三次元フラクタルを3Dプリントする

数学とコンピュータ Advent Calendar 2017 - Qiita
数学とコンピュータアドベントカレンダー23日目の記事です.(30分ほど遅れました.申し訳ありません.)

近年,3Dプリント技術の広まりとともに数学的なオブジェクトを3Dプリントによって造形する活動が広まっています.最近ではVisualizing Mathematics with 3D Printingという本も出ました.三次元フラクタルの3DプリントにはJeremie Brunet氏という先駆者がいます.とてもクオリティの高い作品を数多く公開しており,いつか自分でも試して見たいと思っていました.
www.shapeways.comThe Beauty of Math! These 3D Printed Fractals Will Blow Your Mind | 3DPrint.com | The Voice of 3D Printing / Additive Manufacturing

試行錯誤の結果,最近ようやくそれなりのものを出力できるようになりました.Twitterに投稿したものをまとめたのがこちらになります.
twitter.com


それぞれのオブジェクトの造形には10~20時間ほどの時間がかかりました.これらの3DプリントデータはDMM 3DプリントやShapewaysで公開する予定です.

このフラクタルはQuasi-Fuchsian 3D Fractalと呼ばれるフラクタルです.これに関してはこちらの記事で少しだけ触れました.
自作フラクタルレンダラとシェーダの取り回しについて - 心鏡曼荼羅
この記事ではこれらのフラクタルを3Dプリントするために必要な事項を簡単にまとめます.

モデリング・描画

数学的なオブジェクトをモデリングをするツールは様々にあります.こちらの記事に色々まとまっています.
3Dプリンター用の多彩なデザインの3Dモデルがつくれる計算幾何学アプリが楽しい!
しかし,三次元フラクタルはこういったソフトで計算,出力することは難しいです.多くの三次元フラクタルオブジェクトは,距離関数(Signed Distance Function)で表現されます.距離関数は与えられた点とオブジェクトの表面までの最短距離を返す関数です.たとえば,球の距離関数は次のようになります.

float distSphere(vec3 p){
    return distance(p, centerOfSphere) - r;
}

この関数は与えられた点から球の表面までの距離を返します.また,内側の点に関しては負の値を返します.
距離関数を用いた描画についてはこちらの記事にまとまっています.
全能感UP! GLSLで進めレイマーチング « demoscene.jp

3Dプリントを行うためには,出力したい物体のメッシュデータが必要なので,距離関数で表現されるフラクタルをボリュームデータ(ボクセルデータ)として出力し,これをメッシュ化していくことになります.三次元フラクタルに限らず,距離関数で定義できる(Sphere tracing. ray marhcingで描画できる)形状ならばこれから紹介する方法でメッシュ化してプリントすることができます.先述のJeremie Brunet氏もmandelbulbというソフトウェアで出力したフラクタルのボリュームデータをメッシュ化してプリントしているようです.

メッシュ化

ボクセルデータをメッシュ化する方法,ツールは様々にありますが,今回はParaviewとOpenVDBという二種類の方法を試しました.

ParaViewによるメッシュ化

Paraviewはオープンソースのデータ分析,可視化ソフトウェアです.ボリュームデータを読みこんで可視化,メッシュ化することができます.ボリュームデータはvtkと呼ばれるフォーマットで扱うため,データをこの形式で出力します.

vtkにはレガシーなASCIIフォーマットと,バイナリデータを扱うxmlフォーマットが存在しています.簡単のためにASCIIフォーマットを扱うことにします.vtkファイルに関してはこちらの記事が詳しいです.
ParaViewでVTKレガシーフォーマットを使う その1 - Qiita

フラクタルでない部分とフラクタル内部に分かれた二値のボリュームデータをvtk出力するプログラム例がこちらになります.(ビルドにはnanortが必要です)
Generate vtk file of fractal. · GitHub
テキストファイルなので,格子数によっては巨大ファイルになってしまいますが,一応はフラクタルをボリュームデータにすることができます.

次にParaviewで可視化,メッシュ化してみます.Paraviewを起動し,vtkファイルを投げ込んだらProperties tabからApplyをクリックします.RepresentationのドロップダウンリストからVolumeを選択するとボリュームデータが表示されます.
f:id:soma_arc:20171222124436p:plain

