2.1.

Basic Principle of HiLo

In our original speckle-based HiLo algorithm,⁹ two raw images $I_s(\mathit{\rho })$ and $I_u(\mathit{\rho })$ are acquired with speckle and uniform illumination, respectively ($\mathit{\rho }$ is the 2D spatial coordinate at the camera plane). A final optically sectioned image is synthesized from the fusion of the two, where sectioning of low spatial frequency components is based on the contrast decay of the imaged speckle as a function of defocus. This contrast decay can be accelerated with the help of a wavelet-type filter.^9,28 Specifically, contrast can be defined as

$C_{\delta s}(\mathit{\rho })$ serves as a weighting function that preferentially extracts the optically-sectioned in-focus contributions in $I_u(\mathit{\rho })$. An estimate of these contributions is given by $I_{cu}(\mathit{\rho })=C_{\delta s}(\mathit{\rho })I_u(\mathit{\rho })$; however, this estimate is noisy due to residual speckle noise. The simplest strategy to control this noise is by the use of a low-pass filter, leading to a smoothed estimate of the in-focus contributions for low spatial frequencies, given by $I_{\mathrm{Lo}}(\mathit{\rho })=\mathrm{LP}[I_{cu}(\mathit{\rho })]$. The complementary high-frequency components, which are inherently optically-sectioned, are obtained by high-pass filtering the uniform image: $I_{\mathrm{Hi}}=\mathrm{HP}[I_u]$, where $\mathrm{HP}({\overrightarrow{\kappa }}_{\perp })=1-LP({\overrightarrow{\kappa }}_{\perp })$. The final HiLo image, optically sectioned over the full spatial frequency bandwidth, is synthesized by combining Hi and Lo components with a scaling factor $\eta$,

2.2.

HiLo with Non-local Means Denoising

2.2.1.

$I_{\mathrm{Lo}}$ with NLM denoising

The sectioning capacity of HiLo can be tuned by the cutoff frequency separating Hi and Lo. The larger the frequency, the tighter the optical sectioning, but also the more residual speckle noise becomes apparent in $I_{\mathrm{Lo}}$ because it is only weakly filtered from $I_{cu}$ [see Fig. 1(a)]. To aid this filtering, we incorporate the method of NLM denoising.¹³ The basic idea is to expand the regions over which filtering is performed without lowering the cutoff frequency. This is achieved by searching for similar patches throughout $I_{cu}$ (guided here by the uniform image) and performing a weighted average of these according to their similarity. The net result is that the residual speckle noise is substantially decreased while maintaining tight optical sectioning. Figure 1(a) illustrates this process using a numerically synthesized test target, where the out-of-focus background is simulated from the same target rotated by 180 deg. The algorithm takes two inputs, one $I_u$ from uniform illumination and the other $I_{cu}$ from the pre-processed difference image according to Eq. (1).

Fig. 1

Diagram of NLM-assisted HiLo. (a) For input images $I_u$ and $I_{cu}$, weighted averagings of similar non-local patches obtained from $I_u$ produce denoised ${\tilde{I}}_u$ and ${\tilde{I}}_{cu}$. The intermediate results are then high-pass/low-pass filtered respectively to form $I_{\mathrm{Hi}}$, $I_{\mathrm{Lo}}$, which are then fused to synthesize a final optically sectioned HiLo image $I_{\mathrm{HiLo}}$. (b) Comparison of pixel intensity distributions in the circled regions in $I_{cu}$ and ${\tilde{I}}_{cu}$.

First, similar patches are found in the uniform image $I_u$. The notion of similarity comes from treating a noisy patch as a realization of a random variable of given distribution (e.g., Poisson distribution in the case of shot noise), where the parameters of the distribution are determined by the underlying noise-free patch. A pair of noisy patches are considered similar when they obey the same distribution of common parameters.¹³ Consider a vectorized square patch $\mathit{x}$ containing ${(2N+1)}^2$ pixels centered on pixel $x$, where the center pixel has value $I_u(x)$ and the remaining pixels have value $I_u(x+\mathrm{d}x)$, with $\mathrm{d}x$ being the pixel index across the patch $\mathit{x}$ ($\mathrm{d}x\in [-2N(N+1),2N(N+1)]$). For a pair of patches $({\mathit{x}}_1,{\mathit{x}}_2)$ in the uniform image $I_u$ corrupted by Poisson noise, the patch similarity is calculated based on the log likelihood,

Eq. (3)

$$\mathcal{L}[I_u({\mathit{x}}_1),I_u({\mathit{x}}_2)]=\sum _{\mathrm{d}x=-2N(N+1)}^{2N(N+1)}\log {\mathcal{L}}_G[I_u(x_1+\mathrm{d}x),I_u(x_2+\mathrm{d}x)],$$

where ${\mathcal{L}}_G$ is a pixel-wise implementation of the generalized likelihood ratio,²⁹ which has been demonstrated to provide the best estimate in the case of strong noise levels and reasonably good estimates in the case of medium/low noise levels as compared to other metrics such as Euclidean distance,²⁹

Eq. (4)

$${\mathcal{L}}_G(I_1,I_2)=\frac{(I_1+I_2{)}^{I_1+I_2}}{2^{I_1+I_2}{I}_1^{I_1}{I}_2^{I_2}}.$$

Here, ${\mathcal{L}}_G$ equals 1 for identical pixel values and varies between $(\mathrm{0,1})$ for different pixels.

Following Eq. (3), the weights are calculated by introducing a filtering parameter $h$ that controls the degree of averaging between non-local patches, obtaining

Eq. (5)

$${\omega }_{{\mathit{x}}_1,{\mathit{x}}_2}=\exp \left\{\frac{\mathcal{L}[I_u({\mathit{x}}_1),I_u({\mathit{x}}_2)]}{h}\right\}.$$

For a pixel $x_1$ to be estimated, ${\omega }_{{\mathit{x}}_1,{\mathit{x}}_2}$ determines the contribution of pixel $x_2$ in the weighted average. The number of patches used for calculation is further defined by a large search window $\mathit{y}$ containing ${(2N_y+1)}^2$ pixels [blue rectangle in Fig. 1(a)].

Finally, the same pixel location is found in $I_{cu}$ and the weighted average of $I_{cu}$ is performed for all pixels within $\mathit{y}$,

Eq. (6)

$${\tilde{I}}_{cu}(x_1)=\frac{\sum _{x_2\in \mathit{y}}{\omega }_{{\mathit{x}}_1,{\mathit{x}}_2}{I}_{cu}(x_2)}{\sum _{x_2\in \mathit{y}}{\omega }_{{\mathit{x}}_1,{\mathit{x}}_2}}.$$

We observe that Eq. (6) is essentially an incoherent averaging of speckle realizations in $I_{cu}$, taking into account of sample variation, which is a common practice for suppressing speckle.³⁰ Figure 1(b) compares the intensity distribution in $I_{cu}$ and ${\tilde{I}}_{cu}$ before and after the application of NLM denoising in the circled region in Fig. 1(a) where the intensity is uniform. NLM denoising significantly reduces the standard deviation of the residual speckle contrast, thus improving the local smoothness of $I_{\mathrm{Lo}}$.

3.1.

Multiplane Imaging Setup

A schematic of the experimental setup is shown in Fig. 2. HiLo requires fast switching between speckle and uniform illumination in an otherwise conventional widefield microscope. This was obtained here with two separated beams combined by a polarizing beam splitter (PBS). The first beam was produced by a collimated blue light-emitting diode (LED) (SOLIS-470C); the second beam was from a continuous wave laser (Vortran Stradus, 488 nm) incident on a static diffuser conjugate to the back pupil of the microscope objective (Olympus UMPLFLN 20XW or LUMPLFLN 40XW). Both light sources were triggered by a data acquisition device (DAQ) (National Instruments, NI USB-6356) to produce alternating illumination synchronized to the camera exposure (PCO.edge 4.2). We note that in previous single-beam designs, the uniform illumination image was obtained by rapidly randomizing the speckle illumination patterns with a tilting or rotating diffuser.^8,31 Our current two-beam design avoids two potential drawbacks. First, the toggling between speckle and uniform illumination is no longer slowed down by any need to average over multiple speckle patterns. Second, the speckle illumination pattern in our case remains fixed from frame to frame, meaning that it does not introduce temporal fluctuations from frame to frame that could perturb measurements of relative changes in fluorescence, as utilized, for example, with calcium imaging. Finally, to perform simultaneous multiplane imaging, the emitted fluorescence was directed into a nine-plane z-splitter prism²² and projected onto the single camera sensor with a relay lens pair (demagnification 250/90).

