Chapter 4. Smoothing of Functions

Vladimir Rovenski, Haifa

> restart:

4.1 Method of Least Squares

Example 1. The data from the Example in Section 3.1.

> xx := [2,3,4,5,6,7,8]: yy:=[32,34,36,34,39,40,37]:

> PNTS := plot([xx, yy], style=point, symbol=circle):

> with(stats): P1 := op(2, fit[leastsquare[[x,y]]]([xx,yy]));

> P2 := op(2,fit[leastsquare[[x,y],y=a2*x^2+a1*x+a0]]([xx,yy]));

> P3 := op(2,fit[leastsquare[[x,y],y=a3*x^3+a2*x^2+a1*x+a0]]([xx, yy]));

> PL := plot([P1,P2,P3], x=0..10, color=[blue,green,red]):

> plots[display]([PNTS, PL]);

4.2 Bezier Curves

> xx:=[2,3,4,5,6,7,8]: yy:=[32,34,36,34,39,40,37]:

> PNTS:=plot([seq([xx[i], yy[i]],i=1..7)],style=point,symbol=circle): n:=nops(xx)-1:

> bez:=x->evalm(sum(binomial(n,i)*x^i*(1-x)^(n-i)*[xx[i+1],yy[i+1]],i=0..n)): PBEZ:=plot([bez(x)[1],bez(x)[2],x=0..1]):

> plots[display]([PNTS, PBEZ], axes=framed);

Note that the function bez(x)[1] is linear when the step of its net is constant.

> expand(bez(t)[1]);

Base Bernstein functions are nonnegative with sum equal to 1.

> simplify(sum(binomial(n,i)*x^i*(1-x)^(n-i), i=0..n));

4.3 Rational Bezier Curves

Extend the program for a Bezier curve (for example, = , and other weights are equal to 1).

>

> xx:=[2,3,4,5,6,7,8]: yy:=[32,34,36,34,39,40,37]:

> PNTS:=plot([seq([xx[i], yy[i]],i=1..7)],style=point,symbol=circle): n:=nops(xx)-1:

> w:=[1,1,5,1,5,1,1]: rbez:=x->evalm(sum(binomial(n,i)*x^i*(1-x)^(n-i)*w[i+1]*[xx[i+1], yy[i+1]], i=0..n)/sum(binomial(n,i)*x^i*(1-x)^(n-i)*w[i+1], i=0..n));

> PRBEZ:=plot([rbez(x)[1], rbez(x)[2], x=0..1]):

> plots[display]([PNTS, PRBEZ], axes=framed);

Example 1. Compare the curve PRBEZ given by 8 control points on the perimeter of the square with center

at the origin and the circle CIRC of radius with center at .

> xx:=[0,1,1,1,0,-1,-1,-1,0]: yy:=[-1,-1,0,1,1,1,0,-1,-1]:

> w:=[1,1,1.6,1,1.2,1,1.6,1,1]:

> CIRC:=plot([0.8*cos(t), 0.8*sin(t)-0.2, t=0..2*Pi], scaling=constrained):%;

Solution .

> n:=nops(xx)-1: PNTS:=plot([seq([xx[i], yy[i]],i=1..n+1)],style=point):

> rbez:=x->evalm(sum(binomial(n,i)*x^i*(1-x)^(n-i)*w[i+1]*[xx[i+1], yy[i+1]], i=0..n)/sum(binomial(n,i)*x^i*(1-x)^(n-i)*w[i+1], i=0..n)):

> PRBEZ:=plot([rbez(x)[1], rbez(x)[2], x=0..1]):

> plots[display]([PNTS, PRBEZ,CIRC], axes=framed);