X-ray holography, inverse problem, phase retrieval problem, uniqueness, intensity correlations

1 Introduction↩︎

Holographic X-ray phase contrast computed tomography has become a central tool in biomedical and material sciences [1]–[3]. It provides high-resolution 3D images of large volumes of quasi-transparent specimens such as biological tissues [4].

In this paper, we consider the 2D imaging model \[\label{eq2} I_{f}=|\mathcal{D}(e^{f}u)|^{2},\tag{1}\] where \(f:\mathbb{R}^2 \to \mathbb{C}\) is the unknown quantity of interest, \(u:\mathbb{R}^2 \to \mathbb{C}\) describes the incident beam, \(\mathcal{D}\) is a Fresnel-propagator (a unitary operator defined below), \(|\cdot|^{2}\) has to be understood point-wise, and \(I_{f}:\mathbb{R}^2\to [0,\infty)\) is the observed intensity. In contrast to the standard model, where \(u\) is a deterministic function describing a spatially perfectly coherent incident beam \(\tilde{u}(x_1,x_2,x_3)=u(x_1,x_2) e^{i\kappa x_3}\), we only assume partial coherence and model \(u\) as a complex, circularly symmetric, centered Gaussian process.

The values \(f(x_1,x_2)\) for any \(x_1,x_2\in\mathbb{R}\) are given by line integrals over the perturbation of the refractive index of the sample in the propagation direction \(x_3\) — \(\frac{1}{\kappa}\mathop{\mathrm{Im}}f\) corresponding to the integrated real part of the refractive index (phase contrast) and \(-\frac{1}{\kappa}\mathop{\mathrm{Re}}f\) to the imaginary part (absorption). Under the projection approximation, valid for optically thin samples, \(e^f u\) describes the electromagnetic field in a plane behind the sample. For \(\mathcal{D}=I\) (intensity measurements directly behind the sample) and for a plane incident wave (\(u\equiv 1\)), all information on the phase contrast is lost, and only the absorption contrast can be recovered. To retrieve the phase of \(e^f u\) and thereby phase contrast, the field is propagated to another parallel plane using the Fresnel propagator \(\mathcal{D}\), derived from the Helmholtz equation under the Fresnel approximation. This approach, known as propagation-based or inline holographic phase contrast tomography, can be combined with tomographic techniques to recover the full 3D refractive index of the sample.

Motivation↩︎

In practice, the coherence of an incident beam is never perfect but always partial, and the simplifying assumption of perfect coherence may lead to entirely erroneous or significantly suboptimal reconstruction results. An example where models based on perfect coherence completely fail is the self-amplified spontaneous emission (SASE) process generating X-ray free electron laser (XFEL) pulses [5]. Such ultrashort pulses can even enable time-resolved imaging at extremely small time scales. Another motivation arises from laboratory sources, which exhibit significantly poorer coherence properties compared to synchrotron radiation, yet are indispensable for high-throughput experiments.

If the incident field \(u\) is random, then obviously the measured intensities \(I_{f}\) in 1 are random as well. The primary objective of this work is to investigate and utilize the information contained in intensity correlations \[\begin{align} \label{eq:defi95two95point95corr} \operatorname{cov}[I_{f}](x,y) = \mathop{\mathrm{Cov}}(I_{f}(x),I_{f}(y)), \end{align}\tag{2}\] with a particular focus on the additional information beyond the mean intensities \(\mathop{\mathrm{\mathbb{E}}}[I_{f}(x)] =\sqrt{\operatorname{cov}[I_{f}](x,x)}\) (see eq. 35 ).

Let us first consider the situation that \(u=Z\widehat{u}\) with a complex-valued centered random variable \(Z\) with finite second moments (Gaussianity is not needed here), corresponding to perfect spatial coherence. Then \[\operatorname{cov}[I_{f}](x,y) = \operatorname{Var}(|Z|^2) \widehat{I}_{f}(x)\widehat{I}_{f}(y) =\sqrt{\operatorname{cov}[I_{f}](x,x)\operatorname{cov}[I_{f}](y,y)}, \quad \widehat{I}_{f}:=|\mathcal{D}(e^f\widehat{u})|^2,\] so in this case, intensity correlations and mean intensities uniquely determine each other and hence contain the same information.

To access the covariance function \(\operatorname{cov}[I_{f}]\) in 2 , we assume to have time resolved measurements of intensities \(I_{f,n}=|\mathcal{D}(e^fu_n)|^2\) corresponding to a sequence of fields \(u_1, u_2,\dots, u_N\) with the same distribution as \(u\). Such measurements are possible with modern detectors exploiting the high sensitivity and rapid readout of two-dimensional photon-counting pixel arrays. In other words, we partition the total photon dose into \(N\) fractions. We then have to impose the weak ergodicity assumptions \[\begin{align} \begin{aligned}\label{eq:ergodicity} \mathop{\mathrm{\mathbb{E}}}[I_{f}]&= \lim_{N\to\infty} \overline{I}_{f,N},&& \overline{I}_{f,N}:= \frac{1}{N}\sum_{n=1}^{N}I_{f,n},\\ \operatorname{cov}[I_{f}]&= \lim_{N\to\infty}\widehat{\operatorname{cov}}_{N}[I_{f}],&& \widehat{\operatorname{cov}}_{N}[I_{f}]:= \frac{1}{N}\sum_{n=1}^{N}(I_{f,n}-\overline{I}_{f,N}) (I_{f,n}-\overline{I}_{f,N})^\top. \end{aligned} \end{align}\tag{3}\]

Obviously, 3 is satisfied if the samples \(u_n\) are stochastically independent, which may reasonably be assumed for SASE pulses.

For a coherent beam, i.e., deterministic \(u\) in 1 , it is well known that the intensity \(I_{f}\) cannot be observed directly but only a Poisson process with intensity \(I_{f}\) describing photon counts. If \(u\) is a Gaussian process, then conditioned on \(u\) we also observe a Poisson process. The overall distribution of measured photon count data is then described by a Cox process, which will be considered in our numerical experiments, see 5.1. In such a model \(\operatorname{cov}[I_{f}]\) can be recovered in a joint limit where both \(N\) and the number of photon counts per frame tend to \(\infty\).

