【医学图像预处理过程】

HUI 别摸鱼了

已于 2024-10-08 08:36:29 修改

阅读量4k

点赞数 17

分类专栏：论文笔记文章标签： python

于 2023-06-29 02:32:14 首次发布

本文链接：https://blog.csdn.net/weixin_47244593/article/details/131421465

版权

论文笔记专栏收录该内容

5 篇文章

订阅专栏

注：仅为自用，想到写到哪儿，可能有点乱。

CT图像的预处理操作步骤

step1: 数据加载
step2: 重采样(resampling)
step3: window transform
step4: get mask effective range
step5: generate subimage（或resize图像尺寸）
- 1. 原因
- 2. 方法-裁剪&填充
step6：save the preprocessed data
参考文献

step1: 数据加载

1. 先导知识

1）数字存储

图像宽高： $W, H$ ；
量化级别： $L=2^{k}$ ；
数字存储比特级： $b = W * H * k$
$e g :$
若 $W, H = 384 * 384$ ，256灰度值，则
$b = 384 * 384 * 8$

2）灰度直方图

灰度级的函数；
具有该灰度级的像素个数。
$h(r_k)=n_k$
——灰度级 $r_k$ ；
——该灰度级的像素个数 $n_k$ 。

绘制直方图：直方图与直方图均衡化、图像的直方图

3）Image data types

参考：Image data types

python 读入图像形式为unit8；
PIL读取为unit8； $2^{8}=256$
unit8范围：0~255；
float范围：-1 to 1 or 0 to 1；
图像矩阵形式：0~1。

4）图像读取的宽高（W,H）顺序

在这里插入图片描述

OpenCV的坐标系原点(0,0)是图片的左上角点，这种坐标系在OpenCV的结构体Mat，Rect，Point中都是适用。
OpenCV中坐标系的X轴为图像矩形的上水平线，从左往右；Y轴为图像矩形的左垂直线，从上到下。
注意： $P o in t (x, y)$ 和 $R ec t (x, y)$ 中，第一个参数 $x$ 代表的是元素所在图像的列数cols; 第二个参数 $y$ 代表的是元素所在图像的行数rows。
$x$ 代表width（宽），即column； $y$ 代表height（高），即row。

行row == height == Point.y
列col == width  == Point.x

python中对opencv像素操作，可以认为是对Numpy数组进行操作。

图片原生的表示是(width,height)，这也是电脑编辑图片能看到的分辨率格式。
但由于图片的储存形式为矩阵，而对于矩阵，我们熟悉的都是row first（一般的数学概念上都会是row first），所以说在涉及到具体操作时，还是将图片转化为(height, width)比较方便，因为index时[ $i$ -th row]，[ $j$ -th column]也就对应着(height, width)。
这也是下面例子中，原生格式为(w,h)转化为numpy数组后，就自动变成了(h,w)。

import cv2
from PIL import Image
import numpy as np
from skimage.io import  imread

filepath = './liver.jpg'

cv2_im = cv2.imread(filepath)
print('cv2_im shape ',cv2_im.shape) # (height, width, ch)

im = Image.open(filepath)
print('PIL image size', im.size) # (width, height, ch)
pil_im = np.asarray(im)
print('pil_im shape ',pil_im.shape) # (height, width, ch)

sk_im = imread(filepath)
print('sk_im shape', sk_im.shape) # (height, width, ch)

# pydicom: dcm_slice = pydicom.dcmread(filename)
#          slice_data = dcm_slice.pixel_array   shape:(height, width)
#          若想更改原始像素size或其他操作，必须写回至原始二进制属性PixelData 中
#          dcm_slice.PixelData = dcm_slice.pixel_array.tobytes()
# SimpleITK: itk_im = sitk.ReadImage()  size:(width, height, depth)
#            np_im = sitk.GetArrayFromImage(itk_im) shape:(depth, height, width)

5）人体组织CT范围

人体组织CT范围

MICCAI 2019肾&肿瘤分割挑战赛第一名方案学习笔记

MICCAI BraTS 2020脑肿瘤分割挑战赛冠军方法学习笔记

2. 医学图像格式（nii和dicom）

1）定义

放射生物图像中主要有六种格式–DICOM（医疗中的数字图像和通信），NIFTI（神经影像学信息技术计划），PAR/REC（飞利浦 MRI 扫描格式），ANALYZE（Mayo 医疗成像）以及 NRRD（近乎原始光栅数据）和 MNIC 格式。
其中 DICOM 和 NIFTI 是最常用的格式。
在这里插入图片描述

DICOM格式（全称Digital Imaging and Communications in Medicine，医疗数字影像传输协定），与PNG，JPG，bmp类似，是一种电子文件的储存格式。需要有电脑里专门的软件打开或编辑。
NIfTI 格式的主要特点就是它包含两个能够将每个体素的索引（i,j,k）和它的空间位置（x,y,z）关联起来的仿射坐标。
DICOM 和 NIfTI 这两种格式的主要区别是：NIfTI 中的图像原始数据被存储成了 3 维图像，而 dicom 一些 2 维的图层。

2）医疗数据的组成

医疗数据有四个关键的组成部分–像素深度、光度解释、元数据以及像素数据。这几部分决定了图像的大小和分辨率。

像素深度（Pixel Depth）或者位深度（Bit Depth）或者色深度（Color Depth）就是用来编码每一像素的信息所用的位数。例如，一个 8 位的栅格会拥有从 0 到 255 这 256 种各不相同的数值。
参考：医学图像处理(一)：医学图像格式（nii和dicom）

3. DICOM和nii数据读取

1） DICOM数据读取

用到的包：pydicom；
读取数据：read_file和dcmread两种方式，两种方式读取效果相同。

# 采用pydicom模块
import pydicom

# 数据路径
file_path = r"C:\Users\os\Desktop\data\xxxxx.dcm"

# read_file
data0 = pydicom.read_file(file_path)
# file_data = data0.pixel_array
# print(file_data)
# dcmread
data1 = pydicom.dcmread(file_path)

eg：读取dicom文件

# 采用pydicom模块
import pydicom

label_path = '~/Downloads/3Dircadb/3Dircadb1.1/LABELLED_DICOM'
# 若不遍历读出，则获得的是单张切片，而非三维图片
slices = [pydicom.dcmread(label_path + '/' + s) for s in os.listdir(label_path)]
slices.sort(key = lambda x: int(x.InstanceNumber))
# slice[i]为pydicom.dataset.FileDataset，储存了切片相关的信息
# .pixel_array用来获取像素信息,每张切片都是灰度图
#2019-6-5更新：s.pixel_array.shape:(height,width)
image = np.stack([s.pixel_array for s in slices])
image = image.astype(np.int16)
# 目前还不知道这些label值对应的器官分别是什么？
np.unique(image)
#array([  0,   1,  17,  33,  65,  97, 129, 193, 257, 321, 385, 449, 451,
#       453, 465, 481, 513], dtype=int16)
label_path = '~/Downloads/3Dircadb/3Dircadb1.1/MASKS_DICOM/artery'
slices = [pydicom.dcmread(label_path + '/' + s) for s in os.listdir(label_path)]
slices.sort(key = lambda x: int(x.InstanceNumber))
image = np.stack([s.pixel_array for s in slices])
image = image.astype(np.int16)
np.unique(image)
# array([  0, 255], dtype=int16)

2）nii数据读取

用到的包：SimpleITK
seg = itk.ReadImage()读取出的顺序为(x,y,z)，即(width, height, depth)，可用seg.GetSize() 、seg.GetWidth()、 seg.GetHeight()、 image.GetDepth()来进行验证。
注意！！当将其转化为数组后segimg = sitk.GetArrayFromImage(seg)，通道顺序会反向，变成(z,y,x) ，即(depth, height, width),之所以会变成(z,y,x)，因为python中的numpy默认的order是C，即变换最快的轴向在最后，如改成内部实现改成order=Fortran，则为(x,y,z)，matlab中默认都是order=Fortran。
补充：
读取后立即检查 CT 和分割数组的尺寸。
除非明确需要，否则删除不必要的旋转和变换。
确保以正确的方向写入切片（通常情况下，除非需要，否则不进行轴交换）。

import SimpleITK as sitk

img = sitk.ReadImage(ct_path, sitk.sitkInt16) # 可自行改变存储类型
print(img.GetSize()) # (x,y,z)
img_array =sitk.GetArrayFromImage(img)
print(img_array.shape) # (z,y,x)
# 保存
sitk.WriteImage(out,'simpleitk_save.nii.gz')

将 numpy array 保存为 nii 格式图像之后，有时候会发现使用 itk-snap 却打不开，这是为什么呢？
因为 itk-snap 只接受 int16 类型的数据，所以你需要将 numpy 数组先强制转换成 int16 类型。

array = array.astype(np.int16)
out = sitk.GetImageFromArray(array)
sitk.WriteImage(out,'save.nii.gz')

NIfTI格式医学图像不同方向之间旋转

附： Python的Numpy库中有不少函数都有order参数，比如函数numpy.array()、函数numpy.reshape()。order参数用于控制数据的存储结构(Memory Layout)。

常用的数据存储结构(Memory Layout)有两种：

一是C语言的数据存储结构，即C odrder，C odrder的本质是存储时行优先，即一行一行的存储。C是最后一个索引变化最快。order='C’即按行读取数据并存储到内存中。
二是FORTRAN语言的数据存储结构，即 F order，F order的本质是存储时列优先，即一列一列的存储。F是第一个索引变化最快。order='F’即按列读取数据并存储到内存中。

参考：

eg：公共数据集LITS17

为了方便处理，对其进行了分类处理：

unzip ./*.zip
mkdir ct 
mkdir seg
mv ./volume-*.nii ./ct/
mv ./segmentation-*.nii ./seg/

这样将CT与分割图像分开来，也就有了下面示例代码的路径写法

ct_dir = '～/LITS17/ct/'
seg_dir = '～/LITS17/seg/'
for ct_file in os.listdir(ct_dir):

    # 将CT和金标准入读内存
    ct = sitk.ReadImage(os.path.join(ct_dir, ct_file), sitk.sitkInt16)
    # ct_array:(629, 512, 512)
    # 注意读取出来是z y x，即切片数量在最前
    # 而 origin和position读取出来的是 x y z
    ct_array = sitk.GetArrayFromImage(ct)
    # vol_values=np.unique(ct_array) 有2708个值

    seg = sitk.ReadImage(os.path.join(seg_dir, ct_file.replace('volume', 'segmentation')), sitk.sitkInt8)
    seg_array = sitk.GetArrayFromImage(seg)

3）dcm和nii格式的互相转化

对应的python代码可以转化；
直接将.dcm重命名为.nii后仍然可以读取；
.nii重命名为.nii.gz后仍然可以读取（也可通过gzip压缩后得到），但是反过来.nii.gz重命名为.nii后不可读取。
nii与nii.gz格式的关系

参考：

4. 普通图片（如PNG）读取

用到的包： PIL，cv2, skimage
读取出的格式区别

包	图片通道	格式
PIL	RGB	Image	channel_last
cv2	BGR	array	channel_last
skimage	RGB	array	channel_last

注意：对于opencv来说，无论读取灰度图还是彩图都是(H,W,3)的shape，灰度图的读取会把单通道复制三遍。因此，读取灰度图时
img = cv2.imread('gray.jpg', cv2.CV_LOAD_IMAGE_GRAYSCALE)或者
img = cv2.imread('gray.jpg',0)，其中0代表灰度图，1代表彩图，读出来的shape为（H,W)。

1）skimage读取图片

eg：数据集为超声波灰度图像


import os
import numpy as np

from skimage.io import imsave, imread

data_path = '～/kaggle_Ultrasound'

image_rows = 420
image_cols = 580

train_data_path = os.path.join(data_path, 'train')
images = os.listdir(train_data_path)
# 原图为灰度图
imgs = np.ndarray((total, image_rows, image_cols), dtype=np.uint8)
for image_name in images:
	img = imread(os.path.join(train_data_path, image_name), as_grey=True)
	# 相当于扩维处理 由(420,580) --> (1,420,580)
	img = np.array([img])
	imgs[i] = img

# 可选项
np.save('imgs_train.npy', imgs)
"""
当然，在输入网络前还得进行处理，包括增加通道使格式与网络匹配、归一化处理等等。
"""

2）PIL读取图片（.tif格式）

eg：公共数据集DRIVE，分割眼球血管

import os
import h5py
import numpy as np
from PIL import Image



def write_hdf5(arr,outfile):
  with h5py.File(outfile,"w") as f:
    f.create_dataset("image", data=arr, dtype=arr.dtype)


#------------Path of the images --------------------------------------------------------------
#train
original_imgs_train = "./DRIVE/training/images/"
groundTruth_imgs_train = "./DRIVE/training/1st_manual/"
borderMasks_imgs_train = "./DRIVE/training/mask/"
#test
original_imgs_test = "./DRIVE/test/images/"
groundTruth_imgs_test = "./DRIVE/test/1st_manual/"
borderMasks_imgs_test = "./DRIVE/test/mask/"
#---------------------------------------------------------------------------------------------

Nimgs = 20
channels = 3
height = 584
width = 565
dataset_path = "./DRIVE_datasets_training_testing/"

def get_datasets(imgs_dir,groundTruth_dir,borderMasks_dir,train_test="null"):
    imgs = np.empty((Nimgs,height,width,channels))
    groundTruth = np.empty((Nimgs,height,width))
    border_masks = np.empty((Nimgs,height,width))
    for path, subdirs, files in os.walk(imgs_dir): #list all files, directories in the path
        for i in range(len(files)):
            #original
            print("original image: ",files[i])
            img = Image.open(imgs_dir+files[i])
            imgs[i] = np.asarray(img)
            #corresponding ground truth
            groundTruth_name = files[i][0:2] + "_manual1.gif"
            print("ground truth name: ", groundTruth_name)
            g_truth = Image.open(groundTruth_dir + groundTruth_name)
            groundTruth[i] = np.asarray(g_truth)
            #corresponding border masks
            border_masks_name = ""
            if train_test=="train":
                border_masks_name = files[i][0:2] + "_training_mask.gif"
            elif train_test=="test":
                border_masks_name = files[i][0:2] + "_test_mask.gif"
            else:
                print("specify if train or test!!")
                exit()
            print("border masks name: ", border_masks_name)
            b_mask = Image.open(borderMasks_dir + border_masks_name)
            border_masks[i] = np.asarray(b_mask)

    print("imgs max: ", str(np.max(imgs)))
    print("imgs min: ", str(np.min(imgs)))
    assert(np.max(groundTruth)==255 and np.max(border_masks)==255)
    assert(np.min(groundTruth)==0 and np.min(border_masks)==0)
    print("ground truth and border masks are correctly withih pixel value range 0-255 (black-white)")
    assert(imgs.shape == (Nimgs,height,width,channels))
    groundTruth = np.reshape(groundTruth,(Nimgs,height,width,1))
    border_masks = np.reshape(border_masks,(Nimgs,height,width,1))
    assert(groundTruth.shape == (Nimgs,height,width,1))
    assert(border_masks.shape == (Nimgs,height,width,1))
    return imgs, groundTruth, border_masks

if not os.path.exists(dataset_path):
    os.makedirs(dataset_path)
#getting the training datasets
imgs_train, groundTruth_train, border_masks_train = get_datasets(original_imgs_train,groundTruth_imgs_train,borderMasks_imgs_train,"train")
print("saving train datasets")
write_hdf5(imgs_train, dataset_path + "DRIVE_dataset_imgs_train.hdf5")
write_hdf5(groundTruth_train, dataset_path + "DRIVE_dataset_groundTruth_train.hdf5")
write_hdf5(border_masks_train,dataset_path + "DRIVE_dataset_borderMasks_train.hdf5")

5. 模型数据读取格式npy、png、nii

1）说明

当使用PNG格式作为模型的输入时，需要先将PNG图像读入内存并转换为模型所需的数据类型和形状，这可能会占用较大的内存空间。同时，PNG图像格式的压缩率比较高，读取速度相对较慢。因此，使用PNG格式作为模型输入可能会占用较多的内存空间和时间。
当使用Numpy格式（.npy）作为模型的输入时，可以直接将Numpy数组读入内存，并将其作为模型输入传递。由于Numpy数组已经被预处理为模型所需的数据类型和形状，因此内存占用量相对较小。此外，由于Numpy数组是二进制格式，读取速度也会更快。

综上所述，使用PNG格式作为模型输入可能会占用较多的内存空间和时间，而使用Numpy格式作为模型输入则相对更加高效和方便。

2） png、npy优缺点

保存为常见的jpg格式图像（uint8：80-255）时则会损失数据精度；
- 存储空间方面： cv2.imwrite：得到数据经过压缩，因此占用空间少，同一图像得到的png大小仅4KB；
- 读取速度方面： cv2.imwrite：后续加载数据较慢，同一图像加载用时约0.017s。
保存为.npy文件时则文件大小过大。
- 存储空间方面： np.save：得到的数据占用空间较大，同一图像得到的npy大小为2MB；
- 读取速度方面： np.save：后续加载数据较快，同一图像加载用时约0.003s。

3） png转化npy

在进行深度学习处理时，有些的代码处理的数据格式为npy，但是常常有的数据格式为png，因此下面就用Python实现图像批量png格式转为npy格式，代码进行了详细注释。
参考：png，jpg，bmp等格式图像转.npy文件（附代码）、png转npy

# -- coding: UTF-8 --
'''
Author: Clouds rising
Date: December, 2023
QQ: 504156006
'''

import os
import glob
import numpy as np
from PIL import Image

# 要修改的文件夹路径
input_dir = 'label_png'

# 新文件夹的路径
output_dir = 'label_npy'

# 创建输出文件夹
if not os.path.exists(output_dir):
    os.makedirs(output_dir)

# 遍历文件夹中的所有png图像
for img_path in glob.glob(os.path.join(input_dir, '*.png')):

    # 打开原始的png图像
    img = Image.open(img_path)

    # 转换图像
    img_8bit = img.convert('L') #参数L: 8位像素，黑白

    # 创建一个转换表，将灰度值为255的像素点变为0
    '''
    在这一行代码中，创建了一个包含256个元素的列表（0到255的灰度值）。该列表的每个元素都通过条件判断
    0 if i == 255 else i 进行赋值。这意味着如果灰度值为255，那么赋值为0，否则保持不变。
    这样就形成了一个转换表，用于将灰度值为255的像素点变为0。
    '''
    invert_table = [0 if i == 255 else i for i in range(256)]
    '''
    在这一行代码中，通过point方法应用了之前创建的转换表。
    这个方法会对图像中的每个像素点进行变换，根据转换表将灰度值为255的像素点变为0，而其他灰度值保持不变。
    最终，得到了经过反转处理的图像img_inverted1。
    '''
    img_inverted1 = img_8bit.point(invert_table, 'L')

    # 将图像对象转换为NumPy数组
    img_array = np.asarray(img_inverted1)

    # 生成与输入文件相同的输出文件名
    output_filename = os.path.join(output_dir, os.path.basename(img_path).replace('.png', ''))

    # 自动保存为npy文件
    np.save(output_filename, img_array)

4） npy转化png

5） png标签转nii图像

参考: ※※※【医学影像数据处理】nii 数据格式文件操作汇总、※※nii、npz、npy、dcm、mhd 的数据格式互转，及多目标分割处理、 ※ 将png格式文件转化为nii.gz格式、png标签转nii图像、多张png合成为nii.gz、png图片组合成nii文件

step2: 重采样(resampling)

※仔细检查每一步的结果（只是绘图，如此简单但意义重大）※：Preprocess-of-CT-data
与自然图像不同，在医学影像中，人体部位真实的大小（成像大小）是非常重要的信息。因此例如在CT图像中，存在体素间距（spacing）和体素个数（Voxel count）两个指标：
成像大小 $= s p a c in g * V o x e l co u n t$ ，且成像大小保持不变。

1. 图像分辨率

通常所说的图片分辨率（image resolution）其实是指像素数（pixel count），通常表达为 横向多少个像素 * 纵向多少个像素。像 480x800 这样的表述其实本来应该叫做尺寸（dimensions）的，但是因为数字图片并没有物理的长宽的概念，叫做尺寸反而可能会引起误解。数字图片的"宽"（width）和"高"（height）并非物理意义的长度单位，而是在两个维度上图片包含的像素个数。比如 480x800 这样的图片是由横向 480 个像素、纵向 800 个像素（合计 384000个像素点）构成的。

2. 像素间距（Spacing）

“像素间距”（Pixel Spacing）是指两个像素之间的距离，spacing定义了图像像素的物理大小并且保证了实际距离测量的准确性。比如，如果知道x和y轴的像素间距为 0.4mm，那么在图像中的一条 10 像素的线就会有 4mm的长度。同样，由于知道图像像素中的宽和高(比如对于普通CT来说是 512×512)，就能够找到图像的实际尺寸了： $512 \times 0.4 mm = 204.8 mm$ 。对应的，当图像从2D扩展到3D时，“像素”会扩展成“体素”，“像素间距”（Pixel Spacing）也会扩展成“体素间距”（Voxel Spacing）。

综上，图像分辨率一般是指“像素或体素个数”，而spacing是指“体素间距”，知道了体素数和体素间距，就知道了一张CT的实际“尺寸”。

3. CT图像的重采样

1）原因及相关知识

医院/中心不同、CT设备不同、扫描人员不同、以及层厚的具体参数设置不同都会产生差异，而CNN无法理解体素间距，因此我们需要将所有医学影像的spacing重采样到一致。
影像扫描得到X，Y，Z三个方向的体素间距，X，Y方向的体素间距较小，Z方向的体素间距略大，我们把这种情况称之为各向异性，而重采样至 $1*1*1mm^{3}$ 后，各向异性的图像变为各相同性。

① CT图像里面是保存有图像的pixdim信息的。查看pixdim信息的代码如下：

	#img_name 就是一个nii.gz文件的地址
	img_name = '/media/wanghui/wh2021/KiTS23/volume/case_0000/imaging.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)

结果：
实际spacing： $(0.5 mm, 0.9199219 mm, 0.9199219 mm)$ 。

② 由于不同的病人体型不同，但最后数字成像的分辨率是一样的，这就导致了一定程度的失真变形。但医学图像例如dcm或nii格式，都会带有SliceThickness，PixelSpacing类似的属性，可以利用这些属性去尽量还原真实物体。重采样(resampling)
eg：
若图片分辨率大小（像素pixel）digital_image.shape $= 2 * 2$ ，spacing $= (2, 2)$ ，则数字图像大小digital_image.shape $*$ spacing $= 4 * 4$ ，
现希望new_spacing $= (1, 1)$ ，即当前的像素尺寸refacor_size需要扩大spacing/new_spacing $= (2, 2)$ 倍。
公式：

refacor_size $=$ spacing/new_spacing $=$ new_shape/digital_image.shape
new_shape $=$ round(digital_image.shape $*$ refacor_size)
spacing很多时候都是浮点数，数字图像的分辨率是没有浮点数的，所以需要round处理。

注：
——如果现在希望spacing变为 $(1, 1)$ ，意思是希望数字成像的图片大小就代表原始物体大小。
即new_spacing $= (1, 1)$ ，则refacor_size $=$ spacing，那么 new_shape $=$ digital_image.shape $*$ spacing=数字图像大小。

——若现在希望spacing $= (2, 2)$ ，则表示原始物体大小是数字图像大小的两倍。

2）重采样

重采样的功能：
——将所有病例的体素间距重采样为 3.22 x 1.62 x 1.62。
——矫正图像的direction坐标

使用monai框架中分辨率变换 [Spacingd] ，代码如下：

#对图像重采样到spacing=1.5mm×1.5mm×5.0m
spacing = Spacingd(keys=["image", "label"], pixdim=(1.5, 1.5, 5.0), mode=("bilinear", "nearest"))

注：加载数据有两种方法。使用 Monai 和 PyTorch 预处理 3D Volumes以进行肿瘤分割

1- 第一种是单独加载图像和掩码（如果您想进行图像分类，可以使用这种方式，但它也适用于分割）。

2- 第二种方法是创建一个带有两列的Python字典，一列用于图像路径，一列用于标签路径。然后在每一行输入具有相应掩码的图像的路径。（推荐）

使用SimpleITK库，代码如下：
方法1：

# Step1： 定义转换函数transform
import numpy as np
import SimpleITK as sitk


def transform(image,newSpacing, resamplemethod=sitk.sitkNearestNeighbor):
    # 设置一个Filter
    resample = sitk.ResampleImageFilter()
    # 初始的体素块尺寸
    originSize = image.GetSize()
    # 初始的体素间距
    originSpacing = image.GetSpacing()
    newSize = [
        int(np.round(originSize[0] * originSpacing[0] / newSpacing[0])),
        int(np.round(originSize[1] * originSpacing[1] / newSpacing[1])),
        int(np.round(originSize[2] * originSpacing[2] / newSpacing[2]))
    ]
    print('current size:',newSize)

    # 沿着x,y,z,的spacing（3）
    # The sampling grid of the output space is specified with the spacing along each dimension and the origin.
    resample.SetOutputSpacing(newSpacing)
    # 设置original
    resample.SetOutputOrigin(image.GetOrigin())
    # 设置方向
    resample.SetOutputDirection(image.GetDirection())
    resample.SetSize(newSize)
    # 设置插值方式
    resample.SetInterpolator(resamplemethod)
    # 设置transform
    resample.SetTransform(sitk.Euler3DTransform())
    # 默认像素值   resample.SetDefaultPixelValue(image.GetPixelIDValue())
    return resample.Execute(image)

方法2： 与方法1在最后面有些不同的。
当只改变z轴的sapcing时只需要设置setSpacing（）就可以了。

def resample_image(itk_image, out_spacing,out_direction, is_label=False):
    original_size = itk_image.GetSize()
    original_spacing = itk_image.GetSpacing()
    out_size = [
        int(np.round(original_size[0] * (original_spacing[0] / out_spacing[0]))),
        int(np.round(original_size[1] * (original_spacing[1] / out_spacing[1]))),
        int(np.round(original_size[2] * (original_spacing[2] / out_spacing[2])))
    ]
    resample = sitk.ResampleImageFilter()
    resample.SetOutputSpacing(out_spacing)
    resample.SetSize(out_size)
    resample.SetOutputDirection(out_direction)
    resample.SetOutputOrigin(itk_image.GetOrigin())


    if is_label:
        resample.SetDefaultPixelValue(0)  # 没有图像的地方填充值
        resample.SetInterpolator(sitk.sitkNearestNeighbor)
    else:
        resample.SetDefaultPixelValue(-10)  # -10是我调整的窗宽窗外
        resample.SetInterpolator(sitk.sitkBSpline)

    return resample.Execute(itk_image)

方法3：

def resample(image, scan, new_spacing=[1,1,1]):
    image = imgs_to_process
    scan = patient
    # Determine current pixel spacing
    spacing = map(float, ([scan[0].SliceThickness] + list(scan[0].PixelSpacing)))
    spacing = np.array(list(spacing))

    resize_factor = spacing / new_spacing

    new_real_shape = image.shape * resize_factor
    # 转成浮点数，对给定的数组进行四舍五入。
    new_shape = np.round(new_real_shape)

    real_resize_factor = new_shape / image.shape

    new_spacing = spacing / real_resize_factor
    # 功能是缩放数组，即使用order顺序的样条插值来缩放数组。如果real_resize_factor为float，每个轴的缩放是相同的。
    image = scipy.ndimage.interpolation.zoom(image, real_resize_factor)
    
    return image, new_spacing

print("Shape before resampling\t", imgs_to_process.shape) #Shape before resampling  (129, 512, 512)
imgs_after_resamp, spacing = resample(imgs_to_process, patient, [1,1,1])
print("Shape after resampling\t", imgs_after_resamp.shape) #Shape after resampling   (206, 292, 292)

当然，也可以不进行插值处理，因为插值后会导致图片间的大小不同，无法输入网络，我见到的一个处理是这样的：
只对轴向(z)进行插值，y,x固定大小缩放，
这样对于每一个数据，就只会有切片数量不同。
在输入网络时，每个数据随机取样n个slice数目相同的三维图片，输入网络即可。

上述resampling均定义在图像坐标系下，因为numpy数组完全不会感知到physical world的存在，scipy.ndimage和scikit-image也没有储存图像在世界坐标的几何信息，因此，对于这两种library而言，只有一种坐标系。 而itk和nibabel则对于每个图片都拥有其physical space的信息，一些基本的操作也会在physical space下进行。

# 只对轴向(z)进行插值，y,x固定大小缩放
new_real_shape = image.shape * resize_factor
改成new_real_shape = [image.shape[0] * resize_factor[0],image.shape[1],image.shape[2]]

注意：重采样的插值方法，我试过SimpleITK自带的多种插值方法，线性插值,三次插值以及B样条，比较发现B样条的效果是最好的。
因此，图像image采用sitk.sitkBSpline插值，标签segment采用sitk.sitkNearestNeighbor插值。
如果感兴趣可以自己尝试一下不同的插值方法，或者使用scipy等其他工具包进行重采样。

# Step2：接Step1
data_path = "/root/data/nnUNet_raw_data_base/nnUNet_raw_data/Task040_KiTS/imagesTr"

for path in sorted(os.listdir(data_path)):
    print(path)
    img_path = os.path.join(data_path,path)
    img_itk = sitk.ReadImage(img_path)

    print('origin size:', img_itk.GetSize())
    # B样条
    new_itk = transform(img_itk, [3.22, 1.62, 1.62], sitk.sitkBSpline) # sitk.sitkLinear
    sitk.WriteImage(new_itk, img_path)

print('images is resampled!')
print('-'*20)


label_path = "/root/data/nnUNet_raw_data_base/nnUNet_raw_data/Task040_KiTS/labelsTr"

for path in sorted(os.listdir(label_path)):
    print(path)
    img_path = os.path.join(label_path,path)
    img_itk = sitk.ReadImage(img_path)

    print('origin size:', img_itk.GetSize())
    # NearestNeighbor
    new_itk = transform(img_itk, [3.22, 1.62, 1.62])
    sitk.WriteImage(new_itk, img_path)

print('labels is resampled!')

补充：
当对 CT 图像进行重新取样以标准化体素间距时，切片数会发生变化，尤其是沿 Z 轴（切片轴）的切片数。

为什么在重新取样过程中切片数会发生变化？
CT 扫描沿 x、y 和 z 轴有特定的分辨率（体素间距）。当您对图像进行重新取样以实现统一的体素间距时，原始切片厚度（沿 Z 轴的切片间距）会被调整。下面是它对切片数的影响：

原始间距： CT 扫描的体素间距可能不均匀，尤其是在 Z 轴上（例如，切片厚度可能在 2.5 毫米到 5 毫米之间，而 X 轴和 Y 轴的间距可能在 0.7 毫米到 1 毫米之间）。
重新采样：如果将所有图像重新采样为统一的切片厚度（例如您的情况中为 1 毫米），则沿 Z 轴的切片数可能会发生变化，因为重新采样操作会对原始切片进行插值或压缩，以匹配新的统一间距。

这对切片数有何影响？

如果原始切片厚度大于新的目标厚度（如 5 毫米 -> 1 毫米）：重新取样会在原始切片之间插入新的切片，从而增加 Z 轴上的切片数量。
如果原始切片厚度小于新的目标厚度（例如 0.5mm -> 1mm）：重采样将压缩切片信息，减少 Z 轴上的切片数量。

举例说明：

原始扫描：100 个切片，切片厚度为 5 毫米（沿 Z 轴）。
重新取样为 1 毫米：您将得到大约 500 个切片（因为您实际上沿 Z 轴提高了分辨率）。
相反，如果原始扫描的切片厚度为 0.5 毫米，重新采样到 1 毫米将使切片数量减半。

以您的情况为例：
当重新取样 CT 图像以实现 Z 轴 1 毫米的均匀体素间距时（使用 para.slice_thickness = 1），切片数将根据每个病例的原始切片厚度而发生变化。如果原来沿 Z 轴的间距大于 1 毫米，重采样后切片数将会增加；如果原来的间距小于 1 毫米，切片数将会减少。

ct_array = ndimage.zoom(ct_array, (para.down_scale, para.down_scale, z_spacing_ratio), order=3)

z_spacing_ratio = ct.GetSpacing()[-1] / para.slice_thickness：该因子根据新的期望间距（para.slice_thickness）调整 Z 轴切片数。
该函数将根据原始的 ct.GetSpacing()[-1] 与 para.slice_thickness 的比较结果来改变 Z 轴的切片数。

3）spacing和patch size的权衡

注： target spacing 和patch size 需要不断尝试。
参考（推荐）：
确定input patch size, batch size以及创建U-Net网络架构、
如何设计task specific的3D U-Net结构:nnUNet中结构设计策略整理

We attempted to vary the target spacing for resampling as well as the patch size of our network architecture.

图像的spacing保持一致，图像中体素值的个数（即图像分辨率）却不一定相同。而分割网络的框架是固定的，一般是需要输入图片的分辨率大小是一致的。对于CT或者MRI图像来讲，图像是非常大的，又是一个三维图像，不可能全部输入网络中训练。要么把图像直接Resize到固定的尺寸，要么就是裁剪图像。一般采用裁剪方式，所以通常在训练的时候是用一个固定大小的patch从图像中裁剪采样（在这个过程中原图分辨率不一样是不影响的）。monai框架提供了非常多的裁剪模式，包括中心裁剪，前景裁剪、随机裁剪、滑动窗口遍历整个volume裁剪等等，同时图像不够大的话，也可以进行填充。如：

#patch的size为(192×192×16)
RandCropByPosNegLabeld(keys, label_key=keys[1], spatial_size=(192, 192, 16), num_samples=3),

注意：这些裁剪方式都要求数据格式为通道优先格式（必须有通道维度）。
————————————————
增大CT图像的spacing，其总体的像素/体素个数(pixel/Voxel count)会减少，因此一个patch中获取的细节信息减少，上下文信息增大；而减少spacing，其总体的pixel/Voxel count会增大，一个patch中获取的细节信息会增大，而上下文信息又会减少。因此，优化网络补丁大小中的上下文信息量与图像数据中保留的细节之间的权衡trade-off，对于获得理想的性能至关重要。

Optimizing the trade-off between the amount of contextual information in the networks patch size vs the details retained in the image data is crucial in obtaining ideal performance

step3: window transform

采用不同的窗口宽度和中心作为数据增强。

1. 认知

1）CT图像的Hu值变换

CT值属于医学领域的概念，通常称亨氏单位（hounsfield unit,HU),反映了组织对X射线的吸收程度。黑影表示低吸收区，即低密度区，如含气体多的肺部；白影表示高吸收区，即高密度区，如骨骼。
————————————————
灰度值属于计算机领域的概念，指的是单个像素点的亮度，灰度值越大表示越亮。范围一般从0到255，白色为255，黑色为0，故黑白图片也称为灰度图像。

注：

无论对于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操作。
测试代码：医学图像预处理（三）——windowing(ct对比增强）

2）dicom转化为hu值

公式与代码： $H u = p i x e l v a l u e * s l o p e + in t erce pt$ ，其中slope,intercept可以从元数据中读取。

def get_pixels_hu(scans):
    #type(scans[0].pixel_array)
    #Out[15]: numpy.ndarray
    #scans[0].pixel_array.shape
    #Out[16]: (512, 512)
    # image.shape: (129,512,512)
    image = np.stack([s.pixel_array for s in scans])
    # Convert to int16 (from sometimes int16), 
    # should be possible as values should always be low enough (<32k)
    image = image.astype(np.int16)   # image.shape = (666, 512, 512)

    # Set outside-of-scan pixels to 1
    # The intercept is usually -1024, so air is approximately 0
    # CT扫描边界之外的灰度值是固定的，为-2000，需要把这些值设置为0
    image[image == -2000] = 0
    
    # Convert to Hounsfield units (HU)，转换为HU，就是灰度值*rescaleSlope+rescaleIntercept
    intercept = scans[0].RescaleIntercept
    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)

2. window transform

Hu的范围一般来说很大，这就导致了对比度很差，如果需要针对具体的器官进行处理，效果会不好，于是就有了windowing的方法：

1）认知

图像的亮度取决于window level/center（WL），图像的对比度取决于window width（WW）。
窗口宽度就是一幅CT图片可以显示的CT值范围。CT图像将窗宽范围内的CT值划分为16个灰阶进行显示，例如，CT图像范围为80HU，划分为16个灰阶，则80/16=5HU，在CT图像上，只有CT值相差5HU以上的组织才可以观察到。

A wider window width (2000 HU), therefore, will display a wider range of CT numbers. Consequently, the transition of dark to light structures will occur over a larger transition area to that of a narrow window width (<1000 HU).

窗位大小是窗宽上、下限的平均值。窗口level减少，图片将变亮，level增大，图片变暗。
当给定了window width和window level后，就能计算出窗口的上下界。
超过上界的，是白色；低于下界的，是黑色。
下面是计算方法和举例。

the upper grey level (x) is calculated via $W L + (WW \div 2)$
the lower grey level (y) is calculated via $W L - (WW \div 2)$

eg: a brain is W:80 L:40, therefore, all values above +80 will be white and all values below 0 are black.

2）window transform

强制截断cut off，即Hu值截断来代替此步；
eg: 每个case的强度值范围是[-79，304]。
window transform（推荐），因为无需分析整个数据集获取cut off所必须的截断值。
window transform代码：

def window_transform(self, windowWidth, windowCenter, normal=False):
        """
        注意，这个函数的self.image一定得是float类型的，否则就无效！
        return: trucated image according to window center and window width
        """
        minWindow = float(windowCenter) - 0.5*float(windowWidth)
        newimg = (self.image - minWindow) / float(windowWidth)
        newimg[newimg < 0] = 0  
        newimg[newimg > 1] = 1
        # 将值域转到0-255之间,例如要看头颅时，我们只需将头颅的值域转换到0-255就行了
        if not normal:
            newimg = (newimg * 255).astype('uint8')
        return newimg

结果：

before window transform： -3611.9658 ，1458.9244 （Note：the min,max value of CT data）
After window transform without normalization to [0,1]： -17.759829， 7.594622

注：统计tumor的最大值和最小值，计算窗宽窗位，所得结果最好。代码参考：CT图像之Hu值变换与窗宽窗位调整