Correlated-photon imaging at 10 volumetric images per second

The correlation properties of light provide an outstanding tool to overcome the limitations of traditional imaging techniques. A relevant case is represented by correlation plenoptic imaging (CPI), a quantum-inspired volumetric imaging protocol employing spatio-temporally correlated photons from either entangled or chaotic sources to address the main limitations of conventional light-field imaging, namely, the poor spatial resolution and the reduced change of perspective for 3D imaging. However, the application potential of high-resolution imaging modalities relying on photon correlations is limited, in practice, by the need to collect a large number of frames. This creates a gap, unacceptable for many relevant tasks, between the time performance of correlated-light imaging and that of traditional imaging methods. In this article, we address this issue by exploiting the photon number correlations intrinsic in chaotic light, combined with a cutting-edge ultrafast sensor made of a large array of single-photon avalanche diodes (SPADs). This combination of source and sensor is embedded within a novel single-lens CPI scheme enabling to acquire 10 volumetric images per second. Our results place correlated-photon imaging at a competitive edge and prove its potential in practical applications.


Introduction
The very first studies on the peculiar properties of entangled photons from spontaneous parametric down-conversion (SPDC) 1,2 led to the discovery of quite counter-intuitive correlation imaging modalities 1,3,4 , and to the development of so-called quantum imaging.In this context, spatio-temporal correlations of light have been exploited both toward novel imaging schemes and to overcome the limits of traditional imaging (e.g., in terms of resolution, spectral range, contrast, signal-to-noise ratio) [5][6][7] .On the one hand, quantum correlated photons have been employed to achieve otherwise unattainable quantum effects, such as imaging with undetected photons 8 and sub-shot-noise imaging [9][10][11][12] .On the other hand, it was discovered that specific correlation imaging tasks could still be achieved when replacing SPDC photons with chaotic light sources, such as thermal and pseudo-thermal [13][14][15][16][17] .Such classically correlated sources typically entails a worse signal-to-noise ratio (SNR) than expected with entangled photons 18,19 , but also brings in several practical advantages: light sources are simpler and more feasible (potentially, any available incoherent source could serve the purpose 20,21 ), the image acquisition time can be significantly reduced by avoiding the low production rate of entangled photons, and protocols are potentially insensitive to the deleterious effects of decoherence.A wide range of quantum-inspired imaging protocols based on chaotic light correlations have been demonstrated so far, in different application scenarios: imaging of objects hidden to the main sensor 22,23 , dual wavelength imaging 24,25 (analogous to the implementation with entangled photons 26,27 ), detection of objects surrounded by turbulence [28][29][30] , 3D imaging through computational ghost imaging 31 , refocusing and 3D imaging through correlation plenotpic imaging and microscopy [32][33][34][35][36][37][38] .
However, the variety of advantages entailed by the spatio-temporal correlation properties of light clashes with the most relevant challenge of correlation imaging methods: reducing the acquisition time to make them effectively competitive with the corresponding state-of-the-art traditional techniques.Such a challenge is intrinsic to this approach, since evaluating correlation functions requires sampling a high enough number of statistical realizations.The use of high-resolution CCD and CMOS cameras in correlation imaging setups 6,7,13,16,34,38,39 has enabled to overcome the extremely time-consuming mechanical scanning of the image plane, as performed in pioneering correlation imaging experiments 3,14,17 .There are two time scales that one must consider when dealing with the speed of correlation imaging with a camera 7 .First, if N t camera frames are required, the total acquisition time is equal to T image = N t /R, with R the frame rate; its inverse T −1 image is the rate at which correlation images can be acquired.The acquisition rate can be reduced by i) maximizing R, as in high-speed sensors, and ii) optimizing the trade-off between number of frames and SNR, with the latter expected to depend on √ N t 18, 40-43 .The other crucial time scale to be considered is the gating time, namely, the effective sensitivity window in which a single frame is acquired: if it is larger than the source coherence time, intensity fluctuations in each frame are partially erased 44 , making it more difficult to reconstruct their correlations, and thus increasing the required N t .An innovative class of detectors, made of an array of single-photon avalanche photodiodes (SPAD) [45][46][47][48] , provides one of the most interesting solutions available to take on both temporal issues.SPAD arrays ensure at once fast acquisition rates of up to about 10 5 fps, sub-10 ns gating time, and low noise: SPAD arrays are characterized by the absence of readout noise, whereas dark counts are limited to less than 10 counts per second per pixel.The reduced noise is essential for high-quality sampling of the light statistics with fewer frames.SPAD arrays have thus been extensively employed for correlation measurements, and yet more rarely for correlated-photon imaging (see Ref. 48 for a wide review and detailed bibliography).In this context, despite the development of high-resolution SPAD arrays has permitted to largely cut the acquisition time by exploiting multiple pair coincidences within a single frame, the acquisition times attained in entangled-based high-resolution 2D imaging remain larger than one second 49,50 .In fact, even neglecting the data bandwidth limitation of state-of-the-art high-resolution SPAD arrays, which currently prevents large arrays from being used at full speed for extended periods, the long acquisition time needed by entanglement-based correlation imaging is essentially determined by the low production rate of SPDC.
In this work, we shall demonstrate correlated-photon imaging at a rate of 10 volumetric images per second, with a field of view of 256 × 256 pixels.The key for achieving such performances is the use of: 1) SwissSPAD2 [45][46][47]51 , a SPAD array operated at a 512 × 256-pixel resolution, and capable of acquiring up to 10 5 frames per second; 2) chaotic light illumination, which enables to avoid the low speed related with the low production rates of parametric down-conversion; 3) the principles of correlation plenoptic imaging (CPI) 32,[35][36][37]52 , a technique that enables to acquire information about both the spatial distribution of light and its propagation direction; such a richness of information is used, during data processing, to reconstruct threedimensional snapshots of the scene of interest.In traditional light-field imaging devices [53][54][55][56] , based on the insertion of an array of micro-lenses before the sensor, directional resolution can be gained only at the expense of spatial resolution 57 .Nonetheless, aided by additional post-processing algorithms [58][59][60] , they find applications in the most diverse fields, such as photography 61 and microscopy 62 .Plenoptic or light-field cameras based on intensity measurement can currently operate at the essentially same rates of a standard camera, which evidently outperforms any correlation imaging technique (whether volumetric or not) at similar resolutions.However, many interesting applications of light-field imaging in science, including the neuronal activity detection accomplished in Refs.63,64 , are performed at image rates ranging between 10 and 100 Hz 64,65 .The results presented in this paper thus demonstrate the competitiveness of correlated-photon imaging in these tasks: CPI based on correlated photons from chaotic sources brings a net advantage with respect to state-of-the-art light-field cameras, combining similar time performances with an improved volumetric resolution at analogous numerical aperture 36,66 .To the best of our knowledge, the presented CPI setup is indeed the first case of a feasible light-field device employing a single lens.
The article is organized as follows.In the "Results" section, we summarize the working principle of the novel CPI device and show the possibility to obtain light-field images with a 0.1 s acquisition time.In the "Discussion" section, we highlight the relevance of our results along with their implications for practical applications, and discuss both the current limitations of our experiment and the foreseen future developments.In the "Materials and Methods" section, we present details about the experimental setup.The Supplementary Information document further details on the measured correlation functions and the refocusing process, as well as complementary results on the interdependence between the number of collected frames (hence, the image acquisition time) and the SNR.

