Exact Crystalline Structure Recovery in X-ray Crystallography from Coded Diffraction Patterns

TL;DR

This paper introduces a coded aperture system in X-ray crystallography that guarantees exact structure recovery with fewer measurements, reducing exposure time and preserving crystal integrity.

Contribution

It presents a novel coded diffraction pattern approach with proven exact reconstruction guarantees and an efficient algorithm leveraging sparsity in Fourier domain.

Findings

Exact reconstruction guarantees up to a global phase shift.

Requires 50% fewer measurements than existing methods.

Reduces crystal exposure time, minimizing damage.

Abstract

X-ray crystallography (XC) is an experimental technique used to determine three-dimensional crystalline structures. The acquired data in XC, called diffraction patterns, is the Fourier magnitudes of the unknown crystalline structure. To estimate the crystalline structure from its diffraction patterns, we propose to modify the traditional system by including an optical element called coded aperture which modulates the diffracted field to acquire coded diffraction patterns (CDP). For the proposed coded system, in contrast with the traditional, we derive exact reconstruction guarantees for the crystalline structure from CDP (up to a global shift phase). Additionally, exploiting the fact that the crystalline structure can be sparsely represented in the Fourier domain, we develop an algorithm to estimate the crystal structure from CDP. We show that this method requires 50% fewer measurements to estimate the crystal structure in comparison with its competitive alternatives. Specifically, the proposed method is able to reduce the exposition time of the crystal, implying that under the proposed setup, its structural integrity is less affected in comparison with the traditional. We discuss further implementation of imaging devices that exploits this theoretical coded system.

Equations

$$g_{p,r}^i=|⟨{\mathbf{a}}_{p,i},{\tilde{\mathbf{x}}}_r⟩{|}^2,i=1\cdots ,n,$$

$${\mathbf{a}}_{p,i}=\overline{\hat{\mathbf{D}}}\mathbf{F}\overline{\mathbf{T}}(z_p){\mathbf{f}}_{u_i},$$

$$g_{p,r}^i=|⟨{\mathbf{b}}_{p,i}^r,\mathbf{x}⟩{|}^2,i=1\cdots ,n,$$

$$\underset{\mathbf{x}\in {\mathbb{C}}^n}{\min }h(\mathbf{x})$$

$$\text{ s.t }\Vert \mathbf{Fx}{\Vert }_0$$

$$g_{p,r}^i=({\mathbf{b}}_{p,i}^r{)}^H\mathbf{x}{\mathbf{x}}^H{\mathbf{b}}_{p,i}^r.$$

$$\mathcal{B}(\mathbf{W})={\left[({\mathbf{b}}_{1,1}^1{)}^H\mathbf{W}{\mathbf{b}}_{1,1}^1,\cdots ,({\mathbf{b}}_{P,n}^R{)}^H\mathbf{W}{\mathbf{b}}_{P,n}^R\right]}^T,$$

$$\mathbf{g}=\mathcal{B}(\mathbf{x}{\mathbf{x}}^H).$$

$${\mathcal{T}}_{\mathbf{x}}=\left\{\mathbf{x}{\mathbf{w}}^H+\mathbf{w}{\mathbf{x}}^H|\mathbf{w}\in {\mathbb{C}}^n\right\},$$

$$(1-\delta )\Vert \mathbf{W}{\Vert }_1\le \frac{1}{\beta }\Vert \mathcal{B}(\mathbf{W}){\Vert }_1\le (1+\delta )\Vert \mathbf{W}{\Vert }_1,$$

$$\mathcal{P}\left(\frac{1}{m}\Vert \mathbf{B}{\Vert }_{\infty }^2\le 1+\delta \right)\le 1-2e^{-c_{0}m{ϵ}^2},$$

$$\underset{\mathbf{x}\in {\mathbb{C}}^n}{\min }f(\mathbf{x})$$

$$\text{ s.t }\Vert \mathbf{Fx}{\Vert }_0$$

$$\partial f\left({\mathbf{z}}^{(t)},{\mu }^{(t)}\right)=\frac{2}{m}$$

