| 1 | /// @ref gtx_pca |
| 2 | |
| 3 | #include <algorithm> |
| 4 | #include <utility> |
| 5 | |
| 6 | namespace glm { |
| 7 | |
| 8 | |
| 9 | template<length_t D, typename T, qualifier Q> |
| 10 | GLM_FUNC_QUALIFIER mat<D, D, T, Q> computeCovarianceMatrix(vec<D, T, Q> const* v, size_t n) |
| 11 | { |
| 12 | return computeCovarianceMatrix<D, T, Q, vec<D, T, Q> const*>(v, v + n); |
| 13 | } |
| 14 | |
| 15 | |
| 16 | template<length_t D, typename T, qualifier Q> |
| 17 | GLM_FUNC_QUALIFIER mat<D, D, T, Q> computeCovarianceMatrix(vec<D, T, Q> const* v, size_t n, vec<D, T, Q> const& c) |
| 18 | { |
| 19 | return computeCovarianceMatrix<D, T, Q, vec<D, T, Q> const*>(v, v + n, c); |
| 20 | } |
| 21 | |
| 22 | |
| 23 | template<length_t D, typename T, qualifier Q, typename I> |
| 24 | GLM_FUNC_QUALIFIER mat<D, D, T, Q> computeCovarianceMatrix(I const& b, I const& e) |
| 25 | { |
| 26 | glm::mat<D, D, T, Q> m(0); |
| 27 | |
| 28 | size_t cnt = 0; |
| 29 | for(I i = b; i != e; i++) |
| 30 | { |
| 31 | vec<D, T, Q> const& v = *i; |
| 32 | for(length_t x = 0; x < D; ++x) |
| 33 | for(length_t y = 0; y < D; ++y) |
| 34 | m[x][y] += static_cast<T>(v[x] * v[y]); |
| 35 | cnt++; |
| 36 | } |
| 37 | if(cnt > 0) |
| 38 | m /= static_cast<T>(cnt); |
| 39 | |
| 40 | return m; |
| 41 | } |
| 42 | |
| 43 | |
| 44 | template<length_t D, typename T, qualifier Q, typename I> |
| 45 | GLM_FUNC_QUALIFIER mat<D, D, T, Q> computeCovarianceMatrix(I const& b, I const& e, vec<D, T, Q> const& c) |
| 46 | { |
| 47 | glm::mat<D, D, T, Q> m(0); |
| 48 | glm::vec<D, T, Q> v; |
| 49 | |
| 50 | size_t cnt = 0; |
| 51 | for(I i = b; i != e; i++) |
| 52 | { |
| 53 | v = *i - c; |
| 54 | for(length_t x = 0; x < D; ++x) |
| 55 | for(length_t y = 0; y < D; ++y) |
| 56 | m[x][y] += static_cast<T>(v[x] * v[y]); |
| 57 | cnt++; |
| 58 | } |
| 59 | if(cnt > 0) |
| 60 | m /= static_cast<T>(cnt); |
| 61 | |
| 62 | return m; |
| 63 | } |
| 64 | |
| 65 | namespace _internal_ |
| 66 | { |
| 67 | |
| 68 | template<typename T> |
| 69 | GLM_FUNC_QUALIFIER static T transferSign(T const& v, T const& s) |
| 70 | { |
| 71 | return ((s) >= 0 ? glm::abs(v) : -glm::abs(v)); |
| 72 | } |
| 73 | |
| 74 | template<typename T> |
| 75 | GLM_FUNC_QUALIFIER static T pythag(T const& a, T const& b) { |
| 76 | static const T epsilon = static_cast<T>(0.0000001); |
| 77 | T absa = glm::abs(a); |
| 78 | T absb = glm::abs(b); |
| 79 | if(absa > absb) { |
| 80 | absb /= absa; |
| 81 | absb *= absb; |
| 82 | return absa * glm::sqrt(static_cast<T>(1) + absb); |
| 83 | } |
| 84 | if(glm::equal<T>(absb, 0, epsilon)) return static_cast<T>(0); |
| 85 | absa /= absb; |
| 86 | absa *= absa; |
| 87 | return absb * glm::sqrt(static_cast<T>(1) + absa); |
| 88 | } |
| 89 | |
| 90 | } |
| 91 | |
| 92 | template<length_t D, typename T, qualifier Q> |
| 93 | GLM_FUNC_QUALIFIER unsigned int findEigenvaluesSymReal |
| 94 | ( |
| 95 | mat<D, D, T, Q> const& covarMat, |
| 96 | vec<D, T, Q>& outEigenvalues, |
| 97 | mat<D, D, T, Q>& outEigenvectors |
| 98 | ) |
| 99 | { |
| 100 | using _internal_::transferSign; |
| 101 | using _internal_::pythag; |
| 102 | |
| 103 | T a[D * D]; // matrix -- input and workspace for algorithm (will be changed inplace) |
| 104 | T d[D]; // diagonal elements |
| 105 | T e[D]; // off-diagonal elements |
| 106 | |
| 107 | for(length_t r = 0; r < D; r++) |
| 108 | for(length_t c = 0; c < D; c++) |
| 109 | a[(r) * D + (c)] = covarMat[c][r]; |
| 110 | |
| 111 | // 1. Householder reduction. |
| 112 | length_t l, k, j, i; |
| 113 | T scale, hh, h, g, f; |
| 114 | static const T epsilon = static_cast<T>(0.0000001); |
| 115 | |
| 116 | for(i = D; i >= 2; i--) |
| 117 | { |
| 118 | l = i - 1; |
| 119 | h = scale = 0; |
| 120 | if(l > 1) |
| 121 | { |
| 122 | for(k = 1; k <= l; k++) |
| 123 | { |
| 124 | scale += glm::abs(a[(i - 1) * D + (k - 1)]); |
| 125 | } |
| 126 | if(glm::equal<T>(scale, 0, epsilon)) |
| 127 | { |
| 128 | e[i - 1] = a[(i - 1) * D + (l - 1)]; |
| 129 | } |
| 130 | else |
| 131 | { |
| 132 | for(k = 1; k <= l; k++) |
| 133 | { |
| 134 | a[(i - 1) * D + (k - 1)] /= scale; |
| 135 | h += a[(i - 1) * D + (k - 1)] * a[(i - 1) * D + (k - 1)]; |
| 136 | } |
| 137 | f = a[(i - 1) * D + (l - 1)]; |
| 138 | g = ((f >= 0) ? -glm::sqrt(h) : glm::sqrt(h)); |
| 139 | e[i - 1] = scale * g; |
| 140 | h -= f * g; |
| 141 | a[(i - 1) * D + (l - 1)] = f - g; |
| 142 | f = 0; |
| 143 | for(j = 1; j <= l; j++) |
| 144 | { |
| 145 | a[(j - 1) * D + (i - 1)] = a[(i - 1) * D + (j - 1)] / h; |
| 146 | g = 0; |
| 147 | for(k = 1; k <= j; k++) |
| 148 | { |
| 149 | g += a[(j - 1) * D + (k - 1)] * a[(i - 1) * D + (k - 1)]; |
| 150 | } |
| 151 | for(k = j + 1; k <= l; k++) |
| 152 | { |
| 153 | g += a[(k - 1) * D + (j - 1)] * a[(i - 1) * D + (k - 1)]; |
| 154 | } |
| 155 | e[j - 1] = g / h; |
| 156 | f += e[j - 1] * a[(i - 1) * D + (j - 1)]; |
| 157 | } |
| 158 | hh = f / (h + h); |
| 159 | for(j = 1; j <= l; j++) |
| 160 | { |
| 161 | f = a[(i - 1) * D + (j - 1)]; |
| 162 | e[j - 1] = g = e[j - 1] - hh * f; |
| 163 | for(k = 1; k <= j; k++) |
| 164 | { |
| 165 | a[(j - 1) * D + (k - 1)] -= (f * e[k - 1] + g * a[(i - 1) * D + (k - 1)]); |
| 166 | } |
| 167 | } |
| 168 | } |
| 169 | } |
| 170 | else |
| 171 | { |
| 172 | e[i - 1] = a[(i - 1) * D + (l - 1)]; |
| 173 | } |
| 174 | d[i - 1] = h; |
| 175 | } |
| 176 | d[0] = 0; |
| 177 | e[0] = 0; |
| 178 | for(i = 1; i <= D; i++) |
| 179 | { |
| 180 | l = i - 1; |
| 181 | if(!glm::equal<T>(d[i - 1], 0, epsilon)) |
| 182 | { |
| 183 | for(j = 1; j <= l; j++) |
| 184 | { |
| 185 | g = 0; |
| 186 | for(k = 1; k <= l; k++) |
| 187 | { |
| 188 | g += a[(i - 1) * D + (k - 1)] * a[(k - 1) * D + (j - 1)]; |
| 189 | } |
| 190 | for(k = 1; k <= l; k++) |
| 191 | { |
| 192 | a[(k - 1) * D + (j - 1)] -= g * a[(k - 1) * D + (i - 1)]; |
| 193 | } |
| 194 | } |
| 195 | } |
| 196 | d[i - 1] = a[(i - 1) * D + (i - 1)]; |
| 197 | a[(i - 1) * D + (i - 1)] = 1; |
| 198 | for(j = 1; j <= l; j++) |
| 199 | { |
| 200 | a[(j - 1) * D + (i - 1)] = a[(i - 1) * D + (j - 1)] = 0; |
| 201 | } |
| 202 | } |
| 203 | |
| 204 | // 2. Calculation of eigenvalues and eigenvectors (QL algorithm) |
| 205 | length_t m, iter; |
| 206 | T s, r, p, dd, c, b; |
| 207 | const length_t MAX_ITER = 30; |
| 208 | |
| 209 | for(i = 2; i <= D; i++) |
| 210 | { |
| 211 | e[i - 2] = e[i - 1]; |
| 212 | } |
| 213 | e[D - 1] = 0; |
| 214 | |
| 215 | for(l = 1; l <= D; l++) |
| 216 | { |
| 217 | iter = 0; |
| 218 | do |
| 219 | { |
| 220 | for(m = l; m <= D - 1; m++) |
| 221 | { |
| 222 | dd = glm::abs(d[m - 1]) + glm::abs(d[m - 1 + 1]); |
| 223 | if(glm::equal<T>(glm::abs(e[m - 1]) + dd, dd, epsilon)) |
| 224 | break; |
| 225 | } |
| 226 | if(m != l) |
| 227 | { |
| 228 | if(iter++ == MAX_ITER) |
| 229 | { |
| 230 | return 0; // Too many iterations in FindEigenvalues |
| 231 | } |
| 232 | g = (d[l - 1 + 1] - d[l - 1]) / (2 * e[l - 1]); |
| 233 | r = pythag<T>(g, 1); |
| 234 | g = d[m - 1] - d[l - 1] + e[l - 1] / (g + transferSign(r, g)); |
| 235 | s = c = 1; |
| 236 | p = 0; |
| 237 | for(i = m - 1; i >= l; i--) |
| 238 | { |
| 239 | f = s * e[i - 1]; |
| 240 | b = c * e[i - 1]; |
| 241 | e[i - 1 + 1] = r = pythag(f, g); |
| 242 | if(glm::equal<T>(r, 0, epsilon)) |
| 243 | { |
| 244 | d[i - 1 + 1] -= p; |
| 245 | e[m - 1] = 0; |
| 246 | break; |
| 247 | } |
| 248 | s = f / r; |
| 249 | c = g / r; |
| 250 | g = d[i - 1 + 1] - p; |
| 251 | r = (d[i - 1] - g) * s + 2 * c * b; |
| 252 | d[i - 1 + 1] = g + (p = s * r); |
| 253 | g = c * r - b; |
| 254 | for(k = 1; k <= D; k++) |
| 255 | { |
| 256 | f = a[(k - 1) * D + (i - 1 + 1)]; |
| 257 | a[(k - 1) * D + (i - 1 + 1)] = s * a[(k - 1) * D + (i - 1)] + c * f; |
| 258 | a[(k - 1) * D + (i - 1)] = c * a[(k - 1) * D + (i - 1)] - s * f; |
| 259 | } |
| 260 | } |
| 261 | if(glm::equal<T>(r, 0, epsilon) && (i >= l)) |
| 262 | continue; |
| 263 | d[l - 1] -= p; |
| 264 | e[l - 1] = g; |
| 265 | e[m - 1] = 0; |
| 266 | } |
| 267 | } while(m != l); |
| 268 | } |
| 269 | |
| 270 | // 3. output |
| 271 | for(i = 0; i < D; i++) |
| 272 | outEigenvalues[i] = d[i]; |
| 273 | for(i = 0; i < D; i++) |
| 274 | for(j = 0; j < D; j++) |
| 275 | outEigenvectors[i][j] = a[(j) * D + (i)]; |
| 276 | |
| 277 | return D; |
| 278 | } |
| 279 | |
| 280 | template<typename T, qualifier Q> |
| 281 | GLM_FUNC_QUALIFIER void sortEigenvalues(vec<2, T, Q>& eigenvalues, mat<2, 2, T, Q>& eigenvectors) |
| 282 | { |
| 283 | if (eigenvalues[0] < eigenvalues[1]) |
| 284 | { |
| 285 | std::swap(eigenvalues[0], eigenvalues[1]); |
| 286 | std::swap(eigenvectors[0], eigenvectors[1]); |
| 287 | } |
| 288 | } |
| 289 | |
| 290 | template<typename T, qualifier Q> |
| 291 | GLM_FUNC_QUALIFIER void sortEigenvalues(vec<3, T, Q>& eigenvalues, mat<3, 3, T, Q>& eigenvectors) |
| 292 | { |
| 293 | if (eigenvalues[0] < eigenvalues[1]) |
| 294 | { |
| 295 | std::swap(eigenvalues[0], eigenvalues[1]); |
| 296 | std::swap(eigenvectors[0], eigenvectors[1]); |
| 297 | } |
| 298 | if (eigenvalues[0] < eigenvalues[2]) |
| 299 | { |
| 300 | std::swap(eigenvalues[0], eigenvalues[2]); |
| 301 | std::swap(eigenvectors[0], eigenvectors[2]); |
| 302 | } |
| 303 | if (eigenvalues[1] < eigenvalues[2]) |
| 304 | { |
| 305 | std::swap(eigenvalues[1], eigenvalues[2]); |
| 306 | std::swap(eigenvectors[1], eigenvectors[2]); |
| 307 | } |
| 308 | } |
| 309 | |
| 310 | template<typename T, qualifier Q> |
| 311 | GLM_FUNC_QUALIFIER void sortEigenvalues(vec<4, T, Q>& eigenvalues, mat<4, 4, T, Q>& eigenvectors) |
| 312 | { |
| 313 | if (eigenvalues[0] < eigenvalues[2]) |
| 314 | { |
| 315 | std::swap(eigenvalues[0], eigenvalues[2]); |
| 316 | std::swap(eigenvectors[0], eigenvectors[2]); |
| 317 | } |
| 318 | if (eigenvalues[1] < eigenvalues[3]) |
| 319 | { |
| 320 | std::swap(eigenvalues[1], eigenvalues[3]); |
| 321 | std::swap(eigenvectors[1], eigenvectors[3]); |
| 322 | } |
| 323 | if (eigenvalues[0] < eigenvalues[1]) |
| 324 | { |
| 325 | std::swap(eigenvalues[0], eigenvalues[1]); |
| 326 | std::swap(eigenvectors[0], eigenvectors[1]); |
| 327 | } |
| 328 | if (eigenvalues[2] < eigenvalues[3]) |
| 329 | { |
| 330 | std::swap(eigenvalues[2], eigenvalues[3]); |
| 331 | std::swap(eigenvectors[2], eigenvectors[3]); |
| 332 | } |
| 333 | if (eigenvalues[1] < eigenvalues[2]) |
| 334 | { |
| 335 | std::swap(eigenvalues[1], eigenvalues[2]); |
| 336 | std::swap(eigenvectors[1], eigenvectors[2]); |
| 337 | } |
| 338 | } |
| 339 | |
| 340 | }//namespace glm |