VEXでCatmull-Clark Subdivision Surfaceを実装

はじめに

この記事はこちらの講演に感銘を受けまして自分なりにVEXでCatmull-Clark Subdivision Surfaceを実装してみたという内容です。

あくまで私の実装ですので講演内で紹介されていたものに関してはこちらのチュートリアルを参考にしてください。
Pragmatic VEX

環境

• Houdini 21.0.540

シーンファイル

hipnc

Catmull-Clark Subdivision Surface

アルゴリズムに関してはこちらの動画のシリーズが参考になるかと思います。

式にまとめたものが以下になります。
$F$がポリゴンの中心ポイント、$E$が辺の中心ポイント、$P$が元のポイントです。
$n$は隣接しているポイントまたはポリゴンの数になります。

F=\frac{1}{n}\sum^{n}_{k=1}P_{k}

E=\frac{1}{2}\sum^{2}_{k=1}P_{k}

P_{new}=\frac{1}{n}\left(\frac{1}{n}\sum^{n}_{k=1}F_{k}+\frac{2}{n}\sum^{n}_{k=1}E_{k}+(n-3)P\right)

E_{new}=\frac{1}{2}\left(E+\frac{1}{2}\sum^{2}_{k=1}F_{k}\right)

• 境界の場合
  $P_{neighbor}$は隣接する境界ポイントです。

P_{new}=\frac{1}{8}P_{neighbor1}+\frac{3}{4}P+\frac{1}{8}P_{neighbor2}

E_{new}=\frac{1}{2}\sum^{2}_{k=1}P_{k}

• コーナーの場合

P_{new}=P

既存実装

HoudiniのSOPで使える既存の実装としてSubdivideがあります。
アルゴリズムとしてOpenSubdive Catmull-Clark、Houdini Catmull-Clarkを選択することで対象のジオメトリにCatmull-Clark Subdivision Surfaceを適用することができます。

講演内でも説明がありましたがOpenSubdive Catmull-Clarkは境界の接続処理が行われませんがHoudini Catmull-Clarkでは接続処理が行われます。

実装

前提

今回は閉じたポリゴン以外のプリミティブは考慮しません。
また以下のようなエラーポリゴンの対応も考慮しません。

アルゴリズム

全体のネットワーク

トポロジを変更するノードの色を赤よりの色、それ以外を青よりの色にしています。

1. 分割する面のグループ、元のポイント、元の面のグループを設定
2. uvwの値等必要なアトリビュートを設定
3. 上に同じ
4. 上に同じ
5. ポイント生成に必要なアトリビュートを設定
6. 境界面を設定
7. ポイント生成
8. 生成したポイントと生成元のポイントを紐づけ
9. 面を分割
10. 分割した面のアトリビュートを設定
11. 境界面を再構築
12. 補間用のアトリビュートを設定
13. 上に同じ
14. 上に同じ
15. ポイント、頂点、面それぞれのアトリビュートを補間
16. ポイントを移動
17. 上に同じ
18. 元の面を削除

※コードのコメントに面という単語が出てきますが閉じたポリゴンと同じ意味です。

1. Group

分割部分を指定しているプリミティブグループを_f_subdというグループに名前変更し、元のポイントおよびポリゴンを識別できるように_f_originalと_p_originalという名前のグループを作成します。

識別のためにプリミティブアトリビュートおよびグループに_f_、ポイントアトリビュートおよびグループに_p_、頂点アトリビュートおよびグループに_v_のプリフィックスをつけます。

2. Primitives Wrangle

ポリゴンの重心座標にあたるプリミティブUV値をcentroid_uvwカスタム関数で求め@_f_uvwに記録します。

コード

snippet

// 面の重心にあたるuvwを計算する関数
function vector centroid_uvw(int num_vtx) {
    // 三角ポリゴンの場合
    if (num_vtx == 3) {
        float uv = 1.0 / 3.0;
        vector uvw = set(uv, uv, 0);
        
        return uvw;
    }
    // 四角ポリゴンの場合
    if (num_vtx == 4) {
        vector uvw = set(0.5, 0.5, 0);

        return uvw;
    }

    // それ以外
    vector uvw = set(0, 1, 0);
    
    return uvw;
}

