Let be a metric space, where is a point set and is a metric. Given an open subset , a set is called a tessellation of if for and , where is the closure of . Given a set of points in , the Voronoi cell corresponding to the point is defined as
Elements of are called generators.
Since 1644 (Part III of Principia Phiolosophiae, written by Descartes), Voronoi tessellations had been well studied in the Euclidean space , , in which Voronoi tessellations (as well as their dual structures, well known as Delaunay triangulations) had played a central role as fundamental geometric structures. Voronoi tessellations had also been studied in spaces with non-Euclidean metrics, including spheres , hyperbolic spaces  and the general Riemannian manifolds . Recently, due to the flourishes of big media data from digital sampling, more and more data are appearing in the form of manifold meshes . Quite different from parametrized 2-manifolds or general Riemannian manifolds that are generally smooth (or at least smooth) , manifold meshes are only . In our study, we adopt the discrete geodesic metric .
In this paper, we summarize our recent work on geodesic centroidal Voronoi tessellations (GCVT) — which are provable uniform tessellations on manifold meshes — and we show that they can be used to generate uniform remeshing in computer graphics and build content-sensitive superpixels/supervoxels for images and video in computer vision applications.
Before we introduce GCVTs, we present two close concepts of CVT and RCVT in related work.
2 Related Work
2.1 Centroidal Voronoi Tessellations (CVT)
CVT had been well studied in science and engineering, with a wide range of applications including data compression in digital 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 for an excellent survey.
Let V be a finite region in . The mass centroid of is defined by
where is a density function defined in . Given points in a domain , we can define the Voronoi region corresponding to each based on the Euclidean metric and the Voronoi tessellation of . Let be the mass centroid of each Voronoi region . A Voronoi tessellation is called CVT if all the generators are mass centroids, i.e.
Arbitrarily chosen generators are usually not the mass centroids of their associated Voronoi regions so an arbitrary Voronoi tessellation cannot be a CVT. It can be shown that CVT minimizes the following energy (a.k.a. CVT energy functional):
where is a density function defined in , is an arbitrary tessellation and is any set of points in . To compute CVT, various local methods had been surveyed in . In particular, the Lloyd method  is a simple yet effective local method that iteratively computes mass centroids and Voronoi tessellation.
2.2 Constrained CVT for Surfaces
To extend the domain partitioning of CVT from Euclidean space to general spaces, Du et al.  proposed the constrained CVT (CCVT) that works on a compact and continuous surface defined by
where are some continuous functions.
Given a set of points , the constrained Voronoi region corresponding to each based on the Euclidean metric and restricted in is defined as
The constrained mass centroid of on is defined to be the solution of the following problem:
A tessellation on is called CCVT if and only if the set are both the generators and constrained mass centroids of the tessellation .
3 Geodesic Centroidal Voronoi Tessellations (GCVT)
Given a k-dimensional compact differentiable manifold , a set is called a tessellation of if for and . Denote by the geodesic distance between and in . Given a set of points in , the geodesic Voronoi region corresponding to the point is defined as
Elements of are called generators. The set of geodesic Voronoi regions is called geodesic Voronoi tessellation of . A geodesic Voronoi region is a connected domain and is a non-empty compact set . For each geodesic Voronoi region , the nominal mass centroid of on is defined to be the solution of the following problem
It can be shown that is continuous and the domain is compact so has at least one global minimum in . Therefore, the nominal mass centroid exists.
Given points , we can define the geodesic Voronoi region corresponding to each and geodesic Voronoi tessellation on . For each geodesic Voronoi region , its nominal mass centroid is defined by Eq.(9). We call the geodesic Voronoi tessellation as geodesic centroidal Voronoi tessellation (GCVT) if all the generators are nominal mass centroids, i.e.
We define an energy function from any tessellation of and points :
We call the GCVT energy functional. It can be shown that the necessary condition for being minimized is that is a GCVT  so we can obtain a GCVT by optimizing the GCVT energy functional. The theoretical results for the combinatorial structures of geodesic Voronoi tessellations will be presented in Section 5 and the algorithm for finding a GCVT will be introduced in Section 6.
3.1 Comparison of CVT, CCVT and GCVT
At the first glance, CVT, CCVT and GCVT are all uniform tessellations with similar formulations and we call the Voronoi regions of any tessellation as cells. CVT is defined in Euclidean spaces and its cells are convex polygon/polyhedra that are the intersection of half spaces. Both CCVT and GCVT are defined on manifold domains. However, CCVT needs that the manifold is embedded in a Euclidean space such that the Euclidean metric can be used in Eq.(6). Then we say CCVT is extrinsic, which could also be interpreted as the intersection between the manifold and the Voronoi tessellation in Euclidean space. Therefore, the cells in CCVT may be disconnected or multiply-connected111A 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.; see Figures 1 and 2 for an illustration.
GCVT is defined on a compact differentiable manifold without referring to an embedded Euclidean space. Since GCVT only relies on the geodesic metric, it is intrinsic. Compared to CCVT, all the cells in GCVT are guaranteed to be connected.
Before we present theoretic and algorithmic results of GCVT, we state three representative applications, showing that GCVT is a useful tool in computer vision and graphics.
4.1 Content-Sensitive Superpixels
Image pixels are only the units of image capturing device, but not optimized for image content presentation. Superpixels are a dense over-segmentation of image, which capture well image features and can serve as perceptually meaningful atomic regions for images. Superpixels can be used as a preprocessing for reducing the complexity of subsequent image processing tasks, which includes segmentation, contour closure, object location , object tracking , stereo 3D reconstruction , and many others. See  for a comprehensive survey.
To serve as perceptually meaningful atomic regions, superpixels generally have the following characteristics :
Partition: each pixel in the image is assigned to exactly one superpixel so superpixels are a partition of the image;
Connectivity: each region of superpixel is simply connected;
Compactness: in the non-feature region, superpixels are regular in shape and uniform in size;
Feature preservation: superpixels should adhere well to image boundaries for preserving feature;
Content sensitivity: the density of superpixels is adaptive to the variety of image contents;
Performance: superpixels should be computed in a low cost of time and space.
where is a color image with pixel positions , is the color at the pixel in CIELAB color space, and are global stretching factors. The area elements in the image manifold are a good measure of the content density in the image . Then a uniform tessellation such as GCVT on naturally induce good content-sensitive superpixels in . See Figure 3 for an illustration.
GCVT is a powerful tool for superpixels due to the following reasons:
Partition: GCVT is a tessellation of (and also due to the one-to-one mapping );
Connectivity: each cell in GCVT is guaranteed to be connected;
Compactness: cells in GCVT are regular and uniform in non-feature regions;
Feature preservation: the feature regions (such as object boundary) have a large color variation and therefore lead to a large stretching/area on . The larger the area in , the higher possibility that a cell boundary passes through it;
Content sensitivity: content-dense regions have high intensity or color variation, and then larger area on . Given uniform tessellation on , the superpixels will be smaller in content-dense regions. Similarly, content-sparse regions have large superpixels;
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.
4.2 Content-Sensitive Supervoxels
Akin to superpixels for images, Supervoxels are perceptually meaningful atomic regions in videos, obtained by grouping similar voxels that exhibit coherence in both appearance and motion. Superpixels over-segment a video in the spatiotemporal domain while well preserving its structural content.
where the 3-manifold is embedded in the 6-dimensional combined video and colour space , is a video with voxels, is a voxel with frame index and the pixel position in the frame, is the color of in CIELAB color space , and are global stretching factors.
At the place where the color variation is large in , maps a unit voxel into a large volume in . Therefore, the volume elements in offer a good measure of the content density in . In a similar way to content-sensitive superpixels, a uniform tessellation such as GCVT on 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 6 for an illustration.
4.3 Low-Resolution Remeshing
In computer graphics, 3D shapes are usually represented by triangular 2-manifold. In many engineering applications such as finite element analysis, a high quality mesh with almost congruent triangles is desired. To convert an arbitrary triangular 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 the feature size of the original high-resolution mesh. Due to the existence of thin-shell structures, Voronoi tessellations based 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. A globally optimized GCVT  can generate uniform tessellations with regular cells on a given high-resolution mesh. We propose an efficient sampling criterion such that the intrinsic Delaunay triangulation due to the GVT exists . 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-zero curvatures on curved manifolds , the distance field on characterized by geodesic distances/paths have quite different structures from that in Euclidean space , such as:
Between any two non-duplicated points in , there is one and only one shortest path. However, between any two non-duplicated points on , there may be one, two or infinite shortest paths; See Figure 8 (a-c) for an example. Only in a smooth, simply connected 2-manifold with negative Gaussian curvature everywhere, the geodesic path betwteen any two points on is unique.
Given three points not lying on the same line in , there is a unique circle passing through them. However, given three points not lying on the same geodesic path on , there may be no or more than one geodesic circles passing through them; See Figure 8 (d-e) for an example.
Therefore, the combinatorial structure of geodesic Voronoi tessellations on are quite different from those in . Below we summarize the study of combinatorial structures on 2-manifold meshes in hierarchical way [11, 23, 22].
5.1 Discrete geodesics
Inside every triangle in , geodesic paths are straight line segments;
When crossing a mesh edge , geodesic paths are straight lines if two adjacent faces of are unfolded in the same plane along ;
Starting from the triangle that contain the source point, a visibility wedge (VW) can be initialized and propagated across edges until all edges in are covered;
The mesh vertices lying on any geodesics (i.e., the apexes of any VWs) are called pseudo-sources, which can only be saddle vertices, i.e., the vertices in whose sum of surrounding angles is not smaller than .
The VW structure proposed in  can efficiently answer the single-source-all-destination discrete geodesic problem.
5.2 Iso-contour structure
Given one or more source points , the geodesic distance is , . An iso-contour (a.k.a. level set) of the distance field is the trace of all points on that have the same distance value. Iso-contours had drawn considerable attention in literature. On 2-manifold meshes , their analytical structures are studied in  (ref. Figure 10):
Due to the existance of pseudo-sources, each iso-contour of the distance field on a closed consists of one or more closed curves;
Each closed curve consists of circular arc segments joined at singular points, which are locations where the nearest 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 field function, where a point is a critical point of the distance field function , if the partial derivatives of vanish at . The index of a critical point
is the number of negative eigenvalues of a Hessian matrix ofat .
5.3 Bisector structure
The bisector between any two points on is the trace of all points that have equal geodesic distance to these two points. On 2-manifold meshes , the analytical structure of bisectors are studied in  (ref. Figure 11):
Upon small perturbation on the vertices of , we can assume all vertices do not have the same geodesic distance 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 disjoint closed curves, where is the genus of ;
For each closed curve in a bisector, it can be decomposed at breakpoints, which are the locations where the nearest pseudo-source is changing along the bisector. The bisector is only at breakpoints. Between two adjacent breakpoints, the bisector portion can only be line or hyperbolic segment.
5.4 Geodesic Voronoi Tessellation
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 edge consists of trimmed bisectors. Trimmed bisectors are joined at branch points that are locations on having the same geodesic distance to its three closest generators. One Voronoi edge does not have to contain a branch point.
It is well known that in , given a set of generators, there are at most Voronoi edges and branch points (also called Voronoi vertices) in Voronoi tessellation. The combinatorial structure of geodesic Voronoi tessellation on is studied in :
On a genus- , the number of Voronoi vertices and Voronoi edges is , where is the number of generators;
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- , the combinatorial complexity of geodesic Voronoi tessellation is , where is the number of faces in ;
On a genus- , the number of Voronoi vertices and Voronoi edges is , where is the genus of ;
If the set of generators is dense and the geodesic Voronoi tessellation satisfies the closed ball property , the combinatorial complexity of geodesic Voronoi tessellation on a genus- is .
Some real examples of geodesic Voronoi tessellations are illustrated in Figure 15.
Given the number of point generators on a -manifold mesh, the GCVT can be computed by finding a tessellation that minimizes the GCVT energy in (11). In this section, we present algorithms focusing on . Some methods summarized in this section (e.g., the RCVT approximation method in Section 6.2.2) can be extended to arbitrary 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  that only deals with CVT in and the other is the manifold differential evolution (MDE) method  that is general to deal with GCVT on ( is a special case of
). MCM is a heuristic method without theoretical guarantee, while MDE has a provable probabilistic convergence to the global optimum. In Section6.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 two local methods in Section 6.2.
6.1 MDE Global Method
MDE is a stochastic global optimization method, which extends the classic differential evolution  to the manifold setting.
Classic differential evolution applies agents that have operations of addition, subtraction and scala multiplication, all defined in a vector space. To use these operations on a manifold, MDE assigned an order to the generators in an agent for encoding them into a vector representation such that different agents can be matched akin to matching vectors. The pipeline of MDE is illustrated in Figure16. It initializes a group of agents by random sampling and improves the quality of 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 and , MDE builds a complete bipartite graph whose vertices are and edges weights are the geodesic distances between the pair of any two and . MDE finds a perfect matching in by solving the minimum-weight perfect matching problem using the Hungarian algorithm , which runs in time. An order of generators can be induced by this matching, and thus the subscripts of generators can be corresponding by rearrangement. Let be the -th agent of -th generation, where is the number of generators and is the number of agents.
Mutation operator. MDE produces new mutative agents by the mutation operator. The -th mutative agent can be obtained by randomly selecting three agents and adding a scaled difference between two agents to another agent, i.e.
where outputs a tangent vector at whose direction is the starting direction of the geodesic path from to and whose magnitude is length of the geodesic path, is the scala multiplication in the tangent space and obtains a point on the manifold by (1) parallel transporting to along the geodesic path from the point of to , and denote the resulting vector as ; (2) computing a geodesic path from with a initial direction and the length is equal to , and the end point of this geodesic path is the output.
Crossover operator. For each agent , MDE produces a competitor using the corresponding mutative agent by crossover operator. randomly uses the corresponding component of or , i.e.,
where is crossover rate.
Selection operator. MDE selects a better agent from and and puts it into the new generation, i.e.,
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 parameter specified by user; (2) the solution does not improve in successive several iterations; and (3) the GCVT energy reaches the prescribed value.
It was shown in 
that under some mild assumptions, the MDE solution converges to the global optimization with probability 1.
The globally optimized MDE solution is insensitive to the underlying mesh quality. See Figure 17 for an example.
6.2 Two Local Methods
6.2.1 Approximate nominal mass centroid
The Lloyd method  is a classic algorithm that can be used to efficiently compute the cluster problem, including the construction of CVT and GCVT. The Lloyd method locally minimizes the GCVT energy by iteratively moving the generators 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 corresponding nominal mass centroids and the other is fixing the generators and updating the GVT. Both of the two steps reduce the GCVT energy , 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 the optimization problem in (9). It is difficult or even impossible to solve this problem analytically and thus approximation has to be used. Wang et al.  propose an approximation method that computes the Riemannian center instead of solving the problem (9). Let be the Voronoi vertices of a geodesic Voronoi cell on . The Riemannian center is defined as the local minima of the following function:
Based on the properties studied in [30, 31], an iterative method utilizing the exponential map is proposed in  to quickly find an approximation to the Riemannian center. Another approximation to the nominal mass centroid is proposed in  that makes use of the landmark MDS (LMDS)  to quickly unfold a geodesic Voronoi cell into in a way such that the total distance distortion defined by a graph embedding is minimized. Then the nominal mass centroid is approximated by the mass centroid (ref. Eq.(2)) of unfolded geodesic Voronoi cell in .
6.2.2 RCVT approximation
The two local methods summarized in Section 6.2.1 make use of geodesic distance, which is time-consuming to compute for 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 can only handle the 2-manifold meshes, the RCVT summarized in this section can deal with any -dimensional triangulated meshes, .
Given a -manifold embedded in , , and a set of point generators (not necessary on ), the restricted Voronoi region corresponding to the point is defined as
Elements of are called generators. The set of restricted Voronoi regions forms a restricted Voronoi tessellation of . The mass centroid of is defined by
Similar to point generators , the mass centroid does not need to be on . The restricted Voronoi tessellation is called restricted centroidal Voronoi tessellation(RCVT) if all the generators are mass centroids, i.e.
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 ,
RCVT behaves as a natural bridge between GCVT and CCVT. Finally, RCVT is easy to compute using the Lloyd method with Euclidean distance.
6.3 Algorithm Analysis
In this section, we provide an analysis on the time complexity and space complexity of the above methods for constructing GCVT.
In MDE, the agent matching runs in time and computing the exact geodesic between two points runs in time, where is the number of generators and is the number of vertices. Therefore, generating an agent takes time and MDE runs in , where is population size and is iteration number. It takes space for agents which have generators in total.
IMSLIC  and MSLIC  construct GCVT and RCVT for an image using Lloyd method. It alternately constructs GVT and computes mass centroids. We analyze the complexities of the two steps separately. Fixing generators, IMSLIC can construct a GVT corresponding to the generators using Dijkstra’s algorithm in time and space. Adopting a label correcting method that maintains a bucket data structure can reduce the complexity and it only takes time. Fixing generators, MSLIC can construct a RVT corresponding to the generators by going through all pixels in time and space. Fixing the tessellation, IMSLIC computes approximate nominal mass using LMDS in time and space and MSLIC computes mass centroids by going through all pixels in time and space. Therefore, we can obtain an approximate GCVT in time and space and obtain an RCVT in time and space for an image, where 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 and algorithmic aspects. Their applications in computer graphics and vision are also presented. We hope that this paper can provide some insights for researchers to review the past developments and identify directions for future research on the study of intrinsic geometric structures in intelligent media data processing.
-  Atsuyuki Okabe, Barry Boots, Kokichi Sugihara, and Sung Nok Chiu. Spatial Tessellations: Concept and Applications of Voronoi Diagrams. Wiley, 2000.
-  Jeffrey M Augenbaum and Charles S Peskin. On the construction of the voronoi mesh on a sphere. Journal of Computational Physics, 59(2):177–192, 1985.
-  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.
-  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.
-  H. Sebastian Seung and Daniel D. Lee. The manifold ways of perception. Science, 290(5500):2268–2269, 2000.
-  William M Boothby. An introduction to differentiable manifolds and Riemannian geometry. Academic Press, 1975.
-  Joseph SB Mitchell, David M Mount, and Christos H Papadimitriou. The discrete geodesic problem. SIAM Journal on Computing, 16(4):647–668, 1987.
-  Qiang Du, Vance Faber, and Max Gunzburger. Centroidal voronoi tessellations: Applications and algorithms. Siam Review, 41(4):637–676, 1999.
-  Stuart P. Lloyd. Least squares quantization in pcm’s. IEEE Transactions on Information Theory, 28:129–136, 03 1982.
-  Qiang Du, Max D. Gunzburger, and Lili Ju. Constrained centroidal voronoi tessellations for surfaces. SIAM Journal on Scientific Computing, 24(5):1488–1506, 2003.
-  Liu Yong-Jin, Chen Zhan-Qing, and Tang Kai. Construction of iso-contours, bisectors, and voronoi diagrams on triangulated surfaces. IEEE Transactions on Pattern Analysis & Machine Intelligence, 33(8):1502–1517, 2011.
-  Yong-Jin Liu, Minjing Yu, Bing-Jun Li, and Ying He. Intrinsic manifold SLIC: A simple and efficient method for computing content-sensitive superpixels. IEEE Trans. Pattern Anal. Mach. Intell., 40(3):653–666, 2018.
L. Ming Yu, O. Tuzel, S. Ramalingam, and R. Chellappa.
Entropy rate superpixel segmentation.
Computer Vision & Pattern Recognition, 2011.
-  Alex Levinshtein, Cristian Sminchisescu, and Sven Dickinson. Optimal contour closure by superpixel grouping. In European Conference on Computer Vision, 2010.
-  Brian Fulkerson, Andrea Vedaldi, and Stefano Soatto. Class segmentation and object localization with superpixel neighborhoods. In IEEE International Conference on Computer Vision, 2009.
-  Wang Shu, Huchuan Lu, Yang Fan, and Ming Hsuan Yang. Superpixel tracking. In International Conference on Computer Vision, 2011.
-  Branislav Mičušík and Jana Košecká. Multi-view superpixel stereo in urban environments. International Journal of Computer Vision, 89(1):106–119, 2010.
-  David Stutz, Alexander Hermans, and Bastian Leibe. Superpixels: An evaluation of the state-of-the-art. Computer Vision and Image Understanding, 166:1–27, 2018.
-  Yong-Jin Liu, Chengchi Yu, Minjing Yu, and Ying He. Manifold SLIC: a fast method to compute content-sensitive superpixels. In IEEE Conference on Computer Vision and Pattern Recognition, CVPR ’16, pages 651–659, 2016.
-  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.
-  Yong-Jin Liu, Chun-Xu Xu, Ran Yi, Dian Fan, and Ying He. Manifold differential evolution (MDE): A global optimization method for geodesic centroidal voronoi tessellations on meshes. ACM Transactions on Graphics (SIGGRAPH ASIA 2016), 35(6):243:1–243:10, 2016.
-  Yong-Jin Liu, Dian Fan, Chun-Xu Xu, and Ying He. Constructing intrinsic delaunay triangulations from the dual of geodesic voronoi diagrams. ACM Trans. Graph., 36(2):15:1–15:15, 2017.
-  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.
-  Yong-Jin Liu. Semi-continuity of skeletons in 2-manifold and discrete voronoi approximation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 37(9):1938–1944, 2015.
-  Herbert Edelsbrunner and Nimish R. Shah. Triangulating topological spaces. International Journal of Computational Geometry and Applications, 7(4):365–378, 1997.
-  Lin Lu, Feng Sun, Hao Pan, and Wenping Wang. Global optimization of centroidal voronoi tessellation with monte carlo approach. IEEE Trans. Vis. Comput. Graph., 18(11):1880–1890, 2012.
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.
-  B B. H. Korte and Jens Vygen. Combinatorial Optimization: Theory and Algorithms. Spinger, 2012.
-  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.
-  Xavier Pennec. Intrinsic statistics on riemannian manifolds: Basic tools for geometric measurements. Journal of Mathematical Imaging and Vision, 25(1):127–154, 2006.
-  Raif M. Rustamov. Barycentric coordinates on surfaces. Comput. Graph. Forum, 29(5):1507–1516, 2010.
-  Vin de Silva and Joshua B. Tenenbaum. Global versus local methods in nonlinear dimensionality reduction. In Advances in Neural Information Processing Systems (NIPS ’02), pages 705–712, 2002.