Efficient Hausdorff Distance computation for freeform geometric models in close proximity

Computer-Aided Design 45 (2013) 270–276

Contents lists available at SciVerse ScienceDirect

Computer-Aided Design journal homepage: www.elsevier.com/locate/cad

Efficient Hausdorff Distance computation for freeform geometric models in close proximity Yong-Joon Kim a , Young-Taek Oh b , Seung-Hyun Yoon c , Myung-Soo Kim a,∗ , Gershon Elber d a

School of Computer Science and Engineering, Seoul National University, Seoul 151-744, Republic of Korea

b

Samsung Advanced Institute of Technology, Giheung-gu, Yongin-si, Gyeonggi 446-712, Republic of Korea

c

Department of Multimedia Engineering, Dongguk University, Seoul 100-715, Republic of Korea

d

Computer Science Department, Technion, Haifa 32000, Israel

article

info

Keywords: Hausdorff Distance NURBS surface Coons patch Bounding volume hierarchy (BVH) Hierarchical data structure Geometric algorithm

abstract We present an interactive-speed algorithm for computing the Hausdorff Distance (HD) between two freeform geometric models represented with NURBS surfaces. The algorithm is based on an effective technique for matching a surface patch from one model to the corresponding nearby surface patch on the other model. To facilitate the matching procedure, we employ a bounding volume hierarchy (BVH) for freeform NURBS surfaces, which provides a hierarchy of Coons patches and bilinear surfaces approximating the NURBS surfaces (Kim et al., 2011 [1]). Comparing the local HD upper bound against a global HD lower bound, we can eliminate the majority of redundant surface patches from further consideration. The resulting algorithm and the associated data structures are considerably simpler than the previous BVH-based HD algorithms. As a result, we can compute the HD of two freeform geometric models efficiently and robustly even when the two models are in close proximity. We demonstrate the effectiveness of our approach using several experimental results. © 2012 Published by Elsevier Ltd

1. Introduction The Hausdorff Distance (HD) between two objects measures the maximum deviation of one object from the other, and vice versa. The zero or small HD means that the two objects are identical or similar to each other. The HD measure is thus an essential tool for shape matching and recognition [2], where the majority of applications require an estimation of the HD between two similar objects in close proximity. Other important applications include shape approximation and simplification [3,4], where the error should be bounded within a small HD. Unfortunately, there has been no known efficient algorithm for measuring an extremely small HD between two almost identical objects, overlapping each other everywhere, in particular, when they are non-convex 3D geometric models. In this paper, we address this important issue and propose a simple solution to this non-trivial geometric problem by matching surface patches in close proximity. Fig. 1 shows snapshots from the HD computation between two Utah teapots (shown in red and blue), where the red teapot is translated/rotated from the blue teapot along/about the x, y, zdirections by a small amount. Even though the two teapots overlap each other everywhere, our algorithm can compute their small HD

∗

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

0010-4485/$ – see front matter © 2012 Published by Elsevier Ltd doi:10.1016/j.cad.2012.10.010

in an interactive speed, which is a significant improvement over the conventional algorithms [5–8,3,9,10] including those accelerating the HD computation using GPU. In a recent work, Hanniel et al. [11] independently developed a GPU-based HD algorithm that can handle the near-overlap and near-offset cases, which had been the failure cases in the previous methods. Even though we assume no GPU-acceleration, our interactive-speed HD algorithm can also deal with these failure cases reliably in the current CPU implementation. An important advantage of CPU-basd algorithms is the applicability in a broader class of computing environments including those for mobile computing where resources are severely limited in power, space, bandwidth, etc. The performance improvement of our HD algorithm is based on an effective technique for matching the surface patches in close proximity and then efficiently bounding their maximum distance. For this purpose, we employ the recent result of Kim et al. [1], which generates a hierarchy of simple surfaces (such as Coons patches, bilinear surfaces, and tetrahedra) approximating the given NURBS surfaces in a progressive manner. Fig. 2 illustrates a simple technique for matching a surface patch Ai to a nearby surface patch Bi of similar size from the other surface B. In Fig. 2(a), the four corners of Ai are projected to the nearest Bézier surface B(s, t ). Using the surface parameters of the four footpoints on the surface B, we can extract a surface patch Bi of

Y.-J. Kim et al. / Computer-Aided Design 45 (2013) 270–276

271

