Iterative X-ray spectroscopic ptychography

Chang, H.; Rong, Z.; Enfedaque, P.; Marchesini, S.

doi:10.1107/S1600576720006354

research papers

JOURNAL OF
APPLIED
CRYSTALLOGRAPHY

ISSN: 1600-5767

Volume 53| Part 4| August 2020| Pages 937-948

https://doi.org/10.1107/S1600576720006354

Open

access

Iterative X-ray spectroscopic ptychography¹

Huibin Chang,^a ^* Ziqin Rong,^a Pablo Enfedaque ^b and Stefano Marchesini ^b ^*

^aSchool of Mathematical Sciences, Tianjin Normal University, Tianjin, People's Republic of China, and ^bComputational Research Division, Lawrence Berkeley National Laboratory, Berkeley, CA, USA
^*Correspondence e-mail: changhuibin@gmail.com, smarchesini@lbl.gov

Edited by F. R. N. C. Maia, Uppsala University, Sweden (Received 10 September 2019; accepted 11 May 2020; online 8 July 2020)

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 alternating direction method of multipliers), a novel algorithm to iteratively solve the spectroscopic blind ptychography problem. First, a nonlinear spectro-ptychography model based on Poisson maximum likelihood is designed, and then the proposed method is constructed on the basis of fast iterative splitting operators. SPA can be used to retrieve spectral contrast when considering either 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 with standard state-of-the-art two-step methods. It is demonstrated how SPA can recover accurate chemical maps from Poisson-noised measurements, and its enhanced robustness when reconstructing reduced-redundancy ptychography data using large scanning step sizes is shown.

Keywords: spectromicroscopy; ptychography.

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 (Stöhr, 2013 ; Koningsberger & Prins, 1988 ). 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 measured reference spectra (a dictionary).

Compared with standard lens-based microscopy, X-ray ptychography can provide much finer spatial resolution, while also providing additional phase contrast of the sample (Nellist et al., 1995 ; Chapman, 1996 ; Rodenburg & Faulkner, 2004 ; Rodenburg et al., 2007 ). 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, which also consider a variety of additional experimental challenges (Maiden & Rodenburg, 2009 ; Thibault et al., 2009 ; Thibault & Guizar-Sicairos, 2012 ; Wen et al., 2012 ; Marchesini et al., 2013 ; Horstmeyer et al., 2015 ; Hesse et al., 2015 ; Odstrci et al., 2018 ; Chang et al., 2019a ).

As in standard spectro-microscopy, it is 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 (Beckers et al., 2011 ; Maiden et al., 2013 ; Hoppe et al., 2013 ; Shapiro et al., 2014 ; Farmand et al., 2017 ; Shi et al., 2016 ). 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 spectrum or multivariate analysis (Adams et al., 1986 ; Lerotic et al., 2004 ; Shapiro et al., 2014; Yu et al., 2018 ). More recently, a low-rank constraint (Vaswani et al., 2017 ) 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). Other work proposed a hierarchical model with a Gaussian–Wishart hierarchical prior and developed a variational expectation–maximization algorithm (Liu et al., 2019 ). Also, a matrix-decomposition-based low-rank prior (Chen et al., 2018 ) 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 retrieving 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) (Glowinski & Tallec, 1989 ; Chang et al., 2019a) framework, employing also total variation (TV) regularization (Rudin et al., 1992 ) 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 methods.² The simulation analysis shows the efficient convergence ratio of SPA and demonstrates the increased robustness of the method to large step sizes, 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

The main operators used in this section are given in Table 1.

Table 1
Main variables and operators defined in Section 2

Notations	Explanations
{ I_l} _{l = 0}^L-1	Measured intensities
$[X\in{\bb R}_{+}^{{N,C}}]$	Elemental thickness maps of the sample
$[D\in{\bb C}^{{C,L}}]$	Spectrum dictionary
$[Y\in{\bb C}^{{N,L}}]$	Sample spectral contrast maps
$[\omega\in{\bb C}^{{{\overline m}}}]$	Probe
$[{\cal S}_{j}\in{\bb R}^{{{\overline m},m}}]$	Binary matrix to take image patches
$[{\cal A}:{\bb C}^{{{\overline m}}} \times {\bb C}^{{N}}\rightarrow{\bb C}^{{m}}]$	Forward operator for ptychography
$[{{\cal G}}[{\cal A}(w,u),I_{{l}}]]$	Poisson likelihood estimation
$[D_{\rm r}\in{\bb R}^{{C,L}}]$	Real part of spectrum dictionary
$[{\scr X}]$	Non-negative thickness constraint
$[{\scr W}]$	Normalized constraint of the probe
TV	Total variation regularization

