13472J/1.128J/2158J/16940J COMPUTATIONAL GEOMETRY Lecture 8 N.M. Patrikalakis Massachusetts Institute of Technology Cambridge MA 02139-4307 USA Copyright 2003 Massachusetts Institute of Technology Contents 8 Fitting, Fairing and Generalized Cylinders 8. 1 Least Squares Method of Curve Fitting 8.2 Fairing of Curves and Surfaces 8.2.1 Properties and Definition 8.2.2 Curve Interrogation 224455 8.2.3 Fairing Methods 8.2.4 Surface Fairin 8.3 Generalized Cylinders: Motivation and Definitions 12 8.3.1 Applications 8.3.2 Definition 8.4 Degeneracies of Generalized Cylinders 8.5 Properties of Generalized Cylinders 8.6 Discrete Generalized Cylinders Bibliography 22 Reading in the Textbook Chapter 11, pp. 353-365
13.472J/1.128J/2.158J/16.940J COMPUTATIONAL GEOMETRY Lecture 8 N. M. Patrikalakis Massachusetts Institute of Technology Cambridge, MA 02139-4307, USA Copyright ≥c 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 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . 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 Definitions . . . . . . . . . . . . . . . . 12 8.3.1 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 8.3.2 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 Given N points Pi, i= 1, 2,.,N(N> 4), construct an approximating cubic Bezier curve that interpolates PI and PN(end points Solution 1. Parametrization by chord-length method Let (8.1) where di+1=Pi+1-Pil is the chord length between two consecutive points. The overall chord length is d d (8.2) The parametric value associated with point Pi which is normalized as u; E[0, 1] with u=0 uN 2.L is defined Q(u)=∑QB13(u),0≤u≤1 where Bi3(u) are the cubic Bernstein polynomials
� � Lecture 8 Fitting, Fairing and Generalized Cylinders 8.1 Least Squares Method of Curve Fitting Example problem Given N points Pi, i = 1, 2, ..., N (N � 4), construct an approximating cubic B´ezier curve that interpolates P1 and PN (end points). Solution 1. Parametrization by chord-length method Let uˆ1 = 0; uˆ ˆ i+1 = ui + di+1, i = 1, 2, ..., N − 1 (8.1) where di+1 = |Pi+1 − Pi| is the chord length between two consecutive points. The overall chord length is N d = di (8.2) i=2 The parametric value associated with point Pi ui = uˆi/d (8.3) which is normalized as ui � [0, 1] with u1 = 0 and uN = 1. 2. Linear equations A cubic B´ezier curve is defined as 3 Q(u) = QiBi,3(u), 0 � u � 1 (8.4) i=0 where Bi,3(u) are the cubic Bernstein polynomials. 2
Obviously, the boundary conditions require Qo= P1, Q3= PN. The problem is then represented as a linear system with N-2 equations and 2 unknowns Q: Bi,3(ui)=P,-P1Bo,(ui)-PNB3,3(u i=1 2,3,…N-1 (8.5) or in matrix form B(N-2)×2·q2x1=1(x-2)×1 (8.6) 3. Least Squares Method Define the mean square error as (8.7) then E 1)(B·q-1) BB Bl (8.8) is a function of q and is minimized if we set B Bq-Bl=0 (8.9) *B Bq=B 1(normal equations) (8.10 + q=(BB)-B 1(formal solution) (8.11) The extension to fitting with B-splines is similarly formulated. 1. The choice of internal knots in the B-spline basis should reflect 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 fitting of cubic B-splines [9 (a)Curves: E02BAF (b)Surfaces: E02DAF E02ZAF 4. NAG routines on least square problems provide more flexibility
� Obviously, the boundary conditions require Q0 = P1, Q3 = PN . The problem is then represented as a linear system with N − 2 equations and 2 unknowns: 2 QiBi,3(uj) = Pj − P1B0,3(uj) − PN B3,3(uj) i=1 = Lj , j = 2, 3, ..., N − 1 (8.5) B or in matrix form (N−2)×2 · q2×1 = l(N−2)×1 (8.6) 3. Least Squares Method Define the mean square error as E2 = |B · q − l| 2 (8.7) then E2 = (B · q − l) T (B · q − l) = qT BT Bq − 2qT BT l + l T l (8.8) is a function of q and is minimized if we set θE2 = 0 ≤ BT Bq − BT l = 0 (8.9) θq ≤ BT Bq = BT l (normal equations) (8.10) ≤ q = (BTB) −1 BT l (formal solution) (8.11) The extension to fitting with B-splines is similarly formulated. Notes: 1. The choice of internal knots in the B-spline basis should reflect any knowledge of derivative 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 fitting of cubic B-splines [9] (a) Curves: E02BAF (b) Surfaces: E02DAF & E02ZAF 4. NAG routines on least square problems provide more flexibility. 3
Double knot for cubics Figure 8.1: Set of data reflecting a possible discontinuity of tangent vector 8.2 Fairing of Curves and Surfaces 8.2.1 Properties and Definition 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( C2) 2. Curvature is almost piecewise linear with as few spans as possible For cubics with simple knots, property(1)is automatically satisfied. If R(t)is reasonably parametrized, IR'(t)I is constant, and the curvature k(s)=R()×R"(l R'(t)3 CR(t) Property(2)thus requires that R (t)I be almost piecewise linear. That means that R"(t)needs to be constant, which leads to the following definition of fairness Definition: Q, R are two C cubic splines in t E a, b. Q is fairer than R if for r E a, Q"(r+)-Q(r-)2≤[R 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
Double knot for cubics Figure 8.1: Set of data reflecting a possible discontinuity of tangent vector. 8.2 Fairing of Curves and Surfaces 8.2.1 Properties and Definition 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 ( C2 ). 2. Curvature is almost piecewise linear with as few spans as possible. For cubics with simple knots, property (1) is automatically satisfied. 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 definition of fairness. Q Definition: Q, R are two C2 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 defined as (x2+y2) and it is easy to point out the changes of sign( inflections). The curve is fair if a plot of r(u), hich 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 n'(u)and fair those places(see figure(8.2)) k(u) Fair here Figure 8.2: Plot of curvature k(a) 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 P (a) Obtain data points F/(t),1≤j≤N (b) Fit the points with a spline R(t) (c)Find the knot where macR"(t#)-R"(t )l occurs and attempt fairing at P,where j=J corresponds to the worst jump (d) Using Hermite or Bezier, construct cubic interpolation Q(t), which interpolates R(t1-1)=P R P J+ (8.16) R'(t1-1),R(t+1) (8.17)
� �� 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 defined 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 ( inflections ). 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 figure (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 Pj(tj ), 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 PJ where j = J corresponds to the worst jump. (d) Using Hermite or B´ezier, construct cubic interpolation Q(t), which interpolates R(tJ−1) = PJ−1 (8.15) R(tJ+1) = PJ+1 (8.16) R� (tJ−1) , R� (tJ+1) (8.17) 5
R(t) point to be faired Figure 8.3: Kjellander's fairing method (e) Determine new curve Rnew(t) R(t) if tE[tJ-1, tj+1 Q(t) (8.18) (f)Compute Rnew(tJ). Notice that Rnew(t) is infinitely differentiable there (g) Construct a spline curve on P1,……,PJ-1,Rnen(t),PJ+1,…,PN (h)The resulting curve is usually fairer at tJ Disadvantages ● Global scheme Repeated interpolation 2. Farin's Method 4 Recall boehm's knot insertion method PN4t)=∑PN4() (8. tu, t, tu 20 and the control points are a2P+(1-a)P 21 0<i<l-3 0 l+1<i<n+1 Hence P;=P20≤i<1-3 P1=P-1l+1<i<n+1 P=aP+(1-a)P-1l-2≤i≤l
� � � � � 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 ⇒� [tJ−1, tJ+1] Rnew(t) = (8.18) Q(t) if t � [tJ−1, tJ+1] (f) Compute Rnew(tJ ). Notice that Rnew(t) is infinitely differentiable there. (g) Construct a spline curve on P1, · · ·, PJ−1, Rnew(tJ ), PJ+1, · · ·, PN . (h) The resulting curve is usually fairer at tJ . Disadvantages • Global scheme. • Repeated interpolation. 2. Farin’s Method [4] Recall Boehm’s knot insertion method n n+1 P ˆ ˆ iNi,4(t) = PiNi,4(t) (8.19) i=0 i=0 [t ˆ 0, t1, · · ·] ≤ t0, · · · , tl, t, tl+1, · · · (8.20) and the control points are Pˆi = τiPi + (1 − τi)Pi−1 (8.21) where � ⎧⎧ 1 0 � i � l − 3 τi = 0 l + 1 � i � n + 1 (8.22) ⎧⎧⎠ t ˆ−ti l − 2 � i � l ti+3−ti Hence Pˆi = Pi 0 � i � l − 3 Pˆi = Pi−1 l + 1 � i � n + 1 Pˆi = τiPi + (1 − τi)Pi−1 l − 2 � i � l (8.23) 6
. Idea of farin's method (a) remove a knot first to make curve infinitely differentiable 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=0.1.….l a2P+(1-a)P P -2,l-1 P=Px+1,t=l,l+1,……,n (8.24 Here we have n+2 equations and n+ l unknowns, therefore, an approximate solution can be obtained by the least squares method A sample solution using Farin's method P P 0,1,…,1-3 P l,1+1 825 and from the least squares method, we have the following equations for Pl-2, Pl-1 0 P P-2-(1-a1-2)P 1-al-1 P 26 P PI-agPl+ or in the matrix form which yields P=(AAAf This should be followed by knot insertion to complete the fairing process(if necessary) . Knot insertion T={to,…,t,t,t+1,…]=[1,…,T,T+1,T+2, 29 where the removed knot t is inserted as the knot Ti+1. Hence the control points of the aired curve can be determined by the knot insertion method( Equation 8.23), where t-t Ti+1-Th t1+3-tt Ti+4-T1 t1+2-t1-1T4+3-T7 t-tr (8.30)
• Idea of Farin’s method (a) Remove a knot first to make curve infinitely differentiable 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 = Pi, i = 0, 1, · · ·, l − 3 τ ˆ iPi + (1 − τi)Pi−1 = Pi, i = l − 2, l − 1, l P ˆ i = Pi+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 = Pi, i = 0, 1, · · ·, l − 3 P ˆ i = Pi+1, i = l, l + 1, · · · , n (8.25) and from the least squares method, we have the following equations for Pl−2, Pl−1 ⎨ ⎛ ⎨ ⎛ τl−2 0 ˆ Pl−3 ⎝ ⎞ Pl−2 − (1 − τl−2)ˆ ⎩ ⎩ ⎜ τl−1 ⎜ Pl−2 ⎪ 1 − τ = ⎩ ˆ l−1 � ⎪ P ⎜ (8.26) l−1 � 0 1 − τl Pl−1 Pˆ ˆ l − τlPl+1 or in the matrix form A · p = f ≤ ATAp = AT f (8.27) which yields p = (ATA) −1 AT f (8.28) This should be followed by knot insertion to complete the fairing process (if necessary). • Knot insertion T = [t ˆ 0, · · · , tl, t, tl+1, · · ·] = [T0, · · ·, Tl, Tl+1, Tl+2, · · ·] (8.29) where the removed knot t ˆ is inserted as the knot Tl+1. Hence the control points of the faired curve can be determined by the knot insertion method ( Equation 8.23), where t ˆ− tl Tl+1 − Tl t τl = = l+3 − tl Tl+4 − Tl t ˆ− tl−1 Tl+1 − Tl−1 t τl−1 = = l+2 − tl−1 Tl+3 − Tl−1 t ˆ− tl−2 Tl+1 − Tl−2 t τl−2 = = (8.30) l+1 − tl−2 Tl+2 − Tl−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 R()=∑RN4(t) and the curve after knot removal be Ro(t)=∑QN4(1) where m< n. Then, (n-m) knots removed before are inserted to make Ro(t) have the same knot vector as R(t), such that Ro(t)=∑QN4(t) where Q ,0<is n, are new control points determined by knot insertion, and Ni. 4 (t) are B-spline basis function over the same knot vector as that of r(t) The deviation of two B-spline curves R(t)-Ro(t)l=I2(R-Q:)Ni,4(t) max Imax R -Q remove all the knots Figure 8.5: Deviation of fair curve
� � � 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) = RiNi,4(t) (8.31) i=0 and the curve after knot removal be m R ˜ ˜ 0(t) = QiNi,4(t) (8.32) i=0 where m < n. Then, (n − m) knots removed before are inserted to make R0(t) have the same knot vector as R(t), such that n R0(t) = QiNi,4(t) (8.33) i=0 where Qi, 0 � i � n, are new control points determined by knot insertion, and Ni,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) − R0(t)| = | (Ri − Qi)Ni,4(t)| (8.34) i=0 �n � max i |Ri − Qi| i=0 Ni,4(t) (8.35) = max i |Ri − Qi| (8.36) remove all the knots Figure 8.5: Deviation of fair curve
8.2.4 Surface Fairing Differential Geometry Re Let K1, K2 be principal curvatures, then Gaussian curvature G= K1K2 mean curvature H =-(K1+K2) absolute curvature K =[K1l+kal Local Taylor expansion is z=5(1x2+2y)+hot (8.40) 1. Elliptic case: G=K1K2>0, see Figure 8.6(a) convex surface one side of tangent plane no infection 2. Hyperbolic case: G=K1K2 <0, see Figure 8.6(b) non-convex surface (locally intersection with tangent plane surface inflection Figure 8.6 :(a) Elliptic case; (b) Hyperbolic case Definition 1 Surface inflection at a point P exists iff the surface crosses the tangent plane at P If G<0. surface inflection at P IfG=0 in a region(e. g, a"cylinder"), then surface inflection exists in the region if H int)
8.2.4 Surface Fairing • Differential 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 (ρ1x2 z = + ρ2y2 ) + 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 inflection 2. Hyperbolic case: G = ρ1ρ2 < 0, see Figure 8.6(b). - non-convex surface (locally) - intersection with tangent plane - surface inflection Z κ1 κ2 Y Z Y X X Figure 8.6: (a) Elliptic case; (b) Hyperbolic case. Definition 1 Surface inflection at a point P exists iff the surface crosses the tangent plane at P. If G < 0, surface inflection at P. If G = 0 in a region ( e.g., a “cylinder” ), then surface inflection exists in the region if H changes sign in the region ( not one point)
Definition 2 Surface curve infection exists at a point P of a surface if a planar surface curve through P changes sign of its curvature at P Surface curve inflection Surface inflection Surface Interrogation 1. Gaussian curvature is not useful in cylinders"where G=0 2. Mean curvature is zero for minimal surfaces where K1= K2 3. Absolute curvature does not have such problems. Its curvature plots are useful in study ing surface imperfections( e.g., small oscillations 4. Reflection 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 defined for each of the colur Inns and rows of the control polyhedron
Definition 2 Surface curve inflection exists at a point P of a surface iff a planar surface curve through P changes sign of its curvature at P. Surface curve inflection ≤ Surface inflection • 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 studying surface imperfections ( e.g., small oscillations ). 4. Reflection 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 defined for each of the columns and rows of the control polyhedron