Computational Imaging: a Quest for the Perfect Image Citation Chung, Jaebum (2019) Computational Imaging: a Quest for the Perfect Image. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/8W3A-HE02. https://resolver.caltech.edu/CaltechTHESIS:05202019-151055724 Abstract A physical lens is limited in its ability to capture an image that is both high- resolution and wide-field due to aberrations even with a sophisticated lens design. This thesis explores computational methods that expand on the recently developed Fourier ptychographic microscopy (FPM) to overcome the physical limitations. New algorithms and imaging methods extend the computational aberration correction to more general imaging modalities including fluorescence microscopy and incoherent bright-field imaging so as to allow even a crude lens to perform like an ideal lens. This paradigm shift from the lens design to computational algorithms democratizes high-resolution imaging by making it easier to use and less complicated to build. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: Fourier ptychographic microscopy, FPM, computational aberration correction, deconvolution, phase retrieval, retinal imaging, wide field of view imaging Degree Grantor: California Institute of Technology Division: Engineering and Applied Science Major Option: Electrical Engineering Thesis Availability: Public (worldwide access) Research Advisor(s): Yang, Changhuei Thesis Committee: Wang, Lihong (chair) Yang, Changhuei Shapiro, Mikhail G. Hassibi, Babak Zheng, Guoan Defense Date: 1 May 2019 Non-Caltech Author Email: jaebum (AT) alumni.caltech.edu Funders: Funding Agency Grant Number NIH R21 EY026228A NIH R01 AI096226 Caltech Innovation Initiative (CI2) Research Grant. 25570015 NIH 1DP2OD007307 Department of Defense (DOD) W81XWH-09-1-0050 Record Number: CaltechTHESIS:05202019-151055724 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:05202019-151055724 DOI: 10.7907/8W3A-HE02 Related URLs: URL URL Type Description https://doi.org/10.1117/1.JBO.19.6.066007 DOI Article adapted for ch. 3 https://doi.org/10.1364/BOE.7.000352 DOI Article adapted for ch. 4 https://doi.org/10.1364/OPTICA.6.000647 DOI Article adapted for ch. 5 ORCID: Author ORCID Chung, Jaebum 0000-0003-3932-8428 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 11525 Collection: CaltechTHESIS Deposited By: Jaebum Chung Deposited On: 23 May 2019 16:09 Last Modified: 08 Nov 2023 18:41 Thesis Files Preview PDF - Final Version See Usage Policy. 35MB Repository Staff Only: item control page

