A phase ﬁeld approach to limited-angle tomographic reconstruction

The tomography of an object with limited angle can be addressed through Iterative Reconstruction Reprojection (IRR) procedure, wherein a standard reconstruction procedure is used together with a “ﬁltering” of the image at each iteration. It is here proposed to use as a ﬁlter a phase-ﬁeld — or Cahn-Hilliard — regularization interlaced with a ﬁltered back-projection reconstruction. This reconstruction procedure is tested on a cone-beam tomography of a 3D woven ceramic composite material, and is shown to retrieve a reconstructed volume with very low artifacts in spite of a large missing angle interval (up to 28%).


Introduction
Most of tomographic reconstruction commercial software are based on Filtered Back-Projection (FBP) algorithm (including Feldkamp-Davis-Kress variant for cone-beam geometry used in lab-scale tomographs) [11].Indeed, FBP is efficient and can be implemented on GPU's so that computation is very fast.However, to use this method, projection angles have to sample 360°(or 180°for a parallel beam) with an equi-angular spacing.In some cases, such as transmission electron microscopy tomography [9] or some specific applications of X-ray micro-tomography (X-ray µ-CT) [21], the experimental gantry does not permit the acquisition of a complete set of projection angles, a situation referred to as "limited-angle".This difficulty has led to the development of algebraic algorithms (such as ART, SART, SIRT, MART, etc. [11]).As early as in 1982, it was smartly proposed that FBP algorithms could also be used in such cases, at the expense of using the algorithm iteratively instead of just once when all angles are accessible.This last method was named Iterative Reconstruction-Reprojection (IRR) by Nassi et al. [16].
The spirit of the IRR method is to perform an FBP reconstruction initially based on arbitrary projection data (e.g.blank image) for missing angles.The reconstructed volume is then filtered so as to get rid of the nonphysical data (such as negative absorption coefficients).The resulting volume is then re-projected providing synthetic projections in the missing angle orientations.Using the latter, complemented with the originally acquired projections, as input to FBP allows an enhanced quality reconstruction to be obtained.Through iterations, this procedure provides a more and more accurate reconstruction which by construction matches the recorded projections, and provides plausible projection data for the missing angles.
The importance of taking into account a priori information to clean out nonphysical data in the reconstructions given by FBP, has been noted early.Incorporation of those priors can be considered as a form of regularization that provides a more likely substitute for the missing data and speeds up convergence.
• Medoff et al. [14] give examples of constraints which are split into two categories: (i) constraints on the value of the reconstructed volume data f or of its sinogram s (e.g.positivity of f and s, the fact that all gray levels are in a known range or that a part of f is known), and (ii) constraints on the geometry of the image (e.g. the image is confined within known boundaries).
• Heffernan and Robb [10] have shown on artificial images that a regularization which sets to zero negative values of f and all the pixels of f located outside known boundaries, improves the reconstruction quality.
• Riddell et al. [20] also use constraints on the positivity and on the contour of the image and extrapolate sinogram data to missing angles to initialize IRR.They employ attenuated Radon transform to ensure the quantitative reconstruction.
• Ollinger [17] proposes a regularization based on expectation maximization to improve the robustness of IRR with respect to noise.The Poisson character of the noise affecting the projection is taken into account, thereby giving a fair estimate of the weight given to each data point of the projection data based on its uncertainty.
• Duan et al. [7] add a Total Variation (TV) minimization on the image domain.This procedure introduced in Ref. [5] aims at concentrating gray level variations onto a sparse support, leading to a small number of phases (discrete histogram).The key feature of TV minimization, is that convexity is preserved, and hence the minimizer is unique, and not dependent on the pathway used to the solution, thus allowing for efficient and robust implementation.
Up to now, most studies have been carried out on artificial images (such as the famous Shepp-Logan phantom) in 2D parallel or fan-beam geometry.Using such algorithms on experimental data, requires to account for (and correct) image imperfections and artifacts because the reprojection step has to provide data that is fully consistent with the reconstruction.In the present study, based on lab-scale tomography, the use of wide distribution of energies in the X-ray source produces naturally a significant artifact known as "Beam Hardening" (BH).As the absorption coefficient of materials is generally higher for low energy X-rays, the beam spectrum concentrates on higher energies as it travels through the part, and the Beer-Lambert law (used to interpret the radiography) appears to be only a rough approximation.A usual remedy to account for BH, is to identify an effective Beer-Lambert attenuation law, so that the projection data can be preprocessed with a non-linear grey level correction, from which a standard reconstruction algorithm can be run.
For complex materials and industrial parts, the geometry of the object to scan is often quite complex and generally only a model of the geometry (say a CAD description) is known and, due to the fabrication or elaboration processes, it has to be considered as no more than an approximation.It is often safer to base the regularization on the microstructure, say the homogeneous phases composing the object of interest (which are often just a few).An original regularization approach based on a phase field model is presented in this article.
Phase field models were initially developed to describe the physics of phase separation by free energy minimization [4,1].This global method has also demonstrated its efficiency in image segmentation [12] where the "free energy" now incorporates some known features or properties of the image itself.For instance, the grey level of a pixel can be considered as a parameter which characterizes its phases.In this context, phase field models can be used to very effi-ciently de-noise images, based on a prior knowledge with is included in the very definition of the free-energy.Regularization is expected to lower noise, filtering out high frequencies, however, it is desirable that it does not blur boundaries.Thus ideally, one would anticipate domains of constant grey levels separated by sharp interfaces according to regularization parameters.Those considerations was the motivation for TV filtering, however the latter approach does not exploit the existence of a limited number of phases with say a well defined absorption coefficient.The phase field regularization approach has this potential, from the flexibility offered by the very definition of the "free energy".It is proposed here to use a form of the free energy that includes a phase information enriched by a gradient term for boundaries.Let us stress that such a filtering of images require that no spurious gradients are produced by say BH artifacts, and hence such phenomena have to be cleared from the projection data prior to reconstruction.
If tomography was first developed for medical imagery, today its usage has diffused to a much broader scope of applications, including materials science.X-ray µ-CT is a valuable tool to study the microscopical structure and understand the behaviour of complex materials, such as composites.It is today a technique of choice for non-destructive evaluation.The sample chosen to illustrate this analysis is a Ceramic Matrix Composite (CMC) [15,6] whose outstanding thermo-mechanical properties make it a very promising key material to be used in hot regions of an aircraft engine.Matrix and fibers are made of the same material and to distinguish the texture, an excellent quality tomography is imperatively needed.Moreover, in situ thermo-mechanical testing is extremely desirable to detect the first fore-signs of damage to assess the performance of such parts [13].A mechanical testing machine requires to balance say a tensile force on the sample by an equivalent compressive one that has to be sustained by a rigid frame around the object.For heavy loads, a steel frame is well suited mechanically, but it will partly mask the specimen as it rotates in front of the beam, its shadow on the detector renders inexploitable a domain of the sinogram and thus induces reconstruction issues equivalent to those encountered in "limited-angle" tomography.Therefore, IRR is used to alleviate the difficulty of these missing data, yet imaging faithfully the faint contrast of the material microstructure is very demanding, and thus BH corrections and a rich prior filter (based on phase-field approach) are needed to extract an accurate 3D reconstruction.The present paper aims at exploring the ultimate performance of such an approach.
The principle of IRR implementation is recalled in Section 2 where notations are defined.In Section 3, the sample used as support is presented.Section 4 explains the BH correction strategy.Section 5 details the regularization step based on a phase field model.In Section 6, a reconstruction computed from real experimental data is discussed to illustrate the capacity of the algorithm to provide quality reconstruction with low limited-angle artifacts.Finally, conclusions and further possible improvements are proposed.

