1.

Introduction

One of the most prominent limitations for light microscopy of biological samples is the imaging depth. Scattering in tissue is the limiting factor defining the depth in which biological structures may still be resolved. Larger samples can only be imaged superficially. Deconvolution algorithms for multiview selective plane illumination microscopy (SPIM) have already been shown to increase the ability to image thick samples, provided they are optically accessible from multiple angles.¹

SPIM^2–4 utilizes an extreme darkfield setup by illuminating only the focal plane in the sample by a light sheet. The illumination optics is usually positioned perpendicular to the detection objective. By scanning the sample through the light sheet, a three-dimensional (3-D) image stack of the sample can be acquired in a short period of time. The technique offers a high acquisition rate of entire volumes and low photo damage and bleaching, thus making it an ideal tool for in-vivo imaging of fast processes or developing organisms. In multiview SPIM,^1,5,6 the sample is rotated and multiple image stacks are recorded from different angles. Alternatively, multiple excitation and detection optics allow imaging of the sample from 2 to 4 directions without rotating the sample. The different views are later used to reconstruct a single volume. The most prominent technique to combine the images are blind deconvolution algorithms.

In this work, we introduce a regularizer and a prior to improve the reconstruction and make it more robust. For demonstration purposes, we concentrate on using a single SPIM image stack rather then using multiview data (although the regularizer is applicable to multiview data). This approach offers the advantage of using an SPIM setup without multiple detection optics and a fast aquisition without rotation of the sample. Since the SPIM image stack has to meet only some weak requirements for the reconstruction, it is also possible to improve already existing single view SPIM data.

2.

Reconstruction Algorithm

The scattering influence on intensity within a sample can be treated as a convolution between the real object and an increasingly broad point spread function (PSF). The PSF describes the effect an imaging system has on a point source. For SPIM, the PSF is shaped by three main contributors: (1) the illumination function of the sample. The shape of the light sheet as well as broadening of the light sheet within the sample due to scattering determines the illuminated areas of the sample. Only from these regions we can expect fluorescence. To acquire a full 3-D image of the sample, the light sheet and the sample are moved relative to each other. The direction of movement will be denoted as $z$-direction, while the imaging plane will be labeled the $x-y$ plane. The function describing the illumination will be called $P_{\mathrm{illum}}$. (2) Fluorescence molecules, which are illuminated in the current $z$-position, can be considered point sources. The light from these sources gets scattered on the way to the objective and will contribute to the PSF. This contribution will be denoted as $P_{\mathrm{scatt}}$. (3) Eventually, the detection optics contributes to the PSF. This contribution $P_{\det }$ is sample independent and—for most cases considered here—smaller than the PSF broadening due to scattering. In the following, all three separate contributions are condensed into one single PSF $P_z(x,y)$

2.1.

Spatially Invariant Algorithm

For the sake of clarity, the algorithmic details are first discussed using a spatially invariant model. In each $z$-plane, the detected SPIM image is treated as a convolution between the real object and the PSF. In this case, the real object is the distribution of fluorophores in the sample. An additional noise-term accounts for additive noise influences on the final image introduced by the imaging system

Eq. (1)

$$I_z=O_z\times P_z+N,$$

where $N$ is the noise, $I$ is the image, $O$ is the real object, and $P$ is the PSF.

If $N$ is Gaussian noise, the original image can be estimated by minimizing the least-square error $E$.

Eq. (2)

$$E={\Vert I_z-P_z\times O_z\Vert }_2^2.$$

It should be noted that image noise is typically not Gaussian but Poisson distributed. Nevertheless, the image space reconstruction algorithm (ISRA) which will be introduced below turned out to also be stable in the case of Poisson noise. This reconstruction framework is easy to implement and will be used to demonstrate the function of our regularizer and the prior we use. Equation (2) is typically ill conditioned, even if the PSF $P$ is known. In most cases, multiple views are used to provide the necessary amount of information in order for the algorithm to converge. In our case, the acquired image stack $I$ is the only known variable. Trying to minimize the least-square error will not necessarily converge to one solution. To obtain a converging solution, a prior is used. The main assumption for this prior is that if the step size between two consecutive $z$-planes is small enough so it can be assumed that the same object is illuminated in both positions, while the PSF varies. This idea of using neighboring $z$-planes is depicted in the illustration of the algorithm in Fig. 1(b). Under ideal conditions and under the assumption of multiple observations, $P_i(i=1,2,3)$, which are neighboring $z$-planes, in our case, the coprimeness equation: $P_j\times I_i-P_i\times I_j=0$, $\forall i\ne j$, can be used to determine the PSF. In the presence of noise, this equation becomes an optimization problem

