2017-01-17

OpenVDBのビルドメモ

OpenVDBはボリュームデータのためのフォーマットです．開発元からvdbを扱うツール，ライブラリも公開されているので，今回はそれらをビルドする方法を簡単にまとめておきます．
環境は以下の通りです．

Ubuntu 14.04
OpenVDB 4.0.0

準備

以下のリンクからソース(openvdb-4.0.0.tar.gz)を入手します．
Release v4.0.0 · dreamworksanimation/openvdb · GitHub

OpenVDBには多くの依存ライブラリがありますが，Houdiniがインストールされていると，面倒な依存ライブラリのインストールが必要ありません．OpenVDB側としても，Houdiniをインストールして環境を整えるのが推奨のようです．HoudiniにはApprentice版と呼ばれる無料体験版があり，ほぼすべての機能を試すことができます．これにはHDK(Houdini Development Kit)が付属しており，ここにOpenVDBの殆どの依存ライブラリが含まれています．ユーザー登録をすませたら簡単にインストールすることができます．僕はHoudini15.5をインストールしました．

Makefileの編集

ダウンロードしたソースファイルを解凍したら，openvdb-4.0.0/openvdb/Makefileを編集してビルドを行います．Makefileにも詳しいビルド方法が書いてあるので目を通すと良いでしょう．ちなみにgitリポジトリの方ではCMake対応もされたようですが，僕はうまく使えませんでした．

まず，インストールディレクトリを設定します．デフォルトでは/tmp/OpenVDBです．

#line 84
DESTDIR := /usr/local/OpenVDB

Houdiniに含まれていないライブラリにcppunit, log4cplusがあります．恐らくテスト系のライブラリですが，今回は必要なさそうなのでMakefileに従って無効化します．もしも必要であれば，該当ライブラリをインストール後に，INCL_DIRとLIB_DIRを設定してやります．

#line 129
CPPUNIT_INCL_DIR := #rel/map/generic-2013.22/sys_include

#line 136
LOG4CPLUS_INCL_DIR := #/rel/folio/log4cplus/log4cplus-1.0.3-latest/sys_include

また，openvdb_viewerが必要な場合はGLFWが必要となります．今回は必要ないのでこれも無効化します．

#line 143
GLFW_INCL_DIR := #/rel/third_party/glfw/glfw-3.0.1/include

python module

python moduleを使うにはBoost.pythonが必要です．また，numpyが必要な機能を使う場合にはnumpyもインストールする必要があります．libboost1.55-devとnumpyをaptでインストールし，ヘッダファイルとライブラリへのパスを設定してやります．Houdiniに含まれるPythonは2.6みたいですが，libboostやnumpyは2.7のものでも大丈夫なようです．
以下は僕の環境における設定です．

#line 162
BOOST_PYTHON_LIB_DIR := /usr/lib/x86_64-linux-gnu/
BOOST_PYTHON_LIB := -lboost_python-py27
# The directory containing arrayobject.h
# (leave blank if NumPy is unavailable)
NUMPY_INCL_DIR := /usr/lib/python2.7/dist-packages/numpy/core/include/numpy

もしもpythonモジュールが必要なければ，以下の部分をコメントアウトします．

PYTHON_VERSION := #2.6

その他

epydocやdoxygenがなければこれも無効化しておきます．

#line 169
EPYDOC := #/rel/map/generic_default-2014.24.237/bin/epydoc
#line 180
DOXYGEN := #doxygen

ビルド

Houdiniがインストールされているフォルダへ移動し，houdini_setupをsourceで読み込むと，HDK関連のパスが設定されます．その後，OpenVDBのフォルダへ移動し，makeを走らせます．

cd path/to/houdini
source houdini_setup
cd path/to/openvdb
make -j4
sudo make install

うまくビルドされれば，vdb_printとvdb_renderというバイナリとライブラリができているはずです．
以下のリンクからサンプルの.vdbファイルを入手できます．
OpenVDB - Download
以下の画像はbunny_cloud.vdbをvdb_renderでレンダリングしたものです．
f:id:soma_arc:20170117155917j:plain
ライブラリの使い方に関しては以下のページで簡単なサンプルを見ることができます．

Houdiniをインストールしない方法

もちろん，HoudiniをインストールせずにOpenVDBをビルドすることもできます．libboost（1.55あたりで大丈夫だと思います）, OpenEXR, TBB, python（pythonモジュールが必要であれば）を用意し，ヘッダファイルとライブラリへのパスを設定します．あとはCONCURRENT_MALLOC_LIBをtbbを使うように設定すればビルドできるはずです．

