介紹如何在 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()
![基準影像(fixed image)與調動影像(moving image) 基準影像(fixed image)與調動影像(moving image)](https://officeguide.cc/wp-content/uploads/2020/10/python-itk-linear-translation-registration-tutorial-examples-20201014-01.png)
平移線性對準(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()
![影像對準結果 影像對準結果](https://officeguide.cc/wp-content/uploads/2020/10/python-itk-linear-translation-registration-tutorial-examples-20201014-02.png)
將灰階影像合併成 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()
![影像對準結果 影像對準結果](https://officeguide.cc/wp-content/uploads/2020/10/python-itk-linear-translation-registration-tutorial-examples-20201014-03.png)
如果需要觀察影像對準過程中的 metric 數值收斂狀況,可以透過 observer 設定回呼函數來處理。
參考資料:ITKExamples