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 = ):
        img_bitmap = Image.fromarray(img_array[frame_num])
       return img_bitmap
           

5. 采用CLAHE (Contrast Limited Adaptive Histogram Equalization)技術來優化圖像。

def limitedEqualize(img_array, limit = ):
        img_array_list = []
        for img in img_array:
              clahe = cv2.createCLAHE(clipLimit = limit, tileGridSize = (,))
              img_array_list.append(clahe.apply(img))
       img_array_limited_equalized = np.array(img_array_list)
       return img_array_limited_equalized
           

這一步對于整個圖像處理起到很重要的作用，可以根據不同的原始DICOM圖像的窗位和窗寬來進行動态調整，以達到最佳的識别效果。

最後應用OpenCV的Python架構cv2把每幀圖像組合在一起，生成通用視訊格式。

def writeVideo(img_array):
       frame_num, width, height = img_array.shape
       filename_output = filename.split('.')[] + '.avi'    
       video = cv2.VideoWriter(filename_output, -, , (width, height))    
       for img in img_array:
                  video.write(img)
   video.release()
           

VTK加載DICOM資料

import vtk  
from vtk.util import numpy_support  
import numpy  

PathDicom = "./dir_with_dicom_files/"  
reader = vtk.vtkDICOMImageReader()  
reader.SetDirectoryName(PathDicom)  
reader.Update()  

# Load dimensions using `GetDataExtent`  
_extent = reader.GetDataExtent()  
ConstPixelDims = [_extent[]-_extent[]+, _extent[]-_extent[]+, _extent[]-_extent[]+]  

# Load spacing values  
ConstPixelSpacing = reader.GetPixelSpacing()  

# Get the 'vtkImageData' object from the reader  
imageData = reader.GetOutput()  
# Get the 'vtkPointData' object from the 'vtkImageData' object  
pointData = imageData.GetPointData()  
# Ensure that only one array exists within the 'vtkPointData' object  
assert (pointData.GetNumberOfArrays()==)  
# Get the `vtkArray` (or whatever derived type) which is needed for the `numpy_support.vtk_to_numpy` function  
arrayData = pointData.GetArray()  

# Convert the `vtkArray` to a NumPy array  
ArrayDicom = numpy_support.vtk_to_numpy(arrayData)  
# Reshape the NumPy array to D using 'ConstPixelDims' as a 'shape'  
ArrayDicom = ArrayDicom.reshape(ConstPixelDims, order='F')  
           

PYDICOM加載DICOM資料：

可以在https://github.com/darcymason/pydicom的test裡面看怎麼用代碼。

import dicom  
import os  
import numpy  

PathDicom = "./dir_with_dicom_series/"  
lstFilesDCM = []  # create an empty list  
for dirName, subdirList, fileList in os.walk(PathDicom):  
    for filename in fileList:  
        if ".dcm" in filename.lower():  # check whether the file's DICOM  
            lstFilesDCM.append(os.path.join(dirName,filename))  

# Get ref file  
RefDs = dicom.read_file(lstFilesDCM[])  

# Load dimensions based on the number of rows, columns, and slices (along the Z axis)  
ConstPixelDims = (int(RefDs.Rows), int(RefDs.Columns), len(lstFilesDCM))  

# Load spacing values (in mm)  
ConstPixelSpacing = (float(RefDs.PixelSpacing[]), float(RefDs.PixelSpacing[]), float(RefDs.SliceThickness))  

# The array is sized based on 'ConstPixelDims'  
ArrayDicom = numpy.zeros(ConstPixelDims, dtype=RefDs.pixel_array.dtype)  

# loop through all the DICOM files  
for filenameDCM in lstFilesDCM:  
    # read the file  
    ds = dicom.read_file(filenameDCM)  
    # store the raw image data  
    ArrayDicom[:, :, lstFilesDCM.index(filenameDCM)] = ds.pixel_array  
           

轉換VTK built-in types to SimpleITK/ITK built-ins and vice-versa

import vtk  
import SimpleITK  

dctITKtoVTK = {SimpleITK.sitkInt8: vtk.VTK_TYPE_INT8,  
               SimpleITK.sitkInt16: vtk.VTK_TYPE_INT16,  
               SimpleITK.sitkInt32: vtk.VTK_TYPE_INT32,  
               SimpleITK.sitkInt64: vtk.VTK_TYPE_INT64,  
               SimpleITK.sitkUInt8: vtk.VTK_TYPE_UINT8,  
               SimpleITK.sitkUInt16: vtk.VTK_TYPE_UINT16,  
               SimpleITK.sitkUInt32: vtk.VTK_TYPE_UINT32,  
               SimpleITK.sitkUInt64: vtk.VTK_TYPE_UINT64,  
               SimpleITK.sitkFloat32: vtk.VTK_TYPE_FLOAT32,  
               SimpleITK.sitkFloat64: vtk.VTK_TYPE_FLOAT64}  
dctVTKtoITK = dict(zip(dctITKtoVTK.values(),   
                       dctITKtoVTK.keys()))  

def convertTypeITKtoVTK(typeITK):  
    if typeITK in dctITKtoVTK:  
        return dctITKtoVTK[typeITK]  
    else:  
        raise ValueError("Type not supported")  

def convertTypeVTKtoITK(typeVTK):  
    if typeVTK in dctVTKtoITK:  
        return dctVTKtoITK[typeVTK]  
    else:  
        raise ValueError("Type not supported")  
           

VTK-SimpleITK繪制資料

#!/usr/bin/python  

import SimpleITK as sitk  
import vtk  
import numpy as np  

from vtk.util.vtkConstants import *  

