介紹如何在 Python 中使用 ITK 模組處理 TIFF 檔案序列組成的大型影像。

標準處理方式

如果伺服器的記憶體夠大，足夠容納整個影像的話，可以直接一次將整個影像讀取至記憶體，進行各種處理後，再轉存至輸出檔案中。

ITK 函式庫

這是以標準 ITK 函式庫來處理 TIFF 檔案序列影像降解析度的程式碼，程式碼結構幾乎跟 C++ 版本相同。

#!/usr/bin/python3
import itk
import glob

# 取得 TIFF 影像序列檔案名稱
imgFilenames = sorted(glob.glob("image_folder/*.tif"))

# 設定影像類型
PixelType = itk.ctype('signed short')
Dimension = 3
ImageType = itk.Image[PixelType, Dimension]

# 讀取影像
reader = itk.ImageSeriesReader[ImageType].New()
reader.SetImageIO(itk.TIFFImageIO.New())
reader.SetFileNames(imgFilenames)

# 設定 Voxel Size
changeInfImgFilter = itk.ChangeInformationImageFilter[ImageType].New()
# Voxel Size: 1.82um x 1.82um x 4um
changeInfImgFilter.SetOutputSpacing([1.82, 1.82, 4])
changeInfImgFilter.ChangeSpacingOn()
changeInfImgFilter.SetInput(reader.GetOutput())

# 降低解析度
shrinkFilter = itk.ShrinkImageFilter[ImageType, ImageType].New()
shrinkFilter.SetShrinkFactors([5, 5, 2])
shrinkFilter.SetInput(changeInfImgFilter.GetOutput())

# 輸出影像
WriterType = itk.ImageFileWriter[ImageType]
writer = WriterType.New()
writer.SetFileName("output.nrrd")
writer.SetInput(shrinkFilter.GetOutput())

# 顯示處理進度資訊
watcher1 = itk.SimpleFilterWatcher(reader, "Read Image")
watcher2 = itk.SimpleFilterWatcher(shrinkFilter, "Shrink Image")
watcher3 = itk.SimpleFilterWatcher(writer, "Write Image")

# 實際執行
writer.Update()

SimpleITK 函式庫

這是以 SimpleITK 函式庫來處理 TIFF 檔案序列影像降解析度的程式碼，效果跟上面的 ITK 版本相同，但語法簡潔許多。

#!/usr/bin/python3
import SimpleITK as sitk
import glob

# 取得 TIFF 影像序列檔案名稱
imgFilenames = sorted(glob.glob("image_folder/*.tif"))

# 讀取影像
image = sitk.ReadImage(imgFilenames)

# 設定 Voxel Size
image.SetSpacing([1.82, 1.82, 4])

# 降低解析度
image = sitk.Shrink(image, [5, 5, 2])

# 輸出影像
sitk.WriteImage(image, "output.nrrd")

批次處理單張 2D 影像

當影像的大小過大，無法直接載入至記憶體時，可改用分批處理的方式，一次處理張 TIFF 影像的降解析度，最後再將所有降解析度的 TIFF 影像合併成單一影像檔案。

#!/usr/bin/python3
import SimpleITK as sitk
import glob

# 取得 TIFF 影像序列檔案名稱
imgFilenames = sorted(glob.glob("image_folder/*.tif"))

# 暫存目錄
tmpFolder = "tmp/"

# 降解析度設定
shrinkFactor = [5, 5, 2]

# 原始 Voxel Size
voxelSize = [1.82, 1.82, 4]

# Z 軸降為 1/2 解析度
for imgFilename in imgFilenames[::shrinkFactor[2]]:
    # 讀取影像
    image = sitk.ReadImage(imgFilename)

    # 降低解析度
    image = sitk.Shrink(image, shrinkFactor[:2])

    # 輸出影像
    sitk.WriteImage(image, tmpFolder + imgFilename.replace('/', '_'))

# 取得 TIFF 影像序列檔案名稱
imgFilenames = sorted(glob.glob(tmpFolder + '*.tif'))

# 讀取影像
image = sitk.ReadImage(imgFilenames)

# 設定 Voxel Size
image.SetSpacing([a*b for a,b in zip(voxelSize, shrinkFactor)])

# 輸出影像
sitk.WriteImage(image, "output.nrrd")

壓縮 NRRD 影像檔案

SimpleITK 輸出的 NRRD 影像檔案是未壓縮的，若需要產生有壓縮的 NRRD 檔案，可以改用 nrrd 模組來輸出：

import nrrd

# 輸出壓縮 NRRD 影像
data = sitk.GetArrayViewFromImage(image)
header = {'spacings': [a*b for a,b in zip(voxelSize, shrinkFactor)]}
nrrd.write("output.nrrd", data, header)

參考資料