From dd214962a5dd98658436b5d31160d1d250dfc9ad Mon Sep 17 00:00:00 2001 From: Adrian Kretz Date: Sun, 24 May 2026 09:53:48 +0200 Subject: [PATCH] Merge pull request #29109 from akretz:fix-issue-23644 Better Durand-Kerner Initialization #29109 While investigating issue #23644, I have found [this paper](https://link.springer.com/article/10.1007/BF01935059) which presents a good initialization for the Durand-Kerner algorithm. Basically the idea is to put the initial points equidistantly on a circle on the complex plane. The radius of the circle is computed as image Note that the $a_i$ coefficients in that paper are reversed compared to OpenCV. That's where the `(n - i)` in the code comes from. I have implemented just the mean of the $u_i$'s for the sake of simplicity. That's already enough to make the algorithm converge in all cases I have tested. I have used this to test for convergence for many polynomials of order 2 and 4 and coefficients of different magnitudes: ```cpp TEST(Core_SolvePoly, large_test) { cv::Mat_ coefs3(1,3); cv::Mat_ coefs5(1,5); cv::Mat r; double prec; for (int c0 = -20; c0 <= 20; c0++) { coefs3.at(0) = c0; for (int c1 = -20; c1 <= 20; c1++) { coefs3.at(1) = c1; for (int c2 = -20; c2 <= 20; c2++) { coefs3.at(2) = c2; prec = cv::solvePoly(coefs3, r); EXPECT_LE(prec, 1e-6); } } } for (int c0 = -10; c0 <= 10; c0++) { coefs5.at(0) = c0; for (int c1 = -10; c1 <= 10; c1++) { coefs5.at(1) = c1; for (int c2 = -10; c2 <= 10; c2++) { coefs5.at(2) = c2; for (int c3 = -10; c3 <= 10; c3++) { coefs5.at(3) = c3; for (int c4 = -10; c4 <= 10; c4++) { coefs5.at(4) = c4; prec = cv::solvePoly(coefs5, r); EXPECT_LE(prec, 1e-2); } } } } } for (int i = -10; i < 10; i++) { coefs3.at(0) = pow(2, i); for (int j = -10; j < 10; j++) { coefs3.at(1) = pow(2, j); for (int k = -10; k < 10; k++) { coefs3.at(2) = pow(2, k); prec = cv::solvePoly(coefs3, r); EXPECT_LE(prec, 1e-6); } } } } ``` This test passes, but I have not committed it because it runs for a couple of seconds. This fixes #23644 and replaces #29055. I have checked #29055 and it does not pass the test above. It seems to be optimized to the precise polynomial of #23644. ### Pull Request Readiness Checklist See details at https://github.com/opencv/opencv/wiki/How_to_contribute#making-a-good-pull-request - [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/core/misc/java/test/CoreTest.java | 26 +++++-- modules/core/src/mathfuncs.cpp | 35 ++++++++- modules/core/test/test_math.cpp | 87 +++++++++++++++++++++++ 3 files changed, 140 insertions(+), 8 deletions(-) diff --git a/modules/core/misc/java/test/CoreTest.java b/modules/core/misc/java/test/CoreTest.java index 62d72b706f..8d6b731616 100644 --- a/modules/core/misc/java/test/CoreTest.java +++ b/modules/core/misc/java/test/CoreTest.java @@ -1820,12 +1820,18 @@ public class CoreTest extends OpenCVTestCase { assertGE(1e-6, Math.abs(Core.solvePoly(coeffs, roots))); - truth = new Mat(3, 1, CvType.CV_32FC2) { + List rootsReIm = new ArrayList(); + Core.split(roots, rootsReIm); + Core.sort(rootsReIm.get(0), dst, Core.SORT_EVERY_COLUMN); + + Mat truthRe = new Mat(3, 1, CvType.CV_32F) { { - put(0, 0, 1, 0, 2, 0, 3, 0); + put(0, 0, 1, 2, 3); } }; - assertMatEqual(truth, roots, EPS); + Mat truthIm = Mat.zeros(3, 1, CvType.CV_32F); + assertMatEqual(truthRe, dst, EPS); + assertMatEqual(truthIm, rootsReIm.get(1), EPS); } public void testSolvePolyMatMatInt() { @@ -1836,14 +1842,20 @@ public class CoreTest extends OpenCVTestCase { }; Mat roots = new Mat(); - assertEquals(10.198039027185569, Core.solvePoly(coeffs, roots, 1)); + assertGE(1e-6, Core.solvePoly(coeffs, roots, 10)); - truth = new Mat(3, 1, CvType.CV_32FC2) { + List rootsReIm = new ArrayList(); + Core.split(roots, rootsReIm); + Core.sort(rootsReIm.get(0), dst, Core.SORT_EVERY_COLUMN); + + Mat truthRe = new Mat(3, 1, CvType.CV_32F) { { - put(0, 0, 1, 0, -1, 2, -2, 12); + put(0, 0, 1, 2, 3); } }; - assertMatEqual(truth, roots, EPS); + Mat truthIm = Mat.zeros(3, 1, CvType.CV_32F); + assertMatEqual(truthRe, dst, EPS); + assertMatEqual(truthIm, rootsReIm.get(1), EPS); } public void testSort() { diff --git a/modules/core/src/mathfuncs.cpp b/modules/core/src/mathfuncs.cpp index abcfb54215..c3bd50e0d4 100644 --- a/modules/core/src/mathfuncs.cpp +++ b/modules/core/src/mathfuncs.cpp @@ -1748,7 +1748,40 @@ double cv::solvePoly( InputArray _coeffs0, OutputArray _roots0, int maxIters ) break; } - C p(1, 0), r(1, 1); + // Related issue: https://github.com/opencv/opencv/issues/23644, + // This the initialization scheme of "Initial approximations in Durand-Kerner's root finding method" by Guggenheimer. + // https://link.springer.com/article/10.1007/BF01935059 + // We put the initial points equidistantly on a circle on the complex plane. This code computes the circle radius as in the paper. + Mat absCoeffs(n + 1, 1, CV_64F); + for( i = 0; i <= n; i++ ) + absCoeffs.at(i) = abs(coeffs[i]); + + int nonZeroCoeffs = 0; + Mat u(n, 1, CV_64F, Scalar(0)), v(n, 1, CV_64F, Scalar(0)); + for( i = 0; i <= n; i++ ) + { + double coeff = absCoeffs.at(i); + if( coeff > DBL_EPSILON ) + { + if( i != n ) + u.at(i) = 2.0 * pow(coeff / absCoeffs.at(n), 1.0 / (n - i)); + if( i != 0 ) + v.at(i - 1) = 0.5 * pow(absCoeffs.at(0) / coeff, 1.0 / i); + nonZeroCoeffs++; + } + } + double scale = 1; + if( nonZeroCoeffs > 2 ) + { + Point maxU, minV; + minMaxLoc(u, nullptr, nullptr, nullptr, &maxU); + minMaxLoc(v, nullptr, nullptr, &minV); + u.at(maxU) = 0; + v.at(minV) = 0; + scale = (sum(u).val[0] + sum(v).val[0]) / (2 * nonZeroCoeffs - 2); + } + + C p(scale, 0), r(cos(CV_2PI / n), sin(CV_2PI / n)); for( i = 0; i < n; i++ ) { diff --git a/modules/core/test/test_math.cpp b/modules/core/test/test_math.cpp index 1c9aa81d32..70c244d351 100644 --- a/modules/core/test/test_math.cpp +++ b/modules/core/test/test_math.cpp @@ -2679,6 +2679,93 @@ TEST(Core_SolvePoly, regression_5599) } } +TEST(Core_SolvePoly, regression_23644) +{ + // x^2 - 2x - 3 = 0, roots: 3, -1 + cv::Mat coefs = (cv::Mat_(1,3) << -3, -2, 1 ); + cv::Mat r; + double prec; + prec = cv::solvePoly(coefs, r); + EXPECT_LE(prec, 1e-6); + EXPECT_EQ(2u, r.total()); + ASSERT_EQ(CV_32FC2, r.type()); + checkRoot(r, 3, 0); + checkRoot(r, -1, 0); +} + +TEST(Core_SolvePoly, degree_2_polynomials) +{ + cv::Mat_ coefs(1,3); + cv::Mat r; + double prec; + for (int c0 = -20; c0 <= 20; c0++) + { + coefs.at(0) = c0; + for (int c1 = -20; c1 <= 20; c1++) + { + coefs.at(1) = c1; + for (int c2 = -20; c2 <= 20; c2++) + { + coefs.at(2) = c2; + prec = cv::solvePoly(coefs, r); + EXPECT_LE(prec, 1e-6); + } + } + } +} + +TEST(Core_SolvePoly, degree_4_polynomials) +{ + applyTestTag(CV_TEST_TAG_VERYLONG); + + cv::Mat_ coefs(1,5); + cv::Mat r; + double prec; + for (int c0 = -10; c0 <= 10; c0++) + { + coefs.at(0) = c0; + for (int c1 = -10; c1 <= 10; c1++) + { + coefs.at(1) = c1; + for (int c2 = -10; c2 <= 10; c2++) + { + coefs.at(2) = c2; + for (int c3 = -10; c3 <= 10; c3++) + { + coefs.at(3) = c3; + for (int c4 = -10; c4 <= 10; c4++) + { + coefs.at(4) = c4; + prec = cv::solvePoly(coefs, r); + EXPECT_LE(prec, 1e-3); + } + } + } + } + } +} + +TEST(Core_SolvePoly, different_magnitudes_polynomials) +{ + cv::Mat_ coefs(1,3); + cv::Mat r; + double prec; + for (int i = -10; i < 10; i++) + { + coefs.at(0) = pow(2, i); + for (int j = -10; j < 10; j++) + { + coefs.at(1) = pow(2, j); + for (int k = -10; k < 10; k++) + { + coefs.at(2) = pow(2, k); + prec = cv::solvePoly(coefs, r); + EXPECT_LE(prec, 1e-6); + } + } + } +} + class Core_PhaseTest : public cvtest::BaseTest { int t;