Abstract

Phase contrast tomography is a high sensitivity medical imaging technique. Several regularization methods have been used in the literature to obtain stable solutions for the phase retrieval or the phase contrast tomography problems. Yet, the functional framework and the convergence properties of the methods have not been studied in detail. In this work, the convergence properties of regularization approaches for phase retrieval and phase contrast tomography are investigated.

1. Introduction

X-ray microtomography is a very useful imaging technique in material science [1] and medical imaging [26]. Yet its sensitivity is improved for soft tissues within dense materials with phase contrast imaging obtained with third-generation synchrotrons. X-ray in-line phase tomography is based on coupling of tomography and phase retrieval to reconstruct the real part of the refractive index [7, 8]. For coherent X-rays, phase contrast can be achieved by letting the beam propagate in the free space after interaction with the object and by recording the intensity for one or several propagation distances and for several projection angles (Figure 1). The relationship between the phase shift induced by a sample and the intensity recorded at a given sample-to-detector distance relies on the Fresnel diffraction theory [912]. Several linear algorithms have been proposed for the phase retrieval from the Fresnel diffraction patterns [1, 10, 1318] valid under some restrictive assumptions. Yet, in these works, the convergence of the regularization methods used to obtain stable solutions has not been investigated. The main purpose of this paper is to detail the functional properties of the direct and inverse operators involved in these problems on the basis of classical results of harmonic analysis. Our aim is also to study the convergence of the regularized solutions when the noise and the regularization parameter tend to zero for the phase retrieval and the phase contrast tomography problems in the good functional framework. To prove stability with respect to the noise, we have to define a parameter choice rule for depending on the noise level . This procedure is typical for regularizing ill-posed problems [19] and is presented in detail.

This paper is organized as follows. In Section 2, we present in detail the phase retrieval and phase contrast tomography direct problems. In Section 3, we consider the various regularization methods proposed in the literature for the linearized phase retrieval problem. We investigate the behaviour of the solutions as the regularization parameter tends to zero with the noise level. In Section 4, we study the convergence properties of the solutions of the phase contrast tomography problem obtained with various regularization methods. These methods are based on a combination of inversion schemes for the Radon projector and the phase operator.

2. Phase Retrieval and Phase Contrast Tomography

2.1. The Direct Problem for Phase Contrast Tomography

In the following, we will consider a monochromatic, coherent, parallel X-ray beam. The real and imaginary parts of the complex refractive index to be reconstructed, denoted as and , are defined on a 3D bounded domain () with spatial coordinates . For the sake of simplicity, we assume that is the Cartesian product of bounded intervals on the coordinate axes. We denote by the rotated spatial coordinate system for an angle around the -axis (Figure 1).

The optical properties of an object interacting with coherent X-rays of wavelength are related to its complex refractive index given by where is the refractive index decrement and is the absorption index [20]. For a fixed projection angle , thin objects, and straight line propagation of the beam along the direction, this interaction can be described by a transmittance function of the coordinates : where is the absorption and is the phase shift induced by the object for the projection angle [8]. The phase and the negative logarithm of the absorption, , are the projections of the absorption index and refraction index, respectively: where is the propagation direction of the X-rays. In the framework of the Fresnel diffraction theory, the intensity detected at a distance after the sample is given by the squared modulus of the 2D convolution between the transmittance and the Fresnel propagator for a distance downstream of the object [8, 11, 21]: where

The intensity operator can thus be considered as a nonlinear function of .

2.2. Linear Approaches of the Phase Retrieval Problem

For each value of the projection , the phase retrieval problem to recover the phase from the intensity patterns set in the plane has been linearized in the framework of several models. We first review some widely used relations between the phase and the measured intensity. The Fourier coordinates in the plane will be denoted as . For sufficiently small propagation distance , the function varies in an approximately linear way in and the phase can be written: (see [2224]). In the small limit, the transport of intensity equation (TIE) leads to the linear approximating equation: where is a Laplacian operator, is the intensity at a finite object-detector distance , is the intensity contrast, is the wave number of the incident monochromatic X-rays, and is the phase shift in the exit plane of the object.

In the framework of the contrast transfer function (CTF) model [1418, 25], in the small-relative-phase limit corresponding to the assumption, the relation between the phase and the data can be written in the Fourier space with a linear operator : with and where the 2-dimensional Fourier transform and its inverse are given by This approach is more accurate than the TIE method which relies on a short propagation distance assumption. In the following, we consider only the singularity at the origin because the other singularities can be removed with a combination of several propagation distances.

