简介

即三角网格插值细分

参考文献

http://graphics.stanford.edu/courses/cs468-10-fall/LectureSlides/10_Subdivision.pdf

TIPS

w 采用1/16.

参考示意图

image

code

bug 调了很久，就是 一条半边对应面的方向其实不太稳定。

codeMain

#pragma once
#pragma once
#include "strategy.h"
#include "qwidget.h"
#include <OpenMesh/Core/Mesh/PolyMesh_ArrayKernelT.hh>
#include <OpenMesh/Core/IO/MeshIO.hh>
#include <map>



class Butterfly :public Strategy, public QWidget
{
    Q_OBJECT
private:
    typedef OpenMesh::PolyMesh_ArrayKernelT<> MyMesh;
    MyMesh mesh;
    std::map<int, OpenMesh::Vec3d> facePoints;
    std::map<int, OpenMesh::Vec3d> edgePoints;
    std::map<int, OpenMesh::Vec3d> vertexPoints;
    int times;
    const double w = 1.0 / 16.0;
public:
    Butterfly(Data* data_);
    ~Butterfly();
    void genCube();
    void genFacePoint();
    void genEdgePoint();
    void genVertexPoint();
    void connectPoints();
    void genMesh(std::string name);
    bool Run();
    void getResult();
};

#include "butterfly_surface_subdivision.h"
#include <iostream>
#include <unordered_map>
#include <QInputDialog>

/**
 * @description: 构造函数
 * @param {type}
 * @return {type}
 */
Butterfly::Butterfly(Data* data_) : Strategy(data_) {
    times = 1;
    genCube();
}

/**
 * @description: 析构函数
 * @param {type}
 * @return {type}
 */
Butterfly::~Butterfly() {}

/**
 * @description: 一整套流程
 * @param {type}
 * @return {int} 管线运行是否成功
 */
bool Butterfly::Run() {
    times = QInputDialog::getInt(this, "Surface Mesh", "Please input times",
        1, 1, 1000, 1);
    for (int i = 0; i <= times; i++) {
        if (i == 0) {
            char a[20];
            sprintf(a, "output%d.off", i);
            genMesh(a);
        }
        else {
            genEdgePoint();
            genVertexPoint();
            connectPoints();
            char a[20];
            sprintf(a, "output%d.off", i);
            genMesh(a);
        }
        // 清空变量
        std::cout << "[DEBUG] 迭代了第 " << i << " 次。" << std::endl;
    }
    getResult();

    return true;
}


/**
 * @description: 生成一个立方体（四边形网格）
 * @param {type}
 * @return {type}
 */
