25 #ifndef WTENSORFUNCTIONS_H
26 #define WTENSORFUNCTIONS_H
36 #include <boost/array.hpp>
39 #include <Eigen/Dense>
41 #include "../WAssert.h"
42 #include "../WLimits.h"
43 #include "WCompileTimeFunctions.h"
45 #include "WTensorSym.h"
52 typedef boost::array< std::pair< double, WVector3d >, 3 > RealEigenSystem;
57 typedef boost::array< std::pair< std::complex< double >,
WVector3d >, 3 > EigenSystem;
59 std::ostream& operator<<( std::ostream& os,
const RealEigenSystem& sys );
66 void sortRealEigenSystem( RealEigenSystem* es )
68 if( ( *es )[0].first > ( *es )[2].first )
70 std::swap( ( *es )[0], ( *es )[2] );
72 if( ( *es )[0].first > ( *es )[1].first )
74 std::swap( ( *es )[0], ( *es )[1] );
76 if( ( *es )[1].first > ( *es )[2].first )
78 std::swap( ( *es )[1], ( *es )[2] );
91 template<
typename Data_T >
94 RealEigenSystem& result = *es;
97 ev( 0, 0 ) = ev( 1, 1 ) = ev( 2, 2 ) = 1.0;
108 for(
int i = 0; i < 2; ++i )
110 if( fabs( in( 2, i ) ) > fabs( in( p, q ) ) )
119 if( std::abs( in( 0, 1 ) ) + std::abs( in( 0, 2 ) ) + std::abs( in( 1, 2 ) ) < 1.0e-50 )
121 for(
int i = 0; i < 3; ++i )
123 result[i].first = in( i, i );
124 for(
int j = 0; j < 3; ++j )
126 result[i].second[j] =
static_cast< double >( ev( j, i ) );
129 sortRealEigenSystem( es );
133 Data_T r = in( q, q ) - in( p, p );
134 Data_T o = r / ( 2.0 * in( p, q ) );
137 Data_T signofo = ( o < 0.0 ? -1.0 : 1.0 );
140 t = signofo / ( 2.0 * fabs( o ) );
144 t = signofo / ( fabs( o ) + sqrt( o * o + 1.0 ) );
155 c = 1.0 / sqrt( t * t + 1.0 );
161 while( k == q || k == p )
166 Data_T u = ( 1.0 - c ) / s;
168 Data_T x = in( k, p );
169 Data_T y = in( k, q );
170 in( p, k ) = in( k, p ) = x - s * ( y + u * x );
171 in( q, k ) = in( k, q ) = y + s * ( x - u * y );
174 in( p, p ) = x - t * in( p, q );
175 in( q, q ) = y + t * in( p, q );
176 in( q, p ) = in( p, q ) = 0.0;
178 for(
int l = 0; l < 3; ++l )
180 evp[ l ] = ev( l, p );
181 evq[ l ] = ev( l, q );
183 for(
int l = 0; l < 3; ++l )
185 ev( l, p ) = c * evp[ l ] - s * evq[ l ];
186 ev( l, q ) = s * evp[ l ] + c * evq[ l ];
191 WAssert( iter >= 0,
"Jacobi eigenvector iteration did not converge." );
203 template<
typename Data_T >
206 RealEigenSystem& result = *es;
207 typedef Eigen::Matrix< Data_T, 3, 3 > MatrixType;
210 m << mat( 0, 0 ), mat( 0, 1 ), mat( 0, 2 ),
211 mat( 1, 0 ), mat( 1, 1 ), mat( 1, 2 ),
212 mat( 2, 0 ), mat( 2, 1 ), mat( 2, 2 );
214 Eigen::SelfAdjointEigenSolver< MatrixType > solv( m );
216 for( std::size_t k = 0; k < 3; ++k )
218 result[ k ].first = solv.eigenvalues()[ k ];
219 for( std::size_t j = 0; j < 3; ++j )
221 result[ k ].second[ j ] = solv.eigenvectors()( j, k );
246 template<
template< std::
size_t, std::
size_t,
typename >
class TensorType1,
247 template< std::
size_t, std::
size_t,
typename >
class TensorType2,
251 TensorType2< 2, dim, Data_T >
const& other )
256 for( std::size_t i = 0; i < dim; ++i )
258 for( std::size_t j = 0; j < dim; ++j )
261 for( std::size_t k = 0; k < dim; ++k )
263 res( i, j ) += one( i, k ) * other( k, j );
281 template<
typename Data_T >
287 return gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * tens( 0, 0, 0, 0 )
288 + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 1, 1, 1, 1 )
289 + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 2, 2, 2, 2 )
290 +
static_cast< Data_T
>( 4 ) *
291 ( gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 1 ] * tens( 0, 0, 0, 1 )
292 + gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 2 ] * tens( 0, 0, 0, 2 )
293 + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 0 ] * tens( 1, 1, 1, 0 )
294 + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 0 ] * tens( 2, 2, 2, 0 )
295 + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 2 ] * tens( 1, 1, 1, 2 )
296 + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 1 ] * tens( 2, 2, 2, 1 ) )
297 + static_cast< Data_T >( 12 ) *
298 ( gradient[ 2 ] * gradient[ 1 ] * gradient[ 0 ] * gradient[ 0 ] * tens( 2, 1, 0, 0 )
299 + gradient[ 0 ] * gradient[ 2 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 0, 2, 1, 1 )
300 + gradient[ 0 ] * gradient[ 1 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 0, 1, 2, 2 ) )
301 + static_cast< Data_T >( 6 ) *
302 ( gradient[ 0 ] * gradient[ 0 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 0, 0, 1, 1 )
303 + gradient[ 0 ] * gradient[ 0 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 0, 0, 2, 2 )
304 + gradient[ 1 ] * gradient[ 1 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 1, 1, 2, 2 ) );
317 template<
typename Data_T >
320 return gradient[ 0 ] * gradient[ 0 ] * tens( 0, 0 )
321 + gradient[ 1 ] * gradient[ 1 ] * tens( 1, 1 )
322 + gradient[ 2 ] * gradient[ 2 ] * tens( 2, 2 )
323 +
static_cast< Data_T
>( 2 ) *
324 ( gradient[ 0 ] * gradient[ 1 ] * tens( 0, 1 )
325 + gradient[ 0 ] * gradient[ 2 ] * tens( 0, 2 )
326 + gradient[ 1 ] * gradient[ 2 ] * tens( 1, 2 ) );
335 template<
typename T >
342 eigenEigenvector3D( tens, &sys );
345 if( sys[ 0 ].first <= 0.0 || sys[ 1 ].first <= 0.0 || sys[ 2 ].first <= 0.0 )
347 res( 0, 0 ) = res( 1, 1 ) = res( 2, 2 ) = 1.0;
353 res( 0, 0 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 0 ] * log( sys[ 0 ].first )
354 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 0 ] * log( sys[ 1 ].first )
355 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 0 ] * log( sys[ 2 ].first );
356 res( 0, 1 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 1 ] * log( sys[ 0 ].first )
357 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 1 ] * log( sys[ 1 ].first )
358 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 1 ] * log( sys[ 2 ].first );
359 res( 0, 2 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 2 ] * log( sys[ 0 ].first )
360 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 2 ] * log( sys[ 1 ].first )
361 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 2 ] * log( sys[ 2 ].first );
362 res( 1, 1 ) = sys[ 0 ].second[ 1 ] * sys[ 0 ].second[ 1 ] * log( sys[ 0 ].first )
363 + sys[ 1 ].second[ 1 ] * sys[ 1 ].second[ 1 ] * log( sys[ 1 ].first )
364 + sys[ 2 ].second[ 1 ] * sys[ 2 ].second[ 1 ] * log( sys[ 2 ].first );
365 res( 1, 2 ) = sys[ 0 ].second[ 1 ] * sys[ 0 ].second[ 2 ] * log( sys[ 0 ].first )
366 + sys[ 1 ].second[ 1 ] * sys[ 1 ].second[ 2 ] * log( sys[ 1 ].first )
367 + sys[ 2 ].second[ 1 ] * sys[ 2 ].second[ 2 ] * log( sys[ 2 ].first );
368 res( 2, 2 ) = sys[ 0 ].second[ 2 ] * sys[ 0 ].second[ 2 ] * log( sys[ 0 ].first )
369 + sys[ 1 ].second[ 2 ] * sys[ 1 ].second[ 2 ] * log( sys[ 1 ].first )
370 + sys[ 2 ].second[ 2 ] * sys[ 2 ].second[ 2 ] * log( sys[ 2 ].first );
380 template<
typename T >
387 eigenEigenvector3D( tens, &sys );
391 res( 0, 0 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 0 ] * exp( sys[ 0 ].first )
392 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 0 ] * exp( sys[ 1 ].first )
393 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 0 ] * exp( sys[ 2 ].first );
394 res( 0, 1 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 1 ] * exp( sys[ 0 ].first )
395 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 1 ] * exp( sys[ 1 ].first )
396 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 1 ] * exp( sys[ 2 ].first );
397 res( 0, 2 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 2 ] * exp( sys[ 0 ].first )
398 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 2 ] * exp( sys[ 1 ].first )
399 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 2 ] * exp( sys[ 2 ].first );
400 res( 1, 1 ) = sys[ 0 ].second[ 1 ] * sys[ 0 ].second[ 1 ] * exp( sys[ 0 ].first )
401 + sys[ 1 ].second[ 1 ] * sys[ 1 ].second[ 1 ] * exp( sys[ 1 ].first )
402 + sys[ 2 ].second[ 1 ] * sys[ 2 ].second[ 1 ] * exp( sys[ 2 ].first );
403 res( 1, 2 ) = sys[ 0 ].second[ 1 ] * sys[ 0 ].second[ 2 ] * exp( sys[ 0 ].first )
404 + sys[ 1 ].second[ 1 ] * sys[ 1 ].second[ 2 ] * exp( sys[ 1 ].first )
405 + sys[ 2 ].second[ 1 ] * sys[ 2 ].second[ 2 ] * exp( sys[ 2 ].first );
406 res( 2, 2 ) = sys[ 0 ].second[ 2 ] * sys[ 0 ].second[ 2 ] * exp( sys[ 0 ].first )
407 + sys[ 1 ].second[ 2 ] * sys[ 1 ].second[ 2 ] * exp( sys[ 1 ].first )
408 + sys[ 2 ].second[ 2 ] * sys[ 2 ].second[ 2 ] * exp( sys[ 2 ].first );
418 template< std::
size_t order,
typename Data_T >
420 double threshold,
double minAngularSeparationCosine,
double stepSize,
421 std::vector< WVector3d >
const& startingPoints,
422 std::vector< WVector3d >& maxima )
424 std::vector< double > values;
426 std::vector< WVector3d >::const_iterator it;
427 for( it = startingPoints.begin(); it != startingPoints.end(); ++it )
430 WVector3d gradient = approximateGradient( tensor, position );
433 for( iter = 0; iter < 100; ++iter )
435 WVector3d newPosition = position + stepSize * gradient;
436 normalize( newPosition );
437 WVector3d newGradient = approximateGradient( tensor, newPosition );
439 double cos = dot( newGradient, gradient );
444 else if( cos < 0.866 )
447 if( stepSize < 0.000001 )
454 gradient = newGradient;
455 position = newPosition;
464 std::vector< int > m;
465 m.reserve( maxima.size() );
467 for( std::size_t k = 0; k != maxima.size(); ++k )
470 if( dot( position, maxima[ k ] ) > minAngularSeparationCosine
471 || dot( maxima[ k ], -1.0 * position ) > minAngularSeparationCosine )
474 if( values[ k ] >= newFuncValue )
482 maxima.push_back( position );
483 values.push_back( newFuncValue );
485 for(
int k = static_cast< int >( m.size() - 1 ); k > 0; --k )
487 maxima.erase( maxima.begin() + m[ k ] );
488 values.erase( values.begin() + m[ k ] );
496 for( std::size_t k = 0; k != maxima.size(); ++k )
498 if( values[ k ] > max )
505 while( k < maxima.size() )
507 if( values[ k ] < threshold * max )
509 maxima.erase( maxima.begin() + k );
510 values.erase( values.begin() + k );
519 template< std::
size_t order,
typename Data_T >
525 if( fabs( fabs( pos[ 2 ] ) - 1.0 ) < 0.001 )
533 eXY =
WVector3d( -pos[ 1 ], pos[ 0 ], 0.0 );
536 eZ =
WVector3d( eXY[ 1 ] * pos[ 2 ], -eXY[ 0 ] * pos[ 2 ], eXY[ 0 ] * pos[ 1 ] - eXY[ 1 ] * pos[ 0 ] );
548 double d = 1.0 - pos[ 2 ] * pos[ 2 ];
553 WVector3d res = eZ * dZ + eXY * ( dXY / std::sqrt( d ) );
558 #endif // WTENSORFUNCTIONS_H
Implements a symmetric tensor that has the same number of components in every direction.
A fixed size matrix class.
const double DBL_EPS
Smallest double such: 1.0 + DBL_EPS == 1.0 is still true.
Implements a tensor that has the same number of components in every direction.
Data_T evaluateSphericalFunction(WValue< Data_T > const &gradient) const
Evaluate - for a given gradient - the spherical function represented by this symmetric tensor...
const double MAX_DOUBLE
Maximum double value.