• 最近在研究网格简化，之前转载的一篇博文利用的是QEM度量误差进行网格简化，具体的算法三维网格精简算法（Quadric Error Metrics）附源码（一）一文中说明的已经很清楚了。 在上一篇博文中边收缩的方式是选取三个...
最近在研究网格简化，之前转载的一篇博文利用的是QEM度量误差进行网格简化，具体的算法在三维网格精简算法（Quadric Error Metrics）附源码（一）一文中说明的已经很清楚了。
在上一篇博文中边收缩的方式是选取三个点，分别是边的端点和中点，计算其QEM，从中选择最小的，将边收缩为一点，即为上篇博文中策略一。本文将其扩展为策略二，即每条边都计算使其二次误差矩阵最小的点，如果得不到的话，再运用策略一，从端点和中点中选择二次误差最小的点。即下面的simplification_points_optimals.m函数
其余的函数在上篇博文中均有提到

function [ simV,simF ] = simplification_points_optimals( V,F,Npoints)
%SIMPLIFICATION Summary of this function goes here
%   Detailed explanation goes here
[N] = compute_face_normal(V,F);
N=N';
p = [N, -sum(N .* V(F(:,1),:), 2)];%the constant trem of face according to the point normal form equation
nv = size(V,1); % total vertex number
np = Npoints; % remained vertex number
Q0 = bsxfun(@times, permute(p, [2,3,1]), permute(p, [3,2,1]));

% compute the Q matrices for all the initial vertices.
nf = size(F,1);
Q = zeros(4,4,nv);
valence = zeros(nv,1);
for i = 1:nf
for j = 1:3
valence(F(i,j)) = valence(F(i,j)) + 1;
Q(:,:,F(i,j)) = Q(:,:,F(i,j)) + Q0(:,:,i);
end
end

TR = triangulation(F,V);
E = edges(TR);
% compute Q1+Q2 for each pair
Qbar = Q(:,:,E(:,1)) + Q(:,:,E(:,2));
% a simple scheme: select either v1, v2 or (v1+v2)/2
ne = size(E,1);

costop=1./zeros(ne,1);
vop=zeros(4,1,ne);
for i=1:ne
Qstmp=Qbar(:,:,i);
Qstmp(4,:)=[0 0 0 1];
if (det(Qstmp))
vop(:,:,i)=Qstmp\[0;0;0;1];
costop(i)=vop(:,:,i)'*Qbar(:,:,i)*vop(:,:,i);
else
v1 = permute([V(E(i,1),:),ones(1)], [2,3,1]);
v2 = permute([V(E(i,2),:),ones(1)], [2,3,1]);
vm = 0.5 .* (v1 + v2);
v = [v1, v2, vm];
cost = zeros(1,3);
cost(1) = v1'*Qbar(:,:,i)*v1;
cost(2) = v2'*Qbar(:,:,i)*v2;
cost(3) = vm'*Qbar(:,:,i)*vm;

[min_cost, vidx] = min(cost,[],2);
vop(:,:,i)=v(:,vidx);
costop(i)=min_cost;
end
end

for i = 1:nv-np

[~, k] = min(costop);
e = E(k,:);

% update position for v1
V(e(1),:) = vop(1:3, k)';
V(e(2),:) = NaN;

% update Q for v1
Q(:,:,e(1)) = Q(:,:,e(1)) + Q(:,:,e(2));
Q(:,:,e(2)) = NaN;

% updata face
F(F == e(2)) = e(1);
f_remove = sum(diff(sort(F,2),[],2) == 0, 2) > 0;
F(f_remove,:) = [];

% collapse and delete edge and related edge information
E(E == e(2)) = e(1);
E(k,:) = [];
costop(k) = [];
Qbar(:,:,k) = [];
vop(:,:,k) = [];

% delete duplicate edge and related edge information
[E,ia,ic] = unique(sort(E,2), 'rows'); %#ok<NASGU>
costop = costop(ia);
Qbar = Qbar(:,:,ia);
vop = vop(:,:,ia);

% pairs involving v1
pair = sum(E == e(1), 2) > 0;
npair = sum(pair);

% updata edge information
Qbar(:,:,pair) = Q(:,:,E(pair,1)) + Q(:,:,E(pair,2));

costopt=1./zeros(size(npair,1),1);
pair_vop=zeros(4,1,npair);
indn=0;
for ii=1:size(pair,1)
if pair(ii)
indn=indn+1;
Qstmp=Qbar(:,:,ii);
Qstmp(4,:)=[0 0 0 1];
if (det(Qstmp))
pair_vop(:,:,indn)=Qstmp\[0;0;0;1];
costopt(indn)=pair_vop(:,:,indn)'*Qbar(:,:,ii)*pair_vop(:,:,indn);
else
v1t=permute([V(E(ii,1),:),ones(1)], [2,3,1]);
v2t=permute([V(E(ii,2),:),ones(1)], [2,3,1]);
vmt = 0.5 .* (v1t + v2t);

v = [v1t, v2t, vmt];
cost = zeros(1,3);
cost(1) = v1t'*Qbar(:,:,ii)*v1t;
cost(2) = v2t'*Qbar(:,:,ii)*v2t;
cost(3) = vmt'*Qbar(:,:,ii)*vmt;

