8

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

3
  • 3
    Can you add the expect result and the actual results you got? Aug 18, 2019 at 12:26
  • 1
    @TarunLalwani The results I got are just the Renn's results which are close to the solution.But instead of using the Renn's solution to every axis (X,Y,Z), the expected result is somehow to use a projection on any point in space (x,y,z) in order to see the X-RAY as a 2D Image. and not 3 different slices. Aug 18, 2019 at 15:52
  • "A normal 2D X-ray image is a sum projection through the volume. Send parallel rays through the volume and add up the densities." is not correct in this context. The difference is that that implies a source for every pixel, and their rays travelling parallel to each other. In fact an x ray has a point for a source, and all rays eminate from that point, so a summation over on of the axis will not suffice. That is why your results are slightly different. And that is why you need the position of the source.
    – g_uint
    Aug 23, 2019 at 7:32

3 Answers 3

7
+150

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

3
  • 2
    Hi Renn! Thanks for you approach. I'm not quite sure this is the desired result. I've got this result with my code too, just without summing. I'm trying to see how the end-result should look like from the one who sent me the exercise, I'l update you when he answers.
    – AnaRhisT
    Aug 14, 2019 at 16:13
  • Hi Renn, how can it be upgraded to every point in space and not in X,Y,Z direction? Any ideas? Aug 18, 2019 at 9:28
  • 1
    @RennKane Thanks for your help, although I did successfully figured out to sum on all the axis, but the additional gamma correction was awesome. They accepted my solution :) Aug 26, 2019 at 20:39
3

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

4
  • Like Renn said? Aug 22, 2019 at 15:50
  • Yes, like he said. But the task looks like an exercise, so I don't want to take the chance away from him to learn something valuable! :-)
    – U. W.
    Aug 23, 2019 at 6:55
  • 1
    Not quite what Renn said, the difference is subtle but exactly what @U.W. says. In other words, Renn's method is as if having a source for every pixel, and all rays travelling exactly parallel to each other (along the axis to which is summed), but what is asked is to sum over the direction from the source, which is a point, not a magic xray slab that is just as large as the captured area ;)
    – g_uint
    Aug 23, 2019 at 7:25
  • @g_uint: Exactly. It is the difference between parallel and perspective projection. Have fun implementing it!
    – U. W.
    Aug 23, 2019 at 7:51
1

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! :)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.