Python ITK 線性仿射影像對準 Affine Registration 教學與範例

介紹如何在 Python 中使用 ITK 進行影像的線性仿射影像對準（affine registration）。

準備測試用影像

這裡我們使用一張影像作為基準影像（fixed image），以 AffineTransform 線性仿射轉換套用至基準影像上，產生一張調動影像（moving image）。

import itk
import matplotlib.pyplot as plt

# 影像像素資料類型
PixelType = itk.SS

# 實體座標系統用的資料類型
ScalarType = itk.D

# 影像維度
Dimension = 2

# 影像類型
ImageType = itk.Image[PixelType, Dimension]

# 建立影像 Reader
ReaderType = itk.ImageFileReader[ImageType]
reader = ReaderType.New()
reader.SetFileName("fixed.mhd")
reader.Update()

# 讀取 Fixed 影像
fixedImage = reader.GetOutput()

# 建立 AffineTransform 線性轉換
affineTransform = itk.AffineTransform[ScalarType, Dimension].New()
parameters = affineTransform.GetParameters()
parameters[0] = 0.9 # 旋轉與變形
parameters[1] = 0.1 # 旋轉與變形
parameters[2] = 0.1 # 旋轉與變形
parameters[3] = 0.9 # 旋轉與變形
parameters[4] = 25  # X 軸平移
parameters[5] = -5  # Y 軸平移
affineTransform.SetParameters(parameters)

# 建立 Resampler
resamplerType = itk.ResampleImageFilter[ImageType, ImageType]
resampleFilter = resamplerType.New()

# 設定 Resampler 各項參數
resampleFilter.SetInput(fixedImage)
resampleFilter.SetTransform(affineTransform)
resampleFilter.SetReferenceImage(fixedImage)
resampleFilter.SetUseReferenceImage(True)
resampleFilter.Update()

# 產生 Moving 影像
movingImage = resampleFilter.GetOutput()

# 分開顯示兩張影像
fig, axs = plt.subplots(1, 2)
axs[0].imshow(itk.GetArrayViewFromImage(fixedImage), cmap='gray')
axs[0].set_title('Fixed Image')
axs[1].imshow(itk.GetArrayViewFromImage(movingImage), cmap='gray')
axs[1].set_title('Moving Image')
plt.show()

仿射對準（Affine Registration）

由於為了縮減 Python 模組的大小，Python 中只支援特定幾種影像類型的對準，所以要進行對準之前，必須先把影像轉換為支援的類型，再進行後續的對準動作。

# 影像轉型
InternalImageType = itk.Image[itk.F, Dimension]
fixedCastImageFilter = itk.CastImageFilter[ImageType, InternalImageType].New()
fixedCastImageFilter.SetInput(fixedImage)
movingCastImageFilter = itk.CastImageFilter[ImageType, InternalImageType].New()
movingCastImageFilter.SetInput(movingImage)

# 影像對準所使用的初始 Transform
TransformType = itk.AffineTransform[ScalarType, Dimension]
initialTransform = TransformType.New()
initialTransform.SetIdentity()

# 參數數量
numOfParam = initialTransform.GetNumberOfParameters()
print("Number of parameters: {}".format(numOfParam))

Number of parameters: 6

# 影像對準所使用的 Optimizer
optimizer = itk.RegularStepGradientDescentOptimizerv4.New()
optimizer.SetLearningRate(1)
optimizer.SetMinimumStepLength(0.0001)
optimizer.SetNumberOfIterations(200)
optimizer.SetRelaxationFactor(0.9)

# Optimizer Scales
optimizerScales = itk.OptimizerParameters.D(numOfParam)
optimizerScales.SetElement(0, 1.0)  # 旋轉與變形
optimizerScales.SetElement(1, 1.0)  # 旋轉與變形
optimizerScales.SetElement(2, 1.0)  # 旋轉與變形
optimizerScales.SetElement(3, 1.0)  # 旋轉與變形
optimizerScales.SetElement(4, 0.05) # X 軸平移
optimizerScales.SetElement(5, 0.05) # X 軸平移
optimizer.SetScales(optimizerScales)

# 影像對準所使用的 Metric
#MetricType = itk.MeanSquaresImageToImageMetricv4[InternalImageType, InternalImageType]
MetricType = itk.MattesMutualInformationImageToImageMetricv4[InternalImageType, InternalImageType]
metric = MetricType.New()

# 影像對準
registration = itk.ImageRegistrationMethodv4[InternalImageType, InternalImageType].New(
    FixedImage=fixedCastImageFilter.GetOutput(),
    MovingImage=movingCastImageFilter.GetOutput(),
    Metric=metric,
    Optimizer=optimizer,
    InitialTransform=initialTransform)

