从 3D DICOM 图像重建 2D X 射线

Question

我需要用下面的 Input/Output

编写一个 python 函数或 class

输入：

X射线源的位置（仍然不确定为什么需要它）
板的位置（仍然不确定为什么需要它）
三维 CT 扫描

输出：

A 2D X-ray Scan（模拟 X-ray 扫描，这是一种扫描全身的扫描）

关于我要实现的目标的一些重要说明：

您不需要来自现实世界的额外信息或任何高级知识。
您可以添加您认为合适的任何输入参数。
如果您的方法产生了瑕疵，您可以修复它们。
请解释你的方法的每一步。

到目前为止我做了什么：（添加了.py文件）

我已阅读位于 "Case2" 文件夹中的 .dicom 文件。

这些 .dicom 文件可以从我的 Google 驱动器下载： https://drive.google.com/file/d/1lHoMJgj_8Dt62JaR2mMlK9FDnfkesH5F/view?usp=sharing

我已经按位置对文件进行了排序。

最后，我创建了一个 3D 阵列，并将所有图像添加到该阵列以绘制结果（您可以在添加的图像中看到它们）- 这是 CT 扫描的切片。（参考：https://pydicom.github.io/pydicom/stable/auto_examples/image_processing/reslice.html#sphx-glr-auto-examples-image-processing-reslice-py）

完整代码如下：

import pydicom as dicom
import os
import matplotlib.pyplot as plt
import sys
import glob
import numpy as np
path = "./Case2"
ct_images = os.listdir(path)
slices = [dicom.read_file(path + '/' + s, force=True) for s in ct_images]
slices[0].ImagePositionPatient[2]
slices = sorted(slices, key = lambda x: x.ImagePositionPatient[2])

#print(slices)
# Read a dicom file with a ctx manager
with dicom.dcmread(path + '/' + ct_images[0]) as ds:
    # plt.imshow(ds.pixel_array, cmap=plt.cm.bone)
    print(ds)
    #plt.show()


fig = plt.figure()
for num, each_slice in enumerate(slices[:12]):
    y= fig.add_subplot(3,4,num+1)
    #print(each_slice)
    y.imshow(each_slice.pixel_array)
plt.show()    

for i in range(len(ct_images)):
    with dicom.dcmread(path + '/' + ct_images[i], force=True) as ds:
        plt.imshow(ds.pixel_array, cmap=plt.cm.bone)
        plt.show()       

# pixel aspects, assuming all slices are the same
ps = slices[0].PixelSpacing
ss = slices[0].SliceThickness
ax_aspect = ps[1]/ps[0]
sag_aspect = ps[1]/ss
cor_aspect = ss/ps[0]

# create 3D array
img_shape = list(slices[0].pixel_array.shape)
img_shape.append(len(slices))
img3d = np.zeros(img_shape)

# fill 3D array with the images from the files
for i, s in enumerate(slices):
    img2d = s.pixel_array
    img3d[:, :, i] = img2d