Results
The working principle of the CPI camera is reported in Fig. 1 36 .Two planes O a and O b , chosen arbitrarily within the threedimensional scene of interest, are focused on two high-resolution sensors, D a and D b .Unlike a conventional light-field camera, involving both the usual camera lens and a micro-lens array, our CPI device is realized with a single lens, collecting light from both chosen planes, and focusing them on the two sensors.Light from the scene is chaotic; hence, by computing the equal-time pixel-by-pixel correlation between the number of photons (N a and N b ) detected by the sensors D a and D b , we obtain the correlation function: where ⟨. . .⟩ indicates the averaging process, while ρ ρ ρ a and ρ ρ ρ b are the coordinates identifying pixel positions on the sensors.The correlation function in Eq. ( 1) contains plenoptic information, and thus enables reconstructing features of a 3D object that can be placed both between and beyond the two planes O a and O b imaged on the detectors 35,67 .As explained in detail in the Supplementary Information, Γ(ρ ρ ρ a , ρ ρ ρ b ) encodes a collection of multi-perspective volumetric images; proper processing of these volumetric images provides the refocused image of a specific transverse plane in the scene.Adopting chaotic light illumination entails that the magnitude of the correlation function ( 1) scales like the product on the mean number of photons 44 , thus being crucial in ensuring that correlation measurements provide analogous results both in the case of high-intensity, as in Ref. 38 , and in the single-photon regime, as in the present work.However, a key requirement in this sense, is that the SPAD array works in the linear regime (namely, the probability to detect a photon is proportional to the intensity of the impinging field), far from saturation.
Experimental results are reported in Fig. 2: Both sensors acquire blurred images of three different planar test targets (A, B, and C); the plenoptic information contained in the measured correlation function enables reconstructing the object details, in all three cases.In the panel on the left, we report the out-of-focus images of the test target, which is placed either within (case B) or outside (cases A and C) the volume defined by the two conjugate planes of the detectors; the effective refocusing enabled by CPI is shown in the center panels.The recovery in visibility deriving from refocusing is demonstrated in the right panels: here, we compare the linear images related with both average intensity and CPI, which are obtained by integration along the slit direction.CPI also enables over 10× depth of field (DOF) enhancement at a resolution of 250 µm, and 12× at 160 µm, with respect to a conventional imaging system with the same numerical aperture (NA).This can be seen by considering the curves reported in panel (b) of Fig. 2, together with the axial position and the distance between neighboring slits on the three test targets (A, B, and C) reported in panel (a).The plot shows the expected resolution limit of the refocused images (blue line), with varying axial position z, compared with the analogous limits associated with the conventional images focused on D a (green line) and D b (orange line).In particular, the blue line indicates the object detail size (i.e., the resolution) that can be refocused by our CPI device with 10% visibility, as a function of the longitudinal distance (z) of the object from the lens.The green and red lines represent the natural DOF (as determined by the circle of confusion) of the images separately observed on the two sensors D a and D b , at the given resolution.
The images reported in Fig. 2 have been obtained at an overall acquisition speed of 10 volumetric images per second.This is an unprecedented result in the field of correlation imaging, and indicates the feasibility of correlated-photon imaging at video rate.The noise analysis reported in the Supplementary Information shows the robustness of the developed technique: when the acquisition speed is reduced to 1 image per second by increasing the number of acquired frames by 1 order of magnitude, the SNR increases by nearly 35%, and reaches its maximum value.A comparison between the CPI images acquired at the speed of  1 and 10 volumetric images per second is reported in Fig. S3 of the Supplementary Information.

