Fast and robust computation of the Hausdorff distance between triangle mesh and quad mesh for near-zero cases

Computers & Graphics 81 (2019) 61–72

Contents lists available at ScienceDirect

Computers & Graphics journal homepage: www.elsevier.com/locate/cag

Special Section on CAD & Graphics 2019

Fast and robust computation of the Hausdorff distance between triangle mesh and quad mesh for near-zero cases Yunku Kang a, Seung-Hyun Yoon b, Min-Ho Kyung c, Myung-Soo Kim a,∗ a

Department of Computer Science and Engineering, Seoul National University, Seoul 08826, South Korea Department of Multimedia Engineering, Dongguk Univeristy, Seoul 04620, South Korea c Department of Digital Media, Ajou University, Suwon 16499, South Korea b

a r t i c l e

i n f o

Article history: Received 8 March 2019 Revised 28 March 2019 Accepted 29 March 2019 Available online 16 April 2019 Keywords: Hausdorff distance Shape matching Quad mesh

a b s t r a c t We introduce an algorithm for computing the two-sided Hausdorff distance between a triangle mesh and a quad mesh, guaranteed to be within the given error bound, which can be machine precision-level small. The algorithm expands upon a recent breakthrough that only calculates the one-sided Hausdorff distance from the triangle mesh to the quad mesh using what is called “matching” and “upper bounding” of candidate pieces. We complete the algorithm by accomplishing the computation of the one-sided Hausdorff distance in the opposite direction: from the quad mesh to the triangle mesh. We split each quad into two triangular pieces to simplify the breakdown of matching cases and provide additional matching methods for new cases. By fusing the two one-sided computation algorithms, one can compute the two-sided Hausdorff distance that, for instance, can properly evaluate a quad mesh approximation of a triangle mesh. Experimental results show that our algorithm can handle near-zero Hausdorff distance, which has always been known to be a much difficult task, in an interactive time. Moreover, the improvement in efficiency of the two-sided Hausdorff distance computation over the successive execution of the two one-sided computations is addressed. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Quadrilateral meshes with structural regularity are frequently preferred to triangle meshes in certain applications. They are often better suited to represent principal curvature and feature curves, and parameterization along the surface is much more straightforward. Another significance is that the hexahedral final element mesh, with a quad mesh as its boundary surface, produces less error than the tetrahedral mesh [1]. Thus, conversion of a triangle mesh into a quad mesh has long been a topic of interest [2]. To evaluate a quad mesh approximation, the Hausdorff distance (HD) can be used. The HD between two geometric shapes is a measure that can demonstrate the similarity between the two. However, its computation is notable for its complexity, especially growing substantially more difficult as the HD nears zero—when the two shapes are identical. HD computation for such cases have seldom been explored [3,4], much less to the degree of machine precision.

∗

Corresponding author. E-mail address: [email protected] (M.-S. Kim).

https://doi.org/10.1016/j.cag.2019.03.014 0097-8493/© 2019 Elsevier Ltd. All rights reserved.

Recently, an algorithm has been proposed that quickly calculates the one-sided near-zero HD from a triangle mesh to a quad mesh [5]. It works by iteratively eliminating parts of the triangle mesh (source) that is known to not contain the point where the HD occurs. For every triangle of the source mesh being considered, every point on the triangle is assigned to a point on the destination quad mesh, yielding a “local upper bound” of HD for the piece. If the local upper bound is less than the known lower bound for the HD, the triangle is no longer considered. In this paper, we present an adaptation of the algorithm that computes the HD in the other direction—from quad mesh to triangle mesh. Combining the two results makes possible the computation of the whole two-sided HD between the two meshes, that can be used to evaluate a quad mesh approximation of a triangle mesh. We cover two kinds of interpretation of a quad: two adjoining triangles or a single bilinear surface. The second case is especially challenging because it introduces smoothness into the mesh, affecting our ability to make assumptions about the shape. Nevertheless, we achieve a reasonable speed by devising matching methods for quads. The rest of the paper is organized as follows. In Section 2, we briefly overview recent HD computation algorithms that bear the most significance. Section 3 will illustrate our overall computa-

62

Y. Kang, S.-H. Yoon and M.-H. Kyung et al. / Computers & Graphics 81 (2019) 61–72

