CGALの幾何演算
幾何形状に演算を行うことが可能です．
三角分割(Triangulation)
まあ何をするのにも立方体以上の多面体では何もできない. まず, 領域を三角錐で分解することが重要である．一般の図形はTriangulationにより三角が可能である．とりあえず立方体1個 Q を分解してみよう.
#include <CGAL/Triangulation_3.h>
...
std::vector<CGAL::Point_3<Kernel>> points;
for(auto it=Q.vertices_begin();it!=Q.vertices_end();it++) points.push_back(it->point());
CGAL::Triangulation_3<Kernel> Tri(points.begin(),points.end());
std::cout << "Points= " << Tri.number_of_vertices() << " cells= " << Tri.number_of_cells()
<< " edges= " << Tri.number_of_edges() << " facets= " << Tri.number_of_facets()
<< std::endl;
実行してみると, 立方体1個は セル18個, 辺の総数27, 面の総数36であるそうだ．なんか多いが，これは, 周囲のR3空間も三角化されているためである.
std::cout << "Finite Cells= " << Tri.number_of_finite_cells()
<< " Finite Edges= " << Tri.number_of_finite_edges()
<< " Finite Facets= " << Tri.number_of_finite_facets() << std::endl;
これで, まあ普通のバカが考えるセル数などを取得できる．この場合, セルが6, 辺が19, 面が18だそうである．
バカはセルの一覧が見たくなるし体積が欲しかったり
やはりどんなふうに分割したか，一覧が見たい．どうしても．というわけで
int n=0;
for(auto it=Tri.finite_cells_begin();it!=Tri.finite_cells_end();it++) {
n++;
std::cout << "Cell #" << n;
auto Tra=Tri.tetrahedron(it);
std::cout << "(" << it->vertex(0)->point() << ")(" << it->vertex(1)->point() << ")(" <<
it->vertex(2)->point() << ")(" << it->vertex(3)->point() << ")V= " << Tra.volume() << std::endl;
}
これで調べることができる．そうそう，マニュアルのどこにも書いてない気がするけれど
Cell_handler とは, Cell_iterator と同じものだよ
というのがネットの噂で流れている．実際問題ないみたいだ．上の例では，極めてバカなので三角錐Traを構築した上，その体積も求めている. 結果も見せよう：
Cell #1(1 1 1)(1 1 0)(1 0 0)(0 1 1) with volume= 0.166667
Cell #2(1 0 1)(1 1 1)(1 0 0)(0 0 1) with volume= 0.166667
Cell #3(1 1 0)(0 0 0)(1 0 0)(0 1 1) with volume= 0.166667
Cell #4(1 1 0)(0 1 0)(0 0 0)(0 1 1) with volume= 0.166667
Cell #5(0 0 1)(1 1 1)(1 0 0)(0 1 1) with volume= 0.166667
Cell #6(0 0 0)(0 0 1)(1 0 0)(0 1 1) with volume= 0.166667
つまり，このようになっている：
もともと立方体には12の辺があり, 6つの長方形面が2分割されるので6辺増加, さらに中央面に1辺追加されて, 合計で辺が19個になっている. 三角錐は4面体で6個あるので24面あるが, 元の立方体表面である12面を除けば重複しているので, 12+6=18面存在する．
それはそうと, 上から分かるように, 点をインクリメントして三角化するので, 凸物体しか正しくないのではないか？例えば
こんなのをそのままTriangulationすると, 右のように，穴が埋まってconvex hull化しているぞ．これじゃ体積もまともに出やしねえ．これは後で解決しますよ．
ある点がどのセルにあるか
点を指定して, 分割したセルを特定できます：
Triangulation::Locate_type Type;
int li,lj;
Point here(0,3,0);
Triangulation::Cell_handle CH=Tri.locate(here,Type,li,lj);//場所を指定してType, (li,lj)を取得
switch (Type) {
case Triangulation::VERTEX:
std::cout << "(" << here << ") on the vertex #= " << li << std::endl;break;
case Triangulation::EDGE:
std::cout << "(" << here << ") on the edge between= " << li << " - " << lj <<std::endl;break;
case Triangulation::FACET:
std::cout << "(" << here << ") on the facet #= " << li << std::endl;break;
case Triangulation::CELL:
std::cout << "(" << here << ") in the cell" << std::endl; break;
case Triangulation::OUTSIDE_AFFINE_HULL:
std::cout << "(" << here << ") outside" << std::endl;break;
case Triangulation::OUTSIDE_CONVEX_HULL:
std::cout << "(" << here << ") fucking outside" << std::endl;break;
}
まあ使い方は上のプログラムの通りである．
アフィン変換
あふぃーん！変換が可能である.
- 一般の図形のうち, 三角錐はアフィン変換が可能です．
- NEF図形は，全てにおいてアフィン変換が可能です．
アホくさいのでNEF図形の場合だけ示します．Nef_Affがアフィン写像行列で, もちろん同次座標系を用いるのだから4x4なんですが, 最終行ベクトルは(0,0,0,hw) なので, 要素はhwしか指定できないのよね．
#include <CGAL/Nef_polyhedron_3.h>
#include <CGAL/Aff_transformation_3.h>
using Nef_Polyhedron=CGAL::Nef_polyhedron_3<Kernel>;
using Nef_Aff=Nef_Polyhedron::Aff_transformation_3;
Nef_Polyhedron N(P); //Pは普通のPolyhedronです
//Affine変換行列を1行づつ指定します．hwは同次座標系の1です.
double hw=1.0;
Nef_Aff Mat(v[0],0.,0.,from[0], v[1],width[1],0.,from[1],v[2],0.,width[2],from[2],hw);
N.transform(Mat);
この例では, (y,z(平面の矩形 (0,width[1]) x (0,width[2]) を単位時間内に通過した速度(v[0],v[1],v[2])の分子が存在する領域を作成しています．おや？原点が(from[0],from[1],from[2])だけ平行移動しているな．
凸分解(Convex Decomposition)
まあ何をするのにも凸集合が素晴らしいですよね〜. というわけでぇ，与えられた集合Ωを互いに素な凸集合Ωiの合併で表すのが吉なのよ.
#include <CGAL/Nef_3/SNC_indexed_items.h>
#include <CGAL/convex_decomposition_3.h>
...
using Nef_Polyhedron=CGAL::Nef_polyhedron_3<Kernel>;
using Indexed_Polyhedrons=CGAL::Nef_polyhedron_3<Kernel, CGAL::SNC_indexed_items>;
using Polyhedrons_iterator=Indexed_polyhedrons::Volume_const_iterator;
...
Nef_Polyhedron A;
...
std::vector<Polyhedron> convex_parts;
if ((!A.is_empty())&&(A.is_simple())&&(A.is_bounded())) {
Indexed_polyhedrons N=A;
CGAL::convex_decomposition_3(N);
Polyhedrons_iterator it = ++N.volumes_begin();//最初は領域外部を表すので省略
for(;it != N.volumes_end();++it) {
if (it->mark()) {
Polyhedron P;
N.convert_inner_shell_to_polyhedron(it->shells_begin(), P);
convex_parts.push_back(P);
}
}
}
まあこんな感じで凸集合Ωiのstd::vectorをゲット.
今度はだいじょうぶー
交差
下の図のように図形が交差したりしなかったりすることがありますよね. この交差部分の形状を取得したり, その体積を計算するのはNEF図形ではかーんたんです(ノーマル図形は立体と平面とかしかできねえ・・・ので無視無視〜).
blenderの画面暗いわ．
代表的な論理演算は
| A-B | Aに含まれBに含まれない要素の集合です． |  |
| A*B | AとBの共通部分 A⋂B |  |
| A+B | AとBの合併 A⋃B | |
| A^B | A⋃B - A⋂B |
計算結果はNEF図形ですが，もちろん簡単にノーマル図形に戻せますし，そしたら三角錐分解や体積を求めるのは上と同じで簡単ですね．