ボリュームデータをメッシュに変換します.画面上部のContourツールボタン(もしくはメニューのFilters -> Common -> Contour)をクリックします.
f:id:soma_arc:20171222124944p:plain
画面左のPipeline BrowserにContourというデータができます.これをクリックし,PropertiesタブのApplyをクリックしてやるとメッシュ化されます.あらかじめデータが二値であればデータの境界をメッシュにしてくれます.もしもデータが二値でないのならば,Isosurfacesパラメータを操作する必要があります.
f:id:soma_arc:20171222125119p:plain
フラクタルのボリュームデータをメッシュに変換することができました.データの保存はContourを選択した状態でFile-> Save Data... からファイル形式を.stlに指定して出力してやります.

変換はそれなりに高速で結果も綺麗です.しかし,データサイズが大きくなりがちなので,別のモデリングソフト等でデータを削減する必要があります.この例で使用したデータは比較的小さいですが,簡単に数百MBのデータになってしまうので注意が必要です.

Paraviewによるメッシュ化は非常に手軽にできるので,お勧めです.しかし,なるべくならレンダラにメッシュデータ生成機能を組み込みたいところです.そこで,次に試したのがOpenVDBです.

OpenVDBによるメッシュ化

OpenVDBはボリュームデータを効率よく扱うことのできるオープンソースライブラリです.ボリュームデータからメッシュを生成するような機能も備えています.
OpenVDBでは距離関数によるデータをメッシュ化することができます.OpenVDBのデータ表現に関する情報は次の記事に詳しく書いてあります.
Level sets with OpenVDB. Quick introduction. Part 1 - K. Lykov Blog
データ生成のコード片を置いておきます.メッシュ化処理の呼び出し等を含んだ実際のプログラムはまだ整理できていないため,後で公開します.
genFractalLevelSet.cpp · GitHub

f:id:soma_arc:20171223213220j:plain
Paraviewと比べて遜色ないメッシュデータを作成することができました.また, ボリュームデータフォーマットのvdbファイルはそれなりに小さいので,Houdiniなどの別のレンダラに持って行ってレンダリングするのも面白いかもしれません.Adaptive meshingやsmooth level setなど,試せていない機能もあるので今後の課題です.

OpenVDBの例は少ないですが,このあたりが参考になるかと思います.

メッシュ削減(Mesh Decimation)

生成したメッシュデータ(.obj or .stl)のサイズが大きい場合はメッシュを削減してやる必要があります.多くのモデリングソフトにこのような機能が搭載されています.だいたい頂点数が100万個あたりになるまで減らしています.自分でプリントする分にはデータサイズが大きめでも問題ありませんが,外部のプリントサービスを利用する場合は制限があるので注意が必要です.例えば,DMM3Dプリントだと,100MBを超えるデータはアップロードすることができません

MeshLabによるメッシュ削減

MeshLabはオープンソースのメッシュ処理ソフトです.Jérémie Brunet氏もMeshLabを使用してメッシュの最適化を行っているようです.しかし,今回は処理に時間がかかりすぎて使用を断念しました.データが悪かったのか,処理のパラメータが悪かったのかはわかりませんが,使いこなせればメッシュのクオリティアップに繋がるかもしれません.
MeshLabを使ったメッシュ削減の方法はshapewaysに記事がありました.
Polygon Reduction with Meshlab - Shapeways

Blenderによるメッシュ削減

Blenderオープンソースモデリングソフトです.ModifierのDecimationでメッシュを減らすことができます.

手順は次のようになります.Blenderでメッシュを読み込んだら,そのメッシュを選択し,ModifierタブからAdd Modifier -> Decimationを選択します.
f:id:soma_arc:20171223221751p:plain
0~1の範囲でメッシュを削減する割合を指定します.一度に大きく減らしすぎると,処理に時間がかかって応答しなくなるので,0.2くらいずつ減らしていくと良いかと思います.三種類処理の手法があるようですが,あまり違いはわかっていません.とりあえずCollapseでよいかと思います.
f:id:soma_arc:20171223221818p:plain
ほどほどの数のメッシュになったらApplyで処理を確定し,File->Exportから修正したデータをエクスポートします.ここまでくると,このデータを3Dプリントすることができます.

ZBrushを用いたメッシュ削減

ZBrushはスカルプトモデリングツールと呼ばれる種類のモデリングツールです.粘土をこねるような形でモデリングを行っていくツールで,大量の頂点を効率よく扱うことができます.有償のソフトですが,メッシュの最適化に関する優れた機能をもっています.Blenderでも機能的には十分ですが,やはり時間がかかるのと,後述するUV展開はうまく行えなかったので,ZBrushを使ってみました.ZBrushはなかなかクセのあるソフトですので,作業の流れを紹介するだけにとどめておきます.

ZBrushではDecimation MasterとDynameshという機能を使います.Decimation Masterを使うことで,メッシュの細部を保ったまま,高速にメッシュを減らす事ができます.もしも処理の途中でエラーが発生したらDynameshを使って,メッシュ全体を均一なメッシュに直します.
このあたりの手順はこちらの記事にまとまっています.
デジタル造形の時代〜基本的なワークフローについてわかりやすく解説〜 | 特集 | CGWORLD.jp

ZBrushにも3Dプリント用にデータをエクスポートする機能がついているので,削減した後はデータを出力してプリントできます.

テクスチャデータ生成の試み

ここまでで,3Dプリントデータ生成を見てきました.生成したデータは家庭用(?)の3Dプリンタでプリントする他に3Dプリントサービスを利用してプリントしてもらうことも出来ます.3Dプリントサービスで提供されているプリント素材にはフルカラーのものも存在しています.DMM3Dプリントではフルカラー石膏とフルカラープラスチックという素材を提供しています.カラーリングもフラクタルの見た目において,重要な要素です.フルカラーの3Dプリントデータの作成にも挑戦してみました.

DMM3Dプリントのカラープリントはモデルデータ(.obj)に紐ついた.mtlとテクスチャデータが必要です.そのため,メッシュデータにUV座標をつけ,テクスチャデータを用意する必要があります.

UV展開

まず,三次元空間上の頂点をテクスチャ平面の座標にマッピングする作業である,UV展開を行う必要があります.当初Blenderで試したのですが,メッシュが複雑でデータが大きいためか,うまく行きませんでした.色々試してみましたが,最終的にZBrushで成功しました.

ZBrushではUVMasterという機能でUV展開を行う事ができます.
【UnityAction & ZBrush】UV Masterを使ってみる。: Karasuのアプリ奮闘記
UV展開前のDecimation Masterによるメッシュ削減処理は,UV展開に不適切なメッシュを生成してしまうようです.そのため,Decimation Masterでメッシュをある程度削減した後にZRemesherという機能を使います.ZRemesherを使用すると,UV展開に適した綺麗なポリゴンに変換することができます.

リメッシュ周りの資料です.

テクスチャの生成

メッシュデータにUV座標をつけることができれば,テクスチャデータを生成することは難しくありません.重心座標(Barycentric Coordinates)を計算することで,UV平面の三角形上の点と空間の三角形上の点を対応付けることができます.

愚直な実装ですが,テクスチャを生成するコードを置いておきます.シェーダを使ってレンダリングすればよりクオリティの高いテクスチャを生成できるので,次の課題です.
3d-printing-tools/GenTex at master · soma-arc/3d-printing-tools · GitHub
生成されたテクスチャとマッピングされたオブジェクトは以下のようになります.
f:id:soma_arc:20171223184905p:plain
UV展開の過程でメッシュを削りすぎてしまいましたが,うまくマッピングできました.すぐに3Dプリントサービスに出したいところですが,このままですと,かなり料金がかかります.この状態で見積もりをとった時は十万円を超えてしまいました.かかる料金を減らすため,もうひと手間必要です.

中空化

何らかの3Dプリントサービスに出す場合は,モデルの中身をくり抜いて中空構造にする,肉抜きを行うと,かかる料金を減らすことができます.
【3Dプリンタ】肉抜きをして制作コストを下げよう! | 初心者魂 今更聞けない使い方
BlenderZBrushでも肉抜きを行うことはできますが,meshmixerというソフトウェアを使えば,ほとんど自動で肉抜きを行う事ができます.meshmixerはautodeskが提供するフリーソフトです.

meshmixerを起動し,メッシュを読み込ませたら,編集->中空を選択することで自動的に内側がくり抜かれます.
f:id:soma_arc:20171223173025p:plain
同時に,素材抜き用の穴も生成されます.プリント素材によって素材抜き用の穴の大きさが異なるので注意が必要です.例えば,DMMで石膏素材をプリントする場合,10mm以上の穴が必要になります.

