Multiphase geometric couplings for the segmentation of neural processes

The ability to constrain the geometry of deformable models for image segmentation can be useful when information about the expected shape or positioning of the objects in a scene is known a priori. An example of this occurs when segmenting neural cross sections in electron microscopy. Such images often contain multiple nested boundaries separating regions of homogeneous intensities. For these applications, multiphase level sets provide a partitioning framework that allows for the segmentation of multiple deformable objects by combining several level set functions. Although there has been much effort in the study of statistical shape priors that can be used to constrain the geometry of each partition, none of these methods allow for the direct modeling of geometric arrangements of partitions. In this paper, we show how to define elastic couplings between multiple level set functions to model ribbon-like partitions. We build such couplings using dynamic force fields that can depend on the image content and relative location and shape of the level set functions. To the best of our knowledge, this is the first work that shows a direct way of geometrically constraining multiphase level sets for image segmentation. We demonstrate the robustness of our method by comparing it with previous level set segmentation methods.


Introduction
Deformable models based on level sets have been successfully applied to a variety of computer vision tasks such as image segmentation and video tracking over the last ten years [25,3,7,8].Their success is mostly attributed to their parametrization-free nature, intuitive formulation, and ability to easily adapt to shapes of unknown topology [13].
The problem of image segmentation (partitioning) within this framework is usually cast in a variational formulation; an energy functional is defined on the space of possible contours or image partitions (phases), and the geometric de- Figure 1: Image segmentation of a scene with two objects.Subfigures (a) and (b) show a deformation of Object I. Current multiphase level set methods for image segmentation do not allow for the addition of prior knowledge about the geometric arrangement of the partitioning.However, in this example we might know a priori that region B should always surround region A, and that the partitions C, and D form a different object that should not ever be surrounded by B. This paper introduces a way of grouping and inducing such geometrical arrangements in the partitioning.formable model is then iteratively evolved until an optimal solution is found.
A common approach to constrain the geometry of a deformable model is to build shape priors that are statistically learned from a set of training templates [15,5,10].However, such priors are limited to the subspace of learned deformations from the training set (typically up to an affine or a projective transformation), and they cannot directly model sophisticated geometric arrangements of deformable models.Such ability is important in many imaging applications.For example, as we discuss further in this paper, the segmentation of cellular and intracellular membranes in biomedical imaging often requires partitioning the image into multiple nested contours.Such prior information about the image geometry can be used to avoid undesired segmentations and improve the overall segmentation accuracy.
Here, we introduce a way to directly model geometric objects that are naturally described using multiple level set functions.As opposed to most of the previous work on multi-shape learning for deformable models [20,12], our method does not rely on the statistical inference of a multishape distribution from a set of training samples.We instead provide a way to directly design entire families of partition arrangements using dynamic force fields that can depend on the image content and multiple level set functions.There are a number of important application areas where such an approach is useful.Consider the example of Fig. 1 that shows a scene with two objects and the background.A classical multiphase method for segmentation would partition the scene into several phases according to image features (gradient, statistical variance, etc.), and optionally, according to some statistically learned shape prior for each of the partitions.Such approaches do not allow for the addition of prior knowledge about the relative arrangement of the regions in the partitioning.A real example of this idea can be seen when tracking neural processes (cellular and intra-cellular boundaries) in a sequence of sections from a 3D volume of brain tissue.In each section, cells, mitochondria and other intra-cellular objects have a membrane of homogeneous intensity and varying thickness (see an example in Fig. 2).The tracking of these processes is of important interest in neuroscience, where scientists are aiming at identifying and analyzing the internal wiring of the brain.
The main contributions of our work are first, the definition of elastic couplings between level set functions using dynamic force fields that can be easily integrated in any variational-based multiphase framework, and second, the modeling of ribbon-like deformable models that can be used to segment and track neural processes.To the best of our knowledge none of these issues have been addressed to date.

