3D Snap Rounding

: Let P be a set of n polygons in R 3 , each of constant complexity and with pairwise disjoint interiors. We propose a rounding algorithm that maps P to a simplicial complex Q whose vertices have integer coordinates. Every face of P is mapped to a set of faces (or edges or vertices) of Q and the mapping from P to Q can be done through a continuous motion of the faces such that (i) the L ∞ Hausdorﬀ distance between a face and its image during the motion is at most 3/2 and (ii) if two points become equal during the motion, they remain equal through the rest of the motion. In the worse the size of Q is O ( n 15 ) and the time complexity of the algorithm is O ( n 19 ) but, under reasonable hypotheses, these complexities decrease to O ( n 5 ) and O ( n 6 √ n ) .


Introduction
Rounding 3D polygonal structures is a fundamental problem in computational geometry.Indeed, many implementations dealing with 3D polygonal objects, in academia and industry, require as input pairwise-disjoint polygons whose vertices have coordinates given with fixed-precision representations (usually with 32 or 64 bits).On the other hand, many algorithms and implementations dealing with 3D polygonal objects in computational geometry output polygons whose vertices have coordinates that have arbitrary-precision representations.For instance, when computing boolean operations on polyhedra, some new vertices are defined as the intersection of three faces and their exact coordinates are rational numbers whose numerators and denominators are defined with roughly three times the number of bits used for representing each input coordinate.When applying a rotation to a polyhedron, the new vertices have coordinates that involve trigonometric functions.When sampling algebraic surfaces, the vertices are obtained as solutions of algebraic systems and they may require arbitrary-precision representations since the distance between two solutions may be arbitrarily small (depending on the degree of the surface).
This discrepancy between the precision of the input and output of many geometric algorithms is an issue, especially in industry, because it often prevents the output of one algorithm from being directly used as the input to a subsequent algorithm.
In this context, there exists no solution for rounding the coordinates of 3D polygons with the constraint that their rounded images do not properly intersect and that every input polygon and its rounded image remain close to each other (in Hausdorff distance).In practice, coordinates are often rounded without guarding against changes in topology and there is no guarantee that the rounded faces do not properly intersect one another.
The same problem in 2D for segments, referred to as snap rounding, has been widely studied and admits practical and efficient solutions [1,[5][6][7][8][9][10][11]14].Given a set of possibly intersecting segments in 2D, the problem is to subdivide their arrangement and round the vertices so that no two disjoint segments map to segments that properly intersect.For clarity, all schemes consider that vertices are rounded on the integer grid.It is well known that rounding the endpoints of the edges of the arrangement to their closest integer point is not a good solution because it may map disjoint segments to properly intersecting segments.. Snap rounding schemes propose to further split the edges when they share a pixel (a unit square centered on the integer grid).In such schemes, disjoint edges may collapse but this is inevitable if the rounding precision is fixed and if we bound the Hausdorff distance between the edges and their rounded images.Furthermore, it is NP-hard to determine whether it is possible to round simple polygons with fixed precision and bounded Hausdorff distance, and without changing the topological structure [12].
In dimension three, results are extremely scarse, despite the significance of the problem.Goodrich et al. [5] proposed a scheme for rounding segments in 3D, and Milenkovic [13] sketched a scheme for polyhedral subdivisions but, as pointed out by Fortune [4], both schemes have the property that rounded edges can cross.Fortune [3] suggested a high-level rounding scheme for polyhedra but in a specific setting that does not generalize to polyhedral subdivisions [4].Finally, Fortune [4] proposed a rounding algorithm that maps a set P of n disjoint triangles in R 3 to a set Q of triangles with O(n 4 ) vertices on a discrete grid such that (i) every triangle of P is mapped to a set of triangles in Q at L ∞ Hausdorff distance at most 3  2 from the original face and (ii) the mapping preserves or collapses the vertical ordering of the faces.Unfortunately, this rounding scheme is very intricate and, moreover, it uses a grid precision that depends on the number n of triangles: the vertices coordinates are rounded to integer multiples of about 1  n .The difficulty of snap rounding faces in 3D is described by Fortune [4]: First, it is reasonable to round every vertex to the center of the voxel containing it (a voxel is a unit cube centered on the integer grid).But, by doing so, a vertex may traverse a face and to avoid that, it might be necessary to add beforehand a vertex on the face, which requires triangulating it.Newly formed edges may cross older edges when snapping; to avoid this, new vertices are added to these edges, in turn requiring further triangulating of faces.It is not known whether such schemes terminate.
To better understand the difficulty of the problem, consider the following simple but flawed algorithm.First project all the input faces onto the horizontal plane, subdivide the projected edges as in 2D snap rounding, triangulate the resulting arrangement, lift this triangulation vertically on all faces, and then round all vertices to the centers of their voxels.For an input of size n, this yields an output of size Θ(n 4 ) in the worst case and an L ∞ Hausdorff distance of at most 1 2 between the input faces and their rounded images.Unfortunately, this algorithm does not work in the sense that edges may cross: indeed, consider two almost vertical close triangles whose projections on the horizontal plane are disjoint triangles that are rounded in 2D to the same segment; such triangles in 3D may be rounded into properly overlapping triangles.Fortune [4] solved this problem by using a finer grid to round the vertices and avoiding the faces to round vertically.
Contributions.We present in this paper the first algorithm for rounding a set of interiordisjoint polygons into a simplicial complex whose vertices have integer coordinates and such that the geometry does not change too much: namely, (i) the Hausdorff distance between every input face and its rounded image is bounded by a constant ( 32 for the L ∞ metric) and (ii) the relative positions of the faces are preserved in the sense that there is a continuous motion that deforms all input faces into their rounded images such that if two points collapse at some time, they remain identical up to the end of the motion.This ensure, in particular, that if a line stabs two input faces far enough from their boundaries, the line will stab their rounded images in the same order or in the same point.
The worst-case complexity of our algorithm is polynomial but unsatisfying as our upper bound on the output simplicial complex is O(n 15 ) for an input of size n.However, this upper bound decreases to O(n 5 ) under some reasonable assumption on the input: roughly speaking that the input is a nice discretization of a constant number of surfaces that satisfy some reasonable hypothesis on their curvature.The corresponding time complexity reduces from O(n 19 ) to O(n 6 √ n).Furthermore, it very likely that these bounds are not tight and, in practice on realistic non-pathological data, we anticipate time and space complexities of O(n √ n).
We present the algorithm in Section 3, its proof of correctness in Section 4, and its complexity analysis in Section 5.

Preliminaries
Notation.The coordinates in the Euclidean space R 3 are referred to as x, y, and z and ı, , k is the canonical basis.We use several planes parallel to the axes to project or intersect some faces: the xy-plane is called the floor, the xz-plane is called the back wall and a plane parallel to the yz-plane is called a side wall.Projections on the floor and on the back wall are always considered orthogonal to the plane of projection.
Two polygons, edges, or vertices are said to properly intersect if their intersection is nonempty and not a common face of both.Two polygons (resp.segments) intersect transversally if their relative interiors intersect and if they are not coplanar (resp.collinear).

General position assumption.
For the sake of simplicity, we assume, without loss of generality, some general position on our input set of polygons P. Precisely, we assume: (α) No faces are parallel to the axes of coordinates and no vertices project in direction ±  on an edge (except the endpoints of that edge).

Inria
(β) No supporting plane of a face translated by ±  contains a vertex.
Let I denote the intersection, if not empty, of the supporting plane of a face with the translation by ±  of the supporting plane of another face.By assumption (β), I is a line.
(γ) No vertices project in direction ±  onto such a line I.
(δ) For any point A on a face and with half-integer x and y-coordinates, A ±  does not belong to another face.More generally, no line I crosses any vertical line defined by half-integer x and y-coordinates.
This general position assumption is done with no loss of generality because it can be achieved by a sequence of four symbolic perturbations of decreasing importance: (i) the input faces are translated in the x-direction by 1 , (ii) translated in the y-direction by 2 , (iii) the vector  is scaled by a factor (1 + 3 ), and (iv) the faces are rotated by an angle 4 around a line that is not parallel to the coordinate axes.As shown below, enforcing 1 2 3 4 yields that our perturbation scheme removes all degeneracies.
Consider an intersection I as defined above; I can be a line or a plane.If I is a line L that induces a degeneracy of type (δ), this degeneracy is avoided by a translation (i) in the xdirection if L is not parallel to the yz-plane, and by a translation (ii) in the y-direction, otherwise.Then, perturbations (iii) and (iv) of smaller scales do not reintroduce this degeneracy [2].If the intersection I is a plane, this remains the case after perturbations (i) and (ii), but the intersection becomes empty after a small enough perturbation (iii) and it remains empty after perturbation (iv).Hence, degeneracies of type (δ) are avoided by our scheme of perturbations.
Degeneracies of type (β) and (γ) are not affected by perturbations (i) and (ii), but they are avoided by the scaling of type (iii).Indeed, if I is a line then, viewed in projection on the back wall, the scaling of type (iii) translates the line.Finally, degeneracies of type (α) are not affected by perturbations (i-iii), but they are avoided by a rotation of type (iv).

Algorithm
We first describe the main algorithm in Section 3.1 and then two algorithmic refinements in Section 3.2.that we present separately for clarity.Our algorithm has the following property.
Theorem 1.Given a set P of polygonal faces in 3D in general position and that do not properly intersect, the algorithm outputs a simplicial complex Q whose vertices have integer coordinates and a mapping σ that maps every face F of P onto a set of faces (or edges or vertices) of Q such that there exists a continuous motion that moves every face F into σ(F ) such that (i) the L ∞ Hausdorff distance between F and its image during the motion never exceeds 3  2 and (ii) if two points on two faces become equal during the motion, they remain equal through the rest of the motion.

Main Algorithm
The algorithm is organized in 4 steps.In every step of the algorithm, faces are subdivided and/or modified.We denote by P i the set of faces at the end of Step i and by σ i the mapping from the faces of P i−1 to those of P i (with P 0 = P and P 4 = Q).These mappings are trivial and not explicitly described, except in Step 1.Let σ = σ 4 • • • • • σ 1 be the global mapping from the faces of P to those of the output simplicial complex Q.
(a) The face F j is partially projected onto F i (i < j), i.e.R ji is replaced by the faces R ji , X ji1 and, in (a), X ji2 .(c): F j is also partially projected onto 1. Collapse the faces that are close to one another.Order all the input faces arbitrarily from F1 to Fn .During the process, we modify iteratively the faces.For clarity, we denote by F i the faces that are iteratively modified, which we initially set to F i = Fi for all i.Roughly speaking, for i from 1 to n − 1, we project along y the points of F i+1 , . . ., F n onto F i but only the points that project at distance at most 1.Furthermore, we create, if needed, side walls that connect the boundary of the projected points to their pre-image.Structural properties of the faces at the end of this step are discussed in Section 3.3.Refer to Figure 1.
(a) For i from 1 to n and for j from i+1 to n, modify F j by removing the polygonal region R ji that consists of the points p j ∈ F j whose projection onto F i along the y-direction lies within distance less than 1 from p j , i.e.R ji = {p j ∈ F j | ∃ α ∈ (−1, 1), ∃ p i ∈ F i , p j = p i + α }.Let F1 , . . ., Fn be the resulting faces at the end of the two nested loops and let R ji be the projection of R ji on Fi along the y-direction. 1b) For j from 1 to n, consider on Fj the set of R ji (i < j) and consider their edges, in turn.We define new faces that connect some edges of R ji and R ji , which we refer to as connecting faces (see e.g., faces X jiζ in Figure 1).If edge e is a common edge of R ji and Fj , we define a new face as the convex hull of e and its projection on F i along y.If e is a common edge of R ji and R ji and if e projects (along y) on Fi and on Fi into two distinct segments e i and e i , respectively, we define a new face as the convex hull of e i and e i ; however, if e belongs to that face, we split it in two at e.
(c) For i from 1 to n, subdivide Fi by the arrangement of edges of the R ji , j = i+1, . . ., n.
To summarize, we have removed from every input face Fj the regions R ji , i = 1, . . ., j − 1, we subdivided the resulting faces Fi by the edges of all R ji , j = i + 1, . . ., n, and we created new connecting faces.Finally, we define σ 1 to map every input face Fj to the union of the resulting face Fj , all the regions R ji , i < j (subdivided as in Fi ), and all the connecting faces that are defined by R ji , i < j.

