## 1.

## Introduction

Modern imaging techniques, including two-photon microscopy, produce *in vivo* experimental data of great diversity and volume. New experimental possibilities, such as optogenetics, bring about new scientific questions to address and consequently generate a torrent of data. One of the existing problems is the compromise between a demand for fast sampling to resolve rapid transient processes on one hand and sufficient spatial resolution on the other hand. Although imaging data tend to be collected in high detail, a large portion of it contains background or noise and is discarded downstream in the processing pipeline, while a small fraction is considered the useful “signal.” It is important to robustly identify the useful parts of data in low signal-to-noise ratio (SNR) measurements.

An important feature of analyzed data is its sparseness or compressibility under some transform, i.e., if the data can be decomposed into some basis such as Fourier or wavelet, where only a few coefficients are significantly larger than zero. The notion of sparseness is essential for noise suppression. If transformed data are sparse, one can legitimately assume that only the few large coefficients contain information about the underlying signal, whereas the small-valued coefficients can be attributed to noise. This leads to the idea of thresholding in a transform space for image enhancement.^{1}^{,}^{2} Because wavelet transform provides a nearly sparse representation of piecewise smooth images,^{3} this transform and its variants have become popular in image denoising and signal detection in transform domain,^{4}5.6.7.^{–}^{8} although there are also other transforms that can provide sparse representation, such as singular value decomposition and Radon transform. In short, wavelet transform is a multiscale representation of the input data, obtained via iterative application of band-pass filters. Wavelet coefficients capture the signal features at different locations and hierarchical spatial resolutions. The multiscale property is particularly useful for denoising because a typical useful signal is sparse (concentrated in a few coefficients at several scales), while noise is homogeneously distributed.

Calcium signaling in astrocytes takes versatile forms, one of which is spatially patterned spreading ${\mathrm{Ca}}^{2+}$ signals that pervade astroglial networks. Intercellular ${\mathrm{Ca}}^{2+}$ waves in astrocytic syncytium can be triggered by electrical and mechanical stimulations,^{9} local elevation of extracellular adenosine triphosphate (ATP) level,^{10} or by neuronal activity *in situ*^{11} Spontaneous glial calcium signaling is reported to guide axonal growth and cell migration in the developing brain,^{12}13.^{–}^{14} and calcium waves may represent a reaction to local tissue damage or other pathology. For instance, the incidence of spontaneous ${\mathrm{Ca}}^{2+}$ waves is increased with aging and low-oxygen conditions.^{15} Detection and reconstruction of ${\mathrm{Ca}}^{2+}$ waves in fluorescent imaging recordings pose challenges for data analysis. The key problem is identifying transient low contrast events in large series of images at a low SNR. Conventional data analysis methods can be loosely categorized to region of interest (ROI) type analyses, pixel thresholding, statistical component analyses, and multiscale (usually wavelet based) analyses. ROI analysis and pixel thresholding work particularly well with evoked responses, easily recognizable cellular correlates, relatively low noise, and small datasets. An ROI-based approach is primarily popular because it is often simple to select the ROIs corresponding to neuronal bodies and then process the ROI-averaged fluorescence traces. However, it becomes unwieldy for the analysis of sparse spontaneous events in large datasets and high noise levels or where *a priori* selection of regions of interest is impossible, which is the case for ${\mathrm{Ca}}^{2+}$ waves in the molecular layer of the cerebellum. Independent component analysis (ICA)^{16} is capable of processing large datasets with sparse spontaneous events, but has some limitations. Specifically, the output of ICA relies on the independence of the analyzed signals, and does not preserve the relative amplitude or sign of the detected components; in application to frame series it does not directly take advantage of local correlations in pixel intensities. Spontaneously occurring glial ${\mathrm{Ca}}^{2+}$ waves (GCWs) are hard to detect and quantify by ROI analysis, pixel thresholding, or ICA.

In this paper, we present an example of a multiscale approach to image and signal processing. We discuss a framework for the multiscale vision model (MVM)^{17}^{,}^{18} in application to the problem of detection and reconstruction of spontaneous intercellular glial ${\mathrm{Ca}}^{2+}$ waves. This technique relies on the sparsifying property of the wavelet transform. In short, wavelet transform is a multiscale representation of the input data, made via iterative application of band-pass filters. Wavelet coefficients capture the signal features at different locations and hierarchical spatial resolutions. Due to the nature of the used transform, MVM is most suitable for detection of nearly round faint structures in noisy images or bright blobs in three-dimensional (3-D) data. Because glial calcium waves are manifested as expanding nearly isotropic transient fluorescent elevations, MVM performs well in pursuing such structures frame after frame, and then stitching overlapping objects as snapshots of a single GCW event.

## 2.

## Multiscale Vision Model for Detection of Glial Calcium Waves

Recently, we adapted a wavelet-based framework, multiscale vision model to detect and analyze the spontaneous intercellular calcium waves in mouse cerebellum glial cells *in vivo*.^{18} In short, this framework amounts to noise rejection by thresholding wavelet coefficients of the input image, followed by establishing an interscale relationship between significant coefficients at different scales of decomposition and partial iterative reconstruction of the detected objects. Originally, this framework was suggested for the analysis of astronomy images^{19} and employed redundant wavelet transform with B-spline basis, also called “starlet” transform.^{2} Here, we report two improvements to our previous MVM implementation:^{18} (1) we increase the tolerance to outlier samples by using a mixed multiscale median and the starlet transform, and (2) we implement an iterative partial reconstruction algorithm. To illustrate this improved framework, we use several phantom patterns [Fig. 1(a)] of different shapes, contaminated with additive Gaussian white noise (AGWN) with a few bright “hot-spot” pixels with $50\times $ the noise variance in 0.05% of pixels [Fig. 1(b)]. This can also be described as a small amount of “salt noise” (SN). The framework is summarized in Fig. 4 using an experimental 2-pm record with a spontaneous glial ${\mathrm{Ca}}^{2+}$ wave in the mouse cerebellum as an example. The MVM framework is compared to other denoising techniques in Fig. 6 for phantoms and in Fig. 7 for the experimental recording.

## 2.1.

### Starlet and Multiscale Median Transforms

Starlet transfrom (also known as *à trous* transform) decomposes an original image $I(x,y)$ into a set $\{{w}_{j}(x,y)\}$ representing two-dimensional (2-D) image details at different scales $j$ (wavelet coefficients) and a smoothed approximation ${c}_{N}(x,y)$ at the largest scale:

The low-pass filter is $2\times $ zero-upsampled at each level, leading to interlaced image convolution. In the starlet transform, a discrete filter based on a cubic $B$-spline is used. In one dimension (1-D), this filter takes the form ${h}_{1}=\phantom{\rule{0ex}{0ex}}[(1/16),(1/4),(3/8),(1/4),(1/16)]$. For 2-D images, one can use the outer product of the two one-dimensional filters ${h}_{2}={h}_{1}\otimes {h}_{1}$ or process each dimension separately with a 1-D filter.

A known drawback of the starlet transform is that a very bright point structure will have responses in many scales instead of just being captured by the finest scale.^{2} Such small bright structures can be outliers, like occasional “hot pixels” of the charge coupled device matrix or salt noise. This is illustrated in the starlet decomposition of the phantom image, mixed with AGWN and salt noise [Fig. 1(b)], where the bright outlier points influence coefficients up to the fourth level of decomposition. Such a problem does not exist in the multiscale median transform [Fig. 1(c)]. This transform is organized in the same iterative way as the starlet transform, but the role of low-pass filtering is played by the median filter, with the window size doubling at each level of decomposition. In 2-Ds, the median filter of an image $f$ with window size $L\times L$ can be defined as $\mathrm{Med}(f,L)\text{:}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{f}_{i,k}\to {Q}_{2}\{{f}_{i\pm L,k\pm L}\}$, i.e., each pixel in the image $f$ is replaced by the sample median of pixel values in the $L\times L$ neigborhood of the pixel ${f}_{i,k}$. The window size changes with the decomposition level $j$ as $L=4j+1$. The median filter is nonlinear and provides for robust smoothing, i.e., it attenuates the effect of outlier samples.^{2} Hence, one of the advantages of the multiscale median transform is better separation of structures in the scales, but the starlet transform provides a more robust noise estimation at different scales. In a compromise, one can mix the two transforms and use the median or starlet transform depending on the amplitude of the coefficient. In the original merged median/starlet transform (MST)^{2} at level $j+1$, one first computes the median filtering of the approximation ${c}_{j}$ with a window size $4j+1$:

## 2.2.

### Significant Wavelet Coefficients

Detection of structures of interest that are significantly brighter than the background should be based on the knowledge about the statistical distribution of the wavelet coefficients $\{{w}_{j}(x,y)\}$ in the background. We assumed stationary Gaussian white noise to set significance levels. Statistics of wavelet coefficients at each level were estimated with Monte Carlo simulations: we obtained standard deviations ${\sigma}_{1}(j)$ of noise wavelet coefficients at each level of decomposition $j$ (index 1 in ${\sigma}_{1}$ represents unit variance). The obtained ${\sigma}_{1}(j)$ values were used as weighting factors for arbitrary noise variance in the analyzed images. More details can be found in Ref. 17 and in Chapter 2 of Ref. 20

The knowledge of noise standard deviations at different spatial scales allows the definition of significant coefficients for an image $I(k,l)$ by thresholding at $k\sigma (j)$ at each level. If ${W}_{j}(x,y)>k\sigma (j)$, the coefficient is considered significant and possibly belonging to a bright object. Interested in elevations of ${[{\mathrm{Ca}}^{2+}]}_{i}$, we only test positive coefficients corresponding to the luminous sources. The choice of $k$ can be varied for optimal performance; in the figures below, we used $k=3.3$ which corresponds to the 99.95th percentile of a Gaussian distribution.

## 2.3.

### Interscale Relationship and Object Reconstruction

Image $I(x,y)$ can be modeled as a composition of ${N}_{o}$ objects ${O}_{i}$, smooth background $B$, and noise $N$:

To recover the objects ${O}_{i}(x,y)$ in a given image, we use contiguous regions of significant coefficients at each scale and establish their interscale connectivity relationships. Let a structure ${S}_{j,l}$ be a set of $p$ connected significant wavelet coefficients at scale $j$:

## (5)

$${S}_{j,l}=\{{w}_{j}({x}_{1},{y}_{1}),{w}_{j}({x}_{2},{y}_{2}),\dots ,{w}_{j}({x}_{p},{y}_{p})\},$$When there are several close objects or a smaller object is overlaid on a bigger one, the objects can become entangled in one connectivity tree and should be deblended. A structure ${S}_{j,k}$ will be detached from the tree as a new root node if there exists at least one other structure at the same level belonging to the same tree and the following condition is fulfilled: ${w}_{j-1}^{m}<{w}_{j}^{m}>{w}_{j+1}^{m}$. Here, ${w}_{j}^{m}$ is the maximum wavelet coefficient of ${S}_{j,k}$; ${w}_{j-1}^{m}=\mathrm{max}\{{S}_{j-1,l}\}$, where ${S}_{j-1,l}$ is the structure connected to ${S}_{j,k}$, such that the position of its maximum wavelet coefficient is closest to the position of the maximum of ${S}_{j,k}$, if ${S}_{j,k}$ is not connected to any structure at scale $j-1$ then ${w}_{j-1}^{m}=0$; and ${w}_{j+1}^{m}=\mathrm{max}\text{\hspace{0.17em}}{w}_{j+1}(x,y)|{w}_{j}(x,y)\in {S}_{j,k}$, i.e., the maximum wavelet coefficient at scale $j+1$, such that its location belongs to ${S}_{j,k}$. Recursive application of this deblending procedure to the connectivity tree yields objects as independent structures of significant wavelet coefficients. Partial reconstruction of these objects as images is a nontrivial task. The simplest solution is to perform inverse wavelet transforms for each object, setting all wavelet coefficients not belonging to the object to zero.

In our previos study,^{18} we looked for a compromise between the computational speed and accuracy of reconstruction. For the sake of computational speed, we reconstructed objects simply as the inverse transform ${R}_{w}$ of the wavelet coefficients representing the object. Here, we improve our framework with an iterative reconstruction scheme.^{2} The idea is that we search for an image that will, under the multiscale transform (e.g., starlet), result in coefficients as close as possible to those to ${O}_{i}^{w}$. Thus, we want to minimize $E={\Vert {\mathbf{M}}_{i}({O}_{i}^{w}-{\mathbf{T}}_{w}X)\Vert}_{2}^{2}$ by varying reconstruction image $X$, where ${\mathbf{T}}_{w}$ is the wavelet transform operator, and ${\mathbf{M}}_{i}$ is the multiresolution support of ${O}_{i}^{w}$, i.e., has ones where ${O}_{i}^{w}$ is nonzero and has zeros otherwise. The solution is searched iteratively:

## (7)

$${X}^{n+1}={X}^{n}+\alpha {\mathbf{R}}_{w}{\mathbf{M}}_{i}({O}_{i}^{w}-{\mathbf{T}}_{w}{X}^{n}),$$The reconstruction results for the three decomposition/reconstruction pairs are shown in Fig. 3(d). We should note that in the reconstruction images, each of the pattern was independently reconstructed, but all the individual object reconstructions are summed in one image for the clarity of representation. The starlet-starlet variant retains the point structures along with the patterns, while the schemes involving the median/starlet transform reject this type of noise. Although the MST-MST scheme enhances the object contrast, the residual image (input minus the reconstruction) shows that it does not preserve initial object intensity, while the residual image for the MST-starlet scheme almost exclusively contains noise and little if any traces of the original patterns. This, together with a lower computational load, suggests the MST-starlet decomposition/reconstruction pair as the scheme of choice for the task of glial ${\mathrm{Ca}}^{2+}$ waves detection and reconstruction.

The full MVM framework with the mixed median/starlet decomposition and iterative reconstruction to detect and recover spontaneous glial calcium waves is summarized in Fig. 4. For each frame of normalized imaging data, we perform the median/starlet decomposition and after thresholding find contiguous areas of significant coefficients. These structrures are linked into interscale tree-like connectivity trees (connected structures belonging to one object are shown in the same color in the figure). Based on these connectivity trees, iterative reconstruction using the starlet transform is individually done for each tree according to Eq. (7) (shown in different colors as “significant structures in the ${k}_{\mathrm{th}}$ frame” in the figure). After this the procedure is done for all frames, and the individual objects contained in each frame are stored. The overlapping objects in the neighboring frames are linked into an $XYT$ 3-D representation of an evolution of distinct ${\mathrm{Ca}}^{2+}$ signaling events, shown as 3-D colored blobs in the scheme. Clearly, not all ${\mathrm{Ca}}^{2+}$ signaling events are glial ${\mathrm{Ca}}^{2+}$ waves, but after detection and 3-D reconstruction stages, one can identify the events of different types manually or according to any automated sorting scheme.

## 2.4.

### Comparison to Other Denoising Techniques

We compared the MVM framework with our modifications to other modern denoising techniques, although an exhaustive survey of different denoising and object detection techniques is beyond the scope of this study. Figure 5(a) presents how the PSNR of denoised images changes with the PSNR of noise-corrupted images. Here, we compare the MVM, two algorithms of total variation (TV) denoising, which searches for image approximation with a minimal TV (integral of absolute gradient) and a bilateral filter.^{21} The TV denoising was used in two variants: the Chambolle algorithm^{22} and the split-Bregman algorithm.^{23} The TV and bilateral filtering implementations were used from the open-sourse scikit-image python library.^{24} Because the maximum value of the pure phantom image was 1, the PSNR was calculated as $20\text{\hspace{0.17em}}\mathrm{lg}(1/\sqrt{\mathrm{MSE}})$, where MSE is the mean square error between the tested image $f$ and the pure phantom image $g$, each of size $N\times N$: $\mathrm{MSE}=\phantom{\rule{0ex}{0ex}}(1/{N}^{2})\sum _{k,l}{\Vert f(k,l)-g(k,l)\Vert}^{2}$. It is clear that with AGWN, the MVM resulted in a better PSNR of the recovered images at very low SNRs of the input images and had a similar performance [Fig. 5(a), left pane] to the TV/Chambolle algorithm at a higher PSNR of the input images, although MVM displayed much a higher performance scatter than the other algorithms at the higher PSNR of the input images. The better performance of the MVM at a very low PNSR of the noisy images likely results from background rejection in the MVM and exclusive reconstruction of only the detected structures. The addition of small amounts of salt noise (outlier pixels) did not affect the MVM performance, but deteriorated the performances of the other algorithms [Fig. 5(a), right pane]. It is clear that while the used denoising algorithms are not optimized for outlier measurements, preprocessing of the input images with filters, devoted for removal of salt noise, e.g., adaptive median filter, would have reverted the results to the case of a simple AGWN. Figure 5(b) visualizes the difference between denoising provided by the TV/Chambolle and the MVM approaches at different PSNR values for the degraded phantom image. Comparsion of the images denoised by the two algorithms supports the observation that the higher PSNR values achieved by MVM are due to setting the background to zero.

An important feature of the MVM is the ability to detect and segment structures in an input images, which is built in. Figure 6(a) presents the results of denoising of the test noisy phantom with TV/Chambolle, TV/Bregman, and the sure-let wavelet-based algorithm.^{25} SURE-LET denoising was done using the MATLAB^{®} code provided by Luisier and colleagues.^{26} To detect individual objects in the output of the alternative denoising algorithms, we performed Otsu thresholding^{27} and labeled the resulting contiguous areas as individual objects. Because the given denoising techniques did not effectively suppress the salt noise present in the input phantom, we had to perform adaptive median filtering on the denoised images prior to Otsu segmentation. The results are presented in Fig. 6(b), where each object is shown in a different color. The object separation in MVM matches the ground truth of the object composition, while the alternative methods tend to segment some of the patterns into multiple objects. Some false-positive objects are also detected. It is likely that it is possible to tune the thresholding and segmentation algorithms for the other denoising methods to produce more satisfactory results, but in our view, the power of MVM is that this object separation comes out of the box. The TV/Chambolle produced the best results among the alternative denoising approaches and it was chosen for further comparison of the algorithm performance.

Figure 7 shows the results of the TV/Chambolle denoising and the MVM for the experimental data, displaying several glial ${\mathrm{Ca}}^{2+}$ waves co-occurring in the field of view (FOV). In short, we imaged spontaneous glial ${\mathrm{Ca}}^{2+}$ waves in mouse cerebellum *in vivo* under ketamine anesthesia using two-photon microscope and Oregon green BAPTA-1/AM dye. In the top row, every fourth frame of the original normalized frame sequence is shown (time interval between the shown frames is $\approx 2.9\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$). Fluorescence values in each pixel were normalized to their respective standard deviation after subtraction of the mean value. Experimental details are given in Refs. 15 and 18. Because we do not have ground truth for the events in the experimental data, we compared the denoising and detection results to operator-detected events. Locations of the five glial ${\mathrm{Ca}}^{2+}$ waves in the FOV as detected by the operator in the $\mathrm{\Delta}F/{\sigma}_{F}$ data are shown as circles in the fourth column of images. The TV/Chambolle denoising followed by Otsu thresholding and zeroing pixel values below the threshold (middle row) clearly improves the contrast of the images; however, there are also numerous above-threshold pixels which do not belong to any wave visible in the FOV. In contrast, the MVM (bottom row) provides a cleaner view of the data, recovering areas of elevated ${\mathrm{Ca}}^{2+}$ corresponding to glial ${\mathrm{Ca}}^{2+}$ waves, while at the same time rejecting background and providing for recognition of individual objects.

## 3.

## Conclusion

Our aim in this work was to illustrate the utility of multiscale transforms for the extraction of useful signals from two-photon laser scanning microscopy imaging data. Our main point of application was detection and reconstruction of spontaneous glial ${\mathrm{Ca}}^{2+}$ waves using the multiscale vision model approach. In comparison to our previously reported MVM experience,^{18} we extended the technique by using the more robust mixed multiscale median/startlet transform for image decomposition and implemented an iterative scheme for nonlinear image restoration from the wavelet coefficients representing an object of interest.

With this technique we were able to detect and reconstruct areas of increased ${\mathrm{Ca}}^{2+}$ elevations, corresponding to glial calcium waves and other forms of ${\mathrm{Ca}}^{2+}$ signaling in each of the recorded image frames. Performed in all frames, this procedure allows us to recover the whole dynamics of individual calcium waves by linking overlapping objects together in the neighboring frames. The mixed multiscale median/starlet transform allowed for suppression of both the Gaussian noise and outlying samples such as salt noise, while iterative object reconstruction from the significant wavelet coefficients and multiscale support allowed for a nearly perfect representation of the patterns, although the difference between the source image and the object reconstruction contained almost nothing but noise.

With only minor modifications, MVM should be extendable to other applications, such as quantification of blood vessel cross-sectional areas, since penetrating vessels also appear as nearly circular structures in typical two-photon brain imaging data. There is obviously more to neuroimaging data than roundish calcium waves or vessel cross sections, and some native structures, such as Purkinje cell dendrites, are markedly anisotropic. Representation of such structures should benefit from using, for example, curvelet transforms.^{6} In general, it is probably a winning strategy to not limit the analysis to a single transform, but to combine the significant coefficients from several types of transforms, creating over-redundant dictionaries for morphological component analysis,^{2} which would allow for high flexibility in data representation.

## Acknowledgments

We thank Micael Lønstrup for providing excellent surgical assistance. This study was supported by the NORDEA Foundation/Center for Healthy Aging, the Lundbeck Foundation via the Lundbeck Foundation Center for Neurovascular Signaling (LUCENS), the NOVO-Nordisk Foundation, and the Danish Medical Research Council.

## References

^{2+}oscillations in astrocytes,” J. Neurosci. 27(33), 8957–8966 (2007).JNRSDS0270-6474 http://dx.doi.org/10.1523/JNEUROSCI.2276-07.2007 Google Scholar

## Biography

**Alexey Brazhe** is a senior researcher in the Department of Biophysics, Biological faculty, Moscow State University. His research interests include mathematical modeling and multiscale signal and image analysis tools in application to neuroscience. He received his PhD degree in 2006 from Moscow State University.

**Claus Mathiesen** is a postdoctoral fellow at Translational Neurobiology, University of Copenhagen. He received his PhD degree in 1990 from the University of Copenhagen and studies blood-brain-barrier mechanism, vascular and metabolic coupling to astrocytes, and neuronal activity in the aging brain.

**Barbara Lykke Lind** is a postdoctoral fellow in the Department of Neuroscience and Pharmacology at the University of Copenhagen, Denmark. She is working in Martin Lauritzen’s group studying intracellular Ca^{2+} activities in astrocytes and neurons and the implications for the neurovascular coupling. She completed her PhD degree in 2014 at the University of Copenhagen.

**Andrey Rubin** has been the head of the Department of Biophysics, Biological faculty, Moscow State University since 1976, and chair of the National Committee for Biophysics, Russian Academy of Science.

**Martin Lauritzen** is a professor of clinical neurophysiology and head of the Department of Clinical Neurophysiology at Glostrup Hospital, and professor of translational neurobiology at the Department of Neuroscience and Pharmacology at the University of Copenhagen, Denmark. His research interests focus on the mechanisms coupling neuroglial activity to cerebral blood flow and energy metabolism, the underlying mechanisms for brain imaging, cortical spreading depression, and mechanisms of aging.