Pythonで学ぶ初歩的なITKの利用法

Notification

この記事は自分の勉強のメモとして書きました。
参照元はこのノートブックです。
https://github.com/KitwareMedical/2019-03-13-KRSCourseInBiomedicalImageAnalysisAndVisualization

ITK

画像を可視化する技術に特化したバイオメディカル分野の画像処理に欠かすことのできないライブラリ。
SlicerやImageJで利用されている。

環境

Google Colab
（任意のノートブックを立ち上げておく。）

サンプルデータ

先にサンプルデータをGitHubからsubversionを使ってダウンロード。

!sudo apt install subversion
#データを取得しておきます
#GithubのURLにk含まれる/tree/masterを/trunkへ置き換えて使います
!svn checkout https://github.com/KitwareMedical/2019-03-13-KRSCourseInBiomedicalImageAnalysisAndVisualization/trunk/data

ITKのインストール

#ITKをインストールする
!pip install itk

Image Filtering

ここで学ぶこと

• ITKで使用される画像処理パイプラインモデルを知る。
• ImageやMeshはデータオブジェクトを表すために利用される。
• データオブジェクトを操作するプロセス内で新しいデータオブジェクトが生成される。
• プロセスは、プロセスオブジェクトとして分類され、source、filter、mapperオブジェクトがある。
• source（readerなど）はデータを取得して、オブジェクトをフィルター処理して新しいデータを生成する。
• mapperは、ファイルまたは他のシステムへの出力用にデータを受け入れる。
  

フィルターという用語は、次の3つのタイプとして広く使用される。

• Data pipeline（データパイプライン） 通常、データオブジェクトとプロセスオブジェクトは、SetInput()およびGetOutput()メソッドを使用して相互に変換されます。 何かの処理に対する新しい出力や新しいピクセルデータは、パイプラインの最後（プロセスオブジェクトまたはデータオブジェクト）でUpdate()メソッドが呼び出されるまで発生しません。
• Pipeline updates（パイプラインの更新） 多次元画像データは大きく、さらに大きくすることもできます。 ITKはデータストリーミング機能を介してこの問題に対処します。
• Streaming（ストリーミング） ストリーミングは、イメージをパイプラインの最後で重複しない領域に分割することにより実行されます。その後、RequestedRegionとしてパイプラインを伝播します。

ITKには、3つのImageRegionの種類がある。

• BufferedRegion:メモリに保存されたピクセルの領域
• LargestPossibleRegion: 画像の可能な最大領域
• RequestedRegion: ストリーミング時に単一の処理パスで要求された領域（BufferedRegionとLargestPossibleRegionは、RequestedRegionと同じかそれよりも大きい。）

やってみる。

import numpy as np
import itk
import matplotlib.pyplot as plt

#パックマンで空間フィルタリングを練習
fileName = 'data/PacMan.png'
reader = itk.ImageFileReader.New(FileName=fileName)
#ガウシアンフィルターを使ってみる
smoother = itk.RecursiveGaussianImageFilter.New(Input=reader.GetOutput())

print("reader's Output: %s" % reader.GetOutput())
print("smoother's Output: %s" % smoother.GetOutput())

#処理結果を生成するためにUpdate()を忘れない。
smoother.Update()

print("reader's Output: %s" % reader.GetOutput())
print("smoother's Output: %s" % smoother.GetOutput())

image = reader.GetOutput()#オリジナルの確認
plt.subplot(1,2,1),plt.title("original"),plt.imshow(itk.GetArrayFromImage(image))
smoothed = smoother.GetOutput()#スムージング後
plt.subplot(1,2,2),plt.title("smoothing"),plt.imshow(itk.GetArrayFromImage(smoothed))

（すでに行っているが）パイプラインオブジェクトでアップデートしていれば、処理結果は最新になる

smoother.Update()
# ただし、処理を変更した場合は再度アップデートする必要がある。
smoother.SetSigma(10.0)
smoother.Update()
plt.imshow(itk.GetArrayFromImage(smoothed))

また、あまりないと思うが、readerに途中で変更を加えた場合は、readerを更新することを忘れない。

# reader.Modified()

