# Regularization Methods for Phase Retrieval and Phase Contrast Tomography

Authors: Bruno Sixou
Journal: Abstract and Applied Analysis
 (2015)

---

## 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.


---

## Body

## 1. Introduction

X-ray microtomography is a very useful imaging technique in material science [1] and medical imaging [2–6]. 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 [9–12]. Several linear algorithms have been proposed for the phase retrieval from the Fresnel diffraction patterns [1, 10, 13–18] 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.Figure 1
Experimental set-up in propagation based phase contrast tomography from a single propagation based phase contrast image per projection showing the coordinate system.This paper is organized as follows. In Section2, 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δn and β, are defined on a 3D bounded domain (Σ) with spatial coordinates (x,y,z). For the sake of simplicity, we assume that Σ is the Cartesian product Σx×Σy×Σz of bounded intervals on the coordinate axes. We denote by (xθ,yθ,z) the rotated spatial coordinate system for an angle θ around the z-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
(1)nx,y,z=1-δnx,y,z+iβx,y,z,
where δn 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 yθ direction, this interaction can be described by a transmittance function T of the coordinates Xθ=(xθ,z):
(2)TθδnXθ=exp⁡-BθβXθ+iφθδnXθ=aθ(Xθ)exp⁡(iφθ[δn](Xθ)),
where aθ[β](Xθ) is the absorption and φθ[δn](Xθ) is the phase shift induced by the object for the projection angle θ [8]. The phase and the negative logarithm of the absorption, Bθ[β](Xθ)=-log⁡(aθ[β](Xθ)), are the projections of the absorption index and refraction index, respectively:
(3)BθβXθ=2πλ∫‍βyθ,Xθdyθ,(4)φθδnXθ=2πλ∫‍1-δnyθ,Xθdyθ,
where yθ is the propagation direction of the X-rays. In the framework of the Fresnel diffraction theory, the intensity detected at a distance d after the sample is given by the squared modulus of the 2D convolution between the transmittance Tθ and the Fresnel propagator Pd for a distance d downstream of the object [8, 11, 21]:
(5)Id,θδn=Tθδn*Pd2,
where
(6)Pd(Xθ)=1iλdexp⁡iπλdXθ2.The intensityId 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 Xθ=(xθ,z) 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 Xθ=(xθ,z) plane will be denoted as k→=(ξ,η). For sufficiently small propagation distance d, the function φ varies in an approximately linear way in d and the phase can be written:
(7)φXθ=φd=0Xθ+φ1Xθd
(see [22–24]). In the small d limit, the transport of intensity equation (TIE) leads to the linear approximating equation:
(8)gdXθ=-dkΔφd=0Xθ,
where Δ is a Laplacian operator, Id is the intensity at a finite object-detector distance d, gd(Xθ)=(IdXθ/Id=0)-1 is the intensity contrast, k=2π/λ is the wave number of the incident monochromatic X-rays, and φd=0 is the phase shift in the exit plane of the object.In the framework of the contrast transfer function (CTF) model [14–18, 25], in the small-relative-phase limit corresponding to the assumption,
(9)φd=0X→θ-πdk→k-φd=0X→θ+πdk→k≪1,
the relation between the phase and the data gd can be written in the Fourier space with a linear operator A:
(10)gd=F-12sinuk→2Fφk→=Aφ,
with u=πdλ and where the 2-dimensional Fourier transform and its inverse are given by
(11)f^k→=Fgk→=12π∫x→∈R2‍e-ix→·k→fx→dx→,F-1gx→=12π∫k→∈R2‍eix→.k→fk→dk→.
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 Ω⊂R2 be a bounded open domain, the mathematical model for 2D-tomography is the Radon transform R which maps a function f∈L1(Ω) to its line integrals.Definition 1.
LetL(θ,xθ) be the line defined by L(θ,xθ)={yθθ¯*+xθθ¯:yθ∈R}, with θ¯=(cos⁡(θ),sin(θ)) and θ¯*=(-sin(θ),cos⁡(θ)), and let Ra be the range of the Radon transform and Z=[0,π]×Ra; the Radon transform for f∈L1(Ω) is defined by
(12)Rfθ,xθ=Rθfxθ=∫t∈Lθ,xθ∩Ω‍ftdt.Letg∈L1([0,π]×R) and let U¯=(x,y); the adjoint R*:L2(Z)→L2(Ω) of R is given by
(13)R*gU¯=∫0π‍gθ,U¯·θ¯dθ.
The proof can be found in [26].For parallel beam projection, with a beam parallel to theX=(x,y) plane, f∈L1(Σ) and L(θ,xθ,z) the L(θ,xθ) line for the coordinate z:
(14)Iz⊗Rfθ,xθ,z=Rθfxθ=∫t∈Lθ,xθ,z∩Σ‍ftdt,
where ⊗ denotes a tensor product and Iz is the identity on the z variable. The Sobolev spaces Hα are very useful to characterize the smoothing properties of the Radon transform. The Sobolev norms ·Hα for functions in Rn are defined as
(15)fHα2=∫Rn‍1+ξ2αFfξ2dξ.
Let g∈L1([0,π]×R); the one-dimensional Fourier Fp(g)(θ,τ) is defined by Fp(g)(θ,τ)=(1/2π)∫R‍exp⁡-iτpg(θ,p)dp. Let s∈R, the Sobolev space on [0,π]×R denoted as H0,s([0,π]×R), be the set of all distributions g on ([0,π]×R) such that the Fourier transform Fp(g) is a function and such that the Sobolev norm of g, gH0,s(Z)2, is finite [26]. The norm gH0,s(Z)2 is defined as
(16)gH0,s(Z)=∫0π‍∫R‍1+τ2sFpgθ,τ2dτdθ.The Radon operator is a smoothing operator as detailed in the following proposition.Proposition 2.
Letf be a distribution with a compact support; then, the following bound holds for s≥0:
(17)R(f)H0,s+1/2(Z)≤CfHsR2.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.
LetRaz be the range of the Radon operator for the height z; the operator φ=Iz⊗R[δn]:L2(Σ)→L2([0,π]×Σz×∪Raz) is continuous.Proof.
Fors>0, H0s(R2) is continuously embedded in L2(R2) and thus for compactly supported functions, there is a positive constant C1 such that
(18)fHs(R2)≤C1fL2R2.
Using (18), there exists a positive constant C such that
(19)R(f)H0,0(Z)≤R(f)H0,1/2(Z)≤CfL2R2.
In a Hilbert spaces context, the Radon transform is thus a linear continuous operator fromL2(R2) to L2(Z). For parallel beam projection, let Raz be the range of the Radon operator for the height z. For each value of z, the Radon operator is continuous from L2(Σx×Σy) to L2([0,π]×Raz). By integration on the variable z, the operator φ=Iz⊗R[δn]:L2(Σ)→L2([0,π]×Σz×∪Raz) is continuous.With the Radon projection operator, (4) can be rewritten:
(20)φθδnXθ=2πλIz⊗R1-δnθ,Xθ.


