1. 前言
为了能够在Labelme上对Dicom图像进行编辑,这里对python环境下Dicom文件的读取进行了研究。在Dicom图像中CT的窗宽窗位是一个很重要的概念,但是找了半天在pydicom中没有相关设置函数,这里跟DCMTK还不一样。但是可以根据两个tag得到CT图像的CT值,那就是(0028|1052):rescale intercept和(0028|1053):rescale slope。则按照下面的算子得到CT图像,进而就可以调整窗宽窗位了
Hu = pixel * slope + intercept
至于那个部位的窗宽窗位是多少各位看官就可以自行百度了。
2. 代码实现
# -*- coding=utf-8 -*-import matplotlib.pyplot as pltimport pydicomimport pydicom.uidimport sysimport PIL.Image as Imagefrom PyQt5 import QtGuiimport os
step1:读取Dicom图像数据与得到CT值图像(CT图)
have_numpy = Truetry:import numpyexcept ImportError:have_numpy = Falseraisesys_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):"""Returns-------boolTrue if this pixel data handler might support this transfer syntax.False to prevent any attempt to try to use this handlerto decode the given transfer syntax"""return (dicom_dataset.file_meta.TransferSyntaxUID inNumpySupportedTransferSyntaxes)def needs_to_convert_to_RGB(dicom_dataset):return Falsedef 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------TypeErrorIf there is no Pixel Data or not a supported data type.ImportErrorIf NumPy isn't foundNotImplementedErrorif the transfer syntax is not supportedAttributeErrorif the decoded amount of data does not match the expected amountReturns-------numpy.ndarrayThe contents of the Pixel Data element (7FE0,0010) as an ndarray."""if (dicom_dataset.file_meta.TransferSyntaxUID not inNumpySupportedTransferSyntaxes):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.")# 设置窗宽窗位#dicom_dataset.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 32if dicom_dataset.BitsAllocated == 1:# single bits are used for representation of binary dataformat_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.PixelDataif 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.unpackbitsraise NotImplementedError('Cannot handle BitsAllocated == 1 on this platform')else:pixel_array = numpy.frombuffer(pixel_bytearray, dtype=numpy_dtype)length_of_pixel_array = pixel_array.nbytesexpected_length = dicom_dataset.Rows * dicom_dataset.Columnsif ('NumberOfFrames' in dicom_dataset anddicom_dataset.NumberOfFrames > 1):expected_length *= dicom_dataset.NumberOfFramesif ('SamplesPerPixel' in dicom_dataset anddicom_dataset.SamplesPerPixel > 1):expected_length *= dicom_dataset.SamplesPerPixelif dicom_dataset.BitsAllocated > 8:expected_length *= (dicom_dataset.BitsAllocated // 8)padded_length = expected_lengthif expected_length & 1:padded_length += 1if 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
step2:对于CT图像设置窗宽窗位
# 调整CT图像的窗宽窗位def setDicomWinWidthWinCenter(img_data, winwidth, wincenter, rows, cols):img_temp = img_dataimg_temp.flags.writeable = Truemin = (2 * wincenter - winwidth) / 2.0 + 0.5max = (2 * wincenter + winwidth) / 2.0 + 0.5dFactor = 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 < 0img_temp[min_index] = 0max_index = img_temp > 255img_temp[max_index] = 255return img_temp
step3:获取Dicom中的tag信息
第一种方式:
# 加载Dicom图片中的Tag信息def loadFileInformation(filename):information = {}ds = pydicom.read_file(filename)information['PatientID'] = ds.PatientIDinformation['PatientName'] = ds.PatientNameinformation['PatientBirthDate'] = ds.PatientBirthDateinformation['PatientSex'] = ds.PatientSexinformation['StudyID'] = ds.StudyIDinformation['StudyDate'] = ds.StudyDateinformation['StudyTime'] = ds.StudyTimeinformation['InstitutionName'] = ds.InstitutionNameinformation['Manufacturer'] = ds.Manufacturerprint(dir(ds))print(type(information))return information
第二种方式
dcm = pydicom.dcmread(fileanme) # 加载Dicom数据print(dcm[0x0008, 0x0060])>>(0008, 0060) Modality CS: 'MR'print(dcm[0x0008, 0x0060].VR)>>CSprint(dcm[0x0008, 0x0060].value)>>MR
step4:Dicom图像数据转换为PIL.Image
dcm = pydicom.dcmread(fileanme) # 加载Dicom数据dcm_img = Image.fromarray(img_data) # 将Numpy转换为PIL.Imagedcm_img = dcm_img.convert('L')# 保存为jpg文件,用作后面的生成label用dcm_img.save('temp.jpg')# 显示图像dcm_img.show()
3. 结果展示
调整了窗宽窗位的脑部CT图像: