From d896eba488c793820e129076cab4ec01d64818a8 Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Mon, 6 Apr 2015 21:10:34 -0500 Subject: [PATCH] Filters improved --- .../IterativeGaussianModelEstimator.h | 33 ++- .../IterativeGaussianModelEstimator.hxx | 195 +++++++++++++++--- .../Algorithms/RGBToYPbPrFunction.h | 6 +- 3 files changed, 194 insertions(+), 40 deletions(-) diff --git a/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.h b/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.h index 148fee9..94e5495 100644 --- a/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.h +++ b/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.h @@ -5,6 +5,7 @@ #ifndef __CPPLUGINS__EXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__H__ #define __CPPLUGINS__EXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__H__ +#include #include #include #include @@ -39,15 +40,26 @@ namespace cpPlugins #endif // End concept checking - typedef vnl_matrix< S > TMatrix; + typedef vnl_matrix< S > TMatrix; + typedef std::vector< TMatrix > TMatrices; public: itkNewMacro( Self ); itkTypeMacro( IterativeGaussianModelEstimator, itkObject ); + itkGetConstMacro( MinimumProbability, S ); + itkGetConstMacro( MaximumProbability, S ); + public: - const unsigned long& GetNumberOfSamples( ) const - { return( this->m_N ); } + unsigned long GetNumberOfSamples( ) const + { return( this->m_Samples.size( ) ); } + + void UpdateModel( ); + + template< class V > + S Probability( const V& sample ) const; + + S Probability( const S& s_x, const S& s_y, ... ) const; template< class V, class M > void GetModel( V& m, M& E ) const; @@ -69,9 +81,18 @@ namespace cpPlugins void operator=( const Self& other ); protected: - unsigned long m_N; - TMatrix m_S1; - TMatrix m_S2; + TMatrix m_S1; + TMatrix m_S2; + TMatrices m_Samples; + + bool m_Updated; + TMatrix m_Cov; + TMatrix m_Mean; + TMatrix m_InvCov; + S m_DensityCoeff; + + S m_MinimumProbability; + S m_MaximumProbability; }; } // ecapseman diff --git a/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.hxx b/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.hxx index fa8d8bb..c8c652f 100644 --- a/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.hxx +++ b/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.hxx @@ -5,36 +5,154 @@ #ifndef __CPPLUGINS__EXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__HXX__ #define __CPPLUGINS__EXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__HXX__ +#include #include +#include +#include // ------------------------------------------------------------------------- template< class S, unsigned int D > -template< class V, class M > void cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >:: -GetModel( V& m, M& E ) const +UpdateModel( ) { - TMatrix Es = ( this->m_S1 * this->m_S1.transpose( ) ); - - for( unsigned int i = 0; i < D; ++i ) + if( this->m_Updated ) + return; + + this->m_Cov.set_size( D, D ); + this->m_Mean.set_size( D, 1 ); + + // Compute covariance matrix and mean vector + unsigned long N = this->m_Samples.size( ); + this->m_Mean = this->m_S1; + if( N > 0 ) + this->m_Mean /= S( N ); + if( N > 1 ) { - m[ i ] = this->m_S1[ i ][ 0 ]; - if( this->m_N > 0 ) - m[ i ] /= double( this->m_N ); + this->m_Cov = this->m_S2; + this->m_Cov -= ( this->m_S1 * this->m_S1.transpose( ) ) / S( N ); + this->m_Cov /= S( N - 1 ); + } + else + this->m_Cov.fill( S( 0 ) ); + + /* TODO + TMatrix Es = ( this->m_S1 * this->m_S1.transpose( ) ); + for( unsigned int i = 0; i < D; ++i ) + { + this->m_Mean[ i ][ 0 ] = this->m_S1[ i ][ 0 ]; + if( this->m_N > 0 ) + this->m_Mean[ i ][ 0 ] /= S( this->m_N ); + + for( unsigned int j = 0; j < D; ++j ) + { + this->m_Cov[ i ][ j ] = S( 0 ); + if( this->m_N > 0 ) + this->m_Cov[ i ][ j ] -= Es[ i ][ j ] / S( this->m_N ); + if( this->m_N > 1 ) + { + this->m_Cov[ i ][ j ] += this->m_S2[ i ][ j ]; + this->m_Cov[ i ][ j ] /= S( this->m_N - 1 ); + + } // fi + + } // rof + + } // rof + */ + + // Compute inverse and determinant + S det = vnl_determinant( this->m_Cov ); + if( !( det > S( 0 ) ) ) + { + this->m_InvCov.set_size( D, D ); + this->m_InvCov.fill( S( 0 ) ); + this->m_DensityCoeff = S( 0 ); + } + else + { + this->m_InvCov = vnl_inverse( this->m_Cov ); + double _2piD = std::pow( double( 2 ) * double( vnl_math::pi ), D ); + this->m_DensityCoeff = S( 1 ) / S( std::sqrt( _2piD * double( det ) ) ); - for( unsigned int j = 0; j < D; ++j ) + } // fi + + // Compute minimum and maximum probabilities from input samples + static S sample[ D ]; + for( unsigned long i = 0; i < this->m_Samples.size( ); ++i ) + { + for( unsigned int d = 0; d < D; ++d ) + sample[ d ] = this->m_Samples[ i ][ d ][ 0 ]; + S p = this->Probability( sample ); + if( i == 0 ) { - E[ i ][ j ] = double( 0 ); - if( this->m_N > 0 ) - E[ i ][ j ] -= double( Es[ i ][ j ] / S( this->m_N ) ); - if( this->m_N > 1 ) - { - E[ i ][ j ] += double( this->m_S2[ i ][ j ] ); - E[ i ][ j ] /= double( this->m_N - 1 ); + this->m_MinimumProbability = p; + this->m_MaximumProbability = p; + } + else + { + if( p < this->m_MinimumProbability ) this->m_MinimumProbability = p; + if( this->m_MaximumProbability < p ) this->m_MaximumProbability = p; + + } // fi + + } // rof + this->m_Updated = true; + + /* DEBUG: + std::cout << "--------------------" << std::endl; + std::cout << this->m_Cov << std::endl; + std::cout << "--------------------" << std::endl; + std::cout << this->m_InvCov << std::endl; + std::cout << "--------------------" << std::endl; + std::cout << this->m_Mean << std::endl; + */ +} - } // fi +// ------------------------------------------------------------------------- +template< class S, unsigned int D > +template< class V > +S cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >:: +Probability( const V& sample ) const +{ + TMatrix c( D, 1 ); + for( unsigned int d = 0; d < D; ++d ) + c[ d ][ 0 ] = S( sample[ d ] ) - this->m_Mean[ d ][ 0 ]; + double v = double( ( c.transpose( ) * ( this->m_InvCov * c ) )[ 0 ][ 0 ] ); + v /= double( 2 ); - } // rof + return( S( std::exp( -v ) ) ); +} + +// ------------------------------------------------------------------------- +template< class S, unsigned int D > +S cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >:: +Probability( const S& s_x, const S& s_y, ... ) const +{ + static S sample[ D ]; + + std::va_list args_lst; + va_start( args_lst, s_y ); + sample[ 0 ] = s_x; + sample[ 1 ] = s_y; + for( unsigned int d = 2; d < D; ++d ) + sample[ d ] = S( va_arg( args_lst, double ) ); + va_end( args_lst ); + return( this->Probability( sample ) ); +} + +// ------------------------------------------------------------------------- +template< class S, unsigned int D > +template< class V, class M > +void +cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >:: +GetModel( V& m, M& E ) const +{ + for( unsigned int i = 0; i < D; ++i ) + { + m[ i ] = double( this->m_Mean[ i ][ 0 ] ); + for( unsigned int j = 0; j < D; ++j ) + E[ i ][ j ] = double( this->m_Cov[ i ][ j ] ); } // rof } @@ -45,7 +163,8 @@ void cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >:: Clear( ) { - this->m_N = 0; + this->m_Updated = false; + this->m_Samples.clear( ); this->m_S1.set_size( D, 1 ); this->m_S2.set_size( D, D ); this->m_S1.fill( S( 0 ) ); @@ -60,18 +179,32 @@ void cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >:: AddSample( const V& sample ) { - for( unsigned int i = 0; i < D; ++i ) - { - this->m_S1[ i ][ 0 ] += S( sample[ i ] ); - for( unsigned int j = i; j < D; ++j ) - { - this->m_S2[ i ][ j ] += S( sample[ i ] ) * S( sample[ j ] ); - this->m_S2[ j ][ i ] = this->m_S2[ i ][ j ]; - - } // rof - - } // rof - this->m_N++; + TMatrix s( D, 1 ); + for( unsigned int d = 0; d < D; ++d ) + s[ d ][ 0 ] = S( sample[ d ] ); + + this->m_Samples.push_back( s ); + this->m_S1 += s; + this->m_S2 += s * s.transpose( ); + + /* DEBUG: + std::cout + << sample[ 0 ] << " " << sample[ 1 ] << " " << sample[ 2 ] << std::endl; + */ + + /* TODO + { + this->m_S1[ i ][ 0 ] += S( sample[ i ] ); + for( unsigned int j = i; j < D; ++j ) + { + this->m_S2[ i ][ j ] += S( sample[ i ] ) * S( sample[ j ] ); + this->m_S2[ j ][ i ] = this->m_S2[ i ][ j ]; + + } // rof + + } // rof + */ + this->m_Updated = false; this->Modified( ); } diff --git a/lib/cpPlugins/Extensions/Algorithms/RGBToYPbPrFunction.h b/lib/cpPlugins/Extensions/Algorithms/RGBToYPbPrFunction.h index def48e8..9563a63 100644 --- a/lib/cpPlugins/Extensions/Algorithms/RGBToYPbPrFunction.h +++ b/lib/cpPlugins/Extensions/Algorithms/RGBToYPbPrFunction.h @@ -33,9 +33,9 @@ namespace cpPlugins { static const O M[] = { - O( 0.2126 ), O( 0.7152 ), O( 0.0722 ), - O( -0.2126 ), O( -0.7152 ), O( 0.9278 ), - O( 0.7874 ), O( -0.7152 ), O( -0.0722 ) + O( 0.299 ), O( 0.587 ), O( 0.114 ), + O( -0.169 ), O( -0.331 ), O( 0.500 ), + O( 0.500 ), O( -0.419 ), O( -0.081 ) }; static const vnl_matrix< O > vM( M, 3, 3 ); static const itk::Matrix< O, 3, 3 > iM( vM ); -- 2.47.1