Yuen-Shan Leung,Xiaoning Wang,Ying He,Yong-Jin Liu,and Charlie C.L.Wang

A uni fied framework for isotropic meshing based on narrowband Euclidean distance transformation

Yuen-Shan Leung1,Xiaoning Wang1,Ying He1,Yong-Jin Liu2,and Charlie C.L.Wang3

In this paper,we propose a simpleyet-e ff ective method for isotropic meshing relying on Euclidean distance transformation based centroidal Voronoi tessellation(CVT).Our approach improves the performance and robustness of computing CVT on curved domains while simultaneously providing highquality output meshes.While conventional extrinsic methods compute CVTs in the entire volume bounded by the input model,we restrict the computation to a 3D shell of user-controlled thickness. Taking voxels which contain surface samples as sites,we compute the exact Euclidean distance transform on the GPU.Our algorithm is parallel and memory-efficient, and can construct the shell space for resolutions up to 20483at interactive speed. The 3D centroidal Voronoi tessellation and restricted Voronoi diagrams are also computed efficiently on the GPU.Since the shell space can bridge holes and gaps smaller than a certain tolerance,and tolerate non-manifold edges and degenerate triangles,our algorithm can handle models with such defects,which typically cause conventional remeshing methods to fail.Our method can process implicit surfaces,polyhedral surfaces,and point clouds in a uni fied framework.Computational results show that our GPU-based isotropic meshing algorithm produces results comparable to state-ofthe-art techniques,but is signi fi cantly faster thanconventional CPU-based implementations.

centroidal Voronoi tessellation(CVT); Euclidean distance transformation; GPU;isotropic meshing; polygonal meshes;point clouds;implicit surfaces

1 Introduction

Triangle meshes have found widespread acceptance in computer graphics as a simple,convenient, and versatile representation of surfaces.However, raw meshes obtained from 3D scanners are often unsuitable for subsequent geometric processing,since they may contain holes,gaps,noise,degenerate triangles,and non-manifold edges.

A popular approach to improve mesh quality is to use centroidal Voronoi tessellation(CVT), which can generate a highly regular distribution of sites with respect to a given density function. A typical CVT-based remeshing method iteratively updates the generator of each Voronoi cell until it coincides with its center of mass.An isotropic mesh is then obtained as the dual graph of the computed CVT.A key step in each iteration in CVT computation is to construct a Voronoi diagram (VD).Although this is fairly simple to do in Euclidean space,doing so in curved domains is expensive due to the lack of a closed-form formula for geodesic distance.To tackle this challenge,some researchers parameterize the input models in ?2and computed a 2D CVT whose density function compensates for the parameterization distortion. These parameterization based methods can compute intrinsic CVTs on simple meshes of good quality, but they do not work for point clouds,imperfectmeshes,ormesheswith complicated topology, where parameterization is non-trivial.A practical alternative is to compute a restricted Voronoi diagram(RVD)[1],which is the intersection of the given model and a CVT de fi ned in ?3.RVD methods are easy to implement and can easily handle models with complicated geometry and topology.Although efficient clipping algorithms(e.g.,Ref.[2])have been proposed,their performance is still too low for time-critical applications.Moreover,current RVD methods assume the input is a manifold mesh,and they do not work for other surface representations, such as scanned point samples.

In this paper,we propose a new RVD-based computational framework for isotropic meshing, with the goal of improving performance as well as robustness of computing RVD while simultaneously maintaining the high quality of the output meshes. Rather than computing CVTs in the entire volume bounded by the input model,our idea is to restrict the computation to a 3D shell of user-controlled thickness.Taking voxelswhich contain surface samples as sites,we compute the exact Euclidean distance transform on the GPU.Our algorithm is parallel and memory-efficient,and can construct a shell space with resolution 20483at interactive speed. The 3D centroidal Voronoi tessellation and restricted Voronoi diagrams are also computed efficiently on the GPU.Computational results show that our GPU-enabled isotropic meshing algorithm produces results comparable to state-of-the-art techniques, but runs signi fi cantly faster than conventional CPU-based implementations.It is also worth noting that our method can process various representations, such as implicit surfaces,polyhedral surfaces,and point clouds,in a uni fied framework.Moreover, since the shell space can bridge holes and gaps smaller than a certain tolerance,and also tolerate non-manifold edges and degeneracies,our algorithm works well on imperfect meshes with such defects,for which conventionalremeshingmethods often fail.See Fig.1.

This paper makes the following contributions:

·We give auni fiedframework for constructing isotropic meshes from point clouds,polyhedral surfaces, and implicit surfaces via voxel representation. It avoids computationally expensive components often used in existing methods, such as data fi tting, isosurface extraction,and geodesic distance computation.

·Our framework produces a topologically consistent shell under user control,so it can bridge holes and gaps,and tolerate noise,to some degree.It also works well for models with non-manifold edges and degenerate triangles.

·We give a fast,memory-efficient algorithm for computing narrow-band distance fields on a GPU.Our implementation on an nVidia Quadro K5000 with 4GB RAM can compute theexactEuclidean distance transform at interactive speed.The CVT,RVD,and the dual Delaunay triangulations are also computed in parallel on the GPU.

Fig.1 Isotropic meshing on an imperfect mesh with non-manifold edges,degenerate triangles and holes.

2 Related work

2.1 Geodesic Voronoi diagrams

Although Voronoi diagrams in Euclidean space have been well studied,Voronoi diagrams in 2-manifold meshes(geodesicVoronoi diagrams,GVDs)have received little attention. Liu et al.[3]pointed out the fundamental di ff erence between conventional Voronoi diagrams and GVDs,and they pioneered a practical algorithm for constructing a GVD on a triangle mesh.Given a triangle meshMwithnvertices andm(≤n)sites onM,their algorithm has a theoretical time complexityO(n2logn)and runs empirically in linear time whenmis close ton.Liu and Tang[4]proved that the combinatorial complexity of a GVD isO(m(g+n)),wheregis the genus ofM.Based on a local Voronoi diagram,Xu et al.[5]developed an exact algorithm for computing a polyline-sourced GVD on a triangle mesh.

2.2 Centroidal Voronoi tessellations

A centroidalVoronoitessellation isaVoronoi diagram whose generating points are the centers of mass of the corresponding Voronoi cells.Thanks to its many favorable geometric properties,the CVT has been used in a wide range of applications[6]. It can be de fi ned by the critical points of the CVT energy function.A popular way to compute a CVT is Lloyd’s algorithm[7],which iteratively moves each generator to the corresponding mass center until convergence.Lloyd’s method is conceptually simple and easy to implement,but as a gradient-descent method,it only has linear convergence. Liu et al.[8]proved that the CVT energy function has 2nd order smoothness in most situations encountered in computer graphics,so one can achieve quadratic or super-linear convergence by use of Newton or quasi-Newton optimization methods.Rong et al.[9] developed a GPU-based method for computing a CVT in the plane and observed signi fi cant speedup over CPU methods.

To compute the CVT on a genus-0 surface,Alliez et al.[10]conformally parameterized the surface to a disk and evaluated the centroids over a weighted density function expressed in ?2.Rong et al.[11] extended the CVT energy function to the spherical S2and hyperbolic H2domains,and thus can process surfaces of all topological types.Shuai et al.[12] developed a GPU-enabled algorithm for computing a periodic CVT in H2and used it to construct isotropic meshes for high-genus surfaces.Recently,Wang et al.[13]proposed an intrinsic method for computing the CVT on an arbitrary manifold triangle mesh. Rather than computing the mass centers of Voronoi cells which requires area integration,their algorithm computes the Riemannian centers using exponential maps.These Riemannian centers provide a good approximation of the mass centers if there are sufficient seeds.

Theseparameterization-based and exponential map based CVT algorithms are intrinsic and thereby independent of the embedding space.However,they are computationally expensive and impractical for time-critical applications.Rather than computing CVT on surfaces directly,Yan et al.[1]proposed a novelindirectmethod based on computing arestricted Voronoi diagram,that is,the intersection between the input mesh and a 3D CVT.Using akdtree to quickly identify and compute the intersection of each triangle face with its incident Voronoi cells, their algorithm computes the RVD inO(nlogm) time,wheremis the number of seed points andnis the number of triangles of the input mesh.They also adopted a quasi-Newton method to efficiently compute the 3D CVT in the volume bounded by the input mesh.

Lu et al.[14]computed a CVT using line segments and graphs as generators.CVT can also be de fi ned using anisotropic metrics[15]andLpdistances[16].

Chen et al.[17]proposed an iterative method for generating a constrained centroidal Delaunay mesh(CCDM).With local vertex relation,their method does not require geodesic distances and can guarantee that the computed CCDM has the same topology as the input mesh.However,their method cannot handle non-manifold edges and degenerate triangles.It is also unclear whether their method can be extended to point clouds and implicit functions.

Li et al.[18]presented an elegant method for triangulating the conformal uniformization domain via planar Delaunay re fi nment.They gave explicit estimates for the Hausdor ff distance, normal deviation,and di ff erences in curvature measures between the surface and the mesh.

2.3 Distance computation

Fig.2 Overview of our approach using theSculpturemodel.(a)Input mesh;(b)shell space withd=3;(c)CVT with 3000 seeds; (d)output isotropic mesh.

A key step when intrinsically computing a CVT on a polyhedral surface is to determine geodesic distances,a fundamental problem in computer graphicsand computationalgeometry.Classical computational geometry approaches include the MMP algorithm[19],the CH algorithm[20],and their many variants(e.g.,Refs.[21,22]).Although these methods can compute the exact solution on arbitrary manifold meshes,they are computationally expensive.Partialdi ff erentialequation (PDE) approaches,such as the fast marching method[23] and the heat method[24],are efficient and can fl exibly handle a wide range of geometric domains, including triangle meshes,point clouds,regular grids and volumes.They can also be easily generalized to an anisotropic setting[25].However,they only compute a first-order approximation,and thus the results are sensitive to model tessellation and resolution.Motivated by thelocalstructure of the discrete geodesic problem,Ying et al.[26]proposed the saddle vertex graph(SVG),a sparse graph which encodes geodesic information on a triangle mesh.Computing a geodesic distance is equivalent to fi nding a shortest path in this graph.However, the performance of SVG methods depends greatly on the number and distribution of saddle vertices in the input meshes,and it is unclear whether it can be extended to imperfect meshes,implicit surfaces, or point clouds,as it is difficult to determine saddle vertices in such cases.

Since computing geodesic distances is expensive, many applicationsseek practicalwaysto give approximate solutions with accuracy control.The most widely used approach is to use a distance field,which encodes the distance from each grid node to the surface of the model(uniformly or adaptively)[27].An approximate distance field can be efficiently computed by a distance-transformation method, e.g., based on thechamferdistance transformation(CDT)[28]or thevector-city vector distance transform(VCVDT)[29].

Recent work has focused on theEuclidean distance transformation(EDT).Given ad-dimensional grid withN=ndgrid points,whereSgrid points are colored and the otherN?Sgrid points are not colored,for each grid point,the EDT aims to compute the Euclidean distance to the closest colored grid point.However,this process is typically computationally intensive in 3 or more dimensions. GPU-based algorithms have been developed to provide an efficient solution[30]. Cao et al.[31] presented an efficient algorithm,theparallel banding algorithm(PBA),to compute the exact EDT on a GPU.Their algorithm,however,can quickly exhaust graphics memory:the highest resolution that can be achieved by PBA is 5123on a cutting-edge graphics card,which does not provide sufficient accuracy to guarantee the correctness of the constructed CVT for many real-world models.The new algorithm proposed in this paper is more sophisticated inits memory usage and can achieve 3D EDT for resolutions up to 20483.

2.4 Surface reconstruction and meshing from points

Constructing a high-quality mesh from unorganized point samples plays an important role in computer graphics. Although numerous techniques are available,all tackle the problem in twoseparatestages,that is,constructing an initial surface first, and then improving the mesh quality using the remeshing techniques.Existing surface reconstruction methods are based on either computational geometry or implicit functions.

Computationalgeometry approachesaim to generate a piecewise linear surface,topologically equivalent to the sampled surface and also geometrically close.Representative work includescrust[32],cocone[33],and their many variants. These algorithms provide theoretical guarantees if the input points are densely sampled.However,in practice,this condition may not hold.

In presenceofnoiseand outliers,which is typical for samples obtained by scanning,a common approach is to fi t the samples using the zero set of an implicit function.A popular method fororientedpointsamplesisPoissonsurface reconstruction[34]; an implicit surface is first generated and tessellated into a polygonal mesh later.Although tessellation step can be computed on the GPU, fi tting an implicit surface uses timeconsuming global operators.Sheung and Wang[35] presented a robust mesh reconstruction method forunorientednoisy points.An approximate CVT is computed on the point set to down-sample the input point cloud to a smaller number of pointa. Then,a Voronoi diagram based mesh generation method,tight cocone[36],is employed to generate the mesh.As VD-based surface reconstruction is memory intensive,only datasets with a modest number of points can be handled.Moreover,it is hard to parallelize the computation.

Our proposed method is di ff erent from existing ones in that it performs surface reconstruction and meshing in an integrated manner.Our method ismoreefficientand elegantsinceit can completely bypass implicit fi tting as well as isosurface extraction,making it an ideal tool for processing scanned point samples.

3 Algorithm

LetOdenote the input 3D mesh.We first construct a voxel representationMofOat a given resolutionres.Then we construct a shell spaceˉPconsisting of o ff-surface points,where each point inˉPhas a distancedp≤dto its closest point onM(see Fig.2).The thresholddis speci fied by the user and is model-dependent.

