医学图像处理 pydicom

DICOM

一个 DICOM 文件类似一个类，会有很多属性，属性可能是一个字段，数组或者包含了另外一个类。编程语言会定义基础数据类型，类似的 DICOM 定义很多基础类型叫做 VR - Value Representations 值表现。

DCM 文件中具体数据的处理，说图像可能有些狭隘，广义上还包括波形，如心电图、超声成像的视频等等。

图像重建技术，平面、曲面、表面、三维，已经成为时下放射领域最热门的词。它们也是数字成像对比从前的硬拷贝胶片最独特的优势之一，旧时你根本没法后处理。如果 DICOM 数据对象中没有收集各种信息，那么重建后处理是根本不可能实现的。

考虑到一个简单的数字成像，比如 CT 扫描。除它的像素属性之外，DICOM 会记录所有 有关的长度、三维坐标和方向。

利用 pydicom 库解析 DICOM 文件，自带 dcm 测试文件，各文件信息参考 README 文件。

安装 pydicom 和 numpy 库：

pip install pydicom numpy

单 dcm 文件读取：

import pydicom
import pylab

dcm = 'C:\\Anaconda3\\Lib\\site-packages\\pydicom\\data\\test_files\\RG1_UNCI.dcm'
ds = pydicom.read_file(dcm)
print(ds.dir()) 
print(ds.dir('pat')) 
print(ds.PatientName, ds.PatientSex, ds.StudyDate, ds.PatientBirthDate)

print(ds.data_element('PatientID'))
print(ds.data_element('PatientID').VR, ds.data_element('PatientID').value)

pixel_bytes = ds.PixelData 
pix = ds.pixel_array       
print(pix.shape) 
pylab.imshow(pix, cmap=pylab.cm.bone)
pylab.show()

系列 dicom 文件处理：

import os
import pydicom
import numpy
from matplotlib import pyplot

dcm = 'C:\\Anaconda3\\Lib\\site-packages\\pydicom\\data\\test_files'
lstFilesDCM = []

# 将所有dicom文件读入
for diName, subdirList, fileList in os.walk(dcm):
    for filename in fileList:
        if ".dcm" in filename.lower():
            print(filename)
            lstFilesDCM.append(os.path.join(diName, filename))

## 将第一张图片作为参考图
RefDs = pydicom.read_file(lstFilesDCM[10])
# print(RefDs)
# print(RefDs.pixel_array)
# print(RefDs.PatientPosition)
pyplot.imshow(RefDs.pixel_array, cmap=pyplot.cm.bone)
pyplot.show()

# 建立三维数组，长、宽、层数
ConstPixelDims = (int(RefDs.Rows), int(RefDs.Columns), len(lstFilesDCM))
print(ConstPixelDims)

# 得到 spacing 值 mm 单位
# PixelSpacing - 每个像素点实际的长度与宽度,单位(mm)
# SliceThickness - 每层切片的厚度,单位(mm)
ConstPixelSpacing = (float(RefDs.PixelSpacing[0]), float(RefDs.PixelSpacing[1]), float(RefDs.SliceThickness))

# 三维数据
# 0 到（第一个维数加一*像素间的间隔），步长为 ConstPixelSpacing
x = numpy.arange(0.0, (ConstPixelDims[0] + 1) * ConstPixelSpacing[0], ConstPixelSpacing[0])
y = numpy.arange(0.0, (ConstPixelDims[1] + 1) * ConstPixelSpacing[1], ConstPixelSpacing[1])
z = numpy.arange(0.0, (ConstPixelDims[2] + 1) * ConstPixelSpacing[2], ConstPixelSpacing[2])
print(len(x),"xxxx")

ArrayDicom = numpy.zeros(ConstPixelDims, dtype=RefDs.pixel_array.dtype)

