Polylineのリダクション(Visvalingam-Whyattアルゴリズム)
Visvalingam-Whyatt アルゴリズムは、各頂点が両隣の点と成す三角形の面積を比べて、小さいものから順に消していくポリラインの間引き。入力は残す頂点の割合(0〜1)で、Douglas-Peucker のような距離のしきい値ではないため、モデルの大きさを気にせず使えるのが便利な点。
処理の流れ

両隣のポイントと成す三角形の面積を比較して、小さいものから順に削除していく。 削除したら前後のポイントの三角形の面積を再計算する。両端の点は消さない。

三角形を視覚化したもの。
コード

ポリライン1本を入力する(後段の Python は点番号をそのまま配列の添字にしているので、複数本には対応していない)。初めに初期値として使う各ポイントの三角形の面積を計算する。外積の長さ=ベクトルの成す平行四辺形の面積、を利用する。両端の点は面積を持たない。
// Run Over: Primitives
int pts[] = primpoints(0, @primnum);
for(int i = 1; i < len(pts)-1; i++)
{
vector prev = point(0, "P", pts[i-1]);
vector next = point(0, "P", pts[i+1]);
vector pos = point(0, "P", pts[i]);
float area = length(cross(next-pos, prev-pos))/2;
setpointattrib(0, "area", pts[i], area);
}
この Python SOP ではポイントを削除せず、残すポイントの番号を Detail 属性 result に配列で書き出している。残す割合は null1 ノードの percent パラメータ(0〜1)から読む。毎回全ポイントを走査するので O(n²) だが、ヒープで最小面積を管理すれば O(n log n) にできる。
node = hou.pwd()
geo = node.geometry()
class PT:
def __init__(self, index, pos, area):
self.index = index
self.pos = pos
self.area = area
pts = []
for pt in geo.points():
pts.append(PT(pt.number(), pt.position(), pt.floatAttribValue("area")))
# 残す割合 0-1
percent = `chs("../null1/percent")`
# 頂点が3つ以上を条件とする
if(len(pts) > 2):
# 最小頂点を2に設定する(直線)
num = int(len(pts) * percent)
if(num < 2):
num = 2
# 削減する目標までループ処理
while(len(pts) > num):
# 最小面積のポイントを探す(両端は除く)
minArea = float('inf')
index = -1
for i in range(1, len(pts)-1):
if(pts[i].area <= minArea):
minArea = pts[i].area
index = i
prev = index-1
next = index+1
# 前後の頂点の面積値を更新する
if(prev > 0):
p0 = pts[prev-1].pos
p1 = pts[next].pos
v0 = p0 - pts[prev].pos
v1 = p1 - pts[prev].pos
area = hou.Vector3.cross(v0, v1)
area = hou.Vector3.length(area)/2
pts[prev].area = area
if(next < len(pts)-1):
p0 = pts[prev].pos
p1 = pts[next+1].pos
v0 = p0 - pts[next].pos
v1 = p1 - pts[next].pos
area = hou.Vector3.cross(v0, v1)
area = hou.Vector3.length(area)/2
pts[next].area = area
pts.pop(index)
result = []
for i in range(0, len(pts)):
result.append(pts[i].index)
# 残すポイントの番号を Detail に保存
geo.addArrayAttrib(hou.attribType.Global, 'result', hou.attribData.Int, 1)
geo.setGlobalAttribValue('result', result)
result に入っていないポイントを削除する。
// Run Over: Points
int pts[] = detail(0, "result");
if(len(pts) > 0 && find(pts, @ptnum) < 0)
removepoint(0, @ptnum);