Our isotropic meshing algorithm adopts Lloyd’s framework. Starting withkrandomly generated seeds,it minimizes the CVT energy by iteratively updating the seed positions.In each iteration,it computes a Voronoi diagram con fi ned to the shell space and moves the seeds toward the corresponding mass centers.The algorithm projects the seeds back toˉPif they move outside the shell space.After convergence,it propagates the seed information in the shell space to determine connected Voronoi cells and extracts the dual Delaunay triangulation.Our method is outlined in Algorithm 1,and described in detail in the following.

Algorithm 1:Isotropic meshing based on EDT Input:3D surfaceO,voxel resolutionres,shell space thicknessd,convergence threshold∈,and number of seedsk. Output:Isotropic mesh withkvertices. 1:S←krandom seeds 2:M←Voxelization(O,res)3:ˉP←ShellSpaceConstruction(M,res,d) 4:while convergence not reached do 5:Vk←SearchClosestSeedInShell(ˉP,S) 6:Ck←CenterMass(Vk)7:ˉCk←UpdateSeed(ˉP,Ck) 8:S←ˉC9:end while 10:Vk←ShellFlooding(ˉP,S) 11:return DualTriangulation(Vk)

3.1 Memory-efficient shell space construction

We introduce a memory-efficient way to construct shell spaces in real-time. We extend the parallel banding algorithm(PBA)of Cao et al.[31]to compute the Euclidean distance transform(EDT)in a narrow-band.Their algorithm partitions the input domain into small chunks of equal size,which can be processed in parallel.The results are then mergedconcurrently.Although their method is exact and efficient,it is impractical for large models due to its memory requirements.Our method overcomes this problem by on-the- fl y computation and integrating fast bitmap indexing technology on the GPU.We explain the principle in two dimensions for simplicity; the idea easily extends to three dimensions.

The input image is divided into a virtual grid made up of occupied and unoccupied pixels;the former pixels,thesites∈S,are run-length encoded. Every pixel undergoes a two-step process:(i) fi nd the nearest siteSij,among all sites in rowj;(ii) determine the closest site,among all the nearest sites,in the current columni.During the first step we assign a thread to process each row as it is more efficient to do more computation in a single thread than to repeatedly access global memory with multiple threads.The second step extends the divide-and-merge approach of PBA,by employing warp-vote and warp-shuffle functions in NVIDIA’s CUDA to exchange information between nearest sites within a chunk.(A warp is a pool of threads that executes in parallel.) Every thread in the same warp does the same calculation,avoiding warp variations that often compromise performance.

Figure 3 shows an example.When the threads in a warp reach columni,each row computes the corresponding nearest sitesSij.The nearest site of the current pixel is shown in blue.Note that only some sites(empty blue dots)satisfy the distance constraint.Therefore,sitesSi1of thread 1 can be safely discarded.We set a barrier to ensure that every thread obtains some sites before exchanging information.After synchronization,we sweep eachSijto other threads in the same warp to update the closest site of the current pixel(i,j),based on the distance functiondij. A bitmap stores a Boolean value for every pixel;0 indicates that thecorresponding nearest sitescannotbe the closest sites of the current column.Then the threads repeat the same procedure for the next columni+1,until they reach the last column.In the final step we collect the closest sites with value 1 in di ff erent chunks and merge them to fi nd the pixels that form the shell region. As the nearest sites are computed on-the- fl y and the temporary result only needs indexing by a bitmap,our algorithm requires less memory than the PBA method.

Fig.3 Illustrative example of distance field computation in a narrow band. See the text for an explanation and functionShellSpaceConstructionfor the pseudo code.

3.2 Constructing 3D Voronoi diagrams in shell space

The shell space represented byˉPis used to constrain construction of the Voronoi diagram and update the positions of seeds.Initially,the seedsS={s1,...,sk}are located on the input mesh.We collect points∈ˉPthat share the same closest seeds to build the Voronoi diagram.

We perform a proximity search to fi nd the closest seed for all query points from ˉP.To speed up the process,the seeds are projected onto a uniform gridGwith smaller resolution ofres(e.g.,323),so that, for a query pointq,it just fi nd seeds in the grid cellGqin whichqfalls.If any border of the grid cell is closer toqthan the seed found inGq,query pointqconsiders the neighboring cell ofGq.

3.3 Computing the CVT

Updating the seeds’positions towards a uniform distribution is crucial for constructing the CVT. Optimal positions can be reached by minimizing the following energy function:

whereViis the Voronoi cell of seedsi,p∈ˉP,andρis a non-negative user-de fi ned density function.

In Lloyd’s algorithm,a seedsimoves iteratively toward the corresponding mass centerciof Voronoi diagramViuntil convergence.However,the mass center could be located far from the surface,as it is constructed in the shell space.We determine the following new position to replace the mass center in each iteration:

whereu∈?+is the magnitude of movement of the seeds.We observe that if the seeds move with di ff erent magnitudes,the area of CVT calls will largely vary depending on surface curvature.In order to guarantee topological consistency,the new center is projected back ontoMif it moves outside the shell space,as shown in Fig.4.

3.4 Computing the dual triangulations

Upon convergence,all the generators are uniformly distributed. We next extract the dual Delaunay triangulation.

First,we fi nd thedirectneighbors of all seeds, where direct means two voxels from their Voronoi cells are connected.Adjacent neighbors can be found by fl ooding the information from all seeds to all voxels in the shell∈ˉP.Each voxel stores a hash table to hold the locations of its 26 neighbors.Each propagation step updatesthe current seed information to neighboring voxels until all voxels have been reached.This approach avoids producing the wrong connectivity for seeds that are geometrically close,but topologically far from each other.We then organize direct neighbors in clockwiseorder and finally extract the triangle mesh.

Fig.4 Illustration of the update process in an iteration for two seeds.(a)Yellow dots:seeds of Voronoi cellsVi(in red)andVj(in green).Red and green dots:their respective mass centers.(b)The seeds move along vector?→scto new centers(blue dots).We projectˉcj(light blue dot)to the surface since it is outside the shell region.

4 Experimental results

All tests were performed on a PC with an Intel Xeon E5 2.5GHz CPU and an nVidia Quadro K5000 GPU with 4GB RAM.

4.1 Narrow-banded distance fields

Table 1 lists the time and peak memory requirements for varying parameterdandres(see Fig.5).Clearly, asdincreases,the time increases insigni fi cantly compared to the increased number of nearest sites in the shell.This is due to the low cost of intra-warp communication and the reduced warp divergence in our algorithm.Also,the memory consumption is remarkably small,considering the scene is represented by a high resolution uniform grid.Previous algorithms(e.g.,Ref.[31])typically require at least 10 times more memory than ours. Table 2 compares our algorithm with PBA on the Sculpture model at resolution 5123.Since PBA computes a full distance map,its performance is independent of distanced.Our algorithm consumes signi fi cantly less memory and runs much faster than PBA for a reasonably smalld.

Table 1 Performance of our algorithm for di ff erentdat resolutions 10243and 20483

Table 2 Comparison of our algorithm to PBA at resolution 5123

4.2 CVT computation

Following Ref.[13],we adopted the following criteria to measure the triangle mesh quality:(i)triangle qualityQ(t)de fi ned by,whereandare the inradius and the length of the longest edge of triangletrespectively,(ii)the smallest angleθminand the averageθavgminimal angle for all triangles,(iii)the ratio of singularities,de fi ned byvs/k,wherevsis the number of non-6-valent vertices andkis the number of vertices.

Fig.5 Mesh quality for various shell space parameter valuesd.(a)Triangulation quality measure.(b)Singularity ratio.(c)Our GPU-based Lloyd algorithm usually converges in 100-200 iterations;dhas little impact on the convergence rate.Horizontal axis: iteration number.Vertical axis:normalized CVT energy function.

Flexible inputs.One advantage of our algorithm is that the CVT construction is independent of input topology,thus allowing a wide range of input types such as point clouds,implicit surfaces,and meshes with non-manifold surfaces. Figure 6 shows that our method can be applied to an irregular point cloud.However,the uneven density of points causes failure during triangulation,i.e.,holes and a nonmanifold mesh.The shell space can help to resolve this issue because it covers the uneven space where points cannot be reached.Figure 7 investigates the in fl uence of shell distance on irregular data sets.It shows that the geometric shape and density of the point cloud directly a ff ect the choice of the shell distanced.

Fig.6 Experimental results on scanned point samples.

We have also successfully applied our algorithm to implicit surfaces(see Fig.8)and other implicit representations(see Fig.9).More examples can be found in Fig.12.

Figure 1 shows the result from an imperfect mesh with non-manifold edges,degenerate triangles and holes.Thanks to its defect-tolerant nature,our method can handle these issues easily.

Accuracy and efficiency control.We allow the user to balance accuracy and efficiency by choice of the o ff setd.Figure 5 shows the relationship between the distanced,the number of generators, and the quality of the remeshed surface. As the o ff set increases to 9,with 4000 generators,the mesh quality of the dinosaur model dramatically drops. This also happens when the o ff set decreases to 2 with 1000 generators.Figure 5(b)clearly illustrates the di ff erence in quality for di ff erent values ofd.Along with other examples in Table 3,we can see that the mesh has best quality with o ff set distance in the range 2-6.Figure 10 compares our method with two parameterization-free isotropic meshing methods,anintrinsicCVT method[13]and anextrinsicRVD method[1].

Thanks to its GPU-friendly structure and the computational power of modern GPUs,our method runssigni fi cantly fasterthan theirCPU-based implementations.

