要素間補間による共回転系弾性体シミュレーションの高速化
6
0
0
全文
(2) Vol.2015-CG-160 No.8 2015/8/29. 情報処理学会研究報告 IPSJ SIG Technical Report. すものとする.. 4.1 陰解法 本章では基本的な陰解法の枠組みを復習するとともに, 最適化問題へと置き換える方法を説明する.. n 個の頂点からなる 3 次元の力学系を考える. 時刻 t に おいて頂点の位置と速度を与えるベクトルを, それぞれ. x(t) ∈ R3n , v(t) ∈ R3n とし, ポテンシャルを与える関数を E : R3n → R とすると, 陰解法の位置と速度の更新の式は 以下のように表すことが出来る [1].. x(t+1) = x(t) + hv(t+1). (1). v(t+1) = v(t) − hM−1 ∇E |x(t+1). (2). h は時間の刻み幅であり, M ∈ R3n×3n は質量マトリック スを表している. 代入によって, 式 (15), (2) から速度を消 去して整理すると, 離散化されたニュートンの運動方程式. ( ) M x(t+1) − 2x(t) + x(t−1) = −h2 ∇E |x(t+1). (3). を得る. 式 (3) を解くために, ∇E |x(t+1) の部分を以下のよ うに近似する [1].. ( )( ) ∇E |x(t+1) ≈ ∇E |x(t) + ∇2 E |x(t) x(t+1) − x(t) (4). 図 1 研究概要. Fig. 1 Overview of our method. ここで, 式 (4) を式 (3) に代入して, 整理すると以下の式 (5) を得る.. は McAdams ら [11] や Georgii ら [4] の Multigulid を用い. M (x − y) = −h2 E (x). る手法である. Multigulid を用いるアプローチには解の安 定性が高いという利点もある.. (5). ただし, x = x(t+1) , y = 2x(t) − x(t−1) とした. 式 (5) を解. 3 つ目はシステム行列の扱いを工夫することで高速化を. くことで, 既知の頂点の位置の情報 y から次の時刻におけ. 行うアプローチである. Hecht ら [5] はモデルを領域分割す. る頂点の位置 x を得ることが出来る. ここで, 式 (5) を解. ることで, システム行列を複数のブロックに分割し, 誤差が. いた結果得られる x は, 以下の式 (6) によって与えられる. 大きくなったブロックの成分のみを更新することによりシ. 関数 g(x) の最小値を与える x に等しいため, 式 (5) を解く. ミュレーションを高速化する手法を提案した.. 問題を関数 g(x) の最小化問題として置き換えることが出. 本研究は 1 つ目のアプローチに近い手法である. 本手法. 来る.. ではサンプリングされた代表要素の回転行列の値から残り 全ての要素の回転行列を近似的に補間することで, シミュ レーションの高速化を行う. これは, 捉え方によっては要 素の回転の自由度を減らすことによる高速化と捉えること が出来る.. 3. 背景・表記法 4. 関連研究. g(x) =. 1 T (x − y) M (x − y) + h2 E (x) 2. (6). 実際に, 関数 g(x) の極値を求める式, ∇g(x) = 0 は式 (5) そ のものである. このように最適化問題として置き換えた陰解 法は Martin ら [10] によって提案され, variational implicit. Euler と呼ばれている. 以上の枠組みでは摩擦による減衰などは一切考えていな い. ダンピングの影響を考慮するためには, Liu ら [9] の手. 本稿を通して, 以下のような表記を用いる. Ψ のような. 法に従って, y = x(t) + h˜ v(t) と置き換えればよい. ここで,. 細字の文字はスカラーを, x のようなボールド体の小文字. ˜ (t) はダンピングを考慮した速度であり, 最も単純なモデ v. ˆ はベクトルを, F のようなボールド体の大文字は行列を, q. ˜ (t) = αv(t) ルでは, ダンピング係数 α ∈ [0, 1] を用いて, v. のようなハット付きの文字はクォータニオンをそれぞれ表. と表せる.. ⓒ 2015 Information Processing Society of Japan. 2.
(3) Vol.2015-CG-160 No.8 2015/8/29. 情報処理学会研究報告 IPSJ SIG Technical Report. 4.2 共回転系弾性体. テンソル積である.. 本章では提案手法で用いた弾性体モデルである, 共回転 系弾性体の力学的な取り扱いについて説明する.. よって, 重力などの外力を f ext と書くことにすると, ポ テンシャルエネルギーは以下の式 (12) のように表せる.. 3 次元空間中の弾性体を変形させたときに蓄えられる歪. E (x, r) = xT Lx − 2xT Jr + xT f ext. みエネルギー E は, 弾性体が占める領域を Ω とすると, 式. (7) のように表すことが出来る. ∫ E= Ψ (x, F (x)) dx. (12). 式 (12) を陰解法で最小化すべき関数の式 (6) に代入して整 理すると, 最終的に最小化すべき式 (13) が得られる.. (7). Ω. ここで, F は位置 x における弾性体の変形勾配テンソルで あり, 各成分は Fij = ∂ϕi /∂xj のように表される. また,. ) ( ) 1 T( x M + 2h2 L x − xT My + h2 (2Jr − f ext ) (13) 2 4.4 Local/Global Optimization. Ψ (x, F (x)) 歪みエネルギー密度であり, その具体的な関数. 本章では, 前章で示した式 (13) の最小化問題の解法の 1. 系は採用する弾性体のモデルを定めることで決定される.. つとして, Sokine ら [13] によって提案された Local/Global. 弾性体モデルとして共回転系弾性体を採用した場合, 歪み. approach を用いた最適化の方法を説明する. Local/Global. エネルギー密度は式 (8) のように表すことが出来る.. Ψcor. ( ) λ 2 = µ||F − R||F + tr2 RT F − I 2. approach では Local step と Global step の 2 つのステップ (8). を繰り返すことにより目的関数の最小化を行う.. Local step では頂点ベクトル x を固定し, ベクトル r に. µ, λ はラメ定数であり, R は極分解 F = RS によって得ら. ついて目的関数 g を最小化する. 具体的には, x を固定し. れる回転行列である. 式 (8) において, 解の安定性に大き. た状態で, すべての要素に対して特異値分解 Fi = Ui Σi Qi. く影響するのは µ||F −. 2 R||F. の部分であるので McAdams. ら [11] の手法に従って, 第 2 項を落として近似すると,. Ψcor = µ||F −. 2 R||F. T を行い, Ri = QT i Ui によって g を最小化する回転行列. R1 . . . Rm を計算する. Global step ではベクトル r を固定し, 頂点の位置 x につ. となる.. いて g を最小化する. g の x についての偏微分は. 4.3 離散化 本章では, 4.2 章で説明した式を離散化し, 有限要素から なるモデルに対して適応する方法を説明する.. であるので, 固定された r に対して g を最小化する x は以. ここで, モデルを m 個の要素（提案手法では四面体を用 いた.）に分割して離散化し, 定数項を無視すると, 式 (7) の 歪みエネルギーは以下のような行列を用いた形式で書き表 すことが出来る. m ∑ i=1. 2. µi ||Fi − Ri ||F = xT Lx − 2xT Jr. (9). ル r ∈ R9m は, R = (vec (R1 ) . . . vec (Rm )) であり, 行列. L ∈ R3n×3n と J ∈ R3n×9m はそれぞれ以下の式 (10), (11) のように定義される.. (m ∑. ˜ −1 ˜ −T µi AT i Mi Mi Ai. (15). 化を行い, 最終的な頂点の位置 x を求める.. J=. (m ∑. ˜ −T Ai µi M i. 本研究は 4.4 章で説明した, Local/Global approach の. Local step において回転行列を要素間で補間することに よって計算量を削減し, シミュレーションの高速化を行っ た. 本章ではその具体的な方法について述べる. まず, 5.1 章において特異値分解を行う代表要素のサンプリング方. ) ⊗ I3. (10). i=1. 法について述べる. 次に, 5.2 章において, クォータニオン ブレンディングを用いた回転行列の要素間補間の方法に. ) ⊗ I3. (11). i=1. Ai ∈ R. ( )−1 ( ) x = M + 2h2 L My + h2 (2Jr − f ext ). 5. 提案手法. トルである. ただし, 全頂点数を n とした. また, ベクト. 3×3n. 下の式 (15) によって与えられる.. 以上の 2 つのステップを繰り返し行うことにより g の最小. ここで, x ∈ R3n は弾性体を構成する全頂点の位置ベク. L=. ( ) ( ) ∂g = M + 2h2 L x − My + h2 (2Jr − f ext ) (14) ∂x. は i 番 目 の 要 素 の イ ン デ ッ ク ス 行 列 で,. ついて述べる. 最後に, 5.3 章において, 提案手法を用いて. Local/Global approach により最適化を行う方法について 述べる.. A1i1 = A2i2 = A3i3 = 1, A1i4 = A2i4 = A3i4 = −1 であり, ˜ i ∈ R3×3 は外力が 残りの成分はすべてゼロである. また, M. 5.1 サンプリング 特異値分解を行う代表要素を決定するために, 各要素の. 作用していないときの, i 番目の要素を構成する頂点の位置 ˜ i = (˜ ˜1, . . . , v ˜ 4 とすると, M ˜4, v ˜2 − v ˜4, v ˜3 − v ˜4) をv v1 − v. 重心に対して Schlomer ら [12] の farthest points を用いた. と表せる. また, I3 ∈ R. サンプリングを行った. その結果, 得られた l 個の代表要素. 3×3. は 3 次の単位行列であり, ⊗ は. ⓒ 2015 Information Processing Society of Japan. 3.
(4) Vol.2015-CG-160 No.8 2015/8/29. 情報処理学会研究報告 IPSJ SIG Technical Report. 素法によってポアソン方程式をディリクレ境界条件のもと で解き, 得られた解を用いて補間の重み w を決定した. まず, モデルの頂点座標 v1 , . . . , vn を節点として異なる 境界条件のもとでポアソン方程式を l 回解き, 解となるベ クトル u(1) . . . u(l) を得る. ただし, u(i) ∈ R3n は, 以下の ディリクレ境界条件を課したときのポアソン方程式の解で ある. (i) uj. 図 2 farthest points によるサンプリング. 1 j ∈ e のとき si = 0 j ∈ esk ∧ k ̸= i のとき. (18). Fig. 2 Farthest point sampling. ここで, j ∈ esi は j が要素 esi を構成する頂点のインデッ クスに含まれていることを表している. 更に, これらの解. ˜ (1) . . . u ˜ (l) のそれぞれに対して, 各要素ごとの平均をとり, u を得る. (i). このとき, 要素 ei に対する代表要素 esj の重み wj. を式. (19) によって定めた. u ˜i −1 ) l ( ∑ (k) u ˜i −1 (j). (i). wj = 図 3 ポアソン方程式から計算した重み. (19). k=1. Fig. 3 Blending weights derived from. また, 本手法では式 (16) の計算を高速化するために重み. solutions of Poisson equation. に対して閾値 wthreshold を設定し, 重みの値が閾値以下の. を es1 , . . . esl と表すものとする. また, サンプリングを行. 代表要素は補間される要素の回転行列への寄与が十分小さ. う要素の割合についてはユーザーが指定することが出来る. いとして計算に加えなかった. また, wthreshold の値につい. パラメータであるが, 本稿ではモデルの形状の複雑さに応. てはユーザーが決めることの出来るパラメータであるが,. じて, 全要素の 5% ∼ 15% をサンプリングし, 代表要素と. 本稿では重みの平均値と最小値の中間値とした.. した.. 5.3 最適化 本章では, 提案手法を適用した Local/Global approach. 5.2 回転行列の要素間補間 代表要素に対して計算された回転行列と後述の重み w か. により, 目的関数 g の最適化を行う方法を述べる. 始めに, 事前計算として 5.1 章と 5.2 章において述べた代. ら全要素の回転行列を補間する方法を説明する. サンプリングされた代表要素 es1 , . . . esl に対して, 特異 値分解を行って得られた回転行列 Rs0 , . . . , Rsl をクォー. 表要素 es0 , . . . , esl のサンプリングと重み w の計算を行っ ておく.. ˆ s0 , . . . , q ˆ sl とする. このとき, タニオンに変換したものを, q 要素 ei の回転行列を以下の式 (16) に従って, クォータニオ ンブレンディングによって補間する.. ˆi = q. l ∑. (i). ˆ sj (−1)sj wj q. (16). j=1. ここで, sj は最短経路による補完を行うための整数であり. Kavan ら [8] の手法に従って, 0 ⟨ˆ ˆ sj ⟩ ≥ 0 のとき qs1 , q sj = 1 ⟨ˆ ˆ sj ⟩ < 0 のとき qs1 , q. (17). のように定めるものとする. 式 (16) の計算を代表要素以外 のすべての要素に対して行い, 全要素の回転行列を得る. 次に, 要素間補間の際に用いたクォータニオンブレンディ ングの重み w の求めかたを説明する. 本研究では, 有限要. ⓒ 2015 Information Processing Society of Japan. 図 4 提案手法と Local/Global Optimization の適応結果. Fig. 4 Compare the results. 4.
(5) Vol.2015-CG-160 No.8 2015/8/29. 情報処理学会研究報告 IPSJ SIG Technical Report. 図 5 エネルギー関数の収束. Fig. 5 Our optimization converges.. Local step では, まず, 代表要素 es0 , . . . , esl のそれぞれ. 行列の要素間補間によって, 共回転系弾性体のシミュレー. T に対して特異値分解 Fi = Ui Σi Qi を行い, Ri = QT i Ui. ションを高速化する手法を提案した. これにより, 従来の. によって代表要素の回転行列 R1 . . . Rl を計算する. 次に,. Local/Global approach による最適化と比較して, 各イテ. ˆ s0 , . . . , q ˆ sl に変換し, 式 求めた回転行列をクォータニオン q. レーションに必要な時間を短縮することが可能となり, 結. (16) を用いて残りの要素の回転行列を補間する.. 果として解が収束するまでの時間の短縮を実現した.. Global step では 4.4 章と全く同様に式 (15) によって, g を最小化する x を求める.. 6. 結果 以上の手順によって, 3D モデルに対してシミュレーショ. しかし, 本手法で用いたクォータニオンブレンディングに よる回転行列の要素間補間はあくまで近似計算であるため, 最終的な解の収束点が Newton 法や本来の Local/Global. approach による最適化の結果として得られるものと一致 することは保証されていない.. ンを行った. その結果と, 従来の Local/Global Optimiza-. 今後は代表要素のサンプリングの手法や重み w の計算方. tion を適応したシミュレーションの結果を比較したものを. 法について, より適切な方法を検討していきたい. 特に, 本. 図 5.2 に示す. この結果から, 提案手法を用いても従来の. 手法では重み w をポアソン方程式の解をもとに決定した. Local/Global Optimization を適応した場合とほぼ同じ結. が, 今後は Jacobson ら [7] や Zhang ら [14] の手法を用いた. 果を得ることが可能であることが確認できる.. 場合の結果との比較を行う予定である.. 次に, 2375 個の四面体から構成される, 四角柱のモデル を用いて提案手法と従来の Local/Global Optimization の. 参考文献. エネルギーの相対誤差の減少の速度を比較した結果を図 5. [1]. に示す. ただし, 提案手法の括弧内の値はサンプリングを 行う要素の割合を表している. また, 相対誤差の真値とし てはそれぞれの手法において反復回数を 1000 としたとき の値を用いた. 図 5 から提案手法を用いた場合, 従来の. [2]. Local/Global Optimization を用いた場合よりも収束に要 する時間を短縮出来ることが確認出来る. また, 収束に要する時間の高速化の度合いははサンプ. [3]. リングを行う要素の割合に応じて変化し, サンプリング を行う要素の割合を増やしていくことによって, 従来の. Local/Global Optimization の結果に漸近していくことが. [4]. 分かる. このことから, サンプリングを行う要素の割合を 変化させることで速度と正確さのどちらにどの程度重きを おくか調整をすることが出来る.. [5]. 一方で, 解が収束するまでに要するイテレーションの回 数は提案手法と従来手法で大きな差は見られなかった.. 7. まとめと今後の課題 本稿では, クォータニオンブレンディングを用いた回転. ⓒ 2015 Information Processing Society of Japan. [6]. Baraff, D. and Witkin, A.: Large Steps in Cloth Simulation, Proceedings of the 25th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH ’98, New York, NY, USA, ACM, pp. 43–54 (online), 1998. Barbiˇc, J. and James, D. L.: Real-Time Subspace Integration for St. Venant-Kirchhoff Deformable Models, ACM SIGGRAPH 2005 Papers, SIGGRAPH ’05, New York, NY, USA, ACM, pp. 982–990 (online), 2005. Gao, M., Mitchell, N. and Sifakis, E.: Steklov-Poincar Skinning, Eurographics/ ACM SIGGRAPH Symposium on Computer Animation (Koltun, V. and Sifakis, E., eds.), The Eurographics Association, (online), 2014. Georgii, J. and Westermann, R.: Corotated Finite Elements Made Fast and Stable, Workshop in Virtual Reality Interactions and Physical Simulation ”VRIPHYS” (2008) (Faure, F. and Teschner, M., eds.), The Eurographics Association, (online), 2008. Hecht, F., Lee, Y. J., Shewchuk, J. R. and O’Brien, J. F.: Updated Sparse Cholesky Factors for Corotational Elastodynamics, ACM Trans. Graph., Vol. 31, No. 5, pp. 123:1–123:13 (online), 2012. Jacobson, A., Baran, I., Kavan, L., Popovi´c, J. and Sorkine, O.: Fast Automatic Skinning Transformations, ACM Trans. Graph., Vol. 31, No. 4, pp. 77:1–77:10 (on-. 5.
(6) 情報処理学会研究報告 IPSJ SIG Technical Report. [7]. [8]. [9]. [10]. [11]. [12]. [13]. [14]. Vol.2015-CG-160 No.8 2015/8/29. line), 2012. Jacobson, A., Baran, I., Popovi´c, J. and Sorkine, O.: Bounded Biharmonic Weights for Real-time Deformation, ACM Trans. Graph., Vol. 30, No. 4, pp. 78:1–78:8 (online), 2011. ˇ ara, J. and O’Sullivan, C.: GeKavan, L., Collins, S., Z´ ometric Skinning with Approximate Dual Quaternion Blending, ACM Trans. Graph., Vol. 27, No. 4, pp. 105:1– 105:23 (online), 2008. Liu, T., Bargteil, A. W., O’Brien, J. F. and Kavan, L.: Fast Simulation of Mass-spring Systems, ACM Trans. Graph., Vol. 32, No. 6, pp. 214:1–214:7 (online), 2013. Martin, S., Thomaszewski, B., Grinspun, E. and Gross, M.: Example-based Elastic Materials, ACM Trans. Graph., Vol. 30, No. 4, pp. 72:1–72:8 (online), 2011. McAdams, A., Zhu, Y., Selle, A., Empey, M., Tamstorf, R., Teran, J. and Sifakis, E.: Efficient Elasticity for Character Skinning with Contact and Collisions, ACM Trans. Graph., Vol. 30, No. 4, pp. 37:1–37:12 (online), 2011. Schl¨omer, T., Heck, D. and Deussen, O.: Farthest-point Optimized Point Sets with Maximized Minimum Distance, Proceedings of the ACM SIGGRAPH Symposium on High Performance Graphics, HPG ’11, New York, NY, USA, ACM, pp. 135–142 (online), 2011. Sorkine, O. and Alexa, M.: As-rigid-as-possible Surface Modeling, Proceedings of the Fifth Eurographics Symposium on Geometry Processing, SGP ’07, Aire-la-Ville, Switzerland, Switzerland, Eurographics Association, pp. 109–116 (online), 2007. Zhang, J., Deng, B., Liu, Z., Patan`e, G., Bouaziz, S., Hormann, K. and Liu, L.: Local Barycentric Coordinates, ACM Trans. Graph., Vol. 33, No. 6, pp. 188:1– 188:12 (online), 2014.. ⓒ 2015 Information Processing Society of Japan. 6.
(7)
図
関連したドキュメント
Trichoderma reesei cellobiohydrolase I (TrCel7A) molecules were observed to slide unidirectionally along the crystalline cellulose surface, and the catalytic domain without
In immunostaining of cytokeratin using monoclonal antibodies, the gold particles were scattered in the cytoplasm of the hepatocytes and biliary epithelial cells
Lemma4.1.. This is not true if f is not positively homogeneous as the following example shows.. Let f be positively homogeneous. We shall give an example later to show that
Keywords Poset · Rational function identities · Valuation of cones · Lattice points · Affine semigroup ring · Hilbert series · Total residue · Root system · Weight lattice..
We show that a discrete fixed point theorem of Eilenberg is equivalent to the restriction of the contraction principle to the class of non-Archimedean bounded metric spaces.. We
In particular, we show that the q-heat polynomials and the q-associated functions are closely related to the discrete q-Hermite I polynomials and the discrete q-Hermite II
The oscillations of the diffusion coefficient along the edges of a metric graph induce internal singularities in the global system which, together with the high complexity of
So far, most spectral and analytic properties mirror of M Z 0 those of periodic Schr¨odinger operators, but there are two important differences: (i) M 0 is not bounded from below
