引用来源
参考：Adaptive Subdivision Schemes for Triangular Meshes这篇文章
先按照上面的方法做了线性细分，看起来并没有什么问题，下面详细介绍一下是怎么做的
这篇文章是相当经典的一篇文章，基本上是三角网格自适应细分的先驱，我看到由非常高的引用，该方法主要是基于计算二面角得到的，主要计算一个面和周围的面的夹角来确定这个面是不是应该细分。
算法步骤
算法主要分为四步：
1.对于每个面片计算其法向
2.对于每个面计算其与邻接边的夹角（二面角）
3.如果所有的二面角都小于一个设定的阈值，该面被标记为平坦的flat
4.对于每一个面，与周围三个三角形的二面角是否被标记为flat存在好几种情况，其最小值为0，最大值为3。对于不同的情况，细分情况如下
其中n=1的时候其实有两种不同的细分方案，另外一种我拿橙色的线标记出了，其实两种方法都是可以的，下面就要介绍下怎样进行实现了
MATLAB实现
首先是计算法向的函数compute_face_normal.m
function [normalf] = compute_face_normal(vertex,face)
% compute_normal - compute the normal of a triangulation
%
% [normal,normalf] = compute_normal(vertex,face);
%
% normal(i,:) is the normal at vertex i.
% normalf(j,:) is the normal at face j.
%
% Copyright (c) 2004 Gabriel Peyr?
[vertex,face] = check_face_vertex(vertex,face);
nvert = size(vertex,2);
normal = zeros(3, nvert);
% unit normals to the faces
normalf = crossp( vertex(:,face(2,:))-vertex(:,face(1,:)), ...
vertex(:,face(3,:))-vertex(:,face(1,:)) );
d = sqrt( sum(normalf.^2,1) ); d(d<eps)=1;
normalf = normalf ./ repmat( d, 3,1 );
% unit normal to the vertex
normal = zeros(3,nvert);
% enforce that the normal are outward
v = vertex - repmat(mean(vertex,1), 3,1);
s = sum( v.*normal, 2 );
if sum(s>0)<sum(s<0)
% flip
normalf = -normalf;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z = crossp(x,y)
% x and y are (m,3) dimensional
z = x;
z(1,:) = x(2,:).*y(3,:) - x(3,:).*y(2,:);
z(2,:) = x(3,:).*y(1,:) - x(1,:).*y(3,:);
z(3,:) = x(1,:).*y(2,:) - x(2,:).*y(1,:);
end计算边对应面的函数compute_edge_face_ring.m
输入的是面，输出的是一个稀疏矩阵，是n*n的矩阵，n代表点的个数，第(i,j)位置表示i，j位置的点连线对应的面的索引
function A = compute_edge_face_ring(face)
% compute_edge_face_ring - compute faces adjacent to each edge
%
% e2f = compute_edge_face_ring(face);
%
% e2f(i,j) and e2f(j,i) are the number of the two faces adjacent to
% edge (i,j).
%
% Copyright (c) 2007 Gabriel Peyre
[~,face] = check_face_vertex([],face);
n = max(face(:));
m = size(face,2);
i = [face(1,:) face(2,:) face(3,:)];
j = [face(2,:) face(3,:) face(1,:)];
s = [1:m 1:m 1:m];
% first without duplicate
[~,I] = unique( i+(max(i)+1)*j );
% remaining items
J = setdiff(1:length(s), I);
% flip the duplicates
i1 = [i(I) j(J)];
j1 = [j(I) i(J)];
s = [s(I) s(J)];
% remove doublons
[~,I] = unique( i1+(max(i1)+1)*j1 );
i1 = i1(I); j1 = j1(I); s = s(I);
A = sparse(i1,j1,s,n,n);
% add missing points
I = find( A'~=0 );
I = I( A(I)==0 );
A( I ) = -1;
end
计算一个面邻接面的compute_face_ring.m
同样输入的是面，输出的是1*m个cell，m代表面的个数，每个cell里是该面邻接面的索引
function fring = compute_face_ring(face)
% compute_face_ring - compute the 1 ring of each face in a triangulation.
%
% fring = compute_face_ring(face);
%
% fring{i} is the set of faces that are adjacent
% to face i.
%
% Copyright (c) 2004 Gabriel Peyre
% the code assumes that faces is of size (3,nface)
[~,face] = check_face_vertex([],face);
nface = size(face,2);
nvert = max(max(face));
A = compute_edge_face_ring(face);
[~,~,s1] = find(A); % direct link
[i,j,s2] = find(A'); % reverse link
I = find(i<j);
s1 = s1(I); s2 = s2(I);
fring{nface} = [];
for k=1:length(s1)
if s1(k)>0 && s2(k)>0
fring{s1(k)}(end+1) = s2(k);
fring{s2(k)}(end+1) = s1(k);
end
end
return;
for i=1:nface
face(i,:) = sort(face(i,:));
end
edges = zeros(3*nface,2);
edges( 1:nface, : ) = face(:,1:2);
edges( (nface+1):2*nface, : ) = face(:,2:3);
edges( (2*nface+1):3*nface, : ) = face(:,[1,3]);
for i=1:nface
fring{i} = [];
end
for i=1:3*nface
e = edges(i,:);
I = find(edges(:,1)==e(1) & edges(:,2)==e(2));
if length(I)==2
f1 = mod( I(1)-1, nface)+1;
f2 = mod( I(2)-1, nface)+1;
fring{f1} = [fring{f1}, f2];
fring{f2} = [fring{f2}, f1];
edges(I,:) = rand(length(I),2);
end
end
return;
return;
h = waitbar(0,'Computing Face 1-ring');
for i=1:nface
waitbar(i/nface);
fring{i} = [];
for j=1:3
j1 = j;
j2 = mod(j,3)+1;
v1= face(i,j1);
v2= face(i,j2);
% test if another face share the same vertices
f = [];
for i1=1:nface
for a=1:3
if face(i1,a)==v1 && i1~=i
for b=1:3
if face(i1,b)==v2
% add to the ring
fring{i} = [fring{i}, i1];
end
end
end
end
end
end
end
close(h);
计算二面角的函数compute_adface_angle.m函数
输出值是每个面的flat度和邻接面是否是flat的判断cell，每个cell里是0和1。0表示对应的邻接面是flat的，红色的30可以修改，表示阈值，小于30度标记为flat
function [sum_num_of_flat_faces, num_of_flat_faces]=compute_adface_angle(vertex,face)
%COMPUTE_ADFACE_ANGLE Summary of this function goes here
% Detailed explanation goes here
%the first face is made of 1th and 2nd points
%the second face is made of 1th and 3rd points
%the third face is made of 2nd and 3rd face
fring=compute_face_ring(face);
normalf=compute_face_normal(vertex,face);
nface = size(face,2);
face_angle_cos{nface} = [];
face_angle{nface}=[];
sum_num_of_flat_faces=zeros(1,nface);
num_of_flat_faces{nface}=[];
for i=1:nface
for j=1:3
face_angle_cos{i}(1,j)=dot(normalf(:,i),normalf(:,fring{i}(1,j)));
%compute the angle between the faces
face_angle{i}(1,j)=rad2deg(acos(face_angle_cos{i}(1,j)));
if face_angle{i}(1,j)>30
num_of_flat_faces{i}(1,j)=1;
sum_num_of_flat_faces(i)=sum_num_of_flat_faces(i)+1;
else
num_of_flat_faces{i}(1,j)=0;
end
end
end
end
还有一个辅助函数check_face_vertex.m
function [vertex,face] = check_face_vertex(vertex,face)
% check_face_vertex - check that vertices and faces have the correct size
%
% [vertex,face] = check_face_vertex(vertex,face);
%
% Copyright (c) 2007 Gabriel Peyre
vertex = check_size(vertex,2,4);
face = check_size(face,3,4);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function a = check_size(a,vmin,vmax)
if isempty(a)
return;
end
if size(a,1)>size(a,2)
a = a';
end
if size(a,1)<3 && size(a,2)==3
a = a';
end
if size(a,1)<=3 && size(a,2)>=3 && sum(abs(a(:,3)))==0
% for flat triangles
a = a';
end
if size(a,1)<vmin || size(a,1)>vmax
error('face or vertex is not of correct size');
end
end
最后就是细分函数啦
function [newVertices, newFaces] = linearSubdivision(vertices, faces)
% Mesh subdivision using the Loop scheme.
%
% Dimensions:
% vertices: 3xnVertices
% faces: 3xnFaces
%
% Author: Jesus Mena
global edgeVertex;
global newIndexOfVertices;
newFaces = [];
newVertices = vertices;
nVertices = size(vertices,2);
nFaces = size(faces,2);
edgeVertex = zeros(nVertices, nVertices);
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 = addEdgeVertex(vaIndex, vbIndex);
vqIndex = addEdgeVertex(vbIndex, vcIndex);
vrIndex = addEdgeVertex(vaIndex, vcIndex);
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 = addEdgeVertex(vbIndex, vcIndex);
vrIndex = addEdgeVertex(vaIndex, vcIndex);
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 = addEdgeVertex(vaIndex, vbIndex);
vrIndex = addEdgeVertex(vaIndex, vcIndex);
threeFaces = [vaIndex,vpIndex,vrIndex; vpIndex,vbIndex,vcIndex; vrIndex,vpIndex,vcIndex]';
else
vpIndex = addEdgeVertex(vaIndex, vbIndex);
vqIndex = addEdgeVertex(vbIndex, vcIndex);
threeFaces = [vaIndex,vpIndex,vqIndex; vpIndex,vbIndex,vqIndex; vaIndex,vqIndex,vcIndex]';
end
newFaces = [newFaces, threeFaces];
elseif valence(i)==1
%
% * 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 flat
% vd------vb-----ve vd------vb------ve
% \ / \ / \ / |\ /
% \ / \ / \ / | \ /
% \ / \ / \ / | \ /
% va ----- vc -> va-- vr --vc
% \ / \ /
% \ / \ /
% \ / \ /
% vf vf
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 = addEdgeVertex(vaIndex, vbIndex);
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 = addEdgeVertex(vbIndex, vcIndex);
twoFaces = [vaIndex,vbIndex,vqIndex; vaIndex,vqIndex,vcIndex]';
else
vrIndex = addEdgeVertex(vaIndex, vcIndex);
twoFaces = [vaIndex,vbIndex,vrIndex; vrIndex,vbIndex,vcIndex]';
end
newFaces = [newFaces, twoFaces];
else
oneFace=[vaIndex, vbIndex, vcIndex]';
newFaces = [newFaces, oneFace];
end
end;
% ------------------------------------------------------------------------ %
% positions of the new vertices
for v1=1:nVertices-1
for v2=v1:nVertices
vNIndex = edgeVertex(v1,v2);
if (vNIndex>0)
newVertices(:,vNIndex) = 1/2*(vertices(:,v1)+vertices(:,v2));
end;
end;
end;
end
function vNIndex = addEdgeVertex(v1Index, v2Index)
global edgeVertex;
global newIndexOfVertices;
if (v1Index>v2Index) % setting: v1 <= v2
vTmp = v1Index;
v1Index = v2Index;
v2Index = vTmp;
end;
if (edgeVertex(v1Index, v2Index)==0) % new vertex
newIndexOfVertices = newIndexOfVertices+1;
edgeVertex(v1Index, v2Index) = newIndexOfVertices;
end;
vNIndex = edgeVertex(v1Index, v2Index);
return;
end实现效果
再对之前的200个点的兔子做细分
分别设置阈值为0,30,50,180，也就是二面角小于0,30,50,180度设置为flat的，一次细分得到的分别为
 |  |  |  |
| 0 | 30 | 50 | 180 |
可以很明显的看出来阈值不同对不同地方细分的程度是不同的
扫一扫
4337
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