Fig. 1. Two Utah teapots overlap each other everywhere: (a)–(c) the red teapot is translated from the blue teapot along the x, y, z-directions by a short distance (1% of the model size), and (d)–(f) the red teapot is rotated about the x, y, z-axis by 0.01 radian. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

a

b

Fig. 2. Matching two nearby surface patches: (a) the four corners of Ai are projected to the other surface B, and (b) a surface patch Bi of similar size is extracted from B by bilinearly reparameterizing Bi (u, v) = B(s(u, v), t (u, v)).

similar size to that of Ai by bilinearly reparameterizing Bi (u, v) = B(s(u, v), t (u, v)). The HD between Ai and Bi can be bounded from above by maxu,v ∥Ai (u, v) − Bi (u, v)∥. Based on the surface matching, we maintain a simple priority queue of surface patches (instead of pairs of surface patches), which plays a crucial role in avoiding the memory blow-up usually observed in the conventional HD computation algorithms for two almost identical objects. In computing the one-sided HD from A = ∪Ai to B = ∪Bj , the conventional algorithms maintain a priority queue of node pairs (Ai , Bj ) that have survived upto the current stage of the HD algorithm. When the HD is very small, the conventional algorithms generate many node pairs (Ai , Bj ) that compete each other for the global HD (as the local HDs are about the same as the global HD). As the HD algorithms proceed, the priority queue will blow-up with an exponentially growing number of node pairs. To avoid this problem, we simplify the management scheme by keeping only Ai ’s and Bj ’s in the priority queue and efficiently find their matching Bi ’s and Aj ’s as needed on the fly. Note that we combine the projections from A to B and those from B to A in a single stream of jobs to be processed, which can be effectively managed in a single priority queue for two different types of surface patches Ai ’s and Bj ’s. Instead of computing one-side HD twice and taking the larger one as the final result, we compute the two-sided HD directly. Moreover, this approach is more efficient as the surface patches from the larger one-sided HD will take higher priority in the queue.

Our main contribution can be summarized as follows: – We present an interactive algorithm for computing the HD between two freeform geometric models represented with NURBS surfaces. – Our HD algorithm works efficiently and robustly even for the most challenging case where two almost identical models overlap each other in close proximity. – Using a hierarchy of Coons patches approximating the NURBS surfaces, we have greatly simplified the matching between two nearby freeform surfaces and thus accelerated the elimination of many redundant surface patches from further consideration. The rest of this paper is organized as follows. In Section 2, we briefly review the previous work on the HD computation. In Section 3, we discuss how to bound the HD for two matching surface patches. Section 4 presents the technical details of our HD algorithm. In Section 5, we show the performance of our algorithm using several experimental results. Finally, in Section 6, we conclude this paper. 2. Related work The Hausdorff Distance (HD) computation for polygonal meshes had been a very difficult task to implement for real-time performance [5,7,3]. Using the BVH of polygonal meshes, Tang et al. [10] presented the first real-time algorithm for computing the HD of polygonal meshes. Recently, Krishnamurthy et al. [9] developed a GPU-accelerated real-time HD algorithm for dynamically deformable NURBS surfaces. Nevertheless, it is still a challenging task to develop a real-time HD algorithm for NURBS surfaces without assuming GPU-acceleration. Many applications in computer graphics and geometric modeling require the HD computation between two similar 3D geometric models in close proximity. For example, starting with a dense mesh, simplified meshes are generated for improving the performance of visualization, where the HD computation is needed between two similar meshes in close proximity. The conventional HD algorithms are slow even though they have been developed mainly for the special purpose of measuring a small HD between two similar meshes. Even the recent real-time HD algorithms [9,10] also experience significant reduction in their performance when applied

272

Y.-J. Kim et al. / Computer-Aided Design 45 (2013) 270–276

to two almost identical models overlapping each other in close proximity. In the rest of this section, we review the previous HD algorithms in the respect of their suitability for handling the highly overlapping cases. Early work in computer vision considered efficient HD algorithms for images [12]. Because of the simple data structure of images, it is easy to measure the difference between two almost identical images. Nevertheless, the image resolution is usually insufficient for geometric modeling applications. Moreover, it is rather difficult to extend this approach to 3D because of the considerably larger memory space needed in the higher dimension. The linear-time HD algorithm of Atallah [13] is restricted to the case of two non-overlapping planar convex polygons. The majority of conventional HD algorithms for polygonal meshes approximate the HD between two similar meshes [5,7,3]. Barton et al. [6] compute the precise HD between polygonal meshes. Their performance was far from real-time, in particular since they assume neither GPU-acceleration nor preprocessing of the input mesh models. There are also some previous results for freeform shapes. Jüttler [14] presented a method to estimate bounds on the Hausdorff distance between two freeform/implicit curves. Alt and Scharf [15,16] developed an algorithm for the precise HD between planar rational curves. Elber and Grandine [8] extended this result to a more general problem for rational curves or surfaces. These results are based on solving a system of multivariate polynomial equations, which takes a considerable amount of computing time when dealing with two almost identical freeform curves and surfaces. Recently, Kim et al. [17] developed a real-time algorithm for the precise HD between planar freeform curves using hardware depth buffer. The algorithm works in a fashion similar to the image-based techniques and thus can deal with two almost identical curves. At the end of redundancy elimination, the algorithm switches to the precise HD computation of [8] or to a more accurate approximation of the planar freeform curves using biarcs. Thus the image resolution is not a bottleneck for the HD computation. Nevertheless, it is rather difficult to extend this approach to freeform space curves or surfaces in 3D. There are also two other recent HD algorithms [18,19] for freeform curves. Mainly based on pruning techniques for redundant curve segments, their performance would be slowed down substantially when dealing with two almost identical curves. In a recent work, Hanniel et al. [11] developed a GPU-accelerated algorithm that can compute many point projections from one surface to the other in parallel using the GPU architecure. The tremendous computing power of GPU enabled the sampling-based approach to handle the near-overlap and near-offset cases successfully. On the other hand, we employ only a CPU implementation and yet deal with the near-overlap and near-offset cases in an interactive-speed as demonstrated in Section 4 for the cases of identical objects translated and rotated under many different conditions and also for the HD computation between concentric spheres of different radii. 3. HD upper bound for matching surfaces A freeform surface S (u, v), 0 ≤ u, v ≤ 1, can be approximated by a Coons patch X (u, v) that interpolates the four boundary curves of S (u, v) [20], which provides a high approximation order to the given surface since the interpolation is exact on the surface boundary. In a recent work, Kim et al. [1] proposed a new BVH for NURBS surfaces using Coons approximation, and demonstrated real-time performance of collision detection and minimum distance computation for freeform geometric models. They subdivide a surface S (u, v) recursively until the maximum error: max(u,v) ∥S (u, v) − X (u, v)∥ is reduced within a given tolerance, and record the hierarchy of recursively subdivided Coons patches as a BVH for the surface S (u, v).

Coons patches have a nice property that the whole bivariate surface is completely determined by the boundary information. Thus, approximating the boundary curves of a Coons patch X (u, v), 0 ≤ u, v ≤ 1, using polylines sampled at ui and vj , the Coons patch X (u, v) can be approximated by a piecewise bilinear surface L(u, v) which is again a Coons patch determined by the four boundary polylines. The HD between X (u, v) and L(u, v) is bounded by ϵu + ϵv , the summation of the approximation errors ϵu and ϵv along the boundary curves in the u and v -directions, respectively. The point projection of the corner points of a surface patch Ai (u, v), 0 ≤ u, v ≤ 1, to the other surface B(s, t ) is a special type of the minimum distance computation, which can be computed very efficiently using the Coons BVH structure of Kim et al. [1]. Using the st-coordinates of the four projection footpoints on the surface B(s, t ), we can reparameterize s = s(u, v), t = t (u, v) as bilinear functions of u and v , and extract a surface patch Bi (u, v) = B(s(u, v), t (u, v)), 0 ≤ u, v ≤ 1, which is about the same size as Ai (u, v) in particular when two almost identical surfaces overlap each other. The HD between the two surface patches Ai and Bi can be bounded from above by the maximum size of their bivariate difference surface: max ∥Ai (u, v) − Bi (u, v)∥.

0≤u,v≤1

The reparameterized surface Bi (u, v) has a degree two times higher than that of B(s, t ). Thus the representation for (Ai − Bi )(u, v) requires a degree elevation for Ai (u, v) and a symbolic composition for Bi (u, v) = B(s(u, v), t (u, v)), which is timeconsuming. In this paper, we take a simple sampling approach for estimating an upper bound of the above maximum. Krishnamurthy et al. [9] also took a sampling approach (based on the condition of Filip et al. [21]) for estimating the approximation error of the given C 2 surfaces S (u, v) using surface points uniformly sampled along the u and v parameters. The approximation error can be bounded using the second partial derivatives of S (u, v). We can speed up the sampling error estimation by precomputing the partial derivatives of Ai (u, v) and B(s, t ) and storing them in the Coons BVH of the surfaces. Only the partial derivatives of s(u, v) and t (u, v) need to be computed on the fly as needed. 4. Hausdorff Distance computation Given two surfaces A and B, the one-sided Hausdorff Distance (HD) from A to B:





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

q∈B

measures the maximum deviation of A from B. The (two-sided) Hausdorff Distance is then defined as H (A, B) = max(h(A, B), h(B, A)). In the following discussion, we assume h(A, B) ≥ h(B, A) and compute h(A, B) = H (A, B). The one-sided Hausdorff Distance h(A, B) satisfies the following relation: h(A, B¯ ) ≤ h(A, B) ≤ h(A¯ , B),

¯ For an estimation of the lower for any A ⊂ A ⊂ A¯ and B ⊂ B ⊂ B. bound h = h(A, B¯ ), we construct the subset A using the midpoints pi , for i = 1, . . . , m, of surface patches Ai of A and the superset j B¯ using the bounding volumes Oϵ j (VB ), for j = 1, . . . , n, of surface B

j

patches Bj of B. (The bounding volume Oϵ j (VB ) is an offset of a tetraB

j

hedron VB determined by the four corners of a bilinear surface and j

the offset radius ϵB is an upper bound for the distance between the bilinear surface and the given NURBS surface [1].)

Y.-J. Kim et al. / Computer-Aided Design 45 (2013) 270–276

For an estimation of the upper bound h¯ = h(A¯ , B), we construct the superset A¯ using the bounding volumes Oϵ i (VAi ), for i = 1, A

. . . , m, of surface patches Ai of A, but the subset B is taken as the

union of the matching surface patches Bi of B to the surface patches Ai of A. There is no guarantee that the union ∪Bi completely covers the surface B even though we start with a complete covering of A, i.e., ∪Ai = A. Estimating HD Lower Bound: Each midpoint pi ∈ Ai is first proj jected to the nearest point qi ∈ ∪j Oϵ j (VB ) on the boundary of the B

bounding volumes of B. We take the current lower bound of HD as h = maxi ∥pi − qi ∥. (The point projection is first drawn to the nearj

est point on the tetrahedron VB . The distance is then subtracted by j B

the offset radius ϵ if it is larger than the offset radius; or the distance is set to 0 otherwise.) By recursively going down to the deeper levels of the Coons BVH structure, Kim et al. [1] compute the approximate point projection j qi ∈ Oϵ j (VB ) so that there is one surface point B(si , ti ) within a B

given tolerance δ (e.g., 0.01% of the model size). Thus we have the relation: d(pi , B) − δ < ∥pi − qi ∥ ≤ d(pi , B), where d(pi , B) is the distance between the point pi and the surface B(s, t ). Matching Surface Patches: We first check whether the distance ∥pi − qi ∥ is larger than a certain tolerance (e.g., 1% of the model size). If it is larger, we test if the surface patch Ai is totally contained in the ball of radius h − δ centered at the approximate footpoint qi . When the HD is relatively large, many redundant patches can be eliminated using this simple containment test with respect to a ball and thus it can speed up the HD computation quite dramatically. Otherwise, we proceed to a more expensive test as discussed below. Using the footpoint information qi , we take the nearest Bézier surface on the other NURBS surface B(s, t ). Then we repeat similar point projections of the four corners of Ai to the Bézier surface B(s, t ), and the surface parameters (s, t ) are bilinearly reparameterized to s = s(u, v) and t = t (u, v). The matching surface patch is generated as Bi (u, v) = B(s(u, v), t (u, v)), 0 ≤ u, v ≤ 1, which has about the same size as that of Ai (u, v). Culling Redundant Patches Ai : Using the technique discussed in the previous section, we estimate an upper bound h¯ i for maxu,v ∥Ai (u, v) − Bi (u, v)∥. When the local upper bound hi is smaller than the current global HD lower bound h, we can eliminate the redundant patch Ai from further consideration. Otherwise, the surface patch Ai is inserted into a priority queue, sorted according to the value of h¯ i . The surface patch Ai with the largest value of h¯ i is then dequeued from the priority queue and subdivided into two pieces. The same procedure is then repeated recursively. Estimating HD Upper Bound: The maximum deviation bound h¯ i for the patch Ai just dequeued from the priority queue is the current HD upper bound h¯ since all the remaining Ak ’s in the priority queue have no larger value than h¯ i . If the difference h¯ − h between the upper and lower HD bounds is less than a given tolerance (e.g., 0.1% of the model size), we terminate the HD computation and make a good guess of the exact HD in the interval [h, h¯ ]. Otherwise, we recursively subdivide the subpatch Ai into two pieces and repeat the same procedure recursively. Guessing the Exact HD: The global HD lower bound h is usually closer to the exact HD distance than the HD upper bound h. This is because the HD lower bound h is computed as the maximum of the distances from sample points pi of A to the footpoints on the bounding volumes of B, whereas the HD upper bound h¯ is taken as the maximum among the over-estimated local HD upper bounds h¯ i . Empirically speaking, a more accurate HD distance can be estimated using the footpoints to the bilinear surfaces since they approximate the NURBS surfaces more tightly than the bounding volumes. On the other hand, theoretically speaking, there is no

