0. はじめに

　今回はドロネー三角形分割を Houdini の VEX で実装していきます。
　Houdini 初心者なので間違っている箇所があったら優しく教えてください。

0-1. 動作環境

OS: RedHat Enterprise Linux 9.4
Houdini: 20.5.278 (Production Build) Education Edition

0-2. 免責

この記事をもとに被った不利益について私は一切責任を負いません．

1. アルゴリズムの概要

　ドロネー三角形分割とは、与えられた頂点の集合に対して最小の内角が最大になるように ( 別の言い方をすると prim がなるべく正三角形に近くなるように ) メッシュを張るアルゴリズムです。
　アルゴリズムは以下のように記述できます。

三角形 $T$ を入力ポイント全体の集合 $P$ を覆う三角形で初期化する。
$P$ からランダムなポイント $p$ を選択する。
$p$ を含む $T$ 内の prim $t$ を探索する。
$t$ を $p$ を使って分割する。
違反辺がなくなるまで辺フリップを行う。
$P$　に含まれるすべてのポイントが選択されるまで2-5を繰り返す。
$T$ に属する prim を削除する。

　上記の特徴から地形の計測データの処理などにも使われているらしいです。
　応用すると特徴量を残すメッシュリダクションなどにも使えるのではと思っています。

具体的な処理

2. Houdini での実装

2-1. ヘッダファイルの準備

　処理で使用する関数や構造体は複数の場所で利用する可能性があることとプログラムの保守性から、ヘッダファイルとして $HOUDINI_PATH/vex/include/ ディレクトリに配置します。

$HOUDINI_PATH は Houdini Enviroment に初期化されているシェル上で
hconfig | grep HOUDINI_PATH と入力することで確認できます。

　配置したファイルは Wrangle 上で #include <ファイル名.h> と入力すると、内部の定義を利用できます。
<ファイル名.h> の部分は .txt でも .vex でもファイルの絶対パスでもいいらしいです。(要検証)

2-1-1. エッジを表現する構造体の定義

　このアルゴリズムではエッジに対して様々な処理を行っているが VEX ではエッジに対して操作を行うことができません。
　そこでまず思いつく実装が、関数の引数として primnum を受けて関数内部で処理する頂点の組を解決する方法です。
　しかし、アルゴリズム上操作する対象をエッジにする設計のほうが自然なので、関数の引数はエッジを受けたいです。
そこで今回はエッジを表現する構造体を定義して処理で使用します。

　実際の実装は以下の通りです。

$HOUDINI_PATH/vex/include/edge.h

struct EDGE{
    int first, second;

    int first(){
        return min(first, second);
    }

    int second(){
        return max(first, second);
    }
}

以下は使用例です。詳しくは VEX language reference : : Structs をどうぞ。

// 構造体を作成する。
EDGE e1 = {0, 1};					//	構造体を明示的に作成する。
EDGE e2 = EDGE(hoge, huga);		    //	コンストラクタを使えば変数が使える。

int edgeSecond = e1->second();		//	矢印オペレータで構造体関数をコール。
printf(edgeSecond);					//	これは 1 がプリントされる。

2-1-2. ユーティリティ関数の定義

処理で使う関数を定義します。

$HOUDINI_PATH/vex/include/utility.h

 //    edge と追加した point とチェックする point を受けて、
 //    edge が違反辺なら 1 を、違反辺でないなら 0 を返す関数
function int is_illegal(EDGE e; int pivotPoint, checkPoint) {
   //	外接円の中心を計算する
   int hedge = pointhedge(0, e.first, e.second);
   int prim = hedge_prim(0, hedge);
   int points[] = primpoints(0, prim);

   if(find(points, pivotPoint) < 0) {
       hedge = hedge_nextequiv(0, hedge);
       prim = hedge_prim(0, hedge);
       points = primpoints(0, prim);
   }

   vector pos1 = point(0, "P", points[0]);
   vector pos2 = point(0, "P", points[1]);
   vector pos3 = point(0, "P", points[2]);

   float c = 2.0 * ((pos2.x - pos1.x) * (pos3.z - pos1.z) - (pos2.z - pos1.z) * (pos3.x - pos1.x));
   float x = ((pos3.z - pos1.z) * (pos2.x * pos2.x - pos1.x * pos1.x + pos2.z * pos2.z - pos1.z * pos1.z)
   				+ (pos1.z - pos2.z) * (pos3.x * pos3.x - pos1.x * pos1.x + pos3.z * pos3.z - pos1.z * pos1.z))/c;
   float y = ((pos1.x - pos3.x) * (pos2.x * pos2.x - pos1.x * pos1.x + pos2.z * pos2.z - pos1.z * pos1.z)
   				+ (pos2.x - pos1.x) * (pos3.x * pos3.x - pos1.x * pos1.x + pos3.z * pos3.z - pos1.z * pos1.z))/c;

   vector centerPos = set(x, 0, y);

   //    外接円の半径を求める
   float radius = distance(centerPos, pos1);

   //	チェックする point が外接円の内部にいるかチェックする
   string condition = "@ptnum == " + itoa(checkPoint);
   int nearpoint = nearpoint(0, condition, centerPos, radius);
   int illegal_flg = nearpoint == -1 ? 0 : 1;

   return illegal_flg;
}


