はじめに

当記事では，PythonによるB-spline registrationの使用を目的とします．画像registrationとは，二つの画像間の適切な変形場を求める最適化計算です．

Elastix

B-spline registrationを扱うためのソフトウェアとしてElastix¹が広く知られています．コードは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)

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'
)

登録結果は以下のようになります．

登録前後のLocal NCC²を計算すると，各ピクセルにおいてNCCの絶対値が向上していることが確認できます．

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すると以下のようになります．

画像変形による点群の対応を確認することができました．

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
)

順方向変換が与えられた時，逆方向変換を手に入れられることを確認できました．

Troubleshooting

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

Registration中にエラーが生じる

Too many samples map outside moving image buffer

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

Registrationの結果が歪である

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

FinalGridSpacingInPhysicalUnitsの比較

MaximumStepLengthの比較

まとめ

以上がElastixを用いたB-spline registrationの一連の流れとなります．
registration中のエラーや，結果が不適切であることが繰り返したため，当記事を執筆しました．Hyperparameterの調整は慎重に行う必要があると思います．

elastix: a toolbox for intensity-based medical image registration ↩
itk::ANTSNeighborhoodCorrelationImageToImageMetricv4 ↩