## 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δn and β, are defined on a 3D bounded domain (Σ) with spatial coordinates (x,y,z). For the sake of simplicity, we assume that Σ is the Cartesian product Σx×Σy×Σz of bounded intervals on the coordinate axes. We denote by (xθ,yθ,z) the rotated spatial coordinate system for an angle θ around the z-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
(1)nx,y,z=1-δnx,y,z+iβx,y,z,
where δn 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 yθ direction, this interaction can be described by a transmittance function T of the coordinates Xθ=(xθ,z):
(2)TθδnXθ=exp⁡-BθβXθ+iφθδnXθ=aθ(Xθ)exp⁡(iφθ[δn](Xθ)),
where aθ[β](Xθ) is the absorption and φθ[δn](Xθ) is the phase shift induced by the object for the projection angle θ [8]. The phase and the negative logarithm of the absorption, Bθ[β](Xθ)=-log⁡(aθ[β](Xθ)), are the projections of the absorption index and refraction index, respectively:
(3)BθβXθ=2πλ∫‍βyθ,Xθdyθ,(4)φθδnXθ=2πλ∫‍1-δnyθ,Xθdyθ,
where yθ is the propagation direction of the X-rays. In the framework of the Fresnel diffraction theory, the intensity detected at a distance d after the sample is given by the squared modulus of the 2D convolution between the transmittance Tθ and the Fresnel propagator Pd for a distance d downstream of the object [8, 11, 21]:
(5)Id,θδn=Tθδn*Pd2,
where
(6)Pd(Xθ)=1iλdexp⁡iπλdXθ2.The intensityId 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 Xθ=(xθ,z) 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 Xθ=(xθ,z) plane will be denoted as k→=(ξ,η). For sufficiently small propagation distance d, the function φ varies in an approximately linear way in d and the phase can be written:
(7)φXθ=φd=0Xθ+φ1Xθd
(see [22–24]). In the small d limit, the transport of intensity equation (TIE) leads to the linear approximating equation:
(8)gdXθ=-dkΔφd=0Xθ,
where Δ is a Laplacian operator, Id is the intensity at a finite object-detector distance d, gd(Xθ)=(IdXθ/Id=0)-1 is the intensity contrast, k=2π/λ is the wave number of the incident monochromatic X-rays, and φd=0 is the phase shift in the exit plane of the object.In the framework of the contrast transfer function (CTF) model [14–18, 25], in the small-relative-phase limit corresponding to the assumption,
(9)φd=0X→θ-πdk→k-φd=0X→θ+πdk→k≪1,
the relation between the phase and the data gd can be written in the Fourier space with a linear operator A:
(10)gd=F-12sinuk→2Fφk→=Aφ,
with u=πdλ and where the 2-dimensional Fourier transform and its inverse are given by
(11)f^k→=Fgk→=12π∫x→∈R2‍e-ix→·k→fx→dx→,F-1gx→=12π∫k→∈R2‍eix→.k→fk→dk→.
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 Ω⊂R2 be a bounded open domain, the mathematical model for 2D-tomography is the Radon transform R which maps a function f∈L1(Ω) to its line integrals.Definition 1.
LetL(θ,xθ) be the line defined by L(θ,xθ)={yθθ¯*+xθθ¯:yθ∈R}, with θ¯=(cos⁡(θ),sin(θ)) and θ¯*=(-sin(θ),cos⁡(θ)), and let Ra be the range of the Radon transform and Z=[0,π]×Ra; the Radon transform for f∈L1(Ω) is defined by
(12)Rfθ,xθ=Rθfxθ=∫t∈Lθ,xθ∩Ω‍ftdt.Letg∈L1([0,π]×R) and let U¯=(x,y); the adjoint R*:L2(Z)→L2(Ω) of R is given by
(13)R*gU¯=∫0π‍gθ,U¯·θ¯dθ.
The proof can be found in [26].For parallel beam projection, with a beam parallel to theX=(x,y) plane, f∈L1(Σ) and L(θ,xθ,z) the L(θ,xθ) line for the coordinate z:
(14)Iz⊗Rfθ,xθ,z=Rθfxθ=∫t∈Lθ,xθ,z∩Σ‍ftdt,
where ⊗ denotes a tensor product and Iz is the identity on the z variable. The Sobolev spaces Hα are very useful to characterize the smoothing properties of the Radon transform. The Sobolev norms ·Hα for functions in Rn are defined as
(15)fHα2=∫Rn‍1+ξ2αFfξ2dξ.
Let g∈L1([0,π]×R); the one-dimensional Fourier Fp(g)(θ,τ) is defined by Fp(g)(θ,τ)=(1/2π)∫R‍exp⁡-iτpg(θ,p)dp. Let s∈R, the Sobolev space on [0,π]×R denoted as H0,s([0,π]×R), be the set of all distributions g on ([0,π]×R) such that the Fourier transform Fp(g) is a function and such that the Sobolev norm of g, gH0,s(Z)2, is finite [26]. The norm gH0,s(Z)2 is defined as
(16)gH0,s(Z)=∫0π‍∫R‍1+τ2sFpgθ,τ2dτdθ.The Radon operator is a smoothing operator as detailed in the following proposition.Proposition 2.
Letf be a distribution with a compact support; then, the following bound holds for s≥0:
(17)R(f)H0,s+1/2(Z)≤CfHsR2.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.
LetRaz be the range of the Radon operator for the height z; the operator φ=Iz⊗R[δn]:L2(Σ)→L2([0,π]×Σz×∪Raz) is continuous.Proof.
Fors>0, H0s(R2) is continuously embedded in L2(R2) and thus for compactly supported functions, there is a positive constant C1 such that
(18)fHs(R2)≤C1fL2R2.
Using (18), there exists a positive constant C such that
(19)R(f)H0,0(Z)≤R(f)H0,1/2(Z)≤CfL2R2.
In a Hilbert spaces context, the Radon transform is thus a linear continuous operator fromL2(R2) to L2(Z). For parallel beam projection, let Raz be the range of the Radon operator for the height z. For each value of z, the Radon operator is continuous from L2(Σx×Σy) to L2([0,π]×Raz). By integration on the variable z, the operator φ=Iz⊗R[δn]:L2(Σ)→L2([0,π]×Σz×∪Raz) is continuous.With the Radon projection operator, (4) can be rewritten:
(20)φθδnXθ=2πλIz⊗R1-δnθ,Xθ.


## 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 byf the unknown phase to reconstruct, by g the intensity data, and by A 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 Δ:L2(Ω)→H-2(Ω) is continuous [27]. The phase operators involved in the direct problem can be understood as Fourier multiplication operators, with a Fourier multiplier m(k→) and written as
(21)FAfk→=mk→Ffk→.
The Fourier multiplier of the CTF method is given by
(22)mk→=sinπλdk→2.
It is bounded and thus the direct operator in the CTF method is bounded from L2(R2) to L2(R2) 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 multipliersm(k→) involved in the TIE or CTF approaches m(k→)=1/k→2 and mk→=1/sin(πλdk→2) 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 functionu∈Lloc⁡1(Rn-{0}) and given a function ϕ∈Ck∞, the finite part fp of the integral ∫‍u(y)ϕ(y)dy is defined as
(23)fp∫‍uyϕydy=lim⁡ϵ→0∫|y|>ϵ‍u(y)(ϕ(y)-ϕ(0))dy+∫y>1‍uyϕydy.The following proposition is useful to show that some finite part distributions are well defined.Proposition 5.
Let one assume that the functionu∈Lloc⁡1(Rn-{0}) satisfies the size condition
(24)uy≤Cyn,y≠0;
then, it is possible to define the finite part of u [29–31].Thus,f·p(1/k→2) and f·p(1/sin(πλdk→2)) are well-defined tempered distributions. For data g∈L2(Ω), the product m(k→)×g^ defines a tempered distribution in S′ [29, 30]. The inverse Fourier transform of m(k→)×g^ is well defined. The data function g has a bounded support, g∈L2(Ω)⊂L1(Ω)⊂L1(R2), and the convolution of m and g is well defined [29–31].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 operatorA:X→Y between Hilbert spaces X and Y. 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 A~γ:L2(R2)→L2(Ω) with a bounded kernel in L2. Let A~γ be the approximate inverse operator, depending on the regularization parameter γ. This approximate inverse must verify the pointwise convergence, for all f∈L2(Ω):
(25)lim⁡γ→0A~γAf=f.
Let fγ,δ be the solution obtained from the noisy data gδ, fγ,δ=A~γgδ. We assume that the noisy data gδ and the unnoisy data g are such that
(26)g-gδ≤δ,
where δ is the noise level and · is the L2 norm. We are interested in this work in the convergence of various regularization methods. We will make use of the following inequality for f∈L2(R2):
(27)f-fα,δ≤δA~γ+A~γAf-f.


