Skip to content

Instantly share code, notes, and snippets.

@tnakaicode

tnakaicode/qiita-002.md

Last active August 15, 2020 08:56
Show Gist options
  • Save tnakaicode/f7faa6783664d9f6e5c2324d2bb663e1 to your computer and use it in GitHub Desktop.
Save tnakaicode/f7faa6783664d9f6e5c2324d2bb663e1 to your computer and use it in GitHub Desktop.
Download ZIP
qiita article 002
Raw
qiita-002.md

pyhonocc-core(OpenCASCADE)でSpline Surfaceを作ります。ついでにSTEPで吐き出します。

Surface生成コード

コードはこれです。

赤のSurfaceと青のSurfaceは形状自体は同じものですが、場所をgp_Ax3を使って変えています。

pic

肝になる部分だけ抜粋します。

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として突出した形状にします。

pic

真ん中あたりが異常に飛び上がった形状になります。飛び出した部分に引っ張られて周辺の形状が変形しているのは、Splineを使う以上はどうしても発生します。
こんな異常な形状でも、(正しいかどうかは別として)生成しきってくれるのがOpenCASCADEの強味でもあります。

Raw
qiita_002-1.py
import numpy as np
import matplotlib.pyplot as plt
import sys
import os
import time
from optparse import OptionParser
sys.path.append(os.path.join("../"))
from base import plot2d
import logging
logging.getLogger('matplotlib').setLevel(logging.ERROR)
from OCC.Display.SimpleGui import init_display
from OCC.Core.gp import gp_Pnt, gp_Vec, gp_Dir
from OCC.Core.gp import gp_Ax1, gp_Ax2, gp_Ax3
from OCC.Core.gp import gp_Trsf
from OCC.Core.TColgp import TColgp_Array1OfPnt, TColgp_Array2OfPnt
from OCC.Core.TopoDS import TopoDS_Shape, TopoDS_Compound
from OCC.Core.TopLoc import TopLoc_Location
from OCC.Core.GeomAPI import GeomAPI_PointsToBSplineSurface
from OCC.Core.GeomAbs import GeomAbs_C0, GeomAbs_C1, GeomAbs_C2
from OCC.Core.GeomAbs import GeomAbs_G1, GeomAbs_G2
from OCC.Core.GeomAbs import GeomAbs_Intersection, GeomAbs_Arc
from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_MakeFace
from OCC.Extend.DataExchange import write_step_file
def set_trf(ax1=gp_Ax3(), ax2=gp_Ax3()):
trf = gp_Trsf()
trf.SetTransformation(ax2, ax1)
return trf
def set_loc(ax1=gp_Ax3(), ax2=gp_Ax3()):
trf = set_trf(ax1, ax2)
loc = TopLoc_Location(trf)
return loc
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
if __name__ == '__main__':
argvs = sys.argv
parser = OptionParser()
parser.add_option("--dir", dest="dir", default="./")
parser.add_option("--pxyz", dest="pxyz",
default=[0.0, 0.0, 0.0], type="float", nargs=3)
opt, argc = parser.parse_args(argvs)
print(opt, argc)
display, start_display, add_menu, add_function = init_display()
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
face = spl_face(*mesh, data)
write_step_file(face, "qiita_002-1.stp")
display.DisplayShape(face, color="RED", transparency=0.75)
axs = gp_Ax3()
axs.SetLocation(gp_Pnt(0, 0, 100))
axs.SetDirection(gp_Dir(0.5, 0.5, 1))
face = spl_face(*mesh, data, axs)
display.DisplayShape(face, color="BLUE", transparency=0.75)
display.FitAll()
display.View.Dump("qiita_002-1.png")
start_display()
Raw
qiita_002-2.py
import numpy as np
import matplotlib.pyplot as plt
import sys
import os
import time
from optparse import OptionParser
sys.path.append(os.path.join("../"))
from base import plot2d
import logging
logging.getLogger('matplotlib').setLevel(logging.ERROR)
from OCC.Display.SimpleGui import init_display
from OCC.Core.gp import gp_Pnt, gp_Vec, gp_Dir
from OCC.Core.gp import gp_Ax1, gp_Ax2, gp_Ax3
from OCC.Core.gp import gp_Trsf
from OCC.Core.TColgp import TColgp_Array1OfPnt, TColgp_Array2OfPnt
from OCC.Core.TopoDS import TopoDS_Shape, TopoDS_Compound
from OCC.Core.TopLoc import TopLoc_Location
from OCC.Core.GeomAPI import GeomAPI_PointsToBSplineSurface
from OCC.Core.GeomAbs import GeomAbs_C0, GeomAbs_C1, GeomAbs_C2
from OCC.Core.GeomAbs import GeomAbs_G1, GeomAbs_G2
from OCC.Core.GeomAbs import GeomAbs_Intersection, GeomAbs_Arc
from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_MakeFace
from OCC.Extend.DataExchange import write_step_file
def set_trf(ax1=gp_Ax3(), ax2=gp_Ax3()):
trf = gp_Trsf()
trf.SetTransformation(ax2, ax1)
return trf
def set_loc(ax1=gp_Ax3(), ax2=gp_Ax3()):
trf = set_trf(ax1, ax2)
loc = TopLoc_Location(trf)
return loc
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
if __name__ == '__main__':
argvs = sys.argv
parser = OptionParser()
parser.add_option("--dir", dest="dir", default="./")
parser.add_option("--pxyz", dest="pxyz",
default=[0.0, 0.0, 0.0], type="float", nargs=3)
opt, argc = parser.parse_args(argvs)
print(opt, argc)
display, start_display, add_menu, add_function = init_display()
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
face = spl_face(*mesh, data)
write_step_file(face, "qiita_002-2.stp")
display.DisplayShape(face, color="RED", transparency=0.75)
display.FitAll()
display.View.Dump("qiita_002-2.png")
start_display()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment