異方性カーネル

パーティクルから陰関数を生成する場合,以下のような平滑化カーネルを用いる.

anisotropic_kernel.eq1.gif

ここで,anisotropic_kernel.eq2.gifはスケーリング定数,anisotropic_kernel.eq3.gifはシミュレーション次元, anisotropic_kernel.eq4.gifは距離が遠くなるほど値が滑らかに小さくなるローカルサポートの対称関数である. この定義はいわゆるメタボールと同じようなものであり, 生成された表面は凸凹になってしまう(下図左).このような表面はblobbyであると呼ばれる.

blobbyな表面を改善するために,Yuら (1) は従来の等方性のカーネルの代わりに異方性カーネル(anisotropic kernel)を導入した. つまり,anisotropic_kernel.eq5.gifの代わりに正定値行列anisotropic_kernel.eq6.gifを用いて,

anisotropic_kernel.eq7.gif

とすることで,ベクトルanisotropic_kernel.eq8.gifに回転と伸縮を与え,球の代わりに楕円体形状のカーネルとする. 法線方向に短軸がくるような楕円体を設定することで,滑らかでかつエッジがよりくっきりと出るような表面を得ることができる(下図右). また,表面付近でのパーティクルの不規則な分布を解消するために, パーティクル位置に関しても平滑化を行い更新位置anisotropic_kernel.eq9.gifを計算する. 全パーティクルに対して,anisotropic_kernel.eq10.gifanisotropic_kernel.eq11.gifを計算することで以下の陰関数場を得る.

anisotropic_kernel.eq12.gif
anisotropic_kernel.jpg
左は等方性カーネル,右は異方性カーネルを使った表面

以下で(1)anisotropic_kernel.eq10.gifanisotropic_kernel.eq11.gifの計算方法について述べる.

パーティクル位置更新

不規則なパーティクル配置を修正するために,以下のLaplacianスムージングでパーティクル中心位置を更新する.

anisotropic_kernel.eq13.gif

anisotropic_kernel.eq14.gifは定数で,(1)は0.9から1を推奨している. anisotropic_kernel.eq15.gifは重み関数で,

anisotropic_kernel.eq16.gif

とする.この重み関数は異方性を算出するための共分散行列の計算に用いられる. 異方性をより正確に推定するためにはより多くのパーティクルが必要であるため, 有効半径anisotropic_kernel.eq5.gifはシミュレーションで用いる有効半径の2倍にする.

注意として,anisotropic_kernel.eq9.gifは表面生成のみに用い,シミュレーションには使わない. そのため,更新前の位置anisotropic_kernel.eq17.gifを上書きするのではなく,別の変数に格納するようにする.

共分散行列の計算

anisotropic_kernel.eq6.gifを計算するために,重み付き主成分分析(WPCA:Weighted Principal Component Analysis)(2)を用いる. 主成分分析(PCA)は画像認識などでよく用いられるデータ解析手法で, 分散したデータの特徴をよくとらえた新しい軸(主軸)を設定することでデータ量を減らしたりすることに用いられている. WPCAはPCAに重みをつけることで,外れ値(outlier)やノイズに対してロバストにしたものである.

WPCAではまず,データ点での重みを計算し,次に重みを考慮した共分散行列(covariance matrix)anisotropic_kernel.eq18.gifを算出, 最後にanisotropic_kernel.eq18.gifの固有値解析により固有ベクトルをもとめ,大きな固有値を持つ固有ベクトルを新しい主軸とする. ここでは,WPCAの結果からanisotropic_kernel.eq6.gifを求める.

まず,重み付き平均位置anisotropic_kernel.eq19.gifを計算する.

anisotropic_kernel.eq20.gif

anisotropic_kernel.eq19.gifはパーティクルanisotropic_kernel.eq21.gifの近傍パーティクルの中心を示している. 重みanisotropic_kernel.eq15.gifを使うことで,anisotropic_kernel.eq22.gifから大きく離れた外れ値の影響は小さくなる.

anisotropic_kernel.eq19.gifから各近傍パーティクルの位置anisotropic_kernel.eq23.gifへのベクトルの積をとることで,共分散行列を得る.

anisotropic_kernel.eq24.gif

重み関数anisotropic_kernel.eq15.gifにはパーティクル位置更新で使ったものと同じものを用いる.

共分散行列の特異値分解

楕円の伸縮方向を得るために共分散行列anisotropic_kernel.eq18.gifを特異値分解して,固有値,固有ベクトルを求める.

anisotropic_kernel.eq25.gif

anisotropic_kernel.eq26.gifは固有ベクトルを各列に持つ回転行列,anisotropic_kernel.eq27.gifは対角要素に固有値を持つ対角行列である. 大きな変形や近傍にパーティクルが少ない孤立したパーティクルなどに対応するために以下の処理を行いanisotropic_kernel.eq27.gifを修正する.

  • 固有値をチェックして大きな変形を防ぐ. 固有値は大きい順に並んでいるとして,anisotropic_kernel.eq28.gifとなるanisotropic_kernel.eq29.gifがあったら,anisotropic_kernel.eq30.gifで置き換える(anisotropic_kernel.eq31.gif).
  • 近傍パーティクル数をチェックして,孤立したパーティクルを検出する. もし近傍パーティクル数がanisotropic_kernel.eq32.gifより少ない場合,anisotropic_kernel.eq33.gifとする.ここで,anisotropic_kernel.eq34.gifは単位行列,anisotropic_kernel.eq35.gifはカーネルの大きさを制御するパラメータ.
  • 近傍にパーティクルがフルにある場合の体積を一定にするために,anisotropic_kernel.eq36.gifとなるようなanisotropic_kernel.eq37.gifを導入する. これらを適用して,修正した共分散行列anisotropic_kernel.eq38.gifは以下となる.
    anisotropic_kernel.eq39.gif
anisotropic_kernel.eq40.gif

ここで,anisotropic_kernel.eq41.gifである. (1)は,anisotropic_kernel.eq42.gifを用いている.

異方性カーネル用行列$G$の算出

anisotropic_kernel.eq6.gifanisotropic_kernel.eq38.gifの逆行列にanisotropic_kernel.eq43.gifをかけたものになる.

anisotropic_kernel.eq44.gif

結果

SPH法の実装の結果:SPH法の実装(GPU実装含む)で示した立方体形状の水塊がなだらかなスロープ上に落下するシーンでの結果をしめす.パーティクル数は約25,000,メッシュ化(MC法)のグリッド解像度は256x103x103である.

まず,共分散行列を球に適用することで,パーティクルを楕円体として描画した結果を示す.

dam_break1_ap.jpgdam_break1_np.jpg
左は異方性カーネル,右は従来の等方性カーネルによる結果

メッシュを抽出して,GLSLで屈折表面として描画した結果を以下に示す.

dam_break1_ar.jpgdam_break1_nr.jpg
左は異方性カーネル,右は従来の等方性カーネルによる結果

異方性カーネルで表面が滑らかになるとともに,両サイドの壁に沿ってできた薄い水膜もうまく再現できている. さらにシミュレーションを進めた時の結果を以下に示す.

dam_break2_ap.jpgdam_break2_np.jpg
左は異方性カーネル,右は従来の等方性カーネルによる結果

拡大したものは以下

dam_break3_ap.jpgdam_break3_np.jpg
左は異方性カーネル,右は従来の等方性カーネルによる結果

メッシュを抽出して,GLSLで屈折表面として描画した結果を以下に示す.

dam_break2_ar.jpgdam_break2_nr.jpg
左は異方性カーネル,右は従来の等方性カーネルによる結果

計算時間はだいたい,SPH計算に5ms/frame,異方性カーネルのGの計算に10ms/frame,異方性カーネルを適用した上でのメッシュ化に145ms/frame(従来の等方性カーネルの場合は15ms/frame)であった.ちなみにSPHでの計算,異方性の計算,および,MC法によるメッシュ化はすべてCUDAを使ってGPU上で実装し,GeForceGTX580上でテストした.

パーティクルを直接描画する際に,等方性カーネルはポイントスプライトを使ったが,異方性カーネルの場合はglutSolidSphereに共分散行列を変換行列として設定して楕円体として描画した.そのため,描画に時間がかかっているので注意(私の環境だと異方性カーネルで85ms/frame, ポイントスプライトによる等方性カーネルは15ms/frame).

サンプルコード(GPU実装)

SPH法の実装(GPU実装含む)で示したSPH法のコードに上記の異方性カーネルの実装を追加し,さらにFLTKでGUIを追加したコードを以下に置く.SPH計算,異方性計算,および,MC法によるメッシュ化はすべてCUDAを使って実装した.

  • ビルドするのに必要なライブラリ
    FLTK, freeglut, GLEW, CUDA, boost, FTGL, freetype, GSL
    各ライブラリについてはライブラリのインストールを参照.

