## Abstract

We report a novel generalized optical measurement system and computational approach to determine and correct aberrations in optical systems. The system consists of a computational imaging method capable of reconstructing an optical system’s pupil function by adapting overlapped Fourier coding to an incoherent imaging modality. It recovers 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. We term this 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. We report the demonstration of CACAO-FB with a variety of samples including an *in vivo* imaging experiment on the eye of a rhesus macaque to correct for its inherent aberration in the rendered retinal images. CACAO-FB ultimately allows for an aberrated imaging system to achieve diffraction-limited performance over a wide field of view by casting optical design complexity to computational algorithms in post-processing.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 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. 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) [1]. 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 [2], 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 spatial 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) [3,4]. 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 [5], 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 [6]. Since AO correction is fast (e.g., $<500\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{ms}$ [7]), 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 [8].

Fourier ptychography (FP) circumvents the challenges of adding more optical elements for improving an optical system’s performance by recasting the problem of increasing the system’s spatial 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, FP 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 [9–15]. FP shares its roots with ptychography [16–18] and structured illumination microscopy (SIM) [19–21], which numerically expand the SBP of the imaging system in the spatial and spatial frequency domain, respectively, by capturing multiple images under distinct illumination patterns and computationally synthesizing them into a higher-SBP image. One way to view FP is to note its similarity to synthetic aperture imaging [22–24]. In a standard FP microscope system, images of the target are collected through a low-numerical-aperture (NA) objective with the target illuminated with a series of angularly varied planar or quasi-planar illumination. Viewed in the spatial frequency domain, each image represents a disc of information with its offset from the origin determined by the illumination angle. 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 FP, we employ phase retrieval and the partial information overlap among the image set to converge on the correct phase relationships during the stitching process [9]. At the end of the process, the constituted information in the spatial frequency domain can be Fourier transformed to generate a higher-resolution image of the target that retains the original FOV as set by the objective. It has been demonstrated that a sub-routine can be weaved into the primary FP algorithm that will dynamically determine the pupil function of the imaging system [11]. In fact, the majority of existing FP algorithms incorporate some versions of this aberration determination function to find and subsequently correct out the aberrations from the processed image [25–32]. This particular sub-discipline of FP 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 [33]—this drives home the fact that correcting aberration computationally is a viable alternative to physical correction.

The primary objective of this paper is to report a novel generalized optical measurement system and computational approach to determine and correct aberrations in optical systems. This computational approach is coupled to a general optical scheme designed to efficiently collect the type of images required by the computational approach. Currently, FP’s ability to determine and correct aberration is limited to optical setups with a well-defined, spatially coherent field on the sample plane [9,11,34–42]. We developed a computational imaging method capable of reconstructing an optical system’s pupil function by adapting the FP’s alternating projections as used in overlapped Fourier coding [10] to an incoherent imaging modality, which overcomes the spatial coherence requirement of the original pupil function recovery procedure of FP. 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 [43]. We term this 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 post-processing.

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).

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 [44–49], 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 an intensity-only sample image is extremely prone to noise for any aberrations [45] other than simple ones such as a camera motion blur or a defocus blur [50]. 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 [51–54]. 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 2, we describe the principle of CACAO-FB. In Section 3, we report the demonstration of CACAO-FB with a crude lens and an eye model as imaging systems of interest. Finally, in Section 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.

## 2. PRINCIPLE OF CACAO-FB

