diff --git a/Code/BasicFilters/otbMeanShiftVectorImageFilter.h b/Code/BasicFilters/otbMeanShiftVectorImageFilter.h
new file mode 100644
index 0000000000000000000000000000000000000000..cbf621c67ad1698d450a8ec8f6eca65ba86a8023
--- /dev/null
+++ b/Code/BasicFilters/otbMeanShiftVectorImageFilter.h
@@ -0,0 +1,121 @@
+/*=========================================================================
+
+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 _otbMeanShiftVectorImageFilter_h
+#define _otbMeanShiftVectorImageFilter_h
+
+#include "itkInPlaceImageFilter.h"
+
+namespace otb
+{
+
+ /** \class MeanShiftVectorImageFilter
+ *
+ *
+ *
+ *
+ * \ingroup Streamed
+ * \ingroup Threaded
+ */
+
+template <class TInputImage, class TOutputImage>
+class ITK_EXPORT MeanShiftVectorImageFilter
+ : public itk::InPlaceImageFilter<TInputImage,TOutputImage>
+ {
+ public:
+ /** Standard class typedef */
+ typedef MeanShiftVectorImageFilter Self;
+ typedef itk::InPlaceImageFilter<TInputImage,TOutputImage> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** New and Type macros */
+ itkNewMacro(Self);
+ itkTypeMacro(MeanShiftVectorImageFilter,InPlaceImageFilter);
+
+ /** Template parameters typedefs */
+ typedef TInputImage InputImageType;
+ typedef typename InputImageType::Pointer InputImagePointerType;
+ typedef typename InputImageType::PixelType InputPixelType;
+ typedef typename InputImageType::InternalPixelType InputInternalPixelType;
+ typedef typename InputImageType::PointType PointType;
+ typedef TOutputImage OutputImageType;
+ typedef typename OutputImageType::Pointer OutputImagePointerType;
+ typedef typename OutputImageType::PixelType OutputPixelType;
+ typedef typename OutputImageType::InternalPixelType OutputInternalPixelType;
+ typedef typename OutputImageType::RegionType RegionType;
+ typedef typename RegionType::SizeType SizeType;
+
+ /** Setters / Getters */
+ itkSetMacro(SpatialRadius,double);
+ itkGetMacro(SpatialRadius,double);
+ itkSetMacro(RangeRadius,double);
+ itkGetMacro(RangeRadius,double);
+ itkSetMacro(MaxNumberOfIterations,unsigned int);
+ itkGetMacro(MaxNumberOfIterations,unsigned int);
+ itkGetMacro(ConvergenceDistanceThreshold,double);
+ itkSetMacro(ConvergenceDistanceThreshold,double);
+ itkSetMacro(UseImageSpacing,bool);
+ itkGetMacro(UseImageSpacing,bool);
+ itkBooleanMacro(UseImageSpacing);
+
+ protected:
+ /** This filters use a neighborhood around the pixel, so it needs to redfine the
+ * input requested region */
+ virtual void GenerateInputRequestedRegion();
+
+ /** Threaded generate data */
+ virtual void ThreadedGenerateData(const RegionType & outputRegionForThread,int threadId);
+
+ /** Constructor */
+ MeanShiftVectorImageFilter();
+
+ /** destructor */
+ ~MeanShiftVectorImageFilter(){};
+
+ /**PrintSelf method */
+ virtual void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+ private:
+ MeanShiftVectorImageFilter(const Self&); //purposely not implemented
+ void operator=(const Self&); //purposely not implemented
+
+ /** Spatial radius for mean shift convergence */
+ double m_SpatialRadius;
+ /** Range radius for mean shift convergence */
+ double m_RangeRadius;
+ /** Max number of iterations for convergence */
+ unsigned int m_MaxNumberOfIterations;
+
+ /** Internal radius used by the iterator. The value depends on the spatial radius
+ * and the maximum number of iterations. This value is not intended to be accessed
+ * by users, thus no getter or setter is provided */
+ SizeType m_InternalRadius;
+
+ /** use image spacing */
+ bool m_UseImageSpacing;
+
+ /** Distance threshold for convergence in the spatial domain */
+ double m_ConvergenceDistanceThreshold;
+ };
+}// end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbMeanShiftVectorImageFilter.txx"
+#endif
+
+#endif
diff --git a/Code/BasicFilters/otbMeanShiftVectorImageFilter.txx b/Code/BasicFilters/otbMeanShiftVectorImageFilter.txx
new file mode 100644
index 0000000000000000000000000000000000000000..3f124649bb82c44c0a650851f81f005c7eee8022
--- /dev/null
+++ b/Code/BasicFilters/otbMeanShiftVectorImageFilter.txx
@@ -0,0 +1,263 @@
+/*=========================================================================
+
+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 _otbMeanShiftVectorImageFilter_txx
+#define _otbMeanShiftVectorImageFilter_txx
+
+#include "otbMeanShiftVectorImageFilter.h"
+
+#include "itkConstNeighborhoodIterator.h"
+#include "itkImageRegionIterator.h"
+#include "itkNeighborhoodAlgorithm.h"
+#include "itkZeroFluxNeumannBoundaryCondition.h"
+#include "itkOffset.h"
+#include "itkProgressReporter.h"
+
+#include "otbMacro.h"
+
+namespace otb
+{
+ template <class TInputImage,class TOutputImage>
+ MeanShiftVectorImageFilter<TInputImage,TOutputImage>
+ ::MeanShiftVectorImageFilter()
+ {
+ m_MaxNumberOfIterations = 20;
+ m_ConvergenceDistanceThreshold = 0.5;
+ m_SpatialRadius = 3;
+ m_RangeRadius = 10;
+ m_UseImageSpacing = false;
+ m_InternalRadius.Fill(0);
+ }
+
+
+ template <class TInputImage,class TOutputImage>
+ void
+ MeanShiftVectorImageFilter<TInputImage,TOutputImage>
+ ::GenerateInputRequestedRegion()
+ {
+ // call the superclass' implementation of this method
+ Superclass::GenerateInputRequestedRegion();
+
+ // get pointers to the input and output
+ typename Superclass::InputImagePointer inputPtr =
+ const_cast< TInputImage * >( this->GetInput() );
+ typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
+
+ if ( !inputPtr || !outputPtr )
+ {
+ return;
+ }
+
+ // get a copy of the input requested region (should equal the output
+ // requested region)
+ typename TInputImage::RegionType inputRequestedRegion;
+ inputRequestedRegion = inputPtr->GetRequestedRegion();
+
+ //Compute the radius we will need to execute the filter:
+ if(m_UseImageSpacing)
+ {
+ for(unsigned int dimension = 0; dimension < InputImageType::ImageDimension;++dimension)
+ {
+ m_InternalRadius[dimension] = static_cast<unsigned long>(vcl_floor(m_MaxNumberOfIterations * m_SpatialRadius / inputPtr->GetSpacing()[0]+0.5));
+ }
+ }
+ else
+ {
+ m_InternalRadius.Fill(static_cast<unsigned long>(vcl_floor(m_MaxNumberOfIterations * m_SpatialRadius+0.5)));
+ }
+ std::cout<<"MeanShiftVectorImageFilter: Internal radius = "<<m_InternalRadius<<std::endl;
+
+
+ // pad the input requested region by the operator radius
+ inputRequestedRegion.PadByRadius( m_InternalRadius );
+
+ // crop the input requested region at the input's largest possible region
+ if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
+ {
+ inputPtr->SetRequestedRegion( inputRequestedRegion );
+ return;
+ }
+ else
+ {
+ // Couldn't crop the region (requested region is outside the largest
+ // possible region). Throw an exception.
+
+ // store what we tried to request (prior to trying to crop)
+ inputPtr->SetRequestedRegion( inputRequestedRegion );
+
+ // build an exception
+ itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
+ e.SetLocation(ITK_LOCATION);
+ e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
+ e.SetDataObject(inputPtr);
+ throw e;
+ }
+ }
+
+ template <class TInputImage,class TOutputImage>
+ void
+ MeanShiftVectorImageFilter<TInputImage,TOutputImage>
+ ::ThreadedGenerateData(const RegionType & outputRegionForThread,int threadId)
+ {
+ std::cout<<"Call to threaded generate data, threadId: "<<threadId<<", region: "<<outputRegionForThread<<std::endl;
+ // TODO: update this to handle the UseImageSpacing option.
+
+ // Set up the progress reporter
+ itk::ProgressReporter progress(this,threadId,outputRegionForThread.GetNumberOfPixels());
+
+ // Input and output pointers
+ typename OutputImageType::Pointer outputPtr = this->GetOutput();
+ typename InputImageType::ConstPointer inputPtr = this->GetInput();
+
+ // Iterators
+ itk::ConstNeighborhoodIterator<InputImageType> inputIt(m_InternalRadius,inputPtr,outputRegionForThread);
+ itk::ConstNeighborhoodIterator<InputImageType> internalIt;
+ itk::ImageRegionIterator<OutputImageType> outputIt(outputPtr,outputRegionForThread);
+
+ // Boundary conditions
+ itk::ZeroFluxNeumannBoundaryCondition<InputImageType> nbc;
+ inputIt.OverrideBoundaryCondition(&nbc);
+
+ //TODO: Work with a floating point pixel type
+
+ // local declarations
+ unsigned int nbIterations;
+ typedef typename itk::ConstNeighborhoodIterator<InputImageType>::OffsetType OffsetType;
+
+ OffsetType offset;
+ PointType maxDensityPoint,convergencePoint;
+ InputPixelType maxDensityValue, convergenceValue,diff,current;
+ unsigned int nbPixelsIntoAccount;
+
+ unsigned int nbComponentsPerPixel = inputPtr->GetNumberOfComponentsPerPixel();
+ bool goesOn = true;
+ long startx,stopx,starty,stopy;
+
+ // intialize iterators
+ inputIt.GoToBegin();
+ outputIt.GoToBegin();
+
+ // Walk the images
+ while(!inputIt.IsAtEnd() && !outputIt.IsAtEnd())
+ {
+ nbIterations = 0;
+ goesOn = true;
+ convergencePoint.Fill(0);
+ convergenceValue = inputIt.GetCenterPixel();
+ //std::cout<<"Processing pixel: "<<inputIt.GetIndex()<<", value: "<<convergenceValue<<std::endl;
+
+
+ // While the max number of iterations has not been reached and the convergence index is still moving
+ while(nbIterations < m_MaxNumberOfIterations && goesOn)
+ {
+ //std::cout<<"Loop: "<<nbIterations<<", goesOn: "<<goesOn<<std::endl;
+ maxDensityPoint.Fill(0);
+ maxDensityValue.SetSize(nbComponentsPerPixel);
+ maxDensityValue.Fill(0);
+ nbPixelsIntoAccount = 0;
+
+ startx = vcl_floor(convergencePoint[0]-m_SpatialRadius);
+ stopx = vcl_ceil(convergencePoint[0]+m_SpatialRadius);
+ starty = vcl_floor(convergencePoint[1]-m_SpatialRadius);
+ stopy = vcl_ceil(convergencePoint[1]+m_SpatialRadius);
+
+ // loop on the neighborhood
+ for(long i = startx;i<=stopx;++i)
+ {
+ for(long j =starty;j<=stopy;++j)
+ {
+ // Setting the offset
+ offset[0]=i;
+ offset[1]=j;
+
+ // verify that the current offset is inside or spatial circle
+ if(vcl_sqrt(vcl_pow(static_cast<double>(convergencePoint[0]-offset[0]),2)+vcl_pow(static_cast<double>(convergencePoint[1]-offset[1]),2))<m_SpatialRadius)
+ {
+ //std::cout<<"Point: "<<offset<<" is inside spatial circle"<<std::endl;
+ current = inputIt.GetPixel(offset);
+ diff = convergenceValue - current;
+ //std::cout<<"Difference with current pixel: "<<diff<<std::endl;
+ // Check if we are inside the spectral - spatial circle
+ if(diff.GetSquaredNorm()+vcl_pow(static_cast<double>(offset[0]),2)+vcl_pow(static_cast<double>(offset[1]),2)<nbComponentsPerPixel * vcl_pow(m_RangeRadius,2) + vcl_pow(m_SpatialRadius,2))
+ {
+ //std::cout<<"Point: "<<offset<<" is inside the spatial/spectral disk"<<std::endl;
+ maxDensityPoint[0]+=offset[0];
+ maxDensityPoint[1]+=offset[1];
+ maxDensityValue+=current;
+ ++nbPixelsIntoAccount;
+ }
+ }
+ }
+ }
+ // If no inside pixels
+ if(nbPixelsIntoAccount == 0)
+ {
+ maxDensityPoint.Fill(0);
+ maxDensityValue = inputIt.GetCenterPixel();
+ }
+ else
+ {
+ maxDensityPoint[0]/=nbPixelsIntoAccount;
+ maxDensityPoint[1]/=nbPixelsIntoAccount;
+ maxDensityValue/=nbPixelsIntoAccount;
+ }
+ //std::cout<<"MaxDensityValue: "<<maxDensityValue<<std::endl;
+ convergenceValue = maxDensityValue;
+
+ //std::cout<<"Nb pixels into account: "<<nbPixelsIntoAccount<<", convergence point: "<<convergencePoint<<", maxDensityPoint: "<<maxDensityPoint<<", convergenceValue: "<<convergenceValue<<std::endl;
+
+ // Check if we are still significantly moving
+ double updatedDistance = vcl_sqrt(vcl_pow(static_cast<double>(maxDensityPoint[0]-convergencePoint[0]),2)+vcl_pow(static_cast<double>(maxDensityPoint[1]-convergencePoint[1]),2));
+ if(updatedDistance < m_ConvergenceDistanceThreshold)
+ {
+ goesOn = false;
+ }
+ else
+ {
+ convergencePoint = maxDensityPoint;
+ }
+ ++nbIterations;
+ }
+ // Set the output value
+ // TODO: missing value cast here
+ outputIt.Set(convergenceValue);
+ // Update progress
+ progress.CompletedPixel();
+ // Increment iterators
+ ++inputIt;
+ ++outputIt;
+ }
+ }
+
+ template <class TInputImage,class TOutputImage>
+ void
+ MeanShiftVectorImageFilter<TInputImage,TOutputImage>
+ ::PrintSelf(std::ostream& os, itk::Indent indent) const
+ {
+ Superclass::PrintSelf(os,indent);
+ os<<indent<<"Spatial radius: "<<m_SpatialRadius<<std::endl;
+ os<<indent<<"Range radius: "<<m_RangeRadius<<std::endl;
+ os<<indent<<"Maximum number of iterations: "<<m_MaxNumberOfIterations<<std::endl;
+ os<<indent<<"Internal radius: "<<m_InternalRadius<<std::endl;
+ os<<indent<<"Use image spacing: "<<m_UseImageSpacing<<std::endl;
+ os<<indent<<"Convergence distance threshold: "<<m_ConvergenceDistanceThreshold<<std::endl;
+
+ }
+
+} // end namespace otb
+
+#endif
diff --git a/Testing/Code/BasicFilters/CMakeLists.txt b/Testing/Code/BasicFilters/CMakeLists.txt
index f532cd7bf9057e6e0961b50cfef3b6982a877394..1d261725ae5781bc51100b5728017f8f033a5aab 100755
--- a/Testing/Code/BasicFilters/CMakeLists.txt
+++ b/Testing/Code/BasicFilters/CMakeLists.txt
@@ -883,13 +883,15 @@ ADD_TEST(bfTvContinuousMinimumMaximumImageCalculatorTest ${BASICFILTERS_TESTS9}
# 2 2
# )
-#ADD_TEST(bfTuMeanShiftImageFilterNew ${BASICFILTERS_TESTS9}
-# otbMeanShiftImageFilterNew )
-
-#ADD_TEST(bfTvMeanShiftImageFilter ${BASICFILTERS_TESTS9}
-# otbMeanShiftImageFilter )
-
-
+ADD_TEST(bfTuMeanShiftVectorImageFilterNew ${BASICFILTERS_TESTS9}
+ otbMeanShiftVectorImageFilterNew )
+
+ADD_TEST(bfTvMeanShiftVectorImageFilter ${BASICFILTERS_TESTS9}
+ otbMeanShiftVectorImageFilter
+ ${INPUTDATA}/qb_RoadExtract2sub200x200.tif
+ ${TEMP}/bfTvvMeanShiftVectorImageFilterOutput.tif
+ 3 15 20 0 0.0001
+ )
ADD_TEST(bfTvFunctionToImageFilterNew ${BASICFILTERS_TESTS9}
otbFunctionToImageFilterNew
@@ -1054,8 +1056,8 @@ otbContinuousMinimumMaximumImageCalculatorNew.cxx
otbContinuousMinimumMaximumImageCalculatorTest.cxx
#MeanShift.cxx
#MeanShiftVectorImage.cxx
-#otbMeanShiftImageFilterNew.cxx
-#otbMeanShiftImageFilter.cxx
+otbMeanShiftVectorImageFilterNew.cxx
+otbMeanShiftVectorImageFilter.cxx
otbFunctionToImageFilterNew.cxx
otbFunctionToImageFilter.cxx
otbScalarImageTextureFunctorNew.cxx
diff --git a/Testing/Code/BasicFilters/otbBasicFiltersTests9.cxx b/Testing/Code/BasicFilters/otbBasicFiltersTests9.cxx
index e752c5f9e906a23af98b5bb351e2eb5867e3b90e..b68119fcf2898787f1ff9f26598de30eae01f9fa 100644
--- a/Testing/Code/BasicFilters/otbBasicFiltersTests9.cxx
+++ b/Testing/Code/BasicFilters/otbBasicFiltersTests9.cxx
@@ -31,8 +31,8 @@ REGISTER_TEST(otbChangeLabelImageFilterNew);
REGISTER_TEST(otbChangeLabelImageFilterTest);
REGISTER_TEST(otbContinuousMinimumMaximumImageCalculatorNew);
REGISTER_TEST(otbContinuousMinimumMaximumImageCalculatorTest);
-//REGISTER_TEST(otbMeanShiftImageFilterNew);
-//REGISTER_TEST(otbMeanShiftImageFilter);
+REGISTER_TEST(otbMeanShiftVectorImageFilterNew);
+REGISTER_TEST(otbMeanShiftVectorImageFilter);
//REGISTER_TEST(MeanShiftRef);
//REGISTER_TEST(MeanShiftVectorImageRef);
REGISTER_TEST(otbFunctionToImageFilterNew);
diff --git a/Testing/Code/BasicFilters/otbMeanShiftVectorImageFilter.cxx b/Testing/Code/BasicFilters/otbMeanShiftVectorImageFilter.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..14ac41a5c6aa9361db7fedaaec4151d2648a3c7f
--- /dev/null
+++ b/Testing/Code/BasicFilters/otbMeanShiftVectorImageFilter.cxx
@@ -0,0 +1,67 @@
+/*=========================================================================
+
+ 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 "otbImage.h"
+#include "otbVectorImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+#include "otbMeanShiftVectorImageFilter.h"
+
+int otbMeanShiftVectorImageFilter(int argc, char * argv[])
+{
+ if(argc != 8)
+ {
+ std::cerr<<"Usage: "<<argv[0]<<" infname outfname spatialRadius rangeRadius maxNbIterations useImageSpacing convergenceDistanceThreshold"<<std::endl;
+ return EXIT_FAILURE;
+ }
+
+ const char * infname = argv[1];
+ const char * outfname = argv[2];
+ const double spatialRadius = atof(argv[3]);
+ const double rangeRadius = atof(argv[4]);
+ const unsigned int maxNbIterations = atoi(argv[5]);
+ const bool useImageSpacing = atoi(argv[6]);
+ const double convergenceTol = atof(argv[7]);
+
+ const unsigned int Dimension = 2;
+ typedef double PixelType;
+ typedef otb::VectorImage<PixelType,Dimension> ImageType;
+ typedef otb::ImageFileReader<ImageType> ReaderType;
+ typedef otb::ImageFileWriter<ImageType> WriterType;
+ typedef otb::MeanShiftVectorImageFilter<ImageType,ImageType> FilterType;
+
+ // Instantiating object
+ FilterType::Pointer filter = FilterType::New();
+ ReaderType::Pointer reader = ReaderType::New();
+ WriterType::Pointer writer = WriterType::New();
+
+ reader->SetFileName(infname);
+ writer->SetFileName(outfname);
+
+ filter->SetSpatialRadius(spatialRadius);
+ filter->SetRangeRadius(rangeRadius);
+ filter->SetMaxNumberOfIterations(maxNbIterations);
+ filter->SetUseImageSpacing(useImageSpacing);
+ filter->SetConvergenceDistanceThreshold(convergenceTol);
+
+ filter->SetInput(reader->GetOutput());
+ writer->SetInput(filter->GetOutput());
+
+ writer->Update();
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/BasicFilters/otbMeanShiftVectorImageFilterNew.cxx b/Testing/Code/BasicFilters/otbMeanShiftVectorImageFilterNew.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..8cb9aa85d711640350943b6dad40ed60fa4e7648
--- /dev/null
+++ b/Testing/Code/BasicFilters/otbMeanShiftVectorImageFilterNew.cxx
@@ -0,0 +1,32 @@
+/*=========================================================================
+
+ 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 "otbVectorImage.h"
+#include "otbMeanShiftVectorImageFilter.h"
+
+int otbMeanShiftVectorImageFilterNew(int argc, char * argv[])
+{
+ const unsigned int Dimension = 2;
+ typedef double PixelType;
+ typedef otb::VectorImage<PixelType,Dimension> ImageType;
+ typedef otb::MeanShiftVectorImageFilter<ImageType,ImageType> FilterType;
+
+ // Instantiating object
+ FilterType::Pointer object = FilterType::New();
+ return EXIT_SUCCESS;
+}