00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #ifndef __itkMatrix_h
00018 #define __itkMatrix_h
00019
00020
00021 #include "itkPoint.h"
00022 #include "itkVector.h"
00023 #include "itkCovariantVector.h"
00024 #include "vnl/vnl_matrix_fixed.h"
00025 #include "vnl/algo/vnl_matrix_inverse.h"
00026 #include "vnl/vnl_transpose.h"
00027 #include "vnl/vnl_matrix.h"
00028 #include <vnl/algo/vnl_determinant.h>
00029
00030 namespace itk
00031 {
00032
00041 template<class T, unsigned int NRows=3, unsigned int NColumns=3>
00042 class Matrix
00043 {
00044 public:
00046 typedef Matrix Self;
00047
00049 typedef T ValueType;
00050 typedef T ComponentType;
00051
00053 itkStaticConstMacro(RowDimensions, unsigned int, NRows);
00054 itkStaticConstMacro(ColumnDimensions, unsigned int, NColumns);
00056
00058 typedef vnl_matrix_fixed<T,NRows,NColumns> InternalMatrixType;
00059
00061 Vector<T,NRows> operator*(const Vector<T,NColumns> & vector) const;
00062
00064 Point<T,NRows> operator*(const Point<T,NColumns> & vector) const;
00065
00067 CovariantVector<T,NRows>
00068 operator*(const CovariantVector<T,NColumns> & vector) const;
00069
00071 Self operator*(const Self & matrix) const;
00072
00074 Self operator+(const Self & matrix) const;
00075 const Self & operator+=(const Self & matrix );
00077
00079 Self operator-(const Self & matrix) const;
00080 const Self & operator-=(const Self & matrix );
00082
00084 vnl_matrix<T> operator*(const vnl_matrix<T> & matrix) const;
00085
00087 void operator*=(const Self & matrix);
00088
00090 void operator*=(const vnl_matrix<T> & matrix);
00091
00093 vnl_vector<T> operator*(const vnl_vector<T> & matrix) const;
00094
00096 void operator*=(const T & value)
00097 { m_Matrix *= value; }
00098
00100 Self operator*(const T & value)
00101 {
00102 Self result( *this );
00103 result *= value;
00104 return result;
00105 }
00107
00109 void operator /= (const T & value)
00110 {
00111 m_Matrix /= value;
00112 }
00113
00115 Self operator/(const T & value)
00116 {
00117 Self result( *this );
00118 result /= value;
00119 return result;
00120 }
00122
00124 inline T & operator()( unsigned int row, unsigned int col )
00125 {
00126 return m_Matrix(row,col);
00127 }
00128
00130 inline const T & operator()( unsigned int row, unsigned int col ) const
00131 { return m_Matrix(row,col); }
00132
00134 inline T * operator[]( unsigned int i )
00135 { return m_Matrix[i]; }
00136
00138 inline const T * operator[]( unsigned int i ) const
00139 { return m_Matrix[i]; }
00140
00142 inline InternalMatrixType & GetVnlMatrix( void )
00143 { return m_Matrix; }
00144
00146 inline const InternalMatrixType & GetVnlMatrix( void ) const
00147 { return m_Matrix; }
00148
00150 inline void SetIdentity( void )
00151 { m_Matrix.set_identity(); }
00152
00154 inline void Fill( const T & value )
00155 { m_Matrix.fill( value ); }
00156
00158 inline const Self & operator=( const vnl_matrix<T> & matrix)
00159 {
00160 m_Matrix = matrix;
00161 return *this;
00162 }
00163
00165 inline Matrix(const vnl_matrix<T> & matrix)
00166 {
00167 this->operator=(matrix);
00168 }
00169
00171 inline bool operator==( const Self & matrix) const
00172 {
00173 bool equal = true;
00174 for( unsigned int r=0; r<NRows; r++)
00175 {
00176 for( unsigned int c=0; c<NColumns; c++ )
00177 {
00178 if (m_Matrix(r,c) != matrix.m_Matrix(r,c))
00179 {
00180 equal = false;
00181 break;
00182 }
00183 }
00184 }
00185 return equal;
00186 }
00188
00189 inline bool operator!=( const Self & matrix) const
00190 {
00191 return !this->operator==(matrix);
00192 }
00193
00194 inline const Self & operator=( const InternalMatrixType & matrix )
00195 {
00196 this->m_Matrix = matrix;
00197 return *this;
00198 }
00199
00201 inline explicit Matrix(const InternalMatrixType & matrix)
00202 {
00203 this->operator=(matrix);
00204 }
00205
00207 inline const Self & operator=( const Self & matrix)
00208 {
00209 m_Matrix = matrix.m_Matrix;
00210 return *this;
00211 }
00212
00214 inline vnl_matrix_fixed<T,NColumns,NRows> GetInverse( void ) const
00215 {
00216 if (vnl_determinant(m_Matrix) == 0.0)
00217 {
00218 itkGenericExceptionMacro(<< "Singular matrix. Determinant is 0.");
00219 }
00220 vnl_matrix<T> temp = vnl_matrix_inverse<T>( m_Matrix );
00221 return temp;
00222 }
00224
00226 inline vnl_matrix_fixed<T,NColumns,NRows> GetTranspose( void ) const
00227 {
00228 return m_Matrix.transpose();
00229 }
00230
00232 Matrix() : m_Matrix(NumericTraits<T>::Zero) {};
00233
00235 Matrix(const Self & matrix) : m_Matrix( matrix.m_Matrix ) {};
00236
00237 private:
00238 InternalMatrixType m_Matrix;
00239
00240 };
00241
00242 template< class T, unsigned int NRows, unsigned int NColumns >
00243 ITK_EXPORT std::ostream& operator<<( std::ostream& os, const Matrix<T,NRows,NColumns> & v)
00244 {
00245 os << v.GetVnlMatrix();
00246 return os;
00247 }
00248
00249
00250 }
00251
00252
00253 #define ITK_TEMPLATE_Matrix(_, EXPORT, x, y) namespace itk { \
00254 _(3(class EXPORT Matrix< ITK_TEMPLATE_3 x >)) \
00255 namespace Templates { typedef Matrix< ITK_TEMPLATE_3 x > Matrix##y; } \
00256 }
00257
00258 #if ITK_TEMPLATE_EXPLICIT
00259 # include "Templates/itkMatrix+-.h"
00260 #endif
00261
00262 #if ITK_TEMPLATE_TXX
00263 # include "itkMatrix.txx"
00264 #endif
00265
00266 #endif
00267