計測自動制御学会東北支部 第 279 回研究集会 (2012.3.4)
資料番号 279-5

追尾放射線治療のための高速な肺腫瘍位置変動計測法
High-speed Measurement of Lung Tumor Motion for Tumor
Following Radiotherapy○澁澤直樹 ∗ ，市地慶 ∗∗ ，張暁勇 ∗∗∗ ，本間経康 ∗∗∗ ，髙井良尋 † ，成田雄一郎 † ，石原恵太 ∗∗ 酒井正夫 ‡ ，阿部誠 ∗∗ ，杉田典大 ∗∗ ，吉澤誠 ∗∗∗○ Naoki Shibusawa∗ , Kei Ichiji∗∗ , Xiaoyong Zhang∗∗∗ , Noriyasu Homma∗∗∗ ,Yoshihiro Takai† , Yuichiro Narita† , Keita Ishihara∗∗ , Masao Sakai‡ ,Makoto Abe∗∗ , Norihiro Sugita∗∗ , Makoto Yoshizawa∗∗∗*東北大学工学部，**東北大学大学院工学研究科，***東北大学サイバーサイエンスセンター，† 弘前大学大学院医学研究科，‡ 東北大学教育情報基盤センター*School of Engineering, Tohoku Univ., **Graduate School of Engineering, Tohoku Univ.,***Cyberscience Center, Tohoku Univ., †Hirosaki Univ. Graduate School of Medicine,‡Center for Information Technology in Education, Tohoku Univ.
キーワード : 放射線治療 (radiotherapy)，肺腫瘍 (lung tumor)，マーカーレストラッキング (markerlesstracking)，テンプレートマッチング (template matching)，自己回帰モデル (autoregressive model)

連絡先 : 〒 980-8578 仙台市青葉区荒巻字青葉 6-3 電気・情報系
東北大学サイバーサイエンスセンター 先端情報技術研究部 吉澤・本間研究室
澁澤直樹，Tel.: (022)795-7130, Fax.: (022)795-7129, E-mail: shibusawa@yoshizawa.ecei.tohoku.ac.jp

1.

はじめに位置変動に対応した照射法として，治療ビー

近年，高齢化社会の進行に伴い，癌の罹患数が増加している 1) ．特に高齢者は体力の低下により手術を行うことが困難であることも多く，負担の少ない放射線治療への期待が高まっている．
放射線治療においては，副作用を避けるために健康な組織への照射を避け，治療効果を高めるために腫瘍だけに十分な線量を正確に照射することが求められる．しかし，肺腫瘍のように呼吸や拍動に起因して腫瘍位置が変動する場合にはその位置変動への対策が必要となる．

ムの照射位置自体を移動させて, 腫瘍の位置変動を追尾しながら連続的に照射する動体追尾照射が提唱されている 2) ．このような腫瘍位置変動を追尾する方法では，腫瘍位置をリアルタイムに計測し，その位置を目標値としたサーボ系で照射位置を制御する必要がある．
放射線治療装置の普及機に搭載されている X
線透視装置は Fig. 1 に示すような肺腫瘍位置を直接リアルタイムに撮像することが可能で，現在のところ最も正確なリアルタイム計測を実現–1–置の計測実験を行い，その有効性を実証した．

2.

関連研究

2.1

テンプレートマッチング (TM)TM とは，あらかじめ指定したテンプレートと最も類似度の高い対象画像領域を探索する方法である．類似度の評価指標として画像間の相Fig. 1

kV X 線透視画像例

関係数 (correlation coeﬃcient, CC) などがある．
対象画像の基準位置 (i, j) における画素の輝度できる可能性をもつモダリティである．一方，照値を I(i, j) とすると，サイズ M × N のテンプ

射位置の制御には，やはり普及機に広く搭載されレート T との CC は次式で定義される．ている dynamic Multi-leaf Collimator (dMLC)

CC(x, y)が担うことが可能である．

M−1 N−1∑∑

しかし，これらの装置を用いて動体追尾照射¯ (i, j) − T̄ )(I(x + i, y + j) − I)(Tを行う場合，計測から照射までの間にそれぞれi=0 j=0=M N σI σT

画像計測の計算処理や機械的な要因による時間(1)遅れが内在する．その内訳としては，計算処理に数百 ms, 機械的な処理に数十 ms の時間がかかるという報告もある 3) 4) ．システムの時間遅

ここで I¯，T̄ はそれぞれ対象画像，テンプレート

れが大きい場合，照射位置と腫瘍位置との間に画像，テンプレートの輝度値の標準偏差である．

ずれが生じ, 健康な組織への被曝が懸念される．

ある探索範囲 Ωに対して評価指標の値が最大とこのことは，動体追尾照射の普及機での実現のなる位置 (x∗ , y ∗ ) を対象の位置とする．すなわ妨げとなっている．一方，時間遅れの大部分を占ち，(x∗ , y ∗ ) = arg max CC(x, y) となる．の輝度値の平均であり，σI ，σT はそれぞれ対象(x,y)∈Ω

める画像計測の計算処理の時間を短縮すること

TM は，基本的に対象画像内を全探索し大域が出来れば健康な組織への被曝が軽減できると最適解を得る手法である．したがって，TM の考えられる．たとえば X 線画像計測で対象を追計算量は対象画像のサイズに依存する．具体的跡する目的でよく用いられる手法にテンプレーには対象画像の画素数を X × Y ，テンプレート

トマッチング (template matching, TM) がある
5) ．この手法では，対象が存在する可能性があの画素数をそれぞれ M × N としたとき，(X −M + 1)×(Y − N + 1) 回の評価指標の計算が行る領域を全探索するため，2.1 節で詳しく述べ

われ，サイズが大きいほど計算量も増大してしるように探索領域が広くなると計測時間が膨大まう．になる．
そこで本研究では，対象が存在する可能性がより高い領域のみに探索領域を絞り込むことで，

2.2

Diamond Search

チング法を提案する．提案法の核は，探索領域位置計測処理の高速化を図る先行研究として，
New Diamond Search (DS)6) という手法につい絞り込みを腫瘍位置変動予測に基づいて設定すて説明する．DS は動画圧縮などに用いられる高ることである．実際の臨床画像を用いた腫瘍位速なブロックマッチングアルゴリズムの一種で計測時間を短縮する新しい高速テンプレートマッ–2–

あり，Fig. 2 のような 2 種類の探索パターンを用いる．1 つ目の探索パターンは large diamondsearch pattern (LDSP) と呼ばれ，中心のブロックと周辺の 8 ブロックからなり，ダイヤモンド型の探索パターンを形成する．2 つ目の探索パター
ンは small diamond search pattern (SDSP) と呼ばれ，中心のブロックとその 4 近傍のブロックからなる．

(a) LDSP(b) SDSP

はじめに，動画の前フレームの対象位置を中Fig. 2心とする LDSP の 9 ブロックの評価関数を計算探索パターンする．このとき，評価関数が最大となるブロック

β が LDSP の中心ブロックでなければ，ブロック3.

提案手法

3.1

位置予測を用いた探索範囲限定法

β を中心とする新たな LDSP の評価関数を計算する．LDSP の中心ブロックの評価関数が最大となったら探索パターンを SDSP に変更し，評2.1 節で述べたように TM により位置探索を価関数を計算する．最後に，評価関数が最大と行う場合，計算量は対象画像の大きさに依存しなったブロックが最終的なマッチング位置となて増加する．画像の大きさによらず最終的にマッり，中心ブロックの初期位置とマッチング位置チするのは１箇所のみであるから，対象画像がの座標を結んだ直線を追跡領域の動きベクトル大きいほど無駄な計算が多いということでもあとする．

る．そこで，別な方法，すなわち位置変動予測

以上から明らかなように，DS では対象画像中の限られたブロックのみを評価し，それ以外により，探索範囲を絞り込んで高速化する手法を提案する．の領域を無視することで計算量を抑制している．

具体的には，ある時刻 t において対象画像に対しかし，ピクセル単位での位置計測に用いる場して TM を行う際，時刻 t − 1 までに得られた腫合には最終的にマッチしたブロック内を TM な瘍位置の計測時系列から 1 フレーム先の腫瘍位どによって探索する必要がある．ブロックサイ置を予測し，その予測された腫瘍位置周辺のみをズを M × N とすると，TM の評価関数の計算探索範囲とすることで計算量を減少させる．ここ

回数は M × N となる．つまり，肺腫瘍位置計測を行う際には１フレームごとにブロックサイで，推定位置時系列に自己回帰 (autoregressive,AR) モデル 7) を当てはめ，これによって得らズ分の計算が必要となるので，結局，計算量はれた予測値と予測分散から探索範囲を限定する

ブロックサイズに依存する．ブロックサイズを手法を提案する．

小さく設定することで計算量を減らすことは可能だが，小さくしすぎると腫瘍の全体形状やテ

3.1.1

自己回帰モデル（AR モデル）

クスチャ情報をマッチングに用いることができなくなり，計測誤差が増大してしまう可能性がある．

AR モデルとは時系列を表現するモデルの一つである．m 次の AR モデルは時系列 y1 , y2 · · ·が与えられたとき，次式で定義される．

yt =m∑
i=1–3–

ai yt−i + vt(2)ここで，ai は AR 係数，vt は時刻 t における平均 0，分散 σ 2 の正規分布に従う白色雑音である．
時刻 t までに得られた値より時刻 t + 1 の予測値

ŷt+1 は，以下のような式で表わされる．
ŷt+1 =m−1∑(8)222σ̂m= σ̂m−1{1 − (âmm) }(9)(iv) あらかじめ定めた最大次数 L まで (i) から(iii) を繰り返す．

ai+1 yt−i(3)i=0すなわち，AR モデルによって yt+1 を予測するには，AR 係数 ai を同定する必要がある．同定には，Burg 法 7) を用いた．

3.1.2m−1m−1mwt−m= wt−m− âmm vtBurg 法

また，AR モデルの次数は赤池情報量規準 (AkaikeInformation Criterion, AIC)7) により決定する．

