DIFFUSION RECONSTRUCTION FROM VERY NOISY TOMOGRAPHIC DATA

. As a consequence of very noisy tomographic data the reconstructed images are contaminated by severely ampliﬁed noise. Typically two remedies are considered. Firstly, the data are smoothed, this is called pre-whitening in the engineering literature. The disadvantage here is that the individually treated data sets could become inconsistent. Secondly, the image, reconstructed from the original data sets, is treated by methods of image smoothing. As example diﬀusion ﬁlters are mentioned. In this paper we present a method where the reconstruction of the smoothed image is performed in one step; i.e., we develop special reconstruction kernels, which directly compute the image smoothed by a diﬀusion ﬁlter. Examples from synthetic data are presented.


Introduction
In order to enhance the information content of reconstructed images feature extraction methods are applied.This paper is concerned with presenting a methodology to combine the reconstruction step and the feature extraction in one algorithm.One advantage of this approach is that the involved parameters are jointly optimised, or that, as is the case here for image smoothing, only one set of parameters is needed, resulting in superior reconstructions.This is achieved by precomputing reconstruction kernels which directly result in the enhanced images by computing scalar products of those kernels with the given data.To that aim one selects a suitable mollifier and solves a priori an auxiliary problem for the new reconstruction kernel.In that way, all operations, like differentiation for example, are performed on the analytically defined mollifier, which has to be sufficiently smooth.If it satisfies suitable invariance properties the overall costs for the calculation of the enhanced image are the same as for the original reconstruction step.
The calculation of linear functionals of the solution of inverse problems was mostly used for stabilizing the inversion method.First results are reported by Likht, [7].In [12] it was shown that many linear regularization methods like Tikhonov-Phillips, truncated singular value decomposition or Landweber iteration can be interpreted either as the application of a smoothing operator on the data followed by an inversion or first inversion followed by a smoothing of the reconstruction.Hence, the result can be written as linear functional of the solution (1) with a suitable mollifier e γ .Often it can be represented with the singular value decomposition of the operator.In the case of Tikhonov -Phillips regularization with the minimization of Af − g 2 + γ 2 f 2 the mollifier can be written as where the singular values fulfill Av n = σ n u n .The parameter γ plays the role of a regularization parameter.The method of the approximate inverse directly starts with the selection of such a mollifier, a reconstruction kernel ψ γ is precomputed as solution of A * ψ γ = e γ to represent the regularized solution as (3) which can be written as ( 4) with a smoothing operator M γ .For details the reader is referred to [12].These approaches are concerned with approximating a solution f of an equation Af = g.If special features of the solution are of interest we want to determine Lf where L is a suitable evaluation operator.Here we restrict our consideration to linear operators, both A and L. Instead of first calculating f as a regularized solution of Af = g and then applying the evaluation operator L we aim for determining Lf in just one step from the data g.Lambda tomography may serve as an example where A is the Radon transform and L is the Riesz potential I −1 defined via the Fourier transform as The motivation for this was to find local inversion formulas for the Radon transform and only later it was discovered that the result can be written that way, with the pseudo -differential operator I −1 , see [19,20] for the two -dimensional case.This approach was generalized to 3D X-ray CT, see [9,6] and to other imaging modalities, see [14,18].The combination of A a one-dimensional integral operator and L = d/dx was presented in [11].A general theory for the combination of two operators was developed in [13] and as application an edge detection problem for the Radon transform was analysed; i.e., A the Radon transform in two dimensions and L = ∂/∂x k for k = 1, 2. The combination of an arbitrary linear operator and L the calculation of wavelet coefficients, which means the calculation of wavelet components of f directly from the data was introduced in [10], an application in tomography can be found in [2,16].
The advantage of combining the two steps in one algorithm is, that the different parameters can be chosen optimally and that the overall computational costs is minimized.
In this paper we consider the case of tomography with extremely noisy data.Examples can easily be found, see for example electron microscopy, [5].As a consequence the reconstructions are severely disturbed and methods from image analysis have to be applied to get the searched-for information.When diffusion is used to smooth the reconstruction one considers the noisy image as the initial value of a heat equation.The regularization parameter is the time where the process is stopped.We derive here a reconstruction kernel which starting from the data directly leads to the smoothed reconstruction for a prescribed stopping time.No further regularization for the inversion of the Radon transform is needed.
Finally we can combine the diffusion and the reconstruction of derivatives to get the searched-for edge information.In that case much larger noise levels still can be handled.We conclude with numerical examples.

Approximate Inverse for Combining Reconstruction and Image Analysis
This section is based on [13] where we generalize the method of the approximate inverse as analyzed in [11].Here we present the main ideas, we assume that all the auxiliary problems are solvable.This imposes conditions on the mollifiers e γ , for details see [13].The mollifiers used in this paper fulfill those conditions.
Let A : X → Y be a linear operator between the Hilbert spaces X and Y and L : X → Z be a linear operator between the Hilbert spaces X and Z.As usual we first formulate the reconstruction part ( 5) Next an operation L on the so computed solution f for the image analysis is performed where A † denotes the generalized inverse of A. Now we adapt the concept of approximate inverse, first introduced in [8], where we now compute instead of Lf an approximation (Lf ) γ = Lf, e γ with a prescribed mollifier e γ .We formulate in the following theorem the principle of the reconstruction method.
Theorem 2.1.Let e γ be a suitably chosen mollifier such that the solution ψ L,γ of the auxiliary problem Then the smoothed version of the image analysis operation is directly computed from the given data g as Proof.We write the smoothed version of the image analysis part as Now we use the adjoint operator of L and the auxiliary problem to continue where in the last step we have used the original equation Af = g.
Definition 2.2.The operator S γ : Y → Z defined as is called the approximate inverse of A to compute an approximation of Lf and ψ L,γ is called the reconstruction kernel.
If we know the reconstruction kernel for computing f , then we can solve the above problem for computing Lf in the following way.
Theorem 2.3.Let ψγ be the sufficiently smooth reconstruction kernel for computing f , then the reconstruction kernel ψ L,γ for approximating Lf can be determined as (10) ψ L,γ = L ψγ where L acts on the first variable of ψγ .
Proof.The approximation of Lf is here computed as the application of L on f γ (x) = g, ψγ (x, •) .Interchanging the application of L and the integration, for sufficiently smooth ψγ , gives the result.
If the original mollifier e γ is not smooth enough, then an additional smoothing in the determination of L ψγ may be needed.
It is shown in [13] that S γ is a regularization for computing Lf if the smoothness of e γ is adapted to the smoothing of A and the inverse of L in the following sense (11) lim The computational efficiency of the approximate inverse heavily depends on the use of invariances.We consider again the reconstruction problem in tomography.If we chose for each reconstruction point x a special mollifier, namely e γ (x, •), then the reconstruction kernel also depends on x, the number of values to store is then the number of reconstruction points times the number of data.If we use invariances, for example translation and rotational invariances of the Radon transform and we use these invariances to produce the mollifier we can reduce this number of values to compute and store to just the number of views per direction.The mathematical basis for this can be found in [11].

Diffusion Reconstruction of the Solution in Tomography
The mathematical model of computerized tomography in two dimensions, for the parallel geometry, is the Radon transform, see e.g.[15].It is defined as where θ ∈ S 1 is a unit vector and s ∈ R. In the following we summarize a few results.The central slice theorem, or projection theorem, is nothing but the formal application of the adjoint operator for fixed direction θ on exp(ısσ) (12) Rf (θ, σ) = (2π) 1/2 f (σθ) .
The inversion formula for the two -dimensional Radon transform is Inverse Problems and Imaging Volume X, No. X (200X), X-XX where R * is the backprojection g(θ, x, θ )dθ and the Riesz potential I −1 is defined with the Fourier transform where the Fourier transform acts on the second variable.
To reduce the unavoidable noise in the data typically filters, like Shepp -Logan filter, are applied.For very large noise level this is not sufficient to get the desired information as Figure 1 (left) shows.
For the simulations we used the standard Shepp-Logan phantom with the original densities; i.e., 2 for the bones, and 1 for the brain.The number of data was 720 directions and 1305 rays per view.To the data 12% noise was added.It is clear that then the 'three little tumors' with density differences of 1% are lost in the data.It is known, that this evolution process for the linear heat equation smears out edges to some extent, but nevertheless it gives acceptable results.Nonlinear diffusion would produce better results, see e.g.[17,21], but when large series of images are considered, this may be too time-consuming.
A pre-whitening with a diffusion is in principle possible, but when the stopping of the smoothing is performed adaptively, which means possibly at different times for the projection, then the data become inconsistent which leads to additional artifacts in the reconstructions.
In the following we want to combine as inverse problem with the operator A the Radon transform, A = R, and as evaluation operator L T the time evolution for the heat equation for fixed time T , This operator maps discontinuous functions into infinitely often differentiable functions.So we observe that, due to the injectivity of the Radon transform and the smoothing property of the operator L, the combined problem is well-posed, so no additional regularization is needed.Consequently we choose as mollifier the Delta distribution e γ (x, •) = δ x .The task is to compute a reconstruction kernel which directly results in LR −1 g.We represent the solution of the heat equation with the Green's function as (14) u(T, The combined reconstruction kernel then has to fulfill for each x the equation in y R * ψ T (x, y) = L * T δ x (y) = G(T, x − y) .We solve this auxiliary equation with the following well-known approach, by writing G as R −1 RG getting leading to ψ T = (4π) −1 I −1 RG(T, •) .Using the projection theorem we compute this with the inverse Fourier transform and Ĝ(T, ξ) and integration by parts we find with Formula(7.4.7) from [1] (16) where D is the so-called Dawson integral, see [1], page 295 and 319  The computing time for the combined reconstruction was the same as the reconstruction of the density itself.

Smooth Reconstruction of Derivatives in Tomography from Noisy Data
In this section we aim to calculate derivatives when the data are very noisy.This is part of edge detection algorithms, see e.g.[4].Classically we would compute smoothed derivatives of the smoothed images as shown for example in Figure 1 on the right.In the result with 12% noise, the information is almost completely lost.Hence we repeat here the calculation with half of the noise; i.e., 6% additive noise, see Figure 3.The artifacts outside the object can easily be removed by implementing the support theorem of Boman and Quinto, see [3], which states that the object vanishes on lines parallel to θ not meeting the support of the data.see e.g.[15], and generalizations for higher derivatives.
In combination with the diffusion the mollifier is smooth, hence we also can apply Theorem 2.3.With T ) The implementation of this formula leads both to fast algorithms and also to much better reconstructions of the derivative.In Figure 4 we present the results for the same noise level as above, and also twice the noise level, as in the last section.

Inverse Problems and Imaging
Volume X, No. X (200X), X-XX FIGURE 4: Sum of absolute values of the diffusion reconstruction of the derivatives with respect to the coordinates for 6% noise (left) and 12% noise (right).The image on the left hand side in Figure 4 was produced with the same data as the one in Figure 3.It is obvious that the newly proposed method can handle even larger noise, as shown in Figure4 (right).Hence the optimal combination of the different procedures, here diffusion filter and edge detection, in just one algorithm gives much better results.

FIGURE 1 :
FIGURE 1: Standard Shepp-Logan filter applied to very noisy data ( left ) and after applying a diffusion process( right ) Diffusion is known to smooth those images.The reconstructed image f γ is considered as initial value u(0, •) = f γ of a heat equation u t = ∆u .The time evolution is stopped at a suitable time T to get a smooth image, see Figure 1(right).Here and in the following calculations we used T = 0.1.It is known, that this evolution process for the linear heat equation smears out edges to some extent, but nevertheless it gives acceptable results.Nonlinear diffusion would produce better results, see e.g.[17,21], but when large series of images are considered, this may be too time-consuming.A pre-whitening with a diffusion is in principle possible, but when the stopping of the smoothing is performed adaptively, which means possibly at different times for

( 17 )
D(s) = exp(−s 2 ) s 0 exp(t 2 )dt Inverse Problems and Imaging Volume X, No. X (200X), X-XXAlgorithms are available for the efficient and accurate numerical approximation.The application of this reconstruction kernel directly leads to reconstructions where the noise is smoothed out.

FIGURE 2 :
FIGURE 2: Left: reconstruction and diffusion separately, same as Fig1(right) and on the right the combined reconstruction and diffusion.

FIGURE 3 :
FIGURE 3: Sum of absolute values of the central difference quotient with respect to the coordinates for the image where reconstruction and diffusion are calculated separately (left) and reconstruction with the reconstruction kernel introduced here (right ).Now we change the reconstruction kernel such that we directly calculate the derivatives.Either we use the following result.The Radon transform of a derivative is (18) R ∂ ∂x k f (θ, s) = θ k ∂ ∂s Rf (θ, s) ψ T ( x, θ − s) = θ k ψ T ( x, θ − s)and the derivative of Dawson's integral as (20) D (s) = 1 − 2sD(s) we find the reconstruction kernel ψ T,k for approximating a smoothed version of ∂f /∂x k