几何建模与处理-网格简化的QEM方法老铁，点个在看呗欢迎关注公众号:sumsmile /专注图形渲染、Android开

老铁，点个在看呗 欢迎关注公众号:sumsmile /专注图形渲染、Android开发

基于Carnegie Mellon University的一篇paperQuadric Error Metrics，实现模型简化。准确的说，应该是离散曲面简化。代码参考马玺凯的实现。

实现效果

图中原模型有43243个顶点，86482个面，每次简化减少一半顶点，12次简化后剩22个顶点40个面，当然，只要你愿意，可以一直简化，只是简化的太过分就失去了原模型的特征。

实现原理

原理不复杂，paper讲的非常细致，只要耐心的读读就能理解。

每次减少一条边，即合并两个点v1 v2，如下图所示。

问题的关键有两点：

应该选择哪条边来收缩？
合并后的点的坐标怎么取？

如何选择边

直观理解，应该选择最不重要的边，即减少这条边对模型外观影响最小。

对照图2，假设收缩边 $\overrightarrow{v_1v_2}$ ，即合并顶点v1、v2，用一个新的点 $\overline v$ 替代。

与v1相连有n1个三角面，v_1到这n1个面的距离之和为 $d_1$

计 $d_1$ 为合并v1的影响，同理合并 $v_2$ 的影响为 $d_2$ ，则合并v1、v2对模型整体的影响为 $d_1 + d_2$

单个点到邻接面的距离平方和，用严谨的数学公式表达为：

注意，这里P为平面方程的系数，是列向量，用平方来度量去掉了负号，方便计算。

先暂停下，复习点到平面的距离求解

三维平面可表达为： ax + by + cz + d = 0

为方便计算，要求 a b c是归一化后的系数，即 $a^2 + b^2 + c^2 = 1$

向量 $\overrightarrow n(a, b, c)$ 为该平面的法向量，d为该平面相对原点的位移。则 $\overrightarrow {ov}$ 投影到 $\overrightarrow n$ 上的长度为： $\overrightarrow{n} * \overrightarrow{v}$

还要考虑平面相对原点的平移d。v用齐次坐标表示，则可得v到平面的距离为：

\Delta v = P^T * v = \left[ \begin{matrix} a\\ b\\ c\\ d\\ \end{matrix} \right] * \left[ \begin{matrix} x\\ y\\ z\\ 1\\ \end{matrix} \right]

这是从投影的角度直观的理解点到平面的距离，不是很严谨

回到正题

分别计算出每对点合并后对模型的影响(点到邻接面的平方和)，排序后可得到影响最小的那对点。这个影响用 $\Delta v$ 表示

计算合并后的点的坐标

上面计算 $\Delta v$ 的公式中，v为合并后的新点。继续看 $\Delta v$ 的计算：

其中，

P= \left[ \begin{matrix} a\\ b\\ c\\ d\\ \end{matrix} \right]

则

这个矩阵是这篇paper最核心的内容，每个顶点都可求得这样一个矩阵，确定了合并后的点，乘以这个矩阵就得到一个"影响值"

其中， $Q = \sum_{i = 0} ^m K_p$ ，即点邻接面的度量计算参数之和

看上图， $\Delta v$ 是二次函数，是个凹函数，对其求导，导数为0处对应着最小值 $\overline v$ ，即：

这里可能比较跳跃，自己动手算下，对图6求偏导

如果该矩阵可逆，则该方程有解，可得到 $\overline v$

如果不可逆，求两点平均值： $\overline v = \frac{1}{2} * (v_1 + v_2)$

注意：收缩一条边，实际减少两个点，所以上面Q是收缩前两个点的度量 $\overline Q$ $\overline Q = \overline Q_1 + \overline Q_2$

注意，这么算，会重复考虑了相同的面，不过不要紧，这只是个近似的计算，从结果看影响不大。

代码说明

开发环境：C++ CLion, 代码参考马玺凯同学的实现
依赖库：libigl,用的不是最新的代码(checkout commit 3cb4894e)

代码：
github.com/summer-go/G…

libigl简单而强大，可以处理各种Mesh操作，而且代码文档非常完善，对数学原理也一并做了简要介绍，非常值得学习。 libigl tutorial

主函数逻辑