void Butterfly::genCube()
{
    if (getData()->triMesh.n_vertices() != 0) {
        for (auto it : getData()->triMesh.vertices()) {
            mesh.add_vertex(getData()->triMesh.point(it));
        }
        for (auto it : getData()->triMesh.faces()) {
            std::vector<MyMesh::VertexHandle> face_vhandles;
            for (auto iit : it.vertices()) {
                face_vhandles.push_back(iit);
            }
            mesh.add_face(face_vhandles);
        }
        return;
    }
    MyMesh::VertexHandle vhandle[9];
    vhandle[0] = mesh.add_vertex(MyMesh::Point(-1, -1, 1));
    vhandle[1] = mesh.add_vertex(MyMesh::Point(1, -1, 1));
    vhandle[2] = mesh.add_vertex(MyMesh::Point(1, 1, 1));
    vhandle[3] = mesh.add_vertex(MyMesh::Point(-1, 1, 1));
    vhandle[4] = mesh.add_vertex(MyMesh::Point(-1, -1, -1));
    vhandle[5] = mesh.add_vertex(MyMesh::Point(1, -1, -1));
    vhandle[6] = mesh.add_vertex(MyMesh::Point(1, 1, -1));
    vhandle[7] = mesh.add_vertex(MyMesh::Point(-1, 1, -1));

    std::vector<MyMesh::VertexHandle> face_vhandles;
    face_vhandles.push_back(vhandle[0]);
    face_vhandles.push_back(vhandle[1]);
    face_vhandles.push_back(vhandle[2]);
    mesh.add_face(face_vhandles);
    face_vhandles.clear();
    face_vhandles.push_back(vhandle[0]);
    face_vhandles.push_back(vhandle[2]);
    face_vhandles.push_back(vhandle[3]);
    mesh.add_face(face_vhandles);
    face_vhandles.clear();

    face_vhandles.push_back(vhandle[7]);
    face_vhandles.push_back(vhandle[6]);
    face_vhandles.push_back(vhandle[5]);
    mesh.add_face(face_vhandles);
    face_vhandles.clear();
    face_vhandles.push_back(vhandle[7]);
    face_vhandles.push_back(vhandle[5]);
    face_vhandles.push_back(vhandle[4]);
    mesh.add_face(face_vhandles);
    face_vhandles.clear();

    face_vhandles.push_back(vhandle[1]);
    face_vhandles.push_back(vhandle[0]);
    face_vhandles.push_back(vhandle[4]);
    mesh.add_face(face_vhandles);
    face_vhandles.clear();
    face_vhandles.push_back(vhandle[1]);
    face_vhandles.push_back(vhandle[4]);
    face_vhandles.push_back(vhandle[5]);
    mesh.add_face(face_vhandles);
    face_vhandles.clear();

    face_vhandles.push_back(vhandle[2]);
    face_vhandles.push_back(vhandle[1]);
    face_vhandles.push_back(vhandle[5]);
    mesh.add_face(face_vhandles);
    face_vhandles.clear();
    face_vhandles.push_back(vhandle[2]);
    face_vhandles.push_back(vhandle[5]);
    face_vhandles.push_back(vhandle[6]);
    mesh.add_face(face_vhandles);
    face_vhandles.clear();

    face_vhandles.push_back(vhandle[3]);
    face_vhandles.push_back(vhandle[2]);
    face_vhandles.push_back(vhandle[6]);
    mesh.add_face(face_vhandles);
    face_vhandles.clear();
    face_vhandles.push_back(vhandle[3]);
    face_vhandles.push_back(vhandle[6]);
    face_vhandles.push_back(vhandle[7]);
    mesh.add_face(face_vhandles);
    face_vhandles.clear();

    face_vhandles.push_back(vhandle[0]);
    face_vhandles.push_back(vhandle[3]);
    face_vhandles.push_back(vhandle[7]);
    mesh.add_face(face_vhandles);
    face_vhandles.clear();
    face_vhandles.push_back(vhandle[0]);
    face_vhandles.push_back(vhandle[7]);
    face_vhandles.push_back(vhandle[4]);
    mesh.add_face(face_vhandles);

}


/**
 * @description: 考虑边上和内部
 * @param {type}
 * @return {type}
 */
