diff --git a/CMakeLists.txt b/CMakeLists.txt index 6ea22187e6224a2598ab5b5a8db26e6e0fd54537..9d686e871311e6b9c56962417683fe46871c3f81 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -214,6 +214,9 @@ endif( ) if( ITA_VISTA_BUILD_STATIC ) add_definitions( -DVISTABASE_STATIC -DVISTAMATH_STATIC -DVISTAASPECTS_STATIC -DVISTATOOLS_STATIC -DVISTAINTERPROCCOMM_STATIC ) + if( WIN32 ) + list( APPEND VISTA_USE_PACKAGE_LIBRARIES Ws2_32 ) + endif( ) endif( ) if( NOT WIN32 ) diff --git a/include/ITANumericUtils.h b/include/ITANumericUtils.h index eb180bb8c8769f2f2b70208b0e7baa607bf15108..1f7ae2eae685fc837f6f2fed1d8f5c94a592e49a 100644 --- a/include/ITANumericUtils.h +++ b/include/ITANumericUtils.h @@ -288,6 +288,9 @@ ITA_BASE_API float anglef_proj_0_360_DEG( const float alpha ); */ ITA_BASE_API float anglef_proj_N180_180_DEG( const float alpha ); +//! Project an angle into interval of -90 and 90 degree +ITA_BASE_API float anglef_proj_N90_90_DEG( const float alpha ); + //! (Gerichtete) minimale Winkeldifferenz für zwei Winkel [°] im Intervall [0°,360°) /** * Für zwei gegebene Winkel 0°<= alpha, beta < 360° im Interval [0°,360°) berechnet @@ -324,27 +327,21 @@ ITA_BASE_API float anglef_mindiff_0_360_DEG( const float alpha, const float beta ITA_BASE_API float anglef_mindiff_abs_0_360_DEG( const float alpha, const float beta ); -// Orientierung in Yaw-Pitch-Roll Winkeln in View-Up-Vektor umrechnen (Alle Winkel in Bogenmaß) -ITA_BASE_API void convertYPR2VU( const double yaw, const double pitch, const double roll, - double& vx, double& vy, double& vz, - double& ux, double& uy, double& uz ); +//! Orientierung in Yaw-Pitch-Roll Winkeln in View-Up-Vektor umrechnen (Alle Winkel in Bogenmaß) +ITA_BASE_API void convertYPR2VU( const double yaw, const double pitch, const double roll, double& vx, double& vy, double& vz, double& ux, double& uy, double& uz ); -// Orientierung in Yaw-Pitch-Roll Winkeln in View-Up-Vektor umrechnen (Alle Winkel in Bogenmaß) -ITA_BASE_API void convertYPR2VU( const float yaw, const float pitch, const float roll, - float& vx, float& vy, float& vz, - float& ux, float& uy, float& uz ); +//! Orientierung in Yaw-Pitch-Roll Winkeln in View-Up-Vektor umrechnen (Alle Winkel in Bogenmaß) +ITA_BASE_API void convertYPR2VU( const float yaw, const float pitch, const float roll, float& vx, float& vy, float& vz, float& ux, float& uy, float& uz ); // Orientierung in View-Up-Vektor in Yaw-Pitch-Roll Winkeln umrechnen (Alle Winkel in Bogenmaß) ITA_BASE_API void convertVU2YPR( const double vx, const double vy, const double vz, const double ux, const double uy, const double uz, double& yaw, double& pitch, double& roll ); -// Orientierung in View-Up-Vektor in Yaw-Pitch-Roll Winkeln umrechnen (Alle Winkel in Bogenmaß) -ITA_BASE_API void convertVU2YPR( const float vx, const float vy, const float vz, - const float ux, const float uy, const float uz, - float& yaw, float& pitch, float& roll ); +//! Orientierung in View-Up-Vektor in Yaw-Pitch-Roll Winkeln umrechnen (Alle Winkel in Bogenmaß) +ITA_BASE_API void convertVU2YPR( const float vx, const float vy, const float vz, const float ux, const float uy, const float uz, float& yaw, float& pitch, float& roll ); -// Datenklasse für Fehlerwerte +//! Datenklasse für Fehlerwerte template< typename T > class ITA_BASE_API ErrorValues { @@ -363,7 +360,7 @@ public: {}; }; -// Fehler zwischen zwei Vektoren berechnen +//! Fehler zwischen zwei Vektoren berechnen template< typename Tc, typename Ta, typename Tb > inline ErrorValues computeErrorValues( const Ta* A, const Tb* B, const int size ) { @@ -391,13 +388,8 @@ inline ErrorValues computeErrorValues( const Ta* A, const Tb* B, const int s return v; }; -/* +-----------------------+ - | | - | Listen und Folgen | - | | - +-----------------------+ */ -// Fills a vector with numbers dest = { a+n*s < b | n in N } +//! Fills a vector with numbers dest = { a+n*s < b | n in N } template< typename T > inline void linspace( std::vector< T >& dest, T a, T b, T s = 1 ) { dest.clear(); @@ -425,38 +417,40 @@ template< typename T > inline void linspace( std::vector< T >& dest, T a, T b, T } }; -// Fills a vector with powers of two in the range a=2^i <= 2^j <= 2^k=b -// (Note a and b must be powers of two, otherwise an exception is thrown) +//! Fills a vector with powers of two in the range a=2^i <= 2^j <= 2^k=b +/** + * @note a and b must be powers of two, otherwise an exception is thrown) + */ ITA_BASE_API void pow2space( std::vector< int >& dest, const int a, const int b ); -// returns the elevation angle (0 = frontal direction) in radians from a polar angle theta (0 = above) -ITA_BASE_API double theta2elevation(const double dThetaRAD); +//! Returns the elevation angle (0 = frontal direction) in radians from a polar angle theta (0 = above) +ITA_BASE_API double theta2elevation( const double dThetaRAD ); -// returns the polar angle theta (0 = above) in radians from a given elevation angle (0 = frontal direction) -ITA_BASE_API double elevation2theta(const double dElevationRAD); +//! returns the polar angle theta (0 = above) in radians from a given elevation angle (0 = frontal direction) +ITA_BASE_API double elevation2theta( const double dElevationRAD ); // Calculates the factorial of an positive integer m ITA_BASE_API int factorial( const int m ); -//Calculates the normalizing constant for SHRealvaluedBasefunctions +//! Calculates the normalizing constant for SHRealvaluedBasefunctions ITA_BASE_API double SHNormalizeConst( const int m, const int n ); -//Calculates the Kronecker delta +//! Calculates the Kronecker delta ITA_BASE_API int SHKronecker( const int m ); -// Returns the linear index of a basefunction with degree m and order n, linear indexing starts with 0 +//! Returns the linear index of a basefunction with degree m and order n, linear indexing starts with 0 ITA_BASE_API int SHDegreeOrder2Linear( const int m, const int n ); -// Returns degree and order of a basefunctions from a linear index, linear indexing starts with 0 -ITA_BASE_API void SHLinear2DegreeOrder(const int iLinear, int &m, int &n); +//! Returns degree and order of a basefunctions from a linear index, linear indexing starts with 0 +ITA_BASE_API void SHLinear2DegreeOrder( const int iLinear, int &m, int &n ); -// Calculates the remax weightings up to a given order -ITA_BASE_API std::vector HOARemaxWeights(int iTruncationOrder); +//! Calculates the remax weightings up to a given order +ITA_BASE_API std::vector< double > HOARemaxWeights( int iTruncationOrder ); -//Calculates the realvalued Basefunctions of SH for e.g. Ambisonics -ITA_BASE_API std::vector SHRealvaluedBasefunctions( const double thetaRAD, const double azimuthRAD, const int maxOrder ); +//! Calculates the realvalued Basefunctions of SH for e.g. Ambisonics +ITA_BASE_API std::vector< double > SHRealvaluedBasefunctions( const double thetaRAD, const double azimuthRAD, const int maxOrder ); -//Calculates the associated legendre polynomials -ITA_BASE_API std::vector SHAssociatedLegendre( const int N, const double mu ); +//! Calculates the associated legendre polynomials +ITA_BASE_API std::vector< double > SHAssociatedLegendre( const int N, const double mu ); #endif // INCLUDE_WATCHER_ITA_NUMERIC_UTILS diff --git a/src/ITANumericUtils.cpp b/src/ITANumericUtils.cpp index ebba9cd5655fa2e792fe4ac2e3e4a8a1a2064125..b32fdf4f9b22ba90067a254d1243d81a73b5d966 100644 --- a/src/ITANumericUtils.cpp +++ b/src/ITANumericUtils.cpp @@ -83,12 +83,12 @@ int uprmul( const int x, const int mul ) return ( r == 0 ? x : ( x > 0 ? x + mul - r : x - r ) ); } -unsigned int lwrmulu( const unsigned int x,const unsigned int mul ) +unsigned int lwrmulu( const unsigned int x, const unsigned int mul ) { return x - ( x % mul ); } -unsigned int uprmulu( const unsigned int x, const unsigned int mul ) +unsigned int uprmulu( const unsigned int x, const unsigned int mul ) { unsigned int r = x % mul; // Teilungsrest return ( r == 0 ? x : x + mul - r ); @@ -134,7 +134,7 @@ double grad2rad( const double phi ) float correctAngle180( const float phi ) { - if( fabs( phi ) == 180.0f ) + if( fabs( phi ) == 180.0f ) return phi; const float phi_temp = fmodf( phi, 360.0f ); @@ -242,7 +242,7 @@ void csargpabs( const float in_re, const float in_im, const float dest_arg, floa /* +----------- RAD -----------+ */ float anglef_proj_0_2PI( const float alpha ) { - float alpha_temp= fmodf( alpha, 2 * PI_F ); + float alpha_temp = fmodf( alpha, 2 * PI_F ); if( alpha_temp < 0 ) alpha_temp += 2 * PI_F; return alpha_temp; @@ -254,6 +254,14 @@ float anglef_proj_NPI_PI( const float alpha ) return x; } +float anglef_proj_NPIH_PIH( const float alpha ) +{ + float alpha_temp = fmodf( alpha + PI_F / 2.0f, PI_F ) - PI_F / 2.0f; + if( alpha_temp < -PI_F/2.0f ) + alpha_temp += PI_F; + return alpha_temp; +} + float anglef_mindiff_0_2PI( const float alpha, const float beta ) { float gamma = anglef_proj_0_2PI( beta ) - anglef_proj_0_2PI( alpha ); @@ -278,6 +286,10 @@ float anglef_proj_N180_180_DEG( const float alpha ) { return rad2gradf( anglef_proj_NPI_PI( grad2radf( alpha ) ) ); } +float anglef_proj_N90_90_DEG( const float alpha ) { + return rad2gradf( anglef_proj_NPIH_PIH( grad2radf( alpha ) ) ); +} + float anglef_mindiff_0_360_DEG( const float alpha, const float beta ) { return rad2gradf( anglef_mindiff_0_2PI( grad2radf( alpha ), grad2radf( beta ) ) ); @@ -289,7 +301,7 @@ float anglef_mindiff_abs_0_360_DEG( const float alpha, const float beta ) } -void convertYPR2VU( const float yaw, const float pitch, const float roll, float& vx, float& vy, float& vz, float& ux, float& uy, float& uz ) +void convertYPR2VU( const float yaw, const float pitch, const float roll, float& vx, float& vy, float& vz, float& ux, float& uy, float& uz ) { double vx_, vy_, vz_, ux_, uy_, uz_; convertYPR2VU( double( yaw ), double( pitch ), double( roll ), vx_, vy_, vz_, ux_, uy_, uz_ ); @@ -301,7 +313,7 @@ void convertYPR2VU( const float yaw, const float pitch, const float roll, float& uz = float( uz_ ); } -void convertYPR2VU( const double yaw, const double pitch, const double roll, double& vx, double& vy, double& vz, double& ux, double& uy, double& uz ) +void convertYPR2VU( const double yaw, const double pitch, const double roll, double& vx, double& vy, double& vz, double& ux, double& uy, double& uz ) { /* * Yaw-pitch-roll (YPR) angles rotation order (referring to OpenGL axis) @@ -323,7 +335,7 @@ void convertYPR2VU( const double yaw, const double pitch, const double roll, dou uy = cp*cr; uz = -sy*sr + cy*sp*cr; } -void convertVU2YPR( const float vx, const float vy, const float vz, const float ux, const float uy, const float uz, float& yaw, float& pitch, float& roll ) +void convertVU2YPR( const float vx, const float vy, const float vz, const float ux, const float uy, const float uz, float& yaw, float& pitch, float& roll ) { double y, p, r; convertVU2YPR( double( vx ), double( vy ), double( vz ), double( ux ), double( uy ), double( uz ), y, p, r ); @@ -332,7 +344,7 @@ void convertVU2YPR( const float vx, const float vy, const float vz, const float roll = float( r ); } -void convertVU2YPR( const double vx, const double vy, const double vz, const double ux, const double uy, const double uz, double& yaw, double& pitch, double& roll ) +void convertVU2YPR( const double vx, const double vy, const double vz, const double ux, const double uy, const double uz, double& yaw, double& pitch, double& roll ) { const float EPSILON = 0.0001F; @@ -407,7 +419,7 @@ void convertVU2YPR( const double vx, const double vy, const double vz, const dou // Hint: cos(pitch) = cos( arcsin(vy) ) = sqrt(1-vy^2) double cp = sqrt( 1 - vy*vy ); roll = ( zy <= 0 ? acos( uy / cp ) : -acos( uy / cp ) ); - + /* debug double v_norm = vx*vx + vy*vy + vz*vz; double u_norm = ux*ux + uy*uy + uz*uz; @@ -434,16 +446,32 @@ void pow2space( std::vector& dest, const int a, const int b ) } } +double theta2elevation( const double dThetaRAD ) +{ + if( dThetaRAD<0 || dThetaRAD>ITAConstants::PI_D ) + ITA_EXCEPT1( INVALID_PARAMETER, "thetaRAD is only valid between 0 and PI" ); + + return ITAConstants::PI_D / 2.0f - dThetaRAD; +} + +double elevation2theta( const double dElevationRAD ) +{ + if( abs( dElevationRAD ) < ITAConstants::PI_D / 2.0f ) + ITA_EXCEPT1( INVALID_PARAMETER, "elevationRAD is only valid between -PI/2 and PI/2" ); + + return ITAConstants::PI_D / 2.0f - dElevationRAD; +} + int factorial( const int m ) { - if( m<0 ) + if( m < 0 ) ITA_EXCEPT1( INVALID_PARAMETER, "Factorial can only be calculated for a positive integer" ); return ( m <= 1 ? 1 : m * factorial( m - 1 ) ); } double SHNormalizeConst( const int m, const int n ) { - if( m>n ) + if( m > n ) ITA_EXCEPT1( INVALID_PARAMETER, "Abs(degree) cannot be greater than order" ); double dRes = 0; @@ -452,11 +480,11 @@ double SHNormalizeConst( const int m, const int n ) else dRes = 1; - dRes *= sqrt( ( 2.0f * n + 1 ) * ( 2.0f - SHKronecker( m ) ) * factorial( n - m ) / ( 4.0f * PI_F * factorial( n + m ) ) ); + dRes *= sqrt( ( 2.0f * n + 1 ) * ( 2 - SHKronecker( m ) ) * double( factorial( n - m ) ) / ( 4.0f * PI_F * double( factorial( n + m ) ) ) ); return dRes; } -int SHKronecker( const int m ) +int SHKronecker( const int m ) { if( m == 0 ) return 1; @@ -464,12 +492,13 @@ int SHKronecker( const int m ) return 0; } - -std::vector SHRealvaluedBasefunctions( const double elevation_rad, const double azimuth_rad, const int maxOrder ) +std::vector< double > SHRealvaluedBasefunctions( const double thetaRAD, const double azimuthRAD, const int maxOrder ) { - std::vector Y; - Y.resize( ( maxOrder + 1 )*( maxOrder + 1 ) ); - Y = SHAssociatedLegendre( maxOrder, cos( elevation_rad ) ); + std::vector< double > Y; + int nmax = ( maxOrder + 1 )*( maxOrder + 1 ); + Y.resize( nmax ); + + Y = SHAssociatedLegendre( maxOrder, cos( thetaRAD ) ); for( int n = 0; n <= maxOrder; n++ ) { @@ -477,14 +506,13 @@ std::vector SHRealvaluedBasefunctions( const double elevation_rad, const for( int m = 1; m <= n; m++ ) { double Normalizing = SHNormalizeConst( m, n ); - Y[ SHDegreeOrder2Linear( m, n ) ] *= cos( m*azimuth_rad )*Normalizing; - Y[ SHDegreeOrder2Linear( -m, n ) ] *= sin( m*azimuth_rad )*Normalizing; + Y[ SHDegreeOrder2Linear( m, n ) ] *= cos( m*azimuthRAD )*Normalizing; + Y[ SHDegreeOrder2Linear( -m, n ) ] *= sin( m*azimuthRAD )*Normalizing; } } return Y; } - std::vector SHAssociatedLegendre( const int N, const double mu ) { std::vector< double > P; @@ -516,3 +544,31 @@ int SHDegreeOrder2Linear( const int m, const int n ) return ( n*n + n + m ); } + +std::vector< double > HOARemaxWeights( int iTruncationOrder ) +{ + if( iTruncationOrder > 20 ) + ITA_EXCEPT1( INVALID_PARAMETER, "ReMAx decoding only implemented up to Truncatrion Order 20" ); + + std::vector< double > vdRemax; + double dMaxRootTable[ 21 ] = { 0.993752, 0.993129, 0.992407, 0.991565, 0.990575, 0.989401, 0.987993, 0.986284, 0.984183, 0.978229, 0.978229, 0.973907, 0.96816, 0.96029, 0.949108, 0.93247, 0.90618, 0.861136, 0.774597, 0.57735, 0 }; + double dMaxRoot = dMaxRootTable[ 21 - iTruncationOrder - 1 ]; + std::vector< double > vdAssoLegendre = ( SHAssociatedLegendre( iTruncationOrder, dMaxRoot ) ); + + for( int k = 0; k <= iTruncationOrder; k++ ) + { + vdRemax.push_back( vdAssoLegendre[ SHDegreeOrder2Linear( 0, k ) ] ); + } + + return vdRemax; +} + +void SHLinear2DegreeOrder( const int iLinear, int &m, int &n ) +{ + if( iLinear < 0 ) + ITA_EXCEPT1( INVALID_PARAMETER, "The linear index has to be higher than one" ); + + std::vector nm; + n = ( int ) floor( sqrt( ( double ) iLinear ) ); + m = ( iLinear - n * n - n ); +}