## Abstract

Studies concerning the physiological significance of Ca^{2+} sparks often depend on the detection and measurement of large populations of events in noisy microscopy images. Automated detection methods have been developed to quickly and objectively distinguish potential sparks from noise artifacts. However, previously described algorithms are not suited to the reliable detection of sparks in images where the local baseline fluorescence and noise properties can vary significantly, and risk introducing additional bias when applied to such data sets. Here, we describe a new, conceptually straightforward approach to spark detection in linescans that addresses this issue by combining variance stabilization with local baseline subtraction. We also show that in addition to greatly increasing the range of images in which sparks can be automatically detected, the use of a more accurate noise model enables our algorithm to achieve similar detection sensitivities with fewer false positives than previous approaches when applied both to synthetic and experimental data sets. We propose, therefore, that it might be a useful tool for improving the reliability and objectivity of spark analysis in general, and describe how it might be further optimized for specific applications.

- calcium sparks
- spark analysis
- myocytes
- confocal microscopy
- image processing

the study of brief, localized Ca^{2+} release events, such as Ca^{2+} sparks, using fluorescent indicators and confocal microscopy has afforded new insights into Ca^{2+} signaling in cardiac, skeletal, and smooth muscle (5). The task of analyzing these images is hampered, however, by the high degree of photon (shot) noise normally present. This is a consequence of the low number of photons collected per pixel at the high scan rates required.

Spark detection algorithms seek to speed up analysis and reduce the risk of bias inherent in manually distinguishing Ca^{2+} release events from noise, while identifying release events too small to be seen by eye (6, 12, 22). However, current algorithms were designed primarily only for images with a stable baseline intracellular Ca^{2+} concentration ([Ca^{2+}]_{i}). Here, we consider the effects of baseline variability on detection, and introduce a more general method of spark detection. We focus on linescan images, because these have received more attention in the literature and are better suited than full-frame (XYT) videos for the systematic investigation of the algorithm's detection capabilities. Nevertheless, the similarity of the problems encountered in processing both types of data sets means that the same methods can often be used with minor modification (e.g., Refs. 4, 7, 12, and 23).

Three distinct approaches that completely automate all steps of spark detection in linescans have been detailed in the literature. The most common, referred to as the “conventional algorithm” (4, 22), is most fully described by Cheng et al. (6). This combines median and mean filtering operations to improve the signal-to-noise ratio (SNR) in the image, before excising pixels >1.5 SD above the mean as potential sparks. Resting fluorescence and the noise standard deviation are estimated from the remaining pixels, followed by normalization and hysteresis thresholding to detect the sparks. This strategy is fast, capable of efficient analysis of large data sets, and has been made available both in an IDL implementation (6, 17) and as an ImageJ plugin (15). One important shortcoming is that the noise estimate is not robust: if, for example, no sparks are present, then the removal of the brightest pixels from the computation will lead to underestimating the noise (and potentially more false positive detections), whereas a high true spark density may have the opposite effect. This problem is somewhat overcome by the “à trous wavelet algorithm,” which uses the isotropic undecimated wavelet transform (IUWT) for robust noise estimation, denoising, and final detection, and thereby improves the sensitivity and reliability with which small sparks may be identified (22). The IUWT step combines filtering with subtraction to iteratively smooth an image and extract the difference between consecutive smoothed versions. The original image can be represented by the sum of a single highly smoothed representation and several wavelet (or detail) levels. Larger features are visible with greater contrast on the most highly smoothed wavelet levels. One can effectively isolate objects with properties that resemble sparks by thresholding appropriate wavelet levels. The two-dimensional (2D) IUWT is appropriate when the sparks appear approximately isotropic (rotationally symmetric) in the image space, which implies that the spatial and temporal spread is similar in terms of pixels (22). If this is not the case, the 1D IUWT might be more appropriate. This has previously been used to detect prolonged Ca^{2+} release events dubbed Ca^{2+}-embers (20). Finally, if the properties of the expected sparks are known approximately, the “matched filter detection algorithm” can be used to compute correlations between the image and a template spark (the matched filter) (12). By comparing the correlation coefficients with those obtained by correlating the template with noise generated by randomly reordering the pixels in the original linescan, appropriate thresholds can be determined.

Although each of these algorithms has been extensively tested and validated using synthetic images in which sparks are superimposed on constant backgrounds, none are well suited for processing more complex images in which the baseline fluorescence underlying sparks can change, for example, images acquired under voltage clamp (4, 6). This results from the assumption that the noise can be treated as signal-independent and characterized by the same standard deviation value everywhere. While it has been proposed to apply baseline subtraction before using the conventional algorithm for detecting sparks in such cases (4, 6), this does not take into consideration the signal-dependent nature of the photon noise present in fluorescence confocal images (13). This can cause the noise standard deviation to be several times larger in brighter, as compared with darker, parts of the image. Failing to take this variability into consideration is likely to reduce either the specificity or the sensitivity of spark detection, if not both.