Related work
Variational and energy minimization models for image segmentation based either on level sets or graph cuts are well documented in the literature [25,2,22].For multiobject segmentation, multiway graph cuts and multiphase graph cuts and level sets provide a natural extension of the single object case [22,19,18,21].In the multiphase level sets literature, much attention has gone into the study of topologically constrained flows that avoid vacuum regions and overlap between the different phases [3,17,26,24], and into extending shape priors for single level set evolutions to the multiphase case [6,4,9,12].Two interesting problems have been addressed in this last area, first, identifying which shape prior should be applied to which level set function (this process is also known as shape prior competition) [6,4], and second, applying shape priors while allowing for mergers and splits of multiple phases [9].
Very recently, there has been an effort to employ exter-Figure 2: Example of a cross section of a 3D volume of brain tissue acquired with an electron microscope (1px.= 3nm x 3nm).Given the nature of the images, we can safely make the assumption that sections can be decomposed into a set of "ribbons" of varying thickness.
nal force fields within the multiphase level set framework in a manner which preserves the underlying, assumed known, topology of the problem [7].Our work shares some similarities with this work.Mostly, we both consider the integration of external force fields in multiphase level sets and their application in segmentation.However, our work focuses on how to induce a geometrical arrangement on the different phases, while the work of [7] focuses on preserving the topology of the different phases and avoiding vacuum regions and overlap altogether.

Minimum partition with multiphase level sets
To motivate our work, we consider the classical two dimensional grayscale piecewise-constant segmentation problem in computer vision known as the minimal partition problem or the Mumford and Shah problem.This problem was originally formulated in [16], and since then, several variations have appeared in the literature.In the following we introduce a variation of the so called reduced model [21] of this problem: Let Ω ∈ R 2 be open and bounded, then, given an observed image u, we seek a decomposition G = (Ω 1 , . . ., Ω m ) of Ω and a vector of constants c = (c 1 , . . ., c m ) such that the following energy is minimized: Figure 3: (a): We make the assumption that any cross section of the brain tissue can be decomposed into a set of "ribbons".(b): In the simplified case where we only have one ribbon, the ribbon partitions the image into three different phases.According to our multiphase encoding of Eq. ( 2), we need at least two distance field functions φ 1 and φ 2 to represent three phases.
where Γ i represents the boundary set of each partition Ω i and the tuning parameters λ i ≥ 0 and α i ≥ 0 weigh the relative importance of the different terms in the energy.The above energy favors intensity smoothness along each partition and penalizes the size of their boundary sets.
In order to translate the above problem into a multiphase formulation, we use n = log 2 m level set functions and follow the implementation in [21], which guarantees no vacuum and overlap between phases.In this framework, a vector level set function Φ is defined as Φ = (φ 1 , . . ., φ n ) where each φ i : Ω → R is a level set function.Similarly, a partition Ω i is represented by a binary vector K i of length n.This way, we can define the characteristic function χ i for each partition Ω i as follows: where (H (φ 1 ) , . . ., H (φ n )) is a binary vector that takes different values when evaluated in different partitions of the image.The function H (φ) takes the values H (φ) = 0 when φ ≤ 0 and H (φ) = 1 when φ > 0.
We can now rewrite Eq. ( 1) in terms of the level set functions as: (3) Notice that our level set representation uses the same number of level set functions as [21] and [6], and is different from the work of [3,9,19], where one level set function per partition was used instead, and from [7] where only a total of four functions were needed.The Euler-Lagrange equation corresponding to the gradient descent of the functional in Eq. (3) yields a system of n evolution equations for (φ 1 , . . ., φ n ) (see [21] for an example of a four-phase system).Figure 4: Force field for ribbon consistency.The iso-lines of φ 1 and φ 2 are depicted in blue and red, respectively.Depending on the sign of σ 1 and σ 2 , the force field would attract or repel the boundaries of the ribbon.In the case displayed, both σ 1 and σ 2 are considered to be positive.

