WSymmetricSphericalHarmonic.cpp 16.9 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
//---------------------------------------------------------------------------
//
// Project: OpenWalnut ( http://www.openwalnut.org )
//
// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
// For more information see http://www.openwalnut.org/copying
//
// This file is part of OpenWalnut.
//
// OpenWalnut is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// OpenWalnut is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
//
//---------------------------------------------------------------------------

Mathias Goldau's avatar
Mathias Goldau committed
25 26
#include <stdint.h>

Alexander Wiebel's avatar
[STYLE]  
Alexander Wiebel committed
27 28
#include <vector>

29 30
#include <boost/math/special_functions/spherical_harmonic.hpp>

31
#include "core/common/WLogger.h"
32
#include "../exceptions/WPreconditionNotMet.h"
Stefan Philips's avatar
Stefan Philips committed
33
#include "linearAlgebra/WLinearAlgebra.h"
34 35
#include "WLinearAlgebraFunctions.h"
#include "WMatrix.h"
Mathias Goldau's avatar
Mathias Goldau committed
36
#include "WTensorSym.h"
37
#include "WUnitSphereCoordinates.h"
38 39

#include "WSymmetricSphericalHarmonic.h"
40 41

WSymmetricSphericalHarmonic::WSymmetricSphericalHarmonic() :
42 43
    m_order( 0 ),
    m_SHCoefficients( 0 )
44 45 46
{
}

47
WSymmetricSphericalHarmonic::WSymmetricSphericalHarmonic( const WValue<double>& SHCoefficients ) :
48 49 50 51 52 53 54 55 56 57 58 59 60 61
  m_SHCoefficients( SHCoefficients )
{
  // calculate order from SHCoefficients.size:
  // - this is done by reversing the R=(l+1)*(l+2)/2 formula from the Descoteaux paper
  double q = std::sqrt( 0.25 + 2.0 * static_cast<double>( m_SHCoefficients.size() ) ) - 1.5;
  m_order = static_cast<size_t>( q );
}

WSymmetricSphericalHarmonic::~WSymmetricSphericalHarmonic()
{
}

double WSymmetricSphericalHarmonic::getValue( double theta, double phi ) const
{
62 63 64
  double result = 0.0;
  int j = 0;
  const double rootOf2 = std::sqrt( 2.0 );
65
  for( int k = 0; k <= static_cast<int>( m_order ); k+=2 )
66 67
  {
    // m = 1 .. k
68
    for( int m = 1; m <= k; m++ )
69
    {
70 71 72
      j = ( k*k+k+2 ) / 2 + m;
      result += m_SHCoefficients[ j-1 ] * rootOf2 *
                  static_cast<double>( std::pow( -1.0, m+1 ) ) * boost::math::spherical_harmonic_i( k, m, theta, phi );
73
    }
74 75 76
    // m = 0
    result += m_SHCoefficients[ ( k*k+k+2 ) / 2 - 1 ] * boost::math::spherical_harmonic_r( k, 0, theta, phi );
    // m = -k .. -1
77
    for( int m = -k; m < 0; m++ )
78 79 80 81 82 83
    {
      j = ( k*k+k+2 ) / 2 + m;
      result += m_SHCoefficients[ j-1 ] * rootOf2 * boost::math::spherical_harmonic_r( k, -m, theta, phi );
    }
  }
  return result;
84 85 86 87 88 89 90
}

double WSymmetricSphericalHarmonic::getValue( const WUnitSphereCoordinates& coordinates ) const
{
  return getValue( coordinates.getTheta(), coordinates.getPhi() );
}

91
const WValue<double>& WSymmetricSphericalHarmonic::getCoefficients() const
92 93 94 95
{
  return m_SHCoefficients;
}

96
WValue< double > WSymmetricSphericalHarmonic::getCoefficientsSchultz() const
97
{
98
    WValue< double > res( m_SHCoefficients.size() );
99
    size_t k = 0;
100
    for( int l = 0; l <= static_cast< int >( m_order ); l += 2 )
101
    {
102
        for( int m = -l; m <= l; ++m )
103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
        {
            res[ k ] = m_SHCoefficients[ k ];
            if( m < 0 && ( ( -m ) % 2 == 1 ) )
            {
                res[ k ] *= -1.0;
            }
            else if( m > 0 && ( m % 2 == 0 ) )
            {
                res[ k ] *= -1.0;
            }
            ++k;
        }
    }
    return res;
}