#CONCURRENT_MALLOC_LIB := -ljemalloc
CONCURRENT_MALLOC_LIB := -ltbbmalloc_proxy -ltbbmalloc

一応全部aptで揃うみたいなので，面倒なのはパスの設定ですね．

以上となります．

2017-01-04

blenderでmantaflowで作ったメッシュを表示する

mantaflowは流体シミュレーションのフレームワークです．計算された流体のメッシュをblenderで読み込み，表示する方法のメモです．blenderにも流体シミュレーションを行い，キャッシュを保存しておくbake機能はありますが，他の方法で生成されたメッシュを読み込むきちんとした機能はないようです．しかし，bakeされたファイルのフォーマットに合わせてファイル名を付けてやると，表示することができます．
環境は以下の通りです．

windows10
mantaflow 0.7 prebuilt binary
blender 2.78

mantaflowのインストール

windowsの場合，ビルド済みのバイナリが配布されています．各種dllも同梱されていますが，シーンファイルを読む際にpythonのライブラリがないとエラーが出るみたいです．そのため，python2.7系をインストールしておきます．mantaflowを起動するとシーンファイルの場所を聞かれるので，同梱されている適当なシーンファイルを選び，GUIが表示されれば成功です．

メッシュを作る

同梱のfreesurface.pyからメッシュを作ってみます．メッシュを保存するためには，シーンファイルを少しいじる必要があります．freesurface.pyの62行目あたりにあるmesh.save('phi%04d.bobj.gz' % t)のコメントアウトを外し，ファイル名を'fluidsurface_final_%04d.bobj.gz'に書き換えます．結果としてファイルは以下のようになります．

if (dim==3):
	phi.createMesh(mesh)
	mesh.save('fluidsurface_final_%04d.bobj.gz' % t)

そのあとmantaflowでfreesurface.pyを読み込むと計算されたメッシュが連番の.bobj.gz形式で保存されていきます．ちなみにbobjはバイナリ形式のobjファイルみたいです．
f:id:soma_arc:20170104121203p:plain

blenderで読み込む

blenderを開いたら右側のタブの右端にあるphysicsを選択します．

その他の項目を以下のように設定します．
Enable physics for: -> Fluid
Type: -> Domain
Viewport Display: -> Final
そして，ディレクトリ設定で先ほど保存したメッシュのあるディレクトリを指定すると以下のように読み込まれたメッシュが表示されます．

座標系の違いから横向きになってしまうようです．下部のシークバーでアニメーションできます．
f:id:soma_arc:20170104120438g:plain
以上となります．

References

以下のポストを参考にしました．ありがとうございます．

2016-12-25

GLSLで描くCircle Inversion Fractals

qiita.com
GPUで暖を取りたい人のためのGLSL Advent Calendar 2016，25日目の記事です．

前回の記事で紹介した円の反転という操作を用いてフラクタルを描いてみます．
f:id:soma_arc:20161220183920p:plain
まずはこの画像のような４つの円板とその反転を考えます．それぞれの円の反転は自分以外の円板を自分の内側に移します．よって，それぞれの円の内側には３つの円が移され，１２個の円ができています．次に，ショットキー円の反転を新たに現れた小円にも適用することを繰り替えすと，以下のように円が入れ子状に続いていきます．
f:id:soma_arc:20161220184220p:plain
円が接するときには，このような美しいネックレス上の構造を見ることができます．これらの円列をショットキー円の軌道(Orbit)，円列の極限を極限集合と呼びます．
f:id:soma_arc:20161220184542p:plain
この処理はJavaScriptで書くと以下のような感じになります．

let orbits = [];
orbits.push(schottkyCircles); // level 0
for(let level = 1 ; level < MAX_LEVEL ; level++){
    orbits.push([]);
    for(let i = 0 ; i < orbits[level - 1].length ; i++){
        let orbitCircle = orbits[level-1][i];
        for(let j = 0 ; j < schottkyCircles.length ; j++){
            if(schottkyCircles[j].contains(orbitCircle)) continue;
            let nc = schottkyCircles[j].invert(orbitCircle);
            orbits[level].push(nc);
        }
    }
}
return orbits;

