1 - 九州大学

医用画像処理学（５）
（教科書ｐｐ．123-)
有村秀孝
幾何学的変換
●画像の平行移動、回転、拡大縮小
アフィン変換（線形変換）を用いる
●画像の重ね合わせ（レジストレーション）
ＣＴ画像とＭＲ画像
ＭＲ画像とＰＥＴ画像
•主な画像処理解析技術
１．治療計画
２．臨床的標的体積(CTV)を決める目的で
治療計画CT と他のモダリティ画像（MRI,
PET など）との画像融合（フュージョン）
 レジストレーション
3
CT画像とMRI画像のフュージョンが必要な理由
治療計画CT と他のモダリティ画像（MRI, PET など）とのレジストレーションに
よる画像融合（フュージョン）が重要．
Khoo VS, Joon DL. New developments in MRI for target volume delineation in
radiotherapy. Br J Radiol. 2006 Sep;79 Spec No 1:S2-15.
4
レジストレーション
二つの画像間で対応する特徴点から求めた変換行列によって
画像間の位置合わせを行うこと
固定画像
移動画像
固定画像
特徴点はマニュアルで求める方法と自動的に
求める方法がある．
Tada M, et al. Digital Human Symposium 2009. Zitova B, Flusser J. Image registration methods:
5
a survey. Image and Vision Computing 2003; 21: 977-1000.
レジストレーション
二つの画像間(I1, I2)で対応する特徴点から求めた変換行列Rに
よって画像間の位置合わせを行うこと
●グローバル線形レジストレーション
アフィン変換
位相限定相関法
●局所的な非線形レジストレーション
Free Form Deformation (FFD)法
Thin plate spline (TPS)法
6
グローバル線形レジストレーション (アフィン変換）
C11 C12 C13 xn
Xn
Yn = C21 C22 C23 yn
0 0 1
1
1
X1
Y1
・
・
・
Xn
Yn
=
x1 y1 1
0 0 0
・
・
・
xn yn 1
0 0 0
変換後の座標
0 0 0
x1 y1 1
0 0 0
xn yn 1
変換前の座標
C11
C12
C13
C21
C22
C23
アフィン変換行
列の要素
アフィン変換行列を二つの画像間で対応する特徴点から求める
固定画像
（治療計画DRR画像）
特徴点
移動画像
（ポータル画像）
7
Global registration
Global Registrationでは、アフィン変換により、拡大縮小、
回転、平行移動、せん断の大まかな、位置合わせを行う。
入力：アフィン変換の要素を求める
ために用いる対応するfeature points の組（5組）
変換前の座標から構成されるベクトル A 、
変換後の座標から構成されるベクトル b 、
アフィン変換の要素から構成される行列 ｘ を作る
変換前のfeature pointsを特異値分解し、
A = U・W・VT の形に変換する
A・ｘ = ｂ を解くことによって、
アフィン変換行列の要素を求める
変換前の画像にアフィン変換を行う
その際、線形補間法により、補間する
出力：アフィン変換後の画像
アフィン変換行列 ｘ
C11 C12 C13 xn
Xn
Yn = C21 C22 C23 yn
0 0 1 1
1
変換後の座標
A
X1
Y1
・
・
・
Xn
Yn
=
x1 y 1 1
0 0 0
・
・
・
xn y n 1
0 0 0
変換前の座標
b
0 0 0
x1 y 1 1
0 0 0
xn y n 1
C11
C12
C13
C21
C22
C23
補足： Global Registration
変換前の座標点を特異値分解し、 A = U・W・VT の形に変換する
アフィン変換行列の要素を求める
未知数より、方程式が多いとき、未知数と方程式の数を揃えて、
連立方程式を解くために、特異値分解により、 A = U・W・VTの形に変換する。
[m][n]
[m][1]
[m][1]
[n][1]
A ・ ｘ
=
b
特異値分解
U
[n][n]
[n][n]
[n][1]
W
VT
・ ｘ =
[m][1]
b
A = U・W・VT
未知数ｘを求めるために、両辺に Aの逆行列（U・W・VT）-1 を左から掛ける。
（U・W・VT）-1・（U・W・VT ）・ x = （U・W・VT）-1・ b
x = （U・W・VT）-1・ b
∵直交行列 U, V の逆行列は転置に等しく、
対角行列 W の逆行列は、対角要素を逆数にするだけでよい。
x = V・W-1・UT・ b
つまり、
ｘ =
V
W-1
UT
b
より、
アフィン変換行列の要素（未知数 x）
を求める。
Global registration
1 Input images
2 Change of the matrix size
3 Decide the feature points by manual
4 Affine transform
Local Registration
Local Registrationでは、B-spline を用いた非線形変換を実行する。
概要：
CBF 画像上にcontrol points を定義し、 control points の x 軸方向、
y 軸方向、それぞれのTalairach atlas における座標を求め、その値を基にB-spline
により全ピクセルの座標を変換する。
CBF image の座標を(x,y)とし、 Talairach atlas の座標を(x’,y’) とし、
下式により、x,y方向ごとに、座標を求める。
( x' , y ' )  ( X ( x, y ), Y ( x, y ))
3
3
3
3
 ( Bk ( s ) Bl (t )   X ( i  k ), ( j l ) , Bk ( s ) Bl (t )  Y ( i  k ), ( j l ) )
k 0 l 0
k 0 l 0
この式において、未知数はφ である。これは、control points のもつ値
（ control points の x軸方向、y軸方向それぞれのTalairach atlas における座標）を
示しており、φ を求めなければならない。
グローバル線形レジストレーション (アフィン変換）
1
2
2
1
3
4
治療計画DRR画像
アフィン変換後画像
3
4
ポータル画像 (15°右回転)
重ね合わせ画像
12
（九州大学大学院生 板野航君作成）
FFD法（Free Form Deformation Method）
ー 局所的な非線形レジストレーションの例 ー
移動画像上に制御格子を設定し，格子上の制御点のピクセルと，固定画像
上の類似部位への対応付け（位置合わせ）を， B-spline関数を用いて非線
形に行う変換．
移動画像上の元座標を(x, y)とし、 固定画像上の移動先の座標を
(x’, y’) とする。
変形前
変形後
( x' , y ' )
( x, y )
移動画像の制御格子の変形
Tada M, et al. Digital Human Symposium 2009.
13
局所的な非線形レジストレーション
治療計画DRR画像
ポータル画像
(15°右回転、歪み)
非線型変換後画像
重ね合わせ画像
14
（九州大学大学院 山下泰生君）
脳解剖アトラス（ Talairach ）と T2強調画像 の
レジストレーション
１．Affine 変換による線形なレジストレーション
２．FFD法*を用いた局所的な非線形レジストレーション
脳解剖アトラス
T2強調画像
（九州大学大学院 山下泰生君）
フュージョン画像
非線形レジストレーション
による輪郭抽出支援
（１）正常構造と腫瘍の輪郭テンプレートを
患者のCT画像に大まかに重ねる．
（２）輪郭テンプレートと患者のCT画像の
間の非線形レジストレーションを行う．
（３）医師が自動抽出された輪郭を修正．
その効果は？
計画者間における
輪郭の変動低減
Chao KSC, et al, International journal of radiation
oncology, biology, physics, 68, 1512-1521, 2007.
16
アフィン変換（線形変換）
u  a10 x  a 01 y  a 00
v  b10 x  b01 y  b00
 u   a10 a 01 a 00  x 
  
 
 v    b10 b01 b00  y 
 1   0 0 1  1 
  
 
拡大
 2 0


 0 2
左回転
 cos 

 sin 
反転
1 0 


 0 1
 sin  

cos  
アフィン変換の演習
y
２倍拡大
x軸に関して反転
右回転
2
1
x
0
-1
1
2
方法
治療計画CT
（512*512*306）
PET/CTのCT
（512*512*306）
３Dアフィン変換
PET
（716*716*306）
切り取り 相互情報量を用いた
テンプレートマッチング
変換行列
３Dアフィン変換
PET
（512*512*306）
再補間結果
正しいピクセルサイズが判明
したので、再度cubic補間を
行った
予想値 5.5㎜
↓
実際には、5.46875㎜
512×512に切り取ったPET画像と
PETCTのCT画像を重ね合わせたもの
画像を幾何学的に変換するためのアルゴリズム
（１）変換後の座標系における格子点(x0,y0)に対応する変換前の座標(u0,v0)を
アフィン変換によって求める。しかし、求めた座標(u0,v0)は変換前の座標系で、
格子点とならない。 つまり、整数の座標とならず、非格子点となる。
（２）そこで、周囲の格子点のピクセル値を使って、非格子点の(u0,v0)のピクセル
値を求め、変換後の(x0,y0)におけるピクセル値とする。
(u0,v0)
(x0,y0)
非格子点におけるピクセル値の決定方法
（１）最近傍法(nearest neighbor)
f (u 0, v0)  f ([u 0  0.5],[v0  0.5])
（２）線形補間法(bi-linear interpolation)
f (u 0, v0)  f (u ' , v' )(1   )(1   )  f (u '1, v' ) (1   )
 f (u ' , v'1)(1   )   f (u '1, v'1)
（３）３次補間法(cubic convolution)
f (u 0, v0)   f (uk , vl )C (uk  u 0)C (vl  v0)
k
l
C(x)はsin(πx)/(πx)または３次の近似式（p127の式４．３２）
非格子点におけるピクセル値の決定方法
線形補間法の証明
f 1  (1   ) f (u' , v' )  f (u'1, v' )
f 2  (1   ) f (u' , v'1)  f (u'1, v'1)
f1
f (u 0, v0)  (1   ) f 1  f 2
(1-α)
y  (a  b)( x  x' )  b
 (1   )b  a
(1-β)
b
y
a
α
x’-α
f2
1-α
x’
x’+1-α
最近傍法
線形補間法
問
（１）任意の座標を４５度回転させて、次に２倍拡大する行列を求めなさい。ただ
し、回転行列は、以下の式を使うこと。
cos
-sin
sin
cos
（２）（１）で求めた行列を使って、頂点が(0,0), (1,0), (1,1), (0,1)の正方形を座標
変換した図形を図に描きなさい。
（３）角度によっては、回転した画像の座標は実数になるが、実数の座標のピ
クセル値を求める方法を３つ挙げ、１つだけ簡単に説明しなさい。