Search code examples
pythonprojectionpydicom

2D X-ray reconstruction from 3D DICOM images

I need to write a python function or class with the following Input/Output

Input :

  • The position of the X-rays source (still not sure why it's needed)
  • The position of the board (still not sure why it's needed)
  • A three dimensional CT-Scan

Output :

A 2D X-ray Scan (simulate an X-Ray Scan which is a scan that goes through the whole body)

A few important remarks to what I'm trying to achieve:

  • You don’t need additional information from the real world or any advanced knowledge.
  • You can add any input parameter that you see fit.
  • If your method produces artifacts, you are excepted to fix them.
  • Please explain every step of your method.

What I've done until now: (.py file added)

I've read the .dicom files, which are located in "Case2" folder.

These .dicom files can be downloaded from my Google Drive: https://drive.google.com/file/d/1lHoMJgj_8Dt62JaR2mMlK9FDnfkesH5F/view?usp=sharing

I've sorted the files by their position.

Finally, I've created a 3D array, and added all the images to that array in order to plot the results (you can see them in the added image) - which are slice of the CT Scans. (reference: https://pydicom.github.io/pydicom/stable/auto_examples/image_processing/reslice.html#sphx-glr-auto-examples-image-processing-reslice-py)

Here's the full code:

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

The result isn't what I wanted because:

These are slice of the CT scans. I need to simulate an X-Ray Scan which is a scan that goes through the whole body.

Would love your help to simulate an X-Ray scan that goes through the body.

I've read that it could be done in the following way: "A normal 2D X-ray image is a sum projection through the volume. Send parallel rays through the volume and add up the densities." Which I'm not sure how it's accomplished in code.

References that may help: https://pydicom.github.io/pydicom/stable/index.html

Solution

  • 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