//	prim と edge を受けて、 prim 中で edge に属さない point を返す関数
function int findOppositePoint(int primNum; EDGE e) {
   int points[] = primpoints(0, primNum);
   int oppositePoint = points[0] + points[1] + points[2] - e.first - e.second;
   return oppositePoint;
}


//	prim と edge を受けて、 prim と edge を共有する反対側の prim を返す関数
function int findOppositePrim(int primNum; EDGE e) {
   int hedge = pointhedge(0, e.first, e.second);
   hedge = hedge_prim(0, hedge) == primNum ? hedge_nextequiv(0, hedge) : hedge;
   int oppositePrim = hedge_prim(0, hedge);

   return oppositePrim;
}


//	edge と追加した point とチェックした point を受けて、edge をフリップする関数
function void flipEdge(EDGE e; int pointA, pointB) {
   int hedge = pointhedge(0, e->first(), e->second());
   int prim = hedge_prim(0, hedge);
   removeprim(0, prim, 0);
   hedge = hedge_nextequiv(0, hedge);
   prim = hedge_prim(0, hedge);
   removeprim(0, prim, 0);

   addprim(0, "poly", pointA, e->first(), pointB);
   addprim(0, "poly", pointB, e->second(), pointA);
}

2-2. CVexを使った再帰関数の実装

　処理の2-5は再帰関数を使って実装するわけですが、VEX のユーザ定義関数はコンパイル時にインライン展開されるため、再帰関数を定義することはできません。
$cf.$VEX language reference : : User-defined functions

　変わりに CVex のシェーダコールを使って実装します。

$HOUDINI_PATH/vex/CVex/legalize.vfl

#include <assert.h>
#include <edge.h>
#include <utility.h>
#define println(label, value) printf("%s: %f\n", label, value)

//	再帰的に違反辺をフリップする
cvex Legalize(EDGE edge = {0, 0}; int pivotPoint = 0; export int flipCount = 0) {
	//	oppositePoint を解決する
	int hedge = pointhedge(0, edge->first(), edge->second());
	int checkPrim = hedge_prim(0, hedge);
	int oppositePrim = findOppositePrim(checkPrim, edge);

	int pointA = findOppositePoint(checkPrim, edge);
	int pointB = findOppositePoint(oppositePrim, edge);

//	edge が違反辺なら違反辺が出てこなくなるまでフリップする
	if(checkPrim != oppositePrim && pointA != pointB) {
		if(is_illegal(edge, pointA, pointB) == 1) {
			flipEdge(edge, pointA, pointB);
			flipCount += 1;

			int pointC = pivotPoint != pointA ? pointA : pointB;
			EDGE e1 = EDGE(edge->first(), pointC);
			EDGE e2 = EDGE(edge->second(), pointC);
			int flipCountA = 0;
			int flipCountB = 0;
			Legalize("edge", e1, "pivotPoint", pivotPoint, "flipCount", flipCountA);
			Legalize("edge", e2, "pivotPoint", pivotPoint, "flipCount", flipCountB);
		}
	}
}

引数の export キーワードによるアウトプット指定は必須です。

シェーダをコールする際の引数のキーワード指定は必須です。

　プログラムを作成後、vcc を使ってコンパイルすることでエラーチェックします。まず Houdini Enviroment で初期化されているシェルで $HOUDINI_PATH/vex/CVex ディレクトリに移動します。
　その後、vcc ./legalize.vfl と入力しコンパイルします。
　成功すると、同じディレクトリにLegalize.vex というファイルが作成され、エラーがあると失敗します。