tion algorithm for the one-sided HD from quad mesh to triangle mesh, while Section 4 describes the details in accordance to the differences stemming from the interpretation of the quads. The two-sided HD computation and experimental results on the performance of our algorithm are discussed in Section 5.

However, the lower bound only covers one half of the story. We also wish to guarantee a (global) upper bound for the HD. Consider the quad mesh A being split into multiple pieces, A1 , . . . , An , some of which may contain the p we are looking for. Then the following holds true by definition.

2. Related work

h(A, B ) = max h(Ak , B ) k

The amount of previous research results on Hausdorff distance computation for polygonal mesh models shows its difficulty. Metro [6] and MESH [7], the first two published results, merely return an approximation of the HD, by only taking point samples and calculating their HD. Tang et al. [8] report real time results, but their space complexity tends to blow up for near-zero HD. Meanwhile, Bartonˇ et al. [9] have presented a bulky O(n4 log n) algorithm. It shows how the complexity of the problem is much higher than the HD for NURBS surfaces; the combinatorial possibilities for self-bisectors of the mesh is too numerous. Kang et al. [5] are the only ones to have tackled near-zero HD computation for meshes, but it is limited to the one-sided HD from triangle mesh to quad mesh, although it can also handle the two-sided HD between two triangle meshes. We extend this result to the HD of the opposite direction, fulfilling the two-sided HD computation. The only two other instances where near-zero HD is mentioned compute the HD between NURBS surfaces. Hanniel et al. [3], while presenting a GPU-based sampling approach, show one near-zero example where the algorithm successfully returns, but the difference from the true HD is greater than 10−2 of the model size. On the other hand, Kim et al. [4] guarantee the HD within 10−3 of the model size. They introduce the “matching” of NURBS surfaces as a way to discard irrelevant candidates, also employed by Kang et al. [5], whose error bound is limited only by machine precision. This accuracy is enabled in part by the flatness of meshes, allowing greater certainty when removing unnecessary candidates. But it is nevertheless a notable feat when other factors are considered, such as the number of elements in the models or the treatment of quads as smooth (bilinear) surfaces. 3. Algorithm overview For geometric surfaces A and B, which in our case are a quad mesh and a triangle mesh, the one-sided Hausdorff distance from A to B is defined as follows. (For the most part, we will refer to one-sided HD when we simply say HD.)

h(A, B ) = max min p − q p∈A

q∈B

(1)

For each point in A, we consider the minimum distance to B, and then the HD is the maximum of those distances. Our goal is to find such p, which can obviously be done by examining more and more points in A. In other words, we use the following property of the HD, where A denotes any subset of A.

h(A, B ) ≤ h(A, B )

(2)

In this case, A is the point samples taken from A. As more samples are examined, h(A, B) reaches closer to the HD. We call this the lower bound, denoted by h. As long as the number of samples is finite, h is easily computable, since we only need to perform point projection on B, i.e., finding the minimum distance from a point to B. For this purpose, we utilize the uniform grid data structure [10] that assigns every triangle in B to all cells it crosses. This allows us to quickly search only a few number of triangles. The method fails when there are no triangles nearby, but given our circumstances (near-zero HD), it rarely happens. Yet we must fall back to BVH-based point projection [11] in such cases.

(3)

There is another property of the HD, where B is a subset of B.

h(A, B ) ≤ h(A, B )

(4)

Applying (4) to (3) gives

h(A, B ) ≤ max h(Ak , Bk ), k

(5)

where Bk is Ak ’s match that we have chosen on B. If we can find a local upper bound hk for h(Ak , Bk ), their maximum will act as the global upper bound h. The local upper bound is used in two major ways. Firstly, by continuing to split A into smaller pieces, we gradually obtain smaller hk , also bringing down h closer to the HD. Secondly, any piece with hk that is less than the current known lower bound h is guaranteed to not contain any point where the HD occurs. Therefore, it can safely be removed from the list of candidates, meaning that it does not have to be split any more. It also serves as guidance for point sampling; we do not have to take samples from discarded pieces. While it is preferred to have a match whose local upper bound is as small as possible, it is also important that the calculations are simple, to minimize the HD computation time. So it is key to note that points on Ak are very rarely assigned to the actual closest point on B. To summarize, we compute the HD by iteratively narrowing down the lower bound h and the upper bound h until the difference between the two is less than the desired error bound. First, h is initialized by projecting all vertices of the quad mesh A to the triangle mesh B. Then, we calculate the local upper bound hk for every quad Ak . The ones that satisfy hk < h are thrown away, and the rest are all collected in a queue that keeps the piece with the greatest hk on top. This is done to examine first the piece with high probability of updating h; quick update of h leads to quick elimination of pieces. Now we begin the main iterative process. Every iteration, a piece Ak is popped from the queue, which reduces the h to a smaller or equal value. Ak is split into four smaller pieces, which involves evaluating new points on Ak . These points are projected onto B, updating h. We find the local upper bound for the four new pieces and compare them to h to judge whether to push each into the queue. The cycle is repeated until the queue is empty or the error bound is reached. The exact details of splitting, matching and upper bounding will be discussed in the next section.

