三维网格精简算法（Quadric Error Metrics）附源码（转载）

转载：  https://www.cnblogs.com/shushen/p/5311828.html

在计算机图形应用中，为了尽可能真实呈现虚拟物体，往往需要高精度的三维模型。然而，模型的复杂性直接关系到它的计算成本，因此高精度的模型在几何运算时并不是必须的，取而代之的是一个相对简化的三维模型，那么如何自动计算生成这些三维简化模型就是网格精简算法所关注的目标。

　　[Garland et al. 1997]提出了一种基于二次误差作为度量代价的边收缩算法，其计算速度快并且简化质量较高。该方法在选择一条合适的边进行迭代收缩时，定义了一个描述边收缩代价的变量Δ，具体如下：对于网格中的每个顶点v，我们预先定义一个4×4的对称误差矩阵Q，那么顶点v = [v_x v_y v_z 1]^T的误差为其二次项形式Δ(v) = v^TQv。假设对于一条收缩边(v₁, v₂)，其收缩后顶点变为v_bar，我们定义顶点v_bar的误差矩阵Q_bar为Q_bar = Q₁ + Q₂，对于如何计算顶点v_bar的位置有两种策略：一种简单的策略就是在v₁, v₂和(v₁+ v₂)/2中选择一个使得收缩代价Δ(v_bar)最小的位置。另一种策略就是数值计算顶点v_bar位置使得Δ(v_bar)最小，由于Δ的表达式是一个二次项形式，因此令一阶导数为0，即，该式等价于求解：

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

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

其中p = [a b c d]^T代表平面方程ax + by + cz + d = 0(a² + b² + c² = 1)的系数，K_p为二次基本误差矩阵：

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

%% surface simplification using quadratic error metrics
clc
clear all
close all

load('bunny.mat');

V = surfdata.point;
F = surfdata.face;

nv = size(V,1); % total vertex number
np = 0.1*nv; % remained vertex number

% fundamental error quadric
N = normals(V,F);
N = normalizeVector3d(N);
p = [N, -sum(N .* V(F(:,1),:), 2)];
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

% all valid pairs
E = edges(F);

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

figure;
draw_t = 0;
num = nv;
tic
for i = 1:nv-np
    if (nv - i) < 0.9*num
        num = nv - i;
        
        clf
        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);
    [useless, 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');
    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
', i);
end

clf
drawMesh(V, F, 'facecolor','y', 'edgecolor','k', 'linewidth', 1.2);
view([0 90])
axis equal
axis off
camlight
lighting gouraud
cameramenu