13.472J/1.128J/2.158J/16.940J COMPUTATIONAL GEOMETRY Lecture 20 Dr. K. H. Ko Prof. N. M. Patrikalakis Copyright c≤2003 Massachusetts Institute of Technology Contents 20 Advanced topics in di?erential geometry 2 20.1 Geodesics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 20.1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 20.1.2 De?nition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 20.1.3 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 20.1.4 Two-point boundary value problem . . . . . . . . . . . . . . . . . . . . . . . . . 5 20.1.5 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 20.2 Developable surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 20.2.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 20.2.2 De?nition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 20.2.3 Developable surface in terms of B′ezier surface . . . . . . . . . . . . . . . . . . . 12 20.2.4 Development of developable surface (?attening) . . . . . . . . . . . . . . . . . . 13 20.3 Umbilics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 20.3.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 20.3.2 De?nition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 20.3.3 Computation of umbilical points . . . . . . . . . . . . . . . . . . . . . . . . . . 15 20.3.4 Classi?cation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 20.3.5 Characteristic lines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 20.4 Parabolic, ridge and sub-parabolic points . . . . . . . . . . . . . . . . . . . . . . . . . 21 20.4.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 20.4.2 Focal surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 20.4.3 Parabolic points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 20.4.4 Ridge points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 20.4.5 Sub-parabolic points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Bibliography 25 Reading in the Textbook ? Geodesics : Chapter 10, pp.265 - 291 ? Umbilics : Chapter 9, pp.231 - 264 1 Lecture 20 Advanced topics in di?erential geometry 20.1 Geodesics In this section we study the computation of shortest path between two points on free-form surfaces [14, 11]. 20.1.1 Motivation ? ship design ? robot motion planning ? terrain navigation installation of underwater cable ? 20.1.2 De?nition ? t: Unit tangent vector of C at P n: Unit normal vector of C at P? N: Unit surface normal vector of S at P? ? u: Unit vector perpendicular to t in the tangent plane de?ned by N × t. 2 t N u n k k g k n S P C Figure 20.1: De?nition of geodesic curvature. ? We can decompose the curvature vector k of C into N component k n , which is called normal curvature vector, and u component k g , which is called geodesic curvature vector k = k n + k g = ?ρ n N + ρ g u (20.1) Here ρ n and ρ g are the normal and geodesic curvatures, respectively and de?ned as follows: ρ n = ?k N (20.2)· ρ g = k u (20.3)· ? Consequently, dt ρ g = ds · (N × t) (20.4) ? Geodesic paths are sometimes de?ned as shortest path between points on a surface, however this is not always a satisfactory de?nition. De?nition: Geodesics are curves of zero geodesic curvature [24]. 20.1.3 Governing equations ? The unit tangent vector of the curve C on the surface r is given by dr(u(s), v(s)) du dv t = = r u + r v (20.5) ds ds ds ? Hence using chain rules ? ? 2 ? ? 2 dt du du dv dv d 2 u d 2 v = r uu + 2r uv + r vv + r u + r v (20.6) ds ds ds ds ds ds 2 ds 2 3 ? Consider that the surface normal N has the direction of normal of the geodesic line ±n n r u = 0, n r v = 0. (20.7)· · dt ? Since kn = ds , equation (20.7) can be rewritten as dt dt r u = 0, r v = 0 (20.8) ds · ds · ? By substituting equation (20.6) into equations (20.8) we obtain ? ? 2 ? ? 2 du du dv dv d 2 u d 2 v (r uu · r u ) + 2(r uv · r u ) + (r vv · r u ) + E + F = 0 (20.9) ds ds ds ? ? 2 du du dv ds ds 2 ds 2 ? ? 2 dv d 2 u d 2 v (r uu · r v ) + 2(r uv · r v ) + (r vv · r v ) + F + G = 0 (20.10) ds ds ds ds ds 2 ds 2 v By eliminating d 2 from equation (20.9) using equation (20.10), and eliminating d 2 u from ds 2 ds 2 ? equation (20.10) using equation (20.9) and employing the Christo?el symbols, we obtain ? ? ? 2 d 2 u du ? 2 du dv dv + ? 1 + 2? 1 + ? 1 = 0 (20.11) ds 2 11 22 ds ? d 2 v du ? 2 12 ds ds ds ? ? 2 du dv dv + ? 2 + 2? 2 + ? 2 = 0 (20.12) ds 2 11 22 ds 12 ds ds ds ? jk (i, j, k = 1, 2) are the Christo?el symbols de?ned as follows: Where ? i ? 1 GE u ? 2F F u + F E v ? 2 2EF u ? EE v + F E u = = 11 11 ? 2(EG ? F 2 ) , 2(EG ? F 2 ) 1 GE v ? F G u ? 2 EG u ? F E v = = 12 12 ? 2(EG ? F 2 ) , 2(EG ? F 2 ) (20.13) 1 2GF v ? GG u + F G v ? 2 EG v ? 2F F v + F G u = = 22 22 2(EG ? F 2 ) , 2(EG ? F 2 ) ? These two second order di?erential equations can be rewritten as a system of four ?rst order di?erential equations [6]. du = p (20.14) ds dv = q (20.15) ds dp 2 2 = ?? 1 ? 2? 1 11 p 12 pq ? ? 1 (20.16) 22 q ds dq 2 2 = ?? 2 ? 2? 2 11 p 12 pq ? ? 2 (20.17) 22 q ds 4 ? Euler Lagrange Equations (Calculus of Variations) We can also ?nd this result by means of the general rules of the calculus of variations. We want to minimize ? b ? b I = ds = f(u, v, v˙)du (20.18) a a where dv 2 f(u, v, v˙) = ? E + 2Fv˙ + Gv˙ , v˙ = (20.19) du This leads to the condition πf d πf (20.20) πv ? du πv˙ = 0 from which we can derive the same di?erential equation for geodesics. 20.1.4 Two-point boundary value problem ? We can solve a system of four ?rst order ordinary di?erential equations (20.14) to (20.17) as – Initial-value problem (IVP), where all four boundary conditions are given at one point, or as – Boundary-value problem (BVP), where four boundary conditions are speci?ed at two distinct points. ? The ?rst order di?erential equation for a boundary value problem can be written in vector form as: dy = g(s, y), y(A) = (u A , v A , p A , q A ) T , y(B) = (u B , v B , p B , q B ) T (20.21) ds where p A , p B , q A and q B are unknowns, y = (u, v, p, q) T (20.22) 2 12 pq ? ? 1 2 2 22 q 2 ) T g = (p, q, ?? 1 ? 2? 1 22 q , ?? 2 ? 2? 2 12 pq ? ? 2 (20.23) 11 p 11 p ? There are two commonly used approaches to the numerical solution of BVP. 1. Shooting method: easy to implement but unstable. 2. Relaxation method: more sophisticated but stable. ? Shooting method – We assume values at s = A, which are not given as boundary conditions at s = A and compute the solution of the resulting IVP to s = B. 5 ? ? ? (?) (?) – The computed values of y(B) will not, in general, agree with the corresponding boundary condition at s = B. – Consequently, we need to adjust the initial values and try again. – The process is repeated until the computed values at the ?nal point agree with the boundary conditions and referred as shooting method. – Formulation: Using the ?rst fundamental form, given p A we can obtain q A from 2 ?Fp A ± F 2 p A ? G(Ep 2 q A ? 1) A = . G Thus we assume p A and solve the di?erential equation as IVP using, say Runge- Kutta method. Here we also have to assume the entire arc length of the geodesic path s to stop the integration. Thus the unknowns can be considered as p A and s. ? ), the di?erence can be given B ? , v B If we denote the computed value of (u B , v B ) as (u as (u This can be done by employing the Newton’s method T? )? v B B ? ? u , v B B . We need to adjust p A and s to make the di?erence zero. ?1 ? ? ? B ?s ? B ?u?u ? ? ? ? ? ? ? u B p A p A u ?p A B ?= ? B ?s ? B ?v?v s s ? v B v Bn+1 ?p A n n where ? B ? ( )p A B ? ( + ? ) ?p p u A A B ? B ? ( )s B ? ( + ? ) ?s s u B πu πuu u = , = πP A ?p A πs ?s ? B v ? ( )p A B ? ( + ? ) ?p p v A A B ? B ? ( )s B ? ( + ? ) ?s s v B πv πv v = , = πP A ?p A πs ?s Relaxation method ? dy – The second method is based on a ?nite di?erence approximation to ds on a mesh of points in the interval [A, B]. – This method starts with an initial guess and improves the solution iteratively and referred as, direct method, relaxation method or ?nite di?erence method. – The shooting method is often very sensitive to the unknown initial angles at point A and unless a good initial guess is provided, the integrated path will never reach the other point B, while the relaxation method starts with two end points ?xed and relaxes to the true solution and hence it is much more stable. – Let us consider a mesh of points satisfying A = s 1 < s 2 < . . . < s m = B. We approximate the n ?rst order di?erential equations by the trapezoidal rule [8]. 1Y k ? Y k?1 = [G k + G k?1 ], k = 2, 3, . . ., m (20.24) s k ? s k?1 2 6 a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0 where the n-vectors Y k , G k are meant to approximate y(s k ) and g(s k ) with bound- ary conditions Y 1 = β = (u A , v A , p 1 , q 1 ) T , Y m = ? = (u B , v B , p m , q m ) T (20.25) – Y 1 has n 1 = 2 known components, while Y m has n 2 = n ? n 1 = 4 ? 2 = 2 known components. u2u1 u3 u4 uM?2 uM?1 uM a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a32 a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a32 Boundary Conditions Unknowns a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0 v1 v2 v3 v4 vM?2 vM?1 vM a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a32 a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0p1 p2 p3 p4 pM?2 pM?1 pM a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a32 q1 q2 q3 q4 . . . . . . . . . qM?2 qM?1 qM a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0 s1 s2 s3 s4 sM?2 sM?1 sM Figure 20.2: Mesh points. s – This discrete approximation will be accurate to the order of h 2 (h = max k {s k ? k?1 }). Equation (20.24) forms a system of (m? 1)n nonlinear algebraic equations with (m ? 1)n unknowns. – Let us refer to equation (20.24) as 1 F k = (F 1,k , F 2,k , . . ., F n,k ) T = Y k ? Y k?1 [G k + G k?1 ] = 0 (20.26) s k ? s k?1 ? 2 – and the boundary conditions (20.25) as F F 1 = (F 1,1 , F 2,1 , . . ., F n 1 ,1 ) T = Y 1 ? β = 0, (20.27) m+1 = (F 1,m+1 , F 2,m+1 , . . ., F n 2 ,m+1 ) T = Y m ? ? = 0 – then we have mn nonlinear algebraic equations F = (F T 2 , . . .,F T 1 ,F T m+1 ) T = 0 (20.28) 7 Points Tolerance Iterations Geodesic Distance m κ N L M R L M R 101 1.0E-3 10 1 10 1.661 1.865 1.661 Table 20.1: Numerical conditions and results for the computation of the geodesic path between corner points of the wave-like surface. – and can be computed by quadratically convergent Newton iteration, if a su?ciently accurate starting vector Y (0) = (Y T 2 , . . . , Y T 1 , Y T m ) T is provided. The Newton iter- ation scheme is given by Y (i+1) [J = Y (i) + ?Y (i) (20.29) (i) ]?Y (i) = ?F (i) (20.30) where [J (i) ] is the mn by mn Jacobian matrix of F (i) with respect to Y (i) . x y z Figure 20.3: Geodesic paths on the wave-like bicubic B-spline surface between points of two corners. 20.1.5 Example Bilinear surface r(u, v) = (u, v, uv). E E E = 1 + v 2 , F = uv, G = 1 + u 2 u = 0, F u = v, G u = 2u v = 2v, F v = u, G v = 0 8 x y z Figure 20.4: Convergence of the right geodesic path in Figure 20.3. ? 1 ?2(uv) · v + uv · (2v) = = 0 11 2 2[(1 + v 2 )(1 + u 2 ) ? u v 2 ] ? 2 = ? 1 22 = 0 22 = ? 2 ? 11 1 v = 12 ? u 2 + v 2 + 1 2 u = 12 u 2 + v 2 + 1 Finally the di?erential equations are given by du = p ds dv = q ds dp ?2u = ds u 2 + v 2 + 1 pq dq ?2v = ds u 2 + v 2 + 1 pq. 9 20.2 Developable surface Developable surfaces are a special class of surfaces that can be developed or unfolded onto a plane without stretching or tearing [1, 9, 5, 16, 15]. 20.2.1 Motivation ? In shipbuilding, doubly curved plates are manufactured by roller and line heating pro- cesses, while the singly curved plates (developable surface) are manufactured by roller only. ? The use of developable surfaces has several advantages such as lower labor costs in construction, smaller capital investment in equipment, ease of repair and simple tools for construction. ? In automobile production, body panels, upholstery and window glass are developable surfaces. 20.2.2 De?nition ? A ruled surface is de?ned as a surface generated by the motion of a straight line referred as a generator or ruling [24]. ? The mathematical representation of a ruled surface is given by r(u, v) = r A (u) + vβ(u) (20.31) where r(u) is a 3D curve called the directrix or base curve of the ruled surface and β(u) is a unit vector which gives the direction of the ruling at each point on the directrix see Figure 20.5. ruling or generator directrix α(u) r(u,v)= + rA(u) rA(u) vα(u) Figure 20.5: De?nition of ruled surface 10 ? An alternate expression based on rulings joining corresponding points on two space curves r A (u) and r B (u) is given by r(u, v) = (1 ? v)r A (u) + vr B (u) (20.32) Bilinear surface is a ruled surface. ? r r A (u) = (1 ? u)b 00 + ub 10 B (u) = (1 ? u)b 01 + ub 11 r(u, v) = (1 ? v)r A (u) + vr B (u) = (1 ? u)(1 ? v)b 00 (u) + (1 ? v)ub 10 + v(1 ? u)b 01 + uvb 11 ? A ruled surface is a developable surface if and only if [7] ˙ β) = 0 (20.33)r A · (β × ˙ where × and · are the cross and dot products or (r A ? r B ) · (r˙ A × r˙ B ) = 0 (20.34) ? The following statements are the equivalent necessary and su?cient conditions for a surface to be developable [19]. 1. Gaussian curvature is zero. 2. Geodesics on a developable surface can be mapped onto straight lines in the plane, and the straight lines in plane can map into geodesics on a developable surface. 3. The normal vectors on a developable surface along the ruling are parallel. 4. Developable surfaces possess the same tangent plane at all points of the same gen- erator. ? If the direction of the ruling β(u) is constant, the condition for developability is auto- matically satis?ed since β(u) = 0. This implies that the ruled surface is a cylinder.˙ ? If the direction of the ruling β(u) is given by r A (u) ? a, then the condition (20.33) becomes r˙ A (u) · ((r A ? a) × r˙ A (u)) → 0 and hence the surface is developable. In this case the surface is a cone with apex a. ? Finally, if β(u) is a tangent vector to r A (u), then again the condition (20.33) is satis?ed, and the resulting developable surface is called a tangential developable surface. In this case the base curve coincides with the so-called the edge of regression or cuspidal edge. 11 20.2.3 Developable surface in terms of B′ezier surface ? Let us de?ne a developable B′ezier surface in terms of two B′ezier curves (directrices) and rulings between pairs of points from each curve [1], see Figures 20.6. ? We restrict the two directrices (r A (t), r B (t)) to lie in parallel planes so that r˙ B (t) = ζ(t)r˙ A (t). Here ζ(t) denotes a linear function of t. ? Then the condition for developability becomes (ζ(t) ? 1)r˙ A × (r B ? r A ) · r˙ A 0 and hence→ it is automatically satis?ed. ? Given the control points of the design curve r A (t), we want to determine those of r B (t). Example: r A (t) quadratic and r B (t) cubic a 1 a 2 b 0 b 1 b 2 b 3 generator directrix r B directrix r A a 0 Figure 20.6: A B′ezier Developable Surface with Quadratic r A (t) and Cubic r B (t) Two B′ezier Curves? 2 r A (t) = s a 0 + 2sta 1 + t 2 a 2 2 r B (t) = s 3 b 0 + 3s tb 1 + 3st 2 b 2 + t 3 b 3 where s = 1 ? t ? Tangent vectors r˙ A (t) = 2[s(a 1 ? a 0 ) + t(a 2 ? a 1 )] 2 r˙ B (t) = 3[s 2 (b 1 ? b 0 ) + 2st(b 2 ? b 1 ) + t (b 3 ? b 2 )] Scalar function ? ζ(t) = ζ 0 (1 ? t) + ζ 1 t = ζ 0 s + ζ 1 t 12 ? Substitute r˙ A (t), r˙ B (t) and ζ(t) into r˙ B (t) = ζ(t)r˙ A (t) and collect terms and equating the coe?cients of the independent functions s 2 , 2st and t 2 to zero. ? We get three equations 2ζ 0 (a 1 ? a 0 ) = 3(b 1 ? b 0 ) ζ 0 (a 2 ? a 1 ) + ζ 1 (a 1 ? a 0 ) = 3(b 2 ? b 1 ) 2ζ 1 (a 2 ? a 1 ) = 3(b 3 ? b 2 ) ? 6 scalar equations with 10 scalar unknowns (b 0 , b 1 , b 2 , b 3 , ζ 0 , ζ 1 ) ? We can ?x b 0 and b 3 . 20.2.4 Development of developable surface (?attening) ? In the manufacture of developable surfaces, it is necessary to calculate the plane devel- opment of these surfaces [7]. ? The development is based on the isometric mapping [13]. – The length of any arc on the developable surface is the same as that of on the developed plane. – The coe?cients of the ?rst fundamental forms at the corresponding points are the same. – Isometric surfaces have the same Gaussian curvature at corresponding points. Cor- responding curves on those surfaces have the same geodesic curvature at correspond- ing points. – Every isometric mapping is conformal (the angle of intersection of every arbitrary pair of intersecting arcs are preserved). Procedures ? – Map one of the directrix isometrically onto the plane by integrating a system of ordinary di?erential equation from A to C, see Figure 20.7. d 2 x dy g (s) = 0 (20.35) ds 2 ? ρ ds d 2 y dx + ρ g (s) = 0 (20.36) ds 2 ds – Compute the angle of two isoparametric lines at A. r u (0, 0) r v (0, 0) cos β = (20.37) r u (0, 0) · r v (0, 0)| | | | 13 Generator B C D θ Directrix Directrix Generator A Figure 20.7: Developed surface. – Since we know the length and direction of the generator at A, we may obtain the point B. – Integrate a system of ordinary di?erential equation from B to D. – Connect C and D. 14 ? 20.3 Umbilics 20.3.1 Motivation ? Identi?cation of singular points in the principal direction ?eld ? Fingerprints of shapes ? Object matching ? Object recognition 20.3.2 De?nition An umbilic is a point on a surface where the normal curvatures in all directions are equal and the principal curvature directions are indeterminate. Locally, a surface around an umbilical point can be best approximated by a circle whose radius is equal to a radius of curvature at the umbilical point. 20.3.3 Computation of umbilical points For NURBS surfaces, umbilical points can be calculated by solving a non-linear system of equations derived from the de?nition using the Gaussian (K) and the mean (H) curvatures as follows [24]: ρ 1,2 (u, v) = H(u, v) ± H 2 (u, v) ? K(u, v). (20.38) Let W (u, v) = H 2 ? K. The principal curvatures, ρ 1,2 , are real valued functions so that W √ 0 must hold. From the de?nition of the umbilical point we have W (u, v) = 0. With these two conditions combined, we can infer that at an umbilical point, W (u, v) has a global minimum [17, 18]. Here, we assume that W is at least C 2 smooth. Then, the condition that W has a global minimum at an umbilic implies that ≡W = 0. Therefore, at an umbilic the following equations hold [18]: πW (u, v) πW (u, v) W (u, v) = 0, = 0, = 0. (20.39) πu πv Given a polynomial parametric surface patch such as a rational B′ezier surface patch, we can set W = P N , where P N and P D are polynomials in u and v. With the condition W √ 0, P P D N √ 0 is assured since P D > 0 is always true under the regularity condition of the surface [24]. The equation W = 0 is equivalent to P N = 0. The ?rst derivative of W is ?W = ?x i ( ?P N ?P D D (i = 1, 2), where x 1 = u and x 2 = v, which is reduced to ?W = ( ?P N )/P D ?x i P D ? P N ?x i )/P 2 ?x i ?x i using P N = 0. Therefore, equations (20.39) are reduced to [18] πP N πP N P N (u, v) = 0, = 0, = 0. (20.40) πu πv To locate isolated umbilical points, a set of equations (20.40) can be solved by the rounded interval projected polyhedron algorithm [23, 21]. But when the IPP algorithm encounters re- gions of umbilical points, it slows down dramatically. A di?erent approach, called the adaptive 15 quadtree decomposition which uses the quadtree decomposition and the convex hull properties of Bernstein polynomials, can be adopted in such a case [12]. 20.3.4 Classi?cation Umbilical points can be isolated or form lines or regions. They are classi?ed into two types based on their stability with respect to small perturbations: ? generic ? non-generic Generic umbilical points are stable with respect to small perturbations. Isolated generic umbilical points are further categorized into three types as shown in Figure 20.8 [2]: ? lemon ? star ? monstar Star type umbilical points are further classi?ed into the hyperbolic star and the elliptical star type umbilical points. Several methods are available for the classi?cation of isolated generic umbilical points. What follows is an introduction to techniques of umbilical point classi?cation. Lemon Star Monstar Figure 20.8: Three generic umbilics Index The type of isolated generic umbilical points can be determined by the index. The index is the amount of rotation that a straight line segment tangent to the lines of curvature experiences 16 when rotating in the counterclockwise direction along a small closed path around an umbilic [18]. The equation for direct calculation of the index is given as follows [18]: 1 n ? Index = 2? i=0 ?? i (20.41) ? ? ? ? 1 1 L+?E M +?F where ?? i = ? (i+1) mod n ?? i , ? ? ? ?? i ? ? and ? i = tan ?1 ? , or tan ?1 ? . 2 2 M +?F N +?G Figure 20.9 shows the principal direction ?elds around two distinct types of umbilical points. The index can distinguish the star type umbilical point from the monstar or lemon type umbil- Figure 20.9: Umbilics with principal direction ?elds: 1 ical point. If the index is ? 2 , then the umbilical point is of the star type, whereas the umbilical 1 point is of the monstar or lemon type if the index is 2 . Being topological, this distinction is very robust [10]. Complex θ plane method An umbilical diagram shown in Figure 20.10 [22] is a comprehensive and easy way to distinguish the type of an isolated generic umbilical point. In order to use this diagram, the local surface near an umbilical point has to be represented as a height function or the Monge form with respect to a local coordinate system as follows [18]: r = (x, y, h(x, y)). (20.42) The height function h(x, y) is Taylor expanded at the origin of the local coordinate system. Then we have h(x, y) = ρ (x 2 + y 2 ), (20.43)? 2 + 1 (ax 3 + 3bx 2 y + 3cxy 2 + dy 3 ) + O(4), 6 where ρ is the normal curvature at the umbilical point. Let us set 2 C(x, y) = ax 3 + 3bx 2 y + 3cxy + dy 3 , (20.44) and denote C(x, y) as cubic form. This expression implies that the local structure of a surface near an umbilical point is dominated by the coe?cients of C(x, y), i.e. by (a, b, c, d), which 17 determine the type of the umbilical point [20, 22]. It is convenient to represent the cubic part C(x, y) in the complex plane for analysis purposes. If we set ? = x + iy, then the cubic part C(x, y) becomes ? 3 C(?) = ω? 3 + 3?? 2 ? + 3??? 2 + ω? , (20.45) with 1 ω = 8 [(a ? 3c) + i(d ? 3b)], (20.46) 1 ? = [(a + c) + i(b + d)], 8 where ω = 0. We can express (20.45) in a coordinate system rotated about the normal vector ≥ without losing any essential features to make the coe?cient of ? 3 be equal to 1 [22]. Using ? = ω 1 3 ?, the equation (20.45) becomes 3 C(? ? 3 2 ? ? 2 ) = + 3θ? + 3θ? ? + ? , (20.47) where θ = ?ω ? 1 3 ω ? 2 3 . This means that the cubic part C(x, y) is parametrized with respect to a single complex variable θ [4, 22]. Therefore, the variations of C(x, y) can be mapped onto the complex plane [4, 20, 22]. When ω = 0, the equation (20.45) is reduced to ? C(?) = 3??(?? + ??). (20.48) This equation corresponds to the in?nity in the θ plane [22, 4, 20], which is not considered in this discussion. ? ? Depending on the structure of C(x, y) (in turn C(?)), three characteristic lines are deter- mined as follows [22, 20]: 1 : β ? 1 3 (2e i? + e ?2i? ),? θ = 1, ? | | ? ? 2 : β ? (2e i? + e ?2i? ), where ? 1 and ? 2 are maps from β to the complex θ-plane. They divide the complex plane into sub-regions as shown in Figure 20.10. Each sub-region corresponds to a speci?c type of an umbilical point. In Figure 20.10, ES means the elliptic star, HS the hyperbolic star, MS the monstar and L the lemon. If θ falls on a dividing curve, then the corresponding umbilical point is of non-generic type. The behavior of such an umbilical point can be analyzed with more higher order terms [18]. Using this diagram, the type of an umbilical point is easily determined, see [22, 4, 20]. 20.3.5 Characteristic lines : ? ? 1 3 (2e i? + e ?2i? )? 1 The cubic form (20.47) is parabolic on the deltoid ? 1 [22]. This implies that there are three root lines of (20.47) but two of them coincide. Inside the deltoid, the cubic form (20.47) is 18 ? ? ? 1 Γ ESHS HS HS MS MS MS L L L |ω|=1 2 Γ Figure 20.10: The umbilic diagram adapted from [22] elliptic [22]. It is hyperbolic outside the deltoid. This classi?cation is directly related to the number of ridge lines passing through an umbilical point, and the existence of extrema of the principal curvatures near the umbilical point [4, 10, 18]. Here, a ridge point is de?ned as an extremum point of a line of curvature, and a ridge line is a set of such points [10, 20]. Inside the deltoid, the number of ridge curves is three, the extrema of principal curvatures exist, and an umbilical point is of the elliptical star type [4, 10]. On the other hand, outside the deltoid only one ridge curve passes through an umbilical point, no extremum of a principal curvature exists, and the umbilical point is of the hyperbolic star type [4, 10]. |?| = 1 On the circle θ = 1, the cubic form (20.47) is right-angled [22]. When the root directions of | | the Hessian (20.49) are mutually orthogonal, we call that the cubic form (20.44) is right-angled [22]. Here, the Hessian is de?ned as ? ? 2 C ? 2 C ? ?x?y ? ? H e (x, y) = 1 det ? ? ?x 2 ? 2 C ? . (20.49) 6 ? ? 2 C ? ?x?y ?y 2 This implies that the maximum and minimum lines of curvature are orthogonal at an umbilical point to form approximately a plain rectangular grid pattern [22]. This circle is related to the index [10]. Inside the circle, the index is ? 1 2 , and an umbilical point is of the star type. 19 Outside the circle, the index is 1 2 , and the umbilical point is classi?ed as the lemon or monstar type [10, 18]. On the circle, the index value is zero, and the umbilical point is of non-generic type. This circle also has a relation with the birth/death of generic umbilical points under the evolution of a surface [22]. i? + e ?2i? )? 2 : ? ? (2e Another cubic form, called the Jacobian cubic form is de?ned to explain the deltoid ? 2 as follows: 2 2 3 U (x, y) = bx 3 ? (2c ? a)x y ? (2b ? d)xy ? cy , (20.50) whose root lines are tangent to the lines of curvature near an umbilical point [10]. On the deltoid ? 2 , the cubic form (20.50) becomes parabolic [22]. The Jacobian cubic form (20.50) is related to the number of extrema of the cubic form (20.44) [22]. The cubic form (20.44) can be represented as C(r, β) in polar coordinates with x = r cos β and y = r sin β, and the expression of the direction in which the local extrema of the cubic form C(r, β) occur, i.e. dC(?) = 0 [18] d? is reduced to the Jacobian cubic form (20.50). Inside the deltoid ? 2 , there are three real root lines of the Jacobian cubic form (20.50) or three directions of the extrema of the cubic form (20.44), which implies that three lines of curvature converge to an umbilical point [10, 22]. This umbilical point is classi?ed as the star or monstar type [10, 22]. Outside the deltoid ? 2 , there is one root line of the Jacobian cubic form (20.50), and no extremum of the cubic form (20.44) exists. An umbilical point of this case is of the lemon type [10, 22]. 20 20.4 Parabolic, ridge and sub-parabolic points 20.4.1 Motivation ? Sophisticated classi?cation of umbilics ? Stronger matching conditions than the umbilic matching ? Registration ? Analysis of surface evolution 20.4.2 Focal surfaces For any point on a surface r, two points f 1 and f 2 are called focal points which are de?ned as follows: 1 1 f 1 = r + N, f 2 = r + N (20.51) ρ 1 ρ 2 where N is a normal vector at P and ρ 1 and ρ 2 are the maximum and minimum principal curvatures, respectively. The focal surface is a set of focal points which can be considered as the envelope of the normal to the surface. The focal surface will play an important role in studying ridges, crests and sub-parabolic lines. Figure 20.11: Focal Surfaces 21 20.4.3 Parabolic points Contact function Consider a surface r in Monge form as follows: 2 z = h(x, y) = a 0 x 2 + a 1 xy + a 2 y + H.O.T, (20.52) so that the tangent plane at the origin is the x,y-plane. The parabolic points can be identi?ed 2 by studying the intersection between r and the tangent plane. If a 1 = 4a 0 a 2 then the surface has A 2 contact with the tangent plane. De?nition The parabolic points on a regular surface are points where the tangent plane has specially high ’contact’ with the surface, i.e. A 2 contact [3]. Alternatively, they can be interpreted as the points which separate the elliptic and hyperbolic regions of the surface, namely, points or curves where the Gaussian curvature is zero. 20.4.4 Ridge points Contact functions The sphere of curvature is a sphere centered at one of the centers of principal curvature and having a radius equal to the corresponding radius of curvature. Consider a sphere centered at (0, 0, r) and passing through the origin such that it is tangent to the surface r there. Then the equation of the sphere is 2 2 2 x + y + (z ? r) 2 = r (20.53) The contact function is de?ned with substitution of z = h(x, y). 2 g(x, y) = x 2 + y + (h(x, y) ?r) 2 ?r 2 = 0 (20.54) Expanding g(x, y) as a power series in x and y gives r 2 2 2 g(x, y) = x 2 (1 ? rρ 1 ) + y (1 ? rρ 2 ) ? (b 0 x 3 + 3b 1 x y + 3b 2 xy + b 3 y 3 ) 3 r (c 0 x 4 ρ 2 y 2 ) 2 ? 12 + ···) + ( ρ 1 x 2 + + (20.55) 2 2 ··· De?nition ? The ridge point of a smooth surface is a point where one of the sphere of curvature has more degenerate contact than the usual A 2 contact and the curve of intersection becomes ρ 1 g(x, y) = y 2 (ρ 1 ? ρ 2 ) ? 1 (3b 1 x 2 y + 3b 2 xy 2 + b 3 y 3 ) 3 ρ 3 1 (c 0 x 4 ? 12 + ···) + 1 x 4 + 4 ··· 22 b ? ? 2 1 b 2 b 3 2 x y= (ρ 1 ? ρ 2 ) y ? 2(ρ 1 ? ρ 2 ) 2 ? 2(ρ 1 ? ρ 2 ) xy ? 6(ρ 1 ? ρ 2 ) 1 1 + (ρ 1 ? ρ 2 )(c 0 ? 3ρ 3 1 ))x 4 + (20.56)? 12(ρ 1 ? ρ 2 ) (3b 2 ··· 1 – If r = 1 (or r = ? 2 ) and b 0 = 0 (or b 3 = 0), then the contact function (20.55) ? 1 becomes A 3 singularity. – The coe?cient P 1 = 3b 2 1 ) should not be zero. 1 + (ρ 1 ? ρ 2 )(c 0 ? 3ρ 3 ? A point is a ridge point relative to a principal direction p if and only if the principal curvature ρ 1 corresponding to p has an extremum along the line of curvature in the direction p. ? Ridge points are the pre-image of cuspidal edges on a focal surface of a smooth surface. 20.4.5 Sub-parabolic points ? The point on which one principal curvature has an extremum relative to motion along the other line of curvature is called the sub-parabolic point on the surface. ? It is the pre-image of a parabolic point on the focal surface. 23 Bibliography [1] G. Aumann. Interpolation with developable B′ezier patches. Computer Aided Geometric Design, 8(5):409–420, November 1991. [2] M. V. Berry and J. H. Hannay. Umbilic points on Gaussian random surfaces. Journal of Physics A., 10(11):1809–1821, 1977. [3] J. W. Bruce, P. J. Giblin, and F. Tari. Parabolic curves of evolving surfaces. International Journal of Computer Vision, 17(3):291–306, 1996. [4] J. W. Bruce, P. J. Giblin, and F. Tari. Ridges, crests and sub-parabolic lines of evolving surfaces. International Journal of Computer Vision, 18(3):195–210, 1996. [5] J. S. Chalfant and T. Maekawa. Design for manufacturing using B-spline developable surfaces. Journal of Ship Research, 42(3):207–215, 1998. A. Bj¨[6] G. Dahlquist and ? orck. Numerical Methods. Prentice-Hall, Inc., Englewood Cli?s, NJ, 1974. [7] I. D. Faux and M. J. Pratt. Computational Geometry for Design and Manufacture. Ellis Horwood, Chichester, England, 1981. [8] J. H. Ferziger. Numerical Methods for Engineering Applications. Wiley, 1981. [9] W. H. Frey and D. Bindschadler. Computer-aided design of a class of developable B′ezier surfaces. R&D Publication 8057, General Motors, September 1993. [10] P. W. (editor) Hallinan. Two and Three Dimensional Patterns of the Face. A.K. Peters, Ltd., Wellesley, Massachusetts, 1999. [11] R. Kimmel, A. Amir, and A. M. Bruckstein. Finding shortest paths on surfaces using level sets propagation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 17(6):635–640, June 1995. [12] K. H. Ko, T. Maekawa, N. M. Patrikalakis, H. Masuda, and F.-E. Wolter. Shape intrinsic ?ngerprints for free-form object matching. In G. Elber and V. Shapiro, editors, 8th ACM Symposium on Solid Modeling and Applications, Seattle, WA, June 2003, In press. [13] E. Kreyszig. Di?erential Geometry. University of Toronto Press, Toronto, 1959. [14] T. Maekawa. Computation of shortest paths on free-form parametric surfaces. Journal of Mechanical Design, Transactions of the ASME, 118(4):499–508, December 1996. [15] T. Maekawa and J. S. Chalfant. Computation of in?ection lines and geodesics on devel- opable surfaces. Mathematical Engineering in Industry, 7(2):251–267, 1998. [16] T. Maekawa and J. S. Chalfant. Design and tessellation of B-spline developable surfaces. Journal of Mechanical Design, Transactions of the ASME, 120(3):453–461, September 1998. [17] T. Maekawa and N. M. Patrikalakis. Interrogation of di?erential geometry properties for design and manufacture. The Visual Computer, 10(4):216–237, March 1994. [18] T. Maekawa, F.-E. Wolter, and N. M. Patrikalakis. Umbilics and lines of curvature for shape interrogation. Computer Aided Geometric Design, 13(2):133–161, March 1996. [19] M. J. Mancewicz and W. H. Frey. Developable surfaces: Properties, representations and methods of design. R&D Publication 7637, General Motors, 1992. [20] R. Morris. The sub-parabolic lines of a surface. In G. Mullineux, editor, The Mathematics of Surfaces VI, Proceedings of the 6th IMA Conference on Mathematics of Surfaces VI, pages 79–102, Oxford, UK, 1996. Clarendon Press. [21] N. M. Patrikalakis and T. Maekawa. Shape Interrogation for Computer Aided Design and Manufacturing. Springer-Verlag, Heidelberg, February 2002. [22] I. R. Porteous. Geometric Di?erentiation for the intelligence of curves and surfaces. Cambridge University Press, Cambridge, 1994. [23] E. C. Sherbrooke and N. M. Patrikalakis. Computation of the solutions of nonlinear polynomial systems. Computer Aided Geometric Design, 10(5):379–405, October 1993. [24] D. J. Struik. Lectures on Classical Di?erential Geometry. Addison-Wesley, Cambridge, MA, 1950.