介紹如何使用 SimpleITK 的 IntensityWindowingImageFilter 對影像的像素值進行線性轉換。
IntensityWindowingImageFilter 的作用跟 RescaleIntensityImageFilter 類似，差異在於 IntensityWindowingImageFilter 可以自行指定窗函數的數值區間。
讀取 DICOM 影像
這裡我們使用 pydicom 模組所附帶的 DICOM 影像檔案來做示範。首先透過 pydicom 的 get_testdata_files 取得 CT_small.dcm 這個測試檔案路徑之後，再以 SimpleITK 來讀取檔案中的影像。
import SimpleITK as sitk import matplotlib.pyplot as plt from pydicom.data import get_testdata_files # 取得 Pydicom 附帶的 DICOM 測試影像路徑 filename = get_testdata_files('CT_small.dcm')[0] # 讀取 DICOM 影像 image = sitk.ReadImage(filename)
若想要使用 pydicom 來讀取 DICOM 檔案，可以參考 Python 使用 Pydicom 讀取、編輯 DICOM 影像檔案教學與範例。
檢查原始影像資料型態：
# 查看原始影像資料型態 print(image.GetPixelIDTypeAsString())
32-bit signed integer
這是一張 CT 電腦斷層掃描的影像，在 DICOM 檔案中每個像素的原始數值使用了兩個位元組（bytes）來儲存，但是以 SimpleITK 讀進來之後，變成了四個位元組的整數。
使用 StatisticsImageFilter 計算基本的影像統計資訊，取得像素值的分布範圍：
# 計算影像像素值範圍 statsFilter = sitk.StatisticsImageFilter() statsFilter.Execute(image) print("Minimum:", statsFilter.GetMinimum()) print("Maximum:", statsFilter.GetMaximum())
Minimum: -896.0 Maximum: 1167.0
使用 matplotlib 顯示原始影像的內容與像素分布。
# 顯示影像與像素值分布 nda = sitk.GetArrayViewFromImage(image) fig, axs = plt.subplots(1, 2, figsize=(8, 4)) axs[0].imshow(nda[0,:,:], cmap = 'gray') axs[0].set_title('Image') axs[1].hist(nda.flatten(), bins=30) axs[1].set_title('Histogram') plt.show()
窗函數與線性轉換
如果我們只對影像中部份的像素數值範圍有興趣，可以使用 IntensityWindowingImageFilter 以窗函數將指定的範圍取出來，並以線性轉換轉換至新的範圍，例如將影像數值範圍從 -200 至 300 轉換為 0 至 255：
# 將影像數值範圍從 -200 至 300 轉換為 0 至 255 winFilter = sitk.IntensityWindowingImageFilter() winFilter.SetWindowMinimum(-200) winFilter.SetWindowMaximum(300) winFilter.SetOutputMinimum(0) winFilter.SetOutputMaximum(255) imageRescale = winFilter.Execute(image)
在轉換時若遇到落在範圍以下的數值，會直接轉換為輸出範圍的最小值（OutputMinimum，此處為 0），超過範圍的數值則會轉為輸出範圍的最大值（OutputMaximum，此處為 255）。
計算轉換後的影像像素值範圍：
# 計算影像像素值範圍 statsFilter = sitk.StatisticsImageFilter() statsFilter.Execute(imageRescale) print("Minimum:", statsFilter.GetMinimum()) print("Maximum:", statsFilter.GetMaximum())
Minimum: 0.0 Maximum: 255.0
通常經過轉換之後，會將影像類型轉為較精簡的型態，節省記憶體空間：
# 將影像類型轉換為 8 位元的無號整數
castFilter = sitk.CastImageFilter()
castFilter.SetOutputPixelType(sitk.sitkUInt8)
imageUint8 = castFilter.Execute(imageRescale)
顯示轉換後的影像與像素值分布
# 顯示影像與像素值分布 nda = sitk.GetArrayViewFromImage(imageUint8) fig, axs = plt.subplots(1, 2, figsize=(8, 4)) axs[0].imshow(nda[0,:,:], cmap = 'gray') axs[0].set_title('Image') axs[1].hist(nda.flatten(), bins=30) axs[1].set_title('Histogram') plt.show()