$$\sqrt{g_{p,r}^i}\frac{({\mathbf{b}}_{p,i}^r{)}^H{\mathbf{z}}^{(t)}}{{\phi }_{{\mu }^{(t)}}(|⟨{\mathbf{b}}_{p,i}^r,{\mathbf{z}}^{(t)}⟩|)}){\mathbf{b}}_{p,i}^r.$$

$$\text{relative error :=}\frac{dist(\mathbf{z},\mathbf{x})}{\Vert \mathbf{x}{\Vert }_2},$$

$$\text{dist}(\mathbf{z},\mathbf{x})=\underset{\theta \in [0,2\pi )}{\min }\Vert \mathbf{x}{e}^{-j\theta }-\mathbf{z}{\Vert }_2,$$

$$\Vert \mathcal{B}(\mathbf{W}){\Vert }_1$$

$$\le \sum _{i=1}^n\sum _{r=1}^R\sum _{p=1}^P|{\lambda }_1||⟨{\mathbf{b}}_{p,i}^r,\mathbf{u}⟩{|}^2+|{\lambda }_2||⟨{\mathbf{b}}_{p,i}^r,\mathbf{v}⟩{|}^2$$

$$=|{\lambda }_1|\Vert \mathbf{Bu}{\Vert }_2^2+|{\lambda }_2|\Vert \mathbf{Bv}{\Vert }_2^2\le \Vert \mathbf{W}{\Vert }_1\Vert \mathbf{B}{\Vert }_{\infty }^2,$$

$${\mathbf{B}}^H\mathbf{B}$$

$$=\overline{\hat{\mathbf{D}}}\hat{\mathbf{D}},$$

$$\mathcal{P}\left({\Vert \mathbf{B}\Vert }_{\infty }\ge \sqrt{m}+C\sqrt{s}+t\right)\le 2e^{-c_0{t}^2},$$

$$\mathcal{P}\left(\frac{1}{m}\Vert \mathbf{B}{\Vert }_{\infty }^2\le 1+\delta \right)\le 1-2e^{-c_{0}m{ϵ}^2},$$

$$\frac{1}{m}\Vert \mathcal{B}(\mathbf{W}){\Vert }_1\le (1+\delta )\Vert \mathbf{W}{\Vert }_1,$$

$$\Vert \mathcal{B}(\mathbf{W}){\Vert }_1$$

$$={\lambda }_1\Vert \mathbf{Bu}{\Vert }_2^2+{\lambda }_2\Vert \mathbf{Bv}{\Vert }_2^2$$

$$=\left({\lambda }_1+{\lambda }_2\right)=\Vert \mathbf{W}{\Vert }_1,$$

$$\frac{1}{m}\Vert \mathcal{B}(\mathbf{W}){\Vert }_1\ge \frac{1}{m}\left(1-\delta \right)\Vert \mathbf{W}{\Vert }_1,$$

Full text

††preprint: APS/123-QED

I Introduction

X-ray Crystallography (XC) is known as the leading technique for molecular structure characterization in material analysis [1]. Specifically, XC plays an essential role in fields as biology [2], drug design [3], and material sciences [4], among others. The traditional acquisition system in XC is illustrated in Fig.1, where the acquired data forms “the Ewald sphere” according to the perturbation theory which is mathematically modeled as the Fourier magnitude of the unknown crystalline structure called diffraction patterns [5]. Specifically, in Fig.1 the resulting diffraction patterns, when an X-ray source irradiates a crystal, are recorded while both the sensor and the crystal gradually rotate [6].

