Tomographic Reconstruction Methods for
Decomposing Directional Components
Rasmus Dalgas Kongskov∗

Yiqiu Dong∗

arXiv:1708.06912v1 [cs.CE] 23 Aug 2017

January 25, 2018

Abstract
Decomposition of tomographic reconstructions has many different practical application. We propose two new reconstruction methods that combines the task of tomographic reconstruction with object decomposition.
We demonstrate these reconstruction methods in the context of decomposing directional objects into various directional components. Furthermore we
propose a method for estimating the main direction in a directional object,
directly from the measured computed tomography data. We demonstrate all
the proposed methods on simulated and real samples to show their practical
applicability. The numerical tests show that decomposition and reconstruction can combined to achieve a highly useful fibre-crack decomposition.

1

Introduction

X-ray Computed Tomography (CT) is a highly used non-invasive imaging technique. Applications of this technique ranges from biological and chemical science,
to structural and material science, where the resolution also varies from large scale
(meters) to micro-scale (nano-meters). The technique is based on the attenuation
that X-rays undergoes as they pass trough matter. The attenuated X-ray intensity is measured at a detector. In this paper our CT set-up is 2D and parallel
beam, but the methods introduced here are not limited to this CT-methodology,
but presented in this setting for simplicity. The X-ray source and detector is
rotated around the object, and for each angle of rotation the attenuated X-ray
signal is measured.
Lambert-Beers law relates the object that we are scanning to the measured
intensity. By rewriting Lambert-Beers law, the measured intensity is seen to be
linearly related to the object through line-integrals [4]. If we model the object
as a continuous function the relation between the measured intensity and object
corresponds to the Radon transform. In this paper we focus on the discretized
version of this model where the domain with the object, the reconstruction domain, is split into M ×M equidistant pixels. The discrete model can be expressed
as a simple matrix-vector product as
Az = b.

(1)

∗ Department of Applied Mathematics and Computer Science, Technical University of Denmark, 2800 Kgs. Lyngby, Denmark (rara@dtu.dk, yido@dtu.dk.

1

2

Here the system matrix, A ∈ RN ×M , models the length of each X-ray through
2
each pixel in the reconstruction domain. The reconstruction, z ∈ RM , is the
object, including it’s interior, that we seek to find. The sinogram data, b ∈ RN ,
is the measured data at the detector. The data is measured at a series of discrete
angles, φ ∈ RNφ , where Nφ is the number of discrete angles. For simplicity we
assume that we have parallel beam measurements on a 180◦ degree arch, i.e.
not limited angle measurements. We further assume that we do not have sparse
angle measurements. Here sparse angle measurements refers to the case where
the amount of measurement angles is significantly smaller than the amount of
detector bins Nt . In order make the dimensions fit we have N = Nφ Nt . The
application of A is usually called a forward projection, whereas the application
of AT is a called a backward projection.
The process of constructing z based on known b is usually called a reconstruction. The most widely used reconstruction method for CT-problems was
introduced in [26] and is now known as Filtered Back-Projection (FBP). The
FBP method implicitly assumes to have continuously measured data from the
whole 180◦ angular range. For a discrete problem this assumption is violated,
but with sufficiently many equidistant angular measurements, and a low no noise
level, this method will give results close to the exact object. The main advantage
of FBP is the efficiency since the computational effort of this algorithm relies on
a two Fourier transforms, which numerically can be handled by use of the fast
Fourier transform (fft). The disadvantage of FBP is that the method is sensitive
toward measurement noise and any type of limited data problems.
Variational methods have been used in many types of inverse problems, see
[1, 31] for more on this topic. The main reason for using variational methods,
is that can help us to overcome unwanted the effects of ill-posedness, noise and
limited data. Another advantage of using variational methods is the possibility
for including prior knowledge through regularization. In this paper we will make
use of variational methods for solving the reconstruction problem, which will be
on the form
1
(2)
min kAz − bk22 + R(z).
z
2
The first term, the data-fidelity term, is chosen to be the 2-norm since we model
noise as being Gaussian. The second term is the regularization term and here
we will use several different methods. Typically a regularization parameter, that
balances the two terms, is also present in a variational formulation, here we
include this parameter as part of the regularization term. One of the most widely
used regularization methods for imaging problems is total variation (TV), first
introduced in [29] and used for tomography problems in e.g. [7, 14, 32].
After reconstructing an object one is typically interested in analyzing the
result, which could be aided by the use of segmentation. If we expect that our
object consist of several distinct features that are mixed together, one could
consider to decompose the image, i.e. split the object into different parts. In
order to decompose an object a decomposition model that describes how the
object should be split is needed. In [5] a way to decompose a signal using infimal
convolution is introduced. For two convex functionals J1 and J2 a function is
linearly split into two components as z = z1 + z2 and the infimal convolution of
the two functionals becomes
J(z) = inf {J1 (z1 ) + J2 (z2 )} .
2

This strategy relies on the assumption that each component, z1 and z2 , attains
a smaller functional value for the functional which models the given component.
Meyer further investigated image decomposition in [22], where methods for modelling and extraction of oscillating patterns in images was investigated in great
depth. Theese earlier works sparked research of image decomposition methods in
various aspects, see e.g. [8, 33, 2]. Variational methods for inverse problems is
very closely related to the decomposition methods and in [3, 10] denoising was
combined with image decomposition to texture and structure parts. The idea of
using image decomposition for inverse problems is to split an object linearly into
two components, z = u + w, and let the regularization consist of two different
terms, one for each component. In the more recent work [11] infimal convolution
was shown to give good results when combining several Total Variation (TV) type
functionals.
Solving (2) for z = u + w and different regularization terms for u and w,
will lead to a solution where the two components will follow different priors, and
hence a reconstruction and a decomposition is achieved simultaneously. Other attempts of decomposing CT reconstructions has been through segmentation, either
as post-processing, see e.g. [28], or simultaneous with the reconstruction, see e.g.
[27]. Decomposition of an object by use of spectral CT has also been attempted,
see [20]. In [19] a method for simultaneously decomposing and reconstructing an
object is introduced. In this paper the object is decomposed into a three components: the ground truth object, the limited data artifacts and the measurement
noise, which are assumed to be sparse in the wavelet basis, the discrete cosine
transform basis and in the sinogram domain, respectively.
Directional objects are objects where the texture in the object follows one
main direction. Applications with directional objects are, among other things,
fibres, such as optical fibres, glass fibres, carbon fibres, etc. When analyzing fibre
materials, CT scanners can be used to investigate interior properties, for example
irregularities, see [28, 30, 12]. A specific irregularity which is often sought for in
fibre materials are cracks which are distinctive from the fibres. Both the fibres and
the cracks can be regarded as directional components, that are part of an object.
Recently a regularization method for directional objects has shown promising
results in [16]. The first order regularization method presented in this paper is
based on TV and hence it is called Directional Total Variation (DTV).
In this paper we propose two combined reconstruction and decomposition
methods that benifits from using a variational formulation. The one reconstruction method is based on the infimal convolution decomposition method, where
solving the variational minimization problem decomposes the object into two, or
more, components. The other reconstruction method is motivated from the microlocal analysis results in [25], and utilizes regularization through a variational
formulation to overcome consequential limited data artifacts. We introduce both
reconstruction methods in the light of reconstructing directional objects, and we
demonstrate advantages and disadvantages of the reconstruction methods through
empirical examples. In order to use DTV regularization for the CT reconstruction
problem we need to supply the main object direction. We therefore propose a
method for estimating the main objection direction directly from sinogram data.
A practical application of the directional decomposition methods could be a fibre
material with cracks that occur across the main object direction, i.e. perpendicular to the direction of the fibres.
The paper is organized as follows. In section 2 directional regularization for
3

CT problem is introduced and demonstrated. Furthermore we give a method for
estimating the main direction of an object directly from the measured sinogram
data. In section 3 we introduce the sinogram splitting method, where the entire reconstruction problem is split into two individual problems. In section 4
DTV-decomposition, a variational decomposition model based on two different
DTV regularization term, is introduced. Both sinogram splitting and DTVdecomposition makes use of the direction estimation method and DTV, but one
method is obviously advantageous to the other and we therefore introduce both
method. Numerical experiments are carried out in section 5, where the two introduced decomposition methods are tested and compared. In section 6 conclusions
are drawn.

2

Directional regularization for CT problems

We introduce a regularization method for CT problems where the object of interest is directional, i.e. the texture of the object follows one specific main direction.
Furthermore we propose a method for estimating the main direction of an object
directly from the sinogram data.
The domain consist of M -by-M pixels with equidistant pixel-grid-size 1 × 1.
In this region (i, j) denotes a pixel index with 1 ≤ i, j ≤ M such that ui,j gives
the pixel value at (i, j). Since the scanning is circular, the region of interest, and
our simulated phantoms, are contained in a circle with the diameter M . The
sinogram domain is discretized in Nt -by-Nφ , so we have b ∈ RNt ×Nφ . In the
sinogram domain (l, m) denotes the pixel index such that bl,m is sinogram value
at detector bin tl and measurement angle φm .

2.1

Directional regularization

Regularization methods has been used for CT problems for almost 40 years, see
[18, 21] for some of the first works in this direction. Over time numerous regularization techniques has been adapted, developed and utilized for the objectreconstruction process. Since the publication of the paper that introduced the
ROF-model for imaging problems [29], Total Variation (TV) has been used for
a wide range of imaging applications. In discrete terms the total variation of
z ∈ RM ×M can be expressed as
X
TV(z) =
|(∇z)i,j |2 ,
i,j

for euclidian norm | · |2 and the discrete gradient operator ∇ : RM ×M → R2M ×M .
∇ is defined as


∇x 1 u
∇z =
,
∇x 2 u
where ∇x1 and ∇x2 , for the two dimensions x1 and x2 , are obtained by applying
a forward finite difference scheme with symmetric boundary condition, i.e.,


ui+1,j − ui,j , if i < M,
ui,j+1 − ui,j , if j < M,
+
+
(∇x1 z)i,j =
and (∇x2 z)i,j =
0,
if i = M,
0,
if j = M.

4

A regularization method for processing images while incoporating directional
priors has been introduced for image denoising and image deblurring problems in
[16, 15]. The first order regularization method, DTV, builds on the prior that the
target object is piece-wise constant with the texture along one main direction. In
this paper we will use DTV for regularizing the CT-reconstruction problem. In
discrete terms we can express DTV as
X
DTVθ,a (z) =
|Rθ Λa (∇z)i,j |2 ,
i,j

for scaling and rotation matrices


cos θ − sin θ
Rθ =
cos θ
sin θ


and

Λa =

1
0


0
,
a

(3)

where a ∈ (0, 1] and θ ∈ (0, 2π].
Fibre materials are objects that are often examined using CT, and these fit
well with the prior for DTV, since they are piecewise constant have highly unidirectional texture. We therefore center our examples around fibre materials or
simulated fibre materials. We give an example where the aim is to reconstruct
a simulated directional phantom. In this example we are dealing with a slightly
underdetermined system, 23 , with 1% Gaussian sinogram-noise. We compare the
results produced by FBP, `2 -TV and `2 -DTV by visualizing the results and by
means of the peak signal-to-noise ratio (psnr) measure in fig. 1. By visual and
quantitative comparison we see the advantage of incorporating the additional
prior for fibre objects.

Figure 1: Simulated CT-problem solved using three different reconstruction techniques in order to compare their performance. Regularization parameters for the
`2 -TV and `2 -DTV methods are tuned to maximize peak-signal-to-noise (psnr)
ratio. DTV parameters are chosen as a = 0.15 and θ = π9 (20◦ ).
In order to regularize an inverse problem using DTV, the main direction should
be know. Next we introduce a method for estimating the main direction directly
from the measured CT-data. Furthermore a width-parameter a should be chosen,
but based on the conclusions from empirical tests in [16] we just set a = 0.15 in
this paper.

2.2

Direction estimation from CT-data

In order to use DTV for regularization we have to select the parameters a and θ.
From [16] we know that if the directional prior is met we should choose a ∈]0, 1]
5

