OpenWalnut 1.3.1
WTensorFunctions.h
00001 //---------------------------------------------------------------------------
00002 //
00003 // Project: OpenWalnut ( http://www.openwalnut.org )
00004 //
00005 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
00006 // For more information see http://www.openwalnut.org/copying
00007 //
00008 // This file is part of OpenWalnut.
00009 //
00010 // OpenWalnut is free software: you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as published by
00012 // the Free Software Foundation, either version 3 of the License, or
00013 // (at your option) any later version.
00014 //
00015 // OpenWalnut is distributed in the hope that it will be useful,
00016 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018 // GNU Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public License
00021 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
00022 //
00023 //---------------------------------------------------------------------------
00024 
00025 #ifndef WTENSORFUNCTIONS_H
00026 #define WTENSORFUNCTIONS_H
00027 
00028 #include <algorithm>
00029 #include <cmath>
00030 #include <complex>
00031 #include <iostream>
00032 #include <utility>
00033 #include <vector>
00034 
00035 #include <boost/array.hpp>
00036 
00037 #include "../WAssert.h"
00038 #include "../WLimits.h"
00039 #include "WCompileTimeFunctions.h"
00040 #include "WTensor.h"
00041 #include "WTensorSym.h"
00042 #include "WMatrix.h"
00043 
00044 /**
00045  * An eigensystem has all eigenvalues as well as its corresponding eigenvectors. A RealEigenSystem is an EigenSystem where all
00046  * eigenvalues are real and not complex.
00047  */
00048 typedef boost::array< std::pair< double, WVector3d >, 3 > RealEigenSystem;
00049 
00050 /**
00051  * An eigensystem has all eigenvalues as well its corresponding eigenvectors.
00052  */
00053 typedef boost::array< std::pair< std::complex< double >, WVector3d >, 3 > EigenSystem;
00054 
00055 std::ostream& operator<<( std::ostream& os, const RealEigenSystem& sys );
00056 
00057 namespace
00058 {
00059     void sortRealEigenSystem( RealEigenSystem* es )
00060     {
00061         if( ( *es )[0].first > ( *es )[2].first )
00062         {
00063             std::swap( ( *es )[0], ( *es )[2] );
00064         }
00065         if( ( *es )[0].first > ( *es )[1].first )
00066         {
00067             std::swap( ( *es )[0], ( *es )[1] );
00068         }
00069         if( ( *es )[1].first > ( *es )[2].first )
00070         {
00071             std::swap( ( *es )[1], ( *es )[2] );
00072         }
00073     }
00074 }
00075 
00076 /**
00077  * Compute all eigenvalues as well as the corresponding eigenvectors of a
00078  * symmetric real Matrix.
00079  *
00080  * \note Data_T must be castable to double.
00081  *
00082  * \param[in] mat A real symmetric matrix.
00083  * \param[out] RealEigenSystem A pointer to an RealEigenSystem.
00084  */
00085 template< typename Data_T >
00086 void jacobiEigenvector3D( WTensorSym< 2, 3, Data_T > const& mat, RealEigenSystem* es )
00087 {
00088     RealEigenSystem& result = *es; // alias for the result
00089     WTensorSym< 2, 3, Data_T > in( mat );
00090     WTensor< 2, 3, Data_T > ev;
00091     ev( 0, 0 ) = ev( 1, 1 ) = ev( 2, 2 ) = 1.0;
00092 
00093     int iter = 50;
00094     Data_T evp[ 3 ];
00095     Data_T evq[ 3 ];
00096 
00097     while( iter >= 0 )
00098     {
00099         int p = 1;
00100         int q = 0;
00101 
00102         for( int i = 0; i < 2; ++i )
00103         {
00104             if( fabs( in( 2, i ) ) > fabs( in( p, q ) ) )
00105             {
00106                 p = 2;
00107                 q = i;
00108             }
00109         }
00110 
00111         // Note: If all non diagonal elements sum up to nearly zero, we may quit already!
00112         // Thereby the chosen threshold 1.0e-50 was taken arbitrarily and is just a guess.
00113         if( std::abs( in( 0, 1 ) ) + std::abs( in( 0, 2 ) ) + std::abs( in( 1, 2 ) ) < 1.0e-50 )
00114         {
00115             for( int i = 0; i < 3; ++i )
00116             {
00117                 result[i].first = in( i, i );
00118                 for( int j = 0; j < 3; ++j )
00119                 {
00120                     result[i].second[j] = static_cast< double >( ev( j, i ) );
00121                 }
00122             }
00123             sortRealEigenSystem( es );
00124             return;
00125         }
00126 
00127         Data_T r = in( q, q ) - in( p, p );
00128         Data_T o = r / ( 2.0 * in( p, q ) );
00129 
00130         Data_T t;
00131         Data_T signofo = ( o < 0.0 ? -1.0 : 1.0 );
00132         if( o * o > wlimits::MAX_DOUBLE )
00133         {
00134             t = signofo / ( 2.0 * fabs( o ) );
00135         }
00136         else
00137         {
00138             t = signofo / ( fabs( o ) + sqrt( o * o + 1.0 ) );
00139         }
00140 
00141         Data_T c;
00142 
00143         if( t * t < wlimits::DBL_EPS )
00144         {
00145             c = 1.0;
00146         }
00147         else
00148         {
00149             c = 1.0 / sqrt( t * t + 1.0 );
00150         }
00151 
00152         Data_T s = c * t;
00153 
00154         int k = 0;
00155         while( k == q || k == p )
00156         {
00157             ++k;
00158         }
00159 
00160         Data_T u = ( 1.0 - c ) / s;
00161 
00162         Data_T x = in( k, p );
00163         Data_T y = in( k, q );
00164         in( p, k ) = in( k, p ) = x - s * ( y + u * x );
00165         in( q, k ) = in( k, q ) = y + s * ( x - u * y );
00166         x = in( p, p );
00167         y = in( q, q );
00168         in( p, p ) = x - t * in( p, q );
00169         in( q, q ) = y + t * in( p, q );
00170         in( q, p ) = in( p, q ) = 0.0;
00171 
00172         for( int l = 0; l < 3; ++l )
00173         {
00174             evp[ l ] = ev( l, p );
00175             evq[ l ] = ev( l, q );
00176         }
00177         for( int l = 0; l < 3; ++l )
00178         {
00179             ev( l, p ) = c * evp[ l ] - s * evq[ l ];
00180             ev( l, q ) = s * evp[ l ] + c * evq[ l ];
00181         }
00182 
00183         --iter;
00184     }
00185     WAssert( iter >= 0, "Jacobi eigenvector iteration did not converge." );
00186 }
00187 
00188 /**
00189  * Calculate eigenvectors via the characteristic polynomial. This is essentially the same
00190  * function as in the GPU glyph shaders. This is for 3 dimensions only.
00191  *
00192  * \param m The symmetric matrix to calculate the eigenvalues from.
00193  * \return A std::vector of 3 eigenvalues in descending order (of their magnitude).
00194  */
00195 std::vector< double > getEigenvaluesCardano( WTensorSym< 2, 3 > const& m );
00196 
00197 /**
00198  * Multiply tensors of order 2. This is essentially a matrix-matrix multiplication.
00199  *
00200  * Tensors must have the same data type and dimension, however both symmetric and asymmetric
00201  * tensors (in any combination) are allowed as operands. The returned tensor is always an asymmetric tensor.
00202  *
00203  * \param one The first operand, a WTensor or WTensorSym of order 2.
00204  * \param other The second operand, a WTensor or WTensorSym of order 2.
00205  *
00206  * \return A WTensor that is the product of the operands.
00207  */
00208 template< template< std::size_t, std::size_t, typename > class TensorType1, // NOLINT brainlint can't find TensorType1
00209           template< std::size_t, std::size_t, typename > class TensorType2, // NOLINT
00210           std::size_t dim,
00211           typename Data_T >
00212 WTensor< 2, dim, Data_T > operator * ( TensorType1< 2, dim, Data_T > const& one,
00213                                        TensorType2< 2, dim, Data_T > const& other )
00214 {
00215     WTensor< 2, dim, Data_T > res;
00216 
00217     // calc res_ij = one_ik * other_kj
00218     for( std::size_t i = 0; i < dim; ++i )
00219     {
00220         for( std::size_t j = 0; j < dim; ++j )
00221         {
00222             // components are initialized to zero, so there is no need to zero them here
00223             for( std::size_t k = 0; k < dim; ++k )
00224             {
00225                 res( i, j ) += one( i, k ) * other( k, j );
00226             }
00227         }
00228     }
00229 
00230     return res;
00231 }
00232 
00233 /**
00234  * Evaluate a spherical function represented by a symmetric 4th-order tensor for a given gradient.
00235  *
00236  * \tparam Data_T The integral type used to store the tensor elements.
00237  *
00238  * \param tens The tensor representing the spherical function.
00239  * \param gradient The normalized vector that represents the gradient direction.
00240  *
00241  * \note If the gradient is not normalized, the result is undefined.
00242  */
00243 template< typename Data_T >
00244 double evaluateSphericalFunction( WTensorSym< 4, 3, Data_T > const& tens, WVector3d const& gradient )
00245 {
00246     // use symmetry to reduce computation overhead
00247     // temporaries for some of the gradient element multiplications could further reduce
00248     // computation time
00249     return gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * tens( 0, 0, 0, 0 )
00250          + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 1, 1, 1, 1 )
00251          + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 2, 2, 2, 2 )
00252          + static_cast< Data_T >( 4 ) *
00253          ( gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 1 ] * tens( 0, 0, 0, 1 )
00254          + gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 2 ] * tens( 0, 0, 0, 2 )
00255          + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 0 ] * tens( 1, 1, 1, 0 )
00256          + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 0 ] * tens( 2, 2, 2, 0 )
00257          + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 2 ] * tens( 1, 1, 1, 2 )
00258          + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 1 ] * tens( 2, 2, 2, 1 ) )
00259          + static_cast< Data_T >( 12 ) *
00260          ( gradient[ 2 ] * gradient[ 1 ] * gradient[ 0 ] * gradient[ 0 ] * tens( 2, 1, 0, 0 )
00261          + gradient[ 0 ] * gradient[ 2 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 0, 2, 1, 1 )
00262          + gradient[ 0 ] * gradient[ 1 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 0, 1, 2, 2 ) )
00263          + static_cast< Data_T >( 6 ) *
00264          ( gradient[ 0 ] * gradient[ 0 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 0, 0, 1, 1 )
00265          + gradient[ 0 ] * gradient[ 0 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 0, 0, 2, 2 )
00266          + gradient[ 1 ] * gradient[ 1 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 1, 1, 2, 2 ) );
00267 }
00268 
00269 /**
00270  * Evaluate a spherical function represented by a symmetric 2nd-order tensor for a given gradient.
00271  *
00272  * \tparam Data_T The integral type used to store the tensor elements.
00273  *
00274  * \param tens The tensor representing the spherical function.
00275  * \param gradient The normalized vector that represents the gradient direction.
00276  *
00277  * \note If the gradient is not normalized, the result is undefined.
00278  */
00279 template< typename Data_T >
00280 double evaluateSphericalFunction( WTensorSym< 2, 3, Data_T > const& tens, WVector3d const& gradient )
00281 {
00282     return gradient[ 0 ] * gradient[ 0 ] * tens( 0, 0 )
00283          + gradient[ 1 ] * gradient[ 1 ] * tens( 1, 1 )
00284          + gradient[ 2 ] * gradient[ 2 ] * tens( 2, 2 )
00285          + static_cast< Data_T >( 2 ) *
00286          ( gradient[ 0 ] * gradient[ 1 ] * tens( 0, 1 )
00287          + gradient[ 0 ] * gradient[ 2 ] * tens( 0, 2 )
00288          + gradient[ 1 ] * gradient[ 2 ] * tens( 1, 2 ) );
00289 }
00290 
00291 /**
00292  * Calculate the logarithm of the given real symmetric tensor.
00293  *
00294  * \param tens The tensor.
00295  * \return The logarithm of the tensor.
00296  */
00297 template< typename T >
00298 WTensorSym< 2, 3, T > tensorLog( WTensorSym< 2, 3, T > const& tens )
00299 {
00300     WTensorSym< 2, 3, T > res;
00301 
00302     // calculate eigenvalues and eigenvectors
00303     RealEigenSystem sys;
00304     jacobiEigenvector3D( tens, &sys );
00305 
00306     // first, we check for negative or zero eigenvalues
00307     if( sys[ 0 ].first <= 0.0 || sys[ 1 ].first <= 0.0 || sys[ 2 ].first <= 0.0 )
00308     {
00309         res( 0, 0 ) = res( 1, 1 ) = res( 2, 2 ) = 1.0;
00310         return res;
00311     }
00312 
00313     // this implements the matrix product U * log( E ) * U.transposed()
00314     // note that u( i, j ) = jth value of the ith eigenvector
00315     res( 0, 0 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 0 ] * log( sys[ 0 ].first )
00316                 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 0 ] * log( sys[ 1 ].first )
00317                 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 0 ] * log( sys[ 2 ].first );
00318     res( 0, 1 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 1 ] * log( sys[ 0 ].first )
00319                 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 1 ] * log( sys[ 1 ].first )
00320                 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 1 ] * log( sys[ 2 ].first );
00321     res( 0, 2 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 2 ] * log( sys[ 0 ].first )
00322                 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 2 ] * log( sys[ 1 ].first )
00323                 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 2 ] * log( sys[ 2 ].first );
00324     res( 1, 1 ) = sys[ 0 ].second[ 1 ] * sys[ 0 ].second[ 1 ] * log( sys[ 0 ].first )
00325                 + sys[ 1 ].second[ 1 ] * sys[ 1 ].second[ 1 ] * log( sys[ 1 ].first )
00326                 + sys[ 2 ].second[ 1 ] * sys[ 2 ].second[ 1 ] * log( sys[ 2 ].first );
00327     res( 1, 2 ) = sys[ 0 ].second[ 1 ] * sys[ 0 ].second[ 2 ] * log( sys[ 0 ].first )
00328                 + sys[ 1 ].second[ 1 ] * sys[ 1 ].second[ 2 ] * log( sys[ 1 ].first )
00329                 + sys[ 2 ].second[ 1 ] * sys[ 2 ].second[ 2 ] * log( sys[ 2 ].first );
00330     res( 2, 2 ) = sys[ 0 ].second[ 2 ] * sys[ 0 ].second[ 2 ] * log( sys[ 0 ].first )
00331                 + sys[ 1 ].second[ 2 ] * sys[ 1 ].second[ 2 ] * log( sys[ 1 ].first )
00332                 + sys[ 2 ].second[ 2 ] * sys[ 2 ].second[ 2 ] * log( sys[ 2 ].first );
00333     return res;
00334 }
00335 
00336 /**
00337  * Calculate the exponential of the given real symmetric tensor.
00338  *
00339  * \param tens The tensor.
00340  * \return The exponential of the tensor.
00341  */
00342 template< typename T >
00343 WTensorSym< 2, 3, T > tensorExp( WTensorSym< 2, 3, T > const& tens )
00344 {
00345     WTensorSym< 2, 3, T > res;
00346 
00347     // calculate eigenvalues and eigenvectors
00348     RealEigenSystem sys;
00349     jacobiEigenvector3D( tens, &sys );
00350 
00351     // this implements the matrix product U * exp( E ) * U.transposed()
00352     // note that u( i, j ) = jth value of the ith eigenvector
00353     res( 0, 0 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 0 ] * exp( sys[ 0 ].first )
00354                 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 0 ] * exp( sys[ 1 ].first )
00355                 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 0 ] * exp( sys[ 2 ].first );
00356     res( 0, 1 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 1 ] * exp( sys[ 0 ].first )
00357                 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 1 ] * exp( sys[ 1 ].first )
00358                 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 1 ] * exp( sys[ 2 ].first );
00359     res( 0, 2 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 2 ] * exp( sys[ 0 ].first )
00360                 + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 2 ] * exp( sys[ 1 ].first )
00361                 + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 2 ] * exp( sys[ 2 ].first );
00362     res( 1, 1 ) = sys[ 0 ].second[ 1 ] * sys[ 0 ].second[ 1 ] * exp( sys[ 0 ].first )
00363                 + sys[ 1 ].second[ 1 ] * sys[ 1 ].second[ 1 ] * exp( sys[ 1 ].first )
00364                 + sys[ 2 ].second[ 1 ] * sys[ 2 ].second[ 1 ] * exp( sys[ 2 ].first );
00365     res( 1, 2 ) = sys[ 0 ].second[ 1 ] * sys[ 0 ].second[ 2 ] * exp( sys[ 0 ].first )
00366                 + sys[ 1 ].second[ 1 ] * sys[ 1 ].second[ 2 ] * exp( sys[ 1 ].first )
00367                 + sys[ 2 ].second[ 1 ] * sys[ 2 ].second[ 2 ] * exp( sys[ 2 ].first );
00368     res( 2, 2 ) = sys[ 0 ].second[ 2 ] * sys[ 0 ].second[ 2 ] * exp( sys[ 0 ].first )
00369                 + sys[ 1 ].second[ 2 ] * sys[ 1 ].second[ 2 ] * exp( sys[ 1 ].first )
00370                 + sys[ 2 ].second[ 2 ] * sys[ 2 ].second[ 2 ] * exp( sys[ 2 ].first );
00371     return res;
00372 }
00373 
00374 /**
00375  * Find the maxima of a spherical function represented by a symmetric tensor.
00376  *
00377  *
00378  *
00379  */
00380 template< std::size_t order, typename Data_T >
00381 void findSymmetricSphericalFunctionMaxima( WTensorSym< order, 3, Data_T > const& tensor,
00382                                            double threshold, double minAngularSeparationCosine, double stepSize,
00383                                            std::vector< WVector3d > const& startingPoints,
00384                                            std::vector< WVector3d >& maxima ) //  NOLINT
00385 {
00386     std::vector< double > values;
00387 
00388     std::vector< WVector3d >::const_iterator it;
00389     for( it = startingPoints.begin(); it != startingPoints.end(); ++it )
00390     {
00391         WVector3d position = *it;
00392         WVector3d gradient = approximateGradient( tensor, position );
00393         int iter;
00394 
00395         for( iter = 0; iter < 100; ++iter )
00396         {
00397             WVector3d newPosition = position + stepSize * gradient;
00398             normalize( newPosition );
00399             WVector3d newGradient = approximateGradient( tensor, newPosition );
00400 
00401             double cos = dot( newGradient, gradient );
00402             if( cos > 0.985 ) // less then 10 degree
00403             {
00404                 stepSize *= 2.0;
00405             }
00406             else if( cos < 0.866 ) // more than 30 degree
00407             {
00408                 stepSize *= 0.5;
00409                 if( stepSize < 0.000001 )
00410                 {
00411                     break;
00412                 }
00413             }
00414             if( cos > 0.886 ) // less than 30 degree
00415             {
00416                 gradient = newGradient;
00417                 position = newPosition;
00418             }
00419         }
00420         if( iter != 100 )
00421         {
00422             double newFuncValue = tensor.evaluateSphericalFunction( position );
00423 
00424             bool add = true;
00425 
00426             std::vector< int > m;
00427             m.reserve( maxima.size() );
00428 
00429             for( std::size_t k = 0; k != maxima.size(); ++k )
00430             {
00431                 // find all maxima that are to close to this one
00432                 if( dot( position, maxima[ k ] ) > minAngularSeparationCosine
00433                  || dot( maxima[ k ], -1.0 * position ) > minAngularSeparationCosine )
00434                 {
00435                     m.push_back( k );
00436                     if( values[ k ] >= newFuncValue )
00437                     {
00438                         add = false;
00439                     }
00440                 }
00441             }
00442             if( add )
00443             {
00444                 maxima.push_back( position );
00445                 values.push_back( newFuncValue );
00446 
00447                 for( int k = static_cast< int >( m.size() - 1 ); k > 0; --k )
00448                 {
00449                     maxima.erase( maxima.begin() + m[ k ] );
00450                     values.erase( values.begin() + m[ k ] );
00451                 }
00452             }
00453         }
00454     }
00455 
00456     // remove maxima that are too small
00457     double max = 0;
00458     for( std::size_t k = 0; k != maxima.size(); ++k )
00459     {
00460         if( values[ k ] > max )
00461         {
00462             max = values[ k ];
00463         }
00464     }
00465 
00466     std::size_t k = 0;
00467     while( k < maxima.size() )
00468     {
00469         if( values[ k ] < threshold * max )
00470         {
00471             maxima.erase( maxima.begin() + k );
00472             values.erase( values.begin() + k );
00473         }
00474         else
00475         {
00476             ++k;
00477         }
00478     }
00479 }
00480 
00481 template< std::size_t order, typename Data_T >
00482 WVector3d approximateGradient( WTensorSym< order, 3, Data_T > const& tensor, WVector3d const& pos )
00483 {
00484     WVector3d eXY;
00485     WVector3d eZ;
00486 
00487     if( fabs( fabs( pos[ 2 ] ) - 1.0 ) < 0.001 )
00488     {
00489         eXY = WVector3d( 1.0, 0.0, 0.0 );
00490         eZ = WVector3d( 0.0, 1.0, 0.0 );
00491     }
00492     else
00493     {
00494         // this is vectorProduct( z, pos )
00495         eXY = WVector3d( -pos[ 1 ], pos[ 0 ], 0.0 );
00496         normalize( eXY );
00497         // this is vectorProduct( eXY, pos )
00498         eZ = WVector3d( eXY[ 1 ] * pos[ 2 ], -eXY[ 0 ] * pos[ 2 ], eXY[ 0 ] * pos[ 1 ] - eXY[ 1 ] * pos[ 0 ] );
00499         normalize( eZ );
00500     }
00501 
00502     double dXY = ( tensor.evaluateSphericalFunction( normalize( pos + eXY * 0.0001 ) )
00503                  - tensor.evaluateSphericalFunction( normalize( pos - eXY * 0.0001 ) ) )
00504                  / 0.0002;
00505     double dZ = ( tensor.evaluateSphericalFunction( normalize( pos + eZ * 0.0001 ) )
00506                 - tensor.evaluateSphericalFunction( normalize( pos - eZ * 0.0001 ) ) )
00507                 / 0.0002;
00508 
00509     // std::sqrt( 1.0 - z² ) = sin( acos( z ) ) = sin( theta ) in spherical coordinates
00510     double d = 1.0 - pos[ 2 ] * pos[ 2 ];
00511     if( d < 0.0 ) // avoid possible numerical problems
00512     {
00513         d = 0.0;
00514     }
00515     WVector3d res = eZ * dZ + eXY * ( dXY / std::sqrt( d ) );
00516     normalize( res );
00517     return res;
00518 }
00519 
00520 #endif  // WTENSORFUNCTIONS_H