Iterative X-ray Spectroscopic Ptychography

TL;DR

This paper introduces SPA, a novel iterative algorithm for spectroscopic ptychography that improves chemical composition imaging by leveraging spectral redundancy and handling incomplete spectral dictionaries.

Contribution

The paper presents SPA, an innovative algorithm for spectroscopic ptychography that enhances reconstruction quality and robustness, especially with incomplete spectral information.

Findings

SPA outperforms standard methods in reconstruction quality.

SPA effectively handles Poisson noise in measurements.

SPA demonstrates robustness with reduced redundancy data.

Abstract

Spectroscopic ptychography is a powerful technique to determine the chemical composition of a sample with high spatial resolution. In spectro-ptychography, a sample is rastered through a focused x-ray beam with varying photon energy so that a series of phaseless diffraction data are recorded. Each chemical component in the material under investigation has a characteristic absorption and phase contrast as a function of photon energy. Using a dictionary formed by the set of contrast functions of each energy for each chemical component, it is possible to obtain the chemical composition of the material from high resolution multi-spectral images. This paper presents SPA (Spectroscopic Ptychography with ADMM), a novel algorithm to iteratively solve the spectroscopic blind ptychography problem. We design first a nonlinear spectro-ptychography model based on Poisson maximum likelihood, and construct then the proposed method based on fast iterative splitting operators. SPA can be used to retrieve spectral contrast when considering both a known or an incomplete (partially known) dictionary of reference spectra. By coupling the redundancy across different spectral measurements, the proposed algorithm can achieve higher reconstruction quality when compared to standard state-of-the-art two-step methods. We demonstrate how SPA can recover accurate chemical maps from Poisson-noised measurements, and also show its enhanced robustness when reconstructing reduced redundancy ptychography data using large scanning stepsizes.

Tables

Table 1. Table 1: Main variables and operators defined in section 2

Notations	Explanations
${\left\{I_l\right\}}_{l=0}^{L-1}$	Measured intensities
$X\in {ℝ}_{+}^{N,C}$	Elemental thickness maps of the sample
$D\in {ℂ}^{C,L}$	Spectrum dictionary
$Y\in {ℂ}^{N,L}$	Sample spectral contrast maps
$\omega \in {ℂ}^{\overline{m}}$	Probe
${𝒮}_j\in {ℝ}^{\overline{m},m}$	Binary matrix to take image patches
$𝒜:{ℂ}^{\overline{m}}\times {ℂ}^N\to {ℂ}^m$	Forward operator for ptychography
$𝒢(𝒜(w,u),I_l)$	Poisson likelihood estimation
$D_r\in {ℝ}^{C,L}$	Real part of spectrum dictionary
$𝒳$	Non-negative thickness constraint
$𝒲$	Normalized constraint of the probe
TV	Total variation regularization

Table 2. Table 2: SNR in dB from reconstruction results with different scan stepsizes when using the two-step method, SPA and SPA with TV regularization.

stepsize		36	38	40
SNR	two-step	11.5	7.0	3.4
	SPA	15.1	15.5	12.3
	SPA + TV	16.0	16.0	13.8

Equations

$$I_l=\mathrm{Poi}(|\mathcal{A}(\omega ,Y_l){|}^2)\text{leavevmode}\mathrm{for}\text{leavevmode}l=0,1,\cdots ,L-1,$$

$$\mathcal{A}(\omega,u):=\Big{(}(\mathcal{F}(\omega\circ\mathcal{S}_{0}u))^{T},(\mathcal{F}(\omega\circ\mathcal{S}_{1}u))^{T},\cdots,(\mathcal{F}(\omega\circ\mathcal{S}_{J-1}u))^{T}\Big{)}^{T},$$

$$Y=1+XD,$$

$$\text{To find X and ω, leavevmodenobreak leavevmodenobreak s.t.leavevmodenobreak }|\mathcal{A}(\omega ,Y_l){|}^2\approx I_l,\text{leavevmode}\text{leavevmode}Y=1+XD,\text{leavevmode}\text{leavevmode}X\in \mathcal{X},$$

$$\begin{array}{rl}\mathrm{SP}:& \underset{\omega ,X,Y}{\min }\delta {\sum }_c\mathrm{TV}(X_c)+{\sum }_l\mathcal{G}(\mathcal{A}(\omega ,Y_l);I_l)+{\mathbb{I}}_{\mathcal{X}}(X)+{\mathbb{I}}_{\mathcal{W}}(\omega ),\text{leavevmode}\text{leavevmode}\text{s.t.}\text{leavevmode}\text{leavevmode}Y=1+XD,\end{array}$$

$$\Re (Y)=D_{r}X+1,$$

$$\begin{array}{rl}\mathrm{SPi}:& \underset{\omega ,X,Y}{\min }\delta {\sum }_c\mathrm{TV}(X_c)+{\sum }_l\mathcal{G}(\mathcal{A}(\omega ,Y_l);I_l)+{\mathbb{I}}_{\mathcal{X}}(X)+{\mathbb{I}}_{\mathcal{W}}(\omega ),\text{leavevmode}\text{s.t.}\text{leavevmode}\text{leavevmode}\Re (Y)=1+D_{r}X.\end{array}$$

$$X=(Y-1)\hat{D},$$

