1 Introduction
Xray Computed Tomography (CT) is a highly used noninvasive 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 microscale (nanometers). The technique is based on the attenuation that Xrays undergoes as they pass trough matter. The attenuated Xray intensity is measured at a detector. In this paper our CT setup is 2D and parallel beam, but the methods introduced here are not limited to this CTmethodology, but presented in this setting for simplicity. The Xray source and detector is rotated around the object, and for each angle of rotation the attenuated Xray signal is measured.
LambertBeers law relates the object that we are scanning to the measured intensity. By rewriting LambertBeers law, the measured intensity is seen to be linearly related to the object through lineintegrals [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 matrixvector product as
(1) 
Here the system matrix, , models the length of each Xray 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 CTproblems was introduced in [26] and is now known as Filtered BackProjection (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 illposedness, 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
(2) 
The first term, the datafidelity term, is chosen to be the 2norm 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 postprocessing, 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 DTVdecomposition, 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 by pixels with equidistant pixelgridsize . 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 objectreconstruction process. Since the publication of the paper that introduced the ROFmodel 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 piecewise constant with the texture along one main direction. In this paper we will use DTV for regularizing the CTreconstruction problem. In discrete terms we can express DTV as
for scaling and rotation matrices
(3) 
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 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, , with 1% Gaussian sinogramnoise. We compare the results produced by FBP,  and  by visualizing the results and by means of the peak signaltonoise 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 CTdata. Furthermore a widthparameter should be chosen, but based on the conclusions from empirical tests in [16] we just set in this paper.
2.2 Direction estimation from CTdata
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 basisfunctions with same main direction, will be relatively large in magnitude.
To find the maindirection 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 detectorpixels. 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 basisfunctions, 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.
It should be noted that this angle estimation method is not limited to parallelbeam tomography since a pattern also occurs when the fanbeam geometry is used, although the data should be rebinned 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 realdata 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 noiselevel () and we report the angleestimates in table 1. The two tested samples together with their corresponding noisefree sinograms are shown in fig. 2 and fig. 3. In these two figures we also indicated the angle estimated from the noisefree sinogram and the angledependent 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 
Based on these empirical tests we see that our method is robust up until a noiselevel of at least 20%, which is a high noiselevel for CTdata.
3 Sinogram Splitting
Edges are essential in the analysis of CTscans 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 Xray 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 Xrays 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 sinogramdata is split into two parts, according to the measurement angles, such that each part will include information about singularities along different directions. For a fibrecrack 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 CTmodel in separate parts, where each part is related to an objectcomponent where the directions of the singularities, in the objectdomain, are limited. The simplest split is a twocomponent 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,
(4) 
where and are the object components. and are the splitmatrices. 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 subproblem 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 rangewidthindex , 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 crackcomponent, which could complicate a complete fibrecrack 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 subproblems 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 nonlinear 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 twocomponent 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 0data. 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 noisecorrupted data.
Variational methods has been shown to be advantageous for CTreconstruction problems with limited and/or corrupted data. The motivation for this approach is to solve the subproblems using prior information about the individual components. So for each component we pick a suitable regularization method and solve the corresponding subproblem. 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 fibrecrack decomposition, the prior for the fibrecomponent is that it is piecewise constant and that the texture follows one main direction. The prior for the crackcomponent is that it is independent of the direction, piecewise constant and sparse. Based on these priors we suggest to regularize the fibrecomponent using DTV and the crackcomponent using classical TV together with a norm. In discrete terms we have
Both of these regularization choices are convex and due to the piecewise 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 rangewidthindex should be chosen to match the expected range of the fiber direction, i.e. small for a highly onedirectional 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 backprojection opposed to several hundred, or thousand, forward and backprojections 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 DTVdecomposition 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 texturecartoon decomposition, see [3, 10], we opt to build a variational model where the minimizer will be a decomposed reconstruction.
In many applications fibrestructures 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:
(5) 
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 objectdirection , to be piecewise constant and to be nonnegative since it corresponds to attenuation coefficient values. is assumed to follow the orthogonal main direction , to be piecewise 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 noiserobust algorithm, see section 2.2. The first DTVwidth which we assume to be fixed due to previous empirical tests, see [16]. The balance of DTVterms which should be chosen in combination with the second DTVwidth to achieve sufficient splitting between the two components. The orthogonal main direction which is just orthogonal to the main direction. And the cracksparsity 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
(6) 
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 norminequalities 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 DTVdecomposition, 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 DTVdecomposition where the sum of the components gives the object. The DTVdecomposition method assumes that the crackcomponent follows a specific direction, where the sinogram splitting method includes no such assumption, only an assumption on the fibrecomponent.
For the special case where and 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 DTVdecomposition 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 Xray 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 PrimalDualHybridGradient 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 signaltonoise 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 crackphantom ground truth can be seen in fig. 4. This phantom does not necessarily fit our assumptions, but serves as a stresstest 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 dataset, which is corrupted with 1% Gaussian noise. The data is simulated with detector bins, projections and the reconstruction gridsize 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 fibrecomponent () 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 () reconstructions we see that the noise is very dominating in the FBP reconstruction, whereas the noise is removed by the use of TVregularization. In general we see that all of the cracks are located in the crackcomponents and not present in the fibrecomponents, which is due to a highly directional object and a good choice of the rangewidth 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 edgesingularities will be present in and hence the cracks are present in the reconstruction.
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 crackcomponent, whereas too high a choice of will result in some of the cracks being reconstructed in the fibrecomponent. 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 DTVdecomposition
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 psnrvalue.
We tested the influence of the balance between the DTVterms 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 fibreartifact, whereas a high value will introduce a lot of the cracks in the fibrecomponent and hence be less desirable, from a decomposition point of view. We observe that the visual optimal choice of also coincides with the maximum psnrvalue based on this observation we used the maximum psnrvalue for picking the optimal value.
We demonstrate the improvement of including the sparsity constraint for the DTVdecomposition 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 fibrecomponent is much more accurate and cracks have much sharper edges. The improvement is also reflected by a slight increase of the psnrvalue.
5.3 Sinogram Splitting vs DTVdecomposition
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 crackcomponent is seen to have a nonhomogeneous background, and some ’fibreartifacts’ 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 DTVdecomposition 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 crackedges 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 crackphantom.
Based on the conclusions from the simiulated phantom test we have chosen to compare the sinogram splitting method with the DTVdecomposition 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 200226. 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 fibredirection. Based on the sample we simulate dataacquisition 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.
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.
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 noncracks 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 rangewidth index is much larger than for the crackphantom. This choice is made to avoid having the parts with noncrack singularities, that does not follow the main direction, appearing in the crackcomponent. By choosing a large we force the crackcomponent to have a more homogeneous background.
The splitFBP 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 splitFBP method must be seen as inferior.
The splitDTV result has a sharp edges in the fibrestructure 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 staircasing effects which makes the background highly nonhomogeneous. This makes the crackdetection unreliable, since it not certain whether the reconstructed cracks in this component are indeed cracks.
The DTVDTV reconstruction result has sharp edges for the fibrecomponent and smalldetail parts that could be categorized as cracks are moved to the crack component. The crackcomponent has a homogeneous background and sharp crackedges. Some of the very narrow cracks are not present in the crackcomponent, but can be observed in the fibrecomponent.
For this real sample with fibres along one main direction, cracks mostly perpendicular to the fibredirection and some irregular pieces in bottom right part of sample, the DTVdecomposition has produced the overall best decomposition and reconstruction. The edges are sharper in DTVdecomposition result compared to the splitDTV result, for both components. A disadvantage of the DTVdecomposition is that the narrow crack are not present in the crackcomponent, 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.
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. BlancFé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. Structuretexture image decompositionmodeling, algorithms, and parameter selection. Int. J. Comput. Vis., 67(1):111–136, 2006.
 [4] T. M. Buzug. Computed Tomography : From Photon Statistics to Modern ConeBeam 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 firstorder primaldual 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 edgepreserving regularized reconstruction: An application to limitedangle tomography. IEEE Trans. Image Process., 7(2):204–221, 1998.
 [8] S. Esedoglu and S. J. Osher. Decomposition of images by the anisotropic RudinOsherFatemi 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 unidirectional noncrimp fabric reinforced polyester composite using Xray computed tomography. Compos. Sci. Technol., 136:94–103, 2016.
 [13] K. M. Jespersen, J. Zangenberg, T. Lowe, P. J. Withers, and L. P. Mikkelsen. Xray CT Data: Fatigue Damage in Glass Fibre/polyester Composite Used for Wind Turbine Blades [Dataset], 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 Fouriersynthesis 2D image reconstruction. IEEE Trans. Nucl. Sci., NS26(3):4327–4329, 1979.
 [19] J. Li, C. Miao, Z. Shen, G. Wang, and H. Yu. Robust frame based xray ct reconstruction. Journal of Computational Mathematics, 34(6):683, 2016.
 [20] Y. Long and J. A. Fessler. Multimaterial 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 Xray 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 XRay 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(14):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. Xray 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 conebeam computed tomography by constrained, totalvariation 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 Xray 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.
Comments
There are no comments yet.