CPUで実装すると，変換の木構造を幅優先探索するようなコードになり，円の深さが深くなるに連れて指数オーダーで計算量が増えてしまいます．GLSLをつかったアプローチを用いると，非常に効率よく描くことができます．

アルゴリズム

全ての点に対して，その点がいずれかのショットキー円の内部に存在するときに，点を属しているショットキー円による反転で移します．この操作を変換された点がすべてのショットキー円の外側に落ちるまで続けます．最後に，反転の回数が円の深さになっているので，深さに応じて色付けを行います．このアルゴリズムはIterated Inversion System(IIS)と呼んでいます．
f:id:soma_arc:20161223174333p:plain
このアルゴリズムのコードを以下に示します．

const int ITERATIONS = 30;
int IIS(vec2 pos){
    int loopNum = 0;
    bool loopEnd = true;
    for(int i = 0 ; i < ITERATIONS ; i++){
        loopEnd = true;
        if(distance(pos, c1Pos) < c1Radius){
            pos = circleInvert(pos, c1Pos, c1Radius);
            loopEnd = false;
            loopNum++;
        }else if(distance(pos, c2Pos) < c2Radius){
            pos = circleInvert(pos, c2Pos, c2Radius);
            loopEnd = false;
            loopNum++;
        }else if(distance(pos, c3Pos) < c3Radius){
            pos = circleInvert(pos, c3Pos, c3Radius);
            loopEnd = false;
            loopNum++;
        }else if(distance(pos, c4Pos) < c4Radius){
            pos = circleInvert(pos, c4Pos, c4Radius);
            loopEnd = false;
            loopNum++;
        }
        if(loopEnd) break;
    }
    return loopNum;
}

お気づきかもしれませんが．よく見ると極限集合の内側と外側は黒い基本領域をショットキー円で変換して敷き詰めた形状としても見ることができます．数学的な理由によって，Hyperbolic Tessellationのように極限で内側と外側が分けられるため，最終的に落ちた点が内側か外側かを判定し，処理を切り替えることで，極限のエッジをレンダリングすることができます．
f:id:soma_arc:20161223174902p:plain
f:id:soma_arc:20161223174846p:plain
マンデルブロ集合に似たフラクタル状のエッジが綺麗ですね．
Kissing-schottky Limit set

Bitmap Orbit trapと組み合わせる

Bitmap Orbit trapという手法を用いるとテクスチャの軌道も描くことができます．ループの過程で，あらかじめ置いておいたテクスチャの領域に点が入った時に，テクスチャの色を参照してループを脱します．前々回の記事で猫を書いた時とほぼ同じ手法です．
f:id:soma_arc:20160421173008g:plain
このgifアニメは四か所に猫のテクスチャを配置しています．詳細はコードをご覧ください．
Kleinyan Cat
Bitmap Orbit trapは一年位前に書いた記事でもまとめておきました．
Orbit trapを考える - 心鏡曼荼羅
今読み返すと読みにくいので書き直したいところです…

三次元への拡張

このアルゴリズムを三次元に拡張し，球をベースにしたフラクタルを描くこともできます．今回はアドベントカレンダーの枠と時間がなく紹介することができなかったのでまたの機会に記事を書きたいと思います．興味のある方は僕の公開した論文，スライドやShadertoyのコード等を読んでみてください．
f:id:soma_arc:20161223181740p:plain

備考 - Schottky transformationについて

ショットキーという語はショットキー変換(Schottky transformation)から来ています．ショットキー変換は２つの円をペアとする変換で，相手の円板の外側を自分の内側に，自分の円板の外側を相手の円板の内側へと移す変換です．これは，本質的に円の反転と同じものになることがあります．ショットキー変換を構成する円のことをショットキー円と呼ぶのですが，奇数個の円を扱うことを考えると，厳密にはショットキー円とは呼べません．詳しくは，後述のインドラの真珠を読んでみてください．

おわりに

インタラクティブに円板を配置してフラクタルを描画することのできるWebアプリケーションを開発しています．また，単純な円盤以外の変換や，三次元のものも実装しています．
SchottkyLink
f:id:soma_arc:20161222215658p:plain f:id:soma_arc:20161223181515g:plain
今回のGLSLアドベントカレンダーでは三回にわたって，Tessellationに関連するレンダリング手法を紹介しました．何かの参考になれば幸いです．それでは，来年の皆様のGLSLライフがより良きものになりますように．

2016-12-24

GLSLで描くHyperbolic Tessellation

