渐进插值的LOOP 曲面细分

本文的算法来源于

哈尔滨理工大学孙立镌，鞠志涛的论文——渐进插值的LOOP 曲面细分

摘要

目前很多细分方法都存在不能用同一种方法处理封闭网格和开放网格的问题。对此，一种新的基于插值技术的LOOP 曲面细分方法，其主要思想就是给定一个初始三角网格M，反复生成新的顶点，新顶点是通过其相邻顶点的约束求解得到的，从而构造一个新的控制网格M¯，在取极限的情况下，可以证明插值过程是收敛的;因为生成新顶点使用的是与其相连顶点的约束求解得到的，本质上是一种局部方法，所以，该方法很容易定义。它在本地方法和全局方法中都有优势，能处理任意顶点数量和任意拓扑结构的网格，从而产生一个光滑的曲面并忠实于给定曲面的形状，其控制网格可以是封闭的或者开放的。我们这里主要讨论的是封闭网格

插值细分的几何细分过程

基于渐进插值的LOOP 曲面细分在封闭网格情况下的概念如下: 给定一个三维三角网格M=M0 ，如果初始网格不是三
角网格要首先用Delaunay 三角化算法对网格进行处理。对LOOP 细分曲面的M0插入新顶点，首先要找到一个封闭控制网格M¯，使其LOOP 曲面通过M0 的所有顶点，其次要找到M¯的顶点和M0的顶点的直接关系。本文使用迭代的过程来进行这项工作。
考虑M0 上的LOOP 曲面S0，对于M0 的每个顶点V0，计算该顶点与它在S0上的极限顶点V0∞的距离为D0=V0−V0∞。将该距离与V0相加得到细分一次后的顶点为V1=V0+D0。
所有顶点V1的集合构成一次细分网格M1，接下来计算S1上的极限顶点V1∞与V0的距离为D1=V0−V1∞。将该距离与V1相加得到细分一次后的顶点为V2=V1+D1。
所有V2构成第二次细分网格M2 ，重复该过程可以得到第k次插值后的控制网格Mk，并且Dk=V0－Vk∞。
将这个距离与Vk相加得到Vk+1为Vk+1=Vk+Dk。所有
Vk+1的集合构成Mk+1。
这个过程产生了一系列控制网格Mk和一系列对应LOOP曲面Sk，如果当k 趋近于无穷时，Mk和Sk的距离趋近于0，就能证明该插值过程是收敛的，即Sk 收敛于M0。接下来就证明这
个结论。
注意到上述迭代过程的主要工作是计算Vk∞，Vk∞是Sk上Vk对应的极限顶点。对于LOOP 曲面，通过价为n 的极限控制顶点的约束求解计算，约束方程如下:

V∞=βnV+(1−βn)Q    (1)βn=311−8×(38+(38+14cosπn)2)(2)Q=1n∑i=1nQi    (3)

联立约束方程式(1) (2) (3) 可以求出Vk∞，其中的Qi是V的邻域顶点。由于该计算过程只涉及V 的邻域顶点，该过程是一个局部方法，任何拓扑结构都可以利用该方法，即可以处理任意拓扑结构的网格。这里需要指出的是不必严格地一直按迭代过程计算下去，当M0 和Sk的距离小于给定的一个值时，就可以认为迭代过程已经符合要求。

程序实现及效果

程序实现

我修改了之前的loop细分程序，新的程序命名为adloopsub.m

function [newVertices, newFaces] =  adloopsub(vertices, faces)  % Mesh subdivision using the Loop scheme.  %  %  Dimensions:  %    vertices: 3xnVertices  %    faces:    3xnFaces  %    %  Author: lafengxiaoyu     global edgeVertice;      global newIndexOfVertices;      newVertices = vertices;    [originver, orifaces] = origin();    nVertices = size(vertices,2);      nFaces    = size(faces,2);      edgeVertice = zeros(nVertices, nVertices, 3);      newDeff=vertices;    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));          addEdgeVertice(vaIndex, vbIndex, vcIndex);          addEdgeVertice(vbIndex, vcIndex, vaIndex);          addEdgeVertice(vaIndex, vcIndex, vbIndex);      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;            beta = 3/(11-8*(3/8+((3/8 + 1/4*cos(2*pi/k))^2)));          newVertices(:,v) = (beta)*vertices(:,v) + (1-beta)/k*sum(vertices(:,(adjVertice{v})),2);             newDeff(:,v)=originver(:,v)- newVertices(:,v);          newVertices(:,v)=vertices(:,v)+newDeff(:,v);    end;        newFaces=faces;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  

其中origin.m用来读取数据

function [ vertices, faces ] = origin()%ORIGIN 此处显示有关此函数的摘要%   此处显示详细说明%faces =load('C:\Users\Admin\Documents\MATLAB\face.txt');%faces=faces+1;%faces=faces';%vertices =load('C:\Users\Admin\Documents\MATLAB\ver.txt');%vertices=vertices';[vertices, faces ]= obj__read( 'bunny_200.obj' );%vertices=[-1 -1 0;-1 1 0;1 1 0]';%faces=[1 2 3]';end

效果

执行以下测试代码

[ vertices, faces ]=origin();faces1=faces;vertices1=vertices;%for i=1:5%[vertices, faces] = adloopsub(vertices, faces);%endfor i=1:3[vertices, faces] = loopSubdivision(vertices, faces);endtrimesh(faces', vertices(1,:), vertices(2,:), vertices(3,:),'LineWidth',1,'EdgeColor','k');alpha(0.0); hold on;trimesh(faces1', vertices1(1,:), vertices1(2,:), vertices1(3,:),'EdgeColor','r');alpha(0.0);obj_write('bunny1.obj',vertices,faces);

读入由200个点组成的兔子模型，分别进行loop细分和渐进loop细分，效果如下

loop细分三次后的兔子模型

渐进插值后loop三次的兔子模型

可以明显看得出来效果好一些

参考文献

［1］ LOOP C． Smooth subdivision surfaces based on triangles［D］． Utah:University of Utah，1987．
［2］ HALSTEAD M，KASS M，DEROSE T． Efficient，fair interpolation
using Catmull-Clark surfaces［C］/ /Proc of the 20th SIGGRAPH．1993: 47-61．
［3］ ZHENG Jian-min，CAI Y Yi-yu． Interpolation over arbitrary topology meshes using a two-phase subdivision scheme［J］． IEEE Trans on Visualization and Computer Graphics，2006，12( 3) : 301-310．
［4］ LAI Shu-hua，CHENG F． Similarity based interpolation using Catmull-Clark subdivision surfaces［J］． The Visual Computer，2006，22( 9) : 865-873．
［5］ ZORIN D，SCHRODER P，DEROSE T． Subdivision for modelingand animation［C］/ /Proc of SIGGRAPH． 2006: 30-50．
［6］ CHENG Fu-hua，FAN Feng-tao，LAI Shu-hua，et al． LOOP subdivision surface based progressive interpolation［J］． Journal of Computer Science and Technology，2009，24( 1) : 39-46．