Python处理医学影像学中的DICOM

DICOMDICOM（Digital Imaging and Communications in Medicine）即医学数字成像和通信，是医学图像和相关信息的国际标准（ISO 12052）。它定义了质量能满足临床需要的可用于数据交换的医学图像格式，可用于处理、存储、打印和传输医学影像信息。DICOM可以便捷地交换于两个满足DICOM格式协议的工作站之间。目前该协议标准不仅广泛应用于大型医院，而且已成为小型诊所和牙科诊所医生办公室的标准影像阅读格式。

DICOM被广泛应用于放射医疗、心血管成像以及放射诊疗诊断设备（X射线，CT，核磁共振，超声等），并且在眼科和牙科等其它医学领域得到越来越深入广泛的应用。在数以万计的在用医学成像设备中，DICOM是部署最为广泛的医疗信息标准之一。当前大约有百亿级符合DICOM标准的医学图像用于临床使用。

目前，越来越多的DICOM应用程序和分析软件被运用于临床医学，促使越来越多的编程语言开发出支持DICOM API的框架。今天就让我来介绍一下Python语言下支持的DICOM模块，以及如何完成基本DICOM信息分析和处理的编程方法。

---

Pydicom

Pydicom是一个处理DICOM文件的纯Python软件包。它可以通过非常容易的“Pythonic”的方式来提取和修改DICOM数据，修改后的数据还会借此生成新的DICOM文件。作为一个纯Python包，Pydicom可以在Python解释器下任何平台运行，除了必须预先安装Numpy模块外，几乎无需其它任何配置要求。其局限性之一是无法处理压缩像素图像（如JPEG），其次是无法处理分帧动画图像（如造影电影）。

SimpleITK

Insight Segmentation and Registration Toolkit (ITK)是一个开源、跨平台的框架，可以提供给开发者增强功能的图像分析和处理套件。其中最为著名的就是SimpleITK，是一个简化版的、构建于ITK最顶层的模块。SimpleITK旨在易化图像处理流程和方法。

PIL

Python Image Library (PIL) 是在Python环境下不可缺少的图像处理模块，支持多种格式，并提供强大的图形与图像处理功能，而且API却非常简单易用。

OpenCV

OpenCV的全称是：Open Source Computer Vision Library。OpenCV是一个基于BSD许可（开源）发行的跨平台计算机视觉库，可以运行在Linux、Windows和Mac OS操作系统上。它轻量级而且高效——由一系列 C 函数和少量 C++ 类构成，同时提供了Python、Ruby、MATLAB等语言的接口，实现了图像处理和计算机视觉方面的很多通用算法。

---

下面Python代码来演示如何编程处理心血管冠脉造影DICOM图像信息。

1. 导入主要框架：SimpleITK、pydicom、PIL、cv2和numpy

import SimpleITK as sitk

from PIL import Image

import pydicom

import numpy as np

import cv2

2. 应用SimpleITK框架来读取DICOM文件的矩阵信息。如果DICOM图像是三维螺旋CT图像，则帧参数则代表CT扫描层数；而如果是造影动态电影图像，则帧参数就是15帧/秒的电影图像帧数。

def loadFile(filename):

       ds = sitk.ReadImage(filename)

       img_array = sitk.GetArrayFromImage(ds)

       frame_num, width, height = img_array.shape

      return img_array, frame_num, width, height

3. 应用pydicom来提取患者信息。

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    

      return information

4. 应用PIL来检查图像是否被提取。

def showImage(img_array, frame_num = 0):

        img_bitmap = Image.fromarray(img_array[frame_num])

       return img_bitmap

5. 采用CLAHE (Contrast Limited Adaptive Histogram Equalization)技术来优化图像。

def limitedEqualize(img_array, limit = 4.0):

        img_array_list = []

        for img in img_array:

              clahe = cv2.createCLAHE(clipLimit = limit, tileGridSize = (8,8))

              img_array_list.append(clahe.apply(img))

       img_array_limited_equalized = np.array(img_array_list)

       return img_array_limited_equalized

这一步对于整个图像处理起到很重要的作用，可以根据不同的原始DICOM图像的窗位和窗宽来进行动态调整，以达到最佳的识别效果。

