Single-photon avalanche diodes (SPADs) are highly sensitive optical detectors that can detect and count individual photons. Recently, SPAD arrays, which contain thousands to millions of individual detectors, can be manufactured at scale via current CMOS technology and are entering the larger commercial market [1]. Coupled with a pulsed light source, SPAD array cameras are often implemented in “time-correlated single-photon counting” (TCSPC) schemes that offer picosecond-scale detection resolution for photon time-of-arrival. This has led to impressive 3D imaging capabilities for autonomous vehicles [2], fluorescence lifetime imaging experiments in microscopy [3], as well as computational imaging tasks such as three-dimensional imaging and ranging [4,5], high-speed correction [6], imaging/sensing through scattering media such as fog or tissue [7–9], transient image sensing [10], and forming images of objects hidden by occlusions or around corners [11]. One current limitation of existing SPAD array cameras, as compared to traditional charged-coupled device (CCD) or CMOS camera sensors, is a limited total pixel count. While recent academic work has demonstrated SPAD arrays with up to one million pixels [12,13], commercially available SPAD arrays typically include hundreds to thousands of pixels. CCD or CMOS sensors, on the other hand, are now commonly available with tens of megapixels. If the goal is to form a high-resolution image over a large area with an SPAD array camera, then one typically has to resort to some form of mechanical scanning, either of the image sensor, imaging assembly, or the object of interest, to produce a high pixel count final image [14]. Such mechanical scanning can be slow and lengthens the imaging process.

One alternative method that we study here, termed Fourier ptychography (FP) [15,16], can improve image resolution without requiring any mechanical scanning. FP is a synthetic aperture imaging technique that acquires standard intensity-only images of an object under illumination from different angles, as shown in Figs. 1(a) and 1(b). An iterative phase retrieval algorithm [17] is then used to convert the sequence of captured low-resolution images into a high-resolution image reconstruction, while simultaneously recovering the missing phase of the scattered light. In the past, FP has primarily been implemented with an LED array for illumination [see Fig. 1(a)] to improve the space–bandwidth product [15] of microscopes by a factor of 25× or more, leading to gigapixel-sized high-resolution image reconstructions over large fields of view. FP also recovers the sample’s phase and provides the capability to correct for lens aberrations [18], which has led to low-cost yet high-performance digital microscopes [19,20], and also extends to macroscopic imaging of large objects [21,22]. Over the past several years, FP has been implemented with several types of imaging sensor, including standard CMOS [15], scientific CMOS (sCMOS) [23], intensified CCD (ICCD) cameras [24], and low-cost mobile phone sensors [20]. As SPAD arrays offer a limited pixel count and are also commonly used with rapid active illumination schemes, it appears that FP offers a useful means to improve SPAD array imaging performance in a procedure that is compatible with many use cases.

In this work, we aim to study the feasibility of performing FP using SPAD array cameras. To begin down this path, we must consider and account for a relatively unique characteristic of SPAD array cameras—a significantly limited image bit depth—which is the primary focus of this investigation. There are two general modes of SPAD array cameras, as outlined in Fig. 1(c): (1) the photon-counting mode, which simply counts photon-triggered events, and (2) the TCSPC mode, which records and digitizes time-of-arrival information about photon-triggered events. Both of these modes effectively report low dynamic range or binary spatial intensity information.

There has been prior work considering the dynamic range of images in such computational imaging setups. For example, work within compressive sensing [25–27] has examined the dynamic range problem at the extreme of binary image measurements. In addition, several recent works have considered the general problem of phase retrieval from binary measurements in mathematical context [28–30]. Likewise, work has explored phase retrieval from measured intensities of speckle quantized to 1-bit [31]. However, none of these prior studies consider synthetic aperture imaging-type scenarios such as FP, which uses multiple low-resolution binary images to generate a reconstructed image that simultaneously has a higher resolution and higher dynamic range. Prior work with FP has indirectly accounted for dynamic range challenges with a number of strategies, for example, by varying the image exposure time for each LED [32], ensuring high-dynamic-range capture with an sCMOS camera [23], or employing an high dynamic range (HDR) image capture method with multiple exposures per LED [15]. However, without a priori knowledge of the specimen’s scattering properties, it is generally challenging to select an appropriate set of exposure times for all acquired images (see Section 2.B).

Here, we present a new “quantized” FP image reconstruction procedure that maps low-dynamic-range image data to high-resolution complex-valued image reconstructions. By accounting for each acquired image’s bit depth and exposure duration, we verify both in simulation and experiment that FP can effectively produce a high-dynamic-range reconstruction at expected resolution from raw binary intensity images. We also experimentally demonstrate FP using SPAD array image data, pointing the way toward implementation with novel 3D ranging and time-gating imaging applications in the future. However, our findings here may also be used to improve any general FP imaging setup on standard image sensors as well, as sensor dynamic range is often an important factor that has received limited attention to date.

In FP, we assume that we are imaging a thin sample represented by a complex function s(r→), with amplitude describing the sample’s absorption, and phase describing its imparted phase shift due to the sample’s thickness variation. We treat the illumination generated by each LED at the sample plane as a quasi-monochromatic plane wave of uniform intensity that travels at a particular angle θn: eiϕ(r→,θn), where n denotes light from the n^th LED. Since we have a thin sample, the resultant field emerging from the sample surface is given by the elementwise product s(r→)·eiϕ(r→,θn). This field then propagates through the lens pupil plane (i.e., back focal plane) and then onto the image sensor. Following principles of Fourier optics, we can represent this field directly before the pupil plane as the Fourier transform of the above product: O(k→n)=F[s(r→)·eiϕ(r→,θn)]. Assuming the pupil has a transmission function P(k→), the field immediately after the pupil plane is given by the product of P(k→) and O(k→n). Finally, by modeling propagation to the image sensor with an inverse Fourier transform, we can denote the entire expression for formation of the nth image as zn(r→)=|F−1{F[s(r→).eiϕ(r→,θn)]·P(k→)}|2,(1)where the modulus square is due to the detector’s ability to only detect incident intensity. According to the Fourier shift theorem, the tilted plane wave illumination eiϕ(r→,θn) results in a shift in the Fourier transform of the sample (kn=2πλsinθn, where λ is the wavelength), which allows us to express the intensity observed at the image plane for FP as zn(r→)=|F−1[O(k→−kn)·P(k→)]|2,(2)where O(k→)=F[s(r→)]. This equation represents the standard FP forward model with zn(r→) in units of W/m2 as irradiance immediately before the image detector, which we also assume is spatially discretized to the sensor pixel spacing. A typical digital image is formed via an exposure of some finite duration Δt across a finite pixel area, which converts the irradiance to radiant energy that is detected as an electronic signal. This signal is typically quantized to generate a digital image. To account for the dynamic range of the digital images in FP, we write the quantization process introduced by a general image sensor array as yn(r→)=QM[β·zn(r→)Δtn].(3)