Fig. 2

Experimental setup. Light generated from LED and laser is recombined with a polarizing beamsplitter (PBS), providing uniform and speckle illumination of the sample. A z-splitter prism (here, shown as a three-plane prism) enables simultaneously imaging of multiple depths with a single camera.

3.2.

Sample Preparation

In-vivo imaging was performed of both the heart (isl2b:Gal4 UAS:Dendra, 8dpf) and brain (elavl3:H2B-GCaMP6f, 6dpf) of transgenic zebrafish larvae. To prepare for imaging, larvae were embedded in 1.5% low melting point agarose gelled on a petri dish. Larvae were positioned right-side down for heart imaging and ventral-side down for brain imaging. The solidified agarose and fish were immersed in fish water before being transferred onto the sample stage for imaging.

3.3.

Image Processing

The wavelet filter is defined in frequency space as

The default value of ${\sigma }_w$ was set such that the width of the wavelet filter was about one speckle-grain size. In the case of low signal level, the wavelet size was set to slightly larger than the speckle-grain size to improve SNR. Bias introduced by shot noise and readout noise was subtracted from the measured variance,⁹ presuming a camera gain of $0.46{\mathrm{e}}^{-}/\mathrm{ADU}$ and readout noise of $1.2{\mathrm{e}}^{-}$. The scaling factor $\eta$ was determined either theoretically from system parameters⁹ or empirically based on intensity distributions.

For the implementation of NLM, the patch size was chosen to be $3\times 3\text{pixels}$ in all experiments, corresponding to about one speckle grain per patch, and the search window size was chosen in the range $15\times 15$ to $61\times 61\text{pixels}$ depending on the speckle contrast. A larger search window allows the averaging of more patches (i.e., more uncorrelated speckle) but also takes longer to process. The filtering parameters $h$ for $I_{\mathrm{Lo}}$ and $h^{\prime }$ for $I_{\mathrm{Hi}}$ were empirically set to be around 5 to 20 and 0.4 to 0.6 to make allowances for different speckle and noise variances. The NLM implementation was performed with a compiled $C^{++}$ mex-function in parallel. Other pre/post-processing was performed using MATLAB 2021a.

To characterize the quality of HiLo reconstructions, we made use of the peak SNR (PSNR), defined as

4.1.

Simulation Results

We first simulated widefield imaging of a pollen grain under uniform and speckle illumination. The ground truth was obtained from a confocal image stack of an actual pollen grain obtained with Olympus FV3000 (60×, 1.2 NA). To simulate widefield images, the volume ($323\times 323\times 35\text{pixels}$, $42\times 42\times 19.6\mu {\mathrm{m}}^3$) was convolved with a simulated 3D point spread function and the defocused planes produced the background for each focused slice. Shot noise and readout noise were introduced assuming a camera gain of 2.5 [$\mathrm{ADU}/{\mathrm{e}}^{-}$] and readout noise of $1.6{\mathrm{e}}^{-}$ (rms).

Figure 3 shows the simulated uniform- [Figs. 3(a) and 3(f)] and speckle- [Figs. 3(b) and 3(g)] illumination of raw images at two different depths, where the in-focus information is corrupted by a strong out-of-focus background. HiLo was performed with a wavelet filter size equal to one speckle grain such that the sectioning is approximately one longitudinal speckle length. Our basic HiLo algorithm effectively suppressed background, as expected, but also manifestly suffered from residual low-frequency speckle noise [Figs. 3(c) and 3(h)]. In comparison, NLM HiLo benefited from the incoherent averaging of many more speckle grains, leading to largely denoised images that still preserved the intensity variations of the object itself [Figs. 3(d) and 3(i)], as evidenced by the higher structural similarity index measures (SSIMs). As a result, the fine details of the pollen grain are lost in the maximum-intensity projection (MIP) of the basic HiLo reconstruction [Fig. 3(l)], but they are well preserved in the MIP of the NLM denoised reconstruction [Fig. 3(m)], as verified by comparing with the ground-truth confocal MIP [Fig. 3(n)].

Fig. 3

