Chiebot-Mirror
/
opencv
mirror of https://github.com/opencv/opencv.git
8
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge pull request #11969 from alalek:core_Matx_inv_solve_templates

Browse Source
pull/12029/head
Alexander Alekhin 6 years ago
parent 1dec480974 3c74fde349
commit a5e8ae2183
2 changed files with 96 additions and 49 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 87
      modules/core/include/opencv2/core/operations.hpp
  2. 58
      modules/core/test/test_math.cpp

87
modules/core/include/opencv2/core/operations.hpp
Unescape View File

@ -63,9 +63,9 @@ namespace internal
template<typename _Tp, int m, int n> struct Matx_FastInvOp template<typename _Tp, int m, int n> struct Matx_FastInvOp
{ {
bool operator()(const Matx<_Tp, m, n>&, Matx<_Tp, n, m>&, int) const bool operator()(const Matx<_Tp, m, n>& a, Matx<_Tp, n, m>& b, int method) const
{ {
return false; return invert(a, b, method) != 0;
} }
}; };
@ -73,25 +73,32 @@ template<typename _Tp, int m> struct Matx_FastInvOp<_Tp, m, m>
{ {
bool operator()(const Matx<_Tp, m, m>& a, Matx<_Tp, m, m>& b, int method) const bool operator()(const Matx<_Tp, m, m>& a, Matx<_Tp, m, m>& b, int method) const
{ {
Matx<_Tp, m, m> temp = a; if (method == DECOMP_LU || method == DECOMP_CHOLESKY)
{
Matx<_Tp, m, m> temp = a;
// assume that b is all 0's on input => make it a unity matrix // assume that b is all 0's on input => make it a unity matrix
for( int i = 0; i < m; i++ ) for (int i = 0; i < m; i++)
b(i, i) = (_Tp)1; b(i, i) = (_Tp)1;
if( method == DECOMP_CHOLESKY ) if (method == DECOMP_CHOLESKY)
return Cholesky(temp.val, m*sizeof(_Tp), m, b.val, m*sizeof(_Tp), m); return Cholesky(temp.val, m*sizeof(_Tp), m, b.val, m*sizeof(_Tp), m);
return LU(temp.val, m*sizeof(_Tp), m, b.val, m*sizeof(_Tp), m) != 0; return LU(temp.val, m*sizeof(_Tp), m, b.val, m*sizeof(_Tp), m) != 0;
}
else
{
return invert(a, b, method) != 0;
}
} }
}; };
template<typename _Tp> struct Matx_FastInvOp<_Tp, 2, 2> template<typename _Tp> struct Matx_FastInvOp<_Tp, 2, 2>
{ {
bool operator()(const Matx<_Tp, 2, 2>& a, Matx<_Tp, 2, 2>& b, int) const bool operator()(const Matx<_Tp, 2, 2>& a, Matx<_Tp, 2, 2>& b, int /*method*/) const
{ {
_Tp d = (_Tp)determinant(a); _Tp d = (_Tp)determinant(a);
if( d == 0 ) if (d == 0)
return false; return false;
d = 1/d; d = 1/d;
b(1,1) = a(0,0)*d; b(1,1) = a(0,0)*d;
@ -104,10 +111,10 @@ template<typename _Tp> struct Matx_FastInvOp<_Tp, 2, 2>
template<typename _Tp> struct Matx_FastInvOp<_Tp, 3, 3> template<typename _Tp> struct Matx_FastInvOp<_Tp, 3, 3>
{ {
bool operator()(const Matx<_Tp, 3, 3>& a, Matx<_Tp, 3, 3>& b, int) const bool operator()(const Matx<_Tp, 3, 3>& a, Matx<_Tp, 3, 3>& b, int /*method*/) const
{ {
_Tp d = (_Tp)determinant(a); _Tp d = (_Tp)determinant(a);
if( d == 0 ) if (d == 0)
return false; return false;
d = 1/d; d = 1/d;
b(0,0) = (a(1,1) * a(2,2) - a(1,2) * a(2,1)) * d; b(0,0) = (a(1,1) * a(2,2) - a(1,2) * a(2,1)) * d;
@ -128,10 +135,10 @@ template<typename _Tp> struct Matx_FastInvOp<_Tp, 3, 3>
template<typename _Tp, int m, int l, int n> struct Matx_FastSolveOp template<typename _Tp, int m, int l, int n> struct Matx_FastSolveOp
{ {
bool operator()(const Matx<_Tp, m, l>&, const Matx<_Tp, m, n>&, bool operator()(const Matx<_Tp, m, l>& a, const Matx<_Tp, m, n>& b,
Matx<_Tp, l, n>&, int) const Matx<_Tp, l, n>& x, int method) const
{ {
return false; return cv::solve(a, b, x, method);
} }
}; };
@ -140,12 +147,19 @@ template<typename _Tp, int m, int n> struct Matx_FastSolveOp<_Tp, m, m, n>
bool operator()(const Matx<_Tp, m, m>& a, const Matx<_Tp, m, n>& b, bool operator()(const Matx<_Tp, m, m>& a, const Matx<_Tp, m, n>& b,
Matx<_Tp, m, n>& x, int method) const Matx<_Tp, m, n>& x, int method) const
{ {
Matx<_Tp, m, m> temp = a; if (method == DECOMP_LU || method == DECOMP_CHOLESKY)
x = b; {
if( method == DECOMP_CHOLESKY ) Matx<_Tp, m, m> temp = a;
return Cholesky(temp.val, m*sizeof(_Tp), m, x.val, n*sizeof(_Tp), n); x = b;
if( method == DECOMP_CHOLESKY )
return Cholesky(temp.val, m*sizeof(_Tp), m, x.val, n*sizeof(_Tp), n);
return LU(temp.val, m*sizeof(_Tp), m, x.val, n*sizeof(_Tp), n) != 0; return LU(temp.val, m*sizeof(_Tp), m, x.val, n*sizeof(_Tp), n) != 0;
}
else
{
return cv::solve(a, b, x, method);
}
} }
}; };
@ -155,7 +169,7 @@ template<typename _Tp> struct Matx_FastSolveOp<_Tp, 2, 2, 1>
Matx<_Tp, 2, 1>& x, int) const Matx<_Tp, 2, 1>& x, int) const
{ {
_Tp d = (_Tp)determinant(a); _Tp d = (_Tp)determinant(a);
if( d == 0 ) if (d == 0)
return false; return false;
d = 1/d; d = 1/d;
x(0) = (b(0)*a(1,1) - b(1)*a(0,1))*d; x(0) = (b(0)*a(1,1) - b(1)*a(0,1))*d;
@ -170,7 +184,7 @@ template<typename _Tp> struct Matx_FastSolveOp<_Tp, 3, 3, 1>
Matx<_Tp, 3, 1>& x, int) const Matx<_Tp, 3, 1>& x, int) const
{ {
_Tp d = (_Tp)determinant(a); _Tp d = (_Tp)determinant(a);
if( d == 0 ) if (d == 0)
return false; return false;
d = 1/d; d = 1/d;
x(0) = d*(b(0)*(a(1,1)*a(2,2) - a(1,2)*a(2,1)) - x(0) = d*(b(0)*(a(1,1)*a(2,2) - a(1,2)*a(2,1)) -
@ -210,18 +224,8 @@ template<typename _Tp, int m, int n> inline
Matx<_Tp, n, m> Matx<_Tp, m, n>::inv(int method, bool *p_is_ok /*= NULL*/) const Matx<_Tp, n, m> Matx<_Tp, m, n>::inv(int method, bool *p_is_ok /*= NULL*/) const
{ {
Matx<_Tp, n, m> b; Matx<_Tp, n, m> b;
bool ok; bool ok = cv::internal::Matx_FastInvOp<_Tp, m, n>()(*this, b, method);
if (method == DECOMP_LU || method == DECOMP_CHOLESKY) if (p_is_ok) *p_is_ok = ok;
{
CV_Assert(m == n);
ok = cv::internal::Matx_FastInvOp<_Tp, m, n>()(*this, b, method);
}
else
{
Mat A(*this, false), B(b, false);
ok = (invert(A, B, method) != 0);
}
if( NULL != p_is_ok ) { *p_is_ok = ok; }
return ok ? b : Matx<_Tp, n, m>::zeros(); return ok ? b : Matx<_Tp, n, m>::zeros();
} }
@ -229,18 +233,7 @@ template<typename _Tp, int m, int n> template<int l> inline
Matx<_Tp, n, l> Matx<_Tp, m, n>::solve(const Matx<_Tp, m, l>& rhs, int method) const Matx<_Tp, n, l> Matx<_Tp, m, n>::solve(const Matx<_Tp, m, l>& rhs, int method) const
{ {
Matx<_Tp, n, l> x; Matx<_Tp, n, l> x;
bool ok; bool ok = cv::internal::Matx_FastSolveOp<_Tp, m, n, l>()(*this, rhs, x, method);
if (method == DECOMP_LU || method == DECOMP_CHOLESKY)
{
CV_Assert(m == n);
ok = cv::internal::Matx_FastSolveOp<_Tp, m, n, l>()(*this, rhs, x, method);
}
else
{
Mat A(*this, false), B(rhs, false), X(x, false);
ok = cv::solve(A, B, X, method);
}
return ok ? x : Matx<_Tp, n, l>::zeros(); return ok ? x : Matx<_Tp, n, l>::zeros();
} }

58
modules/core/test/test_math.cpp
Unescape View File

@ -3139,9 +3139,63 @@ TEST(Core_Solve, regression_11888)
cv::Vec<float, 3> b(4, 5, 7); cv::Vec<float, 3> b(4, 5, 7);
cv::Matx<float, 2, 1> xQR = A.solve(b, DECOMP_QR); cv::Matx<float, 2, 1> xQR = A.solve(b, DECOMP_QR);
cv::Matx<float, 2, 1> xSVD = A.solve(b, DECOMP_SVD); cv::Matx<float, 2, 1> xSVD = A.solve(b, DECOMP_SVD);
EXPECT_LE(cvtest::norm(xQR, xSVD, CV_RELATIVE_L2), 0.001); EXPECT_LE(cvtest::norm(xQR, xSVD, NORM_L2 | NORM_RELATIVE), 0.001);
cv::Matx<float, 2, 3> iA = A.inv(DECOMP_SVD); cv::Matx<float, 2, 3> iA = A.inv(DECOMP_SVD);
EXPECT_LE(cvtest::norm(A*iA, Matx<float, 3, 3>::eye(), CV_RELATIVE_L2), 0.6); EXPECT_LE(cvtest::norm(iA*A, Matx<float, 2, 2>::eye(), NORM_L2), 1e-3);
EXPECT_ANY_THROW({
/*cv::Matx<float, 2, 1> xLU =*/ A.solve(b, DECOMP_LU);
std::cout << "FATAL ERROR" << std::endl;
});
}
TEST(Core_Solve, Matx_2_2)
{
cv::Matx<float, 2, 2> A(
2, 1,
1, 1
);
cv::Vec<float, 2> b(4, 5);
cv::Matx<float, 2, 1> xLU = A.solve(b, DECOMP_LU);
cv::Matx<float, 2, 1> xQR = A.solve(b, DECOMP_QR);
cv::Matx<float, 2, 1> xSVD = A.solve(b, DECOMP_SVD);
EXPECT_LE(cvtest::norm(xQR, xSVD, NORM_L2 | NORM_RELATIVE), 1e-3);
EXPECT_LE(cvtest::norm(xQR, xLU, NORM_L2 | NORM_RELATIVE), 1e-3);
cv::Matx<float, 2, 2> iA = A.inv(DECOMP_SVD);
EXPECT_LE(cvtest::norm(iA*A, Matx<float, 2, 2>::eye(), NORM_L2), 1e-3);
}
TEST(Core_Solve, Matx_3_3)
{
cv::Matx<float, 3, 3> A(
2, 1, 0,
0, 1, 1,
1, 0, 1
);
cv::Vec<float, 3> b(4, 5, 6);
cv::Matx<float, 3, 1> xLU = A.solve(b, DECOMP_LU);
cv::Matx<float, 3, 1> xQR = A.solve(b, DECOMP_QR);
cv::Matx<float, 3, 1> xSVD = A.solve(b, DECOMP_SVD);
EXPECT_LE(cvtest::norm(xQR, xSVD, NORM_L2 | NORM_RELATIVE), 1e-3);
EXPECT_LE(cvtest::norm(xQR, xLU, NORM_L2 | NORM_RELATIVE), 1e-3);
cv::Matx<float, 3, 3> iA = A.inv(DECOMP_SVD);
EXPECT_LE(cvtest::norm(iA*A, Matx<float, 3, 3>::eye(), NORM_L2), 1e-3);
}
TEST(Core_Solve, Matx_4_4)
{
cv::Matx<float, 4, 4> A(
2, 1, 0, 4,
0, 1, 1, 3,
1, 0, 1, 2,
2, 2, 0, 1
);
cv::Vec<float, 4> b(4, 5, 6, 7);
cv::Matx<float, 4, 1> xLU = A.solve(b, DECOMP_LU);
cv::Matx<float, 4, 1> xQR = A.solve(b, DECOMP_QR);
cv::Matx<float, 4, 1> xSVD = A.solve(b, DECOMP_SVD);
EXPECT_LE(cvtest::norm(xQR, xSVD, NORM_L2 | NORM_RELATIVE), 1e-3);
EXPECT_LE(cvtest::norm(xQR, xLU, NORM_L2 | NORM_RELATIVE), 1e-3);
cv::Matx<float, 4, 4> iA = A.inv(DECOMP_SVD);
EXPECT_LE(cvtest::norm(iA*A, Matx<float, 4, 4>::eye(), NORM_L2), 1e-3);
} }
softdouble naiveExp(softdouble x) softdouble naiveExp(softdouble x)

Write Preview
Loading…
Cancel
Save