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=[(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.² 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{:}{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.² 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)² at level $j+1$, one first computes the median filtering of the approximation $c_j$ with a window size $4j+1$:

Fig. 2

Illustration of the influence of outlying samples on detection of objects in one-dimensional (1-D) signals. From top to bottom: analyzed signal, containing white noise, a Gaussian with $\sigma =10$ at $t=512$ and an outlier with amplitude $100\times$ higher than the noise standard deviation at $t=485$; starlet coefficients at fifth level of decomposition; multiscale median transform coefficients at fifth level of decomposition; mixed median/starlet transform coefficients at fifth level of decomposition. The outlier distorts object of interest representation for starlet transform, but has no effect on multiscale median and mixed transforms.

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$:

Fig. 3

Iterative reconstruction of the detected structures. (a) Masks for significant (shown in red) starlet/wavelet coefficients for the noisy image shown in Fig. 1, (b) labeled interscale trees of the significant starlet/wavelet coefficients, (c) $L_2$ norm of difference between the significant coefficients and decomposition of the reconstructed image $E={\mathrm{\Sigma }}_i{\Vert {\mathbf{M}}_i(O_i-{\mathbf{T}}_{w}X)\Vert }_2^2$ for decomposition and reconstruction using the starlet transform (black circles), decomposition and reconstruction using the combined median/starlet transform (MST, blue diamonds), and decomposition using MST and reconstruction using the starlet transform (green triangles), (d) noisy image (+AGWN+SN), reconstructions using the different reconstruction variants shown in (c) and the residuals after subtracting the reconstructed image from the input image.

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=\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=\max 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,¹⁸ 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.² 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:

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.

Fig. 4

Scheme of the MVM protocol for a frame sequence: significant bright objects are detected and reconstructed in each frame using 2-D MVM, followed by relation of the objects in the neighboring frames and finally the time evolution of each detected ${\mathrm{Ca}}^{2+}$ signaling event is independently reconstructed.

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.²¹ The TV denoising was used in two variants: the Chambolle algorithm²² and the split-Bregman algorithm.²³ The TV and bilateral filtering implementations were used from the open-sourse scikit-image python library.²⁴ Because the maximum value of the pure phantom image was 1, the PSNR was calculated as $20\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}=(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.

Fig. 5

Performance of the MVM and other state-of-the-art denoising techniques at different PSNR levels of the degraded phantom image shown in Fig. 1(a). (a) Resulting PSNR for images, denoised with the MVM (green) two total variation (TV) techniques with Chambolle (blue) and Bregman (yellow) algorithms, and a bilateral filter algorithm (pink) as a function of PSNR of input (degraded) images; input images are mixed with AGWN (left pane) or AGWN+SN (right pane). Shaded areas around the curves are between 0.05 and 0.95 quantiles of the observed PNSR values and continuous lines are median values (measured after 200 runs), (b) visual comparison between the TV/Chambolle and the MVM denoising performance at different levels of input PSNR (after mixing with AGWN). Grayscale values rescaled to belong to $[\mathrm{0,1}]$ interval for each image independently.

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.²⁵ SURE-LET denoising was done using the MATLAB^® code provided by Luisier and colleagues.²⁶ To detect individual objects in the output of the alternative denoising algorithms, we performed Otsu thresholding²⁷ 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.

Fig. 6

Comparison of performance of MVM and other denoising techniques on the test noisy phantom, shown in Fig. 1(a). Denoising results for two total variation (TV) minimization techniqes with Chambolle and Bregman algorthms, bilateral filter and SURE-LET wavelet-based algorithm (a) and object segmentation and labeling results for the corresponding algorithms (b), segmentation is part of MVM, in other algorithms it was performed by labeling contiguous areas after binarization with automatic Otsu thresholding.

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\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.

Fig. 7

Normalized fluorescence data (top row), results of TV/Chambolle denoising followed by Otsu thresholding and rejecting pixels below threshold (middle row), and results of the MVM reconstruction (bottom row). Objects, overlapping in neighboring frames are attributed to one signaling event and are shown in the same color code. Shown is every fourth frame of the original sequence, same data as in Fig. 4, there are several co-occurring calcium waves in the field of view (FOV). In the MVM.