三維網格精簡算法（Quadric Error Metrics）附源碼

　　在計算機圖形應用中，爲了儘量真實呈現虛擬物體，每每須要高精度的三維模型。然而，模型的複雜性直接關係到它的計算成本，所以高精度的模型在幾何運算時並非必須的，取而代之的是一個相對簡化的三維模型，那麼如何自動計算生成這些三維簡化模型就是網格精簡算法所關注的目標。

　　[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\n', i);
end

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

View Code

 

 

本文爲原創，轉載請註明出處：http://www.cnblogs.com/shushen

 

 

參考文獻：

[1] Michael Garland and Paul S. Heckbert. 1997. Surface simplification using quadric error metrics. In Proceedings of the 24th annual conference on Computer graphics and interactive techniques (SIGGRAPH '97). ACM Press/Addison-Wesley Publishing Co., New York, NY, USA, 209-216.