2.3. The Coupling of Phase Retrieval with the Radon Projection Operator

The direct intensity operator can be rewritten with the Radon projection operator. We first summarize some properties of this projection operator [26]. Let be a bounded open domain, the mathematical model for 2D-tomography is the Radon transform which maps a function to its line integrals.

Definition 1. Let be the line defined by , with and , and let be the range of the Radon transform and ; the Radon transform for is defined by

Let and let ; the adjoint of is given by The proof can be found in [26].

For parallel beam projection, with a beam parallel to the plane, and the line for the coordinate : where denotes a tensor product and is the identity on the variable. The Sobolev spaces are very useful to characterize the smoothing properties of the Radon transform. The Sobolev norms for functions in are defined as Let ; the one-dimensional Fourier is defined by . Let , the Sobolev space on denoted as , be the set of all distributions on such that the Fourier transform is a function and such that the Sobolev norm of , , is finite [26]. The norm is defined as

The Radon operator is a smoothing operator as detailed in the following proposition.

Proposition 2. Let be a distribution with a compact support; then, the following bound holds for :

The proof of this proposition can be found in [26]. It is now possible to use the mapping properties of the Radon transform in the context of phase contrast tomography and in a Hilbert spaces framework.

Proposition 3. Let be the range of the Radon operator for the height z; the operator is continuous.

Proof. For , is continuously embedded in and thus for compactly supported functions, there is a positive constant such that Using (18), there exists a positive constant such that
In a Hilbert spaces context, the Radon transform is thus a linear continuous operator from to . For parallel beam projection, let be the range of the Radon operator for the height . For each value of , the Radon operator is continuous from to . By integration on the variable , the operator is continuous.

With the Radon projection operator, (4) can be rewritten:

3. Regularization for Phase Retrieval

In this section, we investigate the convergence properties of some regularization approaches for phase retrieval. In a generic way, we denote by the unknown phase to reconstruct, by the intensity data, and by the direct operator.

3.1. Preliminaries

Let us assume that the phase is defined on a bounded domain . It is well known that the Laplacian operator of the TIE approach is continuous [27]. The phase operators involved in the direct problem can be understood as Fourier multiplication operators, with a Fourier multiplier and written as The Fourier multiplier of the CTF method is given by It is bounded and thus the direct operator in the CTF method is bounded from to with the Parseval theorem. Some functional properties are lost with the linearization of the direct forward problem. The analysis of the nonlinear direct operator shows that it is a compact operator [28] because of the convolution with the Fresnel propagator. These smoothing properties are lost in the linearized methods presented above. With the linearization, the Laplacian operator in the TIE version or the CTF direct operator is not compact operator.

The inverse Fourier multipliers involved in the TIE or CTF approaches and define tempered distributions with singular kernels. It is necessary to use the finite part to give a sense to these distributions.

Definition 4. Given a function and given a function , the finite part of the integral is defined as

The following proposition is useful to show that some finite part distributions are well defined.

Proposition 5. Let one assume that the function satisfies the size condition then, it is possible to define the finite part of [2931].

Thus, and are well-defined tempered distributions. For data , the product defines a tempered distribution in [29, 30]. The inverse Fourier transform of is well defined. The data function has a bounded support, , and the convolution of and is well defined [2931].

Distribution solutions are obtained because of the singularity of the Fourier multiplier at the origin. The inverse problem formulated with CTF or TIE approach is not a classic linear problem with a compact operator between Hilbert spaces and . The inverse is defined in a distributional sense. A regularized distribution is obtained when the zero frequency component of the test function is suppressed. It is the idea of the “quasiparticule” method detailed below. More generally, it is possible to construct an approximate inverse with a bounded kernel in . Let be the approximate inverse operator, depending on the regularization parameter . This approximate inverse must verify the pointwise convergence, for all : Let be the solution obtained from the noisy data , . We assume that the noisy data and the unnoisy data are such that where is the noise level and is the norm. We are interested in this work in the convergence of various regularization methods. We will make use of the following inequality for :

3.2. Regularization for CTF

In this section, we investigate the convergence properties of some regularization methods for CTF. Several methods of regularization have been tested [1418] to solve this linear inverse problem, like the quadratic Tikhonov regularization used for the classical inverse problems or the “quasiparticule approach.” We detail here the convergence properties for the Tikhonov and quasiparticule approach and for the approximate inverse obtained with mollification.

3.2.1. Tikhonov Regularization for CTF