Simulation of a pollen grain stack. Images of two different depths are selected for presentation. Uniform and speckle illuminated images at $z_1$ (a,b) and $z_2$ (f,g). Corresponding optically-sectioned reconstructions at two depths using basic HiLo (c,h) and NLM HiLo (d,i). (e,j) Confocal images serve as a ground truth. (l)–(n) MIPs and associated SSIMs of image stacks using basic (l) and NLM HiLo (m) compared to ground truth (n). (k) PSNR for different HiLo reconstructions as a function of signal strength. Solid line and shaded area indicate mean and standard deviation of 10 trials. Scale bar: $10\mu \mathrm{m}$.

We further compared the reconstruction performance for different signal levels in Fig. 3(k), where PSNR was used to quantify the similarities between the reconstructions and ground truth. Basic HiLo was compared with NLM HiLo with denoising of $I_{\mathrm{Lo}}$ only (i), denoising of $I_{\mathrm{Hi}}$ only (ii), and denoising of both (iii). At each signal level, the PSNR of the full volume was calculated using the confocal stack as a reference. The resulting PSNR was plotted as a function of the maximum number of photons per pixel in the uniform stack. At low SNR, where shot-noise variance is stronger than speckle variance, denoising of $I_{\mathrm{Hi}}$ plays a more important role than denoising of $I_{\mathrm{Lo}}$, whereas the roles are reversed at high SNR, where speckle noise becomes dominant.

4.2.

Experimental Results with Fixed Slide

Having demonstrated NLM HiLo with simulated data, we next evaluated its performance using our multiplane imaging setup with fixed samples. We imaged a Melittobia digitata slide (Carolina Biological) with a 20×, 0.5 NA objective and compared the multiplane reconstructions from basic HiLo and NLM HiLo with a confocal stack taken of the same volume. The wavelet parameter was set to be ${\sigma }_w=1.5$, meaning that the wavelet filter was slightly larger than a speckle grain in this case. Figures 4(a) and 4(b) show the raw speckle and uniform illumination images of a single plane. The high-frequency images derived from $I_u$ without and with NLM denoising are shown in Figs. 4(c) and 4(f), confirming that NLM effectively attenuates noise while largely preserving image resolution. The low-frequency images derived from $I_{cu}$ and denoised ${\tilde{I}}_{cu}$ are shown in Figs. 4(d) and 4(g). The basic $I_{\mathrm{Lo}}$ was corrupted by strong speckle noise that ultimately propagated to the final HiLo image [Fig. 4(e)], overwhelming the object structure. In comparison, the speckle noise in $I_{\mathrm{Lo}}$ was significantly attenuated by NLM, allowing the object structure to be revealed [Figs. 4(g) and 4(h)]. We also compared multiplane HiLo reconstructions in Figs. 4(i) and 4(j), where the nine-plane projections were color-coded in depth. While both methods achieved optical sectioning comparable to confocal microscopy [Fig. 4(k)], basic HiLo suffered from deleterious speckle noise that became largely attenuated with NLM HiLo.

Fig. 4

Experimental result. Multiplane HiLo imaging of a fixed slide. Widefield imaging of single plane with speckle (a) and uniform (b) illumination. Intermediate $I_{\mathrm{Hi}}$ (c), $I_{\mathrm{Lo}}$ (d), and composite $I_{\mathrm{HiLo}}$ (e) obtained with basic HiLo algorithm. Intermediate denoised $I_{\mathrm{Hi}}$ (f) and $I_{\mathrm{Lo}}$ (g) and composite $I_{\mathrm{HiLo}}$ (h) obtained with NLM HiLo. (i,j) Depth-coded projections across nine planes from basic and NLM HiLo. (k) Depth-coded projections across a confocal stack. Scale bar: $100\mu \mathrm{m}$.

4.3.

Fast Imaging of Beating Zebrafish Heart

The zebrafish larva is a popular animal model in cardiology.³² Real-time contraction measurements provide important information about cardiac function and enable correlation studies with ${\mathrm{Ca}}^{2+}$ transients in the heart.³³ Such measurements benefit from fast, volumetric and high contrast imaging, as enabled here by multiplane HiLo microscopy which provides both optical sectioning and depth information.

A beating zebrafish heart was imaged with nine-plane HiLo at a frame rate of 41 Hz (exposure time 1 ms; total recording time 17 s), using a 40×, 0.8 NA objective that resulted in an axial separation of $5.5\mu \mathrm{m}$ between planes. The results are shown in Fig. 5. Stacks of raw uniform-illumination images and NLM HiLo reconstructions are shown in the left and middle panels of Fig. 5(a). The rightmost panel of Fig. 5(a) shows a zoomed-in region comparing uniform illumination (top), basic HiLo (middle), and NLM HiLo reconstructions (bottom). NLM HiLo effectively attenuates local non-uniformities arising from speckle as well as additional noise from detection (shot and readout noise).