肉抜きしてエクスポートした.objファイルのメッシュグループはバラバラになってしまうようです.
f:id:soma_arc:20171223173950p:plain
利用するサービスによってはひとつにまとめておく必要があるかもしれません.Blenderでは各グループを選択してCtrl + Jで結合することができます.

最終的にDMMで見積もりを取った結果がこちらです.
f:id:soma_arc:20171223175402j:plain
石膏フルカラーで6.7cm x 7.7cm x 5cmで約一万四千円,フルカラープラスチックなら約二万五千円です.流石にこのクオリティでこの値段はしんどいので,メッシュ削減,UV展開の最適化でクオリティを上げ,肉抜きを工夫して値段を減らそうと考えています.

おわりに

ボリュームデータから3Dプリントデータを作成するためのツール群を簡単にまとめました.ボリュームデータさえ用意出来てしまえば,既存のツールの組み合わせで3Dプリントまで持っていくことが可能ですが,根気が必要です.なるべく少ないツール,手間でプリント用データ生成を行えるスキームを確立したいところです.レンダラでレンダリングしたフラクタルの好きな部分を切り取って手軽に3Dプリント用データを生成できるようにする予定です.

かなりの手間がかかる3Dプリントですが,いくつかの場所でデモをした際には,それなりに良い反応をもらえました.やはり,手で触ることのできる立体があることは,アイデアを伝えるために非常に有用であることがわかりました.何か立体化できるようなコンテンツがある方は挑戦してみてはいかがでしょうか.

自作フラクタルレンダラとシェーダの取り回しについて

WebGL Advent Calendar 2017 - Qiita
WebGL Advent Calendar 2017,19日目の記事です.

WebGLを用いたフラクタルレンダラを開発しています.そのコンセプトと中身を簡単に紹介します.今回紹介するのは二種類です.どちらも開発途上で実験段階ではありますが,それなりに面白い図を得ることが出来ます.

フラクタルレンダラ

どちらもWebGL2を要求します.また,新しめのGPUでないと動かないかもしれません.

SchottkyLink

f:id:soma_arc:20171216175254p:plain
f:id:soma_arc:20171216175713p:plain
円や球を元にしたフラクタルを描画することができます.これらの図はリアルタイムに描画され,生成される図を見ながら図形を配置していくことができます.

例えば,下図,左側の図のように,円を配置すると右図のような模様が描かれます.

三次元の場合も同様に,球を配置することでフラクタルが生成されます.
f:id:soma_arc:20171217164052p:plain

しかし,より複雑な絵を描画するためにはアルゴリズム的に困難な部分がでてきてしまい,現在は開発が止まっています.
高速に描画する話や,これを描画するための数学的な話はこちらの論文や記事に書きました.

Experimental Sphairahedron-based Fractal Renderer

f:id:soma_arc:20171216180453p:plain
最近開発を進めているものです.自由度は少ないですが,非常に面白い形状を描画することができます.画像の上部にある二つのパネルのうち,左側に描画されている赤と緑の球によって削られている立体は球面体と呼ばれる立体です.この立体をベースにして右側のパネルに描画されているフラクタルを生成することが出来ます.画像下部の左から二番目のパネルに描画されている赤い点を操作することで,このフラクタルを変形することが出来ます.

このフラクタルの構成方法は以下の動画をみるとなんとなくわかるかもしれません.
www.youtube.com


残念ながら,このフラクタルの描き方に関する文献はほとんどありません.これから論文等を書いていく予定です.また,より自由度を高いパラメータ化も考える予定です.

コンセプト

フラクタル図形の描画というと,マンデルブロ集合をご存知の方は多いかと思います.フラクタルアート界隈ではマンデルブロ集合のようなフラクタルの計算方法*1を用いて三次元のフラクタル図形を描画する式が開発されています.これらのフラクタルのいくつかはFractal Lab - Interactive WebGL Fractal Explorerレンダリングすることができます.
f:id:soma_arc:20171217183255p:plain
これらのフラクタル形状を自分で描画することは,決まった式があるので,さほど難しくはありません.しかし,その式をうまく改変したり,パラメータを操作して自分の思うような形状を得ることは難しいです.

