[DICOM活久见] 序列内部的RescaleIntercept不同导致的问题

本文由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

关于PixelRepresentation这个tag的说明，可以查看网站上的说明:

它主要是指示图像里面像素的值的类型．
如果PixelRepresentation是1, 则像素值可以为负数 (signed);
如果PixelRepresentation是0, 则像素值只能是正数 (unsigned).

而对于CT图像来说，空气的CT值(HU)是负数. 因此，如果是包含了很多空气的肺部，它的HU值也是负数．所以如果一个影像，它的PixelRepresentation是0的话，要想最后的HU值是负数，就需要通过RescaleIntercept是-1024, 才能计算出负数的HU.

这个PixelRepresentation, 一个序列内的影像应该都是相同的．
而RescaleSlope和RescaleIntercept, 也基本是相同的．

因此，之前在创建sitk Image时，便根据第一张图像的RescaleSlope和RescaleIntercept, 来作为整个序列的补偿值．

但正是由于这个假设，导致了上面的问题．

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

通过dcmdump工具，可以查看这个序列下面所有图像的某一个或几个tag值，结果当我查看这个序列的RescaleIntecept时，发现一个包含了134张的序列，它里面的有很多个值．

因此，当我按照一个序列内，所有影像都是按照第一张图像的RescaleIntercept来计算的话，就会导致其他层的图像，在计算HU值时，出现较大的偏差．

所以，最终的修复方案，就是在构建sitk Image时，在构建image array时，每一张图像，都单独计算，单独读取每一张图像的dicom header, 获取它的RescaleSlope和RescaleIntercept, 然后计算单张的image的HU值，最后再构建出整个sitk Image即可．