Constrained Reconstructions in X-ray Phase Contrast Imaging: Uniqueness, Stability and Algorithms

This chapter considers the inverse problem of X-ray phase contrast imaging (XPCI), as introduced in Chap.2. It is analyzed how physical a priori knowledge, e.g. of the approximate size of the imaged sample (support knowledge), affects the inverse problem: uniqueness and—for a linearized model—even well-posedness are shown to hold under support constraints, ensuring stability of reconstruction from real-world noisy data. In order to exploit these theoretical insights, regularized New-tonmethodsareproposedasaclassofreconstructionalgorithmsthatﬂexiblyincorpo-rateconstraintsandaccountfortheinherentnonlinearityofXPCI.AKaczmarz-typevariantoftheapproachisconsideredfor3Dimage-recoveryintomographicXPCI,whichremainsapplicableforlarge-scaledata.Therelevanceofconstraintsandthecapabilitiesoftheproposedalgorithmsaredemonstratedbynumericalreconstructionexamples.

1 Basic physical model of XPCI: incident plane waves scatter on the imaged object, that is parametrized by a spatially varying refractive index n. The resulting diffraction-pattern (hologram) is recorded in the optical near-field at some distance behind the sample

Physical Model and Preliminaries
The basic physical model of XPCI is detailed in Chap. 2 and summarized by Fig. 14.1: incident monochromatic X-rays, modeled by plane waves, are scattered by the imaged sample, that is parametrized by its spatially varying refractive index n(x, z) = 1 − δ(x, z) + iβ(x, z) (δ, β: refractive-and absorption decrement). By the scatteringinteraction, a perturbation (the image) is imprinted upon the transmitted X-ray wavefield. The intensity I of the perturbed wave-field is recorded by a detector placed at a finite distance d > 0 behind the sample.
As derived in Sect. 2.1, the dependence of the hologram-intensities I from the sample-parameters δ and β is given by 2 for all x ∈ R 2 , The phase-and absorption-images φ and μ are 2D-projections of the 3Ddensities δ and β (k: X-ray wavenumber) along the incident z-direction. The Fresnelpropagator D, modeling free-space propagation of the X-rays between object and detector, is defined by Here, F( f )(ξ) := (2π) −m/2 R m exp(−iξ · x) f (x) dx is the Fourier-transform and f = kb 2 /d > 0 is the modified 1 Fresnel-number associated with the physical length b, that is identified with 1 in the chosen dimensionless coordinates.
XPCI-experiments provide intensity data I of the form (14.1) (up to data errors), whereas the images φ, μ are the quantities of interest. Hence, the following principal inverse problem has to be solved: Inverse Problem 1 (XPCI) For some set A, reconstruct a 2D-image h = μ + iφ ∈ A from measured holograms I of the form (14.1).
By rotating the object in Fig. 14.1, holograms I θ j may be acquired for different incident directions θ j ∈ S 2 = {x ∈ R 3 : |x| = 1} of the X-rays onto the sample (in Fig. 14.1, the incident direction coincides with the z-axis). This is the setting of Xray phase contrast tomography (XPCT). A mathematical model will be provided in Sect. 14.1.3. XPCT allows to probe 3D-variations of the parameters δ, β beyond mere projections φ, μ.

A Priori Constraints
The set of admissible images A in inverse Problems 1 and 2 is highly relevant. In order to facilitate and stabilize image reconstruction, the set A should be restricted as far as possible by available physical a priori knowledge: • Support constraints: real-world samples are of finite size. This implies that the functions f ∈ {φ, μ, δ, β} : R m → R have a compact support, i.e. are identically zero outside some bounded object-domain Ω ⊂ R m . • Non-negativity: by the physics of hard X-rays, the decrements δ, β-and thus also φ, μ-are always non-negative. • Pure phase object: especially for biological samples, β and μ are typically orders of magnitude smaller than δ and φ. Assuming a purely shifting-, i.e. non-absorbing object β, μ = 0, is then a good approximation. • Homogeneous objects: as is rigorously true for samples composed of a single material, proportionality of δ and β [φ and μ] may often be assumed. • Regularity: realistic images φ, μ, δ, β are not arbitrarily singular functions, but typically have some characteristic smoothness properties. • Tomographic consistency: Images φ and μ that arise as tomographic projections of one object under different incident directions are correlated.
Focussing on support-knowledge, we study the role of such constraints on inverse Problems 1 and 2 and outline how to exploit them algorithmically.