The current study grew out of a need to efficiently and objectively analyze linescans of multiple smooth muscle cells in retinal arterioles. In these images, sparks may arise from the baseline or be superimposed on prolonged [Ca^{2+}]_{i} elevations, termed Ca^{2+}-oscillations (8). This represents a challenging test of spark detection in the presence of baseline and noise variability. After outlining and validating an appropriate general noise model, we introduce a new algorithm that extends the range of linescan images in which sparks can be automatically detected to include such data sets, while still being suitable for simpler images with constant baselines. By combining variance stabilization with local baseline subtraction, our algorithm effectively adapts detection thresholds automatically to local baseline and noise. This approach reduced false positives without compromising sensitivity. We propose, therefore, that the benefits of this strategy are not limited to complex cases but may extend to spark detection in any linescan.

## METHODS

### Test Images

#### Experimental images.

Linescan images in our laboratory were recorded using a confocal laser-scanning microscope (μRadiance, Bio-Rad) on an inverted microscope (Nikon Eclipse TE300; ×60 oil immersion objective, numerical aperture 1.4). Test solutions and tissues were excited in linescan mode (500 s^{−1}) at 488 nm; emitted light was filtered (530–560 nm) and detected using a photomultiplier tube (PMT). Data acquisition was controlled with Timecourse software (Lasersharp, Bio-Rad). Output images were 8-bit.

To verify the validity of our noise model and acquire spark-free images, we imaged different concentrations of Ca^{2+} (10 nM, 100 nM, 200 nM, 300 nM) and Fluo-4 (5 μM, 10 μM) in a solution containing 100 mM KCl, 10 mM MOPS, and 1 mM EGTA with pH 7.2 at 22°C for various permutations of iris size, gain setting, and laser power. This resulted in 96 linescan stacks, each consisting of 512 × 512 × 8 pixel frames, which were concatenated to form 512 × 4,096 pixel images.

We also obtained linescans of retinal arteriolar smooth muscle from Sprague-Dawley rats loaded with Fluo-4 AM (Molecular Probes). Microvessel segments between 0.2 and 2.0 mm in length, and with an outer diameter of 20–40 μm, were anchored with tungsten wire slips and superfused with bath solution (37°C) containing (in mM) 140 NaCl, 5 KCl, 5 d-glucose; 2 CaCl_{2}, 1.3 MgCl_{2}, and 10 HEPES, pH 7.4 with NaOH. The recording bath was rotated to align the long axis of the arteriole to the *x*-axis of the microscope, so that multiple cells were visible in each linescan. Animals were killed for tissue extraction under UK Home Office Certificate of Designation no. 5012. More detailed descriptions of the protocol are provided in Refs. 8, 16, and 21.

#### Synthetic images.

The algorithm's ability was tested using sparks of known properties embedded in synthetic noisy images (6, 12, 20, 22). Synthetic sparks had an exponential rise and decay and Gaussian spatial profile (11) with a rise time of 10 ms, half-time to decay of 20 ms, and full width at half-maximum (FWHM) of 2 μm (pixel size = 0.15 μm × 2 ms), as described for cardiac sparks (5). Amplitudes, defined as ΔF/F_{0} (F_{0} = the baseline fluorescence), were scaled at 0.05 intervals in the range 0.05–0.6, and randomly positioned in 512 × 2,048 pixel images under the condition that no sparks overlapped with one another or occurred at the image border. Poisson noise was added, and image sets were produced for SNRs of 2, 3, and 4, where the SNR is defined as the square root of F_{0} (12). For each combination of SNR and spark amplitude, five images were generated for spark frequencies of 20, 40, 60, and 80 sparks per image. This resulted in a total of 720 test images, including 36,000 sparks.

In addition, we generated a synthetic linescan resembling arteriolar smooth muscle by superimposing sparks on top of a varying background, which was computed by evaluating a tensor product spline.

### Noise in Confocal Linescans

#### Poisson and Gaussian noise models.

In fluorescence confocal imaging, the “true signal” *λ*_{i} can be viewed as the mean number of photons that would be counted at a certain pixel given sufficient repeated observations, while the recorded signal *f*_{i} is based on a single observation. The primary source of noise in the recorded image is expected to be photon noise, which is due to the inherent randomness of photon emission by the source and follows a Poisson distribution (6, 13, 14). A simple model of the signal is then
*X*_{i} ∼Pois.(*λ*_{i}) is a Poisson distributed random variable, α is a conversion factor (related to the gain of the imaging system), and μ is a constant offset (black level) that may be added to every pixel before digitization. A property of Poisson noise is that the noise variance is equal to the true signal, and the SNR can be defined as the square root of this value (13). However, it is important to note that pixel values in an image do not typically directly equate to photon counts, so that the standard deviation of the noise described by *Eq. 1* is *αλ*_{i}^{1/2}; i.e., it is a function of the unknown true signal and the (often) unknown conversion factor introduced during acquisition.

