core/vnl/vnl_diag_matrix.h
Go to the documentation of this file.
00001 // This is core/vnl/vnl_diag_matrix.h
00002 #ifndef vnl_diag_matrix_h_
00003 #define vnl_diag_matrix_h_
00004 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00005 #pragma interface
00006 #endif
00007 //:
00008 // \file
00009 // \brief Contains class for diagonal matrices
00010 // \author Andrew W. Fitzgibbon (Oxford RRG)
00011 // \date   5 Aug 1996
00012 //
00013 // \verbatim
00014 //  Modifications
00015 //   IMS (Manchester) 16 Mar 2001: Tidied up the documentation + added binary_io
00016 //   Feb.2002 - Peter Vanroose - brief doxygen comment placed on single line
00017 //   Sep.2002 - Peter Vanroose - Added operator+, operator-, operator*
00018 //   Mar.2004 - Peter Vanroose - removed deprecated resize()
00019 //   Oct.2010 - Peter Vanroose - mutators and setters now return *this
00020 //   Jan.2011 - Peter Vanroose - added methods set_diagonal() & get_diagonal()
00021 // \endverbatim
00022 
00023 #include <vcl_cassert.h>
00024 #include <vcl_iosfwd.h>
00025 #include <vnl/vnl_vector.h>
00026 #include <vnl/vnl_matrix.h>
00027 
00028 // forward declarations
00029 template <class T> class vnl_diag_matrix;
00030 template <class T> vnl_vector<T> operator*(vnl_diag_matrix<T> const&, vnl_vector<T> const&);
00031 
00032 //: stores a diagonal matrix as a single vector.
00033 //  vnl_diag_matrix stores a diagonal matrix for time and space efficiency.
00034 //  Specifically, only the diagonal elements are stored, and some matrix
00035 //  operations (currently *, + and -) are overloaded to use more efficient
00036 //  algorithms.
00037 
00038 export
00039 template <class T>
00040 class vnl_diag_matrix
00041 {
00042   vnl_vector<T> diagonal_;
00043 
00044  public:
00045   vnl_diag_matrix() : diagonal_() {}
00046 
00047   //: Construct an empty diagonal matrix.
00048   vnl_diag_matrix(unsigned nn) : diagonal_(nn) {}
00049 
00050   //: Construct a diagonal matrix with diagonal elements equal to value.
00051   vnl_diag_matrix(unsigned nn, T const& value) : diagonal_(nn, value) {}
00052 
00053   //: Construct a diagonal matrix from a vnl_vector.
00054   //  The vector elements become the diagonal elements.
00055   vnl_diag_matrix(vnl_vector<T> const& that): diagonal_(that) {}
00056  ~vnl_diag_matrix() {}
00057 
00058   inline vnl_diag_matrix& operator=(vnl_diag_matrix<T> const& that) {
00059     this->diagonal_ = that.diagonal_;
00060     return *this;
00061   }
00062 
00063   // Operations----------------------------------------------------------------
00064 
00065   //: In-place arithmetic operation
00066   inline vnl_diag_matrix<T>& operator*=(T v) { diagonal_ *= v; return *this; }
00067   //: In-place arithmetic operation
00068   inline vnl_diag_matrix<T>& operator/=(T v) { diagonal_ /= v; return *this; }
00069 
00070   // Computations--------------------------------------------------------------
00071 
00072   vnl_diag_matrix& invert_in_place();
00073   T determinant() const;
00074   vnl_vector<T> solve(vnl_vector<T> const& b) const;
00075   void solve(vnl_vector<T> const& b, vnl_vector<T>* out) const;
00076 
00077   // Data Access---------------------------------------------------------------
00078 
00079   inline T operator () (unsigned i, unsigned j) const {
00080     return (i != j) ? T(0) : diagonal_[i];
00081   }
00082 
00083   inline T& operator () (unsigned i, unsigned j) {
00084     assert(i == j);
00085     return diagonal_[i];
00086   }
00087   inline T& operator() (unsigned i) { return diagonal_[i]; }
00088   inline T const& operator() (unsigned i) const { return diagonal_[i]; }
00089 
00090   inline T& operator[] (unsigned i) { return diagonal_[i]; }
00091   inline T const& operator[] (unsigned i) const { return diagonal_[i]; }
00092 
00093   //: set element with boundary checks.
00094   inline void put (unsigned r, unsigned c, T const& v) {
00095     assert(r == c); assert (r<size()); diagonal_[r] = v;
00096   }
00097 
00098   //: get element with boundary checks.
00099   inline T get (unsigned r, unsigned c) const {
00100     assert(r == c); assert (r<size()); return diagonal_[r];
00101   }
00102 
00103   //: Return a vector (copy) with the content of the (main) diagonal
00104   inline vnl_vector<T> get_diagonal() const { return diagonal_; }
00105 
00106   //: Return diagonal elements as a vector
00107   inline vnl_vector<T> const& diagonal() const { return diagonal_; }
00108 
00109   //: Set all diagonal elements of matrix to specified value.
00110   inline vnl_diag_matrix& fill_diagonal (T const& v) { diagonal_.fill(v); return *this; }
00111 
00112   //: Sets the diagonal elements of this matrix to the specified list of values.
00113   inline vnl_diag_matrix& set_diagonal(vnl_vector<T> const& v) { diagonal_ = v; return *this; }
00114 
00115   // iterators
00116 
00117   typedef typename vnl_vector<T>::iterator iterator;
00118   inline iterator begin() { return diagonal_.begin(); }
00119   inline iterator end() { return diagonal_.end(); }
00120   typedef typename vnl_vector<T>::const_iterator const_iterator;
00121   inline const_iterator begin() const { return diagonal_.begin(); }
00122   inline const_iterator end() const { return diagonal_.end(); }
00123 
00124   inline unsigned size() const { return diagonal_.size(); }
00125   inline unsigned rows() const { return diagonal_.size(); }
00126   inline unsigned cols() const { return diagonal_.size(); }
00127   inline unsigned columns() const { return diagonal_.size(); }
00128 
00129   // Need this until we add a vnl_diag_matrix ctor to vnl_matrix;
00130   inline vnl_matrix<T> asMatrix() const;
00131 
00132   inline vnl_matrix<T> as_ref() const { return asMatrix(); }
00133 
00134   // This is as good as a vnl_diag_matrix ctor for vnl_matrix:
00135   inline operator vnl_matrix<T> () const { return asMatrix(); }
00136 
00137   inline void set_size(int n) { diagonal_.set_size(n); }
00138 
00139   inline void clear() { diagonal_.clear(); }
00140   inline vnl_diag_matrix& fill(T const &x) { diagonal_.fill(x); return *this; }
00141 
00142   //: Return pointer to the diagonal elements as a contiguous 1D C array;
00143   inline T*       data_block()       { return diagonal_.data_block(); }
00144   inline T const* data_block() const { return diagonal_.data_block(); }
00145 
00146   //: Set diagonal elements using vector
00147   inline vnl_diag_matrix& set(vnl_vector<T> const& v)  { diagonal_=v; return *this; }
00148 
00149  private:
00150   #if VCL_NEED_FRIEND_FOR_TEMPLATE_OVERLOAD
00151   friend vnl_vector<T> operator*(vnl_diag_matrix<T> const&,vnl_vector<T> const&);
00152   #endif
00153 };
00154 
00155 //:
00156 // \relatesalso vnl_diag_matrix
00157 template <class T>
00158 vcl_ostream& operator<< (vcl_ostream&, vnl_diag_matrix<T> const&);
00159 
00160 //: Convert a vnl_diag_matrix to a Matrix.
00161 template <class T>
00162 inline vnl_matrix<T> vnl_diag_matrix<T>::asMatrix() const
00163 {
00164   unsigned len = diagonal_.size();
00165   vnl_matrix<T> ret(len, len);
00166   for (unsigned i = 0; i < len; ++i)
00167   {
00168     unsigned j;
00169     for (j = 0; j < i; ++j)
00170       ret(i,j) = T(0);
00171     for (j = i+1; j < len; ++j)
00172       ret(i,j) = T(0);
00173     ret(i,i) = diagonal_[i];
00174   }
00175   return ret;
00176 }
00177 
00178 //: Invert a vnl_diag_matrix in-situ.
00179 // Just replaces each element with its reciprocal.
00180 template <class T>
00181 inline vnl_diag_matrix<T>& vnl_diag_matrix<T>::invert_in_place()
00182 {
00183   unsigned len = diagonal_.size();
00184   T* d = data_block();
00185   T one = T(1);
00186   for (unsigned i = 0; i < len; ++i)
00187     d[i] = one / d[i];
00188   return *this;
00189 }
00190 
00191 //: Return determinant as product of diagonal values.
00192 template <class T>
00193 inline T vnl_diag_matrix<T>::determinant() const
00194 {
00195   T det = T(1);
00196   T const* d = data_block();
00197   unsigned len = diagonal_.size();
00198   for (unsigned i = 0; i < len; ++i)
00199     det *= d[i];
00200   return det;
00201 }
00202 
00203 //: Multiply two vnl_diag_matrices.  Just multiply the diag elements - n flops
00204 // \relatesalso vnl_diag_matrix
00205 template <class T>
00206 inline vnl_diag_matrix<T> operator* (vnl_diag_matrix<T> const& A, vnl_diag_matrix<T> const& B)
00207 {
00208   assert(A.size() == B.size());
00209   vnl_diag_matrix<T> ret = A;
00210   for (unsigned i = 0; i < A.size(); ++i)
00211     ret(i,i) *= B(i,i);
00212   return ret;
00213 }
00214 
00215 //: Multiply a vnl_matrix by a vnl_diag_matrix.  Just scales the columns - mn flops
00216 // \relatesalso vnl_diag_matrix
00217 // \relatesalso vnl_matrix
00218 template <class T>
00219 inline vnl_matrix<T> operator* (vnl_matrix<T> const& A, vnl_diag_matrix<T> const& D)
00220 {
00221   assert(A.columns() == D.size());
00222   vnl_matrix<T> ret(A.rows(), A.columns());
00223   for (unsigned i = 0; i < A.rows(); ++i)
00224     for (unsigned j = 0; j < A.columns(); ++j)
00225       ret(i,j) = A(i,j) * D(j,j);
00226   return ret;
00227 }
00228 
00229 //: Multiply a vnl_diag_matrix by a vnl_matrix.  Just scales the rows - mn flops
00230 // \relatesalso vnl_diag_matrix
00231 // \relatesalso vnl_matrix
00232 template <class T>
00233 inline vnl_matrix<T> operator* (vnl_diag_matrix<T> const& D, vnl_matrix<T> const& A)
00234 {
00235   assert(A.rows() == D.size());
00236   vnl_matrix<T> ret(A.rows(), A.columns());
00237   T const* d = D.data_block();
00238   for (unsigned i = 0; i < A.rows(); ++i)
00239     for (unsigned j = 0; j < A.columns(); ++j)
00240       ret(i,j) = A(i,j) * d[i];
00241   return ret;
00242 }
00243 
00244 //: Add two vnl_diag_matrices.  Just add the diag elements - n flops
00245 // \relatesalso vnl_diag_matrix
00246 template <class T>
00247 inline vnl_diag_matrix<T> operator+ (vnl_diag_matrix<T> const& A, vnl_diag_matrix<T> const& B)
00248 {
00249   assert(A.size() == B.size());
00250   vnl_diag_matrix<T> ret = A;
00251   for (unsigned i = 0; i < A.size(); ++i)
00252     ret(i,i) += B(i,i);
00253   return ret;
00254 }
00255 
00256 //: Add a vnl_diag_matrix to a vnl_matrix.  n adds, mn copies.
00257 // \relatesalso vnl_diag_matrix
00258 // \relatesalso vnl_matrix
00259 template <class T>
00260 inline vnl_matrix<T> operator+ (vnl_matrix<T> const& A, vnl_diag_matrix<T> const& D)
00261 {
00262   const unsigned n = D.size();
00263   assert(A.rows() == n); assert(A.columns() == n);
00264   vnl_matrix<T> ret(A);
00265   T const* d = D.data_block();
00266   for (unsigned j = 0; j < n; ++j)
00267     ret(j,j) += d[j];
00268   return ret;
00269 }
00270 
00271 //: Add a vnl_matrix to a vnl_diag_matrix.  n adds, mn copies.
00272 // \relatesalso vnl_diag_matrix
00273 // \relatesalso vnl_matrix
00274 template <class T>
00275 inline vnl_matrix<T> operator+ (vnl_diag_matrix<T> const& D, vnl_matrix<T> const& A)
00276 {
00277   return A + D;
00278 }
00279 
00280 //: Subtract two vnl_diag_matrices.  Just subtract the diag elements - n flops
00281 // \relatesalso vnl_diag_matrix
00282 template <class T>
00283 inline vnl_diag_matrix<T> operator- (vnl_diag_matrix<T> const& A, vnl_diag_matrix<T> const& B)
00284 {
00285   assert(A.size() == B.size());
00286   vnl_diag_matrix<T> ret = A;
00287   for (unsigned i = 0; i < A.size(); ++i)
00288     ret(i,i) -= B(i,i);
00289   return ret;
00290 }
00291 
00292 //: Subtract a vnl_diag_matrix from a vnl_matrix.  n adds, mn copies.
00293 // \relatesalso vnl_diag_matrix
00294 // \relatesalso vnl_matrix
00295 template <class T>
00296 inline vnl_matrix<T> operator- (vnl_matrix<T> const& A, vnl_diag_matrix<T> const& D)
00297 {
00298   const unsigned n = D.size();
00299   assert(A.rows() == n); assert(A.columns() == n);
00300   vnl_matrix<T> ret(A);
00301   T const* d = D.data_block();
00302   for (unsigned j = 0; j < n; ++j)
00303     ret(j,j) -= d[j];
00304   return ret;
00305 }
00306 
00307 //: Subtract a vnl_matrix from a vnl_diag_matrix.  n adds, mn copies.
00308 // \relatesalso vnl_diag_matrix
00309 // \relatesalso vnl_matrix
00310 template <class T>
00311 inline vnl_matrix<T> operator- (vnl_diag_matrix<T> const& D, vnl_matrix<T> const& A)
00312 {
00313   const unsigned n = D.size();
00314   assert(A.rows() == n); assert(A.columns() == n);
00315   vnl_matrix<T> ret(n, n);
00316   T const* d = D.data_block();
00317   for (unsigned i = 0; i < n; ++i)
00318   {
00319     for (unsigned j = 0; j < i; ++j)
00320       ret(i,j) = -A(i,j);
00321     for (unsigned j = i+1; j < n; ++j)
00322       ret(i,j) = -A(i,j);
00323     ret(i,i) = d[i] - A(i,i);
00324   }
00325   return ret;
00326 }
00327 
00328 //: Multiply a vnl_diag_matrix by a vnl_vector.  n flops.
00329 // \relatesalso vnl_diag_matrix
00330 // \relatesalso vnl_vector
00331 template <class T>
00332 inline vnl_vector<T> operator* (vnl_diag_matrix<T> const& D, vnl_vector<T> const& A)
00333 {
00334   assert(A.size() == D.size());
00335   return element_product(D.diagonal(), A);
00336 }
00337 
00338 //: Multiply a vnl_vector by a vnl_diag_matrix.  n flops.
00339 // \relatesalso vnl_diag_matrix
00340 // \relatesalso vnl_vector
00341 template <class T>
00342 inline vnl_vector<T> operator* (vnl_vector<T> const& A, vnl_diag_matrix<T> const& D)
00343 {
00344   assert(A.size() == D.size());
00345   return element_product(D.diagonal(), A);
00346 }
00347 
00348 #endif // vnl_diag_matrix_h_