B-spline registration 〜非線形画像登録への入門〜

# はじめに
当記事では，PythonによるB-spline registrationの使用を目的とします．画像registrationとは，二つの画像間の適切な変形場を求める最適化計算です．
# Elastix
B-spline registrationを扱うためのソフトウェアとしてElastix[^1]が広く知られています．コードはC++で書かれていますが，wrapされていてPythonライブラリとして取得可能であり，今回はそちらを使用します．ライブラリは以下のコマンドでインストールされます．

```
pip install itk-elastix
```
詳細は以下のリンクをご参照ください．
https://github.com/InsightSoftwareConsortium/ITKElastix
# Data
githubに置かれている2次元CT脳画像を用います．3次元画像を扱う場合も同様の流れになります．

```
fixed_image = itk.imread('data/CT_2D_head_fixed.mha', itk.F)
moving_image = itk.imread('data/CT_2D_head_moving.mha', itk.F)
```
![2D_head_image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520106/e301433c-59e7-eb2a-ff02-36f513f2e0df.png)
# Registration
画像registrationでは，線形変換を行った後，非線形変換を行います．Elastixの線形変換は，

- Translation（平行移動）
- Rigid（回転と平行移動）
- Affine（Rigidと拡大縮小）

が基本的なものであり，非線形変換としては，

- B-spline
- Thin-plate spline 

が選択できます．Elastixでは，`parameter_object`を用意し，`parameter_map`を順次連結することで画像を変形することができます．

```
# Forward Parameter Map
parameter_object = itk.ParameterObject.New()
parameter_map_rigid = parameter_object.GetDefaultParameterMap('rigid')
parameter_object.AddParameterMap(parameter_map_rigid)
parameter_map_affine= parameter_object.GetDefaultParameterMap('affine')
parameter_object.AddParameterMap(parameter_map_affine)
parameter_map_bspline = parameter_object.GetDefaultParameterMap('bspline')
parameter_object.AddParameterMap(parameter_map_bspline)

# Registration
result_image, result_transform_parameters = itk.elastix_registration_method(
    fixed_image, moving_image,
    parameter_object=parameter_object,
    output_directory='exampleoutput/fwd'
)
```
登録結果は以下のようになります．
![2D_head_image_result.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520106/940bf83e-2c66-e9ff-059a-37108f21c727.png)
登録前後のLocal NCC[^2]を計算すると，各ピクセルにおいてNCCの絶対値が向上していることが確認できます．
![2D_head_image_ncc.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520106/8f3d85c3-dfb4-a5c9-1df6-19908088fe89.png)
# Point transform
画像変換と座標変換は逆方向であることに注意します．
moving画像からfixed画像へのregistrationとは，fixed空間のある座標をmoving空間のある座標に移動して，その座標のmoving画像の輝度値を補間することになります．これは，registrationによって疎な空間を作らないためです．

先のregistrationによって得られたmoving→fixed変形場を利用して，fixed空間の点群をmoving空間に移動します．例としてランダムな点を作成して，点群の移動を実行してみます．

```
# Create points
fixed_points = np.random.random([10,2])*200 + 30
np.savetxt("fixed_points.txt", fixed_points, fmt = "%.5f")

# Modify the file
with open("fixed_points.txt", 'r') as f:
    l = f.readlines()

l.insert(0, 'point\n')
l.insert(1, '10\n')

with open("fixed_points.txt", 'w') as f:
    f.writelines(l)
```
Elastixの点群移動関数`itk.transformix_filter`の入力ファイルは，冒頭に`point`ないし`index`と点群数を記入する必要があります．

```
result_image = itk.transformix_filter(
    moving_image, result_transform_parameters,
    fixed_point_set_file_name='fixed_points.txt',
    output_directory = './point_transformed')
```
出力ファイル`outputpoints.txt`の構造は，

```
Point	0	; InputIndex = [ 202 137 ]	; InputPoint = [ 201.934600 136.652650 ]	; OutputIndexFixed = [ 205 138 ]	; OutputPoint = [ 204.842360 138.311854 ]	; Deformation = [ 2.907760 1.659204 ]	; OutputIndexMoving = [ 205 138 ]
Point	1	; InputIndex = [ 212 157 ]	; InputPoint = [ 211.585340 156.863360 ]	; OutputIndexFixed = [ 211 158 ]	; OutputPoint = [ 210.705132 158.448602 ]	; Deformation = [ -0.880208 1.585242 ]	; OutputIndexMoving = [ 211 158 ]
```
のようになっているため，対象となる`OutputPoint`の部分のみを読み込みます．