Contributions↩︎

On the theoretical side, we study the uniqueness of the inverse problem to reconstruct the projected refractive index \(f\) from intensity correlations. We identify three unavoidable sources of non-uniqueness: the first is a global phase shift, and the second are local phase shifts by integer multiples of \(2\pi\). The third is a bit less obvious and is caused by a twisted symmetry of \(f\). To avoid this last type of non-uniqueness, we impose a condition which, roughly speaking, requires the support of \(f\) to be contained in one half of the illuminated area. Then any two \(f_1, f_2\) which lead to identical noise-free intensity correlation data must satisfy \(e^{f_1-f_2} =e^{ic}\) for some real constant \(c\), i.e., \(f\) is identifiable up to the first two sources of non-uniqueness. We also prove a uniqueness theorem for the linearized problem.

The main challenge in the reconstruction process is the huge size of the correlations in high-resolution intensity maps. Note that \(\operatorname{cov}[I_{f}]\) depends on four variables! Such correlations are typically too large to be stored and too expensive to be computed in a pre-processing step. We propose a remedy for the case that the covariance operator of the incident beam \(u\) has low rank which avoids any matrices in the image space of the forward problem. The complexity of our algorithm grows linearly in the number of pixels of the intensity correlation and quadratically in the rank.

Related works↩︎

From the huge literature on phase retrieval problems, we only discuss some works related near-field holographic X-ray imaging as studied in this paper. For perfectly coherent sources (i.e., \(u\equiv 1\) in 1 ), the phase retrieval problem for general compactly supported samples admits a unique solution, provided at least two independent intensity patterns are measured at different object-to-detector distances [6].

The uniqueness of phase contrast imaging with a single intensity measurement has been proven for single material objects if weak absorption and slowly varying phase shifts [7] or small propagation distances, justifying the transport-to-intensity equation [8], [9] are assumed. Moreover, Nugent et al [10] by referring to the phase-vortex counterexample, discussed that two real-valued intensity measurements are not only sufficient but also necessary for unique reconstruction. However, Maretzke in [11] disproved the widely believed existence of ambiguities in single detector near-field phase contrast imaging and proved with the help of the theory of entire functions that a compactly supported complex-valued function can be uniquely determined from intensity measurements only at single detector distance and illumination wavelength. Stability estimates for the linearized problem were derived in [12].

Intensity correlation data were also studied in [13], [14] for a setup with unknown point scatterers where phases can be recovered in a preprocessing step.

Outline↩︎

The rest of the paper is organized as follows: 2 is devoted to the imaging model and the derivation of an explicit formula for the forward operator, for which we then compute the Fréchet derivative and its adjoint in 3. In 4 we discuss sources of non-uniqueness and prove the uniqueness result discussed above. Finally, our numerical algorithm dealing with a realistic noise model based on a Cox process is presented and tested on synthetic data in 5. Moreover, two short appendices are devoted to some supplementary issues.

2 sec:Forward32problem↩︎

We will derive our theoretical results in arbitrary space dimensions \(m\in \mathbb{N}\). The Fresnel transform \(\mathcal{D}:L^{2}(\mathop{\mathrm{\mathbb{R}}}^{m})\longrightarrow L^{2}(\mathop{\mathrm{\mathbb{R}}}^{m})\) can be defined by \[\label{Fresnel32op} (\mathcal{D}f)(x):=\int_{\mathop{\mathrm{\mathbb{R}}}^{m}}k_{\mathfrak f}(x-y)f(y)\,\mathrm{d}y,~~~\text{for all}~x\in\mathop{\mathrm{\mathbb{R}}}^{m}\tag{4}\] (see, e.g., [15]). Here the convolution kernel \(k_{\mathfrak f}:=c_m(\frac{\mathfrak f}{2\pi})^{m/2}n_{\mathfrak f}\) is given by the chirp function \[n_{\mathfrak f}(x):=e^{i\mathfrak f\frac{|x|^{2}}{2}}\] with the dimensionless Fresnel number \(\mathfrak{f}>0\) and the constant \(c_m:=e^{\frac{-im\pi}{4}}\). There exists the following alternative expression for the Fresnel propagator, which will prove valuable in the following: \[\label{Fresnel32op-Alter} (\mathcal{D}f)(x)=c_m\mathfrak{f}^{\frac{m}{2}}n_{\mathfrak{f}}(x)\cdot\mathcal{F}(n_{\mathfrak{f}}\cdot f)(\mathfrak{f}x).\tag{5}\] Here \(\mathcal{F}:L^{2}(\mathop{\mathrm{\mathbb{R}}}^{m})\longrightarrow L^{2}(\mathop{\mathrm{\mathbb{R}}}^{m})\) denotes the Fourier transform with the convention \[\label{FT} (\mathcal{F}f)(\xi):=\frac{1}{(2\pi)^{m/2}}\int_{\mathop{\mathrm{\mathbb{R}}}^{m}}f(x)e^{-i\xi\cdot x}\,\mathrm{d}x, ~~~~~\xi\in\mathop{\mathrm{\mathbb{R}}}^{m}.\tag{6}\] With this definition the Fourier transform is unitary (see, e.g., [16]), which also implies that \(\mathcal{D}\) is unitary.

We assume that the object plane (in dimension \(m=2\)) can be restricted to an open, bounded set \(\mathbb{D}\subset \mathbb{R}^m\). This may be the “illuminated area” \[\label{object32domain} \mathbb{D}:=\{x\in\mathbb{R}^m~:~\operatorname{cov}[u](x,x)>0\}.\tag{7}\] From 1 it is obvious that no information on \(f\) is available in the non-illuminated region \(\mathbb{R}^m \setminus \mathbb{D}\) where \(u=0\) almost surely. As mentioned in the introduction, \(u\) is assumed to be a circularly symmetric Gaussian process with a known covariance \(\operatorname{cov}[u]\). To mitigate the requirement that \(\operatorname{cov}[u]\) is known on \(\mathbb{D}\times \mathbb{D}\), a pinhole may be placed in the object plane. In such a setting \(\mathbb{D}\) describes the shape of the pinhole, which of course should be chosen large enough to contain the support of \(f\).

Recall that circular symmetry of \(u\) means that the distribution of \(e^{i\alpha}u\) is independent of \(\alpha\in\mathbb{R}\) and that the covariance operator \(\mathop{\mathrm{\mathbf{Cov}}}[u]\) of the random process \(u\) can be defined implicitly by \(\langle \mathop{\mathrm{\mathbf{Cov}}}[u]\varphi_{1},\varphi_{2}\rangle:=\mathop{\mathrm{Cov}}(\langle u,\varphi_{1}\rangle,\langle u,\varphi_{2}\rangle)\) for all \(\varphi_{1},\varphi_{2}\in L^2(\mathbb{D})\). For a function \(g\in L_{\mathbb{C}}^{\infty}(\mathbb{D})\) we define the multiplication operator \(M_{g}\in \mathcal{B}(L^2(\mathbb{D}))\) by \(M_{g}\varphi:=g \cdot \varphi\) for \(\varphi\in L^2(\mathbb{D})\) and recall the definition of the Fresnel propagator \(\mathcal{D}\) in 4 . It is straightforward to see that \(v_{f}:=\mathcal{D} M_{e^f}u\) is again a circularly symmetric Gaussian process with covariance operator \[\label{eq:32corr} \mathop{\mathrm{\mathbf{Cov}}}[v_{f}]:=\mathop{\mathrm{\mathbf{Cov}}}[\mathcal{D}M_{e^f}u]=\mathcal{D}M_{e^f} \mathop{\mathrm{\mathbf{Cov}}}[u] M^{\ast}_{e^f} \mathcal{D}^{\ast}.\tag{8}\]

Recall that a compact linear operator \(\mathcal{K}\) on a Hilbert space \(\mathbb{X}\) is called a Hilbert-Schmidt operator if its singular values \(\sigma_n(\mathcal{K})\), \(n\in\mathbb{N}\) are square summable and that the set \(\mathop{\mathrm{\mathcal{HS}}}(\mathbb{X})\) of Hilbert-Schmidt operators on \(\mathbb{X}\) equipped with the norm \(\|\mathcal{K}\|_{\mathop{\mathrm{\mathcal{HS}}}}^2=\sum_{n=1}^{\infty}\sigma_n(K)^2\) is a Hilbert space. For the special case \(\mathbb{X} = L^2(\mathbb{M})\) the kernel-to-operator map defined by \[\begin{align} \label{eq:32HSIO} \begin{aligned} &\mathop{\mathrm{KtO}}:L^{2}(\mathbb{M}\times\mathbb{M})\longrightarrow\mathop{\mathrm{\mathcal{HS}}}(L^{2}(\mathbb{M})),\\ &(\mathop{\mathrm{KtO}}[k]\varphi)(x):=\int_{\mathbb{D}}k(x,y)\varphi(y)\,\mathrm{d}y,~~ \text{for }~x\in\mathbb{M},~ \varphi\in L^2(\mathbb{M}),~k\in L^{2}(\mathbb{M}\times\mathbb{M}) \end{aligned} \end{align}\tag{9}\] is unitary ([17]), so its inverse, the operator-to-kernel map is given by \[\mathop{\mathrm{OtK}}:=\mathop{\mathrm{KtO}}^{-1}=\mathop{\mathrm{KtO}}^{\ast}:\mathop{\mathrm{\mathcal{HS}}}(L^2(\mathbb{M}))\to L^2(\mathbb{M}\times \mathbb{M}).\]

Using the identity \(\mathop{\mathrm{Cov}}(|X|^2,|Y|^2)=|\mathop{\mathrm{Cov}}(X,Y)|^2\) for circularly symmetric complex Gaussian variables \(X,Y\) (see 34 ) with \(X:=v_{f}(x)\) and \(Y:=v_{f}(y)\), we obtain the following relation between the two-point intensity correlations \(\operatorname{cov}[I_{f}]\) in 2 and the two-point correlations \(\operatorname{cov}[v_{f}]\) of the phased fields \(v_{f}\): \[\operatorname{cov}[I_{f}](x,y) = \big|\operatorname{cov}[v_{f}](x,y)\big|^2\quad \text{with}\quad\operatorname{cov}[v_{f}](x,y):=\mathop{\mathrm{Cov}}(v_{f}(x),v_{f}(y)).\] We also introduce an open, bounded subset \(\mathbb{M}\subset \mathbb{R}^m\) in the observation plane (again for \(m=2\)) in which intensity data are available. Using the previous identity and equations 1 , 2 , and 8 , the forward operator \[\tag{10} \begin{equation} F:L_{\mathbb{C}}^{\infty}(\mathbb{D})\longrightarrow L_{\mathbb{R}}^{1}(\mathbb{M}\times\mathbb{M}),~~~ F(f):=\operatorname{cov}[I_{f}] \end{equation} admits the following explicit representation: \begin{equation}\tag{11} F(f)= \big|\operatorname{cov}[v_{f}]\big|^2=\big|\mathop{\mathrm{OtK}}\mathop{\mathrm{\mathbf{Cov}}}[v_{f}]\big|^{2} = \big|\mathop{\mathrm{OtK}}(\mathcal{D}M_{e^f}{\boldsymbol{C}ov}[u]M^{\ast}_{e^f}\mathcal{D}^{\ast})\big|^{2}\,. \end{equation}\] We note that \(\mathop{\mathrm{\mathbf{Cov}}}[v_{f}]\) is a Hilbert-Schmidt operator (see [Adjoin-Frechet-diff], part (i) below) such that the expressions in 10 are well-defined.

