spalldem.html

演示matlab曲线拟和与插直的基本方法

%  knots(1:n+k)                    0     0     0     1     2     3     3     3
%  coefficients(d,n)              -1     0     1     0    -1
%  number n of coefficients  5
%  order k                            3
%  dimension  d  of target     1
%
% However, there is usually NO NEED to know any of these parts. Rather,
% one uses commands like  SP = SPAPI(...)  or  SP = SPAPS(...)  to
% construct the B-form of a spline from some data, then uses commands like
% FNVAL(SP,...), FNPLT(SP,...), FNDER(SP)  etc., to make use of the spline
% constructed, without any need to look at its various parts.
%
% The next slides give more detailed information about the  B-splines, in
% particular about the important role played by knot  MULTIPLICITY .

b = 0:3;
sp = spmak(augknt(b,3),[-1,0,1,0,-1]);
x = linspace(-1,4,51);
plot(x,fnval(sp,x),'x'),
hold on,
axis([-2 5,-1.5,1]),
tx = repmat(b,3,1);
ty = repmat(.15*[1;-1;NaN],1,length(b));
plot(tx(:),ty(:),'-r')
hold off

%%
% Here, for  k=2,3,4, is a B-spline of order  k , and below it its first
% and second (piecewise) derivative, to illustrate the following facts (try
% out the BSPLIGUI if you want to experiment with examples of your own).
%
% 1. The B-spline  B_{j,k} = B( . | t_j, ..., t_{j+k})  is pp of order k
% with breaks at t_j, ..., t_{j+k}  (and nowhere else). Actually, its
% nontrivial polynomial pieces are all of exact degree  k-1 .
%
% E.g., the rightmost B-spline involves  5  knots, hence is of order  4 ,
% i.e., a CUBIC B-spline.  Correspondingly, its second derivative is
% piecewise linear.

cl = ['g','r','b','k','k'];
v = 5.4; d1 = 2.5; d2 = 0; s1 = 1; s2 = .5;
t1 = [0 .8 2];
t2 = [3 4.4 5  6];
t3 = [7  7.9  9.2 10 11];
tt = [t1 t2 t3];
ext = tt([1 end])+[-.5 .5];
plot(ext([1 2]),[v v],cl(5))
hold on
plot(ext([1 2]),[d1 d1],cl(5))
plot(ext([1 2]),[d2 d2],cl(5))
tmp=NaN;
ts = [tt;tt;repmat(NaN,size(tt))];
ty = repmat(.2*[-1;0;NaN],size(tt));
plot(ts(:),ty(:)+v,cl(5))
plot(ts(:),ty(:)+d1,cl(5))
plot(ts(:),ty(:)+d2,cl(5))
b1 = spmak(t1,1);
p1 = [t1;0 1 0];
db1 = fnder(b1);
p11 = fnplt(db1,'j');
p12 = fnplt(fnder(db1));
lw = 2;
plot(p1(1,:),p1(2,:)+v,cl(2),'linew',lw)
plot(p11(1,:),s1*p11(2,:)+d1,cl(2),'linew',lw)
plot(p12(1,:),s2*p12(2,:)+d2,cl(2),'linew',lw)
axis([-1 12 -2 6.5])
b1 = spmak(t2,1);
p1 = fnplt(b1);
db1 = fnder(b1);
p11 = [t2;fnval(db1,t2)];
p12=fnplt(fnder(db1),'j');
plot(p1(1,:),p1(2,:)+v,cl(3),'linew',lw)
plot(p11(1,:),s1*p11(2,:)+d1,cl(3),'linew',lw)
plot(p12(1,:),s2*p12(2,:)+d2,cl(3),'linew',lw)
b1 = spmak(t3,1);
p1 = fnplt(b1);
db1 = fnder(b1);
p11 = fnplt(db1);
p12=[t3;fnval(fnder(db1),t3)];
plot(p1(1,:),p1(2,:)+v,cl(4),'linew',lw)
plot(p11(1,:),s1*p11(2,:)+d1,cl(4),'linew',lw)
plot(p12(1,:),s2*p12(2,:)+d2,cl(4),'linew',lw)
tey = v+1.5;
text(t1(2)-.5,tey,'linear','color',cl(2))
text(t2(2)-.8,tey,'quadratic','color',cl(3))
text(t3(3)-.5,tey,'cubic','color',cl(4))
text(-2,v,'B'), text(-2,d1,'DB'), text(-2,d2,'D^2B')
axis off, hold off

