Automatic neuron segmentation and neural network analysis method for phase contrast microscopy images

Jincheng Pang; Nurdan Özkucur; Michael Ren; David L. Kaplan; Michael Levin; Eric L. Miller

doi:10.1364/BOE.6.004395

1. Introduction

Phase Contrast Microscopy (PCM) is an important technique for the study of living cells. As illustrated in Fig. 1(a), PCM adds an annulus and a phase plate to a conventional light microscope to create a phase shift between light propagating along the “diffracted” and “surround” paths which is converted into intensity variations on the imaging plane [1]. There are two major advantages of PCM due to the fact that specimens are not to be stained as required by fluorescence microscopy. First, researchers can perform long-term observations of unstained living cells without worrying about photobleaching. Second, disruption of the observed cells and thus measurement artifacts caused by introduction of fluorescence can be avoided.

Fig. 1 (a) The optical path of phase contrast microscopy; (b) A cropped region of a PCM neuron image. The blue arrows point to somas; the red arrows point to dendrites; the yellow arrows point to places where dendrites and somas are separated by halos; the purple dash arrows point to shade-offs.

Download Full Size | PDF

However, as seen in the work of [2, 3] and as illustrated in Fig. 1(b), PCM images can suffer from imaging artifacts such as the halo and shade-off due to the imaging physics. Therefore it is challenging to build models for cell intensity distribution and design algorithms to analyze the images. Recently, various algorithms have been developed to help to automatically analyze PCM images for osteosarcoma cells [4], bovine aortic endothelial cells [5], bone marrow stromal cells [6] and fish keratocytes [7]. The above algorithms mainly employ the gradient information among halo, cell and background for segmentation.

Besides the PCM imaging artifacts, the development of analysis algorithms for neurons is even more challenging due to the requirement of identifying both cell bodies (somas) as well as dendrites, two classes of structures with very different geometries and contrasts within the data. Indeed, as illustrated in [8–12], dendrite tracing itself is already a quite challenging problem since these structures can be dim and broken in some cases. Moreover, dendrites and soma of the same neuron are sometimes separated by the halo artifacts in PCM images (see Fig. 1(b)). These gaps need to be filled for the purpose of biological connectivity analysis. Biological connectivity consists of chemical (synapses), electrical (gap junctions), physiological (voltage-sensitive ion channels, synaptic receptors) and anatomical (synaptic clefts) components. In order to study biological connectivity, structural connectivity (identifying somas and connecting corresponding dendrites) needs to first be determined [13,14]. However, computationally bridging this gap is not at all trivial. Currently, the majority of image related research of neurons has been performed on fluorescence microscopy data [15] with very limited work focusing on the development of automatic algorithm for neuron segmentation (including both somas and dendrites) and connectivity analysis for PCM images [16, 17]. The work in [16] mainly considered on dendrite tracing and quantification of the dendrite change while [17] only examined soma segmentation. To the best of our knowledge, the work in this paper is the first to present a method that segments both somas and dendrites simultaneously for PCM neuron images.

Recently, the authors of [2, 3] modeled the imaging process of PCM as a convolution operation. Based on the convolution model, [2] built an optimization function including smoothness and sparseness regularization to restore artifact-free images. Then simple segmentation methods such as thresholding [18] can be used to these restored images.

In this paper, rather than employing a pixel-based parametrization of the problem as in [2], we present a parametric image model consisting of two level set functions to represent the neuron images. Then we formulate an optimization problem merging image restoration and segmentation. Our novel approach has multiple advantages over the methods in [2, 3]. First of all, even after restoration, image segmentation must still be performed on the restored images for the purpose of image analysis in [2, 3]. Our method performs the restoration and segmentation simultaneously, so no other major segmentation operation is needed. Secondly, due to the parametric representation, both dendrites and somas can be segmented simultaneously.

Based on the above segmentation approach, we present the pipeline illustrated in Fig. 2 which performs image analysis for the PCM neuron images. This pipeline is automatic and, aside from the specification of a few parameters, requires no human interaction thereby making it suitable for large scale automated image analysis. As in [2, 3], background bias correction is performed before the essential part of the pipeline as well as the major contributions of this paper which is “variational segmentation.” In this essential part, we build an energy functional of level set functions based on the PCM physical model and image information. Then we solve an optimization problem to obtain a segmenation/restoration of the PCM images. Based on the segmentation of somas and dendrites, some simple morphological operations are performed to refine the results. The feasibility and advantages of our variational segmentation method are demonstrated by synthetic images in Section 5.1. In addition, the whole pipeline has been validated in real PCM neuron images for the neural development study [19]. The remainder of this paper is organized as follows. In Section 2, some background of the level set method is introduced. We propose and illustrate our variational segmentation method in details in Section 3. Then we introduce a set of morphological operations which aim to refine the variational segmentation results in Section 4. We present experimental results and discussions in Section 5. Finally, we draw our conclusion and discuss our future work in Section 6.

Fig. 2 Illustration of the pipeline for PCM neuron image analysis

Download Full Size | PDF

2. The level set method for image segmentation

The level set method is a very powerful and convenient approach to represent boundaries or surfaces for objects [20, 21]. The essential idea of level set methods for segmentation is to represent a curve $\vec{C}$ as the intersection of a horizontal plane z(x) = 0 of an auxiliary level set function ϕ (x); i.e., $\vec{C} = {x | ϕ (x) = 0)}$ .

Image segmentation aims to identify objects of interest and/or their boundaries. From the variational perspective, segmentation is often posed as an optimization problem by constructing a scalar energy functional $E (\vec{C})$ based on $\vec{C}$ , a curve representing the boundary of the object. The energy functional $E (\vec{C})$ , consisting of a data fidelity term and a regularization term, is always constructed in a way that when $\vec{C}$ approaches to the true boundary of the object, the minimum value of $E (\vec{C})$ is achieved.

For example, the well-known Chan-Vese model [22] implicitly assumes that object and background are Gaussian distributed with the same variance, and subsequently seeks a segmentation which minimizes partially a mean square error-type of energy functional. In particular, the Chan-Vese model is defined

E_{c v} (\vec{C}) = \int_{Ω_{1}} {[I (x) - c_{1}]}^{2} d x + \int_{Ω_{1}^{c}} {[I (x) - c_{2}]}^{2} d x + α | \vec{C} |

where I(x) is the image intensity value at location x, Ω₁ is the region inside of the curve

\vec{C}

,

Ω_{1}^{c}

is the region outside of the curve

\vec{C}

and

Ω = Ω_{1} \cup Ω_{1}^{c}

which is the whole image domain. In addition, c₁ and c₂ are the mean image intensity values for the regions Ω₁ and

Ω_{1}^{c}

respectively. Moreover, α is a coefficient,

| \vec{C} |

is the length of the curve

\vec{C}

and the term

α | \vec{C} |

basically is a smoothness regularization term.

By introducing the level set function ϕ such that $\vec{C} = {x | ϕ (x) = 0}$ , eqn. (1) can be converted into a function based on ϕ as