We also note the following alternative expression for the forward operator \[\begin{align} \label{eq:Fourier-PRP} \begin{aligned} F(f)(x,y) &=|\mathcal{F}_{2m}(K_f)(\mathfrak{f}x,-\mathfrak{f}y)|^2,\quad x,y\in\mathbb{M}\\ &\text{with }K_f(p,q):= \mathfrak{f}^{m} n_{\mathfrak{f}}(p)e^{f(p)}K(p,q)e^{\overline{f(q)}}~\overline{n_{\mathfrak{f}}(q)} \end{aligned} \end{align}\tag{12}\] where \(\mathcal{F}_{2m}\) denotes the \(2m\)-dimensional Fourier transform. This follows from the identity \(\operatorname{cov}[v_{f}](x,y)=n_{\mathfrak{f}}(x)(\mathcal{F}_{2m}K_{f})(\mathfrak{f}x,-\mathfrak{f}y)\overline{n_{\mathfrak{f}}(y)}\), \(x,y\in\mathbb{M}\) obtained from the form 5 of the Fresnel propagator. This shows that our inverse problem can be seen as a Fourier phase retrieval problem in even dimensions with a rather special structure.

The complete nonlinear inverse problem can then be formulated as stably reconstructing \(f\) from sample correlations \(\widehat{\operatorname{cov}}_{N}[I_{f}]\) defined in 3 by approximately solving the equation \[\label{eq:nonlinear95inverse95problem} F(f) \approx \widehat{\operatorname{cov}}_{N}[I_{f}].\tag{13}\]

3 Fréchet derivative and its adjoint↩︎

3.1 Fréchet derivative of the forward operator↩︎

The following lemma on well-definedness and Fréchet differentiability of the forward operator is mostly straightforward and mainly formulated to introduce notation:

3.2 Adjoint of Fréchet derivative↩︎