qiita.com
GPUで暖を取りたい人のためのGLSL Advent Calendar 2016，24日目の記事です．

前回の記事で紹介したアルゴリズムを用いて双曲円盤上での敷き詰め模様を描いてみます．前回の例では，ユークリッド平面，つまり私達が普段考える平面での敷き詰めを行いました．敷き詰め模様は無限に続き，果てが見えません．
幾何学を研究する数学者たちは，私達の世界の法則を崩した世界を考えました．双曲幾何学(Hyperbolic Geometry)はそんな世界を考える幾何学の分野の一つです．そして，双曲幾何学の世界を一枚の円盤だとみなしたモデルが存在します．このモデルはポアンカレの円板モデルと呼ばれ，円盤の端がユークリッド平面における無限遠点，世界の端となり，円弧がこの世界における直線だと考えられます．また，三角形の内角の和が180度以下になります．この世界におけるタイルの敷き詰めがHyperbolic Tessellationです．
f:id:soma_arc:20161220204320p:plain
半球の内側にある円板がポアンカレの円板です．ユークリッド平面に対する双曲平面です．直線を半径が無限の円板の円周だとみなすと，あらゆるタイルが曲線で構成されていると考えることができます．円板の内側と外側は円板に関する反転という操作で移り合います．見た目は異なりますが，幾何学的，本質的な意味は変わりません．
画家のエッシャーはこのような敷き詰め模様に着想を得て，Circle Limitや天使と悪魔といった有名な作品を生み出しました．双曲幾何学や，Hyperbolic Tessellationに関して，厳密な話はここでは触れません．教科書等を参照してみてください．

基本タイルと変換を見つける

早速前回同様基本タイルを見つけていきましょう．
f:id:soma_arc:20161220163855p:plain
円板は原点中心，半径１の円とします．中央付近の一番大きなタイルがx軸，y軸と２つの円弧に囲まれているのがわかるでしょうか．Hyperbolic Tessellationでは円に関する反転という操作を使って敷き詰めを行います．今回はx軸対称，y軸対称，２つの円に関する反転をとり続けることで，この敷き詰めを描くことができます．

円の位置を特定する

Hyperbolic Tessellationは，タイルの辺の角度が重要になります．どんな図形でも敷き詰めることができるわけではなく，Hyperbolic Tessellationが成立する条件が存在します．これはポアンカレの貼り合わせ定理として知られています．例えば，双曲平面上での三角形である双曲三角形を考える場合は以下の条件を満たす必要があります．
${ \displaystyle ３内角を\frac{\pi}{p},\frac{\pi}{q},\frac{\pi}{r}とした時，\frac{1}{p} + \frac{1}{q} + \frac{1}{r} < 1 }$
以下の画像は一見正しそうに見えますが，正しい敷き詰めではありません．どこかでタイル同士が重なりあったりしています．
f:id:soma_arc:20161222170252p:plain
双曲四角形は２つの双曲三角形を貼り合わせることで得ることができます．今回のタイルは３つの角がPI/2，残りの角がPI/3の４辺形です．よって軸と垂直に交わり，円同士の角度がPI/3となる２円を求める必要があります．円の位置と角度の関係と双曲直線の条件から連立方程式を解きます．双曲直線の条件は以下の通りです．
円の中心を(x, y)半径をrとします．その円弧が双曲直線であるとき以下の条件を満たします．
${ \displaystyle x^2 + y^2 = r^2 + 1 }$
今回の２円はだいたい中心(1.7317, 0),(0, 1.7317)半径はともに1.413となります．円の座標と半径を求める作業が一番面倒なのですが，時間ができたら他のタイルに関しても計算したいところです．

円に関する反転

円に関する反転は，円の中心を無限遠点に，無限遠点を円の中心にもってくる操作です．この変換には等角性や円円対応といった面白い性質があります．ある点Pを中心C，半径Rの円で反転する式は以下のようになります．
${ \displaystyle \frac{(P - C) * R^2}{distance(P, C)^2} + C }$
f:id:soma_arc:20161222151738p:plain
円の外側に描いた図形を黒の円に関する反転によって変換するとこのような図になります．直線は円に，円は円に移ります．変換の前後で図形の角度は変わりませんが(例えば緑の四角形の内角は変換後も直角を保つ)，面積などは大きく変化します．また，同様にして三次元空間で球の反転を定義することもできます．

Hyperbolic Tessellationを描く

