13472J/1.128J/2158J/16940J COMPUTATIONAL GEOMETRY Lectures 10-12 N.M. Patrikalakis Massachusetts Institute of Technology Cambridge MA 02139-4307 USA Copyright 2003 Massachusetts Institute of Technology Contents 10.1 Overview of intersection problems 10.2 Intersection problem classification 5 10. 2. 1 Classification by dimension 10.2.2 Classification by type of geometr 10.2.3 Classification by number system 6 10.3 Point /point"intersection 7 10.4 Point/curve intersection 10.4.1 Point/Implicit curve intersection 10.4.2 Point /Parametric curve intersection 10.4.3 Point/Procedural parametric (offset, evolute, etc. curve intersection 10.5 Point /surface intersection 10.5.1 Point/Implicit(usually algebraic)surface intersection 10.5.2 Point/Rational polynomial surface intersection 13 10.5.3 Point /Procedural surface intersection 10.6 Curve/curve intersection 20 10.6.1 Case D3: RPP/IA curve intersection 10.6.2 Case D1: RPP/RPP Curve Intersection 10.6.3 Case D2 /D5: RPP/PP and PP/PP Curve Intersections 10.6.4 Case D6: PP/IA Curve Intersection 10.6.5 Case D8: IA/IA Curve Intersection 10.7 Curve/ surface intersection 10.7. 1 Case E3: RPP Curve/IA Surface Intersection 10.7.2 Case El: RPP Curve/RPP Surface Intersection 10.7.3 Case E2/E6: RPP/PP, PP/PP Curve/Surface Intersection 10.7.4 Case E7: PP Curve/IA Surface Intersection
13.472J/1.128J/2.158J/16.940J COMPUTATIONAL GEOMETRY Lectures 10 - 12 N. M. Patrikalakis Massachusetts Institute of Technology Cambridge, MA 02139-4307, USA Copyright c 2003 Massachusetts Institute of Technology Contents 10.1 Overview of intersection problems . . . . . . . . . . . . . . . . . . . . . . . . . 3 10.2 Intersection problem classification . . . . . . . . . . . . . . . . . . . . . . . . . 5 10.2.1 Classification by dimension . . . . . . . . . . . . . . . . . . . . . . . . 5 10.2.2 Classification by type of geometry . . . . . . . . . . . . . . . . . . . . . 5 10.2.3 Classification by number system . . . . . . . . . . . . . . . . . . . . . . 6 10.3 Point/point “intersection” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 10.4 Point/curve intersection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 10.4.1 Point/Implicit curve intersection . . . . . . . . . . . . . . . . . . . . . 8 10.4.2 Point/Parametric curve intersection . . . . . . . . . . . . . . . . . . . . 10 10.4.3 Point/Procedural parametric (offset, evolute, etc.) curve intersection . 12 10.5 Point/surface intersection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 10.5.1 Point/Implicit (usually algebraic) surface intersection . . . . . . . . . . 13 10.5.2 Point/Rational polynomial surface intersection . . . . . . . . . . . . . . 13 10.5.3 Point/Procedural surface intersection . . . . . . . . . . . . . . . . . . . 19 10.6 Curve/curve intersection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 10.6.1 Case D3: RPP/IA curve intersection . . . . . . . . . . . . . . . . . . . 20 10.6.2 Case D1: RPP/RPP Curve Intersection . . . . . . . . . . . . . . . . . 27 10.6.3 Case D2/D5: RPP/PP and PP/PP Curve Intersections . . . . . . . . . 28 10.6.4 Case D6: PP/IA Curve Intersection . . . . . . . . . . . . . . . . . . . . 28 10.6.5 Case D8: IA/IA Curve Intersection . . . . . . . . . . . . . . . . . . . . 29 10.7 Curve/surface intersection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 10.7.1 Case E3: RPP Curve/IA Surface Intersection . . . . . . . . . . . . . . 30 10.7.2 Case E1: RPP Curve/RPP Surface Intersection . . . . . . . . . . . . . 31 10.7.3 Case E2/E6: RPP/PP, PP/PP Curve/Surface Intersection . . . . . . . 31 10.7.4 Case E7: PP Curve/IA Surface Intersection . . . . . . . . . . . . . . . 31 1
10.7.5 Case Ell: IA Curve/IA Surface Intersection 10.7.6 IA Curve/RPP Surface Intersection 10.8 Surface/Surface Intersections 2333 10.8.1 Case F3: RPP/IA Surface Intersection 10.8.2 Case F1: RPP/RPP Surface Intersection 10.8.3 Case F8: IA/IA Surface Intersection 10.9 Nonlinear Solvers 51 10.9.1 Motivation 10.1OLocal Solution Methods 10.11 Classification of Global Solution Methods 10.12Subdivision Method(Projected Polyhedron Method) 10.13Interval Method 10.13.1 Motivation 60 10.132 Definitions 10. 13. 3 Interval Arithmetic 62 10.13.4 Algebraic Properties 10. 13.5 Rounded Interval Arithmetic and its Implementati 10. 14Interval Projected Polyhedron Algorithm 10.15Interval Newton method Bibliography Reading in the Textbook Chapters 4 and 5, pp. 73-160
10.7.5 Case E11: IA Curve/IA Surface Intersection . . . . . . . . . . . . . . . 32 10.7.6 IA Curve/RPP Surface Intersection . . . . . . . . . . . . . . . . . . . . 33 10.8 Surface/Surface Intersections . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 10.8.1 Case F3: RPP/IA Surface Intersection . . . . . . . . . . . . . . . . . . 34 10.8.2 Case F1: RPP/RPP Surface Intersection . . . . . . . . . . . . . . . . . 42 10.8.3 Case F8: IA/IA Surface Intersection . . . . . . . . . . . . . . . . . . . 46 10.9 Nonlinear Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 10.9.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 10.10Local Solution Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 10.11Classification of Global Solution Methods . . . . . . . . . . . . . . . . . . . . . 54 10.12Subdivision Method (Projected Polyhedron Method) . . . . . . . . . . . . . . 55 10.13Interval Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 10.13.1Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 10.13.2Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 10.13.3 Interval Arithmetic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 10.13.4Algebraic Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 10.13.5Rounded Interval Arithmetic and its Implementation . . . . . . . . . . 62 10.14Interval Projected Polyhedron Algorithm . . . . . . . . . . . . . . . . . . . . . 65 10.15Interval Newton method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 Bibliography 72 Reading in the Textbook • Chapters 4 and 5, pp. 73 - 160 2
Lectures 10-12 Intersection Problems, Nonlinear Solvers and robustness issues 10.1 Overview of intersection problems Intersections are fundamental in computational geometry, geometric modeling and design analysis and manufacturing applications. Examples of intersection problems include Shape interrogation(eg. visualization) through contouring (intersection with a series of parallel planes, coaxial cylinders, and cones etc. Numerical control machining (milling) involving intersection of offset surfaces with a series of parallel planes, to create machining paths for ball(spherical) cutters Representation of complex geometries in the "Boundary Representation"scheme; for example, the description of the internal geometry and of structural members cars airplanes, ships, etc, involves Intersections of free-form parametric surfaces with low order algebraic surfaces planes, quadrics, torii Intersections of low order algebraic surfaces in a process called boundary evaluation, in which the Boundary Representation is cre- ated by "evaluating " the Constructive Solid Geometry model of the object. during this process, intersections of the surfaces of primitives(see Figure 10.1)must be found during Boolean operations Boolean opertations on point sets A, B include(see Figure 10.2) Union:A∪B . Intersection: An B. and Difference:A-B 3
Lectures 10 - 12 Intersection Problems, Nonlinear Solvers and Robustness Issues 10.1 Overview of intersection problems Intersections are fundamental in computational geometry, geometric modeling and design, analysis and manufacturing applications. Examples of intersection problems include: • Shape interrogation (eg. visualization) through contouring (intersection with a series of parallel planes, coaxial cylinders, and cones etc.) • Numerical control machining (milling) involving intersection of offset surfaces with a series of parallel planes, to create machining paths for ball (spherical) cutters. • Representation of complex geometries in the “Boundary Representation” scheme; for example, the description of the internal geometry and of structural members of cars, airplanes, ships, etc, involves – Intersections of free-form parametric surfaces with low order algebraic surfaces (planes, quadrics, torii). – Intersections of low order algebraic surfaces. in a process called boundary evaluation, in which the Boundary Representation is created by “evaluating” the Constructive Solid Geometry model of the object. During this process, intersections of the surfaces of primitives (see Figure 10.1) must be found during Boolean operations. Boolean opertations on point sets A, B include (see Figure 10.2) • Union: A ∪ B, • Intersection: A ∩ B, and • Difference: A − B. 3
Figure 10.1: Primitive solids →→→ Figure 10.2: Example of a Boolean operation: union All such operations involve intersections of surfaces to surfaces. In order to solve general surface to surface intersection problems, the following auxiliary intersection problems(similar to distance computation problems used in CAM for inspection of manufactured objects) need to be considered point/point(P/P) point /curve(P/C) point/surface(P/s) curve/ curve(C/c) curve/surface(C/S) All above types of intersection problems are also useful in geometric modeling, robotics collision avoidance. manufacturing simulation scientific visualization etc When the geometric elements involved in intersections are nonlinear(curved), intersection problems typically reduce to solving complex systems of nonlinear equations, which may be either polynomial or more general in character. Solution of nonlinear systems is a very com plex process in general in numerical analysis and there are specialized textbooks on the topic However, geometric modeling applications pose severe robustness, accuracy, automation, and efficiency requirements on solvers of a nonlinear systems as we will see later. Therefore, geo- metric modeling researchers have developed specialized solvers to address these requirements explicitly using geometric formulations
Figure 10.1: Primitive solids. ⇒ ⇒ ⇒ Figure 10.2: Example of a Boolean operation: union. All such operations involve intersections of surfaces to surfaces. In order to solve general surface to surface intersection problems, the following auxiliary intersection problems (similar to distance computation problems used in CAM for inspection of manufactured objects) need to be considered • point/point (P/P) • point/curve (P/C) • point/surface (P/S) • curve/curve (C/C) • curve/surface (C/S) All above types of intersection problems are also useful in geometric modeling, robotics, collision avoidance, manufacturing simulation, scientific visualization, etc. When the geometric elements involved in intersections are nonlinear (curved), intersection problems typically reduce to solving complex systems of nonlinear equations, which may be either polynomial or more general in character. Solution of nonlinear systems is a very complex process in general in numerical analysis and there are specialized textbooks on the topic. However, geometric modeling applications pose severe robustness, accuracy, automation, and efficiency requirements on solvers of a nonlinear systems as we will see later. Therefore, geometric modeling researchers have developed specialized solvers to address these requirements explicitly using geometric formulations. 4
10.2 Intersection problem classification Intersection problems can be classified according to the dimension of the problems and ac- cording to the type of geometric equations involved in defining the various geometric elements (points, curves and surfaces). The solution of intersection problems can also vary according to the number system in which the input is expressed and the solution algorithm is implemented 10.2.1 Classification by dimension /P, P/C, P/S /C, C/S S/S 10.2.2 Classification by type of geometry Parametric Implicit (eg. R=R(t) (eg.f(x,y)=0.,z=0) Rational Procedural Polynomial (algebraic) Figure 10.3: Curve geometry classification 1. Points plicit: R= Ro; R y, Procedural: Intersection of two procedural curves, procedural curve and surface, or three procedure surfaces, eg. offset or blending surfaces Algebraic: f(R)=g(R)=h(r)=0; where f, 9, and h are polynomials A classification of curves is illustrated in Figure 10.3 ● aranetrio R=R(t)A≤t≤B (a) Rational Polynomials(eg: NURBS, rational Bezier) (b)Procedural, eg: offsets, evolutes, ie. the locus of the centers of curvature of Implicit: These require solution of (usually nonlinear)equations
10.2 Intersection problem classification Intersection problems can be classified according to the dimension of the problems and according to the type of geometric equations involved in defining the various geometric elements (points, curves and surfaces). The solution of intersection problems can also vary according to the number system in which the input is expressed and the solution algorithm is implemented. 10.2.1 Classification by dimension • P/P, P/C, P/S • C/C, C/S • S/S 10.2.2 Classification by type of geometry polynomial Procedural Polynomial (algebraic) (eg. f(x, y) = 0, z = 0) Implicit Rational Parametric (eg. R = R(t) ) Figure 10.3: Curve geometry classification 1. Points • Explicit: R = R0; R = [x, y, z] • Procedural: Intersection of two procedural curves, procedural curve and surface, or three procedure surfaces, eg. offset or blending surfaces. • Algebraic: f(R) = g(R) = h(R) = 0; where f, g, and h are polynomials. 2. Curves A classification of curves is illustrated in Figure 10.3. • Parametric: R = R(t) A ≤ t ≤ B (a) Rational Polynomials (eg: NURBS, rational B´ezier). (b) Procedural, eg: offsets, evolutes, ie. the locus of the centers of curvature of a curve. • Implicit: These require solution of (usually nonlinear) equations 5
(a) Algebraic(polynomial) f(R)=g(R)=0 space curve 2=0, f(,y)=0 planar curves (b)Procedural offsets(eg. non-constant distance offsets involving convolution, see Pottmann 1997) 3. Surfaces Parametric R= R(u, v) where u, v vary in some finite domain, the parametric space (a)(Rational) Polynomial(eg: NURBS, Bezier, rational Bezier etc. (b)Procedural offsets blends generalIzed cylinders Implicit: Algebraic f(R)=0 where f is a polynomial 10.2.3 Classification by number system In our discussion of intersection problems, we will refer to various classes of numbers · integer numbers rational numbers, m/n, n+0, where m, n are integers floating point numbers in a computer(which are a subset of the rational numbers) radicals of rational numbers, eg. Vm/n,n+0, where m, n are integers algebraic numbers (roots of polynomials with integer coefficients transcendental (e, T, trigonometric, etc. ) numbers ●re imber interval numbers,(a, b], where a, b are real numbers; rounded interval numbers, [c, d], where c, d are floating point numbers Issues relating to floating point and interval numbers affecting the robustness of intersection Igorithms are addressed in the next section on nonlinear solvers
(a) Algebraics (polynomial) f(R) = g(R) = 0 space curves z = 0, f(x, y) = 0 planar curves (b) Procedural offsets (eg. non-constant distance offsets involving convolution, see Pottmann 1997) 3. Surfaces • Parametric R = R(u, v) where u, v vary in some finite domain, the parametric space. (a) (Rational) Polynomial (eg: NURBS, B´ezier, rational B´ezier etc.) (b) Procedural – offsets – blends – generalized cylinders • Implicit: Algebraics f(R) = 0 where f is a polynomial. 10.2.3 Classification by number system In our discussion of intersection problems, we will refer to various classes of numbers: • integer numbers; • rational numbers, m/n, n 6= 0, where m, n are integers; • floating point numbers in a computer (which are a subset of the rational numbers); • radicals of rational numbers, eg. q m/n, n 6= 0, where m, n are integers; • algebraic numbers (roots of polynomials with integer coefficients); • transcendental (e, π, trigonometric, etc.) numbers. • real numbers; • interval numbers, [a, b], where a, b are real numbers; • rounded interval numbers, [c, d], where c, d are floating point numbers. Issues relating to floating point and interval numbers affecting the robustness of intersection algorithms are addressed in the next section on nonlinear solvers. 6
10.3 Point/point"intersection Check if R1-R2 e, where e represents the maximum allowable tolerance Choice of"tolerances"in a geometric modeller is difficult-an open question Lack of transitivity, see Figure 10.4 R1=R2→R1≠R3 R=Rs R R Figure 10.4: Intersections of points within a tolerance is intransitive ● What should e reflect?
10.3 Point/point “intersection” • Check if |R1 − R2| < , where represents the maximum allowable tolerance. • Choice of “tolerances” in a geometric modeller is difficult–an open question. • Lack of transitivity, see Figure 10.4: R1 = R2 R2 = R3 ⇒ R1 6= R3 R2 R3 . . . R1 Figure 10.4: Intersections of points within a tolerance is intransitive. • What should reflect? 7
10.4 Point/curve intersection 10.4.1 Point/Implicit curve intersection Ro∩{z=0,f(x,y)=0} where f(a, y) is usually a polynomial(and f(a, y)=0 represents an algebraic curve). In an exact arithmetic context, we can substitute Ro in a, f(a, y)=0 and verify if the results are zero. Similarly, we could handle Ro∩{f(R)=9(R)=0} where f(r)=g(r=0 represents an implicit 3D space curve What does verify mean in"Hoating point"arithmetic? · Example A et zo=0 and o, yo satisfy f(xo,y0)<∈<1 where e is a small constant and If(a, y) s l in the domain of interest including(ao, yo) then a"distance"check can be performed by f(xo,3∥d 1 10.2 provided|f(xo,o川≠0. Equation10.1 is called the“ algebraic distance” and Equa tion 10.2 is called the non-algebraic distance". The true geometric distance is given d= min R-Rol; where R=(, y), f(r)=0 (10.3) The true geometric distance is difficult and expensive to compute(particularly for implicit f(r =0 and involves computing the global minimum of R-Rol. Equation 10.2 results from the first order approximation of Equation 10.3 as derived by Taylor expansion and is exact when f(R)is represents a plane
10.4 Point/curve intersection 10.4.1 Point/Implicit curve intersection R0 ∩ {z = 0, f(x, y) = 0} where f(x, y) is usually a polynomial (and f(x, y) = 0 represents an algebraic curve). In an exact arithmetic context, we can substitute R0 in {z, f(x, y) = 0} and verify if the results are zero. Similarly, we could handle: R0 ∩ {f(R) = g(R) = 0} where f(R) = g(R) = 0 represents an implicit 3D space curve. What does verify mean in “floating point” arithmetic? • Example A: Let z0 = 0 and x0, y0 satisfy |f(x0, y0)| < 1 (10.1) where is a small constant and |f(x, y)| ≤ 1 in the domain of interest including (x0, y0), then a “distance” check can be performed by: |f(x0, y0)| | 5 f(x0, y0)| < δ 1 (10.2) provided | 5 f(x0, y0)| 6= 0. Equation 10.1 is called the “algebraic distance” and Equation 10.2 is called the “non-algebraic distance”. The true geometric distance is given by: d = min|R − R0|; where R = (x, y), f(R) = 0 (10.3) The true geometric distance is difficult and expensive to compute (particularly for implicit f(R) = 0 and involves computing the global minimum of |R−R0|. Equation 10.2 results from the first order approximation of Equation 10.3 as derived by Taylor expansion and is exact when f(R) is represents a plane. 8
f(R)=0 R q(R)=0 Figure 10.5: Curves meet at small angle {f(R)=9(R)=0} When curves f=0,g=0 meet at a small angle (h tg= 1), then the condition d61 19l<e and 82=fo<8 (where If1, 191, 81, 82 are evaluated with R= Ro and 6, 8<1) are not enough to guarantee proximity of Ro to the intersection of f, g, see Figure 10.5 vf 83 vg Figure 10.6: Approximate curves with straight lines Using a linear approximation, and letting vf v9 V升| be the angle of intersection as in Figure 10.6 near the intersection point, a better criterion for evaluating if Ro is near the intersection of f and 63=o <6<1
. R0 g(R) = 0 f(R) = 0 Figure 10.5: Curves meet at small angle. • Example B: R0 ∩ {f(R) = g(R) = 0} When curves f = 0, g = 0 meet at a small angle ( 5f |5f| · 5g |5g| ∼= 1), then the condition |f| < and δ1 = |f| |5f| < δ |g| < and δ2 = |g| |5g| < δ (where |f|, |g|, δ1, δ2 are evaluated with R = R0 and , δ 1) are not enough to guarantee proximity of R0 to the intersection of f, g, see Figure 10.5. δ δ δ 1 2 3 φ g f Figure 10.6: Approximate curves with straight lines. Using a linear approximation, and letting φ = cos−1 | 5f | 5 f| · 5g | 5 g| | be the angle of intersection as in Figure 10.6 near the intersection point, a better criterion for evaluating if R0 is near the intersection of f and g is δ3 = φ −1 { |f| | 5 f| + |g| | 5 g| } < δ 1 9
10.4.2 Point /Parametric curve intersection 1. Rational polynomial curves Ro∩R=R(t)A≤t≤B Brute force elementary method We solve each of the following three nonlinear polynomial equations separately and we search for common real roots in a<t< B 0 t1,…,t (t)-30=0→t1,…,tn In principle, this elementary approach is "easy "for polynomials. However, in prad tice, this process is complex and inefficient and prone to numerical inaccuracies Preprocessing and subdivision method Use bounding box of R(t) to eliminate easily resolvable cases, with some level of subdivision(splitting) to reduce bo Concept of subdivision in rational arithmetic: To eliminate numerical error in the subdivision process(which can lead to erroneous decisions), rational arithmetic may be employed (if the input coefficients of R(t) are rational or Hoating point numbers). This can be easily done in object-oriented languages uch as C++ using operator overloading Continue subdivision until box is small -Then, we could use a numerical technique, such as F {R-R()2}t∈D1c[A,B nd use some t from the interval D1 as the initial approximation. Use of the square of the distance function is necessary to avoid possible divergence of the derivative of the distance function, if it approaches zero If the minimization process converges to to and VF(to)<8, t= to is the desired solution Implicitization(perhaps with box preprocessing) such as (o, yo)n f=r0),y=y(t)y Let us consider an example where a(t), y(t) are quadratic polynomials(the curve is a parabola). We will attempt to eliminate t to get a polynomial f(r, y)=0 which the a, y coordinates of all points on the curve satisfy. We start with the system =aot +bot +co aot+ bot t4+bot+co aot+ bot+co
10.4.2 Point/Parametric curve intersection 1. Rational polynomial curves R0 ∩ R = R(t) A ≤ t ≤ B • Brute force elementary method: We solve each of the following three nonlinear polynomial equations separately and we search for common real roots in A ≤ t ≤ B. x(t) − x0 = 0 → t 0 1 , · · · ,t 0 n y(t) − y0 = 0 → t 00 1 , · · · ,t 00 n z(t) − z0 = 0 → t 000 1 , · · · ,t 000 n In principle, this elementary approach is “easy” for polynomials. However, in practice, this process is complex and inefficient and prone to numerical inaccuracies. • Preprocessing and subdivision method – Use bounding box of R(t) to eliminate easily resolvable cases, with some level of subdivision (splitting) to reduce box size. – Concept of subdivision in rational arithmetic: To eliminate numerical error in the subdivision process (which can lead to erroneous decisions), rational arithmetic may be employed (if the input coefficients of R(t) are rational or floating point numbers). This can be easily done in object-oriented languages such as C++ using operator overloading. – Continue subdivision until box is small. – Then, we could use a numerical technique , such as: F(t) = min{|R0 − R(t)| 2 } t ∈ D1 ⊂ [A, B] and use some t from the interval D1 as the initial approximation. Use of the square of the distance function is necessary to avoid possible divergence of the derivative of the distance function, if it approaches zero. – If the minimization process converges to t0 and q F(t0) < δ, t = t0 is the desired solution. • Implicitization (perhaps with box preprocessing) such as (x0, y0) ∩ {x = x(t), y = y(t)} Let us consider an example where x(t), y(t) are quadratic polynomials (the curve is a parabola). We will attempt to eliminate t to get a polynomial f(x, y) = 0 which the x, y coordinates of all points on the curve satisfy. We start with the system x = a0t 2 + b0t + c0 y = a 0 0 t 2 + b 0 0 t + c 0 0 ⇒ a0t 2 + b0t + c0 − x = 0 a 0 0 t 2 + b 0 0 t + c 0 0 − y = 0 10