For iterative regularization methods, we also need the adjoint of the Fréchet derivative. This is mostly standard, except for the adjoint of the mapping \[\label{eq:operatorM} \mathcal{M}:L_{\mathbb{C}}^{\infty}(\mathbb{D})\longrightarrow\mathcal{B}(L^{2}(\mathbb{D})), \qquad h\mapsto M_{h}\tag{14}\] which occurs in \(\mathcal{C}'[f]\). In the discrete setting, it is evident that the adjoint of the mapping \(\mathop{\mathrm{diag}}:\mathbb{C}^{m}\longrightarrow\mathbb{C}^{m\times m}\), defined by \(h\mapsto\mathop{\mathrm{diag}}(h)\) with respect to Euclidean or Frobenius inner products, is the function \(\mathop{\mathrm{Diag}}: \mathbb{C}^{m\times m}\to \mathbb{C}^m\) which maps a square matrix to its diagonal. Formulating a continuous analog is less obvious. Let \(\mathcal{K}\in\mathop{\mathrm{\mathcal{HS}}}(L^{2}(\mathbb{D}))\) be a Hilbert-Schmidt integral operator corresponding to kernel \(K\in L^{2}(\mathbb{D}\times\mathbb{D})\). We are motivated to define \((\mathop{\mathrm{Diag}}\mathcal{K})(x):=K(x,x)\) as the diagonal of the operator kernel \(K\). Since \(K\in L^{2}(\mathbb{D}\times\mathbb{D})\) and the diagonal set \(\lbrace (x,x)\in\mathbb{D}\times\mathbb{D}~:~x\in\mathbb{D}\rbrace\) have zero measure, the restriction of \(K\) to the diagonal is not well-defined. However, for the subspace \(\mathcal{S}_{1}(L^{2}(\mathbb{D}))\) of operators \(\mathcal{K} \in \mathop{\mathrm{\mathcal{HS}}}(L^2(\mathbb{D}))\) for which the singular values of which are not only square summable but even summable, equipped with the norm \(\|\mathcal{K}\|_{\mathcal{S}_1}:=\sum_{n=1}^{\infty}\sigma_n(\mathcal{K})\), two of the authors have recently shown in [18] that there exists a unique bounded linear operator \[\mathop{\mathrm{Diag}}:\mathcal{S}_{1}(L^{2}(\mathbb{D}))\longrightarrow L^{1}(\mathbb{D}) \quad \text{with}\quad \mathop{\mathrm{Diag}}(\mathcal{K})(x)=K(x,x)\] for all \(x\in\mathbb{D}\) and all operators \(\mathcal{K}\in\mathcal{S}_{1}(L^{2}(\mathbb{D}))\) with continuous kernel \(K\). Operators \(\mathcal{K}\in \mathcal{S}_1(L^2(D))\) are called trace class operators since the trace of \(\mathcal{K}\) is well-defined by \(\mathop{\mathrm{tr}}\mathcal{K}:=\sum_{n=1}^{\infty}(\mathcal{K}\varphi_n,\varphi_n)\) for any complete orthonormal system \(\{\varphi_n:n\in\mathbb{N}\}\) in \(L^2(\mathbb{D})\). Moreover, \(\mathop{\mathrm{tr}}(\mathcal{K})=\int_{\mathbb{D}}\mathop{\mathrm{Diag}}(\mathcal{K})\,\mathrm{d}x.\)

The adjoint of \(\mathcal{M}\) in 14 maps from \(\mathcal{B}(L^{2}(\mathbb{D}))'\) to \(L_{\mathbb{C}}^{\infty}(\Omega)'\), and both of these spaces are inconvenient. However, \(\mathcal{B}(L^{2}(\mathbb{D}))\) has a nice predual: we have \(\mathcal{S}_{1}(L^{2}(\mathbb{D})) ' = \mathcal{B}(L^{2}(\mathbb{D}))\) with respect to the dual pairing \(\langle A,B\rangle:=\mathop{\mathrm{tr}}(B^{\ast}A)\) (see, e.g., [19]). It follows that \(\mathcal{S}_{1}(L^{2}(\mathbb{D}))\subset \mathcal{S}_{1}(L^{2}(\mathbb{D}))''=\mathcal{B}(L^{2}(\mathbb{D}))'\) in the sense of the canonical embedding of a Banach space into its bidual. With these preparations, we can formulate the following proposition:

We can now characterize the adjoint of the Fréchet derivative as follows:

4 Uniqueness↩︎

4.1 Unavoidable ambiguities↩︎

A crucial aspect of the mathematical analysis of phase retrieval problems is identifying all potential ambiguities. In this discussion, we will identify unavoidable ambiguities inherent in both linear and non-linear inverse problems, and show that under certain conditions these are all sources of non-uniqueness.

We discuss two further sources of non-uniqueness:

Note that we can avoid this source of non-uniqueness by confining ourselves to continuous functions \(f\).

The non-uniqueness in [rem:nonunique95local95shift] does not occur in the linearized case, but the validity of the linearization becomes questionable for phase jumps in the order of \(\geq 2\pi\). We discuss a third type of non-uniqueness for the nonlinear problem.

4.2 Uniqueness of the nonlinear forward problem↩︎

3 shows that some assumption is needed which excludes symmetry. We will assume that the sample is contained “in one half” of the illuminated area \(\mathbb{D}\) (see 7 ) in the following sense (see Fig. 1 (a)):

Recall that the essential support is defined by \(\mathop{\mathrm{ess\,supp}}f:=\mathbb{D}\setminus\bigcup\{U:U\subset \mathbb{D}\) relatively open, \(f|_U=0\text{ a.e.}\}\). 2 can be checked whenever a bound on the sample’s support is known. However, it has the drawback that, for a given beam, the size of the samples that can be imaged with uniqueness guarantees is reduced by a factor of \(2\).

Note that 3 excludes a nontrivial block decomposition of \(K\) as in 2. Such a block structure is a highly unlikely scenario unless multiple independent beams are employed simultaneously; nevertheless, 3 may still be violated for other reasons. A situation where 3 is highly plausible is the pinhole setting discussed in Section 2 if the coherence length of the incident beam is larger than the diameter of the pinhole.

We are now in a position to state our main theoretical result:

Note that because of the non-uniqueness discussed in [rem:nonunique95local95shift], we cannot conclude that \(f_1-f_2\equiv i(c_1-c_2)\) under the given assumptions. This would be possible, however, if we additionally assume continuity of \(f_1\) and \(f_2\).

As a first step, we establish the following lemma.

Let \(\mathop{\mathrm{ch}}(A)\) denote the convex hull of a set \(A \subset \mathop{\mathrm{\mathbb{R}}}^m\), i.e., the smallest convex set containing \(A\). The essential tool in the proof of the main uniqueness theorem is the following identity for the convex hull of the support of convolutions due to J.-L. Lions ([20]): if \(k_{1}\) and \(k_{2}\) are two distributions with compact support, then \[\begin{align} \label{eq:chull95conv} \mathop{\mathrm{ch}}\mathop{\mathrm{supp}}(k_1*k_2) = \mathop{\mathrm{ch}}\mathop{\mathrm{supp}}(k_1) \oplus \mathop{\mathrm{ch}} \mathop{\mathrm{supp}}(k_2), \end{align}\tag{16}\] where “\(\oplus\)” denotes Minkowski sum and is defined by \(A\oplus B:=\{a+b:a\in A,b\in B\}\) for \(A,B \subset \mathbb{R}^m\).

We also state a corresponding uniqueness result for the linearized problem:

5 sec:Reconstructions↩︎

In this section, our focus is on effective numerical regularization methods to reconstruct jointly the phase and the absorption contrast. Our numerical implementation unfolds in two steps: (i) discretization and noise model, (ii) iterative regularization method.

5.1 Discretization and noise model↩︎

We assume that the illuminated area \(\mathbb{D}\) is contained in the square \([-1,1]^2\) discretize the complex-valued object \(f\in\mathbb{C}^{M_f\times M_f}\) as an equidistant grid dividing it into \(M_f\times M_f\) pixels. Similarly, the measurement domain \(\mathbb{M}\) is divided into \(M_I\times M_I\) equidistant pixels. We implement the Fresnel propagator \(\mathcal{D}\) in the convolution form 4 by FFT using zero-padding of \(\mathbb{D}\) to reduce periodization artifacts and restrictions in Fourier space to adapt the distance and number of detector pixels.

The discrete forward operator is then given by \[\label{GI-discrete} F:\mathbb{C}^{M_f\times M_f}\to\mathbb{R}^{M_I^2\times M_I^2},~~~F(f):=\big| \mathcal{D}\mathop{\mathrm{diag}}(e^{f})\mathop{\mathrm{\mathbf{Cov}}}[u]\mathop{\mathrm{diag}}(e^{\overline{f}})\mathcal{D}^{\ast}\big|^2,\tag{28}\] where \(\mathop{\mathrm{diag}}(e^f)\in \mathbb{C}^{M_f^2\times M_f^2}\) is the diagonal matrix with diagonal \(e^f\).

To generate synthetic data, we draw \(N\) samples \(u_{n}\), \(n=1,\dots,N\) of the Gaussian random field \(u\) and compute the corresponding intensities \[I_{f,n}=|\mathcal{D}\mathop{\mathrm{diag}}(e^{f}) u_{n}|^{2},~~f\in\mathbb{C}^{M_f\times M_f},~~n=1,\dots,N.\] Observed intensities are affected by photon shot noise described by a normalized Poisson process modeled as \[\label{eq:primary95data} I_{f,n}^{\text{obs}}\sim\frac{1}{t}\mathop{\mathrm{Pois}}(tI_{f,n}),\qquad n=1,\dots,N,\tag{29}\] where the parameter \(t > 0\) can be interpreted as the observation time (proportional to the expected number of photon counts) per image. The scaling factor \(\frac{1}{t}\) ensures that \(\mathop{\mathrm{\mathbb{E}}}[I_{f,n}^{\textrm{obs}}|u_n]=I_{f,n}\) (see [21]).

The \(N\) intensity images of size \(M_I\times M_I\) described in 29 are samples of a discrete Cox process and model our primary data. From these data we could in principle compute an estimator of the noise-free intensity correlations \(\operatorname{cov}[I_f]\), which served as input data of the inverse problem in the previous sections by \[\widehat{\operatorname{cov}}_{N,t}[I_{f}]:=\frac{1}{N}\sum_{n=1}^{N}(I_{f,n}^{\text{obs}}-\overline{I}_{f,N}^{\text{obs}})(I_{f,n}^{\text{obs}}-\overline{I}_{f,N}^{\text{obs}})^{\top},~~\text{with}~~\overline{I}_{f,N}^{\text{obs}}:=\frac{1}{N}\sum_{n=1}^{N}I_{f,n}^{\text{obs}}.\] Correction terms on the diagonal due to the Poisson process are discussed in 8, but they are negligible in our setting since they scale like \(\frac{1}{t}\) for large count rates.

Since the size of \(\widehat{\operatorname{cov}}_{N,t}[I_{f}]\) is proportional to \(M_I^4\), the computation of \(\widehat{\operatorname{cov}}_{N,t}[I_{f}]\) severely limits the size of \(M_I\) on available computational resources and impairs the practicality of the method. Therefore, it will be crucial to avoid an explicit computation of \(\widehat{\operatorname{cov}}_{N,t}[I_{f}]\) and only work with the primary data \(I_{f,n}^{\text{obs}}\).

5.2 Iterative regularization method↩︎

We wish to solve the discrete nonlinear inverse problem \[\label{eq:discrete32Inv-Probl} F(f)\approx \widehat{\operatorname{cov}}_{N,t}[I_{f}].\tag{30}\] To improve stability of reconstructions it is typically advantageous to incorporate any available prior information on the solution into the reconstruction process. In the following we will only consider the constraint that phase contrast for x-rays in non-positive and absorption contrast is non-negative [22] , i.e., \(f\) belongs to the set \[\label{non-negativity32cons} \mathbb{K}_{+}:=\{f\in\mathbb{C}^{M_f\times M_f}: -\mathop{\mathrm{Re}}f\geq 0~\text{and}~\mathop{\mathrm{Im}}f\geq 0\;\text{pointwise} \}.\tag{31}\] If additional support constraints are available, they are also easy to implement in this context. To solve 30 both constraints are incorporated to the following Tikhonov regularization problem: \[\label{eq:Tikh-minimization32prob} \Bar{f}\in\mathop{\mathrm{argmin}}\left[\mathcal{S}(F(f))+\alpha \left({\mathpalette\irchi\relax}_{\mathbb{K}_{+}}(f)+\left\| f-f_{0} \right\|_{2}^{2}\right)\right],\qquad \mathcal{S}(g):=\frac{1}{2}\left\| g-\widehat{\operatorname{cov}}_{N,t}[I_{f}] \right\|_{\mathop{\mathrm{\mathcal{HS}}}}^{2}\tag{32}\] Here the norm \(\left\| \cdot \right\|_{\mathop{\mathrm{\mathcal{HS}}}}\) in the data fidelity term denotes the Hilbert-Schmidt (or Frobenius) norm, \(\alpha>0\) signifies the regularization parameter, \(f_{0}\) in the penalty term is an initial guess, and \({\mathpalette\irchi\relax}_{\mathbb{K}_{+}}\) denotes the characteristic function of the set \(\mathbb{K}_{+}\), i.e., \(\mathbb{K}_+(f):=0\) if \(f\in \mathbb{K}_+\) and \(\mathbb{K}_+(f)=\infty\) else.

Avoiding covariance fitting↩︎

If we wish to minimize the Tikhonov functional 32 by gradient methods, we need to compute, in particular, the gradient of the data fidelity term, which is given by \[\mathop{\mathrm{\mathbf{grad}}}(\mathcal{S}\circ F)(f) = F'[f]^{\ast} F(f) - F'[f]^{\ast}\widehat{\operatorname{cov}}_{N,t}[I_{f}].\] If we wish to use Gauß-Newton type methods, i.e., approximate \(\mathcal{S}(F(f))\) by \(\mathcal{S}(F(f_n)+F'[f_n](f-f_n))\) in an iterative process, we additionally need to evaluate terms of the form \(F'[f_n]^{\ast}F'[f_n]h\).

However, evaluating any of these terms in a straightforward manner, in particular pre-processing intensities \(I_{f,n}^{\text{obs}}\) to compute \(\widehat{\operatorname{cov}}_{N,t}[I_{f}]\) prohibitively increases data dimensionality. In X-ray holography imaging, data sets can comprise millions pixels, which would result in the order of \(10^{12}\) independent two-point correlations! This is way too large to be stored in fast memory on most machines.

The crucial point to make this computable for large problem instances is to assume that \(\mathop{\mathrm{\mathbf{Cov}}}[u]\) has small rank \(r\ll m\) and to avoid any \(m\times m\) matrices and in particular any elements of the image space of the forward map. We need at most \(m\times r^2\) or \(r^2\times m\) matrices.

This can be achieved by the following lemma providing a factorization of the pointwise squared modulus of a product of low-rank matrices into a product of low-rank matrices:

To see how Lemma 5 allows to avoid \(m\times m\) matrices, let us introduce the function \[\Theta: \mathop{\mathrm{\mathbb{C}}}^{m\times (r\times r)} \to \mathop{\mathrm{\mathbb{C}}}^{m\times m}, \qquad \Theta(E):=EE^{\ast},\] for the matrix product such that ?? with \(B=C\) becomes \[|BB^{\ast}|^2 = \Theta(\tau(B)).\]

The key point about the formulas in 6 is that thanks to the bracketing in ?? they completely avoid the formation of \(m\times m\) matrices.

If the covariance matrix \(\mathop{\mathrm{\mathbf{Cov}}}[u]\) is of rank \(r\ll m\), then it has a factorization \(\mathop{\mathrm{\mathbf{Cov}}}[u] = VV^{\ast}\) where \(V\) has \(r\) columns and the forward operator corresponding to the correlation data reads as \[F(f) := |B(f)B(f)^{\ast}|^{2} = \Theta(\tau(B(f)))\qquad with \qquad B(f) :=\mathcal{D} \mathop{\mathrm{diag}}(e^f)V,\] and by the chain rule we can compute \(F'[f]^{\ast}F'[f]\) using matrices no larger than \(m\times r^2\).

The forward operator corresponding to mean intensity data is given by \[F_{\mathrm{mean}}(f) := \mathop{\mathrm{Diag}}|B(f)B(f)^{\ast}|^{2} = \Xi(B(f)) \qquad with\qquad \Xi(B):=\left(\sum\nolimits_{j=1}^{r}|B_{kj}|^2\right)_k.\]

:show95measurementsVisualization of the partial coherent incident beam and the fractionated data.

5.3 Numerical implementation↩︎

Software↩︎

All reconstructions were produced by the inverse problems python library “RegPy" [23]. It provides tools to implement custom forward models as well as a variety of regularization methods and stopping rules. Further details can be found on Git https://github.com/regpy/regpy.

Minimization of the Tikhonov functional↩︎

The minimization problem in 32 is solved by the generalized Fast iterative shrinkage-thresholding algorithm (G-FISTA) [24], [25]. To improve convergence and stability, we employ a homotopy (continuation) strategy in the regularization parameter \(\alpha\), where the problem is solved sequentially for decreasing values of \(\alpha\), each initialized with the previous iterate. We choose regularization parameters \(\alpha_{k}=\alpha_{0} c^{k}\) for some \(c\in (0,1)\) to decrease gradually with initial regularization parameter \(\alpha_{0}\) for both intensity correlations and mean intensity.

Reconstruction results↩︎

We numerically analyze the performance of G-FISTA for holographic X-ray phase contrast imaging from both intensity correlations and mean intensity data. To demonstrate the feasibility of this approach, we implement the forward operators corresponding to intensity correlations and mean intensity data for a two-dimensional cell test pattern with \(256 \times 256\) pixels introduced in [26] as contrast \(f\). To produce synthetic data reconstructions, we use a partially coherent beam and \(N=3000\) frames. The regularization parameters are chosen according to \(\alpha\in \{10^{-9}\cdot(\frac{1}{3})^k: k\in\mathbb{N}\}\). For mean intensity, the inversion with G-FISTA is started with an initial guess \(f_{0}=0\), while for intensity correlations, a warm start from mean intensity reconstruction is chosen. For the holographic regime, a Fresnel number \(\mathfrak{f}=\frac{10}{2\pi}\) is chosen for both forward operators.

To construct \(V \in \mathbb{C}^{m\times r}\), we generate \(r=4\) independent samples of smooth random fields as follows: we draw i.i.d. Gaussian Fourier coefficients, multiply them by a rapidly decaying spectral filter of the form \(e^{-\frac{|\xi|^{2}}{\sigma^{2}}}\) with \(\sigma=0.5\), and apply an inverse Fourier transform to obtain spatially smooth realizations. A spatial cutoff is then applied to localize the fields. The resulting samples are orthonormalized via a singular value decomposition to obtain the columns of \(V\).

Throughout the regularization process, we apply the \(L^{2}\)- norm for both data fidelity and penalty term together with a non-negativity constraint for both phase and absorption contrasts. Additionally, the data are polluted with shot noise described by a Cox process with observation time \(t=10^{9}\) photon counts per image.

A comparison of the holographic X-ray phase contrast imaging from intensity correlations and mean intensity data is shown in 2. The joint reconstruction results for mean intensity are achieved using a 5-stage process with 100 G-FISTA iterations per stage, while for intensity correlations, we use a 2-stage process with 50 G-FISTA iterations per stage. The results show that intensity correlations enable simultaneous reconstruction of both phase and absorption contrasts, whereas mean intensities alone do not yield faithful reconstructions. It appears that in the case of mean intensity data, some ghost-like patterns produce artifacts in the reconstructions. The comparison demonstrates that a significant amount of additional information can be achieved by the correlation data compared to the mean intensity data. Numerical simulations, complemented by preliminary theoretical results, suggest that the reconstructions become more accurate as the number of frames, \(N\), increases.

Additionally, we show in [fig:show95measurements] the 6 eigenvectors (or principal components) with the highest eigenvalues. These components correspond to the dominant singular values and capture the principal variations of the input random fields, which is consistent with the stochastic nature of SASE pulses. Each component maximizes the variation of the random fields in a subspace, which is orthonormal to the previous components. According to these results, there is a good agreement between theory and practice, as can be seen by all-at-once phase and absorption reconstruction.

6 Conclusion and outlook↩︎

We have theoretically studied holographic X-ray imaging by exploring and exploiting the information content of the intensity correlations. Specifically, we proved unique identifiability for the nonlinear problem up to unavoidable sources of non-uniqueness.

For large-scale problems, the increased data dimensionality giving rise to a major difficulty when we compute the correlations from the pre-processing intensities. To remedy this we have avoided covariance fitting in the process of regularization by assuming that the covariance matrix is known and has a low rank. By 5 and 6, we showed that the quantitative holographic imaging can then effectively be interpreted as the application of forward-backward propagation, exploiting implicitly the full information content of correlation data. Moreover, the computational cost of the algorithm scales linearly in the pixel number \(m\) and quadratically in the rank \(r\) of \(\mathop{\mathrm{\mathbf{Cov}}}[u]\).

Let us close this manuscript by introducing possible directions for future research: One such direction is to replace the “symmetry-breaking condition” ?? by different conditions that may be better suited for specific applications. Furthermore, one could analyze the analogous, mathematically more challenging problem for the full Helmholtz equation instead of the Fresnel approximation, which might be necessary in different applications with smaller wave numbers of the incident beam. Another natural aim would be to quantify the stability of 13 to find \(f\) from observed correlations depending on the cross-covariance, the Fresnel number, and on the a-priori information. A further direction is to extend our methodology to meet our forward problem when \(\mathop{\mathrm{\mathbf{Cov}}}[u]\) is unknown. This is more applicable in real world problems. The idea is to first estimate the covariance matrix of the empty beam by modal decomposition of SASE pulses and then exploiting the information in an advanced method to get faithful reconstruction.

7 Circularly symmetric Gaussian random vectors↩︎

Recall that a random variable \(Z\) with values in \(\mathbb{C}^m\) is called complex Gaussian if \((\mathop{\mathrm{Re}}Z,\mathop{\mathrm{Im}}Z)\) is multivariate Gaussian in \(\mathbb{R}^{2m}\) and circularly symmetric if \(e^{i\varphi}Z\) has the same distribution as \(Z\) for all \(\varphi\in \mathop{\mathrm{\mathbb{R}}}\). Obviously, if the expectation of a circularly symmetric random vector exists, it must vanish. Let \(Z\) be a circularly symmetric Gaussian random vector. Then \(\mathop{\mathrm{\mathbb{E}}}[ZZ^{\top}]=\mathop{\mathrm{\mathbb{E}}}[e^{i\varphi}Ze^{i\varphi}Z^{\top}] = e^{2i\varphi} \mathop{\mathrm{\mathbb{E}}}[ZZ^{\top}]\), so \[\label{eq:aux95circ95symm} \mathop{\mathrm{\mathbb{E}}}[Z Z^{\top}]= 0\quad if Z is Gaussian and complex symmetric.\tag{33}\] In particular, if \(Z\) is scalar, then taking the imaginary part of \(\mathop{\mathrm{\mathbb{E}}}[Z^2]=0\) shows that \(\mathop{\mathrm{Re}}Z\) and \(\mathop{\mathrm{Im}}Z\) are uncorrelated, which by Gaussianity implies that they are independent.

If \(X\) and \(Y\) are circularly symmetric Gaussian random variables, then Isserlis’ theorem [27] yields the identity \[\mathop{\mathrm{\mathbb{E}}}\left[X\overline{X}Y\overline{Y}\right] = \mathop{\mathrm{\mathbb{E}}}[X\overline{X}]\mathop{\mathrm{\mathbb{E}}}[Y\overline{Y}] + \mathop{\mathrm{\mathbb{E}}}[X\overline{Y}] \mathop{\mathrm{\mathbb{E}}}[\overline{X}Y] + \mathop{\mathrm{\mathbb{E}}}[XY]\mathop{\mathrm{\mathbb{E}}}[\overline{X}\overline{Y}].\] The last term vanishes due to 33 , and the remaining terms can be rearranged to \[\label{eq:Cov95circ95symm} \mathop{\mathrm{Cov}}(|X|^2,|Y|^2)=|\mathop{\mathrm{Cov}}(X,Y)|^2.\tag{34}\] In particular, for \(X=Y\) using \(\mathop{\mathrm{\mathbb{E}}}[X]=0\) we obtain \[\label{eq:E95Var} \boldsymbol{Var}(|X|^2) = \mathop{\mathrm{\mathbb{E}}}[|X|^2]^2.\tag{35}\]

8 Cox processes↩︎

Given a random non-negative function (or measure) \(I\) on a domain \(\mathbb{M}\), a Cox process \(P\) with mean intensity \(I\) is a point process \(P=\sum_{i=1}^N \delta_{x_i}\) which, conditioned on \(I\) is a Poisson process with intensity \(I\). Here both the points \(x_i\in\mathbb{M}\) and the total number of points \(N\) are random. For a more thorough characterization of Cox processes we refer to [28]. For a continuous function \(f:\mathbb{M}\to \mathbb{R}\), let \[\langle P,f\rangle:=\sum\nolimits_{i=1}^N f(x_i).\] If \(c_I(x,y):=\mathop{\mathrm{Cov}}(I(x),I(y))\) and \(f,g:\mathbb{M}\to \mathbb{R}\) are two continuous functions we have \[\mathop{\mathrm{Cov}}\left(\langle P,f\rangle,\langle P,g\rangle\right) = \int_{\mathbb{M}}\int_\mathbb{M}c_I(x,y)f(y)g(x)\,\mathrm{d}y\mathrm{d}x + \int_{\mathbb{M}} (\mathop{\mathrm{\mathbb{E}}}[I])(x) f(x) g(x)\,\mathrm{d}x\] (see [28] for the case of indicator functions \(f,g\) which implies the above formula by density). In this sense we have \[\mathop{\mathrm{Cov}}[P] = \mathop{\mathrm{Cov}}[I] + M_{\mathop{\mathrm{\mathbb{E}}}[I]}.\] Note that \(\mathop{\mathrm{Cov}}[I]\) is quadratic in \(I\) whereas \(\mathop{\mathrm{\mathbb{E}}}[I]\) is linear in \(I\). For the count rates considered in this paper we found the second term to be negligible compared to the first.

Concerning uniqueness the correction term \(M_{\mathop{\mathrm{\mathbb{E}}}[I]}\) does not play a role since its Schwartz kernel is supported only on the diagonal \(\{(x,x):x\in\mathbb{M}\}\), which is a nullset in \(\mathbb{M}\times \mathbb{M}\).

Acknowledgments↩︎

We would like to thank Tim Salditt for many helpful discussions. Financial support by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through Grant 432680300, SFB 1456— Mathematics of Experiments, project C03, is gratefully acknowledged.

References↩︎

---

1. Institute for Numerical and Applied Mathematics, Lotzestraße 16-18, 37083 University of Göttingen, Germany and Max-Planck Institute for Solar System Research, 37077 Göttingen, Germany ()↩︎

2. Corresponding author. Institute for Numerical and Applied Mathematics, Lotzestraße 16-18, 37083 University of Göttingen, Germany ()↩︎

3. Max-Planck Institute for Solar System Research, 37077 Göttingen, Germany ()↩︎