[工具]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;

% ------------------------------------------------------------------------ %

% 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;

% ------------------------------------------------------------------------ %

% 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;

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

% ---------------------------------------------------------------------------- %

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