### 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 [14–18] 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 ratef-fα,δ~δ1/3.Proof.
For the Tikhonov regularization, the approximate inverse can be written:(28)A~γ=γI+A*A-1A*,
where A is a Fourier multiplication operator given in (10). The approximate inverse for CTF is also a Fourier multiplication operator:
(29)A~γf=F-1mk→Ffk→.
With u=πλD, the Fourier multiplier m(k→) is given by
(30)mk→=sin⁡uk→2γ+sinuk→22.
We first obtain an upper bound for the norm of the operatorA~γ:
(31)A~γf=F-1sin(uk→2)γ+sin2(uk→2)Ff=sin(uk→2)γ+sin2(uk→2)Ff≤sup⁡k→sinuk→2γ+sin2uk→2f,
where the first equality results from Plancherel’s theorem.
With the identity(32)xγ+x2≤12γ,
we obtain
(33)A~γ≤C1γ.
Forf∈L2(Ω), we have
(34)I-A~γAf=F-1γγ+sinuk→22Ffk→.
With the change of variable, k→=k→′γ1/4 and with the identity sin(x)≥2x/π, for x∈[0,π/2], we obtain
(35)I-A~γAf2=γγ+sinuk→22Ff(k→)2≤γγ+(2/π)uk→22Ffk→2≤∫‍11+2/πuk→′22Ffk→′γ1/42γ1/2dk→′.
From the Rieman-Lebesgue convergence theorem, the Fourier transform is a continuous function. Thus, there is a constant positive C2 such that
(36)(I-A~γA)f≤C2γ1/4.
Reporting in (27), we obtain for the Tikhonov regularization
(37)f-fα,δ≤δC1γ+C2γ1/4,
and thus the optimal convergence rate is obtained for the a priori choice, γ~δ4/3, and the convergence rate is
(38)f-fα,δ~δ1/3.


#### 3.2.2. Quasiparticle Regularization

In the framework of the quasiparticule approach [22–24], the intensity contrast is filtered in the Fourier space F(gz)(k→) and is replaced by FΘ(gz)(k→):
(39)FΘ(gz)(k→)=Θsinπλdk→2-ϵFgzk→,
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 γ>0 such that
(40)sinπλdγ2=ϵ.
Let Hγ(k→) be the filter function defined by
(41)Hγ(k→)=0,ifk→≤γ,1,ifk→>γ.
The regularized inverse can be written:
(42)A~γf=F-1Hγk→sin⁡uk→2Ffk→,A~γf≤F-1Hγ(k→)sin(uk→2)Ff(k→)≤1sin(uγ2)Hγ(k→)Ff(k→)≤1sin⁡uγ2Ffk→=1sin⁡uγ2f.With the identitysin(x)≥2x/π, for x∈[0,π/2], we obtain that there is a constant C such that
(43)A~γ≤Cγ2.Moreover, we have(44)I-A~γAf=F-11-Hγk→Ff,(I-A~γA)f=1∥k→∥≤γFf.
By the Lebesgue dominated convergence theorem, this term tends to zero as γ→0. Yet it is impossible to derive any convergence rate of the method as γ→0 without any assumption on the position of the spectra of f with respect to the threshold. In the more favorable case, if the spectra of the function f are above γ, this term is zero and
(45)x-xα,δ≤δCγ2.
The convergence is obtained for a choice of the regularization parameter such that δ/γ2→0 as δ→0.


#### 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ϕ∈L1(R2) such that
(46)Fϕ(0)=∫R2‍ϕ(x)dx=1,
and we assume that there exists s>0 such that
(47)1-Fϕk→~k→sask→→0.
Typical examples of mollifier are the Gaussian function
(48)ϕx=2π-1exp⁡-x22
or
(49)ϕν(x)=A1-x2ν,ifx≤1,0,ifx>1,
where A is a normalizing constant.We consider as approximate inverse the operator defined forf∈L2(Ω) by
(50)A~γf=F-11-Fϕ(k→/γ)(sin(uk→2))Ff(k).Proposition 7.
The CTF formulation of the phase retrieval inverse problem can be regularized by the mollification regularization with an exponents=2. The optimal convergence rate is f-fα,δ~δ1/2.Proof.
Let us assumes=2; then,
(51)A~γf=F-11-Fϕ(k→/γ)(sin(uk→2))Ff(k)=1-Fϕ(k→/γ)(sin(uk→2))Ff(k)≤sup⁡k1-Fϕ(k→/γ)sinuk→2Ff(k→)=sup⁡k1-Fϕ(k→/γ)sinuk→2f.
From the Rieman-Lebesque theorem,Fϕ is continuous and decreases at infinity and thus |1-Fϕ(k→/γ)|/|sin(uk→2)| is bounded for the higher value of k→.
Ask→→0,
(52)1-Fϕ(k→/γ)sinuk→2~γ-sk→suk→2;
for s=2, there is a constant C1 such that
(53)1-Fϕ(k→/γ)sinuk→2≤C1γ-2.
Thus, we obtain(54)A~γ≤C1γ-2
with the definition of approximate inverse
(55)(I-A~γA)f=F-1(Fϕk→γFf(k→)).
The mollifier is defined on R2 and thus
(56)Fϕk→γ=γ2Fϕxγ;
we obtain
(57)(I-A~γA)f=γ2F(ϕ(xγ))Ff(k)≤γ2Ffk→≤γ2f.
With (27),
(58)f-fα,δ≤δC1γ-2+γ2f;
the optimal convergence rate obtained with the choice γ~δ1/4 is
(59)f-fα,δ~δ1/2.


## 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 Δ:L2(Ω)→H-2(Ω) is continuous [27]. The phase operators involved in the direct problem can be understood as Fourier multiplication operators, with a Fourier multiplier m(k→) and written as
(21)FAfk→=mk→Ffk→.
The Fourier multiplier of the CTF method is given by
(22)mk→=sinπλdk→2.
It is bounded and thus the direct operator in the CTF method is bounded from L2(R2) to L2(R2) 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 multipliersm(k→) involved in the TIE or CTF approaches m(k→)=1/k→2 and mk→=1/sin(πλdk→2) 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 functionu∈Lloc⁡1(Rn-{0}) and given a function ϕ∈Ck∞, the finite part fp of the integral ∫‍u(y)ϕ(y)dy is defined as
(23)fp∫‍uyϕydy=lim⁡ϵ→0∫|y|>ϵ‍u(y)(ϕ(y)-ϕ(0))dy+∫y>1‍uyϕydy.The following proposition is useful to show that some finite part distributions are well defined.Proposition 5.
Let one assume that the functionu∈Lloc⁡1(Rn-{0}) satisfies the size condition
(24)uy≤Cyn,y≠0;
then, it is possible to define the finite part of u [29–31].Thus,f·p(1/k→2) and f·p(1/sin(πλdk→2)) are well-defined tempered distributions. For data g∈L2(Ω), the product m(k→)×g^ defines a tempered distribution in S′ [29, 30]. The inverse Fourier transform of m(k→)×g^ is well defined. The data function g has a bounded support, g∈L2(Ω)⊂L1(Ω)⊂L1(R2), and the convolution of m and g is well defined [29–31].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 operatorA:X→Y between Hilbert spaces X and Y. 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 A~γ:L2(R2)→L2(Ω) with a bounded kernel in L2. Let A~γ be the approximate inverse operator, depending on the regularization parameter γ. This approximate inverse must verify the pointwise convergence, for all f∈L2(Ω):
(25)lim⁡γ→0A~γAf=f.
Let fγ,δ be the solution obtained from the noisy data gδ, fγ,δ=A~γgδ. We assume that the noisy data gδ and the unnoisy data g are such that
(26)g-gδ≤δ,
where δ is the noise level and · is the L2 norm. We are interested in this work in the convergence of various regularization methods. We will make use of the following inequality for f∈L2(R2):
(27)f-fα,δ≤δA~γ+A~γAf-f.