# plot 3 orthogonal slices
a1 = plt.subplot(2, 2, 1)
plt.imshow(img3d[:, :, img_shape[2]//2])
a1.set_aspect(ax_aspect)

a2 = plt.subplot(2, 2, 2)
plt.imshow(img3d[:, img_shape[1]//2, :])
a2.set_aspect(sag_aspect)

a3 = plt.subplot(2, 2, 3)
plt.imshow(img3d[img_shape[0]//2, :, :].T)
a3.set_aspect(cor_aspect)

plt.show()

结果不是我想要的，因为：

这些是 CT 扫描的切片。我需要模拟 X 射线扫描，这是一种遍历全身的扫描。

希望您能帮助模拟穿过人体的 X 射线扫描。

我读到它可以通过以下方式完成："A normal 2D X-ray image is a sum projection through the volume. Send parallel rays through the volume and add up the densities."我不确定它是如何在代码中完成的。

可能有帮助的参考资料：https://pydicom.github.io/pydicom/stable/index.html

Answer 1

编辑： 如进一步回答所述，此解决方案产生平行投影，而不是透视投影。

根据我对 "A normal 2D X-ray image" 定义的理解，这可以通过对给定方向上投影的每个切片的每个像素的每个密度求和来完成。

对于 3D 体积，这意味着在给定的轴上执行求和，这可以在 numpy 中使用 ndarray.sum(axis) 来完成。

# plot 3 orthogonal slices
a1 = plt.subplot(2, 2, 1)
plt.imshow(img3d.sum(2), cmap=plt.cm.bone)
a1.set_aspect(ax_aspect)

a2 = plt.subplot(2, 2, 2)
plt.imshow(img3d.sum(1), cmap=plt.cm.bone)
a2.set_aspect(sag_aspect)

a3 = plt.subplot(2, 2, 3)
plt.imshow(img3d.sum(0).T, cmap=plt.cm.bone)
a3.set_aspect(cor_aspect)

plt.show()

这会产生以下结果：

对我来说，这看起来像 X 射线图像。

编辑： 结果有点太 "bright"，因此您可能需要应用伽马校正。使用 matplotlib，import matplotlib.colors as colors 并在 plt.imshow 中添加一个 colors.PowerNorm(gamma_value) 作为 norm 参数：

plt.imshow(img3d.sum(0).T, norm=colors.PowerNorm(gamma=3), cmap=plt.cm.bone)

结果：

Answer 2

根据我对任务的理解，您需要编写一个光线跟踪器，跟踪 X 射线从源（这就是您需要其位置的原因）到投影平面（这就是您需要其位置的原因）。

边走边总结这些值，最后映射到允许的灰度值。

查看画线算法以了解如何执行此操作。

这真的不是什么黑魔法，我30多年前就做过这种事。妈的，我老了...

Answer 3

你要的是透视投影，而不是平行投影。为了获得这一点，您需要知道要为投影平面上的每个点求和哪些值。需要牢记多个注意事项：

我们谈论的是体素，因此您需要一种方法来确定 space 中的某个点是否属于您体积中的某个体素。
两点之间的线是直的，但是由于体素是 space 的离散表示，因此确定上述内容的不同方法可能会导致不同的（主要是次要的）结果。根据所使用的算法，这种差异最终也会导致图像略有不同。这是意料之中的。

假设您有一个由 256 个 512x512 像素切片组成的 CT 扫描体积。这将为您提供 512x512x256 体素的体积。对于这些体素中的每一个，您需要知道它们在 x、y、z 坐标中的位置。您可以按如下方式执行此操作：
- 使用 ImagePositionPatient 属性找出给定切片的左上角像素的 x、y、z 坐标（以毫米为单位）。
- 使用 PixelSpacing 属性计算切片中其他像素的 x、y、z 坐标。对所有切片重复

编辑：我刚刚找到了一个针对以下方法的反例，其余的仍然有帮助。会更新

~~现在要找出给定点 (Xa, Ya, Za) 如果源位于 (Xb, Yb, Zb)，需要对哪些体素值求和：~~

找到属于（Xa，Ya，Za）的体素。保留 pixel/voxel 数据。

计算（您可以使用 NumPy 计算）体素 (Xa, Ya, Za) 和 (Xb, Yb, Zb) 之间的距离。这里可以进行优化:)

对于所有直接环绕的体素（即 3x3x3-1 体素的数量）也计算此距离。还可以优化:)

以距离最短的体素为起点进行上述的下一次迭代。添加 pixel/voxel 数据。

~~重复直到超出您的 CT 体积范围。~~

为了获得投影，对投影平面上的所有点重复这些步骤并可视化结果。祝你任务顺利！ :)

从 3D DICOM 图像重建 2D X 射线

2D X-ray reconstruction from 3D DICOM images

python

projection

pydicom