Fig.7 Isotropic meshing on noisy scanned point samples. We observe thatd∈[3,5]helps to bridge holes and tolerate noise,producing fairly good results.However,ifdexceeds that range,there is a large error in the computed CVT,leading to poor results.

Table 3 Model complexity and runtime performance.SS:time(in seconds)for shell construction;m:the number of seeds;T: average time for each Lloyd iteration;n:the number of iterations;I:singularity ratio

Fig.8 Generating isotropic meshes from an implicit function. There are 3000 seeds andd=3.

Fig.9 Generating isotropic meshes with 8000 seeds andd=2 from an implicit representation: layered depth normal image[37].

Fig.10 Comparison with the RVD method[1]and the intrinsic CVT method[13].

Topology consistency. The shell space also guarantees topological consistency between the input and the remeshed surface.Our algorithm constrains the seeds to lie within a shell of width 2dduring the update process,so the output is topologically consistent with respect to the shell. Figure 11 illustrates the importance of this o ff set volume.The bottle model consists of an inner tube which is very close to the outer surface,close enough that only a shell of width 2(d=1)can separate these two parts when the resolution is 10243.However,width 2 is so thin that the seeds can barely move inside,leading to poor remeshing quality(Qavg=0.86,Qmin=0.10).A possible solution would be to increase the grid size.Since our approach is scalable,this could be improved by using a graphic card with larger memory.

Fig.11 Remeshing the genus-2 Knotty Bottle model with 9000 seeds.As its inner tube is very close to the outer surface,d=1 is needed(at resolution 10243)to produce a shell space with the same topology as the mesh,whereasd=2 does not give such a guarantee.

Fig.12 Experimental results.Images are rendered at highresolution,allowing zooming in for examination.

5 Conclusions

This paper has presented a robust, efficient method for constructing isotropic meshes using the Euclidean distance transform.Our algorithm constructsa narrow band enclosing theinput surface,in which 3D centroidal Voronoi tessellations and restricted Voronoi diagrams are computed. Ouralgorithm isfully paralleland memoryefficient,and can construct the shell space with resolution up to 20483at interactive speed. Moreover,our method can process implicit surfaces, polyhedral surfaces and point clouds in a uni fied framework.Computational results show that our GPU-friendly isotropic meshing algorithm produces results comparable to state-of-the-art techniques, but runs signi fi cantly faster than conventional CPU-based implementations.

Our current implementation adopts a constant resolution to construct the shellspace.This, however,is not optimal,since it is pessimistic for the regions with fairly fl at geometry,and it may be inadequate for the highly-curved regions.In the future,we aim to develop a geometry-aware algorithm for constructing in parallel the shell space with adaptive resolution.

Acknowledgements

The NTU authors are partially supported by AcRF RG40/12 and MOE2013-T2-2-011.Y.-J.Liu is partially supported by National Natural Science Foundation of China(Nos.61432003 and 61322206) and the TNList Cross-discipline Foundation.C.C. L.Wang is partially supported by HKSAR Research Grants Council(RGC)General Research Fund (GRF),CUHK/14207414.

Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use,distribution,and reproduction in any medium,provided the original author(s)and the source are credited.

[1]Yan,D.-M.;L′evy,B.;Liu,Y.;Sun,F.;Wang,W. Isotropic remeshing with fast and exact computation of restricted Voronoi diagram.Computer Graphics ForumVol.28,No.5,1445-1454,2009.

[2]Yan,D.;Wang,W.;L′evy,B.;Liu,Y.Efficient computation of clipped Voronoi diagram for mesh generation.Computer-Aided DesignVol.45,No.4, 843-852,2013.

[3]Liu,Y.J.;Chen,Z.;Tang,K.Construction of iso-contours,bisectors,and Voronoi diagrams on triangulated surfaces.IEEE Transactions on PatternAnalysis and Machine IntelligenceVol.33,No.8, 1502-1517,2011.

[4]Liu,Y.-J.;Tang,K.The complexity of geodesic Voronoi diagrams on triangulated 2-manifold surfaces.Information Processing LettersVol.113,No.4,132-136,2013.

[5]Xu,C.;Liu,Y.-J.;Sun,Q.;Li,J.;He,Y.Polylinesourced geodesic Voronoi diagrams on triangle meshes.Computer Graphics ForumVol.33,No.7,161-170, 2014.

[6]Du,Q.;Gunzburger,M.D.;Ju,L.Constrained centroidal Voronoi tessellations for surfaces.SIAM Journal on Scienti fi c ComputingVol.24,No.5,1488-1506,2002.

[7]Lloyd,S.Least squares quantization in PCM.IEEE Transactions on Information TheoryVol.28,No.2, 129-137,1982.