$$\begin{array}{rl}& \underset{\omega ,X,Y}{\min }\delta {\sum }_c\mathrm{TV}(X_c)+{\sum }_l\mathcal{G}(Z_l;I_l)+{\mathbb{I}}_{\mathcal{X}}(X)+{\mathbb{I}}_{\mathcal{W}}(\omega ),\\ & \text{leavevmode}\text{leavevmode}\text{s.t.}\text{leavevmode}\text{leavevmode}{Z}_l=\mathcal{A}(\omega ,Y_l),X=(Y-1)\hat{D}\text{leavevmode}\forall 0\le l\le L-1.\end{array}$$

$$\begin{array}{rl}& \underset{\Lambda ,\Gamma }{\max }\underset{\omega ,X,Y,Z}{\min }{\mathcal{L}}_{\lambda ,\beta }(\omega ,X,Y,Z,\Lambda ,\Gamma )\\ :=& \delta {\sum }_c\mathrm{TV}(X_c)+{\sum }_l\mathcal{G}(Z_l;I_l)+{\mathbb{I}}_{\mathcal{X}}(X)+{\mathbb{I}}_{\mathcal{W}}(\omega ),\\ & +{\sum }_l\left(\lambda \Re ⟨Z_l-\mathcal{A}(\omega ,Y_l),{\Lambda }_l⟩+\frac{\lambda }{2}\Vert Z_l-\mathcal{A}(\omega ,Y_l){\Vert }^2\right),\\ & +\beta \Re ⟨X-(Y-1)\hat{D},\Gamma ⟩+\frac{\beta }{2}\Vert X-(Y-1)\hat{D}{\Vert }^2,\end{array}$$

$$\begin{array}{rl}& {\omega }^{\star }:=\arg \min {\mathcal{L}}_{\lambda ,\beta }(\omega ,X,Y,Z,\Lambda ,\Gamma )\\ & =\arg {\min }_{\omega }\frac{1}{2}{\sum }_l\Vert Z_l+{\Lambda }_l-\mathcal{A}(\omega ,Y_l){\Vert }^2+{\mathbb{I}}_{\mathcal{W}}(\omega )\\ & =\arg {\min }_{\omega \in \mathcal{W}}\frac{1}{2}{\sum }_{l,j}\Vert {\mathcal{F}}^{\ast }(Z_{l,j}+{\Lambda }_{l,j})-\omega \circ {\mathcal{S}}_j{Y}_l{\Vert }^2.\end{array}$$

$$H(\omega ):=\mathrm{diag}({\sum }_{l,j}|{\mathcal{S}}_j{Y}_l{|}^2)\omega -{\sum }_{l,j}{\mathcal{F}}^{\ast }(Z_{l,j}+{\Lambda }_{l,j})\circ {\mathcal{S}}_j{Y}_l^{\ast }.$$

$$\begin{array}{r}{\omega }_{s+1}={\mathrm{Proj}}_{\mathcal{W}}\left({\omega }_s-\frac{H({\omega }_s)}{({\sum }_{l,j}|{\mathcal{S}}_j{Y}_l{|}^2)+{\gamma }_{1}1}\right)\text{leavevmode}\text{leavevmode}s=0,1,\cdots ,\end{array}$$

$$\begin{array}{rl}{X}^{\star }& :=\arg {\min }_X{\mathcal{L}}_{\lambda ,\beta }(\omega ,X,Y,Z,\Lambda ,\Gamma )\\ & =\arg {\min }_X{\sum }_c\left(\frac{\delta }{\beta }\mathrm{TV}(X_c)+\frac{1}{2}\Vert X_c-\Re ((Y-1)\hat{D}-\Gamma {)}_c{\Vert }^2\right)+{\mathbb{I}}_{\mathcal{X}}(X),\end{array}$$

$$X_c^{\star }=\max \{0,{\mathrm{Denoise}}_{\delta /\beta }(\Re (((Y-1)\hat{D}-\Gamma ){)}_c)\}\text{leavevmode}\forall 0\le c\le C-1,$$

$$\begin{array}{rl}& Y^{\star }:=\arg \underset{Y}{\min }{\mathcal{L}}_{\lambda ,\beta }(\omega ,X,Y,Z,\Lambda ,\Gamma )\\ & =\arg {\min }_Y\frac{\lambda }{2}{\sum }_l\Vert \mathcal{A}(\omega ,Y_l)-({\Lambda }_l+Z_l){\Vert }^2+\frac{\beta }{2}\Vert Y\hat{D}-(\Gamma +X+1\hat{D}){\Vert }^2+\frac{{\gamma }_2}{2}\Vert Y-Y_0{\Vert }^2\\ & =\arg {\min }_Y\frac{\lambda }{2}{\sum }_l\Vert \omega \circ {\mathcal{S}}_j{Y}_l-{\mathcal{F}}^{\ast }({\Lambda }_l+Z_l)\Vert +\frac{\beta }{2}\Vert Y\hat{D}-(\Gamma +X+1\hat{D}){\Vert }^2+\frac{{\gamma }_2}{2}\Vert Y-Y_0{\Vert }^2\\ & =\arg {\min }_Y\frac{\lambda }{2}{\sum }_l\Vert {\mathcal{S}}_j^T\omega \circ Y_l-{\mathcal{S}}_j^T{\mathcal{F}}^{\ast }({\Lambda }_l+Z_l)\Vert +\frac{\beta }{2}\Vert Y\hat{D}-(\Gamma +X+1\hat{D})){\Vert }^2+\frac{{\gamma }_2}{2}\Vert Y-Y_0{\Vert }^2,\\ & \end{array}$$