119
WValue< std::complex< double > > WSymmetricSphericalHarmonic::getCoefficientsComplex() const
120
{
121
    WValue< std::complex< double > > res( m_SHCoefficients.size() );
122 123 124 125 126
    size_t k = 0;
    double r = 1.0 / sqrt( 2.0 );
    std::complex< double > i( 0.0, -1.0 );
    for( int l = 0; l <= static_cast< int >( m_order ); l += 2 )
    {
127
        for( int m = -l; m <= l; ++m )
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148
        {
            if( m == 0 )
            {
                res[ k ] = m_SHCoefficients[ k ];
            }
            else if( m < 0 )
            {
                res[ k ] += i * r * m_SHCoefficients[ k - 2 * m ];
                res[ k ] += ( -m % 2 == 1 ? -r : r ) * m_SHCoefficients[ k ];
            }
            else if( m > 0 )
            {
                res[ k ] += i * ( -m % 2 == 0 ? -r : r ) * m_SHCoefficients[ k ];
                res[ k ] += r * m_SHCoefficients[ k - 2 * m ];
            }
            ++k;
        }
    }
    return res;
}

Mathias Goldau's avatar
Mathias Goldau committed
149
double WSymmetricSphericalHarmonic::calcGFA( std::vector< WUnitSphereCoordinates > const& orientations ) const
150 151 152 153 154
{
    double n = static_cast< double >( orientations.size() );
    double d = 0.0;
    double gfa = 0.0;
    double mean = 0.0;
155
    std::vector< double > v( orientations.size() );
156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191

    for( std::size_t i = 0; i < orientations.size(); ++i )
    {
        v[ i ] = getValue( orientations[ i ] );
        mean += v[ i ];
    }
    mean /= n;

    for( std::size_t i = 0; i < orientations.size(); ++i )
    {
        double f = v[ i ];
        // we use gfa as a temporary here
        gfa += f * f;
        f -= mean;
        d += f * f;
    }

    if( gfa == 0.0 || n <= 1.0 )
    {
        return 0.0;
    }
    // this is the real gfa
    gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );

    WAssert( gfa >= -0.01 && gfa <= 1.01, "" );
    if( gfa < 0.0 )
    {
        return 0.0;
    }
    else if( gfa > 1.0 )
    {
        return 1.0;
    }
    return gfa;
}

192
double WSymmetricSphericalHarmonic::calcGFA( WMatrix< double > const& B ) const
193 194
{
    // sh coeffs
195 196 197
    WMatrix< double > s( B.getNbCols(), 1 );
    WAssert( B.getNbCols() == m_SHCoefficients.size(), "" );
    for( std::size_t k = 0; k < B.getNbCols(); ++k )
198 199 200
    {
        s( k, 0 ) = m_SHCoefficients[ k ];
    }
Stefan Philips's avatar
Stefan Philips committed
201
    s = B * s;
202 203
    WAssert( s.getNbRows() == B.getNbRows(), "" );
    WAssert( s.getNbCols() == 1u, "" );
204

205
    double n = static_cast< double >( s.getNbRows() );
206 207 208
    double d = 0.0;
    double gfa = 0.0;
    double mean = 0.0;
209
    for( std::size_t i = 0; i < s.getNbRows(); ++i )
210 211 212 213 214
    {
        mean += s( i, 0 );
    }
    mean /= n;

215
    for( std::size_t i = 0; i < s.getNbRows(); ++i )
216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242
    {
        double f = s( i, 0 );
        // we use gfa as a temporary here
        gfa += f * f;
        f -= mean;
        d += f * f;
    }

    if( gfa == 0.0 || n <= 1.0 )
    {
        return 0.0;
    }
    // this is the real gfa
    gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );

    WAssert( gfa >= -0.01 && gfa <= 1.01, "" );
    if( gfa < 0.0 )
    {
        return 0.0;
    }
    else if( gfa > 1.0 )
    {
        return 1.0;
    }
    return gfa;
}

