医工互联

 找回密码
 注册[Register]

手机动态码快速登录

手机号快速登录

微信登录

微信扫一扫,快速登录

QQ登录

只需一步,快速开始

查看: 219|回复: 0
收起左侧

Python CT图像预处理——nii格式读取、重采样、窗宽窗位设置

[复制链接]

  离线 

发表于 2023-2-23 22:15:25 来自手机 | 显示全部楼层 |阅读模式 <
nii格式CT数据读取

遇到nii格式的CT数据,可以通过nibabel包进行数据的读、写、查看等操作。下面列出常见操作。
读写nii格式文件

nibabel读取数据会将图像旋转九十度,也就是各轴的对应关系为[z,x,y]
  1. import nibabel as  nb
  2. img = nb.load(xxx.nii.gz) #读取nii格式文件
  3. img_affine = img.affine
  4. data = img.get_data()
复制代码
存储nii格式文件,保存时要联通affine矩阵一起保存
  1. import nibabel as  nb
  2. nb.Nifti1Image(data,img_affine).to_filename(xxx.nii.gz)
复制代码
查看

  1. from nibabel.viewers import OrthoSlicer3D
  2. OrthoSlicer3D(data.transpose(1,2,0)).show()
复制代码
重采样

CT图像的重采样是为了使体数据中大小不同的体素变得大小相同。由于CT不同轴的体素尺寸、粗细粒度不同,直接使用不利于进行深度模型的训练和推理。
首先我们可以获取先从nii文件的头文件中获取各轴单位体素对应的空间距离。
  1. header = img.header
  2. print(header)
复制代码
打印出来可见如下信息:
  1. print(header)
  2. <class 'nibabel.nifti1.Nifti1Header'> object, endian='<'
  3. sizeof_hdr      : 348
  4. data_type       : b''
  5. db_name         : b''
  6. extents         : 0
  7. session_error   : 0
  8. regular         : b'r'
  9. dim_info        : 0
  10. dim             : [  3 512 512  81   1   1   1   1]
  11. intent_p1       : 0.0
  12. intent_p2       : 0.0
  13. intent_p3       : 0.0
  14. intent_code     : none
  15. datatype        : int16
  16. bitpix          : 16
  17. slice_start     : 0
  18. pixdim          : [-1.     0.881  0.881  5.     0.     0.     0.     0.   ]
  19. vox_offset      : 0.0
  20. scl_slope       : nan
  21. scl_inter       : nan
  22. slice_end       : 0
  23. slice_code      : unknown
  24. xyzt_units      : 10
  25. cal_max         : 0.0
  26. cal_min         : 0.0
  27. slice_duration  : 0.0
  28. toffset         : 0.0
  29. glmax           : 0
  30. glmin           : 0
  31. descrip         : b'Time=151745.150'
  32. aux_file        : b'Venous_Phase'
  33. qform_code      : scanner
  34. sform_code      : scanner
  35. quatern_b       : 0.0
  36. quatern_c       : 1.0
  37. quatern_d       : 0.0
  38. qoffset_x       : 217.3328
  39. qoffset_y       : -225.04568
  40. qoffset_z       : 1390.0
  41. srow_x          : [ -0.881   -0.       0.     217.3328]
  42. srow_y          : [   0.         0.881      0.      -225.04568]
  43. srow_z          : [   0.    0.    5. 1390.]
  44. intent_name     : b''
  45. magic           : b'n+1'
复制代码
关键字段含义如下:
sizeof_hdr 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)之间的数字,那么这个数据具有相反的字节顺序,所以应该进行字节交换(NIFTI标准没有提供字节顺序的字段,提倡使用dim[0])。
datatype中存储的是数据的类型。
bitpix字段必须与datatype中的代码所对应的bit(s)/pix的大小相等。
slice切片信息
包含字段:slice_start,slice_end, slice_code, slice_duration
slice_duration是存储功能磁共振成像采集的时间相关信息,需要与dim_info字段一起使用。
pixdim体素维度:每个体素维度信息都保存在pixdim中,各自对应dim,但pixdim[0]有特殊意义,其值只能是-1或1。前四个维度将在xyzt_units字段中指定。pixdim[1]对应x轴,pixdim[2]对应y轴,pixdim[3]对应z轴
vox_offset体素偏移量:vox_offset指 单个文件(.nii)图像数据的字节偏移量。
在重采样中主要需要使用到pixdim这一个字段的数据。
  1. pixdim          : [-1.     0.881  0.881  5.     0.     0.     0.     0.   ]
