python医学影像2Ddicom文件转成3Dnii文件(保留原始dicom信息)

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了python医学影像2Ddicom文件转成3Dnii文件(保留原始dicom信息)相关的知识,希望对你有一定的参考价值。

参考技术A

导入需要的包

定义函数

定义待读取文件夹路径和nii文件存储路径

图像识别 | 使用Python对医学Dicom文件的预处理(含代码)


1


前沿

在处理医学图像时,常常会遇到以Dicom格式保存的医学图像,如CT、MRI等。Dicom文件是需要专门的软件或者通过编程,应用相应的库进行处理。为了能够更好地服务下游任务,例如分割或检测腹腔CT图像中某个病灶组织,需要先将Dicom图像进行读取,脱敏,调窗等步骤,以便于后续的编辑。本文使用python对医学Dicom文件进行相应的处理,相比于封装好的软件,笔者认为自己动手的可操作性更强。


2


目录

1 导入相应的包

2 读取Dicom图像数据

3 设置CT图像的窗宽和窗位

4 获取Dicom图像的tag信息

5 结果保存及可视化


3


导入相应的包

# load necessary packages
import matplotlib.pyplot as plt
import pydicom.uid
import sys
from PyQt5 import QtGui
import os
import pydicom
import glob
from PIL import *
import matplotlib.pyplot as plt
from pylab import *
from tkinter.filedialog import *
import PIL.Image as Image

其核心是使用了python中的pydicom库来处理dicom文件。


4


读取Dicom图像数据

have_numpy = True

try:
   import numpy
except ImportError:
   have_numpy = False
   raise

sys_is_little_endian = (sys.byteorder == 'little')

NumpySupportedTransferSyntaxes = [
   pydicom.uid.ExplicitVRLittleEndian,
   pydicom.uid.ImplicitVRLittleEndian,
   pydicom.uid.DeflatedExplicitVRLittleEndian,
   pydicom.uid.ExplicitVRBigEndian,
]

# 支持"传输"语法
def supports_transfer_syntax(dicom_dataset):

   return (dicom_dataset.file_meta.TransferSyntaxUID in
           NumpySupportedTransferSyntaxes)

def needs_to_convert_to_RGB(dicom_dataset):
   return False


def should_change_PhotometricInterpretation_to_RGB(dicom_dataset):
   return False


# 加载 Dicom图像
def get_pixeldata(dicom_dataset):
   """If NumPy is available, return an ndarray of the Pixel Data.
   Raises
   ------
   TypeError
       If there is no Pixel Data or not a supported data type.
   ImportError
       If NumPy isn't found
   NotImplementedError
       if the transfer syntax is not supported
   AttributeError
       if the decoded amount of data does not match the expected amount
   Returns
   -------
   numpy.ndarray
      The contents of the Pixel Data element (7FE0,0010) as an ndarray.
   """
   if (dicom_dataset.file_meta.TransferSyntaxUID not in
           NumpySupportedTransferSyntaxes):
       raise NotImplementedError("Pixel Data is compressed in a "
                                 "format pydicom does not yet handle. "
                                 "Cannot return array. Pydicom might "
                                 "be able to convert the pixel data "
                                 "using GDCM if it is installed.")


   if not have_numpy:
       msg = ("The Numpy package is required to use pixel_array, and "
              "numpy could not be imported.")
       raise ImportError(msg)
   if 'PixelData' not in dicom_dataset:
       raise TypeError("No pixel data found in this dataset.")

   # Make NumPy format code, e.g. "uint16", "int32" etc
   # from two pieces of info:
   # dicom_dataset.PixelRepresentation -- 0 for unsigned, 1 for signed;
   # dicom_dataset.BitsAllocated -- 8, 16, or 32
   if dicom_dataset.BitsAllocated == 1:
       # single bits are used for representation of binary data
       format_str = 'uint8'
   elif dicom_dataset.PixelRepresentation == 0:
       format_str = 'uint{}'.format(dicom_dataset.BitsAllocated)
   elif dicom_dataset.PixelRepresentation == 1:
       format_str = 'int{}'.format(dicom_dataset.BitsAllocated)
   else:
       format_str = 'bad_pixel_representation'
   try:
       numpy_dtype = numpy.dtype(format_str)
   except TypeError:
       msg = ("Data type not understood by NumPy: "
              "format='{}', PixelRepresentation={}, "
              "BitsAllocated={}".format(
                  format_str,
                  dicom_dataset.PixelRepresentation,
                  dicom_dataset.BitsAllocated))
       raise TypeError(msg)

   if dicom_dataset.is_little_endian != sys_is_little_endian:
       numpy_dtype = numpy_dtype.newbyteorder('S')

   pixel_bytearray = dicom_dataset.PixelData

   if dicom_dataset.BitsAllocated == 1:
       # if single bits are used for binary representation, a uint8 array
       # has to be converted to a binary-valued array (that is 8 times bigger)
       try:
           pixel_array = numpy.unpackbits(
               numpy.frombuffer(pixel_bytearray, dtype='uint8'))
       except NotImplementedError:
           # PyPy2 does not implement numpy.unpackbits
           raise NotImplementedError(
               'Cannot handle BitsAllocated == 1 on this platform')
   else:
       pixel_array = numpy.frombuffer(pixel_bytearray, dtype=numpy_dtype)
   length_of_pixel_array = pixel_array.nbytes
   expected_length = dicom_dataset.Rows * dicom_dataset.Columns
   if ('NumberOfFrames' in dicom_dataset and
           dicom_dataset.NumberOfFrames > 1):
       expected_length *= dicom_dataset.NumberOfFrames
   if ('SamplesPerPixel' in dicom_dataset and
           dicom_dataset.SamplesPerPixel > 1):
       expected_length *= dicom_dataset.SamplesPerPixel
   if dicom_dataset.BitsAllocated > 8:
       expected_length *= (dicom_dataset.BitsAllocated // 8)
   padded_length = expected_length
   if expected_length & 1:
       padded_length += 1
   if length_of_pixel_array != padded_length:
       raise AttributeError(
           "Amount of pixel data %d does not "
           "match the expected data %d" %
           (length_of_pixel_array, padded_length))
   if expected_length != padded_length:
       pixel_array = pixel_array[:expected_length]
   if should_change_PhotometricInterpretation_to_RGB(dicom_dataset):
       dicom_dataset.PhotometricInterpretation = "RGB"
   if dicom_dataset.Modality.lower().find('ct') >= 0:  # CT图像需要得到其CT值图像
       pixel_array = pixel_array * dicom_dataset.RescaleSlope + dicom_dataset.RescaleIntercept  # 获得图像的CT值
   pixel_array = pixel_array.reshape(dicom_dataset.Rows, dicom_dataset.Columns*dicom_dataset.SamplesPerPixel)
   return pixel_array, dicom_dataset.Rows, dicom_dataset.Columns

读取到dicom文件中的数据后,实质上是几个图像矩阵,这个过程同时也处理了“脱敏”问题。


5


设置CT图像的窗宽和窗位

def setDicomWinWidthWinCenter(img_data, winwidth, wincenter, rows, cols):
   img_temp = img_data
   img_temp.flags.writeable = True
   min = (2 * wincenter - winwidth) / 2.0 + 0.5
   max = (2 * wincenter + winwidth) / 2.0 + 0.5
   dFactor = 255.0 / (max - min)

   for i in numpy.arange(rows):
       for j in numpy.arange(cols):
           img_temp[i, j] = int((img_temp[i, j]-min)*dFactor)

   min_index = img_temp < min
   img_temp[min_index] = 0
   max_index = img_temp > max
   img_temp[max_index] = 255

   return img_temp

该函数的输入变量winwidth和wincenter即为需要设置的窗宽和窗位,这两个值根据研究的问题(不同的组织器官对应不同的窗宽和窗位,有时候也要根据图像效果进行一定的调整)调整不同的值。网上有很多关于相关的窗位和窗宽对应值,这里给出一些参考资料,如果遇到不确定的,最好借鉴查阅相应领域的论文。

点击链接:https://radiopaedia.org/articles/windowing-ct


6


获取Dicom图像的tag信息

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
   print(dir(ds))
   print(type(information))
   return information

这个步骤视需求而定,如果不需要查看dicom文件的具体tag信息,此步骤可以跳过。


7


结果保存及可视化

可以单张保存,或者批量处理。

读取单张dicom文件

def main_single():
   dcm = dicom.read_file('81228816') # load dicom_file
   # 得到 CT 值,图像的 长, 宽
   pixel_array, dcm.Rows, dcm.Columns = get_pixeldata(dcm)
   # 调整窗位、窗宽
   img_data = pixel_array
   winwidth = 500
   wincenter = 50
   rows = dcm.Rows
   cols = dcm.Columns
   dcm_temp = setDicomWinWidthWinCenter(img_data, winwidth, wincenter, rows, cols)
   # 可视化
   dcm_img = Image.fromarray(dcm_temp)  # 将Numpy转换为PIL.Image
   dcm_img = dcm_img.convert('L')
   # plt.imshow(img, cmap=plt.cm.bone)
   # 保存为jpg文件,用作后面的生成label用
   dcm_img.save('../output/temp.jpg')
   # 显示图像
   dcm_img.show()

同时读取一个文件夹中的 dicom 文件,并处理保存 (写成循环即可)

def main_mulit(path):

   names = os.listdir(path)  # 读取文件夹中的所有文件名
   for i in range(len(names)):
       dicom_name = path+names[i]
       dcm = pydicom.read_file(dicom_name)  # 读取 dicom 文件
       pixel_array, dcm.Rows, dcm.Columns = get_pixeldata(dcm)  # 得到 dicom文件的 CT 值
       img_data = pixel_array
       winwidth = 500
       wincenter = 50
       rows = dcm.Rows
       cols = dcm.Columns
       dcm_temp = setDicomWinWidthWinCenter(img_data, winwidth, wincenter, rows, cols) # 调整窗位、窗宽
       #  可视化
       dcm_img = Image.fromarray(dcm_temp)  # 将 Numpy转换为 PIL.Image
       dcm_img = dcm_img.convert('L')
       # 批量保存
       dcm_img.save("C:/output/%s_%s.png" % (path1, names[i]))

注意,以上都写成函数的形式,运行时需要调用,并注意文件路径的修改。

单张dicom文件处理

main_single()

批量处理

main_mulit(path)


欢迎关注抖音:帕帕科技猫和B站,不时更新数码科技前沿技术。







以上是关于python医学影像2Ddicom文件转成3Dnii文件(保留原始dicom信息)的主要内容,如果未能解决你的问题,请参考以下文章

在ENVI里处理完影像,计算完NDVI,转成tif在arcgis里打开值域不是1到-1

请问在ARCGIS中如何加入影像图?

python包NiBabel对医学影像文件格式进行读写:python包NiBabel简介集示例

ArcGIS中的一些影像格式:.ige .img .rde .rrd

python包NiBabel对医学影像文件格式进行读写并可视化实战:查看和显示.nii.gz.nii文件

怎么把批处理文件导入python