Chapter 14. Spline Curves

Rovenski Vladimir, Haifa

> restart:

14.1 Preliminary Facts and Examples

Example.

1.

> r1:=[4*t,4*t]: r2:=[1+t,1+t]: t1:=0: t2:=1/4: t3:=0: t4:=1:

> p1:=plot([r1[1],r1[2],t=t1..t2]): p2:=plot([r2[1],r2[2],t=t3..t4]):

> c1:=plottools[circle](subs(t=t2, r1), 0.02):

> plots[display]([p1, p2, c1], scaling=constrained);

>

2.

> r1:=[-(1-t)^2, (1-t)^2]: t1:=0: t2:=1: r2:=[t^2, t^2]: t3:=0: t4:=1:

> p1:=plot([r1[1], r1[2], t=t1..t2]): p2:=plot([r2[1], r2[2], t=t3..t4]):

> c1:=plottools[circle](subs(t=t2, r1), 0.02):

> plots[display]([p1, p2, c1], scaling=constrained,axes=framed);

3.

> r1:=[sin(Pi/2*t^2),cos(Pi/2*t^2)]: t1:=0: t2:=1: r2:=[cos(Pi/2*t^2),-sin(Pi/2*t^2)]: t3:=0: t4:=1:

> p1:=plot([r1[1], r1[2], t=t1..t2]): p2:=plot([r2[1], r2[2], t=t3..t4]):

> c1:=plottools[circle](subs(t=t2, r1), 0.02):

> plots[display]([p1, p2, c1], scaling=constrained,axes=framed);

4.

> r1:=[cos(Pi/2*(1-t)^3),sin(Pi/2*(1-t)^3)]: t1:=0: t2:=1: r2:=[2-cos(Pi/2*t^3),-sin(Pi/2*t^3)]: t3:=0: t4:=1:

> p1:=plot([r1[1],r1[2],t=t1..t2]): p2:=plot([r2[1],r2[2],t=t3..t4]):

> c1:=plottools[circle](subs(t=t2, r1), 0.02):

> plots[display]([p1, p2, c1], scaling=constrained);

5.

> r1:=[t^2+t/2-3/2, -t^2+2*t]: t1:=0: t2:=1: r2:=[-t^2+5*t/2, -t^2+1]: t3:=0: t4:=1:

> p1:=plot([r1[1],r1[2],t=t1..t2]): p2:=plot([r2[1],r2[2],t=t3..t4]):

> c1:=plottools[circle](subs(t=t2, r1), 0.02):

> plots[display]([p1, p2, c1], scaling=constrained,axes=framed);

>

14.2 Composed Bezier Curves

14.2.1 Elementary Cubic Bezier Curve

> bez_2d := proc(xp, yp) local s,n,i,p,pp,m; n[0]:=t->(1-t)^3; n[1]:=t->3*t*(1-t)^2; n[2]:=t->3*t^2*(1-t); n[3]:=t->t^3; m:=nops(xp)/2; for s from 1 to 2*m do p[s-1]:=[xp[s], yp[s]] od; p[0]:=2*p[0]-p[1]; p[2*m-1]:=2*p[2*m-1]-p[2*m-2]; pp:=n[0](t-i+1)*(p[2*i-2]+p[2*i-1])/2+n[1](t-i+1)*p[2*i-1]+n[2](t-i+1)*p[2*i]+n[3](t-i+1)*(p[2*i]+p[2*i+1])/2; RETURN(evalm(sum((1+signum(-t+i))*(1+signum(t-i+1))/4*pp, i=1..m-1))) end:

Example 2 .

1. Plot the plane composed Bezier curve using the command bez_2d(xp, yp) .

> xp:=[.056, .287, .655, .716, .228, .269, .666, .929]: yp:=[.820, .202, .202, .521, .521, .820, .820, .227]:

> polygon:=plot([seq([xp[j], yp[j]], j=1..nops(xp))]):

> bez:=bez_2d(xp,yp): curve:=plot([bez[1],bez[2], t=0.001..nops(xp)/2-1.001]): plots[display]([curve, polygon]);

Continuation . We plot the self-repeating composed Bezier curve, where the circled points are derived by the program

> x1:=[op(xp),seq(seq(xp[i]+.61*j,i=3..nops(xp)),j=1..3)]; y1:=[op(yp),seq(seq(yp[i], i=3..nops(yp)),j=1..3)];