原始图像：

经过自动窗位窗宽调节，生成：

再经过CLAHE优化，则生成：

最后应用OpenCV的Python框架cv2把每帧图像组合在一起，生成通用视频格式。

def writeVideo(img_array):

       frame_num, width, height = img_array.shape

       filename_output = filename.split('.')[0] + '.avi'    

       video = cv2.VideoWriter(filename_output, -1, 16, (width, height))    

       for img in img_array:

                  video.write(img)

   video.release()

由于我没有数据，微信图片无法复制过来，只能看看代码了，详细图片请看参考。
======================================================================================================

VTK加载DICOM数据

[cpp] view plain copy

1. import vtk  
2. from vtk.util import numpy_support  
3. import numpy  
4.   
5. PathDicom = "./dir_with_dicom_files/"  
6. reader = vtk.vtkDICOMImageReader()  
7. reader.SetDirectoryName(PathDicom)  
8. reader.Update()  
9.   
10. # Load dimensions using `GetDataExtent`  
11. _extent = reader.GetDataExtent()  
12. ConstPixelDims = [_extent[1]-_extent[0]+1, _extent[3]-_extent[2]+1, _extent[5]-_extent[4]+1]  
13.   
14. # Load spacing values  
15. ConstPixelSpacing = reader.GetPixelSpacing()  
16.   
17. # Get the 'vtkImageData' object from the reader  
18. imageData = reader.GetOutput()  
19. # Get the 'vtkPointData' object from the 'vtkImageData' object  
20. pointData = imageData.GetPointData()  
21. # Ensure that only one array exists within the 'vtkPointData' object  
22. assert (pointData.GetNumberOfArrays()==1)  
23. # Get the `vtkArray` (or whatever derived type) which is needed for the `numpy_support.vtk_to_numpy` function  
24. arrayData = pointData.GetArray(0)  
25.   
26. # Convert the `vtkArray` to a NumPy array  
27. ArrayDicom = numpy_support.vtk_to_numpy(arrayData)  
28. # Reshape the NumPy array to 3D using 'ConstPixelDims' as a 'shape'  
29. ArrayDicom = ArrayDicom.reshape(ConstPixelDims, order='F')  

PYDICOM加载DICOM数据：

可以在https://github.com/darcymason/pydicom的test里面看怎么用代码。

[python] view plain copy

1. import dicom  
2. import os  
3. import numpy  
4.   
5. PathDicom = "./dir_with_dicom_series/"  
6. lstFilesDCM = []  # create an empty list  
7. for dirName, subdirList, fileList in os.walk(PathDicom):  
8.     for filename in fileList:  
9.         if ".dcm" in filename.lower():  # check whether the file's DICOM  
10.             lstFilesDCM.append(os.path.join(dirName,filename))  
11.               
12. # Get ref file  
13. RefDs = dicom.read_file(lstFilesDCM[0])  
14.   
15. # Load dimensions based on the number of rows, columns, and slices (along the Z axis)  
16. ConstPixelDims = (int(RefDs.Rows), int(RefDs.Columns), len(lstFilesDCM))  
17.   
18. # Load spacing values (in mm)  
19. ConstPixelSpacing = (float(RefDs.PixelSpacing[0]), float(RefDs.PixelSpacing[1]), float(RefDs.SliceThickness))  
20.   
21. # The array is sized based on 'ConstPixelDims'  
22. ArrayDicom = numpy.zeros(ConstPixelDims, dtype=RefDs.pixel_array.dtype)  
23.   
24. # loop through all the DICOM files  
25. for filenameDCM in lstFilesDCM:  
26.     # read the file  
27.     ds = dicom.read_file(filenameDCM)  
28.     # store the raw image data  
29.     ArrayDicom[:, :, lstFilesDCM.index(filenameDCM)] = ds.pixel_array  

转换VTK built-in types to SimpleITK/ITK built-ins and vice-versa

[python] view plain copy