Fig. 5

In-vivo imaging of beating larval zebrafish heart. (a) Multiplane stack of uniform illumination (left) and NLM HiLo reconstruction (middle). Comparison of regions outlined by rectangle for uniform illumination (top right), basic HiLo (middle right) and NLM HiLo (bottom right) images. (b) Representative snapshots of projections color-coded in depth showing cardiac contraction. (c) Kymogram obtained from line profile [dashed line in panel (b)] as a function of time. (d) Ventricle diameter as a function of time. Scale bar: $50\mu \mathrm{m}$.

Representative frames showing depth color-coded projections at different times are displayed in Fig. 5(b), where the systole and diastole are easily distinguished by virtue of the fast imaging speed. Each cardiac cycle was sampled by at least 20 images, enabling analysis of cardiac function. Here, a kymogram was determined by tracking the fluorescence signal through the ventricle [dashed line in Fig. 5(b)] and displayed in Fig. 5(c). The local ventricle diameter was calculated based on position differences of the external wall³³ and plotted as a function of time [Fig. 5(d)], illustrating different beating profiles throughout the ventricle. The heart rate was 1.8 Hz on average.

4.4.

In-Vivo Calcium Imaging of Zebrafish Brains

Finally, we illustrate the capabilities of multiplane HiLo when applied to imaging of neuronal activity in zebrafish larvae, which normally can be resolved only with confocal or two-photon. Three-plane widefield imaging of brain dynamics with optical sectioning has been demonstrated with HiLo and an electrically tunable lens.³⁴ Here we make use of a passive z-splitter prism to simultaneously image nine planes with $22\mu \mathrm{m}$ interplane separation. Zebrafish larvae expressing nucleus-localized GCaMP6f were imaged using a 20×, 0.5 NA objective with 40 ms exposure time and 10 Hz frame rate. We recorded for a total duration of 73 seconds. The raw uniform and speckle illumination images were processed with NLM HiLo and the resulting video was post-processed using CaImAn³⁵ for calcium dynamics analysis. The results are shown in Fig. 6.

Fig. 6

In-vivo multiplane calcium imaging of larval zebrafish brain. (a)–(c) Top: temporal MIP of HiLo frames at three different depths. Bottom: enlarged regions delimited by white rectangles showing neurons (left) and their corresponding calcium traces (right). (d) Comparison of uniform illumination (left) and HiLo (right) images for single neuron in panels (a)–(c) when active. (e) SBR for uniform illumination and HiLo images at different planes. (f) Number of identified neurons at different planes. (g) Depth color-coded neuron centroids overlaid on a HiLo MIP image. Scale bar: top panel in (a)–(c),(g): $100\mu \mathrm{m}$; elsewhere $10\mu \mathrm{m}$.

The top panels in Figs. 6(a)–6(c) display temporal MIPs obtained over the full recording duration, showing all active neurons from three central planes. The activity of individual neurons can be resolved in both space and time [Figs. 6(a)–6(c), bottom]. Single-frame widefield and HiLo images obtained from the zoomed-in rectangle regions at three depths are compared in Fig. 6(d). For the uniform illumination image, the in focus neuronal signal was easily overwhelmed by the strong background, resulting in a low signal-to-background ratio (SBR) that hampers the accuracy of neuronal segmentation. In comparison, NLM HiLo images feature lower noise and much higher contrast, by virtue of being optically sectioned. The SBR comparison of uniform illumination and HiLo images at each plane is shown in Fig. 6(e), where SBR was calculated from the ratio of the signal intensity within each neuron while active and the average background intensity within the associated circular surround region (radius $40\mu \mathrm{m}$). HiLo increased SBR by nearly an order of magnitude owing to the rejection of background.

Figure 6(f) shows the total number of active neurons identified per plane over the full recording duration. Most of the neurons were identified in the three central planes ($z=4$ to 6) partly owing to the inherent distribution of neurons and partly owing to the increasing difficulty to resolve neurons at larger depths due to scattering. The identified neuron centroids, color-coded in depth, are shown in Fig. 6(g).