Implementation of IRR
For lab-tomography, radiographies of a specimen rotating for a full 360°rotation between an X-ray source considered as punctual, and a detector are acquired.The intensity of these radiographies is denoted as I(r, θ), where r is the 2D detector coordinates, θ is the sample rotation.The intensity of the beam I 0 (r) without sample is also recorded.From these data, the "projections" defined as can be obtained, and they correspond to the line integral of the specimen X-ray attenuation f (x) along the ray hitting the detector at position r, for a sample rotation of θ.Ideally, all rotation angles, 0 ≤ θ < 2π, and all relevant detector positions r are accessible, so that a reconstruction algorithm R acting on p provides the reconstructed volume f (x), where x stands for 3D coordinates of voxels Mathematically, R is a mere linear operator, which in practice is never written explicitly because of its huge size, and a numerical procedure such as the Filtered Back-Projection algorithm is used as equivalent to R. The inverse operator, is the projection P In the case of obstruction, only some pairs of (r, θ) are accessible forming a set Σ m for which p(r, θ) = p m (r, θ) where the superscript m stands for "measured".
The IRR algorithm is based on the following scheme illustrated in Figure 1: • When the missing data has been complemented with inconsistent data, the resulting f usually contains unphysical features such as negative absorption coefficients.Moreover, the effect of the missing angle data is is spread out over the entire reconstruction domain.A filtering operation F is introduced that makes the data more acceptable, by projecting negative values of f to 0, and cancelling f to 0 outside a known or guessed support.The cleaned reconstruction is denoted as g = F[f ]; • Projecting the cleaned reconstruction provides the computed projection s(r, θ) = P[g]; • Because of the filter, and the fake data used in the missing domain, s does not coincide with p m on the measured points.However, in the complementary domain, s is much more satisfactory than the initial guess.Thus, the complement operation S is introduced so that • As shown in Figure 1, these four operations define a cycle, that produce a physically more consistent missing data in p.This full cycle can thus be repeated as a fixed point algorithm until the filtering or complement operations do not change appreciably their output data.
Iterative Reconstruction-Reprojection algorithm principle: Starting from projection data p(r, θ), the FBP reconstruction algorithm R provides a reconstructed volume f (x).The latter is filtered through a regularization procedure, F, to obtain the volume g(x).The projection operator, P, gives projections s(r, θ).The procedure S restores the projections for the known angles, and leaves s(r, θ) for the missing angles.This loop corresponds to one iteration of the IRR algorithm, which is repeated until convergence to a fixed point For this algorithm to converge, it is essential that the filter F provides a physically realistic information that is consistent with the known projections.The positivity of f or the prescription of its support are very generic filters that were introduced very early.They do not constraint much the reconstruction f and hence the convergence rate is quite slow.
It is worth noting an interesting parallel that can be drawn with regularization.In a seminal paper, Candès et al. [5] proposed to perform tomography with a very coarse angular sampling.They showed that a drastic reduction in the number of projections could be performed without loss on the quality of the reconstruction for the Shepp-Logan phantom.This can be seen also as a kind of missing information on the angles not used as an input.The algorithm they proposed can be cast into the scheme of Figure 1, where the "filter" they proposed was a Total Variation (TV) minimization.The latter promotes uniform absorption coefficient within domains and the sparsity of the support of the non-zero gradient of f is the key underlying assumption to complement the missing information.It is to be noted that the example chosen for illustration was a numerical example deprived of any artifact.
Thus filtering or regularization can be considered as a very powerful lever to extend the performance of tomography, and in the following a novel phasefield regularization will be used.However, it is essential that no model error is present which would render the regularization assumption invalid, ruining the convergence of the IRR procedure.
IRR makes use of two standard primitives for tomography: FBP reconstruction R and (re-)projections P. Here, the free software ASTRA [23,22] is chosen.It offers a python and matlab toolbox with many such flexible primitives wrapped in an efficient GPU-implementation [19].The two transformations R and P are resp.the ASTRA functions FKD3D and FP3D.
A tomographic scan often represents a huge amount of data (more than 10 9 pixels for the full sinogram).To reduce the size of the problem, the spatial resolution of projection can be degraded, so that the computation is faster.Reconstruction and sinogram obtained at the coarse scale is then used as initialization for finer resolution computation.Finally a multiscale algorithm may be designed by a recursive use of this coarsening step, in a pyramidal approach [8].

CMC sample
To illustrate the addressed problem, a SiC-SiC Ceramic Matrix Composite (CMC) sample is chosen.It has a rather faint texture with short wavelength details corresponding to the 3D woven fabric, and some Si inclusions.This rather weak contrast makes the problem quite difficult to handle.Moreover, because the chemical species are light elements the X-ray absorption is low and gives rise to a significant beam-hardening.The sample (Fig. 2a) has a rectangular section with both machined and rough surfaces so that the capacity of the algorithm to reconstruct edges can be demonstrated.For the same reason, a groove was added on one side leaving a very thin wall.This case is thus quite challenging.
The composite sample is placed in a mechanical testing device which presents two columns on the side of the sample (Fig. 3).Those columns (that withstand the compressive load balancing the tensile one applied to the sample) cross the beam line during the tomographic acquisition.Thus angular data are missing to perform the reconstruction.
To quantify the ability of the proposed algorithm results to match the scanned sample, the projections used in the following were taken without the columns.Missing data are numerically suppressed on the measured sinogram to initialize the computation.Hence the ground truth reconstruction is known and errors can be trustfully quantified.

Principle
As earlier mentioned, in a lab-tomograph, the broad spectrum of energy of the X-ray source alters the effective absorption of the beam as it passes through the sample.However, when the nature of the constitutive phases of the studied sample is comparable in terms of attenuation, then the original projection data p can be corrected by a non-linear function of the path length ℓ in the material, and hence of p itself to restore an effective Beer-Lambert law.
Algorithmically, one can infer this correction in the following way.Starting from uncorrected projection data, a first reconstruction is performed and a binary mask h(x) of the studied sample is created (with an Otsu filter [18]) where h = 1 within the material, and 0 elsewhere.The projection of this binary mask P[h] gives precisely the length of the intersection of the beam with the sample ℓ(r, θ).Collecting data for all (r, θ), a large collection of (p, ℓ) data points is obtained from which a non-linear regression can be made to estimate the correspondence p = ψ(ℓ) Then, applying the reverse transformation p eff = ψ −1 (p), restores the validity of Beer-Lambert law for p eff .

