diff --git a/Code/FeatureExtraction/otbEntropyTextureFunctor.h b/Code/FeatureExtraction/otbEntropyTextureFunctor.h index f76cdd22036aadb2b44d90be99d775cb7ace44d1..d0e1fc64cd0b872026aa98ff021a2f6e7a9496f5 100644 --- a/Code/FeatureExtraction/otbEntropyTextureFunctor.h +++ b/Code/FeatureExtraction/otbEntropyTextureFunctor.h @@ -16,7 +16,7 @@ =========================================================================*/ #ifndef __otbEntropyTextureFunctor_h -#define __ootbEntropyTextureFunctor_h +#define __otbEntropyTextureFunctor_h #include "otbTextureFunctorBase.h" @@ -63,82 +63,87 @@ public: typedef std::vector<IntVectorType> IntVectorVectorType; - /** Computes the histogram bins using Scott formula, plus min/max. */ + double GetMaxi(){ return m_Maxi; }; + double GetMini(){ return m_Mini; }; + double GetMaxiOff(){ return m_MaxiOff; }; + double GetMiniOff(){ return m_MiniOff; }; + + /** Computes the histogram bins using Scott formula, plus min/max. */ DoubleVectorType StatComputation(const NeighborhoodType &neigh, const NeighborhoodType &neighOff) - { - DoubleVectorType output; - double binLength = 0; - double binLengthOff =0; - double mean = 0.; - double meanOff = 0.; - RadiusType radius = neigh.GetRadius(); - m_Mini = itk::NumericTraits<double>::max(); - m_MiniOff = itk::NumericTraits<double>::max(); - m_Maxi = itk::NumericTraits<double>::NonpositiveMin(); - m_MaxiOff = itk::NumericTraits<double>::NonpositiveMin(); - double area = static_cast<double>(neigh.GetSize()[0]*neigh.GetSize()[1]); - double areaInv = 1/area; - double scottCoef = 3.5 /(vcl_pow(area, 1/3) ); - - 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; - - 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; - mean += static_cast<double>(neigh[offset]); - meanOff += static_cast<double>(neighOff[offsetOff]); + { + DoubleVectorType output; + double binLength = 0; + double binLengthOff =0; + double mean = 0.; + double meanOff = 0.; + RadiusType radius = neigh.GetRadius(); + m_Mini = itk::NumericTraits<double>::max(); + m_MiniOff = itk::NumericTraits<double>::max(); + m_Maxi = itk::NumericTraits<double>::NonpositiveMin(); + m_MaxiOff = itk::NumericTraits<double>::NonpositiveMin(); + double area = static_cast<double>(neigh.GetSize()[0]*neigh.GetSize()[1]); + double areaInv = 1/area; + double scottCoef = 3.5 /(vcl_pow(area, 1/3) ); + + 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; + + 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; + mean += static_cast<double>(neigh[offset]); + meanOff += static_cast<double>(neighOff[offsetOff]); - m_Mini = std::min(static_cast<double>(neigh[offset]), m_Mini); - m_Maxi = std::max(static_cast<double>(neigh[offset]), m_Maxi); - m_MiniOff = std::min(static_cast<double>(neighOff[offsetOff]),m_MiniOff); - m_MaxiOff = std::max(static_cast<double>(neighOff[offsetOff]),m_MaxiOff); - } - } - mean *= areaInv; - meanOff *= areaInv; + m_Mini = std::min(static_cast<double>(neigh[offset]), m_Mini); + m_Maxi = std::max(static_cast<double>(neigh[offset]), m_Maxi); + m_MiniOff = std::min(static_cast<double>(neighOff[offsetOff]),m_MiniOff); + m_MaxiOff = std::max(static_cast<double>(neighOff[offsetOff]),m_MaxiOff); + } + } + mean *= areaInv; + meanOff *= areaInv; - offsetOff = offsetOffInit; + 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; - binLength += vcl_pow( (mean-static_cast<double>(neigh[offset])), 2); - binLengthOff += vcl_pow( (meanOff-static_cast<double>(neighOff[offsetOff])), 2); - } - } + 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; + binLength += vcl_pow( (mean-static_cast<double>(neigh[offset])), 2); + binLengthOff += vcl_pow( (meanOff-static_cast<double>(neighOff[offsetOff])), 2); + } + } - binLength *= areaInv; - binLength = vcl_sqrt( binLength ); - binLengthOff *= areaInv; - binLengthOff = vcl_sqrt( binLengthOff ); + binLength *= areaInv; + binLength = vcl_sqrt( binLength ); + binLengthOff *= areaInv; + binLengthOff = vcl_sqrt( binLengthOff ); - output.push_back( scottCoef*binLength ); - output.push_back( scottCoef*binLengthOff ); + output.push_back( scottCoef*binLength ); + output.push_back( scottCoef*binLengthOff ); - return output; - } + return output; + } - double ComputeOverSingleChannel(const NeighborhoodType &neigh, const NeighborhoodType &neighOff) + virtual double ComputeOverSingleChannel(const NeighborhoodType &neigh, const NeighborhoodType &neighOff) { DoubleVectorType binsLength = this->StatComputation(neigh, neighOff); @@ -203,74 +208,6 @@ public: out = -(out); return out; - - /* - IntVectorVectorType binsLength = this->StatComputation(it, itOff); - - RadiusType radius = it.GetRadius(); - 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; - - int histoIdX = 0; - int histoIdY = 0; - - //for ( unsigned int i=0; i<nbComp; i++ ) - //{ - IntVectorType histoTemp; - IntVectorVectorType histo; - if (binsLength[0][i] != 0) - histoTemp = IntVectorType( vcl_floor( static_cast<double>(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(static_cast<double>(m_MaxiOff[i]-m_MiniOff[i])/binsLength[1][i])+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; - if ( binsLength[1][i] != 0) - histoIdX = static_cast<int>(vcl_floor( (itOff.GetPixel(offsetOff)[i]-m_MiniOff[i]) / binsLength[1][i] )); - else - histoIdX = 0; - if ( binsLength[0][i] !=0 ) - histoIdY = static_cast<int>(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)); - } - } - - //} - - return outPix; -*/ } diff --git a/Code/FeatureExtraction/otbInverseDifferenceMomentTextureFunctor.h b/Code/FeatureExtraction/otbInverseDifferenceMomentTextureFunctor.h new file mode 100644 index 0000000000000000000000000000000000000000..0f83d2172065f96ddc7e67b5f00d08656a04d9ee --- /dev/null +++ b/Code/FeatureExtraction/otbInverseDifferenceMomentTextureFunctor.h @@ -0,0 +1,137 @@ +/*========================================================================= + + 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 __otbInverseDifferenceMomentTextureFunctor_h +#define __otbInverseDifferenceMomentTextureFunctor_h + +#include "otbEntropyTextureFunctor.h" + +namespace otb +{ +namespace Functor +{ +/** \class InverseDifferenceMomentTextureFunctor + * \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. + * InverseDifferenceMoment 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 InverseDifferenceMomentTextureFunctor : +public EntropyTextureFunctor<TIterInput1, TIterInput2, TOutput> +{ +public: + InverseDifferenceMomentTextureFunctor(){}; + ~InverseDifferenceMomentTextureFunctor(){}; + + 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; + typedef EntropyTextureFunctor<TIterInput1, TIterInput2, TOutput> Superclass; + + + virtual double ComputeOverSingleChannel(const NeighborhoodType &neigh, const NeighborhoodType &neighOff) + { + DoubleVectorType binsLength = Superclass::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]++; + + } + } + + for (unsigned r = 0; r<histo.size(); r++) + { + for (unsigned s = 0; s<histo[r].size(); s++) + { + double p = histo[r][s] * areaInv; + double dist = vcl_pow( ( (static_cast<double>(r)+0.5)*binsLength[1])-((static_cast<double>(s)+0.5)*binsLength[0]), 2); + if (p != 0) + out += ((1/(1+dist)) * p); + } + } + return out; + } + +}; + + + + + } // namespace Functor +} // namespace otb + +#endif + diff --git a/Testing/Code/FeatureExtraction/CMakeLists.txt b/Testing/Code/FeatureExtraction/CMakeLists.txt index 12dac3464b6b877959f5877240dec46534c72494..e5b76257daaa04d54935bb5b2b10545da5de3f93 100644 --- a/Testing/Code/FeatureExtraction/CMakeLists.txt +++ b/Testing/Code/FeatureExtraction/CMakeLists.txt @@ -1045,7 +1045,7 @@ ADD_TEST(feTvEnergyTextureFunctor ${FEATUREEXTRACTION_TESTS11} ${BASELINE}/feTvEnergyTextureFunctor.tif ${TEMP}/feTvEnergyTextureFunctor.tif otbEnergyTextureFunctor - ${INPUTDATA}/poupees.tif + ${INPUTDATA}/poupees_sub.png ${TEMP}/feTvEnergyTextureFunctor.tif 5 # radius -2 # offset[0] @@ -1058,13 +1058,27 @@ ADD_TEST(feTvEntropyTextureFunctor ${FEATUREEXTRACTION_TESTS11} ${BASELINE}/feTvEntropyTextureFunctor.tif ${TEMP}/feTvEntropyTextureFunctor.tif otbEntropyTextureFunctor - ${INPUTDATA}/poupees.tif + ${INPUTDATA}/poupees_sub.png ${TEMP}/feTvEntropyTextureFunctor.tif 5 # radius -2 # offset[0] 2 # offset[1] ) +# ------- otb::InverseDifferenceMomentTextureFunctor ------------- +ADD_TEST(feTvInverseDifferenceMomentTextureFunctor ${FEATUREEXTRACTION_TESTS11} +--compare-image ${EPS} + ${BASELINE}/feTvInverseDifferenceMomentTextureFunctor.tif + ${TEMP}/feTvInverseDifferenceMomentTextureFunctor.tif + otbInverseDifferenceMomentTextureFunctor + ${INPUTDATA}/poupees_sub.png + ${TEMP}/feTvInverseDifferenceMomentTextureFunctor.tif + 5 # radius + -2 # offset[0] + 2 # offset[1] +) + + # ------- otb::EnergyTextureImageFunction ------------- ADD_TEST(feTvEnergyTextureImageFunction ${FEATUREEXTRACTION_TESTS11} --compare-image ${EPS} @@ -1072,7 +1086,7 @@ ADD_TEST(feTvEnergyTextureImageFunction ${FEATUREEXTRACTION_TESTS11} ${TEMP}/feTvEnergyTextureImageFunction.tif otbTextureImageFunction ENJ #energy - ${INPUTDATA}/poupees_1canal.c1.hdr + ${INPUTDATA}/poupees_sub_c1.png ${TEMP}/feTvEnergyTextureImageFunction.tif 5 # radius[0] 5 # radius[1] @@ -1088,7 +1102,7 @@ ADD_TEST(feTvEntropyTextureImageFunction ${FEATUREEXTRACTION_TESTS11} ${TEMP}/feTvEntropyTextureImageFunction.tif otbTextureImageFunction ENT #entropy - ${INPUTDATA}/poupees_1canal.c1.hdr + ${INPUTDATA}/poupees_sub_c1.png ${TEMP}/feTvEntropyTextureImageFunction.tif 5 # radius[0] 5 # radius[1] @@ -1096,6 +1110,23 @@ ADD_TEST(feTvEntropyTextureImageFunction ${FEATUREEXTRACTION_TESTS11} 2 # offset[1] ) + +# ------- otb::InverseDifferenceMomentTextureImageFunction ------------- +ADD_TEST(feTvInverseDifferenceMomentTextureImageFunction ${FEATUREEXTRACTION_TESTS11} +--compare-image ${EPS} + ${BASELINE}/feTvInverseDifferenceMomentTextureImageFunction.tif + ${TEMP}/feTvInverseDifferenceMomentTextureImageFunction.tif + otbTextureImageFunction + IMD #inverse difference momentent + ${INPUTDATA}/poupees_sub_c1.png + ${TEMP}/feTvInverseDifferenceMomentTextureImageFunction.tif + 5 # radius[0] + 5 # radius[1] + -2 # offset[0] + 2 # offset[1] +) + + # A enrichir SET(BasicFeatureExtraction_SRCS1 otbAlignImageToPath.cxx @@ -1228,6 +1259,7 @@ otbCloudDetectionFilter.cxx otbSimplifyManyPathListFilter.cxx otbEnergyTextureFunctor.cxx otbEntropyTextureFunctor.cxx +otbInverseDifferenceMomentTextureFunctor.cxx otbTextureImageFunction.cxx ) diff --git a/Testing/Code/FeatureExtraction/otbFeatureExtractionTests11.cxx b/Testing/Code/FeatureExtraction/otbFeatureExtractionTests11.cxx index e3cd5c1747d8a54c2cba9d18d926c8df43c7d79c..10a02609091d76f6e007ca0ce9e9a303862ac993 100644 --- a/Testing/Code/FeatureExtraction/otbFeatureExtractionTests11.cxx +++ b/Testing/Code/FeatureExtraction/otbFeatureExtractionTests11.cxx @@ -35,5 +35,6 @@ REGISTER_TEST(otbCloudDetectionFilter); REGISTER_TEST(otbSimplifyManyPathListFilter); REGISTER_TEST(otbEnergyTextureFunctor); REGISTER_TEST(otbEntropyTextureFunctor); +REGISTER_TEST(otbInverseDifferenceMomentTextureFunctor); REGISTER_TEST(otbTextureImageFunction); } diff --git a/Testing/Code/FeatureExtraction/otbInverseDifferenceMomentTextureFunctor.cxx b/Testing/Code/FeatureExtraction/otbInverseDifferenceMomentTextureFunctor.cxx new file mode 100644 index 0000000000000000000000000000000000000000..ea15dc4d902f383e348e209b9458e10cb106dced --- /dev/null +++ b/Testing/Code/FeatureExtraction/otbInverseDifferenceMomentTextureFunctor.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 "otbImageFileWriter.h" +#include "otbInverseDifferenceMomentTextureFunctor.h" + + +int otbInverseDifferenceMomentTextureFunctor(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::ImageFileWriter<ImageType> WriterType; + + typedef itk::ConstNeighborhoodIterator<ImageType> IterType;; + typedef otb::Functor::InverseDifferenceMomentTextureFunctor<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/otbTextureImageFunction.cxx b/Testing/Code/FeatureExtraction/otbTextureImageFunction.cxx index c751a4aff750cb4de13ce662f063fb80487bf4a3..50f3422bd2f2ce05cd6b4c8217bff8bac06a9f7e 100644 --- a/Testing/Code/FeatureExtraction/otbTextureImageFunction.cxx +++ b/Testing/Code/FeatureExtraction/otbTextureImageFunction.cxx @@ -29,6 +29,8 @@ // Functors list #include "otbEnergyTextureFunctor.h" #include "otbEntropyTextureFunctor.h" +#include "otbInverseDifferenceMomentTextureFunctor.h" + template<class TInputImage, class TOutputImage, class TFunctor> int generic_TextureImageFunction(int argc, char * argv[]) @@ -99,6 +101,11 @@ int otbTextureImageFunction(int argc, char * argv[]) typedef otb::Functor::EntropyTextureFunctor<IteratorType, IteratorType, VectorType> FunctorType; return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) ); } + else if ( strArgv == "IMD" ) + { + typedef otb::Functor::InverseDifferenceMomentTextureFunctor<IteratorType, IteratorType, VectorType> FunctorType; + return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) ); + } else { return EXIT_FAILURE;