介紹如何使用 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()