> P1:=plot([seq([x1[j], y1[j]], j=1..nops(x1))],scaling=constrained): %;

> bez1:=bez_2d(x1, y1): curve1:=plot([bez1[1], bez1[2], t=0.001..nops(x1)/2-1.001]): plots[display]([curve1, P1]);

3. Plot closed Bezier curves of two types:

(a) with the loop having a singular point, k in N.

> k:=3: xp:=[seq(sin(2*Pi*i/(2*k-1)), i=1..2*k)]: yp:=[seq(cos(2*Pi*i/(2*k-1)), i=1..2*k)]:

> polygon:=plot([seq([xp[j], yp[j]], j=1..nops(xp))]):

> bez:=bez_2d(xp, yp): curve:=plot([bez[1], bez[2], t=0.001..nops(xp)/2-1.001]): plots[display]([curve, polygon]);

(b)

> k:=3: xp:=[seq(cos(2*Pi*i/(2*k)), i=1..2*k+2)]: yp:=[seq(sin(2*Pi*i/(2*k)), i=1..2*k+2)]: xp[1]:=(xp[2]+xp[2*k+1])/2: yp[1]:=(yp[2]+yp[2*k+1])/2: xp[2*k+2]:=xp[1]: yp[2*k+2]:=yp[1]:

> polygon:=plot([seq([xp[j], yp[j]], j=1..nops(xp))]):

> bez:=bez_2d(xp,yp): curve:=plot([bez[1],bez[2], t=0.001..nops(xp)/2-1.001]): plots[display]([curve,polygon]);

14.3 Composed B -Spline Curves

Here is a procedure for deriving and plotting given data in the plane.

> bspl_2d := proc(xp,yp) local s,n,i,p,pp,m; n[0]:=t->(1-t)^3/6; n[1]:=t->(3*t^3-6*t^2+4)/6; n[2]:=t-> (-3*t^3+3*t^2+3*t+1)/6; n[3]:=t->t^3/6; m:=nops(xp)-1; for s from 1 to m+1 do p[s-1]:=[xp[s], yp[s]] od; p[-1]:=p[0]; p[m+1]:=p[m]; p[-2]:=p[0]; p[m+2]:=p[m]; pp:=n[0](t-i+1)*p[i-1]+n[1](t-i+1)*p[i]+n[2](t-i+1)*p[i+1]+n[3](t-i+1)*p[i+2]; RETURN(evalm(sum((1+signum(-t+i))*(1+signum(t-i+1))/4*pp,i=-1..m))) end:

Use the procedure bspl_2d(xp,yp) to plot the plane curve.

> xp:=[.056, .287, .655, .716, .228, .269, .666, .929]: yp:=[.820, .202, .202, .521, .521, .820, .820, .227]:

> polygon:=plot([seq([xp[j], yp[j]], j=1..nops(xp))]):

> bs:= bspl_2d(xp,yp): curve:=plot([bs[1], bs[2], t=-1.999..nops(xp)-1.001]): plots[display]([curve,polygon]);

P rocedure of deriving and plotting B-spline in .

> bspl_3d:=proc(xp,yp,zp) local n,i,s,pp,p; global spl,m,t; n[0]:=t->(1-t)^3/6; n[1]:=t->(3*t^3-6*t^2+4)/6; n[2]:=t->(-3*t^3+3*t^2+3*t+1)/6; n[3]:=t->t^3/6; m:=nops(xp)-1; for s from 1 to m+1 do p[s-1]:=[xp[s], yp[s], zp[s]] od; p[-1]:=p[0]; p[m+1]:=p[m]; p[-2]:=p[0]; p[m+2]:=p[m]; pp:=n[0](t-i+1)*p[i-1]+n[1](t-i+1)*p[i]+n[2](t-i+1)*p[i+1]+n[3](t-i+1)*p[i+2]; spl:=evalm(sum((1+signum(-t+i))*(1+signum(t-i+1))/4*pp, i=-1..m)); RETURN(plots[spacecurve](spl, t=-1.999..m-0.001)) end:

Use the procedure bspl_3d(xp, yp, zp) to plot the space curve.

> xp:=[.056, .287, .655, .716, .228, .269, .666, .929]: yp:=[.820, .202, .202, .521, .521, .820, .820, .227]: zp:=[1/4, 1/4, 0, 0, 0, 0, -1/4, -1/4]:

> polygon:=plots[pointplot3d]([seq([xp[j], yp[j], zp[j]], j=1..nops(xp))],symbol=circle): bs:=bspl_3d(xp,yp,zp):

> plots[display]([bs, polygon], scaling=constrained,axes=boxed);

Plot the knot trefolium

> xp:=[2, 0, -1, 0, 1, 0, -2, -1, 1, 2, 2]: yp:=[0, 0, 2, 4, 2, 0, 0, 2, 2, 0, 0]: zp:=[0, 1, 2, 2, 0, 2, 2, 1, 2, 1, 0]:

> bs:=bspl_3d(xp, yp, zp):%;

> plots[tubeplot]([spl[1], spl[2], spl[3], radius=.2], t=-2..nops(xp)-1, tubepoints=50, style=patchnogrid, color=[1, .8, .1], light=[75,50,1,.9,.5], ambientlight=[.4,.4,.4],axes=framed);

>

14.4 -Spline Curves

14.4.1 Examples of -Spline Curves

Procedure for deriving and plotting in .

> restart:

> beta_2d:=proc(xp,yp,beta1,beta2) local s,i,pp; global d,b1,b2,p,m,n; b1:=beta1; b2:=beta2; d:=evalf(2*b1^3+4*b1^2+4*b1+b2+2); m:=nops(xp)-1; n[0]:=t->2*b1^3*(1-t)^3/d; n[1]:=t->(2*b1^3*t*(t^2-3*t+3)+2*b1^2*(t^3-3*t^2+2)+2*b1*(t^3-3*t+2)+b2*(2*t^3-3*t^2+1))/d; n[2]:=t->(2*b1^2*t^2*(-t+3)+2*b1*t*(-t^2+3)+b2*t^2*(-2*t+3)+2*(-t^3+1))/d; n[3]:=t->2*t^3/d;for s from 1 to m+1 do p[s-1]:=[xp[s], yp[s]] od; p[-1]:=p[0]; p[m+1]:=p[m]; p[-2]:=p[0]; p[m+2]:=p[m]; pp:=n[0](t-i+1)*p[i-1]+n[1](t-i+1)*p[i]+n[2](t-i+1)*p[i+1]+n[3](t-i+1)*p[i+2]; RETURN(evalm(sum((1+signum(-t+i))*(1+signum(t-i+1))/4*pp,i=-1..m))) end:

Use the procedure beta_2d(xp,yp,zp) to plot the curve

> xp:=[.056, .287, .655, .716, .228, .269, .666, .929]: yp:=[.820, .202, .202, .521, .521, .820, .820, .227]:

> P2:=plot([seq([xp[i], yp[i]], i=1..nops(xp))]):

> Bs:=beta_2d(xp,yp,1,7): curve:=plot([Bs[1],Bs[2],t=-2+0.001..nops(xp)-1.001]): plots[display]([curve,P2]);

Change in a Beta-spline when one of its control points (with number i_0) moves.

> i0:=4: delta:=[rand(0..10)()/40, rand(0..10)()/40]:

> char:=i -> (1+signum(-t+i))*(1+signum(t-i+1))/4:

> pnew:=p[i0]+delta: for i from i0-2 to i0+1 do pp:=n[0](t-i+1)*p[i-1]+n[1](t-i+1)*p[i]+n[2](t-i+1)*p[i+1]+n[3](t-i+1)*p[i+2]: old[i]:=pp: new[i]:=subs(p[i0]=evalm(pnew),pp) od: r_d:=sum(char(j)*(new[j]-old[j]),j=i0-2..i0+1): r_new:=evalm(Bs+r_d):

> new_curve:=plot([r_new[1],r_new[2], t=-1.99..m-0.01],thickness=3):%;

> for i from -2 to m+2 do if i<> i0 then q[i]:=p[i] else q[i]:=pnew fi od:

> new_polygon:=plot([seq(q[i], i=0..m)]):%;

> poly2:=plot([[p[i0-1], pnew, p[i0+1]]]):

> plots[display]([new_polygon,P2,new_curve,curve]);

Play with beta_2d using the following data a) - b).

a) An example of random generating of an array of points in the plane.

> mm:=7: xp:=[seq(rand(1..10)(), j=1..mm)]: yp:=[seq(rand(1..10)(), j=1..mm)]: P:=plot([seq([xp[i], yp[i]], i=1..nops(xp))]):

