Tomographic Reconstruction Methods for Decomposing Directional Components

by   Rasmus Dalgas Kongskov, et al.

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.



There are no comments yet.


page 7

page 13

page 14

page 15

page 16

page 18

page 21


Directional Global Three-part Image Decomposition

We consider the task of image decomposition and we introduce a new model...

Space and camera path reconstruction for omni-directional vision

In this paper, we address the inverse problem of reconstructing a scene ...

Kernel Smoothing, Mean Shift, and Their Learning Theory with Directional Data

Directional data consist of observations distributed on a (hyper)sphere,...

General Support-Effective Decomposition for Multi-Directional 3D Printing

We present a method to fabricate general models by multi-directional 3D ...

CodEx: A Modular Framework for Joint Temporal De-blurring and Tomographic Reconstruction

In many computed tomography (CT) imaging applications, it is important t...

The GLD-plot: A depth-based plot to investigate unimodality of directional data

A graphical tool for investigating unimodality of hyperspherical data is...

Image Analysis Using a Dual-Tree M-Band Wavelet Transform

We propose a 2D generalization to the M-band case of the dual-tree decom...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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

equidistant pixels. The discrete model can be expressed as a simple matrix-vector product as


Here the system matrix, , models the length of each X-ray through each pixel in the reconstruction domain. The reconstruction, , is the object, including it’s interior, that we seek to find. The sinogram data, , is the measured data at the detector. The data is measured at a series of discrete angles, , where 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 . In order make the dimensions fit we have . The application of is usually called a forward projection, whereas the application of is a called a backward projection.

The process of constructing based on known 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


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 and a function is linearly split into two components as and the infimal convolution of the two functionals becomes

This strategy relies on the assumption that each component, and , 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, , 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 eq. 2 for and different regularization terms for and , 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 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 DTV-decomposition 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 -by- pixels with equidistant pixel-grid-size . In this region denotes a pixel index with such that gives the pixel value at . Since the scanning is circular, the region of interest, and our simulated phantoms, are contained in a circle with the diameter . The sinogram domain is discretized in -by-, so we have . In the sinogram domain denotes the pixel index such that is sinogram value at detector bin and measurement angle .

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 object-reconstruction 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 can be expressed as

for euclidian norm and the discrete gradient operator . is defined as

where and , for the two dimensions and , are obtained by applying a forward finite difference scheme with symmetric boundary condition, i.e.,

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

for scaling and rotation matrices


where and .

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 uni-directional 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, , with 1% Gaussian sinogram-noise. We compare the results produced by FBP, - and - 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 -TV and -DTV methods are tuned to maximize peak-signal-to-noise (psnr) ratio. DTV parameters are chosen as and .

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 should be chosen, but based on the conclusions from empirical tests in [16] we just set in this paper.

2.2 Direction estimation from CT-data

In order to use DTV for regularization we have to select the parameters and . From [16] we know that if the directional prior is met we should choose as small. Empirical tests in [16] show that 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.

1:Input siogram data and discrete measurement angles .
1D Fourier transform along each angle :
2:Calculate summed magnitude of for each angle and find max response:
return .
Algorithm 1 Main Direction Estimation Method

It should be noted that this angle estimation method is not limited to parallel-beam 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 from algorithm 1.

(%) 0 1 3 5 10 20 30 40
Phantom 20.1 20.1 20.1 20.1 20.1 20.1 20.1 31.7
Real 81.5 81.7 81.5 80.9 81.7 79.5 -1.1 -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 (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 (see 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.

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 and singularities in the corresponding Radon transform is described. The paradigm that is described in [25] is further outlined in [17] as follows:

detects singularities of perpendicular to the line of integration (”visible” singularities), but singularities of in other (”invisible”) directions do not create singularities of 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 consecutive measurement angles, , whereas the second part is based on the remaining angles, . This splitting results in two linear systems,


where and are the object components. and are the split-matrices. And and 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 . Further more is required to have a reasonable 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 and even range-width-index , we can pick .

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 eq. 4. If we use the same linear reconstruction method on both sub-problems we will have have 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 eq. 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 eq. 4 this corresponds to the assumptions: and . 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 . 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

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

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 -norm. In discrete terms we have

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 should be chosen to match the expected range of the fiber direction, i.e. small for a highly one-directional object and larger for an object with fibres within a range of angular directions. The regularization parameters and which control the balance between the fidelity terms and the regularizers. And -sparsity parameter which should be relatively small, compared to , in order to avoid .

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:


In this model represents the fibres and the crack part. The summation of these two components is assumed to model the object, hence the summation in the first term of eq. 5. is assumed to follow the main object-direction , to be piece-wise constant and to be non-negative since it corresponds to attenuation coefficient values. is assumed to follow the orthogonal main direction , to be piece-wise constant and to be sparse.

The parameters of the decomposition model eq. 5 are the following: The regularization parameter , 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 which we assume to be fixed due to previous empirical tests, see [16]. The balance of DTV-terms which should be chosen in combination with the second DTV-width 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 , which should be chosen as a relatively small value, compared to , to avoid forcing to be zero. This leaves four tunable parameter: , , and . In order to avoid an overlap of the feasible sets of the two DTV terms we have the following restriction for :

This bound is obtained by comparing the two DTV functionals. We want to avoid both


for all . Using the matrix definitions in eq. 3 and the fact that we can rewrite the two norm expression above and achieve:

where . We have , since we do not want DTV to become classical TV, i.e. . Based on the two norm-inequalities in eq. 6 we get

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

The model eq. 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 , it also makes eq. 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 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 fibre-component.

For the special case where and we have that the data-fidelity terms for the two methods are the same since

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-Hybrid-Gradient 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 . 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 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.

Figure 4: Fibre crack-phantom with fibres along the direction and cracks in a circular pattern to the left. Simulated sinogram with no noise shown to the right.

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 detector bins, projections and the reconstruction grid-size is . In fig. 5 we compare the two reconstruction methods for a split parameter choice of . The choice of is relatively low, which fits well with the phantom being highly directional. A comparison between , or and ground truth is not suited for picking 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 and based on visual inspection instead. The priority for was to achieve clear edges and an intensity range close to the one in the ground truth. The priority for was to achieve a reconstruction with a homogeneous background and sharp crack edges.

In fig. 5 we see that both fibre-component () 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 crack-component () 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 . 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 and hence the cracks are present in the reconstruction.

Figure 5: Comparison between two reconstruction methods for the sinogram splitting decomposition technique introduced in section 3.

To show the role the parameter we have in fig. 6 compared reconstructions using the variational approach for different values of . The test shows that a too low choice of will result on some fibre elements in the crack-component, whereas too high a choice of 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.

Figure 6: Parameter-test for the sinogram splitting method using DTV-regularized reconstruction. Reconstructions for different choices of the parameter .

5.2 DTV-decomposition

From the definition of DTV and the empirical tests in [16] we know that the choice the width parameter is related to how well the directional prior fits the given function we try to model. When we use DTV on the crack component we do not expect all of the crack to follow the direction completely, so we therefore suggest to chose the parameter higher than , 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 . Based on this choice, together with the already fixed , we get the bound on: . We pick based on a comparison between the sum of the components and the ground truth . 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 , 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.

Figure 7: Parameter test of the DTV balance parameter for the DTV-decomposition method. Reconstructions for different choices of visualized.

We demonstrate the improvement of including the sparsity constraint for the DTV-decomposition method by comparing reconstructions with and . 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.

Figure 8: Comparison of the DTV-decomposition method low emphasis on the sparsity constraint, , and with higher emphasis on the sparsity constraint, . and . Reconstruction results are visualized.

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 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 or . 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 Gaussian noise to the sinogram. The sample is of size and we simulate detector bins and projections.

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

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, and . The reconstructions can be seen in fig. 10.

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

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 component 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 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 we force the crack-component to have a more homogeneous background.

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

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 split-DTV method have clearly marked cracks with sharp edges, but this component suffers a lot from stair-casing effects which makes the background highly non-homogeneous. 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 crack-component, 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 DTV-decomposition is that the narrow crack are not present in the crack-component, 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.


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


  • [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-Fé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.
  • [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., 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.
  • [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 and . SIAM J. Math. Anal., 24(5):1215–1225, 1993.
  • [26] J. Radon. Über die Bestimmung von Funktionen durch ihre Integralwerte lä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.