Application of BH correction
In the studied example, because the material is made mostly of a single phase, SiC, such a beam correction is well justified.Figure 4 shows the raw data points and their fit with the following expression that revealed adequate, ψ(ℓ In the reconstruction, beam hardening appears as an halo around the sample as shown in figure 2a and an overestimate of X-ray absorption near the edge.This artifact is particularly detrimental for image processing based on gray levels such as TV regularization or the approach proposed in the following that would tend to display several phases according to the distance to boundary.After correcting the projection data, the reconstruction shown in (Fig. 2c) is obtained.As compared to the raw reconstruction (Fig. 2a), a much more uniform texture is obtained as could have been anticipated from the large scale homogeneity of the material.

Phase field regularization
In the regularization formulation, prior knowledge about the microstructure is to be inserted to help reaching a physically admissible reconstructed volume.For instance, TV aims at reducing the number of phases (i.e., gray levels in the reconstruction) but only through a penalization of phase boundaries.In fact, no information is provided on the nature of the different phases.Such an approach can easily be generalized through the formulation of a free energy density φ[g(x)], so that the minimization of the total free energy, Φ[g] = φ[g(x)] dx, aims at bringing the reconstruction closer to the actual microstructure.The free energy can be designed so that the microstructure lies at or close to a minimum.Thus for a known initial reconstructed volume, f (x), the filtered image, g(x), is such that it minimizes where the first part of the functional is an "attachment to data" term designed to maintain a close connection between f and g, whereas the second term embodies the prior knowledge on g.For instance, if the free energy density is simply the L 1 norm of the gradient of g, then plain TV regularization is recovered.Another very generic model in this description introduces a double well potential φ(g) = (g − γ 1 ) 2 (g − γ 2 ) 2 together with a L 2 norm of the gradient of g.The former term favors two "pure phases", g = γ 1 or g = γ 2 , and the latter through the gradient term induces a smooth transition (resulting from surface tension) between those values wherever needed.Actually such a form is a standard "phase-field", or Cahn-Hilliard model, introduced in the context of statistical physics [4], and then frequently used to describe phase morphologies [2].It has been exploited more recently very successfully in image processing (for denoising or inpainting) [12,3] that is quite close to the presently proposed usage.
The general form of Φ offers much more flexibility than the above cited examples and allows one to tailor the potential to suit the specimen at hand.A natural strategy is to perform an expansion of φ based on the assumption that the regularity of g is mostly local.At dominant, or 0 th order, φ depends of the local value of g only.At first order, a local nearest neighbor information can be formulated by introducing the gradient of g.Order after order, higher derivatives can be inserted.When inserted into a multiscale approach, non-local information about features having a specific shape can be introduced.

Gray level corrections based on histogram
A very first exploitation of this technique is based on a 0 th order free energy density.In order to limit the introduction of a prior knowledge, (the absolute value of the coefficient of absorption is fragile as discussed in Section 4), a first determination of the reconstructed image f may be used to estimate the values of favored phases through the histogram of gray levels.Let n(f ) be the histogram of f values, a natural idea is to propose a free energy density defined as φ(f where µ is a free parameter to tune the regularization term. Minimizing the free energy defined from this density tends to favor gray values that are already well represented in f .Indeed, all maxima of n will remain invariant under the transformation, and will behave as attractors for the corrected value g, whereas minima of n will be repulsive fixed points (as it can be seen on an example in Fig. 5a).
Because of this observation, the binning of the histogram is very important.For too fine a binning, spurious peaks will hamper the regularization effect, and for a too coarse one, information may be lost by artificially grouping different phases into one.The ideal bin size is given by the noise (or artifact) level which limits the identification of the different phases.
As φ only depends on the gray level value, it is straightforward to solve Eq. 6.In fact, it leads to a correction that only depends on the gray value and not on the voxel location.More precisely, differentiating Eq. 6 leads to The solution to this equation can be written as g = ω(f ) and is seen as a color correction where the initial gray level f is re-encoded with a different value, ω(f ), that solves the above equation.Thereby, a mere look-up table can be written to represent efficiently this filter.In practice, the weight µ is adjusted so that the correspondence ω(f ) fulfills a stable convergence towards the minima of n and hence −1 < dω/ df < 1 at those points, or µ < 2/(max(|n ′′ |)).This regularization actually reinforces the peaks of the histogram and depletes the population of other gray values, thereby selecting stable phases.One may also note that such an effect can easily be strengthened by iterating the gray level mapping, g = ω α (f ), with α > 1.In the limit of a very large α, only "pure" phases are obtained.

Gradient term
The free energy can also be enriched by a gradient term such as the addition of ξ β |∇f | β where ξ and β are tunable parameters.ξ is homogeneous to a length scale, to be interpreted as the width over which the transition between two stable gray levels takes place.A large value of ξ will smooth out (or "blur") boundaries, and a small one will keep transitions rather sharp.The exponent β also plays a role on the sharpness of the transition: When β ≤ 1 discontinuous transitions are encouraged through a localization instability (instability is due to a loss of convexity of this part of the potential when β < 1), while β = 1, aka TV regularization, is just the marginal case, allowing sharp interfaces, but still being convex, a unique combination that made the success of this form).
Here β = 2 is chosen, and equation 8 can be generalized to This equation is to be solved after each reconstruction.Adding this gradient term allows new peaks to emerge in the histogram and to discriminate phases with neighboring gray levels.

Regularized IRR computation
To provide a reference, a first computation is performed on the complete sinogram.Only one IRR loop is performed because, as there is no missing domain, the sinogram is not updated during transformation S (see Fig. 1).Parameters of regularization α and ξ are chosen such that they lower noise without altering the texture of the material.This choice is made by studying the difference between the corrected and the original images, and targeting parameters that give rise to mostly noise in this difference, and no feature that could be read as microstructure.The same parameters will be used to perform IRR on missing-data sinogram so that the results can be compared.Parameters α = 3, µ = 0.2µ max and ξ = 1 voxel are a good compromise for the studied image.

Limited-angle
On the sinogram, the traces of the columns of a mechanical testing machine are now inserted as shown in Figure 6a, and two scenarii are considered to deal with the data they hide.First, referred to as missing data domain A, all radiographs where at least one column can be seen are rejected.The resulting mask is shown in Figure 6b   When compared to the reference volume where no obscuration takes place, it can be observed that gray levels are spatially modulated (left part is slightly darker).Some streaks aligned in the direction of the missing angles are also visible.
The norm of the difference between reconstructed volumes with full or limited angle, (resp.f 0 and f A ), normalized by the norm of the reconstructed volume, is f A − f 0 / f 0 = 5.5%.
These results can also be compared with non-regularized IRR (Fig. 8).In this case convergence is reached in 6 iterations, and the residual relative error is r = 12.7%.Fig. 8b shows that missing angles are responsible for marked streaks that cannot be erased.Because of the filtering stage, each iteration is more costly for the phase field filter, but this is compensated by a faster convergence, and a better quality of the result.

Missing-data domain B
In the considered problem, one can easily exploit all the projection data that is not masked by the columns of the testing machine, as illustrated in Figure 6c.The IRR algorithm is quite similar to the previous case.Only the shape of the sinogram missing-data domain is needed, to define which part is updated and which part, Σ m , is preserved along iterations.
The results of regularized IRR (shown in Fig. 9), are very close to the reconstruction with the complete sinogram data.A relative error of r = 2.8%, is reached after 3 iterations.

Conclusion
A new and efficient regularization, based on a phase field approach, has been implemented in IRR algorithm and has proven its ability to correctly reconstruct volume with large missing-data domains.The few parameters (µ, α, β, ξ) of this phase-field-based regularization are all means to adapt this algorithm to the specificity of the scanned object and can be set to values optimized from a limited volume where they can be tuned at low cost.The counterpart of the missing data is the additional computational time of the regularization step and of the needed iterations.Computation costs were not reported here since the phase field filter was not optimized, but this regularization is essentially a very short range problem so that a GPU implementation should provide very efficient solutions, so that the incurred cost of the missing angle should reach almost the number of iterations times the cost of a full reconstruction.
As explained in part 6.2.2, the proposed formulation for IRR is versatile enough to take into account a missing-data domain of arbitrary shape.For instance, the proposed IRR scheme is directly applicable for any metal artifact corrections.

Figure 2 :Figure 3 :
Figure 2: (a) Sample cross-section as reconstructed from a complete acquisition and no BH corrections; (b) Binary mask obtained with an Otsu filter [18] used to compute the actual length ℓ traversed by each ray; (c) Reconstruction after BH correction where a much more uniform texture is observed; (d) Difference between (c) and (a) showing the net effect of BH corrections: the texture is not much impacted, however the correction produces a better uniformity close to the edges

Figure 4 :
Figure 4: Joint distribution of projection gray level p and path length ℓ plotted as thin black dots in the background.Superimposed on the data points, the bold red curve shows the regression that provides the BH correction ψ −1 function.

Figure 5 :
Figure 5: (a) Example of phase field regularization for α = 1 and µ = (1/2)µ max .The histogram peak sharpens around the maximum.Voxels corresponding to void are not considered in the regularization and hence only pixels with a gray level higher than a threshold (here 0.33) are used to compute the histogram; (b) Initial cross-section image; (c) Effect of the regularization on the image.
, (θ ∈[15°;64°]∪[205.25°;258.75°],or 28% of the data from the complete sinogram are lost).If only the region obscured by the columns is masked, missing data domain B, as shown in Figure6c, only 9% of the data is discarded.Both of these strategies can be handled with IRR, and the results are shown thereafter.

( a )
Sinogram of the sample on the experimental device, shadows of columns are black (b) Missing-data domain A: all radiographs where the columns can be seen are discarded: 72% of the full data range is used (c) Missing-data domain B: Only the region obscured by the columns are discarded: 91% of the full data range is known

Figure 6 :Figure 7
Figure 6: Definition of the missing-data domain

Figure 7 :
Figure 7: Example of IRR results for missing projections.Step 4 is the final step when a stationary solution is reached.The residual is the difference between the considered step of the reconstruction and the one without missing data (after BH correction)

Figure 8 :
Figure 8: Comparison between the IRR results for 28% of missing angles with and without regularization.Convergence is reached in 4 (resp.6) iterations with (resp.without) regularization

Figure 9 :
Figure 9: Example of IRR results for missing-data domain B. Step 3 is the final one