diff --git a/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilter.h b/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilter.h
new file mode 100644
index 0000000000000000000000000000000000000000..c0a5750a9ff705f4cac23c4d757bb97159b594de
--- /dev/null
+++ b/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilter.h
@@ -0,0 +1,228 @@
+/*=========================================================================
+
+Program: ORFEO Toolbox
+Language: C++
+Date: $Date$
+Version: $Revision$
+
+
+Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+See OTBCopyright.txt for details.
+
+Copyright (c) CS Systemes d'information. All rights reserved.
+See CSCopyright.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 _otbImageToSURFKeyPointSetFilter_h
+#define _otbImageToSURFKeyPointSetFilter_h
+
+#include "itkProcessObject.h"
+#include "itkConstNeighborhoodIterator.h"
+#include "otbImageToHessianDeterminantImageFilter.h"
+#include "itkPointSet.h"
+#include <itkRescaleIntensityImageFilter.h>
+#include <otbImageList.h>
+#include "itkResampleImageFilter.h"
+#include "otbImageToPointSetFilter.h"
+
+#include "otbMath.h"
+
+/** \class ImageToSURFKeyPointSetFilter
+ * \brief This class extracts key points from an image through a pyramidal gaussian based decomposition
+ *
+ * This class implements the SURF Key point detector proposed by Tuytelaars and Vangool from the university
+ * of Leuven, 2005
+ *
+ * \li Gaussian Second order derivative Hessian images are computed in each
+ * level and each octave for the input image.
+ * \li For each octave, an extremum search is launched on each 3 adjacent scales.
+ * \li The Key points detected are the ones extremum in the current , previous and
+ * next scale of reserach. 3 scales are needed for the computation (NumberOfScales >=3).
+ * \li Orientation and descriptors are computed for each key point detected.
+ *
+ * Selected Key Points are stored in an itk::PointSet structure.
+ * Points contains the coordinate of the detected point.
+ * DataPoints contain the values of the 64 element descriptor for each key point
+ * detected through the pyramidal analysis.
+ *
+ * Orientation is expressed in degree in the range of [0,360]
+ *
+ * \sa otb::ImageToDeterminantHessianImage
+ */
+
+namespace otb
+{
+ template <class TInputImage , class TOutputPointSet>
+ class ITK_EXPORT ImageToSURFKeyPointSetFilter
+ : public ImageToPointSetFilter<TInputImage, TOutputPointSet>
+ {
+
+ public:
+
+ /** Standard class typedefs. */
+ typedef ImageToSURFKeyPointSetFilter Self;
+ typedef ImageToPointSetFilter<TInputImage,TOutputPointSet> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ /** Run-time type information (and related methods). */
+ itkTypeMacro(ImageToSURFKeyPointSetFilter,ImageToPointSetFilter);
+
+
+ /** Template parameters typedefs*/
+ typedef TInputImage InputImageType;
+ typedef typename InputImageType::Pointer InputImagePointerType;
+ typedef typename InputImageType::IndexType PixelIndex ;
+ typedef typename InputImageType::IndexType IndexType;
+ typedef typename InputImageType::PixelType PixelValue;
+ typedef typename InputImageType::SpacingType SpacingType;
+ typedef typename InputImageType::SizeType SizeType;
+ typedef typename InputImageType::PointType PointImageType;
+
+ typedef TOutputPointSet OutputPointSetType;
+ typedef typename TOutputPointSet::Pointer OutputPointSetPointerType;
+ typedef typename TOutputPointSet::PixelType OutputPixelType;
+ typedef typename TOutputPointSet::PointType OutputPointType;
+ typedef typename TOutputPointSet::PointIdentifier OutputPointIdentifierType;
+
+ /** Set/Get the number of Octaves */
+ itkSetMacro(OctavesNumber,int);
+ itkGetMacro(OctavesNumber,int);
+
+ /** Set/Get the number of scales*/
+ itkSetMacro(ScalesNumber,int);
+ itkGetMacro(ScalesNumber,int);
+
+ /** Get the number of KeyPoints detected*/
+ itkGetMacro(NumberOfPoints,int);
+
+ /** Internal filters typedefs */
+ typedef itk::ConstNeighborhoodIterator<InputImageType> NeighborhoodIteratorType;
+ typedef typename NeighborhoodIteratorType::NeighborhoodType NeighborhoodType;
+
+ typedef otb::ImageToHessianDeterminantImageFilter<InputImageType ,InputImageType > ImageToDetHessianImageType;
+ typedef typename ImageToDetHessianImageType::Pointer DetHessianPointerFilter;
+
+ /** Filter for resampling images for the multi-scale analysis */
+ typedef itk::ResampleImageFilter<InputImageType,InputImageType> ResampleFilterType;
+ typedef typename ResampleFilterType::Pointer ResampleFilterPointerType;
+
+ /** ImageList to store the Hessian determinant image at each scale (sigma width)*/
+ typedef otb::ImageList< InputImageType > ImageListType;
+ typedef typename ImageListType::Pointer ImageListTypePointer;
+
+ /** */
+ typedef std::vector<double> VectorType;
+
+ protected:
+
+ /**
+ * Constructor.
+ */
+ ImageToSURFKeyPointSetFilter();
+ /**
+ * Destructor.
+ */
+ virtual ~ImageToSURFKeyPointSetFilter();
+ /**
+ * Standard PrintSelf method.
+ */
+ virtual void PrintSelf(std::ostream& os, itk::Indent indent) const;
+ /**
+ * Main computation method.
+ */
+ virtual void GenerateData();
+
+
+
+ /** Check local extremum for 8 neighbors (current)
+ *
+ * \param currentScale
+ *
+ * \return true if the Central pixel is extremum
+ */
+ virtual bool IsLocalExtremum(const NeighborhoodType& neigh);
+
+
+
+ /** Check local extremum for 8 neighbors (Previous or Scale)
+ *
+ * \param neigh
+ * \paramCenterValue
+ *
+ * \return true if the Central pixel is extremum
+ */
+ virtual bool IsLocalExtremumAround(const NeighborhoodType& neigh,
+ double CenterValue);
+
+
+ /** AssignOrientation
+ *
+ * \param currentScale neighborhood
+ * \param scale affected to the keypoint
+ *
+ * \return key point orientation
+ */
+ virtual double AssignOrientation(const NeighborhoodType& neigh,
+ double S);
+
+ /** ComputeDescriptor
+ *
+ * \param currentScale Neighboorhood
+ * \param orientation assigned to the key point
+ * \param scale
+ *
+ * \return hsitogram descriptor
+ */
+ virtual VectorType ComputeDescriptor(const NeighborhoodType& neigh,
+ double O,
+ double S);
+ /**
+ * Compute min a b c
+ */
+ virtual int GetMin( int a , int b , int c);
+
+ private:
+
+ ImageToSURFKeyPointSetFilter(const Self&); //purposely not implemented
+ void operator=(const Self&); //purposely not implemented
+
+ /** Number of octaves */
+ int m_OctavesNumber;
+
+ /** Number of scale for each octave */
+ int m_ScalesNumber;
+
+ /** Number of key points detected */
+ int m_NumberOfPoints;
+
+ /** Those images */
+ InputImagePointerType m_determinantImage;
+ InputImagePointerType m_ImageCurrent;
+ InputImagePointerType m_ImageMovedPrev;
+ InputImagePointerType m_ImageMovedNext;
+
+ /** ImageToDeterminantHessianFilter filter */
+ DetHessianPointerFilter m_DetHessianFilter;
+
+ /*Resample Filter*/
+ ResampleFilterPointerType m_ResampleFilter;
+
+ /** ImageList*/
+ ImageListTypePointer m_ImageList;
+ };
+}
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbImageToSURFKeyPointSetFilter.txx"
+#endif
+
+#endif
+
+
diff --git a/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilter.txx b/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilter.txx
new file mode 100644
index 0000000000000000000000000000000000000000..4d58522d7a9f052411fd4ec1461b677ac247af03
--- /dev/null
+++ b/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilter.txx
@@ -0,0 +1,527 @@
+/*=========================================================================
+
+Program: ORFEO Toolbox
+Language: C++
+Date: $Date$
+Version: $Revision$
+
+
+Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+See OTBCopyright.txt for details.
+
+Copyright (c) CS Systemes d'information. All rights reserved.
+See CSCopyright.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 "otbImageToSURFKeyPointSetFilter.h"
+#include "itkCenteredRigid2DTransform.h"
+
+namespace otb
+{
+ /**---------------------------------------------------------
+ * Constructor
+ ----------------------------------------------------------*/
+ template < class TInputImage, class TOutputPointSet >
+ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet >
+ ::ImageToSURFKeyPointSetFilter()
+ {
+ m_OctavesNumber = 1;
+ m_ScalesNumber = 3;
+ m_NumberOfPoints = 0;
+ m_DetHessianFilter = ImageToDetHessianImageType::New();
+ }
+
+
+ /*---------------------------------------------------------
+ * Destructor.c
+ ----------------------------------------------------------*/
+
+ template <class TInputImage, class TOutputPointSet>
+ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
+ ::~ImageToSURFKeyPointSetFilter()
+ {}
+
+ /*-------------------------------------------------------
+ * Generate Data
+ --------------------------------------------------------*/
+ template <class TInputImage, class TOutputPointSet>
+ void
+ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
+ ::GenerateData(void)
+ {
+
+
+ double k;
+ double sigma_in = 2.;
+ SizeType radius;
+
+ /*Output*/
+ OutputPointSetPointerType outputPointSet = this->GetOutput();
+
+ /** Computing the multiplicative factor for scales iteration
+ * scalar used for the wifth of the gaussian
+ */
+ if(m_ScalesNumber > 1 )
+ k = (double)std::pow(2.0,(double)(1/(double)(m_ScalesNumber-1))) ;
+ else
+ k = 3;
+
+
+ /* Computation loop over octaves*/
+ for (int i = 0; i < m_OctavesNumber ; i++ ){
+
+ sigma_in = 2.;
+ m_ImageList = ImageListType::New();
+
+ /*--------------------------------------------------------
+ Octave per octave
+ --------------------------------------------------------*/
+ if(i>0) {
+
+ m_ResampleFilter = ResampleFilterType::New();
+ m_ResampleFilter ->SetInput(this->GetInput(0));
+
+ SizeType size = this->GetInput(0)->GetLargestPossibleRegion().GetSize();
+ for (int k = 0; k < 2; ++k)
+ size[k] = (unsigned int) floor(size[k]/std::pow(2.0,i) );
+ m_ResampleFilter->SetSize( size );
+
+ SpacingType spacing = this->GetInput(0)->GetSpacing();
+ for (int k = 0; k < 2; ++k)
+ spacing[k] = (spacing[k] * std::pow(2.0,i));
+ m_ResampleFilter->SetOutputSpacing( spacing );
+
+ m_ResampleFilter->SetDefaultPixelValue( 0 );
+ m_ResampleFilter->Update();
+ m_determinantImage = m_ResampleFilter->GetOutput();
+
+ otbGenericMsgDebugMacro( <<"ImageToSURFKeyPointSetFilter:: Size of the image at the octave : " \
+ << i << " is " \
+ <<m_determinantImage->GetLargestPossibleRegion().GetSize() );
+
+ }
+
+ for (int j = 0 ; j < m_ScalesNumber; j++ )
+ {
+ /** Incrementation of the gaussian width
+ * the width of the gaussian have to be doubled for
+ * each iteration over octaves
+ */
+
+ if ((i != 0 && j !=0 ) || (i == 0 && (i+j !=0) ) || ( m_ScalesNumber == 1 && i!=0 ))
+ sigma_in *= k;
+
+ /**
+ * For each octave, we serach for the key points
+ */
+
+ /** Hessian Determinant Image */
+ m_DetHessianFilter = ImageToDetHessianImageType::New();
+
+ if ( i == 0 )m_DetHessianFilter->SetInput(this->GetInput(0));
+ else m_DetHessianFilter->SetInput(m_determinantImage );
+
+ m_DetHessianFilter->SetSigma(sigma_in);
+ m_DetHessianFilter->Update();
+ m_determinantImage = m_DetHessianFilter->GetOutput() ;
+
+ if(i+j==0)
+ otbGenericMsgDebugMacro( <<"ImageToSURFKeyPointSetFilter:: Size of the image at the octave : " \
+ << i << " is " \
+ <<m_determinantImage->GetLargestPossibleRegion().GetSize() );
+
+
+ /** For each octave, we fill the imageList for the extremum search*/
+ m_ImageList->PushBack(m_determinantImage);
+ }
+
+ /*----------------------------------------------------*/
+ /* extremum Search over octave's scales */
+ /*----------------------------------------------------*/
+
+ for (int jj = 1 ; jj < (int)(m_ImageList->Size() - 1 ) ; jj++)
+ {
+ m_ImageCurrent = m_ImageList->GetNthElement(jj);
+ m_ImageMovedPrev = m_ImageList->GetNthElement(jj-1);
+ m_ImageMovedNext = m_ImageList->GetNthElement(jj+1);
+
+
+ /** NeighboorhoodIterator parameters*/
+ radius.Fill(1);
+ NeighborhoodIteratorType it(radius, m_ImageCurrent,m_ImageCurrent->GetLargestPossibleRegion());
+ it.GoToBegin();
+
+
+ /* NeighboorhoodIterator Adjacents parameters*/
+ NeighborhoodIteratorType itNeighPrev(radius, m_ImageMovedPrev,m_ImageMovedPrev->GetLargestPossibleRegion());
+ itNeighPrev.GoToBegin();
+
+ NeighborhoodIteratorType itNeighNext(radius, m_ImageMovedNext,m_ImageMovedNext->GetLargestPossibleRegion());
+ itNeighNext.GoToBegin();
+
+ while(!it.IsAtEnd())
+ {
+
+ if(IsLocalExtremum(it.GetNeighborhood())
+ && IsLocalExtremumAround(itNeighPrev.GetNeighborhood(),m_ImageCurrent->GetPixel(it.GetIndex()))
+ && IsLocalExtremumAround(itNeighNext.GetNeighborhood(),m_ImageCurrent->GetPixel(it.GetIndex())) )
+ {
+ OutputPointType keyPoint;
+ itNeighPrev.SetLocation(it.GetIndex());
+ itNeighNext.SetLocation(it.GetIndex());
+
+ keyPoint[0] = it.GetIndex()[0];
+ keyPoint[1] = it.GetIndex()[1];
+
+ //keyPoint[2] = sigma_in/pow(k,(double)jj)*pow(2.,(double)i);
+ double sigmaDetected = sigma_in/pow(k,(double)jj)*pow(2.,(double)i);
+
+ radius.Fill(GetMin((int)(this->GetInput(0)->GetLargestPossibleRegion().GetSize()[0] - keyPoint[0]),
+ (int)(this->GetInput(0)->GetLargestPossibleRegion().GetSize()[1] - keyPoint[1]),
+ (int)(6*keyPoint[2]) ) ) ;
+
+
+
+ /*
+ Computing the orientation of the key point detected
+ */
+ NeighborhoodIteratorType itNeighOrientation(radius,this->GetInput(0) ,
+ this->GetInput(0)->GetLargestPossibleRegion());
+
+ itNeighOrientation.SetLocation(it.GetIndex());
+
+ /** TO DO*/
+ //keyPoint[3] = AssignOrientation( itNeighOrientation.GetNeighborhood() ,keyPoint[2] );
+ double orientationDetected = AssignOrientation( itNeighOrientation.GetNeighborhood() , sigmaDetected );
+
+ /*Filling the Point pointSet Part*/
+ outputPointSet->SetPoint(m_NumberOfPoints, keyPoint);
+
+
+ /*----------------------------------------*/
+ /* Descriptor Computation */
+ /*----------------------------------------*/
+
+ radius.Fill(GetMin((int)(this->GetInput(0)->GetLargestPossibleRegion().GetSize()[0] - keyPoint[0]),
+ (int)(this->GetInput(0)->GetLargestPossibleRegion().GetSize()[1] - keyPoint[1]),
+ (int)(10*sigmaDetected))); // TODO a changer sigmaDetected par Keypoint[2]
+
+ NeighborhoodIteratorType itNeighDescriptor(radius,this->GetInput(0),
+ this->GetInput(0)->GetLargestPossibleRegion());
+ itNeighDescriptor.SetLocation(it.GetIndex());
+ VectorType descriptor;
+ //descriptor = ComputeDescriptor(itNeighDescriptor.GetNeighborhood(),keyPoint[3],keyPoint[2]);
+ descriptor = ComputeDescriptor(itNeighDescriptor.GetNeighborhood(),orientationDetected,sigmaDetected);
+
+
+
+ /*Updating the pointset data with values of the descriptor*/
+ OutputPixelType data;
+ data.SetSize(64);
+
+ unsigned int IndDescriptor = 0;
+
+ typename std::vector<double>::const_iterator itDescriptor =
+ descriptor.begin();
+
+ while(itDescriptor != descriptor.end())
+ {
+ data.SetElement(IndDescriptor, *itDescriptor);
+ IndDescriptor++;
+ itDescriptor++;
+ }
+ outputPointSet->SetPointData(m_NumberOfPoints, data);
+
+ m_NumberOfPoints++;
+ }
+ ++it;
+ ++itNeighPrev;
+ ++itNeighNext;
+
+ }/*End while for extremum search*/
+
+ } /*End Iteration over scales */
+
+ m_ImageList->Clear();
+
+ } /* End Key Points search*/
+
+ otbGenericMsgDebugMacro( <<"ImageToSURFKeyPointSetFilter:: Total Number of Key points "\
+ << m_NumberOfPoints );
+
+ }/** End of GenerateData()*/
+
+ template <class TInputImage, class TOutputPointSet>
+ int
+ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
+ ::GetMin( int a , int b , int c)
+ {
+ return std::min(a,std::min(b,c));
+ }
+
+ /*-----------------------------------------------------------
+ * Find Local Extremum
+ -----------------------------------------------------------*/
+ template <class TInputImage, class TOutputPointSet>
+ bool
+ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
+ ::IsLocalExtremum(const NeighborhoodType& neigh)
+ {
+ int centerIndex = neigh.GetCenterNeighborhoodIndex(), i = 0;
+ double centerValue = neigh[centerIndex];
+ bool max = false, min = false;
+ int flag_min = 0, flag_max = 0;
+
+ while (i!=(int)neigh.Size()){
+ if(i != centerIndex ){
+ if( centerValue> neigh[i] & flag_max == 0) max = true;
+ else { max = false; flag_max = 1; }
+
+ if(centerValue < neigh[i] & flag_min == 0 & centerValue <0) min = true;
+ else { min = false; flag_min = 1; }
+ }
+ ++i;
+ }
+
+ return max || min ;
+ }
+
+ /*-----------------------------------------------------------
+ *Find Local Extremum Around
+ -----------------------------------------------------------*/
+ template <class TInputImage, class TOutputPointSet>
+ bool
+ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
+ ::IsLocalExtremumAround(const NeighborhoodType& neigh , double CenterValue)
+ {
+
+ int i = 0;
+ bool max = false, min = false;
+ int flag_min = 0, flag_max = 0;
+
+ while (i!=(int)neigh.Size()){
+
+ if( CenterValue> neigh[i] & flag_max == 0) max = true;
+ else { max = false; flag_max = 1; }
+
+ if(CenterValue < neigh[i] & flag_min == 0) min = true;
+ else { min = false; flag_min = 1; }
+
+ ++i;
+ }
+
+ return max || min ;
+ }
+
+ /*-----------------------------------------------------------
+ * Compute the orientation of The KeyPoint
+ -----------------------------------------------------------*/
+
+ template <class TInputImage, class TOutputPointSet>
+ double
+ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
+ ::AssignOrientation(const NeighborhoodType& neigh , double S)
+ {
+
+ int i= 0 ;
+ int pas =( (i+S)-(int)(i+S) > 0.5 )?((int)S+1):(int)S ;
+ int Largeur = 2*neigh.GetRadius()[0]+1; // Width & length of a neighborhood
+ int rayon = neigh.GetRadius()[0]; // radius of the neigh
+ int col, raw ;
+ double dist ;
+ double w; // weight of the circular gaussian
+
+ OutputPointType pt ;
+
+ // Gradient orientation histogram
+ double angle;
+ int bin = 0 , Pi = 180;
+ int LengthBin = 60;
+ int NbBins = (2*Pi/LengthBin);
+ double tab[NbBins*2];
+
+
+ while (i < (int)neigh.Size())
+ {
+ col = i%Largeur - rayon ;
+ raw = i/Largeur - rayon ;
+ dist = sqrt(col *col + raw * raw );
+ col +=rayon;
+ raw +=rayon; // Backup to the image coordinate axes
+
+ if(dist < 6*S)
+ {
+ // Haar Wavelets responses accumulated in an histogram with Pi/3 precison
+ if (( col > pas && col < Largeur - pas ) && ( raw > pas && raw < Largeur - pas) )
+ {
+ w = vcl_exp(-((col-rayon)*(col-rayon) + (raw-rayon)*(raw-rayon))/(2*2.5*2.5*S*S) );
+ pt[0] = (neigh[(col+pas) + raw * Largeur] - neigh[(col-pas) + raw *Largeur ]) * w ;
+ pt[1] = (neigh[col + (raw+pas)* Largeur ] - neigh[col + (raw-pas)*Largeur]) * w;
+
+ if (pt[0] + pt[1] != 0)
+ {
+ angle = atan( pt[0]/pt[1] )*( Pi/3.1415);
+ if(angle < 0 )
+ angle += 2*Pi;
+
+ bin = (int)(angle/LengthBin);
+
+ if( bin <= NbBins-1 || bin >= 0 )
+ {
+ tab[2*bin] += pt[0];
+ tab[2*bin+1] += pt[1];
+ }
+ }
+ }
+ }
+ i+= pas;
+ }
+
+ //Find Orientation
+ double indice = 0;
+ double max = 0;
+ double length = 0;
+
+ //Detection de l'orientation du point courant
+ for (int i = 0 ; i < NbBins*2 ; i = i+2){
+ length = sqrt( tab[i]*tab[i] + tab[i+1]*tab[i+1] );
+ if( length > max){
+ max = length ;
+ indice = i/2;
+ }
+ }
+
+ return (indice+0.5)*LengthBin;
+ }
+
+ /*-----------------------------------------------------------
+ * Compute the descriptor of The KeyPoint
+ -----------------------------------------------------------*/
+
+ template <class TInputImage, class TOutputPointSet>
+ typename ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
+ ::VectorType
+ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
+ ::ComputeDescriptor(const NeighborhoodType& neigh , double O , double S )
+ {
+
+ typedef itk::CenteredRigid2DTransform<double> TransformType;
+ TransformType::Pointer eulerTransform = TransformType::New();
+ TransformType::ParametersType ParamVec(5);
+ PointImageType pSrc , pDst;
+ double angle = O * M_PI / 180;
+
+
+
+ int i = 0, col, raw , Nbin, pas = 1 ;
+ double xx = 0, yy = 0;
+ double dx, dy , w ;
+ int Largeur = 2*neigh.GetRadius()[0]+1;
+ double rayon = static_cast<double>(Largeur)/4. ;
+ double r = neigh.GetRadius()[0];
+ double dist = 0;
+ double x0 = neigh.GetCenterNeighborhoodIndex()% Largeur ;
+ double y0 = neigh.GetCenterNeighborhoodIndex()/ Largeur ;
+
+ //std::cout << " x0 " << x0 << " y0 " << y0 << angle << std::endl;
+
+ VectorType descriptorVector;
+ descriptorVector.resize(64);
+
+ /** Parameters of the transformation*/
+ ParamVec[0] = angle;
+ ParamVec[1] = x0;
+ ParamVec[2] = y0;
+ ParamVec[3] = 0;
+ ParamVec[4] = 0;
+ eulerTransform->SetParameters(ParamVec);
+
+ int accumulator = 0;
+ while (i < (int)neigh.Size())
+ {
+ col = i % Largeur ;
+ raw = i / Largeur ;
+
+ if (( col > pas && col < Largeur - pas ) && ( raw > pas && raw < Largeur - pas) )
+ {
+ double distanceX = (raw-r);
+ double distanceY = (col-r);
+ dist = vcl_sqrt(distanceX*distanceX + distanceY*distanceY);
+
+ if(dist <= r )
+ {
+ /* Transform point to compensate the rotation the orientation */
+ pDst[0] = col;
+ pDst[1] = raw;
+ pSrc = eulerTransform->TransformPoint(pDst);
+
+ /** New Coordinates (rotated) */
+ col = static_cast<int>(vcl_floor(pSrc[0]));
+ raw = static_cast<int>(vcl_floor(pSrc[1]));
+
+ if(raw==0) raw =+1;
+ if(col ==0) col +=1;
+
+ xx = static_cast<int> (pSrc[1]/rayon);
+ yy = static_cast<int> (pSrc[0]/rayon);
+
+ //Nbin = static_cast<int> ( vcl_floor(xx + 4 * yy) ) ;
+ Nbin = xx + 4*yy ;
+
+ if( Nbin < 16) //because 64 descriptor length
+ {
+ double distanceXcompensee_2 = (pSrc[0] - r)*(pSrc[0] - r);
+ double distanceYcompensee_2 = (pSrc[1] - r)*(pSrc[1] - r);
+
+ w = vcl_exp(-( distanceXcompensee_2 + distanceYcompensee_2 ) / (2*3.3*3.3*S*S) );
+ //w = vcl_exp(-( (pSrc[0] - r)*(pSrc[0] -r) + (pSrc[1] - r)*(pSrc[1] - r) )/(2*3.3*3.3/**S*S*/));
+
+ dx = 0.5 * (neigh[(col+pas) + raw * Largeur] - neigh[(col-pas) + raw *Largeur]) * w ;
+ dy = 0.5 * (neigh[col + (raw+ pas)* Largeur] - neigh[col + (raw-pas)*Largeur]) * w;
+
+ descriptorVector[4*Nbin ] += dx ;
+ descriptorVector[4*Nbin+1] += dy ;
+ descriptorVector[4*Nbin+2] += vcl_abs(dx) ;
+ descriptorVector[4*Nbin+3] += vcl_abs(dy) ;
+ }
+ else
+ {
+ accumulator++;
+ //std::cout << " xx "<< xx << " yy " << yy << " Nbin " << Nbin << " rayon "<< rayon << " pSrc[1] " << pSrc[1] << " pSrc[0] " <<pSrc[0] << " pDst[0] " << pDst[0] <<" pDst[1] " << pDst[1] << " Largeur "<< Largeur << std::endl;
+ }
+ }
+ }
+ i++;
+ }
+
+ if(accumulator > 0 ) std::cout << " accumulator "<< accumulator << " / "<< i << std::endl;
+
+
+ double accu = 0;
+ for (int i = 0 ; i < 64 ; i++)
+ accu += descriptorVector[i]*descriptorVector[i];
+
+ for (int j = 0 ; j < 64 ; j++)
+ descriptorVector[j] /= vcl_sqrt(accu) ;
+
+ return descriptorVector;
+
+ }
+
+ /*----------------------------------------------------------------
+ PrintSelf
+ -----------------------------------------------------------------*/
+ template <class TInputImage, class TOutputPointSet >
+ void
+ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
+ ::PrintSelf(std::ostream& os, itk::Indent indent) const
+ {
+ Superclass::PrintSelf(os, indent);
+ os << indent << "Number of Key Points " << m_NumberOfPoints << std::endl;
+ }
+}
diff --git a/Testing/Code/FeatureExtraction/CMakeLists.txt b/Testing/Code/FeatureExtraction/CMakeLists.txt
index 16147db51985b2cac75d32dd63de899be25edaba..f3ff009c402da4ef6f240a5750ed5c7762f9d046 100644
--- a/Testing/Code/FeatureExtraction/CMakeLists.txt
+++ b/Testing/Code/FeatureExtraction/CMakeLists.txt
@@ -839,10 +839,36 @@ ADD_TEST(feTvImageToHessianDeterminantImageFilter ${FEATUREEXTRACTION_TESTS9}
${TEMP}/feTvImageToHessianDeterminantImageFilterOutput.tif
1.5
)
-
-# ------- otb::ImageFittingPolygonListFilter -------------
-ADD_TEST(feTuImageFittingPolygonListFilterNew ${FEATUREEXTRACTION_TESTS9}
- otbImageFittingPolygonListFilterNew)
+
+# ------- otb::ImageToSURFKeyPointSetFilterNew -------------
+
+ADD_TEST(feTuImageToSURFKeyPointSetFilterNew ${FEATUREEXTRACTION_TESTS9}
+ otbImageToSURFKeyPointSetFilterNew)
+
+
+ADD_TEST(feTvImageToSURFKeyPointSetFilterSceneOutputImage ${FEATUREEXTRACTION_TESTS9}
+--compare-image ${EPSILON}
+ ${BASELINE}/feTvImageToSURFKeyPointSetFilterSceneImageOutput.png
+ ${TEMP}/feTvImageToSURFKeyPointSetFilterSceneImageOutput.png
+ otbImageToSURFKeyPointSetFilterOutputImage
+ ${INPUTDATA}/scene.png
+ ${TEMP}/feTvImageToSURFKeyPointSetFilterSceneImageOutput.png
+ 1 3
+)
+
+ADD_TEST(feTvImageToSURFKeyPointSetFilterSceneOutputAscii ${FEATUREEXTRACTION_TESTS9}
+--compare-ascii ${EPS}
+ ${BASELINE_FILES}/feTvImageToSURFKeyPointSetFilterSceneKeysOutput.txt
+ ${TEMP}/feTvImageToSURFKeyPointSetFilterSceneKeysOutput.txt
+ otbImageToSURFKeyPointSetFilterOutputAscii
+ ${INPUTDATA}/scene.png
+ ${TEMP}/feTvImageToSURFKeyPointSetFilterSceneKeysOutput.txt
+ 1 3
+)
+
+
+
+
# A enrichir
SET(BasicFeatureExtraction_SRCS1
@@ -948,7 +974,10 @@ otbImageToSIFTKeyPointSetFilterOutputImage.cxx
otbImageToSIFTKeyPointSetFilterOutputAscii.cxx
otbImageToHessianDeterminantImageFilterNew.cxx
otbImageToHessianDeterminantImageFilter.cxx
-otbImageFittingPolygonListFilterNew.cxx
+otbImageToSURFKeyPointSetFilterNew.cxx
+otbImageToSURFKeyPointSetFilterOutputImage.cxx
+otbImageToSURFKeyPointSetFilterOutputAscii.cxx
+#otbImageToSURFKeyPointSetFilterDistanceMap.cxx
)
INCLUDE_DIRECTORIES("${OTBTesting_BINARY_DIR}")
diff --git a/Testing/Code/FeatureExtraction/otbFeatureExtractionTests9.cxx b/Testing/Code/FeatureExtraction/otbFeatureExtractionTests9.cxx
index 9b4221b808cf1c489e92c082fbe1be9076307aeb..359428b9e59a911acbde11ab3faba4cebb9a84f7 100644
--- a/Testing/Code/FeatureExtraction/otbFeatureExtractionTests9.cxx
+++ b/Testing/Code/FeatureExtraction/otbFeatureExtractionTests9.cxx
@@ -35,6 +35,8 @@ REGISTER_TEST(otbImageToSIFTKeyPointSetFilterOutputAscii);
//REGISTER_TEST(otbImageToSIFTKeyPointSetFilterValid);
REGISTER_TEST(otbImageToHessianDeterminantImageFilterNew);
REGISTER_TEST(otbImageToHessianDeterminantImageFilter);
-REGISTER_TEST(otbImageFittingPolygonListFilterNew);
-
+REGISTER_TEST(otbImageToSURFKeyPointSetFilterNew);
+REGISTER_TEST(otbImageToSURFKeyPointSetFilterOutputImage);
+REGISTER_TEST(otbImageToSURFKeyPointSetFilterOutputAscii);
+//REGISTER_TEST(otbImageToSURFKeyPointSetFilterDistanceMap);
}
diff --git a/Testing/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilterDistanceMap.cxx b/Testing/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilterDistanceMap.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..7f39a02e92389ae529b5599456f3af491342ec46
--- /dev/null
+++ b/Testing/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilterDistanceMap.cxx
@@ -0,0 +1,346 @@
+/*=========================================================================
+
+Program: ORFEO Toolbox
+Language: C++
+Date: $Date$
+Version: $Revision$
+
+
+Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+See OTBCopyright.txt for details.
+
+Copyright (c) CS Systemes d'information. All rights reserved.
+See CSCopyright.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 <iostream>
+#include <fstream>
+
+#include "itkPointSet.h"
+#include "itkImageRegionIterator.h"
+#include "itkVariableLengthVector.h"
+#include "itkResampleImageFilter.h"
+#include "itkDanielssonDistanceMapImageFilter.h"
+#include "itkSubtractImageFilter.h"
+#include "itkMinimumMaximumImageCalculator.h"
+#include "itkShrinkImageFilter.h"
+#include "itkExpandImageFilter.h"
+#include "itkPointSetToImageFilter.h"
+#include "itkRescaleIntensityImageFilter.h"
+
+#include "otbImageToSURFKeyPointSetFilter.h"
+#include "otbImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+
+typedef itk::VariableLengthVector<float> RealVectorType;
+typedef itk::PointSet<RealVectorType,2> PointSetType;
+typedef otb::Image<float,2> ImageType;
+typedef otb::Image<unsigned char, 2> OutputImageType;
+
+// PointSet iterator types
+typedef PointSetType::PointsContainer PointsContainerType;
+typedef PointsContainerType::Iterator PointsIteratorType;
+typedef PointSetType::PointDataContainer PointDataContainerType;
+typedef PointDataContainerType::Iterator PointDataIteratorType;
+
+// Filter
+typedef otb::ImageFileReader<ImageType> ReaderType;
+typedef otb::ImageFileWriter<OutputImageType> WriterType;
+typedef otb::ImageFileWriter<ImageType> WriterInputType;
+
+OutputImageType::Pointer surf(ImageType::Pointer input,
+ const unsigned int octaves,
+ const unsigned int scales,
+ const char* surfFileName)
+{
+ typedef otb::ImageToSURFKeyPointSetFilter<ImageType,PointSetType> SurfFilterType;
+ typedef itk::PointSetToImageFilter<PointSetType, OutputImageType> PointSetFilterType;
+
+ SurfFilterType::Pointer surf = SurfFilterType::New();
+ PointSetFilterType::Pointer pointSetFilter = PointSetFilterType::New();
+
+ surf->SetInput(0,input);
+ surf->SetOctavesNumber(octaves);
+ surf->SetScalesNumber(scales);
+
+
+ pointSetFilter->SetInput(surf->GetOutput());
+ pointSetFilter->SetOutsideValue(0);
+ pointSetFilter->SetInsideValue(255);
+ pointSetFilter->SetSize(input->GetLargestPossibleRegion().GetSize());
+ pointSetFilter->SetSpacing(input->GetSpacing());
+ pointSetFilter->SetOrigin(input->GetOrigin());
+ pointSetFilter->Update();
+
+ WriterType::Pointer writer = WriterType::New();
+ writer->SetFileName(surfFileName);
+ writer->SetInput(pointSetFilter->GetOutput());
+ writer->Update();
+
+ return pointSetFilter->GetOutput();
+}
+
+OutputImageType::Pointer ddm( OutputImageType::Pointer input,
+ const char* ddmFileName)
+{
+ typedef itk::DanielssonDistanceMapImageFilter <OutputImageType, OutputImageType> DDMFilterType;
+ DDMFilterType::Pointer ddmFilter = DDMFilterType::New();
+
+ ddmFilter->SetInput(input);
+ ddmFilter->InputIsBinaryOn();
+ ddmFilter->Update();
+
+ WriterType::Pointer writer = WriterType::New();
+ writer->SetFileName(ddmFileName);
+ writer->SetInput(ddmFilter->GetOutput());
+ writer->Update();
+
+ return ddmFilter->GetOutput();
+}
+
+ImageType::Pointer rotate( ImageType::Pointer input,
+ const unsigned int rotation)
+{
+ typedef itk::AffineTransform<double, 2> TransformType;
+ typedef itk::ResampleImageFilter<ImageType, ImageType>
+ ResampleFilterType;
+ ResampleFilterType::Pointer resampler = ResampleFilterType::New();
+
+ TransformType::Pointer transform = TransformType::New();
+ TransformType::OutputVectorType translation1;
+ TransformType::OutputVectorType translation2;
+
+ const ImageType::SpacingType& spacing = input->GetSpacing();
+ const ImageType::PointType& origin = input->GetOrigin();
+ ImageType::SizeType size = input->GetLargestPossibleRegion().GetSize();
+
+ const double imageCenterX = origin[0] + spacing[0] * size[0] / 2.0;
+ const double imageCenterY = origin[1] + spacing[1] * size[1] / 2.0;
+ const double degreesToRadians = atan(1.0) / 45.0;
+ const double angle = rotation * degreesToRadians;
+
+ translation1[0] = -imageCenterX;
+ translation1[1] = -imageCenterY;
+ translation2[0] = imageCenterX;
+ translation2[1] = imageCenterY;
+
+ transform->Translate( translation1 );
+ transform->Rotate2D( -angle, false );
+ transform->Translate( translation2 );
+
+ resampler->SetOutputOrigin( origin );
+ resampler->SetOutputSpacing( spacing );
+ resampler->SetSize( size );
+ resampler->SetTransform(transform);
+ resampler->SetInput(input);
+ resampler->Update();
+ return resampler->GetOutput();
+}
+
+OutputImageType::Pointer invRotate(OutputImageType::Pointer input,
+ const unsigned int rotation)
+{
+ typedef itk::AffineTransform<double, 2> TransformType;
+ typedef itk::ResampleImageFilter<OutputImageType, OutputImageType>
+ ResampleFilterType;
+ ResampleFilterType::Pointer resampler = ResampleFilterType::New();
+
+ TransformType::Pointer transform = TransformType::New();
+ TransformType::OutputVectorType translation1;
+ TransformType::OutputVectorType translation2;
+
+ const ImageType::SpacingType& spacing = input->GetSpacing();
+ const ImageType::PointType& origin = input->GetOrigin();
+ ImageType::SizeType size = input->GetLargestPossibleRegion().GetSize();
+
+ const double imageCenterX = origin[0] + spacing[0] * size[0] / 2.0;
+ const double imageCenterY = origin[1] + spacing[1] * size[1] / 2.0;
+
+ const double degreesToRadians = atan(1.0) / 45.0;
+ const double angle = rotation * degreesToRadians;
+
+ translation1[0] = -imageCenterX;
+ translation1[1] = -imageCenterY;
+ translation2[0] = imageCenterX;
+ translation2[1] = imageCenterY;
+
+ transform->Translate( translation1 );
+ transform->Rotate2D( angle, false );
+ transform->Translate( translation2 );
+
+ resampler->SetOutputOrigin( origin );
+ resampler->SetOutputSpacing( spacing );
+ resampler->SetSize( size );
+ resampler->SetTransform(transform);
+ resampler->SetInput(input);
+ resampler->Update();
+ return resampler->GetOutput();
+}
+
+ImageType::Pointer zoom( ImageType::Pointer input,
+ const unsigned int zoomFactor)
+{
+ typedef itk::ShrinkImageFilter<ImageType, ImageType> ShrinkFilterType;
+ ShrinkFilterType::Pointer shrink = ShrinkFilterType::New();
+
+ shrink->SetInput(input);
+ shrink->SetShrinkFactors(zoomFactor);
+ shrink->Update();
+
+ return shrink->GetOutput();
+}
+
+OutputImageType::Pointer invZoom( OutputImageType::Pointer input,
+ const unsigned int zoomFactor)
+{
+ typedef itk::ExpandImageFilter<OutputImageType, OutputImageType> ExpandFilterType;
+ ExpandFilterType::Pointer expand = ExpandFilterType::New();
+
+ expand->SetInput(input);
+ expand->SetExpandFactors(zoomFactor);
+ expand->Update();
+
+ return expand->GetOutput();
+}
+
+ImageType::Pointer contrast(ImageType::Pointer input,
+ const unsigned int contrastMin,
+ const unsigned int contrastMax)
+{
+ typedef itk::RescaleIntensityImageFilter<ImageType, ImageType> RescaleFilterType;
+ RescaleFilterType::Pointer rescaler= RescaleFilterType::New();
+
+ rescaler->SetInput(input);
+ rescaler->SetOutputMinimum(static_cast<RescaleFilterType::OutputPixelType>(contrastMin));
+ rescaler->SetOutputMaximum(static_cast<RescaleFilterType::OutputPixelType>(contrastMax));
+ rescaler->Update();
+ return rescaler->GetOutput();
+}
+
+OutputImageType::Pointer invContrast( OutputImageType::Pointer input,
+ const unsigned int contrastMin,
+ const unsigned int contrastMax)
+{
+ return input;
+}
+
+void subtract(OutputImageType::Pointer image1,
+ OutputImageType::Pointer image2,
+ const char* subtractFileName)
+{
+ typedef itk::SubtractImageFilter<OutputImageType, OutputImageType, ImageType> SubtractFilterType;
+ typedef itk::MinimumMaximumImageCalculator<ImageType> MaximumCalculatorType;
+ typedef itk::RescaleIntensityImageFilter<ImageType, OutputImageType> OutputRescaleFilterType;
+
+ SubtractFilterType::Pointer subtract = SubtractFilterType::New();
+ MaximumCalculatorType::Pointer maximumCalculator = MaximumCalculatorType::New();
+
+ subtract->SetInput1(image1);
+ subtract->SetInput2(image2);
+ subtract->Update();
+
+ OutputRescaleFilterType::Pointer rescaler = OutputRescaleFilterType::New();
+ rescaler->SetInput(subtract->GetOutput());
+ rescaler->SetOutputMinimum(0);
+ rescaler->SetOutputMaximum(255);
+
+ WriterType::Pointer writer = WriterType::New();
+ writer->SetFileName(subtractFileName);
+ writer->SetInput(rescaler->GetOutput());
+ writer->Update();
+
+ maximumCalculator->SetImage(subtract->GetOutput());
+ maximumCalculator->Compute();
+
+ std::cout << "Mix(sub)= " << maximumCalculator->GetMinimum() << " ";
+ std::cout << "Max(sub)= " << maximumCalculator->GetMaximum() << std::endl;
+}
+
+int otbImageToSURFKeyPointSetFilterDistanceMap(int argc, char * argv[])
+{
+ // Input Image file name
+ const char * infname = argv[1];
+
+ const unsigned int octaves = atoi(argv[2]);
+ const unsigned int scales = atoi(argv[3]);
+
+ // Rotation angle [0,360[
+ const unsigned int rotation = atoi(argv[6]);
+
+ // Zoom factor
+ const unsigned int zoomFactor = atoi(argv[7]);
+
+ // contrast factor
+ const unsigned int contrastMin = atoi(argv[8]);
+ const unsigned int contrastMax = atoi(argv[9]);
+
+ ReaderType::Pointer reader = ReaderType::New();
+ reader->SetFileName(infname);
+ reader->Update();
+
+ ImageType::Pointer _rotated =
+ rotate(reader->GetOutput(), rotation);
+ ImageType::Pointer _zoomed =
+ zoom(reader->GetOutput(), zoomFactor);
+ ImageType::Pointer _contrasted =
+ contrast(reader->GetOutput(), contrastMin, contrastMax);
+
+ ImageType::Pointer _combined1 = zoom(_rotated, zoomFactor);
+ ImageType::Pointer _combined2 = contrast(_combined1, contrastMin, contrastMax);
+
+ OutputImageType::Pointer surf_base =
+ surf(reader->GetOutput(), octaves, scales, "surf_base.png");
+
+ OutputImageType::Pointer _surf_rotated =
+ surf(_rotated, octaves, scales ,"surf_rotated.png");
+
+ OutputImageType::Pointer _surf_zoomed =
+ surf(_zoomed, octaves, scales, "surf_zoomed.png");
+
+ OutputImageType::Pointer _surf_contrasted =
+ surf(_contrasted, octaves, scales, "surf_contrasted.png");
+
+ OutputImageType::Pointer _surf_combined =
+ surf(_combined2, octaves, scales, "surf_combined.png");
+
+ OutputImageType::Pointer surf_rotated =
+ invRotate(_surf_rotated, rotation);
+
+ OutputImageType::Pointer surf_zoomed =
+ invZoom(_surf_zoomed, zoomFactor);
+
+ OutputImageType::Pointer surf_contrasted =
+ invContrast(_surf_contrasted, contrastMin, contrastMax);
+
+ OutputImageType::Pointer _surf_combined1 =
+ invContrast(_surf_combined, contrastMin, contrastMax);
+ OutputImageType::Pointer _surf_combined2 =
+ invRotate(_surf_combined1, rotation);
+ OutputImageType::Pointer surf_combined =
+ invZoom(_surf_combined2, zoomFactor);
+
+ OutputImageType::Pointer ddm_base = ddm(surf_base, "ddm_base.png");
+ OutputImageType::Pointer ddm_rotated = ddm(surf_rotated, "ddm_rotated.png");
+ OutputImageType::Pointer ddm_contrasted = ddm(surf_contrasted, "ddm_contrasted.png");
+ OutputImageType::Pointer ddm_zoomed = ddm(surf_zoomed, "ddm_zoomed.png");
+ OutputImageType::Pointer ddm_combined = ddm(surf_combined, "ddm_combined.png");
+
+ subtract(ddm_base, ddm_rotated,
+ "subtract_rotated.png");
+
+ subtract(ddm_base, ddm_zoomed,
+ "subtract_zoomed.png");
+
+ subtract(ddm_base, ddm_contrasted,
+ "subtract_contrasted.png");
+
+ subtract(ddm_base, ddm_combined,
+ "subtract_combined.png");
+
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilterNew.cxx b/Testing/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilterNew.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..17907339e5c1b16d72e59378465dd9622c6b8c87
--- /dev/null
+++ b/Testing/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilterNew.cxx
@@ -0,0 +1,41 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+ Copyright (c) CS Systemes d'information. All rights reserved.
+ See CSCopyright.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 "otbImageToSURFKeyPointSetFilter.h"
+#include "otbImage.h"
+#include "itkPointSet.h"
+#include "itkVariableLengthVector.h"
+
+
+int otbImageToSURFKeyPointSetFilterNew(int argc, char * argv[])
+{
+ const unsigned int Dimension = 2;
+ typedef double PixelType;
+
+ typedef otb::Image<PixelType,Dimension> ImageType;
+ typedef itk::VariableLengthVector<PixelType> RealVectorType;
+ typedef itk::PointSet<RealVectorType,Dimension> PointSetType;
+ typedef otb::ImageToSURFKeyPointSetFilter<ImageType,PointSetType> FilterType;
+
+ // Instantiating object
+ FilterType::Pointer object = FilterType::New();
+
+
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilterOutputAscii.cxx b/Testing/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilterOutputAscii.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..1397137f01b4c026a3c82ed2aa5a91821acf6db0
--- /dev/null
+++ b/Testing/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilterOutputAscii.cxx
@@ -0,0 +1,102 @@
+/*=========================================================================
+
+Program: ORFEO Toolbox
+Language: C++
+Date: $Date$
+Version: $Revision$
+
+
+Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+See OTBCopyright.txt for details.
+
+Copyright (c) CS Systemes d'information. All rights reserved.
+See CSCopyright.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 "otbImageToSURFKeyPointSetFilter.h"
+#include "otbImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+#include "itkPointSet.h"
+#include "itkVariableLengthVector.h"
+#include "otbRationalQuotientResampleImageFilter.h"
+#include "itkRGBPixel.h"
+#include "itkImageRegionIterator.h"
+
+#include <iostream>
+#include <fstream>
+
+int otbImageToSURFKeyPointSetFilterOutputAscii(int argc, char * argv[])
+{
+
+ if(argc < 5 )
+ {
+ std::cout << " Usage : otbSURFTest imageName FileOutName Octave[int] Level[int]" << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ const char * infname = argv[1];
+ const char * outfname = argv[2];
+
+ const unsigned int octaves = atoi(argv[3]);
+ const unsigned int scales = atoi(argv[4]);
+
+
+ typedef float RealType;
+ const unsigned int Dimension =2;
+
+ typedef otb::Image<RealType,Dimension> ImageType;
+ typedef itk::VariableLengthVector<RealType> RealVectorType;
+ typedef otb::ImageFileReader<ImageType> ReaderType;
+ typedef itk::PointSet<RealVectorType,Dimension> PointSetType;
+ typedef otb::ImageToSURFKeyPointSetFilter<ImageType,PointSetType> ImageToSURFKeyPointSetFilterType;
+
+ // PointSet iterator types
+ typedef PointSetType::PointsContainer PointsContainerType;
+ typedef PointsContainerType::Iterator PointsIteratorType;
+ typedef PointSetType::PointDataContainer PointDataContainerType;
+ typedef PointDataContainerType::Iterator PointDataIteratorType;
+
+ // Instantiating object
+ ReaderType::Pointer reader = ReaderType::New();
+ ImageToSURFKeyPointSetFilterType::Pointer filter = ImageToSURFKeyPointSetFilterType::New();
+
+ reader->SetFileName(infname);
+ filter->SetInput(0,reader->GetOutput());
+ filter->SetOctavesNumber(octaves);
+ filter->SetScalesNumber(scales);
+ filter->Update();
+
+ PointsIteratorType pIt = filter->GetOutput()->GetPoints()->Begin();
+ PointDataIteratorType pDataIt = filter->GetOutput()->GetPointData()->Begin();
+
+ std::ofstream outfile(outfname);
+
+ outfile << "Number of octaves: "<<octaves << std::endl;
+ outfile << "Number of scales: "<<scales << std::endl;
+ outfile << "Number of SURF key points: " << filter->GetNumberOfPoints() << std::endl;
+
+ while( pIt!=filter->GetOutput()->GetPoints()->End() )
+ {
+ outfile << "[" << std::fixed << std::setprecision(2) << pIt.Value()[0] << ", " << std::setprecision(2) << pIt.Value()[1] << "][";
+
+ unsigned int lIterDesc=0;
+
+ while (lIterDesc < pDataIt.Value().Size())
+ {
+ outfile << std::setprecision(3) << pDataIt.Value()[lIterDesc] << " ";
+ lIterDesc++;
+ }
+ outfile << "]" << std::endl;
+ ++pIt;
+ ++pDataIt;
+ }
+
+ outfile.close();
+
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilterOutputImage.cxx b/Testing/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilterOutputImage.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..7cb1117af47e1a2313214760c93f1b6c812bd66e
--- /dev/null
+++ b/Testing/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilterOutputImage.cxx
@@ -0,0 +1,167 @@
+/*=========================================================================
+
+Program: ORFEO Toolbox
+Language: C++
+Date: $Date$
+Version: $Revision$
+
+
+Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+See OTBCopyright.txt for details.
+
+Copyright (c) CS Systemes d'information. All rights reserved.
+See CSCopyright.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 "otbImageToSURFKeyPointSetFilter.h"
+#include "otbImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+#include "itkPointSet.h"
+#include "itkVariableLengthVector.h"
+#include "otbRationalQuotientResampleImageFilter.h"
+#include "itkRGBPixel.h"
+#include "itkImageRegionIterator.h"
+
+#include <iostream>
+#include <fstream>
+
+int otbImageToSURFKeyPointSetFilterOutputImage(int argc, char * argv[])
+{
+ const char * infname = argv[1];
+ const char * outputImageFilename = argv[2];
+
+ const unsigned int octaves = atoi(argv[3]);
+ const unsigned int scales = atoi(argv[4]);
+
+ typedef float RealType;
+ const unsigned int Dimension =2;
+
+ typedef otb::Image<RealType,Dimension> ImageType;
+ typedef itk::VariableLengthVector<RealType> RealVectorType;
+ typedef otb::ImageFileReader<ImageType> ReaderType;
+ typedef itk::PointSet<RealVectorType,Dimension> PointSetType;
+ typedef otb::ImageToSURFKeyPointSetFilter<ImageType,PointSetType> ImageToSURFKeyPointSetFilterType;
+
+ // PointSet iterator types
+ typedef PointSetType::PointsContainer PointsContainerType;
+ typedef PointsContainerType::Iterator PointsIteratorType;
+
+ // Instantiating object
+ ReaderType::Pointer reader = ReaderType::New();
+ ImageToSURFKeyPointSetFilterType::Pointer filter = ImageToSURFKeyPointSetFilterType::New();
+
+ reader->SetFileName(infname);
+ filter->SetInput(0,reader->GetOutput());
+ filter->SetOctavesNumber(octaves);
+ filter->SetScalesNumber(scales);
+
+ filter->Update();
+
+ ImageType::OffsetType t = {{ 0, 1}};
+ ImageType::OffsetType b = {{ 0,-1}};
+ ImageType::OffsetType l = {{ 1, 0}};
+ ImageType::OffsetType r = {{-1, 0}};
+
+ typedef unsigned char PixelType;
+ typedef otb::Image<PixelType,2> UCharImageType;
+
+ typedef itk::RGBPixel<PixelType> RGBPixelType;
+ typedef otb::Image<RGBPixelType, 2> OutputImageType;
+
+ typedef otb::ImageFileWriter<OutputImageType> WriterType;
+ OutputImageType::Pointer outputImage = OutputImageType::New();
+ OutputImageType::RegionType region;
+
+ OutputImageType::SizeType outputSize;
+ outputSize[0] = reader->GetOutput()->GetLargestPossibleRegion().GetSize()[0];
+ outputSize[1] = reader->GetOutput()->GetLargestPossibleRegion().GetSize()[1];
+ region.SetSize(outputSize);
+
+ OutputImageType::IndexType indexStart;
+ indexStart[0] = 0;
+ indexStart[1] = 0;
+ region.SetIndex(indexStart);
+
+ outputImage->SetRegions(region);
+ outputImage->Allocate();
+
+ itk::ImageRegionIterator<OutputImageType> iterOutput(outputImage,
+ outputImage->GetLargestPossibleRegion());
+ itk::ImageRegionIterator<ImageType> iterInput(reader->GetOutput(),
+ reader->GetOutput()->GetLargestPossibleRegion());
+
+ for (iterOutput.GoToBegin(), iterInput.GoToBegin();
+ !iterOutput.IsAtEnd();
+ ++iterOutput, ++iterInput)
+ {
+ OutputImageType::PixelType rgbPixel;
+ rgbPixel.SetRed( static_cast<PixelType>(iterInput.Get()) );
+ rgbPixel.SetGreen( static_cast<PixelType>(iterInput.Get()) );
+ rgbPixel.SetBlue( static_cast<PixelType>(iterInput.Get()) );
+
+ iterOutput.Set(rgbPixel);
+ }
+
+ WriterType::Pointer writerTmp = WriterType::New();
+ writerTmp->SetFileName(outputImageFilename);
+ writerTmp->SetInput(outputImage);
+ writerTmp->Update();
+
+ std::cout << "Copy Input image in Output image" << std::endl;
+
+ PointsIteratorType pIt = filter->GetOutput()->GetPoints()->Begin();
+ ImageType::SpacingType spacing = reader->GetOutput()->GetSpacing();
+ ImageType::PointType origin = reader->GetOutput()->GetOrigin();
+ OutputImageType::SizeType size = outputImage->GetLargestPossibleRegion().GetSize();
+
+ while( pIt != filter->GetOutput()->GetPoints()->End() )
+ {
+ ImageType::IndexType index;
+
+ index[0] = (unsigned int)
+ (vcl_floor
+ ((double)((pIt.Value()[0]-origin[0])/spacing[0]+0.5)));
+
+ index[1] = (unsigned int)
+ (vcl_floor
+ ((double)((pIt.Value()[1]-origin[1])/spacing[1]+0.5)));
+
+ OutputImageType::PixelType keyPixel;
+ keyPixel.SetRed(0);
+ keyPixel.SetGreen(255);
+ keyPixel.SetBlue(0);
+
+ if (outputImage->GetLargestPossibleRegion().IsInside(index))
+ {
+ outputImage->SetPixel(index,keyPixel);
+
+ if (outputImage->GetLargestPossibleRegion().IsInside(index+t))
+ outputImage->SetPixel(index+t,keyPixel);
+
+ if (outputImage->GetLargestPossibleRegion().IsInside(index+b))
+ outputImage->SetPixel(index+b,keyPixel);
+
+ if (outputImage->GetLargestPossibleRegion().IsInside(index+l))
+ outputImage->SetPixel(index+l,keyPixel);
+
+ if (outputImage->GetLargestPossibleRegion().IsInside(index+r))
+ outputImage->SetPixel(index+r,keyPixel);
+ }
+ ++pIt;
+ }
+
+ std::cout << "Copy sift key" << std::endl;
+
+ WriterType::Pointer writer = WriterType::New();
+ writer->SetFileName(outputImageFilename);
+ writer->SetInput(outputImage);
+ writer->Update();
+
+ std::cout << "Write image" << std::endl;
+ return EXIT_SUCCESS;
+}