4. Matching and upper bounding Since the four vertices of a quadrilateral are not guaranteed to lie on the same plane, we cannot assume the quad to be a planar surface either. While the usual solution is to treat one as two triangles sharing an edge, we also consider an alternative where every quad is a bilinear surface connecting its four vertices. Matching has to be done differently depending on the interpretation. One common element is that in the beginning, every quad is split into two triangular pieces for enqueueing. Then during the iteration, a popped piece is again split into four triangular pieces. This is straightforward for the first interpretation—a quad is divided along a predetermined diagonal into two triangles, and then

Y. Kang, S.-H. Yoon and M.-H. Kyung et al. / Computers & Graphics 81 (2019) 61–72

each triangle is repeatedly split into four identical ones, by inserting new vertices at the edges’ midpoints. (These vertices are projected to B to update h.) But for bilinear surfaces, we do something similar to reduce the number of vertices of each piece from four to three. Since we will analyze their projections to determine the matching, we wish to simplify the cases. An analytic description may be necessary for the splitting of bilinear surfaces. Calling its four vertices a, b, c, and d, a quad’s bilinear surface can be defined as

B(u, v ) = a(1 − u )v + b(1 − u )(1 − v ) + cu(1 − v ) + duv,

(6)

where 0 ≤ u, v ≤ 1, so that a = B(0, 1 ), b = B(0, 0 ), c = B(1, 0 ), and d = B(1, 1 ) are satisfied. This can be split into two triangular Bézier surfaces S1 and S2 , by cutting along either of the two diagonals in the parametric domain: u = 1 − v or u = v. To avoid elongated pieces, we choose u = 1 − v if a − c < b − d and u = v otherwise. Splitting along u = 1 − v can be regarded as a reparameterization of B(u, v ). For S1 , we wish to have S1 (1, 0, 0 ) = a, S1 (0, 1, 0 ) = b, and S1 (0, 0, 1 ) = c. We can achieve this by substituting u, 1 − u, v, and 1 − v with w¯ , u¯ + v¯ , u¯ , and v¯ + w¯ , respectively.

B(u, v ) = a(1 − u )v + b(1 − u )(1 − v ) + cu(1 − v ) + duv

63

Fig. 1. If the three corners of Ak are projected onto two triangles sharing an edge, we split it into a triangle and a quadrangle for separate matching.

Fig. 2. One of the three quadrilaterals that form Ak is matched to another quadrilateral on T3 , based on the projections on T3 of the four corners. The projections of the centroid and midpoints are not necessarily the true projections onto Bk .

= a(u¯ + v¯ )u¯ + b(u¯ + v¯ )(v¯ + w¯ ) + cw¯ (v¯ + w¯ ) + dw¯ u¯ = au¯ 2 + bv¯ 2 + cw¯ 2 + (a + b )u¯ v¯ + (b + c )v¯ w¯ + (b + d )w¯ u¯ = S1 (u¯ , v¯ , w¯ ).

(7)

In other words, S1 is a quadratic Bézier triangle defined by the six c d control points a, b, c, a+2 b , b+ , and b+ . Similarly, S2 is defined by 2 2

a d the control points c, d, a, c+2 d , d+ , and b+ . Notice that the new 2 2 control points are the midpoints of the quad’s four edges and the diagonal bd. If we were to choose u = v as the cut, we would need the midpoint of ac. When a Bézier triangle S(u, v, w ) needs to be split, we cut along three lines u = 0.5, v = 0.5, and w = 0.5, resulting in four smaller Bézier triangles. This requires the evaluation of S(0.5, 0.5, 0), S(0, 0.5, 0.5), and S(0.5, 0, 0.5), which we can project onto B to update h, using the BVH introduced by Garland et al. [11] and the uniform grid.