Given L different energies of X-rays going through a sample illuminated by a probe $[\omega \in {\bb C}^{{{\overline 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

$[{I_{l}} = {\rm Poi}[|{\cal A}(\omega,Y_{l})|^{2}]\quad {\rm for}\quad l = 0,1,\ldots,L-1.\eqno(1)]$

Here $[Y = (Y_{0},Y_{1},\ldots,Y_{{L-1}})\in{\bb C}^{{N,L}}]$ is the sample contrast map for each X-ray energy, $[{\cal A}(w,\cdot)]$ is the forward operator for ptychography for a given probe w, Poi denotes the Poisson-noise contamination, and the notations $[|\cdot|,(\cdot)^{2}]$ denote the pointwise absolute and square values of the vector, respectively. Note that the probe ω 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 values.

Specifically, following similar notation to Chang, Enfedaque et al. (2018 ), the bilinear operator $[{\cal A}:{\bb C}^{{{\overline m}}}\times{\bb C}^{{N}}\rightarrow{\bb C}^{{m}}]$ is defined as

$[\eqalignno{{\cal A}(\omega,u) & : = \Big\{[{\cal F}(\omega\circ{\cal S}_{0}u)]^{\rm T},[{\cal F}(\omega\circ{\cal S}_{1}u)]^{\rm T},\ldots,\cr & \quad\quad\ [{\cal F}(\omega\circ{\cal S}_{{J-1}}u)]^{\rm T}\Big\}^{\rm T},&(2)}]$

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

For different energies, assuming that a spectrum dictionary $[D\in{\bb 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)\simeq {\bf 1} + XD]$ ³ as

$[Y = {\bf 1} + XD, \eqno (3)]$

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

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

$[\eqalignno{& {\rm To \, find}\, X \, {\rm and} \, \omega, \, {\rm such \,that}\cr & |{\cal A}(\omega,Y_{l})|^{2} \simeq I_{l}, \quad Y = {\bf 1} + XD,\quad X \in {\scr X},& (4)}]$

with non-negative thickness constraint set $[{\scr X} = \{ X = (X_{{n,c}})\in]$ $[{\bb 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{\scr W}: = \{\omega\in{\bb C}^{{{\overline m}}}: ]$ $[ \|\omega\| = 1\}]$ , the total variation regularized nonlinear optimization model can be established by assuming the piecewise smoothness of the thickness map as

$[\eqalignno{ & {\rm SP}\!: {\mathop{\mathop {\rm min}\limits_{\omega,X,Y}}} \delta \textstyle\sum\limits _{c}{\rm TV}(X_{c})+ \sum\limits _{l}{\cal G}[{\cal A}(\omega,Y_{l})\semi I_{l}]+{\bb I}_{{{\scr X}}}(X)+{\bb I}_{{{\scr W}}}(\omega), \cr & {\rm such \,\,that} \,\, Y = {\bf 1} + XD.& (5)}]$

$[{\cal G}(z;f)]$ : = $[ \textstyle{{1}\over{2}} \sum _{{n = 0}}^{{N-1}} [|z_{n}|^{2}-f_{n}\log(|z_{n}|^{2})] ]$ $[\forall]$ z = $[ (z_{0},z_{1},\ldots,]$ $[z_{{N-1}})^{\rm T}\in]$ $[{\bb C}^{N}]$ , $[f = (f_{0},f_{1}, \ldots, f_{{N-1}})^{\rm T}\in{\bb R}^{N}]$ , derived from the maximum likelihood estimate of Poisson-noised data (Chang, Lou et al., 2018 ), TV denotes the standard total variation semi-norm (Rudin et al., 1992) to enforce the piecewise smooth structure of X_c [the cth column of the mixing matrix (thickness map) defined in equation (3)], and δ is a positive constant to balance the regularization and fitting terms (larger δ produces stronger smoothness). $[{\bb I}_{{{\scr X}}}(X)]$ , $[{\bb I}_{{{\scr W}}}(\omega)]$ denote the indicator functions, with $[{\bb I}_{{{\scr X}}}(X) = 0 ]$ if $[ X\in{\scr X}]$ ; $[{\bb I}_{{{\scr X}}}(X) = +\infty]$ otherwise. We remark that this is a convenient way to enforce hard constraints within an optimization formulation.

Experimentally, only the real-valued part (absorption) of the dictionary D_r is measured. As X is real valued, we consider the following relation:

$[\Re(Y) = D_{\rm r}X + {\bf 1},\eqno(6)]$

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

$[\eqalignno{& {\rm SPi}\!: \mathop{\mathop{\rm min}\limits_{\omega,X,Y }} \delta \textstyle\sum\limits _{c}{\rm TV}(X_{c}) + \textstyle\sum\limits _{l}{\cal G}[{\cal A}(\omega,Y_{l})\semi I_{l}]+ {\bb I}_{{{\scr X}}}(X)+{\bb I}_{{{\scr W}}}(\omega),\cr & {\rm such\,\,that}\,\, \Re(Y) = {\bf 1} + D_{\rm r}X.& (7)}]$

Remark: Rather than 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-{\bf 1}]$ is no greater than that of X.

3. Proposed iterative algorithm

ADMM (Glowinski & Tallec, 1989) is a powerful and flexible tool that has already been applied to both ptychography (Wen et al., 2012; Chang et al., 2019a) and phase tomography problems (Chang et al., 2019b ; Aslan et al., 2019 ). 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

On the basis of the spectro-ptychography model [equation (5)] for a complete dictionary of spectra, we design the proposed SPA algorithm as described below.

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

$[X = (Y-{\bf 1})\hat{D}, \eqno (8)]$

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

$[\eqalignno{& \mathop{\mathop {\rm min } \limits_{\omega,X,Y}} \delta \textstyle\sum\limits _{c}{\rm TV}(X_{c}) + \sum\limits _{l}{\cal G}(Z_{l}\semi I_{l}) + {\bb I}_{{{\scr X}}}(X) + {\bb I}_{{{\scr W}}}(\omega), \cr &{\rm such\,\,that}\,\, Z_{l} = {\cal A}(\omega,Y_{l}),\,\, X = (Y-{\bf 1})\hat{D} \,\, \forall\,\, 0\leq l\leq L-1. & (9)}]$

The benefits of considering equation (8) instead of equation (3) lie 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 with respect to the variable X can be more easily solved.

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

$[\eqalignno { \mathop{\mathop{\rm max}\limits_{\Lambda, \Gamma}}& \mathop{\mathop{\rm min}\limits_{\omega,X,Y,Z}} {\scr L}_{{\lambda,\beta}}(\omega,X,Y,Z,\Lambda,\Gamma)\cr &: = \delta \textstyle\sum\limits _{c}{\rm TV}(X_{c}) + \sum\limits _{l}{\cal G}(Z_{l}\semi I_{l}) + {\bb I}_{{{\scr X}}}(X) + {\bb I}_{{{\scr W}}}(\omega) \cr & + \textstyle\sum\limits _{l}\left[\lambda\Re\langle Z_{l}-{\cal A}(\omega,Y_{l}),\Lambda _{l}\rangle + ({{\lambda}/{2}})\| Z_{l} - {\cal A}(\omega,Y_{l})\|^{2}\right] \cr & + \beta\Re\langle X-(Y - {\bf 1})\hat{D}, \Gamma\rangle + \textstyle({{\beta}/{2}})\| X - (Y - {\bf 1})\hat{D}\|^{2},& (10)}]$

with the multipliers $[\Lambda: = (\Lambda _{0},\ldots,\Lambda _{{L-1}})]$ and Γ, where $[\langle\cdot,\cdot\rangle]$ denotes the inner product of two vectors (or trace norms for two matrices).

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

$[\eqalignno { \omega^{\star}& : = \arg\min{\scr L}_{{\lambda,\beta}}(\omega,X,Y,Z,\Lambda,\Gamma) \cr & = \arg\min\limits _{{\omega}}\textstyle{{1}\over{2}} \sum\limits _{{l}}\| Z_{l}+\Lambda _{l} - {\cal A}(\omega,Y_{l})\|^{2} + {\bb I}_{{{\scr W}}}(\omega) \cr & = \arg\min\limits _{{\omega\in{\scr W}}} \textstyle{{1}\over{2}} \sum\limits _{{l,j}}\|{\cal F}^{*}(Z_{{l,j}} + \Lambda _{{l,j}}) - \omega\circ{\cal S}_{j}Y_{l}\|^{2}.&(11)}]$

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

$[H(\omega): = {\rm diag}\bigg(\textstyle\sum\limits _{{l,j}}|{\cal S}_{j}Y_{l}|^{2}\bigg)\omega-\textstyle\sum\limits _{{l,j}}{\cal F}^{*}(Z_{{l,j}}+\Lambda _{{l,j}})\circ{\cal S}_{j}Y_{l}^{*}.\eqno(12)]$

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

$[{\omega _{{s+1}} = {\rm Proj}_{{{\scr W}}}\left[\omega _{s}-{{H(\omega _{s})} \over {(\textstyle\sum\nolimits _{{l,j}}|{\cal S}_{j}Y_{l}|^{2})+\gamma _{1}{\bf 1}}}\right]}, \quad s = 0,1,\ldots, \eqno (13)]$

with parameter $[\gamma _{1} \,\gt \, 0]$ in order to avoid division by zeros, and $[{\rm Proj}_{{{\scr W}}}(\omega): = {{\omega} /{\|\omega\|}}]$ . Here the parameter $[\gamma _{1}]$ is heuristically set to be a small scalar related to the maximum value of $[\textstyle\sum\nolimits _{{l,j}}|{\cal S}_{j}Y_{l}|^{2}]$ , e.g. $[\gamma _{1} = 0.1\|\textstyle\sum\nolimits _{{l,j}}|{\cal S}_{j}Y^{k}_{l}|^{2}\| _{\infty}]$ .

The X subproblem can be expressed as

$[\eqalignno { X^{\star}&: = \arg\min\limits _{X}{\scr L}_{{\lambda,\beta}}(\omega,X,Y,Z,\Lambda,\Gamma) \cr & = \arg\min\limits _{{X}}\textstyle\sum\limits _{c}\left\{({{\delta}/{\beta}}){\rm TV}(X_{c}) + {{1}\over{2}}\| X_{c}\!-\Re[(Y-{\bf 1})\hat{D}-\Gamma]_{c}\|^{2}\right\} \cr & \quad+ {\bb I}_{{{\scr X}}}(X),&(14)}]$

where $[(\cdot)_{c}]$ denotes the cth column of a matrix. Since it is common practice to solve the total variation denoising problem by using a first-order operator-splitting algorithm (Wu et al., 2011 ; Chambolle & Pock, 2011 ), we directly give the approximate solution below:

$[\eqalignno{& X_{c}^{\star} = \max \{ 0,{\rm Denoise}_{\delta/\beta}\{{\Re}[(Y-{\bf 1})\hat{D}-\Gamma]_{c}\}\}\cr & \forall \,\,0\leq c\leq C-1,&(15)}]$

with $[{\rm Denoise}_{{\nu}}(u_{0}): = \arg\min _{u}\nu\,{\rm TV}(u) + {{1}\over{2}}\| u-u_{0}\|^{2}]$ . Here we remark that, to seek the exact solution with this positivity constraint, one may need more auxiliary variables and inner loops (Chan et al., 2013 ). 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 $[({{\gamma _{2}}/{2}})\| Y-Y_{0}\|^{2}]$ and previous iterative solution Y₀, is expressed as

$[\eqalignno {&Y^{\star} : = \arg\min _{Y}{\scr L}_{{\lambda,\beta}}(\omega,X,Y,Z,\Lambda,\Gamma) \cr & = \arg\min\limits _{{Y}}\,({{\lambda}/{2}})\textstyle\sum\limits _{l}\|{\cal A}(\omega,Y_{l})-(\Lambda _{l} + Z_{l})\|^{2} \cr & \quad+ ({{\beta}/{2}})\| Y\hat{D}-(\Gamma+X+{\bf 1}\hat{D})\|^{2}+ ({{\gamma _{2}}/{2}})\| Y-Y_{0}\|^{2} \cr & = \arg\min\limits _{{Y}}\, ({{\lambda}/{2}})\textstyle\sum\limits _{l}\|\omega\circ{\cal S}_{j}Y_{l}-{\cal F}^{*}(\Lambda _{l}+Z_{l})\| \cr & \quad+ ({{\beta}/{2}})\| Y\hat{D}-(\Gamma+X+{\bf 1}\hat{D})\|^{2}+ ({{\gamma _{2}}/{2}})\| Y-Y_{0}\|^{2} \cr & = \arg\min\limits _{{Y}} \,({{\lambda}/{2}})\textstyle\sum\limits _{l}\|{\cal S}_{j}^{\rm T}\omega\circ Y_{l}-{\cal S}_{j}^{\rm T}{\cal F}^{*}(\Lambda _{l}+Z_{l})\|\cr & \quad +({{\beta}/{2}})\| Y\hat{D}-(\Gamma+X+{\bf 1}\hat{D})\|^{2}+ ({{\gamma _{2}}/{2}})\| Y-Y_{0}\|^{2},&(16)}]$

where $[\gamma _{2}]$ is a positive scalar similar to the parameter $[\gamma _{1}]$ .

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

$[\eqalignno{{\rm diag}&\bigg(\lambda\textstyle\sum\limits _{j}|{\cal S}^{T}_{j}\omega|^{2}+\gamma _{2}{\bf I}\bigg)Y+\beta Y\hat{D}\hat{D}^{*} \cr & = \lambda Q+\gamma _{2}Y_{0}+\beta(\Gamma+X+{\bf 1}\hat{D})\hat{D}^{*},& (17)}]$

with identity operator $[{\bf I}]$ , where $[Q: = (Q_{0},Q_{1},\ldots,Q_{{L-1}})\in]$ $[{\bb C}^{{N,L}}]$ with $[Q_{l}: = \textstyle\sum\nolimits _{j}{\cal S}_{j}^{\rm T}[\omega^{*}\circ{\cal F}^{*}(\Lambda _{l}+Z_{l})]]$ . Equation (17) is actually the Sylvester equation (Sylvester, 1884 ; Simoncini, 2016 ).

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

$[\eqalignno{ {\rm diag}&\bigg(\lambda\textstyle\sum\limits _{j}|{\cal S}^{\rm T}_{j}\omega|^{2}+\gamma _{2}{\bf I}\bigg)\hat{Y}+\beta\hat{Y}{\scr S} \cr & = [\lambda Q+\gamma _{2}Y_{0}+\beta(\Gamma+X+{\bf 1}\hat{D})\hat{D}^{*}]V,&(18)}]$

such that the closed-form solution can be expressed as

$[Y^{\star} = \hat{Y}^{\star}V^{*}, \eqno (19)]$

where $[\hat{Y}^{\star}: = (\hat{Y}^{\star}_{0},\ldots,\hat{Y}^{\star}_{{L-1}})\in{\bb C}^{{N,L}}]$ and

$[\hat{Y}^{\star}_{l} = {{\{[\lambda Q+\gamma _{2}Y_{0}+\beta(\Gamma+X+{\bf 1}\hat{D})\hat{D}^{*}]V\}_{l}}\over{(\lambda\textstyle\sum\nolimits _{j}|{\cal S}^{\rm T}_{j}\omega|^{2}+\gamma _{2}{\bf 1})+\beta{\scr S}_{{l,l}}{\bf 1}}}\,\, \forall \, 0\leq l\leq L-1.\eqno (20)]$

For the Z subproblem, we have (Chang, Lou et al., 2018)

$[Z^{\star}: = \arg\min _{{Z}}\textstyle\sum\limits _{l}{\cal G}(Z_{l}\semi I_{l})+({{\lambda}/{2}})\sum\limits _{l}\| Z_{l}-[{\cal A}(\omega,Y_{l})-\Lambda _{l}]\|^{2},\eqno(21)]$

which gives

$[Z_{l}^{\star} = {{\left[{4(1+\lambda)I_{l}+\lambda^{2}|\hat{Z}_{l}|^{2}} \right]^{1/2}\,+\,\lambda|\hat{Z}_{l}|}\over{2(1+\lambda)}}\circ{\rm sign}(\hat{Z}_{l}),\eqno(22)]$

with $[\hat{Z}_{l}: = {\cal A}(\omega,Y_{l})-\Lambda _{l}]$ . The above calculations and the update of the multipliers form the basis of the baseline SPA algorithm, summarized in Appendix A.

3.2. Incomplete dictionary

A complete dictionary is often difficult to obtain without an independent experiment prior to a spectro-ptychography experiment. The material's 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 (Henke et al., 1993 ). It is more difficult to measure because it requires interferometric or reflectometry measurements rather than simple absorption spectroscopy measurements, and reflectometry experiments are less commonly done. While the Kramers–Kronig relationships relate real and imaginary parts, the relationship requires a 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 equation (7). By assuming that D_r has full row rank, i.e. D_rD_r^* is non-singular, with $[\hat{D}_{\rm r}: =]$ D_r^T(D_rD_r^T)^-1 known in advance, we have

$[X = \Re(Y-{\bf 1})\hat{D}_{\rm r}.\eqno(23)]$

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

$[\eqalignno{& \mathop{\mathop{\rm min}\limits_{\omega,X,Y }} \, \delta \textstyle\sum\limits_{c}{\rm TV}(X_{c}) + \sum\limits_{l}{\cal G}[{\cal A}(\omega,Y_{l})\semi I_{l}]+{\bb I}_{{{\scr X}}}(X)+{\bb I}_{{{\scr W}}}(\omega),\cr & {\rm such\,\,that}\,\, X = \Re(Y-{\bf 1})\hat{D}_{\rm r}. & (24)}]$

Similarly to the previous subsection, introducing the multiplier $[\Gamma _{\rm r}]$ and auxiliary variable Z yields the saddle point problem below, with the help of the augmented Lagrangian of equation (24):

$[\eqalignno { & \mathop{\mathop{\rm max}\limits_{\Lambda, \Gamma_{\rm r}}} \mathop{\mathop{\rm min}\limits _{\omega,X,Y,Z}} {\widetilde{\scr L}}_{\lambda,\beta}(\omega,X,Y,Z,\Lambda,\Gamma _{\rm r}) \cr & : = \delta\textstyle\sum\limits _{c}{\rm TV}(X_{c})+\sum\limits _{l}{\cal G}(Z_{l}\semi I_{l})+{\bb I}_{{{\scr X}}}(X)+{\bb I}_{{{\scr W}}}(\omega) \cr & +\textstyle \sum\limits _{l}\left[\lambda\Re\langle Z_{l}-{\cal A}(\omega,Y_{l}),\Lambda _{l}\rangle+({{\lambda}/{2}})\| Z_{l}-{\cal A}(\omega,Y_{l})\|^{2}\right] \cr & + \beta\langle X-\Re(Y-{\bf 1})\hat{D}_{\rm r},\Gamma _{\rm r}\rangle+ ({{\beta}/{2}})\| X-\Re(Y-{\bf 1})\hat{D}_{\rm r}\|^{2}.&(25)}]$

Below, we focus only on the differences with respect to Algorithm 1 (Appendix A). For the X subproblem, we have

$[\eqalignno{ X^{\star}& : = \arg\min({{\delta}/{\beta}})\textstyle\sum\limits _{c}{\rm TV}(X_{c})+{\bb I}_{{{\scr X}}}(X)\cr & \quad\textstyle+{{1}\over{2}}\| X-[\Re(Y-{\bf 1})\hat{D}_{\rm r}-\Gamma _{\rm r}]\|^{2}.&(26)}]$

Hence we get

$[\eqalignno{ & X_{c}^{\star} = \max\{ 0,{\rm Denoise}_{{\delta/\beta}}\{[\Re(Y-{\bf 1})\hat{D}_{\rm r}-\Gamma _{\rm r}]_{c}\}\} \cr & \forall\,\, 0\leq c\leq C-1.&(27)}]$

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

$[\eqalignno{&Y^{\star} : = \arg\min\limits _{{Y}}({{\lambda}/{2}})\textstyle\sum\limits _{{j,l}}\|{\cal S}_{j}^{\rm T}\omega\circ Y_{l}-{\cal S}_{j}^{\rm T}{\cal F}^{*}(\Lambda _{l}+Z_{l})\|^{2}\cr & \quad +({{\beta}/{2}})\| X-\Re(Y-{\bf 1})\hat{D}_{\rm r}+\Gamma _{\rm r}\|^{2}+({{\hat{\gamma}_{2}}/{2}})\| Y-Y_{0}\|^{2},&(28)}]$

which results in the following equations with respect to the real and imaginary parts, respectively:

$[\eqalign{& {\rm diag}\bigg(\lambda\textstyle\sum\limits _{j}|{\cal S}^{\rm T}_{j}\omega|^{2}+\hat{\gamma}_{2}{\bf 1}\bigg)\Re(Y)+\beta\Re(Y)\hat{D}_{\rm r}\hat{D}_{\rm r}^{\rm T} = \lambda\Re(Q)\cr & \quad\quad+\hat{\gamma}_{2}\Re(Y_{0}) +\beta(\Gamma _{\rm r}+X+{\bf 1}\hat{D}_{\rm r})\hat{D}_{\rm r}^{\rm T} ,\cr & {\rm diag}\bigg(\lambda\textstyle\sum\limits _{j}|{\cal S}^{\rm T}_{j}\omega|^{2}+\hat{\gamma}_{2}{\bf 1}\bigg)\Im(Y) = \lambda\Im(Q)+\hat{\gamma}_{2}\Im(Y_{0}).}\eqno (29)]$

Then, the real part of Y can be solved by equations (19) and (20), while the imaginary part can be simply computed by

$[\Im(Y) = {{\lambda\Im(Q)+\hat{\gamma}_{2}\Im(Y_{0})} \over {\lambda\sum\nolimits _{j}|{\cal S}^{\rm T}_{j}\omega|^{2}+\hat{\gamma}_{2}{\bf 1}}}.\eqno(30)]$

The overall SPA algorithm with an incomplete dictionary is summarized in Appendix B.

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 × 256 pixels), shown in Fig. 1. The real part of the spectrum dictionary [for two different materials, (a) PMMA (polymethyl methacrylate) and (b) PS (polystyrene), plus (c) a constant with respect to ten different energies] was measured at the Advanced Light Source (Yan et al., 2013 ), and the imaginary part was derived using the Kramers–Kronig relations (Kronig, 1926 ). Both real and imaginary part dictionaries are shown in Fig. 2.

Figure 1
Truth for the three different materials: (a) PMMA, (b) PS and (c) constant.

Figure 2
Spectrum dictionaries [three different materials m0 (PMMA), m1 (PS) and m2 (constant)], with real part (a) and imaginary part (b). The x and y axes denote the order of the ten spectra and different energies, respectively.

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 illumination (Fig. 3a) onto the sample after the focused probe (Fig. 3b) has gone through an order-sorting aperture. The zone plate annular aperture is mapped onto the detector by geometric magnification as `outer diameter' (in µm) and corresponds to an annular ring on the detector of dimension (outer diameter/detector pixel size) × (detector distance/focal distance). The illumination probe [Fig. 3(b)] has a beam width (FWHM) of 16 pixels. The relationship between pixels and actual dimensions in the far-field approximation is as follows: illumination pixel (real-space) dimensions = (wavelength × detector distance)/(detector number of pixels × detector pixel size). The zone plate's distance from the sample is assumed to be adjusted proportionally with energy to keep the sample in focus, as is usually done experimentally. We also assume that the detector distance is adjusted to maintain the spatial frequencies on the same detector pixels.

Figure 3
(a) Lens (binary) and (b) probe (64 × 64 pixels).

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

$[{\rm SNR}(X,X_{\rm g}) = -10\log _{{10}}{\| X-X_{\rm g}\|^{2}}/{\| X\|^{2}}, \eqno (31)]$

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 consists of (i) performing ptychography reconstruction using a joint illumination, then (ii) performing spectroscopy analysis with a known dictionary (or known real part), and finally (iii) correcting 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

$[X^{{k+1}} = \max\{ 0,\Re[(Y^{k}-{\bf 1})\hat{D}-\Gamma^{k}]\} \eqno (32)]$

for baseline SPA and

$[X^{{k+1}} = \max\{ 0,\Re(Y^{k}-{\bf 1})\hat{D}_{\rm r}-\Gamma _{\rm r}^{k}\} \eqno (33)]$

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 step size of 32 pixels. Figs. 4 and 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, respectively, of Figs. 4 and 5. Visually, we can see obvious artifacts in the recovered images when using the two-step method (first row of Figs. 4 and 5). Such artifacts are greatly enhanced when the image is 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 Fig. 5. The SNRs of the recovery results parallel the visual analysis. For the completely known dictionary, the two-step method achieves an 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 dB, for the two-step method and SPA with no regularization and regularization, respectively, which achieves a comparative gain of more than 2 dB, similarly to the known dictionary case.

Figure 4
Reconstruction results using a known dictionary of spectra from Poisson-noised data with SNR = 29.2 dB and scan step size = 32. (a)–(c) Standard two-step method; (d)–(f) SPA without regularization; (g)–(i) SPA with TV regularization. The recovered probes are shown in the right column for the two-step method, SPA and SPA with TV (from top to bottom).

Figure 5
Reconstruction results using a partially known dictionary of spectra from Poisson-noised data with SNR = 29.2 dB and scan step size = 32. (a)–(c) Standard two-step method; (d)–(f) SPA without regularization; (g)–(i) SPA with TV regularization. The recovered probes are shown in the right column for the two-step method, SPA and SPA with TV (from top to bottom).

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 an effective phase correction post-process, SPA proves to be more efficient for the simulations performed: higher-quality reconstructions are achieved overall, and there is no need to correct the phase ambiguity because of 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 step sizes. The quantitative results of this simulation are presented in Table 2. The results demonstrate the enhanced robustness of SPA when handling larger step sizes, compared with 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 step size 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.

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

Step size		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

Figure 6
Reconstruction results using a known dictionary of spectra from Poisson-noised data with SNR = 29.0 dB and scan step size = 40. (a)–(c) Standard two-step method; (d)–(f) SPA without regularization; (g)–(i) SPA with TV regularization. The recovered probes are shown in the right column for the two-step method, SPA and SPA with TV (from top to bottom).

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{\bb 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]$ . We know that (1) D₀ = D and (2) the first two rows are exactly the same if s = 1/2 (D_1/2 does not have full row rank). (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 (larger s implies higher similarity). We test the impact of the proposed SPA algorithm by the different similarity of spectral dictionaries (see the SNR changes in Fig. 8 with s $[\in]$ {0, 0.1, 0.2, 0.3, 0.4, 0.45, 0.475, 0.49} and the reconstruction results in Fig. 9 with s = 0.1, 0.3, 0.45). We also know that the quality of reconstruction results obtained by the proposed SPA decays as the spectral dictionaries become similar (the parameter s gets close to 0.5). Hence, to get a better reconstruction, one should design the experiments with little similarity in spectral dictionaries.

Figure 7
Synthetic spectrum dictionaries D_s for different s [three different materials m0 (PMMA), m1 (PS) and m2 (constant)], for real part (left) and imaginary part (right). The x and y axes denote the spectrum and different energies, respectively.

Figure 8
SNR changes versus similarity parameter s.

Figure 9
Reconstruction results by proposed SPA using different dictionaries of spectra (s = 0.1, 0.3, 0.45), from Poisson-noised data with SNR = 29.2 dB and scan step size = 32.

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

Figure 10
Error $[{{\| x^{k}-X^{{k-1}}\|}/{\| X^{k}\|}}]$ variation versus iteration number for SPA without (a) and with TV regularization (b). The x and y axes denote iteration numbers and errors, respectively.

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 a partially known dictionary. 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 (Hirose et al., 2017 ), explore the case using a completely unknown dictionary and further provide software for real experimental data analysis.

APPENDIX A

Algorithm 1: SPA

(0) Initialization: compute the SVD of $[\hat{D}\hat{D}^{*}]$ as $[\hat{D}\hat{D}^{*} = V{\scr S}\,V^{*}]$ . Set $[Y^{0}: = {\bf 1}]$ , $[\omega^{0}: = {\cal F}^{*}[{{\textstyle\sum\nolimits _{{l,j}}({I_{{l,j}}} )^{1/2}}/({L \,J})}]]$ , $[Z^{0}_{{l}}: = {\cal A}(\omega^{0}]$ , Y⁰_l), $[\Lambda = {\bf 0}]$ and $[\Gamma = {\bf 0}]$ . Set k: = 0.

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

$[\eqalignno{\omega^{{k+1}}& = {\rm Proj}_{{{\scr W}}}\Bigg[{{\gamma _{1}\omega^{k}}\over{(\textstyle\sum\nolimits _{{l,j}}|{\cal S}_{j}Y^{k}_{l}|^{2})+\gamma _{1}{\bf 1}}} \cr & + {{\textstyle\sum\nolimits _{{l,j}}{\cal F}^{*}(Z^{k}_{{l,j}}+\Lambda^{k}_{{l,j}})\circ{\cal S}_{j}(Y^{k}_{l})^{*}}\over {(\textstyle\sum\nolimits _{{l,j}}|{\cal S}_{j}Y^{k}_{l}|^{2})+\gamma _{1}{\bf 1}}}\Bigg], &(34)}]$

with $[\gamma _{1} = 0.1\|\sum\nolimits _{{l,j}}|{\cal S}_{j}Y^{k}_{l}|^{2}\| _{\infty}]$ .

(2) Update the thickness function X^k+1 by

$[\eqalignno{ & X_{c}^{{k+1}} = \max\{ 0,{\rm Denoise}_{{\delta/\beta}}\{\Re[(Y^{k}-{\bf 1})\hat{D}-\Gamma^{k}]_{c}\}\}\cr & \forall \,\,0\leq c\leq C-1.& (35)}]$

(3) Update Y^k+1 by

$[\eqalign{ &Y^{{k+1}} = \hat{Y}^{{k+1}}V^{*},\cr& \hat{Y}^{{k+1}}_{l} = {{\{[\lambda Q^{k}+\gamma^{k}_{2}Y^{k}+\beta(\Gamma^{k}+X^{{k+1}}+{\bf 1}\hat{D})\hat{D}^{*}]V\}_{l}}\over{(\lambda\textstyle\sum\nolimits _{j}|{\cal S}^{\rm T}_{j}\omega^{{k+1}}|^{2}+\gamma^{k}_{2}{\bf 1})+\beta{\scr S}_{{l,l}}{\bf 1}}}\cr & \forall\,\, 0\leq l\leq L-1. }\eqno (36)]$

with $[Q^{k}: = (Q^{k}_{0},Q^{k}_{1},\ldots,Q^{k}_{{L-1}})\in{\bb C}^{{N,L}}]$ , $[Q_{l}^{k}: = \sum\nolimits _{j}{\cal S}_{j}^{\rm T}[(\omega^{{k+1}})^{*}\circ]$ $[{\cal F}^{*}(\Lambda _{l}^{k}+Z^{k}_{l})]]$ and $[\gamma _{2}^{k} = 0.1\lambda\|\sum\nolimits _{j}|{\cal S}_{j}\omega^{{k+1}}|^{2}\| _{\infty}]$ .

(4) Update the auxiliary variable Z^k+1 by

$[Z_{l}^{{k+1}} = {{\left[{4(1+\lambda)I_{l}+\lambda^{2}|\hat{Z}^{k}_{l}|^{2}} \right]^{1/2}\,+\,\lambda|\hat{Z}^{k}_{l}|}\over{2(1+\lambda)}} \, \circ{\rm sign}(\hat{Z}^{k}_{l}),\eqno(37)]$

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

(5) Update the multipliers Λ and Γ by

$[\eqalign { \Lambda _{l}& \leftarrow\Lambda _{l}+Z_{l}-{\cal A}(\omega,{\cal T}_{l}Y_{l})\,\, \forall \,\, 0\leq l\leq L-1, \cr \Gamma & \leftarrow\Gamma+X-(Y-{\bf 1})\hat{D}.}\eqno(38)]$

(6) When satisfying the stopping condition, output X^k+1 as the final thickness; otherwise, go to Step 1.

APPENDIX B

Algorithm 2: SPA with incomplete dictionary

(0) Initialization: compute the SVD of $[\hat{D}_{\rm r}\hat{D}_{\rm r}^{*}]$ as $[\hat{D}_{\rm r}\hat{D}_{\rm r}^{*} = V_{\rm r}{\scr S}_{\rm r}V_{\rm r}^{*}.]$ Set $[Y^{0}: = {\bf 1}]$ , $[\omega^{0}: = {\cal F}^{*}[{{\textstyle\sum\nolimits _{{l,j}}({I_{{l,j}}} )^{1/2}}/{(L\, J)}}]]$ , $[Z^{0}_{{l}}: = {\cal A}(\omega^{0}]$ , $[Y^{0}_{{l}}),\Lambda = {\bf 0}]$ and $[\Gamma _{\rm r} = {\bf 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

$[\eqalignno{& X_{c}^{{k+1}} = \max\{ 0,{\rm Denoise}_{{\delta/\beta}}\{[\Re(Y^{k}-{\bf 1})\hat{D}_{\rm r}-\Gamma^{k}_{\rm r}]_{c}\}\} \cr & \forall \,\,0\leq c\leq C-1.& (39)}]$

(3) Update Y^k+1 = Y_r^k+1+i Y_i^k+1 [imaginary unit i: = (-1 )^1/2] with the real part updated by

$[\eqalignno { & Y^{{k+1}}_{\rm r} = \hat{Y}_{\rm r}^{{k+1}}V_{\rm r}^{*}, \cr & (\hat{Y}_{\rm r}^{{k+1}})_{l} = {{\{[\lambda\Re(Q^{k})+\hat{\gamma}^{k}_{2}\Re(Y^{k})+\beta(\Gamma _{\rm r}^{k}+X^{{k+1}}+{\bf 1}\hat{D}_{\rm r})\hat{D}_{\rm r}^{*}]V\}_{l}}\over{(\lambda\textstyle\sum\nolimits _{j}|{\cal S}^{\rm T}_{j}\omega^{{k+1}}|^{2}+\gamma^{k}_{2}{\bf 1})+\beta{\scr S}_{{l,l}}{\bf 1}}}, \cr &\forall \,\, 0\leq l\leq L-1, & (40)}]$

with $[\hat{\gamma}_{2}^{k} = 0.1\lambda \|\textstyle\sum\nolimits _{j}|{\cal S}_{j}\omega^{{k+1}}|^{2}\| _{\infty}]$ , and the imaginary part updated by

$[Y^{{k+1}}_{i} = {{\lambda\Im(Q^{k})+\hat{\gamma}_{2}\Im(Y^{{k}})} \over {\lambda\sum\nolimits _{j}|{\cal S}^{\rm T}_{j}\omega^{{k+1}}|^{2}+\hat{\gamma}_{2}{\bf 1}}}.\eqno(41)]$

(4) Update the auxiliary variable Z^k+1 as Step 4 of Algorithm 1.

(5) Update the multipliers Λ and $[\Gamma _{\rm r}]$ by

$[\eqalign { & \Lambda _{l}\leftarrow\Lambda _{l}+Z_{l}-{\cal A}(\omega,{\cal T}_{l}Y_{l})\,\,\forall \,\,0\leq l\leq L-1, \cr & \Gamma _{\rm r}\leftarrow\Gamma _{\rm r}+X-\Re(Y-{\bf 1})\hat{D}_{\rm r}.}\eqno(42)]$

(6) When satisfying the stopping condition, output X^k+1 as the final thickness; otherwise, go to Step 1.

Footnotes

¹This article will form part of a virtual special issue of the journal on ptychography software and technical developments.

²In the first step of ptychography reconstruction for the two-step method, each sample contrast map is recovered independently, which results in different phase factors for different maps.

³ $[{\bf 1}]$ denotes the matrix with all elements being one, of the same size as Y. $[\exp(\cdot)]$ denotes the pointwise exponentiation of the matrix.

Funding information

The work of HC was partially supported by the National Natural Science Foundation of China (Nos. 11871372 and 11501413), the Natural Science Foundation of Tianjin (No. 18JCYBJC16600), the 2017 Outstanding Young Innovation Team Cultivation Program (No. 043-135202TD1703), the Innovation Project (No. 043-135202XC1605) of Tianjin Normal University, and the Tianjin Young Backbone of Innovative Personnel Training Program and Program for Innovative Research Teams 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 No. DOE-DE-AC03-76SF00098.

References

Adams, J. B., Smith, M. O. & Johnson, P. E. (1986). J. Geophys. Res. 91, 8098–8112. CrossRef Web of Science Google Scholar
Aslan, S., Nikitin, V., Ching, D. J., Bicer, T., Leyffer, S. & Gürsoy, D. (2019). Opt. Express, 27, 9128–9143. Web of Science CrossRef PubMed Google Scholar
Beckers, M., Senkbeil, T., Gorniak, T., Reese, M., Giewekemeyer, K., Gleber, S.-C., Salditt, T. & Rosenhahn, A. (2011). Phys. Rev. Lett. 107, 208101. Google Scholar
Chambolle, A. & Pock, T. (2011). J. Math. Imaging Vis. 40, 120–145. Web of Science CrossRef Google Scholar
Chan, R. H., Tao, M. & Yuan, X. (2013). SIAM J. Imaging Sci. 6, 680–697. Web of Science CrossRef Google Scholar
Chang, H., Enfedaque, P., Lou, Y. & Marchesini, S. (2018). Acta Cryst. A74, 157–169. Web of Science CrossRef IUCr Journals Google Scholar
Chang, H., Enfedaque, P. & Marchesini, S. (2019a). SIAM J. Imaging Sci. 1, 153–185. Web of Science CrossRef Google Scholar
Chang, H., Enfedaque, P. & Marchesini, S. (2019b). IEEE International Conference on Image Processing (ICIP), pp. 2931–2935. IEEE. Google Scholar
Chang, H., Lou, Y., Duan, Y. & Marchesini, S. (2018). SIAM J. Imaging Sci. 11, 24–55. Web of Science CrossRef Google Scholar
Chapman, H. N. (1996). Ultramicroscopy, 66, 153–172. CrossRef CAS Web of Science Google Scholar
Chen, Z., Jagatap, G., Nayer, S., Hegde, C. & Vaswani, N. (2018). IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Vol. I, pp. 6538–6542. IEEE. Google Scholar
Farmand, M., Celestre, R., Denes, P., Kilcoyne, A. D., Marchesini, S., Padmore, H., Tyliszczak, T., Warwick, T., Shi, X., Lee, J., Yu, Y., Cabana, J., Joseph, J., Krishnan, H., Perciano, T., Maia, F. R. N. C. & Shapiro, D. A. (2017). Appl. Phys. Lett. 110, 063101. Google Scholar
Glowinski, R. & Le Tallec, P. (1989). Augmented Lagrangian and Operator-Splitting Methods in Nonlinear Mechanics. Philadelphia: SIAM. Google Scholar
Henke, B. L., Gullikson, E. M. & Davis, J. C. (1993). At. Data Nucl. Data Tables, 54, 181–342. CrossRef CAS Web of Science Google Scholar
Hesse, R., Luke, D. R., Sabach, S. & Tam, M. K. (2015). SIAM J. Imaging Sci. 8, 426–457. Web of Science CrossRef Google Scholar
Hirose, M., Shimomura, K., Burdet, N. & Takahashi, Y. (2017). Opt. Express, 25, 8593–8603. Web of Science CrossRef CAS PubMed Google Scholar
Hoppe, R., Reinhardt, J., Hofmann, G., Patommel, J., Grunwaldt, J.-D., Damsgaard, C. D., Wellenreuther, G., Falkenberg, G. & Schroer, C. (2013). Appl. Phys. Lett. 102, 203104. Google Scholar
Horstmeyer, R., Chen, R. Y., Ou, X., Ames, B., Tropp, J. A. & Yang, C. (2015). New J. Phys. 17, 053044. Google Scholar
Koningsberger, D. & Prins, R. (1988). X-ray Absorption: Principles, Applications, Techniques of EXAFS, SEXAFS and XANES. Chemical Analysis. Wiley-Interscience. Google Scholar
Kronig, R. de L. (1926). J. Opt. Soc. Am. 12, 547–557. CAS Google Scholar
Lerotic, M., Jacobsen, C., Schäfer, T. & Vogt, S. (2004). Ultramicroscopy, 100, 35–57. Web of Science CrossRef PubMed CAS Google Scholar
Liu, K., Wang, J., Xing, Z., Yang, L. & Fang, J. (2019). IEEE Access, 7, 5642–5648. Web of Science CrossRef Google Scholar
Maiden, A., Morrison, G., Kaulich, B., Gianoncelli, A. & Rodenburg, J. (2013). Nat. Commun. 4, 1669. Google Scholar
Maiden, A. M. & Rodenburg, J. M. (2009). Ultramicroscopy, 109, 1256–1262. Web of Science CrossRef PubMed CAS Google Scholar
Marchesini, S., Schirotzek, A., Yang, C., Wu, H.-T. & Maia, F. (2013). Inverse Probl. 29, 115009. Google Scholar
Nellist, P., McCallum, B. & Rodenburg, J. (1995). Nature, 374, 630–632. CrossRef CAS Web of Science Google Scholar
Odstrci, M., Menzel, A. & Guizar-Sicairos, M. (2018). Opt. Express, 26, 3108–3123. Web of Science PubMed Google Scholar
Rodenburg, J. M. & Faulkner, H. M. (2004). Appl. Phys. Lett. 85, 4795–4797. Web of Science CrossRef CAS Google Scholar
Rodenburg, J., Hurst, A., Cullis, A., Dobson, B., Pfeiffer, F., Bunk, O., David, C., Jefimovs, K. & Johnson, I. (2007). Phys. Rev. Lett. 98, 034801. Google Scholar
Rudin, L. I., Osher, S. & Fatemi, E. (1992). Physica D, 60, 259–268. CrossRef Web of Science Google Scholar
Shapiro, D. A., Yu, Y.-S., Tyliszczak, T., Cabana, J., Celestre, R., Chao, W., Kaznatcheev, K., Kilcoyne, A. D., Maia, F., Marchesini, S., Meng, Y. S., Warwick, T., Yang, L. L. & Padmore, H. A. (2014). Nat. Photon. 8, 765–769. Web of Science CrossRef CAS Google Scholar
Shi, X., Fischer, P., Neu, V., Elefant, D., Lee, J., Shapiro, D., Farmand, M., Tyliszczak, T., Shiu, H.-W., Marchesini, S., Roy, S. & Kevan, S. D. (2016). Appl. Phys. Lett. 108, 094103. Google Scholar
Simoncini, V. (2016). SIAM Rev. pp. 377–441. Web of Science CrossRef Google Scholar
Stöhr, J. (2013). NEXAFS Spectroscopy, Springer Series in Surface Sciences, Vol. 25. Berlin, Heidelberg: Springer Science & Business Media. Google Scholar
Sylvester, J. J. (1884). C. R. Acad. Sci. Paris, 99, 67–71. Google Scholar
Thibault, P., Dierolf, M., Bunk, O., Menzel, A. & Pfeiffer, F. (2009). Ultramicroscopy, 109, 338–343. Web of Science CrossRef PubMed CAS Google Scholar
Thibault, P. & Guizar-Sicairos, M. (2012). New J. Phys. 14, 063004. Google Scholar
Vaswani, N., Nayer, S. & Eldar, Y. C. (2017). IEEE Trans. Signal Process. 65, 4059–4074. Web of Science CrossRef Google Scholar
Wen, Z., Yang, C., Liu, X. & Marchesini, S. (2012). Inverse Probl. 28, 115010. Google Scholar
Wu, C., Zhang, J. & Tai, X.-C. (2011). Inverse Probl. 5, 237–261. Web of Science CrossRef Google Scholar
Yan, H., Wang, C., McCarn, A. R. & Ade, H. (2013). Phys. Rev. Lett. 110, 177401. Google Scholar
Yu, Y.-S., Farmand, M., Kim, C., Liu, Y., Grey, C. P., Strobridge, F. C., Tyliszczak, T., Celestre, R., Denes, P., Joseph, J., Krishnan, H., Maia, F. R. N. C., Kilcoyne, A. L. D., Marchesini, S., Leite, T. P. C., Warwick, T., Padmore, H., Cabana, J. & Shapiro, D. A. (2018). Nat. Commun. 9, 921. Google Scholar