2.
Partition the space into slabs.Project all the faces of P 1 on the floor, compute their arrangement, lift all the resulting edges onto all faces of P 1 , and subdivide the faces accordingly.Then, project all these edges on the back wall and compute their arrangement.
The closed region bounded by the two side walls x = c ± 1 2 , c ∈ Z, is called a thin slab S c , if it contains (at least) a vertex of these arrangements.A big slab is a closed region bounded by two consecutive thin slabs.(A big slab may thus be reduced to a plane.) We subdivide all faces of P 1 by intersecting them with the side-wall boundaries of all slabs, resulting in P 2 .
3. Triangulate the faces.We triangulate all the faces of P 2 in every slab in turn.We first consider thin slabs and then big slabs.This order matters because, when triangulating faces in thin slabs, we add vertices on edges, including those on the side-wall boundaries of the thin slabs, which also belong to the adjacent big slabs.
(a) Thin slabs.Project along the x-axis all the faces in a thin slab S c on the side wall x = c and compute the arrangement of the projected edges.In that side wall, denote as hot all the pixels that contain a vertex of that arrangement and split every edge that intersects a hot pixel at its intersection with the pixel boundary.Not considering the dummy vertices, all trapezoids that project on the floor onto one and the same trapezoid are triangulated such that all the diagonals project on the floor onto one and the same diagonal.Trapezoids can only have dummy vertices on the edges on the side walls thus, after splitting a trapezoid in two triangles, each triangle can have dummy vertices on at most one of its edges.For every such triangle, we further triangulate it by adding an edge connecting every dummy vertex to the opposite vertex of the triangle.
4. Snap all vertices to the centers of their voxels.

Algorithm Refinements
We present here two algorithmic refinements that we did not describe above for clarity.

Subdivision of the faces in thin slabs (Steps 2 and 3a).
When snapping all vertices to their voxel centers in Step 4, any planar polygon in a thin slab S c is transformed into a planar polygon in the side-wall x = c.Depending on how the vertices move toward their voxel centers, the polygon may not remain planar during the motion but it is planar at the end of the motion.Hence, the output will be unchanged if, in Step 3a, we avoid triangulating the faces, that is, we avoid triangulating the arrangement in the side wall x = c and only lift the new vertices of the arrangement onto the edges in S c .Still, after snapping the vertices in Step 4, the resulting polygons in the side wall x = c should be triangulated so that the algorithm returns a simplicial complex.Doing so improves the complexity of the algorithm but it is nonetheless convenient for the proof of correctness to consider the triangulations of the faces in Step 3a.
Similarly, in thin slabs, we can avoid subdividing the faces by the vertical projections of edges in Step 2; however, since slabs are defined by the lifted edges, this means in Step 2 to vertically lift the edges on all the faces, without subdividing them, to define the slabs, and to subdivide the faces by the lifted edges only in the big slabs.

Subdivision of the connecting faces (Steps 2, 3a and 3b).
In Step 2, the connecting faces are subdivided by the side-wall boundaries of all slabs.The resulting faces are trapezoids (possibly degenerate to triangles) with two edges of length at most 1 that are parallel to the y-axis.Such trapezoids remain planar when their vertices are moved to their voxel centers in Step 4. Hence, triangulating these trapezoids does not modify the motion of any of its points.Furthermore, subdividing their edges that are parallel to the y-axis does not change the trapezoid motions either.
It follows that, in Steps 2, we do not need to subdivide the connecting faces by the vertical projections of the edges of P 1 , in Step 3a, we do not need to lift any dummy vertices on the connecting-face edges that are parallel to the y-axis and in Step 3b, we do not need to triangulate the connecting faces.Still, we should triangulate the connecting faces at the end of the algorithm in order to obtain a simplicial complex.These triangulations can trivially be done without creating proper intersections between the edges because the connecting faces that need to be triangulated are parallelograms with two edges of length 1 and parallel to the y-axis; edges are not properly intersecting at the end of Step 4, thus the edges that lie in such a parallelogram are all equal to one diagonal and we can triangulate the parallelogram accordingly.(c) The remaining part of F 4 is not projected on any of the other faces.

Properties of the faces of P 1
Observe first that, in Step 1a, R ji is polygonal since its boundary consists of (i) segments of the boundary of F j , (ii) segments of the boundary of F i projected onto F j along the y-direction, and (iii) segments of the intersection of F j and the translated copies of F i by vectors ± .This also implies that, in projection on the back wall, an edge e of R ji lies in an edge e of the boundary of an input face Fu or of Fu ∩ ( Fv ± ) for some u and v.However, the distance between e and e is not bounded from below by a constant, as depicted in Figure 3.
In Step 1b, connecting faces are defined as trapezoids that intersect any side wall in segments of length at most 1 that are parallel to the y-axis.However, it should be noted that connecting faces may overlap, as depicted in Figure 4, and that a connecting face may not contain the edge of R ji that defines it, as depicted in Figure 5.

Proof of Correctness
We prove here Theorem 1.We focus on Step 1 of the algorithm in Section 4.1, on Step 4 in Section 4.2, and we wrap up in Section 4.3.

Step 1
We first prove the main properties of Step 1 and then a technical lemma, which will be used in Lemma 6.
Lemma 2. Every point of the faces of P can be continuously moved so that every face F of P is continuously deformed into σ 1 (F ) such that (i) the L ∞ Hausdorff distance between F and its image during the motion exceeds 1 and (ii) if two points on two faces become equal during the motion, they remain equal through the rest of the motion.
Proof.The motion is decomposed in n successive phases, considering the projection on each F i in turn.For a particular F i , the naive way of moving This can be done so that when two distinct points become equal during the motion, they remain equal through the rest of the motion.The formal proof goes as follows.
Recall that F1 , . . ., Fn denote the input faces of P and that we initially set F i = Fi .In Step 1, for i from 1 to n and for j from i + 1 to n, we modify F j by projecting its subset R ji onto R ji ⊂ F i , and by adding the corresponding connecting faces.When we start the projection on F i , faces F u , for u i, are no longer modified and we have F u = Fu until the end of Step 1.