# いろいろなフィルターを試すことができる。
# https://itk.org/Doxygen/html/group__IntensityImageFilters.html
fileName = 'data/PacMan.png'
reader = itk.ImageFileReader.New(FileName=fileName)
smoother = itk.MedianImageFilter.New(Input=reader.GetOutput())#例えばメディアン
smoother.SetRadius(5)
smoother.Update()
plt.imshow(itk.GetArrayFromImage(smoother.GetOutput()))

補足

# パイプラインの最後にStreamingImageFilterを配置して、パイプラインをストリーミングする。
# smootherは、画像領域のストリーミング分割ごとに、出力を複数回生成する。
# readerはストリーミングできないため、出力は1回しか生成されない。
streamer = itk.StreamingImageFilter.New(Input=smoother.GetOutput())
streamer.SetNumberOfStreamDivisions(4)
reader.Modified()
streamer.Update()
# 書き込みは拡張子を含めて指定するだけでよい
writer = itk.ImageFileWriter.New(Input=smoother.GetOutput())
writer.SetFileName('my_output.mha')
writer.SetNumberOfStreamDivisions(4)
writer.Update()

Image Segmentation

ここで学ぶこと

• ITKで利用可能なセグメンテーション法を知る（まずはConfidenceConnectedImageFilter）
• リージョングローイングとそのパラメータがどのように動作を変えるかを理解する
• レベルセットとそのパラメータがどのように動作を変えるかを理解する

ITKで利用できるセグメンテーションの種類
- Level Set
- Watershed（画像特徴の勾配降下と領域境界に沿った分析）
- Connected component（二値化やラベリングで利用される接続されたオブジェクトを扱うための抽象的な技術）
- Label voting（ラベルとしたピクセル値をカウントしてプロセスの中で利用するためのもの）
- Region Growing
- Classifier（機械学習によるアルゴリズム）

ここでは、RegionGrowingとLevelSetを対象に処理例を紹介する。

Region growing

領域膨張アルゴリズムは、セグメント化されるオブジェクトの内部にあると見なされるシード領域（通常は1つ以上のピクセル）から膨張を開始するのが基本的なアプローチ。
このセグメント化されるオブジェクト領域に隣接するピクセルは、オブジェクトの一部と見なされるべきかどうかを判断するために評価される。
セグメント化される場合、それらは領域に追加され、新しいピクセルが領域に追加される限り、プロセスが続行される。
領域膨張アルゴリズムは、領域にピクセルを含めるかどうかを決定するために使用される基準、近隣を決定するために使用される接続のタイプ、および近隣のピクセルを探索するために使用される処理（またはストラテジー）によって異なる。

ここでは例として次の方法を用いる。
- ConfidenceConnected
- ConnectedThreshold
- IsolatedConnected

まずは、「ConfidenceConnected」の使い方の確認から。
ConfidenceConnectedは、領域境界の上限値と下限値を決めて、この間にある値から標準偏差を求め、この標準偏差の許容範囲（レバレッジのようなもの。SD * Multiplier）を設定する方法。

# BrainProtonDensitySlice.pngを利用する
# 'unsigned char'でロードするのはデータ変換の方法を後で示すため。
input_image = itk.imread('data/BrainProtonDensitySlice.png', itk.ctype('unsigned char'))
# View the input image
plt.imshow(itk.GetArrayFromImage(input_image))

「ConfidenceConnectedImageFilter」クラスの使い方の確認。
使い方が分からないときは調べる。

confidence_connected = itk.ConfidenceConnectedImageFilter.New(input_image)
#pythonのヘルプ機能を使って使い方を確認する
help(confidence_connected)

# confidence connected filter parameterを任意に設定する
confidence_connected.SetMultiplier(2.3)#許容範囲、この例では±SD*2.3
confidence_connected.SetNumberOfIterations(5)
confidence_connected.SetInitialNeighborhoodRadius(3)
confidence_connected.SetReplaceValue(255)

関数の使い方が分からないときは調べる。

# *Multiplier*?
confidence_connected.SetMultiplier?

領域膨張の起点「シードポイント」を設定。ITKのImageオブジェクトの原点は画像座標の左下であることに注意。（NumPyは左上）

# Set the seed points
confidence_connected.SetSeed([100, 110])
#実行して表示
confidence_connected.Update()
plt.imshow(itk.GetArrayFromImage(confidence_connected))

