Add .perpendicular() function in Geometry module (adapted from Eigen1)

Documentation:
 * add an overview for each module.
 * add an example for .all() and Cwise::operator<

Eigen/Array

@@ -5,6 +5,23 @@
 
 namespace Eigen {
 
+/** \defgroup Array Array module
+  * This module provides several handy features to manipulate matrices as simple array of values.
+  * In addition to listed classes, it defines various methods of the Cwise interface
+  * (accessible from MatrixBase::cwise()), including:
+  *  - matrix-scalar sum,
+  *  - coeff-wise comparison operators,
+  *  - sin, cos, sqrt, pow, exp, log, square, cube, reciprocal.
+  *
+  * This module also provides various MatrixBase methods, including:
+  *  - \ref MatrixBase::all() "all", \ref MatrixBase::any() "any",
+  *  - \ref MatrixBase::Random() "random matrix initialization"
+  *
+  * \code
+  * #include <Eigen/Array>
+  * \endcode
+  */
+
 #include "src/Array/CwiseOperators.h"
 #include "src/Array/Functors.h"
 #include "src/Array/AllAndAny.h"

Eigen/Cholesky

@@ -14,6 +14,17 @@
 
 namespace Eigen {
 
+/** \defgroup Cholesky_Module Cholesky module
+  * This module provides two variants of the Cholesky decomposition for selfadjoint (hermitian) matrices.
+  * Those decompositions are accessible via the following MatrixBase methods:
+  *  - MatrixBase::cholesky(),
+  *  - MatrixBase::choleskyNoSqrt()
+  *
+  * \code
+  * #include <Eigen/Cholesky>
+  * \endcode
+  */
+
 #include "src/Cholesky/Cholesky.h"
 #include "src/Cholesky/CholeskyWithoutSquareRoot.h"
 

Eigen/Geometry

@@ -5,7 +5,7 @@
 
 namespace Eigen {
 
-/** \defgroup Geometry
+/** \defgroup Geometry Geometry module
   * This module provides support for:
   *  - fixed-size homogeneous transformations
   *  - 2D and 3D rotations

Eigen/LU

@@ -5,6 +5,17 @@
 
 namespace Eigen {
 
+/** \defgroup LU_Module LU module
+  * This module includes LU decomposition and related notions such as matrix inversion and determinant.
+  * This module defines the following MatrixBase methods:
+  *  - MatrixBase::inverse()
+  *  - MatrixBase::determinant()
+  *
+  * \code
+  * #include <Eigen/LU>
+  * \endcode
+  */
+
 #include "src/LU/Determinant.h"
 #include "src/LU/Inverse.h"
 

Eigen/QR

@@ -15,6 +15,18 @@
 
 namespace Eigen {
 
+/** \defgroup QR_Module QR module
+  * This module mainly provides QR decomposition and an eigen value solver.
+  * This module also provides some MatrixBase methods, including:
+  *  - MatrixBase::qr(),
+  *  - MatrixBase::eigenvalues(),
+  *  - MatrixBase::matrixNorm()
+  *
+  * \code
+  * #include <Eigen/QR>
+  * \endcode
+  */
+
 #include "src/QR/QR.h"
 #include "src/QR/Tridiagonalization.h"
 #include "src/QR/EigenSolver.h"

Eigen/src/Array/AllAndAny.h

@@ -81,7 +81,12 @@ struct ei_any_unroller<Derived, Dynamic>
   * 
   * \returns true if all coefficients are true
   *
-  * \sa MatrixBase::any()
+  * \addexample CwiseAll \label How to check whether a point is inside a box (using operator< and all())
+  *
+  * Example: \include MatrixBase_all.cpp
+  * Output: \verbinclude MatrixBase_all.out
+  *
+  * \sa MatrixBase::any(), Cwise::operator<()
   */
 template<typename Derived>
 bool MatrixBase<Derived>::all(void) const
@@ -105,7 +110,7 @@ bool MatrixBase<Derived>::all(void) const
   * 
   * \returns true if at least one coefficient is true
   *
-  * \sa MatrixBase::any()
+  * \sa MatrixBase::all()
   */
 template<typename Derived>
 bool MatrixBase<Derived>::any(void) const

Eigen/src/Array/CwiseOperators.h

@@ -127,6 +127,8 @@ Cwise<ExpressionType>::cube() const
   * 
   * \returns an expression of the coefficient-wise \< operator of *this and \a other
   *
+  * See MatrixBase::all() for an example.
+  *
   * \sa class CwiseBinaryOp
   */
 template<typename ExpressionType>
@@ -155,6 +157,8 @@ Cwise<ExpressionType>::operator<=(const MatrixBase<OtherDerived> &other) const
   * 
   * \returns an expression of the coefficient-wise \> operator of *this and \a other
   *
+  * See MatrixBase::all() for an example.
+  *
   * \sa class CwiseBinaryOp
   */
 template<typename ExpressionType>
@@ -208,7 +212,8 @@ Cwise<ExpressionType>::operator!=(const MatrixBase<OtherDerived> &other) const
 }
 
 
-/** \returns an expression of \c *this with each coeff incremented by the constant \a scalar */
+/** \array_module
+  * \returns an expression of \c *this with each coeff incremented by the constant \a scalar */
 template<typename ExpressionType>
 inline const typename Cwise<ExpressionType>::ScalarAddReturnType
 Cwise<ExpressionType>::operator+(const Scalar& scalar) const
@@ -216,14 +221,16 @@ Cwise<ExpressionType>::operator+(const Scalar& scalar) const
   return typename Cwise<ExpressionType>::ScalarAddReturnType(m_matrix, ei_scalar_add_op<Scalar>(scalar));
 }
 
-/** \see operator+ */
+/** \array_module
+  * \see operator+() */
 template<typename ExpressionType>
 inline ExpressionType& Cwise<ExpressionType>::operator+=(const Scalar& scalar)
 {
   return m_matrix.const_cast_derived() = *this + scalar;
 }
 
-/** \returns an expression of \c *this with each coeff decremented by the constant \a scalar */
+/** \array_module
+  * \returns an expression of \c *this with each coeff decremented by the constant \a scalar */
 template<typename ExpressionType>
 inline const typename Cwise<ExpressionType>::ScalarAddReturnType
 Cwise<ExpressionType>::operator-(const Scalar& scalar) const
@@ -231,7 +238,8 @@ Cwise<ExpressionType>::operator-(const Scalar& scalar) const
   return *this + (-scalar);
 }
 
-/** \see operator- */
+/** \array_module
+  * \see operator- */
 template<typename ExpressionType>
 inline ExpressionType& Cwise<ExpressionType>::operator-=(const Scalar& scalar)
 {

Eigen/src/Array/PartialRedux.h

@@ -26,7 +26,7 @@
 #ifndef EIGEN_PARTIAL_REDUX_H
 #define EIGEN_PARTIAL_REDUX_H
 
-/** \array_module
+/** \array_module \ingroup Array
   * 
   * \class PartialReduxExpr
   *
@@ -128,7 +128,7 @@ struct ei_member_redux {
   const BinaryOp m_functor;
 };
 
-/** \array_module
+/** \array_module \ingroup Array
   *
   * \class PartialRedux
   *
@@ -141,7 +141,7 @@ struct ei_member_redux {
   * It is the return type of MatrixBase::colwise() and MatrixBase::rowwise()
   * and most of the time this is the only way it is used.
   *
-  * \sa MatrixBase::colwise(), MatrixBase::rowwise()
+  * \sa MatrixBase::colwise(), MatrixBase::rowwise(), class PartialReduxExpr
   */
 template<typename ExpressionType, int Direction> class PartialRedux
 {
@@ -236,9 +236,7 @@ MatrixBase<Derived>::rowwise() const
   return derived();
 }
 
-/** \array_module
-  * 
-  * \returns a row or column vector expression of \c *this reduxed by \a func
+/** \returns a row or column vector expression of \c *this reduxed by \a func
   *
   * The template parameter \a BinaryOp is the type of the functor
   * of the custom redux operator. Note that func must be an associative operator.

Eigen/src/Array/Random.h

@@ -49,7 +49,7 @@ struct ei_functor_traits<ei_scalar_random_op<Scalar> >
   * Example: \include MatrixBase_random_int_int.cpp
   * Output: \verbinclude MatrixBase_random_int_int.out
   *
-  * \sa ei_random(), ei_random(int)
+  * \sa MatrixBase::setRandom(), MatrixBase::Random(int), MatrixBase::Random()
   */
 template<typename Derived>
 inline const CwiseNullaryOp<ei_scalar_random_op<typename ei_traits<Derived>::Scalar>, Derived>
@@ -74,7 +74,7 @@ MatrixBase<Derived>::Random(int rows, int cols)
   * Example: \include MatrixBase_random_int.cpp
   * Output: \verbinclude MatrixBase_random_int.out
   *
-  * \sa ei_random(), ei_random(int,int)
+  * \sa MatrixBase::setRandom(), MatrixBase::Random(int,int), MatrixBase::Random()
   */
 template<typename Derived>
 inline const CwiseNullaryOp<ei_scalar_random_op<typename ei_traits<Derived>::Scalar>, Derived>
@@ -94,7 +94,7 @@ MatrixBase<Derived>::Random(int size)
   * Example: \include MatrixBase_random.cpp
   * Output: \verbinclude MatrixBase_random.out
   *
-  * \sa ei_random(int), ei_random(int,int)
+  * \sa MatrixBase::setRandom(), MatrixBase::Random(int,int), MatrixBase::Random(int)
   */
 template<typename Derived>
 inline const CwiseNullaryOp<ei_scalar_random_op<typename ei_traits<Derived>::Scalar>, Derived>
@@ -110,7 +110,7 @@ MatrixBase<Derived>::Random()
   * Example: \include MatrixBase_setRandom.cpp
   * Output: \verbinclude MatrixBase_setRandom.out
   *
-  * \sa class CwiseNullaryOp, ei_random()
+  * \sa class CwiseNullaryOp, MatrixBase::setRandom(int,int)
   */
 template<typename Derived>
 inline Derived& MatrixBase<Derived>::setRandom()

Eigen/src/Cholesky/Cholesky.h

@@ -25,7 +25,9 @@
 #ifndef EIGEN_CHOLESKY_H
 #define EIGEN_CHOLESKY_H
 
-/** \class Cholesky
+/** \ingroup Cholesky_Module
+  *
+  * \class Cholesky
   *
   * \brief Standard Cholesky decomposition of a matrix and associated features
   *

Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h

@@ -25,7 +25,9 @@
 #ifndef EIGEN_CHOLESKY_WITHOUT_SQUARE_ROOT_H
 #define EIGEN_CHOLESKY_WITHOUT_SQUARE_ROOT_H
 
-/** \class CholeskyWithoutSquareRoot
+/** \ingroup Cholesky_Module
+  *
+  * \class CholeskyWithoutSquareRoot
   *
   * \brief Robust Cholesky decomposition of a matrix and associated features
   *

Eigen/src/Core/Cwise.h

@@ -45,6 +45,8 @@ template<typename Scalar, bool PacketAccess = (int(ei_packet_traits<Scalar>::siz
   * It is the return type of MatrixBase::cwise()
   * and most of the time this is the only way it is used.
   *
+  * Note that some methods are defined in the \ref Array module.
+  *
   * \sa MatrixBase::cwise()
   */
 template<typename ExpressionType> class Cwise

Eigen/src/Core/MatrixBase.h

@@ -31,7 +31,9 @@
   *
   * This class is the base that is inherited by all matrix, vector, and expression
   * types. Most of the Eigen API is contained in this class. Other important classes for
-  * the Eigen API are Matrix, Cwise, and Part.
+  * the Eigen API are Matrix, Cwise, and PartialRedux.
+  * 
+  * Note that some methods are defined in the \ref Array module.
   *
   * \param Derived is the derived type, e.g. a matrix type, or an expression, etc.
   *
@@ -48,7 +50,6 @@
     }
   * \endcode
   *
-  * \nosubgrouping
   */
 template<typename Derived> class MatrixBase
 {
@@ -549,7 +550,7 @@ template<typename Derived> class MatrixBase
     template<typename OtherDerived>
     typename ei_eval<Derived>::type
     cross(const MatrixBase<OtherDerived>& other) const;
-
+    typename ei_eval<Derived>::type perpendicular(void) const;
 };
 
 #endif // EIGEN_MATRIXBASE_H

Eigen/src/Core/arch/SSE/PacketMath.h

@@ -61,7 +61,7 @@ template<> inline __m128i ei_pmul(const __m128i& a, const __m128i& b)
 
 template<> inline __m128  ei_pdiv(const __m128&  a, const __m128&  b) { return _mm_div_ps(a,b); }
 template<> inline __m128d ei_pdiv(const __m128d& a, const __m128d& b) { return _mm_div_pd(a,b); }
-template<> inline __m128i ei_pdiv(const __m128i& a, const __m128i& b)
+template<> inline __m128i ei_pdiv(const __m128i& /*a*/, const __m128i& /*b*/)
 { ei_assert(false && "packet integer division are not supported by SSE"); }
 
 // for some weird raisons, it has to be overloaded for packet integer

Eigen/src/Geometry/Cross.h

@@ -2,6 +2,7 @@
 // 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
@@ -32,6 +33,8 @@ template<typename OtherDerived>
 inline typename ei_eval<Derived>::type
 MatrixBase<Derived>::cross(const MatrixBase<OtherDerived>& other) const
 {
+  EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Derived,3);
+  
   // Note that there is no need for an expression here since the compiler
   // optimize such a small temporary very well (even within a complex expression)
   const typename ei_nested<Derived,2>::type lhs(derived());
@@ -43,4 +46,65 @@ MatrixBase<Derived>::cross(const MatrixBase<OtherDerived>& other) const
   );
 }
 
