--- /dev/null
+/*
+ * Copyright (c) 2011. Philipp Wagner <bytefish[at]gmx[dot]de>.
+ * Released to public domain under terms of the BSD Simplified license.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ * * Neither the name of the organization nor the names of its contributors
+ * may be used to endorse or promote products derived from this software
+ * without specific prior written permission.
+ *
+ * See <http://www.opensource.org/licenses/bsd-license>
+ */
+
+#include "precomp.hpp"
+#include <iostream>
+#include <map>
+#include <set>
+
+namespace cv
+{
+
+// Removes duplicate elements in a given vector.
+template<typename _Tp>
+inline std::vector<_Tp> remove_dups(const std::vector<_Tp>& src) {
+ typedef typename std::set<_Tp>::const_iterator constSetIterator;
+ typedef typename std::vector<_Tp>::const_iterator constVecIterator;
+ std::set<_Tp> set_elems;
+ for (constVecIterator it = src.begin(); it != src.end(); ++it)
+ set_elems.insert(*it);
+ std::vector<_Tp> elems;
+ for (constSetIterator it = set_elems.begin(); it != set_elems.end(); ++it)
+ elems.push_back(*it);
+ return elems;
+}
+
+static Mat argsort(InputArray _src, bool ascending=true)
+{
+ Mat src = _src.getMat();
+ if (src.rows != 1 && src.cols != 1) {
+ String error_message = "Wrong shape of input matrix! Expected a matrix with one row or column.";
+ CV_Error(Error::StsBadArg, error_message);
+ }
+ int flags = SORT_EVERY_ROW | (ascending ? SORT_ASCENDING : SORT_DESCENDING);
+ Mat sorted_indices;
+ sortIdx(src.reshape(1,1),sorted_indices,flags);
+ return sorted_indices;
+}
+
+static Mat asRowMatrix(InputArrayOfArrays src, int rtype, double alpha=1, double beta=0) {
+ // make sure the input data is a vector of matrices or vector of vector
+ if(src.kind() != _InputArray::STD_VECTOR_MAT && src.kind() != _InputArray::STD_VECTOR_VECTOR) {
+ 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<...> >).";
+ CV_Error(Error::StsBadArg, error_message);
+ }
+ // number of samples
+ size_t n = src.total();
+ // return empty matrix if no matrices given
+ if(n == 0)
+ return Mat();
+ // dimensionality of (reshaped) samples
+ size_t d = src.getMat(0).total();
+ // create data matrix
+ Mat data((int)n, (int)d, rtype);
+ // now copy data
+ for(int i = 0; i < (int)n; i++) {
+ // make sure data can be reshaped, throw exception if not!
+ if(src.getMat(i).total() != d) {
+ String error_message = format("Wrong number of elements in matrix #%d! Expected %d was %d.", i, (int)d, (int)src.getMat(i).total());
+ CV_Error(Error::StsBadArg, error_message);
+ }
+ // get a hold of the current row
+ Mat xi = data.row(i);
+ // make reshape happy by cloning for non-continuous matrices
+ if(src.getMat(i).isContinuous()) {
+ src.getMat(i).reshape(1, 1).convertTo(xi, rtype, alpha, beta);
+ } else {
+ src.getMat(i).clone().reshape(1, 1).convertTo(xi, rtype, alpha, beta);
+ }
+ }
+ return data;
+}
+
+static void sortMatrixColumnsByIndices(InputArray _src, InputArray _indices, OutputArray _dst) {
+ if(_indices.getMat().type() != CV_32SC1) {
+ CV_Error(Error::StsUnsupportedFormat, "cv::sortColumnsByIndices only works on integer indices!");
+ }
+ Mat src = _src.getMat();
+ std::vector<int> indices = _indices.getMat();
+ _dst.create(src.rows, src.cols, src.type());
+ Mat dst = _dst.getMat();
+ for(size_t idx = 0; idx < indices.size(); idx++) {
+ Mat originalCol = src.col(indices[idx]);
+ Mat sortedCol = dst.col((int)idx);
+ originalCol.copyTo(sortedCol);
+ }
+}
+
+static Mat sortMatrixColumnsByIndices(InputArray src, InputArray indices) {
+ Mat dst;
+ sortMatrixColumnsByIndices(src, indices, dst);
+ return dst;
+}
+
+
+template<typename _Tp> static bool
+isSymmetric_(InputArray src) {
+ Mat _src = src.getMat();
+ if(_src.cols != _src.rows)
+ return false;
+ for (int i = 0; i < _src.rows; i++) {
+ for (int j = 0; j < _src.cols; j++) {
+ _Tp a = _src.at<_Tp> (i, j);
+ _Tp b = _src.at<_Tp> (j, i);
+ if (a != b) {
+ return false;
+ }
+ }
+ }
+ return true;
+}
+
+template<typename _Tp> static bool
+isSymmetric_(InputArray src, double eps) {
+ Mat _src = src.getMat();
+ if(_src.cols != _src.rows)
+ return false;
+ for (int i = 0; i < _src.rows; i++) {
+ for (int j = 0; j < _src.cols; j++) {
+ _Tp a = _src.at<_Tp> (i, j);
+ _Tp b = _src.at<_Tp> (j, i);
+ if (std::abs(a - b) > eps) {
+ return false;
+ }
+ }
+ }
+ return true;
+}
+
+static bool isSymmetric(InputArray src, double eps=1e-16)
+{
+ Mat m = src.getMat();
+ switch (m.type()) {
+ case CV_8SC1: return isSymmetric_<char>(m); break;
+ case CV_8UC1:
+ return isSymmetric_<unsigned char>(m); break;
+ case CV_16SC1:
+ return isSymmetric_<short>(m); break;
+ case CV_16UC1:
+ return isSymmetric_<unsigned short>(m); break;
+ case CV_32SC1:
+ return isSymmetric_<int>(m); break;
+ case CV_32FC1:
+ return isSymmetric_<float>(m, eps); break;
+ case CV_64FC1:
+ return isSymmetric_<double>(m, eps); break;
+ default:
+ break;
+ }
+ return false;
+}
+
+
+//------------------------------------------------------------------------------
+// cv::subspaceProject
+//------------------------------------------------------------------------------
+Mat LDA::subspaceProject(InputArray _W, InputArray _mean, InputArray _src) {
+ // get data matrices
+ Mat W = _W.getMat();
+ Mat mean = _mean.getMat();
+ Mat src = _src.getMat();
+ // get number of samples and dimension
+ int n = src.rows;
+ int d = src.cols;
+ // make sure the data has the correct shape
+ if(W.rows != d) {
+ 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);
+ CV_Error(Error::StsBadArg, error_message);
+ }
+ // make sure mean is correct if not empty
+ if(!mean.empty() && (mean.total() != (size_t) d)) {
+ String error_message = format("Wrong mean shape for the given data matrix. Expected %d, but was %d.", d, mean.total());
+ CV_Error(Error::StsBadArg, error_message);
+ }
+ // create temporary matrices
+ Mat X, Y;
+ // make sure you operate on correct type
+ src.convertTo(X, W.type());
+ // safe to do, because of above assertion
+ if(!mean.empty()) {
+ for(int i=0; i<n; i++) {
+ Mat r_i = X.row(i);
+ subtract(r_i, mean.reshape(1,1), r_i);
+ }
+ }
+ // finally calculate projection as Y = (X-mean)*W
+ gemm(X, W, 1.0, Mat(), 0.0, Y);
+ return Y;
+}
+
+//------------------------------------------------------------------------------
+// cv::subspaceReconstruct
+//------------------------------------------------------------------------------
+Mat LDA::subspaceReconstruct(InputArray _W, InputArray _mean, InputArray _src)
+{
+ // get data matrices
+ Mat W = _W.getMat();
+ Mat mean = _mean.getMat();
+ Mat src = _src.getMat();
+ // get number of samples and dimension
+ int n = src.rows;
+ int d = src.cols;
+ // make sure the data has the correct shape
+ if(W.cols != d) {
+ 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);
+ CV_Error(Error::StsBadArg, error_message);
+ }
+ // make sure mean is correct if not empty
+ if(!mean.empty() && (mean.total() != (size_t) W.rows)) {
+ String error_message = format("Wrong mean shape for the given eigenvector matrix. Expected %d, but was %d.", W.cols, mean.total());
+ CV_Error(Error::StsBadArg, error_message);
+ }
+ // initialize temporary matrices
+ Mat X, Y;
+ // copy data & make sure we are using the correct type
+ src.convertTo(Y, W.type());
+ // calculate the reconstruction
+ gemm(Y, W, 1.0, Mat(), 0.0, X, GEMM_2_T);
+ // safe to do because of above assertion
+ if(!mean.empty()) {
+ for(int i=0; i<n; i++) {
+ Mat r_i = X.row(i);
+ add(r_i, mean.reshape(1,1), r_i);
+ }
+ }
+ return X;
+}
+
+
+class EigenvalueDecomposition {
+private:
+
+ // Holds the data dimension.
+ int n;
+
+ // Stores real/imag part of a complex division.
+ double cdivr, cdivi;
+
+ // Pointer to internal memory.
+ double *d, *e, *ort;
+ double **V, **H;
+
+ // Holds the computed eigenvalues.
+ Mat _eigenvalues;
+
+ // Holds the computed eigenvectors.
+ Mat _eigenvectors;
+
+ // Allocates memory.
+ template<typename _Tp>
+ _Tp *alloc_1d(int m) {
+ return new _Tp[m];
+ }
+
+ // Allocates memory.
+ template<typename _Tp>
+ _Tp *alloc_1d(int m, _Tp val) {
+ _Tp *arr = alloc_1d<_Tp> (m);
+ for (int i = 0; i < m; i++)
+ arr[i] = val;
+ return arr;
+ }
+
+ // Allocates memory.
+ template<typename _Tp>
+ _Tp **alloc_2d(int m, int _n) {
+ _Tp **arr = new _Tp*[m];
+ for (int i = 0; i < m; i++)
+ arr[i] = new _Tp[_n];
+ return arr;
+ }
+
+ // Allocates memory.
+ template<typename _Tp>
+ _Tp **alloc_2d(int m, int _n, _Tp val) {
+ _Tp **arr = alloc_2d<_Tp> (m, _n);
+ for (int i = 0; i < m; i++) {
+ for (int j = 0; j < _n; j++) {
+ arr[i][j] = val;
+ }
+ }
+ return arr;
+ }
+
+ void cdiv(double xr, double xi, double yr, double yi) {
+ double r, dv;
+ if (std::abs(yr) > std::abs(yi)) {
+ r = yi / yr;
+ dv = yr + r * yi;
+ cdivr = (xr + r * xi) / dv;
+ cdivi = (xi - r * xr) / dv;
+ } else {
+ r = yr / yi;
+ dv = yi + r * yr;
+ cdivr = (r * xr + xi) / dv;
+ cdivi = (r * xi - xr) / dv;
+ }
+ }
+
+ // Nonsymmetric reduction from Hessenberg to real Schur form.
+
+ void hqr2() {
+
+ // This is derived from the Algol procedure hqr2,
+ // by Martin and Wilkinson, Handbook for Auto. Comp.,
+ // Vol.ii-Linear Algebra, and the corresponding
+ // Fortran subroutine in EISPACK.
+
+ // Initialize
+ int nn = this->n;
+ int n1 = nn - 1;
+ int low = 0;
+ int high = nn - 1;
+ double eps = std::pow(2.0, -52.0);
+ double exshift = 0.0;
+ double p = 0, q = 0, r = 0, s = 0, z = 0, t, w, x, y;
+
+ // Store roots isolated by balanc and compute matrix norm
+
+ double norm = 0.0;
+ for (int i = 0; i < nn; i++) {
+ if (i < low || i > high) {
+ d[i] = H[i][i];
+ e[i] = 0.0;
+ }
+ for (int j = std::max(i - 1, 0); j < nn; j++) {
+ norm = norm + std::abs(H[i][j]);
+ }
+ }
+
+ // Outer loop over eigenvalue index
+ int iter = 0;
+ while (n1 >= low) {
+
+ // Look for single small sub-diagonal element
+ int l = n1;
+ while (l > low) {
+ s = std::abs(H[l - 1][l - 1]) + std::abs(H[l][l]);
+ if (s == 0.0) {
+ s = norm;
+ }
+ if (std::abs(H[l][l - 1]) < eps * s) {
+ break;
+ }
+ l--;
+ }
+
+ // Check for convergence
+ // One root found
+
+ if (l == n1) {
+ H[n1][n1] = H[n1][n1] + exshift;
+ d[n1] = H[n1][n1];
+ e[n1] = 0.0;
+ n1--;
+ iter = 0;
+
+ // Two roots found
+
+ } else if (l == n1 - 1) {
+ w = H[n1][n1 - 1] * H[n1 - 1][n1];
+ p = (H[n1 - 1][n1 - 1] - H[n1][n1]) / 2.0;
+ q = p * p + w;
+ z = std::sqrt(std::abs(q));
+ H[n1][n1] = H[n1][n1] + exshift;
+ H[n1 - 1][n1 - 1] = H[n1 - 1][n1 - 1] + exshift;
+ x = H[n1][n1];
+
+ // Real pair
+
+ if (q >= 0) {
+ if (p >= 0) {
+ z = p + z;
+ } else {
+ z = p - z;
+ }
+ d[n1 - 1] = x + z;
+ d[n1] = d[n1 - 1];
+ if (z != 0.0) {
+ d[n1] = x - w / z;
+ }
+ e[n1 - 1] = 0.0;
+ e[n1] = 0.0;
+ x = H[n1][n1 - 1];
+ s = std::abs(x) + std::abs(z);
+ p = x / s;
+ q = z / s;
+ r = std::sqrt(p * p + q * q);
+ p = p / r;
+ q = q / r;
+
+ // Row modification
+
+ for (int j = n1 - 1; j < nn; j++) {
+ z = H[n1 - 1][j];
+ H[n1 - 1][j] = q * z + p * H[n1][j];
+ H[n1][j] = q * H[n1][j] - p * z;
+ }
+
+ // Column modification
+
+ for (int i = 0; i <= n1; i++) {
+ z = H[i][n1 - 1];
+ H[i][n1 - 1] = q * z + p * H[i][n1];
+ H[i][n1] = q * H[i][n1] - p * z;
+ }
+
+ // Accumulate transformations
+
+ for (int i = low; i <= high; i++) {
+ z = V[i][n1 - 1];
+ V[i][n1 - 1] = q * z + p * V[i][n1];
+ V[i][n1] = q * V[i][n1] - p * z;
+ }
+
+ // Complex pair
+
+ } else {
+ d[n1 - 1] = x + p;
+ d[n1] = x + p;
+ e[n1 - 1] = z;
+ e[n1] = -z;
+ }
+ n1 = n1 - 2;
+ iter = 0;
+
+ // No convergence yet
+
+ } else {
+
+ // Form shift
+
+ x = H[n1][n1];
+ y = 0.0;
+ w = 0.0;
+ if (l < n1) {
+ y = H[n1 - 1][n1 - 1];
+ w = H[n1][n1 - 1] * H[n1 - 1][n1];
+ }
+
+ // Wilkinson's original ad hoc shift
+
+ if (iter == 10) {
+ exshift += x;
+ for (int i = low; i <= n1; i++) {
+ H[i][i] -= x;
+ }
+ s = std::abs(H[n1][n1 - 1]) + std::abs(H[n1 - 1][n1 - 2]);
+ x = y = 0.75 * s;
+ w = -0.4375 * s * s;
+ }
+
+ // MATLAB's new ad hoc shift
+
+ if (iter == 30) {
+ s = (y - x) / 2.0;
+ s = s * s + w;
+ if (s > 0) {
+ s = std::sqrt(s);
+ if (y < x) {
+ s = -s;
+ }
+ s = x - w / ((y - x) / 2.0 + s);
+ for (int i = low; i <= n1; i++) {
+ H[i][i] -= s;
+ }
+ exshift += s;
+ x = y = w = 0.964;
+ }
+ }
+
+ iter = iter + 1; // (Could check iteration count here.)
+
+ // Look for two consecutive small sub-diagonal elements
+ int m = n1 - 2;
+ while (m >= l) {
+ z = H[m][m];
+ r = x - z;
+ s = y - z;
+ p = (r * s - w) / H[m + 1][m] + H[m][m + 1];
+ q = H[m + 1][m + 1] - z - r - s;
+ r = H[m + 2][m + 1];
+ s = std::abs(p) + std::abs(q) + std::abs(r);
+ p = p / s;
+ q = q / s;
+ r = r / s;
+ if (m == l) {
+ break;
+ }
+ 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])))) {
+ break;
+ }
+ m--;
+ }
+
+ for (int i = m + 2; i <= n1; i++) {
+ H[i][i - 2] = 0.0;
+ if (i > m + 2) {
+ H[i][i - 3] = 0.0;
+ }
+ }
+
+ // Double QR step involving rows l:n and columns m:n
+
+ for (int k = m; k <= n1 - 1; k++) {
+ bool notlast = (k != n1 - 1);
+ if (k != m) {
+ p = H[k][k - 1];
+ q = H[k + 1][k - 1];
+ r = (notlast ? H[k + 2][k - 1] : 0.0);
+ x = std::abs(p) + std::abs(q) + std::abs(r);
+ if (x != 0.0) {
+ p = p / x;
+ q = q / x;
+ r = r / x;
+ }
+ }
+ if (x == 0.0) {
+ break;
+ }
+ s = std::sqrt(p * p + q * q + r * r);
+ if (p < 0) {
+ s = -s;
+ }
+ if (s != 0) {
+ if (k != m) {
+ H[k][k - 1] = -s * x;
+ } else if (l != m) {
+ H[k][k - 1] = -H[k][k - 1];
+ }
+ p = p + s;
+ x = p / s;
+ y = q / s;
+ z = r / s;
+ q = q / p;
+ r = r / p;
+
+ // Row modification
+
+ for (int j = k; j < nn; j++) {
+ p = H[k][j] + q * H[k + 1][j];
+ if (notlast) {
+ p = p + r * H[k + 2][j];
+ H[k + 2][j] = H[k + 2][j] - p * z;
+ }
+ H[k][j] = H[k][j] - p * x;
+ H[k + 1][j] = H[k + 1][j] - p * y;
+ }
+
+ // Column modification
+
+ for (int i = 0; i <= std::min(n1, k + 3); i++) {
+ p = x * H[i][k] + y * H[i][k + 1];
+ if (notlast) {
+ p = p + z * H[i][k + 2];
+ H[i][k + 2] = H[i][k + 2] - p * r;
+ }
+ H[i][k] = H[i][k] - p;
+ H[i][k + 1] = H[i][k + 1] - p * q;
+ }
+
+ // Accumulate transformations
+
+ for (int i = low; i <= high; i++) {
+ p = x * V[i][k] + y * V[i][k + 1];
+ if (notlast) {
+ p = p + z * V[i][k + 2];
+ V[i][k + 2] = V[i][k + 2] - p * r;
+ }
+ V[i][k] = V[i][k] - p;
+ V[i][k + 1] = V[i][k + 1] - p * q;
+ }
+ } // (s != 0)
+ } // k loop
+ } // check convergence
+ } // while (n1 >= low)
+
+ // Backsubstitute to find vectors of upper triangular form
+
+ if (norm == 0.0) {
+ return;
+ }
+
+ for (n1 = nn - 1; n1 >= 0; n1--) {
+ p = d[n1];
+ q = e[n1];
+
+ // Real vector
+
+ if (q == 0) {
+ int l = n1;
+ H[n1][n1] = 1.0;
+ for (int i = n1 - 1; i >= 0; i--) {
+ w = H[i][i] - p;
+ r = 0.0;
+ for (int j = l; j <= n1; j++) {
+ r = r + H[i][j] * H[j][n1];
+ }
+ if (e[i] < 0.0) {
+ z = w;
+ s = r;
+ } else {
+ l = i;
+ if (e[i] == 0.0) {
+ if (w != 0.0) {
+ H[i][n1] = -r / w;
+ } else {
+ H[i][n1] = -r / (eps * norm);
+ }
+
+ // Solve real equations
+
+ } else {
+ x = H[i][i + 1];
+ y = H[i + 1][i];
+ q = (d[i] - p) * (d[i] - p) + e[i] * e[i];
+ t = (x * s - z * r) / q;
+ H[i][n1] = t;
+ if (std::abs(x) > std::abs(z)) {
+ H[i + 1][n1] = (-r - w * t) / x;
+ } else {
+ H[i + 1][n1] = (-s - y * t) / z;
+ }
+ }
+
+ // Overflow control
+
+ t = std::abs(H[i][n1]);
+ if ((eps * t) * t > 1) {
+ for (int j = i; j <= n1; j++) {
+ H[j][n1] = H[j][n1] / t;
+ }
+ }
+ }
+ }
+ // Complex vector
+ } else if (q < 0) {
+ int l = n1 - 1;
+
+ // Last vector component imaginary so matrix is triangular
+
+ if (std::abs(H[n1][n1 - 1]) > std::abs(H[n1 - 1][n1])) {
+ H[n1 - 1][n1 - 1] = q / H[n1][n1 - 1];
+ H[n1 - 1][n1] = -(H[n1][n1] - p) / H[n1][n1 - 1];
+ } else {
+ cdiv(0.0, -H[n1 - 1][n1], H[n1 - 1][n1 - 1] - p, q);
+ H[n1 - 1][n1 - 1] = cdivr;
+ H[n1 - 1][n1] = cdivi;
+ }
+ H[n1][n1 - 1] = 0.0;
+ H[n1][n1] = 1.0;
+ for (int i = n1 - 2; i >= 0; i--) {
+ double ra, sa, vr, vi;
+ ra = 0.0;
+ sa = 0.0;
+ for (int j = l; j <= n1; j++) {
+ ra = ra + H[i][j] * H[j][n1 - 1];
+ sa = sa + H[i][j] * H[j][n1];
+ }
+ w = H[i][i] - p;
+
+ if (e[i] < 0.0) {
+ z = w;
+ r = ra;
+ s = sa;
+ } else {
+ l = i;
+ if (e[i] == 0) {
+ cdiv(-ra, -sa, w, q);
+ H[i][n1 - 1] = cdivr;
+ H[i][n1] = cdivi;
+ } else {
+
+ // Solve complex equations
+
+ x = H[i][i + 1];
+ y = H[i + 1][i];
+ vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
+ vi = (d[i] - p) * 2.0 * q;
+ if (vr == 0.0 && vi == 0.0) {
+ vr = eps * norm * (std::abs(w) + std::abs(q) + std::abs(x)
+ + std::abs(y) + std::abs(z));
+ }
+ cdiv(x * r - z * ra + q * sa,
+ x * s - z * sa - q * ra, vr, vi);
+ H[i][n1 - 1] = cdivr;
+ H[i][n1] = cdivi;
+ if (std::abs(x) > (std::abs(z) + std::abs(q))) {
+ H[i + 1][n1 - 1] = (-ra - w * H[i][n1 - 1] + q
+ * H[i][n1]) / x;
+ H[i + 1][n1] = (-sa - w * H[i][n1] - q * H[i][n1
+ - 1]) / x;
+ } else {
+ cdiv(-r - y * H[i][n1 - 1], -s - y * H[i][n1], z,
+ q);
+ H[i + 1][n1 - 1] = cdivr;
+ H[i + 1][n1] = cdivi;
+ }
+ }
+
+ // Overflow control
+
+ t = std::max(std::abs(H[i][n1 - 1]), std::abs(H[i][n1]));
+ if ((eps * t) * t > 1) {
+ for (int j = i; j <= n1; j++) {
+ H[j][n1 - 1] = H[j][n1 - 1] / t;
+ H[j][n1] = H[j][n1] / t;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // Vectors of isolated roots
+
+ for (int i = 0; i < nn; i++) {
+ if (i < low || i > high) {
+ for (int j = i; j < nn; j++) {
+ V[i][j] = H[i][j];
+ }
+ }
+ }
+
+ // Back transformation to get eigenvectors of original matrix
+
+ for (int j = nn - 1; j >= low; j--) {
+ for (int i = low; i <= high; i++) {
+ z = 0.0;
+ for (int k = low; k <= std::min(j, high); k++) {
+ z = z + V[i][k] * H[k][j];
+ }
+ V[i][j] = z;
+ }
+ }
+ }
+
+ // Nonsymmetric reduction to Hessenberg form.
+ void orthes() {
+ // This is derived from the Algol procedures orthes and ortran,
+ // by Martin and Wilkinson, Handbook for Auto. Comp.,
+ // Vol.ii-Linear Algebra, and the corresponding
+ // Fortran subroutines in EISPACK.
+ int low = 0;
+ int high = n - 1;
+
+ for (int m = low + 1; m <= high - 1; m++) {
+
+ // Scale column.
+
+ double scale = 0.0;
+ for (int i = m; i <= high; i++) {
+ scale = scale + std::abs(H[i][m - 1]);
+ }
+ if (scale != 0.0) {
+
+ // Compute Householder transformation.
+
+ double h = 0.0;
+ for (int i = high; i >= m; i--) {
+ ort[i] = H[i][m - 1] / scale;
+ h += ort[i] * ort[i];
+ }
+ double g = std::sqrt(h);
+ if (ort[m] > 0) {
+ g = -g;
+ }
+ h = h - ort[m] * g;
+ ort[m] = ort[m] - g;
+
+ // Apply Householder similarity transformation
+ // H = (I-u*u'/h)*H*(I-u*u')/h)
+
+ for (int j = m; j < n; j++) {
+ double f = 0.0;
+ for (int i = high; i >= m; i--) {
+ f += ort[i] * H[i][j];
+ }
+ f = f / h;
+ for (int i = m; i <= high; i++) {
+ H[i][j] -= f * ort[i];
+ }
+ }
+
+ for (int i = 0; i <= high; i++) {
+ double f = 0.0;
+ for (int j = high; j >= m; j--) {
+ f += ort[j] * H[i][j];
+ }
+ f = f / h;
+ for (int j = m; j <= high; j++) {
+ H[i][j] -= f * ort[j];
+ }
+ }
+ ort[m] = scale * ort[m];
+ H[m][m - 1] = scale * g;
+ }
+ }
+
+ // Accumulate transformations (Algol's ortran).
+
+ for (int i = 0; i < n; i++) {
+ for (int j = 0; j < n; j++) {
+ V[i][j] = (i == j ? 1.0 : 0.0);
+ }
+ }
+
+ for (int m = high - 1; m >= low + 1; m--) {
+ if (H[m][m - 1] != 0.0) {
+ for (int i = m + 1; i <= high; i++) {
+ ort[i] = H[i][m - 1];
+ }
+ for (int j = m; j <= high; j++) {
+ double g = 0.0;
+ for (int i = m; i <= high; i++) {
+ g += ort[i] * V[i][j];
+ }
+ // Double division avoids possible underflow
+ g = (g / ort[m]) / H[m][m - 1];
+ for (int i = m; i <= high; i++) {
+ V[i][j] += g * ort[i];
+ }
+ }
+ }
+ }
+ }
+
+ // Releases all internal working memory.
+ void release() {
+ // releases the working data
+ delete[] d;
+ delete[] e;
+ delete[] ort;
+ for (int i = 0; i < n; i++) {
+ delete[] H[i];
+ delete[] V[i];
+ }
+ delete[] H;
+ delete[] V;
+ }
+
+ // Computes the Eigenvalue Decomposition for a matrix given in H.
+ void compute() {
+ // Allocate memory for the working data.
+ V = alloc_2d<double> (n, n, 0.0);
+ d = alloc_1d<double> (n);
+ e = alloc_1d<double> (n);
+ ort = alloc_1d<double> (n);
+ // Reduce to Hessenberg form.
+ orthes();
+ // Reduce Hessenberg to real Schur form.
+ hqr2();
+ // Copy eigenvalues to OpenCV Matrix.
+ _eigenvalues.create(1, n, CV_64FC1);
+ for (int i = 0; i < n; i++) {
+ _eigenvalues.at<double> (0, i) = d[i];
+ }
+ // Copy eigenvectors to OpenCV Matrix.
+ _eigenvectors.create(n, n, CV_64FC1);
+ for (int i = 0; i < n; i++)
+ for (int j = 0; j < n; j++)
+ _eigenvectors.at<double> (i, j) = V[i][j];
+ // Deallocate the memory by releasing all internal working data.
+ release();
+ }
+
+public:
+ EigenvalueDecomposition()
+ : n(0) { }
+
+ // Initializes & computes the Eigenvalue Decomposition for a general matrix
+ // given in src. This function is a port of the EigenvalueSolver in JAMA,
+ // which has been released to public domain by The MathWorks and the
+ // National Institute of Standards and Technology (NIST).
+ EigenvalueDecomposition(InputArray src) {
+ compute(src);
+ }
+
+ // This function computes the Eigenvalue Decomposition for a general matrix
+ // given in src. This function is a port of the EigenvalueSolver in JAMA,
+ // which has been released to public domain by The MathWorks and the
+ // National Institute of Standards and Technology (NIST).
+ void compute(InputArray src)
+ {
+ if(isSymmetric(src)) {
+ // Fall back to OpenCV for a symmetric matrix!
+ cv::eigen(src, _eigenvalues, _eigenvectors);
+ } else {
+ Mat tmp;
+ // Convert the given input matrix to double. Is there any way to
+ // prevent allocating the temporary memory? Only used for copying
+ // into working memory and deallocated after.
+ src.getMat().convertTo(tmp, CV_64FC1);
+ // Get dimension of the matrix.
+ this->n = tmp.cols;
+ // Allocate the matrix data to work on.
+ this->H = alloc_2d<double> (n, n);
+ // Now safely copy the data.
+ for (int i = 0; i < tmp.rows; i++) {
+ for (int j = 0; j < tmp.cols; j++) {
+ this->H[i][j] = tmp.at<double>(i, j);
+ }
+ }
+ // Deallocates the temporary matrix before computing.
+ tmp.release();
+ // Performs the eigenvalue decomposition of H.
+ compute();
+ }
+ }
+
+ ~EigenvalueDecomposition() {}
+
+ // Returns the eigenvalues of the Eigenvalue Decomposition.
+ Mat eigenvalues() { return _eigenvalues; }
+ // Returns the eigenvectors of the Eigenvalue Decomposition.
+ Mat eigenvectors() { return _eigenvectors; }
+};
+
+
+//------------------------------------------------------------------------------
+// Linear Discriminant Analysis implementation
+//------------------------------------------------------------------------------
+
+LDA::LDA(int num_components) : _num_components(num_components) { }
+
+LDA::LDA(InputArrayOfArrays src, InputArray labels, int num_components) : _num_components(num_components)
+{
+ this->compute(src, labels); //! compute eigenvectors and eigenvalues
+}
+
+LDA::~LDA() {}
+
+void LDA::save(const String& filename) const
+{
+ FileStorage fs(filename, FileStorage::WRITE);
+ if (!fs.isOpened()) {
+ CV_Error(Error::StsError, "File can't be opened for writing!");
+ }
+ this->save(fs);
+ fs.release();
+}
+
+// Deserializes this object from a given filename.
+void LDA::load(const String& filename) {
+ FileStorage fs(filename, FileStorage::READ);
+ if (!fs.isOpened())
+ CV_Error(Error::StsError, "File can't be opened for writing!");
+ this->load(fs);
+ fs.release();
+}
+
+// Serializes this object to a given FileStorage.
+void LDA::save(FileStorage& fs) const {
+ // write matrices
+ fs << "num_components" << _num_components;
+ fs << "eigenvalues" << _eigenvalues;
+ fs << "eigenvectors" << _eigenvectors;
+}
+
+// Deserializes this object from a given FileStorage.
+void LDA::load(const FileStorage& fs) {
+ //read matrices
+ fs["num_components"] >> _num_components;
+ fs["eigenvalues"] >> _eigenvalues;
+ fs["eigenvectors"] >> _eigenvectors;
+}
+
+void LDA::lda(InputArrayOfArrays _src, InputArray _lbls) {
+ // get data
+ Mat src = _src.getMat();
+ std::vector<int> labels;
+ // safely copy the labels
+ {
+ Mat tmp = _lbls.getMat();
+ for(unsigned int i = 0; i < tmp.total(); i++) {
+ labels.push_back(tmp.at<int>(i));
+ }
+ }
+ // turn into row sampled matrix
+ Mat data;
+ // ensure working matrix is double precision
+ src.convertTo(data, CV_64FC1);
+ // maps the labels, so they're ascending: [0,1,...,C]
+ std::vector<int> mapped_labels(labels.size());
+ std::vector<int> num2label = remove_dups(labels);
+ std::map<int, int> label2num;
+ for (int i = 0; i < (int)num2label.size(); i++)
+ label2num[num2label[i]] = i;
+ for (size_t i = 0; i < labels.size(); i++)
+ mapped_labels[i] = label2num[labels[i]];
+ // get sample size, dimension
+ int N = data.rows;
+ int D = data.cols;
+ // number of unique labels
+ int C = (int)num2label.size();
+ // we can't do a LDA on one class, what do you
+ // want to separate from each other then?
+ if(C == 1) {
+ String error_message = "At least two classes are needed to perform a LDA. Reason: Only one class was given!";
+ CV_Error(Error::StsBadArg, error_message);
+ }
+ // throw error if less labels, than samples
+ if (labels.size() != static_cast<size_t>(N)) {
+ String error_message = format("The number of samples must equal the number of labels. Given %d labels, %d samples. ", labels.size(), N);
+ CV_Error(Error::StsBadArg, error_message);
+ }
+ // warn if within-classes scatter matrix becomes singular
+ if (N < D) {
+ std::cout << "Warning: Less observations than feature dimension given!"
+ << "Computation will probably fail."
+ << std::endl;
+ }
+ // clip number of components to be a valid number
+ if ((_num_components <= 0) || (_num_components > (C - 1))) {
+ _num_components = (C - 1);
+ }
+ // holds the mean over all classes
+ Mat meanTotal = Mat::zeros(1, D, data.type());
+ // holds the mean for each class
+ std::vector<Mat> meanClass(C);
+ std::vector<int> numClass(C);
+ // initialize
+ for (int i = 0; i < C; i++) {
+ numClass[i] = 0;
+ meanClass[i] = Mat::zeros(1, D, data.type()); //! Dx1 image vector
+ }
+ // calculate sums
+ for (int i = 0; i < N; i++) {
+ Mat instance = data.row(i);
+ int classIdx = mapped_labels[i];
+ add(meanTotal, instance, meanTotal);
+ add(meanClass[classIdx], instance, meanClass[classIdx]);
+ numClass[classIdx]++;
+ }
+ // calculate total mean
+ meanTotal.convertTo(meanTotal, meanTotal.type(), 1.0 / static_cast<double> (N));
+ // calculate class means
+ for (int i = 0; i < C; i++) {
+ meanClass[i].convertTo(meanClass[i], meanClass[i].type(), 1.0 / static_cast<double> (numClass[i]));
+ }
+ // subtract class means
+ for (int i = 0; i < N; i++) {
+ int classIdx = mapped_labels[i];
+ Mat instance = data.row(i);
+ subtract(instance, meanClass[classIdx], instance);
+ }
+ // calculate within-classes scatter
+ Mat Sw = Mat::zeros(D, D, data.type());
+ mulTransposed(data, Sw, true);
+ // calculate between-classes scatter
+ Mat Sb = Mat::zeros(D, D, data.type());
+ for (int i = 0; i < C; i++) {
+ Mat tmp;
+ subtract(meanClass[i], meanTotal, tmp);
+ mulTransposed(tmp, tmp, true);
+ add(Sb, tmp, Sb);
+ }
+ // invert Sw
+ Mat Swi = Sw.inv();
+ // M = inv(Sw)*Sb
+ Mat M;
+ gemm(Swi, Sb, 1.0, Mat(), 0.0, M);
+ EigenvalueDecomposition es(M);
+ _eigenvalues = es.eigenvalues();
+ _eigenvectors = es.eigenvectors();
+ // reshape eigenvalues, so they are stored by column
+ _eigenvalues = _eigenvalues.reshape(1, 1);
+ // get sorted indices descending by their eigenvalue
+ std::vector<int> sorted_indices = argsort(_eigenvalues, false);
+ // now sort eigenvalues and eigenvectors accordingly
+ _eigenvalues = sortMatrixColumnsByIndices(_eigenvalues, sorted_indices);
+ _eigenvectors = sortMatrixColumnsByIndices(_eigenvectors, sorted_indices);
+ // and now take only the num_components and we're out!
+ _eigenvalues = Mat(_eigenvalues, Range::all(), Range(0, _num_components));
+ _eigenvectors = Mat(_eigenvectors, Range::all(), Range(0, _num_components));
+}
+
+void LDA::compute(InputArrayOfArrays _src, InputArray _lbls) {
+ switch(_src.kind()) {
+ case _InputArray::STD_VECTOR_MAT:
+ lda(asRowMatrix(_src, CV_64FC1), _lbls);
+ break;
+ case _InputArray::MAT:
+ lda(_src.getMat(), _lbls);
+ break;
+ default:
+ String error_message= format("InputArray Datatype %d is not supported.", _src.kind());
+ CV_Error(Error::StsBadArg, error_message);
+ break;
+ }
+}
+
+// Projects samples into the LDA subspace.
+Mat LDA::project(InputArray src) {
+ return subspaceProject(_eigenvectors, Mat(), _dataAsRow ? src : src.getMat().t());
+}
+
+// Reconstructs projections from the LDA subspace.
+Mat LDA::reconstruct(InputArray src) {
+ return subspaceReconstruct(_eigenvectors, Mat(), _dataAsRow ? src : src.getMat().t());
+}
+
+}
--- /dev/null
+/*
+ * Copyright (c) 2011. Philipp Wagner <bytefish[at]gmx[dot]de>.
+ * Released to public domain under terms of the BSD Simplified license.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ * * Neither the name of the organization nor the names of its contributors
+ * may be used to endorse or promote products derived from this software
+ * without specific prior written permission.
+ *
+ * See <http://www.opensource.org/licenses/bsd-license>
+ */
+#include "precomp.hpp"
+#include <iostream>
+
+#ifdef _MSC_VER
+#pragma warning( disable: 4305 )
+#endif
+
+namespace cv
+{
+
+static Mat linspace(float x0, float x1, int n)
+{
+ Mat pts(n, 1, CV_32FC1);
+ float step = (x1-x0)/(n-1);
+ for(int i = 0; i < n; i++)
+ pts.at<float>(i,0) = x0+i*step;
+ return pts;
+}
+
+//------------------------------------------------------------------------------
+// cv::sortMatrixRowsByIndices
+//------------------------------------------------------------------------------
+static void sortMatrixRowsByIndices(InputArray _src, InputArray _indices, OutputArray _dst)
+{
+ if(_indices.getMat().type() != CV_32SC1)
+ CV_Error(Error::StsUnsupportedFormat, "cv::sortRowsByIndices only works on integer indices!");
+ Mat src = _src.getMat();
+ std::vector<int> indices = _indices.getMat();
+ _dst.create(src.rows, src.cols, src.type());
+ Mat dst = _dst.getMat();
+ for(size_t idx = 0; idx < indices.size(); idx++) {
+ Mat originalRow = src.row(indices[idx]);
+ Mat sortedRow = dst.row((int)idx);
+ originalRow.copyTo(sortedRow);
+ }
+}
+
+static Mat sortMatrixRowsByIndices(InputArray src, InputArray indices)
+{
+ Mat dst;
+ sortMatrixRowsByIndices(src, indices, dst);
+ return dst;
+}
+
+
+static Mat argsort(InputArray _src, bool ascending=true)
+{
+ Mat src = _src.getMat();
+ if (src.rows != 1 && src.cols != 1)
+ CV_Error(Error::StsBadArg, "cv::argsort only sorts 1D matrices.");
+ int flags = SORT_EVERY_ROW | (ascending ? SORT_ASCENDING : SORT_DESCENDING);
+ Mat sorted_indices;
+ sortIdx(src.reshape(1,1),sorted_indices,flags);
+ return sorted_indices;
+}
+
+template <typename _Tp> static
+Mat interp1_(const Mat& X_, const Mat& Y_, const Mat& XI)
+{
+ int n = XI.rows;
+ // sort input table
+ std::vector<int> sort_indices = argsort(X_);
+
+ Mat X = sortMatrixRowsByIndices(X_,sort_indices);
+ Mat Y = sortMatrixRowsByIndices(Y_,sort_indices);
+ // interpolated values
+ Mat yi = Mat::zeros(XI.size(), XI.type());
+ for(int i = 0; i < n; i++) {
+ int c = 0;
+ int low = 0;
+ int high = X.rows - 1;
+ // set bounds
+ if(XI.at<_Tp>(i,0) < X.at<_Tp>(low, 0))
+ high = 1;
+ if(XI.at<_Tp>(i,0) > X.at<_Tp>(high, 0))
+ low = high - 1;
+ // binary search
+ while((high-low)>1) {
+ c = low + ((high - low) >> 1);
+ if(XI.at<_Tp>(i,0) > X.at<_Tp>(c,0)) {
+ low = c;
+ } else {
+ high = c;
+ }
+ }
+ // linear interpolation
+ yi.at<_Tp>(i,0) += Y.at<_Tp>(low,0)
+ + (XI.at<_Tp>(i,0) - X.at<_Tp>(low,0))
+ * (Y.at<_Tp>(high,0) - Y.at<_Tp>(low,0))
+ / (X.at<_Tp>(high,0) - X.at<_Tp>(low,0));
+ }
+ return yi;
+}
+
+static Mat interp1(InputArray _x, InputArray _Y, InputArray _xi)
+{
+ // get matrices
+ Mat x = _x.getMat();
+ Mat Y = _Y.getMat();
+ Mat xi = _xi.getMat();
+ // check types & alignment
+ CV_Assert((x.type() == Y.type()) && (Y.type() == xi.type()));
+ CV_Assert((x.cols == 1) && (x.rows == Y.rows) && (x.cols == Y.cols));
+ // call templated interp1
+ switch(x.type()) {
+ case CV_8SC1: return interp1_<char>(x,Y,xi); break;
+ case CV_8UC1: return interp1_<unsigned char>(x,Y,xi); break;
+ case CV_16SC1: return interp1_<short>(x,Y,xi); break;
+ case CV_16UC1: return interp1_<unsigned short>(x,Y,xi); break;
+ case CV_32SC1: return interp1_<int>(x,Y,xi); break;
+ case CV_32FC1: return interp1_<float>(x,Y,xi); break;
+ case CV_64FC1: return interp1_<double>(x,Y,xi); break;
+ default: CV_Error(Error::StsUnsupportedFormat, ""); break;
+ }
+ return Mat();
+}
+
+namespace colormap
+{
+
+ class ColorMap {
+
+ protected:
+ Mat _lut;
+
+ public:
+ virtual ~ColorMap() {}
+
+ // Applies the colormap on a given image.
+ //
+ // This function expects BGR-aligned data of type CV_8UC1 or
+ // CV_8UC3. If the wrong image type is given, the original image
+ // will be returned.
+ //
+ // Throws an error for wrong-aligned lookup table, which must be
+ // of size 256 in the latest OpenCV release (2.3.1).
+ void operator()(InputArray src, OutputArray dst) const;
+
+ // Setup base map to interpolate from.
+ virtual void init(int n) = 0;
+
+ // Interpolates from a base colormap.
+ static Mat linear_colormap(InputArray X,
+ InputArray r, InputArray g, InputArray b,
+ int n) {
+ return linear_colormap(X,r,g,b,linspace(0,1,n));
+ }
+
+ // Interpolates from a base colormap.
+ static Mat linear_colormap(InputArray X,
+ InputArray r, InputArray g, InputArray b,
+ float begin, float end, float n) {
+ return linear_colormap(X,r,g,b,linspace(begin,end, cvRound(n)));
+ }
+
+ // Interpolates from a base colormap.
+ static Mat linear_colormap(InputArray X,
+ InputArray r, InputArray g, InputArray b,
+ InputArray xi);
+ };
+
+ // Equals the GNU Octave colormap "autumn".
+ class Autumn : public ColorMap {
+ public:
+ Autumn() : ColorMap() {
+ init(256);
+ }
+
+ Autumn(int n) : ColorMap() {
+ init(n);
+ }
+
+ void init(int n) {
+ float r[] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
+ float g[] = { 0, 0.01587301587301587, 0.03174603174603174, 0.04761904761904762, 0.06349206349206349, 0.07936507936507936, 0.09523809523809523, 0.1111111111111111, 0.126984126984127, 0.1428571428571428, 0.1587301587301587, 0.1746031746031746, 0.1904761904761905, 0.2063492063492063, 0.2222222222222222, 0.2380952380952381, 0.253968253968254, 0.2698412698412698, 0.2857142857142857, 0.3015873015873016, 0.3174603174603174, 0.3333333333333333, 0.3492063492063492, 0.3650793650793651, 0.3809523809523809, 0.3968253968253968, 0.4126984126984127, 0.4285714285714285, 0.4444444444444444, 0.4603174603174603, 0.4761904761904762, 0.492063492063492, 0.5079365079365079, 0.5238095238095238, 0.5396825396825397, 0.5555555555555556, 0.5714285714285714, 0.5873015873015873, 0.6031746031746031, 0.6190476190476191, 0.6349206349206349, 0.6507936507936508, 0.6666666666666666, 0.6825396825396826, 0.6984126984126984, 0.7142857142857143, 0.7301587301587301, 0.746031746031746, 0.7619047619047619, 0.7777777777777778, 0.7936507936507936, 0.8095238095238095, 0.8253968253968254, 0.8412698412698413, 0.8571428571428571, 0.873015873015873, 0.8888888888888888, 0.9047619047619048, 0.9206349206349206, 0.9365079365079365, 0.9523809523809523, 0.9682539682539683, 0.9841269841269841, 1};
+ float b[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+ Mat X = linspace(0,1,64);
+ this->_lut = ColorMap::linear_colormap(X,
+ Mat(64,1, CV_32FC1, r).clone(), // red
+ Mat(64,1, CV_32FC1, g).clone(), // green
+ Mat(64,1, CV_32FC1, b).clone(), // blue
+ n); // number of sample points
+ }
+ };
+
+ // Equals the GNU Octave colormap "bone".
+ class Bone : public ColorMap {
+ public:
+ Bone() : ColorMap() {
+ init(256);
+ }
+
+ Bone(int n) : ColorMap() {
+ init(n);
+ }
+
+ void init(int n) {
+ float r[] = { 0, 0.01388888888888889, 0.02777777777777778, 0.04166666666666666, 0.05555555555555555, 0.06944444444444445, 0.08333333333333333, 0.09722222222222221, 0.1111111111111111, 0.125, 0.1388888888888889, 0.1527777777777778, 0.1666666666666667, 0.1805555555555556, 0.1944444444444444, 0.2083333333333333, 0.2222222222222222, 0.2361111111111111, 0.25, 0.2638888888888889, 0.2777777777777778, 0.2916666666666666, 0.3055555555555555, 0.3194444444444444, 0.3333333333333333, 0.3472222222222222, 0.3611111111111111, 0.375, 0.3888888888888888, 0.4027777777777777, 0.4166666666666666, 0.4305555555555555, 0.4444444444444444, 0.4583333333333333, 0.4722222222222222, 0.4861111111111112, 0.5, 0.5138888888888888, 0.5277777777777778, 0.5416666666666667, 0.5555555555555556, 0.5694444444444444, 0.5833333333333333, 0.5972222222222222, 0.611111111111111, 0.6249999999999999, 0.6388888888888888, 0.6527777777777778, 0.6726190476190474, 0.6944444444444442, 0.7162698412698412, 0.7380952380952381, 0.7599206349206349, 0.7817460317460316, 0.8035714285714286, 0.8253968253968254, 0.8472222222222221, 0.8690476190476188, 0.8908730158730158, 0.9126984126984128, 0.9345238095238095, 0.9563492063492063, 0.978174603174603, 1};
+ float g[] = { 0, 0.01388888888888889, 0.02777777777777778, 0.04166666666666666, 0.05555555555555555, 0.06944444444444445, 0.08333333333333333, 0.09722222222222221, 0.1111111111111111, 0.125, 0.1388888888888889, 0.1527777777777778, 0.1666666666666667, 0.1805555555555556, 0.1944444444444444, 0.2083333333333333, 0.2222222222222222, 0.2361111111111111, 0.25, 0.2638888888888889, 0.2777777777777778, 0.2916666666666666, 0.3055555555555555, 0.3194444444444444, 0.3353174603174602, 0.3544973544973544, 0.3736772486772486, 0.3928571428571428, 0.412037037037037, 0.4312169312169312, 0.4503968253968254, 0.4695767195767195, 0.4887566137566137, 0.5079365079365078, 0.5271164021164021, 0.5462962962962963, 0.5654761904761904, 0.5846560846560845, 0.6038359788359787, 0.623015873015873, 0.6421957671957671, 0.6613756613756612, 0.6805555555555555, 0.6997354497354497, 0.7189153439153438, 0.7380952380952379, 0.7572751322751322, 0.7764550264550264, 0.7916666666666666, 0.8055555555555555, 0.8194444444444444, 0.8333333333333334, 0.8472222222222222, 0.861111111111111, 0.875, 0.8888888888888888, 0.9027777777777777, 0.9166666666666665, 0.9305555555555555, 0.9444444444444444, 0.9583333333333333, 0.9722222222222221, 0.986111111111111, 1};
+ float b[] = { 0, 0.01917989417989418, 0.03835978835978836, 0.05753968253968253, 0.07671957671957672, 0.09589947089947089, 0.1150793650793651, 0.1342592592592592, 0.1534391534391534, 0.1726190476190476, 0.1917989417989418, 0.210978835978836, 0.2301587301587301, 0.2493386243386243, 0.2685185185185185, 0.2876984126984127, 0.3068783068783069, 0.326058201058201, 0.3452380952380952, 0.3644179894179894, 0.3835978835978835, 0.4027777777777777, 0.4219576719576719, 0.4411375661375661, 0.4583333333333333, 0.4722222222222222, 0.4861111111111111, 0.5, 0.5138888888888888, 0.5277777777777777, 0.5416666666666666, 0.5555555555555556, 0.5694444444444444, 0.5833333333333333, 0.5972222222222222, 0.6111111111111112, 0.625, 0.6388888888888888, 0.6527777777777778, 0.6666666666666667, 0.6805555555555556, 0.6944444444444444, 0.7083333333333333, 0.7222222222222222, 0.736111111111111, 0.7499999999999999, 0.7638888888888888, 0.7777777777777778, 0.7916666666666666, 0.8055555555555555, 0.8194444444444444, 0.8333333333333334, 0.8472222222222222, 0.861111111111111, 0.875, 0.8888888888888888, 0.9027777777777777, 0.9166666666666665, 0.9305555555555555, 0.9444444444444444, 0.9583333333333333, 0.9722222222222221, 0.986111111111111, 1};
+ Mat X = linspace(0,1,64);
+ this->_lut = ColorMap::linear_colormap(X,
+ Mat(64,1, CV_32FC1, r).clone(), // red
+ Mat(64,1, CV_32FC1, g).clone(), // green
+ Mat(64,1, CV_32FC1, b).clone(), // blue
+ n); // number of sample points
+ }
+ };
+
+
+
+
+ // Equals the GNU Octave colormap "jet".
+ class Jet : public ColorMap {
+
+ public:
+ Jet() {
+ init(256);
+ }
+ Jet(int n) : ColorMap() {
+ init(n);
+ }
+
+ void init(int n) {
+ // breakpoints
+ Mat X = linspace(0,1,256);
+ // define the basemap
+ float r[] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00588235294117645,0.02156862745098032,0.03725490196078418,0.05294117647058827,0.06862745098039214,0.084313725490196,0.1000000000000001,0.115686274509804,0.1313725490196078,0.1470588235294117,0.1627450980392156,0.1784313725490196,0.1941176470588235,0.2098039215686274,0.2254901960784315,0.2411764705882353,0.2568627450980392,0.2725490196078431,0.2882352941176469,0.303921568627451,0.3196078431372549,0.3352941176470587,0.3509803921568628,0.3666666666666667,0.3823529411764706,0.3980392156862744,0.4137254901960783,0.4294117647058824,0.4450980392156862,0.4607843137254901,0.4764705882352942,0.4921568627450981,0.5078431372549019,0.5235294117647058,0.5392156862745097,0.5549019607843135,0.5705882352941174,0.5862745098039217,0.6019607843137256,0.6176470588235294,0.6333333333333333,0.6490196078431372,0.664705882352941,0.6803921568627449,0.6960784313725492,0.7117647058823531,0.7274509803921569,0.7431372549019608,0.7588235294117647,0.7745098039215685,0.7901960784313724,0.8058823529411763,0.8215686274509801,0.8372549019607844,0.8529411764705883,0.8686274509803922,0.884313725490196,0.8999999999999999,0.9156862745098038,0.9313725490196076,0.947058823529412,0.9627450980392158,0.9784313725490197,0.9941176470588236,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.9862745098039216,0.9705882352941178,0.9549019607843139,0.93921568627451,0.9235294117647062,0.9078431372549018,0.892156862745098,0.8764705882352941,0.8607843137254902,0.8450980392156864,0.8294117647058825,0.8137254901960786,0.7980392156862743,0.7823529411764705,0.7666666666666666,0.7509803921568627,0.7352941176470589,0.719607843137255,0.7039215686274511,0.6882352941176473,0.6725490196078434,0.6568627450980391,0.6411764705882352,0.6254901960784314,0.6098039215686275,0.5941176470588236,0.5784313725490198,0.5627450980392159,0.5470588235294116,0.5313725490196077,0.5156862745098039,0.5};
+ float g[] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.001960784313725483,0.01764705882352935,0.03333333333333333,0.0490196078431373,0.06470588235294117,0.08039215686274503,0.09607843137254901,0.111764705882353,0.1274509803921569,0.1431372549019607,0.1588235294117647,0.1745098039215687,0.1901960784313725,0.2058823529411764,0.2215686274509804,0.2372549019607844,0.2529411764705882,0.2686274509803921,0.2843137254901961,0.3,0.3156862745098039,0.3313725490196078,0.3470588235294118,0.3627450980392157,0.3784313725490196,0.3941176470588235,0.4098039215686274,0.4254901960784314,0.4411764705882353,0.4568627450980391,0.4725490196078431,0.4882352941176471,0.503921568627451,0.5196078431372548,0.5352941176470587,0.5509803921568628,0.5666666666666667,0.5823529411764705,0.5980392156862746,0.6137254901960785,0.6294117647058823,0.6450980392156862,0.6607843137254901,0.6764705882352942,0.692156862745098,0.7078431372549019,0.723529411764706,0.7392156862745098,0.7549019607843137,0.7705882352941176,0.7862745098039214,0.8019607843137255,0.8176470588235294,0.8333333333333333,0.8490196078431373,0.8647058823529412,0.8803921568627451,0.8960784313725489,0.9117647058823528,0.9274509803921569,0.9431372549019608,0.9588235294117646,0.9745098039215687,0.9901960784313726,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.9901960784313726,0.9745098039215687,0.9588235294117649,0.943137254901961,0.9274509803921571,0.9117647058823528,0.8960784313725489,0.8803921568627451,0.8647058823529412,0.8490196078431373,0.8333333333333335,0.8176470588235296,0.8019607843137253,0.7862745098039214,0.7705882352941176,0.7549019607843137,0.7392156862745098,0.723529411764706,0.7078431372549021,0.6921568627450982,0.6764705882352944,0.6607843137254901,0.6450980392156862,0.6294117647058823,0.6137254901960785,0.5980392156862746,0.5823529411764707,0.5666666666666669,0.5509803921568626,0.5352941176470587,0.5196078431372548,0.503921568627451,0.4882352941176471,0.4725490196078432,0.4568627450980394,0.4411764705882355,0.4254901960784316,0.4098039215686273,0.3941176470588235,0.3784313725490196,0.3627450980392157,0.3470588235294119,0.331372549019608,0.3156862745098041,0.2999999999999998,0.284313725490196,0.2686274509803921,0.2529411764705882,0.2372549019607844,0.2215686274509805,0.2058823529411766,0.1901960784313728,0.1745098039215689,0.1588235294117646,0.1431372549019607,0.1274509803921569,0.111764705882353,0.09607843137254912,0.08039215686274526,0.06470588235294139,0.04901960784313708,0.03333333333333321,0.01764705882352935,0.001960784313725483,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+ float b[] = {0.5,0.5156862745098039,0.5313725490196078,0.5470588235294118,0.5627450980392157,0.5784313725490196,0.5941176470588235,0.6098039215686275,0.6254901960784314,0.6411764705882352,0.6568627450980392,0.6725490196078432,0.6882352941176471,0.7039215686274509,0.7196078431372549,0.7352941176470589,0.7509803921568627,0.7666666666666666,0.7823529411764706,0.7980392156862746,0.8137254901960784,0.8294117647058823,0.8450980392156863,0.8607843137254902,0.8764705882352941,0.892156862745098,0.907843137254902,0.9235294117647059,0.9392156862745098,0.9549019607843137,0.9705882352941176,0.9862745098039216,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.9941176470588236,0.9784313725490197,0.9627450980392158,0.9470588235294117,0.9313725490196079,0.915686274509804,0.8999999999999999,0.884313725490196,0.8686274509803922,0.8529411764705883,0.8372549019607844,0.8215686274509804,0.8058823529411765,0.7901960784313726,0.7745098039215685,0.7588235294117647,0.7431372549019608,0.7274509803921569,0.7117647058823531,0.696078431372549,0.6803921568627451,0.6647058823529413,0.6490196078431372,0.6333333333333333,0.6176470588235294,0.6019607843137256,0.5862745098039217,0.5705882352941176,0.5549019607843138,0.5392156862745099,0.5235294117647058,0.5078431372549019,0.4921568627450981,0.4764705882352942,0.4607843137254903,0.4450980392156865,0.4294117647058826,0.4137254901960783,0.3980392156862744,0.3823529411764706,0.3666666666666667,0.3509803921568628,0.335294117647059,0.3196078431372551,0.3039215686274508,0.2882352941176469,0.2725490196078431,0.2568627450980392,0.2411764705882353,0.2254901960784315,0.2098039215686276,0.1941176470588237,0.1784313725490199,0.1627450980392156,0.1470588235294117,0.1313725490196078,0.115686274509804,0.1000000000000001,0.08431372549019622,0.06862745098039236,0.05294117647058805,0.03725490196078418,0.02156862745098032,0.00588235294117645,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+ // now build lookup table
+ this->_lut = ColorMap::linear_colormap(X,
+ Mat(256,1, CV_32FC1, r).clone(), // red
+ Mat(256,1, CV_32FC1, g).clone(), // green
+ Mat(256,1, CV_32FC1, b).clone(), // blue
+ n);
+ }
+ };
+
+ // Equals the GNU Octave colormap "winter".
+ class Winter : public ColorMap {
+ public:
+ Winter() : ColorMap() {
+ init(256);
+ }
+
+ Winter(int n) : ColorMap() {
+ init(n);
+ }
+
+ void init(int n) {
+ float r[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+ float g[] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
+ float b[] = {1.0, 0.95, 0.9, 0.85, 0.8, 0.75, 0.7, 0.65, 0.6, 0.55, 0.5};
+ Mat X = linspace(0,1,11);
+ this->_lut = ColorMap::linear_colormap(X,
+ Mat(11,1, CV_32FC1, r).clone(), // red
+ Mat(11,1, CV_32FC1, g).clone(), // green
+ Mat(11,1, CV_32FC1, b).clone(), // blue
+ n); // number of sample points
+ }
+ };
+
+ // Equals the GNU Octave colormap "rainbow".
+ class Rainbow : public ColorMap {
+ public:
+ Rainbow() : ColorMap() {
+ init(256);
+ }
+
+ Rainbow(int n) : ColorMap() {
+ init(n);
+ }
+
+ void init(int n) {
+ float r[] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0.9365079365079367, 0.8571428571428572, 0.7777777777777777, 0.6984126984126986, 0.6190476190476191, 0.53968253968254, 0.4603174603174605, 0.3809523809523814, 0.3015873015873018, 0.2222222222222223, 0.1428571428571432, 0.06349206349206415, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.03174603174603208, 0.08465608465608465, 0.1375661375661377, 0.1904761904761907, 0.2433862433862437, 0.2962962962962963, 0.3492063492063493, 0.4021164021164023, 0.4550264550264553, 0.5079365079365079, 0.5608465608465609, 0.6137566137566139, 0.666666666666667};
+ float g[] = { 0, 0.03968253968253968, 0.07936507936507936, 0.119047619047619, 0.1587301587301587, 0.1984126984126984, 0.2380952380952381, 0.2777777777777778, 0.3174603174603174, 0.3571428571428571, 0.3968253968253968, 0.4365079365079365, 0.4761904761904762, 0.5158730158730158, 0.5555555555555556, 0.5952380952380952, 0.6349206349206349, 0.6746031746031745, 0.7142857142857142, 0.753968253968254, 0.7936507936507936, 0.8333333333333333, 0.873015873015873, 0.9126984126984127, 0.9523809523809523, 0.992063492063492, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0.9841269841269842, 0.9047619047619047, 0.8253968253968256, 0.7460317460317465, 0.666666666666667, 0.587301587301587, 0.5079365079365079, 0.4285714285714288, 0.3492063492063493, 0.2698412698412698, 0.1904761904761907, 0.1111111111111116, 0.03174603174603208, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+ float b[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.01587301587301582, 0.09523809523809534, 0.1746031746031744, 0.2539682539682535, 0.333333333333333, 0.412698412698413, 0.4920634920634921, 0.5714285714285712, 0.6507936507936507, 0.7301587301587302, 0.8095238095238093, 0.8888888888888884, 0.9682539682539679, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
+ Mat X = linspace(0,1,64);
+ this->_lut = ColorMap::linear_colormap(X,
+ Mat(64,1, CV_32FC1, r).clone(), // red
+ Mat(64,1, CV_32FC1, g).clone(), // green
+ Mat(64,1, CV_32FC1, b).clone(), // blue
+ n); // number of sample points
+ }
+ };
+
+ // Equals the GNU Octave colormap "ocean".
+ class Ocean : public ColorMap {
+ public:
+ Ocean() : ColorMap() {
+ init(256);
+ }
+
+ Ocean(int n) : ColorMap() {
+ init(n);
+ }
+
+ void init(int n) {
+ float r[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.04761904761904762, 0.09523809523809523, 0.1428571428571428, 0.1904761904761905, 0.2380952380952381, 0.2857142857142857, 0.3333333333333333, 0.3809523809523809, 0.4285714285714285, 0.4761904761904762, 0.5238095238095238, 0.5714285714285714, 0.6190476190476191, 0.6666666666666666, 0.7142857142857143, 0.7619047619047619, 0.8095238095238095, 0.8571428571428571, 0.9047619047619048, 0.9523809523809523, 1};
+ float g[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.02380952380952381, 0.04761904761904762, 0.07142857142857142, 0.09523809523809523, 0.119047619047619, 0.1428571428571428, 0.1666666666666667, 0.1904761904761905, 0.2142857142857143, 0.2380952380952381, 0.2619047619047619, 0.2857142857142857, 0.3095238095238095, 0.3333333333333333, 0.3571428571428572, 0.3809523809523809, 0.4047619047619048, 0.4285714285714285, 0.4523809523809524, 0.4761904761904762, 0.5, 0.5238095238095238, 0.5476190476190477, 0.5714285714285714, 0.5952380952380952, 0.6190476190476191, 0.6428571428571429, 0.6666666666666666, 0.6904761904761905, 0.7142857142857143, 0.7380952380952381, 0.7619047619047619, 0.7857142857142857, 0.8095238095238095, 0.8333333333333334, 0.8571428571428571, 0.8809523809523809, 0.9047619047619048, 0.9285714285714286, 0.9523809523809523, 0.9761904761904762, 1};
+ float b[] = { 0, 0.01587301587301587, 0.03174603174603174, 0.04761904761904762, 0.06349206349206349, 0.07936507936507936, 0.09523809523809523, 0.1111111111111111, 0.126984126984127, 0.1428571428571428, 0.1587301587301587, 0.1746031746031746, 0.1904761904761905, 0.2063492063492063, 0.2222222222222222, 0.2380952380952381, 0.253968253968254, 0.2698412698412698, 0.2857142857142857, 0.3015873015873016, 0.3174603174603174, 0.3333333333333333, 0.3492063492063492, 0.3650793650793651, 0.3809523809523809, 0.3968253968253968, 0.4126984126984127, 0.4285714285714285, 0.4444444444444444, 0.4603174603174603, 0.4761904761904762, 0.492063492063492, 0.5079365079365079, 0.5238095238095238, 0.5396825396825397, 0.5555555555555556, 0.5714285714285714, 0.5873015873015873, 0.6031746031746031, 0.6190476190476191, 0.6349206349206349, 0.6507936507936508, 0.6666666666666666, 0.6825396825396826, 0.6984126984126984, 0.7142857142857143, 0.7301587301587301, 0.746031746031746, 0.7619047619047619, 0.7777777777777778, 0.7936507936507936, 0.8095238095238095, 0.8253968253968254, 0.8412698412698413, 0.8571428571428571, 0.873015873015873, 0.8888888888888888, 0.9047619047619048, 0.9206349206349206, 0.9365079365079365, 0.9523809523809523, 0.9682539682539683, 0.9841269841269841, 1};
+ Mat X = linspace(0,1,64);
+ this->_lut = ColorMap::linear_colormap(X,
+ Mat(64,1, CV_32FC1, r).clone(), // red
+ Mat(64,1, CV_32FC1, g).clone(), // green
+ Mat(64,1, CV_32FC1, b).clone(), // blue
+ n); // number of sample points
+ }
+ };
+
+ // Equals the GNU Octave colormap "summer".
+ class Summer : public ColorMap {
+ public:
+ Summer() : ColorMap() {
+ init(256);
+ }
+
+ Summer(int n) : ColorMap() {
+ init(n);
+ }
+
+ void init(int n) {
+ float r[] = { 0, 0.01587301587301587, 0.03174603174603174, 0.04761904761904762, 0.06349206349206349, 0.07936507936507936, 0.09523809523809523, 0.1111111111111111, 0.126984126984127, 0.1428571428571428, 0.1587301587301587, 0.1746031746031746, 0.1904761904761905, 0.2063492063492063, 0.2222222222222222, 0.2380952380952381, 0.253968253968254, 0.2698412698412698, 0.2857142857142857, 0.3015873015873016, 0.3174603174603174, 0.3333333333333333, 0.3492063492063492, 0.3650793650793651, 0.3809523809523809, 0.3968253968253968, 0.4126984126984127, 0.4285714285714285, 0.4444444444444444, 0.4603174603174603, 0.4761904761904762, 0.492063492063492, 0.5079365079365079, 0.5238095238095238, 0.5396825396825397, 0.5555555555555556, 0.5714285714285714, 0.5873015873015873, 0.6031746031746031, 0.6190476190476191, 0.6349206349206349, 0.6507936507936508, 0.6666666666666666, 0.6825396825396826, 0.6984126984126984, 0.7142857142857143, 0.7301587301587301, 0.746031746031746, 0.7619047619047619, 0.7777777777777778, 0.7936507936507936, 0.8095238095238095, 0.8253968253968254, 0.8412698412698413, 0.8571428571428571, 0.873015873015873, 0.8888888888888888, 0.9047619047619048, 0.9206349206349206, 0.9365079365079365, 0.9523809523809523, 0.9682539682539683, 0.9841269841269841, 1};
+ float g[] = { 0.5, 0.5079365079365079, 0.5158730158730158, 0.5238095238095238, 0.5317460317460317, 0.5396825396825397, 0.5476190476190477, 0.5555555555555556, 0.5634920634920635, 0.5714285714285714, 0.5793650793650793, 0.5873015873015873, 0.5952380952380952, 0.6031746031746031, 0.6111111111111112, 0.6190476190476191, 0.626984126984127, 0.6349206349206349, 0.6428571428571428, 0.6507936507936508, 0.6587301587301587, 0.6666666666666666, 0.6746031746031746, 0.6825396825396826, 0.6904761904761905, 0.6984126984126984, 0.7063492063492063, 0.7142857142857143, 0.7222222222222222, 0.7301587301587301, 0.7380952380952381, 0.746031746031746, 0.753968253968254, 0.7619047619047619, 0.7698412698412698, 0.7777777777777778, 0.7857142857142857, 0.7936507936507937, 0.8015873015873016, 0.8095238095238095, 0.8174603174603174, 0.8253968253968254, 0.8333333333333333, 0.8412698412698413, 0.8492063492063492, 0.8571428571428572, 0.8650793650793651, 0.873015873015873, 0.8809523809523809, 0.8888888888888888, 0.8968253968253967, 0.9047619047619048, 0.9126984126984127, 0.9206349206349207, 0.9285714285714286, 0.9365079365079365, 0.9444444444444444, 0.9523809523809523, 0.9603174603174602, 0.9682539682539683, 0.9761904761904762, 0.9841269841269842, 0.9920634920634921, 1};
+ float b[] = { 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4};
+ Mat X = linspace(0,1,64);
+ this->_lut = ColorMap::linear_colormap(X,
+ Mat(64,1, CV_32FC1, r).clone(), // red
+ Mat(64,1, CV_32FC1, g).clone(), // green
+ Mat(64,1, CV_32FC1, b).clone(), // blue
+ n); // number of sample points
+ }
+ };
+
+ // Equals the GNU Octave colormap "spring".
+ class Spring : public ColorMap {
+ public:
+ Spring() : ColorMap() {
+ init(256);
+ }
+
+ Spring(int n) : ColorMap() {
+ init(n);
+ }
+
+ void init(int n) {
+ float r[] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
+ float g[] = { 0, 0.01587301587301587, 0.03174603174603174, 0.04761904761904762, 0.06349206349206349, 0.07936507936507936, 0.09523809523809523, 0.1111111111111111, 0.126984126984127, 0.1428571428571428, 0.1587301587301587, 0.1746031746031746, 0.1904761904761905, 0.2063492063492063, 0.2222222222222222, 0.2380952380952381, 0.253968253968254, 0.2698412698412698, 0.2857142857142857, 0.3015873015873016, 0.3174603174603174, 0.3333333333333333, 0.3492063492063492, 0.3650793650793651, 0.3809523809523809, 0.3968253968253968, 0.4126984126984127, 0.4285714285714285, 0.4444444444444444, 0.4603174603174603, 0.4761904761904762, 0.492063492063492, 0.5079365079365079, 0.5238095238095238, 0.5396825396825397, 0.5555555555555556, 0.5714285714285714, 0.5873015873015873, 0.6031746031746031, 0.6190476190476191, 0.6349206349206349, 0.6507936507936508, 0.6666666666666666, 0.6825396825396826, 0.6984126984126984, 0.7142857142857143, 0.7301587301587301, 0.746031746031746, 0.7619047619047619, 0.7777777777777778, 0.7936507936507936, 0.8095238095238095, 0.8253968253968254, 0.8412698412698413, 0.8571428571428571, 0.873015873015873, 0.8888888888888888, 0.9047619047619048, 0.9206349206349206, 0.9365079365079365, 0.9523809523809523, 0.9682539682539683, 0.9841269841269841, 1};
+ float b[] = { 1, 0.9841269841269842, 0.9682539682539683, 0.9523809523809523, 0.9365079365079365, 0.9206349206349207, 0.9047619047619048, 0.8888888888888888, 0.873015873015873, 0.8571428571428572, 0.8412698412698413, 0.8253968253968254, 0.8095238095238095, 0.7936507936507937, 0.7777777777777778, 0.7619047619047619, 0.746031746031746, 0.7301587301587302, 0.7142857142857143, 0.6984126984126984, 0.6825396825396826, 0.6666666666666667, 0.6507936507936508, 0.6349206349206349, 0.6190476190476191, 0.6031746031746033, 0.5873015873015873, 0.5714285714285714, 0.5555555555555556, 0.5396825396825398, 0.5238095238095238, 0.5079365079365079, 0.4920634920634921, 0.4761904761904762, 0.4603174603174603, 0.4444444444444444, 0.4285714285714286, 0.4126984126984127, 0.3968253968253969, 0.3809523809523809, 0.3650793650793651, 0.3492063492063492, 0.3333333333333334, 0.3174603174603174, 0.3015873015873016, 0.2857142857142857, 0.2698412698412699, 0.253968253968254, 0.2380952380952381, 0.2222222222222222, 0.2063492063492064, 0.1904761904761905, 0.1746031746031746, 0.1587301587301587, 0.1428571428571429, 0.126984126984127, 0.1111111111111112, 0.09523809523809523, 0.07936507936507942, 0.06349206349206349, 0.04761904761904767, 0.03174603174603174, 0.01587301587301593, 0};
+ Mat X = linspace(0,1,64);
+ this->_lut = ColorMap::linear_colormap(X,
+ Mat(64,1, CV_32FC1, r).clone(), // red
+ Mat(64,1, CV_32FC1, g).clone(), // green
+ Mat(64,1, CV_32FC1, b).clone(), // blue
+ n); // number of sample points
+ }
+ };
+
+ // Equals the GNU Octave colormap "cool".
+ class Cool : public ColorMap {
+ public:
+ Cool() : ColorMap() {
+ init(256);
+ }
+
+ Cool(int n) : ColorMap() {
+ init(n);
+ }
+
+ void init(int n) {
+ float r[] = { 0, 0.01587301587301587, 0.03174603174603174, 0.04761904761904762, 0.06349206349206349, 0.07936507936507936, 0.09523809523809523, 0.1111111111111111, 0.126984126984127, 0.1428571428571428, 0.1587301587301587, 0.1746031746031746, 0.1904761904761905, 0.2063492063492063, 0.2222222222222222, 0.2380952380952381, 0.253968253968254, 0.2698412698412698, 0.2857142857142857, 0.3015873015873016, 0.3174603174603174, 0.3333333333333333, 0.3492063492063492, 0.3650793650793651, 0.3809523809523809, 0.3968253968253968, 0.4126984126984127, 0.4285714285714285, 0.4444444444444444, 0.4603174603174603, 0.4761904761904762, 0.492063492063492, 0.5079365079365079, 0.5238095238095238, 0.5396825396825397, 0.5555555555555556, 0.5714285714285714, 0.5873015873015873, 0.6031746031746031, 0.6190476190476191, 0.6349206349206349, 0.6507936507936508, 0.6666666666666666, 0.6825396825396826, 0.6984126984126984, 0.7142857142857143, 0.7301587301587301, 0.746031746031746, 0.7619047619047619, 0.7777777777777778, 0.7936507936507936, 0.8095238095238095, 0.8253968253968254, 0.8412698412698413, 0.8571428571428571, 0.873015873015873, 0.8888888888888888, 0.9047619047619048, 0.9206349206349206, 0.9365079365079365, 0.9523809523809523, 0.9682539682539683, 0.9841269841269841, 1};
+ float g[] = { 1, 0.9841269841269842, 0.9682539682539683, 0.9523809523809523, 0.9365079365079365, 0.9206349206349207, 0.9047619047619048, 0.8888888888888888, 0.873015873015873, 0.8571428571428572, 0.8412698412698413, 0.8253968253968254, 0.8095238095238095, 0.7936507936507937, 0.7777777777777778, 0.7619047619047619, 0.746031746031746, 0.7301587301587302, 0.7142857142857143, 0.6984126984126984, 0.6825396825396826, 0.6666666666666667, 0.6507936507936508, 0.6349206349206349, 0.6190476190476191, 0.6031746031746033, 0.5873015873015873, 0.5714285714285714, 0.5555555555555556, 0.5396825396825398, 0.5238095238095238, 0.5079365079365079, 0.4920634920634921, 0.4761904761904762, 0.4603174603174603, 0.4444444444444444, 0.4285714285714286, 0.4126984126984127, 0.3968253968253969, 0.3809523809523809, 0.3650793650793651, 0.3492063492063492, 0.3333333333333334, 0.3174603174603174, 0.3015873015873016, 0.2857142857142857, 0.2698412698412699, 0.253968253968254, 0.2380952380952381, 0.2222222222222222, 0.2063492063492064, 0.1904761904761905, 0.1746031746031746, 0.1587301587301587, 0.1428571428571429, 0.126984126984127, 0.1111111111111112, 0.09523809523809523, 0.07936507936507942, 0.06349206349206349, 0.04761904761904767, 0.03174603174603174, 0.01587301587301593, 0};
+ float b[] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
+ Mat X = linspace(0,1,64);
+ this->_lut = ColorMap::linear_colormap(X,
+ Mat(64,1, CV_32FC1, r).clone(), // red
+ Mat(64,1, CV_32FC1, g).clone(), // green
+ Mat(64,1, CV_32FC1, b).clone(), // blue
+ n); // number of sample points
+ }
+ };
+
+ // Equals the GNU Octave colormap "hsv".
+ class HSV : public ColorMap {
+ public:
+ HSV() : ColorMap() {
+ init(256);
+ }
+
+ HSV(int n) : ColorMap() {
+ init(n);
+ }
+
+ void init(int n) {
+ float r[] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0.9523809523809526, 0.8571428571428568, 0.7619047619047614, 0.6666666666666665, 0.5714285714285716, 0.4761904761904763, 0.3809523809523805, 0.2857142857142856, 0.1904761904761907, 0.0952380952380949, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.09523809523809557, 0.1904761904761905, 0.2857142857142854, 0.3809523809523809, 0.4761904761904765, 0.5714285714285714, 0.6666666666666663, 0.7619047619047619, 0.8571428571428574, 0.9523809523809523, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
+ float g[] = { 0, 0.09523809523809523, 0.1904761904761905, 0.2857142857142857, 0.3809523809523809, 0.4761904761904762, 0.5714285714285714, 0.6666666666666666, 0.7619047619047619, 0.8571428571428571, 0.9523809523809523, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0.9523809523809526, 0.8571428571428577, 0.7619047619047619, 0.6666666666666665, 0.5714285714285716, 0.4761904761904767, 0.3809523809523814, 0.2857142857142856, 0.1904761904761907, 0.09523809523809579, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+ float b[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.09523809523809523, 0.1904761904761905, 0.2857142857142857, 0.3809523809523809, 0.4761904761904762, 0.5714285714285714, 0.6666666666666666, 0.7619047619047619, 0.8571428571428571, 0.9523809523809523, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0.9523809523809526, 0.8571428571428577, 0.7619047619047614, 0.6666666666666665, 0.5714285714285716, 0.4761904761904767, 0.3809523809523805, 0.2857142857142856, 0.1904761904761907, 0.09523809523809579, 0};
+ Mat X = linspace(0,1,64);
+ this->_lut = ColorMap::linear_colormap(X,
+ Mat(64,1, CV_32FC1, r).clone(), // red
+ Mat(64,1, CV_32FC1, g).clone(), // green
+ Mat(64,1, CV_32FC1, b).clone(), // blue
+ n); // number of sample points
+ }
+ };
+
+ // Equals the GNU Octave colormap "pink".
+ class Pink : public ColorMap {
+ public:
+ Pink() : ColorMap() {
+ init(256);
+ }
+
+ Pink(int n) : ColorMap() {
+ init(n);
+ }
+
+ void init(int n) {
+ float r[] = { 0, 0.1571348402636772, 0.2222222222222222, 0.2721655269759087, 0.3142696805273544, 0.3513641844631533, 0.3849001794597505, 0.415739709641549, 0.4444444444444444, 0.4714045207910317, 0.4969039949999532, 0.5211573066470477, 0.5443310539518174, 0.5665577237325317, 0.5879447357921312, 0.6085806194501846, 0.6285393610547089, 0.6478835438717, 0.6666666666666666, 0.6849348892187751, 0.7027283689263065, 0.7200822998230956, 0.7370277311900888, 0.753592220347252, 0.7663560447348133, 0.7732293307186413, 0.7800420555749596, 0.7867957924694432, 0.7934920476158722, 0.8001322641986387, 0.8067178260046388, 0.8132500607904444, 0.8197302434079591, 0.8261595987094034, 0.8325393042503717, 0.8388704928078611, 0.8451542547285166, 0.8513916401208816, 0.8575836609041332, 0.8637312927246217, 0.8698354767504924, 0.8758971213537393, 0.8819171036881968, 0.8878962711712378, 0.8938354428762595, 0.8997354108424372, 0.9055969413076769, 0.9114207758701963, 0.9172076325837248, 0.9229582069908971, 0.9286731730990523, 0.9343531843023135, 0.9399988742535192, 0.9456108576893002, 0.9511897312113418, 0.9567360740266436, 0.9622504486493763, 0.9677334015667416, 0.9731854638710686, 0.9786071518602129, 0.9839989676081821, 0.9893613995077727, 0.9946949227868761, 1};
+ float g[] = { 0, 0.1028688999747279, 0.1454785934906616, 0.1781741612749496, 0.2057377999494559, 0.2300218531141181, 0.2519763153394848, 0.2721655269759087, 0.2909571869813232, 0.3086066999241838, 0.3253000243161777, 0.3411775438127727, 0.3563483225498992, 0.3708990935094579, 0.3849001794597505, 0.3984095364447979, 0.4114755998989117, 0.4241393401869012, 0.4364357804719847, 0.4483951394230328, 0.4600437062282361, 0.4714045207910317, 0.4824979096371639, 0.4933419132673033, 0.5091750772173156, 0.5328701692569688, 0.5555555555555556, 0.5773502691896257, 0.5983516452371671, 0.6186404847588913, 0.6382847385042254, 0.6573421981221795, 0.6758625033664688, 0.6938886664887108, 0.7114582486036499, 0.7286042804780002, 0.7453559924999299, 0.7617394000445604, 0.7777777777777778, 0.7934920476158723, 0.8089010988089465, 0.8240220541217402, 0.8388704928078611, 0.8534606386520677, 0.8678055195451838, 0.8819171036881968, 0.8958064164776166, 0.9094836413191612, 0.9172076325837248, 0.9229582069908971, 0.9286731730990523, 0.9343531843023135, 0.9399988742535192, 0.9456108576893002, 0.9511897312113418, 0.9567360740266436, 0.9622504486493763, 0.9677334015667416, 0.9731854638710686, 0.9786071518602129, 0.9839989676081821, 0.9893613995077727, 0.9946949227868761, 1};
+ float b[] = { 0, 0.1028688999747279, 0.1454785934906616, 0.1781741612749496, 0.2057377999494559, 0.2300218531141181, 0.2519763153394848, 0.2721655269759087, 0.2909571869813232, 0.3086066999241838, 0.3253000243161777, 0.3411775438127727, 0.3563483225498992, 0.3708990935094579, 0.3849001794597505, 0.3984095364447979, 0.4114755998989117, 0.4241393401869012, 0.4364357804719847, 0.4483951394230328, 0.4600437062282361, 0.4714045207910317, 0.4824979096371639, 0.4933419132673033, 0.5039526306789697, 0.5143444998736397, 0.5245305283129621, 0.5345224838248488, 0.5443310539518174, 0.5539659798925444, 0.563436169819011, 0.5727497953228163, 0.5819143739626463, 0.5909368402852788, 0.5998236072282915, 0.6085806194501846, 0.6172133998483676, 0.6257270902992705, 0.6341264874742278, 0.642416074439621, 0.6506000486323554, 0.6586823467062358, 0.6666666666666666, 0.6745564876468501, 0.6823550876255453, 0.6900655593423541, 0.6976908246297114, 0.7052336473499384, 0.7237468644557459, 0.7453559924999298, 0.7663560447348133, 0.7867957924694432, 0.8067178260046388, 0.8261595987094034, 0.8451542547285166, 0.8637312927246217, 0.8819171036881968, 0.8997354108424372, 0.9172076325837248, 0.9343531843023135, 0.9511897312113418, 0.9677334015667416, 0.9839989676081821, 1};
+ Mat X = linspace(0,1,64);
+ this->_lut = ColorMap::linear_colormap(X,
+ Mat(64,1, CV_32FC1, r).clone(), // red
+ Mat(64,1, CV_32FC1, g).clone(), // green
+ Mat(64,1, CV_32FC1, b).clone(), // blue
+ n); // number of sample points
+ }
+ };
+
+ // Equals the GNU Octave colormap "hot".
+ class Hot : public ColorMap {
+ public:
+ Hot() : ColorMap() {
+ init(256);
+ }
+
+ Hot(int n) : ColorMap() {
+ init(n);
+ }
+
+ void init(int n) {
+ float r[] = { 0, 0.03968253968253968, 0.07936507936507936, 0.119047619047619, 0.1587301587301587, 0.1984126984126984, 0.2380952380952381, 0.2777777777777778, 0.3174603174603174, 0.3571428571428571, 0.3968253968253968, 0.4365079365079365, 0.4761904761904762, 0.5158730158730158, 0.5555555555555556, 0.5952380952380952, 0.6349206349206349, 0.6746031746031745, 0.7142857142857142, 0.753968253968254, 0.7936507936507936, 0.8333333333333333, 0.873015873015873, 0.9126984126984127, 0.9523809523809523, 0.992063492063492, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
+ float g[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.03174603174603163, 0.0714285714285714, 0.1111111111111112, 0.1507936507936507, 0.1904761904761905, 0.23015873015873, 0.2698412698412698, 0.3095238095238093, 0.3492063492063491, 0.3888888888888888, 0.4285714285714284, 0.4682539682539679, 0.5079365079365079, 0.5476190476190477, 0.5873015873015872, 0.6269841269841268, 0.6666666666666665, 0.7063492063492065, 0.746031746031746, 0.7857142857142856, 0.8253968253968254, 0.8650793650793651, 0.9047619047619047, 0.9444444444444442, 0.984126984126984, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
+ float b[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.04761904761904745, 0.1269841269841265, 0.2063492063492056, 0.2857142857142856, 0.3650793650793656, 0.4444444444444446, 0.5238095238095237, 0.6031746031746028, 0.6825396825396828, 0.7619047619047619, 0.8412698412698409, 0.92063492063492, 1};
+ Mat X = linspace(0,1,64);
+ this->_lut = ColorMap::linear_colormap(X,
+ Mat(64,1, CV_32FC1, r).clone(), // red
+ Mat(64,1, CV_32FC1, g).clone(), // green
+ Mat(64,1, CV_32FC1, b).clone(), // blue
+ n); // number of sample points
+ }
+ };
+
+ void ColorMap::operator()(InputArray _src, OutputArray _dst) const
+ {
+ if(_lut.total() != 256)
+ CV_Error(Error::StsAssert, "cv::LUT only supports tables of size 256.");
+ Mat src = _src.getMat();
+ // Return original matrix if wrong type is given (is fail loud better here?)
+ if(src.type() != CV_8UC1 && src.type() != CV_8UC3)
+ {
+ src.copyTo(_dst);
+ return;
+ }
+ // Turn into a BGR matrix into its grayscale representation.
+ if(src.type() == CV_8UC3)
+ cvtColor(src.clone(), src, COLOR_BGR2GRAY);
+ cvtColor(src.clone(), src, COLOR_GRAY2BGR);
+ // Apply the ColorMap.
+ LUT(src, _lut, _dst);
+ }
+
+ Mat ColorMap::linear_colormap(InputArray X,
+ InputArray r, InputArray g, InputArray b,
+ InputArray xi) {
+ Mat lut, lut8;
+ Mat planes[] = {
+ interp1(X, b, xi),
+ interp1(X, g, xi),
+ interp1(X, r, xi)};
+ merge(planes, 3, lut);
+ lut.convertTo(lut8, CV_8U, 255.);
+ return lut8;
+ }
+
+ }
+
+ void applyColorMap(InputArray src, OutputArray dst, int colormap)
+ {
+ colormap::ColorMap* cm =
+ colormap == COLORMAP_AUTUMN ? (colormap::ColorMap*)(new colormap::Autumn) :
+ colormap == COLORMAP_BONE ? (colormap::ColorMap*)(new colormap::Bone) :
+ colormap == COLORMAP_COOL ? (colormap::ColorMap*)(new colormap::Cool) :
+ colormap == COLORMAP_HOT ? (colormap::ColorMap*)(new colormap::Hot) :
+ colormap == COLORMAP_HSV ? (colormap::ColorMap*)(new colormap::HSV) :
+ colormap == COLORMAP_JET ? (colormap::ColorMap*)(new colormap::Jet) :
+ colormap == COLORMAP_OCEAN ? (colormap::ColorMap*)(new colormap::Ocean) :
+ colormap == COLORMAP_PINK ? (colormap::ColorMap*)(new colormap::Pink) :
+ colormap == COLORMAP_RAINBOW ? (colormap::ColorMap*)(new colormap::Rainbow) :
+ colormap == COLORMAP_SPRING ? (colormap::ColorMap*)(new colormap::Spring) :
+ colormap == COLORMAP_SUMMER ? (colormap::ColorMap*)(new colormap::Summer) :
+ colormap == COLORMAP_WINTER ? (colormap::ColorMap*)(new colormap::Winter) : 0;
+
+ if( !cm )
+ CV_Error( Error::StsBadArg, "Unknown colormap id; use one of COLORMAP_*");
+
+ (*cm)(src, dst);
+
+ delete cm;
+ }
+}