Python

Python ITK 2D 影像對準 Hello World 教學與範例

介紹如何在 Python 中使用 ITK 函式庫的 v4 影像對準框架處理二維影像的對準，並提供 hello world 範例程式碼。

ITKv4 影像對準框架

ITKv4 影像對準框架跟上一版比較起來更為彈性，主要差異就是對準過程都在新的虛擬影像（virtual image）中進行，固定影像（fixed image）與調動影像（moving image）都透過轉換（transforms）與內插（interpolators）在虛擬影像空間以 metric 衡量配適程度，然後將結果交給 optimizer 更新轉換參數，重複這個過程直到結果收斂。

ITKv4 影像對準框架

在這張 ITKv4 影像對準框架圖中，虛線的項目代表資料物件（data objects），實線的項目代表處理物件（process objects）。

設定影像類型

在 ITKv4 的架構下要進行影像對準時，首先要決定影像的類型，包含像素的類型以及影像的維度。

import itk
import matplotlib.pyplot as plt

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

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

# 影像維度
Dimension = 2

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

讀取影像

根據影像類型建立影像 reader，並且讀取固定影像與調動影像。

# 建立影像 Reader
FixedReaderType = itk.ImageFileReader[FixedImageType]
fixedImageReader = FixedReaderType.New()
fixedImageReader.SetFileName("BrainProtonDensitySliceBorder20.mhd")
fixedImageReader.Update()

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

# 建立影像 Reader
MovingReaderType = itk.ImageFileReader[MovingImageType]
movingImageReader = MovingReaderType.New()
movingImageReader.SetFileName("BrainProtonDensitySliceShifted13x17y.mhd")
movingImageReader.Update()

# 讀取 Moving 影像
movingImage = movingImageReader.GetOutput()

這裡影像 reader 所指定的影像類型是指讀取影像之後所輸出的影像類型，如果這個類型跟影像檔案中實際的資料不符，影像 reader 會自動處理影像類型的轉換。

若要查看影像內部的詳細資訊，可以直接用 print 輸出影像物件：

# 查看影像詳細資訊
print(fixedImage)

Image (0x25a89c0)
  RTTI typeinfo:   itk::Image<float, 2u>
  Reference Count: 2
  Modified Time: 246
  Debug: Off
  Object Name:
  Observers:
    none
  Source: (0x2960380)
  Source output name: Primary
  Release Data: Off
  Data Released: False
  Global Release Data: Off
  PipelineMTime: 59
  UpdateMTime: 247
  RealTimeStamp: 0 seconds
  LargestPossibleRegion:
    Dimension: 2
    Index: [0, 0]
    Size: [221, 257]
  BufferedRegion:
    Dimension: 2
    Index: [0, 0]
    Size: [221, 257]
  RequestedRegion:
    Dimension: 2
    Index: [0, 0]
    Size: [221, 257]
  Spacing: [1, 1]
  Origin: [0, 0]
  Direction:
1 0
0 1

  IndexToPointMatrix:
1 0
0 1

  PointToIndexMatrix:
1 0
0 1

  Inverse Direction:
1 0
0 1

  PixelContainer:
    ImportImageContainer (0x3dd22d0)
      RTTI typeinfo:   itk::ImportImageContainer<unsigned long, float>
      Reference Count: 1
      Modified Time: 244
      Debug: Off
      Object Name:
      Observers:
        none
      Pointer: 0x3fed990
      Container manages memory: true
      Size: 56797
      Capacity: 56797

使用 matplotlib 的 imshow 顯示兩張影像的內容。

# 分開顯示兩張影像
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()

固定影像與調動影像

影像對準

進行影像對準時，需要定義好 transform（影像轉換，例如線性或非線性轉換）、optimizer（尋找最佳解的方法）、metric（衡量兩張影像距離的方法）、interpolator（計算內差的方法）的類型，接著將各元件串接起來，進行影像對準。

# 影像轉換類型
TransformType = itk.TranslationTransform[ScalarType, Dimension]

# 建立初始轉換
initialTransform = TransformType.New()

# 建立 Optimizer 物件
OptimizerType = itk.RegularStepGradientDescentOptimizerv4[ScalarType]
optimizer = OptimizerType.New()

# 建立 Metric 物件
MetricType = itk.MeanSquaresImageToImageMetricv4[FixedImageType, MovingImageType]
metric = MetricType.New()

# 影像對準物件
RegistrationType = itk.ImageRegistrationMethodv4[FixedImageType, MovingImageType]
registration = RegistrationType.New()

# 設定 Metric
registration.SetMetric(metric)

# 設定 Optimizer
registration.SetOptimizer(optimizer)

# 設定初始轉換
registration.SetInitialTransform(initialTransform)

# 設定輸入影像
registration.SetFixedImage(fixedImage)
registration.SetMovingImage(movingImage)

