Separation of ballistic and diffusive fluorescence photons in confocal Light-Sheet Microscopy of Arabidopsis roots

Image quality in light-sheet fluorescence microscopy is strongly affected by the shape of the illuminating laser beam inside embryos, plants or tissue. While the phase of Gaussian or Bessel beams propagating through thousands of cells can be partly controlled holographically, the propagation of fluorescence light to the detector is difficult to control. With each scatter process a fluorescence photon loses information necessary for the image generation. Using Arabidopsis root tips we demonstrate that ballistic and diffusive fluorescence photons can be separated by analyzing the image spectra in each plane without a priori knowledge. We introduce a theoretical model allowing to extract typical scattering parameters of the biological material. This allows to attenuate image contributions from diffusive photons and to amplify the relevant image contributions from ballistic photons through a depth dependent deconvolution. In consequence, image contrast and resolution are significantly increased and scattering artefacts are minimized especially for Bessel beams with confocal line detection.

A fundamental principle in light microscopy is to uncover the composition of matter by sending photons with defined properties onto a piece of matter and detecting the scattered photons. However, things can become very complicated when the investigated object is large and of complex biological structures and multiple light scattering occurs. Photons that have been scattered too often lose their directional information, i.e. photons are no longer ballistic, but become diffusive. This loss of information does not only occur on the illumination side, but also on the detection side. When arriving on the camera, the origin of the photons is often unknown 1 . In consequence, an image becomes blurred and noisy.
A suitable microscopical method to tackle the problem of light scattering in both the illumination and detection path is light-sheet based microscopy (LSM) 2,3 . Besides advantages such as high acquisition speed, effective sectioning, high contrast and low phototoxicity, this technique is also fascinating because it allows to observe how the illumination light propagates from the side through selected planes of the object. Furthermore, the influence of the detection depth can be observed quite well. In this way, the analysis of both scattered laser light and fluorescence light helps to better understand the formation of useful image data and unwanted imaging artifacts 4 .
Many technical improvements have been achieved in LSM within the last decade -on the detection side, but especially on the illumination side [4][5][6][7] . It turned out that scanning a beam 8 through the focal plane results in a more homogeneous light-sheet 4 than by forming a static light-sheet generated by a cylindrical lens, known as selective plane illumination microscopy (SPIM). Besides the illumination with a scanned Gaussian beam, scanned Bessel beams offer intriguing advantages such as high penetration depths into media due to their self-reconstruction capability 5 . In addition, Bessel beams can generate more imaging contrast and resolution, when combined with confocal line detection 9 , two-photon excitation 10 structured illumination 6,11 or using coherent superposition to form a lattice light-sheet 12 .
Especially the use of beam shaping elements such as spatial light modulators (SLMs) or digital micromirror devices (DMDs) allows to generate nearly arbitrary illumination beams, which, in the optimal case, can be adapted 1 Scientific RepoRts | 6:30378 | DOI: 10.1038/srep30378 specifically to the refractive index inhomogeneities of the object, which would otherwise lead to unwanted beam deflections and distortions 13 . This holographic shaping of the illumination beam requires coherent light, which is characterized by its unique phase dependency between all photons incident onto the object. However, fluorescence light emitted from the fluorophores inside the extended object lack this mutual phase dependency, since the emission of fluorophores is incoherent. Wavefront correction approaches on the detection side can help to improve the image quality 14,15 , but the significant computational efforts and several illumination iterations per image plane set strong limitations concerning the acquisition speed and hamper imaging of dynamic processes. Popular methods to improve the image quality by post-processing such as image deconvolution [16][17][18] often remain unsatisfying because the point-spread function (PSF) at the camera is different for each position inside the 3D object, which is because fluorescence photons are scattered differently at each position. Nevertheless, post-processive deconvolution does not affect imaging speed.
The influence of scattering becomes apparent in an illustrative way when imaging large, scattering structures such as small plants and plant root tips, that consist of regularly shaped cells of variable sizes 19 , but reveal high refractive indices (high polarizabilities) leading to strong deviations of the illumination and detection photons. These effects degrade the image quality of those parts of the plant which are further away from the illumination lens and the detection lens. For example the small and dense pericycle cells within Arabidopsis roots appear blurry 19 . Lateral root initiation originates from the pericycle. Hormone gradients, which are established by polarity in the cell's localized transport proteins of the PIN family, play a crucial role in this process 20 . Detailed real-time analysis of PIN proteins could not be performed so far because of limited sub cellular resolution and contrast. Although illumination from 4-5 sides helps to improve the image quality in the back parts of the object, this comes at the cost of a 4-5 -fold increase in illumination time and bleaching in addition to intensive postprocessing 2 .
Image formation in scattering media has been described by wave-optical models 21,22 and investigated via Monte Carlo simulations 23 . An effective PSF has been introduced 24,25 that allows image estimation 25 and deconvolution 26,27 under different scattering conditions. Image formation in scattering samples can be better understood by quantifying the amount of scattering events a photon has survived and considering this in specific PSFs 28 . However, a practical procedure to separate diffusive fluorescent photons or an analytical model describing their influence on the imaging process is still missing.
In this paper we introduce such a procedure and model. Based on the theoretical description of photon diffusion, we present a method that allows to separate efficiently the parts of the image resulting from the minimally scattered (ballistic) photons and the highly scattered (diffusive) photons. This information can be directly extracted from the image spectra such that no a priori knowledge about the object is required. We demonstrate the principle of separating weakly and strongly scattered fluorescence photons using identical Arabidopsis Thaliana roots, which we image in four different illumination and detection modes. We have developed a concept that quantifies the change of ballistic photons into diffusive photons, which varies with the detection depth and which is implemented in a depth dependent deconvolution of the 3-D images.
The Conceptual Approach. In this study we investigate the influence of three important components of the imaging chain -the illumination beam, the detection method and the image postprocessing -on the 3D image quality. Further, this chain is studied in four different imaging modes. All variations of these components are performed at exactly the same biological object, an Arabidopsis thaliana root, to enable a meaningful comparison of the imaging modes. The roots are about 120 μ m in diameter and scatter strongly both the illumination light and the fluorescent detection light, leading to reduced contrast and resolution, but also to local image artifacts in the back parts of the object.
As the first component of the imaging chain, we use a spatial light modulator (SLM) that can switch plane wise the shape of the illumination beam, which is either a Bessel beam or a Gaussian beam optimized to the extent of the object. The beam is scanned laterally in the plane of focus (see Scheme of Fig. 1). As the second component, we use a camera with a rolling shutter, which enables confocal line detection 9,29 by a narrow width of the slit moving parallel to the illumination beam (see Fig. 1) or conventional widefield detection for a large slit width.
The third component in the imaging chain is a deconvolution based on our model describing the photon diffusion in image formation. The required scattering parameters can be extracted by analyzing the width of the image spectra, which decrease with the detection depth and thereby represent the increase of image blur. We will show that enhancing high frequency components by deconvolution enhances those image components, which mainly originate from ballistic photons. This way, the object and depth dependent deconvolution can be interpreted as post-processive photon separation.

