//
// This file is part of MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2015 Richard S. Smith and collaborators.
//
// If you use MorphoGraphX in your work, please cite:
//   http://dx.doi.org/10.7554/eLife.05864
//
// MorphoGraphX is free software, and is licensed under under the terms of the 
// GNU General (GPL) Public License version 2.0, http://www.gnu.org/licenses.
// 
#ifndef QUATERNION_HPP
#define QUATERNION_HPP

/// \file Quaternion.hpp Implements the quaternion object

#include <Config.hpp>

#include <Geometry.hpp>
#include <Matrix.hpp>
#include <Vector.hpp>

namespace mgx 
{
  /**
   * \class Quaternion Quaternion.hpp <Quaternion.hpp>
   *
   * Implements the quaternion operations
   */
  struct mgx_EXPORT Quaternion : public Point4d {
    /**
     * Default constructor.
     *
     * Provides an identity quaternion FIXME Shouldn't this be 0?
     */
    Quaternion() : Point4d(0, 0, 0, 1) {}
  
    /**
     * Creates a quaternion from a Point4d
     */
    Quaternion(const Point4d &other) : Point4d(other) {}
  
    /**
     * Creates a quaternion specified by its components
     */
    Quaternion(double x, double y, double z, double w) : Point4d(x, y, z, w) {}
  
    /**
     * Copy constructor
     */
    Quaternion(const Quaternion &other) : Point4d(other) {}
  
    /**
     * Creates a Quaternion corresponding to an axis rotation
     *
     * \param axis Axis of the rotation. It needs not be normalized before
     * hand. If it is null, then the Quaternion will correspond to the identity
     * matrix.
     * \param angle Angle of the rotation.
     */
    Quaternion(const Point3d &axis, double angle);
  
    /**
     * Creates the quaternion corresponding to the rotation transforming \p
     * from into \p to.
     */
    Quaternion(const Point3d &from, const Point3d &to);
  
    /**
     * Creates the quaternion corresponding to the rotation matrix \p m
     */
    Quaternion(const Matrix3d &m);
  
    /**
     * Assignment operator for quaternions
     */
    Quaternion &operator=(const Quaternion &other);
  
    /**
     * Set the quaternion to the described rotation
     * \param axis Axis of the rotation
     * \param angle Angle of the rotation
     */
    void setAxisAngle(const Point3d &axis, double angle);
  
    /**
     * Accessing the real part of the quaternion
     */
    double &w() { return elems[3]; }
    /**
     * Accessing the real part of the quaternion
     */
    const double &w() const { return elems[3]; }
  
    /**
     * Quaternion in-place addition
     */
    Quaternion &operator+=(const Quaternion &other)
    {
      for(size_t i = 0; i < 4; ++i)
        elems[i] += other.elems[i];
      return *this;
    }
  
    /**
     * Quaternion addition
     */
    Quaternion operator+(const Quaternion &other) const
    {
      Quaternion tmp(*this);
      tmp += other;
      return tmp;
    }
  
    /**
     * Quaternion multiplication
     */
    Quaternion operator*(const Quaternion &other) const;
    /**
     * Quaternion in-place multiplication
     */
    Quaternion &operator*=(const Quaternion &other);
  
    /**
     * In-place multiplication of a quaternion by a scalar
     */
    Quaternion &operator*=(double s)
    {
      for(size_t i = 0; i < 4; ++i)
        elems[i] *= s;
      return *this;
    }
  
    /**
     * Return the quaternion corresponding to the inverse transform
     */
    Quaternion inverse() const
    {
      double n = mgx::normsq(*this);
      return Quaternion(-x() / n, -y() / n, -z() / n, w() / n);
    }
  
    /**
     * Return the conjugate of the current quaternion
     */
    Quaternion conjugate() const {
      return Quaternion(-x(), -y(), -z(), w());
    }
  
    /**
     * Division of a quaternion by a real number
     */
    Quaternion &operator/=(double v)
    {
      for(size_t i = 0; i < 4; ++i)
        elems[i] /= v;
      return *this;
    }
  
    /**
     * Division of a quaternion by a real number
     */
    Quaternion operator/(double v) const
    {
      Quaternion tmp(*this);
      tmp /= v;
      return tmp;
    }
  
    /**
     * Fill the matrix as argument from the quaternion
     *
     * Multiplying with this matrix is equivalent to performing a rotation with
     * this quaternion.
     */
    void setMatrix(Matrix3d &m) const;
    /**
     * Fill the matrix as argument from the quaternion
     *
     * Multiplying with this matrix is equivalent to performing a rotation with
     * this quaternion.
     */
    void setMatrix(Matrix4d &m) const;
  
    /**
     * Returns the axis of the rotation corresponding to this quaternion
     */
    Point3d axis() const;
  
    /**
     * Returns the angle of the rotation corresponding to this quaternion
     */
    double angle() const;
  
    /**
     * Apply the rotation contained in this quaternion on the vector
     */
    Point3d rotate(const Point3d &v) const;
  
    /**
     * Apply the inverse of the rotation contained in this quaternion on the
     * vector.
     *
     * It is identical to calling \c this->inverse().rotate()
     */
    Point3d inverseRotate(const Point3d &v) const;
  };
  
  /**
   * Perform spherical linear interpolation
   *
   * \relates Quaternion
   */
  mgx_EXPORT Quaternion slerp(const Quaternion &q1, const Quaternion &q2, double t);
  
  inline Quaternion operator*(double s, const Quaternion& q)
  {
    Quaternion tmp(q);
    tmp *= s;
    return tmp;
  }
  
  inline Quaternion operator*(const Quaternion &q, double s)
  {
    Quaternion tmp(q);
    tmp *= s;
    return tmp;
  }
}

namespace std 
{
  /**
   * Power of a quaternion
   *
   * \relates Quaternion
   */
  mgx_EXPORT mgx::Quaternion pow(const mgx::Quaternion &q, double p);
}

#endif
