Generic gaussian model estimator added

diff --git a/appli/examples/CMakeLists.txt b/appli/examples/CMakeLists.txt
index 3a795782084f8aa01e75028c0e9087ea766f22fc..8ef86f17b9fd9488bfeb0e50ff6816f222b7e2ce 100644 (file)
--- 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 (file)
index 0000000..0187f2a
--- /dev/null
+++ b/appli/examples/example_ImageGaussianModelEstimator.cxx
@@ -0,0 +1,85 @@
+#include <cstdlib>
+#include <iostream>
+#include <string>
+
+#include <itkImage.h>
+#include <itkImageFileReader.h>
+#include <itkImageRegionConstIterator.h>
+#include <itkMatrix.h>
+#include <itkRGBPixel.h>
+#include <itkVector.h>
+
+#include <cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.h>
+
+// -------------------------------------------------------------------------
+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 (file)
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 <itkConceptChecking.h>
+#include <itkObject.h>
+#include <vnl/vnl_matrix.h>
+
+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 <cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.hxx>
+
+#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 (file)
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 <cstdarg>
+
+// -------------------------------------------------------------------------
+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$