Computational Imaging: a Quest for the Perfect Image Thesis by Jaebum Chung In Partial Fulfillment of the Requirements for the Degree of Electrical Engineering CALIFORNIA INSTITUTE OF TECHNOLOGY Pasadena, California 2019 Defended May 1, 2019 ii © 2019 Jaebum Chung ORCID: 0000-0003-3932-8428 All rights reserved except where otherwise noted iii ACKNOWLEDGEMENTS I owe my Ph.D. to my advisor, Changhuei Yang, for taking a leap of faith to accept me into his lab even when I did not come from an optics background and could barely grasp what numerical aperture was. I would not know what I know today if it weren’t for his belief in me. I also thank him for setting high standards of work ethic while also providing me with the freedom to pursue an internship opportunity. This pushed me beyond my comfort zone to become a better scientist. The members of our Biophotonics Lab throughout the years have also been essential to my development as a scientist. I would like to give big thanks to Roarke Horstmeyer and Xiaoze Ou for being the sages of Fourier ptychography and beyond to always humble me with their knowledge; Edward Haojiang Zhou, Mooseok Jang, Haowen Ruan, and Atsushi Shibukawa for being a great resource with their extensive knowledge in experimental optics; Joshua Brake for helping me with a broad set of engineering know-hows that extends beyond optics; Jinho Kim and Hangwen Lu for being reliable and insightful research partners; and Daniel Martin, Antony Chan, An Pan, Xiaoyu Lyu, Seung Ah Lee, Benjamin Judkewitz, Jiangtao Huangfu, Chao Han, Yan Liu, Michelle Cua, Jian Xu, Yujia Huang, Changsoon Choi, Cheng Shen, Ruizhi Cao, Ali Kazemi Jahromi, and Baptiste Blochet for inspiring discussions and collaborations. Last but not least, I would like to show my deep appreciation for Anne Sullivan for streamlining the research experience in the lab so that I can do my best work. Most of my work would not have been possible without collaborators outside of the lab: Anthony Williams, Liheng Bian, Srinivas Sadda, Amir Hariri, Karen Lencioni, and Gloria Martinez. A special shoutout to Gloria for letting me torture her weekly with monkeys. I also received guidance from Guoan Zheng, Babak Hassibi, Mikhail Shapiro, Lihong Wang, Amnon Yariv, Christy Jenstad, Alon Greenbaum, Soo-Young Kim, Dierck Hillmann, Martha Neuringer, Laurie Renner, Michael Andrews, and Mark Pennesi. I was also lucky to have known Evan Miyazono who helped me connect with students across various labs at Caltech in the beginning of my Ph.D. career. I cannot imagine my life as a Ph.D. student without Seoyeon Lee, Emily Yun, Donghun Ryu, Serim Ryou, Gahye Jeong, Dong-Gyu Lee, Ryo Adachi, Yu and Chika Horie, Edison Huang, Maiko Obana, Koichiro and Mai Kajikawa, Andrew Cho, iv Tomoyuki Onyama, Sonjong Huang, Peter Lee, Yapeng Su, Xiaozhou Ruan, Qifan Yang, Dongwan Kim, and friends across the globe. I thank them for all the diverse experiences we shared together that made the everyday Ph.D. life unforgettable. Finally, I would like to thank my family and Steph for being the permanent light in my life to guide me through moments of darkness, and cheering me on unconditionally to be the best version of myself. v ABSTRACT A physical lens is limited in its ability to capture an image that is both highresolution and wide-field due to aberrations even with a sophisticated lens design. This thesis explores computational methods that expand on the recently developed Fourier ptychographic microscopy (FPM) to overcome the physical limitations. New algorithms and imaging methods extend the computational aberration correction to more general imaging modalities including fluorescence microscopy and incoherent bright-field imaging so as to allow even a crude lens to perform like an ideal lens. This paradigm shift from the lens design to computational algorithms democratizes high-resolution imaging by making it easier to use and less complicated to build. vi PUBLISHED CONTENT AND CONTRIBUTIONS [1] L. Bian, J. Suo, J. Chung, X. Ou, C. Yang, F. Chen, and Q. Dai. Fourier ptychographic reconstruction using Poisson maximum likelihood and truncated Wirtinger gradient. Nature Scientific Reports, 6:27384, 2016. doi: 10.1364/OPTICA.3.000827. J. Chung participated in the conception of the project, data acquisition, and manuscript writing. [2] J. Chung, X. Ou, R. P. Kulkarni, and C. Yang. Counting White Blood Cells from a Blood Smear Using Fourier Ptychographic Microscopy. PLoS ONE, 10(7):e0133489, 2015. doi: 10.1371/journal.pone.0133489. J. Chung designed and carried out the experiment, processed and analyzed the data, and wrote the manuscript. [3] J. Chung, J. Kim, X. Ou, R. Horstmeyer, and C. Yang. Wide field-of-view fluorescence image deconvolution with aberration-estimation from Fourier ptychography. Biomedical Optics Express, 7(2):352–368, 2016. doi: 10.1364/ BOE.7.000352. J. Chung participated in the conception of the project, designed and carried out the experiment, processed and analyzed the data, and wrote the manuscript. [4] J. Chung, H. Lu, X. Ou, H. Zhou, and C. Yang. Wide-field Fourier ptychographic microscopy using laser illumination source. Biomedical Optics Express, 7(11):4787–4802, 2016. doi: 10.1364/BOE.7.004787. J. Chung participated in the conception of the project, designed and carried out the experiment, processed and analyzed the data, and wrote the manuscript. [5] J. Chung, G. W. Martinez, K. C. Lencioni, S. R. Sadda, and C. Yang. Computational aberration compensation by coded-aperture-based correction of aberration obtained from optical Fourier coding and blur estimation (CACAO-FB). Optica (Accepted), 2019. J. Chung participated in the conception of the project, designed and carried out the experiment, processed and analyzed the data, and wrote the manuscript. [6] R. Horstmeyer, X. Ou, J. Chung, G. Zheng, and C. Yang. Overlapped Fourier coding for optical aberration removal. Optics Express, 22(20):24062–24080, 2014. doi: 10.1364/OE.22.024062. J. Chung participated in processing and analyzing the data, and writing the manuscript. [7] R. Horstmeyer, J. Chung, X. Ou, G. Zheng, and C. Yang. Diffraction tomography with Fourier ptychography. Optica, 3(8):827–835, 2016. doi: 10.1364/OPTICA.3.000827. vii J. Chung participated in the conception of the project, data acquisition and analysis, and manuscript writing. [8] H. Lu, J. Chung, X. Ou, and C. Yang. Quantitative phase imaging and complex field reconstruction by pupil modulation differential phase contrast. Optics Express, 24(22), 2016. doi: 10.1364/OE.24.025345. J. Chung participated in the conception of the project, design of the experiment, and manuscript writing. [9] X. Ou, J. Chung, R. Horstmeyer, and C. Yang. Aperture scanning Fourier ptychographic microscopy. Biomedical Optics Express, 7(8):3140–3150, 2016. doi: 10.1364/BOE.7.003140. J. Chung participated in the conception of the project, data analysis, and manuscript writing. [10] A. Williams, J. Chung, X. Ou, G. Zheng, S. Rawal, Z. Ao, R. Datar, C. Yang, and R. Cote. Fourier ptychographic microscopy for filtration-based circulating tumor cell enumeration and analysis. Journal of Biomedical optics, 19(6): 066007, 2014. doi: 10.1117/1.JBO.19.6.066007. J. Chung participated in the experiment design, carried out the experiment, processed and analyzed the data, and wrote the manuscript. [11] G. Zheng, X. Ou, R. Horstmeyer, J. Chung, and C. Yang. Fourier Ptychographic Microscopy: A Gigapixel Superscope for Biomedicine. Optics and Photonics News, OSA, 25:28–33, 2014. J. Chung participated in writing the manuscript. viii TABLE OF CONTENTS Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v Published Content and Contributions . . . . . . . . . . . . . . . . . . . . . . vi Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Illustrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix Chapter I: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Characteristics of a lens . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 A standard microscope system . . . . . . . . . . . . . . . . . . . . . 4 1.3 Limitations in imaging with physical lenses . . . . . . . . . . . . . . 7 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Chapter II: Fourier ptychographic microscopy . . . . . . . . . . . . . . . . . 12 2.1 Basic principles of FPM . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 Advantages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 Chapter III: High-resolution, wide-field imaging of samples with uneven surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.3 Analysis of CTC images by FPM and comparison with standard microscopy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 Chapter IV: Correcting for aberrations in fluorescence image with FPMrecovered pupil function . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.2 Principle of operation . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.3 Experimental demonstration . . . . . . . . . . . . . . . . . . . . . . 35 4.4 Quantifying the adequate amount of fluorescence signal for deconvolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.6 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 Chapter V: Correcting for aberrations in a general optical system . . . . . . . 54 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.2 Principle of CACAO-FB . . . . . . . . . . . . . . . . . . . . . . . . 56 5.3 Experimental demonstration on artificial samples . . . . . . . . . . . 69 ix 5.4 Adapting CACAO-FB to an in vivo experiment on the eye of a rhesus macaque . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 77 79 84 x LIST OF ILLUSTRATIONS Number Page 1.1 A beam of light can be focused by refraction, reflection, or diffraction. The cross-sectional views of a lens, a GRIN lens, a concave mirror, and a Fresnel zone plate are illustrated as examples. . . . . . . . . . . 2 1.2 A 4f system. The space between the objective lens and the tube lens is referred to as the “infinity space” because a point in the sample plane or a detector plane is mapped to a collimated beam in this region. For a typical microscope, the infinity space may not necessarily be fobj + ftube as it does not have any influence on the intensity image formation on the detector. . . . . . . . . . . . . . . . . . . . . . . . 5 2.1 Experimental setup of FPM. a) An LED matrix is fitted to a standard microscope system, replacing its condenser. An LED matrix placed sufficiently far away (∼8 cm) from the sample provides angularly varying planar illumination. The objective lens’s backfocal plane (i.e., Fourier plane) is indicated by the yellow box. b) Different illumination angles correspond to sampling different spatial frequency information of the sample. c) LEDs are lit up one by one, and the pitch of the LEDs should be small enough to satisfy the overlap requirement of the spatial spectrum regions in the Fourier space. . . . 13 2.2 The embedded pupil function recovery (EPRY) algorithm. EPRY is a phase retrieval process that can synthesize the FPM’s intensity dataset into a complex, high SBP image and simultaneously reconstruct the pupil function of the imaging system. . . . . . . . . . . . . . . . . . 15 3.1 (a) The FPM setup consists of (from the bottom) an LED matrix for sample illumination, a microscope system with a 2x objective, and a camera connected to a computer. (b) The Fourier spectrum point of view. The center red subregion corresponds to the spatial frequency of the low-resolution image captured with plane waves with k x = k y = 0. The off-center red subregion correlates to an oblique angle illumination with wavevector (k x, k y ). (c) Light from an LED at an oblique angle corresponds to a plane wave with a k vector (k x, k y ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 xi 3.2 (a) Full field-of-view color image of the entire microfilter containing captured tumor cells by FPM. Magnified FPM images (b-d1) selected from different areas of the microfilter show detailed morphology of tumor cells, where all sections are well in focus because of the automatic EPRY-FPM program. (d2) A standard microscope (with 40x objective) image shows the corresponding region to (d1), but because of the uneven surface of the microfilter, subregions (d2-1) and (d2-2) cannot be focused simultaneously. Also, its field of view is limited when compared to (d2), as seen by the aperture’s outline at the edge of the image, in contrast to FPM’s wide field-of-view that can provide high-resolution images of the entire micro- filtration area. 3.3 (a) A microfilter’s surface profile characterized by the pupil function recovered by EPRY. The focal plane of the objective is at 150 µm. The maximum difference in-depth across the filter is about 250 µm, which is within FPM’s refocusing capacity of 300 µm. In this case, the captured microfilter image is sectioned into 17 x 17 tiles, and EPRY iteratively characterizes each tile’s defocus level in its highresolution image reconstruction process. (b1)-(c2) Small subregions are extracted from two different surface levels, showing before and after refocusing by EPRY. (b1) is already very close to the focal plane, so there is only a minor improvement after refocusing, as in (b2). (c1) is not in the focal plane and is blurry. (c2) shows the refocused result. 3.4 Graph demonstrating the correlation between tumor cell count by standard microscopy (y-axis) and tumor cell count by FPM (x-axis) in corresponding microfilter samples, where each data point represents a single trial with tumor cells enumerated by both methods. . . . . . 4.1 The experimental setup of joint FPM-fluorescence imaging. The 4f setup has a camera lens, an iris at the back focal plane, and a tube lens. The filter is included only for fluorescence imaging and FPM imaging of the same color channel. First, raw FPM images are captured using varied-angle illumination provided by the LED matrix. Then, a fluorescence image is captured with the illumination from the excitation LED. . . . . . . . . . . . . . . . . . . . . . . . . 26 27 28 36 xii 4.2 Experimental reconstruction of an improved-resolution amplitude and phase image, along with an aberration-corrected fluorescence image, using FPM. (a) Series of low-resolution coherent green images captured with variable-angle illumination from an LED matrix. (b) EPRY is applied on the low-resolution images to generate a highresolution, complex image of the sample, along with the characterization of the microscope’s pupil function. (c) The pupil function is converted into an incoherent PSF, which is the blur kernel induced to a fluorescence image by the imaging system in the same sample region. (d) The PSF is deconvolved from the fluorescence image using Eq. 4.6 to generate an aberration-corrected image. (e) The one-dimensional profile of two beads, in the raw and deconvolved images, shows improved contrast and position accuracy. (f) In the overlay of FPM and fluorescence images, the fluorescence signal is localized with good centroid accuracy after deconvolution. . . . . . . 38 4.3 Comparison between (a) a raw image of a specific sample ROI shifted 4.63 mm off the optical axis (left), its deconvolution result (right), and (b) a raw image captured with the same sample ROI centered on the optical axis, i m (x, y). The PSF recovered from FPM, hm (x, y), is deconvolved from the image to remove aberrations in (a). The result shows close resemblance to the image in (b), i center (x, y), which we assume is minimally impacted by system aberrations. . . . . . . . . . 39 xiii 4.4 Demonstration of FPM and aberration-corrected fluorescence imaging across different regions of a large image FOV. (a)-(c) correspond to regions labeled in the full FOV fluorescence image in (d). With FPM, we improve bright-field image resolution (1st and 2nd column) and characterize the spatially varying pupil functions (6th column). Hot pixels and chromatic aberrations appearing in the low-resolution color images are suppressed and corrected after FPM reconstruction. We correct aberrations in the fluorescence image using each pupil function (3rd and 4th column). The phase gradient images from FPM can be combined with fluorescence data to elucidate the structures and the locations of the nuclei of HeLa cells (5th column). The cell membrane morphology elucidated by the phase gradient can differentiate between cells undergoing cytokinesis (arrow in (b)) and cells in telophase (arrow in (c)), which is otherwise difficult to do in the fluorescence images alone. . . . . . . . . . . . . . . . . . . . . . 41 4.5 (a) Normalized mean square error (NMSE) of the raw, i m (u; t), and deconvolved image, i˜m (u; t), in Fig. 4.3(a) is plotted against the detector exposure time. Higher exposure is equivalent to higher SNR of the captured image. NMSE starts to plateau for both raw and deconvolved images after about 21 min of exposure (raw image SNR = 18.1) as indicated by the green broken line, with the deconvolved image’s NMSE = 0.0057 being 31% lower than the raw image’s NMSE = 0.0083. (b) Example images of the ROI used for this study. The reference image is generated by centering the ROI and capturing with our imaging system. We use the area inside the small yellow box to quantify the SNR. Raw images are captured with varying exposure time (45 sec and 21 min shown here) while the ROI is 4.63 mm away from the center. We use the same small area in each raw image to quantify its SNR. Deconvolved images are generated by applying Eq. 4.6 to the raw images. Longer exposure provides better SNR images and deconvolution results. . . . . . . . . . . . . . . . . . . . 44 xiv 4.6 The variations in experimental setup for comparing the two different pupil function recover methods. (a) Collimated light source illuminates the sample perpendicularly, and the sample is brought to different defocus planes in the defocus diversity-based pupil recovery method. (b) LEDs of varied illumination angles act as the light source in FPM setup. EPRY algorithm jointly solves for the sample spectrum and the pupil function using sample images captured under different LED illumination angles. . . . . . . . . . . . . . . . . . . . 4.7 Full FOV image of a sample of microspheres with associated pupil functions in 3 different sub-regions. The red box indicates the pupil function recovered by EPRY while the blue box stands for defocus diversity-based pupil function recovery method. The two pupil estimation methods correlate well throughout different regions of the FOV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Optical architecture of CACAO-FB. CACAO-FB system consists of 3 tube lenses (L1, L2, and L3) to relay the image from the target system for analysis. The target system consists of an unknown lens and an unknown sample with spatially incoherent field. CACAO-FB system has access to the conjugate plane of the target system’s pupil that can be arbitrarily modulated with binary patterns using a spatial light modulator. The images captured by CACAO-FB system is intensityonly. f0 , f1 , f2 , and f3 are the focal lengths of the unknown lens, L1, L2, and L3, respectively. d is an arbitrary distance smaller than f3 . . . 5.2 Outline of CACAO-FB pipeline. a) The captured images are broken into small tiles of isoplanatic patches (i.e., aberration is spatially invariant within each tile). b) Data acquisition and post-processing for estimating the pupil function, Pt (u, v). Limited-aperture images, i m,t (ξ, η), are captured with small masks, Wm (u, v), applied at the pupil plane. Local PSFs, bm,t (ξ, η), are determined by the blur estimation procedure, Algorithm 1. These PSFs are synthesized into the fullaperture pupil function, Pt (u, v), with Fourier-ptychography-based alternating projections algorithm, Algorithm 2. c) Data acquisition with big masks, An (u, v), at the pupil plane. d) The recovered Pt (u, v) from b) and the big-aperture images φn,t (ξ, η) from c) are used for deconvolution, Algorithm 3, to recover the latent aberration-free intensity distribution of the sample, ot (x, y). . . . . . . . . . . . . . . . 48 49 57 59 xv 5.3 Simulating image acquisition with different small masks at the pupil plane. a) The full pupil function masked by the lens’s NA-limited aperture. Differently masked regions of the pupil, b1-3), give rise to different blur kernels, c1-3), which allows us to capture images of the sample under the influence of different PSFs. Only the phase is plotted for Pt (u, v) and Pm,t (u, v)’s, and their apertures are marked by the black boundaries. W1 (u, v), W45 (u, v), and W52 (u, v) are 3 small masks from a spiraling-out scanning sequence. . . . . . . . . . . . . 5.4 Flowchart of Algorithm 1: blur estimation algorithm for determining local PSFs from images captured with small apertures, Wm,t (u, v) . . . 5.5 Flowchart of Algorithm 2: Fourier-ptychography-based alternating projections algorithm for reconstructing the unknown lens’s pupil function, Pt (u, v) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Simulation of our pupil function recovery procedure and a comparison with blind deconvolution algorithms. a) The Siemens star pattern used in the simulation. b) The system’s pupil function and the associated PSF. c) A series of images, i m,t (ξ, η)’s captured with small masks, Wm (u, v), applied to the pupil function. d) An image captured with the full-pupil-sized mask, An (u, v), on the pupil function, which simulates the general imaging scenario by an aberrated imaging system. e) The system’s pupil function and PSF recovered by our procedure. They show high fidelity to the original functions in b). f) Blur functions recovered by MATLAB’s and Fergus et al.’s blind deconvolution algorithm, respectively. They both show poor reconstructions compared to the recovered PSF in e). . . . . . . . . . 62 62 64 65 xvi 5.7 Simulation that demonstrates the benefit of coded-aperture-based deconvolution. a1-5) Masked pupil functions obtained by masking the same pupil function with the full circular aperture and coded apertures under different rotation angles (0◦ , 45◦ , 90◦ , 135◦ ), their associated OTFs along one spatial frequency axis, and captured images. Each coded aperture is able to shift the null regions of the OTF to different locations. b) Comparison between the OTF of a circular-aperture-masked pupil function and the summed OTFs of the circular- and coded-aperture-masked pupil functions. Null regions in the frequency spectrum are mitigated in the summed OTF, which allows all the frequency content of the sample within the bandlimit to be captured with the imaging system. The OTF of an ideal pupil function is also plotted. c1) Deconvolved image with only a circular aperture shows poor recovery with artifacts corresponding to the missing frequency contents in the OTF’s null regions. c2) A recovered image using one coded aperture only. Reconstruction is better than c1), but still has some artifacts. c3) A recovered image using circular and multiple coded apertures is free of artifacts since it does not have missing frequency contents. . . . . . . . . . . . . . . 67 5.8 Flowchart of Algorithm 3: iterative Tikhonov regularization for recovering the latent sample image, ot (x, y), from aberrated images. Here, Φn,t (u, v) = F{φn,t (ξ, η)}(u, v). . . . . . . . . . . . . . . . . . . 68 5.9 Experimental setup of imaging a sample with a crude lens (i.e., unknown lens). Sample is illuminated by a monochromatic LED (520 nm), and the lens’ surface is imaged onto the SLM by a 1:1 lens relay. The part of light modulated by the SLM is reflected by the PBS and is further filtered by a polarizer to account for PBS’s low extinction ratio in reflection (1:20). Pupil-modulated image of the sample is captured on the sCMOS camera. L: lens; P: polarizer; PBS: polarizing beam splitter. . . . . . . . . . . . . . . . . . . . . . 70 xvii 5.10 Resolution performance measured by imaging a Siemens star target. a) A crude lens has optical aberration that prevents resolving the Siemens star’s features. b) CACAO-FB is able to computationally remove the aberration and resolve 19.6 µm periodicity feature size, which lies between the coherent and incoherent resolution limit given by the focal length of 130 mm, the aperture diameter of 5.5 mm, and the illumination wavelength of 520 nm. c) Pupil function recovered by CACAO-FB used for removing the aberration. d) The PSF associated with the pupil function. e) Intensity values from the circular traces on a) and b) that correspond to the minimum resolvable feature size of 19.6 µm periodicity. The Siemens star’s spokes are not visible in the raw image’s trace, whereas 40 cycles are clearly resolvable in the deconvolved result’s trace. . . . . . . . . . . . . . . . . . . . . . . . 5.11 Spatially varying aberrations compensation result on a grid of USAF target. a) the full FOV captured by our camera with the full circular aperture at 5.5 mm displayed on the SLM. Each small region denoted by b-d1) had a different aberration map as indicated by varying pupil function and PSFs. Spatially varying aberrations are adequately compensated in post-processing as shown by the deconvolution results b-d2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.12 Eye model with a USAF target embedded on the retinal plane. a) A cut-out piece of glass of USAF target is attached on the retina of the eye model. The lid simulates the cornea and also houses a lens element behind it. b) The model is filled with water with no air bubbles in its optical path. c) The water-filled model is secured by screwing it in its case. . . . . . . . . . . . . . . . . . . . . . . . . . 5.13 Experimental setup of imaging an eye model and an in vivo eye. Illumination is provided by a fiber-coupled laser diode (520 nm), and the eye’s pupil is imaged onto the SLM by a 1:1 lens relay. The sample is slightly defocused from the focal length of the crude lens to add additional aberration into the system. Pupil alignment camera provides fiduciary to the user for adequate alignment of the pupil on the SLM. PBS2 helps with removing corneal reflection. Motion-reference camera is synchronized with encoded-image camera to capture images not modulated by the SLM. BS: beam splitter; L: lens; M: mirror; P: polarizer; PBS: polarized beam splitter; QWP: quarter-wave plate. . . 71 72 73 74 xviii 5.14 5.15 5.16 5.17 5.18 5.19 CACAO-FB result of imaging the USAF target in the eye model. a) Raw image (2560 x 1080 pixels) averaged over 12 frames captured with the full circular aperture at 4.5 mm. The pupil function and PSF in each boxed region show the spatially varying aberrations. b-d) Deconvolution results show sharp features of the USAF target. The uneven background is from the rough surface of the eye model’s retina. 75 Showing the importance of masked pupil kernel shape determination for successful deconvolution. a1-3) limited PSFs determined only by considering their centroids. b) Recovered aberration and deconvolution result obtained with centroid-only limited PSFs. Some features of USAF are distorted. c1-3) limited PSFs determined with the blur estimation algorithm. d) Recovered aberration and deconvolution result obtained with the blur-estimated local PSFs. No distortions in the image are present, and more features of the USAF target are resolved. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 CACAO-FB result from imaging an in vivo eye of a rhesus macaque. a) Raw image averaged over 213 frames captured with 4.5 mm full circular aperture. b) Deconvolution result using the c) pupil function reconstructed by CACAO-FB procedure. d) PSF associated with the pupil function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 The influence of the limited masks’ overlap ratio on the pupil function recovery. The recovery becomes poor below 30% overlap. . . . . . . 82 Motion-reference camera images registered for rotation and translation. Only 4 representative frames are shown. . . . . . . . . . . . . . 83 Raw full aperture image. a) 1 frame and b) a sum of 213 frames. . . . 83 xix LIST OF TABLES Number Page 1.1 The resolution and the space-bandwidth product of select Olympus microscope objective lenses. Resolution here is given by Rayleigh criterion, r = 0.61λ NA , where r is the distance from the center of a diffraction limited spot to its first zero for a circular aperture with a given NA. SBP is calculated by dividing the area of the circular FOV (diameter is d = FN/M) by (r/2)2 , which is the pixel area that meets the sampling requirement. The SBP remains constant for low magnification with low NA, but it drastically decreases for high NA objectives due to the difficulty of accounting for increasing aberrations. 8 1 Chapter 1 INTRODUCTION An imaging device allows us to capture and preserve a moment in time. Not only is it a device of sentimental and historical values, but it is also an invaluable tool in scientific discovery, medical diagnosis, military, and infrastructure, as it allows us to see what would otherwise be hard or impossible to observe with the naked eye. To name a few examples, it allows us to observe ever smaller things (e.g., microscopy); see things that are invisible due to their transparency (e.g., phase imaging, infra-red imaging) or they are obscured by another object (e.g., CT imaging, optical phase conjugation); identify chemical structures of material (e.g., fluorescence microscopy, X-ray diffraction); see farther (e.g., astronomy) and identify depth (e.g., LIDAR). Optical imaging systems invariably use lenses to form images. The earliest form of lens dates all the way back to 750 BC [1]. A lens consists of a transparent material with a curved surface that causes the incident light to refract and converge or diverge as it is transmitted through it. Forming an image with lenses refers to using lenses to refract the light originating from points on an object of interest and relay it to points on our detector, which could be a sensor in a camera or the retinal layer of our eye. A lens can also take the form of a flat piece of transparent material whose light propagation speed varies with location (e.g., a gradient-index lens or a spatial light modulator). Some other ways of focusing light include: 1) reflection, where the beam is reflected into a point by a concave mirror; and 2) diffraction, where the wave nature of light is used to diffract the beam into a point. These various focusing methods are illustrated in Fig. 1.1. In this thesis, we are concerned with lenses, which focus light using refraction. In essence, a lens applies different phase delays on the incident light’s wavefront to influence the light’s propagation. A microscope’s objective lens showcases the sophisticated lens design involved to allow for tight focusing over a wide area. A single lens may be sufficient to form a tight focus at a single location the focal distance away from the lens, but its focusing ability degrades for other lateral regions. Multiple lens elements are incorporated into an objective lens to overcome this limitation and allow for tightly focused spots across a large area [9]. The product of the number of distinguishable focus spots and 2 Figure 1.1: A beam of light can be focused by refraction, reflection, or diffraction. The cross-sectional views of a lens, a GRIN lens, a concave mirror, and a Fresnel zone plate are illustrated as examples. the field of view (FOV) is usually referred to as the space-bandwidth product (SBP). There is usually a limit on the maximum achievable SBP by the sophisticated lens design since adding more optical surfaces to reduce aberrations over a larger FOV has diminishing returns. We elaborate on this further in the following sub-sections. Computational imaging methods, riding on the exponentially increasing computation power, are enabling capabilities that were not possible with physical lenses alone. For example, light field cameras [7] and coded aperture photography [6, 15] allow a user to refocus after images have been taken. In ophthalmology, a clear image of retina at the back of the eye can be formed with optimization algorithms [4, 5, 12], overcoming the inherent eye’s aberration computationally. There are imaging modalities that completely bypass use of lens and only use computation to perform imaging (i.e., lensless imaging [11, 13]), but the physical lens still remains the most robust and broadly applicable way to capture images. Not only can leveraging computational methods lead to a simpler setup that can be easily miniaturized and allow for more flexibility (e.g., powerful cellphone camera [8]), it can also be beneficial for certain imaging modalities like fluorescence where phototoxicity and bleaching is of real experimental constraint. This thesis focuses on harnessing and improving a recently developed computationally imaging method, called Fourier ptychographic microscopy (FPM), to surpass the challenges of physical optics. FPM shifts the paradigm of lens design and imaging by accounting for them in the computational realm instead of making sure the perfect image is taken at the moment of acquisition. Not only does it use simple optical components and computational algorithms to achieve much greater SBP than what is allowed conventionally, it also simultaneously corrects for the inherent aber- 3 rations in the captured image due to the imperfections in the system or the imaging scenario. We will detail demonstrations of how we apply this method to solve a common problem of biological specimen imaging, incorporate it into fluorescence imaging to obtain clearer images, and how we are able to modify the optical setup and algorithm further to make the aberration correction method more widely applicable to general optical systems. To describe further, we first need to familiarize ourselves with several concepts and terminologies in the following section. We review the properties of a lens as it constitutes the fundamental building block of an optical imaging system. 1.1 Characteristics of a lens We consider a lens made of a transparent medium, such as glass and plastic. The lens’s medium has a different refractive index, nlens , compared to its surroundings, usually air with nair = 1. A plane wave incident on the lens will experience some phase delay as it propagates through the lens. A convex lens induces phase delays that focuses the incident beam into a point, whereas a concave lens makes the beam diverge. Focal length A convex lens focuses the beam to a point at a certain distance away from the lens. This is called the focal length. The lens’s ability to focus, i.e., focusing power, is higher with more convex surfaces of the lens, and it in turn shortens the focal length. Numerical aperture (NA) Numerical aperture is related to the refractive index of the medium, n, surrounding the lens and the angle, α, subtended by the steepest focused or the collected ray from the lens’s normal. The angular subtense can calculated by the lens’s effective aperture diameter and the focal length. NA is defined as NA = n sin α and it determines the extent of spatial frequency information relayed by the lens [2], as will be described in the following chapter. Fourier transformation For a convex lens of focal length f , the beam profile (i.e., complex field) incident on the lens is Fourier transformed at f away from the lens [2]. The field distribution at f before the lens is exactly Fourier transformed at f behind the lens, while the field at other locations before the lens is also Fourier transformed at f behind the lens albeit with additional phase terms. Thus, the plane at f behind the lens is also 4 referred to as the lens’s Fourier plane. This will be further explored in Chapter 2 and Chapter 5. This property allows us to directly access and manipulate the Fourier spectrum of the incident field, and it is extremely useful for image processing. 1.2 A standard microscope system The most commonly used imaging system in science is a microscope. The series of work illustrated in this thesis is based on an infinity-corrected microscope, which is the standard microscope in use today. In the following section, we outline the building blocks of a standard microscope. 4f arrangement An infinity-corrected microscope refers to a microscope system with its optical elements in a 4f arrangement. In such an arrangement, an objective lens is placed its focal length ( fobj ) away from the sample plane, followed by a tube lens placed at a distance corresponding to the sum of its focal length ( ftube ) and fobj from the objective, followed by a detector plane at ftube away from the tube lens. This is illustrated in Fig. 1.2. Because the sample plane is placed at fobj in front of the objective lens, one has direct access to the sample’s Fourier transform at the objective’s Fourier plane. Placing an element such as a spatial light modulator (SLM) at this plane allows one to directly manipulate the sample’s Fourier spectrum prior to capturing the sample’s image with a detector. Another key characteristic of a 4f system is the fact that a point on the sample plane or the detector plane is mapped to a collimated beam in the space between the objective and the tube lens, often referred to as the “infinity space.” A collimated beam incident on a lens is focused to a point on the lens’s focal plane whose lateral position is defined by the incidence angle, invariant to the incidence location on the lens as shown in Fig. 1.2. Therefore, the length of the infinity space can vary without degrading the image’s intensity formed on the detector plane. The image’s phase information in the case of a coherent imaging scheme would change, as will be described in Chapter 5. Therefore, it is a common practice to add non-refracting optical elements such as a polarizer and a color filter in this space, as the change in the optical path length introduced by these elements in the infinity space will not degrade the final image. Magnification (M) Magnification of a microscope is straight-forward for an infinity-corrected system, and is defined as M = ftube / fobj . As microscopes allows for visualizing small 5 Figure 1.2: A 4f system. The space between the objective lens and the tube lens is referred to as the “infinity space” because a point in the sample plane or a detector plane is mapped to a collimated beam in this region. For a typical microscope, the infinity space may not necessarily be fobj + ftube as it does not have any influence on the intensity image formation on the detector. features, the magnification is typically above 1x, which means the objective lens has a shorter focal length than that of the tube lens. Commercial objective lenses have magnification values associated with them, and they are fundamentally determined by the tube lens for which the manufacturers designed. Tube lenses’ focal lengths are different for different microscope manufacturers. Among the most widely used microscopes, Olympus has ftube = 180 mm while Nikon and Leica has ftube = 200 mm. The focal length of an objective for a given manufacturer can therefore be determined by fobj = ftube /M. Field number (FN) Field number of an objective lens describes the lateral extent of the sample plane transmitted by the objective. It is defined as FN = DFOV M, where DFOV is the diameter of the field of view (FOV) on the sample plane and M is the magnification of the objective. In an objective lens, FOV is determined by a physical or a virtual aperture in the objective’s intermediate image plane. Therefore, the maximum FOV one can capture with a microscope is limited by either the detector size or the objective’s FN, whichever is smaller. Illumination For a non-fluorescent sample, an illumination source is projected onto the sample to carry the sample’s field distribution through the optical system and be captured by the detector. Illumination can be classified into two categories: spatially coherent 6 illumination and spatially incoherent illumination. Spatially coherent illumination can be described as an illumination field originating from a single monochromatic point source. The entire extent of the illumination field incident on the sample can be described by ψcoh (x, y) = a(x, y)e jφ(x,y), (1.1) where (x, y) is the spatial coordinate on the sample plane, j is the complex number, a(x, y) is the amplitude of the field, and φcoh (x, y) is the phase of the field that describes the relative phase delay of the wavefront. Therefore, the fields of the sample at different points can constructively and destructively interfere with each other as they propagate in space. One can manipulate the propagation process by controlling both the amplitude and the phase of the illumination. On the other hand, a spatially incoherent illumination can be considered as an illumination field originating from multiple sources that are not correlated with each other. It can be described by the intensity distribution: ψincoh (x, y) = |a(x, y)| 2 . (1.2) Because it lacks phase, each point on the sample plane propagates independently without interfering with fields originating from other points. For a fluorescent sample, the illumination is provided at the wavelength that can excite the fluorophores in the sample to produce fluorescent signals. Regardless of the coherence of the illumination, the fluorophores’ signals are only dependent on the illumination intensity as the fluorophores are independent from each other. Transfer function The optical signal transmitted by a standard microscope from the sample to the detector plane can be described as a linear process: ψcoh (x, y) = h(x, y) ∗ s(x, y), (1.3) where ψcoh (x, y) is the optical field at the image plane, ∗ is the convolution operator, s(x, y) is the optical field at the sample plane, and h(x, y) is the impulse response of the system, also known as the point-spread function (PSF) [2]. The impulse response is related to the smallest resolvable spot of the microscope. In the spatial frequency domain, the impulse response is described by the transfer function that defines the bandwidth of the system. For a coherent illumination, the transfer function has the shape of the objective’s aperture in its Fourier plane. For a circular aperture, 7 the transfer function is a circular bandpass with a unit amplitude in the 2D spatial frequency space, with the cutoff frequency of NA/λ, where λ is the wavelength of the sample’s optical field. A coherent transfer function is often referred to as the CTF or the pupil function. In the case of an incoherent optical field, the image formed by the microscope is only related to the intensity of the field: ψincoh (x, y) = |h(x, y)| 2 ∗ |s(x, y)| 2 . (1.4) Thus, the PSF of an incoherent imaging system is |h(x, y)| 2 . In the frequency domain, its transfer function is called the optical transfer function (OTF). By the Fourier relationship, OTF can be found by the autocorrelation of CTF: ∫ ∞ OTF(x, y) = (1.5) CTF∗ (x 0, y0)CTF(x + x 0, y + y0)dx 0dy0, −∞ where the superscript ∗ indicates complex conjugation. Due to this property, OTF is twice the width of CTF (i.e., cutoff frequency is 2NA/λ) and has a gradual cut-off in the spatial frequency. Any aberrations due to the imperfections in the objective lens are contained in its pupil function, and the pupil function can also vary across the FOV of the microscope. This will be described further in the next section. 1.3 Limitations in imaging with physical lenses A perfect aberration-free optical lens simply does not exist in reality. As such, all optical imaging systems constructed from a finite number of optical surfaces are going to experience some level of aberration issues. This simple fact underpins the extraordinary amount of optical design efforts that have gone into the design of optical imaging systems. In broad terms, optical imaging system design is largely a complex process by which specialized optical elements and their spatial relationships are chosen in order to minimize aberrations and provide an acceptable image resolution over a desired field of view (FOV) [9]. Its limitations can be summarized as follows: a finite space-bandwidth product, spatially varying aberrations, and chromatic dispersion. Finite space-bandwidth product (SBP) A lens’s diffraction-limited focus spot is governed by the NA. The lens’s aperture diameter can be increased to focus a beam into a tighter spot at the center, but this comes at the expense of increased aberrations in the surrounding area. Multiple optical surfaces can be introduced to mitigate the aberrations. The more optical surfaces available to an optical designer, the greater the extent the aberrations can 8 Objective lens (M/NA/FN) 1.25x/0.04/26.5 2x/0.08/26.5 4x/0.16/26.5 10x/0.3/26.5 20x/0.5/26.5 40x/0.75/26.5 60x/0.9/26.5 100x/1.3/26.5 Resolution at λ = 532 nm 8113 nm 4057 nm 2028 nm 1082 nm 649 nm 433 nm 361 nm 250 nm SBP (megapixels) 21.5 MP 33.5 MP 33.5 MP 18.9 MP 13.1 MP 7.4 MP 4.7 MP 3.5 MP Table 1.1: The resolution and the space-bandwidth product of select Olympus microscope objective lenses. Resolution here is given by Rayleigh criterion, r = 0.61λ NA , where r is the distance from the center of a diffraction limited spot to its first zero for a circular aperture with a given NA. SBP is calculated by dividing the area of the circular FOV (diameter is d = FN/M) by (r/2)2 , which is the pixel area that meets the sampling requirement. The SBP remains constant for low magnification with low NA, but it drastically decreases for high NA objectives due to the difficulty of accounting for increasing aberrations. be minimized. However, this physical system improvement approach for minimizing aberrations has reached a point of diminishing returns in modern optics. Microscope objectives with 15 optical elements have become commercially available in recent years [10], but it is unlikely that another order of magnitude of optical surfaces will be supported within the confines of an objective in the foreseeable future. Moreover, this strategy for minimizing aberration is never expected to accomplish the task of completely zeroing out aberrations. In other words, any optical system’s space-bandwidth product (SBP), which scales as the product of system FOV and inverse resolution, can be expected to remain a design bound dictated by the residual aberrations in the system. As shown in Table 1.1 [14], the decreasing SBPs of a set of Olympus objectives with increasing NA illustrate the greater difficulty in accounting for aberrations for smaller focus spots. FPM overcomes this limit by computationally synthesizing the SBP of multiple images in post processing, as will be described in the following chapter. Spatially varying aberrations Closely related to the finite SBP limit, a lens cannot have a perfectly diffractionlimited spot across its entire FOV. A planar wave incident on a lens’s surface with an angular offset with respect to the lens’s optical axis does not produce a spot with a simple lateral shift on the other side of the lens because refraction is a non-linear 9 process. Refraction is governed by Snell’s law, which relates the light’s incident angle, θi , at the refracting interface and the refractive index of the medium before the interface, ni , to the refracted angle, θr , and the refractive index after the interface, nr , by ni sin θi = nr sin θr . For small angles, we can use the paraxial approximation of Snell’s law, ni θi = nr θr , which linearizes the relationship between the incident and refracted angles such that a plane wave with a tilted illumination angle produces a laterally shifted focal spot on the other side of the lens. However, this approximation immediately breaks down for larger angles and introduces aberrations in the lens’s focal spot [3]. Multiple lens elements can mitigate this effect, but only to a certain extent. The spatially varying aberrations are often represented by the pupil function that depends on different spatial locations. Aberrations can also be caused by the sample. In a typical imaging scenario, the object to be imaged may not have a flat surface and introduce spatially varying defocus aberration in the captured image. In a conventional microscope system, multiple focal planes of the sample need to be captured and patches of different focal depths need to be manually selected to acquire an all-in-focus image. We later show that we can use FPM and modifications of FPM to automatically find the spatially varying pupil function and correct for it in the captured images in post processing. Chromatic dispersion Most optical media have refractive index values that vary with the wavelength. Thus, a polychromatic light beam experiences different refractive angles for its wavelength components upon incidence at an optical interface, resulting in chromatic dispersion. In this thesis, the optical waves of interest have narrow frequency bandwidths (i.e., < 20 nm) such that they are assumed to be monochromatic and free from chromatic dispersion. The rest of the thesis is outlined as follows. In Chapter 2, we introduce the principle of FPM, its advantages, and limitations to lay the foundational knowledge for understanding our endeavor in computational aberration compensation. Chapter 3 illustrates a direct use case of FPM that overcomes the challenges of a standard microscope when capturing a wide FOV image of a sample with uneven surface. In Chapter 4, we demonstrate appending a simple insight to the FPM procedure to correct for aberrations in a fluorescence image with the FPM-reconstructed pupil function. Finally, Chapter 5 describes an FPM-inspired algorithm pipeline that can compensate for a general optical system without spatial coherence requirement of illumination or the sample field. 10 References [1] G. Gasson. The oldest lens in the world: a critical study of the Layard lens. The Ophthalmic Optician, pages 1267–1272, 1972. [2] J. Goodman. Introduction to Fourier Optics. McGraw-Hill, 2008. [3] R. D. Guenther. Modern Optics. John Wiley & Sons, Inc., 1990. [4] D. Hillmann, H. Spahr, C. Hain, H. Sudkamp, G. Franke, C. Pfåffle, C. Winter, and G. Hüttmann. Aberration-free volumetric high-speed imaging of in vivo retina. Scientific Reports, 6:35209, 2016. [5] A. Kumar, D. Fechtig, L. Wurster, L. Ginner, M. Salas, M. Pircher, and R. Leitgeb. Noniterative digital aberration correction for cellular resolution retinal optical coherence tomography in vivo. Optica, 4(8):924–931, 2017. [6] A. Levin, R. Fergus, F. Durand, and W. T. Freeman. Image and depth from a conventional camera with a coded aperture. ACM Transactions on Graphics (TOG), 26(3):70, 2007. [7] M. Levoy, R. Ng, A. Adams, M. Footer, and M. Horowitz. Light Field Microscopy. ACM Transactions on Graphics, 25(3), 2006. [8] W. L. Liang, J. G. Pan, and G. D. J. Su. One-lens camera using a biologically based artificial compound eye with multiple focal lengths. Optica, 6(3):326– 334, 2019. [9] A. W. Lohmann, R. G. Dorsch, D. Mendlovic, Z. Zalevsky, and C. Ferreira. Space-bandwidth product of optical signals and systems. J. Opt. Soc. Am. A, 13(3):470–473, 1996. [10] G. McConnell, J. Trägårdh, R. Amor, J. Dempster, E. Reid, and W. B. Amos. A novel optical microscope for imaging large embryos and tissue volumes with sub-cellular resolution throughout. eLife, 5:e18659, 2016. [11] A. Ozcan and E. McLeod. Lensless imaging and sensing. Annual Review of Biomedical Engineering, 18:77–102, 2016. [12] N. Shemonski, F. South, Y. Z. Liu, S. Adie, S. Carney, and S. Boppart. Computational high-resolution optical imaging of the living human retina. Nature Photonics, 9(7):440–443, 2015. [13] G. Zheng, S. A. Lee, Y. Antebi, M. B. Elowitz, and C. Yang. The ePetri dish, an on-chip cell imaging platform based on subpixel perspective sweeping microscopy (SPSM). Proceedings of the National Academy of Sciences, 108: 16889–16894, 2011. 11 [14] G. Zheng, X. Ou, R. Horstmeyer, J. Chung, and C. Yang. Fourier Ptychographic Microscopy: A Gigapixel Superscope for Biomedicine. Optics and Photonics News, 25:28–33, 2014. [15] C. Zhou and S. Nayar. What are good apertures for defocus deblurring? In Proceedings of IEEE International Conference on Computational Photography, pages 1–8. IEEE, 2009. 12 Chapter 2 FOURIER PTYCHOGRAPHIC MICROSCOPY Fourier ptychographic microscopy (FPM) is a computational coherent imaging method developed by Zheng et al. in 2013 [24]. FPM circumvents the challenges of adding more optical elements for improving an optical system’s performance by recasting the problem of increasing the system’s space-bandwidth product (SBP) as a computational problem that can be solved after image data have been acquired. Rather than striving to get the highest quality images possible through an imaging system, FPM acquires a controlled set of low-SBP images, dynamically determines the system’s aberration characteristics computationally, and reconstitutes a high-SBP, aberration-corrected image from the original controlled image set. FPM shares its roots with ptychography [16], synthetic aperture imaging [20], and structured illumination microscopy (SIM) [3], which numerically expand the SBP of the imaging system in the spatial or spatial frequency domain by capturing multiple images under distinct illumination patterns and computationally synthesizing them into a higher-SBP image. In this chapter, we will provide an overview of the basic principles of FPM’s data acquisition and synthesis process, the advantages offered by FPM over the standard microscope, and its limitations. 2.1 Basic principles of FPM FPM can be divided into the image acquisition process and the computational post-processing. Unlike the standard microscope whose captured image’s quality is dependent on the performance of the optical system at the moment of capture, FPM offloads the stringent requirement on the physical setup to the computational algorithms that can correct for the system’s imperfections with mere numerical operations. FPM’s technology makes use of two key elements: (1) the fact that illuminating a specimen with an oblique plane wave in a microscope setup results in laterally shifting the specimen’s Fourier spectrum in the back focal plane of the objective lens; and (2) the ability to recover the phase information of the specimen from its intensity image by a phase retrieval algorithm [24]. 13 Figure 2.1: Experimental setup of FPM. a) An LED matrix is fitted to a standard microscope system, replacing its condenser. An LED matrix placed sufficiently far away (∼8 cm) from the sample provides angularly varying planar illumination. The objective lens’s backfocal plane (i.e., Fourier plane) is indicated by the yellow box. b) Different illumination angles correspond to sampling different spatial frequency information of the sample. c) LEDs are lit up one by one, and the pitch of the LEDs should be small enough to satisfy the overlap requirement of the spatial spectrum regions in the Fourier space. Image acquisition process In a standard FPM system, images of the sample are collected through a lownumerical-aperture (NA) objective with the target illuminated with a series of angularly varied planar or quasi-planar illumination. The sample is thin such that it can be treated as a 2D object. An example FPM setup, as shown in Fig. 2.1, involves a conventional microscope setup with its light source replaced by an LED matrix. The LEDs are quasi-monochromatic sources with 10 nm bandwidth, which for our purposes can assume to be monochromatic [4]. The light field emitted from one LED can be approximated as a plane wave within a small region at the sample plane because the LED is small (200 µm by 200 µm) and the large distance (∼8 cm) between the LED and the specimen increases the spatial coherence of the LED. The plane wave has vector components (k x, k y ) associated with the oblique illumination of the LED on the sample. Lighting up one LED causes the sample’s Fourier spectrum to be shifted by (k x, k y ) at the objective lens’s back-focal plane due 14 to the Fourier relationship, and the finite numerical aperture (NA) of the objective lens acts as a low-pass filter that transmits only a small subregion at the center of the shifted Fourier spectrum. The low-passed Fourier component is further propagated to the image plane and is captured by a camera. LEDs are lit up sequentially so that the captured images contain partially overlapping Fourier subregions (at least 30% overlap [17] that together span the entire Fourier spectrum of the specimen. The captured images are broken up into smaller tiles of spatially coherent regions, whose dimension is given by the van Cittert-Zernike theorem: L = 0.61λz/a [8], where λ is the LED’s center wavelength, z is the LED-to-sample distance, and a is the radius of the LED’s active area. For example, with the LED dimension of 200 µm x 200 µm, distance between the LED matrix and the sample of 8 cm, and the wavelength of 630 nm, the coherence length is ∼307 µm. Also, within this dimension, the spatially varying aberrations are assumed to be constant. Each tile is considered independently in the later reconstruction process. Phase retrieval algorithm As with synthetic aperture synthesizing, we then stitch the data from the collected series in the spatial frequency domain. Unlike synthetic aperture imaging, we do not have direct knowledge of the phase relationships between each image data set. In FPM, we employ phase retrieval by taking advantage of the partial information overlap amongst the image set to converge on the correct phase relationships during the stitching process [24]. For each tile in the FOV, the phase retrieval algorithm iteratively solves for the expanded Fourier spectrum by using the intensity images as a constraint in the spatial domain and the objective’s finite NA as a shifting lowpass constraint in the Fourier domain. Various phase retrieval algorithms for FPM exist [21], but here we describe embedded pupil function recovery (EPRY) based on the Gerchberg-Saxton algorithm that can simultaneously solve for the objective’s pupil function in the phase retrieval process by interleaving an update step in the algorithm [11, 19]. EPRY’s procedure is described in Fig. 2.2. The images captured with the low NA are labeled as Im (x) (m = 1, 2, 3, ...), the simulated low NA complex field as ψ(x), the synthesized sample spectrum as S(k), and the reconstructed pupil function as P(k) where x = (x, y) and k = (k x, k y ). δ is an arbitrary small value for numerical stability in the division process. In this process, the algorithm is able to not only generate the high-resolution complex amplitude image of the specimen, but also separate the pupil function and the specimen’s complex amplitude function from the image. The pupil function contains 15 Figure 2.2: The embedded pupil function recovery (EPRY) algorithm. EPRY is a phase retrieval process that can synthesize the FPM’s intensity dataset into a complex, high SBP image and simultaneously reconstruct the pupil function of the imaging system. the aberrations of the microscope associated with its lens system and the defocus caused by misalignment of the sample at the focal plane, both of which can be dependent on the tile’s location in the FOV. Therefore, the specimen’s complex amplitude function is an aberration-corrected, diffraction-limited complex image of the specimen. By stitching the tiles together in the spatial domain, we obtain the specimen’s intensity image and phase information separated from aberrations across the entire FOV. 2.2 Advantages Computational postprocessing of FPM presents several advantages over a conventional standard microscope. Increased space-bandwidth product Given the FOV of the conventional objective lens, the resolution of the FPM’s image is no longer limited by the objective’s NA but is given by the summation of the objective’s NA and the illumination NA. For example, an FPM system with the FOV of a 2x objective lens (13 mm diameter) and the resolution typically offered 16 by a 20x objective (0.4 NA) has been achieved [24]. This is extremely valuable for imaging microscope slides (e.g., histology samples) which typically need to be mechanically raster scanned to generate a wide-FOV, high-resolution image for accurate diagnosis. As an example application, we demonstrated that accurate white blood cell counting can be performed on a blood smear slide image captured with FPM [1]. Aberration characterization and correction Automatic aberration compensation is perhaps the quintessential component of FPM for achieving a higher throughput with a given lens. It overcomes the fundamental tradeoff between the achievable FOV and resolution inherent in the physical lens design. Not only can FPM account for the imperfections in the optical system, it can also compensate for the aberrations induced by the sample being imaged by the system. For example, a sample may have an uneven surface profile that places different locations of the sample at different focal planes. We explore this further in Chapter 3 and show that it can reconstruct a high resolution image that is in-focus across the FOV. This particular sub-discipline of FPM has matured to the level that it is even possible to use a very crude lens to obtain high quality images that are typically associated with sophisticated imaging systems [6] – this drives home the fact that correcting aberration computationally is a viable alternative to physical correction. We describe our investigation of extending this aberration compensation ability to imaging applications and to different microscope modalities in the later chapters of this thesis. Phase imaging An image reconstructed by FPM’s phase retrieval is a complex field that contains both amplitude and phase information of the optical field transmitted through the sample. The phase describes the relative optical path length at different locations that can be related to the different refractive indices or thicknesses. This is a very valuable piece of information in imaging biological samples as most cells are transparent. Development of a phase contrast microscope by Frits Zernike [22, 23] and a differential interference contrast microscope by Georges Nomarski [10] revolutionized biological sciences by making cells visible without any staining process. Recently, a digital staining method based on a deep neural network has been demonstrated to successfully convert a label-free sample’s phase information into a stained image from which histologists can more easily diagnose [15]. This 17 eliminates the time- and resource-consuming process of staining the sample prior to imaging. Thus, the high resolution phase imaging capability of FPM can help increase the throughput of medical diagnosis and shed light on valuable structural information of cells and other transparent objects. 2.3 Limitations There are several limitations that currently prevent FPM from replacing the standard microscope. Thin sample limit FPM as implemented by Zheng et al. works for a thin sample with the thickness h within q |k z | < σk 0 , where q sample limit given by  thin  σ = π/(k 0 h) and |k z | = 2 2 2 2 2 2 max k 0 − k0 − (k 0 NAobj ) , k 0 − k 0 − (k0 NAillum ) [12]. The thin sample limit roughly corresponds to the synthesized 3D PSF’s axial resolution, where the 3D PSF can be reconstructed in the 3D spatial frequency domain by synthesizing Ewald caps corresponding to different illumination angles [14]. For a thicker sample, its Fourier spectrum can no longer be assumed to be two-dimensional, and a tilted planar illumination on the sample will not simply shift the sample’s Fourier spectrum in the objective lens’s backfocal plane. Attempting FPM’s reconstruction process on such a sample will lead to severe artifacts that will render unusable images. There have been several developments to address this issue, such as modeling a thick sample as a series of multiple 2D slices [13, 18], or assuming a weakly scattering sample [5]. However, they are not robust as their assumptions are often not applicable in real life experiments. Whereas a standard microscope still allows one to roughly section through a thick sample by changing the focal plane, FPM would simply fail to recover a viable image when those assumptions are not met. Spatial coherence requirement FPM is a coherent imaging modality. The illumination on the sample needs to be spatially coherent for the FPM image acquisition process to capture the correct spatial frequency information of the sample. However, there are many imaging scenarios where providing such illumination is difficult or impossible. These scenarios include: 1) optical systems where the illumination on a sample is provided via a medium with unknown index variations; 2) optical systems where space is so confined that it is not feasible to employ optical propagation to create quasi-planar optical fields as done with LEDs placed at a distance; and 3) optical systems where 18 the optical field at the sample plane is spatially incoherent by nature (e.g., fluorescence emission). This makes FPM quite sensitive to the illumination arrangement, unlike a standard microscope, which only needs a simple condenser with a thermal light source, such that slightly mispositioned LEDs or a small change in the source to sample distance can prevent FPM from rendering a correct reconstructed image. In Chapter 5, we present a method developed off of FPM that overcomes this strict spatial coherence requirement to allow for robust computational aberration correction of a general optical system. Speed As FPM requires multiple images to generate a single image, the acquisition and processing steps present a bottleneck in the imaging speed. Nevertheless, there have been various efforts to address this issue such as by multiplexing the illumination pattern to reduce the number of acquisition [19], replacing the LEDs with a highpower laser source [2, 7], and incorporating a deep neural network to estimate FPM’s reconstructed images from fewer frames [9]. With the rapid increase in computation power and improvement in camera sensors, the issue of speed appears to be solvable in the near future. References [1] J. Chung, X. Ou, R. P. Kulkarni, and C. Yang. Counting White Blood Cells from a Blood Smear Using Fourier Ptychographic Microscopy. PLoS ONE, 10:e0133489, 2015. [2] J. Chung, H. Lu, X. Ou, H. Zhou, and C. Yang. Wide-field Fourier ptychographic microscopy using laser illumination source. Biomedical Optics Express, 7(11):4787–4802, 2016. [3] M. G. L. Gustafsson. Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy. Journal of Microscopy, 198(2): 82–87, 2000. [4] R. Horstmeyer and C. Yang. A phase space model of Fourier ptychographic microscopy. Optics Express, 22(1):338–358, 2014. [5] R. Horstmeyer, J. Chung, X. Ou, G. Zheng, and C. Yang. Diffraction tomography with Fourier ptychography. Optica, 3(8):827–835, 2016. [6] T. Kamal, L. Yang, and W. M. Lee. In situ retrieval and correction of aberrations in moldless lenses using Fourier ptychography. Optics Express, 26(3):2708– 2719, 2018. 19 [7] C. Kuang, Y. Ma, R. Zhou, J. Lee, G. Barbastathis, R. R. Dasari, Z. Yaqoob, and P. T. C. So. Digital micromirror device-based laser-illumination Fourier ptychographic microscopy. Optics Express, 23(21):26999–27010, 2015. [8] L. Mandel and E. Wolf. Optical Coherence and Quantum Optics. Cambridge University Press, 2008. [9] T. Nguyen, Y. Xue, Y. Li, L. Tian, and G. Nehmetallah. Deep learning approach for Fourier ptychography microscopy. Optics Express, 26(20):26470–26484, 2018. [10] G. Nomarski. Interferential polarizing device for study of phase object. U.S. patent 2,924,142 (February 9, 1960). [11] X. Ou, G. Zheng, and C. Yang. Embedded pupil function recovery for Fourier ptychographic microscopy. Optics Express, 22(5):4960–4972, 2014. [12] X. Ou, R. Horstmeyer, G. Zheng, and C. Yang. High numerical aperture Fourier ptychography: principle, implementation and characterization. Optics Express, 23(3):3472–3491, 2015. [13] X. Ou, J. Chung, R. Horstmeyer, and C. Yang. Aperture scanning fourier ptychographic microscopy. Biomedical Optics Express, 7(8):3140–3150, 2016. [14] C. Park, S. Shin, and Y. Park. Generalized quantification of three-dimensional resolution in optical diffraction tomography using the projection of maximal spatial bandwidths. Journal of the Optical Society of America A, 35(11): 1891–1898, 2018. [15] Y. Rivenson, T. Liu, Z. Wei, Y. Zhang, K. de Haan, and A. Ozcan. PhaseStain: the digital staining of label-free quantitative phase microscopy images using deep learning. Light: Science & Applications, 8:23, 2019. [16] J. M. Rodenburg and R. H. T. Bates. The theory of super-resolution electron microscopy via Wigner-distribution deconvolution. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 339 (1655):521–553, 1992. [17] J. Sun, Q. Chen, Y. Zhang, and C. Zuo. Sampling criteria for Fourier ptychographic microscopy in object space and frequency space. Optics Express, 24 (14):15765–15781, 2016. [18] 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. [19] L. Tian, X. Li, K. Ramchandran, and L. Waller. Multiplexed coded illumination for Fourier Ptychography with an LED array microscope. Biomedical Optics Express, 5(7):2376–2389, 2014. 20 [20] T. M. Turpin, L. H. Gesell, J. Lapides, and C. H. Price. Theory of the synthetic aperture microscope. Proc. SPIE, pages 230–240, 1995. [21] L. H. Yeh, J. Dong, J. Zhong, L. Tian, and M. Chen. Experimental robustness of Fourier ptychography phase retrieval algorithms. Optics Express, 23(26): 33214–33240, 2015. [22] F. Zernike. Phase contrast, a new method for the microscopic observation of transparent objects. Physica, 9(7):686–698, 1942. [23] F. Zernike. How I Discovered Phase Contrast. Science, 121(3141):345–349, 1955. [24] G. Zheng, R. Horstmeyer, and C. Yang. Wide-field, high-resolution Fourier ptychographic microscopy. Nature Photonics, 7(9):739–745, 2013. 21 Chapter 3 HIGH-RESOLUTION, WIDE-FIELD IMAGING OF SAMPLES WITH UNEVEN SURFACES High SBP imaging is highly desired in various biomedical applications such as hematology and digital pathology [3]. Typically, it is performed by mechanically scanning the microscope slide with a high-NA objective and digitally stitching the captured images in post-processing to create a wide-field and high-resolution image. There can be an issue in this pipeline if the sample on the microscope slide is not flat so that keeping the objective lens at the same focal plane would lead to some scanned regions being out of focus. A work-around could be to have a microscopist specify preset focal planes for different sample regions, which makes the imaging process extremely laborious, or to capture multiple focal planes at every sample location and pick the images at the right focus, which is highly data-inefficient. Some microscopes incorporate an optical coherence tomography (OCT) system to find the correct depth of the sample prior to capturing an image. In this chapter, we illustrate how FPM can elegantly address the uneven sample issue thanks to its high SBP imaging capability and its ability to correct for spatially varying aberrations. The following is about our work in collaboration with Williams et al. [5] on imaging circulating tumor cells (CTC) captured on a thin filtration membrane. 3.1 Introduction Circulating tumor cells (CTCs) in peripheral blood have emerged in recent years as a valuable biomarker with strong potential to improve prognosis and diagnosis of cancer. Assaying for CTCs requires only a simple, minimally invasive blood draw, providing a unique opportunity for repeated sampling in patients to monitor both metastatic disease as well as therapeutic response in real time. Thus, the enumeration of CTCs with respect to progression-free survival, overall survival, and therapeutic response has been widely reported on in a number of solid tumor malignancies. Among various ways to isolate CTCs, size-based isolation of CTCs from whole blood has been attempted since 1960s [2], and has been revisited more recently. Utilizing the well-known characteristic that the malignant cells are larger than surrounding normal blood cells, CTCs are isolated by using microfilters fabricated with a defined pore size, which allow for the passage of smaller blood 22 cells to pass while capturing larger CTCs [1, 6]. Where the sensitivity and efficiency of affinity-based CTC enrichment strategies rely primarily on tissue- and/or tumorspecific cell surface bio- markers with the potential for highly variable inter tumor expression, size-based enrichment technologies are “antigen expression-agnostic,” allowing analysis of CTCs in tumor types with low or no target antigen expression. With the potential to overcome the limitations accompanying other platforms, sizebased CTC enrichment strategies possess technical limitations of their own. The most significant of these limitations is that following sample processing, the surface of the microfilters becomes uneven with short, microscale modulations. Because captured cells are randomly dispersed throughout the microfilter, often times multiple CTCs are present on different focal planes. This technical limitation requires that the user must constantly change focus while viewing and imaging cells of interest on the microfilter, making sample analysis labor intensive, time consuming, and inefficient. Automated imaging systems have been developed and are widely available on a number of microscopy platforms, and could potentially alleviate these complications. However, such systems cannot readily be employed for CTC analysis using filtration-based technologies due to their inability to focus multiple areas with different focal planes within the same frame. Further complicating this issue, when viewing microfilters under a microscope to identify CTCs, the entire filtration area must be viewed by systematically moving up and down the microfilter in columns (i.e., in y-axis plane), or side to side across the microfilter in rows (i.e., x-axis plane) manually using the stage manipulator knobs. In instances where the observer does not appropriately align adjoining columns and or rows on the microfilter, small areas where tumor cells may reside can escape the field of view (FOV) and fail to be counted, or single events can be counted more than once if adjoining columns or rows are aligned slightly overlapping each other. Without the ability to produce images for CTC analysis in an automated fashion, the potential for inter-operator variability and inconsistent analyses between users and collaborating institutions reviewing the same samples is dramatically increased. These technical limitations, taken together, prevent widespread analysis of CTCs using filtration-based CTC enrichment technologies. To address these challenges, we present the adaptation of Fourier ptychographic microscopy (FPM) for the identification and enumeration of CTCs captured by sizebased enrichment. The use of FPM allows for the rapid generation of continuous, high-resolution images over large areas of interest. Importantly, we also describe 23 the ability to perform digital refocusing of images generated by FPM on a frameby-frame basis, allowing us to create focused images of frames containing cells of interest in multiple focal planes, thus traversing the limitations in automated imaging of uneven surfaces produced by other commercially available technologies. Here, we present an assessment of our ability to analyze CTCs captured on a previously described membrane microfilter device developed by Williams et al. using FPM, evaluating the efficiency of CTC detection as well as the image quality of CTCs generated by FPM relative to standard microscopic analysis. 3.2 Experimental setup Membrane microfilter device for circulating tumor cell capture and characterization Williams et al. have developed a membrane microfilter device for the size-based isolation of CTCs in blood [6]. Microfilters for CTC capture and analysis are fabricated using a precisely defined, stepwise photolithography process on paryleneC. The actual manufacturing process is beyond the scope of this thesis. Parylene-C provides excellent mechanical, optical, and biocompatible properties to be used as the membrane filter. The microfilters have been demonstrated to achieve 90% recovery rate of tumor cell isolation [6]. CTC sample preparation SKBR-3 breast cancer cells were seeded into 5 ml of whole blood from a normal, healthy donor at various concentrations, processed by the microfilter device, and labeled with CK and CK45 by double marker IHC on glass microscope slides. A total of 11 replicates of this experiment were performed. Nine normal donor blood samples were seeded with a range of tumor cells (15 to 800), representing the lower and upper limits of CTCs typically identified from clinical blood samples using the microfilter device. As a negative control, two normal donor blood samples containing no tumor cells were processed by the microfilter device. FPM setup The FPM setup is the same as the one in Chapter 1 and it consists of the following components: a Olympus BX 41 microscope with 0.08 NA using a Plan APO 2x objective lens (Olympus, Center Valley, Pennsylvania), a KAI-29050 interline CCD camera with 5.5 µm pixel size (Eastman Kodak, Rochester, New York) attached to a computer for image capturing and processing, and a square light-emitting 24 Figure 3.1: (a) The FPM setup consists of (from the bottom) an LED matrix for sample illumination, a microscope system with a 2x objective, and a camera connected to a computer. (b) The Fourier spectrum point of view. The center red subregion corresponds to the spatial frequency of the low-resolution image captured with plane waves with k x = k y = 0. The off-center red subregion correlates to an oblique angle illumination with wavevector (k x, k y ). (c) Light from an LED at an oblique angle corresponds to a plane wave with a k vector (k x, k y ). diode (LED) array for illumination. The LED matrix contains 32 x 32 surfacemounted, full-color LEDs and adjacent LEDs are laterally separated by 4 mm. The full-color LED has central wavelengths of 632 nm (red), 532 nm (green), and 472 nm (blue), each offering a spatially coherent quasi-monochromatic source with ∼20 nm bandwidth. When an LED on the matrix is activated, its light field incident on the sample can be approximated as a plane wave due to the large distance (∼8 cm) between the LED and the sample plane. The angular illumination can be characterized by its in-plane wavevector (k x, k y ) within the coordinate system, as depicted in Fig. 3.1. Illuminating a sample with a plane wave of a wavevector in the space domain is equivalent to shifting the center of the sample’s frequency spectrum in the Fourier domain. Because the objective acts as a circular, low-pass filter in Fourier domain, each captured image carries information describing a shifted small subregion, which is geometrically defined by the pupil function of the microscope, of the sample’s frequency spectrum. Images are captured by the camera with single LEDs of one color activated in sequence acting as the light source, and a phase 25 retrieval algorithm is used to stitch the subregions together in the Fourier domain to form a high-resolution complex image, which contains both amplitude and phase information. Red, green, and blue images are acquired separately by altering the color of the LED for each set of acquisitions. Single color channels are used individually for the entire image capturing and reconstruction process. Thus, the procedure is conducted a total of three times for red, green, and blue channels. These images are later combined to create a full-color image, such as the ones shown in Fig. 3.2. In our hands, each sample takes ∼3 min to capture and 10 min to reconstruct each of the three color channels, for a total of 39 min to produce a color image of the entire field. For sample analysis by the observer, the complex whole FOV image created by FPM is refocused and partitioned into a total of ∼300 tiles to be used for CTC identification on the entire field. 3.3 Analysis of CTC images by FPM and comparison with standard microscopy Remarks on FPM-recovered CTC images Because of the uneven surface of the microfilter containing cells of interest, an image taken by a standard microscope suffers from defocus. As shown in Fig. 3.2, when some parts in the FOV are in focus [Figs. 3.2-(d1-1)], other parts can be blurry [Fig. 3.2-(d1- 2)]. FPM with the embedded pupil function recovery (EPRY) algorithm, as described in Chapter 2, can automatically correct for this aberration [4]. When EPRY stitches the subregions in the Fourier domain, it does so by iteratively recovering both the Fourier spectrum of the sample and the pupil function. Because of the pupil function, which contains the aberration of the lens system and the defocus caused by surface unevenness, the Fourier spectrum of the sample is separated during the EPRY process. Performing an inverse Fourier transformation on the sample’s Fourier spectrum results in an aberration-free, flattened image of the microfilter. As shown in Fig. 3.2-(d1), the defocus is corrected automatically and the components of the image shown in Figs. 3.2-(d1-1) and Figs. 3.2-(d1-2) are well focused. Because the algorithm also recovers the pupil function, the sample’s depth information can be obtained and graphed, as shown in Fig. 3.3(a). The graph shows the microfilter’s severe surface unevenness, which requires an imaging device with different focusing levels across the sample. FPM, with its refocusing capacity of up to 300 µm, is able to image the entire microfilter in focus. Fig. 3.3(b1) and Fig. 3.3(c1) show two different areas at different depth levels on the microfilter, with Fig. 3.3(b1) being in focus and Fig. 3.3(c1) being out of focus. After applying FPMs refocusing 26 Figure 3.2: (a) Full field-of-view color image of the entire microfilter containing captured tumor cells by FPM. Magnified FPM images (b-d1) selected from different areas of the microfilter show detailed morphology of tumor cells, where all sections are well in focus because of the automatic EPRY-FPM program. (d2) A standard microscope (with 40x objective) image shows the corresponding region to (d1), but because of the uneven surface of the microfilter, subregions (d2-1) and (d2-2) cannot be focused simultaneously. Also, its field of view is limited when compared to (d2), as seen by the aperture’s outline at the edge of the image, in contrast to FPM’s wide field-of-view that can provide high-resolution images of the entire micro- filtration area. 27 Figure 3.3: (a) A microfilter’s surface profile characterized by the pupil function recovered by EPRY. The focal plane of the objective is at 150 µm. The maximum difference in-depth across the filter is about 250 µm, which is within FPM’s refocusing capacity of 300 µm. In this case, the captured microfilter image is sectioned into 17 x 17 tiles, and EPRY iteratively characterizes each tile’s defocus level in its high-resolution image reconstruction process. (b1)-(c2) Small subregions are extracted from two different surface levels, showing before and after refocusing by EPRY. (b1) is already very close to the focal plane, so there is only a minor improvement after refocusing, as in (b2). (c1) is not in the focal plane and is blurry. (c2) shows the refocused result. algorithm, both areas are brought sharply in focus, as shown in Figs. 3.3(b2) and Fig. 3.3(c2). Comparison with standard microscopy In all samples tested, captured tumor cells were enumerated first by “standard microscopy” – Axio Imager M1 with an Apochromat 20x/0.8 NA objective lens (Carl Zeiss Microscopy LLC, Thornwood, New York) and the tumor cell counts were compared to enumeration done by FPM in corresponding samples. A Pearson’s correlation between cells identified by both technologies on corresponding microfilter samples was conducted to evaluate the consistency of tumor cell identification by FPM relative to standard microscopy. As demonstrated in Fig. 3.4, 28 Figure 3.4: Graph demonstrating the correlation between tumor cell count by standard microscopy (y-axis) and tumor cell count by FPM (x-axis) in corresponding microfilter samples, where each data point represents a single trial with tumor cells enumerated by both methods. the R2 for tumor cells identified in corresponding microfilter samples was 0.99932, indicating a strong correlation between tumor cells counted by both technologies. This demonstrates FPM as a suitable imaging method that provides comparable accuracy in tumor cell detection and enumeration to standard microscopy. References [1] V. De Giorgi, P. Pinzani, F. Salvianti, J. Panelos, M. Paglierani, A. Janowska, M. Grazzini, J. Wechsler, C. Orlando, M. Santucci, T. Lotti, M. Pazzagli, and D. Massi. Application of a Filtration- and Isolation-by-Size Technique for the Detection of Circulating Tumor Cells in Cutaneous Melanoma. Journal of Investigative Dermatology, 130(10):2440–2447, 2010. [2] R. L. Fleischer, H. W. Alter, S. C. Furman, P. B. Price, and R. M. Walker. Particle Track Etching. Science, 178(4058):255–263, 1972. [3] F. Ghaznavi, A. Evans, A. Madabhushi, and M. Feldman. Digital Imaging in Pathology: Whole-Slide Imaging and Beyond. Annual Review of Pathology: Mechanisms of Disease, 8:331–359, 2013. [4] X. Ou, G. Zheng, and C. Yang. Embedded pupil function recovery for Fourier ptychographic microscopy. Optics Express, 22(5):4960–4972, 2014. [5] A. Williams, J. Chung, X. Ou, G. Zheng, S. Rawal, Z. Ao, R. Datar, C. Yang, and R. Cote. Fourier ptychographic microscopy for filtration-based circulating tumor cell enumeration and analysis. J. Biomed. Opt., 19(6):066007, 2014. 29 [6] S. Zheng, H. Lin, J. Q. Liu, M. Balic, R. Datar, R. J. Cote, and Y. C. Tai. Membrane microfilter device for selective capture, electrolysis and genomic analysis of human circulating tumor cells. Journal of Chromatography A, 1162 (2):154–161, 2007. 30 Chapter 4 CORRECTING FOR ABERRATIONS IN FLUORESCENCE IMAGE WITH FPM-RECOVERED PUPIL FUNCTION Due to the limitations of physical lens design, imaging systems tend to exhibit aberrations that vary across their FOV [24]. For wide FOV microscopes (i.e., such as those used for gigapixel imaging via FPM [5, 14, 46, 50]), aberrations significantly deteriorate image quality near the FOV boundary. Thus, aberration correction becomes a critical step to consistently resolve sample features across the entire image plane [54]. As demonstrated in prior work [34] and Chapter 3, the coherent nature of FPM allows us to algorithmically account for the microscope’s complex and spatially varying aberrations using a procedure termed embedded pupil function recovery (EPRY). The output of the EPRY algorithm is both a highresolution estimate of a sample’s amplitude and phase (i.e., the same as FPM), along with an estimate of the microscope’s pupil function. The same algorithm cannot be directly applied to improve fluorescence imaging because both FPM and EPRY only operate in coherent imaging schemes. Fluorescence is incoherent, so it does not change in response to angularly varying illumination. However, some insights from FPM and EPRY do in fact carry over to improve fluorescence resolution using different illumination schemes [6]. Here, we make use of the insight that the aberrations impacting the fluorescence image can be corrected by using the pupil function determined by EPRY. 4.1 Introduction Bright-field imaging is one of the most popular microscope modalities. A brightfield microscope typically illuminates a specimen with white light and captures an image of the transmitted light. Bright-field image contrast offers information about sample structure [30]. In combination with phase contrast methods, one may quantitatively measure sample absorption, thickness, and dispersion [25]. Another increasingly important microscope modality in biology is fluorescence imaging. Fluorescence can help visualize chemical compositions and structures at a molecular level through appropriate labeling with fluorophores [23, 43]. By illuminating a tagged sample at the fluorophores’ excitation light and imaging at the fluorophores’ emission wavelength, biologists can easily identify labeled regions exhibiting chem- 31 ical properties of interest. Combining a bright-field and fluorescence image allows one to locate fluorescing regions relative to the specimen’s underlying structure. In this chapter, we report a method that can simultaneously generate a highresolution, coherent bright-field image, along with an aberration-corrected fluorescence image, across a wide imaging system field-of-view (FOV). Specifically, we first improve bright-field image resolution using Fourier ptychographic microscopy (FPM), a recently developed computational technique that processes a sequence of images acquired under angularly varying illumination [53]. As a phase retrieval method, FPM also recovers the phase information of the sample [33], and can additionally estimate the imaging system’s aberrations [34]. Second, we then capture a fluorescence image using the same microscope setup, and correct its aberrations using the aberration map measured via FPM. Specifically, FPM outputs a complex aberration map in the form of a spatially varying pupil function, with which we compute the microscope’s spatially varying incoherent point spread function (PSF) to use in fluorescence image deconvolution for aberration removal. The aberration removal procedure involves deconvolution with the associated incoherent PSF using Tikhonov regularization. Since both the coherent bright-field and fluorescence images are acquired from the same microscope in quick succession and without any movement, we demonstrate that this direct connection between the pupil function and incoherent PSF is accurate. We note that our deconvolution method is aimed at correcting the aberrations in the fluorescence image caused by imperfections in the imaging system. If done perfectly, it should allow us to render a fluorescence image at the diffraction limit of the system. This deconvolution is distinct from super-resolution deconvolution methods that can reconstruct information outside the passband by incorporating a priori knowledge of sample structure (e.g., one consisting only of lines or points [17, 19, 32, 36]) or an infinite SNR [4, 41]. Here, we are primarily concerned with imaging unknown biological samples, and thus make no a priori assumptions about their structures. The combination of aberration-corrected fluorescence image and a high-resolution bright-field image in a large FOV format makes this method potentially suitable for applications, such as counting fluorescently tagged bacteria [28, 37], studying cell migration dynamics [29], and tracing cell lineage [40]. The remainder of this chapter is outlined in the following manner: in Section 4.2, we describe our algorithm for aberration characterization and aberration removal in incoherent imaging. In Section 4.3, we demonstrate the implementation of our 32 algorithm in a microscope while imaging a sample of fluorescent beads, as well as a fixed, stained, and fluorescently tagged HeLa cell slide. We show that our approach corrects for spatially varying aberrations to enable accurate, wide-field imaging. In Section 4.4, we quantify how much fluorescence signal must be collected for a successful image deconvolution. Finally, we conclude by discussing the implications of our results. 4.2 Principle of operation As noted above, we characterize the aberrations within our FPM setup using the EPRY phase retrieval algorithm [34] as described in Chapter 1. Here, we offer a short summary of it. First, FPM begins by capturing a series of low-resolution images under varying plane wave illumination. In this FPM setup, we use an LED array placed at a large distance behind the sample of interest for our angularly varying illumination source. For N different LEDs, we turn on one LED at a time and capture a unique image. For the remainder of this manuscript, we treat the light emitted from each LED as quasi-monochromatic and spatially coherent across the sample area. A detailed discussion of the requirements and impact of source coherence for FPM is in [20]. Assuming the imaged sample is thin as described in Chapter 1, an angularly varying illumination source creates a laterally shifting sample spectrum (i.e., the Fourier transform of the complex sample) in the back-focal plane of the objective lens. At this plane, the finite extent of the objective lens aperture stop (connected to the lens NA) acts as a low-pass filter. In a conventional infinity-corrected microscope objective lens, the extent of the aperture stop defines its cutoff spatial frequency, which in turn dictates its smallest resolvable feature at the image plane. By laterally shifting the sample spectrum via varied LED illumination and acquiring a sequence of limited resolution images, FPM acquires a data set that contains sufficient information to reconstruct an image with a wider spectrum than offered by a single image. However, since the sample spectrum is complex, and the image sensor can only record the intensity of light, extending the sample spectrum beyond the original objective lens bandpass using the N acquired images is not direct. To solve this inverse problem, FPM applies a phase retrieval reconstruction algorithm [7, 8]. Several algorithms are available [10, 21, 51] including standard non-linear solvers based upon alternating projections, which we utilize in this work, as well as more advanced techniques [1, 42]. With careful initialization, the alternating projection algorithm arrives at the solution with high success rate despite not having convergence guarantees [33, 34, 46, 47, 50, 53]. All 33 techniques require a certain amount of redundancy within the acquired image set. Here, we vary the LED angles such that consecutive images, which correspond to uniquely windowed areas of the sample spectrum, overlap by approximately 65% in the Fourier domain. Independent of the algorithm, we will refer to the application of phase retrieval to expand the aperture stop bandpass as Fourier ptychographic (FP) reconstruction. In the end, we effectively increase the NA of our system by the illumination NA defined via the maximum LED illumination angle [35]. Simultaneous to complex sample reconstruction, it is also possible to use the same set of FPM images to iteratively determine the microscope’s pupil function. We refer to this joint solver as EPRY [34]. The pupil function is a complex function that summarizes a microscope’s aberrations in addition to its spatial frequency bandpass as a variation in amplitude and phase at its aperture plane. In an aberrated microscope, each resulting image is connected to the product of the sample spectrum (shifted laterally for an oblique LED illumination) with the same pupil function, which is nonzero only within the bandpass defined by the microscope’s aperture. Using joint optimization, it is thus possible to simultaneously determine the pupil function within FP reconstruction by splitting the iterative FPM update process into two steps. First, EPRY may use the image data and current pupil function estimate to update the sample spectrum estimate at all even iterations. Then, EPRY may apply the image data and the current spectrum estimate to update the pupil function estimate at all odd iterations. A similar joint process is also used to simultaneously determine the amplitude and phase of the shifting probe field in Xray ptychography [13, 31, 44, 45], leading to a significant boost in solver accuracy. To account for the spatially varying aberrations across the microscope’s FOV, the captured images are segmented into small tiles, for which the aberrations can be considered as spatially invariant within each tile [53, 54]. Here, we select a tile area that is smaller than the entire sample, but larger than approximately 20 x 20 sensor pixels projected onto the image plane, which helps mitigate numerical artifacts during reconstruction. The EPRY process is simply applied separately to distinct image sub-tiles. Image tiling separately considers different areas of sample and image plane, thus allowing EPRY to determine a unique aberration function from each sub-region of the microscope’s FOV. The tile-specific pupil functions obtained from the EPRY algorithm provide an accurate physical model of the microscope’s coherent transfer function, as detailed in the appendix of this chapter, which means that one can also deduce from the pupil functions a correct aberration model for an incoherent imaging scheme. Given 34 that both imaging modalities originate from within the same microscope at similar wavelengths and assuming no inter-image movement, the fluorescence image PSF is simply given by the squared magnitude of the Fourier transform of the pupil function [12]:  2 hm (x, y) = F −1 Pm ( fx, fy ) , (4.1) where (x, y) represent the spatial coordinates of the image plane, ( fx, fy ) represent the coordinates in the aperture plane (i.e., are the Fourier conjugates of the spatial coordinates), hm (x, y) is the incoherent PSF, Pm ( fx, fy ) is the CTF, and F −1 represents an inverse Fourier transform operation. Our goal is to use this information to remove the aberrations from the fluorescence image. We segment the fluorescence image into the same sub-tiles as in the FPM reconstruction process to ensure the aberration map acquired from EPRY corresponds to the same sub-regions in the fluorescence image. We represent one tile of the fluorescent sample as om (x, y), where m stands for the mth tile. When om (x, y) is imaged by an optical system, it is degraded by the possibly spatially variant incoherent PSF, hm (x, y), before reaching the detector. The detected image intensity, i m (x, y), is further corrupted by noise, nm (x, y), originating from the sample’s background signal, shot noise, and detector noise [11]. The imaging process can be represented as i m (x, y) = hm (x, y) ∗ om (x, y) + nm (x, y). (4.2) The goal of the aberration removal in fluorescence imaging is to recover the object om (x, y) from the corrupted image signal, i m (x, y). In the Fourier domain, Eq. 4.2 is represented as follows: Im ( fx, fy ) = Hm ( fx, fy ) · Om ( fx, fy ) + Nm ( fx, fy ), (4.3) where Im ( fx, fy ), Hm ( fx, fy ), Om ( fx, fy ), and Nm ( fx, fy ) are the Fourier transforms of i m (x, y), hm (x, y), om (x, y), and nm (x, y), respectively. Inverting Eq. 4.3 to solve for Om ( fx, fy ) is an ill-posed problem due to the nature of Hm ( fx, fy ), also known as the optical transfer function (OTF). Unlike the coherent transfer function, the OTF can have numerous zeros within its bandpass and its values are very low near the bandpass’s edges [12], which means that the sample information is lost or overwhelmed by noise at these spatial frequencies. Various inversion methods have been developed to account for this information loss, which typically rely upon regularization parameters [39]. Wiener deconvolution [49] is one method 35 that attempts to determine an estimate, Õm ( fx, fy ), for the original object signal, Om ( fx, fy ), by the following: Õm ( fx, fy ) = G m ( fx, fy ) · Im ( fx, fy ), (4.4) where G m ( fx, fy ) is an inverse filter defined as G m ( fx, fy ) = Hm∗ ( fx, fy ) 2 | Nm ( fx, fy ) | Hm ( fx, fy ) + 2 | Om ( fx, fy ) | 2 . (4.5) While it is reasonable to assume Nm ( fx, fy ) is flat white Gaussian noise, Om ( fx, fy ) is hard to determine unless we have some prior knowledge of the sample’s spatial 2 2 distribution. For simplicity, we set Nm ( fx, fy ) / Om ( fx, fy ) to a constant K [11, 38], essentially converting G m ( fx, fy ) to a Tikhonov regularization algorithm [2]. K acts like a regularizer: a smaller K produces sharper details while amplifying noise in the captured image whereas a larger K makes the algorithm more robust to noise at the expense of details. This value is determined visually so that the deconvolution minimizes the background noise while recovering the most detail in the final image [38]. In the end, our final estimate of the original fluorescent object is given by " # ∗(f , f ) H m x y · Im ( fx, fy ) . (4.6) õm (x, y) = F −1 2 Hm ( fx, fy ) + K Equation 4.6 outputs the resulting image after Tikhonov regularization with our assumed noise and sample distribution model. Gaussian noise is a reasonable assumption for the noise in our system because we capture images with long exposures for high signal-to-noise ratio. However, for photon-limited settings where Poisson noise becomes severe and dominant, maximum likelihood deconvolution methods that assume a Poisson imaging process would be more suitable, such as RichardsonLucy deconvolution [22]. We apply Eq. 4.6 for all the sub-tiles to acquire a full FOV aberration-corrected fluorescence image. 4.3 Experimental demonstration Fluorescent microspheres sample To experimentally demonstrate both FPM and fluorescence image enhancement, we use the modified 4f setup diagrammed in Fig. 4.1. Here, the 4f arrangement consists of an f = 200 mm tube lens (ITL200, Thorlabs) and an f = 50 mm Nikon lens (f/1.8D AF Nikkor). We also place an iris at the back-focal plane of the 36 Figure 4.1: The experimental setup of joint FPM-fluorescence imaging. The 4f setup has a camera lens, an iris at the back focal plane, and a tube lens. The filter is included only for fluorescence imaging and FPM imaging of the same color channel. First, raw FPM images are captured using varied-angle illumination provided by the LED matrix. Then, a fluorescence image is captured with the illumination from the excitation LED. f = 50 mm lens, as shown in Fig. 4.1, to allow for a user-controllable NA. The 4f setup has an NA of NAcollection = 0.085 with an M = 3.87 magnification. The sample is placed at the front focal plane of the f = 50 mm lens, and a CCD detector (pixel size 5.5 µm, Prosilica GX6600) is used to capture each image. It is able to capture an area of 6.2 mm by 9.3 mm on the sample plane, defined by the detector’s sensor size in this setup. The setup includes two illumination sources, one for each imaging modality. For high-resolution FPM imaging, we use an LED array placed ∼80 mm behind the sample with 32 x 32 individually addressable elements (pitch size 4 mm), of which we only use a 15x15 segment (225 images captured in total). The increase in NA provided by the LEDs is NAillum = 0.33. The overall system NA is NAsys = NAcollection + NAillum = 0.33 + 0.085 = 0.415, which is a factor of 4.9 resolution gain over the 4f setup with a single plane wave illumination. For fluorescence imaging, we use a separate excitation LED along with an appropriate emission filter that we insert in the optical path behind the iris. FPM raw images and the fluorescence image are captured separately. The emission filter is present only for capturing the fluorescence image and the FPM raw images of the same color channel, which ensures they both image the same spectral range. As a first experimental demonstration of our algorithm pipeline, we image a sample slide containing both 10 µm green-fluorescing microspheres and 15 µm nonfluorescing microspheres (Fisher Scientific). First, we acquire 225 full FOV lowresolution images under variable LED illumination (<3 seconds average exposure 37 time). Only green LEDs are used for this demonstration. We focus our attention on a single image tile in Fig. 4.2, which corresponds to a 360 µm x 360 µm region located 3.8 mm away from the center of the image FOV. From our set of 225 low-resolution images captured under variable green LED illumination (in Fig. 4.2(a)), we apply the EPRY algorithm for 50 iterations (5 seconds per tile) to reconstruct the intensity and phase of the sample, along with the pupil function of the optical system (in Fig. 4.2(b)). Note, this pupil function is only valid for this particular image tile, and the reconstructed complex image is much sharper than the raw images due to resolution improvement and aberration removal. Next, we capture a green fluorescence image by illuminating the sample with our blue excitation LED (470 nm, Thorlabs), in Fig. 4.2(c) (2 minutes exposure time). Both FPM images and the fluorescence image are captured using a green band-pass filter (530 nm, 43 nm pass band) behind the iris. The fluorescence image is severely blurred because it is captured near the edge of the imaging system’s FOV. We can derive the incoherent PSF that characterizes the fluorescence image blur, hm (x, y), from the EPRY algorithm’s computed pupil function, Pm ( fx, fy ), following Eq. 4.1. We then deconvolve the PSF-induced blur from our fluorescence image in Fig. 4.2(c) using Eq. 4.6, which creates the sharp image in Fig. 4.2(d) (< 1 second per tile). Deconvolution removes many of the negative effects of lens-induced aberration. Four fluorescent beads that are challenging to resolve in the raw image are clearly distinguishable after deconvolution. In Fig. 4.2(e), we plot the one dimensional profile through two neighboring beads, both with raw and deconvolved data. This plot highlights that not only does our deconvolution improve image contrast, but it also locates the centroid of each bead with greater accuracy. Specifically, the ratio of the lower peak to the dip between each peak is 0.946 for the raw image and 0.526 for the deconvolved result. The raw image shows the beads’ separation distance is 8.5 µm, which is 15% below the bead’s diameter of 10 µm. The deconvolved image shows the separation of 9.9 µm, which is within the 5% tolerance value indicated within the beads’ manufacturing specifications. We also note that the lateral shift caused by the system’s aberrations is corrected for in the deconvolution result, as indicated by the shift of the inter-peak dip. This is important for identifying the correct spatial correspondence between the high-resolution FPM image and the fluorescent sample image. 38 Figure 4.2: Experimental reconstruction of an improved-resolution amplitude and phase image, along with an aberration-corrected fluorescence image, using FPM. (a) Series of low-resolution coherent green images captured with variable-angle illumination from an LED matrix. (b) EPRY is applied on the low-resolution images to generate a high-resolution, complex image of the sample, along with the characterization of the microscope’s pupil function. (c) The pupil function is converted into an incoherent PSF, which is the blur kernel induced to a fluorescence image by the imaging system in the same sample region. (d) The PSF is deconvolved from the fluorescence image using Eq. 4.6 to generate an aberration-corrected image. (e) The one-dimensional profile of two beads, in the raw and deconvolved images, shows improved contrast and position accuracy. (f) In the overlay of FPM and fluorescence images, the fluorescence signal is localized with good centroid accuracy after deconvolution. 39 Figure 4.3: Comparison between (a) a raw image of a specific sample ROI shifted 4.63 mm off the optical axis (left), its deconvolution result (right), and (b) a raw image captured with the same sample ROI centered on the optical axis, i m (x, y). The PSF recovered from FPM, hm (x, y), is deconvolved from the image to remove aberrations in (a). The result shows close resemblance to the image in (b), i center (x, y), which we assume is minimally impacted by system aberrations. Fluorescently tagged HeLa cells To demonstrate the benefit of combined bright-field and fluorescence imaging with a biological sample, we use our system to image a sample of HeLa cells stained for both fluorescence and bright-field. The sample is a microscope slide with 90% confluent HeLa cells that are fixed, stained with Wright’s stain from Fisher Scientific’s PROTOCOL Hema 3, and fluorescently tagged with DAPI from Life Technologies. The sample is first imaged for FPM reconstruction by sequentially illuminating it with the red, green, and blue channels of the 225 LEDs in our array. We insert a bandpass filter (460 nm, 80 nm band) to capture blue illumination images from within the same spectral range as DAPI’s fluorescence emission. Second, we capture a blue fluorescence image by illuminating the sample with a UV LED (365 nm, Thorlabs) at DAPI’s excitation wavelength. The sample is exposed for 21 minutes, which is the optimal time for the given sample and imaging setup as will be detailed in Section 4.4. We then reconstruct a high-resolution FPM image, along with a spatially varying pupil function, by repeating the EPRY algorithm on different image sub-tiles. Finally, we apply Eq. 4.6 to deconvolve each tile of the fluorescence image. As a qualitative performance test of our deconvolution, we capture two fluorescence images of the same sample region of interest (ROI), shifted to two different locations within the optical system’s FOV. We capture one image with the ROI located at the center of the image FOV (i.e., along the optical axis, shown in Fig. 4.3(b)). In this region, we assume the effects of imaging system aberrations are minimal and treat this image as the ground truth for comparison. We could deconvolve this image with the associated PSF, which could potentially improve the image 40 contrast. However, the result of Tikhonov regularization can at best only be an estimate of the true sample function, as given by the assumption on K in Eq. 4.6. Therefore, we do not perform deconvolution on the center image. We then capture a second image of the same ROI after translating the sample laterally 4.63 mm away from the center, as shown in Fig. 4.3(a). It is clear that lens aberrations more significantly impact image quality off the optical axis, as indicated by the blurry outlines of the DAPI-stained nuclei. The same sub-figure also shows the results of our deconvolution, using the incoherent PSF found via FPM image capture and post-processing of the same image tile. Deconvolution both increases the nuclei’s contrast and highlights features that are not otherwise visible in the raw fluorescence image to the left, such as the shape of the nuclei and their separation gaps. We can confirm the accuracy of our deconvolution result in Fig. 4.3(a) by comparing it with the minimally aberrated image of the same ROI in Fig. 4.3(b). We see close agreement between all deconvolved image features, with few apparent artifacts. Some “hot” pixels, visible as white or black dots within the raw, blurry image in Fig. 4.3(a), lead to ringing artifacts around each after deconvolution. However, their presence does not drastically degrade the quality of the deconvolution image. Fig. 4.4 offers an example of combining both our resolution-improved bright-field imaging and fluorescence imaging capabilities within a single biological application. Here, similar to Fig. 4.3, we process low-resolution images obtained under variableLED illumination using the EPRY algorithm. Prior to reconstruction, the lowresolution color images have severe aberrations because they are captured near the edge of the FOV. We can reconstruct high-resolution, aberration-corrected full-color FPM images to elucidate sharp image features, such as the nucleoli present in the HeLa cells’ nuclei, cell morphology via the reconstructed phase, and boundaries between different cells. Fig. 4.4(d) is the full FOV fluorescence image of the HeLa sample, and Figs. 4.4(a-c) correspond to zoomed-in regions of the image to demonstrate that we are able to obtain high resolution FPM images and aberrationcorrected fluorescence images throughout the entire image FOV. Generating a full FOV image requires applying EPRY and deconvolution algorithms to individual subtiles separately (∼2200 tiles, 5 seconds per tile). We can also combine the fluorescence image with FPM data to generate a phase gradient image overlaid with fluorescence emission. The phase gradient + fluorescence images in Figs. 4.4(a-c) show that, with the help of the cell’s structural information provided through the phase gradient, the fluorescent regions of HeLa cells are indeed the nuclei. It is also possible to distinguish cells in telophase (Fig. 4.4(c) arrow) from the ones in 41 Figure 4.4: Demonstration of FPM and aberration-corrected fluorescence imaging across different regions of a large image FOV. (a)-(c) correspond to regions labeled in the full FOV fluorescence image in (d). With FPM, we improve bright-field image resolution (1st and 2nd column) and characterize the spatially varying pupil functions (6th column). Hot pixels and chromatic aberrations appearing in the low-resolution color images are suppressed and corrected after FPM reconstruction. We correct aberrations in the fluorescence image using each pupil function (3rd and 4th column). The phase gradient images from FPM can be combined with fluorescence data to elucidate the structures and the locations of the nuclei of HeLa cells (5th column). The cell membrane morphology elucidated by the phase gradient can differentiate between cells undergoing cytokinesis (arrow in (b)) and cells in telophase (arrow in (c)), which is otherwise difficult to do in the fluorescence images alone. 42 the cytokinesis stage of cell mitosis (Fig. 4.4(b) arrow), by examining the cell membrane morphology within the phase gradient image. This type of information is missing from fluorescence data alone. It is important to note that, because the aberrations in FPM and fluorescence imaging are corrected by functions derived from the same pupil function, we also observe accurate spatial correspondence (i.e., alignment) between each imaging modality. 4.4 Quantifying the adequate amount of fluorescence signal for deconvolution We observe that the quality of our deconvolution results is proportional to the signalto-noise ratio (SNR) of the input images. As noted previously, Eq. 4.6 is based on Gaussian noise model; thus, applying it on an image dominated by Poisson noise will not result in successful deconvolution. It is important to investigate at what SNR level our method can start yielding accurate deconvolved images. Although we can image the DAPI fluorescence dye (used to stain our HeLa cells in Fig. 4.3 and Fig. 4.4) for an extended period of time without photobleaching to obtain a high SNR, many fluorescence dyes are not as robust [30, 43]. For a weakly fluorescing sample (e.g., low fluorophore concentration), one needs to extend the detector exposure time and runs into the risk of photobleaching the sample before the end of an image capturing process. Thus, the system’s parameters need to be adjusted appropriately to image the sample before photobleaching occurs. Also, to our knowledge, this is the first published work to demonstrate deconvolution of fluorescence images with the aberration function derived from an iterative computational method using coherent microscopy. All aforementioned reasons prompted us to conduct further analysis on our deconvolution method. In the following paragraphs, we systematically analyze how an imaging system’s parameters play a role in the fluorescence image capturing process and quantify how much fluorescence signal is required for our system to offer good deconvolution results. In a CCD detector, photons hitting the sensor generate electrons in each pixel. In the readout process, the number of electrons are converted into a discrete photon count value, S, by a gain factor, g. S includes various kinds of noise, σs , such as dark current, read noise, quantization noise, and shot noise [18]. The per-pixel SNR is then given by SNR = S . σ (4.7) We assume shot noise, which originates from the quantum nature of light, to be the dominant noise and ignore others for simplicity of our analysis. Further assuming p our CCD’s gain to be g = 1, S is directly related to σs by: σs = hSi [18], where 43 hSi is the average pixel count. Thus, SNR increases with an increase in signal value p (i.e., SNR = hSi ). To determine a threshold SNR of fluorescence signal needed for successful deconvolution, we capture multiple fluorescence images of the region in Fig. 4.3 (i.e., 4.63 mm away from center), each with a different exposure time. An increase in exposure time leads to an increase in S registered by our detector’s pixels for the fluorescent nuclei of HeLa cells. Since each HeLa nucleus has a different fluorescence signal level due to different amount of bound DAPI molecules, we select one representative HeLa nucleus within the ROI for quantifying signals in this study. Specifically, we use a small uniform region within the nucleus, as indicated by the yellow box in Fig. 4.5(b), to calculate average SNR values. As we vary the exposure time, t, from 45 seconds to 37.5 minutes, SNR varies from 3.46 to 24.16. We then deconvolve each variable exposure image i m (u; t) via Eq. 4.6 with the same deconvolution filter to create a set of sharpened reconstructions, õm (u; t). We compare each of these sharpened images to a reference image, i center (u), created again by centering the same image ROI along the optical axis. We assume this centered image is a minimally aberrated, and thus a useful comparison benchmark. It has a much greater SNR (SNR = 41.14) compared to the off-center images. Note that õm (u; t) is the estimate of the sample function, om (u), and i center (u) is the convolution of om (u) with the incoherent PSF of the center FOV (which is only minimally aberrated), hcenter (u), as given by Eq. 4.2 (assuming ncenter (u) ∼ 0 for a high SNR). We do not apply Tikhonov regularization to i center (u) to obtain õcenter (u) because õcenter (u) would at best only be an estimation for the true sample function om (u), especially without the knowledge of the sample’s spatial frequency distribution, Om ( fx, fy ), and the exact noise distribution,Nm ( fx, fy ), as indicated in Eq. 4.5. Thus, we treat i center (u) as the ground truth in our analysis. In order to quantify the accuracy of õm (u; t), we convolve õm (u; t) with hcenter (u) to generate i˜m (u; t) for appropriate quantitative comparison between the estimated and the ground truth sample function. We then compare each i˜m (u; t) with i center (u) by computing their normalized mean square error (NMSE) [9], which is a function of exposure time, t: NMSE(t) = 2 Σu i center (u) − αt i˜m (u; t) Σu |i center (u)| 2 . (4.8) Here, u = (x, y), and αt is given by αt = Σui center (u)i˜∗m (u; t) 2 Σu i˜m (u; t) , (4.9) 44 Figure 4.5: (a) Normalized mean square error (NMSE) of the raw, i m (u; t), and deconvolved image, i˜m (u; t), in Fig. 4.3(a) is plotted against the detector exposure time. Higher exposure is equivalent to higher SNR of the captured image. NMSE starts to plateau for both raw and deconvolved images after about 21 min of exposure (raw image SNR = 18.1) as indicated by the green broken line, with the deconvolved image’s NMSE = 0.0057 being 31% lower than the raw image’s NMSE = 0.0083. (b) Example images of the ROI used for this study. The reference image is generated by centering the ROI and capturing with our imaging system. We use the area inside the small yellow box to quantify the SNR. Raw images are captured with varying exposure time (45 sec and 21 min shown here) while the ROI is 4.63 mm away from the center. We use the same small area in each raw image to quantify its SNR. Deconvolved images are generated by applying Eq. 4.6 to the raw images. Longer exposure provides better SNR images and deconvolution results. which allows for NMSE metric to be invariant to a constant multiplicative factor between images being compared. For a quantitative validation of the improvement in image quality after our deconvolution algorithm, we also calculate the NMSE between the raw away-from-center images, i m (u; t), and the reference image, i center (u). As shown in Fig. 4.5, the deconvolved images are able to consistently provide lower NMSE values than do the raw images, with NMSE reduced by 31% in the plateau region (NMSE = 0.0057 versus 0.0087). Both plots start to plateau at about t = 21 min, at which the SNR in the raw captured image is SNR = 18.1. We set this as the adequate SNR value of a fluorescence image for successful deconvolution. We then quantify how much fluorescence a sample needs to emit for our particular imaging system to obtain the adequate SNR image. This quantity can guide the sample preparation process and the imaging system setup so that the sample’s 45 bound fluorophore concentration level and the imaging system’s parameters are made suitable for successful fluorescence image deconvolution. Assuming that the DAPI fluorophores in the nucleus of interest emit n photons per area isotropically (i.e., spherical wavefront from each point on the nucleus), the proportion entering our system, nm, with the given numerical aperture, NA, can be determined by calculating how much fraction of the spherical area is captured by the system. NA limits the light acceptance angle and results in capturing a section of the sphere, a spherical cap [48]. The ratio of the spherical cap to the total spherical area is equal p to m = (1 − 1 − (NA/ni )2 )/2, where ni is the index of refraction of the medium between the sample and the imaging lens. For our system, ni = 1, so it is omitted in subsequent calculations. The photons reach our CCD detector and are converted into electrons with the conversion ratio determined by the quantum efficiency, q (q = 46% at λ = 470 nm for our system). With the magnification of our imaging system, M = 3.87, CCD gain assumed to be g = 1, and the pixel area, p = 30.25 µm2 , the output signal for a pixel is given by p nmpq npq(1 − 1 − NA2 ) S= = . (4.10) M 2M We look at a uniform region within a HeLa nucleus, so S = hSi. Thus, our average SNR is given by s p √ pq(1 − 1 − NA2 ) hSi = n SNR = p . (4.11) 2M hSi Rearranging for n, we find: n = SNR2 2M . p pq(1 − 1 − NA2 ) (4.12) Substituting the values for our system parameters, the number of photons per area required for an adequate SNR in our system is n = 50400 photons/µm2 . This condition allows for our deconvolution algorithm to achieve an accurate result with a low NMSE value. The fluorescent sample needs to be exposed for an adequate amount of time to satisfy the n requirement. For a sample that cannot provide this much signal either due to photobleaching or restrictions on exposure time, the system parameters need to be adjusted appropriately according to Eq. 4.4 to achieve a similar SNR level. 46 4.5 Discussion In summary, we have developed an FPM system capable of generating a wide-field, aberration-corrected fluorescence image with correct spatial correspondence to the high-resolution FPM image. It utilizes the spatially varying pupil function from EPRY algorithm to correct for the aberrations in the fluorescence image. The resulting color intensity image, phase image, and the fluorescence image can be combined to provide multiple layers of information about the sample’s morphology, chemical properties, and functions. The wide FOV fluorescence and high-resolution brightfield image can be beneficial to various biological imaging studies including, but not limited to, cell lineage tracing, counting bacteria, and cell migration. For imaging scenarios requiring the use of a higher NA objective lens, the same algorithm may be used with the FPM setup adjusted accordingly by simply varying the distance of the LED matrix to the sample to allow for steeper angles of illumination, as demonstrated in [35]. The advantages of this aberration-removal method with the same pupil function from FPM over other deconvolution methods are the following: 1) deconvolution by the PSF obtained from FPM allows for proper overlap between the fluorescence image and the FPM image. This is important when overlaying the two images, especially when the pupil function has large asymmetric aberrations in different regions of the FOV. Asymmetric aberrations cause a noticeable lateral shift in the images. Since the presented deconvolution method uses the incoherent PSF determined directly from the pupil function obtained from EPRY, both the fluorescence image and FPM’s high-resolution coherent image can be corrected for aberrations with the same amount of lateral shift. Other deconvolution methods that do not use the pupil function characterized by FPM, such as blind deconvolution, would not be able to generate a corrected image with the proper lateral shift because the methods do not account for the absolute location of the image [26], meaning that the deconvolved image can have any arbitrary lateral shift as its solution. Without the same lateral shift for FPM and fluorescence images, the direct overlay of the images will not have the correct spatial correspondence between them. 2) Because the system can characterize the PSF of the fluorescence image, one can refer to the vast amount of literatures on non-blind deconvolution methods to deconvolve the aberrations from the image. Various high-fidelity and noise-robust methods, such as Ref. [27, 52], can be implemented once the PSF is known. The disadvantage of this method is that it requires a whole set of coherent images 47 to be acquired with an array of coherent light sources of the same wavelength as the fluorescent signals in order to correct for the aberration in the fluorescence image. Not only is this data intensive, but it also requires additional set of light sources for each fluorescence color channel. In Chapter 5, we will discuss another variation of FPM that can address these issues. 4.6 Appendix It is necessary to verify that EPRY is able to recover a physically accurate model of the microscope’s pupil function. Although Ref. [34] shows that EPRY algorithm converges to the actual pupil function in computer simulations, it only demonstrates the algorithm’s robustness in experiments via visual inspection of reconstructed sample images. Here, we proceed to confirm that it converges to the correct pupil function in real imaging situations by comparing the generated pupil function with that obtained from another pupil function recovery method. We crosscheck EPRY with an alternative method capable of characterizing spatially varying aberrations of an imaging system as described in Ref. [54]. Its principle is based on a phaseretrieval technique on a sample via phase diversity introduced by imaging the sample at various defocus planes [3, 15, 16]. The basic procedure is as follows. The sample has identical targets, such as microspheres, distributed uniformly across the imaging FOV. A target at the center of the FOV is considered to be minimally aberrated by the imaging system. Its complex function obtained via the phase retrieval technique is therefore treated as the ground truth of the target’s complex field. The aberration function of any region of the imaging system’s FOV can be recovered using the defocus image stack of a target in that certain area and the target’s ground truth. The aberration function with a given number of parameters can be optimized for those parameters with a pattern search algorithm by minimizing the difference between the measured defocus image stack and a simulated defocus stack. The latter is generated by applying the aberration function on the ground truth target and propagating at the same defocus planes as the measured stack. The aberration function that minimizes the difference provides an accurate description of the pupil function in that region of the imaging system’s FOV. For an experimental implementation of this method, we use the same modified 4f setup diagrammed in Fig. 4.1 without the excitation LED. We replace the LED array with a collimated green LED beam to illuminate the sample at a perpendicular angle for this procedure, as shown in Fig. 4.6(a). 48 Figure 4.6: The variations in experimental setup for comparing the two different pupil function recover methods. (a) Collimated light source illuminates the sample perpendicularly, and the sample is brought to different defocus planes in the defocus diversity-based pupil recovery method. (b) LEDs of varied illumination angles act as the light source in FPM setup. EPRY algorithm jointly solves for the sample spectrum and the pupil function using sample images captured under different LED illumination angles. We use a microscope slide sample consisting of uniformly distributed 3 µm and 15 µm diameter microspheres (Fisher Scientific) immersed in oil and sealed with a coverslip as our target. The sample is brought to 17 different defocus planes, ∆z away from the focal plane of the microscope, where ∆z ranges from -400 µm to +400 µm with 50 µm steps. We select a 3 µm diameter bead at the center of the microscope’s FOV as the ground truth and apply the defocus diversity phase retrieval process to obtain its complex function. We choose one 3 µm bead from different regions in the microscope’s FOV, and apply the aforementioned pattern search optimization algorithm to determine the pupil functions in those regions. We optimize the pupil functions for 8 parameters corresponding to the first, second, third, and fourth order Zernike modes, which are namely x-tilt, y-tilt, x-astigmatism, y-astigmatism, defocus, x-coma, y-coma, and spherical aberration. The experimental setup for EPRY is as follows. The setup is identical as above except for a 32 x 32 individually addressable RGB LED array (4 mm pitch) replacing the collimated light source, as shown in Fig. 4.6(b). We place the LED array ∼80 mm 49 Figure 4.7: Full FOV image of a sample of microspheres with associated pupil functions in 3 different sub-regions. The red box indicates the pupil function recovered by EPRY while the blue box stands for defocus diversity-based pupil function recovery method. The two pupil estimation methods correlate well throughout different regions of the FOV. below the sample, and use a 15 x 15 segment in the acquisition process (225 images captured in total). Only green LEDs were used for this procedure. We segment the full FOV into sub-tiles, and reconstruct both the pupil function and the highresolution sample function specific to each tile. In agreement with our assumption, the central FOV shows minimal aberrations. For accurate comparisons between the two aberration characterization methods, we normalize all the pupil functions obtained by EPRY by the center tile’s pupil function. We compare these to the pupil functions in corresponding regions obtained from the defocus diversity method, of which 3 pairs are plotted in Fig. 4.7. We only include the second, third, and fourth order Zernike modes in this comparison. The first order Zernike modes, namely the x-tilt and y-tilt, merely produce lateral shifts in a captured image, so they are not considered as aberrations. We observe that the recovered pupil functions from the two methods correlate well across the FOV. Slight discrepancies can be attributed to imperfections in the collimation of LED and misalignment between the sample’s defocus axis and the microscope’s optical axis in Fig. 4.6(a)’s setup. We conclude that the pupil function generated by EPRY is an accurate physical description of the imaging system, which 50 can be subsequently appropriated for fluorescence imaging. References [1] H. H. Bauschke, P. L. Combettes, and D. R. Luke. Phase retrieval, error reduction algorithm, and Fienup variants: a view from convex optimization. Journal of the Optical Society of America A, 19(7):1334–1345, 2002. [2] M. Bertero and P. Boccacci. Introduction to Inverse Problems in Imaging. Taylor & Francis, 1998. [3] G. R. Brady and J. R. Fienup. Nonlinear optimization algorithm for retrieving the full complex pupil function. Optics Express, 14(2):474–486, 2006. [4] A. J. den Dekker and A. van den Bos. Resolution: a survey. Journal of the Optical Society of America A, 14(3):547–557, 1997. [5] S. Dong, K. Guo, P. Nanda, R. Shiradkar, and G. Zheng. FPscope: a fieldportable high-resolution microscope using a cellphone lens. Optics Express, 5 (10):3305–3310, 2014. [6] S. Dong, P. Nanda, R. Shiradkar, K. Guo, and G. Zheng. High-resolution fluorescence imaging via pattern-illuminated fourier ptychography. Optics Express, 22(17):20856–20870, 2014. [7] V. Elser. Phase retrieval by iterated projections. Journal of the Optical Society of America A, 20(1):40–55, 2003. [8] J. R. Fienup. Phase retrieval algorithms: a comparison. Applied Optics, 21 (15):2758–2769, 1982. [9] J. R. Fienup. Invariant error metrics for image reconstruction. Applied Optics, 36(32):8352–8357, 1997. [10] P. Godard, M. Allain, V. Chamard, and J. Rodenburg. Noise models for low counting rate coherent diffraction imaging. Optics Express, 20(23):25914– 25934, 2012. [11] R. C. Gonzalez and R. E. Woods. Digital Image Processing. Prentice Hall, 2002. [12] J. Goodman. Introduction to Fourier Optics. McGraw-Hill, 2008. [13] M. Guizar-Sicairos and J. R. Fienup. Phase retrieval with transverse translation diversity: a nonlinear optimization approach. Optics Express, 16(10):7264– 7278, 2008. [14] K. Guo, S. Dong, P. Nanda, and G. Zheng. Optimization of sampling pattern and the design of fourier ptychographic illuminator. Optics Express, 23(5): 6171–6180, 2015. 51 [15] B. M. Hanser, M. G. Gustafsson, D. A. Agard, and J. W. Sedat. Phase retrieval for high-numerical-aperture optical systems. Optics Letters, 28(10):801–803, 2003. [16] B. M. Hanser, M. G. Gustafsson, D. A. Agard, and J. W. Sedat. Phase-retrieved pupil functions in wide-field fluorescence microscopy. Journal of Microscopy, 216(1):32–48, 2004. [17] J. L. Harris. Resolving power and decision theory. Journal of the Optical Society of America, 54(5):606–609, 1964. [18] G. E. Healey and R. Kondepudy. Radiometric CCD camera calibration and noise estimation. IEEE Trans. Pattern Anal. Mach. Intell., 16(3):267–276, 1994. [19] R. Heintzmann. Estimating missing information by maximum likelihood deconvolution. Micron, 38(2):136–144, 2007. [20] R. Horstmeyer and C. Yang. A phase space model of Fourier ptychographic microscopy. Optics Express, 22(1):338–358, 2014. [21] R. Horstmeyer, R. Y. Chen, X. Ou, B. Ames, J. A. Tropp, and C. Yang. Solving ptychography with a convex relaxation. New Journal of Physics, 17(5):053044, 2015. [22] G. M. P. Kempen and V. L. J. Vliet. A quantitative comparison of image restoration methods for confocal microscopy. Journal of Microscopy, 185(3): 354–365, 1997. [23] A. R. Kherlopian, T. Song, Q. Duan, M. A. Neimark, M. J. Po, J. K. Gohagan, and A. F. Laine. A review of imaging techniques for systems biology. BMC Syst. Biol., 2(1):74, 2008. [24] R. Kingslake and R. B. Johnson. Lens Design Fundamentals. SPIE Press, 2010. [25] E. Kissa. Dispersions: Characterization, Testing, and Measurement. Marcel Dekker, Inc, 1999. [26] D. Kundur and D. Hatzinakos. Blind image deconvolution. IEEE Signal Processing Magazine, 13(2):43–64, 1996. [27] J. H. Lee and Y. S. Ho. High-quality non-blind image deconvolution with adaptive regularization. Journal of Visual Communication and Image Representation, 22(7):653–663, 2011. [28] L. G. Leff and A. A. Leff. Use of green fluorescent protein to monitor survival of genetically engineered bacteria in aquatic environments. Applied and Environmental Microbiology, 62(9):3486–3488, 1996. 52 [29] C. C. Liang, A. Y. Park, and J. L. Guan. In vitro scratch assay: a convenient and inexpensive method for analysis of cell migration in vitro. Nature Protocol, 2 (2):329–333, 2007. [30] J. W. Lichtman and J. A. Conchello. Fluorescence microscopy. Nature Methods, 2(12):910–919, 2005. [31] A. M. Maiden and J. M. Rodenburg. An improved ptychographical phase retrieval algorithm for diffractive imaging. Ultramicroscopy, 109(10):1256– 1262, 2009. [32] D. A. Nahrstedt and L. C. Schooley. Alternative approach in decision theory as applied to the resolution of two point images. Journal of the Optical Society of America, 69(6):910–912, 1979. [33] X. Ou, R. Horstmeyer, C. Yang, and G. Zheng. Quantitative phase imaging via Fourier ptychographic microscopy. Optics Letters, 38(22):4845–4848, 2013. [34] X. Ou, G. Zheng, and C. Yang. Embedded pupil function recovery for Fourier ptychographic microscopy. Optics Express, 22(5):4960–4972, 2014. [35] X. Ou, R. Horstmeyer, G. Zheng, and C. Yang. High numerical aperture Fourier ptychography: principle, implementation and characterization. Optics Express, 23(3):3472–3491, 2015. [36] C. Pask. Simple optical theory of super-resolution. Journal of the Optical Society of America, 66(1):68–70, 1976. [37] K. G. Porter and Y. S. Feig. The use of DAPI for identifying and counting aquatic microfloral. Limnology and Oceanography, 25(5):943–948, 1980. [38] J. C. Russ. The Image Processing Handbook. CRC Press, fifth edition, 2007. [39] P. Sarder and A. Nehorai. Deconvolution methods for 3-D fluorescence microscopy images. IEEE Signal Processing Magazine, 23(3):32–45, 2006. [40] T. Schroeder. Long-term single-cell imaging of mammalian stem cells. Nature Methods, 8(4S):S30–S35, 2011. [41] P. J. Sementilli, B. R. Hunt, and M. S. Nadar. Analysis of the limit to superresolution in incoherent imaging. Journal of the Optical Society of America A, 10(11):2265–2276, 1993. [42] Y. Shechtman, Y. C. Eldar, A. Szameit, and M. Segev. Sparsity based subwavelength imaging with partially incoherent light via quadratic compressed sensing. Optics Express, 19(16):14807–14822, 2011. [43] D. J. Stephens and V. J. Allan. Light microscopy techniques for live cell imaging. Science, 300(5616):82–86, 2003. 53 [44] P. Thibault, M. Dierolf, A. Menzel, O. Bunk, C. David, and F. Pfeiffer. Highresolution scanning x-ray diffraction microscopy. Science, 321(5887):379– 382, 2008. [45] P. Thibault, M. Dierolf, O. Bunk, A. Menzel, and F. Pfeiffer. Probe retrieval in ptychographic coherent diffractive imaging. Ultramicroscopy, 109(4):338– 343, 2009. [46] L. Tian, X. Li, K. Ramchandran, and L. Waller. Multiplexed coded illumination for Fourier Ptychography with an LED array microscope. Biomedical Optics Express, 5(7):2376–2389, 2014. [47] L. Tian, Z. Liu, L. H. Yeh, M. Chen, J. Zhong, and L. Waller. Computational illumination for high-speed in vitro Fourier ptychographic microscopy. Optica, 2(10):904–911, 2015. [48] J. Waters and T. Wittmann. Quantitative Imaging in Cell Biology: Methods in Cell Biology. Academic Press, 2014. [49] N. Wiener. The Extrapolation, Interpolation and Smoothing of Stationary Time Series. The MIT Press, 1949. [50] A. Williams, J. Chung, X. Ou, G. Zheng, S. Rawal, Z. Ao, R. Datar, C. Yang, and R. Cote. Fourier ptychographic microscopy for filtration-based circulating tumor cell enumeration and analysis. J. Biomed. Opt., 19(6):066007, 2014. [51] C. Yang, J. Qian, A. Schirotzek, F. Maia, and S. Marchesini. Iterative Algorithms for Ptychographic Phase Retrieval. http://arxiv.org/abs/1105.5628. [52] L. Yuan, J. Sun, L. Quan, and H. Y. Shum. Progressive inter-scale and intrascale non-blind image deconvolution. ACM Transactions on Graphics, 27(3): 1–10, 2008. [53] G. Zheng, R. Horstmeyer, and C. Yang. Wide-field, high-resolution Fourier ptychographic microscopy. Nature Photonics, 7(9):739–745, 2013. [54] G. Zheng, X. Ou, R. Horstmeyer, and C. Yang. Characterization of spatially varying aberrations for wide field-of-view microscopy. Optics Express, 21 (13):15131–15143, 2013. 54 Chapter 5 CORRECTING FOR ABERRATIONS IN A GENERAL OPTICAL SYSTEM Up to this chapter, FPM’s ability to determine and correct aberration has been limited to optical setups with a well-defined, spatially coherent field on the sample plane [4–6, 21, 36, 42, 44–46, 48, 52]. However, the removal of spatial coherence constraints is vitally important in allowing us to apply computational aberration correction to a broader number of imaging scenarios. These scenarios include: 1) optical systems where the illumination on a sample is provided via a medium with unknown index variations; 2) optical systems where space is so confined that it is not feasible to employ optical propagation to create quasi-planar optical fields; 3) optical systems where the optical field at the sample plane is spatially incoherent by nature (e.g., fluorescence emission). In this chapter, we show how we developed a computational imaging method capable of reconstructing an optical system’s pupil function by adapting the FPM’s alternating projections to an incoherent imaging modality. The reconstructed pupil function can then be deconvolved from the aberrated image to recover the latent high-resolution image. 5.1 Introduction A perfect aberration-free optical lens simply does not exist in reality. As such, all optical imaging systems constructed from a finite number of optical surfaces are going to experience some level of aberration issues. This simple fact underpins the extraordinary amount of optical design efforts that have gone into the design of optical imaging systems. Reiterating from Chapter 1, optical imaging system design is largely a complex process by which specialized optical elements and their spatial relationships are chosen in order to minimize aberrations and provide an acceptable image resolution over a desired field of view (FOV) [29]. The more optical surfaces available to the designer, the greater the extent the aberrations can be minimized. However, this physical system improvement approach for minimizing aberrations has reached a point of diminishing returns in modern optics. Microscope objectives with 15 optical elements have become commercially available in recent years [32], but it is unlikely that another order of magnitude of optical surfaces will be supported within the confines of an objective in the foreseeable future. Moreover, this strategy 55 for minimizing aberration is never expected to accomplish the task of completely zeroing out aberrations. In other words, any optical system’s space-bandwidth product (SBP), which scales as the product of system FOV and inverse resolution, can be expected to remain a design bound dictated by the residual aberrations in the system. The issue of aberrations in simpler optical systems with few optical surfaces is, unsurprisingly, more pronounced. The eye is a very good example of such an optical system. While it does a fair job of conveying external scenes onto our retinal layer, its optical quality is actually quite poor. When a clinician desires a high resolution image of the retinal layer itself for diagnostic purposes, the human eye lens and cornea aberrations would have to be somehow corrected or compensated for. The prevalent approach by which this is currently done is through the use of adaptive optics (AO) [12, 49]. This is in effect a sophisticated way of physically correcting aberrations where complex physical optical elements are used to compensate for the aberrations of the lens and cornea. AO forms a guide star on the retina and uses a wavefront detector (e.g., Shack-Hartmann sensor) and a compensation device (e.g., deformable mirror) to correct for the aberrations affecting the guide star and a small region around it as it is under similar aberrations. This region is known as the isoplanatic patch [10] and its size varies depending on the severity of aberrations. To image a larger area beyond the isoplanatic patch, AO needs to be raster-scanned [3]. Since AO correction is fast (e.g., <500 ms [17]), it is still possible to obtain images of multiple isoplanatic patches quickly. However, the AO system can be complicated as it requires the active feedback loop between the wavefront measurement device and the compensation device and needs a separate guide star for the correction process [30]. The primary objective of this chapter is to report a novel generalized optical measurement system and computational approach to determine and correct aberrations in optical systems. It first reconstructs an optical system’s pupil function by adapting the FPM’s alternating projections as used in overlapped Fourier coding [18] to an incoherent imaging modality, which overcomes the spatial coherence requirement of the original pupil function recovery procedure of FPM. It can then recover the high resolution image latent in an aberrated image via deconvolution. The deconvolution is made robust to noise by using coded apertures to capture images [54]. This computational approach is coupled to a general optical scheme designed to efficiently collect the type of images required by the computational approach. We term this 56 method: coded-aperture-based correction of aberration obtained from overlapped Fourier coding and blur estimation (CACAO-FB). It is well-suited for various imaging scenarios where aberration is present and where providing a spatially coherent illumination is very challenging or impossible. CACAO-FB ultimately allows for an aberrated imaging system to achieve diffraction-limited performance over a wide FOV by casting optical design complexity to computational algorithms in postprocessing. CACAO-FB is substantially different from other recent efforts aimed at aberration compensation. Broadly speaking, these efforts can be divided into two major categories: blind and heuristic aberration recovery. Blind recovery minimizes a cost function, typically an image sharpness metric or a maximum-likelihood function, over a search space, usually the coefficient space of Zernike orthonormal basis [1, 9, 16, 39, 40, 43], to arrive at the optimal aberration function. However, blind recovery is prone to converging towards multiple local minima, and requires the aberrated sample image to be a complex field because blind aberration recovery with intensity-only sample image is extremely prone to noise for any aberrations [16] other than simple ones such as a camera motion blur or a defocus blur [25]. Heuristic recovery algorithms rely on several assumptions, such as assuming that the captured complex-field sample image has diffuse distribution in its Fourier spectrum such that each sub-region in the Fourier domain encodes the local aberrated wavefront information [11, 14, 23, 24]. Thus, heuristic methods are limited to specific types of samples and their performance is highly sample dependent. CACAO-FB is capable of achieving a robust aberration recovery performance in a generalized and broadly applicable format. In Section 5.2, we describe the principle of CACAO-FB. In Section 5.3, we report the demonstration of CACAO-FB with a crude lens and an eye model as imaging systems of interest. Finally, in Section 5.4, we demonstrate the potential of using CACAO-FB for retinal imaging in an in vivo experiment on a rhesus macaque’s eye, and discuss the current challenges it needs to address to become a viable alternative to other AO retinal imagers. We summarize our findings and discuss future directions in Section 5.5. 5.2 Principle of CACAO-FB To best understand the overall operation of the CACAO-FB processing, we start by examining the optical scheme (see Fig. 5.1). Suppose we start with an unknown optical system of interest (target system). This target system consists of a lens 57 Figure 5.1: Optical architecture of CACAO-FB. CACAO-FB system consists of 3 tube lenses (L1, L2, and L3) to relay the image from the target system for analysis. The target system consists of an unknown lens and an unknown sample with spatially incoherent field. CACAO-FB system has access to the conjugate plane of the target system’s pupil that can be arbitrarily modulated with binary patterns using a spatial light modulator. The images captured by CACAO-FB system is intensity-only. f0 , f1 , f2 , and f3 are the focal lengths of the unknown lens, L1, L2, and L3, respectively. d is an arbitrary distance smaller than f3 . (unknown lens) placed approximately at its focal length in front of a target sample (unknown sample). The sample is illuminated incoherently. For simplicity in this thought experiment, we will consider the illumination to occur in the transmission mode. The CACAO-FB system collects light from the target system using relay lenses L1, L2, and L3, and an aperture mask in the pupil plane, which is the conjugate plane of the target system’s pupil with coordinates (u, v), that can be modulated into different patterns. Our objective is to resolve the sample at high resolution. It should be clear from this target system description that our ability to achieve the objective is confounded by the presence of the unknown lens and its unknown aberrations. A good example of such a system is the eye: the retinal layer is the unknown sample, and the lens and cornea can be represented by the unknown lens. From this thought experiment, we can see that, to accomplish our objective, we would need to first determine the aberration characteristics of the unknown lens and then use the information to somehow correct out the aberration effects from the final rendered image. CACAO-FB does this by using 3 primary computational imaging algorithmic components that operate in sequence: 1) local aberration recovery with blur estimation, 2) full aberration recovery with FP-based alternating projections algorithm, and 3) latent image recovery by deconvolution with coded apertures. 58 The first two steps determine the target system’s aberrations, and the third step generates an aberration-corrected image. This pipeline is summarized in Fig. 5.2. The sample plane, which has coordinates (x, y), is divided into small tiles within which the aberration can be assumed to be spatially invariant, and CACAO-FB processes each corresponding tile on its image plane, which has coordinates (ξ, η), to recover a high resolution image of the sample tile. In the following analysis, we focus our attention to one tile, t. CACAO-FB begins by capturing a series of images with varying mask patterns in its pupil plane, which has coordinates (u, v). The patterns consist of two kinds: a set of small circular apertures, Wm (u, v), that collectively spans the pupil of the unknown lens; and a set of big apertures, An (u, v), that includes coded apertures and a full circular aperture with their diameters equal to the unknown lens’s pupil’s size. m and n are integers ranging from 1 to the total number of the respective aperture. The images captured with Wm (u, v) are labeled as i m,t (ξ, η), and they encode the local aberration of the unknown lens’s pupil function in their point spread functions (PSF). The blur estimation algorithm extracts these PSFs, bm,t (ξ, η). These intensity values of the spatially filtered pupil function can be synthesized into the full pupil function, Pt (u, v), with FP-based alternating projections algorithm. The images captured with An (u, v), labeled φn,t (ξ, η), are processed with the reconstructed pupil function and the knowledge of the mask patterns to generate the latent, aberration-free image of the sample, ot (x, y). The next 4 sub-sections will explain the mathematical model of the image acquisition process and the 3 imaging algorithmic components in detail. Image acquisition principle of CACAO-FB system We consider a point on the unknown sample, s(x, y), and how it propagates to the camera plane to be imaged. On the sample plane, a unit amplitude point source at (x0, y0 ) can be described by U0 (x, y; x0, y0 ) = δ(x − x0, y − y0 ), (5.1) where U0 (x, y; x0, y0 ) is the complex field of the point on the sample plane, and δ(x − x0, y − y0 ) is the Dirac delta function describing the point located at (x0, y0 ) We then use Fresnel propagation to propagate it to the unknown lens’s plane and apply the phase delay caused by the unknown lens assuming an idealized thin lens with the estimated focal length f0 (Eqs. (5.7-5.8) of Appendix). Any discrepancy from the ideal is incorporated into the pupil function, P(u, v; x0, y0 ), which is usually a circular bandpass filter with a uniform modulus and some phase modulation. Thus, 59 Figure 5.2: Outline of CACAO-FB pipeline. a) The captured images are broken into small tiles of isoplanatic patches (i.e., aberration is spatially invariant within each tile). b) Data acquisition and post-processing for estimating the pupil function, Pt (u, v). Limited-aperture images, i m,t (ξ, η), are captured with small masks, Wm (u, v), applied at the pupil plane. Local PSFs, bm,t (ξ, η), are determined by the blur estimation procedure, Algorithm 1. These PSFs are synthesized into the full-aperture pupil function, Pt (u, v), with Fourier-ptychography-based alternating projections algorithm, Algorithm 2. c) Data acquisition with big masks, An (u, v), at the pupil plane. d) The recovered Pt (u, v) from b) and the big-aperture images φn,t (ξ, η) from c) are used for deconvolution, Algorithm 3, to recover the latent aberration-free intensity distribution of the sample, ot (x, y). 60 the field right after passing through the unknown lens is   2π UF (u, v; x0, y0 ) = P(u, v; x0, y0 ) exp − j (x0 u + y0 v) , λ f0 (5.2) where λ is the wavelength of the field, and (u, v) are the coordinates of both the plane right after the unknown lens and the CACAO-FB system’s aperture plane as these planes are conjugate to each other. Thus, we refer to the aperture plane as the pupil plane. The spatially varying nature of a lens’s aberration is captured by the pupil function’s dependence on (x0, y0 ). We divide our sample into small tiles of isoplanatic patches (e.g., t = 1, 2, 3, ...) and confine our analysis to one tiled region, t, on the sample plane that contains (x0, y0 ) and other points in its vicinity such that the spatially varying aberration can be assumed to be constant in the analysis that follows (i.e., P(u, v; x0, y0 ) = Pt (u, v)). This is a common strategy for processing spatially variant aberration in a wide FOV imaging [15, 53]. We can see from Eq. (5.2) that the field emerging from the unknown lens is essentially its pupil function with additional phase gradient term defined by the point source’s location on the sample plane. At the pupil plane, a user-defined aperture mask, M(u, v), is applied to produce   2π 0 (x0 u + y0 v) , (5.3) UF (u, v; x0, y0 ) = M(u, v)Pt (u, v) exp − j λ f0 where we dropped the constant factor CF (x0, y0 ). After further propagation to the camera plane (Eqs. (5.10-5.14) of Appendix), we obtain the intensity pattern, i PSF,t (ξ, η), that describes the mapping of a point on the sample to the camera plane: i PSF,t (ξ, η; x0, y0 ) = |F {M(u, v)Pt (u, v)} (ξ, η) ∗ δ(ξ + = ht (ξ + x0 y0 ,η + ), λ f0 λ f0 y0 2 x0 ,η + )| λ f0 λ f0 (5.4) where ht (ξ, η) = |F {M(u, v)Pt (u, v)} (ξ, η)| 2 is the intensity of the PSF of the combined system in Fig. 5.1 for a given aperture mask M(u, v) and within the isoplanatic patch t. We observe from Eq. (5.4) that PSFs for different point source locations are related to each other by simple lateral shifts, such that an image, i t (ξ, η), captured by this system of an unknown sample function within the isoplanatic patch, st (x, y), can be represented by i t (ξ, η) = ht (ξ, η) ∗ |st (ξ, η)| 2 = ht (ξ, η) ∗ ot (ξ, η), (5.5) where ot (ξ, η) is the intensity of st (ξ, η), ∗ is the convolution operator, and we ignore the coordinate scaling for simplicity. This equation demonstrates that the image 61 captured by the detector is a convolution of the sample’s intensity field with a PSF associated with the sub-region of the pupil function defined by an arbitrary mask at the pupil plane. This insight allows us to capture images of the sample under the influence of PSFs that originate from different sub-regions of the pupil. We have aperture masks of varying shapes and sizes, mainly categorized into small masks and big masks. Small masks sample small regions of the pupil function to be used for reconstructing the pupil function, as will be described in detail in the following sections. Big masks include a mask corresponding to the full pupil size and several coded apertures that encode the pupil function to assist in the latent image recovery by deconvolution. To avoid confusion, we label the mth small mask, its associated PSF in isoplanatic patch t, and the image captured with it as Wm (u, v), bm,t (ξ, η), and i m,t (ξ, η), respectively; and the nth big mask (coded aperture or a full aperture), its associated PSF, and the image captured with it as An (u, v), hn,t (ξ, η), and φn,t (ξ, η), respectively. CACAO-FB system captures i m,t (ξ, η)’s and φn,t (ξ, η)’s in the data acquisition process, and these data are relayed to post-processing algorithms to recover ot (ξ, η), the underlying aberration-free image of the sample. The algorithm pipeline begins with the blur estimation algorithm using i m,t (ξ, η)’s as described below. In all the following simulations, there are one full aperture and 4 coded apertures, An (u, v), with the diameter of 4.5 mm; 64 small masks, Wm (u, v), with the diameter of 1 mm; an unknown lens, L1, and L2 with the focal length of f0 = f1 = f2 = 100 mm; a tube lens with the focal length of f3 = 200 mm; an image sensor with the pixel size of 6.5 µm (3.25 µm effective pixel size); and a spatially incoherent illumination with the wavelength of 520 nm. Local aberration recovery with blur estimation The blur function, bm,t (ξ, η), associated with the small mask, Wm (u, v), applied to the pupil, Pt (u, v), is also referred to as the local PSF, and it contains valuable information about the target system’s pupil function that we wish to recover. The size of Wm (u, v) is set small enough such that Wm (u, v) applied to a region on Pt (u, v) shows a local phase map that resembles a linear phase gradient, as shown in Fig. 5.3 b1). In such case, the associated bm,t (ξ, η) approximates a diffraction-limited spot with a spatial shift given by the phase gradient. Wm (u, v) applied to other regions on Pt (u, v) may have bm,t (ξ, η) whose shape deviates from a spot if the masked region contains more severe aberrations, as shown in Fig. 5.3 b2-3). In general, the aberration at or near the center of an imaging lens is minimal, and it becomes severe near the edge 62 Figure 5.3: Simulating image acquisition with different small masks at the pupil plane. a) The full pupil function masked by the lens’s NA-limited aperture. Differently masked regions of the pupil, b1-3), give rise to different blur kernels, c1-3), which allows us to capture images of the sample under the influence of different PSFs. Only the phase is plotted for Pt (u, v) and Pm,t (u, v)’s, and their apertures are marked by the black boundaries. W1 (u, v), W45 (u, v), and W52 (u, v) are 3 small masks from a spiraling-out scanning sequence. Figure 5.4: Flowchart of Algorithm 1: blur estimation algorithm for determining local PSFs from images captured with small apertures, Wm,t (u, v) 63 of the aperture because the lens’s design poorly approximates the parabolic shape away from the optical axis [13]. Thus, the image captured with the center mask, i 1,t (ξ, η), is mostly aberration-free with its PSF defined by the diffraction-limited spot associated with the mask’s aperture size. Other i m,t (ξ, η)’s have the same frequency bandlimit as i 1,t (ξ, η), but are under the influence of additional aberration encapsulated by their local PSFs, bm,t (ξ, η)’s. We adopt an image-pair-based blur estimation algorithm widely used in computational photography discipline to determine bm,t (ξ, η). In this algorithm, one of the image pair is assumed to be blur-free while the other is blurred [28, 51]. The blur kernel can be estimated by an iterative PSF estimation method, which is iterative Tikhonov deconvolution [35] in Fourier domain, adopting update scheme in Yuan’s blur estimation algorithm [51] and adjusting the step size to be proportional to |I1,t (u, v)|/|I1,t (u, v)|max for robustness to noise [38], where I1,t (u, v) is the Fourier spectrum of i 1,t (ξ, η). The blur estimation process is described in Algorithm 1 as shown in Fig. 5.4. The recovered bm,t (ξ, η)’s are the intensity information of the different masked pupil regions’ Fourier transforms. They can be synthesized into the full pupil function, Pt (u, v), using FP-based alternating projections algorithm, as described in the following section. Full aberration recovery with Fourier-ptychography-based alternating projections algorithm FP uses alternating projections algorithm to synthesize a sample’s Fourier spectrum from a series of intensity images of the sample captured by scanning an aperture on its Fourier spectrum [18, 52]. In our implementation, the full pupil’s complex field, Pt (u, v), is the desired Fourier spectrum to be synthesized, and the local PSFs, bm,t (ξ, η)’s, are the aperture-scanned intensity images to be used for FP-based alternating projections, as shown in the bottom half of Fig. 5.2 b). Therefore, reconstructing the pupil function from a series of local PSFs’ intensity information in our algorithm is completely analogous to reconstructing the complex spatial spectrum of a sample from a series of its low-passed images. The FP-based alternating projections algorithm is Algorithm 2, and it is described in Fig. 5.5. FP-based alternating projections algorithm requires that the scanned apertures during image acquisition have at least 30% overlap [41] for successful phase retrieval. Thus, the updating bm,t (ξ, η) in Algorithm 2 are ordered in a spiral-out pattern, each 64 Figure 5.5: Flowchart of Algorithm 2: Fourier-ptychography-based alternating projections algorithm for reconstructing the unknown lens’s pupil function, Pt (u, v) having an associated aperture Wm (u, v) that partially overlaps (40% by area) with the previous one’s aperture. The influence of the overlap on the reconstruction is illustrated in Fig. 5.17 of Appendix. For the simulated pupil diameter of 4.5 mm, there are 64 Wm (u, v)’s of 1-mm diameter to span the pupil with 40% overlap. We simulate the image acquisition by an aberrated imaging system and our pupil function reconstruction process in Fig. 5.6. Algorithm 1 and 2 are able to estimate the local PSFs from the 64 images captured with the small masks, Wm (u, v), (Fig. 5.6c)) and reconstruct the complex pupil function, Pt (u, v), successfully (Fig. 5.6e)). A simple Fourier transformation of Pt (u, v) generates the PSF of the aberrated imaging system. On a Macbook Pro with 2.5 GHz Intel Core i7 and 16 GB of RAM, it takes 2 minutes for Algorithm 1 and 20 seconds for Algorithm 2 to operate on the 64 images (1000 by 1000 pixels) taken with Wm (u, v). To gauge our method’s performance among other computational blur estimation methods, we attempt PSF reconstruction with two blind deconvolution algorithms. One is MATLAB’s deconvblind which is a standard blind deconvolution algorithm based on the accelerated, damped Richardson-Lucy algorithm, and the other is the state-of-the-art blind blur kernel recovery method based on variational Bayesian approach by Fergus et al. [8, 27]. They both operate on a single blurred image (Fig. 5.6d)) to simultaneously extract the blur function and the latent image (Fig. 5.6f)). For our purpose, we compare the reconstructed PSFs to gauge the performance. As shown in Fig. 5.6f), the reconstructed blur kernels by MATLAB and Fergus et al. both show poor fidelity to the 65 Figure 5.6: Simulation of our pupil function recovery procedure and a comparison with blind deconvolution algorithms. a) The Siemens star pattern used in the simulation. b) The system’s pupil function and the associated PSF. c) A series of images, i m,t (ξ, η)’s captured with small masks, Wm (u, v), applied to the pupil function. d) An image captured with the full-pupil-sized mask, An (u, v), on the pupil function, which simulates the general imaging scenario by an aberrated imaging system. e) The system’s pupil function and PSF recovered by our procedure. They show high fidelity to the original functions in b). f) Blur functions recovered by MATLAB’s and Fergus et al.’s blind deconvolution algorithm, respectively. They both show poor reconstructions compared to the recovered PSF in e). true PSF. This clearly demonstrates the effectiveness of our algorithm pipeline in reconstructing a complicated PSF which would otherwise be impossible to recover by a blind deconvolution method. The absolute limit of our aberration reconstruction method, assuming an unlimited photon budget, is essentially determined by the number of pixels inside the defined full aperture. However, in real life settings with limited photon budget and a dynamic sample, the smallest subaperture we can use to segment the full aperture is determined by the allowable exposure time and the shot-noise-limited condition of the camera. One has to consider the number of photons required by the camera for the signal to overcome the camera noise and the length of exposure permissible to capture a static image of the sample. Latent image recovery by deconvolution with coded apertures With the knowledge of the pupil function obtained from Algorithm 1 and 2, it is possible to recover ot (x, y) from the aberrated image φt (ξ, η) taken with the full 66 pupil aperture. In the Fourier domain, the image’s spectrum is represented as: Φt (u, v) = Ht (u, v)Ot (u, v), where Ht (u, v) and Ot (u, v) are the spatial spectrum of ht (ξ, η) and ot (x, y), respectively. Ht (u, v) is also called the optical transfer function (OTF) of the optical system and, by Fourier relation, is an auto-correlation of the pupil function, Pt (u, v). In the presence of severe aberrations, the OTF may have values at or close to zero for many spatial frequency regions within the bandpass, as shown in Fig. 5.7. These are due to the phase gradients with opposite slopes found in an aberrated pupil function, which may produce values at or close to zero in the autocorrelation process. Thus, the division of Φt (u, v) by Ht (u, v) during deconvolution will amplify noise at these spatial frequency regions since the information there has been lost in the image acquisition process. This is an ill-posed inverse problem. There are several deconvolution methods that attempt to address the ill-posed problem by using a regularizer [35] or a priori knowledge of the sample, such as by assuming sparsity in its total variation [2, 26]. However, due to their inherent assumptions, these methods work well only on a limited range of samples, and the parameters defining the a priori knowledge need to be manually tuned to produce successful results. Fundamentally, they do not have the information in the spatial frequency regions where the OTF is zero, and the a priori knowledge attempts to fill in the missing gaps. Wavefront coding using a phase mask in the Fourier plane has been demonstrated to remove the null regions in the OTF such that a subsequent deconvolution by the pre-calibrated PSF can recover the latent image [7, 22, 33, 34]. We adopt a similar method called coded aperture proposed by Zhou et al. [54] that uses an amplitude mask in the Fourier domain to achieve the same goal. With the amplitude-modulating SLM already in the optical system, using the amplitude mask over a phase mask is preferred. Combined with the knowledge of the pupil function reconstructed by Algorithm 1 and 2, no a priori knowledge is required to recover the latent image via deconvolution. A coded aperture designed by Zhou et al. at the pupil plane with a defocus aberration can generate a PSF whose OTF does not have zero values within its NA-limited bandpass. The particular coded aperture is generated by a genetic algorithm which searches for a binary mask pattern that maximizes its OTF’s spatial frequency content’s modulus. The optimum aperture’s pattern is different depending on the amount of noise in the imaging condition. We choose the pattern as shown in Fig. 5.7 since it performs well across various noise levels [54]. The pupil function in our imaging scenario does not only consist of defocus, as 67 Figure 5.7: Simulation that demonstrates the benefit of coded-aperture-based deconvolution. a1-5) Masked pupil functions obtained by masking the same pupil function with the full circular aperture and coded apertures under different rotation angles (0◦ , 45◦ , 90◦ , 135◦ ), their associated OTFs along one spatial frequency axis, and captured images. Each coded aperture is able to shift the null regions of the OTF to different locations. b) Comparison between the OTF of a circular-aperture-masked pupil function and the summed OTFs of the circular- and coded-aperture-masked pupil functions. Null regions in the frequency spectrum are mitigated in the summed OTF, which allows all the frequency content of the sample within the bandlimit to be captured with the imaging system. The OTF of an ideal pupil function is also plotted. c1) Deconvolved image with only a circular aperture shows poor recovery with artifacts corresponding to the missing frequency contents in the OTF’s null regions. c2) A recovered image using one coded aperture only. Reconstruction is better than c1), but still has some artifacts. c3) A recovered image using circular and multiple coded apertures is free of artifacts since it does not have missing frequency contents. 68 Figure 5.8: Flowchart of Algorithm 3: iterative Tikhonov regularization for recovering the latent sample image, ot (x, y), from aberrated images. Here, Φn,t (u, v) = F{φn,t (ξ, η)}(u, v). the imaging lenses have severe aberrations. Therefore, our pupil function can have an unsymmetrical phase profile unlike the defocus aberration’s symmetric bullseye phase profile. Thus, rotating the coded aperture can generate PSFs with different spatial frequency distribution, resulting in a different PSF shape beyond the mere rotation of the PSF. Therefore, in the data capturing procedure, we capture a series of images with a sequence of big masks, An (u, v), consisting of 4 coded apertures and a standard circular aperture at the pupil plane, as represented in Fig. 5.2 c). This ensures that we obtain all spatial frequency information within the NA-limited bandpass. The PSF associated with each An (u, v) applied to Pt (u, v) is easily obtained by hn,t (ξ, η) = |F −1 {An (u, v)Pt (u, v)} (ξ, η)| 2 , and its OTF by  Hn,t (u, v) = F hn,t (ξ, η) (u, v). With the measured full and coded aperture images, φn,t (ξ, η)’s, and the knowledge of the OTFs, we perform a combined deconvolution using iterative Tikhonov regularization, similar to Algorithm 1, to recover the object’s intensity distribution, ot (x, y), as described by Algorithm 3 in Fig. 5.8 and represented in Fig. 5.2 d). A simulated deconvolution procedure with the coded aperture on a Siemens star pattern is shown in Fig. 5.7. The combined OTF of a circular aperture and the coded aperture at 4 rotation angles is able to eliminate the null regions found in the circular-aperture-only OTF and thus produce a superior deconvolution result. The deconvolution performance across different frequency components is correlated to 69 the combined OTF’s modulus. We observe that the signal-to-noise ratio (SNR) of at least 40 in the initial aberrated image produces a deconvolution result with minimal artifacts. On a Macbook Pro with 2.5 GHz Intel Core i7 and 16 GB of RAM, it takes 2 seconds for Algorithm 3 to process the 5 images (1000 by 1000 pixels) taken with An (u, v) (1 full aperture, 4 coded apertures) to generate the deconvolution result. 5.3 Experimental demonstration on artificial samples Demonstration of CACAO-FB on a crude lens The CACAO-FB prototype system’s setup is simple, as shown in Fig. 5.9. It consists of a pair of 2-inch, f = 100 mm achromatic doublets (AC508-100-A from Thorlabs) to relay the surface of an imaging lens of interest to the surface of the ferroelectric liquid-crystal-on-silicon (LCOS) spatial light modulator (SLM) (SXGA-3DM-HB from 4DD). A polarized beam splitter (PBS) (PBS251 from Thorlabs) lies in front of the SLM to enable binary modulation of the SLM. A polarizer (LPVISE100-A from Thorlabs) is placed after the PBS to further filter the polarized light to compensate for the PBS’s low extinction ratio in reflection. A f = 200 mm tube lens (TTL200-A from Thorlabs) Fourier transforms the modulated light and images it on a sCMOS camera (PCOedge 5.5 CL from PCO). To determine the orientation of the SLM with respect to the Fourier space in our computational process, we us a phase-only target, such as a microbead sample, illuminated by a collimated laser source to perform an overlapped-Fourier-coding phase retrieval [18]. With the correct orientation, the reconstructed complex field should have the expected amplitude and phase. The imaging system to be surveyed is placed in front of the CACAO-FB system at the first relay lens’s focal length. The imaging system consists of a crude lens and a sample it is supposed to image. The crude lens in our experiment is a +6D trial lens (26 mm diameter, f = 130 mm) from an inexpensive trial lens set (TOWOO TW-104 TRIAL LENS SET). A resolution target is placed at less than the lens’s focal length away to introduce more aberration into the system. The sample is flood-illuminated by a monochromatic LED light source (520 nm, UHP-Microscope-LED-520 from Prizmatix) filtered with a 10-nm-bandpass filter. The relayed lens surface is modulated with various binary patterns by the SLM. The SLM displays a full aperture (5.5 mm diameter), a coded aperture rotated at 0◦ , 45◦ , 90◦ , and 135◦ with the maximum diameter matching the full aperture, and a series of limited apertures (1 mm diameter) shifted to different positions in a spiraling-out pattern within the full aperture dimension. The camera’s exposure is triggered by the SLM for synchronization. Another trigger signal for enabling the camera to begin a 70 Figure 5.9: Experimental setup of imaging a sample with a crude lens (i.e., unknown lens). Sample is illuminated by a monochromatic LED (520 nm), and the lens’ surface is imaged onto the SLM by a 1:1 lens relay. The part of light modulated by the SLM is reflected by the PBS and is further filtered by a polarizer to account for PBS’s low extinction ratio in reflection (1:20). Pupil-modulated image of the sample is captured on the sCMOS camera. L: lens; P: polarizer; PBS: polarizing beam splitter. capture sequence is provided by a data acquisition board (NI ELVIS II from National Instrument) which a user can control with MATLAB. Multiple images for each SLM aperture are captured and summed together to increase their signal-to-noise ratio (SNR). The full-aperture image has SNR=51, with other aperture-scanned images having SNR approximately proportional to the square root of their relative aperture area. SNR is estimated by calculating the mean and variance values in a uniformly bright patch on the image. To quantify the resolution performance of CACAO-FB, we image a Siemens star pattern with the crude +6D lens. The pattern consists of 40 line pairs radiating from the center such that the periodicity along a circle increases with increasing radius. The smallest circle along which the periodic structure is barely resolvable determines the resolution limit of the optical system [19]. For the focal length of 130 mm, the aperture diameter of 5.5 mm, and the illumination wavelength of 520 nm, the expected resolution is between λ/NA = 24.6 µm (coherent) and λ/(2NA) = 12.3 µm (incoherent) periodicity, defined by the spatial frequency cutoff in the coherent and incoherent transfer functions, respectively. As shown in Fig. 5.10, CACAO-FB can resolve features up to 19.6 µm periodicity, which is 71 Figure 5.10: Resolution performance measured by imaging a Siemens star target. a) A crude lens has optical aberration that prevents resolving the Siemens star’s features. b) CACAO-FB is able to computationally remove the aberration and resolve 19.6 µm periodicity feature size, which lies between the coherent and incoherent resolution limit given by the focal length of 130 mm, the aperture diameter of 5.5 mm, and the illumination wavelength of 520 nm. c) Pupil function recovered by CACAO-FB used for removing the aberration. d) The PSF associated with the pupil function. e) Intensity values from the circular traces on a) and b) that correspond to the minimum resolvable feature size of 19.6 µm periodicity. The Siemens star’s spokes are not visible in the raw image’s trace, whereas 40 cycles are clearly resolvable in the deconvolved result’s trace. within the expected resolution limit. The +6D lens is expected to have poor imaging performance that varies across its FOV since it is an inexpensive single element lens. We image a sample slide consisting of an array of USAF targets and select 3 tiled regions, each containing a USAF target pattern, to demonstrate CACAO-FB’s ability to address spatially variant aberration in its latent image recovery procedure. As shown in Fig. 5.11, the recovered PSFs in the three different regions are drastically different, which demonstrates the spatially variant nature of the aberration. Deconvolving each region with the corresponding PSF can successfully recover the latent image. The expected resolution limits as calculated above correspond to a range between Group 5 Element 3 (24.8 µm periodicity) and Group 6 Element 3 (12.4 µm periodicity). Features 72 Figure 5.11: Spatially varying aberrations compensation result on a grid of USAF target. a) the full FOV captured by our camera with the full circular aperture at 5.5 mm displayed on the SLM. Each small region denoted by b-d1) had a different aberration map as indicated by varying pupil function and PSFs. Spatially varying aberrations are adequately compensated in post-processing as shown by the deconvolution results b-d2). up to Group 5 Element 5 (19.68 µm periodicity) are resolved after deconvolution as shown in Fig. 5.11, which matches closely with the resolution determined by the Siemens star pattern. Demonstration of CACAO-FB on an eye model We use an eye model from Ocular Instruments to simulate an in vivo retinal imaging experiment. We embed a cut-out USAF resolution target (2015a USAF from Ready Optics) on the model’s retina, and fill the model’s chamber with de-ionized water, as shown in Fig. 5.12. We make adjustments to our CACAO-FB system as shown in Fig. 5.13 to accommodate the different imaging scenario. First, it has to illuminate the retina in a reflection geometry via the same optical path as imaging. A polarized beam splitter (PBS) is used to provide illumination such that the specular reflection from the eye’s cornea, which mostly maintains the s-polarization from the PBS, is filtered out of the imaging optical path. The scattered light from the retina is depolarized and can partially transmit through the PBS. The light source is a fiber-coupled laser diode (520 nm) (NUGM01T from DTR’s Laser Shop) which is made spatially incoherent by propagating through a 111-meter-long, 600-µm- 73 Figure 5.12: Eye model with a USAF target embedded on the retinal plane. a) A cut-out piece of glass of USAF target is attached on the retina of the eye model. The lid simulates the cornea and also houses a lens element behind it. b) The model is filled with water with no air bubbles in its optical path. c) The water-filled model is secured by screwing it in its case. core-diameter multimode fiber (FP600URT from Thorlabs), following the method in Ref. [37]. The laser diode is triggered such that it is on only during camera exposure. Images are captured at 50 Hz, ensuring that the flashing illumination’s frequency lies outside of the range that can cause photosensitive epilepsy in humans (i.e., between 15 and 20 Hz [31]). We add a pupil camera that outputs the image of the eye’s pupil with fiduciary marks for aligning the eye’s pupil with our SLM. Finally, a motion-reference camera (MRC) that has the identical optical path as the encoded-image camera (EIC) aside from pupil-modulation by SLM is added to the system to account for an in vivo eye’s motion between image frames. The amount of light split between the MRC and EIC can be controlled by the PBS and a quarter-wave plate before the SLM. In Fig. 5.14, the recovered images show severe spatially varying aberrations of the eye model, but good deconvolution performance throughout the FOV, nonetheless. The tile size is set such that it is the biggest tile that could produce an aberration-free image, judged visually. The full aperture in this scenario had 4.5-mm diameter, and its associated aberrated image had SNR of 126. In this imaging scenario, the blur kernels of the limited-aperture images had a significant impact on the deconvolution result, as shown in Fig. 5.15. The aberration of the eye model was severe such that the retrieved blur kernels of the limited-aperture images had distinct shapes in addition to lateral shifts. We observe a much better deconvolution result with the reconstructed pupil that takes blur kernels’ shapes into account compared to the one that does not. The latter is analogous to ShackHartmann wavefront sensing method, which only identifies the centroid of each blur kernel to estimate the aberration. Thus, this demonstrates the importance of the blur 74 Figure 5.13: Experimental setup of imaging an eye model and an in vivo eye. Illumination is provided by a fiber-coupled laser diode (520 nm), and the eye’s pupil is imaged onto the SLM by a 1:1 lens relay. The sample is slightly defocused from the focal length of the crude lens to add additional aberration into the system. Pupil alignment camera provides fiduciary to the user for adequate alignment of the pupil on the SLM. PBS2 helps with removing corneal reflection. Motion-reference camera is synchronized with encoded-image camera to capture images not modulated by the SLM. BS: beam splitter; L: lens; M: mirror; P: polarizer; PBS: polarized beam splitter; QWP: quarter-wave plate. kernel estimation step in our algorithm and the distinct difference of our aberration reconstruction from other wavefront sensing methods. 5.4 Adapting CACAO-FB to an in vivo experiment on the eye of a rhesus macaque The same setup as for imaging the model eye is used for the in vivo experiment on a rhesus macaque’s eye. The animal is anesthetized with 8-10 mg/kg ketamine and 0.02 mg/kg dexdomitor IM. 2 drops of tropicamide (0.5-1%) are placed on the eye to dilate the pupil. To keep the eye open for imaging, a sanitized speculum is be placed between the eyelids. A topical anesthetic (proparacaine 0.5%) is applied to the eye to prevent any irritation from the speculum placement. A rigid gas permeable lens is placed on the eye to ensure that the cornea stays moist throughout imaging. The light intensity is kept below a level of 50 mW/cm2 on the retina in accordance with ANSI recommended safe light dosage. 75 Figure 5.14: CACAO-FB result of imaging the USAF target in the eye model. a) Raw image (2560 x 1080 pixels) averaged over 12 frames captured with the full circular aperture at 4.5 mm. The pupil function and PSF in each boxed region show the spatially varying aberrations. b-d) Deconvolution results show sharp features of the USAF target. The uneven background is from the rough surface of the eye model’s retina. Due to the safety limitation on the illumination power, the captured images of the retina has low SNR (e.g., the full-aperture image has SNR = 7.5). We increase the SNR by capturing multiple redundant frames (213 frames for An (u, v)’s, 4 frames for Wm (u, v)’s) and adding them together (Fig. 5.19 of Appendix). Thus, a long sequence of images (∼45 seconds) has to be captured, and these images with weak retinal signals have to be registered for motion prior to CACAO-FB since the eye has residual motion even under anesthesia. Due to a long averaging window, aberration of high temporal frequency is washed out, but we still expect to be able to resolve the photoreceptors albeit at a lower contrast [17]. Motion registration includes both translation and rotation, and these operations needs to be done such that they do not apply any spatial filter that may alter the images’ spatial frequency spectra. Rotation is performed with fast discrete sinc-interpolation [50] which is a series of Fourier transform operations that can be accelerated by GPU programming. The frames from the motion-reference camera are used for the registration process (Fig. 5.18 of 76 Figure 5.15: Showing the importance of masked pupil kernel shape determination for successful deconvolution. a1-3) limited PSFs determined only by considering their centroids. b) Recovered aberration and deconvolution result obtained with centroid-only limited PSFs. Some features of USAF are distorted. c1-3) limited PSFs determined with the blur estimation algorithm. d) Recovered aberration and deconvolution result obtained with the blur-estimated local PSFs. No distortions in the image are present, and more features of the USAF target are resolved. Appendix). A center region with half the dimensions of the full frame is selected from one of the frames as a template for registration. Normalized cross-correlation (NCC) value is found between the template and each frame [47] for every rotation angle (-1.5 to 1.5 degrees, 0.0015 degrees step size). The set of rotation angle and lateral shift values that produces the maximum NCC value, at a pixel resolution, for each frame corresponds to the motion registration parameters for that frame and the corresponding encoded-image camera’s frame. The registration parameters for all the frames are applied to the images of the encoded-image camera, and the images are grouped by different apertures to be summed together (Fig. 5.19 of Appendix). The deconvolution result is shown in Fig. 5.16. Although the sensor size of 2560 x 2160 pixels is used for capturing raw images, the resultant averaged images are only 1262 x 1614 pixels after the motion registration. The input full-aperture image had SNR=109. Photoreceptors are much better resolved after aberration removal. We expect the entire visible region to have an even spread of photoreceptors, but we observe well-resolved photoreceptors mostly in the brighter regions. This may be 77 Figure 5.16: CACAO-FB result from imaging an in vivo eye of a rhesus macaque. a) Raw image averaged over 213 frames captured with 4.5 mm full circular aperture. b) Deconvolution result using the c) pupil function reconstructed by CACAO-FB procedure. d) PSF associated with the pupil function. due to the lower SNR in the darker regions leading to poorer deconvolution result. We cannot capture more frames of the retina to further increase the SNR because the animal’s eye’s gaze drifts over time and the original visible patch of the retina is shifted out of our system’s FOV. Furthermore, non-uniform specular reflections from the retina add noise to part of the captured data, leading to sub-optimal latent image recovery by CACAO-FB algorithm pipeline. 5.5 Discussion We developed a novel method to characterize the aberration of an imaging system and recover an aberration-free latent image of an underlying sample in post-processing. It does not require the coherent illumination necessary in other computational, aberration-compensating imaging methods. It does not need separate wavefront detection and correction devices found in many conventional adaptive optics systems as the hardware complexities are off-loaded to the software regime that can harness the ever-increasing computational power. Its principle is based on incoherent imaging, which obviates sensitivity issues such as phase fluctuations and incident angles associated with coherent imaging, and allows for characterizing an integrated optical system where the sample plane is only accessible via the target system’s lens. Our demonstrations of CACAO-FB on sub-optimal lenses in benchtop and in vivo experiments show its viability in a broad range of imaging scenarios. Its simple 78 hardware setup is also a key advantage over other aberration-correction methods that may allow for its wide adoption. More severe aberration can be addressed readily by shrinking the scanned aperture size on the SLM so that the aberration within each windowed pupil function remains low-order. This comes at the expense of the acquisition speed as more images need to be captured to cover the same pupil diameter. If the masks on the pupil are shrunk smaller with no overlap, this pupil masking process becomes a Shack-Hartmann (SH) sensing method. This illustrates the key advantages of our scheme over a SH sensor: using bigger masks allows for a fewer number of image acquisitions and increases the images’ SNR. A bigger mask of an aberrated pupil no longer encodes for a simple shifted spot in the spatial domain as would be the case for a SH sensor, but rather a blur kernel as shown in Fig. 5.3. Therefore, reconstructing the blur kernels of the limited aperture images is critical for CACAO-FB’s performance, as is demonstrated in Fig. 5.15. Using an aperture mask in the Fourier plane discards a significant amount of photons in the image acquisition process. One possible way to improve the photon efficiency of our system would be to use a phase mask instead of an amplitude mask to code the Fourier plane as has been demonstrated in Ref. [34] to remove nulls in the OTF of an aberrated imaging system. Although the recovered retinal image in Section 5.4 is not on par with what one can achieve with a typical AO retinal imager, it showcases the proof of concept of using CACAO-FB to correct for aberrations in a general optical system. There are several challenges of imaging an eye in vivo eye that can be addressed in future works to allow CACAO-FB to be a viable alternative to AO retinal imagers. First, increasing the SNR by averaging multiple retinal images of the rhesus macaque’s eye in vivo is challenging as its gaze continues to drift even under general anesthesia. There is a finite number of frames we can capture before the original patch of retina shifts out of our system’s FOV. Imaging a human subject would be less susceptible to this issue as an operator can instruct the subject to focus on a target and maintain the same patch of the retina within the system’s FOV as done in Ref. [17]. The small lateral shifts between captured frames due to the eye’s saccade can be digitally registered prior to averaging. Using a different wavelength invisible to the eye will allow the subject to maintain his/her gaze throughout an extended acquisition time. Second, non-uniform specular reflections from the retinal layer corrupt the captured images. The flood illumination provided on the retina through the pupil does not 79 have sufficient angular coverage to even out the specular reflections such that some images captured with a small aperture contained specular reflections while others did not. A flood illumination with a wider divergence angle can mitigate this problem. On the other hand, our method will be able to provide a readily viable solution for imaging static targets such as wide FOV imaging of fluorescent samples under sample- and system-induced aberrations since the fluorescent signals are inherently spatially incoherent and CACAO-FB principle can be applied to the imaging process. With its simple optical setup, we believe CACAO-FB can be easily incorporated into many existing imaging systems to compensate for the limitations of the physical optics design. 5.6 Appendix Mathematical derivation of capturing images with pupil plane modulation We consider a point on the unknown sample, s(x, y), and how it propagates to the camera plane to be imaged. On the sample plane, a point source at (x0, y0 ) may have an amplitude and phase C, and it can be described by U0 (x, y; x0, y0 ) = Cδ(x − x0, y − y0 ). (5.6) We then use Fresnel propagation: U1 (u, v; x0, y0 ) i h exp j λπf0 (u2 + v 2 ) ∫ ∞ = Cδ(x − x0, y − y0 ) jλ f0 −∞     π 2 2π 2 (x + y ) exp − j (xu + yv) dxdy exp j λ f0 λ f0 i h   exp j λπf0 (u2 + v 2 ) π 2 2 exp j (x + y0 ) =C jλ f0 λ f0 0   2π exp − j (x0 u + y0 v) λ f0 (5.7) and apply the phase delay associated with an idealized thin lens [13] h i having an estimated focal length f0 for the unknown lens, exp − j λπf0 (u2 + v 2 ) , and any discrepancy from the ideal is incorporated into the pupil function, P(u, v; x0, y0 ) = Pt (u, v), 80 where t = 0 is the isoplanatic patch around (x0, y0 ): U2 (u, v; x0, y0 )     C π 2 2π 2 = exp j (x + y0 ) exp − j (x0 u + y0 v) Pt (u, v) jλ f0 λ f0 0 λ f0   2π = C2 (x0, y0 )Pt (u, v) exp − j (x0 u + y0 v) , (5.8) λ f0   where we set C2 (x0, y0 ) = C/( jλ f0 ) exp jπ/(λ f0 )(x02 + y02 ) . Eq. (5.8) is Eq. (2) of the main article. U2 (u, v; x0, y0 ) is relayed to the pupil plane in Fig. 1 without any additional phase term by the 4f system formed by L1 and L2 [13]. The relayed field may be magnified by the factor f2 / f1 , but here we assume no magnification. We apply a mask, M(u, v), to the field: U20 (u, v; x0, y0 )  2π (x0 u + y0 v) = M(u, v)C2 (x0, y0 )Pt (u, v) exp − j λ f0  (5.9) and propagate it by distance d using angular spectrum to the surface of L3:  U3 (s, t; x0, y0 ) = F −1 {F U20 (u, v; x0, y0 )   q 2πd 2 2 exp − j 1 − (λ fx ) − (λ fy ) }(s, t), (5.10) λ where (s, t) are the coordinates on the iL3’s plane. We then apply the phase delay h π associated with L3, exp − j λ f3 (s2 + t 2 ) , and propagate the field by f3 to the camera plane: U4 (ξ, η; x0, y0 ) i h   exp j λπf3 (ξ 2 + η2 ) ∫ ∞ π 2 2 (s + t ) = exp − j jλ f3 λ f3 −∞     π 2 2 2π U3 (s, t; x0, y0 ) exp j (s + t ) exp − j (sξ + tη) dsdt λ f3 λ f3 h i exp j λπf3 (ξ 2 + η2 ) ∫ ∞ = F −1 {F{U20 (u, v; x0, y0 )} jλ f3 −∞   q 2πd exp − j 1 − (λ fx )2 − (λ fx )2 }(s, t) λ   2π exp − j (sξ + tη) dsdt. λ f3 (5.11) 81  Set C4 (ξ 0, η0) = exp jπλ f3 (ξ 02 + η02 ) /( jλ f3 ) and (ξ 0, η0) = (ξ, η)/(λ f3 ):  U40 (ξ 0, η0; x0, y0 ) = C4 (ξ 0, η0)F{F −1 {F{U20 (u, v; x0, y0 )}   q 2πd 2 2 exp − j 1 − (λ fx ) − (λ fy ) }}(ξ 0, η0) λ = C4 (ξ 0, η0)F{U20 (u, v; x0, y0 )∗   q 2πd −1 F {exp − j 1 − (λ fx )2 − (λ fy )2 }}(ξ 0, η0) λ = C4 (ξ 0, η0)F{U20 (u, v; x0, y0 )}(ξ 0, η0)   q 2πd 0 2 0 2 exp − j 1 − (λξ ) − (λη ) λ   q 2πd 0 0 = C4 (ξ , η ) exp − j 1 − (λξ 0)2 − (λη0)2 λ   2π (x0 u + y0 v) }(ξ 0, η0), F{M(u, v)C2 (x0, y0 )Pt (u, v) exp − j λ f0 (5.12) where ∗ is the convolution operator. h i p 0)2 − (λη0)2 C (x , y ) as C (ξ 0, η0; x , y ): Set C4 (ξ 0, η0) exp − j 2πd 1 − (λξ 2 0 0 0 0 5 λ U40 (ξ 0, η0; x0, y0 ) = C5 (ξ 0, η0; x0, y0 )[F{M(u, v)Pt (u, v)}(ξ 0, η0) x0 0 y0 ∗ δ(ξ 0 + ,η + )] λ f0 λ f0 y0 x0 0 ,η + ). = C5 (ξ 0, η0; x0, y0 )F{M(u, v)Pt (u, v)}(ξ 0 + λ f0 λ f0 (5.13) This is the complex field incident on the camera from the point source located at (x0, y0 ). It is the PSF of the system, and we observe that it simply shifts laterally for different (x, y) coordinates. Therefore, the image on the camera sensor can be calculated by a convolution between U40 (ξ 0, η0; x0, y0 ) and the sample field within the isoplanatic patch, st (x, y). However, the phase term in C5 (ξ 0, η0; x0, y0 ) can have a significant impact on the captured images. In our incoherent imaging scenario, the phase relationship between the points on the sample plane during the capturing process is irrelevant. So, we can define an intensity PSF, ht (ξ, η) = |F {M(u, v)Pt (u, v)} (ξ, η)| 2 , 82 Figure 5.17: The influence of the limited masks’ overlap ratio on the pupil function recovery. The recovery becomes poor below 30% overlap. to describe the intensity of the U40 (ξ 0, η0; x0, y0 ) captured by the camera: U40 (ξ 0, η0; x0, y0 ) 2 y0 x0 0 ,η + ) λ f0 λ f0 C y0 x0 0 = 2 ,η + ). ht (ξ 0 + λ f0 λ f0 λ f0 f3 = |C5 (ξ 0, η0; x0, y0 )| 2 ht (ξ 0 + (5.14) The complicated phase fluctuations embedded in C5 (ξ 0, η0; x0, y0 ) become no longer relevant. Dropping the constants and neglecting coordinate scaling, the image of the unknown sample captured by the camera becomes a convolution of ht (ξ, η) with the sample distribution, st , as described by Eq. (5) of the main article. Limited aperture overlap requirement for pupil function reconstruction Given the ground truth pupil function 5 mm in diameter and the limited mask diameter of 1 mm, the area overlap between contiguous masks are varied from 10% to 60% in simulation. In spatial domain, all images are captured satisfying the Nyquist criterion. The simulation is identical to Fig. 6’s data acquisition and reconstruction using Algorithm 1 and 2. As shown in Fig. 5.17, the recovery becomes poor below 30% overlap as observed by disconnected fringe patterns, and eventually contains holes in regions where there was no mask coverage. 83 Figure 5.18: Motion-reference camera images registered for rotation and translation. Only 4 representative frames are shown. Figure 5.19: Raw full aperture image. a) 1 frame and b) a sum of 213 frames. Increasing SNR by averaging over multiple frames In processing the data from the in-vivo experiment, motion-reference camera images are first registered for rotation and translation, as shown in Fig. 5.18, and these registration values are taken into account when summing multiple frames captured with the same SLM aperture pattern. An example of a single frame of full aperture image is shown in Fig. 5.19, and the same aperture image after summing 213 frames. Due to the photon-starved imaging condition, it is imperative to account for the detector noise in the captured images. We use two-point radiometric calibration to account for the fixed pattern noise and inhomogeneous sensitivity of our imaging 84 sensor [20]: I 0(ξ, η) = I(ξ, η) − B(ξ, η) , R(ξ, η) − B(ξ, η) (5.15) where (ξ, η) are the coordinates on the camera sensor’s plane, I 0(ξ, η) is the desired calibrated image, I(ξ, η) is the input image, B(ξ, η) is the dark image captured with the sensor blocked from light, and R(ξ, η) is a reference image captured with the sensor capturing an image of an opal diffuser. References [1] S. Adie, B. Graf, A. Ahmad, S. Carney, and S. Boppart. Computational adaptive optics for broadband optical interferometric tomography of biological tissue. Proceedings of the National Academy of Sciences, 109(19):7175–7180, 2012. [2] J. M. Bioucas-Dias, M. A. T. Figueiredo, and J. P. Oliveira. Total variationbased image deconvolution: a majorization-minimization approach. In Proceedings of IEEE International Conference on Acoustics Speech and Signal Processing, pages 861–864. IEEE, 2006. [3] M. J. Booth. Adaptive optical microscopy: the ongoing quest for a perfect image. Light: Science and Applications, 3:e165, 2014. [4] J. Chung, J. Kim, X. Ou, R. Horstmeyer, and C. Yang. Wide field-of-view fluorescence image deconvolution with aberration-estimation from Fourier ptychography. Biomedical Optics Express, 7(2):352–368, 2016. [5] J. Chung, H. Lu, X. Ou, H. Zhou, and C. Yang. Wide-field Fourier ptychographic microscopy using laser illumination source. Biomedical Optics Express, 7(11):4787–4802, 2016. [6] S. Dong, K. Guo, P. Nanda, R. Shiradkar, and G. Zheng. FPscope: a fieldportable high-resolution microscope using a cellphone lens. Optics Express, 5 (10):3305–3310, 2014. [7] E. Dowski and T. W. Cathey. Extended depth of field through wavefront coding. Applied Optics, 34(11):1859–1866, 1995. [8] R. Fergus, B. Singh, A. Hertzmann, S. T. Roweis, and W. T. Freeman. Removing Camera Shake from a Single Photograph. ACM Transactions on Graphics, 25(3):787–794, 2006. [9] J. R. Fienup and J. J. Miller. Aberration correction by maximizing generalized sharpness metrics. Journal of the Optical Society of America A, 20(4):609–620, 2003. 85 [10] D. L. Fried. Anisoplanatism in adaptive optics. Journal of the Optical Society of America, 72(1):52–61, 1982. [11] L. Ginner, T. Schmoll, A. Kumar, M. Salas, N. Pricoupenko, L. Wurster, and R. Leitgeb. Holographic line field en-face OCT with digital adaptive optics in the retina in vivo. Biomedical Optics Express, 9(2):472–485, 2018. [12] P. Godara, A. Dubis, A. Roorda, J. Duncan, and J. Carroll. Adaptive optics retinal imaging: Emerging clinical applications. Optometry and Vision Science, 87(12):930–941, 2010. [13] J. Goodman. Introduction to Fourier Optics. McGraw-Hill, 2008. [14] G. Gunjala, S. Sherwin, A. Shanker, and L. Waller. Aberration recovery by imaging a weak diffuser. Optics Express, 26(16):21054–21068, 2018. [15] B. K. Gunturk and X. Li. Image Restoration: Fundamentals and Advances, volume 7. CRC Press, 2012. [16] D. Hillmann, H. Spahr, C. Hain, H. Sudkamp, G. Franke, C. Pfåffle, C. Winter, and G. Hüttmann. Aberration-free volumetric high-speed imaging of in vivo retina. Scientific Reports, 6:35209, 2016. [17] H. Hofer, L. Chen, G. Y. Yoon, B. Singer, Y. Yamauchi, and D. R. Williams. Improvement in retinal image quality with dynamic correction of the eye’s aberrations. Optics Express, 8(11):631–643, 2001. [18] R. Horstmeyer, X. Ou, J. Chung, G. Zheng, and C. Yang. Overlapped Fourier coding for optical aberration removal. Optics Express, 22(20):24062–24080, 2014. [19] R. Horstmeyer, R. Heintzmann, G. Popescu, L. Waller, and C. Yang. Standardizing the resolution claims for coherent microscopy. Nature Photonics, 10:68–71, 2016. [20] B. Jähne. Digital Image Processing. Springer, 2005. [21] C. Kuang, Y. Ma, R. Zhou, J. Lee, G. Barbastathis, R. R. Dasari, Z. Yaqoob, and P. T. C. So. Digital micromirror device-based laser-illumination Fourier ptychographic microscopy. Optics Express, 23(21):26999–27010, 2015. [22] K. Kubala, E. R. Dowski, and W. T. Cathey. Reducing complexity in computational imaging systems. Optics Express, 11(18):2102–2108, 2003. [23] A. Kumar, W. Drexler, and R. A. Leitgeb. Subaperture correlation based digital adaptive optics for full field optical coherence tomography. Optics Express, 21 (9):10850–10866, 2013. 86 [24] A. Kumar, D. Fechtig, L. Wurster, L. Ginner, M. Salas, M. Pircher, and R. Leitgeb. Noniterative digital aberration correction for cellular resolution retinal optical coherence tomography in vivo. Optica, 4(8):924–931, 2017. [25] D. Kundur and D. Hatzinakos. Blind image deconvolution. IEEE Signal Processing Magazine, 13(2):43–64, 1996. [26] A. Levin, R. Fergus, F. Durand, and W. T. Freeman. Image and depth from a conventional camera with a coded aperture. ACM Transactions on Graphics (TOG), 26(3):70, 2007. [27] A. Levin, Y. Weiss, F. Durand, and W. T. Freeman. Understanding Blind Deconvolution Algorithms. IEEE Transactions on Pattern Analysis and Machine Intelligence, 33(12):2354–2367, 2011. [28] S. H. Lim and D. A. Silverstein. Method for deblurring an image. U.S. patent 8,654,201 (February 18, 2014). [29] A. W. Lohmann, R. G. Dorsch, D. Mendlovic, Z. Zalevsky, and C. Ferreira. Space-bandwidth product of optical signals and systems. J. Opt. Soc. Am. A, 13(3):470–473, 1996. [30] S. Marcos, J. S. Werner, S. A. Burns, W. H. Merigan, P. Artal, D. A. Atchison, K. M. Hampson, R. Legras, L. Lundstrom, G. Yoon, J. Carroll, S. S. Choi, N. Doble, A. M. Dubis, A. Dubra, A. Elsner, R. Jonnal, D. T. Miller, M. Paques, H. E. Smithson, L. K. Young, Y. Zhang, M. Campbell, J. Hunter, A. Metha, G. Palczewska, J. Schallek, and L. C. Sincich. Vision science and adaptive optics, the state of the field. Vis. Res., 132:3–33, 2017. [31] A. Martins da Silva and B. Leal. Photosensitivity and epilepsy: Current concepts and perspectives–A narrative review. Seizure, 50:209–218, 2017. [32] G. McConnell, J. Trägårdh, R. Amor, J. Dempster, E. Reid, and W. B. Amos. A novel optical microscope for imaging large embryos and tissue volumes with sub-cellular resolution throughout. eLife, 5:e18659, 2016. [33] G. Muyo and A. R. Harvey. Wavefront coding for athermalization of infrared imaging systems. Proc. SPIE, 5612:227–235, 2004. [34] G. Muyo, A. Singh, M. Andersson, D. Huckridge, A. Wood, and A. R. Harvey. Infrared imaging with a wavefront-coded singlet lens. Optics Express, 17(23): 21118–21123, 2009. [35] A. Neumaier. Solving ill-conditioned and singular linear systems: A tutorial on regularization. SIAM Rev., 40(3):636–666, 1998. [36] X. Ou, G. Zheng, and C. Yang. Embedded pupil function recovery for Fourier ptychographic microscopy. Optics Express, 22(5):4960–4972, 2014. 87 [37] J. Rha, R. S. Jonnal, K. E. Thorn, J. Qu, Y. Zhang, and D. T. Miller. Adaptive optics flood-illumination camera for high speed retinal imaging. Optics Express, 14(10):4552–4569, 2006. [38] J. M. Rodenburg and H. M. L. Faulkner. A phase retrieval algorithm for shifting illumination. Appl. Phys. Lett., 85(20):4795–4797, 2004. [39] N. Shemonski, F. South, Y. Z. Liu, S. Adie, S. Carney, and S. Boppart. Computational high-resolution optical imaging of the living human retina. Nature Photonics, 9(7):440–443, 2015. [40] F. Soulez, L. Denis, Y. Tourneur, and É. Thiébaut. Blind deconvolution of 3d data in wide field fluorescence microscopy. In Proceedings of 9th IEEE International Symposium on Biomedical Imaging (ISBI), pages 1735–1738. IEEE, 2012. [41] J. Sun, Q. Chen, Y. Zhang, and C. Zuo. Sampling criteria for Fourier ptychographic microscopy in object space and frequency space. Optics Express, 24 (14):15765–15781, 2016. [42] J. Sun, C. Zuo, L. Zhang, and Q. Chen. Resolution-enhanced Fourier ptychographic microscopy based on high-numerical-aperture illuminations. Scientific Reports, 7:1187, 2017. [43] É. Thiébaut and J. M. Conan. Strict a priori constraints for maximum-likelihood blind deconvolution. J. Opt. Soc. Am. A, 12(3):485–492, 1995. [44] 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. [45] L. Tian, X. Li, K. Ramchandran, and L. Waller. Multiplexed coded illumination for Fourier Ptychography with an LED array microscope. Biomedical Optics Express, 5(7):2376–2389, 2014. [46] L. Tian, Z. Liu, L. H. Yeh, M. Chen, J. Zhong, and L. Waller. Computational illumination for high-speed in vitro Fourier ptychographic microscopy. Optica, 2(10):904–911, 2015. [47] A. R. Wade and F. W. Fitzke. A fast, robust pattern recognition system for low light level image registration and its application to retinal imaging. Optics Express, 3(5):190–197, 1998. [48] A. Williams, J. Chung, X. Ou, G. Zheng, S. Rawal, Z. Ao, R. Datar, C. Yang, and R. Cote. Fourier ptychographic microscopy for filtration-based circulating tumor cell enumeration and analysis. J. Biomed. Opt., 19(6):066007, 2014. [49] D. Williams. Imaging single cells in the living retina. Vision Research, 51 (13):1379–1396, 2011. 88 [50] L. Yaroslavsky. Theoretical Foundations of Digital Imaging Using MATLAB. CRC Press, 2013. [51] L. Yuan, J. Sun, L. Quan, and H. Y. Shum. Image deblurring with blurred/noisy image pairs. ACM Transactions on Graphics, 26(3):1, 2007. [52] G. Zheng, R. Horstmeyer, and C. Yang. Wide-field, high-resolution Fourier ptychographic microscopy. Nature Photonics, 7(9):739–745, 2013. [53] G. Zheng, X. Ou, R. Horstmeyer, and C. Yang. Characterization of spatially varying aberrations for wide field-of-view microscopy. Optics Express, 21 (13):15131–15143, 2013. [54] C. Zhou and S. Nayar. What are good apertures for defocus deblurring? In Proceedings of IEEE International Conference on Computational Photography, pages 1–8. IEEE, 2009.