Adding geometric priors to the multiphase framework
The variational formulation of Eq. ( 1) only uses what are conventionally called internal terms (those associated with the regularity of the interface), and data terms (those associated with the image data).There is nothing in such formulation that constrains the relative shape and positioning of the partitions.In the present section, we present a way of controlling the geometric arrangement of the partitioning by coupling several phases using dynamic force fields.These force fields will generate velocities for each partition that will be added to those generated by the gradient descent of the Mumford Shah (MS) functional.First, we recall the relationship between curve evolution, level sets and their connection with force fields.
The Level Set Method allows to connect the propagation of a 2D front γ to the evolution of the zero level set of a function φ (γ, t).This way, γ propagates with a speed γ t ∈ T γ M , if and only if, by the chain rule, φ propagates according to the Level Set equation: where T γ M ∈ R 2 is the tangent space of the manifold M of closed curves immersed in R 2 , defined at γ [14].
The n evolution equations for (φ 1 , . . ., φ n ) that result from the gradient descent of the MS functional in Eq. ( 3) propagate the zero-level sets of each φ function with speeds γ M S 1t , . . ., γ M S nt due to the connection established by the Level Set Method.Since the encoding of each partition in Eq. ( 2) links the evolution of each level set function with the evolution of each partition, the gradient flow also implicitly propagates each partition (Ω 1 , . . ., Ω m ) with speeds that we denote as Γ M S 1t , . . ., Γ M S mt .These velocities result from the optimization of the minimum partition problem, and therefore have a variational nature.Consider now a vector field v : Ω → R 2 .We can build a ve-Figure 5: σ 1 as a function of φ 2 The graph for σ 2 is obtained by reversing the φ axis (σ 2 (φ) = σ 1 (−φ)).The region labeled as A corresponds to the attraction (σ 1 ≥ 0, v 1 ≥ 0) between the ribbon boundaries (stronger as the boundaries separate from each other).Region B is the region where the ribbon shows a plastic behavior (no resistance towards deformation).Region C corresponds to the repulsion (σ 1 ≤ 0, v 1 ≤ 0) between the ribbon boundaries (stronger as the boundaries come closer to each other).locity field γ v it for the zero-level set of one of our functions Such mapping extracts the normal component of v to the zero level set of φ i .As in the Mumford Shah case, the vector field v implicitly also maps into a vector of velocities (Γ v 1t , . . ., Γ v mt ) for each of the boundary sets of the partitions.
We can extend this concept to build a force field F = (v 1 , . . ., v n ) for each zero-level set γ i .In the more general case, we can generate p force fields that, if designed wisely, could be used to arrange the different partitions (Ω 1 , . . ., Ω m ) on the plane.Since T γ M is a vector space, the velocity fields derived from such force fields can be added to the velocities derived from the MS functional as follows: where γ t = (γ 1t , . . ., γ nt ) is the vector of total velocities of the zero-level sets, γ M S t = γ M S 1t , . . ., γ M S nt is the vector of velocities given by the gradient descent of the MS functional, and is the vector of velocities given by the action of the force field F j = (v 1j , . . ., v nj ) onto each of the φ functions.The parameters µ j determine the strength of the force fields relative to that arising from the MS functional.
It is important to note, however, that since the vector fields v don't have to be irrotational (curl-free), some of them might not equal the gradient of a scalar potential, and therefore it is not always possible to guarantee the existence of an equivalent variational formulation for each vector field.For this reason, we will consider that a solution for the segmentation problem is found (the evolution finishes) when the velocities from the optimization of the MS functional and those from the external force Figure 6: We want the outer side of the active ribbon (φ 2 = 0) to be attracted towards the outer side of the cellular membrane (feature map in green), and the inner side of the ribbon (φ 1 = 0) towards the inner cellular membrane (feature map in pink).fields balance each other and an equilibrium is reached (||(γ 1t , . . ., γ nt )|| L 2 ≤ ).In the best case, the velocities from the MS functional balance their counterparts from the force fields F j : In the next section we show several examples of force fields that can be used to induce geometrical arrangements in the partitioning.

Active Ribbons
Consider again the problem of segmenting and tracking neural processes presented in the introduction.Neural membranes can be composed of myelin in the case of axons, and the analysis of their thickness can reveal important information about the connectivity of biological neural networks and neural-related diseases [1].For this reason it would be useful to be able to segment and extract each of the membranes in an isolated partition.This allows us to use the ideas from the previous section to model such geometric partitioning.
We start by defining an active ribbon as the deformable region between two non-intersecting contours, one contained within the other.Figure 3(a) shows that a single active ribbon yields a partitioning of the image in three regions (inside, ribbon and outside).According to our multiphase encoding of Eq. (2), we need at least two distance field functions φ 1 and φ 2 to represent three regions.
We now consider three force fields that can be used to arrange the image partitions into a set of ribbons.The first two forces control the shape of each ribbon and their ability to find the right cellular boundaries in an image of brain tissue.The third one is required for problems where we wish to track simultaneously multiple neurons and models the interaction between ribbons by controlling their mutual repulsion or attraction.

