diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index d87c6b30b..ea2178fe3 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -387,20 +387,9 @@ template class MatrixBase /////////// Geometry module /////////// - #ifndef EIGEN_PARSED_BY_DOXYGEN - /// \internal helper struct to form the return type of the cross product - template struct cross_product_return_type { - typedef typename ScalarBinaryOpTraits::Scalar,typename internal::traits::Scalar>::ReturnType Scalar; - typedef Matrix type; - }; - #endif // EIGEN_PARSED_BY_DOXYGEN template EIGEN_DEVICE_FUNC -#ifndef EIGEN_PARSED_BY_DOXYGEN - inline typename cross_product_return_type::type -#else - inline PlainObject -#endif + inline typename internal::cross_impl::return_type cross(const MatrixBase& other) const; template diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 503d6511d..8f87c4a2c 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -270,8 +270,10 @@ template class Ho template class JacobiRotation; // Geometry module: +namespace internal { +template::SizeAtCompileTime> struct cross_impl; +} template class RotationBase; -template class Cross; template class QuaternionBase; template class Rotation2D; template class AngleAxis; diff --git a/Eigen/src/Geometry/OrthoMethods.h b/Eigen/src/Geometry/OrthoMethods.h index dc6a762f5..fbf020d30 100644 --- a/Eigen/src/Geometry/OrthoMethods.h +++ b/Eigen/src/Geometry/OrthoMethods.h @@ -15,39 +15,83 @@ namespace Eigen { +namespace internal { + +// Vector3 version (default) +template +struct cross_impl +{ + typedef typename ScalarBinaryOpTraits::Scalar,typename internal::traits::Scalar>::ReturnType Scalar; + typedef Matrix::RowsAtCompileTime,MatrixBase::ColsAtCompileTime> return_type; + + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + return_type run(const MatrixBase& first, const MatrixBase& second) + { + EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Derived,3) + EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived,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) + typename internal::nested_eval::type lhs(first.derived()); + typename internal::nested_eval::type rhs(second.derived()); + return return_type( + numext::conj(lhs.coeff(1) * rhs.coeff(2) - lhs.coeff(2) * rhs.coeff(1)), + numext::conj(lhs.coeff(2) * rhs.coeff(0) - lhs.coeff(0) * rhs.coeff(2)), + numext::conj(lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0)) + ); + } +}; + +// Vector2 version +template +struct cross_impl +{ + typedef typename ScalarBinaryOpTraits::Scalar,typename internal::traits::Scalar>::ReturnType Scalar; + typedef Scalar return_type; + + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + return_type run(const MatrixBase& first, const MatrixBase& second) + { + EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Derived,2); + EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived,2); + typename internal::nested_eval::type lhs(first.derived()); + typename internal::nested_eval::type rhs(second.derived()); + return numext::conj(lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0)); + } +}; + +} // end namespace internal + /** \geometry_module \ingroup Geometry_Module * - * \returns the cross product of \c *this and \a other + * \returns the cross product of \c *this and \a other. This is either a scalar for size-2 vectors or a size-3 vector for size-3 vectors. * - * Here is a very good explanation of cross-product: http://xkcd.com/199/ + * This method is implemented for two different cases: between vectors of fixed size 2 and between vectors of fixed size 3. * - * With complex numbers, the cross product is implemented as - * \f$ (\mathbf{a}+i\mathbf{b}) \times (\mathbf{c}+i\mathbf{d}) = (\mathbf{a} \times \mathbf{c} - \mathbf{b} \times \mathbf{d}) - i(\mathbf{a} \times \mathbf{d} - \mathbf{b} \times \mathbf{c})\f$ + * For vectors of size 3, the output is simply the traditional cross product. + * + * For vectors of size 2, the output is a scalar. + * Given vectors \f$ v = \begin{bmatrix} v_1 & v_2 \end{bmatrix} \f$ and \f$ w = \begin{bmatrix} w_1 & w_2 \end{bmatrix} \f$, + * the result is simply \f$ v\times w = \overline{v_1 w_2 - v_2 w_1} = \text{conj}\left|\begin{smallmatrix} v_1 & w_1 \\ v_2 & w_2 \end{smallmatrix}\right| \f$; + * or, to put it differently, it is the third coordinate of the cross product of \f$ \begin{bmatrix} v_1 & v_2 & v_3 \end{bmatrix} \f$ and \f$ \begin{bmatrix} w_1 & w_2 & w_3 \end{bmatrix} \f$. + * For real-valued inputs, the result can be interpreted as the signed area of a parallelogram spanned by the two vectors. + * + * \note With complex numbers, the cross product is implemented as + * \f$ (\mathbf{a}+i\mathbf{b}) \times (\mathbf{c}+i\mathbf{d}) = (\mathbf{a} \times \mathbf{c} - \mathbf{b} \times \mathbf{d}) - i(\mathbf{a} \times \mathbf{d} + \mathbf{b} \times \mathbf{c})\f$ * * \sa MatrixBase::cross3() */ template template -#ifndef EIGEN_PARSED_BY_DOXYGEN EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE -typename MatrixBase::template cross_product_return_type::type +#ifndef EIGEN_PARSED_BY_DOXYGEN +typename internal::cross_impl::return_type #else -typename MatrixBase::PlainObject +inline std::conditional_t #endif MatrixBase::cross(const MatrixBase& other) const { - EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Derived,3) - EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived,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) - typename internal::nested_eval::type lhs(derived()); - typename internal::nested_eval::type rhs(other.derived()); - return typename cross_product_return_type::type( - numext::conj(lhs.coeff(1) * rhs.coeff(2) - lhs.coeff(2) * rhs.coeff(1)), - numext::conj(lhs.coeff(2) * rhs.coeff(0) - lhs.coeff(0) * rhs.coeff(2)), - numext::conj(lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0)) - ); + return internal::cross_impl::run(*this, other); } namespace internal { diff --git a/doc/QuickReference.dox b/doc/QuickReference.dox index c5dfce421..e96b61767 100644 --- a/doc/QuickReference.dox +++ b/doc/QuickReference.dox @@ -367,7 +367,8 @@ vec2 = vec1.normalized(); vec1.normalize(); // inplace \endcode \link MatrixBase::cross() cross product \endlink \matrixworld\code #include -vec3 = vec1.cross(vec2);\endcode +v3c = v3a.cross(v3b); // size-3 vectors +scalar = v2a.cross(v2b); // size-2 vectors \endcode top diff --git a/doc/TutorialMatrixArithmetic.dox b/doc/TutorialMatrixArithmetic.dox index 5fc569a30..53916c241 100644 --- a/doc/TutorialMatrixArithmetic.dox +++ b/doc/TutorialMatrixArithmetic.dox @@ -158,7 +158,7 @@ For dot product and cross product, you need the \link MatrixBase::dot() dot()\en \verbinclude tut_arithmetic_dot_cross.out -Remember that cross product is only for vectors of size 3. Dot product is for vectors of any sizes. +Cross product is defined in Eigen not only for vectors of size 3 but also for those of size 2, check \link MatrixBase::cross() the doc\endlink for details. Dot product is for vectors of any sizes. When using complex numbers, Eigen's dot product is conjugate-linear in the first variable and linear in the second variable. diff --git a/doc/examples/tut_arithmetic_dot_cross.cpp b/doc/examples/tut_arithmetic_dot_cross.cpp index 5b0fd1eb2..d95e03c9d 100644 --- a/doc/examples/tut_arithmetic_dot_cross.cpp +++ b/doc/examples/tut_arithmetic_dot_cross.cpp @@ -9,5 +9,10 @@ int main() std::cout << "Dot product: " << v.dot(w) << std::endl; double dp = v.adjoint()*w; // automatic conversion of the inner product to a scalar std::cout << "Dot product via a matrix product: " << dp << std::endl; + std::cout << "Cross product:\n" << v.cross(w) << std::endl; + Eigen::Vector2d v2(1,2); + Eigen::Vector2d w2(0,1); + double cp = v2.cross(w2); // returning a scalar between size-2 vectors + std::cout << "Cross product for 2D vectors: " << cp << std::endl; } diff --git a/test/geo_orthomethods.cpp b/test/geo_orthomethods.cpp index 5f7ddb91f..64b392756 100644 --- a/test/geo_orthomethods.cpp +++ b/test/geo_orthomethods.cpp @@ -78,6 +78,36 @@ template void orthomethods_3() VERIFY_IS_APPROX(v2.cross(v1), rv1.cross(v1)); } +template void orthomethods_2() +{ + typedef typename NumTraits::Real RealScalar; + typedef Matrix Vector2; + typedef Matrix Vector3; + + Vector3 v30 = Vector3::Random(), + v31 = Vector3::Random(); + Vector2 v20 = v30.template head<2>(); + Vector2 v21 = v31.template head<2>(); + + VERIFY_IS_MUCH_SMALLER_THAN(v20.cross(v20), Scalar(1)); + VERIFY_IS_MUCH_SMALLER_THAN(v21.cross(v21), Scalar(1)); + VERIFY_IS_APPROX(v20.cross(v21), v30.cross(v31).z()); + + Vector2 v20Rot90(numext::conj(-v20.y()), numext::conj(v20.x())); + VERIFY_IS_APPROX(v20.cross( v20Rot90), v20.squaredNorm()); + VERIFY_IS_APPROX(v20.cross(-v20Rot90), -v20.squaredNorm()); + Vector2 v21Rot90(numext::conj(-v21.y()), numext::conj(v21.x())); + VERIFY_IS_APPROX(v21.cross( v21Rot90), v21.squaredNorm()); + VERIFY_IS_APPROX(v21.cross(-v21Rot90), -v21.squaredNorm()); + + // check mixed product + typedef Matrix RealVector2; + RealVector2 rv21 = RealVector2::Random(); + v21 = rv21.template cast(); + VERIFY_IS_APPROX(v20.cross(v21), v20.cross(rv21)); + VERIFY_IS_APPROX(v21.cross(v20), rv21.cross(v20)); +} + template void orthomethods(int size=Size) { typedef typename NumTraits::Real RealScalar; @@ -119,6 +149,9 @@ template void orthomethods(int size=Size) EIGEN_DECLARE_TEST(geo_orthomethods) { for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( orthomethods_2() ); + CALL_SUBTEST_2( orthomethods_2() ); + CALL_SUBTEST_4( orthomethods_2 >() ); CALL_SUBTEST_1( orthomethods_3() ); CALL_SUBTEST_2( orthomethods_3() ); CALL_SUBTEST_4( orthomethods_3 >() );