A practical local tomography reconstruction algorithm based on known subregion

We propose a new method to reconstruct data acquired in a local tomography setup. This method uses an initial reconstruction and refines it by correcting the low frequency artifacts known as the cupping effect. A basis of Gaussian functions is used to correct the initial reconstruction. The coefficients of this basis are iteratively optimized under the constraint of a known subregion. Using a coarse basis reduces the degrees of freedom of the problem while actually correcting the cupping effect. Simulations show that the known region constraint yields an unbiased reconstruction, in accordance to uniqueness theorems stated in local tomography.


Region of Interest Tomography
Region of Interest (ROI) tomography, also called local tomography, naturally arises when imaging objects that are too large for the detector field of view (FOV), as depicted on Figure 1. It notably occurs in medical imaging, where only a small part of a body is imaged. Local tomography can also originate from a radiation dose concern in medical imaging. Image: (Zeng, 2010) Since the projection data does not cover the entire object, it is said to be truncated with respect to a scan that would cover the entire object. The aim is then to reconstruct the ROI from this "truncated" data.
However, due to the nature of the tomography acquisition, the acquired data is not sufficient to reconstruct the ROI in general: for each angle, rays go through the entire object, not only the ROI. Thus, the data does not only contain information on the ROI, but also contribution from the parts of the object external to the ROI. For example, on Figure 1, the detector gets data from parts of the object located at the left of the ROI. These contributions from the external parts actually preclude from reconstructing exactly the ROI from the acquired data in general.
The problem of reconstructing the interior of an object from truncated data is referred as the interior problem. It is well known that the interior problem does not have a unique solution in general. If P denotes the projection operator, d the acquired data and x a solution of the problem P (x) = d, then x is defined up to a set of ambiguity functions u such that P (x + u) = d. An example is given in (Clackdoyle & Defrise, 2010) where u is non-zero in the ROI, but P (u) = 0 in the detector zone corresponding to the ROI : two solutions differing by u would produce the identical interior data. In (Wang & Yu, 2013), it is emphasized that the ambiguity is an infinitely differentiable function whose variation increases when going outside the ROI. The nonuniqueness of the solution of the interior problem prevents quantitative analysis of the reconstructed slices.
Methods tackling the ROI tomography problem can mainly be classified in two categories. The first category methods aim at completing the data by extrapolating the sinogram. There are often oriented toward easy and practical use, although having no theoretical guarantees. The second category of methods rely on prior knowledge on the object. Many theoretical efforts were made on these methods, providing for example uniqueness and stability results.
Other works use wavelets to localize the Radon transform (Rashid-Farrokhi et al., 1997) (Sastry & Das, 2005) or focus on the detection of discontinuities, the best known being probably Lambda-tomography (Yu & Wang, 2006).

Sinogram extrapolation methods
In a classical tomography acquisition, the whole object is imaged. If nothing is surrounding the object, the rays are not attenuated by the exterior of the object ; thus the sinogram values for each angle go to zero on the left and right parts (after taking IUCr macros version 2.1.10: 2016/01/28 the negative logarithm of the normalized intensity). In a local tomography acquisition, however, the data is "truncated" with respect to what would have been a whole scan.
The incompleteness of the data induces artifacts on the reconstructed image. The first obvious artifact is visible as a bright rim on the exterior of the image. This bright rim is the result of the abrupt transition in the truncated sinogram: the filtration process suffers from a Gibbs phenomenon. Another artifact is referred as the cupping effect: an unwanted background appears in the reconstructed image, which makes further analysis like segmentation challenging. These two artifacts occur simultaneously, but they have different causes. The bright rim comes from the truncation, while the cupping comes from the contribution of the external part. Figure 2 and 3 illustrates these artifacts. A synthetic slice is projected, and the resulting sinogram is truncated to simulate a ROI tomography setup. The filtering step enhances the transition between the ROI and the truncated part which is set to zero. The difference between the filtered whole sinogram and the filtered truncated sinogram also shows the cupping effect, which appears as a low-frequency bias.
is a homogeneous polynomial of degree n in sin θ and cos θ, for all n ≥ 0. An alternative formulation is given by equation (2) : for k > n ≥ 0 and k−n even. In (Van Gompel et al., 2006), (2) is used as a quantitative measure of the sinogram consistency, and is optimized as an objective function.
For many applications, constant extrapolation provides acceptable results (Kyrieleis et al., 2011), although cupping artifact makes the segmentation challenging.

Prior knowledge based interior tomography
It was long believed that ROI tomography cannot be solved exactly, because of the nature of Radon inversion through FBP: the reconstruction of each voxel requires the knowledge of all the (complete) lines passing through this voxel. However, in the last decade, it has been shown that multiple nonequivalent reconstruction formulas allow partial reconstruction from partial data in the 2D case (Clackdoyle & Defrise, 2010).
Alternatively to Filtered Back Projection reconstruction, which requires complete data, Virtual Fan Beam (VFB) and Differentiated Back-Projection (DBP) were developed based on the Hilbert projection equality (Clackdoyle et al., 2009).
Moreover, uniqueness theorems based on analytical continuation of the Hilbert Transform were stated and progressively refined in (Noo et al., 2002), , , (Zou et al., 2005), , (Ye et al., 2007), , (Tang et al., 2012). They ensure an exact and stable reconstruction of the ROI given some assumptions. These assumptions can be of geometric nature, or in the form of a prior knowledge. the ROI on only one side (Ye et al., 2007).
These geometry-based methods do not work, however, when the FOV does not extend the object (Figure 4 (b)). In this case, it has been shown ) that a prior knowledge on the function inside the ROI enables exact and stable reconstruction of the ROI. This knowledge can be in the form of the function values inside a sub-region of the ROI  or can be about the properties of the function to reconstruct, for example sparsity in some domain.
This latest kind of knowledge has led to compressive sensing based ROI tomography.
In (Yang et al., 2010), (Yu & Wang, 2009), (Lee et al., 2015), Total Variation method is used to reconstruct the ROI. In (Niinimäki et al., 2007), the function is assumed to be sparse in the wavelet domain, and a multi-resolution scheme reduces the number of unknown by keeping only fine-scale wavelet coefficients inside the ROI. In (Klann et al., 2015), it is shown that piecewise constant functions are determined everywhere by their ROI data, the underlying hypothesis being formulated as sparsity in the Haar domain.

Low frequencies artifacts correction with Gaussian blobs
Sinogram extrapolation usually copes well with the correction of discontinuities in the truncated sinogram, but does not correct the cupping effect in general. This cupping effect appears as a low frequency bias in the reconstructed image. Sinogram extrapolation and other background correction techniques do not give guarantees that the low frequency bias will actually be removed without distorting the reconstruction.
In this section, we describe a new method using prior knowledge on a subregion of the reconstructed volume to eliminate the low frequency cupping bias. The starting point of this method is an initial reconstruction, hereby denoted x 0 , which can be obtained for example with the padded FBP method. This initial reconstruction is then refined with an additive correction term. This correction term uses the known sub-region as a constraint which should be sufficient, according to uniqueness theorems stated in the references given in 1.3, to accurately reconstruct the region of interest.
As x 0 bears the high frequencies features, the correction term is expressed as a linear combination of Gaussian functions to counterbalance the low frequency artifacts. The coefficients are optimized subject to the knowledge of the subregion, hereby denoted Ω. To constrain all the Gaussian coefficients by the knowledge of the image values in Ω, a reduced set of coefficients is firstly computed inside Ω. Then, the Gaussian coefficients are iteratively optimized to fit the reconstruction error of the whole image, using the coefficients computed inside Ω as a constraint.
Let u 0 denote the "true" object values in the known region Ω. The proposed method can be summarized as follows: • The reconstruction error in Ω, denoted e |Ω = (x 0 ) |Ω − u 0 , is expressed as a linear combination of two dimensional Gaussians. The resulting Gaussian coefficients are denoted g 0 .
• The error in the whole image is iteratively fitted with Gaussians coefficients g, subject to g |Ω = g 0 , to build a consistent reconstruction error in the whole image.
Details of each step are described in the following parts.

