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

时间： 2024-10-01 03:59:41

[工具]Loop subdivision algorithm+matlab code的相关文章

6 Easy Steps to Learn Naive Bayes Algorithm (with code in Python)

6 Easy Steps to Learn Naive Bayes Algorithm (with code in Python) Introduction Here’s a situation you’ve got into: You are working on a classification problem and you have generated your set of hypothesis, created features and discussed the importanc

CEPH CRUSH 算法源码分析 原文CEPH CRUSH algorithm source code analysis

原文地址 CEPH CRUSH algorithm source code analysis http://www.shalandis.com/original/2016/05/19/CEPH-CRUSH-algorithm-source-code-analysis/ 文章比较深入的写了CRUSH算法的原理和过程.通过调试深入的介绍了CRUSH计算的过程.文章中添加了些内容. 写在前面 读本文前,你需要对ceph的基本操作,pool和CRUSH map非常熟悉.并且较深入的读过源码. 分析的方法

支持向量机的smo算法（MATLAB code）

建立smo.m % function [alpha,bias] = smo(X, y, C, tol) function model = smo(X, y, C, tol) % SMO: SMO algorithm for SVM % %Implementation of the Sequential Minimal Optimization (SMO) %training algorithm for Vapnik's Support Vector Machine (SVM) % % This

108节课彻底学通Python-2.2节: 开发工具中的青龙偃月-VS Code

这是一本教同学们彻底学通Python的高质量学习教程,认真地学习每一章节的内容,每天只需学好一节,帮助你成为一名卓越的Python程序员: 本教程面向的是零编程基础的同学,非科班人士,以及有一定编程水平的中高级程序员. 2.2.1 VS Code简介 VS Code的全称是Visual Studio Code,是微软开发的一款垮平台的IDE工具.与“笨重”的PyCharm相比,VS Code非常轻量,启动速度很快. PyCharm之所以“笨重”,是因为它支持更多易于提升开发效率的特性,对初学者比

求平均排序MATLAB code

A0=R(:,1:2:end); for i=1:17 A1=A0(i,:); p=sort(unique(A1)); for j=1:length(p) Rank0(A1==p(j))=j; end Rank(i,:)=Rank0; end RD5=mean(Rank);

word linkage 选择合适的聚类个数matlab code

clear load fisheriris X = meas; m = size(X,2); % load machine % load census % % X = meas; % X=X(1:2000,:); d = pdist(X,'euclidean'); Z = linkage(d,'ward');%Create a hierarchical binary cluster tree using linkage % [Hr,Tr]=dendrogram(Z, 0); %generates

[ML] Gradient Descend Algorithm [Octave code]

function [theta, J_history] = gradientDescentMulti(X, y, theta, alpha, num_iters) m = length(y); % number of training examples J_history = zeros(num_iters, 1); for iter = 1:num_iters theta = theta - alpha * X' * (X * theta - y) / m; iter = iter +1; J

sequential minimal optimization,SMO for SVM, (MATLAB code)

function model = SMOforSVM(X, y, C ) %sequential minimal optimization,SMO tol = 0.001; maxIters = 3000; global i1 i2 K Alpha M1 m1 w b [m, n] = size(X); K = (X*X'); Alpha = zeros(m,1); w = 0; b = 0; flag =1;iters = 1; while flag >0 & iters < max

MFCC matlab code

%function ccc=mfcc(x) %归一化mel滤波器组系数 filename=input('input filename:','s'); [x,fs,bits]=wavread(filename); bank=melbankm(24,256,fs,0,0.5,'m'); bank=full(bank); bank=bank/max(bank(:)); ?T系数,12*24 for k=1:12 n=0:23; dctcoef(k,:)=cos((2*n+1)*k*pi/(2*24))