function SamplPs =CollectSPFergusonCir(N)%N表示样点个数
step =pi/N;%为编程方便,把N作为间隔数
xita =0:step:pi;%等间距取角度向量%根据参数方程计算圆周上的点SamplPs(1,:)=cos(xita);%x的坐标放在矩阵的第一行SamplPs(2,:)=sin(xita);%y的坐标放在矩阵的第二行%根据参数方程计算圆周上的点% plot(SamplPs(1,:),SamplPs(2,:),'k*')%测试程序正确性用的% axis equal;
function ts =ThreeTangEq(SamplPs,t0,tn)%Sampls输入型值点;t0,tn是首末点的切矢量
N =length(SamplPs(1,:));%SamPls(1,:)表示取矩阵Sampls第一行的所有元素%注意矩阵一行或者一列的全部元素形成一个向量%length()就是测一个向量的长度,也就是元素个数%注意矩阵SamPs存储采集的型值点,所有型值点的x坐标放在第一行%参见本文件夹中的函数SamplPs = CollectSPFergusonCir(N)%因此N = length(SamPls(1,:))得到的就是本程序使用的型值点的个数
CoefA =eye(N,N);%建立一个N*N的单位矩阵%对这个矩阵进行更改后就是系数矩阵%更改第2行到第N-1行fori=2:N-1CoefA(i,i-1)=1;CoefA(i,i)=4;CoefA(i,i+1)=1;end
EqB2D =zeros(N,2);%注意每一行表示一个向量,共有N行EqB2D(1,:)= t0;fori=2:N-1EqB2D(i,:)=3*(SamplPs(:,i+1)'-SamplPs(:,i-1)');%V'表示取向量V的转置;SamPls(:,i+1)'就是把列向量表示的点变为行向量表示的点endEqB2D(N,:)= tn;
ts=CoefA\EqB2D;%求解线性方程组得到每个点的切矢量
function DensePs =FergCurvSegDenPs(r0,r1,r0c,r1c)%为绘制参数三次样条曲线段采集密化点%r0,r1,r0c,r1c分别表示首末端点的位置和切矢量%采集绘制曲线段用的密化点
k=1;for u=0:1/100:1DensePs(k,1)=[1 u u*u u*u*u]*[1000;0010;-33-2-1;2-211]*[r0(1)r1(1)r0c(1)r1c(1)]';%参数三次样条曲线段的计算,参见公式(4.3')DensePs(k,2)=[1 u u*u u*u*u]*[1000;0010;-33-2-1;2-211]*[r0(2)r1(2)r0c(2)r1c(2)]';%参数三次样条曲线段的计算,参见公式(4.3')
k=k+1;end%采集绘制曲线段用的密化点
functionRenderFergusonCurve(SamplPs,SamplPTs)%输入型值点SamplPs和型值点处的切矢量SamplPTs%step1 汇总各个曲线段上的密化点
N =length(SamplPs(1,:));%得到密化点数目,参见Test3TEq
N = N-1;%每两个点之间有一个曲线段,这就是曲线段的数目
DensePs =[];%定义空矩阵fori=1:N
r0 =SamplPs(:,i)';r1 = SamplPs(:,i+1)';
r0c =SamplPTs(i,:);r1c =SamplPTs(i+1,:);
Ps =FergCurvSegDenPs(r0,r1,r0c,r1c);%在当前曲线段上构造密化点
DensePs =[DensePs;Ps];%把当前曲线段上采集的密化点数组放入全局数组endplot(DensePs(:,1),DensePs(:,2),'linewidth',2)
hold on
plot(SamplPs(1,:),SamplPs(2,:),'k*')%添加型值点便于观察
N=N+1;fori=1:N%逐个绘制切线段
x =[SamplPs(1,i),SamplPs(1,i)+SamplPTs(i,1)];
y =[SamplPs(2,i),SamplPs(2,i)+SamplPTs(i,2)];%每个切线段的起始点是[SamplPs(1,i),SamplPs(2,i)]%终了点是[SamplPs(1,i) + SamplPTs(i,1),SamplPs(2,i) + SamplPTs(i,2)]plot(x,y,'r','linewidth',1);%把切线绘制为红色,宽度为2end
hold off
axis equal