利用SimpleITK将Numpy转Dcm

目的

在临床和科研中，我们经常需要将numpy数组转换成dicom格式的文件，以方便对这些dicom数据进行处理，例如剂量计算，形变配准或轮廓的映射等等。

代码

话不多说，直接上代码（适合在jupyter-notebook里运行）

import SimpleITK as sitk
import sys, time, os
import numpy as np
#%matplotlib inline
import matplotlib.pyplot as plt
from glob import glob
import shutil
from PIL import Image
from functools import cmp_to_key

# 注意：在运行前，需要将原始dcm文件和欲输出的png文件拷贝到inpupath下的同一个患者id下
inputpath = './'
outpath = './'
subdir=[]
lstDirs=[]
lstDirs = sorted(os.listdir(inputpath))
for name in lstDirs:
    subdir.append(os.path.join(inputpath+name))
print(subdir)


#读取dicom文件名
def read_Dcmfile(path):
    series_reader = sitk.ImageSeriesReader()
    series_IDs = series_reader.GetGDCMSeriesIDs(path)#读取文件夹下的序列，序列通常代表患者的某套影像
    #print(series_IDs)
    series_file_names1 = series_reader.GetGDCMSeriesFileNames(path, series_IDs[0])
    #print(series_file_names)
    return series_file_names1

#读取dicom文件里的tag、pixel value等信息
def read_Dcmimage(filename):
    image_reader = sitk.ImageFileReader()
    image_reader.LoadPrivateTagsOn()
    image_list1 = []
    for file_name in series_file_names:
        #print(file_name)
        image_reader.SetFileName(file_name)
        image_list1.append(image_reader.Execute())
    return image_list1

def compare(a,b):
    if a<b:
        return -1
    elif a==b:
        return 0
    else:
        return 0

def mycmp(x,y):
    s1 = x.split('/')[-1].split('_')
    p1=int(s1[0])
    i1=int(s1[1]) 

    s2= y.split('/')[-1].split('_')
    p2=int(s2[0])
    i2=int(s2[1])

    if p1==p2:
        return compare(i1,i2)
    return compare(p1,p2)


#读入png文件将pixel value读成numpy格式
def read_Pngimage(path):
    pngfilenames = sorted(glob(os.path.join(path,'*_fake_B.png')),key=cmp_to_key(mycmp))
    #print(pngfilenames)
    pngimage_all = []
    for name in pngfilenames:
        #print(name)
        #pngimage = np.array(Image.open(name).convert('L')).astype(np.int16)
        pngimage = np.array(Image.open(name).convert('L'))
        pngimage = pngimage/255*3000-1000
#         #plt.figure(figsize=(10,10))
#         plt.imshow(pngimage, 'gray')#手动勾画的轮廓
#         plt.axis('off')
#         plt.show()
#         plt.close()
        pngimage_all.append(pngimage)
    return pngimage_all

def make_dir(dirpath):
    if not os.path.exists(dirpath):
        os.makedirs(dirpath)



for subdirname in subdir:
    pat_id = subdirname.split("/")[-1]
    dest = os.path.join(outpath,pat_id)
    make_dir(dest)
    
    series_file_names = read_Dcmfile(subdirname)
    #print(series_file_names)
    image_list = read_Dcmimage(series_file_names)
    for i in range(5):
        print(image_list[i].GetOrigin())
    
    pngimage_data = np.array(read_Pngimage(subdirname))
    #print(pngimage_data.shape)
    
    filtered_image = sitk.GetImageFromArray(pngimage_data)
    filtered_image.SetSpacing(image_list[0].GetSpacing())
    filtered_image.SetDirection(image_list[0].GetDirection())
    
    writer = sitk.ImageFileWriter()
    writer.KeepOriginalImageUIDOn()
    modification_time = time.strftime("%H%M%S")
    modification_date = time.strftime("%Y%m%d")
    for i in range(filtered_image.GetDepth()):
        image_slice = filtered_image[:,:,i]
        original_slice = image_list[i]
    # Copy the meta-data except the rescale-intercept, rescale-slope为什么不拷贝？
        for k in original_slice.GetMetaDataKeys():
            if k!="0028|1052" and k!= "0028|1053":
                image_slice.SetMetaData(k, original_slice.GetMetaData(k))
    # Set relevant keys indicating the change, modify or remove private tags as needed
        image_slice.SetMetaData("0008|0031", modification_time)
        image_slice.SetMetaData("0008|0021", modification_date)
        image_slice.SetMetaData("0008|0008", "DERIVED\SECONDARY")
    # Each of the UID components is a number (cannot start with zero) and separated by a '.'
    # We create a unique series ID using the date and time.
        image_slice.SetMetaData("0020|000e", "1.2.826.0.1.3680043.2.1125."+modification_date+".1"+modification_time)# Each of the UID components is a number (cannot start with zero) and separated by a '.'
     # Write to the output directory and add the extension dcm if not there, to force writing is in DICOM format.
        writer.SetFileName(os.path.join(dest, os.path.basename(series_file_names[i])) + ('' if os.path.splitext(series_file_names[i])[1] == '.dcm' else '.dcm'))
        writer.Execute(image_slice)

质控

生成了相应的Dicom文件后，通常需要查验是否和原始的dicom图片是否一样，这时候用到我的第一篇博客里提到的，simpleITK读入image的四个元数据：Origin，Spacing，Size，Direction。
（具体请到第一篇博客里查看）

import SimpleITK as sitk
import numpy as np
filepath = '/'
image = sitk.ReadImage(filepath)#filepath是欲读入的一张图片
print('origin: ' + str(image.GetOrigin()))
print('size: ' + str(image.GetSize()))
print('spacing: ' + str(image.GetSpacing()))
print('direction: ' + str(image.GetDirection()))

后记

此博客参考了https://github.com/zivy/SimpleITK/blob/8e94451e4c0e90bcc6a1ffdd7bc3d56c81f58d80/Examples/DicomSeriesReadModifyWrite/DicomSeriesReadModifySeriesWrite.py。