From 0b03f59ca75a104b37b4c0d006a865a3dfcdecb1 Mon Sep 17 00:00:00 2001
From: Cyrille Valladeau <cyrille.valladeau@c-s.fr>
Date: Fri, 30 Jan 2009 19:46:50 +0100
Subject: [PATCH] ENH : ENTROPY FUNCTOR added
---
...otbFunctionWithNeighborhoodToImageFilter.h | 1 +
...bFunctionWithNeighborhoodToImageFilter.txx | 12 +-
...FunctorNeighborhoodWithOffsetImageFilter.h | 4 +
...nctorNeighborhoodWithOffsetImageFilter.txx | 20 +-
.../otbEnergyTextureFunctor.h | 3 +-
.../otbEntropyTextureFunctor.h | 244 ++++++++++++++++++
Testing/Code/FeatureExtraction/CMakeLists.txt | 13 +
.../otbEntropyTextureFunctor.cxx | 65 +++++
.../otbFeatureExtractionTests11.cxx | 1 +
9 files changed, 358 insertions(+), 5 deletions(-)
create mode 100644 Code/FeatureExtraction/otbEntropyTextureFunctor.h
create mode 100644 Testing/Code/FeatureExtraction/otbEntropyTextureFunctor.cxx
diff --git a/Code/BasicFilters/otbFunctionWithNeighborhoodToImageFilter.h b/Code/BasicFilters/otbFunctionWithNeighborhoodToImageFilter.h
index 968f0e920a..0b37640da2 100644
--- a/Code/BasicFilters/otbFunctionWithNeighborhoodToImageFilter.h
+++ b/Code/BasicFilters/otbFunctionWithNeighborhoodToImageFilter.h
@@ -134,6 +134,7 @@ private:
InputImageSizeType m_Radius;
InputImageOffsetType m_Offset;
+ //std::vector<FunctionType*> m_FunctionVector;
};
} // end namespace otb
diff --git a/Code/BasicFilters/otbFunctionWithNeighborhoodToImageFilter.txx b/Code/BasicFilters/otbFunctionWithNeighborhoodToImageFilter.txx
index 6ea469338e..f51a774807 100644
--- a/Code/BasicFilters/otbFunctionWithNeighborhoodToImageFilter.txx
+++ b/Code/BasicFilters/otbFunctionWithNeighborhoodToImageFilter.txx
@@ -36,6 +36,7 @@ FunctionWithNeighborhoodToImageFilter<TInputImage,TOutputImage,TFunction>
{
m_Radius.Fill(0);
m_Offset.Fill(0);
+ //m_FunctionVector.clear();
}
@@ -47,7 +48,16 @@ FunctionWithNeighborhoodToImageFilter<TInputImage,TOutputImage,TFunction>
Superclass::BeforeThreadedGenerateData();
m_Radius = this->GetFunction()->GetRadius();
- m_Offset = this->GetFunction()->GetOffset();
+ m_Offset = this->GetFunction()->GetOffset();
+ /*
+ for(unsigned int i =0; i<this->GetNumberOfThreads(); i++)
+ {
+ FunctionType * func;
+ func = (this->GetFunction());
+ m_FunctionVector.push_back(func);
+ std::cout<<func<<std::endl;
+ }
+ */
}
template <class TInputImage, class TOutputImage, class TFunction >
diff --git a/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.h b/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.h
index 159f6ce1f9..150046c004 100644
--- a/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.h
+++ b/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.h
@@ -120,6 +120,8 @@ protected:
*/
virtual ~UnaryFunctorNeighborhoodWithOffsetImageFilter() {};
+ virtual void BeforeThreadedGenerateData();
+
/** UnaryFunctorNeighborhoodWithOffsetImageFilter can be implemented as a multithreaded filter.
* Therefore, this implementation provides a ThreadedGenerateData() routine
* which is called for each processing thread. The output image data is
@@ -136,6 +138,7 @@ protected:
* Pad the input requested region by radius
*/
virtual void GenerateInputRequestedRegion(void);
+ std::vector<FunctorType> m_FunctorList;
private:
UnaryFunctorNeighborhoodWithOffsetImageFilter(const Self&); //purposely not implemented
@@ -144,6 +147,7 @@ private:
unsigned int m_Radius;
FunctorType m_Functor;
+
InputImageOffsetType m_Offset;
};
diff --git a/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.txx b/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.txx
index e537531bc2..34b84c801a 100644
--- a/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.txx
+++ b/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.txx
@@ -39,9 +39,24 @@ UnaryFunctorNeighborhoodWithOffsetImageFilter<TInputImage,TOutputImage,TFunction
this->SetNumberOfRequiredInputs( 1 );
m_Radius = 1;
m_Offset.Fill(1);
+ m_FunctorList.clear();
}
+template <class TInputImage, class TOutputImage, class TFunction >
+void
+UnaryFunctorNeighborhoodWithOffsetImageFilter<TInputImage,TOutputImage,TFunction>
+::BeforeThreadedGenerateData()
+{
+ Superclass::BeforeThreadedGenerateData();
+
+ for(unsigned int i =0; i<this->GetNumberOfThreads(); i++)
+ {
+ m_FunctorList.push_back(m_Functor);
+ }
+ }
+
+
template <class TInputImage, class TOutputImage, class TFunction >
void
UnaryFunctorNeighborhoodWithOffsetImageFilter<TInputImage,TOutputImage,TFunction>
@@ -103,7 +118,7 @@ template <class TInputImage, class TOutputImage, class TFunction >
void
UnaryFunctorNeighborhoodWithOffsetImageFilter<TInputImage, TOutputImage, TFunction>
::ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, int threadId)
-{
+{ std::cout<<"threadId : "<<threadId<<std::endl;
itk::ZeroFluxNeumannBoundaryCondition<TInputImage> nbc;
itk::ZeroFluxNeumannBoundaryCondition<TInputImage> nbcOff;
// We use dynamic_cast since inputs are stored as DataObjects. The
@@ -159,7 +174,7 @@ UnaryFunctorNeighborhoodWithOffsetImageFilter<TInputImage, TOutputImage, TFuncti
while ( ! outputIt.IsAtEnd() )
{
- outputIt.Set( m_Functor( neighInputIt, neighInputOffIt) );
+ outputIt.Set( m_FunctorList[threadId]( neighInputIt, neighInputOffIt) );
++neighInputIt;
++neighInputOffIt;
@@ -169,6 +184,7 @@ UnaryFunctorNeighborhoodWithOffsetImageFilter<TInputImage, TOutputImage, TFuncti
++fit;
++fitOff;
}
+std::cout<<"threadIdFIN : "<<threadId<<std::endl;
}
} // end namespace otb
diff --git a/Code/FeatureExtraction/otbEnergyTextureFunctor.h b/Code/FeatureExtraction/otbEnergyTextureFunctor.h
index e04ed396c3..ac19adf18f 100644
--- a/Code/FeatureExtraction/otbEnergyTextureFunctor.h
+++ b/Code/FeatureExtraction/otbEnergyTextureFunctor.h
@@ -82,7 +82,6 @@ class EnergyTextureFunctor
{
offsetOff[0]++;
offsetOff[1] = offsetOffInit[1];
-
offset[0] = l;
ll = l*l;
for( int k = -static_cast<int>(radius[1]); k <= static_cast<int>(radius[1]); k++)
@@ -90,7 +89,7 @@ class EnergyTextureFunctor
offsetOff[1]++;
offset[1] = k;
norm = vcl_pow(static_cast<double>(itOff.GetPixel(offsetOff)[i]-itOff.GetCenterPixel()[i]), 2);
- temp += norm;
+ temp += norm;
}
temp /= area;
outPix[i] = static_cast<OutputPixelType>( vcl_pow(temp, 2) );
diff --git a/Code/FeatureExtraction/otbEntropyTextureFunctor.h b/Code/FeatureExtraction/otbEntropyTextureFunctor.h
new file mode 100644
index 0000000000..65820cf48b
--- /dev/null
+++ b/Code/FeatureExtraction/otbEntropyTextureFunctor.h
@@ -0,0 +1,244 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbEntropyTextureFunctor_h
+#define __ootbEntropyTextureFunctor_h
+
+#include "otbMath.h"
+
+namespace otb
+{
+namespace Functor
+ {
+ /** \class EntropyTextureFunctor
+ * \brief This functor calculates the local energy of an image
+ *
+ * Computes the sqaure gradient which is computed using offset and
+ * angle pixel of the neighborhood.
+ * TIterInput is an ietrator, TOutput is a PixelType.
+ *
+ * \ingroup Functor
+ * \ingroup Statistics
+ */
+template <class TIterInput1, class TIterInput2, class TOutput>
+class EntropyTextureFunctor
+ {
+ public:
+ EntropyTextureFunctor()
+ {
+ m_Offset.Fill(1);
+ //m_BinLength = IntVectyorType(0, 1);
+ };
+ ~EntropyTextureFunctor(){};
+
+ typedef TIterInput1 IterType1;
+ typedef TIterInput2 IterType2;
+ typedef TOutput OutputType;
+ typedef typename IterType1::OffsetType OffsetType;
+ typedef typename IterType1::RadiusType RadiusType;
+ typedef typename OutputType::ValueType OutputPixelType;
+ typedef std::vector<double> DoubleVectorType;
+ typedef std::vector<int> IntVectorType;
+ typedef std::vector<IntVectorType> IntVectorVectorType;
+
+ void SetOffset(OffsetType off){ m_Offset=off; };
+ OffsetType GetOffset(){ return m_Offset; };
+
+ IntVectorVectorType StatComputation(const TIterInput1 &it, const TIterInput2 &itOff)
+ {
+ IntVectorVectorType output;
+ IntVectorType binLength;
+ IntVectorType binLengthOff;
+ RadiusType radius = it.GetRadius();
+ unsigned int nbComp = it.GetCenterPixel().GetSize();
+ m_Mini = DoubleVectorType(nbComp, itk::NumericTraits<double>::max());
+ m_MiniOff = DoubleVectorType(nbComp, itk::NumericTraits<double>::max());
+ m_Maxi = DoubleVectorType(nbComp, itk::NumericTraits<double>::NonpositiveMin());
+ m_MaxiOff = DoubleVectorType(nbComp, itk::NumericTraits<double>::NonpositiveMin());
+ double area = static_cast<double>(radius[0]*radius[1]);
+ double areaInv = 1/area;
+
+ OffsetType offset;
+ offset.Fill(0);
+ OffsetType offsetOff;
+ OffsetType offsetOffInit;
+
+ offsetOffInit[0] = -radius[0]+m_Offset[0]-1;
+ offsetOffInit[1] = -radius[1]+m_Offset[1]-1;
+ IntVectorType mean;
+ IntVectorType meanOff;
+ for( unsigned int i=0; i<nbComp; i++ )
+ {
+ offsetOff = offsetOffInit;
+ double meanCh = 0.;
+ double meanChOff = 0.;
+ for( int l = -static_cast<int>(radius[0]); l <= static_cast<int>(radius[0]); l++ )
+ {
+ offsetOff[0]++;
+ offsetOff[1] = offsetOffInit[1];
+ offset[0] = l;
+ for( int k = -static_cast<int>(radius[1]); k <= static_cast<int>(radius[1]); k++)
+ {
+ //std::cout<<it.GetPixel(offset)[i]<<std::endl;
+ offsetOff[1]++;
+ offset[1] = k;
+ meanCh += static_cast<double>(it.GetPixel(offset)[i]);
+ meanChOff += static_cast<double>(itOff.GetPixel(offsetOff)[i]);
+
+ m_Mini[i] = std::min(static_cast<double>(it.GetPixel(offset)[i]), m_Mini[i]);
+ m_MiniOff[i] = std::min(static_cast<double>(itOff.GetPixel(offsetOff)[i]), m_MiniOff[i]);
+ m_Maxi[i] = std::max(static_cast<double>(it.GetPixel(offset)[i]), m_Maxi[i]);
+ m_MaxiOff[i] = std::max(static_cast<double>(itOff.GetPixel(offsetOff)[i]), m_MaxiOff[i]);
+ }
+ mean.push_back(meanCh * areaInv);
+ meanOff.push_back(meanChOff * areaInv);
+ }
+ }
+
+ for( unsigned int i=0; i<nbComp; i++ )
+ {
+ offsetOff = offsetOffInit;
+ double stdCh = 0.;
+ double stdChOff = 0.;
+ for( int l = -static_cast<int>(radius[0]); l <= static_cast<int>(radius[0]); l++ )
+ {
+ offsetOff[0]++;
+ offsetOff[1] = offsetOffInit[1];
+ offset[0] = l;
+ for( int k = -static_cast<int>(radius[1]); k <= static_cast<int>(radius[1]); k++)
+ {
+ offsetOff[1]++;
+ offset[1] = k;
+ stdCh += vcl_pow( (mean[i]-static_cast<double>(it.GetPixel(offset)[i])), 2);
+ stdChOff += vcl_pow( (meanOff[i]-static_cast<double>(itOff.GetPixel(offsetOff)[i])), 2);
+ }
+ }
+ stdCh *= areaInv;
+ stdCh = vcl_sqrt( stdCh );
+ stdChOff *= areaInv;
+ stdChOff = vcl_sqrt( stdChOff );
+ binLength.push_back( (3.5*stdCh) /(vcl_pow(area, 1/3) ) );
+ binLengthOff.push_back( (3.5*stdChOff) /(vcl_pow(area, 1/3) ) );
+ }
+
+ output.push_back(binLength);
+ output.push_back(binLengthOff);
+
+ return output;
+ }
+
+
+ inline TOutput operator()(const TIterInput1 &it, const TIterInput2 &itOff)
+ {
+ //std::cout<<"operator"<<std::endl;
+ IntVectorVectorType binsLength = this->StatComputation(it, itOff);
+
+ RadiusType radius = it.GetRadius();
+ unsigned int nbComp = it.GetCenterPixel().GetSize();
+ double area = static_cast<double>(radius[0]*radius[1]);
+ double areaInv = 1/area;
+ OutputType outPix;
+ outPix.SetSize( nbComp );
+ outPix.Fill(0);
+
+
+ OffsetType offset;
+ offset.Fill(0);
+ OffsetType offsetOff;
+ OffsetType offsetOffInit;
+
+ offsetOffInit[0] = -radius[0]+m_Offset[0]-1;
+ offsetOffInit[1] = -radius[1]+m_Offset[1]-1;
+
+ int histoIdX = 0;
+ int histoIdY = 0;
+
+ for( unsigned int i=0; i<nbComp; i++ )
+ {
+ //std::cout<<"i :"<<i<<" "<<std::endl;
+
+ IntVectorType histoTemp;
+ IntVectorVectorType histo;
+ if(binsLength[0][i] != 0)
+ histoTemp = IntVectorType( vcl_floor((m_Maxi[i]-m_Mini[i])/binsLength[0][i])+1, 0 );
+ else
+ histoTemp = IntVectorType( 1, 0 );
+ if(binsLength[1][i] != 0)
+ histo = IntVectorVectorType( vcl_floor((m_MaxiOff[i]-m_MiniOff[i])/binsLength[1][i])+1, histoTemp );
+ else
+ histo = IntVectorVectorType( 1, histoTemp );
+
+ //std::cout<<"m_Maxi[i]-m_Mini[i] "<<m_Maxi[i]<<" "<<m_Mini[i]<<std::endl;
+ //std::cout<<"binsLength[0][i] "<<binsLength[0][i]<<std::endl;
+ //std::cout<<histoTemp.size()<<" "<<histo.size()<<std::endl;
+
+ offsetOff = offsetOffInit;
+ for( int l = -static_cast<int>(radius[0]); l <= static_cast<int>(radius[0]); l++ )
+ { //std::cout<<"i :"<<i<<" "<<l<<std::endl;
+ offsetOff[0]++;
+ offsetOff[1] = offsetOffInit[1];
+ offset[0] = l;
+ for( int k = -static_cast<int>(radius[1]); k <= static_cast<int>(radius[1]); k++)
+ {
+ //std::cout<<i<<" "<<l<<" "<<k<<std::endl;
+ offsetOff[1]++;
+ offset[1] = k;
+ if( binsLength[1][i] != 0)
+ histoIdX = vcl_floor( (itOff.GetPixel(offsetOff)[i]-m_MiniOff[i]) / binsLength[1][i] );
+ else
+ histoIdX = 0;
+ if( binsLength[0][i] !=0 )
+ histoIdY = vcl_floor( (it.GetPixel(offset)[i]-m_Mini[i]) / binsLength[0][i] );
+ else
+ histoIdY = 0;
+ histo[histoIdX][histoIdY]++;
+
+ }
+ }
+ for(unsigned r = 0; r<binsLength.size(); r++)
+ {
+ for(unsigned s = 0; s<binsLength[r].size(); s++)
+ {
+ double p = binsLength[r][s] * areaInv;
+ if(p != 0)
+ outPix[i] += static_cast<OutputPixelType>(p * vcl_log(p));
+ }
+ }
+
+ }
+
+ //std::cout<<"operator FIN"<<std::endl;
+ return outPix;
+
+ }
+
+ private:
+ OffsetType m_Offset;
+ DoubleVectorType m_Mini;
+ DoubleVectorType m_MiniOff;
+ DoubleVectorType m_Maxi;
+ DoubleVectorType m_MaxiOff;
+ };
+
+
+
+
+ } // namespace Functor
+} // namespace otb
+
+#endif
+
diff --git a/Testing/Code/FeatureExtraction/CMakeLists.txt b/Testing/Code/FeatureExtraction/CMakeLists.txt
index 7361936bf1..a4dd84ccb5 100644
--- a/Testing/Code/FeatureExtraction/CMakeLists.txt
+++ b/Testing/Code/FeatureExtraction/CMakeLists.txt
@@ -1052,6 +1052,18 @@ ADD_TEST(feTvEnergyTextureFunctor ${FEATUREEXTRACTION_TESTS11}
2 # offset[1]
)
+# ------- otb::EntropyTextureFunctor -------------
+ADD_TEST(feTvEntropyTextureFunctor ${FEATUREEXTRACTION_TESTS11}
+--compare-image ${EPS}
+ ${BASELINE}/feTvEntropyTextureFunctor.tif
+ ${TEMP}/feTvEntropyTextureFunctor.tif
+ otbEntropyTextureFunctor
+ ${INPUTDATA}/poupees.tif
+ ${TEMP}/feTvEntropyTextureFunctor.tif
+ 5 # radius
+ -2 # offset[0]
+ 2 # offset[1]
+)
# ------- otb::EnergyTextureImageFunction -------------
ADD_TEST(feTvEnergyTextureImageFunctionNew ${FEATUREEXTRACTION_TESTS11}
@@ -1203,6 +1215,7 @@ otbSimplifyManyPathListFilter.cxx
otbEnergyTextureFunctor.cxx
otbEnergyTextureImageFunctionNew.cxx
otbEnergyTextureImageFunction.cxx
+otbEntropyTextureFunctor.cxx
)
INCLUDE_DIRECTORIES("${OTBTesting_BINARY_DIR}")
diff --git a/Testing/Code/FeatureExtraction/otbEntropyTextureFunctor.cxx b/Testing/Code/FeatureExtraction/otbEntropyTextureFunctor.cxx
new file mode 100644
index 0000000000..fbaaaf0e8f
--- /dev/null
+++ b/Testing/Code/FeatureExtraction/otbEntropyTextureFunctor.cxx
@@ -0,0 +1,65 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+#include "itkExceptionObject.h"
+
+#include "otbUnaryFunctorNeighborhoodWithOffsetImageFilter.h"
+#include "otbVectorImage.h"
+#include "itkConstNeighborhoodIterator.h"
+#include "otbImageFileReader.h"
+#include "otbStreamingImageFileWriter.h"
+#include "otbEntropyTextureFunctor.h"
+
+
+int otbEntropyTextureFunctor(int argc, char * argv[])
+{
+ const char * inputFileName = argv[1];
+ const char * outputFileName = argv[2];
+
+ typedef double InputPixelType;
+ const int Dimension = 2;
+ typedef otb::VectorImage<InputPixelType,Dimension> ImageType;
+ typedef ImageType::PixelType PixelType;
+ typedef ImageType::OffsetType OffsetType;
+ typedef otb::ImageFileReader<ImageType> ReaderType;
+ typedef otb::StreamingImageFileWriter<ImageType> WriterType;
+
+ typedef itk::ConstNeighborhoodIterator<ImageType> IterType;;
+ typedef otb::Functor::EntropyTextureFunctor<IterType, IterType, PixelType> FunctorType;
+ typedef otb::UnaryFunctorNeighborhoodWithOffsetImageFilter<ImageType, ImageType, FunctorType> UnaryFunctorNeighborhoodImageFilterType;
+
+ // Instantiating object
+ UnaryFunctorNeighborhoodImageFilterType::Pointer object = UnaryFunctorNeighborhoodImageFilterType::New();
+ ReaderType::Pointer reader = ReaderType::New();
+ WriterType::Pointer writer = WriterType::New();
+ reader->SetFileName(inputFileName);
+ writer->SetFileName(outputFileName);
+
+ object->SetInput(reader->GetOutput());
+ object->SetRadius(atoi(argv[3]));
+ OffsetType offset;
+ offset[0] = atoi(argv[4]);
+ offset[1] = atoi(argv[5]);
+
+ object->SetOffset(offset);
+ writer->SetInput(object->GetOutput());
+
+ writer->Update();
+
+
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/FeatureExtraction/otbFeatureExtractionTests11.cxx b/Testing/Code/FeatureExtraction/otbFeatureExtractionTests11.cxx
index cba71ce767..c8e891cc6f 100644
--- a/Testing/Code/FeatureExtraction/otbFeatureExtractionTests11.cxx
+++ b/Testing/Code/FeatureExtraction/otbFeatureExtractionTests11.cxx
@@ -36,4 +36,5 @@ REGISTER_TEST(otbSimplifyManyPathListFilter);
REGISTER_TEST(otbEnergyTextureFunctor);
REGISTER_TEST(otbEnergyTextureImageFunctionNew);
REGISTER_TEST(otbEnergyTextureImageFunction);
+REGISTER_TEST(otbEntropyTextureFunctor);
}
--
GitLab