# 処理結果の改善を試みる。原画像にメディアンフィルタをかけてみる。
median_filtered_image = itk.median_image_filter(input_image, radius=1)
# 再度、パラメータを更新。メディアンフィルタ後の画像を入力。
confidence_connected_image = itk.confidence_connected_image_filter(median_filtered_image,
                                                                   seed=[100,110],
                                                                   replace_value=255,
                                                                   multiplier=3.0,
                                                                   number_of_iterations=5,
                                                                   initial_neighborhood_radius=3)
plt.imshow(itk.GetArrayFromImage(confidence_connected_image))

次は、3つのリージョングローイング処理を比較してみる。
ConnectedThresholdImageFilter vs ConfidenceConnectedmageFilter vs IsolatedConnectedImageFilter

# ConnectedThresholdImageFilter vs ConfidenceConnectedmageFilter vs IsolatedConnectedImageFilter
input_image = itk.imread('data/BrainProtonDensitySlice.png', itk.ctype('unsigned char'))
connected_threshold = itk.ConnectedThresholdImageFilter.New(input_image)
connected_threshold.SetLower(190)
connected_threshold.SetUpper(255)
connected_threshold.SetReplaceValue( 255 )
connected_threshold.SetSeed( [100, 110] );
connected_threshold.Update()

confidence_connected = itk.ConfidenceConnectedImageFilter.New(input_image)
confidence_connected.SetMultiplier(2.3)
confidence_connected.SetSeed( [100, 110] );
confidence_connected.AddSeed( [80, 126] );
confidence_connected.SetNumberOfIterations(5)
confidence_connected.SetInitialNeighborhoodRadius(3)
confidence_connected.SetReplaceValue(255)
confidence_connected.Update()

isolated_connected = itk.IsolatedConnectedImageFilter.New(input_image)
isolated_connected.SetSeed1([98, 112])
isolated_connected.SetSeed2([98, 136])
isolated_connected.SetUpperValueLimit( 245 )
isolated_connected.FindUpperThresholdOff();
isolated_connected.Update()

plt.subplot(131)
plt.title("ConnectedThreshold")
plt.imshow(itk.GetArrayFromImage(connected_threshold))
plt.subplot(132)
plt.title("ConfidenceConnected")
plt.imshow(itk.GetArrayFromImage(confidence_connected))
plt.subplot(133)
plt.title("IsolatedConnected")
plt.imshow(itk.GetArrayFromImage(isolated_connected))
plt.tight_layout()
plt.show()

Fast Marching Segmentation(LevelSet)

次に、LevelSet法の高速マーチング（Fast Marching Segmentation）を試す。
セグメンテーションの対象領域が単純な構造なら、高速マーチングと呼ばれる高速膨張アルゴリズムを使用できる。
空間フィルタを使って画像を単純化してから試してみる。

# smoothing filterが対応できる画像のデータタイプを調べる
itk.CurvatureAnisotropicDiffusionImageFilter.GetTypes()

# サポートされているデータタイプに画像データをキャストする方法
Dimension = input_image.GetImageDimension()
FloatImageType = itk.Image[itk.ctype('float'), Dimension]
caster = itk.CastImageFilter[input_image, FloatImageType].New()
caster.SetInput(input_image)

#スムージングフィルタをかける。
smoothed_image = itk.curvature_anisotropic_diffusion_image_filter(caster, 
                                                                  time_step=0.125, 
                                                                  number_of_iterations=5,
                                                                  conductance_parameter=3.0)

plt.imshow(itk.GetArrayFromImage(smoothed_image))

# そのままエッジ検出フィルタを実行
gradient = itk.gradient_magnitude_recursive_gaussian_image_filter(smoothed_image, sigma=1.0)
plt.imshow(itk.GetArrayFromImage(gradient))

# 続けてsigmoidフィルターを実行

basin_minimum = 2.25
border_minimum = 3.75

alpha = - (border_minimum - basin_minimum) / 3.0
beta = (border_minimum + basin_minimum) / 2.0

sigmoid = itk.sigmoid_image_filter(gradient, output_minimum=0.0, output_maximum=1.0, alpha=alpha, beta=beta)
plt.imshow(itk.GetArrayFromImage(sigmoid))

#高速マーチング
fast_marching = itk.FastMarchingImageFilter.New(sigmoid)

seed_value = 0.0
NodeType = itk.LevelSetNode[itk.ctype('float'), Dimension]
node = NodeType()
node.SetIndex([100,110])
node.SetValue(seed_value)

NodeContainerType = itk.VectorContainer[itk.ctype('unsigned int'), NodeType]
nodes = NodeContainerType.New()
nodes.Initialize()
nodes.InsertElement(0, node)

fast_marching_image = itk.fast_marching_image_filter(sigmoid, trial_points=nodes, stopping_value=80)
plt.imshow(itk.GetArrayFromImage(fast_marching_image))

画像を二値化する。

time_threshold = 100 #閾値を設定
thresholded_image = itk.binary_threshold_image_filter(fast_marching_image,
                                                      lower_threshold = 0,
                                                      upper_threshold=time_threshold,
                                                      outside_value=0,
                                                      inside_value=255)
plt.imshow(itk.GetArrayFromImage(thresholded_image))

ここからShape Detection Level Setを試していく。処理の実行速度は向上するが、曲率はよくないらしい。
パイプラインは次のようになるようだ。

先に図中下段のBalanced画像を作る。

# バイナリ画像はWhiteMatter.pngを使う
binary_mask = itk.imread('data/WhiteMatter.png', itk.ctype('float'))

smoothed_image = itk.smoothing_recursive_gaussian_image_filter(binary_mask)

# 階調処理
rescaler = itk.IntensityWindowingImageFilter[FloatImageType,FloatImageType].New(smoothed_image)
rescaler.SetWindowMinimum(0)
rescaler.SetWindowMaximum(255)
# 入力レベルを内部で負の値を持つように反転
rescaler.SetOutputMinimum(5.0)
rescaler.SetOutputMaximum(-5.0)

rescaler.Update()
plt.imshow(itk.GetArrayFromImage(rescaler))

次に、図中のfeature imageを作っていく。

#上記の手順でロード済みの画像を使って前処理をしていく
gradient = itk.gradient_magnitude_recursive_gaussian_image_filter(caster, sigma=1.0)

basin_minimum = 1
border_minimum = 5.0

alpha = - (border_minimum - basin_minimum) / 3.0
beta = (border_minimum + basin_minimum) / 2.0

# 図中のfeature image
sigmoid = itk.sigmoid_image_filter(gradient, output_minimum=0.0, output_maximum=1.0, alpha=alpha, beta=beta)
plt.imshow(itk.GetArrayFromImage(sigmoid))

#shape_detection_level_set処理を実行
shape_detection_image = itk.shape_detection_level_set_image_filter(rescaler,
                                                                   feature_image=sigmoid,
                                                                   maximum_r_m_s_error=0.001,
                                                                   number_of_iterations=100,
                                                                   propagation_scaling=1.0,
                                                                   curvature_scaling=2.0)
#二値化
thresholded_image = itk.binary_threshold_image_filter(shape_detection_image,
                                                      lower_threshold=-1e7,
                                                      # Cut at the zero set
                                                      upper_threshold=0.0,
                                                      outside_value=0,
                                                      inside_value=255)
plt.imshow(itk.GetArrayFromImage(thresholded_image))

シェイプ検出レベルセットにレベルセットの曲率(Curvature)と伝播(Propagation)の条件に関わるいくつかのパラメータを導入してみる
Min Basin
Min Boundary
Propagation Scaling
Curvature Scaling
Number of Iterations