Force field for ribbon consistency
Consider the force field: where v 1 and v 2 act on φ 1 and φ 2 , respectively, and σ 1 and σ 2 are scalar fields defined on the image plane.The joint action of v 1 and v 2 creates a repulsion or an attraction force between the boundaries of the ribbon depending on the sign of σ 1 and σ 2 (see Fig. 4).Ideally, we want to design the ribbon so that it shows plasticity (no resistance to deformation) when the thickness of the ribbon is within some reasonable range (5-20 nm.).In such cases, the evolution of the ribbons would be mostly driven by the gradient flow of the MS functional Eq. ( 1).On the other hand, we want to trigger the repulsion or attraction between the boundaries of the ribbon when the ribbon has an abnormal thickness.Following this reasoning, we can design σ 1 and σ 2 for φ 1 and φ 2 so that they react to the proximity between the boundaries of the ribbons.We can achieve this effect by setting σ 1 (φ 2 ) and σ 2 (φ 1 ) with a profile such as the one depicted in Fig. 5.A piecewise polynomial approximation of this profile is discussed in Section 6.

Force field for ribbon-cell interaction
The previous force field model adds a purely geometric force to the multiphase partitioning by inducing the creation of elastic ribbons.The combination of such force field with the MS functional will segment the image in homogeneous ribbon-like partitions.However, the model so far does not necessarily guarantee that the two different boundaries of the ribbon will be attracted to different boundaries of the cellular membranes (see Fig. 6).In this section we introduce a second force field that, when added to the previous model, achieves precisely that.
We start by recalling the vector field convolution model The vector field points towards the inner side of the cellular membranes (b): User defined vector kernel for vector field convolution.See [11] for more details on different types of kernels.
introduced in [11].Given a feature map defined on the image e : Ω → R + (i.e. a bitmap) we can build a smooth vector field v (v x , v y ) : Ω → R 2 that points towards the highest values in e as: where k(k x , k y ) is a 2D user-defined vector field kernel such as the one shown in Fig. 8b, and the subscripts refer to the x and y components of each vector field.See [11] for additional notes on kernel selection.Building on this idea, consider the following force field: where I σ is a smoothed version of the image I, and e 1 = ∇I σ • ∇φ 1 and e 2 = −∇I σ • ∇φ 2 are the feature maps for φ 1 and φ 2 , respectively.The above force field is made of a vector field v 1 that acts on φ 1 by pushing its zero-level set towards the inner side of the cellular boundaries, and v 2 that acts on φ 2 by pushing its zero-level set towards their outer side (see Fig. 6).Figures 7 and 8(a) show the feature maps and the resulting force field on examples of real images.

Interaction between ribbons
The previous two force fields control the shape and segmentation of each ribbon in an individual manner, and therefore they do not offer control over the interaction of neighboring structures.However, going back to our application of interest (tracking of neural boundaries), it is known that when moving from one section to the next one through the 3D volume of brain tissue, the cellular boundaries of the neurons should not change their relative position abruptly.That is, if two neurons were adjacent to each other and/or relatively close in one section, they should be so in the next one.We can take advantage of this fact and control the geometric arrangement of the partitioning further to avoid undesired segmentations and speed up the convergence of the multiphase evolution.
Consider that while processing section i of the volume of brain tissue we are given the segmentation results from section i − 1.Since we know the location of each ribbon in section i, we can build a force field F ab (v 1 , v 2 , v 3 , v 4 ) for each pair of ribbons a and b, such that: where the vectors ab i and ab i−1 are the vector that points from the center of mass of a to the one of b, and σ I : R → R is a function that controls the strength of the mutual repulsion (σ I < 0) or attraction (σ I > 0) between the two ribbons (see Fig. 9).Section 6 discusses a piecewise polynomial approximation of σ I .

