1 /* 2 * Copyright (c) 2011. Philipp Wagner <bytefish[at]gmx[dot]de>. 3 * Released to public domain under terms of the BSD Simplified license. 4 * 5 * Redistribution and use in source and binary forms, with or without 6 * modification, are permitted provided that the following conditions are met: 7 * * Redistributions of source code must retain the above copyright 8 * notice, this list of conditions and the following disclaimer. 9 * * Redistributions in binary form must reproduce the above copyright 10 * notice, this list of conditions and the following disclaimer in the 11 * documentation and/or other materials provided with the distribution. 12 * * Neither the name of the organization nor the names of its contributors 13 * may be used to endorse or promote products derived from this software 14 * without specific prior written permission. 15 * 16 * See <http://www.opensource.org/licenses/bsd-license> 17 */ 18 19 #include "precomp.hpp" 20 #include <iostream> 21 #include <map> 22 #include <set> 23 24 namespace cv 25 { 26 27 // Removes duplicate elements in a given vector. 28 template<typename _Tp> 29 inline std::vector<_Tp> remove_dups(const std::vector<_Tp>& src) { 30 typedef typename std::set<_Tp>::const_iterator constSetIterator; 31 typedef typename std::vector<_Tp>::const_iterator constVecIterator; 32 std::set<_Tp> set_elems; 33 for (constVecIterator it = src.begin(); it != src.end(); ++it) 34 set_elems.insert(*it); 35 std::vector<_Tp> elems; 36 for (constSetIterator it = set_elems.begin(); it != set_elems.end(); ++it) 37 elems.push_back(*it); 38 return elems; 39 } 40 41 static Mat argsort(InputArray _src, bool ascending=true) 42 { 43 Mat src = _src.getMat(); 44 if (src.rows != 1 && src.cols != 1) { 45 String error_message = "Wrong shape of input matrix! Expected a matrix with one row or column."; 46 CV_Error(Error::StsBadArg, error_message); 47 } 48 int flags = SORT_EVERY_ROW | (ascending ? SORT_ASCENDING : SORT_DESCENDING); 49 Mat sorted_indices; 50 sortIdx(src.reshape(1,1),sorted_indices,flags); 51 return sorted_indices; 52 } 53 54 static Mat asRowMatrix(InputArrayOfArrays src, int rtype, double alpha=1, double beta=0) { 55 // make sure the input data is a vector of matrices or vector of vector 56 if(src.kind() != _InputArray::STD_VECTOR_MAT && src.kind() != _InputArray::STD_VECTOR_VECTOR) { 57 String error_message = "The data is expected as InputArray::STD_VECTOR_MAT (a std::vector<Mat>) or _InputArray::STD_VECTOR_VECTOR (a std::vector< std::vector<...> >)."; 58 CV_Error(Error::StsBadArg, error_message); 59 } 60 // number of samples 61 size_t n = src.total(); 62 // return empty matrix if no matrices given 63 if(n == 0) 64 return Mat(); 65 // dimensionality of (reshaped) samples 66 size_t d = src.getMat(0).total(); 67 // create data matrix 68 Mat data((int)n, (int)d, rtype); 69 // now copy data 70 for(int i = 0; i < (int)n; i++) { 71 // make sure data can be reshaped, throw exception if not! 72 if(src.getMat(i).total() != d) { 73 String error_message = format("Wrong number of elements in matrix #%d! Expected %d was %d.", i, (int)d, (int)src.getMat(i).total()); 74 CV_Error(Error::StsBadArg, error_message); 75 } 76 // get a hold of the current row 77 Mat xi = data.row(i); 78 // make reshape happy by cloning for non-continuous matrices 79 if(src.getMat(i).isContinuous()) { 80 src.getMat(i).reshape(1, 1).convertTo(xi, rtype, alpha, beta); 81 } else { 82 src.getMat(i).clone().reshape(1, 1).convertTo(xi, rtype, alpha, beta); 83 } 84 } 85 return data; 86 } 87 88 static void sortMatrixColumnsByIndices(InputArray _src, InputArray _indices, OutputArray _dst) { 89 if(_indices.getMat().type() != CV_32SC1) { 90 CV_Error(Error::StsUnsupportedFormat, "cv::sortColumnsByIndices only works on integer indices!"); 91 } 92 Mat src = _src.getMat(); 93 std::vector<int> indices = _indices.getMat(); 94 _dst.create(src.rows, src.cols, src.type()); 95 Mat dst = _dst.getMat(); 96 for(size_t idx = 0; idx < indices.size(); idx++) { 97 Mat originalCol = src.col(indices[idx]); 98 Mat sortedCol = dst.col((int)idx); 99 originalCol.copyTo(sortedCol); 100 } 101 } 102 103 static Mat sortMatrixColumnsByIndices(InputArray src, InputArray indices) { 104 Mat dst; 105 sortMatrixColumnsByIndices(src, indices, dst); 106 return dst; 107 } 108 109 110 template<typename _Tp> static bool 111 isSymmetric_(InputArray src) { 112 Mat _src = src.getMat(); 113 if(_src.cols != _src.rows) 114 return false; 115 for (int i = 0; i < _src.rows; i++) { 116 for (int j = 0; j < _src.cols; j++) { 117 _Tp a = _src.at<_Tp> (i, j); 118 _Tp b = _src.at<_Tp> (j, i); 119 if (a != b) { 120 return false; 121 } 122 } 123 } 124 return true; 125 } 126 127 template<typename _Tp> static bool 128 isSymmetric_(InputArray src, double eps) { 129 Mat _src = src.getMat(); 130 if(_src.cols != _src.rows) 131 return false; 132 for (int i = 0; i < _src.rows; i++) { 133 for (int j = 0; j < _src.cols; j++) { 134 _Tp a = _src.at<_Tp> (i, j); 135 _Tp b = _src.at<_Tp> (j, i); 136 if (std::abs(a - b) > eps) { 137 return false; 138 } 139 } 140 } 141 return true; 142 } 143 144 static bool isSymmetric(InputArray src, double eps=1e-16) 145 { 146 Mat m = src.getMat(); 147 switch (m.type()) { 148 case CV_8SC1: return isSymmetric_<char>(m); break; 149 case CV_8UC1: 150 return isSymmetric_<unsigned char>(m); break; 151 case CV_16SC1: 152 return isSymmetric_<short>(m); break; 153 case CV_16UC1: 154 return isSymmetric_<unsigned short>(m); break; 155 case CV_32SC1: 156 return isSymmetric_<int>(m); break; 157 case CV_32FC1: 158 return isSymmetric_<float>(m, eps); break; 159 case CV_64FC1: 160 return isSymmetric_<double>(m, eps); break; 161 default: 162 break; 163 } 164 return false; 165 } 166 167 168 //------------------------------------------------------------------------------ 169 // cv::subspaceProject 170 //------------------------------------------------------------------------------ 171 Mat LDA::subspaceProject(InputArray _W, InputArray _mean, InputArray _src) { 172 // get data matrices 173 Mat W = _W.getMat(); 174 Mat mean = _mean.getMat(); 175 Mat src = _src.getMat(); 176 // get number of samples and dimension 177 int n = src.rows; 178 int d = src.cols; 179 // make sure the data has the correct shape 180 if(W.rows != d) { 181 String error_message = format("Wrong shapes for given matrices. Was size(src) = (%d,%d), size(W) = (%d,%d).", src.rows, src.cols, W.rows, W.cols); 182 CV_Error(Error::StsBadArg, error_message); 183 } 184 // make sure mean is correct if not empty 185 if(!mean.empty() && (mean.total() != (size_t) d)) { 186 String error_message = format("Wrong mean shape for the given data matrix. Expected %d, but was %d.", d, mean.total()); 187 CV_Error(Error::StsBadArg, error_message); 188 } 189 // create temporary matrices 190 Mat X, Y; 191 // make sure you operate on correct type 192 src.convertTo(X, W.type()); 193 // safe to do, because of above assertion 194 if(!mean.empty()) { 195 for(int i=0; i<n; i++) { 196 Mat r_i = X.row(i); 197 subtract(r_i, mean.reshape(1,1), r_i); 198 } 199 } 200 // finally calculate projection as Y = (X-mean)*W 201 gemm(X, W, 1.0, Mat(), 0.0, Y); 202 return Y; 203 } 204 205 //------------------------------------------------------------------------------ 206 // cv::subspaceReconstruct 207 //------------------------------------------------------------------------------ 208 Mat LDA::subspaceReconstruct(InputArray _W, InputArray _mean, InputArray _src) 209 { 210 // get data matrices 211 Mat W = _W.getMat(); 212 Mat mean = _mean.getMat(); 213 Mat src = _src.getMat(); 214 // get number of samples and dimension 215 int n = src.rows; 216 int d = src.cols; 217 // make sure the data has the correct shape 218 if(W.cols != d) { 219 String error_message = format("Wrong shapes for given matrices. Was size(src) = (%d,%d), size(W) = (%d,%d).", src.rows, src.cols, W.rows, W.cols); 220 CV_Error(Error::StsBadArg, error_message); 221 } 222 // make sure mean is correct if not empty 223 if(!mean.empty() && (mean.total() != (size_t) W.rows)) { 224 String error_message = format("Wrong mean shape for the given eigenvector matrix. Expected %d, but was %d.", W.cols, mean.total()); 225 CV_Error(Error::StsBadArg, error_message); 226 } 227 // initialize temporary matrices 228 Mat X, Y; 229 // copy data & make sure we are using the correct type 230 src.convertTo(Y, W.type()); 231 // calculate the reconstruction 232 gemm(Y, W, 1.0, Mat(), 0.0, X, GEMM_2_T); 233 // safe to do because of above assertion 234 if(!mean.empty()) { 235 for(int i=0; i<n; i++) { 236 Mat r_i = X.row(i); 237 add(r_i, mean.reshape(1,1), r_i); 238 } 239 } 240 return X; 241 } 242 243 244 class EigenvalueDecomposition { 245 private: 246 247 // Holds the data dimension. 248 int n; 249 250 // Stores real/imag part of a complex division. 251 double cdivr, cdivi; 252 253 // Pointer to internal memory. 254 double *d, *e, *ort; 255 double **V, **H; 256 257 // Holds the computed eigenvalues. 258 Mat _eigenvalues; 259 260 // Holds the computed eigenvectors. 261 Mat _eigenvectors; 262 263 // Allocates memory. 264 template<typename _Tp> 265 _Tp *alloc_1d(int m) { 266 return new _Tp[m]; 267 } 268 269 // Allocates memory. 270 template<typename _Tp> 271 _Tp *alloc_1d(int m, _Tp val) { 272 _Tp *arr = alloc_1d<_Tp> (m); 273 for (int i = 0; i < m; i++) 274 arr[i] = val; 275 return arr; 276 } 277 278 // Allocates memory. 279 template<typename _Tp> 280 _Tp **alloc_2d(int m, int _n) { 281 _Tp **arr = new _Tp*[m]; 282 for (int i = 0; i < m; i++) 283 arr[i] = new _Tp[_n]; 284 return arr; 285 } 286 287 // Allocates memory. 288 template<typename _Tp> 289 _Tp **alloc_2d(int m, int _n, _Tp val) { 290 _Tp **arr = alloc_2d<_Tp> (m, _n); 291 for (int i = 0; i < m; i++) { 292 for (int j = 0; j < _n; j++) { 293 arr[i][j] = val; 294 } 295 } 296 return arr; 297 } 298 299 void cdiv(double xr, double xi, double yr, double yi) { 300 double r, dv; 301 if (std::abs(yr) > std::abs(yi)) { 302 r = yi / yr; 303 dv = yr + r * yi; 304 cdivr = (xr + r * xi) / dv; 305 cdivi = (xi - r * xr) / dv; 306 } else { 307 r = yr / yi; 308 dv = yi + r * yr; 309 cdivr = (r * xr + xi) / dv; 310 cdivi = (r * xi - xr) / dv; 311 } 312 } 313 314 // Nonsymmetric reduction from Hessenberg to real Schur form. 315 316 void hqr2() { 317 318 // This is derived from the Algol procedure hqr2, 319 // by Martin and Wilkinson, Handbook for Auto. Comp., 320 // Vol.ii-Linear Algebra, and the corresponding 321 // Fortran subroutine in EISPACK. 322 323 // Initialize 324 int nn = this->n; 325 int n1 = nn - 1; 326 int low = 0; 327 int high = nn - 1; 328 double eps = std::pow(2.0, -52.0); 329 double exshift = 0.0; 330 double p = 0, q = 0, r = 0, s = 0, z = 0, t, w, x, y; 331 332 // Store roots isolated by balanc and compute matrix norm 333 334 double norm = 0.0; 335 for (int i = 0; i < nn; i++) { 336 if (i < low || i > high) { 337 d[i] = H[i][i]; 338 e[i] = 0.0; 339 } 340 for (int j = std::max(i - 1, 0); j < nn; j++) { 341 norm = norm + std::abs(H[i][j]); 342 } 343 } 344 345 // Outer loop over eigenvalue index 346 int iter = 0; 347 while (n1 >= low) { 348 349 // Look for single small sub-diagonal element 350 int l = n1; 351 while (l > low) { 352 s = std::abs(H[l - 1][l - 1]) + std::abs(H[l][l]); 353 if (s == 0.0) { 354 s = norm; 355 } 356 if (std::abs(H[l][l - 1]) < eps * s) { 357 break; 358 } 359 l--; 360 } 361 362 // Check for convergence 363 // One root found 364 365 if (l == n1) { 366 H[n1][n1] = H[n1][n1] + exshift; 367 d[n1] = H[n1][n1]; 368 e[n1] = 0.0; 369 n1--; 370 iter = 0; 371 372 // Two roots found 373 374 } else if (l == n1 - 1) { 375 w = H[n1][n1 - 1] * H[n1 - 1][n1]; 376 p = (H[n1 - 1][n1 - 1] - H[n1][n1]) / 2.0; 377 q = p * p + w; 378 z = std::sqrt(std::abs(q)); 379 H[n1][n1] = H[n1][n1] + exshift; 380 H[n1 - 1][n1 - 1] = H[n1 - 1][n1 - 1] + exshift; 381 x = H[n1][n1]; 382 383 // Real pair 384 385 if (q >= 0) { 386 if (p >= 0) { 387 z = p + z; 388 } else { 389 z = p - z; 390 } 391 d[n1 - 1] = x + z; 392 d[n1] = d[n1 - 1]; 393 if (z != 0.0) { 394 d[n1] = x - w / z; 395 } 396 e[n1 - 1] = 0.0; 397 e[n1] = 0.0; 398 x = H[n1][n1 - 1]; 399 s = std::abs(x) + std::abs(z); 400 p = x / s; 401 q = z / s; 402 r = std::sqrt(p * p + q * q); 403 p = p / r; 404 q = q / r; 405 406 // Row modification 407 408 for (int j = n1 - 1; j < nn; j++) { 409 z = H[n1 - 1][j]; 410 H[n1 - 1][j] = q * z + p * H[n1][j]; 411 H[n1][j] = q * H[n1][j] - p * z; 412 } 413 414 // Column modification 415 416 for (int i = 0; i <= n1; i++) { 417 z = H[i][n1 - 1]; 418 H[i][n1 - 1] = q * z + p * H[i][n1]; 419 H[i][n1] = q * H[i][n1] - p * z; 420 } 421 422 // Accumulate transformations 423 424 for (int i = low; i <= high; i++) { 425 z = V[i][n1 - 1]; 426 V[i][n1 - 1] = q * z + p * V[i][n1]; 427 V[i][n1] = q * V[i][n1] - p * z; 428 } 429 430 // Complex pair 431 432 } else { 433 d[n1 - 1] = x + p; 434 d[n1] = x + p; 435 e[n1 - 1] = z; 436 e[n1] = -z; 437 } 438 n1 = n1 - 2; 439 iter = 0; 440 441 // No convergence yet 442 443 } else { 444 445 // Form shift 446 447 x = H[n1][n1]; 448 y = 0.0; 449 w = 0.0; 450 if (l < n1) { 451 y = H[n1 - 1][n1 - 1]; 452 w = H[n1][n1 - 1] * H[n1 - 1][n1]; 453 } 454 455 // Wilkinson's original ad hoc shift 456 457 if (iter == 10) { 458 exshift += x; 459 for (int i = low; i <= n1; i++) { 460 H[i][i] -= x; 461 } 462 s = std::abs(H[n1][n1 - 1]) + std::abs(H[n1 - 1][n1 - 2]); 463 x = y = 0.75 * s; 464 w = -0.4375 * s * s; 465 } 466 467 // MATLAB's new ad hoc shift 468 469 if (iter == 30) { 470 s = (y - x) / 2.0; 471 s = s * s + w; 472 if (s > 0) { 473 s = std::sqrt(s); 474 if (y < x) { 475 s = -s; 476 } 477 s = x - w / ((y - x) / 2.0 + s); 478 for (int i = low; i <= n1; i++) { 479 H[i][i] -= s; 480 } 481 exshift += s; 482 x = y = w = 0.964; 483 } 484 } 485 486 iter = iter + 1; // (Could check iteration count here.) 487 488 // Look for two consecutive small sub-diagonal elements 489 int m = n1 - 2; 490 while (m >= l) { 491 z = H[m][m]; 492 r = x - z; 493 s = y - z; 494 p = (r * s - w) / H[m + 1][m] + H[m][m + 1]; 495 q = H[m + 1][m + 1] - z - r - s; 496 r = H[m + 2][m + 1]; 497 s = std::abs(p) + std::abs(q) + std::abs(r); 498 p = p / s; 499 q = q / s; 500 r = r / s; 501 if (m == l) { 502 break; 503 } 504 if (std::abs(H[m][m - 1]) * (std::abs(q) + std::abs(r)) < eps * (std::abs(p) 505 * (std::abs(H[m - 1][m - 1]) + std::abs(z) + std::abs( 506 H[m + 1][m + 1])))) { 507 break; 508 } 509 m--; 510 } 511 512 for (int i = m + 2; i <= n1; i++) { 513 H[i][i - 2] = 0.0; 514 if (i > m + 2) { 515 H[i][i - 3] = 0.0; 516 } 517 } 518 519 // Double QR step involving rows l:n and columns m:n 520 521 for (int k = m; k <= n1 - 1; k++) { 522 bool notlast = (k != n1 - 1); 523 if (k != m) { 524 p = H[k][k - 1]; 525 q = H[k + 1][k - 1]; 526 r = (notlast ? H[k + 2][k - 1] : 0.0); 527 x = std::abs(p) + std::abs(q) + std::abs(r); 528 if (x != 0.0) { 529 p = p / x; 530 q = q / x; 531 r = r / x; 532 } 533 } 534 if (x == 0.0) { 535 break; 536 } 537 s = std::sqrt(p * p + q * q + r * r); 538 if (p < 0) { 539 s = -s; 540 } 541 if (s != 0) { 542 if (k != m) { 543 H[k][k - 1] = -s * x; 544 } else if (l != m) { 545 H[k][k - 1] = -H[k][k - 1]; 546 } 547 p = p + s; 548 x = p / s; 549 y = q / s; 550 z = r / s; 551 q = q / p; 552 r = r / p; 553 554 // Row modification 555 556 for (int j = k; j < nn; j++) { 557 p = H[k][j] + q * H[k + 1][j]; 558 if (notlast) { 559 p = p + r * H[k + 2][j]; 560 H[k + 2][j] = H[k + 2][j] - p * z; 561 } 562 H[k][j] = H[k][j] - p * x; 563 H[k + 1][j] = H[k + 1][j] - p * y; 564 } 565 566 // Column modification 567 568 for (int i = 0; i <= std::min(n1, k + 3); i++) { 569 p = x * H[i][k] + y * H[i][k + 1]; 570 if (notlast) { 571 p = p + z * H[i][k + 2]; 572 H[i][k + 2] = H[i][k + 2] - p * r; 573 } 574 H[i][k] = H[i][k] - p; 575 H[i][k + 1] = H[i][k + 1] - p * q; 576 } 577 578 // Accumulate transformations 579 580 for (int i = low; i <= high; i++) { 581 p = x * V[i][k] + y * V[i][k + 1]; 582 if (notlast) { 583 p = p + z * V[i][k + 2]; 584 V[i][k + 2] = V[i][k + 2] - p * r; 585 } 586 V[i][k] = V[i][k] - p; 587 V[i][k + 1] = V[i][k + 1] - p * q; 588 } 589 } // (s != 0) 590 } // k loop 591 } // check convergence 592 } // while (n1 >= low) 593 594 // Backsubstitute to find vectors of upper triangular form 595 596 if (norm == 0.0) { 597 return; 598 } 599 600 for (n1 = nn - 1; n1 >= 0; n1--) { 601 p = d[n1]; 602 q = e[n1]; 603 604 // Real vector 605 606 if (q == 0) { 607 int l = n1; 608 H[n1][n1] = 1.0; 609 for (int i = n1 - 1; i >= 0; i--) { 610 w = H[i][i] - p; 611 r = 0.0; 612 for (int j = l; j <= n1; j++) { 613 r = r + H[i][j] * H[j][n1]; 614 } 615 if (e[i] < 0.0) { 616 z = w; 617 s = r; 618 } else { 619 l = i; 620 if (e[i] == 0.0) { 621 if (w != 0.0) { 622 H[i][n1] = -r / w; 623 } else { 624 H[i][n1] = -r / (eps * norm); 625 } 626 627 // Solve real equations 628 629 } else { 630 x = H[i][i + 1]; 631 y = H[i + 1][i]; 632 q = (d[i] - p) * (d[i] - p) + e[i] * e[i]; 633 t = (x * s - z * r) / q; 634 H[i][n1] = t; 635 if (std::abs(x) > std::abs(z)) { 636 H[i + 1][n1] = (-r - w * t) / x; 637 } else { 638 H[i + 1][n1] = (-s - y * t) / z; 639 } 640 } 641 642 // Overflow control 643 644 t = std::abs(H[i][n1]); 645 if ((eps * t) * t > 1) { 646 for (int j = i; j <= n1; j++) { 647 H[j][n1] = H[j][n1] / t; 648 } 649 } 650 } 651 } 652 // Complex vector 653 } else if (q < 0) { 654 int l = n1 - 1; 655 656 // Last vector component imaginary so matrix is triangular 657 658 if (std::abs(H[n1][n1 - 1]) > std::abs(H[n1 - 1][n1])) { 659 H[n1 - 1][n1 - 1] = q / H[n1][n1 - 1]; 660 H[n1 - 1][n1] = -(H[n1][n1] - p) / H[n1][n1 - 1]; 661 } else { 662 cdiv(0.0, -H[n1 - 1][n1], H[n1 - 1][n1 - 1] - p, q); 663 H[n1 - 1][n1 - 1] = cdivr; 664 H[n1 - 1][n1] = cdivi; 665 } 666 H[n1][n1 - 1] = 0.0; 667 H[n1][n1] = 1.0; 668 for (int i = n1 - 2; i >= 0; i--) { 669 double ra, sa, vr, vi; 670 ra = 0.0; 671 sa = 0.0; 672 for (int j = l; j <= n1; j++) { 673 ra = ra + H[i][j] * H[j][n1 - 1]; 674 sa = sa + H[i][j] * H[j][n1]; 675 } 676 w = H[i][i] - p; 677 678 if (e[i] < 0.0) { 679 z = w; 680 r = ra; 681 s = sa; 682 } else { 683 l = i; 684 if (e[i] == 0) { 685 cdiv(-ra, -sa, w, q); 686 H[i][n1 - 1] = cdivr; 687 H[i][n1] = cdivi; 688 } else { 689 690 // Solve complex equations 691 692 x = H[i][i + 1]; 693 y = H[i + 1][i]; 694 vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q; 695 vi = (d[i] - p) * 2.0 * q; 696 if (vr == 0.0 && vi == 0.0) { 697 vr = eps * norm * (std::abs(w) + std::abs(q) + std::abs(x) 698 + std::abs(y) + std::abs(z)); 699 } 700 cdiv(x * r - z * ra + q * sa, 701 x * s - z * sa - q * ra, vr, vi); 702 H[i][n1 - 1] = cdivr; 703 H[i][n1] = cdivi; 704 if (std::abs(x) > (std::abs(z) + std::abs(q))) { 705 H[i + 1][n1 - 1] = (-ra - w * H[i][n1 - 1] + q 706 * H[i][n1]) / x; 707 H[i + 1][n1] = (-sa - w * H[i][n1] - q * H[i][n1 708 - 1]) / x; 709 } else { 710 cdiv(-r - y * H[i][n1 - 1], -s - y * H[i][n1], z, 711 q); 712 H[i + 1][n1 - 1] = cdivr; 713 H[i + 1][n1] = cdivi; 714 } 715 } 716 717 // Overflow control 718 719 t = std::max(std::abs(H[i][n1 - 1]), std::abs(H[i][n1])); 720 if ((eps * t) * t > 1) { 721 for (int j = i; j <= n1; j++) { 722 H[j][n1 - 1] = H[j][n1 - 1] / t; 723 H[j][n1] = H[j][n1] / t; 724 } 725 } 726 } 727 } 728 } 729 } 730 731 // Vectors of isolated roots 732 733 for (int i = 0; i < nn; i++) { 734 if (i < low || i > high) { 735 for (int j = i; j < nn; j++) { 736 V[i][j] = H[i][j]; 737 } 738 } 739 } 740 741 // Back transformation to get eigenvectors of original matrix 742 743 for (int j = nn - 1; j >= low; j--) { 744 for (int i = low; i <= high; i++) { 745 z = 0.0; 746 for (int k = low; k <= std::min(j, high); k++) { 747 z = z + V[i][k] * H[k][j]; 748 } 749 V[i][j] = z; 750 } 751 } 752 } 753 754 // Nonsymmetric reduction to Hessenberg form. 755 void orthes() { 756 // This is derived from the Algol procedures orthes and ortran, 757 // by Martin and Wilkinson, Handbook for Auto. Comp., 758 // Vol.ii-Linear Algebra, and the corresponding 759 // Fortran subroutines in EISPACK. 760 int low = 0; 761 int high = n - 1; 762 763 for (int m = low + 1; m <= high - 1; m++) { 764 765 // Scale column. 766 767 double scale = 0.0; 768 for (int i = m; i <= high; i++) { 769 scale = scale + std::abs(H[i][m - 1]); 770 } 771 if (scale != 0.0) { 772 773 // Compute Householder transformation. 774 775 double h = 0.0; 776 for (int i = high; i >= m; i--) { 777 ort[i] = H[i][m - 1] / scale; 778 h += ort[i] * ort[i]; 779 } 780 double g = std::sqrt(h); 781 if (ort[m] > 0) { 782 g = -g; 783 } 784 h = h - ort[m] * g; 785 ort[m] = ort[m] - g; 786 787 // Apply Householder similarity transformation 788 // H = (I-u*u'/h)*H*(I-u*u')/h) 789 790 for (int j = m; j < n; j++) { 791 double f = 0.0; 792 for (int i = high; i >= m; i--) { 793 f += ort[i] * H[i][j]; 794 } 795 f = f / h; 796 for (int i = m; i <= high; i++) { 797 H[i][j] -= f * ort[i]; 798 } 799 } 800 801 for (int i = 0; i <= high; i++) { 802 double f = 0.0; 803 for (int j = high; j >= m; j--) { 804 f += ort[j] * H[i][j]; 805 } 806 f = f / h; 807 for (int j = m; j <= high; j++) { 808 H[i][j] -= f * ort[j]; 809 } 810 } 811 ort[m] = scale * ort[m]; 812 H[m][m - 1] = scale * g; 813 } 814 } 815 816 // Accumulate transformations (Algol's ortran). 817 818 for (int i = 0; i < n; i++) { 819 for (int j = 0; j < n; j++) { 820 V[i][j] = (i == j ? 1.0 : 0.0); 821 } 822 } 823 824 for (int m = high - 1; m >= low + 1; m--) { 825 if (H[m][m - 1] != 0.0) { 826 for (int i = m + 1; i <= high; i++) { 827 ort[i] = H[i][m - 1]; 828 } 829 for (int j = m; j <= high; j++) { 830 double g = 0.0; 831 for (int i = m; i <= high; i++) { 832 g += ort[i] * V[i][j]; 833 } 834 // Double division avoids possible underflow 835 g = (g / ort[m]) / H[m][m - 1]; 836 for (int i = m; i <= high; i++) { 837 V[i][j] += g * ort[i]; 838 } 839 } 840 } 841 } 842 } 843 844 // Releases all internal working memory. 845 void release() { 846 // releases the working data 847 delete[] d; 848 delete[] e; 849 delete[] ort; 850 for (int i = 0; i < n; i++) { 851 delete[] H[i]; 852 delete[] V[i]; 853 } 854 delete[] H; 855 delete[] V; 856 } 857 858 // Computes the Eigenvalue Decomposition for a matrix given in H. 859 void compute() { 860 // Allocate memory for the working data. 861 V = alloc_2d<double> (n, n, 0.0); 862 d = alloc_1d<double> (n); 863 e = alloc_1d<double> (n); 864 ort = alloc_1d<double> (n); 865 // Reduce to Hessenberg form. 866 orthes(); 867 // Reduce Hessenberg to real Schur form. 868 hqr2(); 869 // Copy eigenvalues to OpenCV Matrix. 870 _eigenvalues.create(1, n, CV_64FC1); 871 for (int i = 0; i < n; i++) { 872 _eigenvalues.at<double> (0, i) = d[i]; 873 } 874 // Copy eigenvectors to OpenCV Matrix. 875 _eigenvectors.create(n, n, CV_64FC1); 876 for (int i = 0; i < n; i++) 877 for (int j = 0; j < n; j++) 878 _eigenvectors.at<double> (i, j) = V[i][j]; 879 // Deallocate the memory by releasing all internal working data. 880 release(); 881 } 882 883 public: 884 EigenvalueDecomposition() 885 : n(0) { } 886 887 // Initializes & computes the Eigenvalue Decomposition for a general matrix 888 // given in src. This function is a port of the EigenvalueSolver in JAMA, 889 // which has been released to public domain by The MathWorks and the 890 // National Institute of Standards and Technology (NIST). 891 EigenvalueDecomposition(InputArray src) { 892 compute(src); 893 } 894 895 // This function computes the Eigenvalue Decomposition for a general matrix 896 // given in src. This function is a port of the EigenvalueSolver in JAMA, 897 // which has been released to public domain by The MathWorks and the 898 // National Institute of Standards and Technology (NIST). 899 void compute(InputArray src) 900 { 901 if(isSymmetric(src)) { 902 // Fall back to OpenCV for a symmetric matrix! 903 cv::eigen(src, _eigenvalues, _eigenvectors); 904 } else { 905 Mat tmp; 906 // Convert the given input matrix to double. Is there any way to 907 // prevent allocating the temporary memory? Only used for copying 908 // into working memory and deallocated after. 909 src.getMat().convertTo(tmp, CV_64FC1); 910 // Get dimension of the matrix. 911 this->n = tmp.cols; 912 // Allocate the matrix data to work on. 913 this->H = alloc_2d<double> (n, n); 914 // Now safely copy the data. 915 for (int i = 0; i < tmp.rows; i++) { 916 for (int j = 0; j < tmp.cols; j++) { 917 this->H[i][j] = tmp.at<double>(i, j); 918 } 919 } 920 // Deallocates the temporary matrix before computing. 921 tmp.release(); 922 // Performs the eigenvalue decomposition of H. 923 compute(); 924 } 925 } 926 927 ~EigenvalueDecomposition() {} 928 929 // Returns the eigenvalues of the Eigenvalue Decomposition. 930 Mat eigenvalues() { return _eigenvalues; } 931 // Returns the eigenvectors of the Eigenvalue Decomposition. 932 Mat eigenvectors() { return _eigenvectors; } 933 }; 934 935 936 //------------------------------------------------------------------------------ 937 // Linear Discriminant Analysis implementation 938 //------------------------------------------------------------------------------ 939 940 LDA::LDA(int num_components) : _num_components(num_components) { } 941 942 LDA::LDA(InputArrayOfArrays src, InputArray labels, int num_components) : _num_components(num_components) 943 { 944 this->compute(src, labels); //! compute eigenvectors and eigenvalues 945 } 946 947 LDA::~LDA() {} 948 949 void LDA::save(const String& filename) const 950 { 951 FileStorage fs(filename, FileStorage::WRITE); 952 if (!fs.isOpened()) { 953 CV_Error(Error::StsError, "File can't be opened for writing!"); 954 } 955 this->save(fs); 956 fs.release(); 957 } 958 959 // Deserializes this object from a given filename. 960 void LDA::load(const String& filename) { 961 FileStorage fs(filename, FileStorage::READ); 962 if (!fs.isOpened()) 963 CV_Error(Error::StsError, "File can't be opened for writing!"); 964 this->load(fs); 965 fs.release(); 966 } 967 968 // Serializes this object to a given FileStorage. 969 void LDA::save(FileStorage& fs) const { 970 // write matrices 971 fs << "num_components" << _num_components; 972 fs << "eigenvalues" << _eigenvalues; 973 fs << "eigenvectors" << _eigenvectors; 974 } 975 976 // Deserializes this object from a given FileStorage. 977 void LDA::load(const FileStorage& fs) { 978 //read matrices 979 fs["num_components"] >> _num_components; 980 fs["eigenvalues"] >> _eigenvalues; 981 fs["eigenvectors"] >> _eigenvectors; 982 } 983 984 void LDA::lda(InputArrayOfArrays _src, InputArray _lbls) { 985 // get data 986 Mat src = _src.getMat(); 987 std::vector<int> labels; 988 // safely copy the labels 989 { 990 Mat tmp = _lbls.getMat(); 991 for(unsigned int i = 0; i < tmp.total(); i++) { 992 labels.push_back(tmp.at<int>(i)); 993 } 994 } 995 // turn into row sampled matrix 996 Mat data; 997 // ensure working matrix is double precision 998 src.convertTo(data, CV_64FC1); 999 // maps the labels, so they're ascending: [0,1,...,C] 1000 std::vector<int> mapped_labels(labels.size()); 1001 std::vector<int> num2label = remove_dups(labels); 1002 std::map<int, int> label2num; 1003 for (int i = 0; i < (int)num2label.size(); i++) 1004 label2num[num2label[i]] = i; 1005 for (size_t i = 0; i < labels.size(); i++) 1006 mapped_labels[i] = label2num[labels[i]]; 1007 // get sample size, dimension 1008 int N = data.rows; 1009 int D = data.cols; 1010 // number of unique labels 1011 int C = (int)num2label.size(); 1012 // we can't do a LDA on one class, what do you 1013 // want to separate from each other then? 1014 if(C == 1) { 1015 String error_message = "At least two classes are needed to perform a LDA. Reason: Only one class was given!"; 1016 CV_Error(Error::StsBadArg, error_message); 1017 } 1018 // throw error if less labels, than samples 1019 if (labels.size() != static_cast<size_t>(N)) { 1020 String error_message = format("The number of samples must equal the number of labels. Given %d labels, %d samples. ", labels.size(), N); 1021 CV_Error(Error::StsBadArg, error_message); 1022 } 1023 // warn if within-classes scatter matrix becomes singular 1024 if (N < D) { 1025 std::cout << "Warning: Less observations than feature dimension given!" 1026 << "Computation will probably fail." 1027 << std::endl; 1028 } 1029 // clip number of components to be a valid number 1030 if ((_num_components <= 0) || (_num_components > (C - 1))) { 1031 _num_components = (C - 1); 1032 } 1033 // holds the mean over all classes 1034 Mat meanTotal = Mat::zeros(1, D, data.type()); 1035 // holds the mean for each class 1036 std::vector<Mat> meanClass(C); 1037 std::vector<int> numClass(C); 1038 // initialize 1039 for (int i = 0; i < C; i++) { 1040 numClass[i] = 0; 1041 meanClass[i] = Mat::zeros(1, D, data.type()); //! Dx1 image vector 1042 } 1043 // calculate sums 1044 for (int i = 0; i < N; i++) { 1045 Mat instance = data.row(i); 1046 int classIdx = mapped_labels[i]; 1047 add(meanTotal, instance, meanTotal); 1048 add(meanClass[classIdx], instance, meanClass[classIdx]); 1049 numClass[classIdx]++; 1050 } 1051 // calculate total mean 1052 meanTotal.convertTo(meanTotal, meanTotal.type(), 1.0 / static_cast<double> (N)); 1053 // calculate class means 1054 for (int i = 0; i < C; i++) { 1055 meanClass[i].convertTo(meanClass[i], meanClass[i].type(), 1.0 / static_cast<double> (numClass[i])); 1056 } 1057 // subtract class means 1058 for (int i = 0; i < N; i++) { 1059 int classIdx = mapped_labels[i]; 1060 Mat instance = data.row(i); 1061 subtract(instance, meanClass[classIdx], instance); 1062 } 1063 // calculate within-classes scatter 1064 Mat Sw = Mat::zeros(D, D, data.type()); 1065 mulTransposed(data, Sw, true); 1066 // calculate between-classes scatter 1067 Mat Sb = Mat::zeros(D, D, data.type()); 1068 for (int i = 0; i < C; i++) { 1069 Mat tmp; 1070 subtract(meanClass[i], meanTotal, tmp); 1071 mulTransposed(tmp, tmp, true); 1072 add(Sb, tmp, Sb); 1073 } 1074 // invert Sw 1075 Mat Swi = Sw.inv(); 1076 // M = inv(Sw)*Sb 1077 Mat M; 1078 gemm(Swi, Sb, 1.0, Mat(), 0.0, M); 1079 EigenvalueDecomposition es(M); 1080 _eigenvalues = es.eigenvalues(); 1081 _eigenvectors = es.eigenvectors(); 1082 // reshape eigenvalues, so they are stored by column 1083 _eigenvalues = _eigenvalues.reshape(1, 1); 1084 // get sorted indices descending by their eigenvalue 1085 std::vector<int> sorted_indices = argsort(_eigenvalues, false); 1086 // now sort eigenvalues and eigenvectors accordingly 1087 _eigenvalues = sortMatrixColumnsByIndices(_eigenvalues, sorted_indices); 1088 _eigenvectors = sortMatrixColumnsByIndices(_eigenvectors, sorted_indices); 1089 // and now take only the num_components and we're out! 1090 _eigenvalues = Mat(_eigenvalues, Range::all(), Range(0, _num_components)); 1091 _eigenvectors = Mat(_eigenvectors, Range::all(), Range(0, _num_components)); 1092 } 1093 1094 void LDA::compute(InputArrayOfArrays _src, InputArray _lbls) { 1095 switch(_src.kind()) { 1096 case _InputArray::STD_VECTOR_MAT: 1097 lda(asRowMatrix(_src, CV_64FC1), _lbls); 1098 break; 1099 case _InputArray::MAT: 1100 lda(_src.getMat(), _lbls); 1101 break; 1102 default: 1103 String error_message= format("InputArray Datatype %d is not supported.", _src.kind()); 1104 CV_Error(Error::StsBadArg, error_message); 1105 break; 1106 } 1107 } 1108 1109 // Projects samples into the LDA subspace. 1110 Mat LDA::project(InputArray src) { 1111 return subspaceProject(_eigenvectors, Mat(), _dataAsRow ? src : src.getMat().t()); 1112 } 1113 1114 // Reconstructs projections from the LDA subspace. 1115 Mat LDA::reconstruct(InputArray src) { 1116 return subspaceReconstruct(_eigenvectors, Mat(), _dataAsRow ? src : src.getMat().t()); 1117 } 1118 1119 } 1120