医学图像预处理----①小白入门

真的是纯小白，一看代码就头疼，但是任务来了也不能不做，毕竟还是要毕业的。所以呢我就搜啊搜啊搜，终于经过老师指导与查阅文献终于有点处理思路，已完成预处理所有的工作在此做个汇总。本次预处理主要为图像配准时使用

1、dicom格式转Niii。

这一块比较简单，直接上代码了

import numpy as np
import shutil
import os
import SimpleITK as sitk

def dcm2nii_sitk(path_read, path_save):
    reader = sitk.ImageSeriesReader()
    seriesIDs = reader.GetGDCMSeriesIDs(path_read)
    N = len(seriesIDs)
    lens = np.zeros([N])
    for i in range(N):
        dicom_names = reader.GetGDCMSeriesFileNames(path_read, seriesIDs[i])
        lens[i] = len(dicom_names)
    N_MAX = np.argmax(lens)
    dicom_names = reader.GetGDCMSeriesFileNames(path_read, seriesIDs[N_MAX])
    reader.SetFileNames(dicom_names)
    image = reader.Execute()
    if not os.path.exists(path_save):
        os.mkdir(path_save)
    sitk.WriteImage(image, path_save+'/data.nii.gz')

DICOMpath = r"DATA_MD/Changhai_Int_anonymize/dc4f80590509403da78a763bd7f64baa/0a324d770c46429b9e2e08fee5b88761"   #dicom文件夹路径
Midpath = r"SavePng"   #处理中间数据路径
Resultpath = r"SaveRaw"    #保存路径
cases = os.listdir(DICOMpath)  #获取dicom文件夹路径子文件夹名
for c in cases:  #遍历dicom文件夹路径子文件
    path_mid = os.path.join(DICOMpath , c)  #获取dicom文件夹下每一套数据的路径
    dcm2nii_sitk(path_mid , Midpath )  #将dicom转换为nii，并保存在Midpath中
    shutil.copy(os.path.join(Midpath , "data.nii.gz"), os.path.join(Resultpath , c + ".nii.gz"))
    #重新对保存后的nii文件名进行命名，并复制到Resultpath下

2、显示NII格式代码用于整理数据集。（窗体窗位显示的更好）

下面展示一些 第一次显示窗位的代码。

import torchio as tio
import numpy as np
import SimpleITK as sitk
import os
import nibabel as nib
from matplotlib import pyplot as plt
old_path="put/P_06"
cases = os.listdir(old_path)
i=0
def window(img):
    win_min = -300
    win_max = 350
    for i in range(img.shape[0]):
        img[i] = 255.0 * (img[i] - win_min) / (win_max - win_min)
        min_index = img[i] < 0
        img[i][min_index] = 0
        max_index = img[i] > 255
        img[i][max_index] = 255
        img[i] = img[i] - img[i].min()
        if img[i].max()!=0 :
           c = float(255) / img[i].max()
           img[i] = img[i] * c
    return img.astype(np.uint8)
for c in cases:  #遍历p文件子文件
    i = i + 1
    path_mid = os.path.join(old_path ,c)  #获取dicom文件夹下每一套数据的路径
    img = sitk.ReadImage(path_mid)
    img = sitk.GetArrayFromImage(img)
    img = window(img)
    out = sitk.GetImageFromArray(img)
    sitk.WriteImage(out, '{0}.nii.gz'.format(i))
    t1_img = tio.ScalarImage('{0}.nii.gz'.format(i))
    t1_img.plot()
    img = nib.load(path_mid)
    img_hdr = img.header
    voxel = [img_hdr['qoffset_x']]
    voxel.append(img_hdr['qoffset_y'])
    voxel.append(img_hdr['qoffset_z'])
    print("----{0}的头文件信息是,坐标信息是------\n{1}\n{2}".format(c,img,voxel))
print(
i)

下面展示一些 优化显示窗位的代码。在调用函数之前写好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

3、制作表格显示主要信息，便于数据集的使用（多套数据集）

4、对数据集进行切片处理减少干扰（可利用ITK-SNAP查看）

下面展示一些 遍历三维图像找到边界值。

import matplotlib.pyplot as plt
import os
import SimpleITK as sitk
import cv2
old_path="DATA_MD_01/p_18"
cases = os.listdir(old_path)
ori_data=sitk.ReadImage(os.path.join(old_path,cases[0]))
pic_Size=ori_data.GetSize()
data=sitk.GetArrayFromImage(ori_data)

y_list = []
x_list = []
def cv_show(img, name):
    cv2.imshow(name, img)
    cv2.waitKey()
    cv2.destroyAllWindows()
def norm_img(image):  # 归一化像素值到（0，1）之间，且将溢出值取边界值
    image = (image - MIN_BOUND) / (MAX_BOUND - MIN_BOUND)
    image[image > 1] = 1.
    image[image < 0] = 0.
    return image
for i in range(0,pic_Size[2]):
    image=data[i,:,:]
    MIN_BOUND = -1000.0
    MAX_BOUND = 700.0
    image=norm_img(image)
    plt.imshow(image)
    plt.show()
    draw_img = image.copy()
    ret, thresh = cv2.threshold(draw_img, 0.2, 255, cv2.THRESH_BINARY)
    for x in range(0,512):
        for y in range(0,512):
            if thresh[x][y]==255:
                y_list.append(x)
                break
    for y in range(0,512):
        for x in range(0,512):
            if thresh[x][y]==255:
                x_list.append(y)
                break
    x_list.sort()
    y_list.sort()
    x_max_list=x_list[len(x_list)-1]
    y_max_list=y_list[len(y_list)-1]
    x_min=x_list[0]
    y_min=y_list[0]
print(x_min,x_max_list)
print(y_min,y_max_list)

下面展示一些 通过边界值找到切分边界对多余信息切除。

import torchio as tio
import numpy as np
import SimpleITK as sitk
import os

def copy_geometry(image: sitk.Image, ref: sitk.Image):
    image.SetOrigin(ref.GetOrigin())
    image.SetDirection(ref.GetDirection())
    image.SetSpacing(ref.GetSpacing())
    return image
old_path="DATA_MD_01/p_03"
cases = os.listdir(old_path)
print(cases)
i=-1
for c in cases:
    i=i+1
    # print(c)
    path_img=os.path.join(old_path,c)
    img=sitk.ReadImage(path_img)
    # img = sitk.GetArrayFromImage(image)
    # img = sitk.GetImageFromArray(img)
    # img_out_itk = copy_geometry(img, image)
    direction = img.GetDirection()
    pic_Origin=img.GetOrigin()
    pic_Spacing=img.GetSpacing()
    pic_Size=img.GetSize()
    print("原图位置：{}\n原图体素：{}\n原图大小{}".format(pic_Origin,pic_Spacing,pic_Size))
    data_clip = img[10:510,60:390,:]
    data_clip.SetSpacing(pic_Spacing)
    data_clip.SetOrigin(pic_Origin)
    data_clip.SetDirection(direction)
    # print(data_clip)
    Resultpath = r"ab"  # 保存路径
    sitk.WriteImage(data_clip, os.path.join(Resultpath, c + ".gz"))
    t1_img = tio.ScalarImage(os.path.join(Resultpath, c + ".gz"))
    t1_img.plot()
    print("裁剪后位置：{}\n裁剪后体素：{}\n裁剪后大小{}".format(pic_Origin,pic_Spacing,pic_Size))

5、重采样得到训练所需体素大小，归一化处理

下面展示一些 遍历切除后生成的.nii格式文件进行重采样并且通过填充的形式进行调节大小格式一致。
对喽，这个地方的重采样貌似有好多方法，我自己就尝试了三种方法，把感觉最好用的拿出来分享一下

import numpy as np
import os
import SimpleITK as sitk
import torchio as tio


'''调窗'''
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

'''处理-1000——700数据'''
def pretreatmentImage(image):
    image[image < -1000] = -1000
    image[image > 700] = 700
    return image

'''归一化操作'''
def MinMax(image):
    if np.max(image) - np.min(image) != 0:
        image =(image - np.min(image)) / (np.max(image) - np.min(image))  # 归一化
    return image

'''重采样'''
def img_resmaple(ori_data,new_spacing=[2.5, 2.5, 2.5]):

    original_spacing = ori_data.GetSpacing()
    original_size = ori_data.GetSize()
    transform=sitk.Transform()
    transform.SetIdentity()
    new_shape = [
        int(np.round(original_spacing[0] * original_size[0] / new_spacing[0])),
        int(np.round(original_spacing[1] * original_size[1] / new_spacing[1])),
        int(np.round(original_spacing[2] * original_size[2] / new_spacing[2])),
    ]
    # print("新的形状大小为",new_shape)
    resmaple = sitk.ResampleImageFilter()#初始化 #初始化一个图像重采样滤波器resampler
    resmaple.SetInterpolator(sitk.sitkLinear)
    resmaple.SetTransform(transform)
    resmaple.SetDefaultPixelValue(0)
    resmaple.SetOutputSpacing(new_spacing)
    resmaple.SetOutputOrigin(ori_data.GetOrigin())
    resmaple.SetOutputDirection(ori_data.GetDirection())
    resmaple.SetSize(new_shape)
    data = resmaple.Execute(ori_data)
    return data

def copy_geometry(image: sitk.Image, ref: sitk.Image):
    image.SetOrigin(ref.GetOrigin())
    image.SetDirection(ref.GetDirection())
    image.SetSpacing(ref.GetSpacing())
    return image


if __name__ == "__main__":
    '''读取数据'''
    old_path = "ab"
    cases = os.listdir(old_path)
    for c in cases:
        ori_data = sitk.ReadImage(os.path.join(old_path, c))
        data_res = sitk.GetArrayFromImage(ori_data)
        img_pre = pretreatmentImage(data_res)
        img_Min = MinMax(img_pre)
        img = sitk.GetImageFromArray(img_Min)
        img_out_itk = copy_geometry(img, ori_data)
        data_pre=img_resmaple(img_out_itk)
        data=sitk.GetArrayFromImage(data_pre)
        # print(data.shape)
        if data.shape[2]<200 or data.shape[1]<160 or data.shape[0]<168:
            # print("我进入for循环了")
            x_up=int((200-data.shape[2])/2)
            if data.shape[2]%2!=0:
                x_down=x_up+1
            else:
                x_down=x_up
            y_left=int((160-data.shape[1])/2)
            if data.shape[1]%2!=0:
                y_right=y_left+1
            else:
                y_right=x_up
            z_front= int((168 - data.shape[0]) / 2)
            if data.shape[0] % 2 != 0:
                z_after = z_front + 1
            else:
                z_after = z_front
            data_pad=np.pad(data,((z_front,z_after),(y_left,y_right),(x_up,x_down)),'constant', constant_values=(0,0))
        # print('------------------------',data_pad.shape)
        Resultpath = r"SaveRaw"  # 保存路径
        img = sitk.GetImageFromArray(data_pad)
        img_out_itk = copy_geometry(img, data_pre)
        sitk.WriteImage(img_out_itk, os.path.join(Resultpath,c))
        '''验证存储信息是否正确'''
        ori = sitk.ReadImage(os.path.join(Resultpath,c ))
        pic_Origin=ori.GetOrigin()
        pic_Spacing=ori.GetSpacing()
        pic_Size=ori.GetSize()
        print("--------{}文件预处理重采样后的信息------".format(c))
        print("裁剪后位置：{}\n裁剪后体素：{}\n裁剪后大小{}".format(pic_Origin, pic_Spacing, pic_Size))
        img=sitk.GetArrayFromImage(ori)
        # print("数组信息为",img)
        t1_img = tio.ScalarImage(os.path.join(Resultpath, c))
        t1_img.plot()

附加另一种重采样方法zoom或者resize；具体实现自己搞一下很简单

class Resize_img(Base):
    def __init__(self, shape):
        self.shape = shape

    def tf(self, img, k=0):
        if k == 1:
            img = resize(img, (img.shape[0], self.shape[0], self.shape[1], self.shape[2]),
                         anti_aliasing=False, order=0)
        else:
            img = resize(img, (img.shape[0], self.shape[0], self.shape[1], self.shape[2]),
                         anti_aliasing=False, order=3)
        return img

class Resize_img(Base):
    def __init__(self, shape):
        self.shape = shape

    def tf(self, img, k=0):
        if img.shape[-1] == 3:
            img = zoom(img, (1, (self.shape[0] / img.shape[1]), (self.shape[1]/img.shape[2]), 1), order=3)
        else:
            img = zoom(img, (1, (self.shape[0] / img.shape[1]), (self.shape[1] / img.shape[2])), order=3)
        #resize(img, (img.shape[0], self.shape[0], self.shape[1]), anti_aliasing=False, order=3)
        return img

良心点我也给大家附上吧

import numpy as np
import os
import SimpleITK as sitk
from scipy.ndimage import zoom
import collections
from skimage.transform import rescale, resize
from scipy import ndimage, misc
import torchio as tio



'''处理-1000——700数据'''
def pretreatmentImage(image):
    image[image < -1000] = -1000
    image[image > 700] = 700
    return image
'''归一化操作'''
def MinMax(image):
    if np.max(image) - np.min(image) != 0:
        image = (image - np.min(image)) / (np.max(image) - np.min(image))  # 归一化
    return image

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
'''重采样'''
def resize_01(img):
    original_spacing = ori_data.GetSpacing()
    img =zoom(img,(original_spacing[2]/2.5, original_spacing[1]/2.5, original_spacing[0]/2.5),order=3)
    print(img.shape)
    return img
def copy_geometry(image: sitk.Image, ref: sitk.Image):
    image.SetOrigin(ref.GetOrigin())
    image.SetDirection(ref.GetDirection())
    image.SetSpacing([2.5,2.5,2.5])
    return image

if __name__ == "__main__":
    w_width = 400
    w_center = 40
    old_path = "ab"
    cases = os.listdir(old_path)
    ori_data = sitk.ReadImage(os.path.join(old_path, cases[0]))
    data = sitk.GetArrayFromImage(ori_data)
    img_pre = pretreatmentImage(data)
    img_Min = MinMax(img_pre)
    img_re=resize_01(img_Min)
    img = sitk.GetImageFromArray(img_re)
    img_out_itk = copy_geometry(img,ori_data)
    Resultpath = r"SaveRaw"  # 保存路径
    sitk.WriteImage(img_out_itk, os.path.join(Resultpath,cases[0]))
    '''验证存储信息是否正确'''
    ori = sitk.ReadImage(os.path.join(Resultpath, cases[0]))
    print('保存的图像信息', ori)
    img = sitk.GetArrayFromImage(ori)
    # print("数组信息为",img)
    image = adjustMethod1(img, w_width, w_center)
    t1_img = tio.ScalarImage(os.path.join(Resultpath, cases[0]))
    t1_img.plot()

总结一下子，老师是真的牛我是真的菜，一会说一下我检索的来源与出处

如果感兴趣来自哪里就搜一下吧

• 我没用到这个，但是他给我提供了下处理流程的思路： 医学nii图像 预处理——图像裁剪 重采样 灰度区域 归一化 修改图像尺寸
• 我主要用这个调的窗体窗位Python CT图像预处理——nii格式读取、重采样、窗宽窗位设置
• 这是我模仿的这个哥哥的整体思路 CT医学图像的预处理（重采样）
• 这个是主要做输出用0填充背景让大小一致的numpy.pad()函数使用详解
  重采样方法之二ZOOM（当时没看明白多搜了几个自己看看吧）
• scipy.ndimage.zoom矩阵放缩
• Python scipy.ndimage.zoom用法及代码示例
• Python scipy ndimage.zoom()函数
• CT图像预处理之重采样

行了，要是有不对的记得留言告诉我毕竟我也是刚学，那我走了，时间2023年5月8日，当然我会跟着项目尽量到训练模型跑出整个项目出来啊。毕竟论文还是要写的