Automatic regularization for tomographic image reconstruction

Results in Applied Mathematics 6 (2020) 100088

Contents lists available at ScienceDirect

Results in Applied Mathematics journal homepage: www.elsevier.com/locate/results-in-applied-mathematics

Automatic regularization for tomographic image reconstruction ∗

Eduardo Miqueles , Patricio Guerrero Scientific Computing Group, Brazilian Synchrotron Light Laboratory, CNPEM, Brazil

article

info

Article history: Received 6 September 2019 Received in revised form 6 December 2019 Accepted 14 December 2019 Available online xxxx Keywords: Regularization Phase Tomography Synchrotron

a b s t r a c t The phase retrieval process of imaging a sample can be modeled as a simple convolution process. Sometimes, such a convolution depends on physical parameters of the sample which are difficult to estimate a priori. In this case, a blind choice for those parameters usually lead to wrong results, e.g., extracting information from the reconstructed images. In this manuscript, we propose a simple connection between phase-retrieval algorithms and optimization strategies, which lead us to ways of numerically determining the physical parameters. © 2019 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

1. Introduction Automatic segmentation is a well known tool for extracting information from images, and there are many commercial software devoted to this branch of image analysis. Even though it is easy to handle image segmentation using standard tools, there are several examples where such techniques fail. This is the case of soft-tissue samples exposed to a transmission tomographic system. Using standard reconstruction algorithms, and dealing with conventional sinograms for a parallel geometry, it is almost impossible to extract quantitative information, even though it is easy to visually separate phases in the reconstructed slice. Voxels from these regions contain very similar characteristics and usual segmentation methods fail. However, it was previously reported by some authors [1,2] that a regularization on the measured data noticeably improves the contrast in the final reconstruction. As an attempt to differentiate the phases, Paganin’s [3,4] filter can be applied to the raw data. The approach proposed by Paganin et al. to restore the phase shifting of the object assumes that the sample is placed at a certain distance from the source. Even then, there is a major difficulty in defining the parameters {δ, β} for some samples, which compose the refractive index n = 1 − δ + iβ and play a major role on the numerical scheme. A major assumption on phase coefficients {δ, β} is that δ is proportional to β for all projection images. Our contribution on this manuscript is to provide simple forms to obtain the ratio δ/β using a simple one-dimensional optimization routine. This is almost immediate after recognizing the approach proposed by Paganin as a Sobolev regularization [5]. In fact, in our proposed technique, the ratio δ/β comes from a approach similar to the one used in the L-curve [6]. As an applied physics problem, an analytic first order approximated solution of the transport intensity equation (TIE) is obtained under a strong physical assumption on coefficients {δ, β} - which determines the physical parameter ℓ - leading to the effective Paganin equation. From the applied mathematics point of view (ignoring the existence of TIE), we search for the best smooth approximated solution to the measured data. This leads to a regularized least squares approach with ∗ Correspondence to: The Brazilian Center for Research in Energy and Materials (CNPEM), Rua Giuseppe Maximo Scolfaro 10.000, Polo II de Alta Tecnologia, Campinas, São Paulo, Brazil. E-mail address: [email protected] (E. Miqueles). https://doi.org/10.1016/j.rinam.2019.100088 2590-0374/© 2019 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/ licenses/by/4.0/).

2

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

Fig. 1. Fluxogram for the effective Paganin equation.

analytic solution given by the Paganin equation if and only if the regularization parameter is identical to the physical constant ℓ. The regularization approach establishes a different interpretation of TIE, which can be used whenever the sample does not completely satisfy the physical assumption. In fact, for tolerable deviations from the physical hypothesis, the effective Paganin equation can still be used although with a unknown parameter ℓ. A criterion to search for the best ℓ∗ – using the Tikhonov model and the L-curve – is proposed in this manuscript. A complete scenario, containing this discussion, is presented in Fig. 1. The manuscript is organized as follows: Section 2 presents the mathematical description of the effective Paganin equation, Section 3 presents the regularization approach and a motivation of our algorithm obtained from small-angle X-ray scattering. Three different paths to the effective Paganin equation through regularized problems are also provided: (a) Projection regularization, (b) Sinogram regularization and (c) Reconstruction regularization. Section 4 presents the numerical experiments using the proposed algorithm and Section 5 the final concluding remarks. 2. Mathematical description Mathematically, an intensity I is measured at a given pixel y on the plane detector, placed after the sample; see Fig. 2. For a near contact plane at z = 0 the measure is approximately described by Beer–Lambert law. Along the raypath Ω and at plane z = d > 0 - the paraxial wavefront of the scalar electromagnetic field is described by the transport intensity equation (tie) div(I(y , z)∇φ (y , z)) =

∂I (y) ∂z

(1)

where φ indicates the phase and y is a point on the detector plane. A large class of methods for the phase-retrieval problem is summarized in the work of Burvall [1]. It is claimed that a general approximate solution of the tie (1) at the detector plane z = d is proportional to the thickness of the sample along the raypath Ω – say p = p(y) – and must satisfy the following functional relation (−ℓ∇ 2 + 1)e−p(y) = f (y), f (y) :=

I(y) I0 (y)

(2)

with ∇ 2 standing for the Laplacian operator and ℓ a constant depending on the phase coefficients {δ, β} (which compose the refractive index of the sample) and the distance d. Here, I0 is the incident wavefront at the sample. After taking the Fourier transform on the above equation, the approximate solution p becomes p(y) = Pℓ (f ) := − log(Kℓ ⋆ f )(y).

(3)

Operator Pℓ indicates the operator described by Paganin, for a fixed ℓ. In this work, we are particularly interested in the kernel Kℓ described in the frequency domain by

ˆ Kℓ (q) =

1 1 + 4π 2 ℓ∥q∥2

(4)

where ˆ· stands for Fourier transformation fˆ (q) =

∫

f (y)e−iq·y dy R2

(5)

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

3

Fig. 2. Paraxial propagation describing the transport of intensity of a scalar electromagnetic wavefield through a sample.

For example, in the case of an object with a single material, the approximate solution was given by Paganin [3] with ℓ = dδ/µ. Using d = 0 the physical parameter ℓ becomes zero, and (2) determines p as the Radon Transform R. Phase-retrieval discussion starts with the assumption that the integration operation R along the sample s – through the raypath Ω and hitting the detector at point y – can be approximated by a mean value c(y). That is e−R[s](y) ≈ e−c(y)T [s](y)

(6)

where T [s] denotes the support (thickness) of the sample s along the ray parameterized by y. Hence, the general strategy in the literature, as pointed out in the work of Burvall [1], is that p(y) = T [s](y) is the approximate solution of the tie (1) - considering c(y) = 1 - which satisfies the following approximation f (y) = e−R[s](y) ≈ e−p(y) = (Kℓ ⋆ f )(y)





Near contact z=0







Fresnel regime z=d

(7)



First equality on the above equation comes from the standard Beer–Lambert law. Hence, obtaining the map p is a straightforward numerical process by means of the Fast Fourier Transform. Here, the constant ℓ comes from a strong physical assumption on the sample, which either relates to a single or a double material. From (7), the approximation f ≈ e−p indicates that in some sense, f − e−p should be minimized in p. It is obvious, due to the Beer–Lambert law for absorption, that the solution of f = e−p return p as the Radon transform of the attenuation coefficient, which also implies that we are not restoring the phase φ at the plane z = 0. Lemma 1.

The Paganin solution p = Pℓ (f ) - defined in (3) - is positive for all values of ℓ > 0.

Proof. We use the fact that Kℓ and f are positive functions and so is the convolution Kℓ ⋆ f . Since f (y) ≤ 1 - see (2) - and

ˆ Kℓ (0) = 1 we immediately obtain ∫ (Kℓ ⋆ f )(y) = Kℓ (u)f (u − y)du R2 ∫ Kℓ (u)du ≤ max f (y) y ∈R2

(8) (9)

R2

= max f (y) ˆ Kℓ (0) = 1

(10)

y ∈R2

Hence, p(y) = − log(Kℓ ⋆ f )(y) is positive for all values of ℓ.

■

3. Tikhonov regularization versus Paganin equation Before discussing the connection of the regularization approach with the Paganin equation, we briefly review a similar algorithm due to Glatter [7], very popular and extremely used in the small-angle X-ray scattering technique (saxs). Motivation from saxs: A complete description of saxs can be found in [7] and the references therein. The scattering intensity curve Iexp (q) is obtained and compared with the mathematical model I(q) described by the Fourier transform of the distribution function p of the particle under investigation, i.e., I = pˆ . After an orthogonal projection of p on the subspace span{ψ1 , . . . , ψN } (B-splines in the case of [7]), the mathematical problems rely on the determination of the coefficients c = (ck ) ∈ RN in the regularized least squares sense c ∗λ = argmin ∥Iexp − c ∈RN

N ∑ k=1

ˆ k ∥22 + λ∥K c ∥22 ck ψ

4

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

with K a second order finite difference matrix and λ ∈ R+ a regularization parameter. An optimal parameter λ∗ is chosen as the inflexion point of the parametric curve (log ∥K c λ ∥2 , log λ). In fact, a compromise is considered between the jump of the intensity curve I(q) and the euclidean distance of I and Iexp . By jump we refer to the phase shifts of the curve I along each measured radial point q; smoothness on the scattering curve eventually indicate a change of particle geometry, which is undesired for certain experiments. 3.1. Projection regularization In the L2 (R2 ) space equipped with the norm ∥u∥2 = ⟨u, u⟩, we look for the best smooth approximation of the normalized map f . Smoothness in the solution u is evaluated by the use of a ‘‘roughness measurement’’ operator A. This leads to the problem of finding a global minimizer of the following functional

∥u − f ∥2 + L∥Au∥2

minimize

(11)

u∈L2 (R2 )

with a regularization factor L. Indeed, (11) is a standard Tikhonov regularization approach for the smoothed least squares solution u. A typical definition for the functional A, enforcing smoothness on the solution u, is the following

∥Au∥2 =

∫

∥∇ u(y)∥2 dy

(12)

D

The first order optimality condition for the unconstrained problem (11) is determined by the Euler–Lagrange equations, which in the case of (12), coincide with the following normal equations u − f + LA∗ Au = 0 ⇐⇒ (LA∗ A + 1)u = f

(13)

where A∗ is the hermitian-adjoint operator of A = ∇ . Due to the fact that ∇ ∗ = −div and div ∇ = ∇ 2 , the optimality condition (13) becomes (−L∇ 2 + 1)u = f

(14)

which is essentially the same equation as (2), obtained by approximate phase retrieval methods, with L = ℓ. In this case, it means that the optimal regularization parameter λ is the one that obeys some physical assumption. The problem with the above minimization approach is that u is in fact defined by u = e−p , which is not an L2 function. Moreover, u = e−p is bounded within the interval [0, 1] because p is always a positive L2 function. Hence, we have to restrict our optimization problem with the constraint {u ∈ L2 (R2 ): 0 < u(y) < 1}, making the problem more difficult to solve. Therefore, we propose an optimization directly on the L2 space using the variable p through the functional H(p) := ∥e−p − f ∥2 + ℓ

∫

∥∇ e−p(y) ∥2 dy

(15)

R2

and with the additional constraint set defined by C = {p ∈ L2 (R2 ): p(y) ≥ 0}

(16)

which is closed and convex. We note that (a) Paganin solution p¯ = Pℓ (f ) is a feasible point due to Lemma 1, i.e p ∈ C ; (b) The tangential cone at p¯ - denoted as TC (p¯ ) (see [8]) - is identical to L2 since p¯ is an interior point; (c) A second order optimality condition holds due to the convexity of C , [8]: to prove that p¯ is a local minimizer, we have to show that the second Fréchet derivative H ′′ (p¯ ), acting as a bilinear form, satisfies H ′′ (p¯ )(v, v ) ≥ α∥v∥2 , ∀v ∈ TC (p¯ )

(17)

Proposition 1. The in-line phase-Retrieval approximated solution p = p(y) (a projection image) proposed by Paganin et al.– described in (3) and with initial data f = f (y) – is a local minimizer of the regularized constrained problem in L2 (R2 ) (p):

minimize H(p) := ∥e−p − f ∥2 + ℓ p∈C

∫

∥∇ e−p(y) ∥2 dy

(18)

R2

for a fixed parameter ℓ. The analytical solution of problem (p) is determined by (3). Proof. Let us denote the function to be minimized by H(p) = E(p) + ℓR(p), with E(p) = ∥e−p − f ∥2 ,

∫ R(p) = R2

 ( −p(y) )2 ∇ e  dy

(19)

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

5

The Fréchet derivative of nonlinear functionals R and E are presented in the Appendix, using f = e−g . The action of the linear functional H ′ (p) on each h ∈ L2 is determined by H ′ (p)h = 2

∫ =2

∫

e−p (e−p − e−g )h dy − 2ℓ

∫

e−p ∇ 2 e−p h dy

)

(

(20)

e−p [e−p − e−g − ℓ∇ 2 (e−p )]hdy

Hence, H ′ (p¯ ) = 0 if and only if p¯ satisfies the equation

(

¯

e−¯p − e−g − ℓ∇ 2 e−p¯

)

= 0,

(21)

which using u = e−¯p becomes the same as (14). That is, the solution of the Paganin equation is a stationary point of the functional H(p) = E(p) + ℓR(p). The second Fréchet derivative of each functional E and R are presented in Appendix. Form E ′′ (p) evaluated at a pair (v, v ) ∈ L2 × L2 obeys E (p¯ )(v, v ) = −2 ′′

∫

e−¯p [e−g − e−¯p − e−¯p ]v 2 dy

∫

∫ = −2 e−¯p [e−g − e−¯p ]v 2 dy + 2 e−2p¯ v 2 dy ∫ ∫ = 2ℓ e−¯p ∇ 2 (e−¯p )v 2 dy + 2 e−2p¯ v 2 dy

(22)

The second Fréchet derivative of R, at the stationary point p¯ follows R′′ (p¯ )(v, v ) = 2

∫

v e−¯p ∇ 2 (e−¯p v ) + v 2 e−¯p ∇ 2 (e−¯p )dy .

(23)

The bilinear form H ′′ (p) = E ′′ (p) + ℓR′′ (p) at p = p¯ , acting on (v, v ) – after several cumbersome but elementary simplifications – will be given by H (p¯ )(v, v ) = 4ℓ ′

+2ℓ

∫

∫ e

v ∇ (e )dy + 2

−¯p 2

2

−¯p

∫

v 2 e−2p¯ dy (24)

v e−¯p ∇ 2 (e−¯p v )dy

Since p¯ is bounded (due to the compactness of the object) function e−¯p and ∇ 2 (e−¯p ) will also be minorized by some constant, due to the fact that e−¯p ≥ M1 . Let us assume that ∇ 2 (e−¯p ) ≥ M2 . Using this fact, we can minorize (24) H ′ (p¯ )(v, v )∫ ≥ 4ℓM1 M2 ∥v∥2 + 2M12 ∥v∥2

+2ℓ

(25)

v e−¯p ∇ 2 (e−¯p v )dy

To minorize the last integral, we use the fact for all v such that ∥v∥ ≤ ε , v∇ 2 v ∼ 2v 2 and v∇v ∼ v 2 1 with 1 = (1, 1) ∈ R2 . Hence, −¯p v∇ 2 (e = v 2 ∇ 2 (e−¯p ) + 2∇ (e−¯p ) ]· (v∇v ) + e−¯p (v∇ 2 v ) [ 2v )−¯ ∼ ∇ (e p ) + 2∇ (e−¯p ) · 1 + 2e−¯p v 2

Expression between brackets in the above equation has a lower bound M3 since p has bounded variation. Now, (25) becomes H ′ (p¯ )(v, v ) ≥ (4ℓM1 M2 + 2M12 + 2ℓM3 M1 )∥v∥2

(26)

Since H is twice continuously Fréchet-differentiable, the Taylor expansion around p¯ gives 1 ′′ H (p¯ )(p − p¯ , p − p¯ ) + o(∥p − p¯ ∥2 ) 2 for all p such that ∥p − p¯ ∥ < ε . Using (26) and H ′ (p¯ ) = 0 we finally obtain p¯ as the local minimizer of H. H(p) = H(p¯ ) + H ′ (p¯ )(p − p¯ ) +

■

Notice that, as ℓ increases, the retrieved phase p becomes a blurred version of the normalized image f ; on the other hand, taking ℓ = 0 gives us the linear attenuation coefficient. Assuming that approximate phase-retrieval methods give a solution of the optimization problem (11), there are many numerical strategies to find an ideal parameter [2,9] ℓ, e.g., the L-curve [10]. Even though ℓ is strictly related to physical properties of the sample under investigation, the regularized solution (3)pℓ is the one maximizing the contrast of the solution. In fact, from the optimality condition (14), the following decreasing function

∥∇ 2 (e−pℓ )∥2 =

∥e−pℓ − f ∥2 ℓ2

(27)

always attains a maximum curvature at a critical point ℓ > 0. ∗

6

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

Lemma 2. Quantity ∥e−pℓ − f ∥2 /ℓ2 is a decreasing function of ℓ. Proof. Assuming pℓ = − ln(Kℓ ⋆ f ) with Kℓ given by (4), it is easy to note that e−pℓ − f = Kℓ ⋆ f − f . Using Parseval identity we obtain

∥Kℓ ⋆ f − f ∥2 ∥e−pℓ − f ∥2 = 2 2 ℓ ∫ ℓ 2 2 2 1 2 16π ∥q∥ ℓ = 2 |fˆ (q)| dq ℓ R2 [1 + 4πℓ∥q∥2 ]2

(28) (29)

as a monotonic function of ℓ since function 16π 2 ∥q∥2 /[1 + 4πℓ∥q∥2 ]2 is a decreasing function of ℓ, for all points q ∈ R2 . ■ Considering L∗ as the point of maximum curvature of quantity ∥e−pℓ − f ∥2 /ℓ2 defined in (27) is in agreement with the L-curve [6] criterion, which states that a compromise must be taken between the minimization of ∥u − f ∥ and smoothness by operator A through the plot of the parametric curve (∥uℓ − f ∥, ∥Auℓ ∥). Hence, we consider ℓ = ℓ∗ as the point satisfying

ℓ = argmax κf (L);

(30)

L≥0

and κf (ℓ) being the curvature of the function ln ξf (ℓ) with ξf = ∥∇ 2 pℓ ∥2∗ , defined by

κf (ℓ) =

|ξf′′ (ℓ)ξf (ℓ) − (ξf (ℓ)′ )2 | √ [ ′ ]2 3 ξf (ℓ)2 ξ (ℓ) 1 + ξf (ℓ) f 1

(31)

with prime denoting derivative with respect to ℓ (subscript f stands for frame since this is a process done by projections). Computing first and second derivatives from function ξf is a straightforward process due to properties of the Fourier transform. In fact, since e−pℓ = Gℓ with Gℓ = Kℓ ⋆ f from (3), it is easy to obtain n e−pℓ n!(−1)n (4π 2 ∥q∥)2n ∂ˆ ∂n ˆ (q) = e−pℓ (q) = fˆ (q) n n ∂ℓ ∂ℓ (1 + 4π 2 ℓ∥q∥2 )n+1

(32)

which implies that

∂ n [e−pℓ (y) ] = (Kℓ(n) ⋆ f )(y), n = 0, 1, 2 ∂ℓn (1)

(33)

(2)

with kernels {Kℓ , Kℓ } defined in the frequency domain by Eq. (32). Computing the optimal ℓ with an optimization algorithm for one-dimensional functions (such as Newton’s method or bfgs) is fast since the evaluation of the curvature function κ depends only on Fourier transforms. One of the main disadvantages of the Paganin filtering approach, is that each frame has to be processed individually in order to obtain a single filtered block of the raw data. From the new dataset, we extract two or three dimensional reconstructions using some inversion scheme. This is an extremely intensive computational process that can be done with a graphics processing unity to speed up calculations [11]. To avoid nested loops in a programming strategy of in-line phase retrieval method, we propose a filtering strategy directly on the sinogram. In fact, we can use the Lipschitz property of the exponential function

|e−a − e−b | ≤ |a − b|,

(34)

In our case, the unknown p approximates the path integral of the sample – say g = − ln f – (i.e., the Radon transform of µ) on pixel y of the area detector. Therefore, minimization of ∥p − g ∥∗ implies minimization of ∥e−p − e−g ∥ due to (34). Therefore, minimization of ∥p − g ∥2 implies minimization of ∥e−p − e−g ∥2 . This is true because the Fréchet derivative of the function E(p) = ∥e−p − e−g ∥2 is the linear functional E ′ (p): L2 → R determined by 1 2

′

∫

e−p(y) e−g(y) − e−p(y) h(y)dy

(

E (p)h =

)

(35)

y ∈D

Hence, p = g is a critical point of the functional E since E ′ (g) = 0 (first order optimality condition). Since E is convex, p = g must be a global minimizer. A proof of (35) is given in the Appendix. 3.2. Sinogram regularization Another strategy for phase-retrieval, the so-called Fourier method described in [12], claims that a solution p is obtained using p = Tℓ ∗ [− ln f ]

(36)

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

7

Fig. 3. Filtering action by frames (orthoslice F ) or slices (orthoslice S) on the cubic dataset.

with Tℓ a suitable kernel described in [1]. The main difference of (36) and (3) is the filter action, which is performed after the log-normalization. [ In the case ]−1of special Fresnel numbers (lower than one), Fourier methods like Rytov are given by (36), with Tℓ (q) = sin λπ d∥q∥2 . Assuming a cubic dataset with N slices (y2 axis), N rays (y1 axis) and N angles (θ axis), it is easy to note that in order to apply a phase-retrieval strategy by frames on the entire cube, the computational complexity will be bounded by O(N 3 log N), as illustrated by Fig. 3. The disadvantage of using the frame-strategy is the need to restore the cubic dataset prior to reconstruction. A phase-retrieval strategy, by slices, would be easily included in the stacking reconstruction process. In this sense, we look for a smooth sinogram function p which is close to measured sinogram g, i.e., a solution of the following regularized unconstrained optimization problem minimize V (p) := ∥p − g ∥2 + ℓ∥Y p∥2

(37)

p∈L2

with Y a smoothing operator. The Euler–Lagrange equations for the above functional V can be simplified considering the skew-adjoint operator Y = ∂∂ (here, t coincide with the y1 axis, the horizontal axis for a given slice) t

(ℓY ∗ Y + 1)p = g ⇐⇒

( −ℓ

) ∂2 + 1 p=g ∂t2

(38)

After taking Fourier transformations, we arrive at the optimal regularized solution pˆ (σ , θ ) =

(

)

1 1 + 4π 2 ℓσ 2

ˆ σ , θ ) = Tˆℓ (σ )g( ˆ σ , θ) g(

(39)

which is a filtering operation over the sinogram g, i.e., p(t , θ ) = (Tℓ ⋆ g)(t , θ ). The retrieved sinogram p obtained with (39) is similar to the in-line phase retrieval tomography in the Fresnel regime described in [1]. Both assume that phase coefficients {δ, β} are proportional, and also that β ≈ 0 which is true for the fish egg sample. Proposition 2. The phase-retrieval approximated solution p = p(t , θ ) (a sinogram image) proposed by the Fourier method – described in (36), (39) and with initial data g = g(t , θ ) – is a local minimizer of the regularized problem (r):

minimize ∥p − g ∥ + ℓ 2

p∈C

π

∫ ∫ R

0

(

∂ p(t , θ ) ∂t

)2

dtdθ

(40)

Proof. This an immediate consequence of (37) using Y = ∂/∂ t, leading us to the normal equations (38) and thus to a stationary point p determined by (39). It is a local minimizer since the second Fréchet of the objective functional is a positive semi-definite bilinear form, with α = 1 − ℓ (see Eq. (17)). ■ In this particular case, the optimal parameter ℓ can be found as the one that maximizes the curvature κs of the function ln ξs (ℓ) with

ξs (ℓ) = ∥∂t2 pℓ ∥2 ,

pℓ (t) = (Tℓ ⋆ g)(t , θ )

(41)

where ∂t2 stands for the second derivative on the ray axis t (subscript s stands for slice since this is a process done by sinograms). Curvature function κs is the same defined in (31), replacing ξf by ξs . As stated previously, a one-dimensional optimization algorithm can easily take advantage of the fast evaluation of the objection function κs through fast Fourier transforms. For completeness, it is easy to note that

∂ n [pℓ (t , θ )] = (Tℓ(n) ⋆ g)(t , θ ), n = 0, 1, 2 ∂tn

(42)

8

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088 (n)

with Tℓ

being represented in the frequency domain through

ˆ (n)

Tℓ (σ , θ ) =

n!(−1)n (4π 2 σ 2 )2n

(43)

(1 + 4π 2 ℓσ 2 )n+1

Although they never provide identical filtered datasets, the difference is bounded, as we show in Proposition 3. The main difference for the in-line phase retrieval methods described by Proposition 2 relies on the evaluation of the difference

∆ℓ (y) = −(Tℓ ⋆ ln f )(y) + ln (Kℓ ⋆ f ) (y) = (Tℓ ⋆ [− ln f ])(y) + ln (Kℓ ⋆ f ) (y)

(44) (45)

We propose the following result for an upper bound estimate to function ∆. Proposition 3. The difference (44) obtained by in-line phase retrieval methods – described in Proposition 2 – is bounded by ∥ ln f ∥ for all ℓ > 0. For ℓ = 0 the sinograms are equal. 1+4π 2 ℓ Proof. Case ℓ = 0 is trivial since lim Kℓ ⋆ f = f = lim Tℓ ⋆ f

ℓ→0

ℓ→0

⇒ lim |∆(y)| = 0 ℓ→0

Since Kℓ and Tℓ are positive operators, it is true from the triangle inequality that

|∆ℓ (y)| ≤ |(Tℓ ⋆ ln f )(y)| + |ln (Kℓ ⋆ f ) (y)| ≤ (Tℓ ⋆ |ln f |)(y) − ln (Kℓ ⋆ f ) (y)

(46) (47)

Last inequality holds since f < 1 and Kℓ ⋆ f ≤ 1. Due to Jensen’s inequality [13] we can majorize the second convolution integral of (47) to obtain

|∆ℓ (y)| ≤ (Tℓ ⋆ |ln f |)(y) − (Kℓ ⋆ ln f ) (y)

(48)

≤ (Tℓ ⋆ |ln f |)(y) + (Kℓ ⋆ [− ln f ]) (y)

(49)

≤ ∥Tℓ ⋆ |ln f |∥ + ∥Kℓ ⋆ |ln f |∥

(50)

≤ ∥Tℓ ∥2 ∥ ln f ∥ + ∥Kℓ ∥2 ∥ ln f ∥

(51)

Finally, using Parseval’s identity, each L norm ∥Kˆ ℓ ∥ and ∥Tˆℓ ∥ is asymptotically bounded in variable ℓ by the factor because the maximum of Kˆ ℓ and Tˆℓ is one, in variables q and σ , respectively. In this case, 2

max |∆ℓ (y)| ≤ y

2∥ ln f ∥

(52)

1 + 4π 2 ℓ

is the final bound for the difference (44).

1 1+4π 2 ℓ

■

It is important to note that the one-dimensional convolution integral with kernel Tℓ applies over the y 1 axis, keeping y 2 fixed (the slice axis). On the other hand, the two-dimensional integral of the denominator – with kernel Kℓ – applies for the entire domain of y. To overcome the difficulty in comparing kernels Tℓ and Kℓ , let us introduce the operator ∇ϵ2

∇ϵ2 =

∂2 ∂2 +ϵ 2 2 ∂ y1 ∂ y2

(53)

It is obvious that ∇ϵ2 → ∇ 2 as ϵ → 1 and Fourier transforms results in F [∇ϵ2 f ] = −4π 2 [q21 + ϵ q22 ]F [f ]. Returning to starting equation (2) and replacing ∇ 2 by ∇ϵ2 we arrive at an equation similar to (3) with a kernel Kℓ,ϵ → Tℓ as ϵ → 0 defined in the frequency domain by

ˆ K ℓ,ϵ (q) =

1 1 + 4π 2 ℓ[q21 + ϵ q22 ]

(54)

The next Proposition provides a relation between our one-dimensional optimization strategy and the dispersion measure of the input data. The dispersion depends on the cutoff frequency σc defining the interval for integration of the Fourier transform. Proposition 4. The optimal curvature parameter ℓ, which maximize the curvature of function ξs defined in (41), is inversely proportional to the cutoff frequency σc .

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

(k)

(1)

Fig. 4. Curvature of function D(ℓ, σc ) - from Eq. (59) - using σc

9

≤ σc(2) ≤ σc(3) .

Proof. We begin with definition (41) and use the fact that ∂t2 pℓ = (g − pℓ )/ℓ for all θ ∈ [0, π] fixed, but arbitrary. Due to Parseval’s identity we obtain for large values of the cutoff frequency σc

ξs (ℓ) = = =

1

1

1

ℓ

ℓ

ℓ

∥pℓ − g ∥22 = 2

1

∫

ℓ2 1

ℓ2

σc

(

σc −σc

1 1 + 4π 2 ℓσ

−σc

∫

∥pˆℓ − gˆ ∥22 ∼ 2

4π 2 ℓσ 2

(

1 + 4π 2 ℓσ 2

ˆ σ ∗ )2 = 16π 4 g(

∫

σc −σc

ˆ σ ∗ )2 = 32π 4 g(

[

2

ˆ σ ) − g( ˆ σ) g( 2 ˆ σ) g(

)2

σc

∫

ˆ σ ))2 dσ (pˆℓ (σ ) − g(

)2

dσ

(56)

dσ

(57)

σ4 dσ (1 + 4π 2 ℓσ 2 )2

3σc + 8π 2 ℓσc2 32π 4 ℓ2 (4π 2 ℓσc2 + 1)



(58)

3 arctan 2π

(

−

(55)

−σc

√ )] σc

(59)

2(4π 2 ℓ)5/2





D(ℓ,σc )

ˆ σ ∗ )2 for some σ ∗ over the Last equality is obtained using a straightforward integral calculation and a mean value of g( interval (−σc , σc ). It is easy to realize that function D is asymptotically defined as D(ℓ, σc ) ∼ O(σc /ℓ) so that D(·, σc(1) ) ≥ D(·, σc(2) ) ≥ D(·, σc(3) ) (1)

(2)

(3)

for an increasing sequence of cutoff frequencies σc ≤ σc ≤ σc . From the above inequality and since D has an (k) asymptotic behavior with 1/ℓ, the maximum curvature point ℓ∗k from function D(·, σc ) determines an increasing sequence ∗ ∗ ∗ ℓ1 ≤ ℓ2 ≤ ℓ3 , as shown in Fig. 4. ■ A similar discussion applies for the two-dimensional case, changing kernel Tℓ by Kℓ . in the 2D frequency variable q. Further strategies could be used, replacing operator Y in (37) by a generalized linear combination of first and second 2 derivatives on the t-axis in order to obtain smoother sinograms, i.e., Y = a1 ∂∂t + a2 ∂∂t 2 with adjoint given by Y ∗ =

−a1 ∂∂t + a2 ∂∂t 2 . In this case, the resulting differential equation, with his associated Fourier representation is ( ) ˆ σ , θ) ∂2 ∂4 g( ℓ −a21 2 + a22 4 + 1 p = g ⇒ pˆ (σ , θ ) = 2 ∂t ∂t 1 + 4π ℓ(a21 + 4σ 2 a22 )σ 2 2

(60)

3.3. Regularized reconstruction The filtered backprojection inversion algorithm states that a given slice image s can be reconstructed from a measured sinogram g through the inversion formula s(x) = BF [g ](x), x ∈ D. Operator B is the so-called backprojection operator, which is the adjoint of the Radon transform R, i.e., the stacking operator over the family of X-rays passing through the ˆ ˆ σ , θ ), for all θ . The complex sample. Operator F is a low-pass filtering operator, acting on the variable t, i.e., F g(σ ) = |σ |g( computational part of the fbp algorithm is the computation of B for all pixels in the region of interest D. It is an algorithm

10

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

Fig. 5. One-dimensional kernels (a) Tℓ and (b) |σ |Tℓ for the determination of the filtered sinograms using a phase-retrieval filtered backprojection strategy.

with complexity O(N 3 ) with N being the number of pixels for the reconstructed image. It was recently shown [14,15] that B can be computed with complexity O(N 2 log N) using a polar representation. In this case, it is possible to prove [15] that a dual formulation of the Fourier-slice theorem is also valid; it was referred as a backprojection-slice-theorem (bst). The theorem states that the Fourier transform of backprojection of any sinogram image can be computed in polar coordinates of the frequency domain using

ˆ ˆ σ , θ ), B [s](σ cos θ , σ sin θ ) = σ −1 g( with (σ , θ ) ∈ R+ × [−π, π]. Further details about the above formula can be found in [14,15]. The bst formula can be used to obtain an analytical solution of the standard Tikhonov regularization problem minimize ∥RC s − g ∥2 + 4π 2 ℓ∥s∥2

(61)

s∈L2

with C the smoothing self-adjoint operator defined as f (x′ )

∫ C f (x) = R2

∥x′ − x∥

dx′

(62)

In (61), we are looking for the best smooth approximation to the measured data. After exploring the Euler–Lagrange equations for the above optimization problem, we obtain the following representation in the frequency domain for the solution s sˆ(σ cos θ , σ sin θ ) =

(

σ 1 + 4π 2 ℓ2 σ 2

)

ˆ σ , θ) g(

(63)

The proof of the above equation can be found with details in [15], although using C as an identity operator. Since C is selfadjoint, the bst formula with further properties of the Fourier transform gives us (63). We note that (63) is a regularized version of the Fourier-Slice-Theorem and can be used to obtain s explicitly through any gridding strategy [16]. Applying (63) in the Fourier representation of s(x) we finally obtain a new representation for the reconstructed image s,

∫

dσ

s(x) = R

π

∫

dθ

(

0

|σ | 1 + 4π 2 ℓσ 2

)

ˆ σ , θ )|σ |eiσ x·ξθ g(

(64)

Eq. (64) provides exactly the same reconstruction pattern as a typical filtered backprojection reconstruction algorithm, but with the phase-retrieval kernel Tℓ acting on the sinogram g. In fact, we can generalize our regularized strategy in the following representation sℓ (x) = BF 2 [Tℓ g ](x),

x∈D

(65)

Now, {sℓ } is a family of solutions of the optimization problem (61), depending on the regularization parameter ℓ. The filter function Tℓ is defined in (39). Our regularized solution (65) depends explicitly on the computation of the backprojection operator B and either the bst formula, or different strategies could be used. Fig. 5a presents the family filter {Tℓ } for three different values of ℓ, while 5b the lowpass action of the filter F Tℓ on the input sinogram g. 4. Experimental validation Besides representing a major source of animal protein in the Amazon forest, fishes play a critical role in this ecosystem, mainly due to the exceptionally large bowl of the river, which allow them to interact along regional space, across various trophic levels [17]. Freshwater fishes are represented by around 8,500 species (over 40% of all known species of fish), mostly in the vast river systems and tropical lakes [18]. The interaction between the fish fauna with aquatic

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

11

Fig. 6. Slice reconstruction of a fish egg sample. The zoomed regions the difference between phases of the sample. The reconstruction was obtained using the filtered backprojection algorithm on the measured data.

ecosystems and biota takes place through food interrelationships and effects on the chemical composition of the water and sediment. Freshwater fish of the Amazon synchronize their reproductive period with the period of high waters, however, studies presented in recent years in the Intergovernmental Panel on Climate Change have shown global climate change, which has greatly affected the reproduction of this group leading to a decline of fish stocks in the region as well as the associated biotopes and interdependence of these. Also, the group of fish works as an excellent environmental indicator showing responses of changing state of normalcy that are likely to be measured. Thus, it is desirable to develop a protocol for the acquisition of Prochilodus nigricans (popular known as Curimatã in native language) eggs screenshots which emphasize possible responses in the reproductive dynamics of Amazonian fish populations and possible changes due to environmental changes on a global scale. These images could be used in studies eventually leading to proposals that could even steer public policy in the period closed in order to contribute to the conservation of fish stocks and associated biodiversity. Since the sample is mainly composed by soft tissue determined by approximately one single material, all measurements were done using a phase contrast setup. In the present paper, the focus is on the mathematics of the phase contrast tomographic reconstruction, in order to enable the development of such a protocol in the future. To preliminarily test our techniques, we expose a stained [19] sample of fish eggs on the imaging beamline of the Brazilian Synchrotron Laboratory aiming to explore the morphology of the Curimatã sample using conventional micro-tomography. The curimatã fishegg sample was exposed to the imaging beamline from the Brazilian synchrotron light source, using only 200 angles, 300 ms of exposure time and a distance of 27 cm. Fig. 6 presents a slice reconstructed using the filtered backprojection algorithm. It is clear that there is a small absorption of the sample and low contrast. Hence, image segmentation is nearly impossible. Well known image segmentation techniques, commonly implemented in commercial software for 3D visualization and data analysis, grossly failed in this case of study. Region-based segmentation methods [20] use the pixel grayscale to separate the image into different regions and threshold method works well with a bimodal distribution of grayscales, which is clearly not seen in this case, see Fig. 7(a). Watershed, a more robust regionbased segmentation method, comes from the geological idea of valleys and ridges identification, which are related to the gray level of the image. In this specific case, since there is a high grayscale variation observed between neighboring pixels - see Fig. 7(b) - this methodology oversegment the image, creating many small different regions.

12

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

Fig. 7. (a) Histogram for the reconstructed image with normal transmission sinogram data, as shown in Fig. 6. (b) Intensity plot along a line segment through the reconstructed image.

Table 1 Comparison of obtained values for the parameter ℓ (in meters) at a fixed slice.

Frames Slice

Optimal ℓ (meters)

max|∆ℓ (y)|

0.000402175636558 0.000463889940688

0.11386469548 0.11555728758

y

The reconstruction using the filtered backprojection with a filtered sinogram - using kernel Tℓ is presented in Fig. 8, where now contrast is considerably higher if compared with the conventional result of Fig. 6. A plot of the horizontal mark of Fig. 6, as-well as the histogram of the reconstruction after Paganin filter, is presented in Fig. 10. Compared with the histogram of Fig. 7 we see a clear difference of edge enhancement on the reconstructed result. Image segmentation of the new image is now much easier, and is presented in Fig. 9. Here, as the image contrast was enhanced, threshold tool was enough to separate most of the regions (e.g. contour regions 1 and 2 of Fig. 9). However, other regions still present high differences in grayscale (e.g. region 3 of Fig. 9), and a combination of edge detection, based on the image gradient, and fill interior tools had to be applied. The optimal value for ℓ was obtained using the curvature strategy described in Section 3. Two filtered sinograms are depicted in Fig. 11 comparing the regularization strategy by frames and by slice. Although the sinograms are not equal, the critical values of ℓ obtained using the one-dimensional strategy are similar, as shown in Table 1. To validate the numerical curvature strategy, we expose the sample at 29 different distances {d1 , d2 , . . . , d29 }

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

13

Fig. 8. Slice reconstructed of the fish egg sample using the transmission expectation maximization algorithm on the filtered data, using the strategy described in (39) with kernel Tℓ .

Fig. 9. Segmentation of the reconstructed slice of Fig. 8 using the filtering strategy described in (39). Three different regions can be obtained using the segmentation process; see text for details.

14

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

Fig. 10. (a) Histogram for the reconstructed image, using the regularized sinogram, see Fig. 8; (b) Intensity plot along a line segment through the regularized image.

Fig. 11. Sinogram obtained with regularization by (a) frames and (b) slice.

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

15

Fig. 12. Projection at different distances. Higher distance of 282 mm at left, and 4 mm at right.

obtaining a normalized sequence of frames {fk }, fk (y) =

Ik (y) I0 (y)

,

4 mm ≤ dk ≤ 282 mm

A small distance d1 indicates a near-contact propagation while for long distances like d29 , the phase becomes more visible. Each fk satisfies the approximation (2) with a parameter Lk that depends linearly on the distance dk [1,21], i.e., (−Lk ∇ 2 + I)e−p = fk ,

Lk =

δ 4πβ dk , µ = µ λ

(66)

where λ is the wavelength. First and last frame are presented in Fig. 12, where edge enhancement is clearly visible for high distances in the Fresnel regime. According to our numerical scheme, we can restore an approximation ℓk of Lk computing the point of maximum curvature of the function ξf in (30), i.e,

ℓk = argmax κf (ℓ; fk ), ℓ≥0

The sequence {ℓk } as a function of {dk }, obtained for this experiment, is presented in Fig. 13a. Since ℓk is not linear with respect to dk , as predicted by theory, we investigate the effect of changing the cutoff-frequency σc defined as m (67) σc = 2∆ where ∆ is the sampling rate on the square domain for y. Here, m is a positive constant and the mesh for the twodimensional frequency variable q is such that ∥q∥∞ ≤ σc . In this case, the blinded restoration (i.e., without knowing the true value of L) using the optimization strategy presented in the plot of Fig. 13a was obtained using m = 1. In order to investigate the effect of variable m, which defines the cutoff-frequency (67), we investigate our approach using a one-dimensional example. ( ) (Let p(t) ) = rect(t) be the rectangular function representing the electronic density, i.e., p(t) = 21w [sign t + 12 − sign t − 12 ]. We can approximate p using the following approximation for the sign(·)

√

function, sn (t) = t / t 2 + 1n , i.e., p(t) ≈ pn (t) =

1 2w

[ ( sn

x+

1 2

)

( )] 1 − sn x − 2

(68)

We have used w = 300 and n = 104 as depicted in Fig. 14a; in this case, the propagation IL (t) is determined by

( IL (t) =

−L

) ∂2 + 1 e−pn (t) ∂t2

( ) = −ℓ[p′′n (t) + p′n (t)2 ] + 1 e−pn (t)

(69) (70)

with derivatives {pn , pn } known a priori. Function IL is shown in Fig. 14b using value L = dδ/µ = 0.0163522409163 with d = 500 × 103 mm, δ = 1.043 × 10−6 , β = 3.553 × 10−10 , λ = 1.4 × 10−10 and µ = 4πβ/λ. Let us denote the recovered ′′

′

16

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

Fig. 13. Wrong blinded restoration values of ℓ, for an arbitrary value of m = 1 (Eq. (67)). As a function of distance, ℓ is not linear, contradicting the theory. (a) Using real data with 29 distances and (b) a similar simulated one-dimensional example.

Fig. 14. One-dimensional simulated example. (a) Electronic density function; (b) Measured function − log IL with L being fixed.

signal as p(ℓ,m) (t). Some examples with different values of the pair (ℓ, m) are presented in Fig. 15. Columns (a), (b) and (c) are for constant values of ℓ = 0.00163522409163, ℓ = 0.0825788166274 and ℓ = 0.163522409163, respectively. Rows are varying from top to bottom with m = 0.15, m = 0.25 and m = 0.35. Taking the result of Fig. 15, we observe (due to the intermediate value Theorem) that exists a value of ℓ between 0.00163522409163 and 0.163522409163 that restore the signal p. However, the restored value for an arbitrary value of m possibly does not have any relation with the one obtained by the optimization strategy. Let {mj } be the sequence of values for m, which controls the cutoff-frequency (67). The sequence of restored values {ℓ(mj )} ≡ {ℓj } is presented in Fig. 16a, where it is clear that exist a value of m such that ℓj = L. Now, computing the right m for each distance value dk using a force brute algorithm gives us the sequence mk depicted in Fig. 15b. For an increasing value of distance, ℓk behaves linearly as predicted by theory and presented in Fig. 15 and √ 15d. It is important to note that mk versus dk behave as predicted in Proposition 4, i.e., ℓk ∼ 1/ mk . In the last experiment we expose a sample made of a thick kapton tape to the imaging beamline, using 201 angles and a distance of d = 220 × 103 mm, each with an exposure time of 300 ms. The idea was to test the one-dimensional optimization algorithm by frames. A given slice reconstructed with the expectation algorithm maximization is presented in Fig. 17. The contrast is considerably poor inside the sample due to low absorption of the beam on the sample. Applying the one-dimensional strategy by frames with m = 1, we find an optimal ℓ∗ = 0.0013871, resulting in the reconstruction shown in Fig. 18. According to the database from nist [22], since kapton is a polyimide film with that δ = 3.9 × 10−6 and β = 7.02 × 10−9 we obtain the theoretical L given by L = 0.0013617, which is satisfactory to the ℓ∗ obtained by our algorithm.

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

17

Fig. 15. Examples of restored signal p(ℓ,m) with different values of ℓ and m. See text for details.

5. Conclusion In this manuscript we have described a strong relation between phase retrieve methods and a Tikhonov regularization for the imaging of light samples in the Fresnel regime. Considered as a regularized variational problem in an appropriate functional space, we provide a simple one-dimensional algorithm, which aims to restore the physical parameter ℓ, usually given empirically, and needed to restore the phase in the near contact regime. An automatic algorithm capable to restore the right convolution parameter improves quality of image segmentation, which depends on the contrast of the reconstructed image. Other segmentation tools can be used to validate the chosen parameter, such as the collaborative segmentation package proposed in [23]; future developments will incorporate collaborative tools with the reconstruction scheme. The one-dimensional algorithm proposed in this work relates directly to the cutoff-frequency m, also related to the dispersion of the measured data. Even though ℓ and m are also related, we claim that exist an optimal m∗ for which the 1d algorithm provide a best approximation of the appropriate physical parameter ℓ∗ . An algorithm to find the pair (ℓ∗ , m∗ ) was out of the scope of this manuscript. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. CRediT authorship contribution statement Eduardo Miqueles: Conceptualization, Methodology, Software.

18

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

Fig. 16. (a) Relation between the physical parameter L and the sequence {ℓj } for different cutoff frequencies {mj }. (b) Optimal cutoff parameter versus distance between sample and detector. (c) Sequence of restored parameter {ℓk } matching with the original sequence {Lk }. (d) Error between restored sequence {ℓk } and {Lk }.

Acknowledgments We would like to thank Dr. Nathaly L. Archilha and Dr. Thiago Spina for many valuable suggestions regarding the segmentation process of the sample presented in this manuscript. Appendix

A.1. Fréchet derivatives

√

Let ∥ · ∥ = ⟨·, ·⟩ be the usual L2 norm. We compute the Fréchet derivative of the following nonlinear functionals E and R, acting on a function p = p(y) ∈ L2 with y ∈ D: E(p) = ∥e−p − e−g ∥2

(71)

and

∫ R(p) =

 ( −p(y) )2 ∇ e  dy

(72)

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

19

Fig. 17. Reconstructed slice of a thick kapton sample using the expectation maximization algorithm without regularization.

Fig. 18. Reconstructed slice of a thick kapton sample using the expectation maximization algorithm with regularization on the sinogram with kernel Tℓ .

• Function E: First Derivative: For all h ∈ L2 , the error ∆E(p) = E(p + h) − E(p) is given by −p −h 2 −p 2 −p −h −g −p −g ∆E(p) ∫ = ∥e e ∥ − ∥e ∥ − 2∫⟨e e , e ⟩ + 2⟨e , e ⟩

=

e−2p(y) [e−2h(y) − 1]dy − 2

e−p(y) e−g(y) [e−h(y) − 1]dy

Since e−h(y) = 1 − h(y) + o(∥h∥2 ) for y ∈ D we finally get

∆E(p) = 2

∫

e−p(y) [e−g(y) − e−p(y) ]h(y)dy + o(∥h∥2 )

(73)

20

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088

The definition of the Fréchet derivative E ′ (p) comes from (73) since

∆E(p) − E ′ (p)h = lim o(∥h∥2 ) = 0 h→0 h→0 ∥ h∥ lim

Second Derivative: Second Fréchet derivative is a bilinear operator E ′′ (p) acting on a pair (v, w ) ∈ L2 × L2 as E ′′ (p)(v, w ) = lim

R′ (p + t w )v − R′ (p)v t

t →0

(74)

Using the same reasoning, we can easily verify that E ′′ (p)(v, w ) = −2

∫

e−p(y) v (y)[e−p(y) − 2e−p(y) ]w (y)dy

The approximation e−t w = 1 − t w + o(t 2 ) was used to obtain the above formula.

• Function R: First Derivative: Let us consider F the following operator

∫

∥∇ u(y)∥2 dy ,

F (u) = D

acting on functions of bounded variation. The Fréchet derivative of F is given by F ′ (u) = −2∇ 2 u, i.e., the linear operator in the following sense F ′ (u)h = −2

∫

∇ 2 u(y)h(y)dy

Now, our functional R is related to F as

∫ R(p) =

∥∇ (e−p(y) )∥2 dy = F (e−p )

or R(L(u)) = F (u),

L(u) = − ln u

Hence, using the chain rule in the Fréchet sense, R′ (L(u)) ◦ L′ (u) = F ′ (u) ⇐⇒ R′ (− ln u) ◦ L′ (u) = −2∇ 2 u Since L is a also a composition, its easy to realize that L′ (u) = −1/u, therefore ′

(

R (− ln u) ◦

1 u

) I

= −2∇ 2 u

with I the identity operator. Replacing p = − ln u we finally obtain R′ (p) ◦ (ep I) = −2∇ 2 (e−p ) ⇒ R′ (p) = −2e−p ∇ 2 (e−p ) As a linear operator, the action of R′ (p) is given by R′ (p)h = −2

∫

e−p(y) ∇ 2 (e−p(y) )h(y)dy

(75)

Second Derivative: Using (75) and (74) we obtain R′′ (p)(v, w ) = 2

∫

v (y)e−p(y) ∇ 2 w(y) + w(y)e−p(y) ∇ 2 (e−p(y) )v (y) dy

This derivative was obtained using the approximation e−t w = 1 − t w + o(t 2 ) in the definition (74). References [1] Burvall Anna, Lundström Ulf, Takman Per AC, Larsson Daniel H, Hertz Hans M. Phase retrieval in X-ray phase-contrast imaging suitable for tomography. Opt Express 2011;19(11):10359–76. [2] Rositi Hugo, Frindel Carole, Wiart M, Langer M, Olivier Cécile, Peyrin Françoise, et al. Computer vision tools to optimize reconstruction parameters in X-ray in-line phase tomography. Phys Med Biol 2014;59(24):7767. [3] Paganin David, Mayo SC, Gureyev Tim E, Miller Peter R, Wilkins Steve W. Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object. J Microsc 2002;206(1):33–40.

E. Miqueles and P. Guerrero / Results in Applied Mathematics 6 (2020) 100088 [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23]

21

Paganin David. Coherent X-ray optics, no. 6. Oxford University Press on Demand; 2006. Jonas P, Louis AK. A Sobolev space analysis of linear regularization methods for ill-posed problems. J Inverse Ill-Posed Probl 2001;9(1):59–74. Hansen Per Christian. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev 1992;34(4):561–80. Glatter Otto. A new method for the evaluation of small-angle scattering data. J Appl Crystallogr 1977;10(5):415–21. Luenberger David G. Optimization by vector space methods. John Wiley & Sons; 1997. Helou Neto Elias Salomão, Simões Lucas Eduardo Azevedo. A new parameter choice method for inverse problems with poisson noise. In: E-proceedings; 2012. Johnston Peter R, Gulrajani Ramesh M. Selecting the corner in the L-curve approach to Tikhonov regularization. IEEE Trans Biomed Eng 2000;47(9):1293–6. Maia FRNC, MacDowell Alastair, Marchesini Stefano, Padmore Howard A, Parkinson Dula Y, Pien Jack, et al. Compressive phase contrast tomography. In: Image reconstruction from incomplete data VI, vol. 7800. International Society for Optics and Photonics; 2010, p. 78000F. Gureyev Timur E, Davis Timothy J, Pogany Andrew, Mayo Sheridan C, Wilkins Stephen W. Optical phase retrieval by use of first Born-and Rytov-type approximations. Appl Opt 2004;43(12):2418–30. Needham Tristan. A visual explanation of Jensen’s inequality. Amer Math Monthly 1993;100(8):768–71. Miqueles Eduardo X, Helou Elias S. Fast backprojection operator for synchrotron tomographic data. In: European consortium for mathematics in industry. Springer; 2014, p. 243–52. Miqueles Eduardo, Koshev Nikolay, Helou Elias S. A backprojection slice theorem for tomographic reconstruction. IEEE Trans Image Process 2017;27(2):894–906. Marone F, Stampanoni M. Regridding reconstruction algorithm for real-time tomographic imaging. J Synchrotron Radiat 2012;19(6):1029–37. Lowe-McConnell Rosemary H. Estudos ecológicos de comunidades de peixes tropicais. In: Estudos Ecológicos de comunidades de peixes tropicais. 1999 [in Portuguese]. Cohen Daniel M. How many recent fishes are there? California Academy of Sciences; 1970. e Silva Juliana Martins de S, Zanette Irene, Noël Peter B, Cardoso Mateus B, Kimm Melanie A, Pfeiffer Franz. Three-dimensional non-destructive soft-tissue visualization with X-ray staining micro-tomography. Sci Rep 2015;5:14088. Fu King-Sun, Mui JK. A survey on image segmentation. Pattern Recognit 1981;13(1):3–16. Weitkamp T, Haas D, Wegrzynek D, Rack A. Ankaphase: software for single-distance phase retrieval from inline X-ray phase-contrast radiographs. J Synchrotron Radiat 2011;18(4):617–29. Rumble Jr JR, Bickham DM, Powell CJ. The NIST x-ray photoelectron spectroscopy database. Surf Interface Anal 1992;19(1–12):241–6. Spina Thiago V, Stegmaier Johannes, Falcão Alexandre X, Meyerowitz Elliot, Cunha Alexandre. Segment3D: A web-based application for collaborative segmentation of 3D images used in the shoot apical meristem. In: 2018 IEEE 15th international symposium on biomedical imaging. IEEE; 2018, p. 391–5.