CT图像预处理之窗宽窗位调整

孟想成真

目录

1 CT图像的Hu值变换
2 CT图像窗宽窗位调整

1 CT图像的Hu值变换

CT值属于医学领域的概念，通常称亨氏单位（hounsfield unit,HU),反映了组织对X射线的吸收程度。黑影表示低吸收区，即低密度区，如含气体多的肺部；白影表示高吸收区，即高密度区，如骨骼。但是CT 和X线图像相比，CT的密度分辨力高，即有高的密度分辨力（density resolution）。因此，人体软组织的密度差别虽小。吸收系数虽多接近于水，也能形成对比而成像。解剖成像。
灰度值属于计算机领域的概念，指的是单个像素点的亮度，灰度值越大表示越亮。范围一般从0到255，白色为255，黑色为0，故黑白图片也称为灰度图像。
两种最常用的医学图像数据格式：DICOM类型和NIFTI类型
DICOM类型的文件后缀为dcm,常用的python处理库有：SimpleITK,pydicom库
NIFTI类型的文件后缀为nii或nii.gz，常用的python处理库有：SimpleITK,nibabel库
   注：无论对于dcm还是nii格式的图片，只要是ct图，为了更好地观察不同器官，都可以选择将储存的原始数据转化为Hu值，因为Hu值即代表了物体真正的密度。
         对于nii格式的图片，SimpleITK,nibabel中常用的api接口，都会自动的进行上述转化过程，即取出来的值已经是Hu了。除非专门用nib.load('xx').dataobj.get_unscaled()或者itk.ReadImage('xx').GetPixel(x,y,z)才能取得原始数据
         对于dcm格式的图片，SimpleITK,pydicom常用的api接口都不会将原始数据自动转化为Hu。（itk snap软件读入dcm或nii都不会对数据进行scale操作）
  公式与代码：
HU = pixel_val*slope + intercept
其中，slope,intercept可以从元数据中读取

1. def get_pixels_hu(ct_array,slope,intercept):
2.     if slope != 1:
3.         ct_array = slope*ct_array
4.     ct_array += intercept
5.     return ct_array

复制代码

2 CT图像窗宽窗位调整

对于CT图像来说，它的HU值从负几千到正几千都有，我们可以使用ITK依次读取CT图像的每个值，并可以将其每个值保存在Mat中，但这时通过opencv的imshow函数不能够显示出来。这时图像显示和打印面临的一个问题是：图像的亮度和对比度不能充分突出关键部分。这里所指的“关键部分”在 CT 里的例子有软组织、骨头、脑组织、肺、腹部等等。因为显示器往往只有 8-bit， 而数据有 12- 至 16-bits。也就是说，不能直接将 CT 值作为传统意义上的像素值进行输出， 而是要经过一个变换，或者称其为映射。
针对这些问题，研究人员提出一些要求 (requirements)，变换应该尽量满足：
要求一：充分利用 0-255 间的显示有效值域
要求二：尽量减少值域压缩带来的损失
要求三：不能损失应该突出的组织部分
医学图像领域的窗口技术，包括窗宽(window width)和窗位(window center),用于选择感兴趣的CT值范围。因为各种组织结构或病变具有不同的CT值，因此欲显示某一组织结构细节时，应选择适合观察该组织或病变的窗宽和窗位，以获得最佳显示。窗宽和窗位分别对应图像的对比度与亮度。
注：1）窗宽窗位是CT图像特有的概念，MRI图像中可没有此概念2）CT图像必须西安转换成HU值再做窗宽窗位调整
   窗宽是CT图像上显示的CT值范围，在此CT值范围内的组织和病变均以不同的模拟灰度显示。而CT值高于此范围的组织和病变，无论高出程度有多少，均以白影显示，不再有灰度差异; 反之，低于此范围的组织结构，不论低的程度有多少，均以黑影显示，也无灰度差别。增大窗宽，则图像所示CT值范围加大，显示具有不同密度的组织结构增多，但各结构之间的灰度差别减少。减小窗宽，则显示的组织结构减少，然而各结构之间的灰度差别增加。如观察脑质的窗宽常为-15~+85H，即密度在-15 ~+85H范围内的各种结构如脑质和脑脊液间隙均以不同灰度显示。而高于+85H的组织结构如骨质几颅内钙化，其间虽有密度差，但均以白影显示，无灰度差别;而低于-15H组织结构如皮下脂肪及乳突内气体均以黑影显示，其间也无灰度差别。
  窗位是窗的中心位置，同样的窗宽，由于窗位不同，其所包括的CT值范围的CT值也有差异。例如窗宽同为100H，当窗位为0H时，其CT值范围为-50 ~ +50H ; 如窗位为+35H时，则CT值范围为-15~+85H。通常，欲观察某以组织结构及发生的病变，应以该组织的CT值为窗位。例如脑质CT值约为+35H，则观察脑组织及其病变时，选择窗位以+35H为妥。
  设定正确的窗宽窗位
