dicom image display

dicom image display

主要关注的是pixel data到gray data的过程。对于dicom tag的解析此处并不关注。

DICOM image display - DCMTK - three conversions of pixel data (programming.vip)

[Medical] DICOM - 醫學影像的窗位、窗寬、調整斜率與調整截距心得筆記 | 記憶裂縫 - 點部落 (dotblogs.azurewebsites.net)

C.11 Look Up Tables and Presentation States (nema.org)

上面的连接中介绍了dicom影像从裸数据到最后的灰度图的变换过程。整个流程以及相关的dicom tag如下图所示，dicomtag 的具体说明可以见：
VOI LUT Sequence Attribute – DICOM Standard Browser (innolitics.com)

关于上述tag的使用规则如下：

1. Window Center (0028,1050) 与 VOI LUT Sequence (0028,3010) 仅留一个使用
2. Window Center (0028,1050) 需要和 Window Width (0028,1051) 搭配使用
3. Rescale Intercept (0028,1052) 与 Modality LUT Sequence (0028,3000) 仅留一个使用
4. Rescale Intercept (0028,1052) 需要和 Rescale Slope (0028,1053) 搭配使用

Modality LUT 映射

rescale intercept和rescale slope首先对原始数据进行线性变换:

value1 = pixel_data * rescale_slope + rescale_intercept

在不存在rescale intercept字段的情况下，可能会有Modality LUT Sequence。

(0028,3002) LUT Description: 具有三个值，第一个表示lookup table的元素个数。当table的元素个数为2^16，那么这个是为0；第二个表示映射到lut第一个值的pixel value。pixel data中所有小于这个值的都会被映射到table中的第一个元素大小。像素值number of entries + first value mapped - 1会被映射到lut最后一个值。大于这个像素值的也会被映射到lut的最后一个值。最后一个值，表示每个元素占用的bits数量。通常是8或16.

具体的映射关系如下：

num_elements = lut_descriptor[0]
low_pixel_value = lut_descriptor[1]
num_bit = lut_descriptor[2]
lut_data

if (pixel_value <= low_pixel_value)
	value1 = lut_data[0]
else if (pixel_value - low_pixel_value > num_elements - 1)
	value1 = lut_data[num_elements-1]
else
	value1 = lut_data[pixel_value - low_pixel_value]

VOI LUT 映射

VOI LUT Sequence和先前介绍的一致，此处不做额外介绍。

VOI LUT FUNCTION = LINEAR (or default)

当使用Window Center 和 Window Width时。数值的映射是这样的。输出值的范围是ymin, ymax，那么：

if (value1 <= wc - 0.5 - (ww - 1)/2) // wc - 0.5 - ww/2 + 0.5 = wc - ww/2
	value2 = ymin
else if (value1 > wc-0.5 + (ww-1)/2)  // wc - 0.5 + ww/2 - 0.5 = wc + ww/2 - 1
	value2 = ymax
else
	value2 = ((value1 - (wc-0.5))/(ww-1) + 0.5) * (ymax - ymin) + ymin

VOI LUT Function定义了如何利用Window Center和Window Width进行映射，该字段有三种形式，LINEAR，LINEAR_EXACT，SIGMOID。如果没有给出这个字段，则按照WC，WW默认的方式处理(LINEAR)。

VOI LUT Function Attribute – DICOM Standard Browser (innolitics.com)

VOI LUT FUNCTION = SIGMOID

如果VOI LUT FUNCTION等于SIGMOID，那么转换公式如下：

y = (ymax - ymin)/(1+exp(-4*(x-c)/w)) + ymin

• x：是即将执行LUT变换的输入，如Modality LUT变换后的结果，即上述的value1
• c：window center，(0028,1050).
• w：window width，(0028,1051).
• y：转换后的输出；
• ymin：输出的最小值；
• ymax：输出的最大值；

VOI LUT FUNCTION = LINEAR_EXACT

此时转换公式如下：

if (x <= c - w/2)
	y = ymin;
else if (x > c + w/2)
	y = ymax;
else
	y = ((x - c)/w + 0.5) * (ymax - ymin) + ymin

Presentation LUT映射

其中LUT Sequence中存储的内容和上述的LUT一致。如果Presentation LUT Shape为IDENTITY，无需额外处理，如果为INVERSE，需要对结果进行inversion操作。

Python implementation

本节主要利用python实现，对比一下 VOI LUT FUNCTION = LINEAR和LINEAR_EXACT情况下CT影像的差距。

import pydicom
import numpy as np

import matplotlib.pyplot as plt

dcmfile = pydicom.read_file(r"test.dcm")

ww = dcmfile.WindowWidth
wc = dcmfile.WindowCenter

slope = dcmfile.RescaleSlope
intercept = dcmfile.RescaleIntercept

ymin = 0
ymax = 255

pixel_array = dcmfile.pixel_array * slope + intercept

# linear exact array
linear_exact_array = np.zeros(pixel_array.shape)

linear_exact_array_less_idx = pixel_array <= (wc - ww/2)
linear_exact_array_large_idx = pixel_array > (wc + ww/2)
linear_exact_array = ((pixel_array - wc)/ww + 0.5) * (ymax - ymin) + ymin
linear_exact_array[linear_exact_array_less_idx] = ymin
linear_exact_array[linear_exact_array_large_idx] = ymax
linear_exact_array = linear_exact_array.astype(np.uint8)

# linear array
linear_array = np.zeros(pixel_array.shape)
linear_array_less_idx = pixel_array <= (wc - ww/2)
linear_array_large_idx = pixel_array > (wc + ww/2 - 1)
linear_array = ((pixel_array - (wc - 0.5))/(ww - 1) + 0.5) * (ymax - ymin) + ymin
linear_array[linear_array_less_idx] = ymin
linear_array[linear_array_large_idx] = ymax
linear_array = linear_array.astype(np.uint8)

diff_array = linear_array - linear_exact_array

plt.figure()
plt.imshow(linear_exact_array, cmap=plt.cm.gray, interpolation='bilinear')
plt.title("linear_exact_array")
plt.show()

plt.figure()
plt.imshow(linear_array, cmap=plt.cm.gray, interpolation='bilinear')
plt.title("linear_array")
plt.show()

plt.figure()
plt.imshow(diff_array, cmap=plt.cm.gray, interpolation='bilinear')
plt.title("diff_array")
plt.show()

实现结果如下：

针对这种情况，肉眼并不能明显区分两者之间的差别。