The Tikhonov regularization method is a convergent regularization method.

Proposition 6. The Tikhonov regularization method is a convergent regularization method with an optimal convergence rate .

Proof. For the Tikhonov regularization, the approximate inverse can be written: where is a Fourier multiplication operator given in (10). The approximate inverse for CTF is also a Fourier multiplication operator: With , the Fourier multiplier is given by
We first obtain an upper bound for the norm of the operator : where the first equality results from Plancherel’s theorem.
With the identity we obtain
For , we have With the change of variable, and with the identity , for , we obtain From the Rieman-Lebesgue convergence theorem, the Fourier transform is a continuous function. Thus, there is a constant positive such that Reporting in (27), we obtain for the Tikhonov regularization and thus the optimal convergence rate is obtained for the a priori choice, , and the convergence rate is

3.2.2. Quasiparticle Regularization

In the framework of the quasiparticule approach [2224], the intensity contrast is filtered in the Fourier space and is replaced by : where denotes the Heaviside step function and is a threshold for this filter.

Considering only the singularity of the multiplier at the origin, this threshold determines a threshold such that Let be the filter function defined by The regularized inverse can be written:

With the identity , for , we obtain that there is a constant such that

Moreover, we have By the Lebesgue dominated convergence theorem, this term tends to zero as . Yet it is impossible to derive any convergence rate of the method as without any assumption on the position of the spectra of with respect to the threshold. In the more favorable case, if the spectra of the function are above , this term is zero and The convergence is obtained for a choice of the regularization parameter such that as .

3.2.3. Regularization with Mollification

It is not possible to derive any convergence rate for the “quasiparticule” regularization without any assumption about the spectrum of the function to regularize. In this section, we detail the regularization with mollification which is valid without any restrictive assumption.

We consider a mollifier such that and we assume that there exists such that Typical examples of mollifier are the Gaussian function or where is a normalizing constant.

We consider as approximate inverse the operator defined for by

Proposition 7. The CTF formulation of the phase retrieval inverse problem can be regularized by the mollification regularization with an exponent . The optimal convergence rate is .

Proof. Let us assume ; then,
From the Rieman-Lebesque theorem, is continuous and decreases at infinity and thus is bounded for the higher value of .
As , for , there is a constant such that
Thus, we obtain with the definition of approximate inverse The mollifier is defined on and thus we obtain With (27), the optimal convergence rate obtained with the choice is

4. Regularization for Phase Contrast Tomography

In this section, we present some regularization methods for phase contrast tomography which are combination of regularization methods for tomography and phase retrieval. We will use the mollifier regularization for the phase retrieval problem (Section 3.2.3).

4.1. Inversion Formula for Tomography and Regularized Inversion Approaches

We first summarize the well-known filtered backprojection inversion formula for a 2D geometry.

Definition 8. The Riesz potential for is defined by

The Riesz potential is an isomorphism of Hilbert spaces. For , it is a compact operator. The Riesz potential used for tomography reconstruction, , for , is the operator with Fourier multiplier : where the partial Fourier transform of in the variable is defined by the formula

The inversion formula is summarized by the following theorem [26].

Theorem 9. Let , the set of infinitely differentiable functions with 2D compact support; then, can be reconstructed with the formula where is the adjoint of the Radon transform.

This reconstruction formula can be generalized for parallel beam projection with tensor product of identity operators on the coordinate. In this case, is the spatial frequency associated to the axis (Figure 1). For parallel beam projection, the reconstruction formula is

Several regularized inverse operators have been studied in the literature, based on truncated filtered backprojection [26, 32] or on the approximate inverse obtained with mollifiers [33].

In the framework of the truncated filtered backprojection, the approximate inversion formula can be written as , where is the inverse Fourier transform of the truncated multiplier and denotes the convolution in the variable: where is a cutoff frequency. In the following, we denote by the regularized Riesz potential obtained with this kernel, with the power law , where is a positive constant.

In the framework of the approximate inverse theory obtained with mollifiers, the approximate inverse is given by where is the reconstruction kernel and is a mollifier approximating the distribution, centered about 0 and with mean value 1. The reconstruction kernel and the mollifier are related by A family of reconstruction kernel is obtained for as

The theory of the approximate inverse Hilbert space is fully presented in [33, 34]. This method of reconstruction allows for error estimates. Given measured data , assuming that the range of the Radon transform is , the reconstruction formula can be written, for : For a 3D geometry, an approximate inverse of the Radon projector can be obtained with the tensorial product .

