From d72d11d09da7d9f6d77a88bff4ce33334c561159 Mon Sep 17 00:00:00 2001 From: Vadim Pisarevsky Date: Wed, 22 Sep 2010 13:07:51 +0000 Subject: [PATCH] improved matrix expressions efficiency in some cases & simplified the code --- modules/core/src/matop.cpp | 428 +++++++++++++++++---------------------------- 1 file changed, 163 insertions(+), 265 deletions(-) diff --git a/modules/core/src/matop.cpp b/modules/core/src/matop.cpp index 9ec5cde..f249992 100644 --- a/modules/core/src/matop.cpp +++ b/modules/core/src/matop.cpp @@ -59,18 +59,7 @@ public: bool elementWise(const MatExpr& expr) const { return true; } void assign(const MatExpr& expr, Mat& m, int type=-1) const; - - void add(const MatExpr& e1, const MatExpr& e2, MatExpr& res) const; - void add(const MatExpr& e1, const Scalar& s, MatExpr& res) const; - void subtract(const MatExpr& e1, const MatExpr& e2, MatExpr& res) const; - void subtract(const Scalar& s, const MatExpr& expr, MatExpr& res) const; - void multiply(const MatExpr& e1, const MatExpr& e2, MatExpr& res, double scale=1) const; - void multiply(const MatExpr& e1, double s, MatExpr& res) const; - void divide(const MatExpr& e1, const MatExpr& e2, MatExpr& res, double scale=1) const; - void divide(double s, const MatExpr& e, MatExpr& res) const; - - void matmul(const MatExpr& expr1, const MatExpr& expr2, MatExpr& res) const; - + static void makeExpr(MatExpr& res, const Mat& m); }; @@ -85,21 +74,13 @@ public: bool elementWise(const MatExpr& expr) const { return true; } void assign(const MatExpr& expr, Mat& m, int type=-1) const; - void add(const MatExpr& e1, const MatExpr& e2, MatExpr& res) const; void add(const MatExpr& e1, const Scalar& s, MatExpr& res) const; - - void subtract(const MatExpr& e1, const MatExpr& e2, MatExpr& res) const; void subtract(const Scalar& s, const MatExpr& expr, MatExpr& res) const; - - void multiply(const MatExpr& e1, const MatExpr& e2, MatExpr& res, double scale=1) const; void multiply(const MatExpr& e1, double s, MatExpr& res) const; - - void divide(const MatExpr& e1, const MatExpr& e2, MatExpr& res, double scale=1) const; void divide(double s, const MatExpr& e, MatExpr& res) const; void transpose(const MatExpr& e1, MatExpr& res) const; void abs(const MatExpr& expr, MatExpr& res) const; - void matmul(const MatExpr& expr1, const MatExpr& expr2, MatExpr& res) const; static void makeExpr(MatExpr& res, const Mat& a, const Mat& b, double alpha, double beta, const Scalar& s=Scalar()); }; @@ -115,10 +96,7 @@ public: bool elementWise(const MatExpr& expr) const { return true; } void assign(const MatExpr& expr, Mat& m, int type=-1) const; - void multiply(const MatExpr& e1, const MatExpr& e2, MatExpr& res, double scale=1) const; void multiply(const MatExpr& e1, double s, MatExpr& res) const; - - void divide(const MatExpr& e1, const MatExpr& e2, MatExpr& res, double scale=1) const; void divide(double s, const MatExpr& e, MatExpr& res) const; static void makeExpr(MatExpr& res, char op, const Mat& a, const Mat& b, double scale=1); @@ -190,7 +168,6 @@ public: void multiply(const MatExpr& e1, double s, MatExpr& res) const; void transpose(const MatExpr& expr, MatExpr& res) const; - void matmul(const MatExpr& expr1, const MatExpr& expr2, MatExpr& res) const; static void makeExpr(MatExpr& res, const Mat& a, double alpha=1); }; @@ -344,17 +321,34 @@ void MatOp::augAssignXor(const MatExpr& expr, Mat& m) const } -void MatOp::add(const MatExpr& expr1, const MatExpr& expr2, MatExpr& res) const +void MatOp::add(const MatExpr& e1, const MatExpr& e2, MatExpr& res) const { - if( this == expr2.op ) + if( this == e2.op ) { + double alpha = 1, beta = 1; + Scalar s; Mat m1, m2; - expr1.op->assign(expr1, m1); - expr2.op->assign(expr2, m2); - MatOp_AddEx::makeExpr(res, m1, m2, 1, 1); + if( isAddEx(e1) && (!e1.b.data || e1.beta == 0) ) + { + m1 = e1.a; + alpha = e1.alpha; + s = e1.s; + } + else + e1.op->assign(e1, m1); + + if( isAddEx(e2) && (!e2.b.data || e2.beta == 0) ) + { + m2 = e2.a; + beta = e2.alpha; + s += e2.s; + } + else + e2.op->assign(e2, m2); + MatOp_AddEx::makeExpr(res, m1, m2, alpha, beta, s); } else - expr2.op->add(expr1, expr2, res); + e2.op->add(e1, e2, res); } @@ -366,17 +360,34 @@ void MatOp::add(const MatExpr& expr1, const Scalar& s, MatExpr& res) const } -void MatOp::subtract(const MatExpr& expr1, const MatExpr& expr2, MatExpr& res) const +void MatOp::subtract(const MatExpr& e1, const MatExpr& e2, MatExpr& res) const { - if( this == expr2.op ) + if( this == e2.op ) { + double alpha = 1, beta = -1; + Scalar s; Mat m1, m2; - expr1.op->assign(expr1, m1); - expr2.op->assign(expr2, m2); - MatOp_AddEx::makeExpr(res, m1, m2, 1, -1); + if( isAddEx(e1) && (!e1.b.data || e1.beta == 0) ) + { + m1 = e1.a; + alpha = e1.alpha; + s = e1.s; + } + else + e1.op->assign(e1, m1); + + if( isAddEx(e2) && (!e2.b.data || e2.beta == 0) ) + { + m2 = e2.a; + beta = -e2.alpha; + s -= e2.s; + } + else + e2.op->assign(e2, m2); + MatOp_AddEx::makeExpr(res, m1, m2, alpha, beta, s); } else - expr2.op->subtract(expr1, expr2, res); + e2.op->subtract(e1, e2, res); } @@ -388,17 +399,54 @@ void MatOp::subtract(const Scalar& s, const MatExpr& expr, MatExpr& res) const } -void MatOp::multiply(const MatExpr& expr1, const MatExpr& expr2, MatExpr& res, double scale) const +void MatOp::multiply(const MatExpr& e1, const MatExpr& e2, MatExpr& res, double scale) const { - if( this == expr2.op ) + if( this == e2.op ) { Mat m1, m2; - expr1.op->assign(expr1, m1); - expr2.op->assign(expr2, m2); - MatOp_Bin::makeExpr(res, '*', m1, m2, scale); + + if( isReciprocal(e1) ) + { + if( isScaled(e2) ) + { + scale *= e2.alpha; + m2 = e2.a; + } + else + e2.op->assign(e2, m2); + + MatOp_Bin::makeExpr(res, '/', m2, e1.a, scale/e1.alpha); + } + else + { + char op = '*'; + if( isScaled(e1) ) + { + m1 = e1.a; + scale *= e1.alpha; + } + else + e1.op->assign(e1, m1); + + if( isScaled(e2) ) + { + m2 = e2.a; + scale *= e2.alpha; + } + else if( isReciprocal(e2) ) + { + op = '/'; + m2 = e2.a; + scale /= e2.alpha; + } + else + e2.op->assign(e2, m2); + + MatOp_Bin::makeExpr(res, op, m1, m2, scale); + } } else - expr2.op->multiply(expr1, expr2, res, scale); + e2.op->multiply(e1, e2, res, scale); } @@ -410,17 +458,43 @@ void MatOp::multiply(const MatExpr& expr, double s, MatExpr& res) const } -void MatOp::divide(const MatExpr& expr1, const MatExpr& expr2, MatExpr& res, double scale) const +void MatOp::divide(const MatExpr& e1, const MatExpr& e2, MatExpr& res, double scale) const { - if( this == expr2.op ) + if( this == e2.op ) { - Mat m1, m2; - expr1.op->assign(expr1, m1); - expr2.op->assign(expr2, m2); - MatOp_Bin::makeExpr(res, '/', m1, m2, scale); + if( isReciprocal(e1) && isReciprocal(e2) ) + MatOp_Bin::makeExpr(res, '/', e2.a, e1.a, e1.alpha/e2.alpha); + else + { + Mat m1, m2; + char op = '/'; + + if( isScaled(e1) ) + { + m1 = e1.a; + scale *= e1.alpha; + } + else + e1.op->assign(e1, m1); + + if( isScaled(e2) ) + { + m2 = e2.a; + scale /= e2.alpha; + } + else if( isReciprocal(e2) ) + { + m2 = e2.a; + scale /= e2.alpha; + op = '*'; + } + else + e2.op->assign(e2, m2); + MatOp_Bin::makeExpr(res, op, m1, m2, scale); + } } else - expr2.op->divide(expr1, expr2, res, scale); + e2.op->divide(e1, e2, res, scale); } @@ -448,17 +522,46 @@ void MatOp::transpose(const MatExpr& expr, MatExpr& res) const } -void MatOp::matmul(const MatExpr& expr1, const MatExpr& expr2, MatExpr& res) const +void MatOp::matmul(const MatExpr& e1, const MatExpr& e2, MatExpr& res) const { - if( this == expr2.op ) + if( this == e2.op ) { + double scale = 1; + int flags = 0; Mat m1, m2; - expr1.op->assign(expr1, m1); - expr2.op->assign(expr2, m2); - MatOp_GEMM::makeExpr(res, 0, m1, m2); + + if( isT(e1) ) + { + flags = CV_GEMM_A_T; + scale = e1.alpha; + m1 = e1.a; + } + else if( isScaled(e1) ) + { + scale = e1.alpha; + m1 = e1.a; + } + else + e1.op->assign(e1, m1); + + if( isT(e2) ) + { + flags |= CV_GEMM_B_T; + scale *= e2.alpha; + m2 = e2.a; + } + else if( isScaled(e2) ) + { + scale *= e2.alpha; + m2 = e2.a; + } + else + e2.op->assign(e2, m2); + + MatOp_GEMM::makeExpr(res, flags, m1, m2, scale); } else - expr2.op->matmul(expr1, expr2, res); + e2.op->matmul(e1, e2, res); } @@ -1041,11 +1144,6 @@ MatExpr abs(const MatExpr& e) ///////////////////////////////////////////////////////////////////////////////////////////////////// -inline void MatOp_Identity::makeExpr(MatExpr& res, const Mat& m) -{ - res = MatExpr(&g_MatOp_Identity, 0, m, Mat(), Mat(), 1, 0); -} - void MatOp_Identity::assign(const MatExpr& e, Mat& m, int type) const { if( type == -1 || type == e.a.type() ) @@ -1057,91 +1155,11 @@ void MatOp_Identity::assign(const MatExpr& e, Mat& m, int type) const } } -void MatOp_Identity::add(const MatExpr& e1, const MatExpr& e2, MatExpr& res) const -{ - if( isIdentity(e2) ) - { - if( isIdentity(e1) ) - MatOp_AddEx::makeExpr(res, e1.a, e2.a, 1, 1); - else - MatOp::add(e1, e2, res); - } - else - e2.op->add(e1, e2, res); -} - -void MatOp_Identity::add(const MatExpr& e, const Scalar& s, MatExpr& res) const +inline void MatOp_Identity::makeExpr(MatExpr& res, const Mat& m) { - MatOp_AddEx::makeExpr(res, e.a, Mat(), 1, 0, s); + res = MatExpr(&g_MatOp_Identity, 0, m, Mat(), Mat(), 1, 0); } -void MatOp_Identity::subtract(const MatExpr& e1, const MatExpr& e2, MatExpr& res) const -{ - if( isIdentity(e2) ) - { - if( isIdentity(e1) ) - MatOp_AddEx::makeExpr(res, e1.a, e2.a, 1, -1); - else - MatOp::subtract(e1, e2, res); - } - else - e2.op->subtract(e1, e2, res); -} - -void MatOp_Identity::subtract(const Scalar& s, const MatExpr& e, MatExpr& res) const -{ - MatOp_AddEx::makeExpr(res, e.a, Mat(), -1, 0, s); -} - -void MatOp_Identity::multiply(const MatExpr& e1, const MatExpr& e2, MatExpr& res, double scale) const -{ - if( isIdentity(e2) ) - { - if( isIdentity(e1) ) - MatOp_Bin::makeExpr(res, '*', e1.a, e2.a, scale); - else - MatOp::multiply(e1, e2, res, scale); - } - else - e2.op->multiply(e1, e2, res, scale); -} - -void MatOp_Identity::multiply(const MatExpr& e, double s, MatExpr& res) const -{ - MatOp_AddEx::makeExpr(res, e.a, Mat(), s, 0); -} - -void MatOp_Identity::divide(const MatExpr& e1, const MatExpr& e2, MatExpr& res, double scale) const -{ - if( isIdentity(e2) ) - { - if( isIdentity(e1) ) - MatOp_Bin::makeExpr(res, '/', e1.a, e2.a, scale); - else - MatOp::divide(e1, e2, res, scale); - } - else - e2.op->divide(e1, e2, res, scale); -} - -void MatOp_Identity::divide(double s, const MatExpr& e, MatExpr& res) const -{ - MatOp_Bin::makeExpr(res, '/', e.a, Mat(), s); -} - -void MatOp_Identity::matmul(const MatExpr& e1, const MatExpr& e2, MatExpr& res) const -{ - if( isIdentity(e2) ) - { - if( isIdentity(e1) ) - MatOp_GEMM::makeExpr(res, 0, e1.a, e2.a); - else - MatOp::matmul(e1, e2, res); - } - else - e2.op->matmul(e1, e2, res); -} - ///////////////////////////////////////////////////////////////////////////////////////////////////// void MatOp_AddEx::assign(const MatExpr& e, Mat& m, int type) const @@ -1196,41 +1214,13 @@ void MatOp_AddEx::assign(const MatExpr& e, Mat& m, int type) const } -void MatOp_AddEx::add(const MatExpr& e1, const MatExpr& e2, MatExpr& res) const -{ - bool i1 = isIdentity(e1), i2 = isIdentity(e2); - - if( ((isAddEx(e2) && !e2.b.data) || i2) && ((isAddEx(e1) && !e1.b.data) || i1) ) - MatOp_AddEx::makeExpr(res, e1.a, e2.a, i1 ? 1 : e1.alpha, i2 ? 1 : e2.alpha, - (i1 ? Scalar() : e1.s) + (i2 ? Scalar() : e2.s)); - else if( e2.op == this ) - MatOp::add(e1, e2, res); - else - e2.op->add(e1, e2, res); -} - - void MatOp_AddEx::add(const MatExpr& e, const Scalar& s, MatExpr& res) const { res = e; res.s += s; } - -void MatOp_AddEx::subtract(const MatExpr& e1, const MatExpr& e2, MatExpr& res) const -{ - bool i1 = isIdentity(e1), i2 = isIdentity(e2); - - if( ((isAddEx(e2) && !e2.b.data) || i2) && ((isAddEx(e1) && !e1.b.data) || i1) ) - MatOp_AddEx::makeExpr(res, e1.a, e2.a, i1 ? 1 : e1.alpha, i2 ? -1 : -e2.alpha, - (i1 ? Scalar() : e1.s) - (i2 ? Scalar() : e2.s)); - else if( e2.op == this ) - MatOp::subtract(e1, e2, res); - else - e2.op->subtract(e1, e2, res); -} - - + void MatOp_AddEx::subtract(const Scalar& s, const MatExpr& e, MatExpr& res) const { res = e; @@ -1238,21 +1228,7 @@ void MatOp_AddEx::subtract(const Scalar& s, const MatExpr& e, MatExpr& res) cons res.beta = -res.beta; res.s = s - res.s; } - - -void MatOp_AddEx::multiply(const MatExpr& e1, const MatExpr& e2, MatExpr& res, double scale) const -{ - bool i1 = isIdentity(e1), i2 = isIdentity(e2); - - if( (isScaled(e1) || i1) && (isScaled(e2) || i2) ) - MatOp_Bin::makeExpr(res, '*', e1.a, e2.a, scale*(i1 ? 1 : e1.alpha)*(i2 ? 1 : e2.alpha)); - else if( this == e2.op ) - MatOp::multiply(e1, e2, res, scale); - else - e2.op->multiply(e1, e2, res, scale); -} - - + void MatOp_AddEx::multiply(const MatExpr& e, double s, MatExpr& res) const { res = e; @@ -1260,20 +1236,6 @@ void MatOp_AddEx::multiply(const MatExpr& e, double s, MatExpr& res) const res.beta *= s; res.s *= s; } - - -void MatOp_AddEx::divide(const MatExpr& e1, const MatExpr& e2, MatExpr& res, double scale) const -{ - bool i1 = isIdentity(e1), i2 = isIdentity(e2); - - if( (isScaled(e1) || i1) && (isScaled(e2) || i2) ) - MatOp_Bin::makeExpr(res, '/', e1.a, e2.a, scale*(i1 ? 1 : e1.alpha)*(i2 ? 1 : 1./e2.alpha)); - else if( this == e2.op ) - MatOp::divide(e1, e2, res, scale); - else - e2.op->divide(e1, e2, res, scale); -} - void MatOp_AddEx::divide(double s, const MatExpr& e, MatExpr& res) const { @@ -1301,20 +1263,6 @@ void MatOp_AddEx::abs(const MatExpr& e, MatExpr& res) const else MatOp::abs(e, res); } - -void MatOp_AddEx::matmul(const MatExpr& e1, const MatExpr& e2, MatExpr& res) const -{ - bool i1 = isIdentity(e1), i2 = isIdentity(e2); - double alpha1 = i1 ? 1 : e1.alpha, alpha2 = i2 ? 1 : e2.alpha; - - if( (isScaled(e1) || i1) && (isScaled(e2) || i2) ) - MatOp_GEMM::makeExpr(res, 0, e1.a, e2.a, alpha1*alpha2); - else if( this == e2.op ) - MatOp::matmul(e1, e2, res); - else - e2.op->matmul(e1, e2, res); -} - inline void MatOp_AddEx::makeExpr(MatExpr& res, const Mat& a, const Mat& b, double alpha, double beta, const Scalar& s) { @@ -1366,18 +1314,6 @@ void MatOp_Bin::assign(const MatExpr& e, Mat& m, int type) const dst.convertTo(m, type); } -void MatOp_Bin::multiply(const MatExpr& e1, const MatExpr& e2, MatExpr& res, double scale) const -{ - if( isReciprocal(e1) && (isIdentity(e2) || isScaled(e2)) ) - MatOp_Bin::makeExpr(res, '/', e2.a, e1.a, e1.alpha*(isIdentity(e2) ? 1 : e2.alpha)*scale); - else if( isReciprocal(e2) && (isIdentity(e1) || isScaled(e1)) ) - MatOp_Bin::makeExpr(res, '/', e1.a, e2.a, (isIdentity(e1) ? 1 : e1.alpha)*e2.alpha*scale); - else if( e2.op == this ) - MatOp::multiply(e1, e2, res, scale); - else - e2.op->multiply(e1, e2, res, scale); -} - void MatOp_Bin::multiply(const MatExpr& e, double s, MatExpr& res) const { if( e.flags == '*' || e.flags == '/' ) @@ -1389,30 +1325,9 @@ void MatOp_Bin::multiply(const MatExpr& e, double s, MatExpr& res) const MatOp::multiply(e, s, res); } -void MatOp_Bin::divide(const MatExpr& e1, const MatExpr& e2, MatExpr& res, double scale) const -{ - if( isReciprocal(e2) ) - { - if( isIdentity(e1) || isScaled(e1) ) - { - MatOp_Bin::makeExpr(res, '*', e1.a, e2.a, (isIdentity(e1) ? 1 : e1.alpha)*scale/e2.alpha); - return; - } - if( isReciprocal(e1) ) - { - MatOp_Bin::makeExpr(res, '/', e2.a, e1.a, e1.alpha*scale/e2.alpha); - return; - } - } - if( e2.op == this ) - MatOp::divide(e1, e2, res, scale); - else - e2.op->divide(e1, e2, res, scale); -} - void MatOp_Bin::divide(double s, const MatExpr& e, MatExpr& res) const { - if( e.flags == '/' && !e.b.data ) + if( e.flags == '/' && (!e.b.data || e.beta == 0) ) MatOp_AddEx::makeExpr(res, e.a, Mat(), s/e.alpha, 0); else MatOp::divide(s, e, res); @@ -1479,23 +1394,6 @@ void MatOp_T::transpose(const MatExpr& e, MatExpr& res) const MatOp_AddEx::makeExpr(res, e.a, Mat(), e.alpha, 0); } -void MatOp_T::matmul(const MatExpr& e1, const MatExpr& e2, MatExpr& res) const -{ - bool i1 = isIdentity(e1), i2 = isIdentity(e2); - double alpha1 = i1 ? 1 : e1.alpha, alpha2 = i2 ? 1 : e2.alpha; - - if( isT(e1) && (isT(e2) || isScaled(e2) || i2)) - MatOp_GEMM::makeExpr(res, CV_GEMM_A_T + (isT(e2) ? CV_GEMM_B_T : 0), - e1.a, e2.a, e1.alpha*alpha2); - else if( isT(e2) && (isT(e1) || isScaled(e1) || i1)) - MatOp_GEMM::makeExpr(res, CV_GEMM_B_T + (isT(e1) ? CV_GEMM_A_T : 0), - e1.a, e2.a, e2.alpha*alpha1); - else if( this == e2.op ) - MatOp::matmul(e1, e2, res); - else - e2.op->matmul(e1, e2, res); -} - inline void MatOp_T::makeExpr(MatExpr& res, const Mat& a, double alpha) { res = MatExpr(&g_MatOp_T, 0, a, Mat(), Mat(), alpha, 0); -- 2.7.4