int main(int argc, char *argv[])
{
...

    if(!OpenMesh::IO::read_mesh(mesh,"../models/armadillo.obj",opt))
    {
        std::cout << "can't open file"<<std::endl;
    }
    for(MyMesh::VertexIter v_it = mesh.vertices_begin(); v_it != mesh.vertices_end(); ++v_it)
    {
        mesh.point(*v_it) /= 4.0;
    }
    mesh.update_face_normals();
    mesh.update_vertex_normals();
    igl::opengl::glfw::Viewer viewer;

    // 注册按键回调，用来触发简化、还原操作
    viewer.callback_key_down = &key_down;

    // 转成igl需要的数据，点和面用矩阵来存储
    // V F数组长度为12，即缓存了12组简化结果，方便回溯查看
    openMesh_to_igl(mesh,V[0],F[0]);

    std::cout << "vertices : " <<  mesh.n_vertices() << std::endl;
    std::cout << "faces : " << mesh.n_faces() << std::endl;

    // 模型数据设置到igl的viewer中
    viewer.data().set_mesh(V[0], F[0]);
    // 默认展示线框
    viewer.data().show_lines = true;

    // 打开window
    viewer.launch();

}

按键回调操作

bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int modifier)
{
    // 按1 简化模型
    if(key == '1')
    {
        if(ti < cache_size)
        {
            if(!flag[ti])
            {
                // 生成一个简化的模型，0.5表示简化操作保留一半顶点
                Surface_Simplification(mesh,0.5);
                std::cout << "vertices : " << mesh.n_vertices() << std::endl;;
                std::cout << "faces : " << mesh.n_faces() << std::endl;;
                // 转换成igl接受的点面数据
                openMesh_to_igl(mesh,V[ti],F[ti]);
                flag[ti] = true;
            }
            // clear view里的数据
            viewer.data().clear();
            // 更新模型数据
            viewer.data().set_mesh(V[ti],F[ti]);
            // 模型居中
            viewer.core().align_camera_center(V[ti],F[ti]);
            // ti 记录当前简化的次数
            ti++;
        }

    }
    
    // 按2，还原模型，从缓存里读取对应的模型数据
    else if(key == '2')
    {
        if(ti > 1)
        {
            ti--;
            viewer.data().clear();
            viewer.data().set_mesh(V[ti],F[ti]);
            viewer.core().align_camera_center(V[ti],F[ti]);
        }
    }
    return false;
}

模型简化计算