Subdivision of R ji .
For all i and all j ∈ {i + 1, . . ., n}, we project parallel to y the R ji on the back wall, triangulate the resulting arrangement, lift the triangulation back onto the R ji , and subdivide the R ji accordingly.Let T 1 , . . ., T g denote the resulting triangles.
By construction, the boundary of T u ⊂ R ji consists of parts of (i) boundary edges of R ji that define connecting faces (i.e., any boundary edge incident to R ji and some R ji , and that projects along y onto Fi and Fi into two distinct segments), (ii) boundary edges of R ji that do not define connecting faces, and (iii) edges that are induced by the subdivision of the R ki into triangles.Edges of types (ii) and (iii) will be moved onto Fi and they are referred to as red, and edges of type (i) will remain invariant, at first, and they are referred to as blue.A vertex is considered blue if it is the endpoint of a blue edge (among all its incident edges in T 1 , . . ., T g ); otherwise, it is red.

Ordering property.
Consider two triangles T b ⊂ R bi ⊂ Fb and T r ⊂ R ri ⊂ Fr that are on the same side of Fi (with respect to y) and that project onto the same triangle on the back wall.We prove that, if a blue edge e b of T b and a red edge e r of T r coincide in that projection, then T b is farther away than T r from Fi with respect to y (i.e.any line parallel to the y-axis, that intersects these triangles, intersects Fi , T r and T b in that order).
In projection on the back wall, an edge of R ji lies in ∂ Fu or Fu ∩ ( Fv ± ) with u i or u = j, and v i.Indeed, R ji only depends on F i and F j , at the time when R ji is defined in Step 1, and then, F i only depends on , . . ., Fi and F j only depends on F1 , . . ., Fi−1 and Fj .It follows that Fr plays no role in the definition of R bi and thus that the boundary edges of Fr do no overlap the edges of R bi in projection on the back wall, by items (α) and (γ) of the general position assumption.In particular, e r cannot be a boundary edge of Fr .
Thus, since e r is a red edge, either (a) e r lies on a common edge of R ri and some R ri (on Fr ) such that e r projects (along y) on Fi and on Fi onto the same segment e ri = e ri , or (b) e r is induced by the subdivision of the R ki into triangles.Subdivision of T 1 , . . ., T g .
In the rest of the proof, unless otherwise specified, we consider i fixed and we consider among all triangles T 1 , . . ., T g , only those that project in Step 1 on Fi (i.e. that belong to R ji for some j) onto the same triangle, and, furthermore, that are on one same side of Fi (with respect to y).
We choose a parameter small enough and subdivide all these triangles T j as follows.Refer to Figure 7.The general idea is to define a small part of the triangle in the neighborhood of blue edges and vertices.More precisely, for a blue vertex v incident to one or two red edges of the triangle, we define the retractions of v to be the points on these red edges that are, viewed on the back wall at distance λ from v, where λ is the distance between v and its projection on Figure 7: Back wall views of the subdivision of three triangles T j , T k , T that project on Fi into the same triangle and such that Fi , T j , T k , T appear in that order along y.
Fi along y; for a red vertex, its retraction is the vertex itself; and for a blue vertex incident to two blue edges, we define its retraction to be located on the triangle at some relevant place at distance, viewed on the back wall, λ from v.Then, we define the retraction of a triangle T j to be Tj the convex hull of its retracted vertices.
Considering the triangles from the closest to the farthest from their projection of Fi , the ordering property ensures that the blue vertices that are incident to two blue edges can be retracted as defined above such that, on the back wall, the Tj are nested so that the inner of two Tj and Tk is the one that is the farthest away from Fi (along y).Notice that, among all triangles T 1 , . . ., T g , if T j and T k share a red edge (in R 3 ), then Tj and Tk share the same vertices on this red edge.We define the motions of the triangles in four phases.In the first phase (see Figure 8(a)), we infinitesimally deform T j into the surface defined as the union of T j \ Tj , a translated copy of Tj toward Fi along y, and, for each blue edge or isolated blue vertex of T j , a small parallelogram that connects these two parts.This can trivially be done continuously without creating intersection for any two triangles T j and T k that are disjoint or if their projections on the back wall are distinct.On the other hand, if T j and T k share an edge and project (on the back wall) on the same triangle, the shared edge is necessarily red (since input faces have disjoint relative interior).The fact that their subdivision polylines are distinct on the back wall except for a common endpoint on the red edge ensures that the motion can be performed so that the relative interiors of the deforming faces remain disjoint.
In the second phase (see Figure 8(b)), all infinitesimally translated faces Tj ±  are moved to Inria Fi ; this is done trivially my moving each point onto Fi at constant speed parallel to y.The faces T j \ Tj remain invariant.This defines the motions of the vertices of the infinitesimal wall added in the first phase; the points in their relative interiors move accordingly to their barycentric coordinates in their faces.Since the subdivision polylines are disjoint in the back wall except at points where the polylines coincide in R 3 , any two distinct points remain distinct during the motion except possibly at the end of the phase when they may become equal on Fi .
In the third phase (see Figure 8(c)), (i) every face of T j \ Tj is shrunk to the blue edge (or isolated blue vertex) that defines it, (ii) the image of Tj at the end of the second phase, is the projection of Tj onto Fi (parallel to y), is expanded to the projection of T j on Fi , and (iii) the image of the infinitesimal walls of the first phase are moved to the connecting faces (or edges) between T j and its projection onto Fi .
Refer to Figure 8(d) and consider all the triangles T j that project (parallel to y) onto one and the same triangle T on Fi .In the definition of the subdivision of the T j , the Tj are nested in the back wall.Denote by T 1 ⊂ • • • ⊂ T h ⊂ T their projections (parallel to y) on Fi ordered from the innermost to the outermost.For simplicity, define T h+1 = T .All these polygons are close to T (for the Hausdorff distance) and there is thus a natural correspondance between their vertices.
On Fi , we move each vertex of T j at constant speed toward its corresponding vertex on T j+1 .For simplicity, we consider these motions one after the other, for j = 1 to h.During the j-th motion, the points inside T j move accordingly to their barycentric coordinates in T j , and the points in T \ T j+1 are invariant.The points on T k \ Tk move on T k accordingly to the motion of their projection (parallel to y) on Fi .This defines the motion of the points on the boundary of all the faces that are parallel to y and the points in their relative interiors move accordingly to their barycentric coordinates.It is straightforward that any two distinct points remain distinct during this motion except possibly at the end of the phase when some points (viewed on the back wall on the boundary of T ) may become equal.
Finally, the fourth phase goes as follows.Consider, as in Figure 5, a common edge e of R ji and R ji (on Fj ) such that e projects (along y) on Fi and on Fi into two distinct segments e i and e i , respectively, that are on the the same side of Fj with respect to y (i.e., e does not belong to the convex hull of e i and e i ).In Step 1b, we define a connecting face as the convex hull of e i and e i .However, since e is a blue edge and has thus not yet been moved at this stage, R ji and R ji have currently been deformed into a set of faces that contains the convex hulls of e and e i , and the convex hull of e and e i .These two faces overlap between e and, say, e i .In this phase, we simply retract these overlapping parts into segment e i .
Hausdorff distance.The property that the Hausdorff distance between a face and its image during the motion never exceeds 1 (for the L ∞ metric) is straightforward since, (i) any line parallel to y that intersects a triangle T j also intersects its image during the motion, (ii) the image of T j during the motion remains in the convex hull of T j and its projection (parallel to y) on Fi and (iii) all the points of T j are at distance at most 1 along y from Fi (by definition of R ji ).Lemma 3. If a line L parallel to the y-axis intersects the relative interior of a face of P 1 in a single point p then the distance along L from p to any other face of P 1 is at least 1.
Proof.Observe first that any point in a connecting face lies on a segment of length at most 2, parallel to the y-axis and with its endpoints on two non-connecting faces.Indeed, connecting faces are defined either as (i) the convex hull of an edge on Fj and its projection on Fi along y at distance at most 1, or (ii) the convex hulls of an edge on Fj and, respectively, its projections on Fi and Fi in opposite directions along y and each at distance at most 1.
If a line L parallel to the y-axis intersects the relative interior of a face F of P 1 in a single point p, this face is not a connecting face.Assume for a contradiction that L intersects another face F of P 1 at distance less than 1 from p.If F is a connecting face then, by the above observation, L also intersects a non-connecting face at distance at most 1 from p. Furthermore, this distance is exactly 1 only if the above segment has length exactly 2 with p its midpoint; but this is impossible since p is in the interior of F , and faces do not transversally intersect, by Lemma 2. Thus, L intersects a non-connecting face at distance less than 1 from p.In other words, we can assume that F is a non-connecting face.
Since non-connecting faces are not parallel to the y-axis, there exists a line L parallel to the y-axis (and close to L) that intersects the relative interiors of both F and F in two points at distance less than 1.This is impossible after Step 1, which concludes the proof.

Step 4
In the following, we consider in the snapping phase of Step 4 a continuous motion of the vertices such that every vertex moves on a straight line toward the center of its voxel at a speed that is constant for each vertex and so that all vertices start and end their motions simultaneously.The motion of the other points in a face move accordingly to their barycentric coordinates in the face.Note that, in every voxel that contains a vertex, the motion is an homothetic transformation whose factor goes from one to zero.During that motion, we consider that thin and big slabs respectively shrink and expand accordingly.
We recall the standard snap-rounding result for segments in two dimensions.A pixel is called hot if it contains a vertex of the arrangement of segments.
Theorem 4 ([7, Thm.1]).Consider a set of segments in 2D split in fragments at the hot pixels boundaries and a deformation that (i) contracts homothetically all hot pixels at the same speed4 and (ii) moves the fragments outside the hot pixels according to the motions of their endpoints.During the deformation, no fragment endpoint ever crosses over another fragment.Lemma 5. When moving all vertices to the center of their voxels in Step 4, no two faces, edges, or vertices of P 3 properly intersect in thin slabs.
Proof.Consider all the faces of P 3 in a thin slab S c and the arrangement of their projections (along the x-axis) onto the side wall x = c.In that side wall, a pixel that contains a vertex of the arrangement is hot and every edge (in that side wall) that intersects a hot pixel is split at the pixel boundary (Step 3a).By Theorem 4, when moving in that side-wall all the vertices to the centers of their pixels, the topology of the arrangement does not change except possibly at the end of the motion, where edges and vertices may become identical.
It follows that the property that every face of P 3 in S c projects onto a single face of the arrangement in the side wall x = c, which holds by construction at the beginning of the motion (Step 3a), holds during the whole motion of the vertices in 3D and of their projections in the side wall x = c.
Furthermore, the motion preserves the ordering of the x-coordinates of the vertices in S c , until the end when they all become equal to c. Together with the previous property, this implies that, in a thin slab, during the snapping motion, (i) no vertices and edges intersect the relative interior of a face and (ii) if two edges intersect in their relative interior, it is at the end of the motion and they become identical.Furthermore, (iii) no vertices intersect the relative interior of an edge because, in Step 3a, we have split every edge that intersects a hot pixel in projection in the side wall x = c.This concludes the proof.
Triangulated trapezoid at time t = t 1 when edge e (in blue) intersects edge e (and two others; the top edge of T is intersected by e only in projection on the back wall).
Respective positions of M 1 , M 2 and M , (b) inside a column of pixel before snapping the vertices and (c) after snapping the vertices.
In the left side wall: Proof.By construction (Step 2), all the vertices in a big slab are on its side-wall boundaries and, in these side walls, no two edges or vertices properly intersect during the motion, by Lemma 5. Thus, we only have to consider edges that connect the two side-wall boundaries of a big slab and show that such edges do not properly intersect during the snapping Note that input faces are not vertical by assumption and connecting faces are not vertical in big slabs since they are parallel to the y-axis.
Initially, these edges project on the floor onto edges that do not properly intersect pairwise (by definition of the slabs in Step 2).Thus, by Theorem 4, the projections on the floor of two edges either (i) coincide throughout the whole motion, or (ii) they do not properly intersect and do not coincide throughout the whole motion except possibly at the end when they may coincide.In the first case, throughout the whole motion, the edges belong to the same moving vertical plane and they do not properly intersect since they do not initially; indeed, since faces are not vertical, edges may intersect in a vertical plane only if they are boundary edges of trapezoids of P 2 , and such edges do not properly intersect on the back wall by definition of big slabs.Hence, only in the latter case (ii), two edges may properly intersect during the motion, and the first time this may happen is at the end of the motion and then, the two edges belong to the same vertical plane.
Similarly, applying Theorem 4 to the back-wall projection of the boundary edges of the trapezoids, we get that if two boundary edges of trapezoids properly intersect in 3D during the motion, it is at the end and they must coincide in the back-wall projection.Since two edges that coincide in two projections are equal, we get that boundary edges of trapezoids cannot properly intersect throughout the motion.It remains to prove that there is no proper intersections that involve the edges triangulating the trapezoids.
Consider for a contradiction two edges e and e that properly intersect in a vertical plane V at time t 1 , the end of the motion.Since boundary edges of the trapezoids do not properly intersect, we can assume without loss of generality that one of the two edges, say e, is a triangulation edge.Consider at time t 1 , the triangulation of the (deformed) trapezoid F that contains e.We prove that (at time t 1 ) edge e properly intersects (at least) one of the two boundary edges of F .
Assume for a contradiction that e properly intersects none of the two boundary edges of F and refer to Figure 9(a).Consider all the edges of the triangulation of F that are properly intersected by e and the sequence of triangles (of that triangulation) that are incident to these edges; let T and T denote the first and last triangles of that sequence.All these triangles except possibly one, T or T , must be in the vertical plane V ; this is trivial for all triangles but T and T and, if neither T nor T lies in V , then edge e properly intersects the surface formed by these triangles, contradicting the property that t 1 is the first time a proper intersection may occur.As in Figure 9(a), assume without loss of generality that T , the bottommost triangle of the sequence, lies in the vertical plane V at time t 1 .Let M be the endpoint of e that lies in T and let M 1 M 2 be the edge of T that supports M .At time t 1 , M 1 and M 2 are vertically aligned and M is in between them.Thus, before the snapping motion starts, at time t = t 0 , M 1 , M 2 and M must lie in the same vertical column of pixel (in the side wall) and M must be vertically in between M 1 and M 2 (see Figure 9(b)).However, the distance along the y-axis between M and segment M 1 M 2 is at least 1 by Lemma 3. Thus, M and M ±  ∈ M 1 M 2 are initially on opposite sides of the column of pixels and thus that have half-integer x and y-coordinates, which contradicts item (δ) of our general position hypothesis.
Hence, at time t = t 1 , edge e properly intersects one of the boundary edges, say r, of trapezoid F .Since boundary edges do not properly intersect, e must be a triangulation edge of its trapezoid F and we can apply the same argument as above on edges e and r, instead of e and e .We get that r properly intersects a boundary edge r of F , which is a contradiction.

