《Computational Geometry》lecnotes 21 Prof, i .k. atkikalak.pdf_大学文库

13.472J/1.128J/2.158J/16.940J COMPUTATIONAL GEOMETRY Lecture 21 Dr. K. H. Ko Prof. N. M. Patrikalakis Copyright c 2003 Massachusetts Institute of Technology Contents 21 Object Matching 2 21.1 Various matching methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 21.2 Methods through localization/registration . . . . . . . . . . . . . . . . . . . . 2 21.3 Classification of matching methods . . . . . . . . . . . . . . . . . . . . . . . . 2 21.4 Problem statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 21.4.1 Distance metric . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 21.4.2 Distance between a point and a parametric surface . . . . . . . . . . . 4 21.4.3 Distance metric function . . . . . . . . . . . . . . . . . . . . . . . . . . 4 21.5 Matching problems : CGWOS, CPWOS, CGWS or CPWS . . . . . . . . . . . 5 21.5.1 Resolving the scaling effects . . . . . . . . . . . . . . . . . . . . . . . . 5 21.5.2 Rotation and translation . . . . . . . . . . . . . . . . . . . . . . . . . . 5 21.6 Matching problems : IGWOS, IPWOS, IGWS or IPWS . . . . . . . . . . . . . 6 21.6.1 Iterative Closest Point (ICP) algorithm [1] for IGWOS or IPWOS . . . 6 21.6.2 ICP algorithm for scaling effects . . . . . . . . . . . . . . . . . . . . . . 6 21.7 Matching problems : NGWOS or NPWOS . . . . . . . . . . . . . . . . . . . . 7 21.7.1 Search method [2] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 21.7.2 KH method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 21.8 Matching problems : NGWS . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 21.9 Matching problems : NPWS . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 21.9.1 Umbilical point method [7] . . . . . . . . . . . . . . . . . . . . . . . . . 9 21.9.2 Optimization method [7] . . . . . . . . . . . . . . . . . . . . . . . . . . 12 21.10Matching problems : offset method . . . . . . . . . . . . . . . . . . . . . . . . 13 21.10.1Distance function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 21.10.2Objective function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 21.10.3Gradient vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Bibliography 16 1

the matching of part of objects. Matching problems can be further categorized based on the ava ailability of correspondence or initial transformation information between two obj ects a the application of scaling. The classification of matching problems is summarized in Table 21.1. In the table, acronyms are used for simplification as follows C: Correspondence information is provided I: Initial information on correspondence is provided .N: No correspondence information is available ·P: Partial matching ·G: Global matching ●WOS: Without scali ·Ws: With scaling Global matching Partial matching Criteria Without scaling With scaling Without scaling With scaling Correspondence CGWOS CGWS CPWOS CPWS nformation Initial information IGWOS IPWOS IPWS No information NGWOS NGWS NPWOS NPWS Table 21.1: Classification of matching problems I When correspondence information is provided, which is one of the types CGWOS or CP- OS, then a matching problem is simply reduced to calculation of the rigid body trans- formation 3, 4. If no correspondence is known, but a good initial approximation for the transformation is available(IGWOS or IPWOS), then popular iterative schemes such as the Iterative Closest Point(ICP)algorithm [1]can be employed. However, when no prior clue for correspondence or transformation is given(NGWOS or NPWOS), the matching problem be comes more complicated. In this case, the solution process must provide a means to establish such correspondence information such as in 2 Scaling is another factor that needs to be considered separately. If a matching problem involves scaling effects, then direct comparison of quantitative measures cannot be used any longer. For the global matching case, a scaling factor can be estimated by the ratio of surface areas and applied to resolve the scaling transformation. However, when it comes to partial matching, such area information becomes useless for the scaling factor estimation. When the correspondence information between two objects is known(CGWS or CPWS), the scaling factor between the objects can be easily obtained by using the ratio of Euclidean distances between two sets of corresponding points or areas, or the ratio of the principal curvatures. If an initial scaling value as well as a good initial approximation is provided (IGWS or IPWS) the ICP algorithm by Besl [1] or other optimization schemes such as the quasi-Newton method

the matching of part of objects. Matching problems can be further categorized based on the availability of correspondence or initial transformation information between two objects and the application of scaling. The classification of matching problems is summarized in Table 21.1. In the table, acronyms are used for simplification as follows: • C : Correspondence information is provided. • I : Initial information on correspondence is provided. • N : No correspondence information is available. • P : Partial matching. • G : Global matching. • WOS : Without scaling. • WS : With scaling. Global matching Partial matching Criteria Without scaling With scaling Without scaling With scaling Correspondence information CGWOS CGWS CPWOS CPWS Initial information IGWOS IGWS IPWOS IPWS No information NGWOS NGWS NPWOS NPWS Table 21.1: Classification of matching problems When correspondence information is provided, which is one of the types CGWOS or CPWOS, then a matching problem is simply reduced to calculation of the rigid body transformation [3, 4]. If no correspondence is known, but a good initial approximation for the transformation is available (IGWOS or IPWOS), then popular iterative schemes such as the Iterative Closest Point (ICP) algorithm [1] can be employed. However, when no prior clue for correspondence or transformation is given (NGWOS or NPWOS), the matching problem becomes more complicated. In this case, the solution process must provide a means to establish such correspondence information such as in [2]. Scaling is another factor that needs to be considered separately. If a matching problem involves scaling effects, then direct comparison of quantitative measures cannot be used any longer. For the global matching case, a scaling factor can be estimated by the ratio of surface areas and applied to resolve the scaling transformation. However, when it comes to partial matching, such area information becomes useless for the scaling factor estimation. When the correspondence information between two objects is known (CGWS or CPWS), the scaling factor between the objects can be easily obtained by using the ratio of Euclidean distances between two sets of corresponding points or areas, or the ratio of the principal curvatures. If an initial scaling value as well as a good initial approximation is provided (IGWS or IPWS), the ICP algorithm by Besl [1] or other optimization schemes such as the quasi-Newton method 3

[12] can be employed. The problem of NGWS type can be solved by the moment method using the principal moment of inertia and ratio of areas or volumes. No attempt, however, has been made to solve the problems of NPWS type. For a more detailed overview, see Ko [6]. 21.4 Problem statement 21.4.1 Distance metric The Euclidean distance between two points p1 and p2 is defined as de(p1, p2) = |p1 − p2|. (21.1) We also define the minimum distance between a surface r and a point p as follows: dsp(r, p) = min{de(p, pi), ∀pi ∈ r}. (21.2) 21.4.2 Distance between a point and a parametric surface Let us assume that we have a point p and a parametric surface r = r(u, v), 0 ≤ u, v ≤ 1. Then the squared distance function is defined as follows: D(u, v) = |p − r(u, v)| 2 , = (p − r(u, v)) · (p − r(u, v)). (21.3) Finding the minimum distance between p and r is reduced to minimizing (21.3) within the square 0 ≤ u, v ≤ 1. Therefore, the problem needs to be broken up into several sub-problems which consider the behavior of the distance function at the boundary and in the interior of the bound [10]. The sub-problems are summarized as follows: Find the minimum distances (1) in the interior domain, (2) along the boundary curves and (3) from the corner points. Among those minimum distances, the smallest one is chosen as the minimum distance between the point p and the surface r. A robust calculation of the minima of the distance function (21.3) can be achieved by the Interval Projected Polyhedron (IPP) algorithm [13, 10, 14]. 21.4.3 Distance metric function A function can be defined using the squared distance function (21.3) to formulate a matching problem. Suppose that we have a NURBS surface rB and an object rA represented in discrete points or surfaces. Then, a matching problem can be stated as finding the rigid body transformation (a translation vector t and a rotation matrix R) so that a global distance metric function Φ = X ∀p∈rA dsp(rB,(σRp + t)) (21.4) becomes minimal, where σ is a scaling factor. 4

21.5 Matching problems: CGWOS, CPWOS, CGWs or CPWs 21.5.1 Resolving the scaling effects Matching problems of CG WS or CPWS type involve the scaling effects. Therefore, a scaling factor has to be estimated so that the scaling transformation is performed, before calculatin the rigid body transformation. Since correspondence information between two objects are available, a scaling factor can be estimated by using the ratio of the principal curvatures at the corresponding points. After scaling has been resolved, the matching problems of CGWS or CPws type are treated as those of CG wos or CPwos type 21.5.2 Rotation and translation Suppose that we have two 3-tuples m; and n; (i= 1, 2, 3), and the correspondence information for each point. From these points, a translation vector and a rotation matrix can be calculated The translation vector is easily obtained by using the centroids of each 3-tuple. The centroids Cm and cn are given b ni (21.5) and the difference between Cm and cn becomes the translation vector tT Cn-Cm. A rotation matrix consists of three unknown components(the Euler angles). Since the two 3- tuples provide nine constraints, the rotation matrix may be constructed by using some of the constraints. But the results could be different if the remaining constraints are used for the rotation matrix calculation 3. In order to use all the constraints equally, the least squares method may be employed 3]. The basic solution by Horn 3 is described below. Suppose that the translation has been performed. Then what is left is to find the rotation matrix R so that ∑|;-(Rm)|2 n-2∑n;(Rm)+ (21.6) is minimized. Here, D=E_i ni(Rmi has to be maximized to minimize ' The problem can be solved in the quaternion framework. a quaternion can be considered as a vector with four components, i. e. a vector part in 3D and a scalar part. A rotation can be equivalently defined as a unit quaternion q=cos(), sin()ax, sin()ay, sin()a2 which represents a rota tion movement around (az, ay, as)by degree. In the quaternion framework, the problem is reduced to the eigenvalue problem of the 4 x 4 matrix H obtained from the correlation matrix M S11+S22+S H S11-822-S33 S12+s2122-S11-533 S31+S13 S23+S32

21.5 Matching problems : CGWOS, CPWOS, CGWS or CPWS 21.5.1 Resolving the scaling effects Matching problems of CGWS or CPWS type involve the scaling effects. Therefore, a scaling factor has to be estimated so that the scaling transformation is performed, before calculating the rigid body transformation. Since correspondence information between two objects are available, a scaling factor can be estimated by using the ratio of the principal curvatures at the corresponding points. After scaling has been resolved, the matching problems of CGWS or CPWS type are treated as those of CGWOS or CPWOS type. 21.5.2 Rotation and translation Suppose that we have two 3-tuples mi and ni(i = 1, 2, 3), and the correspondence information for each point. From these points, a translation vector and a rotation matrix can be calculated. The translation vector is easily obtained by using the centroids of each 3-tuple. The centroids cm and cn are given by cm = 1 3 X 3 i=1 mi , cn = 1 3 X 3 i=1 ni , (21.5) and the difference between cm and cn becomes the translation vector tT = cn − cm. A rotation matrix consists of three unknown components (the Euler angles). Since the two 3- tuples provide nine constraints, the rotation matrix may be constructed by using some of the constraints. But the results could be different if the remaining constraints are used for the rotation matrix calculation [3]. In order to use all the constraints equally, the least squares method may be employed [3]. The basic solution by Horn [3] is described below. Suppose that the translation has been performed. Then what is left is to find the rotation matrix R so that Φ 0 = X 3 i=1 |ni − (Rmi)| 2 = X 3 i=1 |ni | 2 − 2 X 3 i=1 ni · (Rmi) + X 3 i=1 |Rmi | 2 (21.6) is minimized. Here, D = P3 i=1 ni · (Rmi) has to be maximized to minimize Φ 0 . The problem can be solved in the quaternion framework. A quaternion can be considered as a vector with four components, i.e. a vector part in 3D and a scalar part. A rotation can be equivalently defined as a unit quaternion qˇ = h cos( θ 2 ),sin( θ 2 )ax,sin( θ 2 )ay,sin( θ 2 )az i which represents a rotation movement around (ax, ay, az) by θ degree. In the quaternion framework, the problem is reduced to the eigenvalue problem of the 4 × 4 matrix H obtained from the correlation matrix M: H =      s11 + s22 + s33 s23 − s32 s31 − s13 s12 − s21 s23 − s32 s11 − s22 − s33 s12 + s21 s31 + s13 s31 − s13 s12 + s21 s22 − s11 − s33 s23 + s32 s12 − s21 s31 + s13 s23 + s32 s33 − s22 − s11      , (21.7) 5

where M = X 3 i=1 nimT i =    s11 s12 s13 s21 s22 s23 s31 s32 s33    . (21.8) The eigenvector corresponding to the maximum positive eigenvalue is a quaternion which minimizes the equation (21.6). An orthonormal rotation matrix R can be recovered from a quaternion qˇ = [q0, q1, q2, q3] by R =    q 2 0 + q 2 1 − q 2 2 − q 2 3 2(q1q2 − q0q3) 2(q1q3 + q0q2) 2(q1q2 + q0q3) q 2 0 + q 2 2 − q 2 1 − q 2 3 2(q2q3 − q0q1) 2(q1q3 − q0q2) 2(q2q3 + q0q1) q 2 0 + q 2 3 − q 2 1 − q 2 2    . (21.9) The procedure described above can also be applied to the case where more than three corresponding point pairs are provided. 21.6 Matching problems : IGWOS, IPWOS, IGWS or IPWS 21.6.1 Iterative Closest Point (ICP) algorithm [1] for IGWOS or IPWOS Algorithm • The point set P with Np points {p~i} from the data shape and the model shape X (with Nx supporting geometric primitives: points, lines, or triangles) are given. • The iteration is initialized by setting P0 = P, q~0 = [1, 0, 0, 0, 0, 0, 0]T and k = 0. The registration vectors are defined relative to the initial data set P0 so that the final registration represents the complete transformation. Steps 1,2,3, and 4 are applied until convergence within a tolerance τ . a. Compute the closest points: Yk = C(Pk, X). b. Compute the registration: (~qk, dk) = Q(P0, Yk). c. Apply the registration: Pk+1 = ~qk(P0). d. Terminate the iteration when the change in mean-square error falls below a preset threshold τ > 0 specifying the desired precision of the registration: |dk − dk+1| < τ . 21.6.2 ICP algorithm for scaling effects When initial information on transformation is provided, the ICP method can be extended to resolve scaling effects in the matching problem. A scaling factor σ is included in the objective function (21.4). The scaling transformation is performed at step c in the ICP algorithm. In this case we have to provide seven initial values (three for translation, three for rotation and one for scaling). 6

21.7 Matching problems: NG WOs or NPWOs 21.7.1 Search method 2 Chua and Jarvis 2 developed a method to align two objects through registration assuming no prior knowledge of correspondence between two range data sets. They use a bi-quadratic polynomial to fit data points in the local area in the least squares sense and calculate the principal curvatures and Darboux frames. Then they construct a list of sensed data points based on the fit error. Three seed points are selected to form the list such that the area of the triangle represented by the seed points becomes maximized to reduce any mismatch error. Three constraints(curvature, distance and direction) are imposed to sort out possible corresponding points out of the model data set. Then a list of transformations can be obtained from the candidate points and an optimum transformation is selected. Various searching algorithms are described and demonstrated in[2 21.7.2 Kh method The overall diagram of the KH method [8 is shown in Figure 21.1. The input of the process includes two objects and three pairs of the gaussian and the mean curvatures at three different non-collinear locations. The algorithm yields the minimum value of in the equation(21. 4) and the corresponding rotation matrix R and the translation vector t. Since no scaling effect is involved. we assume that a scaling factor g=1 Step 10 Step 10 is to select three non-collinear points m1, m2 and mg on ri away from the boundary of rI where each point has different, Gaussian K and mean curvature H values. At mi, we have Ki and Hi, where i= 1, 2, 3. Next, subdivide r2 into rational Bezier surface patches B, G=l,.,n)by inserting appropriate knots [ 5, 11]. Then for each rational Bezier surface patch B,, we express K, and H, in the bivariate rational Bernstein polynomial basis using rounded interval arithmetic to formulate the problem. This allows us to use the Interval Projected Polyhedron(IPP) algorithm [10, 13 for solving nonlinear polynomial systems. For ach pair Ki and Hi, we solve the following system of equations by the IPp technique. K(u,v)=K±6 H1(u,v)=H2±6H,(j=1,……, n and i=1,2,3) (21.10) where dk and H represent the uncertainty of estimated curvatures from data points. For each pair of Ki and Hi, a list of roots Li=(uik, Uik),(k=l, .. di) is obtained Step 12(selection process) A simple pruning search based on the Euclidean distance can be applied to the selection process. We have three lists of candidate points, Li=(uik, Vik),(k= 1, . di)from which one 3-tuple(n1, n2, n3) 2(v2k,t2k), (21.11)

21.7 Matching problems : NGWOS or NPWOS 21.7.1 Search method [2] Chua and Jarvis [2] developed a method to align two objects through registration assuming no prior knowledge of correspondence between two range data sets. They use a bi-quadratic polynomial to fit data points in the local area in the least squares sense and calculate the principal curvatures and Darboux frames. Then they construct a list of sensed data points based on the fit error. Three seed points are selected to form the list such that the area of the triangle represented by the seed points becomes maximized to reduce any mismatch error. Three constraints (curvature, distance and direction) are imposed to sort out possible corresponding points out of the model data set. Then a list of transformations can be obtained from the candidate points and an optimum transformation is selected. Various searching algorithms are described and demonstrated in [2]. 21.7.2 KH method The overall diagram of the KH method [8] is shown in Figure 21.1. The input of the process includes two objects and three pairs of the Gaussian and the mean curvatures at three different non-collinear locations. The algorithm yields the minimum value of Φ in the equation (21.4), and the corresponding rotation matrix R and the translation vector t. Since no scaling effect is involved, we assume that a scaling factor σ = 1. Step 10 Step 10 is to select three non-collinear points m1, m2 and m3 on r1 away from the boundary of r1 where each point has different, Gaussian K and mean curvature H values. At mi , we have Ki and Hi , where i = 1, 2, 3. Next, subdivide r2 into rational B´ezier surface patches Bj (j = 1, · · · , n) by inserting appropriate knots [5, 11]. Then for each rational B´ezier surface patch Bj , we express Kj and Hj in the bivariate rational Bernstein polynomial basis using rounded interval arithmetic to formulate the problem. This allows us to use the Interval Projected Polyhedron (IPP) algorithm [10, 13] for solving nonlinear polynomial systems. For each pair Ki and Hi , we solve the following system of equations by the IPP technique. Kj(u, v) = Ki ± δK, Hj(u, v) = Hi ± δH, (j = 1, · · · , n and i = 1, 2, 3), (21.10) where δK and δH represent the uncertainty of estimated curvatures from data points. For each pair of Ki and Hi , a list of roots Li = (uik, vik), (k = 1, · · · , di) is obtained. Step 12 (selection process) A simple pruning search based on the Euclidean distance can be applied to the selection process. We have three lists of candidate points, Li = (uik, vik), (k = 1, · · · , di) from which one 3-tuple (n1, n2, n3) n1 = r2(u1k, v1k), n2 = r2(u2k, v2k), n3 = r2(u3k, v3k), (21.11) 7

Figure 21. 2: Localized surfaces Examples are presented to demonstrate the proposed algorithms. Solids bounded by bicu bic integral B-spline surface patches A and B are used. Solid A is enclosed in a rectangular box of 25mm x 23 48mm x 11mm. Here. the height of solid a is 25mm. Figure 21. 3 shows a sequence of operations for matching of the two surfaces using the principal moments of inertia of input solids. In this example, for clarity, only part of the boundary surfaces of the solids are displayed. The smaller solid has been translated, rotated, uniformly scaled and reparam eterized. In Figure 213-(A), two boundary surfaces of the input solids are shown with their control points. Those two surfaces have similar shape but different numbers of control points and parametrization. Matching the centroids of the two solids is performed by translating the small solid by the position difference between the centroids, which is demonstrated in Figure 213-(B). The orientation of the largest principal moment of inertia of the solid A is aligned to that of the largest one of the solid B. Similarly, the remaining two orientations are aligne based on the values of the principal moments of inertia. After matching the orientations of the principal moment of inertia, the two solids are aligned in their orientations as shown in Figure 213-(C). Figure 21 3-(D)shows that the two solids match after uniform scaling obtained from the ratio between the volumes of the two solids, 4.651, is applied to the small solid 21 9 Matching problems: NPWs 21.9.1 Umbilical point method [7 The correspondence search for matching problems of NGWOS or NPWos type only deals with qualitative aspects. Since the w-plane is not affected by scaling, only qualitative correspondence can be established in the process. This implies that without a scaling factor applied, a rigid body transformation cannot be obtained for aligning two surfaces. Therefore, a scaling factor has to be estimated before any transformation is considered

Figure 21.2: Localized surfaces Examples are presented to demonstrate the proposed algorithms. Solids bounded by bicubic integral B-spline surface patches A and B are used. Solid A is enclosed in a rectangular box of 25mm × 23.48mm × 11mm. Here, the height of solid A is 25mm. Figure 21.3 shows a sequence of operations for matching of the two surfaces using the principal moments of inertia of input solids. In this example, for clarity, only part of the boundary surfaces of the solids are displayed. The smaller solid has been translated, rotated, uniformly scaled and reparameterized. In Figure 21.3-(A), two boundary surfaces of the input solids are shown with their control points. Those two surfaces have similar shape but different numbers of control points and parametrization. Matching the centroids of the two solids is performed by translating the small solid by the position difference between the centroids, which is demonstrated in Figure 21.3-(B). The orientation of the largest principal moment of inertia of the solid A is aligned to that of the largest one of the solid B. Similarly, the remaining two orientations are aligned based on the values of the principal moments of inertia. After matching the orientations of the principal moment of inertia, the two solids are aligned in their orientations as shown in Figure 21.3-(C). Figure 21.3-(D) shows that the two solids match after uniform scaling obtained from the ratio between the volumes of the two solids, 4.651, is applied to the small solid. 21.9 Matching problems : NPWS 21.9.1 Umbilical point method [7] The correspondence search for matching problems of NGWOS or NPWOS type only deals with qualitative aspects. Since the ω-plane is not affected by scaling, only qualitative correspondence can be established in the process. This implies that without a scaling factor applied, a rigid body transformation cannot be obtained for aligning two surfaces. Therefore, a scaling factor has to be estimated before any transformation is considered. 9