应用场景：用神经网络分割liver 和 liver tumor,对训练数据进行预处理。
基本思路：
方法一：
窗位、窗宽可以设置如下：(但有研究者认为该方法没有方法二的量身定做有优势）

1. def get_window_size(organ_name):
2.     if organ_name == 'lung':
3.         # 肺部 ww 1500-2000 wl -450--600
4.         center = -500
5.         width = 2000
6.     elif organ_name == 'abdomen':
7.         # 腹部 ww 300-500 wl 30-50
8.         center = 40
9.         width = 500
10.     elif organ_name == 'bone':
11.         # 骨窗 ww 1000-1500 wl 250-350
12.         center = 300
13.         width = 2000
14.     elif organ_name == 'lymph':
15.         # 淋巴、软组织 ww 300-500 wl 40-60
16.         center = 50
17.         width = 300
18.     elif organ_name == 'mediastinum':
19.         # 纵隔 ww 250-350 wl 30-50
20.         center = 40
21.         width = 350
22.     elif organ_name == 'liver':
23.         #正常肝组织的CT值为50HU左右，
24.         #采用 ww 100-150 wl 70-75，可以检查出常规ww200,wl 50没有检测出的病灶
25.         #两种组织密度差别小的病灶适合用高窗位和低窗宽来观察
26.         center == 70
27.         width == 100
28.    
29.     return center, width

复制代码

方法二：
1）用liver和tumor的mask抠出liver 和 tumor部分
2）统计相应的最大值和最小值，并用其计算窗宽窗位。
优点：每个case的窗宽窗位依据案例量身定做。相较于整个数据集用统一的窗宽窗位，或者采用cut off,有无可比拟的优点。

1. import SimpleITK as sitk
2. import numpy as np
3. import os
5. def saved_preprocessed(savedImg,origin,direction,xyz_thickness,saved_name):
6.     newImg = sitk.GetImageFromArray(savedImg)
7.     newImg.SetOrigin(origin)
8.     newImg.SetDirection(direction)
9.     newImg.SetSpacing((xyz_thickness[0], xyz_thickness[1], xyz_thickness[2]))
10.     sitk.WriteImage(newImg, saved_name)
12. def window_transform(ct_array, windowWidth, windowCenter, normal=False):
13.     """
14.     return: trucated image according to window center and window width
15.     and normalized to [0,1]
16.     """
17.     minWindow = float(windowCenter) - 0.5*float(windowWidth)
18.     newing = (ct_array - minWindow)/float(windowWidth)
19.     newing[newing < 0] = 0
20.     newing[newing > 1] = 1
21.     #将值域转到0-255之间,例如要看头颅时， 我们只需将头颅的值域转换到 0-255 就行了
22.     if not normal:
23.         newing = (newing *255).astype('uint8')
24.     return newing
26. ct_path = 'G:\LITS\Training_dataset'
27. saved_path = 'C:\\Users\\Cerry\\Desktop\\wl'
28. name_list = ['volume-0.nii']
29. for name in name_list:
30.     ct = sitk.ReadImage(os.path.join(ct_path,name))
31.     origin = ct.GetOrigin()
32.     direction = ct.GetDirection()
33.     xyz_thickness = ct.GetSpacing()
34.     ct_array = sitk.GetArrayFromImage(ct)
35.     seg_array = sitk.GetArrayFromImge(sitk.ReadImage(os.path.join(ct_path,name.replace('volume'
36.                                                                        ,'segmentation'))))
38. seg_bg = seg_array == 0
39. seg_liver = seg_array >= 1
40. seg_tumor = seg_array == 2
42. ct_bg = ct_array * seg_bg
43. ct_liver = ct_array*seg_liver
44. ct_tumor = ct_array*seg_tumor
46. liver_min = ct_liver.min(）
47. liver_max = ct_liver.max()
48. tumor_min = ct_tumor.min()
49. tumor_max = ct_tumor.max()
51. #by liver
52. liver_wide = liver_max - liver_min
53. liver_center = (liver_max + liver_min)/2
54. liver_wl = window_transform(ct_array,liver_wide,liver_center,normal=True)
55. saved_name = os.path.join(saved_path,'liver_w1_1.nii')
56. saved_preprocessed(liver_w1,origin,direction,xyz_thickness,saved_name)
58. #by tumor(recommended)
59. tumor_wide = tumor_max - tumor_min
60. tumor_center = (tumor_max + tumor_min) / 2
61. tumor_wl = window_transform(ct_array, tumor_wide, tumor_center, normal=True)
62. saved_name = os.path.join(saved_path, 'tumor_wl_1.nii')
63. saved_preprocessed(tumor_wl, origin, direction, xyz_thickness, saved_name)

复制代码

注意：有研究者认为如果用下面的公式来实现，实验结果会非常的糟糕！

1. #先算图像的最大和最小值
2. min = (2*window_center - window_width)/2.0 + 0.5
3. max = (2*window_center + window_width)/2.0 + 0.5
4. for i in range(nNumPixels):
5.     disp_pixel_val = (pixel_val - min)*255.0/(max - min)
8. #或者：
9. dFactor = 255.0 / (max - min)
10.     d, h, w = np.shape(img)
11.     for n in np.arange(d):
12.         for i in np.arange(h):
13.             for j in np.arange(w):
14.                 img[n, i, j] = int((img[n, i, j] - min) * dFactor)
18. min_index = img < 0
19. img[min_index] = 0
20. max_index = img > 255
21. img[max_index] = 255
23. return img

复制代码

CT图像不同扫描面的像素尺寸、粗细粒度是不同的，这对进行CNN有不好的影响，因此可能需要进行重构采样，这里可参考CT图像相关知识与SimpleITK的使用介绍的12
并且，为了更好地进行网络训练，通常采用数据归一化。
参考
CT图像之Hu值变换与窗宽窗位调整
医学影像“调窗”（window-leveling）的算法
CT图像的相关知识

来源：https://blog.csdn.net/qq_44289607/article/details/123394110
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！