diff options
author | Gael Guennebaud <g.gael@free.fr> | 2008-06-15 11:54:18 +0000 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2008-06-15 11:54:18 +0000 |
commit | 0ee6b08128de893150d7d2d07081bf33b4cb9fb1 (patch) | |
tree | 1d22a187485c99fb9f280d1824dab59c9c1ea06f | |
parent | fbbd8afe3086cc1a1f455cc5264f0b6fc063f8b6 (diff) |
* split Product to a DiagonalProduct template specialization
to optimize matrix-diag and diag-matrix products without
making Product over complicated.
* compilation fixes in Tridiagonalization and HessenbergDecomposition
in the case of 2x2 matrices.
* added an Orientation2D small class with similar interface than Quaternion
(used by Transform to handle 2D and 3D orientations seamlessly)
* added a couple of features in Transform.
-rw-r--r-- | Eigen/Core | 1 | ||||
-rw-r--r-- | Eigen/src/Core/DiagonalProduct.h | 108 | ||||
-rw-r--r-- | Eigen/src/Core/Product.h | 58 | ||||
-rw-r--r-- | Eigen/src/Core/Transpose.h | 2 | ||||
-rw-r--r-- | Eigen/src/Core/util/Constants.h | 4 | ||||
-rw-r--r-- | Eigen/src/Geometry/Transform.h | 239 | ||||
-rwxr-xr-x | Eigen/src/QR/HessenbergDecomposition.h | 3 | ||||
-rwxr-xr-x | Eigen/src/QR/Tridiagonalization.h | 7 | ||||
-rw-r--r-- | test/geometry.cpp | 3 |
9 files changed, 374 insertions, 51 deletions
diff --git a/Eigen/Core b/Eigen/Core index 7129972a9..e08bde4ab 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -41,6 +41,7 @@ namespace Eigen { #include "src/Core/CwiseNullaryOp.h" #include "src/Core/InverseProduct.h" #include "src/Core/Product.h" +#include "src/Core/DiagonalProduct.h" #include "src/Core/Block.h" #include "src/Core/Minor.h" #include "src/Core/Transpose.h" diff --git a/Eigen/src/Core/DiagonalProduct.h b/Eigen/src/Core/DiagonalProduct.h new file mode 100644 index 000000000..5663942ff --- /dev/null +++ b/Eigen/src/Core/DiagonalProduct.h @@ -0,0 +1,108 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2006-2008 Benoit Jacob <jacob@math.jussieu.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_DIAGONALPRODUCT_H +#define EIGEN_DIAGONALPRODUCT_H + +template<typename Lhs, typename Rhs> +struct ei_traits<Product<Lhs, Rhs, DiagonalProduct> > +{ + typedef typename Lhs::Scalar Scalar; + typedef typename ei_nested<Lhs>::type LhsNested; + typedef typename ei_nested<Rhs>::type RhsNested; + typedef typename ei_unref<LhsNested>::type _LhsNested; + typedef typename ei_unref<RhsNested>::type _RhsNested; + enum { + LhsFlags = _LhsNested::Flags, + RhsFlags = _RhsNested::Flags, + RowsAtCompileTime = Lhs::RowsAtCompileTime, + ColsAtCompileTime = Rhs::ColsAtCompileTime, + MaxRowsAtCompileTime = Lhs::MaxRowsAtCompileTime, + MaxColsAtCompileTime = Rhs::MaxColsAtCompileTime, + _RhsVectorizable = (RhsFlags & RowMajorBit) && (RhsFlags & VectorizableBit) + && (ColsAtCompileTime % ei_packet_traits<Scalar>::size == 0), + _LhsVectorizable = (!(LhsFlags & RowMajorBit)) && (LhsFlags & VectorizableBit) + && (RowsAtCompileTime % ei_packet_traits<Scalar>::size == 0), + _LostBits = ~(((RhsFlags & RowMajorBit) && (!_LhsVectorizable) ? 0 : RowMajorBit) + | ((RowsAtCompileTime == Dynamic || ColsAtCompileTime == Dynamic) ? 0 : LargeBit)), + Flags = ((unsigned int)(LhsFlags | RhsFlags) & HereditaryBits & _LostBits) + | (_LhsVectorizable || _RhsVectorizable ? VectorizableBit : 0), + CoeffReadCost = NumTraits<Scalar>::MulCost + _LhsNested::CoeffReadCost + _RhsNested::CoeffReadCost + }; +}; + +template<typename Lhs, typename Rhs> class Product<Lhs, Rhs, DiagonalProduct> : ei_no_assignment_operator, + public MatrixBase<Product<Lhs, Rhs, DiagonalProduct> > +{ + public: + + EIGEN_GENERIC_PUBLIC_INTERFACE(Product) + typedef typename ei_traits<Product>::LhsNested LhsNested; + typedef typename ei_traits<Product>::RhsNested RhsNested; + typedef typename ei_traits<Product>::_LhsNested _LhsNested; + typedef typename ei_traits<Product>::_RhsNested _RhsNested; + + enum { + PacketSize = ei_packet_traits<Scalar>::size + }; + + inline Product(const Lhs& lhs, const Rhs& rhs) + : m_lhs(lhs), m_rhs(rhs) + { + ei_assert(lhs.cols() == rhs.rows()); + } + + private: + + inline int _rows() const { return m_lhs.rows(); } + inline int _cols() const { return m_rhs.cols(); } + + const Scalar _coeff(int row, int col) const + { + int unique = ((Rhs::Flags&Diagonal)==Diagonal) ? col : row; + return m_lhs.coeff(row, unique) * m_rhs.coeff(unique, col); + } + + template<int LoadMode> + const PacketScalar _packetCoeff(int row, int col) const + { + if ((Rhs::Flags&Diagonal)==Diagonal) + { + ei_assert((_LhsNested::Flags&RowMajorBit)==0); + return ei_pmul(m_lhs.template packetCoeff<LoadMode>(row, col), ei_pset1(m_rhs.coeff(col, col))); + } + else + { + ei_assert(_RhsNested::Flags&RowMajorBit); + return ei_pmul(ei_pset1(m_lhs.coeff(row, row)), m_rhs.template packetCoeff<LoadMode>(row, col)); + } + } + + protected: + const LhsNested m_lhs; + const RhsNested m_rhs; +}; + +#endif // EIGEN_DIAGONALPRODUCT_H diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index 6c653a558..6875e3158 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -144,11 +144,12 @@ struct ei_packet_product_impl<false, Index, Dynamic, Lhs, Rhs, PacketScalar> */ template<typename Lhs, typename Rhs> struct ei_product_eval_mode { - enum{ value = Lhs::MaxRowsAtCompileTime >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD - && Rhs::MaxColsAtCompileTime >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD - && Lhs::MaxColsAtCompileTime >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD - && (Rhs::Flags&Diagonal)!=Diagonal - ? CacheFriendlyProduct : NormalProduct }; + enum{ value = ((Rhs::Flags&Diagonal)==Diagonal) || ((Lhs::Flags&Diagonal)==Diagonal) + ? DiagonalProduct + : Lhs::MaxRowsAtCompileTime >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD + && Rhs::MaxColsAtCompileTime >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD + && Lhs::MaxColsAtCompileTime >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD + ? CacheFriendlyProduct : NormalProduct }; }; template<typename T> class ei_product_eval_to_column_major @@ -265,44 +266,25 @@ template<typename Lhs, typename Rhs, int EvalMode> class Product : ei_no_assignm const Scalar _coeff(int row, int col) const { - if ((Rhs::Flags&Diagonal)==Diagonal) - { - return m_lhs.coeff(row, col) * m_rhs.coeff(col, col); - } - else if ((Lhs::Flags&Diagonal)==Diagonal) - { - return m_lhs.coeff(row, row) * m_rhs.coeff(row, col); - } - else - { - Scalar res; - const bool unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT; - ei_product_impl<unroll ? Lhs::ColsAtCompileTime-1 : Dynamic, - unroll ? Lhs::ColsAtCompileTime : Dynamic, - _LhsNested, _RhsNested> - ::run(row, col, m_lhs, m_rhs, res); - return res; - } + Scalar res; + const bool unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT; + ei_product_impl<unroll ? Lhs::ColsAtCompileTime-1 : Dynamic, + unroll ? Lhs::ColsAtCompileTime : Dynamic, + _LhsNested, _RhsNested> + ::run(row, col, m_lhs, m_rhs, res); + return res; } template<int LoadMode> const PacketScalar _packetCoeff(int row, int col) const { - if ((Rhs::Flags&Diagonal)==Diagonal) - { - ei_assert((_LhsNested::Flags&RowMajorBit)==0); - return ei_pmul(m_lhs.template packetCoeff<LoadMode>(row, col), ei_pset1(m_rhs.coeff(col, col))); - } - else - { - const bool unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT; - PacketScalar res; - ei_packet_product_impl<Flags&RowMajorBit ? true : false, Lhs::ColsAtCompileTime-1, - unroll ? Lhs::ColsAtCompileTime : Dynamic, - _LhsNested, _RhsNested, PacketScalar> - ::run(row, col, m_lhs, m_rhs, res); - return res; - } + const bool unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT; + PacketScalar res; + ei_packet_product_impl<Flags&RowMajorBit ? true : false, Lhs::ColsAtCompileTime-1, + unroll ? Lhs::ColsAtCompileTime : Dynamic, + _LhsNested, _RhsNested, PacketScalar> + ::run(row, col, m_lhs, m_rhs, res); + return res; } template<typename Lhs_, typename Rhs_, int EvalMode_, typename DestDerived_, bool DirectAccess_> diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index 1a4843973..23749d67c 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -52,8 +52,6 @@ struct ei_traits<Transpose<MatrixType> > & ~( Like1DArrayBit | LowerTriangularBit | UpperTriangularBit) | (int(_MatrixTypeNested::Flags)&UpperTriangularBit ? LowerTriangularBit : 0) | (int(_MatrixTypeNested::Flags)&LowerTriangularBit ? UpperTriangularBit : 0), - // Note that the above test cannot be simplified using ^ because a diagonal matrix - // has both LowerTriangularBit and UpperTriangularBit and both must be preserved. CoeffReadCost = _MatrixTypeNested::CoeffReadCost }; }; diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index fdde9acd3..163832394 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -80,7 +80,7 @@ const unsigned int Like1DArrayBit = 0x20; /** \ingroup flags * * means all diagonal coefficients are equal to 0 */ -const unsigned int ZeroDiagBit = 0x40; +const unsigned int ZeroDiagBit = 0x40; /** \ingroup flags * @@ -140,7 +140,7 @@ enum { Aligned=0, UnAligned=1 }; enum { ConditionalJumpCost = 5 }; enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight }; enum DirectionType { Vertical, Horizontal }; -enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, LazyProduct}; +enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, DiagonalProduct, LazyProduct}; #endif // EIGEN_CONSTANTS_H diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h index 80f02e71c..f6a525ead 100644 --- a/Eigen/src/Geometry/Transform.h +++ b/Eigen/src/Geometry/Transform.h @@ -25,6 +25,94 @@ #ifndef EIGEN_TRANSFORM_H #define EIGEN_TRANSFORM_H +/** \class Orientation2D + * + * \brief Represents an orientation/rotation 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 representating the rotation angle + * in radian with some additional features such as the conversion from/to + * rotation matrix. Moreover this class aims to provide a similar interface + * to Quaternion in order to facilitate the writting of generic algorithm + * dealing with rotations. + * + * \sa class Quaternion, class Transform + */ +template<typename _Scalar> +class Orientation2D +{ +public: + enum { Dim = 2 }; + /** the scalar type of the coefficients */ + typedef _Scalar Scalar; + typedef Matrix<Scalar,2,2> Matrix2; + +protected: + + Scalar m_angle; + +public: + + inline Orientation2D(Scalar a) : m_angle(a) {} + inline operator Scalar& () { return m_angle; } + inline operator Scalar () const { return m_angle; } + + template<typename Derived> + Orientation2D& fromRotationMatrix(const MatrixBase<Derived>& m); + Matrix2 toRotationMatrix(void) const; + + Orientation2D slerp(Scalar t, const Orientation2D& other) const; +}; + +/** returns the default type used to represent an orientation. + */ +template<typename Scalar, int Dim> +struct ei_get_orientation_type; + +template<typename Scalar> +struct ei_get_orientation_type<Scalar,2> +{ typedef Orientation2D<Scalar> type; }; + +template<typename Scalar> +struct ei_get_orientation_type<Scalar,3> +{ typedef Quaternion<Scalar> type; }; + +/** 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> +Orientation2D<Scalar>& Orientation2D<Scalar>::fromRotationMatrix(const MatrixBase<Derived>& mat) +{ + EIGEN_STATIC_ASSERT(Derived::RowsAtCompileTime==2 && Derived::ColsAtCompileTime==2,you_did_a_programming_error); + 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 Orientation2D<Scalar>::Matrix2 +Orientation2D<Scalar>::toRotationMatrix(void) const +{ + Scalar sinA = ei_sin(m_angle); + Scalar cosA = ei_cos(m_angle); + return Matrix2(cosA, -sinA, sinA, cosA); +} + +/** \returns the spherical interpolation between \c *this and \a other using + * parameter \a t. It is equivalent to a linear interpolation. + */ +template<typename Scalar> +Orientation2D<Scalar> +Orientation2D<Scalar>::slerp(Scalar t, const Orientation2D& other) const +{ + return m_angle * (1-t) + t * other; +} + + /** \class Transform * * \brief Represents an homogeneous transformation in a N dimensional space @@ -32,6 +120,11 @@ * \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 are available if the preprocessor token + * EIGEN_QT_SUPPORT is defined. * */ template<typename _Scalar, int _Dim> @@ -47,6 +140,7 @@ public: typedef Block<MatrixType,Dim,Dim> AffineMatrixRef; typedef Matrix<Scalar,Dim,1> VectorType; typedef Block<MatrixType,Dim,1> VectorRef; + typedef typename ei_get_orientation_type<Scalar,Dim>::type Orientation; protected: @@ -59,13 +153,42 @@ protected: public: + /** Default constructor without initialization of the coefficients. */ + Transform() { } + + inline Transform(const Transform& other) + { m_matrix = other.m_matrix; } + + inline Transform& operator=(const Transform& other) + { m_matrix = other.m_matrix; } + + template<typename OtherDerived> + inline explicit Transform(const MatrixBase<OtherDerived>& other) + { m_matrix = other; } + + template<typename OtherDerived> + inline Transform& operator=(const MatrixBase<OtherDerived>& other) + { m_matrix = other; } + + #ifdef EIGEN_QT_SUPPORT + inline Transform(const QMatrix& other); + inline Transform& operator=(const QMatrix& other); + inline QMatrix toQMatrix(void) const; + #endif + + /** \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 affine (linear) part of the transformation */ inline const AffineMatrixRef affine() const { return m_matrix.template block<Dim,Dim>(0,0); } + /** \returns a writable expression of the affine (linear) part of the transformation */ inline AffineMatrixRef affine() { return m_matrix.template block<Dim,Dim>(0,0); } + /** \returns a read-only expression of the translation vector of the transformation */ inline const VectorRef translation() const { return m_matrix.template block<Dim,1>(0,Dim); } + /** \returns a writable expression of the translation vector of the transformation */ inline VectorRef translation() { return m_matrix.template block<Dim,1>(0,Dim); } template<typename OtherDerived> @@ -78,6 +201,11 @@ public: const typename ProductReturnType<OtherDerived>::Type operator * (const MatrixBase<OtherDerived> &other) const; + /** Contatenates two transformations */ + Product<MatrixType,MatrixType> + operator * (const Transform& other) const + { return m_matrix * other.matrix(); } + void setIdentity() { m_matrix.setIdentity(); } template<typename OtherDerived> @@ -92,13 +220,58 @@ public: template<typename OtherDerived> Transform& pretranslate(const MatrixBase<OtherDerived> &other); + template<typename OtherDerived> + Transform& shear(Scalar sx, Scalar sy); + + template<typename OtherDerived> + Transform& preshear(Scalar sx, Scalar sy); + AffineMatrixType extractRotation() const; AffineMatrixType extractRotationNoShear() const; + template<typename PositionDerived, typename ScaleDerived> + Transform& fromPositionOrientationScale(const MatrixBase<PositionDerived> &position, + const Orientation& orientation, const MatrixBase<ScaleDerived> &scale); + + const Inverse<MatrixType, false> inverse() const + { return m_matrix.inverse(); } + protected: }; +#ifdef EIGEN_QT_SUPPORT +/** Initialises \c *this from a QMatrix assuming the dimension is 2. + */ +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. + */ +template<typename Scalar, int Dim> +Transform<Scalar,Dim>& Transform<Scalar,Dim>::operator=(const QMatrix& other) +{ + EIGEN_STATIC_ASSERT(Dim==2, you_did_a_programming_error); + m_matrix << other.m11(), other.m21(), other.dx(), + other.m12(), other.m22(), other.dy(), + 0, 0, 1; +} + +/** \returns a QMatrix from \c *this assuming the dimension is 2. + */ +template<typename Scalar, int Dim> +QMatrix Transform<Scalar,Dim>::toQMatrix(void) const +{ + EIGEN_STATIC_ASSERT(Dim==2, you_did_a_programming_error); + return QMatrix( other.coeffRef(0,0), other.coeffRef(1,0), + other.coeffRef(0,1), other.coeffRef(1,1), + other.coeffRef(0,2), other.coeffRef(1,2), +} +#endif + template<typename Scalar, int Dim> template<typename OtherDerived> const typename Transform<Scalar,Dim>::template ProductReturnType<OtherDerived>::Type @@ -133,7 +306,7 @@ Transform<Scalar,Dim>::prescale(const MatrixBase<OtherDerived> &other) { EIGEN_STATIC_ASSERT(int(OtherDerived::IsVectorAtCompileTime) && int(OtherDerived::SizeAtCompileTime)==int(Dim), you_did_a_programming_error); - m_matrix.template block<3,4>(0,0) = (other.asDiagonal().eval() * m_matrix.template block<3,4>(0,0)).lazy(); + m_matrix.template block<3,4>(0,0) = (other.asDiagonal() * m_matrix.template block<3,4>(0,0)).lazy(); return *this; } @@ -167,6 +340,39 @@ Transform<Scalar,Dim>::pretranslate(const MatrixBase<OtherDerived> &other) 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> +template<typename OtherDerived> +Transform<Scalar,Dim>& +Transform<Scalar,Dim>::shear(Scalar sx, Scalar sy) +{ + EIGEN_STATIC_ASSERT(int(OtherDerived::IsVectorAtCompileTime) + && int(OtherDerived::SizeAtCompileTime)==int(Dim) && int(Dim)==2, you_did_a_programming_error); + VectorType tmp = affine().col(0)*sy + affine().col(1); + affine() << affine().col(0) + affine().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> +template<typename OtherDerived> +Transform<Scalar,Dim>& +Transform<Scalar,Dim>::preshear(Scalar sx, Scalar sy) +{ + EIGEN_STATIC_ASSERT(int(OtherDerived::IsVectorAtCompileTime) + && int(OtherDerived::SizeAtCompileTime)==int(Dim), you_did_a_programming_error); + m_matrix.template block<3,4>(0,0) = AffineMatrixType(1, sx, sy, 1) * m_matrix.template block<3,4>(0,0); + return *this; +} + /** \returns the rotation part of the transformation using a QR decomposition. * \sa extractRotationNoShear() */ @@ -189,6 +395,22 @@ Transform<Scalar,Dim>::extractRotationNoShear() const .verticalRedux(ei_scalar_sum_op<Scalar>()).cwiseSqrt(); } +/** 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 ScaleDerived> +Transform<Scalar,Dim>& +Transform<Scalar,Dim>::fromPositionOrientationScale(const MatrixBase<PositionDerived> &position, + const Orientation& orientation, const MatrixBase<ScaleDerived> &scale) +{ + affine() = orientation.toRotationMatrix(); + translation() = position; + m_matrix(Dim,Dim) = 1.; + m_matrix.template block<1,Dim>(Dim,0).setZero(); + affine() *= scale.asDiagonal(); +} + //---------- template<typename Scalar, int Dim> @@ -216,12 +438,17 @@ template<typename Other> struct Transform<Scalar,Dim>::ei_transform_product_impl<Other,Dim,1> { typedef typename Transform<Scalar,Dim>::AffineMatrixRef MatrixType; - typedef const CwiseBinaryOp< - ei_scalar_sum_op<Scalar>, - NestByValue<Product<NestByValue<MatrixType>,Other> >, - NestByValue<typename Transform<Scalar,Dim>::VectorRef> > ResultType; + typedef const CwiseUnaryOp< + ei_scalar_multiple_op<Scalar>, + NestByValue<CwiseBinaryOp< + ei_scalar_sum_op<Scalar>, + NestByValue<Product<NestByValue<MatrixType>,Other> >, + NestByValue<typename Transform<Scalar,Dim>::VectorRef> > > + > ResultType; + // FIXME shall we offer an optimized version when the last row is know to be 0,0...,0,1 ? static ResultType run(const Transform<Scalar,Dim>& tr, const Other& other) - { return (tr.affine().nestByValue() * other).nestByValue() + tr.translation().nestByValue(); } + { return ((tr.affine().nestByValue() * other).nestByValue() + tr.translation().nestByValue()).nestByValue() + * (Scalar(1) / ( (tr.matrix().template block<1,Dim>(Dim,0) * other).coeff(0) + tr.matrix().coeff(Dim,Dim))); } }; #endif // EIGEN_TRANSFORM_H diff --git a/Eigen/src/QR/HessenbergDecomposition.h b/Eigen/src/QR/HessenbergDecomposition.h index 8f4710993..880997563 100755 --- a/Eigen/src/QR/HessenbergDecomposition.h +++ b/Eigen/src/QR/HessenbergDecomposition.h @@ -239,7 +239,8 @@ HessenbergDecomposition<MatrixType>::matrixH(void) const // and fill it (to avoid temporaries) int n = m_matrix.rows(); MatrixType matH = m_matrix; - matH.corner(BottomLeft,n-2, n-2).template part<Lower>().setZero(); + if (n>2) + matH.corner(BottomLeft,n-2, n-2).template part<Lower>().setZero(); return matH; } diff --git a/Eigen/src/QR/Tridiagonalization.h b/Eigen/src/QR/Tridiagonalization.h index 1473b5bfa..bcb4a2b59 100755 --- a/Eigen/src/QR/Tridiagonalization.h +++ b/Eigen/src/QR/Tridiagonalization.h @@ -160,8 +160,11 @@ Tridiagonalization<MatrixType>::matrixT(void) const int n = m_matrix.rows(); MatrixType matT = m_matrix; matT.corner(TopRight,n-1, n-1).diagonal() = subDiagonal().conjugate(); - matT.corner(TopRight,n-2, n-2).template part<Upper>().setZero(); - matT.corner(BottomLeft,n-2, n-2).template part<Lower>().setZero(); + if (n>2) + { + matT.corner(TopRight,n-2, n-2).template part<Upper>().setZero(); + matT.corner(BottomLeft,n-2, n-2).template part<Lower>().setZero(); + } return matT; } diff --git a/test/geometry.cpp b/test/geometry.cpp index 9d5a07af8..f41d086de 100644 --- a/test/geometry.cpp +++ b/test/geometry.cpp @@ -116,6 +116,9 @@ template<typename Scalar> void geometry(void) t1.translate(-v0); VERIFY((t0.matrix() * t1.matrix()).isIdentity()); + + t1.fromPositionOrientationScale(v0, q1, v1); + VERIFY_IS_APPROX(t1.matrix(), t0.matrix()); } void test_geometry() |