KITS+肾脏肿瘤预处理+重采样+窗体变换+强度裁剪

摘要

目前正在做肾脏肿瘤的预处理问题，比赛地址。
KiTS19 年比赛的结果。
JunMa大佬在自己的知乎里面给出了自己的看法和对于论文的学习与展望。

我最近正在解决数据预处理的问题，论文的数据预处理分为了重采样，spacing 到 3.22mmx1.62mmx1.62mm，然后将灰度值cut off到[-79,304],然后将cut off 之后的数据进行z-score（subtract 101 and divide by 76.9 to bring the intensity values in a range that is
more easily processed by CNNs）。最后将输入神经网络的数据大小定为80x160x160。

重采样

首先强调一下，我对于论文的复现来自于我对于论文的解读以及自己代码水平，不一定对，但是提供一个思路
我一开始对于spacing到3.22mmx1.62mmx1.62mm表示十分的不理解，像素之间的距离难道不都是1吗？像素距离和实际测量的毫米有什么关系，后来我才知道，CT图像里面是保存有图像的pixdim信息的。具体代码如下

 #img_name 就是一个nii.gz文件的地址img = nibabel.load(img_name)pix_dim = (img.header.structarr['pixdim'][1],img.header.structarr['pixdim'][2],img.header.structarr['pixdim'][3])print(pix_dim)

这样就可以读取pixdim数据我的第一个图片的数据是(0.5mm, 0.9199219mm,0.9199219mm)也就是说它的实际spacing是这个。使用skimage里面的transform.rescale()就可以把数据处理成自己想要得spacing，具体写法如下（对于image的操作）

 from skimage import transformtarget_resolution = (3.22,1.62,1.62)img = np.array(img.get_data())scale_vector = (pix_dim[0]/target_resolution[0],pix_dim[1]/target_resolution[1],pix_dim[2]/target_resolution[2])"""order = 1 样条插值的顺序，必须是在0-5范围内，默认是1，详情参见skimage.transform.warp 这个很重要mode = {'constant','edge','symmetric','reflect','warp'}可选输入边界之外的点根据到给定的模式，模式匹配行为'numpy.pad'cval 与模式constant结合使用，外部值图像边界剪辑 : bool 可选， 是否将输出剪辑到输入值范围，这是默认启用的，因为高阶线性插值可能产生超出给定输入值的范围。preserve_range : bool 可选  是否保持原始值的范围。否则输入图像根据img_as_float的约定进行转换 详情： https://scikit-image.org/docs/dev/user_guide/data_types.htmlmultichannel : bool 可选 图像的最后一个轴是否被解释为多个通道或者空间维度。anti_aliasing : bool 可选 是否用高斯滤波器来平滑图像，缩小比例。将图像下采样到避免解释为多个渠道或者其他空间维度避免混叠伪影anti_aliasing_sigma : {float, tuple of floats} optional 高斯滤波的标准差，以避免混叠伪影，默认情况下，此值选择为(1-2)/2，其中s是缩小比例因子"""yinzi = transform.rescale(img, scale_vector,order=1, mode='constant', cval=0, clip=True,preserve_range=True, multichannel=False,anti_aliasing=True, anti_aliasing_sigma=None)

scale_vector 是 pixdim和target_resolution的比值，对于这个rescale真的有很多参数，如果不想改变原来的值就把order=0，并且把anti_aliasing=False ，并且把preserve_range=True(这个真的特别重要), 因为如果不这样的话，mask里面的数据只有[0,1,2],rescale之后数据就会变得很乱很乱，对于mask的操作如下

 img = nibabel.load(img_name)print(img.shape)pix_dim = (img.header.structarr['pixdim'][1],img.header.structarr['pixdim'][2],img.header.structarr['pixdim'][3])print('pix_dim',pix_dim)img = np.array(img.get_data())print('non_transpos_shape',img.shape)print(np.unique(img))target_resolution = (3.22,1.62,1.62)scale_vector = (pix_dim[0] / target_resolution[0],pix_dim[1] / target_resolution[1],pix_dim[2] / target_resolution[2])img = img.astype(np.float)print('img.dtype', img.dtype)print(np.unique(img))# uint8# 这里把uint8转换成了float64了yinzi = transform.rescale(img, scale_vector,order=0, mode='constant', cval=0, clip=True,preserve_range=True, multichannel=False,anti_aliasing=False, anti_aliasing_sigma=None)

这里比较建议大家去读一读源码里面的注释，然后参考输入输出类型的解释的官方页面进行学习。

窗体变换+强度裁剪

强度裁剪

ct_array[ct_array < -79.0] = -79.0
ct_array[ct_array > 304.0] = 304.0

z-score

newimg = (ct_array - 101) / float(76.9)

crop成（80,160,160）

这个真的是自己瞎编了，我用最传统的赋值方法做的，具体对不对我也不知道了
代码：

def z_crop(img_array,crop_array,z,xy):if z >= 80 and xy >= 160:side = z - 80startz = side//2endnum = side - startzendz = z - endnumstartxy = (xy-160)//2endxy = xy - ((xy-160)-startxy)crop_array[:,:,:] = img_array[startz:endz,startxy:endxy,startxy:endxy]elif z>=80 and xy < 160:side = z - 80startz = side//2endnum = side - startzendz = z - endnumstartxy = (160 - xy) // 2endxy = 160 - ((160 - xy) - startxy)crop_array[:,startxy:endxy,startxy:endxy] = img_array[startz:endz,:,:]elif z < 80 and xy >= 160:side =  80 - zstartz = side//2endnum = side - startzendz = 80 - endnumstartxy = (xy-160)//2endxy = xy - ((xy-160)-startxy)crop_array[startz:endz,:,:] = img_array[:,startxy:endxy,startxy:endxy]elif z < 80 and xy < 160:side =  80 - zstartz = side//2endnum = side - startzendz = 80 - endnumprint(side)print(startz)print(endnum)print(endz)startxy = (160 - xy) // 2endxy = 160 - ((160 - xy) - startxy)crop_array[startz:endz,startxy:endxy,startxy:endxy] = img_array[:,:,:]return crop_array

最后转换成np文件

def preprocessingtonp(pathImg_list):for subsetindex in range(len(pathImg_list)):try:print(pathImg_list[subsetindex])subset_path = KiTs_path + "/" + str(pathImg_list[subsetindex]) + deal_name # 改成计算groundtruthsubmask_path = KiTs_path + "/" + str(pathImg_list[subsetindex]) + dealseg_name#这个地方可以随便写的，如果写斜杠或者反斜杠都可以#print(subset_path)# cropimg_src = sitk.ReadImage(subset_path,sitk.sitkFloat64)img_array = sitk.GetArrayFromImage(img_src)mask_src = sitk.ReadImage(submask_path,sitk.sitkInt8)mask_array = sitk.GetArrayFromImage(mask_src)print(img_array.shape)print(mask_array.shape)print(type(img_array))print(type(mask_array))print(img_array.max())print(img_array.min())print(mask_array.max())print(mask_array.min())file_path = outputImg_path + "/" + str(pathImg_list[subsetindex]) + ".npy"mask_path = outputMask_path + "/" + str(pathImg_list[subsetindex]) + ".npy"## np.save(file_path, img_array)# np.save(mask_path,mask_array)except:subset_path = KiTs_path + "/" + str(pathImg_list[subsetindex]) + image_nameprint(subset_path,"不行")

总结一下

就这不到500行代码写了将近10天，我感觉比较难的地方第一个是重采样部分，要是不是学长给了我一段代码我现在都不知道如何resampling，第二个是transform.rescale，我一开始没有想过后面那些参数这么重要，我一开始就只传入了一个image，一个scale因子的范围，没想到处理mask 一个 order=0 一个mode=“constant” 一个preserve_range=True 一个anti_aliasing=False 一套连招才能处理好mask，读源码的重要性不言而喻。
最后crop我也是瞎几把crop的，我没有根据肾脏的部位来，直接crop中间，我看了看有些数据集还是切到了mask，我也不知道该怎么办了，去他妈的吧，先去训练了。
窗体裁剪也是我自己随着性子写的，网上给了一套方法，我一开始也是按照他的想法写的,虽然总感觉不对劲，但是这篇博客写的基础知识还是很全面的，建议大家可以看一下。

参考资料

https://kits19.grand-challenge.org/
http://results.kits-challenge.org/miccai2019/
https://zhuanlan.zhihu.com/p/76057976
https://scikit-image.org/docs/dev/user_guide/data_types.html
https://blog.csdn.net/normol/article/details/88313888
https://blog.csdn.net/qq_34532604/article/details/104206045
这篇讲重采样不错