OpenWalnut  1.3.1
WTensorFunctions.h
1 //---------------------------------------------------------------------------
2 //
3 // Project: OpenWalnut ( http://www.openwalnut.org )
4 //
5 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
6 // For more information see http://www.openwalnut.org/copying
7 //
8 // This file is part of OpenWalnut.
9 //
10 // OpenWalnut is free software: you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // OpenWalnut is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public License
21 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
22 //
23 //---------------------------------------------------------------------------
24 
25 #ifndef WTENSORFUNCTIONS_H
26 #define WTENSORFUNCTIONS_H
27 
28 #include <algorithm>
29 #include <cmath>
30 #include <complex>
31 #include <iostream>
32 #include <utility>
33 #include <vector>
34 
35 #include <boost/array.hpp>
36 
37 #include "../WAssert.h"
38 #include "../WLimits.h"
39 #include "WCompileTimeFunctions.h"
40 #include "WTensor.h"
41 #include "WTensorSym.h"
42 #include "WMatrix.h"
43 
44 /**
45  * An eigensystem has all eigenvalues as well as its corresponding eigenvectors. A RealEigenSystem is an EigenSystem where all
46  * eigenvalues are real and not complex.
47  */
48 typedef boost::array< std::pair< double, WVector3d >, 3 > RealEigenSystem;
49 
50 /**
51  * An eigensystem has all eigenvalues as well its corresponding eigenvectors.
52  */
53 typedef boost::array< std::pair< std::complex< double >, WVector3d >, 3 > EigenSystem;
54 
55 std::ostream& operator<<( std::ostream& os, const RealEigenSystem& sys );
56 
57 namespace
58 {
59  void sortRealEigenSystem( RealEigenSystem* es )
60  {
61  if( ( *es )[0].first > ( *es )[2].first )
62  {
63  std::swap( ( *es )[0], ( *es )[2] );
64  }
65  if( ( *es )[0].first > ( *es )[1].first )
66  {
67  std::swap( ( *es )[0], ( *es )[1] );
68  }
69  if( ( *es )[1].first > ( *es )[2].first )
70  {
71  std::swap( ( *es )[1], ( *es )[2] );
72  }
73  }
74 }
75 
76 /**
77  * Compute all eigenvalues as well as the corresponding eigenvectors of a
78  * symmetric real Matrix.
79  *
80  * \note Data_T must be castable to double.
81  *
82  * \param[in] mat A real symmetric matrix.
83  * \param[out] RealEigenSystem A pointer to an RealEigenSystem.
84  */
85 template< typename Data_T >
86 void jacobiEigenvector3D( WTensorSym< 2, 3, Data_T > const& mat, RealEigenSystem* es )
87 {
88  RealEigenSystem& result = *es; // alias for the result
91  ev( 0, 0 ) = ev( 1, 1 ) = ev( 2, 2 ) = 1.0;
92 
93  int iter = 50;
94  Data_T evp[ 3 ];
95  Data_T evq[ 3 ];
96 
97  while( iter >= 0 )
98  {
99  int p = 1;
100  int q = 0;
101 
102  for( int i = 0; i < 2; ++i )
103  {
104  if( fabs( in( 2, i ) ) > fabs( in( p, q ) ) )
105  {
106  p = 2;
107  q = i;
108  }
109  }
110 
111  // Note: If all non diagonal elements sum up to nearly zero, we may quit already!
112  // Thereby the chosen threshold 1.0e-50 was taken arbitrarily and is just a guess.
113  if( std::abs( in( 0, 1 ) ) + std::abs( in( 0, 2 ) ) + std::abs( in( 1, 2 ) ) < 1.0e-50 )
114  {
115  for( int i = 0; i < 3; ++i )
116  {
117  result[i].first = in( i, i );
118  for( int j = 0; j < 3; ++j )
119  {
120  result[i].second[j] = static_cast< double >( ev( j, i ) );
121  }
122  }
123  sortRealEigenSystem( es );
124  return;
125  }
126 
127  Data_T r = in( q, q ) - in( p, p );
128  Data_T o = r / ( 2.0 * in( p, q ) );
129 
130  Data_T t;
131  Data_T signofo = ( o < 0.0 ? -1.0 : 1.0 );
132  if( o * o > wlimits::MAX_DOUBLE )
133  {
134  t = signofo / ( 2.0 * fabs( o ) );
135  }
136  else
137  {
138  t = signofo / ( fabs( o ) + sqrt( o * o + 1.0 ) );
139  }
140 
141  Data_T c;
142 
143  if( t * t < wlimits::DBL_EPS )
144  {
145  c = 1.0;
146  }
147  else
148  {
149  c = 1.0 / sqrt( t * t + 1.0 );
150  }
151 
152  Data_T s = c * t;
153 
154  int k = 0;
155  while( k == q || k == p )
156  {
157  ++k;
158  }
159 
160  Data_T u = ( 1.0 - c ) / s;
161 
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 );
166  x = in( p, p );
167  y = in( q, q );
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;
171 
172  for( int l = 0; l < 3; ++l )
173  {
174  evp[ l ] = ev( l, p );
175  evq[ l ] = ev( l, q );
176  }
177  for( int l = 0; l < 3; ++l )
178  {
179  ev( l, p ) = c * evp[ l ] - s * evq[ l ];
180  ev( l, q ) = s * evp[ l ] + c * evq[ l ];
181  }
182 
183  --iter;
184  }
185  WAssert( iter >= 0, "Jacobi eigenvector iteration did not converge." );
186 }
187 
188 /**
189  * Calculate eigenvectors via the characteristic polynomial. This is essentially the same
190  * function as in the GPU glyph shaders. This is for 3 dimensions only.
191  *
192  * \param m The symmetric matrix to calculate the eigenvalues from.
193  * \return A std::vector of 3 eigenvalues in descending order (of their magnitude).
194  */
195 std::vector< double > getEigenvaluesCardano( WTensorSym< 2, 3 > const& m );
196 
197 /**
198  * Multiply tensors of order 2. This is essentially a matrix-matrix multiplication.
199  *
200  * Tensors must have the same data type and dimension, however both symmetric and asymmetric
201  * tensors (in any combination) are allowed as operands. The returned tensor is always an asymmetric tensor.
202  *
203  * \param one The first operand, a WTensor or WTensorSym of order 2.
204  * \param other The second operand, a WTensor or WTensorSym of order 2.
205  *
206  * \return A WTensor that is the product of the operands.
207  */
208 template< template< std::size_t, std::size_t, typename > class TensorType1, // NOLINT brainlint can't find TensorType1
209  template< std::size_t, std::size_t, typename > class TensorType2, // NOLINT
210  std::size_t dim,
211  typename Data_T >
212 WTensor< 2, dim, Data_T > operator * ( TensorType1< 2, dim, Data_T > const& one,
213  TensorType2< 2, dim, Data_T > const& other )
214 {
216 
217  // calc res_ij = one_ik * other_kj
218  for( std::size_t i = 0; i < dim; ++i )
219  {
220  for( std::size_t j = 0; j < dim; ++j )
221  {
222  // components are initialized to zero, so there is no need to zero them here
223  for( std::size_t k = 0; k < dim; ++k )
224  {
225  res( i, j ) += one( i, k ) * other( k, j );
226  }
227  }
228  }
229 
230  return res;
231 }
232 
233 /**
234  * Evaluate a spherical function represented by a symmetric 4th-order tensor for a given gradient.
235  *
236  * \tparam Data_T The integral type used to store the tensor elements.
237  *
238  * \param tens The tensor representing the spherical function.
239  * \param gradient The normalized vector that represents the gradient direction.
240  *
241  * \note If the gradient is not normalized, the result is undefined.
242  */
243 template< typename Data_T >
244 double evaluateSphericalFunction( WTensorSym< 4, 3, Data_T > const& tens, WVector3d const& gradient )
245 {
246  // use symmetry to reduce computation overhead
247  // temporaries for some of the gradient element multiplications could further reduce
248  // computation time
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 ) );
267 }
268 
269 /**
270  * Evaluate a spherical function represented by a symmetric 2nd-order tensor for a given gradient.
271  *
272  * \tparam Data_T The integral type used to store the tensor elements.
273  *
274  * \param tens The tensor representing the spherical function.
275  * \param gradient The normalized vector that represents the gradient direction.
276  *
277  * \note If the gradient is not normalized, the result is undefined.
278  */
279 template< typename Data_T >
280 double evaluateSphericalFunction( WTensorSym< 2, 3, Data_T > const& tens, WVector3d const& gradient )
281 {
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 ) );
289 }
290 
291 /**
292  * Calculate the logarithm of the given real symmetric tensor.
293  *
294  * \param tens The tensor.
295  * \return The logarithm of the tensor.
296  */
297 template< typename T >
298 WTensorSym< 2, 3, T > tensorLog( WTensorSym< 2, 3, T > const& tens )
299 {
301 
302  // calculate eigenvalues and eigenvectors
303  RealEigenSystem sys;
304  jacobiEigenvector3D( tens, &sys );
305 
306  // first, we check for negative or zero eigenvalues
307  if( sys[ 0 ].first <= 0.0 || sys[ 1 ].first <= 0.0 || sys[ 2 ].first <= 0.0 )
308  {
309  res( 0, 0 ) = res( 1, 1 ) = res( 2, 2 ) = 1.0;
310  return res;
311  }
312 
313  // this implements the matrix product U * log( E ) * U.transposed()
314  // note that u( i, j ) = jth value of the ith eigenvector
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 );
333  return res;
334 }
335 
336 /**
337  * Calculate the exponential of the given real symmetric tensor.
338  *
339  * \param tens The tensor.
340  * \return The exponential of the tensor.
341  */
342 template< typename T >
343 WTensorSym< 2, 3, T > tensorExp( WTensorSym< 2, 3, T > const& tens )
344 {
346 
347  // calculate eigenvalues and eigenvectors
348  RealEigenSystem sys;
349  jacobiEigenvector3D( tens, &sys );
350 
351  // this implements the matrix product U * exp( E ) * U.transposed()
352  // note that u( i, j ) = jth value of the ith eigenvector
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 );
371  return res;
372 }
373 
374 /**
375  * Find the maxima of a spherical function represented by a symmetric tensor.
376  *
377  *
378  *
379  */
380 template< std::size_t order, typename Data_T >
381 void findSymmetricSphericalFunctionMaxima( WTensorSym< order, 3, Data_T > const& tensor,
382  double threshold, double minAngularSeparationCosine, double stepSize,
383  std::vector< WVector3d > const& startingPoints,
384  std::vector< WVector3d >& maxima ) // NOLINT
385 {
386  std::vector< double > values;
387 
388  std::vector< WVector3d >::const_iterator it;
389  for( it = startingPoints.begin(); it != startingPoints.end(); ++it )
390  {
391  WVector3d position = *it;
392  WVector3d gradient = approximateGradient( tensor, position );
393  int iter;
394 
395  for( iter = 0; iter < 100; ++iter )
396  {
397  WVector3d newPosition = position + stepSize * gradient;
398  normalize( newPosition );
399  WVector3d newGradient = approximateGradient( tensor, newPosition );
400 
401  double cos = dot( newGradient, gradient );
402  if( cos > 0.985 ) // less then 10 degree
403  {
404  stepSize *= 2.0;
405  }
406  else if( cos < 0.866 ) // more than 30 degree
407  {
408  stepSize *= 0.5;
409  if( stepSize < 0.000001 )
410  {
411  break;
412  }
413  }
414  if( cos > 0.886 ) // less than 30 degree
415  {
416  gradient = newGradient;
417  position = newPosition;
418  }
419  }
420  if( iter != 100 )
421  {
422  double newFuncValue = tensor.evaluateSphericalFunction( position );
423 
424  bool add = true;
425 
426  std::vector< int > m;
427  m.reserve( maxima.size() );
428 
429  for( std::size_t k = 0; k != maxima.size(); ++k )
430  {
431  // find all maxima that are to close to this one
432  if( dot( position, maxima[ k ] ) > minAngularSeparationCosine
433  || dot( maxima[ k ], -1.0 * position ) > minAngularSeparationCosine )
434  {
435  m.push_back( k );
436  if( values[ k ] >= newFuncValue )
437  {
438  add = false;
439  }
440  }
441  }
442  if( add )
443  {
444  maxima.push_back( position );
445  values.push_back( newFuncValue );
446 
447  for( int k = static_cast< int >( m.size() - 1 ); k > 0; --k )
448  {
449  maxima.erase( maxima.begin() + m[ k ] );
450  values.erase( values.begin() + m[ k ] );
451  }
452  }
453  }
454  }
455 
456  // remove maxima that are too small
457  double max = 0;
458  for( std::size_t k = 0; k != maxima.size(); ++k )
459  {
460  if( values[ k ] > max )
461  {
462  max = values[ k ];
463  }
464  }
465 
466  std::size_t k = 0;
467  while( k < maxima.size() )
468  {
469  if( values[ k ] < threshold * max )
470  {
471  maxima.erase( maxima.begin() + k );
472  values.erase( values.begin() + k );
473  }
474  else
475  {
476  ++k;
477  }
478  }
479 }
480 
481 template< std::size_t order, typename Data_T >
482 WVector3d approximateGradient( WTensorSym< order, 3, Data_T > const& tensor, WVector3d const& pos )
483 {
484  WVector3d eXY;
485  WVector3d eZ;
486 
487  if( fabs( fabs( pos[ 2 ] ) - 1.0 ) < 0.001 )
488  {
489  eXY = WVector3d( 1.0, 0.0, 0.0 );
490  eZ = WVector3d( 0.0, 1.0, 0.0 );
491  }
492  else
493  {
494  // this is vectorProduct( z, pos )
495  eXY = WVector3d( -pos[ 1 ], pos[ 0 ], 0.0 );
496  normalize( eXY );
497  // this is vectorProduct( eXY, pos )
498  eZ = WVector3d( eXY[ 1 ] * pos[ 2 ], -eXY[ 0 ] * pos[ 2 ], eXY[ 0 ] * pos[ 1 ] - eXY[ 1 ] * pos[ 0 ] );
499  normalize( eZ );
500  }
501 
502  double dXY = ( tensor.evaluateSphericalFunction( normalize( pos + eXY * 0.0001 ) )
503  - tensor.evaluateSphericalFunction( normalize( pos - eXY * 0.0001 ) ) )
504  / 0.0002;
505  double dZ = ( tensor.evaluateSphericalFunction( normalize( pos + eZ * 0.0001 ) )
506  - tensor.evaluateSphericalFunction( normalize( pos - eZ * 0.0001 ) ) )
507  / 0.0002;
508 
509  // std::sqrt( 1.0 - z² ) = sin( acos( z ) ) = sin( theta ) in spherical coordinates
510  double d = 1.0 - pos[ 2 ] * pos[ 2 ];
511  if( d < 0.0 ) // avoid possible numerical problems
512  {
513  d = 0.0;
514  }
515  WVector3d res = eZ * dZ + eXY * ( dXY / std::sqrt( d ) );
516  normalize( res );
517  return res;
518 }
519 
520 #endif // WTENSORFUNCTIONS_H