dicom檔案處理的一些方式

匯入需要的模組

import os
import SimpleITK
import dicom
import numpy as np
import cv2
import glob
from tqdm import tqdm

首先需要匯入我們需要的處理的dicom檔案，dicom檔案是一組連續的圖片，我們根據圖片中的位置資訊對每張圖片進行間隔計算，然後把結果存到一個列表中，然後將圖片中的畫素資訊進行提取，縮放到1mm1mm1mm的尺度，get_cube_from_img這個函式是從影象中根據座標找到目標的中心，並且切一個包含目標的矩陣，然後把這個三維的矩陣平鋪開成一個64個2維的矩陣並儲存。歸一化的目的是為了加快模型收斂的速度，如果要儲存成灰度圖，需要畫素值乘以255.

def is_dicom_file(filename):
    '''
       判斷某檔案是否是dicom格式的檔案
    :param filename: dicom檔案的路徑
    :return:
    '''
    file_stream = open(filename, 'rb')
    file_stream.seek(128)
    data = file_stream.read(4)
    file_stream.close()
    if data == b'DICM':
        return True
    return False


def
 
 load_patient(src_dir):
    '''
        讀取某資料夾內的所有dicom檔案
    :param src_dir: dicom資料夾路徑
    :return: dicom list
    '''
    files = os.listdir(src_dir)
    slices = []
    for s in files:
        if is_dicom_file(src_dir + '/' + s):
            instance = dicom.read_file(src_dir + '/' + s)
            slices.
 
append(instance)
    slices.sort(key=lambda x: int(x.InstanceNumber))
    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


def get_pixels_hu_by_simpleitk(dicom_dir):
    '''
        讀取某資料夾內的所有dicom檔案,並提取畫素值(-4000 ~ 4000)
    :param src_dir: dicom資料夾路徑
    :return: image array
    '''
    reader = SimpleITK.ImageSeriesReader()
    dicom_names = reader.GetGDCMSeriesFileNames(dicom_dir)
    reader.SetFileNames(dicom_names)
    image = reader.Execute()
    img_array = SimpleITK.GetArrayFromImage(image)
    img_array[img_array == -2000] = 0
    return img_array


def rescale_patient_images(images_zyx, org_spacing_xyz, target_voxel_mm, is_mask_image=False):
    '''
        將dicom影象縮放到1mm:1mm:1mm的尺度
        :param images_zyx: 縮放前的影象(3維)
        :return: 縮放後的影象(3維)
    '''

    print("Spacing: ", org_spacing_xyz)
    print("Shape: ", images_zyx.shape)

    # print "Resizing dim z"
    resize_x = 1.0
    resize_y = float(org_spacing_xyz[2]) / float(target_voxel_mm)
    interpolation = cv2.INTER_NEAREST if is_mask_image else cv2.INTER_LINEAR
    res = cv2.resize(images_zyx, dsize=None, fx=resize_x, fy=resize_y, interpolation=interpolation)
    # print "Shape is now : ", res.shape

    res = res.swapaxes(0, 2)
    res = res.swapaxes(0, 1)
    # print "Shape: ", res.shape
    resize_x = float(org_spacing_xyz[0]) / float(target_voxel_mm)
    resize_y = float(org_spacing_xyz[1]) / float(target_voxel_mm)

    # cv2 can handle max 512 channels..
    if res.shape[2] > 512:
        res = res.swapaxes(0, 2)
        res1 = res[:256]
        res2 = res[256:]
        res1 = res1.swapaxes(0, 2)
        res2 = res2.swapaxes(0, 2)
        res1 = cv2.resize(res1, dsize=None, fx=resize_x, fy=resize_y, interpolation=interpolation)
        res2 = cv2.resize(res2, dsize=None, fx=resize_x, fy=resize_y, interpolation=interpolation)
        res1 = res1.swapaxes(0, 2)
        res2 = res2.swapaxes(0, 2)
        res = np.vstack([res1, res2])
        res = res.swapaxes(0, 2)
    else:
        res = cv2.resize(res, dsize=None, fx=resize_x, fy=resize_y, interpolation=interpolation)

    res = res.swapaxes(0, 2)
    res = res.swapaxes(2, 1)

    print("Shape after: ", res.shape)
    return res


