GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-04 18:54:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

fix cross product for complexes and add support for mixed real-complex cross products

Browse Source
This commit is contained in:
Gael Guennebaud 2011-01-27 11:33:37 +01:00
parent 0bfb78c824
commit e8d6a5ca87
3 changed files with 34 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
Eigen/src/Core/MatrixBase.h
View File

@ -399,8 +399,16 @@ template<typename Derived> class MatrixBase
/////////// Geometry module ///////////
#ifndef EIGEN_PARSED_BY_DOXYGEN
/// \internal helper struct to form the return type of the cross product
template<typename OtherDerived> struct cross_product_return_type {
typedef typename internal::scalar_product_traits<typename internal::traits<Derived>::Scalar,typename internal::traits<OtherDerived>::Scalar>::ReturnType Scalar;
typedef Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime> type;
};
#endif // EIGEN_PARSED_BY_DOXYGEN
template<typename OtherDerived>
PlainObject cross(const MatrixBase<OtherDerived>& other) const;
typename cross_product_return_type<OtherDerived>::type
cross(const MatrixBase<OtherDerived>& other) const;
template<typename OtherDerived>
PlainObject cross3(const MatrixBase<OtherDerived>& other) const;
PlainObject unitOrthogonal(void) const;

28
Eigen/src/Geometry/OrthoMethods.h
View File

@ -35,7 +35,7 @@
*/
template<typename Derived>
template<typename OtherDerived>
inline typename MatrixBase<Derived>::PlainObject
inline typename MatrixBase<Derived>::template cross_product_return_type<OtherDerived>::type
MatrixBase<Derived>::cross(const MatrixBase<OtherDerived>& other) const
{
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Derived,3)
@ -45,10 +45,10 @@ MatrixBase<Derived>::cross(const MatrixBase<OtherDerived>& other) const
// optimize such a small temporary very well (even within a complex expression)
const typename internal::nested<Derived,2>::type lhs(derived());
const typename internal::nested<OtherDerived,2>::type rhs(other.derived());
return typename internal::plain_matrix_type<Derived>::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<OtherDerived>::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<VectorLhs>::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<ExpressionType,Direction>::cross(const MatrixBase<OtherDerived>& 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;
}

11
test/geo_orthomethods.cpp
View File

@ -33,6 +33,7 @@
template<typename Scalar> void orthomethods_3()
{
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar,3,3> Matrix3;
typedef Matrix<Scalar,3,1> Vector3;
@ -44,6 +45,9 @@ template<typename Scalar> 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<typename Scalar> 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<RealScalar, 3, 1> RealVector3;
RealVector3 rv1 = RealVector3::Random();
VERIFY_IS_APPROX(v1.cross(rv1.template cast<Scalar>()), v1.cross(rv1));
VERIFY_IS_APPROX(rv1.template cast<Scalar>().cross(v1), rv1.cross(v1));
}
template<typename Scalar, int Size> 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<float>() );
CALL_SUBTEST_2( orthomethods_3<double>() );
CALL_SUBTEST_4( orthomethods_3<std::complex<double> >() );
CALL_SUBTEST_1( (orthomethods<float,2>()) );
CALL_SUBTEST_2( (orthomethods<double,2>()) );
CALL_SUBTEST_1( (orthomethods<float,3>()) );