4.1. Quads as adjoining triangles Let pi (i = 1, 2, 3) be the three vertices of Ak . To update h, we have already projected them onto B. Based on the relationship between the triangles to which the footpoints ri belong, we categorize Ak into one of three cases and match it accordingly. a. The most basic desirable case is when the three footpoints are on the same triangle T. We match Ak = p1 p2 p3 with r1 r2 r3 , which is a part of B. This means we linearly assign every point on Ak to a point on the match so that pi is assigned to ri . Then the local upper bound can be determined directly from the vertices and footpoints.

ha (p1 p2 p3 , T ) = max pi − ri  i

(8)

In this case, this is also known to be the exact h(Ak , Bk ). b. In many cases, the three vertices may be projected onto two triangles T1 and T2 that share an edge. Without loss of generality, let’s say that r1 and r2 are on T1 , and r3 is on T2 , as in Fig. 1. We choose a point b1 on p1 p3 and b2 on p2 p3 to split Ak into p3 b1 b2 and p1 p2 b2 b1 , so that points on the triangle are roughly closer to T2 than T1 , and points on the quadrangle closer to T1 than T2 ; we are trying to match them

separately. For b1 and b2 , we consider the bisector plane of the two planes containing T1 and T2 , not the correct bisector surface of T1 and T2 themselves. The intersections of the plane and p1 p3 and p2 p3 are set as the new vertices b1 and b2 . In the very rare case where the bisector plane does not meet Ak , we simply choose b1 and b2 to be the midpoints of the edges. We project bj onto T1 and T2 to get footpoints rj1 and rj2 . Matching p3 b1 b2 with r3 r12 r22 , and p1 p2 b2 b1 with r1 r2 r21 r11 , the local upper bound can be set as

hb (p1 p2 p3 , T1 , T2 ) = max(max pi − ri , i

b1 − r11 , b1 − r12 , b2 − r21 , b2 − r22  ).

(9)

c. The last is the most general “fallback” case. ri ’s fall on different triangles Ti ’s, although two of them may be the same. We consider many options for matching and choose the best one—the benefit of minimizing the local upper bound far outweighs the computational cost of exploring the options. The first option is to split Ak into three quadrilaterals by cutting it along the lines cm1 , cm2 , and cm3 , where mi is the midpoint of pi pi+1 and c is the centroid of Ak . We wish to match each pi mi cmi+2 with a quadrilateral on Ti , so we project c, mi , and mi+2 onto Ti , obtaining footpoints ric , ria , and rib ; ri ria ric rib is the match (Fig. 2). We bound each quadrilateral’s HD and take the maximum as a possible local upper bound for Ak .

hc (p1 p2 p3 , T1 , T2 , T3 ) = max max(pi − ri , mi − ria , c − ric , mi+2 − rib  ) i

(10) But we also evaluate the option of matching Ak as a whole with a subtriangle on one of the three Ti ’s, even if the true projection of pi+1 and pi+2 are not on Ti . In such cases, we independently project them onto Ti to calculate ha . (From now on, such occurrences will go without mention.) The final local upper bound becomes the minimum:

hc (p1 p2 p3 , T1 , T2 , T3 )





= min hc (p1 p2 p3 , T1 , T2 , T3 ), min ha (p1 p2 p3 , Ti ) . i

(11)

64

Y. Kang, S.-H. Yoon and M.-H. Kyung et al. / Computers & Graphics 81 (2019) 61–72

4.2. Quads as bilinear surfaces Since we are dealing with quadratic Bézier triangles with three vertices, we can use the same classification of cases as in Section 4.1, although the matching and upper bounding are fairly different. Let qi ’s be the control points that are not pi ’s, i.e., the three vertices, meaning that Ak can be expressed as the following, where 0 ≤ u, v, w ≤ 1 and u + v + w = 1.

(12)

Fig. 3. We can give one-to-one correspondence between Bézier triangle Ak and triangle Bk and quickly bound the norm of the difference using the control points.

a. If all ri ’s are on the same triangle T, we match Ak with r1 r2 r3 . To find a good local upper bound, we represent Bk (a subtriangle of T) as a Bézier triangle.

