better accuracy of _rotatedRectangleIntersection by chacha21 · Pull Request #21677 · opencv/opencv

chacha21 · 2022-03-03T13:46:20Z

Proposal for #21659

Migrating to double precision is a solution, but we can gain some accuracy even with floats by scaling computations according to the actual numeric scale of the current data.

More tests are needed, but it would be a waste of time if we just move to double precision instead.

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

use L2SQR norm to tune the numeric scale

adapt samePointEps with L2 norm

move comment

fix wrong numericalScalingFactor usage

alalek

It makes sense to take a look on existed code coverage and add appropriate tests before refactoring (to avoid regressions).
Test case from the issue should be added too.

alalek · 2022-03-10T19:13:11Z

modules/imgproc/src/intersection.cpp

    rect2.points(pts2);

+    // L2 metric
+    float samePointEps = 1e-6f * (float)std::max(rect1.size.area(), rect2.size.area());


We don't really want to have zero or denormalized zero value here.

What do you mean ? this is the original code, I just moved it.

alalek · 2022-03-10T19:15:02Z

modules/imgproc/src/intersection.cpp

+            vx1 *= numericalScalingFactor;
+            vy1 *= numericalScalingFactor;
+            vx2 *= numericalScalingFactor;
+            vy2 *= numericalScalingFactor;
+
            float det = vx2*vy1 - vx1*vy2;

-            float t1 = (vx2*y21 - vy2*x21) / det;
-            float t2 = (vx1*y21 - vy1*x21) / det;
+            float t1 = (vx2*y21 - vy2*x21)*numericalScalingFactor/det;
+            float t2 = (vx1*y21 - vy1*x21)*numericalScalingFactor/det;


We already have normalization multiplier in the first section. Why do we apply it twice for t1/t2?

This current code is correct : if you look closely, you can see that det is of the form alpha*numericalScalingFactor^2-beta*numericalScalingFactor^2, while (vx2*y21 - vy2*x21) and (vx1*y21 - vy1*x21) are of the form alpha*numericalScalingFactor-beta*numericalScalingFactor
That's why det must be divided by numericalScalingFactor to compute t1 and t2 and get the same result as when not using any numericalScalingFactor.
We could just change the code by float det = (vx2*vy1 - vx1*vy2)/numericalScalingFactor;, but I wanted to keep the initial det formula, and I prefered delaying that division to the computation of t1 and t2.

Agreed.

However more clear code should multiply:

x21 *= numericalScalingFactor; y21 *= numericalScalingFactor;

This code has 2 multiply operations (but could be executed in parallel).

Alternative is to replace 2 mul + 2 div ops => 2 mul + 1 div

float det_inv = numericalScalingFactor/det; float t1 = (...)*det_inv; float t2 = (...)*det_inv;

Bailout check could be changed too:

-if( cvIsInf(t1) || cvIsInf(t2) || cvIsNaN(t1) || cvIsNaN(t2) ) +if (cvIsInf(dev_inv)) continue;

But need to analyze numerical stability of this code first.

BTW, this code section is a generic "Solve for 2x2 Ax=b". Perhaps we could extract its code as a separate function (somewhere into core module).

alalek · 2022-03-10T19:16:00Z

modules/imgproc/src/intersection.cpp

            // 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 A = -vec2[j].y*numericalScalingFactor;
+            float B = vec2[j].x*numericalScalingFactor;
            float C = -(A*pts2[j].x + B*pts2[j].y);


Changing A, B without C doesn't look correct.

C will naturally be scaled since it is of the form C=alpha*A+beta*B, where A and B are already scaled.

alalek

Thank you for update!
Regression test looks good (fails without implementation changes)

alalek · 2022-03-14T13:21:48Z

modules/core/test/test_rotatedrect.cpp


 TEST(Core_RotatedRect, three_point_constructor) { Core_RotatedRectConstructorTest test; test.safe_run(); }

+TEST(Core_RotatedRect, intersection)


This test should be moved to imgproc module, see test_intersection.cpp file

alalek · 2022-03-14T13:24:15Z

modules/core/test/test_rotatedrect.cpp

+    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(intersectionResult, cv::RectanglesIntersectTypes::INTERSECT_NONE);


_EQ

Signature for _EQ is ASSERT_EQ(expected_reference, actual_result);.
Correct order of parameters are required to emit valid error messages.

Reference: https://github.com/opencv/opencv/blob/4.0.0/modules/ts/include/opencv2/ts/ts_gtest.h#L8196-L8200

GTEST_API_ AssertionResult EqFailure(const char* expected_expression, const char* actual_expression, const std::string& expected_value, const std::string& actual_value, bool ignoring_case);

alalek · 2022-03-14T13:30:51Z

modules/imgproc/src/intersection.cpp

+            numericalScalingFactor = !numericalScalingFactor ? 1.f : 1.f/(numericalScalingFactor);
+            if (std::isinf(numericalScalingFactor))
+                numericalScalingFactor = 1.f;


Other OpenCV code uses this pattern without secondary check:

numericalScalingFactor = std::abs(numericalScalingFactor) < 1e-6 ? 1.f : 1.f/(numericalScalingFactor);

std::abs call may be excluded too as we have abs call above.

Could you please to rename variable numerical -> normalize? e.g. normalizationScale or similar

alalek · 2022-03-14T13:44:04Z

modules/imgproc/src/intersection.cpp

+            vx1 *= numericalScalingFactor;
+            vy1 *= numericalScalingFactor;
+            vx2 *= numericalScalingFactor;
+            vy2 *= numericalScalingFactor;
+
            float det = vx2*vy1 - vx1*vy2;

-            float t1 = (vx2*y21 - vy2*x21) / det;
-            float t2 = (vx1*y21 - vy1*x21) / det;
+            float t1 = (vx2*y21 - vy2*x21)*numericalScalingFactor/det;
+            float t2 = (vx1*y21 - vy1*x21)*numericalScalingFactor/det;


Agreed.

However more clear code should multiply:

x21 *= numericalScalingFactor; y21 *= numericalScalingFactor;

This code has 2 multiply operations (but could be executed in parallel).

Alternative is to replace 2 mul + 2 div ops => 2 mul + 1 div

float det_inv = numericalScalingFactor/det; float t1 = (...)*det_inv; float t2 = (...)*det_inv;

Bailout check could be changed too:

-if( cvIsInf(t1) || cvIsInf(t2) || cvIsNaN(t1) || cvIsNaN(t2) ) +if (cvIsInf(dev_inv)) continue;

But need to analyze numerical stability of this code first.

BTW, this code section is a generic "Solve for 2x2 Ax=b". Perhaps we could extract its code as a separate function (somewhere into core module).

alalek · 2022-03-14T13:47:51Z

modules/imgproc/src/intersection.cpp

-    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];


BTW, this code is not performance critical, so probably we could use workaround with double computations inside instead of float.

I think it would be better, but as a second pull request, since the current ~~numericalScaleFactor~~ normalizationScale is a robustness improvement anyway

alalek · 2022-03-14T13:50:34Z

modules/imgproc/src/intersection.cpp

            float vx2 = vec2[j].x;
            float vy2 = vec2[j].y;

+            float numericalScalingFactor = std::min(std::abs(vx1*vx1+vy1*vy1), std::abs(vx2*vx2+vy2*vy2));


Why is std::min?
We need to "eliminate" large numbers.

I think both min and max are valid and (more or less) equivalent (relatively to numerical stability) for this use case.
The idea is to scale the range to be close to 1, so either you do it relatively to the lower bound (with a min), or to the higher bound (with a max); we could even scale relatively to the (logarithmic) mean value

alalek · 2022-03-14T13:52:23Z

modules/imgproc/src/intersection.cpp

+            if (std::isinf(numericalScalingFactor))
+                numericalScalingFactor = 1.f;
+
+            vx1 *= numericalScalingFactor;


BTW, Multiplication besides of 2^n could bring more artifacts.

I think that it would not be worse since a large exponent would certainly already be a problem about the significand digits in the current computations

renaming numericalScaleFctor to normalizationScale refactor some computations more "const"

alalek · 2022-03-16T10:42:29Z

modules/imgproc/src/intersection.cpp

+            vy2 *= normalizationScale;
+
+            const float det = vx2*vy1 - vx1*vy2;
+            const float detInvScaled = !det ? std::numeric_limits<float>::quiet_NaN() : (normalizationScale/det);


!det

Please avoid using bit-exact zero comparison with floating point variables. It could not work as expected.

Use bailout pre-check:

if (std::abs(det) < 1e-6) continue;

alalek

LGTM! Thank you 👍

* 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

chacha21 added 5 commits March 3, 2022 14:43

Update intersection.cpp

4fbb403

use L2SQR norm to tune the numeric scale

Update intersection.cpp

a33f57f

adapt samePointEps with L2 norm

Update intersection.cpp

ca0c915

move comment

Update intersection.cpp

523d4f5

fix wrong numericalScalingFactor usage

alalek reviewed Mar 10, 2022

View reviewed changes

chacha21 added 2 commits March 11, 2022 09:13

added tests

727c7d3

fixed warnings returned by buildbot

a47eccf

asmorkalov removed the pr: needs test New functionality requires minimal tests set label Mar 14, 2022

asmorkalov requested a review from alalek March 14, 2022 11:52

alalek reviewed Mar 14, 2022

View reviewed changes

modifications suggested by reviewer

abd32be

renaming numericalScaleFctor to normalizationScale refactor some computations more "const"

alalek reviewed Mar 16, 2022

View reviewed changes

modifications as suggested by reviewer

ee6bb6c

alalek approved these changes Mar 16, 2022

View reviewed changes

alalek merged commit ef6f421 into opencv:4.x Mar 16, 2022

opencv-pushbot mentioned this pull request Apr 23, 2022

(5.x) Merge 4.x #21903

Merged

Kroppeb mentioned this pull request Apr 25, 2023

rotatedRectangleIntersection does not return all intersection points as expected #23546

Closed

4 tasks


		TEST(Core_RotatedRect, three_point_constructor) { Core_RotatedRectConstructorTest test; test.safe_run(); }

		TEST(Core_RotatedRect, intersection)

Uh oh!

Conversation

chacha21 commented Mar 3, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

alalek left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

chacha21 Mar 10, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

chacha21 Mar 10, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

alalek left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

chacha21 Mar 14, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

alalek left a comment

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

chacha21 commented Mar 3, 2022 •

edited

Loading

alalek left a comment •

edited

Loading

chacha21 Mar 10, 2022 •

edited

Loading

chacha21 Mar 10, 2022 •

edited

Loading

chacha21 Mar 14, 2022 •

edited

Loading