Discussion
We have presented a quantum-inspired imaging system capable of collecting 10 plenoptic images per second.Since plenoptic images are a collection of multi-perspective volumetric images, they enable changing, in post-processing, the focusing plane within the entire axial range enclosed by the blue curve in Fig. 2b).This result entails a large reduction of the gap in timeperformance between CPI and conventional light-field imaging, that is generally performed at a speed between 10 and 100 Hz, in scientific applications 64,65 , but with a significant loss of resolution due to microlens array and intensity measurement.
The key element to achieve such a critical improvement in the acquisition time is the integration of the SwissSPAD2 sensor in a chaotic-light based correlaton plenoptic imaging setup.This SPAD array enables to collect, with single-photon sensitivity, all the frames that contribute to the plenoptic correlation image, at a rate of almost 100.000frames per second.Such a fast rate is combined with both a resolution comparable to that of ordinary detectors and low noise (see Refs. 46,47,51 for a detailed description of the sensor).The low noise of the detector is a key aspect for keeping as low as possible the number of frames N t required for reconstructing light statistics and correlations.Our SwissSPAD2 sensor has an on-board DDR3 memory bank (2 GB) that can be filled with a maximum of 131 072 measured binary frames.By saving the acquired data on the internal memory instead of streaming to an external disk, we were able to exploit the maximum speed of 97.7 kHz at full resolution.However, the limited capacity of the memory has bound us to single-image acquisitions instead of videos.In the future, we shall employ a new generation of SwissSPAD2 capable of streaming data from a 512 × 512 sensor to a workstation at full speed.
The new generation of SwissSPAD2 also provides a relevant improvement in the gating time, which can be reduced down to about 10 ns.This represents an extremely important parameter in correlation imaging based on chaotic light, since reconstruction of light statistics is optimal when the detector exposure time matches the coherence time of light 44 .The possibility to match coherence times as small as hundreds of ps would open the way to CPI with broadband sources, thus leading the way toward passive quantum imaging devices.
It is reasonable to expect the achieved acquisition speed to be further increased through computational techniques enabling to use less frames to achieve a comparable SNR; examples are compressive sensing 68,69 , quantum tomography 70 , and machine learning 71 .All these techniques are currently being developed in the framework of CPI, and we plan to integrate them with our refocusing algorithm 51 .To further increase the SNR while reducing the acquired number of frames, we are also working toward employing, within the data analysis, the statistical properties of the correlation function, in a similar fashion as in 72 .It is interesting to emphasize, however, that all data presented in this work have not been treated with any denoising algorithm, or post-processing method, other than the refocusing algorithm described in the Supplementary Information.
The novelty of the implemented CPI setup also stands in the fact that two arbitrary planes are focused on the two sensors 36 , as opposed to conventional approaches involving imaging of the main lens for retrieving directional information.This approach enables: (i) parallel acquisition of two diffraction limited images within the three-dimensional scene of interest, (ii) single-lens light-field imaging, which is quite significant considering the disadvantages and physical limitations connected with the use of micro-lenses (i.e., resolution loss and reduced 3D imaging capability), (iii) a DOF enhancement by over 1 order of magnitude, without sacrificing diffraction-limited resolution.

