# 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 Cunha^{†}^{†}thanks: 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.

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 .

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

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.

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