nii格式CT数据读取

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

读写nii格式文件

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不同轴的体素尺寸、粗细粒度不同，直接使用不利于进行深度模型的训练和推理。

首先我们可以获取先从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 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这一个字段的数据。

pixdim          : [-1.     0.881  0.881  5.     0.     0.     0.     0.   ]

此外，我们还需要知道各轴的体素总数：

width, height, channel = img.dataobj.shape
print(width, height, channel)

结果为：

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。

### 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

窗宽窗位设置

以下文字内容来源于文章：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为妥。

未设置时，图像显示如下图

以下为两种窗框窗位的设置方式：

方法一：手动设置窗宽窗位

根据常用的预设值设置窗宽窗位。实例中使用的为腹部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

设置窗宽窗位之后，对比度与亮度有了明显的改善。
Python CT图像预处理——nii格式读取、重采样、窗宽窗位设置