PTYCHOGRAPHY IMAGING SYSTEMS AND METHODS WITH CONVEX RELAXATION
Certain aspects pertain to ptychographic imaging systems and methods with convex relaxation. In some aspects, a ptychographic imaging system with convex relaxation comprises one or more electromagnetic radiation sources, a digital radiation intensity detector, and a processor in communication with the digital radiation detector. The electromagnetic radiation provides coherent radiation to a specimen while the digital radiation intensity detector receives light transferred from the sample by diffractive optics and captures intensity distributions for a sequence of low resolution images having diversity. The processor generates a convex problem based on the sequence of low resolution images and optimizes the convex problem to reconstruct a highresolution image of the specimen. In certain aspects, the convex problem is relaxed into a lowrank formulation.
This application is a continuation of U.S. patent application Ser. No. 14/710,947, titled “Ptychography Imaging Systems and Methods with Convex Relaxation” and filed on May 13, 2015, which claims priority to U.S. Provisional Patent Application No. 61/992,505, titled “Solving Conventional and Fourier Ptychography with Convex Optimization” and filed on May 13, 2014; each of these applications is hereby incorporated by reference in its entirety and for all purposes.
BACKGROUNDCertain embodiments described herein are generally related to digital imaging, and more specifically, to standard and Fourier ptychographic imaging systems and methods with convex relaxation.
Ptychography imaging involves collecting lower resolution images and then reconstructing the image data to form a higher resolution image. Over the past two decades, ptychographic imaging has been used in a variety of regimes to produce highresolution, wide fieldofview images of microscopic and nanoscopic phenomena. Whether in the Xray regime at thirdgeneration synchrotron sources, in the electron microscope for atomic scale phenomena, or in the in the optical regime for biological specimens, ptychography has shown an unparalleled ability to acquire hundreds of megapixels of sample information near the diffraction limit. Typically, the underlying operation of ptychography is to sample a series of diffraction patterns from a specimen as it is scanned through a focused beam. These intensityonly measurements are then reconstructed into a complex (i.e. amplitude and phase), highresolution image with more pixels of sample information than any single recorded diffraction pattern.
Most recently, a Fourier ptychographic microscope (FPM) was introduced that uses a Fourier ptychographic technique that can reconstruct gigapixel optical images from a sequence of lower resolution images collected using a low NA objective lens from a conventional microscope. In one example, Fourier ptychographic microscope activates different LEDs in an LED array to illuminate a sample from different directions while the lowresolution images are captured. As in standard ptychography, Fourier ptychography recovers the sample's phase as it merges together the captured image sequence into a highresolution output.
Conventional ptychographic imaging systems can avoid the need for a high NA, wellcorrected objective lens to image at the diffractionlimit by resolving resolutionlimiting factors in their data capture and reconstruction techniques. However, these systems lack stable, robust, and accurate reconstruction methods. For example, conventional ptychographic systems reconstruct the phase of the scattered field from measured intensities using nonconvex algorithms. Most of these conventional systems solve the phase retrieval problem by applying known constraints in an iterative manner using an “alternating projection” (AP) strategy. Reconstruction techniques that use AP strategies tend to converge to incorrect local minima and/or to stagnate.
SUMMARYCertain aspects of this disclosure pertain to standard and Fourier ptychographic imaging systems and methods with convex relaxation.
Certain aspects pertain to ptychographic imaging systems with convex relaxation. In some aspects, a ptychographic imaging system with convex relaxation comprises one or more electromagnetic radiation sources, a digital radiation intensity detector, and a processor in communication with the digital radiation detector. In some cases, the one or more electromagnetic radiation sources configured to provide coherent radiation to a specimen from a plurality of incidence angles at a sequence of sample times. For example, the one or more electromagnetic radiation sources may be an LED array. The digital radiation intensity detector configured to receive light transferred from the sample by diffractive optics. In these cases, the digital radiation intensity detector further configured to capture intensity distributions for a sequence of low resolution intensity images associated with the plurality of incidence angles. The processor in communication with the digital radiation detector to receive image data with intensity distributions for the sequence of low resolution intensity images. The processor configured to generate a convex problem based on the sequence of low resolution intensity images and optimize the convex problem to reconstruct a highresolution image of the specimen. In other cases, the diffractive optics comprises a spatial light modulator configured to provide a pattern at a plurality of locations at a Fourier plane of the specimen. In these cases, the sequence of low resolution images capture by the digital radiation detector is associated with the plurality of locations of the pattern.
Certain aspects pertain to ptychographic imaging methods with convex relaxation. In some aspects, a ptychographic imaging method with convex relaxation comprises collecting a sequence of low resolution images of a specimen, generating a convex problem based on the sequence of low resolution intensity images, and optimizing the convex problem to reconstruct a highresolution image of the specimen. In some cases, collecting a sequence of low resolution images of a specimen comprises providing coherent radiation to a specimen from a plurality of incidence angles at a sequence of sample times, transferring light from the specimen through diffractive optics to a digital radiation intensity detector, and sampling a sequence of low resolution intensity images associated with the plurality of incidence angles. In some cases, generating the convex problem comprises stacking image data from the lowresolution images into a combined image matrix, constructing measurement matrices, and generating the convex problem using convex relaxation.
These and other features are described in more detail below with reference to the associated drawings.
Embodiments of the present invention will be described below with reference to the accompanying drawings. The features illustrated in the drawings may not be to scale.
Certain aspects are directed to ptychographic (standard or Fourier) imaging systems and methods with convex relaxation. Typically, these ptychographic (standard or Fourier) imaging systems comprise coherent electromagnetic (EM) radiation sources for illuminating the specimen, diffraction optics, and a digital radiation intensity detector for taking intensity measurements. The imaging method with convex relaxation starts by providing coherent EM radiation to illuminate the specimen being imaged. The diffractive optics receives light from the specimen and transfers diffraction patterns to the digital radiation intensity detector, which samples intensity measurements at a sequence of sample times. In the case of a standard ptychography imaging system, the diffraction optics are freespace propagation. In the case of a Fourier ptychography imaging system, the diffraction optics are a combination of freespace propagation and a collection optical element. The digital intensity measurements captured at each sample time provide image data of a lowresolution image. The method reconstructs a high resolution image from the image data of the sequence of lowresolution images. This method manipulates the image matrix to form a convex problem by stacking the reduced resolution images in combined matrix, constructing measurement matrices, and using convex relaxation to create the convex program. The method then reconstructs a high resolution image from the convex image matrix by either (Option A) relaxing the convex image matrix into lowrank formulation and then solving the lowrank image matrix using a lowrank ptychography (LRP) process or (Option B) determining the high resolution image directly with a convex lifted ptychography (CLP) technique.
Certain aspects are directed to ptychographic (standard or Fourier) imaging systems and methods with convex relaxation. In some aspects, ptychographic imaging systems and methods with convex relaxation may provide the desired stability, robustness and/or reliability in reconstructing high resolution images from a collection of lower resolution images. In some aspects, these ptychographic imaging systems and methods use reconstruction techniques with convex relaxation that do not have local minima, incorporate noise compensation techniques, and/or use multiple a priori constraints. In one aspect, a reconstruction technique with convex relaxation uses lowrank factorization, whose runtime and memory usage are nearlinear with respect to the size of the output image. In this example, the reconstruction technique may be able to provide 25% lower background variance than conventional ptychographic reconstruction methods.
In certain aspects, ptychographic imaging systems and methods with convex relaxation include an unaided recovery technique that does not use prior sample knowledge or an appropriate heuristic, which may be especially relevant in biological imaging. Moreover, these imaging systems and methods do not have local minima so that a single solution can be found efficiently. In addition, these imaging systems and methods are more noisetolerant than imaging systems that use AP strategies, which makes results more reproducible. Furthermore, a factorization technique can be implemented to obtain solutions at scale. Thus, these aspects provide noisetolerant and efficient reconstruction techniques that may provide for more efficiently and accurately processing multigigapixel high resolution images than conventional systems.
I. Ptychographic Imaging Systems with Convex Relaxation
As mentioned above, the image collection assembly 12 comprises a diversity element 25 that may be a feature(s) of either the one or more coherent EM radiation sources 20 or the diffractive optics 30. The diversity element 25 refers to one or more feature(s) that implement a change between sample times to cause diversity in the captured image data. In some cases, the diversity element 25 may be provided by configuring the one or more coherent EM radiation sources 20 with multiple light sources (i.e. point emitters) providing illumination from different incidence angles to the specimen during sampling. For example, the coherent EM radiation sources 20 may be in the form of a twodimensional LED array (n×m dimensions) of LEDs acting as point emitting light sources at different locations at the illumination plane of the twodimensional LED array. In another example, the diversity element 25 may be provided by mechanically shifting the specimen to different locations at the sample plane using, for example, an XY stage. In another example, the diversity element 25 may be provided by configuring the diffractive optics 30 to comprise a spatial light modulator with its display located at the Fourier plane. The diversity can then be generated during sampling by displaying a pattern at different locations on the spatial light modulator display. In yet another example, the diversity element 25 may be provided by shifting another coded mask around the Fourier plane.
In some cases, the one or more coherent EM radiation sources 20 and diffractive optics 30 are configured to operate in a transillumination mode directing illumination through the specimen and toward a collection element of the diffractive optics 30. In other cases, the one or more coherent EM radiation sources 20 and diffractive optics 30 are configured in epiillumination mode directing illumination toward the specimen and away from a collection element of the diffractive optics 30.
In certain aspects, a digital radiation intensity detector comprises a twodimensional grid of equally spaced discrete elements (e.g., pixels) at a detection plane. At each sample time, each element samples intensity of radiation received. As a group, the grid samples a twodimensional intensity distribution associated with the location of the elements. The digital radiation intensity detector generates a signal(s) with frames of image data of the intensity distribution measured by the grid of radiation detecting elements at the detection plane at each sample time. If visible light radiation is being used to illuminate the specimen, the digital radiation intensity detector may be in the form of a charge coupled device (CCD), a CMOS imaging sensor, an avalanche photodiode (APD) array, a photodiode (PD) array, a photomultiplier tube (PMT) array, or like device. If using THz radiation is used, the digital radiation intensity detector may be, for example, an imaging bolometer. If using microwave radiation, the digital radiation intensity detector may be, for example, an antenna. If Xray radiation is used, the digital radiation intensity detector may be, for example, an xray sensitive CCD. If acoustic radiation is used, the digital radiation intensity detector may be, for example, a piezoelectric transducer array. These examples of digital radiation intensity detectors and others are commercially available. In some aspects, the digital radiation intensity detector may be a color detector e.g., an RGB detector. In other aspects, the digital radiation intensity detector may be a monochromatic detector.
As shown in
The one or more processors 52 may receive instructions stored on the CRM 54 (e.g., memory) and execute those instructions to perform one or more functions of system 10. For example, the processor(s) 52 may execute instructions to perform one or more steps of the imaging with convex relaxation method. For example, the processor(s) 52 may execute instructions stored on the CRM 54 to perform one or more functions of the system 10 such as, for example, 1) interpreting image data, 2) reconstructing a higher resolution image from the image data, and 3) providing display data for displaying one or more images or other output on the display 56. As another example, the processor(s) 52 may provide control instructions for controlling the illumination to the coherent EM radiation source(s) 20. In one case, the processor(s) 52 may provide control instructions to synchronize the illumination by coherent EM radiation source(s) 20 with the sampling times of the digital radiation intensity detector 40. In addition to storing instructions for preforming certain functions of the system 10, the CRM 54 can also store the (lower resolution) intensity and higher resolution image data, and other data produced by the system 10. The display 56 may be a color display or a black and white display. In addition, the display 56 may be a twodimensional display or a threedimensional display. In one embodiment, the display 56 may be capable of displaying multiple views.
More specifically,
To collect lowresolution images using the illustrated Fourier ptychographic imaging system with convex relaxation 300, the low NA lens 330 receives illumination altered by the specimen and filters the illumination based on its low NA. The digital radiation intensity detector 340 receives the filtered illumination from the low NA lens 330 and samples a lowresolution image set b comprising a sequence of j lowresolution images at j sample times during illumination by the j illumination sources. Each lowresolution image is captured at a different sample time during illumination from a different incidence angle. As illustrated by the bottom illustration, the system 300 uses an imaging method with convex relaxation that uses a convex lifted ptychographic (CLP) technique to transform the image set b into a highresolution complex sample image ψ.
To collect lowresolution images using this Fourier ptychographic imaging system with convex relaxation 400, the low NA lens 230 receives illumination altered by the specimen and filters the illumination based on its low NA. The digital radiation intensity detector 240 receives the filtered illumination from the low NA lens 230 and samples a sequence of j lowresolution images at j sample times during illumination by the j illumination sources. Each lowresolution image is captured at a different sample time during illumination from a different incidence angle.
The locations of neighboring apertures have an overlapping area between neighboring apertures such as, for example, the overlapping area 437 between aperture 436(1) and aperture 436(2). When using the reconstruction method with convex relaxation, the overlapping area need only be about 50% or lower of the area of one of the neighboring apertures to converge to a single imaging solution. Conventional FPM reconstruction required a more extensive overlapping area in the range of 80 to 90% in order to converge to an accurate solution. Since more overlap is required, more images and iterations are needed to cover the same area in the conventional systems. Thus, conventional systems required more exposure time and more resources to reconstruct the high resolution image.
II. Ptychographic Imaging Methods with Convex Relaxation
At this juncture in the method, there are two possible options (Option A) go to step 730 or (Option B) go to step 750. Generally speaking, if the convex combined image matrix contains largescale ptychographic data (i.e. the low resolution images have a high number of pixels and/or the sequence has a high number of low resolution images), then Option A may be the more appropriate option. For example, Option A may be used if the number of pixels is more than 50×50 pixels for each low resolution image and/or the number of images is more than 200. If the convex combined image matrix has smaller scale data, Option B may be the appropriate option.
If Option A is used, the method relaxes the convex combined image matrix into lowrank formulation at step 730. Once in lowrank formulation, the method a minima of the new lowrank formulation to determine a highresolution image. The method determines the minima based on a lowrank ptychographic (LRP) technique at step 740.
If Option B is used, the method determines a minima of X to determine a highresolution image at step 750. In this case, the minima is determined using a convex lifted ptychographic (CLP) technique. An example of an appropriate convex solver can be found, for example, at http://cvxr.com/tfocs/.
In the subsections that follow, certain steps described with reference to the flowchart in
A. Collection of Low (Reduced) Resolution Images
This section describes example substeps of the collection of low resolution images step 710 of
With respect to the Fourier ptychographic imaging systems in
Based on this point emitter assumption, the j^{th }light source (e.g., LED illuminated in
s(x,j)=ψ(x)e^{ikxp}^{j} (Eqn. 1)
Where: wavenumber k=2π/λ and

 p_{j}=sin θ_{j }describes the offaxis angle of the j^{th }flight source
At each illumination by a flight source, a j^{th }illuminated sample field s(x,j) is collected by the collection optics having a low NA such as the low NA lens 230 ofFIG. 2 .
 p_{j}=sin θ_{j }describes the offaxis angle of the j^{th }flight source
Neglecting scaling factors and a quadratic phase factor for simplicity, a Fourier optics setup gives the field at the imaging system pupil plane, A(x′), as [s(x,j)]={circumflex over (ψ)}(x′−p_{j}). Here, represents the Fourier transform between conjugate variables x and x′, where {circumflex over (ψ)} is the Fourier transform of ψ, and the Fourier shift property has been applied. The shifted sample spectrum field {circumflex over (ψ)} (x′−p_{j}) is then modulated by the imaging system's aperture function a(x′), which acts as a lowpass filter. In
The spectrum {circumflex over (ψ)} is considered discretized into n pixels with a maximum spatial frequency k. The bandpass cutoff of the aperture function a is denoted as k·m/n, where m is an integer less than n. The modulation of {circumflex over (ψ)} by a results in a field characterized by m discrete samples, which propagates to the imaging plane at g(x) and is sampled by an mpixel digital radiation intensity detector 240. The mpixel digital radiation intensity detector 240 samples intensity distribution related to reducedresolution images that can be combined into a reduced resolution (low resolution) image matrix, g, as:
g(x,j)=[a(x′){circumflex over (ψ)}(x′−p_{j})]^{2} (Eqn. 2)

 Where: g(x,j) is an (m×q) dimensional Fourier ptychographic data matrix.
That is, the j^{th }column contains a lowresolution image of sample intensity while under illumination from the j^{th }optical source.
 Where: g(x,j) is an (m×q) dimensional Fourier ptychographic data matrix.
In the imaging method with convex relaxation, a higher resolution (npixel) complex spectrum {circumflex over (ψ)}(x′) is reconstructed from the plurality of lowresolution (mpixel) intensity measurements contained within the data matrix g. Once {circumflex over (ψ)} is found, an inverseFourier transform yields the desired complex sample reconstruction, ψ.
Standard ptychographic systems resolve the inverse problem using alternating projections (AP) strategies: after initializing a complex sample estimate, ψ_{0}, iterative constraints help force ψ_{0 }to obey all known physical conditions. First, its amplitude is forced to obey the measured intensity measurement set from the detector plane (i.e., the values in g). Second, its spectrum {circumflex over (ψ)}_{0 }is forced to lie within a known support in the plane that is Fourier conjugate to the detector. While these AP strategies are known to converge when each constraint set is convex, the intensity constraint applied at the detector plane is not convex, which may sometimes lead to erroneous solutions and/or stagnation in finding a solution.
In one example, the components of the Fourier ptychographic system shown in
B. Generate Convex Combined Image Matrix
In this subsection, an embodiment of the step 720 of
To solve Eqn. 2 as a convex problem, it is expressed in matrix form. First, the unknown sample spectrum {circumflex over (ψ)} is represented as an (n×1) vector where n is the known sample resolution before being reduced by the finite bandpass of the lens aperture. Next, the j^{th }detected lower resolution image becomes an (m×1) vector g_{j}, where m is the number of pixels in each lowresolution image. The ratio n/m defines the ptychographic resolution improvement factor. It is equivalent to the largest angle of incidence from an offaxis optical source, divided by the acceptance angle of the imaging lens. Third, each lens aperture function a(x+p_{j}) is expressed as an (n×1) discrete aperture vector a_{j}, which modulates the unknown sample spectrum {circumflex over (ψ)}.
To rewrite Eqn. 2 as a matrix product {Aj}_{j=1}^{q }defined as a sequence of (m×n) rectangular matrices that contain a deterministic aperture function a_{j }along a diagonal. For an aberrationfree rectangular aperture, each matrix Aj has a diagonal of ones originating at (0, p′_{j}) and terminating at (m, p′_{j}+m−1), where p′_{j }is now a discretized version of the shift variable, p_{j}. Finally, m×m discrete Fourier transform matrix F^{(m) }is introduced to express the transformation of the lowpass filtered sample spectrum through the fixed imaging system for each lowresolution image g_{j}:
g_{j}F(m)A_{j}{circumflex over (ψ)}^{2}, 1<j<q (Eqn. 3)
The ptychographic system collects a sequence of q lower resolution images, {g_{j}}_{j=1}^{q }such as, for example, at step 710 of
b=FAV{circumflex over (ψ)}^{2}=D{circumflex over (ψ)}^{2} (Eqn. 4)
This is an example of stacking in substep 910 of
The transpose of the i^{th }row of the D matrix is denoted as d_{i}, which is a column vector. The set {d_{i}} forms the measurement vectors. The measured intensity in the i^{th }pixel is the square of the inner product between d_{i }and the spectrum {circumflex over (ψ)}: b_{i}=d_{i},{circumflex over (ψ)}^{2}.
The method “lifts” the solution {circumflex over (ψ)} out of the quadratic relationship in Eqn. 4 and expresses the solution {circumflex over (ψ)} in the space of (n×n) positivesemidefinite matrices:
b_{i}=Tr({circumflex over (ψ)}*d_{i}d_{i}*{circumflex over (ψ)})=Tr(d_{i}d_{i}*{circumflex over (ψ)}{circumflex over (ψ)}*)=Tr(D_{i}X) (Eqn. 6)
where D_{i}=d_{i}d_{i}* is a rank1 measurement matrix constructed from the i^{th }measurement vector d_{i}, X={circumflex over (ψ)}{circumflex over (ψ)}* is an (n×n) rank1 outer product, and 1≤i≤q·m. This is an example of constructing measurement matrices at step 920 of
Eqn. 6 states that quadratic image measurements {b_{i}}_{i=1}^{q·m }are linear transforms of in a higher dimensional space. These q·m linear transforms are combined into a single linear operator to summarize the relationship between the stacked image vector b and the matrix X as, (X)=b.
The phase retrieval problem in ptychography can be posed as the following rank minimization process:
minimize rank(X)
subject to (X)=b, (Eqn. 7)

 X>0, where X>0 denotes that X is positivesemidefinite
However, the minimization problem in Eqn. 7 is not convex. In order to transform the minimization program in Eqn. 7 into a convex problem, convex relaxation is performed on Eqn. 7 by replacing the rank of matrix X with its trace, which generates a convex semidefinite program as follows:
minimize Tr(X)
subject to (X)=b, (Eqn. 8)

 X>0,
To account for the presence of noise, Eqn. 8 may be reformed such that the measured intensities in b are no longer strictly enforced constraints, but instead appear in the objective function as follows:
minimize αTr(X)+½∥A(X)−b∥
subject to X>0 (Eqn. 9)
Here, α is a scalar regularization variable that directly trades off goodness for complexity of fit. Its optimal value depends upon the assumed noise level. Eqn. 9 forms the final convex combined image matrix problem that can be used to recover a resolution improved complex sample ψ from a set of obliquely illuminated images in b.
Both Eqn. 8 and Eqn. 9 use convex relaxation to generate a convex program with a function to minimize (minimization function) and solution constraints. In addition to generating a convex program, Eqn. 9 also accounts for noise. In certain aspects, the minimization function has no local minima.
If Option B is used in
C. Results from Using Convex Lifted Ptychographic (CLP) Technique
In this Section, some simulated results of using Option B of
In certain aspects, using a CLP process returns a lowrank matrix X, with a rapidly decaying spectrum, as the optimal solution of Eqn. 9. The trace term in the CLP objective function is primarily responsible for enforcing the lowrank structure of X. While this trace term also appears like an alternative method to minimize the unknown signal energy, a fair interpretation should consider its effect in a lifted (n×n) solution space. The final complex image estimate ψ can be obtained by first performing a singular value decomposition of X. Given lownoise imaging conditions and spatially coherent illumination, ψ is set to the Fourier transform of the largest resulting singular vector. Viewed as an autocorrelation matrix, useful statistical measurements may also be found within the remaining smaller singular vectors of X. One may also identify X as the discrete mutual intensity matrix of a partially coherent optical field: X={circumflex over (ψ)}{circumflex over (ψ)}*, where denotes an ensemble average. Under this interpretation, Eqn. 9 becomes an alternative solver for the stationary mixed states of a ptychography setup.
Three points distinguish Eqn. 9 from conventional APbased ptychography strategies. First, the convex CLP process has a larger search space. If AP strategies are used to iteratively update an npixel estimate, Eqn. 9 must solve for an n×n positivesemidefinite matrix. Second, this boost in the solution space dimension guarantees that the convex program may find a global optimum with tractable computation. This allows CLP technique to avoid AP's frequent convergence to local minima (i.e., failure to approach the true image). Unlike conventional solvers for the ptychography problem, no local minima exist in the CLP process of embodiments. Finally, Eqn. 9 considers the presence of noise by offering a parameter (a) to tune with an assumed noise level. APbased solvers lack this parameter and can be easily led into incorrect local minima by even low noise levels as discussed in the following section.
CLP Technique Simulations and Noise Performance
In the example results described in this section, the Fourier ptychographic systems of
In the noiseless case, five (5) iterations of (NonConvex) AP reconstruction introduced unpredictable artifacts to both the recovered amplitude and phase images as shown in
The AP and CLP reconstructions were repeated again setting α=0.001 in Eqn. 9 while varying two relevant parameters: the number of captured images q, and their signaltonoise ratio (SNR). In these reconstructions, SNR=10 log_{10}(ψ^{2}/N^{2}), where ψ^{2}θ is the mean sample intensity and N^{2} is the mean intensity of uniform Gaussian noise added to each simulated raw image. To account for the unknown constant phase offset in all phase retrieval reconstructions, the reconstruction meansquared error defined as MSE=Σ_{x}ψ(x)−ρs(x)^{2}/Σ_{x}ψ(x)^{2}, where ρ=Σ_{x}ψ(x)s*(x)/Σ_{x}s(x)^{2 }is a constant phase factor shifting the reconstructed phase to optimally match the known phase of the ground truth sample.
When using the reconstruction method with convex relaxation of certain embodiments, each image spectrum need only overlap with each neighboring image spectrum by about 50% or lower to converge to a single imaging solution. Conventional FPM reconstruction methods require more extensive overlapping, e.g., 8090%, to converge to an accurate solution. Since more overlap is required in conventional systems, more images and iterations are needed to cover the same Fourier space. Thus, conventional systems require more exposure time and more resources to reconstruct the high resolution image.
D. Factorization for LowRank Ptychography (LRP)
Within this subsection, certain details of steps 730 and 740 of
Although generating a convex program based on Eqn. 9 makes reconstruction more efficient, the constraint that X remain positivesemidefinite can be computationally burdensome since each iteration may require a full eigenvalue decomposition of X. In order to process largescale ptychographic data, certain aspects segment each detected image into tiles (e.g., 10^{3 }pixel tile images), and process each tile (segment) separately, and then “tile” the resulting reconstructions back together into a final full resolution image. Tiling parallelization may increase efficiency for processing largescale ptychographic data.
In certain embodiments, the method of convex relaxation provides another process for processing largescale ptychographic data that may be used as an alternative or in conjunction with tiling parallelization. In these embodiments, the method may use Option A of
In these cases, instead of solving for an n×n matrix X, the method can use a lowrank ansatz and factorize the matrix X as X=RR^{T}, where the decision variable R is now an n×r rectangular matrix containing complex entries, with r<n. Inserting this factorization into the optimization problem in Eqn. 8 and writing the constraints in terms of the measurement matrix D_{i}=d_{i}d_{i}^{T }generates the nonconvex program,
Minimize Tr(RR^{T})
subject to Tr(D_{i}RR^{T})=b_{i }for all i. (Eqn. 10)
Besides removing the positive semidefinite constraint in Eqn. 8, the factored form of Eqn. 10 presents two more key adjustments to the original convex formulation. First, using the relationship Tr(RR^{T})=∥R∥_{F}^{2}, where F denotes a Frobenius norm, it is direct to rewrite the objective function and each constraint in Eqn. 10 with just one n×r decision matrix, R. Now instead of storing an n×n matrix like CLP, LRP must only store an n×r matrix. Since most practical applications of ptychography require coherent optics, the desired solution rank r will typically be close to 1, thus significantly relaxing storage requirements (i.e., coherent light satisfies X={circumflex over (ψ)}{circumflex over (ψ)}*, so one would expect R as a column vector and RRT a rank1 outer product). Fixing r at a small value, LRP memory usage now scales linearly instead of quadratically with the number of reconstructed pixels, n. Second, the feasible set of Eqn. 10 is no longer convex and thus, the solution strategy must be shifted away from a simple semidefinite program to use an LRP technique instead.
In one aspect, the method can use the nonconvex program of minimizing Tr(RR^{T})+R_1 is used with the same constraints of Eqn. 10. This latter term is the L1 norm of the unknown sample spectrum matrix, R.
At step 740 of
In certain aspects, to solve Eqn. 10, the following augmented Lagrangian function is minimized:
Where R∈C^{n×r }is the unknown decision variable and the two variables y∈R^{q·m }and σ∈R^{+} are parameters to guide the method to its final reconstruction. Each b_{i }is the intensity measured by one image sensor pixel. The LRP technique iteratively minimizes the function L by sequentially updating R and the parameters Y and σ.
The first term in Eqn. 11 is the objective function from Eqn. 10, indirectly encouraging a lowrank factorized product. This tracks the original assumption of a rank1 solution within a “lifted” solution space. The second term contains the known equality constraints in Eqn. 10 (i.e., the measured intensities), each assigned a weight y_{i}. The third term is a penalized fitting error that is abbreviated with label v. It is weighted by one penalty parameter σ, mimicking the role of a Lagrangian multiplier.
With an appropriate fixed selection of y_{i}'s and σ, the minimization of L(R, y, σ) with respect to R identifies the desired optimum of Eqn. 10. Specifically, if a local minimum of L is identified each iteration (which is nearly always the case in practice), then the minimization sequence accumulation point may be a guaranteed solution. As an unconstrained function, the minimum of L can be found quickly using a quasiNewton approach.
In some aspects, the goal of a lowrank ptychography (LRP) process is to determine a suitable set of (y_{i}, σ) to minimize Eqn. 11 with respect to R, which leads to the desired high resolution image solution. An iterative process is used to sequentially minimize L with respect to R^{k }at iteration k, and then update a new parameter set (y^{k+1}, σ^{k+1}) at iteration k+1. The parameters (y^{k+1}, σ^{k+1}) are updated to ensure their associated term's contribution to the summation forming L is relatively small. This suggests R^{k+1 }is proceeding to a more feasible solution. The relative permissible size of the second and third terms in L are controlled by two important parameters, η<1 and γ>1: if the third term v sufficiently decreases such that v^{k+1}≤ηv^{k}, then multiplier σ is held fixed and the equality constraint multipliers, y_{i}, is updated. Otherwise, σ is increased by a factor γ such that σ^{k+1}=γσ^{k}.
In one example, the LRP process is initialized with an estimate of the unknown highresolution complex sample function ψ_{0}, contained within a lowrank matrix R^{0}. The LRP process terminates if it reaches a sufficient number of iterations or if the minimizer fulfills some convergence criterion. R^{0 }is formed using a spectral method, which can help increase solver accuracy and decrease computation time. Specifically, the r columns of R^{0 }are selected as the leading r eigenvectors of D*diag[b]D, where D is the measurement matrix in Eqn. 4. While this spectral approach works quite well in practice, a random initialization of R^{0 }also often produces an accurate reconstruction.
In certain aspects, the LRP process uses two parameters, γ and η which help guide its solution process. In most of the included experiments, the LRP process set γ=1.5. The most observable consequence of selecting a different value of γ is its influence on the number of iterations needed for a desired level of performance. A larger value for γ will cause a larger change within the augmented Lagrangian function each iteration, and thus a quicker progression to an optimized reconstruction.
In certain aspects, the LRP process sets η=0.5. A different value of η may also significantly alter the quality of reconstruction output with experimental input images (i.e., in the presence of noise).
LRP Simulations and Noise Performance
Following the same procedure used above to simulate the CLP process, MSE performance of the LRP process was tested as a function of SNR and the results are shown in
These simulations arbitrarily fix η and γ at 0.5 and 1.5, respectively, and set the desired rank of the solution, r, to 1. In some embodiments, these free variables may be changed which offers significant freedom to tune the response of LRP process to noise. For example, similar to the noise parameter a in Eqn. 9, the multiplier a (controlled via γ) in Eqn. 11 helps trade off complexity for goodness of fit by reweighting the quadratic fitting error term.
In addition to reducing required memory, the LRP process also improves upon the computational cost of CLP process for large scale data. For an npixel sample reconstruction, the per iteration cost of the CLP algorithm is currently O(n^{3}), using bigO notation. The positivesemidefinite constraint in Eqn. 9, which requires a full eigenvalue decomposition, defines this behavior limit. The periteration cost of the LRP algorithm, on the other hand, is O(n log n). This large periteration cost reduction is the primary source of LRP speedup.
Experimental Results
Experiments were used to verify that the LRP process of certain embodiments can improves the accuracy and noise stability of ptychographic reconstruction by using a Fourier ptychographic (FP) microscope with the configuration of the system in
Quantitative Performance
In these examples, the Fourier ptychographic microscope was comprised of a 15×15 array of surfacemounted LEDs (e.g., model SMD 3528, center wavelength λ=632 nm, 4 mm LED pitch, 150 μm active area diameter), which served as quasicoherent optical sources. The LED array was placed l=80 mm beneath the sample plane, and each LED has an approximate 20 nm spectral bandwidth.
To quantitatively verify resolution improvement, each of the 15×15 LEDs was turned on beneath a U.S. Air Force (USAF) resolution calibration target. A microscope objective (e.g. 2λ e.g., Olympus® objective with apochromatic Plan APO 0.08 NA) transferred each resulting optical field to a CCD detector (e.g., Kodak® KAI29050 detector with 5.5 μm sized pixels), which sampled 225 low resolution images. Using this 0.08 NA microscope objective (5° collection angle) and a 0.35 illumination NA (θ_{max}=20° illumination angle), the Fourier ptychographic microscope provided a total complex field resolution gain of n/m=25. Each image spectrum overlapped by ol≈70% in area with each neighboring image spectrum. For reconstruction, n=25·m was selected and the same aperture parameters were used with both AP and LRP processes to create the highresolution images shown in
Both ˜1 megapixel reconstructions achieved their maximum expected resolving power (i.e., resolved Group 9, Element 3: 1.56 μm line pair spacing). This is approximately 5 times sharper than the smallest resolved feature in one raw image (e.g., Group 7, Element 2). The LRP process avoids certain artifacts that are commonly observed during the nonlinear descent of AP process. Both reconstructions slowly fluctuate in background areas that are expected to be uniformly bright or dark. These fluctuations are caused in part by experimental noise, an imperfect aperture function estimate, and possible misalignments in the LED shift values, pj. In a representative background area marked by a 402 pixel blue box in
The LRP process was then demonstrated on a “highNA” FP microscope configuration comprising a larger 0.5 NA microscope objective lens with a 30° collection angle (e.g., 20× Olympus 0.5 NA UPLFLN). For sample illumination, 28 LEDs are arranged into 3 concentric rings of 8, 8 and 12 evenly spaced light sources (ring radii=16, 32 and 40 mm, respectively). This new light source array was placed 40 mm beneath the sample to create a 0.7 illumination NA with a θ_{max}=45° illumination angle. The synthesized numerical aperture of this FP microscope, computed as the sum of the illumination NA and objective lens NA, is NAs=1.2. With a greaterthanunity synthetic NA, reconstructions can offer oilimmersion quality resolution (˜385 nm smallest resolvable feature spacing), without requiring any immersion medium between the sample and objective lens. A monolayer of polystyrene microspheres (index of refraction n_{m}=1.587) immersed in oil (n_{o}=1.515, both indexes for λ=632 nm light) was imaged by this highNA″ FP microscope configuration.
By using a LRP process of certain embodiments, the process was shown to reconstruct quantitatively accurate phase. Using the same data and parameters for AP and LRP input, the highresolution phase reconstructions were obtained of two adjacent microspheres in
Biological Sample Reconstruction
The third imaging example uses the same highNA FP configuration (i.e. collection angle=30°, θ_{max}=45°) of the embodiments discussed with respect to
The infected blood sample was prepared to maintain erythrocyte asexual stage cultures of the P. falciparum strain 3D7 in culture medium, then smeared, fixed with methanol, and then a Hema 3 stain was applied. An example sample region containing two infected cells, imaged with a conventional highNA oilimmersion microscope (NA=1.25) under Kohler illumination. Twenty eight (28) uniquely illuminated images were captured of these two infected cells using the highNA FP microscope. The top two right images contain an example normally illuminated raw image, which does not clearly resolve the parasite infection. Bottom left six images presents phase retrieval reconstructions using the standard AP algorithm, where m=1202, n=2402, run 6 iterations, and again subtract a constant phase offset. Reconstructions from three data sets were included: images captured with a 1 second exposure (top), a 0.25 second exposure (middle), and 0.1 second exposure (bottom). A shorter exposure time implies increased noise within each raw image. While the 1 sec exposurebased AP reconstruction resolves each parasite, blurred cell boundaries and nonuniform fluctuations in amplitude suggest an inaccurate AP convergence. However, both parasite infections remain visible within the reconstructed phase. The parasites become challenging to resolve within the phase from 0.25 sec exposure data, and are not resolved within the phase from the 0.1 sec exposure data, due to increased image noise. The normalized background variance of each AP amplitude reconstruction, from a representative 402pixel window (marked blue square), is σ2=0.0020, 0.0027, and 0.0059 for the 1 sec, 0.25 sec, and 0.1 sec exposure reconstructions, respectively.
For comparison, reconstructions using the LRP process are shown (sharpest solutions after 15 iterations). For each reconstructed amplitude, the desired solution rank is set to r=3. The 3 modes of the resulting reconstruction are added in an intensity basis to create the displayed amplitude images. For each reconstructed phase, the desired solution matrix rank is set to r=1 and all other parameters are left unchanged. For all three exposure levels, the amplitude of the cell boundaries remains sharper than in the AP images. Both parasite infections are resolvable in either the reconstructed amplitude or phase, or both, for all three exposure levels. The normalized amplitude variances from the same background window are now σ2=0.0016 (1 sec),0.0022 (0.25 sec), and 0.0035 (0.1 sec), an average reduction (i.e., improvement) of 26% with respect to the AP results. The AP reconstructions here offer a generally flatter background phase profile than LRP (i.e., less variation at low spatial frequencies). Without additional filtering or postprocessing, the AP algorithm here might offer superior quantitative analysis during e.g. tomographic cell reconstruction, where loworder phase variations must remain accurate. However, it is clear within
Through the relaxation in Eqn. 8, the traditionally nonlinear phase retrieval process for ptychography is transformed into a convex program. The convex program can be solved with a CLP process if it is a smallscale image set. If it is a largescale image set, the convex program can be relaxed into low rank formulation resulting semidefinite program with an appropriate factorization, and then solved with a LRP process. This method of convex relaxation provides a process that is robust to noise.
Besides removing local minima from the recovery process, perhaps the most significant departure from conventional phase retrieval solvers is a tunable solution rank, r. As noted earlier, r connects to statistical features of the ptychographic experiment, typically arising from the partial coherence of the illuminating field. Coherence effects are significant at thirdgeneration Xray synchrotron sources and within electron microscopes. An appropriately selected r may eventually help LRP process measure the partial coherence of such sources. The solution rank may also help identify setup vibrations, sample autofluorescence, or even 3D sample structure. In some cases, the method of convex relaxation can artificially enlarge the solution rank to encourage the transfer of experimental noise into its smaller singular vectors. Other extensions of LRP include simultaneously solving for unknown aberrations (i.e., the shape of the probe in standard ptychography), systematic setup errors, and/or inserting additional sample priors such as sparsity.
III. SubsystemsThe various components previously described in the Figures may operate using one or more of the subsystems to facilitate the functions described herein. Any of the components in the Figures may use any suitable number of subsystems to facilitate the functions described herein. Examples of such subsystems and/or components are shown in a
In some embodiments, an output device such as the printer 2430 or display 56 of the Fourier camera system can output various forms of data. For example, the Fourier camera system can output 2D color/monochromatic images (intensity and/or phase), data associated with these images, or other data associated with analyses performed by the Fourier camera system.
Modifications, additions, or omissions may be made to any of the abovedescribed embodiments without departing from the scope of the disclosure. Any of the embodiments described above may include more, fewer, or other features without departing from the scope of the disclosure. Additionally, the steps of the described features may be performed in any suitable order without departing from the scope of the disclosure.
It should be understood that the present invention as described above can be implemented in the form of control logic using computer software in a modular or integrated manner. Based on the disclosure and teachings provided herein, a person of ordinary skill in the art will know and appreciate other ways and/or methods to implement the present invention using hardware and a combination of hardware and software.
Any of the software components or functions described in this application, may be implemented as software code to be executed by a processor using any suitable computer language such as, for example, Java, C++ or Perl using, for example, conventional or objectoriented techniques. The software code may be stored as a series of instructions, or commands on a CRM, such as a random access memory (RAM), a read only memory (ROM), a magnetic medium such as a harddrive or a floppy disk, or an optical medium such as a CDROM. Any such CRM may reside on or within a single computational apparatus, and may be present on or within different computational apparatuses within a system or network.
Although the foregoing disclosed embodiments have been described in some detail to facilitate understanding, the described embodiments are to be considered illustrative and not limiting. It will be apparent to one of ordinary skill in the art that certain changes and modifications can be practiced within the scope of the appended claims.
One or more features from any embodiment may be combined with one or more features of any other embodiment without departing from the scope of the disclosure. Further, modifications, additions, or omissions may be made to any embodiment without departing from the scope of the disclosure. The components of any embodiment may be integrated or separated according to particular needs without departing from the scope of the disclosure.
Claims
1. A ptychographic imaging system with convex relaxation, the system comprising:
 one or more electromagnetic radiation sources configured to provide coherent radiation to a specimen;
 a digital radiation intensity detector configured to receive light transferred from the sample by diffractive optics and to capture intensity distributions for a sequence of low resolution images of the specimen; and
 a processor in communication with the digital radiation detector to receive the intensity distributions the sequence of low resolution images, the processor configured to: generate a convex problem based on the sequence of low resolution images; and optimize the convex problem to reconstruct a highresolution image of the specimen.
2. The ptychographic imaging system with convex relaxation of claim 1,
 wherein the one or more electromagnetic radiation sources are configured to provide coherent radiation to a specimen from a plurality of incidence angles at a sequence of sample times; and
 wherein the sequence of low resolution images is associated with the plurality of incidence angles.
3. The ptychographic imaging system with convex relaxation of claim 1,
 wherein the diffractive optics comprises a spatial light modulator configured to provide a pattern at a plurality of locations at a Fourier plane of the specimen; and
 wherein the sequence of low resolution images is associated with the plurality of locations.
4. The ptychographic imaging system with convex relaxation of claim 1, wherein the processor generates the convex problem by:
 stacking image data from the lowresolution images into a combined image matrix;
 constructing measurement matrices; and
 generating the convex problem using convex relaxation.
5. The ptychographic imaging system with convex relaxation of claim 1, where the convex problem comprises a minimization function and solution restraints.
6. The ptychographic imaging system with convex relaxation of claim 5, wherein the minimization function has no local minima.
7. The ptychographic imaging system with convex relaxation of claim 1, wherein the diffractive optics comprises a low NA objective lens.
8. The ptychographic imaging system with convex relaxation of claim 1, wherein the processor optimizes the convex problem by:
 relaxing the convex problem into a lowrank formulation; and
 reconstructing the highresolution imaging of the specimen using a lowrank ptychography process.
9. The ptychographic imaging system with convex relaxation of claim 1, wherein if the convex problem is large scale, then the processor optimizes the convex problem by:
 relaxing the convex problem into a lowrank formulation; and
 reconstructing the highresolution imaging of the specimen using a lowrank ptychography process.
10. The ptychographic imaging system with convex relaxation of claim 1, wherein the convex problem is large scale if the number of low resolution images is above 200 or each low resolution image has a resolution of more than 50×50 pixels.
11. The ptychographic imaging system with convex relaxation of claim 1, wherein the one or more electromagnetic radiation sources are in the form of an LED array.
12. A ptychographic imaging method with convex relaxation, the method comprising:
 collecting a sequence of low resolution images of a specimen;
 generating a convex problem based on the sequence of low resolution images; and
 optimizing the convex problem to reconstruct a highresolution image of the specimen.
13. A ptychographic imaging method with convex relaxation of claim 10, wherein collecting a sequence of low resolution images of a specimen comprises:
 providing coherent radiation to the specimen from a plurality of incidence angles at a sequence of sample times;
 transferring light from the specimen through diffractive optics to a digital radiation intensity detector; and
 sampling a sequence of low resolution images associated with the plurality of incidence angles.
14. A ptychographic imaging method with convex relaxation of claim 10, wherein collecting a sequence of low resolution images of a specimen comprises:
 illuminating the specimen;
 providing a pattern at a plurality of locations on a display of a spatial light modulator located at the Fourier plane of the specimen;
 transferring light from the specimen using diffractive optics with the spatial light modulator to a digital radiation intensity detector; and
 sampling a sequence of low resolution images associated with the plurality of locations of the pattern.
15. The ptychographic imaging method with convex relaxation of claim 10, wherein generating the convex problem comprises:
 stacking image data from the lowresolution images into a combined image matrix;
 constructing measurement matrices; and
 generating the convex problem using convex relaxation.
16. The ptychographic imaging method with convex relaxation of claim 10, where the convex problem comprises a minimization function and solution restraints, wherein the minimization function has no local minima.
17. The ptychographic imaging method with convex relaxation of claim 10, wherein optimizing the convex problem comprises:
 relaxing the convex problem into a lowrank formulation; and
 reconstructing the highresolution imaging of the specimen using a lowrank ptychography process.
18. The ptychographic imaging method with convex relaxation of claim 10, further comprising:
 if the convex problem is large scale, optimizing the convex problem by relaxing the convex problem into a lowrank formulation and reconstructing the highresolution imaging of the specimen using a lowrank ptychography process.
19. The ptychographic imaging method with convex relaxation of claim 15, wherein the convex problem is large scale if the number of low resolution images is above 200 or each low resolution image has a resolution of more than 50×50 pixels.
Type: Application
Filed: Oct 25, 2018
Publication Date: Feb 21, 2019
Inventors: Roarke W. Horstmeyer (Palo Alto, CA), Yuhua Chen (Pasadena, CA), Joel A. Tropp (Los Angeles, CA), Changhuei Yang (South Pasadena, CA)
Application Number: 16/171,270