中央の一番大きな４つのタイルのうち，右上のタイルを基本タイルとした敷き詰めのコードが以下のようになります．

for(int i = 0 ; i < ITERATIONS ; i++){
    fund = true;
    if (pos.x < 0.){
        pos *= vec2(-1, 1);
        invCount++;
        fund = false;
    }
    if(pos.y < 0.){
        pos *= vec2(1, -1);
        invCount++;
        fund = false;
    }
    if(distance(pos, c1Pos) < r1 ){
        pos = circleInvert(pos, c1Pos, r1);
        invCount++;
        fund = false;
    }
    if(distance(pos, c2Pos) < r2 ){
        pos = circleInvert(pos, cPos2, r2);
        invCount++;
        fund = false;
    }
    if(fund)
        break;
}

f:id:soma_arc:20161220201855p:plain
驚くことに，円板の内側の点は，内側の基本タイル(中央，右上の水色のタイル)に，外側の点は外側の基本タイル(右上の最も大きな紫色のタイル)に移されます．円板の内側のピクセルに対して計算を行うことで円板内の敷き詰めのみを描画することができますが，最終的な点の位置によって色付けを変えることで内側と外側を区別することもできます．
Hyperbolic TessellationはCPUベースで実装すると，計算量が多くなるため，端まで描画するのは大変なのですが，この方法を用いると端まできっちりとリアルタイムで描画することができ，GLSLの計算精度の限界まで拡大することができます．
f:id:soma_arc:20161220201953p:plain
f:id:soma_arc:20161223144704g:plain
敷き詰めの過程はこのような感じです．コードはこちら
Process of the Tessellation

おわりに

今回のGLSLで双曲タイリングを行うWebアプリケーションを作ったので試してみてください．対応しているパターンは２種類しかありませんが，タイリングの変形とWebカメラからの入力を敷き詰める機能を実装してあります．
Hyperbolic Tessellator
Shadertoyにもいくつかコードを上げています．

数学の専門分野に入ってしまったうえ，説明もうまくできなかったのでわかりにくかったかもしれません．しかし，少し専門領域に踏み込むことで面白い題材が見つかるのではないでしょうか．明日の記事では今回触れた円に関する反転を用いたフラクタルについて紹介します．

2016-12-23

GLSLで描くTessellation

qiita.com
GPUで暖を取りたい人のためのGLSL Advent Calendar 2016，23日目の記事です．

今回の記事ではGLSLで敷き詰め模様，Tessellationをレンダリングする方法を紹介します. これを読むような人ならばTessellationというと，ポリゴンの三角形分割を思い浮かべる方が多いかもしれませんがこの記事で話題にするのはより広義のTessellation,平面を平面図形で埋め合わせる平面充填のことです．タイリングとよんだりもします．GLSLでなんらかの絵を描いたことがある人ならば一度はチェッカーパターンを描いたことがあると思います．これも敷き詰め模様です．

正方形の敷き詰め

f:id:soma_arc:20161218181520p:plain
まず初めに，ごく簡単な正方形のtessellationを考えてみましょう．画像のように，正方形で平面全体を敷き詰めようとするとき，基本となる正方形を縦と横に動かすことを考えると思います．縦と横の平行移動，全４種類の変換の組み合わせで平面全体を正方形で埋め尽くすことができそうです．
f:id:soma_arc:20161222134855p:plain
GLSLで描画する場合はその逆を行います．すなわち，各ピクセルを敷き詰めを構成する４種類の変換で基本タイルに向かって動かしていきます．そうして，基本タイルに入ったあとで操作の回数によって色付けを行ったり，基本タイル上の色を参照するなどします．
f:id:soma_arc:20161222134927p:plain
以下が実装例です．Shadertoy上の実装です．

vec3 hsv2rgb(vec3 c){
    vec4 K = vec4(1.0, 2.0 / 3.0, 1.0 / 3.0, 3.0);
    vec3 p = abs(fract(c.xxx + K.xyz) * 6.0 - K.www);
    return c.z * mix(K.xxx, clamp(p - K.xxx, 0.0, 1.0), c.y);
}

