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

是根据浙江大学杨勋年老师15年发表的文章Matrix weighted rational curves and surfaces写的博文，这一篇主要是绘制矩阵权（Matrix weighted）有理Bezier曲线。
在上一篇博文讲述了有理Bezier曲线之后，杨老师对其进行了扩展，并赋予其几何意义，将法向信息加以考虑得到矩阵权的有理Bezier曲线，定义为

$$Q(t)=[\sum_{i=0}^n M_iB_{i,n}]^{-1}\sum_{i=0}^nM_iP_iB_{i,n}(t),\ t\in[0,1]$$
其中 $B_{i,n}(t)=\frac{n!}{i!(n-i)!}t^i(1-t)^{n-i}$是Bernstein基
存在很多种定义矩阵权的方法，这里采用一种几何上比较直观的方法定义。假设对于 $\textbf{n}_i,i=0,1,\cdots,n$是一族单位法向量并且 $\omega_i>0,\mu_i>-1,i=0,1,\cdots,n$是一族实数，对于 $Q(\xi)$矩阵权定义如下
$$M_i=\omega_i(I+\mu_i\textbf{n}_i\textbf{n}_i^T),i=0,1,\cdots,n$$
其中I是 $d\times d$的单位阵
这样就可以进行矩阵权有理Bezier曲线的绘制了
这里对文中的图1（a）还原，这里所有的 $\omega_i$都是1，除去 $\mu_2,\mu_5,\mu_8$取20之外所有的 $\mu_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

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