Planar triangulations use SoS (#4435)

* planar_triangulations_use_SoS

* add test

* Update

* add comment

* simplify test

* fix build

* update comment

---------

Co-authored-by: Fedor Chelnokov <fedor.chelnokov@meshinspector.com>

Author: Grantim

source/MRMesh/MR2DContoursTriangulation.cpp

@@ -139,7 +139,7 @@ class SweepLineQueue
 
     bool less_( VertId l, VertId r ) const
     {
-        return std::tuple( pts_[l].x, pts_[l].y, l ) < std::tuple( pts_[r].x, pts_[r].y, r );
+        return smaller( { .id = l,.pt = pts_[l].x }, { .id = r,.pt = pts_[r].x } );
     }
 
 // INITIALIZATION CLASS BLOCK

source/MRMesh/MRHighPrecision.h

@@ -11,9 +11,11 @@ namespace MR
 /// \{
 
 using HighPrecisionInt = boost::multiprecision::checked_int128_t;
+using HighHighPrecisionInt = boost::multiprecision::checked_int256_t;
 
 using Vector2hp = Vector2<HighPrecisionInt>;
 using Vector3hp = Vector3<HighPrecisionInt>;
+using Vector3hhp = Vector3<HighHighPrecisionInt>;
 using Vector4hp = Vector4<HighPrecisionInt>;
 using Matrix3hp = Matrix3<HighPrecisionInt>;
 using Matrix4hp = Matrix4<HighPrecisionInt>;

source/MRMesh/MRPrecisePredicates2.cpp

@@ -1,6 +1,7 @@
 #include "MRPrecisePredicates2.h"
 #include "MRHighPrecision.h"
 #include "MRGTest.h"
+#include "MRPrecisePredicates3.h"
 
 namespace MR
 {
@@ -68,6 +69,22 @@ bool ccw( const PreciseVertCoords2* vs )
     return odd != ccw( vs[order[0]].pt, vs[order[1]].pt, vs[order[2]].pt );
 }
 
+bool inCircle( const std::array<PreciseVertCoords2, 4>& vs )
+{
+    return inCircle( vs.data() );
+}
+
+bool inCircle( const PreciseVertCoords2* vs )
+{
+    PreciseVertCoordsll vs3d[4];
+    for ( int i = 0; i < 4; ++i )
+    {
+        vs3d[i].id = vs[i].id;
+        vs3d[i].pt = Vector3ll( Vector3ll::ValueType( vs[i].pt.x ) * vs[i].pt.x + Vector3ll::ValueType( vs[i].pt.y ) * vs[i].pt.y, vs[i].pt.x, vs[i].pt.y );
+    }
+    return ccw( vs ) == orient3d( vs3d ); // TODO: looks like orient3d is not "honest" enough for this predicate
+}
+
 SegmentSegmentIntersectResult doSegmentSegmentIntersect( const std::array<PreciseVertCoords2, 4> & vs )
 {
     SegmentSegmentIntersectResult res;
@@ -192,4 +209,51 @@ TEST( MRMesh, PrecisePredicates2more )
     EXPECT_TRUE(  ccw( { vs[2],vs[3],vs[0] } ) );
 }
 
+TEST( MRMesh, PrecisePredicates2InCircle )
+{
+    std::array<PreciseVertCoords2, 4> vs =
+    {
+        PreciseVertCoords2{ 3_v, Vector2i{ -1, 2 } },
+        PreciseVertCoords2{ 2_v, Vector2i( 0 , 0 ) },
+        PreciseVertCoords2{ 0_v, Vector2i{ 3, 10 } },
+        PreciseVertCoords2{ 1_v, Vector2i{ 0 , 0 } }
+    };
+    EXPECT_TRUE( ccw( { vs[0],vs[1],vs[2] } ) );
+
+    // These 3 proves that vs[3] is inside vs[0]vs[1]vs[2] triangle
+    EXPECT_TRUE( ccw( { vs[0],vs[1],vs[3] } ) );
+    EXPECT_TRUE( ccw( { vs[1],vs[2],vs[3] } ) );
+    EXPECT_TRUE( ccw( { vs[2],vs[0],vs[3] } ) );
+
+    // Check that vs[3] is inCircle
+    EXPECT_TRUE( inCircle( vs ) );
+}
+
+TEST( MRMesh, PrecisePredicates2InCircle2 )
+{
+    std::array<PreciseVertCoords2, 5> vs =
+    {
+        PreciseVertCoords2{ 0_v, Vector2i{ -106280744 , -1002263723 } },
+        PreciseVertCoords2{ 1_v, Vector2i( -187288916 , -172107608 ) },
+        PreciseVertCoords2{ 2_v, Vector2i{ -25334363 , -1063004405 } },
+        PreciseVertCoords2{ 3_v, Vector2i{ -15200618 , -10122159 } },
+        PreciseVertCoords2{ 4_v, Vector2i{ -106280744 , -1002263723 } }
+    };
+
+    // Prove that 0_v 2_v 4_v circle is in +Y half plane (4_v 2_v is horde in lower part)
+    EXPECT_FALSE( ccw( { vs[2],vs[4],vs[3] } ) ); // 3_v is to the right of 2-4 vec
+    
+    // looks like this should be false
+    EXPECT_TRUE( inCircle( { vs[4],vs[2],vs[0],vs[3] } ) ); // 3_v is in circle
+
+    // prove that 0_v is inside 142 triangle
+    EXPECT_TRUE( ccw( { vs[1],vs[4],vs[0] } ) );
+    EXPECT_TRUE( ccw( { vs[4],vs[2],vs[0] } ) );
+    EXPECT_TRUE( ccw( { vs[2],vs[1],vs[0] } ) );
+    // it means that 142 circle should be larger in +Y half plane and so 3_v should be inside it
+    // but it is not
+    // DISABLED
+    //EXPECT_TRUE( inCircle( { vs[1],vs[4],vs[2],vs[3] } ) );
+}
+
 } //namespace MR