%%
% 2. B_{j,k} is positive on the interval  (t_j .. t_{j+k})  and is zero off
% that interval. It also vanishes at the endpoints of that interval (unless
% the endpoint is a knot of multiplicity  k ; see the rightmost example on
% the next slide).

%%
% 3. Knot   MULTIPLICITY   determines the smoothness with which the two
% adjacent polynomials join across that knot. In shorthand, the rule is:
%
%         knot multiplicity  +  number of smoothness conditions  =  order
%
% Above you see cubic B-splines and, below them, their first two
% derivatives, with a certain knot of multiplicity 1, 2, 3, 4 (as indicated
% by the length of the knot line).
%
% E.g., since the order is  4 , a double knot means just 2 smoothness conditions,
% i.e., just continuity  across that knot of the function and its first derivative.

d2 = -1;
t1=[7  7.9  9.2 10 11]-7;
t2=[7  7.9 7.9  9 10]-2;
t3=[7  7.9 7.9 7.9 9]+2;
t4=[7  7.9 7.9 7.9 7.9]+5;
plot(1,1,'or')
[m,tt]=knt2mlt([t1 t2 t3 t4]);
ext = tt([1 end])+[-.5 .5];
plot(ext,[v v],cl(5)), hold on
plot(ext,[d1 d1],cl(5))
plot(ext,[d2 d2],cl(5))
ts = [tt;tt;repmat(NaN,size(tt))];
ty = .2*[-m-1;zeros(size(m));m];
plot(ts(:),ty(:)+v,cl(5))
plot(ts(:),ty(:)+d1,cl(5))
plot(ts(:),ty(:)+d2,cl(5))
b1 = spmak(t1,1);
p1=fnplt(b1);
db1=fnder(b1);
p11=fnplt(db1);
p12=[t1;fnval(fnder(db1),t1)];
plot(p1(1,:),p1(2,:)+v,cl(1),'linew',lw)
plot(p11(1,:),s1*p11(2,:)+d1,cl(1),'linew',lw)
plot(p12(1,:),s2*p12(2,:)+d2,cl(1),'linew',lw)
text(-2,v,'B'), text(-2,d1,'DB'), text(-2,d2,'D^2B')
b1 = spmak(t2,1);
p1=fnplt(b1);
db1=fnder(b1);
p11=fnplt(db1);
p12=fnplt(fnder(db1),'j');
plot(p1(1,:),p1(2,:)+v,cl(2),'linew',lw)
plot(p11(1,:),s1*p11(2,:)+d1,cl(2),'linew',lw)
plot(p12(1,:),s2*s2*p12(2,:)+d2,cl(2),'linew',lw)
b1 = spmak(t3,1);
p1=fnplt(b1);
db1=fnder(b1);
p11=fnplt(db1,'j');
p12=fnplt(fnder(db1),'j');
plot(p1(1,:),p1(2,:)+v,cl(3),'linew',lw)
plot(p11(1,:),s1*s2*p11(2,:)+d1,cl(3),'linew',lw)
plot(p12(1,:),s2*s2*p12(2,:)+d2,cl(3),'linew',lw)
b1 = spmak(t4,1);
p1=fnplt(b1);
db1=fnder(b1);
p11=fnplt(db1);
p12=fnplt(fnder(db1));
plot(p1(1,:),p1(2,:)+v,cl(4),'linew',lw)
plot(p11(1,:),s2*p11(2,:)+d1,cl(4),'linew',lw)
plot(p12(1,:),s2*s2*p12(2,:)+d2,cl(4),'linew',lw)
text(t2(2)-.5,tey,'2-fold','color',cl(2))
text(t3(2)-.8,tey,'3-fold','color',cl(3))
text(t4(3)-.8,tey,'4-fold','color',cl(4))
axis off, hold off

