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.

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.

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.

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.