From 93d490213fd6520c2af9d2b23f0c99039d355760 Mon Sep 17 00:00:00 2001 From: Pierre Chatelier Date: Tue, 30 May 2023 16:46:39 +0200 Subject: [PATCH] Merge pull request #23690 from chacha21:rotatedRectangleIntersection_precision better accuracy for _rotatedRectangleIntersection() (proposal for #23546) #23690 _rotatedRectangleIntersection() can be (statically) customized to use double instead of float for better accuracy this is a proposal for experimentation around #23546 for better accuracy, _rotatedRectangleIntersection() could use double. It will still return cv::Point2f list for backward compatibility, but the inner computations are controlled by a typedef - [X] I agree to contribute to the project under Apache 2 License. - [X] To the best of my knowledge, the proposed patch is not based on a code under GPL or another license that is incompatible with OpenCV - [X] The PR is proposed to the proper branch - [X] There is a reference to the original bug report and related work - [x] There is accuracy test, performance test and test data in opencv_extra repository, if applicable Patch to opencv_extra has the same branch name. - [x] The feature is well documented and sample code can be built with the project CMake --- modules/imgproc/src/intersection.cpp | 61 +++++++++------------- modules/imgproc/test/test_intersection.cpp | 19 +++++++ 2 files changed, 45 insertions(+), 35 deletions(-) diff --git a/modules/imgproc/src/intersection.cpp b/modules/imgproc/src/intersection.cpp index b9659f666e..858aa3941c 100644 --- a/modules/imgproc/src/intersection.cpp +++ b/modules/imgproc/src/intersection.cpp @@ -47,6 +47,14 @@ namespace cv { +static inline bool _isOnPositiveSide(const Point2f& line_vec, const Point2f& line_pt, const Point2f& pt) +{ + //we are interested by the cross product between the line vector (line_vec) and the line-to-pt vector (pt-line_pt) + //the sign of the only non-null component of the result determining which side of the line 'pt' is on + //the "positive" side meaning depends on the context usage of the current function and how line_vec and line_pt were filled + return (line_vec.y*(line_pt.x-pt.x) >= line_vec.x*(line_pt.y-pt.y)); +} + static int _rotatedRectangleIntersection( const RotatedRect& rect1, const RotatedRect& rect2, std::vector &intersection ) { CV_INSTRUMENT_REGION(); @@ -122,19 +130,10 @@ static int _rotatedRectangleIntersection( const RotatedRect& rect1, const Rotate float vx2 = vec2[j].x; float vy2 = vec2[j].y; - 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 + if (std::abs(det) < 1e-12)//we consider accuracy around 1e-6, i.e. 1e-12 when squared continue; - const float detInvScaled = normalizationScale/det; + const float detInvScaled = 1.f/det; const float t1 = (vx2*y21 - vy2*x21)*detInvScaled; const float t2 = (vx1*y21 - vy1*x21)*detInvScaled; @@ -169,22 +168,19 @@ static int _rotatedRectangleIntersection( const RotatedRect& rect1, const Rotate int posSign = 0; int negSign = 0; - const float x = pts1[i].x; - const float y = pts1[i].y; + const Point2f& pt = pts1[i]; 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 - 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); + // line equation: Ax + By + C = 0 where + // A = -vec2[j].y ; B = vec2[j].x ; C = -(A * pts2[j].x + B * pts2[j].y) + // check which side of the line this point is at + // A*x + B*y + C <> 0 + // + computation reordered for better numerical stability - const float s = A*x + B*y + C; + const bool isPositive = _isOnPositiveSide(vec2[j], pts2[j], pt); - if( s >= 0 ) + if( isPositive ) { posSign++; } @@ -209,24 +205,19 @@ static int _rotatedRectangleIntersection( const RotatedRect& rect1, const Rotate int posSign = 0; int negSign = 0; - const float x = pts2[i].x; - const float y = pts2[i].y; + const Point2f& pt = pts2[i]; for( int j = 0; j < 4; j++ ) { - // line equation: Ax + By + C = 0 - // see which side of the line this point is at - 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); + // line equation: Ax + By + C = 0 where + // A = -vec1[j].y ; B = vec1[j].x ; C = -(A * pts1[j].x + B * pts1[j].y) + // check which side of the line this point is at + // A*x + B*y + C <> 0 + // + computation reordered for better numerical stability - const float s = A*x + B*y + C; + const bool isPositive = _isOnPositiveSide(vec1[j], pts1[j], pt); - if( s >= 0 ) + if( isPositive ) { posSign++; } diff --git a/modules/imgproc/test/test_intersection.cpp b/modules/imgproc/test/test_intersection.cpp index 9ba3bf8125..f2067ba565 100644 --- a/modules/imgproc/test/test_intersection.cpp +++ b/modules/imgproc/test/test_intersection.cpp @@ -486,4 +486,23 @@ TEST(Imgproc_RotatedRectangleIntersection, regression_19824) EXPECT_LE(intersections.size(), (size_t)7); } +TEST(Imgproc_RotatedRectangleIntersection, accuracy_23546) +{ + RotatedRect r1( + Point2f(824.6421183672817f, 280.28737007069833f), + Size2f(565.0f, 140.0f), + -177.80506896972656f); + RotatedRect r2( + Point2f(567.3310438828003f, 270.42527355719545f), + Size2f(275.0f, 50.0f), + 92.19493103027344f); + + std::vector intersection; + double intersectionArea = 0; + int interType = cv::rotatedRectangleIntersection(r1, r2, intersection); + intersectionArea = (intersection.size() <= 2) ? 0. : cv::contourArea(intersection); + EXPECT_EQ(INTERSECT_PARTIAL, interType); + EXPECT_EQ(intersection.size(), (size_t)4); + ASSERT_LE(std::abs(intersectionArea-7000), 1e-1); +} }} // namespace