The crystal structure can be uniquely estimated from the phase of its diffraction patterns [5]. Although, the phase of the diffraction patterns cannot be directly measured, it can be recovered from their phaseless intensity measurements [7], problem that is known as phase retrieval (PR) [7, 8]. Mathematically, the PR problem in XC consists on retrieving a discrete version of the crystalline structure $\mathbf{x}\in\mathbb{C}^{n}$ from the phaseless measurements $y_{k}=\lvert\mathbf{f}_{k}^{H}\mathbf{x}\rvert^{2}$ where $n$ is the discrete size of the crystal, $k=1,\cdots,m$ with $m$ the total number of measurements, $\mathbf{f}_{k}\in\mathbb{C}^{n}$ are the vectors that models the three-dimensional discrete Fourier transform [9], and $(\cdot)^{H}$ represents the conjugate transpose operation. The PR problem has been challenging due to its infinite solutions [7, 10]. Specifically, there are so-called trivial ambiguities that are always present [7]. The following three transformations (or any combination of them) conserve Fourier magnitude: 1) global shift phase: $\mathbf{x}[a]=\mathbf{x}[a]e^{j\phi_{0}}$; 2) conjugate inversion: $\mathbf{x}[a]=\overline{\mathbf{x}[-a]}$; 3) spatial shift: $\mathbf{x}[a]=\mathbf{x}[a+a_{0}]$. Moreover, the number of measurements $m$ must satisfy $m\geq 4n-1$ to accurately retrieve the signal $\mathbf{x}$ (up to trivial ambiguities) [11, 10]. i.e. a 4-fold time exposition X-ray radiation is required, which leads to degradation of the crystalline structure [12, 13, 14, 15, 16, 17]. Additionally, if the crystal deteriorates under this irradiation, material characterization becomes more challenging [14].

Despite such challenges, several methods have been develop to retrieve the crystalline structure. To name a few, the error reduction method [18], dual space programs SnB and ShelxD [19], iterated projections [20], and the charge flipping algorithm [21]. All these methods alternate between real and reciprocal space by the Fourier transform by imposing constraints on the real-space charge density. All these algorithms require to over-expose the crystal to the X-ray source in order to acquired the large enough amount of diffraction patterns needed to estimate the crystalline structure i.e. they work under the $m\geq 4n-1$ regime. Further, these approaches tend to return inaccurate estimates of the crystalline structure when the measurements $y_{k}$ are corrupted by noise [7]. Also, the convergence to the true crystalline structure of the above methods is not guaranteed.

This work proposes to modify the traditional system by including an optical element called coded aperture which modulates the diffracted field to acquire coded diffraction patterns (CDP). The main advantage of the proposed system is that only the global phase shift ambiguity appears, that is, the conjugate inversion and spatial shift ambiguities will not occur [22, 23]. This implies that retrieving the crystalline structure from CDP is much effortless compared with the traditional XC. For the proposed coded system, we derive exact reconstruction guarantees for the crystalline structure from CDP (up to a global shift phase). Additionally, given the fact that the crystalline structure is a periodic element, it can be sparsely represented in the Fourier domain, i.e. the number of non-zero coefficients of the Fourier transformed crystal is much smaller than $n$. Thus, exploiting the sparsity prior of the crystal we develop an algorithm to estimate the crystalline structure from CDP that requires 50% fewer measurements in comparison with its competitive alternatives. Specifically, the proposed method is able to reduce the exposition time of the crystal, implying that under the proposed setup, its structural integrity is less affected in comparison with the traditional. Numerical simulations are conducted to evaluate the performance of the proposed method using synthetic data.

II Problem Formulation

In contrast to the traditional acquisition system illustrated in Fig. 1, a coded diffraction system, as shown in Fig. 2, is proposed. Notice that Fig. 2 includes an optical element known as coded aperture, modeled as $\mathbf{D}$, at distance $z_{p}$ from the sensor. This optical element modules the diffracted field before being measured by the sensor. Besides, changing the distance $z_{p}$, this acquisition system allows acquiring multiple snapshots, where $p=1,\cdots,P$ indexes the sensing distances. The signal $\tilde{\mathbf{x}}\in\mathbb{C}^{n}$ in Fig. 2 corresponds with the Fourier transform of $\mathbf{x}$, that is, $\tilde{\mathbf{x}}=\mathbf{F}\mathbf{x}$, where $\mathbf{F}$ represents the 3D discrete Fourier transform matrix.

Observe that in Fig. 2, we assume that the coded aperture $\mathbf{D}$ is located at the far field from the crystal. Moreover, the distance $z_{p}$ between the coded aperture and the sensor is chosen so that the coded diffracted field, captured at the detector, belongs to the near field. Additionally, we assume that the coded aperture and the sensor jointly perform a raster scanning across “the Ewald sphere”, as illustrated in Fig. 3, while the crystal remains fixed. Thus, in contrast with the traditional architecture, crystal rotation is avoided. Mathematically, a given highlighted region in Fig. 3 is modeled as $\tilde{\mathbf{x}}_{r}=\mathbf{S}_{r}\tilde{\mathbf{x}}$ where $\mathbf{S}_{r}$ is a selection diagonal matrix where $r=1,\cdots,R$ with $R$ as the number of regions.

