Clonon's Web site

Overview

前回は環境を整えたので、今回からは色々プログラミングしていく。

プロジェクト準備

まずソリューション/プロジェクトの準備。
適当なmain.cppを作って、空のフォルダに置く。
同じフォルダに、前回のexampleの各フォルダにあったCMakeLists.txtをコピーして、一部を修正して使いまわす。

# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.

cmake_minimum_required(VERSION 3.1...3.15)
project( Graphics )　======ここを変更。ソリューションの名前になります。

find_package(CGAL QUIET)

//以下略

ビルドディレクトリを作って前回と同じようにcmake。
これでCGALが使える空のソリューションができるので、あとはNuGetでnupenglを入れて準備完了。

プログラミング開始

OpenGLは初めてなので、まずは画面表示やら描写方法から学んだ。
とりあえずなんとか立方体は回転させたりできるようになったので、さっそくやりたいことを初めて見る。

当初のモチベーションとしては、三次元の点群を基にSurfaceを再構成して、ポリゴンを描写してみたいと思っていた。
ルービックキューブのモデリングとして、平面/球面で立体を分割し、パーツごとに移動させたいと考えたから。
例えば分割した各パーツに含まれる微小体積要素をリストアップ、表面にある微小体積要素を基に表面を再構成してパーツを構成。
だけれども、CGALを色々調べてたら3Dポリゴン同士のintersectionも普通に計算できるみたい。
そりゃそうか、出来るわな・・・
ということで必要はなくなりそうだけど、やりたいのでやってみる。

プログラミング開始、再構成まで

まずは点群を生成する。
一番わかりやすい球をつくってみる。
ランダムな球面上の点を乱数で作る。
極座標を単純に乱数で生成すると分布が偏ってしまうので、arccosを使ってz方向の分布を調整し一様にする。

for (int i = 0; i < N; ++i) {
    // dist(mt) : 一様分布乱数 [0,1)
    double t = 2 * M_PI * dist(mt);      // theta
    double p = acos(2 * dist(mt) - 1);   // phi

    double x = sin(p) * cos(t);
    double y = sin(p) * sin(t);
    double z = cos(p);
}

これをCGALに渡してやるわけだが、渡すために型を調整する。
CGALのほうで定義されているPoint_3という型をvectorに詰めて渡す。

typedef CGAL::Exact_predicates_inexact_constructions_kernel::Point_3 Point;
vector<Point> points(N);

for (int i = 0; i < N; ++i) {
    //同上
    points[i] = Point(x, y, z);
}

さて、この点群を再構成するわけで、今回は「Scale-Space Surface Reconstruction」なるものを使ってみた。マニュアルはここ→ https://doc.cgal.org/latest/Scale_space_reconstruction_3/index.html
ここのページは結構親切に書いてくれて助かった。
これによると、ReconstructionオブジェクトにPointのvectorをinsertで突っ込んで関数を呼んでやればいいらしい。
Reconstructionオブジェクトのコンストラクタでvectorのbeginとendを渡すこともできるけど、insert関数を使うのと変わらない模様。 increase_scaleなる関数ではiterationの回数?を指定している。増やすとよりスムーズな表面になるよう複数回iterationするが当然時間はかかるみたい。

typedef CGAL::Scale_space_surface_reconstruction_3 Reconstruction;
typedef Reconstruction::Facet_const_iterator Facet_iterator;

// Construct the mesh in a scale space.
reconstruct.insert(points.begin(), points.end());
reconstruct.increase_scale(1);
reconstruct.reconstruct_surface();

これで再構成完了。お見事。自分でこういうアルゴリズム実装するのはなかなか出来なさそう・・・

描写まで

再構成完了したら、再構成結果を指すiteratorを関数で取得できる。
ここで貰えるiteratorが指すのは、三角形ポリゴンの三頂点を表すインデックス、つまりは点群として渡したvector<Point>での頂点の番号。
型としてはarray<unsigned int, 3>。
これで点群の座標を基に三角形を書いていける。
OpenGLでの光の反射のために、法線ベクトルも計算しておく。
ここでちょっと気づかず困ったのが、CGALで2つのベクトルから法線ベクトルを計算してくれる

CGAL::normal(Point p, Point q, Point r)

の返り値が単位ベクトルでないこと。そういうものか。
これに気付かず、しばらく光の描写がうまくいかなかった・・・
おそらくベクトルの外積をそのまま返してきているので、自前で単位ベクトルに変換しよう。

Vector_3 vec = CGAL::normal(p, q, r);
vec = vec / sqrt(vec.squared_length());

ここまでそろえば後は描くだけ。

typedef Reconstruction::Facet_const_iterator Facet_iterator;

glBegin(GL_TRIANGLES);
for (Facet_iterator it = reconstruct.facets_begin(); it != reconstruct.facets_end(); ++it) {
    array ar = *it;
    Point p,q,r //→ar[0]~ar[2]のPoint
    Vector_3 vec = CGAL::normal(p, q, r);
    vec = vec / sqrt(vec.squared_length());
    glNormal3d(vec.x(), vec.y(), vec.z());
    glVertex3d(p.x(), p.y(), p.z());
    glVertex3d(q.x(), q.y(), q.z());
    glVertex3d(r.x(), r.y(), r.z());
}
glEnd();

結果とパフォーマンスについて

全体として非常に使いやすかった。うちの研究のフレームワークもこれくらい分かりやすく/使いやすく作られてたらなあ・・・
以上の結果を以下に載せておく。
下図左が乱数で作った点群、右が再構成したSurface。
無駄に回転させています。