Blind deconvolution in autocorrelation inversion for multiview light‐sheet microscopy

Abstract Combining the information coming from multiview acquisitions is a problem of great interest in light‐sheet microscopy. Aligning the views and increasing the resolution of their fusion can be challenging, especially if the setup is not fully calibrated. Here, we tackle these issues by proposing a new reconstruction method based on autocorrelation inversion that avoids alignment procedures. On top of this, we add a blind deconvolution step to improve the resolution of the final reconstruction. Our method permits us to achieve inherently aligned, highly resolved reconstructions while, at the same time, estimating the unknown point‐spread function of the system. Research Highlights We tackle the problem of multiview light‐sheet deconvolution with a blind approach of autocorrelation inversion Our method recovers the object and PSF, requires no alignment and calibration, and enhances the reconstruction of the specimen.


| INTRODUCTION
Image formation in microscopy can be described mathematically by the convolution of the perfectly resolved object with the microscope's point spread function (PSF) (Mertz, 2019). The latter represents the impulse response of the system and, ideally, should be as narrow as possible to obtain the best image resolution. When the PSF diverges from the ideal delta function, image-restoration techniques may be used to improve the three-dimensional image quality. Among others, deconvolution (the process of inverting the convolution) aims to retrieve a sharper estimation of the object, reducing optical blur and noise by reassigning the energy distribution of the signal according to the PSF.
In Light Sheet Fluorescence Microscopy (LSFM) (Huisken & Stainier, 2009;Olarte et al., 2018), the illumination and detection optical paths are independent and orthogonal to each other: the PSF depends on both the illumination and detection optics and, typically, results in a function that is more elongated along the detection axis.
For this reason, multi-view acquisitions are frequently adopted as a method to obtain isotropic resolution (Calisesi et al., 2020;Swoger et al., 2007). Moreover, in the presence of slightly absorbing diffusive specimens (Rieckher et al., 2018), multiview acquisition facilitates the reconstruction of the entire imaged volume by accessing favorable views of the object. However, preprocessing steps are necessary to align and fuse all the acquisitions. In a typical multi-view LSFM pipeline, stacks of images are acquired from various angles, they are rotated, aligned, and combined together. Each of these steps requires dedicated care. The alignment procedure involves finding the rigid shift that guarantees the best overlaps of each view against each other. In general, the fusion can be performed by computing simple averages or with more advanced methods, for example, by considering Poisson statistics (Swoger et al., 2007) or optimized Bayesian strategies (Fernandez et al., 2010;Guo et al., 2020;Preibisch et al., 2014).
To complicate the reconstruction further, the knowledge of the PSF is not always accessible (Kundur & Hatzinakos, 1996). In a light sheet microscope, the acquisition of the PSF is performed by acquiring beads placed in the imaging chamber. However, this procedure may be impractical: this is the case in some clearing media (Ueda et al., 2020) where it is difficult to prepare beads phantoms, or in fluidic-based systems (Sala et al., 2020). In these cases, blind deconvolution is the only possible approach to increase imaging resolution (Kundur & Hatzinakos, 1996).
To solve the issues of alignment and fusion, we have proposed the Anchor Update (AU) algorithm . The protocol works out a deconvolved reconstruction from the autocorrelation estimation, exploiting the fact that working in the shift space guarantees implicit alignment. This implies that the solution of the autocorrelation problem always returns reconstructions that are aligned intrinsically at sub-pixel resolution [14].
In the present manuscript, we describe a strategy for tomographic reconstruction in multi-view acquisitions corrupted by unknown blurring, that can be utilized in those cases when the PSF is not easily measurable.
To this end, we design a blind deconvolution approach that acts in the autocorrelation space. This permits the formation of inherently aligned reconstructions, increasing the resolution simultaneously without any prior knowledge of the PSF of the system. We take inspiration from the Richardson-Lucy deconvolution (RLD) (Lucy, 1974;Richardson, 1972) and its blind implementation (Fish et al., 1995), extending the blind approach to the autocorrelation domain. First, we introduce the computational methods by validating the protocol thanks to the reconstruction of synthetic samples. Here, generated images for the object and PSF are used as ground truth to assess the convergence of our algorithm. Once validated, we test our reconstruction method on real data acquired by a custom LSFM setup, with the aim of retrieving both the PSF and the aligned reconstruction. Lastly, we discuss the potential of our method in the field of image processing for optical microscopy applications.