[min_cost, vidx] = min(cost,[],2);
pair_vop(:,:,indx)=v(:,vidx);
costopt(indx)=min_cost;
end
end
end

vop(:,:,pair) = pair_vop;
costop(pair) = costopt;

if(size(unique(sort(F,2),'rows'),1)~=size(F,1))
warning('the %d process wrong, process duplicate faces\n',i);
Fsort=sort(F,2);
[~,m,~]=unique(Fsort,'rows');
a=1:size(F,1);
b=setdiff(a,m);
[~,c]=intersect(Fsort,Fsort(b,:),'rows');
result=[b;c];
F(result,:)=[];
end

end
[ simV,simF ] = rectifyindex( V,F );

end

但是我在这里遇到了一个问题，就是在代码的倒数6-10行

    if(size(unique(sort(F,2),'rows'),1)~=size(F,1))
warning('the %d process wrong, process duplicate faces\n',i);
[~, ia] = unique(sort(F,2),'rows');
F=F(ia,:);
end

这段代码的主要作用是排除掉重面，首先判断是否有重面，如果有的话，排除掉。在这里我也想请教下有没有做网格简化的前辈或者同学，网格简化是否会出现重面。我自己随便想想应该没有重面啊。。但是最后的结果却出现了重面，因此我加了这句话。不过我试了两个例子，只有在进行到某一步的时候出现了重面。去除掉了之后就没有了。
下面这段是新增的，通过调试我发现问题出现在925步
因此跳出，输出网格
显示下一步要收缩的点位于下图中鼻子处，放大得到右图，可以看到拓扑结构

问题出在蓝色的圈中，因为连接两点的1邻域中出现了异常，也就是同时有三个点与要收缩的这条边相连接，蓝色区域类似于四面体的三个面结构。因此如果下一步再次去除这条边的话，将出现重面。
而且这个重面是独立于其他面的（非流形面）。因此这两个面都应该去除

因此后来把倒数6-10行代码改为

    if(size(unique(sort(F,2),'rows'),1)~=size(F,1))
warning('the %d process wrong, process duplicate faces\n',i);
Fsort=sort(F,2);
[~,m,~]=unique(Fsort,'rows');
a=1:size(F,1);
b=setdiff(a,m);
[~,c]=intersect(Fsort,Fsort(b,:),'rows');
result=[b;c];
F(result,:)=[];
end把出现的两个面全部删去，同样的，在上万次的简化操作中，也只是出现了一次这样的异常。不过之前的问题演化成了——这种奇异的结构是否会出现？感觉这个还是有可能的嘿嘿，不过本人没有深入研究过，如果有问题，欢迎批评指正哈

运行的时候执行

name = 'MaxPlanck.obj';
V=OBJ.v;
F=OBJ.f.v;

[simV,simF] = simplification_points_optimals( V,F,2500);
%[simV,simF] = simplification( V,F,2500);
obj_write('MaxPlanck_2500.obj',simV',simF');得到简化网格

原始网格简化网格

展开全文
• 我发现他的一系列文章都挺好，就是总缺点... 本文的算法来源于 Michael Garland在97年的文章 Surface simplification using quadric error metrics 算法介绍在计算机图形应用中，为了尽可能真实呈现虚拟物体，往往
本文转自：http://www.cnblogs.com/shushen/p/5311828.html 我发现他的一系列文章都挺好，就是总缺点东西，所以没法执行

本文的算法来源于 Michael Garland在97年的文章 Surface simplification using quadric error metrics

算法介绍
在计算机图形应用中，为了尽可能真实呈现虚拟物体，往往需要高精度的三维模型。然而，模型的复杂性直接关系到它的计算成本，因此高精度的模型在几何运算时并不是必须的，取而代之的是一个相对简化的三维模型，那么如何自动计算生成这些三维简化模型就是网格精简算法所关注的目标。
[Garland et al. 1997]提出了一种基于二次误差作为度量代价的边收缩算法，其计算速度快并且简化质量较高。该方法在选择一条合适的边进行迭代收缩时，定义了一个描述边收缩代价的变量

Δ

\Delta

，具体如下：对于网格中的每个顶点

v

v

，我们预先定义一个4×4的对称误差矩阵

Q

Q

，那么顶点

v

=

[

v

x

v

y

v

z

1

]

T

v = [v_x\ v_y\ v_z\ 1]^T

的误差为其二次项形式

Δ

(

v

)

=

v

T

Q

v

\Delta(v) = v^TQv

。假设对于一条收缩边

(

v

1

,

v

2

)

(v_1, v_2)

，其收缩后顶点变为

v

b

a

r

v_{bar}

，我们定义顶点

v

b

a

r

v_{bar}

的误差矩阵

Q

b

a

r

Q_{bar}

为

Q

b

a

r

=

Q

1

+

Q

2

Q_{bar} = Q_1 + Q_2

，对于如何计算顶点

v

b

a

r

v_bar

的位置有两种策略：一种简单的策略就是在

v

1

,

v

2

v_1, v_2

和

(

v

1

+

v

2

)

/

2

