📄 spalldem.html
字号:
<html xmlns:mwsh="http://www.mathworks.com/namespace/mcode/v1/syntaxhighlight.dtd"> <head> <meta http-equiv="Content-Type" content="text/html; charset=utf-8"> <!--This HTML is auto-generated from an M-file.To make changes, update the M-file and republish this document. --> <title>Intro to B-form</title> <meta name="generator" content="MATLAB 7.1"> <meta name="date" content="2005-07-27"> <meta name="m-file" content="spalldem"> <link rel="stylesheet" type="text/css" href="../../matlab/demos/private/style.css"> </head> <body> <div class="header"> <div class="left"><a href="matlab:edit spalldem">Open spalldem.m in the Editor</a></div> <div class="right"><a href="matlab:echodemo spalldem">Run in the Command Window</a></div> </div> <div class="content"> <h1>Intro to B-form</h1> <introduction> <p>A quick introduction to the B-form.</p> </introduction> <h2>Contents</h2> <div> <ul> <li><a href="#3">Local partition of unity and convex hull property</a></li> <li><a href="#4">Convex hull property (cont.) and control polygon</a></li> <li><a href="#5">CONTROL POLYGON (cont.)</a></li> <li><a href="#11">REFINEMENT and KNOT INSERTION</a></li> <li><a href="#12">REFINEMENT and KNOT INSERTION (cont.)</a></li> <li><a href="#13">REFINEMENT and KNOT INSERTION (cont.)</a></li> <li><a href="#14">MULTIVARIATE SPLINES</a></li> </ul> </div> <p>In this toolbox, a pp function in B-form is often called a SPLINE.</p> <p>The B-FORM of a (univariate) pp is specified by its (nondecreasing) k n o t sequence t and by its B-spline c o e f f i c i e n t sequence a . </p> <p>SPMAK(t,a) returns the corresponding B-form, for use in FN... commands.</p> <p>The resulting spline is of o r d e r k := LENGTH(t) - SIZE(a,2). This means that all its polynomial pieces have degree < k. </p><pre class="codeinput">t = [.1 .4 .5 .8 .9];fnplt(spmak(t,1),2.5);tmp = repmat(t,3,1);ty = repmat(.1*[1;-1;NaN],1,5);hold <span class="string">on</span>plot(tmp(:),ty(:))text(.65,.5,<span class="string">'B( \cdot | .1, .4, .5, .8, .9)'</span>,<span class="string">'Fontsize'</span>,10)text(.02,1.,<span class="string">'s(x) = \Sigma_j B( x | t_j , \ldots, t_{j+k} ) a(:,j)'</span>,<span class="string">'Fontsize'</span>,18,<span class="string">'color'</span>,<span class="string">'r'</span>)axis([0 1 -.2 1.2])hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="spalldem_01.png"> <p>To say that s has knots t and coefficients a means that</p><pre> ns(x) := sum B_{j,k}(x) * a(:,j) , with B_(j,k) = B( . | t_j, ..., t_{j+k} ) j = 1</pre><p>the j-th B-SPLINE of ORDER k for the given KNOT SEQUENCE t ,</p> <p>i.e., the B-spline with knots t_j, ..., t_{j+k} ; see the sample B-spline above.</p> <h2>Local partition of unity and convex hull property<a name="3"></a></h2> <p>The value of the spline s(x) = sum_j B_{j,k}(x) a(:,j) at any x in the knot interval [t_i .. t_{i+1}] is a CONVEX combination of the k coefficients a(:.i-k+1), ..., a(:,i) since, on that interval, only the k B-splines B_{i-k+1,k}, ..., B_{i,k} are nonzero, and they are nonnegative there and sum to 1. </p> <p>This is often summarized by saying that the B-splines provide a LOCAL (nonnegative) PARTITION of UNITY .</p><pre class="codeinput">k = 3;n = 3;t = [1 1.7 3.2 4.2 4.8 6];tt = (10:60)/10;vals = fnval(spmak(t,eye(k)),tt);plot(tt.',vals.');hold <span class="string">on</span>ind = find(tt>=t(3)&tt<=t(4));plot(tt(ind).',vals(:,ind).',<span class="string">'linewi'</span>,3)plot(t([3 4]),[1 1],<span class="string">'w'</span>,<span class="string">'linewi'</span>,3)axis([-.5 7 -.5 1.5])ty = repmat(.1*[1;-1;NaN],1,6);plot([0 0 -.2 0 0 -.2 0 0],[-.5 0 0 0 1 1 1 1.5],<span class="string">'k'</span>)text(-.5,0,<span class="string">'0'</span>)text(-.5,1,<span class="string">'1'</span>)tmp = repmat(t,3,1);plot(tmp(:),ty(:),<span class="string">'k'</span>);yd = -.25; text(t(1),yd,<span class="string">'t_{i-2}'</span>);text(t(3),yd,<span class="string">'t_i'</span>);text(t(4),yd,<span class="string">'t_{i+1}'</span>);text(1.8,.5,<span class="string">'B_{i-2,3}'</span>);text(5,.45,<span class="string">'B_{i,3}'</span>);axis <span class="string">off</span>hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="spalldem_02.png"> <h2>Convex hull property (cont.) and control polygon<a name="4"></a></h2> <p>When the coefficients are points in the plane and, correspondingly, the spline s(x) = sum_j B_{j,k}(x) a(:,j) traces out a curve, this means that the curve piece { s(x) : t_i <= x <= t_{i+1} } lies in the convex hull (shown in yellow above) of the k points a(:,i-k+1), ... a(:,i). </p> <p>For a quadratic spline (i.e., k = 3), as shown here, it even means that the curve is tangent to the CONTROL POLYGON (shown in black dashes). This is the broken line that connects the coefficients (which are called CONTROL POINTS in this connection and shown here as open circles). </p><pre class="codeinput">t = 1:9;c = [2 1.4;1 .5; 2 -.4; 5 1.4; 6 .5; 5 -.4].';sp = spmak(t,c);fill(c(1,3:5),c(2,3:5),<span class="string">'y'</span>,<span class="string">'edgecolor'</span>,<span class="string">'y'</span>);hold <span class="string">on</span>fnplt(sp,t([3 7]),1.5)fnplt(sp,t([5 6]),3)plot(c(1,:),c(2,:),<span class="string">':ok'</span>)axis([.5 7 -.8 1.8])axis <span class="string">off</span>text(2,-.55,<span class="string">'a(:,i-2)'</span>)text(5,1.6,<span class="string">'a(:,i-1)'</span>)text(6.1,.5,<span class="string">'a(:,i)'</span>)hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="spalldem_03.png"> <h2>CONTROL POLYGON (cont.)<a name="5"></a></h2> <p>We can think of the graph of the SCALAR-valued spline s = sum_j B_{j,k}*a(j) as the curve x |--> (x,s(x)) . Since x = sum_j B_{j,k}(x) t^*_j for x in the interval [t_k .. t_{n+1}] and with t^*_i : = (t_{i+1} + ... + t_{i+k-1})/(k-1), i = 1:n, the knot averages obtainable by AVEKNT(t,k) , the CONTROL POLYGON for a scalar-valued spline is the broken line with vertices (t^*_i, a(i)), i=1:n. </p> <p>The example shows a CUBIC spline (k = 4), with 4-fold end knots, hence t^*_1 = t_1 and t^*_n = t_{n+k} .</p><pre class="codeinput">t = [0 .2 .35 .47 .61 .84 1]*(2*pi);s=t([1 3 4 5 7]);knots = augknt(s,4);sp = spapi(knots,t,sin(t)+1.8);fnplt(sp,2);hold <span class="string">on</span>c = fnbrk(sp,<span class="string">'c'</span>);ts = aveknt(knots,4);plot(ts,c,<span class="string">':ok'</span>)tt = [s;s;repmat(NaN,size(s))];ty = repmat(.25*[-1;1;NaN], size(s));plot(tt(:),ty(:),<span class="string">'r'</span>)plot(ts(1,:),zeros(size(ts)),<span class="string">'*'</span>)text(knots(5),-.5,<span class="string">'t_5'</span>)text(ts(2),-.45,<span class="string">'t^*_2'</span>)text(knots(1)-.28,-.5,<span class="string">'t_1=t_4'</span>)text(knots(end)-.65,-.45,<span class="string">'t_{n+1}=t^*_n=t_{n+4}'</span>)axis([-.72 7 -.5 3.5])axis <span class="string">off</span>hold <span class="string">off</span>savesp = sp;</pre><img vspace="5" hspace="5" src="spalldem_04.png"> <p>The essential parts of the B-form are the knot sequence t and the B-spline coefficient sequence a . Other parts are the NUMBER n of the B-splines or coefficients involved , the ORDER k of its polynomial pieces, and the DIMENSION d of its coefficients. In particular, size(a) equals [d,n] . </p> <p>There is one more part, namely the BASIC INTERVAL , [ t(1) .. t(end) ]. It is used as the default interval when plotting the function. Also, while the spline is taken to be continuous from the right, it is made continuous from the left at the right endpoint of the basic interval, as illustrated above, where values of the spline made by </p><pre>>> sp = spmak(augknt(0:3,3),[-1,0,1,0,-1]);</pre><p>are plotted.</p><pre class="codeinput">b = 0:3;sp = spmak(augknt(b,3),[-1,0,1,0,-1]);x = linspace(-1,4,51);plot(x,fnval(sp,x),<span class="string">'x'</span>)hold <span class="string">on</span>axis([-2 5,-1.5,1])tx = repmat(b,3,1);ty = repmat(.15*[1;-1;NaN],1,length(b));plot(tx(:),ty(:),<span class="string">'-r'</span>)hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="spalldem_05.png"> <p>FNBRK can be used to obtain any or all parts of a B-form. For example, here is the output provided by FNBRK(SP), with SP the B-form of the spline shown above: </p><pre>knots(1:n+k) 0 0 0 1 2 3 3 3coefficients(d,n) -1 0 1 0 -1number n of coefficients 5order k 3dimension d of target 1</pre><p>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. </p> <p>The next slides give more detailed information about the B-splines, in particular about the important role played by knot MULTIPLICITY . </p><pre class="codeinput">b = 0:3;sp = spmak(augknt(b,3),[-1,0,1,0,-1]);x = linspace(-1,4,51);plot(x,fnval(sp,x),<span class="string">'x'</span>),hold <span class="string">on</span>,axis([-2 5,-1.5,1]),tx = repmat(b,3,1);ty = repmat(.15*[1;-1;NaN],1,length(b));plot(tx(:),ty(:),<span class="string">'-r'</span>)hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="spalldem_06.png"> <p>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). </p> <p>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 . </p> <p>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. </p><pre class="codeinput">cl = [<span class="string">'g'</span>,<span class="string">'r'</span>,<span class="string">'b'</span>,<span class="string">'k'</span>,<span class="string">'k'</span>];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 <span class="string">on</span>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,<span class="string">'j'</span>);p12 = fnplt(fnder(db1));lw = 2;plot(p1(1,:),p1(2,:)+v,cl(2),<span class="string">'linew'</span>,lw)plot(p11(1,:),s1*p11(2,:)+d1,cl(2),<span class="string">'linew'</span>,lw)plot(p12(1,:),s2*p12(2,:)+d2,cl(2),<span class="string">'linew'</span>,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),<span class="string">'j'</span>);plot(p1(1,:),p1(2,:)+v,cl(3),<span class="string">'linew'</span>,lw)plot(p11(1,:),s1*p11(2,:)+d1,cl(3),<span class="string">'linew'</span>,lw)plot(p12(1,:),s2*p12(2,:)+d2,cl(3),<span class="string">'linew'</span>,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),<span class="string">'linew'</span>,lw)plot(p11(1,:),s1*p11(2,:)+d1,cl(4),<span class="string">'linew'</span>,lw)plot(p12(1,:),s2*p12(2,:)+d2,cl(4),<span class="string">'linew'</span>,lw)tey = v+1.5;text(t1(2)-.5,tey,<span class="string">'linear'</span>,<span class="string">'color'</span>,cl(2))text(t2(2)-.8,tey,<span class="string">'quadratic'</span>,<span class="string">'color'</span>,cl(3))text(t3(3)-.5,tey,<span class="string">'cubic'</span>,<span class="string">'color'</span>,cl(4))text(-2,v,<span class="string">'B'</span>), text(-2,d1,<span class="string">'DB'</span>), text(-2,d2,<span class="string">'D^2B'</span>)axis <span class="string">off</span>, hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="spalldem_07.png"> <p>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). </p> <p>3. Knot MULTIPLICITY determines the smoothness with which the two adjacent polynomials join across that knot. In shorthand, the rule is: </p><pre> knot multiplicity + number of smoothness conditions = order</pre><p>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). </p> <p>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. </p><pre class="codeinput">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,<span class="string">'or'</span>)[m,tt]=knt2mlt([t1 t2 t3 t4]);ext = tt([1 end])+[-.5 .5];plot(ext,[v v],cl(5)), hold <span class="string">on</span>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),<span class="string">'linew'</span>,lw)⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -