13.472J/1.128J/2.158J/16.940J COMPUTATIONAL GEOMETRY Lecture 8 N. M. Patrikalakis Massachusetts Institute of Technology Cambridge, MA 02139-4307, USA cCopyright ≥2003 Massachusetts Institute of Technology Contents 8 Fitting, Fairing and Generalized Cylinders 2 8.1 Least Squares Method of Curve Fitting . . . . . . . . . . . . . . . . . . . . . . 2 8.2 Fairing of Curves and Surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 8.2.1 Properties and De?nition . . . . . . . . . . . . . . . . . . . . . . . . . . 4 8.2.2 Curve Interrogation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 8.2.3 Fairing Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 8.2.4 Surface Fairing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 8.3 Generalized Cylinders: Motivation and De?nitions . . . . . . . . . . . . . . . . 12 8.3.1 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 8.3.2 De?nition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 8.4 Degeneracies of Generalized Cylinders . . . . . . . . . . . . . . . . . . . . . . . 16 8.5 Properties of Generalized Cylinders . . . . . . . . . . . . . . . . . . . . . . . . 20 8.6 Discrete Generalized Cylinders . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Bibliography 22 Reading in the Textbook ? Chapter 11, pp. 353 - 365 1 ? ? Lecture 8 Fitting, Fairing and Generalized Cylinders 8.1 Least Squares Method of Curve Fitting Example problem Given N points P i , i = 1, 2, ..., N (N ? 4), construct an approximating cubic B′ezier curve that interpolates P 1 and P N (end points). Solution 1. Parametrization by chord-length method Let ?u 1 = 0; ? ?u i+1 = u i + d i+1 , i = 1, 2, ..., N ? 1 (8.1) where d i+1 = |P i+1 ? P i | is the chord length between two consecutive points. The overall chord length is N d = d i (8.2) i=2 The parametric value associated with point P i ?u i = u i /d (8.3) which is normalized as u i ? [0, 1] with u 1 = 0 and u N = 1. 2. Linear equations A cubic B′ezier curve is de?ned as 3 Q(u) = Q i B i,3 (u), 0 ? u ? 1 (8.4) i=0 where B i,3 (u) are the cubic Bernstein polynomials. 2 ? Obviously, the boundary conditions require Q 0 = P 1 , Q 3 = P N . The problem is then represented as a linear system with N ? 2 equations and 2 unknowns: 2 Q i B i,3 (u j ) = P j ? P 1 B 0,3 (u j ) ? P N B 3,3 (u j ) i=1 = L j , j = 2, 3, ..., N ? 1 (8.5) B or in matrix form (N ?2)×2 · q 2×1 = l (N ?2)×1 (8.6) 3. Least Squares Method De?ne the mean square error as E 2 = |B · q ? l| 2 (8.7) then E 2 = (B · q ? l) T (B · q ? l) = q T B T Bq ? 2q T B T l + l T l (8.8) is a function of q and is minimized if we set θE 2 = 0 ≤ B T Bq ? B T l = 0 (8.9) θq ≤ B T Bq = B T l (normal equations) (8.10) ≤ q = (B T B) ?1 B T l (formal solution) (8.11) The extension to ?tting with B-splines is similarly formulated. Notes: 1. The choice of internal knots in the B-spline basis should re?ect any knowledge of deriva- tive discontinuity in the data, as shown in Figure 8.1. 2. Greater density of knots is needed in rapidly changing parts of the shape. 3. NAG routines for approximate ?tting of cubic B-splines [9] (a) Curves: E02BAF (b) Surfaces: E02DAF & E02ZAF 4. NAG routines on least square problems provide more ?exibility. 3 Double knot for cubics Figure 8.1: Set of data re?ecting a possible discontinuity of tangent vector. 8.2 Fairing of Curves and Surfaces 8.2.1 Properties and De?nition Motivation: 1. Spline curves resulting from (a) interpolation of points; (b) manipulation of polygon, usually need fairing. 2. Screen plots ( small resolution ) are misleading concerning curve quality. 3. Full scale plots. 4. Curvature plots are useful as they allow isolation of problem areas on raster devices. Properties of fair curves: [3, 4] 1. Curvature continuity ( C 2 ). 2. Curvature is almost piecewise linear with as few spans as possible. For cubics with simple knots, property (1) is automatically satis?ed. If R(t) is reasonably parametrized, |R ? (t)| is constant, and the curvature ρ(s) = |R ? (t) × R ?? (t)| |R ? (t)| 3 ∈ |R ?? (t)| (8.12) Property (2) thus requires that |R ?? (t)| be almost piecewise linear. That means that R ??? (t) needs to be constant, which leads to the following de?nition of fairness. Q De?nition: Q, R are two C 2 cubic splines in t ? [a, b]. Q is fairer than R if for r ? [a, b], ??? + ) ? Q ??? + ) ?R ??? [ (r (r ? )] 2 ? [R ??? (r (r ? )] 2 (8.13) If r is a knot at which interpolation of data occurs, the above expression means that reducing “shear” forces from supports increases the fairness of splines, if we consider geometric splines as approximations of physical splines. 4 ? ?? 8.2.2 Curve Interrogation We interrogate the fairness of a spline curve by looking into the planar projection of their curvature. The signed curvature is de?ned as x ?? y ? ? x y ρ(u) = (8.14) 3 (x ? 2 + y ? 2 ) 2 and it is easy to point out the changes of sign ( in?ections ). The curve is fair if a plot of ρ(u), which is made up of a few monotone segments, is continuous. The aim of the fairing process is to locate the places of maximum discontinuity of ρ ? (u) and fair those places (see ?gure (8.2)). κ(u) u Fair here Figure 8.2: Plot of curvature ρ(u) along curve as a function of u. 8.2.3 Fairing Methods Assume that we have a spline curve obtained by interpolation of measured data. The curve generally has unwanted behavior. Curve fairing eliminates imperfections by changing data within a measurement tolerance. 1. Kjellander’s Method [7] Procedure (a) Obtain data points P j (t j ), 1 ? j ? N . (b) Fit the points with a spline R(t). + j )?R ??? (t ? (c) Find the knot where max|R ??? (t j )| occurs and attempt fairing at P J where j = J corresponds to the worst jump. (d) Using Hermite or B′ezier, construct cubic interpolation Q(t), which interpolates R(t J ?1 ) = P J ?1 (8.15) R(t J +1 ) = P J +1 (8.16) R ? (t J ?1 ) , R ? (t J +1 ) (8.17) 5 ? ? ? ? ? Pj?1 Pj Pj+1 R(t) point to be faired Rnew(tj) Figure 8.3: Kjellander’s fairing method (e) Determine new curve ? R(t) if t ?? [t J ?1 , t J +1 ] (8.18)R new (t) = Q(t) if t ? [t J ?1 , t J +1 ] (f) Compute R new (t J ). Notice that R new (t) is in?nitely di?erentiable there. (g) Construct a spline curve on P 1 ,··· ,P J ?1 ,R new (t J ),P J +1 ,··· ,P N . (h) The resulting curve is usually fairer at t J . Disadvantages ? Global scheme. ? Repeated interpolation. 2. Farin’s Method [4] Recall Boehm’s knot insertion method n n+1 ? ? P i N i,4 (t) = P i N i,4 (t) (8.19) i=0 i=0 ? [t 0 , t 1 ,···] ≤ t 0 ,··· , t l , t, t l+1 ,··· (8.20) and the control points are ? P i = τ i P i + (1 ? τ i )P i?1 (8.21) where ? ? 1 0 ? i ? l ? 3 ? τ i = 0 l + 1 ? i ? n + 1 (8.22) ? ? t ? ?t i ? l ? 2 ? i ? l t i+3 ?t i Hence ? P i = P i 0 ? i ? l ? 3 ? P i = P i?1 l + 1 ? i ? n + 1 ? P i = τ i P i + (1 ? τ i )P i?1 l ? 2 ? i ? l (8.23) 6 ? Idea of Farin’s method (a) Remove a knot ?rst to make curve in?nitely di?erentiable at that location such that the curve is fairer now in that area. (b) Insert the removed knot back so as to have the same knot vector (if needed; not usually necessary). ? Knot removal Knot removal is an inverse process of knot insertion. From Equations (8.23), we have ? P i = P i , i = 0, 1,··· , l ? 3 ? τ i P i + (1 ? τ i )P i?1 = P i , i = l ? 2, l ? 1, l ? P i = P i+1 , i = l, l + 1,··· , n (8.24) Here we have n + 2 equations and n + 1 unknowns, therefore, an approximate solution can be obtained by the least squares method. A sample solution using Farin’s method: ? P i = P i , i = 0, 1,··· , l ? 3 ? P i = P i+1 , i = l, l + 1,··· , n (8.25) and from the least squares method, we have the following equations for P l?2 ,P l?1 ? ? ? ? ? P l?3 τ l?2 0 ? ? P l?2 ? (1 ? τ l?2 ) ? ? ? ? τ l?1 ? P l?2 = ? ? ? 1 ? τ l?1 ? ? (8.26) ? P l?1 ? 0 1 ? τ l P l?1 ? ? P l ? τ l P l+1 or in the matrix form A · p = f ≤ A T Ap = A T f (8.27) which yields p = (A T A) ?1 A T f (8.28) This should be followed by knot insertion to complete the fairing process (if necessary). ? Knot insertion ? T = [t 0 ,··· , t l , t, t l+1 ,···] = [T 0 ,··· , T l , T l+1 , T l+2 ,···] (8.29) where the removed knot t ? is inserted as the knot T l+1 . Hence the control points of the faired curve can be determined by the knot insertion method ( Equation 8.23), where t ? ? t l T l+1 ? T l t τ l = = l+3 ? t l T l+4 ? T l t ? ? t l?1 T l+1 ? T l?1 t τ l?1 = = l+2 ? t l?1 T l+3 ? T l?1 t ? ? t l?2 T l+1 ? T l?2 t τ l?2 = = (8.30) l+1 ? t l?2 T l+2 ? T l?2 ? ? ? If we remove all or many of the knots as in Figure 8.5, the other constraints (such as deviation) will dominate the problem. Let the curve obtained by either interpolation or approximation be n R(t) = R i N i,4 (t) (8.31) i=0 and the curve after knot removal be m ? ? R 0 (t) = Q i N i,4 (t) (8.32) i=0 where m < n. Then, (n? m) knots removed before are inserted to make R 0 (t) have the same knot vector as R(t), such that n R 0 (t) = Q i N i,4 (t) (8.33) i=0 where Q i , 0 ? i ? n, are new control points determined by knot insertion, and N i,4 (t) are B-spline basis function over the same knot vector as that of R(t). The deviation of two B-spline curves n ? |R(t) ? R 0 (t)| = | (R i ? Q i )N i,4 (t)| (8.34) i=0 n ? ? max i |R i ? Q i | i=0 N i,4 (t) (8.35) = max i |R i ? Q i | (8.36) remove all the knots Figure 8.5: Deviation of fair curve. 8.2.4 Surface Fairing ? Di?erential Geometry Review Let ρ 1 , ρ 2 be principal curvatures, then Gaussian curvature G = ρ 1 ρ 2 (8.37) 1 mean curvature H = 2 (ρ 1 + ρ 2 ) (8.38) absolute curvature K = |ρ 1 | + |ρ 2 | (8.39) Local Taylor expansion is 1 (ρ 1 x 2 z = + ρ 2 y 2 ) + h.o.t. (8.40) 2 1. Elliptic case: G = ρ 1 ρ 2 > 0, see Figure 8.6(a). - convex surface - one side of tangent plane - no in?ection 2. Hyperbolic case: G = ρ 1 ρ 2 < 0, see Figure 8.6(b). - non-convex surface (locally) - intersection with tangent plane - surface in?ection Z κ1 κ2 Y Z Y X X Figure 8.6: (a) Elliptic case; (b) Hyperbolic case. De?nition 1 Surface in?ection at a point P exists i? the surface crosses the tangent plane at P. If G < 0, surface in?ection at P. If G = 0 in a region ( e.g., a “cylinder” ), then surface in?ection exists in the region if H changes sign in the region ( not one point). De?nition 2 Surface curve in?ection exists at a point P of a surface i? a planar surface curve through P changes sign of its curvature at P. Surface curve in?ection ≤ Surface in?ection ? Surface Interrogation 1. Gaussian curvature is not useful in “cylinders” where G = 0. 2. Mean curvature is zero for minimal surfaces where ρ 1 = ρ 2 . 3. Absolute curvature does not have such problems. Its curvature plots are useful in study- ing surface imperfections ( e.g., small oscillations ). 4. Re?ection lines for parallel line light source ( on raster screen, mark points where normal passes through source ) are useful for global imperfection detection. Surface fairing may be performed by fairing auxiliary curves de?ned for each of the columns and rows of the control polyhedron. Generalized Cylinders 8.3 Generalized Cylinders: Motivation and De?nitions Generalized cylinders [12, 13, 15, 8] or sweeps provide greater generality for shape representa- tion than tensor product surfaces. A generalized cylinder is a representation of an elongated object viewed as having a main axis (directrix or spine) and a smoothly varying cross section (generatrix), as de?ned in Figure 8.7. Both directrix and generatrix can be open or closed (periodic) curves. 8.3.1 Applications Some of the more common applications of generalized cylinders are ? Representation of measured data (e.g. CAT scans, deformed solids). ? Representation of manufacturing processes. ? Representation of blends. ? Object recognition and scene interpretation in robotics and computer vision. ? Representation of human and animal shapes. Generatrix Directrix (Spine) Figure 8.7: Generalized cylinder. 8.3.2 De?nition ? Given: 1. A bounded 3-D curve serving as spine. 2. A cross-sectional plane swept along the spine perpendicular to it so that the spine passes through the origin of the 2-D coordinate system on the plane. 3. A cross-sectional curve on the cross-section plane de?ned locally in the cross-section coordinate system, where the size and shape of the curve may vary with the param- eter along the spine curve. ? The surface swept by the curve is a generalized cylinder. ? Some examples of generalized cylinders and the resulting surfaces are: 1. Spine = straight line, generatrix = circle =≤ CYLINDER 2. Spine = circle, generatrix = circle =≤ TORUS 3. Spine = straight line segment, generatrix = linearly tapering circles =≤ CONE Mathematical description of generalized cylinders: ? Directrix (spine): A = A(s), 0 ? s ? 1. ? Generatrix: C = C(t; s) = [x(t; s), y(t; s), 0], 0 ? t ? 1. ? Generalized cylinder surface patch: R(s, t) = A(s) + x(t; s)X(s) + y(t; s)Y(s), where X, Y, Z are orthogonal 3-D unit vectors with the Z tangent to A(s), i.e. A ? (s) Z(s) = |A ? (s)| As an example, X(s), Y(s) could be chosen equal to the normal and binormal vectors of the spine curve A(s) or by rotation of those by some angle, see Figure 8.8. Problems with generalized cylinder representation When A(s) is straight line, X(s), Y(s), Z(s) should be de?ned independent of the Frenet trihedron, eg. using X(s), Y(s) as constants. ? C(t; s) ? k ? Y (s) O ? i ? j ? Z(s) ? X(s) ? A(s) t Figure 8.8: Components of a generalized cylinder. Examples: a. Generalized Cylinders with B-Spline Spine and Generatrix Curves k ? A(s) = A i N i,K (s) i=0 l q ? ? C(t; s) = C ij N i,L (t)N j,Q (s) i=0 j=0 = [x(t; s), y(t; s), 0] where N i,K (s), N i,L (t) and N j,Q (s) are B-spline basis functions or the Bernstein special- izations. b. Pipe Surfaces (See also Section 11.6 of textbook [11]) When the generatrix is a circle, the resulting generalized cylinder is a pipe surface. The pipe surface P C (r) with radius r can be parametrized using the Frenet trihedron (t(t),n(t),b(t)) of the spine curve C(t) as follows: P(t, ?) = C(t) + r[cos ?n(t) + sin ?b(t)] (8.41) where t ? [0, 1] and ? ? [0, 2?]. 8.4 Degeneracies of Generalized Cylinders a. Local Self-Intersection b. Global Self-Intersection Figure 8.9: Types of self-intersection of generalized cylinders. ? X(s) ?n ? Y (s) ? t = ? Z(s) [x(t; s) (t; s), 0] π x i i plane cross-section , y straight line curvature center of , y Figure 8.10: Criterion to avoid local self-intersection of generalized cylinders. There are two types of the degeneracies of generalized cylinders illustrated in Figure 8.9, namely local self-intersection and global self-intersection. A condition to avoid local self- intersection of generalized cylinders is illustrated in Figure 8.10. The condition is 2 max t (x 2 + y ) ? π 2 (s) for all s, where π(s) is the radius of curvature of the spine. ? As a special case, we consider local self-intersection of pipe surfaces (see also Maekawa et al. (1998) for details). The partial derivative of the pipe surface with respect to t is given by P t (t, ?) = C ? (t) + r[cos ?n (t) + sin ?b ? (t)]. (8.42) Equation (8.42) can be rewritten using the Frenet formulae ( n ? (t) = |C ? (t)|(?ρ(t)t + δ(t)b), b ? (t) = ?|C ? (t)|δ(t)n ) as P t (t, ?) = |C ? (t)|(1 ? ρ(t)r cos ?)t ? r|C ? (t)| sin ?δ(t)n + r|C ? (t)| cos ?δ(t)b (8.43) P where ρ(t) and δ(t) are the curvature and torsion of the spine curve. Similarly we can derive ? as P ? (t, ?) = r[? sin ?n(t) + cos ?b(t)]. (8.44) The surface normal of the pipe surface can be obtained by taking the cross product of equations (8.43) and (8.44) yielding P t × P ? = ?|C ? (t)|r[1 ? ρ(t)r cos ?][sin ?b(t) + cos ?n(t)]. (8.45) It is easy to observe that the pipe surface becomes singular (the normal vector vanishes) when 1 ? ρ(t)r cos ? = 0. Since cos ? varies between -1 and 1, there will be no local self-intersection if ρ(t)r < 1. Therefore, to avoid local self-intersection we need to ?nd the largest curvature ρ max of the spine curve and set the radius of the pipe surface such that r < 1/ρ max . Figure 8.11 shows an example of local self-intersection. The curvature ρ(t) of a space curve C(t), is given by |C ? (t) × C ?? (t)| ρ(t) = . (8.46) |C ? (t)| 3 Thus, to ?nd the largest curvature ρ max we need to locate the critical points of ρ(t), i.e. solve the equation ρ ? (t) = 0, and decide whether they are local maxima. Then we compare these local maxima with the curvature at the end points, i.e. ρ(0) and ρ(1), and obtain the global largest curvature. This problem can be solved by elementary calculus. If the spine curve is given by a rational B′ezier curve, equation ρ ? (t) = 0 reduces to a single univariate nonlinear polynomial equation. In the case where the spine curve is a rational B-spline, we can extract the rational B′ezier segments by knot insertion. Cho et al. [1] describe in detail how to obtain ρ ? (t) = 0 for integral B′ezier curves. Global self-intersection of a pipe surface involves the following types of intersections: 1. End circle to end circle: Two end circles of the pipe surface touch each other, see Fig- ure 8.12. 2. Body to body: Two di?erent body portions of the pipe surface touch each other, see Figure 8.13. 3. End circle to body: One of the end circles touches the body, see Figure 8.14. -4 -2 0 2 4 0 2 4 6 8 10 12 14 16 -0.5 0 0.5 Figure 8.11: Local self-intersection The theory of intersections and nonlinear solvers is needed to handle these global intersec- tion problems and we will discuss these later. Let ρ max be the maximum curvature of the spine curve, and r ee , r bb , r eb be the maximum possible upper limit radius of the pipe surface such that it does not globally self-intersect between end circle to end circle, body to body and end circle to body of the pipe surface, respectively. Then we have Theorem Let p(r) be the pipe surface with spine curve c(t) and radius r. Then p(r) is nonsingular if and only if r < ? = min{1/ρ max , r ee , r bb , r eb }. 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 0.5 1 1.5 2 2.5 3 2 2.5 3 3.5 4 Figure 8.12: End circle to end circle global self-intersection -0.2 0 0.2 0.4 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 -0.6 -0.4 -0.2 0 0.2 0.4 Figure 8.13: Body to body tangential intersection and local self-intersection -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 -0.2 0 0.2 0.4 0.6 Figure 8.14: End circle to body tangential global and local self-intersections 8.5 Properties of Generalized Cylinders 1. Unit normal vector R s × R t ?n = |R s × R t | 2. Radius of curvature of spine– A(s) = [x, y, z] |A s | 3 = π(s) |A s × A ss | 2 2 (x 2 + y + z s ) s s = 3 2 [(x s y ss ? y s x ss ) 2 + (y s z ss ? z s y ss ) 2 + (z s x ss ? x s z ss ) 2 ] 1 2 8.6 Discrete Generalized Cylinders Useful in the interpretation of measured data. The constructive de?nition for discrete gener- alized cylinders is: 1. De?ne a piecewise continuous spine. 2. Obtain point measurements on cross-section curves on planes perpendicular to spine at a discrete set of points on spine. 3. Construct a local system of coordinates on each cross-section with origin on spine. 4. Interpolate each cross-section with splines and establish parametric correspondence be- tween cross-sections, see Figure 8.15. 5. Establish an interpolation rule between cross-sections, x(t; s i ), y(t; s i ), z(t; s i ), ?(s i ), see Figure 8.15. b 1 t n 1 1 t b n t b n 2 2 2 3 33 b 2 ’ n 2 ’ θ 2 θ 1 θ 3 = 0 ’ b 3 n ’ 3 Figure 8.15: Cross sections along the spine curve Bibliography [1] W. Cho, T. Maekawa, and N. M. Patrikalakis. Topologically reliable approximation of composite B′ezier curves. Computer Aided Geometric Design, 13(6):497–520, August 1996. [2] Q. Ding and B. J. Davies. Surface Engineering Geometry for Computer-Aided Design and Manufacture. Ellis Horwood, Chichester, UK, 1987. [3] G. Farin. Curves and Surfaces for Computer Aided Geometric Design: A Practical Guide. Academic Press, Boston, MA, 3rd edition, 1993. [4] G. Farin, G. Rein, N. Sapidis, and A. J. Worsey. Fairing cubic B-spline curves. Computer Aided Geometric Design, 4(1-2):91–103, July 1987. [5] I. D. Faux and M. J. Pratt. Computational Geometry for Design and Manufacture. Ellis Horwood, Chichester, England, 1981. [6] J. Hoschek and D. Lasser. Fundamentals of Computer Aided Geometric Design. A. K. Peters, Wellesley, MA, 1993. Translated by L. L. Schumaker. [7] J. A. Kjellander. Smoothing of cubic parametric splines. Computer-Aided Design, 15:175– 179, 1983. [8] T. Maekawa, N. M. Patrikalakis, T. Sakkalis, and G. Yu. Analysis and applications of pipe surfaces. Computer Aided Geometric Design, 15(5):437–458, May 1998. [9] Numerical Algorithms Group, Oxford, England. NAG Fortran Library Introductory Guide, Mark 18 edition, 2000. [10] J. Owen. STEP: An Introduction. Information Geometers, Winchester, UK, 1993. [11] N. M. Patrikalakis and T. Maekawa. Shape Interrogation for Computer Aided Design and Manufacturing. Springer-Verlag, Heidelberg, February 2002. [12] J. Pegna. Variable Sweep Geometric Modeling. PhD thesis, Stanford University, Stanford, CA, 1987. [13] J. Pegna and D. J. Wilde. Spherical and circular blending of functional surfaces. In Proceedings of the 7th International Conference on O?shore Mechanics and Arctic Engi- neering, volume 7, pages 73–82, Houston, Texas, February 1988. New York:ASME, 1988. [14] L. A. Piegl and W. Tiller. The NURBS Book. Springer, New York, 1995. [15] U. Shani and D. H. Ballard. Splines as embeddings for generalized cylinders. Computer Vision, Graphics and Image Processing, 27:129–156, 1984. [16] F. Yamaguchi. Curves and Surfaces in Computer Aided Geometric Design. Springer- Verlag, NY, 1988.