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;endniu(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;endI=[1 0;0 1];for i=1:NumPoint+1M{i}=omega(i)*(I+niu(i)*N(:,i)*N(:,i)');endxx=[];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;    endnewPoint = invMat\Mat;xx=[xx newPoint(1,1)];yy=[yy newPoint(2,1)];endplot(vertices(1,:),vertices(2,:),'b');hold on;for i=1:NumPoint+1quiver(vertices(1,i),vertices(2,i),N(1,i),N(2,i),'r','filled','LineWidth',2);endhold on;grid on;axis tight;axis equal;xlabel('X');ylabel('Y');plot(xx,yy,'r');end

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