* fix Quaternion::setFromTwoVectors (thanks to "benv" from the forum)

* extend PartialRedux::cross() to any matrix sizes with automatic
  vectorization when possible
* unit tests: add "geo_" prefix to all unit tests related to the
  geometry module and start splitting the big "geometry.cpp" tests to
  multiple smaller ones (also include new tests)
This commit is contained in:
Gael Guennebaud 2009-02-17 09:53:05 +00:00
parent 67b4fab4e3
commit e6f1104b57
16 changed files with 445 additions and 157 deletions

View File

@ -1,7 +1,7 @@
// 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) 2008-2009 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
@ -176,7 +176,7 @@ template<typename ExpressionType, int Direction> class PartialRedux
};
typedef typename ExpressionType::PlainMatrixType CrossReturnType;
inline PartialRedux(const ExpressionType& matrix) : m_matrix(matrix) {}
/** \internal */
@ -249,7 +249,7 @@ template<typename ExpressionType, int Direction> class PartialRedux
* \sa MatrixBase::any() */
const typename ReturnType<ei_member_any>::Type any() const
{ return _expression(); }
/** \returns a row (or column) vector expression representing
* the number of \c true coefficients of each respective column (or row).
*
@ -269,8 +269,8 @@ template<typename ExpressionType, int Direction> class PartialRedux
* \sa MatrixBase::prod() */
const typename ReturnType<ei_member_prod>::Type prod() const
{ return _expression(); }
/** \returns a matrix expression
* where each column (or row) are reversed.
*
@ -282,33 +282,9 @@ template<typename ExpressionType, int Direction> class PartialRedux
{
return Reverse<ExpressionType, Direction>( _expression() );
}
/** \returns a 3x3 matrix expression of the cross product
* of each column or row of the referenced expression with the \a other vector.
*
* \geometry_module
*
* \sa MatrixBase::cross() */
template<typename OtherDerived>
const CrossReturnType cross(const MatrixBase<OtherDerived>& other) const
{
EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(CrossReturnType,3,3)
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived,3)
EIGEN_STATIC_ASSERT((ei_is_same_type<Scalar, typename OtherDerived::Scalar>::ret),
YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
if(Direction==Vertical)
return (CrossReturnType()
<< _expression().col(0).cross(other),
_expression().col(1).cross(other),
_expression().col(2).cross(other)).finished();
else
return (CrossReturnType()
<< _expression().row(0).cross(other),
_expression().row(1).cross(other),
_expression().row(2).cross(other)).finished();
}
template<typename OtherDerived>
const CrossReturnType cross(const MatrixBase<OtherDerived>& other) const;
protected:
ExpressionTypeNested m_matrix;

View File