Experiments
In this section we present several experiments that show the robustness of our active ribbons model on real data.First, we compare our model with three other well known level set-based deformable models for the segmentation of single cellular boundaries.We then compare our method with a geometrically-unconstrained multiphase level set model for the segmentation multiple cellular boundaries.Finally, we show how our active ribbon model can efficiently track cellular boundaries on a sequence of crosssections obtained from a volume of brain tissue and give an example of failed segmentation.
Figure 10 shows a comparison of our model with several geometrically-unconstrained deformable models.In the best cases, these models were able to accurately extract Figure 10: First row: Best result obtained using the models of [11] and [23] with the deformable model initialized with a single connected component from outside.Second and third row: Results obtained using the two-phase model of [21] when the deformable model was initialized from outside, and from both inside and outside, respectively; Similar results were obtained with this initialization for the models of [11] and [23].Fourth row: Results obtained with the model of [17] with region descriptors based on the mean and variance for both the foreground and background.Fifth row: Results obtained with our active ribbon model.The columns correspond to iterations 1, 34 and 71. the outer cellular boundaries or the inner cellular boundaries alone, but none of them could segment the cellular membranes in a single isolated partition, making them invalid for the analysis and study of myelin thickness.The parameters chosen for the ribbon-consistency force field in the last row were |b| = 10px and |c| = 25px and a = 10 −9 .
Figure 11 shows a gradual comparison between a classical geometrically-unconstrained multiphase level set model and our active ribbon model.The parameters chosen for the ribbon-to-ribbon interaction force were d = 50px and e = 2px.It is important to note that in our model, by using two level set functions to represent each ribbon, we allow the ribbons to overlap on the image.This overlap is a consequence of our partition encoding of Eq. ( 2), where the zero-level sets of multiple level set functions can intersect.This encoding uses a total of 2m level set functions for m ribbons, but gives direct access to each ribbon via their zero-level sets.Such encoding also guarantees that ribbons can share cellular boundaries and therefore agrees with the biological model of neural membranes in cellular biology.Finally, such overlap facilitates the tracking of neurons throughout the volume of brain tissue, since this way each ribbon can more easily sit on image boundaries.
Figure 12 shows the tracking of a neural process, where in each section, the active ribbon model is initialized with the results obtained in the previous cross-section of the brain tissue.Such approach to tracking only works a neural process is orthogonal to the image plane, since otherwise neural processes would experiment large displacements between consecutive sections.
Finally, Fig. 13 shows an example of a failed segmentation with one ribbon, where the presence of multiple adjacent cellular membranes confused the ribbon to believe it is segmenting a single process.
The active ribbon parameters for the examples of Fig.  of appearance of neural processes were accounted by increasing the overall elasticity of the model, which can be controlled with σ 1 , σ 2 and σ I .

Conclusions and future work
This paper is, to the best of our knowledge, the first work that shows a direct way of geometrically constraining a multiphase level sets flow for image segmentation.This is done using dynamic force fields such as those introduced previously in the literature for active contours for helping them deal with local minima.Our method requires no training set and can be easily combined with other variational level set segmentation models.Future work includes extensions to 3D deformable models, and the study of the ability of dynamic force fields to induce other possible geometric arrangements.

Figure 7 :
Figure 7: (a): Active ribbon (in green) on a sample image.(b): Normalized feature maps e 1 (left) and e 2 (right) for the inner and outer side of the ribbon (blue).

Figure 8 :
Figure 8: (a): Close-up of component v 1 of Force field F 2 generated from the ribbon and sample image in Fig. 7(a).The vector field points towards the inner side of the cellular membranes (b): User defined vector kernel for vector field convolution.See[11] for more details on different types of kernels.

Figure 9 :
Figure 9: (a): A force of mutual repulsion or attraction is created between each pair of ribbons based on the result from previous sections.(b): σ I is chosen to create a elastic repulsion (region A) or attraction force (region C) to keep the distance between the ribbons within some range (region B).

Figure 11 :
Figure 11: Gradual comparison between a classical geometrically-unconstrained multiphase level set model and our active ribbon model.First row: Results obtained with the model of [21] with λ i set to 0 for the background phases in Eq. (3).Second row: Results obtained with our force field for ribbon consistency.Third row: Results obtained with all the force fields discussed in the paper enabled.The columns correspond to iterations 1, 34, 55 and 79.
10 and Fig. 11 were µ 1 = 1.5, µ 2 = µ 3 = 1 and α i = λ i = 1 for the foreground phases and 0 for the background ones.The parameters for Fig. 12 were: µ 1 = 1, µ 2 = 0.35, and same α i 's and λ i 's as in the segmentation examples.The function graphs of Figs. 9 (b) and 5 were approximated using quadratic polynomials that interpolate the point coordinates derived from the horizontal and vertical lines corresponding to the parameters a, b, c, d and e.Large variations

Figure 12 :
Figure 12: Left-to-right and top-to-bottom: Tracking of a cellular membrane that is orthogonal to the scanning plane.The active ribbon is initialized with the results obtained in the previous cross-section of the brain tissue.

Figure 13 :
Figure 13: Example of a failed segmentation with one ribbon.The images correspond to iterations 1, 10, 16 and 30.