b) Data for car profile

> xp:=[226,348,356,400,410,452,470,452,438,496,416,412,430,428,394,344,232,164,58,38,36,98,100,148,158,228,226]:

> yp:=[21,17,5,4,18,24,55,62,82,89,83,79,79,63,68,94,87,62,53,34,21,16,2,1,13,21,21]:

> Bs:=beta_2d(xp,yp,1,7): curve:=plot([Bs[1],Bs[2],t=-2+0.001..nops(xp)-1.001]): plots[display]([curve,P],scaling=constrained);

Analogous procedure for deriving and plotting in .

> beta_3d := proc(xp,yp,zp,beta1,beta2) local s,n,i,p,pp,m; global spl,d,b1,b2; b1:=beta1; b2:=beta2; d:=evalf(2*b1^3+4*b1^2+4*b1+b2+2);n[0]:=t->2*b1^3*(1-t)^3/d;n[1]:=t->(2*b1^3*t*(t^2-3*t+3)+ 2*b1^2*(t^3-3*t^2+2)+2*b1*(t^3-3*t+2)+b2*(2*t^3-3*t^2+1))/d; n[2]:=t->(2*b1^2*t^2*(-t+3)+2*b1*t*(-t^2+3)+b2*t^2*(-2*t+3)+2*(-t^3+1))/d; n[3]:=t->2*t^3/d;m:=nops(xp)-1; for s from 1 to m+1 do p[s-1]:=[xp[s], yp[s], zp[s]] od; p[-1]:=p[0]; p[m+1]:=p[m]; p[-2]:=p[0]; p[m+2]:=p[m]; pp:=n[0](t-i+1)*p[i-1]+n[1](t-i+1)*p[i]+n[2](t-i+1)*p[i+1]+n[3](t-i+1)*p[i+2]; spl:=evalm(sum((1+signum(-t+i))*(1+signum(t-i+1))/4*pp, i=-1..m));RETURN(plots[spacecurve](spl, t=-1.999..m-0.001)) end:

Use the procedure beta_3d(xp, yp, zp) to plot the space curve.

> xp:=[.056, .287, .655, .716, .228, .269, .666, .929]: yp:=[.820, .202, .202, .521, .521, .820, .820, .227]: zp:=[1/4,1/4,0,0,0,0,-1/4,-1/4]:

> polygon:=plots[pointplot3d]([seq([xp[j], yp[j], zp[j]],j=1..nops(xp))],symbol=circle): bs:=beta_3d(xp, yp, zp, 1, 1):

> plots[display]([bs,polygon],scaling=constrained,axes=boxed);

>

14.5 Interpolation Using Cubic Hermite's Curves

14.5.1 Elementary Cubic Hermite's Curve

Procedure for deriving and plotting in .

> herm1_2d:=proc(xp,yp,xq,yq) local s,n,i,m,p,q,pp; global spl; n[0]:=t->2*t^3-3*t^2+1; n[1]:=t->t^2*(-2*t+3); n[2]:=t->t*(t^2-2*t+1); n[3]:=t->t^2*(t-1); m:=nops(xp)-1; for s from 1 to m+1 do p[s-1]:=[xp[s], yp[s]]; q[s-1]:=[xq[s], yq[s]] od; pp:=n[0](t-i+1)*p[i-1]+n[1](t-i+1)*p[i]+n[2](t-i+1)*q[i-1]+n[3](t-i+1)*q[i]; RETURN(evalm(sum((1+signum(-t+i))*(1+signum(t-i+1))/4*pp,i=1..m))) end:

Use the procedure herm1_2d(xp,yp,xq,yq) for ploting the plane curve.

> xp:=[.056, .287, .655, .716, .228, .269, .666, .929]: yp:=[.820, .202, .202, .521, .521, .820, .820, .227]:

> xq:=[.25, .25, -.25, -.25, .25, .25, -.25, -.25]: yq:=[-.25, .25, -.25, .25, -.25, .25, -.25, .25]:

> P:=plot([seq([xp[j], yp[j]], j=1..nops(xp))]):

> vec:=i->plottools[arrow]([xp[i-1], yp[i-1]], [xq[i-1], yq[i-1]], .001,.02,.1): V:=seq(vec(i+1),i=1..nops(xp)):

> H1:=herm1_2d(xp,yp,xq,yq): curve:=plot([H1[1],H1[2],t=.01..nops(xp)-1.01]): plots[display]([curve,P,V],scaling=constrained);

>

14.5.2 Composed Cubic Hermite's Curve

Procedure for deriving and plotting in . The above system is solved in the procedure.

> herm2_2d := proc(xp, yp, xq, yq) local eq,a,b,ss,s,n,i,j,m,p,q,pp; n[0]:=t->2*t^3-3*t^2+1; n[1]:=t->t^2*(-2*t+3); n[2]:=t->t*(t^2-2*t+1); n[3]:=t->t^2*(t-1); m:=nops(xp)-1; p:=array(1..m+1); q:=array(1..m+1); a:=array(1..m-1, 1..m+1); b:=array(1..m-1, 1..m+1); for i from 1 to m-1 do for j from 1 to m+1 do if j=i then a[i,j]:=1; b[i,j]:=-3 elif j=i+1 then a[i,j]:=4; b[i,j]:=0 elif j=i+2 then a[i,j]:=1; b[i,j]:=-3 else a[i,j]:=0; b[i,j]:=0 fi od od; eq:=evalm(a&*q-b&*p); ss:=solve({seq(eq[i]=0, i=1..m-1)},{seq(q[i], i=2..m)}); assign(ss); q[1]:=[xq[1],yq[1]]; q[m+1]:=[xq[2],yq[2]]; for s from 1 to m+1 do p[s]:=[xp[s],yp[s]] od; i:='i'; pp:=n[0](t-i+1)*p[i]+n[1](t-i+1)*p[i+1]+n[2](t-i+1)*q[i]+n[3](t-i+1)*q[i+1]; RETURN(evalm(sum((1+signum(-t+i))*(1+signum(t-i+1))/4*pp,i=1..m))) end:

Plot the plane curve.

> xp:=[.056, .287, .655, .716, .228, .269, .666, .929]: yp:=[.820, .202, .202, .521, .521, .820, .820, .227]:xq:=[.1, -.1]: yq:=[.1, -.1]:

> vecs:=plots[display](plottools[arrow]([xp[1],yp[1]],[xq[1],yq[1]], 0.001, 0.02, 0.1), plottools[arrow]([xp[nops(xp)],yp[nops(yp)]],[xq[2],yq[2]],0.001,0.02,0.1)):

> polygon:=plot([seq([xp[j], yp[j]], j=1..nops(xp))]):

> H2:=herm2_2d(xp,yp,xq,yq): curve:=plot([H2[1],H2[2],t=0.01..nops(xp)-1.01]): plots[display]([vecs,curve,polygon]);

14.6 Composed Catmull-Rom Spline Curves

Procedure for deriving and plotting in .

> restart:

> crom_2d := proc(xp, yp) local s,n,i,p,pp,m; n[0]:=t-> -t*(1-t)^2/2; n[1]:=t-> (3*t^3-5*t^2+2)/2; n[2]:=t-> t*(-3*t^2+4*t+1)/2; n[3]:=t-> t^2*(t-1)/2; m:=nops(xp)-1; for s from 1 to m+1 do p[s-1]:=[xp[s], yp[s]] od; pp:=n[0](t-i+1)*p[i-1]+n[1](t-i+1)*p[i]+n[2](t-i+1)*p[i+1]+n[3](t-i+1)*p[i+2]; RETURN(evalm(sum((1+signum(-t+i))*(1+signum(t-i+1))/4*pp, i=1..m-2))) end:

Plot a plane curve.

> xp:=[.056, .287, .655, .716, .228, .269, .666, .929]: yp:=[.820, .202, .202, .521, .521, .820, .820, .227]:

> P3:=plot([seq([xp[j], yp[j]], j=1..nops(xp))],thickness=2):

> for i from 1 to nops(xp)-2 do v[i-1]:= plots[display](plottools[arrow]([xp[i], yp[i]], [xp[i+2]-xp[i], yp[i+2]-yp[i]], 0.001, 0.01, 0.1)) od:

> vecs:=seq(v[i-1], i=1..nops(xp)-2): cr:=crom_2d(xp,yp):

> curve1:=plot([cr[1], cr[2], t=0.01..nops(xp)-3.01],thickness=2):

> plots[display]([curve1, P3, vecs],scaling=constrained);

>