## 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 [14–18] 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 ratef-fα,δ~δ1/3.Proof.
For the Tikhonov regularization, the approximate inverse can be written:(28)A~γ=γI+A*A-1A*,
where A is a Fourier multiplication operator given in (10). The approximate inverse for CTF is also a Fourier multiplication operator:
(29)A~γf=F-1mk→Ffk→.
With u=πλD, the Fourier multiplier m(k→) is given by
(30)mk→=sin⁡uk→2γ+sinuk→22.
We first obtain an upper bound for the norm of the operatorA~γ:
(31)A~γf=F-1sin(uk→2)γ+sin2(uk→2)Ff=sin(uk→2)γ+sin2(uk→2)Ff≤sup⁡k→sinuk→2γ+sin2uk→2f,
where the first equality results from Plancherel’s theorem.
With the identity(32)xγ+x2≤12γ,
we obtain
(33)A~γ≤C1γ.
Forf∈L2(Ω), we have
(34)I-A~γAf=F-1γγ+sinuk→22Ffk→.
With the change of variable, k→=k→′γ1/4 and with the identity sin(x)≥2x/π, for x∈[0,π/2], we obtain
(35)I-A~γAf2=γγ+sinuk→22Ff(k→)2≤γγ+(2/π)uk→22Ffk→2≤∫‍11+2/πuk→′22Ffk→′γ1/42γ1/2dk→′.
From the Rieman-Lebesgue convergence theorem, the Fourier transform is a continuous function. Thus, there is a constant positive C2 such that
(36)(I-A~γA)f≤C2γ1/4.
Reporting in (27), we obtain for the Tikhonov regularization
(37)f-fα,δ≤δC1γ+C2γ1/4,
and thus the optimal convergence rate is obtained for the a priori choice, γ~δ4/3, and the convergence rate is
(38)f-fα,δ~δ1/3.


### 3.2.2. Quasiparticle Regularization

In the framework of the quasiparticule approach [22–24], the intensity contrast is filtered in the Fourier space F(gz)(k→) and is replaced by FΘ(gz)(k→):
(39)FΘ(gz)(k→)=Θsinπλdk→2-ϵFgzk→,
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 γ>0 such that
(40)sinπλdγ2=ϵ.
Let Hγ(k→) be the filter function defined by
(41)Hγ(k→)=0,ifk→≤γ,1,ifk→>γ.
The regularized inverse can be written:
(42)A~γf=F-1Hγk→sin⁡uk→2Ffk→,A~γf≤F-1Hγ(k→)sin(uk→2)Ff(k→)≤1sin(uγ2)Hγ(k→)Ff(k→)≤1sin⁡uγ2Ffk→=1sin⁡uγ2f.With the identitysin(x)≥2x/π, for x∈[0,π/2], we obtain that there is a constant C such that
(43)A~γ≤Cγ2.Moreover, we have(44)I-A~γAf=F-11-Hγk→Ff,(I-A~γA)f=1∥k→∥≤γFf.
By the Lebesgue dominated convergence theorem, this term tends to zero as γ→0. Yet it is impossible to derive any convergence rate of the method as γ→0 without any assumption on the position of the spectra of f with respect to the threshold. In the more favorable case, if the spectra of the function f are above γ, this term is zero and
(45)x-xα,δ≤δCγ2.
The convergence is obtained for a choice of the regularization parameter such that δ/γ2→0 as δ→0.


### 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ϕ∈L1(R2) such that
(46)Fϕ(0)=∫R2‍ϕ(x)dx=1,
and we assume that there exists s>0 such that
(47)1-Fϕk→~k→sask→→0.
Typical examples of mollifier are the Gaussian function
(48)ϕx=2π-1exp⁡-x22
or
(49)ϕν(x)=A1-x2ν,ifx≤1,0,ifx>1,
where A is a normalizing constant.We consider as approximate inverse the operator defined forf∈L2(Ω) by
(50)A~γf=F-11-Fϕ(k→/γ)(sin(uk→2))Ff(k).Proposition 7.
The CTF formulation of the phase retrieval inverse problem can be regularized by the mollification regularization with an exponents=2. The optimal convergence rate is f-fα,δ~δ1/2.Proof.
Let us assumes=2; then,
(51)A~γf=F-11-Fϕ(k→/γ)(sin(uk→2))Ff(k)=1-Fϕ(k→/γ)(sin(uk→2))Ff(k)≤sup⁡k1-Fϕ(k→/γ)sinuk→2Ff(k→)=sup⁡k1-Fϕ(k→/γ)sinuk→2f.
From the Rieman-Lebesque theorem,Fϕ is continuous and decreases at infinity and thus |1-Fϕ(k→/γ)|/|sin(uk→2)| is bounded for the higher value of k→.
Ask→→0,
(52)1-Fϕ(k→/γ)sinuk→2~γ-sk→suk→2;
for s=2, there is a constant C1 such that
(53)1-Fϕ(k→/γ)sinuk→2≤C1γ-2.
Thus, we obtain(54)A~γ≤C1γ-2
with the definition of approximate inverse
(55)(I-A~γA)f=F-1(Fϕk→γFf(k→)).
The mollifier is defined on R2 and thus
(56)Fϕk→γ=γ2Fϕxγ;
we obtain
(57)(I-A~γA)f=γ2F(ϕ(xγ))Ff(k)≤γ2Ffk→≤γ2f.
With (27),
(58)f-fα,δ≤δC1γ-2+γ2f;
the optimal convergence rate obtained with the choice γ~δ1/4 is
(59)f-fα,δ~δ1/2.


## 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 ratef-fα,δ~δ1/3.Proof.
For the Tikhonov regularization, the approximate inverse can be written:(28)A~γ=γI+A*A-1A*,
where A is a Fourier multiplication operator given in (10). The approximate inverse for CTF is also a Fourier multiplication operator:
(29)A~γf=F-1mk→Ffk→.
With u=πλD, the Fourier multiplier m(k→) is given by
(30)mk→=sin⁡uk→2γ+sinuk→22.
We first obtain an upper bound for the norm of the operatorA~γ:
(31)A~γf=F-1sin(uk→2)γ+sin2(uk→2)Ff=sin(uk→2)γ+sin2(uk→2)Ff≤sup⁡k→sinuk→2γ+sin2uk→2f,
where the first equality results from Plancherel’s theorem.
With the identity(32)xγ+x2≤12γ,
we obtain
(33)A~γ≤C1γ.
Forf∈L2(Ω), we have
(34)I-A~γAf=F-1γγ+sinuk→22Ffk→.
With the change of variable, k→=k→′γ1/4 and with the identity sin(x)≥2x/π, for x∈[0,π/2], we obtain
(35)I-A~γAf2=γγ+sinuk→22Ff(k→)2≤γγ+(2/π)uk→22Ffk→2≤∫‍11+2/πuk→′22Ffk→′γ1/42γ1/2dk→′.
From the Rieman-Lebesgue convergence theorem, the Fourier transform is a continuous function. Thus, there is a constant positive C2 such that
(36)(I-A~γA)f≤C2γ1/4.
Reporting in (27), we obtain for the Tikhonov regularization
(37)f-fα,δ≤δC1γ+C2γ1/4,
and thus the optimal convergence rate is obtained for the a priori choice, γ~δ4/3, and the convergence rate is
(38)f-fα,δ~δ1/3.


