diff --git a/Code/FeatureExtraction/otbHaralickTexturesImageFunction.h b/Code/FeatureExtraction/otbHaralickTexturesImageFunction.h
new file mode 100644
index 0000000000000000000000000000000000000000..e36b5847c1abcc2d6804c80e2e97086ef6391a33
--- /dev/null
+++ b/Code/FeatureExtraction/otbHaralickTexturesImageFunction.h
@@ -0,0 +1,215 @@
+/*=========================================================================
+
+ 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 __otbHaralickTexturesImageFunction_h
+#define __otbHaralickTexturesImageFunction_h
+
+#include "itkImageFunction.h"
+#include "itkFixedArray.h"
+#include "otbMaskedScalarImageToGreyLevelCoocurenceMatrixGenerator.h"
+#include "itkGreyLevelCooccurrenceMatrixTextureCoefficientsCalculator.h"
+
+namespace otb
+{
+/**
+ * \class HaralickTexturesImageFunction
+ * \brief Compute the 8 Haralick texture indices on the neighborhood of the point
+ *
+ * This class computes the 8 Haralick texture indices on the neighborhood of the point
+ * (where \f$ g(i, j) \f$ is the element in
+ * cell i, j of a normalized GLCM):
+ *
+ * "Energy" \f$ = f_1 = \sum_{i,j}g(i, j)^2 \f$
+ *
+ * "Entropy" \f$ = f_2 = -\sum_{i,j}g(i, j) \log_2 g(i, j)\f$, or 0 if \f$g(i, j) = 0\f$
+ *
+ * "Correlation" \f$ = f_3 = \sum_{i,j}\frac{(i - \mu)(j - \mu)g(i, j)}{\sigma^2} \f$
+ *
+ * "Difference Moment" \f$= f_4 = \sum_{i,j}\frac{1}{1 + (i - j)^2}g(i, j) \f$
+ *
+ * "Inertia" \f$ = f_5 = \sum_{i,j}(i - j)^2g(i, j) \f$ (sometimes called "contrast")
+ *
+ * "Cluster Shade" \f$ = f_6 = \sum_{i,j}((i - \mu) + (j - \mu))^3 g(i, j) \f$
+ *
+ * "Cluster Prominence" \f$ = f_7 = \sum_{i,j}((i - \mu) + (j - \mu))^4 g(i, j) \f$
+ *
+ * "Haralick's Correlation" \f$ = f_8 = \frac{\sum_{i,j}(i, j) g(i, j) -\mu_t^2}{\sigma_t^2} \f$
+ * where \f$\mu_t\f$ and \f$\sigma_t\f$ are the mean and standard deviation of the row
+ * (or column, due to symmetry) sums.
+ *
+ * Above, \f$ \mu = \f$ (weighted pixel average) \f$ = \sum_{i,j}i \cdot g(i, j) =
+ * \sum_{i,j}j \cdot g(i, j) \f$ (due to matrix summetry), and
+ *
+ * \f$ \sigma = \f$ (weighted pixel variance) \f$ = \sum_{i,j}(i - \mu)^2 \cdot g(i, j) =
+ * \sum_{i,j}(j - \mu)^2 \cdot g(i, j) \f$ (due to matrix summetry)
+ * This class is templated over the input image type and the
+ * coordinate representation type (e.g. float or double).
+ *
+ * \sa otb::MaskedScalarImageToGreyLevelCooccurrenceMatrixGenerator
+ * \sa itk::GreyLevelCooccurrenceMatrixTextureCoefficientsCalculator
+ *
+ * \ingroup ImageFunctions
+ */
+
+template <class TInputImage, class TCoordRep = double >
+class ITK_EXPORT HaralickTexturesImageFunction :
+ public itk::ImageFunction< TInputImage,
+ itk::FixedArray<
+ ITK_TYPENAME itk::NumericTraits<typename TInputImage::PixelType>::RealType,
+ 8 >,
+ TCoordRep >
+{
+public:
+ /** Standard class typedefs. */
+ typedef HaralickTexturesImageFunction Self;
+ typedef itk::ImageFunction< TInputImage,
+ itk::FixedArray<
+ ITK_TYPENAME itk::NumericTraits<
+ typename TInputImage::PixelType>::RealType,
+ 8 >,
+ TCoordRep > Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Run-time type information (and related methods). */
+ itkTypeMacro(HaralickTexturesImageFunction, ImageFunction);
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ /** InputImageType typedef support. */
+ typedef TInputImage InputImageType;
+ typedef typename Superclass::IndexType IndexType;
+ typedef typename Superclass::ContinuousIndexType ContinuousIndexType;
+ typedef typename Superclass::PointType PointType;
+ typedef typename InputImageType::Pointer InputImagePointerType;
+ typedef typename InputImageType::PixelType InputPixelType;
+ typedef typename InputImageType::RegionType InputRegionType;
+ typedef typename InputRegionType::SizeType SizeType;
+
+ // Textures tools
+ /** Co-occurence matrix and textures calculator */
+ typedef otb::MaskedScalarImageToGreyLevelCooccurrenceMatrixGenerator
+ <InputImageType> CoocurrenceMatrixGeneratorType;
+ typedef typename CoocurrenceMatrixGeneratorType
+ ::Pointer CoocurrenceMatrixGeneratorPointerType;
+ typedef typename CoocurrenceMatrixGeneratorType
+ ::OffsetType OffsetType;
+ typedef typename CoocurrenceMatrixGeneratorType
+ ::HistogramType HistogramType;
+ typedef itk::Statistics::GreyLevelCooccurrenceMatrixTextureCoefficientsCalculator
+ <HistogramType> TextureCoefficientsCalculatorType;
+ typedef typename TextureCoefficientsCalculatorType
+ ::Pointer TextureCoefficientsCalculatorPointerType;
+
+ // Output typedef support
+ typedef typename Superclass::OutputType OutputType;
+ typedef typename OutputType::ValueType ScalarRealType;
+
+ // Coord rep
+ typedef TCoordRep CoordRepType;
+
+ /** Dimension of the underlying image. */
+ itkStaticConstMacro(ImageDimension, unsigned int,
+ InputImageType::ImageDimension);
+
+ /** Evalulate the function at specified index */
+ virtual OutputType EvaluateAtIndex(const IndexType& index) const;
+
+ /** Evaluate the function at non-integer positions */
+ virtual OutputType Evaluate(const PointType& point) const
+ {
+ IndexType index;
+ this->ConvertPointToNearestIndex(point, index);
+ return this->EvaluateAtIndex(index);
+ }
+ virtual OutputType EvaluateAtContinuousIndex(
+ const ContinuousIndexType& cindex) const
+ {
+ IndexType index;
+ this->ConvertContinuousIndexToNearestIndex(cindex, index);
+ return this->EvaluateAtIndex(index);
+ }
+
+ /**
+ * Set the radius of the neighborhood over which the
+ * statistics are evaluated
+ */
+ itkSetMacro( NeighborhoodRadius, unsigned int );
+ /**
+ * Get the radius of the neighborhood over which the
+ * statistics are evaluated
+ */
+ itkGetConstReferenceMacro( NeighborhoodRadius, unsigned int );
+
+ /** Set the offset for co-occurence computation */
+ itkSetMacro(Offset, OffsetType);
+
+ /** Get the offset for co-occurence computation */
+ itkGetMacro(Offset, OffsetType);
+
+ /** Set the number of bin per axis for histogram generation */
+ itkSetMacro(NumberOfBinsPerAxis, unsigned int);
+
+ /** Get the number of bin per axis for histogram generation */
+ itkGetMacro(NumberOfBinsPerAxis, unsigned int);
+
+ /** Set the input image minimum */
+ itkSetMacro(InputImageMinimum, InputPixelType);
+
+ /** Get the input image minimum */
+ itkGetMacro(InputImageMinimum, InputPixelType);
+
+ /** Set the input image maximum */
+ itkSetMacro(InputImageMaximum, InputPixelType);
+
+ /** Get the input image maximum */
+ itkGetMacro(InputImageMaximum, InputPixelType);
+
+protected:
+ HaralickTexturesImageFunction();
+ virtual ~HaralickTexturesImageFunction() {}
+ void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+private:
+ HaralickTexturesImageFunction(const Self &); //purposely not implemented
+ void operator =(const Self&); //purposely not implemented
+
+ /** Radius of the neighborhood over which to compute the textures */
+ unsigned int m_NeighborhoodRadius;
+
+ /** Offset for co-occurences generation */
+ OffsetType m_Offset;
+
+ /** Number of bins per axis for histogram generation */
+ unsigned int m_NumberOfBinsPerAxis;
+
+ /** Input image minimum */
+ InputPixelType m_InputImageMinimum;
+
+ /** Input image maximum */
+ InputPixelType m_InputImageMaximum;
+};
+
+} // namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbHaralickTexturesImageFunction.txx"
+#endif
+
+#endif
+
diff --git a/Code/FeatureExtraction/otbHaralickTexturesImageFunction.txx b/Code/FeatureExtraction/otbHaralickTexturesImageFunction.txx
new file mode 100644
index 0000000000000000000000000000000000000000..54d5a05f3c3c3fa67d301da0febc2d85787ae792
--- /dev/null
+++ b/Code/FeatureExtraction/otbHaralickTexturesImageFunction.txx
@@ -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 __otbHaralickTexturesImageFunction_txx
+#define __otbHaralickTexturesImageFunction_txx
+
+#include "otbHaralickTexturesImageFunction.h"
+#include "itkConstNeighborhoodIterator.h"
+#include "itkNumericTraits.h"
+#include "itkMacro.h"
+#include <complex>
+
+namespace otb
+{
+
+/**
+ * Constructor
+ */
+template <class TInputImage, class TCoordRep>
+HaralickTexturesImageFunction<TInputImage, TCoordRep>
+::HaralickTexturesImageFunction() :
+ m_NeighborhoodRadius(10),
+ m_Offset(),
+ m_NumberOfBinsPerAxis(8),
+ m_InputImageMinimum(0),
+ m_InputImageMaximum(256)
+{}
+
+template <class TInputImage, class TCoordRep>
+void
+HaralickTexturesImageFunction<TInputImage, TCoordRep>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+ this->Superclass::PrintSelf(os, indent);
+ os << indent << " Neighborhood radius value : " << m_NeighborhoodRadius << std::endl;
+}
+
+template <class TInputImage, class TCoordRep>
+typename HaralickTexturesImageFunction<TInputImage,TCoordRep>::OutputType
+HaralickTexturesImageFunction<TInputImage,TCoordRep>
+::EvaluateAtIndex(const IndexType& index) const
+{
+ // Build textures vector
+ OutputType textures;
+
+ // Initialize textures
+ textures.Fill( itk::NumericTraits< ScalarRealType >::Zero );
+
+ // Check for input image
+ if( !this->GetInputImage() )
+ {
+ return textures;
+ }
+
+ // Check for out of buffer
+ if ( !this->IsInsideBuffer( index ) )
+ {
+ return textures;
+ }
+
+ // Build the co-occurence matrix generator
+ CoocurrenceMatrixGeneratorPointerType coOccurenceMatrixGenerator = CoocurrenceMatrixGeneratorType::New();
+ coOccurenceMatrixGenerator->SetInput(this->GetInputImage());
+ coOccurenceMatrixGenerator->SetOffset(m_Offset);
+ coOccurenceMatrixGenerator->SetNumberOfBinsPerAxis(m_NumberOfBinsPerAxis);
+ coOccurenceMatrixGenerator->SetPixelValueMinMax(m_InputImageMinimum, m_InputImageMaximum);
+
+ // Compute the region on which co-occurence will be estimated
+ typename InputRegionType::IndexType inputIndex = index;
+ for(unsigned int dim = 0; dim<InputImageType::ImageDimension;++dim)
+ {
+ inputIndex[dim]-= m_NeighborhoodRadius;
+ }
+ typename InputRegionType::SizeType inputSize;
+ inputSize.Fill(2*m_NeighborhoodRadius+1);
+
+ // Build the input region
+ InputRegionType inputRegion;
+ inputRegion.SetIndex(inputIndex);
+ inputRegion.SetSize(inputSize);
+
+ // Compute the co-occurence matrix
+ coOccurenceMatrixGenerator->SetRegion(inputRegion);
+ coOccurenceMatrixGenerator->SetNormalize(true);
+ coOccurenceMatrixGenerator->Compute();
+
+ // Build the texture calculator
+ TextureCoefficientsCalculatorPointerType texturesCalculator = TextureCoefficientsCalculatorType::New();
+
+ // Compute textures indices
+ texturesCalculator->SetHistogram(coOccurenceMatrixGenerator->GetOutput());
+ texturesCalculator->Compute();
+
+ // Fill the output vector
+ textures[0]=texturesCalculator->GetEnergy();
+ textures[1]=texturesCalculator->GetEntropy();
+ textures[2]=texturesCalculator->GetCorrelation();
+ textures[3]=texturesCalculator->GetInverseDifferenceMoment();
+ textures[4]=texturesCalculator->GetInertia();
+ textures[5]=texturesCalculator->GetClusterShade();
+ textures[6]=texturesCalculator->GetClusterProminence();
+ textures[7]=texturesCalculator->GetHaralickCorrelation();
+
+ // Return result
+ return textures;
+}
+
+} // namespace otb
+
+#endif
diff --git a/Testing/Code/FeatureExtraction/CMakeLists.txt b/Testing/Code/FeatureExtraction/CMakeLists.txt
index 12b13c2b26f4b4e0cb5355e71882cd2c5119a5bb..eb837764ebed4b23525aaced079c244b3032ae54 100644
--- a/Testing/Code/FeatureExtraction/CMakeLists.txt
+++ b/Testing/Code/FeatureExtraction/CMakeLists.txt
@@ -1528,6 +1528,23 @@ ADD_TEST(feTvMetaImageFunction ${FEATUREEXTRACTION_TESTS16}
451846.014047961 5412466.57452216
)
+# ------- otb::HaralickTexturesImageFunction -------------
+ADD_TEST(feTuHaralickTexturesImageFunctionNew ${FEATUREEXTRACTION_TESTS16}
+ otbHaralickTexturesImageFunctionNew
+)
+
+ADD_TEST(feTvHaralickTexturesImageFunction ${FEATUREEXTRACTION_TESTS16}
+--compare-ascii ${EPSILON_8}
+ ${BASELINE_FILES}/feTvHaralickTexturesImageFunction.txt
+ ${TEMP}/feTvHaralickTexturesImageFunction.txt
+ otbHaralickTexturesImageFunction
+ ${INPUTDATA}/ROI_IKO_PAN_LesHalles.tif
+ ${TEMP}/feTvHaralickTexturesImageFunction.txt
+ 451846.014047961 5412466.57452216
+)
+
+
+
# A enrichir
SET(BasicFeatureExtraction_SRCS1
otbFeatureExtractionTests1.cxx
@@ -1724,6 +1741,7 @@ otbLocalHistogramImageFunctionNew.cxx
otbLocalHistogramImageFunctionTest.cxx
otbImageFunctionAdaptor.cxx
otbMetaImageFunction.cxx
+otbHaralickTexturesImageFunction.cxx
)
diff --git a/Testing/Code/FeatureExtraction/otbFeatureExtractionTests16.cxx b/Testing/Code/FeatureExtraction/otbFeatureExtractionTests16.cxx
index 4ce2d662b6c6c19970ed23327ab6483f78fe0892..2dc89f67426caeaf521b69ec651e11a9a504f657 100644
--- a/Testing/Code/FeatureExtraction/otbFeatureExtractionTests16.cxx
+++ b/Testing/Code/FeatureExtraction/otbFeatureExtractionTests16.cxx
@@ -36,4 +36,6 @@ void RegisterTests()
REGISTER_TEST(otbImageFunctionAdaptor);
REGISTER_TEST(otbMetaImageFunctionNew);
REGISTER_TEST(otbMetaImageFunction);
+ REGISTER_TEST(otbHaralickTexturesImageFunctionNew);
+ REGISTER_TEST(otbHaralickTexturesImageFunction);
}
diff --git a/Testing/Code/FeatureExtraction/otbHaralickTexturesImageFunction.cxx b/Testing/Code/FeatureExtraction/otbHaralickTexturesImageFunction.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b2cbc51e7f3f24fd73c0c4f35e108ac1e7e2c58c
--- /dev/null
+++ b/Testing/Code/FeatureExtraction/otbHaralickTexturesImageFunction.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.
+
+=========================================================================*/
+#if defined(_MSC_VER)
+#pragma warning ( disable : 4786 )
+#endif
+
+#include "otbHaralickTexturesImageFunction.h"
+#include "otbImage.h"
+#include "otbImageFileReader.h"
+
+typedef unsigned short InputPixelType;
+const unsigned int Dimension = 2;
+
+typedef otb::Image<InputPixelType, Dimension> InputImageType;
+typedef otb::ImageFileReader<InputImageType> ReaderType;
+
+typedef otb::HaralickTexturesImageFunction<
+ InputImageType,double> HaralickTexturesImageFunctionType;
+typedef HaralickTexturesImageFunctionType::PointType PointType;
+typedef HaralickTexturesImageFunctionType::OutputType OutputType;
+
+int otbHaralickTexturesImageFunctionNew(int argc, char * argv[])
+{
+ HaralickTexturesImageFunctionType::Pointer function = HaralickTexturesImageFunctionType::New();
+
+ return EXIT_SUCCESS;
+}
+
+int otbHaralickTexturesImageFunction(int argc, char * argv[])
+{
+ // Read the input image
+ ReaderType::Pointer reader = ReaderType::New();
+ reader->SetFileName(argv[1]);
+ reader->Update();
+
+ std::ofstream outputStream(argv[2]);
+
+ HaralickTexturesImageFunctionType::Pointer haralick = HaralickTexturesImageFunctionType::New();
+ haralick->SetInputImage(reader->GetOutput());
+ haralick->SetNeighborhoodRadius(10);
+
+ HaralickTexturesImageFunctionType::OffsetType offset;
+ offset.Fill(1);
+ haralick->SetOffset(offset);
+
+ PointType p;
+ p[0] = atof(argv[3]);
+ p[1] = atof(argv[4]);
+
+ OutputType output = haralick->Evaluate(p);
+
+ outputStream<<"Evaluate("<<p<<") = "<<output<<std::endl;
+
+ outputStream.close();
+
+ return EXIT_SUCCESS;
+}