Add option for different quaternion storage layout

Describe the feature you would like to be implemented.

The current implementation of Eigen's quaternion class is constructed using (w, x, y, z) but internally stores the four values as (x, y, z, w).

  EIGEN_DEVICE_FUNC inline Quaternion(const Scalar& w, const Scalar& x, const Scalar& y, const Scalar& z)
      : m_coeffs(x, y, z, w) {}

According to the initial quaternion commit 366971be, this seems to be a rather arbitrary decision but no relevant design issue.

As noted in Section 3 of Joan Sol`a's notes this is not the most common choice. Further, it seems to be very confusing to users that constructing order and storage differs. In fact, when computing Jacobians of expressions involving quaternions this is very annoying since the expressions are not that easy to derive, but more seriously, differ from the formulas found in literature due to the different ordering.

That said, I would propose to add an option to the quaternion class that allows to decide which storage order should be used.

Any hints on how to implement the requested feature?

My suggestion is to add an option like

enum QuaternionStorageOrder {
  /** Storage order (x, y, z, w) */
  RealLast = 0,
  /** Storage order (w, x, y, z) */
  RealFirst = 0x1
};

template <typename Scalar_, int Options_, int StorageOrder_=RealLast>
class Quaternion : public QuaternionBase<Quaternion<Scalar_, Options_> > {
  // ...
}

that specifies the storage ordering of the data. The default case remains the current behavior so no breaking changes are introduced, but users have the possibility to use the other ordering.

At least to my knowledge such an option requires the lines below to be modified

  /** \returns the \c x coefficient */
  EIGEN_DEVICE_FUNC constexpr CoeffReturnType x() const { return this->derived().coeffs().coeff(0); }
  /** \returns the \c y coefficient */
  EIGEN_DEVICE_FUNC constexpr CoeffReturnType y() const { return this->derived().coeffs().coeff(1); }
  /** \returns the \c z coefficient */
  EIGEN_DEVICE_FUNC constexpr CoeffReturnType z() const { return this->derived().coeffs().coeff(2); }
  /** \returns the \c w coefficient */
  EIGEN_DEVICE_FUNC constexpr CoeffReturnType w() const { return this->derived().coeffs().coeff(3); }

  /** \returns a read-only vector expression of the imaginary part (x,y,z) */
  EIGEN_DEVICE_FUNC inline const VectorBlock<const Coefficients, 3> vec() const { return coeffs().template head<3>(); }

  /** \returns a vector expression of the imaginary part (x,y,z) */
  EIGEN_DEVICE_FUNC inline VectorBlock<Coefficients, 3> vec() { return coeffs().template head<3>(); }

  /** \sa QuaternionBase::Identity(), MatrixBase::setIdentity()
   */
  EIGEN_DEVICE_FUNC inline QuaternionBase& setIdentity() {
    coeffs() << Scalar(0), Scalar(0), Scalar(0), Scalar(1);
    return *this;
  }

  /** Constructs and initializes the quaternion \f$ w+xi+yj+zk \f$ from
   * its four coefficients \a w, \a x, \a y and \a z.
   *
   * \warning Note the order of the arguments: the real \a w coefficient first,
   * while internally the coefficients are stored in the following order:
   * [\c x, \c y, \c z, \c w]
   */
  EIGEN_DEVICE_FUNC inline Quaternion(const Scalar& w, const Scalar& x, const Scalar& y, const Scalar& z)
      : m_coeffs(x, y, z, w) {}

  /** Constructs and initializes a quaternion from its real part as a scalar,
   *  and its imaginary part as a 3-vector [\c x, \c y, \c z]
   */
  template <typename Derived>
  EIGEN_DEVICE_FUNC inline Quaternion(const Scalar& w, const Eigen::MatrixBase<Derived>& vec)
      : m_coeffs(vec.x(), vec.y(), vec.z(), w) {
    EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Derived, 3);
  }

All other parts of the file internally use vec.x(), etc. and do not require any changes.

Additional resources

If such a feature is desired and someone can give me a hint how to add this option to the class, I would be glad to create a PR soon.