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 <img width="607" height="178" alt="image" src="https://github.com/user-attachments/assets/ea31b002-c924-4b93-9334-3e59597c896b" /> 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_<float> coefs3(1,3); cv::Mat_<float> coefs5(1,5); cv::Mat r; double prec; for (int c0 = -20; c0 <= 20; c0++) { coefs3.at<float>(0) = c0; for (int c1 = -20; c1 <= 20; c1++) { coefs3.at<float>(1) = c1; for (int c2 = -20; c2 <= 20; c2++) { coefs3.at<float>(2) = c2; prec = cv::solvePoly(coefs3, r); EXPECT_LE(prec, 1e-6); } } } for (int c0 = -10; c0 <= 10; c0++) { coefs5.at<float>(0) = c0; for (int c1 = -10; c1 <= 10; c1++) { coefs5.at<float>(1) = c1; for (int c2 = -10; c2 <= 10; c2++) { coefs5.at<float>(2) = c2; for (int c3 = -10; c3 <= 10; c3++) { coefs5.at<float>(3) = c3; for (int c4 = -10; c4 <= 10; c4++) { coefs5.at<float>(4) = c4; prec = cv::solvePoly(coefs5, r); EXPECT_LE(prec, 1e-2); } } } } } for (int i = -10; i < 10; i++) { coefs3.at<float>(0) = pow(2, i); for (int j = -10; j < 10; j++) { coefs3.at<float>(1) = pow(2, j); for (int k = -10; k < 10; k++) { coefs3.at<float>(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
2026-07-21 19:33:03 +04:00 · 2026-05-24 09:53:48 +02:00
parent 656f3395a7
commit dd214962a5
3 changed files with 140 additions and 8 deletions
@@ -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<Mat> rootsReIm = new ArrayList<Mat>();
+        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<Mat> rootsReIm = new ArrayList<Mat>();
+        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() {
@@ -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<double>(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<double>(i);
+        if( coeff > DBL_EPSILON )
+        {
+            if( i != n )
+                u.at<double>(i) = 2.0 * pow(coeff / absCoeffs.at<double>(n), 1.0 / (n - i));
+            if( i != 0 )
+                v.at<double>(i - 1) = 0.5 * pow(absCoeffs.at<double>(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<double>(maxU) = 0;
+        v.at<double>(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++ )
    {
@@ -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_<float>(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<float>(r, 3, 0);
+    checkRoot<float>(r, -1, 0);
+}
+
+TEST(Core_SolvePoly, degree_2_polynomials)
+{
+    cv::Mat_<float> coefs(1,3);
+    cv::Mat r;
+    double prec;
+    for (int c0 = -20; c0 <= 20; c0++)
+    {
+        coefs.at<float>(0) = c0;
+        for (int c1 = -20; c1 <= 20; c1++)
+        {
+            coefs.at<float>(1) = c1;
+            for (int c2 = -20; c2 <= 20; c2++)
+            {
+                coefs.at<float>(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_<float> coefs(1,5);
+    cv::Mat r;
+    double prec;
+    for (int c0 = -10; c0 <= 10; c0++)
+    {
+        coefs.at<float>(0) = c0;
+        for (int c1 = -10; c1 <= 10; c1++)
+        {
+            coefs.at<float>(1) = c1;
+            for (int c2 = -10; c2 <= 10; c2++)
+            {
+                coefs.at<float>(2) = c2;
+                for (int c3 = -10; c3 <= 10; c3++)
+                {
+                    coefs.at<float>(3) = c3;
+                    for (int c4 = -10; c4 <= 10; c4++)
+                    {
+                        coefs.at<float>(4) = c4;
+                        prec = cv::solvePoly(coefs, r);
+                        EXPECT_LE(prec, 1e-3);
+                    }
+                }
+            }
+        }
+    }
+}
+
+TEST(Core_SolvePoly, different_magnitudes_polynomials)
+{
+    cv::Mat_<float> coefs(1,3);
+    cv::Mat r;
+    double prec;
+    for (int i = -10; i < 10; i++)
+    {
+        coefs.at<float>(0) = pow(2, i);
+        for (int j = -10; j < 10; j++)
+        {
+            coefs.at<float>(1) = pow(2, j);
+            for (int k = -10; k < 10; k++)
+            {
+                coefs.at<float>(2) = pow(2, k);
+                prec = cv::solvePoly(coefs, r);
+                EXPECT_LE(prec, 1e-6);
+            }
+        }
+    }
+}
+
 class Core_PhaseTest : public cvtest::BaseTest
 {
    int t;