Python ITK 影像平移線性對準 Translation Registration 教學與範例

介紹如何在 Python 中使用 ITK 進行 2D 影像平移線性對準（translation registration）。

準備測試用影像

這裡我們使用一張影像作為基準影像（fixed image），以 TranslationTransform 平移轉換套用至基準影像上，產生一張調動影像（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()

# 建立 TranslationTransform 平移轉換
translationTransform = itk.TranslationTransform[ScalarType, Dimension].New()
parameters = translationTransform.GetParameters()
parameters[0] = 29 # X 軸平移
parameters[1] = 14 # Y 軸平移
translationTransform.SetParameters(parameters)

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

# 設定 Resampler 各項參數
resampleFilter.SetInput(fixedImage)
resampleFilter.SetTransform(translationTransform)
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()

平移線性對準（Translation 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.TranslationTransform[ScalarType, Dimension]
initialTransform = TransformType.New()

# 影像對準所使用的 Optimizer
optimizer = itk.RegularStepGradientDescentOptimizerv4.New(
        LearningRate=4,
        MinimumStepLength=0.001,
        RelaxationFactor=0.5,
        NumberOfIterations=200)

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

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

# Moving 影像所使用的初始 Transform
movingInitialTransform = TransformType.New()
initialParameters = movingInitialTransform.GetParameters()
initialParameters[0] = 0
initialParameters[1] = 0
movingInitialTransform.SetParameters(initialParameters)
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()
translationAlongX = finalParameters.GetElement(0)
translationAlongY = finalParameters.GetElement(1)
numberOfIterations = optimizer.GetCurrentIteration()
bestMetricValue = optimizer.GetValue()

print("Registration result: ")
print(" Translation X = " + str(translationAlongX))
print(" Translation Y = " + str(translationAlongY))
print(" Iterations    = " + str(numberOfIterations))
print(" Metric value  = " + str(bestMetricValue))

Registration result: 
 Translation X = -28.999121220546765
 Translation Y = -13.999008455600247
 Iterations    = 27
 Metric value  = 0.00020581286412778255

影像對準所支援的影像類型主要取決於 metric 的支援類型，可用的類型可以使用 GetTypes() 取得，例如：

# 取得 MeanSquaresImageToImageMetricv4 支援的類型
itk.MeanSquaresImageToImageMetricv4.GetTypes()

對於開發者而言，若要將其他的類型加入其中，可以在編譯 ITK 函式庫的時候，將對應類型的 CMake 參數設定為 ON，就可以新增影像的類型，例如（ITK_WRAP_unsigned_short、ITK_WRAP_unsigned_char）；若要新增影像維度，則可將新的維度加入 ITK_WRAP_IMAGE_DIMS 之中。

影像轉換

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

# 建立對準結果的 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()

將灰階影像合併成 RGB 彩色影像，比較影像對準的結果。

# 將像素值型態轉為 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()

如果需要觀察影像對準過程中的 metric 數值收斂狀況，可以透過 observer 設定回呼函數來處理。

參考資料：ITKExamples