GPGPUを活用した実時間脳腫瘍摘出シミュレーション

䢵䣃䢳䢯䣄䢲䢶
GPGPU を活用した実時間脳腫瘍摘出シミュレーション
Real-time Simulation of Brain Tumor Resection Using GPGPU
○ 佐瀬一弥（北大）
正近野敦（北大）
正辻田哲平（東北大／カーネギーメロン大）
福原洸（東北大）陳暁帥（北大）正小水内俊介（北大）
Kazuya SASE, Atsushi KONNO, Hokkaido University,
Teppei TSUJITA, Tohoku University/Carnegie Mellon University
Akira FUKUHARA, Tohoku University,
XiaoShuai CHEN, Shunsuke KOMIZUNAI, Hokkaido University,
[email protected]
This paper proposes a parallelization method for brain tumor resection simulation using General Purpose Graphics Processing Unit (GPGPU).
The simulation is performed by using corotational FEM for fracture handling. However, the calculating cost is too high to execute in real-time.
The bottlenecks of the simulation are solving linear system of equation and assembly procedure of global stiffness matrix. In order to achieve
real-time simulation, we apply the large scale parallelization using GPGPU to the bottlenecks. The proposed method optimizes memory
consumptions and calculations by taking account of symmetry property of the stiffness matrices. This paper shows that the proposed method
enables the simulation to run in around one and half times faster than the case of multi-core CPU parallelization.
Key Words: Surgery simulation, Soft material, GPGPU, Finite element method
1. 緒言
Start
医療技術の高度化によって医師に要求される技能の水準が高
まっており，効率的な手術手技の訓練環境が求められている．
そこで近年，バーチャルリアリティ技術を応用した手術シミュ
レータが注目されている．手術シミュレータは力覚提示装置
をユーザが操作し，計算機によって仮想空間内に構築された
臓器を触った感覚を提示することができる．これによってあ
たかも実際に手術をしているような臨場感を再現することが
期待されており，新しい手術手技訓練環境としての応用研究
が盛んに行われている．脳外科手術の分野では止血や吸引な
どの操作が可能な実時間シミュレータが報告されている [1]．
しかしながら，任意の箇所で剥離操作を行うことはできない．
剥離操作のような破壊を伴う手術シミュレーションは計算コ
ストが高く，さらに不安定化しやすいため手術シミュレー
ションの分野において大きな課題のひとつとなっている．
本研究では脳腫瘍摘出手術を対象とし，幾何学的非線形
性を考慮した有限要素法 (Finite Element Method, FEM) である
Corotational FEM を用い，要素削除に基づく破壊表現を行うこ
とで実時間脳腫瘍摘出手術シミュレーションを実現することを
目的とする．Corotational FEM は線形 FEM やばね質点モデルよ
りも計算コストが高く実時間計算のために並列処理を有効に
用いる必要性が高い．本論文では，全体剛性マトリクス組み
立て，および，連立一次方程式求解の General Purpose Graphics
Processing Unit(GPGPU) 実装について述べ，節点数 10000 規模
である脳腫瘍患部の有限要素モデルに対する実時間計算性能
を明らかにする．
2. 生体組織の変形・破壊計算手法
Collsion detection
Aspirator
(A) Calculate Ke’
(B) Finite element assembly
Apply boundary condition
(C) Solve Avfi+1 = b
Ohter functions
Fig. 1 Simulation of brain tumor resection
Fig. 2 Flowchart
で GPGPU による並列化を実装した部分を示している．以後の
図中に記載した (A)，(B)，(C) の表記は図 2 中の処理と対応して
いることを示す．それぞれの具体的な処理を本章で述べる．
2.2. 要素剛性マトリクスの作成
Corotational FEM では局所的な回転に対する誤差を保障する
ため，要素剛性マトリクスを剛体回転成分を取り除いた座標
系に変換する．要素剛性方程式は次式のように表される．
(
)
f e = Re K e RTe xe − xe0
(1)
ここでf e ∈ R12 , xe ∈ R12 , x0e ∈ R12 , K e ∈ R12×12 はそれぞれ
要素外力ベクトル，要素位置ベクトル，要素初期位置ベクトル，
要素剛性マトリクス，Re ∈ R12×12 は次式により与えられる．
2.1. 概要
本研究では腫瘍摘出過程の剛体回転を含む大変形を計算す
るため，幾何学的非線形性を考慮した変形計算を行う必要が
ある．そこで，生体組織を四面体要素の集合により離散化し，
Corotational FEM [2] を適用し，破壊表現に応力指標に基づく
要素削除を採用する．脳腫瘍摘出シミュレーションの例を図 1
に，フローチャートを図 2 に示す．図中の (A), (B), (C) は本研究
Precomputing (Ke, ....)
Tumor
Re = blockdiag[R R R R]
(2)
ここでR ∈ R3×3 は四面体要素の剛体回転成分を表す回転行列で
あり，四面体要素の変形から剛体回転成分を抽出することで求
められる．式 (1) は連立一次方程式であり，次のように書ける．
f e = K ′e xe + f ′0e
䣐䣱䢰䢢䢳䢶䢯䢴䢢䣒䣴䣱䣥䣧䣧䣦䣫䣰䣩䣵䢢䣱䣨䢢䣶䣪䣧䢢䢴䢲䢳䢶䢢䣌䣕䣏䣇䢢䣅䣱䣰䣨䣧䣴䣧䣰䣥䣧䢢䣱䣰䢢䣔䣱䣤䣱䣶䣫䣥䣵䢢䣣䣰䣦䢢䣏䣧䣥䣪䣣䣶䣴䣱䣰䣫䣥䣵䢮䢢䣖䣱䣻䣣䣯䣣䢮䢢䣌䣣䣲䣣䣰䢮䢢䣏䣣䣻䢢䢴䢷䢯䢴䢻䢮䢢䢴䢲䢳䢶
䢵䣃䢳䢯䣄䢲䢶䢪䢳䢫
(3)
ある．運動方程式は次式となる．
Local index 0
Global index n0
Local index 1
Global index n1
e th element
¨ + C x˙ + K ′ x + f ′0 = f ext
Mx
Connectivity matrix E
E (e, 0) = n0
E (e, 1) = n1
E (e, 2) = n2
E (e, 3) = n3
Local index 3
Global index n3
˙ ，x
¨ はそれぞれxの時間に関する一階微分, 二階微分で
ここで，x
ある．M ∈ R3Nnode ×3Nnode , C ∈ R3Nnode ×3Nnode , f ext ∈ R3Nnode
はそれぞれ質量マトリクス，減衰マトリクス，外力ベクトル
である．運動方程式の境界条件を反映し，時間発展を陰解法
により計算する．この計算は次の連立一次方程式に帰着する．
Local index 2
Global index n2
Fig. 3 Connectivity matrix
0
1
2
Av i+1
=b
f
E(e,3) E(e,1) E(e,2) E(e,0)
E(e,3)
3
E(e,1)
2
E(e,2)
3
Ke
E(e,0)
Keex
ex
Ke ( E(e, 0), E(e, 1) ) = Ke(0,1)
Fig. 4 Exntension of an element stiffness matrix to the dimension of
the global stiffness matrix
ここで，f ′0e ∈ R12 は初期位置によって定義される定数ベクト
ル，外力オフセットである．K ′e ∈ R12×12 は剛体回転成分を
補正した要素剛性行列であり次式により得られる．
K ′e = Re K e RTe
(4)
式 (4) は図 2 中の (A) に該当する．
2.3. 全体剛性マトリクスの組み立て
得られた Ke ′ を全体剛性マトリクスの次元に拡張，総和処
理し，四面体メッシュの全自由度の情報を格納する全体剛性
マトリクス K ′ ∈ R3Nnode ×3Nnode を組み立てる．ここで，Nnode
は節点数である．Ke′ を K ′ の次元に拡張するにはコネクティ
ビティマトリクス E(e, i) を用いる．E(e, i) は四面体要素 e の節
点 i (i = 0, 1, 2, 3) をメッシュ全体における節点番号に変換する
ex
マトリクスである（図 3）．したマトリクス Ke′ は，E(e, i) を
用いて次式のように表される．
Ke′ (E(e, i), E(e, j))block = Ke′ (i, j)block
ex
(5)
ここで，(.)block は3 × 3のブロックマトリクスへのインデックスで
ex
あること示す．式(5)の例としてKe′ (0, 1)をKe′ (E(e, 0), E(e, 1))
ex
′
へ格納する場合の模式図を図 4 に示す．K は Ke′ の総和処理
を行うことによって得られる．
K′ =
Nele −1
∑
Ke′
ex
(6)
e=0
ここで，Nele は四面体要素数である．式 (5), (6) は図 2 中の (B) に
該当する．
2.4. 有限要素方程式の求解
全体剛性方程式は次式となる．
f = K ′ x + f ′0
(9)
ここで，v i+1
∈ R3Nfreenodes は求める外力既知節点の速度ベク
f
トルを表す．A ∈ R3Nfreenodes ×3Nfreenodes ，b ∈ R3Nfreenodes は境界
条件と過去のステップの情報から得られる既知な変数である．
ただし，Nfreenodes は外力既知節点の数である．式 (9) は図 2 中
の (C) に該当する．
0
1
(8)
(7)
ここで，f ∈ R3Nnode , x ∈ R3Nnode , f ′0 ∈ R3Nnode , はそれぞれ全
体外力ベクトル，全体位置ベクトル，全体外力オフセットで
3. GPGPU 実装
3.1. システム概要
数値計算用コンピュータとして，マルチコア CPU (Intel
Core i7 3960X 3.3 [GHz], 6 cores) を 4.8 [GHz] にオーバークロッ
クして利用可能な CIARA 社製ワークステーション KRONOS
S810R を用いた．また，計算用 GPU として，NVIDIA 社製 Tesla
K20c (2496 cores, GPU Clock rate 706 [MHz], Kepler architecture)
を搭載した．数値計算にはマルチコア CPU，および，GPU に
よる並列処理をともに活用した．マルチコア CPU における並
列処理には OpenMP を，GPGPU 実装には NVIDIA 社が提供する
CPU 向け統合開発環境 CUDA を用いた．ワークステーション
は LAN を介してハプティックデバイスと接続された制御コン
ピュータと通信することで力覚提示を可能とした．
図 5 にホスト側のコンピュータとデバイス側の GPU による処
理分担と通信のやり取りを示す．実時間ループの前処理とし
て線形 FEM における要素剛性マトリクスの作成と，後述する
縮約配列を作成し GPU に送信する．実時間ループの中で
CUDA カーネルを呼び出し，実行する．現状の実装では GPU
と CPU の並列実行を行っていない．
3.2. 要素剛性マトリクスの作成
式 (4) より，Ke ′ を得るためには回転行列 Re と線形 FEM の要
素剛性マトリクス Ke が必要である．Ke はシミュレーション
中に変化することはないため，前処理としてあらかじめ計算
し，GPU に送信した．一方，回転行列 R は毎ループで変化す
るため更新が必要である．R の更新はマルチコア CPU によっ
て計算し，GPU に送信する．GPU における演算ではスレッド
あたりにひとつの四面体要素を割り当て，式 (4) を実行する．
この処理においては Ke および Ke′ が対称行列であることを
利用して，メモリ容量，通信量，演算数の削減を行う．Ke は
シミュレーションの前処理において計算され，対角要素を
含む行列の下部を一次元化して送信する．一次元化した
データにおいて (i, j) 成分にアクセスするにはインデックス
(i ∗ (i + 1)/2 + j) を参照する．式 (4) では出力 Ke′ も対称行列と
なることから対角成分を含む下部のみについて演算を行えば
十分であり，これによって演算数を半減できる．
また，式 (4) は 12 × 12 の行列の積となっているが，式 (2) より
Re が対角にブロック行列をもつのみであるため，Re を構築す
る必要はなく，3 × 3 の R を送信し，演算に用いることで演算
数とメモリアクセスの削減を行うことができる．
䣐䣱䢰䢢䢳䢶䢯䢴䢢䣒䣴䣱䣥䣧䣧䣦䣫䣰䣩䣵䢢䣱䣨䢢䣶䣪䣧䢢䢴䢲䢳䢶䢢䣌䣕䣏䣇䢢䣅䣱䣰䣨䣧䣴䣧䣰䣥䣧䢢䣱䣰䢢䣔䣱䣤䣱䣶䣫䣥䣵䢢䣣䣰䣦䢢䣏䣧䣥䣪䣣䣶䣴䣱䣰䣫䣥䣵䢮䢢䣖䣱䣻䣣䣯䣣䢮䢢䣌䣣䣲䣣䣰䢮䢢䣏䣣䣻䢢䢴䢷䢯䢴䢻䢮䢢䢴䢲䢳䢶
䢵䣃䢳䢯䣄䢲䢶䢪䢴䢫
Element data
Host
Device
Reduction array
Precomputing
Reduction array
K’ (csrValues)
Ke
Nonzero entry / thread
Device
Host
Rigid config.
Re
K(CSR)
Real-time
loop
Copy
K’ (csrValues)
Update every simulation loop
K’ (csrRowOffsets)
Invariant during a simulation
K’ (csrColIndices)
Invariant during a simulation
(A) Re Ke ReT
(B) Reduction
Apply
boundary condition
A, b
Fig. 6 Data structure of reduction
(C) Solve Av = b
0
v
1
2
0
3
1
n j
3
E(e, 1) = n
0
Fig. 5 Overview of the global stiffnes matrix assemble process
1
2
n
1
3
0
3
0
1
2
3
0
i
1
2
し，全体剛性マトリクスの組み立てに用いる要素データとし
て格納する．
3
E(e, 3) = i, E(e, 1) = j
1
以上のようにして得られた Ke′ はホスト側に返すことなく
GPU 側で保持する．Ke′ はすべての要素にわたって一次元化
2
2
2
Other function
1
0
0
Global matrix
2
3
3
E(e, 2) = n
Element stifness matices
of tetras which share node n
E(e, 1) = i, E(e, 0) = j
Element stifness matices of tetras
which share an edge (i, j), (i > j)
Fig. 7 Global stiffness matrix and mesh information
3.3. 全体剛性マトリクスの組み立て
全体剛性マトリクスは疎行列であることから，一般的に疎
行列形式で格納される．しかし，疎行列形式ではランダムア
クセスにかかるコストが高い．また，成分の追加を行うため
にメモリ領域の確保が必要になるとさらに効率が低下する．
この処理は有限要素法において連立一次方程式の求解に次い
で計算コストが高く，ボトルネックとなることが多い [3]．
Cecka らはあらかじめ総和処理の入力となる要素データの
インデックスと総和結果の格納先のインデックスを対応付け
るマップ (縮約配列) を作成し，疎行列演算で必要となるイン
デックスの検索やエントリのソートの処理を排除する方法を
提案した [3]．本研究では同手法を採用し Corotaional FEM の
データ構造に適用した．Cecka らの手法はメッシュ構造がシ
ミュレーションにわたって変化しないことを要求する．本研
究では破壊シミュレーションにより要素削除が行われるが，
節点同士の接続関係は不変である．また，四面体が除去され
た場合，要素剛性マトリクス作成部で R = O とすることで要
素剛性マトリクスの寄与を無効化することができるため，同
手法を適用可能である．
図 6 に総和処理におけるデータ構造を示す．要素データに要
素剛性マトリクスの値が格納されている．本実装では K ′ を
CSR(Compressed Sparse Rows) 形式で格納する．メッシュ構造
が不変であれば K ′ の行オフセット配列，列インデックス配列
はシミュレーションにおいて不変である．したがって，K ′ の
値配列のみを更新すれば十分である．同手法では K ′ の値配列
の各エントリ (非ゼロエントリ) に対してスレッドを割り当て
る．非ゼロエントリのインデックスに対応する縮約配列にア
クセスし，総和の入力となる要素データを足し合わせ，最終
的に値配列の該当箇所に総和結果を格納する．このようにす
ることで，各スレッド間に依存関係がなく，アトミック演算
などの同期処理も不要である．
本研究ではさらに，剛性マトリクスの対称性を考慮するこ
とで必要リソースの削減を試みた．縮約配列を作るに当たり，
全体剛性マトリクスの各成分がどの要素剛性マトリクスと関
係しているかを明らかにする．図 7 にメッシュ構造と縮約が行
われる全体剛性マトリクスのインデックスとの対応を示す．
メッシュ構造からある節点 n を共有するすべての四面体要素を
抽出する．節点 n の四面体要素内におけるインデックスを i
(i = 0, 1, 2, 3)とすると四面体要素の(i, i)block 成分のブロックマ
トリクスが全体剛性マトリクスの対角に位置するブロックマト
リクスに対応する．対角に位置するブロックマトリクスは対称
行列であるから，縮約対象は 6 つの成分のみである．辺につい
ても同様であるが，辺に対応するブロックマトリクスは対称行
列ではなく 9 つの成分のすべてを縮約対象とする必要がある．
これらのメッシュ構造と縮約の対応関係をもとに縮約配
列を作成する（図 8）．縮約配列は非ゼロエントリの数だ
け存在する．非ゼロエントリ数はメッシュとの関係性から
6Nnode + 9Nedge である．Cecka ら [3] の方法に習い，ひとつの
縮約配列において最後尾に負の符号をつけた図中の灰色の成
分を配置し，この符号を反転させ縮約先のインデックスとす
る．図中の白色の成分は縮約元の要素データのインデックス
を指す．デバイス側では縮約配列を一次元的に配置する．縮
約配列は個々に長さが異なるため，各配列の先頭インデック
スを格納する配列が必要である．これらのデータ構造を用い
䣐䣱䢰䢢䢳䢶䢯䢴䢢䣒䣴䣱䣥䣧䣧䣦䣫䣰䣩䣵䢢䣱䣨䢢䣶䣪䣧䢢䢴䢲䢳䢶䢢䣌䣕䣏䣇䢢䣅䣱䣰䣨䣧䣴䣧䣰䣥䣧䢢䣱䣰䢢䣔䣱䣤䣱䣶䣫䣥䣵䢢䣣䣰䣦䢢䣏䣧䣥䣪䣣䣶䣴䣱䣰䣫䣥䣵䢮䢢䣖䣱䣻䣣䣯䣣䢮䢢䣌䣣䣲䣣䣰䢮䢢䣏䣣䣻䢢䢴䢷䢯䢴䢻䢮䢢䢴䢲䢳䢶
䢵䣃䢳䢯䣄䢲䢶䢪䢵䢫
160
Reduction array for a nonzero entry
Node 0
Index to csrValues
(negative signed)
9
Number of tetras
which share
the edge i
6
Number of tetras
which share
the node i
Node i
Edge 0
5.4
140 5.1
120
Time [ms]
Indices to element data
(positive signed)
3.3
5.3
60.1
100
30.6
80
32.9
60
40
33.0
44.8
33.7
CPU
(5 threads)
GPU
20
Fracture handling
Calculate stress
(C) Solve Av = b
Apply boundary conditions
Calculate R, (A) Ke, (B) K’
0
Edge j
9
6
Fig. 9 Average calculation time for a simulation loop
Edge Nedge-1
9
50
6
Host
Device
Vectorize
Vectorize
Time [ms]
Node Nnode-1
Reduction array
on the device memory
Indices offsets for
reduction array
40
9.4
30
9.1
Calculate f0
20
10
0 4 8
0
Fig. 8 Data structure of reduction array
8.7
4.1
15.0
Convert K’ to full matrix
Convert K’ to Eigen sparse matrix
9.0
Calculate (A) K’e and (B) K’
3.9
4.1
Calculate R
6.0
6.0
CPU
(5 threads)
GPU
Fig. 10 Calculation time of calculation of R, Ke , K ′
て以下のようにして縮約処理を行う．
1. スレッドは先頭インデックス配列を参照して縮約配列に
アクセスする．
2. 縮約配列を参照して総和対称の要素データにアクセスする．
3. 縮約配列の値が負になるまで要素データの収集を繰り返す．
4. 縮約配列の値が負になったら負号を取り除いて，縮約先
インデックスを取得し，総和の結果をデバイスメモリに
書き込む．
解であり，51 [%] まで削減することができた．図 10 に要素剛
性マトリクスの作成および剛性マトリクスの組み立て部の
詳細な計算時間を示す．このうち GPGPU 実装を行ったのは
“Calculate Ke′ and K ′ ” の部分であり，計算時間を 27 [%] にまで
削減することができた．
以上のように，1 ループ 105 [ms] のほぼリアルタイムシミュ
レーションを実現した．さらなる高速化のために境界条件の適
用部の GPGPU 実装を行うことが効果的であると考えられる．
5. 結言
3.4. 連立一次方程式の求解
連立一次方程式の求解には共役勾配 (Conjugate Gradient, CG)
法を用いた．CG 法で必要となる疎行列と疎ベクトルの積は
CUDA の疎行列ライブラリ CUSPARSE を用いた．
4. 結果と考察
第 3 章で述べた GPGPU 実装を行い，脳腫瘍患部を模擬した
有限要素モデル（節点数 9184，四面体要素数 34050）に対し
腫瘍摘出シミュレーションを行った．図 9 に 1 ループあたりの
計算時間 (100 ループ平均) を示す．計算時間は処理ごとに分け
て表示している．“Calculate R, Ke , K ′ ” は四面体の回転成分抽
出から，全体剛性マトリクス組み立てまでの処理である．
“Apply boundary condition” は境界条件に基づく剛性マトリクス
分解処理であり，GPGPU 実装の対象としなかった．“Solve
Av = b” は連立一次方程式の求解である．“Calculate stress” は
応力計算，“Fracture handling” は破壊処理である．比較のため
マルチコア CPU のみの場合と GPGPU を活用した場合の結果を
示した．全体として GPGPU 活用により計算時間を 71 [%] にま
で削減した．特に効果が大きかったのは連立一次方程式の求
本論文では，剛性マトリクスの対称性を考慮した縮約配列に
よる全体剛性マトリクス組み立て，および，連立一次方程式の
CG 法による求解の GPGPU 実装について述べた．脳腫瘍患部の
有限要素モデルに対して脳腫瘍摘出シミュレーションを行い，
1 ループの計算にマルチコア CPU で 148 [ms] かかっていた処理
を，GPGPU 実装により 105 [ms] まで高速化することができた．
謝辞
本研究は，日本学術振興会最先端・次世代研究開発支援プ
ログラム（課題番号: LR003）の助成を受けて行われた．
文献
[1] Delorme, S. et al., “NeuroTouch: A Physics-Based Virtual Simulator for Cranial
Microneurosurgery Training,” Neurosurgery, vol. 71, pp. 132–142, 2012.
[2] M¨
uller, M. et al., “Interactive Virtual Materials,” Proceedings of 2004 Graphics
Interface Conference, pp. 239–246, 2004.
[3] Cecka, C. et al., “Application of Assembly of Finite Element Methods on Graphics
Processors for Real-time Elastodynamics,” GPU Computing Gems Jade Edition,
pp. 187–205, 2012.
䣐䣱䢰䢢䢳䢶䢯䢴䢢䣒䣴䣱䣥䣧䣧䣦䣫䣰䣩䣵䢢䣱䣨䢢䣶䣪䣧䢢䢴䢲䢳䢶䢢䣌䣕䣏䣇䢢䣅䣱䣰䣨䣧䣴䣧䣰䣥䣧䢢䣱䣰䢢䣔䣱䣤䣱䣶䣫䣥䣵䢢䣣䣰䣦䢢䣏䣧䣥䣪䣣䣶䣴䣱䣰䣫䣥䣵䢮䢢䣖䣱䣻䣣䣯䣣䢮䢢䣌䣣䣲䣣䣰䢮䢢䣏䣣䣻䢢䢴䢷䢯䢴䢻䢮䢢䢴䢲䢳䢶
䢵䣃䢳䢯䣄䢲䢶䢪䢶䢫

Download Report