source/MRMesh/MRPrecisePredicates2.h

@@ -11,6 +11,17 @@ namespace MR
 /// \ingroup MathGroup
 /// \{
 
+struct PreciseVertCoord
+{
+    VertId id;   ///< unique id of the vertex (in both contours)
+    int pt;      ///< coordinate
+};
+
+/// return true if l is smaller then r
+/// uses simulation-of-simplicity to avoid "number are the same"
+inline bool smaller( const PreciseVertCoord& l, const PreciseVertCoord& r )
+{ return l.pt < r.pt || ( l.pt == r.pt && r.id < l.id ); }
+
 /// return true if the smallest rotation from vector (a) to vector (b) is in counter-clock-wise direction;
 /// uses simulation-of-simplicity to avoid "vectors are collinear"
 MRMESH_API bool ccw( const Vector2i & a, const Vector2i & b );
@@ -30,6 +41,10 @@ struct PreciseVertCoords2
 MRMESH_API bool ccw( const std::array<PreciseVertCoords2, 3> & vs );
 MRMESH_API bool ccw( const PreciseVertCoords2* vs );
 
+/// return true if 4th point in array lays inside circumcircle of first 3 points based triangle
+MRMESH_API bool inCircle( const std::array<PreciseVertCoords2, 4>& vs );
+MRMESH_API bool inCircle( const PreciseVertCoords2* vs );
+
 struct SegmentSegmentIntersectResult
 {
     bool doIntersect = false;   ///< whether the segments intersect

source/MRMesh/MRPrecisePredicates3.cpp

@@ -13,39 +13,40 @@ constexpr double cRangeIntMax = 0.99 * std::numeric_limits<int>::max(); // 0.99
 namespace MR
 {
 
-bool orient3d( const Vector3i & a, const Vector3i & b, const Vector3i & c )
+template<typename T, typename hT, typename hhT>
+bool orient3dT( const Vector3<T> & a, const Vector3<T>& b, const Vector3<T>& c )
 {
-    auto vhp = mixed( Vector3hp{ a }, Vector3hp{ b }, Vector3hp{ c } );
+    auto vhp = mixed( Vector3<hhT>{ a }, Vector3<hhT>{ b }, Vector3<hhT>{ c } );
     if ( vhp ) return vhp > 0;
 
-    auto v = cross( Vector2ll{ b.x, b.y }, Vector2ll{ c.x, c.y } );
+    auto v = cross( Vector2<hT>{ b.x, b.y }, Vector2<hT>{ c.x, c.y } );
     if ( v ) return v > 0;
 
-    v = -cross( Vector2ll{ b.x, b.z }, Vector2ll{ c.x, c.z } );
+    v = -cross( Vector2<hT>{ b.x, b.z }, Vector2<hT>{ c.x, c.z } );
     if ( v ) return v > 0;
 
-    v = cross( Vector2ll{ b.y, b.z }, Vector2ll{ c.y, c.z } );
+    v = cross( Vector2<hT>{ b.y, b.z }, Vector2<hT>{ c.y, c.z } );
     if ( v ) return v > 0;
 
-    v = -cross( Vector2ll{ a.x, a.y }, Vector2ll{ c.x, c.y } );
+    v = -cross( Vector2<hT>{ a.x, a.y }, Vector2<hT>{ c.x, c.y } );
     if ( v ) return v > 0;
 
     if ( c.x ) return c.x > 0;
 
     if ( c.y ) return c.y < 0;
 
-    v = cross( Vector2ll{ a.x, a.z }, Vector2ll{ c.x, c.z } );
+    v = cross( Vector2<hT>{ a.x, a.z }, Vector2<hT>{ c.x, c.z } );
     if ( v ) return v > 0;
 
     if ( c.z ) return c.z > 0;
 
 #ifndef NDEBUG
-    v = -cross( Vector2ll{ a.y, a.z }, Vector2ll{ c.y, c.z } );
+    v = -cross( Vector2<hT>{ a.y, a.z }, Vector2<hT>{ c.y, c.z } );
     assert( v == 0 );
     if ( v ) return v > 0;
 #endif
 
-    v = cross( Vector2ll{ a.x, a.y }, Vector2ll{ b.x, b.y } );
+    v = cross( Vector2<hT>{ a.x, a.y }, Vector2<hT>{ b.x, b.y } );
     if ( v ) return v > 0;
 
     if ( b.x ) return b.x < 0;
@@ -57,15 +58,21 @@ bool orient3d( const Vector3i & a, const Vector3i & b, const Vector3i & c )
     return true;
 }
 
-bool orient3d( const std::array<PreciseVertCoords, 4> & vs )
+bool orient3d( const Vector3i& a, const Vector3i& b, const Vector3i& c )
 {
-    return orient3d( vs.data() );
+    return orient3dT<int, long long, HighPrecisionInt>( a, b, c );
 }
 
-bool orient3d( const PreciseVertCoords* vs )
+bool orient3d( const Vector3ll& a, const Vector3ll& b, const Vector3ll& c )
+{
+    return orient3dT<long long, HighPrecisionInt, HighHighPrecisionInt>( a, b, c );
+}
+
+template<typename T>
+bool orient3dT( const T* vs )
 {
     bool odd = false;
-    std::array<int, 4> order = {0, 1, 2, 3};
+    std::array<int, 4> order = { 0, 1, 2, 3 };
 
     for ( int i = 0; i < 3; ++i )
     {
@@ -83,6 +90,26 @@ bool orient3d( const PreciseVertCoords* vs )
     return odd != orient3d( vs[order[0]].pt, vs[order[1]].pt, vs[order[2]].pt, vs[order[3]].pt );
 }
 
+bool orient3d( const std::array<PreciseVertCoords, 4> & vs )
+{
+    return orient3dT( vs.data() );
+}
+
+bool orient3d( const PreciseVertCoords* vs )
+{
+    return orient3dT( vs );
+}
+
+bool orient3d( const std::array<PreciseVertCoordsll, 4>& vs )
+{
+    return orient3dT( vs.data() );
+}
+
+bool orient3d( const PreciseVertCoordsll* vs )
+{
+    return orient3dT( vs );
+}
+
 TriangleSegmentIntersectResult doTriangleSegmentIntersect( const std::array<PreciseVertCoords, 5> & vs )
 {
     TriangleSegmentIntersectResult res;

source/MRMesh/MRPrecisePredicates3.h

@@ -14,21 +14,32 @@ namespace MR
 /// returns true if the plane with orientated triangle ABC has 0 point at the left;
 /// uses simulation-of-simplicity to avoid "0 is exactly on plane"
 MRMESH_API bool orient3d( const Vector3i & a, const Vector3i & b, const Vector3i & c );
+MRMESH_API bool orient3d( const Vector3ll& a, const Vector3ll& b, const Vector3ll& c );
 
 /// returns true if the plane with orientated triangle ABC has D point at the left;
 /// uses simulation-of-simplicity to avoid "D is exactly on plane"
 inline bool orient3d( const Vector3i & a, const Vector3i & b, const Vector3i & c, const Vector3i & d )
     { return orient3d( a - d, b - d, c - d ); }
+inline bool orient3d( const Vector3ll & a, const Vector3ll & b, const Vector3ll & c, const Vector3ll & d )
+    { return orient3d( a - d, b - d, c - d ); }
 
 struct PreciseVertCoords
 {
     VertId id;   ///< unique id of the vertex (in both meshes)
     Vector3i pt; ///< integer coordinates of the vertex
 };
 
+struct PreciseVertCoordsll
+{
+    VertId id;   ///< unique id of the vertex (in both meshes)
+    Vector3ll pt; ///< integer coordinates of the vertex
+};
+
 /// first sorts the indices in ascending order, then calls the predicate for sorted points
 MRMESH_API bool orient3d( const std::array<PreciseVertCoords, 4> & vs );
 MRMESH_API bool orient3d( const PreciseVertCoords* vs );
+MRMESH_API bool orient3d( const std::array<PreciseVertCoordsll, 4>& vs );
+MRMESH_API bool orient3d( const PreciseVertCoordsll* vs );
 
 struct TriangleSegmentIntersectResult
 {