学习资料：
UNet-肝脏肿瘤图像语义分割
医学图像预处理（三）——windowing(ct对比增强）
Python CT图像预处理——nii格式读取、重采样、窗宽窗位设置
医学图像处理步骤为：
1、数据加载
2、CT图像增强
3、直方图均衡化增强
4、获取对应CT、掩膜灰度图并保存
DICOM数据读取：使用pydicom库中的dcmread函数进行读取
image_slices = [pydicom.dcmread(os.path.join(data_path, file_name)) for file_name in os.listdir(data_path)]
上述代码等同于
image_slices_test = []
for file_name_test in os.listdir(data_path):
image_slices_test.append(pydicom.dcmread(os.path.join(data_path, file_name_test)))
作用是读取data_path中的文件，放入到image_slices中，得到一个129个大小的list，每一个list元素表示一个切片的全部信息。
排序：按照InstanceNumber(扫描的EPI时间序列顺序)顺序
img_slices.sort(key=lambda x:x.InstanceNumber)
提取像素值：
img_array = np.array([i.pixel_array for i in img_slices])
对于dicom、nii格式的图片，只要是ct图，就都可以选择将储存的原始数据转化为Hu值，Hu值代表了物体真正的密度。
对于nii格式的图片，在使用nibabel、simpleitk常用的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值和像素值的转化，以dicm数据为例
HU = pixel_val * slope + intercept（如果Slope为1，Intercept为0，则不需要转化的问题，具体Slope和Intercept的值在头文件中都会有说明的）
dicom转化为Hu的代码如下：
def get_pixels_hu(scans):
# np.stack()
# pixel_array 是一个正常的像素矩阵，512x512
image = np.stack([s.pixel_array for s in scans])
image = image.astype(np.int16)
image[image == -2000] = 0
# 缩放斜率 RescaleIntercept
intercept = scans[0].RescaleIntercept
# 截距 RescaleSlope
slope = scans[0].RescaleSlope
if slope != 1:
image = slope * image.astype(np.float64)
image = image.astype(np.int16)
image += np.int16(intercept)
return np.array(image, dtype=np.int16)
根据Hu（CT）值来筛选我们想要的部位的图片，其他部位全部抹黑或者抹白，目的是为了增加图像对比度。使用windowing方法，观察的CT值范围：窗宽。观察的中心CT值即为窗位，然后对肿瘤部分进行二值化处理。
def windowing(img, window_width, window_center):
# params:需要增强的图片, 窗口宽度, 窗中心 通过窗口最小值来线性移动窗口增强
min_windows = float(window_center)-0.5*float(window_width)
new_img = (img-min_windows)/float(window_width)
new_img[new_img<0] = 0 # 二值化处理 抹白
new_img[new_img>1] = 1 # 抹黑
return (new_img * 255).astype('uint8') # 把数据整理成标准图像格式
直方图均衡化
将整个图像分成许多小块，对每个小块进行均衡化。Opencv中将这种方法为：cv2.createCLAHE（）
def clahe_equalized(imgs):
# 输入imgs的形状必须是3维
assert (len(imgs.shape) == 3)
# 定义均衡化函数
clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8, 8))
# 新数组存放均衡化后的数据
img_res = np.zeros_like(imgs)
for i in range(len(imgs)):
img_res[i,:,:] = clahe.apply(np.array(imgs[i,:,:], dtype=np.uint8))
return img_res/255
nii格式CT数据读取
一、读取：使用nibabel包进行读写查看等操作
nibabel读取会将图像旋转九十度，也就是各轴的对应关系为[z, x, y]
import nibabel as nb
img = nb.load(xxx.nii.gz) #读取nii格式文件
img_affine = img.affine
data = img.get_data()
存储nii格式文件，保存时要联通affine矩阵一起保存
import nibabel as nb
nb.Nifti1Image(data,img_affine).to_filename(xxx.nii.gz)
查看
from nibabel.viewers import OrthoSlicer3D
OrthoSlicer3D(data.transpose(1,2,0)).show()
重采样
CT图像的重采样视为了使体数据中大小不同的体素变得大小相同。由于CT不同轴的体素尺寸、粗细粒度不同，直接使用不利于进行深度模型的训练和推理
重采样步骤：
1、获取先从nii文件的头文件中获取各轴单位体素对应的空间距离
header = img.header
print(header)
打印出来可见如下信息：
print(header)
<class 'nibabel.nifti1.Nifti1Header'> object, endian='<'
sizeof_hdr : 348
data_type : b''
db_name : b''
extents : 0
session_error : 0
regular : b'r'
dim_info : 0
dim : [ 3 512 512 81 1 1 1 1]
intent_p1 : 0.0
intent_p2 : 0.0
intent_p3 : 0.0
intent_code : none
datatype : int16
bitpix : 16
slice_start : 0
pixdim : [-1. 0.881 0.881 5. 0. 0. 0. 0. ]
vox_offset : 0.0
scl_slope : nan
scl_inter : nan
slice_end : 0
slice_code : unknown
xyzt_units : 10
cal_max : 0.0
cal_min : 0.0
slice_duration : 0.0
toffset : 0.0
glmax : 0
glmin : 0
descrip : b'Time=151745.150'
aux_file : b'Venous_Phase'
qform_code : scanner
sform_code : scanner
quatern_b : 0.0
quatern_c : 1.0
quatern_d : 0.0
qoffset_x : 217.3328
qoffset_y : -225.04568
qoffset_z : 1390.0
srow_x : [ -0.881 -0. 0. 217.3328]
srow_y : [ 0. 0.881 0. -225.04568]
srow_z : [ 0. 0. 5. 1390.]
intent_name : b''
magic : b'n+1'
关键字段含义如下：
sizeof_hdr：是保存文件的头文件大小，NIFTI-1或者ANALYZE格式的文件sizeof_hdr=348
dim_info：dim_info字段存储着频率编码方向(1, 2, 3)，相位编码方向(1, 2, 3)和采集期间层选择方向(1, 2, 3)，对于径向采集来讲，频率编码和相位编码都设置为0。
dim：short dim[8]保存着前面提到的图像的维度信息。如果第0维不是(1-7)之间的数字，那么这个数据具有相反的字节顺序，所以应该进行字节交换。
datatype中存储的是数据的类型。
bitpix字段必须与datatype中的代码所对应的bit(s)/pix的大小相等。
slice切片信息：
slice_start，slice_end，slice_code，slice_furation
slice_duration是存储功能磁共振成像采集的时间相关信息，需要与dim_info字段一起使用。
pixdim体素维度：每个体素维度信息都保存在pixdim中，各自对应dim，但pixdm[0]有特殊意义，其值只能是-1或1。前四个维度将在xyzt_units字段中指定。pixdim[1]对应x轴，pixdim[2]对应y轴，pixdim[3]对应z轴。
vox_offset体素偏移量：vox_offset指单个文件(.nii)图像数据的字节偏移量。
在重采样中主要需要使用到pixdim这一个字段的数据。
pixdim : [-1. 0.881 0.881 5. 0. 0. 0. 0. ]
此外，还需要知道各轴的体素总数：
width, height, channel = img.dataobj.shape
print(width, height, channel)
结果如下：x轴数据的空间大小为512 * pixdim[1] = 451.07199mm，y轴数据的空间大小为512 * pixdim[2] = 451.07199mm，z轴数据的空间大小为81 * pixdim[3] = 405mm。
512 512 81
根据各轴实际对应的空间大小，可以根据重采样后期望各体素在某个轴上单位体素对应的空间大小，算各种重采样后的体素总数。此处以各轴重采样为1mm spacing为例。使用到的是skimage中的transform函数，调用resize使重采样后的数据具有指定的shape。
### Get original space
width, height, channel = img.dataobj.shape
ori_space = [width,height,channel]
# Resample to have 1.0 spacing in all axes
spacing = [1.0, 1.0, 1.0]
data_resampled = resample(data,ori_space, header, spacing)
OrthoSlicer3D(data_resampled).show()
def resample(data,ori_space, header, spacing):
### Calculate new space
new_width = round(ori_space[0] * header['pixdim'][1] / spacing[0])
new_height = round(ori_space[1] * header['pixdim'][2] / spacing[1])
new_channel = round(ori_space[2] * header['pixdim'][3] / spacing[2])
new_space = [new_width, new_height, new_channel]
data_resampled = skimage.transform.resize(data,new_space,order=1, mode='reflect', cval=0, clip=True, preserve_range=False, anti_aliasing=True, anti_aliasing_sigma=None)
return data_resampled
医学图像领域的窗口技术，包括窗宽(window width) 窗位(window center)，用于选择感兴趣的CT值范围。窗宽和窗位分别对应图像的对比度与亮度。
note：窗宽窗位是CT图像特有的概念，MRI中并没有；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为妥。
窗宽窗位的设置方法：
一、手动设置
根据常用的预设值设置窗宽窗位。实例中使用的为腹部ct，常用的预设值为窗宽 300-500，窗位 30-50。以下示例中使用窗宽400，窗位40。
w_width = 400
w_center = 40
data_adjusted1 = adjustMethod1(data_resampled,w_width,w_center)
def adjustMethod1(data_resampled,w_width,w_center):
val_min = w_center - (w_width / 2)
val_max = w_center + (w_width / 2)
data_adjusted = data_resampled.copy()
data_adjusted[data_resampled < val_min] = val_min
data_adjusted[data_resampled > val_max] = val_max
return data_adjusted
方法二：
根据观察目标的mask，提取对应的目标区域；统计其最大值和最小值，从而计算窗宽窗位。
重采样和窗宽窗位调整参考以下几篇文章：
[https://blog.csdn.net/qq_46511579/article/details/118441519](https://blog.csdn.net/qq_46511579/article/details/118441519)
[https://blog.csdn.net/qq_44289607/article/details/123495983](https://blog.csdn.net/qq_44289607/article/details/123495983)
[https://blog.csdn.net/qq_44289607/article/details/123394110?spm=1001.2014.3001.5502](https://blog.csdn.net/qq_44289607/article/details/123394110?spm=1001.2014.3001.5502)
后续学习：
参考文章：
【医学影像数据处理】nii、npz、npy、dcm、mhd 的数据格式互转，及多目标分割处理汇总
.npy、.mha、.npz
扫一扫
1184
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
