From 7fc2e5dac96eb5b2093e357941f4dcbab234d257 Mon Sep 17 00:00:00 2001
From: Cyrille Valladeau <cyrille.valladeau@c-s.fr>
Date: Mon, 9 Feb 2009 15:12:12 +0100
Subject: [PATCH] ENH : texture formalas correction + add sum entropy
---
.../otbContrastTextureFunctor.h | 2 +-
.../otbDifferenceEntropyTextureFunctor.h | 5 +-
.../otbSumAverageTextureFunctor.h | 5 +-
.../otbSumEntropyTextureFunctor.h | 149 ++++++++++++++++++
Testing/Code/FeatureExtraction/CMakeLists.txt | 32 ++++
.../FeatureExtraction/otbTextureFunctor.cxx | 8 +
.../otbTextureImageFunction.cxx | 7 +-
7 files changed, 201 insertions(+), 7 deletions(-)
create mode 100644 Code/FeatureExtraction/otbSumEntropyTextureFunctor.h
diff --git a/Code/FeatureExtraction/otbContrastTextureFunctor.h b/Code/FeatureExtraction/otbContrastTextureFunctor.h
index 8e508bc525..ad1d97300d 100644
--- a/Code/FeatureExtraction/otbContrastTextureFunctor.h
+++ b/Code/FeatureExtraction/otbContrastTextureFunctor.h
@@ -120,7 +120,7 @@ public:
double rVal = (static_cast<double>(r)+0.5)*binsLength[1];
for (unsigned s = 0; s<histo[r].size(); s++)
{
- if( vcl_abs( rVal - (static_cast<double>(s)+0.5)*binsLength[0] - nCeil) < vcl_abs(binsLength[0]-binsLength[1]) )
+ if( vcl_abs((static_cast<double>(s)+0.5)*binsLength[0] - rVal - nCeil) < vcl_abs(binsLength[0]) )
{
double p = static_cast<double>(histo[r][s])*areaInv;
out += nCeilSquare * p;
diff --git a/Code/FeatureExtraction/otbDifferenceEntropyTextureFunctor.h b/Code/FeatureExtraction/otbDifferenceEntropyTextureFunctor.h
index bf667f36e2..d984f13e17 100644
--- a/Code/FeatureExtraction/otbDifferenceEntropyTextureFunctor.h
+++ b/Code/FeatureExtraction/otbDifferenceEntropyTextureFunctor.h
@@ -115,20 +115,19 @@ public:
{
double Px_y = 0.;
double nCeil = (static_cast<double>(sB)+0.5)*binsLength[0];
- double nCeilSquare = vcl_pow( nCeil, 2);
for (unsigned r = 0; r<histo.size(); r++)
{
double rVal = (static_cast<double>(r)+0.5)*binsLength[1];
for (unsigned s = 0; s<histo[r].size(); s++)
{
- if( vcl_abs( rVal - (static_cast<double>(s)+0.5)*binsLength[0] - nCeil ) < vcl_abs( binsLength[0]+binsLength[1]) )
+ if( vcl_abs((static_cast<double>(s)+0.5)*binsLength[0] - rVal - nCeil) < vcl_abs(binsLength[0]) )
{
Px_y += static_cast<double>(histo[r][s])*areaInv;
}
}
}
if(Px_y != 0.)
- out += Px_y * nCeil * vcl_log(Px_y);
+ out += Px_y * vcl_log(Px_y);
}
if(out != 0)
diff --git a/Code/FeatureExtraction/otbSumAverageTextureFunctor.h b/Code/FeatureExtraction/otbSumAverageTextureFunctor.h
index c9028a7b27..3de106d2f7 100644
--- a/Code/FeatureExtraction/otbSumAverageTextureFunctor.h
+++ b/Code/FeatureExtraction/otbSumAverageTextureFunctor.h
@@ -122,10 +122,11 @@ public:
double sVal = (static_cast<double>(s)+0.5)*binsLength[0];
// In theory don't have the abs but will deals with neighborhood and offset without the same histo
// thus loop over 2*Ng don't have sense
- if( vcl_abs(rVal + sVal - nCeil) < vcl_abs(binsLength[0]+binsLength[1]) || vcl_abs(rVal + sVal - 2*nCeil) < vcl_abs(binsLength[0]+binsLength[1]) )
+ //if( vcl_abs(rVal + sVal - nCeil) < vcl_abs(binsLength[0]+binsLength[1]) || vcl_abs(rVal + sVal - 2*nCeil) < vcl_abs(binsLength[0]+binsLength[1]) )
+ if( vcl_abs(rVal + sVal - nCeil) < vcl_abs(binsLength[0]) || vcl_abs(rVal + sVal - 2*nCeil) < 2*vcl_abs(binsLength[0]) )
{
double p = static_cast<double>(histo[r][s])*areaInv;
- out += sVal * p;
+ out += nCeil * p;
}
}
}
diff --git a/Code/FeatureExtraction/otbSumEntropyTextureFunctor.h b/Code/FeatureExtraction/otbSumEntropyTextureFunctor.h
new file mode 100644
index 0000000000..dad317125a
--- /dev/null
+++ b/Code/FeatureExtraction/otbSumEntropyTextureFunctor.h
@@ -0,0 +1,149 @@
+/*=========================================================================
+
+ 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 __otbSumEntropyTextureFunctor_h
+#define __otbSumEntropyTextureFunctor_h
+
+#include "otbTextureFunctorBase.h"
+
+namespace otb
+{
+namespace Functor
+{
+/** \class SumEntropyTextureFunctor
+ * \brief This functor calculates the inverse difference moment of an image
+ *
+ * Computes joint histogram (neighborhood and offset neighborhood)
+ * which bins are computing using Scott formula.
+ * Computes the probabiltiy p for each pair of pixel.
+ * InverseSumMoment is the sum 1/(1+(pi-poff)²)*p over the neighborhood.
+ * TIterInput is an ietrator, TOutput is a PixelType.
+ *
+ * \ingroup Functor
+ * \ingroup
+ * \ingroup Statistics
+ */
+template <class TIterInput1, class TIterInput2, class TOutput>
+class ITK_EXPORT SumEntropyTextureFunctor :
+public TextureFunctorBase<TIterInput1, TIterInput2, TOutput>
+{
+public:
+ SumEntropyTextureFunctor(){};
+ ~SumEntropyTextureFunctor(){};
+
+ typedef TIterInput1 IterType1;
+ typedef TIterInput2 IterType2;
+ typedef TOutput OutputType;
+ typedef typename IterType1::OffsetType OffsetType;
+ typedef typename IterType1::RadiusType RadiusType;
+ typedef typename IterType1::InternalPixelType InternalPixelType;
+ typedef typename IterType1::ImageType ImageType;
+ typedef itk::Neighborhood<InternalPixelType,::itk::GetImageDimension<ImageType>::ImageDimension> NeighborhoodType;
+ typedef std::vector<double> DoubleVectorType;
+ typedef std::vector<int> IntVectorType;
+ typedef std::vector<IntVectorType> IntVectorVectorType;
+
+
+ virtual double ComputeOverSingleChannel(const NeighborhoodType &neigh, const NeighborhoodType &neighOff)
+ {
+ DoubleVectorType binsLength = this->StatComputation(neigh, neighOff);
+
+ RadiusType radius = neigh.GetRadius();
+ double area = static_cast<double>(neigh.GetSize()[0]*neigh.GetSize()[1]);
+ double areaInv = 1/area;
+ OffsetType offset;
+ offset.Fill(0);
+ OffsetType offsetOff;
+ OffsetType offsetOffInit;
+
+ offsetOffInit[0] = -radius[0]+this->GetOffset()[0]-1;
+ offsetOffInit[1] = -radius[1]+this->GetOffset()[1]-1;
+
+ int histoIdX = 0;
+ int histoIdY = 0;
+ double out = 0.;
+
+ IntVectorType histoTemp;
+ IntVectorVectorType histo;
+ if (binsLength[0] != 0)
+ histoTemp = IntVectorType( vcl_floor( static_cast<double>(this->GetMaxi()-this->GetMini())/binsLength[0])+1., 0);
+ else
+ histoTemp = IntVectorType( 1, 0 );
+
+ if (binsLength[1] != 0)
+ histo = IntVectorVectorType( vcl_floor(static_cast<double>(this->GetMaxiOff()-this->GetMiniOff())/binsLength[1])+1., histoTemp );
+ else
+ histo = IntVectorVectorType( 1, histoTemp );
+
+ offsetOff = offsetOffInit;
+ 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;
+ histoIdX = 0;
+ histoIdY = 0;
+ if ( binsLength[1] != 0)
+ histoIdX = static_cast<int>(vcl_floor( (static_cast<double>(neighOff[offsetOff])-this->GetMiniOff()) / static_cast<double>(binsLength[1]) ));
+ if ( binsLength[0] !=0 )
+ histoIdY = static_cast<int>(vcl_floor( (static_cast<double>(neigh[offset])-this->GetMini()) /static_cast<double>( binsLength[0]) ));
+
+ histo[histoIdX][histoIdY]++;
+
+ }
+ }
+ // loop over bin neighborhood values
+ for (unsigned sB = 0; sB<histo[0].size(); sB++)
+ {
+ double Px_y = 0.;
+ double nCeil = (static_cast<double>(sB)+0.5)*binsLength[0];
+ for (unsigned r = 0; r<histo.size(); r++)
+ {
+ double rVal = (static_cast<double>(r)+0.5)*binsLength[1];
+ for (unsigned s = 0; s<histo[r].size(); s++)
+ {
+ double sVal = (static_cast<double>(s)+0.5)*binsLength[0];
+ if( vcl_abs(rVal + sVal - nCeil) < vcl_abs(binsLength[0]) || vcl_abs(rVal + sVal - 2*nCeil) < vcl_abs(binsLength[0]) )
+ {
+ Px_y += static_cast<double>(histo[r][s])*areaInv;
+ }
+ }
+ }
+ if(Px_y != 0.)
+ out += Px_y * vcl_log(Px_y);
+ }
+
+
+ if(out != 0)
+ out = -out;
+
+
+ return out;
+ }
+
+};
+
+
+} // namespace Functor
+} // namespace otb
+
+#endif
+
diff --git a/Testing/Code/FeatureExtraction/CMakeLists.txt b/Testing/Code/FeatureExtraction/CMakeLists.txt
index 707964e0a8..0d802c625a 100644
--- a/Testing/Code/FeatureExtraction/CMakeLists.txt
+++ b/Testing/Code/FeatureExtraction/CMakeLists.txt
@@ -1187,6 +1187,21 @@ ADD_TEST(feTvDifferenceEntropyTextureFunctor ${FEATUREEXTRACTION_TESTS12}
2 # offset[1]
)
+# ------- otb::SumEntropyTextureFunctor -------------
+ADD_TEST(feTvSumEntropyTextureFunctor ${FEATUREEXTRACTION_TESTS12}
+--compare-image ${EPS}
+ ${BASELINE}/feTvSumEntropyTextureFunctor.tif
+ ${TEMP}/feTvSumEntropyTextureFunctor.tif
+ otbTextureFunctor
+ SEN #difference entropy
+ ${INPUTDATA}/poupees_sub.png
+ ${TEMP}/feTvSumEntropyTextureFunctor.tif
+ 5 # radius
+ -2 # offset[0]
+ 2 # offset[1]
+)
+
+
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ otbFeatureExtractionTests13 ~~~~~~~~~~~~~~~~~~~~~
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1328,6 +1343,23 @@ ADD_TEST(feTvDifferenceEntropyTextureImageFunction ${FEATUREEXTRACTION_TESTS13}
2 # offset[1]
)
+
+# ------- otb::SumEntropyTextureImageFunction -------------
+ADD_TEST(feTvSumEntropyTextureImageFunction ${FEATUREEXTRACTION_TESTS13}
+--compare-image ${EPS}
+ ${BASELINE}/feTuSumEntropyTextureImageFunction.tif
+ ${TEMP}/feTvSumEntropyTextureImageFunction.tif
+ otbTextureImageFunction
+ SEN #difference entropy
+ ${INPUTDATA}/poupees_sub_c1.png
+ ${TEMP}/feTvSumEntropyTextureImageFunction.tif
+ 5 # radius[0]
+ 5 # radius[1]
+ -2 # offset[0]
+ 2 # offset[1]
+ )
+
+
# A enrichir
SET(BasicFeatureExtraction_SRCS1
otbAlignImageToPath.cxx
diff --git a/Testing/Code/FeatureExtraction/otbTextureFunctor.cxx b/Testing/Code/FeatureExtraction/otbTextureFunctor.cxx
index 2ac85ea177..a2714c59f6 100644
--- a/Testing/Code/FeatureExtraction/otbTextureFunctor.cxx
+++ b/Testing/Code/FeatureExtraction/otbTextureFunctor.cxx
@@ -33,6 +33,9 @@
#include "otbContrastTextureFunctor.h"
#include "otbSumAverageTextureFunctor.h"
#include "otbDifferenceEntropyTextureFunctor.h"
+#include "otbSumEntropyTextureFunctor.h"
+
+
template<class TInputImage, class TOutputImage, class TFunctor>
int generic_TextureFunctor(int argc, char * argv[])
@@ -125,6 +128,11 @@ int otbTextureFunctor(int argc, char * argv[])
typedef otb::Functor::DifferenceEntropyTextureFunctor<IteratorType, IteratorType, PixelType> FunctorType;
return( generic_TextureFunctor<ImageType, ImageType, FunctorType>(argc,argv) );
}
+ else if ( strArgv == "SEN" )
+ {
+ typedef otb::Functor::SumEntropyTextureFunctor<IteratorType, IteratorType, PixelType> FunctorType;
+ return( generic_TextureFunctor<ImageType, ImageType, FunctorType>(argc,argv) );
+ }
else
{
return EXIT_FAILURE;
diff --git a/Testing/Code/FeatureExtraction/otbTextureImageFunction.cxx b/Testing/Code/FeatureExtraction/otbTextureImageFunction.cxx
index 45bc396d2a..052cdcabde 100644
--- a/Testing/Code/FeatureExtraction/otbTextureImageFunction.cxx
+++ b/Testing/Code/FeatureExtraction/otbTextureImageFunction.cxx
@@ -36,6 +36,7 @@
#include "otbContrastTextureFunctor.h"
#include "otbSumAverageTextureFunctor.h"
#include "otbDifferenceEntropyTextureFunctor.h"
+#include "otbSumEntropyTextureFunctor.h"
template<class TInputImage, class TOutputImage, class TFunctor>
@@ -98,7 +99,6 @@ int otbTextureImageFunction(int argc, char * argv[])
if(strArgv == "ENJ")
{
- std::cout<<"ENEGRY"<<std::endl;
typedef otb::Functor::EnergyTextureFunctor<IteratorType, IteratorType, VectorType> FunctorType;
return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
}
@@ -142,6 +142,11 @@ int otbTextureImageFunction(int argc, char * argv[])
typedef otb::Functor::DifferenceEntropyTextureFunctor<IteratorType, IteratorType, VectorType> FunctorType;
return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
}
+ else if ( strArgv == "SEN" )
+ {
+ typedef otb::Functor::SumEntropyTextureFunctor<IteratorType, IteratorType, VectorType> FunctorType;
+ return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
+ }
else
{
return EXIT_FAILURE;
--
GitLab