37#ifndef VIGRA_EIGENSYSTEM_HXX
38#define VIGRA_EIGENSYSTEM_HXX
43#include "array_vector.hxx"
44#include "polynomial.hxx"
57template <
class T,
class C1,
class C2>
63 "housholderTridiagonalization(): matrix must be square.");
65 "housholderTridiagonalization(): matrix size mismatch.");
77 for(
int i = n-1;
i > 0; --
i)
83 for(
int k = 0;
k <
i; ++
k)
85 scale = scale + abs(d(
k));
90 for(
int j = 0;
j <
i; ++
j)
101 for(
int k = 0;
k <
i; ++
k)
107 T g = VIGRA_CSTD::sqrt(
h);
114 for(
int j = 0;
j <
i; ++
j)
121 for(
int j = 0;
j <
i; ++
j)
125 g = e(
j) + a(
j,
j) * f;
126 for(
int k =
j+1;
k <=
i-1; ++
k)
134 for(
int j = 0;
j <
i; ++
j)
140 for(
int j = 0;
j <
i; ++
j)
144 for(
int j = 0;
j <
i; ++
j)
148 for(
int k =
j;
k <=
i-1; ++
k)
150 a(
k,
j) -= (f * e(
k) + g * d(
k));
170 d(
k) = a(
k,
i+1) /
h;
177 g += a(
k,
i+1) * a(
k,
j);
201template <
class T,
class C1,
class C2>
207 "tridiagonalMatrixEigensystem(): matrix must be square.");
209 "tridiagonalMatrixEigensystem(): matrix size mismatch.");
221 T
eps = VIGRA_CSTD::pow(2.0,-52.0);
251 T p = (d(
l+1) - g) / (2.0 * e(
l));
257 d(
l) = e(
l) / (p + r);
258 d(
l+1) = e(
l) * (p + r);
287 p = c * d(
i) - s * g;
288 d(
i+1) =
h + s * (c * g + s * d(
i));
295 z(
k,
i+1) = s * z(
k,
i) + c *
h;
296 z(
k,
i) = c * z(
k,
i) - s *
h;
328 std::swap(d(
k), d(
i));
331 std::swap(z(
j,
i), z(
j,
k));
340template <
class T,
class C1,
class C2>
360 scale = scale + abs(
H(
i,
m-1));
373 T g = VIGRA_CSTD::sqrt(
h);
384 for(
int j =
m;
j < n; ++
j)
418 for(
int i = 0;
i < n; ++
i)
420 for(
int j = 0;
j < n; ++
j)
422 V(
i,
j) = (
i ==
j ? 1.0 : 0.0);
459 if(abs(
yr) > abs(
yi))
475template <
class T,
class C>
477 T & p, T &
q, T & r, T & s, T & w, T & x, T & y, T & z)
485 p = (r * s - w) /
H(
m+1,
m) +
H(
m,
m+1);
486 q =
H(
m+1,
m+1) - z - r - s;
488 s = abs(p) + abs(
q) + abs(r);
496 if(abs(
H(
m,
m-1)) * (abs(
q) + abs(r)) <
497 eps * (abs(p) * (abs(
H(
m-1,
m-1)) + abs(z) +
511template <
class T,
class C1,
class C2,
class C3>
529 T
eps = VIGRA_CSTD::pow(2.0,
sizeof(T) ==
sizeof(
float)
533 T p=0,
q=0,r=0,s=0,z=0,t,w,x,y;
545 s = abs(
H(
l-1,
l-1)) + abs(
H(
l,
l));
550 if(abs(
H(
l,
l-1)) <
eps * s)
572 w =
H(n, n-1) *
H(n-1, n);
573 p = (
H(n-1, n-1) -
H(n, n)) / 2.0;
575 z = VIGRA_CSTD::sqrt(abs(
q));
604 r = VIGRA_CSTD::sqrt(p * p+
q *
q);
610 for(
int j = n-1;
j <
nn; ++
j)
613 H(n-1,
j) =
q * z + p *
H(n,
j);
614 H(n,
j) =
q *
H(n,
j) - p * z;
619 for(
int i = 0;
i <= n; ++
i)
622 H(
i, n-1) =
q * z + p *
H(
i, n);
623 H(
i, n) =
q *
H(
i, n) - p * z;
631 V(
i, n-1) =
q * z + p * V(
i, n);
632 V(
i, n) =
q * V(
i, n) - p * z;
662 w =
H(n, n-1) *
H(n-1, n);
670 for(
int i =
low;
i <= n; ++
i)
674 s = abs(
H(n, n-1)) + abs(
H(n-1, n-2));
687 s = VIGRA_CSTD::sqrt(s);
692 s = x - w / ((y - x) / 2.0 + s);
693 for(
int i =
low;
i <= n; ++
i)
707 int m = hessenbergQrDecompositionHelper(
H, n,
l,
eps, p,
q, r, s, w, x, y, z);
708 for(
int i =
m+2;
i <= n; ++
i)
719 for(
int k =
m;
k <= n-1; ++
k)
726 x = abs(p) + abs(
q) + abs(r);
738 s = VIGRA_CSTD::sqrt(p * p +
q *
q + r * r);
762 for(
int j =
k;
j <
nn; ++
j)
767 p = p + r *
H(
k+2,
j);
768 H(
k+2,
j) =
H(
k+2,
j) - p * z;
770 H(
k,
j) =
H(
k,
j) - p * x;
771 H(
k+1,
j) =
H(
k+1,
j) - p * y;
776 for(
int i = 0;
i <= std::min(n,
k+3); ++
i)
778 p = x *
H(
i,
k) + y *
H(
i,
k+1);
781 p = p + z *
H(
i,
k+2);
782 H(
i,
k+2) =
H(
i,
k+2) - p * r;
792 p = x * V(
i,
k) + y * V(
i,
k+1);
795 p = p + z * V(
i,
k+2);
796 V(
i,
k+2) = V(
i,
k+2) - p * r;
798 V(
i,
k) = V(
i,
k) - p;
799 V(
i,
k+1) = V(
i,
k+1) - p *
q;
813 for(n =
nn-1; n >= 0; --n)
824 for(
int i = n-1;
i >= 0; --
i)
828 for(
int j =
l;
j <= n; ++
j)
830 r = r +
H(
i,
j) *
H(
j, n);
858 q = (d(
i) - p) * (d(
i) - p) + e(
i) * e(
i);
859 t = (x * s - z * r) /
q;
863 H(
i+1, n) = (-r - w * t) / x;
867 H(
i+1, n) = (-s - y * t) / z;
874 if((
eps * t) * t > 1)
876 for(
int j =
i;
j <= n; ++
j)
878 H(
j, n) =
H(
j, n) / t;
893 if(abs(
H(n, n-1)) > abs(
H(n-1, n)))
895 H(n-1, n-1) =
q /
H(n, n-1);
896 H(n-1, n) = -(
H(n, n) - p) /
H(n, n-1);
900 cdiv(0.0,-
H(n-1, n),
H(n-1, n-1)-p,
q,
H(n-1, n-1),
H(n-1, n));
904 for(
int i = n-2;
i >= 0; --
i)
909 for(
int j =
l;
j <= n; ++
j)
935 vr = (d(
i) - p) * (d(
i) - p) + e(
i) * e(
i) -
q *
q;
936 vi = (d(
i) - p) * 2.0 *
q;
937 if((
vr == 0.0) && (
vi == 0.0))
940 abs(x) + abs(y) + abs(z));
943 if(abs(x) > (abs(z) + abs(
q)))
945 H(
i+1, n-1) = (-
ra - w *
H(
i, n-1) +
q *
H(
i, n)) / x;
946 H(
i+1, n) = (-
sa - w *
H(
i, n) -
q *
H(
i, n-1)) / x;
950 cdiv(-r-y*
H(
i, n-1),-s-y*
H(
i, n),z,
q,
H(
i+1, n-1),
H(
i+1, n));
956 t = std::max(abs(
H(
i, n-1)),abs(
H(
i, n)));
957 if((
eps * t) * t > 1)
959 for(
int j =
i;
j <= n; ++
j)
961 H(
j, n-1) =
H(
j, n-1) / t;
962 H(
j, n) =
H(
j, n) / t;
979 z = z + V(
i,
k) *
H(
k,
j);
1006template <
class T,
class C1,
class C2,
class C3>
1012 "symmetricEigensystem(): symmetric input matrix required.");
1016 "symmetricEigensystem(): matrix shape mismatch.");
1020 detail::housholderTridiagonalization(
ev,
de);
1021 if(!detail::tridiagonalMatrixEigensystem(
de,
ev))
1030template <
class T,
class C2,
class C3>
1032symmetricEigensystem2x2(T
a00, T
a01, T
a11,
1035 double evec[2]={0,0};
1039 if (fabs(a11)>fabs(a00)){
1045 else if(fabs(a00)>fabs(a11)) {
1052 evec[0]=.5* M_SQRT2;
1053 evec[1]=.5* M_SQRT2;
1059 double temp=a11-a00;
1061 double coherence=sqrt(temp*temp+4*a01*a01);
1063 evec[1]=temp+coherence;
1064 temp=std::sqrt(evec[0]*evec[0]+evec[1]*evec[1]);
1066 evec[0]=.5* M_SQRT2;
1067 evec[1]=.5* M_SQRT2;
1076 ew(0,0)=.5*(a00+a11+coherence);
1077 ew(1,0)=.5*(a00+a11-coherence);
1087template <
class T,
class C2,
class C3>
1089symmetricEigensystem3x3(T a00, T a01, T a02, T a11, T a12, T a22,
1090 MultiArrayView<2, T, C2> & ew, MultiArrayView<2, T, C3> & ev)
1093 &ew(0,0), &ew(1,0), &ew(2,0));
1100 double b1=a00-ew(0,0),
1124 double M = std::max(std::max(abs(ew(0,0)), abs(ew(1,0))), abs(ew(2,0)));
1125 double epsilon = 1e-12*M;
1126 if(l0<epsilon) {
return false; }
1127 if(l1<epsilon) {
return false; }
1128 if(l2<epsilon) {
return false; }
1150template <
class T,
class C1,
class C2,
class C3>
1156 "symmetricEigensystemNoniterative(): symmetric input matrix required.");
1160 detail::symmetricEigensystem2x2(a(0,0), a(0,1), a(1,1),
ew,
ev);
1166 if(detail::symmetricEigensystem3x3(a(0,0), a(0,1), a(0,2), a(1,1), a(1,2), a(2,2),
1172 vigra_precondition(
false,
1173 "symmetricEigensystemNoniterative(): can only handle 2x2 and 3x3 matrices.");
1192template <
class T,
class C1,
class C2,
class C3>
1199 "nonsymmetricEigensystem(): square input matrix required.");
1202 "nonsymmetricEigensystem(): matrix shape mismatch.");
1206 detail::nonsymmetricHessenbergReduction(
H,
ev);
1207 if(!detail::hessenbergQrDecomposition(
H,
ev,
de))
1212 ew(
i,0) = std::complex<T>(
de(
i, 0),
de(
i, 1));
1232template <
class POLYNOMIAL,
class VECTOR>
1235 typedef typename POLYNOMIAL::value_type T;
1236 typedef typename POLYNOMIAL::Real Real;
1237 typedef typename POLYNOMIAL::Complex Complex;
1241 int const degree =
poly.order();
1242 double const eps =
poly.epsilon();
1245 for(
int i = 0;
i < degree; ++
i)
1247 for(
int i = 0;
i < degree - 1; ++
i)
1248 inMatrix(
i + 1,
i) = NumericTraits<T>::one();
1254 for(
int i = 0;
i < degree; ++
i)
1257 vigra::detail::laguerre1Root(
poly,
ew(
i,0), 1);
1258 roots.push_back(vigra::detail::deleteBelowEpsilon(
ew(
i,0),
eps));
1264template <
class POLYNOMIAL,
class VECTOR>
1267 return polynomialRootsEigenvalueMethod(
poly,
roots,
true);
1287template <
class POLYNOMIAL,
class VECTOR>
1290 typedef typename NumericTraits<typename VECTOR::value_type>::ComplexPromote Complex;
1300template <
class POLYNOMIAL,
class VECTOR>
1303 return polynomialRealRootsEigenvalueMethod(p,
roots,
true);
Base class for, and view to, vigra::MultiArray.
Definition multi_array.hxx:705
Class for a single RGB value.
Definition rgbvalue.hxx:128
size_type size() const
Definition tinyvector.hxx:913
iterator end()
Definition tinyvector.hxx:864
iterator begin()
Definition tinyvector.hxx:861
TinyVector & copy(TinyVectorBase< U, USIZE, DATA, DERIVED > const &r)
Definition tinyvector.hxx:1210
Matrix< T, ALLLOC >::NormType norm(const Matrix< T, ALLLOC > &a)
bool symmetricEigensystem(MultiArrayView< 2, T, C1 > const &a, MultiArrayView< 2, T, C2 > &ew, MultiArrayView< 2, T, C3 > &ev)
Definition eigensystem.hxx:1008
MultiArrayIndex columnCount(const MultiArrayView< 2, T, C > &x)
Definition matrix.hxx:684
bool nonsymmetricEigensystem(MultiArrayView< 2, T, C1 > const &a, MultiArrayView< 2, std::complex< T >, C2 > &ew, MultiArrayView< 2, T, C3 > &ev)
Definition eigensystem.hxx:1194
bool polynomialRootsEigenvalueMethod(POLYNOMIAL const &poly, VECTOR &roots, bool polishRoots)
Definition eigensystem.hxx:1233
MultiArrayIndex rowCount(const MultiArrayView< 2, T, C > &x)
Definition matrix.hxx:671
MultiArrayView< 2, T, C > columnVector(MultiArrayView< 2, T, C > const &m, MultiArrayIndex d)
Definition matrix.hxx:727
bool symmetricEigensystemNoniterative(MultiArrayView< 2, T, C1 > const &a, MultiArrayView< 2, T, C2 > &ew, MultiArrayView< 2, T, C3 > &ev)
Definition eigensystem.hxx:1152
bool polynomialRealRootsEigenvalueMethod(POLYNOMIAL const &p, VECTOR &roots, bool)
Definition eigensystem.hxx:1288
bool isSymmetric(const MultiArrayView< 2, T, C > &v)
Definition matrix.hxx:779
R imag(const FFTWComplex< R > &a)
imaginary part
Definition fftw3.hxx:1023
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition fftw3.hxx:1037
R real(const FFTWComplex< R > &a)
real part
Definition fftw3.hxx:1016
void symmetric3x3Eigenvalues(T a00, T a01, T a02, T a11, T a12, T a22, T *r0, T *r1, T *r2)
Compute the eigenvalues of a 3x3 real symmetric matrix.
Definition mathutil.hxx:754
std::ptrdiff_t MultiArrayIndex
Definition multi_fwd.hxx:60
FixedPoint16< IntBits, OverflowHandling > hypot(FixedPoint16< IntBits, OverflowHandling > v1, FixedPoint16< IntBits, OverflowHandling > v2)
Length of hypotenuse.
Definition fixedpoint.hxx:1640