CGALでのメッシュ生成

以下のサイトのGenerating Meshes with CGALに，関数からメッシュを作ってメッシュにアクセスする方法が書いてある．
http://schmid.dk/wiki/index.php/CGAL

この方法で，メッシュ生成すると，ちゃんとメッシュが作れた．
ユーザーの気持ちを知っているすばらしいサイトだと思う．

Facet_iteratorでcomplex上のfaceにアクセスし，アクセスしたfaceのfirst(四面体)からsecond以外の頂点にアクセスする．そのsecond以外の頂点は，complex上のfaceなので，ちゃんと頂点が取れてくる．ということです．
知らないとできないです．こんなの．

C2t3でメッシュ生成後，以下のコード(上のサイトからの転載)で，メッシュにアクセスできます．

    // for each facet (triangle), print the vertices
    for(C2t3::Facet_iterator fi = c2t3.facets_begin();
            fi != c2t3.facets_end(); fi++) {

       C2t3::Cell_handle cell    = fi->first;
       int opposite_vertex_index = fi->second;
       for(int i = 0; i < 4; i++)
           if(i != opposite_vertex_index) {
               std::cout << "(" << cell->vertex(i)->point().x()
                         << "," << cell->vertex(i)->point().y()
                         << "," << cell->vertex(i)->point().z() << ")"
                         << std::endl;
           }
    }

ありがとうございました．

陰関数でのメッシュ生成は多少時間がかかる

陰関数からメッシュ生成することが素晴らしいと思って開発していたのだけれど，計算時間がかかりすぎることを忘れていた．
ひとつの形状を作るのには，優れているかもしれないが，めちゃくちゃ多くのものを作るには，面白くない．
それに，今回の目的に対して陰関数である必要性が，弱い．
他の方法で，簡単にメッシュ作れるならばそっちを適用した方がかっこいいので，そっちに乗り換える．

ドロネー図の点の参照

CGALは，creationの時点で，グラフを作ります．
ですので，宣言して，insertで点を入れていくだけでグラフが形成されます．
問題はどうやってそのグラフを参照するかだったんですが，cellかvertexのiterator経由でハンドルつかんで，点を順に参照していけばよいです．
ドロネー図をFinite_cells_iterator経由で，四面体の4頂点を自作四面体クラス(Tetrahedron)に格納していく例を挙げます．
自分用メモみたいな感じですが，以下コードです．

 Delaunay dt;
 for(int i=0; i < _siteN; i++){
  dt.insert(Point(_site[i][0], _site[i][1], _site[i][2]));
 }
 //std::cout << dt.number_of_finite_cells() << "\n";
 tetrahedraN = dt.number_of_finite_cells();
 tetrahedra = new Tetrahedron*[tetrahedraN];
 for(int i=0; i < tetrahedraN; i++){
  tetrahedra[i] = new Tetrahedron();
 }

 int cnt = 0;
 for(Delaunay::Finite_cells_iterator fci = dt.finite_cells_begin(); fci != dt.finite_cells_end(); fci++){
  for(int j=0;j < 4;j++){
   Vertex_handle v = fci->vertex(j);
   tetrahedra[cnt]->set(j, 0, v->point().x());
   tetrahedra[cnt]->set(j, 1, v->point().y());
   tetrahedra[cnt]->set(j, 2, v->point().z());
  }
  cnt++;
 }

siteは母点です．tetrahedraのset関数は，インデックスと次元，位置を格納します．

CGALのVoronoi図計算実験

CGALで2次元のボロノイ図を計算させてみました．母点をランダムに配置した後，ボロノイ図を計算し，ボロノイエッジを描画しました．上の図のように描画ができますが，ソースコードはいまいちな気がします．

ボロノイ図の計算にあたって，リファレンスの2D Voronoi Diagram Adaptorを読んで使い方を学ぼうと思ったのですが，思考がついていかず，細部の設計が理解できませんでした．ですので，付属のサンプルのソースコードを追っかけてなんとなく使う作戦にしました．
ボロノイ図を計算させるexamples/Voronoi_diagram_2/vd_2_point_location.cppのコンパイルと実行が確認できたので，このソースを元に簡単に変更していきました．
まずは，インクルード周りですが，元にしたソースのままです．

