From d09faa0cdbc33c7871116969b1e2dc02d5ae045e Mon Sep 17 00:00:00 2001 From: Cyrille Valladeau <cyrille.valladeau@c-s.fr> Date: Wed, 14 May 2008 07:36:06 +0000 Subject: [PATCH] Markov code. --- Code/Markov/otbMRFEnergy.h | 115 ++-- .../otbMRFEnergyGaussianClassification.h | 40 +- Code/Markov/otbMRFEnergyPotts.h | 2 +- Code/Markov/otbMRFOptimizer.h | 52 +- Code/Markov/otbMRFOptimizerICM.h | 79 +-- Code/Markov/otbMRFOptimizerMetropolis.h | 64 +- Code/Markov/otbMRFSampler.h | 165 +++-- Code/Markov/otbMRFSamplerMAP.h | 154 ++--- Code/Markov/otbMRFSamplerRandom.h | 137 ++-- Code/Markov/otbMRFSamplerRandomMAP.h | 291 ++++----- Code/Markov/otbMarkovRandomFieldFilter.h | 584 +++++++++--------- Code/Markov/otbMarkovRandomFieldFilter.txx | 304 +++++---- 12 files changed, 1003 insertions(+), 984 deletions(-) diff --git a/Code/Markov/otbMRFEnergy.h b/Code/Markov/otbMRFEnergy.h index 122d7d1993..9fd5d3a7a0 100644 --- a/Code/Markov/otbMRFEnergy.h +++ b/Code/Markov/otbMRFEnergy.h @@ -24,14 +24,21 @@ namespace otb { + /** + * \class MRFEnergy + * \brief This is the base class for energy function used in the MRF framework + * + * Derived class must reimplement the GetSingleValue() method. + */ template< class TInput1, class TInput2 > class ITK_EXPORT MRFEnergy : public itk::Object { public: - typedef MRFEnergy Self; - typedef itk::Object Superclass; - typedef itk::SmartPointer<Self> Pointer; - typedef itk::SmartPointer<const Self> ConstPointer; + typedef MRFEnergy Self; + typedef itk::Object Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + typedef TInput1 InputImageType; typedef TInput2 LabelledImageType; typedef typename InputImageType::PixelType InputImagePixelType; @@ -42,47 +49,58 @@ class ITK_EXPORT MRFEnergy : public itk::Object typedef itk::Array< double > ParametersType; - + /** Le new avait disparu...*/ itkNewMacro(Self); itkTypeMacro(MRFEnergy,itk::Object); + /** Utilisation des accesseur itk*/ + /* + void SetNumberOfParameters(unsigned int nb) + { m_NumberOfParameters = nb; this->Modified(); } unsigned int GetNumberOfParameters(void) const { return m_NumberOfParameters; } - void SetNumberOfParameters(unsigned int nb) - { m_NumberOfParameters = nb; this->Modified(); } + */ + itkSetMacro(NumberOfParameters, unsigned int); + itkGetConstMacro(NumberOfParameters, unsigned int); + // Get the parameters const ParametersType& GetParameters( void ) const { return this->m_Parameters; } + void SetParameters( const ParametersType & parameters ) - { - if( parameters.Size() != m_NumberOfParameters ) + { + if( parameters.Size() != m_NumberOfParameters ) { itkExceptionMacro(<<"Invalid number of parameters"); } m_Parameters = parameters; this->Modified(); - /* - for( unsigned int i=0; i<m_NumberOfParameters; i++ ) - { - if (m_Parameters[i] != parameters[i]) - { - m_Parameters[i] = parameters[i]; - modified = true; - } - } - if (modified) - { - this->Modified(); - } - */ + +// bool modified = false; +// for( unsigned int i=0; i<m_NumberOfParameters; i++ ) +// { +// if (m_Parameters[i] != parameters[i]) +// { +// m_Parameters[i] = parameters[i]; +// modified = true; +// } +// } +// if (modified) +// { +// this->Modified(); +// } } - virtual double GetSingleValue(const InputImagePixelType & value1, const LabelledImagePixelType & value2){return 0.;}; + /** Ne peut allouer un objet de type MRF::Energy parce que les fonctions virtuelles suivantes sont abstraites*/ + virtual double GetSingleValue(const InputImagePixelType & value1, const LabelledImagePixelType & value2) +{ +return 0; +} double GetValue(const InputImagePixelType & value1, const LabelledImagePixelType & value2) { @@ -156,15 +174,22 @@ class ITK_EXPORT MRFEnergy<TInput2,TInput2> : public itk::Object typedef itk::ConstNeighborhoodIterator< LabelledImageType > LabelledNeighborhoodIterator; typedef itk::Array< double > ParametersType; + /** Le new avait disparu...*/ itkNewMacro(Self); itkTypeMacro(MRFEnergy,itk::Object); + /** Utilisation des accesseur itk*/ + /* + void SetNumberOfParameters(unsigned int nb) + { m_NumberOfParameters = nb; this->Modified(); } unsigned int GetNumberOfParameters(void) const { return m_NumberOfParameters; } - void SetNumberOfParameters(unsigned int nb) - { m_NumberOfParameters = nb; this->Modified(); } + */ + itkSetMacro(NumberOfParameters, unsigned int); + itkGetConstMacro(NumberOfParameters, unsigned int); + // Get the parameters const ParametersType& GetParameters( void ) const { @@ -173,32 +198,33 @@ class ITK_EXPORT MRFEnergy<TInput2,TInput2> : public itk::Object void SetParameters( const ParametersType & parameters ) { - if( parameters.Size() != m_NumberOfParameters ) + if( parameters.Size() != m_NumberOfParameters ) { itkExceptionMacro(<<"Invalid number of parameters"); } m_Parameters = parameters; this->Modified(); - /* - bool modified = false; - for( unsigned int i=0; i<(unsigned int)m_NumberOfParameters; i++ ) + /* + bool modified = false; + for( unsigned int i=0; i<m_NumberOfParameters; i++ ) { - if (m_Parameters[i] != parameters[i]) - { - m_Parameters[i] = parameters[i]; - modified = true; - } + if (m_Parameters[i] != parameters[i]) + { + m_Parameters[i] = parameters[i]; + modified = true; + } } - if (modified) + if (modified) { - this->Modified(); - } - */ + this->Modified(); + }*/ } - virtual double GetSingleValue(const LabelledImagePixelType & value1, const LabelledImagePixelType & value2){return 0.;}; - - + /** Ne peut allouer un objet de type MRF::Energy parce que les fonctions virtuelles suivantes sont abstraites*/ + virtual double GetSingleValue(const LabelledImagePixelType & value1, const LabelledImagePixelType & value2) + { + return 0; + } double GetValue(const LabelledImagePixelType & value1, const LabelledImagePixelType & value2) { @@ -230,7 +256,10 @@ class ITK_EXPORT MRFEnergy<TInput2,TInput2> : public itk::Object protected: // The constructor and destructor. - MRFEnergy() {}; + MRFEnergy() + { + m_Parameters=0; + }; virtual ~MRFEnergy() {}; unsigned int m_NumberOfParameters; ParametersType m_Parameters; diff --git a/Code/Markov/otbMRFEnergyGaussianClassification.h b/Code/Markov/otbMRFEnergyGaussianClassification.h index 800b04cd35..d761f8e0d1 100644 --- a/Code/Markov/otbMRFEnergyGaussianClassification.h +++ b/Code/Markov/otbMRFEnergyGaussianClassification.h @@ -43,7 +43,7 @@ namespace otb { public: typedef MRFEnergyGaussianClassification Self; - typedef itk::Object Superclass; + typedef MRFEnergy< TInput1, TInput2> Superclass; typedef itk::SmartPointer<Self> Pointer; typedef itk::SmartPointer<const Self> ConstPointer; @@ -59,24 +59,25 @@ namespace otb itkTypeMacro(MRFEnergyGaussianClassification, MRFEnergy); void SetNumberOfParameters(unsigned int nParameters) - { - this->m_NumberOfParameters=nParameters; - this->m_Parameters.SetSize(this->m_NumberOfParameters); - this->Modified(); - } + { + this->m_NumberOfParameters=nParameters; + this->m_Parameters.SetSize(this->m_NumberOfParameters); + this->Modified(); + } double GetSingleValue(const InputImagePixelType & value1, const LabelledImagePixelType & value2) { - if ((unsigned int)value2 >= this->GetNumberOfParameters()/2) - { - itkExceptionMacro(<<"Number of parameters does not correspond to number of classes" ); - } - + if ((unsigned int)value2 >= this->GetNumberOfParameters()/2) + { + itkExceptionMacro(<<"Number of parameters does not correspond to number of classes" ); + } double val1 = static_cast<double>(value1); - double result = M_SQUARE(val1-this->m_Parameters[2*static_cast<int>(value2)])/(2*M_SQUARE(this->m_Parameters[2*static_cast<int>(value2)+1])) - + vcl_log(vcl_sqrt(2*M_PI)*this->m_Parameters[2*static_cast<int>(value2)+1]); - + double result = M_SQUARE(val1-this->m_Parameters[2*static_cast<int>(value2)]) + / (2*M_SQUARE(this->m_Parameters[2*static_cast<int>(value2)+1])) + + vcl_log(vcl_sqrt(2*M_PI) + *this->m_Parameters[2*static_cast<int>(value2)+1]); + return static_cast<double>( result ); } @@ -84,16 +85,7 @@ namespace otb protected: // The constructor and destructor. - MRFEnergyGaussianClassification() - { - this->SetNumberOfParameters ( 4 ); - this->m_Parameters.SetSize(this->GetNumberOfParameters()); - this->m_Parameters[0]=50.0; //Class 0 mean - this->m_Parameters[1]=10.0; //Class 0 stdev - this->m_Parameters[2]=140.0;//Class 1 mean - this->m_Parameters[3]=10.0; //Class 1 stdev - }; - + MRFEnergyGaussianClassification() {}; virtual ~MRFEnergyGaussianClassification() {}; }; diff --git a/Code/Markov/otbMRFEnergyPotts.h b/Code/Markov/otbMRFEnergyPotts.h index 295fee85fb..0d689e8994 100644 --- a/Code/Markov/otbMRFEnergyPotts.h +++ b/Code/Markov/otbMRFEnergyPotts.h @@ -42,7 +42,7 @@ class ITK_EXPORT MRFEnergyPotts:public MRFEnergy< TInput1, TInput2> { public: typedef MRFEnergyPotts Self; - typedef itk::Object Superclass; + typedef MRFEnergy< TInput1, TInput2> Superclass; typedef itk::SmartPointer<Self> Pointer; typedef itk::SmartPointer<const Self> ConstPointer; diff --git a/Code/Markov/otbMRFOptimizer.h b/Code/Markov/otbMRFOptimizer.h index 0471fce984..eb5ddc0cde 100644 --- a/Code/Markov/otbMRFOptimizer.h +++ b/Code/Markov/otbMRFOptimizer.h @@ -27,6 +27,8 @@ namespace otb /** * \class MRFOptimizer * \brief This is the base class for optimizer used in the MRF framework. + * + * Derived class must reimplement Compute() method. */ class ITK_EXPORT MRFOptimizer : public itk::Object @@ -38,12 +40,17 @@ class ITK_EXPORT MRFOptimizer : public itk::Object typedef itk::SmartPointer<const Self> ConstPointer; typedef itk::Array< double > ParametersType; - itkNewMacro(Self); - + /** Le new avait disparu donc plus de ::Pointer declare...*/ + itkNewMacro(Self); + itkTypeMacro(MRFOptimizer, itk::Object); + /** Utilmisation des accessuers itk*/ + /* unsigned int GetNumberOfParameters(void) const { return m_NumberOfParameters; } + */ + itkGetConstMacro(NumberOfParameters, unsigned int); // Get the parameters const ParametersType& GetParameters( void ) const @@ -59,30 +66,33 @@ class ITK_EXPORT MRFOptimizer : public itk::Object } m_Parameters = parameters; this->Modified(); - /* - bool modified = false; - for( unsigned int i=0; i<m_NumberOfParameters; i++ ) - { - if (m_Parameters[i] != parameters[i]) - { - m_Parameters[i] = parameters[i]; - modified = true; - } - } - if (modified) - { - this->Modified(); - } - */ + +// bool modified = false; +// for( unsigned int i=0; i<m_NumberOfParameters; i++ ) +// { +// if (m_Parameters[i] != parameters[i]) +// { +// m_Parameters[i] = parameters[i]; +// modified = true; +// } +// } +// if (modified) +// { +// this->Modified(); +// } } + /** Ne peut allouer un objet de type MRF::Optimizer parce que les fonctions virtuelles suivantes sont abstraites*/ virtual bool Compute(double deltaEnergy) - { - return false; - } +{ +return 0; +} protected: - MRFOptimizer() {} + MRFOptimizer() + { + m_Parameters=0; + } virtual ~MRFOptimizer() {} unsigned int m_NumberOfParameters; ParametersType m_Parameters; diff --git a/Code/Markov/otbMRFOptimizerICM.h b/Code/Markov/otbMRFOptimizerICM.h index 1902ba4a61..36120fd365 100644 --- a/Code/Markov/otbMRFOptimizerICM.h +++ b/Code/Markov/otbMRFOptimizerICM.h @@ -15,6 +15,7 @@ PURPOSE. See the above copyright notices for more information. =========================================================================*/ + #ifndef _MRFOptimizerICM_h #define _MRFOptimizerICM_h @@ -22,45 +23,47 @@ namespace otb { -/** - * \class MRFOptimizerICM - * \brief This is the optimizer class implementing the ICM algorithm - * - * This is one optimizer to be used in the MRF framework. This optimizer - * follows the ICM algorithm to accept of reject the value proposed by the sampler - */ -class ITK_EXPORT MRFOptimizerICM : public MRFOptimizer - { - public: - - typedef MRFOptimizerICM Self; - typedef itk::Object Superclass; - typedef itk::SmartPointer<Self> Pointer; - typedef itk::SmartPointer<const Self> ConstPointer; - - itkNewMacro(Self); - - itkTypeMacro(MRFOptimizerICM, MRFOptimizer); - - inline bool operator()(double deltaEnergy) - { - if (deltaEnergy < 0) - { - return true; - } - else - { - return false; - } - } - - - protected: - MRFOptimizerICM() {} - virtual ~MRFOptimizerICM() {} + /** + * \class MRFOptimizerICM + * \brief This is the optimizer class implementing the ICM algorithm + * + * This is one optimizer to be used in the MRF framework. This optimizer + * follows the ICM algorithm to accept of reject the value proposed by the sampler + */ + class ITK_EXPORT MRFOptimizerICM : + public MRFOptimizer + { + public: + + typedef MRFOptimizerICM Self; + typedef MRFOptimizer Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + + itkNewMacro(Self); - }; - + itkTypeMacro(MRFOptimizerICM,MRFOptimizer); + + + inline bool Compute(double deltaEnergy) + { + if (deltaEnergy < 0) + { + return true; + } + else + { + return false; + } + } + + + protected: + MRFOptimizerICM() {} + virtual ~MRFOptimizerICM() {} + + }; + } #endif diff --git a/Code/Markov/otbMRFOptimizerMetropolis.h b/Code/Markov/otbMRFOptimizerMetropolis.h index feefc461eb..622ffa72a4 100644 --- a/Code/Markov/otbMRFOptimizerMetropolis.h +++ b/Code/Markov/otbMRFOptimizerMetropolis.h @@ -24,16 +24,22 @@ namespace otb { - /** + /** * \class MRFOptimizerMetropolis * \brief This is the optimizer class implementing the Metropolis algorithm * * This is one optimizer to be used in the MRF framework. This optimizer - * follows the metropolis algorithm to accept of reject the value proposed by the sampler - */ + * follows the metropolis algorithm to accept of reject the value proposed by the sampler. + * + * The MRFOptimizerMetropolis has one parameter corresponding to the temperature T used + * to accept or reject proposed values. The proposed value is accepted with a probability: + * + * \f[ e^{\frac{-\Delta E}{T}} \f] + * + */ class ITK_EXPORT MRFOptimizerMetropolis : public MRFOptimizer - { + { public: typedef MRFOptimizerMetropolis Self; @@ -41,20 +47,10 @@ class ITK_EXPORT MRFOptimizerMetropolis : public MRFOptimizer typedef itk::SmartPointer<Self> Pointer; typedef itk::SmartPointer<const Self> ConstPointer; - typedef Superclass::ParametersType ParametersType; - itkNewMacro(Self); itkTypeMacro(MRFOptimizerMetropolis,MRFOptimizer); - /** Store a value to be used instead of random value.. FOR TEST ONLY*/ - void SetValueInsteadRandom( double val ) - { - m_ValueInsteadRandom = val; - std::cout<<"The m_ValueInsteadRandom varaible has to be used only for tests..."<<std::endl; - this->Modified(); - }; - inline bool Compute(double deltaEnergy) { if (deltaEnergy < 0) @@ -66,38 +62,26 @@ class ITK_EXPORT MRFOptimizerMetropolis : public MRFOptimizer return false; } else - { - double proba = vcl_exp(-(deltaEnergy)/this->m_Parameters[0]); - double val; - if( m_ValueInsteadRandom==itk::NumericTraits<double>::min() ) - { - val = (rand() % 10000); - } - else - { - val = m_ValueInsteadRandom; - } - if ( val < proba*10000) - { - return true; - } - } + { + double proba = vcl_exp(-(deltaEnergy)/this->m_Parameters[0]); + if ( (rand() % 10000) < proba*10000) + { + std::cerr << "Opti true " << std::endl; + return true; + } + } return false; } protected: - MRFOptimizerMetropolis() - { - this->m_NumberOfParameters = 1; - this->m_Parameters.SetSize(this->m_NumberOfParameters); - this->m_Parameters[0]=1.0; - - m_ValueInsteadRandom=itk::NumericTraits<double>::min(); - } + MRFOptimizerMetropolis() { + this->m_NumberOfParameters = 1; + this->m_Parameters.SetSize(this->m_NumberOfParameters); + this->m_Parameters[0]=1.0; + } virtual ~MRFOptimizerMetropolis() {} - /** Store a value to be used instead of random value.. FOR TEST ONLY*/ - double m_ValueInsteadRandom; + double m_Temperature; }; } diff --git a/Code/Markov/otbMRFSampler.h b/Code/Markov/otbMRFSampler.h index f195e35faa..aace58bd1f 100644 --- a/Code/Markov/otbMRFSampler.h +++ b/Code/Markov/otbMRFSampler.h @@ -19,104 +19,89 @@ #define _MRFSampler_h #include "otbMRFEnergy.h" -// #include "otbMRFEnergyFunctorBase.h" #include "itkNeighborhoodIterator.h" namespace otb { -/** - * \class MRFSampler - * \brief This is the base class for sampler methods used in the MRF framework. - */ -template< class TInput1, class TInput2> -class ITK_EXPORT MRFSampler:public itk::Object - { - public: - typedef MRFSampler Self; - typedef itk::Object Superclass; - typedef itk::SmartPointer<Self> Pointer; - typedef itk::SmartPointer<const Self> ConstPointer; - - // typedef itk::NeighborhoodIterator< TInput1 > InputImageNeighborhoodIterator; - typedef itk::ConstNeighborhoodIterator< TInput1 > InputImageNeighborhoodIterator; - typedef itk::NeighborhoodIterator< TInput2 > LabelledImageNeighborhoodIterator; - typedef typename TInput2::PixelType LabelledImagePixelType; - typedef typename TInput1::PixelType InputImagePixelType; - - // typedef MRFEnergyFunctorBase<TInput1, TInput2> EnergyFidelityFunctorType; - // typedef MRFEnergy<EnergyFidelityFunctorType> EnergyFidelityType; - // typedef MRFEnergyFunctorBase<TInput2, TInput2> EnergyRegularizationFunctorType; - // typedef MRFEnergy<EnergyRegularizationFunctorType> EnergyRegularizationType; - typedef MRFEnergy<TInput1, TInput2> EnergyFidelityType; - typedef MRFEnergy<TInput2, TInput2> EnergyRegularizationType; - typedef typename EnergyFidelityType::Pointer EnergyFidelityPointer; - typedef typename EnergyRegularizationType::Pointer EnergyRegularizationPointer; - - itkNewMacro(Self); + /** + * \class MRFSampler + * \brief This is the base class for sampler methods used in the MRF framework. + * + * Derived class must reimplement Compute() method. + * + */ - itkTypeMacro(MRFSampler,itk::Object); - - // Accessors - virtual void SetNumberOfClasses(unsigned int nb){ m_NumberOfClasses = nb; this->Modified(); }; - itkGetMacro(NumberOfClasses, unsigned int); + template< class TInput1, class TInput2> + class ITK_EXPORT MRFSampler:public itk::Object + { + public: + + typedef MRFSampler Self; + typedef itk::Object Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + + typedef itk::ConstNeighborhoodIterator< TInput1 > InputImageNeighborhoodIterator; + typedef typename TInput1::PixelType InputImagePixelType; - itkSetMacro(Lambda, double); - itkGetMacro(Lambda, double); - itkSetMacro(EnergyBefore, double); - itkGetMacro(EnergyBefore, double); - itkSetMacro(EnergyAfter, double); - itkGetMacro(EnergyAfter, double); - itkSetMacro(DeltaEnergy, double); - itkGetMacro(DeltaEnergy, double); - itkSetMacro(EnergyCurrent, double); - itkGetMacro(EnergyCurrent, double); - - itkSetMacro(Value, LabelledImagePixelType); - itkGetMacro(Value, LabelledImagePixelType); - - void SetEnergyRegularization( EnergyRegularizationPointer enerReg ){ m_EnergyRegularization=enerReg; this->Modified(); }; - EnergyRegularizationPointer GetEnergyRegularization(){ return m_EnergyRegularization; }; - void SetEnergyFidelity( EnergyFidelityPointer enerFid ){ m_EnergyFidelity=enerFid; this->Modified(); }; - EnergyFidelityPointer GetEnergyFidelity(){ return m_EnergyFidelity; }; - - virtual int Compute( const InputImageNeighborhoodIterator & itData, const LabelledImageNeighborhoodIterator & itRegul) - { - m_DeltaEnergy = 0; - - return 0; - } - - - protected: - unsigned int m_NumberOfClasses; - double m_EnergyBefore; - double m_EnergyAfter; - double m_DeltaEnergy; - double m_EnergyCurrent; - double m_Lambda; - LabelledImagePixelType m_Value; - EnergyRegularizationPointer m_EnergyRegularization; - EnergyFidelityPointer m_EnergyFidelity; - LabelledImagePixelType m_ValueCurrent; - - - protected: - // The constructor and destructor. - MRFSampler() - { - m_EnergyRegularization = EnergyRegularizationType::New(); - m_EnergyFidelity = EnergyFidelityType::New(); + typedef itk::NeighborhoodIterator< TInput2 > LabelledImageNeighborhoodIterator; + typedef typename TInput2::PixelType LabelledImagePixelType; + + + typedef MRFEnergy<TInput1, TInput2> EnergyFidelityType; + typedef MRFEnergy<TInput2, TInput2> EnergyRegularizationType; + + + typedef typename EnergyFidelityType::Pointer EnergyFidelityPointer; + typedef typename EnergyRegularizationType::Pointer EnergyRegularizationPointer; + + itkTypeMacro(MRFSampler,Object); + + itkSetMacro(NumberOfClasses, unsigned int); + itkGetMacro(NumberOfClasses, unsigned int); + + itkSetMacro(Lambda, double); + itkGetMacro(Lambda, double); + + itkGetMacro(DeltaEnergy, double); + itkGetMacro(Value, LabelledImagePixelType); + + /**Accesseur pour les tests*/ + itkGetConstMacro(EnergyBefore, double); + itkGetConstMacro(EnergyAfter, double); + + void SetEnergyRegularization( EnergyRegularizationPointer enerReg ){ m_EnergyRegularization=enerReg; this->Modified(); }; + void SetEnergyFidelity( EnergyFidelityPointer enerFid ){ m_EnergyFidelity=enerFid; this->Modified(); }; + /* + itkSetObjectMacro( EnergyRegularization, EnergyRegularizationPointer); + itkSetObjectMacro( EnergyFidelity, EnergyFidelityPointer); + */ + virtual int Compute( const InputImageNeighborhoodIterator & itData, + const LabelledImageNeighborhoodIterator & itRegul) = 0; + + protected: + unsigned int m_NumberOfClasses; + double m_EnergyBefore; + double m_EnergyAfter; + double m_DeltaEnergy; + double m_EnergyCurrent; + double m_Lambda; + + LabelledImagePixelType m_Value; + EnergyRegularizationPointer m_EnergyRegularization; + EnergyFidelityPointer m_EnergyFidelity; + LabelledImagePixelType m_ValueCurrent; - m_NumberOfClasses = 1; - m_Lambda = 1.; - m_EnergyCurrent = 1.; - m_DeltaEnergy = 1.; - m_EnergyAfter = 1.; - m_EnergyBefore = 1.; - } - virtual ~MRFSampler() {} + + + protected: + // The constructor and destructor. + MRFSampler() {} + virtual ~MRFSampler() {} - }; + }; + + } #endif diff --git a/Code/Markov/otbMRFSamplerMAP.h b/Code/Markov/otbMRFSamplerMAP.h index 5356fa8a75..d965e9cc91 100644 --- a/Code/Markov/otbMRFSamplerMAP.h +++ b/Code/Markov/otbMRFSamplerMAP.h @@ -15,6 +15,7 @@ PURPOSE. See the above copyright notices for more information. =========================================================================*/ + #ifndef _MRFSamplerMAP_h #define _MRFSamplerMAP_h @@ -22,81 +23,90 @@ namespace otb { -/** - * \class MRFSamplerMAP - * \brief This is the base class for sampler methods used in the MRF framework. - * - * This is one sampler to be used int he MRF framework. This sampler select the - * value which maximizes the apriori probability (minimum energy). - */ - -template< class TInput1, class TInput2> -class ITK_EXPORT MRFSamplerMAP : public MRFSampler< TInput1, TInput2 > -{ - public: - - typedef MRFSamplerMAP Self; - typedef MRFSampler< TInput1, TInput2 > Superclass; - typedef itk::SmartPointer<Self> Pointer; - typedef itk::SmartPointer<const Self> ConstPointer; - - typedef typename Superclass::InputImageNeighborhoodIterator InputImageNeighborhoodIterator; - typedef typename Superclass::LabelledImageNeighborhoodIterator LabelledImageNeighborhoodIterator; - typedef typename Superclass::LabelledImagePixelType LabelledImagePixelType; - typedef typename Superclass::InputImagePixelType InputImagePixelType; - typedef typename Superclass::EnergyFidelityType EnergyFidelityType; - typedef typename Superclass::EnergyRegularizationType EnergyRegularizationType; - typedef typename Superclass::EnergyFidelityPointer EnergyFidelityPointer; - typedef typename Superclass::EnergyRegularizationPointer EnergyRegularizationPointer; - - itkNewMacro(Self); - - itkTypeMacro(MRFSamplerMAP,MRFSampler); - + /** + * \class MRFSamplerMAP + * \brief This is the base class for sampler methods used in the MRF framework. + * + * This is one sampler to be used int he MRF framework. This sampler select the + * value which maximizes the apriori probability (minimum energy). + * + */ + + template< class TInput1, class TInput2> + class ITK_EXPORT MRFSamplerMAP : public MRFSampler< TInput1, TInput2> + { + public: - inline int Compute( const InputImageNeighborhoodIterator & itData, const LabelledImageNeighborhoodIterator & itRegul) - { - this->SetEnergyBefore( this->GetEnergyFidelity()->GetValue(itData, itRegul.GetCenterPixel()) - + this->GetLambda() * this->GetEnergyRegularization()->GetValue(itRegul, itRegul.GetCenterPixel()) ); - - //Try all possible value (how to be generic ?) - this->SetEnergyAfter( this->GetEnergyBefore() ); //default values to current one - this->SetValue( itRegul.GetCenterPixel() ); - - LabelledImagePixelType valueCurrent = 0; - while( valueCurrent<static_cast<LabelledImagePixelType>(this->GetNumberOfClasses()) && valueCurrent != itk::NumericTraits<LabelledImagePixelType>::max() ) - { - this->SetEnergyCurrent( this->GetEnergyFidelity()->GetValue(itData, valueCurrent) - + this->GetLambda() - * this->GetEnergyRegularization()->GetValue(itRegul, valueCurrent) ); - if ( this->GetEnergyCurrent() < this->GetEnergyAfter() ) - { - this->SetEnergyAfter( this->GetEnergyCurrent() ); - this->SetValue( valueCurrent ); - } - valueCurrent++; - } - - // TODO avoir la confirmation cnesienne : premier indince ou dernier - if ( valueCurrent==itk::NumericTraits<LabelledImagePixelType>::max() ) - { - valueCurrent = 0; - } - - this->SetDeltaEnergy( this->GetEnergyAfter() - this->GetEnergyBefore() ); - - return 0; - } - - - protected: - // The constructor and destructor. - MRFSamplerMAP() {} - virtual ~MRFSamplerMAP() {} + typedef MRFSamplerMAP Self; + typedef MRFSampler< TInput1, TInput2> Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + +// typedef itk::ConstNeighborhoodIterator< TInput1 > InputImageNeighborhoodIterator; +// typedef itk::NeighborhoodIterator< TInput2 > LabelledImageNeighborhoodIterator; +// typedef typename TInput2::PixelType LabelledImagePixelType; +// +// typedef MRFEnergy<TInput1, TInput2> EnergyFidelityType; +// typedef MRFEnergy<TInput2, TInput2> EnergyRegularizationType; +// +// typedef typename EnergyFidelityType::Pointer EnergyFidelityPointer; +// typedef typename EnergyRegularizationType::Pointer EnergyRegularizationPointer; + + typedef typename Superclass::InputImageNeighborhoodIterator InputImageNeighborhoodIterator; + typedef typename Superclass::LabelledImageNeighborhoodIterator LabelledImageNeighborhoodIterator; + typedef typename Superclass::LabelledImagePixelType LabelledImagePixelType; + typedef typename Superclass::InputImagePixelType InputImagePixelType; + typedef typename Superclass::EnergyFidelityType EnergyFidelityType; + typedef typename Superclass::EnergyRegularizationType EnergyRegularizationType; + typedef typename Superclass::EnergyFidelityPointer EnergyFidelityPointer; + typedef typename Superclass::EnergyRegularizationPointer EnergyRegularizationPointer; + + itkNewMacro(Self); + + itkTypeMacro(MRFSamplerMAP,MRFSampler); + + + inline int Compute( const InputImageNeighborhoodIterator & itData, + const LabelledImageNeighborhoodIterator & itRegul) + { + this->m_EnergyBefore=this->m_EnergyFidelity->GetValue(itData, itRegul.GetCenterPixel()); + this->m_EnergyBefore += this->m_Lambda + * this->m_EnergyRegularization->GetValue(itRegul, itRegul.GetCenterPixel()); + + //Try all possible value (how to be generic ?) + this->m_EnergyAfter = this->m_EnergyBefore; //default values to current one + this->m_Value = itRegul.GetCenterPixel(); +// for (LabelledImagePixelType valueCurrent = 0; +// valueCurrent< this->m_NumberOfClasses; ++valueCurrent) +// { + LabelledImagePixelType valueCurrent = 0; + while( valueCurrent<static_cast<LabelledImagePixelType>(this->GetNumberOfClasses()) && valueCurrent != itk::NumericTraits<LabelledImagePixelType>::max() ) + { + this->m_EnergyCurrent = this->m_EnergyFidelity->GetValue(itData, valueCurrent); + this->m_EnergyCurrent += this->m_Lambda + * this->m_EnergyRegularization->GetValue(itRegul, valueCurrent); + if (this->m_EnergyCurrent < this->m_EnergyAfter) + { + this->m_EnergyAfter = this->m_EnergyCurrent; + this->m_Value = valueCurrent; + } + valueCurrent++; + } + + this->m_DeltaEnergy= this->m_EnergyAfter - this->m_EnergyBefore; + + return 0; + } + + + protected: + // The constructor and destructor. + MRFSamplerMAP() {} + virtual ~MRFSamplerMAP() {} + + }; -}; - } #endif diff --git a/Code/Markov/otbMRFSamplerRandom.h b/Code/Markov/otbMRFSamplerRandom.h index 2c497edf39..3ec6871812 100644 --- a/Code/Markov/otbMRFSamplerRandom.h +++ b/Code/Markov/otbMRFSamplerRandom.h @@ -19,87 +19,82 @@ #define _MRFSamplerRandom_h #include "otbMRFSampler.h" +#include "itkMersenneTwisterRandomVariateGenerator.h" #include "itkNumericTraits.h" namespace otb { -/** - * \class MRFSamplerRandom - * \brief This is the base class for sampler methods used in the MRF framework. - * - * This is one sampler to be used int he MRF framework. This sampler select the - * value randomly according to a uniform probability. - * - */ - -template< class TInput1, class TInput2> -class ITK_EXPORT MRFSamplerRandom : public MRFSampler< TInput1, TInput2> - { - public: - - typedef MRFSamplerRandom Self; - typedef MRFSampler<TInput1, TInput2> Superclass; - typedef itk::SmartPointer<Self> Pointer; - typedef itk::SmartPointer<const Self> ConstPointer; - - typedef typename Superclass::InputImageNeighborhoodIterator InputImageNeighborhoodIterator; - typedef typename Superclass::LabelledImageNeighborhoodIterator LabelledImageNeighborhoodIterator; - typedef typename Superclass::LabelledImagePixelType LabelledImagePixelType; - typedef typename Superclass::InputImagePixelType InputImagePixelType; - typedef typename Superclass::EnergyFidelityType EnergyFidelityType; - typedef typename Superclass::EnergyRegularizationType EnergyRegularizationType; - typedef typename Superclass::EnergyFidelityPointer EnergyFidelityPointer; - typedef typename Superclass::EnergyRegularizationPointer EnergyRegularizationPointer; - // Now use Superclass tepedefs - //typedef itk::ConstNeighborhoodIterator< TInput1 > InputImageNeighborhoodIterator; - //typedef itk::NeighborhoodIterator< TInput2 > LabelledImageNeighborhoodIterator; - //typedef typename TInput2::PixelType LabelledImagePixelType; - //typedef MRFEnergy<TInput1, TInput2> EnergyFidelityType; - //typedef MRFEnergy<TInput2, TInput2> EnergyRegularizationType; - //typedef typename EnergyFidelityType::Pointer EnergyFidelityPointer; - //typedef typename EnergyRegularizationType::Pointer EnergyRegularizationPointer; - - itkNewMacro(Self); - - itkTypeMacro(MRFSamplerRandom,MRFSampler); + /** + * \class MRFSamplerRandom + * \brief This is the base class for sampler methods used in the MRF framework. + * + * This is one sampler to be used int he MRF framework. This sampler select the + * value randomly according to a uniform probability. + * + */ - /** Store a value to be used instead of random value.. FOR TEST ONLY*/ - void SetValueInsteadRandom( LabelledImagePixelType val ) + template< class TInput1, class TInput2> + class ITK_EXPORT MRFSamplerRandom: public MRFSampler< TInput1, TInput2> { - m_ValueInsteadRandom = val; - std::cout<<"The m_ValueInsteadRandom varaible has to be used only for tests..."<<std::endl; - this->Modified(); - }; + public: + typedef MRFSamplerRandom Self; + typedef otb::MRFSampler< TInput1, TInput2> Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + + typedef typename Superclass::InputImageNeighborhoodIterator InputImageNeighborhoodIterator; + typedef typename Superclass::LabelledImageNeighborhoodIterator LabelledImageNeighborhoodIterator; + typedef typename Superclass::LabelledImagePixelType LabelledImagePixelType; + typedef typename Superclass::InputImagePixelType InputImagePixelType; + typedef typename Superclass::EnergyFidelityType EnergyFidelityType; + typedef typename Superclass::EnergyRegularizationType EnergyRegularizationType; + typedef typename Superclass::EnergyFidelityPointer EnergyFidelityPointer; + typedef typename Superclass::EnergyRegularizationPointer EnergyRegularizationPointer; + + typedef itk::Statistics::MersenneTwisterRandomVariateGenerator RandomGeneratorType; - inline int Compute( const InputImageNeighborhoodIterator & itData, const LabelledImageNeighborhoodIterator & itRegul) - { - this->SetEnergyBefore ( this->GetEnergyFidelity()->GetValue(itData, itRegul.GetCenterPixel()) ); - this->m_EnergyBefore += this->GetLambda() * this->GetEnergyRegularization()->GetValue(itRegul, itRegul.GetCenterPixel()); - LabelledImagePixelType val; - if (m_ValueInsteadRandom == itk::NumericTraits<LabelledImagePixelType>::min()) - { - val = static_cast<LabelledImagePixelType> (rand() % this->GetNumberOfClasses()); - } - else - { - val = m_ValueInsteadRandom; - } + itkNewMacro(Self); - this->SetValue( val ); - this->SetEnergyAfter ( this->GetEnergyFidelity()->GetValue(itData, this->GetValue()) ); - this->m_EnergyAfter += this->GetLambda() * this->GetEnergyRegularization()->GetValue(itRegul, this->GetValue()); - this->SetDeltaEnergy( this->GetEnergyAfter() - this->GetEnergyBefore() ); + itkTypeMacro(MRFSamplerRandom,MRFSampler); - return 0; - } + inline int Compute( const InputImageNeighborhoodIterator & itData, + const LabelledImageNeighborhoodIterator & itRegul) + { + this->m_EnergyBefore = + this->m_EnergyFidelity->GetValue(itData, itRegul.GetCenterPixel()); + this->m_EnergyBefore += this->m_Lambda + * this->m_EnergyRegularization->GetValue(itRegul, itRegul.GetCenterPixel()); + + + this->m_Value = + static_cast<LabelledImagePixelType> + (m_Generator->GetIntegerVariate() % this->m_NumberOfClasses); + + this->m_EnergyAfter = + this->m_EnergyFidelity->GetValue(itData, this->m_Value); + this->m_EnergyAfter += this->m_Lambda + * this->m_EnergyRegularization->GetValue(itRegul, this->m_Value); - protected: - // The constructor and destructor. - MRFSamplerRandom() { m_ValueInsteadRandom=itk::NumericTraits<LabelledImagePixelType>::min(); } - virtual ~MRFSamplerRandom() {} - /** Store a value to be used instead of random value.. FOR TEST ONLY*/ - LabelledImagePixelType m_ValueInsteadRandom; - }; + this->m_DeltaEnergy= this->m_EnergyAfter - this->m_EnergyBefore; + + return 0; + } + + + protected: +// The constructor and destructor. + MRFSamplerRandom() { + m_Generator = + RandomGeneratorType::New(); + m_Generator->SetSeed(); + } + virtual ~MRFSamplerRandom() {} + + private: + RandomGeneratorType::Pointer m_Generator; + }; + + } #endif diff --git a/Code/Markov/otbMRFSamplerRandomMAP.h b/Code/Markov/otbMRFSamplerRandomMAP.h index 58020d6a05..801b2c303b 100644 --- a/Code/Markov/otbMRFSamplerRandomMAP.h +++ b/Code/Markov/otbMRFSamplerRandomMAP.h @@ -15,6 +15,7 @@ PURPOSE. See the above copyright notices for more information. =========================================================================*/ + #ifndef _MRFSamplerRandomMAP_h #define _MRFSamplerRandomMAP_h @@ -22,180 +23,138 @@ namespace otb { -/** - * \class MRFSamplerRandomMAP - * \brief This is the base class for sampler methods used in the MRF framework. - * - * This is one sampler to be used int he MRF framework. This sampler select the - * value randomly according to the apriori probability. - * - * The probability is defined from the energy as: - * - * \f[ P(X=x)= \frac{1}{Z} \exp^{-U(x)} \f] - * - * where \f$ Z = \sum_x \exp^{-U(x)}\f$ - * - */ - -template< class TInput1, class TInput2> -class ITK_EXPORT MRFSamplerRandomMAP : public MRFSampler< TInput1, TInput2> - { - public: - - typedef MRFSamplerRandomMAP Self; - typedef MRFSampler< TInput1, TInput2> Superclass; - typedef itk::SmartPointer<Self> Pointer; - typedef itk::SmartPointer<const Self> ConstPointer; - - typedef typename Superclass::InputImageNeighborhoodIterator InputImageNeighborhoodIterator; - typedef typename Superclass::LabelledImageNeighborhoodIterator LabelledImageNeighborhoodIterator; - typedef typename Superclass::LabelledImagePixelType LabelledImagePixelType; - typedef typename Superclass::InputImagePixelType InputImagePixelType; - typedef typename Superclass::EnergyFidelityType EnergyFidelityType; - typedef typename Superclass::EnergyRegularizationType EnergyRegularizationType; - typedef typename Superclass::EnergyFidelityPointer EnergyFidelityPointer; - typedef typename Superclass::EnergyRegularizationPointer EnergyRegularizationPointer; - - itkNewMacro(Self); - - itkTypeMacro(MRFSamplerRandomMAP, MRFSampler); - - /* - void SetNumberOfClasses(unsigned int arg) - { - otbDebugMacro("setting NumberOfClasses to " << arg); - if (this->m_NumberOfClasses != arg) - { - this->m_NumberOfClasses = arg; - m_RepartitionFunction = (double *) calloc(m_NumberOfClasses, sizeof(double)); - m_Energy = (double *) calloc(m_NumberOfClasses, sizeof(double)); - this->Modified(); - } - } - itkGetMacro(NumberOfClasses, unsigned int); - - itkSetMacro(Lambda, double); - itkGetMacro(Lambda, double); - - itkGetMacro(DeltaEnergy, double); - itkGetMacro(Value, LabelledImagePixelType); - - itkSetObjectMacro( EnergyRegularization, EnergyRegularizationType); - itkSetObjectMacro( EnergyFidelity, EnergyFidelityType); - */ - - inline int Compute( const InputImageNeighborhoodIterator & itData, const LabelledImageNeighborhoodIterator & itRegul) + /** + * \class MRFSamplerRandomMAP + * \brief This is the base class for sampler methods used in the MRF framework. + * + * This is one sampler to be used int he MRF framework. This sampler select the + * value randomly according to the apriori probability. + * + * The probability is defined from the energy as: + * + * \f[ P(X=x)= \frac{1}{Z} \exp^{-U(x)} \f] + * + * where \f$ Z = \sum_x \exp^{-U(x)}\f$ + * + */ + + template< class TInput1, class TInput2> + class ITK_EXPORT MRFSamplerRandomMAP: public MRFSampler< TInput1, TInput2> { - if(this->GetNumberOfClasses() == 0) - { - itkExceptionMacro(<<"m_NumberOfClasses has to be greater than 0."); - } + public: + + typedef MRFSamplerRandomMAP Self; + typedef otb::MRFSampler< TInput1, TInput2> Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + + typedef typename Superclass::InputImageNeighborhoodIterator InputImageNeighborhoodIterator; + typedef typename Superclass::LabelledImageNeighborhoodIterator LabelledImageNeighborhoodIterator; + typedef typename Superclass::LabelledImagePixelType LabelledImagePixelType; + typedef typename Superclass::InputImagePixelType InputImagePixelType; + typedef typename Superclass::EnergyFidelityType EnergyFidelityType; + typedef typename Superclass::EnergyRegularizationType EnergyRegularizationType; + typedef typename Superclass::EnergyFidelityPointer EnergyFidelityPointer; + typedef typename Superclass::EnergyRegularizationPointer EnergyRegularizationPointer; - this->SetEnergyBefore( this->GetEnergyFidelity()->GetValue(itData, itRegul.GetCenterPixel()) - + this->GetLambda() * this->GetEnergyRegularization()->GetValue(itRegul, itRegul.GetCenterPixel()) ); + itkNewMacro(Self); + + itkTypeMacro(MRFSamplerRandomMAP,MRFSampler); + + + void SetNumberOfClasses(const unsigned int nClasses) + { + if (nClasses != this->m_NumberOfClasses) + { + this->m_NumberOfClasses = nClasses; + if (energy != NULL) free(energy); + if (repartitionFunction != NULL) free(repartitionFunction); + energy = (double *) calloc(this->m_NumberOfClasses, sizeof(double)); + repartitionFunction = (double *) calloc(this->m_NumberOfClasses, sizeof(double)); + this->Modified(); + } + } + + inline int Compute( const InputImageNeighborhoodIterator & itData, + const LabelledImageNeighborhoodIterator & itRegul) + { + + this->m_EnergyBefore = this->m_EnergyFidelity->GetValue(itData, itRegul.GetCenterPixel()); + this->m_EnergyBefore += this->m_Lambda + * this->m_EnergyRegularization->GetValue(itRegul, itRegul.GetCenterPixel()); + + //Try all possible value (how to be generic ?) + this->m_EnergyAfter = this->m_EnergyBefore; //default values to current one + this->m_Value = itRegul.GetCenterPixel(); + + //Compute probability for each possibility + double totalProba=0.0; + for (LabelledImagePixelType valueCurrent = 0; + valueCurrent < static_cast<LabelledImagePixelType>(this->m_NumberOfClasses); ++valueCurrent) + { + this->m_EnergyCurrent = this->m_EnergyFidelity->GetValue(itData, valueCurrent); + this->m_EnergyCurrent += this->m_Lambda + * this->m_EnergyRegularization->GetValue(itRegul, valueCurrent); - //Try all possible value (how to be generic ?) - this->SetEnergyAfter( this->GetEnergyBefore() ); //default values to current one - this->SetValue( itRegul.GetCenterPixel() ); + energy[static_cast<unsigned int>(valueCurrent)] = this->m_EnergyCurrent; + repartitionFunction[static_cast<unsigned int>(valueCurrent)] = vcl_exp(-this->m_EnergyCurrent)+totalProba; + totalProba = repartitionFunction[static_cast<unsigned int>(valueCurrent)]; - // otbDebugMacro(<< "Computing MAP for pix " << itData.GetIndex()); - // Compute probability for each possibility - double totalProba=0.0; - for (unsigned int valueCurrent = 0; valueCurrent < this->GetNumberOfClasses(); ++valueCurrent) - { - // otbDebugMacro(<< " --> Proposed value " << static_cast<double>(valueCurrent)); - this->SetEnergyCurrent( this->GetEnergyFidelity()->GetValue(itData, valueCurrent) - + this->GetLambda() * this->GetEnergyRegularization()->GetValue(itRegul, valueCurrent) ); - - m_Energy[valueCurrent] = this->GetEnergyCurrent(); - m_RepartitionFunction[valueCurrent] = vcl_exp(-this->GetEnergyCurrent())+totalProba; - totalProba = m_RepartitionFunction[valueCurrent]; - // otbDebugMacro("valueCurrent, m_RepartitionFunction[valueCurrent] " << (unsigned int) valueCurrent << ", " << m_RepartitionFunction[valueCurrent]); - - } - - //Pick a value according to probability - - double select; - if (m_ValueInsteadRandom == itk::NumericTraits<double>::min()) - { - select = (rand()/(double(RAND_MAX)+1) * totalProba); - } - else - { - select = m_ValueInsteadRandom; - } - // otbDebugMacro("m_RepartitionFunction " << m_RepartitionFunction[0] << " " - // << m_RepartitionFunction[1] << " " << m_RepartitionFunction[2] << " " - // << m_RepartitionFunction[3] << " "); - - // otbDebugMacro("select, totalProba " << select << ", " << totalProba); - unsigned int valueCurrent = 0; - while( valueCurrent<this->GetNumberOfClasses() && m_RepartitionFunction[valueCurrent] <= select) - { - valueCurrent++; - } + } + + //Pick a value according to probability + + double select = (rand()/(double(RAND_MAX)+1) * totalProba); + + /*** POUR EVITER LE SEGFAULT : on sort de la taille du tableau si la condition du break n'est pas respectee...*/ + /* + LabelledImagePixelType valueCurrent = 0; + for (valueCurrent = 0; + valueCurrent < static_cast<LabelledImagePixelType>(this->m_NumberOfClasses); ++valueCurrent) + { + if (repartitionFunction[static_cast<unsigned int>(valueCurrent)] > select) break; + } + */ + unsigned int valueCurrent = 0; + while( valueCurrent<this->GetNumberOfClasses() && repartitionFunction[valueCurrent] <= select) + { + valueCurrent++; + } - // TODO avoir la confirmation cnesienne : premier indince ou dernier - if ( valueCurrent==this->GetNumberOfClasses() ) - { - valueCurrent = 0; - } - - if ( this->GetValue() != static_cast<LabelledImagePixelType>(valueCurrent)) - { - this->SetValue( valueCurrent ); - this->SetEnergyAfter( m_Energy[valueCurrent] ); - } - - this->SetDeltaEnergy( this->GetEnergyAfter() - this->GetEnergyBefore() ); - // otbDebugMacro("Decision " << (unsigned int) valueCurrent); - return 0; - } - - /** Store a value to be used instead of random value.. FOR TEST ONLY*/ - void SetValueInsteadRandom( double val ) - { - m_ValueInsteadRandom = val; - std::cout<<"The m_ValueInsteadRandom varaible has to be used only for tests..."<<std::endl; - this->Modified(); - }; - - void SetNumberOfClasses(unsigned int nb) + // TODO avoir la confirmation cnesienne : premier indince ou dernier + if ( valueCurrent==this->GetNumberOfClasses() ) { - Superclass::SetNumberOfClasses( nb ); - free(m_Energy); - free(m_RepartitionFunction); - m_Energy = (double *) calloc(nb, sizeof(double)); - m_RepartitionFunction = (double *) calloc(nb, sizeof(double)); - - this->Modified(); - }; - - private: - double * m_RepartitionFunction ; - double * m_Energy; - - - - protected: - // The constructor and destructor. - MRFSamplerRandomMAP() - { - m_Energy = (double *) calloc(this->GetNumberOfClasses(), sizeof(double)); - m_RepartitionFunction = (double *) calloc(this->GetNumberOfClasses(), sizeof(double)); - m_ValueInsteadRandom = itk::NumericTraits<double>::min(); - } - virtual ~MRFSamplerRandomMAP() - { - free(m_Energy); - free(m_RepartitionFunction); - } + valueCurrent = 0; + } - /** Store a value to be used instead of random value.. FOR TEST ONLY*/ - double m_ValueInsteadRandom; + if ( this->m_Value != static_cast<LabelledImagePixelType>(valueCurrent)) + { + this->m_Value = static_cast<LabelledImagePixelType>(valueCurrent); + this->m_EnergyAfter = energy[valueCurrent]; + } + + this->m_DeltaEnergy= this->m_EnergyAfter - this->m_EnergyBefore; + + return 0; + } + + private: + double * repartitionFunction; + double * energy; + + protected: + // The constructor and destructor. + MRFSamplerRandomMAP() { + energy=NULL; + repartitionFunction=NULL; + }; + virtual ~MRFSamplerRandomMAP() { + if (energy != NULL) free(energy); + if (repartitionFunction != NULL) free(repartitionFunction); + }; - }; + }; + } diff --git a/Code/Markov/otbMarkovRandomFieldFilter.h b/Code/Markov/otbMarkovRandomFieldFilter.h index 345ea68038..65bfab382d 100644 --- a/Code/Markov/otbMarkovRandomFieldFilter.h +++ b/Code/Markov/otbMarkovRandomFieldFilter.h @@ -38,311 +38,325 @@ #include "otbMRFOptimizer.h" #include "otbMRFSampler.h" -//should be done with generic base class later... -// #include "otbGaussianMRFEnergy.h" -// #include "otbMRFSamplerRandom.h" -// #include "otbMRFOptimizerMetropolis.h" -// #include "otbGaussianClassesMRFEnergy.h" -// #include "otbMRFSamplerMAP.h" -// #include "otbMRFOptimizerICM.h" - - namespace otb { template <class TInputImage, class TClassifiedImage> class ITK_EXPORT MarkovRandomFieldFilter : public itk::ImageToImageFilter<TInputImage,TClassifiedImage> -{ - public: - /** Standard class typedefs. */ - typedef MarkovRandomFieldFilter Self; - typedef itk::ImageToImageFilter<TInputImage,TClassifiedImage> Superclass; - typedef itk::SmartPointer<Self> Pointer; - typedef itk::SmartPointer<const Self> ConstPointer; - typedef typename Superclass::OutputImagePointer OutputImagePointer; - - /** Run-time type information (and related methods). */ - itkTypeMacro(MarkovRandomFieldFilter,itk::ImageToImageFilter); - - /** Method for creation through the object factory. */ - itkNewMacro(Self); - - /** Type definition for the input image. */ - typedef TInputImage InputImageType; - typedef typename TInputImage::Pointer InputImagePointer; - typedef typename TInputImage::ConstPointer InputImageConstPointer; - - /** Type definition for the input image pixel type. */ - typedef typename TInputImage::PixelType InputImagePixelType; - - /** Type definition for the input image region type. */ - typedef typename TInputImage::RegionType InputImageRegionType; - - /** Type definition for the input image region iterator */ - typedef itk::ImageRegionIterator<TInputImage> InputImageRegionIterator; - typedef itk::ImageRegionConstIterator<TInputImage> InputImageRegionConstIterator; - - /** Image dimension */ - itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension); - - /** Type definitions for the training image. */ - typedef TClassifiedImage TrainingImageType; - typedef typename TClassifiedImage::Pointer TrainingImagePointer; - - /** Type definitions for the training image pixel type. */ - typedef typename TClassifiedImage::PixelType TrainingImagePixelType; - - /** Type definitions for the labelled image. - * It is derived from the training image. */ - typedef TClassifiedImage LabelledImageType; - typedef typename TClassifiedImage::Pointer LabelledImagePointer; - /** Type definitions for the classified image pixel type. - * It has to be the same type as the training image. */ - typedef typename TClassifiedImage::PixelType LabelledImagePixelType; - - /** Type definitions for the classified image pixel type. - * It has to be the same type as the training image. */ - typedef typename TClassifiedImage::RegionType LabelledImageRegionType; - - /** Type definition for the classified image index type. */ - typedef typename TClassifiedImage::IndexType LabelledImageIndexType; - typedef typename LabelledImageIndexType::IndexValueType IndexValueType; - - /** Type definition for the classified image offset type. */ - typedef typename TClassifiedImage::OffsetType LabelledImageOffsetType; - - /** Type definition for the input image region iterator */ - typedef itk::ImageRegionIterator<TClassifiedImage> LabelledImageRegionIterator; - - /** Labelled Image dimension */ - itkStaticConstMacro(ClassifiedImageDimension, unsigned int, TClassifiedImage::ImageDimension); - - /** Type definitions for classifier to be used for the MRF lavbelling. */ - typedef itk::ImageClassifierBase<TInputImage,TClassifiedImage> ClassifierType; - - /** Size and value typedef support. */ - typedef typename TInputImage::SizeType SizeType; - - /** Radius typedef support. */ - typedef typename TInputImage::SizeType NeighborhoodRadiusType; - - /** Input image neighborhood iterator and kernel size typedef */ - typedef itk::ConstNeighborhoodIterator< TInputImage > InputImageNeighborhoodIterator; - typedef typename InputImageNeighborhoodIterator::RadiusType InputImageNeighborhoodRadiusType; - typedef itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< TInputImage > InputImageFacesCalculator; - typedef typename InputImageFacesCalculator::FaceListType InputImageFaceListType; - typedef typename InputImageFaceListType::iterator InputImageFaceListIterator; - - /** Labelled image neighborhood interator typedef */ - typedef itk::NeighborhoodIterator< TClassifiedImage > LabelledImageNeighborhoodIterator; - typedef typename LabelledImageNeighborhoodIterator::RadiusType LabelledImageNeighborhoodRadiusType; - typedef itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< TClassifiedImage > LabelledImageFacesCalculator; - typedef typename LabelledImageFacesCalculator::FaceListType LabelledImageFaceListType; - typedef typename LabelledImageFaceListType::iterator LabelledImageFaceListIterator; - - /** Set the pointer to the classifer being used. */ - void SetClassifier( typename ClassifierType::Pointer ptrToClassifier ); - - - /** Set pipeline elements */ - // typedef GaussianMRFEnergy< TInputImage, TClassifiedImage> EnergyType; - typedef MRFEnergy< TClassifiedImage, TClassifiedImage> EnergyRegularizationType; - typedef MRFEnergy< TInputImage, TClassifiedImage> EnergyFidelityType; - typedef typename EnergyRegularizationType::Pointer EnergyRegularizationPointer; - typedef typename EnergyFidelityType::Pointer EnergyFidelityPointer; - - // typedef RandomMRFSampler< TInputImage, TClassifiedImage> SamplerType; - // typedef MRFSamplerMAP< TInputImage, TClassifiedImage> SamplerType; - typedef MRFSampler< TInputImage, TClassifiedImage> SamplerType; - typedef typename SamplerType::Pointer SamplerPointer; - - // typedef MetropolisMRFOptimizer OptimizerType; - // typedef MRFOptimizerICM OptimizerType; - typedef MRFOptimizer OptimizerType; - typedef typename OptimizerType::Pointer OptimizerPointer; - + { + public: + /** Standard class typedefs. */ + typedef MarkovRandomFieldFilter Self; + typedef itk::ImageToImageFilter<TInputImage,TClassifiedImage> Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + typedef typename Superclass::OutputImagePointer OutputImagePointer; + + /** Run-time type information (and related methods). */ + itkTypeMacro(MarkovRandomFieldFilter,itk::ImageToImageFilter); + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + + /** Type definition for the input image. */ + typedef TInputImage InputImageType; + typedef typename TInputImage::Pointer InputImagePointer; + typedef typename TInputImage::ConstPointer InputImageConstPointer; + + /** Type definition for the input image pixel type. */ + typedef typename TInputImage::PixelType InputImagePixelType; + + /** Type definition for the input image region type. */ + typedef typename TInputImage::RegionType InputImageRegionType; + + /** Type definition for the input image region iterator */ + typedef itk::ImageRegionIterator<TInputImage> InputImageRegionIterator; + typedef itk::ImageRegionConstIterator<TInputImage> InputImageRegionConstIterator; + + /** Image dimension */ + itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension); + + /** Type definitions for the training image. */ + typedef TClassifiedImage TrainingImageType; + typedef typename TClassifiedImage::Pointer TrainingImagePointer; + + /** Type definitions for the training image pixel type. */ + typedef typename TClassifiedImage::PixelType TrainingImagePixelType; + + /** Type definitions for the labelled image. + * It is derived from the training image. */ + typedef TClassifiedImage LabelledImageType; + typedef typename TClassifiedImage::Pointer LabelledImagePointer; + + /** Type definitions for the classified image pixel type. + * It has to be the same type as the training image. */ + typedef typename TClassifiedImage::PixelType LabelledImagePixelType; + + /** Type definitions for the classified image pixel type. + * It has to be the same type as the training image. */ + typedef typename TClassifiedImage::RegionType LabelledImageRegionType; + + /** Type definition for the classified image index type. */ + typedef typename TClassifiedImage::IndexType LabelledImageIndexType; + typedef typename LabelledImageIndexType::IndexValueType IndexValueType; - /* + /** Type definition for the classified image offset type. */ + typedef typename TClassifiedImage::OffsetType LabelledImageOffsetType; + + /** Type definition for the input image region iterator */ + typedef itk::ImageRegionIterator<TClassifiedImage> + LabelledImageRegionIterator; + + /** Labelled Image dimension */ + itkStaticConstMacro(ClassifiedImageDimension, unsigned int, + TClassifiedImage::ImageDimension); + + /** Type definitions for classifier to be used for the MRF lavbelling. */ + typedef itk::ImageClassifierBase<TInputImage,TClassifiedImage> ClassifierType; + + /** Size and value typedef support. */ + typedef typename TInputImage::SizeType SizeType; + + /** Radius typedef support. */ + typedef typename TInputImage::SizeType NeighborhoodRadiusType; + + /** Input image neighborhood iterator and kernel size typedef */ + typedef itk::ConstNeighborhoodIterator< TInputImage > + InputImageNeighborhoodIterator; + + + typedef typename InputImageNeighborhoodIterator::RadiusType + InputImageNeighborhoodRadiusType; + + typedef itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< TInputImage > + InputImageFacesCalculator; + + typedef typename InputImageFacesCalculator::FaceListType + InputImageFaceListType; + + typedef typename InputImageFaceListType::iterator + InputImageFaceListIterator; + + /** Labelled image neighborhood interator typedef */ + typedef itk::NeighborhoodIterator< TClassifiedImage > + LabelledImageNeighborhoodIterator; + + typedef typename LabelledImageNeighborhoodIterator::RadiusType + LabelledImageNeighborhoodRadiusType; + + typedef itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< TClassifiedImage > + LabelledImageFacesCalculator; + + typedef typename LabelledImageFacesCalculator::FaceListType + LabelledImageFaceListType; + + typedef typename LabelledImageFaceListType::iterator + LabelledImageFaceListIterator; + + /** Set the pointer to the classifer being used. */ + void SetClassifier( typename ClassifierType::Pointer ptrToClassifier ); + + + /** Set pipeline elements */ + typedef MRFEnergy< TClassifiedImage, TClassifiedImage> EnergyRegularizationType; + typedef MRFEnergy< TInputImage, TClassifiedImage> EnergyFidelityType; + + typedef typename EnergyRegularizationType::Pointer EnergyRegularizationPointer; + typedef typename EnergyFidelityType::Pointer EnergyFidelityPointer; + + typedef MRFSampler< TInputImage, TClassifiedImage> SamplerType; + typedef typename SamplerType::Pointer SamplerPointer; + + typedef MRFOptimizer OptimizerType; + typedef typename OptimizerType::Pointer OptimizerPointer; + + + /** ************ ACCESSORS ************ */ + + void SetEnergyRegularization( EnergyRegularizationPointer enerReg ){ m_EnergyRegularization=enerReg; this->Modified(); }; + EnergyRegularizationPointer GetEnergyRegularization(){ return m_EnergyRegularization; }; + void SetEnergyFidelity( EnergyFidelityPointer enerFid ){ m_EnergyFidelity=enerFid; this->Modified(); }; + EnergyFidelityPointer GetEnergyFidelity(){ return m_EnergyFidelity; }; + /* itkSetObjectMacro( EnergyRegularization, EnergyRegularizationType); itkGetObjectMacro( EnergyRegularization, EnergyRegularizationType); itkSetObjectMacro( EnergyFidelity, EnergyFidelityType); itkGetObjectMacro( EnergyFidelity, EnergyFidelityType); - - itkSetObjectMacro( Sampler, SamplerType); - itkGetObjectMacro( Sampler, SamplerType); - - itkSetObjectMacro( Optimizer, OptimizerType); - itkGetObjectMacro( Optimizer, OptimizerType); - - /** Set/Get the number of classes. */ - itkSetMacro(NumberOfClasses, unsigned int); - itkGetMacro(NumberOfClasses, unsigned int); - - /** Set/Get the number of iteration of the Iterated Conditional Mode - * (ICM) algorithm. A default value is set at 50 iterations. */ - itkSetMacro(MaximumNumberOfIterations, unsigned int); - itkGetMacro(MaximumNumberOfIterations, unsigned int); - - /** Set/Get the error tollerance level which is used as a threshold - * to quit the iterations */ - itkSetMacro(ErrorTolerance, double); - itkGetMacro(ErrorTolerance, double); - - /** Set/Get the degree of smoothing desired - * */ - itkSetMacro(SmoothingFactor, double); - itkGetMacro(SmoothingFactor, double); - - /** Set/Get the temperature - * */ - // itkSetMacro(Temperature, double); - // itkGetMacro(Temperature, double); - - /** Set/Get the regularization coefficient - * */ - itkSetMacro(Lambda, double); - itkGetMacro(Lambda, double); - - - /** Set the neighborhood radius */ - void SetNeighborhoodRadius(const NeighborhoodRadiusType &); - - /** Sets the radius for the neighborhood, calculates size from the - * radius, and allocates storage. */ - - void SetNeighborhoodRadius( const unsigned long ); - void SetNeighborhoodRadius( const unsigned long *radiusArray ); - - /** Get the neighborhood radius */ - const NeighborhoodRadiusType GetNeighborhoodRadius() const - { - NeighborhoodRadiusType m_NeighborhoodRadius; - - for(int i=0; i<InputImageDimension; ++i) - m_NeighborhoodRadius[i] = m_InputImageNeighborhoodRadius[i]; - - return m_NeighborhoodRadius; - } + */ + + + itkSetObjectMacro( Sampler, SamplerType); + itkGetObjectMacro( Sampler, SamplerType); + + itkSetObjectMacro( Optimizer, OptimizerType); + itkGetObjectMacro( Optimizer, OptimizerType); + + + + /** Set/Get the number of classes. */ + itkSetMacro(NumberOfClasses, unsigned int); + itkGetMacro(NumberOfClasses, unsigned int); + + /** Set/Get the number of iteration of the Iterated Conditional Mode + * (ICM) algorithm. A default value is set at 50 iterations. */ + itkSetMacro(MaximumNumberOfIterations, unsigned int); + itkGetMacro(MaximumNumberOfIterations, unsigned int); + + /** Set/Get the error tollerance level which is used as a threshold + * to quit the iterations */ + itkSetMacro(ErrorTolerance, double); + itkGetMacro(ErrorTolerance, double); + + /** Set/Get the degree of smoothing desired + * */ + itkSetMacro(SmoothingFactor, double); + itkGetMacro(SmoothingFactor, double); + + /** Set/Get the regularization coefficient + * */ + itkSetMacro(Lambda, double); + itkGetMacro(Lambda, double); + + + /** Set the neighborhood radius */ + void SetNeighborhoodRadius(const NeighborhoodRadiusType &); + + /** Sets the radius for the neighborhood, calculates size from the + * radius, and allocates storage. */ + + void SetNeighborhoodRadius( const unsigned long ); + void SetNeighborhoodRadius( const unsigned long *radiusArray ); + + /** Get the neighborhood radius */ + const NeighborhoodRadiusType GetNeighborhoodRadius() const + { + NeighborhoodRadiusType m_NeighborhoodRadius; + + for(int i=0; i<InputImageDimension; ++i) + m_NeighborhoodRadius[i] = m_InputImageNeighborhoodRadius[i]; + + return m_NeighborhoodRadius; + } + + + /** Set training image for the starting point. This is not compulsory: + * if the starting image is not specified, a random image will be used + * instead. + */ + virtual void SetTrainingInput( const TrainingImageType * trainingImage); + TrainingImageType* GetTrainingInput(void); + + + //Enum to get the stopping condition of the MRF filter + typedef enum{ + MaximumNumberOfIterations=1, + ErrorTolerance + } StopConditionType; + + /** Get condition that stops the MRF filter (Number of Iterations + * / Error tolerance ) */ + itkGetConstReferenceMacro( StopCondition, StopConditionType ); + + /** Get macro for number of iterations */ + itkGetConstReferenceMacro( NumberOfIterations, unsigned int ); - - /** Set training image */ - // virtual void SetTrainingImage (const TrainingImagePointer _arg) - // { - // itkDebugMacro("setting TrainingImage to " << _arg); - // if (this->m_TrainingImage != _arg) - // { - // // this->SetInput(1,_arg); - // // this->m_TrainingImage = this->GetInput(1); - // this->m_TrainingImage = _arg; - // this->Modified(); - // } - // m_ExternalClassificationSet = true; - // } - // itkSetMacro(TrainingImage,TrainingImagePointer); - // itkGetMacro(TrainingImage,TrainingImagePointer); - // itkSetInputMacro(TrainingImage,TrainingImagePointer,1); - // itkGetInputMacro(TrainingImage,TrainingImagePointer,1); - virtual void SetTrainingInput( const TrainingImageType * trainningImage); - TrainingImageType* GetTrainingInput(void); - - /** Store a value to be used instead of random value.. FOR TEST ONLY*/ - void SetValueInsteadRandom( int val ) - { - m_ValueInsteadRandom = val; - std::cout<<"The m_ValueInsteadRandom varaible has to be used only for tests..."<<std::endl; - this->Modified(); - }; - - //Enum to get the stopping condition of the MRF filter - typedef enum{ MaximumNumberOfIterations=1, ErrorTolerance } StopConditionType; - - /** Get condition that stops the MRF filter (Number of Iterations - * / Error tolerance ) */ - itkGetConstReferenceMacro( StopCondition, StopConditionType ); - - /* Get macro for number of iterations */ - itkGetConstReferenceMacro( NumberOfIterations, unsigned int ); - #ifdef ITK_USE_CONCEPT_CHECKING - /** Begin concept checking */ - itkConceptMacro(UnsignedIntConvertibleToClassifiedCheck, (itk::Concept::Convertible<unsigned int, LabelledImagePixelType>)); - itkConceptMacro(ClassifiedConvertibleToUnsignedIntCheck, (itk::Concept::Convertible<LabelledImagePixelType, unsigned int> )); - itkConceptMacro(ClassifiedConvertibleToIntCheck, (itk::Concept::Convertible<LabelledImagePixelType, int> )); - itkConceptMacro(IntConvertibleToClassifiedCheck, (itk::Concept::Convertible<int, LabelledImagePixelType>)); - itkConceptMacro(SameDimensionCheck, (itk::Concept::SameDimension<InputImageDimension, ClassifiedImageDimension>)); - /** End concept checking */ + /** Begin concept checking */ + itkConceptMacro(UnsignedIntConvertibleToClassifiedCheck, + (itk::Concept::Convertible<unsigned int, LabelledImagePixelType>)); + itkConceptMacro(ClassifiedConvertibleToUnsignedIntCheck, + (itk::Concept::Convertible<LabelledImagePixelType, unsigned int> )); + itkConceptMacro(ClassifiedConvertibleToIntCheck, + (itk::Concept::Convertible<LabelledImagePixelType, int> )); + itkConceptMacro(IntConvertibleToClassifiedCheck, + (itk::Concept::Convertible<int, LabelledImagePixelType>)); + itkConceptMacro(SameDimensionCheck, + (itk::Concept::SameDimension<InputImageDimension, ClassifiedImageDimension>)); + /** End concept checking */ #endif - - - protected: - MarkovRandomFieldFilter(); - ~MarkovRandomFieldFilter(); - void PrintSelf(std::ostream& os, itk::Indent indent) const; - - /** Allocate memory for labelled images. */ - void Allocate(); - - void Initialize() throw (itk::ExceptionObject); - - virtual void ApplyMarkovClassificationFilter(); - virtual void GenerateData(); - virtual void GenerateInputRequestedRegion(); - virtual void EnlargeOutputRequestedRegion( itk::DataObject * ); - virtual void GenerateOutputInformation(); - - - MarkovRandomFieldFilter(const Self&); //purposely not implemented - void operator=(const Self&); //purposely not implemented - - typedef typename TInputImage::SizeType InputImageSizeType; - - - InputImageNeighborhoodRadiusType m_InputImageNeighborhoodRadius; - LabelledImageNeighborhoodRadiusType m_LabelledImageNeighborhoodRadius; - - unsigned int m_NumberOfClasses; - unsigned int m_MaximumNumberOfIterations; - unsigned int m_KernelSize; - unsigned int m_NumberOfIterations; - int m_ErrorCounter; - // int m_NeighborhoodSize; - int m_NeighborhoodRadius; - int m_TotalNumberOfValidPixelsInOutputImage; - int m_TotalNumberOfPixelsInInputImage; - double m_ErrorTolerance; - double m_SmoothingFactor; - double m_Lambda; - // double *m_ClassProbability; //Class liklihood - // double m_Temperature; - bool m_ExternalClassificationSet; - StopConditionType m_StopCondition; - std::vector<double> m_MRFNeighborhoodWeight; - std::vector<double> m_NeighborInfluence; - std::vector<double> m_DummyVector; - - - /** Pointer to different elements */ - // TrainingImagePointer m_TrainingImage; - EnergyRegularizationPointer m_EnergyRegularization; - EnergyFidelityPointer m_EnergyFidelity; - OptimizerPointer m_Optimizer; - SamplerPointer m_Sampler; - - virtual void MinimizeOnce(); - - private: - /** Store a value to be used instead of random value.. FOR TEST ONLY*/ - int m_ValueInsteadRandom; - -}; // class MarkovRandomFieldFilter + + + protected: + MarkovRandomFieldFilter(); + ~MarkovRandomFieldFilter(); + void PrintSelf(std::ostream& os, itk::Indent indent) const; + + /** Allocate memory for labelled images. This is automatically called + * in GenerateData(). + */ + void Allocate(); + + /** Connect the pipeline and propagate the required parameters. This + * is automatically called in GenerateData(). + */ + void Initialize() throw (itk::ExceptionObject); + + + virtual void ApplyMarkovRandomFieldFilter(); + + + virtual void GenerateData(); + virtual void GenerateInputRequestedRegion(); + virtual void EnlargeOutputRequestedRegion( itk::DataObject * ); + virtual void GenerateOutputInformation(); + + + MarkovRandomFieldFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + typedef typename TInputImage::SizeType InputImageSizeType; + + + InputImageNeighborhoodRadiusType m_InputImageNeighborhoodRadius; + LabelledImageNeighborhoodRadiusType m_LabelledImageNeighborhoodRadius; + + unsigned int m_NumberOfClasses; + unsigned int m_MaximumNumberOfIterations; + unsigned int m_KernelSize; + + int m_ErrorCounter; + double m_ImageDeltaEnergy; + + int m_NeighborhoodRadius; + int m_TotalNumberOfValidPixelsInOutputImage; + int m_TotalNumberOfPixelsInInputImage; + double m_ErrorTolerance; + double m_SmoothingFactor; + + unsigned int m_NumberOfIterations; + + + double m_Lambda; + bool m_ExternalClassificationSet; + StopConditionType m_StopCondition; + + TrainingImagePointer m_TrainingImage; + + std::vector<double> m_MRFNeighborhoodWeight; + std::vector<double> m_NeighborInfluence; + std::vector<double> m_DummyVector; + + + /** Pointer to different elements */ + + EnergyRegularizationPointer m_EnergyRegularization; + EnergyFidelityPointer m_EnergyFidelity; + OptimizerPointer m_Optimizer; + SamplerPointer m_Sampler; + + virtual void MinimizeOnce(); + + private: + + }; // class MarkovRandomFieldFilter -} // namespace itk +} // namespace otb #ifndef ITK_MANUAL_INSTANTIATION #include "otbMarkovRandomFieldFilter.txx" diff --git a/Code/Markov/otbMarkovRandomFieldFilter.txx b/Code/Markov/otbMarkovRandomFieldFilter.txx index 7f72ec4958..e5173ff718 100644 --- a/Code/Markov/otbMarkovRandomFieldFilter.txx +++ b/Code/Markov/otbMarkovRandomFieldFilter.txx @@ -22,27 +22,38 @@ namespace otb { - +/** ATTENTION : j'ai passe les initialisation dans le constructeur car certaines plateformes (dont NEMO) n'aiment cette maniere d'initialiser*/ template<class TInputImage, class TClassifiedImage> MarkovRandomFieldFilter<TInputImage,TClassifiedImage> -::MarkovRandomFieldFilter(void): +::MarkovRandomFieldFilter(void) +/* m_NumberOfClasses(0), m_MaximumNumberOfIterations(50), m_ErrorCounter(0), + m_ImageDeltaEnergy(0.0), m_NeighborhoodRadius(1), m_TotalNumberOfValidPixelsInOutputImage(1), m_TotalNumberOfPixelsInInputImage(1), - //m_ErrorTolerance(0.2), - //m_SmoothingFactor(1), - //m_NumberOfIterations(0), - m_ExternalClassificationSet(false), - m_StopCondition(MaximumNumberOfIterations) + m_ErrorTolerance(0.0), + m_SmoothingFactor(1), + m_NumberOfIterations(0), + m_StopCondition(MaximumNumberOfIterations), + m_ExternalClassificationSet(false) + */ { - m_ExternalClassificationSet = false; + m_NumberOfClasses = 0; + m_MaximumNumberOfIterations = 50; + m_ErrorCounter = 0; + m_ImageDeltaEnergy = 0.0; + m_NeighborhoodRadius = 1; + m_TotalNumberOfValidPixelsInOutputImage = 1; + m_TotalNumberOfPixelsInInputImage = 1; + m_ErrorTolerance = 0.0; m_SmoothingFactor = 1; m_NumberOfIterations = 0; - m_ErrorTolerance = 0.2; - this->SetNumberOfRequiredInputs(2); + m_StopCondition = MaximumNumberOfIterations; + m_ExternalClassificationSet = false; + this->SetNumberOfRequiredInputs(1); if( (int)InputImageDimension != (int)ClassifiedImageDimension ) { itk::OStringStream msg; @@ -50,18 +61,14 @@ MarkovRandomFieldFilter<TInputImage,TClassifiedImage> throw itk::ExceptionObject(__FILE__, __LINE__,msg.str().c_str(),ITK_LOCATION); } m_InputImageNeighborhoodRadius.Fill(m_NeighborhoodRadius); - // m_MRFNeighborhoodWeight.resize(0); - // m_NeighborInfluence.resize(0); - // m_DummyVector.resize(0); - // this->SetMRFNeighborhoodWeight( m_DummyVector ); - // this->SetDefaultMRFNeighborhoodWeight(); - - EnergyRegularizationPointer m_EnergyRegularization = EnergyRegularizationType::New();; - EnergyFidelityPointer m_EnergyFidelity = EnergyFidelityType::New();; - OptimizerPointer m_Optimizer = OptimizerType::New();; - SamplerPointer m_Sampler = SamplerType::New(); - - m_ValueInsteadRandom = itk::NumericTraits<int>::min(); + // m_MRFNeighborhoodWeight.resize(0); + // m_NeighborInfluence.resize(0); + // m_DummyVector.resize(0); + // this->SetMRFNeighborhoodWeight( m_DummyVector ); + // this->SetDefaultMRFNeighborhoodWeight(); + + srand((unsigned)time(0)); + } template<class TInputImage, class TClassifiedImage> @@ -76,6 +83,7 @@ void MarkovRandomFieldFilter<TInputImage, TClassifiedImage> ::SetTrainingInput(const TrainingImageType * trainingImage) { + m_ExternalClassificationSet = true; // Process object is not const-correct so the const_cast is required here this->itk::ProcessObject::SetNthInput(1, const_cast< TrainingImageType * >( trainingImage ) ); this->Modified(); @@ -91,7 +99,7 @@ MarkovRandomFieldFilter<TInputImage, TClassifiedImage> { return 0; } - return static_cast<TrainingImageType * > + return static_cast<const TrainingImageType * > (this->itk::ProcessObject::GetInput(1) ); } @@ -102,35 +110,51 @@ MarkovRandomFieldFilter<TInputImage, TClassifiedImage> ::PrintSelf( std::ostream& os, itk::Indent indent ) const { Superclass::PrintSelf(os,indent); + os << indent <<" MRF Image filter object " << std::endl; + os << indent <<" Number of classes: " << m_NumberOfClasses << std::endl; - os << indent <<" Maximum number of iterations: " << m_MaximumNumberOfIterations << std::endl; - os << indent <<" Error tolerance for convergence: " << m_ErrorTolerance << std::endl; - os << indent <<" Size of the MRF neighborhood radius:" << m_InputImageNeighborhoodRadius << std::endl; - os << indent << "StopCondition: " << m_StopCondition << std::endl; - os << indent <<" Number of iterations: " << m_NumberOfIterations << std::endl; + os << indent <<" Maximum number of iterations: " << + m_MaximumNumberOfIterations << std::endl; + + os << indent <<" Error tolerance for convergence: " << + m_ErrorTolerance << std::endl; + + os << indent <<" Size of the MRF neighborhood radius:" << + m_InputImageNeighborhoodRadius << std::endl; + + os << indent << "StopCondition: " + << m_StopCondition << std::endl; + + os << indent <<" Number of iterations: " << + m_NumberOfIterations << std::endl; + + os << indent <<" Lambda: " << + m_Lambda << std::endl; }// end PrintSelf -/* +/** * GenerateInputRequestedRegion method. */ template <class TInputImage, class TClassifiedImage> void MarkovRandomFieldFilter<TInputImage, TClassifiedImage> ::GenerateInputRequestedRegion() -{ +{ + // this filter requires the all of the input images // to be at the size of the output requested region - InputImagePointer inputPtr = const_cast< InputImageType * >( this->GetInput() ); + InputImagePointer inputPtr = + const_cast< InputImageType * >( this->GetInput() ); OutputImagePointer outputPtr = this->GetOutput(); inputPtr->SetRequestedRegion( outputPtr->GetRequestedRegion() ); } -/* +/** * EnlargeOutputRequestedRegion method. */ template <class TInputImage, class TClassifiedImage> @@ -147,7 +171,7 @@ MarkovRandomFieldFilter<TInputImage, TClassifiedImage> -/* +/** * GenerateOutputInformation method. */ template <class TInputImage, class TClassifiedImage> @@ -155,7 +179,7 @@ void MarkovRandomFieldFilter<TInputImage, TClassifiedImage> ::GenerateOutputInformation() { - typename TInputImage::ConstPointer input = this->GetInput(); + typename TInputImage::ConstPointer input = this->GetInput(); typename TClassifiedImage::Pointer output = this->GetOutput(); output->SetLargestPossibleRegion( input->GetLargestPossibleRegion() ); } @@ -168,11 +192,8 @@ void MarkovRandomFieldFilter<TInputImage, TClassifiedImage> ::GenerateData() { - //First run the Gaussian classifier calculator and - //generate the Gaussian model for the different classes - //and then generate the initial labelled dataset. - InputImageConstPointer inputImage = this->GetInput(); +// InputImageConstPointer inputImage = this->GetInput(); //Allocate memory for the labelled images this->Allocate(); @@ -180,15 +201,17 @@ MarkovRandomFieldFilter<TInputImage, TClassifiedImage> //Branch the pipeline this->Initialize(); - this->ApplyMarkovClassificationFilter(); + //Run the Markov random field + this->ApplyMarkovRandomFieldFilter(); + }// end GenerateData -//------------------------------------------------------- -//Set the neighborhood radius -//------------------------------------------------------- +/** +* Set the neighborhood radius from a single value +*/ template<class TInputImage, class TClassifiedImage> void MarkovRandomFieldFilter<TInputImage, TClassifiedImage> @@ -205,7 +228,9 @@ MarkovRandomFieldFilter<TInputImage, TClassifiedImage> }// end SetNeighborhoodRadius - +/** +* Set the neighborhood radius from an array +*/ template<class TInputImage, class TClassifiedImage> void MarkovRandomFieldFilter<TInputImage, TClassifiedImage> @@ -223,8 +248,9 @@ MarkovRandomFieldFilter<TInputImage, TClassifiedImage> - -//Set the neighborhood radius +/** +* Set the neighborhood radius from a radius object +*/ template<class TInputImage, class TClassifiedImage> void MarkovRandomFieldFilter<TInputImage, TClassifiedImage> @@ -242,10 +268,9 @@ MarkovRandomFieldFilter<TInputImage, TClassifiedImage> -//------------------------------------------------------- -//------------------------------------------------------- -//Allocate algorithm specific resources -//------------------------------------------------------- +/** +* Allocate algorithm specific resources +*/ template<class TInputImage, class TClassifiedImage> void MarkovRandomFieldFilter<TInputImage, TClassifiedImage> @@ -256,19 +281,7 @@ MarkovRandomFieldFilter<TInputImage, TClassifiedImage> throw itk::ExceptionObject(__FILE__, __LINE__,"NumberOfClasses <= 0.",ITK_LOCATION); } - InputImageSizeType inputImageSize = this->GetInput()->GetBufferedRegion().GetSize(); - - //--------------------------------------------------------------------- - //Get the number of valid pixels in the output MRF image - //--------------------------------------------------------------------- - int tmp; - for( unsigned int i=0; i < InputImageDimension; i++ ) - { - tmp = static_cast<int>(inputImageSize[i]); - m_TotalNumberOfPixelsInInputImage *= tmp; - m_TotalNumberOfValidPixelsInOutputImage *= ( tmp - 2*m_InputImageNeighborhoodRadius[i] ); - } - + //Set the output labelled and allocate the memory LabelledImagePointer outputPtr = this->GetOutput(); @@ -278,17 +291,18 @@ MarkovRandomFieldFilter<TInputImage, TClassifiedImage> //Copy input data in the output buffer memory or //initialize to random values if not set - LabelledImageRegionIterator outImageIt( outputPtr, outputPtr->GetRequestedRegion() ); - - m_ExternalClassificationSet = true; //TODO switch to random if the ref image is not specified + LabelledImageRegionIterator + outImageIt( outputPtr, outputPtr->GetRequestedRegion() ); + if (m_ExternalClassificationSet) { TrainingImagePointer trainingImage = this->GetTrainingInput(); - LabelledImageRegionIterator trainingImageIt( trainingImage, outputPtr->GetRequestedRegion() ); + LabelledImageRegionIterator + trainingImageIt( trainingImage, outputPtr->GetRequestedRegion() ); while ( !outImageIt.IsAtEnd() ) { - LabelledImagePixelType labelvalue = ( LabelledImagePixelType ) trainingImageIt.Get(); + LabelledImagePixelType labelvalue = static_cast<LabelledImagePixelType> (trainingImageIt.Get()); outImageIt.Set( labelvalue ); ++trainingImageIt; @@ -297,42 +311,53 @@ MarkovRandomFieldFilter<TInputImage, TClassifiedImage> } else //set to random value { - // if it is a test, cancel the rand() - if( m_ValueInsteadRandom == itk::NumericTraits<int>::min() ) - { - srand((unsigned)time(0)); - - while ( !outImageIt.IsAtEnd() ) - { - LabelledImagePixelType randomvalue = static_cast<LabelledImagePixelType>(rand() % static_cast<int>(m_NumberOfClasses)); - outImageIt.Set( randomvalue ); - ++outImageIt; - }// end while - } - else +// srand((unsigned)time(0)); + + while ( !outImageIt.IsAtEnd() ) { - while ( !outImageIt.IsAtEnd() ) - { - LabelledImagePixelType randomvalue = static_cast<LabelledImagePixelType>(m_ValueInsteadRandom % static_cast<int>(m_NumberOfClasses)); - outImageIt.Set( randomvalue ); - ++outImageIt; - } - } + LabelledImagePixelType randomvalue = static_cast<LabelledImagePixelType>( + rand() % static_cast<int>(m_NumberOfClasses) + ); + outImageIt.Set( randomvalue ); + ++outImageIt; + }// end while } }// Allocate - +/** +* Initialize pipeline and values +*/ template<class TInputImage, class TClassifiedImage> void MarkovRandomFieldFilter<TInputImage, TClassifiedImage> ::Initialize() throw (itk::ExceptionObject) -{ - if( m_ValueInsteadRandom == itk::NumericTraits<int>::min() ) +{ + + m_ImageDeltaEnergy=0.0; + + InputImageSizeType inputImageSize = + this->GetInput()->GetBufferedRegion().GetSize(); + + //--------------------------------------------------------------------- + //Get the number of valid pixels in the output MRF image + //--------------------------------------------------------------------- + + m_TotalNumberOfPixelsInInputImage = 1; + m_TotalNumberOfValidPixelsInOutputImage = 1; + + for( unsigned int i=0; i < InputImageDimension; i++ ) { - srand((unsigned)time(0)); - } + m_TotalNumberOfPixelsInInputImage *= static_cast<int>(inputImageSize[i]); + + m_TotalNumberOfValidPixelsInOutputImage *= + ( static_cast<int>(inputImageSize[i]) + - 2*m_InputImageNeighborhoodRadius[i] ); + } + + + srand((unsigned)time(0)); if ( !m_EnergyRegularization ) { @@ -358,50 +383,53 @@ MarkovRandomFieldFilter<TInputImage, TClassifiedImage> m_Sampler->SetEnergyRegularization(m_EnergyRegularization); m_Sampler->SetEnergyFidelity(m_EnergyFidelity); m_Sampler->SetNumberOfClasses(m_NumberOfClasses); - // m_Optimizer->SetTemperature(m_Temperature); } -//------------------------------------------------------- -//------------------------------------------------------- -//Apply the MRF image filter -//------------------------------------------------------- +/** +*Apply the MRF image filter +*/ template<class TInputImage, class TClassifiedImage> void MarkovRandomFieldFilter<TInputImage, TClassifiedImage> -::ApplyMarkovClassificationFilter() +::ApplyMarkovRandomFieldFilter() { - InputImageSizeType inputImageSize = this->GetInput()->GetBufferedRegion().GetSize(); - int totalNumberOfPixelsInInputImage = 1; - for( unsigned int i = 0; i < InputImageDimension; i++ ) - { - totalNumberOfPixelsInInputImage *= static_cast<int>(inputImageSize[ i ]) ; - } - - int maxNumPixelError = (int) ( vnl_math_rnd (m_ErrorTolerance * m_TotalNumberOfValidPixelsInOutputImage) ); + + //Note: error should be defined according to the number of valid pixel in the output + int maxNumPixelError = (int) ( vnl_math_rnd (m_ErrorTolerance * + m_TotalNumberOfPixelsInInputImage) ); + m_NumberOfIterations = 0; - do + m_ErrorCounter = m_TotalNumberOfValidPixelsInOutputImage; + + while(( m_NumberOfIterations < m_MaximumNumberOfIterations ) && + ( m_ErrorCounter >= maxNumPixelError ) ) { - otbMsgDebugMacro(<< "Iteration No." << m_NumberOfIterations); + otbMsgDebugMacro(<< "Iteration No." << m_NumberOfIterations); +// std::cerr << "Iteration No." << m_NumberOfIterations << std::endl; + this->MinimizeOnce(); + otbMsgDebugMacro(<< "m_ErrorCounter/m_TotalNumberOfPixelsInInputImage: " << m_ErrorCounter/((double)(m_TotalNumberOfPixelsInInputImage))); +// std::cerr << "m_ErrorCounter/m_TotalNumberOfPixelsInInputImage: " +// << m_ErrorCounter/((double)(m_TotalNumberOfPixelsInInputImage)) +// << std::endl; +// std::cerr << "m_ImageDeltaEnergy: " << m_ImageDeltaEnergy << std::endl; - m_NumberOfIterations += 1; - // m_ErrorCounter = m_TotalNumberOfValidPixelsInOutputImage - - // totalNumberOfPixelsInInputImage; + ++m_NumberOfIterations; - } - while(( m_NumberOfIterations < m_MaximumNumberOfIterations ) && - ( m_ErrorCounter > maxNumPixelError ) ); + } otbMsgDebugMacro(<< "m_NumberOfIterations: " << m_NumberOfIterations); otbMsgDebugMacro(<< "m_MaximumNumberOfIterations: " << m_MaximumNumberOfIterations); otbMsgDebugMacro(<< "m_ErrorCounter: " << m_ErrorCounter); otbMsgDebugMacro(<< "maxNumPixelError: " << maxNumPixelError); - - +// std::cerr << "m_NumberOfIterations: " << m_NumberOfIterations << std::endl; +// std::cerr << "m_MaximumNumberOfIterations: " << m_MaximumNumberOfIterations << std::endl; +// std::cerr << "m_ErrorCounter: " << m_ErrorCounter << std::endl; +// std::cerr << "maxNumPixelError: " << maxNumPixelError << std::endl; //Determine stop condition if( m_NumberOfIterations >= m_MaximumNumberOfIterations ) @@ -413,34 +441,44 @@ MarkovRandomFieldFilter<TInputImage, TClassifiedImage> m_StopCondition = ErrorTolerance; } -}// ApplyMarkovClassificationFilter - +}// ApplyMarkovRandomFieldFilter +/** +*Apply the MRF image filter on the whole image once +*/ template<class TInputImage, class TClassifiedImage> void MarkovRandomFieldFilter<TInputImage, TClassifiedImage> ::MinimizeOnce() { - LabelledImageNeighborhoodIterator labelledIterator(m_LabelledImageNeighborhoodRadius, this->GetOutput(), this->GetOutput()->GetLargestPossibleRegion() ); - InputImageNeighborhoodIterator dataIterator(m_InputImageNeighborhoodRadius, this->GetInput(), this->GetInput()->GetLargestPossibleRegion() ); + LabelledImageNeighborhoodIterator + labelledIterator(m_LabelledImageNeighborhoodRadius, this->GetOutput(), + this->GetOutput()->GetLargestPossibleRegion() ); + InputImageNeighborhoodIterator + dataIterator(m_InputImageNeighborhoodRadius, this->GetInput(), + this->GetInput()->GetLargestPossibleRegion() ); m_ErrorCounter = 0; - //WARNING, is it the same size ?? - for (labelledIterator.GoToBegin(), dataIterator.GoToBegin(); !labelledIterator.IsAtEnd(); ++labelledIterator, ++dataIterator) - { - LabelledImagePixelType value; - bool changeValue; - m_Sampler->Compute(dataIterator,labelledIterator); - value=m_Sampler->GetValue(); - changeValue= m_Optimizer->Compute(m_Sampler->GetDeltaEnergy()); - if (changeValue) - { - labelledIterator.SetCenterPixel(value); - ++m_ErrorCounter; - } + + for (labelledIterator.GoToBegin(), dataIterator.GoToBegin(); + !labelledIterator.IsAtEnd(); + ++labelledIterator, ++dataIterator){ + + LabelledImagePixelType value; + bool changeValueBool; + m_Sampler->Compute(dataIterator,labelledIterator); + value=m_Sampler->GetValue(); + changeValueBool= m_Optimizer->Compute(m_Sampler->GetDeltaEnergy()); + if (changeValueBool){ + labelledIterator.SetCenterPixel(value); + ++m_ErrorCounter; + m_ImageDeltaEnergy += m_Sampler->GetDeltaEnergy(); } + } + } + } // namespace otb #endif -- GitLab