From 2daf417f13be0c099b128a04ff5e8605002e8810 Mon Sep 17 00:00:00 2001
From: Cyrille Valladeau <cyrille.valladeau@c-s.fr>
Date: Wed, 7 May 2008 13:35:17 +0000
Subject: [PATCH] Modif des Markow : pb d'index et d'nitilaisation de
 m_NumeberOfClasses.

---
 Code/Markov/otbMRFSampler.h          |  3 +-
 Code/Markov/otbMRFSamplerMAP.h       | 23 +++++++++++++-
 Code/Markov/otbMRFSamplerRandomMAP.h | 46 +++++++++++++++++++++++-----
 3 files changed, 62 insertions(+), 10 deletions(-)

diff --git a/Code/Markov/otbMRFSampler.h b/Code/Markov/otbMRFSampler.h
index 720c481f85..f195e35faa 100644
--- a/Code/Markov/otbMRFSampler.h
+++ b/Code/Markov/otbMRFSampler.h
@@ -57,7 +57,7 @@ class ITK_EXPORT MRFSampler:public itk::Object
     itkTypeMacro(MRFSampler,itk::Object);
  
     // Accessors
-    itkSetMacro(NumberOfClasses, unsigned int);
+    virtual void SetNumberOfClasses(unsigned int nb){ m_NumberOfClasses = nb; this->Modified(); };
     itkGetMacro(NumberOfClasses, unsigned int);
 
     itkSetMacro(Lambda, double);
@@ -107,6 +107,7 @@ class ITK_EXPORT MRFSampler:public itk::Object
 	m_EnergyRegularization = EnergyRegularizationType::New();
 	m_EnergyFidelity = EnergyFidelityType::New();
 
+    m_NumberOfClasses = 1;
 	m_Lambda = 1.;
 	m_EnergyCurrent = 1.;
 	m_DeltaEnergy = 1.;
diff --git a/Code/Markov/otbMRFSamplerMAP.h b/Code/Markov/otbMRFSamplerMAP.h
index 0a07dcc476..f050166f9f 100644
--- a/Code/Markov/otbMRFSamplerMAP.h
+++ b/Code/Markov/otbMRFSamplerMAP.h
@@ -74,7 +74,8 @@ class ITK_EXPORT MRFSamplerMAP : public MRFSampler< TInput1, TInput2 >
       //Try all possible value (how to be generic ?)
       this->SetEnergyAfter( this->GetEnergyBefore() ); //default values to current one
       this->SetValue( itRegul.GetCenterPixel() );
-      for (LabelledImagePixelType valueCurrent = 0; valueCurrent< this->m_NumberOfClasses; ++valueCurrent)
+    /*
+      fofor (LabelledImagePixelType valueCurrent = 0; valueCurrent< this->m_NumberOfClasses; ++valueCurrent)
 	{
 	  this->SetEnergyCurrent( this->GetEnergyFidelity()->GetValue(itData, valueCurrent)
 				  + this->GetLambda() * this->GetEnergyRegularization()->GetValue(itRegul, valueCurrent) );  
@@ -85,6 +86,26 @@ class ITK_EXPORT MRFSamplerMAP : public MRFSampler< TInput1, TInput2 >
 	    }
 	  if (valueCurrent == itk::NumericTraits<LabelledImagePixelType>::max()) break;
 	}
+          */
+    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() );
       
