fix(Geometry): now unique points from ColocatedInfo have distances wider than given epsilon from each other

Author: panquez

include/geode/geometry/distance.hpp

@@ -23,6 +23,8 @@
 
 #pragma once
 
+#include <optional>
+
 #include <geode/geometry/common.hpp>
 
 namespace geode

src/geode/geometry/nn_search.cpp

@@ -23,6 +23,7 @@
 
 #include <geode/geometry/nn_search.hpp>
 
+#include <mutex>
 #include <numeric>
 
 #include <absl/algorithm/container.h>
@@ -34,6 +35,25 @@
 #include <geode/basic/logger.hpp>
 #include <geode/basic/pimpl_impl.hpp>
 
+namespace
+{
+    geode::index_t find_min_unmapped_element(
+        absl::Span< const geode::index_t > elements,
+        absl::Span< const geode::index_t > mapping )
+    {
+        geode::index_t min_index{ geode::NO_ID };
+        for( const auto e : elements )
+        {
+            if( mapping[e] != e )
+            {
+                continue;
+            }
+            min_index = std::min( min_index, e );
+        }
+        return min_index;
+    }
+} // namespace
+
 namespace geode
 {
     template < index_t dimension >
@@ -190,16 +210,22 @@ namespace geode
         typename NNSearch< dimension >::ColocatedInfo result;
         std::vector< index_t > mapping( nb_points() );
         absl::c_iota( mapping, 0 );
+        std::mutex mutex;
         async::parallel_for( async::irange( index_t{ 0 }, nb_points() ),
-            [&epsilon, &mapping, this]( index_t p ) {
+            [&epsilon, &mapping, &mutex, this]( index_t p ) {
                 if( mapping[p] == p )
                 {
                     const auto vertices =
                         radius_neighbors( point( p ), epsilon );
-                    const auto min_index = *absl::c_min_element( vertices );
+                    const auto min_index =
+                        find_min_unmapped_element( vertices, mapping );
+                    std::lock_guard< std::mutex > lock( mutex );
                     for( const auto id : vertices )
                     {
-                        mapping[id] = min_index;
+                        if( id == mapping[id] )
+                        {
+                            mapping[id] = min_index;
+                        }
                     }
                 }
             } );

tests/data/surface_degen.og_tsf3d

@@ -23,12 +23,13 @@
 
 #include <geode/basic/logger.hpp>
 
+#include <geode/geometry/distance.hpp>
 #include <geode/geometry/nn_search.hpp>
 #include <geode/geometry/point.hpp>
 
 #include <geode/tests/common.hpp>
 
-void test()
+void first_test()
 {
     const geode::NNSearch2D search{ { geode::Point2D{ { 0.1, 4.2 } },
         geode::Point2D{ { 5.9, 7.3 } }, geode::Point2D{ { 1.8, -5 } },
@@ -70,4 +71,57 @@ void test()
         "[Test] Error in unique points" );
 }
 
+void second_test()
+{
+    geode::Logger::set_level( geode::Logger::LEVEL::debug );
+    std::vector< geode::Point2D > points{
+        geode::Point2D{ { 0, 0.2 } },
+        geode::Point2D{ { 0.25, 0 } },
+        geode::Point2D{ { -0.25, 0 } },
+        geode::Point2D{ { 0, 0.8 } },
+        geode::Point2D{ { 0.25, 1 } },
+        geode::Point2D{ { -0.25, 1 } },
+    };
+    const geode::NNSearch2D colocator{ points };
+    const double DISTANCE{ 0.8 };
+
+    const auto colocated_info = colocator.colocated_index_mapping( DISTANCE );
+    OPENGEODE_EXCEPTION( colocated_info.nb_unique_points() == 2,
+        "[Test] Should be 2 unique points" );
+    for( const auto p : geode::Indices{ points } )
+    {
+        OPENGEODE_EXCEPTION( colocated_info.colocated_mapping[p]
+                                 < colocated_info.unique_points.size(),
+            "[Test] Wrong value of colocated_mapping (bigger than unique "
+            "points size)" );
+        const auto& colocated_point =
+            colocated_info.unique_points[colocated_info.colocated_mapping[p]];
+        OPENGEODE_EXCEPTION(
+            geode::point_point_distance( points[p], colocated_point )
+                <= DISTANCE,
+            "[Test] Colocated point is not close enough to original point" );
+    }
+    for( const auto up0 : geode::Indices{ colocated_info.unique_points } )
+    {
+        for( const auto up1 : geode::Indices{ colocated_info.unique_points } )
+        {
+            if( up1 <= up0 )
+            {
+                continue;
+            }
+            OPENGEODE_EXCEPTION(
+                geode::point_point_distance( colocated_info.unique_points[up0],
+                    colocated_info.unique_points[up1] )
+                    > DISTANCE,
+                "[Test] Colocated points are too close" );
+        }
+    }
+}
+
+void test()
+{
+    first_test();
+    second_test();
+}
+
 OPENGEODE_TEST( "nnsearch" )

tests/geometry/test-nnsearch.cpp

@@ -23,18 +23,23 @@
 
 #include <geode/tests/common.hpp>
 
+#include <async++.h>
+
 #include <geode/basic/assert.hpp>
 #include <geode/basic/logger.hpp>
 
+#include <geode/geometry/nn_search.hpp>
 #include <geode/geometry/point.hpp>
 
 #include <geode/mesh/builder/surface_mesh_builder.hpp>
 #include <geode/mesh/core/surface_mesh.hpp>
+#include <geode/mesh/core/triangulated_surface.hpp>
 #include <geode/mesh/helpers/convert_surface_mesh.hpp>
+#include <geode/mesh/io/triangulated_surface_input.hpp>
+#include <geode/mesh/io/triangulated_surface_output.hpp>
 
-void test()
+void test_create()
 {
-    geode::OpenGeodeMeshLibrary::initialize();
     std::vector< geode::Point2D > points{ geode::Point2D{ { 0, 0 } },
         geode::Point2D{ { 0, 1 } }, geode::Point2D{ { 0, 2 } },
         geode::Point2D{ { 1, 0 } }, geode::Point2D{ { 1, 1 } },
@@ -112,4 +117,62 @@ void test()
         "[Test] Wrong adjacency for { 3, 2 }" );
 }
 
+void test_import()
+{
+    auto surface = geode::load_triangulated_surface< 3 >(
+        absl::StrCat( geode::DATA_PATH, "surface_degen.og_tsf3d" ) );
+    auto bui = geode::TriangulatedSurfaceBuilder3D::create( *surface );
+
+    std::vector< std::reference_wrapper< const geode::SurfaceMesh3D > > meshes{
+        *surface
+    };
+    const auto merged = geode::merge_surface_meshes< 3 >( meshes );
+    geode::save_triangulated_surface(
+        *dynamic_cast< geode::TriangulatedSurface3D* >( merged.get() ),
+        "output.og_tsf3d" );
+
+    std::vector< geode::Point3D > points( merged->nb_vertices() );
+    for( const auto v : geode::Range{ merged->nb_vertices() } )
+    {
+        points[v] = merged->point( v );
+    }
+    geode::NNSearch3D nns( points );
+    const auto mappings = nns.colocated_index_mapping( geode::GLOBAL_EPSILON );
+    OPENGEODE_EXCEPTION( mappings.nb_colocated_points() == 0,
+        "[Test] Should be nomore colocated points" );
+    for( const auto p : geode::Indices{ points } )
+    {
+        OPENGEODE_EXCEPTION(
+            mappings.colocated_mapping[p] < mappings.unique_points.size(),
+            "[Test] Wrong value of colocated_mapping (bigger than unique "
+            "points size)" );
+        const auto& colocated_point =
+            mappings.unique_points[mappings.colocated_mapping[p]];
+        OPENGEODE_EXCEPTION( points[p].inexact_equal( colocated_point ),
+            "[Test] Colocated point is not close enough to original point" );
+    }
+    async::parallel_for(
+        async::irange( geode::index_t{ 0 }, mappings.unique_points.size() ),
+        [&mappings]( geode::index_t up0 ) {
+            for( const auto up1 : geode::Indices{ mappings.unique_points } )
+            {
+                if( up1 <= up0 )
+                {
+                    continue;
+                }
+                OPENGEODE_EXCEPTION( !mappings.unique_points[up0].inexact_equal(
+                                         mappings.unique_points[up1] ),
+                    "[Test] Colocated points are too close" );
+            }
+        } );
+}
+
+void test()
+{
+    geode::OpenGeodeMeshLibrary::initialize();
+    geode::Logger::set_level( geode::Logger::LEVEL::debug );
+    test_create();
+    test_import();
+}
+
 OPENGEODE_TEST( "merge-surface" )