## 3.2.2. Quasiparticle Regularization

In the framework of the quasiparticule approach [22–24], the intensity contrast is filtered in the Fourier space F(gz)(k→) and is replaced by FΘ(gz)(k→):
(39)FΘ(gz)(k→)=Θsinπλdk→2-ϵFgzk→,
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 γ>0 such that
(40)sinπλdγ2=ϵ.
Let Hγ(k→) be the filter function defined by
(41)Hγ(k→)=0,ifk→≤γ,1,ifk→>γ.
The regularized inverse can be written:
(42)A~γf=F-1Hγk→sin⁡uk→2Ffk→,A~γf≤F-1Hγ(k→)sin(uk→2)Ff(k→)≤1sin(uγ2)Hγ(k→)Ff(k→)≤1sin⁡uγ2Ffk→=1sin⁡uγ2f.With the identitysin(x)≥2x/π, for x∈[0,π/2], we obtain that there is a constant C such that
(43)A~γ≤Cγ2.Moreover, we have(44)I-A~γAf=F-11-Hγk→Ff,(I-A~γA)f=1∥k→∥≤γFf.
By the Lebesgue dominated convergence theorem, this term tends to zero as γ→0. Yet it is impossible to derive any convergence rate of the method as γ→0 without any assumption on the position of the spectra of f with respect to the threshold. In the more favorable case, if the spectra of the function f are above γ, this term is zero and
(45)x-xα,δ≤δCγ2.
The convergence is obtained for a choice of the regularization parameter such that δ/γ2→0 as δ→0.


## 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ϕ∈L1(R2) such that
(46)Fϕ(0)=∫R2‍ϕ(x)dx=1,
and we assume that there exists s>0 such that
(47)1-Fϕk→~k→sask→→0.
Typical examples of mollifier are the Gaussian function
(48)ϕx=2π-1exp⁡-x22
or
(49)ϕν(x)=A1-x2ν,ifx≤1,0,ifx>1,
where A is a normalizing constant.We consider as approximate inverse the operator defined forf∈L2(Ω) by
(50)A~γf=F-11-Fϕ(k→/γ)(sin(uk→2))Ff(k).Proposition 7.
The CTF formulation of the phase retrieval inverse problem can be regularized by the mollification regularization with an exponents=2. The optimal convergence rate is f-fα,δ~δ1/2.Proof.
Let us assumes=2; then,
(51)A~γf=F-11-Fϕ(k→/γ)(sin(uk→2))Ff(k)=1-Fϕ(k→/γ)(sin(uk→2))Ff(k)≤sup⁡k1-Fϕ(k→/γ)sinuk→2Ff(k→)=sup⁡k1-Fϕ(k→/γ)sinuk→2f.
From the Rieman-Lebesque theorem,Fϕ is continuous and decreases at infinity and thus |1-Fϕ(k→/γ)|/|sin(uk→2)| is bounded for the higher value of k→.
Ask→→0,
(52)1-Fϕ(k→/γ)sinuk→2~γ-sk→suk→2;
for s=2, there is a constant C1 such that
(53)1-Fϕ(k→/γ)sinuk→2≤C1γ-2.
Thus, we obtain(54)A~γ≤C1γ-2
with the definition of approximate inverse
(55)(I-A~γA)f=F-1(Fϕk→γFf(k→)).
The mollifier is defined on R2 and thus
(56)Fϕk→γ=γ2Fϕxγ;
we obtain
(57)(I-A~γA)f=γ2F(ϕ(xγ))Ff(k)≤γ2Ffk→≤γ2f.
With (27),
(58)f-fα,δ≤δC1γ-2+γ2f;
the optimal convergence rate obtained with the choice γ~δ1/4 is
(59)f-fα,δ~δ1/2.


## 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 (Section3.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Λα:Ht(Rd)→Ht-α(Rd) for α>-d is defined by
(60)FΛαfk→=ξαFfk→.The Riesz potential is an isomorphism of Hilbert spaces. Forα<0, it is a compact operator. The Riesz potential used for tomography reconstruction, Λs, for g∈Cc∞([0,π]×R), is the operator with Fourier multiplier |ξ|:
(61)Λsg=Fs-1ξFsg,
where the partial Fourier transform of g∈Cc∞([0,π]×R) in the s variable is defined by the formula
(62)Fsgϕ,ξ=12π∫-∞∞‍e-isξgϕ,sds.The inversion formula is summarized by the following theorem [26].Theorem 9.
Letf∈Cc∞(R2), the set of infinitely differentiable functions with 2D compact support; then, f can be reconstructed with the formula
(63)f=14πR*ΛsRf,
where R* is the adjoint of the Radon transform.This reconstruction formula can be generalized for parallel beam projection with tensor product of identity operators on thez coordinate. In this case, ξ is the spatial frequency associated to the xθ axis (Figure 1). For parallel beam projection, the reconstruction formula is
(64)fx=14πIz⊗R*Iθ⊗Iz⊗ΛsIz⊗Rfx.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 asf=(1/4π)R*(Π*Rf), where Π is the inverse Fourier transform of the truncated multiplier |ξ| and * denotes the convolution in the s variable:
(65)FsΠ(ξ)=ξ,ifξ≤ω,0,ifξ>ω,
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(66)f*e=R*v*sRf,
where v is the reconstruction kernel and e is a mollifier approximating the δ distribution, centered about 0 and with mean value 1. The reconstruction kernel and the mollifier are related by
(67)v=ΛRe.
A family of reconstruction kernel is obtained for γ>0 as
(68)eγx=γ-2exγ,vγx=γ-2vxγ.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 gδ, assuming that the range of the Radon transform is [-1,1], the reconstruction formula can be written, for y∈Ω:
(69)R~γgδy=f(y)=γ-2∫0π‍∫-11‍gδθ,svθ,s-y,θ¯γdsdθ.
For a 3D geometry, an approximate inverse of the Radon projector can be obtained with the tensorial product R~γ⊗Iz.


### 4.2. The Approximate Inverse for Phase Contrast Tomography

