ITK 計算 3D 二元遮罩影像物件數量與體積

介紹如何在 Python 中使用 ITK 函式庫計算 3D 二元遮罩影像中的物件數量與體積。

讀取原始影像

載入原始的 3D 細胞影像：

import itk
import itkwidgets

# 載入影像
image = itk.imread("3d_monolayer_xy1_ch2.tif")

# 查看影像
itkwidgets.view(image)

這裡所使用的 3D 細胞影像可以從 CellProfiler 的 GitHub 網站上取得。

產生遮罩影像

使用 BinaryThresholdImageFilter 以門檻值方式產生遮罩影像：

# 以門檻值方式產生遮罩影像
binFilter = itk.BinaryThresholdImageFilter[ImageType, ImageType].New()

# 設定門檻值
binFilter.SetLowerThreshold(110)
binFilter.SetUpperThreshold(255)

binFilter.SetInsideValue(255)
binFilter.SetOutsideValue(0)
binFilter.SetInput(image)
binFilter.Update()

# 查看影像
itkwidgets.view(binFilter.GetOutput())

產生物件標註影像

使用 ConnectedComponentImageFilter 將二元影像遮罩轉為物件標註：

# 產生物件標註影像
ccFilter = itk.ConnectedComponentImageFilter[ImageType, CCImageType].New()
ccFilter.SetInput(binFilter.GetOutput())

# 只計算局部區域時可用遮罩控制
#ccFilter.SetMaskImage(anotherMaskImage)

# 計算寬度為 1 像素大小的物件時，需要啟用 Fully Connected 選項
#ccFilter.SetFullyConnected(True)

ccFilter.Update()

# 查看影像
itkwidgets.view(ccFilter.GetOutput())

查看所有物件的數量：

# 物件數量
print("Object Count:", ccFilter.GetObjectCount())

Object Count: 4320

由於影像中有許多雜訊都會被判定為獨立的物件，因此在初步的處理結果中，物件的數量會遠大於實際的物件數量，必須再進一步處理。

保留前 N 個最大的物件

LabelShapeKeepNObjectsImageFilter 可以根據各種物件的屬性（預設是像素數），取出前 N 個最大的物件：

# 保留前 N 個最大的物件
keepNObjFilter = itk.LabelShapeKeepNObjectsImageFilter[CCImageType].New()
keepNObjFilter.SetInput(ccFilter.GetOutput())
keepNObjFilter.SetBackgroundValue(0)

# 保留前 20 個最大的物件
keepNObjFilter.SetNumberOfObjects(20)

# 根據物件的像素數量排序（NUMBER_OF_PIXELS = 100）
keepNObjFilter.SetAttribute(100)

keepNObjFilter.Update()

# 查看影像
itkwidgets.view(keepNObjFilter.GetOutput())

依照物件大小重新編排標註

另一種常見的後處理方式是使用 RelabelComponentImageFilter 依照物件大小重新編排物件標註，並且排除太小的物件：

# 依照物件大小重新編排物件標註
relabelFilter = itk.RelabelComponentImageFilter[CCImageType, CCImageType].New()
relabelFilter.SetInput(ccFilter.GetOutput())

# 設定物件的最低像素數
relabelFilter.SetMinimumObjectSize(1000)

relabelFilter.Update()

# 物件數量
print("Number of Object:", relabelFilter.GetNumberOfObjects())

Number of Object: 18

列出重新編排後的所有物件：

# 列出所有物件的大小
for i in range(1, relabelFilter.GetNumberOfObjects()+1):
    print("Label {}: Size = {}".format(i, relabelFilter.GetSizeOfObjectInPixels(i)))

Label 1: Size = 96923
Label 2: Size = 87725
Label 3: Size = 44710
Label 4: Size = 44152
Label 5: Size = 44129
Label 6: Size = 42240
Label 7: Size = 39144
Label 8: Size = 36027
Label 9: Size = 33964
Label 10: Size = 32491
Label 11: Size = 31588
Label 12: Size = 30674
Label 13: Size = 29314
Label 14: Size = 21026
Label 15: Size = 12549
Label 16: Size = 4779
Label 17: Size = 3656
Label 18: Size = 2733

