00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #ifndef __itkSymmetricEigenAnalysis_h
00018 #define __itkSymmetricEigenAnalysis_h
00019
00020 #include "itkMacro.h"
00021
00022 namespace itk
00023 {
00024
00059 template < typename TMatrix, typename TVector, typename TEigenMatrix=TMatrix >
00060 class SymmetricEigenAnalysis
00061 {
00062 public:
00063 typedef enum {
00064 OrderByValue=1,
00065 OrderByMagnitude,
00066 DoNotOrder
00067 }EigenValueOrderType;
00068
00069 SymmetricEigenAnalysis():
00070 m_Dimension(0),
00071 m_Order(0),
00072 m_OrderEigenValues(OrderByValue)
00073 {};
00074
00075 SymmetricEigenAnalysis( const unsigned int dimension ):
00076 m_Dimension(dimension),
00077 m_Order(dimension),
00078 m_OrderEigenValues(OrderByValue)
00079 {};
00080
00081 ~SymmetricEigenAnalysis() {};
00082
00083 typedef TMatrix MatrixType;
00084 typedef TEigenMatrix EigenMatrixType;
00085 typedef TVector VectorType;
00086
00100 unsigned int ComputeEigenValues(
00101 const TMatrix & A,
00102 TVector & EigenValues) const;
00103
00124 unsigned int ComputeEigenValuesAndVectors(
00125 const TMatrix & A,
00126 TVector & EigenValues,
00127 TEigenMatrix & EigenVectors ) const;
00128
00129
00131 void SetOrder(const unsigned int n)
00132 {
00133 m_Order = n;
00134 }
00135
00139 unsigned int GetOrder() const { return m_Order; }
00140
00144 void SetOrderEigenValues( const bool b )
00145 {
00146 if (b) { m_OrderEigenValues = OrderByValue; }
00147 else { m_OrderEigenValues = DoNotOrder; }
00148 }
00149 bool GetOrderEigenValues() const { return (m_OrderEigenValues == OrderByValue); }
00151
00155 void SetOrderEigenMagnitudes( const bool b )
00156 {
00157 if (b) { m_OrderEigenValues = OrderByMagnitude; }
00158 else { m_OrderEigenValues = DoNotOrder; }
00159 }
00160 bool GetOrderEigenMagnitudes() const { return (m_OrderEigenValues == OrderByMagnitude); }
00162
00165 void SetDimension( const unsigned int n )
00166 {
00167 m_Dimension = n;
00168 if (m_Order == 0 )
00169 {
00170 m_Order = m_Dimension;
00171 }
00172 }
00174
00177 unsigned int GetDimension() const { return m_Dimension; }
00178
00179
00180 private:
00181 unsigned int m_Dimension;
00182 unsigned int m_Order;
00183 EigenValueOrderType m_OrderEigenValues;
00184
00185
00205 void ReduceToTridiagonalMatrix(double *inputMatrix, VectorType &d,
00206 double *e, double *e2) const;
00207
00229 void ReduceToTridiagonalMatrixAndGetTransformation(
00230 double *inputMatrix, VectorType &diagonalElements,
00231 double *subDiagonalElements, double *transformMatrix) const;
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262 unsigned int ComputeEigenValuesUsingQL(
00263 VectorType &d, double *e) const;
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303 unsigned int ComputeEigenValuesAndVectorsUsingQL(
00304 VectorType &d, double *e, double *z) const;
00305
00306 };
00307
00308 template< typename TMatrix, typename TVector, typename TEigenMatrix >
00309 std::ostream & operator<<(std::ostream& os,
00310 const SymmetricEigenAnalysis< TMatrix, TVector, TEigenMatrix > &s)
00311 {
00312 os << "[ClassType: SymmetricEigenAnalysis]" << std::endl;
00313 os << " Dimension : " << s.GetDimension() << std::endl;
00314 os << " Order : " << s.GetOrder() << std::endl;
00315 os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
00316 os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
00317 return os;
00318 }
00319
00320 }
00321
00322
00323 #ifndef ITK_MANUAL_INSTANTIATION
00324 #include "itkSymmetricEigenAnalysis.txx"
00325 #endif
00326
00327 #endif
00328