Here, the product of zn and the exposure Δt follows from the reciprocity property [33] and β is a constant scalar that accounts for factors involved in the mapping of radiant energy per pixel to an electronic signal (pixel fill factor, quantum efficiency, etc.). In the above two equations, we assume that both shot noise and sensor read noise are negligible so that the two scenarios (standard FP and quantized FP models) can be easily compared. QM(x) is a quantization function that maps a continuous quantity to a discrete set of 2M positive integer values within the range [0,2M−1], which corresponds to M bits of information: QM(x)={⌊x·2M⌋ifx<12M−1ifx≥1,(4)where ⌊x⌋ is a flooring function that gives the largest integer less than or equal to x. In other words, the quantization function QM(x) outputs a non-negative integer that is a function of the detected radiant energy with a maximum value 2M. The second line of Eq. (4) typically represents the scenario of image saturation. Finally, adopting a “relative” per-image exposure time, αn=Δtnβ, leads to the following simple expression for the nth detected image: yn(r→)=QM[αn·zn(r→)].(5)

The above forward model provides a straightforward expression for discretized images captured on limited dynamic range sensors. However, it assumes a priori knowledge of an adaptive exposure time, αn, which is typically manipulated to ensure that quantized images formed under illumination from the nth LED are well exposed. For example, With the particular selection of αn=1/max(zn), we can see that the argument of QM cannot exceed 2M. Accordingly, all captured images will be exposed up to (but will not exceed) the maximum dynamic range of the image sensor to produce images without any saturation.

There are three general issues surrounding the experimental specification of αn. First, the maximum intensity reaching the image sensor for each LED, max(zn), is typically an unknown quantity and must be estimated. Second, this quantity will clearly vary as a function of LED illumination angle, as LEDs toward the optical axis will often generate brighter images than those illuminating the sample from higher angles. And third, in many imaging scenarios, image saturation can actually be beneficial. An extreme example of this is in the binary image detection case (e.g.,  with an SPAD array). In this extreme case, from an information collection perspective, it should typically be desirable to select a relative exposure time to ensure that approximately half of the pixels output a quantized measurement of 1, while the remaining pixels output 0. Satisfying this condition will typically require overexposure (i.e., saturation) of a large fraction of pixels. To explore the above issues, we consider four different methods of selecting the relative exposure time αn for each captured image (see Fig. 2).•Uniform: without any a priori knowledge of the specimen, one might ensure that the center-illuminated image (i.e., n=1 LED) is well exposed, and then apply this exposure time to all other captured images such that αn=1/max(z1). We refer to this as the uniform exposure strategy, where the same exposure is used for images taken from different LEDs.•Max: to significantly increase the information preserved through detection and quantization, we can instead apply a strategy that scales the image acquired by each LED by its maximum intensity, such that αn=1/max(zn). Here, each image will be evenly quantized across M bits without saturation. This strategy requires knowledge of the maximum intensity at the image plane produced by each LED for the sample of interest, which is typically not possible in practice.•Informed: to address the experimental limitations of the max method, we propose a third exposure time selection strategy that is practically achievable. We assume that the object of interest can be associated within a general category S of other similar objects (e.g.,  with thin tissue sections or blood smears in microscopic imaging applications) and that we have access to the maximum image intensity per optical source max(zn,m) generated by the other m objects within this category (e.g.,  acquired during prior imaging experiments or in simulation). The relative exposure time per LED is then set at αn,S=1|S|∑m∈S1max(zn,m),(6)where |S| is the cardinality of the image set S. The exposure time is selected as the average per-LED image intensity generated by a set of similar specimens S that have been examined in the past.•Median: one optimal quantization strategy, from an information collection perspective, is to maximize the entropy of each detected image, by making each detected histogram as uniform as possible across all possible quantized values. For the M-bit discretization case, it is generally typically not possible to achieve histogram uniformity without the ability to create a quantization function with non-uniform bin widths, which is not possible with current image sensors. For the binary quantization case, however, it is possible to select a per-image exposure time that generates an even histogram (i.e., pixels of the same number are quantized to 0 as to 1). This is achieved when the normalized exposure time is αn=12×median(zn).(7)

This “median” exposure time method could be implemented with a variety of binary imaging schemes that utilize a separate photodetector to trigger the exposure on/off period, for example.

As shown in Fig. 2, each strategy above will expand or contract the histogram of measured intensities in a unique way. In this example, it can be seen that max and informed methods have a similar effect. The median method shows a significant improvement in relaying specimen information over the max method for binary FP imaging. While we selected to examine these four methods due to their conceptual simplicity and efficiency, more sophisticated exposure strategies are certainly possible and should be explored in future work. We will refer to these four strategies at different points during our simulations and experiments presented below.

In stage 1, the Fourier spectrum estimate from the previous t−1 iteration, Ot−1,n(k→), is shifted according to the illumination angle of the n^th LED as Ot−1,n(k→−kn). It is then low-pass filtered by the pupil function P(k→) and Fourier transformed to generate an estimate of the complex field at the image plane ψ(r→). The intensity zn(r→) of this complex field can then be calculated as |ψ(r→)|2, which corresponds to the forward model shown in Eq. (2) and is an estimate of the high-dynamic-range ground truth image z^n(r→) (at the full precision of the reconstruction software). However, this dynamic range rarely matches that of the image sensor used in experiments, where a quantized low-dynamic-range image y^n(r→), of the ground truth image intensity z^n(r→), is recorded instead. To account for this discrepancy, our proposed algorithm instead computes a quantized image yn(r→) that is produced by the assumed detector, as yn(r→)=QM[αn·|ψ(r→)|2], which corresponds to the forward model in Eq. (5). We hypothesize that this quantized image estimate yn(r→) provides a useful estimate of the experimentally captured and quantized image, y^n(r→). The goal of the quantized FP algorithm is to have yn(r→) match y^n(r→) upon its computation.

In stage 2, the amplitude of the estimated complex field at the image plane, ψ(r→), is updated using experimentally captured data. In standard FP, ψ(r→) is updated by directly replacing the estimated amplitudes with measured amplitudes via the following equation: ψ′(r→)=y^n(r→)exp[i∠ψ(r→)].(8)

Here, ψ′(r→) is the updated complex field in the image plane and ∠ψ(r→) is the phase of ψ(r→). Note that Eq. (8) does not consider image quantization—that is, the measurement y^n(r→) is directly employed for reconstruction update and is not compared to the image estimate in quantized form. In our proposed model, we modify this update step to ensure effective and accurate algorithm convergence with, for example, binary SPAD array data. Specifically, we account for the discrepancy between the current algorithm’s quantized image estimate yn(r→) and the experimentally captured and quantized image y^n(r→) to promote accurate image update and account for update inconsistencies [34]. We introduce an indicator called error map E between actual and estimated discretized images as follows: E={0y^n(r→)=yn(r→)1y^n(r→)≠yn(r→).(9)

If E is zero at a particular pixel, the corresponding estimate |ψ(r→)|’s associated detected and discretized intensity value at that pixel has accurately converged to the correct quantized value. Hence, we assume that particular value for ψ(r→) no longer requires an intensity update. Otherwise, an update is still necessary. To incorporate this observation into ptychography’s iterative phase retrieval cycle, quantized FP replaces Eq. (8) with the following update function: ψ′(r→)=[y^n(r→)αnE+|ψ(r)|(1−E)]exp[i∠ψ(r→)].(10)

In the binary case, this modification ensures that at every iteration the estimated image’s amplitude is either updated with the appropriate value or remains the estimated value, with the latter case being when the estimate’s quantized pixel intensity value agrees with experimental measurement. In stage 3, ψ′(r→) and ψ(r→) are used to update the Fourier spectrum estimate of the sample as follows: Ot,n(k→−k→n)=Ot−1,n(k→−k→n)−δ[P†(k→)]F[ψ′(r→)−ψ(r→)]|P(k→)|max[|P(k→)|2+ξ],(11)where δ is the step size, ξ is a regularization parameter, and P†(k→) is the complex conjugate of P(k→). This updated Fourier spectrum is shifted back according to the LED illumination angle to obtain Ot,n(k→). Once all T iterations are complete, O is inverse Fourier transformed to produce the high-resolution complex image reconstruction. The update function in the above equation is based on the sequential Gauss–Newton method, widely used in FP [17]. Many other similar Fourier spectrum update functions [35] can be used with the proposed quantized error map strategy. Embedded pupil recovery [18], for example, can also be used in our model if pupil aberrations need to be recovered. While this algorithm was developed for the particular case of performing FP with binary imagery, it also has improved performance when using non-binary low bit depth detection schemes as well, as shown in our simulations.

To validate our proposed quantized FP reconstruction algorithm and the adaptive exposure strategy for low-dynamic-range and binary image datasets, we first performed a series of simulations. First, we simulated imaging of standard resolution test targets and natural samples with both 8-bit and 1-bit image sensors to demonstrate the resolution improvement and reconstruction quality. Second, we use several open-source pathology images as samples to compare the average performance of our proposed reconstruction method to standard methods, and to examine the performance of each proposed exposure strategy. Finally, we study the influence of the overlap ratio of Fourier sub-spectra on reconstruction quality.

Simulations assumed an imaging NA of 0.1 with 1.8× magnification and 2.4 μm detector pixel size. We assumed a 0.205 illumination NA provided by a 15×15 LED illumination array (225 low-resolution images captured in total). These parameters are used for all the simulations except for the overlap influence simulations, where the illumination NA is kept constant, but the total number of LEDs and their pitch are varied to provide different overlap values. Poisson noise was added to the low-resolution images generated in the simulation before discretization. To evaluate quality of the complex reconstruction results, we used the following normalized mean square error [36] metric: NMSE=∑r→|S(r→)−γR(r→)|2∑r→|R(r→)|2,γ=∑r→S(r→)R(r→)∑r→|R(r→)|2.(12)

Here, S(r→) is the normalized reference complex image, R(r→) is the normalized reconstructed complex image, and γ is a factor to account for arbitrary constant reconstruction phase offset [36].

Simulation results in Fig. 4 demonstrate the resolution improvement and image quality for cases using 1-bit images captured with the median exposure strategy and 8-bit images using the max exposure strategy. An amplitude-only USAF target and a phase-only Siemens star target were used to verify resolution improvement, where the expected resolution cutoff is 2.04 μm per line pair (lp). In this simulation, group 1 element 1 in the USAF target corresponds to approximately 2 μm per line pair, while group 1 element 2 is 1.78 μm per line pair. We can see that both standard FP and our proposed quantized FP reconstruction using binary data can resolve the former element but not the latter, achieving maximum expected resolution performance. The contrast ratios between binary and 8-bit in Fig. 4(a), columns (1) and (2) are 0.998 and 0.995. A contrast-varying resolution target [Fig. 4(a3)] highlights how the complex reconstruction achieved with an 8-bit image dataset maintains higher resolution in areas of low object contrast (yellow line refers to 19% contrast compared to the maximum sinusoid contrast at right), as compared to the 1-bit image dataset reconstruction. However, the 1-bit dataset reconstruction still offers full resolution improvement to the expected cutoff at slightly higher contrast (34% orange line). Finally, natural sample reconstruction quality looks almost identical in both the 8-bit and binary imaging cases, with similar normalized mean square error (NMSE) values (0.010 and 0.042, respectively). These results first demonstrate that FP is possible with binary measurements, despite significantly altered and deteriorated captured image data [see Fig. 4 column (4)]. Somewhat surprisingly, the quality of high-resolution reconstruction using quantized FP appears comparable to the standard FP case with 8-bit images.

A comparison between the standard FP reconstruction algorithm [Eq. (8)] and the proposed quantized FP algorithm [Eq. (10)] is performed in Fig. 5(a) assuming the use of the max exposure method. We varied the bit depth of captured low-resolution images from 1 to 8 and evaluated the NMSE of the resulting reconstructions from both algorithms (averaged from 30 similar pathology samples). In this plot, quantized FP shows a significant improvement in low-dynamic-range imaging cases, with a small improvement for higher-dynamic-range imaging cases. This trend can also be observed in the sample reconstructions shown in Fig. 5(d). Since the main focus of this work is on FP reconstruction for 1-bit SPAD array data, these results highlight the clear advantage of the proposed algorithm. Reconstruction quality using the four exposure strategies of interest is also compared using quantized FP in Fig. 5(b). As expected, the informed and max strategies trend closely, whereas the uniform strategy falls significantly behind. Also, the median case, which we assume can only be experimentally implemented for 1-bit imaging scenarios currently, shows a dramatic performance improvement. This suggests that the median strategy offers a promising exposure or source brightness selection method for FP implementations with SPAD array cameras.

As mentioned in the previous section, FP algorithm performance for high-resolution reconstruction depends upon various parameters. A crucial parameter is the redundancy in captured data. As a phase retrieval method, FP must sample more pixels than that it reconstructs. While increasing data redundancy is not usually ideal (e.g.,  it increases data capture time and memory requirements), a minimum redundancy criterion should be met for ptychographic image reconstruction strategies [37]. To investigate the impact of this parameter on quantized FP, we varied the number of captured images by varying the number reconstruction pair is shownof assumed illumination sources provided in the illumination array. This allowed us to keep the illumination numerical aperture constant but vary the amount of image dataset redundancy. For each illumination configuration, we then computed FP reconstruction error (NMSE) for a number of different assumed image sensor bit depths. We summarize data redundancy by computing the overlap percentage between adjacent captured areas which are shifted by 2πsinθλ in the Fourier spectrum domain, where θ is the angular shift between adjacent illumination sources. The results of this analysis are plotted in Fig. 6, where each point is an average of NMSE reconstruction across a set of 12 simulated complex samples (from open-source Cell Image Library [38]). Here, it can be seen that the low-bit images (1–2 bits) require much higher overlap and thus greater data redundancy for accurate reconstruction. All bit depths higher than 4 exhibit relatively similar performance. This is partially attributed to our quantized FP algorithm and our use of the proposed max exposure strategy. Reconstruction quality for 1-bit image data improves with larger overlap percentages, but this increase appears insignificant for detected bit depths greater than 3 if >40% overlap is used.

Before testing quantized FP with an SPAD array, we first performed a series of preliminary experiments with a standard CMOS camera (Basler 12-bit CMOS sensor, acA4024-29um) to capture FP datasets, whose image data was manipulated post-capture to exhibit a variety of different bits per pixel. The FP imaging system in these experiments used a 0.1 NA objective lens at 2× magnification with a 0.42 illumination NA. We used 15×15 red LEDs from the Adafruit32×32 RGB LED matrix with a 4 mm separation between adjacent LEDs, placed 60 mm beneath the sample plane. All bright-field images were captured with 40 ms exposure and dark-field images with 200 ms exposure. We normalized each image in the captured dataset with respect to its exposure time to generate a floating point matrix in MATLAB. We then binarized this dataset using the median exposure strategy to create a 1-bit image dataset. This low-resolution binary image dataset was then input into our quantized FP reconstruction algorithm to create a high-resolution complex image estimate. The results (Dataset 1, Ref. [39] and Code 1, Ref. [40]) of this exercise are shown in Fig. 7 for both a resolution test target and a standard biological sample (a peripheral blood smear). The imaging system’s standard resolution is approximately 6 μm/lp. After FP reconstruction, the expected resolution cutoff is 1.2 μm/lp. In this experiment, we were able to demonstrate 1.40 μm/lp resolution using the 12-bit raw image data, and 2.5 μm/lp using the binary image data and quantized FP reconstruction. The reconstruction quality using binary FP images closely matches the standard FP reconstruction algorithm output with 12-bit images, albeit with small added artifacts and a slight decrease in overall resolution, given that 12× less data was utilized. More advanced reconstruction update procedures that better account for noise and experimental calibration errors may improve upon these initial results.

To demonstrate quantized FP with an SPAD array camera, we utilized the PF32 SPAD array camera from Photon Force Ltd., which has 32×32 pixels with a pixel pitch of 50 μm and 1.5% fill factor. A 0.1 NA Olympus plan achromatic objective with a 600 mm focal length tube lens was used for imaging and to provide 18.5× magnification for Nyquist sampling. The field of view (FOV) of the SPAD array camera is extremely small due to the limited number of pixels available, so a regular camera (Basler acA4024-29um) was also used to help with system alignment. A 90:10 beam splitter was inserted after the tube lens and the regular camera was mounted on the 10 percent reflection side. An Adafruit P4 LED array placed 84 mm behind the sample provided illumination. We used a 9×9 632 nm red LED array to provide a 0.18 illumination NA.

SPAD array cameras are typically used with high-powered pulsed optical sources for LIDAR, for example. In our case, we chose LED array illumination for demonstration simplicity. These LEDs are effectively continuous-wave and do not provide the high power required for the optimal SPAD array camera operation, especially at the large source-sample separation that our setup necessitated. Even in a standard FP setup, exposure times up to 1 s are required to achieve sufficient signal-to-noise ratio (SNR) when using such LED arrays [15]. Despite the single photon sensitivity of the SPAD array camera, its low fill factor (1.5%) still required the use of relatively lengthy exposure times (up to 1 s or more for dark-field). This was only possible in the photon counting mode, so we chose this mode for our initial demonstrations in the following experiments.

Before we performed FP with the SPAD array camera, we executed a series of operations to account for the SPAD array’s read noise, which are detailed in Appendix A. First, we captured a sequence of short-exposure images (3 μs) per illumination angle, as opposed to a single long exposure. This provided us with the opportunity to capture more photons and get higher SNR. Second, we also captured a matching dataset without a sample present, which we used for hot-pixel noise subtraction.

Completing these two steps per image capture session resulted in relatively high-quality 32×32×81 SPAD array image datasets that we subsequently processed with our quantized FP algorithm on several different targets (see Fig. 8). First, we imaged a USAF resolution test target to validate the quantitative improvement in SPAD array imaging resolution. Low-resolution HDR images were able to resolve group 7 element 3 (6.2 μm full-pitch resolution). Applying the standard FP reconstruction algorithm with this image data led to image reconstructions in which we could resolve group 8 element 6 (2.2 μm full-pitch resolution), showing an approximate 3× improvement in resolution. This closely matches the expected resolution gain for this FP setup. There are several primarily low-frequency artifacts in the reconstruction due in part to large detector read noise. The quantized FP reconstruction formed by processing binary SPAD array image data is also able to resolve group 8 element 6, albeit at a lower contrast and with larger noise-induced low-frequency artifacts, which is expected following our simulation results.

We also imaged a Siemens star resolution target and performed FP reconstructions in a similar manner, with results (Dataset 1, Ref. [39] and Code 1, Ref. [40]) shown in Figs. 8(c) and 8(d), which demonstrates omnidirectional resolution improvement. While resolution improvement is consistent in all directions, some artifacts still persist. A primarily stained transparent natural sample (a plant root cross section) is also imaged and shown in Figs. 8(e) and 8(f). It can be seen that the reconstructed phases for both reconstructions using the HDR dataset and the binary dataset approximately match. While preliminary, these results successfully verify the ability to improve image resolution using FP on SPAD array cameras. Even in the presence of relatively high amounts of noise, it is still possible to utilize binary image measurements from these single-photon detector arrays to jointly improve image resolution and retrieve the missing object phase. We hope that these results would improve with the ongoing improvements in the SPAD array image sensor technology. We adopted the pre-processing steps of background subtraction and LED position calibration. Other advanced pre-possessing methods, such as addressing aberrations and inter-LED intensity fluctuations [41,42], can be considered in future work.

SPAD array cameras are an exciting new type of single-photon sensitive sensors with high-speed time-gating capabilities. In this work, we studied the feasibility of performing Fourier ptychography with such sensors to jointly improve their imaging resolution and retrieve the unknown phase of the imaged field. To achieve this effectively on these new sensors, we have presented an FP forward model that accounts for the finite bit depth of an image sensor and image exposure selection, and a new quantized FP reconstruction algorithm that accounts for their finite bit depth during iterative reconstruction. This results in the surprising ability to form high-quality complex-valued image reconstructions using binary-valued image data. However, we stopped short of presenting an implementation that utilizes their time-of-flight capability. This extension would most likely require the use of alternative illumination sources (as opposed to LEDs), which would also help to address the imaging noise issues that our demonstration faced to potentially improve resolution enhancement performance. As an alternative to the SPAD array used in this work, large-format quanta image sensors [43] could also be adopted for future implementations of quantized FP. CMOS SPAD quanta image sensors can collect single-photon data with 61% fill factor and 100,000 frames per second [44], for example, which offers a significantly higher fill factor than the SPAD array utilized in our experiments. Such alternative sensors, as well as large-format SPAD arrays with small pixel pitch and high fill factor [45], are currently beginning to enter the commercial market and can also be utilized to collect highly quantized or single-photon image frames for quantized FP. As the pixel count and fill factor of such detectors continue to improve, we are hopeful that future implementations of quantized FP can improve upon our preliminary experiments in this work. There are several direct ways to extend our current demonstration to utilize time-of-flight measurements. Most directly, photons from a particular time range may be selected to form a set of binary images for depth-selective resolution improvement and phase measurement, which may prove particularly helpful when imaging within a scattering environment. Alternatively, one may utilize the time-of-flight image data to improve high-resolution 3D image reconstruction. For this latter goal, one may adopt a 3D imaging forward model for FP, which has been considered in prior work within the context of microscopy [46–48] but not yet within the more macroscopic regime with current SPAD technology’s millimeter-scale depth resolution. Given that FP can currently reconstruct 3D depth information both via phase measurement and by tomographic principles, however, adding in time-of-flight information offers an interesting possibility for future development.

Acknowledgment. The authors would like to thank Kernel Inc. for their kind assistance.