KWStyle - itkSymmetricEigenAnalysis.txx

Matrix View

Description

---

1		/*=========================================================================
2
3		Program:   Insight Segmentation & Registration Toolkit
4		Module:    $RCSfile: itkSymmetricEigenAnalysis.txx.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_txx
18		#define __itkSymmetricEigenAnalysis_txx
19
20		#include "itkSymmetricEigenAnalysis.h"
21		#include "vnl/vnl_math.h"
22
23		namespace itk
24		{
25
26		template< class TMatrix, class TVector, class TEigenMatrix >
27		unsigned int
28		SymmetricEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
29		ComputeEigenValues(const TMatrix  & A,
30		TVector  & D) const
31		{
32		double *workArea1 = new double[ m_Dimension ];
33
34		// Copy the input matrix
35		double *inputMatrix = new double [ m_Dimension * m_Dimension ];
36
37		unsigned int k = 0;
38
39		for( unsigned int row=0; row < m_Dimension; row++ )
40		{
41		for( unsigned int col=0; col < m_Dimension; col++ )
42		{
43		inputMatrix[k++] = A(row,col);
44		}
45		}
46
47		ReduceToTridiagonalMatrix( inputMatrix, D, workArea1, workArea1 );
48		const unsigned int eigenErrIndex =
49		ComputeEigenValuesUsingQL( D, workArea1 );
50
51		delete[] workArea1;
52		delete[] inputMatrix;
53
54		return eigenErrIndex; //index of eigen value that could not be computed
55		}
56
57
58		template< class TMatrix, class TVector, class TEigenMatrix >
59		unsigned int
60		SymmetricEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
61		ComputeEigenValuesAndVectors(
62		const TMatrix  & A,
63		TVector        & EigenValues,
64		TEigenMatrix   & EigenVectors ) const
65		{
66		double *workArea1 = new double[ m_Dimension ];
67		double *workArea2 = new double [ m_Dimension * m_Dimension ];
68
69		// Copy the input matrix
70		double *inputMatrix = new double [ m_Dimension * m_Dimension ];
71
72		unsigned int k = 0;
73
74		for( unsigned int row=0; row < m_Dimension; row++ )
75		{
76		for( unsigned int col=0; col < m_Dimension; col++ )
77		{
78		inputMatrix[k++] = A(row,col);
79		}
80		}
81
82		ReduceToTridiagonalMatrixAndGetTransformation(
83		inputMatrix, EigenValues, workArea1, workArea2 );
84		const unsigned int eigenErrIndex =
85		ComputeEigenValuesAndVectorsUsingQL( EigenValues, workArea1, workArea2 );
86
87		// Copy eigenVectors
88		k = 0;
89		for( unsigned int row=0; row < m_Dimension; row++ )
90		{
91		for( unsigned int col=0; col < m_Dimension; col++ )
92		{
93		EigenVectors[row][col] = workArea2[k++];
94		}
95		}
96
97		delete[] workArea2;
98		delete[] workArea1;
99		delete[] inputMatrix;
100
101		return eigenErrIndex; //index of eigen value that could not be computed
102		}
103
104
105	EML
106	EML
107		template< class TMatrix, class TVector, class TEigenMatrix >
108		void
109		SymmetricEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
110		ReduceToTridiagonalMatrix(double * a, VectorType &d,
111		double *e, double *e2) const
112		{
113		double d__1;
114
115		/* Local variables */
116		double f, g, h;
117		int i, j, k, l;
118		double scale;
119
120		for (i = 0; i < static_cast< int >(m_Order); ++i)
121		{
122		d[i] = a[m_Order-1 + i * m_Dimension];
123		a[m_Order-1 + i * m_Dimension] = a[i + i * m_Dimension];
124		}
125
126
127		for (i = m_Order-1; i >= 0; --i)
128		{
129		l = i - 1;
130		h = 0.;
131		scale = 0.;
132
133		/*     .......... scale row (algol tol then not needed) .......... */
134		for (k = 0; k <= l; ++k)
135		{
136		scale += vnl_math_abs(d[k]);
137		}
138		if (scale == 0.)
139		{
140		for (j = 0; j <= l; ++j)
141		{
142		d[j] = a[l + j * m_Dimension];
143		a[l + j * m_Dimension] = a[i + j * m_Dimension];
144		a[i + j * m_Dimension] = 0.;
145		}
146		e[i] = 0.;
147		e2[i] = 0.;
148		continue;
149		}
150		for (k = 0; k <= l; ++k)
151		{
152		d[k] /= scale;
153		h += d[k] * d[k];
154		}
155
156		e2[i] = scale * scale * h;
157		f = d[l];
158		d__1 = vcl_sqrt(h);
159		g = (-1.0) * vnl_math_sgn0(f) * vnl_math_abs(d__1);
160		e[i] = scale * g;
161		h -= f * g;
162		d[l] = f - g;
163		if (l != 0)
164		{
165
166		/*     .......... form a*u .......... */
167		for (j = 0; j <= l; ++j)
168		{
169		e[j] = 0.;
170		}
171
172		for (j = 0; j <= l; ++j)
173		{
174		f = d[j];
175		g = e[j] + a[j + j * m_Dimension] * f;
176
177		for (k = j+1; k <= l; ++k)
178		{
179		g += a[k + j * m_Dimension] * d[k];
180		e[k] += a[k + j * m_Dimension] * f;
181		}
182		e[j] = g;
183		}
184
185		/*     .......... form p .......... */
186		f = 0.;
187
188		for (j = 0; j <= l; ++j)
189		{
190		e[j] /= h;
191		f += e[j] * d[j];
192		}
193
194		h = f / (h + h);
195		/*     .......... form q .......... */
196		for (j = 0; j <= l; ++j)
197		{
198		e[j] -= h * d[j];
199		}
200
201		/*     .......... form reduced a .......... */
202		for (j = 0; j <= l; ++j)
203		{
204		f = d[j];
205		g = e[j];
206
207		for (k = j; k <= l; ++k)
208		{
209		a[k + j * m_Dimension] = a[k + j * m_Dimension] - f * e[k] - g * d[k];
210		}
211		}
212		}
213
214		for (j = 0; j <= l; ++j)
215		{
216		f = d[j];
217		d[j] = a[l + j * m_Dimension];
218		a[l + j * m_Dimension] = a[i + j * m_Dimension];
219		a[i + j * m_Dimension] = f * scale;
220		}
221		}
222		}
223
224
225	EML
226		template< class TMatrix, class TVector, class TEigenMatrix >
227		void
228		SymmetricEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
229		ReduceToTridiagonalMatrixAndGetTransformation(double * a, VectorType &d,
230		double *e, double *z) const
231		{
232		double d__1;
233
234		/* Local variables */
235		double f, g, h;
236		unsigned int i, j, k, l;
237		double scale, hh;
238
239		for (i = 0; i < m_Order; ++i)
240		{
241		for (j = i; j < m_Order; ++j)
242		{
243		z[j + i * m_Dimension] = a[j + i * m_Dimension];
244		}
245		d[i] = a[m_Order -1 + i * m_Dimension];
246		}
247
248		for (i = m_Order-1; i > 0; --i)
249		{
250		l = i - 1;
251		h = 0.0;
252		scale = 0.0;
253
254		/*     .......... scale row (algol tol then not needed) .......... */
255		for (k = 0; k <= l; ++k)
256		{
257		scale += vnl_math_abs(d[k]);
258		}
259		if (scale == 0.0)
260		{
261		e[i] = d[l];
262
263		for (j = 0; j <= l; ++j)
264		{
265		d[j] = z[l + j * m_Dimension];
266		z[i + j * m_Dimension] = 0.0;
267		z[j + i * m_Dimension] = 0.0;
268		}
269		}
270		else
271		{
272		for (k = 0; k <= l; ++k)
273		{
274		d[k] /= scale;
275		h += d[k] * d[k];
276		}
277
278		f = d[l];
279		d__1 = vcl_sqrt(h);
280		g = (-1.0) * vnl_math_sgn0(f) * vnl_math_abs(d__1);
281		e[i] = scale * g;
282		h -= f * g;
283		d[l] = f - g;
284
285		/*     .......... form a*u .......... */
286		for (j = 0; j <= l; ++j)
287		{
288		e[j] = 0.0;
289		}
290
291		for (j = 0; j <= l; ++j)
292		{
293		f = d[j];
294		z[j + i * m_Dimension] = f;
295		g = e[j] + z[j + j * m_Dimension] * f;
296
297		for (k = j+1; k <= l; ++k)
298		{
299		g += z[k + j * m_Dimension] * d[k];
300		e[k] += z[k + j * m_Dimension] * f;
301		}
302
303		e[j] = g;
304		}
305
306		/*     .......... form p .......... */
307		f = 0.0;
308
309		for (j = 0; j <= l; ++j)
310		{
311		e[j] /= h;
312		f += e[j] * d[j];
313		}
314
315		hh = f / (h + h);
316
317		/*     .......... form q .......... */
318		for (j = 0; j <= l; ++j)
319		{
320		e[j] -= hh * d[j];
321		}
322
323		/*     .......... form reduced a .......... */
324		for (j = 0; j <= l; ++j)
325		{
326		f = d[j];
327		g = e[j];
328
329		for (k = j; k <= l; ++k)
330		{
331		z[k + j * m_Dimension] = z[k + j * m_Dimension] - f * e[k] - g * d[k];
332		}
333
334		d[j] = z[l + j * m_Dimension];
335		z[i + j * m_Dimension] = 0.0;
336		}
337		}
338
339		d[i] = h;
340		}
341
342		/*     .......... accumulation of transformation matrices .......... */
343		for (i = 1; i < m_Order; ++i)
344		{
345		l = i - 1;
346		z[m_Order-1 + l * m_Dimension] = z[l + l * m_Dimension];
347		z[l + l * m_Dimension] = 1.0;
348		h = d[i];
349		if (h != 0.0)
350		{
351		for (k = 0; k <= l; ++k)
352		{
353		d[k] = z[k + i * m_Dimension] / h;
354		}
355
356		for (j = 0; j <= l; ++j)
357		{
358		g = 0.0;
359
360		for (k = 0; k <= l; ++k)
361		{
362		g += z[k + i * m_Dimension] * z[k + j * m_Dimension];
363		}
364
365		for (k = 0; k <= l; ++k)
366		{
367		z[k + j * m_Dimension] -= g * d[k];
368		}
369		}
370		}
371
372		for (k = 0; k <= l; ++k)
373		{
374		z[k + i * m_Dimension] = 0.0;
375		}
376
377		}
378
379		for (i = 0; i < m_Order; ++i)
380		{
381		d[i] = z[m_Order-1 + i * m_Dimension];
382		z[m_Order-1 + i * m_Dimension] = 0.0;
383		}
384
385		z[m_Order-1 + (m_Order-1) * m_Dimension] = 1.0;
386		e[0] = 0.0;
387
388		}
389
390
391	EML
392		template< class TMatrix, class TVector, class TEigenMatrix >
393		unsigned int
394		SymmetricEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
395		ComputeEigenValuesUsingQL(VectorType &d, double *e) const
396		{
397
398		const double c_b10 = 1.0;
399
400		/* Local variables */
401		double c, f, g, h;
402		unsigned int i, j, l, m;
403		double p, r, s, c2, c3=0.0;
404		double s2=0.0;
405		double dl1, el1;
406		double tst1, tst2;
407
408		unsigned int ierr = 0;
409		if (m_Order == 1) {
410	IND	*******return 1;
411	IND	**}
412
413		for (i = 1; i < m_Order; ++i) {
414		e[i-1] = e[i];
415		}
416
417		f = 0.;
418		tst1 = 0.;
419		e[m_Order-1] = 0.;
420
421		for (l = 0; l < m_Order; ++l)
422		{
423		j = 0;
424		h = vnl_math_abs(d[l]) + vnl_math_abs(e[l]);
425		if (tst1 < h)
426		{
427		tst1 = h;
428		}
429		/*     .......... look for small sub-diagonal element .......... */
430		for (m = l; m < m_Order; ++m)
431		{
432		tst2 = tst1 + vnl_math_abs(e[m]);
433		if (tst2 == tst1)
434		{
435		break;
436		}
437		/*     .......... e(n) is always zero, so there is no exit */
438		/*                through the bottom of the loop .......... */
439		}
440
441		if (m != l)
442		{
443
444		do
445		{
446		if (j == 30)
447		{
448		/*     .......... set error -- no convergence to an */
449		/*                eigenvalue after 30 iterations .......... */
450		ierr = l+1;
451		return ierr;
452		}
453		++j;
454		/*     .......... form shift .......... */
455		g = d[l];
456		p = (d[l+1] - g) / (e[l] * 2.);
457		r = vnl_math_hypot(p, c_b10);
458		d[l] = e[l] / (p + vnl_math_sgn0(p) * vnl_math_abs(r));
459		d[l+1] = e[l] * (p + vnl_math_sgn0(p) * vnl_math_abs(r));
460		dl1 = d[l+1];
461		h = g - d[l];
462
463		for (i = l+2; i < m_Order; ++i)
464		{
465		d[i] -= h;
466		}
467
468		f += h;
469		/*     .......... ql transformation .......... */
470		p = d[m];
471		c = 1.;
472		c2 = c;
473		el1 = e[l+1];
474		s = 0.;
475		for (i = m-1; i >= l; --i)
476		{
477		c3 = c2;
478		c2 = c;
479		s2 = s;
480		g = c * e[i];
481		h = c * p;
482		r = vnl_math_hypot(p, e[i]);
483		e[i+1] = s * r;
484		s = e[i] / r;
485		c = p / r;
486		p = c * d[i] - s * g;
487		d[i+1] = h + s * (c * g + s * d[i]);
488		if( i == l )
489		{
490		break;
491		}
492		}
493
494		p = -s * s2 * c3 * el1 * e[l] / dl1;
495		e[l] = s * p;
496		d[l] = c * p;
497		tst2 = tst1 + vnl_math_abs(e[l]);
498		} while (tst2 > tst1);
499		}
500
501		p = d[l] + f;
502
503		if( m_OrderEigenValues == OrderByValue )
504		{
505		// Order by value
506		for (i = l; i > 0; --i)
507		{
508		if (p >= d[i-1])
509	IND	**********break;
510		d[i] = d[i-1];
511		}
512		d[i] = p;
513		}
514		else if( m_OrderEigenValues == OrderByMagnitude )
515		{
516		// Order by magnitude.. make eigen values positive
517		for (i = l; i > 0; --i)
518		{
519		if (vnl_math_abs(p) >= vnl_math_abs(d[i-1]))
520	IND	**********break;
521		d[i] = vnl_math_abs(d[i-1]);
522		}
523		d[i] = vnl_math_abs(p);
524		}
525		else
526		{
527		d[l] = p;
528		}
529		}
530
531		return ierr;    //ierr'th eigen value that couldn't be computed
532
533		}
534
535
536	EML
537		template< class TMatrix, class TVector, class TEigenMatrix >
538		unsigned int
539		SymmetricEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
540		ComputeEigenValuesAndVectorsUsingQL(VectorType &d, double *e, double *z) const
541		{
542
543		const double c_b10 = 1.0;
544
545		/* Local variables */
546		double c, f, g, h;
547		unsigned int i, j, k, l, m;
548		double p, r, s, c2, c3=0.0;
549		double s2=0.0;
550		double dl1, el1;
551		double tst1, tst2;
552
553		unsigned int ierr = 0;
554		if (m_Order == 1)
555		{
556		return 1;
557		}
558
559		for (i = 1; i < m_Order; ++i)
560		{
561		e[i - 1] = e[i];
562		}
563
564		f = 0.0;
565		tst1 = 0.;
566		e[m_Order-1] = 0.;
567
568		for (l = 0; l < m_Order; ++l)
569		{
570		j = 0;
571		h = vnl_math_abs(d[l]) + vnl_math_abs(e[l]);
572		if (tst1 < h)
573		{
574		tst1 = h;
575		}
576
577		/*     .......... look for small sub-diagonal element .......... */
578		for (m = l; m < m_Order; ++m)
579		{
580		tst2 = tst1 + vnl_math_abs(e[m]);
581		if (tst2 == tst1)
582		{
583		break;
584		}
585
586		/*     .......... e(n) is always zero, so there is no exit */
587		/*                through the bottom of the loop .......... */
588		}
589
590		if (m != l)
591		{
592		do
593		{
594
595		if (j == 1000)
596		{
597		/*     .......... set error -- no convergence to an */
598		/*                eigenvalue after 1000 iterations .......... */
599		ierr = l+1;
600		return ierr;
601		}
602		++j;
603		/*     .......... form shift .......... */
604		g = d[l];
605		p = (d[l+1] - g) / (e[l] * 2.);
606		r = vnl_math_hypot(p, c_b10);
607		d[l] = e[l] / (p + vnl_math_sgn0(p) * vnl_math_abs(r));
608		d[l+1] = e[l] * (p + vnl_math_sgn0(p) * vnl_math_abs(r));
609		dl1 = d[l+1];
610		h = g - d[l];
611
612		for (i = l+2; i < m_Order; ++i)
613		{
614		d[i] -= h;
615		}
616
617		f += h;
618		/*     .......... ql transformation .......... */
619		p = d[m];
620		c = 1.0;
621		c2 = c;
622		el1 = e[l+1];
623		s = 0.;
624
625		for (i = m-1; i >= l; --i)
626		{
627		c3 = c2;
628		c2 = c;
629		s2 = s;
630		g = c * e[i];
631		h = c * p;
632		r = vnl_math_hypot(p, e[i]);
633		e[i + 1] = s * r;
634		s = e[i] / r;
635		c = p / r;
636		p = c * d[i] - s * g;
637		d[i + 1] = h + s * (c * g + s * d[i]);
638
639		/*     .......... form vector .......... */
640		for (k = 0; k < m_Order; ++k)
641		{
642		h = z[k + (i + 1) * m_Dimension];
643		z[k + (i + 1) * m_Dimension] = s * z[k + i * m_Dimension] + c * h;
644		z[k + i * m_Dimension] = c * z[k + i * m_Dimension] - s * h;
645		}
646		if( i == l )
647		{
648		break;
649		}
650		}
651
652		p = -s * s2 * c3 * el1 * e[l] / dl1;
653		e[l] = s * p;
654		d[l] = c * p;
655		tst2 = tst1 + vnl_math_abs(e[l]);
656		} while (tst2 > tst1);
657
658		}
659
660		d[l] += f;
661		}
662
663		/*     .......... order eigenvalues and eigenvectors .......... */
664		if( m_OrderEigenValues == OrderByValue )
665		{
666		// Order by value
667		for (i = 0; i < m_Order-1; ++i)
668		{
669		k = i;
670		p = d[i];
671
672		for (j = i+1; j < m_Order; ++j)
673		{
674		if (d[j] >= p)
675		{
676		continue;
677		}
678		k = j;
679		p = d[j];
680	IND	******}
681
682		if (k == i)
683		{
684		continue;
685		}
686		d[k] = d[i];
687		d[i] = p;
688
689		for (j = 0; j < m_Order; ++j)
690		{
691		p = z[j + i * m_Dimension];
692		z[j + i * m_Dimension] = z[j + k * m_Dimension];
693		z[j + k * m_Dimension] = p;
694		}
695		}
696		}
697		else if( m_OrderEigenValues == OrderByMagnitude )
698		{
699		// Order by magnitude
700		for (i = 0; i < m_Order-1; ++i)
701		{
702		k = i;
703		p = d[i];
704
705		for (j = i+1; j < m_Order; ++j)
706		{
707		if (vnl_math_abs(d[j]) >= vnl_math_abs(p))
708		{
709		continue;
710		}
711		k = j;
712		p = d[j];
713	IND	******}
714
715		if (k == i)
716		{
717		continue;
718		}
719		d[k] = vnl_math_abs(d[i]);
720		d[i] = vnl_math_abs(p);
721
722		for (j = 0; j < m_Order; ++j)
723		{
724		p = z[j + i * m_Dimension];
725		z[j + i * m_Dimension] = z[j + k * m_Dimension];
726		z[j + k * m_Dimension] = p;
727		}
728		}
729		}
730
731
732		return ierr;
733		}
734
735
736	EML
737		}  // end namespace itk
738
739		#endif
740
741	EOF

---

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

© Kitware Inc.