E_{c v} (ϕ) = \int_{Ω} {[I (x) - c_{1}]}^{2} H_{ε} (ϕ (x)) d x + \int_{Ω} {[I (x) - c_{2}]}^{2} [1 - H_{ε} (ϕ (x))] d x + α \int_{Ω} | \nabla H_{ε} (ϕ (x)) | d x

where H_ε (·) is an approximate Heaviside function [22], defined as

H_{ε} (z) = \frac{1}{2} (1 + \frac{2}{π} \arctan (\frac{z}{ε})),

and the Heaviside function is defined as

H (z) = {\begin{matrix} 1, & z \geq 0 \\ 0, & e l s e . \end{matrix}

By using calculus of variations [23], the optimization of E_cv can be realized by iteratively updating c₁, c₂ and ϕ using the following Euler-Lagrange equation

\frac{\partial ϕ (x)}{\partial t} = δ_{ε} (ϕ (x)) [α \nabla (\frac{\nabla ϕ (x)}{| \nabla ϕ (x) |}) - (I (x) - c_{1}) + {(I (x) - c_{2})}^{2}]

where

c_{1} = \frac{\int_{Ω} I (x) H_{ε} (ϕ (x)) d x}{\int_{Ω} H_{ε} (ϕ (x)) d x}, c_{2} = \frac{\int_{Ω} I (x) [1 - H_{ε} (ϕ (x))] d x}{\int_{Ω} [1 - H_{ε} (ϕ (x))] d x} .

Here we use δ_ε (·) to represent the derivative of H_ε (·) through the paper. Thus the level set method can be treated as an optimization technique which “evolve” the level set function ϕ as a function of time such that, as t → ∞, $\vec{C} (t)$ approaches the contour of interest in the image.

3. A variational formulation for PCM neuron image segmentation

In principle, phase contrast microscopy converts the phase difference of the surround wave and the diffracted wave which passes through specimen into intensity difference that can be directly measured [1] (see Fig. 1). Yin et al. [2,3] modeled the imaging process of phase contrast microscopy as a convolution operation

I (x) \propto G (x) + [a i r y (x) - δ (x)] * θ (x)

where * represents convolution; I(x) is the observed image; θ (x) is the phase retardation at the location x; δ (x) is the Dirac delta function; G(x) is a bias field caused by uneven illumination and can be estimated using flat-field correction method [24]. The function airy(x) is an obscured airy pattern defined as follow [2, 25, 26]

a i r y (x) = \frac{J_{1} (k_{m} | x |)}{k_{m} | x |} - r_{m} \frac{J_{1} (k_{m} r_{m} | x |)}{k_{m} | x |}, x = (x, y), x, y \in - N_{m}, \dots, - 1, 0, 1, \dots N_{m},

where J₁(·) is the first order Bessel function of the first kind, N_m is the radius of the airy pattern kernel, r_m = 0.85 in our case which is the ratio of the inner radius (R_in) to the outer radius of the phase ring (R_out) in Fig. 1, k_m is a scale parameter related with the properties of microscopy and the wavelength of the imaging light [25, 26]. In addition, the airy pattern is normalized such that 1 * airy(x) = 1 as in [2, 3].

In this section, we present a parametric model to represent θ (x) based on the level set framework. As previously discussed, for our application, a neuron usually consists of a soma and multiple dendrites which need to be segmented [15, 27]. Somas and dendrites differ in two critical ways. First, the shape of a soma is like a blob while the dendrites are typical tubular structures. Second, somas are darker than dendrites in PCM images and halos are more obvious than those of dendrites. This is because somas are usually thicker than dendrites. Thus, in the soma regions, the phase retardation is larger and the PCM image is darker accordingly as shown in Fig. 1(b) and illustrated in Section 5.1. In Section 5.1, a single level set function is demonstrated to have difficulty capturing both somas and dendrites simultaneously. Thus we use two level set functions ϕ₁ and ϕ₂ for somas and dendrites respectively. Specifically, the image to be recovered can be represented as the following Double Level Sets model (DLS)

θ (x) = s_{1} H_{ε} (ϕ_{1} (x)) + s_{2} H_{ε} (ϕ_{2} (x)),

where ϕ₁ and ϕ₂ are the level set functions, s₁ and s₂ are real scalars representing in some sense the average phase retardation values for somas and dendrites respectively, and H_ε (·) is the approximate Heaviside function defined in (3).

Based on the parametric model in (7), we propose the following energy function

E (ϕ_{1}, ϕ_{2}; s_{1}, s_{2}) = E_{p h y} (ϕ_{1}, ϕ_{2}; s_{1}, s_{2}) + λ_{1} E_{l o c} (ϕ_{1}, r_{1}) + λ_{2} E_{l o c} (ϕ_{2}, r_{2}) + λ_{3} E_{w t u b} (ϕ_{2}) .

The energy term E_phy is a data fidelity term based on the physical model as specified in Section 3.1. In addition, E_loc is a local data fidelity term based on localized active contour models [28] as introduced in Section 3.2. Moreover, E_wtub is a spatial weighted version of the geometric regularization term for tubular structure [29] as introduced in Section 3.3. Lastly, λ_i,i = 1 ⋯,3 are real scalars to balance these energy terms.

Here we want to point that, in the case of multiple level set functions, constraints are usually imposed to guarantee each level set function does not overlap with each other [30] etc. In our case for the segmentation of somas and dendrites, the overlap regularization does not apply. A neuron can be divided into a soma and dendrites but no clear boundaries between the soma and dendrites exist. Thus imposing overlap penalty always makes the soma and dendrites in the segmentation separated from each other which complicates the connectivity analysis. Therefore, we do not use regularization to guarantee the non-overlap for somas and dendrites for now. We will develop a regularization method that ensures smooth connection between these two structures in the future.

3.1. A data fidelity term based on the physical model

By compensating for G(x) and the multiplicative factor during preprocessing, and using the sifting property of the Dirac δ function, the imaging model (5) can be written as

I (x) = a i r y (x) * θ (x) - θ (x) .

In order to formulate an optimization problem, the physical model in (5) or (9) are always represented using a matrix form. Following [31], here we define an operator vec: ℂⁿ^×^m → ℂ^nm^×1 such that for a given array w ∈ ℂⁿ^×^m,

vec (w) = {[w_{1, 1}, \dots, w_{n, 1}, w_{1, 2}, \dots, w_{n, 2}, \dots, w_{1, m}, \dots, w_{n, m}]}^{T} .

The operator array defines the inverse of the vec operator such that

a r r a y (vec (w)) = w, vec (a r r a y (w)) = w,

where w ∈ ℂⁿ^×^m and w ∈ ℂ^nm^×1.

Then we can define ϕ₁, ϕ₂ and θ are the vectorized version of ϕ₁, ϕ₂ and θ respectively. Thus the vectorized version of the DLS model (7) is represented by

θ = s_{1} H (ϕ_{1}) + s_{2} H (ϕ_{2})

where ϕ₁ = vec(ϕ₁), ϕ₂ = vec(ϕ₂) and θ = vec(θ). In addition, (9) can be discretized and represented by the following matrix form

I = H θ,

where I is vectorized from I(x) such that I = vec(I). Moreover, H = (K − I_n) where I_n is an identity matrix, K is a matrix with special structure built from the convolution kernel airy(x) [31]. In this paper, we assume θ (x) is periodic so that K is a block-circulant-circulant-block(BCCB) matrix. As we will discuss in Section 3.5, the BCCB assumption enables efficient implementation of the matrix and vector multiplication in (11) through the Fast Fourier Transformation (FFT) [31, 32].

From the physical model (11), a data fitting term E_phy is defined as

E_{p h y} (ϕ_{1}, ϕ_{2}; s_{1}, s_{2}) = \frac{1}{2} {‖ I - H θ ‖}^{2} .

3.2. Localized active contours

In order to deal with the inhomogeneity of images, level set methods driven by local image features were proposed [28, 33]. The general idea of the localized level set method is illustrated by Fig. 3. For each point in the evolving curve $\vec{C}$ such as x, the driven force derived from image at x is only based on its neighborhood B_r(x) rather than the whole image domain such as the Chan-Vese model [22]. Therefore, localized active contours are robust to image heterogeneities and also can capture information in small details such as vessels and dendrites [28, 33].

Fig. 3 An Illustration for the localized level set method. The point x is in the evolving curve $\vec{C}$ and B_r(x) is a neighborhood of x. Note: the neighborhood B_r(x) is divided into two parts by $\vec{C}$ .

Download Full Size | PDF

Specifically, the local image fitting energy function formulation of [28] is defined as follows

E_{l o c} (ϕ; r) = \int δ_{ε} (ϕ (x)) \int B_{r} (x, x^{'}) \cdot F_{l o c} (I (x^{'}), ϕ (x^{'})) d x^{'} d x

where F_loc denotes a metric to evaluate local adherence to a given model. Mathematically, the neighborhood of x with radius equals r is defined by

B_{r} (x, x^{'}) = {\begin{matrix} 1, & {‖ x - x^{'} ‖}_{2} \leq r \\ 0, & e l s e . \end{matrix}

where ∥·∥₂ indicates the l² norm. In this paper, we use a localized version of the mean separation energy [34] for F_loc

F_{l o c} (I (x), ϕ (x)) = {(u_{x} - v x)}^{2},

where u_x and v_x are the mean values of the inside and outside parts of B_r(x) respectively and they can be represented by

u_{x} = \frac{\int B_{r} (x, x^{'}) \cdot H_{ε} (ϕ (x^{'})) \cdot I (x^{'}) d x^{'}}{\int B_{r} (x, x^{'}) \cdot H_{ε} (ϕ (x^{'})) d x^{'}}

v_{x} = \frac{\int B_{r} (x, x^{'}) \cdot [1 - H_{ε} (ϕ (x^{'}))] \cdot I (x^{'}) d x^{'}}{\int B_{r} (x, x^{'}) \cdot [1 - H_{ε} (ϕ (x^{'}))] d x^{'}} .

3.3. A weighted geometric regularization of tubular structure

The most well known geometric regularization is the smoothness constraint in [22, 35] as also defined as the third term in (1) which shrinks the curve in the normal direction. For the PCM neuron image analysis in this paper, we take advantages of a geometric regularization term E_tub for tubular structure proposed by [29]. This tubular regularization is defined on the base of local geometric properties of the level set function. It is similar as the localized level set method introduced in Section 3.2, but the major difference is that only the level set function information rather than image information is used to construct the regularization energy term.

An intuitive illustration of the effect of E_tub can be found in Fig. 4. Rather than shrink the curve into a point as smoothness regularization [22, 35] does, the tubular regularization would extend curve into the major principal direction while keep the other direction unchanged.

Fig. 4 Curve evolution under the geometric tubular structure regularization. From left to right: iteration 0, 20, 40 and 60 respectively.

Download Full Size | PDF

By incorporating this geometric regularization term, most of the gaps between broken dendrites, dendrites and their corresponding somas in the PCM images could be filled for neural connectivity analysis. However, as illustrated later in this section, the energy term E_tub in [29] might impose the regularization in unnecessary regions and thus cause false dendrites. In order to to reduce the false dendrite artifacts, we introduce E_wtub, a spatial weighted version of E_tub as follows.

E_{w t u b} (ϕ) = \int w (x) δ_{ε} (ϕ (x)) f (T r M^{- 1} (x)) | \nabla ϕ (x) | d x

where

M (x) = \int H_{ε} (ϕ (x^{'})) B_{r} (x, x^{'}) \nabla ϕ (x^{'}) \nabla^{t} ϕ (x^{'}) d x^{'},

Tr(·) represents the trace of a matrix, f is a non-negative decreasing function from ℝ⁺ to ℝ⁺. The only difference between E_wtub and E_tub is a spatial varying weight term w(x). Here we define w(x) is static and independent of ϕ as

w (x) = \exp (- d i s t (B W (x)) / a) + b

where a and b are two real positive numbers. BW(x) is a binary image. Here we use the soma initializations as BW(x) which is illustrated in Section 3.4. The function dist(·) calculates the Euclidean distance transform of a binary image. For a zero pixel in BW(x), dist(·) assigns a number that is the distance between that pixel and its nearest nonzero pixel of BW(x). For a nonzero pixel, dist(·) assigns zero to that pixel. Therefore w(x) imposes stronger constraints near the somas and thus reduces false dendrites in the area away from somas. An example of the weight map w(x) is Fig. 5 by setting a = 20 and b = 0.5.

Fig. 5 Left: A binary image; Right: The corresponding weight map w(x) with a = 20 and b = 0.5.

Download Full Size | PDF

The motivation of introducing w(x) is that term E_tub in [29] is independent of image information which does not work well in our case of neuron images. Strong tubular regularization could encourage connections of dendrites to corresponding somas separated by halos. However, if the same weight for the regularization applies over the whole image, false dendrites tend to be generated in areas away from somas as illustrated in Fig. 6. In order to reduce the false created dendrites, w(x) aims to put more tubular regularization near somas and less for regions far from somas. By incorporating the term E_wtub, most of the gaps between broken dendrites, dendrites and their corresponding somas in the PCM images can be filled for neural connectivity analysis.

Fig. 6 Effects of the tubular regularization for segmentation.

Download Full Size | PDF

The effect of the weighting term w(x) as well as the tubular geometric regularization can be found in Fig. 6. In addition, the parameters a and b are set as a = 20 and b = 0.5 for the real PCM images in Section 5.2. In order to show that the E_wtub is not that sensitive to a and b, we consider four combinations of (a,b) values that are both above and below the (20,0.5) values and include the corresponding segmentation results in Fig. 6. From Fig. 6, we find that, without the term E_wtub, some dendrites fail to connect to their corresponding somas as pointed by red arrows in Fig. 6(c). Using a constant weighting term, dendrites are connected to somas with the result that some false dendrites are produced as denoted by light blue arrows in Fig. 6(d). By comparing with Fig. 6(b) to (e–h), we see that varying the a and b values has little impact on the final segmentation. Indeed, there is only a single readily identified difference as denoted by the white arrow in Fig. 6(f)–(h).

3.4. Curve initialization

Curve initialization is very important for level set methods. Good initializations can yield more accurate segmentation results and also make algorithms converge rapidly. Instead of using manual initialization, we propose an automatic initialization method as described in Fig. 7.

Fig. 7 Illustration of the initializations for somas and dendrites. The input image to be processed is in Fig. 2 for curve initializations. (a): local standard deviation (b) Otsu’s thresholding results on (a); (c): closing and hole filling on (b); (d): erosion of (c); (e): steerable filtering enhancement for dendrite initializations; (f) The minimum error thresholding results on (e); (g) The final dendrites initialization.

Download Full Size | PDF

To begin with, two feature maps (Fig. 7(a) and Fig. 7(e)) for somas and dendrites are obtained from the original input image (Fig. 2). For somas, the feature map is the local standard deviation of the image [36]. For dendrites, a steerable filter using the second order Gaussian derivative template [37] is applied to enhance tubular structures. The radius of the local standard deviation in is set to 30 pixels while the standard deviation of the Gaussian template is three pixels in this paper. These two parameters are set based on prior information concerning the approximate sizes of the somas and dendrites respectively. Then by using Otsu’s method [18] and the minimum error thresholding [38], preliminary initializations for somas and dendrites (Fig. 7(b) and (f)) are obtained respectively. The morphological closing and hole filling operations are then performed to make the initialization cover most of the somas as in Fig. 7(c). The structuring element of the closing is a disk with radius equals seven pixels. Next, a morphological erosion using a disk with radius equals ten pixels as the structuring element is used to get the final soma initializations Fig. 7(d). There are two benefits for the erosion operation. Firstly, the initialized soma size is reduced which makes the algorithm converge fast. Secondly, it can keep more dendrites initializations since Fig. 7(g) is obtained by excluding the final soma initializations from the preliminary dendrite initialization.

3.5. Implementation

By combining the gradient flows of terms E_phy, E_loc and E_wtub, the optimization of the energy functional (8) can be achieved using the gradient descent flow method by following Algorithm 1 iteratively.

As summarized by Algorithm 1, there are 4 steps during each iteration. The gradient flows of E_phy, E_loc and E_wtub are defined in (17), (20) and (21) respectively below.

Using the calculus of variation, the gradient flow of E_phy(ϕ₁,ϕ₂,s₁,s₂) with respective to ϕ_k and s_k(k = 1,2) can be obtained as

Algorithm 1. Optimization for energy functional (8) for each iteration

View Table

and s_k(k = 1,2) can be obtained as

\frac{\partial E_{p h y} (ϕ_{1}, ϕ_{2}; s_{1}, s_{2})}{\partial ϕ_{k}} = - s_{k} δ_{ε}^{T} (ϕ_{k}) H^{T} (I - H θ),

\frac{\partial E_{p h y} (ϕ_{1}, ϕ_{2}; s_{1}, s_{2})}{\partial s_{k}} = - H_{ε}^{T} (ϕ_{k}) H^{T} (I - H θ) .

Using the array(·) operator defined in Section 3.1, the array version of (17) is

\frac{\partial E_{p h y} (ϕ_{1}, ϕ_{2}; s_{1}, s_{2})}{\partial ϕ_{k}} = a r r a y (- s_{k} δ_{ε}^{T} (ϕ_{k}) H^{T} (I - H θ))

According to [28], the gradient descent flow of E_loc defined in (12) is

\frac{\partial E_{l o c} (ϕ_{1}, r_{1})}{\partial ϕ_{1}} = - δ_{ε} (ϕ_{1} (x)) \int B_{r_{1}} (x, x^{'}) δ ε (ϕ_{1} (x^{'})) \cdot (\frac{{(I (x^{'}) - u_{x})}^{2}}{A_{u}} - \frac{{(I (x^{'}) - v_{x})}^{2}}{A_{v}}) d x^{'}

where

\begin{matrix} A_{u} = \int B_{r} (x, x^{'}) \cdot H_{ε} (ϕ_{1} (x^{'})) d x^{'} \\ A_{v} = \int B_{r} (x, x^{'}) \cdot [1 - H_{ε} (ϕ_{1} (x^{'}))] d x^{'} . \end{matrix}

For the weighted tubular regularization term E_wtub, according to [29], the gradient flow of (15) is

\frac{\partial E_{w t u b} (ϕ_{2})}{\partial ϕ_{2}} = - w (x) δ_{ε} (ϕ_{2} (x)) {f (T r M^{- 1}) \nabla \cdot (\frac{\nabla ϕ_{2}}{| \nabla ϕ_{2} |}) + \nabla^{t} f (T r M^{- 1}) \cdot \frac{\nabla ϕ_{2}}{| \nabla ϕ_{2} |} + \nabla^{t} ϕ_{2} L \nabla ϕ_{2}}

where

L (x) = \int B (x, x^{'}) δ_{ε} (ϕ_{2} (x^{'})) | \nabla ϕ_{2} (x^{'}) | \dot{f} (x^{'}) M^{- 2} (x^{'}) d (x^{'}),

and

\dot{f}

denotes the derivative of f.

The difficult part of implementing (17) to (21) is the matrix and vector multiplication Hθ. More specifically, as the size of matrix H grows rapidly as the size of image increases making direct, space domain implementation of the model computationally cumbersome. If we assume θ(x), the image to be restored, is n × n, then the size (also the memory cost) of H is n² × n². The time cost of the matrix and vector multiplication operation is O((n²)²) = O(n⁴).

By assuming a BCCB structure, the matrix and and vector multiplication Hθ can be implemented very efficiently through FFT as [31]

H θ = vec (ℱ^{- 1} (ℱ (a i r y (x) - δ (x)) \circ ℱ (θ (x))))

where ℱ(·) and ℱ⁻¹(·) are discrete Fourier transform and its inverse form respectively, ○ is the component-wise multiplication. In addition, we have

ℱ (a i r y (x) - δ (x)) = ℱ (a i r y (x)) - 1.

Thus the implementation of the Dirac δ function can be done in the spatial frequency domain. In addition, from (22), we can see that the matrix and vector multiplication is mainly performed by using several FFT and IFFT operations. Thus the time cost equals O(n² log(n)). Moreover, there is no need to build the huge matrix H and thus the memory cost for the frequency domain implementation is O(n²). Alternative choices for modeling K and the resulting computational methods can be found in [31, 39, 40].

For the two real scalars s₁ and s₂, based on (18), s₁ is updated by numerically solving the following minimization problem (assuming ϕ₁, ϕ₂, s₂ fixed)

\underset{s_{1}}{argmin} | H_{ε}^{T} (ϕ_{1}) H^{T} (I - H θ) |, s u b j e c t t o, 0 < s_{1} < u b .

Here we set ub = 0.2 which is a upper bound of s₁. Next, update s₂ by numerically solving the following minimization problem (assuming ϕ₁, ϕ₂, s₁ fixed)

\underset{s_{2}}{argmin} | H_{ε}^{T} (ϕ_{2}) H^{T} (I - H θ) |, s u b j e c t t o, 0 < s_{2} \leq s_{1} .

By using the constraint 0 < s₂ ≤ s₁, the level set function s₁ would catch the brighter objects as shown in Section 5.1. We update s₁ and s₂ once for every ten updates to the level set functions. This change reduces computation load and makes little difference to the final results.

4. Morphological refinement

An example of the results of our approach to identify somas and dendrites is provided in Fig. 8(b) and (e) respectively. As is typical for problems of this type [27, 41], a degree of post processing can be helpful to refine the accuracy of the basic segmentation methods. Toward this end, the final soma results are obtained by a morphological dilation as seen in Fig. 8(b). The structuring element of the dilation is a disk with radius equals two pixels. To get rid of tubular structures which are not connected to and far away from somas, morphological reconstruction is used. Morphological reconstruction is a useful method for selecting meaningful components from a binary image [42, 43]. Morphological reconstruction eliminates the mask components that do not overlap with the marker. We assume that dendrites always connect to somas. Thus in the neuron images, the mask is the segmented dendrites as in Fig. 8(e) and the marker is the segmented somas as in Fig. 8(c).

Fig. 8 Illustration of morphological operations

Download Full Size | PDF

The effect of morphological reconstruction in our case (Fig. 8) can be seen in by comparing the regions denoted by dash blue rectangles in (e) and (f). By excluding somas (c) from (f), the dendrites are obtained as shown in (g). Lastly the final results with both somas and dendrite traces are obtained by the combination of (c) and the skeletonization of (g). The dilation operation improves the final results in two aspects. Firstly, some dendrites which are just one or two pixels away from somas can be connected to somas. An example is denoted by light blue arrow in (f). Secondly, some false dendrites around soma boundaries can be removed. This effect is visualized by comparing the regions denoted by red dash rectangles in (f) and (g). To summarize, the series of the morphological operations as presented in Fig. 8 have several benefits: firstly, it reduces the false tubular structure; secondly, it guarantees that all dendrites are connected to but not overlap with somas by the exclusion operation.

5. Experiment and discussion

To evaluate the validity of our new method, we conducted two experiments. The first one is based on synthetic images and is intended to demonstrate the feasibility and advantage of our new approach. The second experiment is based on real PCM neuron images of E18 rat cerebral cortical tissue to study the effect of Ivermectin(IVM) to the development of neuron growth [19]. The experiment was approved by Tufts University Institutional Animal Care and Use Committee and complies with the NIH Guide for the Care and Use of Laboratory Animals (IACUC No.B2011 − 45).

There are three parameters, k_m, r_m and N_m, in the airy function formulation (6). All PCM images are collected from rat cortical neuron cultures using a Zeiss Axiovert S100 with an 32× with NA = 0.40 (type LD A-Plan) objective at the resolution of 1728 × 1296. The parameter r_m, the ratio of the inner radius to the outer radius of the phase ring, can be read or measured directly which is 0.85 in our case. The parameters k_m and N_m depend on the hardware specification of the specific microscope and wavelengths of image light [25, 26]. As indicated in [3], wrong parameters generate unsatisfactory restoration results. In addition, restoration results achieve the best using the correct parameters and they are stable around the correct parameters. Thus grid search is used to estimate k_m and N_m by evaluating restoration results on several images manually. In our case, the parameters k_m and N_m are estimated as 0.128 and 15 respectively. The set of parameter values is used for both the synthesized images and the real PCM image data in this section.

5.1. Synthetic image illustration

In this section, we use synthetic images to demonstrate the advantages of our method over Yin’s method in [2, 3] and to justify our use of two level set functions instead of one. Representing both somas and dendrites of the neuron image to be recovered by a Single Level Set (SLS) formulation

θ (x) = s_{1} H_{ε} (ϕ_{1} (x)),

fails to allow for their simultaneous identification. The SLS model can be simplified from the DLS model (7) by assuming s₂ = 0. Thus all calculations and implementations of the SLS model can be obtained from that of the DLS model by setting s₂ = 0.

The synthetic image consists of a bright square and a less bright line as in Fig. 9. The size of square is 30 × 30 and the width of the line is 1 pixel. The square in the true phase retardation image of has the intensity 0.18, the tubular structure has the intensity 0.06 and the standard deviation of the Gaussian noise in the right image of Fig. 9 is 0.02. Here we defined PSNR as

P S N R = 10 \log_{10} \frac{\max (I_{c l e a n})}{M S E (I_{c l e a n}, I_{n o i s y})},

where I_clean and I_noisy are the PCM images with and without noises as displayed in Fig. 9(b) and (c) respectively. In addition, MSE(I_clean, I_noisy) are the mean square error between the noise and clean PCM images.

Fig. 9 Synthesized images. (a): a true phase retardation image; (b): the simulated PCM image without noise; (c): the simulated PCM image with Gaussian noise (PSNR = 15dB)

Download Full Size | PDF

In Fig. 10, the initializations for Yin’s method are obtained directly from the data using the methods in [2, 3] while the initializations for level set formulations are curves obtained using our initialization method in Section 3.4. The radius of the local standard deviation for the square initialization is set to 20 while the standard deviation of the Gaussian template for the steering filter for tubular enhancement ranges from 0.5 to 3 as noise increases. The parameters for weight of the tubular term (16) are set as a = 5, b = 0.2. The parameters in Algorithm 1 are set as follows: r₁ = 15, r₂ = 1, λ₁ = 1, λ₂ = 1 and λ₃ = 0.0015. In order to make a direct comparison, no morphological refinements is used for the synthetic images.

Fig. 10 Results comparison among Yin’s model [2, 3] and the Single Level Set (SLS) and the Double Level Sets (DLS) model for Fig. 9(b) and (c) above. Columns 1-3: Initializations, reconstruction results and segmentation results. Rows 1-3: Yin’s model, SLS and DLS model (top to bottom) for noise free data. Rows 4-6: Yin’s model, SLS and DLS model (top to bottom) for data with noise.

Download Full Size | PDF

For DLS model, the blue curves and the red curves are the initializations for the square and the line structure respectively. For SLS model, all initial curves are blue without distinguishing the square and the line structure and segmentation results only obtain the square. Yin’s method and DLS method both work well in the case without noise. The major advantage of DLS model over Yin’s method is well illustrated for the case with noise. In principle, Yin’s method does an acceptable job for reconstruction by significantly reducing the halo artifacts. In addition, we evaluate the mean square error (MSE) for the reconstructions with respect to different peak signal to noise ratios (PSNR) (20 simulations per PSNR) for Yin’s method, SLS method and DLS method as included in Fig. 11. From Fig. 11, we can see that SLS performs worst for most cases due to the missing of the line structure. In high PSNR cases, both DLS and Yin’s method perform quite well and DLS model outperforms Yin’s method slightly. As PSNR decreases, Yin’s method deteriorates quickly which can also be seen in Fig. 10.

Fig. 11 Mean square error of the reconstruction results

Download Full Size | PDF

For the purpose of image analysis and understanding, the segmentation still needs to be performed on the recovered images for Yin’s method. Moreover, in noisy cases above, simple thresholding skills such as [18] are not capable of recovering the linear structures. While for DLS, no further segmentation is required since the regions {x|ϕ₁(x) > 0} and {x|ϕ₂(x) > 0} already correspond to results for the square and the line respectively.

5.2. Application to real PCM images

In this section, we apply our method on real PCM images of neurons. The parameters for curve initialization for all real PCM images are set as stated in Section 3.4. Through parameters in Algorithm 1, we set the local radii r₁ = 15 and r₂ = 2. The parameters λ₁,λ₂ and λ₃ are set to 5, 80 and 16 respectively. The curve evolution process stops either after 500 iterations or if the region size change of ϕ₁ > 0 and ϕ₂ > 0 is less than 5 pixels between two adjacent iterations. We determine these parameters by manually evaluating the performance of the algorithm using ten images which are not used in the evaluation process below. Finally the parameters for morphological refinements are set as introduced in Section 4. Lastly, the parameters for weight of the tubular term (16) are set as a = 20,b = 0.5.

Results of our proposed method, the approach of Yin et al. in [2, 3] and manual annotation are presented in Fig. 12. In general, Yin’s method (both before and after Otsu-based thresholding [18]) works well in identifying high contrast somas and proves less capable of finding the dendrites. This performance though is to be expected given that the work in [2, 3] was specifically intended to localize cell bodies rather than extended, thin structures such as the dendrites in our data sets. In contrast, our method better recovers the dendrites and yields similar results for the somas. For the dendrites, our method produces more segments than manual results which is not unusual. Actually, even using the same software or algorithm, dendrite tracing results might be different due to different parameter setting or manual interaction. Thus, dendrites segmentation or tracking is usually not evaluated by absolute dendrite length [27, 44]. Later in this section, we present quantitative validation results based on Pearson’s correlation for dendrites as in [27]. In addition, the growth trend of dendrites matches well in the context of the IVM effect evaluation.

Fig. 12 Segmentation results comparison. Columns: different PCM images. Rows (top to bottom): PCM images, Yin’s reconstruction results, Yin’s segmentation results, final results of our method and manual results respectively.

Download Full Size | PDF

In order to quantitatively evaluate our segmentation results, 8 groups of PCM images of E18 rat cerebral cortical tissue were collected. In addition, 4 of these 8 groups are treated with IVM and the other 4 are control groups. Each group consists of 4 images taken at 0 hour, 2 hour, 4 hour, 6 hour and 24 hour for the same sample respectively. One group of the PCM images as well as the segmented results is included in Fig. 13. For the purpose of validation, we manually annotated somas of the 40 PCM images (8 groups multiply 5 images per group) and also dendrites traced manually by using NeuronJ [45]. The manual annotation process takes 5 to 10 minutes to annotate the somas and 30 to 45 minutes to trace the dendrites for one image in a completely hand-operated way. For our approach, no more manual operation is needed except parameter setting. Currently, for a 1728 × 1296 image, our method costs about 30 minutes using unoptimized Matlab codes on a Dell Precision T1600 PC with Quad Core 3.1GHz and 4GB RAM.

Fig. 13 A group of images. Left column: PCM images; Right Column: segmentation results based on our method. Blue areas indicate connected somas, read areas denote isolated somas and yellow lines represent dendrites. Rows (top to bottom): 0, 2, 4, 6 and 24 hour respectively.

Download Full Size | PDF

As in [19], we categorize somas into isolated somas and connected somas in order to study the neuronal growth. The isolated somas are defined as the somas which are not connected to other somas by dendrites and are further than 10 pixels from other somas [19]. As indicated in the right column of Fig. 13, the isolated somas detected by our approach are denoted by red areas. From Fig. 13, we can see that dendrites are growing and connecting more somas as time passes. In addition, the ratio of total soma number to isolated soma number is used as a metric to evaluate weather IVM suppress the neuronal growth as in [19]. Thus we need to perform evaluation on the segmentation of both isolated somas and all somas.

We use several metrics to evaluate soma segmentation. The first metric is called the accuracy (ACC) proposed in [2, 3]. We use P and N to denote soma pixels and background pixels (including dendrites) respectively. The true positive (TP) represents pixels that are both identified manually and our method. The false positive (FP) are pixels labeled by our method but not manually. The false negative (FN) denotes pixels that are segmented manually but not by our method. In addition, |·| represents the number of pixels. Then the accuracy is defined as

A C C = \frac{| T P | + | N | - | F P |}{| P | + | N |} .

The second metric is the Dice coefficient which is defined as

D I C E = \frac{2 | T P |}{(| F P | + | T P |) + (| T P | + | F N |)} .

The ACC and DICE values of all somas and isolated somas for all 40 images are included in Fig. 14. We can see that the accuracies for all somas and the isolated somas are above 0.97 and 0.96 respectively. In addition, most of the dice coefficients range from 0.7 to 0.8 and 0.5 to 0.7 for all somas and isolated somas respectively. We see that the performance for isolated somas is a little worse. This though is not unexpected as identifying isolated somas is dependent on dendrite segmentation. Both false positive dendrites and false negative dendrites might cause incorrect isolated somas classification. Thus incorrect dendrite segmentation might degrade isolated soma segmentation results. Still all these values in Fig. 14 indicate that our method generates rather similar results as manual processing with very little manual effort.

Fig. 14 Quantitative comparison between the proposed method and manual results for somas. Left: ACC results; Right: Dice results.

Download Full Size | PDF

In addition to ACC and DICE which are defined by comparing pixels, we also get the soma numbers for validation as in [27]. From Fig. 15, we can see that our method is consistent with manual results for both isolated somas and total somas at all time points. The total soma numbers almost remain unchanged and the trend of isolated soma numbers is decreasing through different times suggesting that dendrites grow and connect to more somas as time goes by.

Fig. 15 Soma numbers comparison.

Download Full Size | PDF

For dendrites evaluation, we also present dendrite length according to the control and IVM group at different time points in Fig. 16. We can see that our methods yields longer dendrites for both the control group and drug group at all time points. As mentioned in [27], different methods might generated different results for dendrites length. Thus the correlation is always used to validate dendrite segmentation [27]. In addition, the measured dendrites length using the proposed method as well as manual results for all 40 images are included in Fig. 16. The correlation between the two sets of results are 0.94 indicating the consistency of our method and manual results. Moreover, we can find that for both our method and manual results at all time points, the control group has fewer dendrites than that of the drug group suggesting that IVM suppresses the growth of neural cells.

Fig. 16 Quantitative comparison between the proposed method and manual results for dendrites. Left: The x axis represents the manually traced dendrite length; the y axis denotes the segmented dendrites length using the proposed method. Each dot represents a results for each PCM image. Right: Comparison for dendrites at different time points and different groups

Download Full Size | PDF

Finally, we check the ratio of all soma numbers to isolated soma numbers as displayed in Fig. 17. As already indicated in [19], manual results suggest IVM suppresses the neural growth and connectivity. Here we use repeated measures ANOVA test [46] to accommodate the drug and control group and the five time points. The statistical test result indicates that there is a significant (p = 0.0031) difference between the drug and control groups for manual results. For our method, the repeated measures ANOVA test also indicates that the ratios in the control group are significantly larger than that of the drug group (p = 0.0137).

Fig. 17 Ratio of all soma numbers to isolated soma numbers. Left: manual results; Right: proposed method.

Download Full Size | PDF

In summary, experimental results demonstrate that our method produces very similar segmentation results with manual annotation for somas. For dendrites, our method yields consistent results with manual results. In addition, our method can be applied to investigate IVM effect by investigating the significant difference of the ratio of all soma numbers to isolated soma numbers in both control group and drug group as time passes.

6. Conclusion and future work

In this paper, we introduce a parametric image model using the level set framework for the purpose of the PCM image restoration/segmentation. Then we present an optimization problem merging both image restoration and segmentation and give its solutions. Based on the segmentation method just mentioned, we introduce a pipeline for PCM neuron image analysis which is automatic and the first to segment somas and trace dendrites simultaneously for the PCM neuron images. The technical advantages of our method have been well demonstrated by comparing to other state of the art methods using both synthetic and real PCM images. Moreover, our approach is shown to be feasible for the study of IVM effect to neural growth and development.

We are going to build a regularization term to couple the somas and dendrites in the variational formulation (8). Usually for the multi-level set formulation, constraints need to be imposed to avoid overlap [30,47–49]. In the neuron application, no clear boundaries exist between the soma and dendrites for a neuron. Normal overlap penalty term such as [30] would separate somas and dendrites which is not we want for the neural connectivity evaluation. Thus building a regularization term which can avoid overlap between somas and dendrites and also connect them is an interesting problem and also one of the future works.

Currently, the coefficient s₁ and s₂ in (7) are two scalars which means the phase retardation image θ (x) is a piecewise constant model [22]. By extending s₁ and s₂ to be spatial variant fields s₁(x) and s₂(x), θ (x) can be extended to a piecewise smooth model [47] which could better model neural cells over the whole image domain. Moreover, we may also investigate the integration of machine learning techniques such as in [50] in modeling θ (x) within the level set framework.

One important future direction is to develop methods for the optimal selection of the parameters such as λ₁, λ₂ and λ₃ to further reduce manual interaction. Another important area of future work should be devoted to accelerate our method via: using faster programming languages such as C++; using fast level set implementation techniques such as [51]; using parallel computing techniques such as Graphics Processing Unit (GPU).

Phase contrast microscopy belongs to a type of microscopy called partially coherent microscopies [25,26] which does not need fluoresce or staining for specimens. Besides phase contrast microscopy, dark field and differential interference contrast (DIC) microscopy are also examples of spatially partially coherent microscopies. A key issue here is the difference in the point spread functions among these microscopes which in turn give rise to data-domain artifacts that also differ from those seen in this study. While we believe that the basic, two-level set approach pursued here can be used in the context of these other microscopes, changes and adaptation will certainly be required to apply these ideas to other cell lines as well as partially coherent modalities.

Acknowledgments

The authors would like to thank the reviewers for their insightful comments that help improve the manuscript significantly. This study was supported by German Research Council (DFG, OE 541/2-1), National Institutes of Health (NIH) via grants R01 AR005593 and R01 AR061988 and National Science Foundation (NSF) award 1347191. Office of Electricity Delivery and Energy Reliability (OE)

References and links

1. F. Zernike, “How I discovered phase contrast,” Science 121, 345–349 (1955). [CrossRef] [PubMed]

2. Z. Yin, K. Li, T. Kanade, and M. Chen, “Understanding the Optics to Aid Microscopy Image Segmentation,” MICCAI , 6361209–217 (2010). [PubMed]

3. Z. Yin, T. Kanade, and M. Chen, “Understanding the phase contrast optics to restore artifact-free microscopy images for segmentation,” MedIA 16, 1047–1062 (2012).

4. K. Li, E. D. Miller, L. E. Weiss, P. G. Campbell, and T. Kanade, “Online tracking of migrating and proliferating cells imaged with phase-contrast microscopy,” IEEE CVPR workshop, 65 (2006).

5. I. Ersoy, F. Bunyak, K. Palaniappan, M. Sun, and G. Forgacs, “Cell spreading analysis with directed edge profile-guided level set active contours,” MICCAI , 5241376–383 (2008). [PubMed]

6. C. J. Bradhurst, W. Boles, and Y. Xiao, “Segmentation of bone marrow stromal cells in phase contrast microscopy images,” IEEE Inter. Conf. Image Vis. Computing, 1–6 (2008).

7. M. E. Ambuhl, C. Brepsant, J. Meister, A. B. Verkhovsky, and I. F. Sbalzarini, “High-resolution cell outline segmentation and tracking from phase-contrast microscopy images,” J. Microscopy 245, 161–170 (2012). [CrossRef]

8. A. Can, H. Shen, J. N. Turner, H. L. Tanenbaum, and B. Roysam, “Rapid automated tracing and feature extraction from retinal fundus images using direct exploratory algorithms,” IEEE Trans. Inf. Tech. Biomed. 3, 125–138 (1999). [CrossRef]

9. R. F. Frangi, W. J. Niessen, K. L. Vincken, and M. A. Viergever, “Multiscale vessel enhancement filtring,” MICCAI , 1496130–137 (1998).

10. L. M. Lorigo, O. D. Faugeras, W. E. L. Grimson, R. Kervien, R. Kikinis, A. Nabavi, and C. F. Westin, “CURVES: Curve evolution for vessel segmentation,” MedIA 5, 195–206 (2001).

11. K. Krissian, “Model based detection of tubular strucutres in 3D images,” Comput. Vis. Image Und. 80, 130–171 (2000). [CrossRef]

12. A. Vasilevskiy, “Flux maximizing geometric flows,” IEEE Trans. PAMI 24,1565–1578 (2002). [CrossRef]

13. J. E. Vaughn, “Fine structure of synaptogenesis in the vertebrate central nervous system,” Synapse 3, 255–285 (1989). [CrossRef] [PubMed]

14. G. L. Gerstein, P. Bedenbaugh, and M. H. Aertsen, “Neuronal assemblies,” IEEE Trans. Biomed. Eng. 36, 4–14 (1989). [CrossRef] [PubMed]

15. E. Meijering, “Neuron tracing in perspective,” Cytometry Part A 77, 693–704 (2010). [CrossRef]

16. O. Al-Kofahi, R. J. Radke, B. Roysam, and G. Banker, “Automated semantic analysis of changes in image sequences of neurons in culture,” IEEE Trans. on Biomed. Eng. 53, 1109–1123 (2006). [CrossRef]

17. J. S. Vestergaard, A. L. Dahl, P. Holm, and R. Larsen, “Pipeline for tracking neural progenitor cells,” MICCAI Workshop WCV , 7766155–164 (2012).

18. N. Otsu, “A threshold selection method from gray-level histograms,” IEEE Trans. Sys. Man. Cyber. 9, 62–66 (1979). [CrossRef]

19. N. Ozkucur, K. P. Quinn, J. Pang, I. Georgakoudi, E. Miller, M. Levin, and D. L. Kaplan, “Membrane potential depolarization causes alterations in neuron arrangement and connectivity in cocultures,” Brain and Behav. 5, e00295 (2005).

20. S. J. Osher and R. P. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces(Springer-Verlag, 2002).

21. J. A. Sethian, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science(Cambridge University Press, 1999).

22. T. F. Chan and L. A. Vese, “Active contours without edges,” IEEE Trans. Image Process. 10, 266–277 (2001). [CrossRef]

23. H. Sagan, Introduction to the Calculus of Variations (Courier Dover Publications, 2012).

24. D. Murphy, Fundamentals of Light Microscopy and Electronic Imagings(Wiley, 2001).

25. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interfence and Diffraction of Light(Cambridge University Press, 1999). [CrossRef]

26. S. B. Mehta and R. Oldenbourg, “Image simulation for bilogical microscopy: microlith,” Biomed. Opt. Express 5, 1822–1838 (2014). [CrossRef] [PubMed]

27. S-Y Ho, C-Y Chao, H-L Huang, T-W Chiu, P. Charoenkwan, and E. Hwang, “NeurphologyJ: An automatic neuronal morphology quantification method and its application in pharmacological discovery,” BMC Bioinformatics 12, 1471–2105 (2011). [CrossRef]

28. S. Lankton and A. Tannenbaum, “Localizing region-based active contours,” IEEE Trans. Image Process. 17, 2029–2039 (2008). [CrossRef] [PubMed]

29. A. Gooya, H. Liao, K. Matsumiya, K. Masamune, Y. Masutani, and T. Dohi, “A variational method for geometric regularization of vascular segmentation in medical images,” IEEE Trans. Image Process. 17, 1295–1312 (2008). [CrossRef] [PubMed]

30. K. Zimmer and J. Oliver-Marin, “Compled Parametric Active Contours,” IEEE Trans. on PAMI 27, 1838–1842 (2005). [CrossRef]

31. C. R. Vogel, “Computational methods for inverse problems,” SIAM (2002).

32. Y. Wang, J. Yang, W. Yin, and Y. Zhang, “A new alternating minimization algorithm for total vairation image reconstruction,” SIAM J. Imaging Sciences 1, 248–272 (2008). [CrossRef]

33. C. Li, C.Y Kao, J.C. Gore, and Z. Ding, “Minimization of region-scalable fitting energy for image segmentation,” IEEE Trans. Image Process. 17, 1940–1949 (2008). [CrossRef] [PubMed]

34. A. Yezzi, A. Tsai, and A. Willsky, “A fully global approach to image segmenation via coupled curve evolution equation,” J. Vis. Comm. Image Rep. 13, 195–216 (2002). [CrossRef]

35. V. Caselles, R. Kimmel, and G. Sapiro, “Geodesic active contours,” IJCV 22, 61–79 (1997). [CrossRef]

36. K. Wu, D. Gauthier, and M. Levine, “Live cell image segmentation,” IEEE Trans. on Biomed. Eng. 42, 1–12 (1995). [CrossRef]

37. W. Freeman and E. Adelson, “The design and use of steerable filters,” IEEE Trans. PAMI 13, 891–906 (1991). [CrossRef]

38. J. Kittler and J. Illingworth, “Minimum error thresholding,” Pattern Recogn. 19, 41–47 (1986). [CrossRef]

39. M. K. Ng, R. H. Chan, and W. C. Tang, “A fast algorithm for deblurring models with Neuman boundary conditions,” SIAM J. Sci. Comput. 21, 851–8666 (1999). [CrossRef]

40. S. Serra-Capizzano, “A note on Antireflective boundary Conditions and fast deblurring models,” SIAM J. Sci. Comput. 25, 1307–1325 (2004). [CrossRef]

41. D. Nam, J. Mantell, D. Bull, P. Verkade, and A. Achim, “A novel framework for segmentation of secretory granules in electron micrographs,” MedIA 18, 411–424(2014).

42. R. C. Gonzalez and R. E. Woods, Digital Image Processing(Prentice Hall, 2007).

43. R. C. Gonzalez, R. E. Woods, and S. L. Eddins, Digital Image Processing using MATLAB (Gatesmark Publishing, 2009).

44. L. Dehmelt, G. Poplawski, E. Hwang, and S. Halpain, “NeuriteQuant: An open source toolkit for high content screens of neuronal morphogenesis,” BMC Neurosci. 12, 100 (2011). [CrossRef] [PubMed]

45. E. Meijering, M. Jacob, J.-C.F. Sarria, P. Steiner, H. Hirling, and M. Unser, “Design and validation of a tool for neurite tracing and analysis in fluorescence microscopy images,” Cytometry Part A 58A, 167–176 (2004). [CrossRef]

46. C. S. Davis, Statistical Methods for the Analysis of Repeated Measurements(Springer, 2002).

47. L. Vese and T. Chan, “A multiphase level set framework for image segmentation using the Mumford and Shah model,” IJCV 50, 271–293 (2002). [CrossRef]

48. A. Vazquez-Reina, E. Miller, and H. Pfister, “Multiphase geometric couplings for the segmentation of neural processes,” IEEE Conf. CVPR, 2020–2027 (2009).

49. J. Pang, J. B. Driban, T. E. McAlindon, J. G. Tamez-Pena, J. Fripp, and E. L. Miller, “On the use of directional edge information and coupled shape priors for segmentation of magnetic resonance images of the knee,” IEEE J. of Biomed. and Health Inform. 19, 1153–1167 (2015).

50. H. Su, Z. Yin, S. Huh, and T. Kanade, “Cell segmentation in phase contrast microscopy images via semi-supervised clustering over optics-related feataures,” MedIA 17, 746–765 (2013).

51. Y. Shi and W. C. Karl, “A real-time algorithm for the approximation of level-set-based curve evolution,” IEEE Trans. Image Process. 17, 645–656 (2008). [CrossRef] [PubMed]

Automatic neuron segmentation and neural network analysis method for phase contrast microscopy images

Abstract

1. Introduction

2. The level set method for image segmentation

3. A variational formulation for PCM neuron image segmentation

3.1. A data fidelity term based on the physical model

3.2. Localized active contours

3.3. A weighted geometric regularization of tubular structure

3.4. Curve initialization

3.5. Implementation

4. Morphological refinement

5. Experiment and discussion

5.1. Synthetic image illustration

5.2. Application to real PCM images

6. Conclusion and future work

Acknowledgments

References and links

Cited By

Figures (17)

Tables (1)

Equations (39)

Biomedical Optics Express