そこで,開発中のフラクタルレンダラでは円や球,直線,平面といった幾何学形状によって定義される操作を組み合わせることでフラクタル形状を生成するアプローチをとっています.パラメータの操作は,図形の位置や大きさを変化させることで行い,レンダリング結果はリアルタイムで表示します.そうすることである程度複雑なフラクタル形状を直観的に生成することができるようになりました.

最終的にはこれらのフラクタルに関わる数学的な性質を理解し,パラメータを直観的に操作できるようなレンダラを作ることで,自由自在にフラクタル形状を生成することを目指しています.

設計

ここで紹介したレンダラでは,GLSLを用いて図を描画しています.GLSLのフラグメントシェーダを用いてレンダリングを行う際には,シーン毎にシェーダファイルを用意しておく必要があります.例えば,上述のFractal Labでも,各フラクタルのシェーダファイルが予め用意してあり,パラメータがUniform変数で渡されて形状が制御されます.しかし,SchottkyLinkではユーザの操作によって様々なフラクタルのジェネレータが追加されるため,あらかじめシェーダファイルを用意しておくことは難しいです.

そこで,ユーザーの操作に応じてシェーダを動的に生成し,シェーダをコンパイルしなおしています.シェーダの生成にはテンプレートエンジンのNunjucksを使用しています.

例えば,下図のような四つの円から構成されるフラクタルを描画することを考えてみます.

この四つの円による反転という操作を繰り返すことで描画されます.この処理のGLSLコードの例は以下のようになります.

for(int i = 0; i < MAX_ITERATIONS; i++) {
    if(distance(pos, u_circle0.center) < u_circle0.radius){
        pos = circleInvert(pos, u_circle0);
        invNum++;
    } else if(distance(pos, u_circle1.center) < u_circle1.radius){
        pos = circleInvert(pos, u_circle1);
        invNum++;
    } else if(distance(pos, u_circle2.center) < u_circle2.radius){
        pos = circleInvert(pos, u_circle2);
        invNum++;
    } else if(distance(pos, u_circle3.center) < u_circle3.radius){
        pos = circleInvert(pos, u_circle3);
        invNum++;
    }
}

Uniform変数で円や球の位置,半径といったパラメータが制御されます.この円の数はユーザの操作によって減ったり増えたりします.そこで,以下のようなテンプレートを作成しておきます.

for(int i = 0; i < MAX_ITERATIONS; i++) {
{% for n in range(0,  numCircle ) %}
    if(distance(pos, u_circle{{ n }}.center) < u_circle{{ n }}.radius){
        pos = circleInvert(pos, u_circle{{ n }}, dr);
        invNum++;
    }
{% endfor %}
}

numCircleパラメータを渡すことで,このテンプレートから円の数に応じた新たなシェーダが生成できます.実際のシェーダはこのような感じになっています.
SchottkyLink/2dCircles.njk.frag at master · soma-arc/SchottkyLink · GitHub
少々汚いですがマクロ等を使えば綺麗に書けるかと思います.シェーダのモジュール化もこのテンプレートエンジンで行っています.

この方法の欠点はデバッグが面倒なことと,一度に複数のシェーダを再生成するとシェーダのリンクあたりに非常に時間がかかってしまうことです.あらかじめ大きめのUniform配列を用意しておく手もありだとは思いますが,まだ試せてはいません.実行速度にどの程度影響があるかが問題です.

おわりに

フラクタルレンダラとそのシェーダの取り回しについて書きました.
GPUへ依存はありますが)このような複雑なレンダラもブラウザで手軽に動かしてデモできるというのは非常に便利です.特に,こういったある意味マニアックなソフトウェアでも触ってもらいやすいという利点もあります.
完成にはほど遠いですが,数学的な部分も含め,粛々と開発を続けていきます.

*1:Escape-time アルゴリズムとも呼ばれ,特定の式を画面上の点に繰り返し作用させ,その点の収束,発散を見る.

数学ソフトウェアデータベース swMATHの紹介

数学とコンピュータⅡ Advent Calendar 2017 - Qiita
この記事は数学とコンピュータアドベントカレンダーⅡ,14日目の記事です.数学ソフトウェアデータベース,swMATHを紹介します.

f:id:soma_arc:20171211191610p:plain
swMATHは数学ソフトウェア(Mathematical Software)のための情報サービスで,様々な数学ソフトウェアの情報をデータベースとして持っています.これらのソフトウェアの情報は数学の論文の引用情報を元に構築されており,ユーザーは自由にデータベースを検索することができます.

ここでの数学ソフトウェアとは,何らかの数学的な問題を解析したり,解いたり,シミュレーションしたりするソフトや,数学の定理やアルゴリズムに基いたソフトウェア,もしくはライブラリのことを指します.例えば,汎用的で大きなソフトウェアだと,MathematicaMatlab等があります.特定の数学分野に由来するものであると,GeoGebra Geometryや,Desmos,SnapPy等が挙げられるでしょうか.また,様々な計算を行うためのライブラリ(例えばPythonのSciPy, NumPy)も数学ソフトウェアとします.

機能

簡単に機能を紹介します.

適当に検索キーワードを入力して検索します.今回はgeometryで検索してみました.
f:id:soma_arc:20171211192215p:plain
すると,被引用文献数でソートされた検索結果がでてきます.SageMathをクリックしてみます.
f:id:soma_arc:20171211195711p:plain

ソフトウェアの概要とキーワードや関連ソフトが表示されます.ここからさらに検索をすすめる事が可能です.ページ中段に表示されている,ORMSとは別の数学ソフトウェアの情報サービスで,人手で信頼性の高いソフトウェア情報をまとめているサービスです.

ページ下部にはこのソフトウェアを論文中で引用している論文と,そのソフトウェアそのものに関する論文(standard articles)のリストが表示されます.
f:id:soma_arc:20171211195719p:plain
リストのリンクからzbMATHという数学論文データベースに飛ぶことで,その論文の詳細情報にアクセスすることができます.また,画面右側のチェックボックスを用いることで,リストのフィルタを操作することができます.出版時期やMSC(数学の分野の分類タグ)を用いて細かくフィルタできます.

また,各論文出版時のソフトウェアのウェブサイトを見ることもできます.画面右上のVersions: Infoをクリックすると,論文タイトルの横に当時のサイトのWeb Archiveへのリンクが表示されます.
f:id:soma_arc:20171213215825p:plain
f:id:soma_arc:20171213215808p:plain

以上がswMATHの機能の簡単な紹介です.

コンセプト

近年,数学ソフトウェアが数学の研究や教育の分野で非常に重要になっています.しかし,数学の学術論文の引用欄に数学ソフトが並ぶことがありながら,それらの数学ソフトウェアを広くアーカイブするような活動はなかなか行われて来ませんでした.数少ないアーカイブ活動は,特定の数学領域に閉じていたり,メンテナンスされずに古い情報しか残っていないという状態であったようです.

一方で,数学の学術論文アーカイブ活動は昔から行なわれてきています.電子的な情報サービスで有名なものにはzbMATHMathSciNetというプロジェクトが存在しています.こうしたサービスはたくさんの論文情報をアーカイブしています.これらはオープンアクセスではありませんが,レビュワーによる要約,評論もつけられています.

学術論文のアーカイブに対してソフトウェアのアーカイブは難しい問題です.例えば,数学の結果は一度論文として出版されれば,それは永久に残りつづけます.(もちろん間違いが含まれることはありますが)しかし,ソフトウェアの開発に終わりが来ることはありません.出版当時に使用されていたソフトウェアも絶えず改良され,変化していきます.逆にメンテナンスされることがなくなってしまえば,将来的には動かなくなってしまいます.

数学の結果はある程度,その主張が正しいか,正しくないかで評価することができます.しかし,ソフトウェアには数学に比べて,たくさんの評価軸が存在します.また,個人の書き捨てのスクリプトから,大きなソフトウェアプロジェクトまで,大小さまざまな数学ソフトウェアがどんどん増えています.

そのような理由から,ある程度の品質をもった数学ソフトウェアのデータベースを人手で継続的に作っていくのは現実的ではありません.そこでswMATHではPublication-based approachと呼ぶ手法を用いてある程度自動的にデータベースを構築しています.

Publication-based Approach

数学ソフトと数学の出版物は密に繋っています.はじめに出版物に記された数学のアイデアがあり,それがソフトウェアに実装されます.実装されたソフトウェアが新たな数学の結果を生み出します.それらの結果は論文として出版されるわけですが,その論文には使用したソフトウェアが引用されます.swMATHではこの引用情報を利用して,自動的にソフトウェアのデータベースを構築しています.学術論文として出版されている論文に引用されているということで,そのソフトウェアの品質はある程度保証することができます.また,被引用件数がそのままソフトウェアの信頼性の高さを表わす事にもなります.

