OpenWalnut  1.2.5
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
WSymmetricSphericalHarmonic.cpp
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 #include <stdint.h>
26 
27 #include <vector>
28 
29 #include <boost/math/special_functions/spherical_harmonic.hpp>
30 
31 #include "../exceptions/WPreconditionNotMet.h"
32 #include "linearAlgebra/WLinearAlgebra.h"
33 #include "WLinearAlgebraFunctions.h"
34 #include "WMatrix.h"
35 #include "WTensorSym.h"
36 #include "WUnitSphereCoordinates.h"
37 
38 #include "WSymmetricSphericalHarmonic.h"
39 
41  m_order( 0 ),
42  m_SHCoefficients( 0 )
43 {
44 }
45 
47  m_SHCoefficients( SHCoefficients )
48 {
49  // calculate order from SHCoefficients.size:
50  // - this is done by reversing the R=(l+1)*(l+2)/2 formula from the Descoteaux paper
51  double q = std::sqrt( 0.25 + 2.0 * static_cast<double>( m_SHCoefficients.size() ) ) - 1.5;
52  m_order = static_cast<size_t>( q );
53 }
54 
56 {
57 }
58 
59 double WSymmetricSphericalHarmonic::getValue( double theta, double phi ) const
60 {
61  double result = 0.0;
62  int j = 0;
63  const double rootOf2 = std::sqrt( 2.0 );
64  for( int k = 0; k <= static_cast<int>( m_order ); k+=2 )
65  {
66  // m = 1 .. k
67  for( int m = 1; m <= k; m++ )
68  {
69  j = ( k*k+k+2 ) / 2 + m;
70  result += m_SHCoefficients[ j-1 ] * rootOf2 *
71  static_cast<double>( std::pow( -1.0, m+1 ) ) * boost::math::spherical_harmonic_i( k, m, theta, phi );
72  }
73  // m = 0
74  result += m_SHCoefficients[ ( k*k+k+2 ) / 2 - 1 ] * boost::math::spherical_harmonic_r( k, 0, theta, phi );
75  // m = -k .. -1
76  for( int m = -k; m < 0; m++ )
77  {
78  j = ( k*k+k+2 ) / 2 + m;
79  result += m_SHCoefficients[ j-1 ] * rootOf2 * boost::math::spherical_harmonic_r( k, -m, theta, phi );
80  }
81  }
82  return result;
83 }
84 
86 {
87  return getValue( coordinates.getTheta(), coordinates.getPhi() );
88 }
89 
91 {
92  return m_SHCoefficients;
93 }
94 
96 {
98  size_t k = 0;
99  for( int l = 0; l <= static_cast< int >( m_order ); l += 2 )
100  {
101  for( int m = -l; m <= l; ++m )
102  {
103  res[ k ] = m_SHCoefficients[ k ];
104  if( m < 0 && ( ( -m ) % 2 == 1 ) )
105  {
106  res[ k ] *= -1.0;
107  }
108  else if( m > 0 && ( m % 2 == 0 ) )
109  {
110  res[ k ] *= -1.0;
111  }
112  ++k;
113  }
114  }
115  return res;
116 }
117 
119 {
121  size_t k = 0;
122  double r = 1.0 / sqrt( 2.0 );
123  std::complex< double > i( 0.0, -1.0 );
124  for( int l = 0; l <= static_cast< int >( m_order ); l += 2 )
125  {
126  for( int m = -l; m <= l; ++m )
127  {
128  if( m == 0 )
129  {
130  res[ k ] = m_SHCoefficients[ k ];
131  }
132  else if( m < 0 )
133  {
134  res[ k ] += i * r * m_SHCoefficients[ k - 2 * m ];
135  res[ k ] += ( -m % 2 == 1 ? -r : r ) * m_SHCoefficients[ k ];
136  }
137  else if( m > 0 )
138  {
139  res[ k ] += i * ( -m % 2 == 0 ? -r : r ) * m_SHCoefficients[ k ];
140  res[ k ] += r * m_SHCoefficients[ k - 2 * m ];
141  }
142  ++k;
143  }
144  }
145  return res;
146 }
147 
148 double WSymmetricSphericalHarmonic::calcGFA( std::vector< WUnitSphereCoordinates > const& orientations ) const
149 {
150  double n = static_cast< double >( orientations.size() );
151  double d = 0.0;
152  double gfa = 0.0;
153  double mean = 0.0;
154  double v[ 15 ];
155 
156  for( std::size_t i = 0; i < orientations.size(); ++i )
157  {
158  v[ i ] = getValue( orientations[ i ] );
159  mean += v[ i ];
160  }
161  mean /= n;
162 
163  for( std::size_t i = 0; i < orientations.size(); ++i )
164  {
165  double f = v[ i ];
166  // we use gfa as a temporary here
167  gfa += f * f;
168  f -= mean;
169  d += f * f;
170  }
171 
172  if( gfa == 0.0 || n <= 1.0 )
173  {
174  return 0.0;
175  }
176  // this is the real gfa
177  gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
178 
179  WAssert( gfa >= -0.01 && gfa <= 1.01, "" );
180  if( gfa < 0.0 )
181  {
182  return 0.0;
183  }
184  else if( gfa > 1.0 )
185  {
186  return 1.0;
187  }
188  return gfa;
189 }
190 
192 {
193  // sh coeffs
194  WMatrix< double > s( B.getNbCols(), 1 );
195  WAssert( B.getNbCols() == m_SHCoefficients.size(), "" );
196  for( std::size_t k = 0; k < B.getNbCols(); ++k )
197  {
198  s( k, 0 ) = m_SHCoefficients[ k ];
199  }
200  s = B * s;
201  WAssert( s.getNbRows() == B.getNbRows(), "" );
202  WAssert( s.getNbCols() == 1u, "" );
203 
204  double n = static_cast< double >( s.getNbRows() );
205  double d = 0.0;
206  double gfa = 0.0;
207  double mean = 0.0;
208  for( std::size_t i = 0; i < s.getNbRows(); ++i )
209  {
210  mean += s( i, 0 );
211  }
212  mean /= n;
213 
214  for( std::size_t i = 0; i < s.getNbRows(); ++i )
215  {
216  double f = s( i, 0 );
217  // we use gfa as a temporary here
218  gfa += f * f;
219  f -= mean;
220  d += f * f;
221  }
222 
223  if( gfa == 0.0 || n <= 1.0 )
224  {
225  return 0.0;
226  }
227  // this is the real gfa
228  gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
229 
230  WAssert( gfa >= -0.01 && gfa <= 1.01, "" );
231  if( gfa < 0.0 )
232  {
233  return 0.0;
234  }
235  else if( gfa > 1.0 )
236  {
237  return 1.0;
238  }
239  return gfa;
240 }
241 
243 {
244  WAssert( frtMat.getNbCols() == m_SHCoefficients.size(), "" );
245  WAssert( frtMat.getNbRows() == m_SHCoefficients.size(), "" );
246  // Funk-Radon-Transformation as in Descoteaux's thesis
247  for( size_t j = 0; j < m_SHCoefficients.size(); j++ )
248  {
249  m_SHCoefficients[ j ] *= frtMat( j, j );
250  }
251 }
252 
254 {
255  return m_order;
256 }
257 
258 WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< WVector3d >& orientations,
259  int order,
260  double lambda,
261  bool withFRT )
262 {
263  // convert euclid 3d orientations/gradients to sphere coordinates
264  std::vector< WUnitSphereCoordinates > sphereCoordinates;
265  for( std::vector< WVector3d >::const_iterator it = orientations.begin(); it != orientations.end(); it++ )
266  {
267  sphereCoordinates.push_back( WUnitSphereCoordinates( *it ) );
268  }
269  return WSymmetricSphericalHarmonic::getSHFittingMatrix( sphereCoordinates, order, lambda, withFRT );
270 }
271 
272 WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
273  int order,
274  double lambda,
275  bool withFRT )
276 {
278  WMatrix<double> Bt( B.transposed() );
279  WMatrix<double> result( Bt*B );
280  if( lambda != 0.0 )
281  {
283  result += lambda*L;
284  }
285 
286  result = pseudoInverse( result )*Bt;
287 // result *= Bt;
288  if( withFRT )
289  {
291  return ( P * result );
292  }
293  return result;
294 }
295 
296 WMatrix<double> WSymmetricSphericalHarmonic::calcBaseMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
297  int order )
298 {
299  WMatrix<double> B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
300 
301  // calc B Matrix like in the 2007 Descoteaux paper ("Regularized, Fast, and Robust Analytical Q-Ball Imaging")
302  int j = 0;
303  const double rootOf2 = std::sqrt( 2.0 );
304  for(uint32_t row = 0; row < orientations.size(); row++ )
305  {
306  const double theta = orientations[row].getTheta();
307  const double phi = orientations[row].getPhi();
308  for( int k = 0; k <= order; k+=2 )
309  {
310  // m = 1 .. k
311  for( int m = 1; m <= k; m++ )
312  {
313  j = ( k * k + k + 2 ) / 2 + m;
314  B( row, j-1 ) = rootOf2 * static_cast<double>( std::pow( static_cast<double>( -1 ), m + 1 ) )
315  * boost::math::spherical_harmonic_i( k, m, theta, phi );
316  }
317  // m = 0
318  B( row, ( k * k + k + 2 ) / 2 - 1 ) = boost::math::spherical_harmonic_r( k, 0, theta, phi );
319  // m = -k .. -1
320  for( int m = -k; m < 0; m++ )
321  {
322  j = ( k * k + k + 2 ) / 2 + m;
323  B( row, j-1 ) = rootOf2*boost::math::spherical_harmonic_r( k, -m, theta, phi );
324  }
325  }
326  }
327  return B;
328 }
329 
331 WSymmetricSphericalHarmonic::calcComplexBaseMatrix( std::vector< WUnitSphereCoordinates > const& orientations, int order )
332 {
333  WMatrix< std::complex< double > > B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
334 
335  for( std::size_t row = 0; row < orientations.size(); row++ )
336  {
337  const double theta = orientations[ row ].getTheta();
338  const double phi = orientations[ row ].getPhi();
339 
340  int j = 0;
341  for( int k = 0; k <= order; k += 2 )
342  {
343  for( int m = -k; m < 0; m++ )
344  {
345  B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
346  ++j;
347  }
348  B( row, j ) = boost::math::spherical_harmonic( k, 0, theta, phi );
349  ++j;
350  for( int m = 1; m <= k; m++ )
351  {
352  B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
353  ++j;
354  }
355  }
356  }
357  return B;
358 }
359 
361 {
362  size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
363  std::size_t i = 0;
364  WMatrix<double> L( R, R );
365  for( size_t k = 0; k <= order; k += 2 )
366  {
367  for( int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
368  {
369  L( i, i ) = static_cast< double > ( k * k * ( k + 1 ) * ( k + 1 ) );
370  ++i;
371  }
372  }
373  return L;
374 }
375 
377 {
378  size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
379  std::size_t i = 0;
380  WMatrix< double > result( R, R );
381  for( size_t k = 0; k <= order; k += 2 )
382  {
383  double h = 2.0 * piDouble * static_cast< double >( std::pow( static_cast< double >( -1 ), static_cast< double >( k / 2 ) ) ) *
384  static_cast< double >( oddFactorial( ( k <= 1 ) ? 1 : k - 1 ) ) / static_cast<double>( evenFactorial( k ) );
385  for( int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
386  {
387  result( i, i ) = h;
388  ++i;
389  }
390  }
391  return result;
392 }
393 
395  const std::vector< WUnitSphereCoordinates >& orientations )
396 {
397  std::size_t numElements = ( order + 1 ) * ( order + 2 ) / 2;
398  WPrecondEquals( order % 2, 0u );
399  WPrecondLess( numElements, orientations.size() + 1 );
400 
401  // store first numElements orientations as 3d-vectors
402  std::vector< WVector3d > directions( numElements );
403  for( std::size_t i = 0; i < numElements; ++i )
404  {
405  directions[ i ] = orientations[ i ].getEuclidean();
406  }
407 
408  // initialize an index array
409  std::vector< std::size_t > indices( order, 0 );
410 
411  // calc tensor evaluation matrix
412  WMatrix< double > TEMat( numElements, numElements );
413  for( std::size_t j = 0; j < numElements; ++j ) // j is the 'permutation index'
414  {
415  // stores how often each value is represented in the index array
416  std::size_t amount[ 3 ] = { 0, 0, 0 };
417  for( std::size_t k = 0; k < order; ++k )
418  {
419  ++amount[ indices[ k ] ];
420  }
421 
422  // from WTensorSym.h
423  std::size_t multiplicity = calcSupersymmetricTensorMultiplicity( order, amount[ 0 ], amount[ 1 ], amount[ 2 ] );
424  for( std::size_t i = 0; i < numElements; ++i ) // i is the 'direction index'
425  {
426  TEMat( i, j ) = multiplicity;
427  for( std::size_t k = 0; k < order; ++k )
428  {
429  TEMat( i, j ) *= directions[ i ][ indices[ k ] ];
430  }
431  }
432 
433  // from TensorBase.h
434  positionIterateSortedOneStep( order, 3, indices );
435  }
436  directions.clear();
437 
438  // we do not want more orientations than nessessary
439  std::vector< WUnitSphereCoordinates > ori2( orientations.begin(), orientations.begin() + numElements );
440 
441  WMatrix< double > p = pseudoInverse( TEMat );
442 
443  return p * calcBaseMatrix( ori2, order );
444 }
445 
447 {
448  double scale = 0.0;
449  if( m_SHCoefficients.size() > 0 )
450  scale = std::sqrt( 4.0 * piDouble ) * m_SHCoefficients[0];
451  for( size_t i = 0; i < m_SHCoefficients.size(); i++ )
452  {
453  m_SHCoefficients[ i ] /= scale;
454  }
455 }
456 // NOLINT