Issue Details
Surface_mesh_shortest_path::locate() performs an inexact construction under
Exact_predicates_inexact_constructions_kernel, so for a query point that lies exactly on a mesh
vertex it returns barycentric coordinates carrying round-off — typically one component at
-5.6e-17 instead of 0. Surface_mesh_shortest_paths_3::Classify_barycentric_coordinates then
tests those coordinates exactly, so the location is not classified as
BARYCENTRIC_COORDINATES_ON_VERTEX, and expand_root() roots the wavefront somewhere other than
the intended vertex.
The failure is silent — no exception, no assertion, a plausible-looking field — and the affected
vertex reports a distance to its own source point of exactly one edge length instead of zero.
Adding the same vertex via add_source_point(vertex_descriptor) is correct, because
face_location(v) produces exactly (0, 0, 1). So the two documented ways of seeding the same
geometric point disagree, and the composition used by the shortest_path_with_locate example is
the one that is wrong.
It is intermittent with respect to mesh size on identical topology, since it turns on the last bit
of the projection. In the program below, a negative component in the middle slot mis-roots;
n = 8 has a negative in the first slot and survives.
Source Code
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Surface_mesh_shortest_path.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits_3.h>
#include <CGAL/AABB_face_graph_triangle_primitive.h>
#include <iostream>
#include <iomanip>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef CGAL::Surface_mesh<Kernel::Point_3> Mesh;
typedef CGAL::Surface_mesh_shortest_path_traits<Kernel, Mesh> Traits;
typedef CGAL::Surface_mesh_shortest_path<Traits> ShortestPath;
typedef CGAL::AABB_face_graph_triangle_primitive<Mesh> Primitive;
typedef CGAL::AABB_traits_3<Kernel, Primitive> AABBTraits;
typedef CGAL::AABB_tree<AABBTraits> AABBTree;
typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
int main()
{
for (int n : {4, 5, 6, 7, 8, 9})
{
// planar unit square, n x n vertices, two triangles per cell
Mesh mesh;
std::vector<Mesh::Vertex_index> vi;
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j)
vi.push_back(mesh.add_vertex(Kernel::Point_3(double(i)/(n-1), double(j)/(n-1), 0.0)));
for (int i = 0; i < n - 1; ++i)
for (int j = 0; j < n - 1; ++j)
{
int a = i*n+j, b = i*n+j+1, c = (i+1)*n+j, d = (i+1)*n+j+1;
mesh.add_face(vi[a], vi[b], vi[d]);
mesh.add_face(vi[a], vi[d], vi[c]);
}
const vertex_descriptor vd(static_cast<std::size_t>(n - 1)); // corner (0, 1, 0)
ShortestPath by_vertex(mesh);
by_vertex.add_source_point(vd);
const double d_vertex = by_vertex.shortest_distance_to_source_points(vd).first;
ShortestPath by_point(mesh);
AABBTree tree;
by_point.build_aabb_tree(tree);
const auto loc_v = by_vertex.face_location(vd);
const auto loc_p = by_point.locate(mesh.point(vd), tree); // the SAME point
by_point.add_source_point(loc_p);
const double d_point = by_point.shortest_distance_to_source_points(vd).first;
std::cout << "n=" << n
<< " d(vertex source)=" << std::setw(8) << d_vertex
<< " d(point source)=" << std::setw(8) << d_point
<< (d_point > 1e-12 ? " <-- WRONG" : "")
<< "\n face_location bary = " << std::setprecision(17)
<< loc_v.second[0] << ", " << loc_v.second[1] << ", " << loc_v.second[2]
<< "\n locate() bary = "
<< loc_p.second[0] << ", " << loc_p.second[1] << ", " << loc_p.second[2]
<< std::setprecision(6) << "\n";
}
return 0;
}Built header-only:
c++ -std=c++17 -O3 -DNDEBUG -DCGAL_HEADER_ONLY -DCGAL_DISABLE_GMP -DCGAL_USE_BOOST_MP \
-I <cgal>/include -I <boost> -I <eigen> smsp_locate_bug.cpp -o smsp_locate_bug
Actual output
n=4 d(vertex source)= 0 d(point source)=0.333333 <-- WRONG
face_location bary = 0, 0, 1
locate() bary = 0, -5.5511151231257821e-17, 1
n=5 d(vertex source)= 0 d(point source)= 0
face_location bary = 0, 0, 1
locate() bary = 0, 0, 1
n=6 d(vertex source)= 0 d(point source)= 0.2 <-- WRONG
face_location bary = 0, 0, 1
locate() bary = 0, -6.5919492087118657e-17, 1
n=7 d(vertex source)= 0 d(point source)=0.166667 <-- WRONG
face_location bary = 0, 0, 1
locate() bary = 0, -5.5511151231257864e-17, 1
n=8 d(vertex source)= 0 d(point source)=1.38778e-17
face_location bary = 0, 0, 1
locate() bary = -1.1102230246251565e-16, 6.09863722023096e-17, 1
n=9 d(vertex source)= 0 d(point source)= 0
face_location bary = 0, 0, 1
locate() bary = 0, 0, 1
Expected: d(point source) == 0 in every case, since the source point is the queried vertex. The
wrong value is always exactly one edge length, 1/(n-1).
Suggested direction
The composition of an inexact construction with an exact classification is the root of it, so the
repair seems to belong between the two rather than in either alone — for instance normalising a
Face_location onto the simplex (clamping sub-ulp components to zero and renormalising) inside
add_source_point(const Face_location&), or having locate() return coordinates already snapped.
Barycentric coordinates are dimensionless and bounded in [0, 1], so such a threshold is
scale-free rather than a tolerance on a length.
As a downstream workaround we snap components below 8 * DBL_EPSILON to zero and renormalise before
calling add_source_point, which makes all six cases above return zero and leaves face-interior
sources untouched.
Environment
- CGAL 6.1.1 (release tarball, header-only)
- Boost 1.82.0, Eigen 3.4.0
- Apple clang, macOS arm64,
-O3 -DNDEBUG
Issue Details
Surface_mesh_shortest_path::locate()performs an inexact construction underExact_predicates_inexact_constructions_kernel, so for a query point that lies exactly on a meshvertex it returns barycentric coordinates carrying round-off — typically one component at
-5.6e-17instead of0.Surface_mesh_shortest_paths_3::Classify_barycentric_coordinatesthentests those coordinates exactly, so the location is not classified as
BARYCENTRIC_COORDINATES_ON_VERTEX, andexpand_root()roots the wavefront somewhere other thanthe intended vertex.
The failure is silent — no exception, no assertion, a plausible-looking field — and the affected
vertex reports a distance to its own source point of exactly one edge length instead of zero.
Adding the same vertex via
add_source_point(vertex_descriptor)is correct, becauseface_location(v)produces exactly(0, 0, 1). So the two documented ways of seeding the samegeometric point disagree, and the composition used by the
shortest_path_with_locateexample isthe one that is wrong.
It is intermittent with respect to mesh size on identical topology, since it turns on the last bit
of the projection. In the program below, a negative component in the middle slot mis-roots;
n = 8has a negative in the first slot and survives.Source Code
Built header-only:
Actual output
Expected:
d(point source) == 0in every case, since the source point is the queried vertex. Thewrong value is always exactly one edge length,
1/(n-1).Suggested direction
The composition of an inexact construction with an exact classification is the root of it, so the
repair seems to belong between the two rather than in either alone — for instance normalising a
Face_locationonto the simplex (clamping sub-ulp components to zero and renormalising) insideadd_source_point(const Face_location&), or havinglocate()return coordinates already snapped.Barycentric coordinates are dimensionless and bounded in
[0, 1], so such a threshold isscale-free rather than a tolerance on a length.
As a downstream workaround we snap components below
8 * DBL_EPSILONto zero and renormalise beforecalling
add_source_point, which makes all six cases above return zero and leaves face-interiorsources untouched.
Environment
-O3 -DNDEBUG