qhull 組み込みメモ(2) delaunay triangulation

qhullの関数を呼び出してドロネー図を計算する方法が理解できたのでメモ．

Qhull internals http://www.qhull.org/html/qh-in.htm#TOC
上記のサイトにいくつかのヒントが記述されているが，読むには，Cの読解力と幾何学の知識が必要だと思う．初心者にはやさしくない．

この図は，3次元空間内に，ランダムに6点プロットし，その点を母点とし，ドロネー図を計算した結果．
母点は黒色，エッジは青色，面は緑色で透過をかけて描画した．

3次元空間内でのドロネー領域は4点で構成されるため，出力のfacetは母点のインデックス4点を指す．

facet->upperdelaunayフラグが立っているドロネー領域は，実際にドロネー図を構成しているものではないので，排除する．
全facet数であるqh num_facetsを見てみると，ドロネー領域の数より多いことがわかる．
そして，ドロネー領域の数は，qh num_goodに格納されている．上記のとおり，qh num_goodはドロネー領域の数であるとともに，facet->upperdelaunayフラグが立っていないfacetの数でもある．またボロノイ(Voronoi)図の頂点数でもある．

qh_pointid関数は引数に頂点の先頭番地を入れると，その頂点のidを返す．
また，逆にqh_point関数は引数にidを入れると，頂点の先頭番地を返す．

FORALLfacetsはfor文を含んだマクロで，今参照ているfacetをfacet_listの先頭から尻尾まで順に参照する．
FOREACHvertex_(facet->vertices)は(括弧内の)facet下のvertexを先頭から尻尾まで参照する．

メモとして，ソースを載せておく．
私は計算用のメモリと描画用のデータメモリを別に作っているので，以下のソースにはそのメモリの移し替えが含まれている．描画用のメモリは，変数の先頭にg_がついている．

g_vertexに頂点を格納してあり，g_indexにドロネー領域のインデックスを格納している．

#if defined(__cplusplus)
extern "C"
{
#endif
#include 
#include 
#include "qhull/qhull.h"
#include "qhull/mem.h"
#include "qhull/qset.h"
#include "qhull/geom.h"
#include "qhull/merge.h"
#include "qhull/poly.h"
#include "qhull/io.h"
#include "qhull/stat.h"
#if defined(__cplusplus)
}
#endif

int main(int argc, char** argv){
 double *points;           /* array of coordinates for each point */ 
 unsigned int ismalloc;           /* True if qhull should free points in qh_freeqhull() or reallocation */ 
 char flags[]= "qhull d Qbb QJ"; /* option flags for qhull, see qh_opt.htm */
 FILE *outfile= stdout;    /* output from qh_produce_output()   
         use NULL to skip qh_produce_output() */ 
 FILE *errfile= stderr;    /* error messages from qhull code */ 
 int exitcode;             /* 0 if no error from qhull */
 facetT *facet;           /* set by FORALLfacets */
 int curlong, totlong;   /* memory remaining after qh_memfreeshort */

 /* initialize dim, numpoints, points[], ismalloc here */
 g_dim = 3;
 g_vertexN = 6;
 g_vertex = new double*[g_vertexN];
 for(int i = 0; i < g_vertexN; i ++){
  g_vertex[i] = new double[3];
  g_vertex[i][0] = (double)rand()/RAND_MAX - 0.5;
  g_vertex[i][1] = (double)rand()/RAND_MAX - 0.5;
  g_vertex[i][2] = (double)rand()/RAND_MAX - 0.5;
 }
 points = new double[g_dim*g_vertexN];
 for(int i = 0; i < g_vertexN; i ++){
  for(int j=0;j<3;j++){
   points[i*3+j] = g_vertex[i][j];
  }
 }
 ismalloc = 1;

 exitcode= qh_new_qhull (g_dim, g_vertexN, points, ismalloc,
  flags, outfile, errfile);
 if (!exitcode) { /* if no error */ 
  facetT *facet;
  vertexT *vertex, **vertexp;
  g_indexN = qh num_good;
  g_index = new int*[g_indexN];
  for(int i=0;iupperdelaunay) {
    //printf ("%d", qh_setsize (facet->vertices));
    int cnt2 = 0;
    FOREACHvertex_(facet->vertices){
     g_index[cnt][cnt2] = qh_pointid (vertex->point);
     cnt2++;
    }
    cnt++;
   }
  }
 }

 qh_freeqhull(!qh_ALL);  
 qh_memfreeshort (&curlong, &totlong);
 if (curlong || totlong)
  fprintf (errfile, "qhull internal warning (main): did not free %d bytes of long memory (%d pieces)\n", 
  totlong, curlong);
}

測定していないので，どのくらい違いが出ているかわからないけれど，qhullの実行ファイルで標準出力させてプログラムに読み込むよりも実行速度が速く感じる．