| Working principle of blind deconvolution
The goal of any blind deconvolution algorithm is to obtain a deblurred version of the original image without any knowledge of the response of the optical system used for image acquisition. It requires solving the deconvolution problem, o μ ¼ o Ã h, recovering o, the original object, from the blurred measurement o μ without knowing the blurring function h. There are infinite couples of images and PSFs that may reproduce the observed blurred output; thus, it may appear as a seemingly impossible problem to solve. It has been shown, however, that meaningful solutions are achievable by having weak, or even null, knowledge of the object under study and of the PSF (Kundur & Hatzinakos, 1996).
When the PSF is unknown, one of the most effective strategies is provided by a modified blind RLD method, which acts alternatively in the object and PSF domain. This protocol, introduced by Fish et al. (Rieckher et al., 2018), consists in dividing the deconvolution process into two steps, in which the object and the PSF commute their roles.
First, it executes a certain number of RLD iterations to improve the object reconstruction by using a fixed guess for the PSF. Then, the PSF estimation is improved by using the previously reconstructed object as the kernel of another RLD. This cycle is repeated several times until reaching optimal reconstructions.
Since we are interested in aligned multi-view measurements, we propose an extension of this blind method to the autocorrelation space. This approach is aimed at treating alignment and blind deconvolution together. As previously demonstrated by  in the case of known blurring, the AU algorithm inverts the autocorrelation of the measurement, obtaining an aligned and deconvolved reconstruction o rec . In this case, the problem that we try to solve is the inversion of the average autocorrelation, which we call χ μ , estimated from the measured views: In the previous equation, we have denoted the autocorrelation operator with the symbol ? and with m the total number of acquired views. In practice, the problem that we tackle is separating the object from the PSF given the average autocorrelation We address the problem as schematized in Figure 1, implementing the strategy via the iterated execution of two consecutive blocks. We keep the PSF fixed, and we reconstruct the object within the first block. We do the opposite for the second block, fixing the previously obtained object while focusing the reconstruction on the PSF. For the sake of clarity, here, we omit the index k that counts the iterations of the blocks. We start from the autocorrelation of the measurement χ μ by providing two initial guesses: one for the object, We cycle through these AU blocks several times until obtaining accurate reconstructions for the PSF and object. Since we assume that we have no information about the PSF, we will refer to this strategy as the Blind-AU reconstruction protocol. By denoting with the indexes t and t 0 the iterations within AU cycles, and with the subscript k the steps that run through the AU-blocks, formally, the Blind-AU can be written as: As previously mentioned, the object and the PSF are denoted with o and h, whereas we indicate their autocorrelation with the uppercase letters O and H. Here, the tilde specifies that the element is expressed with reversed coordinates: e K ¼ K Àx ð Þ. The code for blind deconvolution is available on GitHub (Ancora & Corbetta, 2021).
In the following section, we test our new Blind-AU method to reconstruct a synthetic sample, mimicking a misaligned acquisition with differently orientated point spread functions.