def get_cube_from_img(img3d, center_x, center_y, center_z, block_size):
    start_x = max(center_x - block_size / 2, 0)
    if start_x + block_size > img3d.shape[2]:
        start_x = img3d.shape[2] - block_size

    start_y = max(center_y - block_size / 2, 0)
    start_z = max(center_z - block_size / 2, 0)
    if start_z + block_size > img3d.shape[0]:
        start_z = img3d.shape[0] - block_size
    start_z = int(start_z)
    start_y = int(start_y)
    start_x = int(start_x)
    res = img3d[start_z:start_z + block_size, start_y:start_y + block_size, start_x:start_x + block_size]
    return res


def normalize_hu(image):
    '''
    將輸入影象的畫素值(-4000 ~ 4000)歸一化到0~1之間
    :param image 輸入的影象陣列
    :return: 歸一化處理後的影象陣列
    '''
    MIN_BOUND = -1000.0
    MAX_BOUND = 400.0
    image = (image - MIN_BOUND) / (MAX_BOUND - MIN_BOUND)
    image[image > 1] = 1.
    image[image < 0] = 0.
    return image


def load_patient_images(src_dir, wildcard="*.*", exclude_wildcards=[]):
    '''
    讀取一個病例的所有png影象，返回值為一個三維影象陣列
    :param image 輸入的一系列png影象
    :return: 三維影象陣列
    '''
    src_img_paths = glob.glob(src_dir + wildcard)
    for exclude_wildcard in exclude_wildcards:
        exclude_img_paths = glob.glob(src_dir + exclude_wildcard)
        src_img_paths = [im for im in src_img_paths if im not in exclude_img_paths]
    src_img_paths.sort()
    images = [cv2.imread(img_path, cv2.IMREAD_GRAYSCALE) for img_path in src_img_paths]
    images = [im.reshape((1,) + im.shape) for im in images]
    res = np.vstack(images)
    return res


def save_cube_img(target_path, cube_img, rows, cols):
    '''
        將3維cube影象儲存為2維影象,方便勘誤檢查
        :param 二維影象儲存路徑, 三維輸入影象
        :return: 二維影象
    '''
    assert rows * cols == cube_img.shape[0]
    img_height = cube_img.shape[1]
    img_width = cube_img.shape[1]
    res_img = np.zeros((rows * img_height, cols * img_width), dtype=np.uint8)

    for row in range(rows):
        for col in range(cols):
            target_y = row * img_height
            target_x = col * img_width
            res_img[target_y:target_y + img_height, target_x:target_x + img_width] = cube_img[row * cols + col]

    cv2.imwrite(target_path, res_img)
if __name__ == '__main__':
    dicom_dir = './data/dicom_demo/'
    # 讀取dicom檔案的元資料(dicom tags)
    slices = load_patient(dicom_dir)
    # 獲取dicom的spacing值
    pixel_spacing = slices[0].PixelSpacing
    pixel_spacing.append(slices[0].SliceThickness)
    print('The dicom spacing : ', pixel_spacing)
    # 提取dicom檔案中的畫素值
    image = get_pixels_hu_by_simpleitk(dicom_dir)
    # 標準化不同規格的影象尺寸, 統一將dicom影象縮放到1mm:1mm:1mm的尺度
    image = rescale_patient_images(image, pixel_spacing, 1.00)
    for i in tqdm(range(image.shape[0])):
        img_path = "./temp_dir/dcm_2_png/img_" + str(i).rjust(4, '0') + "_i.png"
        # 將畫素值歸一化到[0,1]區間
        org_img = normalize_hu(image[i])
        # 儲存影象陣列為灰度圖(.png)
        cv2.imwrite(img_path, org_img * 255)

    # 載入上一步生成的png影象
    pngs = load_patient_images("./temp_dir/dcm_2_png/", "*_i.png")
    # 輸入人工標記的結節位置: coord_x, coord_y, coord_z
    cube_img = get_cube_from_img(pngs, 272, 200, 134, 64)
    print(cube_img)
    save_cube_img('./temp_dir/chapter3_3dcnn_img_X.png', cube_img, 8, 8)