QEM网格简化算法

欢迎关注更多精彩关注我，学习常用算法与数据结构，一题多解，降维打击。qem 全称Quadic Error Metrics 网格简化。它是一种基于二次度量误差的优化算法。边塌缩算法qem算法的基本思想是对某一边进行塌缩，将边的2点合成一个点，如下图。那么要选择哪条边进行塌缩呢。就是基于二次度量误差。上图表示。假设pi,pj点边塌缩前的端点，px为塌缩后的点坐标假设p_i,p_j点边塌缩前的端点，p_

闪电彬彬

10706人浏览 · 2022-05-10 22:57:02

闪电彬彬 · 2022-05-10 22:57:02 发布

欢迎关注更多精彩
关注我，学习常用算法与数据结构，一题多解，降维打击。

qem 全称Quadic Error Metrics 网格简化。
它是一种基于二次度量误差的优化算法。

边塌缩算法

qem算法的基本思想是对某一边进行塌缩，将边的2点合成一个点，如下图。

那么要选择哪条边进行塌缩呢。
就是基于二次度量误差。

上图表示。
$假设p_i,p_j点边塌缩前的端点，p_x为塌缩后的点坐标$

$那么p_x需要到p_i,p_i 各邻面的距离平方和最小$
能量如下：
$E_x = \displaystyle \sum_{P_s \in \Omega(i)}d(p_x,P_s)^2 + \displaystyle \sum_{P_s \in \Omega(j)}d(p_x,P_s)^2$

$d(p_x, P_s) 代表p_x到平面P_s的垂线距离$

$只要向量p_x在P_s法向上的投影减去平面上一点在法向上的投影即可$

$由于p_i是P_s上一点$
$d(p_x, P_s)=p_x\cdot n_s-p_i\cdot n_s$

$E_x = \displaystyle \sum_{P_s \in \Omega(i)}(p_x\cdot n_s-p_i\cdot n_s)^2 + \displaystyle \sum_{P_s \in \Omega(j)}(p_x\cdot n_s-p_j\cdot n_s)^2$

优化上述方程使得E最小就可以得到px坐标。

从矩阵角度来求距离平方和

$上式中n_i, -d_i都是常量，所以Q_i也是常$

再看上式

$令\overline{p_x} = (p_x, 1)$

$\displaystyle \sum_{P_s \in \Omega(i)}d(p_x,P_s)^2$

$=\displaystyle \sum_{P_s \in \Omega(i)}\overline{p_x}^TQ_s \overline{p_x}$

$=\overline{p_x}^T\left[\displaystyle \sum_{P_s \in \Omega(i)}Q_s \right]\overline{p_x}$

从上述可以看出可以先对某点的总和计算。

$Qt_i = \left[\displaystyle \sum_{P_s \in \Omega(i)}Q_s \right]$

$E_x=\overline{p_x}^TQt_i\overline{p_x}+\overline{p_x}^TQt_j\overline{p_x}$

$=\overline{p_x}^T(Qt_i+Qt_j)\overline{p_x}$

$可以看出上只有p_x是变量，是一个二次能量，可以通过优化得到p_x最优解，即为塌缩后的点坐标位置$

$将v看成是一个向量v=[x,y,z,1]，并展开得到下式，然后进行求导，并可以得到Av=b形式, 角得v=A^{-1}b$

算法步骤


for 所有点 i {
	计算点能量总和Qti
}

for 所有边 x{
	得到Ex, 并解出最优点px
}

while(点总数>目标值) {
	先出Ex最小的边
	将该边的两点从网格中删除。
	将两点从网格中删除。
	添加新的点。
	连接受影响的点。
	重新计算新边的能量。
}

代码实现

代码库

//
// Created by chenbinbin on 2022/5/9.
//

#include "hw9.h"
#include "../../PolyMesh/IOManager.h"
#include <Eigen/Dense>
#include <queue>
#include <map>
#include <time.h>
#include <unistd.h>
#include <string>

using namespace std;

#define lambda 1e3

using namespace acamcad;
using namespace polymesh;