273

guarantee that the distance to the bilinear surfaces is even a HD lower bound. Nevertheless, at the end of the HD computation, if the maximum among the distances (from the midpoints pi of Ai ) to the bilinear surfaces (for B(s, t )) is contained in the interval [h, h¯ ], we can take this value as the final HD. Computing Two-sided HD: We cannot assume the relation: h(A, B) ≥ h(B, A) in general, and the previous algorithms repeat the HD computation for h(B, A) one more time. In our approach, the two HD computations for h(A, B) and h(B, A) can be processed in a single stream using the same priority queue for both surface patches Ai and Bj . As the surface patch with the largest HD upper bound is dequeued no matter whether it is Ai or Bj , the basic algorithm works exactly the same way as in the one-side HD computation for h(A, B). 5. Experimental results We have implemented the proposed HD computation algorithm in C++ on an Intel Core i7-2600 3.4 GHz CPU with a 16 GB main memory and an NVIDIA Geforce GTX570. To demonstrate the performance of our approach, we have tested the implemented algorithm for several different freeform geometric models of nontrivial complexity. Tables 1 and 2 report the performance of our algorithm in computing the HD between two identical freeform geometric models, where one model is translated from the other by a given distance (Table 1) or rotated from the other by a given angle (Table 2). In the case of pure translation, the exact HD is the same as the distance of translation. Nevertheless, when we consider rotation, the exact HD is unknown a priori. In the current implementation, we take 0.1% of the model size as the tolerance of terminating the HD computation, i.e., we terminate when the difference h¯ − h between the HD upper and lower bounds is less than 0.1% of the model size. Table 1(a) shows the result for the Utah teapot which is composed of four periodic NURBS surfaces (that contain a total of 160 Bézier surfaces). The first row shows the number of Bézier surfaces in the teapot model. The second row shows the relative distance of the translation with respect to the model size, where the unit length is taken as the largest side length of the minimum bounding box for the teapot and thus 10−3 means 0.1% of the unit length. In this experiment, we take a total of 26 different directions for translating one model from the stationary identical model. The 26 directions are taken as (±1, 0, 0), (0, ±1, 0), (0, 0, ±1), (±1, ±1, 0), (±1, 0, ±1), (0, ±1, ±1), and (±1, ±1, ±1). The third row shows the average of the 26 HD upper bounds thus computed. Similarly, the fifth row shows the average of the 26 HD lower bounds. The fourth row reports the average of the final 26 HDs computed using the footpoints to the bilinear surfaces. (In these three rows, the HD is divided by the translation distance so that the exact HD is rescaled to 1.00; thus the numbers show the relative over and under-estimations of the exact HD. Note that h ∗ 10k = h/10−k , for k = 3, 2, 1, 0.) The last three rows report the time taken in the point projections to the bounding volumes and to the bilinear surfaces, the time taken in the elimination of redundant surface patches, and the total time (in milliseconds) including all the other remaining computations in addition to the two major components. Table 2(a) shows the result for a similar test, but for rotations about the 26 different directions by angle 10−k radian instead of translations, for k = 3, 2, 1, 0. The rest of Tables 1 and 2 reports the results for the three different freeform geometric models shown in Fig. 3(a)–(c). A surface of revolution has singularity at each intersection of the surface with the axis of rotation. Thus in each model of Fig. 3, we made a small hole of size 0.1% of the model size at each singularity. Fig. 3(a) shows two spheres (of the same radius) completely overlapping, but one is rotated from the other. Their

