Segmentation of Cortical Spreading Depression Wavefronts Through Local Similarity Metric

TL;DR

This paper introduces a new region-based segmentation method for cortical spreading depression in microscopy images, improving accuracy and robustness against noise and inhomogeneity.

Contribution

The proposed method uses a local intensity similarity measure with distance maps, achieving higher accuracy and robustness than existing techniques.

Findings

DICE index of 0.9859, 6% higher than previous methods

79.9% reduction in root mean square error

Robust segmentation boundary in noisy, inhomogeneous images

Abstract

In this paper, we present a novel region-based segmentation method for cortical spreading depressions in 2-photon microscopy images. Fluorescent microscopy has become an important tool in neuroscience, but segmentation approaches are challenged by the opaque properties and structures of brain tissue. These challenges are made more extreme when segmenting events such as cortical spreading depressions, where low signal-to-noise ratios and intensity inhomogeneity dominate images. The method we propose uses a local intensity similarity measure that takes advantage of normalized Euclidean and geodesic distance maps of the image. This method provides a smooth segmentation boundary which is robust to the noise and inhomogeneity within cortical spreading depression images. Experimental results yielded a DICE index of 0.9859, an increase of 6% over the current state-of-the-art, and a reduction of root mean square error by 79.9%.

Figures

Equations

$$LSF(x,lc)={\int }_{y\in N_x\ne x}\frac{{\left|I(y)-lc\right|}^2}{d(x,y)}dy$$

$$l{c}_1(x)=\frac{{\int }_{\Omega }M(x,y)I(y)H_{\epsilon }(\varphi (y))dy}{{\int }_{\Omega }M(x,y)H_{\epsilon }(\varphi (y))dy}$$

$$l{c}_2(x)=\frac{{\int }_{\Omega }M(x,y)I(y)(1-H_{\epsilon }(\varphi (y)))dy}{{\int }_{\Omega }M(x,y)(1-H_{\epsilon }(\varphi (y)))dy}.$$

