diff --git a/modules/imgproc/src/intersection.cpp b/modules/imgproc/src/intersection.cpp index 47d3f3f457..b9659f666e 100644 --- a/modules/imgproc/src/intersection.cpp +++ b/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; - float y21 = pts2[j].y - pts1[i].y; + const float x21 = pts2[j].x - pts1[i].x; + 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; - float t2 = (vx1*y21 - vy1*x21) / det; + const float t1 = (vx2*y21 - vy2*x21)*detInvScaled; + 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; - float yi = pts1[i].y + vec1[i].y*t1; + const float xi = pts1[i].x + vec1[i].x*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; - float y = pts1[i].y; + const float x = pts1[i].x; + 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; - float B = vec2[j].x; - float C = -(A*pts2[j].x + B*pts2[j].y); + const float A = -vec2[j].y*normalizationScale ; + const float B = vec2[j].x*normalizationScale ; + 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; - float y = pts2[i].y; + const float x = pts2[i].x; + 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 B = vec1[j].x; - float C = -(A*pts1[j].x + B*pts1[j].y); + float normalizationScale = vec2[j].x*vec2[j].x+vec2[j].y*vec2[j].y; + normalizationScale = (normalizationScale < 1e-12f) ? 1.f : 1.f/normalizationScale; + 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 distPt(N * N); cv::AutoBuffer 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(pt1 - pt0); + const float d2 = normL2Sqr(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; diff --git a/modules/imgproc/test/test_intersection.cpp b/modules/imgproc/test/test_intersection.cpp index c455c439fc..9ba3bf8125 100644 --- a/modules/imgproc/test/test_intersection.cpp +++ b/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 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(