4.2. The Approximate Inverse for Phase Contrast Tomography

The phase contrast tomography problem can be written in a generic way as where denotes in this section the unknown imaginary index part and is the intensity data. The operator , , is linear and continuous. The same result is true for the operator .

The approximate inverse of the direct operator can be expressed as where and are approximate inverses for the Radon projection operator and the phase retrieval operator . We will study successively the convergence of two inversion schemes: a coupling of the CTF approach regularized with mollifiers with the tomography inverse obtained with truncated filtered backprojection or with mollifiers.

4.2.1. Coupling of CTF and Riesz Potential Regularized with a Rectangular Window

We first study the case of a coupling of CTF regularized with mollifiers with a Riesz potential regularized with a rectangular window. In this case, the approximate inverse is given by In this formula, is the Fourier multiplier of the CTF approach regularized with the mollifier approach and and are the Fourier coordinates associated with the and axes. We obtain a formula similar to the one in Gureyev et al.’s paper [25]. We demonstrate in the following that the inversion obtained with this method is convergent.

Proposition 10. Let us assume that, for each z, . The coupling of mollifier regularization for CTF with an inversion method for tomography based on a Riesz potential regularized with a rectangular window is a convergent regularization method. The method is convergent for . The optimal convergence rate is obtained as for and for .

Proof. In order to use a single regularization parameter, we assume that the threshold for the Riesz operator is given by , where is a positive constant. We first obtain a bound for the norm of the truncated Riesz potential: The norm of a tensorial product of operators is the product of the norm of the operators. With the upper bound for the norm of the approximate inverse obtained with the mollifier regularization for the CTF operator (54) and the fact that the operator is bounded, we get Ignoring tensorial products with identity operators for the sake of simplicity, we have the equality we obtain Let us assume that , with a constant ; then, , and (Proposition 2): where is a positive constant. The same bound can be obtained for 3D geometry with parallel beam projection. With the bound obtained for for mollifier regularization, we have
We thus obtain the following bound: The method is convergent for . The optimal convergence rate is obtained as for and the choice and for and the choice .

4.2.2. Coupling of CTF and Mollified Approximate Inverse

In this section, we study the coupling of mollifier regularization for CTF with the inverse Radon transform obtained with the approximate inverse method based on mollifier kernels [33]. Classical results for mollifier show that converge to in as . To have a more precise bound, we will use the following lemma about difference quotients in Sobolev spaces [27].

Definition 11. Given a unit vector and a measurable function , the -translate of in the direction is denoted as and defined by For a measurable function , the difference quotient of in the coordinate direction of length is the measurable function defined by

Proposition 12. Let and let be the weak derivative of u. Then, for any , with and

The convergence results are summarized by the following proposition.

Proposition 13. Let us assume that can be extended to a larger domain such that . The coupling of the mollifier regularization for CTF and of the regularization of the inverse Radon transform with the approximate inverse approach based on mollifiers is a convergent method with a convergence rate .

Proof. We first bound the norm of . We assume that the range of the Radon transform is , with the Young and Schwartz inequalities we obtain, for :
The Fourier slice theorem can be written [32]: where denotes the Fourier transform with respect to the variable and is the 2D Fourier transform; thus, with a polar change of variable. By integration on the variable , there is a positive constant such that We first bound the norm :
With the Hölder inequality, we obtain
By integration, Let us note , the diameter of the domain .
Let us assume that can be extended to a larger bounded domain , such that with . Then, with proposition (8) with , we obtain for a positive constant .
The regularized inverse can be written as . Using (86) and (54) for the mollifiers regularization of CTF, we thus obtain the bound for a positive constant . From the identity and with (90) and (57), we obtain for
With (27), we obtain The optimal convergence rate obtained is with the choice .

It should be noted that the convergence rate obtained is very low. Yet it can not be improved with a Tikhonov regularization for CTF: with the bounds for the norm of for the Tikhonov method (36) and for the norm of   (86), (92) shows that the convergence of towards zero is not achieved.

5. Conclusion

In this paper, we have investigated the convergence of regularization methods for phase retrieval and phase contrast tomography. For the phase retrieval problem, the Tikhonov regularization, quasiparticule approach, and regularization with mollification methods have been studied and the optimal convergence rates have been estimated. The phase contrast tomography problem has also been studied. The inverse of the Radon operator is regularized with a rectangular window for the Riesz kernel or with the approximate inverse method based on mollifiers. The a priori convergence rate for the combined inversion of the Radon and phase retrieval operators is estimated for each approach.

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.