Active Contour Models for Manifold Valued Image Segmentation
Abstract
Image segmentation is the process of partitioning an image into different regions or groups based on some characteristics like color, texture, motion or shape etc. Active contours are a popular variational method for object segmentation in images, in which the user initializes a contour which evolves in order to optimize an objective function designed such that the desired object boundary is the optimal solution.
Recently, imaging modalities that produce Manifold valued images have come up, for example, DTMRI images, vector fields. The traditional active contour model does not work on such images.
In this paper, we generalize the active contour model to work on Manifold valued images. As expected, our algorithm detects regions with similar Manifold values in the image. Our algorithm also produces expected results on usual grayscale images, since these are nothing but trivial examples of Manifold valued images. As another application of our general active contour model, we perform texture segmentation on grayscale images by first creating an appropriate Manifold valued image. We demonstrate segmentation results for manifold valued images and texture images.
1 Introduction
Image segmentation approaches are based on a characteristic property that defines the region of interest to be segmented out of the image, for example color, texture, motion, shape and/or others.
Image segmentation approaches can broadly be categorized into edge based or region based ones. There are various approaches in both categories based on intensity, color, texture and motion using statistical and geometrical framework [5].
Active contour model is a popular segmentation approach, which has both edge based [3] and region based [4] versions.
Active contours, also known as ‘snakes’, are based on evolving an initial contour towards the boundary of an object to be detected [15]. Usually, this evolution equation is a gradient descent for minimizing an appropriate energy functional.
In the Geodesic active contours model [3], the energy functional is written as a length functional with a modified metric such that the minimum, i.e. curve of minimum length corresponds to the object’s boundary. A region based approach was proposed by Chan and Vese [4] where the energy functional was not based on edges. While the traditional approach works with parametric representation of curve, a simpler and efficient representation the level set approach was introduced by Osher and Sethian [20].
Grayscale images can be modeled as functions where is the image domain. These days, there are imaging modalities and other data that can be modeled as a function , where is a Riemannian manifold. We call such images Manifold valued images (henceforth written as MVI). Examples are DTMRI images where , the set of symmetric positive definite matrices, optical flow field images or wind velocity data where . The standard active contour methods will not work on MVIs. In this paper we generalize the active contour model to work on MVI. We adapt both, Geodesic active contours and ChanVese active contour models to work on MVIs. Since is also an example of a manifold, our model also works on grayscale images (modeled as ).
As another application, we pose the texture segmentation problem as a MVI segmentation problem, as follows. Since covariance matrices are often used to characterize textures, we define a covariance matrix of grayvalued images in a neighborhood on every pixel of the given texture image. Covariance matrices are symmetric positive definite (actually semidefinite) matrices of size , the set of these matrices forms a manifold, denoted by . Our algorithm is then able to identify regions of similar covariance matrices, thus segmenting textures in the image. In the rest of the paper denotes the usual Euclidean norm, denoted absolute value of a real number and symbols for other norms used have been defined where used for the first time.
1.1 Related Work
There have been related works on DTMRI image segmentation using active contours, for example work by Lisa Jonasson [13],[14], where the surface is evolved according to
(1) 
where is a speed term proportional to the similarity of the Diffusion tensors of adjacent pixels along the normal direction to the surface, is a curvature based regularization term and is the unit normal to the surface being evolved. The similarity measure between two tensors and used in the above mentioned paper is the Normalized Tensor Scalar Product (NTSP):
A more ChanVese type active contour model was proposed by Wang and Vemuri [29] for DTMRI segmentation, where they define the following energy function on the set of curves :
(2) 
where is the tensor defined at , () is the interior (exterior) of the curve , and
(3)  
(4) 
are the mean tensors in the interior and exterior of the curve . They use the Frobenius norm () in their computations.
In [17], the authors use Geodesic active contour to segment DTMRI images. They use the Euclidean, KLdivergence based and the geodesic distance metric to induce an edgestopping term. They conclude that the geodesic distance based metric yields the best segmentation results.
In [28], the authors use a front propagation scheme to segment DTMRI images.
Sagiv, Sochen and Zeevi[25] propose a ChanVese type active contour model for texture segmentation. Instead of defining the model on intensityvalued texture image, they first compute responses to a manually selected set of Gabor filters :
The scale, orientation and the response which are maximal are included into the feature at every point of the image. This feature space is a D manifold embedded in (they use ). ChanVese active contour is then run on an image where every pixel maps to a point in this feature space. These features can also be used to induce a metric on the D manifold. Lee et.al. [16] use this metric to define a stopping term for Geodesic active contours that segments out texture.
One of the other commonly used feature for characterizing textures is the covariance matrix. In [31], the authors use covariance matrices of the intensity and first order Gaussian derivatives at every pixel over a neighborhood as a feature and segment textures using multiscale Graph cuts. Donoser and Bischof [8] use covariance matrices of intensity and first and second derivatives of intensity values as features to characterize textures. They use the manifold distance on the set of positive semidefinite matrices to extend the ROISEG [7] clustering algorithm to segment textures.
Our contribution is twofold: We give a general active contour model that can be used with intensity valued, vector valued, manifold valued images and as a computer vision application we show how this algorithm can be used to segment textures.
The paper is organized as follows.
In the next section we provide a brief background on active contours, manifold valued images and cite other work on manifold valued data. In section 3, we explain our proposed active contour model for MVI segmentation, followed by implementation details and results.
2 Background
2.1 Active contours and Level sets
In classical active contours [15], the user initializes a curve on an intensity image which evolves and stabilizes on the object boundary. The gradient descent of an energy functional, , given by
(5)  
(6) 
where and are real positive constants, and are first and second derivatives of and is the image gradient, gives us a curve evolution equation. The first two terms are regularizers, while the third term pushes the curve towards the object boundary.
Geodesic active contours [3] are an active contour model where the objective function can be interpreted as the length of a curve in a Riemannian space with metric induced by image intensity. The energy functional for geodesic active contour is given by
where is a positive monotonically decreasing edge detector function. One such choice is . We set to make notations simpler. The curve evolution equation that minimizes this energy is given by
(7) 
where is the inward unit normal and is the curvature of the curve .
A convenient computational procedure for curve evolution is the level set formulation [20, 18]. Here the curve is embedded in the zero level set of a function which evolves so that the corresponding zero level set evolves according to the desired curve evolution equation.
For a curve evolution equation of the form , the corresponding level set evolution is . See Appendix in [3]. In particular the level set evolution for geodesic active contours is given by
(8) 
Another active contour approach was introduced by Chan and Vese [4] where the energy function was based on regional similarity properties of an object, rather than its edges (image gradient). Suppose that is the initial curve defined on the domain of the intensity image . can be divided into two parts, interior (denoted by ) and exterior (denoted by ). Let us represent the mean gray value of the region and by and respectively, then the energy function for which the object boundary is a minimum is given by
(9)  
(10) 
After adding some regularizing terms the energy functional , is given by
(11)  
(12)  
(13) 
where are fixed scalar parameters. The level set evolution equation is given by
(14) 
where is a smooth approximation of the Dirac delta function, and the level set function is assumed to be negative in the interior region and positive in the exterior region of the curve . A nice survey on active contours and level set implementation can be found in [1]. In this paper, we generalize the two active contour models, the Geodesic active contour and ChanVese active contours model to segment objects in MVIs. Before describing our active contour model, we will now introduce and define the kind of images we focus on in this paper.
2.2 Manifold Valued images(MVI)
In this paper, instead of working with intensity images represented as functions , we work on images represented as , where is a Riemannian manifold^{3}^{3}3Manifold is a topological space which is locally Euclidean and a smooth Manifold equipped with a smooth inner product (called Riemannian inner product) on tangent space at every point is called a Riemannian manifold. In what follows we assume the familiarity with basic concepts from differential geometry like geodesics, Exp map, Log map etc. For the sake of completeness, we define these terms in the Appendix. A thorough introduction and discussion can be found in the textbooks [2], [6].
We present some practical instances of MVI on which we propose to segment ”objects” using active contours. Refer to Figure 1 to see examples of the following MVIs.

