使用pydicom和SimpleITK预处理dicom文件

注意我安装了pydicom之后需要安装gdcm依赖,但我不能成功import gdcm,所以在下面的代码中都可能同时(混合)使用了pydicom和SimpleITK包读取的图像数据来做预处理。

可以成功import gdcm的同学请直接移步参考:https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial

(一)简介与可视化

对于医学CT图像来说,每个人有若干张(切片).dcm类型的图像存入同一个文件夹中,分别是不同角度/部位的医学成像,可以用mango/小赛之类的看图软件看一下,也可以用matlab可视化:

clear all;
clc;
close all;
file_path='.\images\001\';
img_path_list=dir(strcat(file_path,'*.dcm'));%获取该文件夹中所有dcm格式的图像
img_num=length(img_path_list);%获取图像总数量
imagename=img_path_list(1).name;%图像名
imagemax=dicominfo(strcat(file_path,imagename));  %读取图像
imweigth=imagemax.Rows;     %DICOM水平像素数量
imhigh=imagemax.Columns;    %DICOM垂直像素数量
imagemax1=dicomread(strcat(file_path,imagename));  %读取图像
I=imagemax1;

for j=2:img_num;  %逐一读取图像
    imagename=img_path_list(j).name;%图像名
    imagemax=dicomread(strcat(file_path,imagename));  %读取图像
    j
            I=imfuse(imagemax,I); %投影图像对应点取两幅图像合成值
end
I=rgb2gray(I);
m=max(max(I));
imshow(I,[]); %显示投影图像
imwrite(I,'001.png','png'); %保存投影图像

(二)读取

3D的CT不像2D图像可以直接读取 预处理 使用,可以借助python中的pydicom或者SimpleITK包来处理dicom。

SimpleITK

首先要将每个人的CT存储为一个大小为N*W*H的3D数组,其中N是这个人的.dcm切片图像数量,W*H是每张图的尺寸。

import SimpleITK as sitk

directorypath = './images/001/'#
def get_img_array(directorypath):
    reader = sitk.ImageSeriesReader()
    reader.MetaDataDictionaryArrayUpdateOn()
    reader.LoadPrivateTagsOn()
    series_IDs = sitk.ImageSeriesReader.GetGDCMSeriesIDs(directorypath)
    series_ID = series_IDs[0]
    dicom_names = reader.GetGDCMSeriesFileNames(directorypath,series_ID)
    reader.SetFileNames(dicom_names)
    image3D = reader.Execute()
    imgArray = sitk.GetArrayFromImage(image3D)
    return image3D, imgArray 

如果想要获取某个切片slice_index的详细信息:

slice_index = 1 # for example
key = reader.GetMetaDataKeys(slice_index)
print(key)
for key in reader.GetMetaDataKeys(slice_index):
    value = reader.GetMetaData(slice_index,key)
    print("({0}) = = \"{1}\"".format(key, value))

pydicom

(1)如果使用pydicom读取目录中的多个文件(它们都是DICOM文件),不能直接:

dicom.read_dicomdir(CT_files_dir)

这样会报错:

PermissionError: [Errno 13] Permission denied:

参考:https://github.com/pydicom/pydicom/issues/322

dirname = '/some/path'
files = os.listdir(dirname)
ds_list = [dicom.read_file(os.path.join(dirname, filename)) for filename in files]

(2)如果缺失提示meta information,读取时使用“,force=True”:

ds_list = [dicom.read_file(os.path.join(dirname, filename), force=True) for filename in files]

(3)如果有的图像缺失部分meta info,则无法获取ImagePositionPatient。

此时可以注释掉第二行#slices.sort(key = lambda x: float(x.ImagePositionPatient[2])):

def load_scan(path):
    slices = [pydicom.read_file(path + '/' + s,force=True) for s in os.listdir(path)]
    #slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))
    try:
        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])
    except:
        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)
    for s in slices:
        s.SliceThickness = slice_thickness     
    return slices

(三)预处理

(1)将dcm图像的值转换为HU值

CT扫描的测量单位是亨斯菲尔德单位(HU),这是一种辐射密度的测量。CT扫描仪经过仔细校准,以准确地测量这一点。

HU(Hounsfiled Unit)值,反映了组织对X射线吸收程度。以水的吸收程度作为参考,即水的HU=0,衰减系数大于水的为正,小于水的为负值。并以骨皮质和空气的HU值为上限和下限。

为了站在医生的角度看问题,所以必须将dcm图像的值转换为HU值。

def get_pixels_hu(slices):
    image = np.stack([s.pixel_array for s in slices])
    # Convert to int16 (from sometimes int16), 
    # should be possible as values should always be low enough (<32k)
    image = image.astype(np.int16)

    # Set outside-of-scan pixels to 0
    # The intercept is usually -1024, so air is approximately 0
    image[image == -2000] = 0
    
    # Convert to Hounsfield units (HU)
    for slice_number in range(len(slices)):
        
        intercept = slices[slice_number].RescaleIntercept
        slope = slices[slice_number].RescaleSlope
        
        if slope != 1:
            image[slice_number] = slope * image[slice_number].astype(np.float64)
            image[slice_number] = image[slice_number].astype(np.int16)
            
        image[slice_number] += np.int16(intercept)
    
    return np.array(image, dtype=np.int16)

因为我安装之后无法成功import gdcm,所以采用迂回战术,即同时使用SimpleITK和pydicom读取的图像sitkarrayslices

def get_pixels_hu(sitkarray, slices):
    #sys.exit(0)
    print(np.array(slices).shape) #maybe different from 'num'
    [num,w,h] = sitkarray.shape
    #print(sitkarray[0])
    image = np.stack([sitkarray[i] for i in range(num)])
    # Convert to int16 (from sometimes int16), 
    # should be possible as values should always be low enough (<32k)
    image = image.astype(np.int16)
    # Set outside-of-scan pixels to 0
    # The intercept is usually -1024, so air is approximately 0
    image[image == -2000] = 0
    # Convert to Hounsfield units (HU)
    for slice_number in range(num): 
        #print(slices[slice_number].RescaleIntercept,slice_number)  
        intercept = slices[slice_number].RescaleIntercept
        slope = slices[slice_number].RescaleSlope 
        if slope != 1:
            image[slice_number] = slope * image[slice_number].astype(np.float64)
            image[slice_number] = image[slice_number].astype(np.int16)    
        image[slice_number] += np.int16(intercept)
    return np.array(image, dtype=np.int16)

(2)重采样

扫描的像素间距如果是[2.5,0.5,0.5],这表示切片之间的距离为2.5毫米。对于不同的扫描,这项参数可能是不同的。

处理这一问题的一种常见方法是将整个数据集重新采样到一定的各向同性分辨率。

重采样到一个同构分辨率1mm *1mm *1mm,以消除扫描仪分辨率的差异:

def resample(sitk_img3D, image, scan, new_spacing=[1,1,1]):
    # Determine current pixel spacing
    #spacing = np.array(sitk_img3D.GetSpacing())
    spacing = np.array([scan[0].SliceThickness] + scan[0].PixelSpacing, dtype=np.float32) 
    #print(spacing)
    #spacing = np.array([scan[0].SliceThickness] + scan[0].PixelSpacing, dtype=np.float32)
    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
    image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')
    return image, new_spacing

(3)归一化

两个优点:

1)加快了梯度下降求最优解的速度;

2)有可能提高精度。

MIN_BOUND = -1000.0
MAX_BOUND = 400.0

def normalize(image):
    image = (image - MIN_BOUND) / (MAX_BOUND - MIN_BOUND)
    image[image>1] = 1.
    image[image<0] = 0.
    return image

(4)中心化(零均值化)

意义:数据中心化和标准化在回归分析中是取消由于量纲不同、自身变异或者数值相差较大所引起的误差。
原理:数据标准化:是指数值减去均值,再除以标准差;

数据中心化:是指变量减去它的均值。

目的:通过中心化和标准化处理,得到均值为0,标准差为1的服从标准正态分布的数据。

PIXEL_MEAN = 0.25 # To determine this mean you simply average all images in the whole dataset. If that sounds like a lot of work, we found this to be around 0.25 in the LUNA16 competition.
def zero_center(image):
    image = image - PIXEL_MEAN
    return image

 

函数调用:

sitk_img3D, sitk_img = get_img_array(directorypath)
dicom_img = load_scan(directorypath)
img_pixels = get_pixels_hu(sitk_img, dicom_img)
pix_resampled, spacing = resample(sitk_img3D, img_pixels, dicom_img, [1,1,1])

参考:

https://www.cnblogs.com/jxblog/p/12010354.html

https://github.com/pydicom/pydicom/issues/322

https://www.jianshu.com/p/f98635abac65

https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial

https://blog.csdn.net/weixin_37536336/article/details/109386431

https://blog.csdn.net/csdnxiekai/article/details/109448128

  • 6
    点赞
  • 50
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
使用SimpleITK库导入DICOM文件并生成STL文件,您可以按照以下步骤操作: ```python import SimpleITK as sitk import vtk def dicom_to_stl(dicom_dir, output_stl): # 读取DICOM序列 reader = sitk.ImageSeriesReader() dicom_names = reader.GetGDCMSeriesFileNames(dicom_dir) reader.SetFileNames(dicom_names) image = reader.Execute() # 使用Marching Cubes算法提取表面 surface_extractor = sitk.MarchingCubes() surface_extractor.SetInput(image) surface_extractor.SetValue(0, 1000) # 根据实际需求设置阈值 # 将表面转换为vtk数据类型 surface = surface_extractor.GetOutput() surface_array = sitk.GetArrayFromImage(surface) vtk_data = vtk.vtkImageData() vtk_data.SetDimensions(surface_array.shape[::-1]) vtk_data.AllocateScalars(vtk.VTK_UNSIGNED_CHAR, 1) vtk_data.GetPointData().SetScalars(vtk.util.numpy_support.numpy_to_vtk(surface_array.ravel(), deep=True)) # 使用vtkSTLWriter将表面写入STL文件 stl_writer = vtk.vtkSTLWriter() stl_writer.SetFileName(output_stl) stl_writer.SetInputData(vtk_data) stl_writer.Write() dicom_dir = "path/to/dicom/files" output_stl = "path/to/output.stl" dicom_to_stl(dicom_dir, output_stl) ``` 请确保您已经安装了SimpleITK和vtk库。在代码中,您需要将`dicom_dir`替换为包含DICOM文件的目录的路径,并将`output_stl`替换为要生成的STL文件的路径。 与之前的代码示例相比,这里使用SimpleITK库而不是pydicom库来读取DICOM文件。然后,我们使用SimpleITK的Marching Cubes算法提取表面,并将表面数据转换为vtk数据类型,最后使用vtkSTLWriter将表面写入STL文件。 请注意,阈值的设置可能需要根据您的DICOM文件的实际情况进行调整。 希望对您有帮助!
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值