const int MAX_ITERATIONS = 100;
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    float ratio = iResolution.x / iResolution.y / 2.0;
    
    vec2 pos = (fragCoord.xy / iResolution.yy ) - vec2(ratio, 0.5);
    pos *= 10.;
    
    bool isFund = true;
    float opCount = 0.;
    float rectSize = 1.;

    for(int i = 0 ; i < MAX_ITERATIONS ; i++){
        isFund = true;
        if(pos.x < 0.){
            pos.x += rectSize;
            opCount++;
            isFund = false;
        }
        if(pos.x > rectSize){
            pos.x -= rectSize;
            opCount++;
            isFund = false;
        }
        if(pos.y < 0.){
            pos.y += rectSize;
            opCount++;
            isFund = false;
        }
        if(pos.y > rectSize){
            pos.y -= rectSize;
            opCount++;
            isFund = false;
        }
        if(isFund) break;
    }   
    fragColor = vec4(hsv2rgb(vec3(opCount * 0.2, 1., 1.)), 1.0 );
}

上の画像は操作の回数によって色をつけています．とても簡単で計算量の少なそうなコードになりました．各ピクセルを変換し続ける処理は，MAX_ITERATIONSによって制限せざるを得ませんが，多くの場合は問題ありません．

moduloによる効率化

多くの方はお気づきだと思いますが，単純な平行移動はmoduloをつかって効率化することができます．

if(pos.x < 0. || rectSize < pos.x){
    opCount += abs(floor(pos.x / rectSize));
    pos.x = mod(pos.x, rectSize);
    isFund = false;
}
if(pos.y < 0. || rectSize < pos.y){
    opCount += abs(floor(pos.y / rectSize));
    pos.y = mod(pos.y, rectSize);
    isFund = false;
}

操作の回数は商の絶対値をとってやれば計算できます．注意が必要なのは，基本タイルの位置が軸に沿っていない時です．正しく商と余りを計算するためには，タイルの辺が軸に沿うように一度平行移動や回転を行う必要があります．

正方形の敷き詰めを行うコードはShadertoyにアップロードしてあります．
https://www.shadertoy.com/view/4t3SDs

鏡映によるTessellation

もう一つ試してみましょう．
f:id:soma_arc:20161218195601p:plain
この模様はどのような基本タイルと変換で敷き詰めているでしょうか．平面は直角二等辺三角形でうめつくされています．また，よく見ると猫が三角形の各辺に対して鏡写しになっているのがわかります．これは三角形の各辺における鏡映変換という操作で移されています．基本としているタイルは(0, 0), (1, 0), (1, 0)を頂点に持つ赤い三角形です．この敷き詰めを描くコードの一部は以下のようになります．
x軸対称，y軸対称，そして斜辺による対称をとり続けることで基本タイルへと点を持っていくことができます．

for(int i = 0 ; i < MAX_ITERATIONS ; i++){
    isFund = true;
        
    if(pos.x < 0.){
        pos.x *= -1.;
        opCount++;
        isFund = false;
    }
    if(pos.y < 0.){
        pos.y *= -1.;
        opCount++;
        isFund = false;
    }
    // the inversion of line (x + y - 1 = 0)
    pos -= vec2(0, 1);
    pos = invRotateM * pos;
    if(pos.x > 0.){
        pos.x *= -1.;
        opCount++;
        isFund = false;
    }
    pos = rotateM * pos;
    pos += vec2(0, 1);
        
    if(isFund) break;
}

今回は，変換の回数で三角形の色を決めるのに加え，基本タイルの上に置いておいた猫のテクスチャを参照しています．完全なコードはこちらです．
https://www.shadertoy.com/view/4l3SDs

まとめ

今回紹介したアルゴリズムをまとめると以下のようになります．

基本タイルを見つける．
基本タイルを敷き詰めるための変換を見つける
各ピクセルを基本タイルに入るまで変換し続ける
変換の回数や，基本タイルの色を使って色をつける

おわりに

GLSLで敷き詰め模様をレンダリングするアルゴリズムを紹介しました．基本がわかれば色々と応用がきくと思います．明日の記事では双曲平面上でのタイルの敷き詰め，Hyperbolic Tessellationについて紹介します．

2016-12-07

数学とアートのカンファレンス，Bridgesの紹介

www.adventar.org
この記事は日曜数学Adbent Calender7日目の記事です．今年の夏に8月9日から8月13日までフィンランドで開催されたBridgesというカンファレンスに参加してきました.今回はこの会議について紹介します．

Bridgesとは