A consequence of this is that the noise varies along with the signal, meaning that the detectability of a spark is a function of the local baseline. Previously, Gaussian noise models have been widely used for spark detection (4, 6, 15, 22), in which the noise is treated as signal independent and characterized by a single standard deviation value. This simplifies detection by allowing the same tests to be applied everywhere and is justifiable when both the SNR is high (greater than about 3.5) and the underlying baseline is constant. However, the Poisson and Gaussian distributions differ appreciably when the SNR is low or if the baseline varies throughout the image.

#### Variance stabilization.

Variance stabilization transforms (VSTs) may be used to convert Poisson data into a form in which the noise may be treated as Gaussian (1). If one has access to the photon counts *X*_{i} (i.e., one has corrected for α and μ in *Eq. 1*), a simple stabilization method is to compute the square root of each count, after which the noise is asymptotically Gaussian with a standard deviation of 1/2 (2). The convergence of the noise can be slightly improved by first adding 3/8 to each count, resulting in the Anscombe transform (1). However, both of these methods become less reliable when the SNR < 5 (18), while the SNR of images used for spark detection is often in the range 2–4 (22). If one first smoothes the data using a filter *h* with a nonzero mean before computing square roots, it has been shown that effective variance stabilization of Poisson data may be achieved at much lower SNRs (around 0.5), assuming that the signal is constant under the support of the filter (24). Both optimal offset values (analogous to the 3/8 in the Anscombe transform) and the final asymptotic noise variance can be calculated from the filter coefficients (24).

This suggests that variance stabilization of filtered data is appropriate for confocal images of sparks. However, one seldom has access to the true photon counts. Determining the effects of gain, digitization, and other noise sources with sufficient accuracy to convert pixel values can be prohibitively difficult in practice. We therefore suggest a simplified form of the transform:

#### Validation of the noise model and variance stabilization technique.

Because a valid understanding of the noise is fundamental to detection, we explored the sufficiency of *Eq. 1* using confocal linescan images of Ca^{2+} solutions. The mean of all pixels in each linescan was taken as a measure of the signal, and the standard deviation of the difference between adjacent pixels, normalized by division by 2^{1/2}, was used as a measure of the noise. This was preferred over simply computing the standard deviation of all pixels because any small unevenness in illumination would inflate the raw standard deviation. To compare results across images, we corrected for μ and α by subtracting a black level value of 4.8 from all pixels, before dividing the result by 1/3 of the value of the gain setting used during image acquisition; these settings were appropriate for our system, for which varying the gain setting has a largely linear effect upon the recorded pixel intensities. A plot of the results is then shown in Fig. 1*A*, from which one can discern a clear relationship between signal and noise. The line *y* = *x*^{1/2} is also superimposed on the plot. This indicates where the data points would be expected to fall if the noise is Poisson and has a standard deviation equal to the square root of the signal.

Figure 1*B* then shows the effect of computing the square root transform after applying a 1D Gaussian filter (FWHM = 12 pixels). The horizontal scale is roughly equivalent to the SNR, since it is the square root of the mean normalized pixel value. Most points in Fig. 1*B* lie on or close to a horizontal line, indicating that the noise has been effectively stabilized, i.e., is independent of the signal. There is some variability, particularly at very low SNRs, probably due to imperfections in estimating μ and the influence of other noise components being more significant when the true signal is very low. However, the stabilization is effective when the SNR is greater than about 2. We therefore conclude that *Eq. 1* is a suitable model that describes the majority of the noise apparent in experimental images, and the square root transform of filtered data is appropriate for noise stabilization at the SNRs of interest in images containing Ca^{2+} sparks, as illustrated for a sample linescan from adjacent myocytes in a retinal arteriole (Fig. 2).

#### Additional noise features.

Two additional localized features of the noise were identified: *1*) isolated bright pixels occur with a frequency greater than would be expected based purely on the Poisson distribution (20), and *2*) there is some correlation between the values of spatially adjacent pixels, i.e., pixels scanned consecutively, even when no Ca^{2+} release event is present.

These observations are linked since several very bright pixels were often seen clustered along the direction of the scan line, as has been observed by others using PMTs (10). If necessary, extreme values may be removed quite easily using 3 × 3 median filtering. To not unduly affect the noise distribution, the median values should only be inserted into the image where bright pixels are found. These locations can be identified, for example, by computing the difference between the raw and median filtered images, and looking for pixels with values > 5× the standard deviation of the difference image.

The correlation between adjacent (consecutively scanned) pixels is a more insidious problem, which may not be initially apparent. Its origins are unclear, but such a relationship may arise as a consequence of the dark noise, minor laser fluctuations, or diffusion of fluorescent molecules when the signal is very weak. We confirmed its existence by calculating Pearson's linear correlation coefficients between adjacent rows and columns respectively of a frame of a representative Ca^{2+} concentration image (selected as having the SNR closest to the median). We found a correlation along the spatial dimension of 0.119 ± 0.046 (means ± SD), which was much greater than along the temporal dimension (0.002 ± 0.043). We conclude then that even when there is only noise present in the image, there can exist a correlation between the pixels adjacent to one another along the direction of scanning. This relationship is sufficiently small that it can normally be disregarded, but becomes important when using sophisticated techniques that assume the full statistical independence of the noise at each pixel in an image, because one risks underestimating the noise, increasing false positives. For this reason, noise is estimated by our algorithm while assuming the independence of the noise at adjacent pixels along the temporal dimension only.

### Spark Detection Using Variance Stabilization

Our spark detection algorithm combines the information from multiple spatial locations by means of a filtering operation, before stabilizing the noise variance using the square root transform and looking for brief fluorescence increases along the temporal dimension. Any background offset (μ in *Eq. 1*) should be subtracted from every pixel before application of the algorithm.

#### Preprocessing and variance stabilization.

Filtering is initially applied along the spatial dimension only, using a 1D Gaussian filter *h*_{s}. When the FWHM of the filter equals that of sparks, this is effectively a matched filter because of the approximately Gaussian spatial profile of sparks (11). If sparks are close together, or if the baseline varies appreciably at different spatial locations (e.g., because of the imaging of multiple cells), a smaller filter FWHM may be preferable to improve variance stabilization or avoid merging events. Each pixel is then replaced by the square root of its filtered value.

#### Noise estimation.

The noise distribution may now be considered Gaussian. Its standard deviation *σ*_{0} can be robustly estimated by taking the median of the absolute difference between all pairs of temporally adjacent pixels in the variance stabilized image, and dividing the result by 0.6745 × 2^{1/2}. This is equivalent to using the Haar filters and applying the median absolute deviation of wavelet coefficients estimation strategy described by Donoho (9).

#### Temporal filtering and baseline subtraction.

The baseline should now be subtracted before detection. A constant background value might be identified as the median of the data, or the iterative spline-fitting approach described by Bray et al. (4) could be used. However, to deal with the possibility of varying baselines while avoiding the high computational requirements of spline fitting, we propose instead to apply two 1D filters along the time dimension of the image: a smaller filter to reduce noise, and a larger filter to suppress brief release events (i.e., sparks) and approximate the baseline. The difference between both filtered images can then be thresholded for detection.

While Gaussian filters might again be used for this purpose, a more efficient algorithm can be defined in terms of the filters commonly used for the 1D IUWT (19). Applied to a signal *c*_{0} = *f*, increasingly smoothed versions are calculated iteratively by convolution with a filter *h*^{↑j}, i.e.,
*h* = [1, 4, 6, 4, 1]/16 is derived from the cubic B-spline, and *h*^{↑j} is the upsampled filter obtained by inserting 2^{j} − 1 zeros between each pair of adjacent coefficients of *h*. The smoothed image *c*_{j} is equivalent to the result of filtering *f* using a single, larger filter *h*^{j}, which may also be calculated iteratively using
*h*^{0} = 1. The iterative algorithm is, however, much faster, because each iteration requires only a filter with 5 non-zero coefficients, whereas *h*^{j} rapidly becomes large and dense as *j* increases. Also, for practical purposes, one must then only choose a number of filtering iterations *m* to apply for noise suppression, and a number of iterations *n* for baseline generation (*n* > *m*). The final image used for detection is the difference between both images, i.e.,
*n* = *m* + 1, then *d* is equivalent to a single set of wavelet coefficients generated using the 1D IUWT. Again, there exists a filter *h*_{t} that could be applied to *f* to generate *d* directly, which is

#### Thresholding.

The mean of *d* is 0, and all non-zero values are either the result of noise or the result of actual variations in fluorescence. If a pixel in *d* has a sufficiently large absolute value relative to the noise, this indicates that the underlying signal is neither constant nor linear under the support of *h*_{t}. Because of the use only of linear filters thus far, the noise distribution remains Gaussian. To determine an appropriate threshold to apply, it is desirable then to know the standard deviation *σ* of the values expected if only noise is present. This is estimated based on *σ*_{0} and the filter *h*_{t} using the equation
*τ* × *σ*, where *τ* is an adjustable threshold parameter. In some cases, particularly when using small filters to detect out-of-focus sparks, several regions might exceed the threshold but relate to the same underlying event. Also, above-threshold regions do not necessarily always contain the peak of the event to which they belong. Therefore it is desirable to expand these regions. This is accomplished in the conventional algorithm by increasing the size of each detected region until it is surrounded by pixels below a low threshold, set at 2 estimated noise standard deviations. This approach could also be used here, but we found that applied globally it could cause very large, in-focus sparks to produce excessively large regions. Therefore, we instead defined a lower threshold separately for each potential spark, so that each above-threshold region expanded until surrounded by pixels less than a user-definable percentage of its amplitude (where the amplitude is defined as the maximum pixel value found within the region). For the results here we used 40%, which ensured that detected regions corresponding to sparks were large enough to contain the peaks. These secondary thresholds were also clipped so that none could be higher than *τ*.

