ITK 使用 ScalarImageKmeansImageFilter 進行 K-Means 分群教學與範例

介紹如何使用 ITK 的 ScalarImageKmeansImageFilter 以 K-means 分群演算法，將影像像素值進行分群（clustering），產生標註（label）影像。

原始影像

載入原始影像。

import itk
import matplotlib.pyplot as plt

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

# 影像維度
Dimension = 2

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

# 建立影像 Reader
ReaderType = itk.ImageFileReader[ImageType]
reader = ReaderType.New()
reader.SetFileName("fixed.mhd")
reader.Update()

# 讀取 Fixed 影像
image = reader.GetOutput()

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

區分前景與背景

使用 ScalarImageKmeansImageFilter 以 K-means 分群演算法，將影像像素值進行分群，產生標註影像。最常見的應用就是將影像的前景與背景分開：

# 建立 ScalarImageKmeansImageFilter 
kmeansFilter = itk.ScalarImageKmeansImageFilter.New()
kmeansFilter.SetInput(image)

# 設定初始群心
kmeansFilter.AddClassWithInitialMean(itk.NumericTraits[PixelType].min())
kmeansFilter.AddClassWithInitialMean(itk.NumericTraits[PixelType].max())

# 執行 KMeans
kmeansFilter.Update()
label = kmeansFilter.GetOutput()

# 顯示結果
plt.imshow(itk.GetArrayViewFromImage(label), cmap='gray')
plt.show()

ScalarImageKmeansImageFilter 所產生的標註影像維度與大小都會跟輸入影像相同，而其中的值會以 0、1、2 等數字代表像素的分群，而若要將這些整數自動擴展至輸出值的整個範圍，可以將 SetUseNonContiguousLabels 設定為 True。

細分區域

如果想要再將影像細分成更多區域，可以搭配繪製像素值分布圖來選擇更多的初始群心值。

# 繪製像素值分布圖
plt.hist(itk.GetArrayViewFromImage(image).flatten(), bins=40)
plt.show()

從像素值分布圖中可以看出像素值有三個主要的群聚區域，將這三個區域的像素值指定為初始群心，即可將影像細分為三個區域。

# 建立 ScalarImageKmeansImageFilter 
kmeansFilter = itk.ScalarImageKmeansImageFilter.New()
kmeansFilter.SetInput(image)

# 設定初始群心
kmeansFilter.AddClassWithInitialMean(0)
kmeansFilter.AddClassWithInitialMean(100)
kmeansFilter.AddClassWithInitialMean(150)

# 執行 KMeans
kmeansFilter.Update()
label = kmeansFilter.GetOutput()

# 顯示結果
plt.imshow(itk.GetArrayViewFromImage(label), cmap='gray')
plt.show()

如果希望影像標註的代碼可以依據影像區域的大小來重新排序，可以使用 RelabelImageFilter 來處理。

套疊原始與標註影像

若想要套疊原始影像與標註影像，可以使用

UCPixelType = itk.ctype('unsigned char')
UCImageType = itk.Image[UCPixelType, Dimension]

# 將標註影像轉為 unsigned char
castImageFilter = itk.CastImageFilter[ImageType, UCImageType].New()
castImageFilter.SetInput(label)

# 將原始影像轉為 unsigned char
intensityRescaler = itk.RescaleIntensityImageFilter[ImageType,
        UCImageType].New(
            Input=image,
            OutputMinimum=itk.NumericTraits[UCPixelType].min(),
            OutputMaximum=itk.NumericTraits[UCPixelType].max())

LabelType = itk.ctype('unsigned long')
LabelObjectType = itk.StatisticsLabelObject[LabelType, Dimension]
LabelMapType = itk.LabelMap[LabelObjectType]

# 將標註影像轉為 LabelMap
imageToLabelMapconverter = itk.LabelImageToLabelMapFilter[UCImageType, LabelMapType].New()
imageToLabelMapconverter.SetInput(castImageFilter.GetOutput())

# 套疊原始影像與 LabelMap
RGBImageType = itk.Image[itk.RGBPixel[itk.UC], Dimension]
overlayFilter = itk.LabelMapOverlayImageFilter[LabelMapType, UCImageType, RGBImageType].New()
overlayFilter.SetInput(imageToLabelMapconverter.GetOutput())
overlayFilter.SetFeatureImage(intensityRescaler.GetOutput())
overlayFilter.SetOpacity(0.5)

# 顯示結果
plt.imshow(itk.GetArrayViewFromImage(overlayFilter.GetOutput()))
plt.show()

參考資料：ITKExamples