ASTRA Toolbox Tutorial: Reconstruction from Projection Images, Part 1

This is part one of a tutorial that describes how to use the ASTRA Toolbox to create a 3D reconstruction from 2D projection images that were taken with a cone-beam CT scanner. In this first part, I’ll create a synthetic test dataset. Part two then starts from this synthetic dataset and creates a reconstruction from it. The motivation for writing this tutorial were several questions on my article on the ASTRA Toolbox. I hope that this tutorial helps you to use the ASTRA Toolbox for your cone-beam reconstructions.

Virtual Scanner

The main difference between the typical reconstruction examples from phantom images and real projection images is that you don’t need to fiddle with the parameters of an actual scanner. However, this also seems to imply that many people have trouble applying the examples that are included with the Toolbox to actual data. In this tutorial, I’ll give “real-world” dimensions of a scanner, albeit a virtual one, so that you can subsitute the given values with the ones from your own scanner. My virtual scanner has the following properties.

  • The X-ray source is a single point, and the detector is a flat panel. The object is on a rotating stage between both. This is a typical cone-beam setup.
  • The origin of the coordinate system of the scanner is the center of the rotating stage.
  • The distance between the X-ray source and the origin, \(d_\mathrm{so}\), is 300 mm.
  • The distance between the origin and the detector, \(d_\mathrm{od}\), is 100 mm.
  • The pixel size of the detector is 1.05 mm in both directions.
  • The detector has 200 rows (vertical direction) and 200 columns (horizontal direction).

Together with the information on the projection angles that I’ll use for the virtual scan, this translates into the following assignment to Python variables.

distance_source_origin = 300  # [mm]
distance_origin_detector = 100  # [mm]
detector_pixel_size = 1.05  # [mm]
detector_rows = 200  # Vertical size of detector [pixels].
detector_cols = 200  # Horizontal size of detector [pixels].
num_of_projections = 180
angles = np.linspace(0, 2 * np.pi, num=num_of_projections, endpoint=False)

The number of projections for the virtual scan is set to 180. Spread out evenly over 360 degrees, this implies an angular step of 2 degrees. Note that you have to use the endpoint=False option of linspace() to avoid setting the last projection to 360 degrees and using the wrong angular step of ≈2.11 degrees.

Test Phantom

As a simple phantom, I’ve created a hollow beam, with a small hole in one of the sides to break the rotational symmetry (Figure 1).

Figure 1. Hollow beam.Figure 1. Hollow beam.

The following code fragment creates a 3D NumPy array with this object and loads it into the ASTRA Toolbox.

vol_geom = astra.creators.create_vol_geom(detector_cols, detector_cols,
                                          detector_rows)
phantom = np.zeros((detector_rows, detector_cols, detector_cols))
hb = 110  # Height of beam [pixels].
wb = 40   # Width of beam [pixels].
hc = 100  # Height of cavity in beam [pixels].
wc = 30   # Width of cavity in beam [pixels].
phantom[detector_rows // 2 - hb // 2 : detector_rows // 2 + hb // 2,
        detector_cols // 2 - wb // 2 : detector_cols // 2 + wb // 2,
        detector_cols // 2 - wb // 2 : detector_cols // 2 + wb // 2] = 1
phantom[detector_rows // 2 - hc // 2 : detector_rows // 2 + hc // 2,
        detector_cols // 2 - wc // 2 : detector_cols // 2 + wc // 2,
        detector_cols // 2 - wc // 2 : detector_cols // 2 + wc // 2] = 0
phantom[detector_rows // 2 - 5 :       detector_rows // 2 + 5,
        detector_cols // 2 + wc // 2 : detector_cols // 2 + wb // 2,
        detector_cols // 2 - 5 :       detector_cols // 2 + 5] = 0
phantom_id = astra.data3d.create('-vol', vol_geom, data=phantom)

With the configuration of the setup given above, these pixel sizes give the example object a physical size of (about) 8.7 cm × 3.2 cm × 3.2 cm. The calculation of the size of a phantom or reconstruction voxel is in part two of this tutorial. The apparent size of the object in the projection images will be larger by a factor of \((d_\mathrm{so}+d_\mathrm{od})/d_\mathrm{so}=4/3\), because of the magnifying effect of the cone-beam projection.

In the virtual scanner, the object will rotate around the vertical axis, which is axis zero of the NumPy array. If you put the object in the scanner with the hole towards the detector, and you plot the so-called central slice, which is the horizontal cut of the object that is projected onto the middle row of the detector, then this looks like Figure 2.

Figure 2. Central slice of hollow beam.Figure 2. Central slice of hollow beam.

With this slice plotted as in Figure 2, the projection from angle zero looks upwards from the bottom of the slice, with the source below the image and the detector above it. This means that the hole in the 3D object is in the side closest to the detector.

Projections

The next step is then to create the projections, i.e., a dataset as it would be recorded by a real scanner. This is achieved by the following code fragment.

proj_geom = \
  astra.create_proj_geom('cone', 1, 1, detector_rows, detector_cols, angles,
                         (distance_source_origin + distance_origin_detector) /
                         detector_pixel_size, 0)
projections_id, projections = \
  astra.creators.create_sino3d_gpu(phantom_id, proj_geom, vol_geom)
projections /= np.max(projections)

How to get to the given values for the parameters of the call to astra.create_proj_geom is explained in part two. This same geometry is also used in the reconstruction process, and it is by far the most confusing part of it, so I wanted to put it in the part of the tutorial that deals with the reconstruction.

Three of the projections are shown in Figure 3, namely the ones at 0, 30, and 90 degrees. The projections are taken from increasing angles, so in a counterclockwise fashion. This is equivalent to the object rotating clockwise, so that the hole that was at the back of the object in the projection at 0 degrees is on the right of the image at 90 degrees.

Figure 3. Projections at 0, 30, and 90 degrees.Figure 3. Projections at 0, 30, and 90 degrees.

The last part of the Python program then applies Poisson noise and saves the projections to files (see below for the full program). I specifically picked 16-bit TIFF for the format, since many scanners produce data in that format. Of course, other formats are possible, although 8-bit ones should probably be avoided.

Complete Python Program

The complete Python program follows below. Note that you need an NVIDIA GPU to run this code, since the 3D algorithms of the ASTRA Toolbox are implemented as CUDA GPU code.

from __future__ import division
 
import numpy as np
from os import mkdir
from os.path import join, isdir
from imageio import get_writer
 
import astra
 
# Configuration.
distance_source_origin = 300  # [mm]
distance_origin_detector = 100  # [mm]
detector_pixel_size = 1.05  # [mm]
detector_rows = 200  # Vertical size of detector [pixels].
detector_cols = 200  # Horizontal size of detector [pixels].
num_of_projections = 180
angles = np.linspace(0, 2 * np.pi, num=num_of_projections, endpoint=False)
output_dir = 'dataset'
 
# Create phantom.
vol_geom = astra.creators.create_vol_geom(detector_cols, detector_cols,
                                          detector_rows)
phantom = np.zeros((detector_rows, detector_cols, detector_cols))
hb = 110  # Height of beam [pixels].
wb = 40   # Width of beam [pixels].
hc = 100  # Height of cavity in beam [pixels].
wc = 30   # Width of cavity in beam [pixels].
phantom[detector_rows // 2 - hb // 2 : detector_rows // 2 + hb // 2,
        detector_cols // 2 - wb // 2 : detector_cols // 2 + wb // 2,
        detector_cols // 2 - wb // 2 : detector_cols // 2 + wb // 2] = 1
phantom[detector_rows // 2 - hc // 2 : detector_rows // 2 + hc // 2,
        detector_cols // 2 - wc // 2 : detector_cols // 2 + wc // 2,
        detector_cols // 2 - wc // 2 : detector_cols // 2 + wc // 2] = 0
phantom[detector_rows // 2 - 5 :       detector_rows // 2 + 5,
        detector_cols // 2 + wc // 2 : detector_cols // 2 + wb // 2,
        detector_cols // 2 - 5 :       detector_cols // 2 + 5] = 0
phantom_id = astra.data3d.create('-vol', vol_geom, data=phantom)
 
# Create projections. With increasing angles, the projection are such that the
# object is rotated clockwise. Slice zero is at the top of the object. The
# projection from angle zero looks upwards from the bottom of the slice.
proj_geom = \
  astra.create_proj_geom('cone', 1, 1, detector_rows, detector_cols, angles,
                         (distance_source_origin + distance_origin_detector) /
                         detector_pixel_size, 0)
projections_id, projections = \
  astra.creators.create_sino3d_gpu(phantom_id, proj_geom, vol_geom)
projections /= np.max(projections)
 
# Apply Poisson noise.
projections = np.random.poisson(projections * 10000) / 10000
projections[projections > 1.1] = 1.1
projections /= 1.1
 
# Save projections.
if not isdir(output_dir):
    mkdir(output_dir)
projections = np.round(projections * 65535).astype(np.uint16)
for i in range(num_of_projections):
    projection = projections[:, i, :]
    with get_writer(join(output_dir, 'proj%04d.tif' %i)) as writer:
        writer.append_data(projection, {'compress': 9})
 
# Cleanup.
astra.data3d.delete(projections_id)
astra.data3d.delete(phantom_id)

Part two of this tutorial creates a reconstruction from the synthetic dataset that is created by the above Python script.

Zhengchun Liu (not verified)

Sat, 10/05/2019 - 00:56

could you create an example that uses an existing sinogram for reconstruction using SIRT? instead of simulate one

I considered that when I wrote the article, but in the end I didn't because these datasets are typically very large (and everyone who wants to run the code would then have to download it).

gideon (not verified)

Fri, 10/11/2019 - 11:30

thank you for the source

Howard (not verified)

Tue, 12/31/2019 - 17:37

Is it possible to place a beam modifier, say, a random binary mask, between the source and the imaging object/phantom? I do not see any examples for this application anywhere. Thanks.

Happy New Year!

Masks are supported by the reconstruction algorithms. See, for example, some general info or the SIRT algorithm (the SinogramMaskId parameter). Masks are also supported by the Python API. I’ve used masks to implement PDART, for example, using SIRT as the “base algorithm” and a mask to remove the dense pixels from the reconstruction.

Howard (not verified)

Fri, 01/03/2020 - 16:19

Thank you for your reply, Tom. I am actually looking for the ability to place the beam modifier between the X-ray source and the imaging object. This beam modifier can be anything like a thin disk to attenuate the beam intensity, say, random thickness disk, not the mask for image reconstruction. I am more interested to generate projection images using this type of beam modifier. I do not see any examples/tutorial on how to place it between X-ray source and the phantom. Is it possible with ASTRA at all? Thanks!

I see what you mean, and I’m afraid that that option is not foreseen. However, given that the ASTRA Toolbox assumes monochromatic rays, the effect of a beam modifier will be something that you can emulate by applying a post-processing step to the projection images, I think. The Toolbox does not support emulating polychromatic (aka realistic) X-ray sources, so the real effect of, say, an aluminum filter, cannot be produced easily. Sadly, using filters probably implies that you want to see the effect of polychromatic rays, so you may be out of luck… If you’re trying to do something like Compressed Sensing, then a post-processing step might be enough…

Thanks Tom for taking the time to answer my question and make suggestions. The reason I am interested in ASTRA for my problem was inspired by the work presented in the paper entitled "Source and coded aperture joint optimization for compressive X-ray tomosynthesis". Fig. 1 shows clearly how they utilized the coded aperture which was placed between the X-ray source and imaging object to better optimize the reconstruction of tomosynthesis. Section 4 on page 10 of their paper described the simulation. I contacted the authors but never received their response. It appears they could manage to perform the simulation using ASTRA with some sort of tricks.

I didn't know that paper, but it actually is about compressed sensing ☺; I didn't know that people were doing that for tomosynthesis or tomography...

I listed the paper to illustrate that the authors managed to simulate what I have been looking for, i.e., place the coded aperture or whatever beam modifier between the X-ray source and the phantom. They used ASTRA so there ought to be a way to implement it. The compressed sensing technique is irrelevant to this issue.

Manju (not verified)

Fri, 03/27/2020 - 11:15

Do you have the virtual scanning codes in MATLAB? I am not familiar with python and so couldn't follow the codes

Sorry, I don't have the MATLAB equivalent of the above code in a readily presentable way…

S Bala (not verified)

Fri, 07/17/2020 - 19:20

Hi, I am new to python. I have difficulty in installing the astra packages using conda install -c astra-toolbox astra-toolbox . Can you help me out?

G.B. (not verified)

Wed, 06/16/2021 - 20:13

Hey, I didn't understand how you directly add the poission noise into the projections? Couldyou give me some details? Many thanks.
# Apply Poisson noise.
projections = np.random.poisson(projections * 10000) / 10000
projections[projections > 1.1] = 1.1
projections /= 1.1

Normally, noise should be added like this:
filename = 'myimage.png'
imagea = (scipy.misc.imread(filename)).astype(float)
poissonNoise = numpy.random.poisson(imagea).astype(float)
noisyImage = imagea + poissonNoise

The code that you show is not correct. You should not create a separate array and add it. This is because Gaussian Noise is Added, Poisson Noise is Applied. Because the projections are in the range [0, 1], the expression np.random.poisson(projections * 10000) / 10000 retuns a noisy image that simulates Poisson noise as if the brightest pixels were hit by, on average, 10000 photons. For example, if a certain pixel has the value 5000, then, for that particular pixel, a value is drawn out of a Poisson distribution with a mean of 5000. This means that the output of np.random.poisson() is immediately the required noisy image.

Anthony Whittemore (not verified)

Tue, 11/16/2021 - 19:27

Tom (or anyone with experience),

Is there a way to visualize your experimental setup? Also, when developing CT simulations what is your typical workflow when using ASTRA?

Thanks!

zhihan (not verified)

Tue, 06/07/2022 - 11:43

Hi, I don't understand how to let astra know the pixel size of input objet. In the case given above, "these pixel sizes give the example object a physical size of (about) 8.7 cm × 3.2 cm × 3.2 cm", how are these numbers calculated? Many thanks.

This is explained in part two of this tutorial, where I determine that the geometry that I picked leads to a pixel size of about 0.79 mm, so that my virtual object of 110×40×40 pixels has a physical size of about 8.7×3.2×3.2 cm.

Hu Zhou (not verified)

Fri, 10/07/2022 - 21:32

May I ask how is the phantom specified? Is each voxel described as a material with density, or as the linear attenuation coefficient (in 1/cm or HU)? Any example?
Thanks.
Hu

Each voxel has either density zero or density one, without referring to any real-world units. So this would be "material with density", from the two options that you give.

Hi Tom,
Thank you very much for this interesting and helpful tutorial.

I am working on PET image reconstruction topic and working with real scanner (Siemens Biograph mCT 1104 scanner) for our clinical data. Since, I want to use GPU accelerated image reconstruction, I need to setup our real scanner geometry in to the ASTRA package. So, my question is that how I can do this. Is there any example to setup a real scanner in to the ASTRA.

Your help and guidance is really appreciated.
Thanks.

Add new comment

The content of this field is kept private and will not be shown publicly.
Spam avoidance measure, sorry for this.

Restricted HTML

  • Allowed HTML tags: <a href hreflang> <em> <strong> <cite> <blockquote cite> <code> <ul type> <ol start type> <li> <dl> <dt> <dd> <h2 id> <h3 id> <h4 id> <h5 id> <h6 id>
  • Lines and paragraphs break automatically.
  • Web page addresses and email addresses turn into links automatically.
Submitted on 3 October 2018