fix(Triangulate): change algorithm for ear clipping

.vscode/settings.json

@@ -5,6 +5,7 @@
         "${workspaceFolder}/build/opengeode",
         "${workspaceFolder}/build/third_party/abseil/install/include",
         "${workspaceFolder}/build/third_party/asyncplusplus/install/include",
+        "${workspaceFolder}/build/third_party/earcut/install/include",
         "${workspaceFolder}/build/third_party/bitsery/install/include",
         "${workspaceFolder}/build/third_party/minizip/install/include",
         "${workspaceFolder}/build/third_party/nanoflann/install/include",

cmake/ConfigureOpenGeode.cmake

@@ -45,7 +45,7 @@ ExternalProject_Add(opengeode
         -DUSE_SUPERBUILD:BOOL=OFF
         -DASYNCPLUSPLUS_INSTALL_PREFIX:PATH=${ASYNCPLUSPLUS_INSTALL_PREFIX}
         -DBITSERY_INSTALL_PREFIX:PATH=${BITSERY_INSTALL_PREFIX}
-        -DFILESYSTEM_INSTALL_PREFIX:PATH=${FILESYSTEM_INSTALL_PREFIX}
+        -DEARCUT_INSTALL_PREFIX:PATH=${EARCUT_INSTALL_PREFIX}
         -DMINIZIP_INSTALL_PREFIX:PATH=${MINIZIP_INSTALL_PREFIX}
         -DNANOFLANN_INSTALL_PREFIX:PATH=${NANOFLANN_INSTALL_PREFIX}
         -DSPDLOG_INSTALL_PREFIX:PATH=${SPDLOG_INSTALL_PREFIX}
@@ -61,6 +61,7 @@ ExternalProject_Add(opengeode
         abseil
         asyncplusplus
         bitsery
+        earcut
         gdal
         minizip
         nanoflann
@@ -80,6 +81,7 @@ add_custom_target(download
         abseil-download
         asyncplusplus-download
         bitsery-download
+        earcut-download
         gdal-download
         minizip-download
         nanoflann-download

cmake/OpenGeode.cmake

@@ -31,6 +31,7 @@ include(cmake/PythonTargets.cmake)
 find_package(absl REQUIRED CONFIG NO_DEFAULT_PATH PATHS ${ABSEIL_INSTALL_PREFIX})
 find_package(Async++ REQUIRED CONFIG NO_DEFAULT_PATH PATHS ${ASYNCPLUSPLUS_INSTALL_PREFIX})
 find_package(Bitsery REQUIRED CONFIG NO_DEFAULT_PATH PATHS ${BITSERY_INSTALL_PREFIX})
+find_package(earcut_hpp REQUIRED CONFIG NO_DEFAULT_PATH PATHS ${EARCUT_INSTALL_PREFIX})
 find_package(minizip-ng REQUIRED CONFIG NO_DEFAULT_PATH PATHS ${MINIZIP_INSTALL_PREFIX})
 find_package(nanoflann REQUIRED CONFIG NO_DEFAULT_PATH PATHS ${NANOFLANN_INSTALL_PREFIX})
 find_package(spdlog REQUIRED CONFIG NO_DEFAULT_PATH PATHS ${SPDLOG_INSTALL_PREFIX})
@@ -52,6 +53,7 @@ install(
 if(NOT BUILD_SHARED_LIBS)
     install(
         DIRECTORY
+            ${EARCUT_INSTALL_PREFIX}/
             ${MINIZIP_INSTALL_PREFIX}/
             ${NANOFLANN_INSTALL_PREFIX}/
             ${SPDLOG_INSTALL_PREFIX}/

src/geode/mesh/CMakeLists.txt

@@ -330,5 +330,6 @@ add_geode_library(
         ${PROJECT_NAME}::geometry
     PRIVATE_DEPENDENCIES
         Async++
+        earcut_hpp::earcut_hpp
         ${PROJECT_NAME}::image
 )

src/geode/mesh/helpers/convert_surface_mesh.cpp

@@ -23,10 +23,13 @@
 
 #include <geode/mesh/helpers/convert_surface_mesh.hpp>
 
+#include <mapbox/earcut.hpp>
+
 #include <geode/basic/attribute_manager.hpp>
 #include <geode/basic/logger.hpp>
 
 #include <geode/geometry/point.hpp>
+#include <geode/geometry/vector.hpp>
 
 #include <geode/mesh/builder/polygonal_surface_builder.hpp>
 #include <geode/mesh/builder/triangulated_surface_builder.hpp>
@@ -39,6 +42,21 @@
 #include <geode/mesh/helpers/detail/surface_merger.hpp>
 #include <geode/mesh/helpers/internal/copy.hpp>
 
+namespace mapbox
+{
+    namespace util
+    {
+        template < std::size_t coord, geode::index_t dimension >
+        struct nth< coord, geode::Point< dimension > >
+        {
+            inline static auto get( const geode::Point< dimension >& point )
+            {
+                return point.value( coord );
+            };
+        };
+    } // namespace util
+} // namespace mapbox
+
 namespace
 {
     template < geode::index_t dimension >
@@ -273,35 +291,47 @@ namespace
     }
 
     template < geode::index_t dimension >
-    void transfer_adjacents( geode::SurfaceMeshBuilder< dimension >& builder,
-        absl::Span< const std::optional< geode::PolygonEdge > > adjacents,
-        absl::Span< const geode::index_t > new_polygons )
+    std::array< absl::FixedArray< geode::Point2D >, 1 > polygon_points(
+        const geode::SurfaceMesh< dimension >& surface,
+        geode::index_t polygon_id,
+        absl::Span< const geode::index_t > vertices );
+
+    template <>
+    std::array< absl::FixedArray< geode::Point2D >, 1 > polygon_points(
+        const geode::SurfaceMesh< 2 >& surface,
+        geode::index_t /*unused*/,
+        absl::Span< const geode::index_t > vertices )
     {
-        if( adjacents.front() )
-        {
-            builder.set_polygon_adjacent(
-                { new_polygons.front(), 0 }, adjacents.front()->polygon_id );
-            builder.set_polygon_adjacent(
-                adjacents.front().value(), new_polygons.front() );
-        }
-        for( const auto v : geode::LRange{ 1, adjacents.size() - 1 } )
+        std::array< absl::FixedArray< geode::Point2D >, 1 > polygons{
+            absl::FixedArray< geode::Point2D >( vertices.size() )
+        };
+        auto& polygon = polygons[0];
+        for( const auto v : geode::LIndices{ vertices } )
         {
-            if( adjacents[v] )
-            {
-                builder.set_polygon_adjacent(
-                    { new_polygons[v - 1], 1 }, adjacents[v]->polygon_id );
-                builder.set_polygon_adjacent(
-                    adjacents[v].value(), new_polygons[v - 1] );
-            }
+            polygon[v] = surface.point( vertices[v] );
         }
-        if( adjacents.back() )
+        return polygons;
+    }
+
+    template <>
+    std::array< absl::FixedArray< geode::Point2D >, 1 > polygon_points(
+        const geode::SurfaceMesh< 3 >& surface,
+        geode::index_t polygon_id,
+        absl::Span< const geode::index_t > vertices )
+    {
+        std::array< absl::FixedArray< geode::Point2D >, 1 > polygons{
+            absl::FixedArray< geode::Point2D >( vertices.size() )
+        };
+        auto& polygon = polygons[0];
+        const auto normal = surface.polygon_normal( polygon_id )
+                                .value_or( geode::Vector3D{ { 0, 0, 1 } } );
+        const auto axis_to_remove = normal.most_meaningful_axis();
+        for( const auto v : geode::LIndices{ vertices } )
         {
-            builder.set_polygon_adjacent(
-                { new_polygons.back(), 2 }, adjacents.back()->polygon_id );
-            builder.set_polygon_adjacent(
-                adjacents.back().value(), new_polygons.back() );
+            polygon[v] =
+                surface.point( vertices[v] ).project_point( axis_to_remove );
         }
-        builder.compute_polygon_adjacencies( new_polygons );
+        return polygons;
     }
 } // namespace
 
@@ -384,20 +414,50 @@ namespace geode
             to_delete[p] = nb_vertices != 3;
             if( nb_vertices > 3 )
             {
-                absl::FixedArray< std::optional< PolygonEdge > > adjacents(
-                    nb_vertices, std::nullopt );
+                using Edge = std::array< index_t, 2 >;
+                absl::flat_hash_map< Edge, PolygonEdge > adjacents;
+                const auto vertices = surface.polygon_vertices( p );
                 for( const auto e : LRange{ nb_vertices } )
                 {
-                    adjacents[e] = surface.polygon_adjacent_edge( { p, e } );
+                    if( const auto adj =
+                            surface.polygon_adjacent_edge( { p, e } ) )
+                    {
+                        adjacents.emplace(
+                            Edge{ vertices[e],
+                                vertices[e + 1 == nb_vertices ? 0 : e + 1] },
+                            adj.value() );
+                    }
                 }
-                absl::FixedArray< index_t > new_polygons( nb_vertices - 2 );
-                const auto vertices = surface.polygon_vertices( p );
-                for( const auto v : LRange{ 2, nb_vertices } )
+                const auto polygons = ::polygon_points( surface, p, vertices );
+                const auto new_triangles =
+                    mapbox::earcut< index_t >( polygons );
+                absl::FixedArray< index_t > new_polygons(
+                    new_triangles.size() / 3 );
+                for( const auto trgl : LIndices{ new_polygons } )
                 {
-                    new_polygons[v - 2] = builder.create_polygon(
-                        { vertices[0], vertices[v - 1], vertices[v] } );
+                    const std::array triangle{
+                        vertices[new_triangles[3 * trgl]],
+                        vertices[new_triangles[3 * trgl + 1]],
+                        vertices[new_triangles[3 * trgl + 2]]
+                    };
+                    new_polygons[trgl] = builder.create_polygon( triangle );
+                    for( const auto e : LRange{ 3 } )
+                    {
+                        const auto vertex0 = triangle[e];
+                        const auto vertex1 = triangle[e == 2 ? 0 : e + 1];
+                        const auto adj_it =
+                            adjacents.find( { vertex0, vertex1 } );
+                        if( adj_it == adjacents.end() )
+                        {
+                            continue;
+                        }
+                        builder.set_polygon_adjacent(
+                            adj_it->second, new_polygons[trgl] );
+                        builder.set_polygon_adjacent( { new_polygons[trgl], e },
+                            adj_it->second.polygon_id );
+                    }
                 }
-                ::transfer_adjacents( builder, adjacents, new_polygons );
+                builder.compute_polygon_adjacencies( new_polygons );
             }
         }
         to_delete.resize( surface.nb_polygons(), false );

tests/mesh/test-convert-surface.cpp

@@ -29,15 +29,17 @@
 #include <geode/geometry/point.hpp>
 #include <geode/geometry/vector.hpp>
 
+#include <geode/mesh/builder/surface_mesh_builder.hpp>
 #include <geode/mesh/core/light_regular_grid.hpp>
+#include <geode/mesh/core/polygonal_surface.hpp>
 #include <geode/mesh/core/triangulated_surface.hpp>
+#include <geode/mesh/io/polygonal_surface_output.hpp>
 #include <geode/mesh/io/triangulated_surface_input.hpp>
 
 #include <geode/mesh/helpers/convert_surface_mesh.hpp>
 
-void test()
+void convert_surface_dimension()
 {
-    geode::OpenGeodeMeshLibrary::initialize();
     const auto surface2d = geode::load_triangulated_surface< 2 >(
         absl::StrCat( geode::DATA_PATH, "3patches.og_tsf2d" ) );
     const auto surface3d =
@@ -70,7 +72,10 @@ void test()
         surface2d->nb_polygons() == converted_surface2d->nb_polygons(),
         "[Test] Number of polygons in re-converted TriangulatedSurface2D "
         "is not correct" );
+}
 
+void convert_grid_to_surface()
+{
     geode::LightRegularGrid2D grid{ geode::Point2D{ { 1, 2 } }, { 5, 6 },
         { 1, 1 } };
     const std::array< geode::index_t, 4 > cells_to_densify{ 5, 11, 12, 13 };
@@ -89,4 +94,32 @@ void test()
         "correct." );
 }
 
+void triangulate_surface()
+{
+    auto surface = geode::PolygonalSurface3D::create();
+    auto builder = geode::SurfaceMeshBuilder3D::create( *surface );
+
+    for( const auto i : geode::Range{ 10 } )
+    {
+        builder->create_point(
+            geode::Point3D{ { i * 0.5, i * 0.5, i * 0.5 } } );
+    }
+    builder->create_point( geode::Point3D{ { 10, 0, 10 } } );
+    std::vector< geode::index_t > polygon( surface->nb_vertices() );
+    absl::c_iota( polygon, 0 );
+    builder->create_polygon( polygon );
+    geode::triangulate_surface_mesh( *surface, *builder );
+    geode::save_polygonal_surface( *surface, "triangulated_surface.og_psf3d" );
+    OPENGEODE_EXCEPTION( surface->nb_polygons() == 9,
+        "[Test] Number of polygons in TriangulatedSurface3D is not correct" );
+}
+
+void test()
+{
+    geode::OpenGeodeMeshLibrary::initialize();
+    convert_surface_dimension();
+    convert_grid_to_surface();
+    triangulate_surface();
+}
+
 OPENGEODE_TEST( "convert-surface" )