通过pydicom在两张dicom图像之间进行线性插值

说明

如果两张dicom空间位置相距较远, 线性插值会有问题, 待研究新的解决方案

里面有将不同序列dicom合并(修改部分dicom标签)的部分代码, 待完善

重要重要重要

有一篇论文叫Nonlinear Image Interpolation using Manifold Learning 如果能实现其中的算法, 插值效果会好很多很多

dicom两层之间进行插值

import pydicom
import numpy as np
import matplotlib.pyplot as plt
import os
from pydicom.uid import ExplicitVRLittleEndian


# 打印dicom的一些tag信息
def print_tag_info(dataset):
    old_media_storage_uid = dataset.file_meta.MediaStorageSOPInstanceUID
    old_sop_uid = dataset.SOPInstanceUID
    print(type(old_sop_uid), old_sop_uid)
    print(type(old_media_storage_uid), old_media_storage_uid)
    print(type(dataset.pixel_array), dataset.pixel_array.dtype)
    print(dataset.pixel_array.shape)
    image_position_patient = dataset.ImagePositionPatient
    slice_thickness = dataset.SliceThickness
    print(type(image_position_patient), image_position_patient)
    print(type(slice_thickness), slice_thickness)
    slice_location = dataset.SliceLocation
    print(type(slice_location), slice_location)
    instance_number = dataset.InstanceNumber
    print(type(instance_number), instance_number)
    print('-'*20)
    return slice_location, slice_thickness


# 保存新的dicom文件
def save_new_dicom(new_ds, slice_number, image_position, out_dir, base_ds):
    # 修改新的dicom 的 tag 信息
    new_ds.file_meta.TransferSyntaxUID = ExplicitVRLittleEndian  # 由于构造的没有压缩的PixelData
    new_ds.ImagePositionPatient = image_position
    new_ds.SliceLocation = image_position[2]
    new_ds.InstanceNumber = slice_number
    old_sop_uid = base_ds.SOPInstanceUID
    old_sop_uid_arr = old_sop_uid.split('.')
    new_uid = ".".join(old_sop_uid_arr[:-1]) + '.' + str(slice_number)
    new_ds.file_meta.MediaStorageSOPInstanceUID = new_uid
    new_ds.SOPInstanceUID = new_uid

    new_ds.SeriesInstanceUID = base_ds.SeriesInstanceUID
    new_ds.SeriesNumber = base_ds.SeriesNumber

    new_ds.InstanceCreationTime = base_ds.InstanceCreationTime
    new_ds.ContentTime = base_ds.ContentTime
    new_ds.SeriesDescription = base_ds.SeriesDescription
    new_ds.DerivationDescription = base_ds.DerivationDescription
    new_ds.ConvolutionKernel = base_ds.ConvolutionKernel
    out_save_name = out_dir + str(slice_number) + '.dcm'
    print(out_save_name)
    new_ds.save_as(out_save_name)


# 在两层dicom之间线性插值
def insert_dicom(dicom_path1, dicom_path2, out_dir):
    ds1 = pydicom.dcmread(dicom_path1)
    ds2 = pydicom.dcmread(dicom_path2)

    slice_location1, slice_thickness1 = print_tag_info(ds1)
    slice_location2, slice_thickness2 = print_tag_info(ds2)
    insert_num = (slice_location2-slice_location1)/slice_thickness1
    insert_num = int(abs(insert_num))
    old_sop_uid = ds1.SOPInstanceUID
    last_name = int(old_sop_uid.split('.')[-1])
    print(type(last_name), last_name, insert_num)
    image_position = ds1.ImagePositionPatient

    for i in range(1, insert_num):
        ds1 = pydicom.dcmread(dicom_path1)
        new_image_position = image_position
        new_image_position[2] = new_image_position[2] - slice_thickness1
        print(new_image_position)
        slice_index = last_name + i
        left_gravity = i/insert_num

        new_pixel_array = ds1.pixel_array*(1-left_gravity)+ds2.pixel_array*left_gravity
        new_pixel_array = new_pixel_array.astype(np.int16)
        new_ds = ds1.copy()
        new_ds.PixelData = new_pixel_array.tobytes()
        print('generate new image...', i, left_gravity)
        save_new_dicom(new_ds, slice_index, new_image_position, out_dir, ds1)
    print('end...')


path1 = 'F:/400.dcm'
path2 = 'F:/650.dcm'
outdir = 'F:/44new/'
# insert_dicom(path1, path2, outdir)