@ -1,7 +1,7 @@
// 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) 2008-2009 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
@ -50,6 +50,42 @@ MatrixBase<Derived>::cross(const MatrixBase<OtherDerived>& other) const
);
}
/** \returns a matrix expression of the cross product of each column or row
* of the referenced expression with the \a other vector.
*
* The referenced matrix must have one dimension equal to 3.
* The result matrix has the same dimensions than the referenced one.
*
* \geometry_module
*
* \sa MatrixBase::cross() */
template<typename ExpressionType, int Direction>
template<typename OtherDerived>
const typename PartialRedux<ExpressionType,Direction>::CrossReturnType
PartialRedux<ExpressionType,Direction>::cross(const MatrixBase<OtherDerived>& other) const
{
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived,3)
EIGEN_STATIC_ASSERT((ei_is_same_type<Scalar, typename OtherDerived::Scalar>::ret),
YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
CrossReturnType res(_expression().rows(),_expression().cols());
if(Direction==Vertical)
{
ei_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);
}
else
{
ei_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);
}
return res;
}
template<typename Derived, int Size = Derived::SizeAtCompileTime>
struct ei_unitOrthogonal_selector
{

View File

@ -345,7 +345,6 @@ inline Quaternion<Scalar>& Quaternion<Scalar>::setFromTwoVectors(const MatrixBas
{
Vector3 v0 = a.normalized();
Vector3 v1 = b.normalized();
Vector3 axis = v0.cross(v1);
Scalar c = v0.dot(v1);
// if dot == 1, vectors are the same
@ -353,7 +352,17 @@ inline Quaternion<Scalar>& Quaternion<Scalar>::setFromTwoVectors(const MatrixBas
{
// set to identity
this->w() = 1; this->vec().setZero();
return *this;
}
// if dot == -1, vectors are opposites
if (ei_isApprox(c,Scalar(-1)))
{
this->vec() = v0.unitOrthogonal();
this->w() = 0;
return *this;
}
Vector3 axis = v0.cross(v1);
Scalar s = ei_sqrt((Scalar(1)+c)*Scalar(2));
Scalar invs = Scalar(1)/s;
this->vec() = axis * invs;

View File

@ -55,7 +55,7 @@ class FancySpheres
mCenters.push_back(Vector3f::Zero());
parents.push_back(-1);
mRadii.push_back(radius);
// generate level 1 using icosphere vertices
radius *= 0.45;
{
@ -89,7 +89,7 @@ class FancySpheres
Vector3f ax1 = ax0.unitOrthogonal();
Quaternionf q;
q.setFromTwoVectors(Vector3f::UnitZ(), ax0);
Transform3f t = Translation3f(c) * q * Scaling3f(mRadii[i]+radius);
Transform3f t = Translation3f(c) * q * Scaling(mRadii[i]+radius);
for (int j=0; j<5; ++j)
{
Vector3f newC = c + ( (AngleAxisf(angles[j*2+1], ax0)
@ -103,14 +103,14 @@ class FancySpheres
start = end;
}
}
void draw()
{
int end = mCenters.size();
glEnable(GL_NORMALIZE);
for (int i=0; i<end; ++i)
{
Transform3f t = Translation3f(mCenters[i]) * Scaling3f(mRadii[i]);
Transform3f t = Translation3f(mCenters[i]) * Scaling(mRadii[i]);
gpu.pushMatrix(GL_MODELVIEW);
gpu.multMatrix(t.matrix(),GL_MODELVIEW);
mIcoSphere.draw(2);
@ -135,7 +135,7 @@ template<typename T> T lerp(float t, const T& a, const T& b)
template<> Quaternionf lerp(float t, const Quaternionf& a, const Quaternionf& b)
{ return a.slerp(t,b); }
// linear interpolation of a frame using the type OrientationType
// linear interpolation of a frame using the type OrientationType
// to perform the interpolation of the orientations
template<typename OrientationType>
inline static Frame lerpFrame(float alpha, const Frame& a, const Frame& b)
@ -171,7 +171,7 @@ public:
Matrix3 m = q.toRotationMatrix();
return *this = m;
}
EulerAngles& operator=(const Matrix3& m)
{
// mat = cy*cz -cy*sz sy
@ -240,7 +240,7 @@ void RenderingWidget::grabFrame(void)
void RenderingWidget::drawScene()
{
static FancySpheres sFancySpheres;
float length = 200;
float length = 50;
gpu.drawVector(Vector3f::Zero(), length*Vector3f::UnitX(), Color(1,0,0,1));
gpu.drawVector(Vector3f::Zero(), length*Vector3f::UnitY(), Color(0,1,0,1));
gpu.drawVector(Vector3f::Zero(), length*Vector3f::UnitZ(), Color(0,0,1,1));
@ -262,12 +262,60 @@ void RenderingWidget::drawScene()
glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, Vector4f(1, 1, 1, 1).data());
glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 64);
glEnable(GL_LIGHTING);
glEnable(GL_LIGHT0);
glEnable(GL_LIGHT1);
// glEnable(GL_LIGHTING);
// glEnable(GL_LIGHT0);
// glEnable(GL_LIGHT1);
glColor3f(0.4, 0.7, 0.4);
sFancySpheres.draw();
{IcoSphere s(5);
float length = 6;
// const std::vector<int>& indices = s.indices(2);
for (unsigned int i=0; i<s.vertices().size() ; ++i)
{
Vector3f n = s.vertices()[i].normalized();
Vector3f p = n*100;
int i,j;
float minc = n.minCoeff(&i);
float maxc = n.maxCoeff(&j);
Vector3f o = n;
o[i] = -n[j];
o[j] = n[i];
//Vector3f u = n.unitOrthogonal().normalized();
Vector3f u = n.cross(o).normalized();
Vector3f v = n.cross(u).normalized();
gpu.drawVector(p, length*u, Color(1,0,0,1));
gpu.drawVector(p, length*v, Color(0,1,0,1));
// gpu.drawVector(p, length*n, Color(0,0,1,1));
}}
for(;;)
{
Vector3f n = Vector3f::Random().normalized();
// Vector3f p = n*100;
//Vector3f u = n.unitOrthogonal().normalized();
// int i,j;
// float minc = n.minCoeff(&i);
// float maxc = n.maxCoeff(&j);
// Vector3f o = n;
// o[i] = -n[j];
// o[j] = n[i];
Vector3f o = Vector3f(-n.y(),n.z(),n.x()).normalized();
Vector3f u = n.cross(o).normalized();
float d = ei_abs(o.dot(n));
if (d>0.9999)
std::cout << d << " " << n.transpose() << "\n";
// Vector3f v = n.cross(u).normalized();
}
// sFancySpheres.draw();
// glVertexPointer(3, GL_FLOAT, 0, mVertices[0].data());
// glNormalPointer(GL_FLOAT, 0, mNormals[0].data());
// glEnableClientState(GL_VERTEX_ARRAY);
@ -275,7 +323,7 @@ void RenderingWidget::drawScene()
// glDrawArrays(GL_TRIANGLES, 0, mVertices.size());
// glDisableClientState(GL_VERTEX_ARRAY);
// glDisableClientState(GL_NORMAL_ARRAY);
glDisable(GL_LIGHTING);
}
@ -318,7 +366,7 @@ void RenderingWidget::animate()
currentFrame.orientation = currentFrame.orientation.inverse();
currentFrame.position = - (currentFrame.orientation * currentFrame.position);
mCamera.setFrame(currentFrame);
updateGL();
}
@ -542,7 +590,7 @@ void RenderingWidget::resetCamera()
// put the camera at that time step:
aux1 = aux0.lerp(duration/2,mInitFrame);
// and make it look at teh target again
// and make it look at the target again
aux1.orientation = aux1.orientation.inverse();
aux1.position = - (aux1.orientation * aux1.position);
mCamera.setFrame(aux1);
@ -647,6 +695,17 @@ QuaternionDemo::QuaternionDemo()
int main(int argc, char *argv[])
{
std::cout << "Navigation:\n";
std::cout << " left button: rotate around the target\n";
std::cout << " middle button: zoom\n";
std::cout << " left button + ctrl quake rotate (rotate around camera position)\n";
std::cout << " middle button + ctrl walk (progress along camera's z direction)\n";
std::cout << " left button: pan (translate in the XY camera's plane)\n\n";
std::cout << "R : move the camera to initial position\n";
std::cout << "A : start/stop animation\n";
std::cout << "C : clear the animation\n";
std::cout << "G : add a key frame\n";
QApplication app(argc, argv);
QuaternionDemo demo;
demo.resize(600,500);

View File

@ -14,7 +14,7 @@ Using STL containers on \ref FixedSizeVectorizable "fixed-size vectorizable Eige
\li A 16-byte-aligned allocator must be used. Eigen does provide one ready for use: aligned_allocator.
\li If you want to use the std::vector container, you need to \#include <Eigen/StdVector>.
These issues arise only with \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types". For other Eigen types, such as Vector3f or MatrixXd, no special care is needed when using STL containers.
These issues arise only with \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types" and \ref StructHavingEigenMembers "structures having such Eigen objects as member". For other Eigen types, such as Vector3f or MatrixXd, no special care is needed when using STL containers.
\section allocator Using an aligned allocator
@ -32,7 +32,7 @@ Note that here, the 3rd parameter "std::less<int>" is just the default value, we
\section vector The case of std::vector
The situation with std::vector was even worse (explanation below) so we had to specialize it for Eigen types. The upside is that our specialization takes care of specifying the aligned allocator, so you don't need to worry about it. All you need to do is to \#include <Eigen/StdVector>.
The situation with std::vector was even worse (explanation below) so we had to specialize it for Eigen types. The upside is that our specialization takes care of specifying the aligned allocator, so you don't need to worry about it. All you need to do is to \#include <Eigen/StdVector> instead of (or before) \#include <vector>.
So as soon as you have
\code

View File

@ -50,7 +50,7 @@ std::map<int, Eigen::Matrix2f> my_map;
then you need to read this separate page: \ref StlContainers "Using STL Containers with Eigen".
Note that here, Eigen::Matrix2f is only used as an example, more generally the issue arises for all \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types".
Note that here, Eigen::Matrix2f is only used as an example, more generally the issue arises for all \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types" and \ref StructHavingEigenMembers "structures having such Eigen objects as member".
\section c3 Cause 3: Passing Eigen objects by value

View File

@ -31,7 +31,7 @@ Fixed-size objects are typically very small, which means that we want to handle
Now, vectorization (both SSE and AltiVec) works with 128-bit packets. Moreover, for performance reasons, these packets need to be have 128-bit alignment.
So it turns out that the only way that fixed-size Eigen objects can be vectorized, is if their size is a multiple of 128 bits, or 16 bytes. Eigen will then request 16-byte alignment for these object, and henceforth rely on these objects being aligned so no runtime check for alignment is performed.
So it turns out that the only way that fixed-size Eigen objects can be vectorized, is if their size is a multiple of 128 bits, or 16 bytes. Eigen will then request 16-byte alignment for these objects, and henceforth rely on these objects being aligned so no runtime check for alignment is performed.
*/

View File

@ -129,10 +129,13 @@ ei_add_test(inverse)
ei_add_test(qr)
ei_add_test(eigensolver " " "${GSL_LIBRARIES}")
ei_add_test(svd)
ei_add_test(geometry)
ei_add_test(hyperplane)
ei_add_test(parametrizedline)
ei_add_test(alignedbox)
ei_add_test(geo_orthomethods)
ei_add_test(geo_quaternion)
ei_add_test(geo_transformations)
ei_add_test(geo_eulerangles)
ei_add_test(geo_hyperplane)
ei_add_test(geo_parametrizedline)
ei_add_test(geo_alignedbox)
ei_add_test(regression)
ei_add_test(stdvector)
ei_add_test(resize)

View File

@ -1,7 +1,7 @@
// 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) 2008-2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@ -65,7 +65,7 @@ template<typename BoxType> void alignedbox(const BoxType& _box)
VERIFY_IS_APPROX(hp1d.template cast<Scalar>(),b0);
}
void test_alignedbox()
void test_geo_alignedbox()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( alignedbox(AlignedBox<float,2>()) );

70
test/geo_eulerangles.cpp Normal file
View File

@ -0,0 +1,70 @@
// 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-2009 Gael Guennebaud <g.gael@free.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/>.
#include "main.h"
#include <Eigen/Geometry>
#include <Eigen/LU>
#include <Eigen/SVD>
template<typename Scalar> void eulerangles(void)
{
typedef Matrix<Scalar,3,3> Matrix3;
typedef Matrix<Scalar,3,1> Vector3;
typedef Quaternion<Scalar> Quaternionx;
typedef AngleAxis<Scalar> AngleAxisx;
Scalar a = ei_random<Scalar>(-Scalar(M_PI), Scalar(M_PI));
Quaternionx q1;
q1 = AngleAxisx(a, Vector3::Random().normalized());
Matrix3 m;
m = q1;
#define VERIFY_EULER(I,J,K, X,Y,Z) { \
Vector3 ea = m.eulerAngles(I,J,K); \
Matrix3 m1 = Matrix3(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z())); \
VERIFY_IS_APPROX(m, Matrix3(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z()))); \
}
VERIFY_EULER(0,1,2, X,Y,Z);
VERIFY_EULER(0,1,0, X,Y,X);
VERIFY_EULER(0,2,1, X,Z,Y);
VERIFY_EULER(0,2,0, X,Z,X);
VERIFY_EULER(1,2,0, Y,Z,X);
VERIFY_EULER(1,2,1, Y,Z,Y);
VERIFY_EULER(1,0,2, Y,X,Z);
VERIFY_EULER(1,0,1, Y,X,Y);
VERIFY_EULER(2,0,1, Z,X,Y);
VERIFY_EULER(2,0,2, Z,X,Z);
VERIFY_EULER(2,1,0, Z,Y,X);
VERIFY_EULER(2,1,2, Z,Y,Z);
}
void test_geo_eulerangles()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( eulerangles<float>() );
CALL_SUBTEST( eulerangles<double>() );
}
}

