Table Of ContentLBNL-1003977
SHARP: a distributed, GPU-based ptychographic solver
Stefano Marchesini,1 Hari Krishnan,1 Benedikt J. Daurer,2 David A.
Shapiro,1 Talita Perciano,1 James A. Sethian,1 and Filipe R.N.C. Maia2
1Lawrence Berkeley National Laboratory, Berkeley, CA, USA
2Uppsala University, Uppsala, Sweden
(Dated: June 21, 2016)
Everbrighterlightsources,fastparalleldetectors,andadvancesinphaseretrievalmethods,have
made ptychography a practical and popular imaging technique. Compared to previous techniques,
ptychographyprovidessuperiorrobustnessandresolutionattheexpenseofmoreadvancedandtime
consuming data analysis. By taking advantage of massively parallel architectures, high-throughput
processing can expedite this analysis and provide microscopists with immediate feedback. These
6
advancesallowreal-timeimagingatwavelengthlimitedresolution,coupledwithalargefieldofview.
1
Here, we introduce a set of algorithmic and computational methodologies used at the Advanced
0
2 Light Source, and DOE light sources packaged as a CUDA based software environment named
SHARP(http://camera.lbl.gov/sharp),aimedatprovidingstate-of-the-arthigh-throughput
n
ptychography reconstructions for the coming era of diffraction limited light sources.
u
J
0 I. INTRODUCTION it has been used in a large array of applications, and
2 shown to be a remarkably robust technique for the char-
acterization of nano materials. A few software imple-
] Reconstructing the 3D map of the scattering poten-
n mentations of the reconstruction algorithm exist such as
tial of a sample from measurements of its far-field scat-
a ptypy (http://ptycho.github.io/ptypy/) and PtychoLib
- tering patterns is an important problem. It arises in a
[24], and a repository for sharing experimental data has
a variety of fields, including optics [1, 2], astronomy [3],
t been established [25].
a X-ray crystallography [4], tomography [5], holography
Ptychography can be used to obtain large high-
d [6, 7] and electron microscopy [8]. As such it has been
. resolution images. It combines the large field of view
s a subject of study for applied mathematicians for over
of a scanning transmission microscope with the resolu-
c a century. The fundamental problem consists of find-
i tion of scattering measurements. In a scanning trans-
s ing the correct phases that go along with the measured
y intensities, such that together they can be Fourier trans- mission microscope, operated in transmission mode, a
h focused beam is rastered across a sample, and the total
formed into the real-space image of the sample. To help
p transmittedintensityisrecordedforeachbeamposition.
recover the correct phases from intensity measurements
[ The pixel positions of the image obtained correspond to
a range of experimental techniques have been proposed
the beam positions used during the scan, and the value
2 along the years, such as interferometry/holography [6],
v of the pixel to the intensity transmitted at that posi-
randomphasemasks[7,9,10],gratings[11]. Avarietyof
8 tion. This limits the resolution of the image to the size
numerical techniques have also been recently developed,
4 of the impinging beam, which is typically limited by the
for example by approximating the problem as a matrix
4 quality of focusing optics and work distance constraints.
1 completion problem [12], or by other convex relaxations
In ptychography, instead of only using the total trans-
0 [13] tractable by semidefinite programming.
mittedintensity, onetypicallyrecordsthedistributionof
.
2 Since its first demonstration [14], progress has been that intensity in the far-field, i.e. the scattering pattern
0 made in solving the phase problem for a single diffrac-
produced by the interaction of the illumination with the
6 tionpatternrecordedfromanon-periodicobject,includ-
sample. Thediffractedsignalcontainsinformationabout
1
ingthedynamicupdateofthesupport[15]andavariety
features much smaller than the size of the x-ray beam,
:
v of projection algortihms [16–18]. Such methods, referred makingitpossibletoachievehigherresolutionsthanwith
i to as coherent diffractive imaging (CDI), attempt to re-
X the scanning techniques. The downside of having to use
cover the complete complex-valued scattering potential
the intensities is that one now has to retrieve the corre-
r
or electron density, and the complex exit wavefront scat-
a spondingphasestobeabletoreconstructanimageofthe
teredfromtheobject,providingphasecontrastaswellas
sample,whichismadeevenmorechallengingbythepres-
a way to overcome depth-of-focus limitations of regular
ence of noise, experimental uncertainties, and perturba-
optical systems.
tionsoftheexperimentalgeometry. Whileitisadifficult
Ptychography, a relatively recent technique, provides problem, it is usually tractable by making use of the re-
the unprecedented capability of imaging macroscopic dundancyinherentinobtainingdiffractionpatternsfrom
specimens in 3D and attain wavelength limited resolu- overlapping regions of the sample. This redundancy also
tion along with chemical specificity [19]. Ptychography permitsthetechniquetoovercomethelackofseveralex-
was proposed in 1969 [20, 21], and later experimentally perimental parameters and measurement uncertainties.
demonstrated [22, 23], with the aim of improving the Forexample, therearemethodstorecoverunknownillu-
resolution in x-ray and electron microscopy. Since then minations [26–29]. As a testament to their success these
2
methods are even used as a way of characterizing high a(1)(q) a(K)(q)
qualityx-rayoptics[30–32],thewavefrontofx-raylasers p p
[33] and EUV lithography tools [34]. y y
p
Ptychographical phasing is a non-linear optimization x
problem [35] still containing many open questions [36].
Several strategies, such as Alternating Directions [37],
projections, gradient [35], conjugate gradient, Newton
[38–40], spectral methods[36, 41] and Monte-carlo [42], x
havebeenproposedtohandlesituationswhenbothsam-
pleandpositions[35,41–43],areunknownparametersin
high dimensions, and to handle experimental situations
FIG. 1. Experimental geometry in ptychography: an un-
such as accounting for noise variance [39, 44], partial co-
known sample with transmission ψ(r) is rastered through an
herence [3, 45–47][41, 46, 48], background[38, 41, 49, 50]
illuminating beam ω(r), and a sequence of diffraction mea-
or vibrations. surements I = |a (q)|2 are recorded on an area detector
(i) (i)
Here, we describe an algorithm approach and software with pixel coordinates p at a distance z from the sample.
D
environment SHARP (Scalable Hetereogeneous Adap-
tiveReal-timePtychography)thatenableshighthrough-
put streaming analysis using computationally efficient ordinates
phase retrieval algorithms. The high performance com-
putational back-end written in C/CUDA and imple- q =kout kin
−
mentedforNVIDIAGPUarchitecuresishiddenfromthe 1 (cid:18) (cid:19)
microscopist, butcanbeaccessedandadaptedtopartic- = (px,py,zD) (0,0,1)
λ √p2+p2+z2 −
x y D
ular needs by using a python interface or by modifying
1
the source code. (p ,p ,0)
x y
Using SHARP we have built an intuitive graphical (cid:39)λzD
user interface that provides visual feedback, of both the where k = (0,0,k) and k = k p are the incident
recordeddiffractiondataaswellasthereconstructedim- in out |p|
and scattered wave vectors that satisfy k = k =
ages, throughout the data aquisition and reconstruction | in| | out|
k = 1/λ, and λ is the wavelength. With a distance p
processes at the Advanced Light Source (ALS). m
from the center to the edge of the detector, the diffrac-
We use a mathematical formulation of ptychography
tion limited resolution (half-period) of the microscope is
which was first introduced in [36–38, 40, 41]. given by the lengthscale r = λzD. As a consequence, the
2pm
coordinates in reciprocal and real space are defined as
II. SHARP SOFTWARE ENVIRONMENT q =(cid:0) µ , ν (cid:1), µ,ν 0,...,m 1
mr mr ∈{ − }
and
A. Forward model
r =(rµ,rν), µ,ν 0,...,m 1 ,
∈{ − }
In a ptychography experiment (see Fig. 1), one per- x =(rµ(cid:48),rν(cid:48)), µ(cid:48),ν(cid:48) 0,...,n m .
(i)
forms a series of diffraction measurement as a sample is ∈{ − }
rastered across an x-ray, electron or visible light beam. Whilex istypicallyrasteredonacoarsergrid,r+x
(i) (i)
The illumination is formed by an x-ray optic such as a spans a finer grid of dimension n n.
zone-plate. Themeasurementisperformedbybrieflyex- Inotherwords,weassumethata×sequenceofK diffrac-
posinganareadetectorsuchasaCCDwhichrecordsthe tion intensity patterns (q) are collected as the posi-
(i)
scattered photons. tion of the object is rasItered on the position x . The
(i)
In a discrete setting, a two-dimensional small beam (cid:112)
simple transform a = (q) is a variance stabiliz-
(i) (i)
withdistributionw(r)ofdimensionm m illuminates I
x y ingtransformforPoissonnoise[51,52]. Therelationship
×
a subregion positioned at x (referred to as frame) of
(i) among the amplitude a (q), the illumination function
(i)
anunknownobjectofinterestψ(r)ofdimensionn n .
x y w(r) and an unknown object ψ(r) to be estimated can
×
Here 0 < m < n, i = 1,...,K and K is the total num-
be expressed as follows:
ber of frames (also referred to as “views” in the litera-
ture). For simplicity we consider square matrices. Gen-
(cid:12) (cid:12)
eralization to non-square matrices is straightforward but a (q)=(cid:12) w(r)ψ(r+x )(cid:12) (1)
(i) (i)
requires more indices and complicates notation. F
The pixel coordinates on a detector placed at a and isthetwo-dimensionaldiscreteFouriertransform,
F
distance z from the sample are described as p =
D (cid:88)
(p ,p ,z ). Underfar-fieldandparaxialapproximations ( f)(q)= √1 e2πiq·rf(r). (2)
x y D F m2
the pixel coordinates are related to reciprocal space co- r
3
where the sum over r is given on all the indices m m B. Phase retrieval
×
of r. We define an operator T , that extracts a frame
(i)
out of an image ψ, and build the illumination operator Projection operators form the basis of every iterative
Q(i), which scales the extracted frame point-wise by the projection and projected gradient algorithms are imple-
illumination function w: mentedinSHARPandaccessiblethroughalibrary. The
projectionP ensuresthattheframesz matchtheexper-
a
Q(i)[ψ](r)=w(r)ψ(r+x(i)), iment, that is, they satisfy Eq. (4), and is referred to as
=w(r)T [ψ](r), data projector:
(i)
=z(i)(r). P z =F∗ Fz a (7)
a
Fz
With the operator Q, eq. (1) can be represented com- | |
pactly as: whiletheprojectionP ontotherangeofQ(seeFig. 2):
Q
(cid:40)
a= Fz , P =Q(Q∗Q)−1Q∗ (8)
a= FQψ∨ , or | | (3) Q
| | z =Qψ∨,
ensures that overlapping frames z are consistent with
each other and satisfy Eq. (5).
wherethesuperscriptψ∨denotesthelinearizedversionof
The projector P is relatively robust to Poisson noise
theimage(thesuperscriptwillbeomittedforsimplicity), a
[51],butweightingfactorstoaccountfornoisypixelscan
and more explicitely as:
be easily added [40].
a∈RKm2 (cid:12)(cid:12)F∈CKm2×Km2 z∈CKm2 (cid:12)(cid:12) Usingrelationship(6),wecanupdatetheimageψfrom
(cid:12) (cid:12) w and frames z:
a (cid:12)F ... 0 z (cid:12)
(1) (cid:12) (1) (cid:12)
... =(cid:12)(cid:12)(cid:12) ... ... ... ... (cid:12)(cid:12)(cid:12), (4) ψ ←QQ∗∗Qz (9)
a (cid:12) 0 ... F z (cid:12) or the illumination w from an image ψ and frames z
(K) (cid:12) (K) (cid:12)
(cid:12) (cid:12) [26, 27] multiplying (Eq. 6) on the left by diagTψ¯ and
(cid:12) (cid:12)
solving for w:
z∈CKm2 Q∈CKm2×n2, ψ∈Cn2
z(1) diag(w)T(1) ψ1 w ←S∗dS∗iaTg|(ψT|ψ2¯)z, (10)
... = ... ... . (5) S where ψ¯ denotes the complex conjugate of ψ. See [29]
z(K) diag(w)T(K) ψn2 for alternative updates, and [28] for convergence theory
behind a similar blockwise optimization strategy. Sev-
where z are K frames extracted from the object ψ and eral possible pathologies need to be accounted for when
multiplied by the illumination function w, and F is the updating both ψ and w:
associated 2D DFT matrix when we write everything in
the stacked form [41]. When both the sample and the il- Combined drift of the illumination and the image
•
luminationareunknown,wecanexpresstherelationship in real space. Drift is eliminated by keeping the
(Eq. 5) between the image ψ, the illumination w, and illumination in the center of the frame by comput-
the frames z in two forms: ingthecenterofmassandcorrectingfordriftsafter
every update of the illumination.
z =Qψ =diag(Sw)Tψ =diag(Tψ)Sw (6)
Fourier space drifts and grid pathologies are su-
•
whereS RKm2×m2 denotestheoperatorthatreplicates pressedbyenforcingeithertheabsolutevalueaw =
∈ w or support m of the Fourier transform of
the illumination w into a stack of K frames, since Qψ = 0 w
|F |
the unknown illumination w .
diag(Sw)Tψ is the entry-wise product of Tψ and Sw. 0
Eq. (6) can be used to find ψ or w from z and the other
Apossibleglobalphasefactorbetweenthesolution
variable. •
and the reconstruction is taken into account in the
The Fourier transform relationship used in equations
error calculation.
(1), (3) and (4) is valid under far-field and paraxial
approximation, which is the focus of the current re- AtypicalreconstructionwithSHARPusesthefollow-
lease of SHARP. For experimental geometries such as ing sequence:
Near Field[53], Fresnel [54], Fourier ptychography [55],
through-focus[56]partiallycoherentmultiplexedgeome-
tries [48, 57, 58], under-sampling conditions [59] and to
1. Input data I(q), translations x. Optional inputs: ini-
account for noise variance [41] , one can substitute the tialimageψ(0),illuminationw(0),illuminationFourier
simple Fourier transform with the appropriate propaga- maskmw andilluminationFourieramplitudesaw.
tor and variance stabilization [38].
4
The initial values for the input data and translations
2. Ifw(0)isnotprovided,initializeilluminationbysetting caneitherbeloadedfromfileorsetbyapythoninterface.
w(0)totheinverseFouriertransformofthesquareroot The starting “zero-th” initial image is loaded from file,
oftheaverageframe.
set to a random image, or taken as a constant image.
3. If ψ(0) is not provided, initialize the image by filling
ψ(0)withrandomnumbersuniformlydrawnfrom[0,1).
C. Computational Methodology
4. Build up Q, Q∗, and (Q∗Q)−1, and frames z(0) =
Qψ(0);
SHARP was developed to achieve the highest perfor-
5. Update the frames z according to [60] using projector
mance, taking advantage of the algorithm described ear-
operatorsdefinedin(Eqs. (7,8))below:
lierandusingadistributedcomputationalbackend. The
z(l):=[2βPQPa+(1−2β)Pa+β(PQ−I)]z(l−1), ptychographic reconstruction algorithm requires one to
compute the product of several linear operators (Q, Q∗,
whereβ∈(0.5,1]isascalarfactorsetbytheuser(set F, F∗, S, S∗) on a set of frames z, an image ψ and an
to0.75bydefault,whichworksinmostcases).
illumination w several times. We use a distributed GPU
6. Updateimageψ((cid:96)) usingusingEq. 9. architecture across multiple nodes for this task (Fig. 2).
Toimplementfastoperators,asetofGPUkernelsand
7. Ifdesired,computeanewilluminationwusingEq. 10.
MPI communication are necessary. The split (Qψ) and
IfmwisgivenapplytheilluminationFouriermaskcon- overlap(Q∗z)kernelsareamongthemostbandwidthde-
straint:
manding kernels and play an important role in the pro-
w((cid:96)):=F−1{(Fw)mw}, cess.
The strategy used to implement those kernels impacts
elseifwI isgivenapplytheilluminationFourierinten- directlytheoverallperformanceofthereconstructional-
sitiesconstraint:
gorithm. To divide the problem among multiple nodes,
w((cid:96)):=F−1(cid:26) Fw aw(cid:27), SHARP initially determines the size of the final image
|Fw| based on the list of translations, frames size, and resolu-
tion. Itsubsequentlyassignsalistoftranslationstoevery
elsesimplykeeptheunconstrainedilluminationw((cid:96)):=
node and loads the corresponding frames onto GPUs.
w.
Thesplit(Qψ)andFFT(F)operationsareeasilypar-
Nowcomputecenterofmassofw((cid:96)) andshiftittofix
allelizedbecauseoftheframewiseintrinsicindependence.
thetranslationoftheobject.
Summing the frames onto an image (Q∗z) requires a re-
8. If desired do background retrieval, that is, estimate duction for every image pixel across neighboring MPI
staticbackgroundandremoveitintheiterationasde-
nodes. WithineachGPUtheimageisdividedintoblocks
scribedin[41](p.7,Eq. 30).
and we first determine which frames contribute to each
9. Iterate from 5 until one of the metrics from Eqs. block. Thecontributingframesaresummedandthenthe
11,12,14 drops below a user defined level or untill a
resultingimageissummedacrossallMPInodes. Weuse
maximum iteration for time-critical applications, and
returnψ((cid:96)) andw. shared memory or constant memory, depending on GPU
compute capability, to store frame translations, and we
usekernelfusiontoreduceaccesstoglobalmemory. The
last step of summing across all MPI nodes does not nec-
The metrics ε , ε ,ε ,ε used to monitor progress
∆ a Q ∆ essarily have to be done at every iteration, at the cost of
amd stagnation are the normalized mean square root er-
slower convergence [61], but that is the default.
ror (nmse) from the corresponding projections of z:
Timing to compute the overlap at each iteration de-
pends on the size of the image and number of frames on
ε (z)= (cid:107)[Pa−I]z(cid:107), (11)
a (cid:107)a(cid:107) top of each pixel, i.e. the density but not the size of the
ε (z)= (cid:107)[PQ−I]z(cid:107), (12) frames.
Q (cid:107)a(cid:107)
In addition to the high performance ptychographic al-
ε∆(cid:16)z(l),z(l−1)(cid:17)= (cid:107)z(l)−(cid:107)az(cid:107)(l−1)(cid:107) (13) gorithm, the SHARP software environment provides a
flexible and modular framework which can be changed
where I is the identity operator, and z(0) =0, and adapted to different needs. Furthermore, the user
For benchmarking purposes, when using a simulation hascontrolofseveraloptionsforthereconstructionalgo-
from a known solution ψ , the following metric can also rithm,whichcanbeusedtoguaranteeabalancebetween
0
be used: performanceandqualityoftheresults. Theseincludethe
choiceofilluminationFouriermask, illuminationFourier
ε0(z)= (cid:107)Q∗1z0(cid:107)mϕin(cid:13)(cid:13)Q∗(eiϕz−z0)(cid:13)(cid:13), (14) intensities and the β parameter, as well as how often
to do different operations such as illumination retrieval,
where ϕ is an arbitrary global phase factor, and z = background retrieval, and synchronization of the differ-
0
Qψ . Notice the additional scaling factor Q∗ used in ε . ent GPUs. For more details we refer the reader to the
0 0
5
squarerooterror10−2 εεεεaQ0∆
mean10−3
Normalized10−4
0 50 100 150 200
Iteration
20
me[s]15
ti
Reconstruction105
FIG. 2. Schematic of the ptychographic reconstruction algo-
0
rithm implemented in SHARP. The iterative reconstruction 0 2 4 6 8 10 12 14 16
Nr.ofGPUs
schemeisshownontheright. Toachievethehighestpossible
throughput and scalability one has to parallelize across mul-
tiple GPUs as shown on the left for the case of 4 GPUs. As FIG. 3. Convergence rate (top) per iteration and timing
mostptychographicscansuseaconstantdensityofscanpoint (bottom)toprocess10,000framesofdimension128×128ex-
acrosstheobject,weexpecttobeabletoachieveaveryeven tractedfromanimageofsize1000×1000asafunctionofthe
division, resulting in good load balancing. SHARP enforces number of nodes. All residuals decrease rapidly; numerical
anoverlapconstraintbetweentheimagesproducedbyeachof precision limits the (weighted) comparison with the known
the GPUs, and also enforces that the illumination recovered solution ε(cid:48)(z). Reconstruction is achieved (ε <5e−4) in
0
on each GPU agree with each other. This is done by default under2secondsusingaclusterwith4computenodeswith4
at every iteration. GTXTitanGPUpernode(16total,43000cores),96GBytes
GPUmemory, 1TByte RAM,and 24TBytes storage, infini-
band. Timingcontributionsforcorrespondingcomputational
kernelsare(Q∗Q)−1Q∗ 30%,F,F∗ 20%,Q20%,S=ast5%,
documentation (http://www.camera.lbl.gov/sharp). elementwise operations 20 %, and residual calculation 5 %.
No illumination retrieval was done, as the exact illumination
wasgiven. Thesimulationwasdoneusingperiodicboundary
conditions to avoid edge effects.
III. APPLICATIONS AND PERFORMANCE
A. Simulations and performance
SHARP enables high-throughput streaming analy-
sis using computationally efficient phase retrieval algo- As a demonstration, we start from a sample that was
rithms. In this section we describe a typical dataset and composed of colloidal gold nanoparticles of 50 nm and
sample that can be collected in less than 1 minute at 10 nm deposited on a transparent silicon nitride mem-
the ALS, and the computational backend to provide fast brane. Anexperimentalimagewasobtainedbyscanning
feedback to the microscopist. electron microscopy, which provides high resolution and
To characterize our performance, we use both simu- contrast but can only view the surface of the sample.
lations and experimental data. We use simulations to We simulate a complex transmission function by scal-
comparetheconvergenceofthereconstructionalgorithm ing the image amplitude from 0 to 50 nm thickness, and
to the “true solution” and characterize the effect of dif- using the complex index of refraction of gold at 750 eV
ferent light sources, contrast, scale, noise, detectors or energy from [henke.lbl.gov]. The illumination is gener-
samples for which no data exists yet. ated by a zone-plate with a diameter of 220 microns and
60 nm outer zone width, discretized into (128 128) pix-
ExperimentaldatafromALSusedtocharacterizebat-
×
els in the far field.
tery materials, green cement, magnetic materials, at dif-
ferent wavelengths and orientation has been successfully
reconstructed[62–66]usingthesoftwaredescribedinthis
article. B. Experimental example
Wealsodescribeastreamingexampleinwhichafront-
end that operates very close to the actual experiment Figure 4 shows ptychographic reconstructions of a
sends the data to the reconstruction backend that runs dataset generated from a sample consisting of gold balls
remotely on a GPU/CPU cluster. Further details about withdiametersof50and10nm. Thedataweregenerated
thestreamingfront-endandprocessingback-endpipeline using 750 eV x-rays at beamline 5.3.2.1 of the Advanced
will be published in an upcoming paper by our group. Light Source, with high stability position control of the
6
Trigger CXI file
Processing / Reconstruction (SHARP)
Data stream Backend
Frontend
CCD
Scan x/y
Framegrabber
FIG. 4. Reconstruction of a test sample consisting of gold Experiment Graphical User Interface
control (GUI)
balls with diameters of 50 and 10 nm. Detector pixel size
30 microns, 1920×960 pixels 80 mm downstream from the
sample, cropped and downsampled to 128, scan of (50× 50)
points,illuminationisgeneratedbyazone-platewithadiam- FIG.5. Overviewofthecomponentsinvolvedinthesoftware
eter of 220 microns and 60 nm outer zone width. A) Phase structureofthestreamingpipeline. Inordertomaximizethe
imagegeneratedbySHARPusingthealgorithmdescribedin performance of this streaming framework, the frontend oper-
sectionIIBapplyingtheilluminationFouriermaskconstraint ates very close to the actual experiment, while the backend
and turning on background retrieval. The red arrow points runs remotely on a powerful GPU/CPU cluster. As soon as
toacollectionof50nmballswhilethebluearrowpointstoa diffractiondataisrecordedbytheCCDcamera,aliveviewof
collectionof10nmballs. Thepixelsizeis10nm. B)Sameas theptychographicreconstructionistransmittedtotheGraph-
(A) except without enforcing the illumination Fourier mask. icalUserInterface,andtheuserisabletocontrolandmonitor
C) Same as (A) but without using the background retrieval (top panel) the current status of the data streams and anal-
algorithm. ysis, (bottom right panel) .
soft x-ray scanning transmission microscope. Exposure This software architecture allows users an intuitive,
time was 1 second and the dataset consists of a square flexible and responsive monitoring and control of their
scan grid with 40 nm spacing ( see [62] for details of the experiments. Suchatightintegrationbetweendataaqui-
experimental setup). The reconstructions consisted of sition and analysis is required to give users the feedback
300 iterations of the RAAR algorithm with a illumina- they expect from a STXM instrument.
tion retrieval and background retrieval step every other
iteration. The initial illumination is generated by (1)
computing the average of the measurements, (2) seting IV. CONCLUSIONS
everything below a threshold to 0, and everything above
a threshold to a constant average value (3) applying an In this paper we described SHARP, a high-
inverseFFT.Theimageisinitializedwithcomplexinde- performance software environment for ptychography re-
pendent identically distributed (i.i.d.) pixels, or a con- constructions, and its application as part of quick feed-
stant average value. back system used by the ptychographic mircoscopes in-
stalled at the Advanced Light Source.
Our software provides a modular interface to the high
C. Interface and Streaming performancecomputationalback-endandcanbeadapted
to different needs. Its fast throughput provides near real
Common processing pipelines used for ptychographic time feedback to microscopists and this also makes it
experiments usually have a series of I/O operations and suitable as a corner stone for demanding higher dimen-
many different components involved. We have devel- sional analysis such as spectro-ptychography or tomo-
oped a streaming pipeline, to be deployed at the COS- ptychography.
MICbeamlineattheALS,whichallowsuserstomonitor With the coming new generation light sources and
andquicklyactuponchangesalongtheexperimentaland faster detectors, the ability to quickly analyse vast
computational pipeline. amountsofdatatoobtainlargehigh-dimensionalimages
Thestreamingpipelineiscomposedofafront-endand will be an enabling tool for science.
aback-end(Fig.5). Thefront-endconsistsofaGraphical
User Interface (see Fig. 5), a worker that grabs frames
fromthedetector,andaninterfacethatmonitorsnetwork V. ACKNOWLEDGMENTS
activity, experimental parameters (position, wavelength,
exposure, etc...), and provides a live view of the ongoing We acknowledge useful discussions with Chao Yang,
reconstruction. H-T Wu, J. Qian and Z. Wen. This work was par-
On the back-end side, the streaming infrastructure is tially funded by the Center for Applied Mathematics
composedofacommunicationhandlerandacollectionof for Energy Research Applications, a joint ASCR-BES
workers addressing different tasks such as dark calibra- funded project within the Office of Science, US Depart-
tion, detector correction, data reduction, ptychographic mentofEnergy, undercontractnumberDOE-DE-AC03-
reconstruction and writing output to file. 76SF00098, by theSwedishResearchCouncil andby the
7
Swedish Foundation for Strategic Research.
[1] J. R. Fienup, Appl. Opt. 21, 2758 (1982). [math.OC].
[2] D. R. Luke, J. V. Burke, and R. G. Lyon, SIAM review [30] C.Kewish,P.Thibault,M.Dierolf,O.Bunk,A.Menzel,
44, 169 (2002). J. Vila-Comamala, K. Jefimovs, and F. Pfeiffer, Ultra-
[3] J. R. Fienup, J. C. Marron, T. J. Schulz, and J. H. microscopy 110, 325 (2010).
Seldin, Appl. Opt. 32, 1747 (1993). [31] S. Ho¨nig, R. Hoppe, J. Patommel, A. Schropp,
[4] M. Eckert, Acta Crystallographica Section A 68, 30 S. Stephan, S. Scho¨der, M. Burghammer, and C. G.
(2012). Schroer, Opt. Express 19, 16324 (2011).
[5] A. Momose, T. Takeda, Y. Itai, and K. Hirano, Nature [32] M. Guizar-Sicairos, S. Narayanan, A. Stein, M. Metzler,
medicine 2, 473 (1996). A. R. Sandy, J. R. Fienup, and K. Evans-Lutterodt,
[6] R. Collier, C. Burckhardt, and L. Lin, Applied Physics Letters 98, 111108 (2011).
Optical holography (Academic Press, New York, [33] A. Schropp, R. Hoppe, V. Meier, J. Patommel, F. Sei-
1971). both, H. J. Lee, B. Nagler, E. C. Galtier, B. Arnold,
[7] S. Marchesini, S. Boutet, and e. a. Sakdinawat, Nature U. Zastrau, et al., Scientific reports 3, 01633 (2013).
Photonics 2, 560 (2008). [34] A. Wojdyla, R. Miyakawa, and P. Naulleau, in
[8] P. W. Hawkes and J. C. H. Spence, eds., SPIE Advanced Lithography (International Society for
Science of microscopy (Springer, New York, 2007). Optics and Photonics, 2014) pp. 904839–904839.
[9] A.FannjiangandW.Liao,J.Opt.Soc.Am.A29,1847 [35] M. Guizar-Sicairos and J. R. Fienup, Opt. Express 16,
(2012). 7264 (2008).
[10] Y. Wang and Z. Xu, ArXiv e-prints 1310.0873 (2013). [36] S. Marchesini, Y.-C. Tu, and H.-t. Wu, Ap-
[11] F.Pfeiffer,T.Weitkamp,O.Bunk, andC.David,Nature plied and Computational Harmonic Analysis ,
physics 2, 258 (2006). 10.1016/j.acha.2015.06.005 (2015).
[12] E.J.Candes,Y.C.Eldar,T.Strohmer, andV.Voronin- [37] Z. Wen, C. Yang, X. Liu, and S. Marchesini, Inverse
ski, SIAM Journal on Imaging Sciences 6, 199 (2013). Problems 28, 115010 (2012).
[13] I. Waldspurger, A. d?Aspremont, and S. Mallat, Math- [38] C. Yang, J. Qian, A. Schirotzek,
ematical Programming (2013). F. Maia, and S. Marchesini,
[14] J.Miao,P.Charalambous,J.Kirz, andD.Sayre,Nature Iterative Algorithms for Ptychographic Phase Retrieval,
400, 342 (1999). Tech. Rep. 4598E, arXiv:1105.5628 (Lawrence Berkeley
[15] S.Marchesini,H.He,H.Chapman,S.Hau-Riege,A.Noy, National Laboratory, 2011).
M. Howells, U. Weierstall, and J. Spence, Phys. Rev. B [39] P. Thibault and M. Guizar-Sicairos, New Journal of
68, 140101 (2003). Physics 14, 063004 (2012).
[16] H.H.Bauschke,P.L.Combettes, andD.R.Luke,Jour- [40] J.Qian,C.Yang,A.Schirotzek,F.Maia, andS.March-
nal of the Optical Society of America. A, Optics, image esini, Inverse Problems and Applications, Contemp.
science, and vision 19, 1334 (2002). Math 615, 261 (2014).
[17] S. Marchesini, Rev Sci Instrum 78, 011301 (2007), [41] S. Marchesini, A. Schirotzek, C. Yang, H.-t. Wu, and
arXiv:physics/0603201. F. Maia, Inverse Problems 29, 115009 (2013).
[18] S. Marchesini, JOSA A 24, 3289 (2007). [42] A. Maiden, M. Humphry, M. Sarahan, B. Kraus, and
[19] J. M. Rodenburg, Advances in Imaging and Electron J. Rodenburg, Ultramicroscopy 120, 64 (2012).
Physics 150, 87 (2008). [43] M. Beckers, T. Senkbeil, T. Gorniak, K. Giewekemeyer,
[20] W. Hoppe, Acta Crystallographica Section A 25, 495 T.Salditt, andA.Rosenhahn,Ultramicroscopy 126,44
(1969). (2013).
[21] R. Hegerl and W. Hoppe, Berichte Der Bunsen- [44] P. Godard, M. Allain, V. Chamard, and J. Rodenburg,
Gesellschaft Fur Physikalische Chemie 74, 1148 (1970). Opt. Express 20, 25914 (2012).
[22] P. D. Nellist, B. C. McCallum, and J. M. Rodenburg, [45] B. Abbey, K. A. Nugent, G. J. Williams, J. N. Clark,
Nature 374, 630 (1995). A.G.Peele,M.A.Pfeifer,M.deJonge, andI.McNulty,
[23] H. N. Chapman, Ultramicroscopy 66, 153 (1996). Nature Physics 4, 394 (2008).
[24] Y.S.G.Nashed,D.J.Vine,T.Peterka,J.Deng,R.Ross, [46] N. C. Jesse and G. P. Andrew, Applied Physics Letters
and C. Jacobsen, Opt. Express 22, 32082 (2014). 99, 154103 (2011).
[25] F. R. Maia, Nature methods 9, 854 (2012). [47] L. W. Whitehead, G. J. Williams, H. M. Quiney, D. J.
[26] P. Thibault, M. Dierolf, A. Menzel, O. Bunk, C. David, Vine, R. A. Dilanian, S. Flewett, K. A. Nugent, A. G.
and F. Pfeiffer, Science 321, 379 (2008). Peele,E.Balaur, andI.McNulty,Phys.Rev.Lett.103,
[27] P. Thibault, M. Dierolf, O. Bunk, A. Menzel, and 243902 (2009).
F. Pfeiffer, Ultramicroscopy 109, 338 (2009). [48] L.Tian,X.Li,K.Ramchandran, andL.Waller,Biomed-
[28] R.Hesse,D.R.Luke,S.Sabach, andM.K.Tam,SIAM ical Optics Express 5, 2376 (2014).
Journal on Imaging Sciences 8, 426 (2015). [49] S. T. Thurman and J. R. Fienup, JOSA A 26, 1008
[29] S. Marchesini and H.-T. Wu, (2009).
Rank-1 accelerated illumination recovery in scanning diffract[i5v0e]imMa.gGinugizbayr-tSriacnasirpoasreanncdyJe.stRim.aFtiieonnu,p, Opt. Express 17,
Tech. Rep. LBNL-6734E, arXiv:1408.1922 (Lawrence 2670 (2009).
Berkeley National Laboratory, 2014) arXiv:1408.1922 [51] F. J. Anscombe, Biometrika 35, 246 (1948).
8
[52] M. Ma¨kitalo and A. Foi, IEEE Transactions on Image [61] J. Liu, S. J. Wright, C. R´e, V. Bittorf, and S. Srid-
Processing 22, 91 (2013). har, The Journal of Machine Learning Research 16, 285
[53] M. Stockmar, P. Cloetens, I. Zanette, B. Enders, (2015).
M. Dierolf, F. Pfeiffer, and P. Thibault, Scientific re- [62] D.A.Shapiro,Y.-S.Yu,T.Tyliszczak,J.Cabana,R.Ce-
ports 3, 01927 (2013). lestre, W. Chao, K. Kaznatcheev, A. Kilcoyne, F. Maia,
[54] D. Vine, G. Williams, B. Abbey, M. Pfeifer, J. Clark, S. Marchesini, et al., Nature Photonics (2014).
M. De Jonge, I. McNulty, A. Peele, and K. Nugent, [63] Y.-S.Yu, C.Kim, D.A.Shapiro, M.Farmand, D.Qian,
Physical Review A 80, 063823 (2009). T.Tyliszczak,A.D.Kilcoyne,R.Celestre,S.Marchesini,
[55] G.Zheng,R.Horstmeyer, andC.Yang,NaturePhoton- J. Joseph, et al., Nano letters 15, 4282 (2015).
ics 7, 739 (2013). [64] S.Bae,R.Taylor,D.Shapiro,P.Denes,J.Joseph,R.Ce-
[56] J. Marrison, L. Ra¨ty, P. Marriott, and P. O’Toole, Sci- lestre,S.Marchesini,H.Padmore,T.Tyliszczak,T.War-
entific reports 3, 02369 (2013). wick,etal.,JournaloftheAmericanCeramicSociety98,
[57] D. J. Batey, D. Claus, and J. M. Rodenburg, Ultrami- 4090 (2015).
croscopy 138, 13 (2014). [65] Y.Li,S.Meyer,J.Lim,S.C.Lee,W.E.Gent,S.March-
[58] S.Dong,R.Shiradkar,P.Nanda, andG.Zheng,Biomed- esini, H. Krishnan, T. Tyliszczak, D. Shapiro, A. L. D.
ical Optics Express 5, 1757 (2014). Kilcoyne, et al., Advanced Materials 27, 6590 (2015).
[59] T.B.Edo, D.J.Batey,A.M.Maiden,C.Rau,U.Wag- [66] X.Shi,P.Fischer,V.Neu,D.Elefant,J.Lee,D.Shapiro,
ner, Z. D. Peˇsi´c, T. A. Waigh, and J. M. Rodenburg, M. Farmand, T. Tyliszczak, H.-W. Shiu, S. Marchesini,
Phys. Rev. A 87, 053850 (2013). et al., Applied Physics Letters 108, 094103 (2016).
[60] R. Luke, Inverse Problems 21, 37 (2005).
