Correcting anisotropic intensity in light sheet images using dehazing and image morphology

Light-sheet fluorescence microscopy (LSFM) provides access to multi-dimensional and multi-scale in vivo imaging of animal models with highly coherent volumetric reconstruction of the tissue morphology, via a focused laser light sheet. The orthogonal illumination and detection LSFM pathways account for minimal photobleaching and deep tissue optical sectioning through different perspective views. Although rotation of the sample and deep tissue scanning constitutes major advantages of LSFM, images may suffer from intrinsic problems within the modality, such as light mismatch of refractive indices between the sample and mounting media and varying quantum efficiency across different depths. To overcome these challenges, we hereby introduce an illumination correction technique integrated with depth detail amelioration to achieve symmetric contrast in large field-of-view images acquired using a low power objective lens. Due to an increase in angular dispersion of emitted light flux with the depth, we combined the dehazing algorithm with morphological operations to enhance poorly separated overlapping structures with subdued intensity. The proposed method was tested on different LSFM modalities to illustrate its applicability on correcting anisotropic illumination affecting the volumetric reconstruction of the fluorescently tagged region of interest.


INTRODUCTION
For understanding dynamic processes taking place at the cellular and molecular level, it is necessary to perform spatiotemporal volumetric analysis. 1 Among non-invasive imaging modalities, light-sheet fluorescence microscopy (LSFM) is emerging as a powerful tool for studying developmental biology due to its highly advantageous capabilities such as short pixel dwell time while still being able to capture a high dynamic range. 2 However, LSFM images with Gaussian illumination may suffer from limitations in the form of stripe effects and nonhomogenous fluorescence intensity due to photon refraction through a heterogeneous scattering medium. 3 Hardware solutions through dual-side and multidirectional selective plane illumination microscopy (mSPIM) have assisted to overcome these issues in conjunction with resonant mirrors. 4 The advent of tissue clearing techniques has also helped to minimize tissue scattering by removing light scattering endogenous pigments and lipid membranes that may otherwise affect the laminar light sheet and cause tissue opacity. 5 In addition, the oblique scanning method facilitated a major advancement in image sub-voxel resolution using a low-power objective lens. 6 Although myriad solutions have been developed, it remains a challenge to obtain an impeccable distinction between the background and foreground due to interference from unfocused fluorescence in the axial direction. Primarily, refractive index mismatch within the tissue and mismatch between the tissue with the embedding medium cause photon scatter outside the focal volume. 7 This leads to sample blurriness axially and laterally, 8 hindering the accuracy of volumetric reconstruction, and is further exacerbated for a low numerical aperture (NA) objective lens that has a smaller solid angle limiting photon capture. 9 As a result, underwater imaging requires pre-processing to maintain spatial homogeneity in the reconstructed volume.
Anisotropic intensity correction through intensity distribution approximation, based on scattering propagation models, has been implemented using conventional techniques such as histogram modeling. 10 Such methods effectively increase the intensity dynamic range of an attenuated image, but they are susceptible to saturating bright pixel neighborhoods where evanescent scattered light flux is on the order of 10-100 greater than the incident light field. 11,12 Biomedical applications require the critical object in the region of interest (ROI) to be aberration free and unaffected by the ambient medium. The dehazing method based on the dark channel prior (DCP) algorithm is an illumination correction model developed for correcting haze affected images captured in an outdoor environment. 13 In this paper, we applied this method to LSFM images by assuming that the light propagation model is based on the underwater medium instead of air. [14][15][16][17] Previous studies have indicated that the scattering coefficient of mounting media (water) and atmospheric light propagation is a function of wavelength. 18 For the underwater condition, scattering depends on the dissolved particulate cross sectional area, 19 whereas the composition of molecular gases determines scattering in the air. 18 Moreover, the dehazing method possesses a high sensitivity to noise as the algorithm corrects attenuation in the aberrated image induced by medium turbidity along the line-of-sight (LOS) propagation. 13 Translation of the mechanical stage along the depth-of-view (DOV) and lens aberration may result in phase delay or a focal point shift along the optical axis, 7 influencing the diffraction-limited point spread function (PSF) resolution. By assuming depth invariancy, 20 we have attempted to apply the dehazing method to recover the asymmetricity of the PSF intensity distribution. Using this method, we were able to resolve sub-pixel details along DOV (axial domain) and fieldof-view (FOV) (lateral domain) without employing a high numerical aperture (NA) lens or using a speckle-based near field aperture projection of the PSF. 21 The use of morphological operations based on isotropic structuring elements has assisted edge restoration in uneven contrast images in imaging modalities, such as x ray and CT (computed topography), which may be affected by different tissue absorption coefficients. 22,23 By utilizing image reconstruction operations such as opening and closing, aberrations affecting the edges such as discontinuities or gaps can be corrected to maintain a uniform contour. 24 The size of the structuring element is the factor that determines the bandpass in the frequency domain and the intensity normalization. 24

