25 #ifndef WTENSORFUNCTIONS_H
26 #define WTENSORFUNCTIONS_H
35 #include <boost/array.hpp>
37 #include "../WAssert.h"
38 #include "../WLimits.h"
39 #include "WCompileTimeFunctions.h"
41 #include "WTensorSym.h"
48 typedef boost::array< std::pair< double, WVector3d >, 3 > RealEigenSystem;
53 typedef boost::array< std::pair< std::complex< double >,
WVector3d >, 3 > EigenSystem;
55 std::ostream& operator<<( std::ostream& os,
const RealEigenSystem& sys );
59 void sortRealEigenSystem( RealEigenSystem* es )
61 if( ( *es )[0].first > ( *es )[2].first )
63 std::swap( ( *es )[0], ( *es )[2] );
65 if( ( *es )[0].first > ( *es )[1].first )
67 std::swap( ( *es )[0], ( *es )[1] );
69 if( ( *es )[1].first > ( *es )[2].first )
71 std::swap( ( *es )[1], ( *es )[2] );
85 template<
typename Data_T >
88 RealEigenSystem& result = *es;
91 ev( 0, 0 ) = ev( 1, 1 ) = ev( 2, 2 ) = 1.0;
102 for(
int i = 0; i < 2; ++i )
104 if( fabs( in( 2, i ) ) > fabs( in( p, q ) ) )
113 if( std::abs( in( 0, 1 ) ) + std::abs( in( 0, 2 ) ) + std::abs( in( 1, 2 ) ) < 1.0e-50 )
115 for(
int i = 0; i < 3; ++i )
117 result[i].first = in( i, i );
118 for(
int j = 0; j < 3; ++j )
120 result[i].second[j] =
static_cast< double >( ev( j, i ) );
123 sortRealEigenSystem( es );
127 Data_T r = in( q, q ) - in( p, p );
128 Data_T o = r / ( 2.0 * in( p, q ) );
131 Data_T signofo = ( o < 0.0 ? -1.0 : 1.0 );
134 t = signofo / ( 2.0 * fabs( o ) );
138 t = signofo / ( fabs( o ) + sqrt( o * o + 1.0 ) );
149 c = 1.0 / sqrt( t * t + 1.0 );
155 while( k == q || k == p )
160 Data_T u = ( 1.0 - c ) / s;
162 Data_T x = in( k, p );
163 Data_T y = in( k, q );
164 in( p, k ) = in( k, p ) = x - s * ( y + u * x );
165 in( q, k ) = in( k, q ) = y + s * ( x - u * y );
168 in( p, p ) = x - t * in( p, q );
169 in( q, q ) = y + t * in( p, q );
170 in( q, p ) = in( p, q ) = 0.0;
172 for(
int l = 0; l < 3; ++l )
174 evp[ l ] = ev( l, p );
175 evq[ l ] = ev( l, q );
177 for(
int l = 0; l < 3; ++l )
179 ev( l, p ) = c * evp[ l ] - s * evq[ l ];
180 ev( l, q ) = s * evp[ l ] + c * evq[ l ];
185 WAssert( iter >= 0,
"Jacobi eigenvector iteration did not converge." );
208 template<
template< std::
size_t, std::
size_t,
typename >
class TensorType1,
209 template< std::
size_t, std::
size_t,
typename >
class TensorType2,
213 TensorType2< 2, dim, Data_T >
const& other )
218 for( std::size_t i = 0; i < dim; ++i )
220 for( std::size_t j = 0; j < dim; ++j )
223 for( std::size_t k = 0; k < dim; ++k )
225 res( i, j ) += one( i, k ) * other( k, j );
243 template<
typename Data_T >
249 return gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * tens( 0, 0, 0, 0 )
250 + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 1, 1, 1, 1 )
251 + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 2, 2, 2, 2 )
252 +
static_cast< Data_T
>( 4 ) *
253 ( gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 1 ] * tens( 0, 0, 0, 1 )
254 + gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 2 ] * tens( 0, 0, 0, 2 )
255 + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 0 ] * tens( 1, 1, 1, 0 )
256 + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 0 ] * tens( 2, 2, 2, 0 )
257 + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 2 ] * tens( 1, 1, 1, 2 )
258 + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 1 ] * tens( 2, 2, 2, 1 ) )
259 + static_cast< Data_T >( 12 ) *
260 ( gradient[ 2 ] * gradient[ 1 ] * gradient[ 0 ] * gradient[ 0 ] * tens( 2, 1, 0, 0 )
261 + gradient[ 0 ] * gradient[ 2 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 0, 2, 1, 1 )
262 + gradient[ 0 ] * gradient[ 1 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 0, 1, 2, 2 ) )
263 + static_cast< Data_T >( 6 ) *
264 ( gradient[ 0 ] * gradient[ 0 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 0, 0, 1, 1 )
265 + gradient[ 0 ] * gradient[ 0 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 0, 0, 2, 2 )
266 + gradient[ 1 ] * gradient[ 1 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 1, 1, 2, 2 ) );
279 template<
typename Data_T >
282 return gradient[ 0 ] * gradient[ 0 ] * tens( 0, 0 )
283 + gradient[ 1 ] * gradient[ 1 ] * tens( 1, 1 )
284 + gradient[ 2 ] * gradient[ 2 ] * tens( 2, 2 )
285 +
static_cast< Data_T
>( 2 ) *
286 ( gradient[ 0 ] * gradient[ 1 ] * tens( 0, 1 )
287 + gradient[ 0 ] * gradient[ 2 ] * tens( 0, 2 )
288 + gradient[ 1 ] * gradient[ 2 ] * tens( 1, 2 ) );
297 template<
typename T >
304 jacobiEigenvector3D( tens, &sys );
307 if( sys[ 0 ].first <= 0.0 || sys[ 1 ].first <= 0.0 || sys[ 2 ].first <= 0.0 )
309 res( 0, 0 ) = res( 1, 1 ) = res( 2, 2 ) = 1.0;
315 res( 0, 0 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 0 ] * log( sys[ 0 ].first )
316 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 0 ] * log( sys[ 1 ].first )
317 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 0 ] * log( sys[ 2 ].first );
318 res( 0, 1 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 1 ] * log( sys[ 0 ].first )
319 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 1 ] * log( sys[ 1 ].first )
320 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 1 ] * log( sys[ 2 ].first );
321 res( 0, 2 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 2 ] * log( sys[ 0 ].first )
322 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 2 ] * log( sys[ 1 ].first )
323 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 2 ] * log( sys[ 2 ].first );
324 res( 1, 1 ) = sys[ 0 ].second[ 1 ] * sys[ 0 ].second[ 1 ] * log( sys[ 0 ].first )
325 + sys[ 1 ].second[ 1 ] * sys[ 1 ].second[ 1 ] * log( sys[ 1 ].first )
326 + sys[ 2 ].second[ 1 ] * sys[ 2 ].second[ 1 ] * log( sys[ 2 ].first );
327 res( 1, 2 ) = sys[ 0 ].second[ 1 ] * sys[ 0 ].second[ 2 ] * log( sys[ 0 ].first )
328 + sys[ 1 ].second[ 1 ] * sys[ 1 ].second[ 2 ] * log( sys[ 1 ].first )
329 + sys[ 2 ].second[ 1 ] * sys[ 2 ].second[ 2 ] * log( sys[ 2 ].first );
330 res( 2, 2 ) = sys[ 0 ].second[ 2 ] * sys[ 0 ].second[ 2 ] * log( sys[ 0 ].first )
331 + sys[ 1 ].second[ 2 ] * sys[ 1 ].second[ 2 ] * log( sys[ 1 ].first )
332 + sys[ 2 ].second[ 2 ] * sys[ 2 ].second[ 2 ] * log( sys[ 2 ].first );
342 template<
typename T >
349 jacobiEigenvector3D( tens, &sys );
353 res( 0, 0 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 0 ] * exp( sys[ 0 ].first )
354 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 0 ] * exp( sys[ 1 ].first )
355 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 0 ] * exp( sys[ 2 ].first );
356 res( 0, 1 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 1 ] * exp( sys[ 0 ].first )
357 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 1 ] * exp( sys[ 1 ].first )
358 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 1 ] * exp( sys[ 2 ].first );
359 res( 0, 2 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 2 ] * exp( sys[ 0 ].first )
360 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 2 ] * exp( sys[ 1 ].first )
361 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 2 ] * exp( sys[ 2 ].first );
362 res( 1, 1 ) = sys[ 0 ].second[ 1 ] * sys[ 0 ].second[ 1 ] * exp( sys[ 0 ].first )
363 + sys[ 1 ].second[ 1 ] * sys[ 1 ].second[ 1 ] * exp( sys[ 1 ].first )
364 + sys[ 2 ].second[ 1 ] * sys[ 2 ].second[ 1 ] * exp( sys[ 2 ].first );
365 res( 1, 2 ) = sys[ 0 ].second[ 1 ] * sys[ 0 ].second[ 2 ] * exp( sys[ 0 ].first )
366 + sys[ 1 ].second[ 1 ] * sys[ 1 ].second[ 2 ] * exp( sys[ 1 ].first )
367 + sys[ 2 ].second[ 1 ] * sys[ 2 ].second[ 2 ] * exp( sys[ 2 ].first );
368 res( 2, 2 ) = sys[ 0 ].second[ 2 ] * sys[ 0 ].second[ 2 ] * exp( sys[ 0 ].first )
369 + sys[ 1 ].second[ 2 ] * sys[ 1 ].second[ 2 ] * exp( sys[ 1 ].first )
370 + sys[ 2 ].second[ 2 ] * sys[ 2 ].second[ 2 ] * exp( sys[ 2 ].first );
380 template< std::
size_t order,
typename Data_T >
382 double threshold,
double minAngularSeparationCosine,
double stepSize,
383 std::vector< WVector3d >
const& startingPoints,
384 std::vector< WVector3d >& maxima )
386 std::vector< double > values;
388 std::vector< WVector3d >::const_iterator it;
389 for( it = startingPoints.begin(); it != startingPoints.end(); ++it )
392 WVector3d gradient = approximateGradient( tensor, position );
395 for( iter = 0; iter < 100; ++iter )
397 WVector3d newPosition = position + stepSize * gradient;
398 normalize( newPosition );
399 WVector3d newGradient = approximateGradient( tensor, newPosition );
401 double cos = dot( newGradient, gradient );
406 else if( cos < 0.866 )
409 if( stepSize < 0.000001 )
416 gradient = newGradient;
417 position = newPosition;
426 std::vector< int > m;
427 m.reserve( maxima.size() );
429 for( std::size_t k = 0; k != maxima.size(); ++k )
432 if( dot( position, maxima[ k ] ) > minAngularSeparationCosine
433 || dot( maxima[ k ], -1.0 * position ) > minAngularSeparationCosine )
436 if( values[ k ] >= newFuncValue )
444 maxima.push_back( position );
445 values.push_back( newFuncValue );
447 for(
int k = static_cast< int >( m.size() - 1 ); k > 0; --k )
449 maxima.erase( maxima.begin() + m[ k ] );
450 values.erase( values.begin() + m[ k ] );
458 for( std::size_t k = 0; k != maxima.size(); ++k )
460 if( values[ k ] > max )
467 while( k < maxima.size() )
469 if( values[ k ] < threshold * max )
471 maxima.erase( maxima.begin() + k );
472 values.erase( values.begin() + k );
481 template< std::
size_t order,
typename Data_T >
487 if( fabs( fabs( pos[ 2 ] ) - 1.0 ) < 0.001 )
495 eXY =
WVector3d( -pos[ 1 ], pos[ 0 ], 0.0 );
498 eZ =
WVector3d( eXY[ 1 ] * pos[ 2 ], -eXY[ 0 ] * pos[ 2 ], eXY[ 0 ] * pos[ 1 ] - eXY[ 1 ] * pos[ 0 ] );
510 double d = 1.0 - pos[ 2 ] * pos[ 2 ];
515 WVector3d res = eZ * dZ + eXY * ( dXY / std::sqrt( d ) );
520 #endif // WTENSORFUNCTIONS_H