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.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)

Liebert (not verified)

Sat, 02/23/2019 - 11:21

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,

Tom

Mon, 02/25/2019 - 11:53

In reply to by Liebert (not verified)

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,

drue (not verified)

Thu, 11/28/2019 - 16:08

In reply to by XavierB (not verified)

thanks for the question, really helped

Tom

Mon, 02/25/2019 - 11:58

In reply to by Liebert (not verified)

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.

Drizzle (not verified)

Sat, 03/09/2019 - 09:56

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!

Sussi (not verified)

Tue, 06/25/2019 - 05:00

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.

Ann (not verified)

Fri, 10/11/2019 - 12:20

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 much more important than the flat-field correction. I'd definitely start with that.

Thank you for your reply, I will verify it.
Have you tried Sophiabeads Datasets? I'm currently challenging these public datasets, but I can't generate reconstructed images well. Following this tutorial, rescaling is done so that the reconstruction voxels get a size of 1x1x1. Is this setting available to all cases when using ASTRA-toolbox? Are there any geometries that cannot be applied?

adriano (not verified)

Thu, 11/28/2019 - 16:04

Hie Tom,
i have tiff images from a micro ct scan amounting to 1000 projections and they are about 7GB in memory size. i wanted to do a reconstruction of the images.

i seem to have difficulties in
1)storing projection data as a 3D array in the matlab environment
2)importing the projection and volume data to the astra environment in matlab

may you or anyone please assist with the neccesary steps as well as a code i can use to perform this

thanks

BIG (not verified)

Mon, 09/07/2020 - 05:10

Hi,tom,thanks for your expertise. my question is about , we have actually dataset matrix(x,y,z ),as you illustrated above'' (detector_rows, num_of_projections, detector_cols)'', then before reconstruction we need to define vol_geom(a,b,c), i want to know the relationship between (x,y,z) and (a,b,c)?

i mean if i set vol_geom according,''vol_geom = astra.creators.create_vol_geom(detector_cols, detector_cols,
detector_rows'', then run the algorithm, there is a error shows did not match the dimension . only when i set vol_geom as vol_geom = astra.creators.create_vol_geom(projections, detector_cols,detector_rows) the algorithm can run but the resolution seems not correct. my question is the any relations between (x,y,z) and (a,b,c);

hi,tom, the two questions above you can ignore them, may be i did not get the point. form my perspective , i want to ask how to compute sinogram from real projections not from simulation volume? thanks.

Hi BIG,

I'm working on the same thing.. Did you figure out how to do it? I'm working with Fluroroscopy Images, and I'm getting an extruded kind of shape instead of a real 3-d shape...

Is there a specific number of images that you need to get, for instance?

Thanks,
SRR

James (not verified)

Wed, 12/16/2020 - 10:39

Thanks for the tutorial. It was incredibly useful.

I am trying to just reconstruct the central slice from the full cone data set, while having the full data set in memory for future reconstruction.. I thought setting the vol_geom(detector_cols,detector_rows,1) would grab it, but doesn't seem to. Any advice?

parth paladiya (not verified)

Wed, 12/30/2020 - 15:53

i have used astra in my 3d construction program where i have 300 images from 1.2 degree angles and i am generating 3d model from that but my images are like gray color type so tomography avoid gray part which is only visible in 15-20 images out of 300 so is there any way i can cover that types of things in model?

Vicky (not verified)

Mon, 04/12/2021 - 16:43

Hi, Tom. It's a pretty nice and clear tutorial for the new developers like me to learn this toolbox.
I have a question about creating the projection geometry data. In the command "astra.data3d.create('-sino', proj_geom, projections)", you choose the object data set as "sino". Does this "sino" means sinogram? If it is means sinogram, should we need to use the command "creators_backprojeciton3d" to change it into sinogram. Then we do the reconstruction.

I ask this question is because I can use your example to test. However, when I use my own dataset, it cannot reconstruct the object which is only zero.

My object is a convex object which has a larger scale on the z-axis. I am wondering if it will generate some problems.

Could you give me some suggestions? Thanks so much!

The '-sino' is indeed short for sinogram. It tells the Toolbox that there is projection data in that particular data structure. The alternative is '-vol', short for volume, for the data structure that will hold the reconstruction.

The term sinogram is used here for a 3D object, although the name originally stems from the 2D sinogram that you get by projecting a single slice. The 3D volume of data is a stack of projection images in one direction and a stack of sinograms in another direction. So you don't need to use astra.creators.create_backprojection3d_gpu.

I'm currently working on a new tutorial that starts from a real dataset. I suspect that that will be useful for you when it is finished (part one is already published, but that's only the introduction...)

Max (not verified)

Thu, 06/24/2021 - 09:58

Hey, thank you for your Tutorial.

I tried it and the Code works perfectly. But with my own data, I got wrong reconstructions and I have no idea why...

I have 500 Projections of 200° then I have distance between the detector and the source of 500mm and between the specimen and the detector of 300mm. So i Thought doing this for my parameter:

distance_source_origin = 300 # [mm]
distance_origin_detector = int(500 - 300) # [mm]
Is that right?

detector_pixel_size = 0.2 # [mm]
detector_rows = 785 # Vertical size of detector [pixels].
detector_cols = 797 # Horizontal size of detector [pixels].
detector_sapcing_x = 159.4
detector_spacing_y = 157
num_of_projections = 500

angles = np.linspace(0, 2 * (200/360) * np.pi, num=num_of_projections, endpoint=False)

I think I have a geomtrical Problem but I can't find it by mysel at the moment.

The data looks like your example data with the Sinogram in X dimension and the Projections in Y dimension.

Thank you in advance

Juan (not verified)

Fri, 08/06/2021 - 16:44

Thanks for the tutorial. It was incredibly useful.

I am using Matlab instead of Python. Therefore I translated the main part of the example (parts I and II) to Matlab, joining the two parts to avoid writing the projections in disk, and also without Poisson noise to try to get a clean, almost perfect, reconstruction.

Here you can download my code: https://pastelink.net/36mdq

But I encounter a problem of artifacts in 3D reconstructions using CUDA with "cone" geometry.
These are the central slices of the reconstruction along the x-z and y-z plane, with unexpected artifacts:
https://ibb.co/wQNRpYd
https://ibb.co/kHRHMhs

I would like to know if these artifacts also appear in the Python implementation.
Could it be because the back-projectors and forward-projectors using GPUs are not matched?

It is a pity that you cannot compare using CPU projectors in 3D, which is not supported by ASTRA toolbox.

Thank you in advance!

These artifacts do not appear in the Python reconstruction, so there must still be something wrong with your Matlab code. Note that the illustration in this article is the central slice from the actual reconstruction made with the Python script.

Thank you for your prompt response!

I have installed AstraToolbox in Python and have found that the artifacts also appear with your code.

The slice to be checked is in the second or third dimension of the reconstructed volume. The first dimension (the one that is being saved is disk) is fine.

I am just adding this:
im = reconstruction[:,100,:] # Central slice, third dimension
import matplotlib.pyplot as plt
plt.imshow(im, cmap='gray')

If you look at the original reconstruction result, before removing values less than zero,
in the line (reconstruction[reconstruction < 0] = 0) the artifact is even more evident.

Without non-negativity constraint: https://ibb.co/vjwV8xX
With non-negativity constraint: https://ibb.co/cLscPN2

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 22 October 2018