### Implementation and Tests

We implemented our algorithm and tests using MATLAB (R2010a, The Mathworks, Natick, MA) and its associated Image Processing Toolbox V7.0 (source code is provided in the Supplemental Material, which is available online at the Journal website). For ease of comparison, we also reimplemented the conventional and à trous algorithms in MATLAB, following the original descriptions and source code given by Cheng et al. (6) and v Wegner et al. (22). The conventional algorithm was used with threshold parameter *Cri* = 3.8, while the wavelet algorithm used denoising and threshold parameters *δ* = 4 and *τ* = 3.75, respectively, with the detection applied to wavelet levels 2–5. These values had been used in the original papers for similarly sized sparks.

Because all tested algorithms output binary images containing spark regions, a detection was defined as a “true positive” if it contained the peak of a real spark and a “false positive” otherwise; if multiple peaks occurred within the region, only one true detection was counted. This necessitated one additional change to be made to the à trous algorithm. A feature of this method is that multiple binary images are produced by thresholding each wavelet level, and previously this information had been combined to give spark regions either by an intersection (i.e., a pixel must be above the threshold on all wavelet levels) (22) or union operation (i.e., a pixel must be above the threshold on any wavelet level) (20). The former approach tends to cause the subdivision of detected sparks into multiple regions—most of which would be considered false positives by our metric—while the reduced requirements of the latter make false detections elsewhere more likely. To improve the fairness of the test, we therefore combined both approaches: spark regions were finally identified using a union operation, but kept only if they contained any pixels that were above the threshold on all wavelet levels. For all algorithms, any detected regions that overlapped image boundaries (where filtered values are less reliable) were always removed first.

## RESULTS

### Algorithm Sensitivity and Positive Predictive Value

To test algorithm performance on synthetic images (stable baseline), we compared the conventional and à trous detection algorithms with two versions of our approach.

#### VST 1.

With the VST1 variant, spatial filter FWHM = 8 pixels, smoothing iterations = 3, baseline iterations = 5, and *τ* = 4.5.

#### VST 2.

With the VST2 variant, spatial filter FWHM = 12 pixels, smoothing iterations = 4, baseline iterations = 8, and *τ* = 4.8.

These variations were chosen to show the effects of using relatively conservative filter sizes, such as may be used when the baseline varies considerably (VST 1), and also larger filters, which are more able to optimize detection when the baseline is approximately constant (VST 2). An example image is shown in Fig. 3, alongside the output from different detection algorithms.

Following previous work, we report the sensitivity and positive predictive value (PPV) in each case for different spark amplitudes at different (fixed) SNRs (4, 6, 15, 17, 20, 22). The sensitivity is defined as the probability that an event was detected, calculated as the ratio of the correctly identified events to the total number of sparks present in the images. The PPV gives the probability that a detection corresponds to a true spark, i.e., the number of true positive detections divided by the total number of detections. The ideal value is therefore 1 in each case.

Figure 4 shows the results, in which curves shifted to the left indicate better performance. From this, one can see that the VST 1 algorithm achieves a sensitivity slightly improved over that of the conventional algorithm, along with a more substantially improved PPV. In fact, the VST 1 leads to fewer false positives than both the conventional and à trous algorithms, but the higher sensitivity of the à trous algorithm also raises its PPV. The overall performance of VST 1 in this test may therefore be considered to fall between that of both previous methods.

The parameters used with the VST 2 variant led to clear improvements in both sensitivity and PPV. Indeed, only 9 false positives were produced by the VST 2 in these tests, compared with a total of 25,710 true detections; the low PPV value when the SNR = 2 and amplitude = 0.05 is a result of a single false detection and no true detections in this particular image set. Apart from this outlier, the PPV is then close to one using VST 2 for all other data points. While the sensitivity of the à trous algorithm is slightly higher than VST 2 for very low spark amplitudes, this coincides with a reduction in PPV, so that the greater number of true detections comes at a cost of false positives (in total, the à trous algorithm produced 25,579 true and 336 false detections). The sensitivity curve of VST 2 shows a steeper ascent, however, and reaches its maximum more quickly, which leads to a slightly higher number of detections. The steeper curve of VST 2 may be considered desirable inasmuch as it indicates a sharper cutoff between detectable and undetectable events, and so a greater consistency of detection for events with fixed properties.

### Spark Detection at Low SNRs

To assess its performance at low SNRs, we retested the VST algorithm using simulated images with a fixed spark amplitude but varying SNRs, allowing direct comparison with the published performance of the matched filter detection algorithm (see Fig. 3 in Ref. 12). A synthetic spark was generated using the IDL code, and default parameters were provided with the matched filter algorithm. Our algorithm was applied using the VST 1 and VST 2 parameters, except that the FWHM of the Gaussian spatial filters was set to 5 pixels. This corresponded more closely with the FWHM of the much narrower synthetic sparks used in this test, allowing a fairer comparison with the matched filter.

Results for sensitivity and PPV for our algorithm are plotted in Fig. 5. SEN_{50} and PPV_{50} are the SNRs at which the algorithm gives half of its maximum sensitivity and PPV, respectively; smaller values are therefore preferable, since this indicates that the algorithm performs better in noisier images. For the conventional algorithm, SEN_{50} and PPV_{50} are reported to be ∼1.2 and ∼1.35, respectively, whereas for the matched filter detection algorithm, these are ∼0.8 and ∼0.5, demonstrating significant improvements in both sensitivity and PPV (12). Using VST 1, SEN_{50} and PPV_{50} were approximately 1.2 and 0.5, indicating similar sensitivity to the conventional algorithm and similar PPV to the matched filter detection algorithm. However, for VST 2 the SEN_{50} is ∼0.7, and no meaningful value for PPV_{50} could be found because the PPV never dips below 50% at any tested SNR. This suggests that the VST algorithm is no less capable of detections at low SNRs compared with the matched filter algorithm when the baseline is constant, while offering a further substantial reduction in the likelihood of false positives.

### False Positive Detections in Experimental Images

We wished to test the robustness of the false positive estimates calculated using simulated images, in case the additional sources of noise in experimental images (see above) might degrade performance. Table 1 shows the result of applying the detection algorithms (with the same settings as above) to the Ca^{2+} concentration images with estimated SNR ≥ 2. These represent experimental linescans in which all detected events are known to be false positives. In all cases, prefiltering (using the simple median-replacing algorithm outlined above) to remove isolated bright pixels helps reduce false detections. It can be seen that the conventional algorithm produces an average of almost one false detection per 512 × 512 pixels, while the à trous algorithm gives slightly more. The VST algorithm is much more robust, with ∼0.1 false positive detections per image when using smaller filters and *τ* = 4.5, which falls to 0.03 with larger filters and higher threshold *τ* = 4.8.

### Application to Ca^{2+}-Linescan Images From Vascular Smooth Muscle

Figure 6 shows the application of spark detection to a linescan of retinal arteriolar smooth muscle. This represents a more complex type of image than those for which previous spark detection algorithms were designed, depicting multiple cells in which sparks occur both from baseline and superimposed on prolonged Ca^{2+}-oscillations (8). Sparks that have been manually identified by both of two experienced independent observers (T. M. Curtis and J. G. McGeown) are shown in Fig. 6*B*. The VST algorithm (using similar settings to VST 1, but with spatial filter FWHM = 12 pixels to better correspond to the FWHM of sparks in these images) detected all nine events that were identified manually, in addition to identifying eight further potential sparks. In the final filtered and baseline subtracted image used for detection (Fig. 6*C*), each pixel has been divided by the estimated noise standard deviation so that it may be effectively seen as a z-score. The detection of an event requires pixels in this image to exceed the threshold, so that from this one can directly ascertain the maximum threshold that would still be low enough to detect any given event. In this case, all manually identified sparks would have been found using *τ* = 7.5, while all the additional events would have been missed. Given the very low number of false positives expected when *τ* = 4.5, this suggests that the algorithm is capable of considerably more sensitive detection than is possible by eye, and the additional detections are likely to correspond to true Ca^{2+} release events. Nevertheless, higher thresholds could be used to increase confidence in the detections yet further, while still maintaining better detection that can be achieved manually.

As previously mentioned, the development of our algorithm was motivated by a need for reliable detection when the baseline varies substantially (spatially and temporally) within an image. We generated a simulated image of arteriolar smooth muscle in which 23 synthetic sparks were added to a varying baseline (Fig. 7*A*). The sparks all had a fixed absolute amplitude, but the varying baseline changed the relative amplitude, and thus the ease of detection, so that some sparks were difficult to identify even by eye. None of the previously described algorithms produce acceptable results. The changing baseline causes high correlations with the matched filter even when no sparks are present and there is no effective control of false positives, while the thresholding of the conventional and wavelet algorithms are unable to identify the sparks. By contrast, the VST algorithm was capable of detecting all the events, and yielded no false positives (Fig. 7*E*).

Finally, the results of applying the VST detection to an experimental linescan are shown in Fig. 8. Sparks are detected throughout the image in different cells, despite differences in resting fluorescence within each cell. In the selected traces, one can also see large sparks summating to form longer [Ca^{2+}] elevations, and further release events can be resolved and detected during the decaying phases. Note that some of the increased variability seen during elevated fluorescence (visible particularly in Fig. 7*C*) is a consequence of the increased noise standard deviation of Poisson noise, and spatiotemporal filtering used to reduce noise might cause these noise variations to resemble potential sparks. For this reason, a detection algorithm that takes into consideration the true noise distribution is required to ensure that only regions in which the fluorescence exceeds this elevated local noise level are detected. The same algorithm parameters were used for Figs. 6–8. Experimental images were preprocessed using the approach described above to remove extreme pixels prior to detection.