1. import vtk  
2. import SimpleITK  
3.   
4. dctITKtoVTK = {SimpleITK.sitkInt8: vtk.VTK_TYPE_INT8,  
5.                SimpleITK.sitkInt16: vtk.VTK_TYPE_INT16,  
6.                SimpleITK.sitkInt32: vtk.VTK_TYPE_INT32,  
7.                SimpleITK.sitkInt64: vtk.VTK_TYPE_INT64,  
8.                SimpleITK.sitkUInt8: vtk.VTK_TYPE_UINT8,  
9.                SimpleITK.sitkUInt16: vtk.VTK_TYPE_UINT16,  
10.                SimpleITK.sitkUInt32: vtk.VTK_TYPE_UINT32,  
11.                SimpleITK.sitkUInt64: vtk.VTK_TYPE_UINT64,  
12.                SimpleITK.sitkFloat32: vtk.VTK_TYPE_FLOAT32,  
13.                SimpleITK.sitkFloat64: vtk.VTK_TYPE_FLOAT64}  
14. dctVTKtoITK = dict(zip(dctITKtoVTK.values(),   
15.                        dctITKtoVTK.keys()))  
16.                          
17. def convertTypeITKtoVTK(typeITK):  
18.     if typeITK in dctITKtoVTK:  
19.         return dctITKtoVTK[typeITK]  
20.     else:  
21.         raise ValueError("Type not supported")  
22.   
23. def convertTypeVTKtoITK(typeVTK):  
24.     if typeVTK in dctVTKtoITK:  
25.         return dctVTKtoITK[typeVTK]  
26.     else:  
27.         raise ValueError("Type not supported")  

VTK-SimpleITK绘制数据

[python] view plain copy

1. #!/usr/bin/python  
2.   
3. import SimpleITK as sitk  
4. import vtk  
5. import numpy as np  
6.   
7. from vtk.util.vtkConstants import *  
8.   
9. def numpy2VTK(img,spacing=[1.0,1.0,1.0]):  
10.     # evolved from code from Stou S.,  
11.     # on http://www.siafoo.net/snippet/314  
12.     importer = vtk.vtkImageImport()  
13.       
14.     img_data = img.astype('uint8')  
15.     img_string = img_data.tostring() # type short  
16.     dim = img.shape  
17.       
18.     importer.CopyImportVoidPointer(img_string, len(img_string))  
19.     importer.SetDataScalarType(VTK_UNSIGNED_CHAR)  
20.     importer.SetNumberOfScalarComponents(1)  
21.       
22.     extent = importer.GetDataExtent()  
23.     importer.SetDataExtent(extent[0], extent[0] + dim[2] - 1,  
24.                            extent[2], extent[2] + dim[1] - 1,  
25.                            extent[4], extent[4] + dim[0] - 1)  
26.     importer.SetWholeExtent(extent[0], extent[0] + dim[2] - 1,  
27.                             extent[2], extent[2] + dim[1] - 1,  
28.                             extent[4], extent[4] + dim[0] - 1)  
29.   
30.     importer.SetDataSpacing( spacing[0], spacing[1], spacing[2])  
31.     importer.SetDataOrigin( 0,0,0 )  
32.   
33.     return importer  
34.   
35. def volumeRender(img, tf=[],spacing=[1.0,1.0,1.0]):  
36.     importer = numpy2VTK(img,spacing)  
37.   
38.     # Transfer Functions  
39.     opacity_tf = vtk.vtkPiecewiseFunction()  
40.     color_tf = vtk.vtkColorTransferFunction()  
41.   
42.     if len(tf) == 0:  
43.         tf.append([img.min(),0,0,0,0])  
44.         tf.append([img.max(),1,1,1,1])  
45.   
46.     for p in tf:  
47.         color_tf.AddRGBPoint(p[0], p[1], p[2], p[3])  
48.         opacity_tf.AddPoint(p[0], p[4])  
49.   
50.     # working on the GPU  
51.     # volMapper = vtk.vtkGPUVolumeRayCastMapper()  
52.     # volMapper.SetInputConnection(importer.GetOutputPort())  
53.   
54.     # # The property describes how the data will look  
55.     # volProperty =  vtk.vtkVolumeProperty()  
56.     # volProperty.SetColor(color_tf)  
57.     # volProperty.SetScalarOpacity(opacity_tf)  
58.     # volProperty.ShadeOn()  
59.     # volProperty.SetInterpolationTypeToLinear()  
60.   
61.     # working on the CPU  
62.     volMapper = vtk.vtkVolumeRayCastMapper()  
63.     compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()  
64.     compositeFunction.SetCompositeMethodToInterpolateFirst()  
65.     volMapper.SetVolumeRayCastFunction(compositeFunction)  
66.     volMapper.SetInputConnection(importer.GetOutputPort())  
67.   
68.     # The property describes how the data will look  
69.     volProperty =  vtk.vtkVolumeProperty()  
70.     volProperty.SetColor(color_tf)  
71.     volProperty.SetScalarOpacity(opacity_tf)  
72.     volProperty.ShadeOn()  
73.     volProperty.SetInterpolationTypeToLinear()  
74.       
75.     # Do the lines below speed things up?  
76.     # pix_diag = 5.0  
77.     # volMapper.SetSampleDistance(pix_diag / 5.0)      
78.     # volProperty.SetScalarOpacityUnitDistance(pix_diag)   
79.       
80.   
81.     vol = vtk.vtkVolume()  
82.     vol.SetMapper(volMapper)  
83.     vol.SetProperty(volProperty)  
84.       
85.     return [vol]  
86.   
87.   
88. def vtk_basic( actors ):  
89.     """ 
90.     Create a window, renderer, interactor, add the actors and start the thing 
91.      
92.     Parameters 
93.     ---------- 
94.     actors :  list of vtkActors 
95.      
96.     Returns 
97.     ------- 
98.     nothing 
99.     """       
100.       
101.     # create a rendering window and renderer  
102.     ren = vtk.vtkRenderer()  
103.     renWin = vtk.vtkRenderWindow()  
104.     renWin.AddRenderer(ren)  
105.     renWin.SetSize(600,600)  
106.     # ren.SetBackground( 1, 1, 1)  
107.    
108.     # create a renderwindowinteractor  
109.     iren = vtk.vtkRenderWindowInteractor()  
110.     iren.SetRenderWindow(renWin)  
111.   
112.     for a in actors:  
113.         # assign actor to the renderer  
114.         ren.AddActor(a )  
115.       
116.     # render  
117.     renWin.Render()  
118.      
119.     # enable user interface interactor  
120.     iren.Initialize()  
121.     iren.Start()  
122.   
123.           
124.   
125. #####  
126.   
127. filename = 'toto.nii.gz'  
128.   
129.   
130. img = sitk.ReadImage( filename ) # SimpleITK object  
131. data = sitk.GetArrayFromImage( img ) # numpy array  
132.   
133. from scipy.stats.mstats import mquantiles  
134. q = mquantiles(data.flatten(),[0.7,0.98])  
135. q[0]=max(q[0],1)  
136. q[1] = max(q[1],1)  
137. tf=[[0,0,0,0,0],[q[0],0,0,0,0],[q[1],1,1,1,0.5],[data.max(),1,1,1,1]]  
138.   
139. actor_list = volumeRender(data, tf=tf, spacing=img.GetSpacing())  
140.   
141. vtk_basic(actor_list)  

下面一个不错的软件：

https://github.com/bastula/dicompyler
还有一个python的库mudicom，https://github.com/neurosnap/mudicom

[python] view plain copy

1. import mudicom  
2. mu = mudicom.load("mudicom/tests/dicoms/ex1.dcm")  
3.   
4. # returns array of data elements as dicts  
5. mu.read()  
6.   
7. # returns dict of errors and warnings for DICOM  
8. mu.validate()  
9.   
10. # basic anonymization  
11. mu.anonymize()  
12. # save anonymization  
13. mu.save_as("dicom.dcm")  
14.   
15. # creates image object  
16. img = mu.image # before v0.1.0 this was mu.image()  
17. # returns numpy array  
18. img.numpy # before v0.1.0 this was mu.numpy()  
19.   
20. # using Pillow, saves DICOM image  
21. img.save_as_pil("ex1.jpg")  
22. # using matplotlib, saves DICOM image  
23. img.save_as_plt("ex1_2.jpg")  

参考：

微信