Merge pull request #1184 from Geode-solutions/fix_removing_apex_refine_in_remesh

fix(Quality): adding function to compute triangle quality, based on angles.

Author: BotellaA

include/geode/geometry/basic_objects/triangle.hpp

@@ -68,6 +68,10 @@ namespace geode
             normal() const;
         template < index_t T = dimension >
         [[nodiscard]]
+        typename std::enable_if< T == 3, std::optional< Vector3D > >::type
+            strict_normal() const;
+        template < index_t T = dimension >
+        [[nodiscard]]
         typename std::enable_if< T == 3, std::optional< Plane > >::type
             plane() const;
         template < index_t T = dimension >
@@ -79,6 +83,10 @@ namespace geode
         typename std::enable_if< T == 3, std::optional< local_index_t > >::type
             pivot() const;
         template < index_t T = dimension >
+        [[nodiscard]] typename std::enable_if< T == 3,
+            std::optional< std::pair< local_index_t, Vector3D > > >::type
+            strict_pivot_and_normal() const;
+        template < index_t T = dimension >
         [[nodiscard]] typename std::enable_if< T == 3,
             std::optional< std::pair< local_index_t, Vector3D > > >::type
             pivot_and_normal() const;

include/geode/geometry/quality.hpp

@@ -28,6 +28,8 @@
 namespace geode
 {
     class Tetrahedron;
+    FORWARD_DECLARATION_DIMENSION_CLASS( Triangle );
+    ALIAS_2D_AND_3D( Triangle );
 } // namespace geode
 
 namespace geode
@@ -39,4 +41,8 @@ namespace geode
 
     [[nodiscard]] double opengeode_geometry_api
         tetrahedron_collapse_aspect_ratio( const Tetrahedron& tetra );
+
+    template < index_t dimension >
+    [[nodiscard]] double triangle_angle_based_quality(
+        const Triangle< dimension >& triangle );
 } // namespace geode

src/geode/geometry/basic_objects/triangle.cpp

@@ -128,6 +128,18 @@ namespace geode
         return std::nullopt;
     }
 
+    template < typename PointType, index_t dimension >
+    template < index_t T >
+    typename std::enable_if< T == 3, std::optional< Vector3D > >::type
+        GenericTriangle< PointType, dimension >::strict_normal() const
+    {
+        if( const auto result = strict_pivot_and_normal() )
+        {
+            return result->second;
+        }
+        return std::nullopt;
+    }
+
     template < typename PointType, index_t dimension >
     template < index_t T >
     typename std::enable_if< T == 3, std::optional< Plane > >::type
@@ -166,6 +178,27 @@ namespace geode
         return std::nullopt;
     }
 
+    template < typename PointType, index_t dimension >
+    template < index_t T >
+    typename std::enable_if< T == 3,
+        std::optional< std::pair< local_index_t, Vector3D > > >::type
+        GenericTriangle< PointType, dimension >::strict_pivot_and_normal() const
+    {
+        const auto result = simple_pivot_and_normal(
+            { vertices_[0], vertices_[1], vertices_[2] } );
+        if( !result )
+        {
+            return std::nullopt;
+        }
+        if( result->pivot != NO_LID )
+        {
+            return std::optional< std::pair< local_index_t, Vector3D > >{
+                std::make_pair( result->pivot, result->normal )
+            };
+        }
+        return std::nullopt;
+    }
+
     template < typename PointType, index_t dimension >
     template < index_t T >
     typename std::enable_if< T == 3,
@@ -373,6 +406,8 @@ namespace geode
 
     template opengeode_geometry_api std::optional< Vector3D >
         GenericTriangle< Point< 3 >, 3 >::normal< 3 >() const;
+    template opengeode_geometry_api std::optional< Vector3D >
+        GenericTriangle< Point< 3 >, 3 >::strict_normal< 3 >() const;
     template opengeode_geometry_api std::optional< Plane >
         GenericTriangle< Point< 3 >, 3 >::plane< 3 >() const;
     template opengeode_geometry_api std::optional< OwnerPlane >
@@ -382,9 +417,14 @@ namespace geode
     template opengeode_geometry_api
         std::optional< std::pair< local_index_t, Vector3D > >
         GenericTriangle< Point< 3 >, 3 >::pivot_and_normal< 3 >() const;