| Blind PSF and object reconstructions of synthetic data
In multi-view light-sheet microscopy, the sample is observed from different angles and always exhibits a point spread function elongated along the detection axis (Swoger et al., 2007). In Figure  F I G U R E 1 Schematics of the blind-AU algorithm. The algorithm starts receiving as input two initial guesses, one for the object and the other for the PSF. In the first step, the autocorrelated signal χ μ is deconvolved and deautocorrelated to improve the object estimation o k . After that, the same signal is used to restore the PSF, maintaining o k unchanged. Upon completion of the outer cycle indicated by arrows, the estimations o k and h k are retrieved, and the process iterated. We choose the number of iterations for each step and the number of times the outer cycle is performed depending on the experimental or simulated conditions F I G U R E 2 Synthetic sample for blind deconvolution. Blurred and noisy orthogonal views generated starting from an image simulating blood vessels. We test our image reconstruction method with this dataset, aiming at obtaining a tomographic slice sharper than the mere average of both images. The original section of the synthetic sample, which is our ground truth, is shown in Figure 3a. Aligning and fusing the orthogonal views leads to the blurred reconstruction in Figure 3b (top-left triangle). The alignment was accomplished by locating the peak of their cross-correlation and compensating for its opposite shift (Guizar-Sicairos et al., 2008). As a result of the fusion, the blurred reconstruction is affected by the combination of the two orthogonal PSFs that generate a star-like point spread function ( Figure 3c).
Our Blind-AU method allows for the reconstruction of both the object, sharper than each single view, and the PSF. Since our algorithm assumes that we do not have any information about the system, we set the initial guess of the PSF to be an isotropic Gaussian with σ ¼ 4 px ( Figure 3d). This initialization is larger than the ground truth PSF in both directions. On the other hand, we keep the aligned average as initial guess for the object. Ideally, our protocol will adjust the guess of the PSF to the actual PSF that blurred the synthetic measurements.
We run a total of 2 Â 10 5 AU iterations to obtain the reconstructions: the protocols are executed in steps of 50 iterations for both the object and the PSF, and the outer blind cycle is repeated 2000 times.
With this choice, the object and the PSF are processed for 10 5 iterations each. The result of the object reconstruction is shown in the bottom-right triangle of Figure 3b. The reconstructed image exhibits higher contrast and enhanced void regions, even when not directly discernible in the blurred version. We use intensity profiles to quantify the resolution increase by comparing the ground truth, the initial guess, and the de-autocorrelated and deconvolved image (Figure 3f).
We observe the profiles along the red line drawn in Figure 3a. Starting from a hump, in which the opposite walls are indistinguishable, the two peaks of the original image are recovered. The deconvolution process transfers the signal energy from the center to the two opposite lobes. Moreover, the recovered PSF reaches a remarkable result: if we start from an isotropic Gaussian function, the kernel quickly converges toward a star-like shape (Figure 3g). For completeness, we report that the cross takes form after a few hundred iterations, and its overall shape continues to refine during the following steps.
The image quality of the reconstruction and the convergence trend of the algorithm can be quantified with the metrics defined in

| Multiview blind deconvolution of zebrafish samples
After validating the method with synthetic samples, we tested it experimentally in a multiview light-sheet microscope. We report the  Figure 4b. For the Blind-AU, we compute the average autocorrelation of each view, as described in the previous section. We recall that this step does not require any alignment process. Then, we feed this estimation to the algorithm described in Equations (1) and (2) and schematized in Figure 1. After running a total of 25 Â 10 4 iterations, we obtain the result displayed in Figure 4c. Our reconstruction is sharper and highly resolved compared to the standard fused data and to the deconvolved one (Figure 4a,b). In particular, we can appreciate the improvement along the tomographic plane (bottom panels of Figure 4), where our method neatly discriminates the two vessels located in opposite regions of the spine. To quantitatively assess the improvement, we examine the profiles of the three volumes along the direction indicated by the red line in Figure 4a. The plot of the profiles is shown in Figure 4d. As expected, the RLD gives results sharper than the initial guess. Our method further surpasses this result by better separating the opposite vessels, reconstructing details sharper than standard approaches. Respectively for RLD and Blind-AU, the convergence MSE is 4:06 Â 10 À16 a:u: and 3:20 Â 10 À15 a:u. The image quality of the direct measurement is quantified by the image SNR and is equal to 63:55 dB. RLD brings the metric of the reconstruction to 72:01 dB, which is further improved by the Blind-AU algorithm to 79:93 dB.
Having assessed the image reconstruction ability, we now focus on the quality of the PSF obtained at the end of the process. We take a slice along the tomographic plane obtained by fusing orthogonal views (Figure 5a). We first average the autocorrelations of all the views. The Blind-AU is initialized with the fused object reconstruction and a completely random guess for the PSF (Figure 5b). The object reconstruction and PSF are optimized alternatively with our method.
After running 5 Â 10 4 iterations, we obtain the reconstructions of the object displayed in Figure 5c, and of the PSF in Figure 5d. graph is plotted in Figure 5f. Our method blindly recovers a PSF close to the measured one, guarantying a resolution increase even when completely ignoring the optical response of the system. Quantitatively, the Blind-AU algorithm improves the image quality of the sample also in the case of a random PSF initial guess: starting from an image SNR of 64:47 dB the value is increased to 76:77 dB after the iterative process. In this case, we assume convergence when the MSE is 1:42 Â 10 À10 a:u.  (Lucy, 1974;Richardson, 1972). On the other hand, additive noise accumulates in autocorrelation space as a delta ing it with the measured optical response along the scanning axis. The setup used in our experiment was a custom-made light-sheet fluorescence microscope with a rotating stage but, our technique is readily usable in any omnidirectional LSFM setups , or in techniques directly based on autocorrelation inversion such as hidden tomography (Ancora et al., 2020).

| CONCLUSIONS
Although it guarantees excellent reconstruction capabilities, our method may have room for further improvements in several aspects.
Concerning the multiview aspect, the objects that appear sharp in one view may display less contrast compared to the view belonging to the opposite sides of the specimen. Even assuming the scattering to be negligible, the presence of tissue absorption weakens the signal coming from deep inside the tissue and makes the fusion problem nontrivial. In turn, this implies that to consider the average of a view may not be the ideal choice to obtain the best reconstruction even if the datasets are perfectly aligned. If the datasets are misaligned (as it typically is), the problem gets even worse and binds to the fusion problem: the optimal fusion does depend upon prior alignment. As thoroughly discussed, our method avoids the alignment procedure by using autocorrelations, and we fuse the views by taking their average. A further in-depth study on the assessment of the optimal choice for fusion might indeed be facilitated with our method, since fusion does not rely on alignment anymore. Our intention is to tackle this problem in follow-up studies, thus, pushing the quality of autocorrelation-based tomography further.  (Guo et al., 2020). This idea is suitable for our implementation and can lead to faster solutions. Since our formulation is structurally similar to RLD, regularizing the solution is an attainable refinement.
Total variation (Laasmaa et al., 2011) and sparsity constrained (Shaked et al., 2011) are two choices we plan to incorporate in our formalism to further enhance the reconstructions. Including priors could also enhance the performance of our method. For example, in multiview geometry, the angle of rotation is generally known.
Designing a functional form for the PSF that accounts for the angular rotation could ease the reconstruction process. Lastly, our algorithm assumes uniform PSF blurring thorough the entire imaged volume. Although this is often considered a good approximation, higher image quality can be reached when considering spatially anisoplanatic deconvolutions (Toader et al., 2021) or mixed optical aberration corrections (Furieri et al., 2022). These considerations fall beyond the scope of the present manuscript, which may open new paths for further studies.

| Synthetic sample generation and measurement
To test the Blind-AU algorithm, we generated an image simulating blood vessels. The simulated multi-view measurement is composed of two orthogonal views of the sample, generated by convolving the original image with an elongated Gaussian PSF (with standard deviation σ z = 2.7 px, σ x = 1 px), oriented along two different axes (θ 1 ¼ 0 and θ 2 ¼ 90 ). Then, images are normalized to 2 16 and Poisson noise is added, with average and variance equal to λ ¼ 2 8 . For N different views, identified by the subscript θ i , each simulated measurement is determined by: Where h θi is the PSF rotated by the angle θ i , and ϵ θi is the Poisson

| Light-sheet fluorescence microscopy
Light-sheet fluorescence microscopy is an optical technique that exploits decoupled illumination and detection, placed orthogonal to each other, to achieve optical sectioning. The illumination arm generates a light sheet (a single plane of illumination) on the sample. The emitted fluorescence is collected and filtered along the detection path, and forms an image on a widefield detector (CMOS camera).
Having independent illumination and detection, the PSF is determined by the product between illumination and detection PSFs and is elongated, typically, along the axial direction (Power & Huisken, 2017).
The system starts with a single-mode blue laser emitting at 473 nm (MBL-FN-473, 50 mW power), operating in CW mode. The laser beam is collimated to a diameter of 4 mm (Thorlabs collimator RC04FC-P01), and its width is reduced to 1 mm by a vertical slit. Then, the light is sent to a galvanometric mirror (Thorlabs GVS001) which is conjugated with the sample position along the optical illumination path.
The galvanometric mirror is controlled by a sinusoidal voltage signal and oscillates with an amplitude of AE1 and frequency of 250 Hz. This pivots the illumination beam to remove shadowing artifacts that are common for side illumination of the sample (Di Battista et al., 2019).
Then, the beam is expanded by a factor of 3 by a telescope This value is also the lateral dimension of the region imaged by every single pixel, by applying a 2 Â 2 binning to the camera.

