2D X-ray reconstruction from 3D DICOM images

EDIT: as further answers noted, this solution yields a parallel projection, not a perspective projection.

From what I understand of the definition of "A normal 2D X-ray image", this can be done by summing each density for each pixel, for each slice of a projection in a given direction.

With your 3D volume, this means performing a sum over a given axis, which can be done with ndarray.sum(axis) in numpy.

# 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()

This yields the following result:

X-ray result

Which, to me, looks like a X-ray image.

EDIT : the result is a bit too "bright", so you may want to apply gamma correction. With matplotlib, import matplotlib.colors as colors and add a colors.PowerNorm(gamma_value) as the norm parameter in plt.imshow:

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

Result:

gamma corrected


The way I understand the task you are expected to write a ray-tracer that follows the X-rays from the source (that's why you need its position) to the projection plane (That's why you need its position).

Sum up the values as you go and do a mapping to the allowed grey-values in the end.

Take a look at line drawing algorithms to see how you can do this.

It is really no black magic, I have done this kind of stuff more than 30 years ago. Damn, I'm old...