CODEWITHSUNDEEP
pythonprojectionpydicom
AnaRhisT
AnaRhisT

Reputation: 121

2D X-ray reconstruction from 3D DICOM images

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

Input :

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:

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

Upvotes: 8

Views: 3046

Answers (3)

Renn Kane
Renn Kane

Reputation: 493

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

Upvotes: 7

g_uint
g_uint

Reputation: 2041

What you want is a perspective projection instead of a parallel projection. In order to obtain this, you need to know which values to sum for each point on the projection plane. There are multiple considerations to keep in mind:

  • We are talking about voxels, so you need to a method to determine whether a certain point in space belongs to a certain voxel in your volume.
  • A line between two points is straight, but because voxels are a discrete representation of space different methods of determining the above can lead to different (mostly minor) results. This difference will ultimately also lead to slightly different images depending on the alogrithms used. This is expected.

Let's say you have a CT scan volume comprising of 256 512x512 pixel slices. This gives you a volume of 512x512x256 voxels. For each of these voxels you need to know what their positions in x,y,z coordinates are. You can do this as follows:
- Use the ImagePositionPatient attribute to find out the x,y,z coordinate of the upper left hand corner pixel in mm for a given slice.
- Use the PixelSpacing attribute to calculate the x,y,z coordinates of the other pixels in your slice. Repeat for all slices

edit: i just found a counterexample against below method, the rest is still helpful. will update

Now to find out for a given point (Xa, Ya, Za) what voxel values need to be summed if the source is at (Xb, Yb, Zb):

  • Find the voxel that belongs to (Xa,Ya, Za). Keep pixel/voxel data.
  • Calculate (you can do this with NumPy) the distance between voxel(Xa, Ya, Za) and (Xb, Yb, Zb). There is an optimalization possible here :)
  • For all directly surrounding voxels (that will be a number of 3x3x3-1 voxels) also calculate this distance. Can also be optimized :)
  • Take the voxel with the shortest distance as the starting point for a next iteration of the above. Add pixel/voxel data.
  • Repeat until out of bounds of you CT volume.

In order to obtain a projection repeat these steps for all points on your projection plane and visualize the result. Good luck with your assignment! :)

Upvotes: 1

U. W.
U. W.

Reputation: 418

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...

Upvotes: 3

Related Questions