# Interpolator 類型
FixedLinearInterpolatorType = itk.LinearInterpolateImageFunction[FixedImageType, ScalarType]
MovingLinearInterpolatorType = itk.LinearInterpolateImageFunction[FixedImageType, ScalarType]

# 建立 Interpolator 物件
fixedInterpolator = FixedLinearInterpolatorType.New()
movingInterpolator = MovingLinearInterpolatorType.New()

# 設定 Interpolator
metric.SetFixedInterpolator(fixedInterpolator)
metric.SetMovingInterpolator(movingInterpolator)

設定調動影像的初始轉換，先取得預設參數，檢查參數數量。

# Moving 影像所使用的初始轉換
movingInitialTransform = TransformType.New()
initialParameters = movingInitialTransform.GetParameters()

# 查看參數數量
print(initialParameters.GetNumberOfElements())

設定初始轉換的參數，此處兩個參數分別代表 X 與 Y 軸方向的位移。

# 設定初始轉換參數
initialParameters[0] = 0
initialParameters[1] = 0
movingInitialTransform.SetParameters(initialParameters)
registration.SetMovingInitialTransform(movingInitialTransform)

固定影像所使用的初始轉換通常都是 identity（維持不變），由於預設值就是 identity，所以亦可省略不寫：

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

optimizer 的 learning rate 是非常重要的參數，其指定 optimizer 每次的步長（step length），不適當的 learning rate 會讓影像對準無法收斂。

learning rate 如果設定太大，會造成 optimizer 不穩定，若設定太小則會讓收斂過程非常緩慢，簡易的設定方式是一開始先設定為較小的數值（例如 1.0 到 5.0 左右），等到其他參數都調整好、整體都可以穩定收斂之後，再慢慢升高 learning rate 直到開始出現收斂不穩定的狀況為止，最佳的 learning rate 就是讓疊代次數最少，但依然維持收斂的穩定性。

learning rate 在實際運算時會跟 metric 的梯度（gradient）相乘產生步長，所以步長也會跟 metric 的數值成正比，在 metric 數值尺度比較大的狀況下，要維持相同的 optimizer 行為，就需要較小的 learning rate。

# 設定 Learning Rate
optimizer.SetLearningRate(4)

每當導數（derivative）方向發生劇烈改變時，通常就是代表剛越過局部極值點（local extrema），此時 optimizer 就會降低 learning，而降低 learning rate 的比例就是由 relaxation factor 參數來控制，這個參數值介於 0 到 1 之間，預設值是 0.5。

# 設定 Relaxation Factor
optimizer.SetRelaxationFactor(0.5)

隨著影像對準的過程，learning rate 會漸漸減小，當其減小到非常小的值的時候，就代表已經收斂了，這個判斷收斂的門檻值可以透過 minimum step length 參數來設定，預設值為 1e-4。

# 設定 Minimum Step Length
optimizer.SetMinimumStepLength(0.001)

在疊代的過程也有可能會出現無法收斂的狀況，因此會設定疊代次數的上限值。

# 設定最高疊代次數
optimizer.SetNumberOfIterations(200)

ITKv4 支援多層級的影像對準（multi-level registration），在不同階段中採用不同解析度的虛擬空間（virtual space），並搭配不同的平滑度（smoothness）的固定影像與調動影像，而這些條件必須在影像對準計算開始之前就設定好，否則就只能採用預設值進行影像對準。

這裡我們採用最簡單的單一層級對準方式，亦不使用任何的降解析度（shrinking）或平滑化（smoothing）。

# 設定 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 = 13.001531179313957
 Translation Y = 16.99948401485064
 Iterations    = 21
 Metric value  = 0.0012263178458750189

TranslationTransform 轉換的參數意義非常單純，兩個參數就是代表 X 軸與 Y 軸的位移。這裡的 Iterations 就是實際計算所使用的疊代次數，這個數字如果很大，可能代表 learning rate 設定太小，會造成計算時間拉長，或是結果尚未收斂。

轉換影像

在影像對準計算完成之後，只是得到了轉換函數的各個參數值，如果要得到調動影像轉換後的結果，必須把轉換函數套用至調動影像上進行轉換。

由於此處我們有指定調動影像的初始轉換，因此這裡必須以 CompositeTransform 將初始轉換與最終計算出來的轉換結合，再套用至調動影像上，才能正確將調動影像轉換至固定影像的空間中。

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

轉換影像的工作必須使用 resampler 搭配轉換函數來處理：

# 建立輸出對準結果用的 Resampler
ResampleFilterType = itk.ResampleImageFilter
resampler = ResampleFilterType[MovingImageType, FixedImageType].New()
resampler.SetInput(movingImage)
resampler.SetTransform(outputCompositeTransform)
resampler.SetUseReferenceImage(True)
resampler.SetReferenceImage(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()