def explore_shape_detection_image_parameters(basin_minimum, boundary_minimum, propagation_scaling,
                                             curvature_scaling, number_of_iterations):

    input_image = itk.imread('data/BrainProtonDensitySlice.png', itk.ctype('unsigned char'))

    # Prepare the initial level set image
    binary_mask = itk.imread('data/WhiteMatter.png', itk.ctype('float'))

    smoothed_image = itk.smoothing_recursive_gaussian_image_filter(binary_mask)

    Dimension = input_image.GetImageDimension()
    FloatImageType = itk.Image[itk.ctype('float'), Dimension]
    caster = itk.CastImageFilter[input_image, FloatImageType].New()
    caster.SetInput(input_image)

    # Prepare the speed image
    gradient = itk.gradient_magnitude_recursive_gaussian_image_filter(caster, sigma=1.0)

    alpha = - (boundary_minimum - basin_minimum) / 3.0
    beta = (boundary_minimum + basin_minimum) / 2.0

    sigmoid = itk.sigmoid_image_filter(gradient, output_minimum=0.0, output_maximum=1.0, alpha=alpha, beta=beta)

    rescaler = itk.IntensityWindowingImageFilter[FloatImageType,FloatImageType].New(smoothed_image)
    rescaler.SetWindowMinimum(0)
    rescaler.SetWindowMaximum(255)
    # Invert the input level set to have negative values inside
    rescaler.SetOutputMinimum(5.0)
    rescaler.SetOutputMaximum(-5.0)

    rescaler.Update()

    shape_detection_image = itk.shape_detection_level_set_image_filter(rescaler,
                                                                       feature_image=sigmoid,
                                                                       maximum_r_m_s_error=0.001,
                                                                       number_of_iterations=number_of_iterations,
                                                                       propagation_scaling=propagation_scaling,
                                                                       curvature_scaling=curvature_scaling)

    thresholded_image = itk.binary_threshold_image_filter(shape_detection_image,
                                                          lower_threshold=-1e7,
                                                          # Cut at the zero set
                                                          upper_threshold=0.0,
                                                          outside_value=0,
                                                          inside_value=255)
    plt.imshow(itk.GetArrayFromImage(thresholded_image))
    plt.title(str(basin_minimum)+"-"+str(boundary_minimum)+"-"+str(propagation_scaling)+"-"+str(curvature_scaling)+"-"+str(number_of_iterations))
    plt.show()
    return

explore_shape_detection_image_parameters(basin_minimum=1, boundary_minimum=5., 
                                         propagation_scaling=1., curvature_scaling=2., number_of_iterations=100)
explore_shape_detection_image_parameters(basin_minimum=1, boundary_minimum=5., 
                                         propagation_scaling=1., curvature_scaling=2., number_of_iterations=500)
explore_shape_detection_image_parameters(basin_minimum=1, boundary_minimum=5., 
                                         propagation_scaling=1., curvature_scaling=4., number_of_iterations=100)
explore_shape_detection_image_parameters(basin_minimum=1, boundary_minimum=20., 
                                         propagation_scaling=1., curvature_scaling=1., number_of_iterations=100)

Image Registration

ここで学ぶこと

• ITKが物理空間で計算を行う方法を知る
• ITKがピクセル空間ではなく物理空間で位置合わせ（レジストレーション）を行う理由を知る
• ITKレジストレーションフレームワークのコンポーネントに触れ、利用可能な値を知る

# 位置合わせの前に、対象のデータ同士が同じ幾何空間にキャリブレーションされているか。リサンプリングが必要ならやる。
# IJのサンプルでやるかな。 future work
# https://imagej.nih.gov/ij/images/
# https://imagej.nih.gov/ij/images/pet-series/
# https://imagej.nih.gov/ij/images/ct-scan.zip

#利用するデータをロード（基準画像）
PixelType = itk.ctype('float')
fixedImage = itk.imread('data/BrainProtonDensitySliceBorder20.png', PixelType)
plt.imshow(itk.array_view_from_image(fixedImage))

#位置がずれている画像
movingImage = itk.imread('data/BrainProtonDensitySliceShifted13x17y.png', PixelType)
plt.imshow(itk.array_view_from_image(movingImage))

# 変換に必要な情報を集める
# 次元数
Dimension = fixedImage.GetImageDimension()
# 画像のタイプ
FixedImageType = type(fixedImage)
MovingImageType = type(movingImage)
#変換するデータのタイプ
TransformType = itk.TranslationTransform[itk.D, Dimension]
#受け皿の初期化
initialTransform = TransformType.New()
#位置合わせ誤差を最小にするように最適化する関数を定義
optimizer = itk.RegularStepGradientDescentOptimizerv4.New(
        LearningRate=4,
        MinimumStepLength=0.001,
        RelaxationFactor=0.5,
        NumberOfIterations=200)
#位置合わせのアルゴリズムを設定（平均自乗誤差）
metric = itk.MeanSquaresImageToImageMetricv4[
    FixedImageType, MovingImageType].New()
#位置合わせオブジェクトを初期化
registration = itk.ImageRegistrationMethodv4.New(FixedImage=fixedImage,
        MovingImage=movingImage,
        Metric=metric,
        Optimizer=optimizer,
        InitialTransform=initialTransform)