(v_1+ v_2)/2

中选择一个使得收缩代价

Δ

(

v

b

a

r

)

\Delta(v_{bar})

最小的位置。另一种策略就是数值计算顶点

v

b

a

r

v_{bar}

位置使得

Δ

(

v

b

a

r

)

\Delta(v_{bar})

最小，由于

Δ

\Delta

的表达式是一个二次项形式，因此令一阶导数为0，即，该式等价于求解：

[

q

11

q

12

q

13

q

14

q

12

q

22

q

23

q

24

q

13

q

23

q

33

q

34

0

0

0

1

]

v

ˉ

=

[

0

0

0

1

]

\left[ \begin{matrix} q_{11}&q_{12}&q_{13}&q_{14}\\ q_{12}&q_{22}&q_{23}&q_{24}\\ q_{13}&q_{23}&q_{33}&q_{34}\\ 0&0&0&1 \end{matrix} \right]\bar{v}=\left[ \begin{matrix} 0\\0\\0\\1 \end{matrix} \right]

其中

q

i

j

q_{ij}

为矩阵

Q

b

a

r

Q_{bar}

中对应的元素。如果系数矩阵可逆，那么通过求解上述方程就可以得到新顶点

v

b

a

r

v_{bar}

的位置，如果系数矩阵不可逆，就通过第一种简单策略来得到新顶点

v

b

a

r

v_{bar}

的位置。根据以上描述，算法流程如下：

对所有的初始顶点计算

Q

Q

矩阵.选择所有有效的边（这里取的是联通的边，也可以将距离小于一个阈值的边归为有效边）对每一条有效边

(

v

1

,

v

2

)

(v_1,v_2)

，计算最优抽取目标

v

ˉ

\bar{v}

.误差

v

ˉ

T

(

Q

1

+

Q

2

)

v

ˉ

\bar{v}^T(Q_1+Q_2)\bar{v}

是抽取这条边的代价（cost）将所有的边按照cost的权值放到一个堆里每次移除代价（cost）最小的边，并且更新包含着

v

1

v_1

的所有有效边的代价

剩下的问题就是如何计算每个顶点的初始误差矩阵

Q

Q

，在原始网格模型中，每个顶点可以认为是其周围三角片所在平面的交集，也就是这些平面的交点就是顶点位置，我们定义顶点的误差为顶点到这些平面的距离平方和：

Δ

(

v

)

=

Δ

(

[

v

x

v

y

v

z

1

]

T

)

=

∑

p

∈

p

l

a

n

e

s

(

v

)

(

p

T

v

)

2

=

∑

p

∈

p

l

a

n

e

s

(

v

)

(

v

T

p

)

(

p

T

v

)

=

∑

p

∈

p

l

a

n

e

s

(

v

)

v

T

(

p

p

T

)

v

=

v

T

(

∑

p

∈

p

l

a

n

e

s

(

v

)

K

p

)

v

\Delta(v)=\Delta([v_x\ v_y\ v_z\ 1]^T)=\sum_{p\in planes(v)}(p^Tv)^2=\sum_{p\in planes(v)}(v^Tp)(p^Tv)=\sum_{p\in planes(v)}v^T(pp^T)v\\ =v^T(\sum_{p\in planes(v)}K_p)v

其中

p

=

[

a

b

c

d

]

T

p = [a\ b\ c\ d]^T

代表平面方程

a

x

+

b

y

+

c

z

+

d

=

0

(

a

2

+

b

2

+

c

2

=

1

)

ax + by + cz + d = 0(a^2 + b^2 + c^2 = 1)

的系数，

K

p

K_p

为二次基本误差矩阵：

K

p

=

p

p

T

[

a

2

a

b

a

c

a

d

a

b

b

2

b

c

b

d

a

c

b

c

c

2

c

d

a

d

b

d

c

d

d

2

]

K_p=pp^T \left[ \begin{matrix} a^2&ab&ac&ad\\ ab&b^2&bc&bd\\ ac&bc&c^2&cd\\ ad&bd&cd&d^2 \end{matrix} \right]

因此原始网格中顶点v的初始误差为

Δ

(

v

)

=

0

\Delta(v) = 0

，当边收缩后，新顶点误差为

Δ

(

v

b

a

r

)

=

v

b

a

r

T

Q

b

a

r

v

b

a

r

\Delta(v_{bar}) = v_{bar}^TQ_{bar}v_{bar}

，我们依次选取收缩后新顶点误差最小的边进行迭代收缩直到满足要求为止。
MATLAB代码
实现代码
V,F是输入的网格，percent代表简化率，simV和simF是输出网格
function [ simV,simF ] = simplification( V,F,percent )
%SIMPLIFICATION Summary of this function goes here
%   Detailed explanation goes here
[N] = compute_face_normal(V,F);
N=N';
p = [N, -sum(N .* V(F(:,1),:), 2)];
nv = size(V,1); % total vertex number
np = percent*nv; % remained vertex number
Q0 = bsxfun(@times, permute(p, [2,3,1]), permute(p, [3,2,1]));

% compute the Q matrices for all the initial vertices.
nf = size(F,1);
Q = zeros(4,4,nv);
valence = zeros(nv,1);
for i = 1:nf
for j = 1:3
valence(F(i,j)) = valence(F(i,j)) + 1;
Q(:,:,F(i,j)) = Q(:,:,F(i,j)) + Q0(:,:,i);
end
end

TR = triangulation(F,V);
E = edges(TR);
% compute Q1+Q2 for each pair
Qbar = Q(:,:,E(:,1)) + Q(:,:,E(:,2));
% a simple scheme: select either v1, v2 or (v1+v2)/2
ne = size(E,1);
v1 = permute([V(E(:,1),:),ones(ne,1)], [2,3,1]);
v2 = permute([V(E(:,2),:),ones(ne,1)], [2,3,1]);
vm = 0.5 .* (v1 + v2);
v = [v1, v2, vm];
cost = zeros(ne,3);
cost(:,1) = sum(squeeze(sum(bsxfun(@times,v1,Qbar),1)).*squeeze(v1),1)';
cost(:,2) = sum(squeeze(sum(bsxfun(@times,v2,Qbar),1)).*squeeze(v2),1)';
cost(:,3) = sum(squeeze(sum(bsxfun(@times,vm,Qbar),1)).*squeeze(vm),1)';

num = nv;
tic
for i = 1:nv-np

[min_cost, vidx] = min(cost,[],2);
[~, k] = min(min_cost);
e = E(k,:);

% update position for v1
V(e(1),:) = v(1:3, vidx(k), k)';
V(e(2),:) = NaN;

% update Q for v1
Q(:,:,e(1)) = Q(:,:,e(1)) + Q(:,:,e(2));
Q(:,:,e(2)) = NaN;

% updata face
F(F == e(2)) = e(1);
f_remove = sum(diff(sort(F,2),[],2) == 0, 2) > 0;
F(f_remove,:) = [];

% collapse and delete edge and related edge information
E(E == e(2)) = e(1);
E(k,:) = [];
cost(k,:) = [];
Qbar(:,:,k) = [];
v(:,:,k) = [];

% delete duplicate edge and related edge information
[E,ia,ic] = unique(sort(E,2), 'rows'); %#ok<NASGU>
cost = cost(ia,:);
Qbar = Qbar(:,:,ia);
v = v(:,:,ia);

% pairs involving v1
pair = sum(E == e(1), 2) > 0;
npair = sum(pair);

% updata edge information
Qbar(:,:,pair) = Q(:,:,E(pair,1)) + Q(:,:,E(pair,2));

pair_v1 = permute([V(E(pair,1),:),ones(npair,1)], [2,3,1]);
pair_v2 = permute([V(E(pair,2),:),ones(npair,1)], [2,3,1]);
pair_vm = 0.5 .* (pair_v1 + pair_v2);
v(:,:,pair) = [pair_v1, pair_v2, pair_vm];

cost(pair,1) = sum(squeeze(sum(bsxfun(@times,pair_v1,Qbar(:,:,pair)),1)).*squeeze(pair_v1),1)';
cost(pair,2) = sum(squeeze(sum(bsxfun(@times,pair_v2,Qbar(:,:,pair)),1)).*squeeze(pair_v2),1)';
cost(pair,3) = sum(squeeze(sum(bsxfun(@times,pair_vm,Qbar(:,:,pair)),1)).*squeeze(pair_vm),1)';

%fprintf('%d\n', i);
end
[ simV,simF ] = rectifyindex( V,F );

end

test.m测试代码
主要是我在调试程序的时候用到，在上面的程序中加入了一段，使得每减少一定的点就显示一遍并且重新写入一遍，为了方便观察-_-
function  test

name = 'bunny_200.obj';
V=OBJ.v;
F=OBJ.f.v;
% [ V, F ] = ply_to_tri_mesh('coww.ply');
% V=V';
% F=F';

[N] = compute_face_normal(V,F);
N=N';
p = [N, -sum(N .* V(F(:,1),:), 2)];
nv = size(V,1); % total vertex number
np = 0.1*nv; % remained vertex number
Q0 = bsxfun(@times, permute(p, [2,3,1]), permute(p, [3,2,1]));

% compute the Q matrices for all the initial vertices.
nf = size(F,1);
Q = zeros(4,4,nv);
valence = zeros(nv,1);
for i = 1:nf
for j = 1:3
valence(F(i,j)) = valence(F(i,j)) + 1;
Q(:,:,F(i,j)) = Q(:,:,F(i,j)) + Q0(:,:,i);
end
end

TR = triangulation(F,V);
E = edges(TR);
% compute Q1+Q2 for each pair
Qbar = Q(:,:,E(:,1)) + Q(:,:,E(:,2));
% a simple scheme: select either v1, v2 or (v1+v2)/2
ne = size(E,1);
v1 = permute([V(E(:,1),:),ones(ne,1)], [2,3,1]);
v2 = permute([V(E(:,2),:),ones(ne,1)], [2,3,1]);
vm = 0.5 .* (v1 + v2);
v = [v1, v2, vm];
cost = zeros(ne,3);
cost(:,1) = sum(squeeze(sum(bsxfun(@times,v1,Qbar),1)).*squeeze(v1),1)';
cost(:,2) = sum(squeeze(sum(bsxfun(@times,v2,Qbar),1)).*squeeze(v2),1)';
cost(:,3) = sum(squeeze(sum(bsxfun(@times,vm,Qbar),1)).*squeeze(vm),1)';

num = nv;
tic
for i = 1:nv-np
if (nv - i) < 0.9*num
num = nv - i;

clf
trimesh(F, V(:,1), V(:,2), V(:,3),'LineWidth',1,'EdgeColor','k');
%drawMesh(V, F, 'facecolor','y', 'edgecolor','k', 'linewidth', 1.2);
view([0 90])
axis equal
axis off
camlight
lighting gouraud
drawnow
end

[min_cost, vidx] = min(cost,[],2);
[~, k] = min(min_cost);
e = E(k,:);

% update position for v1
V(e(1),:) = v(1:3, vidx(k), k)';
V(e(2),:) = NaN;

% update Q for v1
Q(:,:,e(1)) = Q(:,:,e(1)) + Q(:,:,e(2));
Q(:,:,e(2)) = NaN;

% updata face
F(F == e(2)) = e(1);
f_remove = sum(diff(sort(F,2),[],2) == 0, 2) > 0;
F(f_remove,:) = [];

% collapse and delete edge and related edge information
E(E == e(2)) = e(1);
E(k,:) = [];
cost(k,:) = [];
Qbar(:,:,k) = [];
v(:,:,k) = [];

% delete duplicate edge and related edge information
[E,ia,ic] = unique(sort(E,2), 'rows'); %#ok<NASGU>
cost = cost(ia,:);
Qbar = Qbar(:,:,ia);
v = v(:,:,ia);

% pairs involving v1
pair = sum(E == e(1), 2) > 0;
npair = sum(pair);

% updata edge information
Qbar(:,:,pair) = Q(:,:,E(pair,1)) + Q(:,:,E(pair,2));

pair_v1 = permute([V(E(pair,1),:),ones(npair,1)], [2,3,1]);
pair_v2 = permute([V(E(pair,2),:),ones(npair,1)], [2,3,1]);
pair_vm = 0.5 .* (pair_v1 + pair_v2);
v(:,:,pair) = [pair_v1, pair_v2, pair_vm];

cost(pair,1) = sum(squeeze(sum(bsxfun(@times,pair_v1,Qbar(:,:,pair)),1)).*squeeze(pair_v1),1)';
cost(pair,2) = sum(squeeze(sum(bsxfun(@times,pair_v2,Qbar(:,:,pair)),1)).*squeeze(pair_v2),1)';
cost(pair,3) = sum(squeeze(sum(bsxfun(@times,pair_vm,Qbar(:,:,pair)),1)).*squeeze(pair_vm),1)';

end

end


最终得到原来点的1/10 的图
补充代码
读取ply文件来自http://people.sc.fsu.edu/~jburkardt/m_src/ply_io/ply_to_tri_mesh.m 读取obj的我之前的博文MATLAB读取和显示obj文件的数据有，这里这个是简化版的
function obj = readObj(fname)
%
% obj = readObj(fname)
%
% This function parses wavefront object data
% It reads the mesh vertices, texture coordinates, normal coordinates
% and face definitions(grouped by number of vertices) in a .obj file
%
%
% INPUT: fname - wavefront object file full path
%
% OUTPUT: obj.v - mesh vertices
%       : obj.vt - texture coordinates
%       : obj.vn - normal coordinates
%       : obj.f - face definition assuming faces are made of of 3 vertices
%
% Bernard Abayowa, Tec^Edge
% 11/8/07

% set up field types
v = []; vt = []; vn = []; f.v = []; f.vt = []; f.vn = [];

fid = fopen(fname);

% parse .obj file
while 1
tline = fgetl(fid);
if ~ischar(tline),   break,   end  % exit at end of file
ln = sscanf(tline,'%s',1); % line type
%disp(ln)
switch ln
case 'v'   % mesh vertexs
v = [v; sscanf(tline(2:end),'%f')'];
case 'vt'  % texture coordinate
vt = [vt; sscanf(tline(3:end),'%f')'];
case 'vn'  % normal coordinate
vn = [vn; sscanf(tline(3:end),'%f')'];
case 'f'   % face definition
fv = []; fvt = []; fvn = [];
str = textscan(tline(2:end),'%s'); str = str{1};

nf = length(findstr(str{1},'/')); % number of fields with this face vertices

[tok str] = strtok(str,'//');     % vertex only
for k = 1:length(tok) fv = [fv str2num(tok{k})]; end

if (nf > 0)
[tok str] = strtok(str,'//');   % add texture coordinates
for k = 1:length(tok) fvt = [fvt str2num(tok{k})]; end
end
if (nf > 1)
[tok str] = strtok(str,'//');   % add normal coordinates
for k = 1:length(tok) fvn = [fvn str2num(tok{k})]; end
end
f.v = [f.v; fv]; f.vt = [f.vt; fvt]; f.vn = [f.vn; fvn];
end
end
fclose(fid);

% set up matlab object
obj.v = v; obj.vt = vt; obj.vn = vn; obj.f = f;

此外计算面法向的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);

nface = size(face,2);
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);
for i=1:nface
f = face(:,i);
for j=1:3
normal(:,f(j)) = normal(:,f(j)) + normalf(:,i);
end
end
% normalize
d = sqrt( sum(normal.^2,1) ); d(d<eps)=1;
normal = normal ./ repmat( d, 3,1 );

% 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
normal = -normal;
normalf = -normalf;
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,:);

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

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


以及提到的函数rectifyindex.m，主要是用来去除掉网格中明显不对的点（坐标跑到无穷处的点）
function [ recV,recF ] = rectifyindex( V,F )
%RECTIFYINDEX Summary of this function goes here
%   V is nV*3
%   F is nF*3

nV=size(V,1);
nF=size(F,1);

num_of_NaN=zeros(nV,1);
sum=0;
for i=1:nV
if isnan(V(i,1))
sum=sum+1;
end
num_of_NaN(i)=sum;
end

recF=zeros(nF,3);

for i=1:nF
for j=1:3
recF(i,j)=F(i,j)-num_of_NaN(F(i,j));
end
end

recV=zeros(nV-sum,3);
j=1;
for i=1:nV
if ~isnan(V(i,1))
recV(j,:)=V(i,:);
j=j+1;
end
end

end

实现效果
我这里将原来有三万多个点的兔子模型 ~~后来官方把我的模型代码删除了，不知道为什么 ~~ 用里面那个cow.ply也是一样的做简化，得到的效果为 
展开全文
• 三维模型网格简化源码

热门讨论 2014-07-10 12:09:44
This program implements four different mesh simplification algorithms. After loading a mesh, the user can easily remove triangles from the mesh and the results are displayed in real time....
• 三维网格精简算法QEM对原始边的缩减的实现，在Qt中实现。
• 基于八叉树的网格简化算法实现 能正常的运行和对网格进行简化
• Qslim简化流程Qslim简化算法流程大体上可以概括为如下个步骤： 模型的读取 模型顶点Quadrics值的采集 根据顶点Quadrics值进行排序与构造顶点对栈 从顶点对栈中移除顶点对 首先从模型的提取开始说起模型的提取Qslim...
Qslim简化流程
Qslim简化算法流程大体上可以概括为如下三个步骤：
首先从模型的提取开始说起
模型的提取
接上步骤，QSlim程序读取完数据并保存到MxStdModel类后，将会启动Slim类。Slim类分为两种：FSlim与ESlim，分别对应面简化与边简化。他们的差别也会在之后详细解释。
展开全文
• 网格简化 二 、QEM算法

千次阅读 2020-10-16 01:52:46
简化算法的误差测度（度量质量和误差） 误差测度用于度量模型简化的质量和误差，因此它对模型的简化过程和最后的简化结果都具有重要的影响。大多数简化算法采用对象空间(Object-space)的一种或综合几种形式的几何...
简化算法的误差测度（度量质量和误差）
误差测度用于度量模型简化的质量和误差，因此它对模型的简化过程和最后的简化结果都具有重要的影响。大多数简化算法采用对象空间(Object-space)的一种或综合几种形式的几何误差(Geometric errors)作为误差测度，一些视点相关算法通常将对象空间的误差转换为屏幕空间(Screen-space)的误差值为误差测度，有些算法也考虑模型的颜色、法向量和纹理坐标等属性误差(Attribute errors)。几何误差 几何误差测度一般采用欧式空间距离表示。通常有顶点到顶点、顶点到平面和平面到平面的距离等形式。 Hausdorff距离是现有算法中常常用到度量顶点到表面、表面到表面距离的几何误差测度，该距离为两个模型的顶点之间的最小距离中的最大值。 给定欧式空间的两点集 ，Haousdorff距离就是用来衡量这两个点集间的距离。    算法过程：
该算法的时间复杂度是O(n，m)，其中n和m分别为集合A和集合B中的点数。屏幕空间误差计算 视点相关算法常常需要将对象空间误差转换为屏幕空间误差。设对象空间几何误差为e，x为以像素表示的屏幕某方向的分辨率，d是视点到模型对象的距离，θ为视野夹角，则e对应的屏幕空间误差p为：
属性误差（材质、纹理）
网格模型上的三角面片、、法向量、纹理坐标、顶点的颜色是其常见的属性。 网格模型的颜色一般以(r, g, b)三元组形式来表示，各分量分别在[0, 1]中取值。最直接的方法是采用欧式空间距离求解方法来求取颜色的距离。设简化过程的两模型M1、M2的颜色分别表示为(r1, g1, b1)和(r2, g2, b2)，则两模型的颜色距离dc可以表示为：
两个法向量的误差距离dn通常采用角度值进行度量：
多边形表面的纹理坐标用(u, v)坐标对来表示网格模型顶点到二维纹理空间的映射，其中，u，v通常在[0, 1]中取值。一般也是采用欧式空间距离求解方法来计算纹理坐标误差：
简化算法的约束条件或运行条件
模型简化过程中或简化算法运行时往往存在一些限制条件，这些条件也决定了算法采用的技术、算法运行效果和模型简化结构等。细节层次（LOD） 对于各种简化细节层次的LOD模型的管理技术可以分成离散LOD(Discrete LOD)、连续LOD(Continuous LOD)和视点相关LOD技术。 早期简化算法大多采用离散LOD技术。这种技术首先采用离线(offline)方式对原始模型进行预处理，生成一系列不同分辨率的简化模型。在实际运行时，根据需要选择已生成的某个简化模型进行绘制。由于在实时显示绘制时不需要再次进行简化操作，因此该技术具有实时运行速度快、数据存储结构简单等优点。但是因为需要保存多个预处理的中间简化模型，所以占用存储空间大；且在简化预处理时无法考虑视点及实时运行环境因素等要求，只能根据模型本身信息进行简化，因而简化效率不高；同时由于技术限制，预简化生成LOD模型数量不可能过多，粒度不可能太细，因此实时显示绘制时候，在不同简化模型切换过程中会出现画面跳跃、视觉不连续等的效果。离散LOD简化也称为静态简化。 连续LOD技术是对传统的离散LOD技术的改进和发展。与离散LOD技术不同，连续LOD技术的各简化模型不是在预处理中生成，而是通过构造特定的数据结构进行编码存储，在实时显示运行时根据需要生成对应细节层次的简化模型。因此连续LOD技术具有更高的LOD粒度表示，占用空间较小，运行时画面连续性较好等优点；但由于运行时需要进行简化模型的生成处理，因此实时显示速度收到一定影响。连续LOD技术支持多边形网格模型的传输，常被应用于网格模型的各种递进简化算法中。简化模型的拓扑结构保持 简化过程中是否保持网格模型的拓扑结构不变也是区分不同简化算法的一个重要依据。网格模型的拓扑结构通常指构成网格模型的各三角网格之间的连接关系。衡量模型简化算法是否能保持拓扑结构一般是通过判断网格表面的亏格(Genus)和流型(Manifold)是否在简化过程中保持不变来确定。亏格采用网格表面的孔洞数量来计算。 下图为两种类型的顶点对：
模型试验结果：

QEM算法的基本操作基于边折叠，误差测度采用的是二次误差测度。二次误差测度最早是由Garland提出，采用点到平面距离的平方作为误差测度。它的优点是具有较高的计算速度，较小的内存消耗，而且得到的简化网格具有较高质量。它是在速度非常快但简化质量很差、速度很慢但简化质量非常好的两类方法之间的一种折中，是一种兼顾了速度和质量的较理想的误差测度。 在三维欧氏空间中，一个平面可以表示为： 其中 时平面的单位法向量，d时常量。点 到该平面的距离就可以表示为：
表面属性
在计算机图形学中， 三角网格模型最常见的表面属性有颜色、纹理和法线。为了使简化模型同初始模型具有良好的相似性， 必须在保持模型几何信息的同时保留这些属性特征。由于点到平面的距离考虑了简化操作对顶点周围区域属性值变化的影响， 可以比较准确地描述局部属性误差， 同时又比点到表面或表面到表面的距离计算简便快捷。因此， 采用点到平面的距离作为属性误差测度， 将二次误差测度应用到属性误差的计算中。 网格模型的每个顶点除了空间坐标外，还具有描述其属性的数值。在网格模型的三角面上， 属性值根据几何位置插值得到。因此， 三角面上的属性值是连续的， 而且两个属性值之间的距离用欧氏距离来度量。 比如对于颜色属性，可以用三维矢量r,g,bT 来表示(0≤r,g,b≤1 )，所有颜色矢量构成了RGB彩色空间，在RGB彩色空间中点到平面的距离平方同样可以用二次误差Q(v)来计算。边折叠后的新顶点采用子集选择法，不用重新计算顶点的空间位置和属性值，在计算误差的时候不用考虑空间坐标和属性值的相关性，只需分别建立几何二次误差测度和属性二次误差测度，并计算几何和属性误差。
边折叠操作的代价
采用带有颜色属性的模型应用算法，带有其他属性的网格模型可以同理推出。三角网格模型的每个顶点vg=x,y,zT 和vc=r,g,bT (0≤r,g,b≤1 )来表征几何和颜色信息。为每个三角面建立几何二次误差测度Qfg 和颜色二次误差测度Qfc 。各顶点的二次误差测度之和：
当边折叠(v1 , v2 )到顶点v的时候，总的二次误差测度为：
故而边折叠引起的几何误差Eg=Qgvg ，颜色属性误差Ec=Qcvc 。则总的边折叠代价为：
其中α为颜色属性误差在总代价中的影响系数，可以根据实际情况进行调节
参考资料：
Garland Heckbert. 网格简化算法 网格简化 vcglib库 vcglib库文章
展开全文
• 针对三维虚拟场景的物理属性显示需求,提出一种带属性的边折叠的三角形网格简化方法.该算法计算折叠代价时以模型边曲率和边上物理属性的增量以及三角形正则度作为权因子,边上物理属性的增量使简化后的模型很好地保留...
• 网格简化算法研究 摘 要 多边形几何模型变得越来越复杂这无疑给三维物体的实时绘制带来不便 于是网格简化算法成为目前计算机图形学领域的重点研究之一它以算法简化 速度数据结构存储的有效性误差控制等作为衡量标准...
• 在边折叠的网格简化算法的基础上，针对特定三维网格模型—人脸，提出了一种实用的基于特征点的快速模型简化算法。该算法把人脸按特征点的分布进行分块处理，对不同区域采用不同的阈值进行调节。对于需要高细节的区域...
• 随着三维网格在数字游戏、影视动画、虚拟现实等领域的广泛应用,针对三维网格的处理方法也越来越多,包括压缩、简化、嵌入水印和去噪等,这些处理技术均不可避免地导致三维网格的失真。如何更好地评估三维网格的视觉...
• 针对目前三维模型水印算法在将水印嵌入三维模型后，均会使三维模型产生局部失真问题，文章提出一种自适应三维网格水印算法。通过特征点提取算法将模型的特征点提取出来，按其模长进行分组，将每组顶点的模长序列进行...
• 然而，模型的复杂性直接关系到它的计算成本，因此高精度的模型在几何运算时并不是必须的，取而代之的是一个相对简化三维模型，那么如何自动计算生成这些三维简化模型就是网格精简算法所关注的目标。[Garland et al...
• 这次介绍的是一种比较简单的网格简化算法，叫做顶点聚簇。   网格简化 　为了介绍这个算法，首先说明一下网格简化算法。随着计算机绘图在现代科技领域中的广泛应用, 计算机图形在现代制造业中发挥着重要...
• 对水平集网格简化算法和现常用的基于点对收缩的网格简化算法在视觉质量和几何误差方面做了比较和分析，实验表明该方法适用于任意拓扑形状的网格模型，使得模型大规模简化后，在保持较低误差的同时，仍然能够保持相当...
• 提出了一种基于离散曲率的二次误差度量网格简化算法。在代价函数中引入顶点离散曲率,通过将代价函数作为顶点对的权值来控制顶点对合并次序,更好地保留了原模型的细节特征,同时修改模型特征点与特征线的权值,使得简化...
• 基于边界采样和投影方法建立代理曲面,实现了内环三角化细分和三维网格表面的自由变形。通过合并环间区域和引入边缘约束的平滑过程,得到缺失实体孔洞的上缝合面。在此基础上向内扩展,获取孔洞下表面边界轮廓,搜索出...
• 为了解决这个问题，本文作者提出一种改进的基于二次误差测度的网格简化算法．通过对顶点进行分类，在简化过程中更好地保持了模型的细节特征，同时考虑了网格中三角面的分布情况．减小了几何误差．结果表明，算法既...
• 为满足大数据量地学模型可视化的功能需求，实现一种基于拓扑规则和地学规则双重约束的三角网格简化和还原算法，以支持层次细节模型(level of detail，LOD)功能。介绍了实现算法所需的多分辨率数据结构。利用三维地学...
• 公司项目需求所以最近来研究减面，这里主要是大概简单的简述减面的思维等等还有几种常见的算法
• 文中算法首先将三维网格模型进行仿射变换，以获得模型的旋转不变性、缩放不变性；然后，将空间坐标系原点到模型顶点的距离通过一维DCT变化到频域，在频域信号上加入水印，使得嵌入的水印具有不可见性。实验结果表明...
• 基于此,提出一种改进的三维点云重建算法。该方法首先用统计滤波器对点云简化去噪,消除重建表面的锯齿状现象;然后建立点云间拓扑结构,对点云法向量进行法向重定向,以减少法向指向的二义性;最后将具有磁盘拓扑结构的...
• 针对三维模型简化后的精度与效率难以平衡的问题进行研究，提出一种局部特征熵的半边折叠非均匀网格简化算法。采用两次局部区域聚类探测，首先探测三维数据点所在边聚类局部区域，获取该探测区域法向量；其次以三维...
• 三维重建算法综述|传统+深度学习

千次阅读 多人点赞 2020-02-22 13:14:08
来源：基于深度学习的三维重建算法综述 00 前言 01 基于传统多视图几何的三维重建算法 1.1 主动式 （1）结构光 （2）TOF 激光飞行时间法 （3）三角测距法 1.2 被动式 （1）单目视觉 （2）双目/多目视觉 1.3 基于...
• 网格模型简化算法综述 董方敏2 张2 蕊2 刘2 勇2 周学君 三峡大学电气信息学院湖北宜昌2 ! 摘要网格模型的简化是解决复杂三维模型存储传输和实时绘制的有效性与硬件处理能力的局 限性之间矛盾的主要方法1 在分析了...
• 中提出一种带属性的三角网格模型简化的算法, 该算法将表征网格顶点信息的向量由简单的三维几何信息扩展到 包含颜色、纹理等附加信息的多维向量, 在Rn空间中用对称Hausdo rff 距离控制网格简化顺序和精度, 既保证了...
• 实现边坍塌的网格简化方法

千次阅读 2018-02-04 20:47:01
关于边塌陷的网格简化方法，SIGGRAPH 有一篇97年的论文Surface Simplification Using Quadric Error Metrics（Michael&&PaulS）对这个问题进行了深入的探讨。作业里的代码就是基于这篇论文来实现的。最为核心的思想...

...