Commit 31acf5e6 authored by Sebastian Eichelbaum's avatar Sebastian Eichelbaum

[MERGE]

parents aeb66f95 23180c49
......@@ -24,8 +24,6 @@
#include <vector>
#include <Eigen/SVD>
#include "../WAssert.h"
#include "../WLimits.h"
......@@ -292,34 +290,3 @@ bool linearIndependent( const WVector3d& u, const WVector3d& v )
}
return true;
}
void computeSVD( const WMatrix<double>& A,
WMatrix<double>& U,
WMatrix<double>& V,
WValue<double>& S )
{
Eigen::MatrixXd eigenA( A );
Eigen::JacobiSVD<Eigen::MatrixXd> svd( eigenA, Eigen::ComputeFullU | Eigen::ComputeFullV );
U = WMatrix<double>( svd.matrixU() );
V = WMatrix<double>( svd.matrixV() );
S = WValue<double>( svd.singularValues() );
}
WMatrix<double> pseudoInverse( const WMatrix<double>& input )
{
// calc pseudo inverse
WMatrix<double> U( input.getNbRows(), input.getNbCols() );
WMatrix<double> V( input.getNbCols(), input.getNbCols() );
WValue<double> Svec( input.size() );
computeSVD( input, U, V, Svec );
// create diagonal matrix S
WMatrix<double> S( input.getNbCols(), input.getNbCols() );
S.setZero();
for( size_t i = 0; i < Svec.size() && i < S.getNbRows() && i < S.getNbCols(); i++ )
{
S( i, i ) = ( Svec[ i ] == 0.0 ) ? 0.0 : 1.0 / Svec[ i ];
}
return WMatrix<double>( V*S*U.transposed() );
}
......@@ -25,6 +25,7 @@
#ifndef WLINEARALGEBRAFUNCTIONS_H
#define WLINEARALGEBRAFUNCTIONS_H
#include <Eigen/SVD>
#include "WMatrix.h"
#include "linearAlgebra/WLinearAlgebra.h"
......@@ -92,22 +93,60 @@ bool linearIndependent( const WVector3d& u, const WVector3d& v );
*
* A=U*S*V^T
*
* \tparam T The data type.
* \param A Input Matrix
* \param U Output Matrix
* \param S Output of the entries in the diagonal matrix
* \param V Output Matrix
*
*/
void computeSVD( const WMatrix<double>& A, WMatrix<double>& U, WMatrix<double>& V, WValue<double>& S );
template< typename T >
void computeSVD( const WMatrix< T >& A, WMatrix< T >& U, WMatrix< T >& V, WValue< T >& S );
/**
* Calculates for a matrix the pseudo inverse.
*
* \tparam T The data type.
* \param input Matrix to invert
*
* \return Inverted Matrix
*
*/
WMatrix<double> pseudoInverse( const WMatrix<double>& input );
template< typename T >
WMatrix< T > pseudoInverse( const WMatrix< T >& input );
template< typename T >
void computeSVD( const WMatrix< T >& A,
WMatrix< T >& U,
WMatrix< T >& V,
WValue< T >& S )
{
Eigen::Matrix< T, -1, -1 > eigenA( A );
Eigen::JacobiSVD< Eigen::Matrix< T, -1, -1 > > svd( eigenA, Eigen::ComputeFullU | Eigen::ComputeFullV );
U = WMatrix< T >( svd.matrixU() );
V = WMatrix< T >( svd.matrixV() );
S = WValue< T >( svd.singularValues() );
}
template< typename T >
WMatrix< T > pseudoInverse( const WMatrix< T >& input )
{
// calc pseudo inverse
WMatrix< T > U( input.getNbRows(), input.getNbCols() );
WMatrix< T > V( input.getNbCols(), input.getNbCols() );
WValue< T > Svec( input.size() );
computeSVD( input, U, V, Svec );
// create diagonal matrix S
WMatrix< T > S( input.getNbCols(), input.getNbCols() );
S.setZero();
for( size_t i = 0; i < Svec.size() && i < S.getNbRows() && i < S.getNbCols(); i++ )
{
S( i, i ) = ( Svec[ i ] == 0.0 ) ? 0.0 : 1.0 / Svec[ i ];
}
return WMatrix< T >( V*S*U.transposed() );
}
#endif // WLINEARALGEBRAFUNCTIONS_H
......@@ -79,6 +79,20 @@ public:
*/
WMatrix( const Eigen::MatrixXd& newMatrix ); // NOLINT
/**
* Copies the specified Eigen::MatrixXf.
*
* \param newMatrix the Eigen::MatrixXf matrix to copy
*/
WMatrix( const Eigen::MatrixXf& newMatrix ); // NOLINT
/**
* Copies the specified Eigen::MatrixXi.
*
* \param newMatrix the Eigen::MatrixXi matrix to copy
*/
WMatrix( const Eigen::MatrixXi& newMatrix ); // NOLINT
/**
* Makes the matrix contain the identity matrix, i.e. 1 on the diagonal.
* \return Reference to the current matrix which is identity matrix now.
......@@ -130,11 +144,14 @@ public:
operator osg::Matrixd() const;
/**
* Cast this matrix to an Eigen::MatrixXd matrix.
* Cast this matrix to an Eigen::Matrix< EigenDataType, -1, -1 >() matrix.
*
* \tparam EigenDataType Data type of Eigen matrix.
*
* \return casted matrix.
*/
operator Eigen::MatrixXd() const;
template< typename EigenDataType >
operator Eigen::Matrix< EigenDataType, -1, -1 >() const;
/**
* Compares two matrices and returns true if they are equal.
......@@ -210,6 +227,14 @@ public:
protected:
private:
/**
* This function is used by the constructors that have the different Eigen::MatrixX types as parameter.
* \tparam EigenDataType The data type which is used by the Eigen matrix.
* \param newMatrix The source matrix.
*/
template< typename EigenDataType >
void copyFromEigenMatrix( const Eigen::Matrix< EigenDataType, -1, -1 >& newMatrix );
size_t m_nbCols; //!< Number of columns of the matrix. The number of rows will be computed by (size/m_nbCols).
};
......@@ -251,16 +276,20 @@ template< typename T > WMatrix< T >::WMatrix( const WMatrix4d& newMatrix )
template< typename T > WMatrix< T >::WMatrix( const Eigen::MatrixXd& newMatrix )
: WValue< T >( newMatrix.cols() * newMatrix.rows() )
{
m_nbCols = static_cast< size_t >( newMatrix.cols() );
for( int row = 0; row < newMatrix.rows(); ++row )
{
for( int col = 0; col < newMatrix.cols(); ++col )
{
( *this )( row, col ) = static_cast< T >( newMatrix( row, col ) );
}
}
copyFromEigenMatrix< double >( newMatrix );
}
template< typename T > WMatrix< T >::WMatrix( const Eigen::MatrixXf& newMatrix )
: WValue< T >( newMatrix.cols() * newMatrix.rows() )
{
copyFromEigenMatrix< float >( newMatrix );
}
template< typename T > WMatrix< T >::WMatrix( const Eigen::MatrixXi& newMatrix )
: WValue< T >( newMatrix.cols() * newMatrix.rows() )
{
copyFromEigenMatrix< int >( newMatrix );
}
template< typename T > WMatrix< T >::operator WMatrix4d() const
{
......@@ -301,19 +330,21 @@ template< typename T > WMatrix< T >::operator osg::Matrixd() const
}
}
template< typename T > WMatrix< T >::operator Eigen::MatrixXd() const
template< typename T >
template< typename EigenDataType > WMatrix< T >::operator Eigen::Matrix< EigenDataType, -1, -1 >() const
{
Eigen::MatrixXd matrix( this->getNbRows(), this->getNbCols() );
Eigen::Matrix< EigenDataType, -1, -1 > matrix( this->getNbRows(), this->getNbCols() );
for( int row = 0; row < matrix.rows(); ++row )
{
for( int col = 0; col < matrix.cols(); ++col )
{
matrix( row, col ) = ( *this )( row, col );
matrix( row, col ) = static_cast< EigenDataType >( ( *this )( row, col ) );
}
}
return matrix;
}
/**
* Makes the matrix contain the identity matrix, i.e. 1 on the diagonal.
*/
......@@ -494,6 +525,19 @@ template< typename T > bool WMatrix< T >::isIdentity( T delta ) const
return true;
}
template< typename T >
template< typename EigenDataType > void WMatrix< T >::copyFromEigenMatrix( const Eigen::Matrix< EigenDataType, -1, -1 >& newMatrix )
{
m_nbCols = static_cast< size_t >( newMatrix.cols() );
for( int row = 0; row < newMatrix.rows(); ++row )
{
for( int col = 0; col < newMatrix.cols(); ++col )
{
( *this )( row, col ) = static_cast< T >( newMatrix( row, col ) );
}
}
}
template< typename T >
inline std::ostream& operator<<( std::ostream& os, const WMatrix< T >& m )
{
......
......@@ -726,18 +726,18 @@ WMatrix< T > WSymmetricSphericalHarmonic< T >::calcFRTMatrix( size_t order )
template< typename T >
WMatrix< T > WSymmetricSphericalHarmonic< T >::calcSHToTensorSymMatrix( std::size_t order )
{
std::vector< WMatrixFixed< T, 3, 1 > > vertices;
std::vector< WVector3d > vertices;
std::vector< unsigned int > triangles;
// calc necessary tesselation level
int level = static_cast< int >( log( static_cast< float >( ( order + 1 ) * ( order + 2 ) ) ) / 2.0f + 1.0f );
tesselateIcosahedron( &vertices, &triangles, level );
std::vector< WUnitSphereCoordinates< T > > orientations;
for( typename std::vector< WMatrixFixed< T, 3, 1 > >::const_iterator cit = vertices.begin(); cit != vertices.end(); cit++ )
for( typename std::vector< WVector3d >::const_iterator cit = vertices.begin(); cit != vertices.end(); cit++ )
{
// avoid linear dependent orientations
if( ( *cit )[ 0 ] >= 0.0001 )
{
orientations.push_back( WUnitSphereCoordinates< T >( *cit ) );
orientations.push_back( WUnitSphereCoordinates< T >( WMatrixFixed< T, 3, 1 >( *cit ) ) );
}
}
return WSymmetricSphericalHarmonic< T >::calcSHToTensorSymMatrix( order, orientations );
......
......@@ -34,6 +34,8 @@
#include <boost/array.hpp>
#include <Eigen/Dense>
#include "../WAssert.h"
#include "../WLimits.h"
#include "WCompileTimeFunctions.h"
......@@ -187,6 +189,38 @@ void jacobiEigenvector3D( WTensorSym< 2, 3, Data_T > const& mat, RealEigenSystem
WAssert( iter >= 0, "Jacobi eigenvector iteration did not converge." );
}
/**
* Compute all eigenvalues as well as the corresponding eigenvectors of a
* symmetric real Matrix.
*
* \note Data_T must be castable to double.
*
* \param[in] mat A real symmetric matrix.
* \param[out] RealEigenSystem A pointer to an RealEigenSystem.
*/
template< typename Data_T >
void eigenEigenvector3D( WTensorSym< 2, 3, Data_T > const& mat, RealEigenSystem* es )
{
RealEigenSystem& result = *es; // alias for the result
typedef Eigen::Matrix< Data_T, 3, 3 > MatrixType;
MatrixType m;
m << mat( 0, 0 ), mat( 0, 1 ), mat( 0, 2 ),
mat( 1, 0 ), mat( 1, 1 ), mat( 1, 2 ),
mat( 2, 0 ), mat( 2, 1 ), mat( 2, 2 );
Eigen::SelfAdjointEigenSolver< MatrixType > solv( m );
for( std::size_t k = 0; k < 3; ++k )
{
result[ k ].first = solv.eigenvalues()[ k ];
for( std::size_t j = 0; j < 3; ++j )
{
result[ k ].second[ j ] = solv.eigenvectors()( j, k );
}
}
}
/**
* Calculate eigenvectors via the characteristic polynomial. This is essentially the same
* function as in the GPU glyph shaders. This is for 3 dimensions only.
......@@ -303,7 +337,7 @@ WTensorSym< 2, 3, T > tensorLog( WTensorSym< 2, 3, T > const& tens )
// calculate eigenvalues and eigenvectors
RealEigenSystem sys;
jacobiEigenvector3D( tens, &sys );
eigenEigenvector3D( tens, &sys );
// first, we check for negative or zero eigenvalues
if( sys[ 0 ].first <= 0.0 || sys[ 1 ].first <= 0.0 || sys[ 2 ].first <= 0.0 )
......@@ -348,7 +382,7 @@ WTensorSym< 2, 3, T > tensorExp( WTensorSym< 2, 3, T > const& tens )
// calculate eigenvalues and eigenvectors
RealEigenSystem sys;
jacobiEigenvector3D( tens, &sys );
eigenEigenvector3D( tens, &sys );
// this implements the matrix product U * exp( E ) * U.transposed()
// note that u( i, j ) = jth value of the ith eigenvector
......
......@@ -55,7 +55,7 @@ public:
* Constructor for Euclidean coordinates.
* \param vector Euclidean coordinates
*/
explicit WUnitSphereCoordinates( WMatrixFixed< T, 3, 1 > vector );
explicit WUnitSphereCoordinates( const WMatrixFixed< T, 3, 1 >& vector );
/**
* Destructor.
......@@ -95,6 +95,13 @@ public:
*/
WMatrixFixed< T, 3, 1 > getEuclidean() const;
/**
* Returns the stored sphere coordinates as Euclidean coordinates.
*
* \param vector coordinates in euclidean space
*/
void setEuclidean( WMatrixFixed< T, 3, 1 > vector );
protected:
private:
/** coordinate */
......@@ -118,12 +125,9 @@ WUnitSphereCoordinates< T >::WUnitSphereCoordinates( T theta, T phi )
}
template< typename T >
WUnitSphereCoordinates< T >::WUnitSphereCoordinates( WMatrixFixed< T, 3, 1 > vector )
WUnitSphereCoordinates< T >::WUnitSphereCoordinates( const WMatrixFixed< T, 3, 1 >& vector )
{
vector = normalize( vector );
// calculate angles
m_theta = std::acos( vector[2] );
m_phi = std::atan2( vector[1], vector[0] );
setEuclidean( vector );
}
template< typename T >
......@@ -161,4 +165,13 @@ WMatrixFixed< T, 3, 1 > WUnitSphereCoordinates< T >::getEuclidean() const
return WMatrixFixed< T, 3, 1 >( std::sin( m_theta )*std::cos( m_phi ), std::sin( m_theta )*std::sin( m_phi ), std::cos( m_theta ) );
}
template< typename T >
void WUnitSphereCoordinates< T >::setEuclidean( WMatrixFixed< T, 3, 1 > vector )
{
vector = normalize( vector );
// calculate angles
m_theta = std::acos( vector[2] );
m_phi = std::atan2( vector[1], vector[0] );
}
#endif // WUNITSPHERECOORDINATES_H
......@@ -79,16 +79,33 @@ public:
}
/**
* Create a WValue from the given WVector_2.
* \param newValues The WVector_2 with the values..
* Create a WValue from the given Eigen::VectorXd.
* \param newValues The Eigen::VectorXd with the values..
*/
explicit WValue( const WVector_2& newValues )
explicit WValue( const Eigen::VectorXd& newValues )
: m_components( static_cast< std::size_t >( newValues.size() ) )
{
for( std::size_t i = 0; i < m_components.size(); ++i )
{
m_components[ i ] = static_cast< T >( newValues( i ) );
}
copyFromEigenVector( newValues );
}
/**
* Create a WValue from the given Eigen::VectorXf.
* \param newValues The Eigen::VectorXf with the values..
*/
explicit WValue( const Eigen::VectorXf& newValues )
: m_components( static_cast< std::size_t >( newValues.size() ) )
{
copyFromEigenVector( newValues );
}
/**
* Create a WValue from the given Eigen::VectorXi.
* \param newValues The Eigen::VectorXi with the values..
*/
explicit WValue( const Eigen::VectorXi& newValues )
: m_components( static_cast< std::size_t >( newValues.size() ) )
{
copyFromEigenVector( newValues );
}
/**
......@@ -370,6 +387,20 @@ public:
protected:
private:
/**
* This function is used by the constructors that have the different Eigen::MatrixX types as parameter.
* \tparam EigenDataType The data type which is used by the Eigen::VectorX.
* \param newValues The source Eigen::VectorX.
*/
template< typename EigenDataType >
void copyFromEigenVector( const Eigen::Matrix< EigenDataType, -1, 1 >& newValues )
{
for( std::size_t i = 0; i < m_components.size(); ++i )
{
m_components[ i ] = static_cast< T >( newValues( i ) );
}
}
/**
* The components the value is composed of. This contains the actual data
*/
std::vector< T > m_components;
......
......@@ -170,7 +170,7 @@ public:
line.push_back( WPosition( 7, 8, 9 ) );
double expected = length( WPosition( 1, 2, 3 ) - WPosition( 4, 5, 6 ) ) +
length( WPosition( 4, 5, 6 ) - WPosition( 7, 8, 9 ) );
TS_ASSERT_EQUALS( expected, pathLength( line ) );
TS_ASSERT_DELTA( expected, pathLength( line ), 1e-6 );
}
/**
......
//---------------------------------------------------------------------------
//
// Project: OpenWalnut ( http://www.openwalnut.org )
//
// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
// For more information see http://www.openwalnut.org/copying
//
// This file is part of OpenWalnut.
//
// OpenWalnut 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.
//
// OpenWalnut 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 for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
//
//---------------------------------------------------------------------------
#ifndef WUNITSPHERECOORDINATES_TEST_H
#define WUNITSPHERECOORDINATES_TEST_H
#include <string>
#include <vector>
#include <algorithm>
#include <cxxtest/TestSuite.h>
#include "core/common/WLogger.h"
#include "core/common/math/WGeometryFunctions.h"
#include "core/common/math/WSymmetricSphericalHarmonic.h"
#include "../WUnitSphereCoordinates.h"
/**
* Test class for the WUnitSphereCoordinates template.
*/
class WUnitSphereCoordinatesTest : public CxxTest::TestSuite
{
public:
/**
* Test implicit conversion between sphere and euclidean coordinates.
*/
void testSphereToEuclideanConversion()
{
// 1. test case
WVector3d n( 1.0, 0.0, 0.0 );
WUnitSphereCoordinates< double > m( n );
WVector3d n2( m.getEuclidean() );
for( size_t i = 0; i < 3; i++ )
{
TS_ASSERT_DELTA( n[ i ], n2[ i ], 0.001 );
}
// 2. test case
n[ 0 ] = 1.0;
n[ 1 ] = 0.5;
n[ 2 ] = 1.0;
m.setEuclidean( n );
n2 = m.getEuclidean();
n = normalize( n );
for( size_t i = 0; i < 3; i++ )
{
TS_ASSERT_DELTA( n[ i ], n2[ i ], 0.001 );
}
}
private:
};
#endif // WUNITSPHERECOORDINATES_TEST_H
......@@ -151,7 +151,7 @@ public:
void testInsertAlmostMax( void )
{
double max = 10000.000000010001;
WHistogram2D h( -2147483648, max, -2147483648, max, 2, 2 );
WHistogram2D h( -2147483646, max, -2147483646, max, 2, 2 );
h.insert( 10000, 10000 );
h.insert( max - 2.0 * wlimits::FLT_EPS, max - 2.0 * wlimits::FLT_EPS );
TS_ASSERT_EQUALS( h( 1, 1 ), 2 );
......
......@@ -139,7 +139,7 @@ public:
void testInsertAlmostMax( void )
{
double max = 10000.000000010001;
WHistogramBasic h( -2147483648, max );
WHistogramBasic h( -2147483646, max );
h.insert( 10000 );
h.insert( max - 2.0 * wlimits::FLT_EPS );
TS_ASSERT_EQUALS( h[999], 2 );
......
......@@ -25,6 +25,8 @@
#ifndef WTHREADEDTRACKINGFUNCTION_H
#define WTHREADEDTRACKINGFUNCTION_H
#include <stdint.h>
#include <vector>
#include <utility>
......@@ -45,7 +47,11 @@ class WThreadedTrackingFunctionTest; //! forward declaration
namespace wtracking // note that this is not final
{
// an epsilon value for various floating point comparisons
#define TRACKING_EPS 0.0000001
#if INTPTR_MAX == INT32_MAX
#define TRACKING_EPS 0.00001
#else
#define TRACKING_EPS 0.0000001
#endif
/**
* \class WTrackingUtility
......
......@@ -56,13 +56,13 @@ public:
TS_ASSERT_EQUALS( v.getOffsetX(), 1.0 );
TS_ASSERT_EQUALS( v.getOffsetY(), 1.0 );
TS_ASSERT_EQUALS( v.getOffsetZ(), 1.0 );
TS_ASSERT_EQUALS( v.getUnitDirectionX(), WVector3d( 1.0, 0.0, 0.0 ) );
TS_ASSERT_EQUALS( v.getUnitDirectionY(), WVector3d( 0.0, 1.0, 0.0 ) );
TS_ASSERT_EQUALS( v.getUnitDirectionZ(), WVector3d( 0.0, 0.0, 1.0 ) );
TS_ASSERT_EQUALS( v.getDirectionX(), WVector3d( 1.0, 0.0, 0.0 ) );
TS_ASSERT_EQUALS( v.getDirectionY(), WVector3d( 0.0, 1.0, 0.0 ) );
TS_ASSERT_EQUALS( v.getDirectionZ(), WVector3d( 0.0, 0.0, 1.0 ) );
TS_ASSERT_EQUALS( v.getOrigin(), WVector3d( 0.0, 0.0, 0.0 ) );
compareVectors( v.getUnitDirectionX(), WVector3d( 1.0, 0.0, 0.0 ), 0.0001 );
compareVectors( v.getUnitDirectionY(), WVector3d( 0.0, 1.0, 0.0 ), 0.0001 );
compareVectors( v.getUnitDirectionZ(), WVector3d( 0.0, 0.0, 1.0 ), 0.0001 );
compareVectors( v.getDirectionX(), WVector3d( 1.0, 0.0, 0.0 ), 0.0001 );
compareVectors( v.getDirectionY(), WVector3d( 0.0, 1.0, 0.0 ), 0.0001 );
compareVectors( v.getDirectionZ(), WVector3d( 0.0, 0.0, 1.0 ), 0.0001 );
compareVectors( v.getOrigin(), WVector3d( 0.0, 0.0, 0.0 ), 0.0001 );
}
{
TS_ASSERT_THROWS_NOTHING( WGridTransformOrtho v( 2.2, 3.3, -1.0 ) );
......@@ -75,13 +75,13 @@ public:
TS_ASSERT_EQUALS( v.getOffsetX(), 2.2 );
TS_ASSERT_EQUALS( v.getOffsetY(), 3.3 );