二面角自适应的loop细分

最新推荐文章于 2020-06-01 16:53:52 发布

原创最新推荐文章于 2020-06-01 16:53:52 发布 · 957 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#二面角 #matlab #loop-细分 #自适应

CAGD/CG 同时被 2 个专栏收录

24 篇文章

订阅专栏

细分曲线曲面

16 篇文章

订阅专栏

本文介绍了自适应Loop细分算法的实现细节及代码示例，并通过人脸模型对比展示了不同二面角阈值下细分效果的差异。同时探讨了该算法在尖锐和平坦区域的网格密度变化，以及迭代过程中可能遇到的问题。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

跟之前的博文二面角自适应的线性细分是同样的思路，其实和之前我介绍过的线性细分和loop细分之间的文章的差别是几乎一样的

具体实现的代码只需要把之前线性细分的代码替换为如下代码就可以了adloopSubdivision.m

function [newVertices, newFaces] =  adloopSubdivision(vertices, faces)  
% Mesh subdivision using the Loop scheme.  
%  
%  Dimensions:  
%    vertices: 3xnVertices  
%    faces:    3xnFaces  
%    
%  Author: Jesus Mena  
   
    global edgeVertice;  
    global newIndexOfVertices;  
    newFaces = [];  
    newVertices = vertices;  
    
  
    nVertices = size(vertices,2);  
    nFaces    = size(faces,2);  
    edgeVertice = zeros(nVertices, nVertices, 3);  
    newIndexOfVertices = nVertices;
    A = compute_edge_face_ring(faces);
    fring = compute_face_ring(faces);
    [valence, num_of_flat_faces]=compute_adface_angle(vertices,faces);

    
    for i=1:nFaces 
        [vaIndex, vbIndex, vcIndex] = deal(faces(1,i), faces(2,i), faces(3,i)); 
        if valence(i)==3
    % ------------------------------------------------------------------------ %  
    % 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  
  
            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];   
            
            
        elseif valence(i)==2
    % ------------------------------------------------------------------------ %  
    % 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.  
    %    flat
    % vd------vb-----ve        vd------vb------ve             
    %  \     / \     /          \     / |\     /
    %   \   /   \   /            \   /  | vq  /
    %    \ /     \ /              \ /   |/ \ / 
    %     va ----- vc   ->         va-- vr --vc   
    %      \     /                  \      /  
    %       \   /                    \    /  
    %        \ /                      \  /  
    %         vf                       vf     
            for k=1:3
                if num_of_flat_faces{i}(1,k)==0
                the_flat_face=fring{i}(1,k);
                %represent the flat face
                break;
                end
            end
            
            if A(faces(1,i), faces(2,i))==the_flat_face||A(faces(2,i),faces(1,i))==the_flat_face
                vqIndex = addEdgeVertice(vbIndex, vcIndex, vaIndex);  
                vrIndex = addEdgeVertice(vaIndex, vcIndex, vbIndex);
                addoriginEdgeVertice(vaIndex, vbIndex);
                threeFaces = [vaIndex,vbIndex,vrIndex; vrIndex,vbIndex,vqIndex; vrIndex,vqIndex,vcIndex]';  
            elseif A(faces(2,i), faces(3,i))==the_flat_face||A(faces(3,i),faces(2,i))==the_flat_face
                vpIndex = addEdgeVertice(vaIndex, vbIndex, vcIndex);  
                vrIndex = addEdgeVertice(vaIndex, vcIndex, vbIndex);
                addoriginEdgeVertice(vbIndex, vcIndex);
                threeFaces = [vaIndex,vpIndex,vrIndex; vpIndex,vbIndex,vcIndex; vrIndex,vpIndex,vcIndex]';  
            else
                vpIndex = addEdgeVertice(vaIndex, vbIndex, vcIndex);  
                vqIndex = addEdgeVertice(vbIndex, vcIndex, vaIndex);
                addoriginEdgeVertice(vaIndex, vcIndex);
                threeFaces = [vaIndex,vpIndex,vqIndex; vpIndex,vbIndex,vqIndex; vaIndex,vqIndex,vcIndex]';  
            end
            newFaces  = [newFaces, threeFaces];
            
        elseif valence(i)==1
            for k=1:3
                if num_of_flat_faces{i}(1,k)==1
                the_not_flat_face=fring{i}(1,k);
                %represent the flat face
                break;
                end
            end
             
            if     A(faces(1,i), faces(2,i))==the_not_flat_face||A(faces(2,i),faces(1,i))==the_not_flat_face       
                vpIndex = addEdgeVertice(vaIndex, vbIndex, vcIndex);
                addoriginEdgeVertice(vaIndex, vcIndex);
                addoriginEdgeVertice(vbIndex, vcIndex);
                twoFaces = [vaIndex,vpIndex,vcIndex; vpIndex,vbIndex,vcIndex]';  
            elseif A(faces(2,i), faces(3,i))==the_not_flat_face||A(faces(3,i),faces(2,i))==the_not_flat_face
                vqIndex = addEdgeVertice(vbIndex, vcIndex, vaIndex);
                addoriginEdgeVertice(vaIndex, vbIndex);
                addoriginEdgeVertice(vaIndex, vcIndex);
                twoFaces = [vaIndex,vbIndex,vqIndex; vaIndex,vqIndex,vcIndex]';  
            else
                vrIndex = addEdgeVertice(vaIndex, vcIndex, vbIndex);
                addoriginEdgeVertice(vaIndex, vbIndex);
                addoriginEdgeVertice(vbIndex, vcIndex);
                twoFaces = [vaIndex,vbIndex,vrIndex; vrIndex,vbIndex,vcIndex]'; 
            end
            newFaces  = [newFaces, twoFaces]; 
        else
            addoriginEdgeVertice(vaIndex, vbIndex);
            addoriginEdgeVertice(vaIndex, vcIndex);
            addoriginEdgeVertice(vbIndex, vcIndex);
            oneFace=[vaIndex, vbIndex, vcIndex]';
            newFaces  = [newFaces, oneFace];
        end
    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  
                    newVertices(:,vNIndex) = 1/2*(vertices(:,v1)+vertices(:,v2));  
                else  
                    newVertices(:,vNIndex) = 3/8*(vertices(:,v1)+vertices(:,v2)) + 1/8*(vertices(:,vNOpposite1Index)+vertices(:,vNOpposite2Index));  
                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
            newVertices(:,v) = 6/8*vertices(:,v) + 1/8*sum(vertices(:,adjBoundaryVertices),2);
        else
            beta = 1/k*( 5/8 - (3/8 + 1/4*cos(2*pi/k))^2 );
            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  

function addoriginEdgeVertice(v1Index, v2Index)  
    global edgeVertice;  
  
    if (v1Index>v2Index) % setting: v1 <= v2  
        vTmp = v1Index;  
        v1Index = v2Index;  
        v2Index = vTmp;  
    end;  
      
    if (edgeVertice(v1Index, v2Index, 1)==0)  % new vertex  
        edgeVertice(v1Index, v2Index, 1) = -1;
    else  
        edgeVertice(v1Index, v2Index, 3) = -1;   
    end;  
end

代码也是修改了好久才通过测试的，看一个人脸的例子哈，只迭代一次


原始网格	二面角阈值为90度	二面角阈值为50度

很明显在眼睛和耳朵等尖锐的地方有更多的网格，而在脸部等较为平坦的地方网格较少。

对于之前我的疑惑我是这么处理的，还是按照传统的loop细分算法的来计算点的坐标，但是其实这样做会有一定的的问题，比如拿正四面体距离，如果二面角的阈值小于其二面角，就不会增加点，但是，点的坐标，却会变动，不过经过测试，对于正四面体来讲四个点会迭代收敛到一个特定的值附件，不过迭代多少次都是正四面体无疑。

我也请教过杨老师，确实自适应的有时会出现不好的点，就是因为细分格式不是均匀的，对于迭代次数，个人感觉也不宜太高，比如上面的人脸，在迭代次数变多之后有些位置的点会失去特征