さらに,3Dモデルファイルの入出力のrx_model.libも必要(ver0.5以上).

  • 簡単な説明&注意事項
    • Visual Studioから実行する場合は,作業ディレクトリを"../../bin"にしておく.
    • binフォルダに"sph_scene_*.cfg"という名前のファイルがあり,これがシーンを記述するファイル. プログラム起動時にファイルを読み取り,Sceneメニューに列挙する. デフォルトでは6つのシーンを用意している(合計12まで追加可能).
    • rxSPH, rxSPH_GPUクラスのCalAnisotropicKernel関数が行列Gを計算する関数.rxSPHはSPH,行列GともにCPUで,rxSPH_GPUはともにGPU(CUDAで実装)で計算する.
    • CalAnisotropicKernel関数のCPU版は特異値分解にGSLを使っている.GPU版における特異値分解はNumrical Recipesのコードをデバイス関数として移植した.
    • 文字描画にFTGLを用いている.デフォルトではInconsolataフォントを使っている.フォント指定はrx_fltk_glcanvas.cppの103行目あたり.

(1)  J. Yu and G. Turk, Reconstructing Surfaces of Particle-Based Fluids Using Anisotropic Kernels, In Proceedings of the 2010 ACM SIGGRAPH/Eurographics symposium on Computer animation, 2010
(2)  Y. Koren and L. Carmel, Visualization of labeled data using linear transformations, In Proceedings of IEEE Information Visualization, 2003.

添付ファイル: filerx_anisotropic_v1.1.zip 1210件 [詳細] filedam_break3_np.jpg 632件 [詳細] filedam_break2_ap.jpg 615件 [詳細] filedam_break1_np.jpg 570件 [詳細] filecube_np.jpg 604件 [詳細] filecube_ap.jpg 600件 [詳細] filedam_break3_ap.jpg 636件 [詳細] filedam_break2_np.jpg 552件 [詳細] filedam_break2_nr.jpg 535件 [詳細] filedam_break1_nr.jpg 552件 [詳細] filedam_break1_ar.jpg 565件 [詳細] filedam_break1_ap.jpg 556件 [詳細] filedam_break2_ar.jpg 549件 [詳細] fileanisotropic_kernel.jpg 524件 [詳細] fileanisotropic_kernel.eq3.gif 556件 [詳細] fileanisotropic_kernel.eq42.gif 536件 [詳細] fileanisotropic_kernel.eq19.gif 587件 [詳細] fileanisotropic_kernel.eq7.gif 566件 [詳細] fileanisotropic_kernel.eq40.gif 570件 [詳細] fileanisotropic_kernel.eq43.gif 551件 [詳細] fileanisotropic_kernel.eq27.gif 585件 [詳細] fileanisotropic_kernel.eq20.gif 556件 [詳細] fileanisotropic_kernel.eq9.gif 543件 [詳細] fileanisotropic_kernel.eq1.gif 584件 [詳細] fileanisotropic_kernel.eq26.gif 543件 [詳細] fileanisotropic_kernel.eq30.gif 557件 [詳細] fileanisotropic_kernel.eq29.gif 521件 [詳細] fileanisotropic_kernel.eq5.gif 524件 [詳細] fileanisotropic_kernel.eq38.gif 534件 [詳細] fileanisotropic_kernel.eq4.gif 536件 [詳細] fileanisotropic_kernel.eq18.gif 627件 [詳細] fileanisotropic_kernel.eq33.gif 525件 [詳細] fileanisotropic_kernel.eq10.gif 548件 [詳細] fileanisotropic_kernel.eq6.gif 542件 [詳細] fileanisotropic_kernel.eq13.gif 560件 [詳細] fileanisotropic_kernel.eq44.gif 546件 [詳細] fileanisotropic_kernel.eq17.gif 603件 [詳細] fileanisotropic_kernel.eq24.gif 583件 [詳細] fileanisotropic_kernel.eq31.gif 544件 [詳細] fileanisotropic_kernel.eq32.gif 516件 [詳細] fileanisotropic_kernel.eq14.gif 599件 [詳細] fileanisotropic_kernel.eq23.gif 562件 [詳細] fileanisotropic_kernel.eq36.gif 550件 [詳細] fileanisotropic_kernel.eq34.gif 515件 [詳細] fileanisotropic_kernel.eq15.gif 566件 [詳細] fileanisotropic_kernel.eq28.gif 568件 [詳細] fileanisotropic_kernel.eq21.gif 577件 [詳細] fileanisotropic_kernel.eq25.gif 591件 [詳細] fileanisotropic_kernel.eq35.gif 567件 [詳細] fileanisotropic_kernel.eq22.gif 512件 [詳細] fileanisotropic_kernel.eq39.gif 529件 [詳細] fileanisotropic_kernel.eq41.gif 503件 [詳細] fileanisotropic_kernel.eq11.gif 560件 [詳細] fileanisotropic_kernel.eq37.gif 504件 [詳細] fileanisotropic_kernel.eq8.gif 549件 [詳細] fileanisotropic_kernel.eq16.gif 637件 [詳細] fileanisotropic_kernel.eq2.gif 555件 [詳細] fileanisotropic_kernel.eq12.gif 580件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2011-11-28 (月) 13:12:15 (3226d)