Capturing the low frequencies of the error in the known zone
The key assumption of this method is that Ω is large enough to bear sufficient information on the low frequencies artifacts (cupping effect) of a classical reconstruction.
The reconstruction error in Ω is approximated as a linear combination of Gaussian functions. This function is chosen for computational convenience, more details are given in 2.2.3.
Equation (3) gives the expression of the approximation.The error estimation e |Ω is a linear combination of translated Gaussians of weights c i,j , with a spacing s.
For simplicity, all the Gaussians have the same standard deviation σ and all have the same spacing s between them. Their location on the image grid is also fixed, so that the fit turns into a linear inverse problem : where G is the operator taking as an input the coefficients c i , stacked in a vector g ; and producing an image tiled with Gaussians (here in region Ω). The norm · 2 2 is the squared Frobenius norm, that is, the sum of the squares of all components. Choosing the same standard deviation σ and the same spacing s for all Gaussians enables to implement G as a convolution. More precisely, given an image being zero everywhere, coefficients c i are placed every s pixels in Ω and this image is convolved with two dimensional ψ σ .
The reconstruction error in Ω is thus estimated in the least squares sense: g 0 is the vector of Gaussian coefficients giving the best estimation e |Ω of the reconstruction error in the L2 sense. This vector g 0 will be used in the second part of the algorithm.

Correcting the reconstruction error in the image
2.2.1. Overview of the method Outside the known region Ω, the reconstruction error is not known. Like some other methods described in 1.3, this algorithm aims at using the known region information to accurately reconstruct the whole ROI. However, this approach focuses on correcting an initial reconstruction: the reconstruction error in Ω is fitted by as a linear combination of Gaussians, then the whole image is corrected in a coarse Gaussian basis whose coefficients are constrained in the known subregion.
The Filtered Backprojection (FBP) with sinogram extrapolation is widely used in local tomography because it is both simple and gives satisfactory results in general (Kyrieleis et al., 2011). Theoretical investigations found that FBP provides a reconstructed function bearing the same discontinuities as the reference function (Bilgot et al., 2009) , although the cupping effect can make the segmentation challenging. In this method, FBP with padding is used to obtain an initial estimate of the reconstruction ; the aim is to correct the local tomography artifacts on this image using the prior knowledge. Equation (5) gives the expression of the estimate x where x 0 is the initial reconstruction, G is the operator described in 2.1 andĝ is a linear combination of Gaussian functions aiming at counterbalancing the low frequencies artifacts.
x =x 0 + Gĝ The vectorĝ is found by minimizing an objective function which is built as follows. A new image x =x 0 + Gg, containing the initial reconstruction, is created. The imagẽ x 0 is an extension of the initial reconstruction x 0 . This new image is projected with a projector P adapted to the bigger size. To compare with the acquired sinogram d, the computed sinogram P (x 0 + Gg) is truncated by a cropping operator C. The data fidelity is then given by Equation (6).
We emphasize that this approach differs from the full estimation of the ROI based on a subregion. The variables g are in a coarse basis whilex 0 is fixed, which is notably reducing the degrees of freedom of the problem. The operation Gg has two goals: a coarse estimation of the exterior (outside the x 0 support) and a correction of the low frequencies error inside the x 0 support.
As the minimization is on g, the initial estimate of the ROI x 0 is constant, and the data fidelity term (6) can be rewritten as in Equation (7) 1 where d e = d −P x 0 is the difference between the acquired sinogram d and the projection of the initial reconstruction x 0 , andP is the projector adapted to the size of x 0 . The optimization problem is given by Equation (8). This local constraint is propagated in all the variables by the projection operator involved in the optimization process. Uniqueness theorems mentioned in 1.3 state that the knowledge of a subregion of the ROI is sufficient to yield an exact reconstruction.
However, when using a pixel basis without space constraints, the number of degrees of freedom might be too high ; leading to a slow convergence. Using a coarse basis for correcting the low frequencies reduces this number of degrees of freedom.

Details on the involved operations
In this part, more details are given on the different steps of the algorithm. We start by computing a padded FBP reconstruction x 0 which gives an initial estimate of the ROI of size (N, N ). This image is extended to a bigger imagex 0 of size (N 2 , N 2 ) where N 2 > N , and x 0 is placed in the center of the image.
At each iteration k, the image Gg k , where g k is the Gaussian coefficients vector at iteration k, is computed. The operator G consists in placing the coefficients g on a regular grid and convolving with the Gaussian kernel (x, y) → 1 Thus, the operator G can be written G = C σ U s where C σ is the convolution by the aforementioned Gaussian kernel of standard deviation σ, and U s is an operator upsampling an image by a factor of s. As both are linear operators, G is a linear operator and G T = U T s C T σ where U T s is the s-downsampling operator and C T σ is a convolution by the matched Gaussian kernel, which is the same kernel due to symmetry. In our implementation, the Gaussian kernel has a size of 8σ + 1 , i.e the Gaussian is truncated at 4σ. The resulting image is given by Equation (9) x where ψ σ is the discrete Gaussian kernel, z is the image containing the Gaussian coefficients g placed on the grid of size (N 2 , N 2 ) after upsampling, that is, zeros almost everywhere except coefficients every s pixel. The summation in Equation (9) is done on the convolution kernel support. If s < 4σ, the Gaussian functions supports can overlap once placed on the grid. In practice, these Gaussians should overlap to appropriately fit constant regions : for s close to σ, the Gaussians almost yield a partition of unity (Bale et al., 2002).
Once the coefficients are placed on a grid and convolved by the 2D Gaussian function, the image x is projected. The projection operator adapted to the new geometry (the bigger image Gg k ) is denoted P . This is a standard Radon transform. This process is illustrated in Figure 5.
As the new image x is bigger than x 0 , the sinogram P x and the acquired data d cannot be directly compared. The computed sinogram P x is thus cropped to the region corresponding to the ROI. The cropping operator is denoted by C ; it is also a linear operator whose transpose consists in extending the sinogram by inserting zeros on both sides.
The resulting sinogram aims at fitting the error between the acquired sinogram d and the (cropped) projection of the object. As the object is unknown except inside Ω, the reconstruction error is only known in Ω. The Gaussian coefficients g are constrained by those found by fitting the error inside Ω in 2.1. The projection operator involved in the process propagates the constraint to all the other coefficients. Coefficientsĝ from Equation (8) are computed with an iterative solver. As this objective function is quadratic, efficient minimization algorithms like conjugate gradient can be used. The final image is obtained withx = x 0 + Gĝ and is cropped to the region of interest.

Computational aspects
Using Gaussians as functions to iteratively express the error has several computational advantages. Gradient-based algorithms for solving (8) involve the computation of the forward operator P G and its adjoint G T P T . They are usually the computationally expensive steps of iterative solvers. In this case, these operators can be computed in an efficient way.
The Gaussian kernel has an interesting property: it is the only (nonzero) one to be both rotationally invariant and separable (Kannappan & Sahoo, 1992). In our case, the convolution by a Gaussian followed by a projection (forward Radon transform) is equivalent to projecting first and convolving each line of the sinogram by the corresponding one dimensional Gaussian, as illustrated on Figure 5.
The first advantage is of theoretical nature. In many implementations, the projector and backprojector pair are usually not adjoint of each other for performances reasons.
Although giving satisfying results in most practical applications, this raises theoretical issues on convergence of algorithms using iteratively forward and backward operators (Zeng & Gullberg, 2000). By using a point-projector and a point-backprojector implementation, the pair can be exactly adjoint, besides giving more accurate results.
The second advantage is on the computational side. As the operator P G consists in projecting 2D Gaussians disposed on the image, it is equivalent to placing one-pixel coefficients (Dirac functions in the continuous case) in the image on a grid denoted I, projecting the image (with a point-projector) and convolving the sinogram by a one-dimensional Gaussian kernel. The same goes for the adjoint operator G T P T consisting in retrieving the Gaussians coefficients from a sinogram. The standard way to compute this operator would be backprojecting the sinogram, convolving by the two dimensional Gaussian (which is its own matched filter due to symmetry), and sampling the image on the grid I to get the coefficients. Here, the convolution can be first performed in one dimension along the sinogram lines. The sinogram is sampled at locations corresponding to points I in the image domain. The resulting sampled sinogram is then backprojected with a point-backprojector.