# 位置合わせの設定
# 移動する画像の変換タイプ
movingInitialTransform = TransformType.New()
initialParameters = movingInitialTransform.GetParameters()
initialParameters[0] = 0 #x
initialParameters[1] = 0 #y
movingInitialTransform.SetParameters(initialParameters)
registration.SetMovingInitialTransform(movingInitialTransform)
# 基準の変換タイプ
identityTransform = TransformType.New()
identityTransform.SetIdentity()
registration.SetFixedInitialTransform(identityTransform)

registration.SetNumberOfLevels(1) # number of multi-resolution levels
registration.SetSmoothingSigmasPerLevel([0]) # スムージングレベル
registration.SetShrinkFactorsPerLevel([1]) # シュリンク（収縮）レベル
#位置合わせ設定を更新して確定
registration.Update()
# 変換処理に関する結果を取得
transform = registration.GetTransform()
finalParameters = transform.GetParameters()
# 移動距離
translationAlongX = finalParameters.GetElement(0)
translationAlongY = finalParameters.GetElement(1)
# 繰り返し回数
numberOfIterations = optimizer.GetCurrentIteration()
# 手法に関する精度（小さいほどよい）
bestValue = optimizer.GetValue()

print("Result = ")
print(" Translation X = " + str(translationAlongX))
print(" Translation Y = " + str(translationAlongY))
print(" Iterations    = " + str(numberOfIterations))
print(" Metric value  = " + str(bestValue))

# 画像の合成
CompositeTransformType = itk.CompositeTransform[itk.D, Dimension]
outputCompositeTransform = CompositeTransformType.New()
outputCompositeTransform.AddTransform(movingInitialTransform)
outputCompositeTransform.AddTransform(registration.GetModifiableTransform())
# リサンプリング
resampler = itk.ResampleImageFilter.New(Input=movingImage,
        Transform=outputCompositeTransform,
        UseReferenceImage=True,
        ReferenceImage=fixedImage)
#結果表示
resampler.SetDefaultPixelValue(100)

OutputPixelType = itk.ctype('unsigned char')
OutputImageType = itk.Image[OutputPixelType, Dimension]

resampler.Update()

result = resampler.GetOutput()
plt.subplot(1,2,1)
plt.title("registered")
plt.imshow(itk.array_view_from_image(result))
plt.subplot(1,2,2)
plt.title("moving_org")
plt.imshow(itk.array_view_from_image(movingImage))

位置ずれを差分で見てみる

difference = itk.SubtractImageFilter.New(Input1=fixedImage,
        Input2=resampler)
resampler.SetDefaultPixelValue(1)

difference.Update()
plt.title("difference")
plt.imshow(itk.array_view_from_image(difference.GetOutput()))

# レジストレーション前後のずれ
resampler.SetTransform(identityTransform)
difference.Update()
plt.imshow(itk.array_view_from_image(difference.GetOutput()))

次に、異なる種類の画像（マルチモーダル）を位置合わせする。
マルチモーダルな場合の位置合わせによく利用されるMattesMutualInformation（マテス相互情報量）を使う。
剛体変換を使う。

# 基準画像
t1_image = itk.imread('data/DzZ_T1.nrrd', itk.F)
print("Pixel spacing:", t1_image.GetSpacing())
plt.imshow(itk.array_view_from_image(t1_image))

# 回転しているT2強調画像
t2_image = itk.imread('data/DzZ_T2.nrrd', itk.F)
print("Pixel spacing:", t2_image.GetSpacing())
plt.imshow(itk.array_view_from_image(t2_image))

#補足：方向を取得する
pyDirM=itk.GetArrayFromVnlMatrix(t2_image.GetDirection().GetVnlMatrix().as_matrix())
print("Direction matrix:\n", pyDirM)

マルチモーダルな場合の位置合わせによく利用されるMattesMutualInformation（マテス相互情報量）による相関を使って位置合わせする。

T1ImageType = type(t1_image)
T2ImageType = type(t2_image)
assert T1ImageType==T2ImageType, "Images must have the same pixel type!"

#変換タイプを剛体変換として初期化
TransformTypeR = itk.Rigid2DTransform[itk.D]
initialTransform = TransformTypeR.New()

MetricType = itk.MattesMutualInformationImageToImageMetricv4[
        T1ImageType, T2ImageType]
metric = MetricType.New()

scales = itk.OptimizerParameters[itk.D](initialTransform.GetNumberOfParameters())
translation_scale = 1.0 / 1000
scales.SetElement(0, 1.0)
scales.SetElement(1, translation_scale)
scales.SetElement(2, translation_scale)