# 遍历所有的 dicom 文件，读取图像数据，存放在 numpy 数组中
for filenameDCM in lstFilesDCM:
    ds = pydicom.read_file(filenameDCM)
    ArrayDicom[:, :, lstFilesDCM.index(filenameDCM)] = ds.pixel_array


# 轴状面显示
# 每英寸的像素数 dpi - dot per inch 越大图片越清晰。
pyplot.figure(dpi=1000)
# 将坐标轴都变为同等长度
# pyplot.axes().set_aspect('equal', 'datalim')
pyplot.axes().set_aspect('equal')
pyplot.set_cmap(pyplot.gray()) # 灰度化图片

# 第三个维度表示现在展示的是第几层 
pyplot.imshow(ArrayDicom[:, :, 360])
pyplot.show()

# 冠状面显示
pyplot.figure(dpi=100)
pyplot.axes().set_aspect('equal', 'datalim')
pyplot.set_cmap(pyplot.gray())
pyplot.imshow(ArrayDicom[:, 90, :])
pyplot.show()

Pixel Spacing 像素间距属性 Dp 保存在 DICOM(0028,0030)。它定义了图像像素的物理大小并且保证了实际距离测量的准确性。比如，如果你知道 x 和 y 轴的像素间距为 0.4mm，那么在图像中的一条 10 像素的线就会有 4mm 的长度。同 样，由于你知道图像像素中的宽和高，比如对于普通 CT 来说是 512×512，你就能够找到图像的实际尺寸了:512 × 0.4 mm = 204.8 mm。

Image Position 图像位置 DICOM(0020,0032) 属性，Ip 。这表示图像最左上角第一个像素的 x、y、z 坐标， 单位为毫米。这可以让我们了解在三维空间内图像开始的位置。

图像方向 DICOM(0020,0037) 属性。它存储了图像行和列向量 vr 和 vc 在三维方向上的余弦。 这两个向量，源自属性图像位置点 Ip。其在三维空间完整地定义了整 个图像平面。现在，如果我们有一个图像像素 P 在第 r 行、第 c 列，那么我们就能找到它 的三维坐标了，方法如下:

P3D = Ip + r×vr + c×vc

在 x、y、z 三个坐标上分别计算。

层间距 DICOM(0018,0088) 属性，Ds，记录了连续各层图像之间的距离，单位为毫米。它提供了和像素间距 Pixel Spacing (0028,0030) 相同的作用，但是只在 z 方向上。比如， 如果你考虑在第一和第二个 CT 层面上那个相同的像素 (r,c) 位置，那么两层面之间的距离 将等于层间距 (0018,0088) 属性的值。

Tag(7FE0,0010) 即表示 Pixel Data tag 这个元素存储了图像数据

图像相关标记 tag

(0028,0002) Samples per Pixel 像素的取样数，一般 CT，MR，DR 等灰度图像都是 1，而彩超等彩图像都是 3，分别表示R, G, B三个颜色通道。

(0028,0004) Photometric Interpretation：
     Monochrome2 一般的灰度图像都采用这种，Pixel值越大，图像就越白。
     Monochrome1 只有部分CR, DR图像使用，Pixel值越大，图像就越黑。
     Palette Colour 一般用于彩超图像，每个像素占用8位或者16位

(0028,0006) Planar configuration 定义了各个彩色通道值在 Pixel Data 中排列的排列方式。

值为 0 排列为 RGBRGB...。 值为 1 的排列 RRRRR…GGGGG…BBBBB。 
对于多帧图像，它是这样排列的：第一帧 RRR…GGG…BBB…，第二帧 RRR…GGG…BBB…

(0028,0008) Number of Frames 获取帧数，然后，逐帧读取 Pixel Data。 

(0028,0010) Rows 图像的高度
(0028,0011) Columns 图像的宽度

(0028,0030) Pixel Spacing  图像像素间距，主要用于长度测量。
(0028,0100) Bits Allocated 一个像素取样点存储时分配到的 bit 位数。

(0028,0101) Bits Stored 一个像素取样点存储时使用到的 bit 位数。
(0028,0102) High Bit 最高位序号。
(0028,0103) Pixel Representation 值为 0 这表明是无符号类型，1 表明为有符号类型。
(0028,1050) Window Center 窗位
(0028,1051) Window Width 窗宽

(0028,1052) Rescale Intercept 和 
(0028,1053) Rescale Slope 用于根据像素值计算原始值，如 CT 可以用于计算 HU 值。

    HU = Rescale Slope * X + Rescale Intercept.

调色板保存在
(0028,1201) RedPaletteColorLookupTableData, 
(0028,1202) GreenPaletteColorLookupTableData, 
(0028,1203) BluePaletteColorLookupTableData

(0028,2110) Lossy Image Compression 值为 1 表明该图像曾经经过有损压缩处理。
(0028,2112) Lossy Image Compression Ratio 有损压缩压缩率。

色彩模型：

• RGB 这是最常用的彩色图像格式。

• YBR_FULL 另外一种彩色图像格式，存储格式为 Y(Luminance 亮度) B(Blueness 蓝色)， R（Redness, 红色）

• YBR_FULL_422 一般用于 JPG 有损压缩格式的彩色图像，每两个像素共同使用32位，每一个像素都有自己的亮度 Y Luminance，但是共享相同的蓝色 B Blueness，红色 R Redness 两个分量。所以，它的像素值存储方式是：YYBR,YYBT,YYBR

• YBR_RCT 用于JPEG 2000无损压缩彩色图像，Reversible Color Transformation, 可逆色彩变换。

    Y = (R+2G+B)/4，
    CB = B-G ，
    CR = R - G
    G = Y - (CR+CB)/4 ,
    R = CR + G, B = CB + G
  

• YBR_ICT 用于 JPEG 2000有损压缩彩色图像 Irreversible Color Transformation, 不可逆色彩变换。

    Y = + .29900R + .58700G + .11400B
    CB = - .16875R - .33126G + .50000B
    CR = + .50000R - .41869G - .08131B
  

影像的 Transfer Syntax UID (0002,0010) 決定了影像的存储方式。 Pixel data 的存储方式，包括：

• 未压缩：

  • Implicit VR Little Endian:
  • Default Transfer Syntax for DICOM 1.2.840.10008.1.2
  • Explicit VR Little Endian 1.2.840.10008.1.2.1
  • Explicit VR Big Endian 1.2.840.10008.1.2.2

• 无损压缩：

  • JPEG Lossless :
  • Default Transfer Syntax for Lossless JPEG
  • ImageCompression 1.2.840.10008.1.2.4.70
  • JPEG 2000 Image Compression
  • (Lossless Only) 1.2.840.10008.1.2.4.90
  • RLE Lossless 1.2.840.10008.1.2.5

• 有损压缩：

  • Default Transfer Syntax for Lossy
  • JPEG 12 Bit Image Compression

结合开源项目 IMBRA 讲解如何解析一个 DCM 文件。

文件开头会有 128 字节的导言，这部分数据没有内容。接着是 4 字节 DICOM 文件标识，存储这 DICM。然后紧接着就是 dicom 数据数据元素。

旧版本的DCM文件会在开头 8 个字节验证签名，在读完 128 字节后，读取 4 个字节验证 DICOM 标识。

接下来两个字节的 tagId，并且验证大小端。

dicom 数据元素的 tagId 是从大到小的读取的，而最小的 tagId 就是从 0x0002，当 tagId 不是 0x0002 说明已经读取完所以关于 0x0002 的 tagId。从中找出 (0002,0010) 的 tag，这个 tag 设置当前数据是否大端格式，显示 VR 编码还是隐式 VR 编码。 "1.2.840.10008.1.2.2" 表示大端格式，1.2.840.10008.1.2 表示隐式编码。