Materials and Methods
The optical system is aimed at maximizing speed of acquisition and performance in terms of resolution versus DOF trade-off, while guaranteeing flexibility in the focusing capability.The design of the employed CPI setup is oriented to the acquisition of generally demagnified images, as in an ordinary camera.The lack of two synchronized SPAD arrays has imposed using two halves of a single SPAD array as the two sensors D a and D b .This entails some constraints to the setup design: demagnified images can only be obtained if the two CPI paths are separated upstream of the lens, rather than downstream (as reported in Fig. 1, and originally proposed in Ref. 35 ).The experimental setup thus consists of two main parts (Fig. 3): 1. the ultra-fast imaging device, made of a camera lens (Navitar MVL75M1, of focal length 75 mm and focal ratio 2.8) mounted on the SPAD array sensor; 2. a "CPI adapter", represented in Fig. 3 by the dashed gray rectangle.
The CPI adapter endeavors the ultra-fast camera with plenoptic properties by first creating (through the first polarizing beam-splitter) and then recombining (through the second polarizing beam-splitter) two optical paths, which we shall indicate as a (depicted in green) and b (depicted in red).Each optical path contains a delay line, offering the required flexibility for choosing the two arbitrary planes O a and O b , when preparing the acquisition.In our setup the distances between the planes O a,b and the lens are z a = 345 mm and z b = 293 mm, respectively.The delay lines are made by the combined system (QM a,b ) of a quarter-wave plate (QWP) and a mirror.This combined system converts light from H-polarized to V-polarized, and viceversa, so that the beam that is back-reflected by QM a,b is then reflected/transmitted by the corresponding PBS toward the camera lens.Changing the optical path in arms a and b defines the specific plane to be imaged on sensor D a and D b , respectively.In fact, given the lens focal length f and the fixed lens-to-sensor distance z i , the distance z o of the object plane from the lens is uniquely defined by the thin lens equation.Hence, the two planes O a and O b , imaged on D a and D b , respectively, are both placed at an optical distance z o from the lens; however, the actual planes that are imaged on two disjoint halves of the sensor, are determined by length of the delay lines, which enable to arbitrarily choose the distances z a and z b , associated with two different planes within the volume of interest.We should specify that the versatility in choosing the two planes is a useful feature when setting up the acquisition, since it allows the experimenter both to select the planes to be focused and to define the volume that can be refocused (as defined by the blue curve in Fig. 2(b)); the specific choice, however, plays no role during the acquisition itself.Both delay lines are characterized by the same magnifications M = −z i /z o , numerical aperture (NA), and resolution at focus, as defined by the camera lens.The clear aperture of both the polarizing beam splitter, PBS (45 mm), and the optics (2 inches) in the delay lines have been chosen to enable fully exploiting the NA of the camera lens.In order to maximize its fill factor, the SwissSPAD2 sensor is equipped with a microlens array; its NA ≈ 0.25 is larger than the lens NA on the image side (NA i = NA o /|M| = 0.13) and does not limit the NA of the CPI device.However, the pixel size of the sensor is larger than the achievable diffraction limited resolution; hence, the setup has a pixel limited resolution of 95 µm.
In the present experiment, the CPI device was employed to image transmissive planar test targets placed out of focus, as shown in Fig. 2. The targets are illuminated by a chaotic light source of controllable polarization, intensity, and coherence time, made by a green diode laser (Thorlabs CPS532, λ = 532 nm) scattered by a rotating ground glass disk (GGD).At the maximum rotation speed of the GGD (30 Hz), the measured coherence time of the source is t ch ≃ 15 µs.
SwissSPAD2 employs a design that provides one of the largest resolutions (512 × 512 photodiodes operating in Geiger mode) as well as one of the highest sensitivity (50% photon detection probability at 520 nm) and lowest dark count rate (0.26 cps/µm 2 , equivalent to a median value of less than 10 cps per pixel) combinations among SPADs which are built with standard CMOS-process technologies.Its 10.5% native fill factor is improved by 4-5 times, for collimated light, by means of the use of a microlens array.The output of each frame consists of a binary matrix identifying the pixels that have been triggered by at least one photon.Due to the binary nature of the signal, it is of utmost importance for the reconstruction of intensity correlations to work close to the linear regime 45,73 , in which the probability to detect a photon is proportional to the intensity of the impinging electromagnetic field.

Supplementary information S1 Plenoptic properties of the intensity correlation function
As demonstrated in Ref. 36 , plenoptic information is encoded in the correlation function of Eq. ( 1) in the main text, representing the correlation between the intensity fluctuations reaching two points, of which one is placed on the detector D a , and the other on the detector D b .The correlation function reads, up to irrelevant factors that do not depend on either ρ ρ ρ a or ρ ρ ρ b , where A is the aperture function of the object.The function Ψ can be considered as a "second-order point-spread function", which determines the correspondence between object points and detector points.By referring to the parameters in Fig. 3 of the main text, and assuming the source emits chaotic light with an average intensity profile S(ρ ρ ρ s ) and negligible transverse coherence, we have with the functions p j , j = a, b, describing field propagation from the object to the detector.Calling P the pupil function of the lens, we obtain p j (ρ ρ ρ o , ρ ρ ρ j ) = P(ρ ρ ρ ℓ j )e ik z φ j (ρ ρ ρ o ,ρ ρ ρ ℓ j ,ρ ρ ρ j ) d 2 ρ ρ ρ ℓ j (S4) with Notice that the differences in phases φ j with respect to the original results described in Ref. 36 are due to the fact that the optical paths are split upstream of the lens, and not downstream.The plenoptic properties of the correlation function are easily deduced, for example, from the fact that, by applying a stationary-phase approximation to the integrals which defines it, one obtains Therefore, Γ(ρ ρ ρ a , ρ ρ ρ b ) represents a collection of images of the object, whose shifts and scaling depend on the axial position of the latter.The refocusing process, consisting in properly realigning the images contained in the correlation function, is discussed in the next section.

S2 Refocusing and correlation aperture
In a ray-optics approximation, the correlation function that is measured reads 36,74,75 : This equation shows that Γ contains information regarding the object aperture A, and that also the lens aperture P and source intensity profile S have effects on the correlation function.Each aperture defines a region inside the correlation space (ρ ρ ρ a , ρ ρ ρ b ), as shown in Fig. S4: the green and red lines indicate the regions where non-zero correlations can be expected, as defined by the lens; the blue lines, instead, indicate the effects of limited aperture of the source, i.e. the width of the Gaussian profile illuminating the GGD.In fact, correlation can only be measured if the light passes through all the apertures, and this means that the correlation region is the intersection of the three correlation regions defined by P and S. By considering only the correlation can be done by limiting the integration domain to an area defined by the source radius.By comparing Eqs.(S7) and (S11), we see that the coefficients on which the object aperture depends correspond to the first row of the refocusing matrix, while the second row contains the coefficients that appear in the source profile.The first line gives us an insight into what the refocusing algorithm actually does, that is, rescaling the detector coordinates so that, in the transformed plane, the object depends only on the set of two-dimensional coordinates ρ ρ ρ r .It is also clear that the second line, defining the variable ρ ρ ρ s that is integrated in Eq. (S13), plays no relevant role for the object reconstruction, and can be chosen arbitrarily.However, some choices are more convenient than other when defining the integration variable.For example, since the extension of the correlation function is defined by the size of the apertures, a clever choice is to define ρ ρ ρ s as the transverse coordinate of one of the limiting apertures.By doing so, one can integrate only where non-zero correlations are expected, and limit the integration of Eq. (S13) to a much smaller area than what has been measured, so as to speed up calculations quite dramatically.Given the values of the correlation apertures in Eq. (S10), it made sense to us to choose the integration variable as the one defined by the source.If that were not the case, the second line would have been modified with the coefficients of the aperture defining the major limitation on the correlation area.

S3 Data analysis workflow
Data analysis has been performed on MATLAB and consisted of two main parts, that are data reading and correlation, and refocusing.In the first part, the binary frames are read from the disk and the correlation function is evaluated by averaging over all frames.Computation of the correlation function is fast and efficient, so that the speed of the process was mostly defined by the data read rate of the disk; reading ∼ 10 4 frames (512 × 256) and computing the correlation function took about 200 s.After the correlation function is available in the workstation memory as a 4D array, the operation of refocusing took about 14 s per axial coordinate z.

