From d1641ad5bb3587c30c473e78a6ae060a19d7a60e Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Fri, 27 Mar 2015 16:08:23 -0500 Subject: [PATCH] Generic gaussian model estimator added --- appli/examples/CMakeLists.txt | 15 +++ .../example_ImageGaussianModelEstimator.cxx | 85 +++++++++++++ .../IterativeGaussianModelEstimator.h | 84 +++++++++++++ .../IterativeGaussianModelEstimator.hxx | 114 ++++++++++++++++++ 4 files changed, 298 insertions(+) create mode 100644 appli/examples/example_ImageGaussianModelEstimator.cxx create mode 100644 lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.h create mode 100644 lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.hxx diff --git a/appli/examples/CMakeLists.txt b/appli/examples/CMakeLists.txt index 3a79578..8ef86f1 100644 --- a/appli/examples/CMakeLists.txt +++ b/appli/examples/CMakeLists.txt @@ -28,4 +28,19 @@ FOREACH(prog ${EXAMPLES_PROGRAMS}) ) ENDFOREACH(prog) +SET( + NOPLUGINS_EXAMPLES_PROGRAMS + example_ImageGaussianModelEstimator + ) +FOREACH(prog ${NOPLUGINS_EXAMPLES_PROGRAMS}) + ADD_EXECUTABLE( + ${prog} + ${prog}.cxx + ) + TARGET_LINK_LIBRARIES( + ${prog} + ${ITK_LIBRARIES} + ) +ENDFOREACH(prog) + ## eof - $RCSfile$ diff --git a/appli/examples/example_ImageGaussianModelEstimator.cxx b/appli/examples/example_ImageGaussianModelEstimator.cxx new file mode 100644 index 0000000..0187f2a --- /dev/null +++ b/appli/examples/example_ImageGaussianModelEstimator.cxx @@ -0,0 +1,85 @@ +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +// ------------------------------------------------------------------------- +const unsigned int Dim = 2; +typedef double TScalar; +typedef unsigned char TChannel; +typedef itk::RGBPixel< TChannel > TPixel; +typedef itk::Image< TPixel, Dim > TImage; + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + if( argc < 2 ) + { + std::cerr + << "Usage: " << argv[ 0 ] + << " input_image" << std::endl; + return( 1 ); + + } // fi + std::string input_image_file = argv[ 1 ]; + + // Read image + itk::ImageFileReader< TImage >::Pointer reader = + itk::ImageFileReader< TImage >::New( ); + reader->SetFileName( input_image_file ); + try + { + reader->Update( ); + } + catch( itk::ExceptionObject& err ) + { + std::cerr << "Error caught: " << err << std::endl; + return( 1 ); + + } // yrt + TImage::ConstPointer input_image = reader->GetOutput( ); + + // Compute gaussian model + typedef cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< TScalar, 3 > TEstimator; + + TEstimator::Pointer estimator = TEstimator::New( ); + estimator->Clear( ); + + itk::ImageRegionConstIterator< TImage > it( + input_image, + input_image->GetRequestedRegion( ) + ); + for( it.GoToBegin( ); !it.IsAtEnd( ); ++it ) + { + estimator->AddSample( + TScalar( it.Get( ).GetRed( ) ), + TScalar( it.Get( ).GetGreen( ) ), + TScalar( it.Get( ).GetBlue( ) ) + ); + + } // rof + + itk::Vector< TScalar, 3 > mean; + itk::Matrix< TScalar, 3, 3 > cova; + + estimator->GetModel( mean, cova ); + std::cout + << "Covariance = " << std::endl << cova << std::endl; + std::cout + << "Inverse covariance = " + << std::endl << cova.GetInverse( ) << std::endl; + std::cout + << "Mean = " << std::endl << mean << std::endl; + + return( 0 ); +} + +// eof - $RCSfile$ diff --git a/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.h b/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.h new file mode 100644 index 0000000..3e604d2 --- /dev/null +++ b/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.h @@ -0,0 +1,84 @@ +// ------------------------------------------------------------------------- +// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co) +// ------------------------------------------------------------------------- + +#ifndef __CPPLUGINS__EXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__H__ +#define __CPPLUGINS__EXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__H__ + +#include +#include +#include + +namespace cpPlugins +{ + namespace Extensions + { + namespace Algorithms + { + /** + */ + template< class S, unsigned int D > + class IterativeGaussianModelEstimator + : public itk::Object + { + public: + typedef IterativeGaussianModelEstimator Self; + typedef itk::Object Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef S TScalar; + itkStaticConstMacro( Dimension, unsigned int, D ); + + // Begin concept checking +#ifdef ITK_USE_CONCEPT_CHECKING + itkConceptMacro( + ScalarTypeHasFloatResolution, + ( itk::Concept::IsFloatingPoint< S > ) + ); +#endif + // End concept checking + + typedef vnl_matrix< S > TMatrix; + + public: + itkNewMacro( Self ); + itkTypeMacro( IterativeGaussianModelEstimator, itkObject ); + + public: + template< class V, class M > + void GetModel( V& m, M& E ) const; + + void Clear( ); + + template< class V > + void AddSample( const V& sample ); + + void AddSample( const S& s_x, const S& s_y, ... ); + + protected: + IterativeGaussianModelEstimator( ); + virtual ~IterativeGaussianModelEstimator( ); + + private: + // Purposely not implemented + IterativeGaussianModelEstimator( const Self& other ); + void operator=( const Self& other ); + + protected: + unsigned long m_N; + TMatrix m_S1; + TMatrix m_S2; + }; + + } // ecapseman + + } // ecapseman + +} // ecapseman + +#include + +#endif // __CPPLUGINS__EXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__H__ + +// eof - $RCSfile$ diff --git a/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.hxx b/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.hxx new file mode 100644 index 0000000..fa8d8bb --- /dev/null +++ b/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.hxx @@ -0,0 +1,114 @@ +// ------------------------------------------------------------------------- +// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co) +// ------------------------------------------------------------------------- + +#ifndef __CPPLUGINS__EXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__HXX__ +#define __CPPLUGINS__EXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__HXX__ + +#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 +{ + TMatrix Es = ( this->m_S1 * this->m_S1.transpose( ) ); + + for( unsigned int i = 0; i < D; ++i ) + { + m[ i ] = this->m_S1[ i ][ 0 ]; + if( this->m_N > 0 ) + m[ i ] /= double( this->m_N ); + + for( unsigned int j = 0; j < D; ++j ) + { + 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 ); + + } // fi + + } // rof + + } // rof +} + +// ------------------------------------------------------------------------- +template< class S, unsigned int D > +void +cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >:: +Clear( ) +{ + this->m_N = 0; + this->m_S1.set_size( D, 1 ); + this->m_S2.set_size( D, D ); + this->m_S1.fill( S( 0 ) ); + this->m_S2.fill( S( 0 ) ); + this->Modified( ); +} + +// ------------------------------------------------------------------------- +template< class S, unsigned int D > +template< class V > +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++; + this->Modified( ); +} + +// ------------------------------------------------------------------------- +template< class S, unsigned int D > +void +cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >:: +AddSample( const S& s_x, const S& s_y, ... ) +{ + 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 ); + this->AddSample( sample ); +} + +// ------------------------------------------------------------------------- +template< class S, unsigned int D > +cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >:: +IterativeGaussianModelEstimator( ) + : Superclass( ) +{ + this->Clear( ); +} + +// ------------------------------------------------------------------------- +template< class S, unsigned int D > +cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >:: +~IterativeGaussianModelEstimator( ) +{ +} + +#endif // __CPPLUGINS__EXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__HXX__ + +// eof - $RCSfile$ -- 2.45.1