Bridgesは数学とアートやその他の文化に関する国際会議です.
The Bridges Organization - About Bridges
Mathematics, Music, Art, Architecture, Education, Cultureの謳い文句の通り,数学だけではなく,数学とその他の文化的なものとの関わりを育てていくことを目的としています.参加者のバックグラウンドは数学者,教師,アーティスト,音楽家,パフォーマー,プログラマー等と多岐に渡り，学際的な領域にいる人々が集まります.このような領域の会議としては世界最大のものらしく，今年は350人ほどが参加したようです．
会では論文発表やワークショップの他，常設でアート作品展示がされます．展示された作品は投票によって賞が与えられます．また，Public Dayという一般に開かれた日も設けられており，劇や詩の朗読，ジャグリングなどのパフォーマンスが行われます．
その他にも，空いた場所に机をおいて作品販売を行っていたり，ロビーで音楽演奏，ダンスを始める人がいたりと，自由な雰囲気で，会議というよりは祭典という言葉がふさわしいかもしれません．

Call for Submissions

いくつか募集分野について紹介します．2016年12月6日現在，2017年のサイトはあまり準備が整っていないようなので，詳しく知りたい場合は去年のページを読むと良いです．
The Bridges Organization - Bridges 2016

Call for Papers

論文はRegular paperとShort paperが募集されます．それぞれきちんと査読されます．採択されるとRegular paperは３０分，Short paperは１５分のスロットが与えられてプレゼンを行うことになります．過去の論文は以下のリンクから全て読むことができます．様々なネタの宝庫です．
The Bridges Archive
僕はShort paperを出して参加しました．ある特定のフラクタルを高速に描画するためのアルゴリズムについて記述しました．

Art Exihibition

アート作品の展示です．作品は事前に審査されます．論文とは異なり,会場に来れなくても展示料金を払えば展示させてもらえるようなので,信頼できる人に作品を預けて展示してもらうこともできるらしいです.過去の作品はこちらで見ることができます.
2016 Bridges Conference | Mathematical Art Galleries
2015年に行われた展示の様子がアップロードされていました．
vimeo.com

Short Movie Festival

動画作品の募集です．集められた作品はPublic Dayに上映会が行われます．1~5分が望ましい長さとされていますが，厳密な制限はないようです．
2016 Bridges Conference Short Movie Festival | Mathematical Art Galleries

Student Travel Scholarship

海外で行われるカンファレンスですから，結構お金がかかります．学生参加者には奨学金があります．参加料と開催地までの距離に応じて最大1000$までが給付されます．給付の条件は何らかの論文や作品の提出,もしくは学生ボランティアに参加を申し出ることです.今回は僕も奨学金を出していただきました．

終わりに

数学とアートの国際カンファレンスであるBridgesについて紹介しました．来年はカナダのWaterlooで7月27日から31日まで開催されます．普段の活動をこういった学術的，国際的な場所で発表するのも良いのではないでしょうか．それぞれの分野の募集締め切りは2月初めから5月半ばの間に設定されているので確認してみてください．
www.youtube.com
僕も参加のために論文を書いている最中です．Art ExihibitionやShort Movie Festivalにも出したいところですが，あまり目処はたっていません…
それでは,来年の皆様の活動が実りの多きものになりますように.

日曜数学Advent Calendar8日目の記事はToshiki Takahashiさんによる，"ポリゴングラフと経済学"だそうです．

2016-06-04

点列をCatmull Rom曲線を通してBezier曲線列に変換する

前回，レイトレーシングでBezier曲線を描画する方法を紹介しました．しかし，Bezier曲線そのままでは意図した図形を描画することは大変です．今回は，特定の点列を通るCatmull Rom曲線という曲線をBezier曲線に変換する方法を紹介します．
当初，B-spline曲線からBezier曲線に変換する，Boehm's algorithmを使用しようかと考えたのですが，設定等が複雑なことと，点列を必ず通る精密な図形を描きたかったため，点列を必ず通るという特性を持つ．Catmull Rom曲線を使うことにしました．

scheme(Gauche)による実装はこちら
scheme-raytrace/points.scm at master · soma-arc/scheme-raytrace · GitHub
Processingによる実装も公開している人がいました．
Mike's Processing code snippet repository
この方法と，前回紹介したBezier曲線をレイトレする手法を用いて描いた曲線がこちらです．
f:id:soma_arc:20160526164840p:plain

各種曲線の概要はこのサイトが詳しいです．
デジタル・フロンティア-Digital Frontier | DF TALK | スプライン曲線の話
また，変換を理解するにあたって，このサイトが非常に勉強になりました．
t-pot 『3次曲線』
これに従い，順に曲線を見ていくことで，自然と変換の方法がわかります．

Ferguson/Coons曲線

f:id:soma_arc:20160604193159j:plain
基本となる曲線です．三次曲線を表すため，始点と終点，初速度，最終速度の4種のパラメータを用います．これを用いた曲線の計算式はt-potにすべて載っていますが，ここではあまり重要でないので，記述しません．Bezier曲線もCatmull Rom曲線も，この曲線を元に速度の求め方を考えていきます．

Bezier曲線

f:id:soma_arc:20160604193158j:plain
Ferguson/Coons曲線において，速度を直接決めるのは割と面倒ですので，指定された制御点から，速度を求めてやるようにします．これが，三次のBezier曲線です． ${P_0,P_3}$ を始点，終点としたとき，速度 ${v_0, v_1}$ を以下のようにして求めます．
${v_0 = 3(P_1 - P_0)}$
${v_1 = 3(P_3 - P_2)}$
３という係数がBezier曲線の次数ですね．速度が求められれば，Ferguson/Coons曲線の式が使えます．

Catmull Rom曲線

f:id:soma_arc:20160604193201j:plain
Catmull Rom曲線では，点列すべてを通る曲線を定めることができます．点列のセグメントの前後合わせて4つの点を合わせてパラメータを決めます．ある点における速度を前後の移動量を足し合わせたものに ${tightness}$ というパラメータをかけたものとします．この値は通常，0.5が使われるため，基本的には，移動量の平均値をとることになります．
${v_0 = tightness( (P_1 - P_0)+(P_2 - P_1)) = tightness(P_2 - P_0)}$
${v_1 = tightness(P_3 - P_1)}$
ここで注意が必要なのが，点列における始点，終点です．始点，終点では速度を求めるために必要な点がないため，違う速度を使ってあげる必要があります．t-potでは，最初と最後の速度の平均速度を使うとしています．僕は面倒だったので，実装する際には，あらかじめ点を挿入しておいてある，と仮定して無視しました．
${tightness}$ パラメータに関しては，このpdfで知りました．
https://www.cs.cmu.edu/~462/projects/assn2/assn2/catmullRom.pdf

点列からBezier曲線の制御点を求める

さて，Catmull Rom曲線と同じ挙動をとるBezier曲線の制御点を求めてみましょう．Bezier曲線とCarmull Rom曲線の ${v_0, v_1}$ は等しくなることを利用して，方程式を解きます．ベジェ曲線の制御点を ${BP_0, BP_1}$ とすると点列の四点， ${P_0, P_1, P_2, P_3}$ を用いて以下のように求まります．
${BP_0 = \frac{tightness}{3}(P_2 - P_0) + P_1}$
${BP_1 = \frac{tightness}{3}(P_1 - P_3) + P_2}$
先に述べたように， ${tightness}$ は普通0.5です．これで，点列のあるセグメントに関するBezier曲線の制御点を求めることができました．これを点列全体に関して行えば，Bezier曲線列を得ることができます．

おわり

点列をCatmull Rom曲線を用いてBezier曲線に変換する方法を紹介しました．この方法と前回の記事の内容をあわせることで，点列で表現される図形をレイトレで描画することができます．しかし，点列が多い図形だと曲線の数が増えて，どうしても時間がかかってしまいます．BVHなどを用いて効率よくレイとの交差計算を行わなければ，とても使えたものではないでしょう．
Bezier曲線といえば，内分点をとっていく導出しか知りませんでしたが，今回見たような意味づけは知らなかったので勉強になりました．

準備

Makefileの編集

python module

その他

ビルド

Houdiniをインストールしない方法

mantaflowのインストール

メッシュを作る

blenderで読み込む

References

Bitmap Orbit trapと組み合わせる

三次元への拡張

備考 - Schottky transformationについて

Further Reading

おわりに

基本タイルと変換を見つける

円の位置を特定する

円に関する反転

Hyperbolic Tessellationを描く

Further Reading

おわりに

正方形の敷き詰め

moduloによる効率化

鏡映によるTessellation

まとめ

Further Reading

おわりに

Bridgesとは

Call for Submissions

Call for Papers

Art Exihibition

Short Movie Festival

Student Travel Scholarship

終わりに

Ferguson/Coons曲線

Bezier曲線

Catmull Rom曲線

点列からBezier曲線の制御点を求める

おわり