Bk (u, v, w ) = r1 u + r2 v + r3 w = (r1 u + r2 v + r3 w )(u + v + w ) = r1 u2 + r2 v2 + r3 w2 + ( r1 + r2 )uv + (r2 + r3 )vw + (r3 + r1 )wu

(13) r +r

Bk ’s control points become ri ’s and hi ’s, where hi = i 2i+1 . With the convex hull property of Bézier surfaces, we can bound the difference between Ak and Bk by simply taking the maximum among the distances between each of the six control point pairs (Fig. 3). This upper bound will be reused in other cases.

(14)

b. When pi ’s are projected onto two neighboring triangles T1 and T2 , we split Ak into three Bézier triangles (Fig. 4). Let us once again assume that only r3 is on T2 . For some s and t, we cut Ak along two lines in the parameter domain: bs bt and bs b2 , where bs = (s, 0, 1 − s ), bt = (0, t, 1 − t ), and b2 = (0, 1, 0 ). We ¯ ¯ are subdividing two boundary curves p and p 1 p3 2 p3 by the ratios 1 − s : s and 1 − t : t. Without diving into detail, we obtain the following Bézier triangles, each with six control points.

Fig. 4. Ak is subdivided into three Bézier triangles along lines in the parameter domain for separate matching.

We apply ha for each Si and take the maximum.

(17)

c. For the fallback case that can be applied to any situation, we consider three subtriangles of Ak that cover it and overlap; S1 takes the subsurface defined over the range 13 ≤ u ≤ 1, while S2 and S3 do the same for v and w.

(18) (15)

where

pia = where

pib =

ps = p3 ( 1 − s )2 + 2q3 s ( 1 − s ) + p1 s2

qia =

pt = p3 (1 − t )2 + 2q2 t (1 − t ) + p2 t 2 q 1s = q 3 ( 1 − s ) + p1 s

qib =

q 2s = q 2 ( 1 − s ) + q 1 s

qic =

q 3s = p3 ( 1 − s ) + q 3 s q2t = q2 (1 − t ) + p2 t q3t = p3 (1 − t ) + q2 t qst = p3 (1 − s )(1 − t ) + q2 (1 − s )t + q3 s(1 − t ) + q1 st. As for s and t, since we hope ps and pt to be close to the bisector of T1 and T2 , we first measure the distance di from pi to the bisector plane of the two planes containing T1 and T2 . Then we set s and t to be

d3 s= , d1 + d3

d3 t= . d2 + d3

(16)

1 pi + 9 1 pi + 9 1 pi + 3 1 pi + 3 1 pi + 9

4 4 qi + pi+1 9 9 4 4 qi+2 + pi+2 9 9 2 qi 3 2 qi+2 3 2 2 4 qi + qi+2 + qi+1 . 9 9 9

We match Si to Ti and take the maximum of the resulting upper bounds. But similarly to the adjoining triangles case, we also consider matching Ak as a whole to a Ti . This is another scenario where the computational cost is justified by the immense reduction of the number of pieces pushed into the queue.

(19)

Y. Kang, S.-H. Yoon and M.-H. Kyung et al. / Computers & Graphics 81 (2019) 61–72

65

Fig. 5. A visualization of the two-sided HD computation progress for the bunny head example where each quad is a bilinear surface. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Fig. 6. A visualization of the two-sided HD computation progress for the gargoyle example where each quad is a bilinear surface. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

5. Two-sided HD and experimental results The full two-sided HD between a triangle mesh and a quad mesh is the maximum of the two one-sided HDs discussed in Kang et al. [5] and this paper. So while consecutive execution of the two algorithms is a valid option for computing the two-sided HD, blending the two can achieve faster performance, because one of the two directions will benefit from having an increased lower bound, removing candidate pieces even more effectively. The lower bound is initialized by projecting the vertices of the triangle mesh to the quad mesh, and also the vertices of the quad mesh to the triangle mesh. Pieces of both meshes are pushed to the same queue. Depending on whether the popped piece belongs to the triangle mesh or the quad mesh, we apply the appropriate matching and local upper bounding. (An example will be explored later with Fig. 8.) Thus, we have tested the performance of our algorithm for both the one-sided HD from quad mesh to triangle mesh (left halves of Tables A.3–A.6) and the two-sided HD (right halves), computed using the merged algorithm. All tables can be found in the appendix. Table A.1 lists the models used along with the number of vertices and faces. They are the same ten pairs of models as Kang et al. [5], sorted into two groups for tabulating purposes, with the exception of the sphere models, which have been replaced with a more challenging version. The two new sphere models share the same tessellation, where π is divided into 50 segments, but one π . Naturally, the one-sided HD is rotated along the polar axis by 100 in either direction is the same. Other models are all pairs of a triangle mesh and a quad mesh that approximates it. Five pairs are retrieved from the AIM@SHAPE-VISIONAIR Shape Repository, quad mesh approximations using the PGP method by Ray et al. [12]. They are categorized as Group 2 and listed in the latter half of Table A.1, with the rest being Group 1. Tables A.3 and A.5 list

