SimpleITK 計算 2D 二元遮罩影像物件數量與統計值

介紹如何在 Python 中使用 SimpleITK 的 ConnectedComponent 取出二元遮罩影像中的物件，計算物件數量與統計值。

我們希望透過影像分割（segmentation）之後的遮罩影像（mask image），自動計算圖片中的細胞數量（或是各種不相連物件的數量）。

讀取原始影像

首先讀取含有大量細胞的原始影像，並且使用 matplotlib 查看影像內容：

import SimpleITK as sitk
import matplotlib.pyplot as plt

# 讀取原始影像
image = sitk.ReadImage("ExampleHuman/images/AS_09125_050116030001_D03f00d0.tif")

# 顯示影像
plt.imshow(sitk.GetArrayViewFromImage(image), cmap = 'gray')
plt.show()

範例中所使用的圖片可從 CellProfiler 網站下載。

產生遮罩影像

這裡我們只使用最簡單的門檻值方式，產生測試用的遮罩影像，實務上可以使用更進階的影像分割方式產生更高品質的遮罩影像。

# 以門檻值方式產生遮罩影像
maskImage = sitk.BinaryThreshold(image, lowerThreshold=50, upperThreshold=255, insideValue=1, outsideValue=0)

# 顯示影像
plt.imshow(sitk.GetArrayViewFromImage(maskImage), cmap = 'gray')
plt.show()

產生物件標註影像

有了遮罩影像之後，再以 ConnectedComponentImageFilter 每一個沒有連接的物件分別標示不同的 ID 值，產生標註影像：

# 產生物件標註影像
ccFilter = sitk.ConnectedComponentImageFilter()
ccImage = ccFilter.Execute(maskImage)

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

Object Count: 351

# 顯示影像
plt.imshow(sitk.GetArrayViewFromImage(ccImage), cmap = 'gray')
plt.show()

彩色影像顯示標註

由 ConnectedComponentImageFilter 所產生的標註影像主要是給程式處理用的，若要以人眼查看這種影像，可以用 LabelToRGBImageFilter 將原始標註影像轉為 RGB 影像：

# 將標註影像轉為 RGB 彩色影像
rgbCCImage = sitk.LabelToRGB(ccImage)

# 顯示影像
plt.imshow(sitk.GetArrayViewFromImage(rgbCCImage))
plt.show()

計算各物件的影像統計資訊

有了物件標註影像之後，就可以用 LabelIntensityStatisticsImageFilter 根據標註影像計算各物件的大小與像素平均值：

# 計算各物件的影像統計資訊
stats = sitk.LabelIntensityStatisticsImageFilter()
stats.Execute(ccImage, image)

# 列出所有物件統計值
for l in stats.GetLabels():
    print("Label {0}: Mean = {1}, Size = {2}".format(l, stats.GetMean(l), stats.GetPhysicalSize(l)))

Label 1: Mean = 55.08571428571429, Size = 4.35582561728395
Label 2: Mean = 61.5, Size = 0.7467129629629629
Label 3: Mean = 63.85981308411215, Size = 13.316381172839504
Label 4: Mean = 64.84615384615384, Size = 1.6178780864197528
Label 5: Mean = 65.78181818181818, Size = 6.844868827160493
Label 6: Mean = 95.85915492957747, Size = 17.672206790123454
Label 7: Mean = 56.723684210526315, Size = 9.458364197530862
Label 8: Mean = 56.72549019607843, Size = 6.3470601851851844
Label 9: Mean = 56.36842105263158, Size = 9.458364197530862
Label 10: Mean = 52.567567567567565, Size = 4.604729938271604
[略]

取出較大的物件

通常在偵測物件時，會有一些小的雜訊參雜在裡面，若要去除這些雜訊，最簡單的做法就是依照物件的大小排序，只取出體積比較大的物件。

# 依照物件大小重新編排物件標註
relabelFilter = sitk.RelabelComponentImageFilter()
relabelImage = relabelFilter.Execute(ccImage)

# 計算各物件的影像統計資訊
stats2 = sitk.LabelIntensityStatisticsImageFilter()
stats2.Execute(relabelImage, image)

# 列出所有物件統計值
for l in stats2.GetLabels():
    print("Label {0}: Mean = {1}, Size = {2}".format(l, stats2.GetMean(l), stats2.GetPhysicalSize(l)))

Label 1: Mean = 87.69310344827586, Size = 36.09112654320987
Label 2: Mean = 60.44268774703557, Size = 31.486396604938268
Label 3: Mean = 79.988, Size = 31.113040123456784
Label 4: Mean = 105.67782426778243, Size = 29.744066358024686
Label 5: Mean = 72.3225806451613, Size = 27.00611882716049
Label 6: Mean = 88.90760869565217, Size = 22.899197530864193
Label 7: Mean = 69.87078651685393, Size = 22.15248456790123
Label 8: Mean = 69.48795180722891, Size = 20.659058641975307
Label 9: Mean = 87.46060606060605, Size = 20.53460648148148
Label 10: Mean = 92.73291925465838, Size = 20.036797839506168
[略]

參考資料：StackOverflow