vcc で --vex-output オプションを指定しなかったときの出力ファイル名はコンテキスト名( ここでいうcvexキーワードのあとの文字列 -> Legalize)と等しくなります。
詳しくは　Vex compiler を参照してください。

　Wrangle 内で import Legalize; と入力し、シェーダを現行コンテキストに読み込むことでシェーダをコールすることができます。

シェーダはコールするだけでは実行されず、export キーワードで指定した出力変数を評価することで初めて実行されます。

2-3. Houdini 内での処理

2-3-1. ポイントを覆う三角形を作る

　Bound と Wrangle を使って適当に作ります。
　計算は厳密じゃないです。あまり当てにしないでください。

create_huge_tri

vector center = point(1, "P", 0);
vector start = point(2, "P", 0);

float radius = distance(center, start);

int pt1 = addpoint(0, set(center.x - sqrt(3) * radius, 0,  -radius +  (sqrt(3)/5)));
int pt2 = addpoint(0, set(center.x + sqrt(3) * radius, 0,  -radius + (sqrt(2)/5)));
int pt3 = addpoint(0, set(center.x, 0,  radius + (sqrt(3))));

int prim = addprim(0, "poly", pt1, pt2, pt3);

setprimattrib(0, "id", 0, "0", "set");

　あと、Point Cloud の頂点番号をランダムにしておきます。

2-3-2. ポイントを追加分割する

　For each Feedback Loop を使います。Iteration 数はポイント数分です。

作業時はiterations パラメータに $F を入れるとやりやすかったです。
Feedback Loopは累積処理ですが、Single PassはそのIterationしか評価しないので結果が異なり、あまり参考になりません。
iterations パラメータに $F を入れるとタイムスライダを使ってIterationの発展が確認できるので便利です。

　Wrangle やってることは単純で、追加したい prim を削除して張り直しているだけです。

Detail Wrangle : : Annexing_point_to_the_huge_tri

//annexing point
vector pos = point(1, "P", 0);
int pt = addpoint(0, pos);
setpointgroup(0, "__tmp_annexingPoint", pt, 1, "set");

//search prim for annexing
vector tmp;
int deletePrim = -1;
xyzdist(0, pos, deletePrim, tmp);

//restructure mesh
int points[] = primpoints(0, deletePrim);
push(i[]@__targetPrimPoint, points);
removeprim(0, deletePrim, 0);

for(int i = 0; i < len(points); i++){
    int prim = addprim(0, "poly", pt, points[i], points[(i + 1) % 3]);
}

2-3-3. 違反辺をなくす

　まずは３つの辺に対してチェックしたいのでForeach Feedback の中に Wrangle を入れます。 Iterations は 3 です。

Detail Wrangle : : legalize1

import Legalize;
#include <edge.h>
#include <utility.h>
#define println(label, value) printf("%s: %f\n", label, value)
int iter = detail(1, "iteration", 0);
int annexingPoint = expandpointgroup(0, "__tmp_annexingPoint")[0];

int pointA = i[]@__targetPrimPoint[iter];
int pointB = i[]@__targetPrimPoint[(iter + 1) % 3];

setpointgroup(0, "__tmp_edge", pointA, 1);
setpointgroup(0, "__tmp_edge", pointB, 1);

println("Point A--------------", pointA);
println("Point B--------------", pointB);

int flipCountA = 0;

EDGE edgeA = EDGE(pointB, pointA);
Legalize("edge", edgeA, "pivotPoint", annexingPoint, "flipCount", flipCountA);
println("Flip Count A", flipCountA);

#defineでマクロを定義することができます。
#define println(label, value) printf("%s: %f\n", label, value)とすることでややタイプ数を減らせます。
詳しくはVex compiler : : Pre-processorを見てください。

2-3-4. いらないメッシュの削除

delete で hugeTri グループのポイントを削除します。

3. 結果

実装できました。

4. 参考文献

ドロネー三角形分割の期待最速アルゴリズム
 6.886 Delaunay Triangulation

5. おわりに

　CVex を使って Wrangle 内で汎用計算を行う系の情報がほぼなく、ドキュメントにも書いてないハマりポイントで結構時間を消費してしまいました。逆に本体の実装はそこまで時間をかけずに済みました。
　普通のプログラムでやろうとするともうちょっと複雑になってしまうところ、 VEX ならジオメトリを操作するための便利関数が揃っているのでいいですね。
　今度はもう少し複雑な計算幾何系のアルゴリズムを実装してみたいと思います。

　拙い部分もありますが何かの参考になれば幸いです。