as small. Empirical tests in [16] show that a = 0.15 is a reasonable choice for
images, we choose this value for the CT reconstruction problem.
In [15] a method for finding the main direction in corrupted images, using the
Fourier transform, was presented. This method relies on Fourier transforming the
corrupted image and analyzing the magnitude of the coefficients. If the image
texture has a main direction, the coefficients corresponding to basis-functions
with same main direction, will be relatively large in magnitude.
To find the main-direction θ in the object of interest we can utilize that CT
data already consist of angular measurements. Given that we have data measured
from one or more angles relatively close to the main direction of the object, we
can detect this main angle by comparing the individual angular measurements.
Inspired by the method in [15] our method consist of a Fourier transform, though
in 1D, for each angle, i.e. along the detector-pixels. We sum the magnitude of
the Fourier coefficients, for each angle, and find the maximum response. The
algorithm is outlined in algorithm 1. Data measured at the same angle as the
main direction will be oscillatory. An oscillating signal is well represented by the
Fourier basis-functions, and therefore we expect the magnitude of the Fourier
coefficients to be largest at one specific angle, which will therefore become the
estimated main direction.
Algorithm 1 Main Direction Estimation Method
1:

Input siogram data b and discrete measurement angles φ.
1D Fourier transform along each angle φm :
2πιωl
PNt −1
b̂ω,m = l=0
bl,m e− Nt .

2:

Calculate summed magnitude of b̂ for each angle φm and find max response:
P
Wm =
|b̂ω,m |,
ω

h = argmax Wm .
m

3:

θ = φh .
return θ.

It should be noted that this angle estimation method is not limited to parallelbeam tomography since a pattern also occurs when the fan-beam geometry is
used, although the data should be re-binned or the method should be modified
slightly to account for this.
We have tested the angle estimation method on a simulated phantom and
on a real-data object, where we in both cases simulated the projections and the
noise. The real data object is obtained from [12], for more this data see [13].
From the noisy sinogram data we have estimated the main direction of the object
using algorithm 1. We tested the method on the two samples with additive
Gaussian noise of varying noise-level (ρ) and we report the angle-estimates in
table 1. The two tested samples together with their corresponding noise-free
sinograms are shown in fig. 2 and fig. 3. In these two figures we also indicated the
angle estimated from the noise-free sinogram and the angle-dependent summed
magnitudes W from algorithm 1.
Based on these empirical tests we see that our method is robust up until a
noise-level of at least 20%, which is a high noise-level for CT-data.

6

ρ (%)
Phantom
Real

0
20.1
81.5

1
20.1
81.7

3
20.1
81.5

5
20.1
80.9

10
20.1
81.7

20
20.1
79.5

30
20.1
-1.1

40
31.7
-34.9

Table 1: Angle estimation results. See phantoms in fig. 2. For the simulated
phantom the exact angle is 20◦ and for the real data sample we do not have an
exact value.

Figure 2: To the left: Fibre-phantom with estimated main direction marked by
dashed line. To the right: Noise-free sinogram overlayed with a line-plot of W
(see algorithm 1).

Figure 3: To the left: Fibre reconstruction (obtained from [12], see more in [13]),
with estimated main direction marked by dashed line. To the right: Noise-free
sinogram overlayed with a line-plot of W (φ) (see algorithm 1).

7

3

Sinogram Splitting

Edges are essential in the analysis of CT-scans since they indicate material boundaries. An edge is a singularity. A singularity is a point where a function, or its
derivative, is discontinuous. From microlocal analysis we know that a singularity
is only visible in the sinogram, if an X-ray is tangent to the singularity in the
object, see [25].
In directional objects the texture follows the main direction and therefore the
texture edges also follows this direction. The edges along the main direction will
be tangent to the X-rays for a specific angle, or within a limited angular range in
practical applications. Therefore measurements within this limited range of angles
will contain information about the main direction singularities. Based on this
observation, we propose a method where the sinogram-data is split into two parts,
according to the measurement angles, such that each part will include information
about singularities along different directions. For a fibre-crack object, this method
could be used to split the data into two parts: one part with singularities of the
fibres and one for all other singularities, i.e. cracks etc.

3.1

Splitting model

In [25], singularity propagation, from an object through the Radon transform, is
examined using wavefront sets and microlocal analysis. Wavefront sets denotes
points where a given function is not smooth, together with the direction in which
the function is not smooth. In [25] the relation between singularities in a function
f and singularities in the corresponding Radon transform Af is described. The
paradigm that is described in [25] is further outlined in [17] as follows:
A detects singularities of f perpendicular to the line of integration (”visible”
singularities), but singularities of f in other (”invisible”) directions do not
create singularities of Af near the line of integration.
It should be noted here that a singularity perpendicular to the line of integration
corresponds to an edge tangent to the line of integration.
The fact that singularities only propagate when rays are perpendicular to them
inspired us to decompose the linear CT-model in separate parts, where each part
is related to an object-component where the directions of the singularities, in the
object-domain, are limited. The simplest split is a two-component split, which
we will base the following on, but the method can be generalized to any integer
amount of components with different directions. Components are separated according to the measurement angles: the first part of the model is based on K + 1
consecutive measurement angles, φu ∈ RK+1 , whereas the second part is based on
the remaining angles, φw ∈ RN −K−1 . This splitting results in two linear systems,
Au u = bu
2

