KWStyle - itkSymmetricEigenAnalysis.h

Matrix View

Description

---

1		/*=========================================================================
2
3		Program:   Insight Segmentation & Registration Toolkit
4		Module:    $RCSfile: itkSymmetricEigenAnalysis.h.html,v $
5		Language:  C++
6		Date:      $Date: 2006/01/17 19:15:48 $
7		Version:   $Revision: 1.4 $
8
9		Copyright (c) Insight Software Consortium. All rights reserved.
10		See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11
12		This software is distributed WITHOUT ANY WARRANTY; without even
13		the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14		PURPOSE.  See the above copyright notices for more information.
15
16		=========================================================================*/
17		#ifndef __itkSymmetricEigenAnalysis_h
18		#define __itkSymmetricEigenAnalysis_h
19
20		#include "itkMacro.h"
21
22		namespace itk
23		{
24
25		/** \class SymmetricEigenAnalysis
26		* \brief Find Eigen values of a real 2D symmetric matrix. It
27		* serves as a thread safe alternative to the class:
28		* vnl_symmetric_eigensystem, which uses netlib routines.
29		*
30		* The class is templated over the input matrix, (which is expected to provide
31		* access to its elements with the [][] operator), matrix to store eigen
32	LEN	* values, (must provide write operations on its elements with the [] operator),
33	LEN	* EigenMatrix to store eigen vectors (must provide write access to its elements
34		* with the [][] operator).
35		*
36		* The SetOrderEigenValues() method can be used to order eigen values (and their
37		* corresponding eigen vectors if computed) in ascending order. This is the
38		* default ordering scheme. Eigen vectors and values can be obtained without
39		* ordering by calling SetOrderEigenValues(false)
40		*
41		* The SetOrderEigenMagnitudes() method can be used to order eigen values (and
42	LEN	* their corresponding eigen vectors if computed) by magnitude in ascending order.
43		*
44		* The user of this class is explicitly supposed to set the dimension of the
45		* 2D matrix using the SetDimension() method.
46		*
47		* The class contains routines taken from netlib sources. (www.netlib.org).
48		* netlib/tql1.c
49		* netlib/tql2.c
50		* netlib/tred1.c
51		* netlib/tred2.c
52		*
53		* Reference:
54		*     num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
55		*     wilkinson.
56		*     handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
57		*/
58
59		template < typename TMatrix, typename TVector, typename TEigenMatrix=TMatrix >
60		class SymmetricEigenAnalysis
61		{
62		public:
63		typedef enum {
64		OrderByValue=1,
65		OrderByMagnitude,
66		DoNotOrder
67	IND	**}EigenValueOrderType;
68
69		SymmetricEigenAnalysis():
70	IND	******m_Dimension(0),
71	IND	******m_Order(0),
72	IND	******m_OrderEigenValues(OrderByValue)
73	SEM	{} ;
74
75		SymmetricEigenAnalysis( const unsigned int dimension ):
76	IND	******m_Dimension(dimension),
77	IND	******m_Order(dimension),
78	IND	******m_OrderEigenValues(OrderByValue)
79	SEM	{} ;
80
81		~SymmetricEigenAnalysis() {};
82
83		typedef TMatrix      MatrixType;
84		typedef TEigenMatrix EigenMatrixType;
85		typedef TVector      VectorType;
86
87
88	EML
89		/** Compute Eigen values of A
90		* A is any type that overloads the [][] operator and contains the
91		* symmetric matrix. In practice only the upper triangle of the
92		* matrix will be accessed. (Both itk::Matrix and vnl_matrix
93		* overload [][] operator.)
94		*
95		* 'EigenValues' is any type that overloads the []   operator and will contain
96		* the eigen values.
97		*
98		* No size checking is performed. A is expected to be a square matrix of size
99		* m_Dimension.  'EigenValues' is expected to be of length m_Dimension.
100		* The matrix is not checked to see if it is symmetric.
101		*/
102		unsigned int ComputeEigenValues(
103		const TMatrix  & A,
104		TVector        & EigenValues) const;
105
106		/** Compute Eigen values and vectors of A
107		* A is any type that overloads the [][] operator and contains the
108		* symmetric matrix. In practice only the upper triangle of the
109		* matrix will be accessed. (Both itk::Matrix and vnl_matrix
110		* overload [][] operator.)
111		*
112		* 'EigenValues' is any type that overloads the []   operator and will contain
113		* the eigen values.
114		*
115		* 'EigenVectors' is any type that provides access to its elements with the
116		* [][] operator. It is expected be of size m_Dimension * m_Dimension.
117		*
118		* No size checking is performed. A is expected to be a square matrix of size
119		* m_Dimension.  'EigenValues' is expected to be of length m_Dimension.
120		* The matrix is not checked to see if it is symmetric.
121		*
122	LEN	* Each row of the matrix 'EigenVectors' represents an eigen vector. (unlike MATLAB
123		* where the columns of the [EigenVectors, EigenValues] = eig(A) contains the
124		* eigenvectors).
125		*/
126		unsigned int ComputeEigenValuesAndVectors(
127		const TMatrix  & A,
128		TVector        & EigenValues,
129		TEigenMatrix   & EigenVectors ) const;
130
131
132		/** Matrix order. Defaults to matrix dimension if not set **/
133		void SetOrder(const unsigned int n)
134		{
135		m_Order = n;
136		}
137
138		/** Get the Matrix order. Will be 0 unless explicitly set, or unless a
139		* call to SetDimension has been made in which case it will be the
140		* matrix dimension. */
141		unsigned int GetOrder() const { return m_Order; }
142
143		/** Set/Get methods to order the eigen values in ascending order.
144		* This is the default. ie lambda_1 < lambda_2 < ....
145		*/
146		void SetOrderEigenValues( const bool b )
147		{
148		if (b) { m_OrderEigenValues = OrderByValue;     }
149		else   { m_OrderEigenValues = DoNotOrder;       }
150		}
151	LEN	bool GetOrderEigenValues() const { return (m_OrderEigenValues == OrderByValue); }
152
153		/** Set/Get methods to order the eigen value magnitudes in ascending order.
154		* In other words, |lambda_1| < |lambda_2| < .....
155		*/
156		void SetOrderEigenMagnitudes( const bool b )
157		{
158		if (b) { m_OrderEigenValues = OrderByMagnitude; }
159		else   { m_OrderEigenValues = DoNotOrder;       }
160		}
161	LEN	bool GetOrderEigenMagnitudes() const { return (m_OrderEigenValues == OrderByMagnitude); }
162
163		/** Set the dimension of the input matrix A. A is a square matrix of
164		* size m_Dimension. */
165		void SetDimension( const unsigned int n )
166		{
167		m_Dimension = n;
168		if (m_Order == 0 )
169		{
170		m_Order = m_Dimension;
171		}
172		}
173
174		/** Get Matrix dimension, Will be 0 unless explicitly set by a
175		* call to SetDimension. */
176		unsigned int GetDimension() const { return m_Dimension; }
177
178
179		private:
180		unsigned int m_Dimension;
181		unsigned int m_Order;
182		EigenValueOrderType         m_OrderEigenValues;
183
184
185		/** Reduces a real symmetric matrix to a symmetric tridiagonal matrix using
186		*  orthogonal similarity transformations.
187		*  'inputMatrix' contains the real symmetric input matrix. Only the lower
188		*  triangle of the matrix need be supplied. The upper triangle is unaltered.
189		*  'd' contains the diagonal elements of the tridiagonal matrix.
190		*  'e' contains the subdiagonal elements of the tridiagonal matrix in its
191		*  last n-1 positions.  e(1) is set to zero.
192		*  'e2' contains the squares of the corresponding elements of e.
193		*  'e2' may coincide with e if the squares are not needed.
194		*  questions and comments should be directed to burton s. garbow.
195		*  mathematics and computer science div, argonne national laboratory
196		*     this version dated august 1983.
197		*
198		*  Function Adapted from netlib/tred1.c.
199		*  [Changed: remove static vars, enforce const correctness.
200		*            Use vnl routines as necessary].
201		*  Reference:
202		*  num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson.
203		*    handbook for auto. comp., vol.ii-linear algebra, 212-226(1971).    */
204		void ReduceToTridiagonalMatrix(double *inputMatrix, VectorType &d,
205		double *e, double *e2) const;
206
207
208
209		/** Reduces a real symmetric matrix to a symmetric tridiagonal matrix using
210		*  and accumulating orthogonal similarity transformations.
211		*  'inputMatrix' contains the real symmetric input matrix. Only the lower
212		*  triangle of the matrix need be supplied. The upper triangle is unaltered.
213		*  'diagonalElements' will contains the diagonal elements of the tridiagonal
214		*  matrix.
215	LEN	*  'subDiagonalElements' will contain the subdiagonal elements of the tridiagonal
216	LEN	*  matrix in its last n-1 positions.  subDiagonalElements(1) is set to zero.
217		*  'transformMatrix' contains the orthogonal transformation matrix produced
218		*  in the reduction.
219		*
220		*  questions and comments should be directed to burton s. garbow.
221		*  mathematics and computer science div, argonne national laboratory
222		*     this version dated august 1983.
223		*
224		*  Function Adapted from netlib/tred2.c.
225		*  [Changed: remove static vars, enforce const correctness.
226		*            Use vnl routines as necessary].
227		*  Reference:
228		*  num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson.
229		*    handbook for auto. comp., vol.ii-linear algebra, 212-226(1971).    */
230		void ReduceToTridiagonalMatrixAndGetTransformation(
231		double *inputMatrix, VectorType &diagonalElements,
232		double *subDiagonalElements, double *transformMatrix) const;
233
234
235
236		/* Finds the eigenvalues of a symmetric tridiagonal matrix by the ql method.
237		*
238		* On input:
239		* 'd' contains the diagonal elements of the input matrix.
240		* 'e' contains the subdiagonal elements of the input matrix
241		* in its last n-1 positions.  e(1) is arbitrary.
242		* On Output:
243		* 'd' contains the eigenvalues.
244		* 'e' has been destroyed.
245		*
246		* Returns:
247		*          zero       for normal return,
248		*          j          if the j-th eigenvalue has not been
249		*                     determined after 30 iterations.
250		*
251		*
252		* Reference
253		*     this subroutine is a translation of the algol procedure tql1,
254		*     num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
255		*     wilkinson.
256		*     handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
257		*
258		*     Questions and comments should be directed to burton s. garbow,
259		*     mathematics and computer science div, argonne national laboratory
260		*     this version dated august 1983.
261		*
262		*  Function Adapted from netlib/tql1.c.
263		*  [Changed: remove static vars, enforce const correctness.
264		*            Use vnl routines as necessary]                      */
265		unsigned int ComputeEigenValuesUsingQL(
266		VectorType &d, double *e) const;
267
268
269		/* Finds the eigenvalues and eigenvectors of a symmetric tridiagonal matrix
270		* by the ql method.
271		*
272		* On input:
273		* 'd' contains the diagonal elements of the input matrix.
274		* 'e' contains the subdiagonal elements of the input matrix
275		* in its last n-1 positions.  e(1) is arbitrary.
276		* 'z' contains the transformation matrix produced in the reduction by
277		* ReduceToTridiagonalMatrixAndGetTransformation(), if performed. If the
278		* eigenvectors of the tridiagonal matrix are desired, z must contain
279		* the identity matrix.
280
281		* On Output:
282		* 'd' contains the eigenvalues.
283		* 'e' has been destroyed.
284		* 'z' contains orthonormal eigenvectors of the symmetric tridiagonal
285		* (or full) matrix.
286		*
287		* Returns:
288		*          zero       for normal return,
289		*          j          if the j-th eigenvalue has not been
290		*                     determined after 1000 iterations.
291		*
292		* Reference
293		*     this subroutine is a translation of the algol procedure tql1,
294		*     num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
295		*     wilkinson.
296		*     handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
297		*
298		*     Questions and comments should be directed to burton s. garbow,
299		*     mathematics and computer science div, argonne national laboratory
300		*     this version dated august 1983.
301		*
302		*  Function Adapted from netlib/tql2.c.
303		*  [Changed: remove static vars, enforce const correctness.
304		*            Use vnl routines as necessary]
305		*/
306		unsigned int ComputeEigenValuesAndVectorsUsingQL(
307		VectorType &d, double *e, double *z) const;
308
309		};
310
311		template< typename TMatrix, typename TVector, typename TEigenMatrix >
312		std::ostream & operator<<(std::ostream& os,
313		const SymmetricEigenAnalysis< TMatrix, TVector, TEigenMatrix > &s)
314		{
315		os << "[ClassType: SymmetricEigenAnalysis]" << std::endl;
316		os << "  Dimension : " << s.GetDimension() << std::endl;
317		os << "  Order : " << s.GetOrder() << std::endl;
318		os << "  OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
319		os << "  OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
320		return os;
321		}
322
323		} // end namespace itk
324
325
326		#ifndef ITK_MANUAL_INSTANTIATION
327		#include "itkSymmetricEigenAnalysis.txx"
328		#endif
329
330		#endif
331
332	EOF

---

Generated by KWStyle 1.0b on Tuesday January,17 at 02:14:48PM

© Kitware Inc.