```
# for 2D 
result_points = np.loadtxt('point_transformed/outputpoints.txt', dtype='str')[:,27:29].astype('float64')

# for 3D 
# result_points = np.loadtxt('point_transformed/outputpoints.txt', dtype='str')[:,30:33].astype('float64')
```
得られたfixed空間の点群`result_points`をplotすると以下のようになります．
![correspond_points.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520106/33d9346f-2926-3e71-6788-7861eed449df.png)
画像変形による点群の対応を確認することができました．
# Inverse transform
変形の`metric`を`DisplacementMagnitudePenalty`に設定することで，逆変換が求められます．
変換$T_\mu$，座標$x$に対して，`DisplacementMagnitudePenalty`とは，$T_\mu(x)-x$を表す変位に対するペナルティのことです．Registrationにおいては，入力画像と出力画像を同一に設定し，求めた順方向の変換を初期変換に与えます．


```
# Inverse Parameter Map
parameter_object_inv = itk.ParameterObject.New()
parameter_map_rigid = parameter_object_inv.GetDefaultParameterMap('rigid')
parameter_map_rigid['Metric'] = ["DisplacementMagnitudePenalty"]
parameter_object_inv.AddParameterMap(parameter_map_rigid)
parameter_map_affine= parameter_object_inv.GetDefaultParameterMap('affine')
parameter_map_affine['Metric'] = ["DisplacementMagnitudePenalty"]
parameter_object_inv.AddParameterMap(parameter_map_affine)
parameter_map_bspline = parameter_object_inv.GetDefaultParameterMap('bspline')
parameter_map_bspline['FinalGridSpacingInPhysicalUnits'] = np.array(['7']) 
parameter_map_bspline['Metric'] = ["DisplacementMagnitudePenalty"]
parameter_object_inv.AddParameterMap(parameter_map_bspline)

# Registration
result_inverse_image, result_inverse_parameter = itk.elastix_registration_method(
    moving_image,  moving_image,
    initial_transform_parameter_file_name = './exampleoutput/fwd/TransformParameters.2.txt',
    parameter_object=parameter_object_inv,
    output_directory='./exampleoutput/inv/',
    log_to_console=True
)
```
![restore_ncc.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520106/6a694505-cdca-f570-e689-2554767b8973.png)
順方向変換が与えられた時，逆方向変換を手に入れられることを確認できました．

# Troubleshooting
Registrationがうまくいかない場合として考えられる要因がいくつかあります．

- Registration中にエラーが生じる

```
Too many samples map outside moving image buffer
```
というエラーメッセージが度々出てきます．これは初期`Spacing`か位置の大幅なずれがあると生じやすいです．手動で`Spacing`を調整した相似変換ファイルを作成し適用したのち，translationを行うことで良好に作動します．

- Registrationの結果が歪である


いくつかのparameterが原因ですが，B-splineのグリッド幅を決定する`FinalGridSpacingInPhysicalUnits`と更新幅への制約である`MaximumStepLength`が特に重要です．3次元画像に適用するときには，一層気にかける必要があります．

 `FinalGridSpacingInPhysicalUnits`の比較
![FinalGridSpacingInPhysicalUnits_compare.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520106/820b1a6d-7008-e4c5-36ce-7f6e2e56b903.png)
 `MaximumStepLength`の比較
![MaximumStepLength_compare.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520106/60a2881d-1ac2-569d-a4fe-5273309b99e3.png)
# まとめ
以上がElastixを用いたB-spline registrationの一連の流れとなります．
registration中のエラーや，結果が不適切であることが繰り返したため，当記事を執筆しました．Hyperparameterの調整は慎重に行う必要があると思います．

[^1]: <a href = "https://doi.org/10.1109/TMI.2009.2035616">elastix: a toolbox for intensity-based medical image registration</a>

[^2]: <a href = "https://itk.org/Doxygen/html/classitk_1_1ANTSNeighborhoodCorrelationImageToImageMetricv4.html">itk::ANTSNeighborhoodCorrelationImageToImageMetricv4</a>