MATLAB绘制矩阵权（Matrix weighted）有理Bezier曲线

是根据浙江大学杨勋年老师15年发表的文章Matrix weighted rational curves and surfaces写的博文，这一篇主要是绘制矩阵权（Matrix weighted）有理Bezier曲线。

在上一篇博文讲述了有理Bezier曲线之后，杨老师对其进行了扩展，并赋予其几何意义，将法向信息加以考虑得到矩阵权的有理Bezier曲线，定义为

Q(t)=[∑i=0nMiBi,n]−1∑i=0nMiPiBi,n(t), t∈[0,1]

其中Bi,n(t)=n!i!(n−i)!ti(1−t)n−i是Bernstein基

存在很多种定义矩阵权的方法，这里采用一种几何上比较直观的方法定义。假设对于ni,i=0,1,⋯,n是一族单位法向量并且ωi>0,μi>−1,i=0,1,⋯,n是一族实数，对于Q(ξ)矩阵权定义如下

Mi=ωi(I+μininTi),i=0,1,⋯,n

其中I是d×d的单位阵

这样就可以进行矩阵权有理Bezier曲线的绘制了

这里对文中的图1（a）还原，这里所有的ωi都是1，除去μ2,μ5,μ8取20之外所有的μi都是1



整个的MATLAB代码如下

function  normal_bezier
%NORMAL_BEZIER 此处显示有关此函数的摘要
%   此处显示详细说明
N=[2 0.8;1 0;7 -1;10 1;0 1;0 1;0 1;-10 1;-7 -1;-1 0;-2 0.8]';
vertices = [-1 6;-0.75 5;-1.2 2.5;-1.6 1;-1 0;0 0;1 0;1.6 1;1.2 2.5;0.75 5;1 6]';
NumPoint=size(vertices,2)-1;%点的个数
M=cell(1,NumPoint+1);

for i=1:NumPoint+1
niu(i)=1;
omega(i)=1;
end
niu(3)=20;niu(6)=20;niu(9)=20;

for i=1:NumPoint+1;
S=sqrt(N(1,i)^2+N(2,i)^2);
N(1,i)=N(1,i)/S;
N(2,i)=N(2,i)/S;
end

I=[1 0;0 1];
for i=1:NumPoint+1
M{i}=omega(i)*(I+niu(i)*N(:,i)*N(:,i)');
end

xx=[];
yy=[];

for t=0.001:0.001:0.999
Mat=[0;0];
invMat=[0 0;0 0];
for j=0:NumPoint
w=factorial(NumPoint)/(factorial(j)*factorial(NumPoint-j))*(1-t)^(NumPoint-j)*t^(j);
invMat=invMat+M{j+1}*w;
Mat=Mat+M{j+1}*vertices(:,j+1)*w;
end
newPoint = invMat\Mat;
xx=[xx newPoint(1,1)];
yy=[yy newPoint(2,1)];
end

plot(vertices(1,:),vertices(2,:),'b');hold on;
for i=1:NumPoint+1
quiver(vertices(1,i),vertices(2,i),N(1,i),N(2,i),'r','filled','LineWidth',2);
end
hold on;grid on;
axis tight;
axis equal;
xlabel('X');ylabel('Y');
plot(xx,yy,'r');

end

整体效果还是不错。再接再厉^^