Implementation of LSFM
Our tunable LSFM is composed of dual-side illumination pathways integrated with multi-dimensional image fusion [ Fig. 1(a)]. In the illumination path, a cylindrical lens coupled with a 4Â objective lens (4X Plan Apochromat Plan N, Olympus, Tokyo, Japan) generated a collimated Gaussian light-sheet. From the calculation of the Rayleigh length (half of the confocal range), we reduced the mechanical slit to 1.99 mm to generate a 15 mm confocal range to cover the circumferential length of 4 days post-fertilization (dpf) zebrafish. To eliminate stripe effects, we added a galvanometer mirror and optimized the system at a frequency of 2000 Hz [Figs. 1(b) and S1]. When the sample was scanned through the confocal range of the light sheet along the z-direction, the sCMOS camera (ORCA flash 4.0, Hamamatsu, Japan), located at the end of the detection path coupled with a tube lens, recorded a stack of 2D plane images along different z depths.
An attachable oblique scanning adapter could be easily added for applying the voxel super-resolution technique to provide high spatial resolution from a low power objective [ Fig. 1(c)].
Anisotropic illumination correction by integrating the dehazing algorithm and background subtraction However, dehazing based on DCP does not consider endogenous fluorescence in the focal plane due to forward scatter. Due to refractive index mismatch, light scatter takes place through the specimen in a direction perpendicular to the camera. Therefore, the dehazing method is prone to restoring unrequired background autofluorescence Correcting PSF intensity distribution using DCP An improvement in the axial and lateral resolution was observed in the attenuated PSF after dehazing with a shift in the central peak intensity (Fig. S5). We concluded that the skew for the lateral and axial resolution graphs reflects the anisotropic shape of the bead in the xy domain as is reflected on the intensity profile plots (Fig. S6), along with depth variance along the yz domain. As deconvolution is a nonlinear process, we observed salt and pepper noise in the deconvolved images for some cases. To remove any effect of random noise introduced by deconvolution, an edge preserving bilateral smoothening filter was applied. Deconvolution of the background subtracted image using the dehazed PSF shows promising results for symmetric restoration of objects in the focal plane [Figs. 3, 6(e), and S3(c)]. After obtaining the deconvolved image, we applied the top-hat transform to the dehazed image. The top-hat transform is based on image erosion, which is utilized to remove weakly connected regions of interest (ROI) and also image dilation to emphasize and fill-in image boundaries, in that order. Top-hat transform assists by removing any pixel discontinuities induced by background subtraction and effectively isolates foreground objects in the image plane. The bottomhat transform is based on image closing to reinforce depth details, which are blurred out and out-of-focus due to autofluorescence affecting the axial resolution. The final image was reconstructed by adding the image subtraction of the top-hat and bottom-hat transform results to the dehazed image. Isotropic structural integrity was achieved across

Application of the proposed method to various LSFM techniques
We acquired images using different LSFM modalities such as single/dual illumination SPIM, multiview SPIM with dual illumination, and the oblique scanning super-resolution stage integrated with dual illumination. To construct the dual illumination arms orthogonal to the detection camera, a 50:50 beam splitter was used to separate the incoming Gaussian laser spot into two equally intense laser beams For resolving the vasculature of the zebrafish at different scales without loss of spatial homogeneity, we acquired images using an oblique scanning stage to apply the voxel super-resolution algorithm. By applying our enhancement algorithm with the image processing pipeline, minute structures, such as optic arteries and veins on the head, and intersegmental vessels on the tail were reconstructed clearly with a 4Â objective lens [ Fig. 6(b)]. Taking advantage of our processing, we reconstructed a multi-view fused stack from the processed images that were captured by an oblique scanning stage. Multi-view fusion is utilized to overcome opacity or translucency, which affects spatial homogeneity due to light scattering, by adding images from different angles (Fig. S9). The fused image exhibits superior depth reconstruction of different orthogonal perspectives in conjunction with the oblique scanning technique (Fig. S10).

