三维网格精简算法(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]提出了一种基于二次误差作为度量代价的边收缩算法,其计算速度快并且简化质量较高。该方法在选择一条合适的边进行迭代收缩时,定义了一个描述边收缩代价的变量Δ,具体如下:对于网格中的每个顶点v,我们预先定义一个4×4的对称误差矩阵Q,那么顶点v=[vx vy vz 1]T的误差为其二次项形式Δ(v)=vTQv。假设对于一条收缩边(v1,v2),其收缩后顶点变为vbar,我们定义顶点vbar的误差矩阵QbarQbar=Q1+Q2,对于如何计算顶点vbar的位置有两种策略:一种简单的策略就是在v1,v2(v1+v2)/2中选择一个使得收缩代价Δ(vbar)最小的位置。另一种策略就是数值计算顶点vbar位置使得Δ(vbar)最小,由于Δ的表达式是一个二次项形式,因此令一阶导数为0,即,该式等价于求解:

q11q12q130q12q22q230q13q23q330q14q24q341v¯=0001

其中qij为矩阵Qbar中对应的元素。如果系数矩阵可逆,那么通过求解上述方程就可以得到新顶点vbar的位置,如果系数矩阵不可逆,就通过第一种简单策略来得到新顶点vbar的位置。根据以上描述,算法流程如下:

  1. 对所有的初始顶点计算Q矩阵.
  2. 选择所有有效的边(这里取的是联通的边,也可以将距离小于一个阈值的边归为有效边)
  3. 对每一条有效边(v1,v2),计算最优抽取目标v¯.误差v¯T(Q1+Q2)v¯是抽取这条边的代价(cost)
  4. 将所有的边按照cost的权值放到一个堆里
  5. 每次移除代价(cost)最小的边,并且更新包含着v1的所有有效边的代价

剩下的问题就是如何计算每个顶点的初始误差矩阵Q,在原始网格模型中,每个顶点可以认为是其周围三角片所在平面的交集,也就是这些平面的交点就是顶点位置,我们定义顶点的误差为顶点到这些平面的距离平方和:

Δ(v)=Δ([vx vy vz 1]T)=pplans(v)(pTv)2=pplans(v)(vTp)(pTv)=pplans(v)vT(ppT)v=vT(pplans(v)Kp)v

其中p=[a b c d]T代表平面方程ax+by+cz+d=0(a2+b2+c2=1)的系数,Kp为二次基本误差矩阵:

Kp=ppTa2abacadabb2bcbdacbcc2cdadbdcdd2

因此原始网格中顶点v的初始误差为Δ(v)=0,当边收缩后,新顶点误差为Δ(vbar)=vTbarQbarvbar,我们依次选取收缩后新顶点误差最小的边进行迭代收缩直到满足要求为止。

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也是一样的做简化,得到的效果为
这里写图片描述

原创粉丝点击