243
void WSymmetricSphericalHarmonic::applyFunkRadonTransformation( WMatrix< double > const& frtMat )
244
{
245 246
    WAssert( frtMat.getNbCols() == m_SHCoefficients.size(), "" );
    WAssert( frtMat.getNbRows() == m_SHCoefficients.size(), "" );
247
    // Funk-Radon-Transformation as in Descoteaux's thesis
248
    for( size_t j = 0; j < m_SHCoefficients.size(); j++ )
249 250 251
    {
        m_SHCoefficients[ j ] *= frtMat( j, j );
    }
252 253 254 255 256 257
}

size_t WSymmetricSphericalHarmonic::getOrder() const
{
  return m_order;
}
258

259
WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< WVector3d >& orientations,
Alexander Wiebel's avatar
[STYLE]  
Alexander Wiebel committed
260 261 262
                                                                        int order,
                                                                        double lambda,
                                                                        bool withFRT )
263 264
{
  // convert euclid 3d orientations/gradients to sphere coordinates
Mathias Goldau's avatar
Mathias Goldau committed
265
  std::vector< WUnitSphereCoordinates > sphereCoordinates;
266
  for( std::vector< WVector3d >::const_iterator it = orientations.begin(); it != orientations.end(); it++ )
267
  {
Mathias Goldau's avatar
Mathias Goldau committed
268
    sphereCoordinates.push_back( WUnitSphereCoordinates( *it ) );
269 270 271 272
  }
  return WSymmetricSphericalHarmonic::getSHFittingMatrix( sphereCoordinates, order, lambda, withFRT );
}

273
WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
274 275 276
                                                                 int order,
                                                                 double lambda,
                                                                 bool withFRT )
277
{
278 279 280
  WMatrix<double> B( WSymmetricSphericalHarmonic::calcBaseMatrix( orientations, order ) );
  WMatrix<double> Bt( B.transposed() );
  WMatrix<double> result( Bt*B );
281
  if( lambda != 0.0 )
282
  {
283
    WMatrix<double> L( WSymmetricSphericalHarmonic::calcSmoothingMatrix( order ) );
Alexander Wiebel's avatar
[STYLE]  
Alexander Wiebel committed
284
    result += lambda*L;
285
  }
Stefan Philips's avatar
Stefan Philips committed
286

Mathias Goldau's avatar
Mathias Goldau committed
287
  result = pseudoInverse( result )*Bt;
288
  if( withFRT )
289
  {
290
    WMatrix<double> P( WSymmetricSphericalHarmonic::calcFRTMatrix( order ) );
291 292 293 294 295
    return ( P * result );
  }
  return result;
}

296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346
WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrixForConstantSolidAngle( const std::vector< WVector3d >& orientations,
                                                                                      int order,
                                                                                      double lambda )
{
  // convert euclid 3d orientations/gradients to sphere coordinates
  std::vector< WUnitSphereCoordinates > sphereCoordinates;
  for( std::vector< WVector3d >::const_iterator it = orientations.begin(); it != orientations.end(); it++ )
  {
    sphereCoordinates.push_back( WUnitSphereCoordinates( *it ) );
  }
  return WSymmetricSphericalHarmonic::getSHFittingMatrixForConstantSolidAngle( sphereCoordinates, order, lambda );
}

WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrixForConstantSolidAngle(
    const std::vector< WUnitSphereCoordinates >& orientations,
    int order,
    double lambda )
{
  WMatrix<double> B( WSymmetricSphericalHarmonic::calcBaseMatrix( orientations, order ) );
  WMatrix<double> Bt( B.transposed() );
  WMatrix<double> result( Bt*B );

  // regularisation
  if( lambda != 0.0 )
  {
    WMatrix<double> L( WSymmetricSphericalHarmonic::calcSmoothingMatrix( order ) );
    result += lambda*L;
  }

  result = pseudoInverse( result )*Bt;

  // multiply with eigenvalues of coefficients / Laplace-Beltrami operator
  WMatrix<double> LB( WSymmetricSphericalHarmonic::calcMatrixWithEigenvalues( order ) );
  wlog::debug( "" ) << "LB: " << LB;
  result = LB*result;
  wlog::debug( "" ) << "LB*result: " << result;

  // apply FRT
  WMatrix<double> P( WSymmetricSphericalHarmonic::calcFRTMatrix( order ) );
  result = P * result;
  wlog::debug( "" ) << "P: " << P;
  wlog::debug( "" ) << "P*result: " << result;

  // correction factor
  double correctionFactor = 1.0 / ( 16.0 * std::pow( piDouble, 2 ) );
  result *= correctionFactor;

  return result;
}