Evaluation parameter of image restoration
For restoration quality assessment, the structural similarity index (SSIM) full reference method was used by using a histogram adjusted raw image as reference [ Fig. S11(d)]. The lower/upper 10% of intensity distribution was saturated in order to enhance under exposed structures. Using the reference image, we compared the loss/degradation by our proposed method using dehazing followed by background  Adjusting parameters for either algorithms yielded similar global SSIM values with the lowest score being þ0.1 and the best score being 0.9. We also calculated the peak signal to noise ratio (PSNR) to understand improvement in terms of noise reduction during each step of the processing. For PSNR evaluation, an oblique scanned raw image was used for comparison (Fig. S12).

DISCUSSION
In this paper, we demonstrated the applicability of the DCP algorithm for restoring isotropic intensity in the FOV and morphological operations for correcting diffraction-limited volumes with nonhomogenous spatial homogeneity. We observed some intrinsic limitations in the single/dual side LSFM images acquired at the detection focal plane that affected isotropic reconstruction of the z stack. The images composed of local pixel patches were innately high in photon count with reduced visibility [ Fig. 6(a)]. For this study, we assumed that scattered wavelengths in these regions are longer than the excitation wavelengths as seen in Raman Stokes scattering and are nonpropagating. 11 This light scatter spreads laterally across the tissue but  APL Bioengineering ARTICLE scitation.org/journal/apb does not participate in the far field image formation. It decays evanescently in a direction perpendicular to the surface with the intensity distribution being amplified by an order of 10-100 in the evanescent field region. 11 Coupled with a laser intensity beyond a certain threshold, 7 it can cause optical distortions for any particulates that are less than half the input wavelength. As a result, out of focus tissue planes in the evanescent light field are illuminated. Due to this limitation, dull/low light patches with relatively uniform colocalized zero intensity crossings appear saturated due to over estimation of haze thickness (Figs. S13 and S14). The capability of an interrogating light sheet to illuminate the sample and the ability of the emitted light to escape the sample and reach the detection camera constitute important factors for determining image resolution. Using a low NA focusing objective can compromise the excitation efficiency of the focused beam waist 7,18 (Fig. S15). As LSFM involves light propagation through media with differing refractive indices, it is also necessary to develop or select intensity correction algorithms that account for these changes in optical properties of the media. Even when refractive index mismatch between the FEP tube and mounting media is neglected, light scattering occurs due to the variation in tissue refractive indices. [25][26][27][28][29] Tissue birefringence due to the tissue refractive index is a function of the density, resulting in a gradient refractive index. 30 Tissue transparency in animals is a function of surface irregularities at the cellular/molecular level, which affect the potential of tissue to reflect light. Using the Weber formula for defining contrast (Cd) of a transparent animal, we get 28 where Ld is the radiance of the object viewed at a distance in the underwater medium and Lb is the background radiance, 28 Ld ¼ ðLoe ðÀcdÞ þ Lbð1 À e ðÀcdÞ ÞÞ; where Lo is the inherent radiance at zero viewing distance, c is the beam attenuation coefficient (wavelength is a function of the depth), and d is the distance of the object to the viewer. Here, the detected object radiance (Ld) depends on object irradiance at the point of focus and the effect of atmospheric light scatter along the LOS. As we are using along working distance air lens with low NA, we can infer from the above equation that as the viewing distance is increased, inherent object radiance will decrease and the scattered light (Lb) will dominate. For a transparent animal, the inherent radiance consists of two parts, namely, the background radiance that is transmitted through the animal and the light scattered toward the viewer; hence, 28 Ld ¼ ðLbTe ðÀcdÞ þ LdSe ðÀcdÞ þ Lbð1 À e ðÀcdÞ ÞÞ; (3) where T represents transmission through the animal and S is a fraction of the environmental light scattered to the observer. Comparing equation (3) with the Dehazing equation, 13,14 we get where I(x) is the attenuated image recovered using the camera (Fig.  S15), J(x) is the original scene irradiance that has a multiplicative effect with the transmission distance, t(x) (Fig. S15), and A is the ambient/ background light. Therefore, we infer that the DCP method does not take into account the background radiance, Lb, which is transmitted through the illuminated portion, T, of the object. To counter this problem, we performed the background subtraction operation on the dehazed image to remove innately present autofluorescence in the surrounding tissue. Fluorescent images captured at higher frame rates suffer from poor distinction between the fluorescently tagged marker and background noise. To further elucidate depth details and separate overlapping edges, we used top-hat and bottom-hat morphological transforms based on flat structuring elements. As morphological operations are non-linear operations based on spatial ordering of the pixel coordinates and not the pixel intensity values, these methods can be applied efficaciously to grayscale turbid media images. To maintain isotropic kernel convolution values in all directions, we used a sphereshaped structuring element that was greater than the smallest regionof-interest (ROI) for all the morphological operations.
As the PSF of diffraction-limited microscopy is depth and space invariant with respect to the detection optic axis, we applied the dehazing algorithm as a PSF intensity correction model. An intensity correction of the recorded emission intensity for a fixed dipole, or fluorescent bead, was implemented by assuming that the intensity will vary according to the random orientation of the bead (Fig. S9). Assuming anisotropic photon collection efficiency, the intensity will experience a higher degree of exponential decay for long working distance objective lenses.
In summary, with the help of our image processing, we were able to overcome intrinsic LSFM limitations and improve the spatial resolution of the turbid media image. The dehazing model was used to enhance subdued intensity peaks and thus restore attenuation caused due to propagation through an opaque medium or light sheet skew. Using morphological operators, we were able to localize peaks and valleys in the local image structure in order to resolve diffraction limited structures affected by light scatter or stripe artifacts. At the same time, we were able to capitalize on the advantages provided by the LSFM modality, namely, the wide FOV provided by the low NA infinity corrected objective lens, and acquire zebrafish z-stacks with short acquisition time. This research will benefit in vivo developmental biology studies, visualized using fluorescence modalities that require accurate volumetric reconstruction for further analysis.

