

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:
[20200417 18:10:31.988] [cpf] [info] Brute force id = 386 [20200417 18:10:31.988] [cpf] [info] CGAL id = 84 [20200417 18:10:31.988] [cpf] [info] Brute force distance = 0.0 [20200417 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?


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:
>
>
> [20200417 18:10:31.988] [cpf] [info] Brute force id = 386
>
> [20200417 18:10:31.988] [cpf] [info] CGAL id = 84
>
> [20200417 18:10:31.988] [cpf] [info] Brute force distance = 0.0
>
> [20200417 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 cgaldiscuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgaldiscuss


Sure, find attached a minimal example and the data I am using. El vie., 17 abr. 2020 a las 18:41, Sebastien Loriot (GeometryFactory) (< [hidden email]>) 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:
>
>
> [20200417 18:10:31.988] [cpf] [info] Brute force id = 386
>
> [20200417 18:10:31.988] [cpf] [info] CGAL id = 84
>
> [20200417 18:10:31.988] [cpf] [info] Brute force distance = 0.0
>
> [20200417 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 cgaldiscuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgaldiscuss


On Fri, 17 Apr 2020, Juan Jose Casafranca wrote:
> Sure, find attached a minimal example and the data I am using.
It is far from minimal, there are several unnecessary #includes at least.
Your property map looks suspicious, with operator[] that returns a value
but reference that is a typedef for a reference type, and a category that
implies a reference.
Did you try fsanitize=address?

Marc Glisse

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


So if I understand correctly, the problem here is that the property map must model an LValuePropertyMap and that concept requires th operator[] to return a reference, right? I have tried storing a vector of K::Point_3 points and works fine. Is there anyway I can build the K::Point_3 when needed, instead, similarly to how I am doing it right now?
Thanks El vie., 17 abr. 2020 a las 21:03, Marc Glisse (< [hidden email]>) escribió: On Fri, 17 Apr 2020, Juan Jose Casafranca wrote:
> Sure, find attached a minimal example and the data I am using.
It is far from minimal, there are several unnecessary #includes at least.
Your property map looks suspicious, with operator[] that returns a value
but reference that is a typedef for a reference type, and a category that
implies a reference.
Did you try fsanitize=address?

Marc Glisse

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


Hi Are You using any query algorithm for find the closest neighbor point. Thanks On Sat, Apr 18, 2020 at 1:08 AM Juan Jose Casafranca < [hidden email]> wrote: So if I understand correctly, the problem here is that the property map must model an LValuePropertyMap and that concept requires th operator[] to return a reference, right? I have tried storing a vector of K::Point_3 points and works fine. Is there anyway I can build the K::Point_3 when needed, instead, similarly to how I am doing it right now?
Thanks
El vie., 17 abr. 2020 a las 21:03, Marc Glisse (< [hidden email]>) escribió: On Fri, 17 Apr 2020, Juan Jose Casafranca wrote:
> Sure, find attached a minimal example and the data I am using.
It is far from minimal, there are several unnecessary #includes at least.
Your property map looks suspicious, with operator[] that returns a value
but reference that is a typedef for a reference type, and a category that
implies a reference.
Did you try fsanitize=address?

Marc Glisse

You are currently subscribed to cgaldiscuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgaldiscuss
 With Warm regards
Lokesh kumar


In reply to this post by Juan Jose Casafranca
On Fri, 17 Apr 2020, Juan Jose Casafranca wrote:
> So if I understand correctly, the problem here is that the property map
> must model an LValuePropertyMap and that concept requires th operator[] to
> return a reference, right? I have tried storing a vector of K::Point_3
> points and works fine. Is there anyway I can build the K::Point_3 when
> needed, instead, similarly to how I am doing it right now?
I am not sure, but LValuePropertyMap might be ok with some proxy type
(like vector<bool> uses), or not. Also, the doc says LValuePropertyMap is
required, but the code actually contains a fallback for more general
property maps, so you could try that.
It may be easier to work with a single point type. If you cannot use
CGAL's points everywhere, you could instead provide traits that define
suitable functions for MyPoint.

Marc Glisse

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


I have finally decide to use an LValuePropertyMap and create a std::vector<K::Point_3>, at least for now. Maybe on the future I will create the appropriate traits for the other Point type I am using. Thanks, it would had been really complex to debug this without your help. El sáb., 18 abr. 2020 a las 9:05, Marc Glisse (< [hidden email]>) escribió: On Fri, 17 Apr 2020, Juan Jose Casafranca wrote:
> So if I understand correctly, the problem here is that the property map
> must model an LValuePropertyMap and that concept requires th operator[] to
> return a reference, right? I have tried storing a vector of K::Point_3
> points and works fine. Is there anyway I can build the K::Point_3 when
> needed, instead, similarly to how I am doing it right now?
I am not sure, but LValuePropertyMap might be ok with some proxy type
(like vector<bool> uses), or not. Also, the doc says LValuePropertyMap is
required, but the code actually contains a fallback for more general
property maps, so you could try that.
It may be easier to work with a single point type. If you cannot use
CGAL's points everywhere, you could instead provide traits that define
suitable functions for MyPoint.

Marc Glisse

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

