介紹如何在 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