#スケールを自動推定することもできるが、MattesMutualInformationでは使えないようだ。
#ScalesEstimatorType = itk.RegistrationParameterScalesFromPhysicalShift[MetricType]
#scalesEstimator = ScalesEstimatorType.New(Metric=metric)

optimizer = itk.RegularStepGradientDescentOptimizerv4.New(
        Scales=scales,
        #ScalesEstimator=scalesEstimator,
        #LearningRate=1.0,
        MinimumStepLength=0.001,
        # RelaxationFactor=0.8,
        NumberOfIterations=300)

registration = itk.ImageRegistrationMethodv4.New(FixedImage=t1_image,
        MovingImage=t2_image,
        Metric=metric,
        MetricSamplingPercentage=0.1, # 10%
        Optimizer=optimizer,
        InitialTransform=initialTransform)

initialParameters = initialTransform.GetParameters()
# オフセットの設定
initialParameters[0] = 0.0 # rotation in radians
initialParameters[1] = -40.0 # x translation in millimeters
initialParameters[2] = +30.0 # y translation in millimeters

initialTransform.SetParameters(initialParameters)

resampler = itk.ResampleImageFilter.New(Input=t2_image,
        Transform=initialTransform,
        UseReferenceImage=True,
        ReferenceImage=t1_image)
resampler.SetDefaultPixelValue(0)
resampler.Update()

plt.subplot(131)
plt.imshow(itk.GetArrayFromImage(resampler.GetOutput()))
plt.subplot(132)
plt.imshow(itk.GetArrayFromImage(t1_image))

結構ずれる。パラメータを調整すれば良いのかも。
めげずに、位置合わせ後の画像を使ってもう一回繰り返してみる。

registration.SetMovingInitialTransform(initialTransform)

Dimension = t1_image.GetImageDimension()
identityTransform = TransformTypeR.New()
identityTransform.SetIdentity()
registration.SetFixedInitialTransform(identityTransform)

registration.SetNumberOfLevels(1)
registration.SetSmoothingSigmasPerLevel([0])
registration.SetShrinkFactorsPerLevel([1])

try:
    registration.Update()
except RuntimeError as exc:
    print("Exception ocurred:\n", exc, "\n\n")
finally:
    print("Registration finished")

CompositeTransformType = itk.CompositeTransform[itk.D, Dimension]
outputCompositeTransform = CompositeTransformType.New()
outputCompositeTransform.AddTransform(registration.GetTransform())
outputCompositeTransform.AddTransform(initialTransform)

resamplerResult = itk.ResampleImageFilter.New(Input=t2_image,
        Transform=outputCompositeTransform,
        UseReferenceImage=True,
        ReferenceImage=t1_image)
resamplerResult.SetDefaultPixelValue(0)
resamplerResult.Update()
# final position of T2 image when resampled into T1 image's pixel grid
plt.imshow(itk.GetArrayFromImage(resamplerResult.GetOutput()))

位置合わせの情報を出力する。

transform = registration.GetTransform()
finalParameters = transform.GetParameters()
rotationInRadians = finalParameters.GetElement(0)
translationAlongX = finalParameters.GetElement(1)
translationAlongY = finalParameters.GetElement(2)
rotationInDegrees = rotationInRadians * 360 / 3.141592 # radian to degree
numberOfIterations = optimizer.GetCurrentIteration()

bestValue = optimizer.GetValue()

print("Result = ")
print(" Rotation degr = " + str(rotationInDegrees))
print(" Translation X = " + str(translationAlongX))
print(" Translation Y = " + str(translationAlongY))
print(" Rotation rad. = " + str(rotationInRadians))
print(" Iterations    = " + str(numberOfIterations))
print(" Metric value  = " + str(bestValue))

差分を見てみる。

difference = itk.AddImageFilter.New(Input1=t1_image,
        Input2=resamplerResult)
difference.Update()
plt.imshow(itk.GetArrayFromImage(difference.GetOutput()))

以上です。

Reference

• https://github.com/KitwareMedical/2019-03-13-KRSCourseInBiomedicalImageAnalysisAndVisualization
• https://itk.org/ITKSoftwareGuide/html/Book2/ITKSoftwareGuide-Book2ch4.html#x37-1770004