and

2

Aw w = bw ,

(4)
2

where u ∈ RM and w ∈ RM are the object components. Au ∈ R(K+1)Nt ×M
2
and Aw ∈ R(N −K−1)×M are the split-matrices. And bu ∈ R(K+1)Nt and bw ∈
R(N −K−1) are the split measurements. This splitting method is not so suitable for
limited data problems, such as limited angle and sparse angle problems, since for
such problems we might not have measurements of the singularities that we need
in order to be determine φh . Further more K > 0 is required to have a reasonable
8

splitting, otherwise the data is simply to limited for the one sub-problem and the
reconstruction will be unreliable.
For the example of a split between fibre and crack components, we can split the
data and the matrices according to the main object direction θ. For main direction
θ = φh and even range-width-index K, we can pick φu = (φh− K , . . . , φh , . . . , φh+ K ).
2
2
This splitting has the downside that objects that do not follow the main
direction are very likely to end up in the crack-component, which could complicate
a complete fibre-crack segmentation.
Next we introduce two different ways of solving the linear systems in (4). If
we use the same linear reconstruction method on both sub-problems we will have
have z = u + w and we can therefore analyze both the combined reconstruction
and each component individually. On the other hand we can choose a non-linear
reconstruction method and use the flexibility of this to impose prior information
on each component.

3.2

Reconstruction methods

We introduce two different ways to solve the two linear systems (4). The first
solution strategy is to use FBP on each system. For the second solution strategy
we focus on the variational methods, due to their applicability for inverse problems. This will not necessarily give us a linear reconstruction method, but on
the other hand, it will allow us to easily impose some prior information through
regularization.
When using standard FBP for limited angle problems, the data from the missing measurement angles is implicitly assumed to be 0, since line integration over
0 will give no impact on the reconstruction. For the two-component split (4) this
corresponds to the assumptions: Au w = 0 and Aw u = 0. This will result in
artifacts which comes from the artificial singularities in the sinogram at the transition between the measured data and the assumed 0-data. For more on limited
angle artifacts see [23, 9]. Since FBP is linear, the limited angle artifacts will
cancel when summing the two components u + w. One can therefore examine
the separate components with limited angle artifacts or summed reconstruction
without the limited angle artifacts. It should be noted, as also mentioned previously, that the FBP method is not suited for handling underdetermined systems
or noise-corrupted data.
Variational methods has been shown to be advantageous for CT-reconstruction
problems with limited and/or corrupted data. The motivation for this approach
is to solve the sub-problems using prior information about the individual components. So for each component we pick a suitable regularization method and solve
the corresponding sub-problem. Formulating the two problems (4) as variational
problems and incorporating regularization we get the two optimization problems
min kAu u − bu k2 + Ru (u)
u≥0

and

min kAw w − bw k2 + Rw (w).

w≥0

Both of problems are still limited angle problems, but now the assumptions are
only imposed through the regularizers Ru and Rw . Since each component is
solved individually, it is possible to incorporate the most suitable regularization
technique for each of them.

9

For the fibre-crack decomposition, the prior for the fibre-component is that it is
piecewise constant and that the texture follows one main direction. The prior for
the crack-component is that it is independent of the direction, piece-wise constant
and sparse. Based on these priors we suggest to regularize the fibre-component
using DTV and the crack-component using classical TV together with a `1 -norm.
In discrete terms we have
Ru (u) = λu DTV(u)

and Rw (w) = λw TV(w) + βkwk1 .

Both of these regularization choices are convex and due to the piece-wise constant
prior they will diminish the undesirable effects of noise.
Solving the sinogram splitting problem using the above suggested regularizers
will leave us in total with four tunable parameters: The range-width-index K
should be chosen to match the expected range of the fiber direction, i.e. small K
for a highly one-directional object and larger K for an object with fibres within
a range of angular directions. The regularization parameters λu > 0 and λw > 0
which control the balance between the fidelity terms and the regularizers. And
w-sparsity parameter β > 0 which should be relatively small, compared to λw , in
order to avoid w = 0.
A general advantage of using FBP instead of a variational method is the
efficiency, since FBP only requires one single back-projection opposed to several hundred, or thousand, forward- and back-projections when the optimization
problem is solved iteratively. General advantages for using a variational method
instead of FBP is that artifacts and noise effects can be diminished or removed,
while desired features, such as piecewise constancy, can be enhanced.

4

DTV-decomposition model

We wish to impose different priors through regularization on each component in
a decomposition model, and hence split the object. This strategy relies on the
assumption that the desired components attains a smaller functional value for a
desired regularizer opposed to an undesired one.
For an object which consist of some regions that fulfill the directional prior
and other regions that do not fulfill the prior, we want to develop a model where
the two regions are reconstructed separately as different components. Inspired by
previous works within texture-cartoon decomposition, see [3, 10], we opt to build
a variational model where the minimizer will be a decomposed reconstruction.
In many applications fibre-structures are analyzed with the aim to detect
cracks and/or other types of deterioration. Whereas the texture of the fibre material follows one main direction, the deteriorated parts are mainly perpendicular,
or close to perpendicular, to the main direction. Moreover the deteriorated parts
are much less dominating in the object, opposed to the fibre materials. Based on
these observation we propose the following decomposition model:


1
kA(u + w) − bk22 + λ DTVθ,au (u) + αDTVθ⊥ ,aw (w) + βkwk`1 . (5)
min
u≥0,w 2
2

2

In this model u ∈ RM represents the fibres and w ∈ RM the crack part. The
summation of these two components is assumed to model the object, hence the
summation in the first term of (5). u is assumed to follow the main objectdirection θ, to be piece-wise constant and to be non-negative since it corresponds
10

to attenuation coefficient values. w is assumed to follow the orthogonal main
direction θ⊥ , to be piece-wise constant and to be sparse.
The parameters of the decomposition model (5) are the following: The regularization parameter λ > 0, which controls the trade off between the regularizer
and the fidelity term. The main direction angle θ which can be estimated with a
fast and noise-robust algorithm, see section 2.2. The first DTV-width au which
we assume to be fixed due to previous empirical tests, see [16]. The balance of
DTV-terms α > 0 which should be chosen in combination with the second DTVwidth aw ∈ (0, 1] to achieve sufficient splitting between the two components. The
orthogonal main direction θ⊥ which is just orthogonal to the main direction. And
the crack-sparsity parameter β > 0, which should be chosen as a relatively small
value, compared to λ, to avoid forcing w to be zero. This leaves four tunable
parameter: λ, α, aw and β. In order to avoid an overlap of the feasible sets of
the two DTV terms we have the following restriction for α:
au < α <

1
.
aw

This bound is obtained by comparing the two DTV functionals. We want to avoid
both
|Rθ Λau y|2 < α|Rθ Λaw y|2

and α|Rθ Λaw y|2 < |Rθ Λau y|2 ,

(6)

for all y ∈ R2 . Using the matrix definitions in (3) and the fact that θ⊥ = θ +
we can rewrite the two norm expression above and achieve:
p
(y1 cos θ + y2 sin θ)2 + a2u (−y1 sin θ + y2 cos θ)2
|Rθ Λau y|2 =
p
α|Rθ Λaw y|2 = (αaw )2 (y1 cos θ + y2 sin θ)2 + α2 (−y1 sin θ + y2 cos θ)2

π
2

