Python 使用 SimpleElastix 處理 3D 影像對準教學與範例

介紹如何使用 Python 與 SimpleElastix 處理 3D 影像對準，並顯示處理結果。

Allen Mouse Brain Common Coordinate Framework

第三代的 Allen Mouse Brain Common Coordinate Framework（簡稱 CCFv3）提供了高解析度的老鼠腦圖譜，這裡我們使用 CCFv3 的腦影像作為基準，將自己的腦影像與 CCFv3 的腦影像進行影像對準。

CCFv3 的高解析度標準腦影像可以從 Allen Brain Altas 的網頁下載。

# USHORT(16 bit) anatomical template of CCFv3 - a shape and intensity average of 1675 specimen brains
wget http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/average_template/average_template_10.nrrd
wget http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/average_template/average_template_25.nrrd
wget http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/average_template/average_template_50.nrrd
wget http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/average_template/average_template_100.nrrd

# FLOAT(32 bit) reconstructed Allen Reference Atlas Nissl deformably registered to the anatomical template of CCFv3
wget http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/ara_nissl/ara_nissl_10.nrrd
wget http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/ara_nissl/ara_nissl_25.nrrd
wget http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/ara_nissl/ara_nissl_50.nrrd
wget http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/ara_nissl/ara_nissl_100.nrrd

# UINT (32 bit) structure gray matter and fiber tract annotation of CCFv3 (May 2015)
wget http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/annotation/ccf_2015/annotation_10.nrrd
wget http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/annotation/ccf_2015/annotation_25.nrrd
wget http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/annotation/ccf_2015/annotation_50.nrrd
wget http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/annotation/ccf_2015/annotation_100.nrrd

CCFv3 所提供的影像都是 nrrd 格式，可以使用 Python 的 nrrd 模組來讀取，使用前先安裝該模組：

# 安裝 pynrrd 模組
sudo pip3 install pynrrd

讀取影像之後取出切面影像觀看：

import nrrd
import matplotlib.pyplot as plt

# 讀取 NRRD 影像（平均標準腦）
avgImg, avgMeta = nrrd.read('average_template_25.nrrd');

# 抽取影像切面
sliceImg = avgImg[264,:,:]

# 顯示 2D 影像
plt.imshow(sliceImg, cmap = 'gray')
plt.show()

顯示腦區標注影像：

anoImg, anoMeta = nrrd.read('annotation_25.nrrd');
sliceImg = anoImg[264,:,:].astype(float)
plt.imshow(sliceImg, vmin = 0, vmax = 1000, cmap = 'jet')
plt.show()

基本影像比對

比對 fixed 與 moving 兩張影像，以旋轉方式修正影方向：

import SimpleITK as sitk
import matplotlib.pyplot as plt
from matplotlib.pyplot import figure
import numpy as np

# 讀取 CCFv3 標準腦與自己的腦影像
fixedITKImage = sitk.ReadImage("average_template_25.nrrd")
movingITKImage = sitk.ReadImage("my_image.nrrd")

# 轉為 NumPy 陣列
fixedImage = sitk.GetArrayViewFromImage(fixedITKImage)
movingImage = sitk.GetArrayViewFromImage(movingITKImage)

# 旋轉影像
movingImageRotation = np.rot90(movingImage, 1, (1, 2))
movingImageRotation = np.rot90(movingImageRotation, 2, (0, 1))

print("Fixed Image: {}".format(fixedImage.shape))
print("Moving Image(Rotation): {}".format(movingImageRotation.shape))