Diffusion tensor magnetic resonance imaging (DTMRI) produces diffusion tensors corresponding to diffusion of white matter. The diffusion tensors produced are found to be symmetric and positive definite matrices of size , which forms a Riemannian manifold. A detailed analysis of Diffusion tensor data from DTMRI can be found in the work by Fletcher and Joshi [9]. In order to display such images in this paper, we represent the positive definite matrix by a planar projection of an ellipsoid with major axis given by its eigenvectors and length of the axis being proportional to the corresponding eigenvalues.

A variety of vector field images are produced in many application like directional field of wind flow, optical flow [12], and others. The manifold in this case could be one of the following: , , , or .

Given a texture intensity image , we compute a covariance matrix at every point in the following manner. Let denote the vector containing intensity values of the neighborhood of the pixel , arranged in a specific order. Let be the covariance matrix of patches in some neighborhood defined by , where denotes transpose of . Covariance matrices are symmetric positive definite^{4}^{4}4actually positive semidefinite, but we assume independant random variables giving us positive definiteness. matrices which form a manifold, denoted by . Thus we get a manifold valued map .
In the last decade or so, a lot of work has been done on manifold data. We mention a few here. Computing nonlinear statistics on such data was proposed by Fletcher and Joshi [11] and by Xavier Pennec [22]. Weickert and Brox [30], Tschumperlè and Deriche[26] and Rosman et al. [24] have all proposed regularization schemes for vector valued and/or matrix valued images using PDE’s. Tuzel et al. [27] presented methods using Lie group modeling for tracking objects.
In the next section, we describe our active contour model for MVI segmentation.
3 Adapting active contour model for MVI segmentation
We generalize the Geodesic active contours and ChanVese active contours for MVIs.
3.1 Geodesic active contours(GAC)
In edge dependent active contour models, for example GAC, the curve evolution is made to stop at the object boundary by defining a speed function (or a metric) that is inversely proportional to an edge detection function which is itself a function of the image gradient (Refer to Section 2). The gradient of a differentiable intensity function is a vector along which the function increases the most and whose length is that amount of increase. With Euclidean inner product on , the gradient at a point can be computed as . Note that in geodesic active contours, the image gradient magnitude (not the image gradient vector itself) plays a significant role: the edge detector function is a function of the image gradient magnitude. Image gradient magnitude can be interpreted as the maximum rate of change in the value of the function at a point .
(15) 
where is the differential of the function at point , and are evaluated at . Since this map is linear, the magnitude of maximum increase is the same as the magnitude of maximum decrease in the function value. The function increases the maximum along and decreases the maximum along while the magnitude of rate of change in both cases is equal to . We generalize this concept for active contours to work on MVIs. Note that MVI are maps of the kind . It does not make sense to talk about gradient of such maps, but we define an analogous concept to the gradient magnitude of intensity images for MVIs. The differential of the MVI function^{5}^{5}5In what follows, we assume that is differentiable in the sense described in the book[2]. The definition of Differential of a smooth map between manifolds is given in the Appendix 6. at point is given by
(16)  
(17) 
where is the tangent space to at the point , and . Moreover
(18) 
where is the Riemannian norm on the tangent space . We define the gradient magnitude of at (denoted by ) as the maximum possible value of the norm of the differential:
(19) 
Let the inner product on be defined as
(20) 
where is a symmetric positive definite matrix that varies smoothly over . Let
(21) 
Then it can be shown that
(22) 
where is the maximum singular value of the matrix . Finally we explain how we compute and on discretized MVIs. In case of intensity images, the following approximation is frequently used:
(23)  
(24) 
Subtraction in vector space can be reinterpreted on Riemannian manifolds as the Riemannian Log map[21]. Using this, the analogous approximations to for MVIs , are given by
(25)  
(26) 
where is the Riemannian Log map at point for point . Having defined the required gradient magnitude for MVIs, the level set evolution equation corresponding to geodesic active contours is given as
(27) 
where
(28) 
3.2 ChanVese active contour model
In ChanVese active contour model the energy function is
(29)  
(30)  
(31) 
where and are mean over the interior and exterior of the contour , respectively.
For MVIs, are the intrinsic means of the manifold data in the interior and exterior of the curve . Intrinsic mean of a collection of points is the minimizer of the sum of squared Riemannian distances from each of the given points:
(32) 
where is the Riemannian distance on and is a generalization of the Euclidean distance to Riemannian manifolds [21]. This is computed using a gradient descent approach. For a detailed explanation refer to [11], and refer to Section 4 for a computational algorithm to compute the intrinsic mean. The corresponding energy function for ChanVese active contours on MVI is thus given by
(33) 
where and are intrinsic means of manifold points lying in the regions and respectively. The level set evolution equation comes out to be
(34) 
where, following the ChanVese model, we assume that the level set function is negative in the interior and positive in the exterior of the curve . In the next section, we implement both GAC and ChanVese active contour model on several examples of MVIs. We also provide computational algorithms for computing Riemannian Exp map, Log map and intrinsic mean on the corresponding manifolds. The Riemannian Exp map is the inverse of Riemannian Log map, and intuitively can be interpreted as vector addition in Riemannian manifolds [21].
4 Experiments
Our Active contour model requires computing the Exp map, Log map and intrinsic mean for each manifold. For short definitions of these terms refer to the Appendix 6. A summary of algorithms to compute these maps on (D sphere), (D rotation matrices) and manifolds are provided next, while detailed derivations can be found in one or more of these papers: [10, 11, 23, 19].
4.1 Computation on different manifolds


Exp map
(35) where and .

Log map
(36) where .



Exp map
(37) (38) where and is a representation of in .

Log map
(39) (40) where , and .



Exp map
Input : Initial point
Tangent vector
Output :
Let diagonal)
Let diagonal)

Log map
Input : Initial point
End point
Output :
Let diagonal)
Let diagonal)

We have used the following algorithm to find the intrinsic mean on all manifolds, details of which can be found in [11].
With all the required computational machinery set, we next demonstrate our segmentation results for both Geodesic active contours and ChanVese active contours on various manifold valued images.
4.2 Results
Results on , , and manifold valued images for Geodesic active contours as well as ChanVese active contours are shown in Figure 2. Result on Real DTMRI data^{6}^{6}6from http://www.cise.ufl.edu/~abarmpou/lab/fanDTasia/ is shown in Figure 3.
Texture segmentation problem can be posed as an MVI segmentation problem, as already explained in Section 2. We obtain a valued image from a texture image over which our algorithm successfully segments different textures^{7}^{7}7from http://sipi.usc.edu and http://www.nada.kth.se/cvap/databases/kthtips as shown in Figure 4. Since texture boundaries are not defined based on grayvalue image edges, the geodesic active contour model does not yield appropriate segmentation results. The results shown are for ChanVese active contour model. We have used a covariance matrix to characterize a texture and it is computed over a neighborhood of pixels.
5 Conclusion and future scope
In this paper we have generalized the active contours for MVI segmentation. We provide several such examples which can be dealt under our general framework. The drawbacks and benefits of the proposed model are same as those of the usual active contour models. As a computer vision application we pose the texture segmentation problem as an MVI segmentation problem and demonstrate some texture segmentation results using our algorithm. We take a neighborhood of size pixels in an intensity image to form covariance matrices over a larger neighborhood of size at every pixel to get a valued image . This fixes the scale of features in which we are interested. Some textures may need multi scale information to be properly characterized. Simply increasing the neighborhood size is not going to help in dealing with large scale texture as it may prevent proper localization of the texture boundary. One needs to incorporate a mechanism that automatically detects the scale of the given texture.
The computational cost for the model extended from ChanVese active contour is high since one needs to compute the Riemannian Exp and Log map for every manifold value in the image. This cost will of course vary from manifold to manifold. The texture segmentation algorithm is computationally very expensive: To generate results on images of size pixels, it took about minutes on a GHz Intel laptop. This is primarily due to repeated Exp and Log map computation on manifold. One may use multigrid approaches to speed up the evolution.
References
 [1] Gilles Aubert and Pierre Kornprobst. Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations (Applied Mathematical Sciences). SpringerVerlag New York, Inc., Secaucus, NJ, USA, 2006.
 [2] W.M. Boothby. An introduction to differentiable manifolds and Riemannian geometry. Pure and Applied Mathematics. Elsevier Science, 1975.
 [3] Vicent Caselles, Ron Kimmel, and Guillermo Sapiro. Geodesic active contours. International Journal of Computer Vision, 22(1):61–79, feb 1997.
 [4] Tony F. Chan, B. Yezrielev Sandberg, and Luminita A. Vese. Active contours without edges for vectorvalued images. Journal of Visual Communication and Image Representation, 11:130–141, 2000.
 [5] Daniel Cremers, Mikael Rousson, and Rachid Deriche. A review of statistical approaches to level set segmentation: Integrating color, texture, motion and shape. International Journal of Computer Vision, 72:215, 2007.
 [6] M.P. do Carmo. Differential geometry of curves and surfaces. PrenticeHall, 1976.
 [7] Michael Donoser and Horst Bischof. Roiseg: Unsupervised color segmentation by combining differently focused sub results. In CVPR, 2007.
 [8] Michael Donoser and Horst Bischof. Using covariance matrices for unsupervised texture segmentation. In ICPR, pages 1–4, 2008.
 [9] P. Thomas Fletcher and Sarang Joshi. Riemannian geometry for the statistical analysis of diffusion tensor data. Signal Process., 87(2):250–262, feb 2007.
 [10] P. Thomas Fletcher and Sarang C. Joshi. Principal geodesic analysis on symmetric spaces: Statistics of diffusion tensors. In ECCV Workshops CVAMIA and MMBIA, pages 87–98, 2004.
 [11] P. Thomas Fletcher, Conglin Lu, Stephen M. Pizer, and Sarang Joshi. Principal geodesic analysis for the study of nonlinear statistics of shape. IEEE Transactions on Medical Imaging, 23(8):995–1005, 2004.
 [12] Berthold K. P. Horn and Brian G. Schunck. Determining optical flow. Artificial Intelligence, 17:185–203, 1981.
 [13] Lisa Jonasson, Xavier Bresson, Patric Hagmann, Olivier Cuisenaire, Reto Meuli, and JeanPhilippe Thiran. White matter fiber tract segmentation in dtmri using geometric flows. Medical Image Analysis, 9(3):223–236, 2005.
 [14] Lisa Jonasson, Patric Hagmann, Claudio Pollo, Xavier Bresson, Cecilia Richero Wilson, Reto Meuli, and JeanPhilippe Thiran. A level set method for segmentation of the thalamus and its nuclei in dtmri. Signal Processing, 87(2):309–321, 2007.
 [15] Michael Kass, Andrew Witkin, and Demetri Terzopoulos. Snakes: Active contour models. International Journal of Computer Vision, 1(4):321–331, 1988.
 [16] SangMook Lee, A. Lynn Abbott, Neil A. Clark, and Philip A. Araman. Active contours on statistical manifolds and texture segmentation. In ICIP (3), pages 828–831, 2005.
 [17] Christophe Lenglet, Mikaël Rousson, and Rachid Deriche. Dti segmentation by statistical surface evolution. IEEE Trans. Med. Imaging, 25(6):685–700, 2006.
 [18] Ravikanth Malladi, James A. Sethian, and Baba C. Vemuri. Shape modeling with front propagation: A level set approach. IEEE Transactions on Pattern Analysis and Machine Intelligence, 17:158–175, 1995.
 [19] Maher Moakher and Mourad Zéraï. The riemannian geometry of the space of positivedefinite matrices and its application to the regularization of positivedefinite matrixvalued data. Journal of Mathematical Imaging and Vision, 40(2):171–187, 2011.
 [20] Stanley Osher and James A. Sethian. Fronts propagating with curvature dependent speed: Algorithms based on hamiltonjacobi formulations. Journal of computational Physics, 79(1):12–49, 1988.
 [21] X. Pennec. Statistical computing on manifolds: From riemannian geometry to computational anatomy. In Frank Nielsen, editor, Emerging Trends in Visual Computing, volume 5416 of LNCS, pages 347–386, 2008.
 [22] Xavier Pennec. Intrinsic statistics on Riemannian manifolds: Basic tools for geometric measurements. Journal of Mathematical Imaging and Vision, 25(1):127–154, July 2006. A preliminary appeared as INRIA RR5093, January 2004.
 [23] Quentin Rentmeesters and P.A. Absil. Algorithm comparison for karcher mean computation of rotation matrices and diffusion tensors. In Proceedings of the 19th European Signal Processing Conference (EUSIPCO 2011), pages 2229–2233. EURASIP, 2011.
 [24] Guy Rosman, Yu Wang, XueCheng Tai, Ron Kimmel, and Alfred M. Bruckstein. Fast regularization of matrixvalued images. In ECCV (3), pages 173–186, 2012.
 [25] Chen Sagiv, Nir A. Sochen, and Yehoshua Y. Zeevi. Integrated active contours for texture segmentation. IEEE Transactions on Image Processing, 15(6):1633–1646, 2006.
 [26] David Tschumperlé and R. Deriche. Vectorvalued image regularization with pdes: A common framework for different applications. In IEEE Transactions on Pattern Analysis and Machine Intelligence, pages 506–517, 2003.
 [27] O. Tuzel, F. Porikli, and P. Meer. Learning on lie groups for invariant detection and tracking. In Computer Vision and Pattern Recognition, 2008. CVPR 2008. IEEE Conference on, pages 1–8, 2008.
 [28] Guo W., Chen Y., and Zeng Q. Dti segmentation using an information theoretic tensor dissimilarity measure. Philos Trans A Math Phys Eng Sci., 366(10):2279–2292, 2008.
 [29] Zhizhou Wang and Baba C. Vemuri. Dti segmentation using an information theoretic tensor dissimilarity measure. IEEE Trans. Med. Imaging, 24(10):1267–1277, 2005.
 [30] Joachim Weickert and Thomas Brox. Diffusion and regularization of vector and matrixvalued images, 2002.
 [31] Horst Wildenauer, Branislav Mičušík, and Markus Vincze. Efficient texture representation using multiscale regions. In Proceedings of the 8th Asian conference on Computer vision  Volume Part I, ACCV’07, pages 65–74, Berlin, Heidelberg, 2007. SpringerVerlag.