$$\begin{array}{r}\mathrm{diag}(\lambda {\sum }_j|{\mathcal{S}}_j^T\omega {|}^2+{\gamma }_2\mathbf{I})Y+\beta Y\hat{D}{\hat{D}}^{\ast }=\lambda Q+{\gamma }_2{Y}_0+\beta (\Gamma +X+1\hat{D}){\hat{D}}^{\ast },\end{array}$$

$$\begin{array}{r}\mathrm{diag}(\lambda {\sum }_j|{\mathcal{S}}_j^T\omega {|}^2+{\gamma }_2\mathbf{I})\hat{Y}+\beta \hat{Y}\mathcal{S}=(\lambda Q+{\gamma }_2{Y}_0+\beta (\Gamma +X+1\hat{D}){\hat{D}}^{\ast })V,\end{array}$$

$$Y^{\star }={\hat{Y}}^{\star }{V}^{\ast },$$

$${\hat{Y}}_l^{\star }=\frac{((\lambda Q+{\gamma }_2{Y}_0+\beta (\Gamma +X+1\hat{D}){\hat{D}}^{\ast })V{)}_l}{(\lambda {\sum }_j|{\mathcal{S}}_j^T\omega {|}^2+{\gamma }_{2}1)+\beta {\mathcal{S}}_{l,l}1}\text{leavevmode}\text{leavevmode}\forall \text{leavevmode}0\le l\le L-1.$$

$$Z^{\star }:=\arg \underset{Z}{\min }{\sum }_l\mathcal{G}(Z_l;I_l)+\frac{\lambda }{2}{\sum }_l\Vert Z_l-(\mathcal{A}(\omega ,Y_l)-{\Lambda }_l){\Vert }^2,$$

$$Z_l^{\star }=\frac{\sqrt{4(1+\lambda )I_l+{\lambda }^2|{\hat{Z}}_l{|}^2}+\lambda |{\hat{Z}}_l|}{2(1+\lambda )}\circ \mathrm{sign}({\hat{Z}}_l),$$

$$X=\Re (Y-1){\hat{D}}_r.$$

$$\begin{array}{r}\underset{\omega ,X,Y}{\min }\delta {\sum }_c\mathrm{TV}(X_c)+{\sum }_l\mathcal{G}(\mathcal{A}(\omega ,Y_l);I_l)+{\mathbb{I}}_{\mathcal{X}}(X)+{\mathbb{I}}_{\mathcal{W}}(\omega ),\text{leavevmode}\text{s.t.}\text{leavevmode}\text{leavevmode}X=\Re (Y-1){\hat{D}}_r.\end{array}$$

$$\begin{array}{rl}& \underset{\Lambda ,{\Gamma }_r}{\max }\underset{\omega ,X,Y,Z}{\min }{\tilde{\mathcal{L}}}_{\lambda ,\beta }(\omega ,X,Y,Z,\Lambda ,{\Gamma }_r):=\delta {\sum }_c\mathrm{TV}(X_c)+{\sum }_l\mathcal{G}(Z_l;I_l)+{\mathbb{I}}_{\mathcal{X}}(X)+{\mathbb{I}}_{\mathcal{W}}(\omega )\\ & +{\sum }_l\left(\lambda \Re ⟨Z_l-\mathcal{A}(\omega ,Y_l),{\Lambda }_l⟩+\frac{\lambda }{2}\Vert Z_l-\mathcal{A}(\omega ,Y_l){\Vert }^2\right)\\ & +\beta ⟨X-\Re (Y-1){\hat{D}}_r,{\Gamma }_r⟩+\frac{\beta }{2}\Vert X-\Re (Y-1){\hat{D}}_r{\Vert }^2.\end{array}$$

$$\begin{array}{rl}& X^{\star }:=\arg \min \frac{\delta }{\beta }{\sum }_c\mathrm{TV}(X_c)+{\mathbb{I}}_{\mathcal{X}}(X)+\frac{1}{2}\Vert X-(\Re (Y-1){\hat{D}}_r-{\Gamma }_r){\Vert }^2.\end{array}$$

$$X_c^{\star }=\max \{0,{\mathrm{Denoise}}_{\delta /\beta }((\Re (Y-1){\hat{D}}_r-{\Gamma }_r{)}_c)\}\text{leavevmode}\forall 0\le c\le C-1.$$

$$\begin{array}{rl}& Y^{\star }:=\arg {\min }_Y\frac{\lambda }{2}\sum _{j,l}\Vert {\mathcal{S}}_j^T\omega \circ Y_l-{\mathcal{S}}_j^T{\mathcal{F}}^{\ast }({\Lambda }_l+Z_l){\Vert }^2+\frac{\beta }{2}\Vert X-\Re (Y-1){\hat{D}}_r+{\Gamma }_r{\Vert }^2+\frac{{\hat{\gamma }}_2}{2}\Vert Y-Y_0{\Vert }^2,\end{array}$$