# 重新保存指定文件夹下的dicom
def re_save_dicom(dicom_dir, new_dir):
    for root, dirs, files in os.walk(dicom_dir):
        for file in files:
            file_name = os.path.join(root, file)
            new_file_dir = new_dir + root[9:]
            if not os.path.exists(new_file_dir):
                os.makedirs(new_file_dir)
                print("Folder created successfully")
            dicom_image = pydicom.dcmread(file_name)
            new_file_name = os.path.join(new_file_dir, file)
            dicom_image.save_as(new_file_name)
            print('save succeed...', new_file_name)

合并多个dicom序列为一个序列

g_slice_index = 1


# 修改dicom标签并重新保存
def change_tag_and_re_save(base_dicom_path, dicom_files_dataset, out_dir):
    global g_slice_index
    base_ds = pydicom.dcmread(base_dicom_path)
    for file_dataset in dicom_files_dataset:
        new_ds = file_dataset
        image_position = base_ds.ImagePositionPatient
        image_position[2] = g_slice_index
        new_ds.ImagePositionPatient = image_position
        new_ds.SliceLocation = g_slice_index
        new_ds.InstanceNumber = g_slice_index

        old_sop_uid = base_ds.SOPInstanceUID
        old_sop_uid_arr = old_sop_uid.split('.')
        new_uid = ".".join(old_sop_uid_arr[:-1]) + '.' + str(g_slice_index)
        new_ds.file_meta.MediaStorageSOPInstanceUID = new_uid
        new_ds.SOPInstanceUID = new_uid

        new_ds.PatientName = base_ds.PatientName
        new_ds.PatientID = base_ds.PatientID
        new_ds.StudyID = base_ds.StudyID
        new_ds.AccessionNumber = base_ds.AccessionNumber
        new_ds.StudyInstanceUID = base_ds.StudyInstanceUID
        new_ds.SeriesInstanceUID = base_ds.SeriesInstanceUID
        new_ds.SeriesNumber = base_ds.SeriesNumber
        new_ds.IrradiationEventUID = base_ds.IrradiationEventUID
        new_ds.FrameOfReferenceUID = base_ds.FrameOfReferenceUID

        new_ds.StudyTime = base_ds.StudyTime
        new_ds.SeriesTime = base_ds.SeriesTime
        new_ds.AcquisitionTime = base_ds.AcquisitionTime
        new_ds.InstanceCreationTime = base_ds.InstanceCreationTime
        new_ds.ContentTime = base_ds.ContentTime
        new_ds.SeriesDescription = base_ds.SeriesDescription
        new_ds.DerivationDescription = base_ds.DerivationDescription
        new_ds.ConvolutionKernel = base_ds.ConvolutionKernel

        new_file_name = os.path.join(out_dir, str(g_slice_index)+'.dcm')
        g_slice_index = g_slice_index+1
        new_ds.save_as(new_file_name)
        print('save succeed...', new_file_name)


# 对指定文件夹下的dicom进行排序
def sort_dicom_dir(dicom_dir):
    all_file_dataset = {}

    for root, dirs, files in os.walk(dicom_dir):
        for file in files:
            file_name = os.path.join(root, file)
            dcm_dataset = pydicom.dcmread(file_name)
            if not hasattr(dcm_dataset, 'SliceLocation'):
                continue

            if dcm_dataset.SeriesInstanceUID in all_file_dataset:
                all_file_dataset[dcm_dataset.SeriesInstanceUID].append(dcm_dataset)
            else:
                all_file_dataset[dcm_dataset.SeriesInstanceUID] = []
                all_file_dataset[dcm_dataset.SeriesInstanceUID].append(dcm_dataset)

    for key in all_file_dataset:
        print(key, len(all_file_dataset[key]))
        new_slices = sorted(all_file_dataset[key], key=lambda s: s.SliceLocation)
        all_file_dataset[key] = new_slices
    return all_file_dataset


# 合并多个序列为一个序列
def merge_multi_series(base_path, multi_series, out_dir):
    if not os.path.exists(out_dir):
        os.makedirs(out_dir)

    for series in multi_series:
        print(series)
        all_file_dataset = sort_dicom_dir(series)
        for key in all_file_dataset:
            change_tag_and_re_save(base_path, all_file_dataset[key], out_dir)


base_dcm_path = r'F:\1.dcm'
series_dirs = [r'F:\only']
one_dir = r'F:\one'
merge_multi_series(base_dcm_path, series_dirs, one_dir)