fix(Geometry): implements a new computation of aspect ratio for tetrahedron

Author: panquez

include/geode/geometry/quality.hpp

@@ -36,4 +36,7 @@ namespace geode
         const Tetrahedron& tetra );
     [[nodiscard]] double opengeode_geometry_api
         tetrahedron_volume_to_edge_ratio( const Tetrahedron& tetra );
+
+    [[nodiscard]] double opengeode_geometry_api
+        tetrahedron_collapse_aspect_ratio( const Tetrahedron& tetra );
 } // namespace geode

src/geode/geometry/quality.cpp

@@ -21,18 +21,115 @@
  *
  */
 
-// Implementation inspired from https://github.com/sandialabs/verdict
-
 #include <geode/geometry/quality.hpp>
 
 #include <limits>
 
 #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>
 
+namespace
+{
+
+    std::array< geode::local_index_t, 3 > lu_decomposition(
+        geode::SquareMatrix3D& matrix )
+    {
+        std::array< geode::local_index_t, 3 > indices{ 0, 0, 0 };
+        double inverse_biggest[3];
+        for( const auto i : geode::LRange{ 3 } )
+        {
+            double biggest{ 0.0 };
+            for( const auto j : geode::LRange{ 3 } )
+            {
+                const auto t = std::fabs( matrix.value( i, j ) );
+                if( biggest < t )
+                    biggest = t;
+            }
+            if( biggest == 0.0 )
+            {
+                return indices;
+            }
+            inverse_biggest[i] = 1.0 / biggest;
+            indices[i] = i;
+        }
+        geode::local_index_t pivot_id{ 0 };
+        for( const auto k : geode::LRange{ 2 } )
+        {
+            double biggest{ 0.0 };
+            for( const auto i : geode::LRange{ k, 3 } )
+            {
+                const double t = std::fabs( matrix.value( indices[i], k ) )
+                                 * inverse_biggest[indices[i]];
+                if( biggest < t )
+                {
+                    biggest = t;
+                    pivot_id = i;
+                }
+            }
+            if( biggest == 0.0 )
+            {
+                return indices;
+            }
+            if( pivot_id != k )
+            {
+                const int j = indices[k];
+                indices[k] = indices[pivot_id];
+                indices[pivot_id] = j;
+            }
+            const auto pivot = matrix.value( indices[k], k );
+            for( const auto i : geode::LRange{ k + 1, 3 } )
+            {
+                const double m = matrix.value( indices[i], k ) / pivot;
+                matrix.set_value( indices[i], k, m );
+                if( m != 0.0 )
+                {
+                    for( const auto j : geode::LRange{ k + 1, 3 } )
+                    {
+                        matrix.set_value( indices[i], j,
+                            matrix.value( indices[i], j )
+                                - m * matrix.value( indices[k], j ) );
+                    }
+                }
+            }
+        }
+        return indices;
+    }
+    std::array< double, 3 > lu_solving( const geode::SquareMatrix3D& lu_matrix,
+        const std::array< double, 3 >& rhs,
+        const std::array< geode::local_index_t, 3 >& indices )
+    {
+        std::array< double, 3 > result{ 0., 0., 0. };
+        for( const auto i : geode::LRange{ 3 } )
+        {
+            double s{ 0.0 };
+            for( const auto j : geode::LRange{ 0, i } )
+            {
+                s += lu_matrix.value( indices[i], j ) * result[j];
+            }
+            result[i] = rhs[indices[i]] - s;
+        }
+        for( const auto i :
+            geode::LReverseRange{ 3 } ) // for(  int i = 2; i >= 0; --i )
+        {
+            double s{ 0.0 };
+            for( const auto j : geode::LRange{ i + 1, 3 } )
+            {
+                s += lu_matrix.value( indices[i], j ) * result[j];
+            }
+            result[i] = ( result[i] - s ) / lu_matrix.value( indices[i], i );
+        }
+
+        return result;
+    }
+
+} // namespace
+
 namespace geode
 {
+    // Implementation inspired from https://github.com/sandialabs/verdict
     double tetrahedron_aspect_ratio( const Tetrahedron& tetra )
     {
         const auto& vertices = tetra.vertices();
@@ -87,4 +184,57 @@ namespace geode
         const auto l_rms = std::sqrt( sq_len / 6 );
         return 6 * std::sqrt( 2 ) * signed_volume / ( l_rms * l_rms * l_rms );
     }
+
+    double tetrahedron_collapse_aspect_ratio( const Tetrahedron& tetra )
+    {
+        if( geode::tetrahedron_volume( tetra ) < geode::GLOBAL_EPSILON )
+        {
+            return std::numeric_limits< double >::max();
+        }
+
+        const auto& vertices = tetra.vertices();
+        const Vector3D edge_ab{ vertices[0], vertices[1] };
+        const Vector3D edge_ac{ vertices[0], vertices[2] };
+        const Vector3D edge_ad{ vertices[0], vertices[3] };
+        const Vector3D edge_bc{ vertices[1], vertices[2] };
+        const Vector3D edge_bd{ vertices[1], vertices[3] };
+        const Vector3D edge_cd{ vertices[2], vertices[3] };
+        const auto edge_ab_l2 = edge_ab.length2();
+        const auto edge_bc_l2 = edge_bc.length2();
+        const auto edge_ac_l2 = edge_ac.length2();
+        const auto edge_ad_l2 = edge_ad.length2();
+        const auto edge_bd_l2 = edge_bd.length2();
+        const auto edge_cd_l2 = edge_cd.length2();
+        const auto longest_edge_length = std::sqrt( std::max( { edge_ab_l2,
+            edge_bc_l2, edge_ac_l2, edge_ad_l2, edge_bd_l2, edge_cd_l2 } ) );
+
+        geode::SquareMatrix3D matrix{ { edge_ad, edge_bd, edge_cd } };
+        const auto row_indices = lu_decomposition( matrix );
+
+        std::array< Vector3D, 4 > N;
+        for( const auto j : geode::LRange{ 3 } )
+        {
+            std::array< double, 3 > b{ 0., 0., 0. };
+            b[j] = 1.0; // Positive means the inside direction
+            N[j] = Vector3D{ lu_solving( matrix, b, row_indices ) };
+        }
+        N[3] = ( N[0] + N[1] + N[2] ) * -1.;
+
+        std::array< double, 4 >
+            H; // H[i] = inverse of the height of its corresponding face
+        for( const auto i : geode::LRange{ 4 } )
+        {
+            H[i] = N[i].length();
+        }
+
+        auto heightinv = H[0];
+        for( const auto i : geode::LRange{ 1, 4 } )
+        { // Get the biggest H[i] (corresponding to the smallest height)
+            if( H[i] > heightinv )
+            {
+                heightinv = H[i];
+            }
+        }
+        return ( longest_edge_length * heightinv ) / std::sqrt( 3. / 2. );
+    }
 } // namespace geode

