Python

SimpleITK 使用 IntensityWindowingImageFilter 轉換影像數值範圍教學與範例

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