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