the results for Group 1, while A.4 and A.6 list for Group 2. For Tables A.3 and A.4, each quad is treated as two adjoining triangles, while for A.5 and A.6, quads are regarded as bilinear surfaces. Table A.2 covers triangle mesh-to-quad mesh HD computation for the new sphere models, that is missing from Kang et al. [5]. For each configuration, we compute the HD with a very high precision and list some statistics as we follow the gradual progress of the algorithm. The HD and error bound are all listed as relative values compared to the diagonal length of the models’ bounding box. The first column lists the achieved error bound in logarithmic steps. Every time h − h first reaches 10−l , where l is a natural number, we show (in order) the known lower and upper bounds at the moment, the time taken in milliseconds, the number of iterations taken, and the size of the queue at the time. The first row always shows 0 iterations, so the time listed represents the time it took to project all vertices and sort through all faces once. We also list the preprocessing time required for uniform grid generation, and the maximum queue size on the top. All times are measured as an average taken from a 100 executions on an Intel i7-8700 machine. But before we look into the tables, let us look over what is not on the tables and examine some interesting characteristics. Figs. 5 and 6 visualize how parts of the input meshes are discarded during HD computation, specifically, the two-sided HD computation for the bilinear surface case. The triangle mesh is rendered in red, while the quad mesh is rendered in blue. There are some noticeable differences between the two examples. Even before iteration begins, most of the two gargoyle meshes, especially the quad mesh, has been eliminated, which is not the case for the bunny head example, which retains much of the input halfway through the iteration. This means that the gargoyle example shows relatively larger variation in how much points on one mesh are far from (or close to) the other mesh; points on the quad mesh are close to the triangle mesh, but not vice versa. Also, the yellow dots

66

Y. Kang, S.-H. Yoon and M.-H. Kyung et al. / Computers & Graphics 81 (2019) 61–72

Fig. 7. Quad mesh-to-triangle mesh HD computation progress for the bunny head example where each quad is a bilinear surface.

show where the last known lower bound has occurred, while the green dots mark the final position of the reported HD in advance. For the bunny head example, the lower bound is not updated for nearly the entire duration of the algorithm, only to increase a few times in the end. But sometimes, the yellow dot jumps around as the algorithm discovers greater lower bounds. To be clear, if very sharp features such as spikes are present in the meshes, the distance distribution would have peaks much more prominent than the gargoyle example. Even if a large number of spikes all have very similar distance from their target mesh, spikes actually tend to allow the algorithm to easily localize the candidate area and more quickly find the HD. Such cases are therefore not suitable for detailed demonstration of topics other than that fact. Meanwhile, since all elements in the meshes are tested individually, inverted elements and self-intersections do not impact the algorithm in any way. Now, the very first table from Table A.5, minus the last two columns, is plotted in Fig. 7 to demonstrate the typical HD computation progress. The rows are placed along the x-axis, so the algorithm proceeds from left to right. The known bounds for the HD at each stage are represented as green lines, placed on a logarithmic scale. Observe that they approach each other, sandwiching the true HD. The red line shows the time it took to reach each error bound. Our algorithm calculates the HD within an interactive time, a few seconds at worst for the weightiest models, comparable to the results of Kang et al. [5]. And although we have limited the listing up to 10−8 , it can be noted that after some bottleneck, nearly infi-