$$\{\begin{array}{ll}& \mathrm{diag}(\lambda {\sum }_j|{\mathcal{S}}_j^T\omega {|}^2+{\hat{\gamma }}_{2}1)\Re (Y)+\beta \Re (Y){\hat{D}}_r{\hat{D}}_r^T=\lambda \Re (Q)+{\hat{\gamma }}_2\Re (Y_0)+\beta ({\Gamma }_r+X+1{\hat{D}}_r){\hat{D}}_r^T\\ & \mathrm{diag}(\lambda {\sum }_j|{\mathcal{S}}_j^T\omega {|}^2+{\hat{\gamma }}_{2}1)\Im (Y)=\lambda \Im (Q)+{\hat{\gamma }}_2\Im (Y_0).\end{array}$$

$$\Im (Y)=\frac{\lambda \Im (Q)+{\hat{\gamma }}_2\Im (Y_0)}{\lambda {\sum }_j|{\mathcal{S}}_j^T\omega {|}^2+{\hat{\gamma }}_{2}1}.$$

Full text

\journalcode

J

1 Introduction

X-ray spectro-microscopy is a powerful technique to study the chemical and morphological structure of a material at high resolution. The contrast of the material under study is recorded as a function of photon energy, and this spectral absorption contrast can later be used to reveal details about its chemical, orbital or magnetic state [stohr2013nexafs, koningsberger1988x]. The idea is that because different chemical components interact differently with the beam at different energies, the composition map of a sample can be solved by using a measured reference of spectra (dictionary).

Compared to standard lens-based microscopy, x-ray ptychography can provide much finer spatial resolution, while also providing additional phase contrast of the sample [nellist1995resolution, chapman1996phase, rodenburg2004phase, rodenburg2007hard]. Ptychography is based on retrieving the phase of diffraction data recorded to a numerical aperture that is far larger than what x-ray optics can technically achieve. In ptychography, the probe (illumination) is almost never completely known, so a joint recovery problem (sample and probe) is typically considered, referred to as blind ptychography. Several algorithms to solve both standard and blind ptychography problems have been published in the literature to solve both standard and blind ptychography problems, which also consider a variety of additional experimental challenges [maiden2009improved, thibault2009probe, thibault2012maximum, wen2012, marchesini2013augmented, horstmeyer2015solving, hesse2015proximal, odstrvcil2018iterative, chang2018Blind].

As in standard spectro-microscopy, it is also possible to perform spectroscopic ptychography by recording diffraction data at different x-ray photon energies. In recent years, spectro-ptychography has become an increasingly popular chemical analysis technique [beckers2011chemical, maiden2013soft, hoppe2013high, shapiro2014chemical, farmand2017near, shi2016soft]. However, the standard methodology involves independent ptychographic reconstructions for each energy, followed by component analysis, i.e. spectral imaging analysis based on a known reference spectra or multivariate analysis [adams1986spectral, lerotic2004cluster, shapiro2014chemical, yu2018three]. More recently, a low rank constraint [vaswani2017low] for multi-channel samples was proposed together with a gradient descent algorithm with spectral initialization to recover the higher dimension phase retrieval problem (without component analysis). Another work proposed a hierarchical model with Gaussian-Wishart hierarchical prior and developed a variational expectation-maximization algorithm [liu2019low]. Also, matrix decomposition based low rank prior [chen2018low] has been exploited to reconstruct dynamic, time-varying targets in Fourier ptychographic imaging.

In this paper we propose a novel technique to solve the blind x-ray spectro-ptychography problem based on coupling the diffraction data from each photon energy and iteratively retrieve the chemical map of the sample. The proposed algorithm, referred to as SPA (Spectroscopic Ptychography with ADMM), works with both completely and partially known reference spectra. The method is designed using the Alternating Direction Method of Multipliers (ADMM) [glowinski1989augmented, chang2018Blind] framework, employing also total variation (TV) regularization [rudin1992nonlinear] on the chemical map. Compared with the standard two-step methods, the proposed joint reconstruction algorithm can generate much higher quality results without presenting the phase ambiguity problem inherent to two-step methods111 Since in the first step for ptychography reconstruction for the two-step method, each sample contrast map is recovered independently, that causes to different phase factors for different maps.. The simulation analysis shows the efficient convergence ratio of SPA and demonstrates the increased robustness of the method to large stepsizes, being able to retrieve features lost when using standard two-step methods. The algorithm is described and analyzed with and without TV regularization for both partially and completely known dictionary cases.

2 Spectroscopic ptychography model

Given $L$ different energies of x-rays going through a sample illuminated by a probe $\omega\in\mathbb{C}^{\bar{m}}$, a collection of phaseless intensities $\{I_{l}\}_{l=0}^{L-1}$ are measured in the far field, such that with Poisson fluctuation caused by photon counting, we have

[TABLE]

where $Y=(Y_{0},Y_{1},\ldots,Y_{L-1})\in\mathbb{C}^{N,L}$ is the sample contrast map for each X-ray energy, $\mathcal{A}(w,\cdot)$ is the forward operator for ptychography for a given probe $w$, $\mathrm{Poi}$ denotes the Poisson-noise contamination, and the notations $|\cdot|,(\cdot)^{2}$ denoting the pointwise absolute and square values of the vector respectively. Note that the probe $\omega$, and each column of contrast maps $\{Y_{l}\}$ are all 2D images, written as vectors by a lexicographical order. The relationship between the contrast map $Y_{l}$ observed by ptychography at each energy and an unknown sample elemental map $X$ made of $C$ elements is governed by the spectral contrast of each element, stored in a “dictionary” $D$ of known tabulated values.

Specifically, following similar notation as in [chang2018partially], the bilinear operator $\mathcal{A}:\mathbb{C}^{\bar{m}}\times\mathbb{C}^{N}\rightarrow\mathbb{C}^{m}$ is defined as

[TABLE]

where $\mathcal{F}$ denotes the discrete Fourier transform, $\circ$ denotes the Hadamard product (pointwise multiplication) of two vectors, and $\mathcal{S}_{j}\in\mathbb{R}^{\bar{m},m}$ is a binary matrix that defines a small window with the index $j$ and size $\bar{m}$ over the entire image $u$ (taking small patches out of the entire image).

For different energies, assuming that a spectrum dictionary $D\in\mathbb{C}^{C,L}$ (or its absorption part) is measured in advance, having $C$ components for different materials or particles, and given a sufficiently thin specimen, the sample contrast maps can be approximated by first-order Taylor expansion $\exp(XD)\approx\bm{1}+XD$ 222$\bm{1}$ denotes the matrix with all elements being one, of the same size of $Y$. $\exp(\cdot)$ denotes the pointwise exponentiation of the matrix. as

[TABLE]

with $X=(X_{0},X_{1},\ldots,X_{C-1})\in\mathbb{R}_{+}^{N,C}$ being the elemental thickness map of the sample (each column of thickness map denotes the thickness of each components in the object).

To determine the thickness map $X$, with a completely known spectrum $D$, one has to solve the following problem :

[TABLE]

with non-negative thickness constraint set $\mathscr{X}=\{X=(X_{n,c})\in\mathbb{R}^{N,C}:X_{n,c}\geq 0,\ 0\leq n\leq N-1,0\leq c\leq C-1\}.$ Letting the illumination be normalized, i.e. $\omega\in\mathscr{W}:=\{\omega\in\mathbb{C}^{\bar{m}}:\leavevmode\nobreak\ \|\omega\|=1\}$, the total variation regularized nonlinear optimization model can be established by assuming the piecewise smoothness of the thickness map as:

[TABLE]

where \mathcal{G}(z;f):=\tfrac{1}{2}\sum_{n=0}^{N-1}\big{(}|z_{n}|^{2}-f_{n}\log(|z_{n}|^{2})\big{)}\leavevmode\nobreak\ \forall z=(z_{0},z_{1},\cdots,z_{N-1})^{T}\in\mathbb{C}^{N}, $f=(f_{0},f_{1},\cdots,f_{N-1})^{T}\in\mathbb{R}^{N},$ derived from the maximum likelihood estimate of Poisson noised data [chang2016Total], $\mathrm{TV}$ denotes the standard total variation semi-norm [rudin1992nonlinear] to enforce the piecewise smooth structure of $X_{c}$ (the $c^{th}$ column of the mixing matrix (thickness map) denoted in (1)), and $\delta$ is a positive constant to balance the regularization and fitting terms (Bigger $\delta$ produces stronger smoothness). Here $\mathbb{I}_{\mathscr{X}}(X),\mathbb{I}_{\mathscr{W}}(\omega)$ denotes the indicator functions, with $\mathbb{I}_{\mathscr{X}}(X)=0\leavevmode\nobreak\ \text{if\leavevmode\nobreak\ }X\in\mathscr{X};\mathbb{I}_{\mathscr{X}}(X)=+\infty,\leavevmode\nobreak\ \text{\leavevmode\nobreak\ otherwise.}$ We remark that it is a convenient way to enforce hard constraints within an optimization formulation.

Experimentally, only the real values part (absorption) of the dictionary $D_{r}$ are measured. As $X$ is real-valued, we consider the following relation:

[TABLE]

where $D_{r}:=\Re(D)$. Similarly, we derive the following spectroscopic ptychography with incomplete dictionary (SPi) as

[TABLE]

** Remark 2.1**

Unlike solving the ptychography imaging independently for each energy, we use the low-rank structure of the recovery results of different energies, i.e. the rank of the matrix $Y-\bm{1}$ is no greater than that of $X$.

3 Proposed iterative algorithm

ADMM [glowinski1989augmented] is a powerful and flexible tool that has already been applied to both ptychography [wen2012, chang2018Blind] and phase tomography problems [chang2019iterative, aslan2019joint]. In this work we also adopt the ADMM framework to design an iterative joint spectro-ptychography solution. We construct the proposed algorithm considering both complete and incomplete dictionary cases.

3.1 Complete dictionary

Based on the spectro-ptychography model (3) for a complete dictionary of spectra, we design the proposed SPA algorithm (Spectroscopic Ptychography with ADMM) as described below.

Let $DD^{*}\in\mathbb{C}^{C,C}$ be non-singular, where $D^{*}$ denotes the Hermitian matrix of $D$, i.e. $D^{*}:=\mathrm{conj}(D^{T})$. Considering the constraint in (1), the following equivalent form can be derived

[TABLE]

with $\hat{D}:=D^{*}(DD^{*})^{-1}\in\mathbb{C}^{L,C}$. Accordingly, the following equivalent model can be considered, by introducing auxiliary variables $\{Z_{l}\}$:

[TABLE]

The benefit of considering (5) instead of (1) lies in the fact that: (i) the multiplier will be a low-dimensional variable, since the dimension of $Y$ is much higher than that of $X$, and (ii) the subproblem w.r.t. the variable $X$ can be more easily solved.

An equivalent saddle point problem for (6), based on the augmented Lagrangian, can be derived as:

[TABLE]

with the multipliers $\Lambda:=(\Lambda_{0},\cdots,\Lambda_{L-1})$ and $\Gamma$, where $\langle\cdot,\cdot\rangle$ denotes the inner product of vectors and matrices (trace norms) respectively.

The above saddle point problem can be solved by alternating minimization and update of the multipliers. We first define each sub-minimization problem. The $\omega-$subproblem, with the additional proximal term, can be expressed as

[TABLE]

The first-order gradient of the above least squares problem (without constraint) is given as

[TABLE]

Consequently, the projected gradient descent scheme with preconditioning can be derived as

[TABLE]

with parameter $\gamma_{1}>0$ in order to get rid of division by zeros, and $\mathrm{Proj}_{\mathscr{W}}(\omega):=\frac{\omega}{\|\omega\|}$. Here the parameter $\gamma_{1}$ heuristically set to be a small scalar related with the maximum value of the $\textstyle\sum\nolimits_{l,j}|\mathcal{S}_{j}Y_{l}|^{2}$, e.g. $\gamma_{1}=0.1\times\|\sum\nolimits_{l,j}|\mathcal{S}_{j}Y^{k}_{l}|^{2}\|_{\infty}.$

The $X-$subproblem can be expressed as

[TABLE]

where $(\cdot)_{c}$ denotes the $c^{th}$ column of a matrix. Since it is common practise to solve the total variation denoising problem by using a first-order operator-splitting algorithm [wu2011augmented, chambolle2011first], we directly give the approximate solution below:

[TABLE]

with $\mathrm{Denoise}_{\nu}(u_{0}):=\arg\min_{u}\nu\mathrm{TV}(u)+\tfrac{1}{2}\|u-u_{0}\|^{2}.$ Here we remark that to seek for exact solution with this positivity constraint, one may need more auxiliary variables and inner loops as [Chan2013]. For simplicity, we did not exactly solve the constraint problem, and instead, the above approximation is derived by the standard TV-L2 denoising without constraint and then a projection to the positivity constraint set.

The $Y-$subproblem, with additional proximal term $\tfrac{\gamma_{2}}{2}\|Y-Y_{0}\|^{2}$ and previous iterative solution $Y_{0}$, is expressed as

[TABLE]

where $\gamma_{2}$ is a positive scalar similarly to the parameter $\gamma_{1}$.

By calculating the first-order gradient of the above least squares problem, one has

[TABLE]

with identity operator $\mathbf{I},$ where $Q:=(Q_{0},Q_{1},\ldots,Q_{L-1})\in\mathbb{C}^{N,L}$ with $Q_{l}:=\textstyle\sum\nolimits_{j}\mathcal{S}_{j}^{T}(\omega^{*}\circ\mathcal{F}^{*}(\Lambda_{l}+Z_{l})),$ which is actually the Sylvester equation [Sylvester, Simoncini2016].

Assuming that the positive Hermitian $\hat{D}\hat{D}^{*}$ has the singular-value-decomposition (SVD) as $\hat{D}\hat{D}^{*}=V\mathscr{S}V^{*},$ with diagonal matrix (diagonal elements are singular values) $\mathscr{S}\in\mathbb{R}^{L,L}$ and unitary matrix $V\in\mathbb{C}^{L,L}$, by introducing $\hat{Y}:=YV,$ we derive:

[TABLE]

such that the closed form solution can be expressed as

[TABLE]

where $\hat{Y}^{\star}:=(\hat{Y}^{\star}_{0},\cdots,\hat{Y}^{\star}_{L-1})\in\mathbb{C}^{N,L},$ and

[TABLE]

For the $Z-$subproblem, we have [chang2016Total]:

[TABLE]

which gives

[TABLE]

with $\hat{Z}_{l}:=\mathcal{A}(\omega,Y_{l})-\Lambda_{l}$. Based on the above calculations and the update of multipliers, the baseline SPA algorithm is summarized in the appendix.

3.2 Incomplete dictionary

A complete dictionary is often difficult to obtain without an independent experiment prior to a spectro-ptychography experiment. The materials components and their chemical states are often not known in advance. Moreover the real part of the refractive index component is often not well known [henke1993x]. It is more difficult to measure because it requires an interferometric or reflectometry measurements rather than simple absorption spectroscopy measurements and reflectometry experiments are less commonly done. While the Cramers-Kronig relationships relate real and imaginary parts, the relationship requires an spectral measurement from 0 to infinity which is not possible to measure in finite time; Standard techniques to extend absorption spectra can only produce approximate values in the imaginary component. Hence, it is attractive in practice to provide a version working with the real part only.

In this subsection, we propose a variation of the SPA algorithm to solve the joint spectro-ptychography problem when the dictionary of spectra is only partially know, based on the model proposed in (4). By assuming that $D_{r}$ has full row-rank, i.e. $D_{r}D_{r}^{*}$ is non-singular, with known $\hat{D}_{r}:=D_{r}^{T}(D_{r}D_{r}^{T})^{-1}$ in advance, we have:

[TABLE]

Consequently, the following equivalent problem can be solved instead of (4):

[TABLE]

Similarly to the previous subsection, introducing the multiplier $\Gamma_{r}$ and auxiliary variable $Z$ yields the saddle point problem below, with the help of the augmented Lagrangian of (10):

[TABLE]

Below, we focus only on the differences with respect to the Algortihm 1. For the $X-$subproblem, we have

[TABLE]

Hence we get

[TABLE]

For the $Y-$subproblem with proximal terms $\|Y-Y_{0}\|^{2}$ , we have

[TABLE]

which results in the following equations w.r.t. the real and imaginary parts, respectively:

[TABLE]

Then, the real part of $Y$ can be solved by (8)-(9), while the imaginary part can be simply computed by

[TABLE]

The overall SPA algorithm with an incomplete dictionary is summarized in the appendix.

4 Simulation and Reconstruction Results

In the simulation analysis of the proposed algorithms we consider the synthetic thickness maps of three different materials, extracted from three (RGB) channels of a natural color image (after thresholding and shift, consisting of $256\times 256$ pixels), shown in Fig. 1. The real part of the spectrum dictionary (for two different materials, PMMA (polymethyl methacrylate) and PS (polystyrene), plus a constant w.r.t. ten different energies) was measured at the Advanced Light Source [yan2013accurate], and the imaginary part was derived using Kramers-Kronig relations [kronig1926theory]. Both real and imaginary part dictionaries are shown in Fig. 3.

The ptychography measurements are simulated with Poisson noise contamination, using a single grid scan at each energy.

A standard zone-plate with annular shape diffracts an image (Fig. 2 ) onto the detector after going through an order-sorting aperture. The zone plate annular aperture is mapped onto the detector by geometric magnification as outer-diameter (in microns) corresponds to an annular ring on the detector of dimension outer-diameter/ detector-pixel-size $\times$ detector-distance / focal-distance . The illumination probe (Fig. 2b) has beam width (full width-half max, FWHM) of 16 pixels. The relationship between pixels and actual dimensions in the Far-Field approximation is as follows: illumination pixel (real-space) dimensions = (wavelenght $\times$ detector-distance)/ detector-number-of-pixels $\times$ detector-pixel-size. The zone plate distance to the sample is assumed to be adjusted proportionally with energy to keep the sample in focus as usually done experimentally. We also assume that the detector distance is adjusted to maintain the spatial frequencies on the same detector pixels.

In order to evaluate the recovered results, the signal-to-noise ratio (SNR) in dB is used, which is denoted below:

[TABLE]

where $X_{g}$ corresponds to the ground truth thickness.

We compare the proposed iterative SPA algorithm with the standard two-step method. The two-step method consist of (i) perform ptychography reconstruction using a joint illumination, then (ii) perform spectroscopy analysis with a known dictionary (or known real part), and finally (iii) correct the phase ambiguity for different energies.

When assessing the performance of SPA, we consider both with and without regularization cases, where we simply set the regularization parameter $\delta=0$ and slightly adjust the algorithm by replacing Step 2 with

[TABLE]

for baseline SPA and

[TABLE]

for the incomplete dictionary case.

4.1 Reconstruction quality

The first simulation assesses the reconstruction quality achieved by the proposed SPA algorithm, compared with the two-step method, when using a scan stepsize of 32 pixels. Figs. 4-5 depict the reconstructed images when using complete and incomplete dictionaries, respectively. The SPA simulations are performed without and with regularization, in rows 2 and 3 of Figs. 4-5, respectively. Visually, we can see obvious artifacts in the recovered images when using the two-step method (first row of Figs. 4-5). Such artifacts are greatly enhanced when reconstructed using SPA. Specifically, clear improvements can be identified in the regions corresponding to the red and blue circles for all three materials in both Fig. 4 and 5. The SNRs of the recovery results parallel the visual analysis. For the completely known dictionary, the two-step method achieves a SNR of 14.0 dB for the above simulation, whereas SPA achieves 18.1 dB (no regularization), and 18.8 dB (regularization). In the partially known dictionary simulation, the SNRs are 13.8, 15.8, and 16.7 dBs, for the two-step method and SPA with no-regularization and regularization, respectively, which achieves a comparative gain of more than 2 dB, similarly as in the known dictionary case.

The phase ambiguity is an inherent problem of the two-step method that causes a loss in reconstruction accuracy. For example, for the simulation shown in Fig. 4, the SNR without phase correction is only 12.3 dB, reaching 14.0 dB after applying correction. Even when using effective phase correction post-process, SPA proves to be more efficient for the simulations performed: higher quality reconstruction are achieved overall, and there is no need to correct the phase ambiguity due to the iterative reconstruction exploiting the low-rank structure and positivity constraint of the thickness function.

4.2 Robustness and convergence

The following simulation assesses the robustness of the proposed algorithm when varying the scanning stepsizes. The quantitative results of this simulation are presented in Table 2. The results demonstrate the enhanced robustness of SPA when handling bigger stepsizes, compared to the reference two-step method, achieving up to 10 dB increase in SNR. To permit a better visual analysis, we provide the reconstruction results of the three algorithms with 40 pixels stepsize in Fig. 6. The figure highlights the dramatic improvement achieved by SPA compared with the standard two-step method when reconstructing low-redundancy ptychographic data. Specifically, we can see how the features within the blue and red circles are almost lost in the two-step reconstruction, while they can be clearly observed when reconstructing using SPA.

Generally speaking, to make the proposed algorithms work, the basic condition is to assume $D$ has full row rank such that $DD^{*}$ is non-singular. However, the performance should also rely on the similarity of spectral elements. Here we introduce a factor $s\in[0,0.5]$ to generate a new dictionary $D_{s}\in\mathbb{C}^{3\times 10}$ with $D_{s}[1,k]=(1-s)D[1,k]+sD[2,k],D_{s}[2,k]=(1-s)D[2,k]+sD[1,k]\forall 1\leq k\leq 10.$ Readily one knows that (1) $D_{0}=D$; and (2) the first two rows are exactly the same if $s=1/2$ ($D_{1/2}$ does not have full row rank). Please see Fig. 7 for the dictionaries with $s=0.1,0.3,$ and $0.45$. Therefore, the parameter $s$ can be used to control the similarity of the new spectral dictionary (bigger $s$ implies higher similarity). We test the impact of proposed algorithm SPA by the different similarity of spectral dictionaries, and please see reconstruction results in Fig. 9 with $s=0.1,0.3,0.45$, and SNR changes in Fig. 8 with $s\in\{0,0.1,0.2,0.3,0.4,0.45,0.475,0.49\}$. One readily knows that the quality of reconstruction results by proposed SPA decays as the spectral becomes similar (the parameter $s$ gets close to $0.5$). Hence, to get better reconstruction, one should design the experiments with little similarity in spectral dictionaries.

The last simulation depicts the error curve achieved by the SPA algorithm, shown in Fig. 10. The results demonstrate an steady decrease of the successive errors of the proposed algorithm, both with and without regularization.

5 Conclusions

This paper presents the first iterative spectroscopy ptychography solution. The proposed SPA algorithm is based on a novel spectro-ptychography model and it is constructed considering both a completely known and partially known dictionaries. Numerical simulations show that SPA produces more accurate results with clearer features compared with the standard two-step method. In the future, we will extend our work to thicker samples, where the first-order Taylor expansion is not sufficiently accurate. We also plan to investigate the use of Kramers-Kronig relationships [hirose2017use], explore the case using a completely unknown dictionary, and further provide software for real experimental data analysis.

6 Acknowledgments

This work of the first author was partially supported by National Natural Science Foundation of China (Nos.11871372, 11501413), Natural Science Foundation of Tianjin (No.18JCYBJC16600), 2017-Outstanding Young Innovation Team Cultivation Program (No.043-135202TD1703) and Innovation Project (No.043-135202XC1605) of Tianjin Normal University, Tianjin Young Backbone of Innovative Personnel Training Program and Program for Innovative Research Team in Universities of Tianjin (No.TD13-5078). This work was also partially funded by the Advanced Light Source and the Center for Advanced Mathematics for Energy Research Applications, a joint ASCR-BES funded project within the Office of Science, US Department of Energy, under contract number DOE-DE-AC03-76SF00098.

**Appendix

**

Algorithm 1: SPA

Initialization: Compute the SVD of $\hat{D}\hat{D}^{*}$ as $\hat{D}\hat{D}^{*}=V\mathscr{S}V^{*}.$ Set $Y^{0}:=\mathbf{1},\omega^{0}:=\mathcal{F}^{*}\left(\tfrac{\textstyle\sum\nolimits_{l,j}\sqrt{I_{l,j}}}{L\times J}\right),Z^{0}_{l}:=\mathcal{A}(\omega^{0},Y^{0}_{l}),\Lambda=\bm{0},$ and $\Gamma=\bm{0}.$ Set $k:=0.$

• 1.

Update the probe $\omega^{k+1}$ by one-step projected gradient descent method

[TABLE]

with $\gamma_{1}=0.1\times\|\sum\nolimits_{l,j}|\mathcal{S}_{j}Y^{k}_{l}|^{2}\|_{\infty}.$

• 2.

Update the thickness function $X^{k+1}$ by

[TABLE]

• 3.

Update $Y^{k+1}$ by

[TABLE]

with $Q^{k}:=(Q^{k}_{0},Q^{k}_{1},\ldots,Q^{k}_{L-1})\in\mathbb{C}^{N,L},$ $Q_{l}^{k}:=\sum\nolimits_{j}\mathcal{S}_{j}^{T}((\omega^{k+1})^{*}\circ\mathcal{F}^{*}(\Lambda_{l}^{k}+Z^{k}_{l})),$ and $\gamma_{2}^{k}=0.1\times\lambda\|\sum\nolimits_{j}|\mathcal{S}_{j}\omega^{k+1}|^{2}\|_{\infty}.$

• 4.

Update the auxiliary variable $Z^{k+1}$ by

[TABLE]

with $\hat{Z}^{k}_{l}:=\mathcal{A}(\omega^{k+1},Y^{k+1}_{l})-\Lambda^{k}_{l}$.

• 5.

Update the multipliers $\Lambda$ and $\Gamma$ by

[TABLE]

• 6.

When satisfying the stopping condition, output $X^{k+1}$ as the final thickness, otherwise, go to Step 1.

Algorithm 2: SPA with incomplete dictionary

Initialization: Compute the SVD of $\hat{D}_{r}\hat{D}_{r}^{*}$ as $\hat{D}_{r}\hat{D}_{r}^{*}=V_{r}\mathscr{S}_{r}V_{r}^{*}.$ Set $Y^{0}:=\mathbf{1},\omega^{0}:=\mathcal{F}^{*}\left(\tfrac{\textstyle\sum\nolimits_{l,j}\sqrt{I_{l,j}}}{L\times J}\right),Z^{0}_{l}:=\mathcal{A}(\omega^{0},Y^{0}_{l}),\Lambda=\bm{0},$ and $\Gamma_{r}=\bm{0}.$ Set $k:=0.$

• 1.

Update the probe $\omega^{k+1}$ as Step 1 of Algorithm 1.

• 2.

Update the thickness function $X^{k+1}$ by

[TABLE]

• 3.

Update $Y^{k+1}=Y_{r}^{k+1}+\mathrm{i}\times Y_{i}^{k+1}$ (Imaginary unit $\mathrm{i}:=\sqrt{-1}$) with real part update by

[TABLE]

with $\hat{\gamma}_{2}^{k}=0.1\times\lambda\|\sum\nolimits_{j}|\mathcal{S}_{j}\omega^{k+1}|^{2}\|_{\infty},$ and imaginary part update by

[TABLE]

• 4.

Update the auxiliary variable $Z^{k+1}$ as Step 4 of Algorithm 1.

• 5.

Update the multipliers $\Lambda$ and $\Gamma_{r}$ by

[TABLE]

• 6.

When satisfying the stopping condition, output $X^{k+1}$ as the final thickness, otherwise, go to Step 1.