diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index f318bfd5d..096438f96 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -399,8 +399,16 @@ 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 internal::scalar_product_traits::Scalar,typename internal::traits::Scalar>::ReturnType Scalar; + typedef Matrix type; + }; + #endif // EIGEN_PARSED_BY_DOXYGEN template - PlainObject cross(const MatrixBase& other) const; + typename cross_product_return_type::type + cross(const MatrixBase& other) const; template PlainObject cross3(const MatrixBase& other) const; PlainObject unitOrthogonal(void) const; diff --git a/Eigen/src/Geometry/OrthoMethods.h b/Eigen/src/Geometry/OrthoMethods.h index 64ca52e16..52b469881 100644 --- a/Eigen/src/Geometry/OrthoMethods.h +++ b/Eigen/src/Geometry/OrthoMethods.h @@ -35,7 +35,7 @@ */ template template -inline typename MatrixBase::PlainObject +inline typename MatrixBase::template cross_product_return_type::type MatrixBase::cross(const MatrixBase& other) const { EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Derived,3) @@ -45,10 +45,10 @@ MatrixBase::cross(const MatrixBase& other) const // optimize such a small temporary very well (even within a complex expression) const typename internal::nested::type lhs(derived()); const typename internal::nested::type rhs(other.derived()); - return typename internal::plain_matrix_type::type( - lhs.coeff(1) * rhs.coeff(2) - lhs.coeff(2) * rhs.coeff(1), - lhs.coeff(2) * rhs.coeff(0) - lhs.coeff(0) * rhs.coeff(2), - lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0) + return typename cross_product_return_type::type( + internal::conj(lhs.coeff(1) * rhs.coeff(2) - lhs.coeff(2) * rhs.coeff(1)), + internal::conj(lhs.coeff(2) * rhs.coeff(0) - lhs.coeff(0) * rhs.coeff(2)), + internal::conj(lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0)) ); } @@ -62,9 +62,9 @@ struct cross3_impl { run(const VectorLhs& lhs, const VectorRhs& rhs) { return typename internal::plain_matrix_type::type( - lhs.coeff(1) * rhs.coeff(2) - lhs.coeff(2) * rhs.coeff(1), - lhs.coeff(2) * rhs.coeff(0) - lhs.coeff(0) * rhs.coeff(2), - lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0), + internal::conj(lhs.coeff(1) * rhs.coeff(2) - lhs.coeff(2) * rhs.coeff(1)), + internal::conj(lhs.coeff(2) * rhs.coeff(0) - lhs.coeff(0) * rhs.coeff(2)), + internal::conj(lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0)), 0 ); } @@ -121,16 +121,16 @@ VectorwiseOp::cross(const MatrixBase& ot if(Direction==Vertical) { eigen_assert(CrossReturnType::RowsAtCompileTime==3 && "the matrix must have exactly 3 rows"); - res.row(0) = _expression().row(1) * other.coeff(2) - _expression().row(2) * other.coeff(1); - res.row(1) = _expression().row(2) * other.coeff(0) - _expression().row(0) * other.coeff(2); - res.row(2) = _expression().row(0) * other.coeff(1) - _expression().row(1) * other.coeff(0); + res.row(0) = (_expression().row(1) * other.coeff(2) - _expression().row(2) * other.coeff(1)).conjugate(); + res.row(1) = (_expression().row(2) * other.coeff(0) - _expression().row(0) * other.coeff(2)).conjugate(); + res.row(2) = (_expression().row(0) * other.coeff(1) - _expression().row(1) * other.coeff(0)).conjugate(); } else { eigen_assert(CrossReturnType::ColsAtCompileTime==3 && "the matrix must have exactly 3 columns"); - res.col(0) = _expression().col(1) * other.coeff(2) - _expression().col(2) * other.coeff(1); - res.col(1) = _expression().col(2) * other.coeff(0) - _expression().col(0) * other.coeff(2); - res.col(2) = _expression().col(0) * other.coeff(1) - _expression().col(1) * other.coeff(0); + res.col(0) = (_expression().col(1) * other.coeff(2) - _expression().col(2) * other.coeff(1)).conjugate(); + res.col(1) = (_expression().col(2) * other.coeff(0) - _expression().col(0) * other.coeff(2)).conjugate(); + res.col(2) = (_expression().col(0) * other.coeff(1) - _expression().col(1) * other.coeff(0)).conjugate(); } return res; } diff --git a/test/geo_orthomethods.cpp b/test/geo_orthomethods.cpp index 9152d99c8..020ae7103 100644 --- a/test/geo_orthomethods.cpp +++ b/test/geo_orthomethods.cpp @@ -33,6 +33,7 @@ template void orthomethods_3() { + typedef typename NumTraits::Real RealScalar; typedef Matrix Matrix3; typedef Matrix Vector3; @@ -44,6 +45,9 @@ template void orthomethods_3() // cross product VERIFY_IS_MUCH_SMALLER_THAN(v1.cross(v2).dot(v1), Scalar(1)); + VERIFY_IS_MUCH_SMALLER_THAN(v1.dot(v1.cross(v2)), Scalar(1)); + VERIFY_IS_MUCH_SMALLER_THAN(v1.cross(v2).dot(v2), Scalar(1)); + VERIFY_IS_MUCH_SMALLER_THAN(v2.dot(v1.cross(v2)), Scalar(1)); Matrix3 mat3; mat3 << v0.normalized(), (v0.cross(v1)).normalized(), @@ -68,6 +72,12 @@ template void orthomethods_3() v40.w() = v41.w() = v42.w() = 0; v42.template head<3>() = v40.template head<3>().cross(v41.template head<3>()); VERIFY_IS_APPROX(v40.cross3(v41), v42); + + // check mixed product + typedef Matrix RealVector3; + RealVector3 rv1 = RealVector3::Random(); + VERIFY_IS_APPROX(v1.cross(rv1.template cast()), v1.cross(rv1)); + VERIFY_IS_APPROX(rv1.template cast().cross(v1), rv1.cross(v1)); } template void orthomethods(int size=Size) @@ -115,6 +125,7 @@ void test_geo_orthomethods() for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1( orthomethods_3() ); CALL_SUBTEST_2( orthomethods_3() ); + CALL_SUBTEST_4( orthomethods_3 >() ); CALL_SUBTEST_1( (orthomethods()) ); CALL_SUBTEST_2( (orthomethods()) ); CALL_SUBTEST_1( (orthomethods()) );