diff --git a/Code/Common/otbHistogramStatisticFunction.h b/Code/Common/otbHistogramStatisticFunction.h index cfc177a6b37a7bf2851ad53314260c3f9b5083ea..14e4971ab01d03bbc07bcfe5f5a95ff45e5979f6 100644 --- a/Code/Common/otbHistogramStatisticFunction.h +++ b/Code/Common/otbHistogramStatisticFunction.h @@ -60,20 +60,38 @@ public: /** Returns the entropy value */ OutputType GetEntropy(); - + + /** Stores the histogram pointer */ + void SetInputHistogram( const TInputHistogram * histogram ) + { + if ( m_InputHistogram != histogram ) + { + m_InputHistogram = histogram ; + this->Modified() ; + m_IsModified = true; + } + } + protected: + HistogramStatisticFunction() ; virtual ~HistogramStatisticFunction() {} void PrintSelf(std::ostream& os, itk::Indent indent) const; /** Calculates the thresholds and save them */ - void GenerateData() ; + void GenerateData(); + + /** Calculate the entropy value */ + void CalculateEntropy(); private: OutputType m_entropy ; - + bool m_IsModified; + /** Target histogram data pointer */ + typename TInputHistogram::ConstPointer m_InputHistogram ; + } ; // end of class } // end of namespace otb diff --git a/Code/Common/otbHistogramStatisticFunction.txx b/Code/Common/otbHistogramStatisticFunction.txx index 20c44cdecacfc3da38fac53d24183f3b7e8c55c4..9b7e9d0a5562f9d2db65fbb6d127333779b9fa59 100644 --- a/Code/Common/otbHistogramStatisticFunction.txx +++ b/Code/Common/otbHistogramStatisticFunction.txx @@ -28,14 +28,27 @@ template< class TInputHistogram, class TOutput > HistogramStatisticFunction< TInputHistogram, TOutput> ::HistogramStatisticFunction() { + m_IsModified = true; } - + template< class TInputHistogram, class TOutput > typename HistogramStatisticFunction< TInputHistogram, TOutput>::OutputType HistogramStatisticFunction< TInputHistogram, TOutput> ::GetEntropy() { - typename TInputHistogram::ConstPointer histogram = this->GetInputHistogram(); + if(m_IsModified == true) + { + this->Update(); + } + return m_entropy; +} + +template< class TInputHistogram, class TOutput > +void +HistogramStatisticFunction< TInputHistogram, TOutput> +::CalculateEntropy() +{ + typename TInputHistogram::ConstPointer histogram = m_InputHistogram; // TODO: as an improvement, the class could accept multi-dimensional histograms // and the user could specify the dimension to apply the algorithm to. @@ -50,29 +63,35 @@ HistogramStatisticFunction< TInputHistogram, TOutput> EntropyType entropy = itk::NumericTraits<EntropyType>::Zero; FrequencyType globalFrequency = histogram->GetTotalFrequency(); - + if(globalFrequency == 0) + { + itkExceptionMacro(<<"Histogram must contain at least 1 element."); + } while (iter != end) { - EntropyType Proba = static_cast<EntropyType>(iter.GetMeasurementVector()[0]); - Proba /= static_cast<EntropyType>(globalFrequency); - - entropy += Proba * log (Proba); + EntropyType Proba = static_cast<EntropyType>(iter.GetFrequency()); + Proba /= static_cast<EntropyType>(globalFrequency); + if(Proba !=0.0) + { + entropy -= Proba * log(Proba); + } ++iter ; } m_entropy = static_cast<OutputType>(entropy); - return m_entropy; } - template< class TInputHistogram, class TOutput > void HistogramStatisticFunction< TInputHistogram, TOutput> ::GenerateData() { - m_entropy=GetEntropy(); + + CalculateEntropy(); + m_IsModified = false; } + template< class TInputHistogram, class TOutput > void HistogramStatisticFunction< TInputHistogram, TOutput> diff --git a/Testing/Code/Common/CMakeLists.txt b/Testing/Code/Common/CMakeLists.txt index 6fc17db6839f44067afcc4f749a35799ce7b28fb..20d154ba1445a34e0bc606e0602ec9bcf913ac21 100644 --- a/Testing/Code/Common/CMakeLists.txt +++ b/Testing/Code/Common/CMakeLists.txt @@ -224,7 +224,7 @@ ADD_TEST(coTuPathListToHistogramGenerator ${COMMON_TESTS} # ------- otb::HistogramStatisticFunction --------------------------- ADD_TEST(coTuHistogramStatisticFunction ${COMMON_TESTS} - otbHistogramStatisticFunction) + otbHistogramStatisticFunction 100) # ------- Fichiers sources CXX ----------------------------------- SET(BasicCommon_SRCS diff --git a/Testing/Code/Common/otbHistogramStatisticFunction.cxx b/Testing/Code/Common/otbHistogramStatisticFunction.cxx index fed608975a378abdbf6a84f9adaa68d5528925af..b0736a685c0702f03a4d6219fc17f504a1d921a3 100644 --- a/Testing/Code/Common/otbHistogramStatisticFunction.cxx +++ b/Testing/Code/Common/otbHistogramStatisticFunction.cxx @@ -26,29 +26,31 @@ #include "itkHistogram.h" #include "otbHistogramStatisticFunction.h" -int otbHistogramStatisticFunction(int, char*[]) +int otbHistogramStatisticFunction(int argc, char* argv[]) { + unsigned int NbOfBins((unsigned int)::atoi(argv[1])); + typedef float MeasurementType ; typedef itk::Statistics::Histogram< MeasurementType, 1 > HistogramType ; HistogramType::Pointer histogram = HistogramType::New() ; // initialize histogram HistogramType::SizeType size; - size.Fill(64) ; + size.Fill(NbOfBins) ; HistogramType::MeasurementVectorType lowerBound ; HistogramType::MeasurementVectorType upperBound ; lowerBound[0] = 0.0 ; - upperBound[0] = 64.0 ; + upperBound[0] = NbOfBins ; histogram->Initialize(size, lowerBound, upperBound ) ; - // create histogram + // create histogram with same value for each frequency for (HistogramType::Iterator iter = histogram->Begin(); iter != histogram->End(); ++iter) { iter.SetFrequency(1.0); } - + typedef otb::HistogramStatisticFunction<HistogramType,MeasurementType> HistogramStatisticFunctionType; HistogramStatisticFunctionType::Pointer HistogramStatisticFunction = HistogramStatisticFunctionType::New(); @@ -57,11 +59,34 @@ int otbHistogramStatisticFunction(int, char*[]) MeasurementType Entropy; Entropy = HistogramStatisticFunction->GetEntropy(); - std::cout << "Entropy : " << Entropy << std::endl; + std::cout << "Entropy 1 : " << Entropy << std::endl; + + if(fabs(Entropy-log(NbOfBins))>0.00001 ) + { + std::cout << "Error in entropy estimation" << std::endl; + return EXIT_FAILURE; + } - if(fabs(Entropy)>0.0 ) + // create histogram just all value equal to zero except the first one + for (HistogramType::Iterator iter = histogram->Begin(); iter != histogram->End(); ++iter) + { + if(iter == histogram->Begin()) + { + iter.SetFrequency(1.0); + } + else + { + iter.SetFrequency(0.0); + } + } + + HistogramStatisticFunction->Update(); + Entropy = HistogramStatisticFunction->GetEntropy(); + std::cout << "Entropy 2 : " << Entropy << std::endl; + if( Entropy!=0.0 ) { - return EXIT_FAILURE; + std::cout << "Error in entropy estimation" << std::endl; + return EXIT_FAILURE; } return EXIT_SUCCESS;