Results
3D imaging of root tips. Arabidopsis thaliana root tip is the preferred model for root development in plants due to its small size and simplicity. With conventional light-sheet microscopy only outer, but not inner cell layers were accessible for microscopic analysis with reasonable resolution 19 .
Here, we imaged Arabidopsis thaliana roots expressing the plasma membrane located protein LTi6b fused to GFP homogenously in all cells 30 . Seeds of 35s::LTi6b-eGFP) were germinated for 7 days on solid medium in the light. Excised roots were embedded in 1.5% low-melting agarose.
It is well known that Gaussian and Bessel beams exhibit different scattering properties in inhomogeneous media and thereby generate different light-sheets and images of different quality in LSM. The usage of a spatial light modulator (SLM) in the illumination path allows to shape nearly arbitrary illumination beams, which will then be scanned laterally across the object as indicated in Fig. 1. This offers the opportunity to directly compare the images acquired by using e.g. Gaussian or Bessel beam illumination 31 .
Scientific RepoRts | 6:30378 | DOI: 10.1038/srep30378 The concentric ring system around the thin main lobe of a Bessel beam enables the superior propagation stability and the self-reconstruction capability of the beam. However, it is also the ring system that generates the low image contrast because it excites fluorescence also out of the focal plane. This problem can be well solved by using a confocal line detection scheme, which mainly detects the fluorescence generated by the thin main lobe of the Bessel beam and generates unreached contrast and resolution 9 . Synchronizing the light-sheet mode of the camera (ORCA-Flash4.0, Hamamatsu) with the scanned illumination beam allows confocal line detection at high framerates 29 and excellent sectioning. In this way 3D image stacks including image aberrations can be obtained with high sampling also in a plane vertical to the light-sheet revealing so far unknown imaging peculiarities on a sub-micron scale.
In total four 3D image stacks each consisting of 260 planes separated by Δ y = 0.5 μ m were recorded from a single root tip. To allow a fair and precise comparison without bleaching or drift effects, each plane was imaged successively with the above mentioned four different imaging modes. Both the Gaussian (axial full width at half maximum (FWHM) Δ z = 300 μ m) and the Bessel beam (axial FWHM Δ z = 300 μ m; numerical aperture NA = 0.15) illumination have been used to image in conventional widefield mode (detection slit width d slit = M T · 50 μ m) and in confocal line detection mode (detection slit width d slit = 8 · 6.5 μ m = M T · 1.3 μ m corresponding to 8 pixels on the camera, at magnification M T = 40). Figure 2a,b show image slices of an Arabidopsis root tip parallel to the light-sheet with Gaussian and Bessel illumination and conventional detection. Whereas the Gaussian beam reveals the better image contrast, the Bessel beam image hardly shows artifacts caused by scattering of the illumination light as visible with the Gaussian beam marked with a white arrow. In the top right corner the Fourier transforms of the images are displayed (in the spatial frequency coordinates k x and k y ) revealing a broader image spectrum (better transfer of high frequencies) for the conventional Gaussian beam.
The combination with confocal line detection results in a distinct contrast improvement for both the Gaussian beam and the Bessel beam as can be seen in Fig. 2c,d and also in the broadened image spectra. Although the image contrast is still better with the Gaussian beam, the scattering artifacts in the image become more pronounced especially in the back part of the image, since the Gaussian beams are stronger deflected and scattered than the Bessel beams 9 and thereby do not propagate straight along the detection slit of the camera. In consequence, significantly less ballistic photons are collected and the image reveals dark stripes. On the other side the image obtained with Bessel beam illumination combined with confocal line detection reveals a good contrast with much less artifacts since fluorescence from the ring system is efficiently blocked.
For a complete comparison of the beam propagation and imaging properties of all four imaging modes image slices perpendicular to the light-sheet have to be inspected.  Depth dependent image contrast. It was shown by the yz-slices of Fig. 2 that for all imaging modes a strong decrease in image contrast is going along with an increasing detection depth (y-axis). One method to provide a quantitative measure of contrast is the "useful contrast" introduced by Truong et al. 32 . A similar procedure is applied in our study and is further described in the methods section. This approach determines the contrast coefficient Q, which quantifies the ratio of the high spatial frequency image components (HSF, k ⊥ > k c ) and the low spatial frequency image components (LSF, k ⊥ ≤ k c ) (see eq. (11)). The corner-frequency k c = 2π /d cell used for the spatial filtering is defined by the largest object structure, which is the longer cell length of d cell ≈ 20 μ m (see inset of Fig. 3b).
Therefore, a high value of Q refers to high contrast in an xz-image slice. The contrast coefficient is plotted in Fig. 3a for all imaging modes over the detection depth y 0 (y 0 is defined in Fig. 2f). As expected, a strong decrease in contrast is visible with increasing detection depth. The drop in contrast Q caused by a few 10 micrometers propagation through the object is even stronger than the contrast difference caused by the imaging mode. This effect points out that the compensation for the detection depth induced loss in contrast is of great importance and might be differently effective for each imaging mode.
For a quantitative comparison of the different imaging modes, the contrast is normalized by the contrast coefficient of the standard imaging mode (Gaussian illumination and conventional detection). The result is plotted in Fig. 3b. It is obvious that confocal line detection improves the contrast. At low detection depth (y 0 = 10 μ m) the contrast is improved roughly by 70% for Bessel illumination (compare the bright and dark blue lines) and by 30% for Gaussian illumination (the red line above the green line). The contrast improvement is further enhanced for higher detection depths. For e.g. y 0 = 100 μ m the contrast coefficient is improved by more than 150% for both illumination modes. For higher detection depth the amount of diffusive photons involved in the image formation increases and image quality drops down. Since diffusive photons are multiply scattered on the way to the detector and thus are displaced in the image plane, they have a higher probability to be blocked by the confocal slit. This effect leads to the enhanced contrast improvement by confocal detection for higher detection depths. Since the slit is parallel to the z-axis, only photons displaced in x-direction are blocked thereby increasing optical resolution along x. This effect is shown in Fig. 3c, where the 1/e width of the spectra in all directions of the k x k z -plane are plotted relative to the spectrum obtained by the standard imaging mode (see the green circle as a reference). The elongation of the confocal image spectra in k x -direction demonstrates that the improved contrast enhancement at high detection depths is caused by the confocal slit, which works in in x-, but not in z-direction. A detailed discussion on the effect of confocal detection with special attention on the difference between Bessel and Gaussian illumination can be found in Supplementary Text 2.
Image formation with ballistic and diffusive photons. The image is formed by ballistic (hardly scattered) and diffusive (multiply scattered) photons. Due to the fact that diffusive photons carry only little high-frequency information about the object, the image quality decreases if the amount of diffusive photons involved in the imaging process increases. Consequently, the image quality is increased if the percentage of ballistic photons is enhanced. This process is called gating 33,34 . Confocal detection is a first step in this direction since the detection slit blocks photons, which are displaced by multiple scattering. A more advanced method is separating ballistic and diffusive photons by their time of flight 35,36 but technical complexity has prevented a combination with LSM. However, an alternative, which needs no further hardware, is separating the effect of the ballistic and diffusive photons on the image in a post processing step. The defined suppression of low-frequency information from the diffusive photons and the defined enhancement of high-frequency information from the ballistic photons is nothing else than a deconvolution with a PSF describing both the optical response and the detection depth dependent object response. The contrast coefficient is normalized to that of the standard Gaussian illumination with conventional detection. An enhanced contrast improvement by confocal detection is observed for high detection depths. The first inset illustrates the separation of useful high spatial frequency (HSF) information and low spatial frequency (LSF) background in Fourier space. The ratio of the mean value of the HSF and LSF components in Fourier domain gives the contrast coefficient. The second inset shows the relevant coordinate definitions in Fourier space. (c) Width of image spectra from 2D slices at different detection depths normalized to the Gaussian illumination with conventional detection. The 1/e width in all directions is obtained by exponential fitting in all directions. The elliptical shape of the confocal spectra for high detection depths indicates the contrast improvement perpendicular to the detection slit (x-direction).
Scientific RepoRts | 6:30378 | DOI: 10.1038/srep30378 A diffusive photon differs from a ballistic photon by its propagation angle θ . Any angular change during the propagation of a fluorescence photon on its way to the detector will result in a wrong position on the camera and to a wrong image contribution. The image is described by where h ill (r) and h det (r) are the purely optical response functions for light-sheet illumination and detection. (Further explanation of the mathematical description of the image process can be found in the methods section and in Supplementary Text 6) The contribution from both the ballistic and diffusive photons is defined by the scattering properties of the object f(r) and can be accounted for quantitatively by an object response function h obj (r, y 0 ), which consists of two parts. One part for diffusive photons, which is nearly of Gaussian shape and broadens with the detection depth y 0 and one part describing the ballistic photons by a delta-point function. This is explained in detail in the methods section. The degree of the linearly increasing influence of the object, i.e. the broadening of h obj (r, y 0 ) is expressed by the scattering parameter γ, which is a material constant specific for the object. The corresponding object transfer function H obj (k r , y 0 ), i.e. the Fourier transform of the response function, describes the signal loss for each spatial frequency k r . It consists of a Gaussian like part for the diffusive photons that narrows with y 0 and a constant offset for the ballistic photons that drops off exponentially with y 0 : (1) The shape and behavior of the object transfer function H obj (k ⊥ , y 0 ) in lateral direction is shown in Fig. 4a for three different detection depths y 0 . The experimental data is obtained through angular averaging of the normalized image spectra. Figure 4b,c illustrate photon propagation through the object, where a small fraction c 0 of the fluorescent photons leaves the object unscattered, and fractions c j are scattered j times. The parameter μ sca describes the probability that a photon is scattered. Each scattering event results in a random change of the propagation direction θ . The probability for a directional change follows a Gaussian distribution p(θ ), characterized by the material parameter γ (see Fig. 4d). In detail explanation of eq. (1) can be found in the method section.
Depth dependent frequency transfer. The image spectra shown as insets in Fig. 3a for y 0 = 10 μ m, 55 μ m and 100 μ m reveal that the spectral widths (xz-slices) decrease for higher detection depths. This means that a reduced amount of high frequency information is transferred through the scattering sample. To extract this loss in high frequencies caused by the object, the spectra must be normalized by a reference spectrum at the bottom of the object (see eq. (12)), where the influence of the object is assumed to be negligibly small. Division by this reference spectrum eliminates the effect of the microscope and the frequency spectrum of the object. Hereby we assume that the object frequencies do not change significantly with the detection depth. Deviations from this assumption are elaborated in the discussion. The obtained frequency transfer is plotted as scattered markers for three different detection depths (y 0 = 10 μ m, 30 μ m and 50 μ m in Fig. 4b). The data has been extracted from the 3D image data captured with Gaussian illumination and conventional detection. The frequency transfer is independent on the direction and can be plotted over the radial spatial frequency = + ⊥ k k k x z 2 2 . As expected from theory (see method section), the frequency transfer consists of a Gaussian like part centered around k ⊥ = 0 describing the diffusive photons and a constant frequency transfer. Since ballistic photons are not affected by the object, their spatial frequencies do not change, which results in a constant frequency transfer function. The width of the Gaussian part reduces with increasing detection depth y 0 and the constant level drops exponentially with y 0 as expected from Lambert-Beer's law.
Extracting the object transfer function allows to compensate for the reduced high frequency transfer through the object by applying a deconvolution according to our model for an effective system PSF. Therefore, the scattering parameters μ sca and γ which describe the scattering behavior averaged over the 3D image have to be known. If no such parameters are available for example for heterogeneous materials like Arabidopsis roots the parameters can be extracted from the images. This procedure is described in the methods part and in Supplementary Text 3. By fitting H obj (k ⊥ , y 0 ) to the frequency transfer data, we found μ sca = 50 mm −1 and γ = 22.4. The model-fitted curve is plotted in Fig. 4a with solid lines. Supplementary Movie 1 shows the fitting result over the whole range of the detection depth. As mentioned above, the frequency transfer has been calculated for the conventional detection mode since it is independent of the direction in the k x k z -plane. This is different for confocal detection. Nevertheless, it is possible to extract the scattering parameters by only considering the frequency transfer in k z -direction (See Supplementary Text 3 for more details).