void Butterfly::genEdgePoint()
{
    edgePoints.clear();
    for (auto e_it = mesh.edges_begin(); e_it != mesh.edges_end(); ++e_it)
    {
        // 得到边所代表的半边
        OpenMesh::HalfedgeHandle heh1 = mesh.halfedge_handle(*e_it, 0); // 默认一个方向的半边
        OpenMesh::Vec3d edgePoint(0, 0, 0);
        int edgePointsNumber = 0;
        OpenMesh::DefaultTraits::Point pointV = mesh.point(mesh.from_vertex_handle(heh1)); // 这条（半）边的起点
        OpenMesh::DefaultTraits::Point pointW = mesh.point(mesh.to_vertex_handle(heh1));   // 这条（半）边的终点
        

        OpenMesh::FaceHandle fh1 = mesh.face_handle(heh1);
        OpenMesh::HalfedgeHandle heh = mesh.opposite_halfedge_handle(heh1);
        OpenMesh::FaceHandle fh2 = mesh.face_handle(heh);
        std::vector<OpenMesh::VertexHandle> ss1; // d3 d4
        std::vector<OpenMesh::VertexHandle> ss2; // d5 d6 d7 d8
        std::vector<OpenMesh::VertexHandle> ss_all;
        // 得到这两个面的所有顶点
        for (auto fv_it : mesh.fv_range(fh1)) {
            if (fv_it.idx() != mesh.from_vertex_handle(heh1).idx() &&
                fv_it.idx() != mesh.to_vertex_handle(heh1).idx())
                ss1.push_back(fv_it);
        }
        for (auto fv_it : mesh.fv_range(fh2)) {
            if (fv_it.idx() != mesh.from_vertex_handle(heh1).idx() &&
                fv_it.idx() != mesh.to_vertex_handle(heh1).idx())
                ss1.push_back(fv_it); 
        }
        if (ss1.size() != 2) {
            std::cout << "[ERROR] " << ss1.size() << std::endl;
        }
        ss_all.assign(ss1.begin(), ss1.end());
        ss_all.push_back(mesh.from_vertex_handle(heh1));
        ss_all.push_back(mesh.to_vertex_handle(heh1));
        for (const auto& fe : mesh.fe_range(fh1)) {
            OpenMesh::HalfedgeHandle fe_heh1 = mesh.halfedge_handle(fe, 0);
            //OpenMesh::FaceHandle fe_fh1 = mesh.face_handle(fe_heh1);
            OpenMesh::HalfedgeHandle fe_heh = mesh.opposite_halfedge_handle(fe_heh1);
            OpenMesh::FaceHandle fe_fh2 = mesh.face_handle(fe_heh);
            if (fe_fh2 == fh1 || fe_fh2 == fh2) {
                OpenMesh::FaceHandle fe_fh1 = mesh.face_handle(fe_heh1);
                if (fe_fh1 == fh1 || fe_fh1 == fh2) continue;
                for (const auto& ffv : mesh.fv_range(fe_fh1)) {
                    if (std::find(ss_all.begin(), ss_all.end(), ffv) == ss_all.end()) {
                        ss2.push_back(ffv);
                    }
                }
            }
            else {
                for (const auto& ffv : mesh.fv_range(fe_fh2)) {
                    if (std::find(ss_all.begin(), ss_all.end(), ffv) == ss_all.end()) {
                        ss2.push_back(ffv);
                    }
                }
            }
        }
        //std::cout << "=1=\n";
        //for (const auto& ffv : mesh.fv_range(fh1)) {
        //    std::cout << mesh.point(ffv)[0] << " " << mesh.point(ffv)[1] << " " << mesh.point(ffv)[2] << "\n";
        //}
        //for (const auto& ffv : mesh.fv_range(fh2)) {
        //    std::cout << mesh.point(ffv)[0] << " " << mesh.point(ffv)[1] << " " << mesh.point(ffv)[2] << "\n";
        //}


        //std::cout << mesh.point(ss_all[0])[0] << " " << mesh.point(ss_all[0])[1] << " " << mesh.point(ss_all[0])[2] << "\n";
        //std::cout << mesh.point(ss_all[1])[0] << " " << mesh.point(ss_all[1])[1] << " " << mesh.point(ss_all[1])[2] << "\n";
        //std::cout << mesh.point(ss_all[2])[0] << " " << mesh.point(ss_all[2])[1] << " " << mesh.point(ss_all[2])[2] << "\n";
        //std::cout << mesh.point(ss_all[3])[0] << " " << mesh.point(ss_all[3])[1] << " " << mesh.point(ss_all[3])[2] << "\n";
        //std::cout << "==\n";
        for (const auto& fe : mesh.fe_range(fh2)) {
            OpenMesh::HalfedgeHandle fe_heh1 = mesh.halfedge_handle(fe, 0);
            //OpenMesh::FaceHandle fe_fh1 = mesh.face_handle(fe_heh1);
            OpenMesh::HalfedgeHandle fe_heh = mesh.opposite_halfedge_handle(fe_heh1);
            OpenMesh::FaceHandle fe_fh2 = mesh.face_handle(fe_heh);
            if (fe_fh2 == fh1 || fe_fh2 == fh2) {
                OpenMesh::FaceHandle fe_fh1 = mesh.face_handle(fe_heh1);
                if (fe_fh1 == fh1 || fe_fh1 == fh2) continue;
                for (const auto& ffv : mesh.fv_range(fe_fh1)) {
                    if (std::find(ss_all.begin(), ss_all.end(), ffv) == ss_all.end()) {
                        ss2.push_back(ffv);
                    }
                }
            }
            else {
                for (const auto& ffv : mesh.fv_range(fe_fh2)) {
                    if (std::find(ss_all.begin(), ss_all.end(), ffv) == ss_all.end()) {
                        ss2.push_back(ffv);
                    }
                }
            }
        }

        if (ss2.size() != 4) {
            std::cout << "[ERROR] " << ss2.size() << std::endl;
        }

        edgePoints[heh1.idx()] = 1.0 / 2.0 * (pointV + pointW) + w * (mesh.point(ss1[0]) + mesh.point(ss1[1]))
                -w/2.0 *(mesh.point(ss2[0]) + mesh.point(ss2[1]) + mesh.point(ss2[2]) + mesh.point(ss2[3]));

    }
}

