[工具]Loop subdivision algorithm+matlab code
引言：分析Loop subdivision algorithm的原理及代码实现。该算法可以用在对3D网格的细分上。同时，在形状检索领域，经常需要选取视点（viewpoint）来对模型进行绘制（render），例如，在zhouhui lian的CM-BOF方法中，对正八面体进行细化分解得到均匀的视点分布。
关于原理，有一些中文博客，但是讲的不是太详细，对照代码可能可以看的更清楚。
直接上代码吧。中文注释是我注释的，代码是在mathwork上看到的。
function [newVertices, newFaces] = loopSubdivision(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;
% ------------------------------------------------------------------------ %
% 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)%v1，v2相连，vNIndex表示v1，v2之间的顶点序列号
vNOpposite1Index = edgeVertice(v1,v2,2);
vNOpposite2Index = edgeVertice(v1,v2,3);
if (vNOpposite2Index==0) % boundary case当v1v2确定的边为边界时，新坐标即为v1，v2中点
newVertices(:,vNIndex) = 1/2*(vertices(:,v1)+vertices(:,v2));
else %否则新顶点坐标为3/8*(v1+v2)+1/8*(v3+v4)，v3，v4表示的是v1,v2两边的顶点。如上图求vr，则v1，v2分别指代va，vc；v3，v4分别指代vb，vd
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)
%如果v和vTmp是邻接的，那么v的邻接表就指向vTmp
adjVertice{v}(end+1) = vTmp;
end;
end;
end;
% ------------------------------------------------------------------------ %
% new positions of the original vertices
for v=1:nVertices
k = length(adjVertice{v}); %k为邻接顶点个数
adjBoundaryVertices = [];
for i=1:k
vi = adjVertice{v}(i);
if (vi>v) && (edgeVertice(v,vi,3)==0) || (vi<v) && (edgeVertice(vi,v,3)==0)
%边界的情况，将vi放在边界表里
adjBoundaryVertices(end+1) = vi;
end;
end;
if (length(adjBoundaryVertices)==2) % boundary case
%如果边界顶点v在边界上的两相邻顶点为v0,v1，则新顶点的坐标为6/8*v+1/8*v0+1/8*v1
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
%如果v1和v2之间没有新的顶点（也就是说v1，v2确定的边为图形的边界），则生成新的顶点，将新顶点的序号赋给v1，v2之间的新顶点
%同时将v3的序号赋给v1，v2之间的第一相反顶点（first opposite vertex）
newIndexOfVertices = newIndexOfVertices+1;
edgeVertice(v1Index, v2Index, 1) = newIndexOfVertices;
edgeVertice(v1Index, v2Index, 2) = v3Index;
else
%否则，将v3赋给第二相反顶点（second opposite vertex）
edgeVertice(v1Index, v2Index, 3) = v3Index;
end;
%返回v1，v2之间的顶点
vNIndex = edgeVertice(v1Index, v2Index, 1);
return;
end
难理解的部分是，这里使用了edgeVertice，adjVertice，adjBoundaryVertices三个变量。
其中edgeVertice是一个全局变量，大小为nVertex*nVertex*3。
edgeVertice（x，y，1）记录了顶点x，y之间的顶点序列，默认没有顶点，值为0.
edgeVertice（x，y，2）记录了顶点x，y确定的边的两边的顶点序列。
edgeVertice（x，y，3）同2，但是在该边为边界时，该值为0.
其中adjVertice保存顶点的邻接顶点。
其中adjBoundaryVertices保存为边界的顶点。
下面是MATLAB的测试与显示代码
function plotMesh(vertices, faces)
hold on;
trimesh(faces', vertices(1,:), vertices(2,:), vertices(3,:));
colormap gray(1);
axis tight;
axis square;
axis off;
view(3);
end
% Test: Mesh subdivision using the Loop scheme.
%
% Author: Jesus Mena
% Example: Box
vertices = [10 10 10; -10 10 10; 10 -10 10; -10 -10 10; 10 10 -10; -10 10 -10; 10 -10 -10; -10 -10 -10]';
faces = [1 2 3; 4 3 2; 1 3 5; 7 5 3; 1 5 2; 6 2 5; 8 6 7; 5 7 6; 8 7 4; 3 4 7; 8 4 6; 2 6 4]';
figure(1);
subplot(1,4,1);
plotMesh(vertices, faces);
for i=2:4
subplot(1,4,i);
[vertices, faces] = loopSubdivision(vertices, faces);
plotMesh(vertices, faces);
end
% Example: Tetrahedron
vertices = [10 10 10; -100 10 -10; -100 -10 10; 10 -10 -10]';
faces = [1 2 3; 1 3 4; 1 4 2; 4 3 2]';
figure(2);
subplot(1,4,1);
plotMesh(vertices, faces);
for i=2:4
subplot(1,4,i);
[vertices, faces] = loopSubdivision(vertices, faces);
plotMesh(vertices, faces);
end
% Example: Cylinder
vertices = [0 -5 0; 0 5 0; 10 -5 0; 9.65 -5 2.58; 8.66 -5 5; 7.07 -5 7.07; 5 -5 8.66; 2.58 -5 9.65; 0 -5 10; -2.58 -5 9.65; -5 -5 8.66; -7.07 -5 7.07; -8.66 -5 5; -9.65 -5 2.58; -10 -5 0; -9.65 -5 -2.58; -8.66 -5 -5; -7.07 -5 -7.07; -5 -5 -8.66; -2.58 -5 -9.65; -0 -5 -10; 2.58 -5 -9.65; 5 -5 -8.66; 7.07 -5 -7.07; 8.66 -5 -5; 9.65 -5 -2.58; 10 5 0 ; 9.65 5 2.58; 8.66 5 5; 7.07 5 7.07; 5 5 8.66; 2.58 5 9.65; 0 5 10; -2.58 5 9.65; -5 5 8.66; -7.07 5 7.07; -8.66 5 5; -9.65 5 2.58; -10 5 0; -9.65 5 -2.58; -8.66 5 -5; -7.07 5 -7.07; -5 5 -8.66; -2.58 5 -9.65; -0 5 -10; 2.58 5 -9.65; 5 5 -8.66; 7.07 5 -7.07; 8.66 5 -5; 9.65 5 -2.58]';
faces = [1 3 4; 1 4 5; 1 5 6; 1 6 7; 1 7 8; 1 8 9; 1 9 10; 1 10 11; 1 11 12; 1 12 13; 1 13 14; 1 14 15; 1 15 16; 1 16 17; 1 17 18; 1 18 19; 1 19 20; 1 20 21; 1 21 22; 1 22 23; 1 23 24; 1 24 25; 1 25 26; 1 26 3; 2 28 27; 2 29 28; 2 30 29; 2 31 30; 2 32 31; 2 33 32; 2 34 33; 2 35 34; 2 36 35; 2 37 36; 2 38 37; 2 39 38; 2 40 39; 2 41 40; 2 42 41; 2 43 42; 2 44 43; 2 45 44; 2 46 45; 2 47 46; 2 48 47; 2 49 48; 2 50 49; 2 27 50; 3 27 28; 3 28 4; 4 28 29; 4 29 5; 5 29 30; 5 30 6; 6 30 31; 6 31 7; 7 31 32; 7 32 8; 8 32 33; 8 33 9; 9 33 34; 9 34 10; 10 34 35; 10 35 11; 11 35 36; 11 36 12; 12 36 37; 12 37 13; 13 37 38; 13 38 14; 14 38 39; 14 39 15; 15 39 40; 15 40 16; 16 40 41; 16 41 17; 17 41 42; 17 42 18; 18 42 43; 18 43 19; 19 43 44; 19 44 20; 20 44 45; 20 45 21; 21 45 46; 21 46 22; 22 46 47; 22 47 23; 23 47 48; 23 48 24; 24 48 49; 24 49 25; 25 49 50; 25 50 26; 26 50 27; 26 27 3]';
figure(3);
subplot(1,4,1);
plotMesh(vertices, faces);
for i=2:4
subplot(1,4,i);
[vertices, faces] = loopSubdivision(vertices, faces);
plotMesh(vertices, faces);
end
% Example: Grid
vertices = [-4 -4 0; -2 -4 0; 0 -4 0; 2 -4 0; 4 -4 0; -4 -2 0; -2 -2 0; 0 -2 0; 2 -2 0; 4 -2 0; -4 0 0; -2 0 0; 0 0 0; 2 0 0; 4 0 0; -4 2 0; -2 2 0; 0 2 0; 2 2 0; 4 2 0; -4 4 0; -2 4 0; 0 4 0; 2 4 0; 4 4 0]';
faces = [7 2 1; 1 6 7; 8 3 2; 2 7 8; 9 4 3; 3 8 9; 10 5 4; 4 9 10; 12 7 6; 6 11 12; 13 8 7; 7 12 13; 14 9 8; 8 13 14; 15 10 9; 9 14 15; 17 12 11; 11 16 17; 18 13 12; 12 17 18; 19 14 13; 13 18 19; 20 15 14; 14 19 20; 22 17 16; 16 21 22; 23 18 17; 17 22 23; 24 19 18; 18 23 24; 25 20 19; 19 24 25]';
figure(4);
subplot(1,4,1);
plotMesh(vertices, faces);
for i=2:4
subplot(1,4,i);
[vertices, faces] = loopSubdivision(vertices, faces);
plotMesh(vertices, faces);
end