Preparation of zebrafish imaging
The experiments were performed in compliance with the approval from the UT Arlington Institutional Animal Care and Use Committee (IACUC) protocol (#A17.014). The transgenic lines, Tg(flk:gfp), were raised at the UT Arlington Aquatic Animal Core Facility. The flk promoter-driven green fluorescent protein (GFP) from the Tg(flk:gfp) zebrafish line was expressed in vascular endothelial and endocardial cells. 31 To maintain the transparency of zebrafish embryos, the medium was supplemented with 0.0025% phenylthiourea (PTU) to suppress pigmentation at 20 h post fertilization (hpf). At 4 dpf, live transgenic zebrafish embryos were anesthetized in 0.05% tricaine and immersed in 37 C, 0.5%, low-melt agarose for imaging. Prior to agarose solidification, the embryos with low-melt agarose were transferred to a fluorinated ethylene propylene (FEP) tube to match the refractive index (1.34) to that of water (1.33).

LSFM setup
Our tunable LSFM was set up to be able to use single side illumination, dual-side illumination, multi-view fusion, and oblique scanning for the VSR (Voxel super resolution) technique. To cover the entire width of the zebrafish, we calculated the confocal range from double of the Rayleigh length (z R ), which is where k is the wavelength of our laser (473 nm) to illuminate the green fluorescent protein in the transgenic zebrafish line and w 0 is the beam waist, which is the thinnest focal point of the Gaussian beam. The width of the light-sheet wðzÞ along the z-axis can be calculated as follows: We also adjusted the slit width to control the light-sheet thickness at the sides of the sample of ffiffi ffi 2 p wðzÞ, which is the thickness of the lightsheet under the Rayleigh length.

Voxel super-resolution
The technique of super-resolution by using oblique scanning follows the principle of the function of pixel super-resolution techniques in a 2D spatial domain. Super-resolved estimate I is most consistent with multiple measurements P k after a series of degradation operators. 32 Here, we solved I via minimizing the following cost function: where q is the difference between the model and measurements. S k is the geometric motion operator between HR (high-resolution) estimate I and k th LR (low-resolution) P k , the point-spread-function of the oblique scanning LSFM system is modeled by the blur operator O, and D k represents the decimation operator that models digital sampling of the camera. Theoretically, the computation estimates a unique HR image, which has maximum likelihood to the LR inputs after given degradations S k , O k , and D k are applied. Practically, a steepest descent method is provided to iteratively approach a converged super-resolved solution at high efficiency, where S T k ; O T k ; D T k represent the inverse operation of S k ; O k ; D k , respectively. This maximum likelihood estimation (MLE) solving process is further described in the form of a block diagram. A detailed description can be found in our previous research.

Image processing pipeline
The proposed method was implemented and tested in MATLAB R2019B and ImageJ. Observing better lateral resolution as compared to axial resolution, the processing was applied to every individual image slice in the z-stack (Fig. S16). The dehazing, deconvolution, and morphological operations were performed using the imreducehaze (dehazing operation), deconvlucy (point spread function deconvolution), and imtophat and imbothat (morphological operations) functions, respectively. We performed the background subtraction method in Fiji and used Amira to reconstruct the z-stack by specifying the voxel size parameters.
The dehazing method tends to amplify the autofluorescence of tissue outside the focal plane and as a result hampers the intensity restoration [ Fig. 6(c)]. To avoid this artifact, we performed background subtraction on the dehazed image [ Fig. 6(d)]. Furthermore, we deconvolved the background subtracted image using a dehazed point spread function [Figs. 6(e), 6(e 0 ), and 6(e 00 )]. The deconvolution result is a function of the number of iterations and hence can introduce random noise. We smoothened the image after deconvolution using a Gaussian smoothening filter in MATLAB using the imgaussfilt function or alternatively an edge preserving bilateral smoothening filter in ImageJ (Plugins > Process > Bilateral Filter). Both methods produced similar results without any significant change. It is important to adjust the parameters of PSF dehazing and select the PSF along LOS propagation unaffected by spherical aberrations such that it does not lead to PSF oversampling causing unwanted saturation. The final image [Figs. 6(b), 6(b 0 ), and 6(b 000 )] was obtained through an image arithmetic operation as follows (Fig. S16): (background subtracted dehazed image) -[tophat(background subtracted dehazed image) þ bottomhat(dehazed image deconvolved with the dehazed point spread function)].
After reconstructing the z-stack by specifying the required voxel arguments, we were able to visualize the 3D local morphology across different scales by adjusting the transparency of the volrenGreen color map and the histogram to saturate critical objects in the 3D FOV.
(1) Dehazing: The raw images were contrast-corrected using the dehazing method based on the dark channel prior (DCP) algorithm. Image formation achieved by the dehazing algorithm is based on the following model: 32 where I(x) is the attenuated image, J(x) is the dehazed image at the original scene of irradiance captured using the camera over LOS propagation, A is the ambient light, and t(x) is the light transmission that is unaffected by interference by the ambient surroundings.
(2) Background subtraction: The sliding paraboloid background subtraction filter in Imagej based on the rolling ball averaging algorithm 33 was used to remove the effect of tissue autofluorescence and lateral light scatter. It is based on image subtraction of an averaged background value for each pixel from the original image by sliding a paraboloid estimated by four parabolas in four directions: x,y and two 45 directions. The recovered image contains a gradient refractive index of overlapping tissue represented as an unsigned integer in the 2D image vector. These overlapping spatial variations represented as intensity variations in local pixel neighborhoods can be effectively corrected using background subtraction as it is based on the assumption that each pixel intensity on the XY spatial image plane can be imagined as a third dimension with respect to the intensity values of other pixels in the local neighborhood. This is analogous to using non-flat structural elements for the mentioned morphological operations. (3) Gaussian blurring filter: The deconvolution process may result in over-emphasized intensity peaks or unwanted noise amplification. To reduce this effect, we applied the following Gaussian smoothing function: where x denotes the distance from the horizontal axis, y is the distance from the vertical axis, and r (sigma) is the standard deviation for a 2D isotropic Gaussian convolution kernel. We specifically used the Gaussian smoothing filter, as other filters that replace the pixel values with the weighted average of neighboring pixel values may pass intensity values that may be beyond the passband (ringing oscillations). This is also because convolution with a smoothing kernel in the spatial image domain is analogous to multiplication in the frequency domain with a low pass frequency filter. Since the value of sigma (rÞ determines the degree of smoothing, we used a smaller r value to maintain a tighter roll-off for the Gaussian function and to avoid excess blurring. The top hat transform is based on image subtraction of the original image by the opened image. This morphological operation was utilized to restore contrast in uneven illumination along the line of sight (LOS) propagation in this case. Image opening can be defined as image erosion followed by image dilation. The bottom hat transform was used to emphasize poorly illuminated depth details, and the bottom-hat image was obtained by dilating the image and eroding it. The top hat transform was employed to maintain the size of larger objects and remove smaller discontinuities, whereas the bottom hat transform was used to fill-in smaller regions. For images that were relatively unaffected by any intensity aberrations, we used the opening operation as an alternative to the top hat operation and closing operation instead of bottom hat filtering.

SUPPLEMENTARY MATERIAL
See the supplementary material for understanding the image restoration methods in the implemented pipeline.