From 00b43e66efdbe42e0b95264e241780b00dff7a14 Mon Sep 17 00:00:00 2001
From: Emmanuel Christophe <emmanuel.christophe@orfeo-toolbox.org>
Date: Tue, 29 Apr 2008 16:02:35 +0000
Subject: [PATCH] Ajout d'un sampler et energy
---
Code/Markov/otbMRFEnergyEdgeFidelity.h | 54 ++++++
Code/Markov/otbMRFSamplerRandomMAP.h | 154 ++++++++++++++++++
.../Markov/MarkovClassification1Example.cxx | 2 +-
.../Markov/MarkovClassification2Example.cxx | 8 +-
4 files changed, 213 insertions(+), 5 deletions(-)
create mode 100644 Code/Markov/otbMRFEnergyEdgeFidelity.h
create mode 100644 Code/Markov/otbMRFSamplerRandomMAP.h
diff --git a/Code/Markov/otbMRFEnergyEdgeFidelity.h b/Code/Markov/otbMRFEnergyEdgeFidelity.h
new file mode 100644
index 0000000000..bf3501a736
--- /dev/null
+++ b/Code/Markov/otbMRFEnergyEdgeFidelity.h
@@ -0,0 +1,54 @@
+
+#ifndef _otbMRFEnergyEdgeFidelity_h
+#define _otbMRFEnergyEdgeFidelity_h
+
+#include "otbMRFEnergy.h"
+#define M_SQUARE(x) ((x)*(x))
+
+namespace otb
+{
+ /**
+ * \class MRFEnergyEdgeFidelity
+ * \brief This is the implementation of an edge preserving model for Markov denoising.
+ *
+ * This is the implementation of an edge fidelity model for Markov denoising, to be used for
+ * regularization. Energy is:
+ * \f[ \sum_{t \in \mathcal{V}_s} U(x_s,x_t) = \Phi(x_s-x_t) \f]
+ * with
+ * - \f$ x_s \f$ the label on site s
+ * - \f$ x_t \f$ the label on site t, a neighbor of s
+ * - \f$ \Phi \f$ an edge preserving function:
+ \f[ \Phi(u) = \frac{u^2}{1+u^2} \f]
+ */
+
+ template< class TInput1, class TInput2>
+ class ITK_EXPORT MRFEnergyEdgeFidelity : public MRFEnergy< TInput1, TInput2>
+ {
+ public:
+ typedef MRFEnergyEdgeFidelity Self;
+ typedef itk::Object Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ typedef itk::ConstNeighborhoodIterator< TInput1 > NeighborhoodIterator;
+ typedef typename TInput2::PixelType LabelledImagePixelType;
+
+ itkNewMacro(Self);
+
+// itkSetMacro(Beta, double);
+ itkTypeMacro(MRFEnergy,Object);
+
+ double GetSingleValue(const InputImagePixelType & value1, const LabelledImagePixelType & value2)
+ {
+ return M_SQUARE((value1 - value2))/double(1+M_SQUARE((value1 - value2)));
+ }
+
+
+ protected:
+ // The constructor and destructor.
+ MRFEnergyEdgeFidelity() {};
+ virtual ~MRFEnergyEdgeFidelity() {};
+ };
+}
+
+#endif
diff --git a/Code/Markov/otbMRFSamplerRandomMAP.h b/Code/Markov/otbMRFSamplerRandomMAP.h
new file mode 100644
index 0000000000..4b5aaa3803
--- /dev/null
+++ b/Code/Markov/otbMRFSamplerRandomMAP.h
@@ -0,0 +1,154 @@
+
+#ifndef _MRFSamplerRandomMAP_h
+#define _MRFSamplerRandomMAP_h
+
+// #include "otbMRFEnergy.h"
+#include "otbMRFSampler.h"
+
+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 itk::Object Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ typedef itk::ConstNeighborhoodIterator< TInput1 > InputImageNeighborhoodIterator;
+// typedef itk::NeighborhoodIterator< 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(MRFSamplerRandomMAP,Object);
+
+// void SetNumberOfClasses(unsigned int arg)
+// {
+// otbDebugMacro("setting NumberOfClasses to " << arg);
+// if (this->m_NumberOfClasses != arg)
+// {
+// this->m_NumberOfClasses = arg;
+// repartitionFunction = (double *) calloc(m_NumberOfClasses, sizeof(double));
+// 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)
+ {
+ this->m_EnergyBefore=this->m_EnergyFidelity->GetValue(itData, itRegul.GetCenterPixel())
+ + 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();
+// otbDebugMacro(<< "Computing MAP for pix " << itData.GetIndex());
+ //Compute probability for each possibility
+ double totalProba=0.0;
+ for (LabelledImagePixelType valueCurrent = 0;
+ valueCurrent < this->m_NumberOfClasses; ++valueCurrent)
+ {
+// otbDebugMacro(<< " --> Proposed value " << static_cast<double>(valueCurrent));
+ this->m_EnergyCurrent = this->m_EnergyFidelity->GetValue(itData, valueCurrent)
+ + this->m_Lambda * this->m_EnergyRegularization->GetValue(itRegul, valueCurrent);
+
+ energy[valueCurrent] = this->m_EnergyCurrent;
+ repartitionFunction[valueCurrent] = vcl_exp(-this->m_EnergyCurrent)+totalProba;
+ totalProba = repartitionFunction[valueCurrent];
+// otbDebugMacro("valueCurrent, repartitionFunction[valueCurrent] " << (unsigned int) valueCurrent << ", " << repartitionFunction[valueCurrent]);
+
+ }
+
+ //Pick a value according to probability
+
+ double select = (rand()/(double(RAND_MAX)+1) * totalProba);
+// otbDebugMacro("repartitionFunction " << repartitionFunction[0] << " "
+// << repartitionFunction[1] << " " << repartitionFunction[2] << " "
+// << repartitionFunction[3] << " ");
+
+// otbDebugMacro("select, totalProba " << select << ", " << totalProba);
+ LabelledImagePixelType valueCurrent = 0;
+ for (valueCurrent = 0;
+ valueCurrent < this->m_NumberOfClasses; ++valueCurrent)
+ {
+ if (repartitionFunction[valueCurrent] > select) break;
+ }
+
+ if ( this->m_Value != valueCurrent)
+ {
+ this->m_Value = valueCurrent;
+ this->m_EnergyAfter = energy[valueCurrent];
+ }
+
+ this->m_DeltaEnergy= this->m_EnergyAfter - this->m_EnergyBefore;
+// otbDebugMacro("Decision " << (unsigned int) valueCurrent);
+ return 0;
+ }
+
+
+ private:
+// unsigned int m_NumberOfClasses;
+// double m_EnergyBefore;
+// double m_EnergyAfter;
+// double m_DeltaEnergy;
+// LabelledImagePixelType m_Value;
+// EnergyRegularizationPointer m_EnergyRegularization;
+// EnergyFidelityPointer m_EnergyFidelity;
+// LabelledImagePixelType valueCurrent;
+// double m_EnergyCurrent;
+ double * repartitionFunction ;
+ double * energy;
+// double m_Lambda;
+
+
+ protected:
+ // The constructor and destructor.
+ MRFSamplerRandomMAP() {
+ energy = (double *) calloc(this->m_NumberOfClasses, sizeof(double));
+ repartitionFunction = (double *) calloc(this->m_NumberOfClasses, sizeof(double));
+ }
+ virtual ~MRFSamplerRandomMAP() {}
+
+ };
+
+
+}
+
+#endif
diff --git a/Examples/Markov/MarkovClassification1Example.cxx b/Examples/Markov/MarkovClassification1Example.cxx
index 497b880d8a..46c5167940 100644
--- a/Examples/Markov/MarkovClassification1Example.cxx
+++ b/Examples/Markov/MarkovClassification1Example.cxx
@@ -39,7 +39,7 @@
//
// This example applies the \doxygen{otb}{MarkovClassificationFilter} to
// classify an image into four classes defined by their mean and variance. The
-// optimization is done using an Metropolis algorithm with a MAP estimator. The
+// optimization is done using an Metropolis algorithm with a random sampler. The
// regularization energy is defined by a Potts model and the fidelity by a
// Gaussian model.
//
diff --git a/Examples/Markov/MarkovClassification2Example.cxx b/Examples/Markov/MarkovClassification2Example.cxx
index 7c25c723a9..0d5bf64e78 100644
--- a/Examples/Markov/MarkovClassification2Example.cxx
+++ b/Examples/Markov/MarkovClassification2Example.cxx
@@ -39,7 +39,7 @@
//
// This example applies the \doxygen{otb}{MarkovClassificationFilter} to
// classify an image into four classes defined by their mean and variance. The
-// optimization is done using an Metropolis algorithm with a MAP estimator. The
+// optimization is done using an ICM algorithm with a MAP estimator. The
// regularization energy is defined by a Potts model and the fidelity by a
// Gaussian model.
//
@@ -66,8 +66,8 @@
// #include "otbMRFEnergyPotts.h"
// #include "otbMRFEnergyGaussianClassification.h"
// #include "otbMRFOptimizerICM.h"
-#include "otbMRFSamplerMAP.h"
-// #include "otbMRFSamplerRandomMAP.h"
+// #include "otbMRFSamplerMAP.h"
+#include "otbMRFSamplerRandomMAP.h"
#include "otbMRFOptimizerICM.h"
// #include "otbMRFSamplerRandom.h"
// Software Guide : EndCodeSnippet
@@ -161,7 +161,7 @@ int main(int argc, char* argv[] )
// Software Guide : BeginCodeSnippet
- typedef otb::MRFSamplerMAP< InputImageType, LabelledImageType> SamplerType;
+ typedef otb::MRFSamplerRandomMAP< InputImageType, LabelledImageType> SamplerType;
// Software Guide : EndCodeSnippet
--
GitLab