*異方性カーネル [#n8bdf672]
パーティクルから陰関数を生成する場合，以下のような平滑化カーネルを用いる．
#ref(anisotropic_kernel.eq1.gif,nolink,70%)
ここで，&ref(anisotropic_kernel.eq2.gif,nolink,70%);はスケーリング定数，&ref(anisotropic_kernel.eq3.gif,nolink,70%);はシミュレーション次元，
&ref(anisotropic_kernel.eq4.gif,nolink,70%);は距離が遠くなるほど値が滑らかに小さくなるローカルサポートの対称関数である．
この定義はいわゆるメタボールと同じようなものであり，
生成された表面は凸凹になってしまう(下図左)．このような表面はblobbyであると呼ばれる．
blobbyな表面を改善するために，Yuら
¬e{Yu2010: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};
は従来の等方性のカーネルの代わりに異方性カーネル(anisotropic kernel)を導入した．
つまり，&ref(anisotropic_kernel.eq5.gif,nolink,70%);の代わりに正定値行列&ref(anisotropic_kernel.eq6.gif,nolink,70%);を用いて，
#ref(anisotropic_kernel.eq7.gif,nolink,70%)
とすることで，ベクトル&ref(anisotropic_kernel.eq8.gif,nolink,70%);に回転と伸縮を与え，球の代わりに楕円体形状のカーネルとする．
法線方向に短軸がくるような楕円体を設定することで，滑らかでかつエッジがよりくっきりと出るような表面を得ることができる(下図右)．
また，表面付近でのパーティクルの不規則な分布を解消するために，
パーティクル位置に関しても平滑化を行い更新位置&ref(anisotropic_kernel.eq9.gif,nolink,70%);を計算する．
全パーティクルに対して，&ref(anisotropic_kernel.eq10.gif,nolink,70%);と&ref(anisotropic_kernel.eq11.gif,nolink,70%);を計算することで以下の陰関数場を得る．
#ref(anisotropic_kernel.eq12.gif,nolink,70%)
|#ref(anisotropic_kernel.jpg,nolink,100%)|
|CENTER:左は等方性カーネル，右は異方性カーネルを使った表面|
以下で¬e{Yu2010};の&ref(anisotropic_kernel.eq10.gif,nolink,70%);，&ref(anisotropic_kernel.eq11.gif,nolink,70%);の計算方法について述べる．
**パーティクル位置更新 [#da028bd8]
不規則なパーティクル配置を修正するために，以下のLaplacianスムージングでパーティクル中心位置を更新する．
#ref(anisotropic_kernel.eq13.gif,nolink,70%)
&ref(anisotropic_kernel.eq14.gif,nolink,70%);は定数で，¬e{Yu2010};は0.9から1を推奨している．
&ref(anisotropic_kernel.eq15.gif,nolink,70%);は重み関数で，
#ref(anisotropic_kernel.eq16.gif,nolink,70%)
とする．この重み関数は異方性を算出するための共分散行列の計算に用いられる．
異方性をより正確に推定するためにはより多くのパーティクルが必要であるため，
有効半径&ref(anisotropic_kernel.eq5.gif,nolink,70%);はシミュレーションで用いる有効半径の2倍にする．
注意として，&ref(anisotropic_kernel.eq9.gif,nolink,70%);は表面生成のみに用い，シミュレーションには使わない．
そのため，更新前の位置&ref(anisotropic_kernel.eq17.gif,nolink,70%);を上書きするのではなく，別の変数に格納するようにする．
***共分散行列の計算 [#u6e862f8]
&ref(anisotropic_kernel.eq6.gif,nolink,70%);を計算するために，重み付き主成分分析(WPCA:Weighted Principal Component Analysis)¬e{Koren2003:Y. Koren and L. Carmel, Visualization of labeled data using linear transformations, In Proceedings of IEEE Information Visualization, 2003.};を用いる．
主成分分析(PCA)は画像認識などでよく用いられるデータ解析手法で，
分散したデータの特徴をよくとらえた新しい軸(主軸)を設定することでデータ量を減らしたりすることに用いられている．
WPCAはPCAに重みをつけることで，外れ値(outlier)やノイズに対してロバストにしたものである．
WPCAではまず，データ点での重みを計算し，次に重みを考慮した共分散行列(covariance matrix)&ref(anisotropic_kernel.eq18.gif,nolink,70%);を算出，
最後に&ref(anisotropic_kernel.eq18.gif,nolink,70%);の固有値解析により固有ベクトルをもとめ，大きな固有値を持つ固有ベクトルを新しい主軸とする．
ここでは，WPCAの結果から&ref(anisotropic_kernel.eq6.gif,nolink,70%);を求める．
まず，重み付き平均位置&ref(anisotropic_kernel.eq19.gif,nolink,70%);を計算する．
#ref(anisotropic_kernel.eq20.gif,nolink,70%)
&ref(anisotropic_kernel.eq19.gif,nolink,70%);はパーティクル&ref(anisotropic_kernel.eq21.gif,nolink,70%);の近傍パーティクルの中心を示している．
重み&ref(anisotropic_kernel.eq15.gif,nolink,70%);を使うことで，&ref(anisotropic_kernel.eq22.gif,nolink,70%);から大きく離れた外れ値の影響は小さくなる．
&ref(anisotropic_kernel.eq19.gif,nolink,70%);から各近傍パーティクルの位置&ref(anisotropic_kernel.eq23.gif,nolink,70%);へのベクトルの積をとることで，共分散行列を得る．
#ref(anisotropic_kernel.eq24.gif,nolink,70%)
重み関数&ref(anisotropic_kernel.eq15.gif,nolink,70%);にはパーティクル位置更新で使ったものと同じものを用いる．
***共分散行列の特異値分解 [#e0594448]
楕円の伸縮方向を得るために共分散行列&ref(anisotropic_kernel.eq18.gif,nolink,70%);を特異値分解して，固有値，固有ベクトルを求める．
#ref(anisotropic_kernel.eq25.gif,nolink,70%)
&ref(anisotropic_kernel.eq26.gif,nolink,70%);は固有ベクトルを各列に持つ回転行列，&ref(anisotropic_kernel.eq27.gif,nolink,70%);は対角要素に固有値を持つ対角行列である．
大きな変形や近傍にパーティクルが少ない孤立したパーティクルなどに対応するために以下の処理を行い&ref(anisotropic_kernel.eq27.gif,nolink,70%);を修正する．
-固有値をチェックして大きな変形を防ぐ．
固有値は大きい順に並んでいるとして，&ref(anisotropic_kernel.eq28.gif,nolink,70%);となる&ref(anisotropic_kernel.eq29.gif,nolink,70%);があったら，&ref(anisotropic_kernel.eq30.gif,nolink,70%);で置き換える(&ref(anisotropic_kernel.eq31.gif,nolink,70%);)．
-近傍パーティクル数をチェックして，孤立したパーティクルを検出する．
もし近傍パーティクル数が&ref(anisotropic_kernel.eq32.gif,nolink,70%);より少ない場合，&ref(anisotropic_kernel.eq33.gif,nolink,70%);とする．ここで，&ref(anisotropic_kernel.eq34.gif,nolink,70%);は単位行列，&ref(anisotropic_kernel.eq35.gif,nolink,70%);はカーネルの大きさを制御するパラメータ．
-近傍にパーティクルがフルにある場合の体積を一定にするために，&ref(anisotropic_kernel.eq36.gif,nolink,70%);となるような&ref(anisotropic_kernel.eq37.gif,nolink,70%);を導入する．
これらを適用して，修正した共分散行列&ref(anisotropic_kernel.eq38.gif,nolink,70%);は以下となる．
#ref(anisotropic_kernel.eq39.gif,nolink,70%)
#ref(anisotropic_kernel.eq40.gif,nolink,70%)
ここで，&ref(anisotropic_kernel.eq41.gif,nolink,70%);である．
¬e{Yu2010};は，&ref(anisotropic_kernel.eq42.gif,nolink,70%);を用いている．
***異方性カーネル用行列$G$の算出 [#md43c27c]
&ref(anisotropic_kernel.eq6.gif,nolink,70%);は&ref(anisotropic_kernel.eq38.gif,nolink,70%);の逆行列に&ref(anisotropic_kernel.eq43.gif,nolink,70%);をかけたものになる．
#ref(anisotropic_kernel.eq44.gif,nolink,70%)
*結果 [#ta3ab463]
[[SPH法の実装の結果:SPH法の実装(GPU実装含む)#k4a9f2ef]]で示した立方体形状の水塊がなだらかなスロープ上に落下するシーンでの結果をしめす．パーティクル数は約25,000，メッシュ化(MC法)のグリッド解像度は256x103x103である．
まず，共分散行列を球に適用することで，パーティクルを楕円体として描画した結果を示す．
|&ref(dam_break1_ap.jpg,nolink,100%);|&ref(dam_break1_np.jpg,nolink,100%);|
|>|CENTER:左は異方性カーネル，右は従来の等方性カーネルによる結果|
メッシュを抽出して，GLSLで屈折表面として描画した結果を以下に示す．
|&ref(dam_break1_ar.jpg,nolink,100%);|&ref(dam_break1_nr.jpg,nolink,100%);|
|>|CENTER:左は異方性カーネル，右は従来の等方性カーネルによる結果|
異方性カーネルで表面が滑らかになるとともに，両サイドの壁に沿ってできた薄い水膜もうまく再現できている．
さらにシミュレーションを進めた時の結果を以下に示す．
|&ref(dam_break2_ap.jpg,nolink,100%);|&ref(dam_break2_np.jpg,nolink,100%);|
|>|CENTER:左は異方性カーネル，右は従来の等方性カーネルによる結果|
拡大したものは以下
|&ref(dam_break3_ap.jpg,nolink,100%);|&ref(dam_break3_np.jpg,nolink,100%);|
|>|CENTER:左は異方性カーネル，右は従来の等方性カーネルによる結果|
メッシュを抽出して，GLSLで屈折表面として描画した結果を以下に示す．
|&ref(dam_break2_ar.jpg,nolink,100%);|&ref(dam_break2_nr.jpg,nolink,100%);|
|>|CENTER:左は異方性カーネル，右は従来の等方性カーネルによる結果|
計算時間はだいたい，SPH計算に5ms/frame，異方性カーネルのGの計算に10ms/frame，異方性カーネルを適用した上でのメッシュ化に145ms/frame(従来の等方性カーネルの場合は15ms/frame)であった．ちなみにSPHでの計算，異方性の計算，および，MC法によるメッシュ化はすべてCUDAを使ってGPU上で実装し，GeForceGTX580上でテストした．
パーティクルを直接描画する際に，等方性カーネルはポイントスプライトを使ったが，異方性カーネルの場合はglutSolidSphereに共分散行列を変換行列として設定して楕円体として描画した．そのため，描画に時間がかかっているので注意(私の環境だと異方性カーネルで85ms/frame, ポイントスプライトによる等方性カーネルは15ms/frame)．
*サンプルコード(GPU実装) [#b03ec21e]
[[SPH法の実装(GPU実装含む)]]で示したSPH法のコードに上記の異方性カーネルの実装を追加し，さらにFLTKでGUIを追加したコードを以下に置く．SPH計算，異方性計算，および，MC法によるメッシュ化はすべてCUDAを使って実装した．
#ref(rx_anisotropic.zip);
#ref(rx_anisotropic_v1.1.zip);
-ビルドするのに必要なライブラリ
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:http://www.slis.tsukuba.ac.jp/~fujis/cgi-bin/wiki/index.php?GSL]]を使っている．GPU版における特異値分解はNumrical Recipesのコードをデバイス関数として移植した．
--文字描画に[[FTGL:http://homepages.paradise.net.nz/henryj/code/]]を用いている．デフォルトではInconsolataフォントを使っている．フォント指定はrx_fltk_glcanvas.cppの103行目あたり．
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
