Abstract. Point pattern matching is an important problem in computational geometry, with ... algebraic manipulation of high-degree manifolds in a robust fashion.
Combinatorial and Experimental Methods for Approximate Point Pattern Matching∗ Martin Gavrilov
Piotr Indyk
Rajeev Motwani
Suresh Venkatasubramanian
Abstract Point pattern matching is an important problem in computational geometry, with applications in areas like computer vision, object recognition, molecular modelling, and image registration. Traditionally, it has been studied in an exact formulation, where the input point sets are given with arbitrary precision. This leads to algorithms that typically have running times of the order of high degree polynomials, and require robust calculations of intersection points of high degree surfaces. We study approximate point pattern matching, with the goal of developing algorithms that are more efficient and more practical than exact algorithms. Our work is motivated by the observation that in practice, data sets that form instances of pattern matching problems are noisy, and so approximate formulations are more appropriate. We present new and efficient algorithms for approximate point pattern matching in two and three dimensions, based on approximate combinatorial distance bounds on sets of points, and via the use of methods from combinatorial pattern matching. We also present an average case analysis and a detailed empirical study of our methods.
1
Introduction
Given point sets P and Q, the problem of pattern matching is to determine whether P is similar to some portion of Q. Similarity is typically defined in terms of a distance measure, and P and Q may undergo transformations like rotations, translations and scaling before the distance between them is computed. This general formulation captures a wide variety of problems, arising in image registration [33], model-based object recognition [4, 21], computer vision [35], molecular modelling [27], and many other areas. In the domain of computational geometry, these problems were mostly studied in an exact formulation: point coordinates are given with infinite precision, and two points match if their coordinates are identical. The disadvantage with this formulation is two-fold: firstly, this is an inaccurate representation of reality, where the input points are typically noisy, and any match requires some notion of a tolerance parameter. Secondly, the algorithms designed to compute similarity were complex, with running times in the order of large polynomials, and required sophisticated algebraic manipulation of high-degree manifolds in a robust fashion. This provides the motivation for approximate pattern matching. In an approximate setting, one might replace points by balls of radius ² (where ² is a tolerance parameter) and would be willing to compute approximate similarity estimates rather than exact answers. As in most kinds of optimization, the hope is that by doing so one can design algorithms that are significantly faster (and simpler) than the exact counterparts, and thus are more practical as well.
1.1
Our Work
In this paper, we consider problems in approximate point set matching. We make use of two distinct techniques to design algorithms for approximate point set matching. The first is a combinatorial method based on the estimation of the number of approximately identical copies of a set of (two or three) points in a given point set. As we shall see, this quantity can be used to estimate the running times of simple alignment-based algorithms. Our second technique is a novel connection to combinatorial pattern matching (or string matching). Techniques like convolutions have long been used in combinatorial pattern matching to obtain near-linear time algorithms for ∗
This paper expands and combines the contents of papers presented at the 10th ACM-SIAM Conference on Discrete Algorithms [31] and the 15th ACM Symposium on Computational Geometry [23]
various pattern matching problems. We show that geometric matching problems can be reduced to combinatorial pattern matching problems, thus obtaining fast approximations. This idea was independently used by Cardoze and Schulman [8] to provide a near-linear scheme for pattern matching; we discuss their work in relation to ours later. Subsequent to the first presentation of this work, other problems in geometric matching have been solved via a reduction to combinatorial pattern matching [17, 32]. Diameter Dependence It has been observed in computational geometry that the parameter n (defining the number of geometric objects in the input) is often not the most useful parameter to determine the running time of geometric algorithms. Examples of properties of a point set that influence running time include fatness [16] and aspect ratio [42]. The running times of our algorithms contain terms involving ∆, the diameter of a point set. More formally, we denote ∆ to be the ratio of the largest distance between pairs of points to the smallest distance between pairs of points. This quantity, also called the spread of a point set, has been exploited in other work [22] to obtain fast algorithms for various problems via a careful analysis of the dependency on the diameter. Experimental Evaluation We conduct a detailed experimental evaluation of the schemes that we present. The asymptotic bounds presented hide several constants that reflect different degrees of dependence on the tolerance parameters, and thus an empirical study provides deeper insight into the efficacy of our methods in practice. Much of the recent empirical work in geometric pattern matching [29, 25, 33] has examined algorithms based on an exhaustive search of the transformation space (with clever pruning techniques). Our evaluation naturally complements this work. However, an exact comparison of our running times with those obtained in these works is not very meaningful, because the problems being solved are different. Huttenlocher and Rucklidge [29] consider the problem of computing the Hausdorff distance from P to Q under translation (and scaling), and Mount et al. [33] consider the problem of image registration in which the goal is to determine whether a large overlap exists between P and Q. In the case of [25], the distance measure used is different. Outline Of The Paper Section 2 contains basic definitions and we survey related work in Section 3. In Sections 4,5, we discuss combinatorial bounds for approximate distances in planar point sets and present an algorithm 4=3 1=3 ˜ running in time O(kn ∆ ) for matching two point sets P, Q where |P| = k, |Q| = n. Sections 6,7 provide similar results for point sets in three dimensions. In Section 9 we present a reduction from geometric matching to combinatorial pattern matching, yielding an algorithm running in time O(n(n + ∆) for pattern matching in two dimensions. In Section 8 we consider pattern matching under alternate measures and show how our techniques can be used to solve these problems. In the rest of the paper (Sections 10-13), we present a detailed empirical evaluation of the above schemes, and also provide an average case analysis (under probabilistic input models).
2
Definitions and Preliminaries
In this paper, P and Q refer to point sets in Rd , where |P| = k, |Q| = n ≥ k and d = 2, 3. Unless otherwise specified, the underlying metric on points in Rd is that induced by the `2 norm. Let T be a transformation operating on points in Rd , and let T (P) denote the action of T on the set P. Definition 2.1 (Hausdorff distance). Given two point sets P and Q, the Hausdorff distance from P to Q is dH (P, Q) = max min d(p, q) p∈P q∈Q
Definition 2.2 (Pattern Matching (PM)). Given point sets P, Q and space of transformations T dPM (P, Q) = min dH (T (P), Q) T ∈T
2
Definition 2.3 (Congruence (CONG)). Given point sets P and Q and space of transformations T , dC = max(dPM (P, Q), dPM (Q, P)) Definition 2.4 (Largest Common Point Set (LCP)). Given point sets P and Q, space of transformations T , and parameter ², determine A ⊆ P, B ⊆ Q such that min dH (T (A), B) ≤ ²
T ∈T
and |A| is maximized. Definition 2.5 (β-approximate pattern matching PM(², β) [26]). Given ², β > 0 and point sets P and Q, • If ²∗ ≤ ², then return a transformation T ∈ T such that dH (T (P), Q) ≤ (1 + β)²; • If ²∗ > (1 + β)², then return NONE. • Otherwise, when ²∗ ∈ (², (1 + β)²], return any transformation T ∈ T . In the above definition, ²∗ = dPM (P, Q). In essence, the above formulation captures the intuition that determining an answer close to optimal is harder than determining an answer that is further away. Hence, if ²∗ lies in the given approximation range (², (1 + β)²], the algorithm is not guaranteed to return a near-optimal answer, but if not, it gives an answer close to the prescribed error bound ². In general, one can define the following promise problem: Definition 2.6. Given sets P and Q, approximation parameters α and β, distance measure d, space of transformations T , and numbers l, r, find T ∈ T such that 1. There exists A ⊆ P, B ⊆ Q such that d(T (A), B) ≤ (1 + β)r 2. |A| ≥ (1 − α)l Here l and r are the promise parameters: we are guaranteed that there exists a subset A ⊆ P of size l that matches some B ⊆ Q to within distance r, and we wish to approximate this solution. We call a solution to this problem an (α, β)-approximation. Heffernan and Schirra provided a very useful technique to compute such approximations efficiently by discretizing the space of transformations, which can then be explored in a systematic manner. Implicit in their work is a general reduction technique from the β-approximate formulation of pattern matching under the Hausdorff measure to pattern matching under the exact measure, denoted by dE (P, Q) = 0 if P ⊆ Q, and 1 otherwise. Theorem 2.1 ([26]). Let there be an algorithm for PM(dE , T ) running in time O(f(n). Then there exists an algorithm for PM(dH , T , ², β) running in time O(f(n)/βO(d) ), where d is the dimension of the space of transformations T. In the sequel, we will address the approximation versions of CONG, PM and LCP under the Hausdorff measure and under different groups of transformations. Following the notation of [25], let Td denote the space of translations in Rd . Let SOd denote the group of orthogonal transforms, and Md = Td · SOd denote the group of all rigid transformations. 3
3
Related Work
Early work in this area was done by Huttenlocher, Kedem and Sharir in 1993 [28], where they presented an algorithm running in time O(kn(k + n)α(kn) log(k + n)) for PM(dH , T2 , ²) This was generalized to rigid motions ˜ 3 n2 ). Briefly, the algorithms proceed by computing a partition in [9], yielding an algorithm running in time O(k of transformation space such that each cell of the partition corresponds to a fixed cost assignment of points in P to points in Q. The best such transformation is then found by traversing this partition and evaluating the cost of each cell. For a more detailed overview of the techniques, the reader is referred to the survey by Alt and Guibas [3]. 2 ˜ In the approximate setting, Goodrich, Mitchell and Orletsky [24] developed a very simple O(kn ) algorithm for 2 PM(dH , M , ², 2 + γ) (for any γ ≥ 0). Their scheme is based on the alignment method, a popular heuristic for point set matching that we outline in Algorithm Naive Alignment below. Prior to our work, this was the best known algorithm for approximate pattern matching in the plane. Distance Distributions The study of combinatorial distance bounds was initiated by Paul Erd˝os [18] in 1946, when he posed the question “Given a set of n points in the plane, how many pairs of points can be a unit distance apart?”. Erd˝os showed an upper bound of O(n3=2 ) and a lower bound of Ω(n1+c= log log n ). Szemeredi and Trotter [38] reduced the upper bound to O(n4=3 ). Later results [10, 37] reduced the constants and simplified the proof considerably. The lower bound remains unimproved and is conjectured to be tight. Akutsu, Tamaki, and Tokuyama [2] and Boxer [6], made the key observation that these bounds could be used to improve the analysis of simple point set matching algorithms. Consider the following naive algorithm to solve PM(dE , M2 ) i.e to determine whether a k-point set P ⊂ R2 is contained in an n-point set Q ⊂ R2 . Firstly, fix a pair p1 , p2 ∈ P. Then, for each pair q1 , q2 ∈ Q, determine the (unique up to reflection) transformation T mapping p1 → q1 , p2 → q2 , if one exists. For each T check if every point in T (P) matches a point in Q; if so, return T . We refer to this algorithm as Naive Alignment. Algorithm Naive Alignment: 1. Fix a pair p1 , p2 ∈ P, 2. For each pair q1 , q2 ∈ Q, (a) Determine the (unique up to reflection) transformation T mapping p1 → q1 , p2 → q2 , if one exists. (b) Compute T (P). If every point in T (P) matches a point in Q, return T . This algorithm has running time O(kn2 log n), where a factor of k log n comes from the cost of transforming P and determining the nearest neighbor in Q for each point in T (P) [13]. However, given a pair (p1 , p2 ), we need only consider pairs (q1 , q2 ) such that |p1 − p2 | = |q1 − q2 |. By the combinatorial distance bounds, we know that there are O(n4=3 ) such pairs, and Agarwal, Aronov, Sharir, and Suri [1] have given an effective way of extracting these distances in time O(n4=3 log n). This immediately improves the running time of the naive algorithm to O(kn4=3 log n). Using more sophisticated ideas with a similar flavor, Akutsu et al. presented an algorithm for PM(dE , M2 ) running in time O(k1:43 n0:77 ). This technique can be generalized to higher dimensions, where the combinatorial object of interest is now a simplex in Rd . Akutsu et. al provide results for PM in three and higher dimensions using such bounds. More recently, Braβ [7] has used similar methods for exact point set matching in three dimensions. Interestingly, a simple construction presented by Erd˝os in 1967 [19] shows that for any d ≥ 4, there exists a set of n points in Rd having Θ(n2 ) pairs of points a unit distance apart. Approximate Distance Bounds A natural question would be whether we could extend the above methods to deal with measures like dH . For this purpose, we need a noisy version of the above combinatorial problem. Such a 4
question was posed and solved (in the plane) by Erd˝os, Makai, Pach and Spencer in 1991: Theorem 3.1 ([20]). Given t > 0, and a set of n points in the plane ¥with minimum distance at least 1, the number ¦ 2 of pairs (i, j) such that d(pi , pj ) lies in the range [t, t + 1] is at most n /4 , provided n is sufficiently large. They also show that this bound is tight for sufficiently large t. However, this requires t = Ω(n2 ). Notice that the minimum interpoint requirement cannot be removed; if this condition were not present, we could ensure that all n(n − 1)/2 distances were in the range by clustering all the points together. In itself, this bound provides no asymptotic improvement on the worst case. However, the bounds are tight only for sets with large diameter. In the sequel, we extend and generalize the above bound by explicitly establishing the dependence of the bound on the diameter of the point set. The connection between the diameter of a point set and such approximate combinatorial quantities has been exploited before. Valtr and others [14, 39, 40, 41] investigated bounds on the number of approximate incidences between families of well-separated lines and points; assuming that the point sets are dense, i.e., their diameter √ is O( n), he showed that several key results developed for the exact incidence model have counterparts in the approximate setting.
4
A Bound On Approximate Distances
In the sequel we use B(p, r) to denote a 2 or 3-dimensional ball of radius r centered at p. We also use A(p, r1 , r2 ) to denote a 2 or 3-dimensional annulus B(p, r2 ) − B(p, r1 ). For a point p, let Nt (p) = P ∩ A(p, t, t + ²). We denote a set of points as δ-separated if the minimum inter-point distance in this set is δ. Fact 4.1. For any n, there exists √ a set of n 1-separated points in the plane such that the number of distance pairs in the interval [t, t + 1] is Ω(n t). Proof. Consider two rows of n/2 points that are t units apart. In each row, the points are spaced at unit intervals. If we consider the annulus of inner radius t and outer √ radius t + 1 swept out by any point√in the first row, it will intersect a line segment on the second row of length t, which implies that it intersects Ω( t) points. Lemma 4.1. Given ², l > 1 and two points p and q in the plane such that d(p, q) = d and 1 ≤ d < 2l, the number ofp 1-separated points whose distance from p and q lies in the interval [l, l + ²] is O(²2 l/d) when d ≤ l, and O(²2 l/(2l − d)) when l < d ≤ 2l. ^ = d/². The minimum inter-point separation is now Proof. Let us scale all distances by ², setting ^l = l/² and d 1/². Construct around each point concentric circles of radius ^l and ^l + 1. We need to estimate the diameter of the common area. To this end, we observe that the four circles intersect in exactly eight points (in two symmetric parts). Consider the top four intersection points A, B, C, and D (see Figure 1). We wish to estimate the number of points lying in the region ABCD. It is easy to see that this can be bounded by ²2 max(d(A, C), (B, D)). If we set the origin ^ 0), the coordinates of A, B, C, D are as follows: at p, and let q lie at (d, q ^ − (2^l + 1)/2d, ^ ^l2 − (d/2 ^ − (2^l + 1)/2d) ^ 2 ), A : (d/2 q ^ ^ 2 ), B: (d/2, (^l + 1)2 − (d/2) q ^ + (2^l + 1)/2d, ^ ^l2 − (d/2 ^ − (2^l + 1)/2d) ^ 2 ), C : (d/2 q ^ ^l2 − (d/2) ^ 2 ), D: (d/2, q q ^ 2 − ^l2 − (d/2) ^ 2 ). It is easy to verify (^l + 1)2 − (d/2) p that for d ≤ l, this gives a bound of O(²2 l/d); for l < d ≤ 2l, the bound is O(²2 l/(2l − d)).
^ which yields max(d(A, C), d(B, D))=max((2^l + 1)/d,
5
B A
C D
p.
q. d
Figure 1: The four points of intersection of the two annuli Armed with this combinatorial bound, we can now prove our main lemma, which bounds the number of distance intervals induced by a single point. Lemma 4.2. Let P be a 1-separated set of points in the plane. Given ² > 1 and t, there exists a point p ∈ P such that Nt (p) = O((²2 nt)1=3 log n). Proof. Let f be a parameter such that for all points p ∈ P, Nt (p) ≥ f. We show that f = O(²(nt)1=3 log n). Take any point p ∈ P. As |Nt (p)| ≥ f, there must be a segment Sp of A(p, t, t + ²) subtending an angle π/3 such that |Sp ∩ P| = u ≥ f/6. Let q1 , q2 , . . . , qu be the members of Sp , ordered in clockwise fashion around p. Notice that the distance between any pair qi and qj is at most t. Partition this sequence into subsequences S1 , S2 , . . . , Sk of size l. Let q10 , q20 , . . . , qk0 be the subsequence formed by taking the first member of each Si , for 1 ≤ i ≤ k, and let Ni = Nt (qi0 ). Note that lk ≥ f
(4.1)
By Fact 4.1, two points at distance l ≤ t can have at most O(²2 t/l) common neighbors in the distance range [t, t + ²]. We can approximate arc lengths of a circle by a straight line with only a constant factor error, which implies that for some c and for 1 ≤ j < i ≤ k, |Nj ∩ Ni | ≤ c²2 t/(l(i − j)). Now, i-1 X
|Nj ∩ Ni | ≤
j=1
i-1 X j=1
c
²2 t ²2 t log n ²2 t =c Hi ≤ c (i − j)l l l
where Hi is the ith harmonic number. The final inequality follows from the fact that i ≤ n. 3 The total contribution of ∪j 1, t > 1, and a set of n 1-separated points in the plane, the number of pairs (i, j) such that d(pi , pj ) lies in the range [t, t + ²] is O((n4 t²2 )1=3 ). Proof. Let p be the point having the property described by Lemma 4.2. We can delete p from P, thus removing O(²(nt)1=3 log n) approximate distance pairs. Repeating for each point, the result follows. Notice that for any reasonably sparse point set (∆ = o(n2 )), this yields a bound that is asymptotically better than the worst-case of Θ(n2 ).
5
An Algorithm For Pattern Matching in 2D
We now have the main component of an algorithm for PM(dH , M2 , ², β). We can choose a pair of points (p1 , p2 ) from P, determine all pairs of points (q1 , q2 ) in Q having a separation approximately equal to that of (p1 , p2 ), and then map P to Q using the transformation defined by these pairs. Theorem 4.1 bounds the number of such pairs that we need to examine. We thus need to provide a scheme for extracting these pairs efficiently, and a way of discretizing transformation space to achieve a (², β)-approximation. Let p1 , p2 , q1 , q2 be any points. Then for any ² > 0, we define T (p1 , p2 , q1 , q2 ) = {T |dH ({T (p1 ), T (p2 )}, {q1 , q2 }) ≤ ²} S (p1 , p2 , q1 , q2 ) = {(T (p1 ), T (p2 ))|T ∈ T (p1 , p2 , q1 , q2 )} T is an infinite subset of the transformation space and S is the space of images of P under T . The following proposition (in the spirit of [26]) shows that this space can be approximated by a small set of transformations. Proposition 5.1. For any S = S (p1 , p2 , q1 , q2 ) and γ > 0, there exists a set S∗ ⊂ S of cardinality O(1/γ3 ) such that for any pair (s1 , s2 ) ∈ S, there exists a pair (s∗1 , s∗2 ) ∈ S∗ with dH ({s1 , s2 }, {s∗1 , s∗2 }) ≤ γ². Such a set S∗ is denoted by S∗ ; (p1 , p2 , q1 , q2 ). Let d = dH , and let γ = β/4. The algorithm is as follows: Algorithm 2D Alignment: 1. Compute a diameter pair p1 , p2 of P. 2. Find the set A of all pairs (q1 , q2 ) of points in Q such that kp1 − p2 k − 2² ≤ kq1 − q2 k ≤ kp1 − p2 k + 2². 3. Construct T = ∪(q1 ;q2 )∈A {T0 (p1 , p2 , s1 , s2 )|(s1 , s2 ) ∈ S∗ ; (p1 , p2 , q1 , q2 )}. 4. Search for T ∗ ∈ T such that dH (T ∗ (P), Q) ≤ (1 + 3γ)². Lemma 5.1. Algorithm 2D Alignment is a β-approximation algorithm for PM. Proof. Assume that there exists a transformation T such that dH (T (P), Q) ≤ ². Let p1 , p2 be the diameter points found by the algorithm. Let s1 = T (p1 ) and s2 = T (p2 ) and let q1 , q2 ∈ Q be points such that d(si , qi ) ≤ ² for i = 1, 2. By the definition of S∗ = S∗ ; (p1 , p2 , q1 , q2 ) we know that there exists (s∗1 , s∗2 ) ∈ S∗ such that dH ({s∗1 , s∗2 }, {q1 , q2 }) ≤ γ². Let T ∗ be the transformation which maps (p1 , p2 ) to (s∗1 , s∗2 ). Consider any point p ∈ P. We need to estimate d(T (p), T ∗ (p)). We represent T ∗ as R ◦ I ◦ T , where I is a translation which moves s1 to s∗1 , and R is a rotation centered at s∗1 which moves I(s2 ) to s∗2 . Then d(T (p), T ∗ (p)) ≤ d(T (p), (I ◦ T )(p)) + d((I ◦ T )(p), (R ◦ I ◦ T )(p)) ≤ d(s1 , I(s1 )) + d(I(s2 ), s∗2 )
≤ d(s1 , s∗1 ) + (d(I(s2 ), s2 ) + d(s2 , s∗2 )) ≤ 3γ² ≤ β². 7
Therefore d(T ∗ (p), Q) ≤ d(T (p), Q)+d(T (p), T ∗ (p)) ≤ ²+β² = (1+β)², thus T ∗ gives the desired approximate mapping. It is straightforward to check that our algorithm finds such T ∗ . Note that it is crucial that the pair p1 , p2 found is a diameter pair. This ensures that the second inequality in the equation is satisfied. Lemma 5.2. Algorithm 2D Alignment can be implemented to run in time à õ ¶ !! ² 2 min(²2=3 n4=3 ∆1=3 log n, n2 ) O k · min , log n β3 β for any ² ≥ 1 and β ≤ 1, where ∆ is the diameter of P. For ² < 1 the running time is as for ² = 1; similarly for β > 1. Proof. Consider the complexity of each step of the algorithm: Step 1: Can be performed in time O(n log n) using a standard algorithm [34]. Step 2: Let r = d(p1 , p2 ). We need to find A = ∪q1 ∈Q {(q1 , q2 )|q2 ∈ Q ∩ A(q1 , r − 2², 4²)}. If ∆ ≤ n2 then we employ the following bucketing scheme. Let G denote a planar grid consisting of rectangular cells of width √ √ √ . Notice, that 4² and height r. For each i = 0, . . . , r we define Gi to be the grid G rotated by an angle 2i r √ for any point q the ring A(q, r − 2², 4²) can be covered by O( r) grid cells from G0 , . . . , Gr (this can be proved similarly to Fact 4.1). During the preprocessing for all G0 , . . . , G√r and all q ∈ Q we store q at the grid cell of all Gi containing √ q (using hashing). Then, in order to compute Q ∩ A(q, r − 2², 4²) we retrieve the contents of the O( r) buckets covering A(q, r − 2², 4²). The total complexity of this step can be bounded by the total number of 4=3 1=3 2 the points retrieved (which √ is O(min(²n ∆ log n, n )) by Theorem 4.1) plus the√total number of buckets accessed (which is O(n ∆)) plus the cost of the preprocessing (which is again O(n ∆)). If ∆ > n2 , then we can compute A by computing all n2 pairwise distances in time O(n2 ). Therefore, Step 2 √ can be implemented to run in O(min(n2 , ²2=3 n4=3 ∆1=3 log n + n ∆)) time, which is O(min(n2 , ²2=3 n4=3 ∆1=3 log n)). Steps 3/4: In these steps, we perform at most O(min(²n4=3 ∆1=3 log n, n2 )/β3 ) computations of dH (T (P), Q). Each computation requires k comparisons of dH (T (p), Q) to (1 + 3γ)². It is easy to perform this comparison in O(log n) using a Voronoi diagram. One can however achieve constant time per query as follows. Impose a uniform grid of side γ². For each q ∈ Q we store it in all grid cells intersecting B(q, (1 + 3γ)²). This takes O(n/γ2 ) time which is subsumed by the complexity of searching T . In order to verify if dH (p, Q) ≤ (1+β)², it is sufficient if the grid containing p contains any q.
6
Bounds For Noisy Triangles
Having shown combinatorial bounds for approximate distances in the plane, a natural question is whether we can show an analogous result in three dimensions, where the corresponding combinatorial object is a triangle. Let ∼ d 4 0 if and only if dH (4, 4 0 ) ≤ ². For a 4 = (p1 , p2 , p3 ) and 4 0 = (p10 , p20 , p30 ) be two triangles. Then 4 = point set P where |P| = n, we define HP (4), the multiplicity of 4, as the following; ∼ 4, pi ∈ P, i = 1, 2, 3}| HP (4) = |{4 0 = (p1 , p2 , p3 )|4 0 = Let Hn (4) = maxjPj=n HP (4). In the sequel, we drop n from Hn if its value is clear from the context. √ Theorem 6.1. For any triangle 4 of largest side l, H (4) = O(min(²2:5 n2:25 l, n3 )).
8
Proof. Take any 1-separated point set P. Let 4 = (p1 , p2 , p3 ) and assume kp1 − p2 k = l. Let a = a(4) be the point where the line p1 − p2 intersects the altitude from p3 to the base p1 − p2 ; define r = ka − p3 k. For convenience, we assume that α = kp1l-ak ≥ kp2l-ak ; it can be verified that the other case can be treated similarly. For 0 ≤ β ≤ 1, we define J (p, p 0 ) = (1 − β)p + βp 0 . ∼ (p 0 , p 0 , p 0 ) where p1 corresponds to p 0 and Let C4 (p10 , p20 ) denote the set of all points p 0 in R3 such that 4 = 1 2 1 0 p2 corresponds to p2 . For any point p, vector v and reals ρ and δ we also define R(p, v, ρ, δ) to be the torus centered at p orthogonal to the vector v, with outer radius ρ and cross section δ. The following fact allows us to approximate C4 (p10 , p20 ) by a torus. Fact 6.1. There exists constant c such that C4 (p10 , p20 ) ⊂ R(J1= (p10 , p20 ), p20 − p10 , r, c²). ∼ (p, s, p 0 ), For any point p define Bp = A(p, l − 2², l + 2²) ∩ P. Note that Bp contains all points s such that 4 = for some p 0 ∈ P. Each s ∈ Bp defines a torus Rps = R(J1= (p, s), p − s, r, c²). Let Rp = {Rps | s ∈ Bp }. The following inequality allows us to obtain a bound for H (4): H (4) ≤
X X
|R ∩ P| ≤ n · max p∈P
p∈P R∈Rp
X
|R ∩ P|.
R∈Rp
P Therefore, it suffices to determine an upper bound on the expression I(Rp , P) = R∈Rp |R ∩ P| for any fixed p ∈ P. We recall that the points in Bp lie between two spheres Sp (resp. Sp0 ) of radius l − 2² (resp. l + 2²). Therefore the centers of all R ∈ Rp (which we we denote by ρ(R)) are contained in a ring A = R(p, (l − 2²)/α, (l + 2²)/α). It is also easy to see that the centers of all R ∈ Rp are 1/α-separated, as all points in Bp are 1-separated. We split 1 . In the following, each A into disjoint concentric rings A1 , . . . , Au for u = O(²), such that each ring is of width 2 we estimate the value of I(Rp ∩ Aj ) for a fixed j; when multiplied by ², this gives the desired bound. Clearly, if the centers ρ(R), for R ∈ Rp ∩ Aj , are projected along a line joining them to p onto B(p, l/α), they are still c 0 /α separated for some c 0 > 0. Notice that by perturbing them in this way we only change the bound by a constant. Therefore, we assume that all centers of tori from Rp lie on the sphere Cp of B(p, l/α) and are 1/α-separated. Let N(R) be the number of points of P inside torus R. Fix a number f ≤ n. We proceed as in the 2-dimensional case. Stage 1: Delete all tori R for which N(R) ≤ f. Stage 2: Now all tori have at least f points in them. In the following, we upper bound f. 0 0 PDefine φ(R, R0 ) to be the number of points in P that an torus R can share with another torus R . Let Φ(R) = R 0 ∈Rp φ(R, R ). We now bound Φ(R). The proof considers two cases: when r > cl/α ( for some c > 0) and when r < cl/α. Case 1 [r > cl/α, for c > 0]. First, we partition the sphere B(p, l/α) into a constant number of cones, such that each cone spans an angle less than π/3. Clearly, it is sufficient to bound the above sum only for tori R 0 in a fixed cone (and multiply the result by a constant). Therefore, we assume that all tori centers lie in a fixed cone C. In this case, we have the following lemma whose proof is similar to that of Lemma 4.1. , for some c > 0. Lemma 6.1. Let R 0 be a torus such that |ρ(R) − ρ(R 0 )| = u. Then φ(R, R 0 ) ≤ c²3 l= u Now we proceed as follows: Φ(R) =
X R 0 ∈Rp
φ(R, R 0 ) ≤
²3 l X φ(ρ(R), ρ(R 0 )). α 0 R ∈Rp
9
As the centers ρ(R 0 ) are 1/α-separated and lie on a sphere, we know that the number of R 0 ’s such that |ρ(R)−ρ(R 0 )| ∈ [i/α, (i + 1)/α] is O(i). Therefore, the above sum can be estimated by √
n
√ ²3 l X 1 i = ²3 l n. α i/α i=1
Case 2 [r ≤ cl/α]. As before, we scale so that the inter-sphere and torus width is 1. Consider one hemisphere of Cp (we lose only a constant factor by this). Project this hemisphere onto the plane (again losing only a constant factor). By Lemma 4.1 in Section 4, the fact that no torus has more than f points contained in it, and applying a volume scaling of ², ± min(²3 q r/d(R, R 0 ), f) d(R, R 0 ) ≤ r 0 φ(R, R ) = r min(²3 r-d(R;R 2r > d(R, R 0 ) > r 0 )=2 , f) Since φ(R, R 0 ) is unimodal with a minimum at d(R, R 0 ) = r, we maximize Φ(R) by locating R 0 maximally (with respect to φ) around R. 0 Let R1 = {R 0 ∈ r} and R2 = {RP0 ∈ Rp | r ≤ d(R, R 0 ) < 2r}, and then Φ(R) = Φ1 (R) + Φ2 (R) PRp | d(R, R ) ≤ where Φ1 (R) = R 0 ∈R1 φ(R, R 0 ) and Φ2 (R) = R 0 ∈R2 φ(R, R 0 ). Note that if d(R, R 0 ) > 2r, then R and R 0 do not intersect. Let η = maxR 0 ∈R1 d(R, R 0 ); similarly, let λ = maxR 0 ∈R2 2r − d(R, R 0 ). √ Lemma 6.2. η = O(min(r, n/α)) and λ = O(η2 /r). Proof. Clearly, η ≤ r. All tori centers are 1/α−separated, hence the number of tori R ∈ R1 is O(η2 α2 ). But √ |R1 | ≤ n, implying that η ≤ n/α. Similarly, λ ≤ r. The area of the torus of inner radius 2r − λ and outer radius 2r is O(rλ). As before, since tori centers are 1/α−separated, this implies that |R2 | = O(rλα2 ) ≤ n which yields λ = O(n/rα2 ). Consider any interval I = [i, i + 1]. All the tori centers are separated by 1/α, therefore |{R 0 | d(R, R 0 ) ∈ I}| is O(iα2 ). Using this, and substituting for φ(R, R 0 ) we can write: X
Φ1 (R) ≤
iα2 min(²3 r/i, f),
i=1
Φ2 (R) ≤
r
2r -1 X
2
iα min(f, ²
3
i=2r-
¯ Lemma 6.3. Φ(R) =
r ). r − i/2
O(α2 fη2 ) r > ηf/²3 O(α2 rη²3 ) r ≤ ηf/²3
Proof. Case 1 [r ≤ f/²3 ]. Φ1 (R) ≤
X
α2 r²3 = α2 rη²3
i=1
Φ2 (R) ≤
r
2r -1 X
2 3
iα ²
i=2r-
10
r r − i/2
Let j = r − i/2. =2 √ 2 3 X r−j √ Φ2 (R) ≤ 2 rα ² j j=1=2
=2 X
≤ 2r3=2 α2 ²3
p 1/ j
j=1=2
= O(r
3=2
2 3
α ²
Z =2
dj √ ) j
1=2
√ = O(r α ² λ) √ = O(α2 ²3 r rλ) 3=2
2 3
= O(α2 rη²3 ) from Lemma 6.2 Case 2 [ηf/²3 ≥ r > f/²3 ]. Let γ = ²3 r/f. From above, we know that 1 < γ ≤ η.
X
Φ1 (R) ≤
2
α if +
i=1 2
X
α2 ²3 r
i= 2
≤ α [γ f/2 + ²3 (ηr − γr)] ≤ α2 ²3 (ηr − γr/2) since γf = ²3 r. = O(α2 ηr²3 ) as γ ≤ η Applying the change of variable j = r − i/2 as above, we obtain Φ2 (R) ≤ 2α2
=2 X
r (r − j) min(f, ²3
j= 1
r
=2
2
≤ α r
X
3
min(f, ²
j=1
r ) j
r ) j
Let γ = ²6 r/f2 . We get Φ2 (R) ≤ α2 r[γf + ²3
=2 r X r j=
j
]
Z =2 r
r dj]) j
√ √ = O(α2 r[γf + ²3 ( λr − γr)]) √ = O(α2 r(γf + 2²3 (η − γr))) from Lemma 6.2 √ = O(α2 rη²3 ) as γf = ²3 γr 2
3
= O(α r[γf + ²
3 Case 3 [ηf/² ≤ r]. Since ηf/²3 ≤ r, min(f, ²3 r/i) = f ∀i ≤ η. Hence, Φ1 (R) = O(α2 fη2 ). In addition, q r min(f, ²3 j ) = f ∀j ≤ η2 /r = Ω(λ). Therefore, Φ2 (R) = O(α2 fη2 ).
11
Putting it all together: Φ(R) is the total number of points that torus R can share with all other tori. Assume that k such tori suffice to ensure that all tori contain f points each. This implies that k maxR∈R Φ(R) ≥f n by the Pigeon-hole Principle. In addition, we know that kf ≤ n, since each torus has f points in it. These two inequalities give us an upper bound on f. √ √ Case 1 p [r > cl/α]. Φ(R) = O(l n²3 ). Therefore, we have kl nn²3 /n ≥ f and kf ≤ n. The two inequalities √ yield f ≤ l²3 n. This yields I(Rp , P) = nf = O(n1:25 l0:5 ²3 , which gives H (4) = O(n2:25 l0:5 ²4 ) Case 2 [r ≤ cl/α]. There are two sub-cases to consider: p Sub-case 1 [ηf/²3 ≥ r]. f ≤ k²3 α2 rη/n, and f ≤ n/k. Combining these yields f ≤ α rη²3 . Sub-case 2 [ηf/²3 < r]. As before, the calculation yields f ≤ α2 η2 . This yields p ηf/²3 ≥ r ⇒ f ≤ α rη²3 ηf/²3 < r ⇒ f ≤ α2 η2 √ where η = min(r, n/α) from Lemma 6.2. √ (a). r ≤ n/α, η = r √ In this case, the two conditions collapse to f ≤ αr²3=2 ≤ n²3=2 . √ √ (b) r > n/α, η = n/α √ Let y = αr²3 . Substituting η = n/α in the above conditions, and replacing αr²3 by y, we obtain q √ √ f ≥ y/ n ⇒ f ≤ y n, √ f < y/ n ⇒ f ≤ n. p √ p √ √ √ For y ≤ n3=2 , we have y/ n ≤ y n ≤ p n. For y ≥ n3=2 we have n ≤ y n ≤ y/ n. Therefore, √ from the above we obtain f = min( y n, n). Since r ≤ cl/α, y = O(l²3 ) which implies that p conditions √ f = O(min( l²3 n, n)). Thus, we obtain H (4) = O(min(max(n2:25 l0:5 , n2:5 )²5=2 , n3 )). Fact 6.2. For any l > 1 such√that l = O(n2 ), there exists a point set P and a triangle 4 with the largest side length l such that H1P (4) = Ω(n2 l).
7
A 3D Pattern Matching Algorithm
The above combinatorial result can be translated into an algorithm for PM(dH , M3 , ², β). As in Section 5, we first state a proposition pertaining to the discretization of the transformation space. Let p1 , p2 , p3 and q1 , q2 , q3 be any points. Then for any ² > 0, we define T (p1 , p2 , p3 , q1 , q2 , q3 ) = {T |dH ({T (p1 ), T (p2 ), T (p3 )}, {q1 , q2 , q3 }) ≤ ²} S (p1 , p2 , p3 , q1 , q2 , q3 ) = {(T (p1 ), T (p2 ), T (p3 ))|T ∈ T (p1 , p2 , p3 , q1 , q2 , q3 )} Proposition 7.1. For any S = S (p1 , p2 , p3 , q1 , q2 , q3 ) and γ > 0, there exists a set S∗ ⊂ S of cardinality O(1/γ6 ) such that for any pair (s1 , s2 , s3 ) ∈ S, there exists a triplet (s∗1 , s∗2 , s∗3 ) ∈ S∗ with dH ({s1 , s2 , s3 }, {s∗1 , s∗2 , s∗3 }) ≤ γ². Such a set S∗ is denoted by S∗ ; (p1 , p2 , p3 , q1 , q2 , q3 ). 12
The algorithm itself is quite simple. Following the alignment paradigm outlined in Section 5, we perform the following four steps: Algorithm 3D Alignment: 1. Compute a diameter pair p1 , p2 of P. Choose a third point p3 in P arbitrarily. Let 4 denote the triangle (p1 , p2 , p3 ). 2. Find the set A of all triplets (q1 , q2 , q3 ) of points in Q such that dH (4, (q1 , q2 , q3 )) ≤ 2². 3. Construct the set T = ∪(q1 ;q2 ;q3 )∈A {T0 (p1 , p2 , p3 , s1 , s2 , s3 )|(s1 , s2 , s3 ) ∈ S∗ ; (p1 , p2 , p3 , q1 , q2 , q3 )} 4. Search for T ∗ ∈ T such that dH (T ∗ (P), Q) ≤ (1 + 3γ)². The following lemma is immediate: Lemma 7.1. Algorithm 3D Alignment is a β-approximation algorithm for PM. Let l, r, s denote the sides of the triangle 4 = (p1 , p2 , p3 ) selected in Step 1 of the algorithm, where l ≥ r ≥ s is the length of the longest side. Let (q1 , q2 ) be a pair of points in Q such that dH ({p1 , p2 }, {q1 , q2 }) ≤ 2². Following the notation of Section 6, let C4 (q1 , q2 ) denote the set of all points q3 such that (q1 , q2 , q3 ) ∈ A. From Fact 6.1, we know that C4 (q1 , q2 ) can be approximated by a torus of outer radius l 0 ≤ l and width O(²). Let this torus be denoted by R4 (q1 , q2 ). Consider the ball B(q, r) centered at q. For any q 0 ∈ Q − {q}, R4 (q, q 0 ) is contained within this ball. Let ~v denote the unit vector directed from q to the center of the torus R4 (q, q 0 ). Consider the projection of R4 (q, q 0 ) onto the surface of B(q, r), where the image of p ∈ R4 (q, q 0 ) is the point of intersection of B(q, r) with a ray originating from p in direction ~v. Note that the projections change all dimensions by at most a constant factor. The image of R4 (q, q 0 ) under this projection is an annulus of width O(²) and radius O(l). Let B(q, r) be decomposed into a constant number of uniformly-sized caps, and consider the projection of such a cap into R2 . Each such cap is partitioned by √ (elliptical) sections of the annuli, each of which (by arguments similar to those in Section 5) can be covered by O( l) rectangles of width O(²). √ Thus given the collection of tori Rq = ∪q 0 ∈Q-fqg R4 (q, q 0 ), we can construct O( l) grids such that each torus in this √ √ collection is covered by at most O(1) cells from each grid (where each cell is a rectangular pipe of length O( l) having as cross-section a square of side O(²)). Repeating this for each point q ∈ Q, we obtain O(n l) such grids. In a preprocessing step, for each pair of points q1 , q2 ∈ Q, the grid cells covering the torus R4 (q1 , q2 )) are √ 2 labelled with the pair (q√ l) overall. In the next step, for every point q ∈ Q, the grid 1 , q2 ). This step takes O(n cell in each of the O(n l) grids corresponding to q is located. Each label (q1 , q2 ) present in one of these cells ∼ 4. The total time taken to extract all such triangles is bounded by the number of yields a triangle (q1 , q2 , q) = √ cells and the number of triangles obtained. √ Summing over all q ∈ Q, the first number is O(n2 l) and the second number (from Theorem 6.1) is O(min(n2:25 l, n3 )). For each such triangle, an alignment can be performed in time ˜ O(k) using methods similar to those in Section 5. Therefore, we have (omitting terms in ², β): Theorem 7.1. Algorithm 3D Alignment can be implemented to run in time à õ ¶ !! √ min(²2:5 n2:25 ∆, n3 ) ² 3 O( k · min , log n β6 β where ∆ is the diameter of P.
13
8
Other Alignment-Based Methods
The techniques above imply algorithms for the corresponding CONG formulations of the matching problem; however they do not yield any results for the LCP formulation. This is because the alignment-based methods rely on the pair of points from P being a diameter pair. Even if we could choose a pair of points from P that is guaranteed (possibly with high probability) to be in the common overlap, we cannot guarantee that this pair is a diameter pair (or even has a distance within a constant factor of the diameter). In order to solve LCP we employ a voting scheme similar to the Hough transform [5]. First, we discretize the rotation space, obtaining O(∆) different angles. Then, for each angle we enumerate all pairs p ∈ P and q ∈ Q. Each such pair “votes” for a translation moving p to q. At the end we count the number of votes for each (rotation,translation) pair and output all pairs with more than k votes. This yields an algorithm for LCP(dH , M2 , ², β) running in time1 O(∆kn). In three dimensions, the bounds are obtained in a similar way, yielding an algorithm running in time O(∆3 kn) In Table 1, we summarize our approximation algorithms for the problems described above. For clarity, we omit terms that are polynomial in 1/² and 1/β. Our algorithms are enumerative, in that they can also enumerate all transformations satisfying the approximation criteria. It should be noted that all the algorithms are deterministic.
LCP PM
2D ˜ O(∆kn) 4 ˜ O(k(n ∆)1=3 )
3D ˜ 3 kn) O(∆ √ ˜ min(n2:25 ∆, n3 )) O(k
Table 1: Results for two and three dimensions Sum Measures. A variant of the Hausdorff distance, called the Σ−Hausdorff distance, can be defined as follows: Definition 8.1 (Σ-Hausdorff distance). Given point sets P and Q, the Σ-Hausdorff distance from P to Q is X d min d(p, q) H (P, Q) = p∈P
q∈Q
Consider any transformation T . Suppose that for any ² we know the number f(²) of points from T (P) which are within distance ² of Q. Assume for simplicity that all distances from p ∈ P to Q are distinct. Define f-1 (i) to be the “inverse” of f, i.e., the function which when given a number i ∈ {0, . . . , n} returns the smallest ² such that f(²) = i. P Then it is easy to see that the Σ-Hausdorff distance between T (P) and Q is equal to ki=1 f-1 (i). Similarly, we ˜ can compute an approximation of that distance if we know an approximation of f, i.e., a function f˜ such that f(²) ∈ {f(²), . . . , f((1 + β)²) }, for some β > 0. Such an approximation can be obtained by solving a logarithmic number of instances of approximate LCP between T (P) and Q under Hausdorff measure for ² = 1, (1 + γ), (1 + γ)2 , . . . , ∆ and γ = Θ(²). Thus, we can obtain an approximation of the Σ-Hausdorff for a fixed transformation T . However, our Hough-transform-based algorithm computes the LCP for the whole (properly discretized) space of transformations. Therefore we can choose the one which minimizes the distance value, obtaining an approximate solution for PM under the Σ-Hausdorff measure. An algorithm for LCP under the same measure can be similarly obtained.
9
Combinatorial Pattern Matching
The algorithms of the previous sections make use of approximate distance distributions on point sets to achieve their bounds. We now present a completely different approach to approximate pattern matching that exploits techniqes from combinatorial pattern matching. The main result of this section is an algorithm running in time O(n(n + ∆)) for PM(dH , M2 , ², β) that employs a novel reduction to the problem of subset matching: 1
We omit polylogarithmic factors in 1= and 1= .
14
Definition 9.1 (Subset Matching). Given a text string t[0, . . . , n − 1] and a pattern string p[0, . . . , m − 1] over a finite alphabet Σ such that for all i and j, t[i], p[j] ⊆ Σ, compute a binary array o[0, . . . , n − 1] such that o[i] = 1 if and only if p[j] ⊂ t[i + j mod n] for all j ∈ {0, . . . , m − 1}. An example instance of a subset matching problem is shown in Figure 2. In this example, the first three positions of the pattern match the corresponding positions in the text.
Text
x p q r
Pattern
a
c
f
1
g
d
f d
3 t p4 y
3 4
t
q
r
x
q
Figure 2: An instance of subset matching
Recent work of Cole and Hariharan [11] showed that this problem can be solved in (randomized) O(n log3 m) time. Subsequently, Cole, Hariharan, and Indyk [12] have shown that the above algorithm can be derandomized to run in O(n log3 m) time. We will solve PM(dH , M2 , ², β) by reducing it to multiple instances of the subset matching problem. The basic idea of the algorithm is as follows. From earlier remarks we can assume that the rigid transformation to be computed can be expressed as a composition of a translation and a rotation. 1. Choose an arbitrary point p in P and translate it to all points in B(q, ²), for all q ∈ Q. Using the same technique as in Section 5, we can replace this infinite space of translations by a discrete set. 2. For each such alignment, subdivide the plane into concentric rings around p, of increasing radius. Each ring will yield one instance of a subset matching problem. 3. Further subdivide each ring concentrically and radially into cells of fixed size. Since the point set is 1separated, this implies that each cell contains at most one point, for appropriate choice of constants. 4. For any ring, all the cells at a fixed distance from the center have the same identifier. Construct a text string made up of identifiers of cells that contain points of Q, where each position in the string corresponds to the set of all such cells that lie along the same radial line. Similarly construct a pattern string using points of P. Each match of text and pattern corresponds to a rotational shift of the corresponding ring. It is now easy to see that there exists a rotation matching P approximately into Q iff there exists a single rotational shift that corresponds to a pattern match for each ring. We now present the algorithm in more detail. (1) Choose an arbitrary point (say p) from P. For each q ∈ Q align p to a set of O(1/β2 ) points spaced uniformly inside the disk of radius ² centered at q; this is done similarly as in Algorithm 2d Alignment. For each such translation T , move the points P according to T ; for simplicity we still refer to T (P) as P. (2) In the next step, split the plane into l = O(log ∆) concentric rings R1 , . . . , Rl centered at p (see Figure 3); note that R1 is a full disk. The inner radius of the i-th ring (for i ≥ 2) is equal to ²2i-1 , the outer radius (for i ≥ 1) is ri = ²2i . For each ring Ri , let Ii be the (infinite) set of rotations which brings P ∩ Ri within Hausdorff distance ² of Q ∩ (Ri-1 ∪ Ri ∪ Ri+1 ). Clearly, P is within Hausdorff distance ² of Q under rotation iff I1 ∩ . . . ∩ Il 6= ∅; this follows from the fact that a disk of radius ² intersects at most two rings. Therefore it is sufficient to compute the sets I1 , . . . , Il and find their intersection. 15
p
γ
ε2
i
ε2
i+1
Figure 3: The grid decomposition Each Ii can be represented as a finite union of intervals. In order to solve the approximate matching problem, we approximate the Ii as follows: the endpoints of the intervals in each Ii restricted to multiples of αi = γ/2i , for γ > 0 specified later. Denote the approximate intervals by I˜ i . Include an interval [jαi , (j + 1)αi ] in I˜ i if Ii ∩ [jαi , (j + 1)αi ] 6= ∅; we call such an interval important. Notice that rotation of P ∩ Ri by angle αi changes the position of each point in P ∩ Ri by at most 2i αi = γ. Therefore, we can ensure that all important intervals are contained in I˜ i by including all intervals [jαi , (j + 1)αi ] such that rotating Pi = P ∩ Ri by both angles jαi and (j + 1)αi results in a set which is within Hausdorff distance ² + γ from Qi = Q ∩ (Ri-1 ∪ Ri ∪ Ri+1 ). This admits some false matches, but as we explain below, the error they induce is bounded. (3) To find all such intervals, check for each angle jαi if Pi0 ⊂ Si , where the set Pi0 is obtained by rotating Pi by jαi and Si = ∪q 0 ∈Qi B(q 0 , ² + γ). Partition each Ri into 2π/αi sectors, the sectors being partitioned further by 2i /γ uniformly placed concentric circles (Figure 3). We denote the set of grid cells obtained from the ring Ri by Gi ; the union of all Gi ’s (i.e., the whole partition) is denoted by G. For any point x, let G(x) be the cell of G to which x belongs (ties are broken arbitrarily); the function G can be extended to sets of points in a natural way. Each grid cell √ has diameter cγ, for c > 0 with value near 2. (4) Now, for each angle jαi , check if G(Pi0 ) ⊂ G(Si ) for the set Pi0 obtained by rotating Pi by an angle jαi . This condition clearly implies that Pi0 ⊂ Si ; although it also introduces false matches, the error can be bounded as claimed. Notice however, that rotation of Pi by αi results in “shifting” the set G(Pi ) by one position in Gi . Therefore, it is sufficient to check (for all shifts) if the shifted set G(Pi ) belongs to G(Si ). This can be done using the subset matching algorithm as follows. Define a signature of a grid cell to be its distance to the origin point p; note that all grids from one sector of a ring Ri have different signatures, while the signatures of cells from different sectors can be equal. Define the pattern p to be a sequence of sets p[0], p[1], . . . such that each set p[j] contains signatures of grid cells from G(Pi ) belonging to the jth sector of Ri (the first sector is chosen arbitrarily). The text t is constructed analogously to G(Si ). It is easy to check that the subset matching algorithm finds the desired shift. By the above discussion it follows that whenever there is an angle δ such that rotating P by δ results in a set P 0 within distance ² to Q, the intersection I˜ 1 . . . I˜ l will be non-empty. On the other hand, we know that for any q 0 ∈ Q the distance from q 0 to any point in G(B(q 0 , ² + γ)) is at most ² + (1 + c)γ. Therefore, a false match rotation occurs only if it brings P with distance ² + (1 + c)γ to Q; for sufficient γ, this distance can be made less than β². This scheme constructs l = log ∆ instances of the subset matching problem. Each instance is drawn from an alphabet of size O(∆/²β). The pattern (and text) lengths are O(2i /γ) = O(∆/²β). From this, and the above algorithm, we obtain the following theorem: 16
Theorem 9.1. 2D Concentric Pattern Matching is a β-approximation algorithm for the PM problem, for suitable γ = Θ(β) . It can be implemented in O(n( + n2 ) log3 n · log ∆) time, for any β ≤ 1. For β > 1 the running time is as for β = 1.
10
An Experimental Study
In Section 5, we studied the alignment paradigm in some detail, presenting an algorithm for PM(dH , M2 , ², β) 4 ˜ that runs in time O(k(n ∆)1=3 ). In addition, Section 9 also presented a new paradigm that obtains an efficient algorithm for the problem by reducing it to a combinatorial pattern matching problem; the algorithm running time is ˜ O(n(n + ∆)), a significant improvement over the alignment-based method. We now undertake a detailed empirical study of these schemes. In particular, we investigate the performance of the following algorithms for PM(dH , M2 , ², β): BA: An alignment scheme based on the work of Goodrich et al [24]. MGRID: The alignment-based scheme from Section 5. GRID: A simple scheme that modifies MGRID by performing a preliminary filtering. CPM: The algorithm based on combinatorial pattern matching from Section 9. For all schemes we give both analytical and experimental estimations of their performance. Specifically, we analyze the running times of the first three algorithms for randomly distributed input; the (quadratic) running time of CPM does not depend on the point set distribution. We show that when the number of matches is small, BA runs in ˜ 2 ) time, while both MGRID and GRID run in time O(n ˜ 1:5 ). On the other hand, when the number of matches is O(n √ large, the above bounds are multiplied by kn. Additionally, we show that for arbitrary input instances, the running 4=3 1=3 ˜ 6=5 k3=5 ∆3=5 ); this is close to the bound of O(kn ˜ time of GRID is O(n ∆ ) for MGRID given in [31]. These analytical results suggest an interesting tradeoff between alignment-based schemes and CPM parameterized in terms of the number of matches and the pattern/image sizes, at least in the asymptotic sense. However, CPM is more complex and therefore has larger constant factors in the running time bounds, which could potentially affect the running time significantly for medium-size input instances. Therefore, we perform several experimental studies of these schemes to measure their real performance. The conclusions from our study can be roughly summarized as follows: • For typical values of ², GRID works best. This is due to the fact that it is simple (and thus has small big-oh constants) as opposed to CPM and MGRID; at the same time, it filters matches effectively, which allows it to beat BA. • For small ², the running time of CPM is order(s) of magnitude larger than the running time of the other algorithms. On the other hand, this time does not increase with ² (unlike the other schemes); thus, for large values of ² CPM outperforms the other algorithms. • For small ², the running time of GRID does not increase with the size of the pattern, as the worst-case bounds would suggest. What is probably even more counterintuitive is that the running time of GRID actually decreases as pattern size increases. This surprising behavior is a consequence of the fact that as the pattern size becomes larger, its diameter increases as well. As a result, if the pattern and the image have comparable diameters, the filtering mechanism of GRID performs extremely well.
17
11
The Algorithms
We describe implementation details for the four algorithms that we study. The first three are alignment-based algorithms that all conform to the same basic template (recall Algorithm 2D Alignment from Section 5): Algorithm ALIGNMENT: 1. Compute a diameter pair p1 , p2 of P. This is a pair p1 , p2 such that diam(P) = kp1 − p2 k = ∆. 2. Find the set A of all pairs (q1 , q2 ) of points in Q such that ∆ − 2² ≤ kq1 − q2 k ≤ ∆ + 2². 3. Construct a set T of (possibly many) transformations mapping (p1 , p2 ) close to (q1 , q2 ) (the details will be given shortly) 4. Search for T ∗ ∈ T such that dH (T ∗ (P), Q) ≤ (1 + β)². Step 1 of the algorithm can be performed in O(n log n) time [34]. We now discuss in detail the implementation of Step 2-4.
11.1
Step 2: Extracting Candidate Pairs
The three schemes differ in their implementation of Step 2. Basic Alignment (BA): This is essentially the scheme employed by Goodrich et al. ’[24] that yields a constant factor approximation to dH (P, Q). We enumerate all pairs (q, q 0 ), q 6= q 0 , retaining only those pairs that satisfy the condition of Step 2. Multiple Grids (MGRID): This is the scheme outlined in Section 5. For clarity of presentation, we recall the basic idea here. √ √ Let G denote a planar grid consisting of rectangular cells of width 4² and height ∆. For each i = 0, . . . , ∆ we i define Gi to be the grid G rotated by an angle 2√ . Notice, that for any point q the ring R(q, ∆ − 2², ∆ + 2²) can √ be covered by O( ∆) grid cells taken from G0 , . . . , G . During the preprocessing step, we store q at the grid cell of all Gi containing q. Now, to compute the set A satisfying the conditions of Step 2, we must compute Q ∩ R(q, ∆ − 2², √ ∆ + 2²) for each point q in Q. In order to compute Q ∩ R(q, ∆ − 2², ∆ + 2²) we retrieve the contents of the O( ∆) buckets covering R(q, ∆ − 2², ∆ + 2²), eliminating all q 0 such that q 0 ∈ / R(q, ∆ − 2², ∆ + 2²). The total complexity of this step (over all q ∈ Q) can √ be bounded by the total number of the points retrieved√plus the total number of buckets accessed (which is O(n ∆)) plus the cost of preprocessing (which is again O(n ∆)). Single Grid (GRID): In this scheme, we use a single square grid, with cells of size µ². We place all points in Q in this grid, and with every point q ∈ Q we store the set of all grid cells covering R(q, ∆ − 2², ∆ + 2²). As before, to compute the set A in Step 2 of the algorithm, we traverse this set of cells, eliminating points as before. As we 6=5 shall see in Section 12, this step runs in time O( nk2=5 ∆3=5 ) (ignoring polynomial terms in 1/²).
11.2
Step 3: Constructing The Transformation Set
In order to solve the β-approximate version of the problem, we can construct a set of O(1/β3 ) transformations mapping the pair (p1 , p2 ) close to (q1 , q2 ), as described in [26]. For simplicity, we merely compute the single transformation T 0 that maps p1 to q1 and aligns p1~p2 to q1~q2 . Since p1 , p2 is a diameter pair, the error incurred in doing this is no more than 2², implying that the approximation is within a factor of 2 of the minimum value. 18
11.3
Step 4: Verifying The Match
In Step 4 of the algorithm we need to verify if for a given transformation T the point set T (P) is close to Q. This is done by checking for each p ∈ P if there is a point in Q within ² from T (p). One √way of implementing this check is as described in Section 5. We subdivide the plane into a grid of cell width β²/ 2. For each grid cell, we make a list of all points of Q contained within it. Now, for each point p ∈ T (P), we find the cell c containing it, and check whether this cell (and all adjacent cells) are empty or not. As it turns out, a more efficient implementation can be obtained by observing that we only need to test whether a grid cell is empty or not. Therefore, instead of maintaining (for each cell) the set of points contained within it, we need only maintain one bit for each cell which is set if there is some point within distance ² from it. This bitmap can be constructed in advance, by spreading each point q ∈ Q to all cells that intersect B(q, ²). If the cell size is α², α √ ≤ 1, then given a point p ∈ T (P), this procedure will return YES if there is a point of Q within distance ²(1 + 2α) of p. This scheme yields significant improvements in running time over the first method.
11.4
The Combinatorial Pattern Matching Scheme
We use the algorithm of Section 9. The main subroutine is the subset matching algorithm, that involves computing convolutions of binary vectors over (·, +) [30]. We use the NTL package developed by Victor Shoup [36] that implements various operations on finite fields, including fast polynomial multiplication (which we use to implement convolution).
12
Running Time Analysis
We present two analytical estimations of the running time of GRID: a worst case bound in terms of the ratio between the diameter and the closest pair, and an upper bound under probabilistic assumptions about the input.
12.1
Worst-case Analysis
Recall that the algorithm GRID proceeds as follows. First a grid of a specified cell side δ = µ², µ ≥ 1 is imposed on the image. Then the algorithm computes a diameter pair p, p 0 of P; let r = |p − p 0 |. For each point q ∈ Q we select the cells Cq which intersect with an annulus centered at q with inner radius r − 2² and outer radius r + 2². Finally, for all pairs (q, q 0 ), q 0 ∈ Cq we align (p, p 0 ) to (q, q 0 ) and verify the match. Assume that the minimum interpoint distance is 1. Let ∆ ≥ 1 denote the diameter. ˜ 6=5 k3=5 ∆3=5 ). Theorem 12.1. GRID runs in time O(n ˜ 4=3 (∆δ2 )1=3 ). Also, the Proof. By the analysis of Section 4 we know that the total number of pairs found is O(n 2=5 (and assuming this number is greater total number of cells visited is O(n∆/δ), as r ≤ ∆. By setting δ = k3=5 n1=5 6=5 3=5 3=5 ˜ than ²) we get the total running time bound of O(n k ∆ ).
12.2
Probabilistic Analysis
We consider two following random input models. In both cases the pattern is generated by choosing k points uniformly at random from the square [0, t]2 , for t such that t2 = nk ; note that t ≤ 1. This choice guarantees that the density of the points in the image and the pattern is roughly similar. The models are differentiated by the way we select the image. In the first case (which we call model A) the image consists of n points chosen independently at random from the unit square [0, 1]2 . In the second case (called model B) the first n − k points are chosen at random from the unit square. The last k points are chosen by picking a random vector from [1 − t]2 , translating P by this vector and adding the resulting point set into Q. The first model has the advantage of being simple: the pattern and the image are chosen independently, which simplifies the analysis. On the other hand, for small values of ² the probability of an existence of a match is very 19
small; in such a situation it is not very realistic. The model B does not have this drawback (by definition); however, it is more difficult to deal with due to the fact that the pattern and the image are not independent. The estimations which we give below are valid for both models. However, for the sake of simplicity, we give proofs only for model A. We show that if the number of occurrences of the √ pattern in the image is small, then the expected running time of GRID is subquadratic (more specifically, O(n k)). On the other hand if the number of matches is large, then GRID runs in superquadratic time (i.e. roughly O(kn1:5 ). Thus depending on the characteristic of the input data, the algorithm is expected to run either faster or slower than CPM, which runs in quadratic time independently of the input properties. q √ Theorem 12.2. If ² = O( √1n ), then GRID runs in expected time O(n k). For ² > c lognn for large enough c > 1, then the expected running time of GRID is O(kn1:5 ).
Proof. In order to estimate the running time of the procedure, we first observe that the total area Aq of the union of cells from Cq is O(∆δ), which is at most atδ for some constant a > 0. Consider now a specific point q ∈ Q. From the way we generate the random points it is easy to see that the expected number of points in Q falling into the cells from Cq is Aq n ≤ atδn. On the other hand, we know that |Cq | ≤ t . Thus if we denote the expected cost of verifying an alignment by C, then the expected cost CA of the algorithm is the order of CA = n · (t/δ + Ctδn) = tn(1/δ + Cδn). If ² is small (formally O( √1n )), then it is not difficult to verify that C = O(1); note that in this case the number of matches is 0 with high probability. Then the cost of the algorithm is bounded by √ Csmall = O(nt(1/δ + nδ)) = O(n3=2 t) = O(n k) q 1 √ for δ = n . On the other hand when ² is large (formally greater than c lognn for large enough c > 1), then the algorithm cost is Clarge = O(nt(1/δ + knδ)) = O(kn3=2 log n) for δ = ² (since tδ ≤ 1). Notice that in the latter case many matches exist. 2 1:5 Note that a similar analysis to the one above yields an expected running time √ of O(n ) for BA and O(n ) for 1 MGRID when ² ≤ O( √n ). For large ² the above bounds are multiplied by kn.
13
Experiments
We now describe the experiments that we performed on these algorithms. All of the above algorithms were implemented using SGI C++ without any optimization flags. The machine used was an SGI Indigo running IRIX 6.2 with a 195 MHz MIPS R10000 processor, 384 MB RAM and a 32 KB data cache.
13.1
Data Sets
Random Data: The first data set consists of random sets of points. Specifically, given two parameters n, ∆, we generate a random image consisting of n points inside a bounding box of side ∆. It is easy to see that for such an image, the expected ratio of the largest distance to the smallest distance is O(n). To generate the pattern, we use a parameter α, and perturbation parameters ρ and θ. We place a bounding box of side α∆ randomly in the image, and extract the set P0 of all points within. We then perturb each point in P0 by a uniformly distributed random value (in the range [0..ρ]) and then rotate P0 by θ to obtain P. This construction ensures dH (P, Q) ≤ ρ. In all the experiments that follow, we will set θ = 50◦ . We will also set ρ = 2 in all cases where a fixed perturbation is called for. 20
2500
16 CPM MGrid Grid BA
Running Time (secs)
Running Time (secs)
2000
MGrid Grid BA
14
1500
1000
500
12 10 8 6 4 2
0
0 0.5
1
2 3 Size of image (n) * 1000
4
5
Figure 4: Running time vs n
0.5
1
2 3 Size of image (n) * 1000
4
5
Figure 5: Running times vs n (Alignment)
Satellite Data: Mount et al [33] consider data sets obtained by feature extraction from satellite images. We consider two of these sets, one drawn from a satellite image taken over Haifa, Israel, and the other from an image taken over South Africa. We construct patterns from these image sets in the same way as above. The first set (S1) has 1020 points. Its diameter is 173.9. The second set (S2) has size 927 and diameter 188.6. It is interesting to note that the behavior of the algorithms on these sets closely mirrors their behavior on random data sets. As we shall see, these data sets share certain key properties which enable such behavior.
13.2
Running Times
Our first suite of experiments studies the variation in running time of the algorithms as various input parameters are changed. For this paper, we maintain β to be fixed, so the four parameters that determine the running time of the algorithms are k (the pattern size), n (the image size), ², and ∆. In reality, n, ² and ∆ are not really independent of each other. Although ² can in principle vary arbitrarily, typical values of ² will be close to the average closest-point √ distance, which for random data sets can be shown to be Θ(∆/ n). Running time vs image size The first experiment compares the running time of the four schemes as a function of n, the size of the image. For this experiment, we used random instances of cardinality varying from 100 to 1000 √ points. The diameters of the sets were chosen in order to maintain the ratio ∆/ n to be constant. To generate the patterns, we set α = 0.6 in the choice of the random window, and used perturbation parameters ρ = 2 and θ = 50◦ to perturb the points. For each instance, we set ² = 2 and β = 1. For GRID, we recorded the optimal running time achieved by varying the mesh cell size µ². Figure 4 shows how the four schemes compare. Notice that CPM performs orders of magnitude worse than the alignment-based schemes. Hence, in Figure 5, we focus on the three alignment-based schemes. We see that BA performs the worst of the three as n increases, while GRID consistently beats the other two schemes.
GRID BA MGRID CPM
S1 0.412 0.815 0.869 119.792
S2 0.577 0.860 1.265 87.510
Table 2: Times for Satellite data (in seconds)
In Table 2 we show the running times for S1 and S2 using the same pattern extraction parameters as above. As 21
700
3.5 CPM MGrid Grid BA
500
MGrid Grid BA
3 Running Time (secs)
Running Time (secs)
600
400 300 200 100
2.5 2 1.5 1 0.5
0
0 0
5
10
15
20
25
1
2
Epsilon
Figure 6: Running time vs noise
3
4 Epsilon
5
6
7
Figure 7: Running time vs noise (Alignment)
we can see, the relative ordering of the algorithms is the same as above. Running time vs distance threshold Next, we investigate the behavior of the schemes as we vary the noise parameter ². We do this by varying the perturbation parameter ρ, while keeping the input sizes fixed. For this experiment, we fix the (random) image size at 1000, and select the pattern using α = 0.5. As ² increases, the number of solutions increases, and this is reflected in the increased running time of all the alignment-based schemes, as they generate solutions one at a time. However, as the asymptotic complexity of CPM is independent of the number of solutions, its performance is relatively independent of ². This enables it to outperform GRID (and BA) for large ². In Figure 7, we compare the three alignment based schemes. Once again, we observe that GRID outperforms the other two algorithms. Note that as ² increases, more and more candidate pairs (q1 , q2 ) are chosen in Step 3 of the basic alignment algorithm, and so the performance of GRID will tend towards that of BA (which checks all pairs).
13.3
The Quality Of Filtering
The previous suite of experiments establishes the superior performance of GRID for most values of ² and n. The difference between GRID and the other schemes is in its filtering procedure. In this section, we study the filtering performed by GRID more closely. Distance pairs vs distance One of the key parameters governing the running time of GRID is the number of candidate pairs (q1 , q2 ) retrieved in Step 3 of Algorithm 11. For a given point set, we calculate the number of pairs whose distance lies in the range [i, i + 1] for all i < ∆. We scale the ranges such that ∆ = 100, and then normalize each such plot so that the area under each curve is 1. In Figure 8, we plot the resulting distance distributions for a random point set having 1000 points and for S1 and S2. One can immediately observe similarities between the distributions. One common property they all share is that the number of large pairs is very small. This would imply that if the pattern and image sizes are comparable, GRID should perform very well. Another observation is that the distribution curves are quite smooth, without significant peaks at any specific distance. This means that even for smaller patterns the filtering mechanism of GRID should work well. Our next graph demonstrates this behavior. We fix a point set and using the random-window approach described earlier, we generate patterns of varying diameters. We then plot the running time of GRID against the ratio of pattern diameter to image diameter. In Figure 9 we present plots for a random set of 1000 points and diameter 179.2, and S1 and S2. In all cases, the running time decreases as pattern diameter increases. 22
0.025 Random S1 S2
Number of pairs
0.02
0.015
0.01
0.005
0 0
20
40 60 80 Edge length (normalized)
100
Figure 8: Distance distributions
4.5
1.1 S0 S1 S2
4
Running Time (secs)
Running Time (s)
3.5 3 2.5 2 1.5 1
0.9 0.8 0.7 0.6 0.5 0.4
0.5 0 0.2
S0 S1 S2
1
0.3 0.3
0.4 0.5 0.6 0.7 0.8 0.9 Pattern Diameter/Text Diameter
Figure 9: Running time vs. Pattern size
1.0
0
20
40 60 80 Mesh Size / Epsilon
100
Figure 10: Running time vs Mesh size
23
120
Mesh cell size dependence Recall that GRID uses a uniform grid of side µ², where we refer to µ as the mesh factor. Notice that as µ tends to ∆/², GRID will behave more and more like BA, becoming identical to it when µ² = ∆. Therefore, the effectiveness of GRID can be seen by examining the value of µ which yields the best performance. We consider three point sets, S1, S2 and a random set of size 1000. For each set we extract a pattern (where α = 0.5), perturbing it as before. We now run GRID on this instance with varying values of µ, keeping all other parameters fixed. In Figure 10, we plot running time against mesh factor for each of the point sets. Notice that in each case, there is a well-defined non-trivial point at which GRID performs fastest.
13.4
Summary Of Results
The main conclusion of this section is that for typical values of ² GRID is the best choice. Its simplicity allows it to beat the (much more complex) CPM, while its filtering mechanism of matches avoids the quadratic behavior of BA. Moreover, for small ² the match verification runs effectively in constant time and thus the pattern size has an insignificant effect on the running time of the algorithm, making it subquadratic (or even close to linear when the diameter of the pattern is close to the image diameter). On the other hand the CPM scheme, although better in the worst case, is order(s) of magnitude slower than other algorithms for small values of ². However, this ratio changes when the value of ² increases. This is due to a rapid increase of the number of matches or potential matches, which makes the match verification of alignment schemes more costly. Thus for sufficiently large ² CPM achieves the lowest running time. This preliminary study seems to suggest that algorithms based on a reduction to combinatorial pattern matching are primarily of theoretical interest. However, it would be premature to dismiss these methods entirely. It is conceivable that methods employing simpler reductions might well perform as well as purely geometric algorithms.
14
Conclusions
This paper presents schemes for approximate pattern matching in two and three dimensions. We use two main ideas. The first one is a combinatorial analysis of (approximate) distance distributions for point sets in the plane and in three dimensions. We show how these distributions can be used to develop efficient algorithms for point set matching. The algorithms are very simple, drawing on the popular alignment heuristic. Our second contribution is a reduction of geometric pattern matching to combinatorial matching, which allows us to exploit the power of convolutions to obtain near-quadratic time algorithms for geometric pattern matching in two dimensions. In all of the above, a key idea is establishing a connection between the diameter of the point sets and the running time of the algorithms, thus allowing us to escape the pathological settings that drive traditional worst-case analyses of such schemes. We also demonstrate the empirical efficiency of our methods with a detailed study of their behaviour in practice and on average case input. Future Directions In Section 10 we briefly surveyed other empirical methods for pattern matching that are based (primarily) on branch-and-bound search of the transformation space for (near-)optimal solutions. The only direct comparison of these methods against alignment-based schemes was performed by Goodrich et al. [24], who observed that in most cases, a scheme like BA outperformed branch-and-bound methods. The work by Mount et al. attempts to exploit this by employing a hybrid scheme that uses branch-and-bound to quickly eliminate large portions of the transformation space and uses alignment at lower levels of the search procedure. To our knowledge, no rigorous analysis of branch-and-bound schemes has been performed in order to shed light on their behaviour. In addition, our experiments above (especially those addressing the algorithm behaviour as the pattern size tends to the image size) suggest that a way of understanding the behaviour of these methods might come from making assumptions about the range of permissible transformations. In the same spirit, the following problem was suggested by Efrat [15]: 24
Definition 14.1 (Dog in the Field). Given two point sets P, Q in the plane and ² > 0 such that there exists a unique transformation T such that d(T (P), Q) ≤ ², determine T in time O(n + kc ), where |P| = k, |Q| = n, c ≥ 1. The formulation is motivated by the following idea: suppose we have a picture of a dog in a field. Given a picture of the dog, we would like to identify the location of the dog in the picture. There is only one dog, and the rest of the picture can be considered to be noise. Is there some way of exploiting this information to obtain an algorithm that is ˜ 2 k3 ) runs in time linear in the noise ? Notice that known techniques will give us an algorithm running in time O(n (for d = dH ), and so under the assumption that k ¿ n, a solution satisfying the desired form will yield tremendous savings over the worst-case bound. We would like to thank David Mount for providing us with the data that was used in [33]. We also would like to thank Victor Shoup for providing us with the NTL package for finite field arithmetic, Dragomir Angelov for helpful discussions, and the anonymous reviewers for helping clarify the presentation of the paper.
References [1] AGARWAL , P., A RONOV, B., S HARIR , M., AND S URI , S. Selecting distances in the plane. Algorithmica 9, 5 (1993), 495–514. [2] A KUTSU , T., TAMAKI , H., AND T OKUYAMA , T. Distribution of distances and triangles in a point set and algorithms for computing the largest common point set. In Proc. 13th Annual ACM Symposium on Computational Geometry (1997), pp. 314–323. [3] A LT, H., AND G UIBAS , L. Discrete geometric shapes: Matching, interpolation, and approximation. A survey. In Handbook of Computational Geometry, J.-R. Sack and J. Urrutia, Eds. Elsevier Science Publishers B.V. North-Holland, 1999, pp. 121–153. [4] BAIRD , H. S. Model-based image matching using location. Distinguised Dissertation Series. MIT Press, 1985. [5] BALLARD , D. Generalizing the Hough Transform to detect arbitrary shapes. Pattern Recognition 13, 2 (1981), 111–122. [6] B OXER , L. Point set pattern matching in 3-D. Pattern Recognition Letters 17, 12 (1996), 1293–1297. [7] B RAβ, P. Exact point pattern matching and the number of congruent triangles in a three dimensional point set. In ESA 2000 — European Symposium on Algorithms (2000), M. Paterson, Ed., vol. 1879 of Lecture Notes in Computer Science, Springer-Verlag, pp. 112–119. [8] C ARDOZE , D., AND S CHULMAN , L. Pattern matching for spatial point sets. In Proc. 39th Annual Symposium on the Foundations of Computer Science (November 1998), IEEE, pp. 156–165. [9] C HEW, L., G OODRICH , M., H UTTENLOCHER , D., K EDEM , K., K LEINBERG , J., AND K RAVETS , D. Geometric pattern matching under euclidean motion. In Proc. Fifth Canadian Conference on Computational Geometry (1993), pp. 151–156. [10] C LARKSON , K., E DELSBRUNNER , H., G UIBAS , L., S HARIR , M., AND W ELZL , E. Combinatorial complexity bounds for arrangements of curves and spheres. Discrete Computational Geometry 5 (1990), 99–160. [11] C OLE , R., AND H ARIHARAN , R. Tree pattern matching and subset matching in randomized O(n log3 m) time. In Proc. 29th Annual ACM Symposium on Theory of Computing (1997), pp. 66–75. [12] C OLE , R., H ARIHARAN , R., AND I NDYK , P. Tree pattern matching and subset matching in deterministic O(n log3 m) Time. In Proc. 10th Annual ACM-SIAM Symposium on Discrete Algorithms (1999), pp. 245– 254. 25
[13] D OBKIN , D., 181–186.
AND
L IPTON , R. Multidimensional search problems. SIAM Journal on Computing 5 (1976),
[14] E DELSBRUNNER , H., VALTR , P., AND W ELZL , E. Cutting dense point sets in half. Discrete Comp. Geom 17 (1997), 243–255. [15] E FRAT, A. The dog in the field problem. Personal communication, 1998. [16] E FRAT, A. The complexity of the union of (alpha, beta)-covered objects. In Proc. 15th ACM Symposium on Computational Geometry (1999), pp. 134–142. [17] E FRAT, A., I NDYK , P., AND V ENKATASUBRAMANIAN , S. Pattern matching for sets of segments. In Proc. 12th ACM-SIAM Symposium on Discrete Algorithms (2001), pp. 295–304. ˝ , P. On sets of distances of n points. American Mathematical Monthly 53 (1946), 248–250. [18] E RD OS ˝ , P. On some applications of graph theory to geometry. Canadian J. Mathematics 19 (1967), 968–971. [19] E RD OS ˝ , P., M AKAI , E., PACH , J., AND S PENCER , J. Gaps in difference sets, and the graph of nearly [20] E RD OS equal distances. Applied Geometry and Discrete Mathematics: The Victor Klee Festschrift (P. Gritzmann and B. Sturmfels, eds.) 4 (1991), 265–273. [21] E RIC , W., AND G RIMSON , L. Object Recognition by Computer: The Role of Geometric Constraints. MIT Press, 1990. [22] E RICKSON , J. Dense point sets have sparse delaunay triangulations. In Proc. 13th ACM-SIAM Symposium on Discrete Algorithms (2002), pp. 96–105. [23] G AVRILOV, M., I NDYK , P., M OTWANI , R., AND V ENKATASUBRAMANIAN , S. Geometric pattern matching: A performance study. In Proc. 15th ACM Symposium on Computational Geometry (1999), pp. 79–85. [24] G OODRICH , M., M ITCHELL , J., AND O RLETSKY, M. Practical methods for approximate geometric pattern matching under rigid motions. In Proc. 10th Annual ACM Symposium on Computational Geometry (1994), pp. 103–113. [25] H AGEDOORN , M., AND V ELTKAMP, R. C. Reliable and efficient pattern matching using an affine invariant metric. International Journal of Computer Vision 31, 2/3 (1999), 203–225. [26] H EFFERNAN , P. J., AND S CHIRRA , S. Approximate decision algorithms for point set congruence. Computational Geometry: Theory and Applications 4, 3 (1994), 137–156. [27] H OLM , L., AND S ANDER , C. Mapping the protein universe. Science 273, 5275 (August 1996), 595–602. [28] H UTTENLOCHER , D. P., K EDEM , K., AND S HARIR , M. The upper envelope of Voronoi surfaces and its applications. Discrete Comput. Geom. 9 (1993), 267–291. [29] H UTTENLOCHER , D. P., AND RUCKLIDGE , W. T. A multi-resolution technique for comparing images using the hausdorff distance. In Proc. IEEE Conference on Computer Vision and Pattern Recognition (June 1993), IEEE, pp. 705–706. [30] I NDYK , P. Faster algorithms for string matching problems: matching the convolution bound. In Proc. 38th IEEE Symposium on Foundations of Computer Science (1998), pp. 166–173. [31] I NDYK , P., M OTWANI , R., AND V ENKATASUBRAMANIAN , S. Geometric matching under noise: Combinatorial bounds and algorithms. In Proc. 10th Annual SIAM-ACM Symposium on Discrete Algorithms (1999), pp. 457–465. 26
[32] I NDYK , P., AND V ENKATASUBRAMANIAN , S. Approximate congruence in near-linear time. International Journal of Computational Geometry and Applications (2002). To Appear. [33] M OUNT, D., N ETANYAHU , N., AND L E M OIGNE , J. Improved algorithms for robust point pattern matching and applications to image registration. In Proc. 14th ACM Symposium on Computational Geometry (June 1998), pp. 155–164. [34] P REPARATA , F. P., York., 1985.
AND
S HAMOS , M. I. Computational Geometry: An Introduction. Springer-Verlag, New
[35] ROSENFELD , A. Image analysis and computer vision: 1998. Computer Vision and Image Understanding 74, 1 (1998), 36–95. [36] S HOUP, V. NTL: A Library for doing Number Theory. http://www.cs.wisc.edu/˜shoup/ntl/. [37] S Z E´ KELY, L. Crossing numbers and hard Erd˝os problems in discrete geometry. Combinatorics, Probability and Computing 6, 3 (1997), 353–358. [38] S ZEMEREDI , E., 381–392.
AND
T ROTTER , W. T. Extremal problems in discrete geometry. Combinatorica 3 (1983),
[39] VALTR , P. Planar point sets with bounded ratios of distances. PhD thesis, Fachbereich Mathematik, Freie Universit¨at Berlin, Berlin, Germany, 1994. [40] VALTR , P. Lines, line-point incidences and crossing families in dense sets. Combinatorica 16 (1996), 269–294. [41] V ERBARG , K. Approximate center points in dense sets. In Abstracts 13th European Workshop Comp. Geom. (1997), pp. 55–56. [42] Z HOU , Y., 833–857.
AND
S URI , S. Analysis of a bounding box heuristic for object intersection. J. ACM 46 (1999),
27
