Merge pull request #21677 from chacha21:rectangle_intersection

* better accuracy of _rotatedRectangleIntersection

instead of just migrating to double-precision (which would work), some computations are scaled by a factor that depends on the length of the smallest vectors.
There is a better accuracy even with floats, so this is certainly better for very sensitive cases

* Update intersection.cpp

use L2SQR norm to tune the numeric scale

* Update intersection.cpp

adapt samePointEps with L2 norm

* Update intersection.cpp

move comment

* Update intersection.cpp

fix wrong numericalScalingFactor usage

* added tests

* fixed warnings returned by buildbot

* modifications suggested by reviewer

renaming numericalScaleFctor to normalizationScale
refactor some computations
more "const"

* modifications as suggested by reviewer

modules/imgproc/src/intersection.cpp

@@ -51,15 +51,15 @@ static int _rotatedRectangleIntersection( const RotatedRect& rect1, const Rotate
 {
     CV_INSTRUMENT_REGION();
 
-    // L2 metric
-    const float samePointEps = std::max(1e-16f, 1e-6f * (float)std::max(rect1.size.area(), rect2.size.area()));
-
     Point2f vec1[4], vec2[4];
     Point2f pts1[4], pts2[4];
 
     rect1.points(pts1);
     rect2.points(pts2);
 
+    // L2 metric
+    float samePointEps = 1e-6f * (float)std::max(rect1.size.area(), rect2.size.area());
+
     int ret = INTERSECT_FULL;
 
     // Specical case of rect1 == rect2
@@ -99,14 +99,22 @@ static int _rotatedRectangleIntersection( const RotatedRect& rect1, const Rotate
         vec2[i].y = pts2[(i+1)%4].y - pts2[i].y;
     }
 
+    //we adapt the epsilon to the smallest dimension of the rects
+    for( int i = 0; i < 4; i++ )
+    {
+        samePointEps = std::min(samePointEps, std::sqrt(vec1[i].x*vec1[i].x+vec1[i].y*vec1[i].y));
+        samePointEps = std::min(samePointEps, std::sqrt(vec2[i].x*vec2[i].x+vec2[i].y*vec2[i].y));
+    }
+    samePointEps = std::max(1e-16f, samePointEps);
+
     // Line test - test all line combos for intersection
     for( int i = 0; i < 4; i++ )
     {
         for( int j = 0; j < 4; j++ )
         {
             // Solve for 2x2 Ax=b
-            float x21 = pts2[j].x - pts1[i].x;
+            const float x21 = pts2[j].x - pts1[i].x;
-            float y21 = pts2[j].y - pts1[i].y;
+            const float y21 = pts2[j].y - pts1[i].y;
 
             float vx1 = vec1[i].x;
             float vy1 = vec1[i].y;
@@ -114,10 +122,22 @@ static int _rotatedRectangleIntersection( const RotatedRect& rect1, const Rotate
             float vx2 = vec2[j].x;
             float vy2 = vec2[j].y;
 
-            float det = vx2*vy1 - vx1*vy2;
+            float normalizationScale  = std::min(vx1*vx1+vy1*vy1, vx2*vx2+vy2*vy2);//sum of squares : this is >= 0
+            //normalizationScale is a square, and we usually limit accuracy around 1e-6, so normalizationScale should be rather limited by ((1e-6)^2)=1e-12
+            normalizationScale  = (normalizationScale < 1e-12f) ? 1.f : 1.f/normalizationScale;
+
+            vx1 *= normalizationScale;
+            vy1 *= normalizationScale;
+            vx2 *= normalizationScale;
+            vy2 *= normalizationScale;
+
+            const float det = vx2*vy1 - vx1*vy2;
+            if (std::abs(det) < 1e-12)//like normalizationScale, we consider accuracy around 1e-6, i.e. 1e-12 when squared
+              continue;
+            const float detInvScaled = normalizationScale/det;
 
-            float t1 = (vx2*y21 - vy2*x21) / det;
+            const float t1 = (vx2*y21 - vy2*x21)*detInvScaled;
-            float t2 = (vx1*y21 - vy1*x21) / det;
+            const float t2 = (vx1*y21 - vy1*x21)*detInvScaled;
 
             // This takes care of parallel lines
             if( cvIsInf(t1) || cvIsInf(t2) || cvIsNaN(t1) || cvIsNaN(t2) )
@@ -127,8 +147,8 @@ static int _rotatedRectangleIntersection( const RotatedRect& rect1, const Rotate
 
             if( t1 >= 0.0f && t1 <= 1.0f && t2 >= 0.0f && t2 <= 1.0f )
             {
-                float xi = pts1[i].x + vec1[i].x*t1;
+                const float xi = pts1[i].x + vec1[i].x*t1;
-                float yi = pts1[i].y + vec1[i].y*t1;
+                const float yi = pts1[i].y + vec1[i].y*t1;
 
                 intersection.push_back(Point2f(xi,yi));
             }
@@ -149,18 +169,20 @@ static int _rotatedRectangleIntersection( const RotatedRect& rect1, const Rotate
         int posSign = 0;
         int negSign = 0;
 
-        float x = pts1[i].x;
+        const float x = pts1[i].x;
-        float y = pts1[i].y;
+        const float y = pts1[i].y;
 
         for( int j = 0; j < 4; j++ )
         {
+            float normalizationScale  = vec2[j].x*vec2[j].x+vec2[j].y*vec2[j].y;
+            normalizationScale  = (normalizationScale < 1e-12f) ? 1.f : 1.f/normalizationScale;
             // line equation: Ax + By + C = 0
             // see which side of the line this point is at
-            float A = -vec2[j].y;
+            const float A = -vec2[j].y*normalizationScale ;
-            float B = vec2[j].x;
+            const float B = vec2[j].x*normalizationScale ;
-            float C = -(A*pts2[j].x + B*pts2[j].y);
+            const float C = -(A*pts2[j].x + B*pts2[j].y);
 
-            float s = A*x+ B*y+ C;
+            const float s = A*x + B*y + C;
 
             if( s >= 0 )
             {
@@ -187,18 +209,22 @@ static int _rotatedRectangleIntersection( const RotatedRect& rect1, const Rotate
         int posSign = 0;
         int negSign = 0;
 
-        float x = pts2[i].x;
+        const float x = pts2[i].x;
-        float y = pts2[i].y;
+        const float y = pts2[i].y;
 
         for( int j = 0; j < 4; j++ )
         {
             // line equation: Ax + By + C = 0
             // see which side of the line this point is at
-            float A = -vec1[j].y;
+            float normalizationScale  = vec2[j].x*vec2[j].x+vec2[j].y*vec2[j].y;
-            float B = vec1[j].x;
+            normalizationScale  = (normalizationScale < 1e-12f) ? 1.f : 1.f/normalizationScale;
-            float C = -(A*pts1[j].x + B*pts1[j].y);
+            if (std::isinf(normalizationScale ))
+                normalizationScale  = 1.f;
+            const float A = -vec1[j].y*normalizationScale ;
+            const float B = vec1[j].x*normalizationScale ;
+            const float C = -(A*pts1[j].x + B*pts1[j].y);
 
-            float s = A*x + B*y + C;
+            const float s = A*x + B*y + C;
 
             if( s >= 0 )
             {
@@ -223,7 +249,7 @@ static int _rotatedRectangleIntersection( const RotatedRect& rect1, const Rotate
     }
 
     // Get rid of duplicated points
-    int Nstride = N;
+    const int Nstride = N;
     cv::AutoBuffer<float, 100> distPt(N * N);
     cv::AutoBuffer<int> ptDistRemap(N);
     for (int i = 0; i < N; ++i)
@@ -233,7 +259,7 @@ static int _rotatedRectangleIntersection( const RotatedRect& rect1, const Rotate
         for (int j = i + 1; j < N; )
         {
             const Point2f pt1 = intersection[j];
-            float d2 = normL2Sqr<float>(pt1 - pt0);
+            const float d2 = normL2Sqr<float>(pt1 - pt0);
             if(d2 <= samePointEps)
             {
                 if (j < N - 1)
@@ -252,10 +278,10 @@ static int _rotatedRectangleIntersection( const RotatedRect& rect1, const Rotate
         float minD = distPt[1];
         for (int i = 0; i < N - 1; ++i)
         {
-            float* pDist = distPt.data() + Nstride * ptDistRemap[i];
+            const float* pDist = distPt.data() + Nstride * ptDistRemap[i];
             for (int j = i + 1; j < N; ++j)
             {
-                float d = pDist[ptDistRemap[j]];
+                const float d = pDist[ptDistRemap[j]];
                 if (d < minD)
                 {
                     minD = d;

modules/imgproc/test/test_intersection.cpp

@@ -366,6 +366,84 @@ TEST(Imgproc_RotatedRectangleIntersection, regression_12221_2)
     EXPECT_LE(intersections.size(), (size_t)8);
 }
 
+TEST(Imgproc_RotatedRectangleIntersection, accuracy_21659)
+{
+    float scaleFactor = 1000;//to challenge the normalizationScale in the algorithm
+    cv::RectanglesIntersectTypes intersectionResult = cv::RectanglesIntersectTypes::INTERSECT_NONE;
+    std::vector<cv::Point2f> intersection;
+    double intersectionArea = 0;
+    cv::RotatedRect r1 = cv::RotatedRect(cv::Point2f(.5f, .5f)*scaleFactor, cv::Size2f(1.f, 1.f)*scaleFactor, 0);
+    cv::RotatedRect r2;
+
+    r2 = cv::RotatedRect(cv::Point2f(-2.f, -2.f)*scaleFactor, cv::Size2f(1.f, 1.f)*scaleFactor, 0);
+    intersectionResult = (cv::RectanglesIntersectTypes) cv::rotatedRectangleIntersection(r1, r2, intersection);
+    intersectionArea = (intersection.size() <= 2) ? 0. : cv::contourArea(intersection);
+    ASSERT_EQ(cv::RectanglesIntersectTypes::INTERSECT_NONE, intersectionResult);
+    ASSERT_LE(std::abs(intersectionArea-0), 1e-1);
+
+    r2 = cv::RotatedRect(cv::Point2f(1.5f, .5f)*scaleFactor, cv::Size2f(1.f, 2.f)*scaleFactor, 0);
+    intersectionResult = (cv::RectanglesIntersectTypes) cv::rotatedRectangleIntersection(r1, r2, intersection);
+    intersectionArea = (intersection.size() <= 2) ? 0. : cv::contourArea(intersection);
+    ASSERT_EQ(cv::RectanglesIntersectTypes::INTERSECT_PARTIAL, intersectionResult);
+    ASSERT_LE(std::abs(intersectionArea-0), 1e-1);
+
+    r2 = cv::RotatedRect(cv::Point2f(1.5f, 1.5f)*scaleFactor, cv::Size2f(1.f, 1.f)*scaleFactor, 0);
+    intersectionResult = (cv::RectanglesIntersectTypes) cv::rotatedRectangleIntersection(r1, r2, intersection);
+    intersectionArea = (intersection.size() <= 2) ? 0. : cv::contourArea(intersection);
+    ASSERT_EQ(cv::RectanglesIntersectTypes::INTERSECT_PARTIAL, intersectionResult);
+    ASSERT_LE(std::abs(intersectionArea-0), 1e-1);
+
+    r2 = cv::RotatedRect(cv::Point2f(.5f, .5f)*scaleFactor, cv::Size2f(1.f, 1.f)*scaleFactor, 0);
+    intersectionResult = (cv::RectanglesIntersectTypes) cv::rotatedRectangleIntersection(r1, r2, intersection);
+    intersectionArea = (intersection.size() <= 2) ? 0. : cv::contourArea(intersection);
+    ASSERT_EQ(cv::RectanglesIntersectTypes::INTERSECT_FULL, intersectionResult);
+    ASSERT_LE(std::abs(intersectionArea-r2.size.area()), 1e-1);
+
+    r2 = cv::RotatedRect(cv::Point2f(.5f, .5f)*scaleFactor, cv::Size2f(.5f, .5f)*scaleFactor, 0);
+    intersectionResult = (cv::RectanglesIntersectTypes) cv::rotatedRectangleIntersection(r1, r2, intersection);
+    intersectionArea = (intersection.size() <= 2) ? 0. : cv::contourArea(intersection);
+    ASSERT_EQ(cv::RectanglesIntersectTypes::INTERSECT_FULL, intersectionResult);
+    ASSERT_LE(std::abs(intersectionArea-r2.size.area()), 1e-1);
+
+    r2 = cv::RotatedRect(cv::Point2f(.5f, .5f)*scaleFactor, cv::Size2f(2.f, .5f)*scaleFactor, 0);
+    intersectionResult = (cv::RectanglesIntersectTypes) cv::rotatedRectangleIntersection(r1, r2, intersection);
+    intersectionArea = (intersection.size() <= 2) ? 0. : cv::contourArea(intersection);
+    ASSERT_EQ(cv::RectanglesIntersectTypes::INTERSECT_PARTIAL, intersectionResult);
+    ASSERT_LE(std::abs(intersectionArea-500000), 1e-1);
+
+    r2 = cv::RotatedRect(cv::Point2f(.5f, .5f)*scaleFactor, cv::Size2f(1.f, 1.f)*scaleFactor, 45);
+    intersectionResult = (cv::RectanglesIntersectTypes) cv::rotatedRectangleIntersection(r1, r2, intersection);
+    intersectionArea = (intersection.size() <= 2) ? 0. : cv::contourArea(intersection);
+    ASSERT_EQ(cv::RectanglesIntersectTypes::INTERSECT_PARTIAL, intersectionResult);
+    ASSERT_LE(std::abs(intersectionArea-828427), 1e-1);
+
+    r2 = cv::RotatedRect(cv::Point2f(1.f, 1.f)*scaleFactor, cv::Size2f(1.f, 1.f)*scaleFactor, 45);
+    intersectionResult = (cv::RectanglesIntersectTypes) cv::rotatedRectangleIntersection(r1, r2, intersection);
+    intersectionArea = (intersection.size() <= 2) ? 0. : cv::contourArea(intersection);
+    ASSERT_EQ(cv::RectanglesIntersectTypes::INTERSECT_PARTIAL, intersectionResult);
+    ASSERT_LE(std::abs(intersectionArea-250000), 1e-1);
+
+    //see #21659
+    r1 = cv::RotatedRect(cv::Point2f(4.48589373f, 12.5545063f), cv::Size2f(4.0f, 4.0f), 0.0347290039f);
+    r2 = cv::RotatedRect(cv::Point2f(4.48589373f, 12.5545235f), cv::Size2f(4.0f, 4.0f), 0.0347290039f);
+    intersectionResult = (cv::RectanglesIntersectTypes) cv::rotatedRectangleIntersection(r1, r2, intersection);
+    intersectionArea = (intersection.size() <= 2) ? 0. : cv::contourArea(intersection);
+    ASSERT_EQ(cv::RectanglesIntersectTypes::INTERSECT_PARTIAL, intersectionResult);
+    ASSERT_LE(std::abs(intersectionArea-r1.size.area()), 1e-3);
+
+    r1 = cv::RotatedRect(cv::Point2f(4.48589373f, 12.5545063f + 0.01f), cv::Size2f(4.0f, 4.0f), 0.0347290039f);
+    r2 = cv::RotatedRect(cv::Point2f(4.48589373f, 12.5545235f), cv::Size2f(4.0f, 4.0f), 0.0347290039f);
+    intersectionResult = (cv::RectanglesIntersectTypes) cv::rotatedRectangleIntersection(r1, r2, intersection);
+    intersectionArea = (intersection.size() <= 2) ? 0. : cv::contourArea(intersection);
+    ASSERT_LE(std::abs(intersectionArea-r1.size.area()), 1e-1);
+
+    r1 = cv::RotatedRect(cv::Point2f(45.0715866f, 39.8825722f), cv::Size2f(3.0f, 3.0f), 0.10067749f);
+    r2 = cv::RotatedRect(cv::Point2f(45.0715866f, 39.8825874f), cv::Size2f(3.0f, 3.0f), 0.10067749f);
+    intersectionResult = (cv::RectanglesIntersectTypes) cv::rotatedRectangleIntersection(r1, r2, intersection);
+    intersectionArea = (intersection.size() <= 2) ? 0. : cv::contourArea(intersection);
+    ASSERT_LE(std::abs(intersectionArea-r1.size.area()), 1e-3);
+}
+
 TEST(Imgproc_RotatedRectangleIntersection, regression_18520)
 {
     RotatedRect rr_empty(