research papers
A practical local tomography reconstruction algorithm based on a known subregion
^{a}ESRF, 71 avenue des Martyrs, 38000 Grenoble, France, ^{b}Université de Grenoble, GIPSALab, 11 Rue des Mathématiques, 38400 SaintMartind'Hères, France, and ^{c}GIPSALab, Grenoble Images Parole Signal Automatique, Institut Polytechnique de Grenoble, France
^{*}Correspondence email: pierre.paleo@esrf.fr, mirone@esrf.fr
A new method to reconstruct data acquired in a local tomography setup is proposed. This method uses an initial reconstruction and refines it by correcting the lowfrequency artifacts, known as the cupping effect. A basis of Gaussian functions is used to correct the initial reconstruction. The coefficients of this basis are found by optimizing iteratively a fidelity term under the constraint of a known subregion. Using a coarse basis reduces the
of the problem while actually correcting the cupping effect. Simulations show that the known region constraint yields an unbiased reconstruction, in accordance with uniqueness theorems stated in local tomography.Keywords: tomography; reconstruction; algorithm; interior problem.
1. Introduction
In this section, we briefly recall the region of interest tomography problem and review the related work.
1.1. 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). 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.
Since the projection data do 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 the `truncated' data.
However, due to the nature of the tomography acquisition, the acquired data are not sufficient to reconstruct the ROI in general: for each angle, rays go through the entire object, not only the ROI. Thus, the data do not only contain information on the ROI but also contributions from the parts of the object external to 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 to 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 by Clackdoyle & Defrise (2010), where u is nonzero in the ROI but P(u) = 0 in the detector zone corresponding to the ROI: two solutions differing by u would produce identical interior data. Wang & Yu (2013) emphasize that the ambiguity is an infinitely differentiable function whose variation increases when moving 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 into two categories. The first category methods aim at completing the data by extrapolating the sinogram. These are often oriented toward easy and practical use, although having no theoretical guarantees. The second category of methods rely on prior knowledge of the object. Many theoretical efforts have been made on these methods, providing for example uniqueness and stability results.
Other works use wavelets to localize the Radon transform (RashidFarrokhi et al., 1997; Sastry & Das, 2005) or focus on the detection of discontinuities, the best known being probably Lambdatomography (Yu & Wang, 2007).
1.2. 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 the negative logarithm of the normalized intensity). In a local tomography acquisition, however, the data are `truncated' with respect to what would have been a whole scan.
The incompleteness of the data induces artifacts in 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 to as the cupping effect: an unwanted background appears in the reconstructed image, which makes further analysis such as 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.
Figs. 1 and 2 illustrate 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 lowfrequency bias.
Sinogram extrapolation methods primarily aim at eliminating the bright rim resulting from the truncation by ensuring a
between the ROI and the external part. In addition, efforts have been put into estimating the missing data in order to reduce the cupping effect. These techniques are referred to as sinogram extrapolation methods: the external part is estimated from the truncated data with some extrapolating function.The extrapolating function can be, for example, constant (the outermost left/right values are replicated), polynomial or cos^{2}. Shuangren Zhao & Yang (2011) used a mixture of exponential and quadratic functions to estimate the external part, possibly iteratively. The projection of a circle, for which a closedform formula is known, can also be used (Van Gompel, 2009). A common approach is using the values of the left/right part of the sinogram to estimate the external part, that is, replicating the borders values.
In general, sinogram extrapolation methods do not take into account the sinogram theoretical properties. For example, given an object being nonzero only inside a circle of a given radius, the sinogram decreases to zero at the left and right boundaries. Generally speaking, a sinogram of complete measurements satisfies the Helgason–Ludwig consistency conditions (1) (Van Gompel, 2009):
is a homogeneous polynomial of degree n in and , for all n 0. An alternative formulation is given by equation (2):
for k > n 0 and kn even. Van Gompel et al. (2006) used (2) 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 the cupping artifact makes the segmentation challenging.
1.3. Priorknowledgebased interior tomography
It was long believed that ROI tomography could be solved exactly, because of the nature of the Radon inversion through filtered backprojection (FBP): the reconstruction of each voxel requires 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 twodimensional case (Clackdoyle & Defrise, 2010). Alternatively to FBP reconstruction, which requires complete data, virtual fan beam (VFB) and differentiated backprojection (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 (Noo et al., 2002, 2004; Clackdoyle et al., 2004; Zou et al., 2005; Defrise et al., 2006; Ye et al., 2007; Kudo et al., 2008; Tang et al., 2012). These authors ensure an exact and stable reconstruction of the ROI given some assumptions. These assumptions can be of geometric nature, or in the form of prior knowledge.
Geometrybased prior knowledge is related to the acquisition geometry. For example, in DBPbased reconstruction, a point can be reconstructed if it lies on a line segment extending outside the object on both sides, and all lines crossing the segment are measured (Clackdoyle & Defrise, 2010). Similar results were obtained under less restrictive assumptions, for example the FOV extending the ROI on only one side (Ye et al., 2007).
These geometrybased methods do not work, however, when the FOV does not extend the object. In this case, it has been shown (Courdurier et al., 2008) that a prior knowledge of 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 subregion of the ROI (Kudo et al., 2008) 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. Yang et al. (2010), Yu & Wang (2009) and Lee et al. ( 2015) used the total variation method to reconstruct the ROI. Niinimäki et al. (2007) assumed the function to be sparse in the wavelet domain, and that a multiresolution scheme reduces the number of unknowns by keeping only finescale wavelet coefficients inside the ROI. Klann et al. (2015) showed that piecewise constant functions are determined everywhere by their ROI data, the underlying hypothesis being formulated as sparsity in the Haar domain.
2. Lowfrequencies 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 lowfrequency bias in the reconstructed image. Sinogram extrapolation and other background correction techniques do not give guarantees that the lowfrequency bias will actually be removed without distorting the reconstruction.
In this section, we describe a new method using prior knowledge of a subregion of the reconstructed volume to eliminate the lowfrequency cupping bias. The starting point of this method is an initial reconstruction, hereby denoted x_{0}, which can be obtained for example using the padded FBP method. This initial reconstruction is then refined with an additive correction term. This correction term uses the known subregion as a constraint which should be sufficient, according to uniqueness theorems stated in the references given in §1.3, to accurately reconstruct the ROI.
As x_{0} bears the highfrequencies features, the correction term is expressed as a linear combination of Gaussian functions to counterbalance the lowfrequency artifacts. The representation in a basis of Gaussian blobs in tomography was first introduced by Lewitt (1992); we use it here for representing coarse features. The coefficients are optimized subject to knowledge of the subregion, hereby denoted Ω. To constrain all the Gaussian coefficients by 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:
(i) The reconstruction error in Ω, denoted = , is expressed as a linear combination of twodimensional Gaussians. The resulting Gaussian coefficients are denoted g_{0}.
(ii) The error in the whole image is iteratively fitted with Gaussians coefficients g, subject to = g_{0}, to build a consistent reconstruction error in the whole image.
Details of each step are described in the following parts.
2.1. Capturing the low frequencies of the error in the known zone
The reconstruction error in Ω is first 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 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 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 twodimensional .
The reconstruction error in Ω is thus estimated in the leastsquares sense: g_{0} is the vector of Gaussian coefficients giving the best estimation of the reconstruction error in the L2 sense. This vector g_{0} will be used in the second part of the algorithm.
2.2. 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 as a linear combination of Gaussians, then the whole image is corrected on a coarse Gaussian basis whose coefficients are constrained in the known subregion.
The 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., 2011), 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 lowfrequencies artifacts,
The vector is found by minimizing an objective function which is built as follows. A new image x = , containing the initial reconstruction, is created. The image 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 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 while is fixed, which is notably reducing the 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 lowfrequencies 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),
where d_{e} = is the difference between the acquired sinogram d and the projection of the initial reconstruction x_{0}, and is the projector adapted to the size of x_{0}. The optimization problem is given by equation (8),
g_{0} is the vector of Gaussian coefficients found in (4), such that Gg_{0} approximates the error in the known zone. The set denotes the subset of the Gaussian basis corresponding to Ω in the pixel basis: if a coefficient c_{i} of g lies in in the Gaussian basis, then (Gg)_{i} lies in Ω in the pixel basis. Equation (8) boils down to finding coefficients g minimizing the reconstruction error in the whole image, under the constraint that Gg should give the (known) reconstruction error in Ω. Alternatively, the coefficients g can be optimized with the constraint expressed in the pixel domain, yielding equation (9),
The local constraint is propagated in all of the variables by the projection operator involved in the optimization process. Uniqueness theorems mentioned in §1.3 state that 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 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.
2.2.2. Details of the involved operations
In this part, more details are given about 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 image 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 of placing the coefficients g on a regular grid and convolving with the Gaussian kernel (x,y) . Thus, the operator G can be written G = where 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} = where U_{s}^{ T} is the sdownsampling operator and 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 , i.e. the Gaussian is truncated at . The resulting image is given by equation (10),
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 (10) is done on the convolution kernel support. If s < , 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 twodimensional 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 Fig. 3.
As the new image x is bigger than x_{0}, the sinogram Px and the acquired data d cannot be directly compared. The computed sinogram Px 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 of 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 with = and is cropped to the ROI.
2.2.3. Computational aspects
Using Gaussians as functions to iteratively express the error has several computational advantages. Gradientbased algorithms for solving (8) involve the computation of the forward operator PG 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 onedimensional Gaussian, as illustrated in Fig. 3.
The first advantage is of a 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 pointprojector and a pointbackprojector implementation, the pair can be exactly adjoint, besides giving more accurate results.
The second advantage is on the computational side. As the operator PG consists of projecting twodimensional Gaussians disposed on the image, it is equivalent to placing onepixel coefficients (Dirac functions in the continuous case) in the image on a grid denoted , projecting the image (with a pointprojector) and convolving the sinogram by a onedimensional Gaussian kernel. The same goes for the adjoint operator G^{ T} P^{ T} consisting of retrieving the Gaussians coefficients from a sinogram. The standard way to compute this operator would be backprojecting the sinogram, convolving by the twodimensional Gaussian (which is its own matched filter due to symmetry), and sampling the image on the grid to obtain 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 in the image domain. The resulting sampled sinogram is then backprojected with a pointbackprojector.
2.3. 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 §2.1, the second builds the resulting image as described in §2.2.
A complete implementation of the proposed method is available (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 pointprojector 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.
The location of the known zone Ω can be simply implemented as a tuple of pixels (i_{0},j_{0}) and a radius r for a circular zone.
In practice, the final image is cropped to the ROI. 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. Another approach, expressing the known zone constraint in the pixel domain [equation (9)], only requires Algorithm 2 to be run.
3. 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 ROI 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 synthetic sinogram is then truncated to the radius of the ROI. The truncated sinogram is the input of the methods. The proposed method is compared with 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 in the final image.
The first test involves the standard Shepp–Logan phantom (Fig. 4), 256 × 256 pixels. The ROI is embedded inside the `absorbing outer material' (ellipse with the highest gray values) to simulate a local tomography acquisition. For an easier interpretation of the line profiles in the final reconstructed images, the values of the standard phantom are multiplied by 250 so that all the values are between 0 and 500. The width of the extended for reconstruction image is N_{2} = 260.
This phantom is the `original' Shepp–Logan phantom, not the `modified Shepp–Logan phantom' where the contrast is improved. In our case, low contrast is important for the tests, as the cupping effect is stronger and directly visible in the reconstructed slice. For highcontrast images, the cupping effect only affects few lowfrequency components, and is thus less detrimental to the reconstruction quality.
Fig. 5 shows the result of the reconstruction with padded FBP and with the proposed method. The Gaussian coefficients were computed with = 4 on a grid of spacing s = 6. The known region radius is R = 20 pixels, and the extended image width is N_{2} = 260 pixels. On the padded FBP reconstruction, the cupping effect is relatively important due to the low contrast of the data.
By visual inspection, this method does not induce new artifacts in the reconstruction. Fig. 6 shows a line profile of this reconstruction. The cupping effect is visible for the padded FBP, and it has been removed using the proposed method. More importantly, the average reconstructed values are distributed around the true interior values. This provides an illustration of the uniqueness theorem: knowing the values of a subregion of the ROI enables the ROI to be exactly reconstructed (up to numerical errors). The reconstruction with the proposed method bears the same high frequencies as the FBP with full data, which is a good indication that this method does not induce new artifacts. The fact that the reconstruction has the same mean values as the true interior could enable quantitative analysis of the reconstructed volume, which is not easily achievable in local tomography.
Fig. 7 shows the difference between the reconstructions and the interior values (denoted x^{#}). As expected, the cupping effect is visible for padded FBP, while being almost entirely suppressed in the reconstruction with the proposed method.
The known zone constraint is important to remove the mean bias. As can be seen in Figs. 8 and 9, the cupping effect remains if no constraint is applied. In this case, the uniqueness theorem does not apply when there is no constraint, hence there are no guarantees that the method converges to an acceptable solution.
The second test involves the test image `Lena', 512 ×512 pixels, bearing both smooth regions and highfrequencies textures. Ellipses with strong intensity values have been superimposed in the exterior of the ROI to simulate absorbing material.
Fig. 10 shows the test setup. The known region has been 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.
Fig. 11 shows the difference between the true interior and the reconstruction with the proposed method, with and without the known zone constraint. In this case, the differences between the proposed reconstruction and the true interior are not significant enough to display the curves of absolute line profiles, so only the difference profiles are shown. The parameter = s = 3 has been used for these reconstructions. Fig. 12 shows the reconstructed images. As the contrast is high, the cupping effect is not prominent on the reconstructions; however, an excessive can be seen on the bottom side.
It is also interesting to visualize the reconstruction of the whole extended image. As can be seen in Fig. 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.
The third test involves the image of a pencil resulting from a scan at the ESRF ID19 beamline, 512 ×512 pixels, shown on Fig. 14. The width of the extended image is N_{2} = 520.
Fig. 15 shows profiles of the difference between the reconstruction and the true interior. On this image, a greater radius also improves the cupping removal. The profile of a line through the reconstructed image is depicted in Fig. 16.
As a last remark, Fig. 17 shows the result of this method without using the known zone constraint, that is, without applying the constraint = g_{0} in (8). As expected, there is a notnull mean bias, even if it has been reduced with respect to padded FBP.
Beside visual inspection, reconstructions can be quantitatively compared with the true interior of the test image. Table 1 shows a comparison with two image metrics: peak signaltonoise 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 yields better overall reconstruction quality than padded FBP. In particular, it does not induce spurious distortion in the reconstruction. For the `Lena' test case, a similar reconstruction quality was obtained with = (4, 6) with respect to = (3, 3), which indicates that a coarser basis does not always yield worse reconstruction results.
For all the reconstructions with the proposed methods, the conjugate gradient optimization algorithm was used. The convergence is reached within 400 iterations. The prototype method given by Paleo (2016) takes the following execution times on a machine with an Intel Xeon CPU E51607 3.00 GHz CPU, and a Nvidia GTX 750 Ti GPU: 7 s for a 260 ×260 extended slice, 32 s for a 520 ×520 extended slice, and 152 s for a 1040 ×1040 extended slice. A highperformance implementation of this method, where the data are kept on GPU and the projector is implemented as described in §2.2.3, can probably achieve much shorter execution times.
The proposed method depends on some parameters. The first is the size of the extended image, which should be big enough to model the contribution of the external part. The other parameters are the Gaussian standard deviation σ and the spacing s of the grid. Both are related in a way that the Gaussian functions should slightly overlap to approximate constant functions: if the value of s is high, then σ should also be high, and vice versa. These parameters essentially tune how coarse the Gaussian basis is: high values would yield fast convergence but coarse result, while small values would yield slow convergence and fine result.
Using a Gaussian basis does not yield an exact correction of the error, as Gaussian functions defined in equation (3) do not form a basis. For example, Gaussian functions do not yield a partition of unity, although a very close approximation of this property 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 seen in Figs. 6 and 16.
4. Advantages against the existing exact approaches
As mentioned in the Introduction, there are two main types of exact methods for local tomography. The first requires a special geometry: the FOV has to extend the object on at least one side (Clackdoyle & Defrise, 2010). This approach is not suited for interior tomography that we are focusing on. In the second type, a prior knowledge is used, either on a subregion or on the slice properties.
In the work of Kudo et al. (2008), an algorithm for exact local reconstruction is proposed, along with proof of uniqueness. This algorithm is based on differentiated backprojection and projection onto convex sets to invert the truncated Hilbert transform. Although this algorithm seems very promising, its implementation is not straightforward; for example, convex sets may be cumbersome to define, and projections on these sets are not easy to compute. To our knowledge, no standard implementation of this method can be found.
On the other hand, several compressedsensing methods were developed for interior tomography (Niinimäki et al., 2007; Chaves Brandao~dos Santos et al., 2014; Klann et al., 2015). As these methods boil down to solving an optimization problem, the known region constraint can be integrated in these cases.
To highlight the advantage of the coarse basis representation with respect to a standard pixel basis, the proposed method is tested against a method solving
for x in the pixel space, meaning that all the variables are free, except those in the known region. The term is the total variation (TV) seminorm which helps to reduce the illposedness of problem (11). This problem is a simple instance of the aforementioned compressed sensing methods, where the TV expresses the prior knowledge that the image is piecewiseconstant. Related work on this method is given in §3 of Rashed & Kudo (2010). The proposed methods have two fundamental differences from this approach: firstly, the variables are the coefficients of a coarse basis; secondly, the variables are optimized with respect to the reconstruction error.
Equation (11) was tested on the 256 ×256 Shepp–Logan phantom (Fig. 4) with the same conditions as previously. As may be seen from Figs. 18 and 19, it yields a nearperfect reconstruction (PSNR = 58.63, SSIM = 0.9999), as TV is adapted for this piecewiseconstant phantom. However, it took of the order of 5000 iterations of the preconditioned Chambolle–Pock optimization algorithm (Pock & Chambolle, 2011) to yield this result, as the convergence is extremely slow. An approximate result, with a slight cupping remaining, can be obtained with fewer iterations (Fig. 18b). This test case of 256 ×256 took approximately 67 s. The execution time for a 512 ×512 image, with 15000 iterations to converge, was measured to be 430 s. Thus, this method becomes impracticable, even with an efficient implementation, as it would take dozens of minutes for a single slice of modern datasets (2048 ×2048 or 4096 ×4096 pixels).
The advantage of the proposed method is to reduce the number of unknowns by modifying the forward model, expressing the variables in a coarse basis. On the implementation side, operators can be efficiently implemented according to §2.2.3, which is crucial when dealing with large image sizes. Results presented in this work were obtained with a prototype, which suggests that shorter execution times can be expected with an efficient implementation.
5. Conclusion
We have presented a new technique of local tomography reconstruction based on the knowledge of a zone of the ROI. This technique corrects the cupping effect in an initial reconstruction by expressing the error in a coarse basis of Gaussian functions. In accordance with local tomography uniqueness theorems, this method yields 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 a prototype of this method is given by Paleo (2016).
References
Aarle, W. van, Palenstijn, W. J., De Beenhouwer, J., Altantzis, T., Bals, S., Batenburg, K. J. & Sijbers, J. (2015). Ultramicroscopy, 157, 35–47. Web of Science PubMed Google Scholar
Bale, R. A., Grossman, J. P., Margrave, G. F. & Lamoureux, M. P. (2002). Multidimensional partitions of unity and Gaussian terrains. CREWES Research Report, Vol. 14. Google Scholar
Bilgot, A., Desbat, L. & Perrier, V. (2011). 2011 IEEE Nucl. Sci. Symp. Conf. Rec. pp. 40804085. CrossRef Google Scholar
Chaves Brandao dos Santos, L., Gouillart, E. & Talbot, H. (2014). 2014 IEEE International Conference on Image Processing (ICIP), pp. 1768–1772. Google Scholar
Clackdoyle, R. & Defrise, M. (2010). IEEE Signal Process. Mag. 27, 60–80. CrossRef Google Scholar
Clackdoyle, R., Noo, F., Junyu Guo, & Roberts, J. (2004). IEEE Trans. Nucl. Sci. 51, 2570–2578. Google Scholar
Clackdoyle, R., Roy, D., Defrise, M., Mennessier, C. & Mohamed, M.S. (2009). 2009 IEEE Nuclear Science Symposium Conference Record (NSS/MIC), pp. 3367–3371. IEEE. Google Scholar
Courdurier, M., Noo, F., Defrise, M. & Kudo, H. (2008). Inverse Probl. 24, 065001. CrossRef Google Scholar
Defrise, M., Noo, F., Clackdoyle, R. & Kudo, H. (2006). Inverse Probl. 22, 1037–1053. Web of Science CrossRef Google Scholar
Kannappan, P. & Sahoo, P. K. (1992). SIAM J. Math. Anal. 23, 1342–1351. CrossRef Google Scholar
Klann, E., Quinto, E. T. & Ramlau, R. (2015). Inverse Probl. 31, 025001. CrossRef Google Scholar
Kudo, H., Courdurier, M., Noo, F. & Defrise, M. (2008). Phys. Med. Biol. 53, 2207–2231. Web of Science CrossRef PubMed Google Scholar
Kyrieleis, A., Titarenko, V., Ibison, M., Connolley, T. & Withers, P. (2011). J. Microsc. 241, 69–82. CrossRef CAS Google Scholar
Lee, M., Han, Y., Ward, J. P., Unser, M. & Ye, J. C. (2015). SIAM J. Imaging Sci. 8, 2452–2486. CrossRef Google Scholar
Lewitt, R. M. (1992). Phys. Med. Biol. 37, 705–716. CrossRef CAS Google Scholar
Niinimäki, K., Siltanen, S. & Kolehmainen, V. (2007). Phys. Med. Biol. 52, 6663–6678. Google Scholar
Noo, F., Clackdoyle, R. & Pack, J. D. (2004). Phys. Med. Biol. 49, 3903–3923. Web of Science CrossRef PubMed Google Scholar
Noo, F., Defrise, M., Clackdoyle, R. & Kudo, H. (2002). Phys. Med. Biol. 47, 2525–2546. CrossRef Google Scholar
Palenstijn, W., Batenburg, K. & Sijbers, J. (2011). J. Struct. Biol. 176, 250–253. Web of Science CrossRef CAS PubMed Google Scholar
Paleo, P. (2016). Implementation of a new method for local tomography reconstruction, https://github.com/pierrepaleo/localtomo. Google Scholar
Pock, T. & Chambolle, A. (2011). International Conference on Computer Vision, pp. 1762–1769. IEEE. Google Scholar
Rashed, E. A. & Kudo, H. (2010). Mathematical Programming in the 21st Century: Algorithms and Modeling, 1676, 145–156. Google Scholar
RashidFarrokhi, F., Liu, K., Berenstein, C. A. & Walnut, D. (1997). IEEE Trans. Image Processing, 6, 1412–1430. CAS Google Scholar
Sastry, C. S. & Das, P. (2005). ANZIAM J. 46, 341–360. CrossRef Google Scholar
Shuangren Zhao, K. Y. & Yang, X. (2011). J. Xray Sci. Technol. 19, 155. Google Scholar
Tang, S., Yang, Y. & Tang, X. (2012). J. Xray Sci. Technol. 20, 405. Google Scholar
Van Gompel, G. (2009). Phd thesis, University of Antwerp, Antwerp, Belgiium. Google Scholar
Van Gompel, G., Defrise, M. & Van Dyck, D. (2006). Proc. SPIE, 6142, 61424B. CrossRef Google Scholar
Wang, G. & Yu, H. (2013). Phys. Med. Biol. 58, R161–R186. CrossRef Google Scholar
Yang, J., Yu, H., Jiang, M. & Wang, G. (2010). Inverse Probl. 26, 035013. Web of Science CrossRef Google Scholar
Ye, Y., Yu, H., Wei, Y. & Wang, G. (2007). Int. J. Biomed. Imaging, 2007, 63634. PubMed Google Scholar
Yu, H. & Wang, G. (2007). Int. J. Biomed. Imaging. 2007, 95295. Google Scholar
Yu, H. & Wang, G. (2009). Phys. Med. Biol. 54, 2791–2805. Web of Science CrossRef PubMed Google Scholar
Zeng, G. L. & Gullberg, G. T. (2000). IEEE Trans. Med. Imaging, 19, 548–555. CrossRef CAS Google Scholar
Zou, Y., Pan, X. & Sidky, E. Y. (2005). Phys. Med. Biol. 50, 13–27. CrossRef Google Scholar
This is an openaccess article distributed under the terms of the Creative Commons Attribution (CCBY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.