%% REFINEMENT and KNOT INSERTION
% Any B-form can be REFINED, i.e., converted, by KNOT INSERTION, into the
% B-form for the same function, but for a finer knot sequence. The finer
% the knot sequence, the closer is the control polygon to the function
% being represented.
%
% The refined control polygon shown above in red is obtained for the cubic
% spline  SP  of a previous picture by the following commands:
%
%  >> b = knt2brk(fnbrk(sp,'k')); spref = fnrfn( sp, (b(2:end)+b(1:end-1))/2);
%  >> cr = fnbrk(spref,'c'); plot( aveknt( fnbrk(spref,'knots'), 4), cr,':r*')

sp = savesp;
fnplt(sp,2.5);
hold on
c = fnbrk(sp,'c');
plot(aveknt(fnbrk(sp,'k'),4),c,':ok')
b=knt2brk(fnbrk(sp,'k'));
spref = fnrfn(sp,(b(2:end)+b(1:end-1))/2);
cr = fnbrk(spref,'c');
plot(aveknt(fnbrk(spref,'knots'),4),cr,':*r')
axis([-.72 7 -.5 3.5])
axis off
hold off

%% REFINEMENT and KNOT INSERTION (cont.)
% In this second example, we start with the vertices of the standard
% diamond as our control points, but run through the sequence twice:
%
%  >> ozmz = [1 0 -1 0]; c = [ozmz ozmz 1; 0 ozmz ozmz]; circle = spmak(-4:8,c);
%  >> fnplt(circle), hold on, plot(c(1,:),c(2,:),':ok'), hold off
%
% However, when we plot the resulting spline, we get a curve that begins
% and ends at the origin,  due to the fact that we chose to make the knot
% sequence simple, hence our spline vanishes at the endpoints of its basic
% interval, [-4 .. 8] . We really only want the part corresponding to the
% interval [0 .. 4], as can be seen by plotting just that part, more
% boldly:   >> fnplt(circle,[0 4],2)

ozmz = [1 0 -1 0];
c = [ozmz ozmz 1; 0 ozmz ozmz];
circle = spmak(-4:8,c);
fnplt(circle)
hold on
plot(c(1,:),c(2,:),':ok')
axis(1.1*[-1 1 -1 1])
axis equal
axis off
fnplt(circle,[0 4],2)
hold off

%% REFINEMENT and KNOT INSERTION (cont.)
% To get just the circle, we restrict our spline to the interval [0 .. 4].
% We do this by converting to ppform, restricting to [0 .. 4], then
% converting to B-form.
%
%  >> circ = fn2fm( fnbrk( fn2fm(circle,'pp'), [0 4]), 'B-'); fnplt(circ,2.5)
%  >> cc = fnbrk(circ,'c');  hold on, plot(cc(1,:),cc(2,:),':ok')
%
% Refinement of the resulting knot sequence leads to a control polygon much
% closer to the circle:
%
%  >> ccc = fnbrk(fnrfn(circ,0.5:4),'c');   plot(ccc(1,:),ccc(2,:),'-r*'), hold off

circ = fn2fm(fnbrk(fn2fm(circle,'pp'),[0 4]),'B-');
fnplt(circ,2.5)
hold on
axis(1.1*[-1 1 -1 1])
axis equal
axis off
cc = fnbrk(circ,'c');
plot(cc(1,:),cc(2,:),':ok')
ccc = fnbrk(fnrfn(circ,.5:4),'c');
plot(ccc(1,:),ccc(2,:),'-r*')
hold off

%% MULTIVARIATE SPLINES
% A spline in this toolbox can also be multivariate, namely the tensor
% product of univariate splines. The B-form for such a function is only
% slightly more complicated, with the knots now a cell array containing the
% various univariate knot sequences, and the coefficient array suitably
% multidimensional.
%
% For example, the above random spline surface results from the following:
%
%  >> fnplt( spmak( {augknt(0:4,4),augknt(0:4,3)}, rand(7,8) ) )
%
% This spline is cubic in the first variable (there are 11 knots and 7
% coefficients in that variable), but only piecewise constant in the second
% variable ((2+5+2)-8 = 1).

fnplt( spmak( {augknt(0:4,4),augknt(0:4,3)}, rand(7,8) ) )


displayEndOfDemoMessage(mfilename)

##### SOURCE END #####
-->
   </body>
</html>