三维网格精简算法(Quadric Error Metrics)附源码(一)
来源:互联网 发布:宁波行知小学总课表 编辑:程序博客网 时间:2024/03/28 20:34
本文转自:http://www.cnblogs.com/shushen/p/5311828.html
我发现他的一系列文章都挺好,就是总缺点东西,所以没法执行
本文的算法来源于 Michael Garland在97年的文章
Surface simplification using quadric error metrics
算法介绍
在计算机图形应用中,为了尽可能真实呈现虚拟物体,往往需要高精度的三维模型。然而,模型的复杂性直接关系到它的计算成本,因此高精度的模型在几何运算时并不是必须的,取而代之的是一个相对简化的三维模型,那么如何自动计算生成这些三维简化模型就是网格精简算法所关注的目标。
[Garland et al. 1997]提出了一种基于二次误差作为度量代价的边收缩算法,其计算速度快并且简化质量较高。该方法在选择一条合适的边进行迭代收缩时,定义了一个描述边收缩代价的变量
其中
- 对所有的初始顶点计算
Q 矩阵.- 选择所有有效的边(这里取的是联通的边,也可以将距离小于一个阈值的边归为有效边)
- 对每一条有效边
(v1,v2) ,计算最优抽取目标v¯ .误差v¯T(Q1+Q2)v¯ 是抽取这条边的代价(cost)- 将所有的边按照cost的权值放到一个堆里
- 每次移除代价(cost)最小的边,并且更新包含着
v1 的所有有效边的代价
剩下的问题就是如何计算每个顶点的初始误差矩阵
其中
因此原始网格中顶点v的初始误差为
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 numbernp = percent*nv; % remained vertex numberQ0 = 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:nffor j = 1:3valence(F(i,j)) = valence(F(i,j)) + 1;Q(:,:,F(i,j)) = Q(:,:,F(i,j)) + Q0(:,:,i);endendTR = triangulation(F,V);E = edges(TR);% compute Q1+Q2 for each pairQbar = Q(:,:,E(:,1)) + Q(:,:,E(:,2));% a simple scheme: select either v1, v2 or (v1+v2)/2ne = 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;ticfor 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 testname = 'bunny_200.obj';OBJ=readObj(name);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 numbernp = 0.1*nv; % remained vertex numberQ0 = 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:nffor j = 1:3valence(F(i,j)) = valence(F(i,j)) + 1;Q(:,:,F(i,j)) = Q(:,:,F(i,j)) + Q0(:,:,i);endendTR = triangulation(F,V);E = edges(TR);% compute Q1+Q2 for each pairQbar = Q(:,:,E(:,1)) + Q(:,:,E(:,2));% a simple scheme: select either v1, v2 or (v1+v2)/2ne = 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;ticfor 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 cameramenu 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)';endend
最终得到原来点的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 typesv = []; 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]; endendfclose(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 facesnormalf = 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 vertexnormal = zeros(3,nvert);for i=1:nface f = face(:,i); for j=1:3 normal(:,f(j)) = normal(:,f(j)) + normalf(:,i); endend% normalized = sqrt( sum(normal.^2,1) ); d(d<eps)=1;normal = normal ./ repmat( d, 3,1 );% enforce that the normal are outwardv = 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) dimensionalz = 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 Peyrevertex = check_size(vertex,2,4);face = check_size(face,3,4);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function a = check_size(a,vmin,vmax)if isempty(a) return;endif size(a,1)>size(a,2) a = a';endif size(a,1)<3 && size(a,2)==3 a = a';endif size(a,1)<=3 && size(a,2)>=3 && sum(abs(a(:,3)))==0 % for flat triangles a = a';endif 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*3nV=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;endrecF=zeros(nF,3);for i=1:nF for j=1:3 recF(i,j)=F(i,j)-num_of_NaN(F(i,j)); endendrecV=zeros(nV-sum,3);j=1;for i=1:nV if ~isnan(V(i,1)) recV(j,:)=V(i,:); j=j+1; endendend
实现效果
我这里将原来有三万多个点的兔子模型
~~ 后来官方把我的模型代码删除了,不知道为什么~~
用里面那个cow.ply也是一样的做简化,得到的效果为
- 三维网格精简算法(Quadric Error Metrics)附源码(一)
- 三维网格精简算法(Quadric Error Metrics)附源码(二)
- 网格细分算法(Catmull-Clark subdivision & Loop subdivision)附源码
- 网格形路径(算法源码)
- 网格聚类算法(一)
- Java服务器性能监控(一) Metrics
- jre精简详细教程(附精简工具)
- jre精简详细教程(附精简工具)
- GL游戏算法(附源码)
- GL 游戏算法(附源码)
- 浅析KMP算法(附C++源码)
- SM4密码算法(附源码)
- SM4密码算法(附源码)
- SM4密码算法(附源码)
- tinyweb: C语言 + libuv 开发的最精简的WebServer (附源码)
- tinyweb: C语言 + libuv 开发的最精简的WebServer (附源码)
- Keras/Python深度学习中的网格搜索超参数调优(附源码)
- Keras/Python深度学习中的网格搜索超参数调优(附源码)
- 软件光栅化渲染器(九)
- Windows上部署Jenkins遇到的问题
- byte为何范围是-128~127
- 编程第五十六天
- 【Leetcode】389. Find the Difference
- 三维网格精简算法(Quadric Error Metrics)附源码(一)
- 双队列=>栈
- 记录一下这几天eclipse建立maven工程遇到的问题
- Tensorflow学习笔记--使用迁移学习做自己的图像分类器(Inception v3)
- 1.2.1.1 GNU Debugger
- 从LINQ开始之LINQ to Objects(上)
- 1.2.1.1 GNU调试器
- Problem A: 你会定义类吗?
- 基于主分量和支持向量机的人脸识别