// 面の重心にあたるuvw
v@_f_uvw = centroid_uvw(@numvtx);

3. Vertices Wrangle

後続の処理で使用するため現在の頂点の接続状態を記録します。

コード

snippet

// 現在の面頂点番号
int vtx_idx = vertexprimindex(0, @vtxnum);

int nxt_vtx_idx = vtx_idx + 1; // 次の面頂点番号
int prv_vtx_idx = vtx_idx - 1; // 前の面頂点番号

// 面頂点番号の範囲外処理
if (nxt_vtx_idx == @numvtx) nxt_vtx_idx = 0;
if (prv_vtx_idx < 0) prv_vtx_idx = @numvtx - 1;

i@_v_nxtvtx = vertexindex(0, @primnum, nxt_vtx_idx); // 次の面頂点
i@_v_prvvtx = vertexindex(0, @primnum, prv_vtx_idx); // 前の面頂点

i@_v_nxtpt = vertexpoint(0, i@_v_nxtvtx); // 次の面頂点のポイント
i@_v_prvpt = vertexpoint(0, i@_v_prvvtx); // 前の面頂点のポイント

4. Vertices Wrangle

頂点に優先順位、境界の判定、辺の中心座標のプリミティブuvを記録します。

コード

snippet

// ポリゴンの辺の外周に沿ったuvwを取得する関数
function vector perimeter_uvw(float u; int num_vtx, closed) {
    // uの値の範囲外処理を行う
    if (u > num_vtx || u < 0) 
        u = closed ? u % num_vtx : clamp(u, 0, float(num_vtx));

    // 三角形ポリゴンの処理
    if (num_vtx == 3) {
        if (u < 1.0) return set(u, 0, 0);
        if (u > 2.0) return set(0, 3.0 - u, 0);
        return set(2.0 - u, u - 1.0, 0);
    }
    
    // 四角形ポリゴンの処理
    if (num_vtx == 4) {
        if (u < 1.0) return set(0, u, 0);
        if (u < 2.0) return set(u - 1.0, 1, 0);
        if (u < 3.0) return set(1, 3.0 - u, 0);
        return set(4.0 - u, 0, 0);
    }
    
     return set(u / num_vtx, 0, 0);
}

// 閉じたポリゴン専用のオーバーロード
function vector perimeter_uvw(float u; int num_vtx) {
    return perimeter_uvw(u, num_vtx, 1);
}

int prims[] = array(@primnum);

// 自身とポイントを共有する頂点を走査
foreach (int vtx; pointvertices(0, i@_v_nxtpt)) {
    int prim = vertexprim(0, vtx);

    // 自身が所属する面もしくは分割対象の面でなければ次へ
    if (!inprimgroup(0, "_f_subd", prim) || prim == @primnum) continue;

    int nxt_pt = vertex(0, "_v_nxtpt", vtx); // 次の面頂点のポイント
    int prv_pt = vertex(0, "_v_prvpt", vtx); // 前の面頂点のポイント

    // 自身と紐づいているポイントと同一であれば次へ
    if (nxt_pt != @ptnum && prv_pt != @ptnum) continue;
    
    append(prims, prim);
}

// 次の面頂点との辺の中心に相当する値
float perimeter_u = float(vertexprimindex(0, @vtxnum)) + 0.5;

i@_v_isprimary = min(prims) == @primnum; // 分割対象の面の中で最小の@primnumが自身が所属する面の@primnumと一致するなら優先頂点としてマーク
v@_v_euvw = perimeter_uvw(perimeter_u, @numvtx); // 辺中心のuvw値
i@_v_isboundary = len(prims) == 1; // 分割対象の面が一つしかなければ境界頂点としてマーク

ハーフエッジを使わない理由

単純にアクセスするだけでも高コストです。
おそらく本来データとしてハーフエッジ構造を持っているわけではなくハーフエッジ関連の関数を使用した際にデータを新たにデータを構築していると思われます。
辺を共有する頂点に優先順位をつけるという処理を行うだけでもVertices Wrangleの組み合わせを用いたほうが高速になります。