$$E_{\text{coprime}}=\sum _{i\ne j}{\Vert P_j*I_i-P_i*I_j\Vert }_2^2.$$

The effectiveness of algorithms based on this prior has already been shown by multiple groups.^7–9 However, they use either different views or time series of the same view with slightly varying PSF. Due to the physical properties of the imaging system, it can additionally be assumed that PSF $P\ge 0$, as well as $\sum _{x,y}P(x,y)=1$. The sum constraint is only valid if absorption in the sample can be neglected, which is reasonable for adjacent $z$-planes. The positive constraint is also true for the object $O$. Furthermore, the support of the PSF should be smaller than the image $I$. These constraints together with the assumptions of a constant object over neighboring $z$-planes [whereas the PSF $P_i(i=1,2,3)$ changes] lead to the following minimization problem:

$$\underset{P_i\ge 0,\sum P_i=1,O\ge 0}{\min }\sum _i{\Vert I_i-P_i\times O\Vert }_2^2+\sum _{i\ne j}{\Vert P_j\times I_i-P_i\times I_j\Vert }_2^2.$$

Fig. 1

(a) Schematic of the SPIM setup and (b) illustration of the algorithm with the quadratically increasing regularization (graph to the left), the coprimness prior (to neighboring slices are used to reconstruct the red plane) and the iteration scheme to the right.

The last term, the coprimeness regularization mentioned above, helps to stabilize the deconvolution problem and can even lead to a unique solution.⁷

To ensure the convergence of both the PSF $P$ and the object $O$, an alternating optimization scheme is used, where the optimization iterates between $P$ and $O$. In this work, the optimization is done with a modified nonnegative least-square estimator, the ISRA.^10–12 In the first step, the PSF is updated by applying the ISRA

$$P_k(t+1)=P_k(t)\frac{O\times O^{*}\times I_k+{\sum }_{j\ne k}{I}_j\times I_j^{*}\times P_j(t)\times I_k}{(O\times O^{*}+{\sum }_{(j\ne k)}{\mathrm{I}}_{\mathrm{j}}\times {\mathrm{I}}_{\mathrm{j}}^{*})\times P_k(t)}.$$

In the second iteration step, the PSF is fixed and the object is updated in the same way. The iteration scheme is illustrated in Fig. 1(b) on the right. To ensure convergence of both the PSF and the object, a proximal minimization method is applied to the algorithm.^13,14 Alternatively to ISRA, a Richardson–Lucy deconvolution could be used for deconvolution which should result in a slightly longer computation time.

With increasing imaging depth, the images will suffer from increased scattering. The idea in this work is to compensate for this increased scattering by adjusting the algorithm along with the imaging depth. A simple approximation for light scattering in tissue is found in Refs. 15 and 16. With increasing depth, the size of the PSF grows quadratically [see Fig. 1(b), where the increase in regularization with depth is schematically plotted]. This simple model was implemented by the use of the Tikhonov regularization.¹⁷ A higher Tikhonov regularization leads to an increased PSF size in the reconstruction. This is not an intrinsic effect of the regularization itself, but a combination effect with the constant intensity constraint $\sum P_i=1$. The ratio between the l1-norm (sum constraint) and the l2-norm (Tikhonov regularization) is a measure of signal sparsity with higher values describing sparser signals. Since the Tikhonov term reduces the l2-norm of the solution (this is penalized in the final cost function), the regularized solution automatically has a lower sparsity, i.e., it has a larger support (broader blur).

To reflect the increase in PSF size, the regularization parameter was increased quadratically, which was empirically found to be an effective way to model the increased scattering. However, no assumptions about the exact shape of the PSF are made. As a result of this, even elliptical-shaped PSF (due to aberrations or motion) can be corrected.

3.

Results

We demonstrate the effect of the algorithm on a real dataset of a common model organism for biological research, a transgenic zebrafish embryo. The 3-day-old embryo of the type Fli1A had its vasculature labeled with enhanced green fluorescent protein (EGFP). Imaging was done ex vivo with SPIM. The microscope setup, which is schematically shown in Fig. 1(a), was built by ZEISS, (Carl Zeiss AG, Jena, Germany) with a custom-made Apo Z $1.5\times /0.37$ objective (free working distance 30 mm). Two cylindrical lenses ($f=30\mathrm{mm}$) are used to achieve a double-sided illumination with a light sheet thickness of $16\mu \mathrm{m}$ (full width at half maximum). For detection, an sCMOS camera (pco.egde, PCO, Kelheim, Germany) is mounted on the microscope and the raw data are saved as a tiff image stack. The fish was imaged sideways with a step size of $5\mu \mathrm{m}$ between adjacent $z$-planes. Two complete 3-D image stacks of the zebrafish embryo were acquired, and the sample was rotated by 180 deg between the image stacks. While reconstruction was only performed on the first image stack (S1), the second image stack (S2) was used as a control image. Since the imaging was done from two opposed sites, structures which appear deep and scattered in S1 are more visible in S2. This circumstance is used to validate the reconstruction results from S1. To be able to use the second dataset (S2) for validation, a region from the tail of the fish with a width of about $175\mu \mathrm{m}$ was selected. Figures 2(a)–2(c) show image slices from three different depths: 105, 130, and $155\mu \mathrm{m}$, respectively. For this sample, the depicted depths are just in the range where scattering starts to drastically reduce the image quality. In the left column, the original slice from dataset S1 is presented, the column in the middle shows the same slice after processing with the introduced algorithm, and the column to the right shows the corresponding slice from dataset S2.

Fig. 2

The images show details of a zebrafish from the genetic line Fli1A 3 days post fertilization. The vasculature is labeled with EGFP. (a)–(c) show the same region from the tail of a zebrafish at different depths (105, 130, $155\mu \mathrm{m}$). The column on the left shows the original dataset (S1), the column in the middle shows the reconstruction of the dataset to the left, and the column on the right side shows a second dataset (S2), which was recorded from the opposite side of the fish to validate the reconstruction results. (d) MIP of the original dataset. The red-colored box indicated the region which can be seen in (e) and (f). In (e), the MIP of the reconstructed image stack can be seen, in (f), the MIP of the original image stack. The red arrows indicate areas where the improvement can clearly be seen. (g) and (h) show cuts through the slices in (a) and (b), as indicated by the dashed lines. The blue graph corresponds to the original dataset S1, red to our reconstruction, and green represents the control dataset S2. Scale bar indicates $100\mu \mathrm{m}$.

The section shown was reconstructed using 25 windows. The red arrows indicate regions where the reconstructed dataset shows a clear improvement in image quality compared to the original dataset S1. The structures which can be seen in $105\text{-}\mu \mathrm{m}$ depth agree well with the structures in the second dataset S2. The effect is still visible at a depth of $130\mu \mathrm{m}$. The good agreement with the control dataset S2 can also be seen in Figs. 2(g) and 2(h), where line profiles through the slices at 105- [Fig. 2(g)] and $130\text{-}\mu \mathrm{m}$ depth [Fig. 2(h)] are plotted. The original dataset S1 is plotted in blue, red represents our reconstruction, and green corresponds to the control dataset S2. The position of the chosen cut is indicated by a white dashed line in Figs. 2(a) and 2(b). Deeper inside the sample, the ability of the algorithm to correct the image diminishes. At a depth of $155\mu \mathrm{m}$, the image improvement is no longer able to reveal detailed structures.

The adaptive regularization prevents overregularization on image planes closer to the surface. Images from all planes benefit from the used prior, since the algorithm can correct for small motion artifacts which alter the PSF shape of neighboring planes. A region of interest covering part of the head and the eye of the zebrafish is marked red in Fig. 2(d). Figure 2(f) shows a projection of the original image stack. A projection of the reconstructed image stack is depicted in Fig. 2(e). The overall thickness of the head is around $320\mu \mathrm{m}$. The blur due to scattering and other effects is clearly reduced in the reconstructed image and virtually no artifacts can be seen. Two red arrows mark small details which are not/only partly visible in the original maximum intensity projection (MIP) but are clearly visible in the reconstructed MIP.

