Python 處理醫學影像學中的DICOM
阿新 • • 發佈：2019-02-01
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 cv22. 應用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, height3. 應用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 information4. 應用PIL來檢查影象是否被提取。
def showImage(img_array, frame_num = 0):
img_bitmap = Image.fromarray(img_array[frame_num])
return img_bitmap5. 採用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影象的窗位和窗寬來進行動態調整，以達到最佳的識別效果。
最後應用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資料
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[1]-_extent[0]+1, _extent[3]-_extent[2]+1, _extent[5]-_extent[4]+1]
# 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()==1)
# Get the `vtkArray` (or whatever derived type) which is needed for the `numpy_support.vtk_to_numpy` function
arrayData = pointData.GetArray(0)
# Convert the `vtkArray` to a NumPy array
ArrayDicom = numpy_support.vtk_to_numpy(arrayData)
# Reshape the NumPy array to 3D using 'ConstPixelDims' as a 'shape'
ArrayDicom = ArrayDicom.reshape(ConstPixelDims, order='F') 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[0])
# 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[0]), float(RefDs.PixelSpacing[1]), 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=[1.0,1.0,1.0]):
# 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(1)
extent = importer.GetDataExtent()
importer.SetDataExtent(extent[0], extent[0] + dim[2] - 1,
extent[2], extent[2] + dim[1] - 1,
extent[4], extent[4] + dim[0] - 1)
importer.SetWholeExtent(extent[0], extent[0] + dim[2] - 1,
extent[2], extent[2] + dim[1] - 1,
extent[4], extent[4] + dim[0] - 1)
importer.SetDataSpacing( spacing[0], spacing[1], spacing[2])
importer.SetDataOrigin( 0,0,0 )
return importer
def volumeRender(img, tf=[],spacing=[1.0,1.0,1.0]):
importer = numpy2VTK(img,spacing)
# Transfer Functions
opacity_tf = vtk.vtkPiecewiseFunction()
color_tf = vtk.vtkColorTransferFunction()
if len(tf) == 0:
tf.append([img.min(),0,0,0,0])
tf.append([img.max(),1,1,1,1])
for p in tf:
color_tf.AddRGBPoint(p[0], p[1], p[2], p[3])
opacity_tf.AddPoint(p[0], p[4])
# 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(600,600)
# 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(),[0.7,0.98])
q[0]=max(q[0],1)
q[1] = max(q[1],1)
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]]
actor_list = volumeRender(data, tf=tf, spacing=img.GetSpacing())
vtk_basic(actor_list) 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") 本文轉自：