#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$