Wrap up, proof of Theorem 1
First, by construction, the algorithm outputs faces that have integer coordinates.
Second, there is a continuous motion of every input face F into σ(F ) so that the Hausdorff distance between F and its image during the motion never exceeds 3  2 for the L ∞ metric.Indeed, by Lemma 2, the Hausdorff distance never exceeds 1 between F and its image during the motion Step 1; in Steps 2 and 3 the faces are only subdivided; and the Hausdorff distance between any face of P 3 and and its image during the motion of Step 4 clearly never exceeds 1  2 since vertices are moved to the centers of their respective voxels.
Third, if two points on two faces become equal during the motion, they remain equal through the rest of the motion.This is proved in Lemma 2 for the motion of Step 1 and this also holds for the motion of Step 4 since, by Lemmas 5 and 6, if two faces, edges or vertices intersect during this motion, they share a common face of both, whose motion is uniquely defined by its vertices (actually, we show in the proofs of Lemmas 5 and 6 that no two distinct points become equal except possibly at the end of the motion).
Finally, the algorithm outputs a simplicial complex by Lemmas 5 and 6.

Complexity
We first analyse in Section 5.1 the complexity of our algorithm in the worst case (Proposition 7) and then, in Section 6, its complexity under some reasonable assumptions on the input (Proposition 10).We finally argue in Remark 12 that we can anticipate time and space complexities of O(n √ n) in practice on realistic non-pathological data.

Worst-Case Complexity
We define the z-cylinder of a face F as the volume defined by all the lines parallel to the z-axis that intersect F ; similarly for x and y-cylinders.Over all input faces F , let f d be the maximum number of input faces that are (i) intersected by one such cylinder of F and (ii) at distance at Inria most d from F .Denote f = f ∞ the maximum number of faces intersected by one such cylinder.Let w x be the the maximum number of input faces that are intersected by any side wall x = c.
Proposition 7. Given a set of O(n) polygons, each of constant complexity, the algorithm outputs a simplicial complex of complexity O(nw 2 x f 7 f5 1 ) in time O(n 2 w x f 9 f 7 1 ).
Proof.For analyzing the complexities of the algorithm and of its output, we bound the number of edges that we consider for splitting the input faces.However, for each face, we first count the number of edges without considering any intersection; in other words, for each face, we bound the number of lines supporting the edges instead of the complexity of the induced arrangement.
To avoid confusion, we refer to these edges as unsplit edges.).)The number of thin and big slabs is thus O(nf 5 f 2 1 ).Complexity in thin slabs.For any thin slab S c , we analyze the complexity of the output in the side wall x = c.We do not consider here any triangulation edge inside the faces as discussed in Section 3.2.1.Thus, the edges in S c are subdivisions of the edges of P 1 and subdivisions of the edges defined as the intersection of the faces of P 1 and the side walls x = c ± 1 2 .More precisely, these edges are pieces of either (a.i) edges of P 1 that lie on the input faces, (a.ii) edges of intersection of an input face with a side wall x = c ± 1 2 , (b.i) edges of connecting faces that are parallel to the y-axis, or (b.ii) edges of intersection of a connecting face with a side wall x = c ± 1 2 .Note that, although we do not consider the faces triangulated in S c , the edges are subdivided according the arrangement of their projections on the side wall x = c; however, we consider them unsplit for now.
As mentioned above, every input face supports O(f f 1 ) unsplit edges at the end of Step 1.Thus, if F Sc denotes the number of input faces that are intersected by the thin slab S c , there are O(F Sc f f 1 ) edges of types (a) in S c .On the other hand, the O(f f 1 ) unsplit edges on the back wall yield an arrangement of size O(f 2 f 2 1 ); viewed on the back wall, edges of type (b) are incident to the vertices of this arrangement, so the number of edges of type (b) in S c is bounded by O(F Sc f 2 f 2 1 ).We bound the number of intersection between these edges at the end of the algorithm.However, as noted in Section 3.2.2,we do not consider intersections that involve edges of type (b) because such intersections necessarily belong to the hot pixels defined by the edge endpoints.We thus consider here only edges of type (a).Every edge on a given input face F may only intersect, after projection on the side wall x = c, edges that lie on faces that intersect S c and the x-cylinder of F ; 5 there are O(f 1 ) such faces and O(f f 1 ) edges on each of them, thus every edge may intersect O(f f 2 1 ) edges on the side wall x = c.There are O(F Sc f f 1 ) edges of type (a) in S c , thus, there are O(F Sc f 2 f 3 1 ) intersections between edges of type (a).
In addition to these O(F Sc f 2 f 3 1 ) intersections, there are O(F Sc f 2 f 2 1 ) edges (and thus edge endpoints) in S c .The number of hot pixels in the side wall x = c is thus H c = O(F Sc f 2 f 3 1 ) at the end of Step 3a.
The number of vertices (dummy or not) in S c at the end of Step 3a is the number of hot pixels, H c , times the number of edges in S c , O(F Sc f 2 f 2 1 ).However, the size of the arrangement in the side wall x = c after snapping the vertices to their voxel centers in Step 4 is of the order of the number hot pixels, H c , and it remains as such after triangulating the faces (see Section 3.2.1).
Complexity in big slabs.The number of edges in a big slab after the triangulation of Step 3b is (i) the number of edges of connecting faces in the slab, that is O(w x f f 1 ) (the number O(f f 1 ) of unsplit edges on input faces times the number O(w x ) of input faces that intersect the slab) plus (ii) the number O(w x ) of input faces that intersect the slab times the number H c1 + H c2 of hot pixels in the adjacent thin slabs S c1 and S c2 .This sums up to O(w x (H c1 + H c2 )).
Complexity of the output.The total complexity of the output is thus O(w x ) times the sum over all thin slabs of the number of hot pixels H c in S c .Denoting by ). Times w x gives the complexity of the output: O(nw 2 x f 7 f 5 1 ).Time complexity.All the arrangements and triangulations performed by the algorithm can be done in time complexities that match their worst sizes.However, this does not match the complexity of the output because the complexity of P 3 may be larger before than after snapping.As noted above, the complexity of P 3 in a thin slab is O(H c F Sc f 2 f 2 1 ).On the other hand, the complexity of P 3 in a big slab is, as above, O(w x (H c1 + H c2 )).The sum of these complexities over all slabs is bounded by nf 2 f 2 1 times the total number O(nw x f 7 f 5 1 ) of hot pixels, that is O(n 2 w x f 9 f 7 1 ).