Pseudocode of the proposed method
In this section, the different steps of the proposed method are summarized in two algorithms. The first performs the fitting of the error in the known zone as described in section 2.1, the second builds the resulting image as described in section 2.2.
A complete implementation of the proposed method is available at (Paleo, 2016).
It contains comments on the different steps and can be tuned for various setups. This implementation relies on the ASTRA toolbox (Palenstijn et al., 2011) (van Aarle et al., 2015, the point-projector scheme described in 2.2.3 is not implemented for readability ; but this approach would be more suited to a production reconstruction algorithm where performances are an issue. In practice, the final imagex is cropped to the region of interest. In algorithm 2, the optimization (line 5) can be done with a gradient algorithm, as differentiating the quadratic error term requires only the operators and their adjoints.

Results and discussion
In this section, results and discussions on three test cases are presented. Synthetic sinograms are generated by projecting an object and truncating the sinogram to the radius of a given region of interest in the image.
The following notations are used: σ is the standard deviation of the Gaussians of the basis, s is the grid spacing, N 2 is the size (width or height in pixels) of the extended image and R is the radius (in pixels) of the known region. In practice, the size of the "original image" (which corresponds to the size of an image that would contain the whole object in practice) is unknown, hence N 2 is always chosen different from the width of the original test image.
In all cases, the input image is projected with a projector covering the entire object.
The resulting sinogram is then truncated to the radius of the region of interest. The truncated sinogram is the input of the methods. The proposed method is compared to the padded FBP. As the padded FBP is used as an initial reconstruction by the proposed method, the benchmark is mainly about checking that the cupping effect is actually removed, and that the correction does not induce distortion to the final image.
The first test involves the standard Shepp-Logan phantom (Figure 6   The fact that the reconstruction has the same mean values than the true interior could enable quantitative analysis of the reconstructed volume, which is not easily achievable in local tomography. The second test involves the test image "Lena", 512 × 512 pixels, bearing both smooth regions and high frequencies textures. Figure 10 shows the test setup. The known region has be chosen as slowly varying as possible, as in real acquisitions the known region is likely to be air or coarse features. The width of the extended image is N 2 = 520.  Figure 11 shows the difference between the true interior and the reconstruction with the proposed method, with varying values of the radius R of the known region. The parameters σ = s = 3 has been used for these reconstructions. As it can be expected, the cupping effect removal is better when the known region is wide. It is also interesting to visualize the reconstruction of the whole extended image. As it can be seen on Figure 13, the Gaussian basis even yields an approximation of the exterior. This approximation is actually important for modeling the contribution of the external part in the acquired sinogram. The bias correction is thus closely related to the modeling of the external part.  As a last remark, Figure 17 shows the result of this method without using the known zone constraint, that is, without applying the constraint g |Ωg = g 0 in (8). As expected, there is a not-null mean bias, even if it has been reduced with respect to padded FBP. Beside visual inspection, reconstructions can be quantitatively compared to the true interior of the test image. Table 1 shows the comparison with two image metrics: peak signal to noise ratio (PSNR) and the structural similarity index (SSIM). As these metrics are indicators of an average distance between two images, we believe it is well suited for this purpose of evaluating how the low frequencies are corrected by the proposed method. These results suggest that the proposed method yield better overall reconstruction quality than padded FBP. In particular, it does not induce critical distortion in the reconstruction.
The following advantages of this method can be highlighted. Using a Gaussian basis allows for efficient computation: the correction can be implemented as a convolution with a simple kernel. This basis can also be extended to a multi-resolution grid, where big Gaussian functions are placed outside the ROI and small Gaussian functions are placed inside the ROI, to reduce the degrees of freedom even further. This approach would be similar to the method proposed in (Niinimäki et al., 2007), although here only the correction is expressed in a multi resolution grid, not the image variables.
What can also be noted is that no assumption is done on the shape or location of the ROI and the known region: the known region can for example be several regions of various shapes, corresponding to pores in a sample. Finally, by using the correction in the forward model (8) can be achieved (Bale et al., 2002). Thus, the final reconstruction cannot be exact due to the basis coarseness, but can provide results quite close to FBP with full data as seeing Figures 8, 12 and 16.
The fact that only a known zone of the ROI is enough to guarantee an almost-exact reconstruction might be counterintuitive, especially in our case where this constraint is expressed in a coarse basis. In the Gaussian basis, local constraints are propagated to the global image by the projection and backprojection operators involved in the process. Using a coarse basis greatly reduces the degrees of freedom of problem (8).
The classical tomographic reconstruction problem P x = d turned into a least squares optimization argmin x P x − d 2 2 is ill-posed, even for complete data. In a local tomography setup, the ill-posedness is even worse (Clackdoyle & Defrise, 2010). Iterative solvers dealing with problem (10) argmin for x in the pixel space, have very slow convergence in general due to the high number of degrees of freedom, even with spatial constraints. The importance of reducing the degrees of freedom of (10) is highlighted for example in (Niinimäki et al., 2007) and

Conclusion
We presented a new technique of local tomography reconstruction based on the knowledge of a zone of the region of interest. This technique corrects the cupping effect in an initial reconstruction by expressing the error in a coarse basis of Gaussian functions.
In accordance to local tomography uniqueness theorems, this method yield almost exact reconstructions, in spite of being only a correction with a coarse basis. Besides, practical considerations are given for an efficient implementation suitable for reconstruction of real data. A commented implementation of this method can be found at (Paleo, 2016).