swMATHでは,zbMATHを中心として数学論文データベースに登録されている論文の題名や,引用欄などから自動的にソフトウェア名,そのソフトの分類を検出し,データベースを構築しています.ソフトウェアの重要な情報元となるウェブページは,WebArchiveへのリンクを用いることで保存しています.なるべく人手をかけずにデータを収集し,メンテナンスコストがかかるような機能の実装は避けて,今日まで運用されています.

おわりに

数学ソフトウェアの情報サービス,swMATHを紹介しました.swMathは2011年あたりからプロジェクトがスタートしたようなのですが,利用者数,登録されているソフトウェア,論文数は増え続けているようです.ちなみに,データベースに登録されていないソフトウェアは登録を申請することもできます.
f:id:soma_arc:20171211193330p:plain
グラフはswMATH-2chartsより

数学に関連した情報サービスとして,とても面白いサービスだと思います.一度チェックしてみてはいかがでしょうか.

参考文献

  • Sebastian B¨onisch, Michael Brickenstein, Hagen Chrapary, Gert-Martin Greuel, and Wolfram Sperber.swMATH – A New Information Service for Mathematical Software
  • Gert-Martin Greuel and Wolfram Sperber. swMATH - An Information Service for Mathematical Software. In ICMS2014 Conference Proceedings
  • Hagen Chrapary and Yue Ren. The Software Portal swMATH: A State of the Art Report and Next Steps. In ICMS2016 Conference Proceedings
  • Helge Holzmann, Mila Runnwerth, Wolfram Sperber. Linking Mathematical Software in Web Archive. In ICMS2016 Conference Proceedings (arXiv)

Tokyo Demo Fest 2017に参加しました

2月18日と2月19日に3331 Arts Chiyodaで開催されたTokyo Demo Fest(TDF)2017に参加しました.
tokyodemofest.jp
去年に引き続き,二回目の参加になります.去年の参加記事はこちら.
TokyoDemoFest2016に参加しました - 心鏡曼荼羅
TDFは日本で唯一のデモパーティです.詳しくはTDF公式サイトを参照してもらいたいのですが,すごく簡単にいうとデモと呼ばれるリアルタイムに映像や音楽を生成するプログラムを発表したり,参加者同士で交流したりする場です.作品発表の様子はYouTubeでもみることができます.
Tokyo Demo Fest 2017 - Combined Demo Compo [Live footage] - YouTube

今回はCombined Graphics CompoとGLSL Graphics Compoに作品を投稿し,Combined Graphicsで二位,GLSL Graphicsで一位をいただきました.



あまり制作に時間は取れなかったのですが,評価していただけて嬉しいです.簡単に作品の紹介と解説をつけておきます.

TokyoDemoFesTessellation

f:id:soma_arc:20170224173532p:plain
https://www.shadertoy.com/view/MsscR4
Tokyo Demo Festの頭文字TDFを使って平面を敷き詰めた作品です.テセレーションというとポリゴン分割のことを考える人が多いかと思いますが,平面や空間を図形で敷き詰めることもテセレーションといいます.
以前から普通の双曲タイリングは描いていましたが,こういったデザインのものを作ったのは初めてです.本当は回転を含むような複雑なパターンにしたかったのですが,三種類のアルファベット展開というのはなかなかに難しいです.あちらをたてればこちらがたたず…割と単純なものになってしまいました.アイデアメモの一部を置いておきます.

グリッドで考えて塗りつぶすような実装にしました.

また,GLSLで描いたのでこちらの記事([メガデモ]VisualStudio2013で4kbのexeを作る | notargs.com)を参考にexecutableにしてPC 4k Graphicsとして投稿してみました.本番では解像度の違いから図が左に寄ってしまうというミスがありましたが…
f:id:soma_arc:20170225184307p:plain
元々小さなファイルサイズというのはあまり興味がなかったのですが,案外面白かったので,次回もPC 4k Graphicsに出してみたいところです.

4kは手間がかかりそうだと思っていたのですが,こちらの発表資料(デモシーンへようこそ-4KBで映像作品をつくる技術-.pptx - Google ドライブ)や,i_saintさんのテンプレート(GitHub - i-saint/CEDEC2016_4kintro)で捗りそうです.