$$M(x,y)\{\begin{array}{cc}1& d(x,y)<r\\ 0& else\end{array}$$

$$\begin{array}{rl}& E^{RLSF}(x,\varphi (x))=\\ & ={\lambda }_1{\int }_{\Omega }LS{F}_1(x,l{c}_1(x))H_{ϵ}(\varphi (x))dx\\ & +{\lambda }_2{\int }_{\Omega }LS{F}_2(x,l{c}_2(x))(1-H_{ϵ}(\varphi (x)))dx\\ & +\mu \int {\delta }_{ϵ}(\varphi (x))|▽\varphi (x)|dx\end{array}$$

$$D_I(x)=(2H_{ϵ}({\varphi }_0(x))-1)\frac{d_0(x)g_0(x)}{ma{x}_{y\in \Omega }\left\{d_0(y)g_0(y)\right\}}$$

$$\varphi (x,T)=\{\begin{array}{cc}1& D_I(x)\le T\\ -1& else.\end{array}$$

$$\begin{array}{rl}& E^{LSM}(T)=\\ & {\lambda }_1{\int }_{\Omega }LS{F}_1(x,l{c}_1(x))H_{ϵ}(\varphi (x,T))dx\\ & +{\lambda }_2{\int }_{\Omega }LS{F}_2(x,l{c}_2(x))(1-H_{ϵ}(\varphi (x,T)))dx.\end{array}$$

$$\frac{\partial T}{\partial t}={\lambda }_1{\int }_{N_T}LSF(x,l{c}_1)dx-{\lambda }_2{\int }_{N_T}LSF(x,l{c}_2)dx,$$

$$T_{n+1}=T_n+\Delta t\Delta T_n$$

$$\Delta t=\frac{1}{Nma{x}_{x\in N_T}(LS{F}_1(x,l{c}_1)-LS{F}_2(x,l{c}_2))}$$

Full text

**Index Terms— ** Biomedical Image Analysis, Image and Video Segmentation, Microscopy, Neuroimaging

1 Introduction

In vivo imaging of the central nervous system has become an increasingly prevalent approach to neuroscience research. Advances in transgenic animal models and two-photon microscopy have allowed researchers to directly observe the development and behavior of many cells and structures within the brain which hitherto have only been approximated through in vitro studies and indirect observation [1, 2, 3]. However, due to the complexities of the tissues being imaged, this technique has several challenges which affect image quality. The primary challenge is that tissue is adept at scattering light, which results in noise and decreasing fluorescence intensity as the imaging depth increases [4]. This problem of noise and clarity is compounded by the fluorescence of tissue outside the focal plane, as well as partial or total occlusion of objects in the image by structures that lie above the focal plane, such as blood vessels. Because of these difficulties, the development of robust image segmentation strategies, both for general and specific applications, is critical to successful analysis.

Traditional approaches to segmentation of biological images such as active contours [5, 6], level set methods [7], and watershed techniques [8] fail due to poor contrast and inhomogeneity of the observed intensity. In this paper, we propose a method for segmenting the wavefront of cortical spreading depressions (CSDs) in two-photon calcium images of mouse brains. This method is inspired by the region-based, noise-robust method proposed by Nui et al. [9] and converts a level set method to a threshold-based method that is essentially based on fast marching [10]. The primary motivation of the proposed method is to conserve the shape of the segmentation boundary between iterations, preventing the final segmentation from being overly sensitive to the discontinuities in the CSD wavefront that occur due to noise and occlusion. Furthermore, the proposed method seeks not to segment an object in an image from the background in the traditional sense, but rather to separate the image functionally into two regions: one where the neurons have depolarized in response to the CSD and one where they have yet to depolarize.

2 Background

2.1 Cortical Spreading Depressions

CSDs are a pathological neurological phenomena associated with traumatic brain injury, stroke, and migraines. They are characterized by a slowly-propagating (1.7-9.2mm/min) wave of near-complete neuron depolarization, followed by a period of suppressed neural activity. In traumatic brain injuries and stroke, CSDs occur repeatedly with increasing intensity over time, damaging healthy neurons and exacerbating the injury [11, 12]. Because of this, research into the causes of CSDs, the neurological responses they induce, and methods of inhibiting their generation and propagation are important to improving patient outcome.

While CSDs and their role in brain injury have been investigated in previous studies, these primarily relied on the use of electrodes for in vivo recording of the generation and propagation of CSDs [11]. In recent years, researchers have begun to use two-photon microscopy to image CSDs [2]. These imaging studies are particularly useful in investigating the suspected role of non-neuronal cells, such as microglia, in the propagation of CSDs and the resulting cellular response [2, 13, 14]. In these studies, the propagation of the CSD is detected through calcium indicators (such as GCaMP) which fluoresce in response to increased intracellular calcium [15]. In this paradigm, the CSD appears as a region of increased fluorescence as neuronal dendrites depolarize in the focal plane, which rapidly propagates as a wavefront across the image (see Fig. 1).

Accurately segmenting the boundary of the CSD can provide valuable information about which regions have undergone depolarization in a given image, as well as measuring the direction and speed of the CSD. However, these imaging sequences are often characterized by a low signal-to-noise ratio (SNR), intensity inhomogeneity, and discontinuities in the visible wavefront as a result of light scattering and occlusion. These obstacles make delineating the wavefront difficult for many segmentation methods which rely on edge detection or global information. This is because the low SNR can make it challenging to accurately calculate edges within the image, while the intensity inhomogeneity and discontinuities make it challenging to calculate a global statistic that properly distinguishes the depolarized regions of the image from the rest of the image. This results in segmentations where the contour boundary fails to advance towards the CSD wavefront, or ones which ”sink into” the discontinuous or inhomogeneous regions of the wavefront, resulting in erroneous segmentation. Because of this, we seek to implement a region-based method which is robust to such noise and inhomogeneity, thus resulting in a segmentation that accurately converges to the wavefront boundary in a smooth and continuous manner.

2.2 Local Similarity

The inspiration behind the method proposed in this paper is the local similarity factor (LSF). First introduced by Nui et al. [9], the LSF provides a distance-weighted, region-based measurement of the similarity between the intensities of pixels within a region and the region’s mean intensity. For a given pixel, $x$, within an image, $I$, the LSF is defined as:

[TABLE]

where $N_{x}$ is a square-shaped window defining the local region, $d(x,y)$ is the Euclidean distance between pixels $x$ and $y$, and $lc$ is the local mean intensity. The LSF metric has two major advantages that are valuable for the problem addressed in this paper. First, the method does not require preprocessing prior to segmentation, in the form of noise reduction or contrast enhancement, for the metric to be effective. Second, the LSF-based model the authors present is robust to high levels of noise as well as intensity inhomogeneity, making LSF an appealing metric to apply to the CSD images this paper focuses on [9].

Two different LSF values are computed per pixel for each iteration of the algorithm. $LSF_{1}$ compares pixels that are inside the current segmentation boundary while $LSF_{2}$ compares pixels that are outside the current segmentation boundary. As a result, two local mean intensities, $lc_{1}$ and $lc_{2}$, are calculated when evaluated over a level set $\phi(\cdot)$:

[TABLE]

and

[TABLE]

$M(x,y)$ is a mask of the local region defined as

[TABLE]

and $H_{\varepsilon}(\cdot)$ is the regularized Heaviside function.

By combining equations (1-4), Nui et al. produce the following energy functional, called the Region-based model via Local Similarity Factor (RLSF):

[TABLE]

where $\delta_{\epsilon}(\cdot)$ is the regularized Dirac delta function and $\lambda_{1}$, $\lambda_{2}$, and $\mu$ are weighting terms. The last integral of the energy functional serves as a smoothing parameter, with larger values of $\mu$ resulting in a smoother contour.

3 Method: Local Similarity Metric

3.1 Distance map-based segmentation with LSM

In order to compensate for the high level of noise, intensity inhomogeneity, and wavefront occlusion present in confocal images of CSDs, we propose an modified version of the RLSF model which recontextualizes it from an active contour to a threshold-based approach, which we call the Local Similarity Metric (LSM).

In the LSM model, the boundary of a given threshold is defined according to a distance map representation of the image rather than the pixel intensities of the image. This normalized distance map, $D_{I}(x)$, is defined as:

[TABLE]

where $d_{0}$ and $g_{0}$ are normalized Euclidean and geodesic distance maps of the image, respectively, and $\phi_{0}(\cdot)$ is the user-defined initial contour. The initial contour’s boundary pixels, excluding those on the edges of the image, serve as the ”zero-points” of both distance maps, meaning that the value of a pixel $x$ is the shortest distance from it to one of the boundary pixels. In order to encourage a smooth curve shape over various threshold values, a median filter is applied to the map.

Using (6), we can define a pseudo-level set of the image, with respect to a given threshold $T$, as such:

[TABLE]

Using this psuedo-level set, we define the LSM energy functional as follows:

[TABLE]

Note that the smoothing term has been dropped from the original RLSF model (5), as the contour shape is strictly defined by $D_{I}$ and the choice of $T$. A smooth contour is enforced through the use of Euclidean distance in generating $D_{I}$ and the median filter applied to it. Furthermore, both $lc_{1}$ and $lc_{2}$, in (2-3), also use $\phi(x,T)$ instead of $\phi(x)$. For the sake of brevity, we will not explicitly redefine them here.

3.2 Threshold optimization through gradient descent

Once the initial contour has been defined and the appropriate distance map has been generated, the optimal threshold value for the image is calculated through the use of a gradient descent algorithm with a fast marching method-like implementation.

The partial differential equation used for the gradient descent algorithm is as follows:

[TABLE]

where $N_{T}$ is a subset of the image containing pixels that are within a specified distance from the current threshold boundary. However, the curved geometry of the CSD wavefront often leads to an imbalance within $N_{T}$ between the number of pixels outside versus inside the current threshold, which biases one side over the other. We equalize this imbalance by removing pixels on the overrepresented side from $N_{T}$.

In order to minimize (8), the following gradient descent formulation is utilized:

[TABLE]

where $\Delta T_{n}$ is a numerical approximation of (9). $\Delta t$ is the step-size, defined as:

[TABLE]

where $N$ is the number of pixels within $N_{T}$.

4 Experimental Results

4.1 Dataset and parameter selection

In order to evaluate the performance of our proposed method, we compiled a dataset of twenty 512x512 pixel images of CSDs taken from two-photon confocal calcium imaging sequences of mouse brains. CSDs were induced in the mice either through surgical occlusion of the middle cerebral artery (MCA) or a 1 M injection of KCl directly into the cortex. Images from these sequences were selected such that the dataset contained a variety of possible noise and contrast levels, as well as various shapes that the CSD wavefront can assume. While we would have liked to have had a larger dataset for this analysis, the novelty of this imaging technique and the overall cost of such experiments limits the quantity and quality of sequences at our disposal.

As stated before, our algorithm requires several parameters to be specified by the user prior to segmentation. Given the resolution of our imaging instrument, for the local window, $N_{x}$, we selected a size of 17x17 pixels, while the regional mask, $M$, had a radius of 13 pixels. For the region used to calculate the gradient descent, $N_{T}$, we considered pixels within 7 pixels of the current boundary. Finally, the weighting terms were set as $\lambda_{1}=\lambda_{2}=1$. The geodesic distance map used in creating $D_{I}$ was generated using the MATLAB (MathWorks, CA) command immsegfmm.

4.2 Performance Evaluation

To evaluate the performance of the LSM method, we compared the results of our approach with that of Chan-Vese [16] and Mukherjee et al. [17]. The implementation of Mukherjee et al. method that we used in this evaluation utilizes a fast-marching method approach which prevents contour edges that lie on the edges of an image from evolving. Chan-Vese segementation was performed with MATLAB’s built in implementation, activecontour. Both Mukherjee et al. and Chan-Vese were executed for 1000 iterations, while LSM ran for only 50 iterations.

Performance was measured using two different metrics. The first was the DICE index, which measures the overlap between the segmentation result and the ground truth [18]. This is calculated as $DICE(R_{1},R_{2})=\frac{2Area(R_{1}\cap R_{2})}{Area(R_{1})+Area(R_{2})}$ where $R_{1}$ is the ground truth and $R_{2}$ is the segmentation. A DICE index closer to 1 indicates superior performance. However, due to the large area of the image being segmented, a poor segmentation can still result in a large DICE index. Therefore, we used the root mean square error (RMSE) as a second performance metric. Here, the error being measured is the Euclidean distance between a given point on the segmentation boundary and the nearest point on the ground truth boundary. For the initialization, an initial contour was hand-drawn such that the only edge of the contour which did not lie on the image boundary was a rough approximation of the CSD wavefront, displaced approximately 40 pixels away from the wavefront boundary in the image. These initial segments had a mean DICE index of .9213 and a mean RMSE of 23.00 pixels.

Over the entirety of the dataset, our proposed method demonstrated superior performance compared to the other methods examined (Fig. 2). The LSM yeilded a mean DICE index of .9859, compared to that of .9308 and .8136 for Mukherjee et al. and Chan-Vese, respectively. Likewise, LSM had a mean RMSE of 4.52 pixels, compared to 22.57 pixels and 30.19 pixels for Mukherjee et al. and Chan-Vese, respectively. Failures of Mukherjee et al.’s method were largely characterized by the contour’s inability to evolve towards the wavefront. By contrast, intensity inhomogeneity was the primary cause of failure for Chan-Vese. The active contour would converge around the brightest portions of the depolarized region, often excluding regions where noise and occlusion obscured the increased fluorescence.

5 Conclusion

In this paper, a novel approach to segmenting wavefront boundaries in two-photon calcium images of CSDs is presented. Qualitative analysis of the performance for LSM indicates improvements in comparison to the state of the art. In addition to accurate segmentation results, the LSM method has a fast rate of convergence, requiring less than 50 iterations. The most significant downside of this method is its dependence on the user’s a priori knowledge of the wavefront’s shape, which limits throughput on larger datasets. Inclusion of the LSM into models which can incorporate a priori shape information in an adaptable manner is an appealing solution for future work.