diff options
Diffstat (limited to 'eigen/Eigen/src/Eigen2Support/Geometry')
| -rw-r--r-- | eigen/Eigen/src/Eigen2Support/Geometry/AlignedBox.h | 159 | ||||
| -rw-r--r-- | eigen/Eigen/src/Eigen2Support/Geometry/All.h | 115 | ||||
| -rw-r--r-- | eigen/Eigen/src/Eigen2Support/Geometry/AngleAxis.h | 214 | ||||
| -rw-r--r-- | eigen/Eigen/src/Eigen2Support/Geometry/CMakeLists.txt | 6 | ||||
| -rw-r--r-- | eigen/Eigen/src/Eigen2Support/Geometry/Hyperplane.h | 254 | ||||
| -rw-r--r-- | eigen/Eigen/src/Eigen2Support/Geometry/ParametrizedLine.h | 141 | ||||
| -rw-r--r-- | eigen/Eigen/src/Eigen2Support/Geometry/Quaternion.h | 495 | ||||
| -rw-r--r-- | eigen/Eigen/src/Eigen2Support/Geometry/Rotation2D.h | 145 | ||||
| -rw-r--r-- | eigen/Eigen/src/Eigen2Support/Geometry/RotationBase.h | 123 | ||||
| -rw-r--r-- | eigen/Eigen/src/Eigen2Support/Geometry/Scaling.h | 167 | ||||
| -rw-r--r-- | eigen/Eigen/src/Eigen2Support/Geometry/Transform.h | 786 | ||||
| -rw-r--r-- | eigen/Eigen/src/Eigen2Support/Geometry/Translation.h | 184 |
12 files changed, 2789 insertions, 0 deletions
diff --git a/eigen/Eigen/src/Eigen2Support/Geometry/AlignedBox.h b/eigen/Eigen/src/Eigen2Support/Geometry/AlignedBox.h new file mode 100644 index 0000000..2e4309d --- /dev/null +++ b/eigen/Eigen/src/Eigen2Support/Geometry/AlignedBox.h @@ -0,0 +1,159 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +// no include guard, we'll include this twice from All.h from Eigen2Support, and it's internal anyway + +namespace Eigen { + +/** \geometry_module \ingroup Geometry_Module + * \nonstableyet + * + * \class AlignedBox + * + * \brief An axis aligned box + * + * \param _Scalar the type of the scalar coefficients + * \param _AmbientDim the dimension of the ambient space, can be a compile time value or Dynamic. + * + * This class represents an axis aligned box as a pair of the minimal and maximal corners. + */ +template <typename _Scalar, int _AmbientDim> +class AlignedBox +{ +public: +EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim==Dynamic ? Dynamic : _AmbientDim+1) + enum { AmbientDimAtCompileTime = _AmbientDim }; + typedef _Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar,AmbientDimAtCompileTime,1> VectorType; + + /** Default constructor initializing a null box. */ + inline AlignedBox() + { if (AmbientDimAtCompileTime!=Dynamic) setNull(); } + + /** Constructs a null box with \a _dim the dimension of the ambient space. */ + inline explicit AlignedBox(int _dim) : m_min(_dim), m_max(_dim) + { setNull(); } + + /** Constructs a box with extremities \a _min and \a _max. */ + inline AlignedBox(const VectorType& _min, const VectorType& _max) : m_min(_min), m_max(_max) {} + + /** Constructs a box containing a single point \a p. */ + inline explicit AlignedBox(const VectorType& p) : m_min(p), m_max(p) {} + + ~AlignedBox() {} + + /** \returns the dimension in which the box holds */ + inline int dim() const { return AmbientDimAtCompileTime==Dynamic ? m_min.size()-1 : AmbientDimAtCompileTime; } + + /** \returns true if the box is null, i.e, empty. */ + inline bool isNull() const { return (m_min.cwise() > m_max).any(); } + + /** Makes \c *this a null/empty box. */ + inline void setNull() + { + m_min.setConstant( (std::numeric_limits<Scalar>::max)()); + m_max.setConstant(-(std::numeric_limits<Scalar>::max)()); + } + + /** \returns the minimal corner */ + inline const VectorType& (min)() const { return m_min; } + /** \returns a non const reference to the minimal corner */ + inline VectorType& (min)() { return m_min; } + /** \returns the maximal corner */ + inline const VectorType& (max)() const { return m_max; } + /** \returns a non const reference to the maximal corner */ + inline VectorType& (max)() { return m_max; } + + /** \returns true if the point \a p is inside the box \c *this. */ + inline bool contains(const VectorType& p) const + { return (m_min.cwise()<=p).all() && (p.cwise()<=m_max).all(); } + + /** \returns true if the box \a b is entirely inside the box \c *this. */ + inline bool contains(const AlignedBox& b) const + { return (m_min.cwise()<=(b.min)()).all() && ((b.max)().cwise()<=m_max).all(); } + + /** Extends \c *this such that it contains the point \a p and returns a reference to \c *this. */ + inline AlignedBox& extend(const VectorType& p) + { m_min = (m_min.cwise().min)(p); m_max = (m_max.cwise().max)(p); return *this; } + + /** Extends \c *this such that it contains the box \a b and returns a reference to \c *this. */ + inline AlignedBox& extend(const AlignedBox& b) + { m_min = (m_min.cwise().min)(b.m_min); m_max = (m_max.cwise().max)(b.m_max); return *this; } + + /** Clamps \c *this by the box \a b and returns a reference to \c *this. */ + inline AlignedBox& clamp(const AlignedBox& b) + { m_min = (m_min.cwise().max)(b.m_min); m_max = (m_max.cwise().min)(b.m_max); return *this; } + + /** Translate \c *this by the vector \a t and returns a reference to \c *this. */ + inline AlignedBox& translate(const VectorType& t) + { m_min += t; m_max += t; return *this; } + + /** \returns the squared distance between the point \a p and the box \c *this, + * and zero if \a p is inside the box. + * \sa exteriorDistance() + */ + inline Scalar squaredExteriorDistance(const VectorType& p) const; + + /** \returns the distance between the point \a p and the box \c *this, + * and zero if \a p is inside the box. + * \sa squaredExteriorDistance() + */ + inline Scalar exteriorDistance(const VectorType& p) const + { return ei_sqrt(squaredExteriorDistance(p)); } + + /** \returns \c *this with scalar type casted to \a NewScalarType + * + * Note that if \a NewScalarType is equal to the current scalar type of \c *this + * then this function smartly returns a const reference to \c *this. + */ + template<typename NewScalarType> + inline typename internal::cast_return_type<AlignedBox, + AlignedBox<NewScalarType,AmbientDimAtCompileTime> >::type cast() const + { + return typename internal::cast_return_type<AlignedBox, + AlignedBox<NewScalarType,AmbientDimAtCompileTime> >::type(*this); + } + + /** Copy constructor with scalar type conversion */ + template<typename OtherScalarType> + inline explicit AlignedBox(const AlignedBox<OtherScalarType,AmbientDimAtCompileTime>& other) + { + m_min = (other.min)().template cast<Scalar>(); + m_max = (other.max)().template cast<Scalar>(); + } + + /** \returns \c true if \c *this is approximately equal to \a other, within the precision + * determined by \a prec. + * + * \sa MatrixBase::isApprox() */ + bool isApprox(const AlignedBox& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const + { return m_min.isApprox(other.m_min, prec) && m_max.isApprox(other.m_max, prec); } + +protected: + + VectorType m_min, m_max; +}; + +template<typename Scalar,int AmbiantDim> +inline Scalar AlignedBox<Scalar,AmbiantDim>::squaredExteriorDistance(const VectorType& p) const +{ + Scalar dist2(0); + Scalar aux; + for (int k=0; k<dim(); ++k) + { + if ((aux = (p[k]-m_min[k]))<Scalar(0)) + dist2 += aux*aux; + else if ( (aux = (m_max[k]-p[k]))<Scalar(0)) + dist2 += aux*aux; + } + return dist2; +} + +} // end namespace Eigen diff --git a/eigen/Eigen/src/Eigen2Support/Geometry/All.h b/eigen/Eigen/src/Eigen2Support/Geometry/All.h new file mode 100644 index 0000000..e0b00fc --- /dev/null +++ b/eigen/Eigen/src/Eigen2Support/Geometry/All.h @@ -0,0 +1,115 @@ +#ifndef EIGEN2_GEOMETRY_MODULE_H +#define EIGEN2_GEOMETRY_MODULE_H + +#include <limits> + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +#if EIGEN2_SUPPORT_STAGE < STAGE20_RESOLVE_API_CONFLICTS +#include "RotationBase.h" +#include "Rotation2D.h" +#include "Quaternion.h" +#include "AngleAxis.h" +#include "Transform.h" +#include "Translation.h" +#include "Scaling.h" +#include "AlignedBox.h" +#include "Hyperplane.h" +#include "ParametrizedLine.h" +#endif + + +#define RotationBase eigen2_RotationBase +#define Rotation2D eigen2_Rotation2D +#define Rotation2Df eigen2_Rotation2Df +#define Rotation2Dd eigen2_Rotation2Dd + +#define Quaternion eigen2_Quaternion +#define Quaternionf eigen2_Quaternionf +#define Quaterniond eigen2_Quaterniond + +#define AngleAxis eigen2_AngleAxis +#define AngleAxisf eigen2_AngleAxisf +#define AngleAxisd eigen2_AngleAxisd + +#define Transform eigen2_Transform +#define Transform2f eigen2_Transform2f +#define Transform2d eigen2_Transform2d +#define Transform3f eigen2_Transform3f +#define Transform3d eigen2_Transform3d + +#define Translation eigen2_Translation +#define Translation2f eigen2_Translation2f +#define Translation2d eigen2_Translation2d +#define Translation3f eigen2_Translation3f +#define Translation3d eigen2_Translation3d + +#define Scaling eigen2_Scaling +#define Scaling2f eigen2_Scaling2f +#define Scaling2d eigen2_Scaling2d +#define Scaling3f eigen2_Scaling3f +#define Scaling3d eigen2_Scaling3d + +#define AlignedBox eigen2_AlignedBox + +#define Hyperplane eigen2_Hyperplane +#define ParametrizedLine eigen2_ParametrizedLine + +#define ei_toRotationMatrix eigen2_ei_toRotationMatrix +#define ei_quaternion_assign_impl eigen2_ei_quaternion_assign_impl +#define ei_transform_product_impl eigen2_ei_transform_product_impl + +#include "RotationBase.h" +#include "Rotation2D.h" +#include "Quaternion.h" +#include "AngleAxis.h" +#include "Transform.h" +#include "Translation.h" +#include "Scaling.h" +#include "AlignedBox.h" +#include "Hyperplane.h" +#include "ParametrizedLine.h" + +#undef ei_toRotationMatrix +#undef ei_quaternion_assign_impl +#undef ei_transform_product_impl + +#undef RotationBase +#undef Rotation2D +#undef Rotation2Df +#undef Rotation2Dd + +#undef Quaternion +#undef Quaternionf +#undef Quaterniond + +#undef AngleAxis +#undef AngleAxisf +#undef AngleAxisd + +#undef Transform +#undef Transform2f +#undef Transform2d +#undef Transform3f +#undef Transform3d + +#undef Translation +#undef Translation2f +#undef Translation2d +#undef Translation3f +#undef Translation3d + +#undef Scaling +#undef Scaling2f +#undef Scaling2d +#undef Scaling3f +#undef Scaling3d + +#undef AlignedBox + +#undef Hyperplane +#undef ParametrizedLine + +#endif // EIGEN2_GEOMETRY_MODULE_H diff --git a/eigen/Eigen/src/Eigen2Support/Geometry/AngleAxis.h b/eigen/Eigen/src/Eigen2Support/Geometry/AngleAxis.h new file mode 100644 index 0000000..af598a4 --- /dev/null +++ b/eigen/Eigen/src/Eigen2Support/Geometry/AngleAxis.h @@ -0,0 +1,214 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +// no include guard, we'll include this twice from All.h from Eigen2Support, and it's internal anyway + +namespace Eigen { + +/** \geometry_module \ingroup Geometry_Module + * + * \class AngleAxis + * + * \brief Represents a 3D rotation as a rotation angle around an arbitrary 3D axis + * + * \param _Scalar the scalar type, i.e., the type of the coefficients. + * + * The following two typedefs are provided for convenience: + * \li \c AngleAxisf for \c float + * \li \c AngleAxisd for \c double + * + * \addexample AngleAxisForEuler \label How to define a rotation from Euler-angles + * + * Combined with MatrixBase::Unit{X,Y,Z}, AngleAxis can be used to easily + * mimic Euler-angles. Here is an example: + * \include AngleAxis_mimic_euler.cpp + * Output: \verbinclude AngleAxis_mimic_euler.out + * + * \note This class is not aimed to be used to store a rotation transformation, + * but rather to make easier the creation of other rotation (Quaternion, rotation Matrix) + * and transformation objects. + * + * \sa class Quaternion, class Transform, MatrixBase::UnitX() + */ + +template<typename _Scalar> struct ei_traits<AngleAxis<_Scalar> > +{ + typedef _Scalar Scalar; +}; + +template<typename _Scalar> +class AngleAxis : public RotationBase<AngleAxis<_Scalar>,3> +{ + typedef RotationBase<AngleAxis<_Scalar>,3> Base; + +public: + + using Base::operator*; + + enum { Dim = 3 }; + /** the scalar type of the coefficients */ + typedef _Scalar Scalar; + typedef Matrix<Scalar,3,3> Matrix3; + typedef Matrix<Scalar,3,1> Vector3; + typedef Quaternion<Scalar> QuaternionType; + +protected: + + Vector3 m_axis; + Scalar m_angle; + +public: + + /** Default constructor without initialization. */ + AngleAxis() {} + /** Constructs and initialize the angle-axis rotation from an \a angle in radian + * and an \a axis which must be normalized. */ + template<typename Derived> + inline AngleAxis(Scalar angle, const MatrixBase<Derived>& axis) : m_axis(axis), m_angle(angle) {} + /** Constructs and initialize the angle-axis rotation from a quaternion \a q. */ + inline AngleAxis(const QuaternionType& q) { *this = q; } + /** Constructs and initialize the angle-axis rotation from a 3x3 rotation matrix. */ + template<typename Derived> + inline explicit AngleAxis(const MatrixBase<Derived>& m) { *this = m; } + + Scalar angle() const { return m_angle; } + Scalar& angle() { return m_angle; } + + const Vector3& axis() const { return m_axis; } + Vector3& axis() { return m_axis; } + + /** Concatenates two rotations */ + inline QuaternionType operator* (const AngleAxis& other) const + { return QuaternionType(*this) * QuaternionType(other); } + + /** Concatenates two rotations */ + inline QuaternionType operator* (const QuaternionType& other) const + { return QuaternionType(*this) * other; } + + /** Concatenates two rotations */ + friend inline QuaternionType operator* (const QuaternionType& a, const AngleAxis& b) + { return a * QuaternionType(b); } + + /** Concatenates two rotations */ + inline Matrix3 operator* (const Matrix3& other) const + { return toRotationMatrix() * other; } + + /** Concatenates two rotations */ + inline friend Matrix3 operator* (const Matrix3& a, const AngleAxis& b) + { return a * b.toRotationMatrix(); } + + /** Applies rotation to vector */ + inline Vector3 operator* (const Vector3& other) const + { return toRotationMatrix() * other; } + + /** \returns the inverse rotation, i.e., an angle-axis with opposite rotation angle */ + AngleAxis inverse() const + { return AngleAxis(-m_angle, m_axis); } + + AngleAxis& operator=(const QuaternionType& q); + template<typename Derived> + AngleAxis& operator=(const MatrixBase<Derived>& m); + + template<typename Derived> + AngleAxis& fromRotationMatrix(const MatrixBase<Derived>& m); + Matrix3 toRotationMatrix(void) const; + + /** \returns \c *this with scalar type casted to \a NewScalarType + * + * Note that if \a NewScalarType is equal to the current scalar type of \c *this + * then this function smartly returns a const reference to \c *this. + */ + template<typename NewScalarType> + inline typename internal::cast_return_type<AngleAxis,AngleAxis<NewScalarType> >::type cast() const + { return typename internal::cast_return_type<AngleAxis,AngleAxis<NewScalarType> >::type(*this); } + + /** Copy constructor with scalar type conversion */ + template<typename OtherScalarType> + inline explicit AngleAxis(const AngleAxis<OtherScalarType>& other) + { + m_axis = other.axis().template cast<Scalar>(); + m_angle = Scalar(other.angle()); + } + + /** \returns \c true if \c *this is approximately equal to \a other, within the precision + * determined by \a prec. + * + * \sa MatrixBase::isApprox() */ + bool isApprox(const AngleAxis& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const + { return m_axis.isApprox(other.m_axis, prec) && ei_isApprox(m_angle,other.m_angle, prec); } +}; + +/** \ingroup Geometry_Module + * single precision angle-axis type */ +typedef AngleAxis<float> AngleAxisf; +/** \ingroup Geometry_Module + * double precision angle-axis type */ +typedef AngleAxis<double> AngleAxisd; + +/** Set \c *this from a quaternion. + * The axis is normalized. + */ +template<typename Scalar> +AngleAxis<Scalar>& AngleAxis<Scalar>::operator=(const QuaternionType& q) +{ + Scalar n2 = q.vec().squaredNorm(); + if (n2 < precision<Scalar>()*precision<Scalar>()) + { + m_angle = 0; + m_axis << 1, 0, 0; + } + else + { + m_angle = 2*std::acos(q.w()); + m_axis = q.vec() / ei_sqrt(n2); + } + return *this; +} + +/** Set \c *this from a 3x3 rotation matrix \a mat. + */ +template<typename Scalar> +template<typename Derived> +AngleAxis<Scalar>& AngleAxis<Scalar>::operator=(const MatrixBase<Derived>& mat) +{ + // Since a direct conversion would not be really faster, + // let's use the robust Quaternion implementation: + return *this = QuaternionType(mat); +} + +/** Constructs and \returns an equivalent 3x3 rotation matrix. + */ +template<typename Scalar> +typename AngleAxis<Scalar>::Matrix3 +AngleAxis<Scalar>::toRotationMatrix(void) const +{ + Matrix3 res; + Vector3 sin_axis = ei_sin(m_angle) * m_axis; + Scalar c = ei_cos(m_angle); + Vector3 cos1_axis = (Scalar(1)-c) * m_axis; + + Scalar tmp; + tmp = cos1_axis.x() * m_axis.y(); + res.coeffRef(0,1) = tmp - sin_axis.z(); + res.coeffRef(1,0) = tmp + sin_axis.z(); + + tmp = cos1_axis.x() * m_axis.z(); + res.coeffRef(0,2) = tmp + sin_axis.y(); + res.coeffRef(2,0) = tmp - sin_axis.y(); + + tmp = cos1_axis.y() * m_axis.z(); + res.coeffRef(1,2) = tmp - sin_axis.x(); + res.coeffRef(2,1) = tmp + sin_axis.x(); + + res.diagonal() = (cos1_axis.cwise() * m_axis).cwise() + c; + + return res; +} + +} // end namespace Eigen diff --git a/eigen/Eigen/src/Eigen2Support/Geometry/CMakeLists.txt b/eigen/Eigen/src/Eigen2Support/Geometry/CMakeLists.txt new file mode 100644 index 0000000..c347a8f --- /dev/null +++ b/eigen/Eigen/src/Eigen2Support/Geometry/CMakeLists.txt @@ -0,0 +1,6 @@ +FILE(GLOB Eigen_Eigen2Support_Geometry_SRCS "*.h") + +INSTALL(FILES + ${Eigen_Eigen2Support_Geometry_SRCS} + DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Eigen2Support/Geometry + ) diff --git a/eigen/Eigen/src/Eigen2Support/Geometry/Hyperplane.h b/eigen/Eigen/src/Eigen2Support/Geometry/Hyperplane.h new file mode 100644 index 0000000..b95bf00 --- /dev/null +++ b/eigen/Eigen/src/Eigen2Support/Geometry/Hyperplane.h @@ -0,0 +1,254 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +// no include guard, we'll include this twice from All.h from Eigen2Support, and it's internal anyway + +namespace Eigen { + +/** \geometry_module \ingroup Geometry_Module + * + * \class Hyperplane + * + * \brief A hyperplane + * + * A hyperplane is an affine subspace of dimension n-1 in a space of dimension n. + * For example, a hyperplane in a plane is a line; a hyperplane in 3-space is a plane. + * + * \param _Scalar the scalar type, i.e., the type of the coefficients + * \param _AmbientDim the dimension of the ambient space, can be a compile time value or Dynamic. + * Notice that the dimension of the hyperplane is _AmbientDim-1. + * + * This class represents an hyperplane as the zero set of the implicit equation + * \f$ n \cdot x + d = 0 \f$ where \f$ n \f$ is a unit normal vector of the plane (linear part) + * and \f$ d \f$ is the distance (offset) to the origin. + */ +template <typename _Scalar, int _AmbientDim> +class Hyperplane +{ +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim==Dynamic ? Dynamic : _AmbientDim+1) + enum { AmbientDimAtCompileTime = _AmbientDim }; + typedef _Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar,AmbientDimAtCompileTime,1> VectorType; + typedef Matrix<Scalar,int(AmbientDimAtCompileTime)==Dynamic + ? Dynamic + : int(AmbientDimAtCompileTime)+1,1> Coefficients; + typedef Block<Coefficients,AmbientDimAtCompileTime,1> NormalReturnType; + + /** Default constructor without initialization */ + inline Hyperplane() {} + + /** Constructs a dynamic-size hyperplane with \a _dim the dimension + * of the ambient space */ + inline explicit Hyperplane(int _dim) : m_coeffs(_dim+1) {} + + /** Construct a plane from its normal \a n and a point \a e onto the plane. + * \warning the vector normal is assumed to be normalized. + */ + inline Hyperplane(const VectorType& n, const VectorType& e) + : m_coeffs(n.size()+1) + { + normal() = n; + offset() = -e.eigen2_dot(n); + } + + /** Constructs a plane from its normal \a n and distance to the origin \a d + * such that the algebraic equation of the plane is \f$ n \cdot x + d = 0 \f$. + * \warning the vector normal is assumed to be normalized. + */ + inline Hyperplane(const VectorType& n, Scalar d) + : m_coeffs(n.size()+1) + { + normal() = n; + offset() = d; + } + + /** Constructs a hyperplane passing through the two points. If the dimension of the ambient space + * is greater than 2, then there isn't uniqueness, so an arbitrary choice is made. + */ + static inline Hyperplane Through(const VectorType& p0, const VectorType& p1) + { + Hyperplane result(p0.size()); + result.normal() = (p1 - p0).unitOrthogonal(); + result.offset() = -result.normal().eigen2_dot(p0); + return result; + } + + /** Constructs a hyperplane passing through the three points. The dimension of the ambient space + * is required to be exactly 3. + */ + static inline Hyperplane Through(const VectorType& p0, const VectorType& p1, const VectorType& p2) + { + EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(VectorType, 3) + Hyperplane result(p0.size()); + result.normal() = (p2 - p0).cross(p1 - p0).normalized(); + result.offset() = -result.normal().eigen2_dot(p0); + return result; + } + + /** Constructs a hyperplane passing through the parametrized line \a parametrized. + * If the dimension of the ambient space is greater than 2, then there isn't uniqueness, + * so an arbitrary choice is made. + */ + // FIXME to be consitent with the rest this could be implemented as a static Through function ?? + explicit Hyperplane(const ParametrizedLine<Scalar, AmbientDimAtCompileTime>& parametrized) + { + normal() = parametrized.direction().unitOrthogonal(); + offset() = -normal().eigen2_dot(parametrized.origin()); + } + + ~Hyperplane() {} + + /** \returns the dimension in which the plane holds */ + inline int dim() const { return int(AmbientDimAtCompileTime)==Dynamic ? m_coeffs.size()-1 : int(AmbientDimAtCompileTime); } + + /** normalizes \c *this */ + void normalize(void) + { + m_coeffs /= normal().norm(); + } + + /** \returns the signed distance between the plane \c *this and a point \a p. + * \sa absDistance() + */ + inline Scalar signedDistance(const VectorType& p) const { return p.eigen2_dot(normal()) + offset(); } + + /** \returns the absolute distance between the plane \c *this and a point \a p. + * \sa signedDistance() + */ + inline Scalar absDistance(const VectorType& p) const { return ei_abs(signedDistance(p)); } + + /** \returns the projection of a point \a p onto the plane \c *this. + */ + inline VectorType projection(const VectorType& p) const { return p - signedDistance(p) * normal(); } + + /** \returns a constant reference to the unit normal vector of the plane, which corresponds + * to the linear part of the implicit equation. + */ + inline const NormalReturnType normal() const { return NormalReturnType(*const_cast<Coefficients*>(&m_coeffs),0,0,dim(),1); } + + /** \returns a non-constant reference to the unit normal vector of the plane, which corresponds + * to the linear part of the implicit equation. + */ + inline NormalReturnType normal() { return NormalReturnType(m_coeffs,0,0,dim(),1); } + + /** \returns the distance to the origin, which is also the "constant term" of the implicit equation + * \warning the vector normal is assumed to be normalized. + */ + inline const Scalar& offset() const { return m_coeffs.coeff(dim()); } + + /** \returns a non-constant reference to the distance to the origin, which is also the constant part + * of the implicit equation */ + inline Scalar& offset() { return m_coeffs(dim()); } + + /** \returns a constant reference to the coefficients c_i of the plane equation: + * \f$ c_0*x_0 + ... + c_{d-1}*x_{d-1} + c_d = 0 \f$ + */ + inline const Coefficients& coeffs() const { return m_coeffs; } + + /** \returns a non-constant reference to the coefficients c_i of the plane equation: + * \f$ c_0*x_0 + ... + c_{d-1}*x_{d-1} + c_d = 0 \f$ + */ + inline Coefficients& coeffs() { return m_coeffs; } + + /** \returns the intersection of *this with \a other. + * + * \warning The ambient space must be a plane, i.e. have dimension 2, so that \c *this and \a other are lines. + * + * \note If \a other is approximately parallel to *this, this method will return any point on *this. + */ + VectorType intersection(const Hyperplane& other) + { + EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(VectorType, 2) + Scalar det = coeffs().coeff(0) * other.coeffs().coeff(1) - coeffs().coeff(1) * other.coeffs().coeff(0); + // since the line equations ax+by=c are normalized with a^2+b^2=1, the following tests + // whether the two lines are approximately parallel. + if(ei_isMuchSmallerThan(det, Scalar(1))) + { // special case where the two lines are approximately parallel. Pick any point on the first line. + if(ei_abs(coeffs().coeff(1))>ei_abs(coeffs().coeff(0))) + return VectorType(coeffs().coeff(1), -coeffs().coeff(2)/coeffs().coeff(1)-coeffs().coeff(0)); + else + return VectorType(-coeffs().coeff(2)/coeffs().coeff(0)-coeffs().coeff(1), coeffs().coeff(0)); + } + else + { // general case + Scalar invdet = Scalar(1) / det; + return VectorType(invdet*(coeffs().coeff(1)*other.coeffs().coeff(2)-other.coeffs().coeff(1)*coeffs().coeff(2)), + invdet*(other.coeffs().coeff(0)*coeffs().coeff(2)-coeffs().coeff(0)*other.coeffs().coeff(2))); + } + } + + /** Applies the transformation matrix \a mat to \c *this and returns a reference to \c *this. + * + * \param mat the Dim x Dim transformation matrix + * \param traits specifies whether the matrix \a mat represents an Isometry + * or a more generic Affine transformation. The default is Affine. + */ + template<typename XprType> + inline Hyperplane& transform(const MatrixBase<XprType>& mat, TransformTraits traits = Affine) + { + if (traits==Affine) + normal() = mat.inverse().transpose() * normal(); + else if (traits==Isometry) + normal() = mat * normal(); + else + { + ei_assert("invalid traits value in Hyperplane::transform()"); + } + return *this; + } + + /** Applies the transformation \a t to \c *this and returns a reference to \c *this. + * + * \param t the transformation of dimension Dim + * \param traits specifies whether the transformation \a t represents an Isometry + * or a more generic Affine transformation. The default is Affine. + * Other kind of transformations are not supported. + */ + inline Hyperplane& transform(const Transform<Scalar,AmbientDimAtCompileTime>& t, + TransformTraits traits = Affine) + { + transform(t.linear(), traits); + offset() -= t.translation().eigen2_dot(normal()); + return *this; + } + + /** \returns \c *this with scalar type casted to \a NewScalarType + * + * Note that if \a NewScalarType is equal to the current scalar type of \c *this + * then this function smartly returns a const reference to \c *this. + */ + template<typename NewScalarType> + inline typename internal::cast_return_type<Hyperplane, + Hyperplane<NewScalarType,AmbientDimAtCompileTime> >::type cast() const + { + return typename internal::cast_return_type<Hyperplane, + Hyperplane<NewScalarType,AmbientDimAtCompileTime> >::type(*this); + } + + /** Copy constructor with scalar type conversion */ + template<typename OtherScalarType> + inline explicit Hyperplane(const Hyperplane<OtherScalarType,AmbientDimAtCompileTime>& other) + { m_coeffs = other.coeffs().template cast<Scalar>(); } + + /** \returns \c true if \c *this is approximately equal to \a other, within the precision + * determined by \a prec. + * + * \sa MatrixBase::isApprox() */ + bool isApprox(const Hyperplane& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const + { return m_coeffs.isApprox(other.m_coeffs, prec); } + +protected: + + Coefficients m_coeffs; +}; + +} // end namespace Eigen diff --git a/eigen/Eigen/src/Eigen2Support/Geometry/ParametrizedLine.h b/eigen/Eigen/src/Eigen2Support/Geometry/ParametrizedLine.h new file mode 100644 index 0000000..9b57b7e --- /dev/null +++ b/eigen/Eigen/src/Eigen2Support/Geometry/ParametrizedLine.h @@ -0,0 +1,141 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +// no include guard, we'll include this twice from All.h from Eigen2Support, and it's internal anyway + +namespace Eigen { + +/** \geometry_module \ingroup Geometry_Module + * + * \class ParametrizedLine + * + * \brief A parametrized line + * + * A parametrized line is defined by an origin point \f$ \mathbf{o} \f$ and a unit + * direction vector \f$ \mathbf{d} \f$ such that the line corresponds to + * the set \f$ l(t) = \mathbf{o} + t \mathbf{d} \f$, \f$ l \in \mathbf{R} \f$. + * + * \param _Scalar the scalar type, i.e., the type of the coefficients + * \param _AmbientDim the dimension of the ambient space, can be a compile time value or Dynamic. + */ +template <typename _Scalar, int _AmbientDim> +class ParametrizedLine +{ +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim) + enum { AmbientDimAtCompileTime = _AmbientDim }; + typedef _Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar,AmbientDimAtCompileTime,1> VectorType; + + /** Default constructor without initialization */ + inline ParametrizedLine() {} + + /** Constructs a dynamic-size line with \a _dim the dimension + * of the ambient space */ + inline explicit ParametrizedLine(int _dim) : m_origin(_dim), m_direction(_dim) {} + + /** Initializes a parametrized line of direction \a direction and origin \a origin. + * \warning the vector direction is assumed to be normalized. + */ + ParametrizedLine(const VectorType& origin, const VectorType& direction) + : m_origin(origin), m_direction(direction) {} + + explicit ParametrizedLine(const Hyperplane<_Scalar, _AmbientDim>& hyperplane); + + /** Constructs a parametrized line going from \a p0 to \a p1. */ + static inline ParametrizedLine Through(const VectorType& p0, const VectorType& p1) + { return ParametrizedLine(p0, (p1-p0).normalized()); } + + ~ParametrizedLine() {} + + /** \returns the dimension in which the line holds */ + inline int dim() const { return m_direction.size(); } + + const VectorType& origin() const { return m_origin; } + VectorType& origin() { return m_origin; } + + const VectorType& direction() const { return m_direction; } + VectorType& direction() { return m_direction; } + + /** \returns the squared distance of a point \a p to its projection onto the line \c *this. + * \sa distance() + */ + RealScalar squaredDistance(const VectorType& p) const + { + VectorType diff = p-origin(); + return (diff - diff.eigen2_dot(direction())* direction()).squaredNorm(); + } + /** \returns the distance of a point \a p to its projection onto the line \c *this. + * \sa squaredDistance() + */ + RealScalar distance(const VectorType& p) const { return ei_sqrt(squaredDistance(p)); } + + /** \returns the projection of a point \a p onto the line \c *this. */ + VectorType projection(const VectorType& p) const + { return origin() + (p-origin()).eigen2_dot(direction()) * direction(); } + + Scalar intersection(const Hyperplane<_Scalar, _AmbientDim>& hyperplane); + + /** \returns \c *this with scalar type casted to \a NewScalarType + * + * Note that if \a NewScalarType is equal to the current scalar type of \c *this + * then this function smartly returns a const reference to \c *this. + */ + template<typename NewScalarType> + inline typename internal::cast_return_type<ParametrizedLine, + ParametrizedLine<NewScalarType,AmbientDimAtCompileTime> >::type cast() const + { + return typename internal::cast_return_type<ParametrizedLine, + ParametrizedLine<NewScalarType,AmbientDimAtCompileTime> >::type(*this); + } + + /** Copy constructor with scalar type conversion */ + template<typename OtherScalarType> + inline explicit ParametrizedLine(const ParametrizedLine<OtherScalarType,AmbientDimAtCompileTime>& other) + { + m_origin = other.origin().template cast<Scalar>(); + m_direction = other.direction().template cast<Scalar>(); + } + + /** \returns \c true if \c *this is approximately equal to \a other, within the precision + * determined by \a prec. + * + * \sa MatrixBase::isApprox() */ + bool isApprox(const ParametrizedLine& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const + { return m_origin.isApprox(other.m_origin, prec) && m_direction.isApprox(other.m_direction, prec); } + +protected: + + VectorType m_origin, m_direction; +}; + +/** Constructs a parametrized line from a 2D hyperplane + * + * \warning the ambient space must have dimension 2 such that the hyperplane actually describes a line + */ +template <typename _Scalar, int _AmbientDim> +inline ParametrizedLine<_Scalar, _AmbientDim>::ParametrizedLine(const Hyperplane<_Scalar, _AmbientDim>& hyperplane) +{ + EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(VectorType, 2) + direction() = hyperplane.normal().unitOrthogonal(); + origin() = -hyperplane.normal()*hyperplane.offset(); +} + +/** \returns the parameter value of the intersection between \c *this and the given hyperplane + */ +template <typename _Scalar, int _AmbientDim> +inline _Scalar ParametrizedLine<_Scalar, _AmbientDim>::intersection(const Hyperplane<_Scalar, _AmbientDim>& hyperplane) +{ + return -(hyperplane.offset()+origin().eigen2_dot(hyperplane.normal())) + /(direction().eigen2_dot(hyperplane.normal())); +} + +} // end namespace Eigen diff --git a/eigen/Eigen/src/Eigen2Support/Geometry/Quaternion.h b/eigen/Eigen/src/Eigen2Support/Geometry/Quaternion.h new file mode 100644 index 0000000..4b6390c --- /dev/null +++ b/eigen/Eigen/src/Eigen2Support/Geometry/Quaternion.h @@ -0,0 +1,495 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +// no include guard, we'll include this twice from All.h from Eigen2Support, and it's internal anyway + +namespace Eigen { + +template<typename Other, + int OtherRows=Other::RowsAtCompileTime, + int OtherCols=Other::ColsAtCompileTime> +struct ei_quaternion_assign_impl; + +/** \geometry_module \ingroup Geometry_Module + * + * \class Quaternion + * + * \brief The quaternion class used to represent 3D orientations and rotations + * + * \param _Scalar the scalar type, i.e., the type of the coefficients + * + * This class represents a quaternion \f$ w+xi+yj+zk \f$ that is a convenient representation of + * orientations and rotations of objects in three dimensions. Compared to other representations + * like Euler angles or 3x3 matrices, quatertions offer the following advantages: + * \li \b compact storage (4 scalars) + * \li \b efficient to compose (28 flops), + * \li \b stable spherical interpolation + * + * The following two typedefs are provided for convenience: + * \li \c Quaternionf for \c float + * \li \c Quaterniond for \c double + * + * \sa class AngleAxis, class Transform + */ + +template<typename _Scalar> struct ei_traits<Quaternion<_Scalar> > +{ + typedef _Scalar Scalar; +}; + +template<typename _Scalar> +class Quaternion : public RotationBase<Quaternion<_Scalar>,3> +{ + typedef RotationBase<Quaternion<_Scalar>,3> Base; + +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,4) + + using Base::operator*; + + /** the scalar type of the coefficients */ + typedef _Scalar Scalar; + + /** the type of the Coefficients 4-vector */ + typedef Matrix<Scalar, 4, 1> Coefficients; + /** the type of a 3D vector */ + typedef Matrix<Scalar,3,1> Vector3; + /** the equivalent rotation matrix type */ + typedef Matrix<Scalar,3,3> Matrix3; + /** the equivalent angle-axis type */ + typedef AngleAxis<Scalar> AngleAxisType; + + /** \returns the \c x coefficient */ + inline Scalar x() const { return m_coeffs.coeff(0); } + /** \returns the \c y coefficient */ + inline Scalar y() const { return m_coeffs.coeff(1); } + /** \returns the \c z coefficient */ + inline Scalar z() const { return m_coeffs.coeff(2); } + /** \returns the \c w coefficient */ + inline Scalar w() const { return m_coeffs.coeff(3); } + + /** \returns a reference to the \c x coefficient */ + inline Scalar& x() { return m_coeffs.coeffRef(0); } + /** \returns a reference to the \c y coefficient */ + inline Scalar& y() { return m_coeffs.coeffRef(1); } + /** \returns a reference to the \c z coefficient */ + inline Scalar& z() { return m_coeffs.coeffRef(2); } + /** \returns a reference to the \c w coefficient */ + inline Scalar& w() { return m_coeffs.coeffRef(3); } + + /** \returns a read-only vector expression of the imaginary part (x,y,z) */ + inline const Block<const Coefficients,3,1> vec() const { return m_coeffs.template start<3>(); } + + /** \returns a vector expression of the imaginary part (x,y,z) */ + inline Block<Coefficients,3,1> vec() { return m_coeffs.template start<3>(); } + + /** \returns a read-only vector expression of the coefficients (x,y,z,w) */ + inline const Coefficients& coeffs() const { return m_coeffs; } + + /** \returns a vector expression of the coefficients (x,y,z,w) */ + inline Coefficients& coeffs() { return m_coeffs; } + + /** Default constructor leaving the quaternion uninitialized. */ + inline Quaternion() {} + + /** 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] + */ + inline Quaternion(Scalar w, Scalar x, Scalar y, Scalar z) + { m_coeffs << x, y, z, w; } + + /** Copy constructor */ + inline Quaternion(const Quaternion& other) { m_coeffs = other.m_coeffs; } + + /** Constructs and initializes a quaternion from the angle-axis \a aa */ + explicit inline Quaternion(const AngleAxisType& aa) { *this = aa; } + + /** Constructs and initializes a quaternion from either: + * - a rotation matrix expression, + * - a 4D vector expression representing quaternion coefficients. + * \sa operator=(MatrixBase<Derived>) + */ + template<typename Derived> + explicit inline Quaternion(const MatrixBase<Derived>& other) { *this = other; } + + Quaternion& operator=(const Quaternion& other); + Quaternion& operator=(const AngleAxisType& aa); + template<typename Derived> + Quaternion& operator=(const MatrixBase<Derived>& m); + + /** \returns a quaternion representing an identity rotation + * \sa MatrixBase::Identity() + */ + static inline Quaternion Identity() { return Quaternion(1, 0, 0, 0); } + + /** \sa Quaternion::Identity(), MatrixBase::setIdentity() + */ + inline Quaternion& setIdentity() { m_coeffs << 0, 0, 0, 1; return *this; } + + /** \returns the squared norm of the quaternion's coefficients + * \sa Quaternion::norm(), MatrixBase::squaredNorm() + */ + inline Scalar squaredNorm() const { return m_coeffs.squaredNorm(); } + + /** \returns the norm of the quaternion's coefficients + * \sa Quaternion::squaredNorm(), MatrixBase::norm() + */ + inline Scalar norm() const { return m_coeffs.norm(); } + + /** Normalizes the quaternion \c *this + * \sa normalized(), MatrixBase::normalize() */ + inline void normalize() { m_coeffs.normalize(); } + /** \returns a normalized version of \c *this + * \sa normalize(), MatrixBase::normalized() */ + inline Quaternion normalized() const { return Quaternion(m_coeffs.normalized()); } + + /** \returns the dot product of \c *this and \a other + * Geometrically speaking, the dot product of two unit quaternions + * corresponds to the cosine of half the angle between the two rotations. + * \sa angularDistance() + */ + inline Scalar eigen2_dot(const Quaternion& other) const { return m_coeffs.eigen2_dot(other.m_coeffs); } + + inline Scalar angularDistance(const Quaternion& other) const; + + Matrix3 toRotationMatrix(void) const; + + template<typename Derived1, typename Derived2> + Quaternion& setFromTwoVectors(const MatrixBase<Derived1>& a, const MatrixBase<Derived2>& b); + + inline Quaternion operator* (const Quaternion& q) const; + inline Quaternion& operator*= (const Quaternion& q); + + Quaternion inverse(void) const; + Quaternion conjugate(void) const; + + Quaternion slerp(Scalar t, const Quaternion& other) const; + + template<typename Derived> + Vector3 operator* (const MatrixBase<Derived>& vec) const; + + /** \returns \c *this with scalar type casted to \a NewScalarType + * + * Note that if \a NewScalarType is equal to the current scalar type of \c *this + * then this function smartly returns a const reference to \c *this. + */ + template<typename NewScalarType> + inline typename internal::cast_return_type<Quaternion,Quaternion<NewScalarType> >::type cast() const + { return typename internal::cast_return_type<Quaternion,Quaternion<NewScalarType> >::type(*this); } + + /** Copy constructor with scalar type conversion */ + template<typename OtherScalarType> + inline explicit Quaternion(const Quaternion<OtherScalarType>& other) + { m_coeffs = other.coeffs().template cast<Scalar>(); } + + /** \returns \c true if \c *this is approximately equal to \a other, within the precision + * determined by \a prec. + * + * \sa MatrixBase::isApprox() */ + bool isApprox(const Quaternion& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const + { return m_coeffs.isApprox(other.m_coeffs, prec); } + +protected: + Coefficients m_coeffs; +}; + +/** \ingroup Geometry_Module + * single precision quaternion type */ +typedef Quaternion<float> Quaternionf; +/** \ingroup Geometry_Module + * double precision quaternion type */ +typedef Quaternion<double> Quaterniond; + +// Generic Quaternion * Quaternion product +template<typename Scalar> inline Quaternion<Scalar> +ei_quaternion_product(const Quaternion<Scalar>& a, const Quaternion<Scalar>& b) +{ + return Quaternion<Scalar> + ( + a.w() * b.w() - a.x() * b.x() - a.y() * b.y() - a.z() * b.z(), + a.w() * b.x() + a.x() * b.w() + a.y() * b.z() - a.z() * b.y(), + a.w() * b.y() + a.y() * b.w() + a.z() * b.x() - a.x() * b.z(), + a.w() * b.z() + a.z() * b.w() + a.x() * b.y() - a.y() * b.x() + ); +} + +/** \returns the concatenation of two rotations as a quaternion-quaternion product */ +template <typename Scalar> +inline Quaternion<Scalar> Quaternion<Scalar>::operator* (const Quaternion& other) const +{ + return ei_quaternion_product(*this,other); +} + +/** \sa operator*(Quaternion) */ +template <typename Scalar> +inline Quaternion<Scalar>& Quaternion<Scalar>::operator*= (const Quaternion& other) +{ + return (*this = *this * other); +} + +/** Rotation of a vector by a quaternion. + * \remarks If the quaternion is used to rotate several points (>1) + * then it is much more efficient to first convert it to a 3x3 Matrix. + * Comparison of the operation cost for n transformations: + * - Quaternion: 30n + * - Via a Matrix3: 24 + 15n + */ +template <typename Scalar> +template<typename Derived> +inline typename Quaternion<Scalar>::Vector3 +Quaternion<Scalar>::operator* (const MatrixBase<Derived>& v) const +{ + // Note that this algorithm comes from the optimization by hand + // of the conversion to a Matrix followed by a Matrix/Vector product. + // It appears to be much faster than the common algorithm found + // in the litterature (30 versus 39 flops). It also requires two + // Vector3 as temporaries. + Vector3 uv; + uv = 2 * this->vec().cross(v); + return v + this->w() * uv + this->vec().cross(uv); +} + +template<typename Scalar> +inline Quaternion<Scalar>& Quaternion<Scalar>::operator=(const Quaternion& other) +{ + m_coeffs = other.m_coeffs; + return *this; +} + +/** Set \c *this from an angle-axis \a aa and returns a reference to \c *this + */ +template<typename Scalar> +inline Quaternion<Scalar>& Quaternion<Scalar>::operator=(const AngleAxisType& aa) +{ + Scalar ha = Scalar(0.5)*aa.angle(); // Scalar(0.5) to suppress precision loss warnings + this->w() = ei_cos(ha); + this->vec() = ei_sin(ha) * aa.axis(); + return *this; +} + +/** Set \c *this from the expression \a xpr: + * - if \a xpr is a 4x1 vector, then \a xpr is assumed to be a quaternion + * - if \a xpr is a 3x3 matrix, then \a xpr is assumed to be rotation matrix + * and \a xpr is converted to a quaternion + */ +template<typename Scalar> +template<typename Derived> +inline Quaternion<Scalar>& Quaternion<Scalar>::operator=(const MatrixBase<Derived>& xpr) +{ + ei_quaternion_assign_impl<Derived>::run(*this, xpr.derived()); + return *this; +} + +/** Convert the quaternion to a 3x3 rotation matrix */ +template<typename Scalar> +inline typename Quaternion<Scalar>::Matrix3 +Quaternion<Scalar>::toRotationMatrix(void) const +{ + // NOTE if inlined, then gcc 4.2 and 4.4 get rid of the temporary (not gcc 4.3 !!) + // if not inlined then the cost of the return by value is huge ~ +35%, + // however, not inlining this function is an order of magnitude slower, so + // it has to be inlined, and so the return by value is not an issue + Matrix3 res; + + const Scalar tx = Scalar(2)*this->x(); + const Scalar ty = Scalar(2)*this->y(); + const Scalar tz = Scalar(2)*this->z(); + const Scalar twx = tx*this->w(); + const Scalar twy = ty*this->w(); + const Scalar twz = tz*this->w(); + const Scalar txx = tx*this->x(); + const Scalar txy = ty*this->x(); + const Scalar txz = tz*this->x(); + const Scalar tyy = ty*this->y(); + const Scalar tyz = tz*this->y(); + const Scalar tzz = tz*this->z(); + + res.coeffRef(0,0) = Scalar(1)-(tyy+tzz); + res.coeffRef(0,1) = txy-twz; + res.coeffRef(0,2) = txz+twy; + res.coeffRef(1,0) = txy+twz; + res.coeffRef(1,1) = Scalar(1)-(txx+tzz); + res.coeffRef(1,2) = tyz-twx; + res.coeffRef(2,0) = txz-twy; + res.coeffRef(2,1) = tyz+twx; + res.coeffRef(2,2) = Scalar(1)-(txx+tyy); + + return res; +} + +/** Sets *this to be a quaternion representing a rotation sending the vector \a a to the vector \a b. + * + * \returns a reference to *this. + * + * Note that the two input vectors do \b not have to be normalized. + */ +template<typename Scalar> +template<typename Derived1, typename Derived2> +inline Quaternion<Scalar>& Quaternion<Scalar>::setFromTwoVectors(const MatrixBase<Derived1>& a, const MatrixBase<Derived2>& b) +{ + Vector3 v0 = a.normalized(); + Vector3 v1 = b.normalized(); + Scalar c = v0.eigen2_dot(v1); + + // if dot == 1, vectors are the same + if (ei_isApprox(c,Scalar(1))) + { + // set to identity + this->w() = 1; this->vec().setZero(); + return *this; + } + // if dot == -1, vectors are opposites + if (ei_isApprox(c,Scalar(-1))) + { + this->vec() = v0.unitOrthogonal(); + this->w() = 0; + return *this; + } + + Vector3 axis = v0.cross(v1); + Scalar s = ei_sqrt((Scalar(1)+c)*Scalar(2)); + Scalar invs = Scalar(1)/s; + this->vec() = axis * invs; + this->w() = s * Scalar(0.5); + + return *this; +} + +/** \returns the multiplicative inverse of \c *this + * Note that in most cases, i.e., if you simply want the opposite rotation, + * and/or the quaternion is normalized, then it is enough to use the conjugate. + * + * \sa Quaternion::conjugate() + */ +template <typename Scalar> +inline Quaternion<Scalar> Quaternion<Scalar>::inverse() const +{ + // FIXME should this function be called multiplicativeInverse and conjugate() be called inverse() or opposite() ?? + Scalar n2 = this->squaredNorm(); + if (n2 > 0) + return Quaternion(conjugate().coeffs() / n2); + else + { + // return an invalid result to flag the error + return Quaternion(Coefficients::Zero()); + } +} + +/** \returns the conjugate of the \c *this which is equal to the multiplicative inverse + * if the quaternion is normalized. + * The conjugate of a quaternion represents the opposite rotation. + * + * \sa Quaternion::inverse() + */ +template <typename Scalar> +inline Quaternion<Scalar> Quaternion<Scalar>::conjugate() const +{ + return Quaternion(this->w(),-this->x(),-this->y(),-this->z()); +} + +/** \returns the angle (in radian) between two rotations + * \sa eigen2_dot() + */ +template <typename Scalar> +inline Scalar Quaternion<Scalar>::angularDistance(const Quaternion& other) const +{ + double d = ei_abs(this->eigen2_dot(other)); + if (d>=1.0) + return 0; + return Scalar(2) * std::acos(d); +} + +/** \returns the spherical linear interpolation between the two quaternions + * \c *this and \a other at the parameter \a t + */ +template <typename Scalar> +Quaternion<Scalar> Quaternion<Scalar>::slerp(Scalar t, const Quaternion& other) const +{ + static const Scalar one = Scalar(1) - machine_epsilon<Scalar>(); + Scalar d = this->eigen2_dot(other); + Scalar absD = ei_abs(d); + + Scalar scale0; + Scalar scale1; + + if (absD>=one) + { + scale0 = Scalar(1) - t; + scale1 = t; + } + else + { + // theta is the angle between the 2 quaternions + Scalar theta = std::acos(absD); + Scalar sinTheta = ei_sin(theta); + + scale0 = ei_sin( ( Scalar(1) - t ) * theta) / sinTheta; + scale1 = ei_sin( ( t * theta) ) / sinTheta; + if (d<0) + scale1 = -scale1; + } + + return Quaternion<Scalar>(scale0 * coeffs() + scale1 * other.coeffs()); +} + +// set from a rotation matrix +template<typename Other> +struct ei_quaternion_assign_impl<Other,3,3> +{ + typedef typename Other::Scalar Scalar; + static inline void run(Quaternion<Scalar>& q, const Other& mat) + { + // This algorithm comes from "Quaternion Calculus and Fast Animation", + // Ken Shoemake, 1987 SIGGRAPH course notes + Scalar t = mat.trace(); + if (t > 0) + { + t = ei_sqrt(t + Scalar(1.0)); + q.w() = Scalar(0.5)*t; + t = Scalar(0.5)/t; + q.x() = (mat.coeff(2,1) - mat.coeff(1,2)) * t; + q.y() = (mat.coeff(0,2) - mat.coeff(2,0)) * t; + q.z() = (mat.coeff(1,0) - mat.coeff(0,1)) * t; + } + else + { + int i = 0; + if (mat.coeff(1,1) > mat.coeff(0,0)) + i = 1; + if (mat.coeff(2,2) > mat.coeff(i,i)) + i = 2; + int j = (i+1)%3; + int k = (j+1)%3; + + t = ei_sqrt(mat.coeff(i,i)-mat.coeff(j,j)-mat.coeff(k,k) + Scalar(1.0)); + q.coeffs().coeffRef(i) = Scalar(0.5) * t; + t = Scalar(0.5)/t; + q.w() = (mat.coeff(k,j)-mat.coeff(j,k))*t; + q.coeffs().coeffRef(j) = (mat.coeff(j,i)+mat.coeff(i,j))*t; + q.coeffs().coeffRef(k) = (mat.coeff(k,i)+mat.coeff(i,k))*t; + } + } +}; + +// set from a vector of coefficients assumed to be a quaternion +template<typename Other> +struct ei_quaternion_assign_impl<Other,4,1> +{ + typedef typename Other::Scalar Scalar; + static inline void run(Quaternion<Scalar>& q, const Other& vec) + { + q.coeffs() = vec; + } +}; + +} // end namespace Eigen diff --git a/eigen/Eigen/src/Eigen2Support/Geometry/Rotation2D.h b/eigen/Eigen/src/Eigen2Support/Geometry/Rotation2D.h new file mode 100644 index 0000000..19b8582 --- /dev/null +++ b/eigen/Eigen/src/Eigen2Support/Geometry/Rotation2D.h @@ -0,0 +1,145 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +// no include guard, we'll include this twice from All.h from Eigen2Support, and it's internal anyway + +namespace Eigen { + +/** \geometry_module \ingroup Geometry_Module + * + * \class Rotation2D + * + * \brief Represents a rotation/orientation in a 2 dimensional space. + * + * \param _Scalar the scalar type, i.e., the type of the coefficients + * + * This class is equivalent to a single scalar representing a counter clock wise rotation + * as a single angle in radian. It provides some additional features such as the automatic + * conversion from/to a 2x2 rotation matrix. Moreover this class aims to provide a similar + * interface to Quaternion in order to facilitate the writing of generic algorithms + * dealing with rotations. + * + * \sa class Quaternion, class Transform + */ +template<typename _Scalar> struct ei_traits<Rotation2D<_Scalar> > +{ + typedef _Scalar Scalar; +}; + +template<typename _Scalar> +class Rotation2D : public RotationBase<Rotation2D<_Scalar>,2> +{ + typedef RotationBase<Rotation2D<_Scalar>,2> Base; + +public: + + using Base::operator*; + + enum { Dim = 2 }; + /** the scalar type of the coefficients */ + typedef _Scalar Scalar; + typedef Matrix<Scalar,2,1> Vector2; + typedef Matrix<Scalar,2,2> Matrix2; + +protected: + + Scalar m_angle; + +public: + + /** Construct a 2D counter clock wise rotation from the angle \a a in radian. */ + inline Rotation2D(Scalar a) : m_angle(a) {} + + /** \returns the rotation angle */ + inline Scalar angle() const { return m_angle; } + + /** \returns a read-write reference to the rotation angle */ + inline Scalar& angle() { return m_angle; } + + /** \returns the inverse rotation */ + inline Rotation2D inverse() const { return -m_angle; } + + /** Concatenates two rotations */ + inline Rotation2D operator*(const Rotation2D& other) const + { return m_angle + other.m_angle; } + + /** Concatenates two rotations */ + inline Rotation2D& operator*=(const Rotation2D& other) + { return m_angle += other.m_angle; return *this; } + + /** Applies the rotation to a 2D vector */ + Vector2 operator* (const Vector2& vec) const + { return toRotationMatrix() * vec; } + + template<typename Derived> + Rotation2D& fromRotationMatrix(const MatrixBase<Derived>& m); + Matrix2 toRotationMatrix(void) const; + + /** \returns the spherical interpolation between \c *this and \a other using + * parameter \a t. It is in fact equivalent to a linear interpolation. + */ + inline Rotation2D slerp(Scalar t, const Rotation2D& other) const + { return m_angle * (1-t) + other.angle() * t; } + + /** \returns \c *this with scalar type casted to \a NewScalarType + * + * Note that if \a NewScalarType is equal to the current scalar type of \c *this + * then this function smartly returns a const reference to \c *this. + */ + template<typename NewScalarType> + inline typename internal::cast_return_type<Rotation2D,Rotation2D<NewScalarType> >::type cast() const + { return typename internal::cast_return_type<Rotation2D,Rotation2D<NewScalarType> >::type(*this); } + + /** Copy constructor with scalar type conversion */ + template<typename OtherScalarType> + inline explicit Rotation2D(const Rotation2D<OtherScalarType>& other) + { + m_angle = Scalar(other.angle()); + } + + /** \returns \c true if \c *this is approximately equal to \a other, within the precision + * determined by \a prec. + * + * \sa MatrixBase::isApprox() */ + bool isApprox(const Rotation2D& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const + { return ei_isApprox(m_angle,other.m_angle, prec); } +}; + +/** \ingroup Geometry_Module + * single precision 2D rotation type */ +typedef Rotation2D<float> Rotation2Df; +/** \ingroup Geometry_Module + * double precision 2D rotation type */ +typedef Rotation2D<double> Rotation2Dd; + +/** Set \c *this from a 2x2 rotation matrix \a mat. + * In other words, this function extract the rotation angle + * from the rotation matrix. + */ +template<typename Scalar> +template<typename Derived> +Rotation2D<Scalar>& Rotation2D<Scalar>::fromRotationMatrix(const MatrixBase<Derived>& mat) +{ + EIGEN_STATIC_ASSERT(Derived::RowsAtCompileTime==2 && Derived::ColsAtCompileTime==2,YOU_MADE_A_PROGRAMMING_MISTAKE) + m_angle = ei_atan2(mat.coeff(1,0), mat.coeff(0,0)); + return *this; +} + +/** Constructs and \returns an equivalent 2x2 rotation matrix. + */ +template<typename Scalar> +typename Rotation2D<Scalar>::Matrix2 +Rotation2D<Scalar>::toRotationMatrix(void) const +{ + Scalar sinA = ei_sin(m_angle); + Scalar cosA = ei_cos(m_angle); + return (Matrix2() << cosA, -sinA, sinA, cosA).finished(); +} + +} // end namespace Eigen diff --git a/eigen/Eigen/src/Eigen2Support/Geometry/RotationBase.h b/eigen/Eigen/src/Eigen2Support/Geometry/RotationBase.h new file mode 100644 index 0000000..b1c8f38 --- /dev/null +++ b/eigen/Eigen/src/Eigen2Support/Geometry/RotationBase.h @@ -0,0 +1,123 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +// no include guard, we'll include this twice from All.h from Eigen2Support, and it's internal anyway + +namespace Eigen { + +// this file aims to contains the various representations of rotation/orientation +// in 2D and 3D space excepted Matrix and Quaternion. + +/** \class RotationBase + * + * \brief Common base class for compact rotation representations + * + * \param Derived is the derived type, i.e., a rotation type + * \param _Dim the dimension of the space + */ +template<typename Derived, int _Dim> +class RotationBase +{ + public: + enum { Dim = _Dim }; + /** the scalar type of the coefficients */ + typedef typename ei_traits<Derived>::Scalar Scalar; + + /** corresponding linear transformation matrix type */ + typedef Matrix<Scalar,Dim,Dim> RotationMatrixType; + + inline const Derived& derived() const { return *static_cast<const Derived*>(this); } + inline Derived& derived() { return *static_cast<Derived*>(this); } + + /** \returns an equivalent rotation matrix */ + inline RotationMatrixType toRotationMatrix() const { return derived().toRotationMatrix(); } + + /** \returns the inverse rotation */ + inline Derived inverse() const { return derived().inverse(); } + + /** \returns the concatenation of the rotation \c *this with a translation \a t */ + inline Transform<Scalar,Dim> operator*(const Translation<Scalar,Dim>& t) const + { return toRotationMatrix() * t; } + + /** \returns the concatenation of the rotation \c *this with a scaling \a s */ + inline RotationMatrixType operator*(const Scaling<Scalar,Dim>& s) const + { return toRotationMatrix() * s; } + + /** \returns the concatenation of the rotation \c *this with an affine transformation \a t */ + inline Transform<Scalar,Dim> operator*(const Transform<Scalar,Dim>& t) const + { return toRotationMatrix() * t; } +}; + +/** \geometry_module + * + * Constructs a Dim x Dim rotation matrix from the rotation \a r + */ +template<typename _Scalar, int _Rows, int _Cols, int _Storage, int _MaxRows, int _MaxCols> +template<typename OtherDerived> +Matrix<_Scalar, _Rows, _Cols, _Storage, _MaxRows, _MaxCols> +::Matrix(const RotationBase<OtherDerived,ColsAtCompileTime>& r) +{ + EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(Matrix,int(OtherDerived::Dim),int(OtherDerived::Dim)) + *this = r.toRotationMatrix(); +} + +/** \geometry_module + * + * Set a Dim x Dim rotation matrix from the rotation \a r + */ +template<typename _Scalar, int _Rows, int _Cols, int _Storage, int _MaxRows, int _MaxCols> +template<typename OtherDerived> +Matrix<_Scalar, _Rows, _Cols, _Storage, _MaxRows, _MaxCols>& +Matrix<_Scalar, _Rows, _Cols, _Storage, _MaxRows, _MaxCols> +::operator=(const RotationBase<OtherDerived,ColsAtCompileTime>& r) +{ + EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(Matrix,int(OtherDerived::Dim),int(OtherDerived::Dim)) + return *this = r.toRotationMatrix(); +} + +/** \internal + * + * Helper function to return an arbitrary rotation object to a rotation matrix. + * + * \param Scalar the numeric type of the matrix coefficients + * \param Dim the dimension of the current space + * + * It returns a Dim x Dim fixed size matrix. + * + * Default specializations are provided for: + * - any scalar type (2D), + * - any matrix expression, + * - any type based on RotationBase (e.g., Quaternion, AngleAxis, Rotation2D) + * + * Currently ei_toRotationMatrix is only used by Transform. + * + * \sa class Transform, class Rotation2D, class Quaternion, class AngleAxis + */ +template<typename Scalar, int Dim> +static inline Matrix<Scalar,2,2> ei_toRotationMatrix(const Scalar& s) +{ + EIGEN_STATIC_ASSERT(Dim==2,YOU_MADE_A_PROGRAMMING_MISTAKE) + return Rotation2D<Scalar>(s).toRotationMatrix(); +} + +template<typename Scalar, int Dim, typename OtherDerived> +static inline Matrix<Scalar,Dim,Dim> ei_toRotationMatrix(const RotationBase<OtherDerived,Dim>& r) +{ + return r.toRotationMatrix(); +} + +template<typename Scalar, int Dim, typename OtherDerived> +static inline const MatrixBase<OtherDerived>& ei_toRotationMatrix(const MatrixBase<OtherDerived>& mat) +{ + EIGEN_STATIC_ASSERT(OtherDerived::RowsAtCompileTime==Dim && OtherDerived::ColsAtCompileTime==Dim, + YOU_MADE_A_PROGRAMMING_MISTAKE) + return mat; +} + +} // end namespace Eigen diff --git a/eigen/Eigen/src/Eigen2Support/Geometry/Scaling.h b/eigen/Eigen/src/Eigen2Support/Geometry/Scaling.h new file mode 100644 index 0000000..b8fa6cd --- /dev/null +++ b/eigen/Eigen/src/Eigen2Support/Geometry/Scaling.h @@ -0,0 +1,167 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +// no include guard, we'll include this twice from All.h from Eigen2Support, and it's internal anyway + +namespace Eigen { + +/** \geometry_module \ingroup Geometry_Module + * + * \class Scaling + * + * \brief Represents a possibly non uniform scaling transformation + * + * \param _Scalar the scalar type, i.e., the type of the coefficients. + * \param _Dim the dimension of the space, can be a compile time value or Dynamic + * + * \note This class is not aimed to be used to store a scaling transformation, + * but rather to make easier the constructions and updates of Transform objects. + * + * \sa class Translation, class Transform + */ +template<typename _Scalar, int _Dim> +class Scaling +{ +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Dim) + /** dimension of the space */ + enum { Dim = _Dim }; + /** the scalar type of the coefficients */ + typedef _Scalar Scalar; + /** corresponding vector type */ + typedef Matrix<Scalar,Dim,1> VectorType; + /** corresponding linear transformation matrix type */ + typedef Matrix<Scalar,Dim,Dim> LinearMatrixType; + /** corresponding translation type */ + typedef Translation<Scalar,Dim> TranslationType; + /** corresponding affine transformation type */ + typedef Transform<Scalar,Dim> TransformType; + +protected: + + VectorType m_coeffs; + +public: + + /** Default constructor without initialization. */ + Scaling() {} + /** Constructs and initialize a uniform scaling transformation */ + explicit inline Scaling(const Scalar& s) { m_coeffs.setConstant(s); } + /** 2D only */ + inline Scaling(const Scalar& sx, const Scalar& sy) + { + ei_assert(Dim==2); + m_coeffs.x() = sx; + m_coeffs.y() = sy; + } + /** 3D only */ + inline Scaling(const Scalar& sx, const Scalar& sy, const Scalar& sz) + { + ei_assert(Dim==3); + m_coeffs.x() = sx; + m_coeffs.y() = sy; + m_coeffs.z() = sz; + } + /** Constructs and initialize the scaling transformation from a vector of scaling coefficients */ + explicit inline Scaling(const VectorType& coeffs) : m_coeffs(coeffs) {} + + const VectorType& coeffs() const { return m_coeffs; } + VectorType& coeffs() { return m_coeffs; } + + /** Concatenates two scaling */ + inline Scaling operator* (const Scaling& other) const + { return Scaling(coeffs().cwise() * other.coeffs()); } + + /** Concatenates a scaling and a translation */ + inline TransformType operator* (const TranslationType& t) const; + + /** Concatenates a scaling and an affine transformation */ + inline TransformType operator* (const TransformType& t) const; + + /** Concatenates a scaling and a linear transformation matrix */ + // TODO returns an expression + inline LinearMatrixType operator* (const LinearMatrixType& other) const + { return coeffs().asDiagonal() * other; } + + /** Concatenates a linear transformation matrix and a scaling */ + // TODO returns an expression + friend inline LinearMatrixType operator* (const LinearMatrixType& other, const Scaling& s) + { return other * s.coeffs().asDiagonal(); } + + template<typename Derived> + inline LinearMatrixType operator*(const RotationBase<Derived,Dim>& r) const + { return *this * r.toRotationMatrix(); } + + /** Applies scaling to vector */ + inline VectorType operator* (const VectorType& other) const + { return coeffs().asDiagonal() * other; } + + /** \returns the inverse scaling */ + inline Scaling inverse() const + { return Scaling(coeffs().cwise().inverse()); } + + inline Scaling& operator=(const Scaling& other) + { + m_coeffs = other.m_coeffs; + return *this; + } + + /** \returns \c *this with scalar type casted to \a NewScalarType + * + * Note that if \a NewScalarType is equal to the current scalar type of \c *this + * then this function smartly returns a const reference to \c *this. + */ + template<typename NewScalarType> + inline typename internal::cast_return_type<Scaling,Scaling<NewScalarType,Dim> >::type cast() const + { return typename internal::cast_return_type<Scaling,Scaling<NewScalarType,Dim> >::type(*this); } + + /** Copy constructor with scalar type conversion */ + template<typename OtherScalarType> + inline explicit Scaling(const Scaling<OtherScalarType,Dim>& other) + { m_coeffs = other.coeffs().template cast<Scalar>(); } + + /** \returns \c true if \c *this is approximately equal to \a other, within the precision + * determined by \a prec. + * + * \sa MatrixBase::isApprox() */ + bool isApprox(const Scaling& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const + { return m_coeffs.isApprox(other.m_coeffs, prec); } + +}; + +/** \addtogroup Geometry_Module */ +//@{ +typedef Scaling<float, 2> Scaling2f; +typedef Scaling<double,2> Scaling2d; +typedef Scaling<float, 3> Scaling3f; +typedef Scaling<double,3> Scaling3d; +//@} + +template<typename Scalar, int Dim> +inline typename Scaling<Scalar,Dim>::TransformType +Scaling<Scalar,Dim>::operator* (const TranslationType& t) const +{ + TransformType res; + res.matrix().setZero(); + res.linear().diagonal() = coeffs(); + res.translation() = m_coeffs.cwise() * t.vector(); + res(Dim,Dim) = Scalar(1); + return res; +} + +template<typename Scalar, int Dim> +inline typename Scaling<Scalar,Dim>::TransformType +Scaling<Scalar,Dim>::operator* (const TransformType& t) const +{ + TransformType res = t; + res.prescale(m_coeffs); + return res; +} + +} // end namespace Eigen diff --git a/eigen/Eigen/src/Eigen2Support/Geometry/Transform.h b/eigen/Eigen/src/Eigen2Support/Geometry/Transform.h new file mode 100644 index 0000000..fab60b2 --- /dev/null +++ b/eigen/Eigen/src/Eigen2Support/Geometry/Transform.h @@ -0,0 +1,786 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +// no include guard, we'll include this twice from All.h from Eigen2Support, and it's internal anyway + +namespace Eigen { + +// Note that we have to pass Dim and HDim because it is not allowed to use a template +// parameter to define a template specialization. To be more precise, in the following +// specializations, it is not allowed to use Dim+1 instead of HDim. +template< typename Other, + int Dim, + int HDim, + int OtherRows=Other::RowsAtCompileTime, + int OtherCols=Other::ColsAtCompileTime> +struct ei_transform_product_impl; + +/** \geometry_module \ingroup Geometry_Module + * + * \class Transform + * + * \brief Represents an homogeneous transformation in a N dimensional space + * + * \param _Scalar the scalar type, i.e., the type of the coefficients + * \param _Dim the dimension of the space + * + * The homography is internally represented and stored as a (Dim+1)^2 matrix which + * is available through the matrix() method. + * + * Conversion methods from/to Qt's QMatrix and QTransform are available if the + * preprocessor token EIGEN_QT_SUPPORT is defined. + * + * \sa class Matrix, class Quaternion + */ +template<typename _Scalar, int _Dim> +class Transform +{ +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Dim==Dynamic ? Dynamic : (_Dim+1)*(_Dim+1)) + enum { + Dim = _Dim, ///< space dimension in which the transformation holds + HDim = _Dim+1 ///< size of a respective homogeneous vector + }; + /** the scalar type of the coefficients */ + typedef _Scalar Scalar; + /** type of the matrix used to represent the transformation */ + typedef Matrix<Scalar,HDim,HDim> MatrixType; + /** type of the matrix used to represent the linear part of the transformation */ + typedef Matrix<Scalar,Dim,Dim> LinearMatrixType; + /** type of read/write reference to the linear part of the transformation */ + typedef Block<MatrixType,Dim,Dim> LinearPart; + /** type of read/write reference to the linear part of the transformation */ + typedef const Block<const MatrixType,Dim,Dim> ConstLinearPart; + /** type of a vector */ + typedef Matrix<Scalar,Dim,1> VectorType; + /** type of a read/write reference to the translation part of the rotation */ + typedef Block<MatrixType,Dim,1> TranslationPart; + /** type of a read/write reference to the translation part of the rotation */ + typedef const Block<const MatrixType,Dim,1> ConstTranslationPart; + /** corresponding translation type */ + typedef Translation<Scalar,Dim> TranslationType; + /** corresponding scaling transformation type */ + typedef Scaling<Scalar,Dim> ScalingType; + +protected: + + MatrixType m_matrix; + +public: + + /** Default constructor without initialization of the coefficients. */ + inline Transform() { } + + inline Transform(const Transform& other) + { + m_matrix = other.m_matrix; + } + + inline explicit Transform(const TranslationType& t) { *this = t; } + inline explicit Transform(const ScalingType& s) { *this = s; } + template<typename Derived> + inline explicit Transform(const RotationBase<Derived, Dim>& r) { *this = r; } + + inline Transform& operator=(const Transform& other) + { m_matrix = other.m_matrix; return *this; } + + template<typename OtherDerived, bool BigMatrix> // MSVC 2005 will commit suicide if BigMatrix has a default value + struct construct_from_matrix + { + static inline void run(Transform *transform, const MatrixBase<OtherDerived>& other) + { + transform->matrix() = other; + } + }; + + template<typename OtherDerived> struct construct_from_matrix<OtherDerived, true> + { + static inline void run(Transform *transform, const MatrixBase<OtherDerived>& other) + { + transform->linear() = other; + transform->translation().setZero(); + transform->matrix()(Dim,Dim) = Scalar(1); + transform->matrix().template block<1,Dim>(Dim,0).setZero(); + } + }; + + /** Constructs and initializes a transformation from a Dim^2 or a (Dim+1)^2 matrix. */ + template<typename OtherDerived> + inline explicit Transform(const MatrixBase<OtherDerived>& other) + { + construct_from_matrix<OtherDerived, int(OtherDerived::RowsAtCompileTime) == Dim>::run(this, other); + } + + /** Set \c *this from a (Dim+1)^2 matrix. */ + template<typename OtherDerived> + inline Transform& operator=(const MatrixBase<OtherDerived>& other) + { m_matrix = other; return *this; } + + #ifdef EIGEN_QT_SUPPORT + inline Transform(const QMatrix& other); + inline Transform& operator=(const QMatrix& other); + inline QMatrix toQMatrix(void) const; + inline Transform(const QTransform& other); + inline Transform& operator=(const QTransform& other); + inline QTransform toQTransform(void) const; + #endif + + /** shortcut for m_matrix(row,col); + * \sa MatrixBase::operaror(int,int) const */ + inline Scalar operator() (int row, int col) const { return m_matrix(row,col); } + /** shortcut for m_matrix(row,col); + * \sa MatrixBase::operaror(int,int) */ + inline Scalar& operator() (int row, int col) { return m_matrix(row,col); } + + /** \returns a read-only expression of the transformation matrix */ + inline const MatrixType& matrix() const { return m_matrix; } + /** \returns a writable expression of the transformation matrix */ + inline MatrixType& matrix() { return m_matrix; } + + /** \returns a read-only expression of the linear (linear) part of the transformation */ + inline ConstLinearPart linear() const { return m_matrix.template block<Dim,Dim>(0,0); } + /** \returns a writable expression of the linear (linear) part of the transformation */ + inline LinearPart linear() { return m_matrix.template block<Dim,Dim>(0,0); } + + /** \returns a read-only expression of the translation vector of the transformation */ + inline ConstTranslationPart translation() const { return m_matrix.template block<Dim,1>(0,Dim); } + /** \returns a writable expression of the translation vector of the transformation */ + inline TranslationPart translation() { return m_matrix.template block<Dim,1>(0,Dim); } + + /** \returns an expression of the product between the transform \c *this and a matrix expression \a other + * + * The right hand side \a other might be either: + * \li a vector of size Dim, + * \li an homogeneous vector of size Dim+1, + * \li a transformation matrix of size Dim+1 x Dim+1. + */ + // note: this function is defined here because some compilers cannot find the respective declaration + template<typename OtherDerived> + inline const typename ei_transform_product_impl<OtherDerived,_Dim,_Dim+1>::ResultType + operator * (const MatrixBase<OtherDerived> &other) const + { return ei_transform_product_impl<OtherDerived,Dim,HDim>::run(*this,other.derived()); } + + /** \returns the product expression of a transformation matrix \a a times a transform \a b + * The transformation matrix \a a must have a Dim+1 x Dim+1 sizes. */ + template<typename OtherDerived> + friend inline const typename ProductReturnType<OtherDerived,MatrixType>::Type + operator * (const MatrixBase<OtherDerived> &a, const Transform &b) + { return a.derived() * b.matrix(); } + + /** Contatenates two transformations */ + inline const Transform + operator * (const Transform& other) const + { return Transform(m_matrix * other.matrix()); } + + /** \sa MatrixBase::setIdentity() */ + void setIdentity() { m_matrix.setIdentity(); } + static const typename MatrixType::IdentityReturnType Identity() + { + return MatrixType::Identity(); + } + + template<typename OtherDerived> + inline Transform& scale(const MatrixBase<OtherDerived> &other); + + template<typename OtherDerived> + inline Transform& prescale(const MatrixBase<OtherDerived> &other); + + inline Transform& scale(Scalar s); + inline Transform& prescale(Scalar s); + + template<typename OtherDerived> + inline Transform& translate(const MatrixBase<OtherDerived> &other); + + template<typename OtherDerived> + inline Transform& pretranslate(const MatrixBase<OtherDerived> &other); + + template<typename RotationType> + inline Transform& rotate(const RotationType& rotation); + + template<typename RotationType> + inline Transform& prerotate(const RotationType& rotation); + + Transform& shear(Scalar sx, Scalar sy); + Transform& preshear(Scalar sx, Scalar sy); + + inline Transform& operator=(const TranslationType& t); + inline Transform& operator*=(const TranslationType& t) { return translate(t.vector()); } + inline Transform operator*(const TranslationType& t) const; + + inline Transform& operator=(const ScalingType& t); + inline Transform& operator*=(const ScalingType& s) { return scale(s.coeffs()); } + inline Transform operator*(const ScalingType& s) const; + friend inline Transform operator*(const LinearMatrixType& mat, const Transform& t) + { + Transform res = t; + res.matrix().row(Dim) = t.matrix().row(Dim); + res.matrix().template block<Dim,HDim>(0,0) = (mat * t.matrix().template block<Dim,HDim>(0,0)).lazy(); + return res; + } + + template<typename Derived> + inline Transform& operator=(const RotationBase<Derived,Dim>& r); + template<typename Derived> + inline Transform& operator*=(const RotationBase<Derived,Dim>& r) { return rotate(r.toRotationMatrix()); } + template<typename Derived> + inline Transform operator*(const RotationBase<Derived,Dim>& r) const; + + LinearMatrixType rotation() const; + template<typename RotationMatrixType, typename ScalingMatrixType> + void computeRotationScaling(RotationMatrixType *rotation, ScalingMatrixType *scaling) const; + template<typename ScalingMatrixType, typename RotationMatrixType> + void computeScalingRotation(ScalingMatrixType *scaling, RotationMatrixType *rotation) const; + + template<typename PositionDerived, typename OrientationType, typename ScaleDerived> + Transform& fromPositionOrientationScale(const MatrixBase<PositionDerived> &position, + const OrientationType& orientation, const MatrixBase<ScaleDerived> &scale); + + inline const MatrixType inverse(TransformTraits traits = Affine) const; + + /** \returns a const pointer to the column major internal matrix */ + const Scalar* data() const { return m_matrix.data(); } + /** \returns a non-const pointer to the column major internal matrix */ + Scalar* data() { return m_matrix.data(); } + + /** \returns \c *this with scalar type casted to \a NewScalarType + * + * Note that if \a NewScalarType is equal to the current scalar type of \c *this + * then this function smartly returns a const reference to \c *this. + */ + template<typename NewScalarType> + inline typename internal::cast_return_type<Transform,Transform<NewScalarType,Dim> >::type cast() const + { return typename internal::cast_return_type<Transform,Transform<NewScalarType,Dim> >::type(*this); } + + /** Copy constructor with scalar type conversion */ + template<typename OtherScalarType> + inline explicit Transform(const Transform<OtherScalarType,Dim>& other) + { m_matrix = other.matrix().template cast<Scalar>(); } + + /** \returns \c true if \c *this is approximately equal to \a other, within the precision + * determined by \a prec. + * + * \sa MatrixBase::isApprox() */ + bool isApprox(const Transform& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const + { return m_matrix.isApprox(other.m_matrix, prec); } + + #ifdef EIGEN_TRANSFORM_PLUGIN + #include EIGEN_TRANSFORM_PLUGIN + #endif + +protected: + +}; + +/** \ingroup Geometry_Module */ +typedef Transform<float,2> Transform2f; +/** \ingroup Geometry_Module */ +typedef Transform<float,3> Transform3f; +/** \ingroup Geometry_Module */ +typedef Transform<double,2> Transform2d; +/** \ingroup Geometry_Module */ +typedef Transform<double,3> Transform3d; + +/************************** +*** Optional QT support *** +**************************/ + +#ifdef EIGEN_QT_SUPPORT +/** Initialises \c *this from a QMatrix assuming the dimension is 2. + * + * This function is available only if the token EIGEN_QT_SUPPORT is defined. + */ +template<typename Scalar, int Dim> +Transform<Scalar,Dim>::Transform(const QMatrix& other) +{ + *this = other; +} + +/** Set \c *this from a QMatrix assuming the dimension is 2. + * + * This function is available only if the token EIGEN_QT_SUPPORT is defined. + */ +template<typename Scalar, int Dim> +Transform<Scalar,Dim>& Transform<Scalar,Dim>::operator=(const QMatrix& other) +{ + EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE) + m_matrix << other.m11(), other.m21(), other.dx(), + other.m12(), other.m22(), other.dy(), + 0, 0, 1; + return *this; +} + +/** \returns a QMatrix from \c *this assuming the dimension is 2. + * + * \warning this convertion might loss data if \c *this is not affine + * + * This function is available only if the token EIGEN_QT_SUPPORT is defined. + */ +template<typename Scalar, int Dim> +QMatrix Transform<Scalar,Dim>::toQMatrix(void) const +{ + EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE) + return QMatrix(m_matrix.coeff(0,0), m_matrix.coeff(1,0), + m_matrix.coeff(0,1), m_matrix.coeff(1,1), + m_matrix.coeff(0,2), m_matrix.coeff(1,2)); +} + +/** Initialises \c *this from a QTransform assuming the dimension is 2. + * + * This function is available only if the token EIGEN_QT_SUPPORT is defined. + */ +template<typename Scalar, int Dim> +Transform<Scalar,Dim>::Transform(const QTransform& other) +{ + *this = other; +} + +/** Set \c *this from a QTransform assuming the dimension is 2. + * + * This function is available only if the token EIGEN_QT_SUPPORT is defined. + */ +template<typename Scalar, int Dim> +Transform<Scalar,Dim>& Transform<Scalar,Dim>::operator=(const QTransform& other) +{ + EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE) + m_matrix << other.m11(), other.m21(), other.dx(), + other.m12(), other.m22(), other.dy(), + other.m13(), other.m23(), other.m33(); + return *this; +} + +/** \returns a QTransform from \c *this assuming the dimension is 2. + * + * This function is available only if the token EIGEN_QT_SUPPORT is defined. + */ +template<typename Scalar, int Dim> +QTransform Transform<Scalar,Dim>::toQTransform(void) const +{ + EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE) + return QTransform(m_matrix.coeff(0,0), m_matrix.coeff(1,0), m_matrix.coeff(2,0), + m_matrix.coeff(0,1), m_matrix.coeff(1,1), m_matrix.coeff(2,1), + m_matrix.coeff(0,2), m_matrix.coeff(1,2), m_matrix.coeff(2,2)); +} +#endif + +/********************* +*** Procedural API *** +*********************/ + +/** Applies on the right the non uniform scale transformation represented + * by the vector \a other to \c *this and returns a reference to \c *this. + * \sa prescale() + */ +template<typename Scalar, int Dim> +template<typename OtherDerived> +Transform<Scalar,Dim>& +Transform<Scalar,Dim>::scale(const MatrixBase<OtherDerived> &other) +{ + EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived,int(Dim)) + linear() = (linear() * other.asDiagonal()).lazy(); + return *this; +} + +/** Applies on the right a uniform scale of a factor \a c to \c *this + * and returns a reference to \c *this. + * \sa prescale(Scalar) + */ +template<typename Scalar, int Dim> +inline Transform<Scalar,Dim>& Transform<Scalar,Dim>::scale(Scalar s) +{ + linear() *= s; + return *this; +} + +/** Applies on the left the non uniform scale transformation represented + * by the vector \a other to \c *this and returns a reference to \c *this. + * \sa scale() + */ +template<typename Scalar, int Dim> +template<typename OtherDerived> +Transform<Scalar,Dim>& +Transform<Scalar,Dim>::prescale(const MatrixBase<OtherDerived> &other) +{ + EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived,int(Dim)) + m_matrix.template block<Dim,HDim>(0,0) = (other.asDiagonal() * m_matrix.template block<Dim,HDim>(0,0)).lazy(); + return *this; +} + +/** Applies on the left a uniform scale of a factor \a c to \c *this + * and returns a reference to \c *this. + * \sa scale(Scalar) + */ +template<typename Scalar, int Dim> +inline Transform<Scalar,Dim>& Transform<Scalar,Dim>::prescale(Scalar s) +{ + m_matrix.template corner<Dim,HDim>(TopLeft) *= s; + return *this; +} + +/** Applies on the right the translation matrix represented by the vector \a other + * to \c *this and returns a reference to \c *this. + * \sa pretranslate() + */ +template<typename Scalar, int Dim> +template<typename OtherDerived> +Transform<Scalar,Dim>& +Transform<Scalar,Dim>::translate(const MatrixBase<OtherDerived> &other) +{ + EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived,int(Dim)) + translation() += linear() * other; + return *this; +} + +/** Applies on the left the translation matrix represented by the vector \a other + * to \c *this and returns a reference to \c *this. + * \sa translate() + */ +template<typename Scalar, int Dim> +template<typename OtherDerived> +Transform<Scalar,Dim>& +Transform<Scalar,Dim>::pretranslate(const MatrixBase<OtherDerived> &other) +{ + EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived,int(Dim)) + translation() += other; + return *this; +} + +/** Applies on the right the rotation represented by the rotation \a rotation + * to \c *this and returns a reference to \c *this. + * + * The template parameter \a RotationType is the type of the rotation which + * must be known by ei_toRotationMatrix<>. + * + * Natively supported types includes: + * - any scalar (2D), + * - a Dim x Dim matrix expression, + * - a Quaternion (3D), + * - a AngleAxis (3D) + * + * This mechanism is easily extendable to support user types such as Euler angles, + * or a pair of Quaternion for 4D rotations. + * + * \sa rotate(Scalar), class Quaternion, class AngleAxis, prerotate(RotationType) + */ +template<typename Scalar, int Dim> +template<typename RotationType> +Transform<Scalar,Dim>& +Transform<Scalar,Dim>::rotate(const RotationType& rotation) +{ + linear() *= ei_toRotationMatrix<Scalar,Dim>(rotation); + return *this; +} + +/** Applies on the left the rotation represented by the rotation \a rotation + * to \c *this and returns a reference to \c *this. + * + * See rotate() for further details. + * + * \sa rotate() + */ +template<typename Scalar, int Dim> +template<typename RotationType> +Transform<Scalar,Dim>& +Transform<Scalar,Dim>::prerotate(const RotationType& rotation) +{ + m_matrix.template block<Dim,HDim>(0,0) = ei_toRotationMatrix<Scalar,Dim>(rotation) + * m_matrix.template block<Dim,HDim>(0,0); + return *this; +} + +/** Applies on the right the shear transformation represented + * by the vector \a other to \c *this and returns a reference to \c *this. + * \warning 2D only. + * \sa preshear() + */ +template<typename Scalar, int Dim> +Transform<Scalar,Dim>& +Transform<Scalar,Dim>::shear(Scalar sx, Scalar sy) +{ + EIGEN_STATIC_ASSERT(int(Dim)==2, YOU_MADE_A_PROGRAMMING_MISTAKE) + VectorType tmp = linear().col(0)*sy + linear().col(1); + linear() << linear().col(0) + linear().col(1)*sx, tmp; + return *this; +} + +/** Applies on the left the shear transformation represented + * by the vector \a other to \c *this and returns a reference to \c *this. + * \warning 2D only. + * \sa shear() + */ +template<typename Scalar, int Dim> +Transform<Scalar,Dim>& +Transform<Scalar,Dim>::preshear(Scalar sx, Scalar sy) +{ + EIGEN_STATIC_ASSERT(int(Dim)==2, YOU_MADE_A_PROGRAMMING_MISTAKE) + m_matrix.template block<Dim,HDim>(0,0) = LinearMatrixType(1, sx, sy, 1) * m_matrix.template block<Dim,HDim>(0,0); + return *this; +} + +/****************************************************** +*** Scaling, Translation and Rotation compatibility *** +******************************************************/ + +template<typename Scalar, int Dim> +inline Transform<Scalar,Dim>& Transform<Scalar,Dim>::operator=(const TranslationType& t) +{ + linear().setIdentity(); + translation() = t.vector(); + m_matrix.template block<1,Dim>(Dim,0).setZero(); + m_matrix(Dim,Dim) = Scalar(1); + return *this; +} + +template<typename Scalar, int Dim> +inline Transform<Scalar,Dim> Transform<Scalar,Dim>::operator*(const TranslationType& t) const +{ + Transform res = *this; + res.translate(t.vector()); + return res; +} + +template<typename Scalar, int Dim> +inline Transform<Scalar,Dim>& Transform<Scalar,Dim>::operator=(const ScalingType& s) +{ + m_matrix.setZero(); + linear().diagonal() = s.coeffs(); + m_matrix.coeffRef(Dim,Dim) = Scalar(1); + return *this; +} + +template<typename Scalar, int Dim> +inline Transform<Scalar,Dim> Transform<Scalar,Dim>::operator*(const ScalingType& s) const +{ + Transform res = *this; + res.scale(s.coeffs()); + return res; +} + +template<typename Scalar, int Dim> +template<typename Derived> +inline Transform<Scalar,Dim>& Transform<Scalar,Dim>::operator=(const RotationBase<Derived,Dim>& r) +{ + linear() = ei_toRotationMatrix<Scalar,Dim>(r); + translation().setZero(); + m_matrix.template block<1,Dim>(Dim,0).setZero(); + m_matrix.coeffRef(Dim,Dim) = Scalar(1); + return *this; +} + +template<typename Scalar, int Dim> +template<typename Derived> +inline Transform<Scalar,Dim> Transform<Scalar,Dim>::operator*(const RotationBase<Derived,Dim>& r) const +{ + Transform res = *this; + res.rotate(r.derived()); + return res; +} + +/************************ +*** Special functions *** +************************/ + +/** \returns the rotation part of the transformation + * \nonstableyet + * + * \svd_module + * + * \sa computeRotationScaling(), computeScalingRotation(), class SVD + */ +template<typename Scalar, int Dim> +typename Transform<Scalar,Dim>::LinearMatrixType +Transform<Scalar,Dim>::rotation() const +{ + LinearMatrixType result; + computeRotationScaling(&result, (LinearMatrixType*)0); + return result; +} + + +/** decomposes the linear part of the transformation as a product rotation x scaling, the scaling being + * not necessarily positive. + * + * If either pointer is zero, the corresponding computation is skipped. + * + * \nonstableyet + * + * \svd_module + * + * \sa computeScalingRotation(), rotation(), class SVD + */ +template<typename Scalar, int Dim> +template<typename RotationMatrixType, typename ScalingMatrixType> +void Transform<Scalar,Dim>::computeRotationScaling(RotationMatrixType *rotation, ScalingMatrixType *scaling) const +{ + JacobiSVD<LinearMatrixType> svd(linear(), ComputeFullU|ComputeFullV); + Scalar x = (svd.matrixU() * svd.matrixV().adjoint()).determinant(); // so x has absolute value 1 + Matrix<Scalar, Dim, 1> sv(svd.singularValues()); + sv.coeffRef(0) *= x; + if(scaling) + { + scaling->noalias() = svd.matrixV() * sv.asDiagonal() * svd.matrixV().adjoint(); + } + if(rotation) + { + LinearMatrixType m(svd.matrixU()); + m.col(0) /= x; + rotation->noalias() = m * svd.matrixV().adjoint(); + } +} + +/** decomposes the linear part of the transformation as a product rotation x scaling, the scaling being + * not necessarily positive. + * + * If either pointer is zero, the corresponding computation is skipped. + * + * \nonstableyet + * + * \svd_module + * + * \sa computeRotationScaling(), rotation(), class SVD + */ +template<typename Scalar, int Dim> +template<typename ScalingMatrixType, typename RotationMatrixType> +void Transform<Scalar,Dim>::computeScalingRotation(ScalingMatrixType *scaling, RotationMatrixType *rotation) const +{ + JacobiSVD<LinearMatrixType> svd(linear(), ComputeFullU|ComputeFullV); + Scalar x = (svd.matrixU() * svd.matrixV().adjoint()).determinant(); // so x has absolute value 1 + Matrix<Scalar, Dim, 1> sv(svd.singularValues()); + sv.coeffRef(0) *= x; + if(scaling) + { + scaling->noalias() = svd.matrixU() * sv.asDiagonal() * svd.matrixU().adjoint(); + } + if(rotation) + { + LinearMatrixType m(svd.matrixU()); + m.col(0) /= x; + rotation->noalias() = m * svd.matrixV().adjoint(); + } +} + +/** Convenient method to set \c *this from a position, orientation and scale + * of a 3D object. + */ +template<typename Scalar, int Dim> +template<typename PositionDerived, typename OrientationType, typename ScaleDerived> +Transform<Scalar,Dim>& +Transform<Scalar,Dim>::fromPositionOrientationScale(const MatrixBase<PositionDerived> &position, + const OrientationType& orientation, const MatrixBase<ScaleDerived> &scale) +{ + linear() = ei_toRotationMatrix<Scalar,Dim>(orientation); + linear() *= scale.asDiagonal(); + translation() = position; + m_matrix.template block<1,Dim>(Dim,0).setZero(); + m_matrix(Dim,Dim) = Scalar(1); + return *this; +} + +/** \nonstableyet + * + * \returns the inverse transformation matrix according to some given knowledge + * on \c *this. + * + * \param traits allows to optimize the inversion process when the transformion + * is known to be not a general transformation. The possible values are: + * - Projective if the transformation is not necessarily affine, i.e., if the + * last row is not guaranteed to be [0 ... 0 1] + * - Affine is the default, the last row is assumed to be [0 ... 0 1] + * - Isometry if the transformation is only a concatenations of translations + * and rotations. + * + * \warning unless \a traits is always set to NoShear or NoScaling, this function + * requires the generic inverse method of MatrixBase defined in the LU module. If + * you forget to include this module, then you will get hard to debug linking errors. + * + * \sa MatrixBase::inverse() + */ +template<typename Scalar, int Dim> +inline const typename Transform<Scalar,Dim>::MatrixType +Transform<Scalar,Dim>::inverse(TransformTraits traits) const +{ + if (traits == Projective) + { + return m_matrix.inverse(); + } + else + { + MatrixType res; + if (traits == Affine) + { + res.template corner<Dim,Dim>(TopLeft) = linear().inverse(); + } + else if (traits == Isometry) + { + res.template corner<Dim,Dim>(TopLeft) = linear().transpose(); + } + else + { + ei_assert("invalid traits value in Transform::inverse()"); + } + // translation and remaining parts + res.template corner<Dim,1>(TopRight) = - res.template corner<Dim,Dim>(TopLeft) * translation(); + res.template corner<1,Dim>(BottomLeft).setZero(); + res.coeffRef(Dim,Dim) = Scalar(1); + return res; + } +} + +/***************************************************** +*** Specializations of operator* with a MatrixBase *** +*****************************************************/ + +template<typename Other, int Dim, int HDim> +struct ei_transform_product_impl<Other,Dim,HDim, HDim,HDim> +{ + typedef Transform<typename Other::Scalar,Dim> TransformType; + typedef typename TransformType::MatrixType MatrixType; + typedef typename ProductReturnType<MatrixType,Other>::Type ResultType; + static ResultType run(const TransformType& tr, const Other& other) + { return tr.matrix() * other; } +}; + +template<typename Other, int Dim, int HDim> +struct ei_transform_product_impl<Other,Dim,HDim, Dim,Dim> +{ + typedef Transform<typename Other::Scalar,Dim> TransformType; + typedef typename TransformType::MatrixType MatrixType; + typedef TransformType ResultType; + static ResultType run(const TransformType& tr, const Other& other) + { + TransformType res; + res.translation() = tr.translation(); + res.matrix().row(Dim) = tr.matrix().row(Dim); + res.linear() = (tr.linear() * other).lazy(); + return res; + } +}; + +template<typename Other, int Dim, int HDim> +struct ei_transform_product_impl<Other,Dim,HDim, HDim,1> +{ + typedef Transform<typename Other::Scalar,Dim> TransformType; + typedef typename TransformType::MatrixType MatrixType; + typedef typename ProductReturnType<MatrixType,Other>::Type ResultType; + static ResultType run(const TransformType& tr, const Other& other) + { return tr.matrix() * other; } +}; + +template<typename Other, int Dim, int HDim> +struct ei_transform_product_impl<Other,Dim,HDim, Dim,1> +{ + typedef typename Other::Scalar Scalar; + typedef Transform<Scalar,Dim> TransformType; + typedef Matrix<Scalar,Dim,1> ResultType; + static ResultType run(const TransformType& tr, const Other& other) + { return ((tr.linear() * other) + tr.translation()) + * (Scalar(1) / ( (tr.matrix().template block<1,Dim>(Dim,0) * other).coeff(0) + tr.matrix().coeff(Dim,Dim))); } +}; + +} // end namespace Eigen diff --git a/eigen/Eigen/src/Eigen2Support/Geometry/Translation.h b/eigen/Eigen/src/Eigen2Support/Geometry/Translation.h new file mode 100644 index 0000000..2b9859f --- /dev/null +++ b/eigen/Eigen/src/Eigen2Support/Geometry/Translation.h @@ -0,0 +1,184 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +// no include guard, we'll include this twice from All.h from Eigen2Support, and it's internal anyway + +namespace Eigen { + +/** \geometry_module \ingroup Geometry_Module + * + * \class Translation + * + * \brief Represents a translation transformation + * + * \param _Scalar the scalar type, i.e., the type of the coefficients. + * \param _Dim the dimension of the space, can be a compile time value or Dynamic + * + * \note This class is not aimed to be used to store a translation transformation, + * but rather to make easier the constructions and updates of Transform objects. + * + * \sa class Scaling, class Transform + */ +template<typename _Scalar, int _Dim> +class Translation +{ +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Dim) + /** dimension of the space */ + enum { Dim = _Dim }; + /** the scalar type of the coefficients */ + typedef _Scalar Scalar; + /** corresponding vector type */ + typedef Matrix<Scalar,Dim,1> VectorType; + /** corresponding linear transformation matrix type */ + typedef Matrix<Scalar,Dim,Dim> LinearMatrixType; + /** corresponding scaling transformation type */ + typedef Scaling<Scalar,Dim> ScalingType; + /** corresponding affine transformation type */ + typedef Transform<Scalar,Dim> TransformType; + +protected: + + VectorType m_coeffs; + +public: + + /** Default constructor without initialization. */ + Translation() {} + /** */ + inline Translation(const Scalar& sx, const Scalar& sy) + { + ei_assert(Dim==2); + m_coeffs.x() = sx; + m_coeffs.y() = sy; + } + /** */ + inline Translation(const Scalar& sx, const Scalar& sy, const Scalar& sz) + { + ei_assert(Dim==3); + m_coeffs.x() = sx; + m_coeffs.y() = sy; + m_coeffs.z() = sz; + } + /** Constructs and initialize the scaling transformation from a vector of scaling coefficients */ + explicit inline Translation(const VectorType& vector) : m_coeffs(vector) {} + + const VectorType& vector() const { return m_coeffs; } + VectorType& vector() { return m_coeffs; } + + /** Concatenates two translation */ + inline Translation operator* (const Translation& other) const + { return Translation(m_coeffs + other.m_coeffs); } + + /** Concatenates a translation and a scaling */ + inline TransformType operator* (const ScalingType& other) const; + + /** Concatenates a translation and a linear transformation */ + inline TransformType operator* (const LinearMatrixType& linear) const; + + template<typename Derived> + inline TransformType operator*(const RotationBase<Derived,Dim>& r) const + { return *this * r.toRotationMatrix(); } + + /** Concatenates a linear transformation and a translation */ + // its a nightmare to define a templated friend function outside its declaration + friend inline TransformType operator* (const LinearMatrixType& linear, const Translation& t) + { + TransformType res; + res.matrix().setZero(); + res.linear() = linear; + res.translation() = linear * t.m_coeffs; + res.matrix().row(Dim).setZero(); + res(Dim,Dim) = Scalar(1); + return res; + } + + /** Concatenates a translation and an affine transformation */ + inline TransformType operator* (const TransformType& t) const; + + /** Applies translation to vector */ + inline VectorType operator* (const VectorType& other) const + { return m_coeffs + other; } + + /** \returns the inverse translation (opposite) */ + Translation inverse() const { return Translation(-m_coeffs); } + + Translation& operator=(const Translation& other) + { + m_coeffs = other.m_coeffs; + return *this; + } + + /** \returns \c *this with scalar type casted to \a NewScalarType + * + * Note that if \a NewScalarType is equal to the current scalar type of \c *this + * then this function smartly returns a const reference to \c *this. + */ + template<typename NewScalarType> + inline typename internal::cast_return_type<Translation,Translation<NewScalarType,Dim> >::type cast() const + { return typename internal::cast_return_type<Translation,Translation<NewScalarType,Dim> >::type(*this); } + + /** Copy constructor with scalar type conversion */ + template<typename OtherScalarType> + inline explicit Translation(const Translation<OtherScalarType,Dim>& other) + { m_coeffs = other.vector().template cast<Scalar>(); } + + /** \returns \c true if \c *this is approximately equal to \a other, within the precision + * determined by \a prec. + * + * \sa MatrixBase::isApprox() */ + bool isApprox(const Translation& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const + { return m_coeffs.isApprox(other.m_coeffs, prec); } + +}; + +/** \addtogroup Geometry_Module */ +//@{ +typedef Translation<float, 2> Translation2f; +typedef Translation<double,2> Translation2d; +typedef Translation<float, 3> Translation3f; +typedef Translation<double,3> Translation3d; +//@} + + +template<typename Scalar, int Dim> +inline typename Translation<Scalar,Dim>::TransformType +Translation<Scalar,Dim>::operator* (const ScalingType& other) const +{ + TransformType res; + res.matrix().setZero(); + res.linear().diagonal() = other.coeffs(); + res.translation() = m_coeffs; + res(Dim,Dim) = Scalar(1); + return res; +} + +template<typename Scalar, int Dim> +inline typename Translation<Scalar,Dim>::TransformType +Translation<Scalar,Dim>::operator* (const LinearMatrixType& linear) const +{ + TransformType res; + res.matrix().setZero(); + res.linear() = linear; + res.translation() = m_coeffs; + res.matrix().row(Dim).setZero(); + res(Dim,Dim) = Scalar(1); + return res; +} + +template<typename Scalar, int Dim> +inline typename Translation<Scalar,Dim>::TransformType +Translation<Scalar,Dim>::operator* (const TransformType& t) const +{ + TransformType res = t; + res.pretranslate(m_coeffs); + return res; +} + +} // end namespace Eigen |