Curves

1. Curves • Path - trajectory - motion • Parameterization • Designing polyloops • Velocity, tangent, curvature, jerk • Smoothing polyloops • B-spline, Bezier and 4-points subdivisions • The Jarek subdivision schemes for polyloops • Arc-length parameterization • Tangent and curvature estimation • Curvature-driven adaptive parameterization Updated January 1, 2020

2. How to specify a curve? • Implicit (equation): x2+y2=r2 • Not convenient for walking on the curve • Parametric (function): P(s) = ( r cos(s) , r sin(s) ) • Polyloop: cycles of edges joining consecutive points (vertices) • Not smooth,unless you use lots of vertices • Subdivision: result to which a subdivision process converges • Procedural: snowflake, fractal, subdivision… V2 V3 V1 V5 V4 V0

3. It is all about curves Points, vectors curves Shapes Motions Animated shapes

4. Outline • We will talk about polygonal loops (polyloops) in the plane • We will discuss algorithms to transform them into smoothcurves • We want to support an intuitive design paradigm: • A designer adjust the vertices of the control polygon (polyloop) • The corresponding smooth curve is regenerated automatically • The effect of an adjustment should be local and predictible • We will look at several subdivision schemes for converting a control polyloop into a smooth curve • We will study their properties and look at several applications • Smooth rendering, culling and clipping • Computing intersections, areas • We will extend these results to • 3D curves • Surfaces • Animations

5. Geometry review • The constructions and algorithms that follow are formulated in terms of simple geometric constructions that involve points, vectors, and operators. • We will first • provide an overview of these fundamental concepts • propose a few convenient operations • define our notation

6. P=(4,3) J O I Points • A point is a location in space (zero dimensional entity) • Initially space will the plane (later we will go 3D) • How to represent a point P? • Assume a known origin O and two orthogonal directions I and J (unit basis vectors) • Represent P by its two coordinatesP=(x,y) where P=O+xI+yJ • To get to P: • start at O, • walk x units along the direction I, • then y units along the direction J • REMEMBER THIS EXPLANATION! • It is the key to understanding • Coordinate systems • Graphic transformations • Animations

7. Vectors • A vector is either a direction or a displacement between two points (locations). • We will also use vectors to represent forces and velocities. • A vector is often represented by its coordinates: V=(x,y) • Please do not confuse points and vectors! • A direction is a unit vector. It may represent • the direction of a ray from the eye through a pixel • the direction towards the light from a point on a surface • the tangent to a curve at a specific point • a basis vector for expressing the coordinates of points • Displacement vectors may be computed as point subtractions • Let V be the displacement vector from A=(xA,yA) to B=(xB,yB) • V may be denoted B–A, or for short AB • Its coordinates may be computed as V =(xB–xA, yB–yA) • Hence, given O, a point P corresponds to the vector OP • Assuming that the coordinates of the origin O are (0,0)

8. V+V’ V’ V –V Operations on vectors • Consider two vectors: V=(x,y) and V’=(x’,y’) • Scaling vectors sV= (sx,sy) • Reversing a vector: –V=(–x,–y) • Adding vectors V+V’= (x+x’,y+y’) • Subtracting vectors V–V’= (x–x’,y–y’) • What do the coordinate mean? • Consider a basis of two orthogonal unit vectors I and J • V=(x,y) means that V=xI+yJ • Walk x units along I and y units along J

9. U U.away(V) U.along(V) V “Advanced” vector operators and notation • Let V=(x,y) and U=(a,b) • Length of V is computed by V.norm=sqrt(x2+y2) • Unit length vector obtained by normalizing V: V.unit=(1/V.norm)V, • Orthogonal vector (rotating V ccw by 90 degrees) V.left=(–y,x) • Cascading “.” notation (left-to-right): V.left.norm=(V.left).norm • Dot product: UV=ax+by (a scalar!) • When V.norm=U.norm=1, then UV is the projection of U onto V • Proof: write U=cV+dV.left=(cx–dy,dx+cy), solving for c with x2+y2=1 yields c=ax+by • Otherwise: UV is the projection of U onto V scaled by the length of V • UV = U.norm*V.norm*cos(angle(U,V)) • Properties of dot product • UV = VU • (–U)V = –(VU) • UV = 0 implies that they are orthogonal or one has zero length • Tangential and orthogonal decomposition U with respect to a vector V: • along V: U.along(V):=(UV.unit)V.unit • awayfrom V: U.away(V):=U–U.along(V)

10. Walking along an edge • Consider an edge from point A to point B: denoted Edge(A,B) • You can use a parameter, s, in [0,1] to identify any point P(s) on it P(s)=A+sAB • This is the proper notation: start at A and move by s along vector AB • For convenience, we often use another equivalent notation: P = A+s(B–A) = (1–s)A+sB • Note that P(0)=A and P(1)=B 1–s B s P(0.6) A

11. P P AP B AP H B H A A s=1.4 s=0.6 Closest point on an edge • Consider Edge(A,B) and a point P in the plane • Compute the point C of Edge(A,B) that is the closest to P • First compute the parameter s of the orthogonal projection H of P on Edge(A,B) • Then test it against [0,1] and decide Procedure project (A,B,P) { s:=APAB/ABAB; # assuming AB IF 0<s<1 THEN RETURN A+sAB ELSE IF s<0 THEN RETURN A ELSE RETURN B } To compute closest point to P on polyloop: Compute it for all edges and take closest

12. Arclength of an edge • Write the procedure L(A,B) that returns the arclength of edge (A,B)

13. Total arclength of a polyloop • Write the procedure L(C) that returns the total arclength of polyloop with vertices C[i], for i = 0, 1, … n

14. Terminology (for lack of a better one) • A path is a parameterized curve (2D or 3D) • A point (of an object) may follow a curve • The path indicates where we are along the curve at each moment of the animation • The orientation of the object is not defined (could be pure translation?) • A motion defines the evolution of a pose • Position and orientation of the local CS or object(s) defined in it • A trajectory is a motion defined by a path • The orientation is defined by the tangent to the curve • …and bi-tangent in 3D…

15. Curve and path Given a path P(t), what is the curve traced by P(t)? Union of all points P(t) What if P(t) starts at A, moves a few times back and forth along the segment (A,B) then finishes at B? The path is complicated The curve is trivial

16. Motivation • Assume that you have a polyloop P path defined by sampling positions Pi or produced by subdividing a control polygon L • In general, the vertices of P are not uniformly spaced • Their spacing may reflect the animator’s desire to slow down near turns or at specific places • or it may be the undesired side effect of data acquisition or creation • We may want to re-parameterize P to obtain a path P’ that moves at constant speed if you use s(P’i,t,P’i+1) on each segment • Then, you can warp time (as in project P1a) to impose your own speed profile, independently on the original sampling of P

17. 4-pt subdivision slower faster resampling constant speed Example where resampling is needed Subdivision of a bad control polygon

18. P(1)=P(0) t P(.5) Arclength reparameterization • Write the procedure for computing a point on the polyloop that corresponds to parameter t in [0,1] • Demonstrate it by animating a point at constant speed

19. Research challenge!!! Given a polyloop P, place an ordered set of k points along P so that each point is at the same Euclidean (shortest) distance to its two neighbors. Is is always possible? Is there a deterministic algorithm? Would an iterative algorithms always converge?

20. Polygon area calculation: 2 methods • Sum of signed areas of triangles, each joining an arbitrary origin o to a different edge(a,b) • SUM (oaob) for each edge (a,b) • uv is the dot product of v with the the vector obtained by rotating u by 90° • Sum of signed areas of trapezoids between each oriented edge and its orthogonal projection (shadow) on the x-axis y by b ay a x ax bx area(a,b):=(ay+by)(bx–ax)/2

21. Area implementation float area () { float A=0; for (int i=0; i<vn; i++) A+=trapezeArea(P[i],P[this.in(i)]); return(A); } float trapezeArea(pt A, pt B) {return((B.x-A.x)*(B.y+A.y)/2.);}

22. Centroid of polygonal region • Weighted sum of centroids of trapezoids • Weight = signed area of trapezoid / total area

23. Centroid of trapezoid Any line through centroid cuts shape in 2 parts with same moment x = 0∫1 sy(s)ds / 0∫1 y(s)ds, with y(s)=a+s(b–a) x = 0∫1 (b–a)s2 + as ds / 0∫1 (b–a)s +a ds x=(a+2b)/(3(a+b)) Similar derivation: y= (a2+ab+b2)/(3(a+b)) y b a (x,y) x 0 1

24. Implementation of centroid pt barycenter () {G.setTo(0,0); for (int i=0; i<vn; i++) G.addScaledPt( trapezeArea(P[i],P[this.in(i)]),trapezeCenter(P[i],P[this.in(i)])); G.scaleBy(1./this.area()); return(G); } pt trapezeCenter(pt A, pt B) { return(new pt( A.x+(B.x-A.x)*(A.y+2*B.y)/(A.y+B.y)/3., (A.y*A.y+A.y*B.y+B.y*B.y)/(A.y+B.y)/3.) ); }

25. Principal directions of polyloop • For a set of points sampled along the curve • Not the same as principal axes (second order inertia moments)

26. Implementation of principal direction void showAxis() {pt G=this.barycenter(); stroke(black); G.show(10); float xx=0, xy=0, yy=0, px=0, py=0, mx=0, my=0; for (int i=0; i<vn; i++) {xx+=(P[i].x-G.x)*(P[i].x-G.x); xy+=(P[i].x-G.x)*(P[i].y-G.y); yy+=(P[i].y-G.y)*(P[i].y-G.y);}; float a = atan(2*xy/(xx-yy))/2.; vec V= new vec(50,0); V.rotateBy(a); vec U = V.makeTurnedLeft(); for (int i=0; i<vn; i++) {vec W=G.makeVecTo(P[i]); float vd=dot(W,V); float ud=dot(W,U); if(vd>0) px+=vd; else mx-=vd; if(ud>0) py+=ud; else my-=ud; }; if(mx>max(px,py,my)) V.rotateBy(PI); if(py>max(px,mx,my)) V.rotateBy(PI/2.); if(my>max(px,mx,py)) V.rotateBy(-PI/2.); strokeWeight(1); stroke(black); V.showArrowAt(G); stroke(black,60); V.turnLeft(); V.showAt(G); V.turnLeft(); V.showAt(G); V.turnLeft(); V.showAt(G); }

27. Questions • Describe two techniques for testing whether a point q lies inside a polygon P • Describe two techniques for computing the area of a polygon P • Assume that you have a function fill_tri(a,b,c) that will visit all pixels whose center falls inside the triangle with vertices (a,b,c). Assume that each time the function visits a pixel, it toggles its status. Initially, all the status of each pixel is OFF. Explain how you could turn to ON the status of all pixels whose center falls inside a polygon P.

28. Motivation • Smooth curves and surfaces are used for aesthetic, manufacturing, and analysis applications where discontinuities due to piecewise linear approximations would create misleading artifacts. • Designers like to specify such curves by adjusting a control network with only a few vertices • How to go from a crude 2D control polyloop (black) to a smooth curve (red)? • Suggestions?

29. compromise Different requirements smooth interpolating

30. Operators on polyloops Simplification Subsampling Smoothing Denoising Subdivision Refinement Adding noise

31. Motivation for curve smoothing • Smooth curves and surfaces are used for aesthetic, manufacturing, and analysis applications where discontinuities due to triangulated approximations would create misleading artifacts. • Smooth curves may be used to design trajectories and paths for movingobjects and cameras in animations. • Designers like to specify such curves by adjusting a crude control polygon with only a few vertices • How to go from a crude control polygon (black) to a smooth curve (red)? • Suggestions?

32. Local curve analysis • We first investigate local differential properties of a polyloop L • to measure forces, compare schemes, adjust speeds… • … but these are trivial for a polyloop • straight line segments (zero curvature) • not differentiable at vertices (infinite curvature) • So… we really want the properties of the smooth curve or path J defined by L • tangent, normal, acceleration, curvature, jerk • … but, we usually do not know J • it may be defined as the limit of a smoothing process • Therefore, we use discreteestimators for these properties

33. AB B V(B) L C A V(B) Velocity at a vertex: V(B)=(AB+BC)/2 • Let A,B,C,D,E… be consecutive vertices in L • AB, BC,… are the edges of L • J is the unknown smooth curve passing through the vertices • When traveling along J from A to B, it takes 1 sec. • Hence the average velocity vector for that segment is AB. • Hence the velocity at B is estimated as V(B)=(AB+BC)/2

34. AB B V(B) L C A V(B) Implementing V(B)=(AB+BC)/2 class vec { float x,y; vec average(vec U, vec V) {return(new vec((U.x+V.x)/2,(U.y+V.y)/2)); }; void unit() {float n=sqrt(sq(x)+sq(y)); if (n>0.000001) {x/=n; y/=n;};}; … class pt { float x,y; vec vecTo(pt P) {return(new vec(P.x-x,P.y-y)); }; Write a method for the pt class that computes the velocity at point B vec velocity (pt A, pt C) {…}; tangent T(B) = V(B).unit()

35. N(B) B V(B) C A Normal at a vertex: N(B)=AC.left().unit() • The normal N(B) at B is the unit vector orthogonal to AC. • We pick the one pointing left: N(B)=AC.left().unit() class vec { float x,y; void left() {float w=x; x=-y; y=w;}; void unit() {float n=sqrt(sq(x)+sq(y)); if (n>0.000001) {x/=n; y/=n;};}; Write a method for the pt class that computes the normal at point B vec normal (pt A, pt C) {…};

36. B V(B) C g(B) A Acceleration at a vertex: g(B)=BC–AB • The acceleration at B is the velocity change g(B)=BC–AB • It is useful for computing forces during animation (dynamics) Write a method for the pt class that computes the acceleration at point B vec acceleration (pt A, pt C) {…};

37. Radius of curvature: r(B)=V2/(2ABN(B)) • The radius r(B) of curvature (sharpness of turn in the path) at B could be estimated as the radius of the circle through A, B, and C, but this approach yields unexpected results when the angle at b is acute. • I prefer to use the parabolic curvature, r(B)=V2/(2h), where V is the magnitude of V(B) and where h is the distance from B to the Line(A,C). • h=ABN(B) is the normal component of g(B) • It measures the centrifugal force • Curvature = 1/r(B) Write a method for the pt class that computes the r(B) at point B vec radiusOfCurvature (pt A, pt C) {…};

38. j(B,C) B C –g(B) g(C) g(B) A D Jerk: j(B,C)=AD+3CB It is the change of acceleration between B and C. It measures the change in the force felt by a person traveling along the curve j(B,C) = g(C)–g(B) = (CD–BC)–(BC–AB) = CD+CB+CB+AB = AB+BC+CD + CB+CB+CB = AD+3CB To show the second degree discontinuities in the curve, draw the normal component of j(B,C) at (B+C)/2

39. smoothing Smoothing a polyloop Given a polyloop L, we wish to move its vertices to produce a curve J that resembles L, but is “smoother” (less jerk?) Smoothing does not change the number of vertices, it only moves them

40. tuck(.5) tuck(.5) tuck(.5) tuck(.5) Tuck smoothing tuck(s) moves the vertices by executing the following two steps: • For each vertex B, compute M(B)=(BA+BC)/2 • A (r. C) is the previous (r. next) vertex in the cyclic order • Moving B to B+M(B) brings it to the average of its neighbors • For each vertex B, move B to B+sM(B) Repeating tuck(s) tends to eliminate sharp features, but shrinks the curve …

41. Avoid oscillations Keep s in [0,2/3] to avoid oscillations tuck(1) tuck(1)

42. Avoid shrinking: Tuck-tuck smoothing • Alternate tuck(s), which shrinks, and tuck(–s), which expands • untuck cannot recover the sharp features lost during the tuck

43. dual dual Tuck(1/2) = dual of dual • Tuck(.5) is the dual of the dual • refine introduces a new vertex in the middle of each edge • kill deletes the old vertices • dual = (refine, kill) • tuck(.5) = (dual, dual) refine kill

44. C(0) C(1) C(-1) M B D C(2) C(-2) A E M’ smoothing bulge M’M/3 Cubic fit tuck-tuck: (–A+4B+4D–E)/6 Compute a cubic polynomial curve, C(t)=at3+bt2+ct+d, through vertices A, B, D, E satisfying: C(–2)=A, C(2)=E, C(–1)=B, and C(1)=D We want to insert a point C=C(0)=d between B and E. Solve for d C(–2)=–8a+4b–2c+d=A and C(2)=8a+4b+2c+d=E yields A+E=8b+2d C(–1)=–a+b–c+d=B, and C(1)=a+b+c+d=D yields A+E=8b+2d A+E=8b+2d and A+E=8b+2d yields 6d=4(B+D)–(A+E) d = (–A+4B+4D–E)/6 = M+M’M/3, with M=(B+D)/2, M’=(A+E)/2) To insert such a cubic fit for each edge: tuck(2/6), tuck(–2/6) Note that 2/6=0.82, which may be unstable Move towards C

45. Smoothing vs subdivision • Smoothing should be performed when the desired curve has been densely sampled but the resulting polyloop has undesiredsharp features (acquisition noise, too many details, sharp changes in direction or curvature) • Subdivision refines the polyloop (adds more vertices). It should be performed when the polyloop is a crude outline of the desired curve. subdivision smoothing

46. Polyloop refine/smooth operators We will use the following operators: • r (refine) introduces a new vertex in the middle of each edge • k (kill) deletes the old vertices (those present before r) • d (dual) performs (r, k) • t (tuck) performs tuck(1/2), which is (d, d) • us (parameterized untuck) performs tuck(–s) • u (untuck) performs tuck(–1/2), same as u.5 Write a polyloops smoothing scheme using these

47. Quadratic uniform B-spline (r,d) Repeating (r,d) converges to a piecewise quadratic B-spline curve that interpolates the mid-edge points r rd rdrd rdrdrdrd

48. Cubic uniform B-spline (r,t) Repeating (r,t), which is (r,d,d) converges to a piecewise cubic B-spline curve that has second degree continuity (no jerk) rt rtr rtrt rtrtrt

49. 2 B A 3 1 4 8 5 7 C D 6 Split&tweak B-spline subdivision • Given a control polygon, for example (A,B,C,D), repeat the following sequence of two steps, until all consecutive 4-tuples of control points are nearly co-linear. 1. “Split”: insert a new control point in the middle of each edge (2,4,6,8) 2. “Tweak”: move the old control points half-way towards the average of their new neighbors (1,3,5,7) • The polyloop converges rapidly to a smooth curve, which happens to be a cubic B-spline curve.

50. Can the limit curve self-intersect? • Can it self-intersect when P does not?