/**
 * @description: 生成新的顶点
 * @param {type}
 * @return {type}
 */
void Butterfly::genVertexPoint()
{
    vertexPoints.clear();
    // 原始点接触的面的所有的面点的均值
    for (auto v_it = mesh.vertices_begin(); v_it != mesh.vertices_end(); v_it++)
    {
        if (mesh.is_boundary(v_it)) { // 如果这个顶点处于边界
            OpenMesh::Vec3d point(0, 0, 0);
            for (auto vv_it = mesh.vv_begin(v_it); vv_it.is_valid(); vv_it++) {
                if (mesh.is_boundary(vv_it)) {
                    point += mesh.point(vv_it);
                }
            }
            vertexPoints[(*v_it).idx()] = 3.0 / 4.0 * mesh.point(v_it) + 1.0 / 8.0 * (point);
        }
        else {
            // 计算度 n, 度
            int n = 0;
            for (auto vf_it = mesh.vv_begin(v_it); vf_it.is_valid(); ++vf_it) {
                n++;
            }
            // 计算β
            double beta = 1.0 / n * (5.0 / 8.0 - pow(3.0 / 8.0 + 1.0 / 4.0 * cos(2 * M_PI / n), 2));
            // 下面的需要么??

            if (n < 3)
                beta = 3.0 / 16.0;


            OpenMesh::Vec3d point(0, 0, 0);
            for (auto vv_it = mesh.vv_begin(v_it); vv_it.is_valid(); ++vv_it) {
                point += beta * mesh.point(vv_it);
            }
            vertexPoints[(*v_it).idx()] = (1 - n * beta) * mesh.point(v_it) + point;
        }
    }
}

/**
 * @description: 连接面点和边点
 * @param {type}
 * @return {type}
 */