nite precision can be reached in almost no time, only hindered by machine precision, which is a meaningful and surprising result. The performance gain of the merged two-sided HD algorithm can be seen when we add the time/iteration measurements of the left-hand side tables to those of the tables from Kang et al. [5] (or Table A.2), and compare them to those of the right-hand side tables. The number of iterations taken in the end for the merged algorithm is always less than the number it would have taken if each one-sided HD were computed separately, the only exception being the sphere example, where the two lower bounds are the same, in which case the number of iterations remain unaffected. Fig. 8, converted from the gargoyle tables of Table A.6, along with the corresponding triangle mesh-to-quad mesh HD (h(B, A)) table from Kang et al. [5], illustrates one example. The initial lower bound for h(B, A), 0.0084811 ≈ 10−2.07 , is greater than that for h(A, B), 0.0 0 09520 ≈ 10−3.02 . That value is also used as the initial lower bound for two-sided HD, so some of the pieces of A that are not immediately removed during the computation of h(A, B) are eliminated instantly for the two-sided HD computation, reducing the latter’s computation time. On the other hand, since h(A, B) is greater than h(B, A), which is why the two-sided HD still converges toward h(A, B), some B pieces are certainly discarded earlier during the two-sided HD computation than they are during h(B, A) computation. Fig. 9 depicts the time comparison more clearly for the sphere, John, and bust examples. Quads are interpreted as two adjoined triangles in the sphere example, while each quad is a bilinear surface in the other two. The figure traces only the time column from the tables, but juxtaposes the three HD computation scenarios, along with a theoretical fourth one. The purple and yellow lines show the time for h(A, B) and h(B, A) computation, respectively. The solid red line is for the merged two-sided HD algorithm, while the dotted line is simply the sum of the first two lines, representing the time it would take if the two one-sided HD algorithms are performed one after the other. The latter two examples portray typical time-saving of the merged algorithm. The advantage is especially noticeable for the John example, where the computation time has decreased by over 35% by merging the two algorithms. This is due to the relatively huge difference between the two one-sided HDs, leaving much to gain. However, for the sphere example, where all lower bounds are equal throughout, the advantage disappears, and the computation time has rather increased, likely owing to the management of a larger queue. Nevertheless, the loss is minimal, and it would be wise to use the merged algorithm when the HDs are unknown beforehand. Meanwhile, the memory requirement is tame, with the queue size only reaching six digits. An interesting observation is that for seven model pairs, the queue for the quad mesh-to-triangle mesh HD computation starts much larger than the triangle mesh-to-quad

Fig. 8. HD computation progress for the gargoyle example where each quad is a bilinear surface.

Y. Kang, S.-H. Yoon and M.-H. Kyung et al. / Computers & Graphics 81 (2019) 61–72

67

Fig. 9. HD computation time comparison for three examples. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

mesh HD or the two-sided HD case and does not shrink easily. This is because the initial h is exceptionally small, since—the quad mesh being an approximation of the triangle mesh—the vertices of the quad mesh were generated very close to the triangle mesh. In fact, a vertex being the occurrence of the HD, denoted by boldfaced model names in the tables, happens for none of the nine quad mesh approximations. However, regardless of the queue size, the algorithm still returns in a quick manner, because only a relatively small number of iterations are required. While it is possible to check and purge the unnecessarily inserted pieces, it would only add to the computation time and rather cause slowdown.

Appendix A. Experimental results

Table A.1 The number of vertices and faces of our example models.

6. Conclusion We have introduced the first non-GPU interactive-speed algorithm for the computation of the near-zero Hausdorff distance from quad mesh to triangle mesh, which is robustly guaranteed within any desired error bound down to machine precision. Combined with the previous result of Kang et al. [5], complete two-sided Hausdorff distance computation is possible, paving the way for fast evaluation of quad mesh approximations of triangle meshes.

Table A.2 Triangle mesh-to-quad mesh HD computation for the sphere models. (a) Each quad as two triangles.

Declaration of interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. The authors declare the following financial interests/personal relationships which may be considered as potential competing interests.

Spheres 10

−l

10−2 10−3 10−4 10−5 10−6 10−7 10−8

Prep(ms) 18

Max 18951

h

h

Time

Iter

Size

0.0 0 01426 0.0 0 01426 0.0 0 01426 0.0 0 01426 0.0 0 01426 0.0 0 01426 0.0 0 01426

0.0090776 0.0011424 0.0 0 02422 0.0 0 01510 0.0 0 01430 0.0 0 01427 0.0 0 01426

11 72 114 149 159 161 162

0 13300 23048 32302 34802 35402 35602

9600 18800 11851 2900 800 200 0

(b) Each quad as a bilinear surface.

