Biomedical Physics & Engineering
Express
PAPER •
OPEN ACCESS
TIGRE: a MATLAB-GPU toolbox for CBCT image
reconstruction
To cite this article: Ander Biguri et al 2016 Biomed. Phys. Eng. Express 2 055010
View the article online for updates and enhancements.
Related content
GPU-based iterative CBCT reconstruction
using tight frame regularization
Xun Jia, Bin Dong, Yifei Lou et al.
-
Parameter selection in limited data cone-
beam CT reconstruction using edge-
preserving total variation algorithms
Manasavee Lohvithee, Ander Biguri and
Manuchehr Soleimani
-
A general method for motion
compensation in x-ray computed
tomography
Ander Biguri, Manjit Dosanjh, Steven
Hancock et al.
-
Recent citations
Time separation technique: Accurate
solution for 4D C-Arm-CT perfusion
imaging using a temporal decomposition
model
Sebastian Bannasch et al
-
Wei Zhang and Yan Kang-
Non-rigid CT/CBCT to CBCT registration
for online external beam radiotherapy
guidance
Cornel Zachiu et al
-
This content was downloaded from IP address 188.184.3.52 on 26/04/2018 at 09:55
Biomed. Phys. Eng. Express 2 (2016) 055010 doi:10.1088/2057-1976/2/5/055010
PAPER
TIGRE: a MATLAB-GPU toolbox for CBCT image reconstruction
Ander Biguri
1,3
, Manjit Dosanjh
2
, Steven Hancock
2
and Manuchehr Soleimani
1
1
Engineering Tomography Lab (ETL), Electronic and Electrical Engineering, University of Bath, Bath, UK
2
CERN, Geneva, Switzerland
3
Author to whom any correspondence should be addressed.
E-mail:
a.biguri@bath.ac.uk and m.soleimani@bath.ac.uk
Keywords: cone beam CT, image reconstruction, tomography software, GPU
Abstract
In this article the Tomographic Iterative GPU-based Reconstruction (TIGRE) Toolbox, a MATLAB/
CUDA toolbox for fast and accurate 3D x-ray image reconstruction, is presented. One of the key
features is the implementation of a wide variety of iterative algorithms as well as FDK, including a
range of algorithms in the SART family, the Krylov subspace family and a range of methods using total
variation regularization. Additionally, the toolbox has GPU-accelerated projection and back
projection using the latest techniques and it has a modular design that facilitates the implementation
of new algorithms. We present an overview of the structure and techniques used in the creation of the
toolbox, together with two usage examples. The TIGRE Toolbox is released under an open source
licence, encouraging people to contribute.
1. Introduction
Among the techniques for x-ray computed tomogra-
phy (CT) in widespread use, cone beam (CB) geometry
is getting increasing attention nowadays, from medical
imaging to material science. The possibility of recon-
structing full 3D images using a reduced x-ray radia-
tion dose is an important feature for CBCT
development in medicine and it has led to high-quality
3D reconstruction in micro-CT [
1]. Applications
include maxillofacial imaging [
2], guidance for radia-
tion therapy in oncology [
3], insect imaging [4] and
material science [5].
In all applications of CBCT, the working principle
is the same: 2D x-ray images of the ‘sample’ are
obtained from different angles and a tomographic
reconstruction algorithm is used to create an image
from the data. The fact that in circular CBCT the origi-
nal image is mathematically impossible to obtain [
6, 7]
and other factors, such as the high dimensionality of
the problem or the inconsistency created by different
physical effects with photons, make the image recon-
struction problem what mathematicians define as ill-
posed. Advanced mathematics is needed to generate a
solution. This has led to extended research in image
reconstruction algorithms, with a wide range of
published approaches that give differing results. And it
remains a hot topic.
While the use of CBCT is being increasingly exten-
ded to cover different imaging fields and research on
reconstruction algorithms still sees new methods pub-
lished, the end users of the images, both in medicine
and microtomography, mainly use the simplest recon-
struction algorithm, FDK [8]. This is worrying
because, while FDK produces satisfactory images for
good quality, full-projection, noiseless data, it per-
forms poorly in less ideal scenarios. It has been repeat-
edly demonstrated that iterative algorithms [
9]
outperform FDK [10–14].
There are a few factors that influence the lack of
connection between mathematics and usage. The
main one is computation time. Most, if not all, alter-
native algorithms are iterative. They need to recom-
pute repeatedly operations that are very memory- and
computationally expensive, while FDK needs less time
than a single such iteration. This is an important point
especially in medical applications, where an image is
needed rapidly as medical decisions are taken on the
basis of what can be read from it. The time scales for
iterative algorithms to run on a modern computer
CPU are of the order of hours, days or even weeks for
the most complex algorithms and bigger data. Another
OPEN ACCESS
RECEIVED
3 April 2016
REVISED
29 June 2016
ACCEPTED FOR PUBLICATION
8 July 2016
PUBLISHED
8 September 2016
Original content from this
work may be used under
the terms of the
Creative
Commons Attribution 3.0
licence
.
Any further distribution of
this work must maintain
attribution to the
author(s) and the title of
the work, journal citation
and DOI.
© 2016 IOP Publishing Ltd
factor is the lack of easy-to-use and free-distribution
iterative algorithms. While some of the most recent
toolboxes (presented later) do include some iterative
algorithms, the vast majority of these algorithms have
been completely ignored by both open source and
commercial image reconstruction software. This lack,
in conjunction with the fact that the research on
reconstruction algorithms requires a deep under-
standing of such fields of mathematics as linear algebra
and inverse problems, makes iterative algorithms
somewhat out of reach for the end users of the recon-
structed images.
In order to reduce the gap between algorithm
research and end use, we have developed the Tomo-
graphic Iterative GPU-based Reconstruction (TIGRE)
Toolbox, a MATLAB/GPU toolbox featuring a wide
range of iterative algorithms. It uses the higher level
abstraction of MATLAB with the lower (hardware-
specific) performance of CUDA in order to make it
fast and easy to use. In an attempt to bring the different
fields together, we addressed the computation-time
problem using the latest technologies in GPU comput-
ing with massive parallelization and memory manage-
ment efficiency. Only the main computationally
expensive blocks have been parallelized, with a mod-
ular design, allowing algorithm researchers easily to
plug new methods together with the GPU blocks pro-
vided. Additionally, the algorithms can be used as
single-line functions, giving total abstraction to
researchers who are only interested in the resultant
images, rather than in algorithm development.
Before explaining the specifics of the TIGRE tool-
box, it is worth mentioning some other toolboxes that
are also available. There are several commercial and
free software packages for FDK reconstruction,
including (but not exhaustively) CoBRA [
15] , Ultra-
fast CB CT reconstruction [
16] , OSCaR [17], Accel-
erating ConebeamCT [
18]. Additionally, some more
advanced toolboxes that include one or two iterative
reconstruction algorithms (SIRT and/or CGLS) are
also available, such as ASTRA [
19], RTK [20] and 3D
CB CT MATLAB [
21]. Of these, ASTRA and RTK are
the toolboxes that are most complete, however their
infrastructure in low-level programming languages
make them less suitable to work with when developing
new algorithms.
In this paper we briefly describe the CBCT image
reconstruction problem and some of the many ways to
solve it. Thereafter, we give an overview of the struc-
ture of the TIGRE Toolbox and show some perfor-
mance results and some reconstructed images. Finally,
we discuss the future vision of the toolbox.
2. Methods
2.1. CBCT geometry
The geometry of CBCT can be represented as in
figure
1. An x-ray source,
S
, is located at distance
DSO
from a centre of rotation
O
, where the origin of a
cartesian coordinate system is located. The x-ray
source irradiates a cone-shaped region containing the
image volume
and a detector
measures the
intensity of the photons impinging on it, photons that
have been attenuated following the Beer-Lamber law.
The image is centred at position
¢
O
, which is displaced
by
¾
V
orig
from the coordinate system origin. The
detector, located at distance
DSD
from the source and
centred at
¢D
, has an offset of
¾
V
det
from
D
, which is a
point lying in the xy-plane at distance
-DSD DSO
from the origin. A projection coordinate system uv is
defined centred at the lower left corner of the detector.
During the measurement acquisition, the source and
the detector rotate around the z-axis at an angle of α
from their initial position.
The geometric variables described above are used
in the TIGRE Toolbox to perform the necessary opera-
tions for image reconstruction, as shown in code snip-
pet
1. It is worth mentioning that both
¾
V
det
and
¾
V
orig
are vectors that define a single offset per projection.
Code Snippet 1. Geometry definition in TIGRE.
%% Geometry structure definition.
% Distances
geo.DSD = 1536; % Distance Source Detector
geo.DSO = 1000; % Distance Source Origin
% Detector parameters
geo.nDetector=[512; 512]; % number of pixels
geo.dDetector=[0.8; 0.8]; % size in mm of each pixel
geo.sDetector=geo.nDetector.
*
geo.dDetector; % total size of the detector in mm
% Image parameters
geo.nVoxel=[512;512;512]; % number of voxels in the image
geo.sVoxel=[256;256;256]; % total size of the image in mm
geo.dVoxel=geo.sVoxel./geo.nVoxel; % size in mm of each voxel
% Offsets
geo.offOrigin=[0; 0; 0]; % V_orig
geo.offDetector=[0; 0]; % V_det
2
Biomed. Phys. Eng. Express 2 (2016) 055010 A Biguri et al
2.2. Image reconstruction problem
For a given geometry, image reconstruction can be
described by two different approaches. The first one is
by solving the inverse Radon transform, a mathema-
tical tool to describe the integral of a function over
straight lines. The solution of the inverse Radon
transform is well known in tomography and Feld-
kamp, Davis, and Kress modified it for CB geometries.
Their solution is known as the FDK algorithm [
22].
While the FDK algorithm will produce good quality
images in an ideal scenario, it copes poorly with
common unaccounted sources of error, such as beam
hardening or photon scattering [
23].
Alternatively, the image reconstruction problem
has been described as a minimization one as in
equation (
1), where b are the projection data, x is the
image and A is a matrix describing the intersections of
x-rays and voxels in the image. In this equation, G(x) is
an optional term that describes a regularization func-
tional
(·
)G
. This functional can be used to introduce
additional constraints to the image reconstruction
algorithm
=-+
ˆ
() ()xbAxGxargmin . 1
x
2
While the minimization formulation allows the
use of advanced linear algebra techniques, there is a
significant complication: the size of the matrix A. This
is especially important as most, if not all, iterative tech-
niques to solve equation (
1) use one Ax and one A
T
b
matrix-vector multiplication. As an example, matrix A
for the geometry described in code snippet
1 for 360
projections has 94 371 840×134 217 728 elements
with a sparsity index of 0.0017%. This requires 320 Gb
of memory even using optimized sparse memory
methods.
In order to cope with a problem of this scale, the
most common approach is to substitute the matrix-
vector multiplications Ax and A
T
b by operators A(x)
and
()
A
b
T
, recomputing the relevant matrix values
whenever necessary. While computationally very
expensive to perform, the operators have an
advantage: the values are completely independent of
each other, making them suitable for parallel comput-
ing. In the TIGRE Toolbox, these two operators have
been implemented using CUDA in a GPU capable of
computing over 60K floating-point operations simul-
taneously
4
. This has resulted in speed-ups of up to
1400 times compared to matrix-based methods.
2.3. Toolbox structure
In this section an overview of the structure of the
toolbox is given (see figure
2). As mentioned in the
previous section, the main building blocks of any
iterative algorithm are the so-called projection (A(x))
and back projection (
()
A
b
T
) operators. In the TIGRE
Toolbox, these two blocks have been optimized for
GPU computing using CUDA. They lie in the lowest
layer of the toolbox design and are constantly used by
the other layers. The algorithms themselves lie in the
topmost layer and are all coded in MATLAB, which
provides the power and flexibility of a high-level
language. To be able to communicate between the
low-level, hardware-oriented CUDA and the high-
level, design-oriented MATLAB, a set of the so-called
MEX functions are needed. The toolbox has been
designed not to have any specific data types or classes.
Instead, it comprises only the basic MATLAB types,
such as matrices and structures.
2.3.1. Projection and back projection
As already mentioned, the main building blocks of the
toolbox are the CUDA/C++ implementations of the
projection and back projection operators. Concep-
tually, the matrix A is a linearization of the model that
describes the x-ray attenuation measured over a given
domain and several different approaches to compute
this may be found in the literature. Similarly, back
projection is a ‘smearing’ over the domain with a
weighting applied. Without explaining in detail all the
Figure 1. Geometry of CBCT.
4
In specific GPU models.
3
Biomed. Phys. Eng. Express 2 (2016) 055010 A Biguri et al
methods available, we briefly describe those used in
the toolbox.
For the projection, two approaches have been
implemented: the voxel-ray intersection approach and
the interpolation approach. The first of these uses the
Siddon ray-tracing algorithm [
24] with optimized
operations [
25]. This algorithm computes the distance
between a given voxel and an infinitesimally narrow
x-ray beam and multiplies that by the voxel intensity.
This approach is the fastest way of computing the pro-
jector. However, it is known to introduce discretiza-
tion square-block artefacts due to the finite size of the
voxels, artefacts which become more significant the
bigger the voxel size. To avoid this problem, a trilinear
interpolation approach has been implemented where
the path integral is evaluated every fixed
D
l
and image
values are interpolated using advanced texture mem-
ory. To implement this an additional variable is added
to the geometric definition of the problem:
=geo.accuracy 0.5
, which defines
D
l
as a fraction of
the voxel size. This fraction is best chosen to be 0.5 or
lower, as Jia et al [
26] demonstrated.
For the back projection, two different approaches
based on the same concept are used. Initially a ray is
linked from the source location to the desired voxel,
and extended to the detector. There, using bilinear
interpolation, a value is read and added to that voxel
with a weight. The difference between the two back
projections is in this weight. One of them implements
the FDK weight. However, this makes the back projec-
tion unmatched, i.e., it makes the back projection
operator not equivalent to the transpose of the pro-
jector. While not important for most algorithms, this
is crucial for Krylov subspace methods. In order to
change that, a matched weight as described by Jia et al
[
27] is used. While not completely matched, they claim
that it is above 99% similar to the transpose of matrix
A. Both back projectors perform similarly.
2.4. Algorithms
One of the key features that we wish to introduce with
the TIGRE Toolbox is algorithm variety. The field of
image reconstruction has seen the development of a
wide variety of methods to solve equation (
1) using
different solvers and regularization techniques. There
are four main families of reconstruction algorithms
present in the current implementation of the toolbox:
the filtered back projection family, the SART-type
family, the Krylov subspace method family and the
total variation regularization family. A brief descrip-
tion of each algorithm subgroup follows, together with
which algorithms are included in the toolbox.
The filtered back projection family is a set of algo-
rithms based on solving the inverse of the Radon
transform. Different variations of the algorithm have
been proposed in the literature, but the toolbox con-
tains just the standard FDK implementation with a
small choice of filtering kernels
5
.
The SART-type family [
28] is set of algorithms that
derives originally from the Kaczmarz method and is
adapted to work projection by projection instead of
row by row. This family of algorithms follows
equation
l=+ -
+
() ()xxVAWbAx,2
kkkT k1
where V and W are weight matrices based on ray
length. The algorithms of this family mainly differ by
the number of projections used simultaneously. In the
TIGRE Toolbox, SIRT, OS-SART [
29] and SART are
implemented, where the image is updated using all
projections, subsets of projections or projection by
projection, respectively. Additionally, the toolbox
provides different options for tuning the algorithms.
For example different initialization techniques are
implemented, such as FDK, multi-grid, or user-
specified image. The main difference between the
performance of the algorithms in this family is in
convergence versus speed. The more data used in one
update, the faster the algorithm will be per iteration,
but slower (in number of iterations) to converge. For a
more accurate solution, SART is suggested, while a
faster result is obtained using SIRT, and with OS-
SART somewhere in between.
Krylov subspace methods consitute a set of faster
algorithms for solving linear equations. They iterate
through Krylov subspaces, minimizing the eigenvec-
tors of the residuals in descending order and so have
Figure 2. Structure of the TIGRE Toolbox.
5
FDK adapted from 3D CB CT MATLAB [21], with permission.
4
Biomed. Phys. Eng. Express 2 (2016) 055010 A Biguri et al