openmesh 用矩阵法映射到圆盘

转自：https://blog.csdn.net/a605907914/article/details/51416732

#include<iostream>  
#include <OpenMesh\Core\IO\MeshIO.hh>  
#include <OpenMesh\Core\Mesh\TriMesh_ArrayKernelT.hh>  
#include <OpenMesh\Core\Utils\Property.hh>
#include <Eigen/Dense>    
#include <Eigen/Sparse>
#include <map>
struct MyTraits : public OpenMesh::DefaultTraits
{
    VertexTraits
    {
        int some_additional_index;
    };
};

typedef OpenMesh::TriMesh_ArrayKernelT<MyTraits> MyMesh;

using namespace std;
using namespace Eigen;

const double inf = 0.0000001;
bool judge(float xx, float yy){
    //cout << xx - yy << endl;   
    if ((xx - yy) > -1 * inf && (xx - yy) < inf)return true;
    return false;
}
double get_cot(MyMesh::Point t1, MyMesh::Point t2){
    double p1 = t1.data()[0];
    double p2 = t1.data()[1];
    double p3 = t2.data()[0];
    double p4 = t2.data()[1];
    double f1 = (p1*p3 + p2*p4);
    double len_a = powf(t1.data()[0], 2) + powf(t1.data()[1], 2);
    double len_b = powf(t2.data()[0], 2) + powf(t2.data()[1], 2);
    double COS = f1 / sqrtf(len_a*len_b);
    double SIN=sqrtf(1.0-COS*COS);
    return COS / SIN;
}

bool solve_eigen(Eigen::SparseMatrix<double>&A, VectorXd&b, VectorXd&x){
    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
    solver.compute(A);
    if (solver.info() != Success){ return false; }
    x = solver.solve(b);
    return solver.info() == Success;
}

