pyhonocc-core(OpenCASCADE)でSpline Surfaceを作ります。ついでにSTEPで吐き出します。
Surface生成コード
コードはこれです。
赤のSurfaceと青のSurfaceは形状自体は同じものですが、場所をgp_Ax3を使って変えています。
肝になる部分だけ抜粋します。
def spl_face(px, py, pz, axs=gp_Ax3()):
nx, ny = px.shape
pnt_2d = TColgp_Array2OfPnt(1, nx, 1, ny)
for row in range(pnt_2d.LowerRow(), pnt_2d.UpperRow() + 1):
for col in range(pnt_2d.LowerCol(), pnt_2d.UpperCol() + 1):
i, j = row - 1, col - 1
pnt = gp_Pnt(px[i, j], py[i, j], pz[i, j])
pnt_2d.SetValue(row, col, pnt)
#print (i, j, px[i, j], py[i, j], pz[i, j])
api = GeomAPI_PointsToBSplineSurface(pnt_2d, 3, 8, GeomAbs_G2, 0.001)
api.Interpolate(pnt_2d)
face = BRepBuilderAPI_MakeFace(api.Surface(), 1e-6).Face()
face.Location(set_loc(gp_Ax3(), axs))
return face
px, py, pzは(n x n)のnp.arrayをnp.meshgridで生成して、適当な曲面のGridデータをpzとして入力します。
あとは、TColgp_Array2OfPntというOpenCASCADEのgp_Pntの2次元配列に代入していき(for loop部分)、GeomAPI_PointsToBSplineSurfaceに食わせます。
GeomAPI_PointsToBSplineSurfaceが実質的なSpline Surfaceを生成する関数ですが、詳しくはこちらを読んでください。
TColgp_Array2OfPntのあとの2つの数値(3,8)はSplineを生成する近似多項式の際の最小次数と最大次数です。
GeomAbs_G2は滑らかさを定義しています。ほかにもGeomAbs_C0, GeomAbs_C1, GeomAbs_C2, GeomAbs_G1などがあります。G2とういのは2階微分可能ということになります。
最後の0.001は近似の際の最大誤差で、入力の点群とSurfaceの誤差がこの値を下回るようにSurfaceを生成します。
1点だけ異常な値を入れてみる
きれいな点群データを入れても面白くないので、生成したデータに異常値を入れてみます。
コードはこちらです。
データ生成部分だけ抜粋します。
px = np.linspace(-1, 1, 100) * 100 + 50
py = np.linspace(-1, 1, 200) * 100 - 50
mesh = np.meshgrid(px, py)
data = mesh[0]**2 / 1000 + mesh[1]**2 / 2000
data[100, 50] = 100.0
(200 x 100)の2次元arrayの大体真ん中あたり(100,50)を100.0として突出した形状にします。
真ん中あたりが異常に飛び上がった形状になります。飛び出した部分に引っ張られて周辺の形状が変形しているのは、Splineを使う以上はどうしても発生します。
こんな異常な形状でも、(正しいかどうかは別として)生成しきってくれるのがOpenCASCADEの強味でもあります。