struct Edge_priority
{
    double cost;
    MEdge* eh;
    int state;
    MPoint3 NewPoint;
};

struct cmp
{
    bool operator()(const Edge_priority &a, const Edge_priority &b)
    {
        return a.cost > b.cost;
    }
};

std::priority_queue<Edge_priority, std::vector<Edge_priority>, cmp> Cost;
std::map<MVert*, Eigen::Matrix4d> Qv;
std::map<MEdge*, int> State;


void cal_Q(MVert* vh, PolyMesh* mesh)
{
    MVector3 p_vh = vh->position();
    Eigen::Matrix4d Q_temp;
    Q_temp.setZero();
    MVector3 face_nor;
    double a, b, c, d;
    Eigen::Matrix<double, 4, 1> p;
    /*for (auto vf_it = mesh->vf_iter(vh); vf_it.isValid(); ++vf_it)
    {
        face_nor = (*vf_it)->normal();
        a = face_nor[0], b = face_nor[1], c = face_nor[2], d = -dot(face_nor, p_vh);
        p = { a,b,c,d };
        Q_temp += p * p.transpose();
    }
    Qv[vh] = Q_temp;*/
    if (!mesh->isBoundary(vh))
    {
        for (auto vf_it = mesh->vf_iter(vh); vf_it.isValid(); ++vf_it)
        {
            face_nor = (*vf_it)->normal();
            a = face_nor[0], b = face_nor[1], c = face_nor[2], d = -dot(face_nor, p_vh);
            p = { a,b,c,d };
            Q_temp += p * p.transpose();
        }
        Qv[vh] = Q_temp;
    }
    else
    {
        /*for (auto vf_it = mesh->vf_iter(vh); vf_it.isValid(); ++vf_it)
        {
            face_nor = (*vf_it)->normal();
            a = face_nor[0], b = face_nor[1], c = face_nor[2], d = -dot(face_nor, p_vh);
            p = { a,b,c,d };
            Q_temp += p * p.transpose();
            for (auto fhh_it = mesh->fhe_iter(*vf_it); fhh_it.isValid(); ++fhh_it)
            {
                if ((*fhh_it)->pair()->isBoundary())
                {
                    MVector3 vir_face_nor = (cross((*fhh_it)->pair()->tangent(), face_nor)).normalized();
                    a = vir_face_nor[0], b = vir_face_nor[1], c = vir_face_nor[2], d = -dot(vir_face_nor, p_vh);
                    p = { a,b,c,d };
                    Q_temp += lambda * (p * p.transpose());
                    break;
                }
            }
        }*/

        for (auto vhh_it = mesh->voh_iter(vh); vhh_it.isValid(); ++vhh_it)
        {
            if (!(*vhh_it)->isBoundary())
            {
                MPolyFace* vf = (*vhh_it)->polygon();
                MVector3 face_nor = vf->normal();
                double a = face_nor[0], b = face_nor[1], c = face_nor[2], d = -dot(face_nor, p_vh);
                Eigen::Matrix<double, 4, 1> p = { a,b,c,d };
                Q_temp += p * p.transpose();
            }
            else
            {
                MHalfedge* hh_boundary = *vhh_it;
                MHalfedge* prev_hh_boundary = hh_boundary->prev();
                MPolyFace* face_hh = hh_boundary->pair()->polygon(), *face_prev_hh = prev_hh_boundary->pair()->polygon();
                MVector3 face_hh_nor = face_hh->normal(), face_prev_hh_nor = face_prev_hh->normal();
                MVector3 vir_face_hh_nor = cross(hh_boundary->tangent(), face_hh_nor).normalized(), vir_face_prev_hh_nor = cross(prev_hh_boundary->tangent(), face_prev_hh_nor).normalized();
                double a = vir_face_hh_nor[0], b = vir_face_hh_nor[1], c = vir_face_hh_nor[2], d = -dot(vir_face_hh_nor, p_vh);
                Eigen::Matrix<double, 4, 1> p = { a,b,c,d };
                Q_temp += lambda * (p * p.transpose());
                a = vir_face_prev_hh_nor[0], b = vir_face_prev_hh_nor[1], c = vir_face_prev_hh_nor[2], d = -dot(vir_face_prev_hh_nor, p_vh);
                p = { a,b,c,d };
                Q_temp += lambda * (p * p.transpose());
            }
        }
        Qv[vh] = Q_temp;
    }
}

void cal_Cost(MEdge* eh)
{
    Edge_priority temp;
    MHalfedge* hh = eh->halfEdge();
    MVert* v_from = hh->fromVertex(), *v_to = hh->toVertex();
    Eigen::Matrix4d Q_plus = Qv[v_from] + Qv[v_to], Q_solve = Q_plus;
    Q_solve(3, 0) = 0.0, Q_solve(3, 1) = 0.0, Q_solve(3, 2) = 0.0, Q_solve(3, 3) = 1.0;
    MPoint3 new_point;
    Eigen::Vector4d new_vec;
    if (Q_solve.determinant() == 0)
    {
        MVector3 temp = 0.5*(v_from->position() + v_to->position());
        new_point = { temp[0], temp[1], temp[2] };
        new_vec = { new_point[0], new_point[1], new_point[2], 1.0 };
    }
    else
    {
        Eigen::Vector4d temp = { 0.0,0.0,0.0,1.0 };
        new_vec = Q_solve.inverse()*temp;
        new_point = { new_vec[0], new_vec[1], new_vec[2] };
    }
    temp.cost = new_vec.transpose()*Q_plus*new_vec;
    temp.eh = eh;
    temp.state = State[eh];
    temp.NewPoint = new_point;
    Cost.push(temp);
}

void update_Q_and_Cost(MVert* vh, PolyMesh* mesh)
{
    cal_Q(vh, mesh);
    for (auto vv_it = mesh->vv_iter(vh); vv_it.isValid(); ++vv_it)
    {
        cal_Q(*vv_it, mesh);
    }
    for (auto ve_it = mesh->ve_iter(vh); ve_it.isValid(); ++ve_it)
    {
        State[*ve_it]++;
        cal_Cost(*ve_it);
    }
}

bool QEM_collapse(const Edge_priority &temp_edge, PolyMesh* mesh)
{
    bool is_collapse = false;
    MEdge* eh = temp_edge.eh;
    MHalfedge* hh = eh->halfEdge(), *hh_oppo = hh->pair();
    MVert *v_from = hh->fromVertex(), *v_to = hh->toVertex();
    MVert* vh;
//    if (mesh->is_collapse_ok(hh))
	if (mesh->is_collapse_ok_Triangle(hh))
    {
        v_to->setPosition(temp_edge.NewPoint);
        //mesh->collapse(hh);
        mesh->collapseTriangle(hh);
        is_collapse = true;
        vh = v_to;
    }
//    else if (mesh->is_collapse_ok(hh_oppo))
	else if (mesh->is_collapse_ok_Triangle(hh_oppo))
    {
        v_from->setPosition(temp_edge.NewPoint);
        //mesh->collapse(hh_oppo);
        mesh->collapseTriangle(hh_oppo);
        is_collapse = true;
        vh = v_from;
    }
    if (is_collapse)
    {
        update_Q_and_Cost(vh, mesh);
    }
    return is_collapse;
}

void QEM(PolyMesh* mesh)
{
    if (mesh->numVertices() == 3)
        return;
    /// Initial Qv
    for (MVert* vh : mesh->vertices())
    {
        MVector3 p_vh = vh->position();
        Eigen::Matrix4d Q_temp;
        Q_temp.setZero();
        MVector3 face_nor;
        double a, b, c, d;
        Eigen::Matrix<double, 4, 1> p;
        if (!mesh->isBoundary(vh))
        {
            for (auto vf_it = mesh->vf_iter(vh); vf_it.isValid(); ++vf_it)
            {
                face_nor = (*vf_it)->normal();
                a = face_nor[0], b = face_nor[1], c = face_nor[2], d = -dot(face_nor, p_vh);
                p = { a,b,c,d };
                Q_temp += p * p.transpose();
            }
            Qv.insert(std::make_pair(vh, Q_temp));
        }
        else
        {
            for (auto vhh_it = mesh->voh_iter(vh); vhh_it.isValid(); ++vhh_it)
            {
                if (!(*vhh_it)->isBoundary())
                {
                    MPolyFace* vf = (*vhh_it)->polygon();
                    MVector3 face_nor = vf->normal();
                    double a = face_nor[0], b = face_nor[1], c = face_nor[2], d = -dot(face_nor, p_vh);
                    Eigen::Matrix<double, 4, 1> p = { a,b,c,d };
                    Q_temp += p * p.transpose();
                }
                else
                {
                    MHalfedge* hh_boundary = *vhh_it;
                    MHalfedge* prev_hh_boundary = hh_boundary->prev();
                    MPolyFace* face_hh = hh_boundary->pair()->polygon(), *face_prev_hh = prev_hh_boundary->pair()->polygon();
                    MVector3 face_hh_nor = face_hh->normal(), face_prev_hh_nor = face_prev_hh->normal();
                    MVector3 vir_face_hh_nor = cross(hh_boundary->tangent(), face_hh_nor).normalized(), vir_face_prev_hh_nor = cross(prev_hh_boundary->tangent(), face_prev_hh_nor).normalized();
                    double a = vir_face_hh_nor[0], b = vir_face_hh_nor[1], c = vir_face_hh_nor[2], d = -dot(vir_face_hh_nor, p_vh);
                    Eigen::Matrix<double, 4, 1> p = { a,b,c,d };
                    Q_temp += lambda * (p * p.transpose());
                    a = vir_face_prev_hh_nor[0], b = vir_face_prev_hh_nor[1], c = vir_face_prev_hh_nor[2], d = -dot(vir_face_prev_hh_nor, p_vh);
                    p = { a,b,c,d };
                    Q_temp += lambda * (p * p.transpose());
                }
            }
            Qv.insert(std::make_pair(vh, Q_temp));
        }

        /*for (auto vf_it = mesh->vf_iter(vh); vf_it.isValid(); ++vf_it)
        {
            face_nor = (*vf_it)->normal();
            double a = face_nor[0], b = face_nor[1], c = face_nor[2], d = -dot(face_nor, p_vh);
            Eigen::Matrix<double, 4, 1> p = { a,b,c,d };
            Q_temp += p * p.transpose();
        }
        Qv.insert(std::make_pair(vh, Q_temp));*/
    }

    ///Initial Cost
    for (MEdge* eh : mesh->edges())
    {
        State.insert(std::make_pair(eh, 0));
        cal_Cost(eh);
    }

    /// simplification
    int N_V = mesh->numVertices();
    int target_num = std::min((int)(0.5*N_V), 1000);
    while(N_V > target_num)
    {
        Edge_priority temp_edge = Cost.top();
        Cost.pop();
        MEdge* eh = temp_edge.eh;
        if (temp_edge.state == State[eh])
        {
            if (eh->index() != -1)
            {
                if (QEM_collapse(temp_edge, mesh))
                {
                    N_V--;
                }
            }
        }
    }
}

void qem_simplification()
{
    char buffer[500];
    getcwd(buffer, 500);
    printf("The current directory is: %s/../\n", buffer);
    string mesh_path = buffer;
    mesh_path += "/../src/hw9/dragon.obj";
    PolyMesh* mesh = new PolyMesh();
    loadMesh(mesh_path, mesh);

    // read input mesh
    //PolyMesh* mesh = new PolyMesh();
    //loadMesh("cat_open.obj", mesh);

    clock_t start, end;
    std::cout << "Begin QEM" << std::endl;
    start = clock();
    QEM(mesh);
    end = clock();
    std::cout << "time: " << (double)(end - start) / CLOCKS_PER_SEC << "s" << std::endl;
    writeMesh("dragon-simplification-tri.obj", mesh);
}