TooN 2.1
|
00001 // -*- c++ -*- 00002 00003 // Copyright (C) 2005,2009 Tom Drummond (twd20@cam.ac.uk), 00004 // Ed Rosten (er258@cam.ac.uk), Gerhard Reitmayr (gr281@cam.ac.uk) 00005 // 00006 // This file is part of the TooN Library. This library is free 00007 // software; you can redistribute it and/or modify it under the 00008 // terms of the GNU General Public License as published by the 00009 // Free Software Foundation; either version 2, or (at your option) 00010 // any later version. 00011 00012 // This library is distributed in the hope that it will be useful, 00013 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 // GNU General Public License for more details. 00016 00017 // You should have received a copy of the GNU General Public License along 00018 // with this library; see the file COPYING. If not, write to the Free 00019 // Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, 00020 // USA. 00021 00022 // As a special exception, you may use this file as part of a free software 00023 // library without restriction. Specifically, if other files instantiate 00024 // templates or use macros or inline functions from this file, or you compile 00025 // this file and link it with other files to produce an executable, this 00026 // file does not by itself cause the resulting executable to be covered by 00027 // the GNU General Public License. This exception does not however 00028 // invalidate any other reasons why the executable file might be covered by 00029 // the GNU General Public License. 00030 00031 00032 #ifndef TOON_INCLUDE_HELPERS_H 00033 #define TOON_INCLUDE_HELPERS_H 00034 00035 #include <TooN/TooN.h> 00036 #include <TooN/gaussian_elimination.h> 00037 #include <cmath> 00038 #include <functional> 00039 #include <utility> 00040 00041 #ifndef M_PI 00042 #define M_PI 3.14159265358979323846 00043 #endif 00044 00045 #ifndef M_SQRT1_2 00046 #define M_SQRT1_2 0.707106781186547524401 00047 #endif 00048 00049 #ifdef WIN32 00050 namespace std { 00051 inline int isfinite( const double & v ){ return _finite(v); } 00052 inline int isfinite( const float & v ){ return _finite(v); } 00053 inline int isnan( const double & v ){ return _isnan(v); } 00054 inline int isnan( const float & v ){ return _isnan(v); } 00055 } 00056 #endif 00057 00058 namespace TooN { 00059 00060 ///Invert a matrix. 00061 /// 00062 ///For sizes other than 2x2, @link gDecomps decompositions @endlink provide a suitable solition. 00063 /// 00064 ///@param m Matrix to invert. 00065 ///@ingroup gDecomps 00066 inline Matrix<2> inv(const Matrix<2>& m) 00067 { 00068 double d = 1./(m[0][0]*m[1][1] - m[1][0]*m[0][1]); 00069 00070 return Data( 00071 m[1][1]*d, -m[0][1]*d, 00072 -m[1][0]*d, m[0][0]*d 00073 ); 00074 } 00075 00076 00077 00078 ///\deprecated 00079 ///@ingroup gLinAlg 00080 template<int Size, class Precision, class Base> TOON_DEPRECATED void Fill(Vector<Size, Precision, Base>& v, const Precision& p) 00081 { 00082 for(int i=0; i < v.size(); i++) 00083 v[i]= p; 00084 } 00085 00086 ///\deprecated 00087 ///@ingroup gLinAlg 00088 template<int Rows, int Cols, class Precision, class Base> TOON_DEPRECATED void Fill(Matrix<Rows, Cols, Precision, Base>& m, const Precision& p) 00089 { 00090 for(int i=0; i < m.num_rows(); i++) 00091 for(int j=0; j < m.num_cols(); j++) 00092 m[i][j] = p; 00093 } 00094 00095 ///Compute the \f$L_2\f$ norm of \e v 00096 ///@param v \e v 00097 ///@ingroup gLinAlg 00098 template<int Size, class Precision, class Base> inline Precision norm(const Vector<Size, Precision, Base>& v) 00099 { 00100 using std::sqrt; 00101 return sqrt(v*v); 00102 } 00103 00104 ///Compute the \f$L_2^2\f$ norm of \e v 00105 ///@param v \e v 00106 ///@ingroup gLinAlg 00107 template<int Size, class Precision, class Base> inline Precision norm_sq(const Vector<Size, Precision, Base>& v) 00108 { 00109 return v*v; 00110 } 00111 00112 ///Compute the \f$L_1\f$ norm of \e v 00113 ///@param v \e v 00114 ///@ingroup gLinAlg 00115 template<int Size, class Precision, class Base> inline Precision norm_1(const Vector<Size, Precision, Base>& v) 00116 { 00117 using std::abs; 00118 Precision n = 0; 00119 for(int i=0; i < v.size(); i++) 00120 n += abs(v[i]); 00121 return n; 00122 } 00123 00124 ///Compute the \f$L_\infty\f$ norm of \e v 00125 ///@param v \e v 00126 ///@ingroup gLinAlg 00127 template<int Size, class Precision, class Base> inline Precision norm_inf(const Vector<Size, Precision, Base>& v) 00128 { 00129 using std::abs; 00130 using std::max; 00131 Precision n = 0; 00132 n = abs(v[0]); 00133 00134 for(int i=1; i < v.size(); i++) 00135 n = max(n, abs(v[i])); 00136 return n; 00137 } 00138 00139 ///Compute the \f$L_2\f$ norm of \e v. 00140 ///Synonym for norm() 00141 ///@param v \e v 00142 ///@ingroup gLinAlg 00143 template<int Size, class Precision, class Base> inline Precision norm_2(const Vector<Size, Precision, Base>& v) 00144 { 00145 return norm(v); 00146 } 00147 00148 00149 00150 00151 ///Compute a the unit vector \f$\hat{v}\f$. 00152 ///@param v \e v 00153 ///@ingroup gLinAlg 00154 template<int Size, class Precision, class Base> inline Vector<Size, Precision> unit(const Vector<Size, Precision, Base> & v) 00155 { 00156 using std::sqrt; 00157 return TooN::operator*(v,(1/sqrt(v*v))); 00158 } 00159 00160 //Note because of the overload later, this function will ONLY receive sliced vectors. Therefore 00161 //a copy can be made, which is still a slice, so operating on the copy operates on the original 00162 //data. 00163 ///Normalize a vector in place 00164 ///@param v Vector to normalize 00165 ///@ingroup gLinAlg 00166 template<int Size, class Precision, class Base> inline void normalize(Vector<Size, Precision, Base> v) 00167 { 00168 using std::sqrt; 00169 v /= sqrt(v*v); 00170 } 00171 00172 //This overload is required to operate on non-slice vectors 00173 /** 00174 \overload 00175 */ 00176 template<int Size, class Precision> inline void normalize(Vector<Size, Precision> & v) 00177 { 00178 normalize(v.as_slice()); 00179 } 00180 00181 ///For a vector \e v of length \e i, return \f$[v_1, v_2, \cdots, v_{i-1}] / v_i \f$ 00182 ///@param v \e v 00183 ///@ingroup gLinAlg 00184 // Don't remove the +0 in the first template parameter of the return type. It is a work around for a Visual Studio 2010 bug: 00185 // https://connect.microsoft.com/VisualStudio/feedback/details/735283/vc-2010-parse-error-in-template-parameters-using-ternary-operator 00186 template<int Size, typename Precision, typename Base> inline Vector<(Size==Dynamic?Dynamic:Size-1)+0, Precision> project( const Vector<Size, Precision, Base> & v){ 00187 static const int Len = (Size==Dynamic?Dynamic:Size-1); 00188 return TooN::operator/(v.template slice<0, Len>(0 , v.size() - 1) , v[v.size() - 1]); 00189 } 00190 00191 //This should probably be done with an operator to prevent an extra new[] for dynamic vectors. 00192 ///For a vector \e v of length \e i, return \f$[v_1, v_2, \cdots, v_{i}, 1]\f$ 00193 ///@param v \e v 00194 ///@ingroup gLinAlg 00195 // Don't remove the +0 in the first template parameter of the return type. It is a work around for a Visual Studio 2010 bug: 00196 // https://connect.microsoft.com/VisualStudio/feedback/details/735283/vc-2010-parse-error-in-template-parameters-using-ternary-operator 00197 template<int Size, typename Precision, typename Base> inline Vector<(Size==Dynamic?Dynamic:Size+1)+0, Precision> unproject( const Vector<Size, Precision, Base> & v){ 00198 Vector<(Size==Dynamic?Dynamic:Size+1), Precision> result(v.size()+1); 00199 static const int Len = (Size==Dynamic?Dynamic:Size); 00200 result.template slice<0, Len>(0, v.size()) = v; 00201 result[v.size()] = 1; 00202 return result; 00203 } 00204 00205 /** 00206 \overload 00207 */ 00208 template<int R, int C, typename Precision, typename Base> inline Matrix<R-1, C, Precision> project( const Matrix<R,C, Precision, Base> & m){ 00209 Matrix<R-1, C, Precision> result = m.slice(0,0,R-1,m.num_cols()); 00210 for( int c = 0; c < m.num_cols(); ++c ) { 00211 result.slice(0,c,R-1,1) /= m[R-1][c]; 00212 } 00213 return result; 00214 } 00215 00216 template<int C, typename Precision, typename Base> inline Matrix<-1, C, Precision> project( const Matrix<-1,C, Precision, Base> & m){ 00217 Matrix<-1, C, Precision> result = m.slice(0,0,m.num_rows()-1,m.num_cols()); 00218 for( int c = 0; c < m.num_cols(); ++c ) { 00219 result.slice(0,c,m.num_rows()-1,1) /= m[m.num_rows()-1][c]; 00220 } 00221 return result; 00222 } 00223 00224 /** 00225 \overload 00226 */ 00227 template<int R, int C, typename Precision, typename Base> inline Matrix<R+1, C, Precision> unproject( const Matrix<R, C, Precision, Base> & m){ 00228 Matrix<R+1, C, Precision> result; 00229 result.template slice<0,0,R,C>() = m; 00230 result[R] = Ones; 00231 return result; 00232 } 00233 00234 template<int C, typename Precision, typename Base> inline Matrix<-1, C, Precision> unproject( const Matrix<-1, C, Precision, Base> & m){ 00235 Matrix<-1, C, Precision> result( m.num_rows()+1, m.num_cols() ); 00236 result.slice(0,0,m.num_rows(),m.num_cols()) = m; 00237 result[m.num_rows()] = Ones; 00238 return result; 00239 } 00240 00241 /// Frobenius (root of sum of squares) norm of input matrix \e m 00242 ///@param m \e m 00243 ///@ingroup gLinAlg 00244 template <int R, int C, typename P, typename B> 00245 P inline norm_fro( const Matrix<R,C,P,B> & m ){ 00246 using std::sqrt; 00247 P n = 0; 00248 for(int r = 0; r < m.num_rows(); ++r) 00249 for(int c = 0; c < m.num_cols(); ++c) 00250 n += m[r][c] * m[r][c]; 00251 00252 return sqrt(n); 00253 } 00254 00255 /// \e L<sub>∞</sub> (row sum) norm of input matrix m 00256 /// computes the maximum of the sums of absolute values over rows 00257 ///@ingroup gLinAlg 00258 template <int R, int C, typename P, typename B> 00259 P inline norm_inf( const Matrix<R,C,P,B> & m ){ 00260 using std::abs; 00261 using std::max; 00262 P n = 0; 00263 for(int r = 0; r < m.num_rows(); ++r){ 00264 P s = 0; 00265 for(int c = 0; c < m.num_cols(); ++c) 00266 s += abs(m(r,c)); 00267 n = max(n,s); 00268 } 00269 return n; 00270 } 00271 00272 /// \e L<sub>1</sub> (col sum) norm of input matrix m 00273 /// computes the maximum of the sums of absolute values over columns 00274 ///@ingroup gLinAlg 00275 template <int R, int C, typename P, typename B> 00276 P inline norm_1( const Matrix<R,C,P,B> & m ){ 00277 using std::abs; 00278 using std::max; 00279 P n = 0; 00280 for(int c = 0; c < m.num_cols(); ++c){ 00281 P s = 0; 00282 for(int r = 0; r < m.num_rows(); ++r) 00283 s += abs(m(r,c)); 00284 n = max(n,s); 00285 } 00286 return n; 00287 } 00288 00289 namespace Internal { 00290 ///@internal 00291 ///@brief Exponentiate a matrix using a the Taylor series 00292 ///This will not work if the norm of the matrix is too large. 00293 template <int R, int C, typename P, typename B> 00294 inline Matrix<R, C, P> exp_taylor( const Matrix<R,C,P,B> & m ){ 00295 TooN::SizeMismatch<R, C>::test(m.num_rows(), m.num_cols()); 00296 Matrix<R,C,P> result = TooN::Zeros(m.num_rows(), m.num_cols()); 00297 Matrix<R,C,P> f = TooN::Identity(m.num_rows()); 00298 P k = 1; 00299 while(norm_inf((result+f)-result) > 0){ 00300 result += f; 00301 f = (m * f) / k; 00302 k += 1; 00303 } 00304 return result; 00305 } 00306 00307 ///@internal 00308 ///@brief Logarithm of a matrix using a the Taylor series 00309 ///This will not work if the norm of the matrix is too large. 00310 template <int R, int C, typename P, typename B> 00311 inline Matrix<R, C, P> log_taylor( const Matrix<R,C,P,B> & m ){ 00312 TooN::SizeMismatch<R, C>::test(m.num_rows(), m.num_cols()); 00313 Matrix<R,C,P> X = m - TooN::Identity * 1.0; 00314 Matrix<R,C,P> F = X; 00315 Matrix<R,C,P> sum = TooN::Zeros(m.num_rows(), m.num_cols()); 00316 P k = 1; 00317 while(norm_inf((sum+F/k)-sum) > 0){ 00318 sum += F/k; 00319 F = -F*X; 00320 k += 1; 00321 } 00322 return sum; 00323 } 00324 00325 }; 00326 00327 /// computes the matrix exponential of a matrix m by 00328 /// scaling m by 1/(powers of 2), using Taylor series and 00329 /// squaring again. 00330 /// @param m input matrix, must be square 00331 /// @return result matrix of the same size/type as input 00332 /// @ingroup gLinAlg 00333 template <int R, int C, typename P, typename B> 00334 inline Matrix<R, C, P> exp( const Matrix<R,C,P,B> & m ){ 00335 using std::max; 00336 SizeMismatch<R, C>::test(m.num_rows(), m.num_cols()); 00337 const P l = log2(norm_inf(m)); 00338 const int s = max(0,(int)ceil(l)); 00339 Matrix<R,C,P> result = Internal::exp_taylor(m/(1<<s)); 00340 for(int i = 0; i < s; ++i) 00341 result = result * result; 00342 return result; 00343 } 00344 00345 /// computes a matrix square root of a matrix m by 00346 /// the product form of the Denman and Beavers iteration 00347 /// as given in Chen et al. 'Approximating the logarithm of a matrix to specified accuracy', 00348 /// J. Matrix Anal Appl, 2001. This is used for the matrix 00349 /// logarithm function, but is useable by on its own. 00350 /// @param m input matrix, must be square 00351 /// @return a square root of m of the same size/type as input 00352 /// @ingroup gLinAlg 00353 template <int R, int C, typename P, typename B> 00354 inline Matrix<R, C, P> sqrt( const Matrix<R,C,P,B> & m){ 00355 SizeMismatch<R, C>::test(m.num_rows(), m.num_cols()); 00356 Matrix<R,C,P> M = m; 00357 Matrix<R,C,P> Y = m; 00358 Matrix<R,C,P> M_inv(m.num_rows(), m.num_cols()); 00359 const Matrix<R,C,P> id = Identity(m.num_rows()); 00360 do { 00361 M_inv = gaussian_elimination(M, id); 00362 Y = Y * (id + M_inv) * 0.5; 00363 M = 0.5 * (id + (M + M_inv) * 0.5); 00364 } while(norm_inf(M - M_inv) > 0); 00365 return Y; 00366 } 00367 00368 /// computes the matrix logarithm of a matrix m using the inverse scaling and 00369 /// squaring method. The overall approach is described in 00370 /// Chen et al. 'Approximating the logarithm of a matrix to specified accuracy', 00371 /// J. Matrix Anal Appl, 2001, but this implementation only uses a simple 00372 /// taylor series after the repeated square root operation. 00373 /// @param m input matrix, must be square 00374 /// @return the log of m of the same size/type as input 00375 /// @ingroup gLinAlg 00376 template <int R, int C, typename P, typename B> 00377 inline Matrix<R, C, P> log( const Matrix<R,C,P,B> & m){ 00378 int counter = 0; 00379 Matrix<R,C,P> A = m; 00380 while(norm_inf(A - Identity*1.0) > 0.5){ 00381 ++counter; 00382 A = sqrt(A); 00383 } 00384 return Internal::log_taylor(A) * pow(2.0, counter); 00385 } 00386 00387 /// Returns true if every element is finite 00388 ///@ingroup gLinAlg 00389 template<int S, class P, class B> bool isfinite(const Vector<S, P, B>& v) 00390 { 00391 using std::isfinite; 00392 for(int i=0; i < v.size(); i++) 00393 if(!isfinite(v[i])) 00394 return 0; 00395 return 1; 00396 } 00397 00398 /// Returns true if any element is NaN 00399 ///@ingroup gLinAlg 00400 template<int S, class P, class B> bool isnan(const Vector<S, P, B>& v) 00401 { 00402 using std::isnan; 00403 for(int i=0; i < v.size(); i++) 00404 if(isnan(v[i])) 00405 return 1; 00406 return 0; 00407 } 00408 00409 /// Symmetrize a matrix 00410 ///@param m \e m 00411 ///@return \f$ \frac{m + m^{\mathsf T}}{2} \f$ 00412 ///@ingroup gLinAlg 00413 template<int Rows, int Cols, typename Precision, typename Base> 00414 void Symmetrize(Matrix<Rows,Cols,Precision,Base>& m){ 00415 SizeMismatch<Rows,Cols>::test(m.num_rows(), m.num_cols()); 00416 for(int r=0; r<m.num_rows()-1; r++){ 00417 for(int c=r+1; c<m.num_cols(); c++){ 00418 const Precision temp=(m(r,c)+m(c,r))/2; 00419 m(r,c)=temp; 00420 m(c,r)=temp; 00421 } 00422 } 00423 } 00424 00425 /// computes the trace of a square matrix 00426 ///@ingroup gLinAlg 00427 template<int Rows, int Cols, typename Precision, typename Base> 00428 Precision trace(const Matrix<Rows, Cols, Precision, Base> & m ){ 00429 SizeMismatch<Rows,Cols>::test(m.num_rows(), m.num_cols()); 00430 Precision tr = 0; 00431 for(int i = 0; i < m.num_rows(); ++i) 00432 tr += m(i,i); 00433 return tr; 00434 } 00435 00436 /// creates an returns a cross product matrix M from a 3 vector v, such that for all vectors w, the following holds: v ^ w = M * w 00437 /// @param vec the 3 vector input 00438 /// @return the 3x3 matrix to set to the cross product matrix 00439 ///@ingroup gLinAlg 00440 template<int Size, class P, class B> inline TooN::Matrix<3, 3, P> cross_product_matrix(const Vector<Size, P, B>& vec) 00441 { 00442 SizeMismatch<Size,3>::test(vec.size(), 3); 00443 00444 TooN::Matrix<3, 3, P> result; 00445 00446 result(0,0) = 0; 00447 result(0,1) = -vec[2]; 00448 result(0,2) = vec[1]; 00449 result(1,0) = vec[2]; 00450 result(1,1) = 0; 00451 result(1,2) = -vec[0]; 00452 result(2,0) = -vec[1]; 00453 result(2,1) = vec[0]; 00454 result(2,2) = 0; 00455 00456 return result; 00457 } 00458 00459 namespace Internal { 00460 template<int Size, typename Precision, typename Base, typename Func, typename Ret> inline Ret accumulate( const Vector<Size, Precision, Base> & v ) { 00461 Func func; 00462 if( v.size() == 0 ) { 00463 return func.null(); // What should we return, exception? 00464 } 00465 func.initialise( v[0], 0 ); 00466 for( int ii = 1; ii < v.size(); ii++ ) { 00467 func( v[ii], ii ); 00468 } 00469 return func.ret(); 00470 } 00471 00472 template<int R, int C, typename Precision, typename Base, typename Func, typename Ret> inline Ret accumulate( const Matrix<R, C, Precision, Base> & m ) { 00473 Func func; 00474 if( m.num_rows() == 0 || m.num_cols() == 0) { 00475 return func.null(); // What should we return, exception? 00476 } 00477 func.initialise( m[0][0], 0, 0 ); 00478 for(int r=0; r<m.num_rows(); r++){ 00479 for(int c=0; c<m.num_cols(); c++){ 00480 func( m[r][c], r, c ); 00481 } 00482 } 00483 return func.ret(); 00484 } 00485 template<int R, int C, typename Precision, typename Base, typename Func, typename Ret> inline Ret accumulate_horizontal( const Matrix<R, C, Precision, Base> & m ) { 00486 Func func( m.num_rows() ); 00487 if( m.num_cols() == 0 || m.num_rows() == 0 ) { 00488 func.null(); // What should we return, exception? 00489 } 00490 for(int r=0; r<m.num_rows(); r++){ 00491 func.initialise( m[r][0], r, 0 ); 00492 for(int c=1; c<m.num_cols(); c++){ 00493 func( m[r][c], r, c ); 00494 } 00495 } 00496 return func.ret(); 00497 } 00498 template<int R, int C, typename Precision, typename Base, typename Func, typename Ret> inline Ret accumulate_vertical( const Matrix<R, C, Precision, Base> & m ) { 00499 Func func( m.num_cols() ); 00500 if( m.num_cols() == 0 || m.num_rows() == 0 ) { 00501 func.null(); // What should we return, exception? 00502 } 00503 for(int c=0; c<m.num_cols(); c++){ 00504 func.initialise( m[0][c], 0, c ); 00505 for(int r=1; r<m.num_rows(); r++){ 00506 func( m[r][c], r, c ); 00507 } 00508 } 00509 return func.ret(); 00510 } 00511 00512 template<typename Precision, typename ComparisonFunctor> 00513 class accumulate_functor_vector { 00514 Precision bestVal; 00515 public: 00516 Precision null() { 00517 return 0; 00518 } 00519 void initialise( Precision initialVal, int ) { 00520 bestVal = initialVal; 00521 } 00522 void operator()( Precision curVal, int ) { 00523 if( ComparisonFunctor()( curVal, bestVal ) ) { 00524 bestVal = curVal; 00525 } 00526 } 00527 Precision ret() { return bestVal; } 00528 }; 00529 template<typename Precision, typename ComparisonFunctor> 00530 class accumulate_element_functor_vector { 00531 Precision bestVal; 00532 int nBestIndex; 00533 public: 00534 std::pair<Precision,int> null() { 00535 return std::pair<Precision,int>( 0, 0 ); 00536 } 00537 void initialise( Precision initialVal, int nIndex ) { 00538 bestVal = initialVal; 00539 nBestIndex = nIndex; 00540 } 00541 void operator()( Precision curVal, int nIndex ) { 00542 if( ComparisonFunctor()( curVal, bestVal ) ) { 00543 bestVal = curVal; 00544 nBestIndex = nIndex; 00545 } 00546 } 00547 std::pair<Precision,int> ret() { 00548 return std::pair<Precision,int>( bestVal, nBestIndex ); 00549 } 00550 }; 00551 template<typename Precision, typename ComparisonFunctor> 00552 class accumulate_functor_matrix { 00553 Precision bestVal; 00554 public: 00555 Precision null() { 00556 return 0; 00557 } 00558 void initialise( Precision initialVal, int, int ) { 00559 bestVal = initialVal; 00560 } 00561 void operator()( Precision curVal, int, int ) { 00562 if( ComparisonFunctor()( curVal, bestVal ) ) { 00563 bestVal = curVal; 00564 } 00565 } 00566 Precision ret() { return bestVal; } 00567 }; 00568 template<typename Precision, typename ComparisonFunctor> 00569 class accumulate_element_functor_matrix { 00570 Precision bestVal; 00571 int nBestRow; 00572 int nBestCol; 00573 public: 00574 std::pair<Precision,std::pair<int,int> > null() { 00575 return std::pair<Precision,std::pair<int,int> >( 0, std::pair<int,int>( 0, 0 ) ); 00576 } 00577 void initialise( Precision initialVal, int nRow, int nCol ) { 00578 bestVal = initialVal; 00579 nBestRow = nRow; 00580 nBestCol = nCol; 00581 } 00582 void operator()( Precision curVal, int nRow, int nCol ) { 00583 if( ComparisonFunctor()( curVal, bestVal ) ) { 00584 bestVal = curVal; 00585 nBestRow = nRow; 00586 nBestCol = nCol; 00587 } 00588 } 00589 std::pair<Precision,std::pair<int,int> > ret() { 00590 return std::pair<Precision,std::pair<int,int> >( bestVal, 00591 std::pair<int,int>( nBestRow, nBestCol ) ); 00592 } 00593 }; 00594 template<typename Precision, typename ComparisonFunctor> 00595 class accumulate_vertical_functor { 00596 Vector<Dynamic,Precision>* bestVal; 00597 public: 00598 accumulate_vertical_functor() { 00599 bestVal = NULL; 00600 } 00601 accumulate_vertical_functor( int nNumCols ) { 00602 bestVal = new Vector<Dynamic,Precision>( nNumCols ); 00603 } 00604 Vector<Dynamic,Precision> null() { 00605 return Vector<Dynamic,Precision>( 0 ); 00606 } 00607 void initialise( Precision initialVal, int, int nCol ) { 00608 (*bestVal)[nCol] = initialVal; 00609 } 00610 void operator()( Precision curVal, int, int nCol ) { 00611 if( ComparisonFunctor()( curVal, (*bestVal)[nCol] ) ) { 00612 (*bestVal)[nCol] = curVal; 00613 } 00614 } 00615 Vector<Dynamic,Precision> ret() { 00616 if( bestVal == NULL ) { 00617 return null(); 00618 } 00619 Vector<Dynamic,Precision> vRet = *bestVal; 00620 delete bestVal; 00621 return vRet; 00622 } 00623 }; 00624 template<typename Precision, typename ComparisonFunctor> 00625 class accumulate_element_vertical_functor { 00626 Vector<Dynamic,Precision>* bestVal; 00627 Vector<Dynamic,Precision>* bestIndices; 00628 public: 00629 accumulate_element_vertical_functor() { 00630 bestVal = NULL; 00631 bestIndices = NULL; 00632 } 00633 accumulate_element_vertical_functor( int nNumCols ) { 00634 bestVal = new Vector<Dynamic,Precision>( nNumCols ); 00635 bestIndices = new Vector<Dynamic,Precision>( nNumCols ); 00636 } 00637 std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> > null() { 00638 Vector<Dynamic,Precision> vEmpty( 0 ); 00639 return std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> >( vEmpty, vEmpty ); 00640 } 00641 void initialise( Precision initialVal, int nRow, int nCol ) { 00642 (*bestVal)[nCol] = initialVal; 00643 (*bestIndices)[nCol] = nRow; 00644 } 00645 void operator()( Precision curVal, int nRow, int nCol ) { 00646 if( ComparisonFunctor()( curVal, (*bestVal)[nCol] ) ) { 00647 (*bestVal)[nCol] = curVal; 00648 (*bestIndices)[nCol] = nRow; 00649 } 00650 } 00651 std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> > ret() { 00652 if( bestVal == NULL ) { 00653 return null(); 00654 } 00655 std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> > vRet = 00656 std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> > (*bestVal, *bestIndices ); 00657 delete bestVal; bestVal = NULL; 00658 delete bestIndices; bestIndices = NULL; 00659 return vRet; 00660 } 00661 }; 00662 template<typename Precision, typename ComparisonFunctor> 00663 class accumulate_horizontal_functor { 00664 Vector<Dynamic,Precision>* bestVal; 00665 public: 00666 accumulate_horizontal_functor() { 00667 bestVal = NULL; 00668 } 00669 accumulate_horizontal_functor( int nNumRows ) { 00670 bestVal = new Vector<Dynamic,Precision>( nNumRows ); 00671 } 00672 Vector<Dynamic,Precision> null() { 00673 return Vector<Dynamic,Precision>( 0 ); 00674 } 00675 void initialise( Precision initialVal, int nRow, int ) { 00676 (*bestVal)[nRow] = initialVal; 00677 } 00678 void operator()( Precision curVal, int nRow, int ) { 00679 if( ComparisonFunctor()( curVal, (*bestVal)[nRow] ) ) { 00680 (*bestVal)[nRow] = curVal; 00681 } 00682 } 00683 Vector<Dynamic,Precision> ret() { 00684 if( bestVal == NULL ) { 00685 return null(); 00686 } 00687 Vector<Dynamic,Precision> vRet = *bestVal; 00688 delete bestVal; bestVal = NULL; 00689 return vRet; 00690 } 00691 }; 00692 template<typename Precision, typename ComparisonFunctor> 00693 class accumulate_element_horizontal_functor { 00694 Vector<Dynamic,Precision>* bestVal; 00695 Vector<Dynamic,Precision>* bestIndices; 00696 public: 00697 accumulate_element_horizontal_functor() { 00698 bestVal = NULL; 00699 bestIndices = NULL; 00700 } 00701 accumulate_element_horizontal_functor( int nNumRows ) { 00702 bestVal = new Vector<Dynamic,Precision>( nNumRows ); 00703 bestIndices = new Vector<Dynamic,Precision>( nNumRows ); 00704 } 00705 std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> > null() { 00706 Vector<Dynamic,Precision> vEmpty( 0 ); 00707 return std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> >( vEmpty, vEmpty ); 00708 } 00709 void initialise( Precision initialVal, int nRow, int nCol ) { 00710 (*bestVal)[nRow] = initialVal; 00711 (*bestIndices)[nRow] = nCol; 00712 } 00713 void operator()( Precision curVal, int nRow, int nCol ) { 00714 if( ComparisonFunctor()( curVal, (*bestVal)[nRow] ) ) { 00715 (*bestVal)[nRow] = curVal; 00716 (*bestIndices)[nRow] = nCol; 00717 } 00718 } 00719 std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> > ret() { 00720 if( bestVal == NULL ) { 00721 return null(); 00722 } 00723 std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> > vRet = 00724 std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> >( *bestVal, *bestIndices ); 00725 delete bestVal; bestVal = NULL; 00726 delete bestIndices; bestIndices = NULL; 00727 return vRet; 00728 } 00729 }; 00730 } 00731 00732 00733 /// Finds the minimal value of a vector. 00734 /// @param v a vector 00735 /// @return the smallest value of v 00736 template<int Size, typename Precision, typename Base> inline Precision min_value( const Vector<Size, Precision, Base>& v) { 00737 typedef Internal::accumulate_functor_vector<Precision, std::less<Precision> > vector_accumulate_functor; 00738 return Internal::accumulate<Size,Precision,Base, 00739 vector_accumulate_functor, Precision >( v ); 00740 } 00741 /// Finds the largest value of a vector. 00742 /// @param v a vector 00743 /// @return the largest value of v 00744 template<int Size, typename Precision, typename Base> inline Precision max_value( const Vector<Size, Precision, Base>& v) { 00745 typedef Internal::accumulate_functor_vector<Precision, std::greater<Precision> > vector_accumulate_functor; 00746 return Internal::accumulate<Size,Precision,Base, 00747 vector_accumulate_functor, Precision >( v ); 00748 } 00749 00750 /// Finds the smallest value of a matrix. 00751 /// @param m a matrix 00752 /// @return the smallest value of m 00753 template<int R, int C, typename Precision, typename Base> inline Precision min_value( const Matrix<R, C, Precision, Base> & m) { 00754 typedef Internal::accumulate_functor_matrix<Precision, std::less<Precision> > matrix_accumulate_functor; 00755 return Internal::accumulate<R,C,Precision,Base, 00756 matrix_accumulate_functor, Precision>( m ); 00757 } 00758 /// Finds the largest value of a matrix. 00759 /// @param m a matrix 00760 /// @return the largest value of m 00761 template<int R, int C, typename Precision, typename Base> inline Precision max_value( const Matrix<R, C, Precision, Base> & m) { 00762 typedef Internal::accumulate_functor_matrix<Precision, std::greater<Precision> > matrix_accumulate_functor; 00763 return Internal::accumulate<R,C,Precision,Base, 00764 matrix_accumulate_functor, Precision>( m ); 00765 } 00766 /// Finds the smallest values of each column of a matrix. 00767 /// @param m a matrix 00768 /// @return a vector of size C 00769 template<int R, int C, typename Precision, typename Base> inline Vector<Dynamic,Precision> min_value_vertical( const Matrix<R, C, Precision, Base> & m) { 00770 typedef Internal::accumulate_vertical_functor<Precision,std::less<Precision> > matrix_accumulate_vertical_functor; 00771 return Internal::accumulate_vertical<R,C,Precision,Base, 00772 matrix_accumulate_vertical_functor, Vector<Dynamic,Precision> >( m ); 00773 } 00774 /// Finds the largest values of each column of a matrix. 00775 /// @param m a matrix 00776 /// @return a vector of size C 00777 template<int R, int C, typename Precision, typename Base> inline Vector<Dynamic,Precision> max_value_vertical( const Matrix<R, C, Precision, Base> & m) { 00778 typedef Internal::accumulate_vertical_functor<Precision,std::greater<Precision> > matrix_accumulate_vertical_functor; 00779 return Internal::accumulate_vertical<R,C,Precision,Base, 00780 matrix_accumulate_vertical_functor, Vector<Dynamic,Precision> >( m ); 00781 } 00782 /// Finds the smallest values of each row of a matrix. 00783 /// @param m a matrix 00784 /// @return a vector of size R 00785 template<int R, int C, typename Precision, typename Base> inline Vector<Dynamic,Precision> min_value_horizontal( const Matrix<R, C, Precision, Base> & m) { 00786 typedef Internal::accumulate_horizontal_functor<Precision,std::less<Precision> > matrix_accumulate_horizontal_functor; 00787 return Internal::accumulate_horizontal<R,C,Precision,Base, 00788 matrix_accumulate_horizontal_functor, Vector<Dynamic,Precision> >( m ); 00789 } 00790 /// Finds the largest values of each row of a matrix. 00791 /// @param m a matrix 00792 /// @return a vector of size R 00793 template<int R, int C, typename Precision, typename Base> inline Vector<Dynamic,Precision> max_value_horizontal( const Matrix<R, C, Precision, Base> & m) { 00794 typedef Internal::accumulate_horizontal_functor<Precision,std::greater<Precision> > matrix_accumulate_horizontal_functor; 00795 return Internal::accumulate_horizontal<R,C,Precision,Base, 00796 matrix_accumulate_horizontal_functor, Vector<Dynamic,Precision> >( m ); 00797 } 00798 /// Finds the smallest value of a vector and its index. 00799 /// @param v a vector 00800 /// @return a pair containing the smallest value and its index 00801 template<int Size, typename Precision, typename Base> inline std::pair<Precision,int> min_element( const Vector<Size, Precision, Base>& v) { 00802 typedef Internal::accumulate_element_functor_vector<Precision, std::less<Precision> > vector_accumulate_functor; 00803 return Internal::accumulate<Size,Precision,Base, 00804 vector_accumulate_functor, std::pair<Precision,int> >( v ); 00805 } 00806 /// Finds the largest value of a vector and its index. 00807 /// @param v a vector 00808 /// @return a pair containing the largest value and its index 00809 template<int Size, typename Precision, typename Base> inline std::pair<Precision,int> max_element( const Vector<Size, Precision, Base>& v) { 00810 typedef Internal::accumulate_element_functor_vector<Precision, std::greater<Precision> > vector_accumulate_functor; 00811 return Internal::accumulate<Size,Precision,Base, 00812 vector_accumulate_functor, std::pair<Precision,int> >( v ); 00813 } 00814 /// Finds the smallest value of a matrix and its row and column. 00815 /// @param m a matrix 00816 /// @return a pair containing the smallest value and a pair 00817 /// containing its row and column 00818 template<int R, int C, typename Precision, typename Base> inline std::pair<Precision,std::pair<int,int> > min_element( const Matrix<R, C, Precision, Base> & m) { 00819 typedef Internal::accumulate_element_functor_matrix<Precision, std::less<Precision> > matrix_accumulate_functor; 00820 typedef std::pair<Precision,std::pair<int,int> > Ret; 00821 return Internal::accumulate<R,C,Precision,Base, 00822 matrix_accumulate_functor, Ret>( m ); 00823 } 00824 /// Finds the largest value of a matrix and its row and column. 00825 /// @param m a matrix 00826 /// @return a pair containing the largest value and a pair 00827 /// containing its row and column 00828 template<int R, int C, typename Precision, typename Base> inline std::pair<Precision,std::pair<int,int> > max_element( const Matrix<R, C, Precision, Base> & m) { 00829 typedef Internal::accumulate_element_functor_matrix<Precision, std::greater<Precision> > matrix_accumulate_functor; 00830 typedef std::pair<Precision,std::pair<int,int> > Ret; 00831 return Internal::accumulate<R,C,Precision,Base, 00832 matrix_accumulate_functor, Ret>( m ); 00833 } 00834 /// Finds the smallest values of each column of a matrix and their 00835 /// indices. 00836 /// @param m a matrix 00837 /// @return a pair of vectors of size C containg the values and 00838 /// their indices 00839 template<int R, int C, typename Precision, typename Base> inline std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> > min_element_vertical( const Matrix<R, C, Precision, Base> & m) { 00840 typedef Internal::accumulate_element_vertical_functor<Precision,std::less<Precision> > matrix_accumulate_vertical_functor; 00841 typedef std::pair<Vector< Dynamic, Precision >,Vector< Dynamic, Precision > > Ret; 00842 return Internal::accumulate_vertical<R,C,Precision,Base, 00843 matrix_accumulate_vertical_functor, Ret >( m ); 00844 } 00845 /// Finds the largest values of each column of a matrix and their 00846 /// indices. 00847 /// @param m a matrix 00848 /// @return a pair of vectors of size C containg the values and 00849 /// their indices 00850 template<int R, int C, typename Precision, typename Base> inline std::pair<Vector< Dynamic, Precision >,Vector< Dynamic, Precision > > max_element_vertical( const Matrix<R, C, Precision, Base> & m) { 00851 typedef Internal::accumulate_element_vertical_functor<Precision,std::greater<Precision> > matrix_accumulate_vertical_functor; 00852 typedef std::pair<Vector< Dynamic, Precision >,Vector< Dynamic, Precision > > Ret; 00853 return Internal::accumulate_vertical<R,C,Precision,Base, 00854 matrix_accumulate_vertical_functor, Ret >( m ); 00855 } 00856 /// Finds the smallest values of each row of a matrix and their 00857 /// indices. 00858 /// @param m a matrix 00859 /// @return a pair of vectors of size R containg the values and 00860 /// their indices 00861 template<int R, int C, typename Precision, typename Base> inline std::pair<Vector< Dynamic, Precision >,Vector< Dynamic, Precision > > min_element_horizontal( const Matrix<R, C, Precision, Base> & m) { 00862 typedef Internal::accumulate_element_horizontal_functor<Precision,std::less<Precision> > matrix_accumulate_vertical_functor; 00863 typedef std::pair<Vector< Dynamic, Precision >,Vector< Dynamic, Precision > > Ret; 00864 return Internal::accumulate_horizontal<R,C,Precision,Base, 00865 matrix_accumulate_vertical_functor, Ret >( m ); 00866 } 00867 /// Finds the largest values of each row of a matrix and their 00868 /// indices. 00869 /// @param m a matrix 00870 /// @return a pair of vectors of size R containg the values and 00871 /// their indices 00872 template<int R, int C, typename Precision, typename Base> inline std::pair<Vector< Dynamic, Precision >,Vector< Dynamic, Precision > > max_element_horizontal( const Matrix<R, C, Precision, Base> & m) { 00873 typedef Internal::accumulate_element_horizontal_functor<Precision,std::greater<Precision> > matrix_accumulate_vertical_functor; 00874 typedef std::pair<Vector< Dynamic, Precision >,Vector< Dynamic, Precision > > Ret; 00875 return Internal::accumulate_horizontal<R,C,Precision,Base, 00876 matrix_accumulate_vertical_functor, Ret >( m ); 00877 } 00878 } 00879 #endif