diff --git a/Code/BasicFilters/otbFunctionWithNeighborhoodToImageFilter.txx b/Code/BasicFilters/otbFunctionWithNeighborhoodToImageFilter.txx index 5fd3b39ebef082b5994f5fcc309de34c37482036..299c4a8b3cc737db826a15170ef5af937aa31855 100644 --- a/Code/BasicFilters/otbFunctionWithNeighborhoodToImageFilter.txx +++ b/Code/BasicFilters/otbFunctionWithNeighborhoodToImageFilter.txx @@ -49,15 +49,6 @@ FunctionWithNeighborhoodToImageFilter<TInputImage,TOutputImage,TFunction> m_Radius = this->GetFunction()->GetRadius(); 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 > @@ -137,7 +128,6 @@ FunctionWithNeighborhoodToImageFilter<TInputImage,TOutputImage,TFunction> while ( !inputIt.IsAtEnd() ) { - //std::cout<<"MOTHER "<<threadId<<" : "<<inputIt.GetIndex()<<" "<<inputIt.Get()<<std::endl; outputIt.Set( static_cast<OutputImagePixelType>(this->GetFunction()->EvaluateAtIndex(inputIt.GetIndex())) ); ++inputIt; ++outputIt; diff --git a/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.txx b/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.txx index 7a2d52092baa3ec69e629134da72c8ce3e27bad3..4476bad01421fcf8aa5f0b5fe00172e8296a4589 100644 --- a/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.txx +++ b/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.txx @@ -119,7 +119,6 @@ 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 @@ -185,7 +184,6 @@ 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 95c49a17f6676decb96c362496053f6926c08178..b6990c60f2fac34c642fba933a4955f97919d288 100644 --- a/Code/FeatureExtraction/otbEnergyTextureFunctor.h +++ b/Code/FeatureExtraction/otbEnergyTextureFunctor.h @@ -16,9 +16,11 @@ =========================================================================*/ #ifndef __otbEnergyTextureFunctor_h -#define __ootbEnergyTextureFunctor_h +#define __otbEnergyTextureFunctor_h + + +#include "otbTextureFunctorBase.h" -#include "otbMath.h" namespace otb { @@ -29,94 +31,67 @@ namespace Functor * * 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 TextureFunctorBase * \ingroup Statistics */ template <class TIterInput1, class TIterInput2, class TOutput> -class EnergyTextureFunctor +class ITK_EXPORT EnergyTextureFunctor : +public TextureFunctorBase<TIterInput1, TIterInput2, TOutput> { public: - EnergyTextureFunctor() - { - m_Offset.Fill(1); - }; - ~EnergyTextureFunctor() {}; - - typedef TIterInput1 IterType1; - typedef TIterInput2 IterType2; - typedef TOutput OutputType; - typedef typename IterType1::OffsetType OffsetType; - typedef typename IterType1::RadiusType RadiusType; - typedef typename OutputType::ValueType OutputPixelType; - - void SetOffset(OffsetType off) + EnergyTextureFunctor(){}; + ~EnergyTextureFunctor(){}; + + 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; + + double ComputeOverSingleChannel(const NeighborhoodType &neigh, const NeighborhoodType &neighOff) { - m_Offset=off; - }; - OffsetType GetOffset() - { - return m_Offset; - }; - - inline TOutput operator()(const TIterInput1 &it, const TIterInput2 &itOff) - { - RadiusType radius = it.GetRadius(); - unsigned int nbComp = it.GetCenterPixel().GetSize(); - double area = static_cast<double>(radius[0]*radius[1]); - OutputType outPix; - outPix.SetSize( nbComp ); - outPix.Fill(0); - + RadiusType radius = neigh.GetRadius(); + double area = static_cast<double>(neigh.GetSize()[0]*neigh.GetSize()[1]); 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; - - //std::cout<<offsetOffInit<<std::endl; + offsetOffInit[0] = -radius[0]+this->GetOffset()[0]-1; + offsetOffInit[1] = -radius[1]+this->GetOffset()[1]-1; double temp = 0.; - int ll = 0; double norm = 0.; - for ( unsigned int i=0; i<nbComp; i++ ) - { - offsetOff = offsetOffInit; - temp = 0.; - for ( int l = -static_cast<int>(radius[0]); l <= static_cast<int>(radius[0]); l++ ) + double output = 0.; + + 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; - ll = l*l; for ( int k = -static_cast<int>(radius[1]); k <= static_cast<int>(radius[1]); k++) - { - offsetOff[1]++; - offset[1] = k; - norm = vcl_pow(static_cast<double>(itOff.GetPixel(offsetOff)[i]-itOff.GetCenterPixel()[i]), 2); - temp += norm; - } - temp /= area; - outPix[i] = static_cast<OutputPixelType>( vcl_pow(temp, 2) ); + { + offsetOff[1]++; + offset[1] = k; + norm = vcl_pow(static_cast<double>( ( neigh[offset] - neighOff[neighOff.GetCenterNeighborhoodIndex()] ) ), 2); + temp += norm; + } } - } - - - return outPix; + temp /= area; + return vcl_pow(temp, 2); } -private: - OffsetType m_Offset; - }; - - } // namespace Functor } // namespace otb diff --git a/Code/FeatureExtraction/otbEnergyTextureImageFunction.txx b/Code/FeatureExtraction/otbEnergyTextureImageFunction.txx index 65f77cc36983be6b047aeb435c0e75dd82471f0d..a7df15db8ac49c8052f4ba2687df61be62fbd2e3 100644 --- a/Code/FeatureExtraction/otbEnergyTextureImageFunction.txx +++ b/Code/FeatureExtraction/otbEnergyTextureImageFunction.txx @@ -21,7 +21,7 @@ #include "otbEnergyTextureImageFunction.h" #include "itkConstNeighborhoodIterator.h" #include "otbEnergyTextureFunctor.h" - +#include "itkVariableLengthVector.h" namespace otb { @@ -72,6 +72,7 @@ EnergyTextureImageFunction<TInputImage,TCoordRep> } typedef itk::ConstNeighborhoodIterator<InputImageType> IterType; + typedef typename IterType::NeighborhoodType NeighborhoodType; IterType it(m_Radius, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion()); it.SetLocation(index); SizeType radiusOff; @@ -80,52 +81,13 @@ EnergyTextureImageFunction<TInputImage,TCoordRep> IterType itOff(radiusOff, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion()); itOff.SetLocation(index); - /* - typedef typename Functor::EnergyTextureFunctor<IterType, IterType, RealType> FunctType; + // false third template... : double has no type ValueType but no impact + typedef itk::VariableLengthVector<double> VecType; + typedef typename Functor::EnergyTextureFunctor<IterType, IterType, VecType/*RealType*/> FunctType; FunctType funct; - funct.SetRadius(m_Radius); funct.SetOffset(m_Offset); - funct(it, itOff); - //FunctType::operator()(it, itOff); - - return 0; - */ - double area = static_cast<double>(m_Radius[0]*m_Radius[1]); - RealType output; - - OffsetType offset; - offset.Fill(0); - OffsetType offsetOff; - OffsetType offsetOffInit; - - offsetOffInit[0] = -m_Radius[0]+m_Offset[0]-1; //++ in for - offsetOffInit[1] = -m_Radius[1]+m_Offset[1]-1; //++ in for - - double temp = 0.; - int ll = 0; - double norm = 0.; - offsetOff = offsetOffInit; - temp = 0.; - for ( int l = -static_cast<int>(m_Radius[0]); l <= static_cast<int>(m_Radius[0]); l++ ) - { - offsetOff[0]++; - offset[0] = l; - ll = l*l; - offsetOff[1] = offsetOffInit[1]; - for ( int k = -static_cast<int>(m_Radius[1]); k <= static_cast<int>(m_Radius[1]); k++) - { - offsetOff[1]++; - offset[1] = k; - norm = vcl_pow(static_cast<double>(itOff.GetPixel(offsetOff)-itOff.GetCenterPixel()), 2); - temp += norm; - } - temp /= area; - output = static_cast<RealType>( vcl_pow(temp, 2) ); - } - return output; - - - + + return static_cast<RealType>(funct.ComputeOverSingleChannel( it.GetNeighborhood(), itOff.GetNeighborhood()) ); } diff --git a/Code/FeatureExtraction/otbEntropyTextureFunctor.h b/Code/FeatureExtraction/otbEntropyTextureFunctor.h index c0f9e4381e97ab3143adbb310890da86484f41a7..8b78d0320e065b35381beb3b279770cb26ef6907 100644 --- a/Code/FeatureExtraction/otbEntropyTextureFunctor.h +++ b/Code/FeatureExtraction/otbEntropyTextureFunctor.h @@ -18,7 +18,7 @@ #ifndef __otbEntropyTextureFunctor_h #define __ootbEntropyTextureFunctor_h -#include "otbMath.h" +#include "otbTextureFunctorBase.h" namespace otb { @@ -37,47 +37,46 @@ namespace Functor * \ingroup Statistics */ template <class TIterInput1, class TIterInput2, class TOutput> -class EntropyTextureFunctor +class ITK_EXPORT EntropyTextureFunctor : +public TextureFunctorBase<TIterInput1, TIterInput2, TOutput> { public: EntropyTextureFunctor() - { - m_Offset.Fill(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; - }; + { + m_Mini = itk::NumericTraits<double>::max(); + m_MiniOff = itk::NumericTraits<double>::max(); + m_Maxi = itk::NumericTraits<double>::NonpositiveMin(); + m_MaxiOff = itk::NumericTraits<double>::NonpositiveMin(); + }; + ~EntropyTextureFunctor(){}; + + 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; + /** Computes the histogram bins using Scott formula, plus min/max. */ - IntVectorVectorType StatComputation(const TIterInput1 &it, const TIterInput2 &itOff) + DoubleVectorType StatComputation(const NeighborhoodType &neigh, const NeighborhoodType &neighOff) { - 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]); + 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) ); @@ -86,85 +85,130 @@ public: OffsetType offsetOff; OffsetType offsetOffInit; - offsetOffInit[0] = -radius[0]+m_Offset[0]-1; - offsetOffInit[1] = -radius[1]+m_Offset[1]-1; - DoubleVectorType mean; - DoubleVectorType 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++ ) + 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; - 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); + { + 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; + + 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); + } } - } + + binLength *= areaInv; + binLength = vcl_sqrt( binLength ); + binLengthOff *= areaInv; + binLengthOff = vcl_sqrt( binLengthOff ); + + output.push_back( scottCoef*binLength ); + output.push_back( scottCoef*binLengthOff ); - 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++ ) + return output; + } + + + 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>(m_Maxi-m_Mini)/binsLength[0])+1., 0); + else + histoTemp = IntVectorType( 1, 0 ); + if (binsLength[1] != 0) + histo = IntVectorVectorType( vcl_floor(static_cast<double>(m_MaxiOff-m_MiniOff)/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]; + offsetOff[1] = offsetOffInit[1]; offset[0] = l; - for( int k = -static_cast<int>(radius[1]); k <= static_cast<int>(radius[1]); k++) + 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); + histoIdX = 0; + histoIdY = 0; + if ( binsLength[1] != 0) + histoIdX = static_cast<int>(vcl_floor( (static_cast<double>(neighOff[offsetOff])-m_MiniOff) / static_cast<double>(binsLength[1]) )); + if ( binsLength[0] !=0 ) + histoIdY = static_cast<int>(vcl_floor( (static_cast<double>(neigh[offset])-m_Mini) /static_cast<double>( binsLength[0]) )); + + histo[histoIdX][histoIdY]++; + } } - - stdCh *= areaInv; - stdCh = vcl_sqrt( stdCh ); - stdChOff *= areaInv; - stdChOff = vcl_sqrt( stdChOff ); - binLength.push_back( scottCoef*stdCh ); - binLengthOff.push_back( scottCoef*stdChOff ); - - } - output.push_back(binLength); - output.push_back(binLengthOff); - - return output; - } + for (unsigned r = 0; r<histo.size(); r++) + { + for (unsigned s = 0; s<histo[r].size(); s++) + { + double p = histo[r][s] * areaInv; + if (p != 0) + out += (p * vcl_log(p)); + } + } + + return out; - inline TOutput operator()(const TIterInput1 &it, const TIterInput2 &itOff) - { + /* 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; @@ -176,8 +220,8 @@ public: int histoIdX = 0; int histoIdY = 0; - for ( unsigned int i=0; i<nbComp; i++ ) - { + //for ( unsigned int i=0; i<nbComp; i++ ) + //{ IntVectorType histoTemp; IntVectorVectorType histo; if (binsLength[0][i] != 0) @@ -221,20 +265,21 @@ public: } } - } + //} return outPix; - +*/ } + + private: - OffsetType m_Offset; - /** Stores min/max neighborhood area values */ - DoubleVectorType m_Mini; - DoubleVectorType m_Maxi; - /** Stores min/max neighborhood+offset values */ - DoubleVectorType m_MiniOff; - DoubleVectorType m_MaxiOff; + /** Stores min/max neighborhood area values */ + double m_Mini; + double m_Maxi; + /** Stores min/max neighborhood+offset values */ + double m_MiniOff; + double m_MaxiOff; }; diff --git a/Code/FeatureExtraction/otbEntropyTextureImageFunction.h b/Code/FeatureExtraction/otbEntropyTextureImageFunction.h new file mode 100644 index 0000000000000000000000000000000000000000..d57100306d5b8f5a77bcabd9950f8a9fbf89457f --- /dev/null +++ b/Code/FeatureExtraction/otbEntropyTextureImageFunction.h @@ -0,0 +1,118 @@ +/*========================================================================= + + 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 __otbEntropyTextureImageFunction_h +#define __otbEntropyTextureImageFunction_h + +#include "itkImageFunction.h" +#include "itkNumericTraits.h" + + +namespace otb +{ + +/** + * \class EntropyImageFunction + * \brief Calculate the energy in the neighborhood of a pixel + * + * This class is templated over the input image type and the + * coordinate representation type (e.g. float or double ). + * + * \ingroup ImageFunctions + */ +template <class TInputImage, class TCoordRep = float > +class ITK_EXPORT EntropyTextureImageFunction : + public itk::ImageFunction< TInputImage, ITK_TYPENAME itk::NumericTraits<typename TInputImage::PixelType>::RealType, TCoordRep > +{ + public: + /** Standard class typedefs. */ + typedef EntropyTextureImageFunction Self; + typedef itk::ImageFunction<TInputImage, ITK_TYPENAME itk::NumericTraits<typename TInputImage::PixelType>::RealType, + TCoordRep > Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + + /** Run-time type information (and related methods). */ + itkTypeMacro(EntropyTextureImageFunction, itk::ImageFunction); + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** typedef support. */ + typedef TInputImage InputImageType; + typedef typename InputImageType::OffsetType OffsetType; + typedef typename InputImageType::SizeType SizeType; + typedef typename InputImageType::PixelType PixelType; + typedef typename Superclass::OutputType OutputType; + typedef typename Superclass::IndexType IndexType; + typedef typename Superclass::ContinuousIndexType ContinuousIndexType; + typedef typename Superclass::PointType PointType; + typedef typename itk::NumericTraits<typename InputImageType::PixelType>::RealType RealType; + + /** Dimension of the underlying image. */ + itkStaticConstMacro(ImageDimension, unsigned int,InputImageType::ImageDimension); + + + /** Evalulate the function at specified index */ + virtual RealType EvaluateAtIndex( const IndexType& index ) const; + + /** Evaluate the function at non-integer positions */ + virtual RealType Evaluate( const PointType& point ) const + { + IndexType index; + this->ConvertPointToNearestIndex( point, index ); + return this->EvaluateAtIndex( index ); + } + virtual RealType EvaluateAtContinuousIndex( + const ContinuousIndexType& cindex ) const + { + IndexType index; + this->ConvertContinuousIndexToNearestIndex( cindex, index ); + return this->EvaluateAtIndex( index ) ; + } + + /** Get/Set the radius of the neighborhood over which the + statistics are evaluated */ + itkSetMacro( Radius, SizeType); + itkGetMacro( Radius, SizeType); + itkSetMacro( Offset, OffsetType); + itkGetMacro( Offset, OffsetType ); + + protected: + EntropyTextureImageFunction(); + ~EntropyTextureImageFunction() {}; + void PrintSelf(std::ostream& os, itk::Indent indent) const; + +private: + EntropyTextureImageFunction( const Self& ); //purposely not implemented + void operator=( const Self& ); //purposely not implemented + + SizeType m_Radius; + OffsetType m_Offset; + +}; + +} // end namespace otb + + + +#ifndef OTB_MANUAL_INSTANTIATION +# include "otbEntropyTextureImageFunction.txx" +#endif + +#endif + diff --git a/Code/FeatureExtraction/otbEntropyTextureImageFunction.txx b/Code/FeatureExtraction/otbEntropyTextureImageFunction.txx new file mode 100644 index 0000000000000000000000000000000000000000..4d24a06bc8ae34f37ffc454b61d5c4a920c7c6f3 --- /dev/null +++ b/Code/FeatureExtraction/otbEntropyTextureImageFunction.txx @@ -0,0 +1,97 @@ +/*========================================================================= + + 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 __otbEntropyTextureImageFunction_txx +#define __otbEntropyTextureImageFunction_txx + +#include "otbEntropyTextureImageFunction.h" +#include "itkConstNeighborhoodIterator.h" +#include "otbEntropyTextureFunctor.h" +#include "itkVariableLengthVector.h" + +namespace otb +{ + +/** + * Constructor + */ +template <class TInputImage, class TCoordRep> +EntropyTextureImageFunction<TInputImage,TCoordRep> +::EntropyTextureImageFunction() +{ + m_Radius.Fill(0); + m_Offset.Fill(0); +} + + +/** + * + */ +template <class TInputImage, class TCoordRep> +void +EntropyTextureImageFunction<TInputImage,TCoordRep> +::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + this->Superclass::PrintSelf(os,indent); + os << indent << "Radius: " << m_Radius << std::endl; + os << indent << "Offset: " << m_Offset << std::endl; +} + + +/** + * + */ +template <class TInputImage, class TCoordRep> +typename EntropyTextureImageFunction<TInputImage,TCoordRep> +::RealType +EntropyTextureImageFunction<TInputImage,TCoordRep> +::EvaluateAtIndex(const IndexType& index) const +{ + if ( !this->GetInputImage() ) + { + return ( itk::NumericTraits<RealType>::max() ); + } + + if ( !this->IsInsideBuffer( index ) ) + { + return ( itk::NumericTraits<RealType>::max() ); + } + + typedef itk::ConstNeighborhoodIterator<InputImageType> IterType; + typedef typename IterType::NeighborhoodType NeighborhoodType; + IterType it(m_Radius, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion()); + it.SetLocation(index); + SizeType radiusOff; + radiusOff[0] = m_Radius[0] + vcl_abs(m_Offset[0]); + radiusOff[1] = m_Radius[1] + vcl_abs(m_Offset[1]); + IterType itOff(radiusOff, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion()); + itOff.SetLocation(index); + + // false third template.. but no impact + typedef itk::VariableLengthVector<double> VecType; + typedef typename Functor::EntropyTextureFunctor<IterType, IterType, VecType /*RealType*/> FunctType; + FunctType funct; + funct.SetOffset(m_Offset); + + return static_cast<RealType>(funct.ComputeOverSingleChannel( it.GetNeighborhood(), itOff.GetNeighborhood()) ); + +} + + +} // end namespace otb + +#endif diff --git a/Code/FeatureExtraction/otbTextureFunctorBase.h b/Code/FeatureExtraction/otbTextureFunctorBase.h new file mode 100644 index 0000000000000000000000000000000000000000..19980ad08750f2e9fdfd257d7d242405e9d7ff82 --- /dev/null +++ b/Code/FeatureExtraction/otbTextureFunctorBase.h @@ -0,0 +1,124 @@ +/*========================================================================= + + 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 __otbTextureFunctorBase_h +#define __otbTextureFunctorBase_h + +#include "otbMath.h" +#include "itkNeighborhood.h" + +namespace otb +{ +namespace Functor +{ +/** \class TextureFunctorBase + * \brief This functor is the basic class for texture computation. + * + * Implements the offset setting and the computation loop over each channel + * TIterInput1 and TIterInput2 are an iterators. The first one for the neighborhood, the second for the considered offset. + * TOutput is a PixelType. + * + * \ingroup Functor + * \ingroup Statistics + */ + +template <class TIterInput1, class TIterInput2, class TOutput> +class TextureFunctorBase +{ +public: + TextureFunctorBase() + { + m_Offset.Fill(1); + }; + ~TextureFunctorBase() {}; + + 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 typename IterType1::InternalPixelType InternalPixelType; + typedef typename IterType1::ImageType ImageType; + typedef itk::Neighborhood<InternalPixelType,::itk::GetImageDimension<ImageType>::ImageDimension> NeighborhoodType; + + void SetOffset(OffsetType off) + { + m_Offset=off; + }; + OffsetType GetOffset() + { + return m_Offset; + }; + + + inline TOutput operator()(const IterType1 &it, const IterType2 &itOff) + { + RadiusType radius = it.GetRadius(); + RadiusType radiusOff = itOff.GetRadius(); + OutputType outPix; + outPix.SetSize( it.GetCenterPixel().GetSize() ); + outPix.Fill(0); + OffsetType offset; + offset.Fill(0); + + for ( unsigned int i=0; i<it.GetCenterPixel().GetSize(); i++ ) + { + NeighborhoodType inNeigh; + inNeigh.SetRadius(it.GetRadius()); + NeighborhoodType offNeigh; + offNeigh.SetRadius(itOff.GetRadius()); + + for ( int l = -static_cast<int>(radius[0]); l <= static_cast<int>(radius[0]); l++ ) + { + offset[0] = l; + for ( int k = -static_cast<int>(radius[1]); k <= static_cast<int>(radius[1]); k++) + { + offset[1] = k; + inNeigh[offset] = it.GetPixel(offset)[i]; + } + } + offset.Fill(0); + for ( int l = -static_cast<int>(radiusOff[0]); l <= static_cast<int>(radiusOff[0]); l++ ) + { + offset[0] = l; + for ( int k = -static_cast<int>(radiusOff[1]); k <= static_cast<int>(radiusOff[1]); k++) + { + offset[1] = k; + offNeigh[offset] = itOff.GetPixel(offset)[i]; + } + } + outPix[i] = static_cast<OutputPixelType>( this->ComputeOverSingleChannel(inNeigh, offNeigh) ); + } + return outPix; + } + + virtual double ComputeOverSingleChannel(const NeighborhoodType &neigh, const NeighborhoodType &neighOff) = 0; + + private: + OffsetType m_Offset; + +}; + + + + +} // namespace Functor +} // namespace otb + +#endif + diff --git a/Testing/Code/FeatureExtraction/CMakeLists.txt b/Testing/Code/FeatureExtraction/CMakeLists.txt index 3262c5858e5a09e9df19bfc8651a4f383d2a33a5..2f728ba6f85c7b514a3ab089e573672123c32ad1 100644 --- a/Testing/Code/FeatureExtraction/CMakeLists.txt +++ b/Testing/Code/FeatureExtraction/CMakeLists.txt @@ -1082,6 +1082,24 @@ ADD_TEST(feTvEnergyTextureImageFunction ${FEATUREEXTRACTION_TESTS11} 2 # offset[1] ) + +# ------- otb::EntropyTextureImageFunction ------------- +ADD_TEST(feTuEntropyTextureImageFunctionNew ${FEATUREEXTRACTION_TESTS11} + otbEntropyTextureImageFunctionNew +) +ADD_TEST(feTvEntropyTextureImageFunction ${FEATUREEXTRACTION_TESTS11} +--compare-image ${EPS} + ${BASELINE}/feTvEntropyTextureImageFunction.tif + ${TEMP}/feTvEntropyTextureImageFunction.tif + otbEntropyTextureImageFunction + ${INPUTDATA}/poupees_1canal.c1.hdr + ${TEMP}/feTvEntropyTextureImageFunction.tif + 5 # radius[0] + 5 # radius[1] + -2 # offset[0] + 2 # offset[1] +) + # A enrichir SET(BasicFeatureExtraction_SRCS1 otbAlignImageToPath.cxx @@ -1216,6 +1234,8 @@ otbEnergyTextureFunctor.cxx otbEnergyTextureImageFunctionNew.cxx otbEnergyTextureImageFunction.cxx otbEntropyTextureFunctor.cxx +otbEntropyTextureImageFunctionNew.cxx +otbEntropyTextureImageFunction.cxx ) INCLUDE_DIRECTORIES("${OTBTesting_BINARY_DIR}") diff --git a/Testing/Code/FeatureExtraction/otbEnergyTextureImageFunction.cxx b/Testing/Code/FeatureExtraction/otbEnergyTextureImageFunction.cxx index 5df9617113bf22f759702b5f2e87b43a4ee3f6b7..f8b9e319147c18a393dc6af9d3a85f5311541fcc 100644 --- a/Testing/Code/FeatureExtraction/otbEnergyTextureImageFunction.cxx +++ b/Testing/Code/FeatureExtraction/otbEnergyTextureImageFunction.cxx @@ -18,13 +18,13 @@ #include "itkExceptionObject.h" #include "otbFunctionWithNeighborhoodToImageFilter.h" -#include "otbEnergyTextureImageFunction.h" +#include "otbEntropyTextureImageFunction.h" #include "otbImage.h" #include "otbImageFileReader.h" #include "otbStreamingImageFileWriter.h" -int otbEnergyTextureImageFunction(int argc, char * argv[]) +int otbEntropyTextureImageFunction(int argc, char * argv[]) { const char * inputFileName = argv[1]; const char * outputFileName = argv[2]; @@ -38,7 +38,7 @@ int otbEnergyTextureImageFunction(int argc, char * argv[]) typedef otb::ImageFileReader<ImageType> ReaderType; typedef otb::StreamingImageFileWriter<ImageType> WriterType; - typedef otb::EnergyTextureImageFunction<ImageType> FunctionType; + typedef otb::EntropyTextureImageFunction<ImageType> FunctionType; typedef otb::FunctionWithNeighborhoodToImageFilter<ImageType, ImageType, FunctionType> FilterType; FunctionType::Pointer energyFunction = FunctionType::New(); diff --git a/Testing/Code/FeatureExtraction/otbEnergyTextureImageFunctionNew.cxx b/Testing/Code/FeatureExtraction/otbEnergyTextureImageFunctionNew.cxx index 1ae8ad7d7194519a677797ad672da15a621ae0f6..ba481f761117554ed01eead976d820e54d911eb2 100644 --- a/Testing/Code/FeatureExtraction/otbEnergyTextureImageFunctionNew.cxx +++ b/Testing/Code/FeatureExtraction/otbEnergyTextureImageFunctionNew.cxx @@ -17,51 +17,23 @@ =========================================================================*/ #include "itkExceptionObject.h" -//#include "otbFunctionWithNeighborhoodToImageFilter.h" -#include "otbEnergyTextureImageFunction.h" + +#include "otbEntropyTextureImageFunction.h" #include "otbImage.h" -//#include "otbImageFileReader.h" -//#include "otbImageFileWriter.h" -int otbEnergyTextureImageFunctionNew(int argc, char * argv[]) +int otbEntropyTextureImageFunctionNew(int argc, char * argv[]) { //const char * inputFileName = argv[1]; //const char * outputFileName = argv[2]; - typedef double InputPixelType; + typedef double InputPixelType; const int Dimension = 2; typedef otb::Image<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::EnergyTextureFunctor<IterType, IterType, PixelType> FunctorType; - typedef otb::EnergyTextureImageFunction<ImageType> EnergyTextureImageFunctionType; - - EnergyTextureImageFunctionType::Pointer energy = EnergyTextureImageFunctionType::New(); - - // 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]); + typedef otb::EntropyTextureImageFunction<ImageType> EntropyTextureImageFunctionType; - object->SetOffset(offset); - writer->SetInput(object->GetOutput()); + EntropyTextureImageFunctionType::Pointer entropy = EntropyTextureImageFunctionType::New(); - writer->Update(); - */ return EXIT_SUCCESS; } diff --git a/Testing/Code/FeatureExtraction/otbEntropyTextureImageFunction.cxx b/Testing/Code/FeatureExtraction/otbEntropyTextureImageFunction.cxx new file mode 100644 index 0000000000000000000000000000000000000000..5df9617113bf22f759702b5f2e87b43a4ee3f6b7 --- /dev/null +++ b/Testing/Code/FeatureExtraction/otbEntropyTextureImageFunction.cxx @@ -0,0 +1,72 @@ +/*========================================================================= + + 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 "otbFunctionWithNeighborhoodToImageFilter.h" +#include "otbEnergyTextureImageFunction.h" +#include "otbImage.h" +#include "otbImageFileReader.h" +#include "otbStreamingImageFileWriter.h" + + +int otbEnergyTextureImageFunction(int argc, char * argv[]) +{ + const char * inputFileName = argv[1]; + const char * outputFileName = argv[2]; + + typedef double InputPixelType; + const int Dimension = 2; + typedef otb::Image<InputPixelType,Dimension> ImageType; + typedef ImageType::SizeType SizeType; + typedef ImageType::OffsetType OffsetType; + + typedef otb::ImageFileReader<ImageType> ReaderType; + typedef otb::StreamingImageFileWriter<ImageType> WriterType; + + typedef otb::EnergyTextureImageFunction<ImageType> FunctionType; + typedef otb::FunctionWithNeighborhoodToImageFilter<ImageType, ImageType, FunctionType> FilterType; + + FunctionType::Pointer energyFunction = FunctionType::New(); + FilterType::Pointer filter = FilterType::New(); + + + // Instantiating object + ReaderType::Pointer reader = ReaderType::New(); + WriterType::Pointer writer = WriterType::New(); + reader->SetFileName(inputFileName); + writer->SetFileName(outputFileName); + + filter->SetInput(reader->GetOutput()); + + SizeType radius; + radius[0] = atoi(argv[3]); + radius[1] = atoi(argv[4]); + energyFunction->SetRadius(radius); + OffsetType offset; + offset[0] = atoi(argv[5]); + offset[1] = atoi(argv[6]); + energyFunction->SetOffset(offset); + + filter->SetFunction(energyFunction); + writer->SetInput(filter->GetOutput()); + writer->SetNumberOfStreamDivisions(1); + writer->Update(); + + + return EXIT_SUCCESS; +} diff --git a/Testing/Code/FeatureExtraction/otbEntropyTextureImageFunctionNew.cxx b/Testing/Code/FeatureExtraction/otbEntropyTextureImageFunctionNew.cxx new file mode 100644 index 0000000000000000000000000000000000000000..7397bf4bacd3765e90f199fb42b547850ee15fb3 --- /dev/null +++ b/Testing/Code/FeatureExtraction/otbEntropyTextureImageFunctionNew.cxx @@ -0,0 +1,44 @@ +/*========================================================================= + + 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 "otbFunctionWithNeighborhoodToImageFilter.h" +#include "otbEnergyTextureImageFunction.h" +#include "otbImage.h" +//#include "otbImageFileReader.h" +//#include "otbImageFileWriter.h" + + +int otbEnergyTextureImageFunctionNew(int argc, char * argv[]) +{ + //const char * inputFileName = argv[1]; + //const char * outputFileName = argv[2]; + + typedef double InputPixelType; + const int Dimension = 2; + typedef otb::Image<InputPixelType,Dimension> ImageType; + typedef ImageType::PixelType PixelType; + typedef ImageType::OffsetType OffsetType; + + typedef otb::EnergyTextureImageFunction<ImageType> EnergyTextureImageFunctionType; + + EnergyTextureImageFunctionType::Pointer energy = EnergyTextureImageFunctionType::New(); + + + return EXIT_SUCCESS; +} diff --git a/Testing/Code/FeatureExtraction/otbFeatureExtractionTests11.cxx b/Testing/Code/FeatureExtraction/otbFeatureExtractionTests11.cxx index 33ff10183eef4366287d8cf2a6f61aa2567aa9e0..80e54914627ffe770cb33dd70b0bd22b0e629932 100644 --- a/Testing/Code/FeatureExtraction/otbFeatureExtractionTests11.cxx +++ b/Testing/Code/FeatureExtraction/otbFeatureExtractionTests11.cxx @@ -27,14 +27,16 @@ void RegisterTests() { - REGISTER_TEST(otbCloudEstimatorFilterNew); - REGISTER_TEST(otbCloudEstimatorDefaultFilter); - REGISTER_TEST(otbCloudEstimatorFilter); - REGISTER_TEST(otbCloudDetectionFilterNew); - REGISTER_TEST(otbCloudDetectionFilter); - REGISTER_TEST(otbSimplifyManyPathListFilter); - REGISTER_TEST(otbEnergyTextureFunctor); - REGISTER_TEST(otbEnergyTextureImageFunctionNew); - REGISTER_TEST(otbEnergyTextureImageFunction); - REGISTER_TEST(otbEntropyTextureFunctor); +REGISTER_TEST(otbCloudEstimatorFilterNew); +REGISTER_TEST(otbCloudEstimatorDefaultFilter); +REGISTER_TEST(otbCloudEstimatorFilter); +REGISTER_TEST(otbCloudDetectionFilterNew); +REGISTER_TEST(otbCloudDetectionFilter); +REGISTER_TEST(otbSimplifyManyPathListFilter); +REGISTER_TEST(otbEnergyTextureFunctor); +REGISTER_TEST(otbEnergyTextureImageFunctionNew); +REGISTER_TEST(otbEnergyTextureImageFunction); +REGISTER_TEST(otbEntropyTextureFunctor); +REGISTER_TEST(otbEntropyTextureImageFunctionNew); +REGISTER_TEST(otbEntropyTextureImageFunction); }