之前 MATLAB绘制矩阵权（Matrix weighted）有理Bezier曲线提到了矩阵权的方法，现在我将其用到loop细分上，实现矩阵权的loop细分
loop细分的算法在我之前的博客中已经多次提到了，下面将其推广到矩阵权的loop细分上
浙江大学 杨勋年老师论文——Matrix weighted rational curves and surfaces
loop细分规则
1.网格内部V-顶点位置：
设内部顶点v0的相邻点为v1、v2，⋯,vn，则该顶点更新后位置为v=(1−nβ)v0+β∑ni=1vi，其中β=1n[58−(38+14cos2πn)2]。
2.网格边界V-顶点位置：
设边界顶点v0的相邻点为v1、v2，则该顶点更新后位置为v=3/4∗v0+1/8∗(v1+v2)。
3.网格内部E-顶点位置：
设内部边的两个端点为v0、v1，相对的两个顶点为v2、v3，则新增加的顶点位置为v=3/8∗(v0+v1)+1/8∗(v2+v3)。
4.网格边界E-顶点位置：
设边界边的两个端点为v0、v1，则新增加的顶点位置为v=1/2∗(v0+v1)。
矩阵权的loop细分
权矩阵的计算
1.网格内部V-顶点：
设内部顶点v0的相邻点的权矩阵为M1、M2，⋯,Mn，则该顶点更新后权矩阵为Mvi=(1−nβ)M0+β∑ni=1Mi，其中β=1n[58−(38+14)cos2πn)2]。
2.网格边界V-顶点：
设边界顶点v0的相邻点为M1、M2，则该顶点更新后权矩阵为Mve=3/4∗M0+1/8∗(M1+M2)。
3.网格内部E-顶点：
设内部边的两个端点的权矩阵为M0、M1，相对的两个顶点为的权矩阵为M2、M3，则新增加的顶点权矩阵为Mei=3/8∗(M0+M1)+1/8∗(M2+M3)。
4.网格边界E-顶点：
设边界边的两个端点的权矩阵为M0、M1，则新增加的顶点权矩阵为Mee=1/2∗(M0+M1)。
细分点的计算
1.网格内部V-顶点位置：
该顶点更新后位置为v=M−1vi(M0(1−nβ)v0+β∑ni=1Mivi)，其中β=1n[58−(38+14)cos2πn)2]。
2.网格边界V-顶点位置：
该顶点更新后位置为v=M−1ve(3/4∗M0v0+1/8∗(M1v1+M2v2))。
3.网格内部E-顶点位置：
新增加的顶点位置为v=M−1ei(3/8∗(M0v0+M1v1)+1/8∗(M2v2+M3v3))。
4.网格边界E-顶点位置：
新增加的顶点位置为v=M−1ee(1/2∗(M0v0+M1v1))。
实现代码
用MATLAB修改之前的loop细分代码，如下
function [newVertices, newFaces,newM] = mwloopSubdivision(vertices, faces,M)
% Mesh subdivision using the Loop scheme.
%
% Dimensions:
% vertices: 3xnVertices
% faces: 3xnFaces
%
% Author: Jesus Mena
global edgeVertice;
global newIndexOfVertices;
newFaces = [];
newVertices = vertices;
newM=M;
nVertices = size(vertices,2);
nFaces = size(faces,2);
edgeVertice = zeros(nVertices, nVertices, 3);
newIndexOfVertices = nVertices;
% ------------------------------------------------------------------------ %
% create a matrix of edge-vertices and the new triangulation (newFaces).
% computational complexity = O(3*nFaces)
%
% * edgeVertice(x,y,1): index of the new vertice between (x,y)
% * edgeVertice(x,y,2): index of the first opposite vertex between (x,y)
% * edgeVertice(x,y,3): index of the second opposite vertex between (x,y)
%
% 0riginal vertices: va, vb, vc, vd.
% New vertices: vp, vq, vr.
%
% vb vb
% / \ / \
% / \ vp--vq
% / \ / \ / \
% va ----- vc -> va-- vr --vc
% \ / \ /
% \ / \ /
% \ / \ /
% vd vd
for i=1:nFaces
[vaIndex, vbIndex, vcIndex] = deal(faces(1,i), faces(2,i), faces(3,i));
vpIndex = addEdgeVertice(vaIndex, vbIndex, vcIndex);
vqIndex = addEdgeVertice(vbIndex, vcIndex, vaIndex);
vrIndex = addEdgeVertice(vaIndex, vcIndex, vbIndex);
fourFaces = [vaIndex,vpIndex,vrIndex; vpIndex,vbIndex,vqIndex; vrIndex,vqIndex,vcIndex; vrIndex,vpIndex,vqIndex]';
newFaces = [newFaces, fourFaces];
end;
% ------------------------------------------------------------------------ %
% positions of the new vertices
for v1=1:nVertices-1
for v2=v1:nVertices
vNIndex = edgeVertice(v1,v2,1);
if (vNIndex~=0)
vNOpposite1Index = edgeVertice(v1,v2,2);
vNOpposite2Index = edgeVertice(v1,v2,3);
if (vNOpposite2Index==0) % boundary case
%矩阵权内容
ME=1/2*(M{v1}+M{v2});
newVertices(:,vNIndex) = (ME)\(1/2*(M{v1}*vertices(:,v1)+M{v2}*vertices(:,v2)));
newM{vNIndex}=ME;
else
%矩阵权内容
ME=3/8*(M{v1}+M{v2})+1/8*(M{vNOpposite1Index}+M{vNOpposite2Index});
newVertices(:,vNIndex) = (ME)\(3/8*(M{v1}*vertices(:,v1)+M{v2}*vertices(:,v2)) + 1/8*(M{vNOpposite1Index}*vertices(:,vNOpposite1Index)+M{vNOpposite2Index}*vertices(:,vNOpposite2Index)));
newM{vNIndex}=ME;
end;
end;
end;
end;
% ------------------------------------------------------------------------ %
% adjacent vertices (using edgeVertice)
adjVertice{nVertices} = [];
for v=1:nVertices
for vTmp=1:nVertices
if (v<vTmp && edgeVertice(v,vTmp,1)~=0) || (v>vTmp && edgeVertice(vTmp,v,1)~=0)
adjVertice{v}(end+1) = vTmp;
end;
end;
end;
% ------------------------------------------------------------------------ %
% new positions of the original vertices
for v=1:nVertices
k = length(adjVertice{v});
adjBoundaryVertices = [];
for i=1:k
vi = adjVertice{v}(i);
if (vi>v) && (edgeVertice(v,vi,3)==0) || (vi<v) && (edgeVertice(vi,v,3)==0)
adjBoundaryVertices(end+1) = vi;
end;
end;
if (length(adjBoundaryVertices)==2) % boundary case
MP=6/8*M{v}+1/8*(M{adjBoundaryVertices(1)}+M{adjBoundaryVertices(2)});
newVertices(:,v) = MP\(6/8*M{v}*vertices(:,v) + 1/8*(M{adjBoundaryVertices(1)}*vertices(:,adjBoundaryVertices(1))+M{adjBoundaryVertices(2)}*vertices(:,adjBoundaryVertices(2))));
newM{v}=MP;
else
beta = 1/k*( 5/8 - (3/8 + 1/4*cos(2*pi/k))^2 );
MP=(1-k*beta)*M{v};
Msum=(1-k*beta)*M{v}*vertices(:,v);
for i=1:k
MP=MP+ beta*M{adjVertice{v}(i)};
Msum=Msum+beta*M{adjVertice{v}(i)}*vertices(:,(adjVertice{v}(i)));
end
newVertices(:,v) =MP\(Msum);
newM{v}=MP;
% newVertices(:,v) = (1-k*beta)*vertices(:,v) + beta*sum(vertices(:,(adjVertice{v})),2);
end;
end;
end
% ---------------------------------------------------------------------------- %
function vNIndex = addEdgeVertice(v1Index, v2Index, v3Index)
global edgeVertice;
global newIndexOfVertices;
if (v1Index>v2Index) % setting: v1 <= v2
vTmp = v1Index;
v1Index = v2Index;
v2Index = vTmp;
end;
if (edgeVertice(v1Index, v2Index, 1)==0) % new vertex
newIndexOfVertices = newIndexOfVertices+1;
edgeVertice(v1Index, v2Index, 1) = newIndexOfVertices;
edgeVertice(v1Index, v2Index, 2) = v3Index;
else
edgeVertice(v1Index, v2Index, 3) = v3Index;
end;
vNIndex = edgeVertice(v1Index, v2Index, 1);
return;
end调用的方法和之前的loop细分是一样的
演示及对比
显示的代码为
function display_origin_modol_with_normal
%DISPLAY_ORIGIN_MODOL_WITH_NORMAL 此处显示有关此函数的摘要
% 此处显示详细说明
name = 'bunny_200.obj';
OBJ=readObj(name);vertices=OBJ.v;faces=OBJ.f.v;
% trimesh(fac, OBJ.v(:,1), OBJ.v(:,2), OBJ.v(:,3),'LineWidth',1,'EdgeColor','k');alpha(0.0); hold on;
% for i=1:1289
% quiver3(ver(i,1),ver(i,2),ver(i,3),nor(i,1)*0.1,nor(i,2)*0.1,nor(i,3)*0.1,'r','filled','LineWidth',0.1);hold on;
% end %
%options.name = name;
[normal] = compute_normal(vertices,faces);
%options.normal = normal;
%clf; plot_mesh(vertices,faces,options);
%shading interp; %axis tight;
[normal_weight]=compute_normal_weight(vertices,faces);
vertices=vertices';faces=faces';NumPoint=size(vertices,2);
%点的个数
M=cell(1,NumPoint);
M_weight=cell(1,NumPoint);
%对其单位化
N=zeros(3,NumPoint);
N_weight=zeros(3,NumPoint);
for i=1:NumPoint
S=sqrt(normal(1,i)^2+normal(2,i)^2+normal(3,i)^2);
N(1,i)=normal(1,i)/S;N(2,i)=normal(2,i)/S;N(3,i)=normal(3,i)/S;
S_weight=sqrt(normal_weight(1,i)^2+normal_weight(2,i)^2+normal_weight(3,i)^2);
N_weight(1,i)=normal_weight(1,i)/S_weight;N_weight(2,i)=normal_weight(2,i)/S_weight;N_weight(3,i)=normal_weight(3,i)/S_weight;
end
omega=ones(1,NumPoint);
%mu=ones(1,NumPoint);
mu=compute_mu(vertices,faces,N);
mu_weight=compute_mu(vertices,faces,N_weight);
%计算矩阵权
I=eye(3);
for i=1:NumPoint
M{i}=omega(i)*(I+mu(i)*N(:,i)*N(:,i)');
M_weight{i}=omega(i)*(I+mu_weight(i)*N_weight(:,i)*N_weight(:,i)');
end
faces1=faces;vertices1=vertices;verticesl=vertices;facesl=faces;
for i=1:2
[vertices, faces,M] = mwloopSubdivision(vertices, faces,M);
[verticesl, facesl,M_weight] = mwloopSubdivision(verticesl,facesl,M_weight);
end
figure(2);
trimesh(faces', vertices(1,:), vertices(2,:), vertices(3,:),'LineWidth',1,'EdgeColor','k');alpha(0.0); hold on;axis off;trimesh(faces1', vertices1(1,:), vertices1(2,:), vertices1(3,:),'EdgeColor','r');alpha(0.0);
obj_write('bunny_200_matrix.obj',vertices, faces);
figure(3)
trimesh(facesl', verticesl(1,:), verticesl(2,:), verticesl(3,:),'LineWidth',1,'EdgeColor','k');alpha(0.0); hold on;axis off;trimesh(faces1', vertices1(1,:), vertices1(2,:), vertices1(3,:),'EdgeColor','r');alpha(0.0);
obj_write('bunny_200_loop.obj',verticesl, facesl);
end另附readObj.m文件
function obj = readObj(fname)
%
% obj = readObj(fname)
%
% This function parses wavefront object data
% It reads the mesh vertices, texture coordinates, normal coordinates
% and face definitions(grouped by number of vertices) in a .obj file
%
%
% INPUT: fname - wavefront object file full path
%
% OUTPUT: obj.v - mesh vertices
% : obj.vt - texture coordinates
% : obj.vn - normal coordinates
% : obj.f - face definition assuming faces are made of of 3 vertices
%
% Bernard Abayowa, Tec^Edge
% 11/8/07
% set up field types
v = []; vt = []; vn = []; f.v = []; f.vt = []; f.vn = [];
fid = fopen(fname);
% parse .obj file
while 1
tline = fgetl(fid);
if ~ischar(tline), break, end % exit at end of file
ln = sscanf(tline,'%s',1); % line type
%disp(ln)
switch ln
case 'v' % mesh vertexs
v = [v; sscanf(tline(2:end),'%f')'];
case 'vt' % texture coordinate
vt = [vt; sscanf(tline(3:end),'%f')'];
case 'vn' % normal coordinate
vn = [vn; sscanf(tline(3:end),'%f')'];
case 'f' % face definition
fv = []; fvt = []; fvn = [];
str = textscan(tline(2:end),'%s'); str = str{1};
nf = length(findstr(str{1},'/')); % number of fields with this face vertices
[tok str] = strtok(str,'//'); % vertex only
for k = 1:length(tok) fv = [fv str2num(tok{k})]; end
if (nf > 0)
[tok str] = strtok(str,'//'); % add texture coordinates
for k = 1:length(tok) fvt = [fvt str2num(tok{k})]; end
end
if (nf > 1)
[tok str] = strtok(str,'//'); % add normal coordinates
for k = 1:length(tok) fvn = [fvn str2num(tok{k})]; end
end
f.v = [f.v; fv]; f.vt = [f.vt; fvt]; f.vn = [f.vn; fvn];
end
end
fclose(fid);
% set up matlab object
obj.v = v; obj.vt = vt; obj.vn = vn; obj.f = f;其中两种表示分别用上一篇博客的两种方法计算法向得到的不同模型。
原始模型（200点） | loop细分模型 | 不按面积加权法向 | 按面积加权法向 |
|
|
|
|
通过上图可以看出，在点较多的地方效果是不错的，几乎可以插值控制点，但是在尖锐的地方效果比较差，这是因为法向的计算经常是不准确的，所以兔子耳朵出现了比较大的变形。
而对于面积进行加权之后效果略微有所改善，说明按照面积加权计算法向相对来说比较准确^^如果能有所改进就更好了