In order to easily model the acquired data in Fig. 2, $\hat{\mathbf{D}}\in\mathbb{C}^{n\times n}$ is defined as a diagonal matrix whose entries are the elements of $\mathbf{D}$. Thus, from the traditional diffraction theory [23, 24], the measurements being captured at the detector for the $r$-th region and at a distance $z_{p}$, are given by

[TABLE]

where $\mathbf{a}_{p,i}$ are the sampling vectors defined as

[TABLE]

where $u_{i}=(i-1)\mod n+1$, $\mathbf{T}(z_{p})$ is the spatial frequency transfer function which depends on the distance $z_{p}$ [23, 24], and $\overline{(\cdot)}$ represents the conjugate operation. Now, considering the definition of $\tilde{\mathbf{x}}_{r}$ we have that (1) can be equivalently expressed as

[TABLE]

where $\mathbf{b}^{r}_{p,i}=\mathbf{F}^{H}\mathbf{S}_{r}\mathbf{a}_{p,i}$. Observe that under the proposed acquisition system, the maximum number of measurements is given by $m\leq nRP$. Additionally, we remark that the acquired measurements ${g^{i}_{p,r}}$ in Fig. 2 are different than ${y_{k}}$ in Fig. 1 due to the inclusion of the coded aperture.

This paper deals with the problem of recovering $\mathbf{x}\in\mathbb{C}^{n}$ from the phaseless measurements $\{g^{i}_{p,r}\}$. Additionally, we assume that the entries of a coded aperture $\mathbf{D}$ are $i.i.d$ copies of a discrete random variable $d$ obeying $|d|\leq 1$. This assumption over the coded aperture is important because under a random sensing process uniqueness is guaranteed (up to a global phase shift) as will be shown in Section III. This implies that retrieving the crystalline structure from CDP is much easier compared with the traditional XC since the conjugate inversion and spatial shift ambiguities will not occur [22, 23]. Moreover, exploiting the fact that the crystalline structure is periodic, it can be sparsely represented in the Fourier domain, i.e. $\lVert\mathbf{Fx}\rVert_{0}=s\ll n$ where $\lVert\cdot\rVert_{0}$ is the $\ell_{0}$ pseudo-norm, we formulate the following optimization problem

[TABLE]

To solve (4) we propose a gradient descend method that will be described in Section IV.

III Exact Recovery Guarantees

Observe that (3) can be equivalently expressed as

[TABLE]

Now, consider the linear operator $\mathcal{B}:\mathcal{S}^{n\times n}\rightarrow\mathbb{R}^{m}$($\mathcal{S}^{n\times n}$ is the space of self-adjoint matrices) defined as

[TABLE]

and stacking the measurements $\{g^{i}_{p,r}\}$ as $\mathbf{g}:=[g^{1}_{1,1},\cdots,g^{n}_{P,R}]^{T}$, then we have that

[TABLE]

Thus, to prove that the signal $\mathbf{x}$ can be exactly recovered from the measurements $g^{i}_{p,r}$ in (5), from (7) we have that $\mathcal{B}(\cdot)$ must be injective [25, 26]. More precisely, this work follows the strategy in [27] that considers

[TABLE]

as the tangent space of the manifold of all rank-1 Hermitian matrices at the point $\mathbf{x}\mathbf{x}^{H}$. Thus, if the operator $\mathcal{B}$ satisfies the following condition, which is proved in Theorem 1, one can guarantee recovery with high probability[25].

Condition 1**.**

For any $\delta\in(0,1)$ and some constant $\beta>0$ the linear operator $\mathcal{B}$ satisfies

[TABLE]

for all matrices $\mathbf{W}\in\mathcal{T}_{\mathbf{x}}$, where $\|\mathbf{W}\|_{1}=\sum_{i}\sigma_{i}(\mathbf{W})$ with $\sigma_{i}(\mathbf{W})$ as the $i$-th singular value of $\mathbf{W}$.

Theorem 1**.**

Fix any $\delta\in(0,1)$ and the coded aperture $\hat{\mathbf{D}}\in\mathbb{C}^{n\times n}$, with $i.i.d$ copies of a random variable $d$ such that $\lvert d\rvert\leq 1$. Then, we have that

