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

: Phase Contrast Microscopy (PCM) is an important tool for the long term study of living cells. Unlike ﬂuorescence methods which suffer from photobleaching of ﬂuorophore or dye molecules, PCM image contrast is generated by the natural variations in optical index of refraction. Unfortunately, the same physical principles which allow for these studies give rise to complex artifacts in the raw PCM imagery. Of particular interest in this paper are neuron images where these image imperfections manifest in very different ways for the two structures of speciﬁc interest: cell bodies (somas) and dendrites. To address these challenges, we introduce a novel parametric image model using the level set framework and an associated variational approach which simultaneously restores and segments this class of images. Using this technique as the basis for an automated image analysis pipeline, results for both the synthetic and real images validate and demonstrate the advantages of our approach.


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. 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][9][10][11][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 (voltagesensitive 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.

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 C as the intersection of a horizontal plane z(x) = 0 of an auxiliary level set function φ (x); i.e., 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( C) based on C, a curve representing the boundary of the object. The energy functional E( C), consisting of a data fidelity term and a regularization term, is always constructed in a way that when C approaches to the true boundary of the object, the minimum value of E( 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 where I(x) is the image intensity value at location x, Ω 1 is the region inside of the curve C, Ω c 1 is the region outside of the curve C and Ω = Ω 1 ∪ Ω c 1 which is the whole image domain. In addition, c 1 and c 2 are the mean image intensity values for the regions Ω 1 and Ω c 1 respectively. Moreover, α is a coefficient, | C| is the length of the curve C and the term α| C| basically is a smoothness regularization term.
By introducing the level set function φ such that C = {x|φ (x) = 0}, eqn. (1) can be converted into a function based on φ as where H ε (·) is an approximate Heaviside function [22], defined as and the Heaviside function is defined as By using calculus of variations [23], the optimization of E cv can be realized by iteratively updating c 1 , c 2 and φ using the following Euler-Lagrange equation where 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 → ∞, C(t) approaches the contour of interest in the image.

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 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] airy where J 1 (·) 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 φ 1 and φ 2 for somas and dendrites respectively. Specifically, the image to be recovered can be represented as the following Double Level Sets model (DLS) where φ 1 and φ 2 are the level set functions, s 1 and s 2 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 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.

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 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: C n×m → C nm×1 such that for a given array w ∈ C n×m , vec(w) = [w 1,1 , · · · , w n,1 , w 1,2 , · · · , w n,2 , · · · , w 1,m , · · · , w n,m ] T .
The operator array defines the inverse of the vec operator such that array(vec(w)) = w, vec(array(w)) = w, where w ∈ C n×m and w ∈ C nm×1 .
Then we can define φ φ φ 1 , φ φ φ 2 and θ θ θ are the vectorized version of φ 1 , φ 2 and θ respectively. Thus the vectorized version of the DLS model (7) is represented by where φ φ φ 1 = vec(φ 1 ), φ φ φ 2 = vec(φ 2 ) and θ θ θ = vec(θ ). In addition, (9) can be discretized and represented by the following matrix form 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-circulantblock(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

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 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]. Specifically, the local image fitting energy function formulation of [28] is defined as follows 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 where || · || 2 indicates the l 2 norm. In this paper, we use a localized version of the mean separation energy [34] for F loc 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

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. 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. where Tr(·) represents the trace of a matrix, f is a non-negative decreasing function from R + to R + . 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 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. 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.
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).

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

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 (φ 1 , φ 2 , s 1 , s 2 ) with respective to φ k
Step 4: Update s 2 by numerically solving (24). and s k (k = 1, 2) can be obtained as Using the array(·) operator defined in Section 3.1, the array version of (17) is According to [28], the gradient descent flow of E loc defined in (12) is For the weighted tubular regularization term E wtub , according to [29], the gradient flow of (15) is andḟ 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 2 × n 2 .
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 2 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 2 ). Alternative choices for modeling K and the resulting computational methods can be found in [31,39,40]. For the two real scalars s 1 and s 2 , based on (18), s 1 is updated by numerically solving the following minimization problem (assuming φ 1 , φ 2 , s 2 fixed) arg min Here we set ub = 0.2 which is a upper bound of s 1 . Next, update s 2 by numerically solving the following minimization problem (assuming φ 1 , φ 2 , s 1 fixed) arg min By using the constraint 0 < s 2 ≤ s 1 , the level set function s 1 would catch the brighter objects as shown in Section 5.1. We update s 1 and s 2 once for every ten updates to the level set functions. This change reduces computation load and makes little difference to the final results.

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

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.

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 fails to allow for their simultaneous identification. The SLS model can be simplified from the DLS model (7) by assuming s 2 = 0. Thus all calculations and implementations of the SLS model can be obtained from that of the DLS model by setting s 2 = 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 PSNR = 10 log 10 max(I clean ) MSE(I clean , I noisy ) , 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. 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 1 = 15, r 2 = 1, λ 1 = 1, λ 2 = 1 and λ 3 = 0.0015. In order to make a direct comparison, no morphological refinements is used for the synthetic images.
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. 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|φ 1 (x) > 0} and {x|φ 2 (x) > 0} already correspond to results for the square and the line respectively.

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 1 = 15 and r 2 = 2. The parameters λ 1 , λ 2 and λ 3 are set to 5, 80 and 16 respectively. The curve evolution process stops either after 500 iterations or if the region size change of φ 1 > 0 and φ 2 > 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] 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. 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.
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 The second metric is the Dice coefficient which is defined as .
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.
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.
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 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).
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.

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][48][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 1 and s 2 in (7) are two scalars which means the phase retardation image θ (x) is a piecewise constant model [22]. By extending s 1 and s 2 to be spatial variant fields s 1 (x) and s 2 (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 λ 1 , λ 2 and λ 3 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.