mirror/opencv
1
0
Fork 0
mirror of https://github.com/opencv/opencv.git synced 2024-11-25 19:50:38 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity

squash

Browse Source
This commit is contained in:
hoangviet1985 2015-11-20 14:48:29 -05:00
parent 6e2a68fdf3
commit b96def885f
2 changed files with 91 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
modules/core/src/mathfuncs.cpp
View File

@ -2085,12 +2085,53 @@ double cv::solvePoly( InputArray _coeffs0, OutputArray _roots0, int maxIters )
{
p = roots[i];
C num = coeffs[n], denom = coeffs[n];
int num_same_root = 1;
for( j = 0; j < n; j++ )
{
num = num*p + coeffs[n-j-1];
if( j != i ) denom = denom * (p - roots[j]);
if( j != i )
{
if ( (p - roots[j]).re != 0 || (p - roots[j]).im != 0 )
denom = denom * (p - roots[j]);
else
num_same_root++;
}
}
num /= denom;
if( num_same_root > 1)
{
double old_num_re = num.re;
double old_num_im = num.im;
int square_root_times = num_same_root % 2 == 0 ? num_same_root / 2 : num_same_root / 2 - 1;
for( j = 0; j < square_root_times; j++)
{
num.re = old_num_re*old_num_re + old_num_im*old_num_im;
num.re = sqrt(num.re);
num.re += old_num_re;
num.im = num.re - old_num_re;
num.re /= 2;
num.re = sqrt(num.re);
num.im /= 2;
num.im = sqrt(num.im);
if( old_num_re < 0 ) num.im = -num.im;
}
if( num_same_root % 2 != 0){
Mat cube_coefs(4, 1, CV_64FC1);
Mat cube_roots(3, 1, CV_64FC2);
cube_coefs.at<double>(3) = -(pow(old_num_re, 3));
cube_coefs.at<double>(2) = -(15*pow(old_num_re, 2) + 27*pow(old_num_im, 2));
cube_coefs.at<double>(1) = -48*old_num_re;
cube_coefs.at<double>(0) = 64;
solveCubic(cube_coefs, cube_roots);
if(cube_roots.at<double>(0) >= 0) num.re = pow(cube_roots.at<double>(0), 1./3);
else num.re = -pow(-cube_roots.at<double>(0), 1./3);
num.im = sqrt(pow(num.re, 2) / 3 - old_num_re / (3*num.re));
}
}
roots[i] = p - num;
maxDiff = std::max(maxDiff, cv::abs(num));
}

49
modules/core/test/test_math.cpp
View File

@ -2349,6 +2349,55 @@ void Core_SolvePolyTest::run( int )
}
}
template<typename T>
static void checkRoot(Mat& r, T re, T im)
{
for (int i = 0; i < r.cols*r.rows; i++)
{
Vec<T, 2> v = *(Vec<T, 2>*)r.ptr(i);
if (fabs(re - v[0]) < 1e-6 && fabs(im - v[1]) < 1e-6)
{
v[0] = std::numeric_limits<T>::quiet_NaN();
v[1] = std::numeric_limits<T>::quiet_NaN();
return;
}
}
GTEST_NONFATAL_FAILURE_("Can't find root") << "(" << re << ", " << im << ")";
}
TEST(Core_SolvePoly, regression_5599)
{
// x^4 - x^2 = 0, roots: 1, -1, 0, 0
cv::Mat coefs = (cv::Mat_<float>(1,5) << 0, 0, -1, 0, 1 );
{
cv::Mat r;
double prec;
prec = cv::solvePoly(coefs, r);
EXPECT_LE(prec, 1e-6);
EXPECT_EQ(4, r.total());
//std::cout << "Preciseness = " << prec << std::endl;
//std::cout << "roots:\n" << r << "\n" << std::endl;
ASSERT_EQ(CV_32FC2, r.type());
checkRoot<float>(r, 1, 0);
checkRoot<float>(r, -1, 0);
checkRoot<float>(r, 0, 0);
checkRoot<float>(r, 0, 0);
}
// x^2 - 2x + 1 = 0, roots: 1, 1
coefs = (cv::Mat_<float>(1,3) << 1, -2, 1 );
{
cv::Mat r;
double prec;
prec = cv::solvePoly(coefs, r);
EXPECT_LE(prec, 1e-6);
EXPECT_EQ(2, r.total());
//std::cout << "Preciseness = " << prec << std::endl;
//std::cout << "roots:\n" << r << "\n" << std::endl;
ASSERT_EQ(CV_32FC2, r.type());
checkRoot<float>(r, 1, 0);
checkRoot<float>(r, 1, 0);
}
}
class Core_PhaseTest : public cvtest::BaseTest
{
public: