Commit 39bb9b90 authored by Sebastian Eichelbaum's avatar Sebastian Eichelbaum
Browse files

[MERGE]

parents 18e9f453 262a61c9
......@@ -27,6 +27,8 @@
#include <cmath>
#include <boost/math/constants/constants.hpp>
#if defined ( _MSC_VER )
#include "float.h"
#endif
......@@ -36,6 +38,13 @@
*/
namespace wmath
{
// Pi constants - we dont use the macro M_PI, because it is not part of the C++-standard.
// ref.: http://stackoverflow.com/questions/1727881/how-to-use-the-pi-constant-in-c
/** the pi constant in float format */
const float piFloat = boost::math::constants::pi<float>();
/** the pi constant in double format */
const double piDouble = boost::math::constants::pi<double>();
/**
* Tests whether the number stored in the parameter is finite.
* \param number the number to be tested
......@@ -52,6 +61,30 @@ inline int myIsfinite( double number )
WAssert( false, "isfinite not provided on this platform or platform not known." );
#endif
}
/**
* Calculates the odd factorial. This means 1*3*5* ... * border if border is odd, or 1*3*5* ... * (border-1) if border is even.
* \param border the threshold for the factorial calculation.
*/
inline unsigned int oddFactorial( unsigned int border )
{
unsigned int result = 1;
for ( unsigned int i = 3; i <= border; i+=2 )
result *= i;
return result;
}
/**
* Calculates the even factorial. This means 2*4*6 ... * \param border if border is even, or 2*4*6* ... * ( \param border - 1 ) if border is odd.
* \param border the threshold for the factorial calculation.
*/
inline unsigned int evenFactorial( unsigned int border )
{
unsigned int result = 1;
for ( unsigned int i = 2; i <= border; i+=2 )
result *= i;
return result;
}
}
#endif // WMATH_H
//---------------------------------------------------------------------------
//
// 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/>.
//
//---------------------------------------------------------------------------
#include <boost/math/special_functions/spherical_harmonic.hpp>
#include "WUnitSphereCoordinates.h"
#include "WValue.h"
#include "WSymmetricSphericalHarmonic.h"
using wmath::WSymmetricSphericalHarmonic;
WSymmetricSphericalHarmonic::WSymmetricSphericalHarmonic() :
m_order( 0 ),
m_SHCoefficients( 0 )
{
}
WSymmetricSphericalHarmonic::WSymmetricSphericalHarmonic( const WValue<double>& SHCoefficients ) :
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
{
double result = 0.0;
int j = 0;
for ( int k = 0; k <= static_cast<int>( m_order ); k+=2 )
{
// m = 1 .. k
for ( int m = 1; m <= k; m++ )
{
j = ( k*k+k+2 ) / 2 + m;
result += m_SHCoefficients[ j-1 ] * std::sqrt( 2.0 ) * boost::math::spherical_harmonic_i( k, m, theta, phi );
}
// m = 0
result += m_SHCoefficients[ ( 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;
result += m_SHCoefficients[ j-1 ] * std::sqrt( 2.0 ) * boost::math::spherical_harmonic_r( k, m, theta, phi );
}
}
return result;
}
double WSymmetricSphericalHarmonic::getValue( const WUnitSphereCoordinates& coordinates ) const
{
return getValue( coordinates.getTheta(), coordinates.getPhi() );
}
const wmath::WValue<double>& WSymmetricSphericalHarmonic::getCoefficients() const
{
return m_SHCoefficients;
}
void WSymmetricSphericalHarmonic::applyFunkRadonTransformation()
{
// Funk-Radon-Transformation like in the 2007 Descoteaux paper
// TODO(philips): implement
WAssert( false, "implement!" );
int lj = 0.0;
double factor = 0.0;
for ( size_t j = 0; j < 2*m_SHCoefficients.size(); j++ )
{
lj = 2 * ( ( j == 1 ) ? 0 : ( size_t )( j+2 )/4 );
// lj is always even!
factor = 2.0 * wmath::piDouble * static_cast<double>( std::pow( -1, ( lj / 2 ) ) ) *
static_cast<double>( wmath::oddFactorial( lj ) ) / static_cast<double>( wmath::evenFactorial( lj ) );
// factor = (double) std::pow( -1, ( lj / 2 ) ) * (double)wmath::oddFaculty( lj ) / (double)wmath::evenFaculty( lj );
// std::cerr << "factor: " << factor << std::endl;
// m_SHCoefficients[ j ] = m_SHCoefficients[ j ] * factor;
}
}
size_t WSymmetricSphericalHarmonic::getOrder() const
{
return m_order;
}
//---------------------------------------------------------------------------
//
// 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/>.
//
//---------------------------------------------------------------------------
#ifndef WSYMMETRICSPHERICALHARMONIC_H
#define WSYMMETRICSPHERICALHARMONIC_H
#include <vector>
#include "WMath.h"
#include "WUnitSphereCoordinates.h"
#include "WValue.h"
namespace wmath
{
/**
* Class for symmetric spherical harmonics
* The index scheme of the coefficients/basis values is like in the Descoteaux paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging".
*/
class WSymmetricSphericalHarmonic
{
// TODO(all): implement test
// friend class WSymmetricSphericalHarmonicTest;
public:
/**
* Default constructor.
*/
WSymmetricSphericalHarmonic();
/**
* Constructor.
* \param SHCoefficients the initial coefficients (stored like in the mentioned Descoteaux paper).
*/
explicit WSymmetricSphericalHarmonic( const WValue<double>& SHCoefficients );
/**
* Destructor.
*/
virtual ~WSymmetricSphericalHarmonic();
/**
* Return the value on the sphere.
* \param theta angle for the position on the unit sphere
* \param phi angle for the position on the unit sphere
*/
double getValue( double theta, double phi ) const;
/**
* Return the value on the sphere.
* \param coordinates for the position on the unit sphere
*/
double getValue( const WUnitSphereCoordinates& coordinates ) const;
/**
* Returns the used coefficients (stored like in the mentioned Descoteaux paper).
*/
const wmath::WValue<double>& getCoefficients() const;
/**
* Applies the Funk-Radon-Transformation.
*/
void applyFunkRadonTransformation();
/**
* Return the order of the spherical harmonic.
*/
size_t getOrder() const;
protected:
private:
/** order of the spherical harmonic */
size_t m_order;
/** coefficients of the spherical harmonic */
WValue<double> m_SHCoefficients;
};
}
#endif // WSYMMETRICSPHERICALHARMONIC_H
//---------------------------------------------------------------------------
//
// 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/>.
//
//---------------------------------------------------------------------------
#include <cmath>
#include "WUnitSphereCoordinates.h"
using wmath::WUnitSphereCoordinates;
WUnitSphereCoordinates::WUnitSphereCoordinates() :
m_theta( 0.0 ),
m_phi( 0.0 )
{
}
WUnitSphereCoordinates::WUnitSphereCoordinates( double theta, double phi ) :
m_theta( theta ),
m_phi( phi )
{
}
WUnitSphereCoordinates::WUnitSphereCoordinates( wmath::WVector3D vector )
{
// normalize vector
vector *= vector.norm();
// calculate angles
m_theta = std::acos( vector[2] );
m_phi = std::atan2( vector[1], vector[0] );
}
WUnitSphereCoordinates::~WUnitSphereCoordinates()
{
}
double WUnitSphereCoordinates::getTheta() const
{
return m_theta;
}
double WUnitSphereCoordinates::getPhi() const
{
return m_phi;
}
void WUnitSphereCoordinates::setTheta( double theta )
{
m_theta = theta;
}
void WUnitSphereCoordinates::setPhi( double phi )
{
m_phi = phi;
}
wmath::WVector3D WUnitSphereCoordinates::getEuclidean() const
{
return wmath::WVector3D( std::sin( m_theta )*std::cos( m_phi ), std::sin( m_theta )*std::sin( m_phi ), std::cos( m_theta ) );
}
//---------------------------------------------------------------------------
//
// 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/>.
//
//---------------------------------------------------------------------------
#ifndef WUNITSPHERECOORDINATES_H
#define WUNITSPHERECOORDINATES_H
#include <vector>
#include "../../common/math/WVector3D.h"
namespace wmath
{
/**
* This class stores coordinates on the unit sphere.
*/
class WUnitSphereCoordinates
{
// TODO(all): implement test
// friend class WUnitSphereCoordinatesTest;
public:
/**
* Default constructor.
*/
WUnitSphereCoordinates();
/**
* Constructor for unit sphere angles.
* \param theta coordinate
* \param phi coordinate
*/
WUnitSphereCoordinates( double theta, double phi );
/**
* Constructor for euclidean coordinates.
* \param vector euclidean coordinates
*/
explicit WUnitSphereCoordinates( wmath::WVector3D vector );
/**
* Destructor.
*/
virtual ~WUnitSphereCoordinates();
/**
* Return the theta angle.
*/
double getTheta() const;
/**
* Return the phi angle.
*/
double getPhi() const;
/**
* Set theta angle.
* \param theta Value for theta.
*/
void setTheta( double theta );
/**
* Set phi angle.
* \param phi Value for phi.
*/
void setPhi( double phi );
/**
* Returns the stored sphere coordinates as euclidean coordinates.
*/
wmath::WVector3D getEuclidean() const;
protected:
private:
/** coordinate */
double m_theta;
/** coordinate */
double m_phi;
};
}
#endif // WUNITSPHERECOORDINATES_H
//---------------------------------------------------------------------------
//
// 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/>.
//
//---------------------------------------------------------------------------
#include <stdint.h>
#include <string>
#include <vector>
#include "../common/WAssert.h"
#include "WDataSetSingle.h"
#include "WDataSetSphericalHarmonics.h"
// prototype instance as singleton
boost::shared_ptr< WPrototyped > WDataSetSphericalHarmonics::m_prototype = boost::shared_ptr< WPrototyped >();
WDataSetSphericalHarmonics::WDataSetSphericalHarmonics( boost::shared_ptr< WValueSetBase > newValueSet,
boost::shared_ptr< WGrid > newGrid ) :
WDataSetSingle( newValueSet, newGrid )
{
m_valueSet = boost::shared_dynamic_cast< WValueSet<double> >( newValueSet );
WAssert( newValueSet, "No value set given." );
WAssert( newGrid, "No grid given." );
WAssert( m_valueSet, "No WValueSet<double given." );
}
WDataSetSphericalHarmonics::WDataSetSphericalHarmonics()
: WDataSetSingle()
{
}
WDataSetSphericalHarmonics::~WDataSetSphericalHarmonics()
{
}
boost::shared_ptr< WPrototyped > WDataSetSphericalHarmonics::getPrototype()
{
if ( !m_prototype )
{
m_prototype = boost::shared_ptr< WPrototyped >( new WDataSetSphericalHarmonics() );
}
return m_prototype;
}
wmath::WSymmetricSphericalHarmonic WDataSetSphericalHarmonics::interpolate( const wmath::WPosition& pos, bool* success ) const
{
boost::shared_ptr< WGridRegular3D > grid = boost::shared_dynamic_cast< WGridRegular3D >( m_grid );
*success = grid->encloses( pos );
if( !*success )
{
return wmath::WSymmetricSphericalHarmonic();
}
// ids of vertices for interpolation
std::vector< size_t > vertexIds = grid->getCellVertexIds( grid->getCellId( pos ) );
wmath::WPosition localPos = pos - grid->getPosition( vertexIds[0] );
double lambdaX = localPos[0] / grid->getOffsetX();
double lambdaY = localPos[1] / grid->getOffsetY();
double lambdaZ = localPos[2] / grid->getOffsetZ();
std::vector< double > h( 8 );
// lZ lY
// | /
// | 6___/_7
// |/: /|
// 4_:___5 |
// | :...|.|
// |.2 | 3
// |_____|/ ____lX
// 0 1
h[0] = ( 1 - lambdaX ) * ( 1 - lambdaY ) * ( 1 - lambdaZ );
h[1] = ( lambdaX ) * ( 1 - lambdaY ) * ( 1 - lambdaZ );
h[2] = ( 1 - lambdaX ) * ( lambdaY ) * ( 1 - lambdaZ );
h[3] = ( lambdaX ) * ( lambdaY ) * ( 1 - lambdaZ );
h[4] = ( 1 - lambdaX ) * ( 1 - lambdaY ) * ( lambdaZ );
h[5] = ( lambdaX ) * ( 1 - lambdaY ) * ( lambdaZ );
h[6] = ( 1 - lambdaX ) * ( lambdaY ) * ( lambdaZ );
h[7] = ( lambdaX ) * ( lambdaY ) * ( lambdaZ );
// take
wmath::WValue< double > interpolatedCoefficients( m_valueSet->dimension() );
for( size_t i = 0; i < 8; ++i )
{
interpolatedCoefficients += h[i] * m_valueSet->getWValue( vertexIds[i] );
}
*success = true;
return wmath::WSymmetricSphericalHarmonic( interpolatedCoefficients );
}
wmath::WSymmetricSphericalHarmonic WDataSetSphericalHarmonics::getSphericalHarmonicAt( size_t index ) const
{
if ( index < m_valueSet->size() ) return wmath::WSymmetricSphericalHarmonic( m_valueSet->getWValue( index ) );
return wmath::WSymmetricSphericalHarmonic();
}
const std::string WDataSetSphericalHarmonics::getName() const
{
return "WDataSetSphericalHarmonics";
}
const std::string WDataSetSphericalHarmonics::getDescription() const
{
return "Contains factors for spherical harmonics.";
}
//---------------------------------------------------------------------------
//
// 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/>.
//
//---------------------------------------------------------------------------
#ifndef WDATASETSPHERICALHARMONICS_H
#define WDATASETSPHERICALHARMONICS_H
#include <string>
#include <vector>
#include "../common/math/WSymmetricSphericalHarmonic.h"
#include "WValueSet.h"
#include "WDataSetSingle.h"