Complexity under some Assumptions
We consider, in Proposition 10, the complexity of our algorithm for approximations of "nice" surfaces, defined as follows.
Definition 8.A (ε, ζ)-sampling of a surface S is a set of vertices on S so that there is at least 1 and at most ζ vertices strictly inside any ball of radius ε centered on S. It is straightforward that a (ε, ζ)-sampling of a compact surface has Θ(n) vertices with n = 1 ε 2 .The Delaunay triangulation of a set of points P restricted to a surface S is the set of simplices of the Delaunay triangulation of P whose dual Voronoi faces intersect S. If P ⊂ S, we simply refer to the restricted Delaunay triangulation of P on S.
A surface S is k-monotone (with respect to z) if every line parallel to the z-axis intersects S in at most k points.Let ∆ and k be any two positive constants.A surface S is nice if it is a compact smooth k-monotone surface such that the Gaussian curvature of S is larger than a positive constant in a ball of radius ∆ centered at any point p ∈ S where the tangent plane to S is vertical.
The following lemma is a technical though rather straightforward result.Remark 11.It should be stressed that the above complexities do not depend on the grid size, including in the hidden constants.This is important because for a fixed grid size, the number of voxels that are intersected by a given compact surface (or its convex hull) is a constant C hence, independently of n, the complexity of the snap-rounded simplicial complex will be in O(C) = O(1).It follows that, for a fixed grid size, it would be pointless to consider the asymptotic complexity of the simplicial complex in terms of n.However, this makes sense if this complexity is independent from the grid size, which is the case in the above proposition.
Remark 12.In practice on realistic non-pathological data, one can anticipate a better complexity of O(n √ n) for the size of the output and for the time complexity.
Indeed, the x, y or z-cylinders of most faces will intersect a constant number of other faces and we can expect that the few that intersect a non-constant number of faces will not impact the final complexities; hence we can consider that f 1 < f = O(1) in Proposition 7.Moreover, considering the proof of Proposition 7 with f 1 < f = O(1), we get that the number of hot pixels in a thin slab S c is H c = O(F Sc ), that the number of edges in S c is O(F Sc ), and that the number of edges in a big slab is O(w x ).In practice, this should yield, before snapping, arrangements of complexity O(w x + F Sc ) in a big and thin slab instead of the worst-case complexities O(w x F Sc ) and O(F 2 Sc ) for big and thin slabs, respectively.Summing over all slabs, this yields an anticipated practical complexity of O(nw x ) where w x = O( √ n) as in Lemma 9.

Conclusion
Natural questions arise.The algorithm we presented is simple enough that it should be implementable without particular difficulties.However, its worst-case complexity, even under reasonable assumptions (Propositions 7 and 10), is prohibitive in practice.Hence, the question of whether our estimated practical complexity of O(n √ n) is correct on real data is crucial for applications.Furthermore, it would be interesting to design heuristics for improving the practical efficiency of the algorithm.Another issue is that faces can drift arbitrarily far when the snap rounding scheme is applied repeatedly.Several approaches were presented to address this issue in 2D [8,10,14] but the problem in 3D is naturally entirely open.

Figure 1 :
Figure 1: View inside a side wall x = cst.(a-d):The face F j is partially projected onto F i (i < j), i.e.R ji is replaced by the faces R ji , X ji1 and, in (a), X ji2 .(c): F j is also partially projected ontoF k (i < k < j).(d): If instead i < j < k, it is F k that is partially projected onto F j .

2
Triangulate the resulting arrangement3 and lift it back (still along the x-axis) onto all the faces in the slab and subdivide them accordingly.All these new vertices are referred to as dummy vertices to differentiate them from the other vertices.(b) Big slabs.Refer to Figure2.Not considering the dummy vertices of Step 3a, all faces are trapezoids (possibly degenerated to triangles) such that each of the parallel edges lie on each of the side-wall boundaries of the big slab; any two trapezoids are either identical, disjoint, or share exactly one edge or vertex, and the same holds for their projections on the floor.The dummy vertices lie on the trapezoid edges that lie on the side-wall boundaries of the big slab.
Triangulation of trapezoids not considering dummy vertices.

