在医学影像的数据存储领域，是存在一定的行业标准的。X光、CT机器等等医疗器械等生产企业，会依据行业标准，对采集的数据进行规范化的存储。
这里面就包括了大名鼎鼎的DICOM 3.0协议，上述的摄影形式大部分也都是以这种形式进行存储和传播的。
但是呢，在医学领域进行数据处理的时候，经常会遇到除DICOM外其他的数据形式，比如常见的 .nii、.npy、.npz、.mhd、.raw等等格式，千奇百怪，五花八门。
还有就是常见的日常图像存储的图像格式，比如 .jpeg、.png 、.bmp等等格式，这种格式下我们可以轻而易举的用电脑中的图像软件打开，这也是最方便我们查看的图像存储格式。
本系列就留作一个笔记本，主要记录上述医学影像的数据存储形式，以及数据间的转换，希望对初次接触的你，有所帮助。本文就对dicom文件，进行数据读取、存储等等的操作，都会汇总到这篇博客下。
一、导言
DICOM（Digital Imaging and Communications in Medicine）数据格式， 数据存储形式是一个存储信息比较丰富的文件格式，也是医学领域专业的数据格式。存储的不仅仅有图像数组信息，还有其他相关的信息。比如拍摄日期、拍摄条件、患者信息等等。
二、dicom文件读取
dicom文件，多以dcm作为文件的后缀形式。如果你拿到的数据，明知道是dicom标准存储的文件，但是确没有.dcm后缀，多半是导出时候给忘记加上这个字段了。
此时，你只需要将.dcm的后缀，加到对应的文件后面即可。为了验证这个文件是不是dicom文件，可以采用microDicom软件进行打开。可以一次双击打开一个文件，也可以鼠标右键存储很多dcm文件的文件夹，点击打开，就可以一次性打开整个文件夹的图像了，这个场景对于查看CT图像很有用。
与此同时，microDicom软件是公开的，可以直接去官网下载安装就可以了。软件内容很丰富，可以查看图像，导出图像，调整窗位、窗宽帮助我们查看图像都可以。这里我就不过多的赘述了，网上资料较多。
下面就讲重点，就是dicom如何用Python代码，批量的读取和处理，转成方便我们能查看的png文件呢？
2.1、pydicom库帮助打开dicom文件
Pydicom是一个用于处理DICOM格式文件的Python包。支持DICOM格式的读取：可以将dicom文件读入python结构，同时支持修改后的数据集可以再次写入DICOM格式文件。但需要注意，它不是被设计为查看图像，主要是用来操作DICOM文件的各种数据元素。
导入PyDicom包
import pydicom
打开dicom文件，并查看dicom文件内记录的信息
ds = pydicom.dcmread(file_path, force=True) # 读取dcm
patient_id = str(ds.PatientID)
patient_sex = str(ds.PatientSex)
patient_age = str(ds.PatientAge)
此时，可以采用打印的形式，将上面三个患者信息，打印到terminal上
print('patient info:', patient_id, patient_sex, patient_age)
在读取dicom文件时候，除了采用pydicom.dcmread()进行dicom文件的读取外，也可以采用pydicom.read_file()进行读取，两者之间没有区别。
2.2、dicom文件转成png文件
在2.1中展示了dicom如何读取，以及对dicom文件内存储的信息进行打印。那么，储存图像信息的在哪里呢？需要怎么转化为png图片呢？
首先，图像的像素信息，是都储存在pixel_array里面的，可以采用ds.pixel_array去除信息，如下：
ori_img = np.array(ds.pixel_array)
要显示DICOM格式的医学图像，必须将原始图像数据经过一系列的转换才能得到可直接在显示设备上显示的数据(称之为P Values)。DICOM医学图像显示需要经过Modality LUT、VOI LUT、Presentation LUT三个转换过程，最终输出的P Values才是可以直接显示的图像数据。
其中，影响显示效果的，还有窗位和窗宽，其中：
- 窗位(Window Center/Window Level)：代表可视范围(或是感兴趣区域)的 CT 值范围中心
- 窗宽(Window Width)：可视范围大小
对于这部分的详细内容，可以参考我的另一篇文件，【DICOM医学影像1】数据格式存储于显示，基本知识科普指南
接下来，就对转图部分进行详细的介绍，和代码展示。在DICOM文件的tag文件里面，有一个和显示转换相关的参数，关系到如何转图，能够看到和dicom文件读取器显示一样的图像。那就是VOI LUT Function定义了如何利用Window Center和Window Width进行映射，该字段有三种形式，
LINEAR；LINEAR_EXACT；SIGMOID。
如果没有给出这个字段，则按照WC，WW默认的方式处理(LINEAR)。处理方式如下：
import pydicom, cv2
import numpy as np
import matplotlib.pyplot as plt
class ConvertClass(object):
def __init__(self, ds):
super(ConvertClass, self).__init__()
self.pixel_array = ds.pixel_array
if type(ds.WindowCenter).__name__ == 'DSfloat':
self.wc = int(ds.WindowCenter)
self.ww = int(ds.WindowWidth)
elif type(ds.WindowCenter).__name__ == 'MultiValue':
self.wc = int(ds.WindowCenter[0])
self.ww = int(ds.WindowWidth[0])
self.ymax = 255
self.ymin = 0
if hasattr(ds, 'RescaleSlope'):
self.slope = ds.RescaleSlope
else:
self.slope = 1
if hasattr(ds, 'RescaleIntercept'):
self.intercept = ds.RescaleIntercept
else:
self.intercept = 0
if hasattr(ds, 'PhotometricInterpretation'):
self.PhotometricInterpretation = ds.PhotometricInterpretation
else:
self.PhotometricInterpretation = 'other'
if hasattr(ds, 'VOILUTFunction'):
self.VOILUTFunction = ds.VOILUTFunction
else:
self.VOILUTFunction = 'LINEAR'
# VOILUTFunction linear_exact
def linear_exact(self):
pixel_array = self.pixel_array * self.slope + self.intercept
linear_exact_array_less_idx = pixel_array <= (self.wc - self.ww / 2)
linear_exact_array_large_idx = pixel_array > (self.wc + self.ww / 2)
linear_exact_array = ((self.pixel_array - self.wc) / self.ww + 0.5) * (self.ymax - self.ymin) + self.ymin
linear_exact_array[linear_exact_array_less_idx] = self.ymin
linear_exact_array[linear_exact_array_large_idx] = self.ymax
linear_exact_array = linear_exact_array.astype(np.uint8)
if self.PhotometricInterpretation == 'MONOCHROME1':
linear_exact_array = 255 - linear_exact_array
return linear_exact_array
# VOILUTFunction sigmoid
def sigmoid(self):
pixel_array = self.pixel_array * self.slope + self.intercept
sigmoid_array = (self.ymax - self.ymin)/(1+np.exp(-4*(pixel_array-self.wc)/self.ww)) + self.ymin
sigmoid_array = sigmoid_array.astype(np.uint8)
if self.PhotometricInterpretation == 'MONOCHROME1':
sigmoid_array = 255 - sigmoid_array
return sigmoid_array
# VOILUTFunction linear
def linear(self):
pixel_array = self.pixel_array * self.slope + self.intercept
linear_array_less_idx = pixel_array <= (self.wc - self.ww/2)
linear_array_large_idx = pixel_array > (self.wc + self.ww/2 - 1)
linear_array = ((pixel_array - (self.wc - 0.5))/(self.ww - 1) + 0.5) * (self.ymax - self.ymin) + self.ymin
linear_array[linear_array_less_idx] = self.ymin
linear_array[linear_array_large_idx] = self.ymax
linear_array = linear_array.astype(np.uint8)
if self.PhotometricInterpretation == 'MONOCHROME1':
linear_array = 255 - linear_array
return linear_array
def noWW(self):
min_value = np.min(self.pixel_array)
max_value = np.max(self.pixel_array)
pixel_array = (self.pixel_array - min_value) / (max_value - min_value) * 255
pixel_array = pixel_array.astype(np.uint8)
if self.PhotometricInterpretation == 'MONOCHROME1':
pixel_array = 255 - pixel_array
return pixel_array
ds = pydicom.read_file(r"01.dcm")
convert = ConvertClass(ds)
noWW_arr = convert.noWW()
sigmoid_array = convert.sigmoid()
linear_exact_array = convert.linear_exact()
linear_array = convert.linear()
cv2.imwrite(r'01_no.png', noWW_arr)
cv2.imwrite(r'01_l_e.png', linear_exact_array)
cv2.imwrite(r'01_sig.png', sigmoid_array)
cv2.imwrite(r'01_l.png', linear_array)
同一个dicom文件，经过不同的处理方式，最终转换出来的图像的差异性，如下：
上述代码实现的计算公式部分，详情可以参考文章：dicom image display
将VOILUTFunction加入进来，汇总后的代码（可以直接使用下面这个代码），如下：
import pydicom, cv2
import numpy as np
import matplotlib.pyplot as plt
class ConvertClass(object):
def __init__(self, ds):
super(ConvertClass, self).__init__()
self.pixel_array = ds.pixel_array
if type(ds.WindowCenter).__name__ == 'DSfloat':
self.wc = int(ds.WindowCenter)
self.ww = int(ds.WindowWidth)
elif type(ds.WindowCenter).__name__ == 'MultiValue':
self.wc = int(ds.WindowCenter[0])
self.ww = int(ds.WindowWidth[0])
self.ymax = 255
self.ymin = 0
if hasattr(ds, 'RescaleSlope'):
self.slope = ds.RescaleSlope
else:
self.slope = 1
if hasattr(ds, 'RescaleIntercept'):
self.intercept = ds.RescaleIntercept
else:
self.intercept = 0
if hasattr(ds, 'PhotometricInterpretation'):
self.PhotometricInterpretation = ds.PhotometricInterpretation
else:
self.PhotometricInterpretation = 'other'
if hasattr(ds, 'VOILUTFunction'):
self.VOILUTFunction = ds.VOILUTFunction
else:
self.VOILUTFunction = 'LINEAR'
# VOILUTFunction linear_exact
def linear_exact(self):
pixel_array = self.pixel_array * self.slope + self.intercept
linear_exact_array_less_idx = pixel_array <= (self.wc - self.ww / 2)
linear_exact_array_large_idx = pixel_array > (self.wc + self.ww / 2)
linear_exact_array = ((self.pixel_array - self.wc) / self.ww + 0.5) * (self.ymax - self.ymin) + self.ymin
linear_exact_array[linear_exact_array_less_idx] = self.ymin
linear_exact_array[linear_exact_array_large_idx] = self.ymax
linear_exact_array = linear_exact_array.astype(np.uint8)
if self.PhotometricInterpretation == 'MONOCHROME1':
linear_exact_array = 255 - linear_exact_array
return linear_exact_array
# VOILUTFunction sigmoid
def sigmoid(self):
pixel_array = self.pixel_array * self.slope + self.intercept
sigmoid_array = (self.ymax - self.ymin)/(1+np.exp(-4*(pixel_array-self.wc)/self.ww)) + self.ymin
sigmoid_array = sigmoid_array.astype(np.uint8)
if self.PhotometricInterpretation == 'MONOCHROME1':
sigmoid_array = 255 - sigmoid_array
return sigmoid_array
# VOILUTFunction linear
def linear(self):
pixel_array = self.pixel_array * self.slope + self.intercept
linear_array_less_idx = pixel_array <= (self.wc - self.ww/2)
linear_array_large_idx = pixel_array > (self.wc + self.ww/2 - 1)
linear_array = ((pixel_array - (self.wc - 0.5))/(self.ww - 1) + 0.5) * (self.ymax - self.ymin) + self.ymin
linear_array[linear_array_less_idx] = self.ymin
linear_array[linear_array_large_idx] = self.ymax
linear_array = linear_array.astype(np.uint8)
if self.PhotometricInterpretation == 'MONOCHROME1':
linear_array = 255 - linear_array
return linear_array
def noWW(self):
min_value = np.min(self.pixel_array)
max_value = np.max(self.pixel_array)
pixel_array = (self.pixel_array - min_value) / (max_value - min_value) * 255
pixel_array = pixel_array.astype(np.uint8)
if self.PhotometricInterpretation == 'MONOCHROME1':
pixel_array = 255 - pixel_array
return pixel_array
def convert(self):
if self.VOILUTFunction=='SIGMOID':
sigmoid_array = self.sigmoid()
return sigmoid_array
elif self.VOILUTFunction=='LINEAR_EXACT':
linear_exact_array = self.linear_exact()
return linear_exact_array
else:
linear_array = self.linear()
return linear_array
ds = pydicom.read_file(r"01.dcm")
convert = ConvertClass(ds)
image_arr = convert.convert()
cv2.imwrite(r'image.png', image_arr)
根据不同dicom在数据采集阶段，设定的参数不同，采用针对性的转图方式，尤为重要。尤其是在为了保持转图不失真的情况下，必须这么操作。
2.3、转储图像做直方图均衡化
在上文2.2中，已经可以将我们看到，显示在microDicom软件的图像，转化为了png图像，存储到了本地。但是，由于拍摄环境等等的原因，会导致最后呈现的图像发生整体偏暗、或者整体偏亮，对最后的诊断产生不必要的干扰。
外部原因可能是曝光过高、或者曝光多低，会导致这个问题。一般医院遇着这种较为严重的情况，会采取重新拍摄的方式，而对于不那么严重的，则会继续诊断。毕竟再拍摄一次，会增加一次患者射线辐射。
此时，如果能够对这种质量不那么好的片子，采用后处理的方式，降低影响，是最好的了。于是，采用直方图均衡化，对这种像素分布过于集中的片子进行后处理，很有必要。
完整代码，如下所示
import pydicom
import skimage.transform as transform
import copy
import os
import numpy as np
import cv2
def over_lap(x1,y1,x2,y2):
if max(x1,y1) <= min(x2,y2) or min(x1,y1) >= max(x2,y2):
return 0
else:
delta = max(max(x1,y1) - min(x1,y1), max(x2,y2) - min(x2,y2))
op = min(abs(y2 - x1) + abs(y1 - x2), abs(y2 - y1) + abs(x2 - x1))
if (delta - op)/delta < 0.2:
return 0
else:
return 1
def dcm2png(file):
print('FIle:', file)
ds = pydicom.dcmread(file, force=True)
ds.file_meta.TransferSyntaxUID = pydicom.uid.ImplicitVRLittleEndian
ori_img = np.array(ds.pixel_array)
sharp = ori_img.shape
_h = sharp[0]
_w = sharp[1]
if len(sharp) == 3:
ori_img = ori_img[:, :, 0]
img = transform.resize(ori_img, (256, 256))
img1024 = transform.resize(ori_img, (_h, _w))
start = img.min()
end = img.max()
print("img shape:", img.shape)
vals = img.flatten()
counts, bins = np.histogram(vals, bins=256) # 一个参数是数据集，第二个参数是划分数据的bins个数或边缘值
temp_counts = copy.copy(counts) # 第一个频数的列表，第二个元素是数据的范围。
print(len(counts), len(bins))
hist_accum = np.zeros(256)
temp_sum = 0
for i in range(2, len(counts) - 2): # 计算累积直方图
temp_sum = temp_sum + counts[i]
hist_accum[i] = temp_sum / (65536.0 - counts[0] - counts[1] - counts[-2] - counts[-1])
# print(i, ':', hist_accum[i])
# 找到直方图中分布在0.3~0.7之间的像素值范围
try:
count_start = np.where(hist_accum > 0.3)[0][0]
count_end = np.where(hist_accum > 0.7)[0][0]
except Exception as e:
print("get histgram start and end error:", e)
# 将temp_counts中分布在count_start和count_end之外的像素值频数设为0
temp_counts[0:count_start] = 0
temp_counts[count_end:256] = 0
max_v = temp_counts.max() # 获取temp_counts中的最大值
loc = np.where(counts == max_v) # 获取max_v在counts中的位置
hist_sum_left = sum(counts[0: loc[0][0] + 1])
hist_sum_right = 65536 - hist_sum_left
# 根据左右两侧像素值的频数比例，确定start和end的值
for inx in range(loc[0][0], 1, -1):
left_ratio = hist_sum_left / 65536
if counts[inx] + counts[inx - 1] < 50 and left_ratio < 0.05:
start = (bins[inx] + bins[inx + 1]) / 2.0
break
hist_sum_left = hist_sum_left - counts[inx]
for inx in range(loc[0][0] + 1, 255):
right_ratio = hist_sum_right / 65536
if counts[inx] + counts[inx + 1] < 50 and right_ratio < 0.05:
end = (bins[inx] + bins[inx + 1]) / 2.0
break
hist_sum_right = hist_sum_right - counts[inx]
try:
print('WC type:', type(ds.WindowCenter).__name__)
if type(ds.WindowCenter).__name__ == 'DSfloat':
wc = int(ds.WindowCenter)
wl = int(ds.WindowWidth)
elif type(ds.WindowCenter).__name__ == 'MultiValue':
wc = int(ds.WindowCenter[0])
wl = int(ds.WindowWidth[0])
low = int(wc - wl / 2)
up = int(wc + wl / 2)
print("low,up:", low, up)
# 如果直方图分布范围太小，则使用窗位和窗宽参数代替start和end的值
if count_end - count_start < 4:
start = low
end = up
# 如果start和end与窗位和窗宽参数有重叠，则使用窗位和窗宽参数代替start和end的值
elif over_lap(low, up, start, end):
start = low
end = up
except Exception as e:
print(e, " get window error")
# 将img1024中小于start的像素值设为start，大于end的像素值设为end，并将像素值归一化到[0,255]范围内
img1024[img1024 < start] = start
img1024[img1024 > end] = end
img1024 = np.array((img1024 - start) * 255.0 / (end - start))
if hasattr(ds, 'PhotometricInterpretation'):
if ds.PhotometricInterpretation == 'MONOCHROME1':
img1024 = 255 - img1024
img1024 = img1024.astype(np.uint8)
cv2.imwrite(r'F:\huangao.jpeg', img1024)
dcm2png('huangao.dcm')
展示均衡化调整后的对比，如下：
可以看出，直方图均衡化在该片中，处理后的效果并不好，甚至有些曝光过度了，尤其是腕骨部分。
所以说，并不是所有的情况都是适合采用该方法的，需要针对性的进行调整，一般状况下可不调整。
三、生成dicom文件
既然dicom文件可以另存为普通的png文件，那么，png文件是不是也可以存储为dicom文件呢？通过前面对raw和mhd文件分别存储信息的学习，知道肯定是可以的。
在读取dicom文件阶段，先取出整体，然后采用后缀的形式，不断的取出需要的属性。再生成dicom阶段，恰好是解析dicom的反过程。他需要先将需要的各个属性，加入到Dataset，然后再储存，即可生成想要的dicom了。
对于dicom内部的这些参数都是什么意思感兴趣的，可以参考这篇文章：DICOM内部信息详解篇
3.1、jpeg图像生成dicom文件
通过读取一种jpeg图片，把它转成dicom文件，然后用microDicom查看，完整的代码如下所示：
import os
import datetime
import numpy as np
from PIL import Image
import pydicom
from pydicom.dataset import Dataset, FileDataset
filename_little_endian ='bb.dcm'
filename_big_endian = 'testr.dcm'
def save_to_dcm(image_path):
img = np.asarray(Image.open(image_path).convert("L"), dtype=np.uint8)
print(img.shape)
print("Setting file meta information...")
# Populate required values for file meta information
file_meta = Dataset()
file_meta.MediaStorageSOPClassUID = '1.2.840.10008.5.1.4.1.1.2'
file_meta.MediaStorageSOPInstanceUID = "1.2.3"
file_meta.ImplementationClassUID = "1.2.3.4"
print("Setting dataset values...")
# Create the FileDataset instance (initially no data elements, but file_meta
# supplied)
ds = FileDataset(filename_little_endian, {},
file_meta=file_meta, preamble=b"\0" * 128)
# # Write as a different transfer syntax XXX shouldn't need this but pydicom
# # 0.9.5 bug not recognizing transfer syntax
ds.file_meta.TransferSyntaxUID = pydicom.uid.ImplicitVRLittleEndian
# Set creation date/time
dt = datetime.datetime.now()
ds.ContentDate = dt.strftime('%Y%m%d')
timeStr = dt.strftime('%H%M%S.%f') # long format with micro seconds
ds.ContentTime = timeStr
ds.Modality = "US"
ds.PatientName = "ZhangSan"
ds.patientAge = "105"
ds.SeriesInstanceUID = pydicom.uid.generate_uid()
ds.StudyInstanceUID = pydicom.uid.generate_uid()
ds.FrameOfReferenceUID = pydicom.uid.generate_uid()
ds.BitsStored = 8
ds.BitsAllocated = 8
ds.SamplesPerPixel = 3
ds.HighBit = 7
ds.PlanarConfiguration = 0
ds.InstanceNumber = 1
ds.ImagePositionPatient = r"0\0\1"
ds.ImageOrientationPatient = r"1\0\0\0\-1\0"
ds.ImagesInAcquisition = "1"
ds.Rows, ds.Columns = img.shape[:2]
ds.PixelSpacing = [1, 1]
ds.PixelRepresentation = 0
ds.PixelData = img.tobytes()
# import pdb;pdb.set_trace()
# ds.PhotometricInterpretation = 'RGB'
print("Writing file", 'test.dcm')
ds.save_as('test.dcm')
print("File saved.")
if __name__ == '__main__':
image_path = r"00d35ceb3c14c3318f4f3ec998879827.jpeg"
save_to_dcm(image_path)
3.2、操作dicom文件图像重写入
dicom文件是一个规范性的存储格式。那能不能通过操作，将一些附加的信息，添加到dicom文件内呢？简言之就是实现了对dicom文件的修改。
下面展示了如何对dicom文件的图像部分进行修改。当然，也可以对附加信息部分进行修改，比如改写studyDate。后半部分比较简单，这里不再赘述了。着重介绍前半部分修改影像吧。下面的演示就是在图像中添加两条直线，然后再把dicom重新保存。代码如下：
# 获取dcm文件的pixel_array
filename_dcm = "a.dcm"
ds = pydicom.dcmread(os.path.join(dcm_Dir, filename_dcm)) # 读取dcm
img_dcm = np.array(ds.pixel_array)
# 采用OpenCV方式，对pixel_array进行划直线操作
cv2.line(img_dcm, (x1, line1_y), (raw_w-x1, line1_y), (0, 0, 255))
cv2.line(img_dcm, (x1, line2_y), (raw_w-x1, line2_y), (0, 0, 255))
# 将修改后的pixel_array，重新保存到dcm文件内，其他内容不发生改变
ds.PixelData = img_dcm.tobytes()
ds.Rows, ds.Columns = img_dcm.shape
ds.save_as('a_r.dcm')
print("File saved.")
这种方式需要对dicom操作的内容比较少，原dicom的属性都还是在的，只是对dicom内的pixel_array进行了操作，然后再重新还给dicom文件，保存好即可。
microDicom软件打开，前后的展示如下：
四、总结
本文基本上围绕着dicom文件格式，前前后后都介绍到了。尤其是一些关于图像操作的部分，想必对你的操作会有较大的帮助的。对于文中的一些参考链接，也是一个非常好的学习工具，会介绍的更加的详细，建议可以着重看看。
