Commit 673bb7af by Stefan Philips

[CHANGE] Renamed function calcBMatrix to calcBaseMatrix

parent c48a0694
......@@ -22,16 +22,14 @@
//
//---------------------------------------------------------------------------
#include<vector>
#include <boost/math/special_functions/spherical_harmonic.hpp>
#include "stdint.h"
#include "WLinearAlgebraFunctions.h"
#include "WUnitSphereCoordinates.h"
#include "WValue.h"
#include "WMatrix.h"
#include "WValue.h"
#include "WSymmetricSphericalHarmonic.h"
......@@ -119,6 +117,7 @@ const wmath::WValue<double>& WSymmetricSphericalHarmonic::getCoefficients() cons
void WSymmetricSphericalHarmonic::applyFunkRadonTransformation()
{
// Funk-Radon-Transformation like in the 2007 Descoteaux paper
// TODO(philips): implement
size_t lj = 0;
// double ljDbl = 0.0;
double factor = 0.0;
......@@ -148,9 +147,9 @@ size_t WSymmetricSphericalHarmonic::getOrder() const
void WSymmetricSphericalHarmonic::calcLj( size_t order )
{
m_lj.clear();
for ( size_t k = 0; k <= order; k += 2 )
for ( size_t k=0; k<=order; k+=2 )
{
for ( int m = -k; m <= static_cast<int>( k ); m++ )
for ( int m=-k; m<=static_cast<int>( k ); m++ )
{
m_lj.push_back( k );
}
......@@ -158,10 +157,7 @@ void WSymmetricSphericalHarmonic::calcLj( size_t order )
}
wmath::WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< wmath::WVector3D >& orientations,
int order,
double lambda,
bool withFRT )
wmath::WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< wmath::WVector3D >& orientations, int order, double lambda, bool withFRT )
{
// convert euclid 3d orientations/gradients to sphere coordinates
std::vector< wmath::WUnitSphereCoordinates > sphereCoordinates;
......@@ -172,18 +168,15 @@ wmath::WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const st
return WSymmetricSphericalHarmonic::getSHFittingMatrix( sphereCoordinates, order, lambda, withFRT );
}
wmath::WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< wmath::WUnitSphereCoordinates >& orientations,
int order,
double lambda,
bool withFRT )
wmath::WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< wmath::WUnitSphereCoordinates >& orientations, int order, double lambda, bool withFRT )
{
wmath::WMatrix<double> B( WSymmetricSphericalHarmonic::calcBMatrix( orientations, order ) );
wmath::WMatrix<double> B( WSymmetricSphericalHarmonic::calcBaseMatrix( orientations, order ) );
wmath::WMatrix<double> Bt( B.transposed() );
wmath::WMatrix<double> result( Bt*B );
if ( lambda != 0.0 )
{
wmath::WMatrix<double> L( WSymmetricSphericalHarmonic::calcSmoothingMatrix( order ) );
result += lambda*L;
result += lambda*L;
}
result = wmath::pseudoInverse( result )*Bt;
// result *= Bt;
......@@ -195,9 +188,9 @@ wmath::WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const st
return result;
}
wmath::WMatrix<double> WSymmetricSphericalHarmonic::calcBMatrix( const std::vector< wmath::WUnitSphereCoordinates >& orientations, int order )
wmath::WMatrix<double> WSymmetricSphericalHarmonic::calcBaseMatrix( const std::vector< wmath::WUnitSphereCoordinates >& orientations, int order )
{
wmath::WMatrix<double> B( orientations.size(), ( ( order+1 ) * ( order+2 ) ) / 2 );
wmath::WMatrix<double> B( orientations.size(), ((order+1)*(order+2))/2 );
// calc B Matrix like in the 2007 Descoteaux paper ("Regularized, Fast, and Robust Analytical Q-Ball Imaging")
int j = 0;
......@@ -212,15 +205,14 @@ wmath::WMatrix<double> WSymmetricSphericalHarmonic::calcBMatrix( const std::vect
for ( int m = 1; m <= k; m++ )
{
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 );
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 );
}
// m = 0
B( row, ( k*k+k+2 ) / 2 - 1 ) = boost::math::spherical_harmonic_r( k, 0, theta, phi );
B( row, (k*k+k+2)/2-1 ) = boost::math::spherical_harmonic_r( k, 0, theta, phi );
// m = -k .. -1
for ( int m = -k; m < 0; m++ )
{
j = ( k*k+k+2 ) / 2 + m;
j = (k*k+k+2)/2+m;
B( row, j-1 ) = rootOf2*boost::math::spherical_harmonic_r( k, -m, theta, phi );
}
}
......@@ -228,7 +220,7 @@ wmath::WMatrix<double> WSymmetricSphericalHarmonic::calcBMatrix( const std::vect
return B;
// old implementation
// wmath::WMatrix<double> B( orientations.size(), ((order+1)*(order+2))/2 );
//
//
// // calc B Matrix like in the 2007 Descoteaux paper ("Regularized, Fast, and Robust Analytical Q-Ball Imaging")
// int j = 0;
// for (uint row = 0; row < orientations.size(); row++ )
......@@ -261,14 +253,14 @@ wmath::WMatrix<double> WSymmetricSphericalHarmonic::calcSmoothingMatrix( size_t
{
std::vector<int> lj;
lj.clear();
for ( size_t k = 0; k <= order; k += 2 )
for ( size_t k=0; k<=order; k+=2 )
{
for ( int m =- k; m <= static_cast<int>( k ); m++ )
for ( int m=-k; m<=static_cast<int>( k ); m++ )
{
lj.push_back( k );
}
}
size_t R = ( ( order+1 ) * ( order+2 ) ) / 2;
size_t R = ( (order+1) * (order+2) ) / 2;
wmath::WMatrix<double> L( R, R );
for ( size_t i = 0; i < R; i++ ) L( i, i ) = static_cast<double> ( lj[i] * lj[i] * ( lj[i] + 1) * ( lj[i] + 1) );
return L;
......@@ -278,8 +270,9 @@ wmath::WMatrix<double> WSymmetricSphericalHarmonic::calcFRTMatrix( size_t order
{
size_t dim = ( order + 1 ) * ( order + 2 ) / 2;
if ( WSymmetricSphericalHarmonic::m_lj.size() < dim ) WSymmetricSphericalHarmonic::calcLj( order );
wmath::WMatrix<double> result( dim, dim );
double factor = 0.0;
for ( size_t j = 0; j < dim; j++ )
{
size_t lj = WSymmetricSphericalHarmonic::m_lj[j];
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment