使用pydicom和SimpleITK预处理dicom文件

梦里梦外; 2023-01-10 11:45 187阅读 0赞

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

可以成功import gdcm的同学请直接移步参考：[https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial][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'); %保存投影图像

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L01JS0FTQTM_size_16_color_FFFFFF_t_70][]

## （二）读取 ##

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][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读取的图像**sitkarray**和**slices**：

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_www.cnblogs.com_jxblog_p_12010354.html]

[https://github.com/pydicom/pydicom/issues/322][https_github.com_pydicom_pydicom_issues_322]

[https://www.jianshu.com/p/f98635abac65][https_www.jianshu.com_p_f98635abac65]

[https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial][https_www.kaggle.com_gzuidhof_full-preprocessing-tutorial]

[https://blog.csdn.net/weixin\_37536336/article/details/109386431][https_blog.csdn.net_weixin_37536336_article_details_109386431]

[https://blog.csdn.net/csdnxiekai/article/details/109448128][https_blog.csdn.net_csdnxiekai_article_details_109448128]

[https_www.kaggle.com_gzuidhof_full-preprocessing-tutorial]: https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L01JS0FTQTM_size_16_color_FFFFFF_t_70]: /images/20221119/dafd79d8734044f49d98f0b429aff56a.png
[https_github.com_pydicom_pydicom_issues_322]: https://github.com/pydicom/pydicom/issues/322
[https_www.cnblogs.com_jxblog_p_12010354.html]: https://www.cnblogs.com/jxblog/p/12010354.html
[https_www.jianshu.com_p_f98635abac65]: https://www.jianshu.com/p/f98635abac65
[https_blog.csdn.net_weixin_37536336_article_details_109386431]: https://blog.csdn.net/weixin_37536336/article/details/109386431
[https_blog.csdn.net_csdnxiekai_article_details_109448128]: https://blog.csdn.net/csdnxiekai/article/details/109448128