CGAL locate with AABB tree returning wrong results

Question

I am using CGAL to perform some location queries. In particular, I want to verify in which triangle or a CGAL::Surface_mesh a point is, if any.

To do that, I am using the polygon mesh processing package. I have built an AABB tree to accelerate this queries and depending on the result of the query, I return an optional with the index of the face or nullopt.

class PMPTriangleIntersection
{
public:
    using Mesh = CGAL::Surface_mesh<VNCS::Sim2D::Kernel::Point>;
    using Point = VNCS::Sim2D::Kernel::Point;
    using Point3 = VNCS::Sim2D::Kernel::K::Point_3;
    using Real = VNCS::Sim2D::Kernel::Real;

    struct Point3VPM {
        using key_type = Mesh::Vertex_index;
        using value_type = Point3;
        using reference = value_type;
        using category = boost::readable_property_map_tag;

        Point3VPM();

        Point3VPM(const Mesh &m);

        friend Point3VPM::value_type get(const Point3VPM &map, Mesh::Vertex_index idx)
        {
            Point p = map.meshPtr->point(idx);
            return {p[0], p[1], 0};
        }

        const Mesh *meshPtr;
    };

    using AABBTreeTraits = CGAL::AABB_traits<VNCS::Sim2D::Kernel::K,  //
                                             CGAL::AABB_face_graph_triangle_primitive<Mesh, Point3VPM>>;
    using AABBTree = CGAL::AABB_tree<AABBTreeTraits>;

    PMPTriangleIntersection(const Mesh &mesh);

    std::optional<Mesh::Face_index> query(const Point &p) const;

private:
    std::reference_wrapper<const Mesh> m_mesh;
    AABBTree m_tree;
};

namespace PMP = CGAL::Polygon_mesh_processing;

VNCS::Sim2D::PMPTriangleIntersection::PMPTriangleIntersection(const Mesh &mesh)
    : m_mesh(std::cref(mesh))
{
    Point3VPM vpm(m_mesh.get());

    PMP::build_AABB_tree(m_mesh.get(), m_tree, CGAL::parameters::vertex_point_map(vpm));
}

std::optional<VNCS::Sim2D::PMPTriangleIntersection::Mesh::Face_index> VNCS::Sim2D::PMPTriangleIntersection::query(
    const Point &p) const
{
    auto location = PMP::locate_with_AABB_tree(p, m_tree, m_mesh.get());
    std::cout << p << "\n";
    std::cout << location.first.idx() << "\n";
    std::cout << location.second[0] << "\n";
    std::cout << location.second[1] << "\n";
    std::cout << location.second[2] << "\n";

    // Build triangle
    const auto vIndices = m_mesh.get().vertices_around_face(m_mesh.get().halfedge(location.first));
    const auto v =
        vIndices | ranges::views::transform([this](const auto v) { return m_mesh.get().point(v); }) | ranges::to_vector;

    std::cout << v[0] << "\n";
    std::cout << v[1] << "\n";
    std::cout << v[2] << "\n";
    std::cout << std::endl;

    if (ranges::all_of(location.second, [](const auto v) { return v > 0.0; }))
        return location.first;
    return {};
}

VNCS::Sim2D::PMPTriangleIntersection::Point3VPM::Point3VPM()
    : meshPtr(nullptr)
{
}

VNCS::Sim2D::PMPTriangleIntersection::Point3VPM::Point3VPM(const Mesh &m)
    : meshPtr(&m)
{
}

The issue is that for a particular point that is outside the triangle mesh, it is returning wrong results. In particular, the output of all the couts is the following:

-0.760221 -0.145442
17
0.5
0.5
0
-0.5 -0.25
-0.0334291 -0.0134272
-0.634804 0.105392

I have verified that my mesh and the query point are correct in paraview

Face:

Point:

I dont understand why I get this result for the barycentric coordinates. I would expect on of them to be negative or above 1.0.

Any hints?

Answer 1

Since point are expected to be closed to the input mesh (on the mesh up to rounding errors), locate function returns the location of the closest point onto the mesh from the query. Since you have a 2D mesh and you take a query outside of the mesh, then the closest feature is an edge. So locate will return a face on the boundary and one out of the three barycentric coordinates will be 0.