5. Points Wrangle

Point Generateで用いるアトリビュートを設定します。
また生成されたポイントを格納する_p_newptsアトリビュートも作成します。

コード

snippet

// 分割する面に含まれるポイントをグループ化
i@group__p_subd = 1;

i[]@_p_vtxs = array(); // 優先頂点の配列
i[]@_p_vtxprims = array(); // 優先頂点の所属する面の配列
i[]@_p_vtxisboundary = array(); // 優先頂点が境界かどうかの配列
i[]@_p_prims = array(); // 面頂点番号が0の面の配列

int num_boundaries = 0; // 境界頂点の総数
int num_prims = 0; // 面の総数

// 自身に紐づいている頂点を走査
foreach (int vtx; pointvertices(0, @ptnum)) {
    int prim = vertexprim(0, vtx);

    // 頂点の面が分割対象でないなら次へ
    if (!inprimgroup(0, "_f_subd", prim)) continue;
    
    int prv_vtx = vertex(0, "_v_prvvtx", vtx); // 前の面頂点
    int is_primary = vertex(0, "_v_isprimary", vtx); // 頂点が優先か
    int is_boundary = vertex(0, "_v_isboundary", vtx); // 頂点が境界か
    int prv_is_boundary = vertex(0, "_v_isboundary", prv_vtx); // 前の面頂点が優先か

    // 現在の頂点もしくは次の面頂点番号の頂点が境界であればカウントする
    num_boundaries += is_boundary || prv_is_boundary;

    // 面の数をカウントする
    num_prims += 1;

    // 現在の頂点が最初の面頂点番号であればその面を収集
    if (0 == vertexprimindex(0, vtx)) append(i[]@_p_prims, prim);

    // 現在の頂点が優先であればポイントを作成する
    if (is_primary) {
        append(i[]@_p_vtxs, vtx);
        append(i[]@_p_vtxisboundary, is_boundary);
        append(i[]@_p_vtxprims, prim);
    }
}

// 境界タイプ
i@_p_boundarytype = -1;
if (num_prims == 1) {
    i@_p_boundarytype = 1; // コーナー境界
} else if (num_boundaries > 0) {
    i@_p_boundarytype = 2; // コーナー以外の境界
}
// 境界ポイントをグループ化
i@group__p_boundary = i@_p_boundarytype > 0;

i@_p_num = len(i[]@_p_vtxs) + len(i[]@_p_prims); // 生成するポイントの数
i[]@_p_newpts = array(); // 生成されたポイントを紐づけるための配列

6. Primitives Wrangle

境界ポリゴンを設定します。

コード

int num_pts = 0;
int pts[] = primpoints(0, @primnum);

foreach (int pt; pts) {
    if (inpointgroup(0, "_p_subd", pt)) ++num_pts;
}

// 分割エリアに含まれるポイントが一つでもあれば境界面としてグループ化
i@group__f_boundary = num_pts > 1;

7. Point Generate

ポリゴンの中心ポイント、辺の中心ポイントを作成します。

VEXのaddpointを使わない理由

ポイントの作成に関してはPoint Generateが一番速いと思われます。
作成と同時に元になった@ptnumおよび何番目かのアトリビュートを設定してくれる機能も非常に強力です。
1,494,008のポイントにそれぞれ4個のポイントを作成しアトリビュートを設定するという処理をPoint WrangleとPoint Generateで行うと以下のような速度差になります。

8. Points Wrangle

生成したポイントを生成元のポイントの_p_newptsまたはポリゴンの_f_newptに追加します。

コード

snippet

int vtxs[] = point(1, "_p_vtxs", i@_p_ptnum);
int vtx_prims[] = point(1, "_p_vtxprims", i@_p_ptnum);
int vtx_is_boundary[] = point(1, "_p_vtxisboundary", i@_p_ptnum);

int num_vtxs = len(vtxs);

i@_p_nxtpt = -1;
i@_p_prvpt = -1;

