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);

← 記事一覧へ