+template<typename Derived, int Size = Derived::SizeAtCompileTime>
+struct ei_perpendicular_selector
+{
+  typedef typename ei_eval<Derived>::type VectorType;
+  typedef typename ei_traits<Derived>::Scalar Scalar;
+  inline static VectorType run(const Derived& src)
+  {
+    VectorType perp;
+    /* Let us compute the crossed product of *this with a vector
+     * that is not too close to being colinear to *this.
+     */
+
+    /* unless the x and y coords are both close to zero, we can
+     * simply take ( -y, x, 0 ) and normalize it.
+     */
+    if((!ei_isMuchSmallerThan(src.x(), src.z()))
+    || (!ei_isMuchSmallerThan(src.y(), src.z())))
+    {
+      Scalar invnm = Scalar(1)/src.template start<2>().norm();
+      perp.template start<3>() << -ei_conj(src.y())*invnm, ei_conj(src.x())*invnm, 0;
+    }
+    /* if both x and y are close to zero, then the vector is close
+     * to the z-axis, so it's far from colinear to the x-axis for instance.
+     * So we take the crossed product with (1,0,0) and normalize it.
+     */
+    else
+    {
+      Scalar invnm = Scalar(1)/src.template end<2>().norm();
+      perp.template start<3>() << 0, -ei_conj(src.z())*invnm, ei_conj(src.y())*invnm;
+    }
+    if (Derived::SizeAtCompileTime>3
+    || (Derived::SizeAtCompileTime==Dynamic && src.size()>3))
+      perp.end(src.size()-3).setZero();
+
+    return perp;
+   }
+};
+
+template<typename Derived>
+struct ei_perpendicular_selector<Derived,2>
+{
+  typedef typename ei_eval<Derived>::type VectorType;
+  inline static VectorType run(const Derived& src)
+  { return VectorType(-ei_conj(src.y()), ei_conj(src.x())).normalized(); }
+};
+
+/** \Returns an orthogonal vector of \c *this
+  *
+  * The size of \c *this must be at least 2. If the size is exactly 2,
+  * then the returned vector is a counter clock wise rotation of \c *this, \ie (-y,x).
+  *
+  * \sa cross()
+  */
+template<typename Derived>
+typename ei_eval<Derived>::type
+MatrixBase<Derived>::perpendicular() const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived);
+  return ei_perpendicular_selector<Derived>::run(derived());
+}
+
 #endif // EIGEN_CROSS_H

