Geodesic Centroidal Voronoi Tessellations: Theories, Algorithms and Applications
GG EODESIC C ENTROIDAL V ORONOI T ESSELLATIONS :T HEORIES , A
LGORITHMS AND A PPLICATIONS
A P
REPRINT
Zipeng Ye, Ran Yi, Minjing Yu, Yong-Jin Liu ∗ Department of Computer Science and Technology,Tsinghua University, China
Ying He
School of Computer Engineering,Nanyang Technological University, SingaporeJuly 2, 2019 A BSTRACT
Nowadays, big data of digital media (including images, videos and 3D graphical models) arefrequently modeled as low-dimensional manifold meshes embedded in a high-dimensional featurespace. In this paper, we summarized our recent work on geodesic centroidal Voronoi tessellations(GCVTs), which are intrinsic geometric structures on manifold meshes. We show that GCVT canfind a widely range of interesting applications in computer vision and graphics, due to the efficiencyof search, location and indexing inherent in these intrinsic geometric structures. Then we present thechallenging issues of how to build the combinatorial structures of GCVTs and establish their timeand space complexities, including both theoretical and algorithmic results. K eywords Voronoi tessellation · Geodesic · Computational Geometry
Let ( X, d ) be a metric space, where X is a point set and d : X × X → R is a metric. Given an open subset Ω ⊆ X , aset { V i } ki =1 is called a tessellation of Ω if V i ∩ V j = ∅ for i (cid:54) = j and ∪ ki =1 V i = Ω , where A is the closure of A . Given aset of points { g i } ki =1 in Ω , the Voronoi cell corresponding to the point g i is defined as ˆ V i = { x ∈ Ω (cid:12)(cid:12) d ( x, g i ) < d ( x, g j ) for j = 1 , . . . , k, j (cid:54) = i } . (1)Elements of { g i } ki =1 are called generators .Since 1644 (Part III of Principia Phiolosophiae, written by Descartes), Voronoi tessellations had been well studiedin the Euclidean space R n , n ∈ Z + [1], in which Voronoi tessellations (as well as their dual structures, well knownas Delaunay triangulations) had played a central role as fundamental geometric structures. Voronoi tessellations hadalso been studied in spaces with non-Euclidean metrics, including spheres [2], hyperbolic spaces [3] and the generalRiemannian manifolds [4]. Recently, due to the flourishes of big media data from digital sampling, more and more dataare appearing in the form of manifold meshes [5]. Quite different from parametrized 2-manifolds or general Riemannianmanifolds that are generally C ∞ smooth (or at least C smooth) [6], manifold meshes are only C . In our study, weadopt the discrete geodesic metric [7].In this paper, we summarize our recent work on geodesic centroidal Voronoi tessellations (GCVT) — which areprovable uniform tessellations on manifold meshes — and we show that they can be used to generate uniform remeshingin computer graphics and build content-sensitive superpixels/supervoxels for images and video in computer visionapplications.Before we introduce GCVTs, we present two close concepts of CVT and RCVT in related work. ∗ Corresponding author a r X i v : . [ c s . G R ] J u l PREPRINT - J
ULY
2, 2019
CVT had been well studied in science and engineering, with a wide range of applications including data compression indigital image processing, optimal quadrature in numerical methods, quantization and clustering in machine learning,finite difference methods in solid mechanics and fluid dynamics, distribution of resources in operational research,cellular patterns in biology, and the territorial behavior of animals; see [8] for an excellent survey.Let V be a finite region in R n . The mass centroid m of V is defined by m (cid:44) (cid:82) x ∈ V xρ ( x ) dx (cid:82) x ∈ V ρ ( x ) dx , (2)where ρ is a density function defined in V . Given k points { g i } ki =1 in a domain Ω ⊆ R n , we can define the Voronoiregion ˆ V i corresponding to each g i based on the Euclidean metric d E and the Voronoi tessellation { ˆ V i } ki =1 of Ω . Let m i be the mass centroid of each Voronoi region ˆ V i . A Voronoi tessellation is called CVT if all the generators are masscentroids, i.e. g i = m i , i = 1 , . . . , k. (3)Arbitrarily chosen generators are usually not the mass centroids of their associated Voronoi regions so an arbitraryVoronoi tessellation cannot be a CVT. It can be shown that CVT minimizes the following energy (a.k.a. CVT energyfunctional): ε E ( { ( p i , V i ) } ki =1 ) = k (cid:88) i =1 (cid:90) x ∈ V i ρ ( x ) d E ( x, p i ) dx, (4)where ρ is a density function defined in Ω , { V i } ki =1 is an arbitrary tessellation and { p i } ki =1 is any set of k points in Ω .To compute CVT, various local methods had been surveyed in [8]. In particular, the Lloyd method [9] is a simple yeteffective local method that iteratively computes mass centroids and Voronoi tessellation. To extend the domain partitioning of CVT from Euclidean space to general spaces, Du et al. [10] proposed theconstrained CVT (CCVT) that works on a compact and continuous surface S ⊂ R N defined by S (cid:44) { x ∈ R N : g ( x ) = 0 and g i ( x ) ≤ , for i = 1 , , · · · , m } (5)where { g i } mi =0 are some continuous functions.Given a set of points { g i } ki =1 ∈ S , the constrained Voronoi region ˆ V i corresponding to each g i based on the Euclideanmetric d E and restricted in S is defined as ˆ V i (cid:44) { x ∈ S (cid:12)(cid:12) d E ( x, g i ) < d E ( x, g j ) for j = 1 , . . . , k, j (cid:54) = i } . (6)The constrained mass centroid m i of ˆ V i on S is defined to be the solution of the following problem: min x ∈ S ε i ( z ) , where ε i ( z ) = (cid:90) x ∈ ˆ V i ρ ( y ) d E ( x, z ) dx (7)A tessellation on S is called CCVT if and only if the set { g i } ki =1 ∈ S are both the generators and constrained masscentroids of the tessellation { ˆ V i } ki =1 .The applications of CCVT including polynomial interpolation and numerical integration on the sphere are illustrated in[10]. Given a k-dimensional compact differentiable manifold M , a set { V i } ki =1 is called a tessellation of M if V i ∩ V j = ∅ for i (cid:54) = j and ∪ ki =1 V i = M . Denote by d g ( a, b ) the geodesic distance between a and b in M . Given a set of points { g i } ki =1 in M , the geodesic Voronoi region corresponding to the point g i is defined as ˆ V i = { x ∈ M (cid:12)(cid:12) d g ( x, g i ) < d g ( x, g j ) for j = 1 , . . . , k, j (cid:54) = i } . (8)2 PREPRINT - J
ULY
2, 2019 (a) A 2-manifold M (b) A constrained Voronoi tessellation (c) A geodesic Voronoi tessellation V ( g k ) Figure 1: (b) shows a constrained Voronoi tessellation of 30 point generators on the 2-manifold (a) and its Voronoi cell V ( g k ) is disconnected; i.e., it consists of three disjoint components (shown in red areas). Given the same set of pointgenerators, each cell in geodesic Voronoi tessellation (c) is guaranteed to be connected.Figure 2: A geodesic Voronoi tessellation on a genus-2 manifold whose cells are multiply connected.Elements of { g i } ki =1 are called generators . The set of geodesic Voronoi regions { ˆ V i } ki =1 is called geodesic Voronoitessellation of M . A geodesic Voronoi region is a connected domain and is a non-empty compact set [11]. For eachgeodesic Voronoi region ˆ V i , the nominal mass centroid m i of ˆ V i on M is defined to be the solution of the followingproblem min z ∈ M ε i ( z ) , where ε i ( z ) = (cid:90) x ∈ ˆ V i ρ ( x ) d g ( x, z ) dx. (9)It can be shown that ε i ( z ) is continuous and the domain M is compact so ε i ( z ) has at least one global minimum in M .Therefore, the nominal mass centroid exists.Given k points { g i } ki =1 , we can define the geodesic Voronoi region ˆ V i corresponding to each g i and geodesic Voronoitessellation { ˆ V i } ki =1 on M . For each geodesic Voronoi region ˆ V i , its nominal mass centroid m i is defined by Eq.(9).We call the geodesic Voronoi tessellation as geodesic centroidal Voronoi tessellation (GCVT) if all the generators arenominal mass centroids, i.e. g i = m i , i = 1 , . . . , k. (10)We define an energy function ε g from any tessellation { V i } ki =1 of M and k points { p i } ki =1 ∈ M : ε g ( { ( g i , V i ) } ki =1 ) = k (cid:88) i =1 (cid:90) x ∈ V i ρ ( x ) d g ( x, g i ) dx. (11)We call ε g the GCVT energy functional . It can be shown that the necessary condition for ε g being minimized is that { ( g i , ˆ V i ) } ki =1 is a GCVT [12] so we can obtain a GCVT by optimizing the GCVT energy functional. The theoreticalresults for the combinatorial structures of geodesic Voronoi tessellations will be presented in Section 5 and the algorithmfor finding a GCVT will be introduced in Section 6. At the first glance, CVT, CCVT and GCVT are all uniform tessellations with similar formulations and we call the Voronoiregions of any tessellation as cells . CVT is defined in Euclidean spaces and its cells are convex polygon/polyhedra thatare the intersection of half spaces. Both CCVT and GCVT are defined on manifold domains. However, CCVT needsthat the manifold is embedded in a Euclidean space such that the Euclidean metric d E can be used in Eq.(6). Thenwe say CCVT is extrinsic , which could also be interpreted as the intersection between the manifold and the Voronoi3 PREPRINT - J
ULY
2, 2019 :𝐼 → 𝑅 𝐼 𝑅 (𝐼) 𝑅 A pixel 𝑝 (𝑝) (a) Stretching map (a) Map a grey image I to image manifold M R (b) Computing GCVT on M and content sensitive superpixels on I GCVT - (b) Content-sensitive superpixels via image manifold Figure 3: (a) The map Φ that stretches an image I ⊂ R to an image manifold M ⊂ R . (b) A uniform tessellation on M such as GCVT naturally induces content-sensitive superpixels in I via Φ − ; for an easy illustration, a gray image isused for a mapping in R .tessellation in Euclidean space. Therefore, the cells in CCVT may be disconnected or multiply-connected ; see Figures1 and 2 for an illustration.GCVT is defined on a compact differentiable manifold without referring to an embedded Euclidean space. Since GCVTonly relies on the geodesic metric, it is intrinsic . Compared to CCVT, all the cells in GCVT are guaranteed to beconnected. Before we present theoretic and algorithmic results of GCVT, we state three representative applications, showing thatGCVT is a useful tool in computer vision and graphics.
Image pixels are only the units of image capturing device, but not optimized for image content presentation. Superpixelsare a dense over-segmentation of image, which capture well image features and can serve as perceptually meaningfulatomic regions for images. Superpixels can be used as a preprocessing for reducing the complexity of subsequentimage processing tasks, which includes segmentation[13], contour closure[14], object location [15], object tracking[16], stereo 3D reconstruction [17], and many others. See [18] for a comprehensive survey.To serve as perceptually meaningful atomic regions, superpixels generally have the following characteristics [12]:(1)
Partition : each pixel in the image is assigned to exactly one superpixel so superpixels are a partition of the image;(2)
Connectivity : each region of superpixel is simply connected;(3)
Compactness : in the non-feature region, superpixels are regular in shape and uniform in size;(4)
Feature preservation : superpixels should adhere well to image boundaries for preserving feature;(5)
Content sensitivity : the density of superpixels is adaptive to the variety of image contents;(6)
Performance : superpixels should be computed in a low cost of time and space.In [19, 12], we propose an image manifold that maps a color image I from R to a 2-manifold M embedded in the5-dimensional combined image and colour space R : Φ( u, v ) = ( λ p, λ c ) = ( λ u, λ v, λ l, λ a, λ b ) , (12)where I ( u, v ) is a color image with pixel positions p = ( u, v ) , c ( p ) = ( l ( u, v ) , a ( u, v ) , b ( u, v )) is the color at the pixel p in CIELAB color space, λ and λ are global stretching factors. The area elements in the image manifold M are agood measure of the content density in the image I . Then a uniform tessellation such as GCVT on M naturally inducegood content-sensitive superpixels in I . See Figure 3 for an illustration.GCVT is a powerful tool for superpixels due to the following reasons:(1) Partition : GCVT is a tessellation of M (and also I due to the one-to-one mapping Φ );(2) Connectivity : each cell in GCVT is guaranteed to be connected; A region is simply connected if any simple closed curve in it can be continuously shrunk into a point without leaving the region.A connected region that is not simply connected is called multiply connected . PREPRINT - J
ULY
2, 2019
Original image MSLIC IMSLIC
Figure 4: Some examples of two content-sensitive superpixels (MSLIC [19] and IMSLIC [12]) based on the uniformtessellation on the image manifold. : → 𝑅 𝑅 ( ) 𝑅 Figure 5: The stretching map
Φ : Υ → M ⊂ R maps a video Υ ⊂ R into a 3-manifold M ⊂ R .(3) Compactness : cells in GCVT are regular and uniform in non-feature regions;(4)
Feature preservation : the feature regions (such as object boundary) have a large color variation and therefore lead toa large stretching/area on M . The larger the area in M , the higher possibility that a cell boundary passes throughit;(5) Content sensitivity : content-dense regions have high intensity or color variation, and then larger area on M . Givenuniform tessellation on M , the superpixels will be smaller in content-dense regions. Similarly, content-sparseregions have large superpixels;(6) Performance : we propose efficient computation methods in Section x that can quickly approximate GCVTs.See Figure 4 for some qualitative results of content-sensitive superpixels.
Akin to superpixels for images, Supervoxels are perceptually meaningful atomic regions in videos, obtained bygrouping similar voxels that exhibit coherence in both appearance and motion. Superpixels over-segment a video in thespatiotemporal domain while well preserving its structural content.To compute content-sensitive supvoxels, Yi et al. [20] extend the image manifold concept to the video manifold usingthe stretching map
Φ : Υ → M ⊂ R (Figure 5): Φ(Υ) = Φ( u, v, t ) = ( λ u, λ v, λ t, λ c ) , (13)5 PREPRINT - J
ULY
2, 2019 F r a m e F r a m e F r a m e Original frame Supervoxels
Figure 6: Examples of content-sensitive supervoxels [20] based on the uniform tessellation on the video manifold. Foran easy illustration, supervoxels are clipped in each frame and shown as cross-sectional superpixels. (a) GCVT (b) Dual IDT
Figure 7: (a) GCVT with a small number of generators [21]. (b) Its dual intrinsic Delaunay triangulation (IDT) [22].where the 3-manifold M is embedded in the 6-dimensional combined video and colour space R , Υ is a videowith N voxels, υ ( u, v, t ) ∈ Υ is a voxel with frame index t and the pixel position ( u, v ) in the frame, c ( υ ) =( l ( u, v, t ) , a ( u, v, t ) , b ( u, v, t )) is the color of υ ( u, v, t ) in CIELAB color space λ , λ and λ are global stretchingfactors.At the place where the color variation is large in Υ , Φ maps a unit voxel into a large volume in M . Therefore,the volume elements in M offer a good measure of the content density in Υ . In a similar way to content-sensitivesuperpixels, a uniform tessellation such as GCVT on M will naturally induce content-sensitive supervoxels in Υ , i.e.,supervoxels are typically larger and longer in content-sparse regions (i.e., with homogeneous appearance and motion),and smaller and shorter in content-dense regions (i.e., with high variation of appearance and/or motion). See Figure 6for an illustration. In computer graphics, 3D shapes are usually represented by triangular 2-manifold. In many engineering applicationssuch as finite element analysis, a high quality mesh with almost congruent triangles is desired. To convert an arbitrarytriangular mesh into a high quality mesh while preserving geometric shapes, remeshing techniques are developed.Low-resolution remeshing is to generate a mesh with a small number of vertices and the vertex size approaching thefeature size of the original high-resolution mesh. Due to the existence of thin-shell structures, Voronoi tessellationsbased on Euclidean metric frequently results in disjoint fragments in a Voronoi cell. See Figure 1b for an illustration.GCVT can guarantee that each Voronoi cell is connected and thus is suitable for this low-resolution remeshing task. Aglobally optimized GCVT [21] can generate uniform tessellations with regular cells on a given high-resolution mesh.6
PREPRINT - J
ULY
2, 2019 p p p (a) (b) (c) (d) (e)Geodesic paths Geodesic circles Figure 8: Geodesic paths and circles on curved manifolds. There are two different geodesic circles passing throughthree points p , p and p on M : one is in the front view (d) and the other is in the back view (e).Figure 9: A visibility wedge (VW) on the bottom mesh edge.We propose an efficient sampling criterion such that the intrinsic Delaunay triangulation due to the GVT exists [22].Therefore, it provides an efficient solution to low-resolution remeshing. See Figure 7 for an illustration. Geodesic paths are locally shortest paths between any two points on the manifold. Due to the bending of non-zerocurvatures on curved manifolds M , the distance field on M characterized by geodesic distances/paths have quitedifferent structures from that in Euclidean space R n , such as: • Between any two non-duplicated points in R n , there is one and only one shortest path. However, between anytwo non-duplicated points on M , there may be one, two or infinite shortest paths; See Figure 8 (a-c) for anexample. Only in a smooth, simply connected 2-manifold with negative Gaussian curvature everywhere, thegeodesic path betwteen any two points on M is unique. • Given three points not lying on the same line in R n , there is a unique circle passing through them. However,given three points not lying on the same geodesic path on M , there may be no or more than one geodesiccircles passing through them; See Figure 8 (d-e) for an example.Therefore, the combinatorial structure of geodesic Voronoi tessellations on M are quite different from those in R n .Below we summarize the study of combinatorial structures on 2-manifold meshes M in hierarchical way [11, 23, 22]. Mitchell et al. [7] establish the discrete geodesic structure on M (ref. Figure 9): • Inside every triangle in M , geodesic paths are straight line segments; • When crossing a mesh edge e , geodesic paths are straight lines if two adjacent faces of e are unfolded in thesame plane along e ; 7 PREPRINT - J
ULY
2, 2019
Singular points
An iso-contourSource An obstacle on M Pseudo-sources (a) Isocontour structure
SourceSingular points Iso-contours shown in black or red curves (b) Real example
Figure 10: Each iso-contour of the distance field on a closed M consists of one or more closed curves and each closedcurve consists of circular arc segments joined at singular points. The obstacle in (a) can be a mountain shape with asufficient height. Source p Source q An obstacle on M Break points B( p , q )Pseudo-sources (a) Bisector structure (b) Real example Figure 11: (a) If all vertices do not have the same geodesic distance to a given pair of source points, their bisectorconsists of 1D curve segments, C jointed at breakpoints. Between two breakpoints, the bisector portion is a line orhyperbolic segment. (b) A real example: the geodesic distance field (left) and bisectors between two points (right). • Starting from the triangle that contain the source point, a visibility wedge (VW) can be initialized andpropagated across edges until all edges in M are covered; • The mesh vertices lying on any geodesics (i.e., the apexes of any VWs) are called pseudo-sources , which canonly be saddle vertices, i.e., the vertices in M whose sum of surrounding angles is not smaller than π .The VW structure proposed in [7] can efficiently answer the single-source-all-destination discrete geodesic problem. Given one or more source points P = { p i } Ki =1 , the geodesic distance is d g ( x ) = min i { d g ( x, p i ) , p i ∈ P } , ∀ x ∈ M .An iso-contour (a.k.a. level set) of the distance field d g ( x ) is the trace of all points on M that have the same distancevalue. Iso-contours had drawn considerable attention in literature. On 2-manifold meshes M , their analytical structuresare studied in [11] (ref. Figure 10): • Due to the existance of pseudo-sources, each iso-contour of the distance field on a closed M consists of oneor more closed curves; • Each closed curve consists of circular arc segments joined at singular points, which are locations where thenearest pseudo-source is changing from one to another; • The number of closed curves in an isocontour depends on the indices of critical points of the distance fieldfunction, where a point c ∈ M is a critical point of the distance field function D , if the partial derivatives of D vanish at c . The index of a critical point c is the number of negative eigenvalues of a Hessian matrix of D at c . The bisector between any two points on M is the trace of all points that have equal geodesic distance to these twopoints. On 2-manifold meshes M , the analytical structure of bisectors are studied in [11] (ref. Figure 11):8 PREPRINT - J
ULY
2, 2019
A singular vertex v s 𝑝 𝑞 Figure 12: If a saddle vertex v s lies on the bisector of two points p and q , i.e., d g ( p, v s ) = d g ( q, v s ) , this bisectorcontains a 2D regions (shaded area).Figure 13: The geodesic Voronoi tessellation of the generator set { s i } mi =1 . The geodesic Voronoi cell of s m has m − closed Voronoi edges, each of which is a bisector between s m and s i , i = 1 , , · · · , m − . • The bisector between any two points p and q on M may not be 1D. If a saddle vertex lies on the bisector, thenthis bisector contains a 2D region on M . See Figure 12. Another example can be found in [24] (Figure 4); • Upon small perturbation on the vertices of M , we can assume all vertices do not have the same geodesicdistance to a given pair of source points, and then their bisector consists of 1D curve segments; • If the bisector of two points is 1D, this bisector can have at most g + 1 disjoint closed curves, where g is thegenus of M ; • For each closed curve in a bisector, it can be decomposed at breakpoints, which are the locations where thenearest pseudo-source is changing along the bisector. The bisector is only C at breakpoints. Between twoadjacent breakpoints, the bisector portion can only be line or hyperbolic segment. Given a set of generators, the trimmed bisectors among them partition M into geodesic Voronoi cells. Their analyticalstructures are studied in [11] (ref. Figures 13 and 14): • Each geodesic Voronoi cell is connected, but may not be singly connected; • Each geodesic Voronoi cell is bounded by one or more closed curves called Voronoi edges. Each Voronoi edgeconsists of trimmed bisectors. Trimmed bisectors are joined at branch points that are locations on M havingthe same geodesic distance to its three closest generators. One Voronoi edge does not have to contain a branchpoint.It is well known that in R , given a set of n generators, there are at most n − Voronoi edges and n − branch points(also called Voronoi vertices ) in Voronoi tessellation. The combinatorial structure of geodesic Voronoi tessellation on M is studied in [23]: 9 PREPRINT - J
ULY
2, 2019Figure 14: A real example of geodesic Voronoi tessellation of 12 point generators on a 2-manifold mesh with 2022faces. Left is the geodesic Voronoi tessellation on the mesh. Right shows the boundary of each geodesic Voronoi cellwith transparent surface rendering.Figure 15: Some examples of geodesic Voronoi tessellations. • On a genus- M , the number of Voronoi vertices and Voronoi edges is O ( m ) , where m is the number ofgenerators; • The combinatorial complexity of geodesic Voronoi tessellation is defined to be the total number of breakpoints,Voronoi vertices, Voronoi edges and Voronoi cells. On a genus- M , the combinatorial complexity of geodesicVoronoi tessellation is O ( mk ) , where k is the number of faces in M ; • On a genus- g M , the number of Voronoi vertices and Voronoi edges is O ( m + g ) , where g is the genus of M ; • If the set of generators is dense and the geodesic Voronoi tessellation satisfies the closed ball property [25], thecombinatorial complexity of geodesic Voronoi tessellation on a genus- g M is O (( n + g ) k ) .Some real examples of geodesic Voronoi tessellations are illustrated in Figure 15. Given the number N of point generators on a d -manifold mesh, the GCVT can be computed by finding a tessellationthat minimizes the GCVT energy in (11). In this section, we present algorithms focusing on d = 2 . Some methodssummarized in this section (e.g., the RCVT approximation method in Section 6.2.2) can be extended to arbitrary d ≥ dimensions.The energy in (11) can be minimized globally or locally. There are two existing global optimization methods for CVT:one is the Monte Carlo with minimization (MCM) framework [26] that only deals with CVT in R and the other is themanifold differential evolution (MDE) method [21] that is general to deal with GCVT on M ( R is a special case of M ). MCM is a heuristic method without theoretical guarantee, while MDE has a provable probabilistic convergence tothe global optimum. In Section 6.1, we summarize the MDE global method.Although the global method can achieve high-quality GCVT results and is insensitive to the initial position of generators,it is very time-consuming. Therefore, several fast, local optimization methods had proposed and we summarize twolocal methods in Section 6.2. 10 PREPRINT - J
ULY
2, 2019 X i ,0 X i ,rand0 X i ,rand2 V UDifference of two agents Add difference to an agent
Crossover
Randomly select three agentsInitial population X i ,rand1 X i ,0 Selection X i ,0 X i ,rand0 X i ,rand1 X i ,rand2 D = X i ,rand0 ⊖ X i ,rand1 V = D ⊕ X i ,rand2 V U = Crossover ( V ,X i ,0 ) X i +1,0 = Selection ( U , X i ,0 ) X i ,0 Figure 16: MDE pipeline.
MDE is a stochastic global optimization method, which extends the classic differential evolution [27] to the manifoldsetting.Classic differential evolution applies agents that have operations of addition, subtraction and scala multiplication, alldefined in a vector space. To use these operations on a manifold, MDE assigned an order to the generators in an agentfor encoding them into a vector representation such that different agents can be matched akin to matching vectors. Thepipeline of MDE is illustrated in Figure 16. It initializes a group of agents by random sampling and improves the qualityof the agents iteratively. There are three steps in each iteration of MDE, i.e. mutation, crossover and selection.
Vector representation and agent matching.
Given any two agents G = { g i } ni =1 and G (cid:48) = { g (cid:48) i } ni =1 , MDE builds acomplete bipartite graph K n,n whose vertices are G ∪ G (cid:48) and edges weights are the geodesic distances between the pairof any two g i and g (cid:48) j . MDE finds a perfect matching in K n,n by solving the minimum-weight perfect matching problemusing the Hungarian algorithm [28], which runs in O ( n ) time. An order of generators can be induced by this matching,and thus the subscripts of generators can be corresponding by rearrangement. Let X k,j = { x k,j,i } Ni =1 , j = 1 . . . M bethe j -th agent of k -th generation, where N is the number of generators and M is the number of agents. Mutation operator.
MDE produces M new mutative agents by the mutation operator. The j -th mutative agent can beobtained by randomly selecting three agents and adding a scaled difference between two agents to another agent, i.e. V k,j = X k,rand ⊕ ( λ ⊗ ( X k,rand (cid:9) X k,rand )) , (14)where (cid:9) ( y, x ) : M × M → T x M outputs a tangent vector at x whose direction is the starting direction of the geodesicpath from x to y and whose magnitude is length of the geodesic path, ⊗ is the scala multiplication in the tangent spaceand ⊕ ( z, v ) : M × T ( M ) → M obtains a point on the manifold by (1) parallel transporting v to z along the geodesicpath from the point of v to z , and denote the resulting vector as v (cid:48) ; (2) computing a geodesic path from z with a initialdirection v (cid:48) and the length is equal to v (cid:48) , and the end point of this geodesic path is the output. Crossover operator.
For each agent X k,j , MDE produces a competitor U k,j = { u k,j,i } Ni =1 using the correspondingmutative agent V k,j by crossover operator. u k,j,i randomly uses the corresponding component of X k,j or V k,j , i.e., u k,j,i = (cid:26) v k,j,i rand (0 , < C r x k,j,i otherwise , (15)where C r is crossover rate . Selection operator.
MDE selects a better agent from X k,j and U k,j and puts it into the new generation, i.e., X k +1 ,j = (cid:26) X k,j ε ( X k,j ) < ε ( U k,j ) U k,j otherwise , (16)where ε is the GCVT energy of the agent computed by Equation 11.The terminate condition of MDE is meeting one of the three conditions: (1) the iteration number exceeds the parameterspecified by user; (2) the solution does not improve in successive several iterations; and (3) the GCVT energy reachesthe prescribed value. 11 PREPRINT - J
ULY
2, 2019 (a) Regular mesh (b) Highly irregular mesh
Figure 17: The regular mesh in (a) left and the irregular mesh in (b) left represent the same geometry. By applying theglobally optimized MDE solution, the GCVTs on both meshes with the same number of point generators are the same.It was shown in [21] that under some mild assumptions, the MDE solution converges to the global optimization withprobability 1.The globally optimized MDE solution is insensitive to the underlying mesh quality. See Figure 17 for an example.
The Lloyd method [9] is a classic algorithm that can be used to efficiently compute the cluster problem, including theconstruction of CVT and GCVT. The Lloyd method locally minimizes the GCVT energy by iteratively moving thegenerators to the corresponding nominal mass centroids and updating the GVT of these generators. I.e., in each iteration,there are two steps in Lloyd method: one is fixing the tessellation and the moving generators to the correspondingnominal mass centroids and the other is fixing the generators and updating the GVT. Both of the two steps reduce theGCVT energy [12], which ensures the convergence of the Lloyd algorithm.The bottleneck of this Lloyd method lies on the computation of the nominal mass centroids, which requires to solve theoptimization problem in (9). It is difficult or even impossible to solve this problem analytically and thus approximationhas to be used. Wang et al. [29] propose an approximation method that computes the Riemannian center instead ofsolving the problem (9). Let v , v , · · · , v m be the Voronoi vertices of a geodesic Voronoi cell on M . The Riemanniancenter is defined as the local minima of the following function: E ( x ) = m (cid:88) i =1 d g ( x, v i ) (17)Based on the properties studied in [30, 31], an iterative method utilizing the exponential map is proposed in [29] toquickly find an approximation to the Riemannian center. Another approximation to the nominal mass centroid isproposed in [12] that makes use of the landmark MDS (LMDS) [32] to quickly unfold a geodesic Voronoi cell into R in a way such that the total distance distortion defined by a graph embedding is minimized. Then the nominal masscentroid is approximated by the mass centroid (ref. Eq.(2)) of unfolded geodesic Voronoi cell in R . The two local methods summarized in Section 6.2.1 make use of geodesic distance, which is time-consuming to computefor updating GVT in each Lloyd iteration.Restricted centroidal Voronoi tessellation(RCVT) — which utilizes Euclidean distance in embedded Euclidean space —is proposed in [19, 20] as a fast approximation to GCVT. Different from the two local methods in Section 6.2.1 that canonly handle the 2-manifold meshes, the RCVT summarized in this section can deal with any d -dimensional triangulatedmeshes, d ≥ . 12 PREPRINT - J
ULY
2, 2019Given a d -manifold M embedded in R n , d < n , and a set of point generators { g i } ki =1 ∈ R n (not necessary on M ), the restricted Voronoi region corresponding to the point g i is defined as ˆ V i = { x ∈ M (cid:12)(cid:12) d E ( x, g i ) < d E ( x, g j ) for j = 1 , . . . , k, j (cid:54) = i } . (18)Elements of { g i } ki =1 are called generators . The set of restricted Voronoi regions { ˆ V i } ki =1 forms a restricted Voronoitessellation of M . The mass centroid m i of ˆ V i is defined by m i = (cid:82) x ∈ ˆ V i xdx (cid:82) x ∈ ˆ V i dx . (19)Similar to point generators { g i } ki =1 , the mass centroid m i does not need to be on M . The restricted Voronoi tessellationis called restricted centroidal Voronoi tessellation(RCVT) if all the generators are mass centroids, i.e. g i = m i , i = 1 , . . . , k. (20)Due to the following properties: • Compared to GCVT, RCVT uses Euclidean distance and • Compared to CCVT, the mass centroids do not need to be on M ,RCVT behaves as a natural bridge between GCVT and CCVT. Finally, RCVT is easy to compute using the Lloydmethod with Euclidean distance. In this section, we provide an analysis on the time complexity and space complexity of the above methods forconstructing GCVT.In MDE, the agent matching runs in O ( m ) time and computing the exact geodesic between two points runs in O ( n log n ) time, where m is the number of generators and n is the number of vertices. Therefore, generating anagent takes O ( m + n log n ) time and MDE runs in O ( N p K ( m + n log n )) , where N p is population size and K isiteration number. It takes O ( N p m ) space for N p agents which have N p m generators in total.IMSLIC [12] and MSLIC [19] construct GCVT and RCVT for an image using Lloyd method. It alternately constructsGVT and computes mass centroids. We analyze the complexities of the two steps separately. Fixing generators, IMSLICcan construct a GVT corresponding to the generators using Dijkstra’s algorithm in O ( n log n ) time and O ( n ) space.Adopting a label correcting method that maintains a bucket data structure can reduce the complexity and it only takes O ( n ) time. Fixing generators, MSLIC can construct a RVT corresponding to the generators by going through all pixelsin O ( n ) time and O ( n ) space. Fixing the tessellation, IMSLIC computes approximate nominal mass using LMDS in O ( n ) time and O ( n ) space and MSLIC computes mass centroids by going through all pixels in O ( n ) time and O (1) space. Therefore, we can obtain an approximate GCVT in O ( nK ) time and O ( n ) space and obtain an RCVT in O ( nK ) time and O ( n ) space for an image, where K is the number of iteration. Geodesic centroidal Voronoi tessellations (GCVTs) are intrinsic geometric structure inherent in manifolds. In this paper,we summarized our recent work on GCVTs on triangulated manifold. We show our results from both theoretical andalgorithmic aspects. Their applications in computer graphics and vision are also presented. We hope that this paper canprovide some insights for researchers to review the past developments and identify directions for future research on thestudy of intrinsic geometric structures in intelligent media data processing.
References [1] Atsuyuki Okabe, Barry Boots, Kokichi Sugihara, and Sung Nok Chiu.
Spatial Tessellations: Concept andApplications of Voronoi Diagrams . Wiley, 2000.[2] Jeffrey M Augenbaum and Charles S Peskin. On the construction of the voronoi mesh on a sphere.
Journal ofComputational Physics , 59(2):177–192, 1985.[3] Kensuke Onishi and Nobuki Takayama. Construction of voronoi diagram on the upper half-plane.
IEICE Trans.Fundamentals of Electronics, Communications and Computer Sciences , E79-A(4):533–539, 1996.13
PREPRINT - J
ULY
2, 2019[4] Greg Leibon and David Letscher. Delaunay triangulations and voronoi diagrams for riemannian manifolds. In
Proceedings of the Sixteenth Annual Symposium on Computational Geometry , SCG ’00, pages 341–349. ACM,2000.[5] H. Sebastian Seung and Daniel D. Lee. The manifold ways of perception.
Science , 290(5500):2268–2269, 2000.[6] William M Boothby.
An introduction to differentiable manifolds and Riemannian geometry . Academic Press,1975.[7] Joseph SB Mitchell, David M Mount, and Christos H Papadimitriou. The discrete geodesic problem.
SIAMJournal on Computing , 16(4):647–668, 1987.[8] Qiang Du, Vance Faber, and Max Gunzburger. Centroidal voronoi tessellations: Applications and algorithms.
Siam Review , 41(4):637–676, 1999.[9] Stuart P. Lloyd. Least squares quantization in pcm’s.
IEEE Transactions on Information Theory , 28:129–136, 031982.[10] Qiang Du, Max D. Gunzburger, and Lili Ju. Constrained centroidal voronoi tessellations for surfaces.
SIAMJournal on Scientific Computing , 24(5):1488–1506, 2003.[11] Liu Yong-Jin, Chen Zhan-Qing, and Tang Kai. Construction of iso-contours, bisectors, and voronoi diagrams ontriangulated surfaces.
IEEE Transactions on Pattern Analysis & Machine Intelligence , 33(8):1502–1517, 2011.[12] Yong-Jin Liu, Minjing Yu, Bing-Jun Li, and Ying He. Intrinsic manifold SLIC: A simple and efficient method forcomputing content-sensitive superpixels.
IEEE Trans. Pattern Anal. Mach. Intell. , 40(3):653–666, 2018.[13] L. Ming Yu, O. Tuzel, S. Ramalingam, and R. Chellappa. Entropy rate superpixel segmentation. In
ComputerVision & Pattern Recognition , 2011.[14] Alex Levinshtein, Cristian Sminchisescu, and Sven Dickinson. Optimal contour closure by superpixel grouping.In
European Conference on Computer Vision , 2010.[15] Brian Fulkerson, Andrea Vedaldi, and Stefano Soatto. Class segmentation and object localization with superpixelneighborhoods. In
IEEE International Conference on Computer Vision , 2009.[16] Wang Shu, Huchuan Lu, Yang Fan, and Ming Hsuan Yang. Superpixel tracking. In
International Conference onComputer Vision , 2011.[17] Branislav Miˇcušík and Jana Košecká. Multi-view superpixel stereo in urban environments.
International Journalof Computer Vision , 89(1):106–119, 2010.[18] David Stutz, Alexander Hermans, and Bastian Leibe. Superpixels: An evaluation of the state-of-the-art.
ComputerVision and Image Understanding , 166:1–27, 2018.[19] Yong-Jin Liu, Chengchi Yu, Minjing Yu, and Ying He. Manifold SLIC: a fast method to compute content-sensitivesuperpixels. In
IEEE Conference on Computer Vision and Pattern Recognition , CVPR ’16, pages 651–659, 2016.[20] Ran Yi, Yong-Jin Liu, and Yu-Kun Lai. Content-sensitive supervoxels via uniform tessellations on video manifolds.In
IEEE Conference on Computer Vision and Pattern Recognition , CVPR ’18, pages 646–655, 2018.[21] Yong-Jin Liu, Chun-Xu Xu, Ran Yi, Dian Fan, and Ying He. Manifold differential evolution (MDE): A globaloptimization method for geodesic centroidal voronoi tessellations on meshes.
ACM Transactions on Graphics(SIGGRAPH ASIA 2016) , 35(6):243:1–243:10, 2016.[22] Yong-Jin Liu, Dian Fan, Chun-Xu Xu, and Ying He. Constructing intrinsic delaunay triangulations from the dualof geodesic voronoi diagrams.
ACM Trans. Graph. , 36(2):15:1–15:15, 2017.[23] Yong-Jin Liu and Kai Tang. The complexity of geodesic voronoi diagrams on triangulated 2-manifold surfaces.
Information Processing Letters , 113(4):132–136, 2013.[24] Yong-Jin Liu. Semi-continuity of skeletons in 2-manifold and discrete voronoi approximation.
IEEE Transactionson Pattern Analysis and Machine Intelligence , 37(9):1938–1944, 2015.[25] Herbert Edelsbrunner and Nimish R. Shah. Triangulating topological spaces.
International Journal of Computa-tional Geometry and Applications , 7(4):365–378, 1997.[26] Lin Lu, Feng Sun, Hao Pan, and Wenping Wang. Global optimization of centroidal voronoi tessellation withmonte carlo approach.
IEEE Trans. Vis. Comput. Graph. , 18(11):1880–1890, 2012.[27] Swagatam Das and Ponnuthurai Nagaratnam Suganthan. Differential evolution: A survey of the state-of-the-art.
IEEE Transactions on Evolutionary Computation , 15(1):4–31, 2011.[28] B B. H. Korte and Jens Vygen.
Combinatorial Optimization: Theory and Algorithms . Spinger, 2012.14
PREPRINT - J
ULY
2, 2019[29] Xiaoning Wang, Xiang Ying, Yong-Jin Liu, Shi-Qing Xin, Wenping Wang, Xianfeng Gu, Wolfgang Müller-Wittig,and Ying He. Intrinsic computation of centroidal voronoi tessellation (CVT) on meshes.
Computer-Aided Design ,58:51–61, 2015.[30] Xavier Pennec. Intrinsic statistics on riemannian manifolds: Basic tools for geometric measurements.
Journal ofMathematical Imaging and Vision , 25(1):127–154, 2006.[31] Raif M. Rustamov. Barycentric coordinates on surfaces.
Comput. Graph. Forum , 29(5):1507–1516, 2010.[32] Vin de Silva and Joshua B. Tenenbaum. Global versus local methods in nonlinear dimensionality reduction. In