Depth dependent deconvolution.
With knowledge of the scattering parameters μ sca and γ, a depth dependent deconvolution is possible. Here, we used a Wiener filter (see methods section) instead of more advanced deconvolution algorithms, to make the effect of photon separation better visible. The deconvolution intensifies the influence of ballistic photons on the imaging process at high detection depths, such that it can be classified as a gating technique. Figure 5a-d show how contrast of xz-slices is increased by a depth-dependent deconvolution. It displays how the enhancement of high frequencies is adapted to the blur which increases with the detection depth. Figure 5e-h show yz-slices through 3D stacks, which have been processed by a depth-dependent deconvolution. 12 different y 0 -position and 12 different PSFs (see methods section) have been used. The images reveal a strong increase of contrast compared to the unprocessed images displayed in Fig. 2e-h. Some overshooting is visible in the left side of Fig. 5h. This is due to the fact that the round shape of the root is not considered and the effective detection depth in the area, where overshooting appears, is shorter than expected. Additional scattering of the illumination light, which is not considered, leads to a widening of the illumination beam from left to right, causing a gradient of image quality. Compensating for the scattering of illumination light needs a separate discussion, because of the coherent nature of the illumination beam.

Discussion
3D imaging of root tips with and without confocal detection. The advantages of confocal line detection using a camera with a moving slit (rolling shutter) become clearly visible in Fig. 2. The image contrast improves, the lateral image spectra become broader (see insets) and especially the resolution in detection direction is improved (Fig. 2g,h). However, this comes at the cost of stronger stripe artifacts for Gaussian beam illumination, which do not exhibit the directional stability of Bessel beams. Gaussian beams are deflected during propagation and after some distance do not propagate parallel to the detection slit. In addition, it becomes apparent by the yz-slices of Fig. 2 that Gaussian beams cannot illuminate entirely the cell membranes oriented parallel to their propagation direction. However, this effect is not further investigated here. Bessel beams on the other side provide good contrast over the extent of the 100 μ m thick root without producing any significant artifacts, which is the consequence of their propagation stability and their self-reconstruction capability.
A recent study using structured light-sheet illumination reveals a strong image darkening in the central part of the root tips, since the illumination modulation frequency, does not match the PSF broadening. In our study, however, the small central cells can be reasonably resolved 19 . In structured illumination, the frequency of the illumination pattern has to be adapted to the PSF width and thereby also to the detection depth in order to exploit its full power, which would be possible with our model.
The gating effect of a confocal slit, i.e. the blocking of diffusive photons, has been analyzed also theoretically. It was shown by Supplementary Figs 3 and 4 that contrast improvement by confocal detection using Gaussian illumination is mainly caused by blocking diffusive photons, but not by blocking fluorescence excited out of the focal plane. In other words, confocal detection in weakly scattering media with Gaussian illumination has almost Scientific RepoRts | 6:30378 | DOI: 10.1038/srep30378 no effect. In contrast Bessel illumination with confocal line detection benefits from the true confocal effect and gating further improves contrast.
Separating fluorescence photons in the image. Ideal, diffraction limited imaging is only possible by ballistic photons, which transfer information from the focal plane inside the object to the detector. Since the number of ballistic photons at the detector decreases exponentially with the detection depth y 0 , the useful image signal decays in the same way. The image background, defined by the diffusive photons, will increase exponentially with the detection depth. The separation between ballistic and diffusive photons, i.e. a strong enhancement in contrast, is possible, as long as the amount of ballistic photons is distinguishable from the noise level.
To extract the frequency transfer through the object each image spectrum at y 0 is divided by the image spectrum at y 0 ≈ 0, which is mainly defined by the imaging optics. In principle, the normalized spectrum can be obtained experimentally without the fit functions derived in our model. However, the fit functions provide an additional control and deliver the material parameters μ sca and γ. Figure 4 summarizes the change of the object transfer function for different detection depths y 0 . Both the experimentally and theoretically obtained H obj (k r , y 0 ) reveal a quite good coincidence in shape and thereby demonstrate the ongoing transfer of ballistic photons to diffusive photons with increasing y 0 . Deviations between the measurement data and fit curves can be first a consequence of the inhomogeneous scattering properties of the root tip and second, the characteristic structure size and thus the object spectrum changes with detection depth, such that the effect of the object spectrum cannot be completely eliminated by scaling with a reference spectrum. For example, cells in the middle of the root (Endodermis, Pericycle and Stele) have a size of 5 to 7 μ m and are smaller than the cells in the two outer cell layers of the root (Epidermis and Protoxylem), which measure 10 to 20 μ m. Since the reference spectrum was obtained from an outer cell layer, frequencies from 0.05 to 0.1 μ m −1 are more pronounced than in a layer at 50 μ m detection depth (green curve in Fig. 2b)). Thus, in the normalized spectrum at y 0 = 50 μ m frequencies from 0.05 to 0.1 μ m −1 are reduced and frequencies from 0.14 to 0.2 μ m −1 occur amplified. The fitting algorithm turned out to be robust against these deviations. In practice, there are several objects offering depth independent object spectra (e.g. Zebrafish with labeled nuclei) and thus are even more suitable for the presented method. Nevertheless, by imaging membrane labeled root tips, we demonstrate that the method is not limited to such samples.

Depth dependent deconvolution.
Deconvolution with a single response function generates images with either too little or too much contrast (effect of overshooting) depending on the choice of the PSF. However, by adapting the object PSF to the corresponding detection depth, resolution and contrast in the image can be optimized over the complete 3D volume. For illustration of this effect, see Supplementary Text 4. The depth dependent deconvolution algorithm was tested successfully also at a dense bead cluster. Results are shown and discussed in Supplementary Text 5.
Our results prove that scattering of fluorescent light can be well approximated by a PSF even for inhomogeneous samples like Arabidopsis roots. In comparison to other deconvolution techniques, which utilize point like structures in the object to extract the PSFs 37 our technique can be applied to any structure without inserting beads. Recently developed blind deconvolution algorithms are based on a depth-variant PSF model 38 . Here, spherical aberrations from refractive index steps between immersion medium and sample were considered multiple scattering a main source for image degradation is neglected.
A perfect depth-dependent deconvolution would require a complete 3D deconvolution for all N xz-planes. A solution with less computational effort is described in the methods part, where twelve 3D deconvolutions for the whole stack have been performed. However, it is not necessary to compute twelve 3D deconvolutions for the whole stack, since only a limited number of planes around the specific plane at y 0 is of interest. A further reduction of computational efforts is possible, if each of the twelve 3D deconvolution is performed on a section, which is six times thinner than the whole stack. In the end, the computational effort will increase only by a factor of two relative to a standard 3D deconvolution with a single PSF.
Furthermore, all image processing operations can be outsourced directly to camera systems exploiting fast and flexible power of GPUs or FPGAs.
By the strong increase in contrast it becomes apparent that Bessel illumination with confocal detection benefits most from the deconvolution and leads to a nearly artefact-free image. The already mentioned artefacts caused by Gaussian illumination are even more obvious after deconvolution. Only without deconvolution, one may argue whether the high contrast achieved by Gaussian illumination or the elimination of artefacts by Bessel illumination leads to the best image quality. In general, this question has to be answered individually for each application. After deconvolution the best image quality is clearly achieved by Bessel illumination. Fig. 4a that the model-based

Extracting scattering parameters out of the image. It was shown in
sca introduced in eq. (1) fits reasonably well to the frequency transfer extracted from the experimental data. The decrease of the offset level allowed to extract μ sca , whereas the decrease in Gaussian widths allowed to extract γ. We tested the universality of our approach by using the two different illumination and detection modes, to extract the same material parameters of the same root tip. Because of the gating effect by confocal detection the influence of the various scattering orders on image formation is changed. Thus the weighting by the photon fractions in the superposition of eq. (1) has to be modified in order to enable the fitting of the frequency transfer extracted from confocal images.
Alternative methods like collimated transmission measurements and goniometry 39 require thin samples. Our model is also applicable to thick samples, which can be described neither by the quasi ballistic regime nor by the diffusive regime 40 .
The results are summarized in the following table (Table 1) and indicate that the extracted material parameters μ sca and γ agree with each other reasonably. Further details and illustrating movies can be found in the Supplementary Text 3. In Supplementary Text 5, the scattering parameters of a bead cluster were extracted. It was shown that the scattering parameters are close to the values predicted by Mie theory.
It will be interesting in biology to analyze both scattering parameters over a longer time, e.g. during root tip development and during responses to external stimuli, since they provide biophysical, structural and possibly mechanical information about the state of the cellular compound. This information is usually not visible with standard fluorescence imaging. It will be an interesting future task to extract this information spatially resolved over areas of different cellular morphology within the root tip.

Conclusion
A severe and unsolved problem is the multiple scattering of fluorescent light on the way to the detector, thereby blurring the images. This effect is the stronger, the larger the object and the longer the way of the photons through the scattering object. We have presented an effective and elegant solution to this problem, leading to a significant improvement in image quality. We separated the image contributions of the diffusive and the ballistic photons by postprocessing. Our approach, which does not require a priori information, neither slows down the acquisition speed of light sheet microscopy, nor are any iterative image acquisitions required, as this is the case with adaptive optical approaches. Our algorithm is applicable to most new and existing 3D data sets suffering from image blur due to scattering of fluorescence light. It should be beneficial to every light-sheet microscopist and could in principle be applied to other microscopy techniques.

Methods
Theory of image formation. Image generation in light-sheet microscopy. The formation of a 3D image p(r) can be described by a convolution (symbol * ) of a 3D object distribution f(r), e.g. the fluorescence emission distribution, with the system PSF h sys (r): sys The impulse response h sys (r) of the imaging system is assumed to be shift invariant and contains the information on resolution and contrast. In LSM an effective systems-PSF h sys (r) can be defined, which is given by , with T being the exposure time of the camera and d slit the slit width, if confocal detection is applied. For better readability the transversal magnification of the detection optics is set to |M| = 1. Thus the illumination PSF effective for one pixel line during T reads Image generation in combination with an object response function. If LSM is used to capture a 3D image stack of strongly scattering objects a significant decrease in image quality with increasing detection depth y 0 can be seen. This loss of image quality with y 0 corresponds to a low pass filtering, which can be modeled by a convolution with an additional object response function h obj (r, y 0 ) describing the influence of the object on the detection process. Consequently, hdet(r) is replaced by h det (r) * h obj (r, y 0 ) and the effective system PSF can be written as The fact that h sys is now a function of y 0 indicates that the image process is no longer a 3D shift invariant problem and h sys (r, y 0 ) should be used for one specific detection depth y 0 . Thus, the final 2D-image obtained from position y 0 is given by Eq. (5) shows that the influence of h obj (r, y 0 ) on the image process is highly anisotropic and depends on h sheet (r). A good light-sheet has only a small extent in detection direction y. Due to the multiplication in eq. (5), the spread caused by h obj (r, y 0 ) is limited to the light-sheet thickness. The same is true in scanning direction x, if confocal detection is applied (d slit → 0, h ill (r) → h SB (r)). This results in a gating effect.
Depth specific scattering of fluorescence photons. To understand the depth dependence of h obj (r, y 0 ), one has to consider the exponential decrease of the number of ballistic photons and the increase in diffusive photons that travel over the distance y 0 through the object towards the detector. The larger the detection depth y 0 , the larger is the average of the scattering order j, which indicates the number of scattering events a fluorescence photon involved in the imaging process has survived. In this way, c j (y 0 ) denotes the percentage of j times scattered photons at depth y 0 , such that c j=0 (y 0 ) is the fraction of ballistic photons contributing to the image. Since every scattering order carries a different amount of information, the image process for each order j is described by a specific object response h obj,j (r, y 0 ). Then the overall object response h obj (r, y 0 ) is given by the weighted superposition j j j obj 0 sca 0 0 sca obj, 0 As indicated by eq. (7), the fraction c j depends on the scattering coefficient μ sca , which describes the probability for a photon to be scattered while propagating a certain distance through the object. By neglecting absorption, one finds (see Supplementary Text 7) extends the Lambert-Beer law to higher scattering orders. Figure 6 shows an example of the distribution of fluorescent photons over the scattering orders given by eq. (8).
The second term in eq. (7) describes the change in shape of the specific object response h obj,j (r, y 0 ), which depends on the scattering order and can be approximated by obj,j 0 where A is an unimportant normalization factor. A detailed derivation for eq. (9) can be found in Supplementary Text 8. The broader the specific object response function, the more high-frequency information is lost. This loss in information increases with every scattering event and with the detection depth y 0 to the focal plane. It is important to note that for ballistic photons h obj,j=0 (r, y 0 ) becomes a δ -function. This is in agreement with the fact that these photons are not affected by the object. h obj,j (r, y 0 ) also depends on the material specific unit γ, which describes the angular width of the scattering phase function  Fig. 4b,d)). The relation to the anisotropy factor g HG , known from e.g. the Henyey-Greenstein function 42 , is given by = − γ γ ( ) g coth HG 2 2 2 2 with coth being the hyperbolic cotangent. We find γ → 0 for Rayleigh scatterers (PDF θ (θ) = 1) and γ → ∞ for forward scatterers (PDF θ (θ) = δ(θ)).
In Fourier domain one finds the depth dependent fluorescence photon object transfer function H obj (k r , y 0 )  contribution of the ballistic photons in real space is described by an infinitely small Gaussian function identical to a Dirac delta function. Accordingly in Fourier domain ballistic photons are considered by an infinitely wide Gaussian function or a constant offset. By use of the power series representation of the exponential function the object transfer function can be rewritten as follows: Here + = ∆ + ⌈ ⌉ n y y 1 / 1 obj = 12 is the number of sections of thickness Δ y and y obj is the object dimension in y-direction. The Wiener filter is given by