precise predicate segmentIntersectionOrder for 2D (#4790)

Fedr · web-flow · commit f2896730f288 · 2025-07-04T16:54:53.000+03:00
diff --git a/source/MRMesh/MRMesh.vcxproj b/source/MRMesh/MRMesh.vcxproj
@@ -237,6 +237,7 @@
     <ClInclude Include="MRSeparationPoint.h" />
     <ClInclude Include="MRSharedThreadSafeOwner.h" />
     <ClInclude Include="MRSolarRadiation.h" />
+    <ClInclude Include="MRSparsePolynomial.h" />
     <ClInclude Include="MRSphere.h" />
     <ClInclude Include="MRStacktrace.h" />
     <ClInclude Include="MRSurfaceLineOffset.h" />
diff --git a/source/MRMesh/MRMesh.vcxproj.filters b/source/MRMesh/MRMesh.vcxproj.filters
@@ -1337,6 +1337,9 @@
     <ClInclude Include="MREndMill.h">
       <Filter>Source Files\Gcode</Filter>
     </ClInclude>
+    <ClInclude Include="MRSparsePolynomial.h">
+      <Filter>Source Files\Math</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="MRParallelProgressReporter.cpp">
diff --git a/source/MRMesh/MRPrecisePredicates2.cpp b/source/MRMesh/MRPrecisePredicates2.cpp
@@ -1,10 +1,73 @@
 #include "MRPrecisePredicates2.h"
 #include "MRHighPrecision.h"
 #include "MRPrecisePredicates3.h"
+#include "MRSparsePolynomial.h"
 
 namespace MR
 {
 
+namespace
+{
+
+struct PointDegree
+{
+    Vector2i pt;
+    int d = 0; // degree of epsilon for pt.y
+};
+
+std::array<PointDegree, 6> getPointDegrees( const std::array<PreciseVertCoords2, 6> & vs )
+{
+    struct VertN
+    {
+        VertId v;
+        int n = 0;
+    };
+    std::array<VertN, 6> as;
+    for ( int i = 0; i < 6; ++i )
+        as[i] = { vs[i].id, i };
+    std::sort( begin( as ), end( as ), []( const auto & a, const auto & b ) { return a.v < b.v; } );
+
+    std::array<PointDegree, 6> res;
+    int d = 1;
+    for ( int i = 0; i < 6; ++i )
+    {
+        assert( i == 0 || as[i-1].v < as[i].v ); // no duplicate vertices are permitted
+        const auto n = as[i].n;
+        res[n] = { vs[n].pt, d };
+        d *= 9;
+    }
+    return res;
+}
+
+SparsePolynomial<FastInt128> ccwPoly(
+    const PointDegree & a,
+    const PointDegree & b,
+    const PointDegree & c,
+    int db ) // degree.x = degree.y * db
+{
+    using Poly = SparsePolynomial<FastInt128>;
+
+    const Poly xx( a.pt.x - c.pt.x, a.d * db, 1, c.d * db, -1 );
+    const Poly xy( a.pt.y - c.pt.y, a.d     , 1, c.d     , -1 );
+    const Poly yx( b.pt.x - c.pt.x, b.d * db, 1, c.d * db, -1 );
+    const Poly yy( b.pt.y - c.pt.y, b.d     , 1, c.d     , -1 );
+    auto det = xx * yy;
+    det -= xy * yx;
+    return det;
+}
+
+bool isPositive( const SparsePolynomial<FastInt128>& poly )
+{
+    const auto & mapDegToCf = poly.get();
+    if ( !mapDegToCf.empty() )
+        return mapDegToCf.begin()->second > 0;
+
+    assert (false);
+    return false;
+}
+
+} // anonymous namespace
+
 // see https://arxiv.org/pdf/math/9410209 Table 4-i:
 // a=(pi_i,1, pi_i,2)
 // b=(pi_j,1, pi_j,2)
@@ -206,6 +269,32 @@ SegmentSegmentIntersectResult doSegmentSegmentIntersect( const std::array<Precis
     return res;
 }
 
+bool segmentIntersectionOrder( const std::array<PreciseVertCoords2, 6> & vs )
+{
+    // s=01, sa=23, sb=45
+    assert( doSegmentSegmentIntersect( { vs[0], vs[1], vs[2], vs[3] } ) );
+    assert( doSegmentSegmentIntersect( { vs[0], vs[1], vs[4], vs[5] } ) );
+    const auto ds = getPointDegrees( vs );
+
+    // res = ( ccw(sa,s[0])*ccw(sb,ds[1])   -   ccw(sb,s[0])*ccw(sa,ds[1]) ) /
+    //       ( ccw(sa,s[0])-ccw(sa,ds[1]) ) * ( ccw(sb,s[0])-ccw(sb,ds[1]) )
+    const auto polySaOrg  = ccwPoly( ds[2], ds[3], ds[0], 3 );
+    const auto polySaDest = ccwPoly( ds[2], ds[3], ds[1], 3 );
+    assert( isPositive( polySaOrg ) != isPositive( polySaDest ) );
+
+    const auto polySbOrg  = ccwPoly( ds[4], ds[5], ds[0], 3 );
+    const auto polySbDest = ccwPoly( ds[4], ds[5], ds[1], 3 );
+    assert( isPositive( polySbOrg ) != isPositive( polySbDest ) );
+
+    auto nom = polySaOrg * polySbDest;
+    nom -= polySbOrg * polySaDest;
+
+    bool res = isPositive( nom );
+    if ( isPositive( polySaOrg ) != isPositive( polySbOrg ) ) // denominator is negative
+        res = !res;
+    return res;
+}
+
 Vector2i findSegmentSegmentIntersectionPrecise(
     const Vector2i& ai, const Vector2i& bi, const Vector2i& ci, const Vector2i& di )
 {
diff --git a/source/MRMesh/MRPrecisePredicates2.h b/source/MRMesh/MRPrecisePredicates2.h
@@ -76,14 +76,19 @@ struct CoordinateConverters2
     ConvertToFloatVector2 toFloat{};
 };
 
+/// given line segment s=01 and two other segments sa=23, sb=45 known to intersect it, finds the order of intersection using precise predicates:
+/// true:  s[0], s ^ sa, s ^ sb, s[1]
+/// false: s[0], s ^ sb, s ^ sa, s[1]
+[[nodiscard]] MRMESH_API bool segmentIntersectionOrder( const std::array<PreciseVertCoords2, 6> & vs );
+
 /// finds intersection precise, using high precision int inside
 /// this function input should have intersection
 [[nodiscard]] MRMESH_API Vector2i findSegmentSegmentIntersectionPrecise(
     const Vector2i& a, const Vector2i& b, const Vector2i& c, const Vector2i& d );
 
 /// finds intersection precise, using high precision int inside
 /// this function input should have intersection
-[[nodiscard]] MRMESH_API Vector2f findSegmentSegmentIntersectionPrecise( 
+[[nodiscard]] MRMESH_API Vector2f findSegmentSegmentIntersectionPrecise(
     const Vector2f& a, const Vector2f& b, const Vector2f& c, const Vector2f& d,
     CoordinateConverters2 converters );
 
diff --git a/source/MRMesh/MRSparsePolynomial.h b/source/MRMesh/MRSparsePolynomial.h
@@ -0,0 +1,140 @@
+#pragma once
+
+#include <map>
+
+namespace MR
+{
+
+template <typename C, typename D>
+class SparsePolynomial;
+template <typename C, typename D>
+SparsePolynomial<C,D> operator *( const SparsePolynomial<C,D>& a, const SparsePolynomial<C,D>& b );
+
+/// The class to store a polynomial with a large number of zero coefficient (only non-zeros are stored in std::map)
+/// \tparam C - type of coefficients
+/// \tparam D - type of degrees
+template <typename C, typename D = int>
+class SparsePolynomial
+{
+public:
+    /// constructs zero polynomial
+    SparsePolynomial() = default;
+
+    /// takes existing coefficients in ownership
+    SparsePolynomial( std::map<D, C> && );
+
+    /// constructs polynomial c0 + c1*x^d1
+    SparsePolynomial( C c0, D d1, C c1 );
+
+    /// constructs polynomial c0 + c1*x^d1 + c2*x^d2
+    SparsePolynomial( C c0, D d1, C c1, D d2, C c2 );
+
+    /// gets read-only access to all not-zero coefficients
+    const std::map<D, C> & get() const { return map_; }
+
+    SparsePolynomial& operator +=( const SparsePolynomial& b );
+    SparsePolynomial& operator -=( const SparsePolynomial& b );
+    friend SparsePolynomial operator *<>( const SparsePolynomial& a, const SparsePolynomial& b );
+
+private:
+    std::map<D, C> map_; // degree -> not-zero coefficient
+};
+
+template <typename C, typename D>
+SparsePolynomial<C,D>::SparsePolynomial( std::map<D, C> && m ) : map_( std::move( m ) )
+{
+#ifndef NDEBUG
+    for ( const auto & [deg, cf] : map_ )
+        assert( cf != 0 );
+#endif
+}
+
+template <typename C, typename D>
+SparsePolynomial<C,D>::SparsePolynomial( C c0, D d1, C c1 )
+{
+    assert( c1 != 0 );
+    assert( d1 != 0 );
+    if ( c0 != 0 )
+        map_[D(0)] = c0;
+    map_[d1] = c1;
+}
+
+template <typename C, typename D>
+SparsePolynomial<C,D>::SparsePolynomial( C c0, D d1, C c1, D d2, C c2 )
+{
+    assert( c1 != 0 );
+    assert( d1 != 0 );
+    assert( c2 != 0 );
+    assert( d2 != 0 );
+    assert( d1 != d2 );
+    if ( c0 != 0 )
+        map_[D(0)] = c0;
+    map_[d1] = c1;
+    map_[d2] = c2;
+}
+
+template <typename C, typename D>
+SparsePolynomial<C,D>& SparsePolynomial<C,D>::operator +=( const SparsePolynomial& b )
+{
+    for ( const auto & [degB, cfB] : b.map_ )
+    {
+        assert( cfB != 0 );
+        auto [it,inserted] = map_.insert( { degB, cfB } );
+        if ( !inserted )
+        {
+            const auto sum = it->second + cfB;
+            if ( sum != 0 )
+                it->second = sum;
+            else
+                map_.erase( it );
+        }
+    }
+    return * this;
+}
+
+template <typename C, typename D>
+SparsePolynomial<C,D>& SparsePolynomial<C,D>::operator -=( const SparsePolynomial& b )
+{
+    for ( const auto & [degB, cfB] : b.map_ )
+    {
+        assert( cfB != 0 );
+        auto [it,inserted] = map_.insert( { degB, -cfB } );
+        if ( !inserted )
+        {
+            const auto sum = it->second - cfB;
+            if ( sum != 0 )
+                it->second = sum;
+            else
+                map_.erase( it );
+        }
+    }
+    return * this;
+}
+
+template <typename C, typename D>
+SparsePolynomial<C,D> operator *( const SparsePolynomial<C,D>& a, const SparsePolynomial<C,D>& b )
+{
+    std::map<D,C> resMap;
+    for ( const auto & [degA, cfA] : a.map_ )
+    {
+        assert( cfA != 0 );
+        for ( const auto & [degB, cfB] : b.map_ )
+        {
+            assert( cfB != 0 );
+            const auto deg = degA + degB;
+            const auto cf = cfA * cfB;
+            auto [it,inserted] = resMap.insert( { deg, cf } );
+            if ( !inserted )
+            {
+                const auto sum = it->second + cf;
+                if ( sum != 0 )
+                    it->second = sum;
+                else
+                    resMap.erase( it );
+            }
+        }
+    }
+    return resMap;
+}
+
+} //namespace MR
diff --git a/source/MRTest/MRPrecisePredicatesTests.cpp b/source/MRTest/MRPrecisePredicatesTests.cpp
@@ -109,6 +109,10 @@ TEST( MRMesh, PrecisePredicates2More )
 
     // segments 03 and 12 intersect one with another
     EXPECT_TRUE( doSegmentSegmentIntersect( { vs[0], vs[3], vs[1], vs[2] } ).doIntersect );
+
+    // intersection of 45 and 03 is closer to 4 than intersection of 45 and 12
+    EXPECT_TRUE(  segmentIntersectionOrder( { vs[4], vs[5], vs[0], vs[3], vs[1], vs[2] } ) );
+    EXPECT_FALSE( segmentIntersectionOrder( { vs[5], vs[4], vs[0], vs[3], vs[1], vs[2] } ) );
 }
 
 TEST( MRMesh, PrecisePredicates2FullDegen )
@@ -192,6 +196,33 @@ TEST( MRMesh, PrecisePredicates2InCircle2 )
     EXPECT_FALSE( inCircle( { vs[1],vs[4],vs[2],vs[3] } ) );
 }
 
+TEST( MRMesh, PreciseSegmentIntersectionOrder2 )
+{
+    PreciseVertCoords2 vs[6] =
+    {
+        // s:
+        PreciseVertCoords2{ 0_v, Vector2i( 0, 0 ) },
+        PreciseVertCoords2{ 1_v, Vector2i( 3, 0 ) },
+        // sa:
+        PreciseVertCoords2{ 2_v, Vector2i( 1,-1 ) },
+        PreciseVertCoords2{ 3_v, Vector2i( 1, 1 ) },
+        // sb:
+        PreciseVertCoords2{ 5_v, Vector2i( 2,-1 ) },
+        PreciseVertCoords2{ 6_v, Vector2i( 2, 1 ) }
+    };
+
+    EXPECT_TRUE(  segmentIntersectionOrder( { vs[0], vs[1], vs[2], vs[3], vs[4], vs[5] } ) );
+    EXPECT_TRUE(  segmentIntersectionOrder( { vs[0], vs[1], vs[3], vs[2], vs[4], vs[5] } ) );
+    EXPECT_TRUE(  segmentIntersectionOrder( { vs[0], vs[1], vs[3], vs[2], vs[5], vs[4] } ) );
+    EXPECT_FALSE( segmentIntersectionOrder( { vs[1], vs[0], vs[3], vs[2], vs[5], vs[4] } ) );
+
+    // swapped sa and sb
+    EXPECT_FALSE( segmentIntersectionOrder( { vs[0], vs[1], vs[4], vs[5], vs[2], vs[3] } ) );
+    EXPECT_FALSE( segmentIntersectionOrder( { vs[0], vs[1], vs[4], vs[5], vs[3], vs[2] } ) );
+    EXPECT_FALSE( segmentIntersectionOrder( { vs[0], vs[1], vs[5], vs[4], vs[3], vs[2] } ) );
+    EXPECT_TRUE(  segmentIntersectionOrder( { vs[1], vs[0], vs[5], vs[4], vs[3], vs[2] } ) );
+}
+
 TEST( MRMesh, PrecisePredicates3 )
 {
     const std::array<PreciseVertCoords, 5> vs = 