Figure 2 :
Figure 2: Triangulations of trapezoids in a big slab.

R 87 Figure 3 :
Figure 3: View in a side wall: the boundary edge e of R can be far way from the edge e = F1 ∩ ( F2 − ) that coincides with e on the back wall.

Figure 4 :
Figure 4: Example of overlapping connecting faces X 31 and X , viewed inside a side wall.(a): The four input faces; (b) F 3 is projected onto F 1 , then (c) F 4 onto F 1 , and then (d) F 3 onto F 2 .The remaining part of F 4 is not projected on any of the other faces.

Figure 5 :
Figure 5: A connecting face X induced by the common edge e of R 31 and R 32 , and that does not contain e.View inside a side wall: (a) The three input faces; (b) F 2 is projected onto F 1 , then (c) F 3 is projected onto F 1 and F 2 .
by moving each point of R ji along the y-direction and at constant speed to R ji , and to transform edge e ζ into face X ζ for every edge e ζ of the boundary of R ji that defines a connecting face X ζ .However, this does not define a function since segments are mapped to faces.The definition of a continuous motion requires a bit of technicality but the straightforward underlying idea is to subdivide R ji by considering, for each edge e ζ , a tiny quadrilateral bounded by e ζ .We then transform continuously each tiny quadrilateral bounded by e ζ into the connecting face X ζ and move the complement of these quadrilaterals, which is a slightly shrunk version of R ji , into R ji .

FiFigure 6 :
Figure 6: Counter example for proof of the ordering property.Refer to Figure 6.In case (a), since input are interior disjoint, the segment e ri = e ri is on the common boundary of the input faces Fi and Fi .Hence, e r is a common edge of T r ⊂ R ri and some T r ⊂ R ri such that they respectively project in Step 1 onto two triangles T i ⊂ Fi and T i ⊂ Fi that share edge e ri = e ri .Assume for contradiction that edge e b is in between e r and e ri (i.e. in their convex hull).Since e b is a blue edge, the triangle T b ⊂ Fb that share edge e b with T b projects in Step 1 onto a triangle T b = T i .Consider any line parallel to the y-axis that intersects the interior of the four triangles T b , T i , T b and T r in points p b , p i , p b and p r , respectively.By definition segment p r p i has length less than 1 and contains p b .Thus p b p i < 1 and since T b projects in Step 1 onto T b = T i , we have b < i .Furthermore, i < r by definition of R ri .Thus b < i < r and since T r projects on T i instead of T b , p r p b > 1 while p b p b < 1.Hence, p b , p i , p b and p r appear in that order on their supporting line.It follows that p i p b < 1 and thus T i should have been projected on T b in Step 1, which contradicts the definition of R ri , and thus e b is not in between e r and e ri .Since e r and e b are on the same side of Fi , e b is farther away than e r from Fi .Case (b) is similar: the only difference is that T i lies in Fi instead of Fi but T i and T i are still incident.The ordering property follows.

Figure 8 :
Figure 8: Motion of T j . of the T .We define the motions of the triangles in four phases.In the first phase (see Figure8(a)), we infinitesimally deform T j into the surface defined as the union of T j \ Tj , a translated copy of Tj toward Fi along y, and, for each blue edge or isolated blue vertex of T j , a small parallelogram that connects these two parts.This can trivially be done continuously without creating intersection for any two triangles T j and T k that are disjoint or if their projections on the back wall are distinct.On the other hand, if T j and T k share an edge and project (on the back wall) on the same triangle, the shared edge is necessarily red (since input faces have disjoint relative interior).The fact that their subdivision polylines are distinct on the back wall except for a common endpoint on the red edge ensures that the motion can be performed so that the relative interiors of the deforming faces remain disjoint.In the second phase (see Figure8(b)), all infinitesimally translated faces Tj ±  are moved to

Figure 9 :
Figure 9: For the proof of Lemma 6.

Lemma 9 . 1 4Proposition 10 .
The restricted Delaunay triangulation T of a (ε, ζ)-sampling on a compact surface has complexity O(n) = O( 1 ε 2 ).Any plane x = c intersects at most O( √ n) = O( 1 ε ) faces of T .Furthermore, for any face f of T , the set of vertical lines through f intersects at most O(n ) = O( 1 √ ε ) faces of T .Inria hence the area of S in the cylinder is O(ε √ ε).The number of sampling points in the cylinder is thus O( Given the arrangement of the restricted Delaunay triangulations of the (ε, ζ)samplings of a constant number of nice surfaces, the algorithm outputs a simplicial complex of complexity O(n 5 ) in time O(n 6 √ n).Proof.Consider the restricted Delaunay triangulations of the (ε, ζ)-samplings of two surfaces.By definition of (ε, ζ)-samplings, it is straightforward that any triangle of one triangulation intersects a constant number of triangles of the other triangulation.Hence, the complexities of Lemma 9 hold for the arrangement of the two triangulations, and similarly for a constant number of triangulations.Thus, for the arrangement of triangulations, Lemma 9 yields w x = O( √ n) and f 1 < f = O( 4 √ n) (as defined in Section 5.1) and plugging these values in the complexities of Proposition 7 yields the result.
Number of slabs.After projection on the back wall, the edges of R ji and R ji , in Step 1, are pieces of the boundary edges of the input faces Fk and of the segments of intersection Fk ∩( F ± ).In a y-cylinder of a face F , only f faces project on the back wall and thus there are O(f f 1 ) such edges.Thus, at the end of Step 1, every input face ends up supporting O(f f 1 ) unsplit edges.InStep 2, we thus lift O(f 2 f 1 ) unsplit edges onto every face.Every unsplit edge on a given face F may only intersect, after projection on the back wall, edges that lie on the faces that intersect the y-cylinder of F ; there are O(f ) such faces and O(f 2 f 1 ) unsplit edges on each of them, thus every unsplit edge may intersect O(f 3 f 1 ) edges on the back wall.There are O(nf 2 f 1 ) unsplit edges in total, hence, in Step 2, the back wall arrangement has complexity O(nf 5 f 2 1 ).(Similarly, the floor arrangement has complexity O(nf 3 f 2 1 Sc the number of input faces that are entirely inside S c , we have that F Sc Sc + 2w x and that the sum over all thin slabs of Sc is O(n).Hence, the sum over all N = O(nf 5 f 2 1 ) thin slabs of the numbers of hot pixels H • F • F • F