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.

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.path import join 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. 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)

Hi,

Thank you very much for the example above. I have used it to test in my datasets successfully in the way that I managed to get the reconstruction. However some more issues are getting me stuck.

First, even having 8 GB of memory in the GPU, I get the following error for datasets of about 300 MB:

Error: CUDA error 2: out of memory.

Error: Failed to allocate 672x356x768 GPU buffer

Error: CUDA error 11: invalid argument.

The other point are:

How do I control some other parameters, such as the center of rotation, number of iterations, output histogram limits?

I will appreciate your help ...

Kind regards,

The errors that you get seem to point to insufficient memory, even though your dataset is only 300 MB. You have to take into account that a full reconstruction volume (typically n×n×n floats for n×n projection images) has to fit in the memory of the GPU card, in addition to your dataset that takes up a volume of n×n×p floats, with p the number of projections (even if your images are, e.g., 16-bit compressed TIFFs). I would start by downscaling the projection images and creating a smaller reconstruction first. Another option is to crop the projection images, although it is essential that you do not cut off parts of the object in any of the projections.

Hi Tom,

Thank you for this tutorial.

I have the same problem as Liebert. I want to use the code for a real CBCT case: n projections of 1840x1456 pixels. In my first test, the reconstruction volume is defined as 438x438x308 voxels. When n is less than 49, everything works. But when n = 50 or more, the call to astra_mex_algorithm_run() (when called astra_mex_algorithm('iterate', alg_id, this.IterNum)) crashes my matlab session with the following system error: "the display pilot was no longer responding and was recovered".

Looking at the GPU RAM consumption, I see that the crash occurs when the GPU RAM limit is reached (3GB) and 37 MB is used for each projection. So, for a complete projection of 1400 images, I will need a 51GB GPU RAM.

Does the Astra toolbox handle the usual CBCT cases, or is it limited to small datasets?

Thanks in advance for your answer.

Kind regards,

Your other questions are more diverse and more difficult to answer. To correct the center of rotation, have a look at the

`cone_vec`

geometry. The`FDK_CUDA`

algorithm of the example is not iterative, so there is no number of iterations. See the documentation for algorithms such as`SIRT_CUDA`

that are iterative.This is a wonderful tutorial and thank you very much for sharing.

I notice that these are Python codes and I'm using MATLAB to reconstruct 3D images. So could you please tell me how to obtain the projection id in MATLAB? I mean, in Python we can use "projections_id = astra.data3d.create('-sino', proj_geom, projections)", however I used "proj_id=astra_mex_data3d('create', '-proj3d', proj_geom, proj_data);" in MATLAB and got a wrong reconstruction result.

Thank you in advance!

I've solved this problem. I used "proj_id=astra_mex_data3d('create', '-sino', proj_geom, proj_data);" and got a correct result finally. It seems that changing '-sino' to '-proj3d' doesn't influence the result. The reason I got the wrong result before was that the angle of my observation is not correct. I mean, it was the correct result indeed. Thank you again for sharing!

Thank you very much Tom! It's a very useful tutorial, and amazing website! I've learned a lot from your fabulous work.

I have a similar problem as Liebert, it's also a memory issue. I just want to make sure, if my x-ray image is 2000*2000, is it necessary that my volume matrix must be 2000*2000*2000, that will be a huge RAM consumption. I saw you mentioned about downscaling the projection images, any suggestion how to do that? Sorry I'm very new to tomography.

Much appreciate your kind help!

Sussi

If the object does not reach the edges of your projection images, then you could first of all crop them. However, you should carefully check that no part of the object is cut off in any of the images.

If you have placed the projections virtually in the center of the reconstruction volume, as I do in the tutorial, then you indeed need a reconstruction volume that is the size of the projection images. However, you could start by reconstructing only the central (2D) slice. This will surely fit in the GPU memory (if you also adapt the projections), and it'll give you an idea of how a full size reconstruction would look.

Downscaling the projection images is simply making them smaller, e.g., convert your 2000x2000 images to 500x500 images). Of course, you have to adapt the geometry accordingly.

I'll update this answer with more information on reconstructions that don't fit in the memory of the GPU in the future, but I'll have to ask for your patience for now.

Hi,

Thank you very much for sharing this example! I have successfully used it and I have a question about the use of Optomo object with the astra toolbox. I know this was not covered in your tutorial but maybe you are also familiar with it. After using astra.create_projector(..), I created an OpTomo object with astra.OpTomo(..) and I would like to get the dense projection matrix (on a toy example so that I won't run out of memory). Do you know if it's possible with this package and how to do it?

Thank you in advance,

Kind regards

I have very little experience with the OpTomo class, but I don't think that it's possible to get the projection matrix out of it. However, using the normal Python API you can use

`astra.projector.matrix()`

and`astra.matrix.get()`

to get the projection matrix. I've used that successfully (for small problems, of course), albeit using the MATLAB API of the Toolbox.Thank you very much for sharing your amazing tutorial. I have a question about loading datasets and pre-processing.

Currently, I'm working on a Cone beam CT reconstruction for a sample tube containing a certain animal's tissue. The outline of the tube has been reconstructed almost clearly, but the contents of the tube has not been reconstructed. Nothing is painted.

I think this is a pre-processing issue. When the projection images are loaded in float format, the maximum value is about 60000, but the minimum value differed from 10000 to 30000. I think some rescale and pre-processing are necessary.

Should the minimum value is 0 as in this case?, Or Should I apply pre-processing technique such as minus-log using tomopy ? Would you please advise me about something?

If your pixels have the values that you mention, then I would guess that your images are still raw X-ray images (with the dense parts of the sample dark and the background white). These cannot be used for tomography; they must be, as you mention, "log-transformed". I explain this on Tomography, Part 1: Projections, but maybe I need to add an extra tutorial, because this questions pops regularly… You can estimate the value of I(0) from the images by averaging background pixels.

Thank you for your message. After checking, as you say, there is a high probability that my images still raw X-ray images. I'll try to apply the log-transformed or flat-field correction(?). Thank you very much.

Yes, you need to apply both corrections. However, the "log-transform" is

muchmore important than the flat-field correction. I'd definitely start with that.## Add new comment