Main Page   Groups   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Concepts

vnl_diag_matrix.h

Go to the documentation of this file.
00001 #ifndef vnl_diag_matrix_h_
00002 #define vnl_diag_matrix_h_
00003 
00004 //:
00005 //  \file
00006 //  \brief Contains class for diagonal matrices
00007 //  \author Andrew W. Fitzgibbon (Oxford RRG) 5/8/96
00008 //
00009 //
00010 // \verbatim
00011 //  Modifications
00012 //  IMS (Manchester) 16/03/2001: Tidied up the documentation + added binary_io
00013 // \endverbatim
00014 
00015 
00016 
00017 //  forward declare friend functions
00018 //  template <class T> class vnl_diag_matrix;
00019 //  template<class T> bool epsilon_equals
00020 
00021 //        (vnl_diag_matrix<T> const& m1, vnl_diag_matrix<T> const& m2, double alt_epsilon = 0);
00022 
00023 
00024 #include <vcl_cassert.h>
00025 #include <vnl/vnl_vector.h>
00026 #include <vnl/vnl_matrix.h>
00027 
00028 //: stores a diagonal matrix as a single vector
00029 //  vnl_diag_matrix stores a diagonal matrix for time and space efficiency.
00030 //  Specifically, only the diagonal elements are stored, and some matrix
00031 //  operations (currently *, + and -) are overloaded to use more efficient
00032 //  algorithms.
00033 
00034 export
00035 template <class T>
00036 class vnl_diag_matrix {
00037 public:
00038   vnl_diag_matrix() {}
00039 
00040   //: Construct an empty diagonal matrix.
00041   vnl_diag_matrix(unsigned nn) : diagonal_(nn) {}
00042 
00043   //: Construct a diagonal matrix with diagonal elements equal to value.
00044   vnl_diag_matrix(unsigned nn, T const& value) : diagonal_(nn, value) {}
00045 
00046   //: Construct a diagonal matrix from a Vector.  The vector elements become
00047   // the diagonal elements.
00048   vnl_diag_matrix(vnl_vector<T> const& that): diagonal_(that) {}
00049  ~vnl_diag_matrix() {}
00050 
00051 
00052   vnl_diag_matrix& operator=(vnl_diag_matrix<T> const& that) {
00053   this->diagonal_ = that.diagonal_;
00054   return *this;
00055     }
00056 
00057   // Operations----------------------------------------------------------------
00058   //: In-place arithmetic operations
00059   vnl_diag_matrix<T>& operator*=(T v) { diagonal_ *= v; return *this; }
00060   //: In-place arithmetic operations
00061   vnl_diag_matrix<T>& operator/=(T v) { diagonal_ /= v; return *this; }
00062 
00063   // Computations--------------------------------------------------------------
00064   void invert_in_place();
00065   T determinant() const;
00066   vnl_vector<T> solve(vnl_vector<T> const& b);
00067   void solve(vnl_vector<T> const& b, vnl_vector<T>* out);
00068 
00069   // Data Access---------------------------------------------------------------
00070   T operator () (unsigned i, unsigned j) const {
00071     return (i != j) ? T(0) : diagonal_[i];
00072   }
00073 
00074   T& operator () (unsigned i, unsigned j) {
00075     assert(i == j);
00076     return diagonal_[i];
00077   }
00078   T& operator() (unsigned i) { return diagonal_[i]; }
00079   T const& operator() (unsigned i) const { return diagonal_[i]; }
00080 
00081   // iterators
00082   typedef typename vnl_vector<T>::iterator iterator;
00083   inline iterator begin() { return diagonal_.begin(); }
00084   inline iterator end() { return diagonal_.end(); }
00085   typedef typename vnl_vector<T>::const_iterator const_iterator;
00086   inline const_iterator begin() const { return diagonal_.begin(); }
00087   inline const_iterator end() const { return diagonal_.end(); }
00088 
00089   unsigned size() const { return diagonal_.size(); }
00090   unsigned n() const { return diagonal_.size(); } // ** deprecated ? **
00091   unsigned rows() const { return diagonal_.size(); }
00092   unsigned cols() const { return diagonal_.size(); }
00093   unsigned columns() const { return diagonal_.size(); }
00094 
00095   // Need this until we add a vnl_diag_matrix ctor to vnl_matrix;
00096   inline vnl_matrix<T> asMatrix() const;
00097 
00098   void resize(int n) { diagonal_.resize(n); }
00099   void clear() { diagonal_.clear(); }
00100 
00101   //: Return pointer to the diagonal elements as a contiguous 1D C array;
00102   T*       data_block()       { return diagonal_.data_block(); }
00103   T const* data_block() const { return diagonal_.data_block(); }
00104 
00105   //: Return diagonal elements as a vector
00106   vnl_vector<T> const& diagonal() const { return diagonal_; }
00107 
00108   //: Set diagonal elements using vector
00109   void set(vnl_vector<T> const& v)  { diagonal_=v; }
00110 
00111 protected:
00112   vnl_vector<T> diagonal_;
00113 
00114 private:
00115   //  // This is private because it's not really a matrix operation.
00116   //  T operator()(unsigned i) const {
00117   //    return diagonal_[i];
00118   //  }
00119 
00120   #if VCL_NEED_FRIEND_FOR_TEMPLATE_OVERLOAD
00121   friend vnl_vector<T> operator*(vnl_diag_matrix<T> const&,vnl_vector<T> const&);
00122   #endif
00123 };
00124 
00125 template <class T>
00126 vcl_ostream& operator<< (vcl_ostream&, vnl_diag_matrix<T> const&);
00127 
00128 
00129 
00130 
00131 // Define this now,
00132 // #define IUE_DEFINED_vnl_diag_matrix
00133 
00134 template <class T> vcl_ostream& operator<< (vcl_ostream&, vnl_diag_matrix<T> const&);
00135 
00136 //: Convert a vnl_diag_matrix to a Matrix.
00137 template <class T>
00138 inline vnl_matrix<T> vnl_diag_matrix<T>::asMatrix() const
00139 {
00140   unsigned len = diagonal_.size();
00141   vnl_matrix<T> ret(len, len);
00142   for(unsigned i = 0; i < len; ++i)
00143     for(unsigned j = 0; j < len; ++j)
00144       ret(i,j) = ((i == j) ? diagonal_[i] : T(0));
00145   return ret;
00146 }
00147 
00148 //: Invert a vnl_diag_matrix in-situ.
00149 // Just replaces each element with its reciprocal.
00150 template <class T>
00151 inline void vnl_diag_matrix<T>::invert_in_place()
00152 {
00153   unsigned len = diagonal_.size();
00154   T* d = data_block();
00155   T one = T(1);
00156   for(unsigned i = 0; i < len; ++i)
00157     d[i] = one / d[i];
00158 }
00159 
00160 //: Return determinant as product of diagonal values.
00161 template <class T>
00162 inline T vnl_diag_matrix<T>::determinant() const
00163 {
00164   T det = T(1);
00165   T const* d = data_block();
00166   unsigned len = diagonal_.size();
00167   for(unsigned i = 0; i < len; ++i)
00168     det *= d[i];
00169   return det;
00170 }
00171 
00172 //: Multiply a Matrix by a vnl_diag_matrix.  Just scales the columns - mn flops
00173 template <class T>
00174 inline vnl_matrix<T> operator* (vnl_matrix<T> const& A, vnl_diag_matrix<T> const& D)
00175 {
00176   vnl_matrix<T> ret(A.rows(), A.columns());
00177   for(unsigned i = 0; i < A.rows(); ++i)
00178     for(unsigned j = 0; j < A.columns(); ++j)
00179       ret(i,j) = A(i,j) * D(j,j);
00180   return ret;
00181 }
00182 
00183 //: Multiply a vnl_diag_matrix by a Matrix.  Just scales the rows - mn flops
00184 template <class T>
00185 inline vnl_matrix<T> operator* (vnl_diag_matrix<T> const& D, vnl_matrix<T> const& A)
00186 {
00187   vnl_matrix<T> ret(A.rows(), A.columns());
00188   T const* d = D.data_block();
00189   for(unsigned i = 0; i < A.rows(); ++i)
00190     for(unsigned j = 0; j < A.columns(); ++j)
00191       ret(i,j) = A(i,j) * d[i];
00192   return ret;
00193 }
00194 
00195 //: Add a vnl_diag_matrix to a Matrix.  n adds, mn copies.
00196 template <class T>
00197 inline vnl_matrix<T> operator + (vnl_matrix<T> const& A, vnl_diag_matrix<T> const& D)
00198 {
00199   unsigned n = D.n();
00200   vnl_matrix<T> ret(A);
00201   T const* d = D.data_block();
00202   for(unsigned j = 0; j < n; ++j)
00203     ret(j,j) += d[j];
00204   return ret;
00205 }
00206 
00207 //: Add a Matrix to a vnl_diag_matrix.  n adds, mn copies.
00208 template <class T>
00209 inline vnl_matrix<T> operator + (vnl_diag_matrix<T> const& D, vnl_matrix<T> const& A)
00210 {
00211   return A + D;
00212 }
00213 
00214 //: Subtract a vnl_diag_matrix from a Matrix.  n adds, mn copies.
00215 template <class T>
00216 inline vnl_matrix<T> operator - (vnl_matrix<T> const& A, vnl_diag_matrix<T> const& D)
00217 {
00218   unsigned n = D.n();
00219   vnl_matrix<T> ret(A);
00220   T const* d = D.data_block();
00221   for(unsigned j = 0; j < n; ++j)
00222     ret(j,j) -= d[j];
00223   return ret;
00224 }
00225 
00226 //: Subtract a Matrix from a vnl_diag_matrix.  n adds, mn copies.
00227 template <class T>
00228 inline vnl_matrix<T> operator - (vnl_diag_matrix<T> const& D, vnl_matrix<T> const& A)
00229 {
00230   unsigned n = D.n();
00231   vnl_matrix<T> ret(n, n);
00232   T const* d = D.data_block();
00233   for(unsigned i = 0; i < n; ++i)
00234     for(unsigned j = 0; j < n; ++j)
00235       if (i == j)
00236     ret(i,j) = d[j] - A(i,j);
00237       else
00238     ret(i,j) = -A(i,j);
00239   return ret;
00240 }
00241 
00242 //: Multiply a vnl_diag_matrix by a Vector.  n flops.
00243 template <class T>
00244 inline vnl_vector<T> operator* (vnl_diag_matrix<T> const& D, vnl_vector<T> const& A)
00245 {
00246   return element_product(D.diagonal(), A);
00247 }
00248 
00249 //: Multiply a Vector by a vnl_diag_matrix.  n flops.
00250 template <class T>
00251 inline vnl_vector<T> operator* (vnl_vector<T> const& A, vnl_diag_matrix<T> const& D)
00252 {
00253   return element_product(D.diagonal(), A);
00254 }
00255 
00256 
00257 #endif // vnl_diag_matrix_h_

Generated at Fri May 21 01:15:47 2004 for ITK by doxygen 1.2.15 written by Dimitri van Heesch, © 1997-2000