CGAL users discussion list

[Juan Jose Casafranca <>]

Sure, find attached a minimal example and the data I am using.

Thanks

El vie., 17 abr. 2020 a las 18:41, Sebastien Loriot (GeometryFactory) (<>) escribió:

Could you please provide a minimal example that we could compile, run
and that would be showing the error?

At first sight, I don't see anything suspicious.

Thanks,

Sebastien.

On 4/17/20 6:25 PM, Juan Jose Casafranca wrote:
> I am trying to find the closest neighbor in a point cloud to a given
> point. I have followed the documentation example in which the points are
> stored in a user vector and the Tree only contains indices. Since the
> results were not what I was expecting, I am also computing the closest
> neighbor by brute force. Ill explain my code:
>
> I have a vector of points (coarsePoints) and I have a vector of indices
> (xC_vertices). The indices are the points in coarsePoints that can be
> marked as closest neighbor.
>
> This is how I compute the closest point using a brute force algorithm
> ```
>
> for (int i = 0; i < d_vertices.size(); ++i) {
>
> Real minDistance = std::numeric_limits<Real>::infinity();
>
> int vertexId = -1;
>
> const auto &dPoint = highResPoints[d_vertices[i]];
>
> const CPF::Point_3 pD{dPoint[0], dPoint[1], dPoint[2]};
>
> for (int j = 0; j < xC_vertices.size(); ++j) {
>
> const auto &xPoint = coarsePoints[xC_vertices[j]];
>
> const CPF::Point_3 pC{xPoint[0], xPoint[1], xPoint[2]};
>
>
> if (CGAL::squared_distance(pC, pD) < minDistance) {
>
> vertexId = j;
>
> minDistance = CGAL::squared_distance(pC, pD);
>
> }
>
> }
>
> ```
>
>
> When using CGAL, I have the following (based on the documentation example):
>
>
> ```
>
> class My_point_property_map
>
> {
>
> const sofa::helper::vector<CPF::SofaTypes::Point> &points;
>
> const std::vector<CPF::SofaTypes::Tetra::value_type> &xC_vertices;
>
> public:
>
> typedef CPF::Point_3 value_type;
>
> typedef value_type &reference;
>
> typedef std::size_t key_type;
>
> typedef boost::lvalue_property_map_tag category;
>
> My_point_property_map(const sofa::helper::vector<CPF::SofaTypes::Point>
> &pts,
>
> const std::vector<CPF::SofaTypes::Tetra::value_type> &xC_vertices_)
>
> : points(pts)
>
> , xC_vertices(xC_vertices_)
>
> {
>
> }
>
> value_type operator[](key_type k) const
>
> {
>
> const CPF::SofaTypes::Point &sPoint = points[xC_vertices[k]];
>
> return CPF::Point_3{sPoint[0], sPoint[1], sPoint[2]};
>
> }
>
> friend value_type get(const My_point_property_map &ppmap, key_type i) {
> return ppmap[i]; }
>
> };
>
>
> using SearchTraitsBase = CGAL::Search_traits_3<CPF::K>;
>
> using SearchTraits = CGAL::Search_traits_adapter<std::size_t,
> My_point_property_map, SearchTraitsBase>;
>
> using NeighborSearch = CGAL::Orthogonal_k_neighbor_search<SearchTraits>;
>
> using Tree = NeighborSearch::Tree;
>
> using Splitter = typename Tree::Splitter;
>
> using SearchDistance = typename NeighborSearch::Distance;
>
>
> int searchNeighbor() {
>
> My_point_property_map ppmap(coarsePoints, xC_vertices);
>
> Tree searchTree(boost::counting_iterator<int>(0),
>
> boost::counting_iterator<int>(xC_vertices.size()),
>
> Splitter(),
>
> SearchTraits(ppmap));
>
> SearchDistance trDistance(ppmap);
>
>
> for (int i = 0; i < d_vertices.size(); ++i) {
>
> NeighborSearch search(searchTree, CPF::Point_3{dPoint[0], dPoint[1],
> dPoint[2]}, 1, 0, true, trDistance);
>
> for (const auto &searchResult : search) {
>
> if (searchResult.first != vertexId) {
>
> spdlog::get("cpf")->info("Brute force id = {}", vertexId);
>
> spdlog::get("cpf")->info("CGAL id = {}", searchResult.first);
>
> spdlog::get("cpf")->info("Brute force distance = {}", minDistance);
>
> spdlog::get("cpf")->info("CGAL distance = {}", searchResult.second);
>
> }
>
> }
>
> } }
>
> ```
>
>
> This return a los of vertices that are actually different:
>
>
> [2020-04-17 18:10:31.988] [cpf] [info] Brute force id = 386
>
> [2020-04-17 18:10:31.988] [cpf] [info] CGAL id = 84
>
> [2020-04-17 18:10:31.988] [cpf] [info] Brute force distance = 0.0
>
> [2020-04-17 18:10:31.988] [cpf] [info] CGAL distance = 0.041339367896971506
>
>
> I am not sure if I am using incorrectly the NeighborSearch or what is
> happening. Any ideas what might be wrong?
>
>
>
>

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/Search_traits_3.h>
#include <CGAL/Search_traits_adapter.h>
#include <boost/iterator/counting_iterator.hpp>
#include <CGAL/AABB_face_graph_triangle_primitive.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/Polygon_mesh_processing/locate.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>

#include <fstream>

using Real = double;
using K = CGAL::Simple_cartesian<Real>;
using Point_3 = K::Point_3;
using MyPoint = std::array<Real, 3>;
using PointId = unsigned int;

class My_point_property_map
{
    const std::vector<MyPoint> &points;
    const std::vector<PointId> &xC_vertices;

public:
    typedef K::Point_3 value_type;
    typedef value_type &reference;
    typedef std::size_t key_type;
    typedef boost::lvalue_property_map_tag category;
    My_point_property_map(const std::vector<MyPoint> &pts, const std::vector<PointId> &xC_vertices_)
        : points(pts)
        , xC_vertices(xC_vertices_)
    {
    }
    value_type operator[](key_type k) const
    {
        const MyPoint &sPoint = points[xC_vertices[k]];
        return Point_3{sPoint[0], sPoint[1], sPoint[2]};
    }
    friend value_type get(const My_point_property_map &ppmap, key_type i) { return ppmap[i]; }
};

using SearchTraitsBase = CGAL::Search_traits_3<K>;
using SearchTraits = CGAL::Search_traits_adapter<std::size_t, My_point_property_map, SearchTraitsBase>;
using NeighborSearch = CGAL::Orthogonal_k_neighbor_search<SearchTraits>;
using Tree = NeighborSearch::Tree;
using Splitter = typename Tree::Splitter;
using SearchDistance = typename NeighborSearch::Distance;

int main()
{
    std::vector<PointId> d_vertices;
    std::vector<PointId> xC_vertices;

    std::vector<MyPoint> coarsePoints;
    std::vector<MyPoint> highResPoints;

    std::ifstream str;
    str.open("/tmp/test.bin", std::istream::binary);
    std::size_t coarsePointsSize;
    std::size_t highResPointsSize;
    std::size_t xC_verticesSize;
    std::size_t d_verticesSize;
    str >> coarsePointsSize;
    coarsePoints.resize(coarsePointsSize);
    str.read(reinterpret_cast<char *>(coarsePoints.data()), sizeof(MyPoint) * coarsePoints.size());
    str >> highResPointsSize;
    highResPoints.resize(highResPointsSize);
    str.read(reinterpret_cast<char *>(highResPoints.data()), sizeof(MyPoint) * highResPoints.size());
    str >> xC_verticesSize;
    xC_vertices.resize(xC_verticesSize);
    str.read(reinterpret_cast<char *>(xC_vertices.data()), sizeof(unsigned int) * xC_vertices.size());
    str >> d_verticesSize;
    d_vertices.resize(d_verticesSize);
    str.read(reinterpret_cast<char *>(d_vertices.data()), sizeof(unsigned int) * d_vertices.size());

    My_point_property_map ppmap(coarsePoints, xC_vertices);
    Tree searchTree(boost::counting_iterator<int>(0),
                    boost::counting_iterator<int>(xC_vertices.size()),
                    Splitter(),
                    SearchTraits(ppmap));

    SearchDistance trDistance(ppmap);

    for (int i = 0; i < d_vertices.size(); ++i) {
        Real minDistance = std::numeric_limits<Real>::infinity();
        int vertexId = -1;
        const auto &dPoint = highResPoints[d_vertices[i]];
        const Point_3 pD{dPoint[0], dPoint[1], dPoint[2]};
        for (int j = 0; j < xC_vertices.size(); ++j) {
            const auto &xPoint = coarsePoints[xC_vertices[j]];
            const Point_3 pC{xPoint[0], xPoint[1], xPoint[2]};

            const auto distance = trDistance.transformed_distance(pD, j);
            if (distance < minDistance) {
                vertexId = j;
                minDistance = distance;
            }
        }

        NeighborSearch search(searchTree, Point_3{dPoint[0], dPoint[1], dPoint[2]}, 1, 0.0001, true, trDistance);
        for (const auto &searchResult : search) {
            if (searchResult.first != vertexId) {
                std::cout << "Brute force id = " << vertexId << "\n";
                std::cout << "CGAL id = " << searchResult.first << "\n";
                std::cout << "Brute force distance = " << minDistance << "\n";
                std::cout << "CGAL distance = " << searchResult.second << "\n";
            }
        }
    }
}

Attachment: test.bin
Description: Binary data