6 Appendix
We give some basic definitions from Differential geometry required for our paper. For a thorough explanation, we refer the reader to the books [2], [6].

Differentiable Manifolds:
A differentiable manifold of dimension is a set and a family of injective mappings = of open sets of into such that
, i.e. the open sets cover .

for any pair with , the mapping is differentiable.

The family is maximal, which means that if , is such that: for each element of , with implies that is a diffeomorphism, then in fact .


Differential of a smooth map between differential manifolds:
Let be a smooth map between two differentiable manifolds and . Given a point , the differential of at is a linear mapfrom the tangent space of at to the tangent space of at .

Riemannian Metric:
A Riemannian metric on a manifold is a correspondence which associates to each point an inner product on the tangent space , which varies smoothly. In terms of local coordinates, the metric at each point is given by a matrix, , where are tangent vectors to at , and it varies smoothly with . A Geodesic curve is a local minimizer of arclength computed with a Riemannian metric. 
Geodesics:
A parameterized curve is a geodesic if , where is called the Covariant derivative and intuitively is the orthogonal projection of the usual derivative to the tangent space. Geodesics are also local minimizers of arclength. 
Exponential map:
The exponential map is a map , that maps for , to a point on obtained by going out the length equal to , starting from , along a geodesic which passes through with velocity equal to . The geodesic starting at with initial velocity can thus be parametrized as 
Log map: For in a sufficiently small neighborhood of , the length minimizing curve joining and is unique as well. Given and , the direction in which to travel geodesically from in order to reach is given by the result of the logarithm map . We get the corresponding geodesic as the curve . In other words, is the inverse of in the neighborhood.