347
WMatrix<double> WSymmetricSphericalHarmonic::calcBaseMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
Alexander Wiebel's avatar
[STYLE]  
Alexander Wiebel committed
348
                                                                    int order )
349
{
350
  WMatrix<double> B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
351 352 353 354

  // calc B Matrix like in the 2007 Descoteaux paper ("Regularized, Fast, and Robust Analytical Q-Ball Imaging")
  int j = 0;
  const double rootOf2 = std::sqrt( 2.0 );
355
  for(uint32_t row = 0; row < orientations.size(); row++ )
356 357 358
  {
    const double theta = orientations[row].getTheta();
    const double phi = orientations[row].getPhi();
359
    for( int k = 0; k <= order; k+=2 )
360 361
    {
      // m = 1 .. k
362
      for( int m = 1; m <= k; m++ )
363
      {
Alexander Wiebel's avatar
[STYLE]  
Alexander Wiebel committed
364 365 366
        j = ( k * k + k + 2 ) / 2 + m;
        B( row, j-1 ) = rootOf2 * static_cast<double>( std::pow( static_cast<double>( -1 ), m + 1 ) )
                        * boost::math::spherical_harmonic_i( k, m, theta, phi );
367 368
      }
      // m = 0
Alexander Wiebel's avatar
[STYLE]  
Alexander Wiebel committed
369
      B( row, ( k * k + k + 2 ) / 2 - 1 ) = boost::math::spherical_harmonic_r( k, 0, theta, phi );
370
      // m = -k .. -1
371
      for( int m = -k; m < 0; m++ )
372
      {
Alexander Wiebel's avatar
[STYLE]  
Alexander Wiebel committed
373
        j = ( k * k + k + 2 ) / 2 + m;
374 375 376 377 378 379 380
        B( row, j-1 ) = rootOf2*boost::math::spherical_harmonic_r( k, -m, theta, phi );
      }
    }
  }
  return B;
}

381
WMatrix< std::complex< double > >
Mathias Goldau's avatar
Mathias Goldau committed
382
WSymmetricSphericalHarmonic::calcComplexBaseMatrix( std::vector< WUnitSphereCoordinates > const& orientations, int order )
383
{
384
    WMatrix< std::complex< double > > B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
385 386 387 388 389 390 391 392 393

    for( std::size_t row = 0; row < orientations.size(); row++ )
    {
        const double theta = orientations[ row ].getTheta();
        const double phi = orientations[ row ].getPhi();

        int j = 0;
        for( int k = 0; k <= order; k += 2 )
        {
394
            for( int m = -k; m < 0; m++ )
395 396 397 398 399 400 401 402 403 404 405 406 407 408 409
            {
                B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
                ++j;
            }
            B( row, j ) = boost::math::spherical_harmonic( k, 0, theta, phi );
            ++j;
            for( int m = 1; m <= k; m++ )
            {
                B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
                ++j;
            }
        }
    }
    return B;
}
410