void Butterfly::connectPoints() {
    if (!mesh.has_vertex_status())  mesh.request_vertex_status();
    if (!mesh.has_face_status())    mesh.request_face_status();
    if (!mesh.has_edge_status())    mesh.request_edge_status();
    // 先将旧顶点移动到新顶点
    std::unordered_map<int, OpenMesh::VertexHandle> vertexHandle;
    for (auto v_it = mesh.vertices_begin(); v_it != mesh.vertices_end(); v_it++) {
        //mesh.set_point(*v_it, (OpenMesh::DefaultTraits::Point)vertexPoints[(*v_it).idx()]);
        vertexHandle[(*v_it).idx()] = *v_it;
    }
    std::vector<MyMesh::VertexHandle> facePointsHandle;
    std::vector<MyMesh::FaceHandle> faceHandle;
    std::unordered_map<int, MyMesh::VertexHandle> edgesHandle;
    // 再加入所有的要加入的边点
    for (const auto& eh : mesh.edges()) {
        edgesHandle[mesh.halfedge_handle(eh, 0).idx()] = mesh.add_vertex(MyMesh::Point(edgePoints[mesh.halfedge_handle(eh, 0).idx()]));
    }
    for (const auto& fh : mesh.faces()) {
        faceHandle.push_back(fh);
    }
    std::vector<std::vector<MyMesh::VertexHandle>> v;
    for (const auto& fh : mesh.faces()) {
        std::vector<MyMesh::VertexHandle> handlerVector(6);
        int pointsNumber = 0;
        for (const auto& fvh : mesh.fv_range(fh)) {
            handlerVector[pointsNumber] = vertexHandle[fvh.idx()];
            pointsNumber++;
        }
        for (const auto& feh : mesh.fe_range(fh)) {
            auto fehh = mesh.halfedge_handle(feh, 0);
            handlerVector[pointsNumber] = edgesHandle[fehh.idx()];
            pointsNumber++;
        }
        v.push_back(handlerVector);
    }

    for (int i = 0; i < v.size(); i++) {
        mesh.delete_face(faceHandle[i], true);
        std::vector<MyMesh::VertexHandle> handlerVector;
        for (int j = 0; j < v[i].size(); j++) {
            handlerVector.push_back(v[i][j]);
        }

        std::vector<MyMesh::VertexHandle> face_vhandles;
        face_vhandles.push_back(handlerVector[3]);
        face_vhandles.push_back(handlerVector[0]);

        face_vhandles.push_back(handlerVector[4]);
        mesh.add_face(face_vhandles);
        face_vhandles.clear();

        face_vhandles.push_back(handlerVector[4]);
        face_vhandles.push_back(handlerVector[1]);

        face_vhandles.push_back(handlerVector[5]);
        mesh.add_face(face_vhandles);
        face_vhandles.clear();

        face_vhandles.push_back(handlerVector[3]);
        face_vhandles.push_back(handlerVector[4]);
        face_vhandles.push_back(handlerVector[5]);
        mesh.add_face(face_vhandles);
        face_vhandles.clear();

        face_vhandles.push_back(handlerVector[2]);
        face_vhandles.push_back(handlerVector[3]);
        face_vhandles.push_back(handlerVector[5]);
        mesh.add_face(face_vhandles);
        face_vhandles.clear();
    }

    mesh.garbage_collection();
    if (mesh.has_vertex_status())  mesh.release_vertex_status();
    if (mesh.has_face_status())    mesh.release_face_status();
    if (mesh.has_edge_status())    mesh.release_edge_status();
}



/**
 * @description: 输出新网格
 * @param {string} name 文件名称 默认 output.obj
 * @return {type}
 */
void Butterfly::genMesh(std::string name) {
    if (name == "") {
        name = "output.off";
    }
    try {
        if (!OpenMesh::IO::write_mesh(mesh, name)) {
            std::cerr << "Cannot write mesh to file 'output.off'" << std::endl;
            return;
        }
    }
    catch (std::exception& e) {
        std::cerr << e.what() << std::endl;
        return;
    }

}

void Butterfly::getResult() {
    //if (getData()->edges.size() == 0) {
    //	getData()->edges.push_back(std::vector<V3f>());
    //}
    getData()->edges.clear();
    getData()->edges.push_back(std::vector<V3f>());
    for (auto e_it = mesh.edges_begin(); e_it != mesh.edges_end(); ++e_it)
    {
        // 得到边所代表的半边
        OpenMesh::HalfedgeHandle heh1 = mesh.halfedge_handle(*e_it, 0); // 默认一个方向的半边
        OpenMesh::Vec3d edgePoint(0, 0, 0);
        int edgePointsNumber = 0;
        OpenMesh::DefaultTraits::Point pointV = mesh.point(mesh.from_vertex_handle(heh1)); // 这条（半）边的起点
        OpenMesh::DefaultTraits::Point pointW = mesh.point(mesh.to_vertex_handle(heh1));   // 这条（半）边的终点
        getData()->edges[0].push_back({ pointV[0], pointV[1], pointV[2] });
        getData()->edges[0].push_back({ pointW[0], pointW[1], pointW[2] });
    }
}