本文由Markdown语法编辑器编辑完成.
1. 背景:
本文记录在工作中遇到的一些比较罕见的dicom图像.
这对于在未来工作中, 处理图像时, 需要考虑方案的完整性, 会有很大的帮助.
本文介绍的, 是目前我工作10年来, 头一次见到的一个CT序列, 它的序列内的RescaleIntercept值, 不是完全相同的. 这会导致如果在生成体数据时, 如果按照第一层的RescaleIntercept的值, 来作为整个体数据的RescaleIntercept时, 会引起很大偏差的问题.
2. 问题描述:
在医学图像处理中, 我们通用会将一个序列内的所有dicom图像, 生成一个包含该序列内所有影像像素信息的体数据. 这样无论对于后期的读取, 计算, 都可以基于这个体数据来进行.
生成体数据, 主要是基于simpleITK提供的方法来进行, 大致逻辑如下:
def get_sitk_image_from_series(series_dataset_list: list):
"""
根据series内每一张图像的信息, 得到包含一个序列内所有图像信息的np array.
Args:
series_dataset_list (str): 根据ImagePositionPatient的z值递增排序后的数组.
"""
bottom_dicom = series_dataset_list[0]
try:
image_count = len(series_dataset_list)
rows, cols = bottom_dicom.get("Rows"), bottom_dicom.get("Columns")
rescale_intercept, rescale_slope = None, None
if not bottom_dicom.get("RescaleIntercept") or not bottom_dicom.get("RescaleSlope"):
rescale_intercept = 0.0
rescale_slope = 1.0
else:
rescale_intercept, rescale_slope = bottom_dicom.get("RescaleIntercept"), bottom_dicom.get("RescaleSlope")
series_array = np.zeros((image_count, rows, cols), dtype=np.int16)
slice_location_list = list()
for index, img_dataset in enumerate(series_dataset_list):
img_pixel_array = img_dataset.pixel_array * rescale_slope + rescale_intercept
slice_location_list.append(img_dataset.get("ImagePositionPatient")[-1])
series_array[index] = img_pixel_array
# origin direcation spacing
top_dicom = series_dataset_list[-1]
direction = cal_direction_with_orientation(bottom_dicom.get("ImageOrientationPatient", [1, 0, 0, 0, 1, 0]))
sitk_img = sitk.GetImageFromArray(series_array)
sitk_img.SetDirection(direction)
sitk_img.SetOrigin(bottom_dicom.get("ImagePositionPatient"))
z_spacing = (top_dicom.get("ImagePositionPatient")[-1] - bottom_dicom.get("ImagePositionPatient")[-1]) / (
len(series_dataset_list) - 1
)
pixel_spacing = list(copy.deepcopy(img_dataset.get("PixelSpacing")))
pixel_spacing.append(z_spacing)
sitk_img.SetSpacing(pixel_spacing)
except Exception:
return None
return sitk_img
def cal_direction_with_orientation(orientation):
row_direction = np.array(orientation[:3])
column_direction = np.array(orientation[3:])
slice_direction = np.cross(row_direction, column_direction)
row_direction /= np.linalg.norm(row_direction)
column_direction /= np.linalg.norm(column_direction)
slice_direction /= np.linalg.norm(slice_direction)
return list(row_direction) + list(column_direction) + list(slice_direction)
这一段逻辑, 主要是将将一个序列内读入的, image_dataset的list, 构建出一个sitkImg, 然后再将sitkImg, 转存成为一个.mha, .niig.gz等体数据的文件.
这样就实现了, 将一系列二维的图像, 生成一个三维的体数据文件, 便于后续的读取和操作.
那么如何再将.mha的体数据, 还原为原始的dicom文件呢?
主要就是从.mha中, 读取出特定层数的像素信息, 然后再加上之前的dicom header信息, 就可以组合为一张一张dicom图像.
因为有时候, 也是需要读取单张dicom图像的.
那么现在问题来了. 在处理一套病例数据时, 发现经过这样的处理后, 返回的dicom文件, 出现了白图的情况. 但是如果只是查看.mha文件, 又是正常的.
那么究竟是哪里出了问题呢?
2.1 问题表现
通过从浏览器，可以查看到，通过.mha和dicom header,组合后生成的dicom图像．
第一幅图像，是经过还原后的异常图像．图像出现了很严重的反白现象；而第二幅图像，则是原始的dicom图像，非常正常．
那么我们再看一下，通过原始的dicom图像，生成的.mha是否正常．将生成的.mha文件，拖入slicer中，可以看到:
从slicer中看，体数据也是"没问题"的．
至少从大体上看，它是能够看到各个部位的轮廓的．
为了对比，我再把原始序列的dicoms文件，也加载到slicer里面看一下．
但是仔细观察，slicer加载: .mha和原始dicom后，图像的亮度是有差异的．
那么，观察左侧的图像的灰度范围，可以看出两者的差异:
| 图像类型 | 灰度下限 | 灰度上限 |
|---|---|---|
| .mha | -3896 | 9772 |
| 原始dicom | -3986 | 9811 |
可以看出，两者之间的差异，主要体现在了灰度的最大值不同，差了大概40多．虽然不是很多，但是也是有差异的．
按理来说，这里的.mha文件，是原始的dicom文件，在仅保留像素的情况下，生成的．所以，生成的.mha文件，和原始的dicoms文件，加载到slicer里面，应该是完全一致，灰度的范围也是要完全一致，才可以的．
那么，现在的问题，就是需要定位，这相差的40的灰度值，是如何引起的．
2.2 问题定位:
通过查看生成.mha体数据的逻辑．
生成.mha文件前，需要先基于读取的dicoms的灰度信息，构建出一个sitk的Image, 然后将这个Image, 写成一个.mha文件．
那么问题的关键，就是构建这个sitk Image的逻辑了．
图像原始的灰度信息，是存在了:
pixel_array中．而构建sitk Image时，需要将pixel_array中的灰度，经过RescaleSlope和RescaleIntercept的调节，生成HU值．
计算的方式为:
HU = pixel_array * RescaleSlope + RescaleIntercept
对于CT图像来说，一般是两类型的值．而这个值的取值，又依赖于PixelRepresentation的值.
| PixelRepresentation | RescaleSlope | RescaleIntercept |
|---|---|---|
| 1 | 1.0 | 0 |
| 0 | 1.0 | -1024 |