tests/geometry/test-quality.cpp

@@ -25,7 +25,9 @@
 
 #include <limits>
 
+#include <geode/geometry/basic_objects/segment.hpp>
 #include <geode/geometry/basic_objects/tetrahedron.hpp>
+#include <geode/geometry/basic_objects/triangle.hpp>
 #include <geode/geometry/distance.hpp>
 #include <geode/geometry/point.hpp>
 #include <geode/geometry/quality.hpp>
@@ -43,8 +45,14 @@ void test_perfect_tetrahedron()
         bary + geode::Point3D{ { 0, 0, std::sqrt( 1 - length * length ) } };
     const geode::Tetrahedron tetra{ point0, point1, point2, point3 };
     const auto quality = geode::tetrahedron_aspect_ratio( tetra );
-    OPENGEODE_EXCEPTION( quality - 1 < geode::GLOBAL_EPSILON,
+    const auto quality2 = geode::tetrahedron_volume_to_edge_ratio( tetra );
+    const auto quality3 = geode::tetrahedron_collapse_aspect_ratio( tetra );
+    OPENGEODE_EXCEPTION( std::fabs( quality - 1 ) < geode::GLOBAL_EPSILON,
         "[Test] Wrong aspect ratio for perfect tetrahedron ", quality );
+    OPENGEODE_EXCEPTION( std::fabs( quality2 - 1 ) < geode::GLOBAL_EPSILON,
+        "[Test] Wrong aspect ratio 2 for perfect tetrahedron ", quality2 );
+    OPENGEODE_EXCEPTION( std::fabs( quality3 - 1 ) < geode::GLOBAL_EPSILON,
+        "[Test] Wrong aspect ratio 3 for perfect tetrahedron ", quality3 );
 }
 
 void test_regular_tetrahedron()
