Computational Microscopy: Turning Megapixels into Gigapixels Citation Horstmeyer, Roarke William (2016) Computational Microscopy: Turning Megapixels into Gigapixels. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/Z95Q4T1W. https://resolver.caltech.edu/CaltechTHESIS:10202015-173005082 Abstract The layout of a typical optical microscope has remained effectively unchanged over the past century. Besides the widespread adoption of digital focal plane arrays, relatively few innovations have helped improve standard imaging with bright-field microscopes. This thesis presents a new microscope imaging method, termed Fourier ptychography, which uses an LED to provide variable sample illumination and post-processing algorithms to recover useful sample information. Examples include increasing the resolution of megapixel-scale images to one gigapixel, measuring quantitative phase, achieving oil-immersion quality resolution without an immersion medium, and recovering complex three dimensional sample structure. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: Computational imaging, biophotonics, microscopy, optics Degree Grantor: California Institute of Technology Division: Engineering and Applied Science Major Option: Electrical Engineering Awards: Everhart Distinguished Graduate Student Lecturer Award, 2015. Charles and Ellen Wilts Prize, 2016. Thesis Availability: Public (worldwide access) Research Advisor(s): Yang, Changhuei Thesis Committee: Yang, Changhuei (chair) Tropp, Joel A. Vaidyanathan, P. P. Hassibi, Babak Zheng, Guoan Defense Date: 15 October 2015 Non-Caltech Author Email: roarke.horstmeyer (AT) gmail.com Funders: Funding Agency Grant Number National Institutes of Health 1DP2OD007307-01 National Institutes of Health R01AI096226-01 Record Number: CaltechTHESIS:10202015-173005082 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:10202015-173005082 DOI: 10.7907/Z95Q4T1W ORCID: Author ORCID Horstmeyer, Roarke William 0000-0002-2480-9141 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 9231 Collection: CaltechTHESIS Deposited By: Roarke Horstmeyer Deposited On: 12 Nov 2015 22:11 Last Modified: 08 Nov 2023 18:41 Thesis Files Preview PDF - Final Version See Usage Policy. 38MB Repository Staff Only: item control page

Computational microscopy: turning megapixels into gigapixels Thesis by Roarke Horstmeyer In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy California Institute of Technology Pasadena, California 2016 (Defended October 15, 2015) ii c 2015 Roarke Horstmeyer All Rights Reserved iii Acknowledgments This work would not be possible without the help of the Biophotonics group at Caltech. First, I would like to thank Guoan Zheng for initiating this field of research, Xiaoze Ou for all his experimental knowhow, and Jinho Kim, Jaebum Chung and Joshua Brake for their additional help with making Fourier ptychography into a practical tool. My knowledge of optics would not be what it is today without many lessons from other group members, past and present: Benjamin Judkewitz, Ivo M. Vellekoop, Ying Min Wang, Mooseok Jang, Hao Wen, Edward Zhou, Chao Han, Sean Peng, Seung Ah Lee, and Jian Ren, among many others who I am likely leaving out. Finally, a big thanks to Changhuei Yang for all of his support and guidance throughout my time at Caltech. I was lucky enough to receive a lot of help and guidance from outside of the Biophotonics group. Those from Caltech include Richard Chen, Joel Tropp, Sid Assawaworrarit, Mark Harfouche, Dongwan Kim, Christos Santis, Alon Greenbaum, Fabio Arai, P. P. Vaidyanathan, Babak Hassibi, Christos Thramboulidis, Ramya Vinayak, Wael Halbawi, Kishore Jaganathan and all of the others in the EE department who answered many of my pestering questions. I also received a lot of guidance from outside of Caltech. Thank you to Ravi Athale, Gary Euliss, Michael Stenner, Chuck DiMarzio, Joseph Hollmann, George Barbastathis, Laura Waller, Lei Tian, Salman Asif, Jason Holloway, Ashok Veeraraghavan, Ollie Cossairt, Ioannis Papadopoulos, Mike Henninger, Ed Boyden, Ramesh Raskar, Otkrist Gupta, Nikhil Naik, Doug Lanman, Ankit Mohan, Tom Cuypers, and the rest of the early members of the Camera Culture group, as well as all of the others at the Media Lab. And finally, a special thank you to Se Baek Oh, who helped me through my early years of studying optics at MIT, and planted in my head the seed of the idea to use phase retrieval as a means to improve microscope image resolution. As for financial support, I would first like to thank the National Defense Science and Engineering Graduate Fellowship (NDSEG) for providing me with the freedom to hop around projects, until I was able to identify some really great experiments. Thanks are also due to the National Institutes of Health (grant numbers 1DP2OD007307-01, R01AI096226-01), Clearbridge Biophotonics Pte Ltd., Singapore (Agency Award no. Clearbridge 1), the Office of Naval Research (award N00014-11-1002), a Sloan Research Fellowship, and the Moore Foundation. Finally, I’d like to thank my family and friends for, you know, everything. Sorry for not listing all of you out, but you know who you are. And of course Margie, ’cause she’s the best! iv Abstract The layout of a typical optical microscope has remained effectively unchanged over the past century. Besides the widespread adoption of digital focal plane arrays, relatively few innovations have helped improve standard imaging with bright-field microscopes. This thesis presents a new microscope imaging method, termed Fourier ptychography, which uses an LED to provide variable sample illumination and post-processing algorithms to recover useful sample information. Examples include increasing the resolution of megapixel-scale images to one gigapixel, measuring quantitative phase, achieving oil-immersion quality resolution without an immersion medium, and recovering complex three dimensional sample structure. v Contents Acknowledgments iii Abstract iv 1 Introduction 1 1.1 The simple microscope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.1 The illumination plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.2 The sample plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.3 The aperture plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.4 The image plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.2 Relevant definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3 Challenges for the standard microscope . . . . . . . . . . . . . . . . . . . . . . . . . 10 Bibliography 14 2 Fourier ptychography for gigapixel imaging 17 2.1 Resolution improvement through shifting illumination . . . . . . . . . . . . . . . . . 17 2.2 Mathematical model of Fourier ptychography . . . . . . . . . . . . . . . . . . . . . . 19 2.3 Reconstruction goals and extensions . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 Bibliography 24 3 Fourier ptychographic image reconstruction 26 3.1 The phase retrieval algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.2 The Fourier ptychographic microscopy (FPM) algorithm . . . . . . . . . . . . . . . . 28 3.3 Experimental demonstration of FPM . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 Bibliography 36 4 Fourier ptychography for quantitative phase measurement 37 4.1 Quantitative phase for sample thickness . . . . . . . . . . . . . . . . . . . . . . . . . 37 vi 4.1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.1.2 Verification with microspheres . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.1.3 Verification with human blood cells . . . . . . . . . . . . . . . . . . . . . . . . 41 4.2 FPM phase measurements for digital pathology . . . . . . . . . . . . . . . . . . . . . 42 4.3 Improving reconstruction quality with FPM phase measurements . . . . . . . . . . . 46 4.3.1 Digital refocusing for enhanced depth-of-field . . . . . . . . . . . . . . . . . . 46 4.3.2 Simultaneously removal of microscope aberrations . . . . . . . . . . . . . . . 47 Bibliography 49 5 Modeling ptychography in phase space 52 5.1 The conventional ptychography (CP) setup . . . . . . . . . . . . . . . . . . . . . . . 53 5.1.1 Phase space representation of CP . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.2 Mathematical representation of Fourier ptychography (FP) . . . . . . . . . . . . . . 56 5.3 Visualizing connections between both ptychographic domains . . . . . . . . . . . . . 60 5.4 A complete statistical model with partially coherent light . . . . . . . . . . . . . . . 62 5.4.1 Partially coherent source description . . . . . . . . . . . . . . . . . . . . . . . 63 5.4.2 CP with partially coherent light . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.4.3 FP with partially coherent light . . . . . . . . . . . . . . . . . . . . . . . . . . 64 Case study: CP and FP under partially coherent illumination . . . . . . . . . . . . . 66 5.5.1 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.5.2 Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.5 Bibliography 70 6 Solving ptychography with a convex relaxation 72 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 6.2 Convex Lifted Ptychography (CLP) . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 6.2.1 Mathematical fundamentals . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 6.2.2 The CLP solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 6.2.3 CLP simulations and noise performance . . . . . . . . . . . . . . . . . . . . . 79 Factorization for Low-Rank Ptychography (LRP) . . . . . . . . . . . . . . . . . . . . 80 6.3.1 The LRP solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 6.3.2 LRP simulations and noise performance . . . . . . . . . . . . . . . . . . . . . 82 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 6.4.1 Quantitative performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 6.4.2 Biological sample reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . 88 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 6.3 6.4 6.5 vii Bibliography 92 7 Fourier ptychography in a conventional camera 95 7.1 Introduction and Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 7.2 Principle of operation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 7.2.1 Image acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 7.2.2 Aberration-free OFC reconstruction . . . . . . . . . . . . . . . . . . . . . . . 100 7.2.3 Aberration-free OFC simulation 7.3 . . . . . . . . . . . . . . . . . . . . . . . . . 101 OFC with simulated annealing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 7.3.1 Characterization and removal of low-order aberrations . . . . . . . . . . . . . 103 7.3.2 Characterization and removal of geometric distortion . . . . . . . . . . . . . . 106 7.3.3 The complete OFC algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 7.4 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 7.5 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 Bibliography 115 8 Diffraction tomography with Fourier ptychography 118 8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 8.2 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 8.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 8.4 8.5 8.3.1 Image formation in FPT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 8.3.2 FPT reconstruction algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 8.4.1 Quantitative verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 8.4.2 Biological experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 Conclusion Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 136 1 Chapter 1 Introduction The microscope is an invaluable tool for scientific discovery. The fundamental aim of magnifying small objects dates back at least several thousand years, to the invention of the simple lens [1]. The finding that two lenses, when placed in sequence, can create an image with a very large magnification likely dates to the 16th century. Early “compound” microscopes invented during this era truly opened up a new visual world to the curious eye. Astonishing views of cells and bacteria were achieved around the same time that the telescope was offering astronomers their first glimpses of our neighboring planets. Over the past century, a number of findings have pushed the microscope into new realms. Frits Zernike’s insights during the 1930’s and 1940’s led both to the development of aberration theory [2], as well as the creation of the phase contrast microscope [3, 4], for which he was awarded the Nobel Prize. The concept that light travels as a wave, with a defined amplitude and phase, was wellknown for many years prior to Zernike’s work. However, the insight that it might be possible to capture this phase information, within an intensity-only image, led to significant breakthroughs for measurements in vivo. A related technique, termed differential interference contrast, was developed around the same time by Nomarski [5]. Both methods effectively mix an optical field’s phase into the field’s amplitude, which can suddenly reveal the three-dimensional structure of cells. A decade later, Marvin Minsky’s invention of the confocal microscope [6] further transformed imaging at the micro-scale, most notably within biology. The confocal microscope aims to improve image contrast by removing the negative effects of scattered light. It applies a simple insight: all rays that do not originate from a desired focal point of interest are blocked by a pinhole. As a result, biological imaging could extend beyond imaging only the surface of organisms and materials, and now could peer into them as well. Confocal microscopes resolve sharp features from below the superficial layers of e.g., tissue, and from defined volumes within otherwise murky and opaque samples that are up to hundreds of micrometers thick. Finally, over the past several decades, three primary insights have kept the optical microscope at the forefront of biological and chemical science. First, a large collection of recently developed 2 (a) Compound microscope! Illumination! Sample! Microscope objective! Eyepiece! Intermediate image! x! z! y! (b) Simple microscope! Illumination! Detector! (c) Detailed microscope lens system! Sample! Aperture! a(x',y')! ψ(x,y)! Image! g(M x, M y)! Obj. lens! Back focal plane! Tube lens! a(x',y')! w! δx   l do! di! Figure 1.1: (a) Diagram of a compound microscope. (b) Diagram of a “simple” microscope, offering a mathematically simpler and effectively equivalent description of microscope image formation. optical methods can now both probe and activate biochemical content with visible light. Example applications include the photochemical activation of drugs [7], the photorelease of biomolecules [8], stimulation of neural activity through optogenetic tags [9], and imaging with fluorescent markers [10]. Second, projectors and advanced optical sources now allow one to pattern the illumination incident upon a microscope sample. Examples include using a digital projector or a pulsed laser to achieve a desired illumination field shape (e.g., for stimulation emission-depletion microscopy [11]). Third, and perhaps more importantly, digital detector arrays (i.e., the charged-coupled device [CCD] and complementary metal-oxide semiconductor [CMOS] pixel arrays) now directly connect the microscope to the computer. Acquired images no longer have to appear sharp and crisp, but instead can be post-processed, while taking into account knowledge of the behavior of the optical system, to extract additional information. The combination of designed illumination, fluorescent labeling and computational recovery formed the basis for the 2015 Nobel Prize in Chemistry, this year. This thesis focuses on a new microscope technique, termed Fourier ptychography (FP) [32], which uses insights from both illumination design and computational post-processing to push microscopy to new heights. Examples of what the FP microscope can achieve include transforming megapixel images into gigapixel maps, acquiring the surface profile of a sample to nanometer accuracy, resolving sub-wavelength phenomena without the need for oil immersion, and producing volumetric tomograms of thick samples. In the remainder of this section, we provide a quick review of background information that is relevant to fully understanding the new capabilities of Fourier ptychography. 3 1.1 The simple microscope A schematic of a standard compound microscope is shown in Fig. 1.1(a). Here, light originates at an illumination plane and propagates to a sample plane, where we place our object of interest for inspection. The sample is imaged via a microscope objective lens to an intermediate plane. This intermediate plane is subsequently imaged by an eyepiece lens to a digital detector. For the majority of this thesis, we will neglect the effects of the microscope eyepiece lens, since it is rarely used in digital microscope imaging. We will instead work with a simplified microscope setup, shown in Fig. 1.1(b), which we term a “simple” microscope. In practice, most microscopes in the lab use an infinity-corrected objective lens, along with a tube lens of fixed focal length (180 mm), to form an image onto the detector. For diagrammatic simplicity, we will summarize the effect of both these lenses during microscope image formation as a single focusing element, as shown in the box in Fig. 1.1(b). In practice, the infinity-corrected objective lens and tube lens form a modified 4f imaging system, whose Fourier plane is at the back focal plane of the objective lens (Fig. 1.1(c)). In our simplified diagram, we will draw the Fourier plane as a plane that is internal to our single focusing element. While this single-lens simplification neglects various minor effects, like vignetting and possible aberrations (which might be specific to the multi-lens diagram), we maintain that our simplification of the compound microscope into a single-lens schematic offers a clear and mathematically accurate description of wave-based image formation. Unless otherwise stated, this simplified picture will guide our primary mathematical models of microscope operation in this thesis. Next, we examine in detail the four primary planes of interest that compose our simple microscope. 1.1.1 The illumination plane Illumination is a critical component of any microscope. The majority of current bright-field microscopes use Kohler illumination, as diagrammed in Fig. 1.2(a). Here, we show a trans-illumination Kohler scheme, where light passes through the sample to the image plane. However, the same principle also extends to epi-illumination. In either case, Kohler designs use a set of lenses (a collector and condenser lens) to spread incoherent light as evenly as possible across the sample plane. Typically, a thermal source, such as a light bulb, creates the incoherent light. The benefits of using incoherent illumination in a microscope are two-fold. First, a wide range of optical frequencies (i.e., colors) reach the sample and subsequently transmit information about possible absorption, excitation, and fluorescence at different energy levels to the detector. Second, incoherent illumination offers a slight benefit in terms of maximum achievable image resolution, as detailed later in this thesis. An alternative illumination geometry using an incoherent source is shown in Fig. 1.2(b). Here, a mask is placed between the collector and condenser lens to create a “dark-field” image. The mask Different  types  of  microscope  illumina/on   4 (a) Kohler illumination! Bulb! Collector! Condenser! (c) Laser illumination! Sample! (b) Dark-field illumination! Patch stop! Coherent source! Sample! Wavefront! (d) LED array illumination! LEDs! Sample! Figure 1.2: Four different trans-illumination schemes that are common within microscopy. (a) Kohler illumination provides even incoherent light across the sample. (b) Dark-field illumination blocks the central rays in a Kohler setup before they reach the sample. (c) Coherent illumination is helpful, e.g., in digital holographic microscopy. (d) An LED array provides a computationally addressable, spatially coherent illumination source. prevents all rays traveling at small angles with respect to the optical axis from reaching the sample. As a result, only light traveling at large angles will hit the sample. This angled incident light will not enter the objective lens unless it is diffracted by the sample. In other words, the incident light is traveling at an angle that is too large to pass directly into the objective lens in the absence of a sample, and must somehow interact with the sample to deflect into the lens acceptance angle. Dark-field imaging is particularly helpful at highlighting fine sample features which diffract incident light into a wide range of angles. We will return to this fundamental principle in Chapter 2 during our explanation of Fourier ptychography. Instead of using an incoherent light source, it is also possible to illuminate a sample of interest with spatially and temporally coherent light, such as light from a laser. Spatially and temporally coherent light can be modeled as a wave, with a defined amplitude and phase across space. Sufficiently spatially coherent light helps to preserve any phase information contained within the optical field exiting the surface of a sample (see Fig. 1.2(c)). It is possible to measure this exiting field’s phase using digital holography or via alternative computational methods [13]. As we will detail, knowledge of this phase is also helpful when reconstructing thick, three-dimensional samples during tomographic measurements. A final type of microscope illumination that will play a major role in this thesis originates from an LED array (Figure 1.2(d)). Light from an LED lies somewhere between that originating from a thermal bulb or a laser source. Typical LED sources emit a narrow range of wavelengths (5-10 5 nm spectral bandwidth). Furthermore, the active area of an LED is quite compact, typically on the order of several hundred microns in diameter. Thus, it is often possible to consider light from one distant LED as originating from a small point source that is effectively spatially coherent and quasi-monochromatic [23] (see definition below). An array of individually addressable LEDs allows one to turn on multiple individually coherent, yet mutually incoherent sources at a time. Or, one may effectively shift one coherent source to different spatial locations along the illumination plane. Note that while certainly possible, this thesis does not consider placing any optics between the LED array and sample, or curving the LED array plane. This variable source of LED illumination is closely related to other “structured illumination” techniques used in microscopy. Specifically, a structured illumination setup creates and shines a specific pattern of optical intensity onto the sample. Examples include sinusoidal stripes [17, 18] or random speckle [19, 20]. Since structured illumination is almost always used with the goal of fluorescent imaging in mind, the phase of the illumination light is rarely manipulated. The LED array in Fig. 1.2(d) can be thought of as a method to provide structured illumination with a uniform intensity across the surface of the sample, but a spatially varying phase, whose profile depends upon the distance of the activated LED from the optical axis. 1.1.2 The sample plane The sample plane is located directly above the illumination plane. We denote its spatial coordinates as (x, y), which are perpendicular to the axis of propagation, z (see labels in Fig. 1.1). For the majority of this thesis, we will consider imaging thin samples. The thin sample condition holds if 2 /πλ, where δres is the sampling resolution and λ the maximum sample thickness t obeys t << 4δres is the illuminating light’s central wavelength [12]. If a sample satisfies this condition, then we may completely summarize its interactions with light using a complex two-dimensional function, ψ(x, y). Specifically, the amplitude of ψ(x, y) defines the amount of light absorbed by the sample at each spatial location (x, y) when illuminated with a uniform plane wave. Likewise, the phase of ψ(x, y) defines the spatially varying phase delay imparted by the sample to the incident plane wave. If the sample does not strictly satisfy the above thickness condition, we often find that a 2D function still offers a very useful sample description, up to a thickness of approximately 50 µm. Chapter 8 details how FP operates with thick samples. 1.1.3 The aperture plane Next, the optical field exiting the surface of the sample propagates into our microscope. Here, we adopt the common convention of treating the microscope as a generalized “black box” imaging system [10]. Under this interpretation, we may summarize all properties of the complicated system 6 of microscope lenses using just an entrance pupil and an exit pupil, which are each images of the same limiting aperture. This limiting aperture typically includes a physical stop within the system of lenses to minimize the effect of aberrations, and to block any stray light. We define the plane of this limiting aperture as our “aperture plane”. It may be described by an aperture function, a(x0 , y 0 ). Here, we define (x0 , y 0 ) as the spatial coordinates at the aperture plane perpendicular to the optical axis. Following the simple analysis in [10], it is direct to show that the coordinates (x0 , y 0 ) are the Fourier conjugate coordinates of (x, y). In practice, when using an infinity-corrected objective lens and tube lens to form an image, the Fourier conjugate aperture plane is located at the microscope objective lens back focal plane. Equivalent descriptions of microscopes sometimes rely upon the image of the aperture plane from the point of view of the image plane (i.e., the exit pupil). We define the image of the aperture function a(x0 , y 0 ) from the point of view of the image plane, as the pupil function, p(x0 , y 0 ). Typically, p(x0 , y 0 ) is equivalent to or a scaled version of the aperture function, a(x0 , y 0 ). In this thesis, we treat the aperture and pupil functions as equivalent. Both are Fourier conjugate to the sample and image planes. However, care should be taken in actual system analysis to ensure equality when appropriate. 1.1.4 The image plane After passing through the microscope, our optical field of interest will terminate at the image plane. In all imaging setups considered here, our image plane contains a digital detector array with a finite pixel size. Each pixel will sample the intensity of the incoming optical field. This sampling process can be described as a convolution between the incident optical field, which is band-limited, and a pixel sampling function [16]. Through Shannon’s sampling theorem, it is thus possible to exactly describe the intensity of the incoming optical field with the set of discrete measurements from the detector array. A detailed analysis of the effects of pixel sampling on image formation is given in Chapter 7 of [16]. Typical CCD and CMOS detectors contain approximately 106 -107 pixels, with pixel sizes that are in the range of δx = 1.5 − 10 µm. To avoid aliasing, we will always assume that the pixel size δx exceeds the maximum spatial frequency of incoming (coherent) light. Unless otherwise mentioned, our experiments use a 20 megapixel Kodak KAI-29050 CCD detector with a pixel size of 5.5 µm. 1.2 Relevant definitions This thesis repeatedly characterizes the performance of microscope imaging using several related parameters. In this section, we briefly define these common parameters for quick reference. Most of these parameters are also diagrammed in Fig. 1.3. 7 Diagram  to  show  FOV,  NA  and  PSF  (and  magnifica/on?)   (a) Geometric optics parameters! LEDs! h! Sample ! Aperture! Image! FOV! Field   number   θa! θi! l do! di! (b) Wave optics parameters! a(x') = rect(x'/w)! h(x) = sinc(xŸλ/NA)! CTF! PSF! Coherence length Lc ! LED area a! x'! x! z! x! z! z! Figure 1.3: (a) Relevant geometric optics parameters labeled for the simple microscope. (b) Wave optics parameters. Within the framework of wave optics, we may define a coherence length for our LED source, as well as a coherent transfer function (CTF) and point-spread function (PSF) for our microscope system. 8 • Magnification: In an infinity-corrected microscope setup, such as those considered in this thesis, the image magnification M is typically defined as, M = fo /ft , where fo is the objective lens focal length and ft = 180 mm is the tube lens focal length. To ensure that this standard is followed by our simple microscope, and that the well-known single lens magnification definition M = di /do also holds true, we simply set do = fo and di = fl = 180 mm. We maintain these equalities for the entire thesis. • Field number and field-of-view: The field number (FN) of a microscope is defined as the diameter of the measured optical field at the intermediate image plane of a compound microscope. In most stand-up microscopes, this corresponds to the plane of the digital detector. We assume the optical field at the image plane extends a large distance along x and y, such that the FN in our simple microscope is limited by the width of the digital detector. A typical field number is FN= 26.5 mm. The field-of-view (FOV, also called field size) is the field number de-magnified, which is its spatial extent at the image plane: FOV=FN/M . • Numerical aperture (NA): The microscope objective numerical aperture is NA = n·sin(θa ), where n is the refractive index of the medium of propagation (n = 1 in air), and θa is the maximum acceptance half-angle of the microscope objective. In our simple microscope, we define the maximum lens acceptance half-angle θa along the optical axis. We note that our LED array also has an effective illumination numerical aperture, N Ai , defined using the maximum possible angle θi of a plane wave generated by the most laterally displaced LED: NA = n·sin(θi ) • Impulse response and point-spread function: The impulse response of a coherent imaging system, h(x, y), defines its complex response to an idealized point source placed along the optical axis at the sample plane. In the absence of aberrations, the shape of h(x, y) remains shift-invariant across the entire image plane. The function h(x, y) is given by the inverse Fourier transform of the CTF, defined below. Typically, the CTF is a circular function in two dimensions, and we find h(x, y) = Jinc(xλ/N A). Here, the “Jinc” function is defined as Jinc(x) = J1 (2πx)/x, where J1 is a Bessel function of the first kind, order-1. A good rule of thumb is to assume the width of the impulse response as approximately λ/N A. As we detail in Chapter 2, this width approximates the smallest resolvable sample feature when imaging with a conventional microscope. Finally, the point-spread function of an imaging system is the 2 squared magnitude of the impulse response, |h(x, y)| . • Coherent transfer function (CTF): In a coherent imaging system, such as those considered in a large part of this thesis, the CTF defines the imaging system’s response in the spatial frequency domain. Specifically, given a sample function ψ(x, y), we may de- fine its spatial frequency domain representation as ψ̂(x0 , y 0 ), where again (x0 , y 0 ) are the 9 Fourier conjugate variables of (x, y), and the hat denotes a two-dimensional Fourier transform. Following principles from Fourier optics, the spatial frequency representation of the optical field at the image plane (i.e., its Fourier transform, or spectrum), ĝ(x0 , y 0 ), is given as, ĝ(x0 , y 0 ) = ψ̂(x0 , y 0 )ĥ(x0 , y 0 ). Here, ĥ(x0 , y 0 ) is the imaging system CTF. A fundamental insight from Fourier optics states that the CTF is simply a scaled version of the lens aperture function: ĥ(x0 , y 0 ) = a(λdi x0 , λdi y 0 ) [10]. Typically, we will ignore constant coordinate scaling factors when examining the behavior of our Fourier ptychographic microscope, for simplicity. This allows us to define the coordinates of the aperture plane as simply (x0 , y 0 ), the Fourier conjugate variables of the coordinates at the sample plane. Furthermore, we’ll also ignore any scaling effects between the sample and image plane and will label both using (x, y), as shown if Fig. 1.3. • Space-bandwidth product: The space-bandwidth product of an imaging system (SP) is the total number of resolvable features it can capture. Mathematically, the SP is approximately given by the imaging system field-of-view divided by the average width of its impulse response, h(x, y). If the width of the impulse response varies across the image plane, as is typically the case in most microscope objectives (due to the influence of aberrations), then this variation must be taken into account when computing the SP. This thesis will also occasionally use a related alternative definition of the SP of a band-limited optical system, given as the product of the spatial extent and spatial frequency range that it can fully capture [22]. • Quasi-monochromatic: As first noted in [23], the quasi-monochromatic condition must be met if one wishes to neglect the effect of the finite spectral bandwidth of an optical field. Specifically, if we wish to accurately assume that an optical source used within an imaging experiment is one frequency ν, then its spectral bandwidth ∆ν must satisfy the following inequality: ν/∆ν > nx . Here, nx is the number of pixels along one axis of the digital detector. If a source meets this condition, then one can neglect the effects of its spectral bandwidth on the intensity values within each detected image. In our setup, a ∆ν of several nanometers is required to fulfill the quasi-monochromatic condition, which a highly temporally coherent LED may satisfy. Unless otherwise stated, this thesis treats each LED in our illumination array as quasi-monochromatic. • Coherence length: While the above quasi-monochromatic assumption offers an accurate model for the spectral response of our LED array microscope, the assumption that each LED is an ideal point source is not accurate, in practice. Instead, we must take into account the finite width wl of each LED, which we treat as square light source that is fully incoherent within its photon generating region. From the Van Cittert-Zernike theorem, a fully incoherent source of finite width wl will emit a statistical field that will gain coherence upon propagation. Space-bandwidth product, Olympus microscope objectives! 10 Objective Resolution Space-Bandwidth Product (SP) Magnification/NA/Field number 532 nm incident wavelength (μm) 1.25X/0.04/26.5 8.12 21.5 MP 2X/0.08/26.5 4.06 33.5 MP 4X/0.16/26.5 2.03 33.5 MP 10X/0.3/26.5 1.08 18.9 MP 20X/0.5/26.5 0.65 13.1 MP 40X/0.75/26.5 0.43 7.4 MP 60X/0.9/26.5 0.36 4.7 MP 100X/1.3/26.5 0.25 3.5 MP megapixels Figure 1.4: Space-bandwidth product (SP) of various microscope objectives. Although maximum resolutions vary significantly, all lenses exhibit a SP less than 50 megapixels It is useful to define a measure of the field’s coherence [24], given by the coherence length Lc a distance z away from the LED, as Lc = λz/wl . Within this coherence length, a partially spatially coherent field remains effectively correlated with itself, and can thus be approximated as a coherent field. In other words, a two-slit experiment will produce visible fringes up to a slit separation of approximately Lc , after which the resulting fringes will decrease in contrast to zero. 1.3 Challenges for the standard microscope To capture a standard microscope image, one selects any of the illumination schemes from Fig. 1.2, illuminates the sample, forms an image of the illuminated sample on the detector, and records the optical intensity for a finite exposure time. While commercial systems can form very sharp images across a variety of magnifications, these snapshots still lack several key properties desired by the experimentalist. Below, we outline several of these key properties that Fourier ptychography can help solve, as well as some of the prior work that attempts to achieve a similar goal: • A high space-bandwidth product (SP): As detailed in [2], the SP of an imaging system is primarily influenced by two phenomena: diffraction and aberrations. Diffraction effects are minimized by using a larger lens (i.e., the PSF width scales inversely with the lens diameter and numerical aperture). Unfortunately, the size of aberrations also increase linearly with size of the lens. Thus, a tradeoff space emerges, where lenses must be highly optimized to both achieve a sharp PSF and also provide minimal aberrations across as large a FOV as possible. This optimization process is the job of the lens designer, who combines up to 20 individual 11 glass elements to form one microscope objective. Due to this lens scaling law, a tradeoff space eventually emerges between image resolution and field-of-view. Microscope objectives designed to capture very high resolution images suffer from a narrow field-of-view, and objectives designed to image a wide field-of-view exhibit poor resolution. Thus, all microscope systems are currently limited to a SP of approximately 50 megapixels (or less). Fig. 1.4 lists a variety of different objective lenses and the total number of pixels (i.e., the SP) that each can capture [26]. As described in Chapter 2, Fourier ptychography overcomes this tradeoff by first capturing a sequence of SP-limited images. Then, it fuses them together to reconstruct an image with an effective SP of 1 gigapixel. This lens scaling law is not limited to microscopes, but applies to all generalized imaging systems. Several recent efforts use the principle of multi-scale lens design to overcome this scaling law [27]. Unlike such prior work, Fourier ptychography requires variable illumination to achieve its resolution gain, and is thus best suited to increase the SP within microscopes, where one typically has control over the illumination source. • Quantitative phase: Light, as an electromagnetic field, has both an amplitude and phase. All optical detectors can only directly measure the amplitude of the field (specifically, the intensity) and not the phase. As mentioned in the introduction, optical phase can be extremely helpful. However, these phase-sensitive techniques from the early years of microscope design (e.g., phase contrast and differential interference contrast) do not offer quantitative measurements, but instead provide indirect evidence of phase variation. Quantitative phase directly measures variations in sample thickness and index of refraction. It also enables direct removal of system aberrations, and the ability to digitally refocus an image formed with spatially coherent light (see below). Measuring quantitative phase has a long history. Holography inherently relies upon recording information connected to the phase of the sample of interest. Digital holography, like well known in-line phase shifting technique [28], allows exact recovery of quantitative phase (up to a constant unknown phase offset) through a sequence of four images. Computational methods, such as phase retrieval [2,4], offer a means to estimate the phase of a sample without requiring a reference path and wavelength-precise shifting, but cannot always guarantee a completely accurate solution. Finally, the transport of intensity equation (TIE) offers another useful means to estimate phase [9, 31, 32]. However, TIE requires motion of the sample or within the imaging system, and typically operates with transparent samples. As we show in Chapters 2-4, Fourier ptychography acquires the quantitative phase of the optical field emerging from a sample, simultaneous to improving the imaging system spatial resolution. It may also simultaneously measure the quantitative phase deviations across the aperture plane (i.e., the 12 microscope aberrations), which can be removed to form a much sharper image. • Digital refocusing: In a standard microscope, focusing to the ideal plane of interest can be a challenge. Often, the sample is not perfectly flat at different areas across the sample plane, causing the image to come into focus at different depths along the optical axis, z. Frequent users of microscopes, such as pathologists, are well aware of the challenges surrounding manual focusing. One computational imaging technique, termed light field microscopy, replaces manual refocusing with limited digital image refocusability, post-capture [5, 34]. However, light field images offer a significantly limited resolution, and require the insertion of a microlens array near the image plane, thus necessitating a customized microscope body. Once one acquires both the amplitude and phase of the optical field at the detector, it is simple to digitally refocus the field through a wide range of axial planes, by computational propagation (e.g., using the angular spectrum method [10]). Computational propagation enables direct refocusing of FP images while maintaining their sub-micron resolution. In addition, it is common for a sample to occupy an unknown defocus plane, or perhaps lie at an unknown tilt. During FP reconstruction, we can additionally solve for these unknown placements and tilts to offer a sharp image across a significantly extended depth-of-field (up to 75 times a comparable objective lens, as discussed in Chapter 3). • Aberration removal: As noted above, all lenses exhibit aberrations, which deteriorate image quality. Over the past half-century, many unique aberration characterization methods have been reported [36–38]. These methods typically attempt to estimate the phase deviations or the frequency response of the optical system under testing. Several relatively simple noninterferometric procedures utilize a Shack-Hartmann wavefront sensor [38], consisting of an array of microlenses that each focus light onto a detector. Despite offering high accuracy, measuring aberrations with a Shack-Hartmann sensor often requires considerable modification to an existing optical setup. Alternatively, wavefront aberrations can be inferred directly from intensity measurements, by relying upon phase retrieval procedures [24]. Such computational methods typically require capturing multiple images, while inducing some unknown change to the optical system between each capture (i.e., while applying “measurement diversity” [24,40]). As we will show in this thesis, Fourier ptychography provides a very useful form measurement diversity through its variable illumination source. This enables robust computational estimation of system aberrations and their subsequent removal, which significantly improves the effective resolution of our final image reconstructions. • 3D structural information: Typical microscope images are two-dimensional. As we detail in Chapter 8, it is possible to computationally invert an FP dataset into a three dimensional 13 representation of a thick sample. This inversion process closely resembles diffraction tomography (DT) [14]. Unlike in all prior DT setups, FP does not measure the phase of the optical field. This allows one to use a standard microscope, outfitted with an LED array and computer, to computationally recover a quantitative measurement of the complex index of refraction of a thick sample of interest (t > 100µm), from throughout its volume, at a micrometer-scale resolution. There is one key item that many current microscope experiments currently focus on that is missing from the above list: improving fluorescent image capture. For example, recent fluorescent particle localization techniques, such as PALM and STORM, can exceed standard microscope resolution limits [42]. Likewise, structured illumination methods may yield similarly unbounded resolution gains [18]. Standard Fourier ptychography assumes all optical interactions are coherent (i.e., a phase relationship is preserved, at least locally). Fluorescent excitation does not preserve phase and thus will not obey the interactions FP predicts. While it is possible to adapt the principle of Fourier ptychography to achieve fluorescent resolution enhancement [43], this thesis will not discuss in detail this possible direction. We hope the reader keeps this important point in mind throughout the following discussions. Here is an outline for the rest of this thesis. In Chapter 2, we will overview the principle of FP in an optical microscope, from a Fourier optics perspective. In Chapter 3, we will discuss the process of Fourier ptychographic image reconstruction using a phase retrieval algorithm. In Chapter 4, we will examine the quantitative accuracy of Fourier ptychographic phase reconstruction. We also briefly visit several applications of measuring high-resolution phase, including the ability to digitally refocus images, remove aberrations, and determine structural information regarding biological tissue. In Chapter 5, we will connect Fourier ptychography to its “standard” counterpart in X-ray imaging, ptychography, using a unique mathematical phase-space model. This model highlights the role of partial coherence during data acquisition and image reconstruction. In Chapter 6, we will present a convex approach to process both standard and Fourier ptychographic data, which performs better in the presence of noise than prior reconstruction techniques. In Chapter 7, we implement Fourier ptychography within a conventional camera setup, which helps remove unknown aberrations from a final reconstructed image. Finally, in Chapter 8, we apply the principle of ptychography to reconstruct thick samples, using a process that resembles diffraction tomography but does not require the measurement of optical phase. 14 Bibliography [1] R. M. Allen, The Microscope (Chapman and Hall Limited, 1940). [2] F. Zernike, “Beugungstheorie des Schneidenverfahrens und seiner verbesserten Form, der Phasenkontrastmethode,” Physica 1, 689–704 (1934) [3] F. Zernike, “How I Discovered Phase Contrast”. Science 121, 345–349 (1955). [4] F. Zernike, “Phase contrast, a new method for the microscopic observation of transparent objects part I”, Physica 9 (7), 686–698 (1942). [5] G. Nomarski, “Interferential polarizing device for study of phase objects,” US Patent US2924142A (1960). [6] M. Minsky, “Microscopy apparatus,” US Patent US-3013467A (1961). [7] Z. Huang, “A Review of Progress in Clinical Photodynamic Therapy,” Technology in Cancer Research and Treatment, 4(3), 283–293 (2005). [8] G. C. R. Ellis-Davies, “Caged compounds: photorelease technology for control of cellular chemistry and physiology,” Nat. Meth. 4(8), 619–628 (2007). [9] F. Zhang et al., “Optogenetic interrogation of neural circuits: technology for probing mammalian brain structures,” Nat. Protocols 5(3), 439–456 (2010). [10] J. W. Lichtman and J. A. Conchello, “Fluorescence microscopy,” Nat. Meth. 2(12), 910–919 (2005). [11] S. W. Hell and J. Wichmann, “Breaking the diffraction resolution limit by stimulated emission: Stimulated-emission-depletion fluorescence microscopy,” Opt. Lett. 19(11), 780-782 (1994). [12] G. Zheng, R. Horstmeyer and C. Yang, “Wide-field, high-resolution Fourier ptychographic microscopy,” Nature Photon. 7, 739–745 (2013). [13] P. Ferraro, A. Wax and Z. Zalevsky, Coherent light microscopy: imaging and quantitative phase analysis (Springer, New York, 2011). 15 [14] K. Nugent, “Coherent methods in the X-ray sciences,” Adv. Phys. 59(1), 1–99 (2010). [15] J. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1996). [16] D. Brady, Optical Imaging and Spectroscopy (John Wiley & Sons, 2009). [17] M. Gustafsson, “Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy,” J. Microsc. 198, 82–87 (2000). [18] M. G. L. Gustafsson, “Nonlinear structured-illumination microscopy: wide-field fluorescence imaging with theoretically unlimited resolution,” Proc. Nat. Acad. Sci. 102, 13081–13086 (2005). [19] J. Garcia, Z. Zalevsky and D. Fixler, “Synthetic aperture superresolution by speckle pattern projection,” Opt. Express 13, 6073–6078 (2005). [20] E. Mundry et al., “Structured illumination microscopy using unknown speckle patterns,” Nat. Photonics 6, 312–315 (2012). [21] P. D. Nellist, B. C. McCallum and J. M. Rodenburg, “Resolution beyond the ‘infromation limit’ in transmission electron microscopy,” Nature 374, 630–632 (1995). [22] A. W. Lohmann, R.G. Dorsch, D. Mendlovic, Z Zalevsky, and C. Ferreira, “Spacebandwidth product of optical signals and systems,” J. Opt. Soc. Am. A 13, 470473 (1996). [23] A. W. Lohmann, “Matched Filtering with Self-Luminous Objects, Appl. Opt. 7, 561563 (1968). [24] J. Goodman, Statistical Optics (John Wiley and Sons, 2000). [25] A. W. Lohmann, “Scaling laws for lens systems,” Appl. Opt. 28, 4996-4998 (1989) [26] G. Zheng, X. Ou, R. Horstmeyer, J. Chung and C. Yang, “Fourier ptychographic microscopy: creating a gigapixel microscope for biomedicine,” Optics and Photonics News 25(4), 26–33 (2014). [27] D. Brady et al., “Multiscale gigapixel photography”. Nature 486, 386389 (2012). [28] R. Gerchberg, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik (Stuttg.) 35, 237 (1972). [29] J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21(15), 2758–2769 (1982). [30] I. Yamaguchi and T. Zhang, “Phase-shifting digital holography,” Opt. Lett. 22(16) 1268–1270 (1997). [31] N. Streibl, “Phase imaging by the transport equation of intensity,” Opt. Commun. 49(1), 6–10 (1984). 16 [32] T. E. Gureyev and K. A. Nugent, “Rapid quantitative phase imaging using the transport of intensity equation,” Opt. Commun. 133, 339–346 (1997). [33] L. Waller, L. Tian, and G. Barbastathis, “Transport of Intensity phase-amplitude imaging with higher order intensity derivatives,” Opt. Express 18(12), 12552–12561 (2010). [34] M. Levoy, R. Ng, A. Adams, M. Footer and M. Horowitz, “Light field microscopy,” ACM Trans. Graphics 25, 924–934 (2006). [35] M. Levoy, Z. Zhang and I. McDowall, “Recording and controlling the 4D light field in a microscope using microlens arrays,”. J. Microsc. 235, 144–162 (2009). [36] M. Takeda and S. Kobayashi, “Lateral aberration measurements with a digital Talbot interferometer,” Appl. Opt. 23(11), 1760–1764 (1984). [37] S. Gong and S. S. Hsu, “Aberration measurement using axial intensity,” Opt. Eng. 33(4), 1176– 1186 (1994). [38] J. L. Beverage, R. V. Shack, and M. R. Descour, “Measurement of the three-dimensional microscope point spread function using a Shack-Hartmann wavefront sensor,” J. Microsc. 205(1), 61–75 (2002). [39] R. G. Paxman, T. J. Schulz, and J. R. Fienup, “Joint estimation of object and aberrations by using phase diversity,” JOSA A 9(7), 10721085 (1992). [40] M. Guizar-Sicairos and J. R. Fienup, “Phase retrieval with transverse translation diversity: a nonlinear optimization approach,” Opt. Express 16, 7264–7278 (2008). [41] V. Lauer, “New approach to optical diffraction tomography yielding a vector equation of diffraction tomography and a novel tomographic microscope,” Journal of Microscopy 205, 165–176 (2002). [42] M.A. Thompson, M.D. Lew, and W.E. Moerner, “Extending microscopic resolution with singlemolecule imaging and active control,” Annu Rev Biophys 41, 321–42 (2012). [43] S. Dong, P. Nanda, K. Guo, and G. Zheng, “High-resolution fluorescence imaging via patternilluminated Fourier ptychography,” Opt. Express, 22 (17), 20856–20870 (2014). 17 Chapter 2 Fourier ptychography for gigapixel imaging In this chapter, we use our model of a simple microscope from Chapter 1 to explain the principle of Fourier ptychography (FP). The goal of FP is to computationally resolve an image containing one billion pixels (i.e., a gigapixel) from a sequence of low-resolution measurements. To begin, we will first motivate the concept of improving resolution by means of a shifting illumination source. Then we will present a mathematical model to explain the process of Fourier ptychographic data capture. We will limit our discussion here to a two-dimensional optical geometry (i.e., one spatial axis that is orthogonal to the axis of optical propagation). At the end of this chapter, we will discuss extension to a three-dimensional setup. 2.1 Resolution improvement through shifting illumination As noted in Chapter 1, a microscope’s numerical aperture (NA) defines its resolution performance. Specifically, the minimum resolvable feature within a standard microscope is approximated as λ/NA = λ/ sin θ. Here, we assume the standard microscope is under coherent illumination (from a distant point source) and operates within air, which sets the index of refraction term within the NA equal to 1. As diagrammed in Fig.2.1(a), the half-angle θ denotes the largest cone of wavevectors (i.e., rays) that can pass from the sample and into the imaging lens, and defines the system’s maximum resolvable feature. Most biological samples of interest contain many micrometer-scale features that diffract light into a large cone of wavevectors. In Fig.2.1(a), this larger green cone subtends an angle ω > θ and is thus not completely captured by the microscope. To resolve smaller features at the image plane, we would ideally like to capture all of the rays within this large green cone. Specifically, we would like to somehow keep the low-NA objective lens in Fig.2.1 for imaging, since it offers a very large field of view (FOV), but still collect the entire green cone of rays, which will simultaneously offer us 18 (a) Direct illumination! (b) Angled illumination! Microscope lens (red)! Illumination source! Shift illumination! Tilted ray cone! (c) LED illumination for multiple angles! Illuminate n LEDs! Capture n images! Detector! Sample! ω   θ" shift! ω   φ" φ" θ" Figure 2.1: Capturing a wider cone of wavevectors. (a) A standard microscope with a limited NA = n sin θ captures only the red cone of rays emerging at half-angle θ from the sample. An ideal microscope can capture the wide green cone of rays emerging at angle ω. (b) By shifting the illumination source, we can rotate the green cone of rays emerging from the thin sample surface. Now, a different segment of the green cone will pass through the fixed “red” aperture. (c) FP uses an array of n LEDs to shift the wide green cone of rays n times. Each time the cone is shifted, a unique image is captured. From the set of n images, the computational goal of FP is to synthesize an image with a larger NA, as if it originated from a lens that could originally capture the wide green cone of rays. high resolution. The combination of both a wide image FOV and high resolution will lead to a large space-bandwidth product (SP) system. To capture the large green cone of rays through our low-NA microscope objective lens, Fourier ptychography shifts its illumination source, as diagrammed in Fig.2.1(b). A shifting illumination source will illuminate the sample from a variable angle, φ. If the sample is thin, this will cause the green cone of rays emerging from the sample surface to also rotate, by the same angle φ. At the microscope lens, the green cone of rays will subsequently shift laterally by a finite distance. As is clear in Fig.2.1(b), the shifted source will cause a different segment of the green cone to pass through our fixed microscope objective lens (denoted by the red circle). As we shift the illumination source across many different positions, then the entire green cone of rays, albeit at different points in time, will pass through the fixed objective lens and subsequently propagate to the detector. If we capture an image of the optical field that passes through the objective lens at each unique location of the shifted optical source, then we have conceptually captured information corresponding to the entire green cone of rays. The goal of Fourier ptychography is to use this set of captured images to computationally recover a single image that appears to have passed through a “synthetic” lens, whose effective size extends across the entire green cone of rays (i.e., a new lens with an acceptance half-angle ω). This synthesized image will now have a much higher resolution – its smallest resolvable element will now approach λ/ sin ω. Instead of actually shifting an illumination source laterally, Fourier ptychography typically uses an LED to create its angular illumination (Fig.2.1(c)). Next, we develop a mathematical model to 19 (a) FP setup! LEDs! Sample ψ (x)! Aperture a (x')! Image g (x)! jth source! θa" θj" r" l L(x')! S(x)! A(x')! Sample! Shifted spectrum! exp(ikxpj)!Ÿ! ψ(x)! ψ(x-pj)! D(x)! (b) FP, wave optics! jth source! r" Plane wave! Image! pj   F [a(x)ψ(x-pj)]! 2! Figure 2.2: Optical diagram for our simple 2D FP setup, with (a) planes and distances of interest labeled, and (b) wave optics descriptions of the resulting optical field at each plane of interest. The final image is both a function of space (x) and incident illumination angle (pj ). Each captured image thus forms one column of our FP data matrix, g(x, pj ). solidify the above concept of Fourier ptychographic data capture, before explaining how we transform the data into a high resolution image, in Chapter 3. 2.2 Mathematical model of Fourier ptychography We mathematically model our microscope imaging system in two dimensions (x, z), for simplicity. We assume that a distant plane L(x0 ) contains q different quasi-monochromatic optical sources (central wavelength λ) evenly distributed along x0 with a spacing r. We assume each optical source acts as an effective point emitter that illuminates a sample ψ(x) at a plane S(x) a large distance l away from L(x0 ). Under this assumption, the jth source illuminates the sample with a spatially coherent plane wave at angle θj = tan−1 (jr/l), where −q/2 ≤ j ≤ q/2. Next, we additionally assume the sample ψ(x) is thin (see definition in Chapter 1). Under this assumption, we may express the optical 20 field exiting the thin sample as the sample-plane wave product, s(x, pj ) = ψ(x)eikxpj , (2.1) where the wavenumber k = 2π/λ and pj = sin θj describes the off-axis angle of the jth optical source. The jth illuminated sample field s(x, pj ) then enters an imaging system with a low numerical aperture (NA). Neglecting scaling factors and a quadratic phase factor for simplicity, Fourier optics gives the field at the imaging system aperture plane, A(x0 ), as F [s(x, pj )] = ψ̂(x0 − pj ). (2.2) Here, F represents the Fourier transform between conjugate variables x and x0 , and ψ̂ is the Fourier transform of ψ. Eq. 2.2 uses the Fourier shift property to show that the spectrum of a thin sample, when illuminated by a plane wave at an angle given by the inverse sine of pj , will shift by pj laterally across the imaging system’s aperture plane. The shifted spectrum field ψ̂(x0 − pj ) is then modulated by the aperture function of the imaging system, a(x0 ), which acts as a low-pass filter. As outlined in Chapter 1, the aperture function a(x0 ) is typically a physical stop placed at the back focal plane of the microscope objective. It usually exhibits a circular shape in two dimensions. In the conceptual schematic of FP in Fig. 2.2, a(x0 ) is the limited width of the transparent lens area at the aperture plane, which may be mathematically expressed as a rect function. The shape of a(x0 ) is directly proportional to the imaging system coherent transfer function (CTF), which in turn defines the lens cutoff spatial frequency, its maximum acceptance angle, and thus also its NA. It is now useful to consider the sample spectrum ψ̂ discretized into n pixels with a maximum spatial frequency k. Following the Nyquist-Shannon sampling theorem, since our optical signal is now bandlimited to the aperture area a(x0 ), this discretization process is exact. We denote the bandpass cutoff of the aperture function a as k · m/n, where m is an integer less than n. The modulation of ψ̂ by a results in a field characterized by m discrete samples, which propagates to the camera imaging plane, D(x). Here, the digital detector critically samples the incident field using an m-pixel digital detector (e.g., a CCD or CMOS detector, which we assume has a perfect fill factor). This forms a reduced-resolution intensity image, g(x, pj ): h i2 g(x, pj ) = F a(x0 )ψ̂(x0 − pj ) . (2.3) Here, we use the subscript j to denote that this image is formed by the jth LED. By turning on one LED at a time and saving the resulting image, for −q/2 ≤ j ≤ q/2, we will compile an (m×q) Fourier 21 Shift jyŸr LED Array θj Shift jxŸr Sample ψ (x,y) Aperture a (x',y' ) Image g (x, y, pjx, pjy ) L(x',y') S(x,y) A(x',y') Fourier Ptychography Setup Detector D(x,y) High-resolution complex sample ψ Low-resolution image set g Fourier Ptychography … 0! 2π! intensity phase … intensity Figure 2.3: Optical diagram for an example FP setup in three dimensions, including all of the same parameters and planes as shown in Fig. 2.2. (bottom) Actual raw data and resulting reconstruction of a resolution target, demonstrating a 25X increase in microscope resolution and the simultaneous acquisition of the optical phase from the sample plane. ptychography data matrix, g(x, pj ). The jth column of this data matrix contains a low-resolution image of the sample intensity while it is under illumination from the jth optical source. The goal of Fourier ptychographic post-processing is to reconstruct a high-resolution (n-pixel) complex spectrum, ψ̂(x0 ), from the multiple low-resolution (m-pixel) intensity measurements contained within the data matrix g(x, pj ). Once ψ̂ is found, an inverse-Fourier transform will yield the desired complex sample reconstruction, ψ. As we will explain in the next chapter, most current ptychography setups solve this inverse problem using an alternating projections (AP) algorithm: 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 set of measured intensities from the detector plane (i.e., the values in g). Second, its spectrum ψ̂0 is forced to lie within a known support in the plane that is Fourier conjugate to the detector. Different projection operators and update rules are available, but are closely related [4, 5, 46]. It is also possible to solve the inverse ptychography problem with a convex program, which this thesis examines in detail in Chapter 6. Before examining possible methods to reconstruct the sample at high resolution from the data 22 Synthesized  spectrum   es2mate,    ψ0 (x')   Image 1! Detect 3 complex images:! g' (x',1-3)! Image 2! Image 3! Inverse Fourier transform! LEDs! Sample! Aperture! Figure 2.4: Conceptual outline of using a synthetic aperture to improve image resolution. If each detected image in the sequence of low-resolution measurements also captures the phase of the optical field, then reconstruction is direct. Each detected image is inverse Fourier transformed to the aperture plane (to form the values of the optical field within each colored circle). These three Fourier transforms are then stitched together, side-by-side (to form the large green circle). An inverse Fourier transform of this synthesized spectrum yields the desired high-resolution complex image. Unfortunately, the phase of light is not detected by standard optical pixel arrays, which requires us to use a slightly more involved high-resolution image reconstruction procedure, as outlined in the next chapter. matrix in Eq. 6.2, we note that the above discussion easily generalizes to a three-dimensional geometry. This extension is diagrammed in Fig. 2.3, where all functions of x and x0 are now also functions of y and y 0 , respectively. LED scanning is performed along two dimensions separately. Combined with two-dimensional images, this creates a four-dimensional data matrix, g(x, y, pxj , pyj ), where now r/2 ≤ jy ≤ r/2 is the counter for LEDs along the second axis, where a total of q × r LEDs are in the possibly rectangular array. 2.3 Reconstruction goals and extensions The goal of ptychographic recovery is to convert the acquired data matrix into a complex sample reconstruction with improved resolution. It will be useful at several points throughout this thesis to return to the above data matrix picture, as several different algorithms are discussed. Here, we briefly outline one related strategy to help connect the above concepts to our final goal of resolution improvement. First, let us assume for simplicity that our digital detector can also detect the phase of incident light (note this is typically not true in practice, especially within a standard microscope). However, 23 if the phase were detected, then our data matrix would no longer be squared: h i g 0 (x, pj ) = F a(x0 )ψ̂(x0 − pj ) . (2.4) Furthermore, let us assume that each captured image originated from a unique but adjacent (i.e., nonoverlapping) region of the spectrum, ψ̂(x0 ). In other words, we assume that our LEDs are arranged such that the angular shift they impart to the spectrum at each step, pj − pj−1 , is equal to the total extent of the aperture in the Fourier domain, k ·m/n. This ideal condition is diagrammed in Fig. 2.4. Under these two assumptions, the recovery of a high resolution image from our measurements in the data matrix would require three straightforward steps: • Inverse Fourier transform each image in the complex data matrix in Eq. 2.4 to recover each associated spectrum “tile”, ĝ 0 (x, pj ) = a(x0 )ψ̂(x0 − pj ). • Form one long spectrum estimate vector, ψ̂0 (x0 ), by arranging each spectrum tile from the data matrix adjacent to one another. For all j, take the spectrum tile a(x0 )ψ̂(x0 − pj ), and place it in ψ̂0 starting at entry (j − 1) · k · m/n, and ending at entry j · k · m/n. This spectrum synthesis process is diagrammed in Fig. 2.4. • After all spectrum tiles have been concatenated together to form ψ̂0 (x0 ), take the inverse Fourier transform of ψ̂0 (x0 ) to recover the high resolution complex sample estimate, ψ0 (x). This simple inversion process, in the absence of noise, allows for exact recovery of a high resolution image from a sequence of low-resolution complex image measurements. The entire process of low resolution capture and high resolution synthesis is commonly referred to as “synthetic aperture” imaging. Many examples of synthetic aperture imaging arise in the area of remote sensing, which often operates within the microwave and radio regimes, where the phase of the incident field is easily measured on a digital detector. It is also possible to implement the same concept in optical setups, using holographic techniques to determine the phase of each individual low-resolution image [4, 5, 7, 8, 8–16]. Since it is not possible to directly measure phase within a standard microscope, synthetic aperture imaging will not help us directly achieve our goal of gigapixel image formation. In the next chapter, we detail how Fourier ptychography applies a phase retrieval algorithm to overcome this barrier and offer high-resolution, complex sample maps from a series of low-resolution intensity measurements. 24 Bibliography [1] J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21(15), 2758–2769 (1982). [2] S. Marchesini, “A unified evaluation of iterative projection algorithms for phase retrieval,” Rev. Sci. Instrum. 78, 011301 (2007) [3] P. Godard, M. Allain, V. Chamard and J. Rodenburg, “Noise models for low counting rate coherent diffraction imaging,” Opt. Express 20 25914–25934 (2014). [4] T. Turpin, L. Gesell, J. Lapides and C. Price, “Theory of the synthetic aperture microscope”. Proc. SPIE 2566, 230–240 (1995). [5] J. Di et al., “High resolution digital holographic microscopy with a wide field of view based on a synthetic aperture technique and use of linear CCD scanning,” Appl. Opt. 47, 5654–5659 (2008). [6] T. R. Hillman, T. Gutzler, S. A. Alexandrov and D. D. Sampson, “High-resolution, wide-field object reconstruction with synthetic aperture Fourier holographic optical microscopy”. Opt. Express 17, 7873–7892 (2009). [7] L. Granero, V. Mico, Z. Zalevsky and J Garcia, “Synthetic aperture superresolved microscopy in digital lensless Fourier holography by time and angular multiplexing of the object information,” Appl. Opt. 49, 845–857 (2010). [8] M. Kim et al., “High-speed synthetic aperture microscopy for live cell imaging,” Opt. Lett. 36, 148–150 (2011). [9] C. J. Schwarz, Y. Kuznetsova and S. Brueck, “Imaging interferometric microscopy,” Opt. Lett. 28, 1424–1426 (2003). [10] P. Feng, X. Wen and R. Lu, “Long-working-distance synthetic aperture Fresnel off-axis digital holography,” Opt. Express 17, 5473–5480 (2009). [11] V. Mico, Z. Zalevsky, P. Garcia-Martinez and J. Garcia, “Synthetic aperture superresolution with multiple off-axis holograms,” J. Opt. Soc. Am. A 23, 3162–3170 (2006). 25 [12] C. Yuan, H. Zhai, and H. Liu, “Angular multiplexing in pulsed digital holography for aperture synthesis,” Opt. Lett. 33, 2356–2358 (2008). [13] V. Mico, Z. Zalevsky, and J. Garcia, “Synthetic aperture microscopy using off-axis illumination and polarization coding,” Opt. Commun. 276, 209–217 (2007). [14] A. E. Tippie, A. Kumar and J. R. Fienup, “High-resolution synthetic-aperture digital holography with digital phase and pupil correction,” Opt. Express 19, 12027–12038 (2011). [15] T. Gutzler, T. R. Hillman, S. A. Alexandrov and D. D. Sampson, “Coherent aperture-synthesis, wide-field, high-resolution holographic microscopy of biological tissue,” Opt. Lett. 35, 1136–1138 (2010). [16] S. A. Alexandrov, T. R. Hillman, T. Gutzler and D. D. Sampson, “Synthetic aperture Fourier holographic optical microscopy.,” Phys. Rev. Lett. 97, 168102 (2006). 26 Chapter 3 Fourier ptychographic image reconstruction In Chapter 2, we derived an expression for the data matrix in Fourier ptychography (see Eq. 6.2). We also saw that if the phase of the data matrix elements is known, then a synthetic aperture technique can recover a high-resolution microscope image from its set of low-resolution measurements. Since the phase of each low-resolution measurement is not known, the inverse problem of resolution enhancement becomes more challenging. In this chapter, we detail how to use an iterative phase retrieval algorithm to accurately solve this inverse problem. 3.1 The phase retrieval algorithm The problem of recovering a discretized complex signal (i.e., its amplitudes and phases) from knowledge of just its amplitudes has a long history. Within the context of optics, one of the first successful attempts to solve this problem was initiated by Gerchburg and Saxton in the 1970’s [2]. In this section, we outline their general computational strategy, described within the context of our microscope. We refer to this general class of algorithm as a “phase retrieval” method. Here, we explain the simple “error reduction” (ER) phase retrieval algorithm [4], although one of many related strategies may also be used [46]. Phase retrieval often iteratively projects an initial estimate of the unknown complex sample, which we will define as hk (x, y) = h0 (x, y), onto two constraints in two different domains. Here, the subscript k indicates the estimate is in the kth iteration of our iterative solution process. One constraint is always the measured signal amplitudes, which in the case of our microscope is at the image plane (i.e., the square root of the measured intensities in the spatial domain). To begin, let us consider just one image from our microscope, |g(x, y)|, which results from illuminating the sample from directly below with the centered LED. To determine the sample’s unknown phase values, it is common in coherent imaging to use as a second constraint the assumption of a finite sample spectrum support. This second constraint 27 Conventional phase retrieval! Aperture plane! Initialize:! (1)   Support constraint, A! h0(x,y)! hk(x',y')! (4)   Spectrum Image plane! Fourier transform! (2)   F [hk]! Inverse Fourier transform! estimate! h'k(x,y)! Image estimate! hk(x,y)! (3)  Amplitude constraint, |!g!|! F-1 [h'k]! h'k(x,y)! Figure 3.1: The ER phase retrieval algorithm for coherent imaging. Each step is detailed in the text. constitutes additional a-priori signal knowledge, and is needed for our algorithm to accurately solve the inverse phase retrieval problem. When imaging, one often knows the shape and extent of the imaging system aperture function, a(x0 , y 0 ), which is almost always zero past a certain cutoff frequency. Thus, one has a-priori knowledge that the Fourier transform of the unknown signal estimate, ĥk (x0 , y 0 ) must be zero outside a certain region, or in other words must be a band-limited signal. Occasionally (but not always), knowledge of the signal amplitudes and the extent of its band-limited support is sufficient to accurately recover the unknown optical phase at the image plane. ER phase retrieval uses both the first “amplitude” and the second “support” constraint to iteratively encourage an initial signal estimate to converge to a solution containing the correct phase values (see diagram in Fig. 3.1). First, ER initializes an estimate of the complex sample spectrum, ĥ0 (x0 , y 0 ), at the aperture plane. Second, ER digitally propagates this estimate to the image plane. Following the description of our microscope from Chapter 2, propagation to the image plane is given by a Fourier transform: hk (x, y) = F[ĥk (x0 , y 0 )]. Again, k here denotes the kth iterative loop, for 0 ≤ k ≤ n iterations, and we continue to use the “hat” to denote the signal spectrum at the aperture plane. Next, ER enforces the amplitude constraint at the image plane. It replaces the amplitudes of hk (x, y) with the experimentally measured amplitudes at the detector, |g(x, y)|: 0 hk (x, y) = |g(x, y)| hk (x, y) . |hk (x, y)| (3.1) We may equivalently write our estimate in terms of its amplitudes and phases, hk (x, y) = α(x, y)eiφ(x,y) , 28 0 and then represent the estimate update step in Eq. 3.1 as, hk (x, y) = |g(x, y)| eiφk (x,y) , which makes clear that the original phase values of the estimate, φk (x, y), remain unchanged. Third, ER 0 propagates this amplitude-constrained hologram estimate back to the sample plane: ĥk (x0 , y 0 ) = h 0 i F −1 hk (x, y) . Fourth, ER applies the sample support constraint. It leaves unchanged all values within a defined subset of pixels, A, representing the nonzero locations within the aperture function bandpass, a(x0 , y 0 ). However, it assumes that outside of this area of interest the spectrum must be zero: 0 ĥk+1 (x, y) =   ĥ0 (x, y), k (x, y) ∈ A  0, (x, y) ∈ / S. (3.2) In this last step, ER also increments the iteration counter value, k, to begin the next iteration. It again transforms the estimate to the image plane, constrain its amplitudes, etc. Typically, the above alternating projection loop runs for a fixed number of iterations q, or until some convergence criteria is satisfied. The complex algorithm output, hq (x, y), typically offers an accurate estimate of the amplitude and phase of the original optical field at the image plane. 3.2 The Fourier ptychographic microscopy (FPM) algorithm The reconstruction algorithm used by Fourier ptychographic microscopy (FPM) effectively combines the synthetic aperture strategy outlined at the end of Chapter 2, with the phase retrieval algorithm outlined above. A detailed presentation of this algorithm and its implementation within a microscope is presented in [32]. Here, we offer an alternative summary, now examining just a one-dimensional signal (along x) for simple notation. After capturing the data matrix in Eq. 6.2, we start the FPM recovery algorithm by initiating a high-resolution complex spectrum estimate, ŝ(x0 ). Here, the structure of ŝ(x0 ) is not critical: it may be all zeros, or proportional to the Fourier transform of one of the raw images. The resolution of ŝ(x0 ) must match the final desired high-resolution after reconstruction, n, which may be computed via knowledge of the number of LEDs and known discrete distance of aperture shift between each LED. Specifically, if each low-resolution image contains m pixels, then we may use the numerical aperture of the objective lens (NAo ) and the illumination (NAi ) to find that, n=m× N Ao + N Ai . N Ao (3.3) It is worth noting here that a critical component of the success of Fourier ptychography lies within the correct selection of n with respect to the number of captured images, q. Unlike the synthetic aperture scenario, FPM requires that the regions of the spectrum that correspond to each captured image overlap to a certain extent in Fourier space. To understand this in a bit more detail, it is helpful to rewrite the jth field at the aperture plane as a fixed spectrum, at full-resolution and 29 FPM recovery algorithm! Image plane! Aperture plane! Initialize:! Select j th aperture region:! s0(x)! sk(x') a(x'-pj)! Fourier transform! F [skaj]! Sub-Image estimate! hk(x,j)! Increment j! Replace within aperture region! sk(x') a(x'-pj) = !h'k(x,j)! Update of subspectrum! h'k(x,j)! Inverse Fourier transform! F-1 [h'k]! Amplitude constraint, |!g!|! h'k(x,j)! Figure 3.2: The FPM algorithm. Each step is detailed in the text. containing n pixels, windowed by a shifted aperture: a(x0 + pj )ψ̂(x0 ). Note this is equivalent to the actual physical scenario of a shifted spectrum, reflected in Eq. 2.2. If the finite bandpass of a(x0 ) is m pixels wide, then FPM requires that pj ≤ m/2, such that a(x0 + pj−1 ) and a(x0 + pj ) each mask a similar finite portion of ψ̂(x0 ), or in other words overlap with one another. We define the percentage of aperture overlap as, ol = 1 − (n − m)/qm, where one may interpret m as the width of a and n as the width of ŝ. Typically, an overlap percentage of ol ≥ 60% is needed for successful FPM algorithm convergence. Given we satisfy the above overlap criterion, by carefully selecting the separation distance of our LEDs, the first step of FPM recovery starts with j = 0 and runs until j = q (the total number of captured images, see Fig. 3.2). We select the jth region of the spectrum, corresponding to the jth image, by computing the product, ŝk (x0 ) × a(x0 + pj ). Here, we have again shifted the aperture function a across the spectrum by a distance pj from the origin, which is mathematically equivalent to a shifted spectrum ŝ across a fixed aperture (i.e., what our FPM implements via the LED array). Then, we Fourier transform this limited region of the spectrum to the image plane: hk (x, j) = F [ŝk (x0 ) × a(x0 + pj )] . (3.4) Eq. 3.4 uses our high-resolution estimate to simulate the formation of an microscope image under illumination from the jth LED. We label hk also with the variable j to maintain its association with the jth image. Second, we update the amplitude of this simulated image with the known amplitude from the jth measured image, g(x, pj ). This is the same as the amplitude constraint used in the ER 30 Spectrum estimate, hk(x',y') F [Ÿ]   Fourier transform of region j Ÿ eiϕ(x,y)! y'! j-1 …   α(x,y)! j …   F-1 [Ÿ]   Replace amplitudes with j th image Ÿ eiϕ(x,y)! x'! Ÿ eiϕ(x,y)! α(x,y)! Figure 3.3: Diagram of FPM operation. Each local sub-spectrum region is Fourier transformed to the image plane (to the right, one at a time). The jth image is then updated with the measured intensities from the image corresponding to illumination from the jth LED. After update, it is inverse Fourier transformed back and inserted into the sub-spectrum region it originated from. phase retrieval algorithm: 0 hk (x, j) = |g(x, pj )| hk (x, j) . |hk (x, j)| (3.5) Here, the phase of hk is left unchanged. Third, we inverse Fourier transform this updated image h 0 i 0 back to the aperture plane: ĥk (x0 , j) = F −1 hk (x, j) . This forms an update of the local region of the spectrum that we assume passed through the aperture when we turned on the jth LED. Fourth, we locally update our spectrum estimate ŝk (x0 ) with the values in the jth updated 0 spectrum region, ĥk (x0 , j). During this update, we only modify spectrum estimate values that lie within the finite bandpass (that is, support) of the jth shifted aperture. Mathematically, this support update takes the form, 0 ŝk (x0 ) = (1 − a(x0 + pj )) × ŝk (x0 ) + a(x0 + pj ) × hk (x0 , j) (3.6) We assume in Eq. 3.6 that a(x0 + pj ) is 1 within its finite passband and 0 outside. If this is not the case, a similar expression may still yield a correct update. Finally, we increment our counter j and select a new spectrum region, corresponding to the values of the spectrum that passed through the microscope aperture when illuminating the (j + 1)th LED. This begins our process of sub-spectrum region update again, which proceeds until j = q. Then, we increment k, again looping through all sub-spectrum regions and updating them with the amplitudes within each image. This process is terminated typically after several iterations of k, or until reaching some threshold error metric value. A diagram of this update process for a 2D spectrum is in Fig. 3.3. 31 3.3 Experimental demonstration of FPM To demonstrate FPM’s ability to improve microscope resolution beyond the standard cutoff frequency, we construct the experimental microscope setup shown in Fig. 3.4 [32]. This FP microscope is simply a conventional microscope body, modified by placing a 15 × 15 array of surface-mounted LEDs (model SMD 3528, center wavelength λ = 632 nm, 4 mm LED pitch, 150 µm active area diameter), which serve as our quasi-coherent optical sources. The array is placed l =80 mm beneath the sample plane, and each LED has an approximate 20 nm spectral bandwidth. Prior work establishes that the impact of non-ideal source coherence is gradual [28], which we examine in further detail in Chapter 5. To quantitatively verify resolution improvement, we turn on each of the 15 × 15 LEDs beneath a U.S. Air Force (USAF) resolution calibration target. A 2X Olympus microscope objective (apochromatic Plan APO, NAo = 0.08) transfers each resulting optical field to a CCD detector (Kodak KAI-29050, 5.5 µm pixels), creating q = 137 low resolution images. Using this 0.08 NA microscope objective (5◦ collection angle) and with illumination NAi = 0.35 illumination NA (20◦ maximum p illumination angle), our FP microscope offers a total complex field resolution gain of n/m = 5 in one dimension, and n/m = 25 for both dimensions together, which follows from Eq. 3.3. Each image spectrum overlaps by ol = 70% in area with each neighboring image spectrum. For reconstruction, we select n = 25 × m and run the above FPM algorithm on our captured data matrix. For computational efficiency, we sequentially process small tiles (500 × 500 pixels per tile) of each image within the data matrix. That is, after capturing all 137 images, we first consider the upper 5002 pixels in each image, and only process this limited image region. Then, we shift our region of interest by 490 pixels to the right, and process a new but slightly overlapping 5002 pixel region of the image set. After processing these limited-FOV image regions, we tile the resulting complex field solutions together, blending at the 10 pixel-wide borders, to form a full-FOV high-resolution image output. The results of our first validation experiment are in Fig. 3.4. Fig. 3.4(c) displays a single image of an Air Force resolution target captured by our microscope under illumination from the central LED using a 2X microscope objective. Here, the small aperture of the 2X objective lens limits the image resolution to approximately 3.2 µm (half-pitch width, the width of one bar in Group 7 Element 3). In Fig. 3.4(d), we show the result of an FPM reconstruction of the same resolution target after acquiring 137 different images, each under unique illumination from a single LED in the array. Resolution improvement over a single image in Fig. 3.4(c) is clear (from 3.2 µm to 0.7 µm half-pitch, now resolving Group 9 Element 3). Across the entire image FOV, FPM has increased the total number of resolvable pixels (i.e., the SP) with respect to a single image from approximately 23 megapixels to 230 megapixels, which equates to 4.6 × 108 individual measurements of amplitude 32 a detector g(x) h(x'-pj) MO lens sample LED Array h(x) pj 10 µm Figure 3.4: The FPM setup (figure adapted from [32]). (a) An LED array sequentially illuminates a sample from different directions, which is then imaged by a 2X microscope objective (MO) lens. (b) The actual FPM setup, showing the LED array and an inset of a single color LED. (c) A single image of a resolution target with this 2X objective offers a wide FOV, but cannot resolve group 8 (approximately 6 µm half-pitch resolution). (d) An FPM-reconstructed image resolves group 9 ( approximately 0.75 µm half-pitch resolution). 33 Low  res  version   One raw data frame 137 uniquely illuminated ~20 megapixel images 0.23 Gigapixel Image (0.78 µm2 pixels, 120 mm2) FP  algorithm   20um   Reconstruction 1 mm! 100 μm! 20um   Figure 3.5: FPM reconstruction of a single layer of stained human blood cells. (middle) Our full image reconstruction, demonstrating the ability to see tens of thousands of individual cells across our FPM setup’s 120 mm2 FOV. Enlarging just one small segment (200×200 pixels, shown in lower right) of this final image (15000×15000 pixels) highlights our ability to clearly resolve each red and white blood cell, whereas the individual cells are hardly resolved in the single raw image. and phase. This enhancement in information throughput can come of great aid during acquisition in many microscope imaging scenarios. To demonstrate this benefit in a biological imaging scenario, we run the same procedure outlined above while imaging a slide of red and white blood cells spread out into a single layer. The results of this experiment are shown in Fig. 3.5. We again transform q = 137 low resolution images into a 0.23 gigapixel image reconstruction. Zooming into a small segment of this final image reconstruction in the lower right of Fig. 3.5, we see that our FPM system clearly resolves each red blood cell, and resolve the fine intercellular features of white blood cells. The raw image segment shown above, which is what a conventional microscope system captures over the same image FOV, cannot clearly resolve any of the cells. Finally, we conclude this section by demonstrating how FPM can also improve the resolution of a pathology slide, enabling whole-slide imaging at approximately 1 µm resolution without any moving parts. In Fig. 3.6, we image a stained pathology sample of human breast tissue exhibiting adenocarcinoma (cancer).The FOV of each image captured by our microscope, based on a 2X objective lens, is nearly as wide as the entire microscope slide. Note that unlike other whole-slide imaging setups for digital pathology, FPM images the entire tissue section at high resolution (comparable to a 20X objective lens) without physically scanning the lens or sample, but instead by simply illuminating different LEDs. Fig. 3.6(b)-(d) zooms into different segments of this image to show our achievable resolution is at the cellular level. We compare what our system can resolve in Fig. 3.6(c1) to the resolution of a conventional microscope using a 20X objective lens in Fig. 3.6(c2), to see each offers comparable performance. However, the FOV of the 20X objective lens setup is approximately 100 times less than our system’s FOV, as indicated by the blue circle marked in Fig. 3.6(a). Beyond 34 20X FOV! Figure 3.6: FPM reconstruction of a pathology slide (figure adapted from [32]). (a) Full color reconstruction (0.23 gigapixels, 120 mm2 FOV). (b)-(d) Different regions of the reconstruction, zoomed-in to highlight the FPM system resolution. Our reconstruction (c1) offers approximately the same resolution as a 20X 0.5 NA objective lens, which we used to capture the same region as shown in (c2). However, a 20X objective lens can only see the FOV denoted by the blue circle in (a) and marked on the microscope slide inset. FPM increases the total number of resolved pixels by approximately 25 in this example, seen by comparing our reconstruction in (c1) to the raw image data in (c3). 35 direct resolution improvement of intensity imagery, our FPM system may also capture quantitative optical phase, extend the microscope depth-of-field, and help estimate and computationally remove intrinsic system aberrations. We detail these extensions in the next chapter. 36 Bibliography [1] R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of the phase from image and diffraction plane pictures,” Optik 35, 237–246 (1972). [2] J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Optics 21, 2758–2769 (1982). [3] S. Marchesini, “A unified evaluation of iterative projection algorithms for phase retrieval,” Rev. Sci. Instrum. 78, 011301 (2007) [4] G. Zheng, R. Horstmeyer and C. Yang, “Wide-field, high-resolution Fourier ptychographic imaging,” Nat. Photonics 7, 739–745 (2013). [5] R. Horstmeyer and C. Yang, “A phase space model of Fourier ptychographic microscopy,” Opt. Express 22, 338–358 (2014). 37 Chapter 4 Fourier ptychography for quantitative phase measurement In Chapter 3, we experimentally verified that Fourier ptychographic microscopy (FPM) improves the resolution of a standard microscope beyond its spatial frequency cutoff limit. In this chapter, we show that FPM also accurately measures the optical phase emerging from thin samples. We demonstrate that this phase information quantitatively measures of sample optical thickness, helps to determine sample scattering properties, and also helps form sharper images through system aberration removal. 4.1 Quantitative phase for sample thickness 4.1.1 Background There are many methods to measure optical phase. Examples include on and off-axis digital holography and interferometry. These interference-based techniques require a highly coherent source and a somewhat complicated optical setup, which is prone to alignment errors, the presence of speckle noise, and must be phase-stable. As we saw in Chapter 3, FPM applies a phase retrieval algorithm to transform optical intensities, from a conventional microscope, into amplitudes and phases. Here, we show these measured phases are accurate and robust. Other strategies exist to computationally determine optical phase without using interference. The Gerchberg-Saxton (GS) algorithm [2], developed for electron microscopy in the 1960s, is one of the earliest strategies for recovering the phase of a specimen from intensity measurements. As outlined in Chapter 3, this iterative procedure alternatively constrains an initial complex estimate to conform to the measured intensity data in the spatial domain, and to obey a known constraint in the Fourier domain. It is closely related to the ER algorithm detailed previously [4]. While proven to weakly converge, stagnation and local minima issues (especially in the presence of noise) limit the applicability of the simple GS/ER algorithms [4, 8]. Gonsalves [4] and Fienup [4, 24] both recognized that applying multiple unique intensity measurement constraints, as opposed to a single 38 intensity constraint, helps prevent stagnation and greatly improves convergence speed. This type of phase diversity procedure has branched into many forms over the past several years, and includes variants based on translational diversity [6], defocus diversity [7], wavelength diversity [8, 9], and sub-aperture piston diversity [10]. Each unique approach here relies on a different technique to vary the experimental system between sequential image capture, before running a variant of the GS/ER iterative algorithms across the entire captured data set. Phase retrieval schemes based on translational-diversity (i.e., moving the sample laterally) are the closest relative to FPM. Phase retrieval with translation diversity is also commonly referred to as ptychography [8, 9, 11]. Ptychographic phase retrieval has found successful use in X-ray [3] and electron microscopy [2]. While setups exist in many flavors [4, 12, 17, 18, 20, 21, 29], the general ptychographic approach consists of three major steps: 1) illuminating a sample with a spatially confined probe beam and capturing an image of its far-field diffraction pattern, 2) mechanically translating the sample to multiple unique spatial locations (i.e., applying translational diversity) while repeating step 1, and 3) using the set of captured images as constraints in an iterative phase retrieval algorithm. While a useful review of details regarding ptychography’s operation are available in [8, 12], it is important to note that the sample or the probe beam must be physically scanned over the desired field-of-view, and the object support is imposed by the shifting probe in the spatial domain. The ability of traditional ptychographic procedures to obtain quantitatively accurate phase measurements has been demonstrated previously with both electron [2] and x-ray [4, 12, 17, 18, 20, 21, 29] illumination, and recently extended to the optical regime [6, 6, 7]. FPM [32] requires no mechanical movement and facilitates microscopic imaging well beyond the cutoff frequency defined by the employed imaging optics, as we saw in the previous chapter. However, it is not directly obvious that FPM, like conventional ptychography, might also capture and postprocess data with its phase retrieval algorithm to construct an accurate solution of the optical phase at the sample plane. Here, we show that although its multiple illumination sources offer a limited spatial coherence and much of the images’ redundant information is utilized to improve spatial resolution, the FPM algorithm also converges to an accurate phase measurement. More details regarding this insight are available in [33]. 4.1.2 Verification with microspheres To verify that FPM measures quantitatively accurate phase, we use a similar experimental setup as outlined in Chapter 3. Our microscope system consists of a conventional microscope with a 15×15 red LED matrix (center wavelength 633 nm, 12 nm spectrum bandwidth, approximately 150 µm size) as the illumination source. This setup is diagrammed in Fig. 6.1, where we label the jth twodimensional illumination plane wavevector as (kxj , kyj ). We have also labeled the two-dimensional Fourier plane coordinates as wavevectors (kx , ky ), which are directly proportional to the spatial 39 j j Image   spectrum   Image   spectrum   j j j j Figure 4.1: FPM setup for quantitative phase measurement (figure adapted from [33]). (a) The full setup, (b) diagram of capture at the aperture plane, (c) close up of plane wave illumination and (d) LED scanning. 40 Raw data, FPM! Phase from FPM! Phase from holography! π Phase (radians)! -π Raw data, FPM! Phase from FPM! Phase from holography! π Phase (radians)! -π Figure 4.2: Quantitative measurement of microsphere and red blood cell thickness (figure adapted from [33]). Here, we also compare FPM phase reconstructions to digital holographic and theoretical data. FPM transforms low-resolution intensity images from a 2X objective (a1) into a high-resolution phase map (a2) of different-sized polystyrene microbeads, as compared with a DH reconstruction (a3) using a 40X objective. (b) A similar image sequence highlights FPM phase imaging capabilities on a human blood smear. (c) Line traces through the microspheres and a red blood cell demonstrate quantitative agreement with expected phase performance. frequencies (x0 , y 0 ) used in previous chapters. In this coordinate system, the center of each shifted spectrum lies at (kxj , kyj ) in the Fourier domain, as shown in Fig. 6.1(b). We insert our sample at the focal plane of the microscope’s 2X, 0.08 NA objective lens. Then, we collect a sequence of q = 225 low-resolution intensity images as the sample is successively illuminated by each of the 225 LEDs in the array. These images are input to the FPM phase retrieval algorithm from Chapter 3, which reconstructs a high-resolution map of the complex field at the sample plane, ψ(x, y). For example, the 500×500 pixel quantitative phase map in Fig. 4.2(a2) is generated from a sequence of 225 different 50×50 pixel cropped low-resolution images, an example of which is displayed in Fig. 4.2(a1). To first verify the accuracy of FPM phase measurement, we image a sample containing polystyrene microspheres in oil (3 µm and 6.5 µm diameter, noil = 1.48, nsphere = 1.6). One raw intensity image is shown in Fig. 4.2(a). It is challenging to resolve each sphere within this raw image data. 41 After running the FPM phase retrieval process, we obtain the phase map in 4.2(a2). Besides clearly resolving each sphere, we may also use this phase map to measure each sphere’s optical thickness. Specifically, the measured phase at each pixel φ and the sample’s thickness t are related via t = k∆φ(nsphere − noil )−1 , where k is the average wavenumber and ∆φ = φ − φ0 is the reconstructed phase minus a constant phase offset. Unwrapped line traces of the optical phase shift induced by two different-sized spheres lead to estimated microsphere thickness curves in Fig. 4.2(c1)-(c2), exhibiting close agreement with theory. The theoretical curves are defined as the length of the vertical chords connecting the top and bottom arcs of a 3 and 6.5 µm circle, respectively (i.e., the diameter of each sphere, which we know a-priori). The root mean-squared error (RMSE) between the experimentally determined and theoretically calculated thickness is 0.25 µm for the 3 µm sphere and 0.33 µm for the 6.5 µm sphere. As a second verification, we also compare our FPM measured phase to that found using a phaseshifting digital holographic (DH) microscope. This microscope consists of a 40X objective lens for clear imaging of each bead. Our DH setup splits a solid-state 532 nm laser into a sample and reference arm (both spatially filtered and collimated). The sample arm passes through the sample. The reference arm passes through an electro-optic phase modulator (Thorlabs EO-PM-NR-C1). Both beams are then recombined and imaged (Prosilica GX 1920, 4.54 µm pixels) via the 40X objective lens (0.65 NA Nikon Plan N) and a tube lens. 4 images are captured with a π/2 phase shift added to the reference between each image. Sample phase is calculated from the 4 images via the phase recovery equation [28]. The resulting phase map is in Fig. 4.2(a3). A RMSE of 0.41 µm and 0.30 µm for the 3 µm and 6.5 µm line traces also offers close agreement between the DH experimental measurements and theory (see Fig. 4.2(c1)-(c2)). 4.1.3 Verification with human blood cells We additionally verify the accuracy of the optical phase measured through a biological sample, in Fig. 4.2(b). Here, we image a human blood smear immersed in oil (noil = 1.51). While a theoretical benchmark is not available, the FPM image in (b2) and the ground-truth DH phase map in (b3) closely match. A trace through the phase of the same red blood cell, noted by the dotted line in (b2) and (b3), is shown in Fig. 4.2(c3). Here, the MSE between the two phase measurements is 0.58 µm, but both closely match the average red blood cell thickness of approximately 2 µm. Sources of error for the FPM setup include the inclusion of slight aberrations by the objective lens, effects of a partially coherent illumination source, and the influence of noise within the iterative reconstruction scheme. The primary source of error in the DH data is speckle noise caused by a coherent illumination source. Artifacts also result from the removal of a global low-frequency background phase fluctuation, caused by an imperfect reference arm. In general, the FPM data tends towards a smoother phase profile, in part a consequence of its LEDs’ partially coherent illumination, 42 Figure 4.3: Using the recovered sample phase, it is direct to simulate image formation under alternative microscope configurations. Here, we show computed phase gradient images in the x direction (a) and y direction (b) from the human blood smear phase map in Fig. 4.2, as well as (c) the phase gradient magnitude to highlight sample boundaries (figure adapted from [33]). which avoids high-frequency spatial fluctuations (i.e., the coherent speckle artifacts). We examine the effects of partially coherent illumination in more detail in the next chapter. Finally, we highlight the benefits of measuring the quantitative phase of blood cells in Fig. 4.3. Here, we use the measured phase to compute phase-gradient images along both the x and y directions. The resulting gradient maps are similar in appearance to that which would be observed through a differential-interference-contrast microscope, where edge-based features are more clearly visible than within a conventional microscope. However, we note that such computational post-processing does not produce new information about the sample beyond what is contained within the recovered amplitude and phase. 4.2 FPM phase measurements for digital pathology While currently of limited use in clinical pathology, quantitative phase measurement may play a important role in future pathology systems that image with digital sensors. Accurately measured phase can reveal both useful sample features and relax digital microscope design constraints. In this section, we briefly outline possible applications of optical phase measurements to the pathologist. A more detailed discussion of this topic may be found in [29]. First, an n-pixel complex image contains up to two times the amount of information than an npixel intensity image. This additional information lies within the image’s complex phase. As is well known with other phase imaging microscope setups (e.g., differential interference and phase contrast microscopy), optical phase may be manipulated to improve the visibility of primarily transparent samples, and provides a direct indication of the optical path length difference (i.e., optical thickness) between adjacent sample regions, which we already demonstrated earlier in this chapter. This 43 Intensity from FPM! Phase from FPM! π Phase (radians)! -π Figure 4.4: Segments of FPM intensity and phase images of a stained tissue (i.e., pathology) sample. As indicated by the red arrow, some cell features are transparent in the intensity image, but clearly visible in the phase image, indicating the utility of phase in revealing otherwise hidden features (adapted from [33]). additional information may help pathologists reach a diagnosis that was otherwise not possible with intensity-only information. An example of this is shown in Fig. 4.4, where the FPM quantitative phase image reveals an otherwise invisible structure contained within a pathology sample (e.g., see the red arrow, or bright highlights in phase where the intensity is low). Such a primarily transparent sample region exhibits a varying thickness that would be quite challenging to detect via direct intensity-only observation. A phase map captured by FPM also contains embedded information about sample scattering properties (e.g., scattering from a histological tissue section). The scattering properties of a sample of tissue are directly tied to its spatially varying refractive index profile. Such refractive index maps can successfully distinguish between healthy and cancerous cells in digitized histology slide images [30]. Cancer can manifest itself within the cell nucleus as a slight index of refraction shift, which is both difficult for the clinical pathologist to detect upon direct observation, and also difficult for a computer algorithm to recognize from an intensity-only image. Furthermore, a phase-derived refractive index map maintains its utility even in the presence of uneven histology sample staining [31], thus suggesting it as a robust computational aid that can fit within pathology’s well-established clinical diagnostic workflow. Next, we briefly demonstrate how to transform an FPM phase image into a spatial scattering coefficient map, which may in turn be used as an aid for diagnosis. Recent work has revealed two close links between the statistics of a sample’s phase map and its scattering parameters [32, 33]. First, the spatial variance across the phase map is linearly related to the sample scattering coefficient, µs . Second, the variance of its phase gradient is related to the reduced scattering coefficient of the sample, µ0s . Denoting the output phase map of FPM as φ(x, y), we may follow the simple set of calculations in [32, 33] to transform φ(x, y) into an estimate of the 44 (a) FPM reconstructed phase! 20  μm   (b) Phase (magnified)! ϕ(x,y) (c) Phase gradient ! π 1  μm   1  μm   ϕ(x,y)| 0 Δ|   -π 1  μm   5.5 π 1  μm   2.5 ϕ(x,y)| Δ|   -π 0 Figure 4.5: (a) An FPM phase map of microspheres, allowing us to verify the accuracy of our procedure to estimate tissue scattering parameters (figure adapted from [29]). (b) Cropped regions of the phase image around two different-sized microspheres. (c) Gradient of the phase in (b). Variances from the regions in (b) and (c) are used to determine the scattering and reduced scattering coefficients of the microspheres following Eq. 4.1 and Eq. 4.2, respectively. scattering and reduced scattering coefficients: µs (x0 , y0 ) = E 1D 2 (φ(x, y) − µφ ) L ∆x,∆y (4.1) µ0s (x0 , y0 ) = E 1 D 2 |∇φ(x, y)| . 2k 2 L ∆x,∆y (4.2) In Eq. 4.1, L is the thickness of the sample, µφ = hφ(x, y)i is the mean of the phase map, and h·i(∆x,∆y) denotes an expectation value over a finite spatial window, (∆x, ∆y). In short, Eq. 4.1 indicates that the scattering coefficient at a particular location (x0 , y0 ) is given by examining the variance of the phase map φ(x, y) within a spatial window (∆x, ∆y) centered on pixel (x0 , y0 ). In Eq. 4.2, ∇ represents the 2D gradient operator. It indicates that a spatial average over a similar finite window of the square of the phase gradient yields the reduced scattering coefficient, µ0s . Both above equations assume the sample under consideration is thin. First, we may verify the accuracy of Eq. 4.1 and Eq. 4.2 through an experiment using a phase map measured by FPM. Using the same experimental setup described above, we first image a slide of microspheres. Given the known size and geometry of the microspheres, we may compute µs and µ0s for direct comparison to our phase-derived expressions, offering a quantitative verification. Second, we will apply the same computations to a histology slide with unknown µs and µ0s . Doing so, we aim to demonstrate how FPM’s acquired quantitative phase map can offer extra information that is otherwise missing from intensity-only imagery. Fig. 4.5 displays an example 500×500 pixel FPM-recovered phase map of a monolayer of two 45 (a) FPM Reconstructed Intensity! Reconstructed Phase! Derived scattering coefficient images! radians! π 40  μm   40  μm   40  μm   1  mm   -π (d) Log(µs' (x,y)) (c) µs (x,y) (b) ϕ(x,y) 0 A.U.! 0.4 0 A.U.! 0.25 Figure 4.6: Recovering a spatial map of tissue scattering parameters from a histology slide (figure adapted from [29]). (a) Full color FPM gigapixel image of a histology slide. (b) The recovered phase from a 0.2 mm2 region at the sample center. (c) The scattering coefficient µs for this area of tissue, computed using Eq. 4.1. (d) The reduced scattering coefficient µ0s for the same area of tissue from Eq. 4.2, here shown on log scale. sizes of microspheres (r = 3 µm and 6.5 µm, n = 1.6) immersed in oil (n = 1.48). This limited FOV image is one cropped area of a larger image. Selecting one 6.5 µm microsphere of interest, as shown in (b) top, we first apply Eq. 4.1 to find that µs =0.696 µm−1 and Eq. 4.2 to find that µ0s =0.0239 µm−1 . Here, we have used λ = 632 nm and a sample thickness L = 6.5 µm. These two values may be combined to find the anisotropy factor g = 1 − (µ0s )µs = 0.966. We can compare these scattering coefficients to predictions from Mie theory code [34] for a sparse set of spheres of similar size and refractive index, which yields µs =0.658 µm−1 , µ0s =0.0241 µm−1 , and g = 0.963. These predicted values match our experimental measurements closely, thus confirming the accuracy of our quantitative phase map as well as verifying the validity of Eq. 4.1 and Eq. 4.2. The same set of calculations applied to one selected 3 µm microsphere (shown in Fig. 4.5, bottom row) offers similarly accurate results. Given a verifiably accurate process of measuring the scattering and reduced scattering coefficients from the FPM phase map, we now test a biological sample. Fig. 4.6(a) includes an example FPM histology slide image, shown here in color after combining three separate reconstructions under illumination from a set of red, green, and blue LEDs. Fig. 4.6(b) displays the phase map of a 1200×1200 pixel area of interest, from which we calculate the phase and phase gradient variance for each image pixel following Eq. 4.1 and Eq. 4.2. As with the microspheres, we again assume (∆x, ∆y) is a 22-pixel window and the sample thickness is 5 µm (an estimation). The result of these calculations are the spatial maps of µs and µ0s shown in Fig. 4.6(c)-(d), where regions of homogenous scattering and specific points of inhomogeneity are clearly observable. A pathologist may find the larger µ0s in the upper left corner of the image segment a helpful indication of a different tissue region. This type of sectioning is otherwise challenging to achieve with intensity-only imagery. Alternatively, the large peaks in the µ0s map might indicate cells that exhibit unique (i.e., highly 46 (a) Z = -150 microns! (b) Z = 0 microns! (c) Z = +150 microns! (d) Trace through Group 9 El. 3! Figure 4.7: Example of digital refocusing with FPM data (figure adapted from [32]). (a) Raw image recorded from the center LED at -150 µm defocus, showing clear defocus effects at top, and FPM reconstruction after digital refocusing at bottom. (b) Same without defocus, and (c) same with defocus +150 µm in the opposite direction. FPM reconstruction, after digital refocusing, maintains resolution of Group 9, Element 3, thus offering a depth of field of approximately 300 µm, as compared to the 80 µm depth of field of the 2X objective lens during standard operation. scattering) structural properties. We hope that these initial findings encourage the future integration of FPM into the diagnosis pipeline of digital pathology. 4.3 Improving reconstruction quality with FPM phase measurements Finally, we briefly note here two applications of FPM-derived phase that additionally improve the resolution and fidelity of the reconstruction. These insights may also help relax current microscope design constraints (e.g., increase the working distance, extend the depth of field, and/or reduce the objective lens size and complexity). However, since they do not form the main focus of this thesis, we only mention them here briefly, and point the interested reader to the relevant sections of Refs. [32] and [35] for more details. 4.3.1 Digital refocusing for enhanced depth-of-field The phase acquired by FPM can be used to digitally refocus images into sharp focus. We achieve digital refocusing by propagating the complex field solution towards or away from the image plane via a propagation simulator. Without an accurate measurement of sample phase at the image plane, this type of digital propagation is not possible. Since digital refocusing can be applied to FPM data at any point after image acquisition, we anticipate that it may help improve the efficiency of a pathologist’s diagnosis pipeline. For example, a slide that was originally slightly tilted or a histology sample that was not perfectly flattened between 47 Aberrations, red channel! Aberrations, green channel! Aberrations, blue channel! -π π Region c! Region b! Region a! Intensity! Figure 4.8: Removing aberrations from FPM (figure adapted from [35]). By slightly modifying the FPM iterative phase retrieval procedure, it is possible to simultaneously update an estimate of the complex sample as well as an estimate of the objective lens aberrations. This simultaneous update effectively separates the effect of lens aberrations from the recovered sample image. (left) An example recovered complex sample with aberrations removed. (a)-(c) Magnified views of labeled regions, where column (2) shows the associated phase aberrations removed from the red color channel, column (3) from the green color channel, and (4) from the blue color channel. It is clear this technique can remove aberrations that are both a function of wavelength and spatial location. its slide and cover glass will not have to be re-imaged. Instead, a technician can computationally correct minor slide positioning errors at any later date using the acquired sample phase information. Whats more, by digitally refocusing any sample region not in sharp focus to its correct focal plane, FPM can significantly extend a microscope’s effective depth-of-field, e.g., from 80 µm to 300 µm for the setup demonstrated in [32, 33]. An example of using digital refocusing to bring a resolution target axially offset from the focal plane by +/− 150 µm is in Fig. 4.7. One may either use an a-priori estimate of the required defocus distance, as in the example in Fig. 4.7, or rely upon an image sharpness metric, to establish an unknown defocus distance. Additional details regarding FPM digital refocusing are available in [32]. 4.3.2 Simultaneously removal of microscope aberrations The FPM phase retrieval algorithm outlined in Chapter 3 is only a starting point for more involved recovery procedures. As previously explored within the realm of conventional ptychography, more advanced computational methods can simultaneously estimate any errors shared between captured images of a sequentially illuminated sample. Typically, these fixed errors arise from a non-uniform aperture function. For FPM, this corresponds to the optical aberration map of the microscope objective lens (i.e., the phase deviations across the back focal plane). As demonstrated in [35], it is possible to simultaneously estimate and remove this optical aberration map during FPM reconstruction. Iterative removal of the aberration map significantly sharpens 48 the spatial resolution of the final FPM reconstruction, which is especially beneficial in off-axis regions of the image plane (i.e., in the corners of the image). Because both the sample and aberration estimates are complex (i.e., include optical phase), this removal is not ill-conditioned (e.g., like point spread function deconvolution from intensity-only images). What’s more, it requires no initial calibration data. An example of FPM aberration removal is included in Fig. 4.8. Further details regarding this correction procedure may be found in [35]. 49 Bibliography [1] R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of the phase from image and diffraction plane pictures,” Optik 35, 237 (1972). [2] J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21 (1982). [3] J. R. Fienup and C. C. Wackerman, “Phase-retrieval stagnation problems and solutions,” J. Opt. Soc. Am. A 3, 1897–1907 (1986). [4] R.A. Gonsalves, “Phase retrieval and diversity in adaptive optics,” Opt. Eng. 21, (1982). [5] R. G. Paxman, T. J. Schulz and J.R. Fienup, “Joint estimation of object and aberrations by using phase diversity,” J. Opt. Soc. Am. A 9 (1992). [6] M. Guizar-Sicairos and J.R. Fienup, “Phase retrieval with transverse translation diversity: a nonlinear optimization approach,” Opt. Express 16 (2008). [7] R. G. Paxman and J. R. Fienup, “Optical misalignment sensing and image reconstruction using phase diversity,” J. Opt. Soc. Am. A 5, 914–923 (1988). [8] P. Bao, F. Zhang, G. Pedrini and W. Osten, “Phase retrieval using multiple wavelengths,” Opt. Lett. 33(4), (2008). [9] L. Waller, S. S. Kou, C. J. R. Sheppard, and G. Barbastathis, “Phase from chromatic aberrations,” Opt. Express 18(22), (2010) [10] M. R. Bolcar and J. R. Fienup, “Method of phase diversity in multi-aperture systems utilizing individual sub-aperture control,” Proc. SPIE 5896 (2005). [11] W. Hoppe, “Diffraction in inhomogeneous primary wave fields. 1. Principle of phase determination from electron diffraction interference,” Acta Crystallogr. A25, 495–501 (1969). [12] J.M. Rodenburg and R.H.T. Bates, “The theory of super-resolution electron microscopy via Wigner-distribution deconvolution,” Phil. Trans. R. Soc. Lond. A 339, 521–553 (1992). [13] H. N. Chapman, “Phase retrieval x-ray microscopy by Wigner distribution deconvolution,” Ultramicroscopy 66, 153 (1996). 50 [14] J. M. Rodenburg et al., “Hard-X-ray lensless imaging of extended objects,” Phys. Rev. Lett 98, 034801 (2007). [15] F. Hue et al., “Wavefront phase retrieval in transmission electron microscopy via ptychography,” Phys. Rev. B 121415(R) (2010). [16] K. Nugent, “Coherent methods in the x-ray sciences,” Advances in Physics 59(1), 1–99 (2010). [17] I. Johnson et al., “Coherent diffractive imaging using phase front modifications,” Phys. Rev. Lett. 100, 155503 (2008). [18] C.T . Putkunz et al., Phase-diverse coherent diffractive imaging: high sensitivity with low dose, Phys. Rev. Lett. 106, 013903 (2011). [19] P. Thibault et al., “High-resolution scanning x-ray diffraction microscopy,” Science 321, 379–382 (2008). [20] K. Giewekemeyer et al., “Quantitative biological imaging by ptychographic x-ray diffraction microscopy,” Proc. Nat. Acad. Sci. 107 (2), 529–534 (2010). [21] Y. Takahashi et al., “Towards high-resolution ptychographic x-ray diffraction microscopy,” Phys. Rev. B 83, 214109 (2011). [22] P. Thibault and A. Menzel, “Reconstructing state mixtures from diffraction measurements,” Nature 494, 7435 (2013). [23] H. M. L. Faulkner and J. M. Rodenburg, “Movable aperture lensless transmission microscopy: A novel phase retrieval algorithm,” Phys. Rev. Lett. 93, 023903 (2004). [24] A. M. Maiden, J. M. Rodenburg and M. J. Humphry, “Optical ptychography: a practical implementation with useful resolution,” Opt. Lett. 35(15), 2585–2587 (2010). [25] A. M. Maiden, M. J. Humphry, F. Zhang and J. M. Rodenburg, “Superresolution imaging via ptychography,” J. Opt. Soc. Am. A 28(4), 604“612 (2011). [26] G. Zheng, R. Horstmeyer and C. Yang, “Wide-field, high-resolution Fourier ptychographic microscopy,” Nat. Photonics 7, 739–745 (2013). [27] X. Ou, R. Horstmeyer, G. Zheng and C. Yang, “Quantitative phase imaging via Fourier ptychographic microscopy,” Opt. Lett. 38, 4845–4848 (2013). [28] I. Yamaguchi and T. Zhang, “Phase-shifting digital holography,” Opt. Lett. 22(16), 1268–1270 (1997). 51 [29] R. Horstmeyer, X. Ou, G. Zheng, P. Willems and C. Yang, “Digital pathology with Fourier ptychography,” Computerized Medical Imaging and Graphics 42, 38–43 (2015). [30] Z. Wang, H. Ding and G. Popescu, “Scattering phase theorem,” Optics Letters 36(7), 1215–1217 (2011). [31] M. Xu, “Scattering phase theorem: anomalous diffraction by forward-peaked scattering media,” Optics Express 19(22), 21643–21651 (2011). [32] P. Wang et al., “Nanoscale nuclear architecture for cancer diagnosis beyond pathology via spatial-domain low-coherence quantitative phase microscopy,” Journal of Biomedical Optics15, 066028 (2010). [33] S. Uttam, R. .K. Bista, D. J. Hartman, R. E. Brand and Y. Liu, “Correction of stain variations in nuclear refractive index of clinical histology specimens,” Journal of Biomedical Optics 16, 116103 (2011). [34] S. Prahl, Interactive Mie Calculator. Oregon Medical Laser Center (2012), accessed Jan. 14, 2014. http://omlc.ogi.edu/calc/mie_calc.html [35] X. Ou, G. Zheng and C. Yang, “Embedded pupil function recovery for Fourier ptychographic microscopy,” Optics Express 22, 4960–72 (2014). 52 Chapter 5 Modeling ptychography in phase space In ptychographic imaging, also commonly referred to as scanning diffraction microscopy or more broadly as coherent diffractive imaging, a sample is shifted across a narrow illumination beam and a series of diffraction intensity patterns are recorded. The acquired image data is then computationally processed into an improved-resolution estimate of the sample’s amplitude and phase transmittance. The unique procedure of ptychography has recently lead to the generation of many impressive X-ray and electron microscope images that defy the conventional resolution limitations of their detectors and focusing elements [1–4,38]. Such resolution enhancement has also spread to optical imaging [6,7]. This multi-image procedure closely matches our description of FP thus far. In this chapter, we first develop a mathematical model for “conventional” ptychography (CP) based in phase space. Second, we use the same model to explain how Fourier ptychography (FP) operates. We show how these two otherwise unique optical setups – one capturing the diffracted light from a moving sample, and the other capturing microscope images of a fixed sample illuminated by an array of sources – create nearly identical datasets. Finally, we derive how the partial coherence of the optical sources in each experimental setup impact their measurements, and test this derivation with a simple experiment. By mathematically connecting FP with the well-established method CP, we hope to encourage a cross-pollination of ideas between the two experimental disciplines. As we will see in the next three chapters, many of our new insights concerning FP may port directly over to the discipline of CP. While this has yet to be attempted, we believe that both our convex algorithm and unique tomographic reconstruction technique (the subjects of the next two chapters) may be of significant benefit to the X-ray and electron microscope communities. Conventional Ptychography Setup! Source! 53 Aperture a(x )! Sample ψ(x'-p )! Detector! shift position p" Image g(x)! w" δx   l L(x') ! f" A(x) ! d" S(x') ! D(x) ! Figure 5.1: The optical setup for conventional ptychography (CP). A sample ψ (in green) is shifted through many positions p as the intensity of the probe light it diffracts is recorded at a far-field detector. In a typical visible light setup, the lens at planeA(x) is a multi-element system containing the aperture stop function, a(x), at some intermediate plane, as diagrammed (figure adapted from [28]). 5.1 The conventional ptychography (CP) setup Our first steps toward a common mathematical framework are to outline the standard elements of a CP setup, model how light passes through it, and then convert our findings into a suitable phase space representation. The basic setup, notations, and derivations used here closely follow those previously employed in [8, 9]. Unlike this prior work, our final expression demonstrates a unique convolution relationship that helps directly connect the parameters in CP with FP. Reciprocal space coordinates will be designated with a prime, and reciprocal space functions will include a tilde (e.g., the Fourier transform of a(x) is ã(x0 )). Note that in practice, both x and x0 will have units of meters, since they represent the spatial axis of an imaging system’s two Fourier conjugate planes. A schematic diagram of a scanning CP setup containing two sets of such planes is in Fig. 5.1. While deviations exist, most recent ptychographic experiments generally follow the optical outline in Fig. 5.1. The following analysis considers a two-dimensional imaging geometry, for simplicity. Extension to three dimensions is direct. A standard CP setup first focuses light from an illumination plane L(x0 ) onto a shifting sample and records a series of far-field diffraction patterns. We assume L(x0 ) contains an ideal point light source that produces a quasi-monochromatic plane wave (wavelength λ) propagating parallel to the optical axis at a large distance `. The case of a non-ideal point source will be considered later, when partial coherence is added. At distance ` is an aperture plane A(x) containing a lens of focal length f . Immediately past this plane, the optical field may be described across all space simply as a(x), the aperture transmission function. This incident plane wave, confined to a(x), is focused by the lens to a small area at the sample plane, S(x0 ). Note here the sample plane in conventional ptychography is considered as lying in 54 reciprocal space coordinates, which will allow us to directly connect it to FP shortly. Under the Fresnel approximation, the shape of the focal spot before hitting the sample is proportional to the scaled Fourier transform of the field at aperture transmission function [10]: s+ (x0 ) = exp  jk 02 2f x jλf  Z   jk a(x)exp − x0 · x dx ≈ Fxx0 [a(x)] = ã(x0 ), f (5.1) where Fxx0 is the Fourier transform operator from x to x0 , s+ (x0 ) is the field directly before the sample, and the approximation assumes the phase pre-factor is unity. All integrals are assumed to extend from negative to positive infinity. The above expression ignores a constant coordinate scaling factor: ã(x0 ) should actually be written as ã(x0 /λf ). For clarity, we will generally neglect constant scaling factors. Details of scaling effects may be found in Appendix A of [28]. ã(x0 ) typically takes the form of a sinc function, as shown in Fig. 5.2, but may be arbitrarily shaped in general. Independent of its specific distribution, the confined beam ã(x0 ) then interacts with a shifted sample ψ to produce an exiting optical field, s(x0 ). We assume the effect of sample thickness upon diffraction is negligible, allowing us to define the optical field s(x0 ) directly past the sample as a multiplication of ã(x0 ) and the sample transmission function ψ: s(x0 ) = ã(x0 )ψ(x0 − p). (5.2) Here, p is the sample’s shift distance perpendicular to the optical axis. The thin object approxi2 mation holds if the maximum sample thickness t obeys t << 4δres /πλ, where δres is the sampling resolution [12]. s(x0 ) then propagates a large distance d to far-field detector plane D(x), where (as a first approximation) the intensity of the Fourier transform of s(x0 ) is measured: g(p, x) = |Fx0 x [ã(x0 )ψ(x0 − p)] |2 . (5.3) Here, g(p, x) is a two dimensional function of probe shift distance (p) and space (x), and comprises our data matrix. In experiment, g(p, x) is filled up, column-by-column, with discretized diffraction images captured at the detector for many shift distances p (see example in Fig. 5.2). For twodimensional images, g(p, x) is a four-dimensional function. Note that throughout the rest of this chapter, we write the data matrix with images in its columns, instead of in its rows (i.e., as g(x, p), which is the style used in previous chapters). We choose to follow this convention to keep with the correct notation of phase space functions using primal variables along the horizontal axis, as will become clear shortly. 55 PTY! Plane S(x')! Sample ψ(x'-p )! Probe a(x ')! pth diffraction image ! F  [ψ(x'-p) Ÿ a(x')]   Intensity (AU)! optical field! Shift! Ÿ   .1! -.2! Wa (x', u)! x'-p (mm)! -.2! .2! u (m-1) ! u (m-1) ! x (mm)! .2! g (p, x)! Wψ (x', u)! 1! x (m)! x' (mm)! -.1! WDF model! Plane D(x)! Data  Matrix     1! =   *  2D   0! -.2! x' (m)! x' (m)! p (m) ! Figure 5.2: Conventional ptychography (CP) data acquisition. A chirped amplitude grating (400 µm wide, 4µm minimum pitch) serves as our sample ψ(x0 ). It is shifted and illuminated by a probe function ã(x0 ), which is here a sinc function from a rectangular-shaped focusing element. At detector plane D, the diffracted light’s intensity is recorded. (Bottom) Corresponding probe and sample Wigner functions, whose two-dimensional convolution creates the CP data matrix, g(p, x) (figure adapted from [28]). 5.1.1 Phase space representation of CP The structure of CP’s data matrix, g(p, x), reveals information about the spatial frequency content of the sample ψ along the x dimension, thanks to the Fourier transform in Eq. (5.3). Likewise, since the probe only hits a narrow segment of the sample ψ at a given scan position, the spatial structure of ψ is also partially resolved along the scanning dimension p. This joint preservation of both the spatial and spatial frequency sample information within g is a property held by any optical phase space function [13], as first explored early on in [8]. A clear connection between g(p, x) and optical phase space is found by applying a few mathematical transformations to Eq. (5.3). First, expanding it into integral form produces, ZZ g(p, x) = ã(x01 )ã∗ (x02 )ψ(x01 − p)ψ ∗ (x02 − p)exp [−jkx · (x01 − x02 )] dx01 dx02 , (5.4) where the double integral over new spatial variables (x01 , x02 ) results from measurement of intensity at the detector, and ∗ denotes complex conjugate. From here, straightforward manipulations produce an expression for the data matrix g as a convolution of two functions: ZZ g(p, x) = Wψ (x0 − p, u)Wã (x0 , x − u) dx0 du, (5.5) 56 where constant pre-integral multipliers are neglected for clarity. The function W applied to ψ takes the form, Wψ (x0 , u) = Z ∞  y ∗  0 y ψ x0 + ψ x − exp(−jkyu) dy 2 2 −∞ (5.6) and is known as the Wigner distribution function (WDF) of ψ. Eq. (5.5) describes the set of diffraction intensity images in CP as a convolution of two functions. Each function is solely related to the shape of the sample and the probe beam, respectively (i.e., the WDF separates the sample transmission function and probe beam into a linear expression). This is graphically depicted in Fig. 5.2. Note that while not explicitly included in this paper, the interested reader is invited to use the derivation steps in Appendix B of [28] to help create Eq. (5.5) from Eq. (5.4). The WDF is a well-studied phase space distribution that is often used to analyze optical imaging setups [13–15]. Like the Fourier transform, it transfers a function of one “primal” variable, here x0 , into a new space. Unlike the Fourier transform, which offers a one-to-one mapping between the primal variable x0 and its conjugate u (here a mapping between the space and spatial frequency coordinates at the sample plane S(x0 )), this new space is two-dimensional. The WDF is a joint function of both the primal spatial variable x0 and the conjugate spatial frequency variable u. Although defined in a higher-dimensional space, Wψ maintains a one-to-one relationship with the complex function ψ (apart from a constant phase shift). While not always exact, it is convenient to connect the value of W (x00 , u0 ) to the amount of optical power at point x00 propagating in direction u0 . However, while the WDF is real-valued it is not necessarily non-negative, which requires this interpretation to be taken loosely. The goal of ptychography’s many post-processing algorithms is to recover the complex sample function ψ, which has a one-to-one relationship with Wψ , from its recorded dataset g. This goal is computationally related to deconvolving the effect of the aperture a, described by Wã , from g(p, x) in Eq. (5.5). After deconvolution, one has an estimate for Wψ , which can be converted into an estimate for ψ. Deconvolution is often indirectly achieved through a phase retrieval algorithm [6], very similar to the post-processing of FP, which we now summarize using the same mathematical methods. 5.2 Mathematical representation of Fourier ptychography (FP) FP also acquires a sequence of images that are compiled into a data matrix (here labeled gF ) but does so using the unique optical setup in Fig. 5.3. Two primary experimental differences set FP apart from the CP setup outlined above: an array of n LEDs now occupy the illumination plane L(x0 ), and the locations of the sample and aperture planes are effectively switched. Instead of recording the diffraction pattern from a small illuminated sample region, FP images the entire sample under illumination from different directions. Fourier Ptychography Setup! LED Array! 57 Sample ψ(x )! Aperture a(x' )! Detector! shift illumination angle hi" Image gF (x )! ! w" δx   l L(x') ! do" S(x) ! di" A(x') ! D(x) ! Figure 5.3: Fourier ptychography (FP) optical setup. An LED array replaces the single illumination source for CP in Fig. 5.1, and planes S(x) and A(x0 ) have switched places along the optical axis. Each LED sequentially illuminates the sample, now with coordinate x, from a different angle (figure adapted from [28]). Again, we begin by assuming each LED in the array at illumination plane L(x0 ) emits a quasimonochromatic and spatially coherent field at wavelength λ. Each LED sequentially illuminates the entire sample plane S(x), a distance ` away with an angled plane wave. Next, the illuminated sample is imaged by a lens of focal length f , which is located at the aperture plane A(x0 ). At detector plane D(x), a pixel array samples the image intensity at spacing δx = λw/2f (to avoid aliasing issues). From Eq. (5.1), we note that the optical field at A(x0 ) is proportional to the Fourier transform of the field both at the sample plane S(x) and the image plane D(x), a feature that distinguishes FP from CP and lends to its name. Again applying the thin object approximation, the optical field s(x) directly past the sample plane may be written as a multiplication between the incident field and the sample transmission function ψ as, s(x) = ψ(x)ejkxpj . (5.7) Here, pj represents the sine of the plane wave generated by the j th LED, located a distance hj away q from the optical axis: pj = hj / h2j + `2 , with ` the distance between the LED array and the sample. The optical field s(x) continues to propagate to aperture plane A(x0 ), mathematically represented through the scaled Fourier transform in Eq. (5.1). The field is attenuated at A(x0 ) by the aperture transmittance function a(x0 ) (i.e., the shape of the MO pupil plane), creating the optical field, Fxx0 [s(x)] a(x0 ) = s̃(x0 )a(x0 ) = ψ̃(x0 − p)a(x0 ). (5.8) Again, we’ve neglected coordinate scaling factors for clarity (see Appendix A). Finally, this attenuated field propagates to image plane D(x), represented through a scaled Fourier transform. At 58 optical field! FPM! Plane S(x)! Plane A(x')! Plane D(x )! Sample ψ(x )! Aperture a(x' )! pth low-res. image ! amplitude! Intensity (AU)! .2! u (m-1) ! WDF model! Sample Wψ (u, x')! -35! x' (mm)! 35! Aperture Wa (u, x')! x (mm)! .2! gF (p, x)! Data  Matrix     =   *   x' (m)! -.2! x (m)! x (mm)! u (m-1) ! -.2! x' (m)! p (m) ! Figure 5.4: FP data acquisition diagram. (Top) The same grating sample ψ(x0 ) used in Fig. 5.2 is sequentially illuminated by tilted plane waves, adding a different linear phase ∝ p to each image (tilted green line). At plane A(x0 ), the aperture a(x0 ) limits the extent of the field before the sample is imaged to detector plane D(x) at low resolution. (Bottom) Corresponding WDFs and their convolution, representing the FP data matrix. Note each WDF here is rotated by 90 degrees with respect to its typical orientation as well as the orientation in Fig. 5.2, due to the need to swap x0 and u to correctly model FP. Color maps here follow those included in Fig. 5.2 (figure adapted from [28]). 59 D(x), the digital pixel array detects the field’s intensity gF (p, x): h i gF (p, x) = |Fx0 x ψ̃(x0 − p)a(x0 ) |2 . (5.9) As in Eq. (5.3), the p-variable in gF (p, x) also represents a shift parameter. Here, instead of describing the sample’s lateral displacement, it connects each image to the j th illumination LED. Specifically, it informs us how much the sample spectrum has shifted across our fixed imaging aperture. Each column of FP’s data matrix in Eq. (5.9) contains a 1D image captured under a unique illumination direction from one of the j ∈ {1, . . . , n} LEDs in the array, and thus a unique spectrum position. The simulated FP data matrix in Fig. 5.4 is visually quite similar to CP’s, shown in Fig. 5.2. The only mathematical difference between the two data matrices, expressed compactly in Eq. (5.3) and Eq. (5.9), is whether the aperture function a or the sample function ψ is Fourier-transformed. To more explicitly connect CP and FP, we can expand Eq. (5.9) into, ZZ gF (p, x) = ψ̃(x01 − p)ψ̃ ∗ (x02 − p)a(x01 )a∗ (x02 )exp (−jkx · (x01 − x02 )) dx01 dx02 . (5.10) A straightforward derivation (detailed in Appendix B of [28]) leads to the more compact representation, ZZ gF (p, x) = Wψ (−u − p, x0 )Wã (u, x0 − x) du dx0 . (5.11) Here, the functions Wψ and Wã are again the WDF’s of the sample and the Fourier transform of the aperture, respectively, as included in our phase space model for CP in Eq. (5.5). However, the u and x0 variables within each Wigner function have traded places. This leads to a 90 degree rotation of each phase space function. Directly comparing Eq. (5.5) and Eq. (5.11) yields the following simple relationship between the data collected by CP and FP: gF (p, x) = g(x, −p). (5.12) Fourier ptychographic microscopy’s data matrix is simply a rotated version of the data recorded by a conventional ptychography setup (i.e, trading the x and p variables is equivalent to a 90◦ matrix rotation). Since the data collected by each procedure is related through an isomorphic transformation, various post-processing algorithms developed specifically for CP will, in theory, work equally well with FP, and vice-versa. Besides allowing CP and FP to exchange data-processing techniques, Eq. (5.12) also suggests that alternative optical setups can capture the data matrix g under different linear transformations (e.g., a matrix rotation that is not 90◦ , or another isomorphic transform besides rotation). These alternative setups will most likely offer application-specific advantages. For example, one could imagine jointly shifting the sample across a limited range and using a small 60 number of illumination sources to increase collection efficiency of a setup that only tolerates minimal movement. Of course, many practical considerations can influence a selection between different optical setups that produce mathematically equivalent data. In the next section, we will examine how such experimental requirements and practical sampling conditions manifest themselves within our mathematical framework, before deriving a more detailed model including the effects of partially coherent light. 5.3 Visualizing connections between both ptychographic domains The phase space model in the above section offers an excellent visualization of the close link between the data collected by CP and FP. This connection will prove highly advantageous as research surrounding each setup continues to progress, because both systems may draw benefits from new computational and experimental developments in either area. However, it is not correct to assume the exact linear relationship in Eq. (5.12) implies that CP and FP are experimentally identical – a number of system-specific factors may influence each data matrix uniquely. The first goal of this section is to use our phase space model to visualize how experimental factors impact data collection, as Fig. 5.5 outlines. This type of visualization can help inform future decisions about which optical applications each setup is best suited for. At the same time, ensuring the two setups produce data exactly following Eq. (5.12)’s rotation relationship is not particularly challenging. The second goal of the following discussion is to identify a set of carefully chosen setup parameters that lead to such an exact relationship, which we will use in a comparison in a subsequent section. Most experimental aspects of CP and FP fit nicely into one of four categories describing a particular data matrix property: 1. Scaling along the optical axis: For both ptychographic procedures, distances between the optical source, sample, detector, and the lens focal length will lead to constant scaling variations along x and p in their respective data matrices. Details of these scaling relationships are presented in Appendix A of [28]. 2. Sampling along x: The digital detector’s sampling conditions for CP and FP both manifest themselves along their corresponding data matrices’ x axis (Fig. 5.5, green text). For CP, the detector width must match the aperture’s maximum transmitted spatial frequency. This width defines the resolution limit of a final reconstructed image. The detector size and distance together define a geometric NA, which much match the detector pixel size to avoid aliasing [9]. For FP, sampling along the x axis follows a typical imaging setup: the detector width is paired to the imaging lens FOV, and the detector pixel size matches the imaging optics’ point-spread function (PSF) width to 61 CP! FP! Data matrix! Image  i   x! g (p,x) Image  i   Detector width!   x! g (p,x) F Image FOV!   Pixel/ PSF size! Pixel size! p! p! LED pitch! Scan step! Max. LED angle! Scan range! (Max. scan range/ accepted kx)! Finite probe (overlap)! u! W (x', u)! a Blur Kernel! u! W (x', u)! ã Max. spatial frequency! x'! Image PSF! x'! Figure 5.5: The experimental factors influencing CP and FP data matrices. (top) Geometrical factors define the data matrix scaling and sampling, while (bottom) parameters specific to the focusing/imaging lens define the function responsible for data matrix blurring for both setups (figure adapted from [28]). avoid aliasing. 3. Scanning along p: Sampling along the data matrix p-dimension is tied to the operation of the illumination in each setup (Fig. 5.5, blue text). In CP, the probe beam’s total scanning distance sets the maximum extent along p, which also defines the final reconstructed image’s FOV. In FP, however, the maximum extent along p is set by the maximum LED-sample illumination angle. This in turn defines the final reconstructed image’s maximum resolution, as opposed to FOV. This outstanding feature of FP allows for the extension of a lens’s typical resolution cutoff by simply illuminating the sample from large off-axis angles. The sampling rate along the p-dimension of both data matrices is set by the number of captured images. Our above model assumes the WDF is ideally discretized, requiring the number of detector pixels along x to match the number of collected images along p. In practice, accurate high-resolution sample reconstruction does not require full population of g(p, x) or gF (p, x) along p over a fixed scanning range [50]. Under-sampling remains an effective strategy because the WDF is a redundant 2D representation of a complex 1D signal. Phase retrieval algorithms exploit this redundancy to arrive at accurate sample and probe reconstructions, as shown in [18]. 4. Data matrix blur kernel: The finite probe width in CP causes blurring between images, and the finite extent of its aperture will typically define the maximum spatial frequency cutoff for 62 each image. These limiting effects manifest themselves, respectfully, along the x0 and u dimensions of CP’s aperture WDF, Wã (x0 , u), shown in the bottom of Fig. 5.5. Convolution with Wã (x0 , u) in Eq. (5.5) describes how sample information is blurred during the detection process. Since it is zero beyond a certain cutoff value along u, Wã (x0 , u) removes from the data matrix any sample information above this associated spatial frequency range. FP’s rotated blur kernel Wã (u, x0 ) is defined by its imaging aperture. It also blurs and cutoffs sample information from the data matrix in a similar manner as CP’s blur kernel. The simulations presented later in this chapter use a fixed set of example setup parameters, to ensure that the data matrices from the example CP and FP setups only vary by a rotation. For CP, we assume a lens (diameter w = 37.5 mm, focal length f = 105 mm) creates a sinc of estimated width 18 µm (peak-to-zero) at the sample from an LED located ` = 300 mm away. The sample plane contains a grating with 4 µm minimum feature size that is shifted in 4 µm steps. In Fig. 5.2 and Fig. 5.4, we simulate a grating that is 0.4 mm wide. We assume a 4 mm-wide detector, containing 4 µm pixels and with full factor, captures its diffraction pattern, which approximately requires d = 30 mm assuming free space propagation. For FP, we assume a similar lens (diameter w = 37.5 mm, do = 300 mm and di = 105 mm) images the sample onto an identical detector. FP’s LED array is fixed at a distance l = 100 mm and illuminates the same sample. The array extends across a total distance h = 24 mm perpendicular to the optical axis, yielding a 240 µm pitch for Fig. 5.2 and Fig. 5.4. One important parameter still missing from the above analysis is the light’s coherence state, connected to the active area of each optical source. We will now extend our phase space model to account for this critical effect. 5.4 A complete statistical model with partially coherent light In practice, the illumination sources used by each form of ptychography exhibit a limited spatial and temporal coherence. The rarity of ideally coherent electron and X-ray sources has led to the theoretical and experimental examination of coherence effects in CP setups [19, 29]. Here, we reexamine the impact of partial coherence in CP, and introduce for the first time a thorough analysis of partial coherence for FP. We use our phase space model to show that in either optical setup, partially coherent illumination does not limit the ability to recover an exact sample amplitude and phase estimate. We conclude that while partial coherence impacts CP and FP performance differently, it remains a mathematical separable expression that can be removed by computational post-processing. This finding motivates future investigations to examine the potential for a variety of limited-coherence, high-throughput sources to improve the acquisition rates and noise performance for both CP and FP. 63 5.4.1 Partially coherent source description To accurately model experimentally realistic optical sources, we must introduce a statistical measure of spatial coherence into our phase space descriptions of CP in Eq. (5.5) and FP in Eq. (5.11). We achieve this by treating the field emitted by the optical source, U (x0 , t), as a temporally stationary stochastic process and examining its correlation across space and time: hU (x01 , t1 )U ∗ (x02 , t2 )i = Γ̃(x01 , x02 , τ ). Here, Γ̃ is the light’s mutual coherence, τ = t2 − t1 is a constant time difference, and the expectation value is performed over time. Furthermore, we model the source using the conjugate coordinate x0 , since this is the plane where we expect to find the illumination source. From the Weiner-Khinchine theorem, the cross-spectral density (CSD) of this stochastic process is R defined through a Fourier transform of the mutual coherence, Γ(x01 , x02 , ω) = Γ̃(x01 , x02 , τ )e−jωτ dτ . The spectral density C(x0 , ω) = Γ(x0 , x0 , ω) represents the intensity of light at location x0 at a certain frequency ω. We will assume our illumination sources are fully spatially incoherent within their photon-generating area, leading to a CSD function at source plane L, ΓL (x01 , x02 , ω) = γ 2 C(x01 , ω)δ(x01 − x02 ), (5.13) where C represents the geometric shape of the source intensity for each frequency ω (typically a circfunction in two dimensions), γ is its spatial coherence cross section, and δ is a Dirac delta function. For the remainder of this section, we will drop our interest in the spectral variable ω for simplicity, assuming a notch filter is used in the experiment to effectively isolate a narrow spectrum from the source. Although not detailed here, effects of a spectrally broad (i.e., temporally incoherent) source are an important consideration and may be included through incoherent superposition of the following equations. The Van Cittert-Zernike theorem relates the CSD of the source, ΓL in Eq. (5.13), to the CSD a distance z away, Γz : Γz (∆x) = e −jkq 2z Z jk 0 C(x0 )e z (x ∆x) dx0 ≈ C̃(x0 ), (5.14) where a constant multiplier is neglected for simplicity, ∆x = x1 − x2 and q = x21 − x22 . Assuming (x21 − x22 )/λz << 1 allows us to neglect the phase factor up front. With this assumption, we arrive at an approximate scaled Fourier relationship between the shape of an incoherent illumination source, C, and the CSD function, Γz , at any subsequent plane a large distance z from this source. 5.4.2 CP with partially coherent light In conventional ptychography, the first distant plane the source’s light interacts with is the aperture plane A(x). Here, the light’s CSD function Γ` (x1 − x2 ) is given by Eq. (5.14), with z = `. The aperture a(x) then modulates Γ` (x1 − x2 ) before the light is focused by the lens to the sample plane. 64 As developed in [28], it is direct to show that the CSD at the sample plane, ΓãS , is given by, ΓãS (x01 , x02 ) = Z C(q)ã (x01 − q) ã∗ (x02 − q) dq. (5.15) Here, we used the variable replacement q = x0 for notational clarity. We can update our original expression for the intensity at the detector g(p, x) in Eq. (5.4) to reflect our new partially coherent probe beam with a simple replacement. Instead of multiplying the sample ψ with coherent probe wave ã, we multiply ψ with the probe wave CSD in Eq. (5.15): ZZ g(p, x) = ΓãS (x01 , x02 )ψ(x01 − p)ψ ∗ (x02 − p)exp [−jkx · (x01 − x02 )] dx01 dx02 . (5.16) Plugging Eq. (5.15) into Eq. (5.16) and performing several straightforward manipulations (outlined in Appendix C of [28]) produces the following mathematical description of the CP data matrix g(p, x) in terms of the aperture’s WDF, the sample’s WDF, and the illumination source’s geometric shape C: ZZZ g(p, x) = C(q)Wψ (x0 − p, u)Wã (x0 − q, x − u) dx0 du dq. (5.17) Partially coherent light alters CP’s data matrix with an additional convolution along the scan variable p (Fig. 5.6(a)). The goal of ptychographic data post-processing under partially coherent illumination is to recover a complex description of the sample Wψ from data matrix g(p, x). This is achieved by deconvolving the effect of both Wã and C(q). This is identical to the coherent case, but with an additional (yet still separable) blurring term. 5.4.3 FP with partially coherent light Unlike CP, FP uses an array of spatially offset and partially coherent LEDs at its illumination plane. We now use p, our shift variable, to represent the distance from a given LED to the optical axis. The CSD of one LED may be expressed by modifying Eq. (5.13) to incorporate a spatial offset by p: ΓL (x01 , x02 ) = γ 2 C(x01 − p)δ(x01 − x02 ). This LED’s shifted source light first illuminates the sample at plane S(x). We can write the CSD at the sample plane, ΓS , in terms of ΓL as, Z ΓS (ρ1 − ρ2 ) = 0 C(x0 − p)ejkx (ρ1 −ρ2 ) dx0 = C̃(ρ1 − ρ2 )exp(−jkp(ρ1 − ρ2 )), (5.18) where (ρ1 , ρ2 ) replace (x1 , x2 ) as the spatial coordinates at the sample plane S(x), for notational clarity. As developed in [28], we may transform this CSD through the imaging lens to the detector plane to write, ZZ ΓD (x1 , x2 ) = ΓS (ρ1 − ρ2 )ψ(ρ1 )ψ ∗ (ρ2 )ã(ρ1 − x1 )ã∗ (ρ2 − x2 ) dρ1 dρ2 , (5.19) 65 (a) Conventional ptychography ! g (p, x)   g0 (p, x)   x! x! C(q)! q! blur! =   *   p! p! p! (b) Fourier ptychography! gF0 (p, x)   gF (p, x)   x! x! C(q)! q! 1! blur! =   *   p! p! p! 0! Figure 5.6: Partially coherent light manifests itself as an additional convolution along the data matrix scan dimension p for both (a) CP and (b) FP. The convolution is one-dimensional, as indicated by the vertical bar. With matrices rotated by 90◦ with respect to one another, this convolution will mix the data from each respective setup in a unique manner. For this simulation, we used the same setup parameters as for Fig. 5.2 and Fig. 5.4, but assumed each illumination source C(q) (i.e., LED) is a square with a width of 200 µm (figure adapted from [28]). where ΓD (x1 , x2 ) is the CSD of partially coherent light at the detector. Finally, we may express the measured data matrix in FP in terms of the aperture WDF, sample WDF, and LED source geometry now with, ZZZ gF (p, x) = C(q)Wψ (q − u − p, x0 )Wã (u, x0 − x) dx0 du dq. (5.20) Details of this last step are in Appendix D of [28]. Comparing Eq. (5.20) to Eq. (5.11)’s coherent description of FP, we see that partial coherence manifests itself again as an additional convolution along the data matrix p-dimension (Fig. 5.6(b)). Practically, this indicates that each FP image, captured from a different LED and compiled along p, will begin to look increasingly similar with increasingly incoherent illumination. In the limit of a completely incoherent source, spatial shifting will leave all image features nearly unchanged. Since this blur remains a separable function, it is still possible to deconvolve the effects of both C and Wa to obtain an accurate sample estimate Wψ . Put simply, using a partially coherent source in a CP setup blurs together the sample’s spatial information within its recorded data matrix. In FP, using an array of partially coherent sources instead blurs the sample’s spatial frequency content, as Fig. 5.6 clearly depicts. 66 Deconvolved! Blurred: C = 100μm ! (c) Ptychography and FP recovery error vs. LED size! (a) CP! g (p, x)   RMSEp  =  0.065   x! x! W-1(C, σ)! p! (b) FP! p! RMSEF  =  0.029     gF (p, x)   (d) Input data   x! x! W-1(C, σ)! p! g0(p, x)! 1! x! ψ (x' )! p! -.7! x' (mm)! .7! p! 0! Figure 5.7: Simulation of partially coherent effects produce blurred (a) CP and (b) FP data matrices of an example grating. A Wiener filter can approximately recover the coherent data matrix for each setup, from which an accurate sample reconstruction is direct. (c) Reconstruction error as a function of LED diameter (i.e., blur kernel width) increases for both CP and FP, although FP’s error is consistently lower. (d) The chirped grating sample and its coherent CP data matrix, for comparison (figure adapted from [28]). 5.5 Case study: CP and FP under partially coherent illumination To briefly demonstrate the validity and utility of our phase space model, we now attempt to measure and remove the effects of partial coherence with example CP and FP data matrices. We attempt this both in simulation and experiment. This exercise will allow us to verify the accuracy of Eq. (5.17) and Eq. (5.20). Additionally, it will demonstrate how both setups can successfully operate with low-coherence illumination, assuming that accurate coherence function characterization is possible and minimal noise is introduced. For both simulation and experiment, we carefully designed the scaling and distance parameters to match those listed at the end of Section 3 for three purposes. First, these optimized parameters ensure both data matrices g and gF match, after a rotation. Second, the listed parameters require both setups to use the same lens numerical aperture, detector pixel size and count, and nearly the same total optical path length, offering as even a comparison as possible. Third, the parameters correspond closely with previous optical CP [6, 7] and FP experimental testing platforms. One exception to this close match is the width of the CP’s probe beam at the sample plane, which is typically allowed to be several times wider than what we simulate to allow for under-sampling along p by a similar factor. 67 5.5.1 Simulation In our first investigation, we simulate the partially coherent imaging performance of CP and FP as a function of LED size. Both systems capture 350 one-dimensional images containing 103 pixels each, which combine to form each data matrix. Note that all figures display the central 350-pixel area of each captured image to aid in visualization. As in Fig. 5.2 and Fig. 5.4, our sample here is a chirped grating with minimum feature size of 4 µm. Unlike previous simulations, the grating is now 1.33 mm-wide and is of a slightly different structure to match our experimental sample (see Fig. 5.7(d)). We first apply a Fresnel-based propagation simulation to create this grating’s CP and FP data matrices under partially coherent illumination, as in Fig. 5.7(a)-(b). We then numerically compute Eq. (5.17) and Eq. (5.20) using the same grating function, ψ. In doing so, we find agreement up to an average error of < 1% caused by numerical approximation, which verifies our phase space formulation. Given a valid model, we next test if partial coherence effects can be effectively removed from CP and FP. Successful digital removal of the blurring effects caused by a finite source shape C will allow both setups to maintain high-resolution imaging performance using larger, brighter optical sources (i.e., with higher photon throughput). As a standard benchmark, we apply the well-known Wiener filter in our deconvolution attempt. Previously used to recover complex sample data in [8, 9], it has since been replaced by more advanced phase retrieval-based algorithms [6, 29]. The example blurred CP and FP data matrix inputs in the left of Fig. 5.7(a)-(b) assume quasimonochromatic illumination from sources with 100 µm-diameter active area (0.11◦ angular extent). The associated Wiener deconvolution outputs are shown directly to the right. Gaussian noise (normalized variance of 10−3 ) was added to the data before deconvolution. Noise variance and source size were assumed as prior knowledge to assure optimal filter performance. Fig. 5.7(c) plots the average root-mean-squared error (RMSE) of recovered data matrices as a function of source diameter after Wiener deconvolution. Each point in this plot is an average over 10 experiments with noise variances ranging evenly from 10−2 to 10−4 . The linear process of recovering a sample estimate from its coherent data matrix ensures sample reconstruction RMSE will follow a similar curve. Three important trends are worth noting. First, RMSE increases as a function of LED diameter, but accurate sample recovery is still possible up to quite large-diameter sources. In the tested setup, an angular source extending up to a 0.5◦ maintained manageable error after deconvolution (under realistic noise assumptions). Second, it is easier to globally remove the effects of partial coherence from FP’s data matrix than from CP’s. This key conclusion is a result of the direction of features within the data matrix for this particular simulated object. Blurring occurs along the chirped grating ridges for FP, while it blurs the ridges together for CP, which is harder to invert. Since we expect intensities will vary more quickly along a biological sample’s spatial dimension as opposed to its spatial frequency, the trend of superior FP performance should hold for most samples 68 C = 250 μm! 1! x (mm)! 0! 0! 1.5! x (mm)! (b) Experiment! C = 1000 μm! Intensity (AU)! (a) Simulation! C = 150 μm! 1.5! 0! 0! p (mm)! 2.5! Figure 5.8: (a) Simulated and (b) experimental FP data matrices with varying degrees of partially coherent illumination. The experimental sample closely matches the distribution of ψ(x0 ) in Fig. 5.7(d). C at top indicates the LED active area diameter used in each experiment (figure adapted from [28]). of interest (biological samples tend towards sparse spatial distributions). Third, the blurring induced here by conventional LEDs assumes a rectangularly shaped blur kernel. Due to zeros in its transfer function, this blur kernel is quite challenging to computationally remove. As is well-known within linear filter design, modifying the LED shape and/or apodizing its emission area can remove zeros from its filter and improve deconvolution performance. In general, a broad range of computational imaging systems using either coherent or incoherent active illumination may draw common benefits from considering how coherence alters acquired data, with many computational techniques more advanced than the Wiener filter currently available to improve system design and performance. 5.5.2 Experiment To experimentally verify the findings of the simulation in Fig. 5.7, we constructed a simplified FP setup with an illumination system to scan along one dimension. Experimental parameters closely match the parameters used in simulation (see Section 5.3). Our experimental setup exhibits two primary differences from the diagram in Fig. 5.3. First, a single LED on a motorized linear stage (Newport ESP301) was used instead of a fixed LED array at illumination plane L to facilitate easy variation of LED coherence area. This variation was achieved by placing pinholes of different diameter (100 µm-1000 µm) directly in front of the active area of a 532 nm central-wavelength diode. Second, an f = 50 mm, w = 50 mm collection lens was inserted 50 mm in front of the LED source to assure uniform illumination of the sample. We experimentally determined this lens has minimal effect on the coherence area at the detector plane. Our imaging setup used an f = 105 mm, w = 37.5 mm compound lens (Nikon Micro-Nikkor f/2.8G) positioned do = 300 mm from the sample that 69 imaged onto a 4.54µm pixel CMOS array (Prosilica-GX 1920). Fig. 5.8 displays an example set of simulated and experimental data matrices of the same chirped grating sample in Fig. 5.7 under three different illumination coherence states. Each data matrix was compiled by scanning the LED-pinhole unit at 250 µm steps across 25 mm, capturing an image of the linear grating at each step p, and selecting a single row of the CMOS detector array to form data matrix column p. Each image’s maximum pixel value is scaled to 1 (i.e., each data matrix column in Fig. 5.8 is normalized to it’s maximum value), which enhances the appearance of noise in low-intensity areas but aids with visualization of coherence effects. Fig. 5.8 highlights three primary effects of illumination coherence on FP data. First, the striped “diffraction cone” within each matrix gF (p, x) broadens along the x-dimension when using a largerdiameter source, as the convolution relationship in Eq. (5.20) predicts. Conceptually, an increasingly incoherent source will extend the lens’s coherent spatial frequency cutoff at, k·NA, to its incoherent spatial frequency cutoff at 2k·NA, hence broadening what is captured along p. Second, Eq. (5.20)’s convolution also predicts features along p to blur with increased incoherence, which is clearly observed at the edge of the diffraction cone. Finally, increasingly incoherent illumination still allows the FP setup to acquire high-frequency sample information that otherwise would not be captured by a conventional imaging setup. This is indicated by the dark “tails” at the bottom of each data matrix, which represent the high-frequency grating information that is diffracted into the imaging lens from an off-axis LED, which would otherwise be cutoff from a single image. The density of this high-frequency information tail decreases with increasingly incoherent illumination. However, it is still clearly present with a low-coherence source, thus allowing computational improvement of a reconstructed image’s resolution beyond the conventional imaging lens NA cutoff. This informationpreserving feature of ptychography in the presence of incoherent light is a very powerful tool that has yet to be studied in full, and is the main conclusion of this experiment. 70 Bibliography [1] P. D. Nellist, B. C. McCallum and J. M. Rodenburg, “Resolution beyond the ‘information limit’ in transmission electron microscopy,” Nature 374, 630–632 (1995). [2] F. Hue, J. M. Rodenburg, A. M. Maiden, F. Sweeney and P. A. Midgley, “Wave-front phase retrieval in transmission electron microscopy via ptychography,” Phys. Rev. B 82, 121415(R) (2010). [3] J. M. Rodenburg et al., “Hard-X-ray lensless imaging of extended objects,” PRL 98, 034801 (2007). [4] P. Thibault, M. Dierolf, A. Menzel, O. Bunk, C. David and F. Pheiffer, “High-resolution scanning X-ray diffraction microscopy,” Science 321, 379–382 (2008). [5] M. Dierolf, A. Menzel, P. Thibault, P. Schneider, C. M. Kewish, R. Wepf, O. Bunk and F. Pheiffer, “Ptychographic X-ray computed tomography at the nanoscale,” Nature 467, 437–439 (2010). [6] A. M. Maiden, J. M. Rodenburg and M. J. Humphry, “Optical ptychography: a practical implementation with useful resolution,” Opt. Lett. 35(15), 2585–2587 (2010). [7] A. M. Maiden, M. J. Humphry, F. Zhang and J. M. Rodenburg, “Superresolution imaging via ptychography,” JOSA A 28(4), 604–612 (2011). [8] J. M. Rodenburg and R. H. T. Bates, “The theory of super-resolution electron microscopy via Wigner-distribution deconvolution,” Phil. Trans. R. Soc. Lond. A 339, 521–553 (1992). [9] H. N. Chapman, “Phase retrieval x-ray microscopy by Wigner distribution deconvolution,” Ultramicroscopy 66, 153 (1996). [10] J. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1996). [11] R. Horstmeyer and C. Yang, “A phase space model of Fourier ptychographic microscopy,” Opt. Express 22(1), 338–358 (2014). [12] K. Nugent, “Coherent methods in the X-ray sciences,” Adv. Phys. 59(1), 1–99 (2010). 71 [13] M. Testorf, B. M. Hennelly and J. Ojeda-Castaneda, Phase-Space Optics: Fundamentals and Applications (McGraw-Hill, 2010). [14] M. J. Bastiaans, “Application of the Wigner distribution function to partially coherent light,” JOSA A 3(8), 1227–1238 (1986). [15] R. Horstmeyer, S. B. Oh and R. Raskar, “Iterative aperture mask design in phase space using a rank constraint,” Opt. Express 18(21), 22545–22555 (2010). [16] H. M. L. Faulkner and J. M. Rodenburg, “Movable aperture lensless transmission microscopy: A novel phase retrieval algorithm,” Phys. Rev. Lett. 93, 023903 (2004). [17] M. Bunk, M. Dierolf, S. Kynde, I. Johnson, O. Marti and F. Pfeiffer, “Influence of the overlap parameter on the convergence of the ptychographical iterative engine,” Ultramicroscopy 108, 481–487 (2008). [18] C. Teale, D. Adams, M. Murnane, H. Kapteyn and D. J. Kane, “Imaging by integrating stitched spectrograms,” Opt. Express 21(6), 6783–6793 (2012). [19] J. N. Clark, X. Huang, R. Harder and I. K. Robinsion, “High-resolution three-dimensional partially coherent diffraction imaging,” Nat. Commun. 3, 993 (2012). [20] P. Thibault and A. Menzel, “Reconstructing state mixtures from diffraction measurements,” Nature 494, 68–71 (2013). 72 Chapter 6 Solving ptychography with a convex relaxation As we have demonstrated in the past several chapters, ptychography is a powerful computational imaging technique for transforming a collection of low-resolution images into a high-resolution reconstruction. Unfortunately, algorithms that currently solve this reconstruction problem lack stability, robustness, and theoretical guarantees. In this chapter, we present a convex formulation of the ptychography problem. This formulation has no local minima, it can be solved using a wide range of algorithms, it can incorporate appropriate noise models, and it can include multiple a priori constraints. As also detailed in [1], the new algorithm developed in this chapter relies upon lowrank factorization. Experiments demonstrate that this approach offers a near-linear runtime, linear memory usage and a 25% lower background variance on average than alternating projections, the ptychographic reconstruction algorithm that is currently in widespread use. 6.1 Introduction As we saw in the previous chapter, standard and Fourier ptychography both avoid the need for a large, well-corrected lens to image at the diffraction-limit. Instead, they shift the majority of resolution-limiting factors into the computational realm. Unfortunately, an accurate and reliable solver does not yet exist. As we learned in Chapter 3 for Fourier ptychography, and Chapter 5 for conventional ptychography, high-resolution reconstruction is only achieved if the phase of the scattered field is also determined. Given that both techniques only measure intensity, this presents a challenging problem. To date, most ptychography algorithms recover the unknown phase (i.e., solve the “phase retrieval” problem) by applying known constraints in an iterative manner. We categorize this class of algorithm as an “alternating projection” (AP) strategy. The simplest examples of AP strategies are the Gerchburg-Saxton and error reduction algorithms [2, 4], as reviewed in Chapter 3. Our AP 73 category also includes the iterative projection and gradient search techniques reviewed by Fienup [4] and Marchesini [46], which map to analogous procedures in ptychography [5]. All AP strategies, including several related variants [6, 7], often converge to incorrect local minima or can stagnate [8]. Few guarantees exist regarding convergence, let alone convergence to a reasonable solution. Despite these shortcomings, many authors have pushed beyond the basic algorithms [6] to account for unknown system parameters [22, 52], improve outcomes by careful initialization [12], perform multiplexed acquisition [13], and attempt three-dimensional imaging [34, 36]. Indeed, the images included thus far in this thesis have all resulted from an AP-based recovery scheme. In this chapter, we formulate a convex program for the ptychography problem, which offers an alternative approach to obtain a reliable image reconstruction using efficient computational methods. Convex optimization has recently matured into a powerful computational tool that now solves a variety of challenging problems [16]. However, many sub-disciplines of imaging, especially those involving phase retrieval, have been slow to reap its transformative benefits. Several prior works [17– 21] have connected convex optimization with abstract phase retrieval problems. Here, we consider how convex optimization may improve the quickly growing field of high-resolution ptychography. While it is possible in some experiments to improve reconstruction performance using prior sample knowledge or an appropriate heuristic, we consider here the general case of unaided recovery, which is especially relevant in biological imaging. Under these circumstances, we will show that our convex optimization approach to ptychographic reconstruction has many advantages over AP. Our formulation has no local minima, so we can always obtain a solution with minimum cost. The methodology is significantly more noise-tolerant than AP, and the results are more reproducible. There are also opportunities to establish theoretical guarantees using machinery from convex analysis. Furthermore, there are many efficient algorithms for our convex formulation of the ptychography problem. To obtain solutions at scale, we apply a factorization method due to Burer and Monteiro [22, 23]. This method avoids the limitations of earlier convex algorithms for abstract phase retrieval, whose storage and complexity scale cubically in the number of reconstructed pixels [21]. Moreover, recent results establish that this factorization technique converges globally under certain conditions [24], offering a promising theoretical guarantee. The end result is a new, noise-tolerant algorithm for ptychographic reconstruction that is efficient enough to process the multi-gigapixel images that future applications will demand. Here is an outline for this chapter. First, we manipulate a linear algebra framework for ptychography to pose its sample recovery problem as a convex program. This initial algorithm, termed “convex lifted ptychography” (CLP), supports a-priori knowledge of noise statistics to significantly increase the accuracy of image reconstruction in the presence of noise. Second, we build upon research in low-rank semidefinite programming [22, 23] to develop a second non-convex algorithm, 74 LED Array θj Shift jŸr Sample ψ (x) Aperture a (x' ) Image g (x ) L(x') S(x) A(x') Fourier Ptychography Setup Detector High-resolution complex sample ψ Low-resolution image set g Convex Lifted Ptychography (CLP) … 0! 2π! intensity phase … intensity Figure 6.1: Diagram of the Fourier ptychography setup (top), where we use an LED array to illuminate a sample from different directions and acquire an image set g (bottom). This chapter introduces a convex phase retrieval algorithm to transform this image set into a high-resolution complex sample estimate, ψ. Included image set and reconstructed resolution target are experimental results (figure adapted from [1]). called “low-rank ptychography” (LRP), which improves on the computational profile of CLP. Finally, we explore the performance of LRP in both simulation and experiment to demonstrate how it may be of great use in reducing the image capture time of Fourier ptychography. 6.2 Convex Lifted Ptychography (CLP) 6.2.1 Mathematical fundamentals We derived the Fourier ptychography data matrix in Chapter 2. Here, we present a brief review of this model, which we then convert into a new framework using linear algebra. Following Fig. 6.1, but restricting our attention to just the x dimension to start (extension to two dimensions is direct), we assume that each LED is an effective point emitter that illuminates a sample ψ(x) at a plane S(x) a large distance l away. The jth source illuminates the sample with a spatially coherent plane wave at angle θj = tan−1 (jr/l), where −q/2 ≤ j ≤ q/2. Additionally assuming the sample ψ(x) is 75 thin, we may express the optical field exiting the thin sample as the product, s(x, pj ) = ψ(x)eikxpj . (6.1) Here, the wavenumber k = 2π/λ and pj = sin θj describes the off-axis angle of the jth optical source. This pj th field at the sample plane is then transformed through the microscope to the detector plane, where its intensity is recorded on the digital detector containing m pixels. As we showed in Chapter 2, the resulting image captured at the detector plane can be expressed as, h i2 g(x, pj ) = F a(x0 )ψ̂(x0 − pj ) . (6.2) Here, a(x0 ) is the imaging system aperture function, ψ̂(x0 ) is the Fourier transform of ψ(x), and g(x, pj ) is the (m × q) Fourier ptychography data matrix. The goal of Fourier ptychographic post-processing is to reconstruct a high-resolution (n-pixel) complex spectrum ψ̂(x0 ), from the multiple low-resolution (m-pixel) intensity measurements contained within the data matrix g. Once ψ̂ is found, an inverse-Fourier transform will yield the desired complex sample reconstruction, ψ. As discussed in Chapter 3, most current ptychography setups solve this inverse problem using alternating projections (AP): 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 set from the detector plane (i.e., the values in g). Second, its spectrum ψ̂0 is forced to lie within a known support in the plane that is Fourier conjugate to the detector. Different projection operators and update rules are available, but are closely related [4,5,46]. While these projection strategies are known to converge when each constraint set is convex, the intensity constraint applied at the detector plane is not convex [25], leading to erroneous solutions [26] and possible stagnation [8]. 6.2.2 The CLP solver We begin the process of solving equation 6.2 as a convex program by expressing it in matrix form. First, we represent the unknown sample spectrum ψ̂ as an (n × 1) vector. Again, n is the known sample resolution before it is reduced by the finite bandpass of the lens aperture. Second, the jth detected image becomes an (m × 1) vector gj , where again m is the number of pixels in each low-resolution image. The ratio n/m defines the ptychographic resolution improvement factor. It is equivalent to the largest angle of incidence from an off-axis optical source, divided by the acceptance angle of the imaging lens. Third, we express each lens aperture function a(x + pj ) as an (n × 1) discrete aperture vector aj , which modulates the unknown sample spectrum ψ̂. Note we now shift a by pj , as opposed to ψ by −pj , as the two are mathematically equivalent. 76 Image set b ! 2m pixels! Block-diagonal F! g1(x)   A1! F! =! Sparse A! g2(x)   F! A2! ψ(x')! 2! p1   n pixels! p2   shift! 2m x n! 2m x 2m! Figure 6.2: A set of images captured by Fourier ptychography stack together into a long data vector, b. Each associated matrix transform is similarly stacked and combined to form our final measurement matrix, D = FA. Here, we show stacking of just two images for simplicity. Typically, over 200 images are stacked (figure adapted from [1]). To rewrite equation 6.2 as a matrix product, we define {Aj }qj=1 to be the sequence of (m × n) rectangular matrices that contain a deterministic aperture function aj along a diagonal. For an aberration-free rectangular aperture, each matrix Aj has a diagonal of ones originating at (0, p0j ) and terminating at (m, p0j + m − 1), where p0j is now a discretized version of our shift variable pj . Finally, we introduce an m×m discrete Fourier transform matrix F(m) to express the transformation of the low-pass filtered sample spectrum through our fixed imaging system for each low-resolution image gj : 2 gj = F(m) Aj ψ̂ , 1 ≤ j ≤ q. (6.3) Ptychography acquires a series of q images, {gj }qj=1 . We combine this image set into a single vector by “stacking” all images in equation 6.3: 2 b = FAψ̂ 2 = Dψ̂ . (6.4) Here, b is {g} expressed as a (q · m × 1) stacked image vector (see Fig. 6.2). In addition, we define D = FA, where F is a (q · m × q · m) block diagonal matrix containing q copies of the low-resolution DFT matrices F(m) in its diagonal blocks, and A has size (q · m × n) and is formed by vertically stacking each aperture matrix Aj :  F(m)   . F =  ..  0 ··· .. . ···  A1     .    , A =  ..  .    F(m) Aq 0 .. .   (6.5) We denote the transpose of the ith row of D as di , which is a column vector. The set {di } forms our measurement vectors. The measured intensity in the ith pixel is the square of the inner product 77 D E between di and the spectrum ψ̂: bi = | di , ψ̂ |2 . Next, we “lift” the solution ψ̂ out of the quadratic relationship in equation 6.4. As suggested in [19], we may instead express it in the space of (n × n) positive-semidefinite matrices:     bi = Tr ψ̂ ∗ di d∗i ψ̂ = Tr di d∗i ψ̂ ψ̂ ∗ = Tr (Di X) , (6.6) where Di = di d∗i is a rank-1 measurement matrix constructed from the ith measurement vector di , X = ψ̂ ψ̂ ∗ is an (n × n) rank-1 outer product, and 1 ≤ i ≤ q · m. Equation 6.6 states that our quadratic image measurements {bi }q·m i=1 are linear transforms of ψ̂ in a higher dimensional space. We may combine these q · m linear transforms into a single linear operator A to summarize the relationship between the stacked image vector b and the matrix X as, A (X) = b. One can now pose the phase retrieval problem in ptychography as a rank minimization procedure: minimize rank(X) subject to A (X) = b, (6.7) X  0, where X  0 denotes X is positive-semidefinite. This rank minimization problem is not convex and is a computational challenge. Instead, adapting ideas from [18], we form a convex relaxation of equation 6.7 by replacing the rank of matrix X with its trace. This creates a convex semidefinite program: minimize Tr(X) subject to A (X) = b, (6.8) X  0. Several recent results establish that the relaxation in equation 6.8 is equivalent to equation 6.7 under certain conditions on the operator A [27, 28]. Although not necessarily equivalent in general, this relaxation consistently offers us highly accurate experimental performance. To account for the presence of noise, we may reform equation 6.8 such that the measured intensities in b are no longer strictly enforced constraints, but instead appear in the objective function: minimize 1 α Tr(X) + kA (X) − bk 2 (6.9) subject to X  0. Here, α is a scalar regularization variable that directly trades off goodness for complexity of fit. Its optimal value depends upon the assumed noise level. Equation 6.9 forms our final convex problem to recover a resolution-improved complex sample ψ from a set of obliquely illuminated images in b. Many standard tools are available to solve this semidefinite program (see Appendix A). Its solution 78 (a) Sample, simulated! (c) CLP Reconstruction! (d) AP Reconstruction! (non-convex)! 1! 1! 1! 0! 0! 0! π! π! π! -π! -π! -π! Phase! Amplitude! (convex)! (b) Simulated raw image set! (12 x 12 pixels)! …! …! Figure 6.3: Simulation of the CLP algorithm. (a) An n = 36 × 36 pixel complex sample (simulated) consisting of absorptive microspheres modulated with an independent quadratic phase envelope. (b) Sequence of low-resolution simulated intensity measurements (m = 12 × 12 pixels each), serving as algorithm input. (c)-(d) Example CLP and AP reconstructions, where CLP is successful but AP converges to an incorrect local minimum. Here we use q = 82 images to achieve a resolution gain of 3 along each spatial dimension and simultaneously acquire phase (figure adapted from [1]). defines our Convex Lifted Ptychography (CLP) approach. In practice, CLP returns a low-rank matrix X, with a rapidly decaying spectrum, as the optimal solution of equation 6.9. The trace term in the CLP objective function is primarily responsible for enforcing the low-rank structure of X. While this trace term also appears like an alternative method to minimize the unknown signal energy, we caution that a fair interpretation must consider its effect in a lifted (n × n) solution space. We obtain our final complex image estimate ψ by first performing a singular value decomposition of X. Given low-noise imaging conditions and spatially coherent illumination, we set ψ to the Fourier transform of the largest resulting singular vector. Viewed as an autocorrelation matrix, we may also find useful statistical measurements within the remaining smaller singular vectors of X. We note that one may also identify X as the discrete mutual intensity D E matrix of a partially coherent optical field: X = ψ̂ ψ̂ ∗ , where hi denotes an ensemble average [29]. Under this interpretation, equation 6.9 becomes an alternative solver for the stationary mixed states of a ptychography setup [29]. Without any further modification, three points distinguish equation 6.9 from existing AP-based ptychography solvers. First, the convex solver has a larger search space. If AP is used to iteratively update an n-pixel estimate, equation 6.9 must solve for an n × n positive-semidefinite matrix. Second, this boost in the solution space dimension guarantees the convex program may find its global optimum with tractable computation. This allows CLP to avoid AP’s frequent convergence to local 79 minima (i.e., failure to approach the true image). Unlike prior solvers for the ptychography problem, no local minima exist in the CLP approach. However, CLP cannot yet claim a single global minimum, since it is not necessarily a strictly convex program. Finally, equation 6.9 implicitly considers the presence of noise by offering a parameter (α) to tune with an assumed noise level. AP-based solvers lack this parameter and can be easily led into incorrect local minima by even low noise levels, which we demonstrate next. 6.2.3 CLP simulations and noise performance We simulate Fourier ptychography following the setup in Fig. 6.1. We capture multiple two- dimensional images in (x, y) from a three-dimensional optical geometry. The simulated FP setup contains a detector with m = 122 pixels that are each 4 µm wide, a 0.1 NA lens at plane A(x0 , y 0 ) (6◦ collection angle, unity magnification), and an array of spatially coherent optical sources at plane L(x0 , y 0 ) (632 nm center wavelength, 10 nm spectral bandwidth). The array is designed to offer an illumination NA of 0.2 (θmax = 11.5◦ maximum illumination angle). Together, the lens and illumination NAs define the reconstructed resolution of our complex sample as n = 362 pixels, increasing the pixel count of one raw image by a factor n/m = 9. Figure 6.3(b) shows example simulated raw images from a sample of absorptive microspheres modulated by a quadratic phase envelope. Within each raw image, the set of microspheres is not clearly resolved. Here, we simulate the capture of q = 82 low resolution images, each uniquely illuminated from one of q = 82 optical sources in the square array. We then input this image set into both the standard AP algorithm (i.e., the PIE strategy) [6], as well as CLP in equation 6.9, to recover a high resolution (36 × 36 pixel) complex sample. Here, we select the PIE strategy as our comparison benchmark for two reasons. First, it is one of the most widely used ptychography algorithms. Second, similar to CLP, its structure implicitly assumes a Gaussian noise model [5]. Even in the noiseless case, five iterations of nonlinear AP introduces unpredictable artifacts to both the recovered amplitude and phase (Fig. 6.3(d)), while CLP offers near perfect recovery (Fig. 6.3(c)). A constant phase offset is subtracted from both reconstructions for fair comparison, and we selected α = .001. Next, we quantify AP and CLP performance. We repeat the reconstructions in Fig. 6.3, again setting α = .001 in equation 6.9 while varying two relevant parameters: the number of captured images q, and their signal-to-noise ratio (SNR). We define the SNR as, SNR = 10 log10 ( |ψ|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, we follow [52] and define our reconstruction mean-squared error as P P P 2 P 2 2 MSE = x |ψ(x) − ρs(x)| / x |ψ(x)| , where ρ = x ψ(x)s∗ (x)/ x |s(x)| is a constant phase factor shifting our reconstructed phase to optimally match the known phase of the ground truth 80 Log MSE (dB)! Reconstruction MSE vs. SNR (simulation with quadratic phase)! AP   LRP  (a v g.)   CLP   Log SNR (dB)! Figure 6.4: Reconstruction MSE versus signal to noise ratio (SNR) of CLP and AP (log scale, dB). Each curve represents reconstruction with a different number of captured images, q, corresponding to a different percentage of spectrum overlap (ol, noted in legend). Each point is an average over five independent algorithm runs with unique additive noise. Also included is the average performance of our LRP algorithm over the same three spectrum overlap settings (see Section 6.4, figure adapted from [1]). sample. Figure 6.4 plots MSE as a function of SNR for this large set of CLP and AP reconstructions. Each of the algorithms’ 3 independent curves simulates reconstruction using a different number of captured images, q. We summarize q by defining a Fourier spectrum overlap percentage: ol = 1 − (n − m)/qm. Each of the 6 points within one curve simulates a different level of additive measurement noise. Each point is an average over 5 independent trials. Since AP tends not to converge in the presence of noise, we represent each AP trial with the reconstruction that offers the lowest MSE across all iteration steps (up to 20 iterations). All CLP reconstructions improve linearly as SNR increases, while AP performance fluctuates unpredictably. For both algorithms, performance improves with increased spectrum overlap ol, and reconstruction fidelity quickly deteriorates and then effectively fails when ol drops below ∼60%. 6.3 Factorization for Low-Rank Ptychography (LRP) Posing the inverse problem of ptychography as a semidefinite program (equation 6.9) is a good first step towards a more tractable solver. However, the constraint that X remain positive-semidefinite is computationally demanding: each iteration typically requires a full eigenvalue decomposition of X. As the size of X scales with n2 , processable image sizes are limited to an order of 104 pixels 81 on current desktop machines. This scaling limit does not prevent large-scale CLP processing of ptychography data. It is common practice to segment each detected image into as few as 103 pixels, process each segment separately, and then “tile” the resulting reconstructions back together into a final full-resolution solution [32]. CLP may also parallelize its computation with this strategy. 6.3.1 The LRP solver While such tiling parallelization offers significant speedup, a simple observation helps avoid the poor scaling of CLP altogether: the desired solution of the ptychography problem in equation 6.7 is lowrank. Instead of solving for an n × n matrix X, it is thus natural to adopt a low-rank ansatz and factorize the matrix X as X = RRT , where our decision variable R is now an n×r rectangular matrix containing complex entries, with r < n [22, 23]. Inserting this factorization into our optimization problem in equation 6.8 and writing the constraints in terms of the measurement matrix Di = di dTi creates the non-convex program, minimize Tr(RRT ) subject to Tr(Di RRT ) = bi (6.10) for all i. Besides removing the positive semidefinite constraint in equation 6.8, the factored form of equation 6.10 presents two more key adjustments to our original convex formulation. First, using the relationship Tr(RRT ) = kRk2F , where F denotes a Frobenius norm, it is direct to rewrite the objective function and each constraint in equation 6.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 = ψ̂ ψ̂ ∗ , so we expect R as a column vector and RRT a rank-1 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 equation 6.10 is no longer convex. We thus must shift our solution strategy away from a simple semidefinite program. Prior work in [22, 23] suggests that an efficient and practically successful method of solving equation 6.10 is to minimize the following augmented Lagrangian function: L(R, y, σ) = Tr(RRT ) − X
σ X
2 yi · Tr(Di RRT ) − bi + · Tr(Di RRT ) − bi , i i 2 (6.11) where R ∈ Cn×r is the unknown decision variable and the two variables y ∈ Rq·m and σ ∈ R+ are new parameters to help guide our algorithm to its final reconstruction. The first term in equation 6.11 is the objective function from equation 6.10, indirectly encouraging a low-rank factorized product. 82 This tracks our original assumption of a rank-1 solution within a “lifted” solution space. The second term contains the known equality constraints in equation 6.10 (i.e., the measured intensities), each assigned a weight yi . The third term is a penalized fitting error that we abbreviate with label v. It is weighted by one penalty parameter σ, mimicking the role of a Lagrangian multiplier. With an appropriate fixed selection of yi ’s and σ, the minimization of L(R, y, σ) with respect to R identifies our desired optimum of equation 6.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 is a guaranteed solution [23]. As an unconstrained function, the minimum of L is quickly found via standard tools (e.g., a quasi-Newton approach such as the LBFGS algorithm [32]), as previously demonstrated across a wide range of applications and experiments [22]. The goal of our low-rank ptychography (LRP) algorithm thus reduces to the following task: determine a suitable set of (yi , σ) such that we may minimize equation 6.11 with respect to R, which leads to our desired solution. We use the iterative algorithm suggested in [22] to sequentially minimize L with respect to Rk at iteration k, and then update a new parameter set (y k+1 , σ k+1 ) at iteration k + 1. We update parameters (y k+1 , σ k+1 ) to ensure their associated term’s contribution to the summation forming L is relatively small. This suggests Rk+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 we hold its multiplier σ fixed and update the equality constraint multipliers, yi . Otherwise, we increase σ by a factor γ such that σ k+1 = γσ k . A detailed discussion of these algorithmic steps is in [22, 23]. We initialize the LRP algorithm with an estimate of the unknown high-resolution complex sample function ψ0 , contained within a low-rank matrix R0 . We terminate the algorithm either if it reaches a sufficient number of iterations, or if the minimizer fulfills some convergence criterion. We form R0 using a spectral method, which can help increase solver accuracy and decrease computation time [31]. Specifically, we select the r columns of R0 as the leading r eigenvectors of D∗ diag[b]D, where D is the measurement matrix in equation 6.4. While this spectral approach works quite well in practice, a random initialization of R0 also often produces an accurate reconstruction. 6.3.2 LRP simulations and noise performance Following the same procedure used to simulate the CLP algorithm, we test the MSE performance of the LRP algorithm as a function of SNR in Fig. 6.5. We again add different amounts of uncorrelated Gaussian noise to each simulated raw image set and compare the LRP reconstruction with ground truth. This simulated sample is the experimentally obtained amplitude and phase of a human blood smear. It is qualitatively similar to the sample used in Fig. 6.3. Unlike with the simulations in Figs. 6.3–6.4, the AP algorithm no longer malfunctions at lower spectrum overlap percentages (i.e., 83 (a) Example simulation: LRP vs. AP! Sample, simulated! LRP Reconstruction! AP Reconstruction! Amplitude! 1! 0! Phase! π! -π! Log MSE (dB)! (b) Reconstruction MSE vs. SNR (red blood cell simulation)! ! AP   LRP   Log SNR (dB)! Figure 6.5: Simulation of the LRP and AP algorithms using the same parameters as for Figs. 6.3–6.4, but now with a different “red blood cell” sample. (a) Using 82 simulated intensity measurements as input (SNR=19, 122 pixels each), both algorithms successfully recover each cell’s phase, but AP is less accurate. (b) MSE versus SNR plot with varying amounts of noise added to the same data set. The MSE for LRP is ∼5-10 dB lower than for AP across all noise levels and aperture overlap settings (each point from 5 independent trials, figure adapted from [1]). lower values of ol). Despite this apparent success, the MSE of the LRP minimizer is still ∼5-10 dB better than the MSE of the AP minimizer, across all levels of SNR. This reduced LRP reconstruction error follows without any parameter optimization or explicit noise modeling. In these simulations, we somewhat arbitrarily fix η and γ at 0.5 and 1.5, respectively, and set the desired rank of the solution, r, to 1. In practice, these free variables offer significant freedom to tune the response of LRP to noise. For example, similar to the noise parameter α in equation 6.9, the multiplier σ (controlled via γ) in equation 6.11 helps trade off complexity for goodness of fit by re-weighting the quadratic fitting error term. In addition to reducing required memory, the LRP algorithm also improves upon the computational cost of CLP. For an n-pixel sample reconstruction, the per iteration cost of the CLP algorithm Third  Version?  -­‐48x48x3x3:   84 (a) AP reconstruction, full! (b) LRP reconstruction, full! 20 μm! 20 μm! 1! Intensity (AU)! 0! 8 μm! 8 μm! (c) Raw data! 8 μm! Figure 6.6: Experimental reconstruction of a USAF target, where the number of resolved pixels is increased by a factor of 25. We test two different ptychography algorithms: (a) AP and (b) LRP. Here we only show reconstructed intensity. LRP avoids artifacts (e.g., boxed in green) commonly encountered in the AP approach. Cited variances are measured in blue boxes (top). (c) Same cropped region of one low-resolution raw image, for comparison (figure adapted from [1]). is currently O(n3 ), using big-O notation. The positive-semidefinite constraint in equation 6.9, which requires a full eigenvalue decomposition, defines this behavior limit. The per-iteration cost of the LRP algorithm, on the other hand, is O(n log n). This large per-iteration cost reduction is the primary source of LRP speedup. For example, LRP required ∼21 seconds to complete an average simulation of the example in Fig. 6.3, while CLP required ∼170 minutes and AP required 1 second on the same desktop machine. 6.4 Experiments We experimentally verify how LRP improves the accuracy and noise stability of ptychographic reconstruction using a Fourier ptychographic (FP) microscope. Our experimental procedure closely follows the protocol in [32]. While we demonstrate at optical wavelengths, it is straightforward to acquire a Fourier ptychographic data set in an X-ray or electron microscope (e.g., with a tilting source [34]). Alternatively, two trivial changes within equation 6.10 directly prepares standard ptychographic data for LRP processing (see end of section 2). Given its removal of local minima and improved treatment of noise, we expect our strategy will benefit both experimental arrangements. In this section, we first quantitatively verify that LRP accurately measures high resolution and sample phase. Compared with AP reconstructions, our LRP algorithm generates fewer undesirable artifacts in experiment. Second, we will compare the AP and LRP reconstructions of a biological 85 sample, which establishes the improved noise stability of our new algorithm. 6.4.1 Quantitative performance Our FP microscope consists of a 15×15 array of surface-mounted LEDs (model SMD 3528, center wavelength λ=632 nm, 4 mm LED pitch, 150 µm active area diameter), which serve as our quasicoherent optical sources. The array is placed l=80 mm beneath the sample plane, and each LED has an approximate 20 nm spectral bandwidth. Prior work establishes that the impact of non-ideal source coherence is gradual [33]. While negligible in these experiments, we may eventually account for source statistics in the multi-rank structure of the LRP optimizer R. To quantitatively verify resolution improvement, we turn on each of the 15 × 15 LEDs beneath a U.S. Air Force (USAF) resolution calibration target. A 2X Olympus microscope objective (apochromatic Plan APO 0.08 NA) transfers each resulting optical field to a CCD detector (Kodak KAI-29050, 5.5 µm pixels), creating 225 low resolution images. Using this 0.08 NA microscope objective (5◦ collection angle) and a 0.35 illumination NA (θmax = 20◦ illumination angle), our FP microscope offers a total complex field resolution gain of n/m = 25. Each image spectrum overlaps by ol ≈ 70% in area with each neighboring image spectrum. For reconstruction, we select n = 25 · m and use the same aperture parameters with AP and LRP to create the high-resolution images in Fig. 6.6. For computational efficiency, we segment each lowresolution image into 3×3 tiles (n=4802 per tile) and process the tiles in parallel, as performed in [32]. We determine the optimal number of AP and LRP algorithm iterations as 6 and 15, respectively, and fixed this for each tile (and all subsequent reconstructions). We typically initialize LRP with the following parameters: γ=1.5, η=0.3, y 0 =10 and σ 0 =10. We determine the microscope aperture function with an iterative procedure [30] before each experiment and fix it for each algorithm trial. Both ∼1 megapixel reconstructions achieve their maximum expected resolving power (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 (Group 7, Element 2 in Fig. 6.6(c)). Our new LRP algorithm avoids certain artifacts that are commonly observed during the nonlinear descent of AP (boxed in green). Both reconstructions slowly fluctuate in background areas that we expect 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 Fig. 6.6, AP and LRP exhibit normalized 2 2 background amplitude variances of σA = 5.4 × 10−4 and σL = 5.0 × 10−4 , respectively. Accounting for experimental uncertainty in the aperture function a and shifts pj (e.g., following [30, 33]) may reduce this error in both algorithms. To verify that our LRP solver reconstructs quantitatively accurate phase, we next image a monolayer of polystyrene microspheres (index of refraction nm = 1.587) immersed in oil (no = 1.515, both 86 π! -π! 1 μm! (c) Measured thickness! (b) LRP phase map! 1 μm! 6! Thickness (μm)! (a) AP phase map! 4! 2! 0! -5! Width (μm)! 5! Figure 6.7: Experimental measurement of the quantitative optical phase emerging from two polystyrene microspheres. Both (a) AP and (b) LRP reconstruct phase maps that appear qualitatively similar, although the AP phase map flattens at the sphere’s center. Variances measured in red boxes. (c) Plot of microsphere thickness from a trace through the center of the large sphere (dashed line) demonstrates close agreement between LRP and ground truth (GT). Figure adapted from [1]. indexes for λ = 632 nm light). To demonstrate the LRP algorithm easily generalizes to any ptychographic arrangement, we perform this experiment on a new “high-NA” FP microscope. The high-NA setup uses a larger 0.5 NA microscope objective lens with a 30◦ collection angle (20X Olympus 0.5 NA UPLFLN). For sample illumination, we now arrange 28 LEDs into three concentric rings of 8, 8, and 12 evenly spaced light sources (ring radii=16, 32 and 40 mm, respectively). We place this new light source array 40 mm beneath the sample to create a 0.7 illumination NA with a θmax = 45◦ illumination angle. The synthesized numerical aperture of this new FP microscope, computed as the sum of the illumination NA and objective lens NA, is NAs = 1.2. With a greater-than-unity synthetic NA, our reconstructions can offer oil-immersion quality resolution (∼385 nm smallest resolvable feature spacing [36]), without requiring any immersion medium between the sample and objective lens. Using the same data and parameters for AP and LRP input, we obtain the high-resolution phase reconstructions of two adjacent microspheres in Fig. 6.7 (3 µm and 6 µm diameters). For this reconstruction, we set m=1602 and n=3202 . We have subtracted a constant phase offset from the LRP solution in (b) to allow for direct comparison to the AP solution in (a). The two reconstructions appear qualitatively similar except at the center of the 6 µm sphere, where the AP phase profile unexpectedly flattens. We highlight this flattening by selecting phase values along each marked dashed line to plot the resulting sample thickness in Fig. 6.7(c). Phase φ and sample thickness t are related via t = k∆φ(nm − no )−1 , where k is the average wavenumber and ∆φ = φ − φ0 is the reconstructed phase minus a constant offset. LRP closely matches the optical thickness of a groundtruth sphere (GT, black curve): the length of the vertical chord connecting the top and bottom arcs of a 6 µm diameter circle. The normalized amplitude variances from a 402 -pixel background area 87 AP reconstruction! LRP reconstruction! Full FOV! 1 second ! 1/3 second ! Intensity variance vs. Exposure! 1/10 second! Figure 6.8: Experimental reconstructions of a resolution target imaged by our high-NA FP microscope. Results from three different exposure times are shown. (lower left) Plot of the normalized background variance of the intensity of our reconstructed image, as a function of exposure time for 5 different data sets. The background variance of the LRP reconstruction is consistently lower than the variance of the AP reconstruction. The difference between the background variance from each algorithm increases at shorter exposure times, where LRP outperforms AP by roughly a factor of 2. 2 2 = 5.8 × 10−4 , respectively. This again supports our observation that = 9.2 × 10−4 and σL are σA the high resolution reconstructions formed by LRP are more accurate than those formed by AP. For a third quantitative experiment, we use the same high-NA FP setup (collection angle= 30◦ , θmax = 45◦ ) to image a resolution target with varying exposure times. The resulting reconstructions are shown in Fig. 6.8. We assume a shorter exposure time implies increased noise within each raw image. We include reconstructions from three data sets: images captured with a 1 second exposure (top), a 0.25 second exposure (middle), and 0.1 second exposure (bottom). It is clear that the LRP solver maintains a smoother solution as additional noise is introduced. This observation is quantitatively verified by selecting a 20 by 20 pixel background area where the resolution target should be a constant value, and computing the normalized variance for reconstructions based on 5 different exposure times, ranging from 0.05 to 1 second. The results are plotted in the lower left. Here, we see that LRP maintains a consistently lower background variance, which becomes much lower than the result of AP (roughly a factor of 2) at low exposure times (high noise). However, we note that the LRP result still maintains a slowly varying, uneven background intensity. The source 88 (a) Oil immersion (1.25 NA) image (b) LRP, full-FOV (c) Raw image 4 μm 4 μm 2 μm (e) LRP reconstructions 2 μm 2 μm 0.25 seconds 2 μm 2 μm 2 μm 2 μm 1 second (d) AP reconstructions 0.1 seconds 2 μm 0 Amp. 1 0 Phase π 0 Amp. 1 0 Phase π Figure 6.9: Experimental reconstruction of malaria-infected human red blood cells. (a) Oil immersion microscope image (1.25 NA) identifies two infected cells of interest (marked with arrows). (b) Example LRP reconstruction (area of interest in red box). (c) One example raw image used for AP and LRP algorithm input. (d) AP-reconstructed amplitude and phase from three different 29-image data sets, using 1 sec (top), 0.25 sec (middle) and 0.1 sec (bottom) exposure times for all images in each set. Variances measured within green boxes. Increased noise within short-exposure images deteriorates reconstruction quality until both parasites are not resolved. (e) LRP reconstructions using the same three data sets. Both parasites are clearly resolved in the reconstructed phase for all three exposure levels (figure adapted from [1]). of this background term is the result of error in reconstructing lower image spatial frequencies. 6.4.2 Biological sample reconstruction Our fourth imaging example uses our high-NA FP configuration to resolve a biological phenomenon: the infectious spread of malaria in human blood. The early stages of a Plasmodium falciparum infection in erythrocytes (i.e., red blood cells) includes the formation of small parasitic “rings”. It is challenging to resolve these parasites under a microscope without using an immersion medium, even after appropriate staining. Oil-immersion is required for an accurate diagnosis of infection [37]. We use FP to resolve Plasmodium falciparum-infected cells with a 0.5 NA objective lens and using no oil in Fig. 6.9. We first prepare an infected blood sample following the protocol in [38]: we 89 maintain erythrocyte asexual stage cultures of the P. falciparum strain 3D7 in culture medium, then we smear, fix with methanol and apply a Hema 3 stain. An example sample region containing two infected cells, imaged with a conventional high-NA oil-immersion microscope (NA = 1.25) under Kohler illumination, is in Fig. 6.9(a). Next, we capture 28 uniquely illuminated images of these two infected cells using our high-NA FP microscope. Figure 6.9(c) contains an example normally illuminated raw image, which does not clearly resolve the parasite infection. Figure 6.9(d) presents phase retrieval reconstructions using the standard AP algorithm, where we set m=1202 , n=2402 , run six iterations, and again subtract a constant phase offset. We include reconstructions from three data sets: images captured with a 1 second exposure (top), a 0.25 second exposure (middle), and 0.1 second exposure (bottom). As with our third example, a shorter exposure time implies increased noise within each raw image. While the 1 sec exposure-based AP reconstruction resolves each parasite, blurred cell boundaries and non-uniform 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 402 -pixel window (marked blue square), is 2 = .0020, .0027, and .0059 for the 1 sec, 0.25 sec, and 0.1 sec exposure reconstructions, respectively. σA For comparison, reconstructions using our LRP algorithm are shown in Fig. 6.9(e) (sharpest solutions after 15 iterations). For each reconstructed amplitude, we set the desired solution rank to r = 3. We add the 3 modes of the resulting reconstruction in an intensity basis to create the displayed amplitude images. For each reconstructed phase, we set the desired solution matrix rank to r = 1 and leave all other parameters 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 2 = .0016 (1 sec), .0022 (0.25 amplitude variances from the same background window are now σL sec), and .0035 (0.1 sec), an average reduction (i.e., improvement) of 26% with respect to the AP results. While not observed within our previous simulations or experiments, 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 post-processing, the AP algorithm here might offer superior quantitative analysis during, e.g., tomographic cell reconstruction, where low-order phase variations must remain accurate. However, it is clear within Fig. 6.9 that LRP better resolves the fine structure of each infection, which is critical during malaria diagnosis. A shorter image exposure time (i.e., up to 10 times shorter) may still enable accurate infection diagnosis when using LRP, as opposed to the standard AP approach. 90 6.5 Discussion and Conclusion Through the relaxation in equation 6.8, we first transform the traditionally nonlinear phase retrieval process for ptychography into a convex program. We may now use well-established optimization tools to find the ptychography problem’s global minimum. Then, we suggest a practically efficient algorithm to solve the resulting semidefinite program with an appropriate factorization. The result is a new ptychographic image recovery algorithm that is robust to noise. We demonstrate its successful performance in three unique experiments, concluding with a practical biological imaging scenario: the identification of malaria infection without using an oil immersion medium and under shortexposure imaging conditions. Much future work remains to fully explore the specific benefits of our problem reformulation. Besides removing local minima from the recovery process, perhaps the most significant departure from prior 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 third-generation X-ray synchrotron sources and within electron microscopes. An appropriately selected r may eventually help LRP measure the partial coherence of such sources, as outlined in [29]. The solution rank may also help identify setup vibrations, sample auto-fluorescence, or even 3D sample structure. As in prior work with low-rank matrix optimization, we may also artificially enlarge our 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 inserting additional sample priors such as sparsity. These refinements are currently a critical component of ptychographic recovery in the fields of X-ray and electron microscopy, and will also improve our optical results. Along with algorithm refinement, a detailed comparison between LRP and various other recovery methods, especially under different sources of noise and error, will help identify the experimental conditions under which our strategy is of greatest benefit. What’s more, as a particular solution to the general problem of phaseless measurement, our findings can also inform a wide variety of coherent diffractive imaging techniques. Regardless of the specific experimental application, convex analysis will continue to provide useful theoretical guarantees regarding phase retrieval algorithm performance, a crucial feature missing from current nonlinear solvers. Appendix A. Computational specifics We performed all processing on a high-end desktop containing two Intel Xeon 2.0 GHz CPUs and two 3GB GeForce GTX GPUs. Code was written in Matlab with built-in GPU acceleration. We 91 solved our CLP semidefinite program using the TFOCS code package [39]. Our LRP algorithm borrows concepts from the LBFGS solver in [32] for one specific minimization step. LRP’s total recovery time for the 1 megapixel example in Fig. 6.6 was approximately 130 seconds, while AP completed in approximately 15 seconds on the same desktop. 92 Bibliography [1] R. Horstmeyer, R. Y. Chen, X. Ou, B. Ames, J. A. Tropp and C. Yang, “Solving ptychography with a convex relaxation,” New. J. Phys. 17, 053004 (2015). [2] R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik 35, 237–246 (1972). [3] J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 22, 2758–2769 (1982). [4] S. Marchesini, “A unified evaluation of iterative projection algorithms for phase retrieval,” Rev. Sci. Instrum. 78, 011301 (2007). [5] P. Godard, M. Allain, V. Chamard and J. Rodenburg, “Noise models for low counting rate coherent diffraction imaging,” Opt. Express 20(23), 25914–25934 (2012). [6] V. Elser, I. Rankenburg and P. Thibault, “Searching with iterated maps,” Proc. Natl. Acad. Sci. 104, 418–423 (2007). [7] P. Thibault and M. Guizar-Sicairos, “Maximum-likelihood refinement for coherent diffractive imaging,” New. J. Phys. 14, 063004 (2012). [8] J. R. Fienup, C. C. Wackerman, “Phase-retrieval stagnation problems and solutions,” J. Opt. Soc. Am. A 3, 1897–1907 (1986). [9] H. M. L. Faulkner and J. M. Rodenburg, “Movable aperture lensless transmission microscopy: A novel phase retrieval algorithm,” Phys. Rev. Lett. 93, 023903 (2004). [10] A. M. Maiden and J. M. Rodenburg, “An improved ptychographical phase retrieval algorithm for diffractive imaging,” Ultramicroscopy 109, 1256–1562 (2009). [11] A. M. Maiden, M. J. Humphry, M. C. Sarahan, B. Kraus and J. M. Rodenburg, “An annealing algorithm to correct positioning errors in ptychography,” Ultramicroscopy 120, 64–72 (2012). [12] S. Marchesini, A. Schirotzek, C. Yang, H. Wu and F. Maia, “Augmented projections for ptychographic imaging,” Inverse Probl. 29, 115009 (2013). 93 [13] L. Tian, X. Li, K. Ramchandran and L. Waller, ”Multiplexed coded illumination for Fourier ptychography with an LED array microscope,” Biomed. Opt. Express 5, 2376–2389 (2014). [14] T. M. Godden et al., “Ptychographic microscope for three-dimensional imaging,” Opt. Express 22, 12513–12524 (2014). [15] L. Tian and L. Waller, “3D intensity and phase imaging from light field measurements in an LED array microscope,” Optica 2, 104-111 (2015). [16] S. Boyd and L. Vandenberghe, Convex Optimization. Cambridge University Press (2009). [17] H. H. Bauschke, P. L. Combettes and D. R. Luke, “Phase retrieval, error reduction algorithm, and Fienup variants: a view from convex optimization,” J. Opt. Soc. Am. A 19, 1334–1345 (2002). [18] M. Fazel, Matrix rank minimization with applications. PhD thesis, Stanford University, Palo Alto, CA (2002). [19] R. Balan, B. G. Bodmann, P. G. Casazza and D. Edidin, “Painless reconstruction from magnitudes of frame coefficients,”. J. Fourier Anal. Appl. 15, 488–501 (2009). [20] Y. Shechtman, Y. C. Eldar, A. Szameit and M. Segev, “Sparsity based sub-wavelength imaging with partially incoherent light via quadratic compressed sensing,” Opt Express 19, 14807–14822 (2011). [21] E. J. Candes, Y. C. Eldar, T. Strohmer and V. Voroninsk, “Phase retrieval via matrix completion,” SIAM Journal on Imaging Sciences 6, 199–225 (2012) [22] S. Burer and R. D. C. Monteiro, “A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization,” Math. Program., Ser. B 95, 329–357 (2003). [23] S. Burer and R. D. C. Monteiro, “Local minima and convergence in low-rank semidefinite programming,” Math. Program., Ser. A 103, 427–444 (2005). [24] C. De Sa, K. Olukotun and C. Re, “Global convergence of stochastic gradient descent for some non convex matrix problems,” arxiv:1411.1134v2 (2014). [25] H. H. Bauschke, P. L. Combettes and D. R. Luke, “Hybrid projection-reflection method for phase retrieval,” J. Opt. Soc. Am. A 20, 1025–1024 (2003). [26] J. L. C. Sanz, T. S. Huang and T. Wu, “A note on iterative Fourier transform phase reconstruction from magnitude,” IEEE Trans Acoust Speech Signal Processing 32, 1251–1254 (1984). [27] B. Recht, M. Fazel and P. A. Parrilo, “Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization,” SIAM Rev. 54, 471–501 (2010). 94 [28] E. J. Candes, T. Strohmer and V. Voroninski, “PhaseLift: exact and stable signal recovery from magnitude measurements via convex programming,” Comm. Pure Appl. Math. 66, 1241–1274 (2013). [29] H. M. Ozaktas, S. Yuksel and M. A. Kutay, “Linear algebraic theory of partial coherence: discrete fields and measures of partial coherence,” J. Opt. Soc. Am. A 19, 1563–1571 (2002). [30] P. Thibault and A. Menzel, “Reconstructing state mixtures from diffraction measurements,” Nature 494, 68–71 (2013). [31] E. J. Candes, X. Li and M. Soltanolkotabi, “Phase retrieval via Wirtinger flow: theory and algorithms,” IEEE Trans. Inf. Theory 61, 19852007 (2014). [32] M. Schmidt, MicFunc. http://www.cs.ubc.ca/˜schmidtm/Software/minFunc.html. Minfunc, Accessed Nov 2014 (2005) [33] R. Horstmeyer, X. Ou, G. Zheng and C. Yang, “Overlapped Fourier coding for optical aberration removal,” Opt. Express 22, 24062–24080 (2014). [34] S. J. Haigh, H. Sawada and A. I. Kirkland, “Atomic structure imaging beyond conventional resolution limits in the transmission electron microscope,” Phys. Rev. Lett 103, 126101 (2009). [35] X. Ou, G. Zheng and C. Yang, “Embedded pupil function recovery for Fourier ptychographic microscopy,” Opt. Express 22, 4960–4972 (2014). [36] X. Ou, R. Horstmeyer, G. Zheng and C. Yang, “High numerical aperture Fourier ptychography: principle, implementation and characterization,” Opt. Express 23, 3472–3491 (2015). [37] World Health Organization, Malaria microscopy quality assurance manual, version 1. World Health Organization (2009). [38] D. L. van de Hoef et al., “Plasmodium falciparum-deprived uric acid precipitates incude maturation of dendritic cells,” PLoS One 8(2). 55584 (2013). [39] S. Becker, E. J. Candes and M. Grant, “Templates for convex cone problems with applications to sparse signal recovery,” Technical report, Department of Statistics, Stanford University (2010). 95 Chapter 7 Fourier ptychography in a conventional camera So far, our various investigations with Fourier ptychography have used a microscope. In this chapter, we implement Fourier ptychography within a conventional camera. This shift from a microscope to a conventional camera requires both new hardware and post-processing software. This alternative computational imaging technique is termed “overlapped Fourier coding” [1]. We use it to estimate and remove the intrinsic aberrations and misalignments within an imaging system. 7.1 Introduction and Background As we have seen, the space-bandwidth product of a standard microscope is limited to approximately 10 megapixels. Through the use of active illumination, Fourier ptychography can improve this limit to hundreds of megapixels. Besides microscopes, standard cameras, ranging from cell phone-sized cameras to large satellite based imagers, are also limited in the total number of features they can resolve. This physical resolution limit scales with the dimensions of the lens, but is often on the order of 10-100 megapixels for most applications. As detailed in [2], the unavoidable influence of lens aberrations is the root cause of this upper limit. Often, additional optical elements are used to correct for the aberrations that limit a camera’s space-bandwidth product (SP). These extra elements simultaneously increase the size, weight, and complexity of the physical setup. In this chapter, we explore an alternative approach to improving a camera’s resolution performance. Instead of adding optical components, we instead use a computationally based capture and post-processing scheme to remove aberrations. This procedure, termed overlapped Fourier coding (OFC), relaxes the complexity of the optics and digital sensor while maintaining a large SBP. Our correction scheme also acquires the phase of the optical field exiting the sample plane, which may in turn be used to digitally refocus the sample post-capture to ensure all image segments are in sharp focus. 96 (a)   Detector! SLM! Sample! (b) SLM aperture! c   1! 2! l f! f! f! (b) Aberrated Images Ij (x,y )! (c) High Res. Output sq(x,y )! ! 2! 1! Intensity (AU)! 0! Amplitude! 1! n 2! 2! π! Phase map! f! …   n L! Image I (x)! S (kx)! Wj(kx)! s (x )! …   System aberrations Aq(kx,ky)! π! -π! ! -π! Figure 7.1: Outline of the OFC procedure. (a) We place a transmissive SLM in the Fourier plane of a 4f system to digitally create different sub-aperture functions. (b) We capture a sequence of aberrated images while the SLM displays a sub-aperture shifted to a unique location between each snapshot. (c) We computationally transform the captured image set into a high-resolution amplitude and phase map, as well as an estimation of the camera’s low-order aberrations (figure adapted from [1]). 97 Just like Fourier ptychography, OFC first acquires a sequence of images. Instead of using active illumination via an LED array, it modulates the camera’s Fourier plane with a shifting aperture function between each snapshot. Here, the shifting aperture function we use is a small square pupil generated on an amplitude-modulating spatial light modulator (SLM), although other aperture function types, including phase-only modulation, may benefit alternative applications. It is important that the shifted sub-apertures spatially overlap with one another in the pupil plane by a certain degree. Second, OFC uses a unique algorithm to digitally synthesize the sub-aperture images into a high-resolution estimate of the complex optical field at the image plane. This algorithm extends prior work with FP [32], so it is also directly related to the well-known methods of phase-retrieval [4, 5] and ptychography [6, 52]. Unlike other holographic [8, 9] or non-holographic [10–12, 32] resolution-enhancing optical setups that measure phase, this work does not use external illumination. So while OFC cannot extend an optical system’s resolution beyond its aperture-defined cutoff, it can increase its SBP via removal of undesired aberrations and misalignments. This sets our goal as distinct from prior methods employing SLM modulation to acquire phase [13–16] or the related incoherent light field [17, 18]. Since the sample’s complex field is simultaneously acquired, removal of aberrative effects is not an ill-posed intensity PSF deconvolution as employed in fluorescent [19] or light-field [6] microscopy. Our procedure relies upon simulated annealing to iteratively reduce discrepancies between the image’s digital Fourier representation and the camera system’s physical Fourier plane. Annealing is a well-studied optimization procedure [21] that has been previously applied within the area of Xray ptychography to correct for an optical probe’s unknown shifted location [22, 23]. The proposed algorithm is closely related to this prior work, but instead searches over either a set of unknown Zernike modes, or a space of Fourier plane misalignments, or both, to increase a camera’s SBP. While many previous systems [24–26] can pre-calibrate for optical aberrations for later removal (assuming they also obtain a measure of the sample’s quantitative phase), OFC can continually update an improved estimate of a system’s imperfections from each image it captures. This may prove especially important in microscopy applications where objective lenses are often shifted or replaced, or in conventional cameras where aberrations are a function of object depth and zoom lens position, for example. In a broader sense, this situates the OFC scheme close to the realm of adaptive optics, which utilizes correction schemes targeting and removing optical distortions that can change with each acquired image. The remainder of this chapter is outlined as follows. First, we present the OFC imaging strategy and discuss its recovery algorithm in the context of an aberration-free setup. Second, we introduce the experimentally realistic situation of a camera containing low-order Zernike aberrations. We also demonstrate how OFC uses simulated annealing to estimate and remove these aberrations from its final reconstruction. Third, we show how the simulated annealing approach may additionally 98 account for aberrations that do not fall within the typical Zernike model, such as distortion and lens misalignment. Finally, we experimentally demonstrate OFC with an SLM-based shifting subaperture to recover a sample’s image and phase with resolution at the camera’s maximum cutoff spatial frequency. Our annealing algorithm likewise helps remove included aberrations to increase the setup’s SBP. 7.2 Principle of operation In this section, we first develop a mathematical model of OFC’s imaging routine. Then, we detail OFC’s reconstruction algorithm to recover a sample’s amplitude and phase at high-resolution from a series of low-resolution measurements, assuming aberration-free optics. 7.2.1 Image acquisition As diagrammed in Fig. 8.1, OFC is implemented in a conventional 4f imaging system with a simple modification: a transmissive amplitude SLM inserted into the camera’s Fourier plane (i.e., aperture plane). A 4f imaging system is not a requirement for effective operation – placing the SLM at any plane conjugate to the sample plane (i.e., the aperture plane of a compound lens system, or the back-focal plane of a microscope objective) will function similarly. We will assume the sample s(x, y) at the object plane is illuminated with spatially coherent, quasi-monochromatic light. OFC directly extends to operate with partially coherent illumination. We also assume that our imaging setup contains pixels that are matched to the maximum cutoff frequency of its unmodulated Fourier plane (i.e., each pixel matches the system’s minimum PSF width). This enables us to effectively verify the accuracy of our reconstruction through a direct comparison with an unmodulated image. A more practical OFC setup should match the pixel size to the NA of a single sub-aperture image, which will lead to a final reconstructed resolution that surpasses the detector pixel count. The square SLM of width L is configured to display a small square sub-aperture (i.e., optically transparent area) of width and height `. This sub-aperture is shifted to n2 different equally spaced locations arranged on a 2D rectilinear grid. The j th sub-aperture will be displaced from the origin by a two-dimensional vector cj = (cxj , cyj ) for 1 ≤ j < n2 . By setting n > L/`, we ensure that each sub-aperture overlaps partially with its neighboring sub-apertures. In practice, it is useful to select n such that the sub-aperture images overlap by approximately 70%, following a similar parameter optimized for ptychographic imaging [50]. All simulations and experiments in this work use n = 9, requiring a total of n2 = 81 images per capture sequence. We set L/` = 2.5 to ensure each sub-aperture window overlaps with its neighbors by 75%. An ideal binary amplitude SLM will completely block light when its pixels are switched to opaque and pass 100% of any incident light when switched to transparent. In a practical setup, a 99 (b) Window! Wj-1   Wj   SŸWj   ky! (d) Intensity Measurements Ij(x,y)! (c) Constraint! F I1   …   S'ŸWj   y! pj(x,y)   pj'(x,y)   F -1! Cutoff! 1! (a) Spectrum estimate S0(kx, ky)! Eq.(3)! x! Ij   kx! 0! …   Figure 7.2: One stage of the basic OFC algorithm. For each window position Wj , a segment of the spectrum estimate Ŝ0 (kx , ky ) is extracted (corresponding to the j th position of the SLM’s shifting sub-aperture). In the spatial domain, the amplitude associated with this windowed spectrum is constrained with the measured image Ij (x, y) to form p0 (x, y). Example measured images are in (d). The result is Fourier transformed back to the aperture plane, where it is used to update the spectrum estimate, Ŝ0 (kx , ky ). Figure adapted from [1]. transmissive SLM will have a finite optical density b, which we include in our model by defining the j th SLM sub-aperture’s transmission function Wj as, Wj (kx − cxj , ky − cyj ) =   1, |kx | ≤ ` and |ky | ≤ ` 2 2  b, (7.1) otherwise, where (kx , ky ) represent spatial frequencies, the spatial coordinates at the aperture plane. Note that we use (kx , ky ) here to label to aperture plane, instead of (x0 , y 0 ) as in previous chapters, to distinguish this system from the microscope case. Eq. (7.1) states the SLM’s modulation is a biased rect function of width ` and center cj in the Fourier plane. Modeling the optical field emerging from the sample surface as s(x, y), we can write the complex field directly before the aperture (i.e., SLM) plane as F[s(x, y)] = Ŝ(kx , ky ), where F is the Fourier transform operation, and we will refer to Ŝ as the sample’s spectrum. After passing through the j th sub-aperture window, Wj , the exiting field is the product Ŝ(kx , ky )Wj (kx , ky ) assuming the SLM is thin. This modulated field then propagates to a detector at the image plane, where to a first-order approximation the j th image measures Ij (x, y) = |F[Wj (kx , ky )Ŝ(kx , ky ))]|2 (7.2) for 1 ≤ j ≤ n2 . The first goal of OFC post-processing is to recover a high-resolution complex estimate of the sample S from the above set of n2 intensity measurements, before attempting aberration estimation and removal. 100 7.2.2 Aberration-free OFC reconstruction The basic OFC post-processing algorithm employs the well-known strategy of alternating projections, similar to FP and other prior phase retrieval methods [4, 5]. Our measured images constrain our sample estimate to a known set of amplitudes in the spatial domain, while our shifting SLM sub-aperture uniquely constrains its support in the Fourier domain. Assuming an aberration-free setup, the OFC algorithm outputs a reconstruction of the complex field exiting the sample plane at the employed imaging system’s maximum resolution. Unlike prior formulations of the related FP procedure in previous chapters [32], the sample here may be thick, optically reflective, and under arbitrary quasi-coherent illumination. The aperture scanning procedure detailed in [16] also follows similar steps to both FP and OFC, and the interested reader should refer to this related work for a comparison of aperture scanning versus full camera translation for ptychographic resolution improvement. OFC’s SLM modulator requires a modified algorithm in comparison to these other techniques. The following explains one iteration of the OFC process. Before initiating the loop, we generate a first guess of a high-resolution complex sample spectrum, Ŝ0 (kx , ky ), as shown in Fig. 7.2(a). A good starting point is an up-sampled version of the low-resolution image recorded with a centered scanning aperture, modulated by a randomly assigned phase. First, we compute a masked spectrum estimate by modulating Ŝ0 with one of our SLM’s sub-apertures, from Eq. (7.1). Starting with j = 0 and continuing until j = n2 − 1, we compute the spectrum product P̂j (kx , ky ) = Ŝ0 (kx , ky )Wj (kx , ky ). For an SLM with perfect contrast, this corresponds to selecting a windowed region of Ŝ0 , as in Fig. 7.2(b). Second, we compute the Fourier transform of this spectrum product to simulate light’s propagation to the image plane: pj (x, y) = F [Pj ((kx , ky )]. Third, we replace the amplitude of the p resulting Fourier transform pj with the known amplitude from the j th detected image Ij to form p0j : p0j (x, y) = q Ij (x, y) pj (x, y) . |pj (x, y)| (7.3) Like other phase retrieval strategies, this leaves the estimated field’s phase unchanged. Fourth, we   inverse-Fourier transform p0j to create a new spectrum: P̂j0 (kx , ky ) = F −1 p0j (x, y) . Finally, we update our sample spectrum estimate with P̂j0 . While we accounted for the SLM’s imperfect modulation when extracting P̂j from Ŝ0 , we do not use the same support function when re-inserting P̂j0 in this update. Instead, we only update the aperture area that in an ideal case is optically transparent, leaving all other areas unchanged: Ŝ0 (kx − cjx , ky − cjy ) =   P̂ 0 (kx − cj , ky − cj ), |kx | ≤ ` and |ky | ≤ ` x y j 2 2  Ŝ0 (kx − cj , ky − cj ), x y otherwise. (7.4) 101 Target  phase        U(x,y)! ! π! -π! Reconstructed  |sq(x,y)|   Reconstructed        sq(x,y)   Phase! Amplitude! π! -π! Reconstruction error Eq vs. iteration! Per pixel NMSE Eq! Phase! Amplitude! Target  amplitude  |U(x,y)|   Iteration # q! Figure 7.3: OFC Simulation for an aberration-free noisy imaging system. Both the sample amplitude and phase are well-reconstructed. Reconstruction error grows linearly with system noise while the OFC algorithm continues to converge (figure adapted from [1]). This selective support constraint is unique to OFC. However, it shares close parallels with the methods of hybrid input-output [4] and iterated projections [5] that do not strictly enforce an optical system’s support constraint every iteration. While alternative update strategies are certainly possible, we empirically determined that Eq. (7.4) leads to quick and accurate algorithm convergence. After updating all n2 overlapping sub-aperture areas of Ŝ0 , we proceed to a second iteration by advancing our spectrum estimate: Ŝ1 ← Ŝ0 and resetting j = 0. After q iterations, we inverse Fourier transform the final recovered spectrum Ŝq (kx , ky ) to recover our high-resolution sample estimate, sq (x, y). 7.2.3 Aberration-free OFC simulation Fig. 7.3 demonstrates the OFC algorithm’s successful convergence to an accurate amplitude and phase image in simulation. With an ideal camera (no aberrations or misalignments, but including Gaussian noise), recovery at the system’s native resolution proceeds directly following the five update steps in Section 7.2.2. Our ideal simulated 4f camera has a maximum F-number of 7.5 (PSF width= 5µm width at the detector assuming λ = 632 nm illumination), and its detector pixels are 5 µm to match this unmodulated PSF. We shift a square sub-aperture that is ` = 4 mm wide across a total distance of L = 10 mm, where each shift is ∆c = 1 mm along one dimension. In two dimensions, this leads to a total of 81 images and an aperture overlap percentage of 75%. Simulated intensity images are N = 10002 pixels and have a 12-bit well depth. These setup parameters closely match our experimental setup (see Section 7.5) and are used for all subsequent simulations. The unknown complex sample U (x, y) in Fig. 7.3 is an Air Force resolution target with a 2.5 µm minimum feature size (i.e., half a pixel) with a multiplied cubic phase function: ∠U (x, y) = α(x3 + y 3 ), where ∠ indicates phase and α = 20π. We add 2% random Gaussian noise to each subimage after detection of its intensity to simulate detector-induced noise. Note that the reconstructed 102 image is at the full resolution of aberration-free imaging system – approximately three times larger than the resolution of each single sub-aperture image (e.g., compare Fig. 7.3 and Fig. 7.2(d)). This resolution increase follows expectations, as our final synthesized aperture is roughly three times wider than a single sub-aperture. The reconstructed phase closely matches the target phase but is impacted in quality more than the target amplitude, especially at lower spatial frequencies, and is offset by a constant global phase shift. The normalized mean-squared error (NMSE) Eq between the amplitude of the target sample U (x, y) and reconstruction sq (x, y) offers a useful metric to examine algorithm performance: Eq = sX X x (|sq (x, y)| − |U (x, y)|)2 /N. (7.5) y This metric is also used in [4]. Here N is the total number of image pixels and q indicates iteration number. In simulation, the original input amplitude |U (x, y)| is known a-priori. Fig. 7.3 plots the NMSE as a function of 15 iterations, where in one iteration we update all 9 x 9 overlapping sub-aperture regions. Here, each curve represents a different amount of zero-mean Gaussian noise: no noise, with variance σ 2 = 0.5%, and with variance σ 2 = 2% of the maximum signal value. As expected, the algorithm is error-reducing and convergence accuracy decreases with increased noise. A perfect optical system is never experimentally realizable. Inherent aberrations eventually limit the performance of all lenses. Next, we modify the OFC approach to account for and computationally remove aberrations. First, we show how to correct for an imaging setup’s unknown low-order aberrations. Second, we demonstrate removal of geometric distortions. Both of these steps increase the SBP of simple imaging lenses. 7.3 OFC with simulated annealing By repeatedly taking pictures through an overlapping sub-aperture, OFC captures slightly redundant data. This redundancy not only allows us to accurately extract sample phase, but also helps us compute a low number of unknown variables that influence the imaging process. A well-known method of searching over a space of unknown variables, to maximize or minimize a particular function of merit (i.e., error function), is termed simulated annealing (SA). Instead of exhaustively searching through all possible unknown variable configurations, SA takes an iterative approach. At one iteration, annealing first randomly searches through a small number of different configurations and selects the configuration that minimizes its error function. Then, it uses this configuration as starting point for the search during the next iteration. As iteration continues, the algorithm slowly reduces the range over which it randomly searches for error-minimizing states. For many problems, this type of iterative local search is very efficient at seeking out global minima of nonlinear functions [21]. 103 Simulated  Annealing   Simple aberration A(kx, ky)! π! Intensity Measurements! Spectrum Sj(kx, ky)! Wj   1! SjŸWjŸAj! F ψj   I1   Ÿ   -π! 0! Sj'ŸWjŸAj! Eq. (11)!   ψj'   F -1! …   Update Sj   Ij   Update Aj! …   Figure 7.4: Schematic of the OFC algorithm with simulated annealing (SA-OFC). We use the same PR loop outlined in Fig. 7.2 with three additions. First, an estimate of the systems aberrations at the Fourier plane Aj is now multiplied with the spectrum estimate Ŝ at each PR loop. Second, SA is used to compare T different perturbed versions of the j th windowed spectrum with the j th intensity measurement. Third, the error-minimizing aberration perturbation Atjmin and corresponding sample estimate ψjtmin are calculated through Eq. (7.9)–Eq. (7.10), which then update Aj+1 and Ŝj+1 via Eq. (7.11). Figure adapted from [1]. OFC’s space of unknown variables includes any aberrations or misalignments within the optical system. At each sub-aperture location, our function to minimize is the mean-squared error difference between the recorded intensity image, Ij (x, y), and the corresponding image that would result from detecting our complex windowed spectrum estimate, P̂ (kx , ky ). For each iteration, we make several guesses about how the optical system might be aberrated or misaligned, compute the resulting image captured through each imperfect camera, and then select the imperfections that yield the closest image to our observed data as the starting point for the next iteration. As we will demonstrate, this process is both effective at recovering the correct imperfections and robustly removing these imperfections from a final sample solution. However, we must assume a-priori knowledge about which subspace of imperfections to search through, as search time will scale linearly with the dimensions of this search space. 7.3.1 Characterization and removal of low-order aberrations We may account for the wavefront-based aberrations of our camera system using a multiplicative phase function, A(kx , ky ), at its Fourier plane. It is common to decompose A into a sum of weighted P Zernike polynomials on the unit disk: A(ρ, θ) = exp (2πi l al Wl (ρ, θ)), where Wl are the orthogonal Zernike polynomials and al are the associated weights. We will begin our explanation by only considering the most significant single Zernike phase aberration, defocus, which in Cartesian coordinates 104 takes the form, 2 2 A(kx , ky ) = e(id(kx +ky )) , (7.6) where d is a defocus aberration weight. Selecting Eq. (7.6) here as our aberration model calls upon prior knowledge that this simple example system is primarily susceptible to defocus. In the following section, we discuss how to correct for multiple higher-order Zernike aberrations. In addition to an estimated spectrum Ŝ0 , simulated annealing OFC (SA-OFC) also initializes an estimated aberration map A0 , as shown in Fig. 7.4. In all tests here, we initialize with A0 (kx , ky ) = 1. SA-OFC begins with j = 1 and computes T different candidate aberration functions at the j th loop: 2 2 Atj (kx , ky ) = Aj (kx , ky )e(i∆t (kx +ky )) (7.7) for 1 ≤ t ≤ T . Here, ∆t is a number selected randomly from a uniform distribution on (−ra , ra ), where ra is the annealing search radius. Each candidate Atj can be thought of as a random perturbation of the current iteration’s estimated aberration function, which is restricted to a limited search distance, ra . Next, we compute T candidate aberrated spectra, Ψ̂tj = Ŝj Atj . We then window these spectra with the sub-aperture function Wj and Fourier transform the result to form a set of T uniquely aberrated low-resolution image fields: ψjt (x, y) = F[Wj (kx , ky )Ψ̂tj (kx , ky )], (7.8) The annealing ends with identification of the candidate aberration perturbation Atjmin (kx , ky ) that minimizes the MSE between the set of candidate images ψjt (x, y) and our measured image through the j th sub-aperture, Ij (x, y): tmin = arg min t XX x |ψjt (x, y)| − q 2 ! Ij (x, y) . (7.9) y Note that since we only measure intensity, the merit function in Eq. (7.9) only considers the amplitude of each ψjt . We use the error-minimizing aberration perturbation, Atjmin , as our annealing search’s starting point for the next (j + 1)th sub-aperture image: Aj+1 (kx , ky ) = Atjmin (kx , ky ). (7.10) We also re-insert Atjmin into Eq. (7.8) to find ψjtmin , our optimal aberrated image field estimate. p ψtmin 0 We then constrain ψjtmin with our measured intensities following Eq. (7.3): ψj = Ij jtmin . After 0 0 |ψj | Fourier transforming ψj into Ψj , we are then ready to update our unaberrated sample spectrum estimate, Ŝj . To remove the effects of aberrations, we adopt a strategy common to prior algorithms 105 (a) OFC (no annealing)! A! ϕ! π! Ground truth aberration A(kx, ky)! π! -π! -π! (b) SA-OFC (with annealing)! A! ϕ! π! Recovered aberration Aq(kx, ky)! π! -π! -π! (c) Example detected image: centered aperture! (d) Example detected image: diagonal edge! Figure 7.5: Simulation results of the SA-OFC algorithm. (a) Without simulated annealing, OFC cannot combine a set of aberrated sub-images into an accurate full-resolution complex field estimate. The induced defocus aberration is show to right. (b) The annealing OFC algorithm accurately recovers both the sample’s amplitude and phase and the 4f setup’s aberration map Aq (kx , ky ). Here, q = 20. (c)-(d) Example low-resolution sub-aperture images, for comparison (figure adapted from [1]). like ePIE [52] and EPRY [30] and effectively divide out the aberration function estimate Aj from 0 Ψj : Ŝj+1 (0) = Ŝj (0) + A∗j 0 (Ψ̂ − Ψ̂j ). |Aj | j (7.11) Here, we have modified our notation to include the iteration number in parenthesis and indicate update by the j th sub-image with a subscript. This random search and update process is repeated for all 0 ≤ j ≤ n2 −1 sub-apertures to complete one iteration. To encourage the SA-OFC algorithm’s convergence, we linearly reduce the search radius ra at the completion of each iteration: ra ← ra −αa , where typically αa ≈ ra /q. Iteration continues for q loops to form our final camera aberration map, A(q), and sample spectrum solution, Ŝ(q). Fig. 7.5 presents an example simulation of the SA-OFC algorithm using the same 4f setup from Section 7.2. The target sample matches Fig. 7.3(a)’s, but now with α = 5 defining its cubic phase envelope. We add defocus aberration by setting d = 200 in Eq. (7.6) and multiplying the resulting A(kx , ky ) with each sub-aperture spectrum during image simulation. To first test the effect of aberrations without simulated annealing, we run the basic OFC algorithm from Section 7.2 to recover the amplitude and phase maps shown in Fig. 7.5(a). Because the algorithm incorrectly updates our sample estimate with aberrated low-resolution intensity images, the final solution does 106 not closely resemble the target sample. Switching to the SA-OFC algorithm in Section 7.3.1 enables simultaneous aberration map estimation and removal, as shown in Fig. 7.5(b). For this example, we use the same initial conditions and noise as Section 7.2, setting T = 8, q = 20, ra = d/2, and the initial aberration estimate to a constant function: A0 (x, y) = 1. The recovered aberration map’s defocus coefficient, dq , differs from ground truth by ∆d = (dq − d)2 = 25.71. Annealing generally performs better with more candidate search functions (i.e., a larger value of T ), but will require linearly more computation [22]. The algorithm is also sensitive to the selected search radius ra . While all simulations here used ra = d/2 and linearly decreased it to zero by the q th iteration by setting αa = ra /(q + 1), it may prove useful to additionally optimize over this free parameter. We plot the SA-OFC algorithm’s error reduction versus iteration number, using the error metric from Eq. (7.5), in Fig. 7.6(a). Here, we compare simulated annealing (blue curve) with two alternative post-processing options. All parameters here match those for the data in Fig. 7.5 (here now each curve is averaged over 5 runs with σ 2 = 2% noise). First, no annealing (“no SA”, red) leads to a significantly higher sample recovery error (see Fig. 7.5(a)) but continues to offer errordecreasing performance, demonstrating algorithm stability. Second, we assume a-priori knowledge of the 4f setup’s aberration map A0 (kx , ky ) (“known”, green) to recover an almost exact solution. Here, instead of following Eq. (7.10) to estimate the aberration function each iteration, we simply set Aj+1 (kx , ky ) = A0 (kx , ky ), the actual aberration map in Eq. (7.6), for all j. The SA-OFC algorithm’s error falls somewhere between. While recovery is not perfect, there is certainly enough redundancy within the captured dataset to significantly improve an image’s SBP via removal of the negative defocus aberrations. In Fig. 7.6(b), we again repeat Fig. 7.5’s simulation, but now vary the amount of defocus aberration d in the ground-truth aberration, A(kx , ky ). Again, each plot point is an average over five independent tests with σ 2 = 2% noise. As expected, SA-OFC remains bounded below by the case of knowing and directly accounting for a complex aberration map within this redundant dataset. However, it achieves much lower error than an annealing-free algorithm, even for significantly defocused image sets. 7.3.2 Characterization and removal of geometric distortion A number of camera imperfections that negatively impact image fidelity cannot be summarized as a phase-only modification to the Fourier plane (i.e., a Zernike aberration). Examples of such imperfections include unknown changes in magnification, image distortion, and vignetting caused by system misalignment. In this section, we explore how to computationally correct for these undesirable effects. We focus our attention on one specific form of misalignment that significantly impacts our experiments–the displacement of the Fourier plane from its assumed location, which is directly connected to image distortion. We outline how the SA-OFC algorithm accounts for this 107 (a) Error Eq vs. iteration q, comparison! (b) Error Eq vs. aberration size d! x 10-4! Error Eq (AU)! ! Error Eq (AU)! 7! 50! 6! 30! 5! 10! 50! Iteration # q! 200! 100! 300! Aberration estimate error Δd! 70! 8! 400! Aberration size d! Figure 7.6: (a) SA-OFC is error reducing (blue) and exhibits much lower error than the regular OFC algorithm (red) in the presence of aberrations (assuming defocus d = 200). The ideal case of performing SA-OFC with an a-priori known aberration, initialized and enforced each iteration, is plotted in green. (b) As the aberration size increases, SA-OFC, OFC, and the ideal case all slowly decrease in MSE performance. SA-OFC exhibits an aberration recovery mean-square error ∆d scaling roughly as 10-15% of d, which may be improved with additional fine-tuning of annealing parameters (figure adapted from [1]). type of unknown parameter by updating the estimated size and location of each SLM sub-aperture, leading to an improved-resolution complex sample reconstruction. It is challenging to construct an imaging system with a perfectly flat, centered Fourier conjugate plane (i.e., aperture plane or back focal plane). Slight curvature across the plane is often encountered, even with the aid of advanced lens design software [31]. Element misalignment during system assembly will additionally shift the Fourier plane away from its ideal location. Misalignments may point perpendicular to the optical axis, which will cause the Fourier plane to laterally shift and add a linear phase ramp across the image plane. Or, they may point axially, which will primarily cause the Fourier plane to scale in magnification. Depending upon the size and linearity of this scaling, the field at the image plane can become magnified, geometrically distorted, and also defocused. The SA-OFC algorithm can help measure and remove this type of misalignment. First, let us assume that the location of each SLM pixel is known a-priori and we can accurately measure the optical field associated with one sub-aperture image, Ij (x, y). This field’s digital inverse Fourier transform should create a masked spectrum with a clearly visible window function, Wj (kx , ky ), modulating the sample spectrum via Eq. (7.2). Any deviation in the position and size of the computed window Wj (kx , ky ) from the known SLM pixel map will inform us of how the optical setup’s Fourier transformation differs from an exact digital computation. We may computationally account for any such deviation (i.e., system misalignment) by digitally updating the window’s assumed position and size with a more accurate reflection of the imperfect optical setup. In practice, although we assume accurate knowledge of each sub-aperture’s centered location
cxj , cyj and size `j when we extract and insert the updated spectra (e.g., Eq. (7.4)), these variables are not known a-priori. Moreover, each sub-aperture image only measures the amplitude of the 108 (a) Geometric distortion annealing! (b) Misalignment correction! ky! ky! Error Eq (AU)! Wj! (c) Error Eq vs. iteration! Random shifting! kx! kx! Wj(kx, ky)! start! δ   Wj(kx-δx, ky-δy)! Iteration # q! 10th iteration! ! actual! Figure 7.7: SA-OFC removes geometric distortions from within a camera. (a) SLM sub-aperture centers at the Fourier plane (blue dots) may be distorted by optical misalignments to unknown positions. During algorithm iteration, each estimated center is randomly perturbed by vector (δx , δy ) as we generate and compare trial images to measured data. (b) Ground-truth sub-aperture centers from the simulation in Fig. 7.5 are here radially offset, which simulates a Fourier plane displaced by 100 µm axially. The actual center of each displaced sub-aperture is marked with an ‘x’. Initializing each sub-aperture center on an assumed rectilinear grid (blue dots), the simulated annealing process draws these estimates close to the actual centers after 10 iterations. (c) Error decreases with iteration, similar to Fig. 7.6(a). Figure adapted from [1]. complex field, thus preventing us from adopting the above simple misalignment correction scheme. However, just like the determination of the unknown aberration function A(kx , ky ) in the last section, we can determine the 3 · n representative variables cjx , cjy and `j from our intensity measurements using a simulated annealing algorithm. Our “geometric distortion” annealing process proceeds as follows. At the j th iteration, we construct T candidate window functions Wjt , defined for 1 ≤ t ≤ T as, Wjt (kx − cjx − δxt , ky − cjy − δyt ) =  t t  1, |kx | ≤ δa ` and |ky | ≤ δa ` 2 2  b, (7.12) otherwise. Here, (δxt , δyt ) are random perturbations selected from the uniform distribution [−rx , rx ] that modify each window function’s center, and δat is a random perturbation from the uniform distribution [−rw , rw ] that randomly scales the aperture size. Just like ra , rx and rw are annealing search radii that we reduce by an α factor each iteration, again using αx = rx /(q + 1) and αw = rw /(q + 1) in all demonstrations. Next, we create T different windowed spectra, each of which we Fourier transform into a simulated image: ψjt (x, y) = F[Wjt (kx , ky )Ŝj (kx , ky )], (7.13) As with the aberration annealing process, we compare each simulated image ψjt with the j th measured image Ij to find the error-minimizing perturbed window function, following Eq. (7.9). We save a map of all updated centers, which serve as the starting window positions and sizes during the next 109 cycle through all n windowed images (i.e., as the algorithm proceeds from the 1st to qth iteration). Fig. 7.7 presents a simulation of the geometric distortion SA-OFC process. Here, we again use our familiar noisy optical system, now free of Zernike aberrations but geometrically distorted in a radial manner. Each sub-aperture center is displaced from its original location by δcj = 0.05cj , suggesting a Fourier plane that is axially offset from its true location by 250 µm. Fig. 7.7(b) illustrates this distortion, where the blue circles denote the distortion-free rectilinear centers cj (i.e., the algorithm’s assumed centers) and the x’s mark their actual locations, cj + δcj . Iterating q = 20 times accurately identifies the unknown misaligned sub-aperture centers. Using parameters T = 8, rx = 10 pixels and rw = 0 pixels, we reduce an initial average center offset from 125 µm to 14 µm. Removal of these geometric effects qualitatively improves the final image similar to Fig. 7.5. As with using SA to correct for misalignments in ptychography [22], we also empirically find that holding (δxt , δYt ) to zero for the first few iterations helps improve convergence. The simulation in Fig. 7.7 holds the search radius to 0 for the first 5 iterations, as the blue error curve for SA indicates in Fig. 7.7(c). As with our Zernike aberration example’s error in Fig. 7.6, this curve’s final error is much lower than OFC without simulated annealing (red) and is bounded below by running OFC with the correct geometric distortions known a-priori (green). 7.3.3 The complete OFC algorithm We may also simultaneously determine and remove low-order Zernike aberrations and geometric misalignments in one SA loop. By searching over both candidate aberrations with Eq. (7.7) and sub-aperture positions with Eq. (7.12) in a parallel manner, we may use the simulated annealing process summarized in Eq. (7.7) – (7.11) to search over a larger sub-space of unknown camera configurations. However, as we increase the dimension of this parameter search space, computation time will grow exponentially. For example, to search over both T different possible defocus settings and T different sub-aperture positions with the same amount of coverage, each loop must now test T 2 candidate states. If we would also like to search for the effects of x and y astigmatism, the number of candidate tests jumps to T 4 . For a p-dimensional parameter search space, the number of candidate tests becomes T p . As our constructed optical setup was primarily influenced by subaperture position shift and low-order aberrations, this exponential scaling does not become a major concern, as we experimentally demonstrate in the following section. However, for systems suffering from many possible equally-weighted aberrations and geometric imperfections, alternative iterative strategies, such as conjugate gradient descent [23], will prove more efficient than simulated annealing. 110 (a) Experimental sub-image! (b) OFC experiment! (c) SA-OFC experiment! Intensity (AU)! 1! 0! (d) Experimental traces! 0! 50  µm   50  µm   pixel! 150μm! 50  µm   Figure 7.8: Experimental results from an OFC setup. We first image an Air Force resolution target to test the algorithm’s spatial resolution performance. (a) A single sub-aperture image exhibits low spatial resolution. (b) The OFC algorithm without annealing recovers a sharper image, but still contains artifacts. (c) The SA-OFC algorithm further improves the output field sq (x, y)’s spatial resolution (see text for algorithm parameters), as highlighted by traces in (d). Figure adapted from [1]. 7.4 Experimental results Our experimental OFC setup closely matches the optical conditions used for the simulations thus far. We created a 4f camera using two bi-convex lenses (diameter= 25.4 mm, f= 75 mm, Thorlabs AC254075). Without any additional corrective elements, it is easy to confirm that such a single-element large F-number lens (˜F/3) exhibits significant off-axis aberrations. We use an amplitude-modulating SLM (Epson HDTV LCD, BBS Bildsysteme 1920×1080 pixels placed between two crossed polarizers) to create the Fourier plane’s shifting sub-aperture masks. The SLM pixel size is 23 µm. Only the central 1202 pixels of the SLM are varied, defining the camera’s full aperture width as 2.76 mm. Outside of this finite aperture range, significant aberrations from each singlet lens proved challenging to correct for accurately and were thus blocked. At the image plane, we place a CMOS detector (Prosilica GX-1920 with 4.54 µm pixels). We illuminate our sample with a collimated 632 nm laser beam (no spatial filtering used). Unlike FP, the illumination field’s specific shape and phase is not critical and can remain unknown. The spatial coherence length of the illuminating beam must be as wide as the sample at the object plane, and it should be relatively narrow-band (≤ 20 nm spectral bandwidth [16, 28]). This may be achieved with an LED placed sufficient far behind the sample. Alternative algorithms exist to additionally incorporate the effects of a partially coherent source to further reduce these illumination requirements [29]. 111 The sub-aperture displayed on the SLM is a square with side length ` = 0.92 mm, which we sequentially shift laterally n = 9 times in x and y across an L × L =2.76 mm× 2.76 mm square Fourier plane area. In each step, we shift the aperture by 0.28 mm, raster scanning it until it passes through all n2 = 81 unique aperture locations. This leads to a sub-aperture overlap of approximately 70% and an expected resolution gain of 3. To ensure that each of the 81 captured images is properly exposed, we take 3 snapshots at each sub-aperture position with an exposure sequence of .01, .1 and 1 second, and combine these three images via high-dynamic range (HDR) processing [32]. The HDR process may be omitted by using a high bit-depth detector or if larger sub-apertures are utilized. Finally, we experimentally measure the optical density of the SLM as 22 (i.e., the central sub-aperture region is 22 times as bright as the outer apodizing region), making b = 0.045 in Eq. (7.1). Fig. 7.8(a) displays an example sub-aperture image captured with the aperture centered on the optical axis (i.e., cj = 0). First, we attempt OFC reconstruction without simulated annealing as outlined in Section 7.2. The reconstructed image intensity after q = 20 iterations is in Fig. 7.8(b). Note the resolution is improved significantly as compared to that of a single sub-aperture image, but artifacts remain primarily due to a misaligned Fourier plane. To correct these misalignments, we then implement OFC with simulated annealing. We use Section 7.3.3’s “complete” SA-OFC algorithm to jointly search over aberrations and geometric misalignments. We adopt an aberration model that contains the first five primary aberrations beyond linear tilt in x and y – defocus, x/y   P8 astigmatism and x/y coma: A(ρ, θ) = exp 2πi l=4 al Wl (ρ, θ) . We apply Eq. (7.7) to search now over these 5 orthogonal aberration modes, in addition to Eq. (7.12) to correct geometric sub-aperture errors. This creates a total search space dimension of p = 6 (T 6 candidate tests). Setting the SA-OFC algorithm parameters to q = 20 iterations, T = 2, ra = 10, rx = 8, rw = 1, and with linearly reducing α’s, we obtain the intensity image in Fig. 7.8(c). Feature sharpness increases most notably along the vertical direction. Fig. 7.8(d) plots a line trace through each image’s group 4, element 6. Both the sub-aperture image (a) and uncorrected output (b) fail to resolve this feature (17.54 µm-width). The aberration-corrected image (expected spatial frequency cutoff at 17.2 µm−1 ) does resolve this group, confirming the ability to restore a potentially misaligned camera system close to its native resolution with SA-OFC. More advanced models, which account for increasingly complex optical aberrations and misalignments, may achieve a sharper resolution, but will require additional computation. The experimental system aberration map Aq (kx , ky ) and Fourier plane misalignments, simultaneously recovered with the image in Fig. 7.8(c), are shown in Fig. 7.9(a) and Fig. 7.9(b), respectively. We note the modal weights of the final aberration function decrease with mode number, as expected for most lens systems. Linear x and y-tilt aberrations are partially corrected for by the geometric Fourier plane realignments and are thus left out of Aq (kx , ky ) for computational efficiency. Fig. 7.9(b) 112 (a) Experimental aberration! ky! (c) Experimental error vs. iteration! (b) Experimental misalignment! 6! .02! kx! -π! ky! Weight! Avg. δcy! 10! x! x! .01! 8! 1! -6! Mode #! 6! 4! .015! .005! π! 0! Error Eq (AU)! Avg. δcx! -6! kx! start! 3! 6! 9! 12! Iteration # q! 10th iteration! ! Figure 7.9: (a) OFC-computed aberration map, recovered simultaneously with the image in Fig. 7.8(c). (b) Geometric misalignments of our 4f setup’s Fourier plane also simultaneously recovered by the SA-OFC algorithm. Average shifts for each row/column are plotted on the side/top. (c) Plot demonstrating the algorithm’s error reduction with iteration. After 10 iterations, SA-OFC’s NMSE is 3 times lower than the direct OFC algorithm without annealing (figure adapted from [1]). shows that geometric realignments are primarily a linear function along the vertical direction, as highlighted by its plots of average shift along kx and ky (average δcx and δcy ). This matches the significantly sharpened vertical image features in Fig. 7.8(c). Finally, we plot the algorithm’s NMSE, Eq , in Fig. 7.9(c) (calculated via Eq. (7.5)). The addition of annealing (SA-OFC) brings our image much closer to expected measured intensities. We also plot Eq after reducing the geometric misalignment search radius to rx = 4, offering an idea of algorithm sensitivity to this type of free parameter. Next, to demonstrate that the OFC procedure can also acquire quantitative phase, we image a monolayer of polystyrene microspheres (diameter = 117 µm, index nsphere = 1.594 at λ = 632 nm) coated on a microscope slide and immersed in oil (index noil = 1.5915). Using the same SA-OFC algorithm parameters outlined above, we converge upon a high-resolution sample amplitude and phase reconstruction. A phase map containing five microspheres is in Fig. 7.10(a). We take a line trace through one microsphere’s phase (dashed line) and plot its optical thickness h, computed from λ the measured phase map ∆φ(x) via the following equation: h = 2π ∆φ(nsphere − noil )−1 . This curve in Fig. 7.10(b) closely matches the shape and optical thickness of a perfect sphere. Fig. 7.10(d) includes the microsphere slide’s computed aberration map, indicating successful convergence following two observations. First, the aberration map exhibits a quite similar structure to that for the resolution target shown in Fig. 7.9(a), apart from a constant phase offset. Second, microspheres recovered without aberration correction in Fig. 7.9(e) do not closely match the expected thickness profile, thus pointing towards the necessity of an aberration and misalignment correction strategy to ensure ptychography-based recovery schemes like OFC remain quantitatively accurate. 113 (a) Microsphere phase image (b) Microsphere thickness Thickness (µm) 2π 0 100  µm   Width (µm) (c) Examples from Low-res. sequence cj = (.56,.56)! 1! cj = (1.12,1.12)! Intensity (AU)! cj = (0,0)! 0! π! ky! kx! -π! (e) No annealing (d) Aberrations 100  µm   π! 100  µm   -π! Figure 7.10: SA-OFC for quantitative phase. (a) Recovered phase map of 5 microspheres with aberration correction. (b) Line trace through one sphere demonstrates quantitative accuracy. (c) Example captured images. (d) Simultaneously recovered aberration map exhibits a similar structure as Fig. 7.9(a), as expected for the same optical setup. (e) Phase map recovered by the OFC algorithm without simulated annealing, where low-order aberrations clearly impact the fidelity of the reconstructed phase map (figure adapted from [1]). 114 7.5 Discussion and Conclusion We have demonstrated how a sequence of low-resolution images can be computationally processed into a full-resolution amplitude and phase image, while simultaneously extracting camera aberrations. While our experimental imaging setup offers proof-of-concept aberration removal from a series of simple lenses, several steps may help move the OFC concept towards a more practical device. First, the transmissive SLM used by this setup exhibits an imperfect pixel fill-factor, which diffracts light incident upon the Fourier plane. An alternative optical modulator such as liquid crystal-onsilicon can help avoid this undesired effect. Second, due to the limited sub-aperture size, required exposure times were quite long. Moving to an alternative aperture coding strategy to increase light throughput, while maintaining sufficient overlap, can immediately address this shortcoming. Adopting a phase-only modulation strategy is another possible direction to decrease total capture time. Third, this work does not consider the effects of a pixel-limited optical system. A fully optimized OFC setup will match each sub-aperture image PSF to the detector pixel size, leading to a 3X resolution improvement above the sensor pixel count with a very simple modification to our current hardware. Furthermore, two primary shortcomings currently limit SA-OFC algorithm performance. First, as already noted above, computational scaling issues require the annealing search to consider only a small number of aberrations (i.e., Zernike coefficients) and misalignment parameters. Second, instabilities are introduced when solving for multiple modes that may not be orthogonal. Although each Zernike mode is mutually orthogonal and thus does not confuse the annealing’s movement towards minimal error, the geometric alignments and aberrations are not necessarily orthogonal, which might lead to algorithm stagnation. Future work will examine more robust methods to search within a high-dimensional unknown parameter space. Possibilities include gradient descent or a maximum likelihood-based solver, where the log of the aberration function presents itself as a simple linear sum. These alternative strategies may additionally benefit from a modified aperture coding strategy. 115 Bibliography [1] R. Horstmeyer, X. Ou, J. Chung, G. Zheng and C. Yang, “Overlapped Fourier coding for optical aberration removal,” Opt. Express 22(20), 24062–24080 (2014). [2] A. Lohmann, “Scaling laws for lens systems,” Appl. Optics 28, 4996–4998 (1989). [3] G. Zheng, R. Horstmeyer and C. Yang, “Wide-field, high-resolution Fourier ptychographic microscopy,” Nature Photon. 7, 739–745 (2013). [4] J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Optics 21, 2758–2769 (1982). [5] V. Elser, “Phase retrieval by iterated projections,” J. Opt. Soc. Am. A 20, 40-56 (2003). [6] H. M. L. Faulkner and J. M. Rodenburg, “Movable aperture lensless transmission microscopy: A novel phase retrieval algorithm,” Phys. Rev. Lett. 93, 023903 (2004). [7] A. M. Maiden and J. M. Rodenburg, “An improved ptychographical phase retrieval algorithm for diffractive imaging,” Ultramicroscopy 109, 1256–1562 (2009). [8] C. J. Schwarz, Y. Kuznetsova and S. R. J. Brueck, “Imaging interferometric microscopy,” Opt. Lett. 28, 1424–1426 (2003). [9] T. R. Hillman, T. Gutzler, S. A. Alexandrov and D. D. Sampson, “High-resolution, wide-field object reconstruction with synthetic aperture Fourier holographic optical microscopy,” Opt. Express 17, 7873–7892 (2009). [10] S. Chowdhury and J. Izatt, “Structured illumination diffraction phase microscopy for broadband, subdiffraction resolution, quantitative phase imaging,” Opt. Lett. 39, 1015–1018 (2014). [11] R. Gao, G. Pedrini and W. Osten, “Phase retrieval with resolution enhancement by using structured illumination,” Opt. Lett. 38, 5204–5207 (2013). [12] D. J. Lee and A. M. Weiner, “Optical phase imaging using a synthetic aperture phase retrieval technique,” Opt. Express 22, 9380–9394 (2014). [13] L. Camacho, V. Mico, Z. Zalevsky and J. Garcia, “Quantitative phase microscopy using defocusing by means of a spatial light modulator,” Opt. Express 18, 6755–6766 (2010). 116 [14] M. P. Lee, G. M. Gibson, R. Bowman, S. Bernet, M. Ritsch-Marte, D. B. Phillips and M. J. Padgett, “A multi-modal stereo microscope based on a spatial light modulator,” Opt. Express 21, 16541–16551 (2013). [15] B. Bhaduri, H. Pham, M. Mir and G. Popescu, “Diffraction phase microscopy with white light,” Opt. Lett. 37, 1094–1096 (2012). [16] S. Dong, R. Horstmeyer, R. Shiradkar, K. Guo, X. Ou, Z. Bian, H. Xin and G. Zheng, “Aperturescanning Fourier ptychography for 3D refocusing and super-resolution macroscopic imaging,” Opt. Express 22, 13586–13599 (2014). [17] C. Liang, T. Lin, B. Wong, C. Liu and H. Chen, “Programmable aperture photography: multiplexed light field acquisition,” ACM Trans. Graph. 27(3), 55 (2008). [18] K. Marwah, G. Wetzstein, Y. Bando and R. Raskar, “Compressive light field photography using overcomplete dictionaries and optimized projections,” ACM Trans. Graph. 32(4), 46 (2013). [19] S. Sarder and A. Nehorai, “Deconvolution methods for 3-D fluorescence microscopy images,” IEEE Sig. Proc. Mag. 23(3), 32–45 (2006). [20] M. Broxton, L. Grosenick, S. Yang, N. Cohen, A. Andalman, K. Deisseroth and M. Levoy, “Wave optics theory and 3-D deconvolution for the light field microscope,” Opt. Express 21(21), 25418–25439 (2013). [21] S. Kirkpatrick, C. D. Gelatt and M. P. Vecchi, “Optimization by simulated annealing,” Science 220, 671–680 (1983). [22] A. M. Maiden, M. J. Humphry, M. C. Sarahan, B. Kraus and J. M. Rodenburg, “An annealing algorithm to correct positioning errors in ptychography,” Ultramicroscopy 120, 64–72 (2012). [23] A. Tripathi, I. McNulty and O. G. Shpyrko, “Ptychographic overlap constraint errors and the limits of their numerical recovery using conjugate gradient descent methods,” Opt. Express 22, 1452–1466 (2013). [24] R. G. Paxman, T. J. Schulz, and J. R. Fienup, “Joint estimation of object and aberrations by using phase diversity,” JOSA A 9(7), 10721085 (1992). [25] B. M. Hanser, M. G. Gustafsson, D. A. Agard, and J. W. Sedat, “Phase-retrieved pupil functions in wide-field fluorescence microscopy,” J. Microsc. 216(1), 3248 (2004). [26] G. Zheng, X. Ou, R. Horstmeyer and C. Yang, “Characterization of spatially varying aberrations for wide field-of-view microscopy,” Opt. Express 21, 15131–15143 (2013). 117 [27] O. Bunk, M. Dierolf, S. Kynde, I. Johnson, O. Marti and F. Pfeiffer, “Influence of the overlap parameter on the convergence of the ptychographical iterative engine,” Ultramicroscopy 108, 481–487 (2008). [28] R. Horstmeyer and C. Yang, “A phase space model of Fourier ptychographic microscopy,” Opt. Express 22, 338–358 (2014). [29] P. Thibault and A. Menzel, “Reconstructing state mixtures from diffraction measurements,” Nature 494, 68–71 (2013). [30] X. Ou, G. Zheng and C. Yang, “Embedded pupil function recovery for Fourier ptychographic microscopy,” Opt. Express 22, 4960–4972 (2014). [31] R. Kingslake and R. B. Johnson, Lens Design Fundamentals: Second Edition (Elsevier and SPIE Press, 2010). [32] P. E. Debevec and J. Malik, “Recovering high dynamic range radiance maps from photographs,” In Proc. SIGGRAPH 97, ACM SIGGRAPH/Addison Wesley Computer Graphics Proceedings, Annual Conference Series, 369–378 (1997). 118 Chapter 8 Diffraction tomography with Fourier ptychography In this chapter, we modify the Fourier ptychography reconstruction algorithm to determine the complex index of refraction, across a three-dimensional volume, of a thick sample. We term this procedure Fourier ptychographic tomography (FPT), which conceptually shares the sample goal as diffraction tomography (DT). In DT, one typically illuminates a sample from a variety of directions, captures each resulting complex diffracted field with holographic techniques, and applies a computational procedure to reconstruct the three-dimensional sample structure, assuming minimal inter-sample scattering. With FPT, we no longer need to measure both the amplitude and phase of the diffracted fields, but instead only their intensity. Thus, holography is not needed, and we can now obtain limited-angle complex tomographic reconstructions with micrometer-scale resolution, using a standard microscope. 8.1 Introduction It is challenging to image thick samples with a standard microscope. High-resolution objective lenses offer a shallow depth-of-field, which require one to axially scan through the sample to visualize three-dimensional shape. Unfortunately, refocusing does not remove light from areas above and below the plane of interest. This longstanding problem has inspired a number of solutions, the most widespread being confocal designs, two-photon excitation methods, light sheet microscopy, and optical coherence tomography. These above methods “gate out” light from sample areas away from the point of interest. They offer excellent signal enhancement, especially for thick, fluorescent samples [1]. Such gating techniques also encounter several problems. First, they typically must scan out each image, which might require physical movement, and can be time consuming. Second, the available signal (i.e., the number of ballistic photons) decreases exponentially with depth. To overcome this 119 limit, one must use a high NA lens, which provides a proportionally smaller image field-of-view (FOV). Finally, little light is backscattered when imaging non-fluorescent samples that are primarily transparent, such as commonly seen in embryology, in model organisms such as zebrafish, and after the application of recent tissue-clearing [2] and expansion [3] techniques. Instead of capturing just the ballistic photons emerging from the sample, one might instead image the entire optical field, including the scattered components. This avoids point scanning, and allows one to record a very wide image FOV in a single snapshot. Several techniques have been proposed to enable depth selectivity after full-field capture. First, one might perform optical sectioning by capturing a focal stack, and then attempting digital deconvolution [4]. A second related example is light-field imaging [5,6]. Point-spread function engineering is a third possibility [7], but this typically requires a sparse sample. All three of these methods primarily operate with incoherent light, e.g. from fluorescent samples. They are thus not ideal tools for obtaining the complex refractive index distribution of a primarily transparent medium. To do so, it is useful to use coherent illumination. For example, the amplitude and phase of a digital hologram may be computationally propagated to different depths within a thick sample, much like refocusing a microscope. However, it is still influenced by the field at out-of-focus planes. Several techniques have improved upon depth selectivity with quasi-coherent illumination, based upon the acquisition of multiple images [8–12]. A very useful framework to summarize how coherent light scatters through thick samples is diffraction tomography (DT) [13]. This framework connects the optical fields that diffract from a sample, under arbitrary illumination, to its 3D complex refractive index. In a typical DT experiment, one illuminates a sample of interest with a series of tilted plane waves and measures the resulting complex diffraction patterns in the far field. These measurements may then be combined with a suitable algorithm into a tomographic reconstruction. As a synthetic aperture technique, DT comes with the additional benefit of improving the limited resolution of an imaging element beyond its traditional diffraction cutoff [14]. Thus, it appears a well-suited method for the study of thick, transparent samples at high resolution. However, as a technique that models both amplitude and phase of a coherent field, nearly all prior implementations of DT required a reference beam and holographic measurement, or some sort of phase-stable interference (including SLM coding strategies, e.g. as in [22]). Reference fields require sub-micrometer stability in terms of both motion and phase drift, which has thus far limited DT to well-controlled, customized setups. While several prior works have considered solving the DT problem from intensity-only measurements from a theoretical perspective [23–27], none have implemented a DT system within a standard microscope, or connected their reconstruction attempts to ptychography. Here, we perform DT based upon intensity images from variable LED illumination with an 120 (a) Setup! (b) Raw image set! gn! g2! g1! Detector! Intensities g(x,j)! …   j = 2! j = 1! (c) 3D k-space scattering potential ! Microscope objective! Spectrum Uj (k)! input! kzmax! kz! FPT algorithm! -kzmax! kx! kmax! kmax! Thick Sample! j = q! ky! output! Scattered field Uj (r)! (d) 3D reconstruction of complex refractive index! 30! 15! j th LED illumination ! LED Matrix! Figure 8.1: Fourier ptychographic tomography (FPT) setup. (a) A labeled diagram of the FPT microscope, with optical fields of interest labeled. (b) Multiple images are acquired under varied LED illumination. (c) A ptychography-inspired algorithm combines these images within a 3D kspace representation of the complex sample of interest. (d) FPT outputs a 3D tomographic map of the complex index of refraction of the sample. Included images are experimental measurements (starfish embryo). array source. Our technique, termed Fourier ptychographic tomography (FPT), captures a sequence of images while changing the light pattern displayed on the LED array. Then, it combines these images using a phase retrieval-based ptychographic reconstruction algorithm, which computationally (as opposed to physically) rejects light from all areas above and below each plane of interest. FPT also improves the lateral image resolution beyond the standard cutoff of the objective lens used for imaging. The end result is a quantitatively accurate three-dimensional map of the complex index of refraction of a volumetric sample, obtained directly from a sequence of standard microscope images. 8.2 Related Work The theoretical foundations of DT were first developed by Wolf [13]. A number of implementations based upon holography have followed. An early demonstration by Lauer is a good example [14]. Prior methods have also implemented tomography within a microscope-like setup, but required the addition of a phase-stable reference beam. The first results operated under the projection approximation, which models light as a ray [15]. Subsequent work has taken the effects of diffraction into account [16–18]. 121 Instead of relying upon holography, this work measures intensity images and computationally recovers the missing phase. As mentioned above, a few prior works consider the reconstruction problem from detected optical intensities, but must either move the focal plane between measurements, or assume a sample support constraint. They do not attempt ptychographic phase retrieval. Tomographic phase retrieval based on lifting [19] is another related strategy. The connection between the first Born approximation and phase retrieval has also been explored within the context of volume hologram design [20]. Related efforts to reconstruct volumetric samples from wide-field intensity-only measurements outside of the realm of DT include lensless on-chip setups [28, 29], lensless techniques that assume an appropriate linearization [30], and methods relying upon effects like defocus (e.g., the transport of intensity equation [31]) or spectral variations [21]. None of these techniques fit within a standard microscope setup, nor offer the ability to simultaneously improve spatial resolution. Using very similar hardware, Fourier ptychography (FP) [32] can simultaneously improve image resolution and measure quantitative phase [33]. However, it is restricted to thin samples. FPT effectively extends prior developments of FP into the third dimension. One recent work also examined the problem of 3D imaging from intensities in an LED microscope [34]. This recent example borrows its 3D reconstruction technique from the related field of 3D ptychography [35, 36], where the sample under examination is split up into a specified number of infinitesimally thin slices, and the beam propagation method (i.e., multi-slice approximation) is used [37]. Unlike the multi-slice approach, which works well with distinctly separated absorbing layers, FPT is best suited for continuous, primarily transparent samples. A number of related methods to perform 3D X-ray ptychography have also been proposed [38–40]. However, none seem to apply DT under the first Born or Rytov approximation, to the best of our knowledge. A popular technique appears to use standard 2D ptychographic solvers to determine the complex field for individual projections of a slowly rotated sample, which are subsequently combined using standard DT techniques, as shown with both crystallographic [41] and unordered specimens [42]. Here, we first outline a solid foundation for the application of ptychographic phase retrieval to DT. Unlike approaching the problem from a projection-based or multi-slice perspective, the framework of DT (under the first Born approximation) follows directly from the scalar wave equation. It offers a clear picture of achievable resolution in 3D, spells out sampling and data redundancy requirements for an accurate reconstruction, and presents a clear path forward for future extensions to account for multiple scattering [43]. Furthermore, our method does not require the arbitrary assignment of the number slices in the 3D volume, or their location, or for us to select a particular order in which to address each slice as iterations proceed. Instead, it simply inserts the measured 3D data into its appropriate location in Fourier space and ensures phase consistency between each measured image, given a sufficient amount of data redundancy (just like ptychography). From the initial 122 (b) k-space representation! (a) Image formation! Detector! Aperture! g(x,j)   -kz! Uj(kx)Ÿa(kx)   |k|   Scattering potential! V(kx, kz)! kx! kj! Uj(x')   Thick sample! θj   Apertures shift along illumin.! sphere ! V(x,z)   (c) Reconstructed bandpass volume! LEDs! θj   -kz! z! 2π/Δz! Ve(kx, kz)! x! kx! 2π/Δx! Figure 8.2: Mathematical summary of FPT. (a) Diagram of the FPT setup in 2D. The field from the jth LED scatters through the sample and exits its top surface as Uj (x0 ). This field then propagates to form Ûj (kx ) at the microscope back focal plane, where it is band-limited by the finite microscope k(z)   then propagates to the aperture, a(kk(z)   k(z)  image plane, where its intensity is x ). This band-limited field sampled to form the jth image. (b) Under the first Born approximation, each detected image is the squared magnitude of the Fourier transform of one ”shell” in (kx , kz ) space. (c) By filling in this space with a ptychographic phase retrieval algorithm, FPT reconstructs the complex values within k(y)   k(x)   the finite bandpass Fourier transform of k(y)   this reconstruction yields k(y)   our complex sample k(x)   volume. The k(x)   index of refraction map. starting point of solving for the first term in the Born expansion, we aim this approach as a general framework to eventually solve the challenging problem of forming tomographic maps of volumetric samples, at sub-micrometer resolution, in the presence of significant scattering. 8.3 Methods In this section, we develop a mathematical expression for our image measurements using the FPT framework, and then summarize our reconstruction algorithm. We describe our setup and reconstruction in 3D with the vector r = (rx , ry , rz ) defining the sample coordinates and the vector k = (kx , ky , kz ) defining the k-space (wavevector) coordinates (see Fig. 8.1). 123 8.3.1 Image formation in FPT It is helpful to begin our discussion by introducing a quantity termed the scattering potential, which contains the complex index of refraction of an arbitrarily thick volumetric sample: V (r) =
k n2 (r) − n2b . 4π (8.1) Here, n(r) is the spatially varying and complex refractive index profile of the sample, nb is the index of refraction of the background (which we assume is constant), and k = 2π/λ is the wavenumber in vaccuum. It is informative to point out that n(r) = nr (r) + ini (r), where nr is associated with the sample’s refractive index and ni is associated with its absorptivity. We typically neglect the dependence of n on λ since we illuminate with quasi-monochromatic light. This dependence cannot be neglected when imaging with polychromatic light. Next, to understand what happens to light when it passes through this volumetric sample, we define the complex field that results from illuminating the thick sample, U (r), as a sum of two fields: U (r) = Ui (r) + Us (r). Here, Ui (r) is the field incident upon the sample (i.e., from one LED) and Us (r) is the resulting field that scatters off of the sample. As detailed in [13], we may insert this decomposition into the scalar wave equation for light propagating through an inhomogeneous medium. Using Green’s theorem, we may determine the total field scattered by the medium as, 0 0 U (r ) = Ui (r ) + Z G(|r0 − r|)V (r)U (r)dr. (8.2) Here, G(|r0 − r|) is the Green’s function connecting light scattered from various sample locations, denoted by r, to an arbitrary location r0 . V (r) is the scattering potential from Eq. 8.1. Since U (r) is unknown at all sample locations, it is challenging to solve Eq. 8.2. Instead, it is helpful to apply the first Born approximation, which replaces U (r) in the integrand with Ui (r). This approximation assumes that Ui (r)  Us (r). It is the first term in the Born expansion that describes the scattering response of an arbitrary sample [13]. Our system sequentially illuminates the sample with an LED array, which contains q = qx × qy sources positioned a large distance l from the sample (in a uniform grid, with inter-LED spacing c, see Fig. 8.1). It is helpful to label each LED with a 2D counter variable (jx , jy ), where −qx /2 ≤ jx ≤ qx /2 and −qy /2 ≤ jy ≤ qy /2, as well as a single counter variable 1 ≤ j ≤ q. Assuming each LED acts as a spatially coherent and quasi-monochromatic source (central wavelength λ), the incident field takes the form of a plane wave traveling at a variable angle such that θjx = tan−1 (jx · c/l) and θjy = tan−1 (jy · c/l) with respect to the x and y axes, respectively. We may express the jth field incident upon the sample as, (j) Ui (r) = exp(ikj · r), (8.3) 124  where kj = k sin θjx , k sin θjy , k q 1 − sin2 θjx − sin2 θjy  is the wavevector of the jth LED plane wave. As θjx and θjy vary, kj will always assume values along a spherical shell in 3D (kx , ky , kz ) space (i.e., the Ewald sphere), since the value of kz is a deterministic function of kx and ky . (j) After replacing U (r) in Eq. 8.2 with Ui (r) in Eq. 8.3, and additionally approximating the Green’s function G as a far field response, the following relationship emerges between the scattering potential V and the jth scattered field, Û (j) (r0 ), in the far field [13]: Û (j) (r0 ) = exp(ikj · r0 ) + V̂ (k − kj ) (8.4) Here, V̂ (k) is the three-dimensional Fourier transform of V (r), which we refer to as the k-space scattering potential, and k denotes the scattered wavevector. The total field scattered through the sample and viewed at a distant plane, Û (j) (r0 ), is given as the original incident LED plane wave plus the values along a specific manifold, or spherical “shell”, of the k-space scattering potential, here written as V̂ (k − kj ). We illustrate the geometric connection between V̂ (k) and Û (j) (r0 ) for a 2D optical geometry in Fig. 8.2(b). The center of the jth shell is defined by the incident wavevector, kj . For a given shell center, the radial distance to each value of interest is given by |k| = k (see multi-colored arcs in Fig. 8.2(b)). As kj varies with the changing LED illumination, the shell center shifts along a second shell with similar radius (since kj is itself constrained to lie on an Ewald sphere, see gray circle in Fig. 8.2(b)). The goal of DT is to determine all complex values within the volume V̂ , from a set of q scattered fields, {Û }qj=1 . It is common to measure these scattered fields holographically [14, 18]. Each 2D holographic measurement maps to the complex values of V̂ that lie along one 2D shell. Values from multiple measurements (i.e., the multiple shells in Fig. 8.2(b)) can be combined to form a k-space scattering potential estimate, V̂e [44]. Nearly all stationary optical setups will yield only an estimate, since it is challenging to measure data from the entire k-space scattering potential without rotating the sample. Fig. 8.1(c) and Fig. 8.2(c) each display a typical measurable volume, also termed a bandpass, from a limited-angle illumination and detection setup. Once sampled, an inverse 3D Fourier transform of the band-limited V̂e (k) yields the desired complex scattering potential estimate, Ve (r), from which the quantitative index of refraction is directly obtained. In FPT, we do not measure the scattered fields holographically. Instead, we use a standard microscope to detect image intensities and apply a ptychographic phase retrieval algorithm to solve for the unknown complex potential. The scattered fields in Eq. 8.4 are defined at the microscope objective back focal plane (i.e., its Fourier plane), where the coordinate r0 is conjugate to the sample plane coordinate r. We may thus replace r0 with the 2D k-space coordinate system k2D = (kx , ky ) of the back focal plane. After neglecting the effect of the background plane wave term, jth shifted field at our microscope back focal plane is simply Û (j) (k2D ) = V̂ (k − kj ). 125 Each shifted scattered field is then bandlimited by the microscope aperture function, a(k2D ), before propagating to the image plane. The limited extent of a(k2D ) (defined by the imaging system NA) sets the maximum extent of each shell along kx and ky . The jth intensity image acquired by the detector is given by the squared Fourier transform of the bandlimited field at the microscope back focal plane: h i2 g(x, y, j) = F V̂ (k − kj ) · a(k2D ) . (8.5) Here, F denotes a 2D Fourier transform and we neglect the effects of magnification, for simplicity, by assuming the image plane coordinates match the (x, y) coordinates at the sample plane. The goal of FPT is to determine the complex 3D function V̂ (k) from the real, positive data matrix, g(x, y, j). A final 3D Fourier transform of V̂ (k) yields the desired scattering potential, and subsequently the refractive index distribution, of the thick sample. 8.3.2 FPT reconstruction algorithm Eq. 8.5 closely resembles the data matrix measured by Fourier ptychography (FP). Now, however, intensities are sampled from a volumetric function along shells in a 3D space (i.e., the curves in Fig. 8.2). We use an iterative reconstruction procedure, mirroring that from FP [32], to “fill in” the k-space scattering potential with data from each recorded intensity image. Just like FP needs a certain amount of data redundancy (i.e., overlap in k-space) to recover the unknown optical phase, FPT also requires overlap between shell regions in 3D k-space. Since our discretized k-space now has an extra dimension, overlap is less frequent and more images are required for successful algorithm convergence. We may encourage overlap with increased discretization, a smaller LED array pitch and/or a larger array-sample distance along z. As we demonstrate experimentally, several hundred images are sufficient for a complex reconstruction that offers approximately 35 unique slices along the axial dimension. It is important to select the correct limits and discretization of 3D k-space (i.e., the FOV and resolution of the complex sample reconstruction). The maximum resolvable wavevector along kx and ky is proportional to k(NAo + NAi ), where NAo is the objective NA and NAi is maximum NA of LED illumination. This lateral spatial resolution limit matches FP [45]. The maximum resolvable wavevector range along kz is also determined   as a function of the objective and illumination NA q q as, kzmax = k 2 − 1 − NA2o − 1 − NA2i . This relationship is easily derived from the geometry of the k-space bandpass volume in Fig. 8.2, as shown in [14]. We typically specify the maximum imaging range along the axial dimension, zmax , to equal approximately twice the expected sample thickness. This then sets the discretization level along kz : ∆kz = 2π/zmax . The total number of resolved slices along z is set by the ratio kzmax /∆kz . We now summarize the FPT reconstruction algorithm in the following 5 steps: 126 1. Initialize a discrete estimate of the unknown k-space scattering potential, V̂e (r0 ), selecting the appropriate 3D array size following the discussion above. Either a single raw image may be padded along all three dimensions and then Fourier transformed for this initialization, or the raw intensity imagery may be used to form a refocused light field [34]. A constant matrix is also often an adequate initialization. 2. For j = 1 to q images, compute the center coordinate, kj , and select its associated shell (radius k, maximum width 2k · NAo ). This selection process samples a discrete 2D function, dˆj (kx , ky ) from the 3D k-space volume. The selected voxels must partially overlap with voxels from adjacent shells. Currently, no interpolation is used to map voxels from the discrete shell to pixels within dˆj (kx , ky ). 3. Fourier transform dˆj (kx , ky ) to the image plane, and constraint its amplitudes to match the measured amplitudes from the jth image. For example, the update may take the simple form, p d0j (x, y) = g(x, y, j)·dj (x, y)/|dj (x, y)|. More advanced alternating projection-based updates are also available [46]. 4. Inverse 2D Fourier transform the image plane update, d0j (x, y), back to 2D k-space to form dˆ0j (kx , ky ). Replace the voxel values of V̂e (r0 ) at locations where voxel values were extracted in step 2. Use the values of dˆ0j (kx , ky ) for replacement. 5. Repeat this select, update and replace process for all j = 1 to q images. This completes one iteration of the FPT algorithm. Continue for a fixed number of iterations, or until satisfying some error metric. At the end, 3D inverse Fourier transform V̂e (k) to recover the complex scattering potential, Ve (r). In practice, we also implement a pupil function recovery procedure [47] as we update each extracted shell from k-space. This allows us to simultaneously estimate and remove possible aberrations present in the microscope back focal plane. 8.4 Results We experimentally verify our reconstruction technique using a standard microscope outfitted with an LED array. The microscope uses an infinity corrected objective lens (NAo = 0.4, Olympus MPLN, 20X), to image onto a digital detector containing 4.54 µm pixels (Prosilica GX 1920, 1936×1456 pixel count). The LED array contains 31 × 31 surface-mounted elements (model SMD3528, center wavelength λ =632 nm, 4 mm LED pitch, 150 µm active area diameter). For this first demonstration, we position the LED array approximately 135 mm beneath the sample to create a maximum illumination NA of NAi =0.41. This leads to an effective lateral NA of NAo + NAi = 0.81, and a 127 (a) Raw image! 3 µm! (b) FPT reconstruction, Δz=0! 1 µm! 0! 3 µm! 1 µm! 0.8 µm! +.06! Re[Δn]! Intensity (AU)! 0.8 µm! y! 1! -.04! x! Figure 8.3: FPT improves the lateral resolution of a standard microscope. (a) A single raw image of a layer of 0.8 µm microspheres immersed in oil, where beads within each cluster are not resolved. (b) The real component of the index of refraction from one slice out of thirty for our FPT reconstruction Other   ok  color   aps:   (∆z = 0 slice), which clearly resolves each m microsphere. lateral resolution gain along (x, y) of slightly over a factor of 2 (from a 1.6 µm minimum resolved spatial period in the raw images to a 0.78 µm minimum resolved spatial period in the reconstruction). The associated axial resolution is computed at 3.7 µm, and we reconstruct quantitative sample information across a depth range of approximately zmax =110 µm (approximately 20 times larger than the stated objective lens DOF of 5.8 µm). For most of the reconstructions presented below, we capture n = 675 images from the same fixed pattern of LEDs, which are then input into the FPT algorithm. We typically use the following parameters for reconstruction: each raw image is cropped to 1000 × 1000 pixels, the reconstruction voxel size is set at 0.39 × 0.39 × 3.7 µm3 for sampling at the Nyquist-Shannon rate, the reconstruction array contains approximately 2100 × 2100 × 30 voxels, and the algorithm runs for 5 iterations. 8.4.1 Quantitative verification We include three different quantitative verifications of FPT performance using polystyrene microspheres as reference targets. First, we verify the ability of FPT to improve lateral image resolution. This matches the goal of FP for thin 2D samples. Here, our sample consists of 800 nm diameter microspheres (index of refraction ns = 1.59) immersed in oil (index of refraction no = 1.515). We highlight a small group of these microspheres in Fig. 8.3. First, we show a single raw image in (a) (generated from the center LED), where the individual spheres, gathered in small clusters, are not resolved at all. Based upon the coherent Sparrow limit for resolving two points (0.68λ/NAo ) this raw image cannot resolve points that are closer than 1.1 µm. After FPT reconstruction, we obtain the complex index of refraction in Fig. 8.3(b), where here we show the real component of the recovered index. The FPT reconstruction along the z = 0 plane clearly resolves the spheres within each cluster. This distance is close to the new expected Sparrow limit of 0.68λ/ (NAo + NAi ) = 0.54µm, thereby verifying lateral resolution performance close to, but not exactly at, the theoretical limit. 128 (a) 3D quantitative refractive index! Refractive index Δnr (x, y, z=0)! +.01! Refractive index Δnr (x, 0, z)! z! z! y! x! 5 µm! Index along x! 0! Index along z! .012! Δnr! -.005! x! 5 µm! -.005! 0! x (µm)! 15! 0! z (µm)! 20! (b) Comparing reconstructions along z! z! z! 5 µm! Light field (intensity)! 2D FP with refocus! 3D FP! z! x! x! z! x! x! Figure 8.4: FPT quantitatively measures the complex index of refraction of samples in 3D. (a) Tomographic reconstruction of 12 µm microspheres immersed in oil, where we show a lateral (∆z = 0) slice on the left, an axial (∆y = 25 µm) slice on the right, and 1D plots of the index shift along both x and z, demonstrating quantitative performance. (b) We use the same dataset to obtain an FP reconstruction and propagate the result along z (middle), and also perform light field refocusing (right). Our FPT reconstruction (left) offers the closest match to the expected axial profile of a spherical bead. 129 Second, we check the quantitative accuracy of FPT by imaging microspheres that extend across more than just a few reconstruction voxels. Fig. 8.4 displays a reconstruction of 12 µm diameter microspheres (index of refraction ns = 1.59) immersed in oil (index of refraction no = 1.58). We use the same data capture and post-processing steps as in Fig. 8.3. The full reconstructed scattering potential, cropped to 200×200×15 voxels, is shown in Fig. 8.4. We again display the real (nonabsorptive) component of the recovered index across both a lateral slice (along the z = 0 plane), a vertical slice (along the y = 25 µm plane). Detailed 1D traces along the center of both of these slices are also included. Three observations are noteworthy regarding this experiment. First, the measured index shift approximately matches the expected shift of ∆n = ns − no = 0.01 across the entire bead, thus demonstrating quantitatively accurate performance. Second, for each 1D trace through the center of each microsphere, we would ideally expect a perfect rect function (from ∆n = 0 to ∆n = 0.01 and then back down). This is unlike 2D FP, which reconstructs the phase delay though each sphere, leading to a parabolic function (due their varying thickness along the optical axis). While the system can resolve an approximate step function through the center of the sphere along the lateral (x) dimension, it is not a step function function along the axial (z) direction. This is caused by the limited extent of the measurable volume of 3D k-space (i.e., the limited bandpass). The “missing cone” of information, primarily surrounding the kz axis, creates a noticeably wide pointspread function along z, which leads to its distinct sinc shape. Various methods are available to computationally fill in the missing cone using prior sample information [48, 49]. For our third observation, we compare FPT with two alternative techniques for 3D imaging in Fig. 8.4(c). First, we use the same dataset to perform 2D FP, and then attempt to holographically refocus its complex 2D solution. We obtain this solution using the same number of images (q = 675), with the procedure in [32], after focusing the objective lens at the axial center of the 12 µm microspheres. The “out of focus noise” above and below the plane of the microsphere, created by digital propagation of the complex field via the angular spectrum method, quite noticeably hides its spherical shape. Second, we interpret the same raw image set as a light field and perform light field refocusing [5]. While the refocused light field clearly resolves the outline of microsphere along the z-dimension, it does not offer a quantitative picture of the sample interior, nor a measure of its complex index of refraction. The areas above the microsphere are very bright due to its lensing effect (i.e., the light field displays the optical intensity at each plane, and thus displays high energy where the microsphere focuses light). Compared to these two alternatives, FPT more accurately measures the microsphere’s ground-truth 3D shape. For our third and final quantitative test, we verify the axial resolution of FPT along z. Here, we prepare a sample containing two closely separated layers. Each layer contains 2 µm microspheres (ns = 1.59) distributed across the surface of a glass slide, which we sandwich together, with oil 130 (a) Sample geometry! 2 µm! Standard microscope! (b) Raw image, Δz = 0! (c) Re-focused +1.9 µm! ! 2 µm! 1! Inten. (AU)! 0! x! z! 2  µm   3.9 µm! 2  µm   3D FPT! (d) Δz = -1.9 µm! (e) Δz = 0 µm! (f) Δz = +1.9 µm! .09! Abs[Δn]! 0! 2  µm   2  µm   2  µm   Figure 8.5: Experimentally measuring the axial resolution of FPT. (a) The reconstructed sample contains two layers of microspheres separated by a thin layer of oil. Raw images (b) focused at the center of the two layers and (c) on the top layer do not clearly resolve overlapping microspheres (e.g., in red box). (d)-(e) Slices of the FPT tomographic reconstruction, showing |∆n|, clearly resolve the two individual sphere layers. 131 between (no = 1.515). The separation between the two microsphere layers, measured from the center of each sphere along z, is 3.9 µm (i.e., the separation between the microscope slide surfaces is 5.9 µm, as diagrammed in Fig. 8.5(a)). The 3.9 µm center-to-center distance is close to the expected axial resolution limit of 3.7 µm for the FPT microscope. Conventional microscope images of the sample, using the center LED for illumination, are in Fig. 8.5(b)-(c), where we focus to the center of the two layers (∆z = 0) and the top layer (∆z = 1.9 µm), in an attempt to distinguish each microsphere layer. At the top of each image, where two layers of microspheres overlap, it is especially hard to resolve each sphere, or determine which sphere is in a particular layer. Next, we return the focus to the ∆z = 0 plane and implement FPT. We display three slices of our reconstructed scattering potential in Fig. 8.5(d)-(f). Here, we show the absolute value of the potential near the plane of the top layer, at the center, and near the plane of the bottom layer. The originally indistinguishable spheres within the top and bottom layers are now clearly resolved in each z-plane. Due to the system’s limited axial resolution, the reconstruction at the middle plane (∆z = 0) still shows the presence of spheres from both the top and bottom layers. Comparing Fig. 8.5(b)-(c) with Fig. 8.5(e)-(f) clearly maintains that the axial resolution of FPT is sharper than manual refocusing. Not only is each layer clearly distinguishable (as predicted theoretically), but we now also have quantitative information about the complex refractive index. 8.4.2 Biological experiments Next, we use FPT to reconstruct the 3D complex refractive index distribution of two different thick, biological specimens. First, we tomographically reconstruct a 3D volume containing a trichinella spiralis parasite (Fig. 8.6). Here, since the worm extends along a larger distance than the width of our detector, we performed FPT twice, shifting the FOV between to capture the left and right side of the worm with 10% overlap between. We then merged each tomographic reconstruction together with a simple averaging operation (matching that from FP [32]). The total imaging volume here is 0.8 mm×0.4 mm×110 µm. Note that if we were to replace our current digital detector with one that occupied the entire microscope FOV, we would increase our fixed imaging volume to 2.2 mm×2.2 mm×110 µm, and form tomograms that each contain approximately 109 voxels. We display a thresholded 3D scattering potential reconstruction of the parasite at the top of Fig. 8.6 (real component, threshold applied at Re[∆n] > 0.7 after |∆n| normalized to 1, undersampled for clarity). Its 3D curved trajectory is especially clear in the 3 separate z-slices of the reconstructed tomogram in Fig. 8.6(a). The two downward bends in the parasite body are lower than the upward bend in the middle, as well as at its front and back ends. It is very challenging to resolve these depth-dependent sample features by simply refocusing a standard microscope. Fig. 8.6(b) displays such an attempt, where the same three z planes are brought into focus manually. Since 132 z (mm)! 0.1! 1! Re(Δn), norm.! 0! 0.8! 0! 0.7! 0! Δz = -7.4 µm   x (mm)! Δz = 0 µm   0.4! Δz = +7.4 µm   1! (a) 3D FPT! Re(Δn), norm.! 50 µm! Δz = 0 µm   Δz = +7.4 µm   1! Intensity (AU)! (b) Manual refocus! Δz = -7.4 µm   0.7! 50 µm! 0! Figure 8.6: Tomographic reconstruction of a trichinella spiralis parasite. (a) The worm’s curved trajectory is clearly resolved within the various z-planes. (b) Refocusing the same distance to each respective plane does not clearly distinguish each in-focus worm segment (marked by white arrows). Since the worm is primarily transparent, in-focus worm sections exhibit minimal intensity contrast, presenting significant challenges for segmentation. FPT, on the other hand, exhibits maximal contrast at each voxel containing the worm. 133 the sample is primarily transparent, in-focus areas in each standard image actually exhibit minimal contrast (see marked arrows in Fig. 8.6(b)), making any attempt at depth segmentation (e.g., deconvolution of a focal stack [4]) nearly impossible. Since FPT is a phase contrast technique, points along the parasite offer maximal contrast in each respective voxel, enabling direct segmentation via thresholding, as achieved in the top plot. For our second 3D biological example, we tomographically reconstruct a starfish embryo at its larvae stage (Fig. 8.7). Here, we again show three different closely spaced z-slices of the reconstructed scattering potential (Re(∆n), no thresholding applied). Each z-slice contains sample features that are not present in the adjacent z-slices. For example, the large oval structure in the upper left of the ∆z = 0 plane, which is a developing stomach, nearly completely disappears in the ∆z = −3.3 µm plane. Now at this z-slice, however, small developmental structures appear in the lower right (marked as feature 1). Both the particular plane of the developing stomach and even the presence of feature 1 are completely missing from the refocused images. This is due to the inability of the standard microscope to segment each particular plane of interest, the inability to accurately reconstruct transparent structures without a phase contrast mechanism, and an inferior lateral (x, y) resolution with respect to FPT. 8.5 Conclusion The FPT method performs diffraction tomography using intensity measurements, captured with a standard microscope and an LED illuminator. Its reconstruction algorithm extends previous work with FP to now operate in 3D. The current system offers a lateral resolution of approximately 400 nm (550 nm at the Sparrow limit, 800 nm full period) and an axial resolution of 3.7 µm at the Nyquist-Shannon sampling limit. The maximum axial extent attempted thus far was 110 µm along z, which leads to approximately one giga-voxel of complex sampling points per acquisition if imaging over the total microscope FOV (2.2 × 2.2 mm). We demonstrated quantitative measurement of the real and imaginary terms of the complex index of refraction within thick biological specimens. We believe that FPT can be significantly improved with additional experimental development. First, an improved LED array geometry will enable a higher angle of illumination to improve resolution. Second, we set the number of captured images here to match previously determined data redundancy requirements [50]. However, we have observed that reconstructions are successful with much fewer images than otherwise expected. Along with using a multiplexed illumination strategy [51], this may help significantly speed up tomogram capture time. Third, we set our reconstruction range along the z-axis somewhat arbitrarily at 110 µm. We expect the ability to further extend this axial range in future setups. Alternative computational approaches may also improve FPT. Here, we list a number of possible 134 (a) 3D FPT! Δz = -3.7 µm! Δz = 0! Δz = +7.4 µm! .04! Re(Δn)! 10 µm ! -.04! Protocoel! 3 µm ! Cell wall! Hydropore! Feature 1! Feature 2! (b) Manual refocus! Δz = -3.7 µm! Δz = 0! Δz = +7.4 µm! 1! Intensity (AU)! 10 µm ! 10 µm ! 3 µm ! 0! 3 µm ! Figure 8.7: 3D reconstruction of a starfish embryo at the larvae stage. (a) Three different axial planes of the FPT tomogram show significantly different features within the larvae. For example, the protocoel is completely missing from the ∆z = −3.7 µm plane. Likewise, the developmental “feature 1” (see lower right) is only visible in the ∆z = −3.7 µm plane. (b) This type of axial information, and even certain structures (like feature 1 and feature 2, marked in (a)) are completely missing from standard microscope images after manually refocusing to each plane of interest. 135 directions. First, we adopted the well-known alternating projections method (i.e., the ePIE algorithm [52]) for ptychographic update. Other methods, such as convex-based approaches [53], offer better performance in the presence of noise. Second, alternative approximations are also available to solve the first Born approximation [54]. Third, a big detriment to resolution is currently the missing cone in 3D k-space, and various methods are available to fill this cone in, e.g., by assuming the sample is positive-only, sparse, or of a finite spatial support [48, 49]. Finally, there are already suggested methods to solve for the full Born series, which take into account the effects of multiple scattering [43]. Connections between this type of multiple scattering solver, recent methods applying the multi-slice approximation [34,55], and FPT may lead to successful reconstruction of increasingly turbid samples. 136 Bibliography [1] V. Ntziachristos, “Going deeper than microscopy: the optical imaging frontier in biology,” Nature Methods 7, 603–611 (2010). [2] K. Chung et al., “Structural and molecular interrogation of intact biological systems,” Nature 497, 332–337 (2013). [3] F. Chen, P. W. Tillberg and E. S. Boyden, “Expansion microscopy,” Science 347(6221), 543-548 (2015). [4] D. A. Agard, “Optical sectioning microscopy: cellular architecture in three dimensions,” Rev. Biophys. Bioeng. 13, 191–219 (1984). [5] M. Levoy, Z. Zhang and I. McDowall, “Recording and controlling the 4D light field in a microscope,” J. of. Microscopy 435, 144-162 (2009). [6] M. Broxton et al., “Wave optics theory and 3D deconvolution for the light field microscope,” Opt. Express 21(21), 25418–25439 (2013). [7] S. R. P. Pavani et al., “Three-dimensional single-molecule fluorescence imaging beyond the diffraction limit using a double-helix point spread function,” Proc. Nat. Acad. Sci. 106, 2995 (2009). [8] A. Dubois, L. Vabre, A. C. Boccara, and E. Beaurepaire, “High-resolution full-field optical coherence tomography with a Linnik microscope,” Appl. Opt. 41, 805–812 (2002). [9] S. G. Adie et al., “Computational adaptive optics for broadband optical interferometric tomography of biological tissue,” Proc. Nat. Acad. Sci. 109, 7175–7180 (2012). [10] S. Kang et al., “Imaging deep within a scattering medium using collective accumulation of single-scattered waves,” Nature Photon. 9, 253–258 (2015). [11] T. E. Matthews et al., “Deep tissue imaging using spectroscopic analysis of multiply scattered light,” Optica 1, 105–111 (2014). 137 [12] N. Streibl, “Three-dimensional imaging by a microscope,” J. Opt. Soc. Am. A 2(2), 121-127 (1985). [13] E. Wolf, “Three-dimensional structure determination of semi-transparent objects from holographic data,” Opt. Comm. 1(4), 153–156 (1969). [14] V. Lauer, “New approach to optical diffraction tomography yielding a vector equation of diffraction tomography and a novel tomographic microscope,” J. Microscopy 205, 165–176 (2002). [15] Choi et al., “Tomographic phase microscopy,” Nature Methods 4(9), 717-719 (2007). [16] M. Debailleul et al., “Holographic microscopy and diffractive microtomography of transparent samples,” Meas. Sci. Technol. 19 074009 (2008). [17] M. Debailleul et al., “High-resolution three-dimensional tomographic diffractive microscopy of transparent inorganic and biological samples,” Opt. Lett. 24, 79–81 (2009). [18] Y. Sung et al., “Optical diffraction tomography for high resolution live cell imaging,” Opt. Express 1(17), 266–277 (2009). [19] A. Chai, M. Moscoso and G. Papanicolaou, “Array imaging using intensity-only measurements,” Inv. Prob. 27 015005 (2011). [20] T. D. Gerke and R. Piestun, “Aperiodic volume optics,” Nature Photon. 10, 1–6 (2010). [21] T. Kim et al., “White light diffraction tomography of unlabeled live cells,” Nature Photon. 8, 256–263 (2014). [22] K. Kim et al., “Diffraction optical tomography using a quantitative phase imaging unit,” Opt. Lett. 39(24), 6935-6938 (2014). [23] M. H. Maleki and A. J. Devaney, “Phase-retrieval and intensity-only reconstruction algorithms for optical diffraction tomography,” J. Opt. Soc. Am. A 10(5), 1086–1092 (1993). [24] T. C. Wedberg. and J. J. Stamnes, “Comparison of phase retrieval methods for optical diffraction tomography,” Pure Appl. Opt. 4, 39–54 (1995). [25] G. Gbur and E. Wolf, “Diffraction tomography without phase information,” Opt. Lett. 27(21), 1890–1892 (2002). [26] M. A. Anastasio, D. Shi, Y. Huang and G. Gbur, “Image reconstruction in spherical-wave intensity diffraction tomography,” J. Opt. Soc. Am. A 22(12), 2651–2661 (2005). [27] Y. Huang and M. A. Anastasio, “Statistically principled use of in-line measurements in intensity diffraction tomography,” J. Opt. Soc. Am. A 24(3), 626–642 (2007) 138 [28] S. O. Isikman et al., “Lens-free optical tomographic microscope with a large imaging volume on a chip,” Proc. Nat. Acad. Sci. 108(18), 7296–7301 (2011). [29] C. Zuo, J. Sun. J. Zhang, Y. Hu and Q. Chen, “Lensless phase microscopy and diffraction tomography with multi-angle and multi-wavelength illuminations using a LED matrix,” Opt. Express 23(11), 14314–14328 (2015). [30] T. E. Gureyev et al., “Phase-and-amplitude computer tomography,” Appl. Phys. Lett. 89, 034102 (2006). [31] A. V. Bronnikov, “Theory of quantitative phase-contrast computed tomography,” J. Opt. Soc. Am. A 19(3), 472–480 (2002). [32] G. Zheng, R. Horstmeyer and C. Yang, “Wide-field, high-resolution Fourier ptychographic microscopy,” Nat. Photonics 7, 739–745 (2013). [33] X. Ou, R. Horstmeyer, G. Zheng and C. Yang, “Quantitative phase imaging via Fourier ptychographic microscopy,” Opt. Lett. 38, 4845–4848 (2013). [34] L. Tian and L. Waller, “3D intensity and phase imaging from light field measurements in an LED array microscope,” Optica 2(2), 104-111 (2015). [35] A. M. Maiden, M. J. Humphry and J. M. Rodenburg, “Ptychographic trans- mission microscopy in three dimensions using a multi-slice approach,” J. Opt. Soc. Am. A 29(8), 1606–1614 (2012). [36] T. M. Godden, R. Suman, M. J. Humphry, J. M. Rodenburg and A. M. Maiden, “Ptychographic microscope for three-dimensional imaging,” Opt. Express 22(10), 12513–12523 (2014). [37] J. V. Roey, J. V. Donk, and P. E. Lagasse, “Beam-propagation method: analysis and assessment,” J. Opt. Soc. Am. 71(7), 803–810 (1981). [38] M. Dierolf et al., “Ptychographic X-ray computed tomography at the nanoscale,” Nature 467, 436-440 (2010). [39] I. Peterson et al., “Nanoscale Fresnel coherent diffraction imaging tomography using ptychography,” Opt. Express 20(22), 24678–24685 (2012). [40] K. Shimomura, A. Suzuki, M. Hirose and Y. Takahashi, “Precession x-ray ptychography with multislice approach,” Phys. Rev. B 91, 214114 (2015). [41] S. O. Hruszkewycz, M. V. Holt, A. Tripathi, J. Maser and P. H. Fuoss, “Framework for threedimensional coherent diffraction imaging by focused beam x-ray Bragg ptychography,” Opt. Lett. 36(12), 2227-2229 (2011). 139 [42] C. T. Putkunz et al., “Fresnel coherent diffraction tomography,” Opt. Express 18(11), 11746– 11753 (2010). [43] Y. M. Wang and W. C. Chew, “An Iterative Solution of the Two-Dimensional Electromagnetic Inverse Scattering Problem,” Int. J. Imag. Syst. Tech. 1(1), 100–108 (1989). [44] O. Haeberle, K. Belkebir, H. Giovaninni and A. Sentenac, “Tomographic diffractive microscopy: basics, techniques and perspectives,” J. Mod. Opt. 57(9), 686–699 (2010). [45] X. Ou, R. Horstmeyer, G. Zheng and C. Yang, “High numerical aperture Fourier ptychography: principle, implementation and characterization,” Opt. Express 23(3), 3472–3491 (2015). [46] S. Marchesini, “A unified evaluation of iterative projection algorithms for phase retrieval,” Rev. Sci. Instrum. 78, 011301 (2007). [47] X. Ou, G. Zheng and C. Yang, “Embedded pupil function recovery for Fourier ptychographic microscopy,” Opt. Express 22(5), 4960–4972 (2014). [48] K. C. Tam and V. Perezmendez, “Tomographical imaging with limited-angle input,” J. Opt. Soc. Am. 71, 582–592 (1981). [49] B. P. Medoff, W. R. Brody, M. Nassi, and A. Macovski, “Iterative convolution backprojection algorithms for image-reconstruction from limited data,” J. Opt. Soc. Am. 73, 1493–1500 (1983). [50] O. Bunk, M. Dierolf, S. Kynde, I. Johnson, O. Marti and F. Pfeiffer, “Influence of the overlap parameter on the convergence of the ptychographical iterative engine,” Ultramicroscopy 108, 481–487 (2008). [51] L. Tian, X. Li, K. Ramchandran and L. Waller, “Multiplexed coded illumination for Fourier Ptychography with an LED array microscope,” Biomed. Opt. Express 19, 2376–2389 (2014). [52] A. M. Maiden and J. M. Rodenburg, “An improved ptychographical phase retrieval algorithm for diffractive imaging,” Ultramicroscopy 109, 1256–1262 (2009). [53] R. Horstmeyer, R. C. Chen, X. Ou, B. Ames, J. A. Tropp and C. Yang, “Solving ptychography with a convex relaxation,” New. J. Phys. 17 053044 (2015). [54] A. J. Devaney, ”Inverse-scattering theory within the Rytov approximation,” Opt. Lett. 6, 374– 376 (1981). [55] U. S. Kamilov et al., “Learning approach to optical tomography,” Optica 2(6), 517–522 (2015).