# Moving 影像所使用的初始 Transform
movingInitialTransform = TransformType.New()
movingInitialTransform.SetIdentity()
registration.SetMovingInitialTransform(movingInitialTransform)

# Fixed 影像所使用的初始 Transform
identityTransform = TransformType.New()
identityTransform.SetIdentity()
registration.SetFixedInitialTransform(identityTransform)

# 設定 Levels 相關參數
registration.SetNumberOfLevels(1)
registration.SetSmoothingSigmasPerLevel([0])
registration.SetShrinkFactorsPerLevel([1])

# 進行影像對準
registration.Update()

# 取得影像對準結果
transform = registration.GetTransform()
finalParameters = transform.GetParameters()
numberOfIterations = optimizer.GetCurrentIteration()
bestMetricValue = optimizer.GetValue()

print("Registration result: ")
print(" Iterations    = " + str(numberOfIterations))
print(" Metric value  = " + str(bestMetricValue))

Registration result: 
 Iterations    = 200
 Metric value  = -0.6965207175899221

影像轉換

做完影像對準之後，接著進行影像轉換：

# 建立對準結果的 Transform
CompositeTransformType = itk.CompositeTransform[ScalarType, Dimension]
outputCompositeTransform = CompositeTransformType.New()
outputCompositeTransform.AddTransform(movingInitialTransform)
outputCompositeTransform.AddTransform(registration.GetModifiableTransform())

# 建立輸出對準結果用的 Resampler
resampler = itk.ResampleImageFilter.New(Input=movingImage,
        Transform=outputCompositeTransform,
        UseReferenceImage=True,
        ReferenceImage=fixedImage)
resampler.SetDefaultPixelValue(100)

# 分開顯示兩張影像
fig, axs = plt.subplots(1, 2)
axs[0].imshow(itk.GetArrayViewFromImage(fixedImage), cmap='gray')
axs[0].set_title('Fixed Image')
axs[1].imshow(itk.GetArrayViewFromImage(resampler.GetOutput()), cmap='gray')
axs[1].set_title('Result Image')
plt.show()

# 將像素值型態轉為 unsigned char
OutputPixelType = itk.ctype('unsigned char')
OutputImageType = itk.Image[OutputPixelType, Dimension]
fixedIntensityRescaler = itk.RescaleIntensityImageFilter[ImageType,
        OutputImageType].New(
            Input=fixedImage,
            OutputMinimum=itk.NumericTraits[OutputPixelType].min(),
            OutputMaximum=itk.NumericTraits[OutputPixelType].max())
movingIntensityRescaler = itk.RescaleIntensityImageFilter[ImageType,
        OutputImageType].New(
            Input=movingImage,
            OutputMinimum=itk.NumericTraits[OutputPixelType].min(),
            OutputMaximum=itk.NumericTraits[OutputPixelType].max())
resultIntensityRescaler = itk.RescaleIntensityImageFilter[ImageType,
        OutputImageType].New(
            Input=resampler.GetOutput(),
            OutputMinimum=itk.NumericTraits[OutputPixelType].min(),
            OutputMaximum=itk.NumericTraits[OutputPixelType].max())

# 合併為 RGB 影像
ComponentType = itk.UC
RGBPixelType = itk.RGBPixel[ComponentType]
RGBImageType = itk.Image[RGBPixelType, Dimension]

composeFilter = itk.ComposeImageFilter[OutputImageType, RGBImageType].New()
composeFilter.SetInput(0, fixedIntensityRescaler.GetOutput())
composeFilter.SetInput(1, movingIntensityRescaler.GetOutput())
composeFilter.SetInput(2, fixedIntensityRescaler.GetOutput())

composeFilter2 = itk.ComposeImageFilter[OutputImageType, RGBImageType].New()
composeFilter2.SetInput(0, fixedIntensityRescaler.GetOutput())
composeFilter2.SetInput(1, resultIntensityRescaler.GetOutput())
composeFilter2.SetInput(2, fixedIntensityRescaler.GetOutput())

# 分開顯示兩張影像
fig, axs = plt.subplots(1, 2)
axs[0].imshow(itk.GetArrayViewFromImage(composeFilter.GetOutput()))
axs[0].set_title('Initial Image')
axs[1].imshow(itk.GetArrayViewFromImage(composeFilter2.GetOutput()))
axs[1].set_title('Result Image')
plt.show()

參考資料：ITKExamples、Qiita、Biomedical Image Analysis and Visualization: ITK