[8]Liu,Y.;Wang,W.;L′evy,B.;Sun,F.;Yan,D.-M.; Lu,L.;Yang,C.On centroidal Voronoi tessellation—Energysmoothnessand fastcomputation.ACM Transactions on GraphicsVol.28,No.4,Article No. 101,2009.

[9]Rong,G.;Liu,Y.;Wang,W.;Yin,X.;Gu,X.D.;Guo, X.GPU-assisted computation of centroidal Voronoi tessellation.IEEE Transactions on Visualization and Computer GraphicsVol.17,No.3,345-356,2011.

[10]Alliez,P.;de Verdi`ere,′E.C.;Devillers,O.;Isenburg, M.Isotropic surface remeshing.In:Proceedings of the Shape Modeling International,49,2003.

[11]Rong,G.;Jin,M.;Shuai,L.;Guo,X.Centroidal Voronoi tessellation in universal covering space of manifold surfaces.Computer Aided Geometric DesignVol.28,No.8,475-496,2011.

[12]Shuai,L.;Guo,X.;Jin,M.GPU-based computation of discrete periodic centroidal Voronoi tessellation in hyperbolic space.Computer-Aided DesignVol.45,No. 2,463-472,2013.

[13]Wang,X.;Ying,X.;Liu,Y.-J.;Xin,S.-Q.;Wang, W.;Gu,X.;Mueller-Wittig,W.;He,Y.Intrinsic computation of centroidal Voronoi tessellation(CVT) on meshes.Computer-Aided DesignVol.58,51-61, 2015.

[14]Lu,L.;L′evy,B.;Wang,W.Centroidal Voronoi tessellation of line segments and graphs.Computer Graphics ForumVol.31,No.2,775-784,2012.

[15]L′evy,B.;Bonneel,N.Variational anisotropic surface meshing with Voronoi parallel linear enumeration. In:Proceedings of the 21st International Meshing Roundtable,349-366,2013.

[16]L′evy,B.;Liu,Y.Lpcentroidal Voronoi tessellation and its applications.ACM Transactions on GraphicsVol.29,No.4,Article No.119,2010.

[17]Chen,Z.;Cao,J.;Wang,W.Isotropic surface remeshing using constrained centroidalDelaunay mesh.Computer Graphics ForumVol.31,No.7,2077-2085,2012.

[18]Li,H.;Zeng,W.;Morvan,J.M.;Chen,L.;Gu,X. Surface meshing with curvature convergence.IEEETransactions on Visualization and Computer GraphicsVol.20,No.6,919-934,2014.

[19]Mitchell,J.S.B.;Mount,D.M.;Papadimitriou,C. H.The discrete geodesic problem.SIAM Journal on ComputingVol.16,No.4,647-668,1987.

[20]Chen,J.;Han,Y.Shortest paths on a polyhedron. In:Proceedings of the 6th Annual Symposium on Computational Geometry,360-369,1990.

[21]Xu,C.;Wang,T.;Liu,Y.;Liu,L.;He,Y.Fast wavefront propagation(FWP)for computing exact discrete geodesics on meshes.IEEE Transactions on Visualization and Computer GraphicsVol.21,No.7, 822-834,2015.

[22]Ying,X.;Xin,S.-Q.;He,Y.Parallel chen-han(PCH) algorithm for discrete geodesics.ACM Transactions on GraphicsVol.33,No.1,Article No.9,2014.

[23]Kimmel,R.;Sethian,J.A.Computing geodesic paths on manifolds.Proceedings of the National Academy of Sciences of the United States of AmericaVol.95,8431-8435,1998.

[24]Crane,K.;Weischedel,C.;Wardetzky,M.Geodesics in heat:A new approach to computing distance based on heat fl ow.ACM Transactions on GraphicsVol.32, No.5,Article No.152,2013.

[25]Campen,M.;Heistermann,M.;Kobbelt,L.Practical anisotropic geodesy.Computer Graphics ForumVol. 32,No.5,63-71,2013.

[26]Ying,X.;Wang,X.;He,Y.Saddle vertex graph (SVG):A novel solution to the discrete geodesic problem.ACM Transactions on GraphicsVol.32,No. 6,Article No.170,2013.

[27]Jones,M.W.;Baerentzen,J.A.;Sramek,M. 3D distance fields: A survey of techniques and applications.IEEE Transactions on Visualization and Computer GraphicsVol.12,No.4,581-599,2006.

[28]Marchand-Maillet,S.;Sharaiha,Y.M.Euclidean ordering via chamfer distance calculations.Computer Vision and Image UnderstandingVol.73,No.3,404-413,1999.

[29]Satherley,R.;Jones,M.W.Vector-city vector distance transform.Computer Vision and Image UnderstandingVol.82,No.3,238-254,2001.

[30]Ho ffIII,K.E.;Keyser,J.;Lin,M.;Manocha,D.; Culver,T.Fast computation of generalized Voronoi diagrams using graphics hardware.In:Proceedings of the 26th Annual Conference on Computer Graphics and Interactive Techniques,277-286,1999.

[31]Cao,T.-T.;Tang,K.;Mohamed,A.;Tan,T.-S. Parallel banding algorithm to compute exact distance transform with the GPU.In:Proceedings of the 2010 ACM SIGGRAPH Symposium on Interactive 3D Graphics and Games,83-90,2010.

[32]Amenta,N.;Bern,M.;Kamvysselis,M.A new Voronoi-based surface reconstruction algorithm.In: Proceedingsofthe25th AnnualConferenceon Computer Graphics and Interactive Techniques,415-421,1998.

[33]Amenta,N.;Choi,S.;Dey,T.K.;Leekha,N.A simple algorithm for homeomorphic surface reconstruction. In:Proceedings of the 16th Annual Symposium on Computational Geometry,213-222,2000.

[34]Kazhdan,M.M.;Bolitho,M.;Hoppe,H.Poisson surface reconstruction.In: Proceedings of the 4th Eurographics Symposium on Geometry Processing, 61-70,2006.

[35]Sheung, H.; Wang, C.C.L.Robustmesh reconstruction from unoriented noisy points.In: ProceedingsofSIAM/ACMJointConferenceon Geometric and Physical Modeling,13-24,2009.

[36]Dey,T.K.;Goswami,S.Tight cocone:A water-tight surface reconstructor.In:Proceedings of the 8th ACM Symposium on Solid Modeling and Applications,127-134,2003.

[37]Wang,C.C.L.;Leung,Y.-S.;Chen,Y.Solid modeling of polyhedral objects by layered depth-normal images on the GPU.Computer-Aided DesignVol.42,No.6, 535-544,2010.

Yuen-Shan Leung received her B.Eng.degree and her Ph.D.degree from the Chinese University of Hong Kong.She is currently a research fellow at the School of Computer Engineering, Nanyang Technological University, Singapore. Her research interests include geometric& solid modeling, computer-aided design&manufacturing,and computer graphics.

Xiaoning Wang received the B.Eng. degree from Tianjin University,China, in 2011.He iscurrently a Ph.D. candidate in the School of Computer Engineering, Nanyang Technological Univeristy, Singapore. His research interestsare geometry analysisand computer graphics.

Ying He is currently an associate professor at the School of Computer Engineering, Nanyang Technological University,Singapore.He received the B.S.and M.S.degreesin electrical engineering from Tsinghua University, China, and the Ph.D.degree in computer science from Stony Brook University,USA.Hisresearch interestsfallinto the general area of visual computing,and he is particularly interested in the problems which require geometric analysis and computation. For more information,please visit http://www.ntu.edu.sg/home/yhe.

Yong-Jin Liu received his B.Eng. degree from Tianjin University,China, in 1998,and his M.Phil.and Ph.D. degrees from Hong Kong University of Science and Technology(HKUST), Hong Kong,China, in 2000 and 2004,respectively. Heisnow an associateprofessoratthe Tsinghua National Laboratory for Information Science and Technology, the Department of Computer Science and Technology,Tsinghua University, China. His research interests include computational geometry, multimedia, computer graphics, and computeraided design. For more information, please visit http://cg.cs.tsinghua.edu.cn/people/∽Yongjin/yongjin.htm.

Charlie C.L.Wang is a Fellow of the American Society of Mechanical Engineers(ASME)with expertise in geometric computing, design, and manufacturing.He received his Ph.D. in mechanical engineering from HKUST in 2002. After that,he joined the Chinese University ofHong Kong in 2003,and is now a full professor of mechanical and automation engineering. His research interests include geometric computing,computer-aided design,advanced manufacturing,and computational physics.

Other papers from this open access journal are available free of charge from http://www.springer.com/journal/41095. To submit a manuscript,please go to https://www. editorialmanager.com/cvmj.

1School of Computer Engineering,Nanyang Technological University,Singapore.E-mail:Y.-S.Leung,ysleung@ ntu.edu.sg;X.Wang,wang0684@e.ntu.edu.sg;Y.He, yhe@ntu.edu.sg.

2DepartmentofComputerScienceand Technology, Tsinghua University, China.E-mail:liuyongjin@ tsinghua.edu.cn.

3Department of Mechanical and Automation Engineering, The Chinese University of Hong Kong,Hong Kong, China.E-mail:cwang@mae.cuhk.edu.hk.

Manuscript received:2015-08-17;accepted:2015-08-29