View File

@ -127,7 +127,7 @@ template<typename Scalar> void lines()
}
}
void test_hyperplane()
void test_geo_hyperplane()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( hyperplane(Hyperplane<float,2>()) );

109
test/geo_orthomethods.cpp Normal file
View File

@ -0,0 +1,109 @@
// 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-2009 Gael Guennebaud <g.gael@free.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/>.
#include "main.h"
#include <Eigen/Geometry>
#include <Eigen/LU>
#include <Eigen/SVD>
/* this test covers the following files:
Geometry/OrthoMethods.h
*/
template<typename Scalar> void orthomethods_3()
{
typedef Matrix<Scalar,3,3> Matrix3;
typedef Matrix<Scalar,3,1> Vector3;
Vector3 v0 = Vector3::Random(),
v1 = Vector3::Random(),
v2 = Vector3::Random();
// cross product
VERIFY_IS_MUCH_SMALLER_THAN(v1.cross(v2).dot(v1), Scalar(1));
Matrix3 mat3;
mat3 << v0.normalized(),
(v0.cross(v1)).normalized(),
(v0.cross(v1).cross(v0)).normalized();
VERIFY(mat3.isUnitary());
// colwise/rowwise cross product
mat3.setRandom();
Vector3 vec3 = Vector3::Random();
Matrix3 mcross;
int i = ei_random<int>(0,2);
mcross = mat3.colwise().cross(vec3);
VERIFY_IS_APPROX(mcross.col(i), mat3.col(i).cross(vec3));
mcross = mat3.rowwise().cross(vec3);
VERIFY_IS_APPROX(mcross.row(i), mat3.row(i).cross(vec3));
}
template<typename Scalar, int Size> void orthomethods(int size=Size)
{
typedef Matrix<Scalar,Size,1> VectorType;
typedef Matrix<Scalar,3,Size> Matrix3N;
typedef Matrix<Scalar,Size,3> MatrixN3;
typedef Matrix<Scalar,3,1> Vector3;
VectorType v0 = VectorType::Random(size),
v1 = VectorType::Random(size),
v2 = VectorType::Random(size);
// unitOrthogonal
VERIFY_IS_MUCH_SMALLER_THAN(v0.unitOrthogonal().dot(v0), Scalar(1));
VERIFY_IS_APPROX(v0.unitOrthogonal().norm(), Scalar(1));
// colwise/rowwise cross product
Vector3 vec3 = Vector3::Random();
int i = ei_random<int>(0,size-1);
Matrix3N mat3N(3,size), mcross3N(3,size);
mat3N.setRandom();
mcross3N = mat3N.colwise().cross(vec3);
VERIFY_IS_APPROX(mcross3N.col(i), mat3N.col(i).cross(vec3));
MatrixN3 matN3(size,3), mcrossN3(size,3);
matN3.setRandom();
mcrossN3 = matN3.rowwise().cross(vec3);
VERIFY_IS_APPROX(mcrossN3.row(i), matN3.row(i).cross(vec3));
}
void test_geo_orthomethods()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( orthomethods_3<float>() );
CALL_SUBTEST( orthomethods_3<double>() );
CALL_SUBTEST( (orthomethods<float,2>()) );
CALL_SUBTEST( (orthomethods<double,2>()) );
CALL_SUBTEST( (orthomethods<float,3>()) );
CALL_SUBTEST( (orthomethods<double,3>()) );
CALL_SUBTEST( (orthomethods<float,7>()) );
CALL_SUBTEST( (orthomethods<double,8>()) );
CALL_SUBTEST( (orthomethods<float,Dynamic>(36)) );
CALL_SUBTEST( (orthomethods<double,Dynamic>(35)) );
}
}

