三维DICOM图像的二维x射线重建

Question

我需要用下面的输入/输出编写一个python函数或类

输入：

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

输出：

二维X射线扫描（模拟X射线扫描，即扫描整个身体）

关于我要实现的目标，我要说几句重要的话：

• 你不需要来自真实世界的额外信息或任何先进的知识。在
• 您可以添加任何您认为合适的输入参数。在
• 如果您的方法产生了工件，那么您必须修复它们。在
• 请解释你方法的每一步。在

到目前为止我所做的事情：（.py文件已添加）

我已经阅读了.dicom文件，它们位于“Case2”文件夹中。在

这些.dicom文件可以从我的谷歌硬盘下载： 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光扫描。在

我读过这样的书：“一个普通的2dx射线图像是通过体积的和投影。把平行光穿过体积，再把密度加起来。”我不知道这是如何用代码实现的。在

可能有帮助的引用：https://pydicom.github.io/pydicom/stable/index.html

Answer 1

我理解这个任务的方式，你需要写一个光线跟踪器，它跟随X射线从源（这就是为什么你需要它的位置）到投影平面（这就是为什么你需要它的位置）。在

总结你的价值，并做一个映射到允许的灰色值在最后。在

看看画线算法，看看你是如何做到这一点的。在

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

Answer 2

你想要的是透视投影而不是平行投影。要获得此值，需要知道投影平面上每个点的和值。有多个注意事项需要记住：

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

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

编辑：我刚刚找到了一个反例来反对下面的方法，其余的仍然有用。将更新

现在要找出给定点（Xa，Ya，Za）在（Xb，Yb，Zb）处需要求和哪些体素值：

• 找到属于（Xa，Ya，Za）的体素。保留像素/体素数据。在
• 计算（你可以用NumPy）体素（Xa，Ya，Za）和（Xb，Yb，Zb）之间的距离。这里有一个可能的优化：）
• 对于所有直接围绕的体素（即3x3x3-1体素）也计算此距离。也可以优化：）
• 以距离最短的体素作为上一次迭代的起点。添加像素/体素数据。在
• 重复，直到超出CT容积的范围。在

为了获得投影，对投影平面上的所有点重复这些步骤，并将结果可视化。祝你的作业顺利！：）

Answer 3

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

根据我对“普通2D X射线图像”定义的理解，这可以通过求和给定方向上投影的每个像素的每个密度来实现。在

对于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光照片。在

编辑：结果有点太“亮”，所以您可能需要应用gamma校正。使用matplotlib，import matplotlib.colors as colors并在plt.imshow中添加一个colors.PowerNorm(gamma_value)作为norm参数：

^{pr2}$

结果：