| PSF characterization
The point-spread function (PSF), representing the intensity impulse response of the microscope, describes the image degradation process.
The PSF was estimated experimentally by measuring a sample composed of sparse fluorescent nanobeads, much smaller than the microscope resolution. They act as individual point sources, whose recorded image describes the impulse response of the system. Beads should be measured under the same experimental conditions as the sample, to avoid variations due to wavelength changes, and retrieve information about the alignment of the optical elements (Sibarita, 2005). We used fluorescent nanobeads (Estapor, XC dye, 100 nm size) absorbing between 460 and 500 nm and emitting in the 515 À 570 nm range. They are embedded in phytogel with a concentration of 1 : 10 5 and inserted in a cylindrical FEP tube with an internal diameter of 2 mm and an external one of 4 mm. To characterize the 3D PSF, we acquired an entire volume of beads at the same experimental conditions as the sample; then, we selected one fluorescent bead located in the center of the FOV and applied a Gaussian fit along the three dimensions to retrieve an analytical expression. The groundtruth PSF of the 2D reconstruction ( Figure 5) (Kaufmann et al., 2012). The tube is mounted in the microscope chamber, filled with a stock solution composed of fish water with 0:016% tricaine concentration.

| Richardson-Lucy deconvolution
Richardson-Lucy deconvolution (RLD) algorithm is an iterative numerical method for image restoration, stable in presence of high noise levels (Lucy, 1974;Richardson, 1972). It is derived by interpreting images and PSFs as probabilities and applying Bayes theorem. In microscopy, the system's point spread function, h, is not ideal; thus, the measurement f μ of an object o can be described by the convolu- In this description, we can state that pixel-value of f μ is determined by the re-distribution of the energy of o accordingly to h: Assuming that the PSF is shift-invariant (isoplanatic condition), the RLD iterative process can be written in a convolutional form: With the notation e h, we indicate the PSF with reversed coordinates: e h ¼ h Àx ð Þ, and the index t identifies iterations. Starting from an initial guess o 0 ≥ 0, nonnegativity is preserved during the iteration process, without the need for positivity constraints (Fish et al., 1995).

| Anchor-update deconvolved deautocorrelation
Our blind method strongly relies on the Anchor-update (AU) algorithm  that, here, we recall for completeness. The goal is to recover a sharp object o rec from the measurement of the autocorrelation χ μ blurred by a given function H so that: Since in our specific case, H ¼ h ? h. The solution with respect to o rec can be approached via an iterative fixed-point multiplicative method similar to RL deconvolution: where K t ¼ o t ? H is an effective kernel continuously updated as the reconstruction progresses. Equation (6) defines the AU iteration as a function of the iteration step t.

| Metrics to assess the image quality
We monitored the convergence trend by computing the mean squared error (MSE) that estimates the similarity between the final reconstruction blurred by the PSF of the system and the initial measurement: Here we denote with o rec the reconstruction obtained with the iterative algorithm, with h being the PSF of the system, and with s μ the initial measurement. If this metric is applied to the RLD algorithm, h is the known PSF of the optical system. Instead, when considering the AU algorithm, h is the PSF estimated by the blind iterative process. The measurement s μ and the blurred image o rec Ã h are energynormalized.
On the other hand, we assess the image quality by calculating the signal-to-noise ratio (SNR) of the original measurements and the reconstructions, according to the expression: where s max is the peak intensity of the image, and ϵ is mean value of the background noise estimated by in a region of the image where the sample is not present.

| Preprocessing of experimental data
In the present manuscript, a multiview measurement is composed of four orthogonal views, acquired by rotating the sample by steps of 90 .
First, we subtract the background signal from each matrix by measuring a full stack without laser illumination. Second, we rotate each stack against the reference view: the passage is performed in multiples of 90 .
The measure passed to the algorithm lies in the autocorrelation space, and it is computed by autocorrelating each stack and averaging them. It guarantees intrinsic subpixel alignment . Here, autocorrelation is carried out as multiplication in the Fourier domain, according to the cross-correlation theorem. To avoid numerical errors before the fusion, we take the absolute value of each autocorrelation.
To determine an initial guess for the sample image, we decided to align the views. This step is carried out by setting the first view as a reference, then computing the cross-correlation between this stack and the others.
For each couple, the position of the cross-correlation maximum represents the translational shift to apply in the three spatial coordinates.
Finally, the initial guess is computed as the average of all the aligned views. In addition, we removed all the null elements by substituting a small value, taken as 1/100 of the image average. Otherwise, those pixels will remain zero during deconvolution. It is worth noticing that the alignment procedure is not strictly necessary to run the algorithm but ensures convergence in case the object is truncated in the image plane.