OpenWalnut 1.3.1
|
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