矩阵权（Matrix weighted）Bezier三角（曲面）片

参考文献仍然是杨老师的这篇

Matrix weighted rational curves and surfaces

文章，结合上篇博文中的算法，再次将其引申到三角面片上

算法描述

假设Pi,j,k,0≤i,j,k≤n,i+j+k=n是给定的控制点，并且Mijk是对应的加权矩阵。矩阵权有理三角面片由下式定义

Q3(u,v,w)=[∑i+j+k=nMijkBni,j,k(u,v,w)]−1∑i+j+k=nMijkPi,j,kBni,j,k(u,v,w)，(u,v)∈[uk1−1,um+1]×[vk2−1,vn+1]

其中Bni,j,k(u,v,w)=n!i!j!k!是定义在三角域上的Berstein基函数。

代码实现

用matlab编写的如下程序为实现函数wm_tri.m

function trepBer = wm_tri( controls,num,M,N )%WM_TRI Summary of this function goes here%   Detailed explanation goes herepoints=cell(num,num);for i=1:num    for j=1:i        points{i,j}=[(j-1)*(1/(num-1)),1-(i-1)*(1/(num-1)),(i-j)*(1/(num-1))];    endendn=size(controls,1);Bpoints=cell(num,num);for i=1:num    for j=1:i        Bpoints{i,j}=[0;0;0];        tempM=zeros(3);        for u=1:n            for v=1:u                %i=v-1;j=n-u;k=u-v                tempM=tempM+M{u+v-1}*factorial(n-1)/(factorial(v-1)*factorial(n-u)*factorial(u-v))*(points{i,j}(1)^(v-1)*points{i,j}(2)^(n-u)*points{i,j}(3)^(u-v));                Bpoints{i,j}=Bpoints{i,j}+M{u+v-1}*controls{u,v}'*factorial(n-1)/(factorial(v-1)*factorial(n-u)*factorial(u-v))*(points{i,j}(1)^(v-1)*points{i,j}(2)^(n-u)*points{i,j}(3)^(u-v));            end        end        Bpoints{i,j}=(tempM\Bpoints{i,j})';    endendtriConPoi=zeros(n*(n+1)/2,3);temp=1;for i=1:n    for j=1:i       triConPoi(temp,:)=controls{i,j};       temp=temp+1;    endendtriConSur=zeros((n-1)^2,3);temp=1;for i=2:n%the first in the i-th line is n(n-1)/2    for j=1:i-1        %Counter clockwise order,up triangle        triConSur(temp,:)=[i*(i-1)/2+j,i*(i-1)/2+j+1,i*(i-1)/2+j+1-i];        temp=temp+1;    endendfor i=2:n-1%the first in the i-th line is n(n-1)/2    for j=1:i-1        %Counter clockwise order,down triangle        triConSur(temp,:)=[i*(i-1)/2+j+1,i*(i-1)/2+j,i*(i-1)/2+j+1+i];        temp=temp+1;    endendtrepCon=triangulation(triConSur,triConPoi); triBezPoi=zeros(n*(n+1)/2,3);temp=1;for i=1:num    for j=1:i       triBezPoi(temp,:)=Bpoints{i,j};       temp=temp+1;    endendtriBerSur=zeros((num-1)^2,3);temp=1;for i=2:num%the first in the i-th line is n(n-1)/2    for j=1:i-1        %Counter clockwise order,up triangle        triBerSur(temp,:)=[i*(i-1)/2+j,i*(i-1)/2+j+1,i*(i-1)/2+j+1-i];        temp=temp+1;    endendfor i=2:num-1%the first in the i-th line is n(n-1)/2    for j=1:i-1        %Counter clockwise order,down triangle        triBerSur(temp,:)=[i*(i-1)/2+j+1,i*(i-1)/2+j,i*(i-1)/2+j+1+i];        temp=temp+1;    endendtrepBer=triangulation(triBerSur,triBezPoi); trisurf(trepCon,'edgecolor','k','FaceColor', 'interp');alpha(.2);axis equal;hold on;trisurf(trepBer,'edgecolor','none','FaceColor', 'interp');alpha(.9);hold on;for i=1:3quiver3(triConPoi(i,1),triConPoi(i,2),triConPoi(i,3),N(i,1),N(i,2),N(i,3),'r','filled','LineWidth',2);hold on;endend

示例

下面绘制一个瞅瞅哈

controls = cell(2,2);controls{1,1}=[0,0,0];controls{2,1}=[0.5,1,-0.2];controls{2,2}=[-0.4,0.9,0.1];N=[0.1,-0.3,0.6;0.4,0.3,0.6;-0.3,0.4,0.4];omega=ones(3,1);niu=ones(3,1);M=cell(3,1);for i=1:3S=sqrt(N(i,1)^2+N(i,2)^2+N(i,3)^2);N(i,1)=N(i,1)/S;N(i,2)=N(i,2)/S;N(i,3)=N(i,3)/S;endI=[1 0 0;0 1 0;0 0 1];for i=1:3M{i}=omega(i)*(I+niu(i)*N(i,:)'*N(i,:));endtrepBer = wm_tri( controls,30,M,N);

结果如下图

再来一个例子

来一整个模型瞅瞅，比如说人脸的例子

原本的网格 Bezier曲面片

老铁没毛病