// 辺の中心ポイント
if (i@_p_ptidx < num_vtxs) {
    int vtx = vtxs[i@_p_ptidx];
    int prim = vtx_prims[i@_p_ptidx];

    // 辺中心のuvw
    vector edge_uvw = vertex(1, "_v_euvw", vtx);

    // 接続予定のポイントをそれぞれ設定
    i@_p_nxtpt = vertex(1, "_v_nxtpt", vtx);
    i@_p_prvpt = vertexpoint(1, vtx);
    
    // アトリビュート補間用アトリビュート
    i@_p_srcprim = prim;
    v@_p_uvw = edge_uvw;

    // 境界だった場合にグループ化
    i@group__p_boundary = vtx_is_boundary[i@_p_ptidx];

    // 生成元のポイントに自身の番号を格納する
    setpointattrib(geoself(), "_p_newpts", i@_p_ptnum, @ptnum, "append");
} else { // 面の中心ポイント
    int prims[] = point(1, "_p_prims", i@_p_ptnum);
    int prim = prims[i@_p_ptidx - len(vtxs)];

    // 面中心のuvw
    vector uvw = prim(1, "_f_uvw", prim);

    // アトリビュート補間用アトリビュート
    i@_p_srcprim = prim;
    v@_p_uvw = uvw;

    // 面の中心のポイントとしてグループ化
    i@group__p_center = 1;

    // 面の中心になる新規ポイントはその面に自身の番号を格納する
    setprimattrib(geoself(), "_f_newpt", prim, @ptnum, "set");
}

9. Primitices Wrangle

ポリゴンの頂点数だけポリゴンを新規に作成します。

コード

snippet

int new_pts[] = array();
int new_pts_start[] = array(0);
int cnt_new_pt = prim(0, "_f_newpt", @primnum);

// 既存のポイントから新規ポイントを全て収集
foreach (int pt; primpoints(0, @primnum)) {
    int _new_pts[] = point(0, "_p_newpts", pt);
    
    append(new_pts, _new_pts);
    append(new_pts_start, new_pts_start[-1] + len(_new_pts));
}

// 元のポリゴンの頂点を走査
for (int vtx_idx = 0; vtx_idx < @numvtx; ++vtx_idx) {
    // 次の面頂点番号、前の面頂点番号
    int nxt_vtx_idx = vtx_idx + 1;
    int prv_vtx_idx = vtx_idx - 1;

    // 面頂点番号の範囲外処理
    if (nxt_vtx_idx == @numvtx) nxt_vtx_idx = 0;
    if (prv_vtx_idx < 0) prv_vtx_idx = @numvtx - 1;

    // ポイントに紐づけた新規ポイント
    int crt_new_pts[] = slice(new_pts, new_pts_start[vtx_idx], new_pts_start[vtx_idx+1]);
    int nxt_new_pts[] = slice(new_pts, new_pts_start[nxt_vtx_idx], new_pts_start[nxt_vtx_idx+1]);
    int prv_new_pts[] = slice(new_pts, new_pts_start[prv_vtx_idx], new_pts_start[prv_vtx_idx+1]);

    int crtpt = primpoint(0, @primnum, vtx_idx); // 現在の面頂点のポイント
    int nxtpt = primpoint(0, @primnum, nxt_vtx_idx); // 次の面頂点のポイント
    int prvpt = primpoint(0, @primnum, prv_vtx_idx); // 前の面頂点のポイント

    // 接続先ポイント
    int nxt_new_pt = -1;
    int prv_new_pt = -1;

    // 現在の面頂点のポイントからの頂点から接続先のポイントを検索
    foreach (int new_pt; crt_new_pts) {
        int nxt_pt = point(0, "_p_nxtpt", new_pt);
        
        if (nxtpt == nxt_pt) nxt_new_pt = new_pt;
        if (prvpt == nxt_pt) prv_new_pt = new_pt;
    }

    // 見つからなかった場合、次の面頂点のポイントから接続先のポイントを検索
    if (nxt_new_pt < 0) {
        foreach (int new_pt; nxt_new_pts) {
            int nxt_pt = point(0, "_p_nxtpt", new_pt);
            
            if (crtpt == nxt_pt) {
                nxt_new_pt = new_pt;
                break;
            }
        }
    }

    // 見つからなかった場合、前の面頂点のポイントから接続先のポイントを検索
    if (prv_new_pt < 0) {
        foreach (int new_pt; prv_new_pts) {
            int nxt_pt = point(0, "_p_nxtpt", new_pt);
            
            if (crtpt == nxt_pt) {
                prv_new_pt = new_pt;
                break;
            }
        }
    }

    // 新たなポリゴンを生成
    int pr = addprim(geoself(), "poly", crtpt, nxt_new_pt, cnt_new_pt, prv_new_pt);
}

10. Primitices Wrangle

分割したポリゴンに元の頂点に対応した番号を_f_primidxに設定します。

コード

snippet

// 分割した面としてグループ化
i@group__f_divided = 1;

// 新規に作成した面はすべて四角形かつ面頂点番号2は面中心のポイント
int pt = primpoint(0, @primnum, 2);

i@_f_srcprim = point(0, "_p_srcprim", pt);
i@_f_primidx = findsorted(pointprims(0, pt), @primnum);

11. Primitives Wrangle

境界ポリゴンを再構築します。

コード

snippet

int prim_pts[] = primpoints(0, @primnum);
int new_pts[] = array();

// 元の頂点の配列
i[]@_f_srcvtxs = array();

foreach (int vtx_idx; int crt_pt; prim_pts) {
    append(i[]@_f_srcvtxs, primvertex(0, @primnum, vtx_idx));
    append(new_pts, crt_pt);

    // ポイントに紐づけた新規ポイント
    int crt_new_pts[] = point(0, "_p_newpts", crt_pt);

    // 現在の面頂点のポイント
    int nxt_pt = prim_pts[(vtx_idx + 1) % @numvtx];

    // 接続先ポイント
    int nxt_new_pt = -1;

    // 現在の面頂点のポイントからの頂点から接続先のポイントを検索
    foreach (int new_pt; crt_new_pts) {
        int _nxt_pt = point(0, "_p_nxtpt", new_pt);
        
        if (nxt_pt == _nxt_pt) {
            nxt_new_pt = new_pt;
            break;
        }
    }

    // 見つからなかった場合、次の面頂点のポイントから接続先のポイントを検索
    if (nxt_new_pt < 0) {
        int nxt_new_pts[] = point(0, "_p_newpts", nxt_pt);
        
        foreach (int new_pt; nxt_new_pts) {
            int _nxt_pt = point(0, "_p_nxtpt", new_pt);
            
            if (crt_pt == _nxt_pt) {
                nxt_new_pt = new_pt;
                break;
            }
        }
    }
    
    if (nxt_new_pt < 0) continue;
    
    append(i[]@_f_srcvtxs, -1); // 新規ポイントの頂点は-1として記録
    append(new_pts, nxt_new_pt);
}

int num_vtx = @numvtx;

// ポリゴンを再構築
foreach (int vtx_idx; int new_pt; new_pts) {
    // 頂点があるところは再接続、それ以上は新規追加
    if (vtx_idx < num_vtx) {
        setvertexpoint(geoself(), @primnum, vtx_idx, new_pt);
    } else {
        addvertex(geoself(), @primnum, new_pt);
    }
}

12. Vertices Wrangle

境界以外の新規頂点の補間用アトリビュートを設定します。

コード

snippet

// 新規頂点としてグループ化
i@group__v_new = 1;

int src_prim = prim(0, "_f_srcprim", @primnum);
int prim_idx = prim(0, "_f_primidx", @primnum);

i@_v_srcprim = src_prim; // 元の面
i[]@_v_srcvtxs = array(); // 補間元の頂点の配列
f[]@_v_srcweights = array(); // 補間元のウェイトの配列

int vtx_idx = vertexprimindex(0, @vtxnum); // 面頂点番号
int src_vtxs[] = primvertices(1, src_prim); // 元の頂点配列

if (vtx_idx == 0) {
    // コーナー頂点は元の頂点をそのまま引き継ぐ
    append(i[]@_v_srcvtxs, src_vtxs[prim_idx]);
    append(i[]@_v_srcweights, 1.0);
} else if (vtx_idx == 2) {
    append(i[]@_v_srcvtxs, src_vtxs);

    // 面の中心ポイントの頂点は元の面の頂点の平均にする
    float weight = 1.0 / primvertexcount(1, src_prim);
    foreach (int src_vtx; src_vtxs) append(i[]@_v_srcweights, weight);
} else {
    int src_vtx = src_vtxs[prim_idx];

    int neighbour_vtx = vtx_idx == 1 ? 
        vertex(1, "_v_nxtvtx", src_vtx) : 
        vertex(1, "_v_prvvtx", src_vtx);

    // 辺の中心ポイントの頂点は近傍頂点の平均にする
    append(i[]@_v_srcvtxs, src_vtx);
    append(i[]@_v_srcvtxs, neighbour_vtx);
    i[]@_v_srcweights = {0.5, 0.5};
}

13. Vertices Wrangle

境界の新規頂点の補間用アトリビュートを設定します。

コード

snippet

// 新規頂点としてグループ化
i@group__v_new = 1;

int src_prim = prim(0, "_f_srcprim", @primnum);

i@_v_srcprim = src_prim; // 元の面
i[]@_v_srcvtxs = array(); // 補間元の頂点の配列
f[]@_v_srcweights = array(); // 補間元のウェイトの配列

int src_vtxs[] = prim(0, "_f_srcvtxs", @primnum);
int vtx_idx = vertexprimindex(0, @vtxnum);
int src_vtx = src_vtxs[vtx_idx];

if (src_vtx < 0) {
    // 次の面頂点、前の面頂点
    int nxt_vtx_idx = vtx_idx + 1;
    int prv_vtx_idx = vtx_idx - 1;

    // 面頂点番号の範囲外処理
    if (nxt_vtx_idx == @numvtx) nxt_vtx_idx = 0;
    if (prv_vtx_idx < 0) prv_vtx_idx = @numvtx - 1;

    // 辺の中心ポイントの頂点は近傍頂点の平均にする
    append(i[]@_v_srcvtxs, src_vtxs[nxt_vtx_idx]);
    append(i[]@_v_srcvtxs, src_vtxs[prv_vtx_idx]);
    i[]@_v_srcweights = {0.5, 0.5};
} else {
    append(i[]@_v_srcvtxs, src_vtx);
    append(f[]@_v_srcweights, 1.0);
}

14. Primitives Wrangle

新規ポリゴンの補間用アトリビュートを設定します。

コード

snippet

i[]@_f_srcprims = array(i@_f_srcprim); // 補間元の面の配列
f[]@_f_srcweights = array(1.0); // 補間元のウェイトの配列

15. Attribute Interpolate

新規ポイント、新規頂点、新規ポリゴンのアトリビュートを補間します。
ポイントのみプリミティブuvwを用いその他はウェイトを用いた補間を選択しています。

16. Points Wrangle

新規で作成したポイントを移動します。
境界となるポイント、ポリゴンの中心ポイントは移動させる必要がないため省いています。

コード

snippet

int pts[] = array();
vector pos = 0;

foreach (int prim; pointprims(0, @ptnum)) {
    int pt = primpoint(0, prim, 2);
    if (find(pts, pt) < 0) append(pts, pt);
}

foreach (int pt; pts) {
    pos += point(0, "P", pt);
}

pos /= len(pts);

// 辺中心のポイントに式を適用させて移動
@P = lerp(@P, pos, 0.5);

17. Points Wrangle

元のポイントを移動します。

コード

snippet

if (inpointgroup(0, "_p_boundary", @ptnum)) {
    // コーナーではない境界ポイントの処理
    if (i@_p_boundarytype == 2) {
        int vtx_idxs[] = {1, 3}; // 分割した面の面頂点番号1と3は辺の中心ポイント
        int pts[] = array();        

        // 接続している面を走査
        foreach (int prim; pointprims(0, @ptnum)) {
            // 分割した面でなければ次へ
            if (!inprimgroup(0, "_f_divided", prim)) continue;
            
            foreach (int vtx_idx; vtx_idxs) {
                int pt = primpoint(0, prim, vtx_idx);

                // 境界ポイントでなければ次へ
                if (!inpointgroup(0, "_p_boundary", pt)) continue;
                
                int nxt_pt = point(0, "_p_nxtpt", pt); // 次の境界ポイント
                int prv_pt = point(0, "_p_prvpt", pt); // 前の境界ポイント
                
                if (nxt_pt != @ptnum && find(pts, nxt_pt) < 0) append(pts, nxt_pt);
                if (prv_pt != @ptnum && find(pts, prv_pt) < 0) append(pts, prv_pt);
            }
        }

        // 隣接ポイントの平均座標
        vector pos = 0;
        
        foreach (int pt; pts) {
            pos += point(1, "P", pt);
        }

        // 境界の式を適用させて移動
        @P = ((pos / len(pts)) + (@P * 3)) / 4;
    }
} else {
    int num_prim = 0;

    vector cnt_pos = 0; // 面中心ポイントの重心座標
    vector edge_pos = 0; // 辺中心ポイントの重心座標

    // 接続している面を走査
    foreach (int prim; pointprims(0, @ptnum)) {
        if (!inprimgroup(0, "_f_divided", prim)) continue;
        
        int cnt_pt = primpoint(0, prim, 2); // 新規で作成した面の中心ポイント
        int edge_pt = primpoint(0, prim, 1); // 新規で作成した辺の中心ポイント
        
        cnt_pos += point(1, "P", cnt_pt);
        edge_pos += point(1, "P", edge_pt);
        
        num_prim += 1;
    }

    cnt_pos /= num_prim;
    edge_pos /= num_prim;

    // 式を適用させて移動
    @P = (cnt_pos + (2 * edge_pos) + ((num_prim - 3) * @P)) / num_prim;
}

18. Primitives Wrangle

分割前のポリゴンを削除します。

コード

snippet

// 最後の引数に0指定で面のみ削除
removeprim(geoself(), @primnum, 0);

removeprimを用いたほうが若干高速だったためこちらを採用していますがBlastで問題ないと思います。

以下にSubdivide - OpenSubdiv Catmull-Clarkと本実装を適用させたジオメトリを並べています。

重ねるとポイントの位置が一致しています。

1~18まで順に適用している様子です。

RubberToyの一部のみを分割してみました。

境界ポリゴンもしっかり接続されているのがわかります。

処理速度

4種類のジオメトリをそれぞれSubdivide - OpenSubdiv Catmull-Clark、 Subdivide - Houdini Catmull-Clark および本実装で分割しその際の速度をパフォーマンスモニターで計測しました。

	ポイント数	ポリゴン数
A	2,886	2,889
B	12,874	12,854
C	90,085	89,942
D	1,494,008	2,241,391

ジオメトリ	Subdivide OpenSubdiv Catmull-Clark	Subdivide Houdini Catmull-Clark	本実装
A	0.008s	0.029s	0.028s
B	0.036s	0.149s	0.073s
C	0.265s	1.164s	0.461s
D	4.097s	11.477s	7.8s

おおむねSubdivide - OpenSubdiv Catmull-Clarkの200%程度の速度で動作しています。

Dを本実装で分割した際のパフォーマンスモニタの結果が以下になります。

トポロジの変更の処理が特に重く、次いで頂点への操作を行う処理が重いです。

おわりに

Catmull-Clark Subdivision Surfaceそのものの実装より面分割の実装に苦労しました。
現在の実装ではクリーズやuvのバイリニア以外の変位ができていないため今後実装していきたいです。

参考文献

• Implementing Catmull-Clark Subdivision - Introduction & Rules
• Catmull-Clark subdivision surface のアルゴリズム解説
• Subdivision Exterior Calculus for Geometry Processing Supplemental Material: Subdivision Rules