Reconstruction of a Real Dataset, Part 4: The Actual Reconstruction

This is part four of a tutorial that describes how to create a reconstruction of a cone-beam computed tomography (CT) dataset (also read part one, part two, and part three). In this part, I’ll show how to create the actual reconstruction from the projections that were created in part three. To keep you in suspense no longer, Figure 1 shows the central slice of the reconstruction.

Figure 1. Central slice of the reconstruction.

Python Code

The actual reconstruction code is very close to the code from ASTRA Toolbox Tutorial: Reconstruction from Projection Images, Part 2, so I won’t go into a lot of detail here, and go straight to the actual script! A few remarks follow.

This script imports the file config.py from part one.
As for the code in ASTRA Toolbox Tutorial: Reconstruction from Projection Images, Part 2, you need an NVIDIA GPU to run this code, since the 3D algorithms of the ASTRA Toolbox are implemented as CUDA GPU code (for the previous parts of the current tutorial, a GPU or the ASTRA Toolbox were not required yet).
As is often done, I have put the detector virtually in the origin of the coordinate system, so that it is in the center of the object. This means that the size of the detector pixels has to be adjusted, which is done in the computation of detector_pixel_size_in_origin.
A funny statement is im = np.roll(im, cfg.horizontal_shift, axis=1). This shifts each projection image two pixels to the left (the value of cfg.horizontal_shift is −2), to compensate for the axis of rotation not being entirely in the middle of the detector. This implementation is not entirely kosher, because the roll operation adds the two leftmost columns of the image to the right of the output image. This doesn’t make sense, but I can get away with it because there is only background in those columns. The ASTRA Toolbox allows you to do a direct adjustment of the geometry, which then also supports shifting by a fraction of a pixel, but I’ll keep that for a later post.
There is some postprocessing after the reconstruction. First I zero negative pixels. Then I compute the maximum using the central 200 slices. As for the code in part part three, the maximum is needed to create nicely scaled TIFF files. I use only the central 200 slices to avoid some bright spots in slices further away from the middle.

#filename: reconstruct.py
from __future__ import print_function
from __future__ import division
 
from os import makedirs
from os.path import join, exists
import numpy as np
from imageio import imread, get_writer
 
import astra
import config as cfg
 
def save_tiff(name, im_in):
    # Expects image as floating point in range [0, 1].
    im = im_in.copy()
    im = np.round(im * 65535).astype(np.uint16)
    with get_writer(name) as writer:
        writer.append_data(im, {'compress': 0})
 
if not exists(cfg.reco_dir):
    makedirs(cfg.reco_dir)
 
detector_pixel_size_in_origin = \
  cfg.detector_pixel_size * cfg.distance_source_origin / \
  (cfg.distance_source_origin + cfg.distance_origin_detector)
 
# Determine dimensions of projection images.
im = imread(join(cfg.proj_dir, 'proj0000.tif'))
dims = im.shape
# dims[0]: Number of rows in the projection image, i.e., the height of the
#          detector. This is Y in the Cartesian coordinate system.
# dims[1]: Number of columns in the projection image, i.e., the width of the
#          detector. This is X in the Cartesian coordinate system.
detector_rows = dims[0]
detector_columns = dims[1]
 
# Load projection images.
projections = np.zeros((detector_rows, cfg.num_projections, detector_columns))
for proj in range(cfg.num_projections):
    im = imread(join(cfg.proj_dir, 'proj{:04d}.tif'.format(proj))).astype(float)
    im /= 65535
    im = np.roll(im, cfg.horizontal_shift, axis=1)
    projections[:, proj, :] = im
 
# Copy projection images into ASTRA Toolbox.
proj_geom = astra.create_proj_geom('cone', 1, 1, detector_rows,
                                   detector_columns, cfg.angles,
                                   cfg.distance_source_origin /
                                   detector_pixel_size_in_origin, 0)
projections_id = astra.data3d.create('-sino', proj_geom, projections)
 
# Create reconstruction.
vol_geom = astra.creators.create_vol_geom(detector_columns, detector_columns,
                                          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)
print('Reconstructing... ', end='')
astra.algorithm.run(algorithm_id)
print('done')
reconstruction = astra.data3d.get(reconstruction_id)
 
# Limit and scale reconstruction.
reconstruction[reconstruction < 0] = 0
middle = detector_rows // 2
M = np.max(reconstruction[middle - 100 : middle + 100, :, :])
reconstruction /= M
reconstruction[reconstruction > 1] = 1
 
# Save reconstruction.
for proj in range(detector_rows):
    print('Saving slice %d' % proj)
    slice = reconstruction[proj, :, :]
    save_tiff(join(cfg.reco_dir, 'reco%04d.tif' % proj), slice)
 
# Cleanup.
astra.algorithm.delete(algorithm_id)
astra.data3d.delete(reconstruction_id)
astra.data3d.delete(projections_id)

Download the Apple Dataset

If you think it would help you in understanding the code and applying it to your own datasets, I can actually provide you with the dataset that I’ve used in this series of articles! The nice people at UGCT, who created this dataset, have given me permission to make it available for download.

Simply shoot me a message and I’ll provide you with a download link (I don’t dare to put an actual download link here, since that might use a lot of bandwidth if it is linked from somewhere popular).