@@ -55,8 +63,15 @@ void test_regular_tetrahedron()
     const geode::Point3D point3{ { 0, 0, 1 } };
     const geode::Tetrahedron tetra{ point0, point1, point2, point3 };
     const auto quality = geode::tetrahedron_aspect_ratio( tetra );
-    OPENGEODE_EXCEPTION( quality - 1.36603 < geode::GLOBAL_EPSILON,
+    const auto quality2 = geode::tetrahedron_volume_to_edge_ratio( tetra );
+    const auto quality3 = geode::tetrahedron_collapse_aspect_ratio( tetra );
+    OPENGEODE_EXCEPTION(
+        std::fabs( quality - 1.366026 ) < geode::GLOBAL_EPSILON,
         "[Test] Wrong aspect ratio for regular tetrahedron ", quality );
+    OPENGEODE_EXCEPTION( std::fabs( quality2 - 0.7698 ) < geode::GLOBAL_EPSILON,
+        "[Test] Wrong aspect ratio 2 for regular tetrahedron ", quality2 );
+    OPENGEODE_EXCEPTION( std::fabs( quality3 - 2 ) < geode::GLOBAL_EPSILON,
+        "[Test] Wrong aspect ratio 3 for regular tetrahedron ", quality3 );
 }
 
 void test_sliver_tetrahedron()
@@ -67,15 +82,45 @@ void test_sliver_tetrahedron()
     const geode::Point3D point3{ { 1, 1, 0 } };
     const geode::Tetrahedron tetra{ point0, point1, point2, point3 };
     const auto quality = geode::tetrahedron_aspect_ratio( tetra );
+    const auto quality2 = geode::tetrahedron_volume_to_edge_ratio( tetra );
+    const auto quality3 = geode::tetrahedron_collapse_aspect_ratio( tetra );
     OPENGEODE_EXCEPTION( quality == std::numeric_limits< double >::max(),
-        "[Test] Wrong aspect ratio for regular tetrahedron ", quality );
+        "[Test] Wrong aspect ratio for sliver tetrahedron ", quality );
+    OPENGEODE_EXCEPTION( std::fabs( quality2 ) < geode::GLOBAL_EPSILON,
+        "[Test] Wrong aspect ratio 2 for sliver tetrahedron ", quality2 );
+    OPENGEODE_EXCEPTION( quality3 == std::numeric_limits< double >::max(),
+        "[Test] Wrong aspect ratio 3 for sliver tetrahedron ", quality3 );
+}
+
+void test_other_tetrahedron()
+{
+    const geode::Point3D point0{ { 0, 0, 0 } };
+    const geode::Point3D point1{ { 4, 0, 0 } };
+    const geode::Point3D point2{ { 2, 1, 0 } };
+    const geode::Point3D point3{ { 2, 1, 1 } };
+    const geode::Tetrahedron tetra{ point0, point1, point2, point3 };
+    const auto quality = geode::tetrahedron_aspect_ratio( tetra );
+    const auto quality2 = geode::tetrahedron_volume_to_edge_ratio( tetra );
+    const auto quality3 = geode::tetrahedron_collapse_aspect_ratio( tetra );
+
+    OPENGEODE_EXCEPTION(
+        std::fabs( quality - 2.884068 ) < geode::GLOBAL_EPSILON,
+        "[Test] Wrong aspect ratio for random tetrahedron ", quality );
+    OPENGEODE_EXCEPTION(
+        std::fabs( quality2 - 0.341354 ) < geode::GLOBAL_EPSILON,
+        "[Test] Wrong aspect ratio 2 for random tetrahedron ", quality2 );
+    OPENGEODE_EXCEPTION(
+        std::fabs( quality3 - 4.618802 ) < geode::GLOBAL_EPSILON,
+        "[Test] Wrong aspect ratio 3 for random tetrahedron ", quality3 );
 }
 
 void test()
 {
+    geode::Logger::set_level( geode::Logger::LEVEL::debug );
     test_perfect_tetrahedron();
     test_regular_tetrahedron();
     test_sliver_tetrahedron();
+    test_other_tetrahedron();
 }
 
 OPENGEODE_TEST( "quality" )