## DISCUSSION

### Comparison With Other Algorithms

In both the current and conventional algorithms the aim is to reduce noise before thresholding the resultant image to detect sparks. We have attempted to retain the principal advantages of the conventional algorithm while improving its performance and generality. This was achieved by replacing nonlinear normalization and filtering, which can be unduly affected by baseline variability and differences in spark density, with more predictable variance stabilization and linear filtering operations, based on a more robust noise estimate.

Our proposed detection method also shares some similarities with the à trous algorithm (inasmuch as both are based on separable, Gaussian-like filtering) (22) and the matched filter detection algorithm (in that filter sizes should be appropriate for the size of events to be detected) (12). However, in addition to differences in implementation, one may appreciate a more fundamental distinction in these approaches by considering the surface plots of the filters shown in Fig. 9. The matched filter naturally looks most “sparklike” (Fig. 9*A*). When the filtering is applied, this leads to high correlations with any similar sparks found within the image. However, the shape and relatively large, flat background of the filter mean that it also correlates with non-sparklike changes in fluorescence. This makes it insufficiently selective to cope with background changes. By contrast, the filters associated with the IUWT produce zero correlations not only when the signal is constant under the support of the filter but also if it is linear, and so IUWT filters of an appropriate size for spark detection are less prone to give non-zero correlations because of a more slowly varying baseline. This means that any correlation significantly higher than that expected from noise alone is much more likely to indicate a spark, rather than a baseline change. The à trous algorithm uses the 2D IUWT, which effectively convolves an image with a filter that has a positive, Gaussian-like central component, surrounded by negative values (Fig. 9*B*). The filter associated with the VST algorithm is based on the 1D IUWT (albeit extended into 2D by the separate Gaussian spatial filter). This has a similar positive central component but is surrounded by negative values only along the temporal dimension (Fig. 9*C*). In both cases, a significant correlation occurs whenever the weighted sum of all the pixels under the positive central region is sufficiently large relative to the weighted sum of all other surrounding pixels under the filter support. This means that a fluorescence increase along any dimension could potentially lead to a significant correlation using the 2D IUWT, whereas the increase must occur along the time dimension using the 1D IUWT. Consequently, if the baseline is known to be flat, then the 2D IUWT can be preferable because it makes use of the fact that a spark peak is completely surrounded by lower values, but the 1D IUWT is needed to process images such as ours to prevent cell boundaries (e.g., Fig. 6*A*) from being misclassified as [Ca^{2+}] release events. Nevertheless, we have shown that appropriate settings allow the VST algorithm to achieve similar sensitivity to the à trous algorithm when applied to constant baseline images, while retaining an improved PPV (Fig. 4).

Event detection using our algorithm then results from a rejection of the hypothesis that the (square root of the) signal is either locally constant, or increasing or decreasing linearly along the temporal dimension under the filter support. Clearly, as filter sizes increase this criterion becomes harder to meet if the baseline is not truly constant or many other sparks occur within the filter support, so that when using very large filters there is an increased risk of detecting non-spark-like fluorescence increases (such as the peaks of more prolonged waves or oscillations). Consequently, a balance must be struck between sensitivity and specificity. This means that smaller filters are to be preferred whenever they are applied to images that are very complex, while larger filters can improve sensitivity in simpler images. In any case, the need for a locally constant or linear baseline is clearly less stringent than the requirement of the conventional and matched filter algorithms for a constant baseline everywhere.

### Algorithm Parameters and Refinements

Suitable parameters to use with the VST algorithm can be inferred from knowledge of the image and can generally be kept constant for a set of similar linescans. The threshold parameter *τ* is the simplest to determine. At *τ* = 5, the likelihood of false positive detections is already extremely low, and so higher values are only really suitable if one observes other phenomena in the image (such as brief, subtle variations in fluorescence or pattern noise) that are smaller than one wishes to detect. The full-width at half maximum (FWHM) of the spatial filter should be less than or equal to the FWHM of a typical spark. The most difficult choice is then the selection of temporal filters, for which several combinations may need to be tried, although the parameter space using our proposed iterative algorithm is rather small. We found that no more than five smoothing iterations should be applied to approximate the baseline for microvascular images, because more resulted in the misclassification of the peaks of longer events as sparks. Although we have used no fewer than three iterations for noise suppression here so as to improve the detection sensitivity of very small events in very noisy images, in practice, two iterations may be preferable for images in which sparks may be seen occurring very close to one another in time. For example, in Fig. 8*D* a separate, brief rise in fluorescence can be distinguished immediately before the second detected spark, but this was not recognized as a separate event because it was merged with the much larger, subsequent spark during detection. By using only two smoothing iterations instead of three, the events would have been separated.