void Surface_Simplification(MyMesh &mesh, float ratio)
{
    assert(ratio >=0 && ratio <= 1);
    int it_num = (1.0f - ratio) * mesh.n_vertices();

    //1. Compute the Q matrices for all the initial vertices
    auto Q = OpenMesh::makeTemporaryProperty<OpenMesh::VertexHandle, Eigen::Matrix4f>(mesh);
    auto v = OpenMesh::makeTemporaryProperty<OpenMesh::VertexHandle, Eigen::Vector4f>(mesh);
    auto flag = OpenMesh::makeTemporaryProperty<OpenMesh::VertexHandle, int>(mesh); //与vto_flag vfrom_flag配合判断点对是否还有效
    auto p = OpenMesh::makeTemporaryProperty<OpenMesh::FaceHandle, Eigen::Vector4f>(mesh);

    // 求每个面的表达式，遍历面，取出法线作为平面方程系数
    // 平面：ax + by + cz + d = 0, [a, b, c]为该平面的法线， d 为平面从原点(0, 0, 0) 移动的距离
    for(MyMesh::FaceIter fit = mesh.faces_begin(); fit != mesh.faces_end(); ++fit)
    {
        float a, b, c, d;

        // 获取法线，得到参数 a b c
        a = mesh.normal(*fit)[0];
        b = mesh.normal(*fit)[1];
        c = mesh.normal(*fit)[2];

        // 代入平面上任意一点，求得 d
        MyMesh::Point tp = mesh.point(fit->halfedge().to());
        d = - (a * tp[0] + b * tp[1] + c * tp[2]);
        p[*fit][0] = a;
        p[*fit][1] = b;
        p[*fit][2] = c;
        p[*fit][3] = d;
    }

    // 求每个点的Q.遍历点，求出每个点到面的"距离矩阵"
    for(MyMesh::VertexIter vit = mesh.vertices_begin(); vit != mesh.vertices_end(); ++vit)
    {
        Eigen::Matrix4f mat;
        mat.setZero();
        for(MyMesh::VertexFaceIter vf_it = mesh.vf_iter(*vit); vf_it.is_valid(); ++vf_it)
        {
            mat += p[*vf_it] * p[*vf_it].transpose();
        }
        Q[*vit] = mat;
        v[*vit][0] = mesh.point(*vit)[0];
        v[*vit][1] = mesh.point(*vit)[1];
        v[*vit][2] = mesh.point(*vit)[2];
        v[*vit][3] = 1.0f;
        flag[*vit] = 0;
    }

    // 使用优先级队列，实现一个小顶堆，从小到大排序
    // 每次都用新的q
    std::priority_queue <edge_Collapse_structure, std::vector<edge_Collapse_structure>, std::less<edge_Collapse_structure>> q;
    for(MyMesh::EdgeIter eit = mesh.edges_begin(); eit != mesh.edges_end(); ++eit)
    {
        Eigen::Matrix4f newQ = Q[eit->v0()] + Q[eit->v1()];
        Eigen::Matrix4f tQ = newQ;
        Eigen::Vector4f b(0.0f,0.0f,0.0f,1.0f);

        // 对Q求偏导得到的矩阵和Q区别在于：偏导后最后一行为[0, 0, 0, 1]
        tQ(3,0) = 0.0f;
        tQ(3,1) = 0.0f;
        tQ(3,2) = 0.0f;
        tQ(3,3) = 1.0f;
        Eigen::FullPivLU<Eigen::Matrix4f> lu(tQ);
        Eigen::Vector4f vnew;
        // if is invertible, solve the linear equation
        if(lu.isInvertible())
        {
            vnew = tQ.inverse() * b;
        }
        // else take the midpoint
        else
        {
            vnew = (v[eit->v0()] + v[eit->v1()]) / 2.0f;
        }
        //std::cout << vnew << std::endl;
        edge_Collapse_structure ts;
        ts.hf = eit->halfedge(0);
        ts.cost = vnew.transpose() * newQ * vnew;
        MyMesh::Point np(vnew[0], vnew[1], vnew[2]);
        ts.np = np;
        ts.vto = eit->halfedge(0).to();
        ts.vfrom = eit->halfedge(0).from();
        ts.Q_new = newQ;
        //
        q.push(ts);
    }
    mesh.request_vertex_status();
    mesh.request_edge_status();
    mesh.request_face_status();

    int i = 0;
    
    // 上面都是准备工作，
    // 真正进入到边收缩阶段，迭代it_num次
    while( i < it_num)
    {
        edge_Collapse_structure s = q.top();
        q.pop();
        if(mesh.status(s.vfrom).deleted() || mesh.status(s.vto).deleted())
            continue;
        if(s.vto_flag != flag[s.vto] || s.vfrom_flag != flag[s.vfrom])
            continue;
        MyMesh::VertexHandle tvh;
        if(mesh.is_collapse_ok(s.hf))
        {
            mesh.collapse(s.hf);
            tvh = s.vto;
            flag[s.vto] ++;
            flag[s.vfrom] ++;
        }

        else if(mesh.is_collapse_ok(mesh.opposite_halfedge_handle(s.hf)))
        {
            mesh.collapse(mesh.opposite_halfedge_handle(s.hf));
            tvh = s.vfrom;
            flag[s.vto] ++;
            flag[s.vfrom] ++;
        }
        else
        {
            continue;
        }

        mesh.set_point(tvh,s.np);
        Q[tvh] = s.Q_new;
        v[tvh][0] = s.np[0];
        v[tvh][1] = s.np[1];
        v[tvh][2] = s.np[2];
        v[tvh][3] = 1.0f;
        for(MyMesh::VertexOHalfedgeIter vh_it = mesh.voh_iter(tvh); vh_it.is_valid(); ++vh_it)
        {
            MyMesh::VertexHandle tt = vh_it->to();
            Eigen::Matrix4f newQ = s.Q_new + Q[tt];
            Eigen::Matrix4f tQ = newQ;
            Eigen::Vector4f b(0.0f,0.0f,0.0f,1.0f);
            tQ(3,0) = 0.0f;
            tQ(3,1) = 0.0f;
            tQ(3,2) = 0.0f;
            tQ(3,3) = 1.0f;
            Eigen::FullPivLU<Eigen::Matrix4f> lu(tQ);
            Eigen::Vector4f vnew;
            // if is invertible, solve the linear equation
            if(lu.isInvertible())
            {
                vnew = tQ.inverse() * b;
            }
                // else take the midpoint
            else
            {
                vnew = (v[tvh] + v[tt]) / 2.0f;
            }
            //std::cout << vnew << std::endl;
            edge_Collapse_structure ts;
            ts.hf = *vh_it;
            ts.cost = vnew.transpose() * newQ * vnew;
            MyMesh::Point np(vnew[0], vnew[1], vnew[2]);
            ts.np = np;
            ts.vto = tt;
            ts.vto_flag = flag[tt];
            ts.vfrom = tvh;
            ts.vfrom_flag = flag[tvh];
            ts.Q_new = newQ;
            q.push(ts);
        }
        i++;

    }
    mesh.garbage_collection();
}

老铁，点个在看呗 欢迎关注公众号:sumsmile /专注图形渲染、Android开发

参考资料

[1] Quadric Error Metrics: mgarland.org/files/paper…

[2] libigl tutorial: libigl.github.io/tutorial/