The phase contrast tomography problem can be written in a generic way as(70)Uf=A⊗IθIz⊗Rf=g,
where f denotes in this section the unknown imaginary index part δn and g is the intensity data. The operator Iz⊗R, φ=Iz⊗R[δn]:L2(Σ)→L2([0,π]×Σz×∪Raz), is linear and continuous. The same result is true for the operator A⊗Iθ:L2([0,π]×Σz×∪Raz)→L2([0,π]×R2).The approximate inverse of the direct operator can be expressed as(71)U~=R~γ⊗IzA~γ⊗Iθ,
where A~γ and R~γ are approximate inverses for the Radon projection operator and the phase retrieval operator A. 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 inverseU~ is given by
(72)U~=(Iz⊗R*)(Iθ⊗Iz⊗Λγ)(Iθ⊗A~)=(Iz⊗R*)(Iθ⊗Iz⊗F1-1ξ1ξ≤ωF1)·Iθ⊗F2-1mξ,ηF2.
In this formula, m(ξ,η) is the Fourier multiplier of the CTF approach regularized with the mollifier approach and ξ and η are the Fourier coordinates associated with the xθ and z 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,f(·,·,z)∈H0,1/2+ϵ(R2). 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 τ<4/3. The optimal convergence rate is obtained as f-fα,δ~δτϵ/(τϵ+τ+2) for τϵ<2-τ and f-fα,δ~δ(2-τ)/4 for τϵ>2-τ.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:
(73)Λγf=F-11|ξ|≤ωξFf=1ξ≤ωξFf≤ωFf=γ-τf.
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 R* is bounded, we get
(74)U~≤C1γ-τ-2.
Ignoring tensorial products with identity operators for the sake of simplicity, we have the equality
(75)I-U~U=I-R*Λγ(A~A-I)R+R*ΛR+R*(Λγ-Λ)R=-R*ΛγA~A-IR+R*Λγ-ΛR;
we obtain
(76)(I-U~U)f≤R*ΛγA~A-IRf+R*Λγ-ΛRf.
Let us assume that f∈H1/2+ϵ(R2), with a constant ϵ>0; then, g=Rf∈H0,1+ϵ(Z), and |ξ|1+ϵFg∈L2 (Proposition 2):
(77)Λγ-ΛRf=1|ξ|≥ωξFRf=1|ξ|≥ωξ1+ϵξ-ϵFRf≤1|ξ|≥ωξ1+ϵξ-ϵFRf≤sup⁡|ξ|≥ωξ-ϵξ1+ϵFRf≤γτϵξ1+ϵFRf≤C2γτϵ,
where C2 is a positive constant. The same bound can be obtained for 3D geometry with parallel beam projection. With the bound obtained for (A~A-I)f for mollifier regularization, we have
(78)(I-U~U)f≤R*Λr(A~A-I)Rf+R*(Λr-Λ)Rf≤C2max⁡(γ2-τ,γτϵ)≤C2γmin⁡⁡2-τ,τϵ.
We thus obtain the following bound:(79)f-fα,δ≤C1δγ-τ-2+C2γmin⁡⁡2-τ,τϵ.
The method is convergent for τ<2. The optimal convergence rate is obtained as f-fα,δ~δτϵ/(τϵ+τ+2) for τϵ<2-τ and the choice γ=δ1/(τϵ+τ+2) and f-fα,δ~δ(2-τ)/4 for τϵ>2-τ and the choice γ=δ1/4.


#### 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 e*f converge to f in L2 as γ→0. To have a more precise bound, we will use the following lemma about difference quotients in Sobolev spaces [27].Definition 11.
Given a unit vectorv and a measurable function u, the h-translate of u in the direction v is denoted as τvhu and defined by
(80)τvhux=ux+vh.
For a measurable function u, the difference quotient of u in the coordinate direction ei of length h≠0 is the measurable function defined by
(81)Δihu=τeihu-uh.Proposition 12.
Letu∈W1,p(Ω′) and let Du be the weak derivative of u. Then, Δhu∈Lp(Ω) for any Ω⊂Ω′, with h<dist⁡(Ω,∂Ω′) and
(82)Δhup;Ω≤Dup;Ω′.The convergence results are summarized by the following proposition.Proposition 13.
Let us assume thatf can be extended to a larger domain Ω′ such that f∈W1,2(Ω′). 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 f-fα,δ~δ1/8.Proof.
We first bound the norm ofR~γ. We assume that the range of the Radon transform is V=[-1,1], with the Young and Schwartz inequalities we obtain, for y∈Ω:
(83)R~γgy2=γ-2∫0π‍∫-11‍gθ,svθ,s-y,θ¯γdsdθ2≤∫0π‍gθ,sL2Vvγθ,sL2Vdθ2≤∫0π‍gθ,sL2V2dθ∫0π‍vγθ,sL2V2dθ.
The Fourier slice theorem can be written [32]:
(84)12πFsRfθ,σ=Ffσθ¯=f^σθ¯,
where Fs denotes the Fourier transform with respect to the s variable and Ff=f^ is the 2D Fourier transform; thus,
(85)ΛsRθeγs=2πFs-1σe^γσθ¯,ΛsRθeγL2V2=2πσe^γ(σθ¯)L2(V)2=2π∫‍σ2e^γσθ¯2dσ=2π∫‍σ2e^γσθ¯2dσ=2πγ-3∫‍σ2e^σθ¯2dσ,∫0π‍ΛsRθeγL2V2dθ=γ-3∫0π‍2π∫‍σ2e^σθ¯2dσdθ=γ-3∫R2‍k→e^k→2dk→≤γ-3∫R2‍1+k→212e^k→2dk→=γ-3e1/2,
with a polar change of variable. By integration on the variable y, there is a positive constant C such that
(86)R~γ≤Cγ-3/2.
We first bound the norm eγ*f-f:
(87)eγ*f(x)-f(x)=∫‍eγ(y)(f(x-y)-f(x))dy=γ-2∫‍eyγfx-y-fxdy=∫‍ezfx-γz-fxdz.
With the Hölder inequality, we obtain(88)eγ*fx-fx≤∫‍ezfx-γz-fx2dz1/2∫‍ezdz1/2≤∫‍ezfx-γz-fx2dz1/2.
By integration,(89)eγ*f-f2≤∫‍ezfx-γz-fx2dz.
Let us note d(Ω), the diameter of the domain Ω.
Let us assume thatf can be extended to a larger bounded domain Ω′, such that f∈W1,p(Ω′) with γd(Ω)<dist⁡(Ω,∂Ω′). Then, with proposition (8) with p=2, we obtain
(90)eγ*f-f≤C2γ,
for a positive constant C2.
The regularized inverse can be written asU~=(Iz⊗R~γ)(Iθ⊗A~γ). Using (86) and (54) for the mollifiers regularization of CTF, we thus obtain the bound
(91)U~≤C1γ-7/2,
for a positive constant C1. From the identity
(92)I-U~U=I-R~γR-RγAγA-IR,
and with (90) and (57), we obtain for C2≥0(93)(I-U~U)f≤C2max⁡(γ,γ1/2)≤C2γ1/2.
With (27), we obtain
(94)f-fα,δ≤C1δγ-7/2+C2γ1/2.
The optimal convergence rate obtained is f-fα,δ~δ1/8 with the choice γ=δ1/8.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 ofA~γA-I for the Tikhonov method (36) and for the norm of Rγ  (86), (92) shows that the convergence of I-U~U towards zero is not achieved.


## 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Λα:Ht(Rd)→Ht-α(Rd) for α>-d is defined by
(60)FΛαfk→=ξαFfk→.The Riesz potential is an isomorphism of Hilbert spaces. Forα<0, it is a compact operator. The Riesz potential used for tomography reconstruction, Λs, for g∈Cc∞([0,π]×R), is the operator with Fourier multiplier |ξ|:
(61)Λsg=Fs-1ξFsg,
where the partial Fourier transform of g∈Cc∞([0,π]×R) in the s variable is defined by the formula
(62)Fsgϕ,ξ=12π∫-∞∞‍e-isξgϕ,sds.The inversion formula is summarized by the following theorem [26].Theorem 9.
Letf∈Cc∞(R2), the set of infinitely differentiable functions with 2D compact support; then, f can be reconstructed with the formula
(63)f=14πR*ΛsRf,
where R* is the adjoint of the Radon transform.This reconstruction formula can be generalized for parallel beam projection with tensor product of identity operators on thez coordinate. In this case, ξ is the spatial frequency associated to the xθ axis (Figure 1). For parallel beam projection, the reconstruction formula is
(64)fx=14πIz⊗R*Iθ⊗Iz⊗ΛsIz⊗Rfx.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 asf=(1/4π)R*(Π*Rf), where Π is the inverse Fourier transform of the truncated multiplier |ξ| and * denotes the convolution in the s variable:
(65)FsΠ(ξ)=ξ,ifξ≤ω,0,ifξ>ω,
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(66)f*e=R*v*sRf,
where v is the reconstruction kernel and e is a mollifier approximating the δ distribution, centered about 0 and with mean value 1. The reconstruction kernel and the mollifier are related by
(67)v=ΛRe.
A family of reconstruction kernel is obtained for γ>0 as
(68)eγx=γ-2exγ,vγx=γ-2vxγ.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 gδ, assuming that the range of the Radon transform is [-1,1], the reconstruction formula can be written, for y∈Ω:
(69)R~γgδy=f(y)=γ-2∫0π‍∫-11‍gδθ,svθ,s-y,θ¯γdsdθ.
For a 3D geometry, an approximate inverse of the Radon projector can be obtained with the tensorial product R~γ⊗Iz.