Eigen/src/QR/EigenSolver.h

@@ -25,7 +25,9 @@
 #ifndef EIGEN_EIGENSOLVER_H
 #define EIGEN_EIGENSOLVER_H
 
-/** \class EigenSolver
+/** \ingroup QR_Module
+  *
+  * \class EigenSolver
   *
   * \brief Eigen values/vectors solver for non selfadjoint matrices
   *

Eigen/src/QR/HessenbergDecomposition.h

@@ -25,7 +25,9 @@
 #ifndef EIGEN_HESSENBERGDECOMPOSITION_H
 #define EIGEN_HESSENBERGDECOMPOSITION_H
 
-/** \class HessenbergDecomposition
+/** \ingroup QR_Module
+  *
+  * \class HessenbergDecomposition
   *
   * \brief Reduces a squared matrix to an Hessemberg form
   *

Eigen/src/QR/QR.h

@@ -25,7 +25,9 @@
 #ifndef EIGEN_QR_H
 #define EIGEN_QR_H
 
-/** \class QR
+/** \ingroup QR_Module
+  *
+  * \class QR
   *
   * \brief QR decomposition of a matrix
   *

Eigen/src/QR/SelfAdjointEigenSolver.h

@@ -25,7 +25,7 @@
 #ifndef EIGEN_SELFADJOINTEIGENSOLVER_H
 #define EIGEN_SELFADJOINTEIGENSOLVER_H
 
-/** \qr_module
+/** \qr_module \ingroup QR_Module
   *
   * \class SelfAdjointEigenSolver
   *

Eigen/src/QR/Tridiagonalization.h

@@ -25,7 +25,9 @@
 #ifndef EIGEN_TRIDIAGONALIZATION_H
 #define EIGEN_TRIDIAGONALIZATION_H
 
-/** \class Tridiagonalization
+/** \ingroup QR_Module
+  *
+  * \class Tridiagonalization
   *
   * \brief Trigiagonal decomposition of a selfadjoint matrix
   *

doc/snippets/MatrixBase_all.cpp

@@ -0,0 +1,7 @@
+Vector3f boxMin(Vector3f::Zero()), boxMax(Vector3f::Ones());
+Vector3f p0 = Vector3f::Random(), p1 = Vector3f::Random().cwise().abs();
+// let's check if p0 and p1 are inside the axis aligned box defined by the corners boxMin,boxMax:
+cout << "Is (" << p0.transpose() << ") inside the box: "
+     << ((boxMin.cwise()<p0).all() && (boxMax.cwise()>p0).all()) << endl;
+cout << "Is (" << p1.transpose() << ") inside the box: "
+     << ((boxMin.cwise()<p1).all() && (boxMax.cwise()>p1).all()) << endl;

test/geometry.cpp

@@ -45,10 +45,23 @@ template<typename Scalar> void geometry(void)
   Vector3 v0 = Vector3::Random(),
     v1 = Vector3::Random(),
     v2 = Vector3::Random();
+  Vector2 u0 = Vector2::Random();
   Matrix3 matrot1;
 
   Scalar a = ei_random<Scalar>(-M_PI, M_PI);
 
+  // cross product
+  VERIFY_IS_MUCH_SMALLER_THAN(v1.cross(v2).dot(v1), Scalar(1));
+  Matrix3 m;
+  m << v0.normalized(),
+      (v0.cross(v1)).normalized(),
+      (v0.cross(v1).cross(v0)).normalized();
+  VERIFY(m.isUnitary());
+
+  // perpendicular
+  VERIFY_IS_MUCH_SMALLER_THAN(u0.perpendicular().dot(u0), Scalar(1));
+  VERIFY_IS_MUCH_SMALLER_THAN(v0.perpendicular().dot(v0), Scalar(1));
+
   q1 = AngleAxis(ei_random<Scalar>(-M_PI, M_PI), v0.normalized());
   q2 = AngleAxis(ei_random<Scalar>(-M_PI, M_PI), v1.normalized());
 
@@ -82,14 +95,6 @@ template<typename Scalar> void geometry(void)
   VERIFY_IS_APPROX(q1 * (q1.inverse() * v1), v1);
   VERIFY_IS_APPROX(q1 * (q1.conjugate() * v1), v1);
 
-  // cross product
-  VERIFY_IS_MUCH_SMALLER_THAN(v1.cross(v2).dot(v1), Scalar(1));
-  Matrix3 m;
-  m << v0.normalized(),
-      (v0.cross(v1)).normalized(),
-      (v0.cross(v1).cross(v0)).normalized();
-  VERIFY(m.isUnitary());
-
   // AngleAxis
   VERIFY_IS_APPROX(AngleAxis(a,v1.normalized()).toRotationMatrix(),
     Quaternion(AngleAxis(a,v1.normalized())).toRotationMatrix());