ASTRA Toolbox Tutorial: Reconstruction from Projection Images, Part 2

This is part two 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 part one, I’ve created a synthetic dataset. In this part, I’ll use that dataset to create a reconstruction. Note that, if you start from a real dataset instead of from this simulated one, you might still have to convert your X-ray images into projection images. This is explained in part 1 of my previous series of articles on tomography. However, apart from that, these projection images are completely equivalent to real ones, albeit that they look as if taken with a mechanically perfect (but not noiseless!) scanner.

With the risk of repeating a few things from part one, I’ll list all steps to take to go from the projection images to the final reconstruction.

Scan Parameters

The first thing to do is to gather all relevant information about the scan. This information will be in some sort of logfile that comes with the scan, since the parameters can be adjusted for each scan, even if the same scanner is used. For example, the position of the rotating stage with the sample can be shifted closer to or farther away from the detector, to change the magnification. For the purpose of this article, I’ll use the Python fragment from part one of this tutorial, which is repeated below.

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)

This “logfile” provides enough information to be able to create a reconstruction from the projection images. Of course, you also have to be sure that the projections were actually taken with a cone-beam scanner.

The origin of the coordinate system of the scanner is the center of the object (assuming that the object is in the center of the rotating stage, of course). From the Python code, you can see that the distance between the source and the origin, let’s call that \(d_\mathrm{so}\), is 300 mm, and that the distance between the origin and the detector, \(d_\mathrm{od}\), is 100 mm. This means that the magnification is \((d_\mathrm{so}+d_\mathrm{sd})/d_\mathrm{so}=4/3\) in this case. The detector has square pixels with size \(s_\mathrm{det}=1.05\,\mathrm{mm}\), and has 200 rows in the vertical direction and 200 columns in the horizontal direction. The angles are equally spaced in 180 steps over 360 degrees, which means that a projection was made every two degrees.

You should be able to adapt these values by looking at the logfiles of your own scans. Sometimes the distances are given as source–origin and source–detector, so you might to have to adapt for that. In the Python code that follows, I assume that your projections are in a directory called dataset and that your files are named proj####.tif.

Loading the Dataset

Loading the dataset is done in two steps. The first step is reading the dataset from disk and putting it in a 3D NumPy array, as follows.

projections = np.zeros((detector_rows, num_of_projections, detector_cols))
for i in range(num_of_projections):
    im = imread(join(input_dir, 'proj%04d.tif' % i)).astype(float)
    im /= 65535
    projections[:, i, :] = im

This is straighforward, but note that it is axis one that increases with each projection. The code assumes 16-bit TIFF files. This is quite common, but you might have to adapt this if your projection images are in a different format. I assume here that the projections use the full range of 0 to 65535. This wil normally never be the case with actual data, so you will have to rescale your data. The actual maximum is not important, because we are working with relative numbers here. I have rescaled for a maximum of one in this example. The minimum, however, should be zero after you have converted the X-ray images into projection images (again, see my existing series of articles on tomography). This is because, in places where there is no object in between the source and the detector, the measured density of the (non-existant) material should be exactly zero.

Creating the Projection Geometry

The second step is then to put the NumPy variable projections into the toolbox by using the function astra.data3d.create(). To do that, you first have to describe the geometry of your projection data by creating a projection geometry through calling astra.create_proj_geom().

The creation of the projection geometry is one of the crucial steps that you need to understand well to successfully apply the ASTRA Toolbox. There is one important property of the geometries of the ASTRA Toolbox, which is that the voxels of the reconstruction area have a size of 1 (so 1×1×1 in the case of the 3D reconstruction that we are working on). This implies that you have to scale the rest of your parameters to make this so. A little “trick” to make this easier to handle is to put the detector virtually in the middle of the object. This is not mandatory, but it simplifies thinking about it. Or, at least, that’s what I think.

The parameters of the scanner as given in the Python code fragment above are the following.

\[\begin{align}
d_\mathrm{so}&=300\,\mathrm{mm}\\
d_\mathrm{od}&=100\,\mathrm{mm}\\
s_\mathrm{det}&=1.05\,\mathrm{mm}
\end{align}\]

Shifting the detector to the origin amounts to setting the origin–detector distance to zero and rescaling the size of the detector pixels to compensate for that. This results in a new set of parameters:

\[\begin{align}
d_\mathrm{so}'&=d_\mathrm{so}=300\,\mathrm{mm}\\
d_\mathrm{od}'&=0\,\mathrm{mm}\\
s_\mathrm{det}'&=\frac{d_\mathrm{so}}{d_\mathrm{so}+d_\mathrm{od}}s_\mathrm{det}\approx 0.79\,\mathrm{mm}
\end{align}\]

After shifting the detector to the origin, the parameter \(s_\mathrm{det}'\) gives you the actual size of the reconstruction voxels in real-world units. In the example scanner, each voxel in the reconstructed object will correspond to a cube of 0.79×0.79×0.79 mm. This is important information if you need to know or measure the actual dimensions of the scanned object in real-world units.

The next step is then to scale the parameters so that the reconstruction voxels get a size of 1×1×1. This results in the final set of parameters:

\[\begin{align}
d_\mathrm{so}''&=\frac{d_\mathrm{so}'}{s_\mathrm{det}'}=\frac{d_\mathrm{so}+d_\mathrm{od}}{s_\mathrm{det}}\approx 381\\
d_\mathrm{od}''&=0\\
s_\mathrm{det}''&=1
\end{align}\]

I’ve left out the units in this last set of expressions, since it seems confusing to keep including them after this scaling. In Python, this results in the following code.

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 = astra.data3d.create('-sino', proj_geom, projections)

The result of this step is that you now have a projections_id handle to the projection data in the Toolbox that you need to provide to subsequent calls to the Toolbox.

Reconstruction

The reconstruction itself is normally not the difficult part. You create a volume geometry and then create a 3D data object inside of the Toolbox to hold the reconstruction. This provides you with a reconstruction_id handle. You then choose an algorithm and specify its parameters. For a cone-beam dataset, the FDK_CUDA algorithm is the obvious one to start with. The resulting Python code is shown below.

vol_geom = astra.creators.create_vol_geom(detector_cols, detector_cols,
                                          detector_rows)
reconstruction_id = astra.data3d.create('-vol', vol_geom, data=0)
alg_cfg = astra.astra_dict('FDK_CUDA')
alg_cfg['ProjectionDataId'] = projections_id
alg_cfg['ReconstructionDataId'] = reconstruction_id
algorithm_id = astra.algorithm.create(alg_cfg)
astra.algorithm.run(algorithm_id)
reconstruction = astra.data3d.get(reconstruction_id)

After this part of the code has run, the NumPy array reconstruction contains the 3D reconstructed object. You could then cut off negative values, which should be due to noise if the scaling of your dataset was done correctly. To save the reconstruction as a series of individual slices, you can convert the reconstruction volume to an 8-bit (or 16-bit) unsigned int and save the slices as PNG or TIFF files:

reconstruction[reconstruction < 0] = 0
reconstruction /= np.max(reconstruction)
reconstruction = np.round(reconstruction * 255).astype(np.uint8)
 
for i in range(detector_rows):
    im = reconstruction[i, :, :]
    im = np.flipud(im)
    imwrite(join(output_dir, 'reco%04d.png' % i), im)

Note the essential call to flipud(). This necessary because of the inverted Y-axis that you get when you convert from a NumPy array to an image. The resulting PNG for the so-called central slice, which is slice 100 in this case, is shown in Figure 1.

Figure 1. Central slice of the reconstructed object.

I hope that this extensive tutorial article helps you in applying the ASTRA Toolbox to your own cone-beam CT datasets!

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 imread, imwrite
 
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)
input_dir = 'dataset'
output_dir = 'reconstruction'
 
# Load projections.
projections = np.zeros((detector_rows, num_of_projections, detector_cols))
for i in range(num_of_projections):
    im = imread(join(input_dir, 'proj%04d.tif' % i)).astype(float)
    im /= 65535
    projections[:, i, :] = im
 
# Copy projection images into ASTRA Toolbox.
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 = astra.data3d.create('-sino', proj_geom, projections)
 
# Create reconstruction.
vol_geom = astra.creators.create_vol_geom(detector_cols, detector_cols,
                                          detector_rows)
reconstruction_id = astra.data3d.create('-vol', vol_geom, data=0)
alg_cfg = astra.astra_dict('FDK_CUDA')
alg_cfg['ProjectionDataId'] = projections_id
alg_cfg['ReconstructionDataId'] = reconstruction_id
algorithm_id = astra.algorithm.create(alg_cfg)
astra.algorithm.run(algorithm_id)
reconstruction = astra.data3d.get(reconstruction_id)
 
# Limit and scale reconstruction.
reconstruction[reconstruction < 0] = 0
reconstruction /= np.max(reconstruction)
reconstruction = np.round(reconstruction * 255).astype(np.uint8)
 
# Save reconstruction.
if not isdir(output_dir):
    mkdir(output_dir)
for i in range(detector_rows):
    im = reconstruction[i, :, :]
    im = np.flipud(im)
    imwrite(join(output_dir, 'reco%04d.png' % i), im)
 
# Cleanup.
astra.algorithm.delete(algorithm_id)
astra.data3d.delete(reconstruction_id)
astra.data3d.delete(projections_id)