CT影像文件格式 dicom pydicom

CT影像文件格式

转自  https://blog.csdn.net/Acmer_future_victor/article/details/106428407

CT图像的文件格式是 dicom 格式，可以用 pydicom 进行处理，其含有许多的DICOM Tag信息。查看一些tag信息的代码实现如下所示。

    # __author: Y
    # date: 2019/12/10
     
    import pydicom
    import numpy as np
    import matplotlib
    import pandas
    import SimpleITK as sitk
    import cv2
    from PIL import Image
     
    # 应用pydicom来提取患者信息
    def loadFile(filename):
        ds = sitk.ReadImage(filename)
        image_array = sitk.GetArrayFromImage(ds)
        frame_num, width, height = image_array.shape
        print('frame_num:%s, width:%s, height:%s'%(frame_num, width, height))
        return image_array, frame_num, width, height
     
     
    def loadFileInformation(filename):
        information = {}
        ds = pydicom.read_file(filename)
        information['PatientID'] = ds.PatientID
        information['PatientName'] = ds.PatientName
        information['PatientBirthDate'] = ds.PatientBirthDate
        information['PatientSex'] = ds.PatientSex
        information['StudyID'] = ds.StudyID
        information['StudyDate'] = ds.StudyDate
        information['StudyTime'] = ds.StudyTime
        information['InstitutionName'] = ds.InstitutionName
        information['Manufacturer'] = ds.Manufacturer
        information['NumberOfFrames'] = ds.NumberOfFrames
        print(information)
        return information
     
    loadFile('../000000.dcm')
    loadFileInformation('abdominallymphnodes-26828')

 
CT图像是根据人体不同组织器官对X射线的吸收能力不同扫描得到的，由许多轴向切片组成三维图像，从三个方向观察可以分为三个视图，分别是轴状图、冠状图和矢状图。运用pydicom读取dcm格式的CT图像切片的代码实现如下所示。

    def load_scan(path):
        # 获取切片
        slices = [pydicom.read_file(path + '/' + s) for s in os.listdir(path)]
        # 按ImagePositionPatient[2]排序，否则得到的扫描面是混乱无序的
        slices.sort(key=lambda x: int(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

 
为了更好地观察不同器官，需要将像素值转换为CT值，单位为HU。计算方法为HU=pixel*rescale slope+rescale intercept。其中，rescale slope和rescale intercept是dicom图像文件的两个tag信息。代码实现如下所示

    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)  # image.shape = (666, 512, 512)
        # Set outside-of-scan pixels to 0
        # 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
        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)

 
将像素值转换为CT值之后，可以设置窗宽、窗位来更好地观察不同组织、器官。每种组织都有一定的CT值或CT值范围，如果想观察这一特定组织，就将窗位设置为其对应的CT值，而窗宽是CT图像可以显示的CT值范围，窗位大小是窗宽上、下限的平均值。CT图像将窗宽范围内的CT值划分为16个灰阶进行显示，例如，CT图像范围为80HU，划分为16个灰阶，则80/16=5HU，在CT图像上，只有CT值相差5HU以上的组织才可以观察到。设置窗位、窗宽的代码实现如下所示。

    def get_window_size(organ_name):
        if organ_name == 'lung':
            # 肺部 ww 1500-2000 wl -450--600
            center = -500
            width = 2000
        elif organ_name == 'abdomen':
            # 腹部 ww 300-500 wl 30-50
            center = 40
            width = 500
        elif organ_name == 'bone':
            # 骨窗 ww 1000-1500 wl 250-350
            center = 300
            width = 2000
        elif organ_name == 'lymph':
            # 淋巴、软组织 ww 300-500 wl 40-60
            center = 50
            width = 300
        elif organ_name == 'mediastinum':
            # 纵隔 ww 250-350 wl 250-350
            center = 40
            width = 350
        return center, width
     
     
    def setDicomCenWid(slices, organ_name):
        img = slices
        center, width = get_window_size(organ_name)
        min = (2 * center - width) / 2.0 + 0.5
        max = (2 * center + width) / 2.0 + 0.5
     
        dFactor = 255.0 / (max - min)
        d, h, w = np.shape(img)
        for n in np.arange(d):
            for i in np.arange(h):
                for j in np.arange(w):
                    img[n, i, j] = int((img[n, i, j] - min) * dFactor)
     
        min_index = img < 0
        img[min_index] = 0
        max_index = img > 255
        img[max_index] = 255
     
        return img

 
CT图像不同扫描面的像素尺寸、粗细粒度是不同的，这对进行CNN有不好的影响，因此需要进行重构采样，将图像重采样为[1,1,1]的代码实现如下所示

    def resample(image, slice, new_spacing=[1, 1, 1]):
        spacing = map(float, ([slice.SliceThickness] + [slice.PixelSpacing[0], slice.PixelSpacing[1]]))
        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
        image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')
     
        return image, new_spacing

 
为了更好地进行网络训练，通常进行标准化，有min-max标准化和0-1标准化。
————————————————

原文链接：https://blog.csdn.net/Acmer_future_victor/article/details/106428407