int main()
{
    MyMesh mesh;
    MyMesh mymesh;

    //mesh.request_vertex_texcoords2D();  

    if (!OpenMesh::IO::read_mesh(mesh, "meshes/alex.off")){
        cerr << "can't open the file." << endl;
        exit(1);
    }

    vector<MyMesh::HalfedgeHandle>h_handle;
    vector<MyMesh::VertexHandle>v_handle;
    vector<MyMesh::VertexHandle>v_in_handle;
    map<MyMesh::VertexHandle, int>V_Map;
    OpenMesh::VPropHandleT<bool>VH_Bool;
    OpenMesh::VPropHandleT<double>VH_Double;
    OpenMesh::VPropHandleT<MyMesh::TexCoord2D>VH_TC;
    OpenMesh::HPropHandleT<bool>HE_bool;
    //OpenMesh::HPropHandleT<double>
    mesh.add_property(VH_Bool, "VH_Bool");
    mesh.add_property(VH_Double, "VH_Double");
    mesh.add_property(HE_bool, "HE_bool");
    V_Map.clear();
    float PAI = acosf(-1.0);
    float r = 1.0;
    MyMesh::VertexHandle temp_vh_it;
    MyMesh::VertexHandle temp_mark;
    //int dd = 0;  
    for (MyMesh::VertexIter vh_it = mesh.vertices_begin(); vh_it != mesh.vertices_end(); ++vh_it){
        if (mesh.is_boundary(*vh_it)){
            temp_vh_it = *vh_it;
            break;
        }
    }
    //cout << temp_vh_it.idx() << endl;  

    //cout << "111111" << endl;  
    int uu = 0;
    int break_ll = 0;
    temp_mark = temp_vh_it;
    while (1){
        //if (break_ll == 1)break;  
        for (MyMesh::VertexOHalfedgeCWIter VOH_it = mesh.voh_cwbegin(temp_vh_it); VOH_it != mesh.voh_cwend(temp_vh_it); ++VOH_it){
            if (mesh.is_boundary(*VOH_it)){
                MyMesh::VertexHandle temp_in_v = mesh.to_vertex_handle(*VOH_it);
                if (mesh.property(VH_Bool, temp_in_v)){ break_ll = 1; break; }
                h_handle.push_back(*VOH_it);
                //v_handle.push_back(temp_in_v);
                
                temp_vh_it = temp_in_v;
                mesh.property(VH_Bool, temp_in_v) = true;
                mesh.property(VH_Bool).set_persistent(true);
                //uu++;  
                //cout << uu << endl;  
                //cout << temp_mark.idx() << " " << temp_vh_it.idx() << endl;  
                break;
            }


        }
        if (temp_mark == temp_vh_it)break;

        //cout <<"IDX = "<<temp_vh_it.idx() << endl;  
        //if (uu == 466)break;  
    }

    cout << "222222" << endl;  
    float sum = 0.0;
    float temp_sum = 0.0;
    for (auto it = h_handle.begin(); it != h_handle.end(); it++){
        MyMesh::VertexHandle VH_it_to = mesh.to_vertex_handle(*it);
        MyMesh::VertexHandle VH_it_from = mesh.from_vertex_handle(*it);
        MyMesh::Point P_it_to = mesh.point(VH_it_to);
        MyMesh::Point P_it_from = mesh.point(VH_it_from);
        temp_sum = powf(P_it_from.data()[0] - P_it_to.data()[0], 2.0) + powf(P_it_from.data()[1] - P_it_to.data()[1], 2.0) + powf(P_it_from.data()[2] - P_it_to.data()[2], 2.0);
        temp_sum = sqrt(temp_sum);
        sum += temp_sum;
        mesh.property(VH_Double, VH_it_to) = sum;
    }
    cout << "333333" << endl;  

    mesh.request_vertex_texcoords2D();
    MyMesh::TexCoord2D;


    for (auto it = mesh.vertices_begin(); it != mesh.vertices_end(); it++){
        if (mesh.is_boundary(mesh.halfedge_handle(*it))){
            float arf = (mesh.property(VH_Double, *it)*PAI * 2 / sum);
            float point_x = cosf(arf)*r;
            float point_y = sinf(arf)*r;
            mesh.set_texcoord2D(*it, OpenMesh::Vec2f(point_x, point_y));
        }
        else{
            v_in_handle.push_back(*it);
            V_Map[*it] = v_in_handle.size() - 1;
            mesh.set_texcoord2D(*it, OpenMesh::Vec2f(0.0, 0.0));
        }
        //cout << (*it).idx() << endl;  
    }
    
    int pp = 0;
    cout << "Begin computing!!!!" << endl;  

    int in_size = v_in_handle.size();
    cout << "In_Point = " << in_size << endl;
    vector<Eigen::Triplet<double>>A_coefficient;
    Eigen::SparseMatrix<double> A(in_size, in_size);
    Eigen::VectorXd U(in_size);
    Eigen::VectorXd V(in_size);
    for (int ii = 0; ii < v_in_handle.size(); ii++){
        MyMesh::VertexHandle temp_v = v_in_handle[ii];
        double Kij = 0.0;
        double Kij_U = 0.0;
        double Kij_V = 0.0;
        for (auto it = mesh.voh_ccwbegin(temp_v); it != mesh.voh_ccwend(temp_v); it++){                                                 //out in ?
            
            //MyMesh::HalfedgeHandle temp_he_handle = mesh.halfedge_handle(*it, 0);
            MyMesh::HalfedgeHandle temp_he_op_handle = mesh.opposite_halfedge_handle(*it);
            MyMesh::HalfedgeHandle temp_he_handle_next = mesh.next_halfedge_handle(*it);
            MyMesh::HalfedgeHandle temp_he_op_handle_next = mesh.next_halfedge_handle(temp_he_op_handle);

            MyMesh::VertexHandle temp_is_to = mesh.to_vertex_handle(*it);
            MyMesh::VertexHandle temp_is_from = mesh.from_vertex_handle(*it);
            MyMesh::VertexHandle temp_is_third = mesh.to_vertex_handle(temp_he_handle_next);
            

            MyMesh::VertexHandle temp_is_op_to = mesh.to_vertex_handle(temp_he_op_handle);
            MyMesh::VertexHandle temp_is_op_from = mesh.from_vertex_handle(temp_he_op_handle);
            MyMesh::VertexHandle temp_is_op_third = mesh.to_vertex_handle(temp_he_op_handle_next);

            MyMesh::Point temp1_point = mesh.point(temp_is_to) - mesh.point(temp_is_third);
            MyMesh::Point temp2_point = mesh.point(temp_is_from) - mesh.point(temp_is_third);
            double cot1 = get_cot(temp1_point, temp2_point);
            MyMesh::Point temp1_op_point = mesh.point(temp_is_op_to) - mesh.point(temp_is_op_third);
            MyMesh::Point temp2_op_point = mesh.point(temp_is_op_from) - mesh.point(temp_is_op_third);
            double cot2 = get_cot(temp1_op_point, temp2_op_point);
            double temp_Kij = (cot1 + cot2)/2 ;
            Kij -= temp_Kij;
            if (mesh.is_boundary(mesh.to_vertex_handle(*it))){
                Kij_U -= temp_Kij*mesh.texcoord2D(mesh.to_vertex_handle(*it)).data()[0];
                Kij_V -= temp_Kij*mesh.texcoord2D(mesh.to_vertex_handle(*it)).data()[1];
            }
            else{
                A_coefficient.push_back(Eigen::Triplet<double>(ii, V_Map[mesh.to_vertex_handle(*it)], temp_Kij));
            }
        }
        U(ii) = Kij_U;
        V(ii) = Kij_V;
        A_coefficient.push_back(Eigen::Triplet<double>(ii, ii, Kij));
    }
    A.setFromTriplets(A_coefficient.begin(), A_coefficient.end());
    Eigen::VectorXd x(in_size);
    Eigen::VectorXd y(in_size);
    solve_eigen(A,U,x);
    solve_eigen(A,V,y);

    for (int ii = 0; ii < x.size(); ii++){
        mesh.set_texcoord2D(v_in_handle[ii],OpenMesh::Vec2f(x[ii],y[ii]));
    }

    
    /*
    while (1){
        int mark_break = 0;
        pp = 0;
        for (auto it = mesh.vertices_begin(); it != mesh.vertices_end(); it++){
            if (mesh.is_boundary(*it))continue;
            else {
                float xx_sum = 0.0;
                float yy_sum = 0.0;
                float xx_pre = 0.0;
                float yy_pre = 0.0;
                float num = 0.0;
                for (auto v_it = mesh.vv_begin(*it); v_it != mesh.vv_end(*it); v_it++){
                    float xx = mesh.texcoord2D(*v_it).data()[0];
                    float yy = mesh.texcoord2D(*v_it).data()[1];
                    xx_sum += xx;
                    yy_sum += yy;
                    num++;
                }
                if (judge(num, 0.0))continue;
                xx_pre = mesh.texcoord2D(*it).data()[0];
                yy_pre = mesh.texcoord2D(*it).data()[1];
                yy_sum /= num;
                xx_sum /= num;
                float f_sum = (xx_pre - xx_sum)*(xx_pre - xx_sum) + (yy_pre - yy_sum)*(yy_pre - yy_sum);
                f_sum = sqrtf(f_sum);
                if (!judge(f_sum, 0.0)){ mark_break = 1; pp++; }
                mesh.set_texcoord2D(*it, OpenMesh::Vec2f(xx_sum, yy_sum));
            }
        }
        if (!mark_break)break;
        //if (pp%1000==0)  
        //cout << pp << endl;  
        //if (pp <= 63000)cout << pp << endl;  
    }
    */
    


    cout << "Output Answer!!!!" << endl;

    mesh.property(VH_Bool).set_persistent(true);
    mesh.property(VH_Double).set_persistent(true);
    //mesh.release_vertex_texcoords2D();  
    OpenMesh::IO::Options opt = OpenMesh::IO::Options::VertexTexCoord;

    if (!OpenMesh::IO::write_mesh(mesh, "gg.off", opt)){
        cerr << "Can't write the mesh." << endl;
        cin.get();

    }
    cout << "YOU GET THE ANSWER!!!!" << endl;

    cin.get();
    return 0;
}