decomposition.hpp
Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 2012. Philipp Wagner <bytefish[at]gmx[dot]de>.
00003  * Released to public domain under terms of the BSD Simplified license.
00004  *
00005  * Redistribution and use in source and binary forms, with or without
00006  * modification, are permitted provided that the following conditions are met:
00007  *   * Redistributions of source code must retain the above copyright
00008  *     notice, this list of conditions and the following disclaimer.
00009  *   * Redistributions in binary form must reproduce the above copyright
00010  *     notice, this list of conditions and the following disclaimer in the
00011  *     documentation and/or other materials provided with the distribution.
00012  *   * Neither the name of the organization nor the names of its contributors
00013  *     may be used to endorse or promote products derived from this software
00014  *     without specific prior written permission.
00015  *
00016  *   See <http://www.opensource.org/licenses/bsd-license>
00017  */
00018 #ifndef __DECOMPOSITION_HPP__
00019 #define __DECOMPOSITION_HPP__
00020 
00021 #include "opencv2/core/core.hpp"
00022 
00023 using namespace cv;
00024 using namespace std;
00043 namespace cv
00044 {
00045 
00046 class EigenvalueDecomposition
00047 {
00048 private:
00049 
00050         // Holds the data dimension.
00051         int n;
00052 
00053         // Stores real/imag part of a complex division.
00054         double cdivr, cdivi;
00055 
00056         // Pointer to internal memory.
00057         double *d, *e, *ort;
00058         double **V, **H;
00059 
00060         // Holds the computed eigenvalues.
00061         Mat _eigenvalues;
00062 
00063         // Holds the computed eigenvectors.
00064         Mat _eigenvectors;
00065 
00066         // Allocates memory.
00067         template<typename _Tp>
00068         _Tp *alloc_1d(int m)
00069         {
00070                 return new _Tp[m];
00071         }
00072 
00073         // Allocates memory.
00074         template<typename _Tp>
00075         _Tp *alloc_1d(int m, _Tp val)
00076         {
00077                 _Tp *arr = alloc_1d<_Tp>(m);
00078                 for (int i = 0; i < m; i++)
00079                         arr[i] = val;
00080                 return arr;
00081         }
00082 
00083         // Allocates memory.
00084         template<typename _Tp>
00085         _Tp **alloc_2d(int m, int n)
00086         {
00087                 _Tp **arr = new _Tp*[m];
00088                 for (int i = 0; i < m; i++)
00089                         arr[i] = new _Tp[n];
00090                 return arr;
00091         }
00092 
00093         // Allocates memory.
00094         template<typename _Tp>
00095         _Tp **alloc_2d(int m, int n, _Tp val)
00096         {
00097                 _Tp **arr = alloc_2d<_Tp>(m, n);
00098                 for (int i = 0; i < m; i++)
00099                 {
00100                         for (int j = 0; j < n; j++)
00101                         {
00102                                 arr[i][j] = val;
00103                         }
00104                 }
00105                 return arr;
00106         }
00107 
00108         void cdiv(double xr, double xi, double yr, double yi)
00109         {
00110                 double r, d;
00111                 if (std::abs(yr) > std::abs(yi))
00112                 {
00113                         r = yi / yr;
00114                         d = yr + r * yi;
00115                         cdivr = (xr + r * xi) / d;
00116                         cdivi = (xi - r * xr) / d;
00117                 }
00118                 else
00119                 {
00120                         r = yr / yi;
00121                         d = yi + r * yr;
00122                         cdivr = (r * xr + xi) / d;
00123                         cdivi = (r * xi - xr) / d;
00124                 }
00125         }
00126 
00127         // Nonsymmetric reduction from Hessenberg to real Schur form.
00128 
00129         void hqr2()
00130         {
00131 
00132                 //  This is derived from the Algol procedure hqr2,
00133                 //  by Martin and Wilkinson, Handbook for Auto. Comp.,
00134                 //  Vol.ii-Linear Algebra, and the corresponding
00135                 //  Fortran subroutine in EISPACK.
00136 
00137                 // Initialize
00138                 int nn = this->n;
00139                 int n = nn - 1;
00140                 int low = 0;
00141                 int high = nn - 1;
00142                 double eps = pow(2.0, -52.0);
00143                 double exshift = 0.0;
00144                 double p = 0, q = 0, r = 0, s = 0, z = 0, t, w, x, y;
00145 
00146                 // Store roots isolated by balanc and compute matrix norm
00147 
00148                 double norm = 0.0;
00149                 for (int i = 0; i < nn; i++)
00150                 {
00151                         if (i < low | i > high)
00152                         {
00153                                 d[i] = H[i][i];
00154                                 e[i] = 0.0;
00155                         }
00156                         for (int j = max(i - 1, 0); j < nn; j++)
00157                         {
00158                                 norm = norm + std::abs(H[i][j]);
00159                         }
00160                 }
00161 
00162                 // Outer loop over eigenvalue index
00163                 int iter = 0;
00164                 while (n >= low)
00165                 {
00166 
00167                         // Look for single small sub-diagonal element
00168                         int l = n;
00169                         while (l > low)
00170                         {
00171                                 s = std::abs(H[l - 1][l - 1]) + std::abs(H[l][l]);
00172                                 if (s == 0.0)
00173                                 {
00174                                         s = norm;
00175                                 }
00176                                 if (std::abs(H[l][l - 1]) < eps * s)
00177                                 {
00178                                         break;
00179                                 }
00180                                 l--;
00181                         }
00182 
00183                         // Check for convergence
00184                         // One root found
00185 
00186                         if (l == n)
00187                         {
00188                                 H[n][n] = H[n][n] + exshift;
00189                                 d[n] = H[n][n];
00190                                 e[n] = 0.0;
00191                                 n--;
00192                                 iter = 0;
00193 
00194                                 // Two roots found
00195 
00196                         }
00197                         else if (l == n - 1)
00198                         {
00199                                 w = H[n][n - 1] * H[n - 1][n];
00200                                 p = (H[n - 1][n - 1] - H[n][n]) / 2.0;
00201                                 q = p * p + w;
00202                                 z = sqrt(std::abs(q));
00203                                 H[n][n] = H[n][n] + exshift;
00204                                 H[n - 1][n - 1] = H[n - 1][n - 1] + exshift;
00205                                 x = H[n][n];
00206 
00207                                 // Real pair
00208 
00209                                 if (q >= 0)
00210                                 {
00211                                         if (p >= 0)
00212                                         {
00213                                                 z = p + z;
00214                                         }
00215                                         else
00216                                         {
00217                                                 z = p - z;
00218                                         }
00219                                         d[n - 1] = x + z;
00220                                         d[n] = d[n - 1];
00221                                         if (z != 0.0)
00222                                         {
00223                                                 d[n] = x - w / z;
00224                                         }
00225                                         e[n - 1] = 0.0;
00226                                         e[n] = 0.0;
00227                                         x = H[n][n - 1];
00228                                         s = std::abs(x) + std::abs(z);
00229                                         p = x / s;
00230                                         q = z / s;
00231                                         r = sqrt(p * p + q * q);
00232                                         p = p / r;
00233                                         q = q / r;
00234 
00235                                         // Row modification
00236 
00237                                         for (int j = n - 1; j < nn; j++)
00238                                         {
00239                                                 z = H[n - 1][j];
00240                                                 H[n - 1][j] = q * z + p * H[n][j];
00241                                                 H[n][j] = q * H[n][j] - p * z;
00242                                         }
00243 
00244                                         // Column modification
00245 
00246                                         for (int i = 0; i <= n; i++)
00247                                         {
00248                                                 z = H[i][n - 1];
00249                                                 H[i][n - 1] = q * z + p * H[i][n];
00250                                                 H[i][n] = q * H[i][n] - p * z;
00251                                         }
00252 
00253                                         // Accumulate transformations
00254 
00255                                         for (int i = low; i <= high; i++)
00256                                         {
00257                                                 z = V[i][n - 1];
00258                                                 V[i][n - 1] = q * z + p * V[i][n];
00259                                                 V[i][n] = q * V[i][n] - p * z;
00260                                         }
00261 
00262                                         // Complex pair
00263 
00264                                 }
00265                                 else
00266                                 {
00267                                         d[n - 1] = x + p;
00268                                         d[n] = x + p;
00269                                         e[n - 1] = z;
00270                                         e[n] = -z;
00271                                 }
00272                                 n = n - 2;
00273                                 iter = 0;
00274 
00275                                 // No convergence yet
00276 
00277                         }
00278                         else
00279                         {
00280 
00281                                 // Form shift
00282 
00283                                 x = H[n][n];
00284                                 y = 0.0;
00285                                 w = 0.0;
00286                                 if (l < n)
00287                                 {
00288                                         y = H[n - 1][n - 1];
00289                                         w = H[n][n - 1] * H[n - 1][n];
00290                                 }
00291 
00292                                 // Wilkinson's original ad hoc shift
00293 
00294                                 if (iter == 10)
00295                                 {
00296                                         exshift += x;
00297                                         for (int i = low; i <= n; i++)
00298                                         {
00299                                                 H[i][i] -= x;
00300                                         }
00301                                         s = std::abs(H[n][n - 1]) + std::abs(H[n - 1][n - 2]);
00302                                         x = y = 0.75 * s;
00303                                         w = -0.4375 * s * s;
00304                                 }
00305 
00306                                 // MATLAB's new ad hoc shift
00307 
00308                                 if (iter == 30)
00309                                 {
00310                                         s = (y - x) / 2.0;
00311                                         s = s * s + w;
00312                                         if (s > 0)
00313                                         {
00314                                                 s = sqrt(s);
00315                                                 if (y < x)
00316                                                 {
00317                                                         s = -s;
00318                                                 }
00319                                                 s = x - w / ((y - x) / 2.0 + s);
00320                                                 for (int i = low; i <= n; i++)
00321                                                 {
00322                                                         H[i][i] -= s;
00323                                                 }
00324                                                 exshift += s;
00325                                                 x = y = w = 0.964;
00326                                         }
00327                                 }
00328 
00329                                 iter = iter + 1; // (Could check iteration count here.)
00330 
00331                                 // Look for two consecutive small sub-diagonal elements
00332                                 int m = n - 2;
00333                                 while (m >= l)
00334                                 {
00335                                         z = H[m][m];
00336                                         r = x - z;
00337                                         s = y - z;
00338                                         p = (r * s - w) / H[m + 1][m] + H[m][m + 1];
00339                                         q = H[m + 1][m + 1] - z - r - s;
00340                                         r = H[m + 2][m + 1];
00341                                         s = std::abs(p) + std::abs(q) + std::abs(r);
00342                                         p = p / s;
00343                                         q = q / s;
00344                                         r = r / s;
00345                                         if (m == l)
00346                                         {
00347                                                 break;
00348                                         }
00349                                         if (std::abs(H[m][m - 1]) * (std::abs(q) + std::abs(r)) < eps * (std::abs(p) * (std::abs(H[m - 1][m - 1]) + std::abs(z) + std::abs(H[m + 1][m + 1]))))
00350                                         {
00351                                                 break;
00352                                         }
00353                                         m--;
00354                                 }
00355 
00356                                 for (int i = m + 2; i <= n; i++)
00357                                 {
00358                                         H[i][i - 2] = 0.0;
00359                                         if (i > m + 2)
00360                                         {
00361                                                 H[i][i - 3] = 0.0;
00362                                         }
00363                                 }
00364 
00365                                 // Double QR step involving rows l:n and columns m:n
00366 
00367                                 for (int k = m; k <= n - 1; k++)
00368                                 {
00369                                         bool notlast = (k != n - 1);
00370                                         if (k != m)
00371                                         {
00372                                                 p = H[k][k - 1];
00373                                                 q = H[k + 1][k - 1];
00374                                                 r = (notlast ? H[k + 2][k - 1] : 0.0);
00375                                                 x = std::abs(p) + std::abs(q) + std::abs(r);
00376                                                 if (x != 0.0)
00377                                                 {
00378                                                         p = p / x;
00379                                                         q = q / x;
00380                                                         r = r / x;
00381                                                 }
00382                                         }
00383                                         if (x == 0.0)
00384                                         {
00385                                                 break;
00386                                         }
00387                                         s = sqrt(p * p + q * q + r * r);
00388                                         if (p < 0)
00389                                         {
00390                                                 s = -s;
00391                                         }
00392                                         if (s != 0)
00393                                         {
00394                                                 if (k != m)
00395                                                 {
00396                                                         H[k][k - 1] = -s * x;
00397                                                 }
00398                                                 else if (l != m)
00399                                                 {
00400                                                         H[k][k - 1] = -H[k][k - 1];
00401                                                 }
00402                                                 p = p + s;
00403                                                 x = p / s;
00404                                                 y = q / s;
00405                                                 z = r / s;
00406                                                 q = q / p;
00407                                                 r = r / p;
00408 
00409                                                 // Row modification
00410 
00411                                                 for (int j = k; j < nn; j++)
00412                                                 {
00413                                                         p = H[k][j] + q * H[k + 1][j];
00414                                                         if (notlast)
00415                                                         {
00416                                                                 p = p + r * H[k + 2][j];
00417                                                                 H[k + 2][j] = H[k + 2][j] - p * z;
00418                                                         }
00419                                                         H[k][j] = H[k][j] - p * x;
00420                                                         H[k + 1][j] = H[k + 1][j] - p * y;
00421                                                 }
00422 
00423                                                 // Column modification
00424 
00425                                                 for (int i = 0; i <= min(n, k + 3); i++)
00426                                                 {
00427                                                         p = x * H[i][k] + y * H[i][k + 1];
00428                                                         if (notlast)
00429                                                         {
00430                                                                 p = p + z * H[i][k + 2];
00431                                                                 H[i][k + 2] = H[i][k + 2] - p * r;
00432                                                         }
00433                                                         H[i][k] = H[i][k] - p;
00434                                                         H[i][k + 1] = H[i][k + 1] - p * q;
00435                                                 }
00436 
00437                                                 // Accumulate transformations
00438 
00439                                                 for (int i = low; i <= high; i++)
00440                                                 {
00441                                                         p = x * V[i][k] + y * V[i][k + 1];
00442                                                         if (notlast)
00443                                                         {
00444                                                                 p = p + z * V[i][k + 2];
00445                                                                 V[i][k + 2] = V[i][k + 2] - p * r;
00446                                                         }
00447                                                         V[i][k] = V[i][k] - p;
00448                                                         V[i][k + 1] = V[i][k + 1] - p * q;
00449                                                 }
00450                                         } // (s != 0)
00451                                 } // k loop
00452                         } // check convergence
00453                 } // while (n >= low)
00454 
00455                 // Backsubstitute to find vectors of upper triangular form
00456 
00457                 if (norm == 0.0)
00458                 {
00459                         return;
00460                 }
00461 
00462                 for (n = nn - 1; n >= 0; n--)
00463                 {
00464                         p = d[n];
00465                         q = e[n];
00466 
00467                         // Real vector
00468 
00469                         if (q == 0)
00470                         {
00471                                 int l = n;
00472                                 H[n][n] = 1.0;
00473                                 for (int i = n - 1; i >= 0; i--)
00474                                 {
00475                                         w = H[i][i] - p;
00476                                         r = 0.0;
00477                                         for (int j = l; j <= n; j++)
00478                                         {
00479                                                 r = r + H[i][j] * H[j][n];
00480                                         }
00481                                         if (e[i] < 0.0)
00482                                         {
00483                                                 z = w;
00484                                                 s = r;
00485                                         }
00486                                         else
00487                                         {
00488                                                 l = i;
00489                                                 if (e[i] == 0.0)
00490                                                 {
00491                                                         if (w != 0.0)
00492                                                         {
00493                                                                 H[i][n] = -r / w;
00494                                                         }
00495                                                         else
00496                                                         {
00497                                                                 H[i][n] = -r / (eps * norm);
00498                                                         }
00499 
00500                                                         // Solve real equations
00501 
00502                                                 }
00503                                                 else
00504                                                 {
00505                                                         x = H[i][i + 1];
00506                                                         y = H[i + 1][i];
00507                                                         q = (d[i] - p) * (d[i] - p) + e[i] * e[i];
00508                                                         t = (x * s - z * r) / q;
00509                                                         H[i][n] = t;
00510                                                         if (std::abs(x) > std::abs(z))
00511                                                         {
00512                                                                 H[i + 1][n] = (-r - w * t) / x;
00513                                                         }
00514                                                         else
00515                                                         {
00516                                                                 H[i + 1][n] = (-s - y * t) / z;
00517                                                         }
00518                                                 }
00519 
00520                                                 // Overflow control
00521 
00522                                                 t = std::abs(H[i][n]);
00523                                                 if ((eps * t) * t > 1)
00524                                                 {
00525                                                         for (int j = i; j <= n; j++)
00526                                                         {
00527                                                                 H[j][n] = H[j][n] / t;
00528                                                         }
00529                                                 }
00530                                         }
00531                                 }
00532 
00533                                 // Complex vector
00534 
00535                         }
00536                         else if (q < 0)
00537                         {
00538                                 int l = n - 1;
00539 
00540                                 // Last vector component imaginary so matrix is triangular
00541 
00542                                 if (std::abs(H[n][n - 1]) > std::abs(H[n - 1][n]))
00543                                 {
00544                                         H[n - 1][n - 1] = q / H[n][n - 1];
00545                                         H[n - 1][n] = -(H[n][n] - p) / H[n][n - 1];
00546                                 }
00547                                 else
00548                                 {
00549                                         cdiv(0.0, -H[n - 1][n], H[n - 1][n - 1] - p, q);
00550                                         H[n - 1][n - 1] = cdivr;
00551                                         H[n - 1][n] = cdivi;
00552                                 }
00553                                 H[n][n - 1] = 0.0;
00554                                 H[n][n] = 1.0;
00555                                 for (int i = n - 2; i >= 0; i--)
00556                                 {
00557                                         double ra, sa, vr, vi;
00558                                         ra = 0.0;
00559                                         sa = 0.0;
00560                                         for (int j = l; j <= n; j++)
00561                                         {
00562                                                 ra = ra + H[i][j] * H[j][n - 1];
00563                                                 sa = sa + H[i][j] * H[j][n];
00564                                         }
00565                                         w = H[i][i] - p;
00566 
00567                                         if (e[i] < 0.0)
00568                                         {
00569                                                 z = w;
00570                                                 r = ra;
00571                                                 s = sa;
00572                                         }
00573                                         else
00574                                         {
00575                                                 l = i;
00576                                                 if (e[i] == 0)
00577                                                 {
00578                                                         cdiv(-ra, -sa, w, q);
00579                                                         H[i][n - 1] = cdivr;
00580                                                         H[i][n] = cdivi;
00581                                                 }
00582                                                 else
00583                                                 {
00584 
00585                                                         // Solve complex equations
00586 
00587                                                         x = H[i][i + 1];
00588                                                         y = H[i + 1][i];
00589                                                         vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
00590                                                         vi = (d[i] - p) * 2.0 * q;
00591                                                         if (vr == 0.0 & vi == 0.0)
00592                                                         {
00593                                                                 vr = eps * norm * (std::abs(w) + std::abs(q) + std::abs(x) + std::abs(y) + std::abs(z));
00594                                                         }
00595                                                         cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi);
00596                                                         H[i][n - 1] = cdivr;
00597                                                         H[i][n] = cdivi;
00598                                                         if (std::abs(x) > (std::abs(z) + std::abs(q)))
00599                                                         {
00600                                                                 H[i + 1][n - 1] = (-ra - w * H[i][n - 1] + q * H[i][n]) / x;
00601                                                                 H[i + 1][n] = (-sa - w * H[i][n] - q * H[i][n - 1]) / x;
00602                                                         }
00603                                                         else
00604                                                         {
00605                                                                 cdiv(-r - y * H[i][n - 1], -s - y * H[i][n], z, q);
00606                                                                 H[i + 1][n - 1] = cdivr;
00607                                                                 H[i + 1][n] = cdivi;
00608                                                         }
00609                                                 }
00610 
00611                                                 // Overflow control
00612 
00613                                                 t = max(std::abs(H[i][n - 1]), std::abs(H[i][n]));
00614                                                 if ((eps * t) * t > 1)
00615                                                 {
00616                                                         for (int j = i; j <= n; j++)
00617                                                         {
00618                                                                 H[j][n - 1] = H[j][n - 1] / t;
00619                                                                 H[j][n] = H[j][n] / t;
00620                                                         }
00621                                                 }
00622                                         }
00623                                 }
00624                         }
00625                 }
00626 
00627                 // Vectors of isolated roots
00628 
00629                 for (int i = 0; i < nn; i++)
00630                 {
00631                         if (i < low | i > high)
00632                         {
00633                                 for (int j = i; j < nn; j++)
00634                                 {
00635                                         V[i][j] = H[i][j];
00636                                 }
00637                         }
00638                 }
00639 
00640                 // Back transformation to get eigenvectors of original matrix
00641 
00642                 for (int j = nn - 1; j >= low; j--)
00643                 {
00644                         for (int i = low; i <= high; i++)
00645                         {
00646                                 z = 0.0;
00647                                 for (int k = low; k <= min(j, high); k++)
00648                                 {
00649                                         z = z + V[i][k] * H[k][j];
00650                                 }
00651                                 V[i][j] = z;
00652                         }
00653                 }
00654         }
00655 
00656         // Nonsymmetric reduction to Hessenberg form.
00657         void orthes()
00658         {
00659                 //  This is derived from the Algol procedures orthes and ortran,
00660                 //  by Martin and Wilkinson, Handbook for Auto. Comp.,
00661                 //  Vol.ii-Linear Algebra, and the corresponding
00662                 //  Fortran subroutines in EISPACK.
00663                 int low = 0;
00664                 int high = n - 1;
00665 
00666                 for (int m = low + 1; m <= high - 1; m++)
00667                 {
00668 
00669                         // Scale column.
00670 
00671                         double scale = 0.0;
00672                         for (int i = m; i <= high; i++)
00673                         {
00674                                 scale = scale + std::abs(H[i][m - 1]);
00675                         }
00676                         if (scale != 0.0)
00677                         {
00678 
00679                                 // Compute Householder transformation.
00680 
00681                                 double h = 0.0;
00682                                 for (int i = high; i >= m; i--)
00683                                 {
00684                                         ort[i] = H[i][m - 1] / scale;
00685                                         h += ort[i] * ort[i];
00686                                 }
00687                                 double g = sqrt(h);
00688                                 if (ort[m] > 0)
00689                                 {
00690                                         g = -g;
00691                                 }
00692                                 h = h - ort[m] * g;
00693                                 ort[m] = ort[m] - g;
00694 
00695                                 // Apply Householder similarity transformation
00696                                 // H = (I-u*u'/h)*H*(I-u*u')/h)
00697 
00698                                 for (int j = m; j < n; j++)
00699                                 {
00700                                         double f = 0.0;
00701                                         for (int i = high; i >= m; i--)
00702                                         {
00703                                                 f += ort[i] * H[i][j];
00704                                         }
00705                                         f = f / h;
00706                                         for (int i = m; i <= high; i++)
00707                                         {
00708                                                 H[i][j] -= f * ort[i];
00709                                         }
00710                                 }
00711 
00712                                 for (int i = 0; i <= high; i++)
00713                                 {
00714                                         double f = 0.0;
00715                                         for (int j = high; j >= m; j--)
00716                                         {
00717                                                 f += ort[j] * H[i][j];
00718                                         }
00719                                         f = f / h;
00720                                         for (int j = m; j <= high; j++)
00721                                         {
00722                                                 H[i][j] -= f * ort[j];
00723                                         }
00724                                 }
00725                                 ort[m] = scale * ort[m];
00726                                 H[m][m - 1] = scale * g;
00727                         }
00728                 }
00729 
00730                 // Accumulate transformations (Algol's ortran).
00731 
00732                 for (int i = 0; i < n; i++)
00733                 {
00734                         for (int j = 0; j < n; j++)
00735                         {
00736                                 V[i][j] = (i == j ? 1.0 : 0.0);
00737                         }
00738                 }
00739 
00740                 for (int m = high - 1; m >= low + 1; m--)
00741                 {
00742                         if (H[m][m - 1] != 0.0)
00743                         {
00744                                 for (int i = m + 1; i <= high; i++)
00745                                 {
00746                                         ort[i] = H[i][m - 1];
00747                                 }
00748                                 for (int j = m; j <= high; j++)
00749                                 {
00750                                         double g = 0.0;
00751                                         for (int i = m; i <= high; i++)
00752                                         {
00753                                                 g += ort[i] * V[i][j];
00754                                         }
00755                                         // Double division avoids possible underflow
00756                                         g = (g / ort[m]) / H[m][m - 1];
00757                                         for (int i = m; i <= high; i++)
00758                                         {
00759                                                 V[i][j] += g * ort[i];
00760                                         }
00761                                 }
00762                         }
00763                 }
00764         }
00765 
00766         // Releases all internal working memory.
00767         void release()
00768         {
00769                 // releases the working data
00770                 delete[] d, e, ort;
00771                 for (int i = 0; i < n; i++)
00772                 {
00773                         delete[] H[i];
00774                         delete[] V[i];
00775                 }
00776                 delete[] H, V;
00777         }
00778 
00779         // Computes the Eigenvalue Decomposition for a matrix given in H.
00780         void compute()
00781         {
00782                 // Allocate memory for the working data.
00783                 V = alloc_2d<double>(n, n, 0.0);
00784                 d = alloc_1d<double>(n);
00785                 e = alloc_1d<double>(n);
00786                 ort = alloc_1d<double>(n);
00787                 // Reduce to Hessenberg form.
00788                 orthes();
00789                 // Reduce Hessenberg to real Schur form.
00790                 hqr2();
00791                 // Copy eigenvalues to OpenCV Matrix.
00792                 _eigenvalues.create(1, n, CV_64FC1);
00793                 for (int i = 0; i < n; i++)
00794                 {
00795                         _eigenvalues.at<double>(0, i) = d[i];
00796                 }
00797                 // Copy eigenvectors to OpenCV Matrix.
00798                 _eigenvectors.create(n, n, CV_64FC1);
00799                 for (int i = 0; i < n; i++)
00800                         for (int j = 0; j < n; j++)
00801                                 _eigenvectors.at<double>(i, j) = V[i][j];
00802                 // Deallocate the memory by releasing all internal working data.
00803                 release();
00804         }
00805 
00806 public:
00807         EigenvalueDecomposition() :
00808                 n(0)
00809         {
00810         }
00811 
00812         // Initializes & computes the Eigenvalue Decomposition for a general matrix
00813         // given in src. This function is a port of the EigenvalueSolver in JAMA,
00814         // which has been released to public domain by The MathWorks and the
00815         // National Institute of Standards and Technology (NIST).
00816         EigenvalueDecomposition(InputArray src)
00817         {
00818                 compute(src);
00819         }
00820 
00821         // This function computes the Eigenvalue Decomposition for a general matrix
00822         // given in src. This function is a port of the EigenvalueSolver in JAMA,
00823         // which has been released to public domain by The MathWorks and the
00824         // National Institute of Standards and Technology (NIST).
00825         void compute(InputArray src);
00826 
00827         ~EigenvalueDecomposition()
00828         {
00829         }
00830 
00831         // Returns the eigenvalues of the Eigenvalue Decomposition.
00832         Mat eigenvalues()
00833         {
00834                 return _eigenvalues;
00835         }
00836         // Returns the eigenvectors of the Eigenvalue Decomposition.
00837         Mat eigenvectors()
00838         {
00839                 return _eigenvectors;
00840         }
00841 };
00842 
00843 } // namespace
00844 #endif


cob_people_detection
Author(s): Richard Bormann , Thomas Zwölfer
autogenerated on Fri Aug 28 2015 10:24:12