Compressed imaging with focused light

The flowchart in figure 1(b) summarises the modelling sequence implemented for two-dimensional (2D) computational super-resolution imaging based on compressed imaging and deconvolution. First, we generated different illumination bases ill: focus illumination, speckle illumination, and array illumination. Focus illumination consisted in a sequence of single diffraction-limited focus at different lateral positions in the sample, which is equivalent to standard point-scanning imaging (supplementary video 1). Speckle illumination consisted in a sequence of unique speckle-like intensity distribution covering the surface of the image (supplementary video 2). Array illumination consisted in a sequence of multiple diffraction-limited foci at different lateral position in the sample following a structured arrangement (supplementary video 3). All illumination bases were generated by modelling propagation through a MMF. In some cases, perturbations were introduced in the illumination ill_per in the form of optical aberrations or bending of the fiber. For the latter, it was necessary to calculate new sets of illumination bases. However, ill_per was equal to ill in most cases.

In parallel, we also generated artificial 2D bead samples (Object A). Then, we generated the noisy signal S_σ . For compressed imaging, each measurement S_i in the signal S was calculated by integrating all the elements in the matrix resulting from the element-wise multiplication between A and ill_per. For deconvolution, A was convolved with a single focus from the focus illumination to generate S. S_σ was generated by adding varying amounts of zero-mean Gaussian white noise with a specified variance σ.

Next, we performed computational super-resolution imaging using CS and deconvolution algorithms. The super-resolution image B was obtained based on two primary inputs: the noisy signal S_σ and the unperturbed illumination ill. In other words, the algorithms had no knowledge of any perturbations, noise level, or ground truth A. Compressed image reconstruction was performed using the basis pursuit algorithm. Deconvolution was performed using a modified Lucy-Richardson (LR) algorithm with total variation regularisation. Finally, the performance of the SR image reconstruction process was assessed quantitatively by comparing the SR image B to the ground truth object A by calculating the correlation coefficient R.

Using the methodology outlined above, the performance of compressed imaging with the different illumination bases and deconvolution was assessed for varying amounts of noise. The use of focus and speckle illumination was also compared for compressed imaging in the presence of perturbation. The three-dimensional (3D) analysis is described in details below.

Synthetic imaging data was generated by simulating 2D and 3D bead images. The 2D bead images were simulated by producing an array of 2D points at random locations using a salt and pepper noise function with a density of 0.001 on a $120\, \times\, 120$ null matrix. The beads were then restricted within a given radius from the centre of the image (53 pixels). Finally, the intensity of the beads was randomised following a Gaussian distribution (σ = 0.5) and all beads were given the same diameter using Gaussian filtering (σ = 1.33). A total of 21 2D images were generated and the width, given as the full width at half maximum, of each bead objects was 1.44 µm. Similarly, 3D bead images were simulated by producing an array of 3D points at random locations using a salt and pepper noise function with a density of 0.0015 on a $40 \times 40 \times Z$ null matrix, where Z is the number of axial planes (Z = 60). The beads were then restricted within a radius of 16 pixels from the image centre and between the 8th and 52nd axial planes. Finally, the beads were given the same diameter using Gaussian filtering (σ = 0.93). Five 3D images were generated.

To facilitate the comparison between speckle and focus illumination, a MMF-based endo-microscope system was simulated. This system enabled the generation of both speckles and foci in the sample by giving different shapes to the wavefront entering a straight fiber at the proximal end. The original MATLAB code that was modified for this work can be accessed at [24] and is based on the model described in Plöschner et al [25].

Foci were generated in the sample by first evaluating the transmission matrix using a sequence of proximal foci and then calculating its regularised inversion. A proximal grid of $49{\,}\times{\,}49$ pixels and a distal grid of $120{\,}\times{\,}120$ pixels were used for the simulation. The resulting illumination matrix was then sub-sampled to select 2401 foci uniformly covering the sample over the surface of the core. The focus illumination basis is shown in supplementary video 1.

Speckle patterns were generated in the sample by focusing light at different location of the MMF proximal end. The resulting speckle patterns were determined also by using the transmission matrix. A proximal and distal grid of $120 \times 120$ pixels were used for speckle simulation. The resulting illumination was then sub-sampled to 2401 illuminations uniformly covering the proximal core area. The speckle illumination basis is shown in supplementary video 2.

For all 2D analyses, light at 532 nm was modelled going through a MMF having a diameter of 50 µm, a numerical aperture of 0.22, and a length of 1 m. Normalisation was applied such that the total intensity was the same for both illumination strategies. For bending, a single fiber segment was considered with a curvature defined by a bending radius $r_\textrm{bending}$ [24].

For compressed imaging with either illumination basis, the signal was calculated by integrating the matrix resulting from the element-wise multiplication of the 2D illumination and 2D sample for each illumination composing the basis. For deconvolution, the object was convolved with the central focus. For all reconstruction methods, a zero-mean Gaussian white noise with a specified variance σ was superimposed onto the signal when needed. To characterise the quality of the reconstruction process, the correlation coefficient R, the mutual coherence µ, and the error probability $P_\textrm{e,worst}$ were calculated as described in the supplementary information.

For bending, the signal was calculated by pixel-wise multiplication of the 2D bent illumination and 2D sample for an illumination in the sequence and then integrated. The SR image reconstruction processed was agnostic to the bending and assumed illuminations from an straight fiber. Information regarding the modelling with optical aberrations can be found in the supplementary information.

Analyses for volumetric imaging were implemented for focus illumination only. Distal foci in 3D were created by first generating 2D foci as described in the previous section and then propagating them axially using a defocus function. A proximal grid of $49{\,}\times{\,}49$ pixels and a distal grid of $180 \times 180$ pixels were used for the simulation. The resulting illumination matrix was then sub-sampled by a factor of 2 to select foci uniformly covering the sample over the surface of the core. For all 3D analyses, light at 930 nm was modelled going through a MMF having a diameter of 50 µm, a numerical aperture of 0.37, and a length of 1 m. The signal was calculated by integrating the matrix resulting from the element-wise multiplication of the 3D illumination and 3D sample for each illumination composing the focus illumination basis. The outcome of SR image reconstruction was evaluated by calculating the axial and lateral full width half maximum (FWHM) of beads. The lateral and axial profiles were manually obtained for two randomly selected beads in each of the 3D reconstructions. The FWHM was then calculated by finding the four points neighbouring 50$\%$ of the peak intensity value on each side of the peak and performing a linear interpolation between the two points to determine the positional values of both half maximum [26].

For compressed image reconstruction, we used the basis pursuit (BP) algorithm [8] with a maximum number of iterations set to 3000 and a modified calculation of pseudo-inverse calculation for additional precision [27]. More information about the theory of compressed sensing is available elsewhere [7–11]. The augmented Lagrangian and over-relaxation parameters were both set to 1.0. For deconvolution, we used the modified LR algorithm with total variation regularisation [5, 22, 28, 29]. A uniform point-response was assumed across the surface of the fiber, as obtained in the focus illumination basis. The algorithm performed nine iterations, which was empirically determined to provide the best performance.

The performance of the three different SR image reconstruction methods was first compared in the noiseless scenario with the 2D bead data. Figure 2(a) shows an example set including the ground truth object, the diffraction-limited image, and the reconstructed SR images using the BP algorithm with speckle and focus illuminations and using the LR deconvolution algorithm. Extended spatial resolution was achieved with all SR image reconstruction methods when compared to the diffraction-limited image (figure 2(a)). Use of the BP algorithm with any illumination resulted in no visually perceptible difference between the SR image and the ground truth object (figures 2(a) and (b)). When using LR deconvolution, small errors in the reconstruction were visible when tracing line profiles but overall the reconstruction quality was still high, which is expected for a noise-free signal (figure 2(b)).

Zoom In Zoom Out Reset image size

The correlation coefficient R between the ground truth object and the SR images was calculated for the analysis to be quantitative. For the original bead density (0.01), R was above 0.998 with the BP algorithm and below 0.997 with the LR deconvolution algorithm in accordance with the above observations. A denser (0.1, 10×) and sparser (0.001, 0.1×) dataset were also generated to assess the effect of sparsity on the quality of the reconstruction. As expected for compressed imaging, R increased with increased sparsity (figure 2(c)). The same dependence was observed for the LR algorithm. However, while both illumination strategies yielded similar R value for all bead densities with the BP algorithm (figure 2(c)), the LR algorithm appeared to be more sensitive to sparsity as exemplified by the much lower R value of 0.9816 ± 0.0016 for the denser objects. In summary, all SR image reconstruction methods were successful when no noise was present in reconstructing images with increased spatial resolution, but use of the BP algorithm provided the optimal results regardless of the chosen illumination pattern. The results from figure 1(c) also confirm that the samples were sufficiently sparse in the spatial domain to meet the sparsity condition of CS algorithm and as such no basis transformation was necessary.

Next, the impact of noise for the three different SR image reconstruction methods was assessed by processing the dataset with varying amount of noise. The three methods had strikingly different noise dependence as shown by the graph of R with respect to the ground truth object as a function of noise (figure 3(a)). By fitting a logarithmic logistic function, the 50%-R value, denoted τ_noise, was found to be $5.2 \times 10^{-15}$, $2.6 \times 10^{-11}$ and $4.9 \times 10^{-6}$ for the LR deconvolution algorithm, BP algorithm with speckle illumination and BP algorithm with focus illumination, respectively. It is of note that the BP algorithm with focus illumination utilised the same type of illumination but only a fraction of the measurements compared to LR deconvolution, yet outperformed the latter. To explore this result further, the BP algorithm was used to reconstruct SR images with varying noise using the same illumination, including the number of measurements (no compression), as for LR deconvolution (figure 3(a)). A τ_noise value of $1.2 \times 10^{-12}$ was obtained, which is larger than for the LR algorithm on the same data. This observation suggests that CS algorithms might be preferable over LR deconvolution even when there is no data compression.

Zoom In Zoom Out Reset image size

More importantly, figure 3(a) shows that, for a given CS algorithm, the robustness to noise will depend on the illumination. Indeed, unlike for the noiseless case, the robustness of the BP reconstruction algorithm to noise was better for focus than speckle illumination, given equivalent data compression (figures 3(b) and (c)). In addition, focus illumination with compression outperformed the no compression case when the total incident power was equal. Together, these results suggest that the structure of the illumination and the efficiency of spatial sampling should both be considered as key factors. The comparative noise analysis between focus and speckle illumination using the BP reconstruction algorithm was repeated on images of extended objects. For this purpose, 21 images of hand-written numbers from the MNIST dataset were processed; no parameters were changed. The τ_noise was found to be $1.2 \times 10^{-9}$ and $4.8 \times 10^{-5}$ for the BP algorithm with speckle illumination and BP algorithm with focus illumination, respectively (figure 4(a)). This result, further illustrated by example SR images in figures 4(b) and (c), confirms that focus illumination provides more robustness to noise than speckle illumination when used with a CS algorithm, even when applied to more complex images.

Zoom In Zoom Out Reset image size

The source of the dependence on the illumination when processing noisy data with CS algorithms was investigated by using a third type of illumination with the BP algorithm. This third illumination consisted of structured arrays of foci (supplementary video 3). The number of foci per array was varied such that a range of τ_noise values, between that of single foci and speckle illuminations, were obtained when processing the data and fitting the logarithmic logistic function (figure 5(a)). The mutual coherence µ was then calculated to quantify the incoherence condition. This condition dictates that a maximally incoherent sampling provides an optimal sampling of the space and, accordingly, a low mutual coherence is preferred for the illumination [12]. In line with this principle and the above observations, ${\mu_{\textrm{BP}}}_{\textrm{focus}}$ was lower than ${\mu_{\textrm{BP}}}_{\textrm{speckle}}$ (figure 5(b)). Sampling incoherence is not necessarily the only factor at play in explaining the difference between focus and speckle illuminations because array illumination with 8 or more foci per illumination have a higher µ than speckle illumination, yet their τ_noise is higher (figures 5(a) and (b)).

Zoom In Zoom Out Reset image size

Another important factor to consider is the probability of making an error during the reconstruction. This probability is dependent on the illumination as well as the object and is not easily calculated. The higher bound of the error probability $P_\textrm{e,worst}$ serves as a good proxy metric as its calculation is tractable and can be expressed in an object-independent manner for the relative comparison of different illuminations [30]. Figure 5(c) indicates that focus illumination has a lower $P_\textrm{e,worst }$ than speckle illumination and that as the number of foci in array illumination increases $P_\textrm{e,worst,array}$ approaches $P_\textrm{e,worst,speckle }$. This can be intuitively understood by considering that a fixed amount of noise is added to a decreasing peak signal as the number of foci in the array illumination increases (figure 5(d)). In brief, several factors may contribute to the increased robustness to noise with focus illumination compared to speckle illumination, including improved space sampling and decreased error probability.

Several approaches have been implemented to achieve 3D compressed imaging with the goal of minimising photo-bleaching and increasing imaging speed [31–33]. One approach for 3D compressed imaging consists of acquiring data at different axial planes with a series of known 3D speckle patterns and then reconstructing a 3D image with a CS algorithm. The robustness to noise would again be improved with focus illumination as the reconstruction does not depend on the image dimensionality. Indeed, the signal is treated as a one-dimensional vector by the CS algorithm, even when objects are 2D and 3D. The computational requirement of this approach scales with the cube of the number of axial planes as the number of required measurements increases linearly and the size of the illumination matrix with the square of the number of axial planes. The entire volume also needs to be sampled, even when capturing a single axial plane is desired.

SR compressed imaging combined with optical sectioning would be a time- and computation-efficient strategy for 3D imaging or for rapidly imaging a single axial plane. Only the axially-restricted illumination pattern in the plane of interest would need to be known and probed for 2D image processing [33]. A 3D image can then be constructed from a series of 2D images, individually processed. We implemented this approach with focused light as it enables the use of two-photon excited fluorescence, which intrinsically provides optical sectioning. An axial-lateral projection shows that optical sectioning can indeed be achieved (figures 6(a) and (b)). Analysis of bead profiles by FWHM measurements indicated that an extended spatial resolution is attained along the lateral dimension ($L_\textrm{lat,SR} = 1.05~\pm~0.05$) and diffraction-limited performance along the axial dimension ($L_\textrm{axial,SR} = 3.91~\pm~0.11$) when compared to the diffraction-limited data ($L_\textrm{lat,DL} = 1.41~\pm~0.04$, $L_\textrm{axial, DL} = 4.39~\pm~0.02$) and reference object ($L_\textrm{lat,Obj} = 1.005~\pm~0.003$, $L_\textrm{axial, Obj} = 1.005~\pm~0.003$, figures 6(c) and (d)).

Zoom In Zoom Out Reset image size

In summary, fast volumetric imaging with extended lateral resolution can be achieved with two-photon excited fluorescence, focused light, a.nd SR CS algorithm.

The above analyses are valid for any optical system in which focusing is possible. In this section, we focus on the specific issue of imaging through a bent MMF. For imaging through a MMF, the wavefront entering the proximal end can be shaped such that a focus is formed in the sample [34–36]. Determination of the optimal wavefront requires a calibration, which is only valid for the mechanical state in which the calibration was performed, typically in a straight MMF. If the MMF is bent, the calibration is no longer valid and the focus quality decreases [25, 37, 38]. The same is true for speckle patterns; for a given input, the output will change as the MMF is bent.

The resilience to bending of SR compressed imaging through a MMF was evaluated for speckle and focus illuminations. Illumination changes were simulated with different bending radius $r_\textrm{bending}$ and the signal processed with the expected and known illumination patterns (that of a straight MMF). It was observed that using the BP algorithm with speckle illumination was less resilient to bending then using the same algorithm with focused light (figures 7(a) and (b)). Indeed, the reconstruction quality, as assessed by R with respect to the ground truth, dropped off at a much larger $r_\textrm{bending}$ for speckle than for focus illumination (figure 7(c)).

Zoom In Zoom Out Reset image size

Here again the sampling incoherence and error probability may play a role. However, the evolution of the illumination with bending is likely a more important factor. The correlation coefficient between the straight and bent illumination $R_\textrm{illum}$ was therefore calculated as a function of $r_\textrm{bending}$ (figure 7(d)). The illumination remained more similar with focus than speckle illumination during bending. This difference can be understood as follows: while the intensity distribution changes uniformly across the area illuminated with a speckle pattern with bending, the fraction of light contained in the focus is still much higher than the background outside the focus. $R_\textrm{illum}$ only partially captures this idea. In sum, the overall illumination changes less rapidly with bending for focused light than for speckles patterns and, as such, focus illumination offers more resilience to bending.