#ifndef __JEIGEN3D__ #define __JEIGEN3D__ #include #include #include #include #include "TMatrixDSymEigen.h" #include "JLang/JException.hh" #include "JMath/JMath.hh" #include "JMath/JLegendre.hh" #include "JGeometry3D/JVector3D.hh" #include "JGeometry3D/JQuaternion3D.hh" #include "JGeometry3D/JGeometry3DToolkit.hh" /** * \author mdejong */ namespace JGEOMETRY3D {} namespace JPP { using namespace JGEOMETRY3D; } namespace JGEOMETRY3D { /** * Eigen values in 3D. */ struct JEigenValues3D : public std::map > { /** * Constructor. * * \param __begin begin of data * \param __end end of data */ template JEigenValues3D(T __begin, T __end) { const JCenter3D center(__begin, __end); // RMS matrix TMatrixDSym A(3); A = 0.0; for (T i = __begin; i != __end; ++i) { const double dx = center.getX() - i->getX(); const double dy = center.getY() - i->getY(); const double dz = center.getZ() - i->getZ(); A(0,0) += (dx * dx); A(0,1) += (dx * dy); A(0,2) += (dx * dz); A(1,1) += (dy * dy); A(1,2) += (dy * dz); A(2,2) += (dz * dz); } A(1,0) = A(0,1); A(2,0) = A(0,2); A(2,1) = A(1,2); const TMatrixDSymEigen result(A); const TVectorD& u = result.GetEigenValues(); const TMatrixD& V = result.GetEigenVectors(); for (Int_t i = 0; i != u.GetNoElements(); ++i) { (*this)[u[i]] = JVector3D(V(0,i), V(1,i), V(2,i)); } } }; } namespace JMATH { using JGEOMETRY3D::JQuaternion3D; using JLANG::JDivisionByZero; /** * Template spacialisation for averaging quaternions. */ template<> struct JAverage { /** * Numerical precision for eigen value evaluation. */ static double MINIMAL_EIGEN_VALUE; /** * Default constructor. */ JAverage() : A(4), W(0.0) {} /** * Constructor. * * \param __begin begin of data * \param __end end of data */ template JAverage(T __begin, T __end) : A(4), W(0.0) { for (T i = __begin; i != __end; ++i) { this->put(*i); } } /** * Reset. */ void reset() { A = 0.0; W = 0.0; } /** * Type conversion operator. * * \return quaternion */ operator JQuaternion3D() const { JQuaternion3D Q; const TMatrixDSymEigen result(A); const TVectorD& u = result.GetEigenValues(); const TMatrixD& V = result.GetEigenVectors(); if (u.GetNoElements() != 0 && fabs(u[0]) > MINIMAL_EIGEN_VALUE) { Q = JQuaternion3D(V(0,0), V(1,0), V(2,0), V(3,0)); if (Q.getA() < 0.0) { Q.negate(); } Q.pow(u[0] / W); } else { //THROW(JDivisionByZero, "Invalid eigen value."); } return Q; } /** * Put quaternion. * * \param Q quaternion * \param w weight */ void put(const JQuaternion3D& Q, const double w = 1.0) { TMatrixDSym a(4); a(0,0) = Q.getA()*Q.getA(); a(0,1) = Q.getA()*Q.getB(); a(0,2) = Q.getA()*Q.getC(); a(0,3) = Q.getA()*Q.getD(); a(1,0) = a(0,1); a(1,1) = Q.getB()*Q.getB(); a(1,2) = Q.getB()*Q.getC(); a(1,3) = Q.getB()*Q.getD(); a(2,0) = a(0,2); a(2,1) = a(1,2); a(2,2) = Q.getC()*Q.getC(); a(2,3) = Q.getC()*Q.getD(); a(3,0) = a(0,3); a(3,1) = a(1,3); a(3,2) = a(2,3); a(3,3) = Q.getD()*Q.getD(); a *= w; A += a; W += w*w; } private: TMatrixDSym A; double W; }; /** * Numerical precision for eigen value evaluation. */ double JAverage::MINIMAL_EIGEN_VALUE = 1.0e-8; /** * Template specialisation for function evaluation of Legendre polynome of quaternions for undefined number of degrees. */ template<> struct JLegendre : public JLegendre_t { /** * Default constructor. */ JLegendre() {} /** * Constructor. * * \param xmin minimal abscissa value * \param xmax maximal abscissa value */ JLegendre(const double xmin, const double xmax) { this->set(xmin, xmax); } /** * Function value. * * \param x abscissa value * \return function value */ virtual JQuaternion3D getValue(const double x) const override { const double z = this->getX(x); JQuaternion3D y = zero; for (size_t n = 0; n != this->size(); ++n) { y *= pow((*this)[n], legendre(n,z)); } return y.normalise(); } }; /** * Template specialisation for function evaluation of Legendre polynome of quaternions for defined number of degrees. */ template struct JLegendre : JLegendre { using JLegendre_t::set; /** * Default constructor. */ JLegendre() { this->JLegendre::resize(N+1); } /** * Constructor. * * \param xmin minimal abscissa value * \param xmax maximal abscissa value */ JLegendre(const double xmin, const double xmax) { this->JLegendre::resize(N+1); this->set(xmin, xmax); } /** * Constructor. * * The template argument T refers to an iterator of a data structure which should have the following data members: * - double %first; // abscissa * - JQuaternion3D %second; // ordinate * * which conforms with std::pair. * * \param __begin begin of data * \param __end end of data */ template JLegendre(T __begin, T __end) { this->JLegendre::resize(N+1); this->set(__begin, __end); } /** * Set Legendre polynome of quaternions. * * The template argument T refers to an iterator of a data structure which should have the following data members: * - double %first; // abscissa * - JQuaternion3D %second; // ordinate * * which conforms with std::pair. * * \param __begin begin of data * \param __end end of data */ template void set(T __begin, T __end) { this->set(__begin, __end, __end); } /** * Set Legendre polynome of quaternions. * * The template argument T refers to an iterator of a data structure which should have the following data members: * - double %first; // abscissa * - JQuaternion3D %second; // ordinate * * which conforms with std::pair. * * \param __begin begin of data * \param __not not in data * \param __end end of data */ template void set(T __begin, T __not, T __end) { // factor to be applied as power to eigen value obtained with JAverage when degree larger than zero. static const double factor = 1.0 / (PI * 45.0 / 180.0); for (size_t n = 0; n != N + 1; ++n) { (*this)[n] = JQuaternion3D(); } this->xmin = std::numeric_limits::max(); this->xmax = std::numeric_limits::lowest(); for (T i = __begin; i != __end; ++i) { if (i != __not) { if (i->first < this->xmin) { this->xmin = i->first; } if (i->first > this->xmax) { this->xmax = i->first; } } } for (size_t n = 0; n != N + 1; ++n) { JAverage Q; for (T i = __begin; i != __end; ++i) { if (i != __not) { const double x = i->first; const JQuaternion3D& y = i->second; const double z = this->getX(x); const double w = legendre(n, z); const JQuaternion3D q = this->getValue(x).getConjugate() * y; Q.put(q, w); } } (*this)[n] = Q; if (n != 0) { (*this)[n].pow(factor); } } } /** * Read Legendre polynome from input. * * \param in input stream * \param object Legendre polynome * \return input stream */ friend inline std::istream& operator>>(std::istream& in, JLegendre& object) { for (typename JLegendre::iterator i = object.begin(); i != object.end(); ++i) { in >> *i; } return in; } private: void clear(); void resize(); void erase(); void push_back(); void pop_back(); }; } #endif