diff --git a/src/rawtoaces_core/mathOps.h b/src/rawtoaces_core/mathOps.h index eb9b8eab..3f96764f 100644 --- a/src/rawtoaces_core/mathOps.h +++ b/src/rawtoaces_core/mathOps.h @@ -20,17 +20,6 @@ namespace core { // Non-class functions -inline double invertD( double val ) -{ - assert( fabs( val - 0.0 ) >= DBL_EPSILON ); - - return 1.0 / val; -}; - -template T clip( T val, T target ) -{ - return std::min( val, target ); -}; template int isSquare( const vector> &vm ) { @@ -126,16 +115,6 @@ template vector invertV( const vector &vMtx ) return result; }; -template vector> diagVM( const vector &vct ) -{ - assert( vct.size() != 0 ); - vector> vctdiag( vct.size(), vector( vct.size(), T( 0.0 ) ) ); - - FORI( vct.size() ) vctdiag[i][i] = vct[i]; - - return vctdiag; -}; - template vector diagV( const vector &vct ) { assert( vct.size() != 0 ); @@ -205,55 +184,6 @@ template void scaleVector( vector &vct, const T scale ) return; }; -template void scale_vector_max( vector &vector ) -{ - Eigen::Matrix column_vector; - column_vector.resize( vector.size(), 1 ); - - FORI( vector.size() ) - { - column_vector( i, 0 ) = vector[i]; - } - column_vector *= ( 1.0 / column_vector.maxCoeff() ); - - FORI( vector.size() ) - { - vector[i] = column_vector( i, 0 ); - } - - return; -}; - -template void scale_vector_min( vector &vector ) -{ - Eigen::Matrix column_vector; - column_vector.resize( vector.size(), 1 ); - - FORI( vector.size() ) - { - column_vector( i, 0 ) = vector[i]; - } - column_vector *= ( 1.0 / column_vector.minCoeff() ); - - FORI( vector.size() ) - { - vector[i] = column_vector( i, 0 ); - } - - return; -}; - -template void scaleVectorD( vector &vct ) -{ - Eigen::Matrix v; - v.resize( vct.size(), 1 ); - - FORI( v.rows() ) v( i, 0 ) = vct[i]; - FORI( v.rows() ) vct[i] = v.maxCoeff() / vct[i]; - - return; -}; - template vector mulVectorElement( const vector &vct1, const vector &vct2 ) { @@ -275,21 +205,6 @@ vector mulVectorElement( const vector &vct1, const vector &vct2 ) return vct3; }; -template -vector divVectorElement( const vector &vct1, const vector &vct2 ) -{ - assert( vct1.size() == vct2.size() ); - - vector vct2D( vct2.size(), T( 1.0 ) ); - FORI( vct2.size() ) - { - assert( vct2[i] != T( 0.0 ) ); - vct2D[i] = T( 1.0 ) / vct2[i]; - } - - return mulVectorElement( vct1, vct2D ); -}; - template vector mulVector( vector vct1, vector vct2, int k = 3 ) { @@ -368,73 +283,6 @@ vector mulVector( const vector &vct1, const vector> &vct2 ) return mulVector( vct2, vct1 ); }; -template -T *mulVectorArray( - T *data, - const uint32_t total, - const uint8_t dim, - const vector> &vct ) -{ - assert( vct.size() == dim && isSquare( vct ) ); - - /** - // new implementation based on Eigen::Eigen::Matrix (Slow...) - - Eigen::Matrix MI, mvct; - MI.resize(total/dim, dim); - mvct.resize(dim, dim); - FORIJ(MI.rows(), MI.cols()) MI(i,j) = data[i*dim+j]; - FORIJ(dim, dim) mvct(i,j) = static_cast(vct[i][j]); - - Eigen::Matrix MR(MI * (mvct.transpose())); - FORI(total) data[i] = MR(i); - */ - - if ( dim == 3 || dim == 4 ) - { - for ( uint32_t i = 0; i < total; i += dim ) - { - T temp[4]; - - for ( uint8_t j = 0; j < dim; j++ ) - { - temp[j] = 0; - - for ( uint8_t k = 0; k < dim; k++ ) - temp[j] += vct[j][k] * data[i + k]; - } - - for ( uint8_t j = 0; j < dim; j++ ) - data[i + j] = temp[j]; - } - } - - return data; -}; - -template -vector> -solveVM( const vector> &vct1, const vector> &vct2 ) -{ - - Eigen::Matrix m1, m2, m3; - m1.resize( vct1.size(), vct1[0].size() ); - m2.resize( vct2.size(), vct2[0].size() ); - - FORIJ( vct1.size(), vct1[0].size() ) - m1( i, j ) = vct1[i][j]; - FORIJ( vct2.size(), vct2[0].size() ) - m2( i, j ) = vct2[i][j]; - - // colPivHouseholderQr() - m3 = m1.jacobiSvd( Eigen::ComputeThinU | Eigen::ComputeThinV ).solve( m2 ); - - vector> vct3( m3.rows(), vector( m3.cols() ) ); - FORIJ( m3.rows(), m3.cols() ) vct3[i][j] = m3( i, j ); - - return vct3; -}; - /// Calculate the Sum of Squared Errors (SSE) between two vectors. /// The SSE measures how well the calculated values (tcp) match the reference values (src). /// Formula: Σ((tcp[i] / src[i] - 1)²) diff --git a/tests/testMath.cpp b/tests/testMath.cpp index fd4b1be4..eb027711 100644 --- a/tests/testMath.cpp +++ b/tests/testMath.cpp @@ -13,30 +13,6 @@ using namespace rta::core; -void test_InvertD() -{ - double a = 1.0; - OIIO_CHECK_EQUAL_THRESH( invertD( a ), 1.0, 1e-9 ); - - double b = 1000.0; - OIIO_CHECK_EQUAL_THRESH( invertD( b ), 0.001, 1e-9 ); - - double c = 1000000.0; - OIIO_CHECK_EQUAL_THRESH( invertD( c ), 0.000001, 1e-9 ); -}; - -void test_Clip() -{ - double a = 254.9; - OIIO_CHECK_EQUAL_THRESH( clip( a, 255.0 ), a, 1e-5 ); - - double b = 255.1; - OIIO_CHECK_EQUAL_THRESH( clip( b, 255.0 ), 255.0, 1e-5 ); - - double c = 63355.0; - OIIO_CHECK_EQUAL_THRESH( clip( c, 63355.0 ), c, 1e-5 ); -}; - void test_IsSquare() { vector> a; @@ -126,24 +102,6 @@ void test_InvertV() FORI( 9 ) OIIO_CHECK_EQUAL_THRESH( V_Inverse[i], MV_Inverse[i], 1e-5 ); }; -void test_DiagVM() -{ - double M[3][3] = { { 1.0, 0.0, 0.0 }, - { 0.0, 2.0, 0.0 }, - { 0.0, 0.0, 3.0 } }; - - double vd[3] = { 1.0, 2.0, 3.0 }; - vector MV( vd, vd + 3 ); - vector> MVD = diagVM( MV ); - - FORI( 3 ) - { - OIIO_CHECK_EQUAL_THRESH( MVD[i][0], M[i][0], 1e-5 ); - OIIO_CHECK_EQUAL_THRESH( MVD[i][1], M[i][1], 1e-5 ); - OIIO_CHECK_EQUAL_THRESH( MVD[i][2], M[i][2], 1e-5 ); - } -}; - void test_DiagV() { double v[3] = { 1.0, 2.0, 3.0 }; @@ -200,41 +158,6 @@ void test_SumVector() OIIO_CHECK_EQUAL_THRESH( sum, 55.0000, 1e-5 ); }; -void test_ScaleVectorMax() -{ - double M[10] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 }; - double M_Scaled[10] = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 }; - vector MV( M, M + 10 ); - - scale_vector_max( MV ); - FORI( MV.size() ) - OIIO_CHECK_EQUAL_THRESH( M_Scaled[i], MV[i], 1e-5 ); -}; - -void test_ScaleVectorMin() -{ - double M[10] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 }; - vector MV( M, M + 10 ); - - scale_vector_min( MV ); - FORI( MV.size() ) - OIIO_CHECK_EQUAL_THRESH( M[i], MV[i], 1e-5 ); -}; - -void test_scaleVectorD() -{ - double M[10] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 }; - double M_Scaled[10] = { 10.0000000000, 5.0000000000, 3.3333333333, - 2.5000000000, 2.0000000000, 1.6666666667, - 1.4285714286, 1.2500000000, 1.1111111111, - 1.0000000000 }; - vector MV( M, M + 10 ); - - scaleVectorD( MV ); - FORI( MV.size() ) - OIIO_CHECK_EQUAL_THRESH( MV[i], M_Scaled[i], 1e-5 ); -}; - void test_MulVectorElement() { double M1[10] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 }; @@ -249,19 +172,6 @@ void test_MulVectorElement() OIIO_CHECK_EQUAL_THRESH( MV3[i], 10.0000000000, 1e-5 ); }; -void test_DivVectorElement() -{ - double M1[10] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 }; - double M2[10] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 }; - - vector MV1( M1, M1 + 10 ); - vector MV2( M2, M2 + 10 ); - - vector MV3 = divVectorElement( MV1, MV2 ); - FORI( MV3.size() ) - OIIO_CHECK_EQUAL_THRESH( MV3[i], 1.0000000000, 1e-5 ); -}; - void test_MulVector1() { double M1[3][3] = { { 1.0, 0.0, 0.0 }, @@ -312,61 +222,6 @@ void test_MulVector2() OIIO_CHECK_EQUAL_THRESH( MV3[i], M2[i], 1e-5 ); }; -void test_MulVectorArray() -{ - double data[9] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0 }; - double M[3][3] = { { 1.0000000000, 0.1000000000, 0.0100000000 }, - { 0.1000000000, 2.0000000000, 0.0100000000 }, - { 0.1000000000, 0.0100000000, 3.0000000000 } - - }; - - double data_test[9] = { 1.2300000000, 4.13000000000, 9.12000000000, - 4.5600000000, 10.4600000000, 18.4500000000, - 7.8900000000, 16.7900000000, 27.7800000000 }; - - vector> MV( 3, vector( 3 ) ); - FORIJ( 3, 3 ) - MV[i][j] = M[i][j]; - - mulVectorArray( data, 9, 3, MV ); - FORI( 9 ) - OIIO_CHECK_EQUAL_THRESH( data[i], data_test[i], 1e-5 ); -}; - -void test_SolveVM() -{ - double M1[3][3] = { { 1.0000000000, 0.0000000000, 0.0000000000 }, - { 0.0000000000, 2.0000000000, 0.0000000000 }, - { 0.0000000000, 0.0000000000, 3.0000000000 } - - }; - double M2[3][3] = { { 1.0000000000, 0.0000000000, 0.0000000000 }, - { 0.0000000000, 1.0000000000, 0.0000000000 }, - { 0.0000000000, 0.0000000000, 1.0000000000 } - - }; - - double M3_test[3][3] = { { 1.0000000000, 0.0000000000, 0.0000000000 }, - { 0.0000000000, 0.5000000000, 0.0000000000 }, - { 0.0000000000, 0.0000000000, 0.3333333333 } - - }; - - vector> MV1( 3, vector( 3 ) ); - vector> MV2( 3, vector( 3 ) ); - - FORIJ( 3, 3 ) - { - MV1[i][j] = M1[i][j]; - MV2[i][j] = M2[i][j]; - } - - vector> MV3 = solveVM( MV1, MV2 ); - FORIJ( 3, 3 ) - OIIO_CHECK_EQUAL_THRESH( MV3[i][j], M3_test[i][j], 1e-5 ); -}; - void test_FindIndexInterp1() { int M[100]; @@ -909,27 +764,18 @@ void test_GetCalcXYZt() int main( int, char ** ) { - test_InvertD(); - test_Clip(); test_IsSquare(); test_AddVectors(); test_SubVectors(); test_Cross2(); test_InvertVM(); test_InvertV(); - test_DiagVM(); test_DiagV(); test_TransposeVec(); test_SumVector(); - test_ScaleVectorMax(); - test_ScaleVectorMin(); - test_scaleVectorD(); test_MulVectorElement(); - test_DivVectorElement(); test_MulVector1(); test_MulVector2(); - test_MulVectorArray(); - test_SolveVM(); test_FindIndexInterp1(); test_Interp1DLinear(); testIDT_XytoXYZ();