From: Vadim Pisarevsky Date: Tue, 29 Jun 2010 14:52:43 +0000 (+0000) Subject: added "small matrix" class Matx X-Git-Tag: accepted/2.0/20130307.220821~4811 X-Git-Url: http://review.tizen.org/git/?a=commitdiff_plain;h=10b5a51731610a1db2f550c4da9e9bb98c584c69;p=profile%2Fivi%2Fopencv.git added "small matrix" class Matx --- diff --git a/modules/core/include/opencv2/core/core.hpp b/modules/core/include/opencv2/core/core.hpp index fe76ded..a51e3c7 100644 --- a/modules/core/include/opencv2/core/core.hpp +++ b/modules/core/include/opencv2/core/core.hpp @@ -78,15 +78,51 @@ using std::string; template class CV_EXPORTS Size_; template class CV_EXPORTS Point_; template class CV_EXPORTS Rect_; - +template class CV_EXPORTS Vec; +template class CV_EXPORTS Matx; + typedef std::string String; typedef std::basic_string WString; +class Mat; +class MatND; +template class CV_EXPORTS MatExpr_Base_; +typedef MatExpr_Base_ MatExpr_Base; +template class MatExpr_; +template class MatExpr_Op1_; +template class MatExpr_Op2_; +template class MatExpr_Op3_; +template class MatExpr_Op4_; +template class MatExpr_Op5_; +template class CV_EXPORTS MatOp_DivRS_; +template class CV_EXPORTS MatOp_Inv_; +template class CV_EXPORTS MatOp_MulDiv_; +template class CV_EXPORTS MatOp_Repeat_; +template class CV_EXPORTS MatOp_Set_; +template class CV_EXPORTS MatOp_Scale_; +template class CV_EXPORTS MatOp_T_; + +template class CV_EXPORTS MatIterator_; +template class CV_EXPORTS MatConstIterator_; +template class CV_EXPORTS MatCommaInitializer_; + CV_EXPORTS string fromUtf16(const WString& str); CV_EXPORTS WString toUtf16(const string& str); CV_EXPORTS string format( const char* fmt, ... ); + +// matrix decomposition types +enum { DECOMP_LU=0, DECOMP_SVD=1, DECOMP_EIG=2, DECOMP_CHOLESKY=3, DECOMP_QR=4, DECOMP_NORMAL=16 }; +enum { NORM_INF=1, NORM_L1=2, NORM_L2=4, NORM_TYPE_MASK=7, NORM_RELATIVE=8, NORM_MINMAX=32}; +enum { CMP_EQ=0, CMP_GT=1, CMP_GE=2, CMP_LT=3, CMP_LE=4, CMP_NE=5 }; +enum { GEMM_1_T=1, GEMM_2_T=2, GEMM_3_T=4 }; +enum { DFT_INVERSE=1, DFT_SCALE=2, DFT_ROWS=4, DFT_COMPLEX_OUTPUT=16, DFT_REAL_OUTPUT=32, + DCT_INVERSE = DFT_INVERSE, DCT_ROWS=DFT_ROWS }; + + /*! The standard OpenCV exception class. Instances of the class are thrown by various functions and methods in the case of critical errors. @@ -404,6 +440,7 @@ public: Vec(_Tp v0, _Tp v1, _Tp v2, _Tp v3, _Tp v4, _Tp v5, _Tp v6, _Tp v7); //!< 8-element vector constructor Vec(_Tp v0, _Tp v1, _Tp v2, _Tp v3, _Tp v4, _Tp v5, _Tp v6, _Tp v7, _Tp v8); //!< 9-element vector constructor Vec(_Tp v0, _Tp v1, _Tp v2, _Tp v3, _Tp v4, _Tp v5, _Tp v6, _Tp v7, _Tp v8, _Tp v9); //!< 10-element vector constructor + explicit Vec(const _Tp* values); Vec(const Vec<_Tp, cn>& v); static Vec all(_Tp alpha); @@ -421,10 +458,14 @@ public: template operator Vec() const; //! conversion to 4-element CvScalar. operator CvScalar() const; + + Matx<_Tp, 1, cn> t() const; /*! element access */ const _Tp& operator [](int i) const; _Tp& operator[](int i); + const _Tp& operator ()(int i) const; + _Tp& operator ()(int i); _Tp val[cn]; //< vector elements }; @@ -461,6 +502,158 @@ typedef Vec Vec4d; typedef Vec Vec6d; +////////////////////////////// Small Matrix /////////////////////////// + +/*! + A short numerical vector. + + This template class represents short numerical vectors (of 1, 2, 3, 4 ... elements) + on which you can perform basic arithmetical operations, access individual elements using [] operator etc. + The vectors are allocated on stack, as opposite to std::valarray, std::vector, cv::Mat etc., + which elements are dynamically allocated in the heap. + + The template takes 2 parameters: + -# _Tp element type + -# cn the number of elements + + In addition to the universal notation like Vec, you can use shorter aliases + for the most popular specialized variants of Vec, e.g. Vec3f ~ Vec. + */ + +struct CV_EXPORTS Matx_AddOp {}; +struct CV_EXPORTS Matx_SubOp {}; +struct CV_EXPORTS Matx_ScaleOp {}; +struct CV_EXPORTS Matx_MulOp {}; +struct CV_EXPORTS Matx_MatMulOp {}; +struct CV_EXPORTS Matx_TOp {}; + +template class CV_EXPORTS Matx : public Vec<_Tp, m*n> +{ +public: + typedef _Tp value_type; + typedef Vec<_Tp, m*n> base_type; + typedef Vec<_Tp, MIN(m, n)> diag_type; + typedef Matx<_Tp, m, n> mat_type; + enum { depth = DataDepth<_Tp>::value, rows = m, cols = n, channels = rows*cols, + type = CV_MAKETYPE(depth, channels) }; + + //! default constructor + Matx(); + + Matx(_Tp v0); //!< 1x1 matrix + Matx(_Tp v0, _Tp v1); //!< 1x2 or 2x1 matrix + Matx(_Tp v0, _Tp v1, _Tp v2); //!< 1x3 or 3x1 matrix + Matx(_Tp v0, _Tp v1, _Tp v2, _Tp v3); //!< 1x4, 2x2 or 4x1 matrix + Matx(_Tp v0, _Tp v1, _Tp v2, _Tp v3, _Tp v4); //!< 1x5 or 5x1 matrix + Matx(_Tp v0, _Tp v1, _Tp v2, _Tp v3, _Tp v4, _Tp v5); //!< 1x6, 2x3, 3x2 or 6x1 matrix + Matx(_Tp v0, _Tp v1, _Tp v2, _Tp v3, _Tp v4, _Tp v5, _Tp v6); //!< 1x7 or 7x1 matrix + Matx(_Tp v0, _Tp v1, _Tp v2, _Tp v3, _Tp v4, _Tp v5, _Tp v6, _Tp v7); //!< 1x8, 2x4, 4x2 or 8x1 matrix + Matx(_Tp v0, _Tp v1, _Tp v2, _Tp v3, _Tp v4, _Tp v5, _Tp v6, _Tp v7, _Tp v8); //!< 1x9, 3x3 or 9x1 matrix + Matx(_Tp v0, _Tp v1, _Tp v2, _Tp v3, _Tp v4, _Tp v5, _Tp v6, _Tp v7, _Tp v8, _Tp v9); //!< 1x10, 2x5 or 5x2 or 10x1 matrix + Matx(_Tp v0, _Tp v1, _Tp v2, _Tp v3, + _Tp v4, _Tp v5, _Tp v6, _Tp v7, + _Tp v8, _Tp v9, _Tp v10, _Tp v11); //!< 1x12, 2x6, 3x4, 4x3, 6x2 or 12x1 matrix + Matx(_Tp v0, _Tp v1, _Tp v2, _Tp v3, + _Tp v4, _Tp v5, _Tp v6, _Tp v7, + _Tp v8, _Tp v9, _Tp v10, _Tp v11, + _Tp v12, _Tp v13, _Tp v14, _Tp v15); //!< 1x16, 4x4 or 16x1 matrix + explicit Matx(const _Tp* vals); //!< initialize from a plain array + + Matx(const base_type& v); + static Matx all(_Tp alpha); + static Matx zeros(); + static Matx ones(); + static Matx eye(); + static Matx diag(const Vec<_Tp, MIN(m,n)>& d); + static Matx randu(_Tp a, _Tp b); + static Matx randn(_Tp m, _Tp sigma); + + //! convertion to another data type + template operator Matx() const; + + //! change the matrix shape + template Matx<_Tp, m1, n1> reshape() const; + + //! extract part of the matrix + template Matx<_Tp, m1, n1> minor(int i, int j) const; + + //! extract the matrix row + Matx<_Tp, 1, n> row(int i) const; + + //! extract the matrix column + Vec<_Tp, m> col(int i) const; + + //! extract the matrix diagonal + Vec<_Tp, MIN(m,n)> diag() const; + + //! transpose the matrix + Matx<_Tp, n, m> t() const; + + //! invert matrix the matrix + Matx<_Tp, n, m> inv(int method=DECOMP_LU) const; + + //! solve linear system + template Matx<_Tp, n, l> solve(const Matx<_Tp, m, l>& rhs, int flags=DECOMP_LU) const; + Vec<_Tp, n> solve(const Vec<_Tp, m>& rhs, int method) const; + + //! multiply two matrices element-wise + Matx<_Tp, m, n> mul(const Matx<_Tp, m, n>& a) const; + + //! element access + const _Tp& operator ()(int i, int j) const; + _Tp& operator ()(int i, int j); + + //! 1D element access + const _Tp& operator ()(int i) const; + _Tp& operator ()(int i); + + Matx(const Matx<_Tp, m, n>& a, const Matx<_Tp, m, n>& b, Matx_AddOp); + Matx(const Matx<_Tp, m, n>& a, const Matx<_Tp, m, n>& b, Matx_SubOp); + template Matx(const Matx<_Tp, m, n>& a, _T2 alpha, Matx_ScaleOp); + Matx(const Matx<_Tp, m, n>& a, const Matx<_Tp, m, n>& b, Matx_MulOp); + template Matx(const Matx<_Tp, m, l>& a, const Matx<_Tp, l, n>& b, Matx_MatMulOp); + Matx(const Matx<_Tp, n, m>& a, Matx_TOp); +}; + + +typedef Matx Matx12f; +typedef Matx Matx12d; +typedef Matx Matx13f; +typedef Matx Matx13d; +typedef Matx Matx14f; +typedef Matx Matx14d; +typedef Matx Matx16f; +typedef Matx Matx16d; + +typedef Matx Matx21f; +typedef Matx Matx21d; +typedef Matx Matx31f; +typedef Matx Matx31d; +typedef Matx Matx41f; +typedef Matx Matx41d; +typedef Matx Matx61f; +typedef Matx Matx61d; + +typedef Matx Matx22f; +typedef Matx Matx22d; +typedef Matx Matx23f; +typedef Matx Matx23d; +typedef Matx Matx32f; +typedef Matx Matx32d; + +typedef Matx Matx33f; +typedef Matx Matx33d; + +typedef Matx Matx34f; +typedef Matx Matx34d; +typedef Matx Matx43f; +typedef Matx Matx43d; + +typedef Matx Matx44f; +typedef Matx Matx44d; +typedef Matx Matx66f; +typedef Matx Matx66d; + //////////////////////////////// Complex ////////////////////////////// /*! @@ -767,7 +960,7 @@ public: is that each of them is basically a tuple of numbers of the same type. Each "scalar" can be represented by the depth id (CV_8U ... CV_64F) and the number of channels. OpenCV matrices, 2D or nD, dense or sparse, can store "scalars", - as long as the number of channels does not exceed CV_MAX_CN (currently set to 32). + as long as the number of channels does not exceed CV_CN_MAX. */ template class DataType { @@ -1046,45 +1239,13 @@ protected: //////////////////////////////// Mat //////////////////////////////// -class Mat; -class MatND; -template class CV_EXPORTS MatExpr_Base_; -typedef MatExpr_Base_ MatExpr_Base; -template class MatExpr_; -template class MatExpr_Op1_; -template class MatExpr_Op2_; -template class MatExpr_Op3_; -template class MatExpr_Op4_; -template class MatExpr_Op5_; -template class CV_EXPORTS MatOp_DivRS_; -template class CV_EXPORTS MatOp_Inv_; -template class CV_EXPORTS MatOp_MulDiv_; -template class CV_EXPORTS MatOp_Repeat_; -template class CV_EXPORTS MatOp_Set_; -template class CV_EXPORTS MatOp_Scale_; -template class CV_EXPORTS MatOp_T_; - typedef MatExpr_ >, Mat> MatExpr_Initializer; - -template class CV_EXPORTS MatIterator_; -template class CV_EXPORTS MatConstIterator_; -template class CV_EXPORTS MatCommaInitializer_; - + enum { MAGIC_MASK=0xFFFF0000, TYPE_MASK=0x00000FFF, DEPTH_MASK=7 }; static inline size_t getElemSize(int type) { return CV_ELEM_SIZE(type); } -// matrix decomposition types -enum { DECOMP_LU=0, DECOMP_SVD=1, DECOMP_EIG=2, DECOMP_CHOLESKY=3, DECOMP_QR=4, DECOMP_NORMAL=16 }; -enum { NORM_INF=1, NORM_L1=2, NORM_L2=4, NORM_TYPE_MASK=7, NORM_RELATIVE=8, NORM_MINMAX=32}; -enum { CMP_EQ=0, CMP_GT=1, CMP_GE=2, CMP_LT=3, CMP_LE=4, CMP_NE=5 }; -enum { GEMM_1_T=1, GEMM_2_T=2, GEMM_3_T=4 }; -enum { DFT_INVERSE=1, DFT_SCALE=2, DFT_ROWS=4, DFT_COMPLEX_OUTPUT=16, DFT_REAL_OUTPUT=32, - DCT_INVERSE = DFT_INVERSE, DCT_ROWS=DFT_ROWS }; - /*! The matrix class. @@ -1322,8 +1483,12 @@ public: Mat(const IplImage* img, bool copyData=false); //! builds matrix from std::vector with or without copying the data template explicit Mat(const vector<_Tp>& vec, bool copyData=false); - //! builds matrix from cv::Vec; the data is copied - template explicit Mat(const Vec<_Tp, n>& vec); + //! builds matrix from cv::Vec; the data is copied by default + template explicit Mat(const Vec<_Tp, n>& vec, + bool copyData=true); + //! builds matrix from cv::Matx; the data is copied by default + template explicit Mat(const Matx<_Tp, m, n>& mtx, + bool copyData=true); //! builds matrix from a 2D point template explicit Mat(const Point_<_Tp>& pt); //! builds matrix from a 3D point @@ -2116,7 +2281,8 @@ public: Mat_(const MatExpr_Base& expr); //! makes a matrix out of Vec, std::vector, Point_ or Point3_. The matrix will have a single column explicit Mat_(const vector<_Tp>& vec, bool copyData=false); - template explicit Mat_(const Vec<_Tp, n>& vec); + template explicit Mat_(const Vec<_Tp, n>& vec, bool copyData=true); + template explicit Mat_(const Matx<_Tp, m, n>& mtx, bool copyData=true); explicit Mat_(const Point_<_Tp>& pt); explicit Mat_(const Point3_<_Tp>& pt); explicit Mat_(const MatCommaInitializer_<_Tp>& commaInitializer); @@ -2359,20 +2525,29 @@ public: void assignTo(Mat& m, int type=-1) const; }; -#if 0 -template class VectorCommaInitializer_ + +template class CV_EXPORTS VecCommaInitializer { public: - VectorCommaInitializer_(vector<_Tp>* _vec); - template VectorCommaInitializer_<_Tp>& operator , (T2 val); - operator vector<_Tp>() const; - vector<_Tp> operator *() const; + VecCommaInitializer(Vec<_Tp, n>* _vec); + template VecCommaInitializer<_Tp, n>& operator , (T2 val); + Vec<_Tp, n> operator *() const; - vector<_Tp>* vec; + Vec<_Tp, n>* vec; int idx; }; -#endif - + + +template class CV_EXPORTS MatxCommaInitializer : + public VecCommaInitializer<_Tp, m*n> +{ +public: + MatxCommaInitializer(Matx<_Tp, m, n>* _mtx); + template MatxCommaInitializer<_Tp, m, n>& operator , (T2 val); + Matx<_Tp, m, n> operator *() const; +}; + + /*! Automatically Allocated Buffer Class @@ -3917,7 +4092,7 @@ public: SeqIterator<_Tp> end() const; //! returns the number of elements in the sequence size_t size() const; - //! returns the type of sequence elements (CV_8UC1 ... CV_64FC(CV_MAX_CN) ...) + //! returns the type of sequence elements (CV_8UC1 ... CV_64FC(CV_CN_MAX) ...) int type() const; //! returns the depth of sequence elements (CV_8U ... CV_64F) int depth() const; diff --git a/modules/core/include/opencv2/core/internal.hpp b/modules/core/include/opencv2/core/internal.hpp index 5e83efd..5f4d164 100644 --- a/modules/core/include/opencv2/core/internal.hpp +++ b/modules/core/include/opencv2/core/internal.hpp @@ -663,7 +663,7 @@ CvFuncTable; typedef struct CvBigFuncTable { - void* fn_2d[CV_DEPTH_MAX*CV_CN_MAX]; + void* fn_2d[CV_DEPTH_MAX*4]; } CvBigFuncTable; diff --git a/modules/core/include/opencv2/core/mat.hpp b/modules/core/include/opencv2/core/mat.hpp index 8691ab4..e5635b4 100644 --- a/modules/core/include/opencv2/core/mat.hpp +++ b/modules/core/include/opencv2/core/mat.hpp @@ -231,13 +231,39 @@ template inline Mat::Mat(const vector<_Tp>& vec, bool copyData) } -template inline Mat::Mat(const Vec<_Tp, n>& vec) +template inline Mat::Mat(const Vec<_Tp, n>& vec, bool copyData) : flags(MAGIC_VAL | DataType<_Tp>::type | CV_MAT_CONT_FLAG), rows(0), cols(0), step(0), data(0), refcount(0), datastart(0), dataend(0) { - create(n, 1, DataType<_Tp>::type); - memcpy(data, &vec[0], n*sizeof(_Tp)); + if( !copyData ) + { + rows = n; + cols = 1; + step = sizeof(_Tp); + data = datastart = (uchar*)vec.val; + dataend = datastart + rows*step; + } + else + Mat(n, 1, DataType<_Tp>::type, vec.val).copyTo(*this); +} + + +template inline Mat::Mat(const Matx<_Tp,m,n>& M, bool copyData) + : flags(MAGIC_VAL | DataType<_Tp>::type | CV_MAT_CONT_FLAG), + rows(0), cols(0), step(0), data(0), refcount(0), + datastart(0), dataend(0) +{ + if( !copyData ) + { + rows = m; + cols = n; + step = sizeof(_Tp); + data = datastart = (uchar*)M.val; + dataend = datastart + rows*step; + } + else + Mat(m, n, DataType<_Tp>::type, (uchar*)M.val).copyTo(*this); } @@ -619,13 +645,17 @@ template inline Mat_<_Tp>::Mat_(const Mat_& m, const Range& rowRan template inline Mat_<_Tp>::Mat_(const Mat_& m, const Rect& roi) : Mat(m, roi) {} -template template inline Mat_<_Tp>::Mat_(const Vec<_Tp, n>& vec) - : Mat(n, 1, DataType<_Tp>::type) +template template inline + Mat_<_Tp>::Mat_(const Vec<_Tp, n>& vec, bool copyData) + : Mat(vec, copyData) { - _Tp* d = (_Tp*)data; - for( int i = 0; i < n; i++ ) - d[i] = vec[i]; } + +template template inline + Mat_<_Tp>::Mat_(const Matx<_Tp,m,n>& M, bool copyData) + : Mat(M, copyData) +{ +} template inline Mat_<_Tp>::Mat_(const Point_<_Tp>& pt) : Mat(2, 1, DataType<_Tp>::type) diff --git a/modules/core/include/opencv2/core/operations.hpp b/modules/core/include/opencv2/core/operations.hpp index b0b56bb..e762756 100644 --- a/modules/core/include/opencv2/core/operations.hpp +++ b/modules/core/include/opencv2/core/operations.hpp @@ -101,6 +101,8 @@ #endif +#include + namespace cv { @@ -284,6 +286,12 @@ template inline Vec<_Tp, cn>::Vec(_Tp v0, _Tp v1, _Tp v2, for(int i = 10; i < cn; i++) val[i] = _Tp(0); } + +template inline Vec<_Tp, cn>::Vec(const _Tp* values) +{ + for( int i = 0; i < cn; i++ ) val[i] = values[i]; +} + template inline Vec<_Tp, cn>::Vec(const Vec<_Tp, cn>& v) { @@ -304,6 +312,7 @@ template inline _Tp Vec<_Tp, cn>::dot(const Vec<_Tp, cn>& return s; } + template inline double Vec<_Tp, cn>::ddot(const Vec<_Tp, cn>& v) const { double s = 0; @@ -311,11 +320,20 @@ template inline double Vec<_Tp, cn>::ddot(const Vec<_Tp, c return s; } + template inline Vec<_Tp, cn> Vec<_Tp, cn>::cross(const Vec<_Tp, cn>& v) const { - return Vec<_Tp, cn>(); // for arbitrary-size vector there is no cross-product defined + CV_Error(CV_StsError, "for arbitrary-size vector there is no cross-product defined"); + return Vec<_Tp, cn>(); +} + + +template inline Matx<_Tp, 1, cn> Vec<_Tp, cn>::t() const +{ + return (const Matx<_Tp, 1, cn>&)*this; } + template template inline Vec<_Tp, cn>::operator Vec() const { @@ -333,9 +351,30 @@ template inline Vec<_Tp, cn>::operator CvScalar() const return s; } -template inline const _Tp& Vec<_Tp, cn>::operator [](int i) const { return val[i]; } -template inline _Tp& Vec<_Tp, cn>::operator[](int i) { return val[i]; } +template inline const _Tp& Vec<_Tp, cn>::operator [](int i) const +{ + CV_DbgAssert( (unsigned)i < (unsigned)cn ); + return val[i]; +} + +template inline _Tp& Vec<_Tp, cn>::operator [](int i) +{ + CV_DbgAssert( (unsigned)i < (unsigned)cn ); + return val[i]; +} + +template inline const _Tp& Vec<_Tp, cn>::operator ()(int i) const +{ + CV_DbgAssert( (unsigned)i < (unsigned)cn ); + return val[i]; +} +template inline _Tp& Vec<_Tp, cn>::operator ()(int i) +{ + CV_DbgAssert( (unsigned)i < (unsigned)cn ); + return val[i]; +} + template static inline Vec<_Tp1, cn>& operator += (Vec<_Tp1, cn>& a, const Vec<_Tp2, cn>& b) { @@ -439,6 +478,7 @@ Vec& operator += (Vec& a, const Vec& b) return a; } + template static inline Vec& operator += (Vec& a, const Vec& b) { @@ -449,6 +489,7 @@ Vec& operator += (Vec& a, const Vec& b) return a; } + template static inline double norm(const Vec& a) { @@ -457,6 +498,31 @@ double norm(const Vec& a) s += (double)a.val[i]*a.val[i]; return std::sqrt(s); } + + +template static inline +double norm(const Vec& a, int normType) +{ + if( normType == NORM_INF ) + { + T1 s = 0; + for( int i = 0; i < n; i++ ) + s = std::max(s, std::abs(a.val[i])); + return s; + } + + if( normType == NORM_L1 ) + { + T1 s = 0; + for( int i = 0; i < n; i++ ) + s += std::abs(a.val[i]); + return s; + } + + CV_DbgAssert( normType == NORM_L2 ); + return norm(a); +} + template static inline bool operator == (const Vec& a, const Vec& b) @@ -471,6 +537,698 @@ bool operator != (const Vec& a, const Vec& b) { return !(a == b); } + +template static inline +VecCommaInitializer<_Tp, n> operator << (const Vec<_Tp, n>& vec, _T2 val) +{ + VecCommaInitializer<_Tp, n> commaInitializer((Vec<_Tp, n>*)&vec); + return (commaInitializer, val); +} + +template inline +VecCommaInitializer<_Tp, n>::VecCommaInitializer(Vec<_Tp, n>* _vec) + : vec(_vec), idx(0) +{} + +template template inline +VecCommaInitializer<_Tp, n>& VecCommaInitializer<_Tp, n>::operator , (_T2 value) +{ + CV_DbgAssert( idx < n ); + vec->val[idx++] = saturate_cast<_Tp>(value); + return *this; +} + +template inline +Vec<_Tp, n> VecCommaInitializer<_Tp, n>::operator *() const +{ + CV_DbgAssert( idx == n ); + return *vec; +} + +//////////////////////////////// Matx ///////////////////////////////// + + +template Matx<_Tp,m,n>::Matx() +{} + +template +inline Matx<_Tp,m,n>::Matx(_Tp v0) +: Matx<_Tp,m,n>::base_type(v0) +{} + +template +inline Matx<_Tp,m,n>::Matx(_Tp v0, _Tp v1) +: Matx<_Tp,m,n>::base_type(v0, v1) +{} + +template +inline Matx<_Tp,m,n>::Matx(_Tp v0, _Tp v1, _Tp v2) +: Matx<_Tp,m,n>::base_type(v0, v1, v2) +{} + +template Matx<_Tp,m,n>::Matx(_Tp v0, _Tp v1, _Tp v2, _Tp v3) +: Matx<_Tp,m,n>::base_type(v0, v1, v2, v3) +{} + +template +inline Matx<_Tp,m,n>::Matx(_Tp v0, _Tp v1, _Tp v2, _Tp v3, _Tp v4) +: Matx<_Tp,m,n>::base_type(v0, v1, v2, v3, v4) +{} + +template +inline Matx<_Tp,m,n>::Matx(_Tp v0, _Tp v1, _Tp v2, _Tp v3, _Tp v4, _Tp v5) +: Matx<_Tp,m,n>::base_type(v0, v1, v2, v3, v4, v5) +{} + +template +inline Matx<_Tp,m,n>::Matx(_Tp v0, _Tp v1, _Tp v2, _Tp v3, _Tp v4, _Tp v5, _Tp v6) +: Matx<_Tp,m,n>::base_type(v0, v1, v2, v3, v4, v5, v6) +{} + +template +inline Matx<_Tp,m,n>::Matx(_Tp v0, _Tp v1, _Tp v2, _Tp v3, + _Tp v4, _Tp v5, _Tp v6, _Tp v7) +: Matx<_Tp,m,n>::base_type(v0, v1, v2, v3, v4, v5, v6, v7) +{} + +template +inline Matx<_Tp,m,n>::Matx(_Tp v0, _Tp v1, _Tp v2, + _Tp v3, _Tp v4, _Tp v5, + _Tp v6, _Tp v7, _Tp v8) +: Matx<_Tp,m,n>::base_type(v0, v1, v2, v3, v4, v5, v6, v7, v8) +{} + +template +inline Matx<_Tp,m,n>::Matx(_Tp v0, _Tp v1, _Tp v2, _Tp v3, _Tp v4, + _Tp v5, _Tp v6, _Tp v7, _Tp v8, _Tp v9) +: Matx<_Tp,m,n>::base_type(v0, v1, v2, v3, v4, v5, v6, v7, v8, v9) +{} + +template +inline Matx<_Tp,m,n>::Matx(_Tp v0, _Tp v1, _Tp v2, _Tp v3, + _Tp v4, _Tp v5, _Tp v6, _Tp v7, + _Tp v8, _Tp v9, _Tp v10, _Tp v11) +{ + assert(m*n == 12); + this->val[0] = v0; this->val[1] = v1; this->val[2] = v2; this->val[3] = v3; + this->val[4] = v4; this->val[5] = v5; this->val[6] = v6; this->val[7] = v7; + this->val[8] = v8; this->val[9] = v9; this->val[10] = v10; this->val[11] = v11; +} + +template +inline Matx<_Tp,m,n>::Matx(_Tp v0, _Tp v1, _Tp v2, _Tp v3, + _Tp v4, _Tp v5, _Tp v6, _Tp v7, + _Tp v8, _Tp v9, _Tp v10, _Tp v11, + _Tp v12, _Tp v13, _Tp v14, _Tp v15) +{ + assert(m*n == 16); + this->val[0] = v0; this->val[1] = v1; this->val[2] = v2; this->val[3] = v3; + this->val[4] = v4; this->val[5] = v5; this->val[6] = v6; this->val[7] = v7; + this->val[8] = v8; this->val[9] = v9; this->val[10] = v10; this->val[11] = v11; + this->val[12] = v12; this->val[13] = v13; this->val[14] = v14; this->val[15] = v15; +} + + +template +inline Matx<_Tp,m,n>::Matx(const _Tp* vals) +: Matx<_Tp,m,n>::base_type(vals) +{} + + +template + inline Matx<_Tp,m,n>::Matx(const Matx<_Tp,m,n>::base_type& v) +: base_type(v) +{ +} + +template inline +Matx<_Tp,m,n> Matx<_Tp,m,n>::all(_Tp alpha) +{ + base_type v = base_type::all(alpha); + return (Matx<_Tp,m,n>&)v; +} + +template inline +Matx<_Tp,m,n> Matx<_Tp,m,n>::zeros() +{ + return all(0); +} + +template inline +Matx<_Tp,m,n> Matx<_Tp,m,n>::ones() +{ + return all(1); +} + +template inline +Matx<_Tp,m,n> Matx<_Tp,m,n>::eye() +{ + Matx<_Tp,m,n> M; + for(int i = 0; i < MIN(m,n); i++) + M(i,i) = 1; + return M; +} + + +template inline +Matx<_Tp,m,n> Matx<_Tp,m,n>::diag(const Vec<_Tp,MIN(m,n)>& d) +{ + Matx<_Tp,m,n> M; + for(int i = 0; i < MIN(m,n); i++) + M(i,i) = d[i]; +} + +template inline +Matx<_Tp,m,n> Matx<_Tp,m,n>::randu(_Tp a, _Tp b) +{ + Matx<_Tp,m,n> M; + Mat matM(M, false); + cv::randu(matM, Scalar(a), Scalar(b)); + return M; +} + +template inline +Matx<_Tp,m,n> Matx<_Tp,m,n>::randn(_Tp a, _Tp b) +{ + Matx<_Tp,m,n> M; + Mat matM(M, false); + cv::randn(matM, Scalar(a), Scalar(b)); + return M; +} + +template template +inline Matx<_Tp, m, n>::operator Matx() const +{ + Vec v = *this; + return (Matx&)v; +} + + +template template inline +Matx<_Tp, m1, n1> Matx<_Tp, m, n>::reshape() const +{ + CV_DbgAssert(m1*n1 == m*n); + return (const Matx<_Tp, m1, n1>&)*this; +} + + +template +template inline +Matx<_Tp, m1, n1> Matx<_Tp, m, n>::minor(int i, int j) const +{ + CV_DbgAssert(0 <= i && i+m1 <= m && 0 <= j && j+n1 <= n); + Matx<_Tp, m1, n1> s; + for( int di = 0; di < m1; di++ ) + for( int dj = 0; dj < n1; dj++ ) + s(di, dj) = (*this)(i+di, j+dj); + return s; +} + + +template inline +Matx<_Tp, 1, n> Matx<_Tp, m, n>::row(int i) const +{ + CV_DbgAssert((unsigned)i < (unsigned)m); + return (Matx<_Tp, 1, n>&)(*this)(i,0); +} + + +template inline +Vec<_Tp, m> Matx<_Tp, m, n>::col(int j) const +{ + CV_DbgAssert((unsigned)j < (unsigned)n); + Vec<_Tp, m> v; + for( int i = 0; i < m; i++ ) + v[i] = (*this)(i,j); + return v; +} + + +template inline +Vec<_Tp, MIN(m,n)> Matx<_Tp, m, n>::diag() const +{ + diag_type d; + for( int i = 0; i < MIN(m, n); i++ ) + d[i] = (*this)(i,i); + return d; +} + + +template inline +const _Tp& Matx<_Tp, m, n>::operator ()(int i, int j) const +{ + CV_DbgAssert( (unsigned)i < (unsigned)m && (unsigned)j < (unsigned)n ); + return this->val[i*n + j]; +} + + +template inline +_Tp& Matx<_Tp, m, n>::operator ()(int i, int j) +{ + CV_DbgAssert( (unsigned)i < (unsigned)m && (unsigned)j < (unsigned)n ); + return this->val[i*n + j]; +} + + +template inline +const _Tp& Matx<_Tp, m, n>::operator ()(int i) const +{ + CV_DbgAssert( (m == 1 || n == 1) && (unsigned)i < (unsigned)(m+n-1) ); + return this->val[i]; +} + + +template inline +_Tp& Matx<_Tp, m, n>::operator ()(int i) +{ + CV_DbgAssert( (m == 1 || n == 1) && (unsigned)i < (unsigned)(m+n-1) ); + return this->val[i]; +} + + +template static inline +Matx<_Tp1, m, n>& operator += (Matx<_Tp1, m, n>& a, const Matx<_Tp2, m, n>& b) +{ + for( int i = 0; i < m*n; i++ ) + a.val[i] = saturate_cast<_Tp1>(a.val[i] + b.val[i]); + return a; +} + + +template static inline +Matx<_Tp1, m, n>& operator -= (Matx<_Tp1, m, n>& a, const Matx<_Tp2, m, n>& b) +{ + for( int i = 0; i < m*n; i++ ) + a.val[i] = saturate_cast<_Tp1>(a.val[i] - b.val[i]); + return a; +} + + +template inline +Matx<_Tp,m,n>::Matx(const Matx<_Tp, m, n>& a, const Matx<_Tp, m, n>& b, Matx_AddOp) +{ + for( int i = 0; i < m*n; i++ ) + this->val[i] = saturate_cast<_Tp>(a.val[i] + b.val[i]); +} + + +template inline +Matx<_Tp,m,n>::Matx(const Matx<_Tp, m, n>& a, const Matx<_Tp, m, n>& b, Matx_SubOp) +{ + for( int i = 0; i < m*n; i++ ) + this->val[i] = saturate_cast<_Tp>(a.val[i] - b.val[i]); +} + + +template template inline +Matx<_Tp,m,n>::Matx(const Matx<_Tp, m, n>& a, _T2 alpha, Matx_ScaleOp) +{ + for( int i = 0; i < m*n; i++ ) + this->val[i] = saturate_cast<_Tp>(a.val[i] * alpha); +} + + +template inline +Matx<_Tp,m,n>::Matx(const Matx<_Tp, m, n>& a, const Matx<_Tp, m, n>& b, Matx_MulOp) +{ + for( int i = 0; i < m*n; i++ ) + this->val[i] = saturate_cast<_Tp>(a.val[i] * b.val[i]); +} + + +template template inline +Matx<_Tp,m,n>::Matx(const Matx<_Tp, m, l>& a, const Matx<_Tp, l, n>& b, Matx_MatMulOp) +{ + for( int i = 0; i < m; i++ ) + for( int j = 0; j < n; j++ ) + { + _Tp s = 0; + for( int k = 0; k < l; k++ ) + s += a(i, k) * b(k, j); + this->val[i*n + j] = s; + } +} + + +template inline +Matx<_Tp,m,n>::Matx(const Matx<_Tp, n, m>& a, Matx_TOp) +{ + for( int i = 0; i < m; i++ ) + for( int j = 0; j < n; j++ ) + this->val[i*n + j] = a(j, i); +} + + +template static inline +Matx<_Tp, m, n> operator + (const Matx<_Tp, m, n>& a, const Matx<_Tp, m, n>& b) +{ + return Matx<_Tp, m, n>(a, b, Matx_AddOp()); +} + + +template static inline +Matx<_Tp, m, n> operator - (const Matx<_Tp, m, n>& a, const Matx<_Tp, m, n>& b) +{ + return Matx<_Tp, m, n>(a, b, Matx_SubOp()); +} + + +template static inline +Matx<_Tp, m, n>& operator *= (Matx<_Tp, m, n>& a, int alpha) +{ + for( int i = 0; i < m*n; i++ ) + a.val[i] = saturate_cast<_Tp>(a.val[i] * alpha); +} + +template static inline +Matx<_Tp, m, n>& operator *= (Matx<_Tp, m, n>& a, float alpha) +{ + for( int i = 0; i < m*n; i++ ) + a.val[i] = saturate_cast<_Tp>(a.val[i] * alpha); +} + +template static inline +Matx<_Tp, m, n>& operator *= (Matx<_Tp, m, n>& a, double alpha) +{ + for( int i = 0; i < m*n; i++ ) + a.val[i] = saturate_cast<_Tp>(a.val[i] * alpha); +} + +template static inline +Matx<_Tp, m, n> operator * (const Matx<_Tp, m, n>& a, int alpha) +{ + return Matx<_Tp, m, n>(a, alpha, Matx_ScaleOp()); +} + +template static inline +Matx<_Tp, m, n> operator * (const Matx<_Tp, m, n>& a, float alpha) +{ + return Matx<_Tp, m, n>(a, alpha, Matx_ScaleOp()); +} + +template static inline +Matx<_Tp, m, n> operator * (const Matx<_Tp, m, n>& a, double alpha) +{ + return Matx<_Tp, m, n>(a, alpha, Matx_ScaleOp()); +} + +template static inline +Matx<_Tp, m, n> operator * (int alpha, const Matx<_Tp, m, n>& a) +{ + return Matx<_Tp, m, n>(a, alpha, Matx_ScaleOp()); +} + +template static inline +Matx<_Tp, m, n> operator * (float alpha, const Matx<_Tp, m, n>& a) +{ + return Matx<_Tp, m, n>(a, alpha, Matx_ScaleOp()); +} + +template static inline +Matx<_Tp, m, n> operator * (double alpha, const Matx<_Tp, m, n>& a) +{ + return Matx<_Tp, m, n>(a, alpha, Matx_ScaleOp()); +} + +template static inline +Matx<_Tp, m, n> operator - (const Matx<_Tp, m, n>& a) +{ + return Matx<_Tp, m, n>(a, -1, Matx_ScaleOp()); +} + + +template static inline +Matx<_Tp, m, n> operator * (const Matx<_Tp, m, l>& a, const Matx<_Tp, l, n>& b) +{ + return Matx<_Tp, m, n>(a, b, Matx_MatMulOp()); +} + + +template static inline +Vec<_Tp, m> operator * (const Matx<_Tp, m, n>& a, const Vec<_Tp, n>& b) +{ + return Matx<_Tp, m, 1>(a, (const Matx<_Tp, n, 1>&)b, Matx_MatMulOp()); +} + + +template inline +Matx<_Tp, m, n> Matx<_Tp, m, n>::mul(const Matx<_Tp, m, n>& a) const +{ + return Matx<_Tp, m, n>(*this, a, Matx_MulOp()); +} + + +CV_EXPORTS int LU(float* A, int m, float* b, int n); +CV_EXPORTS int LU(double* A, int m, double* b, int n); +CV_EXPORTS bool Cholesky(float* A, int m, float* b, int n); +CV_EXPORTS bool Cholesky(double* A, int m, double* b, int n); + + +template struct CV_EXPORTS Matx_DetOp +{ + double operator ()(const Matx<_Tp, m, m>& a) const + { + Matx<_Tp, m, m> temp = a; + double p = LU(temp.val, m, 0, 0); + if( p == 0 ) + return p; + for( int i = 0; i < m; i++ ) + p *= temp(i, i); + return p; + } +}; + + +template struct CV_EXPORTS Matx_DetOp<_Tp, 1> +{ + double operator ()(const Matx<_Tp, 1, 1>& a) const + { + return a(0,0); + } +}; + + +template struct CV_EXPORTS Matx_DetOp<_Tp, 2> +{ + double operator ()(const Matx<_Tp, 2, 2>& a) const + { + return a(0,0)*a(1,1) - a(0,1)*a(1,0); + } +}; + + +template struct CV_EXPORTS Matx_DetOp<_Tp, 3> +{ + double operator ()(const Matx<_Tp, 3, 3>& a) const + { + return a(0,0)*(a(1,1)*a(2,2) - a(2,1)*a(1,2)) - + a(0,1)*(a(1,0)*a(2,2) - a(2,0)*a(1,2)) + + a(0,2)*(a(1,0)*a(2,1) - a(2,0)*a(1,1)); + } +}; + +template static inline +double determinant(const Matx<_Tp, m, m>& a) +{ + return Matx_DetOp<_Tp, m>()(a); +} + + +template static inline +double trace(const Matx<_Tp, m, n>& a) +{ + _Tp s = 0; + for( int i = 0; i < std::min(m, n); i++ ) + s += a(i,i); + return s; +} + + +template inline +Matx<_Tp, n, m> Matx<_Tp, m, n>::t() const +{ + return Matx<_Tp, n, m>(*this, Matx_TOp()); +} + + +template struct CV_EXPORTS Matx_FastInvOp +{ + bool operator()(const Matx<_Tp, m, m>& a, Matx<_Tp, m, m>& b, int method) const + { + Matx<_Tp, m, m> temp = a; + + // assume that b is all 0's on input => make it a unity matrix + for( int i = 0; i < m; i++ ) + b(i, i) = (_Tp)1; + + if( method == DECOMP_CHOLESKY ) + return Cholesky(temp.val, m, b.val, m); + + return LU(temp.val, m, b.val, m) != 0; + } +}; + + +template struct CV_EXPORTS Matx_FastInvOp<_Tp, 2> +{ + bool operator()(const Matx<_Tp, 2, 2>& a, Matx<_Tp, 2, 2>& b, int) const + { + _Tp d = determinant(a); + if( d == 0 ) + return false; + d = 1/d; + b(1,1) = a(0,0)*d; + b(0,0) = a(1,1)*d; + b(0,1) = -a(0,1)*d; + b(1,0) = -a(1,0)*d; + return true; + } +}; + + +template struct CV_EXPORTS Matx_FastInvOp<_Tp, 3> +{ + bool operator()(const Matx<_Tp, 3, 3>& a, Matx<_Tp, 3, 3>& b, int) const + { + _Tp d = determinant(a); + if( d == 0 ) + return false; + d = 1/d; + b(0,0) = (a(1,1) * a(2,2) - a(1,2) * a(2,1)) * d; + b(0,1) = (a(0,2) * a(2,1) - a(0,1) * a(2,2)) * d; + b(0,2) = (a(0,1) * a(1,2) - a(0,2) * a(1,1)) * d; + + b(1,0) = (a(1,2) * a(2,0) - a(1,0) * a(2,2)) * d; + b(1,1) = (a(0,0) * a(2,2) - a(0,2) * a(2,0)) * d; + b(1,2) = (a(0,2) * a(1,0) - a(0,0) * a(1,2)) * d; + + b(2,0) = (a(1,0) * a(2,1) - a(1,1) * a(2,0)) * d; + b(2,1) = (a(0,1) * a(2,0) - a(0,0) * a(2,1)) * d; + b(2,2) = (a(0,0) * a(1,1) - a(0,1) * a(1,0)) * d; + return true; + } +}; + + +template inline +Matx<_Tp, n, m> Matx<_Tp, m, n>::inv(int method) const +{ + Matx<_Tp, n, m> b; + bool ok; + if( method == DECOMP_LU || method == DECOMP_CHOLESKY ) + ok = Matx_FastInvOp<_Tp, m>()(*this, b, method); + else + { + Mat A(*this, false), B(b, false); + ok = invert(A, B, method); + } + return ok ? b : Matx<_Tp, n, m>::zeros(); +} + + +template struct CV_EXPORTS Matx_FastSolveOp +{ + 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, m> temp = a; + x = b; + if( method == DECOMP_CHOLESKY ) + return Cholesky(temp.val, m, x.val, n); + + return LU(temp.val, m, x.val, n) != 0; + } +}; + + +template struct CV_EXPORTS Matx_FastSolveOp<_Tp, 2, 1> +{ + bool operator()(const Matx<_Tp, 2, 2>& a, const Matx<_Tp, 2, 1>& b, + Matx<_Tp, 2, 1>& x, int method) const + { + _Tp d = determinant(a); + if( d == 0 ) + return false; + d = 1/d; + x(0) = (b(0)*a(1,1) - b(1)*a(0,1))*d; + x(1) = (b(1)*a(0,0) - b(0)*a(1,0))*d; + return true; + } +}; + + +template struct CV_EXPORTS Matx_FastSolveOp<_Tp, 3, 1> +{ + bool operator()(const Matx<_Tp, 3, 3>& a, const Matx<_Tp, 3, 1>& b, + Matx<_Tp, 3, 1>& x, int method) const + { + _Tp d = determinant(a); + if( d == 0 ) + return false; + d = 1/d; + x(0) = d*(b(0)*(a(1,1)*a(2,2) - a(1,2)*a(2,1)) - + a(0,1)*(b(1)*a(2,2) - a(1,2)*b(2)) + + a(0,2)*(b(1)*a(2,1) - a(1,1)*b(2))); + + x(1) = d*(a(0,0)*(b(1)*a(2,2) - a(1,2)*b(2)) - + b(0)*(a(1,0)*a(2,2) - a(1,2)*a(2,0)) + + a(0,2)*(a(1,0)*b(2) - b(1)*a(2,0))); + + x(2) = d*(a(0,0)*(a(1,1)*b(2) - b(1)*a(2,1)) - + a(0,1)*(a(1,0)*b(2) - b(1)*a(2,0)) + + b(0)*(a(1,0)*a(2,1) - a(1,1)*a(2,0))); + return true; + } +}; + + +template template inline +Matx<_Tp, n, l> Matx<_Tp, m, n>::solve(const Matx<_Tp, m, l>& rhs, int method) const +{ + Matx<_Tp, n, l> x; + bool ok; + if( method == DECOMP_LU || method == DECOMP_CHOLESKY ) + ok = Matx_FastSolveOp<_Tp, m, 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(); +} + + +template inline +Vec<_Tp, n> Matx<_Tp, m, n>::solve(const Vec<_Tp, m>& rhs, int method) const +{ + return solve(Matx<_Tp, m, 1>(rhs), method); +} + +template static inline +MatxCommaInitializer<_Tp, m, n> operator << (const Matx<_Tp, m, n>& mtx, _T2 val) +{ + MatxCommaInitializer<_Tp, m, n> commaInitializer((Matx<_Tp, m, n>*)&mtx); + return (commaInitializer, val); +} + +template inline +MatxCommaInitializer<_Tp, m, n>::MatxCommaInitializer(Matx<_Tp, m, n>* _mtx) + : VecCommaInitializer<_Tp, m*n>((Vec<_Tp,m*n>*)_mtx) +{} + +template template inline +MatxCommaInitializer<_Tp, m, n>& MatxCommaInitializer<_Tp, m, n>::operator , (_T2 value) +{ + return (MatxCommaInitializer<_Tp, m, n>&)VecCommaInitializer<_Tp, m*n>::operator ,(value); +} + +template inline +Matx<_Tp, m, n> MatxCommaInitializer<_Tp, m, n>::operator *() const +{ + CV_DbgAssert( this->idx == n*m ); + return (Matx<_Tp, m, n>&)*(this->vec); +} //////////////////////////////// Complex ////////////////////////////// @@ -1481,35 +2239,6 @@ inline Point LineIterator::pos() const p.x = ((ptr - ptr0) - p.y*step)/elemSize; return p; } - -#if 0 - template inline VectorCommaInitializer_<_Tp>:: - VectorCommaInitializer_(vector<_Tp>* _vec) : vec(_vec), idx(0) {} - - template template inline VectorCommaInitializer_<_Tp>& - VectorCommaInitializer_<_Tp>::operator , (T2 val) - { - if( (size_t)idx < vec->size() ) - (*vec)[idx] = _Tp(val); - else - vec->push_back(_Tp(val)); - idx++; - return *this; - } - - template inline VectorCommaInitializer_<_Tp>::operator vector<_Tp>() const - { return *vec; } - - template inline vector<_Tp> VectorCommaInitializer_<_Tp>::operator *() const - { return *vec; } - - template static inline VectorCommaInitializer_<_Tp> - operator << (const vector<_Tp>& vec, T2 val) - { - VectorCommaInitializer_<_Tp> commaInitializer((vector<_Tp>*)&vec); - return (commaInitializer, val); - } -#endif /////////////////////////////// AutoBuffer //////////////////////////////////////// diff --git a/modules/core/include/opencv2/core/types_c.h b/modules/core/include/opencv2/core/types_c.h index 345001c..159a120 100644 --- a/modules/core/include/opencv2/core/types_c.h +++ b/modules/core/include/opencv2/core/types_c.h @@ -518,7 +518,7 @@ IplConvKernelFP; * Matrix type (CvMat) * \****************************************************************************************/ -#define CV_CN_MAX 64 +#define CV_CN_MAX 1024 #define CV_CN_SHIFT 3 #define CV_DEPTH_MAX (1 << CV_CN_SHIFT) diff --git a/modules/core/src/lapack.cpp b/modules/core/src/lapack.cpp index 86e1671..134080b 100644 --- a/modules/core/src/lapack.cpp +++ b/modules/core/src/lapack.cpp @@ -57,7 +57,162 @@ namespace cv { + +/****************************************************************************************\ +* LU & Cholesky implementation for small matrices * +\****************************************************************************************/ + +template static inline int LUImpl(_Tp* A, int m, _Tp* b, int n) +{ + int i, j, k, p = 1; + + for( i = 0; i < m; i++ ) + { + k = i; + + for( j = i+1; j < m; j++ ) + if( std::abs(A[j*m + i]) > std::abs(A[k*m + i]) ) + k = j; + + if( std::abs(A[k*m + i]) < std::numeric_limits<_Tp>::epsilon() ) + return 0; + + if( k != i ) + { + for( j = i; j < m; j++ ) + std::swap(A[i*m + j], A[k*m + j]); + if( b ) + for( j = 0; j < n; j++ ) + std::swap(b[i*n + j], b[k*n + j]); + p = -p; + } + + _Tp d = -1/A[i*m + i]; + + for( j = i+1; j < m; j++ ) + { + _Tp alpha = A[j*m + i]*d; + + for( k = i+1; k < m; k++ ) + A[j*m + k] += alpha*A[i*m + k]; + + if( b ) + for( k = 0; k < n; k++ ) + b[j*n + k] += alpha*b[i*n + k]; + } + + A[i*m + i] = -d; + } + + if( b ) + { + for( i = m-1; i >= 0; i-- ) + for( j = 0; j < n; j++ ) + { + _Tp s = b[i*n + j]; + for( k = i+1; k < m; k++ ) + s -= A[i*m + k]*b[k*n + j]; + b[i*n + j] = s*A[i*m + i]; + } + } + + return p; +} + + +int LU(float* A, int m, float* b, int n) +{ + return LUImpl(A, m, b, n); +} + + +int LU(double* A, int m, double* b, int n) +{ + return LUImpl(A, m, b, n); +} + +template static inline bool CholImpl(_Tp* A, int m, _Tp* b, int n) +{ + _Tp* L = A; + int i, j, k; + double s; + + for( i = 0; i < m; i++ ) + { + for( j = 0; j < i; j++ ) + { + s = A[i*m + j]; + for( k = 0; k < j; k++ ) + s -= L[i*m + k]*L[j*m + k]; + L[i*m + j] = (_Tp)(s*L[j*m + j]); + } + s = A[i*m + i]; + for( k = 0; k < j; k++ ) + { + double t = L[i*m + k]; + s -= t*t; + } + if( s < std::numeric_limits<_Tp>::epsilon() ) + return 0; + L[i*m + i] = (_Tp)(1./std::sqrt(s)); + } + + if( !b ) + return false; + + // LLt x = b + // 1: L y = b + // 2. Lt x = y + + /* + [ L00 ] y0 b0 + [ L10 L11 ] y1 = b1 + [ L20 L21 L22 ] y2 b2 + [ L30 L31 L32 L33 ] y3 b3 + + [ L00 L10 L20 L30 ] x0 y0 + [ L11 L21 L31 ] x1 = y1 + [ L22 L32 ] x2 y2 + [ L33 ] x3 y3 + */ + + for( i = 0; i < m; i++ ) + { + for( j = 0; j < n; j++ ) + { + s = b[i*n + j]; + for( k = 0; k < i; k++ ) + s -= L[i*m + k]*b[k*n + j]; + b[i*n + j] = (_Tp)(s*L[i*m + i]); + } + } + + for( i = m-1; i >= 0; i-- ) + { + for( j = 0; j < n; j++ ) + { + s = b[i*n + j]; + for( k = m-1; k > i; k-- ) + s -= L[k*m + i]*b[k*n + j]; + b[i*n + j] = (_Tp)(s*L[i*m + i]); + } + } + + return true; +} + + +bool Cholesky(float* A, int m, float* b, int n) +{ + return CholImpl(A, m, b, n); +} + +bool Cholesky(double* A, int m, double* b, int n) +{ + return CholImpl(A, m, b, n); +} + /****************************************************************************************\ * Determinant of the matrix * \****************************************************************************************/ diff --git a/modules/core/src/matrix.cpp b/modules/core/src/matrix.cpp index c15361b..8c93b9e 100644 --- a/modules/core/src/matrix.cpp +++ b/modules/core/src/matrix.cpp @@ -2344,7 +2344,7 @@ SparseMat::Hdr::Hdr( int _dims, const int* _sizes, int _type ) dims = _dims; valueOffset = (int)alignSize(sizeof(SparseMat::Node) + - sizeof(int)*(dims - CV_MAX_DIM), CV_ELEM_SIZE1(_type)); + sizeof(int)*std::max(dims - CV_MAX_DIM, 0), CV_ELEM_SIZE1(_type)); nodeSize = alignSize(valueOffset + CV_ELEM_SIZE(_type), (int)sizeof(size_t));