DOPの描画(2) -Qhullの利用-

前回に引き続きDOPの描画です。多少の変更はあったものの、前回考えていた通り、「DOPを構成する平面の交点を列挙し、その凸包を作る」という方法で無事に描画することができました。ただし、凸包で思いのほか苦労してしまいました。Qhullというライブラリを使ったのですが、その使い方を間違えていて一日中頭を悩ませていました。今後もまた使う機会があると思うので、自分用にメモを残しておきます。知ったかで書いているので、間違っている点などがあったら教えてください。

Qhullについて

Qhullは、凸包(convex hull)、ドローネ三角形分割（Delaunay triangulation）、ボロノイ図（Voronoi diagram）などを作るためのソフトウェアです。その名前からもわかるとおり、凸包のアルゴリズムとしてquickhullを採用しています。3次元以上の凸包もできるようです。Qhull自体は入力された点群に対して処理を行うソフトウェアなのですが、ソースが公開されており、多くのプラットフォームで動くことからライブラリとしても使われています。

使ってみての注意点など

DownloadのページからWindows用のものを落としてきて使いました。バージョンは 2003.1 です。コンパイルの仕方やQhullのソフトの使い方などはReadmeに書かれているので割愛します。
以下、ひっかかったところを箇条書きにしておきます。

・C++で使う場合は"qhull_interface.cpp"を参考

そのままリンクすると関数が見つからないなどのエラーが出る。

・サンプルプログラムは"user_eg.c"と"user_eg2.c"

user_eg.cが基本。user_eg2.cはより複雑だが、細かく使えるらしい。それから、"unix.c"もサンプルかな。

・"qh "は"qh_qh->"か"qh_qh."を表すマクロ

"qh num_vertices"とか"qh num_facets"なんていうのがよく出てくるけど、これはグローバル変数であるqh_qh参照のためのマクロ。これ、気持ち悪い。

・全ての処理は"qh_new_qhull()"から

凸包はもちろん、その他の処理も全て"qh_new_qhull()"で行える。引数には頂点の次元数や頂点の数を渡すが、特に重要なのは”コマンド”。必ず"qhull "で始まる文字列を渡す。コンソールからqhullプログラムを実行する時のオプションと同じものを指定できる。今回は次のようにした。qh_new_qhull( 3, verts.size(), verts[0].vertex, False, "qhull ", NULL, stderr)

・FORALLpointsやFOREACHpoint_で頂点を参照

頂点やインデックスは入り組んでいるので、アクセスするにはFORALL*マクロやFOREACH*_マクロを使う。これらはfor分を隠蔽したマクロ。イテレータとなる変数は前もって宣言しておく必要があるので注意。

・REALfloatを1にする

デフォルトではdoubleで計算される。"user.h"にあるREALfloatを1にすればfloatになる。制度が必要な計算などでは警告が出る。

・qh_ORIENTclockを1にする

デフォルトでは反時計回りらしいんだけど、僕がやったときは逆な感じがした。"user.h"にあるqh_ORIENTclockを1にすれば逆転する。あとで詳しく調べておこう。

・qh_facet3vertex()を使用する

qh_ORIENTclockで指定した通りの頂点インデックスが欲しい場合は、qh_facet3vertex()を使う必要がある。そうじゃないと、順番がばらばらだったりする。

・こんな風に使った

// 凸包の実行
if ( qh_new_qhull( 3, verts.size(), verts[0].vertex, False, "qhull ", NULL, stderr) )
	return;	// 凸包に失敗
// 頂点を取得
pointT *point, *pointtemp;
FORALLpoints
	dop8.verts.push_back( Vector3d(point[0],point[1],point[2]) );
// 面を取得
facetT *facet;
vertexT *vertex, **vertexp;
FORALLfacets
{
	setT *vertices = qh_facet3vertex(facet);	// 逆時計回り順にする
	dop8.idxs.push_back(qh_setsize(vertices));
	FOREACHvertex_(vertices)
		dop8.idxs.push_back( qh_pointid(vertex->point) );
	qh_settempfree(&vertices);
}
qh_freeqhull(!qh_ALL);                   // free long memory
int curlong, totlong;	    // memory remaining after qh_memfreeshort
qh_memfreeshort (&curlong, &totlong);    // free short memory and memory allocator
if (curlong || totlong) 
{
	fprintf (stderr, "qhull internal warning (user_eg, #1): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
	return;
}