diff --git a/Code/Markov/otbMRFSamplerRandomMAP.h b/Code/Markov/otbMRFSamplerRandomMAP.h
index 375d7f1ffa..08c80d5f4b 100644
--- a/Code/Markov/otbMRFSamplerRandomMAP.h
+++ b/Code/Markov/otbMRFSamplerRandomMAP.h
@@ -97,12 +97,20 @@ class ITK_EXPORT MRFSamplerRandomMAP : public MRFSampler< TInput1, TInput2>
     
     inline int Compute( const InputImageNeighborhoodIterator & itData, const LabelledImageNeighborhoodIterator & itRegul)             
       {
+        if(this->GetNumberOfClasses() == 0)
+        {
+            itkExceptionMacro(<<"m_NumberOfClasses has to be greater than 0.");
+        }
+        std::cout<<this->GetEnergyAfter()<<std::endl;
+        std::cout<<this->GetEnergyBefore()<<std::endl;
+        
 	this->SetEnergyBefore( this->GetEnergyFidelity()->GetValue(itData, itRegul.GetCenterPixel())
 			       + this->GetLambda() * this->GetEnergyRegularization()->GetValue(itRegul, itRegul.GetCenterPixel()) );
-	
+		std::cout<<this->GetEnergyBefore()<<std::endl;
 	//Try all possible value (how to be generic ?)
 	this->SetEnergyAfter( this->GetEnergyBefore() ); //default values to current one
 	this->SetValue( itRegul.GetCenterPixel() );
+	std::cout<<this->GetEnergyAfter()<<std::endl;
 	// otbDebugMacro(<< "Computing MAP for pix " << itData.GetIndex());
 	// Compute probability for each possibility
 	double totalProba=0.0;
@@ -111,8 +119,12 @@ class ITK_EXPORT MRFSamplerRandomMAP : public MRFSampler< TInput1, TInput2>
 	    // otbDebugMacro(<< " --> Proposed value " << static_cast<double>(valueCurrent)); 
 	    this->SetEnergyCurrent( this->GetEnergyFidelity()->GetValue(itData, valueCurrent)
 	      + this->GetLambda() * this->GetEnergyRegularization()->GetValue(itRegul, valueCurrent) );
-	    
+          
+        std::cout<<this->GetEnergyFidelity()->GetValue(itData, valueCurrent)<<" ��� "<<this->GetLambda()<<" ��� "<<this->GetEnergyRegularization()->GetValue(itRegul, valueCurrent)<<std::endl;
+	    std::cout<<this->GetEnergyCurrent()<<std::endl;
+        std::cout<<valueCurrent<<" ### "<<m_Energy[valueCurrent] <<std::endl;
 	    m_Energy[valueCurrent] = this->GetEnergyCurrent();
+        std::cout<<valueCurrent<<" ### "<<m_Energy[valueCurrent] <<std::endl;
 	    m_RepartitionFunction[valueCurrent] = vcl_exp(-this->GetEnergyCurrent())+totalProba;
 	    totalProba = m_RepartitionFunction[valueCurrent];
 	    // otbDebugMacro("valueCurrent, m_RepartitionFunction[valueCurrent] " << (unsigned int)  valueCurrent << ", " << m_RepartitionFunction[valueCurrent]);
@@ -136,15 +148,22 @@ class ITK_EXPORT MRFSamplerRandomMAP : public MRFSampler< TInput1, TInput2>
 	
 	// otbDebugMacro("select, totalProba " <<  select << ", " << totalProba);
 	unsigned int valueCurrent = 0;
-	for (valueCurrent = 0; valueCurrent < this->m_NumberOfClasses; ++valueCurrent)
-	  {
-	    if (m_RepartitionFunction[valueCurrent] > select) break;
-	  }
-	
-	if ( this->GetValue() != static_cast<LabelledImagePixelType>(valueCurrent))
+    while( valueCurrent<this->GetNumberOfClasses() && m_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] );
+        	std::cout<<valueCurrent<<"##### "<<m_Energy[valueCurrent] <<std::endl;
 	  }
 	
 	this->SetDeltaEnergy( this->GetEnergyAfter() - this->GetEnergyBefore() );
@@ -160,6 +179,17 @@ class ITK_EXPORT MRFSamplerRandomMAP : public MRFSampler< TInput1, TInput2>
 	this->Modified();
       };
 
+    void SetNumberOfClasses(unsigned int nb)
+        {
+            Superclass::SetNumberOfClasses( nb );
+            free(m_Energy);
+            free(m_RepartitionFunction);
+            m_Energy = (double *) calloc(nb, sizeof(double));
+            m_RepartitionFunction = (double *) calloc(nb, sizeof(double));
+                std::cout<<m_Energy[0]<<std::endl;
+            this->Modified(); 
+        };
+      
   private:
     double * m_RepartitionFunction ;   
     double * m_Energy; 
-- 
GitLab