def numpy2VTK(img,spacing=[,,]):  
    # evolved from code from Stou S.,  
    # on http://www.siafoo.net/snippet/314  
    importer = vtk.vtkImageImport()  

    img_data = img.astype('uint8')  
    img_string = img_data.tostring() # type short  
    dim = img.shape  

    importer.CopyImportVoidPointer(img_string, len(img_string))  
    importer.SetDataScalarType(VTK_UNSIGNED_CHAR)  
    importer.SetNumberOfScalarComponents()  

    extent = importer.GetDataExtent()  
    importer.SetDataExtent(extent[], extent[] + dim[] - ,  
                           extent[], extent[] + dim[] - ,  
                           extent[], extent[] + dim[] - )  
    importer.SetWholeExtent(extent[], extent[] + dim[] - ,  
                            extent[], extent[] + dim[] - ,  
                            extent[], extent[] + dim[] - )  

    importer.SetDataSpacing( spacing[], spacing[], spacing[])  
    importer.SetDataOrigin( ,, )  

    return importer  

def volumeRender(img, tf=[],spacing=[,,]):  
    importer = numpy2VTK(img,spacing)  

    # Transfer Functions  
    opacity_tf = vtk.vtkPiecewiseFunction()  
    color_tf = vtk.vtkColorTransferFunction()  

    if len(tf) == :  
        tf.append([img.min(),,,,])  
        tf.append([img.max(),,,,])  

    for p in tf:  
        color_tf.AddRGBPoint(p[], p[], p[], p[])  
        opacity_tf.AddPoint(p[], p[])  

    # working on the GPU  
    # volMapper = vtk.vtkGPUVolumeRayCastMapper()  
    # volMapper.SetInputConnection(importer.GetOutputPort())  

    # # The property describes how the data will look  
    # volProperty =  vtk.vtkVolumeProperty()  
    # volProperty.SetColor(color_tf)  
    # volProperty.SetScalarOpacity(opacity_tf)  
    # volProperty.ShadeOn()  
    # volProperty.SetInterpolationTypeToLinear()  

    # working on the CPU  
    volMapper = vtk.vtkVolumeRayCastMapper()  
    compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()  
    compositeFunction.SetCompositeMethodToInterpolateFirst()  
    volMapper.SetVolumeRayCastFunction(compositeFunction)  
    volMapper.SetInputConnection(importer.GetOutputPort())  

    # The property describes how the data will look  
    volProperty =  vtk.vtkVolumeProperty()  
    volProperty.SetColor(color_tf)  
    volProperty.SetScalarOpacity(opacity_tf)  
    volProperty.ShadeOn()  
    volProperty.SetInterpolationTypeToLinear()  

    # Do the lines below speed things up?  
    # pix_diag = 5.0  
    # volMapper.SetSampleDistance(pix_diag / 5.0)      
    # volProperty.SetScalarOpacityUnitDistance(pix_diag)   


    vol = vtk.vtkVolume()  
    vol.SetMapper(volMapper)  
    vol.SetProperty(volProperty)  

    return [vol]  


def vtk_basic( actors ):  
    """ 
    Create a window, renderer, interactor, add the actors and start the thing 

    Parameters 
    ---------- 
    actors :  list of vtkActors 

    Returns 
    ------- 
    nothing 
    """       

    # create a rendering window and renderer  
    ren = vtk.vtkRenderer()  
    renWin = vtk.vtkRenderWindow()  
    renWin.AddRenderer(ren)  
    renWin.SetSize(,)  
    # ren.SetBackground( 1, 1, 1)  

    # create a renderwindowinteractor  
    iren = vtk.vtkRenderWindowInteractor()  
    iren.SetRenderWindow(renWin)  

    for a in actors:  
        # assign actor to the renderer  
        ren.AddActor(a )  

    # render  
    renWin.Render()  

    # enable user interface interactor  
    iren.Initialize()  
    iren.Start()  



#####  

filename = 'toto.nii.gz'  


img = sitk.ReadImage( filename ) # SimpleITK object  
data = sitk.GetArrayFromImage( img ) # numpy array  

from scipy.stats.mstats import mquantiles  
q = mquantiles(data.flatten(),[,])  
q[]=max(q[],)  
q[] = max(q[],)  
tf=[[,,,,],[q[],,,,],[q[],,,,],[data.max(),,,,]]  

actor_list = volumeRender(data, tf=tf, spacing=img.GetSpacing())  

vtk_basic(actor_list)  
           

下面一個不錯的軟體：

https://github.com/bastula/dicompyler

還有一個python的庫mudicom，https://github.com/neurosnap/mudicom

import mudicom  
mu = mudicom.load("mudicom/tests/dicoms/ex1.dcm")  

# returns array of data elements as dicts  
mu.read()  

# returns dict of errors and warnings for DICOM  
mu.validate()  

# basic anonymization  
mu.anonymize()  
# save anonymization  
mu.save_as("dicom.dcm")  

# creates image object  
img = mu.image # before v0.1.0 this was mu.image()  
# returns numpy array  
img.numpy # before v0.1.0 this was mu.numpy()  

# using Pillow, saves DICOM image  
img.save_as_pil("ex1.jpg")  
# using matplotlib, saves DICOM image  
img.save_as_plt("ex1_2.jpg")  
           

本文轉自：

http://mp.weixin.qq.com/s?__biz=MzAxOTk4NTIwMw==&mid=2247483968&idx=1&sn=2844930d81f3e1f45260338dae21d8ea&chksm=9c3fe41cab486d0aa9d03a6494865c0ffdbb0e878f4804227a73d2d403382432fe22e9c03b71&mpshare=1&scene=23&srcid=0122xHKRmyLclSvbojaS1ICX##

http://blog.csdn.net/langb2014/article/details/54667905#

http://www.jianshu.com/p/df64088e9b6b