#include < CGAL/basic.h >

// includes for defining the Voronoi diagram adaptor
#include < CGAL/Exact_predicates_inexact_constructions_kernel.h >
#include < CGAL/Delaunay_triangulation_2.h >
#include < CGAL/Voronoi_diagram_2.h >
#include < CGAL/Delaunay_triangulation_adaptation_traits_2.h >
#include < CGAL/Delaunay_triangulation_adaptation_policies_2.h >

// typedefs for defining the adaptor
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Delaunay_triangulation_2 < K > DT;
typedef CGAL::Delaunay_triangulation_adaptation_traits_2 < DT > AT;
typedef CGAL::Delaunay_triangulation_caching_degeneracy_removal_policy_2 < DT > AP;
typedef CGAL::Voronoi_diagram_2 < DT,AT,AP > VD;

// typedef for the result type of the point location
typedef AT::Site_2 Site_2;
typedef AT::Point_2 Point_2;

typedef VD::Locate_result Locate_result;
typedef VD::Vertex_handle Vertex_handle;
typedef VD::Face_handle Face_handle;
typedef VD::Halfedge_handle Halfedge_handle;
typedef VD::Ccb_halfedge_circulator Ccb_halfedge_circulator;

次にボロノイ図の計算周りですが，ランダムな点を母点とし，母点の位置にある要素をとってきて描画しました．

 double w = 2.5;
 double h = 2.0;

 int n = 100;
 double** sitePos = new double*[n];
 for(int i = 0; i < n; i++){
  sitePos[i] = new double[2];
  sitePos[i][0] = math::myRandom(-0.50)*2.0*w;
  sitePos[i][1] = math::myRandom(-0.50)*2.0*h;
 }

 //2D Voronoi adaptor
 VD vd;
 for(int i = 0; i < n; i++){
  //insertion with insert() func
  vd.insert(Site_2(sitePos[i][0], sitePos[i][1]));
 }

 for(int i = 0; i < n; i++){
  //母点の指定
  Point_2 p(sitePos[i][0], sitePos[i][1]);

  //指定母点上のボロノイ図の要素
  Locate_result lr = vd.locate(p);

  //Face_handleを受け取る
  Face_handle* f = boost::get< Face_handle >(&lr);

  //Face_handleからhalfedgeを受け取る
  Ccb_halfedge_circulator ec_start = (*f)->outer_ccb();

  //Voronoi edgeのスタート地点をecに渡す
  Ccb_halfedge_circulator ec = ec_start;

  do {
   //Voronoi edgeの有無
   if(ec->has_source()){
    //voronoi edgeの端点その1の位置
    double vxs = ec->source()->point().x();
    double vys = ec->source()->point().y();
    if(ec->has_target()){
     //Voronoi edgeの端点その2の位置
     double vxt = ec->target()->point().x();
     double vyt = ec->target()->point().y();
     //描画
     glVertex2d(vxs,vys);
     glVertex2d(vxt,vyt);
    }
   }
   else{
    //無限遠
    std::cout << "inf" << std::endl;
   }
  } while ( ++ec != ec_start );//スタート地点に戻ってくるまでのdo while文
 }

 for(int i = 0; i < n; i++)
  delete[] sitePos[i];
 delete[] sitePos;

ボロノイエッジの描画だけならもっとさらっとしたコードになりそうですが，今はよくわからないことが多いので，これが今の実力では最善です．

ISVD2010気になる

ボロノイ図に関する国際会議ISVD2010の論文募集がcg-mlで流れていた．
ボロノイ図がらみというだけでマニア心をくすぐる．
ボロノイ図がらみならば，募集の領域をCGや計算幾何に限定していない，というのも気になる．

委員を見てたら某JAIST某浅野先生がいた．気になる．