From 1ab75b5b7bc202e2cde741dde7496de23836f5ca Mon Sep 17 00:00:00 2001 From: Andy Maloney Date: Sat, 25 Apr 2020 22:11:06 -0400 Subject: [PATCH] Add CCCoreLib namespace to three files Addresses part of #26 --- include/Garbage.h | 161 +++++----- include/Jacobi.h | 716 ++++++++++++++++++++++---------------------- include/RayAndBox.h | 185 ++++++------ 3 files changed, 533 insertions(+), 529 deletions(-) diff --git a/include/Garbage.h b/include/Garbage.h index 1eaae75..b4500b6 100644 --- a/include/Garbage.h +++ b/include/Garbage.h @@ -9,98 +9,97 @@ //STL #include -//! Garbage container (automatically deletes pointers when destroyed) -template class Garbage +namespace CCCoreLib { -public: - //! Puts an item in the trash - inline void add(C* item) + //! Garbage container (automatically deletes pointers when destroyed) + template class Garbage { - try + public: + //! Puts an item in the trash + inline void add(C* item) { - m_items.insert(item); + try + { + m_items.insert(item); + } + catch (const std::bad_alloc&) + { + //what can we do?! + } } - catch (const std::bad_alloc&) + + //! Removes an item from the trash + /** \warning The item won't be destroyed! **/ + inline void remove(C* item) { - //what can we do?! + m_items.erase(item); } - } - - //! Removes an item from the trash - /** \warning The item won't be destroyed! - **/ - inline void remove(C* item) - { - m_items.erase(item); - } - - //! To manually delete an item already in the trash - inline void destroy(C* item) - { - m_items.erase(item); - delete item; - } - - //! Destructor - /** Automatically deletes all items - **/ - ~Garbage() - { - //dispose of left over - for (auto it = m_items.begin(); it != m_items.end(); ++it) - delete *it; - m_items.clear(); - } - - //! Items to delete - std::unordered_set m_items; -}; - -//! Speciailization for ScalarFields -template <> class Garbage -{ -public: - //! Puts an item in the trash - inline void add(CCCoreLib::ScalarField* item) - { - try + + //! To manually delete an item already in the trash + inline void destroy(C* item) { - m_items.insert(item); + m_items.erase(item); + delete item; } - catch (const std::bad_alloc&) + + //! Destructor + /** Automatically deletes all items **/ + ~Garbage() { - //what can we do?! + //dispose of left over + for (auto it = m_items.begin(); it != m_items.end(); ++it) + delete *it; + m_items.clear(); } - } - - //! Removes an item from the trash - /** \warning The item won't be destroyed! - **/ - inline void remove(CCCoreLib::ScalarField* item) + + //! Items to delete + std::unordered_set m_items; + }; + + //! Specialization for ScalarFields + template <> class Garbage { - m_items.erase(item); - } - - //! Manually deltes an item already in the trash - inline void destroy(CCCoreLib::ScalarField* item) - { - m_items.erase(item); - item->release(); - } - - //! Destructor - /** Automatically deletes all items - **/ - ~Garbage() - { - //dispose of left over - for (auto item : m_items) + public: + //! Puts an item in the trash + inline void add(ScalarField* item) + { + try + { + m_items.insert(item); + } + catch (const std::bad_alloc&) + { + //what can we do?! + } + } + + //! Removes an item from the trash + /** \warning The item won't be destroyed! **/ + inline void remove(ScalarField* item) { + m_items.erase(item); + } + + //! Manually deltes an item already in the trash + inline void destroy(ScalarField* item) + { + m_items.erase(item); item->release(); } - m_items.clear(); - } - - //! Items to delete - std::unordered_set m_items; -}; + + //! Destructor + /** Automatically deletes all items **/ + ~Garbage() + { + //dispose of left over + for (auto item : m_items) + { + item->release(); + } + m_items.clear(); + } + + //! Items to delete + std::unordered_set m_items; + }; +} diff --git a/include/Jacobi.h b/include/Jacobi.h index 8768ee5..1c126eb 100644 --- a/include/Jacobi.h +++ b/include/Jacobi.h @@ -6,447 +6,449 @@ //Local #include "SquareMatrix.h" -#define ROTATE(a,i,j,k,l) { Scalar g = a[i][j]; h = a[k][l]; a[i][j] = g-s*(h+g*tau); a[k][l] = h+s*(g-h*tau); } - -//! Jacobi eigen vectors/values decomposition -template class Jacobi +namespace CCCoreLib { -public: - - using SquareMatrix = CCCoreLib::SquareMatrixTpl; - using EigenValues = std::vector; - - //! Computes the eigenvalues and eigenvectors of a given square matrix - /** It uses Rutishauser's modfications of the classical Jacobi rotation method with threshold pivoting. - +#define ROTATE(a,i,j,k,l) { Scalar g = a[i][j]; h = a[k][l]; a[i][j] = g-s*(h+g*tau); a[k][l] = h+s*(g-h*tau); } + + //! Jacobi eigen vectors/values decomposition + template class Jacobi + { + public: + + using SquareMatrix = SquareMatrixTpl; + using EigenValues = std::vector; + + //! Computes the eigenvalues and eigenvectors of a given square matrix + /** It uses Rutishauser's modfications of the classical Jacobi rotation method with threshold pivoting. + Note: this code is inspired from John Burkardt's 'jacobi_eigenvalue' code See https://people.sc.fsu.edu/~jburkardt/cpp_src/jacobi_eigenvalue/jacobi_eigenvalue.cpp - + \warning DGM: this method gives strange results in some particular cases!!! - + \param[in] matrix input square matrix \param[out] eigenVectors eigenvectors (as a square matrix) \param[out] eigenValues eigenvalues \param[in] maxIterationCount max number of iteration (optional) \return success **/ - static bool ComputeEigenValuesAndVectors2(const SquareMatrix& matrix, SquareMatrix& eigenVectors, EigenValues& eigenValues, unsigned maxIterationCount = 50) - { - if (!matrix.isValid()) - { - return false; - } - - unsigned n = matrix.size(); - - //duplicate the current matrix (as it may be modified afterwards!) - SquareMatrix a(matrix); - - //output eigen vectors matrix - eigenVectors = SquareMatrix(n); - if (!eigenVectors.isValid()) - { - return false; - } - eigenVectors.toIdentity(); - - std::vector bw, zw; - try + static bool ComputeEigenValuesAndVectors2(const SquareMatrix& matrix, SquareMatrix& eigenVectors, EigenValues& eigenValues, unsigned maxIterationCount = 50) { - eigenValues.resize(n); - zw.resize(n, 0); - bw.resize(n); - - //initialize bw and d to the diagonal of this matrix - for (unsigned i = 0; i < n; i++) + if (!matrix.isValid()) { - bw[i] = eigenValues[i] = a.m_values[i][i]; + return false; } - } - catch (const std::bad_alloc&) - { - //not enough memory - return false; - } - - for (unsigned it = 0; it < maxIterationCount; ++it) - { - //The convergence threshold is based on the size of the elements in - //the strict upper triangle of the matrix. - double sum = 0.0; - for (unsigned j = 0; j < n; ++j) + + unsigned n = matrix.size(); + + //duplicate the current matrix (as it may be modified afterwards!) + SquareMatrix a(matrix); + + //output eigen vectors matrix + eigenVectors = SquareMatrix(n); + if (!eigenVectors.isValid()) { - for (unsigned i = 0; i < j; ++i) + return false; + } + eigenVectors.toIdentity(); + + std::vector bw, zw; + try + { + eigenValues.resize(n); + zw.resize(n, 0); + bw.resize(n); + + //initialize bw and d to the diagonal of this matrix + for (unsigned i = 0; i < n; i++) { - double val = a.m_values[i][j]; - sum += val * val; + bw[i] = eigenValues[i] = a.m_values[i][i]; } } - - Scalar thresh = static_cast(sqrt(sum) / (4 * n)); - if (thresh == 0) + catch (const std::bad_alloc&) { - break; + //not enough memory + return false; } - - for (unsigned p = 0; p < n; ++p) + + for (unsigned it = 0; it < maxIterationCount; ++it) { - for (unsigned q = p + 1; q < n; ++q) + //The convergence threshold is based on the size of the elements in + //the strict upper triangle of the matrix. + double sum = 0.0; + for (unsigned j = 0; j < n; ++j) { - Scalar gapq = 10 * std::abs(a.m_values[p][q]); - Scalar termp = gapq + std::abs(eigenValues[p]); - Scalar termq = gapq + std::abs(eigenValues[q]); - - //Annihilate tiny offdiagonal elements - if (it > 4 - && termp == std::abs(eigenValues[p]) - && termq == std::abs(eigenValues[q])) + for (unsigned i = 0; i < j; ++i) { - a.m_values[p][q] = 0; + double val = a.m_values[i][j]; + sum += val * val; } - //Otherwise, apply a rotation - else if (thresh <= std::abs(a.m_values[p][q])) + } + + Scalar thresh = static_cast(sqrt(sum) / (4 * n)); + if (thresh == 0) + { + break; + } + + for (unsigned p = 0; p < n; ++p) + { + for (unsigned q = p + 1; q < n; ++q) { - Scalar h = eigenValues[q] - eigenValues[p]; - Scalar term = std::abs(h) + gapq; - - Scalar t = 0; - if (term == std::abs(h)) + Scalar gapq = 10 * std::abs(a.m_values[p][q]); + Scalar termp = gapq + std::abs(eigenValues[p]); + Scalar termq = gapq + std::abs(eigenValues[q]); + + //Annihilate tiny offdiagonal elements + if (it > 4 + && termp == std::abs(eigenValues[p]) + && termq == std::abs(eigenValues[q])) { - t = a.m_values[p][q] / h; + a.m_values[p][q] = 0; } - else + //Otherwise, apply a rotation + else if (thresh <= std::abs(a.m_values[p][q])) { - Scalar theta = h / (2 * a.m_values[p][q]); - t = 1 / (std::abs(theta) + sqrt(1 + theta*theta)); - if (theta < 0) + Scalar h = eigenValues[q] - eigenValues[p]; + Scalar term = std::abs(h) + gapq; + + Scalar t = 0; + if (term == std::abs(h)) { - t = -t; + t = a.m_values[p][q] / h; + } + else + { + Scalar theta = h / (2 * a.m_values[p][q]); + t = 1 / (std::abs(theta) + sqrt(1 + theta*theta)); + if (theta < 0) + { + t = -t; + } + } + Scalar c = 1 / sqrt(1 + t*t); + Scalar s = t * c; + Scalar tau = s / (1 + c); + h = t * a.m_values[q][p]; + + //Accumulate corrections to diagonal elements. + zw[p] -= h; + zw[q] += h; + eigenValues[p] -= h; + eigenValues[q] += h; + + a.m_values[p][q] = 0; + + // Rotate, using information from the upper triangle of A only. + unsigned j; + for (j = 0; j < p; ++j) + { + Scalar g = a.m_values[j][p]; + Scalar h = a.m_values[j][q]; + a.m_values[j][p] = g - s * (h + g * tau); + a.m_values[j][q] = h + s * (g - h * tau); + } + + for (j = p + 1; j < q; ++j) + { + Scalar g = a.m_values[p][j]; + Scalar h = a.m_values[j][q]; + a.m_values[p][j] = g - s * (h + g * tau); + a.m_values[j][q] = h + s * (g - h * tau); + } + + for (j = q + 1; j < n; ++j) + { + Scalar g = a.m_values[p][j]; + Scalar h = a.m_values[q][j]; + a.m_values[p][j] = g - s * (h + g * tau); + a.m_values[q][j] = h + s * (g - h * tau); + } + + //Accumulate information in the eigenvector matrix. + for (j = 0; j < n; ++j) + { + Scalar g = eigenVectors.m_values[j][p]; + Scalar h = eigenVectors.m_values[j][q]; + eigenVectors.m_values[j][p] = g - s * (h + g * tau); + eigenVectors.m_values[j][q] = h + s * (g - h * tau); } - } - Scalar c = 1 / sqrt(1 + t*t); - Scalar s = t * c; - Scalar tau = s / (1 + c); - h = t * a.m_values[q][p]; - - //Accumulate corrections to diagonal elements. - zw[p] -= h; - zw[q] += h; - eigenValues[p] -= h; - eigenValues[q] += h; - - a.m_values[p][q] = 0; - - // Rotate, using information from the upper triangle of A only. - unsigned j; - for (j = 0; j < p; ++j) - { - Scalar g = a.m_values[j][p]; - Scalar h = a.m_values[j][q]; - a.m_values[j][p] = g - s * (h + g * tau); - a.m_values[j][q] = h + s * (g - h * tau); - } - - for (j = p + 1; j < q; ++j) - { - Scalar g = a.m_values[p][j]; - Scalar h = a.m_values[j][q]; - a.m_values[p][j] = g - s * (h + g * tau); - a.m_values[j][q] = h + s * (g - h * tau); - } - - for (j = q + 1; j < n; ++j) - { - Scalar g = a.m_values[p][j]; - Scalar h = a.m_values[q][j]; - a.m_values[p][j] = g - s * (h + g * tau); - a.m_values[q][j] = h + s * (g - h * tau); - } - - //Accumulate information in the eigenvector matrix. - for (j = 0; j < n; ++j) - { - Scalar g = eigenVectors.m_values[j][p]; - Scalar h = eigenVectors.m_values[j][q]; - eigenVectors.m_values[j][p] = g - s * (h + g * tau); - eigenVectors.m_values[j][q] = h + s * (g - h * tau); } } } + + for (unsigned i = 0; i < n; ++i) + { + bw[i] += zw[i]; + eigenValues[i] = bw[i]; + zw[i] = 0.0; + } } - - for (unsigned i = 0; i < n; ++i) - { - bw[i] += zw[i]; - eigenValues[i] = bw[i]; - zw[i] = 0.0; - } - } - - return true; - } - - //! Computes eigen vectors (and values) with the Jacobian method - /** See the Numerical Recipes. - **/ - static bool ComputeEigenValuesAndVectors( const SquareMatrix& matrix, - SquareMatrix& eigenVectors, - EigenValues& eigenValues, - bool absoluteValues = true, - unsigned maxIterationCount = 50) - { - if (!matrix.isValid()) - { - return false; - } - - unsigned n = matrix.size(); - unsigned matrixSquareSize = n * n; - - //output eigen vectors matrix - eigenVectors = SquareMatrix(n); - if (!eigenVectors.isValid()) - { - return false; - } - eigenVectors.toIdentity(); - - std::vector b, z; - try - { - eigenValues.resize(n); - b.resize(n); - z.resize(n); - } - catch (const std::bad_alloc&) - { - //not enough memory - return false; + + return true; } - Scalar* d = eigenValues.data(); - - //init + + //! Computes eigen vectors (and values) with the Jacobian method + /** See the Numerical Recipes. **/ + static bool ComputeEigenValuesAndVectors( const SquareMatrix& matrix, + SquareMatrix& eigenVectors, + EigenValues& eigenValues, + bool absoluteValues = true, + unsigned maxIterationCount = 50) { - for (unsigned ip = 0; ip < n; ip++) + if (!matrix.isValid()) { - b[ip] = d[ip] = matrix.m_values[ip][ip]; //Initialize b and d to the diagonal of a. - z[ip] = 0; //This vector will accumulate terms of the form tapq as in equation (11.1.14) + return false; } - } - - for (unsigned i = 1; i <= maxIterationCount; i++) - { - //Sum off-diagonal elements - Scalar sm = 0; + + unsigned n = matrix.size(); + unsigned matrixSquareSize = n * n; + + //output eigen vectors matrix + eigenVectors = SquareMatrix(n); + if (!eigenVectors.isValid()) { - for (unsigned ip = 0; ip < n - 1; ip++) - { - for (unsigned iq = ip + 1; iq < n; iq++) - sm += std::abs(matrix.m_values[ip][iq]); - } + return false; } - - if (sm == 0) //The normal return, which relies on quadratic convergence to machine underflow. + eigenVectors.toIdentity(); + + std::vector b, z; + try { - if (absoluteValues) - { - //we only need the absolute values of eigenvalues - for (unsigned ip = 0; ip < n; ip++) - d[ip] = std::abs(d[ip]); - } - - return true; + eigenValues.resize(n); + b.resize(n); + z.resize(n); } - - Scalar tresh = 0; - if (i < 4) + catch (const std::bad_alloc&) { - tresh = sm / static_cast(5 * matrixSquareSize); //...on the first three sweeps. + //not enough memory + return false; } - - for (unsigned ip = 0; ip < n - 1; ip++) + Scalar* d = eigenValues.data(); + + //init + { + for (unsigned ip = 0; ip < n; ip++) + { + b[ip] = d[ip] = matrix.m_values[ip][ip]; //Initialize b and d to the diagonal of a. + z[ip] = 0; //This vector will accumulate terms of the form tapq as in equation (11.1.14) + } + } + + for (unsigned i = 1; i <= maxIterationCount; i++) { - for (unsigned iq = ip + 1; iq < n; iq++) + //Sum off-diagonal elements + Scalar sm = 0; { - Scalar g = std::abs(matrix.m_values[ip][iq]) * 100; - //After four sweeps, skip the rotation if the off-diagonal element is small. - if (i > 4 - && static_cast(std::abs(d[ip]) + g) == static_cast(std::abs(d[ip])) - && static_cast(std::abs(d[iq]) + g) == static_cast(std::abs(d[iq]))) + for (unsigned ip = 0; ip < n - 1; ip++) { - matrix.m_values[ip][iq] = 0; + for (unsigned iq = ip + 1; iq < n; iq++) + sm += std::abs(matrix.m_values[ip][iq]); } - else if (std::abs(matrix.m_values[ip][iq]) > tresh) + } + + if (sm == 0) //The normal return, which relies on quadratic convergence to machine underflow. + { + if (absoluteValues) { - Scalar h = d[iq] - d[ip]; - Scalar t = 0; - if (static_cast(std::abs(h) + g) == static_cast(std::abs(h))) - { - t = matrix.m_values[ip][iq] / h; - } - else - { - Scalar theta = h / (2 * matrix.m_values[ip][iq]); //Equation (11.1.10). - t = 1 / (std::abs(theta) + sqrt(1 + theta*theta)); - if (theta < 0) - t = -t; - } - - Scalar c = 1 / sqrt(t*t + 1); - Scalar s = t*c; - Scalar tau = s / (1 + c); - h = t * matrix.m_values[ip][iq]; - z[ip] -= h; - z[iq] += h; - d[ip] -= h; - d[iq] += h; - matrix.m_values[ip][iq] = 0; - - //Case of rotations 1 <= j < p - { - for (unsigned j = 0; j + 1 <= ip; j++) - ROTATE(matrix.m_values, j, ip, j, iq) - } - //Case of rotations p < j < q - { - for (unsigned j = ip + 1; j + 1 <= iq; j++) - ROTATE(matrix.m_values, ip, j, j, iq) - } - //Case of rotations q < j <= n + //we only need the absolute values of eigenvalues + for (unsigned ip = 0; ip < n; ip++) + d[ip] = std::abs(d[ip]); + } + + return true; + } + + Scalar tresh = 0; + if (i < 4) + { + tresh = sm / static_cast(5 * matrixSquareSize); //...on the first three sweeps. + } + + for (unsigned ip = 0; ip < n - 1; ip++) + { + for (unsigned iq = ip + 1; iq < n; iq++) + { + Scalar g = std::abs(matrix.m_values[ip][iq]) * 100; + //After four sweeps, skip the rotation if the off-diagonal element is small. + if (i > 4 + && static_cast(std::abs(d[ip]) + g) == static_cast(std::abs(d[ip])) + && static_cast(std::abs(d[iq]) + g) == static_cast(std::abs(d[iq]))) { - for (unsigned j = iq + 1; j < n; j++) - ROTATE(matrix.m_values, ip, j, iq, j) + matrix.m_values[ip][iq] = 0; } - //Last case + else if (std::abs(matrix.m_values[ip][iq]) > tresh) { - for (unsigned j = 0; j < n; j++) - ROTATE(eigenVectors.m_values, j, ip, j, iq) + Scalar h = d[iq] - d[ip]; + Scalar t = 0; + if (static_cast(std::abs(h) + g) == static_cast(std::abs(h))) + { + t = matrix.m_values[ip][iq] / h; + } + else + { + Scalar theta = h / (2 * matrix.m_values[ip][iq]); //Equation (11.1.10). + t = 1 / (std::abs(theta) + sqrt(1 + theta*theta)); + if (theta < 0) + t = -t; + } + + Scalar c = 1 / sqrt(t*t + 1); + Scalar s = t*c; + Scalar tau = s / (1 + c); + h = t * matrix.m_values[ip][iq]; + z[ip] -= h; + z[iq] += h; + d[ip] -= h; + d[iq] += h; + matrix.m_values[ip][iq] = 0; + + //Case of rotations 1 <= j < p + { + for (unsigned j = 0; j + 1 <= ip; j++) + ROTATE(matrix.m_values, j, ip, j, iq) + } + //Case of rotations p < j < q + { + for (unsigned j = ip + 1; j + 1 <= iq; j++) + ROTATE(matrix.m_values, ip, j, j, iq) + } + //Case of rotations q < j <= n + { + for (unsigned j = iq + 1; j < n; j++) + ROTATE(matrix.m_values, ip, j, iq, j) + } + //Last case + { + for (unsigned j = 0; j < n; j++) + ROTATE(eigenVectors.m_values, j, ip, j, iq) + } } } } - } - - //update b, d and z - { - for (unsigned ip = 0; ip < n; ip++) + + //update b, d and z { - b[ip] += z[ip]; - d[ip] = b[ip]; - z[ip] = 0; + for (unsigned ip = 0; ip < n; ip++) + { + b[ip] += z[ip]; + d[ip] = b[ip]; + z[ip] = 0; + } } } + + //Too many iterations! + return false; } - - //Too many iterations! - return false; - } - - //! Sorts the eigenvectors in the decreasing order of their associated eigenvalues - /** \param eigenVectors eigenvectors (as a square matrix) + + //! Sorts the eigenvectors in the decreasing order of their associated eigenvalues + /** \param eigenVectors eigenvectors (as a square matrix) \param eigenValues eigenvalues \return success **/ - static bool SortEigenValuesAndVectors(SquareMatrix& eigenVectors, EigenValues& eigenValues) - { - if (!eigenVectors.isValid() - || eigenVectors.size() < 2 - || eigenVectors.size() != eigenValues.size()) + static bool SortEigenValuesAndVectors(SquareMatrix& eigenVectors, EigenValues& eigenValues) { - assert(false); - return false; - } - - unsigned n = eigenVectors.size(); - for (unsigned i = 0; i < n - 1; i++) - { - unsigned maxValIndex = i; - for (unsigned j = i + 1; j eigenValues[maxValIndex]) - maxValIndex = j; - - if (maxValIndex != i) + if (!eigenVectors.isValid() + || eigenVectors.size() < 2 + || eigenVectors.size() != eigenValues.size()) { - std::swap(eigenValues[i], eigenValues[maxValIndex]); - for (unsigned j = 0; j < n; ++j) - std::swap(eigenVectors.m_values[j][i], eigenVectors.m_values[j][maxValIndex]); + assert(false); + return false; } + + unsigned n = eigenVectors.size(); + for (unsigned i = 0; i < n - 1; i++) + { + unsigned maxValIndex = i; + for (unsigned j = i + 1; j eigenValues[maxValIndex]) + maxValIndex = j; + + if (maxValIndex != i) + { + std::swap(eigenValues[i], eigenValues[maxValIndex]); + for (unsigned j = 0; j < n; ++j) + std::swap(eigenVectors.m_values[j][i], eigenVectors.m_values[j][maxValIndex]); + } + } + + return true; } - - return true; - } - - //! Returns the given eigenvector - /** \param eigenVectors eigenvectors (as a square matrix) + + //! Returns the given eigenvector + /** \param eigenVectors eigenvectors (as a square matrix) \param index requested eigenvector index (< eigenvectors matrix size) \param eigenVector output vector (size = matrix size) \return success **/ - static bool GetEigenVector(const SquareMatrix& eigenVectors, unsigned index, Scalar eigenVector[]) - { - if (eigenVector && index < eigenVectors.size()) + static bool GetEigenVector(const SquareMatrix& eigenVectors, unsigned index, Scalar eigenVector[]) { - for (unsigned i = 0; i < eigenVectors.size(); ++i) + if (eigenVector && index < eigenVectors.size()) { - eigenVector[i] = eigenVectors.m_values[i][index]; + for (unsigned i = 0; i < eigenVectors.size(); ++i) + { + eigenVector[i] = eigenVectors.m_values[i][index]; + } + return true; + } + else + { + assert(false); + return false; } - return true; - } - else - { - assert(false); - return false; } - } - - //! Returns the biggest eigenvalue and its associated eigenvector - /** \param eigenVectors eigenvectors (as a square matrix) + + //! Returns the biggest eigenvalue and its associated eigenvector + /** \param eigenVectors eigenvectors (as a square matrix) \param eigenValues eigenvalues \param maxEigenValue biggest eigenvalue \param maxEigenVector eigenvector vector corresponding to the biggest eigenvalue \return success **/ - static bool GetMaxEigenValueAndVector(const SquareMatrix& eigenVectors, const EigenValues& eigenValues, Scalar& maxEigenValue, Scalar maxEigenVector[]) - { - if (!eigenVectors.isValid() - || eigenVectors.size() < 2 - || eigenVectors.size() != eigenValues.size()) + static bool GetMaxEigenValueAndVector(const SquareMatrix& eigenVectors, const EigenValues& eigenValues, Scalar& maxEigenValue, Scalar maxEigenVector[]) { - assert(false); - return false; + if (!eigenVectors.isValid() + || eigenVectors.size() < 2 + || eigenVectors.size() != eigenValues.size()) + { + assert(false); + return false; + } + + unsigned maxIndex = 0; + for (unsigned i = 1; i eigenValues[maxIndex]) + maxIndex = i; + + maxEigenValue = eigenValues[maxIndex]; + return GetEigenVector(eigenVectors, maxIndex, maxEigenVector); } - - unsigned maxIndex = 0; - for (unsigned i = 1; i eigenValues[maxIndex]) - maxIndex = i; - - maxEigenValue = eigenValues[maxIndex]; - return GetEigenVector(eigenVectors, maxIndex, maxEigenVector); - } - - //! Returns the smallest eigenvalue and its associated eigenvector - /** \param eigenVectors eigenvectors (as a square matrix) + + //! Returns the smallest eigenvalue and its associated eigenvector + /** \param eigenVectors eigenvectors (as a square matrix) \param eigenValues eigenvalues \param minEigenValue smallest eigenvalue \param minEigenVector eigenvector vector corresponding to the smallest eigenvalue \return success **/ - static bool GetMinEigenValueAndVector(const SquareMatrix& eigenVectors, const EigenValues& eigenValues, Scalar& minEigenValue, Scalar minEigenVector[]) - { - if (!eigenVectors.isValid() - || eigenVectors.size() < 2 - || eigenVectors.size() != eigenValues.size()) + static bool GetMinEigenValueAndVector(const SquareMatrix& eigenVectors, const EigenValues& eigenValues, Scalar& minEigenValue, Scalar minEigenVector[]) { - assert(false); - return false; + if (!eigenVectors.isValid() + || eigenVectors.size() < 2 + || eigenVectors.size() != eigenValues.size()) + { + assert(false); + return false; + } + + unsigned minIndex = 0; + for (unsigned i = 1; i < eigenVectors.size(); ++i) + if (eigenValues[i] < eigenValues[minIndex]) + minIndex = i; + + minEigenValue = eigenValues[minIndex]; + return GetEigenVector(eigenVectors, minIndex, minEigenVector); } - - unsigned minIndex = 0; - for (unsigned i = 1; i < eigenVectors.size(); ++i) - if (eigenValues[i] < eigenValues[minIndex]) - minIndex = i; - - minEigenValue = eigenValues[minIndex]; - return GetEigenVector(eigenVectors, minIndex, minEigenVector); - } -}; + }; +} diff --git a/include/RayAndBox.h b/include/RayAndBox.h index 3cce819..bf09861 100644 --- a/include/RayAndBox.h +++ b/include/RayAndBox.h @@ -5,97 +5,100 @@ #include "CCGeom.h" -//! Simple Ray structure -template struct Ray +namespace CCCoreLib { - Ray(const Vector3Tpl& rayAxis, const Vector3Tpl& rayOrigin) - : dir(rayAxis) - , origin(rayOrigin) - , invDir(0,0,0) - , sign(0,0,0) + //! Simple Ray structure + template struct Ray { - dir.normalize(); - invDir = Vector3Tpl(1/rayAxis.x, 1/rayAxis.y, 1/rayAxis.z); // +/-infinity is acceptable here because we are mainly interested in the sign - sign = Tuple3i(invDir.x < 0, invDir.y < 0, invDir.z < 0); - } - - inline double radialSquareDistance(const Vector3Tpl& P) const - { - Vector3Tpl OP = P - origin; - return OP.cross(dir).norm2d(); - } - - inline double squareDistanceToOrigin(const Vector3Tpl& P) const - { - Vector3Tpl OP = P - origin; - return OP.norm2d(); - } - - inline void squareDistances(const Vector3Tpl& P, double& radial, double& toOrigin) const - { - Vector3Tpl OP = P - origin; - radial = OP.cross(dir).norm2d(); - toOrigin = OP.norm2d(); - } - - Vector3Tpl dir, origin; - Vector3Tpl invDir; - Tuple3i sign; -}; - -//! Simple axis aligned box structure -template struct AABB -{ - AABB(const Vector3Tpl& minCorner, const Vector3Tpl& maxCorner) + Ray(const Vector3Tpl& rayAxis, const Vector3Tpl& rayOrigin) + : dir(rayAxis) + , origin(rayOrigin) + , invDir(0,0,0) + , sign(0,0,0) + { + dir.normalize(); + invDir = Vector3Tpl(1/rayAxis.x, 1/rayAxis.y, 1/rayAxis.z); // +/-infinity is acceptable here because we are mainly interested in the sign + sign = Tuple3i(invDir.x < 0, invDir.y < 0, invDir.z < 0); + } + + double radialSquareDistance(const Vector3Tpl& P) const + { + Vector3Tpl OP = P - origin; + return OP.cross(dir).norm2d(); + } + + double squareDistanceToOrigin(const Vector3Tpl& P) const + { + Vector3Tpl OP = P - origin; + return OP.norm2d(); + } + + void squareDistances(const Vector3Tpl& P, double& radial, double& toOrigin) const + { + Vector3Tpl OP = P - origin; + radial = OP.cross(dir).norm2d(); + toOrigin = OP.norm2d(); + } + + Vector3Tpl dir, origin; + Vector3Tpl invDir; + Tuple3i sign; + }; + + //! Simple axis aligned box structure + template struct AABB { - assert(minCorner.x <= maxCorner.x); - assert(minCorner.y <= maxCorner.y); - assert(minCorner.z <= maxCorner.z); - - corners[0] = minCorner; - corners[1] = maxCorner; - } - - /* - * Ray-box intersection using IEEE numerical properties to ensure that the - * test is both robust and efficient, as described in: - * - * Amy Williams, Steve Barrus, R. Keith Morley, and Peter Shirley - * "An Efficient and Robust Ray-Box Intersection Algorithm" - * Journal of graphics tools, 10(1):49-54, 2005 - * - */ - bool intersects(const Ray &r, T* t0 = 0, T* t1 = 0) const - { - T tmin = (corners[ r.sign.x].x - r.origin.x) * r.invDir.x; - T tmax = (corners[1-r.sign.x].x - r.origin.x) * r.invDir.x; - T tymin = (corners[ r.sign.y].y - r.origin.y) * r.invDir.y; - T tymax = (corners[1-r.sign.y].y - r.origin.y) * r.invDir.y; - - if (tmin > tymax || tymin > tmax) - return false; - if (tymin > tmin) - tmin = tymin; - if (tymax < tmax) - tmax = tymax; - - T tzmin = (corners[ r.sign.z].z - r.origin.z) * r.invDir.z; - T tzmax = (corners[1-r.sign.z].z - r.origin.z) * r.invDir.z; - - if (tmin > tzmax || tzmin > tmax) - return false; - if (tzmin > tmin) - tmin = tzmin; - if (tzmax < tmax) - tmax = tzmax; - - if (t0) - *t0 = tmin; - if (t1) - *t1 = tmax; - - return true; - } - - Vector3Tpl corners[2]; -}; + AABB(const Vector3Tpl& minCorner, const Vector3Tpl& maxCorner) + { + assert(minCorner.x <= maxCorner.x); + assert(minCorner.y <= maxCorner.y); + assert(minCorner.z <= maxCorner.z); + + corners[0] = minCorner; + corners[1] = maxCorner; + } + + /* + * Ray-box intersection using IEEE numerical properties to ensure that the + * test is both robust and efficient, as described in: + * + * Amy Williams, Steve Barrus, R. Keith Morley, and Peter Shirley + * "An Efficient and Robust Ray-Box Intersection Algorithm" + * Journal of graphics tools, 10(1):49-54, 2005 + * + */ + bool intersects(const Ray &r, T* t0 = 0, T* t1 = 0) const + { + T tmin = (corners[ r.sign.x].x - r.origin.x) * r.invDir.x; + T tmax = (corners[1-r.sign.x].x - r.origin.x) * r.invDir.x; + T tymin = (corners[ r.sign.y].y - r.origin.y) * r.invDir.y; + T tymax = (corners[1-r.sign.y].y - r.origin.y) * r.invDir.y; + + if (tmin > tymax || tymin > tmax) + return false; + if (tymin > tmin) + tmin = tymin; + if (tymax < tmax) + tmax = tymax; + + T tzmin = (corners[ r.sign.z].z - r.origin.z) * r.invDir.z; + T tzmax = (corners[1-r.sign.z].z - r.origin.z) * r.invDir.z; + + if (tmin > tzmax || tzmin > tmax) + return false; + if (tzmin > tmin) + tmin = tzmin; + if (tzmax < tmax) + tmax = tzmax; + + if (t0) + *t0 = tmin; + if (t1) + *t1 = tmax; + + return true; + } + + Vector3Tpl corners[2]; + }; +}