复制代码
此外,我们还需要知道各轴的体素总数:
  1. width, height, channel = img.dataobj.shape
  2. print(width, height, channel)
复制代码
结果为:
  1. 512 512 81
复制代码
即x轴数据的空间大小为512*pixdim[1]=451.07199mm,y轴数据的空间大小为512*pixdim[2]=451.07199mm,z轴数据的空间大小为81*pixdim[3]=405mm。
获得了各轴实际对应的空间大小后,就可以根据重采样后期望各体素在某个轴上单位体素对应的空间大小,算各轴重采样后的体素总数。
此处以各轴重采样为1mm spacing为例。这里使用到的是skimage的transform函数,调用resize使重采样后的数据具有指定的shape。
  1. ### Get original space
  2. width, height, channel = img.dataobj.shape
  3. ori_space = [width,height,channel]
  4. # Resample to have 1.0 spacing in all axes
  5. spacing = [1.0, 1.0, 1.0]
  6. data_resampled = resample(data,ori_space, header, spacing)
  7. OrthoSlicer3D(data_resampled).show()
  8. def resample(data,ori_space, header, spacing):
  9.     ### Calculate new space
  10.     new_width = round(ori_space[0] * header['pixdim'][1] / spacing[0])
  11.     new_height = round(ori_space[1] * header['pixdim'][2] / spacing[1])
  12.     new_channel = round(ori_space[2] * header['pixdim'][3] / spacing[2])
  13.     new_space = [new_width, new_height, new_channel]
  14.     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)
  15.     return data_resampled
复制代码
窗宽窗位设置

以下文字内容来源于文章:https://blog.csdn.net/qq_44289607/article/details/123394110?spm=1001.2014.3001.5502
CT值属于医学领域的概念,通常称亨氏单位(hounsfield unit,HU),反映了组织对X射线的吸收程度。黑影表示低吸收区,即低密度区,如含气体多的肺部;白影表示高吸收区,即高密度区,如骨骼。
对于CT图像来说,它的HU值从负几千到正几千都有,直接显示图像的亮度和对比度不能充分突出关键部分。这里所指的“关键部分”在 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为妥。
  未设置时,图像显示如下图
1.png

以下为两种窗框窗位的设置方式:
方法一:手动设置窗宽窗位

根据常用的预设值设置窗宽窗位。实例中使用的为腹部ct,常用的预设值为窗宽 300-500,窗位 30-50。以下示例中使用窗宽400,窗位40。
  1. w_width = 400
  2. w_center = 40
  3. data_adjusted1 = adjustMethod1(data_resampled,w_width,w_center)
  4. def adjustMethod1(data_resampled,w_width,w_center):
  5.     val_min = w_center - (w_width / 2)
  6.     val_max = w_center + (w_width / 2)
  7.     data_adjusted = data_resampled.copy()
  8.     data_adjusted[data_resampled < val_min] = val_min
  9.     data_adjusted[data_resampled > val_max] = val_max
  10.     return data_adjusted
复制代码
设置窗宽窗位之后,对比度与亮度有了明显的改善。
2.png

方法二:

如果存在目标组织的mask,我们可以将对应的组织提取出来,统计相应的最大值和最小值,并用其计算窗宽窗位。
当然该方法是在mask存在的情况下才能实现。
优点:每个case的窗宽窗位依据案例量身定做。相较于整个数据集用统一的窗宽窗位,或者采用cut off,有无可比拟的优点。
重采样和窗宽窗位调整参考以下几篇文章:
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/123394110?spm=1001.2014.3001.5502

来源:https://blog.csdn.net/lavinia_chen007/article/details/125389503
免责声明:如果侵犯了您的权益,请联系站长,我们会及时删除侵权内容,谢谢合作!
回复

使用道具 举报

提醒:禁止复制他人回复等『恶意灌水』行为,违者重罚!
您需要登录后才可以回帖 登录 | 注册[Register] 手机动态码快速登录 微信登录

本版积分规则

发布主题 快速回复 收藏帖子 返回列表 客服中心 搜索
简体中文 繁體中文 English 한국 사람 日本語 Deutsch русский بالعربية TÜRKÇE português คนไทย french

QQ|RSS订阅|小黑屋|处罚记录|手机版|联系我们|Archiver|医工互联 |粤ICP备2021178090号 |网站地图

GMT+8, 2024-11-5 12:30 , Processed in 0.282384 second(s), 64 queries .

Powered by Discuz!

Copyright © 2001-2023, Discuz! Team.

快速回复 返回顶部 返回列表