Additional Notation
We study inverse Problems 1 and 2 in spaces of square-integrable functions: 3) The focus lies on functions f ∈ L 2 (R m ) that have compact support supp( f ), i.e. that vanish outside some bounded domain Ω ⊂ R m : Furthermore, we define spaces of real-valued L 2 -functions: where Re(·), Im(·) denote the real-and imaginary parts, respectively.

Forward Operators for XPCI
Based on Sect. 14.1.1, we introduce forward maps F : X → Y modeling different settings of XPCI. Note that we define the maps in arbitrary dimensions m ∈ {1, 2, 3, . . .} although the natural case are images and holograms in m = 2 dimensions. The benefit of this will be seen in Sect. 14.3.4.1.

General Nonlinear Forward Operator
The most general (and most challenging) XPCI-setting is the reconstruction of both phase φ and absorption μ from a single hologram. According to (14.1), this setting is modeled by the forward map for complex-valued images h = μ + iφ. Note that the constant background intensity 1 has been subtracted, such that N (0) = 0. As a benefit, N can be analyzed as an operator on L 2 -spaces: for any bounded Ω ⊂ R m , is essentially 2 a well-defined, nonlinear operator. Moreover, it can be shown [1,2] that N is continuously Fréchet-differentiable, i.e. sufficiently smooth to admit local linear approximations. The derivative is given by

Linearized Forward Map and Contrast-Transfer-Functions
The nonlinearity of the forward map N causes difficulties in both analysis and practical image reconstruction. It is therefore standard [3][4][5][6][7] to resort to a linearization valid for weakly scattering samples (see e.g. [7] for details on the regime-ofvalidity): the idea is that the image f is sufficiently "small" so that higher-order terms are negligible: The linearized forward map T = N [0] is also known as the contrast-transferfunction-(CTF-)model, which refers to the following alternate form (compare with Sect. 2.2): According to (14.11), the linearized contrast in Fourier-space is given by a superposition of the Fourier-transforms of phase-and absorption-image φ, μ modulated by the oscillatory CTFs s 0 and c 0 , respectively.
Such a homogeneity-constraint may be incorporated into the general forward model, via a modified forward map 12) The linearized model under a homogeneity constraint may be expressed via a single CTF s ν (ξ) := sin |ξ| 2 /(2f) + ν : for ϕ ∈ L 2 (R m , R), it holds that Although (14.13) only holds for real-valued ϕ, we define S ν : L 2 (R m ) → L 2 (R m ) ( S ν = 2) on general L 2 -spaces. For its properties, it is widely irrelevant if real-or complex-valued functions are considered, as S ν commutes with the pointwise real-part:

Multiple Holograms
In order to obtain richer data in XPCI, it is standard to acquire multiple holograms I 1 , I 2 , . . . , I at several object-to-detector-distances, corresponding to different Fresnel-numbers f 1 , f 2 , . . . , f . This may be modeled by combining the forward maps for the individual holograms F j : ν } to a "vector-valued" operator: (14.14)

Forward Operators for XPCT
In X-ray phase contrast tomography (XPCT), holograms are measured under different incident directions θ ∈ S 2 . According to the basic model (14.1), the resulting intensities I θ are then given by where P θ is the parallel-beam projector along θ (θ ⊥ n x ⊥ n y ⊥ θ): According to the standard theory of computed tomography, projection-data {P θ ( f )} θ∈Θ for a suitable set of incident-directions Θ allows to reconstruct the underlying 3D-function f : R 3 → C. Analogously, the goal of XPCT is to reconstruct 3D-variations of the decrements δ and β of the sample's refractive index from a tomographic series of holograms {I θ } θ∈Θ .

Uniqueness Theory
In practice, it is highly relevant whether the measured intensity data I uniquely determines the sought image h = μ + iφ (or f = kβ + ikδ in XPCT). Otherwise, it might happen that two structurally different samples are indistinguishable by the imaging method, which is not desirable. (Non-)uniqueness of an inverse problem is equivalent to (non-)injectivity of the governing forward operator F : X → Y . Hence, it depends on different aspects: 1. The richness of the data, i.e. the size of the data-space Y : for example, it is commonly argued that measuring several holograms I 1 , I 2 , . . . at different Fresnelnumbers (see Sect. 14.1.2.4) helps to ensure uniqueness in XPCI. 2. Available a priori knowledge, i.e. the size of the object-space X : the smaller X the more likely it is that any two images In addition, it may happen that the nonlinear forward model is unique but its linearization is non-unique or vice verser. Accordingly, the different forward models from Sect. 14.1.2 have to be investigated individually.