A benefit of the strategy we have outlined is that it results in a single image to be thresholded for detection, which can be viewed and one can determine by visual inspection whether appropriate settings have been used to detect the type of events that are of interest (see Fig. 3, *E* and *F*, and Fig. 6*C*). However, improved detection of events with more widely varying properties might be achieved by replacing the temporal filtering step with a discrete undecimated or continuous wavelet transform, which effectively allows the detection to be applied using different levels of filtering. Because a 1D transform is applied along the time dimension, one is not limited to the use of the IUWT, the benefits of which are most apparent for the detection of isotropic objects in higher dimensions. The primary cost of this would be an increase in the complexity of the algorithm both in terms of its implementation and use, because multiple images would be used for detection and the information must somehow be recombined to give the final result.

### The Control of False Positives

The results in Table 1 show that, when applied to real linescans, one can be more confident that a detected region given by the VST algorithm corresponds to a real spark when compared with a region output by the conventional and à trous algorithms, with the number of false positives reduced by ∼90% or more when using the VST approach. One explanation is that the VST algorithm is more tractable. Because it involves no nonlinear operations other than the square root transform for variance stabilization, the eventual threshold is more simply defined in terms of true noise standard deviations. The algorithm assumes that the noise is essentially Poisson in nature and that the signal is constant (or linear) under the support of the filters used when no event is observed. Under these conditions, event detection requires a pixel to have a value greater than 4.5–5 times the noise standard deviation, which is expected to be a very rare event. On the other hand, the nonlinear normalization and median filtering steps used by the other algorithms, in addition to possible wavelet denoising, mean that the eventual noise distribution is less predictable and the threshold less interpretable.

The high numbers of false positives reported here when using the conventional and à trous algorithm can be reduced in practice by subsequently filtering detections based upon measurements (22). This has not been used here primarily because the outcome of such postprocessing is therefore dependent on the measurement strategy employed [e.g., filtering and threshold crossing (6, 22), or Gaussian and exponential fitting (11)], and, in any case, such postprocessing may be applied to the output of any detection algorithm, including the VST, to improve its performance. Because measurement strategies do not generally take into consideration the noise properties, we consider the goal of the detection algorithm to be to identify regions in which an increase in signal is present, such that the null hypothesis that it is due to noise is rejected. The limitations of low SNRs mean it is likely that events that are too small to accurately measure using current methods can nonetheless be reliably detected.

### Detection Bias and Variance Stabilization

Previously, it has been suggested that one may apply spark detection in complex images using global thresholds after baseline subtraction (4, 6). The definition of a global threshold in terms of spark amplitude (ΔF/F_{0}), rather than noise, was justified on the grounds that widely different criteria would be used to identify sparks within the same image if the threshold were adapted to local noise (6). However, using an absolute threshold during detection not only implies that many events can be missed when the baseline fluorescence is low (since a high threshold is required to be suitable when the baseline is elevated), but it affords room to introduce more bias. This occurs because the probability of obtaining a false positive is dependent on the size of the threshold criterion relative to the noise, and so also increases with the baseline. Consequently, one potentially underestimates and overestimates spark frequencies for low and high baselines respectively. Using variance stabilization to effectively adapt thresholds greatly reduces this risk, so that at least the likelihood of false positives remains similar everywhere.

Furthermore, although we have used a straightforward Poisson noise model for simplicity, the method of variance stabilization employed here has been shown to be suitable when an additional source of Gaussian noise is also present in an image (25). Therefore, our algorithm is able to deal with additional read noise, which may occur when using PMTs, but which is more commonly associated with full-frame imaging using conventional charge-coupled device cameras (14). Nevertheless, if one is confident that the noise in an image is similar everywhere, the variance stabilization step can be easily omitted from our algorithm by simply not computing square roots. This may even lead to a slightly improved sensitivity in some cases (since spark profiles are modified by the transform, which may have an effect on their ease of detection). One should, however, apply caution, and possibly raise thresholds accordingly, since the very good PPV reported here is assisted by the variance stabilization, and so any apparent increase in detection numbers from omitting this step should be carefully checked for false positives.

### Conclusions

Here, we have described in some detail the issues encountered in automating spark detection, and demonstrated that improvements in detection sensitivity, reliability, and generality do not necessarily require a complex algorithm. We have introduced a detection algorithm that foregoes nonlinear normalization and filtering in favor of a strategy based on separable linear filtering and variance stabilization. By using both synthetic and experimentally acquired images, we have shown that our algorithm allows one not only to achieve similar or better sensitivity to current automated methods, but also to extend the range of images that can be processed automatically while reducing the risk of false positive detections in all cases.

## GRANTS

This work was supported by the European Social Fund (to P. Bankhead).

Present addresses: P. Bankhead, Nikon Imaging Center, Bioquant BQ 0004, Im Neuenheimer Feld 267, D-69120 Heidelberg, Germany; C. Norman, Department of Physiology, Naresuan University, Phitsanulok 65000, Thailand.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

- Copyright © 2011 the American Physiological Society