Acknowledgments We would like to thank our anonymous reviewers for their precious feedback regarding irregular cases and the illustration of the algorithm. This work was funded in part by the MSIT/IITP of Korea (No. 2017-0-00367), and in part by the National Research Foundation of Korea (No. NRF-2018R1D1A1B07048036 and NRF2019R1A2C1003490).

Spheres 10

−l

10−1 10−2 10−3 10−4 10−5 10−6 10−7 10−8

Prep(ms) 40

Max 140 0 0

h

h

Time

Iter

Size

0.0 0 01426 0.0 0 01426 0.0 0 01426 0.0 0 01426 0.0 0 01426 0.0 0 01426 0.0 0 01426 0.0 0 01426

0.0124781 0.0085280 0.0010670 0.0 0 02414 0.0 0 01523 0.0 0 01435 0.0 0 01427 0.0 0 01426

19 108 432 665 761 789 797 863

0 9200 49200 76800 88400 91600 92600 1010 0 0

9600 13800 13600 11600 40 0 0 1400 400 0

68

Y. Kang, S.-H. Yoon and M.-H. Kyung et al. / Computers & Graphics 81 (2019) 61–72

Table A.3 HD computation progress for models in Group 1 where each quad is two adjoining triangles.

Y. Kang, S.-H. Yoon and M.-H. Kyung et al. / Computers & Graphics 81 (2019) 61–72

Table A.4 HD computation progress for models in Group 2 where each quad is two adjoining triangles.

69

70

Y. Kang, S.-H. Yoon and M.-H. Kyung et al. / Computers & Graphics 81 (2019) 61–72

Table A.5 HD computation progress for models in Group 1 where each quad is a bilinear surface.

Y. Kang, S.-H. Yoon and M.-H. Kyung et al. / Computers & Graphics 81 (2019) 61–72

Table A.6 HD computation progress for models in Group 2 where each quad is a bilinear surface.

71

72

Y. Kang, S.-H. Yoon and M.-H. Kyung et al. / Computers & Graphics 81 (2019) 61–72

References [1] Cottrell JA, Hughes TJR, Bazilevs Y. Isogeometric analysis: toward integration of CAD and FEA. John Wiley & Sons; 2009. [2] Bommes D, Lévy B, Pietroni N, Puppo E, Silva C, Tarini M, et al. Quad-mesh generation and processing: a survey. Comput Gr Forum 2013;32(6):51–76. [3] Hanniel I, Krishnamurthy A, McMains S. Computing the Hausdorff distance between NURBS surfaces using numerical iteration on the GPU. Gr Models 2012;74(4):255–64. [4] Kim Y-J, Oh Y-T, Yoon S-H, Kim M-S, Elber G. Efficient Hausdorff distance computation for freeform geometric models in close proximity. Comput Aided Des 2013;45(2):270–6. [5] Kang Y, Kyung M-H, Yoon S-H, Kim M-S. Fast and robust Hausdorff distance computation from triangle mesh to quad mesh in near-zero cases. Comput Aided Geom Des 2018;62:91–103. [6] Cignoni P, Rocchini C, Scopigno R. Metro: measuring error on simplified surfaces. Comput Gr Forum 1998;17(2):167–74.

[7] Aspert N, Santa-Cruz D, Ebrahimi T. MESH: measuring errors between surfaces using the Hausdorff distance. In: Proceedings of the IEEE International Conference on Multimedia and Expo. IEEE; 2002. p. 705–8. [8] Tang M, Lee M, Kim YJ. Interactive Hausdorff distance computation for general polygonal models. ACM Trans Gr 2009;28(3) 74:1–74:10. [9] Bartonˇ M, Hanniel I, Elber G, Kim M-S. Precise Hausdorff distance computation between polygonal meshes. Comput Aided Geom Des 2010;27(8):580–91. [10] Franklin WR. Nearest point query on 184M points in E3 with a uniform grid.. In: Proceedings of the seventeenth Canadian conference on computational geometry; 2005. p. 239–42. [11] Garland M, Willmott A, Heckbert PS. Hierarchical face clustering on polygonal surfaces. In: Proceedings of the 2001 Symposium on Interactive 3D Graphics. ACM; 2001. p. 49–58. [12] Ray N, Li WC, Lévy B, Sheffer A, Alliez P. Periodic global parameterization. ACM Trans Gr 2006;25(4):1460–85.