Preliminary Results and Counter-Examples
We first review some known results on (non-)uniqueness of XPCI. Firstly, image reconstruction from a single hologram is generally non-unique: • Linearized model: • Nonlinear model (example from [8]): Images h ± : R 2 \ {0} → C; x → a(|x|) ± iν arctan2(x) for ν ∈ N and smooth functions a : R ≥0 → R give rise to so-called phase-vortices in the wave-field. The sign of the vortex is not determined by Fresnel-intensities (A := exp(−a)): |D(exp(−h + ))| 2 = |D(A · exp(−iν arctan2(·))| 2 = |D(A · exp(iν arctan2(·))| 2 = |D(exp(−h − ))| 2 (14.19) Based on these negative results, it is typically argued that at least two holograms and/or a homogeneity-constraint are required for uniqueness. Indeed, the situation improves substantially in the latter settings: • Uniqueness under homogeneity-constraints (linear): the operator S ν : L 2 (R m ) → L 2 (R m ) from Sect. 14.1.2.3 is injective, as the zero-manifolds of the Fouriermultiplier s ν are sets of the Lebesgue-measure 0 in R m . • Uniqueness for two holograms (linear): in [9], it is shown by a similar argument based on the CTF-representation (14.11) that also the operator Moreover, it is argued in [9] that both results carry over to the nonlinear model, provided that the image h is compactly supported. Indeed, a much stronger uniqueness result holds true under such an assumption, as will be seen in the following.

Sources of Non-uniqueness-The Phase Problem
According to the basic physical model (14.1), image-formation mathematically amounts to three operations: pointwise exponential, h → exp(−h), Fresnelpropagation, exp(−h) → D(exp(−h)), and computation of the pointwise squared modulus, D(exp(−h)) → |D(exp(−h))| 2 . Among those, D is an invertible operation, i.e. does not destroy information. This is not true for the other two operations, which give rise to different sources of non-uniqueness: • Phase-wrapping: The exponential is 2π-periodic in the imaginary-part of its argument. Hence, the phase-image φ = Im(h) may only be determined by the data up to increments by multiples of 2π. • Phase problem: The squared modulus, arising from the restriction of X-ray detectors to measuring intensities, eliminates the phase-information.
The first aspect is simpler to analyze and often turns out to be of lesser practical impact in XPCI: for moderately strongly scattering samples, φ is a priori known to assume values within [0; 2π), so that non-uniqueness due to phase-wrapping is not an issue. In the following, we therefore focus on possible ambiguities due to the phase problem.

Relation to Classical Phase Retrieval Problems
Up to possibly remaining phase-wrapping ambiguities, the image reconstruction problem in XPCI may be rephrased as follows: Such settings are known as phase retrieval problems as recovering O is equivalent to retrieving the missing phase of D(O) (and then inverting D). Uniqueness of phase retrieval has been extensively studied ever since the pioneering works of Walther [10] and Akutowicz [11,12], primarily for the case where D is replaced by the Fouriertransform F, i.e. for the reconstruction from phaseless Fourier-data. We refer to [13][14][15][16][17] for reviews.
Indeed, Fresnel-data may be readily reduced to the classical Fourier-setting, by rewriting the Fresnel-propagator in the form then the holograms in XPCI provide Fourier-data forÕ: Based on the identification in (14.21), uniqueness results for Fourier-phase retrieval may be adapted to the Fresnel-regime. Notably, however, most of such uniqueness theorems assume a compact support of the objective. Importantly, this is not justified in the setting of XPCI: The OTF O is not a compactly supported function in any realistic setting. Only the contrast o := O − 1 typically has compact support.

Holographic Nature of Phase Retrieval in XPCI
In order to emphasize the structural difference to classical phase retrieval problems, it is illustrative to rewrite the XPCI-model in the form where it has been used that D maps constant functions onto themselves. 3 According to the physical model from Sect. 14.1.1, the summands on the r.h.s. of (14.22) can be interpreted in terms of the scattered-and transmitted parts of the X-ray wave-field: the constant 1 is the intensity of the incident plane wave and the last summand that of the waves scattered by the object, whereas the second term describes the interference of these two wave-field components on the detector. Formula (14.22) places the inverse problem of XPCI in the realm of holographic phase retrieval problems, i.e. reconstruction in the presence of a reference signalhere provided by the unscattered part of the incident X-rays. Several theoretical and practical works have shown that such a holographic reference facilitates phase retrieval, see e.g. [18][19][20][21].

General Uniqueness Under Support Constraints
According to (14.22), image reconstruction in XPCI is equivalent to retrieving o = exp(−h) − 1 from data of the form (14.22) (up to possible phase-wrapping). By invertibility of the Fresnel-propagator D, uniqueness thus holds if it is possible to disentangle the summands D(o), D(o), and |D(o)| 2 . As shown in [22] using the theory of entire functions, the latter is indeed possible whenever o is known to have compact support, which is true for any sample of finite size. The principal result reads as follows: For any such K and Ω ⊂ R m bounded, Importantly, Theorem 14.1 establishes uniqueness in the most challenging setting of XPCI: single hologram, no homogeneity-constraint. The result trivially extends to every less difficult case with more data or additional constraints. However, note that the extension of uniqueness to restricted measurements I | K is based on analytic continuation of the data-a very unstable procedure in practice.

Uniqueness for the Linearized Model
Uniqueness for the linearized XPCI-model has to be shown individually. According to Sect. 14.1.2.2, it corresponds to data of the form (14.22), merely the quadratic term |D(o)| 2 is omitted and o = exp(−h) − 1 is replaced by −h (note that this rules out phase-wrapping!). Hence, the principal uniqueness argument from [22] remains valid: the summands D(h) and D(h) may be disentangled owing to their different "finger-prints" as entire functions: Corollary 14.1 (Uniqueness of linearized XPCI [22]) For any bounded domain Ω ⊂ R m and any K ⊂ R m that contains an open set, the linearized forward operator

Uniqueness for XPCT
By combining with standard results on uniqueness of tomographic reconstruction described by the theory of the Radon transform, the uniqueness theorems may be easily extended to XPCT. We refer to [22] for details.

Stability Theory
The uniqueness results of the preceding Sect. 14.2, suprisingly strong though they are, do not guarantee that accurate images may actually be reconstructed from holograms acquired in real-world XPCI-setups. Experimental data always contains errors due to noise and/or inaccuracies of the physical model. As detailed in Chap. 5 such data errors may lead to arbitrarily strongly corrupted images due to the phenomenon of ill-posedness: even if a forward model F : X → Y is injective, its inverse F −1 : F(X ) → X may be discontinuous such that small perturbations in the data g obs = F( f ) + may be arbitrarily amplified in the reconstruction The aim of this section is thus to supplement the uniqueness results with an analysis of stability, exploring how susceptible image reconstruction is to data errors. Thereby, it sheds a light on the question which reconstructions are feasible in practice. Due to difficulties arising from nonlinearity, the stability analysis is restricted to the linearized forward models.

Lipschitz-Stability and its Meaning
Although other (weaker) concepts of stability are common in the field of inverse problems, the notion of Lipschitz-stability turns out to be most suitable for XPCI: a forward map F : X → Y between normed spaces X, Y is said to be Lipschitz-stable if a stability estimate of the form holds for some constant C stab > 0. In this case, F has a Lipschitz-continuous inverse Notably, this implies robustness to data errors: given measurements g = F( f † ) + ∈ F(X ), the resulting reconstruction-error is bounded by The bound (14.25) states that data errors manifest at most amplified by a finite factor C −1 stab in the recovered object. Therefore C stab should be as large as possible: if C stab 1, the error-amplification predicted by (14.25) may be too large to guarantee accurate reconstructions at realistic noise-levels Y . Notably, for linear forward models F : X → Y , (14.24) is equivalent to

Stability for General Objects and one Hologram
Firstly, we consider the most challenging setting of reconstructing arbitrary phaseand absorption-images φ, μ from a single (linearized) hologram I ≈ 1 + T (h). Stable inversion of the forward map T is commonly argued to be infeasible. Indeed, as seen in Sect. 14.2.1, the forward model is not even unique for general images φ, μ ∈ L 2 (R m ), but only if φ, μ are compactly supported. Accordingly, we assume a support contraint in the following:

Analytical Approach
Our approach to analyzing stability is ultimately based on the principle of holographic reconstruction [23], that earned Dennis Gabor the Nobel Prize in physics in 1971.
The idea is to rewrite the forward map in the form This is Gabor's original idea of holographic reconstruction, which is illustrated in the first and second panel of Fig. 14.2 for a real-data example.
For stability analysis, we use the idea in a converse manner. By the constraint h ∈ L 2 (Ω), the sharp twin-image (the valuable part ins Gabor's eyes!) can be eliminated from (14.29) by restricting to the complement Ω c of Ω: By (14.30), we are left with incomplete (but phased!) Fresnel-data D 2 (h)| Ω c . Notably, to this point, only stable operations have been applied to the XPCI-data, which do not amplify data errors in L 2 -norm: for any g ∈ L 2 (R m ), it holds that D(g)| Ω c L 2 ≤ g L 2 . When applied to (14.30) this bound yields Finally, by employing the alternate form (14.20) of D, the Fresnel-data on the r.h.s. of (14.31) may be identified with incomplete Fourier-data:

Stability Bound
Since h L 2 = h L 2 , the bound (14.32) can be regarded as a relative stability estimate: recovering an image h ∈ L 2 (Ω) from XPCI-data T (h) is at least as stable as the reconstruction ofh ∈ L 2 (Ω) from Fourier-data outside the domain Ω f . Reconstruction from incomplete Fourier-data in turn is a well-studied problem: an uncertainty principle from [27] implies that Lipschitz-stability holds, F(h)| Ω c f ≥ C gen stab h for some C gen stab > 0, provided that Ω is bounded along at least one dimension. For rectangular domains Ω, the stability-constant C gen stab may be expressed in terms of the principal eigenvalue of a compact selfadjoint operator, for which asymptotics are derived in [28]. Via (14.32), these results yield stability estimates for linearized XPCI: While Theorem 14.2 only gives a worst-case bound on the data-contrast T (h) / h over all images h, the result may be sharpened considerably, as detailed in [24]: for any h ∈ L 2 (Ω), an individual lower bound for T (h) may be given based on the eigenvalues from [28] and the images that minimize T (h) / h may be characterized in terms of the associated eigenmodes.

Stability in a Practical Sense?
Numerical computations in [24] indicate that the bound (14.34) is quite sharp. While this is good news for a (pure) mathematician, it is bad news from an applied perspective: the predicted (quasi) exponential decay C gen stab (Ω, f) ∼ exp(f/8) implies that the constant quickly becomes very small for larger values of f, e.g. C gen stab (Ω, f) 10 −5 , for f 100. Notably, f = kb 2 /d is the modified Fresnel-number associated with the width of the support-domain Ω, i.e. with the diameter of the imaged sample. 4 In typical XPCI-experiments at synchrotrons, one has 10 2 f 10 5 , so that Theorem 14.2 only guarantees stability in practice for imaging settings at the lower end of typical Fresnel-numbers.
Notably, this is in line with empirics: after all, independent reconstruction of phase-and absorption-image φ and μ from a single hologram, as analyzed here, is widely considered as infeasible by practioners. It is thus highly surprising in the first place that the problem is technically well-posed at all.

Extension to Other Domains
Theorem 14.2 seemingly only applies to a very particular choice of the domain Ω ⊂ R m . Yet, it may be readily generalized via the following properties: • Translation-and rotation-invariance: As the map T is invariant under shifts and/or rotations of the coordinates, it holds that C
According to Sect. 14.2.1, uniqueness then also holds without support constraints, but image reconstruction is still ill-posed in general: the associated forward maps S ν : L 2 (R m ) → L 2 (R m ) and T (f 1 ,...,f ) : L 2 (R m ) → L 2 (R m ) do not have a bounded inverse due to zeros of the CTFs.
When both homogeneity-and support constraints can be assumed, well-posedness holds true with an improved stability constant compared to (14.34):

Theorem 14.3 (Stability estimate for homogeneous objects [24]) Let
for some constants c j > 0 that depend only on m.
A similar improvement applies for the reconstruction of general objects (no homogeneity-constraint) from two holograms: Theorem 14.4 (Stability estimate for two holograms [24]) Let Then In particular, for any support-domain Ω ⊂ R m , the following stability estimate holds true:

Order-Optimality
For ν = 0, it can be shown that the ∼ f −1 order of the decay in Theorem 14.3 cannot be improved: for a fixed bounded domain Ω ⊂ R m with non-empty interior, there exists a constant c max (Ω) > 0 such that This is a consequence of the bound S 0 (ϕ) Notably, better rates do not even hold in a setting with multiple holograms: for any Fresnel-numbers f 1 , . . . , f , it holds by a similar argument that (14.39) The reason for this surprising negative result on the benefit of multiple holograms is that the CTFs s (f i ) corresponds to the well-known low-frequency instability of XPCI that gives rise to the proven f −1 -rates of the stability constant.

Numerical Stability Computations
Other than for the setting in Theorem 14.2, the prediction (14.35) for the stability constant C hom stab (and thus for C two stab ) is far from sharp if the analytical bounds on the constants c j from [24] are inserted. Sharp values of C hom stab may however be computed numerically by approximating the minimum singular value of the operator S ν via techniques presented in [29,Sect. 3.4].

Phase Contrast Tomography
Although the physical setting of XPCI corresponds to m = 2 dimensions, the stability results in Theorems 14.2 to 14.4 have been formulated for arbitrary m. As a benefit, stability may be readily extended to XPCT: for the considered linearized forward models, XPCT data is of the form I θ − 1 = T (P θ ( f )) for T ∈ {T , S ν } and incident directions θ ∈ Θ ⊂ S 2 , compare Sect. 14.1.3. As noted in [30,31], the order of the projector P θ and T may be interchanged: (14.40) where } is the equivalent of T in m = 3 dimensions. As detailed in [29,Sect. 3.3], the relation (14.40) allows to express stability of linearized XPCT via known results for tomographic reconstruction, combined with stability bounds for T (3d) , S (3d) where Ω ⊂ R 3 . Stability then depends on a three-dimensional support constraint supp(β + iδ) ⊂ Ω ⊂ R 3 for the imaged sample's refractive index.

Imaging with Finite Detectors
There are a number of idealizing assumptions underlying to the obtained stability estimates: in addition to the neglected nonlinearity and idealizations in the basic physical model such as full coherence, it has also been assumed that the hologram I is measured in the whole detector-plane in Fig. 14.1. Due to the finite size of realworld detectors (and-more fundamentally-the finite width of illuminating X-ray beams), however, only restrictions I | K to some bounded domain K are available in practice.
According to Theorem 14.1, such restricted data has no impact on uniqueness (if K contains an open set). The situation is quite different in terms of stability, as analyzed in [32]: for any bounded K ⊂ R m -however large-the inverse problem of XPCI becomes severely ill-posed, i.e. Lipschitz-stability is lost so that data errors may severely corrupt the reconstructed images. Yet, it is also proven in [32] that the situation may be repaired by restricting to images h = μ + iφ of finite resolution (smoothness constraint in the sense of Sect. 14.1.1.1): by imposing that the h are B-splines on a Cartesian grid of sufficiently large spacing r (Ω, K , f) > 0 (i.e. pixelated images in some sense), Lipschitz-stability can be restored in the finite-detector setting. Physically, the necessity of such a restriction corresponds to a resolution limit that arises due to the finite numerical aperture associated with the detector size.

Regularized Newton Methods for XPCI
The following section considers regularized Newton-type methods for image reconstruction in XPCI. The proposed algorithm is motivated by the theoretical insights gained from Sects. 14.2 and 14.3.

Significance of Constraints
The stability results of Sect. 14.3 heavily rely on support constraints-without such, XPCI is ill-posed or even non-unique. To guarantee stability in practice, image reconstruction methods must thus be able to exploit support-knowledge. Also other types of a priori knowledge (see Sect. 14.1.1.1) are known to be beneficial. In particular, imposing non-negativity often has a similar stabilizing effect as support constraints.

Necessity of Iterative Methods
By far the most commonly used reconstruction method for XPCI at synchrotrons is direct CTF-inversion, as presented in Sect. 2.3. Within the notation of this chapter, the approach corresponds to quadratic Tikhonov regularization applied to the linearized forward maps S (f 1 ,...,f ) ν or T (f 1 ,...,f ) . Owing to the linearity and translationinvariance of these maps, the reconstruction may be implemented via a multiplication in Fourier-space (deconvolution), which renders the approach computationally fast.
However, direct CTF-inversion is incompatible with the above constraints: • Support constraints supp(h) ⊂ Ω for Ω R m break translation-invariance • Non-negativity is a nonlinear constraint: any reconstruction imposing it depends nonlinearly on the data I − 1-even for linear forward models!
In either case, reconstruction may thus no longer be achieved by deconvolution. Thus, iterative algorithms have to be applied to impose support-and/or nonnegativityconstraints in lack of efficient direct reconstruction formulas.

XPCI Beyond Linear Models
Although the linear CTF-model of XPCI has a surprisingly large regime-of-validity, there are settings where linear image reconstruction induces severe artifacts arising from the neglected nonlinearity, as demonstrated in Sect. 13.3. Reconstruction algorithms based on the full nonlinear XPCI-model are thus preferable in principle. The main obstacle in using such is that direct inversion formulas for the nonlinear model are not known. However, when iterative methods are needed anyway (Sect. 14.4.1.2), nonlinear forward maps cause little additional difficulty.

Reconstruction Method
In the following, we propose a reconstruction algorithm that meets the requirements discussed in Sect. 14.4.1. Details can be found in [33]. By choosing F : for Ω ⊂ R m , optional homogeneity-and/or support constraints are incorporated in the forward operator F. Consequently, such constraints are imposed automatically if image reconstruction in XPCI is performed by inverting F via any generic regularization method for inverse problems, see Chap. 5. In order to exploit Fréchet-differentiability of N (see Sect. 14.1.2.1) and the comparably moderate nonlinearity of XPCI, we choose regularized Newton methods as introduced in Chap. 5: for k = 0, 1, . . . , k stop , with initial guess h 0 ∈ X (usually h 0 = 0), observed (noisy) hologram(s) I obs and regularization parameters α k > 0. Note that we use a standard squared L 2 -norm as a data-fidelity term in (14.41), in lack of an accurate model for the data error statistics in flat-field corrected holograms.
s ) imposes tunable (by the choice of s ≥ 0) smoothness of the iterates h k and acts as a regularizer. Finally, R ≥0 (h, h k ) is a quadratic penalty term that is designed to correct negative values of Re(h k ) or Im(h k ) in the subsequent iterate h k+1 , see [33] for details.
In the numerical algorithm, a discretized analogue of the quadratic minimization problem in (14.41) is solved for images h * ∈ C N , data I obs ∈ R M and forward map F dis : C N → R M , via a conjugate-gradient method. The α k and k stop are chosen in a widely automated fashion, as detailed in [33].

Reconstruction Example
We assess the capabilities of the proposed method by reconstructing phase φ and absorption μ as independent parameters from a single simulated noisy hologram, which is shown in Fig. 14.3a. The considered test case is detailed in [33], where also a real-data example is considered for an analogous setting.
The true phase-image φ (Fig. 14.3b) is given by a bulk disk of magnitude 0.2, whereas the true absorption-image 0 ≤ μ ≤ 0.02 shows a logo-structure ( Fig. 14.3c). Accordingly, no homogeneity-constraint is applicable so that the test-case is situated in the most challenging, unstable setting of XPCI, which has been analyzed in Sect. 14.3.2. In particular, recall that image reconstruction is non-unique without exploiting further constraints.
The data is reconstructed using the regularized Newton method from Sect. 14.4.2, imposing non-negativity of φ and μ as well as support constraint, allowing nonzero values of φ, μ only within the circular region marked by the blue dashed line in Fig. 14.3b, c. The reconstructed images in Fig. 14.3d, e show that the proposed method correctly attributes the disk-structure to the phase-image φ and the logo-pattern to μ, without visible signs of "mixing things up". The overall lower reconstructionquality in μ compared to φ is due to the lower signal-to-noise in this parameter, as a realistically low absorption-refraction-ratio β/δ ≤ 0.1 has been assumed in the test case.
Now why does reconstruction of both φ and μ from a single hologram work here, contrary to the usual experience? The diameter of the circular support corresponds to a relatively low (modified) Fresnel-number f ≈ 87. According to the analysis in Sect. 14.3.2, this ensures stability of image reconstruction, as is discussed to greater detail in [33] and [24,Sect. 6]. By its ability to impose support constraints (and non-negativity), the proposed Newton-type method allows to exploit this theoretical stability in practice.

Regularized Newton-Kaczmarz-SART for XPCT
In the final section, we present a Newton-type reconstruction method for X-ray phase contrast tomography (XPCT) that is a compromise between flexibility w.r.t. a priori constraints and computational performance. We note that the method is an all-at-once approach, as also proposed in [30,31,34]: the 3D-object parameters δ, β are recovered directly from the full tomographic hologram-series, instead of first reconstructing 2D-images φ, μ for each hologram individually. Thereby, tomographic consistency is imposed as an additional constraint in image reconstruction, compare Sect. 14.1.1.1.
By replacing F ∈ {N , N ν } with the corresponding tomographic forward operator from Sect. 14.1.3, F PCT : X → L 2 (R 2 ) t ; f → (F(P θ j ( f ))) t j=1 with X = L 2 (Ω, (R)), Ω ⊂ R 3 , the regularized Newton method from Sect. 14.4.2 may be readily adapted to solve the inverse problem of XPCT: This is done in [2]. Yet, typical problem-sizes in XPCT with ∼ 10 9 dimensions of the discretized object-and data-space, are too large for this approach to be competitive in terms of computation times and memory requirements. As a remedy, we supplement the approach with a Kaczmarz-type strategy that exploits the block-structure of the XPCT-problem (14.42). The idea is to cyclically perform regularized Newton-steps w.r.t. the small sub-problems I obs θ j − 1 ≈ F(P θ j ( f )) defined by the measured holograms I obs θ j under the different tomographic incident directions θ j : for k = 0, 1, . . . , tn stop − 1 with n stop ∈ N. The parameters α > 0, 0 ≤ γ ≤ 1 control the regularization and smoothing w.r.t. the preceding iterate f k .
Iterations of the form (14.43) are known as regularized Newton-Kaczmarz [35]. The advantage compared to bulk (i.e. non-Kaczmarz-)methods is that the operatorblocks f → F(P θ j ( f )) require much less computations to evaluate than the total XPCT operator F PCT , which permits efficient computation of the iterates (14.41). Moreover, Kaczmarz-type methods often exhibit fast initial convergence, typically reaching a good reconstruction already after one or two cycles over the data, i.e. for n stop ∈ {1, 2}. To promote convergence, the processing order { j 1 , j 2 , . . .} ⊂ {1, . . . , t} of the data-blocks should be chosen such that subsequently fitted directions θ j k , θ j k+1 differ as strongly as possible, which we achieve by following a "multi-levelscheme" from [36].

Efficient Computation by Generalized SART
Although the processed data-size is reduced by the Kaczmarz-strategy, the iterates (14.41) still involve a minimization problem on a high-dimensional space of 3Dobjects f . Moreover, if the minimization is performed iteratively, each iteration requires evaluations of the (discretized) projector P θ j k and its adjoint P * θ j k , the back-projector, both of which typically amount to much higher computational costs than evaluating the XPCI forward map F. 5 Both computational issues can be resolved by computing the iterates (14.41) via a generalized SART 6 (GenSART-) scheme, as introduced in [38] for a much more general class of tomographic Kaczmarz-iterations: GenSART for Newton-Kaczmarz-iterations: 1. Forward-projection: p k := P θ j k ( f k ) 2. Optimization in projection-space (u j = P θ j (1 Ω )): The main benefit of the approach is that the required minimization is cast to projection-space, i.e. no longer needs to be solved on a high-dimensional space of 3D-objects but merely on 2D-images. Moreover, the whole scheme requires only a single evaluation of P θ j k (1.) and its adjoint P * θ j k (3.), whereas the optimization (2.) does not involve any of these costly operations anymore.

Parallelization and Large-Scale Implementation
Regularized Newton-Kaczmarz, computed via GenSART-schemes, is well-suited for large-scale computations and can be efficiently implemented in a parallelized manner. While we refer to [29,Sect. 6.3] for a detailed discussion, we mention the most important aspects here: • Low memory requirements: if the back-projection update (3.) (as well as the optional non-negativity projection) is implemented as an in-place operation, only a single 3D-array (storing f 0 , f 1 , ( f 1,≥0 , ) f 2 , . . .) needs to be kept in memory throughout the whole Newton-Kaczmarz-reconstruction. • Parallelized optimization: as the optimization-step (2.) works on 2D-images only, its memory-requirements are low enough to be performed on a single graphical processing unit (GPU) even for large-scale data. This permits efficient parallized implementation of this step. • Parallelized 3D-computations: The only operations on the 3D-objects f k are forward-and back-projections P θ j k , P * θ j k and pointwise arithmetics. All of these can be easily parallelized at low communication requirements between the different processors. In fact, it is possible to implement GenSART-schemes in a distributed manner: the object-iterates f k may be split into chunks, that are stored and managed by dedicated machines throughout the whole reconstruction. This property allows to run Newton-Kaczmarz reconstructions efficiently on multiple GPUs.
The biological sample constitutes a pure phase object to good approximation, i.e. vanishing absorption β = 0 may be assumed. Moreover, the sample is localized in a small subdomain of the imaged 2048 × 2048-sized field-of-view, as can be seen from Fig. 14.4a-c, i.e. support constraints may be imposed.
For comparison, we reconstruct the XPCT-data with different methods: • The additional constraints exploited in "Linear Kaczmarz" compared to "CTF+FBP" widely eliminate low-frequency background-artifacts (compare Fig. 14.4e-h) and thereby enable quantitatively correct reconstructions δ. • Though the sample-induced phase shifts are moderate, φ θ = kP θ (δ) 1, going over to the nonlinear XPCI-model has significant effects: especially in Fig. 14.4h, it can be seen that using the linearized model causes artificial distortions in the recovered object-density compared to the nonlinear Newton-Kaczmarz-reconstruction in Fig. 14.4i-l.
Accordingly, both the nonlinearity and the ability to exploit a priori constraints of the proposed Newton-Kaczmarz method turn out to be vital here to accurately reconstruct the anticipated 3D structure of the imaged bacteria 7 : cytoplasm with blob-shaped inclusions containing the DNA, where each of the two compounds is of approximately uniform density. The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.