To further highlight the advantage over deconvolution algorithms without the suggested prior and regularization, two slices from the head region of the zebrafish are shown in Fig. 3. In the upper row, Figs. 3(a)–3(d) depict a plane from a depth of around $140\mu \mathrm{m}$. The lower row in Figs. 3(f)–3(i) show an image plane from the eye close to the surface at a depth of 5 to $10\mu \mathrm{m}$. The original image slices can be seen in Figs. 3(a) and 3(f), and the result from our algorithm in Figs. 3(b) and 3(g). The two planes in Figs. 3(c) and 3(h) were reconstructed without the coprimeness prior and with a fixed PSF. Instead of iteratively reconstructing both the PSF and the image, a fixed (simulated) Gaussian PSF was used. This kind of correction is widely used in microscopic imaging to correct the optical aberrations of the microscope. Figures 3(d) and 3(i) show a modified version of our algorithm where we kept the prior but used a constant regularization. The amount of regularization was optimized for a depth of around $140\mu \mathrm{m}$. To visualize the difference between the three reconstructions, Figs. 3(e) and 3(j) show a cut through the images. The location of this cut though is indicated by a white dashed line in the original image. In Fig. 3(e), a cut through the plane of $140\mu \mathrm{m}$ depth is shown. The original dataset is plotted in blue, the reconstruction with fixed PSF in green, and the results from our algorithm in red. The results from the constant regularization are identical to the results from our algorithm since the regularization was optimized for the depth of $140\mu \mathrm{m}$. As expected, our algorithm enhances structural contrast and is reducing the structure size to correct the scattering, while the reconstruction with the fixed PSF falls behind. In this depth, the PSF broadening due to scattering is much larger than the microscope PSF, therefore, the PSF size in the fixed reconstruction is not sufficient to yield a better result.

Fig. 3

Comparison between different deconvolution methods. The head region of the zebrafish shown in Figure 2 is chosen and two different slices from the volume are picked. The upper row (a)–(d) shows a plane from a depth of $\sim 140\mu \mathrm{m}$. The lower row (f)–(i) shows a plane from the surface (of the eye). In (e) and (j), a cut through the image planes is plotted. The location of the cut is indicated by the white dotted line in the original images. In the plot, blue color corresponds to the original dataset, red to the results from our algorithm, green to the reconstruction with fixed PSF, and cyan to the reconstruction with constant regularization. Artifacts produced by the constant regularization are marked with a red arrow. Scale bar indicates $100\mu \mathrm{m}$

The effect of scattering can be neglected at a depth of only $5\mu \mathrm{m}$. Here, we assume the original image Fig. 3(f) to be very close to the ground truth.

Applying excessive regularization can lead to artifacts in planes close to the surface. Figure 3(j) shows a cut through to one of these planes at a depth of $5\mu \mathrm{m}$. The reconstruction with constant regularization is plotted in cyan. One artifact of overregularization could be observed in our simulation and which is also visible in Fig. 3(j) is the formation of oscillations near plateaus or structures of high intensity. The plot in Fig. 3(j) shows an example of such an oscillation at a position of around $50\mu \mathrm{m}$ (marked the red arrow). Our algorithm and the fixed PSF reconstruction are much closer to the original dataset at this point. In the image, the oscillations might appear as a “shadow” around bright structures and spatially overlapping objects of different intensity are artificial separated.

The selected step size of $5\mu \mathrm{m}$ was sufficiently small for a stable reconstruction. In the shown datasets, most observed features were at least three times larger than the step size. For samples with smaller labeled structures, a smaller step size together with a thinner light sheet would probably be beneficial.

The results from the imaged zebrafish show that the algorithm works on $z$-stacks of SPIM images. However, this does not mean that the algorithm is limited to this imaging modality. The algorithm should be able to work with other optical microscopy methods, as well as other imaging methods which are limited by the scattering of light in tissue. There are three requirements which must be met: a full 3-D image stack is available, the step size between two adjacent planes is small compared to the observed structures, and the scattering in the sample can sufficiently be described by the introduced model. While a single dataset was used for this demonstration, the algorithm can be extended to handle multiple datasets to form one sample, either with varying illumination or even with varying views of the sample.

All animal experiments have been carried out according to national and institutional guidelines of the Government of Upper Bavaria and the highest standards of safety have been maintained during the course of the project.

4.

Conclusion

The acquired results show that the regularizer and prior presented in this work are able to partially correct for scattering within the sample and enhance the image quality and depth. The adaptive regularizer effectively avoids overregularization of the reconstruction. In the specific case of the zebrafish sample, the increase in imaging depth was about $30\mu \mathrm{m}$. For this correction, only information from a single image stack was used, together with a simple model for the dispersive broadening of light in tissue. Tikhonov regularization was used efficiently as a model to account for an increased scattering, which comes along with an increased imaging depth. The advantage of this approach is that neither multiple image stacks have to be acquired nor adaptations on the hardware level are required and the algorithm can be applied to already acquired data. The algorithm only requires a volumetric image stack with a small step size between two adjacent imaging planes, thus the algorithm should also be able to work on other imaging techniques which provide a 3-D image volume.