こういった敷き詰め作品を作りたい方には,こちらの本がおすすめです.
www.amazon.co.jp

また,GLSLで描く方法については少し前に記事を書きました.
GLSLで描くTessellation - 心鏡曼荼羅
GLSLで描くHyperbolic Tessellation - 心鏡曼荼羅

Schottky Waltz

f:id:soma_arc:20170224173544p:plain
https://www.shadertoy.com/view/XslyzH


これはCircle inversion fractalやSelf inversion fractalと呼ばれる種類のフラクタルです.ショットキー変換と呼ばれる,円のペアから定義される変換の名前に由来しています.

このフラクタルは円の外側と内側を入れ替える"反転"という操作を繰り返すことで描かれます.例えば,下図のような4つの円の反転で構成されるものを考えてみましょう.

1つの円の反転によって,他の3つの円は反転を行った円の内側に移されます.そのため,4つの円で反転を行うと,それらの内側に12個の小円ができます.

新たにできた12個の小円に対しても4つの円の反転を作用させると小円の下にさらに円ができます.

これを更に繰り返していくことで,円が無限に連なる図を得ることができます.

これをGLSLで描画する詳しいアルゴリズムに興味がある方はこちらの記事(GLSLで描くCircle Inversion Fractals - 心鏡曼荼羅)をご覧ください.

アルゴリズムは単純なので,面白い図を生成する円の配置を見つけて,動きをつけます.円の配置を試す際には拙作のSchottkyLinkを用いました.円や球で構成されるフラクタルを直観的に構成することができるウェブアプリケーションです.普段からこういったフラクタルの可視化等に取り組んでいるのでこういった形で役に立つのは嬉しいですね.まだバグも多いのですが,興味があればお試しください.
f:id:soma_arc:20170225211055p:plain

動きをつける際に参考にさせていただいたのはこちらの記事(2Dの小技 動くお絵かき - Qiita)と各種のeasingが紹介されており,実際に動きをみることができるこちらのサイト(http://gizma.com/easing/)です.前々からこのフラクタルに小気味良いモーションを付けたいと思っていたので実現できてよかったです.

去年投稿した作品はこのフラクタルの三次元版でした.形状は気に入っていたのですが,シェーディングがわりとお粗末だったのがマイナスポイントでした.三次元は魅せ方が難しいうえ,時間もなかったので,二次元にしてわかりやすい形にしたのは良い選択だったと思います.

おわりに

今回のTDFではやりたいと思っていたことを試す非常に良い機会になりました.他の作品にも創作意欲を刺激されたので来年に備えたいですね.また,各部門にもっと作品が増えて欲しいと思っておりますので,興味のある方は是非参加,作品発表をしていただければと思います.
それでは,Demoscenerの皆様の一年がより良きものでありますように.

OpenVDBのビルドメモ

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

準備

以下のリンクからソース(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, pythonpythonモジュールが必要であれば)を用意し,ヘッダファイルとライブラリへのパスを設定します.あとはCONCURRENT_MALLOC_LIBをtbbを使うように設定すればビルドできるはずです.

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

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

以上となります.

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
以上となります.

GLSLで描くCircle Inversion Fractals

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

前回の記事で紹介した円の反転という操作を用いてフラクタルを描いてみます.
f:id:soma_arc:20161220183920p:plain
まずはこの画像のような4つの円板とその反転を考えます.それぞれの円の反転は自分以外の円板を自分の内側に移します.よって,それぞれの円の内側には3つの円が移され,12個の円ができています.次に,ショットキー円の反転を新たに現れた小円にも適用することを繰り替えすと,以下のように円が入れ子状に続いていきます.
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)から来ています.ショットキー変換は2つの円をペアとする変換で,相手の円板の外側を自分の内側に,自分の円板の外側を相手の円板の内側へと移す変換です.これは,本質的に円の反転と同じものになることがあります.ショットキー変換を構成する円のことをショットキー円と呼ぶのですが,奇数個の円を扱うことを考えると,厳密にはショットキー円とは呼べません.詳しくは,後述のインドラの真珠を読んでみてください.

Further Reading

今回紹介したようなフラクタルや数学的背景に興味のある方は必読です.

このアルゴリズムを発表した論文とスライドです.三次元への拡張に関する事項も書いてあります.

おわりに

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