+    template opengeode_geometry_api
+        std::optional< std::pair< local_index_t, Vector3D > >
+        GenericTriangle< Point< 3 >, 3 >::strict_pivot_and_normal< 3 >() const;
 
     template opengeode_geometry_api std::optional< Vector3D >
         GenericTriangle< RefPoint< 3 >, 3 >::normal< 3 >() const;
+    template opengeode_geometry_api std::optional< Vector3D >
+        GenericTriangle< RefPoint< 3 >, 3 >::strict_normal< 3 >() const;
     template opengeode_geometry_api std::optional< Plane >
         GenericTriangle< RefPoint< 3 >, 3 >::plane< 3 >() const;
     template opengeode_geometry_api std::optional< OwnerPlane >
@@ -394,4 +434,8 @@ namespace geode
     template opengeode_geometry_api
         std::optional< std::pair< local_index_t, Vector3D > >
         GenericTriangle< RefPoint< 3 >, 3 >::pivot_and_normal< 3 >() const;
+    template opengeode_geometry_api
+        std::optional< std::pair< local_index_t, Vector3D > >
+        GenericTriangle< RefPoint< 3 >, 3 >::strict_pivot_and_normal< 3 >()
+            const;
 } // namespace geode

src/geode/geometry/quality.cpp

@@ -25,14 +25,35 @@
 
 #include <limits>
 
+#include <geode/basic/logger.hpp>
+
 #include <geode/geometry/basic_objects/tetrahedron.hpp>
 #include <geode/geometry/basic_objects/triangle.hpp>
 #include <geode/geometry/mensuration.hpp>
 #include <geode/geometry/square_matrix.hpp>
 #include <geode/geometry/vector.hpp>
 
+#include <geode/basic/detail/disable_debug_logger.hpp>
+
 namespace
 {
+    template < geode::index_t dimension >
+    double compute_angle( const geode::Point< dimension >& point,
+        const geode::Point< dimension >& point_prev,
+        const geode::Point< dimension >& point_next )
+    {
+        const auto prev =
+            geode::Vector< dimension >{ point, point_prev }.normalize();
+        const auto next =
+            geode::Vector< dimension >{ point, point_next }.normalize();
+        const auto dot = std::clamp( prev.dot( next ), -1.0, 1.0 );
+        const auto angle = std::acos( dot );
+        if( std::isnan( angle ) )
+        {
+            return 0;
+        }
+        return angle;
+    }
 
     std::array< geode::local_index_t, 3 > lu_decomposition(
         geode::SquareMatrix3D& matrix )
@@ -241,4 +262,45 @@ namespace geode
         }
         return ( longest_edge_length * heightinv ) / std::sqrt( 3. / 2. );
     }
+
+    template < index_t dimension >
+    double triangle_angle_based_quality( const Triangle< dimension >& triangle )
+    {
+        try
+        {
+            const auto& vertices = triangle.vertices();
+            std::array< double, 3 > sinus;
+            for( const auto v : LRange{ 3 } )
+            {
+                const auto point = vertices[v].get();
+                const auto point_prev = vertices[( v + 2 ) % 3].get();
+                const auto point_next = vertices[( v + 1 ) % 3].get();
+                const auto angle =
+                    compute_angle( point, point_prev, point_next );
+                DEBUG( angle );
+                sinus[v] = std::sin( angle );
+                DEBUG( sinus[v] );
+            }
+            const auto denominator = sinus[0] + sinus[1] + sinus[2];
+            if( denominator == 0 )
+            {
+                return 0;
+            }
+            const auto quality =
+                4 * sinus[0] * sinus[1] * sinus[2] / ( denominator );
+            DEBUG( quality );
+            return quality;
+        }
+        catch( ... )
+        {
+            return 0;
+        }
+    }
+
+    template opengeode_geometry_api double triangle_angle_based_quality< 2 >(
+        const Triangle2D& );
+    template opengeode_geometry_api double triangle_angle_based_quality< 3 >(
+        const Triangle3D& );
 } // namespace geode
+
+#include <geode/basic/detail/enable_debug_logger.hpp>