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
core: fix solveCubic numerical instability via coefficient normalization (fixes#27748) #28117
Summary
This PR fixes numerical instability in `cv::solveCubic` when the leading coefficient `a` is non-zero but extremely small relative to other coefficients (Issue #27748).
It introduces a **normalization step** that scales all coefficients by their maximum magnitude before solving. This ensures robust detection of when the equation should degenerate to a quadratic solver, without breaking valid cubic equations that happen to have small coefficients (e.g., scaled by 1e-9).
The Problem (Issue #27748)
The previous implementation checked `if (a == 0)` to decide whether to use the cubic or quadratic formula.
- When `a` is extremely small (e.g., 1e-17) but not exactly zero, and other coefficients are normal (e.g., 5.0), the standard cubic formula suffers from catastrophic cancellation and overflow, producing incorrect roots (e.g., 1e14).
The Fix
1. Normalization: The solver now finds `max_coeff = max(|a|, |b|, |c|, |d|)` and scales all coefficients by `1.0 / max_coeff`.
2. Relative Threshold: It then checks `if (abs(a) < epsilon)` on the *normalized* coefficients.
Why this is better than previous attempts
In a previous attempt (PR #28057), a simple absolute check `abs(a) < epsilon` was proposed. That approach was rejected because it failed for scaled equations.
Fixes#27748
Improve solveCubic accuracy #27347
### Pull Request Readiness Checklist
Fix#27323
```
2e-13 * x^3 + x^2 - 2 * x + 1 = 0 -> x^3 + 5e12 * x^2 - 1e13 * x + 5e12 = 0
```
The problem that coefficients have quite big magnitudes and current calculations are subject to round-off error
```
Q = (a1 * a1 - 3 * a2) * (1./9)
R = (2 * a1 * a1 * a1 - 9 * a1 * a2 + 27 * a3) * (1./54)
Qcubed = Q * Q * Q = a1^6/729 - (a1^4 a2)/81 + (a1^2 a2^2)/27 - a2^3/27
R * R = R^2 = a1^6/729 - (a1^4 a2)/81 + (a1^2 a2^2)/36 + (a1^3 a3)/27 - (a1 a2 a3)/6 + a3^2/4
d = Qcubed - R * R
```
Let `a1`, `a2`, `a3` have quite big same magnitudes, then we see that `Qcubed` and `R * R` have same terms `a1^6/729` and `-(a1^4 a2)/81` (which will be reduced in `d`), but they level out the other terms (these terms have `6`th and `5`th degree and other terms - less or equal than `4`th degree).
So, if these terms will participate in the calculation, this will lead to a huge round-off error.
But if we expand the expression, then round-off error should be less
```
d = Qcubed - R * R = 1/108 (a1^2 a2^2 - 4 a2^3 - 4 a1^3 a3 + 18 a1 a2 a3 - 27 a3^2)
```
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
Add tests for solveCubic #27331
### Pull Request Readiness Checklist
Related to #27323
I found only randomized tests with number of roots always equal to `1` or `3`, `x^3 = 0` and some simple test for Java and Swift.
Obviously, they don't cover all cases (implementation has strong branching and number of roots can be equal to `-1`, `0` and `2` additionally).
So, I think it will be useful to try explicitly cover more cases (and implementation branches correspondingly)
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
Add a basic sanity test to verify the rounding functions
work as expected.
Likewise, extend the rounding performance test to cover the
additional float -> int fast math functions.
* rewrote Mat::convertTo() and convertScaleAbs() to wide universal intrinsics; added always-available and SIMD-optimized FP16<=>FP32 conversion
* fixed compile warnings
* fix some more compile errors
* slightly relaxed accuracy threshold for int->float conversion (since we now do it using single-precision arithmetics, not double-precision)
* fixed compile errors on iOS, Android and in the baseline C++ version (intrin_cpp.hpp)
* trying to fix ARM-neon builds
* trying to fix ARM-neon builds
* trying to fix ARM-neon builds
* trying to fix ARM-neon builds
* core:OE-27 prepare universal intrinsics to expand (#11022)
* core:OE-27 prepare universal intrinsics to expand (#11022)
* core: Add universal intrinsics for AVX2
* updated implementation of wide univ. intrinsics; converted several OpenCV HAL functions: sqrt, invsqrt, magnitude, phase, exp to the wide universal intrinsics.
* converted log to universal intrinsics; cleaned up the code a bit; added v_lut_deinterleave intrinsics.
* core: Add universal intrinsics for AVX2
* fixed multiple compile errors
* fixed many more compile errors and hopefully some test failures
* fixed some more compile errors
* temporarily disabled IPP to debug exp & log; hopefully fixed Doxygen complains
* fixed some more compile errors
* fixed v_store(short*, v_float16&) signatures
* trying to fix the test failures on Linux
* fixed some issues found by alalek
* restored IPP optimization after the patch with AVX wide intrinsics has been properly tested
* restored IPP optimization after the patch with AVX wide intrinsics has been properly tested
- 'if' logic is moved into templates.
- removed unnecessary cv::Mat objects creation.
- fixed inv() test (invA * A == eye)
- added more Matx tests to cover all defined template specializations
- removed tr1 usage (dropped in C++17)
- moved includes of vector/map/iostream/limits into ts.hpp
- require opencv_test + anonymous namespace (added compile check)
- fixed norm() usage (must be from cvtest::norm for checks) and other conflict functions
- added missing license headers
* add accuracy test and performance check for matmul
* add performance tests for transform and dotProduct
* add test Core_TransformLargeTest for 8u version of transform
* remove raw SSE2/NEON implementation from matmul.cpp
* use universal intrinsic instead of raw intrinsic
* remove unused templated function
* add v_matmuladd which multiply 3x3 matrix and add 3x1 vector
* add v_rotate_left/right in universal intrinsic
* suppress intrinsic on some function and platform
* add pure SW implementation of new universal intrinsics
* add test for new universal intrinsics
* core: prevent memory access after the end of buffer
* fix perf tests
Alexander Smorkalov
Adrian Kretz
Kumataro
Vincent Rabaud
satyam yadav
Maxim Smolskiy
Alexander Smorkalov
Maksim Shabunin
Pierre Chatelier
Alexander Alekhin
Rostislav Vasilikhin
Suleyman TURKMEN
yuki takehara
Paul E. Murphy
Vadim Pisarevsky
Alexander Alekhin
Vitaly Tuzov
Vitaly Tuzov
Rostislav Vasilikhin
Tomoaki Teshima
Vladislav Sovrasov
LaurentBerger
Alexander Alekhin
Maksim Shabunin
hoangviet1985
Adil Ibragimov
Roman Donchenko