Geometric Median Shapes

Geometric Median Shapes

Abstract

We present an algorithm to compute the geometric median of shapes which is based on the extension of median to high dimensions. The median finding problem is formulated as an optimization over distances and it is solved directly using the watershed method as an optimizer. We show that computing the geometric median of shapes is robust in the presence of outliers and it is superior to the mean shape which can easily be affected by the presence of outliers. The geometric median shape thus faithfully represents the true central tendency of the data, contaminated or not. Our approach can be applied to manifold and non manifold shapes, with connected or disconnected shapes. The application of distance transforms and watershed algorithm, two well established constructs of image processing, lead to an algorithm that can be quickly implemented to generate fast solutions with linear storage requirements. We demonstrate our methods in synthetic and natural shapes and compare median and mean results under increasing contamination by strong outliers.

Geometric Median Shapes

Alexandre Cunhathanks: We are grateful for funding to the Center for Advanced Methods in Biological Image Analsysis from the Beckman Institute at Caltech, and support from the Gordon and Betty Moore Foundation through grant GBMF3406. Correspondence should be addressed to cunha@caltech.edu.
Center for Advanced Methods in Biological Image Analysis
Center for Data Driven Discovery
California Institute of Technology, Pasadena, CA, USA


Index Terms—  Shape Analysis, Geometric Median, Median Shapes, Image Segmentation, Collaborative Segmentation

1 Introduction

The computation of an average shape from a pool of candidate shapes is an important tool when classifying, comparing, and searching for objects in a database or an image. It is helpful, for example, when morphology is key in the identification of healthy and diseased cells, when differentiating wild from mutated cells, and when combining different candidate segmentations for the same image. The average usually takes the form of the mean or the median, the latter being preferred due to its robustness when data is potentially corrupted by outliers. Shape analysis is a vast field of research and many solutions have been proposed to determine mean shapes – see, for example, [1] and citations therein – with applications ranging from cell biology to space exploration and methods in the fields of elastic shape analysis [2, 3] and most recently in deep learning [4, 5].

We propose here a method that quickly build median shapes from contours representing shapes in images. Our median is computed using the notion of geometric median, an area of active research and promising results. The geometric median, also known as the spatial median, Fermat-Weber point, and median, is an extension of the median of numbers to points in higher dimensional spaces [6]. It has the property that it minimizes the sum of its Euclidean distances to points ,

(1)

and it is robust in the presence of outliers. The geometric median has a breakdown point of 50%, i.e. the data can be corrupted with up to 50% and it still gives a good prediction representing the sample [7]. Unlike the mean of points, which can be found directly using a analytical expression, the geometric median can only be approximated by numerical algorithms. Researchers have developed fast numerical solutions for the convex problem in eq.(1). Weiszfeld developed in 1937 (translated reprint [8]) at the age of 16 years old, an algorithm which is still quite competitive (see e.g. [9]). Recent developments claim faster algorithms that can handle large data sets [10, 11]. The geometric median lies within the convex hull of the given points (see e.g. [12]), a property we will explore in our approach. The median point always exist and it is unique – due to convexity of eq.(1) – provided the points are not collinear [13].

Comparatively, very few have proposed methods for computing medians of shapes. Fletcher et al. [14] developed a method to compute the geometric median on Riemannian shape manifolds showing robustness to outliers on the Kendall shape space. In [15] the authors formulate the median finding problem as a variational problem taking into account the area of the symmetric difference of shapes, a model explored much earlier in [16]. The proposed method is extremely slow, taking minutes to compute a median for a few small images. More recently, the work in [17] represents shapes as currents and consider median finding using a variational formulation which the authors claim lead to an efficient linear program in practice. They observed that smooth shapes do not necessarily lead to smooth medians, an aspect we also verified in our approach.

2 Methods

2.1 A distance for planar shapes

We define here the expression to measure the distance between planar shapes represented by their contours. It is this metric that we will be optimizing when computing the geometric median of a set of shapes. Shapes that are near (far) in space express lower (higher) distance values. We assume contours are rectifiable curves, but do not restrict their topology to planar manifolds, crossings are allowed.

Let give the distance of a point to a curve ,

(2)

where the infimum is attained for one or more points closest to , and the norm gives the Euclidean distance between points (see Fig.1). Note that . We then define the distance from a rectifiable curve , i.e. , to another curve using a line integral over all points ,

(3)

where is given by eq(2), is the arc length differential, and is the length of curve . By integrating distances for all points along the curve and dividing by its length we have a measure that tells us how far, in average, is the source curve away from the target curve . If is a parameterization of on the unit interval, for , we then have

(4)

This distance is independent of curve parameterization. It follows that and . If and are parallel curves for which , then , which is what we intuitively expect. As an example, for two concentric circles with radii , we have . Another property of is its invariance under isometries (rigid transformations including reflections), owing to the intrinsic property of curve length and to the preservation of point distances under such transformations.

Fig. 1: Illustration of distances between curves. The left figure shows the distances and from points and , respectively, on curves and . The figure on the right illustrates a case where a curve , horizontal black line, is equally distant to three different other curves, . The symmetric distance captures our two-way expectations: . Curves only differ at the bulgy part. Slanted arrows on the right figure are directed to the farthest points on from . To compute , curves are 8-connected in 429x295 images.

 

In general, this distance is asymmetric, , unless the curves themselves are symmetric with respect to an axis. Similar to the Hausdorff distance, we combine and to obtain a symmetric distance ,

(5)

possessing the following properties: {enumerate*}[label= (0) ]

(non–negativity);

(identity);

(symmetry);

(triangle inequality). Properties (1) to (3) are due to the Euclidean norm used in , while the triangle inequality property requires formal proof, a task for future work. All properties have been experimentally validated.

2.2 Geometric median of shapes

We formulate the problem of finding the geometric median of shapes as an optimization problem over distances, similar to the high dimensional geometric median of points. Find the optimal contour that minimizes the sum of distances to a given set of contours ,

(6)

where is given by eq.(3) and the solution lies in the convex hull of , . When each curve degenerates to a point, with , then the problem is reduced to finding the geometric median of points. From eq.(6) we then have

(7)

where and commute due to Tonelli’s theorem for non-negative functions [18] – we have – and is the closest point to , as in eq.(2) and Fig.1. It should be clear that the location of each depends on the position of in , thus . Let’s for a moment ignore above and consider an approximation of by an ordered set of points . At the solution , the term in parentheses on the last equation in (7) is optimized for each which implies that any point in is in fact the geometric median of points , each from a single ,

(8)

Since we know , see e.g. [12], we can then show by contradiction that , as stated above. This is a property that will help us in the optimization. Directly computing each from points is not straightforward due to the nonlinear dependency of on . It is also not clear how to properly select a set of points from the curves to produce each median point in .

By discarding the contribution of the curve length and optimizing solely based on distances, i.e. minimizing the following functional

(9)

we end up with a simpler problem while, surprisingly, still solving the original optimization problem in eq.(6). Our results consistently demonstrate we are also optimizing for the curve length albeit not explicitly using it in the optimization. When we use a squared norm in eq.(9) and elsewhere, we say is a mean shape.

2.3 Optimization with watershed transform

We look for the contour having overall the smallest sum of distances from its points to all given contours . We solve this problem as an optimal path finding problem where the watershed with markers method [19] is applied over the accumulated distance transforms of the given curves to generate the median contour . Others have applied the watershed method to solve optimization problems on either a graph or a discrete grid. For example, Cuprie et al. [20] used the Power Watershed [21] in their formulation of the surface reconstruction from point cloud problem to minimize a weighted total variation functional where weights are proportional to Euclidean distances between points.

We abuse the notation here and use to designate both the 8-connected contour and the image containing it. We represent each contour as the zero level set of a unsigned distance transform function , such that for , for , and gives the shortest distance from to curve . Thus, from eq.(2), . If we take the sum for every point in the domain of the image, the median contour is a minimal path in the image .

To create markers for the watershed, we take the union of the basins in (the inverted ) with the complement of the convex hull of contours , . The marker imposed outside the convex hull prevents the creation of false positive regions due to eventual basins placed close to the boundary of the image. Points in regional minima are those farthest from the contours, and hence serve well as markers. The watershed lines computed from these markers form . The GEMS algorithm below generates .

1:procedure GEMS()
2:     Input: (minima height), (contour images)
3:     Output: (geometric median)
4:     
5:     for  do
6:          
7:                
8:     end
9:     
10:     
11:     
12:     
13:     
Algorithm 1 GEMS – Geometric Median of Shapes
Our current implementation of GEMS is a script that uses distance transforms, height minima, minima, convex hull, and watershed implementations from the pink [22] package, with image I/O and operations complemented by the gmic [23] tool. Note that memory space is linear in image size and does not grow with the number of contours . The only parameter that needs adjustments is the height used for the detection of watershed basis.

 

When the contours are highly dissimilar in shape and location then a shape representing the central tendency of the data is either too weak or non existent. This is reflected in the value of the height parameter : the smaller the value needed to create , the less cohesive is the data. For example, we could not obtain a median for the colored outliers in Fig.4, even after taking , as they are very distinct.

When contours reach the boundaries of the image, GEMS results need to be augmented. This is, for example, the case when building a consensus out of many binary contour segmentations for tiles of an image – see Fig.5. Disconnected basins are formed at the boundary of the image and we no longer have the larger marker outside the convex hull to eliminate them. This in turn leads to the creation of extra regions not representing the actual data. We repair this over segmentation, whenever present, by removing edges from according to the following criteria: {enumerate*}[label=(0) ]

Let . Edges in not present in are removed: if the edge has no candidate from the pool of segmentations, there is no justification to keep it; this eliminates the great majority of extraneous edges in .

Edges in which have a small count (currently 20%) in the pool and an average intensity close to background are eliminated; this is done to remove edges introduced only by a few users that can be potentially annotation mistakes.

Our experiments have consistently demonstrated that the median shapes obtained using GEMS have a cost lower than any other contour in its vicinity, , and elsewhere. One might argue that the cost could be reduced, according to eq.(3), by picking a longer curve. But this is not the case. Choosing a contour longer than , a few pixels away from the optimal configuration, increases the sum of distances at a greater rate, as the accumulated point distances grows faster than the curve length . Conversely, choosing a curve with fewer points to decrease does not reduce the cost as now the curve length is smaller and although we have now fewer points each has a higher accumulated distance. Evidences suggest then that GEMS not only finds the contour with lowest accumulated distances, , but it also picks the path with the right length. This indicates that optimizes distances , and consequently . We report values for in our experiments. q

3 Experimental Results

 

     

Fig. 2: Median and mean shapes for Ginkgo biloba leaves. Sixteen leaf curves are ordered according to their distances to their geometric median. They were midvein aligned and scaled keeping their original aspect ratio and shape, all shown on the curves panel. Their respective geometric median and mean shapes are the black curves on the medians and means panels. Median and mean are pixel units distant. Note that the petiole (lower tip) is better captured by the median. The gray mean shapes were obtained after tainting data with two (black and brown curves) and five outliers, shown on the rightmost outliers panel, resulting in poor average shapes distant and from the true mean. We can no longer obtain a simple mean curve when these five outliers are considered. These outliers barely affect the median. The gray median, shown on the medians panel, is only pixel units apart from the true median. It was obtained after contaminating data with ten outliers, twice as much used for the mean, thus illustrating the superiority of the geometric median. Leaf images have 800x800 pixels.

We demonstrate GEMS on a set of planar shapes ranging from plant leaves to hand traced circles and non-manifold contours from interactive segmentation of biological cells. The examples in the figures are self contained and address particular cases of interest.

In the Ginkgo biloba leaves experiment, Fig.2, we compute median and mean shapes for sixteen aligned contours, with and without outliers. Median and mean are initially pretty close, , but their distances increase with the introduction of outliers. The influence of data alignment for the median computation is presented in Fig.3, where tulip leaves are registered in different ways leading to distinct but close in shape medians. We illustrate in the circles example, Fig4, a case that the geometric median holds the central tendency even when the contamination with outliers is beyond 50%. This is not the case for the mean shape which deteriorates quickly as we increase the number of random outliers tainting the data. In the last experiment, Fig.5, we use the augmented GEMS described in section 2.3, to create a consensus segmentation for a small set of interactively generated segmentation results. These segmentations are for small tiles cropped from large cell images. We compare our results with those created by the Staple method [24] and show that our approach is comparable, and sometimes better, and can produce results much faster.

In all these experiments, the fast generation of geometric medians has been robust to outliers and to spurious data frequently present in biological images, showing that the proposed approach is competitive and suitable for averaging large data sets, contaminated or not.

   

Fig. 3: Alignment and medians. Sixteen tulip tree (Liriodendron tulipifera) leaves, shown on the left panel, are aligned in different ways leading to different median shapes (bottom right curves) conforming to the alignment. We overlay the medians, shown as thick red contours, onto the aligned leaf curves to visually demonstrate their placement: in the convex hull of the given shapes and not biased towards the extremes – see e.g. the top right case. Registrations were: no registration (top,left), registration to a line passing through side lobe corners (top right), and rigid registration, all done using the same chosen target image from from the pool.

 

    #outliers 2 4 6 8 10 12 14 median 0.52 0.88 1.15 1.48 1.88 2.47 3.04 mean 9.91 20.47 27.72 35.13 40.01 43.70 47.69  

Fig. 4: Robustness to outliers. Thirteen noisy circles hand traced by different users, shown on top row, were collected using our Collaborative Segmentation web application. We progressively added randomly generated outliers to contaminate the circle data with up to fourteen outliers, 50%+1 contamination (colored rubber bands on outliers panel, one color per outlier, some disconnected). These are strong outliers as they greatly differ in size, shape, topology, and location when compared to the circles – see circles+outliers panel. For each additional outlier added to the pool of circles, we compute new median and mean shapes and then measure their distance to, respectively, the median and mean shapes obtained without outliers. These numbers are shown on the table above for seven out of fourteen contaminations. The geometric median is marginally affected – is at most 3 pixel units – while the mean significantly deteriorates, early and quickly. Panels medians and means show respectively geometric median and mean shapes as outliers progress, from one to fourteen. All processed images have 400x400 pixels.

 

Fig. 5: Segmentation consensus. Results of interactive segmentation (columns second to sixth) are shown for four small tiles (first column) of much larger images. Although users create different segmentations for the same image, they nevertheless exhibit some agreement most prominent in clear, unambiguous parts of the image. Using our geometric median to combine users’ results into a single one eliminates outliers and spurious annotations, e.g. dangling scribbles and speckles, as shown in the last column. Our solution is on par or better (third row) than results obtained with the Staple method [24], shown on the panel before last, and orders of magnitude faster (e.g. 1.06s against 237.35s for sixteen 800x800 images) when compared to the expectation maximization approach adopted by Staple. Users’ segmentation results are from our Collaborative Segmentation web application. .

 

4 Conclusions

We presented an algorithm, GEMS, to compute the geometric median of shapes which is simple and straightforward to implement and requires adjustment of a single parameter. The approach leads to median shapes that are robust to outliers supporting a contamination of over 50% while still capturing the central tendency of given shapes. We adopted the watershed transform as the optimizer to identify paths of minimal cost in the accumulated distance transform image. Our experiments have shown that the proposed symmetric distance is optimized for the geometric median contour even though the curve length is not explicitly used in the optimization. From our initial experimentation, our approach compares well with the Staple method in the segmentation consensus problem while being faster by many orders of magnitude.