[TABLE]

for some constant $c_{0}>0$ provided that $m\geq Cs$ where $s$ is the sparsity of the crystalline structure with $C>0$ and $\mathbf{B}$ is given by $\mathbf{B}=\left[\mathbf{b}^{r}_{1,1},\cdots,\mathbf{b}^{R}_{P,n}\right]^{H}$. Also, Condition 1 is satisfied with the same probability taking $\beta=m$, where $\lVert\mathbf{B}\rVert_{\infty}$ denotes the spectral norm of $\mathbf{B}$.

The proof is deferred to Appendix A. Theorem 1 establishes two aspects: first, the only ambiguity that appears in the proposed system is the global phase shift. Second, the condition $m\geq 4n-1$ imposed by the traditional system in XC can be defeated by the proposed architecture since now the number of measurements $m$ depends on the sparsity ($s\ll n$) of the crystal, i.e. $m\geq Cs$ for some constant $C>0$. The value of $C>0$ will be numerically estimated in Section V.

IV Crystalline Structure Reconstruction Methodology

In order to solve (4), this work adapted the Sparse Phase Retrieval via Smoothing Function (SPRSF) method introduced in [28]. This algorithm claims to have better performance in terms of the number of measurements compared with recent approaches in the state-of-the-art. SPRSF uses a special mapping called smoothing function, which is useful to eliminate the non-smoothness of $h(\cdot)$. Specifically, the optimization problem solved by SPRSF is formulated as

[TABLE]

where $\varphi_{\mu}(w)$ is defined as $\varphi_{\mu}(w)=\sqrt{w^{2}+\mu^{2}}$. The proposed reconstruction algorithm consists on two stages as follows:

• •

*Initialization step: * this procedure consists on estimating the crystalline structure $\mathbf{x}$ as the leading eigenvector of a carefully designed matrix.

• •

*Refining step: * the outcome of the first step is refined upon a sequence of Wirtinger gradient iterations.

These two stages are summarized in Algorithm 1.

Algorithm 1 requires the sampling vectors and the acquired coded diffraction patterns as modeled in (3) (Line 1). The initialization step is presented in Lines 2-3. Further, a thresholding step is calculated in Line 5, where the operators $\mathcal{H}_{s}(\mathbf{w})$ set all the entries in the vector $\mathbf{w}\in\mathbb{C}^{n}$ to zero, except its $s$ largest absolute values. Additionally if the condition in Line 7 is not satisfied the smoothing parameter $\mu$ is updated in Line 9 to obtain a new point. Remark that each vector $\partial f(\mathbf{z}^{(t)},\mu^{(t)})$ in Algorithm 1 is calculated using the Wirtinger derivative [29] as

[TABLE]

The theoretical guarantees of convergence of this algorithm are demonstrated in [28].

V Simulations and Results

The performance of the proposed algorithm is presented. Three different tests are performed: first, the empirical success of the proposed method is analyzed, among 100 trial runs. The second examines the reconstruction performance of Algorithm 1 for recovering the crystalline structure. The third test studies the stability behavior of the recovery algorithm under additive noise.

The default values of the parameters of Algorithm 1 were determined using a cross-validation strategy. They were fixed as $\tau=0.3$, $\gamma=0.8$, $\gamma_{1}=0.5$, $\mu^{(0)}=60$ and $T=800$. The performance metric used is

[TABLE]

where $dist(\mathbf{z},\mathbf{x})$ is defined as

[TABLE]

with $j=\sqrt{-1}$. The simulated coded aperture was a block-unblock ensemble as used in our previous works [30, 23]. All simulations were implemented in Matlab R2019a on an Intel Core i7 3.41Ghz CPU with 32 GB RAM.

To examine the recovery success rate against the level of sparsity $s$ and the number of measurements $m$ under a noiseless scenario, we randomly simulate 1000 crystalline structures. The success rate is determined over 100 trial runs for each crystalline structure when a relative error of $10^{-5}$ is reached. The results are shown in Fig. 4, where the success rate is plotted in gray scale when white and black represent 100% and 0% probability of success, respectively.

From Fig. 4, it can be concluded that the proposed algorithm is able to estimate the crystalline structure when $m/n\geq 2$ for all the tested sparsity levels. This result implies that the proposed acquisition system, along with its reconstruction approach, is able to estimate the crystal structure using $50\%$ less measurements than the state-of-the-art methods that require $m/n\geq 4$ [11, 10]. Thus, considering the fact that $s\ll n$, then the constant $C$ of Theorem 1 is bounded from above as $C\leq 2$. Observe that in practice we can achieve this reduction by choosing the number of sensing distances as $P=2$.

Finally, the reconstruction accuracy of Algorithm 1 for both noiseless and noisy scenarios is analyzed. The tested crystalline structure is the sodium chloride (NaCl) as illustrated in Fig. 5. Specifically the green spheres in Fig. 5 models the Cl atoms and the yellow ones the sodium atoms. For the noisy case, we study the stability behavior of Algorithm 1 under additive white noise with a Signal-to-Noise-Ratio (SNR) fixed as $SNR=30dB$ where SNR$=20\log(\lVert\mathbf{g}\rVert_{2}/(m\sigma))$ with $\sigma$ the variance of the noise. The number of measurements $m$ used for this experiments satisfies $m/n=2$. The attained reconstructions are shown in Fig. 5 suggesting the effectiveness of the proposed method to estimate the crystalline structure from both noiseless and noisy measurements.

VI Conclusion and Discussion

This work presented an algorithm to recover the 3D structure of a crystal, under a system that records coded diffraction patterns. Our approach takes advantage of the fact that the crystalline structure can be sparsely represented in the Fourier domain to significantly reduce the number of measurements. Simulations show that our approach is able to reconstruct the crystalline structure with up to $50\%$ less amount of measurements compared with the traditional methods in the state-of-the-art. Also, the results suggest that the proposed approach allows reconstructing the crystalline structure even in noisy scenarios.

To implement the proposed acquisition system in Fig. 2 a block-unblock coded aperture is feasible [30, 23]. Specifically, the blocking elements of these coded apertures can be fabricated using tungsten, since this material can stop an x-ray beam, resulting in low fabrication costs [31, 32, 33].

Appendix A Proof of Theorem 1

To prove Theorem 1 we divide it into two parts. First, we prove the right inequality in 9, and then as a second part we prove the left inequality. Thus, let $\mathbf{W}\in\mathcal{T}_{\mathbf{x}}$. As $\mathbf{W}$ has rank at most two, we can choose normalized vectors $\mathbf{u},\mathbf{v}\in\mathbb{C}^{n}$ such that $\mathbf{W}=\lambda_{1}\mathbf{u}\mathbf{u}^{H}+\lambda_{2}\mathbf{v}\mathbf{v}^{H}$. Then considering the definition of the linear map $\mathcal{B}$ in (6) we have that

[TABLE]

in which the first and second inequalities are obtained using the triangular inequality, and matrix $\mathbf{B}$ as defined in Theorem 1. Further, considering definition of matrix $\mathbf{B}$ we have that

[TABLE]

since $\mathbf{F}^{H}\mathbf{F}=\mathbf{F}\mathbf{F}^{H}=\mathbf{I}$, $\mathbf{T}(z_{p})$ is a diagonal orthogonal matrix, and $\sum_{r=1}^{R}\mathbf{S}_{r}\mathbf{S}_{r}=\mathbf{I}$. Thus, using the fact that the admissible random variable $d$ is assumed $|d|\leq 1$ we have from (14) that $\mathbf{B}$ is an isotropic subgaussian matrix [34]. Then, from Theorem 5.39 in [34] we have that

[TABLE]

for constants $c_{0},C>0$ and any $t>0$. Then, taking $m\geq C^{2}\epsilon^{-2}s$ and $t=\sqrt{m}\epsilon$ for any $\epsilon\in(0,1/2)$, we have from (15) that

[TABLE]

for $\delta=2\epsilon$. Thus, we conclude that

[TABLE]

for any $\delta\in(0,1)$.

On the other hand, from (13) we can also conclude that

[TABLE]

in which the third equality comes from observation in (14), using that $\mathbf{W}$ is assumed to be positive semidefinite. Thus, we have that

[TABLE]

for any $\delta\in(0,1)$. Thus, combining (17) and (19) the result holds.