## 4.2. The Approximate Inverse for Phase Contrast Tomography

The phase contrast tomography problem can be written in a generic way as(70)Uf=A⊗IθIz⊗Rf=g,
where f denotes in this section the unknown imaginary index part δn and g is the intensity data. The operator Iz⊗R, φ=Iz⊗R[δn]:L2(Σ)→L2([0,π]×Σz×∪Raz), is linear and continuous. The same result is true for the operator A⊗Iθ:L2([0,π]×Σz×∪Raz)→L2([0,π]×R2).The approximate inverse of the direct operator can be expressed as(71)U~=R~γ⊗IzA~γ⊗Iθ,
where A~γ and R~γ are approximate inverses for the Radon projection operator and the phase retrieval operator A. 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 inverseU~ is given by
(72)U~=(Iz⊗R*)(Iθ⊗Iz⊗Λγ)(Iθ⊗A~)=(Iz⊗R*)(Iθ⊗Iz⊗F1-1ξ1ξ≤ωF1)·Iθ⊗F2-1mξ,ηF2.
In this formula, m(ξ,η) is the Fourier multiplier of the CTF approach regularized with the mollifier approach and ξ and η are the Fourier coordinates associated with the xθ and z 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,f(·,·,z)∈H0,1/2+ϵ(R2). 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 τ<4/3. The optimal convergence rate is obtained as f-fα,δ~δτϵ/(τϵ+τ+2) for τϵ<2-τ and f-fα,δ~δ(2-τ)/4 for τϵ>2-τ.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:
(73)Λγf=F-11|ξ|≤ωξFf=1ξ≤ωξFf≤ωFf=γ-τf.
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 R* is bounded, we get
(74)U~≤C1γ-τ-2.
Ignoring tensorial products with identity operators for the sake of simplicity, we have the equality
(75)I-U~U=I-R*Λγ(A~A-I)R+R*ΛR+R*(Λγ-Λ)R=-R*ΛγA~A-IR+R*Λγ-ΛR;
we obtain
(76)(I-U~U)f≤R*ΛγA~A-IRf+R*Λγ-ΛRf.
Let us assume that f∈H1/2+ϵ(R2), with a constant ϵ>0; then, g=Rf∈H0,1+ϵ(Z), and |ξ|1+ϵFg∈L2 (Proposition 2):
(77)Λγ-ΛRf=1|ξ|≥ωξFRf=1|ξ|≥ωξ1+ϵξ-ϵFRf≤1|ξ|≥ωξ1+ϵξ-ϵFRf≤sup⁡|ξ|≥ωξ-ϵξ1+ϵFRf≤γτϵξ1+ϵFRf≤C2γτϵ,
where C2 is a positive constant. The same bound can be obtained for 3D geometry with parallel beam projection. With the bound obtained for (A~A-I)f for mollifier regularization, we have
(78)(I-U~U)f≤R*Λr(A~A-I)Rf+R*(Λr-Λ)Rf≤C2max⁡(γ2-τ,γτϵ)≤C2γmin⁡⁡2-τ,τϵ.
We thus obtain the following bound:(79)f-fα,δ≤C1δγ-τ-2+C2γmin⁡⁡2-τ,τϵ.
The method is convergent for τ<2. The optimal convergence rate is obtained as f-fα,δ~δτϵ/(τϵ+τ+2) for τϵ<2-τ and the choice γ=δ1/(τϵ+τ+2) and f-fα,δ~δ(2-τ)/4 for τϵ>2-τ and the choice γ=δ1/4.


### 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 e*f converge to f in L2 as γ→0. To have a more precise bound, we will use the following lemma about difference quotients in Sobolev spaces [27].Definition 11.
Given a unit vectorv and a measurable function u, the h-translate of u in the direction v is denoted as τvhu and defined by
(80)τvhux=ux+vh.
For a measurable function u, the difference quotient of u in the coordinate direction ei of length h≠0 is the measurable function defined by
(81)Δihu=τeihu-uh.Proposition 12.
Letu∈W1,p(Ω′) and let Du be the weak derivative of u. Then, Δhu∈Lp(Ω) for any Ω⊂Ω′, with h<dist⁡(Ω,∂Ω′) and
(82)Δhup;Ω≤Dup;Ω′.The convergence results are summarized by the following proposition.Proposition 13.
Let us assume thatf can be extended to a larger domain Ω′ such that f∈W1,2(Ω′). 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 f-fα,δ~δ1/8.Proof.
We first bound the norm ofR~γ. We assume that the range of the Radon transform is V=[-1,1], with the Young and Schwartz inequalities we obtain, for y∈Ω:
(83)R~γgy2=γ-2∫0π‍∫-11‍gθ,svθ,s-y,θ¯γdsdθ2≤∫0π‍gθ,sL2Vvγθ,sL2Vdθ2≤∫0π‍gθ,sL2V2dθ∫0π‍vγθ,sL2V2dθ.
The Fourier slice theorem can be written [32]:
(84)12πFsRfθ,σ=Ffσθ¯=f^σθ¯,
where Fs denotes the Fourier transform with respect to the s variable and Ff=f^ is the 2D Fourier transform; thus,
(85)ΛsRθeγs=2πFs-1σe^γσθ¯,ΛsRθeγL2V2=2πσe^γ(σθ¯)L2(V)2=2π∫‍σ2e^γσθ¯2dσ=2π∫‍σ2e^γσθ¯2dσ=2πγ-3∫‍σ2e^σθ¯2dσ,∫0π‍ΛsRθeγL2V2dθ=γ-3∫0π‍2π∫‍σ2e^σθ¯2dσdθ=γ-3∫R2‍k→e^k→2dk→≤γ-3∫R2‍1+k→212e^k→2dk→=γ-3e1/2,
with a polar change of variable. By integration on the variable y, there is a positive constant C such that
(86)R~γ≤Cγ-3/2.
We first bound the norm eγ*f-f:
(87)eγ*f(x)-f(x)=∫‍eγ(y)(f(x-y)-f(x))dy=γ-2∫‍eyγfx-y-fxdy=∫‍ezfx-γz-fxdz.
With the Hölder inequality, we obtain(88)eγ*fx-fx≤∫‍ezfx-γz-fx2dz1/2∫‍ezdz1/2≤∫‍ezfx-γz-fx2dz1/2.
By integration,(89)eγ*f-f2≤∫‍ezfx-γz-fx2dz.
Let us note d(Ω), the diameter of the domain Ω.
Let us assume thatf can be extended to a larger bounded domain Ω′, such that f∈W1,p(Ω′) with γd(Ω)<dist⁡(Ω,∂Ω′). Then, with proposition (8) with p=2, we obtain
(90)eγ*f-f≤C2γ,
for a positive constant C2.
The regularized inverse can be written asU~=(Iz⊗R~γ)(Iθ⊗A~γ). Using (86) and (54) for the mollifiers regularization of CTF, we thus obtain the bound
(91)U~≤C1γ-7/2,
for a positive constant C1. From the identity
(92)I-U~U=I-R~γR-RγAγA-IR,
and with (90) and (57), we obtain for C2≥0(93)(I-U~U)f≤C2max⁡(γ,γ1/2)≤C2γ1/2.
With (27), we obtain
(94)f-fα,δ≤C1δγ-7/2+C2γ1/2.
The optimal convergence rate obtained is f-fα,δ~δ1/8 with the choice γ=δ1/8.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 ofA~γA-I for the Tikhonov method (36) and for the norm of Rγ  (86), (92) shows that the convergence of I-U~U towards zero is not achieved.