274

Y.-J. Kim et al. / Computer-Aided Design 45 (2013) 270–276

Table 1 HD computation for two identical NURBS surfaces, where one surface is translated along 26 different directions by distance 10−k , for k = 3, 2, 1, 0; the average of 26 HD computation results is reported, where h ∗ 10k means the relative HD with respect to the translation distance 10−k ; and the computing times are reported in miliseconds. (a)

Table 2 HD computation for two identical NURBS surfaces, where one surface is rotated about 26 different directions by angle 10−k radian, for k = 3, 2, 1, 0; the average of 26 HD computation results is reported, where h ∗ 10k means the relative HD with respect to the rotation angle 10−k ; and the computing times are reported in miliseconds. (a)

Teapot

#Bézier(160)

10−3

10−2

10−1

1.00

h ∗ 10 h ∗ 10k h ∗ 10k

1.18 0.39 0.30

0.52 0.43 0.43

0.44 0.43 0.43

0.29 0.29 0.29

Proj. Cull

5 50

11 50

4 49

2 15

Total (ms)

56

62

53

17

10−3

10−2

10−1

1.00

h ∗ 10 h ∗ 10k h ∗ 10k

1.14 0.42 0.30

0.43 0.41 0.38

0.03 0.02 0.02

0.00 0.00 0.00

Proj. Cull

1 22

0 21

6 23

27 24

Total (ms)

24

22

29

51

10−3

10−2

10−1

1.00

h ∗ 10 h ∗ 10k h ∗ 10k

1.05 0.49 0.39

0.53 0.45 0.44

0.46 0.45 0.45

0.32 0.32 0.32

Teapot

10−3

10−2

10−1

1.00

h ∗ 10k h ∗ 10k h ∗ 10k

1.82 1.00 0.95

1.09 1.00 0.99

1.01 1.00 1.00

1.00 1.00 1.00

Proj. Cull

2 47

10 40

19 30

7 4

Total (ms)

49

50

49

11

(b)

k

(b) Sphere

Sphere

#Bézier(8)

k

10−3

10−2

10−1

1.00

h ∗ 10k h ∗ 10k h ∗ 10k

1.68 1.00 0.91

1.04 1.00 0.99

1.00 1.00 1.00

1.00 1.00 1.00

Proj. Cull

2 13

3 24

1 14

1 6

Total (ms)

16

27

16

7

(c)

(c) Bishop k

Bishop

#Bézier(75) 10−3

10−2

10−1

1.00

h ∗ 10 h ∗ 10k h ∗ 10k

1.94 1.00 0.95

1.10 1.00 1.00

1.01 1.00 1.00

1.00 1.00 1.00

Proj. Cull

6 16

6 19

5 12

4 6

Total (ms)

23

25

17

10

Proj. Cull

13 20

17 26

8 17

5 9

(d)

Total (ms)

35

45

26

15

10−3

10−2

10−1

1.00

h ∗ 10 h ∗ 10k h ∗ 10k

1.52 0.59 0.54

0.61 0.52 0.51

0.51 0.50 0.50

0.32 0.32 0.32

Proj. Cull

59 67

69 80

74 45

11 2

Total (ms)

172

189

146

15

k

Bear k

(d) Bear

#Bézier(1393) 10−3

10−2

10−1

1.00

h ∗ 10 h ∗ 10k h ∗ 10k

1.95 1.01 0.96

1.09 1.00 0.99

1.01 1.00 1.00

1.00 1.00 1.00

Proj. Cull

75 54

87 86

86 11

15 1

Total (ms)

172

208

122

19

k

HD is realized at each pole as the distance from the pole to the boundary circle of the missing hole. It is interesting to note that the most difficult case for each test is not always the one for the shortest translation distance or the smallest rotation angle. For models of relatively small size, at the shortest distance 10−3 , the surface matching produces many pairs of similar surface patches to be compared and then immediately discarded as redundant because of their relatively small bounds of the maximum difference. On the other hand, for the same models at a larger distance 10−2 , the surface matching produces many pairs of different shapes and thus the surface patches should be further subdivided and compared. For the tests using different angles of rotation, three models (except the sphere) also experience the worst performance at the rotation angle of 10−2 rad. The rotations of a sphere generate almost identical spheres except the small holes (punctured at the north and south poles) located at various different places. The size of holes becomes relatively smaller as the angle of rotation gets larger; consequently, the relative HD also gets reduced. Though it is unclear from Tables 1 and 2, where each item reports the average of 26 computation results, the HD computation

Fig. 3. HD computation for sphere, bishop, bear, and bunny models represented with NURBS surfaces.

Y.-J. Kim et al. / Computer-Aided Design 45 (2013) 270–276

HD computation for the translation of distance 10−3 along the up/down direction takes considerably more computing than along the front/back direction. Nevertheless, for translations of large distances, the relative difference is not that significant. The eagle model also shows a similar pattern of performance variations. Our HD algorithm is optimized for the case of two almost identical NURBS surfaces in close proximity. Nevertheless, the basic algorithm works reasonably well for the general case of measuring the HD for two different surfaces as well as for two identical surfaces at large distances. Fig. 5 shows the HD computation between a stationary Utah teapot and a moving teapot under a continuous translation and rotation. The HD computation time increases sharply as their HD approaches to zero, reaches the peak at the overlap configuration of the two models, and then immediately drops as the two models are separated. Fig. 6 is a similar scenario but the HD is measured between two different models. Because of the interactions among various sharp features of the two relatively complex models (i.e., the bear and bunny models consisting of 1393 and 1012 Bézier surfaces, respectively), their HD computing time shows a rather fluctuating fashion locally. Nevertheless, the global change in the HD computation time still follows a gross pattern in which the computing time increases as the HD decreases. The HD computation for two concentric spheres of different radii is a challenging task, where the HD is realized almost everywhere and the local HDs are the same as the global HD. Note that our sphere model is not the complete sphere as we have made two small holes at the poles. The global HD for our concentric sphere models is often realized at the projection of a pole onto the other sphere. The global HD is thus slightly larger than the local HD which is the same as the difference between the two sphere radii. Fig. 7 shows the HD computation results for the case of continuously enlarging the radius of one concentric sphere from the unit length

Fig. 4. B58 and eagle models with many flat NURBS surfaces.

for a certain direction of translation/rotation takes longer than for other directions. The relative difference can be quite large for translations in short distance. The translation of distance 10−3 along a direction orthogonal to large flat surfaces takes a much longer computing time than translations to other directions or in large distances. In particular, for the B58 model in Fig. 4, the

a

275

b

Fig. 5. HD computation between a stationary teapot and a moving Utah teapot: (a) one teapot is moving under a continuous translation and rotation, and (b) the HD computation time increases as the HD approaches to zero.

a

b

Fig. 6. HD computation between a stationary bear model and a moving bunny model: (a) the bunny model is moving under a continuous translation and rotation, and (b) the HD computation time increases as their HD decreases.

276

Y.-J. Kim et al. / Computer-Aided Design 45 (2013) 270–276

a

b

Fig. 7. HD computation between two concentric spheres: (a) the spheres share the same axis of rotation, and (b) their axes of rotation make angle 60°.

(which is the radius of the stationary concentric sphere) to twice that. Fig. 7(a) is for the case where the two concentric spheres share the same axis of rotation. Thus the HD is realized between the two circles bounding the missing holes of the two concentric spheres. On the other hand, in Fig. 7(b), the growing sphere has the axis of rotation that makes an angle of 60° with the axis of the stationary sphere. Comparing the two results, we can notice that the case (b) is in general more difficult to compute than the case of (a), in particular when the difference in their radii is small. Moreover, the HD computation becomes more difficult as the HD increases as the relative size of the missing holes gets smaller and thus the problem becomes more and more singular near the poles of the concentric spheres. 6. Conclusions We have presented an efficient algorithm for computing the HD between two freeform NURBS models in close proximity, which is a significant improvement over the previous HD algorithms. Though the algorithm demonstrates an interactive speed, the current implementation employs a relatively large tolerance bound, i.e., 0.1% of the model size. Nevertheless, using footpoints to the bilinear surfaces, we have estimated the exact HD considerably more precisely than the tolerance bound. In future work, we plan to improve the performance of our algorithm by exploring the GPU acceleration of the maximum distance computation for each pair of matching surface patches. Acknowledgments The authors would like to thank Mr. Jeeho Lee for modeling various freeform NURBS models that have been used in our experiments. This research was partly supported by the Israel Science Foundation (Grant No. 346/07), in part by the Israeli Ministry of Science (Grant No. 3-8273), and in part by the NRF of Korea (Grants Nos. 2011-00356, and 2012-0005037, and 2012R1A1A1008908). References [1] Kim Y-J, Oh Y-T, Yoon S-H, Kim M-S, Elber G. Coons BVH for freeform geometric models. ACM Trans Graphics 30(6) Article 169, (SIGGRAPH Asia 2011).

[2] Alt H, Guibas L. Discrete geometric shapes: matching, interpolation, and approximation. In: Handbook of computational geometry. Elsevier Science; 1999. [3] Guthe M, Borodin P, Klein R. Fast and accurate Hausdorff distance calculation between meshes. J WSCG 2005;13:41–8. V. Skala (ed.). [4] Varadhan G, Manocha D. Accurate Minkowski sum approximation of polyhedral models. Graph Models 2006;68(4):343–55. [5] Aspert N, Santa-Cruz D, Ebrahimi T. MESH: measuring errors between surfaces using the Hausdorff distance. In: Proc. of the IEEE int’l conf. in multimedia and expo (ICME). vol. 1. 2002. p. 705–8. [6] Barton M, Hanniel I, Elber G, Kim M-S. Precise Hausdorff distance computation between polygonal meshes. Comput Aided Geom Design 2010;27(8):580–91. [7] Cignoni P, Rocchini C, Scopigno R. Metro: measuring error on simplipied surfaces. Comput Graph Forum 1998;17(2):167–74. [8] Elber G, Grandine T. Hausdorff and minimal distances between parametric freeforms in R2 and R3 . In: Chen F, Jüttler B, editors. Advances in geometric modeling and processing. Lecture notes in computer science, vol. 4975. Springer; 2008. p. 191–204. Procs. of the 5th int’l conf., GMP 2008. [9] Krishnamurthy A, McMains S, Hanniel I. GPU-accelerated Hausdorff distance computation between dynamically deformable NURBS surfaces. Comput Aided Design 2011;43(11):1370–9. [10] Tang M, Lee M, Kim YJ. Interactive Hausdorff distance computation for general polygonal models. ACM Trans Graphics 28(3) Article 74, (SIGGRAPH 2009). [11] Hanniel I, Krishnamurthy A, McMains S. Computing the Hausdorff distance between NURBS surfaces using numerical iteration on the GPU. Graph Models 2012;74(4):255–64. (The special issue of Geometric Modeling and Processing (GMP2012), June 18–22, 2012, Huangshan, China). [12] Rucklidge W. Efficient Visual Recognition using the Hausdorff Distance. Lecture Notes in Comput Sci 1996;1173. [13] Atallah M. A linear time algorithm for the Hausdorff distance between convex polygons. Inform Process Lett 1983;17:207–9. [14] Jüttler B. Bounding the Hausdorff distance of implicitly defined and/or parametric curves. Math Methods in CAGD 2000;1–10. [15] Alt H, Scharf L. Computing the Hausdorff distance between sets of curves. In: Procs of the 20th European workshop on computational geometry (EWCG). 2004. p. 233–6. [16] Alt H, Scharf L. Computing the Hausdorff distance between curved objects. J Comput Geom Appl 2008;18(4):307–20. [17] Kim Y-J, Oh Y-T, Yoon S-H, Kim M-S, Elber G. Precise Hausdorff distance computation for planar freeform curves using biarcs and depth buffer. Vis Comput 2010;26(6–8):1007–16. [18] Bai Y-B, Yong J-H, Liu C-Y, Liu X-M, Meng Y. Polyline approach for approximating the Hausdorff distance between two planar free-form curves. Comput Aided Design 2011;43(6):687–98. [19] Chen X-D, Ma W, Xu G, Paul J-C. Computing the Hausdorff distance between two B-spline curves. Comput Aided Design 2010;42(12):1197–206. [20] Coons S. Surfaces for computer-aided design, Technical report, MIT, 1964. Available as AD 663 504 from the national technical information service. Springfield, VA, 22161. [21] Filip D, Magedson R, Markot R. Surface approximations using bounds on derivatives. Comput Aided Geom Design 1986;3(4):295–311.