矩阵权（Matrix weighted）的LOOP细分实现

之前 MATLAB绘制矩阵权（Matrix weighted）有理Bezier曲线提到了矩阵权的方法，现在我将其用到loop细分上，实现矩阵权的loop细分
loop细分的算法在我之前的博客中已经多次提到了，下面将其推广到矩阵权的loop细分上

浙江大学 杨勋年老师论文——Matrix weighted rational curves and surfaces

loop细分规则

1.网格内部V-顶点位置：
设内部顶点$v_0$的相邻点为$v_1、v_2，\cdots ,v_n$，则该顶点更新后位置为$v=(1-n\beta)v_0+\beta\sum_{i=1}^n v_i$，其中$\beta=\frac{1}{n}[\frac{5}{8}-(\frac{3}{8}+\frac{1}{4}\cos\frac{2\pi}{n})^2]$。

2.网格边界V-顶点位置：
　　设边界顶点$v_0$的相邻点为$v_1、v_2$，则该顶点更新后位置为$v = 3/4*v_0 + 1/8*(v_1 + v_2)$。

3.网格内部E-顶点位置：
　　设内部边的两个端点为$v_0、v_1$，相对的两个顶点为$v_2、v_3$，则新增加的顶点位置为$v = 3/8*(v_0 + v_1) + 1/8*(v_2 + v_3)$。

4.网格边界E-顶点位置：
　　设边界边的两个端点为$v_0、v_1$，则新增加的顶点位置为$v = 1/2*(v_0 + v_1)$。

矩阵权的loop细分

权矩阵的计算

1.网格内部V-顶点：
设内部顶点$v_0$的相邻点的权矩阵为$M_1、M_2，\cdots ,M_n$，则该顶点更新后权矩阵为$M_{vi}=(1-n\beta)M_0+\beta\sum_{i=1}^n M_i$，其中$\beta=\frac{1}{n}[\frac{5}{8}-(\frac{3}{8}+\frac{1}{4})\cos\frac{2\pi}{n})^2]$。

2.网格边界V-顶点：
　　设边界顶点$v_0$的相邻点为$M_1、M_2$，则该顶点更新后权矩阵为$M_{ve} = 3/4*M_0 + 1/8*(M_1 + M_2)$。

3.网格内部E-顶点：
　　设内部边的两个端点的权矩阵为$M_0、M_1$，相对的两个顶点为的权矩阵为$M_2、M_3$，则新增加的顶点权矩阵为$M_{ei} = 3/8*(M_0 + M_1) + 1/8*(M_2 + M_3)$。

4.网格边界E-顶点：
　　设边界边的两个端点的权矩阵为$M_0、M_1$，则新增加的顶点权矩阵为$M_{ee} = 1/2*(M_0 + M_1)$。

细分点的计算

1.网格内部V-顶点位置：
该顶点更新后位置为$v=M_{vi}^{-1}(M_0(1-n\beta)v_0+\beta\sum_{i=1}^n M_iv_i)$，其中$\beta=\frac{1}{n}[\frac{5}{8}-(\frac{3}{8}+\frac{1}{4})\cos\frac{2\pi}{n})^2]$。

2.网格边界V-顶点位置：
该顶点更新后位置为$v = M_{ve}^{-1}(3/4*M_0v_0 + 1/8*(M_1v_1 + M_2v_2))$。

3.网格内部E-顶点位置：
新增加的顶点位置为$v = M_{ei}^{-1}(3/8*(M_0v_0 + M_1v_1) + 1/8*(M_2v_2 + M_3v_3))$。

4.网格边界E-顶点位置：
新增加的顶点位置为$v =M_{ee}^{-1}( 1/2*(M_0v_0 + M_1v_1))$。

实现代码

用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细分模型	不按面积加权法向	按面积加权法向

通过上图可以看出，在点较多的地方效果是不错的，几乎可以插值控制点，但是在尖锐的地方效果比较差，这是因为法向的计算经常是不准确的，所以兔子耳朵出现了比较大的变形。
而对于面积进行加权之后效果略微有所改善，说明按照面积加权计算法向相对来说比较准确^^如果能有所改进就更好了