S4 Study of the SNR dependence on the number of frames
We report in Fig. S5 the full analysis of the dependence of the SNR on the number of frames N t .The analysis has been performed in the experimental conditions corresponding to case "C" of Fig. 2 (main text).The experimental points indicates that the SNR tends to saturate as the number of frames increases.This behavior deviates from the √ N t scaling, expected if the 9/13 which accounts for both the √ N t behavior at a small number of frames and the saturation at a high number of frames.The best fit provides a = 3.60 × 10 −2 and b = 3.18 × 10 2 .The plot in Fig. S5 indicates that variations in the number of acquired frames around N t = 4.0 × 10 5 does not yield an appreciable improvement in image quality.Actually, a reduction of N t by a factor close to 4 (N t = 9.8 × 10 4 ) with respect to the whole dataset, which enables to collect each CPI image in one second (T CPI = 1 s), produces a decrease in the SNR by less than 5%.Moreover, a further reduction in the number of frames to N t = 9.8 × 10 3 , which enables to collect 10 CPI in one second (T CPI = 0.1 s), Fig. S6 shows a comparison between these three cases, with the corresponding estimated values of the SNR.The SNR in the region of interest enclosed in the cyan rectangle highlighted in Fig. S6, has been estimated as where Σ in represents the average value of the signal, in the image, in correspondence of the transmissive parts of the considered five-slit group; the denominator ∆ in Σ represents the standard deviation of the values contributing to the numerator.Such a definition relies on the realistic assumption that the statistical distributions of the signal in the refocused image, in correspondence of the transmissive parts, are identical.

Figure 1 .
Figure 1.Working principle of the developed CPI protocol.O a and O b are the two conjugate planes of the high-resolution sensors D a and D b ; they are placed at the generic distances z a and z b from the lens, respectively.BS is a beam splitter sending light from the lens toward the two sensors.Pixel-by-pixel correlations between photon number fluctuations are evaluated by software and employed to reconstruct the volumetric image of the scene.See "Materials and Methods" for the detailed experimental setup.

Figure 2 .
Figure 2. (a) Comparison between the images directly retrieved by the sensors, through their average intensity (left panels), and the images refocused by means of CPI (central panel); the plots in the right panel are obtained by integrating the highlighted rectangles in the corresponding average intensities (gray) and refocused images (blue) along the slit direction.All the reported data have been obtained by acquiring N t = 9.8 × 10 3 frames, at full resolution, at ∼ 9.8 × 10 4 frames per second, resulting in an overall acquisition speed of 10 volumetric images per second.(b) Comparison of the resolution achieved by conventional imaging (orange for D b , green for D a ) and CPI (blue) as a function of the axial distance from the lens.A (z = 275 mm), B (z = 319 mm) and C (z = 373 mm) indicate the three different masks shown in panels (a); the plot shows their axial position and the distances between neighboring slits within the masks.

Figure 3 .
Figure 3. Technical scheme of the developed CPI setup.The chaotic source is made of a diode laser illuminating a rotating ground glass disk.Two planes, O a and O b , arbitrarily chosen in the surrounding of the scene of interest, are imaged by a unique lens onto two disjoint high-resolution detectors D a and D b , which are practically implemented by using two halves of the same SwissSPAD2 sensor.Two optical paths, one for each detector, are realized by means of two polarizing beams splitters (PBS), two quarter-wave plates (QWP) and four mirrors; each pair of QWP and mirror is mounted on a translation stage, which offers flexibility in the choice of the two planes O a and O b , when preparing the acquisition.

Figure S5 .
Figure S5.Behavior of the SNR in the image refocused by CPI (referred to case C shown in Fig.3of the main text), as a function of the number of collected frames N t .Blue circles represent the SNR evaluated on the experimental refocused images, while the solid black line corresponds to the fitting curve (a + b/N t ) −1/2 .The points discussed in the text, and corresponding to N t = 4.0 × 10 5 , N t = 9.8 × 10 4 , and N t = 9.8 × 10 3 , are highlighted in black, magenta, and green, respectively.The red dotted lines correspond to N t = 8 × 10 4 , where we estimated a decrease of the SNR by 5% of its maximum.

(a) N t = 9 . 8 × 3 Figure S6 .
Figure S6.Comparison of the refocused images acquired at 10 (a) and 1 (b) volumetric images per second.Panel (c) shows the refocused image using the whole acquired dataset.The cyan rectangle indicates the area where the SNR has been evaluated on the three images.