where y = (y1 , y2 ). We have au , aw ∈]0, 1[, since we do not want DTV to become
classical TV, i.e. a = 1. Based on the two norm-inequalities in (6) we get
1 < αaw < α

∨

au < α

& αaw < 1

∨

α < au < 1.

Two of the above inequalities contradict each other and the remaining two gives
the bound for α
au < α

& αaw < 1 ⇐⇒ au < α <

1
.
aw

The model (5) is convex, which is desirable when we want to find a solution
to the minimization problem. Furthermore the sparsity constraint is not only a
reasonable regularization method for w, it also makes (5) strictly convex, i.e. the
minimizer will be unique.
Having introduced two different methods for combined decomposition and reconstruction, sinogram splitting and DTV-decomposition, we sum up the advantages and disadvantages of the two methods: Both methods include four tunable
parameters, so in this regard they are similar. The sinogram splitting method
risks, by design, to reconstruct incorrect attenuation coefficient values. If we
choose the DTV reconstruction method for the sinogram splitting, the components are not summable. On the other hand, if we choose the FBP reconstruction
11

method the components are summable, just as for DTV-decomposition where the
sum of the components gives the object. The DTV-decomposition method assumes that the crack-component follows a specific direction, where the sinogram
splitting method includes no such assumption, only an assumption on the fibrecomponent.
For the special case where Au w = 0 and Aw u = 0 we have that the datafidelity terms for the two methods are the same since
 
 
Au w = 0
Au
bu
⇐⇒
(u + w) =
⇐⇒ A(u + w) = b ⇐⇒ Az = b.
Aw
bw
Aw u = 0
The target for using regularization is also quite different for the two methods.
For both methods the regularization is used to suppress noise. The sinogram
splitting method furthermore uses the regularization to overcome the limited angle
artifact that occur due to the splitting. The DTV-decomposition method utilizes
the regularization directly for decomposing the components.
Based on the advantages and disadvantages above, it is not obvious that one
method should be superior to the other in all cases, which is why we empirically
examine both method in the next section.

5

Numerical Experiments

In this section we demonstrate the performance of the methods introduced in
section 3 and 4. We do this for a range of simulated X-ray CT problems and
we compare the performance of the proposed methods. In order to set the stage
for the numerical experiments we first give some discretization and experiment
details which are valid for the following tests.
We solve the variational optimization problems using the Primal-Dual-HybridGradient method from [6]. We stop the algorithm when the relative change of the
objective function goes below a threshold tollerance. In all the experiments this
tollerance is set to 10−5 . All of the algorithms are implemented in Matlab, where
we use the parallel beam GPU code described in [24] from the ASTRA toolbox,
see [35, 34], to calculate forward and backward projections. In all the numerical
tests we set the tolerance to 10−5 and we pick the optimal regularization by
ground truth comparison based on the peak signal-to-noise ratio (psnr), unless
otherwise stated. We use a simulated crack phantom, which has cracks in a 360◦
circular pattern, to demonstrate the performance of the previously introduced
reconstruction methods. The crack-phantom ground truth can be seen in fig. 4.
This phantom does not necessarily fit our assumptions, but serves as a stress-test
of the introduced methods.

5.1

Sinogram splitting

For the sinogram splitting method we compare the two reconstruction techniques
presented in section 3.2, namely FBP and the variational method approach. Both
reconstruction methods are tested on the same data-set, which is corrupted with
1% Gaussian noise. The data is simulated with Nt = 256 detector bins, Nφ = 171
projections and the reconstruction grid-size is M = 256. In fig. 5 we compare
the two reconstruction methods for a split parameter choice of K = 10. The
choice of K = 10 is relatively low, which fits well with the phantom being highly
12

1

0.8

0.6

0.4

0.2

0

Figure 4: Fibre crack-phantom with fibres along the direction 20◦ and cracks in a
circular pattern to the left. Simulated sinogram with no noise shown to the right.
directional. A comparison between u, w or u+w and ground truth z is not suited
for picking λw since similarity between these signals is not we the model tries to
achieve, nor will this yield the most desirable results. When using the variational
reconstruction method for sinogram splitting we therefore pick the optimal λu
and λw based on visual inspection instead. The priority for u was to achieve
clear edges and an intensity range close to the one in the ground truth. The
priority for w was to achieve a reconstruction with a homogeneous background
and sharp crack edges.
In fig. 5 we see that both fibre-component (u) reconstructions are visually similar, but the colorbar shows that the intensity level in the FBP reconstruction has
an offset of arround 0.5. If we examine the details, the variational reconstruction
has much sharper edges and less artifacts than the FBP result. For the crackcomponent (w) reconstructions we see that the noise is very dominating in the
FBP reconstruction, whereas the noise is removed by the use of TV-regularization.
In general we see that all of the cracks are located in the crack-components and
not present in the fibre-components, which is due to a highly directional object
and a good choice of the range-width index K = 10. We even see that the
cracks along the main direction are present in crack components. This is due
to the edges of the cracks that are perpendicular to the main direction. These
edge-singularities will be present in bw and hence the cracks are present in the
reconstruction.
To show the role the parameter K we have in fig. 6 compared reconstructions
using the variational approach for different values of K. The test shows that a too
low choice of K will result on some fibre elements in the crack-component, whereas
too high a choice of K will result in some of the cracks being reconstructed in
the fibre-component. The choice of this parameter should be chosen according
to prior knowledge about the object, i.e. if the object is highly directional a
relatively low value will be sufficient.

5.2

DTV-decomposition

From the definition of DTV and the empirical tests in [16] we know that the
choice the width parameter a is related to how well the directional prior fits the
given function we try to model. When we use DTV on the crack component
w we do not expect all of the crack to follow the direction completely, so we
13

0.5

1

0.8

0.6
0
0.4

0.2

-0.5

0

1

0.8
0.7
0.6

0.5
0.5
0.4
0.3

0

0.2
0.1
-0.5

0

Figure 5: Comparison between two reconstruction methods for the sinogram splitting decomposition technique introduced in section 3.
therefore suggest to chose the parameter aw higher than au , but still lower than
1, to avoid ending up with TV. The value should be chosen according the prior
given about the object of interest and for the crack phantom in fig. 4 we set
aw = 0.5. Based on this choice, together with the already fixed au = 0.15, we get
the bound on: 0.15 < α < 2. We pick λ based on a comparison between the sum
of the components u + w and the ground truth z. We choose the value of λ that
maximizes the psnr-value.
We tested the influence of the balance between the DTV-terms α and demonstrated some results in fig. 7. To avoid that the sparsity constraint will influence
the results we fix β = 10−6 , which is relatively low, in this test. The visual
interpretation of the results in fig. 7 tells us that a low α value will result in
more cracks in the the crack component, but also start to introduce noise and
fibre-artifact, whereas a high α-value will introduce a lot of the cracks in the
fibre-component and hence be less desirable, from a decomposition point of view.
We observe that the visual optimal choice of α also coincides with the maximum
psnr-value based on this observation we used the maximum psnr-value for picking
the optimal α-value.
We demonstrate the improvement of including the sparsity constraint for the
DTV-decomposition method by comparing reconstructions with β = 10−6 and
β = 10−4 . We visualize the the results in fig. 8. From the results we see a clear
improvement of both components. The intensity range for the fibre-component
is much more accurate and cracks have much sharper edges. The improvement is
also reflected by a slight increase of the psnr-value.

5.3

Sinogram Splitting vs DTV-decomposition

The demonstration of the two reconstruction methods in the two previous sections
shows that the sinogram splitting method delivers a much more complete split
14

1

1

1

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0

0.8

0.8

0.8

0.7
0.6

0.6
0.5

0.6

0.4
0.4

0.4

0.2

0.3
0.2

0.2

0

0.1
0

Figure 6: Parameter-test for the sinogram splitting method using DTVregularized reconstruction. Reconstructions for different choices of the parameter
K.
between the fibres and the cracks along any given direction. For the sinogram
splitting method the crack-component is seen to have a non-homogeneous background, and some ’fibre-artifacts’ are present in the regularized reconstruction.
Furthermore the cracks appear to be a bit wider than they should be, which is
due to the missing data and the strong emphasis which put on the regularizer to
remove limited data artifacts and noise. The DTV-decomposition method is best
at splitting the fibres from the cracks that are in the range [θ + 45◦ , θ + 135◦ ] or
[θ − 45◦ , θ − 135◦ ]. When the sparsity constraint is enforced the we see that the
background of the crack component is highly homogeneous, while the appearing
crack-edges are still sharp. Depending of the object and the analysis task on
method does not seem superior to the other when testing with this simulated
crack-phantom.
Based on the conclusions from the simiulated phantom test we have chosen to
compare the sinogram splitting method with the DTV-decomposition method on
a real sample object. The carbon fibre sample we have chosen for the comparison
can be seen in fig. 9. The sample is taken from [28] where it is analyzed, see
page 200-226. We have chosen this sample since it is unidirectional in most of the
sample, and that most of the cracks are perpendicular to the main fibre-direction.
Based on the sample we simulate data-acquisition as previously and we further
simulate noise by adding 1% Gaussian noise to the sinogram. The sample is of
size M = 426 and we simulate Nt = 426 detector bins and Nφ = 284 projections.
In order to have a common reference for the two decomposition methods we
have reconstructed the sample using three method that does not decompose the
object: FBP, `2 − TV and `2 − DTV. The reconstructions can be seen in fig. 10.
For the sinogram splitting we have tested both FBP and DTV for the reconstruction. We show both the fibre and the crack component for each reconstruction in fig. 11. For all three reconstruction methods we tuned the parameters
by visual inspection, where we prioritize a decomposition of cracks in one com15

1

1

1

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0

0.3
0.2

0.4
0.2
0.1

0.2

0.1

0
0

0

-0.1
-0.2

-0.1

-0.2
-0.3
-0.2
-0.4

-0.4

Figure 7: Parameter test of the DTV balance parameter α for the DTVdecomposition method. Reconstructions for different choices of α visualized.
ponent and non-cracks in the other. Furthermore we tuned the parameters to
achieve a homogeneous background for the crack component without smoothing
away small detail of the cracks. For the sinogram splitting reconstructions the
range-width index K is much larger than for the crack-phantom. This choice is
made to avoid having the parts with non-crack singularities, that does not follow
the main direction, appearing in the crack-component. By choosing a large K we
force the crack-component to have a more homogeneous background.
The split-FBP result is seen to be highly influenced by noise in both components and also limited angle artifacts occur. Compared with the two other
reconstruction results, the split-FBP method must be seen as inferior.
The split-DTV result has a sharp edges in the fibre-structure along the main
direction, while other edges are blurry, especially the ones perpendicular to the
main direction. On the other hand the cracks in the central part of the object are
not present in the fibre component as desired. The crack component for the splitDTV method have clearly marked cracks with sharp edges, but this component
suffers a lot from stair-casing effects which makes the background highly nonhomogeneous. This makes the crack-detection unreliable, since it not certain
whether the reconstructed cracks in this component are indeed cracks.
The DTV-DTV reconstruction result has sharp edges for the fibre-component
and small-detail parts that could be categorized as cracks are moved to the
crack component. The crack-component has a homogeneous background and
sharp crack-edges. Some of the very narrow cracks are not present in the crackcomponent, but can be observed in the fibre-component.
For this real sample with fibres along one main direction, cracks mostly perpendicular to the fibre-direction and some irregular pieces in bottom right part of
sample, the DTV-decomposition has produced the overall best decomposition and
reconstruction. The edges are sharper in DTV-decomposition result compared
to the split-DTV result, for both components. A disadvantage of the DTVdecomposition is that the narrow crack are not present in the crack-component,

16

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0
0

0.3
-0.1
0.2
0.1

-0.2

0

-0.3

-0.1

-0.4

-0.2

-0.5

-0.3

-0.6

-0.4

-0.7

Figure 8: Comparison of the DTV-decomposition method low emphasis on the
sparsity constraint, β = 10−6 , and with higher emphasis on the sparsity constraint, β = 10−4 . α = 0.7 and λ = 0.0038. Reconstruction results are visualized.
but these were already difficult to reconstruct from the underdetermined system
with noisy data, cf. fig. 10.

6

Conclusions

We have proposed two new tomographic reconstruction methods that makes utilize variational formulations. Both reconstruction methods combine reconstruction and decomposition into one problem and aim to solve both simultaneously.
We compare the two method by discussing their theoretical differences. We have
also proposed and demonstrated a new method for estimating the main object
direction directly from measured computed tomography data. The combined reconstruction and decomposition methods are compared empirically for both a
simulated and real data sample. The simulated phantom tests serves as a general
performance tests of the methods. In these tests we demonstrate what can be
achieved with the proposed methods. The real data sample tests show how well
these methods perform in practice. In order achieve even better reconstruction results we consider if the two combined decomposition and reconstruction methods
can be combined into a single method.

Acknowledgements
The work was supported by Advanced Grant 291405 from the European Research
Council.

17

Figure 9: Carbon fibre sample. Reconstruction from CT scan of carbon fibre.
Taken from [28].

Figure 10: Simulated CT-problem based on real data object. Solved using three
different reconstruction techniques.

References
[1] G. Aubert and P. Kornprobst. Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations. Springer
Science + Business Media, 2006.
[2] J. F. Aujol, G. Aubert, L. Blanc-Féraud, and A. Chambolle. Image decomposition into a bounded variation component and an oscillating component.
J. Math. Imaging Vis., 22(1):71–88, 2005.
[3] J. F. Aujol, G. Gilboa, T. Chan, and S. Osher. Structure-texture image
decomposition-modeling, algorithms, and parameter selection. Int. J. Comput. Vis., 67(1):111–136, 2006.
[4] T. M. Buzug. Computed Tomography : From Photon Statistics to Modern
Cone-Beam CT. Springer, 2008.
[5] A. Chambolle and P.-L. Lions. Image recovery via total variation minimization and related problems. Numer. Math., 76(2):167–188, 1997.

18

[6] A. Chambolle and T. Pock. A first-order primal-dual algorithm for convex
problems with applications to imaging. J. Math. Imaging Vis., 40(1):120–
145, dec 2011.
[7] A. H. Delaney and Y. Bresler. Globally convergent edge-preserving regularized reconstruction: An application to limited-angle tomography. IEEE
Trans. Image Process., 7(2):204–221, 1998.
[8] S. Esedoglu and S. J. Osher. Decomposition of images by the anisotropic
Rudin-Osher-Fatemi model. Commun. Pure Appl. Math., 57(12):1609–1626,
2004.
[9] J. Frikel and E. T. Quinto. Characterization and reduction of artifacts in
limited angle tomography. Inverse Probl., 29(12):125007, 2013.
[10] J. Gilles. Noisy image decomposition: A new structure, texture and noise
model based on local adaptivity. J. Math. Imaging Vis., 28(3):285–295, 2007.
[11] M. Holler and K. Kunisch. On infimal convolution of total tariation
type functionals and applications. SIAM J. Imaging Sci. J. Imaging Sci.,
7(4):2258–2300, 2014.
[12] K. M. Jespersen, J. Zangenberg, T. Lowe, P. J. Withers, and L. P. Mikkelsen.
Fatigue damage assessment of uni-directional non-crimp fabric reinforced
polyester composite using X-ray computed tomography. Compos. Sci. Technol., 136:94–103, 2016.
[13] K. M. Jespersen, J. Zangenberg, T. Lowe, P. J. Withers, and L. P. Mikkelsen.
X-ray CT Data: Fatigue Damage in Glass Fibre/polyester Composite Used
for Wind Turbine Blades [Data-set], 2016.
[14] E. Jonsson, T. Chan, and S.-C. Huang. Total variation regularization in
Positron Emission Tomography. Technical report, Dept. Mathematics, Univ.
California., Los Angeles., 1998.
[15] R. D. Kongskov and Y. Dong. Directional total generalized variation regularization for impulse noise removal Rasmus. Scale Sp. Var. Methods Comput.
Vis. 2017, 6667:221–231, 2017.
[16] R. D. Kongskov, Y. Dong, and K. Knudsen. Directional total generalized
variation regularization. http://arxiv.org/abs/1701.02675, 2017.
[17] V. P. Krishnan and E. T. Quinto. Handbook of Mathematical Methods in
Imaging. Springer Science + Business Media, 2015.
[18] E. Levitan, J. Degani, and J. Zak. Regularization in Fourier-synthesis 2-D
image reconstruction. IEEE Trans. Nucl. Sci., NS-26(3):4327–4329, 1979.
[19] J. Li, C. Miao, Z. Shen, G. Wang, and H. Yu. Robust frame based x-ray ct
reconstruction. Journal of Computational Mathematics, 34(6):683, 2016.
[20] Y. Long and J. A. Fessler. Multi-material decomposition using statistical
image reconstruction for spectral CT. IEEE Trans. Med. Imag., 33(8):1614–
26, Aug. 2014.

19

[21] A. K. Louis and F. Natterer. Mathematical problems of computerized tomography. Proc. IEEE, 71(3):379–389, 1983.
[22] Y. Meyer. Oscillating Patterns in Image Processing and Nonlinear Evolution
Equations, volume 22. American Mathematical Society, 2001.
[23] F. Natterer. The Mathematics of Computerized Tomography, volume 29.
John Wiley & Sons Ltd. and B G Teubner, Stuttgart, 1986.
[24] W. J. Palenstijn, K. J. Batenburg, and J. Sijbers. Performance improvements
for iterative electron tomography reconstruction using graphics processing
units (GPUs). J. Struct. Biol., 176(2):250–253, 2011.
[25] E. T. Quinto. Singularities of the X-ray transform and limited data tomography in R2 and R3 . SIAM J. Math. Anal., 24(5):1215–1225, 1993.
[26] J. Radon. Über die Bestimmung von Funktionen durch ihre Integralwerte
längs gewisser Mannigfaltigkeiten. Akad. Wiss., 69:262–277, 1917.
[27] M. Romanov, A. B. Dahl, Y. Dong, and P. C. Hansen. Simultaneous tomographic reconstruction and segmentation with class priors. Inverse Problems
in Science and Engineering, 24(8):1432–1453, 2016.
[28] J. E. Rouse. Characterisation of Impact Damage in Carbon Fibre Reinforced
Plastics by 3D X-Ray Tomography. PhD thesis, University of Manchester,
2012.
[29] L. I. Rudin, S. Osher, and E. Fatemi. Nonlinear total variation based noise
removal algorithms. Phys. D Nonlinear Phenom., 60(1-4):259–268, 1992.
[30] S. R. Sandoghchi, G. T. Jasion, N. V. Wheeler, S. Jain, Z. Lian, J. P.
Wooler, R. P. Boardman, N. K. Baddela, Y. Chen, J. R. Hayes, E. N. Fokoua,
T. Bradley, D. R. Gray, S. M. Mousavi, M. N. Petrovich, F. Poletti, and D. J.
Richardson. X-ray tomography for structural analysis of microstructured and
multimaterial optical fibers and preforms. Opt. Express, 22(21):26181, 2014.
[31] O. Scherzer. Handbook of Mathematical Methods in Imaging. Springer Science + Business Media, 2010.
[32] E. Y. Sidky and X. Pan. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization. Phys. Med.
Biol., 53(17):4777–4807, sep 2008.
[33] J. L. Starck, M. Elad, and D. L. Donoho. Image Decomposition via the
combination of sparse representation and a variational approach. IEEE
Trans.Image Process., 14(10):1570–1582, 2005.
[34] W. van Aarle, W. J. Palenstijn, J. Cant, E. Janssens, F. Bleichrodt,
A. Dabravolski, J. De Beenhouwer, K. Joost Batenburg, and J. Sijbers. Fast
and flexible X-ray tomography using the ASTRA toolbox. Opt. Express,
24(22):25129, 2016.
[35] W. van Aarle, W. J. Palenstijn, J. De Beenhouwer, T. Altantzis, S. Bals,
K. J. Batenburg, and J. Sijbers. The ASTRA Toolbox: A platform for
advanced algorithm development in electron tomography. Ultramicroscopy,
157:35–47, 2015.

20

0.8
0.6

0.6
0.4

0.4
0.2
0.2
0
0
-0.2
-0.2
-0.4
-0.4
1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

1
0.1
0.8

0.6

0
-0.1
-0.2

0.4
-0.3
0.2

-0.4
-0.5

0

Figure 11: Simulated CT-problem based on real data object. Solved and simultaneously decomposed using three different reconstruction techniques.

21

