Barycentric mapping

[petrasvestartas]

Is there any way to speed up the barycentric mapping?

For one mesh this method runs for around 3.5 seconds.

std::vector<std::vector<int>> map_mesh_cropped(
    Eigen::Ref<const compas::RowMatrixXd> v, 
    Eigen::Ref<const compas::RowMatrixXi> f, 
    Eigen::Ref<const compas::RowMatrixXd> uv,
    Eigen::Ref<compas::RowMatrixXd> pattern_v, 
    const std::vector<std::vector<int>>& pattern_f, 
    Eigen::Ref<const compas::RowMatrixXd> pattern_uv)
{
    // Use regular MatrixXd to avoid type conflicts with AABB functions
    Eigen::MatrixXd V_uv = uv;
    Eigen::MatrixXi F_faces = f;
    
    igl::AABB<Eigen::MatrixXd, 2> tree;
    tree.init(V_uv, F_faces);

    Eigen::MatrixXd C;
    Eigen::VectorXi I;
    Eigen::VectorXd sqrD;
    Eigen::MatrixXd pattern_uv_eigen = pattern_uv;
    
    // Find closest points on the mesh for all pattern vertices
    tree.squared_distance(V_uv, F_faces, pattern_uv_eigen, sqrD, I, C);
    
    // Time how long the barycentric mapping takes
    auto bary_start = std::chrono::high_resolution_clock::now();
    
    // Map each pattern vertex to 3D using barycentric coordinates
    for(int id = 0; id < pattern_uv.rows(); id++)
    {
        // Prepare point for barycentric coordinate calculation
        Eigen::MatrixXd P(1, 3);
        P << pattern_uv(id, 0), pattern_uv(id, 1), 0;
        
        // Inline get_triface for UV coordinates
        int faceID = I(id);
        Eigen::MatrixXd UV_A(1, 3), UV_B(1, 3), UV_C(1, 3);
        
        if(0 <= faceID && faceID < f.rows()) {
            // UV coordinates
            UV_A << uv(f(faceID, 0), 0), uv(f(faceID, 0), 1), 0;
            UV_B << uv(f(faceID, 1), 0), uv(f(faceID, 1), 1), 0;
            UV_C << uv(f(faceID, 2), 0), uv(f(faceID, 2), 1), 0;
        }
        
        // Calculate barycentric coordinates
        Eigen::MatrixXd L;
        igl::barycentric_coordinates(P, UV_A, UV_B, UV_C, L);
        
        // Inline get_triface for 3D vertex positions
        Eigen::MatrixXd A(1, 3), B(1, 3), C(1, 3);
        
        if(0 <= faceID && faceID < f.rows()) {
            // 3D coordinates
            A << v(f(faceID, 0), 0), v(f(faceID, 0), 1), v(f(faceID, 0), 2);
            B << v(f(faceID, 1), 0), v(f(faceID, 1), 1), v(f(faceID, 1), 2);
            C << v(f(faceID, 2), 0), v(f(faceID, 2), 1), v(f(faceID, 2), 2);
        }
        
        // Update pattern vertex position through barycentric interpolation
        pattern_v.row(id) = A * L(0, 0) + B * L(0, 1) + C * L(0, 2);
    }
    
    
    return pattern_f;
}