## 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 inverseU~ is given by
(72)U~=(Iz⊗R*)(Iθ⊗Iz⊗Λγ)(Iθ⊗A~)=(Iz⊗R*)(Iθ⊗Iz⊗F1-1ξ1ξ≤ωF1)·Iθ⊗F2-1mξ,ηF2.
In this formula, m(ξ,η) is the Fourier multiplier of the CTF approach regularized with the mollifier approach and ξ and η are the Fourier coordinates associated with the xθ and z 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,f(·,·,z)∈H0,1/2+ϵ(R2). 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 τ<4/3. The optimal convergence rate is obtained as f-fα,δ~δτϵ/(τϵ+τ+2) for τϵ<2-τ and f-fα,δ~δ(2-τ)/4 for τϵ>2-τ.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:
(73)Λγf=F-11|ξ|≤ωξFf=1ξ≤ωξFf≤ωFf=γ-τf.
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 R* is bounded, we get
(74)U~≤C1γ-τ-2.
Ignoring tensorial products with identity operators for the sake of simplicity, we have the equality
(75)I-U~U=I-R*Λγ(A~A-I)R+R*ΛR+R*(Λγ-Λ)R=-R*ΛγA~A-IR+R*Λγ-ΛR;
we obtain
(76)(I-U~U)f≤R*ΛγA~A-IRf+R*Λγ-ΛRf.
Let us assume that f∈H1/2+ϵ(R2), with a constant ϵ>0; then, g=Rf∈H0,1+ϵ(Z), and |ξ|1+ϵFg∈L2 (Proposition 2):
(77)Λγ-ΛRf=1|ξ|≥ωξFRf=1|ξ|≥ωξ1+ϵξ-ϵFRf≤1|ξ|≥ωξ1+ϵξ-ϵFRf≤sup⁡|ξ|≥ωξ-ϵξ1+ϵFRf≤γτϵξ1+ϵFRf≤C2γτϵ,
where C2 is a positive constant. The same bound can be obtained for 3D geometry with parallel beam projection. With the bound obtained for (A~A-I)f for mollifier regularization, we have
(78)(I-U~U)f≤R*Λr(A~A-I)Rf+R*(Λr-Λ)Rf≤C2max⁡(γ2-τ,γτϵ)≤C2γmin⁡⁡2-τ,τϵ.
We thus obtain the following bound:(79)f-fα,δ≤C1δγ-τ-2+C2γmin⁡⁡2-τ,τϵ.
The method is convergent for τ<2. The optimal convergence rate is obtained as f-fα,δ~δτϵ/(τϵ+τ+2) for τϵ<2-τ and the choice γ=δ1/(τϵ+τ+2) and f-fα,δ~δ(2-τ)/4 for τϵ>2-τ and the choice γ=δ1/4.


## 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 e*f converge to f in L2 as γ→0. To have a more precise bound, we will use the following lemma about difference quotients in Sobolev spaces [27].Definition 11.
Given a unit vectorv and a measurable function u, the h-translate of u in the direction v is denoted as τvhu and defined by
(80)τvhux=ux+vh.
For a measurable function u, the difference quotient of u in the coordinate direction ei of length h≠0 is the measurable function defined by
(81)Δihu=τeihu-uh.Proposition 12.
Letu∈W1,p(Ω′) and let Du be the weak derivative of u. Then, Δhu∈Lp(Ω) for any Ω⊂Ω′, with h<dist⁡(Ω,∂Ω′) and
(82)Δhup;Ω≤Dup;Ω′.The convergence results are summarized by the following proposition.Proposition 13.
Let us assume thatf can be extended to a larger domain Ω′ such that f∈W1,2(Ω′). 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 f-fα,δ~δ1/8.Proof.
We first bound the norm ofR~γ. We assume that the range of the Radon transform is V=[-1,1], with the Young and Schwartz inequalities we obtain, for y∈Ω:
(83)R~γgy2=γ-2∫0π‍∫-11‍gθ,svθ,s-y,θ¯γdsdθ2≤∫0π‍gθ,sL2Vvγθ,sL2Vdθ2≤∫0π‍gθ,sL2V2dθ∫0π‍vγθ,sL2V2dθ.
The Fourier slice theorem can be written [32]:
(84)12πFsRfθ,σ=Ffσθ¯=f^σθ¯,
where Fs denotes the Fourier transform with respect to the s variable and Ff=f^ is the 2D Fourier transform; thus,
(85)ΛsRθeγs=2πFs-1σe^γσθ¯,ΛsRθeγL2V2=2πσe^γ(σθ¯)L2(V)2=2π∫‍σ2e^γσθ¯2dσ=2π∫‍σ2e^γσθ¯2dσ=2πγ-3∫‍σ2e^σθ¯2dσ,∫0π‍ΛsRθeγL2V2dθ=γ-3∫0π‍2π∫‍σ2e^σθ¯2dσdθ=γ-3∫R2‍k→e^k→2dk→≤γ-3∫R2‍1+k→212e^k→2dk→=γ-3e1/2,
with a polar change of variable. By integration on the variable y, there is a positive constant C such that
(86)R~γ≤Cγ-3/2.
We first bound the norm eγ*f-f:
(87)eγ*f(x)-f(x)=∫‍eγ(y)(f(x-y)-f(x))dy=γ-2∫‍eyγfx-y-fxdy=∫‍ezfx-γz-fxdz.
With the Hölder inequality, we obtain(88)eγ*fx-fx≤∫‍ezfx-γz-fx2dz1/2∫‍ezdz1/2≤∫‍ezfx-γz-fx2dz1/2.
By integration,(89)eγ*f-f2≤∫‍ezfx-γz-fx2dz.
Let us note d(Ω), the diameter of the domain Ω.
Let us assume thatf can be extended to a larger bounded domain Ω′, such that f∈W1,p(Ω′) with γd(Ω)<dist⁡(Ω,∂Ω′). Then, with proposition (8) with p=2, we obtain
(90)eγ*f-f≤C2γ,
for a positive constant C2.
The regularized inverse can be written asU~=(Iz⊗R~γ)(Iθ⊗A~γ). Using (86) and (54) for the mollifiers regularization of CTF, we thus obtain the bound
(91)U~≤C1γ-7/2,
for a positive constant C1. From the identity
(92)I-U~U=I-R~γR-RγAγA-IR,
and with (90) and (57), we obtain for C2≥0(93)(I-U~U)f≤C2max⁡(γ,γ1/2)≤C2γ1/2.
With (27), we obtain
(94)f-fα,δ≤C1δγ-7/2+C2γ1/2.
The optimal convergence rate obtained is f-fα,δ~δ1/8 with the choice γ=δ1/8.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 ofA~γA-I for the Tikhonov method (36) and for the norm of Rγ  (86), (92) shows that the convergence of I-U~U 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.



---

*Source: 943501-2015-02-04.xml*