411
WValue<double> WSymmetricSphericalHarmonic::calcEigenvalues( size_t order )
412
{
413 414
    size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
    std::size_t i = 0;
415
    WValue<double> result( R );
416
    for( size_t k = 0; k <= order; k += 2 )
417
    {
418
        for( int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
419
        {
420
            result[ i ] = -static_cast< double > ( k * ( k + 1 ) );
421 422
            ++i;
        }
423
    }
424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445
    return result;
}

WMatrix<double> WSymmetricSphericalHarmonic::calcMatrixWithEigenvalues( size_t order )
{
    WValue<double> eigenvalues( WSymmetricSphericalHarmonic::calcEigenvalues( order ) );
    WMatrix<double> L( eigenvalues.size(), eigenvalues.size() );
    for( std::size_t i = 0; i < eigenvalues.size(); ++i )
    {
        L( i, i ) = eigenvalues[ i ];
    }
    return L;
}

WMatrix<double> WSymmetricSphericalHarmonic::calcSmoothingMatrix( size_t order )
{
    WValue<double> eigenvalues( WSymmetricSphericalHarmonic::calcEigenvalues( order ) );
    WMatrix<double> L( eigenvalues.size(), eigenvalues.size() );
    for( std::size_t i = 0; i < eigenvalues.size(); ++i )
    {
        L( i, i ) = std::pow( eigenvalues[ i ], 2 );
    }
446
    return L;
447 448
}

449
WMatrix<double> WSymmetricSphericalHarmonic::calcFRTMatrix( size_t order )
450
{
451 452
    size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
    std::size_t i = 0;
453
    WMatrix< double > result( R, R );
454
    for( size_t k = 0; k <= order; k += 2 )
455
    {
Mathias Goldau's avatar
Mathias Goldau committed
456 457
        double h = 2.0 * piDouble * static_cast< double >( std::pow( static_cast< double >( -1 ), static_cast< double >( k / 2 ) ) ) *
                    static_cast< double >( oddFactorial( ( k <= 1 ) ? 1 : k - 1 ) ) / static_cast<double>( evenFactorial( k ) );
458
        for( int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
459 460 461 462 463 464
        {
            result( i, i ) = h;
            ++i;
        }
    }
    return result;
ledig's avatar
ledig committed
465
}
466

467 468
WMatrix< double > WSymmetricSphericalHarmonic::calcSHToTensorSymMatrix( std::size_t order,
                                                                               const std::vector< WUnitSphereCoordinates >& orientations )
469 470 471 472 473 474
{
    std::size_t numElements = ( order + 1 ) * ( order + 2 ) / 2;
    WPrecondEquals( order % 2, 0u );
    WPrecondLess( numElements, orientations.size() + 1 );

    // store first numElements orientations as 3d-vectors
475
    std::vector< WVector3d > directions( numElements );
476 477 478 479 480 481 482 483 484
    for( std::size_t i = 0; i < numElements; ++i )
    {
        directions[ i ] = orientations[ i ].getEuclidean();
    }

    // initialize an index array
    std::vector< std::size_t > indices( order, 0 );

    // calc tensor evaluation matrix
485
    WMatrix< double > TEMat( numElements, numElements );
486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511
    for( std::size_t j = 0; j < numElements; ++j ) // j is the 'permutation index'
    {
        // stores how often each value is represented in the index array
        std::size_t amount[ 3 ] = { 0, 0, 0 };
        for( std::size_t k = 0; k < order; ++k )
        {
            ++amount[ indices[ k ] ];
        }

        // from WTensorSym.h
        std::size_t multiplicity = calcSupersymmetricTensorMultiplicity( order, amount[ 0 ], amount[ 1 ], amount[ 2 ] );
        for( std::size_t i = 0; i < numElements; ++i ) // i is the 'direction index'
        {
            TEMat( i, j ) = multiplicity;
            for( std::size_t k = 0; k < order; ++k )
            {
                TEMat( i, j ) *= directions[ i ][ indices[ k ] ];
            }
        }

        // from TensorBase.h
        positionIterateSortedOneStep( order, 3, indices );
    }
    directions.clear();

    // we do not want more orientations than nessessary
Mathias Goldau's avatar
Mathias Goldau committed
512
    std::vector< WUnitSphereCoordinates > ori2( orientations.begin(), orientations.begin() + numElements );
513

514
    WMatrix< double > p = pseudoInverse( TEMat );
515 516 517

    return p * calcBaseMatrix( ori2, order );
}
Mathias Goldau's avatar
Mathias Goldau committed
518

519 520 521
void WSymmetricSphericalHarmonic::normalize()
{
  double scale = 0.0;
522
  if( m_SHCoefficients.size() > 0 )
523
    scale = std::sqrt( 4.0 * piDouble ) * m_SHCoefficients[0];
524
  for( size_t i = 0; i < m_SHCoefficients.size(); i++ )
525 526 527 528
  {
    m_SHCoefficients[ i ] /= scale;
  }
}
529
// NOLINT