自転車生活♪: pythonでSTL形式のデータを作ってみた

gnuplotでは、以下のようなスクリプトで球面を描画したりします。

set parametric # 媒介変数モード
set urange[0:2*pi] # uの範囲を設定
set vrange[0:2*pi] # vの範囲を設定
set ticslevel 0
set hidden3d
set isosample 50
set view equal xyz
splot cos(u/2.0)*cos(v),sin(u/2.0)*cos(v),sin(v) w l

このようなものを3Dモデルとして出力したいということで、 pythonで(u, v)の関数をSTL形式にして出力するスクリプトを書いてみました。ほかの言語で書き換えて遊ぶのにも、ちょうどよい内容と長さかと思います。

import math

def Sphere(u, v):
    x = math.cos(u/2.0)*math.cos(v)
    y = math.sin(u/2.0)*math.cos(v)
    z = math.sin(v)
    return x, y, z

def VectorSub(v1, v2):
    return [e1-e2 for e1, e2 in zip(v1, v2)] 

def VectorDiv(v, a):
    return [e/a for e in v] 

def VectorCross(v1, v2):
    x1, y1, z1 = v1
    x2, y2, z2 = v2
    return [y1*z2-z1*y2, z1*x2-x1*z2, x1*y2-y1*x2]

def VectorAbs(v):
    return math.sqrt(sum([e**2 for e in v]))

def NormalVector(p1, p2, p3):
    v1 = VectorSub(p2, p1)
    v2 = VectorSub(p3, p2)
    vc = VectorCross(v1, v2)
    a = VectorAbs(vc)
    v = VectorDiv(vc, a) if a!=0.0 else False
    return v

def STLfacet(apolygon, inverse=False):
    if inverse:
        points = apolygon[::-1]
    else:
        points = apolygon[:]

    nv = NormalVector(*points)

    if not nv:
        return ""

    l = "facet normal %g %g %g\n" % tuple(nv)
    l = l + "outer loop\n"
    for p in points:
        l = l + "vertex %g %g %g\n" % tuple(p)
    l = l + "endloop\n"
    l = l + "endfacet\n"
    return l

def PolygonsToSTL(name, polygons, inverse=False):
    l = "solid %s\n" % name
    for p in polygons:
        l = l + STLfacet(p, inverse)
    l = l+"endsolid\n"
    return l

def Polygons_of_Function_UV(func):
    nu = 90 
    nv = 90 
    du = 2.0 * math.pi / nu
    dv = 2.0 * math.pi / nv

    u = [du * i for i in range(nu+1)]
    v = [dv * i for i in range(nv+1)]

    polygons = []
    for i in range(nu):
        for j in range(nv):
            p1 = func(u[i], v[j])
            p2 = func(u[i+1], v[j])
            p3 = func(u[i+1], v[j+1])
            p4 = func(u[i], v[j+1])
            polygons.append((p1, p2, p4))
            polygons.append((p4, p2, p3))
    return polygons

def OutputSTL(name, polygons):
    with open('%s.stl' % name, "w") as fp:
        fp.write(PolygonsToSTL(name, polygons, False))

if __name__ == '__main__':
    OutputSTL("Sphere", Polygons_of_Function_UV(Sphere))

出力されたSTLデータを見るには、GLC Playerが便利ですね。また、Windows10には、3D Builderがついているので、そのままでSTLデータを見ることができます。GLC Playerでスナップショットを取ったものが下の図です。

関数を次のものにしてみます。

def Dounut(u, v):
    x = math.cos(u)*(1+0.5*math.cos(v))
    y = math.sin(u)*(1+0.5*math.cos(v))
    z = 0.5*math.sin(v)
    return x, y, z

ドーナツの形になりました。

def Helical(u, v):
    n = 10.0
    a = 0.6
    b = 0.3
    rm = (lambda v: a*b/math.sqrt(b*b*math.cos(v)**2 + a*a*math.sin(v)**2))
    r = (lambda u, v: 1 + rm(v)*math.cos(v+n*u/2.0))
    x = r(u, v) * math.cos(u)
    y = r(u, v) * math.sin(u)
    z = rm(v) * math.sin(v + n*u/2.0)
    return x, y, z