JEigen3D.hh

Go to the documentation of this file.

#ifndef __JEIGEN3D__
#define __JEIGEN3D__

#include <cmath>
#include <map>
#include <limits>
#include <functional>

#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<double, JVector3D, std::greater<double> >
  {
    /**
     * Constructor.
     *
     * \param  __begin        begin of data
     * \param  __end          end   of data
     */
    template<class T>
    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<JQuaternion3D> {

    /**
     * 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<class T>
    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<JQuaternion3D>::MINIMAL_EIGEN_VALUE = 1.0e-8;
  

  /**
   * Template specialisation for function evaluation of Legendre polynome of quaternions for undefined number of degrees.
   */
  template<>
  struct JLegendre<JQuaternion3D, (size_t) -1> :
    public JLegendre_t<JQuaternion3D>
  {
    /**
     * 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<size_t N>
  struct JLegendre<JQuaternion3D, N> :
    JLegendre<JQuaternion3D, (size_t) -1>
  {
    using JLegendre_t<JQuaternion3D>::set;


    /**
     * Default constructor.
     */
    JLegendre()
    {
      this->JLegendre<JQuaternion3D>::resize(N+1);
    }


    /**
     * Constructor.
     *
     * \param  xmin             minimal abscissa value
     * \param  xmax             maximal abscissa value
     */
    JLegendre(const double xmin,
              const double xmax)
    {
      this->JLegendre<JQuaternion3D>::resize(N+1);
      this->set(xmin, xmax);
    }

    
    /**
     * Constructor.
     *
     * The template argument <tt>T</tt> 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<class T>
    JLegendre(T __begin, T __end)
    {
      this->JLegendre<JQuaternion3D>::resize(N+1);
      this->set(__begin, __end);
    }


    /**
     * Set Legendre polynome of quaternions.
     *
     * The template argument <tt>T</tt> 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<class T>
    void set(T __begin, T __end)
    {
      this->set(__begin, __end, __end);
    }


    /**
     * Set Legendre polynome of quaternions.
     *
     * The template argument <tt>T</tt> 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<class T>
    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<double>::max();
      this->xmax = std::numeric_limits<double>::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<JQuaternion3D> 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