www.gusucode.com > curvefit 案例源码程序 matlab代码 > curvefit/spcrvdem.m
%% Splines in the Plane % % This example shows how to use the |spmak|, |spcrv|, |cscvn| and |rscvn| % commands from Curve Fitting Toolbox(TM) to construct spline curves in the % plane. This includes plotting tangents and computing the area enclosed by a % curve. % Copyright 1987-2012 The MathWorks, Inc. %% A Simple Spline Curve % Curve Fitting Toolbox can handle _vector-valued_ splines. A % d-vector-valued univariate spline provides a curve in d-space. In this % mode, |d = 2| is most common, as it gives plane curves. % % Here is an example, in which a spline with 2-dimensional coefficients is % constructed and plotted. knots = [1,1:9,9]; curve = spmak( knots, repmat([ 0 0; 1 0; 1 1; 0 1 ], 2,1).' ); t = linspace(2,8,121); values = fnval(curve,t); plot(values(1,:),values(2,:),'LineWidth',2); axis([-.4 1.4 -.2 1.2]), axis equal title('A Spline Curve'); hold on %% A Word of Caution % You may have noticed that this example did not use |fnplt| to plot the % curve, but instead plotted some points on the curve obtained by |fnval|. % Here is the code again: % % t = linspace(2,8,121); % values = fnval(curve,t); % plot(values(1,:),values(2,:),'LineWidth',2) % % Using |fnplt| directly with this particular spline curve gives the red curve % in the figure below. fnplt(curve,'r',.5); title('The Full Spline Curve, in Red') %% % The explanation? % % The spline is of order 4, yet the end knots in the knot sequence knots %% % only have multiplicity 2. Therefore, all the B-splines of order 4 for this % knot sequence are 0 at the endpoints of the basic interval. This makes the % curve start and stop at (0,0). %% A Remedy % Since, in this case, we are really interested only in the curve segment % corresponding to the parameter interval [3 .. 7], we can use |fnbrk| to % extract that part, and then have no difficulty plotting it, in yellow, with % |fnplt|. mycv = fnbrk(curve,[3 7]); fnplt(mycv,'y',2.5); title('The Spline Curve of Interest, in Yellow') %% The Area Enclosed By This Curve % Since you now have a spline, namely |mycv|, that describes the curve % (and nothing else), you can easily compute the area enclosed by this % closed curve, as follows. area = diff(fnval(fnint( ... fncmb(fncmb(mycv,[0 1]),'*',fnder(fncmb(mycv,[1 0]))) ... ),fnbrk(mycv,'interval'))) %% % With a little effort, you can recognize this as the value of the % integral % % int y(t) d(x(t)) = int y(t) Dx(t) dt % % over the basic interval of the spline |mycv|, with |(x(t),y(t)) := % fnval(mycv,t)| the point on the curve corresponding to the parameter value % |t|. Here, |fncmb(mycv,[1,0])|, |fncmb(mycv,[0,1])| describe the two % components of the spline curve, i.e., the scalar-valued splines |x| and |y|. % % Also, the curve is roughly a circle with radius 1/2. Hence, you would % expect an area of, roughly, disp(pi/4) %% % But why is the computed area _negative_? Because the area enclosed by the % curve lies to the left as one travels on the curve with increasing |t|. % To verify this, we draw some tangent vectors. %% Add Some Tangent Vectors % We redraw the curve and also draw the tangent vector to the curve at % some points. hold off fnplt(mycv,'y',2.5); hold on t = 3:.4:6.2; cv = fnval(curve, t); cdv = fnval(fnder(curve), t); quiver(cv(1,:),cv(2,:), cdv(1,:),cdv(2,:)); title('A Spline Curve With Some Tangents') axis([-.4 1.4 -.2 1.2]), axis equal %% The Intersection of the Curve With a Straight Line % If you wanted to determine the points of intersection of this spline % curve with the straight line |y = x|, the following code would give them % to you, and plot the segment of that straight line inside the curve: cuts = fnval(mycv, ... mean(fnzeros(fncmb(fncmb(mycv,[0,1]),'-',fncmb(mycv,[1,0]))))); plot(cuts(1,:), cuts(2,:),'y','LineWidth',2.5) hold off title('A Spline Curve With Some Tangents and a Cut Across') %% SPCRV: The Control Polygon and the Corresponding Spline Curve % Spline curves are used extensively in the generation of illustrations in % which nothing more than a smooth curve of a certain roughly imagined % shape is required. For this, Curve Fitting Toolbox contains a special % command, |spcrv|, which can be used independently of the rest of the % toolbox. % % Given a sequence of points in the plane and, optionally, an order |k|, % |spcrv| generates, by repeated midpoint knot insertion, the spline curve of % order |k| whose control polygon is specified by the given sequence. % % The figure below shows such a control polygon, and the corresponding spline % curve of order 3. points = [0 0; 1 0; 1 1; 0 2; -1 1; -1 0; 0 -1; 0 -2].'; values = spcrv(points,3); plot(points(1,:),points(2,:),'k'); axis([-2 2.25 -2.1 2.2]); hold on plot(values(1,:),values(2,:),'r','LineWidth',1.5); legend({'Control Polygon' 'Quadratic Spline Curve'}, 'location','SE'); %% % Notice that the curve touches each segment of the control polygon at its % midpoint, and follows the shape outlined by the control polygon. %% Raising the Order % Raising the order |k| will pull the curve away from the control polygon and % make it smoother, but also shorter. Here, we have added the corresponding % spline curve of order 4. value4 = spcrv(points,4); plot(value4(1,:),value4(2,:),'b','LineWidth',2); legend({'Control Polygon' 'Quadratic Spline Curve' ... 'Cubic Spline Curve'}, 'location','SE'); %% CSCVN % On the other hand, to obtain an interpolating curve, you could use the % |cscvn| command, which provides a parametric `natural' cubic spline curve. fnplt(cscvn(points), 'g',1.5); legend({'Control Polygon' 'Quadratic Spline Curve' ... 'Cubic Spline Curve' 'Interpolating Spline Curve'}, ... 'location','SE'); %% % By adding the point (.95,-.05) near the second control point, (1,0), we % can create an interpolating spline curve that turns faster there. np = size(points, 2); fnplt( cscvn([ points(:,1) [.95; -.05] points(:,2:np) ]), 'm',1.5); plot(.95,-.05,'*'); legend({'Control Polygon' 'Quadratic Spline Curve' ... 'Cubic Spline Curve' 'Interpolating Spline Curve' ... 'Faster Turning Near (1,0)'}, ... 'location','SE'); hold off %% RSCVN % You can also obtain a tangent-continuous curve composed of circular arcs % that passes through a given sequence of points in the plane and, % optionally, is orthogonal to given normal directions at the points. The % command |rscvn| provides such a curve. % % For example, the following generates a circle c = rscvn([-1 1 -1;0 0 0],[1 1;0 0]); %% % as its plot shows. fnplt(c); axis([-1.05 1.05 -1.05 1.05]), axis equal, axis off %% % |c| is a quadratic rational spline consisting of just two pieces, % as the following commands make clear. [form, order, breaks] = fnbrk(c,'f','o','b') %% % It is easy to generate striking patterns with this tool using just a few % data points. For example, here is a version of the design on the % Bronze Triskele Medallion in the Ulster Museum in Belfast, supposedly % done by pieces of circular arcs a long time ago. pp =[zeros(1,7); 5.4, 3, 6.9, 2.75, 2.5, .5, 5]; alpha = 2*pi/3; ca = cos(alpha); sa = sin(alpha); c = [ca sa;-sa ca]; d = [0 0 .05 -.05;1 -1 .98 .98]; d = [d c*d]; yin = rscvn([pp(:,[7,1:3]),c*pp(:,3:4),pp(:,3)], d(:,[1 2 1 4 7 5 1])); fnplt(yin), hold on, fnplt(fncmb(yin,c)), fnplt(fncmb(yin,c')) yang = rscvn([pp(:,6),-pp(:,6),pp(:,5),c*pp(:,4)],[d(:,[2 1 1]),c(:,2)]); fnplt(yang), fnplt(fncmb(yang,c)), fnplt(fncmb(yang,c')) axis([-7.2 7.2 -7.2 7.2]), axis equal, axis off, hold off displayEndOfDemoMessage(mfilename)