fix(Geometry): tier 2 robustness fixes for degenerate inputs

include/geode/geometry/basic_objects/tetrahedron.hpp

@@ -43,12 +43,12 @@
 namespace geode
 {
     template < typename PointType >
     class GenericTetrahedron
     {
     public:
         static constexpr auto dim = 3;
         static constexpr inline std::array< std::array< index_t, 3 >, 4 >
             tetrahedron_facet_vertex{ { { { 1, 3, 2 } }, { { 0, 2, 3 } },
                 { { 3, 1, 0 } }, { { 0, 1, 2 } } } };
 
         GenericTetrahedron( PointType point0,
@@ -66,6 +66,10 @@
         void set_point( local_index_t vertex, PointType point );
         [[nodiscard]] const std::array< PointType, 4 >& vertices() const;
         [[nodiscard]] BoundingBox3D bounding_box() const;
+        /*!
+         * A tetrahedron is degenerated if all its vertices lie in a slab
+         * (thick plane) of thickness GLOBAL_EPSILON.
+         */
         [[nodiscard]] bool is_degenerated() const;
         [[nodiscard]] std::string string() const;
 
@@ -73,7 +77,7 @@
         std::array< PointType, 4 > vertices_;
     };
 
     class opengeode_geometry_api OwnerTetrahedron
         : public GenericTetrahedron< Point3D >
     {
         using Base = GenericTetrahedron< Point3D >;
@@ -84,14 +88,14 @@
             Point3D point2,
             Point3D point3 ) noexcept;
 
         OwnerTetrahedron( const Tetrahedron& other ) noexcept;
         OwnerTetrahedron( const OwnerTetrahedron& other ) noexcept;
         OwnerTetrahedron& operator=( const OwnerTetrahedron& other ) noexcept;
         OwnerTetrahedron( OwnerTetrahedron&& other ) noexcept;
         OwnerTetrahedron& operator=( OwnerTetrahedron&& other ) noexcept;
     };
 
     class opengeode_geometry_api Tetrahedron
         : public GenericTetrahedron< RefPoint3D >
     {
         using Base = GenericTetrahedron< RefPoint3D >;
@@ -103,7 +107,7 @@
             const Point3D& point3 ) noexcept;
 
         Tetrahedron( const Tetrahedron& other ) noexcept;
         Tetrahedron( const OwnerTetrahedron& other ) noexcept;
         Tetrahedron& operator=( const Tetrahedron& other ) noexcept;
         Tetrahedron( Tetrahedron&& other ) noexcept;
         Tetrahedron& operator=( Tetrahedron&& other ) noexcept;

include/geode/geometry/basic_objects/triangle.hpp

@@ -56,7 +56,7 @@
     } // namespace internal
 
     template < typename PointType, index_t dimension >
     class GenericTriangle
     {
     public:
         static constexpr auto dim = dimension;
@@ -98,6 +98,11 @@
         [[nodiscard]] local_index_t longest_edge_index() const;
         [[nodiscard]] local_index_t smallest_edge_index() const;
         [[nodiscard]] double minimum_height() const;
+        /*!
+         * A triangle is degenerated if all its vertices lie in a strip
+         * (thick line) of thickness GLOBAL_EPSILON, i.e. if its minimum
+         * height is below GLOBAL_EPSILON.
+         */
         [[nodiscard]] bool is_degenerated() const;
         [[nodiscard]] std::string string() const;
 
@@ -108,7 +113,7 @@
     };
 
     template < index_t dimension >
     class OwnerTriangle
         : public GenericTriangle< Point< dimension >, dimension >
     {
         using Base = GenericTriangle< Point< dimension >, dimension >;
@@ -127,7 +132,7 @@
     };
     ALIAS_2D_AND_3D( OwnerTriangle );
     template < index_t dimension >
     class Triangle : public GenericTriangle< RefPoint< dimension >, dimension >
     {
         using Base = GenericTriangle< RefPoint< dimension >, dimension >;
 

include/geode/geometry/internal/predicates.hpp → include/geode/geometry/detail/predicates.hpp

@@ -561,32 +561,32 @@ namespace GEO
 {
     namespace PCK
     {
-        SIGN orient_2d( const geode::Point2D& p0,
+        opengeode_geometry_api SIGN orient_2d( const geode::Point2D& p0,
             const geode::Point2D& p1,
             const geode::Point2D& p2 );
 
-        SIGN orient_3d( const geode::Point3D& p0,
+        opengeode_geometry_api SIGN orient_3d( const geode::Point3D& p0,
             const geode::Point3D& p1,
             const geode::Point3D& p2,
             const geode::Point3D& p3 );
 
-        SIGN det_3d( const geode::Vector3D& p0,
+        opengeode_geometry_api SIGN det_3d( const geode::Vector3D& p0,
             const geode::Vector3D& p1,
             const geode::Vector3D& p2 );
 
-        bool aligned_3d( const geode::Point3D& p0,
+        opengeode_geometry_api bool aligned_3d( const geode::Point3D& p0,
             const geode::Point3D& p1,
             const geode::Point3D& p2 );
 
-        SIGN dot_3d( const geode::Point3D& p0,
+        opengeode_geometry_api SIGN dot_3d( const geode::Point3D& p0,
             const geode::Point3D& p1,
             const geode::Point3D& p2 );
 
-        SIGN dot_2d( const geode::Point2D& p0,
+        opengeode_geometry_api SIGN dot_2d( const geode::Point2D& p0,
             const geode::Point2D& p1,
             const geode::Point2D& p2 );
 
-        void initialize();
+        opengeode_geometry_api void initialize();
     } // namespace PCK
 } // namespace GEO
 

include/geode/geometry/distance.hpp

@@ -43,6 +43,7 @@ namespace geode
     class Plane;
     class Tetrahedron;
     class Circle;
+    class Cylinder;
     template < index_t dimension >
     using Ball = Sphere< dimension >;
     using Disk = Circle;
@@ -261,6 +262,19 @@ namespace geode
     [[nodiscard]] std::tuple< double, Point3D > opengeode_geometry_api
         point_disk_distance( const Point3D& point, const Disk& disk );
 
+    /*!
+     * Compute the smallest distance between a point and a cylinder surface,
+     * i.e. its wall and its two end disks
+     * @return a tuple containing:
+     * - the smallest distance.
+     * - the closest point on the cylinder surface.
+     * @details Result is always positive or null, even if the point is inside
+     * the cylinder.
+     */
+    [[nodiscard]] std::tuple< double, Point3D >
+        opengeode_geometry_api point_cylinder_distance(
+            const Point3D& point, const Cylinder& cylinder );
+
     /*!
      * Compute the smallest distance between a point and an ellipse
      * @return a tuple containing:

include/geode/geometry/internal/intersection_from_sides.hpp

@@ -23,7 +23,7 @@
 
 #pragma once
 
-#include <geode/geometry/internal/predicates.hpp>
+#include <geode/geometry/detail/predicates.hpp>
 
 namespace geode
 {

include/geode/geometry/internal/position_from_sides.hpp

@@ -23,7 +23,7 @@
 
 #pragma once
 
-#include <geode/geometry/internal/predicates.hpp>
+#include <geode/geometry/detail/predicates.hpp>
 
 namespace geode
 {

include/geode/geometry/is_point_inside.hpp

@@ -37,7 +37,8 @@ namespace geode
 {
     /*!
      * Find if point is inside a polygon.
-     * @param[in] point The point to rotate.
+     * Points on the polygon boundary are considered inside.
+     * @param[in] point The point to test.
      */
     [[nodiscard]] bool opengeode_geometry_api is_point_inside_polygon(
         const Point2D& point, const Polygon2D& polygon );

include/geode/geometry/nn_search.hpp

@@ -106,7 +106,7 @@ namespace geode
         /*!
          * Get the closest neighbor from the given point
          * @param[in] point The requested point
-         * @return the index of the closest point
+         * @return the index of the closest point, NO_ID if there is no point
          */
         [[nodiscard]] index_t closest_neighbor(
             const Point< dimension >& point ) const;

src/geode/geometry/CMakeLists.txt

@@ -54,7 +54,7 @@ add_geode_library(
         "perpendicular.cpp"
         "points_sort.cpp"
         "position.cpp"
-        "internal/predicates.cpp"
+        "detail/predicates.cpp"
         "point.cpp"
         "projection.cpp"
         "quality.cpp"
@@ -105,11 +105,12 @@ add_geode_library(
     ADVANCED_HEADERS
         "detail/aabb_impl.hpp"
         "detail/bitsery_archive.hpp"
+        "detail/predicates.hpp"
     INTERNAL_HEADERS
         "internal/intersection_from_sides.hpp"
         "internal/position_from_sides.hpp"
-        "internal/predicates.hpp"
     SKIP_UNITY
+        "internal/predicates.cpp"
         "position.cpp"
     PUBLIC_DEPENDENCIES
         Bitsery::bitsery
@@ -120,3 +121,13 @@ add_geode_library(
         earcut_hpp::earcut_hpp
         nanoflann::nanoflann
 )
+
+# Exact predicates rely on exact floating-point rounding of each operation:
+# forbid contraction of a*b+c into fused multiply-add instructions.
+if(CMAKE_CXX_COMPILER_ID MATCHES "GNU|Clang")
+    set_source_files_properties(
+        "${PROJECT_SOURCE_DIR}/src/geode/geometry/internal/predicates.cpp"
+        PROPERTIES
+            COMPILE_OPTIONS "-ffp-contract=off"
+    )
+endif()

src/geode/geometry/basic_objects/circle.cpp

@@ -60,11 +60,7 @@ namespace geode
     {
         const auto sin_angle = [this]( Vector3D direction ) {
             const auto dP = plane_.normal().dot( direction );
-            if( dP <= -1.0 )
-            {
-                return M_PI;
-            }
-            if( dP >= 1.0 )
+            if( std::fabs( dP ) >= 1.0 )
             {
                 return 0.;
             }

src/geode/geometry/basic_objects/polygon.cpp

@@ -304,17 +304,18 @@ namespace geode
             geode::Segment< dimension >{
                 vertices_[max_length_edge], vertices_[next] }
         };
-        index_t opposite_vertex{ 0 };
+        double max_distance{ 0. };
         for( const auto vertex : geode::Range{ nb_vertices } )
         {
             if( vertex == max_length_edge || vertex == next )
             {
                 continue;
             }
-            opposite_vertex = vertex;
+            const geode::Point< dimension >& point = vertices_[vertex];
+            max_distance = std::max(
+                max_distance, geode::point_line_distance( point, line ) );
         }
-        const geode::Point< dimension >& point = vertices_[opposite_vertex];
-        return geode::point_line_distance( point, line );
+        return max_distance;
     }
 
     template < typename PointType, index_t dimension >

src/geode/geometry/basic_objects/tetrahedron.cpp

@@ -27,6 +27,7 @@
 #include <geode/geometry/basic_objects/triangle.hpp>
 #include <geode/geometry/bounding_box.hpp>
 #include <geode/geometry/distance.hpp>
+#include <geode/geometry/vector.hpp>
 
 namespace geode
 {
@@ -85,18 +86,30 @@ namespace geode
     template < typename PointType >
     bool GenericTetrahedron< PointType >::is_degenerated() const
     {
-        const Point3D& point0 = vertices_.at( 0 );
-        const Point3D& point1 = vertices_.at( 1 );
-        const Point3D& point2 = vertices_.at( 2 );
-        const Triangle3D triangle{ point0, point1, point2 };
-        if( triangle.is_degenerated() )
-        {
-            return true;
-        }
-        const Point3D& point3 = vertices_.at( 3 );
-        return std::get< 0 >(
-                   point_plane_distance( point3, Plane{ triangle } ) )
-               <= GLOBAL_EPSILON;
+        // The tetrahedron is degenerated if all its vertices lie in a slab of
+        // thickness GLOBAL_EPSILON. The thinnest slab is either bounded by a
+        // facet plane (width = 6 * volume / facet double area) or by two
+        // parallel planes containing opposite edges (width = 6 * volume /
+        // |edge cross product|): the largest of these denominators gives the
+        // minimum width.
+        const Point3D& point0 = vertices_[0];
+        const Point3D& point1 = vertices_[1];
+        const Point3D& point2 = vertices_[2];
+        const Point3D& point3 = vertices_[3];
+        const Vector3D edge01{ point0, point1 };
+        const Vector3D edge02{ point0, point2 };
+        const Vector3D edge03{ point0, point3 };
+        const Vector3D edge12{ point1, point2 };
+        const Vector3D edge13{ point1, point3 };
+        const Vector3D edge23{ point2, point3 };
+        const auto six_volume = edge01.dot( edge02.cross( edge03 ) );
+        const auto max_cross2 = std::max( { edge01.cross( edge02 ).length2(),
+            edge01.cross( edge03 ).length2(), edge02.cross( edge03 ).length2(),
+            edge12.cross( edge13 ).length2(), edge01.cross( edge23 ).length2(),
+            edge02.cross( edge13 ).length2(),
+            edge03.cross( edge12 ).length2() } );
+        return six_volume * six_volume
+               <= GLOBAL_EPSILON * GLOBAL_EPSILON * max_cross2;
     }
 
     template < typename PointType >

src/geode/geometry/basic_objects/triangle.cpp

@@ -35,6 +35,8 @@
 #include <geode/geometry/basic_objects/segment.hpp>
 #include <geode/geometry/bounding_box.hpp>
 #include <geode/geometry/distance.hpp>
+#include <geode/geometry/perpendicular.hpp>
+#include <geode/geometry/vector.hpp>
 
 namespace
 {
@@ -73,6 +75,19 @@ namespace
         }
         return result;
     }
+
+    double squared_double_area(
+        const geode::Vector2D& edge0, const geode::Vector2D& edge1 )
+    {
+        const auto double_area = geode::dot_perpendicular( edge0, edge1 );
+        return double_area * double_area;
+    }
+
+    double squared_double_area(
+        const geode::Vector3D& edge0, const geode::Vector3D& edge1 )
+    {
+        return edge0.cross( edge1 ).length2();
+    }
 } // namespace
 
 namespace geode
@@ -295,16 +310,18 @@ namespace geode
     template < typename PointType, index_t dimension >
     bool GenericTriangle< PointType, dimension >::is_degenerated() const
     {
-        const Point< dimension >& point0 = vertices_.at( 0 );
-        const Point< dimension >& point1 = vertices_.at( 1 );
-        const Segment< dimension > edge{ point0, point1 };
-        if( edge.is_degenerated() )
-        {
-            return true;
-        }
-        const Point< dimension >& point2 = vertices_.at( 2 );
-        return point_line_distance( point2, InfiniteLine< dimension >{ edge } )
-               <= GLOBAL_EPSILON;
+        // The triangle is degenerated if its minimum height, i.e. twice its
+        // area divided by its longest edge length, is below GLOBAL_EPSILON.
+        const Point< dimension >& point0 = vertices_[0];
+        const Point< dimension >& point1 = vertices_[1];
+        const Point< dimension >& point2 = vertices_[2];
+        const Vector< dimension > edge01{ point0, point1 };
+        const Vector< dimension > edge02{ point0, point2 };
+        const Vector< dimension > edge12{ point1, point2 };
+        const auto max_length2 = std::max(
+            { edge01.length2(), edge02.length2(), edge12.length2() } );
+        return squared_double_area( edge01, edge02 )
+               <= GLOBAL_EPSILON * GLOBAL_EPSILON * max_length2;
     }
 
     template < typename PointType, index_t dimension >