AIC は次式によって計算される．
2AICm = T (log 2πσ̂m+ 1) + 2(m + 1） (10)最適な次数の AR モデルを用いるため各 AR モ

時系列 y１ , . . . , yT から AR モデルの同定を行う．まず，vtm−1 を時刻 t における m − 1 次の前向きの AR モデルの予測誤差，wtm−1 を m−1 次

デルの AIC，すなわち AIC0 , . . . , AICL を計算し，その中で最小の値をとる AR モデルを選択する．の後向きの AR モデルの予測誤差とする．yt の3.1.3

平均を 0 とすると vt0 = wt0 = yt ，n = 1, . . . , Tとおける．そのときの予測誤差分散は1∑ 2ytT

時刻 t + 1 において新たな画像が得られた際，Tσ̂02 =予測値に基づく探索範囲の設定

(4)

t=1

提案法では，まず推定位置時系列 yt−T +1 , . . . , ytを用いて時刻 t + 1 の腫瘍位置の予測値 ŷt+1 と予測分散 σ̂ 2 を求める．次に予測値 ŷt+1 を中心となる．
次数 m の AR モデルの係数は，以下の手順でとし, 探索範囲を± 2σ̂ とすることにより, 腫瘍位置計測の計算量を予測精度に応じて適応的に求めることができる．

減少させる．

(i) AR 係数 âmm を次式より求める．
2âmm =T∑4.

m−1
vtm−1 wt−mt=m+1T∑m−1 2(wt−m) +t=m+1T∑検証実験

位置計測性能と計算速度の検証のため，実際(vtm−1 )2の臨床画像を用いて，全探索，先験情報を用いt=m+1(5)て探索範囲を限定する ROI 探索，DS，提案法の 4 手法による肺腫瘍位置計測誤差と計算速度(ii) (i) で求めた âmm より，m 次 AR モデルのm他の係数 âm1 , . . . , âm−1 を次式より求める．
m−1m−1âm− âmj = âjm âm−j(6)の比較検証を行った．なお本研究では，以下の環境下で実験を行った．CPU: Intel Xeon CPU

E3-1225 V2 @ 3.20GHz，OS: Windows 7 Professional 64 bit，RAM: 16 GB, 4 GB modules(iii) m 次 AR モデルの前向き予測誤差 vtm ，後mおよび予測誤差分散
向き予測誤差 wt−mx4. PC3-12800 (800 MHz)，Software: MATLAB R2012b．

2 をそれぞれ以下の式より計算する．
σˆmm−1
vtm = vtm−1 − âmm wt−m(7)–4–

評価指標

4.1

4.3

実験結果

位置計測性能の指標として x 座標，y 座標そFig. 4 に各手法の計測誤差 RMSE の一例を，れぞれに対して平均二乗平方根誤差 (root mean

Table 2 に全症例，全手法の計測誤差をまとめる．

square error, RMSE) を用いた．

これらの結果より，計測誤差は x 方向は ROI 探

また，計算速度の評価指標として 1 秒当たり索，y 方向は DS が最小だったが，手動作成した正解の標準偏差を考えると，提案法の誤差も同の処理枚数を用いた．

等程度である．すなわち提案法は従来法と同程度以上の計測性能を達成できている．

データセット4.2

今回の検証に用いた臨床データは４症例である．用いた画像シーケンスは各症例 300 フレー
ムで構成され，画像の解像度は 0.259 mm/pixelである．今回実験に用いた症例の腫瘍はＸ線透視画像上で 80 × 80 pixel のブロックに収まる大きさであったので，DS で用いるブロックサイズを 80 × 80 pixel にした．

4.2.1

Fig. 4

症例 3 の計測誤差

Table 2

位置計測誤差 (mm)正解位置

臨床画像の場合，腫瘍位置の真値は不明である．そこで，実験評価のために医師らの指導を
症例 1
症例 2
症例 3
症例 4
平均

受けた 2 名が腫瘍位置時系列を手動で作成し，その平均を正解位置とした．作成した正解位置時系列の例を，各作成者により作成された位置

全探索
ROI 探索DS
提案法
RMSEx RMSEy RMSEx RMSEy RMSEx RMSEy RMSEx RMSEy
0.67
1.37
0.18
0.53
0.170.210.58
1.34
0.2
0.64
0.11
0.53
0.13
0.41
0.14
0.59
0.95
0.95
0.69
0.83
0.47
0.62
0.79
0.77
1.93
2.8
0.88
1.01
2.15
1.58
1.24
1.71
0.94
1.44
0.47
0.73
0.73
0.71
0.69
1.10とあわせて Fig. 3 に示す．また，Table 1 に各症例毎に作成者間の標準偏差をまとめる．

また，300 フレームの処理にかかった時間から 1 秒当たりに処理したフレーム数を算出した計算速度の結果を Table 3 にまとめる．実際の放射線治療に用いられる X 線透視装置のリアルタイム画像取得レートは最低 10 fps 程度である
ため，リアルタイム位置計測の最低条件として1 秒間に 10 フレーム以上の処理が求められる．
Fig. 3 症例 1 の正解位置 x 方向（上），y 方向（下）Table 3 より，その条件を達成しているのは提案法のみであり，十分高速な計測が可能であることがわかる．
これに対し，計測誤差が小さかった DS は計Table 1

正解作成者間の標準偏差 (mm)

標準偏差

症例 10.85

症例 20.33

症例 30.70

症例 40.62

平均0.62

算速度が不十分である一方，同様に誤差が小さかった ROI 探索は提案法を除いた各手法よりも高速であり，探索範囲の絞り込みが計測時間の短縮に有効であることが確認できる．今回は–5–

ROI 探索の先験情報として，事前に腫瘍位置の参考文献正解から移動振幅を求め，探索範囲を限定したが，臨床上は，事前に未来の腫瘍位置変動の振幅を正確に知ることはできないので，今回比較した ROI 探索は実際には用いることはできず，今回の実験環境ではリアルタイム計測の最低条件をクリアすることはできなかった．また仮に，事前情報として腫瘍の振幅の情報を得ることができたとしても計測中にその範囲を逸脱しないという保証はなく，実際そのような逸脱は臨床上よく見られる 8) ．