# 顯示各方向影像切面
fig, axs = plt.subplots(2, 3, figsize=(15, 10))
axs[0,0].imshow(fixedImage[fixedImage.shape[0]//2,:,:], cmap = 'jet')
axs[0,0].set_title('fixedImage(X)')
axs[0,1].imshow(fixedImage[:,fixedImage.shape[1]//2,:], cmap = 'jet')
axs[0,1].set_title('fixedImage(Y)')
axs[0,2].imshow(fixedImage[:,:,fixedImage.shape[2]//2], cmap = 'jet')
axs[0,2].set_title('fixedImage(Z)')

windowSize = 300
sliceX = movingImageRotation[movingImageRotation.shape[0]//2,:,:]
axs[1,0].imshow(sliceX, vmin = sliceX.min(), vmax = sliceX.min() + windowSize, cmap = 'jet')
axs[1,0].set_title('movingImage(X)')
sliceY = movingImageRotation[:,movingImageRotation.shape[1]//2,:]
axs[1,1].imshow(sliceY, vmin = sliceY.min(), vmax = sliceY.min() + windowSize, cmap = 'jet')
axs[1,1].set_title('movingImage(Y)')
sliceZ = movingImageRotation[:,:,movingImageRotation.shape[2]//2]
axs[1,2].imshow(sliceZ, vmin = sliceZ.min(), vmax = sliceZ.min() + windowSize, cmap = 'jet')
axs[1,2].set_title('movingImage(Z)')

plt.show()

# 輸出影像
image = sitk.GetImageFromArray(movingImageRotation)
vs = movingITKImage.GetSpacing()
image.SetSpacing([vs[0], vs[2], vs[1]])
sitk.WriteImage(image, "my_image_rotation.nrrd")

增加 Fixed 影像邊緣

增加 Fixed 影像邊緣大小：

import SimpleITK as sitk
import matplotlib.pyplot as plt

# 讀取 CCFv3 標準腦影像
fixedITKImage = sitk.ReadImage("average_template_25.nrrd")

# 取得原始影像大小
fixedImageSize = fixedITKImage.GetSize()

# 增添空白框邊
padding = 100

# 建立空白影像
expand = sitk.Image(fixedImageSize[0] + padding * 2,
                    fixedImageSize[1] + padding * 2,
                    fixedImageSize[2] + padding * 2,
                    sitk.sitkUInt16)

# 貼上原始影像
result = sitk.Paste(expand, fixedITKImage, fixedImageSize, (0, 0, 0), (padding, padding, padding))

# 轉為 NumPy 陣列
nda = sitk.GetArrayViewFromImage(result)

# 顯示影像切面
plt.imshow(nda[nda.shape[0]//2,:,:])
plt.show()

# 設定 Voxel Size
result.SetSpacing(fixedITKImage.GetSpacing())

# 輸出影像
sitk.WriteImage(result, "average_template_25_expand.nrrd")

影像對準

使用 SimpleElastix 進行影像對準：

import SimpleITK as sitk
import matplotlib.pyplot as plt
from matplotlib.pyplot import figure
import numpy as np

# 讀取 CCFv3 標準腦與自己的腦影像
fixedITKImage = sitk.ReadImage("average_template_25_expand.nrrd")
movingITKImage = sitk.ReadImage("my_image_rotation.nrrd")

# 影像資訊
print("Fixed Image:")
print("  Image Size: {}".format(fixedITKImage.GetSize()))
print("  Voxel Size: {}".format(fixedITKImage.GetSpacing()))
print("Moving Image:")
print("  Image Size: {}".format(movingITKImage.GetSize()))
print("  Voxel Size: {}".format(movingITKImage.GetSpacing()))

# 影像對準
elastixImageFilter = sitk.ElastixImageFilter()
elastixImageFilter.SetFixedImage(fixedITKImage)
elastixImageFilter.SetMovingImage(movingITKImage)
elastixImageFilter.LogToConsoleOn()

# 參數組合
#parameterMap = sitk.GetDefaultParameterMap('translation')
parameterMap = sitk.ReadParameterFile("Parameters_Rigid.txt")
elastixImageFilter.SetParameterMap(parameterMap)

# 參數組合
#parameterMap = sitk.GetDefaultParameterMap('affine')
parameterMap = sitk.ReadParameterFile("Parameters_Affine.txt")
elastixImageFilter.AddParameterMap(parameterMap)

# 參數組合
#parameterMap = sitk.GetDefaultParameterMap('bspline')
parameterMap = sitk.ReadParameterFile("Parameters_BSpline.txt")
elastixImageFilter.AddParameterMap(parameterMap)

# 執行影像對準
elastixImageFilter.Execute()

# 取得影像對準結果
resultITKImage = elastixImageFilter.GetResultImage()

# 轉為 NumPy 陣列
fixedImage = sitk.GetArrayViewFromImage(fixedITKImage)
resultImage = sitk.GetArrayViewFromImage(resultITKImage)

fig, axs = plt.subplots(2, 3, figsize=(15, 10))
axs[0,0].imshow(fixedImage[fixedImage.shape[0]//2,:,:], cmap = 'jet')
axs[0,0].set_title('fixedImage(X)')
axs[0,1].imshow(fixedImage[:,fixedImage.shape[1]//2,:], cmap = 'jet')
axs[0,1].set_title('fixedImage(Y)')
axs[0,2].imshow(fixedImage[:,:,fixedImage.shape[2]//2], cmap = 'jet')
axs[0,2].set_title('fixedImage(Y)')

windowSize = 300
sliceX = resultImage[resultImage.shape[0]//2,:,:]
axs[1,0].imshow(sliceX, vmin = sliceX.min(), vmax = sliceX.min() + windowSize, cmap = 'jet')
axs[1,0].set_title('resultImage(X)')
sliceY = resultImage[:,resultImage.shape[1]//2,:]
axs[1,1].imshow(sliceY, vmin = sliceY.min(), vmax = sliceY.min() + windowSize, cmap = 'jet')
axs[1,1].set_title('resultImage(Y)')
sliceZ = resultImage[:,:,resultImage.shape[2]//2]
axs[1,2].imshow(sliceZ, vmin = sliceZ.min(), vmax = sliceZ.min() + windowSize, cmap = 'jet')
axs[1,2].set_title('resultImage(Z)')

plt.show()

# 將像素值型態轉為 UInt8
fixedImageUint8 = sitk.Cast(sitk.IntensityWindowing(fixedITKImage, windowMinimum = 0, windowMaximum = 258, outputMinimum = 0.0, outputMaximum = 255.0), sitk.sitkUInt8)
resultImageUint8 = sitk.Cast(sitk.IntensityWindowing(resultITKImage, windowMinimum = 80, windowMaximum = 300, outputMinimum = 0.0, outputMaximum = 255.0), sitk.sitkUInt8)

# 合併影像
rgbITKImage = sitk.Cast(sitk.Compose(resultImageUint8, fixedImageUint8, resultImageUint8), sitk.sitkVectorUInt8)
rgbImage = sitk.GetArrayViewFromImage(rgbITKImage)

# 顯示 RGB 影像
fig, axs = plt.subplots(3, 1, figsize=(16, 48))

sliceX = rgbImage[rgbImage.shape[0]//2,:,:]
axs[0].imshow(sliceX)
axs[0].set_title('resultImage(X)')
sliceY = rgbImage[:,rgbImage.shape[1]//2,:]
axs[1].imshow(sliceY)
axs[1].set_title('resultImage(Y)')
sliceZ = rgbImage[:,:,rgbImage.shape[2]//2]
axs[2].imshow(sliceZ)
axs[2].set_title('resultImage(Z)')

plt.show()

# 設定影像參數
rgbITKImage.SetOrigin(fixedITKImage.GetOrigin())
rgbITKImage.SetSpacing(fixedITKImage.GetSpacing())

# 儲存檔案
sitk.WriteImage(rgbITKImage, "result_rgb.mha")

將 fixed 與 moving 兩張影像合併為一張 RGBA 彩色影像：

# 建立透明度
alpha = sitk.MaximumImageFilter().Execute(fixedImageUint8, resultImageUint8)

# 合併為 RGBA 影像
rgbaITKImage = sitk.Cast(sitk.Compose(resultImageUint8, fixedImageUint8, resultImageUint8, alpha), sitk.sitkVectorUInt8)

# 設定影像參數
rgbaITKImage.SetOrigin(fixedITKImage.GetOrigin())
rgbaITKImage.SetSpacing(fixedITKImage.GetSpacing())

# 儲存檔案
sitk.WriteImage(rgbaITKImage, "result_rgba.mha")

參考資料：Allen Institute、Office 指南、Office 指南、Office 指南、SimpleITK Notebooks