python - 从3D DICOM图像进行2D X射线重建

我需要使用以下输入/输出编写python函数或类

输入:

X射线源的位置(仍然不确定为什么需要这样做)

板的位置(仍然不确定为什么需要这样做)

三维CT扫描

输出:

2D X射线扫描(模拟X射线扫描，它是遍及整个 body 的扫描)

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

您不需要来自现实世界的任何其他信息或任何高级知识。

您可以添加您认为合适的任何输入参数。

如果您的方法产生工件，则无法修复它们。

请说明您的方法的每个步骤。

到目前为止，我所做的事情:(添加了.py文件)

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

这些.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射线扫描，它是遍及整个 body 的扫描。

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

我已经读到它可以通过以下方式完成:“正常的2D X射线图像是通过体积的总投影。通过体积发送平行光线并累加密度。”我不确定它是如何在代码中完成的。

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

最佳答案

编辑:作为进一步的答案指出，此解决方案产生一个平行投影，而不是透视投影。

根据我对“正常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并在colors.PowerNorm(gamma_value)中添加norm作为plt.imshow参数:

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

结果:

关于python - 从3D DICOM图像进行2D X射线重建，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/57497695/