Tomographic Reconstruction Methods for Decomposing Directional Components

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 (1) * Department of Applied Mathematics and Computer Science, Technical University of Denmark, 2800 Kgs. Lyngby, Denmark (rara@dtu.dk, yido@dtu.dk.
Here the system matrix, A ∈ R N ×M 2 , models the length of each X-ray through each pixel in the reconstruction domain. The reconstruction, z ∈ R M 2 , is the object, including it's interior, that we seek to find. The sinogram data, b ∈ R N , is the measured data at the detector. The data is measured at a series of discrete angles, φ ∈ R N φ , 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 N t . In order make the dimensions fit we have N = N φ N t . The application of A is usually called a forward projection, whereas the application of A T 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 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 J 1 and J 2 a function is linearly split into two components as z = z 1 + z 2 and the infimal convolution of the two functionals becomes This strategy relies on the assumption that each component, z 1 and z 2 , 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 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.

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 u i,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 N t -by-N φ , so we have b ∈ R Nt×N φ . In the sinogram domain (l, m) denotes the pixel index such that b l,m is sinogram value at detector bin t l and measurement angle φ m .

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 ∈ R M ×M can be expressed as for euclidian norm | · | 2 and the discrete gradient operator ∇ : where ∇ x1 and ∇ x2 , for the two dimensions x 1 and x 2 , 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 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, 2 3 , 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. 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.

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] 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: Calculate summed magnitude ofb for each angle φ m and find max response: 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.     [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).

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.

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 ∈ R K+1 , whereas the second part is based on the remaining angles, φ w ∈ R N −K−1 . This splitting results in two linear systems, where u ∈ R M 2 and w ∈ R M 2 are the object components. A u ∈ R (K+1)Nt×M 2 and A w ∈ R (N −K−1)×M 2 are the split-matrices. And b u ∈ R (K+1)Nt and b w ∈ 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 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 . 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.

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: A u w = 0 and A w 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 Both of problems are still limited angle problems, but now the assumptions are only imposed through the regularizers R u and R w . 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 1 -norm. In discrete terms we have R u (u) = λ u DTV(u) and R w (w) = λ w TV(w) + β w 1 .
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.

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 u ∈ R M 2 represents the fibres and w ∈ R M 2 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 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 a u 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 a w ∈ (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: λ, α, a w 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 |R θ Λ au y| 2 < α|R θ Λ aw y| 2 and α|R θ Λ aw y| 2 < |R θ Λ au y| 2 , for all y ∈ R 2 . Using the matrix definitions in (3) and the fact that θ ⊥ = θ + π 2 we can rewrite the two norm expression above and achieve: |R θ Λ au y| 2 = (y 1 cos θ + y 2 sin θ) 2 + a 2 u (−y 1 sin θ + y 2 cos θ) 2 α|R θ Λ aw y| 2 = (αa w ) 2 (y 1 cos θ + y 2 sin θ) 2 + α 2 (−y 1 sin θ + y 2 cos θ) 2 where y = (y 1 , y 2 ). We have a u , a w ∈]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 Two of the above inequalities contradict each other and the remaining two gives the bound for α 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 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 A u w = 0 and A w u = 0 we have that the datafidelity 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.

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

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 N t = 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  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 b w 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.

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 therefore suggest to chose the parameter a w higher than a u , 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 a w = 0.5. Based on this choice, together with the already fixed a u = 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.

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 [θ + 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 N t = 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 com-  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 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 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, but these were already difficult to reconstruct from the underdetermined system with noisy data, cf. fig. 10.

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.   Figure 11: Simulated CT-problem based on real data object. Solved and simultaneously decomposed using three different reconstruction techniques.