Table 3
症例 1
症例 2
症例 3
症例 4
平均

5.

全探索
0.076
0.074
0.074
0.086
0.078

計算速度 (frame/s)
ROI 探索
7.75
6.47
1.36
1.49
4.27

DS
0.73
0.79
0.50
0.82
0.71

提案法
20.3
19.5
12.6
21.3
18.41) 独 立 行 政 法 人 国 立 が ん 研 究 セ ン タ ー がん 対 策 情 報 セ ン タ ー：が ん の 統 計’12，
http://ganjoho.jp/data/public/
statistics/backnumber/2012/files/
cancer_statistics_2012.pdf
2) 上坂充, 中川恵一, 西尾禎治, 金井達明 (監修): 医学物理の理工学-上巻-，64/65，養賢堂 (2012)3) A. Sawant, R. L. Smith, R. B. Venkat, L.
Santanam, B. Cho, H. Cattell, L. J. Newell,
P. Parikh, and P. J. Keall: Toward submillimeter accuracy in the management ofintrafraction motion: The integration ofreal-time internal position monitoring andmultileaf collimator target tracking, Int.
J. Radiation Oncology Biol. Phys., 74-2,575/582 (2009)4) P. R. Poulsen, B. Cho, A. Sawant, D. Ruan,
and P. J. Keall: Detailed analysis of latencies in image-based dynamic MLC tracking, Med. Phys., 37-9, 4998/5005 (2010)

おわりに本研究では，過去の計測結果から現在の腫瘍位置を予測することにより，探索範囲を適応的に狭める高速な画像位置計測法を提案した．これにより X 線透視画像での肺腫瘍位置計測時間を，位置計測性能を犠牲にすることなく短縮することができ，リアルタイム計測の条件である毎秒 10 フレーム以上の処理性能を達成した．
今後より多くの症例において提案法の信頼性を詳細に検証をする必要があるが，今回実証された提案法の性能は非常に有望なものであり，動体追尾照射の実現に向けての貢献が期待される．

謝辞
本研究は，一部，栢森情報科学振興財団研究助成 K24 研 X VII 第 389 号により行われた．こ
こに記して謝意を表する．

–6–5) J. Meyer, A. Richter, K. Baier, J. Wilbert,
M. Guckenberger, M. Flentje: Trackingmoving objects with megavoltage portal
imaging: a feasibility study, Med. Phys.,
33-5, 1275/1280 (2006)
6) Shan Zhu and Kai-Kuang Ma: A new
diamond search algorithm for fast blockmatching motion estimation, IEEE Transactions on Image Processing, 9-2, 287/290(2000)
7) 北 川 源 四 郎:時 系 列 解 析 入 門,
(2005)

岩波出版

8) A. Richter, J. Wilbert, K. Baier, M. Flentje and M. Guckenberger: Feasibility Studyfor Markerless Tracking of Lung Tumorsin Stereotactic Body Radiotherapy, Int.
J. Radiation Oncology Biol. Phys., 78-2,618/627 (2010)