diff --git a/Code/FeatureExtraction/otbContrastTextureFunctor.h b/Code/FeatureExtraction/otbContrastTextureFunctor.h index 8e508bc525b542f77a4ba2abb75265d3d915e23a..ad1d97300d9ce6cdb309321339fe2e55ebc48688 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 bf667f36e22bb7374ce1cab45ce7c0202929c8c4..d984f13e173f22babbb83429f1c3f37ad91e6b9c 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 c9028a7b273677322805bd2f6c4c412f18b36a5e..3de106d2f793dd0c425f7ec8d673b0e7066702c5 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 0000000000000000000000000000000000000000..dad317125a8372d09b9e9b02c4ffa3117cfe4ed8 --- /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 707964e0a8a05902f96a48b47fce1e2ad6c8272d..0d802c625a825129ca9d67dc66163ffec4a04ae2 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 2ac85ea17750682eba037522bbd6da7fb7d004cb..a2714c59f653b05c2f353f46f7cf9f55f56493dd 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 45bc396d2a1d624a559f51e9c53b5af598b1184f..052cdcabdea520dd463481dc465dab6e75d40c17 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;