View File

@ -66,7 +66,7 @@ template<typename LineType> void parametrizedline(const LineType& _line)
VERIFY_IS_APPROX(hp1d.template cast<Scalar>(),l0);
}
void test_parametrizedline()
void test_geo_parametrizedline()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( parametrizedline(ParametrizedLine<float,2>()) );

116
test/geo_quaternion.cpp Normal file
View File

@ -0,0 +1,116 @@
// 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-2009 Gael Guennebaud <g.gael@free.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/>.
#include "main.h"
#include <Eigen/Geometry>
#include <Eigen/LU>
#include <Eigen/SVD>
template<typename Scalar> void quaternion(void)
{
/* this test covers the following files:
Quaternion.h
*/
typedef Matrix<Scalar,3,3> Matrix3;
typedef Matrix<Scalar,3,1> Vector3;
typedef Quaternion<Scalar> Quaternionx;
typedef AngleAxis<Scalar> AngleAxisx;
Scalar largeEps = test_precision<Scalar>();
if (ei_is_same_type<Scalar,float>::ret)
largeEps = 1e-3f;
Scalar eps = ei_random<Scalar>() * 1e-2;
Vector3 v0 = Vector3::Random(),
v1 = Vector3::Random(),
v2 = Vector3::Random(),
v3 = Vector3::Random();
Scalar a = ei_random<Scalar>(-Scalar(M_PI), Scalar(M_PI));
// Quaternion: Identity(), setIdentity();
Quaternionx q1, q2;
q2.setIdentity();
VERIFY_IS_APPROX(Quaternionx(Quaternionx::Identity()).coeffs(), q2.coeffs());
q1.coeffs().setRandom();
VERIFY_IS_APPROX(q1.coeffs(), (q1*q2).coeffs());
q1 = AngleAxisx(a, v0.normalized());
q2 = AngleAxisx(a, v1.normalized());
// angular distance
Scalar refangle = ei_abs(AngleAxisx(q1.inverse()*q2).angle());
if (refangle>Scalar(M_PI))
refangle = Scalar(2)*Scalar(M_PI) - refangle;
if((q1.coeffs()-q2.coeffs()).norm() > 10*largeEps)
{
VERIFY(ei_isApprox(q1.angularDistance(q2), refangle, largeEps));
}
// rotation matrix conversion
VERIFY_IS_APPROX(q1 * v2, q1.toRotationMatrix() * v2);
VERIFY_IS_APPROX(q1 * q2 * v2,
q1.toRotationMatrix() * q2.toRotationMatrix() * v2);
VERIFY( (q2*q1).isApprox(q1*q2, largeEps)
|| !(q2 * q1 * v2).isApprox(q1.toRotationMatrix() * q2.toRotationMatrix() * v2));
q2 = q1.toRotationMatrix();
VERIFY_IS_APPROX(q1*v1,q2*v1);
// angle-axis conversion
AngleAxisx aa = q1;
VERIFY_IS_APPROX(q1 * v1, Quaternionx(aa) * v1);
VERIFY_IS_NOT_APPROX(q1 * v1, Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1);
// from two vector creation
VERIFY_IS_APPROX( v2.normalized(),(q2.setFromTwoVectors(v1, v2)*v1).normalized());
VERIFY_IS_APPROX( v1.normalized(),(q2.setFromTwoVectors(v1, v1)*v1).normalized());
VERIFY_IS_APPROX(-v1.normalized(),(q2.setFromTwoVectors(v1,-v1)*v1).normalized());
v3 = v1.cwise()+eps;
VERIFY_IS_APPROX( v3.normalized(),(q2.setFromTwoVectors(v1, v3)*v1).normalized());
VERIFY_IS_APPROX(-v3.normalized(),(q2.setFromTwoVectors(v1,-v3)*v1).normalized());
// inverse and conjugate
VERIFY_IS_APPROX(q1 * (q1.inverse() * v1), v1);
VERIFY_IS_APPROX(q1 * (q1.conjugate() * v1), v1);
// test casting
Quaternion<float> q1f = q1.template cast<float>();
VERIFY_IS_APPROX(q1f.template cast<Scalar>(),q1);
Quaternion<double> q1d = q1.template cast<double>();
VERIFY_IS_APPROX(q1d.template cast<Scalar>(),q1);
}
void test_geo_quaternion()
{
for(int i = 0; i < g_repeat; i++) {
// CALL_SUBTEST( quaternion<float>() );
CALL_SUBTEST( quaternion<double>() );
}
}