To best understand the overall operation of CACAO-FB processing, we start by examining the optical scheme (see Fig. 1). Suppose we start with an unknown optical system of interest (target system). This target system consists of a lens (unknown lens) placed approximately at its focal length in front of a target sample (unknown sample). The sample is illuminated incoherently. For the sake of 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 three primary computational imaging algorithmic components that operate in sequence: 1) local aberration recovery with blur estimation; 2) full aberration recovery with a FP-based alternating projections algorithm; and 3) latent image recovery by deconvolution with coded apertures. 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. 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 $(\xi ,\eta )$, 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 ${W}_{m}(u,v)$, which collectively spans the pupil of the unknown lens, and a set of big apertures ${A}_{n}(u,v)$, which 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 ${W}_{m}(u,v)$ are labeled as ${i}_{m,t}(\xi ,\eta )$, and they encode the local aberration of the unknown lens’s pupil function in their point spread functions (PSF). The blur estimation algorithm (Section 2.B) extracts these PSFs; ${b}_{m,t}(\xi ,\eta )$. These intensity values of the spatially filtered pupil function can be synthesized into the full pupil function ${P}_{t}(u,v)$ with an FP-based alternating projections algorithm (Section 2.C). The images captured with ${A}_{n}(u,v)$, labeled ${\varphi}_{n,t}(\xi ,\eta )$, are processed with the reconstructed pupil function and the knowledge of the mask patterns to generate the latent, aberration-free image of the sample ${o}_{t}(x,y)$ (Section 2.D).

The next four sub-sections will explain the mathematical model of the image acquisition process and the three imaging algorithmic components in detail.

#### A. 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 $({x}_{0},{y}_{0})$ can be described by

where ${U}_{0}(x,y;{x}_{0},{y}_{0})$ is the complex field of the point on the sample plane, and $\delta (x-{x}_{0},y-{y}_{0})$ is the Dirac delta function describing the point located at $({x}_{0},{y}_{0})$.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 ${f}_{0}$ [Eqs. (S2) and (S3) of Supplement 1]. Any discrepancy from the ideal is incorporated into the pupil function $P(u,v;{x}_{0},{y}_{0})$, which is usually a circular bandpass filter with a uniform modulus and some phase modulation. Thus, the field right after passing through the unknown lens is

At the pupil plane, a user-defined aperture mask $M(u,v)$ is applied to produce

The CACAO-FB system captures ${i}_{m,t}(\xi ,\eta )$s and ${\varphi}_{n,t}(\xi ,\eta )$s in the data acquisition process, and these data are relayed to post-processing algorithms to recover ${o}_{t}(\xi ,\eta )$, the underlying aberration-free image of the sample. The algorithm pipeline begins with the blur estimation algorithm using ${i}_{m,t}(\xi ,\eta )$s as described below. In all the following simulations, there is one full aperture and four coded apertures ${A}_{n}(u,v)$ with the diameter of 4.5 mm; 64 small masks ${W}_{m}(u,v)$ with the diameter of 1 mm; an unknown lens L1 and L2 with the focal length of ${f}_{0}={f}_{1}={f}_{2}=100\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$; a tube lens with the focal length of ${f}_{3}=200\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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.

#### B. Local Aberration Recovery with Blur Estimation

The blur function ${b}_{m,t}(\xi ,\eta )$ associated with the small mask ${W}_{m}(u,v)$ applied to the pupil ${P}_{t}(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 ${W}_{m}(u,v)$ is set small enough such that ${W}_{m}(u,v)$ applied to a region on ${P}_{t}(u,v)$ shows a local phase map that resembles a linear phase gradient, as shown in Fig. 3(b1). In such case, the associated ${b}_{m,t}(\xi ,\eta )$ approximates a diffraction-limited spot with a spatial shift given by the phase gradient. ${W}_{m}(u,v)$ applied to other regions on ${P}_{t}(u,v)$ may have ${b}_{m,t}(\xi ,\eta )$, whose shape deviates from a spot if the masked region contains more severe aberrations, as shown in Figs. 3(b2) and 3(b3). In general, the aberration at or near the center of an imaging lens is minimal, and it becomes severe near the edge of the aperture because the lens’s design poorly approximates the parabolic shape away from the optical axis [57]. Thus, the image captured with the center mask ${i}_{1,t}(\xi ,\eta )$ is mostly aberration free with its PSF defined by the diffraction-limited spot associated with the mask’s aperture size. Other ${i}_{m,t}(\xi ,\eta )$s have the same frequency band limit as ${i}_{1,t}(\xi ,\eta )$ but are under the influence of additional aberration encapsulated by their local PSFs ${b}_{m,t}(\xi ,\eta )$s.

We adopt an image-pair-based blur estimation algorithm widely used in computational photography discipline to determine ${b}_{m,t}(\xi ,\eta )$. In this algorithm, one of the image pairs is assumed to be blur free while the other is blurred [58,59]. The blur kernel can be estimated by an iterative PSF estimation method, which is iterative Tikhonov deconvolution [60] in the Fourier domain, adopting the update scheme in Yuan’s blur estimation algorithm [58] and adjusting the step size to be proportional to $|{I}_{1,t}(u,v)|/{|{I}_{1,t}(u,v)|}_{\mathrm{max}}$ for robustness to noise [61], where ${I}_{1,t}(u,v)$ is the Fourier spectrum of ${i}_{1,t}(\xi ,\eta )$. The blur estimation process is described in Algorithm 1 as shown in Fig. 4.

The recovered ${b}_{m,t}(\xi ,\eta )$s are the intensity information of the different masked pupil regions’ Fourier transforms. They can be synthesized into the full pupil function ${P}_{t}(u,v)$ using an FP-based alternating projections algorithm, as described in the following section.

#### C. Full Aberration Recovery with Fourier-Ptychography-based Alternating Projections Algorithm

FP uses the 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 [9,10]. In our implementation, the full pupil’s complex field ${P}_{t}(u,v)$ is the desired Fourier spectrum to be synthesized, and the local PSFs ${b}_{m,t}(\xi ,\eta )$s are the aperture-scanned intensity images to be used for FP-based alternating projections, as shown in the bottom half of Fig. 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.

The FP-based alternating projections algorithm requires that the scanned apertures during image acquisition have at least 30% overlap [62] for successful phase retrieval. Thus, the updating ${b}_{m,t}(\xi ,\eta )$s in Algorithm 2 are ordered in a spiral-out pattern, each having an associated aperture ${W}_{m}(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. S1 of Supplement 1. For the simulated pupil diameter of 4.5 mm, there are 64 ${W}_{m}(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. 6. Algorithms 1 and 2 are able to estimate the local PSFs from the 64 images captured with the small masks ${W}_{m}(u,v)$ [Fig. 6(c))], and they reconstruct the complex pupil function ${P}_{t}(u,v)$ successfully [Fig. 6(e)]. A simple Fourier transformation of ${P}_{t}(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 min for Algorithm 1 and 20 s for Algorithm 2 to operate on the 64 images (1000 by 1000 pixels) taken with ${W}_{m}(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.* [63,64]. They both operate on a single blurred image [Fig. 6(d)] to simultaneously extract the blur function and the latent image [Fig. 6(f)]. For our purpose, we compare the reconstructed PSFs to gauge the performance. As shown in Fig. 6(f), the reconstructed blur kernels by MATLAB and Fergus *et al.* both show poor fidelity to the 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.

#### D. Latent Image Recovery by Deconvolution with Coded Apertures

With the knowledge of the pupil function obtained from Algorithms 1 and 2, it is possible to recover ${o}_{t}(x,y)$ from the aberrated image ${\varphi}_{t}(\xi ,\eta )$ taken with the full pupil aperture. In the Fourier domain, the image’s spectrum is represented as ${\mathrm{\varphi}}_{t}(u,v)={H}_{t}(u,v){O}_{t}(u,v)$, where ${H}_{t}(u,v)$ and ${O}_{t}(u,v)$ are the spatial spectrum of ${h}_{t}(\xi ,\eta )$ and ${o}_{t}(x,y)$, respectively. ${H}_{t}(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 ${P}_{t}(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. 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 auto-correlation process. Thus, the division of ${\mathrm{\varphi}}_{t}(u,v)$ by ${H}_{t}(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 [60] or *a priori* knowledge of the sample, such as by assuming sparsity in its total variation [65,66]. 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 [67–70]. We adopt a similar method called coded aperture proposed by Zhou and Nayar [43] 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 Algorithms 1 and 2, no *a priori* knowledge is required to recover the latent image via deconvolution. A coded aperture designed by Zhou and Nayar 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 that 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. 7 since it performs well across various noise levels [43].

The pupil function in our imaging scenario does not only consist of defocus, as 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 ${A}_{n}(u,v)$ consisting of four coded apertures and a standard circular aperture at the pupil plane, as represented in Fig. 2(c). This ensures that we obtain all spatial frequency information within the NA-limited bandpass. The PSF associated with each ${A}_{n}(u,v)$ applied to ${P}_{t}(u,v)$ is easily obtained by ${h}_{n,t}(\xi ,\eta )={|{\mathcal{F}}^{-1}\{{A}_{n}(u,v){P}_{t}(u,v)\}(\xi ,\eta )|}^{2}$ and its OTF by ${H}_{n,t}(u,v)=\mathcal{F}\{{h}_{n,t}(\xi ,\eta )\}(u,v)$. With the measured full and coded aperture images ${\varphi}_{n,t}(\xi ,\eta )$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 ${o}_{t}(x,y)$, as described by Algorithm 3 in Fig. 8 and represented in Fig. 2(d).

A simulated deconvolution procedure with the coded aperture on a Siemens star pattern is shown in Fig. 7. The combined OTF of a circular aperture and the coded aperture at four 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 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 s for Algorithm 3 to process the 5 images (1000 by 1000 pixels) taken with ${A}_{n}(u,v)$ (one full aperture, four coded apertures) to generate the deconvolution result.

## 3. EXPERIMENTAL DEMONSTRATION ON ARTIFICIAL SAMPLES

#### A. Demonstration of CACAO-FB on a Crude Lens

The CACAO-FB prototype system’s setup is simple, as shown in Fig. 9. It consists of a pair of 2 inch, $f=100\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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 [10]. 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 $+6\mathrm{D}$ trial lens (26 mm diameter, $f=130\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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 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 $\mathrm{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 $+6\mathrm{D}$ 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 [71]. 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 $\lambda /\mathrm{NA}=24.6\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$ (coherent) and $\lambda /(2\mathrm{NA})=12.3\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$ (incoherent) periodicity, defined by the spatial frequency cutoff in the coherent and incoherent transfer functions, respectively. As shown in Fig. 10, CACAO-FB can resolve features up to 19.6 μm periodicity, which is within the expected resolution limit.

The $+6\mathrm{D}$ 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 three 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. 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 up to Group 5 Element 5 (19.68 μm periodicity) are resolved after deconvolution as shown in Fig. 11, which matches closely with the resolution determined by the Siemens star pattern.

#### B. 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. 12. We make adjustments to our CACAO-FB system as shown in Fig. 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 m long, 600 μm core diameter multimode fiber (FP600URT from Thorlabs), following the method in Ref. [72]. 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 [73]). 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. 14, the recovered images show severe spatially varying aberration 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 a 4.5 mm diameter, and its associated aberrated image had a 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. 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 Shack–Hartmann wavefront sensing method, which only identifies the centroid of each blur kernel to estimate the aberration. Thus, this demonstrates the importance of the blur kernel estimation step in our algorithm and the distinct difference of our aberration reconstruction from other wavefront sensing methods.

## 4. ADAPTING CACAO-FB TO AN *IN VIVO* EXPERIMENT ON THE EYE OF A RHESUS MACAQUE

The same setup as in Section 3.B 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. Two 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 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\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mW}/\mathrm{cm}}^{2}$ on the retina in accordance with ANSI recommended safe light dosage.

Due to the safety limitation on the illumination power, the captured images of the retina have low SNR (e.g., the full-aperture image has $\mathrm{SNR}=7.5$). We increase the SNR by capturing multiple redundant frames [213 frames for ${A}_{n}(u,v)$s, 4 frames for ${W}_{m}(u,v)$s] and adding them together (Fig. S3 of Supplement 1). Thus, a long sequence of images ($\sim 45\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$) 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 [7]. Motion registration includes both translation and rotation, and these operations need 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 [74], 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. S2 of Supplement 1 and Visualization 1). A center region with half the dimensions of the full frame is selected from one of the frames as a template for registration. The normalized cross-correlation (NCC) value is found between the template and each frame [75] for every rotation angle ($-1.5$ to 1.5 deg, 0.0015 deg 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. S3 of Supplement 1).

The deconvolution result is shown in Fig. 16. Although the sensor size of $2560\times 2160$ pixels is used for capturing raw images, the resultant averaged images are only $1262\times 1614$ pixels after the motion registration. The input full-aperture image had $\mathrm{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 due to the lower SNR in the darker regions leading to a 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 the CACAO-FB algorithm pipeline.

## 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, which 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 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 fewer 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. 3. Therefore, reconstructing the blur kernels of the limited aperture images is critical for CACAO-FB’s performance, as is demonstrated in Fig. 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. [70] to remove nulls in the OTF of an aberrated imaging system.

Although the recovered retinal image in Section 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 *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. [7]. 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 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 the 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.

## Funding

National Institutes of Health (NIH) (R21 EY026228A).

## Acknowledgment

We thank Amir Hariri for assisting with the experiment; Soo-Young Kim, Dierck Hillmann, Martha Neuringer, Laurie Renner, Michael Andrews and Mark Pennesi for being of tremendous help over email regarding the experimental setup and *in vivo* eye imaging; and Mooseok Jang, Haowen Ruan, Edward Haojiang Zhou, Joshua Brake, Michelle Cua, Hangwen Lu, and Yujia Huang for helpful discussions.

See Supplement 1 for supporting content.

## REFERENCES

**1. **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**, 470–473 (1996). [CrossRef]

**2. **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). [CrossRef]

**3. **P. Godara, A. Dubis, A. Roorda, J. Duncan, and J. Carroll, “Adaptive optics retinal imaging: Emerging clinical applications,” Optom. Vis. Sci. **87**, 930–941 (2010). [CrossRef]

**4. **D. Williams, “Imaging single cells in the living retina,” Vis. Res. **51**, 1379–1396 (2011). [CrossRef]

**5. **D. L. Fried, “Anisoplanatism in adaptive optics,” J. Opt. Soc. Am. **72**, 52–61 (1982). [CrossRef]

**6. **M. J. Booth, “Adaptive optical microscopy: the ongoing quest for a perfect image,” Light Sci. Appl. **3**, e165 (2014). [CrossRef]

**7. **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,” Opt. Express **8**, 631–643 (2001). [CrossRef]

**8. **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). [CrossRef]

**9. **G. Zheng, R. Horstmeyer, and C. Yang, “Wide-field, high-resolution Fourier ptychographic microscopy,” Nat. Photonics **7**, 739–745 (2013). [CrossRef]

**10. **R. Horstmeyer, X. Ou, J. Chung, G. Zheng, and C. Yang, “Overlapped Fourier coding for optical aberration removal,” Opt. Express **22**, 24062–24080 (2014). [CrossRef]

**11. **X. Ou, G. Zheng, and C. Yang, “Embedded pupil function recovery for Fourier ptychographic microscopy,” Opt. Express **22**, 4960–4972 (2014). [CrossRef]

**12. **X. Ou, R. Horstmeyer, C. Yang, and G. Zheng, “Quantitative phase imaging via Fourier ptychographic microscopy,” Opt. Lett. **38**, 4845–4848 (2013). [CrossRef]

**13. **Z. Bian, S. Dong, and G. Zheng, “Adaptive system correction for robust Fourier ptychographic imaging,” Opt. Express **21**, 32400–32410 (2013). [CrossRef]

**14. **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,” Sci. Rep. **6**, 27384 (2016). [CrossRef]

**15. **L. Bian, J. Suo, G. Zheng, K. Guo, F. Chen, and Q. Dai, “Fourier ptychographic reconstruction using Wirtinger flow optimization,” Opt. Express **23**, 4856–4866 (2015). [CrossRef]

**16. **J. M. Rodenburg and R. H. T. Bates, “The theory of super-resolution electron microscopy via Wigner-distribution deconvolution,” Philos. Trans. R. Soc. A **339**, 521–553 (1992). [CrossRef]

**17. **H. M. L. Faulkner and J. M. Rodenburg, “Movable aperture lensless transmission microscopy: a novel phase retrieval algorithm,” Phys. Rev. Lett. **93**, 023903 (2004). [CrossRef]

**18. **A. Pan and B. Yao, “Three-dimensional space optimization for near-field ptychography,” Opt. Express **27**, 5433–5446 (2019). [CrossRef]

**19. **M. G. L. Gustafsson, “Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy,” J. Microsc. **198**, 82–87 (2000). [CrossRef]

**20. **M. G. L. Gustafsson, “Nonlinear structured-illumination microscopy: wide-field fluorescence imaging with theoretically unlimited resolution,” Proc. Natl. Acad. Sci. USA **102**, 13081–13086 (2005). [CrossRef]

**21. **J. Qian, S. Dang, Z. Wang, X. Zhou, D. Dan, B. Yao, Y. Tong, H. Yang, Y. Lu, Y. Chen, X. Yang, M. Bai, and M. Lei, “Large-scale 3D imaging of insects with natural color,” Opt. Express **27**, 4845–4857 (2019). [CrossRef]

**22. **T. M. Turpin, L. H. Gesell, J. Lapides, and C. H. Price, “Theory of the synthetic aperture microscope,” Proc. SPIE **2566**, 230–240 (1995). [CrossRef]

**23. **J. Di, J. Zhao, H. Jiang, P. Zhang, Q. Fan, and W. Sun, “High resolution digital holographic microscopy with a wide field of view based on a synthetic aperture technique and use of linear CCD scanning,” Appl. Opt. **47**, 5654–5659 (2008). [CrossRef]

**24. **T. R. Hillman, T. Gutzler, S. A. Alexandrov, and D. D. Sampson, “High-resolution, wide-field object reconstruction with synthetic aperture Fourier holographic optical microscopy,” Opt. Express **17**, 7873–7892 (2009). [CrossRef]

**25. **L. H. Yeh, J. Dong, J. Zhong, L. Tian, and M. Chen, “Experimental robustness of Fourier ptychography phase retrieval algorithms,” Opt. Express **23**, 33214–33240 (2015). [CrossRef]

**26. **L. Bian, J. Suo, G. Situ, G. Zheng, F. Chen, and Q. Dai, “Content adaptive illumination for Fourier ptychography,” Opt. Lett. **39**, 6648–6651 (2014). [CrossRef]

**27. **J. Sun, Q. Chen, Y. Zhang, and C. Zuo, “Efficient positional misalignment correction method for Fourier ptychographic microscopy,” Biomed. Opt. Express **7**, 1336–1350 (2016). [CrossRef]

**28. **Y. Zhang, W. Jiang, L. Tian, L. Waller, and Q. Dai, “Self-learning based Fourier ptychographic microscopy,” Opt. Express **23**, 18471–18486 (2015). [CrossRef]

**29. **P. Li, D. J. Batey, T. B. Edo, and J. M. Rodenburg, “Separation of three-dimensional scattering effects in tilt-series Fourier ptychography,” Ultramicroscopy **158**, 1–7 (2015). [CrossRef]

**30. **R. Horstmeyer, J. Chung, X. Ou, G. Zheng, and C. Yang, “Diffraction tomography with Fourier ptychography,” Optica **3**, 827–835 (2016). [CrossRef]

**31. **A. Pan, Y. Zhang, K. Wen, M. Zhou, J. Min, M. Lei, and B. Yao, “Subwavelength resolution Fourier ptychography with hemispherical digital condensers,” Opt. Express **26**, 23119–23131 (2018). [CrossRef]

**32. **A. Pan, Y. Zhang, T. Zhao, Z. Wang, D. Dan, M. Lei, and B. Yao, “System calibration method for Fourier ptychographic microscopy,” J. Biomed. Opt. **22**, 096005 (2017). [CrossRef]

**33. **T. Kamal, L. Yang, and W. M. Lee, “In situ retrieval and correction of aberrations in moldless lenses using Fourier ptychography,” Opt. Express **26**, 2708–2719 (2018). [CrossRef]

**34. **J. Chung, J. Kim, X. Ou, R. Horstmeyer, and C. Yang, “Wide field-of-view fluorescence image deconvolution with aberration-estimation from Fourier ptychography,” Biomed. Opt. Express **7**, 352–368 (2016). [CrossRef]

**35. **L. Tian, X. Li, K. Ramchandran, and L. Waller, “Multiplexed coded illumination for Fourier Ptychography with an LED array microscope,” Biomed. Opt. Express **5**, 2376–2389 (2014). [CrossRef]

**36. **J. Chung, H. Lu, X. Ou, H. Zhou, and C. Yang, “Wide-field Fourier ptychographic microscopy using laser illumination source,” Biomed. Opt. Express **7**, 4787–4802 (2016). [CrossRef]

**37. **L. Tian and L. Waller, “3D intensity and phase imaging from light field measurements in an LED array microscope,” Optica **2**, 104–111 (2015). [CrossRef]

**38. **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**, 904–911 (2015). [CrossRef]

**39. **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**, 066007 (2014). [CrossRef]

**40. **S. Dong, K. Guo, P. Nanda, R. Shiradkar, and G. Zheng, “FPscope: a field-portable high-resolution microscope using a cellphone lens,” Opt. Express **5**, 3305–3310 (2014). [CrossRef]

**41. **J. Sun, C. Zuo, L. Zhang, and Q. Chen, “Resolution-enhanced Fourier ptychographic microscopy based on high-numerical-aperture illuminations,” Sci. Rep. **7**, 1187 (2017). [CrossRef]

**42. **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,” Opt. Express **23**, 26999–27010 (2015). [CrossRef]

**43. **C. Zhou and S. Nayar, “What are good apertures for defocus deblurring?” in *IEEE International Conference on Computational Photography* (IEEE, 2009), pp. 1–8.

**44. **J. R. Fienup and J. J. Miller, “Aberration correction by maximizing generalized sharpness metrics,” J. Opt. Soc. Am. A **20**, 609–620 (2003). [CrossRef]

**45. **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,” Sci. Rep. **6**, 35209 (2016). [CrossRef]

**46. **F. Soulez, L. Denis, Y. Tourneur, and E. Thiébaut, “Blind deconvolution of 3D data in wide field fluorescence microscopy,” in *9th IEEE International Symposium on Biomedical Imaging (ISBI)* (IEEE, 2012), pp. 1735–1738.

**47. **E. Thiébaut and J.-M. Conan, “Strict a priori constraints for maximum-likelihood blind deconvolution,” J. Opt. Soc. Am. A **12**, 485–492 (1995). [CrossRef]

**48. **S. Adie, B. Graf, A. Ahmad, S. Carney, and S. Boppart, “Computational adaptive optics for broadband optical interferometric tomography of biological tissue,” Proc. Natl. Acad. Sci. USA **109**, 7175–7180 (2012). [CrossRef]

**49. **N. Shemonski, F. South, Y.-Z. Liu, S. Adie, S. Carney, and S. Boppart, “Computational high-resolution optical imaging of the living human retina,” Nat. Photonics **9**, 440–443 (2015). [CrossRef]

**50. **D. Kundur and D. Hatzinakos, “Blind image deconvolution,” IEEE Signal Process. Mag. **13**(3), 43–64 (1996). [CrossRef]

**51. **A. Kumar, W. Drexler, and R. A. Leitgeb, “Subaperture correlation based digital adaptive optics for full field optical coherence tomography,” Opt. Express **21**, 10850–10866 (2013). [CrossRef]

**52. **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**, 924–931 (2017). [CrossRef]

**53. **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*,” Biomed. Opt. Express **9**, 472–485 (2018). [CrossRef]

**54. **G. Gunjala, S. Sherwin, A. Shanker, and L. Waller, “Aberration recovery by imaging a weak diffuser,” Opt. Express **26**, 21054–21068 (2018). [CrossRef]

**55. **G. Zheng, X. Ou, R. Horstmeyer, and C. Yang, “Characterization of spatially varying aberrations for wide field-of-view microscopy,” Opt. Express **21**, 15131–15143 (2013). [CrossRef]

**56. **B. K. Gunturk and X. Li, *Image Restoration: Fundamentals and Advances* (CRC Press, 2012), Vol. 7.

**57. **J. Goodman, *Introduction to Fourier Optics* (McGraw-Hill, 2008).

**58. **L. Yuan, J. Sun, L. Quan, and H.-Y. Shum, “Image deblurring with blurred/noisy image pairs,” ACM Trans. Graph. **26**, 1 (2007). [CrossRef]

**59. **S. H. Lim and D. A. Silverstein, “Method for deblurring an image,” U.S. patent 8,654,201 (February 18, 2014).

**60. **A. Neumaier, “Solving ill-conditioned and singular linear systems: a tutorial on regularization,” SIAM Rev. **40**, 636–666 (1998). [CrossRef]

**61. **J. M. Rodenburg and H. M. L. Faulkner, “A phase retrieval algorithm for shifting illumination,” Appl. Phys. Lett. **85**, 4795–4797 (2004). [CrossRef]

**62. **J. Sun, Q. Chen, Y. Zhang, and C. Zuo, “Sampling criteria for Fourier ptychographic microscopy in object space and frequency space,” Opt. Express **24**, 15765–15781 (2016). [CrossRef]

**63. **R. Fergus, B. Singh, A. Hertzmann, S. T. Roweis, and W. T. Freeman, “Removing camera shake from a single photograph,” ACM Trans. Graph. **25**, 787–794 (2006). [CrossRef]

**64. **A. Levin, Y. Weiss, F. Durand, and W. T. Freeman, “Understanding blind deconvolution algorithms,” IEEE Trans. Pattern Anal. Mach. Intell. **33**, 2354–2367 (2011). [CrossRef]

**65. **J. M. Bioucas-Dias, M. A. T. Figueiredo, and J. P. Oliveira, “Total variation-based image deconvolution: a majorization-minimization approach,” in *IEEE International Conference on Acoustics Speech and Signal Processing* (IEEE, 2006), pp. 861–864.

**66. **A. Levin, R. Fergus, F. Durand, and W. T. Freeman, “Image and depth from a conventional camera with a coded aperture,” ACM Trans. Graph. **26**, 70 (2007). [CrossRef]

**67. **E. Dowski and T. W. Cathey, “Extended depth of field through wavefront coding,” Appl. Opt. **34**, 1859–1866 (1995). [CrossRef]

**68. **K. Kubala, E. R. Dowski, and W. T. Cathey, “Reducing complexity in computational imaging systems,” Opt. Express **11**, 2102–2108(2003). [CrossRef]

**69. **G. Muyo and A. R. Harvey, “Wavefront coding for athermalization of infrared imaging systems,” Proc. SPIE **5612**, 227–235 (2004). [CrossRef]

**70. **G. Muyo, A. Singh, M. Andersson, D. Huckridge, A. Wood, and A. R. Harvey, “Infrared imaging with a wavefront-coded singlet lens,” Opt. Express **17**, 21118–21123 (2009). [CrossRef]

**71. **R. Horstmeyer, R. Heintzmann, G. Popescu, L. Waller, and C. Yang, “Standardizing the resolution claims for coherent microscopy,” Nat. Photonics **10**, 68–71 (2016). [CrossRef]

**72. **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,” Opt. Express **14**, 4552–4569 (2006). [CrossRef]

**73. **A. Martins da Silva and B. Leal, “Photosensitivity and epilepsy: current concepts and perspectives-a narrative review,” Seizure **50**, 209–218 (2017). [CrossRef]

**74. **L. Yaroslavsky, *Theoretical Foundations of Digital Imaging Using MATLAB* (CRC Press, 2013).

**75. **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,” Opt. Express **3**, 190–197 (1998). [CrossRef]