這裡透過 GetSizeOfObjectInPixels 所取得的體積其單位是 voxel，若要獲取物件實際的體積，可以改用 GetSizeOfObjectInPhysicalUnits。

串流處理

對於大型影像，可以使用 StreamingImageFilter 串流分批處理的方式處理。

import itk
import itkwidgets

PixelType = itk.UC
Dimension = 3
ImageType = itk.Image[PixelType, Dimension]
CCImageType = itk.Image[itk.US, Dimension]

# 載入影像
imageFilename = "3d_monolayer_xy1_ch2.tif"
reader = itk.ImageFileReader[ImageType].New(FileName=imageFilename)

# 以門檻值方式產生遮罩影像
binFilter = itk.BinaryThresholdImageFilter[ImageType, ImageType].New()
binFilter.SetLowerThreshold(110)
binFilter.SetUpperThreshold(255)
binFilter.SetInsideValue(255)
binFilter.SetOutsideValue(0)
binFilter.SetInput(reader.GetOutput())

# 產生物件標註影像
ccFilter = itk.ConnectedComponentImageFilter[ImageType, CCImageType].New()
ccFilter.SetInput(binFilter.GetOutput())

# 依照物件大小重新編排物件標註
relabelFilter = itk.RelabelComponentImageFilter[CCImageType, CCImageType].New()
relabelFilter.SetInput(ccFilter.GetOutput())
relabelFilter.SetMinimumObjectSize(1000)

# 設定分割數量
numberOfSplits = 4

# 建立 PipelineMonitorImageFilter 查看 Streaming 每次處理的範圍
monitorFilter = itk.PipelineMonitorImageFilter[CCImageType].New()
monitorFilter.SetInput(relabelFilter.GetOutput())

# 建立 StreamingImageFilter 將資料以 Streaming 分批處理
streamingFilter = itk.StreamingImageFilter[CCImageType, CCImageType].New()
streamingFilter.SetInput(monitorFilter.GetOutput())
streamingFilter.SetNumberOfStreamDivisions(numberOfSplits)

# 實際執行
streamingFilter.Update()

# 輸出影像完整的 Region 資訊
print('The output LargestPossibleRegion is: ' +
      str(streamingFilter.GetOutput().GetLargestPossibleRegion()))

# 輸出 Streaming 的分批 Region 資訊
updatedRequestedRegions = monitorFilter.GetUpdatedRequestedRegions()
print("Updated RequestedRegion's:")
for r in range(len(updatedRequestedRegions)):
    print('  ' + str(updatedRequestedRegions[r]))

The output LargestPossibleRegion is: itkImageRegion3([0, 0, 0], [256, 256, 60])
Updated RequestedRegion's:
  itkImageRegion3([0, 0, 0], [256, 256, 15])
  itkImageRegion3([0, 0, 15], [256, 256, 15])
  itkImageRegion3([0, 0, 30], [256, 256, 15])
  itkImageRegion3([0, 0, 45], [256, 256, 15])

# 物件數量
print("Number of Object:", relabelFilter.GetNumberOfObjects())

Number of Object: 18

# 列出所有物件的大小
for i in range(1, relabelFilter.GetNumberOfObjects()+1):
    print("Label {}: Size = {}".format(i, relabelFilter.GetSizeOfObjectInPixels(i)))

Label 1: Size = 96923
Label 2: Size = 87725
Label 3: Size = 44710
Label 4: Size = 44152
Label 5: Size = 44129
Label 6: Size = 42240
Label 7: Size = 39144
Label 8: Size = 36027
Label 9: Size = 33964
Label 10: Size = 32491
Label 11: Size = 31588
Label 12: Size = 30674
Label 13: Size = 29314
Label 14: Size = 21026
Label 15: Size = 12549
Label 16: Size = 4779
Label 17: Size = 3656
Label 18: Size = 2733