View File

@ -1,7 +1,7 @@
// 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) 2008-2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@ -27,7 +27,7 @@
#include <Eigen/LU>
#include <Eigen/SVD>
template<typename Scalar> void geometry(void)
template<typename Scalar> void transformations(void)
{
/* this test covers the following files:
Cross.h Quaternion.h, Transform.cpp
@ -56,32 +56,10 @@ template<typename Scalar> void geometry(void)
v1 = Vector3::Random(),
v2 = Vector3::Random();
Vector2 u0 = Vector2::Random();
Matrix3 matrot1;
Matrix3 matrot1, m;
Scalar a = ei_random<Scalar>(-Scalar(M_PI), Scalar(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());
// Quaternion: Identity(), setIdentity();
Quaternionx q1, q2;
q2.setIdentity();
VERIFY_IS_APPROX(Quaternionx(Quaternionx::Identity()).coeffs(), q2.coeffs());
q1.coeffs().setRandom();
VERIFY_IS_APPROX(q1.coeffs(), (q1*q2).coeffs());
// unitOrthogonal
VERIFY_IS_MUCH_SMALLER_THAN(u0.unitOrthogonal().dot(u0), Scalar(1));
VERIFY_IS_MUCH_SMALLER_THAN(v0.unitOrthogonal().dot(v0), Scalar(1));
VERIFY_IS_APPROX(u0.unitOrthogonal().norm(), Scalar(1));
VERIFY_IS_APPROX(v0.unitOrthogonal().norm(), Scalar(1));
VERIFY_IS_APPROX(v0, AngleAxisx(a, v0.normalized()) * v0);
VERIFY_IS_APPROX(-v0, AngleAxisx(Scalar(M_PI), v0.unitOrthogonal()) * v0);
VERIFY_IS_APPROX(ei_cos(a)*v0.squaredNorm(), v0.dot(AngleAxisx(a, v0.unitOrthogonal()) * v0));
@ -89,30 +67,11 @@ template<typename Scalar> void geometry(void)
VERIFY_IS_APPROX(Matrix3::Identity(), m * AngleAxisx(a, v0.normalized()));
VERIFY_IS_APPROX(Matrix3::Identity(), AngleAxisx(a, v0.normalized()) * m);
Quaternionx q1, q2;
q1 = AngleAxisx(a, v0.normalized());
q2 = AngleAxisx(a, v1.normalized());
// angular distance
Scalar refangle = ei_abs(AngleAxisx(q1.inverse()*q2).angle());
if (refangle>Scalar(M_PI))
refangle = Scalar(2)*Scalar(M_PI) - refangle;
if((q1.coeffs()-q2.coeffs()).norm() > 10*largeEps)
{
VERIFY(ei_isApprox(q1.angularDistance(q2), refangle, largeEps));
}
// rotation matrix conversion
VERIFY_IS_APPROX(q1 * v2, q1.toRotationMatrix() * v2);
VERIFY_IS_APPROX(q1 * q2 * v2,
q1.toRotationMatrix() * q2.toRotationMatrix() * v2);
VERIFY( (q2*q1).isApprox(q1*q2, largeEps) || !(q2 * q1 * v2).isApprox(
q1.toRotationMatrix() * q2.toRotationMatrix() * v2));
q2 = q1.toRotationMatrix();
VERIFY_IS_APPROX(q1*v1,q2*v1);
matrot1 = AngleAxisx(Scalar(0.1), Vector3::UnitX())
* AngleAxisx(Scalar(0.2), Vector3::UnitY())
* AngleAxisx(Scalar(0.3), Vector3::UnitZ());
@ -126,14 +85,6 @@ template<typename Scalar> void geometry(void)
VERIFY_IS_APPROX(q1 * v1, Quaternionx(aa) * v1);
VERIFY_IS_NOT_APPROX(q1 * v1, Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1);
// from two vector creation
VERIFY_IS_APPROX(v2.normalized(),(q2.setFromTwoVectors(v1,v2)*v1).normalized());
VERIFY_IS_APPROX(v2.normalized(),(q2.setFromTwoVectors(v1,v2)*v1).normalized());
// inverse and conjugate
VERIFY_IS_APPROX(q1 * (q1.inverse() * v1), v1);
VERIFY_IS_APPROX(q1 * (q1.conjugate() * v1), v1);
// AngleAxis
VERIFY_IS_APPROX(AngleAxisx(a,v1.normalized()).toRotationMatrix(),
Quaternionx(AngleAxisx(a,v1.normalized())).toRotationMatrix());
@ -382,11 +333,6 @@ template<typename Scalar> void geometry(void)
DiagonalMatrix<double,3> sc1d; sc1d = (sc1.template cast<double>());
VERIFY_IS_APPROX(sc1d.template cast<Scalar>(),sc1);
Quaternion<float> q1f = q1.template cast<float>();
VERIFY_IS_APPROX(q1f.template cast<Scalar>(),q1);
Quaternion<double> q1d = q1.template cast<double>();
VERIFY_IS_APPROX(q1d.template cast<Scalar>(),q1);
AngleAxis<float> aa1f = aa1.template cast<float>();
VERIFY_IS_APPROX(aa1f.template cast<Scalar>(),aa1);
AngleAxis<double> aa1d = aa1.template cast<double>();
@ -397,47 +343,12 @@ template<typename Scalar> void geometry(void)
VERIFY_IS_APPROX(r2d1f.template cast<Scalar>(),r2d1);
Rotation2D<double> r2d1d = r2d1.template cast<double>();
VERIFY_IS_APPROX(r2d1d.template cast<Scalar>(),r2d1);
m = q1;
// m.col(1) = Vector3(0,ei_random<Scalar>(),ei_random<Scalar>()).normalized();
// m.col(0) = Vector3(-1,0,0).normalized();
// m.col(2) = m.col(0).cross(m.col(1));
#define VERIFY_EULER(I,J,K, X,Y,Z) { \
Vector3 ea = m.eulerAngles(I,J,K); \
Matrix3 m1 = Matrix3(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z())); \
VERIFY_IS_APPROX(m, Matrix3(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z()))); \
}
VERIFY_EULER(0,1,2, X,Y,Z);
VERIFY_EULER(0,1,0, X,Y,X);
VERIFY_EULER(0,2,1, X,Z,Y);
VERIFY_EULER(0,2,0, X,Z,X);
VERIFY_EULER(1,2,0, Y,Z,X);
VERIFY_EULER(1,2,1, Y,Z,Y);
VERIFY_EULER(1,0,2, Y,X,Z);
VERIFY_EULER(1,0,1, Y,X,Y);
VERIFY_EULER(2,0,1, Z,X,Y);
VERIFY_EULER(2,0,2, Z,X,Z);
VERIFY_EULER(2,1,0, Z,Y,X);
VERIFY_EULER(2,1,2, Z,Y,Z);
// colwise/rowwise cross product
mat3.setRandom();
Vector3 vec3 = Vector3::Random();
Matrix3 mcross;
int i = ei_random<int>(0,2);
mcross = mat3.colwise().cross(vec3);
VERIFY_IS_APPROX(mcross.col(i), mat3.col(i).cross(vec3));
mcross = mat3.rowwise().cross(vec3);
VERIFY_IS_APPROX(mcross.row(i), mat3.row(i).cross(vec3));
}
void test_geometry()
void test_geo_transformations()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( geometry<float>() );
CALL_SUBTEST( geometry<double>() );
CALL_SUBTEST( transformations<float>() );
CALL_SUBTEST( transformations<double>() );
}
}

View File

@ -34,7 +34,6 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType&
double density = std::max(8./(rows*cols), 0.01);
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
typedef Matrix<Scalar,Dynamic,1> DenseVector;
Scalar eps = 1e-6;
// test matrix-matrix product
{