![MathJax Logo](/templates/jsp/_style2/_tandf/pb2/images/math-jax.gif)
Abstract
A new nonlinear method to reconstruct the complex refractive index distribution with in-line phase tomography measurements is presented. The inverse problem is regularized with the Tikhonov smoothing norm of the index. The original nonlinear iterative approach is based on the Fréchet derivative of the intensity recorded at a single propagation distance and on a Simultaneous Algebraic Reconstruction technique. The reconstruction method requires no a priori knowledge about the materials. The algorithm is successfully applied to some simulated data with noise.
1 Introduction
X-ray microtomography is a very useful imaging technique in material science [Citation1] and medical imaging.[Citation2–Citation6] Yet, its sensitivity is now largely 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 a coupling of tomography and phase retrieval to reconstruct the real part of the refractive index.[Citation7, Citation8] 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 ). 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.[Citation9–Citation12] Several linear algorithms have been proposed for the phase retrieval from the Fresnel diffraction patterns [Citation1, Citation10, Citation13–Citation18] valid under some restrictive assumptions. More recently, some nonlinear approaches have been developed [Citation12, Citation19, Citation20] that improve the reconstruction results.
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.
![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.](/cms/asset/51a34f0c-a531-4da8-b659-45cc4c3db9b4/gipe_a_947480_f0001_b.gif)
In order to reconstruct the refractive index in phase contrast tomography, the prior that the ratio of the imaginary to the real part of the complex refractive index is piecewise constant can be included.[Citation8, Citation11, Citation21] If the attenuation image can be measured, it is possible to use a scaled version of this image to construct an a prior knowledge of the phase for multimaterial objects.[Citation8, Citation11] A linear algorithm has also been proposed for multimaterial objects to recover the distance travelled in each material and the real and complex part of the refractive index.[Citation21, Citation22]
All these schemes rely on restrictive assumptions and formulate the phase contrast tomography as a linear inverse problem. The reconstruction images obtained need improvements. The simultaneous nonlinear reconstruction of the real and of the imaginary part of the index from noisy intensity measurements for different projection angles in in-line phase tomography has not been theoretically investigated. The coupling of phase retrieval with tomography sets an ill-posed nonlinear inverse problem. Several regularization methods can be investigated. In this work, we investigate the functional properties of the direct operator and we present new reconstruction results obtained for a single propagation distance with a nonlinear Tikhonov regularization scheme coupled with Simultaneous Algebraic Reconstruction technique.
An extensive comparison with the linear state of the art algorithms in the literature on experimental data is beyond the scope of this study. The new nonlinear reconstruction method is tested on a simulated multi-material object, with an object of a given material embedded in a matrix of a second material. In this case, the complex refractive index takes two distinct values.
We first study the direct problem of the image formation. Then we investigate the functional properties of the direct intensity operator and detail the nonlinear Tikhonov regularization method used. The numerical results obtained with a coupling of a modified Landweber method with Simultaneous Algebraic Reconstruction techniques (ART) on noisy simulated multimaterial data are presented and discussed before concluding remarks.
2 The direct problem of the image formation
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
be the rotated spatial coordinate system for an angle
around the z axis (Figure ).
2.1 The Fresnel intensity
The optical properties of an object interacting with coherent X-rays of wavelength are related to its complex refractive index given by:
1
1 where
is the refractive index decrement and
is the absorption index for the spatial coordinates
in the object domain.[Citation23] 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
:
2
2 where
is the absorption and
the phase shift induced by the object for the projection angle
.[Citation8] The phase and the negative logarithm of the absorption,
, are the projections of the absorption index and refraction index, respectively:
3
3
4
4 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 [Citation8, Citation11, Citation24]:
5
5 where
6
6 The intensity
operator can thus be considered as a nonlinear function of
and
.
2.2 The coupling with Radon projection operator
The direct intensity operator can be rewritten with the Radon projection operator. We first summarize some properties of this projection operator.[Citation25] 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 2.1
Let the line defined by
, with
and
,
the range of the Radon transform and
, the Radon transform for
is defined by:
7
7
Proposition 2.2
Let and
, the adjoint
of
is given by:
8
8
Proof
The proof can be found in [Citation25].
For parallel beam projection, with a beam parallel to the plane,
and
the
line for the coordinate
.
9
9
Definition 2.3
Let , and the one-dimensional Fourier
defined by
. Let
, the Sobolev space on
denoted as
, is the set of all distributions g on (
such that the Fourier transform
is a function and such that the Sobolev norm of g,
, is finite.[Citation25] The norm
is defined as:
10
10
The Radon operator is a smoothing operator as detailed in the following proposition.
Proposition 2.4
Let be a distribution with a compact support, then the following bound holds for
:
11
11 where the constant
depend on the size of the support, and
is the classical Sobolev space of order s of functions such that all partial derivatives of f of order up to s are square integrable.
Proof
The proof of this proposition can be found in [Citation25].
It is now possible to use the mapping properties of the Radon transform in the context of phase contrast tomography and in an Hilbert spaces framework.
Proposition 2.5
Let be range of the Radon operator for the height z, the operators
is continuous. The same result is true for the operator
.
Proof
For ,
is continuously embedded in
and thus for compactly supported functions, there is a positive constant
such that:
12
12 Using Equation Equation11
11
11 , there exists a positive constant
such that:
13
13 In an Hilbert spaces context, the Radon transform is thus a linear continuous operator from
to
. For parallel beam projection, let
be range of the Radon operator for the height z. For each value of z, the Radon operator is continuous from
to
. By integration on the variable
, the operator
is continuous. The same result is holds for the operator
.
With the Radon projection operator, Equations Equation33
3 and Equation4
4
4 can be rewritten:
14
14
15
15 Let
the multiplication operator,
the convolution with
and
defined by
. With the above-defined operators and Equations Equation2
2
2 , Equation5
5
5 , Equation13
13
13 and Equation14
14
14 , we obtain:
Proposition 2.6
The nonlinear intensity operator admits the following decomposition, for
:
16
16
3 The Fréchet derivative of the Fresnel intensity operator and its adjoint
The minimization of the regularization functional can be performed with iterative methods and it is useful to calculate the Fréchet derivative of the intensity operator and its adjoint.
Proposition 3.1
Let , the Fréchet derivative of the operator
at the point
is the operator
defined by:
17
17
Proof
The Fréchet derivative can be calculated with the decomposition given Equation Equation1313
13 and follows from the chain rule. The operators
and
are linear and continuous operators and they are identical to their Fréchet derivatives. The Fréchet derivative of the operator
at the point
is given for
by:
18
18 with
.
Proposition 3.2
The adjoint is thus defined by
with:
19a
19a
19b
19b
Proof
This proposition follows from the definition of the Hilbert space adjoint of the Fresnel diffraction intensity Fréchet derivative defined for all , and for all
by:
20
20
4 Nonlinear Tikhonov regularization for phase contrast tomography
4.1 Nonlinear Tikhonov regularization
Our aim is thus to estimate the two components , from the noisy intensity
measurements such that
21
21 where
is the noise free intensity data, and
the noise level.
The Radon operator and the convolution by the Fresnel operator are compact operators. Our inverse problem is ill-posed and must be regularized. In order to obtain stable solutions, we have used a nonlinear Tikhonov regularization, which is a widely used approach for solving nonlinear inverse problems in a stable way. The aim of the nonlinear Tikhonov regularization is to minimize the functional of the two components
:
22
22 where
is a regularization parameter. In the following, we will make use of the stability and convergence result of Engl et al. [Citation26–Citation28] for regularization methods in Hilbert spaces. We first detail the functional properties of the forward intensity operator for phase tomography.
Proposition 4.1
The Fresnel intensity operator is continuous, weakly sequentially closed and Fréchet differentiable with a Lipschitz continuous Fréchet derivative on its domain.
Proof
In [Citation20], it was shown that the operator is continuous, weakly sequentially closed and Fréchet differentiable with a Lipschitz continuous Fréchet derivative on its domain. In Section 2, it was shown that the Radon operator
is linear and continuous between well-defined Hilbert spaces and thus the former properties are satisfied for
.
Assume that a minimum norm solution exists with
, then the following convergence and stability results hold:
Theorem 4.2
Assume , and
for all k. Then
has a convergent subsequence. Every convergent subsequence converges to a minimizer of
.
Theorem 4.3
Assume that satisfies
and
as
. Let
be a sequence of positive numbers converging to 0 such that the data satisfy
. Let
. Then
has a convergent subsequence. The limit is a
-minimal norm solution .
Proof
The operator is continuous and weakly sequentially closed and the result follows from the classical Tikhonov regularization in Hilbert spaces.[Citation27, Citation28]
Some convergence rate can also be derived under suitable source conditions.[Citation27, Citation28]
4.2 Implementation of the nonlinear regularization method
In order to minimize the Tikhonov regularization functional, many iterative methods have been investigated based on Landweber or Gauss–Newton approaches.[Citation27, Citation28] In the field of discrete ART for tomography, several schemes have been studied which are some variants of the Landweber methods.[Citation25, Citation29, Citation30] In this work, we use a variant of the Landweber method and of the Simultaneous Algebraic Reconstruction Technique where the update is performed per projection angle.
The first order optimality condition for the minimization of the regularization functional can be written:23
23
Table 1. Values of the and
values for the materials in the object, for 24 keV X-rays, from http://henke.lbl.gov/optical_constants.
Figure 5. Diagonal profiles of the map:intial map (blue), true map (red), reconstructed map (Black) (PPSNR
24 dB).
![Figure 5. Diagonal profiles of the β map:intial map (blue), true map (red), reconstructed map (Black) (PPSNR = 24 dB).](/cms/asset/dae41261-a4df-43d4-9113-d4e70982bcf3/gipe_a_947480_f0005_oc.gif)
Figure 6. Diagonal profiles of the map:intial map (blue), true map (red), reconstructed map (black) (PPSNR
24 dB).
![Figure 6. Diagonal profiles of the δ map:intial map (blue), true map (red), reconstructed map (black) (PPSNR = 24 dB).](/cms/asset/7140f07a-1547-40e8-8d20-32da84b4146a/gipe_a_947480_f0006_oc.gif)
Figure 7. Evolution of the partial data term corresponding to the selected projection angle as a function of the iterations number without noise (full line) and with noise (dotted line) for a dB.
![Figure 7. Evolution of the partial data term corresponding to the selected projection angle as a function of the iterations number without noise (full line) and with noise (dotted line) for a PPSNR=24 dB.](/cms/asset/bcd63902-1454-4e51-8dc9-0ec8ccdc7ac3/gipe_a_947480_f0007_b.gif)
Figure 8. Evolution of the normalized mean square error on as a function of the iterations number without noise (full line) and with noise (dotted line) for a PPSNR
24 dB.
![Figure 8. Evolution of the normalized mean square error on δ as a function of the iterations number without noise (full line) and with noise (dotted line) for a PPSNR = 24 dB.](/cms/asset/6cc3ac68-370c-47a8-aad8-0ad242699c46/gipe_a_947480_f0008_b.gif)
Figure 9. Evolution of the normalized mean square error on as a function of the iterations number without noise (full line) and with noise (dotted line) for a PPSNR
24 dB.
![Figure 9. Evolution of the normalized mean square error on β as a function of the iterations number without noise (full line) and with noise (dotted line) for a PPSNR = 24 dB.](/cms/asset/e2753278-2b86-4a06-a631-2e2990e57494/gipe_a_947480_f0009_b.gif)
The update formula in our Landweber type method becomes:24
24 In our implementation, to avoid large data storage, one projection angle is first iteratively chosen at random for the calculation of the adjoint of the Radon projection operator. The step length parameter
is chosen with a backtracking strategy in order to obtain the best decrease of the Tikhonov’s functional along the descent direction
, with
. This type of method has been used for linear ill-posed problems.[Citation31] It corresponds to a modification of the standard Landweber method where the step length is constant. In the nonlinear case, the convergence results are established for a small constant step length.[Citation32] Yet, in our simulations the best numerical results were achieved with an optimized descent parameter. The choice of the regularization parameter
and the stopping index of the iterations
are based on the Morozov’s discrepancy principle.[Citation27] The regularization parameter is chosen such that the discrepancy term satisfy the condition:
25
25 for a constant
.
5 Results and discussion
5.1 Simulation details
The nonlinear regularization method used for the recovery of the real and imaginary part of the refraction index has been tested on simulated data with an homogeneous material embedded in another one. The simulated object consists of an Al cylinder of in diameter and
in height embedded in another Poly(methyl-methacrylate (
) cylinder with diameter
. Let
, the
and
values for the materials in the object, for 24 keV X-rays, are listed in Table . The
and
values were discretized on a regular grid with a pixel size of
. The cylinders are included in a rectangular volume of size
with
and
used for the simulations. Horizontal sections of the original
and
maps to be retrieved are displayed in Figure .
The discrete approximation of the Radon transform is the projection operator implemented in the Mablab Toolbox. The nonlinear algorithm was implemented with 400 projections angles. The Fresnel diffraction image was simulated with a hard X-ray beam of wavelength (24 keV) at one sample-to-detector distances
. An initialization index map was obtained by smoothing the true index map. Horizontal sections of the initial
and
functions are displayed in Figure . The intensity data were corrupted with additive Gaussian white noise with a peak-to-peak signal-to-noise ratio (PPSNR) between 0 and 24 dB.
5.2 Results
Figure displays the horizontal section of the real and imaginary part of the reconstructed index map for a PPSNR of 24 dB. The diagonal profiles for this section for and
are compared in Figures and . The reconstruction errors have been significantly reduced. Similar results are obtained for the other sections. The evolution of the partial data term
corresponding to the selected projection angle
is displayed as a function of the number of iterations in Figure without noise and for a PPSNR of 24 dB. Let
and
be the real and imaginary parts of the complex refractive index to be retrieved. In order to have a more quantitative information about the convergence of the method, the evolution of the normalized least square error using the
norm
and
is displayed as a function of the number of iterations without noise on the projections and for a PPSNR of 24 dB on Figures and . The normalized mean square errors on the two components of the refractive index are much decreased. It should be noted that the method developed does not require any a priori knowledge about the complex refractive index of the material and a priori about the total projected thickness of the sample for each orientation contrary to the other approaches found in the literature.[Citation8, Citation11, Citation21, Citation22] Others regularization schemes and prior information may be introduced in the object domain.
6 Conclusion
In this work, we have investigated the new nonlinear inverse problem associated to the reconstruction of the real and complex part of the refractive index in phase contrast tomography. The functional properties of the direct operator are investigated. The regularization scheme investigated is a classical nonlinear Tikhonov scheme. Our iterative algorithm to minimize the regularization functional is based on the use of the adjoint of the Fréchet derivative of the intensity operator, a modified Landweber method and the Simultaneous Algebraic Reconstruction technique. The proposed regularization method is tested on a multimaterial phantom. The reconstruction errors are much decreased. The proposed nonlinear method will be compared with other algorithms in future studies. This method opens promising perspectives to process experimental data in various applications.
References
- Momose A, Takeda TT, Tai Y, Yoneyama A, Hirano K. Phase-contrast tomographic imaging using an X-ray interferometer. J. Synchrotron. Rad. 1998;5:309–314.
- Davis GR, Wong SL. X-ray microtomography of bones and teeth. Physiol. Meas. 1996;17:121–146.
- Salome M, Peyrin F, Cloetens P, Odet C, Laval-Jeantet AM, Baruchel J, Spanne P. A synchrotron radiation microtomography system for the analysis of trabecular bone samples. Med. Phys. 1999;26:2194–2204.
- Nuzzo S, Peyrin F, Cloetens P, Baruchel J, Boivin G. Quantification of the degree of mineralization of bone in three dimensions using synchrotron radiation microtomograph. Med. Phys. 2002;29:2672–2681.
- Bayat S, Apostol L, Boller E, Brochard T, Peyrin F. In vivo imaging of bone micro-architecture in mice with 3D synchrotron radiation micro-tomography. Nucl. Instrum. Methods Phys. Res. 2005;548:247–252.
- Chappard C, Basillais A, Benhamou L, Bonassie A, Bonnet N, Brunet-Imbault B, Peyrin F. Comparison of synchrotron radiation and conventional X-ray microcomputed tomography for assessing trabecular bone microarchitecture of human femoral heads. Med. Phys. 2006;33:3568–3577.
- Cloetens P, Ludwig PW, Baruchel J, Van Dyck D, Van Landuyt J, Guigay JP, Schlenker M. Holotomography: quantitative phase tomography with micrometer resolution using hard synchrotron radiation X rays. Appl. Phys. Lett. 1999;75:2912–2914.
- Langer M, Cloetens P, Pacureanu A, Peyrin F. X-ray in-line phase tomography of multimaterial objects. Opt. Lett. 2012;37:2151–2154.
- Cloetens P, Barrett R, Baruchel J, Guigay JP, Schlenker M. Phase objects in synchrotron radiation hard X-ray imaging. J. Phys. D. 1996;29:133–146.
- Langer M, Cloetens P, Guigay JP, Peyrin F. Quantitative comparison of direct phase retrieval algorithms in in-line phase tomography. Med. Phys. 2008;35:4556–4565.
- Langer M, Cloetens P, Peyrin F. Regularization of phase retrieval with phase attenuation duality prior for 3D holotomography. IEEE Trans. Image Process. 2010;19:2425–2436.
- Davidoiu V, Sixou B, Langer M, Peyrin F. Nonlinear iterative phase retrieval based on Fréchet derivative. Opt. Express. 2011;23:22809–22819.
- Guigay JP, Langer M, Boistel R, Cloetens P. A mixed contrast transfer and transport of intensity approach for phase retrieval in the Fresnel region. Opt. Lett. 2007;32:1617–1629.
- Gureyev TE, Nugent KA. Phase retrieval with the transport of intensity equation: orthogonal series solution for non uniform illumination. Opt. Commun. 1996;13:1670–1682.
- Gureyev TE. Composite techniques for phase retrieval in the Fresnel region. Opt. Commun. 2003;220:49–58.
- Nugent KA. Coherent methods in the X-rays science. Adv. Phys. 2010;59:1–99.
- Paganin D, Mayo SC, Gureyev TE, Miller PR, Wilkins SW. Simulateneous phase and amplitude extraction from a single defocused image of a homogeneous object. J. Microsc. 2002;206:33–40.
- Paganin D. Coherent X-ray optics. New York (NY): Oxford University Press; 2006.
- Moosmann J, Hofmann R, Bronnikov A, Baumbach T. Non-linear phase retrieval from single-distance radioagraph. Opt. Lett. 2010;18:25771–25785.
- Sixou B, Davidoiu V, Langer M, Peyrin F. Absorption and phase retrieval in phase contrast imaging with nonlinear Tikhonov regularization and joint sparsity constraint regularization. Inverse Prob. Imaging. 2013;7:267–282.
- Beltran MA, Paganin DM, Uesugi Kand Kitchen MJ. 2D and 3D X-ray phase retrieval of multi-material objects using a single defocus distance. Opt. Express. 2010;18:6423–6436.
- Beltran MA, Paganin DM, Siu KKW, Fouras A, Hooper SB, Reser DH, Kitchen MJ. Interface-specific X-ray phase retrieval tomography of complex biological organs. Phys. Med. Biol. 2011;56:7353–7369.
- Born M, Wolf E. Principles of optics. Cambridge: Cambridge University Press; 1997.
- Goodman JW. Intoduction fo Fourier optics. Greenwood Village (CO): Roberts; 2005.
- Natterer F. The mathematics of computerized tomography. New York (NY): Wiley; 1986.
- Engl H, Kunisch K, Neubauer A. Convergence rates for Tikhonov regularization of nonlinear ill-posed problem. Inverse Prob. 1989;5:523–540.
- Engl H, Hanke M, Neubauer A. Regularization of inverse problems and its applications. Dordrecht: Kluwer; 1996.
- Scherzer O, Grasmair M, Grossauer H, Haltmeier M, Lenzen F. Variational methods in imaging. New York (NY): Springer; 2008.
- Herman GT. Image reconstruction from projections: the fundamentals of computerized tomography. New York (NY): Academic Press; 1980.
- Kak AC, Slaney M. Principles of computerized tomographic imaging. New York (NY): IEEE Press; 1988.
- Daubechies I, Fornasier M, Loris I. Accelerated projected gradient method for linear inverse problems with sparsity constraints. JFAA. 2008;14:764–792.
- Hanke M, Neubauer A, Scherzer O. A convergence analysis of the Landweber iteration for nonlinear ill-posed problems. Numer. Math. 1995;72:21–37.