References

  • [1] Luciano da Fontoura Costa and Roberto Marcond Cesar Jr, Shape analysis and classification: theory and practice, CRC press, 2010.
  • [2] Washington Mio, Anuj Srivastava, and Shantanu Joshi, “On shape of plane elastic curves,” International Journal of Computer Vision, vol. 73, no. 3, pp. 307–324, 2007.
  • [3] Anuj Srivastava, Eric Klassen, Shantanu H Joshi, and Ian H Jermyn, “Shape analysis of elastic curves in euclidean spaces,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 33, no. 7, pp. 1415–1428, 2011.
  • [4] Davide Boscaini, Jonathan Masci, Emanuele Rodolà, and Michael Bronstein, “Learning shape correspondence with anisotropic convolutional neural networks,” in Advances in Neural Information Processing Systems, 2016, pp. 3189–3197.
  • [5] Ayan Sinha, Jing Bai, and Karthik Ramani, “Deep learning 3D shape surfaces using geometry images,” in European Conference on Computer Vision. Springer, 2016, pp. 223–240.
  • [6] JBS Haldane, “Note on the median of a multivariate distribution,” Biometrika, vol. 35, no. 3-4, pp. 414–417, 1948.
  • [7] Daniel Gervini, “Robust functional estimation using the median and spherical principal components,” Biometrika, vol. 95, no. 3, pp. 587–600, 2008.
  • [8] Endre Weiszfeld and Frank Plastria, “On the point for which the sum of the distances to n given points is minimum,” Annals of Operations Research, vol. 167, no. 1, pp. 7–41, 2009.
  • [9] Amir Beck and Shoham Sabach, “Weiszfeld’s method: old and new results,” Journal of Optimization Theory and Applications, vol. 164, no. 1, pp. 1–40, 2015.
  • [10] Yehuda Vardi and Cun-Hui Zhang, “The multivariate L1-median and associated data depth,” Proceedings of the National Academy of Sciences, vol. 97, no. 4, pp. 1423–1426, 2000.
  • [11] Hervé Cardot, Peggy Cénac, Pierre-André Zitt, et al., “Efficient and fast estimation of the geometric median in Hilbert spaces with an averaged stochastic gradient algorithm,” Bernoulli, vol. 19, no. 1, pp. 18–43, 2013.
  • [12] Stanislav Minsker, “Geometric median and robust estimation in Banach spaces,” Bernoulli, vol. 21, no. 4, pp. 2308–2335, 2015.
  • [13] J. H. B. Kemperman, “The median of a finite measure on a Banach space,” Statistical Data Analysis Based on the L 1 -norm and Related Methods (Neuchtel, 1987), pp. 217–230, 1987.
  • [14] P Thomas Fletcher, Suresh Venkatasubramanian, and Sarang Joshi, “Robust statistics on riemannian manifolds via the geometric median,” in Computer Vision and Pattern Recognition, 2008. CVPR 2008. IEEE, 2008, pp. 1–8.
  • [15] Benjamin Berkels, Gina Linkmann, and Martin Rumpf, “An SL (2) invariant shape median,” Journal of Mathematical Imaging and Vision, vol. 37, no. 2, pp. 85–97, 2010.
  • [16] Alexandre Cunha, A fully eulerian method for shape optimization, with application to Navier-Stokes flows, chapter Shape Matching, pp. 9–43, PhD thesis, Carnegie Mellon University, September 2004.
  • [17] Yunfeng Hu, Matthew Hudelson, Bala Krishnamoorthy, Altansuren Tumurbaatar, and Kevin R Vixie, “Median shapes,” arXiv preprint arXiv:1802.04968, 2018.
  • [18] Gerald B Folland, Real analysis: modern techniques and their applications, John Wiley & Sons, 2013.
  • [19] Serge Beucher and Fernand Meyer, “The morphological approach to segmentation: the watershed transformation,” Optical Engineering-New York-Marcel Dekker Incorporated-, vol. 34, pp. 433–433, 1992.
  • [20] Camille Couprie, Xavier Bresson, Laurent Najman, Hugues Talbot, and Leo Grady, “Surface reconstruction using power watershed,” in International Symposium on Mathematical Morphology and Its Applications to Signal and Image Processing. Springer, 2011, pp. 381–392.
  • [21] Camille Couprie, Leo Grady, Laurent Najman, and Hugues Talbot, “Power watershed: a unifying graph-based optimization framework,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 33, no. 7, pp. 1384–1399, 2011.
  • [22] Michel Couprie, Laszlo Marak, Hugues Talbot, et al., “The Pink Image Processing Library,” https://perso.esiee.fr/~coupriem/Pink/doc/html/.
  • [23] David Tschumperlé, “G’MIC, GREYC’s Magic for Image Computing,” https://gmic.eu.
  • [24] Simon K Warfield, Kelly H Zou, and William M Wells, “Simultaneous truth and performance level estimation (STAPLE): an algorithm for the validation of image segmentation,” IEEE Transactions on Medical Imaging, vol. 23, no. 7, pp. 903–921, 2004.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
313642
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description