diff --git a/Code/BasicFilters/otbBinaryImageToDensityImageFilter.txx b/Code/BasicFilters/otbBinaryImageToDensityImageFilter.txx
index 76cb3fc4b1ae8a34f2dc9f4a3ffe4baa57868ab4..2a95dd2579a77402a9dc662c3eb61f11f77598cf 100644
--- a/Code/BasicFilters/otbBinaryImageToDensityImageFilter.txx
+++ b/Code/BasicFilters/otbBinaryImageToDensityImageFilter.txx
@@ -20,6 +20,7 @@
 
 #include "otbBinaryImageToDensityImageFilter.h"
 #include "itkImageRegionIterator.h"
+#include "itkImageRegionConstIterator.h"
 #include "otbMacro.h"
 
 namespace otb
@@ -53,7 +54,7 @@ BinaryImageToDensityImageFilter<TInputImage, TOutputImage, TCountFunction>
 
   m_CountFunction->SetInputImage(input);
 
-  itk::ImageRegionIterator<InputImageType>   it(input,outputRegionForThread );
+  itk::ImageRegionConstIterator<InputImageType>   it(input,outputRegionForThread );
   itk::ImageRegionIterator<OutputImageType>   itOut(output,outputRegionForThread );
   
   it.GoToBegin();
diff --git a/Code/BasicFilters/otbEdgeDensityImageFilter.h b/Code/BasicFilters/otbEdgeDensityImageFilter.h
index 7abdd35848fff1204ec0ff9fa006c3e40cd1b969..e2d6f6ee270a3aaa6a09bcea5245c6c9dfe65c5a 100644
--- a/Code/BasicFilters/otbEdgeDensityImageFilter.h
+++ b/Code/BasicFilters/otbEdgeDensityImageFilter.h
@@ -36,7 +36,7 @@ PURPOSE.  See the above copyright notices for more information.
 
 namespace otb
 {
-template <class TInputImage , class TEdgeDetector, class TDensityCount , class TOutputImage>
+template <class TInputImage, class TOutputImage , class TEdgeDetector, class TDensityCount>
 class ITK_EXPORT EdgeDensityImageFilter
       : public itk::ImageToImageFilter<TInputImage, TOutputImage>
 {
diff --git a/Code/BasicFilters/otbEdgeDensityImageFilter.txx b/Code/BasicFilters/otbEdgeDensityImageFilter.txx
index bb169c19c61e3e3bc01dc4d747ca812352397c32..806cad874a311c2840927e50c9df4624878b866f 100644
--- a/Code/BasicFilters/otbEdgeDensityImageFilter.txx
+++ b/Code/BasicFilters/otbEdgeDensityImageFilter.txx
@@ -26,8 +26,8 @@ namespace otb
 /**---------------------------------------------------------
  * Constructor
  ----------------------------------------------------------*/
-template <class TInputImage , class TEdgeDetector, class TDensityCount , class TOutputImage>
-EdgeDensityImageFilter<TInputImage, TEdgeDetector, TDensityCount, TOutputImage>
+template <class TInputImage , class TOutputImage, class TEdgeDetector, class TDensityCount>
+EdgeDensityImageFilter<TInputImage, TOutputImage, TEdgeDetector, TDensityCount>
 ::EdgeDensityImageFilter()
 {
   this->SetNumberOfRequiredInputs( 1 );
@@ -41,8 +41,8 @@ EdgeDensityImageFilter<TInputImage, TEdgeDetector, TDensityCount, TOutputImage>
 /*---------------------------------------------------------
  * Destructor.c
  ----------------------------------------------------------*/
-template <class TInputImage , class TEdgeDetector, class TDensityCount , class TOutputImage>
-EdgeDensityImageFilter<TInputImage, TEdgeDetector, TDensityCount, TOutputImage>
+template <class TInputImage , class TOutputImage, class TEdgeDetector, class TDensityCount>
+EdgeDensityImageFilter<TInputImage, TOutputImage, TEdgeDetector, TDensityCount>
 ::~EdgeDensityImageFilter()
 {}
 
@@ -53,9 +53,9 @@ EdgeDensityImageFilter<TInputImage, TEdgeDetector, TDensityCount, TOutputImage>
 /**
 * ThreadedGenerateData Performs the pixel-wise addition
 */
-template <class TInputImage , class TEdgeDetector, class TDensityCount , class TOutputImage>
+template <class TInputImage , class TOutputImage, class TEdgeDetector, class TDensityCount>
 void
-EdgeDensityImageFilter<TInputImage, TEdgeDetector, TDensityCount, TOutputImage>
+EdgeDensityImageFilter<TInputImage, TOutputImage, TEdgeDetector, TDensityCount>
 //::GenerateData()
 ::ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, int threadId )
 {
@@ -82,9 +82,9 @@ EdgeDensityImageFilter<TInputImage, TEdgeDetector, TDensityCount, TOutputImage>
 /**
  * Set Detector
  */
-template <class TInputImage , class TEdgeDetector, class TDensityCount , class TOutputImage>
+template <class TInputImage , class TOutputImage, class TEdgeDetector, class TDensityCount>
 void
-EdgeDensityImageFilter<TInputImage, TEdgeDetector, TDensityCount, TOutputImage>
+EdgeDensityImageFilter<TInputImage, TOutputImage, TEdgeDetector, TDensityCount>
 ::SetDetector(DetectorType* detector)
 {
   m_Detector = detector;
@@ -94,10 +94,10 @@ EdgeDensityImageFilter<TInputImage, TEdgeDetector, TDensityCount, TOutputImage>
 /**
  * Get Detector
  */
-template <class TInputImage , class TEdgeDetector, class TDensityCount , class TOutputImage>
-typename EdgeDensityImageFilter<TInputImage, TEdgeDetector, TDensityCount, TOutputImage> 
+template <class TInputImage , class TOutputImage, class TEdgeDetector, class TDensityCount>
+typename EdgeDensityImageFilter<TInputImage, TOutputImage, TEdgeDetector, TDensityCount> 
 ::DetectorType *
-EdgeDensityImageFilter<TInputImage, TEdgeDetector, TDensityCount, TOutputImage> 
+EdgeDensityImageFilter<TInputImage, TOutputImage, TEdgeDetector, TDensityCount> 
 ::GetDetector()
 {
   return m_Detector;
@@ -108,9 +108,9 @@ EdgeDensityImageFilter<TInputImage, TEdgeDetector, TDensityCount, TOutputImage>
 /*----------------------------------------------------------------
   PrintSelf
   -----------------------------------------------------------------*/
-template <class TInputImage , class TEdgeDetector, class TDensityCount , class TOutputImage>
+template <class TInputImage , class TOutputImage, class TEdgeDetector, class TDensityCount>
 void
-EdgeDensityImageFilter<TInputImage, TEdgeDetector, TDensityCount, TOutputImage>  
+EdgeDensityImageFilter<TInputImage, TOutputImage, TEdgeDetector, TDensityCount>  
 ::PrintSelf(std::ostream& os, itk::Indent indent) const
 {
   Superclass::PrintSelf(os, indent);
diff --git a/Code/BasicFilters/otbKeyPointDensityImageFilter.h b/Code/BasicFilters/otbKeyPointDensityImageFilter.h
index efac0877f37dc8e01c5f80b77989117c4baf7da1..c306fb6fb4536bc7b6d3323f1ee610fc4788a5b8 100644
--- a/Code/BasicFilters/otbKeyPointDensityImageFilter.h
+++ b/Code/BasicFilters/otbKeyPointDensityImageFilter.h
@@ -35,7 +35,7 @@ PURPOSE.  See the above copyright notices for more information.
 
 namespace otb
 {
-template <class TInputImage , class TDetector, class TOutputImage>
+template <class TInputImage, class TOutputImage , class TDetector>
 class ITK_EXPORT KeyPointDensityImageFilter
       : public itk::ImageToImageFilter<TInputImage, TOutputImage>
 {
diff --git a/Code/BasicFilters/otbKeyPointDensityImageFilter.txx b/Code/BasicFilters/otbKeyPointDensityImageFilter.txx
index c16649ad4c21d49e312f2dc44548c296cce65911..7b9d1bd056b0b9ae3500f5116be63848d08382c7 100644
--- a/Code/BasicFilters/otbKeyPointDensityImageFilter.txx
+++ b/Code/BasicFilters/otbKeyPointDensityImageFilter.txx
@@ -27,8 +27,8 @@ namespace otb
 /**---------------------------------------------------------
  * Constructor
  ----------------------------------------------------------*/
-template <class TInputImage , class TDetector, class TOutputImage>
-KeyPointDensityImageFilter<TInputImage , TDetector,  TOutputImage>
+template <class TInputImage , class TOutputImage, class TDetector>
+KeyPointDensityImageFilter<TInputImage , TOutputImage, TDetector>
 ::KeyPointDensityImageFilter()
 {
   this->SetNumberOfRequiredInputs( 1 );
@@ -41,8 +41,8 @@ KeyPointDensityImageFilter<TInputImage , TDetector,  TOutputImage>
 /*---------------------------------------------------------
  * Destructor.c
  ----------------------------------------------------------*/
-template <class TInputImage , class TDetector, class TOutputImage>
-KeyPointDensityImageFilter<TInputImage, TDetector, TOutputImage >
+template <class TInputImage , class TOutputImage, class TDetector>
+KeyPointDensityImageFilter<TInputImage, TOutputImage, TDetector >
 ::~KeyPointDensityImageFilter()
 {}
 
@@ -53,9 +53,9 @@ KeyPointDensityImageFilter<TInputImage, TDetector, TOutputImage >
 /**
 * ThreadedGenerateData Performs the pixel-wise addition
 */
-template <class TInputImage , class TDetector,  class TOutputImage>
+template <class TInputImage , class TOutputImage, class TDetector>
 void
-KeyPointDensityImageFilter<TInputImage, TDetector, TOutputImage >
+KeyPointDensityImageFilter<TInputImage, TOutputImage, TDetector >
 ::GenerateData()
 //::GenerateData( const OutputImageRegionType &outputRegionForThread, int threadId )
 {
@@ -83,9 +83,9 @@ KeyPointDensityImageFilter<TInputImage, TDetector, TOutputImage >
 /**
  * Set Detector
  */
-template <class TInputImage , class TDetector, class TOutputImage>
+template <class TInputImage , class TOutputImage, class TDetector>
 void
-KeyPointDensityImageFilter<TInputImage, TDetector, TOutputImage >
+KeyPointDensityImageFilter<TInputImage, TOutputImage, TDetector >
 ::SetDetector(DetectorType* detector)
 {
   m_Detector = detector;
@@ -95,10 +95,10 @@ KeyPointDensityImageFilter<TInputImage, TDetector, TOutputImage >
 /**
  * Get Detector
  */
-template <class TInputImage , class TDetector, class TOutputImage>
-typename KeyPointDensityImageFilter< TInputImage ,  TDetector,  TOutputImage >
+template <class TInputImage , class TOutputImage, class TDetector>
+typename KeyPointDensityImageFilter< TInputImage ,  TOutputImage, TDetector >
 ::DetectorType *
-KeyPointDensityImageFilter< TInputImage , TDetector, TOutputImage >
+KeyPointDensityImageFilter< TInputImage , TOutputImage, TDetector >
 ::GetDetector()
 {
   return m_Detector;
@@ -109,9 +109,9 @@ KeyPointDensityImageFilter< TInputImage , TDetector, TOutputImage >
 /*----------------------------------------------------------------
   PrintSelf
   -----------------------------------------------------------------*/
-template <class TInputImage , class TDetector, class TOutputImage>
+template <class TInputImage , class TOutputImage, class TDetector>
 void
-KeyPointDensityImageFilter< TInputImage ,  TDetector,  TOutputImage >
+KeyPointDensityImageFilter< TInputImage ,  TOutputImage, TDetector >
 ::PrintSelf(std::ostream& os, itk::Indent indent) const
 {
   Superclass::PrintSelf(os, indent);
diff --git a/Code/BasicFilters/otbProlateInterpolateImageFunction.h b/Code/BasicFilters/otbProlateInterpolateImageFunction.h
index b3435b5691b4c32186964ab8b2d3722889069c17..686c422c6b2c5adcdee9cb2a151f230d2f94a444 100644
--- a/Code/BasicFilters/otbProlateInterpolateImageFunction.h
+++ b/Code/BasicFilters/otbProlateInterpolateImageFunction.h
@@ -62,7 +62,7 @@ public:
   {
     return m_OriginalProfile;
   }
-  double ComputeEnergy(double resampleRatio);
+  double ComputeEnergy(double resampleRatio) const;
 
   inline TOutput operator()( const TInput & A ) const
   {
diff --git a/Code/BasicFilters/otbProlateInterpolateImageFunction.txx b/Code/BasicFilters/otbProlateInterpolateImageFunction.txx
index 19b0f0eca8aaf46611ebad35aa477af59dea5936..86b95ff465d46ad68d7875bb1b10b3089fe8cbcb 100644
--- a/Code/BasicFilters/otbProlateInterpolateImageFunction.txx
+++ b/Code/BasicFilters/otbProlateInterpolateImageFunction.txx
@@ -101,7 +101,7 @@ ProlateFunction<TInput, TOutput>
 template<class TInput, class TOutput>
 double
 ProlateFunction<TInput, TOutput>
-::ComputeEnergy(double resampleRatio)
+::ComputeEnergy(double resampleRatio) const
 {
   vnl_vector<vcl_complex<double> > resampledProfile(1024);
   resampledProfile.fill(0);
diff --git a/Code/Common/otbDrawLineSpatialObjectListFilter.h b/Code/Common/otbDrawLineSpatialObjectListFilter.h
index a186bc8f370e67895d66f33eb0557f2dba407a0e..5ada193180a4917557b2d59c236d0fa61a6fc1b2 100644
--- a/Code/Common/otbDrawLineSpatialObjectListFilter.h
+++ b/Code/Common/otbDrawLineSpatialObjectListFilter.h
@@ -94,12 +94,31 @@ protected:
   void PrintSelf(std::ostream& os, itk::Indent indent) const;
 
   virtual void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, int threadId ) ;
+  
+  /**
+   * compute the intersection of the segment to draw with the region
+   */
+  virtual void CropSegment(OutputIndexType *indexToCrop,OutputIndexType *otherIndex, const OutputImageRegionType *outputRegionForThread )const;
+  
+  virtual void CropRightSegment(OutputIndexType *indexToCrop,OutputIndexType *otherIndex, const OutputImageRegionType *outputRegionForThread )const;
+  
+
+  virtual bool IsUpsideTheRegion(OutputIndexType *indexToCrop, const OutputImageRegionType *outputRegionForThread ) const;
+
+  virtual bool IsDownsideTheRegion(OutputIndexType *indexToCrop, const OutputImageRegionType *outputRegionForThread )const;
+  
+  virtual bool IsDownsideTheImage(OutputIndexType *indexToCrop )const;
+    
+  
+  virtual bool IsColumnOutsideOfTheRegion(OutputIndexType *indexToCheck, OutputIndexType *otherToCheck, const OutputImageRegionType *outputRegionForThread)const;
 
 private:
   DrawLineSpatialObjectListFilter(const Self&); //purposely not implemented
   void operator=(const Self&); //purposely not implemented
 
   OutputPixelType       m_Value;
+  int          m_Length;
+  int          m_Width;
 
 };
 } // end namespace otb
diff --git a/Code/Common/otbDrawLineSpatialObjectListFilter.txx b/Code/Common/otbDrawLineSpatialObjectListFilter.txx
index 305de1436a49e1cd877b76a83f14827de80a8398..f171477f71f7b932a360a5c3c474df1465694b36 100644
--- a/Code/Common/otbDrawLineSpatialObjectListFilter.txx
+++ b/Code/Common/otbDrawLineSpatialObjectListFilter.txx
@@ -59,7 +59,7 @@ DrawLineSpatialObjectListFilter<TInputImage, TOutputImage>
 ::GetInputLineSpatialObjectList(void)
 {
   //ROMAIN
-  return static_cast</*const */LinesListType *>
+  return static_cast</*const*/ LinesListType* >
          (this->ProcessObjectType::GetInput(1) );
 }
 
@@ -71,7 +71,7 @@ DrawLineSpatialObjectListFilter<TInputImage, TOutputImage>
 {
   typename InputImageType::ConstPointer input  = this->GetInput();
   typename OutputImageType::Pointer output  = this->GetOutput();
-  typename LinesListType::Pointer    list = this->GetInputLineSpatialObjectList();
+  typename LinesListType::Pointer    list = const_cast<LinesListType*>(this->GetInputLineSpatialObjectList());
        
   /** Copy the input requested region in the output requested region*/
   typedef itk::ImageRegionIterator< OutputImageType >      OutputIteratorType;
@@ -91,9 +91,11 @@ DrawLineSpatialObjectListFilter<TInputImage, TOutputImage>
   OutputIndexType                                  indexBeginLine, indexEndLine;
   LineListIterator                                 itList = list->begin();
 
-  /** Dimensions of the buffered region*/
-  //   typename OutputImageRegionType::SizeType size = outputRegionForThread.GetSize();
-  //   typename OutputImageRegionType::IndexType start = outputRegionForThread.GetIndex();
+  typename InputImageType::SizeType size = input->GetLargestPossibleRegion().GetSize();
+  m_Length = size[1];
+  m_Width  = size[0];
+  
+  
   
   while(itList != list->end())
     {
@@ -107,23 +109,216 @@ DrawLineSpatialObjectListFilter<TInputImage, TOutputImage>
 
       indexEndLine[0] = static_cast<unsigned int>((*itPoints).GetPosition()[0]);
       indexEndLine[1] = static_cast<unsigned int>((*itPoints).GetPosition()[1]);
+      
 
-      /** Instanciation of the line iterator with begin and ending index*/
-      LineIteratorFilter   itLine(output,indexBeginLine ,indexEndLine  );
+      /** Crop the segment if it is outside the region in the left*/
       
-      /** Iteration over the line and writing white lines */
-      while(!itLine.IsAtEnd())
+      if( !(this->IsColumnOutsideOfTheRegion(&indexBeginLine,&indexEndLine,&outputRegionForThread) && this->IsColumnOutsideOfTheRegion(&indexEndLine,&indexBeginLine,&outputRegionForThread)))
 	{
-	if(outputRegionForThread.IsInside(itLine.GetIndex()))
-	  itLine.Set(m_Value);
-	++itLine;
+	  if(indexEndLine[0] >=static_cast< int>(size[0]))
+	    this->CropRightSegment(&indexEndLine,&indexBeginLine, &outputRegionForThread);
+	  
+	  if( indexBeginLine[0] >= static_cast< int>(size[0]) )
+	    this->CropRightSegment(&indexBeginLine,&indexEndLine, &outputRegionForThread);
 	}
+
+      /** 
+       * If an extremity is under the region 
+       * Technically, the X component of the index is inside the image 
+       */
+      if(this->IsDownsideTheImage(&indexBeginLine)  && input->GetLargestPossibleRegion().IsInside(indexEndLine))
+	this->CropSegment(&indexBeginLine,&indexEndLine, &outputRegionForThread);
+      
+      if(this->IsDownsideTheImage(&indexEndLine) && input->GetLargestPossibleRegion().IsInside(indexBeginLine))
+	this->CropSegment(&indexEndLine,&indexBeginLine, &outputRegionForThread);
+    
+
+      /** If the segments are not in the region (upside or downside the region)*/
+      if(!(this->IsUpsideTheRegion(&indexBeginLine,&outputRegionForThread)   && this->IsUpsideTheRegion(&indexEndLine,&outputRegionForThread)) &&
+	 !(this->IsDownsideTheRegion(&indexBeginLine,&outputRegionForThread) && this->IsDownsideTheRegion(&indexEndLine,&outputRegionForThread)) &&
+	 !(this->IsColumnOutsideOfTheRegion(&indexBeginLine,&indexEndLine, &outputRegionForThread) && this->IsColumnOutsideOfTheRegion(&indexEndLine,&indexBeginLine,&outputRegionForThread)) 
+	 )
+	{
+
+	  /** Instanciation of the line iterator with begin and ending index*/
+	  LineIteratorFilter   itLine(output,indexBeginLine ,indexEndLine  );
       
+	  /** Iteration over the line and writing white lines */
+	  while(!itLine.IsAtEnd())
+	    {
+	      if(outputRegionForThread.IsInside(itLine.GetIndex()))
+		itLine.Set(m_Value);
+	      ++itLine;
+	    }
+	}
+   
       ++itList;
     }
 }
 
 
+/**
+ * Compute the intersection between the segment to draw and the region belonging to the current thread
+ * It is used if the segment abcisse is greater than the width of the image
+ *
+ */
+template <class TInputImage, class TOutput>
+void
+DrawLineSpatialObjectListFilter<TInputImage, TOutput>
+::CropRightSegment(OutputIndexType *indexToCrop,OutputIndexType *otherIndex, const OutputImageRegionType *outputRegionForThread ) const
+
+{
+
+ /** Dimensions of the buffered region*/
+ typename OutputImageRegionType::SizeType  size  = outputRegionForThread->GetSize();
+ typename OutputImageRegionType::IndexType start = outputRegionForThread->GetIndex();
+ 
+ /** Equation of the line (Begin, End)*/
+ double lengthSegment = -(*otherIndex)[1] + (*indexToCrop)[1];
+ double slope         =  lengthSegment/(  (*indexToCrop)[0]  - (*otherIndex)[0]);
+ double origin        =  (*otherIndex)[1] - (slope * (*otherIndex)[0]);
+
+ (*indexToCrop)[0] = static_cast<unsigned int>(size[0]-1);
+ (*indexToCrop)[1] = static_cast<unsigned int>(slope *(*indexToCrop)[0] + origin +0.5) ;
+}
+
+/**
+ * Define if the point is upside the region or not
+ *
+ */
+template <class TInputImage, class TOutput>
+bool
+DrawLineSpatialObjectListFilter<TInputImage, TOutput>
+::IsUpsideTheRegion(OutputIndexType *indexToCrop, const OutputImageRegionType *outputRegionForThread ) const
+
+{
+  /** Dimensions of the buffered region*/
+  typename OutputImageRegionType::SizeType  size  = outputRegionForThread->GetSize();
+  typename OutputImageRegionType::IndexType start = outputRegionForThread->GetIndex();
+
+  return (*indexToCrop)[1] < static_cast<unsigned int>(start[1]);
+}
+
+/**
+ * Define if the point is upside the region or not
+ *
+ */
+template <class TInputImage, class TOutput>
+bool
+DrawLineSpatialObjectListFilter<TInputImage, TOutput>
+::IsDownsideTheRegion(OutputIndexType *indexToCrop, const OutputImageRegionType *outputRegionForThread ) const
+
+{
+  /** Dimensions of the buffered region*/
+  typename OutputImageRegionType::SizeType  size  = outputRegionForThread->GetSize();
+  typename OutputImageRegionType::IndexType start = outputRegionForThread->GetIndex();
+
+  return (*indexToCrop)[1] >= static_cast< int>(start[1]+size[1]); //The down limit of the region in the Y direction
+}
+
+/**
+ * Define if the point is outside the hole image
+ *
+ */
+template <class TInputImage, class TOutput>
+bool
+DrawLineSpatialObjectListFilter<TInputImage, TOutput>
+::IsDownsideTheImage(OutputIndexType *indexToCrop) const
+
+{
+  return (*indexToCrop)[1] >= static_cast<int>(m_Length); //The down limit of the Image in the Y direction
+}
+
+
+/**
+ * 
+ *
+ */
+template <class TInputImage, class TOutput>
+bool
+DrawLineSpatialObjectListFilter<TInputImage, TOutput>
+::IsColumnOutsideOfTheRegion(OutputIndexType *indexToCheck, OutputIndexType *otherToCheck, const OutputImageRegionType *outputRegionForThread) const
+{
+  /** Dimensions of the buffered region*/
+  typename OutputImageRegionType::SizeType size = outputRegionForThread->GetSize();
+  bool res = false, res1= false , res2 = false;
+  
+  if (  ((*indexToCheck)[0]>=static_cast< int>(size[0])) && ((*otherToCheck)[0]>=static_cast< int>(size[0]) ))
+    res  = true;
+  
+  if((*indexToCheck)[0]>=static_cast< int>(size[0]) && this->IsUpsideTheRegion(otherToCheck,outputRegionForThread))  
+    res1 = true;
+  
+  if((*indexToCheck)[0]>=static_cast< int>(size[0]) && this->IsDownsideTheRegion(otherToCheck,outputRegionForThread) )
+    res2 = true;
+
+
+  
+  return res || res1 || res2;
+
+}
+
+
+/**
+ * Compute the intersection between the segment to draw and the region belonging to the current thread
+ *
+ */
+template <class TInputImage, class TOutput>
+void
+DrawLineSpatialObjectListFilter<TInputImage, TOutput>
+::CropSegment(OutputIndexType *indexToCrop,OutputIndexType *otherIndex,  const OutputImageRegionType *outputRegionForThread) const
+
+{
+  OutputIndexType tempIndex;
+  tempIndex = *indexToCrop;
+  
+  /** Dimensions of the buffered region*/
+  typename OutputImageRegionType::SizeType  size = outputRegionForThread->GetSize();
+  typename OutputImageRegionType::IndexType start = outputRegionForThread->GetIndex();
+
+  /** Equation of the line (Begin, End)*/
+  double slope = 0.;
+  double lengthSegment =0.;
+  double origin = 0.;
+  double tempOtherIndexX= 0.;
+  
+  /** Equation of the first Line*/
+
+  if(vcl_abs( (*otherIndex)[0] -(*indexToCrop)[0] ) <1e-4)
+    tempOtherIndexX= 0.000001;
+  else
+    tempOtherIndexX = static_cast<double>((*otherIndex)[0]);
+
+
+  if( (*indexToCrop)[0] < (*otherIndex)[0])
+    lengthSegment = (*otherIndex)[1] -(*indexToCrop)[1];
+  else 
+    lengthSegment = (*indexToCrop)[1] -(*otherIndex)[1];
+  
+  slope = lengthSegment/( tempOtherIndexX - static_cast<double>((*indexToCrop)[0]));
+  origin = (*indexToCrop)[1] - (slope * static_cast<double>((*indexToCrop)[0]));
+  
+
+  if((*indexToCrop)[1] < 0)
+    {
+      unsigned int Y = 0;
+      tempIndex[1] = Y;
+      tempIndex[0] = static_cast<unsigned int>((Y-origin) / slope);  // X = (Y-B)/A
+    }
+  
+  if(this->IsDownsideTheImage(indexToCrop ))
+    {
+      double Y = static_cast<double>(m_Length-1)/*tstart[1]+size[1]-1*/;
+      tempIndex[1] = static_cast<unsigned int>(Y);
+      tempIndex[0] = static_cast<unsigned int>((Y-origin) / slope);  // X = (Y-B)/A
+    } 
+
+  (*indexToCrop)[0] = tempIndex[0];
+  (*indexToCrop)[1] = tempIndex[1];
+}
+
+
+
 /**
  * Standard "PrintSelf" method
  */
diff --git a/Code/Common/otbObjectList.h b/Code/Common/otbObjectList.h
index 66366d4d0a642b92e4e8fe8c7bd98794df6700e7..1d8f4a56727716c83d6e21075afa55fc0080dac8 100644
--- a/Code/Common/otbObjectList.h
+++ b/Code/Common/otbObjectList.h
@@ -117,6 +117,20 @@ public:
   friend class ReverseIterator;
   friend class ReverseConstIterator;
 
+  /**
+   * Insert an element at a given position 
+   * \param position A random access iterator
+   * \return An iterator that points to the newly insereted element.
+   */
+  Iterator Insert ( Iterator position, ObjectPointerType element );
+  /**
+   * Insert an element at a given position 
+   * \param position A reverse iterator
+   * \return A reverse iterator that points to the newly insereted element.
+   */
+  ReverseIterator Insert ( ReverseIterator position, ObjectPointerType element );
+
+
   /** \class Iterator
    *  \brief Iterator of the object list.
    */
@@ -142,6 +156,13 @@ public:
     {
       return (*m_Iter);
     };
+    /**
+     * Set the current object
+     */
+    void Set ( ObjectPointerType element )
+    {
+      (*m_Iter) = element;
+    }
     /**
        * Increment.
        */
@@ -333,6 +354,13 @@ public:
     {
       return (*m_Iter);
     };
+    /**
+     * Set the current object
+     */
+    void Set ( ObjectPointerType element )
+    {
+      (*m_Iter) = element;
+    }
     /**
        * Increment.
        */
@@ -379,6 +407,14 @@ public:
       m_Iter=it.m_Iter;
     };
 
+    /**
+      * Get the current internal iterator
+    */
+    InternalReverseIteratorType & GetIter(void)
+    {
+      return(m_Iter);
+    }
+
   private:
     // Internal iterator.
     InternalReverseIteratorType m_Iter;
diff --git a/Code/Common/otbObjectList.txx b/Code/Common/otbObjectList.txx
index 3d337bcd8c3aa2622daf44e722b0578e973202cd..45d27e7364c87079b6ae198f26a45a38d256af94 100644
--- a/Code/Common/otbObjectList.txx
+++ b/Code/Common/otbObjectList.txx
@@ -166,7 +166,29 @@ ObjectList<TObject>
   m_InternalContainer.clear();
 }
 
+/**
+ * Insert an element
+ */
+template <class TObject>
+typename ObjectList<TObject>::Iterator
+ObjectList<TObject>
+::Insert(Iterator position, ObjectPointerType element)
+{
+  Iterator iter ( m_InternalContainer.insert( position.GetIter(), element ) );
+  return iter;
+}
 
+/**
+ * Insert an element
+ */
+template <class TObject>
+typename ObjectList<TObject>::ReverseIterator
+ObjectList<TObject>
+::Insert(ReverseIterator position, ObjectPointerType element)
+{
+  ReverseIterator iter ( m_InternalContainer.insert( position.GetIter(), element ) );
+  return iter;
+}
 
 /**
  * Get an Iterator that points to the beginning of the container.
diff --git a/Code/Common/otbTestMain.h b/Code/Common/otbTestMain.h
index 530d64baecff5bd619bbff565be1e5855fbbd42b..f0fc741078015707bd79a808186f106c4f5fd7b8 100644
--- a/Code/Common/otbTestMain.h
+++ b/Code/Common/otbTestMain.h
@@ -668,11 +668,11 @@ bool isHexaPointerAddress(std::string str)
   return result;
 }
 
-std::string VectorToString(otb::ImageBase::VectorType vector)
+std::string VectorToString(otb::MetaDataKey::VectorType vector)
 {
   otb::StringStream oss;
   oss.str("");
-  otb::ImageBase::VectorType::iterator it = vector.begin();
+  otb::MetaDataKey::VectorType::iterator it = vector.begin();
   oss<<"[";
   while (it!=vector.end())
   {
diff --git a/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.txx b/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.txx
index 263567b5978fd0c8b228095821ffeb1115f23ea2..ac9a423ed5dfae397aef09156bb5aaa6daf055b1 100644
--- a/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.txx
+++ b/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.txx
@@ -159,7 +159,7 @@ UnaryFunctorNeighborhoodWithOffsetImageFilter<TInputImage, TOutputImage, TFuncti
     while ( ! outputIt.IsAtEnd() )
     {
 
-      outputIt.Set( m_FunctorList[threadId]( neighInputOffIt ) );
+      outputIt.Set( m_FunctorList[threadId]( neighInputOffIt.GetNeighborhood() ) );
 
       ++neighInputOffIt;
       ++outputIt;
diff --git a/Code/Common/otbVectorDataExtractROI.h b/Code/Common/otbVectorDataExtractROI.h
new file mode 100644
index 0000000000000000000000000000000000000000..bbeeef948bdacf87997e460132c4ce0b65376500
--- /dev/null
+++ b/Code/Common/otbVectorDataExtractROI.h
@@ -0,0 +1,106 @@
+/*=========================================================================
+
+  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 __otbVectorDataExtractROI_h
+#define __otbVectorDataExtractROI_h
+
+#include "otbVectorDataToVectorDataFilter.h"
+#include "itkMacro.h"
+
+namespace otb
+{
+
+/** \class VectorDataExtractROI
+ * \brief Extract a subset of a Vector Data.
+ *
+ * \note Parameter to this class for input and outputs are vectorData
+ *
+ */
+template <class TVectorData>
+class ITK_EXPORT VectorDataExtractROI:
+      public VectorDataToVectorDataFilter<TVectorData ,TVectorData >
+{
+public:
+  /** Standard class typedefs. */
+  typedef VectorDataExtractROI                                     Self;
+  typedef VectorDataToVectorDataFilter<TVectorData ,TVectorData >  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(VectorDataExtractROI,VectorDataToVectorDataFilter);
+
+  /** Image type information. */
+  typedef TVectorData                                        VectorDataType;
+  typedef typename VectorDataType::DataTreeType               DataTreeType;
+  
+  /** Get/Set Macro for ROI Column size */
+  itkGetMacro(SizeX, unsigned int);
+  itkSetMacro(SizeX, unsigned int);
+  
+  /** Get/Set Macro for ROI Lines size */
+  itkGetMacro(SizeY, unsigned int);
+  itkSetMacro(SizeY, unsigned int);
+  
+  /** Get/Set Macro for ROI Start Point Coordinate */
+  itkGetMacro(StartX, unsigned int);
+  itkSetMacro(StartX, unsigned int);
+  
+  /** Get/Set Macro for  ROI Start Point Abcisse */
+  itkGetMacro(StartY, unsigned int);
+  itkSetMacro(StartY, unsigned int);  
+
+  
+  /** Prototype of the generate data method*/
+  void GenerateData(void );
+
+  
+  
+
+protected:
+  VectorDataExtractROI();
+  ~VectorDataExtractROI() {};
+  void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+  /** VectorDataExtractROI
+   *
+   * \sa VectorDataExtractROIBase::GenerateOutputInformaton()  */
+  //virtual void GenerateOutputInformation();
+
+private:
+  VectorDataExtractROI(const Self&); //purposely not implemented
+  void operator=(const Self&); //purposely not implemented
+
+  unsigned int m_SizeX;
+  unsigned int m_SizeY;
+  unsigned int m_StartX;
+  unsigned int m_StartY;
+  
+
+};
+
+
+} // end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbVectorDataExtractROI.txx"
+#endif
+
+#endif
diff --git a/Code/Common/otbVectorDataExtractROI.txx b/Code/Common/otbVectorDataExtractROI.txx
new file mode 100644
index 0000000000000000000000000000000000000000..f545a6c08dc41dbaae1cbbf2a95b11cb0f047b1a
--- /dev/null
+++ b/Code/Common/otbVectorDataExtractROI.txx
@@ -0,0 +1,136 @@
+/*=========================================================================
+
+  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 __otbVectorDataExtractROI_txx
+#define __otbVectorDataExtractROI_txx
+
+#include "otbVectorDataExtractROI.h"
+#include "itkImageRegion.h"
+#include "itkPreOrderTreeIterator.h"
+#include "itkIndex.h"
+#include "itkSize.h"
+
+namespace otb
+{
+
+/**
+ *
+ */
+template <class TVectorData>
+VectorDataExtractROI<TVectorData>
+::VectorDataExtractROI()
+{
+  /** Global Variables initialisation*/
+  m_SizeX = 1;
+  m_SizeY = 1;
+  m_StartX = 0;
+  m_StartY = 0;
+
+  
+  
+}
+
+/**
+ *
+ */
+template <class TVectorData>
+void
+VectorDataExtractROI<TVectorData>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+  Superclass::PrintSelf(os,indent);
+}
+
+
+/**
+ * VectorDataExtractROI can produce an output vector data which is different depth
+ * than the input vectorData
+ * As such, VectorDataExtractROI have to provide an implementation for GenerateOutputInformation
+ * in order to inform the pipeline execution model. 
+ * \sa ProcessObject::GenerateOutputInformaton()
+ */
+// template <class TVectorData>
+// void
+// VectorDataExtractROI<TVectorData>
+// ::GenerateOutputInformation()
+// {
+//   // Call to the base class method
+//   Superclass::GenerateOutputInformation();
+// }
+
+template <class TVectorData>
+void
+VectorDataExtractROI<TVectorData>
+::GenerateData(void)
+{
+  /** Get The input and the outptut*/
+  typename VectorDataType::ConstPointer   input = this->GetInput();
+  if(!input)
+    std::cout << " Probleme avec la recuperation du input"<<std::endl;
+  
+  /** Create a region with the right size*/
+  const   unsigned int                        Dimension = 2;
+  typedef itk::ImageRegion<Dimension>         ImageRegionType;
+  typedef itk::Index<Dimension>               IndexType;
+  typedef itk::Size<Dimension>                SizeType;
+  
+  /** */
+  ImageRegionType                             roi;
+  IndexType                                   index;
+  SizeType                                    size;
+
+  /** Update the region information*/
+  index[0] = m_StartX;
+  index[1] = m_StartY;
+  size[0]  = m_SizeX;
+  size[1]  = m_SizeY;
+
+  /** Create the region*/
+  roi.SetSize(size);
+  roi.SetIndex(index);
+  
+  /** Loop in the vectorData file*/
+  typedef itk::PreOrderTreeIterator<DataTreeType>                 TreeIteratorType;
+  TreeIteratorType                                                it(input->GetDataTree());
+
+  it.GoToBegin();
+  while (!it.IsAtEnd())
+  {
+    
+    itk::PreOrderTreeIterator<DataTreeType> itParent = it;
+    bool goesOn = true;
+
+    if (it.Get()->IsPolygonFeature())
+      std::cout << " C'est un polygone et les coordonnees " << it.Get()->GetPolygonExteriorRing()->GetVertexList()->GetElement(0)  << std::endl;
+    
+    if (it.Get()->IsLineFeature())
+	std::cout << "Vertex List Size " << it.Get()->GetLine()->GetVertexList()->Size() <<std::endl;
+
+    if (it.Get()->IsPointFeature())
+      std::cout << " C'est un point" << std::endl;
+        
+
+    ++it;
+  }
+
+
+}
+
+
+} // end namespace otb
+
+#endif
diff --git a/Code/FeatureExtraction/otbEnergyTextureFunctor.h b/Code/FeatureExtraction/otbEnergyTextureFunctor.h
index 6086c9c5ab874e4f72ddc092097842e817a9f0bb..31467044b469953152e018e16220865c3934cfd6 100755
--- a/Code/FeatureExtraction/otbEnergyTextureFunctor.h
+++ b/Code/FeatureExtraction/otbEnergyTextureFunctor.h
@@ -37,21 +37,21 @@ namespace Functor
  *  \ingroup Functor
  *  \ingroup Statistics
  */
-template <class TIterInput, class TOutput>
+template <class TScalarInputPixelType, class TScalarOutputPixelType>
 class ITK_EXPORT EnergyTextureFunctor :
-public TextureFunctorBase<TIterInput, TOutput>
+public TextureFunctorBase<TScalarInputPixelType, TScalarOutputPixelType>
 {
 public:
   EnergyTextureFunctor(){};
   virtual ~EnergyTextureFunctor(){};
 
-  typedef TIterInput                           IterType;
-  typedef TOutput                              OutputType;
-  typedef typename IterType::OffsetType        OffsetType;
-  typedef typename IterType::RadiusType        RadiusType;
-  typedef typename IterType::InternalPixelType InternalPixelType;
-  typedef typename IterType::ImageType         ImageType;
-  typedef itk::Neighborhood<InternalPixelType,::itk::GetImageDimension<ImageType>::ImageDimension>    NeighborhoodType;
+  typedef TScalarInputPixelType                  InputScalarType;
+  typedef TScalarOutputPixelType                 OutputScalarType;
+  typedef TextureFunctorBase<TScalarInputPixelType, TScalarOutputPixelType> Superclass;
+  typedef typename Superclass::OffsetType        OffsetType;
+  typedef typename Superclass::RadiusType        RadiusType;
+  typedef typename Superclass::NeighborhoodType  NeighborhoodType;
+ 
 
   double ComputeOverSingleChannel(const NeighborhoodType &neigh, const NeighborhoodType &neighOff)
   {
diff --git a/Code/FeatureExtraction/otbEntropyTextureFunctor.h b/Code/FeatureExtraction/otbEntropyTextureFunctor.h
index 1f8175abd654c2ab6eaae15ba7bf5d6cedf69bff..b4163202025f2737194b29b7fc153e7badb910ba 100755
--- a/Code/FeatureExtraction/otbEntropyTextureFunctor.h
+++ b/Code/FeatureExtraction/otbEntropyTextureFunctor.h
@@ -37,20 +37,27 @@ namespace Functor
  *  \ingroup Statistics
  */
 
-template <class TIterInput, class TOutput>
+template <class TScalarInputPixelType, class TScalarOutputPixelType>
 class ITK_EXPORT EntropyTextureFunctor :
-public TextureFunctorBase<TIterInput, TOutput>
+public TextureFunctorBase<TScalarInputPixelType, TScalarOutputPixelType>
 {
 public:
   EntropyTextureFunctor(){};
   virtual ~EntropyTextureFunctor(){};
 
+  typedef TScalarInputPixelType                  InputScalarType;
+  typedef TScalarOutputPixelType                 OutputScalarType;
+  typedef TextureFunctorBase<TScalarInputPixelType, TScalarOutputPixelType> Superclass;
+  typedef typename Superclass::OffsetType        OffsetType;
+  typedef typename Superclass::RadiusType        RadiusType;
+  typedef typename Superclass::NeighborhoodType  NeighborhoodType;
+  /*
   typedef TIterInput                            IterType;
   typedef TOutput                               OutputType;
   typedef typename IterType::InternalPixelType  InternalPixelType;
   typedef typename IterType::ImageType          ImageType;
   typedef itk::Neighborhood<InternalPixelType,::itk::GetImageDimension<ImageType>::ImageDimension>    NeighborhoodType;
-
+  */
 
   virtual double ComputeOverSingleChannel(const NeighborhoodType &neigh, const NeighborhoodType &neighOff)
   {
diff --git a/Code/FeatureExtraction/otbLineDirectionFunctor.h b/Code/FeatureExtraction/otbLineDirectionFunctor.h
index e261f78c9cb42dda0bd7dcdcf485ccdc03bff801..3714ca558825b0d56e24942fe58160d7dcdbcdae 100644
--- a/Code/FeatureExtraction/otbLineDirectionFunctor.h
+++ b/Code/FeatureExtraction/otbLineDirectionFunctor.h
@@ -26,10 +26,28 @@
 namespace otb
 {
 /** \class LineDirectionFunctor
- *  \brief This functor first computes the spectral angle according to a reference pixel.
- *  \brief Then multiplies the result by a gaussian coefficient
- *  \brief And reverse the pixel values.
+ *  \brief This functor computes textures based on line direction analysis through the central pixel.
+ * 
+ *  Directions are computed using NumberOfDirection, used to compute a constant step angle.
+ *  A direction is defined as : $\mathit{d_{i} = \sqrt{(m^{e1}-m{e2})^{2}+(n^{e1}-n{e2})^{2}}}$
+ *  From  $\mathit{d_{i}}, histograms are defined :
+ *  $\mathit{H(c) : \{c \in I \mid \lbrack d_{1}(c), \ldots , d_{i}(c), \ldots , d_{D}(c)\rbrack  \}}$
+ *  Thus, 6 textures are defined :
+ *  $\mathit{length = \max_{i \in \lbrack1;D\rbrack}(d_{i}(c)}$
+ *  $\mathit{width = \min_{i \in \lbrack1;D\rbrack}(d_{i}(c)}$
+ *  $\mathit{PSI = \frac{1}{D}\sum_{1=1}^{D}d_{i}(c)}$
+ *  $\mathit{\omega-mean = \frac{1}{D}\sum_{1=1}^{D}\frac{\alpha.(k_{i}-1)}{st_{i}}d_{i}(c)}$
+ *  $\mathit{ratio = \arctan{\frac{\sum_{j=1}^{n}{sort_{min}^{j}(H(c))}}{\sum_{j=1}^{n}{sort_{max}^{j}(H(c))}}}}$
+ *  $\mathit{SD = \frac{1}{D-1}\sqrt{\sum_{1=1}^{D}(d_{i}(c)-PSI)^{2}}}$
+ *
+ *  For more details, please refer to refer to Xin Huang, Liangpei Zhang and Pingxiang Li publication,
+ *  Classification and Extraction of Spatial Features in Urban Areas
+ *  Using High-Resolution Multispectral Imagery.
+ *  IEEE Geoscience and Remote Sensing Letters,
+ *  vol. 4, n. 2, 2007, pp 260-264
  */
+
+
 namespace Functor
 {
 template<class TIter,class TOutputValue>
@@ -43,15 +61,16 @@ public:
     m_RatioMaxConsiderationNumber = 5;
     m_Alpha = 1;
     this->SetNumberOfDirections(20); // set the step too
-  };
-
+    m_SelectedTextures = std::vector<bool>(6, 1);
+  }
   ~LineDirectionFunctor() {};
 
   typedef typename TIter::InternalPixelType InternalPixelType;
   typedef typename TIter::SizeType          SizeType;
   typedef typename TIter::IndexType         IndexType;
   typedef typename TIter::OffsetType        OffsetType;
-  typedef typename TOutputValue::ValueType  InternalOutputPixelType;
+  typedef TOutputValue                      OutputValueType;
+  typedef std::vector<OutputValueType>      OutputType;
  
   void SetSpatialThreshold( unsigned int thresh ){ m_SpatialThreshold=thresh; };
   void SetSpectralThreshold( InternalPixelType thresh ){ m_SpectralThreshold=thresh; };
@@ -63,39 +82,43 @@ public:
       m_DirectionStep = 2*M_PI/static_cast<double>(D);
     };
   void SetDirectionStep( double step ){ m_DirectionStep = step; };
+  void SetSelectedTextures( std::vector<bool> vect )
+    { 
+      m_SelectedTextures.clear();
+      m_SelectedTextures = vect;
+    };
+  void SetTextureStatus( unsigned int id, bool isSelected ){ m_SelectedTextures[id] = isSelected; };
+
   unsigned int GetSpatialThreshold(){ return m_SpatialThreshold; };
   InternalPixelType GetSpectralThreshold(){ return m_SpectralThreshold; };  
   unsigned int GetRatioMaxConsiderationNumber(){ return m_RatioMaxConsiderationNumber; }; 
   double  GetAlpha(){ return m_Alpha; }; 
   unsigned int GetNumberOfDirections(){ return m_NumberOfDirections(); };
+  std::vector<bool> GetTexturesStatus(){ return m_SelectedTextures; };
 
-
-  inline TOutputValue operator()(const TIter& it)
+  inline OutputType operator()(const TIter& it)
   {
     double length = itk::NumericTraits<double>::NonpositiveMin();
     double width = itk::NumericTraits<double>::max();
-    double sum = 0.;
-    std::vector<double> di(m_NumberOfDirections, 0.);
-    std::vector<double> minSorted(m_RatioMaxConsiderationNumber, width);
-    std::vector<double> maxSorted(m_RatioMaxConsiderationNumber, length);
-    std::vector<double> sti(m_NumberOfDirections, 0.);
+    double SpatialThresholdDouble = static_cast<double>(m_SpatialThreshold);
+    double NumberOfDirectionsDouble = static_cast<double>(m_NumberOfDirections);
+    double dist     = 0.;
+    double angle    = 0.;
+    double sdiVal   = 0.;
+    double sumWMean = 0.;
+    double sum      = 0.;
+    std::vector<double>       di(m_NumberOfDirections, 0.);
+    std::vector<double>       minSorted(m_RatioMaxConsiderationNumber, width);
+    std::vector<double>       maxSorted(m_RatioMaxConsiderationNumber, length);
+    std::vector<double>       sti(m_NumberOfDirections, 0.);
     std::vector<unsigned int> lengthLine(m_NumberOfDirections, 0);
 
     std::vector<double>::iterator itVector;
-    TOutputValue out;
-    out.SetSize(6);
-    out.Fill(0);
-    //std::vector<double> radius(2, static_cast<double>(it.GetRadius()[0]));
-    //radius[1] = static_cast<double>(it.GetRadius()[0]);
- 
+    OutputType out(6, 0);
+  
     OffsetType off;
     off.Fill(0);
-    double SpatialThresholdDouble = static_cast<double>(m_SpatialThreshold);
-    double NumberOfDirectionsDouble = static_cast<double>(m_NumberOfDirections);
-    double dist = 0.;
-    double angle = 0.;
-    double sdiVal = 0.;
-    double sumWMean = 0.;
+  
 
     for( unsigned int d = 0; d<m_NumberOfDirections; d++ )
       {
@@ -103,87 +126,106 @@ public:
 	angle = m_Alpha*static_cast<double>(d);
 
 	// last offset in the diraction respecting spatial threshold
-	off[0] = vcl_floor(SpatialThresholdDouble*vcl_cos( angle ) + 0.5);
-	off[1] = vcl_floor(SpatialThresholdDouble*vcl_sin( angle ) + 0.5);
+	off[0] = static_cast<unsigned int>(vcl_floor(SpatialThresholdDouble*vcl_cos( angle ) + 0.5));
+	off[1] = static_cast<unsigned int>(vcl_floor(SpatialThresholdDouble*vcl_sin( angle ) + 0.5));
 	// last indices in the diration respecting spectral threshold
 	OffsetType offEnd = this->FindLastOffset( it, off );
 	// computes distance = dist between the 2 segment point. One of them is the center pixel -> (0,0)
   	dist = vcl_sqrt( vcl_pow(static_cast<double>(offEnd[0]), 2 ) + vcl_pow(static_cast<double>(offEnd[1]), 2 ) );
 
 	// for length computation
-	if( dist>length)
-	  length = dist;
+	if( m_SelectedTextures[0] == true )
+	  if( dist>length )
+	    length = dist;
 
-	// for width computation 
-	if( dist<width)
-	  width = dist;
+	// for width computation
+ 	if( m_SelectedTextures[1] == true )
+	  if( dist<width)
+	    width = dist;
 
 	// for PSI computation
-	sum += dist;
+	if( m_SelectedTextures[2] == true || m_SelectedTextures[5] == true )
+	  sum += dist;
 
 	// for w-mean computation
-	sdiVal = this->ComputeSDi(it, offEnd);
+	if( m_SelectedTextures[3] == true )
+	  sdiVal = this->ComputeSDi(it, offEnd);
 
 	// for Ratio computation
-	bool doo = false;
-	itVector = maxSorted.begin();
-	while( itVector != maxSorted.end() && doo==false )
+	if( m_SelectedTextures[4] == true )
 	  {
-	    if( dist>(*itVector) )
+	    bool doo = false;
+	    itVector = maxSorted.begin();
+	    while( itVector != maxSorted.end() && doo==false )
 	      {
-		maxSorted.insert(itVector, dist);
-		maxSorted.pop_back();
-		doo=true;
+		if( dist>(*itVector) )
+		  {
+		    maxSorted.insert(itVector, dist);
+		    maxSorted.pop_back();
+		    doo=true;
+		  }
+		++itVector;
 	      }
-	    ++itVector;
-	  }
-	doo = false;
-	itVector = minSorted.begin();
-	while( itVector != minSorted.end() && doo==false )
-	  {
-	    if( dist<(*itVector) )
+	    doo = false;
+	    itVector = minSorted.begin();
+	    while( itVector != minSorted.end() && doo==false )
 	      {
-		minSorted.insert(itVector, dist);
-		minSorted.pop_back();
-		doo=true;
+		if( dist<(*itVector) )
+		  {
+		    minSorted.insert(itVector, dist);
+		    minSorted.pop_back();
+		    doo=true;
+		  }
+		++itVector;
 	      }
-	    ++itVector;
 	  }
+
 	di[d] = dist;
-	lengthLine[d] = dist;//static_cast<unsigned int>( vcl_sqrt(vcl_pow(static_cast<double>(offEnd[0]), 2) + vcl_pow(static_cast<double>(offEnd[1]), 2)) );
-	sti[d] = sdiVal;
-	if(sdiVal!=0.)
-	  sumWMean += (m_Alpha*(dist-1)*dist/*lengthLine[n]*di[n]*/)/sdiVal;
+	if( m_SelectedTextures[3] == true )
+	  {
+	    lengthLine[d] = static_cast<unsigned int>(dist);//static_cast<unsigned int>( vcl_sqrt(vcl_pow(static_cast<double>(offEnd[0]), 2) + vcl_pow(static_cast<double>(offEnd[1]), 2)) );
+	    sti[d] = sdiVal;
+	    if(sdiVal!=0.)
+	      sumWMean += (m_Alpha*(dist-1)*dist/*lengthLine[n]*di[n]*/)/sdiVal;
+	  }
       }
 
     /////// FILL OUTPUT
     // length
-    out[0] = length;
+    if( m_SelectedTextures[0] == true )
+      out[0] = static_cast<OutputValueType>(length);
     // width
-    out[1] = width;
+    if( m_SelectedTextures[1] == true )
+      out[1] = static_cast<OutputValueType>(width);
     // PSI
-    out[2] = sum/NumberOfDirectionsDouble;
+    if( m_SelectedTextures[2] == true )
+      out[2] = static_cast<OutputValueType>(sum/NumberOfDirectionsDouble);
     // w-mean
-    out[3] = sumWMean/NumberOfDirectionsDouble;
+    if( m_SelectedTextures[3] == true )
+      out[3] = static_cast<OutputValueType>(sumWMean/NumberOfDirectionsDouble);
     // ratio
-    double sumMin = 0;
-    double sumMax = 0;
-    for(unsigned int t=0; t<m_RatioMaxConsiderationNumber ; t++)
+    if( m_SelectedTextures[4] == true )
       {
-	sumMin += minSorted[t];
-	sumMax += maxSorted[t];
+	double sumMin = 0;
+	double sumMax = 0;
+	for(unsigned int t=0; t<m_RatioMaxConsiderationNumber ; t++)
+	  {
+	    sumMin += minSorted[t];
+	    sumMax += maxSorted[t];
+	  }
+	if (sumMax != 0.)
+	  out[4] = static_cast<OutputValueType>(vcl_atan(sumMin/sumMax));
+	else if (sumMax == 0. && sumMin == 0.)
+	  out[4] = static_cast<OutputValueType>(1.);
       }
-    if (sumMax != 0.)
-      out[4] = vcl_atan(sumMin/sumMax);
-    else if (sumMax == 0. && sumMin == 0.)
-      out[4] = 1.;
-    
     // SD
-    double sumPSI = 0;
-    for(unsigned int n=0; n<di.size(); n++)
-      sumPSI += vcl_pow(di[n] - out[3] , 2);
-    out[5] = vcl_sqrt(sumPSI)/(NumberOfDirectionsDouble-1.);
-
+    if( m_SelectedTextures[4] == true )
+      {
+	double sumPSI = 0;
+	for(unsigned int n=0; n<di.size(); n++)
+	  sumPSI += vcl_pow(di[n] - sumWMean/NumberOfDirectionsDouble , 2);
+	out[5] = static_cast<OutputValueType>(vcl_sqrt(sumPSI)/(NumberOfDirectionsDouble-1.));
+      }
    
     return out;
 
@@ -197,16 +239,13 @@ public:
   OffsetType FindLastOffset( const TIter & it, const OffsetType & stopOffset )
     {	
       bool res = true;
+      int signX = this->ComputeStep( stopOffset[0] );
+      int signY = this->ComputeStep( stopOffset[1] );
+      
       OffsetType currentOff;
-      currentOff.Fill(0);
-      int signX = 1;
-      if(stopOffset[0]<0)
-	signX = -1;
-      int signY = 1;
-      if(stopOffset[1]<0)
-	signY = -1;
-
+      currentOff.Fill(0);       
       currentOff[0]=signX;
+
       double slop = 0.;
       if(stopOffset[0]!=0)
 	slop = static_cast<double>(stopOffset[1] / static_cast<double>(stopOffset[0]) );
@@ -214,23 +253,16 @@ public:
       bool isInside = true;
       while( isInside == true && res == true)	
 	{
-	  // get the nearest point : Bresenham algo
-	  if(stopOffset[0]!=0)
-	    currentOff[1] = static_cast<int>( vcl_floor( slop*static_cast<double>(currentOff[0]) + 0.5 ));
-	  else
-	    currentOff[1] += signY;
+	  this->ComputePointLine( currentOff, slop, signY, stopOffset[0] );
 
-	  if( vcl_abs(it.GetPixel(currentOff)[0]-it.GetCenterPixel()[0]) > m_SpectralThreshold )
+	  if( vcl_abs(it.GetPixel(currentOff)-it.GetCenterPixel()) > m_SpectralThreshold )
 	    {
 	      res = false;
 	    }
 	  else
 	    currentOff[0]+=signX;
        
-	  if( signX*currentOff[0]>signX*stopOffset[0] && stopOffset[0]!=0)
-	    isInside = false;
-	  else if( signY*currentOff[1]>signY*stopOffset[1] && stopOffset[1] != 0 )
-	    isInside = false;
+	  isInside = this->CheckIsInside(signX, signY, currentOff, stopOffset);
 	}
 
       return currentOff;
@@ -240,46 +272,36 @@ public:
   /** Computes SD in the ith direction */
   double ComputeSDi( const TIter & it, const OffsetType & stopOffset)
     {
-      double SDi = 0.;
-      bool canGo = true;
-      OffsetType currentOff;
-      currentOff.Fill(0);
+      bool         canGo = true;
       unsigned int nbElt = 0;
-      double mean = 0.;
-      double slop = 0.;
-      if(stopOffset[0]!=0)
+      double       SDi   = 0.;  
+      double       mean  = 0.;
+      double       slop  = 0.;
+      if(stopOffset[0] != 0)
 	slop = static_cast<double>(stopOffset[1] / static_cast<double>(stopOffset[0]) );
       
-      int signX = 1;
-      if(stopOffset[0]<0)
-	signX = -1;
-      int signY = 1;
-      if(stopOffset[1]<0)
-	signY = -1;
+      int signX = this->ComputeStep( stopOffset[0] );
+      int signY = this->ComputeStep( stopOffset[1] );
+      
+      OffsetType currentOff;
+      currentOff.Fill(0);   
       currentOff[0]=signX;
 
       bool isInside = true;
       // First compute mean
       while( isInside == true && canGo == true )
 	{
-	  // get the nearest point : Bresenham algo
-	  if(stopOffset[0]!=0)
-	    currentOff[1] = static_cast<int>( vcl_floor( slop*static_cast<double>(currentOff[0]) + 0.5 ));
-	  else
-	    currentOff[1]+=signY;
- 
-	  mean += static_cast<double>(it.GetPixel(currentOff)[0]);
+	  this->ComputePointLine( currentOff, slop, signY, stopOffset[0] );
+
+	  mean += static_cast<double>(it.GetPixel(currentOff));
 	  nbElt++;
 
-	  if( vcl_abs(it.GetPixel(currentOff)[0]-it.GetCenterPixel()[0]) >= m_SpectralThreshold )
+	  if( vcl_abs(it.GetPixel(currentOff)-it.GetCenterPixel()) >= m_SpectralThreshold )
 	      canGo = false;
 	  else
 	    currentOff[0]+=signX;
- 
-	  if( signX*currentOff[0]>signX*stopOffset[0] && stopOffset[0]!=0)
-	    isInside = false;
-	  else if( signY*currentOff[1]>signY*stopOffset[1] && stopOffset[1] != 0 )
-	    isInside = false;
+	  
+	  isInside = this->CheckIsInside(signX, signY, currentOff, stopOffset);
 	}
       
       mean /= static_cast<double>(nbElt);
@@ -289,26 +311,59 @@ public:
 
       while( isInside == true && canGo == true )
 	{
-	  // get the nearest point : Bresenham algo
-	  if(stopOffset[0]!=0)
-	    currentOff[1] = static_cast<int>( vcl_floor( slop*static_cast<double>(currentOff[0]) + 0.5 ));
-	  else
-	    currentOff[1]+=signY;
-	  SDi += vcl_pow((static_cast<double>(it.GetPixel(currentOff)[0]) - mean), 2);
-	  if( vcl_abs(it.GetPixel(currentOff)[0]-it.GetCenterPixel()[0]) >= m_SpectralThreshold )
+	  this->ComputePointLine( currentOff, slop, signY, stopOffset[0] );
+
+	  SDi += vcl_pow((static_cast<double>(it.GetPixel(currentOff)) - mean), 2);
+	  if( vcl_abs(it.GetPixel(currentOff)-it.GetCenterPixel()) >= m_SpectralThreshold )
 	      canGo = false;
 	  else
 	    currentOff[0]+=signX;
 
-	       
-	  if( signX*currentOff[0]>signX*stopOffset[0] && stopOffset[0]!=0)
-	    isInside = false;
-	  else if( signY*currentOff[1]>signY*stopOffset[1] && stopOffset[1] != 0 )
-	    isInside = false;
+	  isInside = this->CheckIsInside(signX, signY, currentOff, stopOffset);
+
 	}
       return vcl_sqrt(SDi);
     }
 
+
+  /** Check if the current offset is inside the stop one. */
+  bool CheckIsInside(const int & signX, const int & signY, const OffsetType & currentOff, const OffsetType & stopOffset)
+    {
+      bool isInside = true;
+      if( signX*currentOff[0]>=signX*stopOffset[0] && stopOffset[0]!=0)
+	isInside = false;
+      else if( signY*currentOff[1]>=signY*stopOffset[1] && stopOffset[1] != 0 )
+	isInside = false;
+
+      return isInside;
+    }
+
+  /** Compute the y coordinate according to a given x coordinate.
+   *  Use the Bresenham algo if the line is not horizontal (stopOffsetX==0).
+   *  Otherwise, it increments of 1 Y.
+   */
+  void ComputePointLine( OffsetType & currentOff, const double & slop, const int & signY, const int & stopOffsetX )
+    {
+      if(stopOffsetX!=0)
+	currentOff[1] = static_cast<int>( vcl_floor( slop*static_cast<double>(currentOff[0]) + 0.5 ));
+      else
+	currentOff[1]+=signY;
+    }
+
+
+  /** Compute the to give at x for line computation.
+   *  according to the signof stopOffset.
+   */
+  int ComputeStep( const int & stopOffset )
+    {
+      int sign = 1;
+      if(stopOffset<0)
+	sign = -1;
+
+      return sign;
+    }
+
+
 protected:
   /** spectral threshold conditon*/
   InternalPixelType m_SpectralThreshold;
@@ -322,6 +377,16 @@ protected:
   unsigned int m_NumberOfDirections;
   /** Angular step between 2 directions (between  and 2*pi). */
   double m_DirectionStep;
+  /** List of wanted textures :
+   *  0: length
+   *  1: width
+   *  2: PSI
+   *  3: w-mean
+   *  4: ratio
+   *  5: SD
+   *  Set to 1 means the texture will be computed.
+   **/
+  std::vector<bool> m_SelectedTextures;
 };
 
 } // end namespace functor
diff --git a/Code/FeatureExtraction/otbLineDirectionImageFilter.h b/Code/FeatureExtraction/otbLineDirectionImageFilter.h
index 5eba6ee8b03aa9f31d4c9da44cc54c9ad8ee27d3..bbdf04230011aaaedef542e26872d17ea3ad541c 100644
--- a/Code/FeatureExtraction/otbLineDirectionImageFilter.h
+++ b/Code/FeatureExtraction/otbLineDirectionImageFilter.h
@@ -15,52 +15,96 @@ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 PURPOSE.  See the above copyright notices for more information.
 
 =========================================================================*/
-
 #ifndef __otbLineDirectionImageFilter_h
 #define __otbLineDirectionImageFilter_h
 
-#include "otbMath.h"
 #include "otbLineDirectionFunctor.h"
-#include "otbUnaryFunctorNeighborhoodImageFilter.h"
+#include "itkImageToImageFilter.h"
+#include "itkImageRegionIteratorWithIndex.h"
 #include "itkConstNeighborhoodIterator.h"
 
 
 namespace otb
 {
-  /** \class LineDirectionImageFilter
-   *  \brief This functor computes water, ndvi and spectral index of an image
-   */
+/** \class LineDirectionImageFilter
+ *  \brief This functor computes the texture describes in the following publication
+ *  It is based on line direction estimation.
+ *
+ * Please refer to Xin Huang, Liangpei Zhang and Pingxiang Li publication,
+ * Classification and Extraction of Spatial Features in Urban Areas
+ * Using High-Resolution Multispectral Imagery.
+ * IEEE Geoscience and Remote Sensing Letters,
+ * vol. 4, n. 2, 2007, pp 260-264
+ *
+ * The texture is computaed for each pixel using its neighborhood.
+ * User can set the spatial threshold taht is the max line length, the spectral threshold
+ * that is the max difference authorized between a pixel of the line and the center pixel
+ * of the current neighborhood. Alpha and RationMaxConsideration are used to compute 
+ * the \omega -mean value. Finally, The number of direction can be precised with 
+ * NumberOfDirections.
+ * You can choose the computed textures using SetTextureStatus method (1:length, 2:width,
+ * 3:PSI, 4:w-mean, 5:ratio, 6:SD).
+ *
+ * \sa LineDirectionFunctor
+*/
 
 template <class TInputImage, class TOutputImage>
 class ITK_EXPORT LineDirectionImageFilter :
-public UnaryFunctorNeighborhoodImageFilter< TInputImage,
-				            TOutputImage, 
-				            Functor::LineDirectionFunctor< ITK_TYPENAME itk::ConstNeighborhoodIterator<TInputImage>, 
-						           		   ITK_TYPENAME TOutputImage::PixelType> >
+public itk::ImageToImageFilter<TInputImage,TOutputImage>
 {
 public:
   /** Standard class typedefs. */
   typedef LineDirectionImageFilter                           Self;
   typedef TInputImage                                        InputImageType;
   typedef TOutputImage                                       OutputImageType;
-  typedef typename InputImageType::PixelType                 InputPixelType;
-  typedef typename InputImageType::InternalPixelType         InputInternalPixelType;
-  typedef typename OutputImageType::PixelType                OutputPixelType;
-  typedef typename otb::UnaryFunctorNeighborhoodImageFilter< TInputImage,
-                                                             TOutputImage, 
-               				                     Functor::LineDirectionFunctor< typename itk::ConstNeighborhoodIterator<TInputImage>, 
-                                                                                            OutputPixelType > > Superclass;
-  typedef itk::SmartPointer<Self>        Pointer;
-  typedef itk::SmartPointer<const Self>  ConstPointer;
+  typedef itk::ImageToImageFilter<TInputImage,TOutputImage >  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(UnaryFunctorNeighborhoodImageFilter,ImageToImageFilter);
+  
+  /** Some convenient typedefs. */
+  typedef typename InputImageType::ConstPointer         InputImagePointerType;
+  typedef typename InputImageType::RegionType           InputImageRegionType;
+  typedef typename InputImageType::PixelType            InputImagePixelType;
+  typedef typename InputImageType::SizeType             InputImageSizeType;
+  typedef typename OutputImageType::Pointer             OutputImagePointerType;
+  typedef typename OutputImageType::RegionType          OutputImageRegionType;
+  typedef typename OutputImageType::PixelType           OutputImagePixelType;
+  typedef itk::ConstNeighborhoodIterator<TInputImage>   NeighborhoodIteratorType;
+  typedef typename NeighborhoodIteratorType::RadiusType RadiusType;
+  typedef Functor::LineDirectionFunctor< NeighborhoodIteratorType,OutputImagePixelType  > FunctorType; 
+  typedef typename FunctorType::OutputType              FunctorOutputType;
+  typedef itk::ProcessObject                            ProcessObjectType;
+
+  /**Set/Get the radius of neighborhood.*/
+  itkGetMacro(Radius,unsigned int);
+
+  /** Functor accessors */
+  FunctorType& GetFunctor()
+  {
+    return m_Functor;
+  }
+
+  const FunctorType& GetFunctor() const
+  {
+    return m_Functor;
+  };
+  void SetFunctor(const FunctorType& functor)
+  {
+    m_Functor = functor;
+    this->Modified();
+  }
+
   /** Spatial Threshold accessor */
   void SetSpatialThreshold( unsigned int thresh )
     { 
       this->GetFunctor().SetSpatialThreshold( thresh ); 
-      this->SetRadius(thresh);
+      m_Radius = thresh;
       this->Modified();
     };
   unsigned int GetSpatialThreshold()
@@ -68,13 +112,13 @@ public:
       return this->GetFunctor().GetSpatialThreshold(); 
     };
   /** Spectral Threshold accessor */
-  void SetSpectralThreshold( InputInternalPixelType thresh )
+  void SetSpectralThreshold( InputImagePixelType thresh )
     { 
       this->GetFunctor().SetSpectralThreshold( thresh ); 
     };
-  InputInternalPixelType GetSpectralThreshold()
+  InputImagePixelType GetSpectralThreshold()
     { 
-      return this->GetFunctor().GSetSpectralThreshold(); 
+      return this->GetFunctor().GetSpectralThreshold(); 
     };
   /** RatioMaxConsiderationNumber accessor */
   void SetRatioMaxConsiderationNumber( unsigned int value )
@@ -106,29 +150,85 @@ public:
   unsigned int GetNumberOfDirections()
     { 
       return this->GetFunctor().GetNumberOfDirections(); 
-    }; 
-
+    };
 
- virtual void GenerateOutputInformation()
+  /** Texture selection accessors 
+   *  1: length
+   *  2: width
+   *  3: PSI
+   *  4: w-mean
+   *  5: ratio
+   *  6: SD
+   *  Set to 1 means the texture will be computed.
+   **/
+  void SetTextureStatus( unsigned int id, bool isSelected )
     {
-      Superclass::GenerateOutputInformation();
-      this->GetOutput()->SetNumberOfComponentsPerPixel(6);
+      if ( id>this->GetTexturesStatus().size() || id == 0 )
+	{
+	  itkExceptionMacro(<<"Invalid texture index "<<id<<", must be in [1;"<<this->GetTexturesStatus().size()<<"]");
+	}
+      else
+	{
+	  this->GetFunctor().SetTextureStatus( id-1, isSelected );
+	}
     }
+  std::vector<bool> GetTexturesStatus()
+    {
+      return this->GetFunctor().GetTexturesStatus();
+    }
+
+  /** Return output length image */
+  const OutputImageType * GetLengthOutput() const;
+  OutputImageType * GetLengthOutput();
+
+   /** Return output width image */
+  const OutputImageType * GetWidthOutput() const;
+  OutputImageType * GetWidthOutput();
+
+  /** Return output PSI image */
+  const OutputImageType * GetPSIOutput() const;
+  OutputImageType * GetPSIOutput();
+  
+  /** Return output WMean image */
+  const OutputImageType * GetWMeanOutput() const;
+  OutputImageType * GetWMeanOutput();
+  
+  /** Return output ratio image */
+  const OutputImageType * GetRatioOutput() const;
+  OutputImageType * GetRatioOutput();
+  
+  /** Return output SD image */
+  const OutputImageType * GetSDOutput() const;
+  OutputImageType * GetSDOutput();
+
+
+ 
+  virtual void GenerateOutputInformation();
+  std::vector<FunctorType> m_FunctorList;
 
 protected:
-  LineDirectionImageFilter()
-    {
-      this->SetRadius(this->GetSpatialThreshold());
-    };
+  LineDirectionImageFilter();
   virtual ~LineDirectionImageFilter(){};
+  void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+  virtual void BeforeThreadedGenerateData();
+  virtual void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId );
+  /** Pad the input requested region by radius */
+  virtual void GenerateInputRequestedRegion(void);
 
 private:
   LineDirectionImageFilter(const Self&); //purposely not implemented
   void operator=(const Self&); //purposely not implemented
-  
+
+  unsigned int m_Radius;
+  FunctorType m_Functor;
+
 };
 
 } // end namespace otb
 
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbLineDirectionImageFilter.txx"
+#endif
 
 #endif
diff --git a/Code/FeatureExtraction/otbLineDirectionImageFilter.txx b/Code/FeatureExtraction/otbLineDirectionImageFilter.txx
new file mode 100644
index 0000000000000000000000000000000000000000..0cca3fa0e41f9add5d764b8a0f954683e87fbd6b
--- /dev/null
+++ b/Code/FeatureExtraction/otbLineDirectionImageFilter.txx
@@ -0,0 +1,451 @@
+/*=========================================================================
+
+  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 __otbLineDirectionImageFilter_txx
+#define __otbLineDirectionImageFilter_txx
+
+#include "otbLineDirectionImageFilter.h"
+
+#include "itkProgressReporter.h"
+#include "itkImageRegionIterator.h"
+#include "itkNeighborhoodAlgorithm.h"
+#include "otbMath.h"
+
+namespace otb
+{
+
+template <class TInputImage, class TOutputImage>
+LineDirectionImageFilter<TInputImage,TOutputImage>
+::LineDirectionImageFilter()
+{
+  this->SetNumberOfInputs(1);
+  this->SetNumberOfRequiredInputs( 1 );
+  this->SetNumberOfOutputs(6);
+  this->SetNumberOfRequiredOutputs(1);
+  
+  this->SetNthOutput(0,OutputImageType::New());
+  this->SetNthOutput(1,OutputImageType::New());
+  this->SetNthOutput(2,OutputImageType::New());
+  this->SetNthOutput(3,OutputImageType::New());
+  this->SetNthOutput(4,OutputImageType::New());
+  this->SetNthOutput(5,OutputImageType::New());
+
+  m_Radius = this->GetSpatialThreshold();
+  m_FunctorList.clear();
+}
+/************************************************************
+ *
+ *              OUTPUTS MANIPULATION
+ *
+ ************************************************************/
+// Return output length image
+template <class TInputImage, class TOutputImage>
+const typename LineDirectionImageFilter<TInputImage, TOutputImage>::OutputImageType *
+LineDirectionImageFilter<TInputImage, TOutputImage>
+::GetLengthOutput() const
+{
+  if (this->GetNumberOfOutputs() < 1)
+  {
+    return 0;
+  }
+  if (this->GetTexturesStatus()[0] == false)
+  {
+    itkExceptionMacro(<<"Impossible to create length image : texture not selected");
+  }
+  return static_cast<const OutputImageType * > (this->itk::ProcessObject::GetOutput(0) );
+}
+template <class TInputImage, class TOutputImage>
+typename LineDirectionImageFilter<TInputImage, TOutputImage>::OutputImageType *
+LineDirectionImageFilter<TInputImage, TOutputImage>
+::GetLengthOutput()
+{
+  if (this->GetNumberOfOutputs() < 1)
+  {
+    return 0;
+  }
+  if (this->GetTexturesStatus()[0] == false)
+  {
+    itkExceptionMacro(<<"Impossible to create length image : texture not selected");
+  }
+  return static_cast<OutputImageType * >(this->itk::ProcessObject::GetOutput(0) );
+}
+
+// Return output width image
+template <class TInputImage, class TOutputImage>
+const typename LineDirectionImageFilter<TInputImage, TOutputImage>::OutputImageType *
+LineDirectionImageFilter<TInputImage, TOutputImage>
+::GetWidthOutput() const
+{
+  if (this->GetNumberOfOutputs() < 2)
+  {
+    return 0;
+  }
+ if (this->GetTexturesStatus()[1] == false)
+  {
+     itkExceptionMacro(<<"Impossible to create width image : texture not selected");
+  }
+  return static_cast<const OutputImageType * > (this->itk::ProcessObject::GetOutput(1) );
+}
+template <class TInputImage, class TOutputImage>
+typename LineDirectionImageFilter<TInputImage, TOutputImage>::OutputImageType *
+LineDirectionImageFilter<TInputImage, TOutputImage>
+::GetWidthOutput()
+{
+  if (this->GetNumberOfOutputs() < 2)
+  {
+    return 0;
+  }
+ if (this->GetTexturesStatus()[1] == false)
+  {
+    itkExceptionMacro(<<"Impossible to create width image : texture not selected");
+  }
+  return static_cast<OutputImageType * >(this->itk::ProcessObject::GetOutput(1) );
+}
+
+//Return output PSI image 
+template <class TInputImage, class TOutputImage>
+const typename LineDirectionImageFilter<TInputImage, TOutputImage>::OutputImageType *
+LineDirectionImageFilter<TInputImage, TOutputImage>
+::GetPSIOutput() const
+{
+  if (this->GetNumberOfOutputs() < 3)
+  {
+    return 0;
+  }
+ if (this->GetTexturesStatus()[2] == false)
+  {
+    itkExceptionMacro(<<"Impossible to create PSI image : texture not selected");
+  }
+  return static_cast<const OutputImageType * > (this->itk::ProcessObject::GetOutput(2) );
+}
+template <class TInputImage, class TOutputImage>
+typename LineDirectionImageFilter<TInputImage, TOutputImage>::OutputImageType *
+LineDirectionImageFilter<TInputImage, TOutputImage>
+::GetPSIOutput()
+{
+  if (this->GetNumberOfOutputs() < 3)
+  {
+    return 0;
+  }
+  if (this->GetTexturesStatus()[2] == false)
+  {
+    itkExceptionMacro(<<"Impossible to create PSI image : texture not selected");
+  }
+
+  return static_cast<OutputImageType * >(this->itk::ProcessObject::GetOutput(2) );
+}
+
+// Return output WMean image
+template <class TInputImage, class TOutputImage>
+const typename LineDirectionImageFilter<TInputImage, TOutputImage>::OutputImageType *
+LineDirectionImageFilter<TInputImage, TOutputImage>
+::GetWMeanOutput() const
+{
+  if (this->GetNumberOfOutputs() < 4)
+  {
+    return 0;
+  }
+  if (this->GetTexturesStatus()[3] == false)
+  {
+    itkExceptionMacro(<<"Impossible to create W-Mean image : texture not selected");
+  }
+  return static_cast<const OutputImageType * > (this->itk::ProcessObject::GetOutput(3) );
+}
+template <class TInputImage, class TOutputImage>
+typename LineDirectionImageFilter<TInputImage, TOutputImage>::OutputImageType *
+LineDirectionImageFilter<TInputImage, TOutputImage>
+::GetWMeanOutput()
+{
+  if (this->GetNumberOfOutputs() < 4)
+  {
+    return 0;
+  }
+ if (this->GetTexturesStatus()[3] == false)
+  {
+    itkExceptionMacro(<<"Impossible to create W-Mean image : texture not selected");
+  }
+  return static_cast<OutputImageType * >(this->itk::ProcessObject::GetOutput(3) );
+}
+
+// Return output ratio image
+template <class TInputImage, class TOutputImage>
+const typename LineDirectionImageFilter<TInputImage, TOutputImage>::OutputImageType *
+LineDirectionImageFilter<TInputImage, TOutputImage>
+::GetRatioOutput() const
+{
+  if (this->GetNumberOfOutputs() < 5)
+  {
+    return 0;
+  }
+ if (this->GetTexturesStatus()[4] == false)
+  {
+    itkExceptionMacro(<<"Impossible to create Ratio image : texture not selected");
+  }
+  return static_cast<const OutputImageType * > (this->itk::ProcessObject::GetOutput(4) );
+}
+template <class TInputImage, class TOutputImage>
+typename LineDirectionImageFilter<TInputImage, TOutputImage>::OutputImageType *
+LineDirectionImageFilter<TInputImage, TOutputImage>
+::GetRatioOutput()
+{
+  if (this->GetNumberOfOutputs() < 5)
+  {
+    return 0;
+  }
+  if (this->GetTexturesStatus()[4] == false)
+  {
+    itkExceptionMacro(<<"Impossible to create Ratio image : texture not selected");
+  }
+  return static_cast<OutputImageType * >(this->itk::ProcessObject::GetOutput(4) );
+}
+
+// Return output SD image
+template <class TInputImage, class TOutputImage>
+const typename LineDirectionImageFilter<TInputImage, TOutputImage>::OutputImageType *
+LineDirectionImageFilter<TInputImage, TOutputImage>
+::GetSDOutput() const
+{
+  if (this->GetNumberOfOutputs() < 6)
+  {
+    return 0;
+  }
+ if (this->GetTexturesStatus()[5] == false)
+  {
+    itkExceptionMacro(<<"Impossible to create SD image : texture not selected");
+  }
+  return static_cast<const OutputImageType * > (this->itk::ProcessObject::GetOutput(5) );
+}
+template <class TInputImage, class TOutputImage>
+typename LineDirectionImageFilter<TInputImage, TOutputImage>::OutputImageType *
+LineDirectionImageFilter<TInputImage, TOutputImage>
+::GetSDOutput()
+{
+  if (this->GetNumberOfOutputs() < 6)
+  {
+    return 0;
+  }
+  if (this->GetTexturesStatus()[5] == false)
+  {
+    itkExceptionMacro(<<"Impossible to create SD image : texture not selected");
+  }
+
+  return static_cast<OutputImageType * >(this->itk::ProcessObject::GetOutput(5) );
+}
+
+
+
+template <class TInputImage, class TOutputImage>
+void
+LineDirectionImageFilter<TInputImage, TOutputImage>
+::BeforeThreadedGenerateData()
+{
+  Superclass::BeforeThreadedGenerateData();
+  if(this->GetSpatialThreshold() < this->GetRatioMaxConsiderationNumber())
+    {
+      itkExceptionMacro(<<"Spatial Threshold ("<<this->GetSpatialThreshold()
+			<<") is lower than Ration Max Consideration Number ("
+			<<this->GetRatioMaxConsiderationNumber()<<") what is not allowed.");
+    }
+  for (int i =0; i<this->GetNumberOfThreads(); i++)
+    {
+      m_FunctorList.push_back(m_Functor);
+    }
+}
+
+template <class TInputImage, class TOutputImage>
+void
+LineDirectionImageFilter<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 outputPtr1 = this->GetOutput(0);//this->GetLengthOutput();
+  typename Superclass::OutputImagePointer outputPtr2 = this->GetOutput(1);//this->GetWidthOutput();
+  typename Superclass::OutputImagePointer outputPtr3 = this->GetOutput(2);//this->GetPSIOutput();
+  typename Superclass::OutputImagePointer outputPtr4 = this->GetOutput(3);//this->GetWMeanOutput();
+  typename Superclass::OutputImagePointer outputPtr5 = this->GetOutput(4);//this->GetRatioOutput();
+  typename Superclass::OutputImagePointer outputPtr6 = this->GetOutput(5);//this->GetSDOutput();
+
+  if ( !inputPtr || !outputPtr1 || !outputPtr2 || !outputPtr3 || !outputPtr4 || !outputPtr5 || !outputPtr6 )
+  {
+    return;
+  }
+  // get a copy of the input requested region (should equal the output
+  // requested region)
+  typename TInputImage::RegionType inputRequestedRegion;
+  inputRequestedRegion = inputPtr->GetRequestedRegion();
+
+  // pad the input requested region by the operator radius
+  InputImageSizeType rad;
+  rad[0] = m_Radius;
+  rad[1] = m_Radius;
+  inputRequestedRegion.PadByRadius( rad );
+
+  // 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__);
+    itk::OStringStream msg;
+    msg << this->GetNameOfClass()
+    << "::GenerateInputRequestedRegion()";
+    e.SetLocation(msg.str().c_str());
+    e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
+    e.SetDataObject(inputPtr);
+    throw e;
+  }
+}
+
+
+
+template <class TInputImage, class TOutputImage>
+void
+LineDirectionImageFilter<TInputImage, TOutputImage>
+::GenerateOutputInformation()
+{
+  Superclass::GenerateOutputInformation();
+  //this->GetOutput()->SetNumberOfComponentsPerPixel(6);
+}
+
+
+template <class TInputImage, class TOutputImage>
+void
+LineDirectionImageFilter<TInputImage, TOutputImage>
+::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
+{
+  itk::ZeroFluxNeumannBoundaryCondition<TInputImage> nbc;
+
+// We use dynamic_cast since inputs are stored as DataObjects.  The
+  // ImageToImageFilter::GetInput(int) always returns a pointer to a
+  // TInputImage so it cannot be used for the second input.
+  InputImagePointerType inputPtr = dynamic_cast<const TInputImage*>(ProcessObjectType::GetInput(0));
+  OutputImagePointerType outputPtr1 = this->GetOutput(0);
+  OutputImagePointerType outputPtr2 = this->GetOutput(1);
+  OutputImagePointerType outputPtr3 = this->GetOutput(2);
+  OutputImagePointerType outputPtr4 = this->GetOutput(3);
+  OutputImagePointerType outputPtr5 = this->GetOutput(4);
+  OutputImagePointerType outputPtr6 = this->GetOutput(5);
+
+  RadiusType r;
+  r.Fill(this->GetRadius());
+  NeighborhoodIteratorType neighInputIt;
+
+  itk::ImageRegionIterator<TOutputImage> outputIt1, outputIt2, outputIt3, outputIt4, outputIt5, outputIt6;
+  FunctorOutputType outputFunctor;
+
+  // Find the data-set boundary "faces"
+  typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<TInputImage>::FaceListType faceList;
+  typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<TInputImage> bC;
+  faceList = bC(inputPtr, outputRegionForThread, r);
+
+
+  typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<TInputImage>::FaceListType::iterator fit;
+
+  // support progress methods/callbacks
+  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
+
+  // Process each of the boundary faces.  These are N-d regions which border
+  // the edge of the buffer.
+
+  std::vector<bool> textStatus = this->GetTexturesStatus();
+  for (fit=faceList.begin(); fit != faceList.end(); ++fit)
+  {
+    neighInputIt = itk::ConstNeighborhoodIterator<TInputImage>(r, inputPtr, *fit);
+    
+    outputIt1 = itk::ImageRegionIterator<TOutputImage>(outputPtr1, *fit);
+    outputIt2 = itk::ImageRegionIterator<TOutputImage>(outputPtr2, *fit);
+    outputIt3 = itk::ImageRegionIterator<TOutputImage>(outputPtr3, *fit);
+    outputIt4 = itk::ImageRegionIterator<TOutputImage>(outputPtr4, *fit);
+    outputIt5 = itk::ImageRegionIterator<TOutputImage>(outputPtr5, *fit);
+    outputIt6 = itk::ImageRegionIterator<TOutputImage>(outputPtr6, *fit);
+    
+    std::vector< itk::ImageRegionIterator<TOutputImage> *> outItList;
+    outItList.push_back(&outputIt1);
+    outItList.push_back(&outputIt2);
+    outItList.push_back(&outputIt3);
+    outItList.push_back(&outputIt4);
+    outItList.push_back(&outputIt5);
+    outItList.push_back(&outputIt6);
+    
+    neighInputIt.OverrideBoundaryCondition(&nbc);
+    neighInputIt.GoToBegin();
+    
+    for(unsigned int i = 0; i<outItList.size(); i++)
+      {
+	(*outItList[i]).GoToBegin();
+      }
+    
+    while ( !outputIt1.IsAtEnd() )
+    {
+
+      outputFunctor = m_FunctorList[threadId]( neighInputIt);    
+      for(unsigned int i = 0; i<outItList.size(); i++)
+	{
+	  if( textStatus[i]==true )
+	    (*outItList[i]).Set( outputFunctor[i] );
+	}	
+      
+      ++neighInputIt;
+      for(unsigned int i = 0; i<outItList.size(); i++)
+	{
+	  ++(*outItList[i]);
+	}
+      
+      progress.CompletedPixel();
+    }
+  }
+}
+
+
+/**
+ * Standard "PrintSelf" method
+ */
+template <class TInputImage, class TOutputImage>
+void
+LineDirectionImageFilter<TInputImage, TOutputImage>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+  Superclass::PrintSelf( os, indent );
+
+  
+  //os << indent << "Spatial Threshold             : "  << this->GetSpatialThreshold() << std::endl;
+  //os << indent << "Spectral Threshold            : "  << this->GetSpectralThreshold() << std::endl;
+  //os << indent << "Ratio Max Consideration Number: "  << this->GetRatioMaxConsiderationNumber() << std::endl;
+  //os << indent << "Alpha                         : "  << this->GetAlpha() << std::endl;
+  //os << indent << "Number Of Directions          : "  << this->GetNumberOfDirections() << std::endl;
+  
+}
+
+
+} // end namespace otb
+
+
+#endif
diff --git a/Code/FeatureExtraction/otbLineSegmentDetector.h b/Code/FeatureExtraction/otbLineSegmentDetector.h
index b36f9f0cd4544f6d05c67c47e85fe83a46b34987..e17c88a72b9d6aed55c6f25e2ba1acc557cd5031 100644
--- a/Code/FeatureExtraction/otbLineSegmentDetector.h
+++ b/Code/FeatureExtraction/otbLineSegmentDetector.h
@@ -60,8 +60,9 @@ namespace Functor
 	
 	inline TOutputPixel operator()(const TInputPixel& input)
 	  {
-	    TOutputPixel resp = vcl_atan2(input[0],-input[1]);
- 	    if (resp<0)
+	    TOutputPixel resp = static_cast<TOutputPixel>(vcl_atan2(input[0],-input[1]));
+	    
+ 	    if (resp< itk::NumericTraits<TOutputPixel>::Zero)
  	      {
  		resp = -resp;
  	      }
@@ -79,7 +80,7 @@ namespace Functor
  * 
  */
 
-template <class TInputImage,class TPrecision>
+template <class TInputImage,class TPrecision = double>
 class ITK_EXPORT LineSegmentDetector :
       public otb::ImageToLineSpatialObjectListFilter< TInputImage >
 {
@@ -130,20 +131,20 @@ public:
   typedef typename DirectionVectorType::iterator                        DirectionVectorIteratorType; 
 
   /** */ 
-  typedef itk::GradientRecursiveGaussianImageFilter<InputImageType > GradientFilterType;
+  typedef itk::GradientRecursiveGaussianImageFilter<OutputImageType > GradientFilterType;
   //typedef itk::GradientImageFilter<InputImageType > GradientFilterType;
   typedef typename GradientFilterType::Pointer GradientFilterPointerType;
   typedef typename GradientFilterType::OutputImageType GradientOutputImageType;
 
 
-  typedef itk::UnaryFunctorImageFilter<GradientOutputImageType,InputImageType,
-  Functor::MagnitudeFunctor<typename GradientOutputImageType::PixelType,typename  InputImageType::PixelType> > MagnitudeFilterType;
+  typedef itk::UnaryFunctorImageFilter<GradientOutputImageType,OutputImageType,
+  Functor::MagnitudeFunctor<typename GradientOutputImageType::PixelType,TPrecision> > MagnitudeFilterType;
   typedef typename MagnitudeFilterType::Pointer                                   MagnitudeFilterPointerType;
   typedef typename MagnitudeFilterType::OutputImageType::PixelType                MagnitudePixelType;
   typedef typename MagnitudeFilterType::OutputImageType                           MagnitudeImageType;
             
-  typedef itk::UnaryFunctorImageFilter<GradientOutputImageType,InputImageType,
-  Functor::OrientationFunctor<typename GradientOutputImageType::PixelType,typename InputImageType::PixelType> > OrientationFilterType;
+  typedef itk::UnaryFunctorImageFilter<GradientOutputImageType,OutputImageType,
+  Functor::OrientationFunctor<typename GradientOutputImageType::PixelType,TPrecision> > OrientationFilterType;
   typedef typename OrientationFilterType::Pointer OrientationFilterPointerType;
   typedef typename OrientationFilterType::OutputImageType                         OutputImageDirType;
   typedef typename OutputImageDirType::RegionType                                 OutputImageDirRegionType ;
@@ -153,13 +154,11 @@ public:
   typedef typename LabelImageType::Pointer                          LabelImagePointerType;
   
   /** Vector to store the rectangle characteization  center, width, orientation ,( begin ,end ) of the central line*/
-  typedef std::vector<float>                                        RectangleType;
+  typedef std::vector<double>                                       RectangleType;
   typedef typename RectangleType::iterator                          RectangleIteratorType;
   typedef std::vector< RectangleType>                               RectangleListType;
   typedef typename RectangleListType::iterator                      RectangleListTypeIterator; 
 
-  //typedef otb::ImageFileWriter<MagnitudeImageType  >                writerType;
-  //typedef typename writerType::Pointer                              writerPointerType;
 
 protected:
   LineSegmentDetector();
@@ -188,38 +187,38 @@ protected:
   virtual void GrowRegion(InputIndexType  index);
   
   /** Define if two are aligned */
-  virtual bool IsAligned(float Angle, float regionAngle, float prec);
+  virtual bool IsAligned(double Angle, double regionAngle, double prec);
   
   /** For each region of the region List it builds a rectangle */
   virtual int ComputeRectangles();
   
   /** */
-  virtual void Region2Rect(IndexVectorType  region , float angleRegion);
+  virtual void Region2Rect(IndexVectorType  region , double angleRegion);
   
   /** */
-  virtual float ComputeRegionOrientation(IndexVectorType  region , float x, float y , float angleRegion);
+  virtual double ComputeRegionOrientation(IndexVectorType  region , double x, double y , double angleRegion);
 
   /** */
-  virtual float angle_diff(float a, float b);
+  virtual double angle_diff(double a, double b);
 
   /**  Compute the Number Of False Alarm for a rectangle*/
-  virtual float ComputeRectNFA(RectangleType  rec);
+  virtual double ComputeRectNFA(RectangleType  rec);
 
   /** */
-  virtual float ImproveRectangle(RectangleType  * rec);
+  virtual double ImproveRectangle(RectangleType  * rec);
 
   /** NFA For a rectangle*/
-  virtual float NFA(int n, int k, double p, double logNT);
+  virtual double NFA(int n, int k, double p, double logNT);
   
   /** Create a copy of a rectangle*/
   virtual void CopyRectangle(RectangleType * rDst , RectangleType  *rSrc );
   
 
   /** Rutines from numerical recipes*/
-  virtual float betacf(float a, float b, float x);
-  virtual float gammln(float xx);
-  virtual float betai(float a, float b, float x);
-  virtual float factln(int n);
+  virtual double betacf(double a, double b, double x);
+  virtual double gammln(double xx);
+  virtual double betai(double a, double b, double x);
+  virtual double factln(int n);
 
   /** Printself method*/
   void PrintSelf(std::ostream& os, itk::Indent indent) const;
@@ -233,14 +232,14 @@ private:
   LabelImagePointerType             m_UsedPointImage;
   RectangleListType                 m_RectangleList;
   
-  float                             m_Threshold;
-  float                             m_Prec;
-  float                             m_DirectionsAllowed;
-  unsigned int                      m_MinimumRegionSize;
-  unsigned int                      m_NumberOfImagePixels;
-
-  int                      m_Length;
-  int                      m_Width;
+  double                             m_Threshold;
+  double                             m_Prec;
+  double                             m_DirectionsAllowed;
+  unsigned int                       m_MinimumRegionSize;
+  unsigned int                       m_NumberOfImagePixels;
+
+  int                                m_Length;
+  int                                m_Width;
   
   /** Gradient filter */
   GradientFilterPointerType m_GradientFilter;
diff --git a/Code/FeatureExtraction/otbLineSegmentDetector.txx b/Code/FeatureExtraction/otbLineSegmentDetector.txx
index 30c52aa4e971d9ec353d7ea34d822a3c31c1c63e..960e024397786c7f1865250b05f154b92ae06e36 100644
--- a/Code/FeatureExtraction/otbLineSegmentDetector.txx
+++ b/Code/FeatureExtraction/otbLineSegmentDetector.txx
@@ -25,6 +25,7 @@
 #include "itkConstNeighborhoodIterator.h"
 #include "itkNeighborhoodIterator.h"
 #include "otbPolygon.h"
+#include "itkCastImageFilter.h"
 
 #include "otbMath.h"
 
@@ -81,15 +82,20 @@ LineSegmentDetector<TInputImage,TPrecision >
   m_UsedPointImage->Allocate();
   m_UsedPointImage->FillBuffer(0);
 
+  /** Cast the MagnitudeOutput Image in */
+  typedef itk::CastImageFilter<InputImageType, OutputImageType>      castFilerType;
+  typename castFilerType::Pointer  castFilter =  castFilerType::New();
+  castFilter->SetInput(this->GetInput());
+  
   /** Compute the modulus and the orientation gradient image*/
-  m_GradientFilter->SetInput(this->GetInput());
+  m_GradientFilter->SetInput(castFilter->GetOutput());
   m_GradientFilter->SetSigma(1.1);
   m_MagnitudeFilter->SetInput(m_GradientFilter->GetOutput());
-    
   m_OrientationFilter->SetInput(m_GradientFilter->GetOutput());
+  
   m_MagnitudeFilter->Update();
   m_OrientationFilter->Update();
-  
+
   /** Comupute the seed histogram to begin the search*/
   CoordinateHistogramType   CoordinateHistogram;
   CoordinateHistogram = this->SortImageByModulusValue(m_MagnitudeFilter->GetOutput());
@@ -126,15 +132,15 @@ LineSegmentDetector<TInputImage,TPrecision >
   /** 
    *  Compute the minimum region size 
    */
-  float logNT = 5.*(vcl_log10(static_cast<float>(m_Width)) + vcl_log10(static_cast<float>(m_Length)))/2.;
-  float log1_p = vcl_log10(m_DirectionsAllowed);
-  float rapport = logNT/log1_p;
+  double logNT = 5.*(vcl_log10(static_cast<double>(m_Width)) + vcl_log10(static_cast<double>(m_Length)))/2.;
+  double log1_p = vcl_log10(m_DirectionsAllowed);
+  double rapport = logNT/log1_p;
   m_MinimumRegionSize = -1*static_cast<unsigned  int>(rapport);
 
   
   /** Definition of the min & the max of an image*/
-  OutputPixelType   min = itk::NumericTraits<MagnitudePixelType >::Zero;
-  OutputPixelType   max = itk::NumericTraits<MagnitudePixelType >::Zero;
+  OutputPixelType   min = itk::NumericTraits</*MagnitudePixelType*/TPrecision>::Zero;
+  OutputPixelType   max = itk::NumericTraits</*MagnitudePixelType*/TPrecision>::Zero;
   
   /** Computing the min & max of the image*/
   typedef  itk::MinimumMaximumImageCalculator<OutputImageType>  MinMaxCalculatorFilter;
@@ -151,14 +157,13 @@ LineSegmentDetector<TInputImage,TPrecision >
 
   /** Computing the length of the bins*/  
   unsigned int NbBin = 10;  
-  MagnitudePixelType lengthBin = (max - min)/static_cast<MagnitudePixelType>(NbBin-1) ;
+  double lengthBin = static_cast<double>((max - min))/static_cast<double>(NbBin-1) ;
   CoordinateHistogramType  tempHisto(NbBin);  /** Initializing the histogram */
 
   
   itk::ImageRegionIterator<OutputImageType> it(modulusImage, 
 					       modulusImage->GetRequestedRegion());
-    
-  //std::cout << "min "<< min<<  " max " << max  << " lengthBin " << lengthBin <<  " et thresh" << m_Threshold << std::endl;
+
   it.GoToBegin();
   while(!it.IsAtEnd())
     {
@@ -241,7 +246,7 @@ LineSegmentDetector<TInputImage, TPrecision>
   RectangleListTypeIterator            itRec = m_RectangleList.begin();
   while(itRec != m_RectangleList.end())
     {
-      float NFA = this->ImproveRectangle(&(*itRec) );
+      double NFA = this->ImproveRectangle(&(*itRec) );
       /**
        * Here we start building the OUTPUT :a LineSpatialObjectList. 
        */
@@ -251,9 +256,9 @@ LineSegmentDetector<TInputImage, TPrecision>
 	  PointListType pointList;
 	  PointType     point;
 	  
-	  point.SetPosition((*itRec)[0],(*itRec)[1]);
+	  point.SetPosition(static_cast<TPrecision>((*itRec)[0]),static_cast<TPrecision>((*itRec)[1]));
 	  pointList.push_back(point);
-	  point.SetPosition((*itRec)[2],(*itRec)[3]);
+	  point.SetPosition(static_cast<TPrecision>((*itRec)[2]),static_cast<TPrecision>((*itRec)[3]));
 	  pointList.push_back(point);
 	  
 	  typename LineSpatialObjectType::Pointer line = LineSpatialObjectType::New();
@@ -295,17 +300,17 @@ LineSegmentDetector<TInputImage, TPrecision>
  * the components of the rectangle
  */
 template <class TInputImage, class TPrecision  >
-float
+double
 LineSegmentDetector<TInputImage, TPrecision>
 :: ImproveRectangle(RectangleType  * rec)
 {
   int n = 0;
-  float nfa_new;
-  float delta = 0.5;
-  float delta_2 = delta / 2.0;
+  double nfa_new;
+  double delta = 0.5;
+  double delta_2 = delta / 2.0;
   RectangleType r;
 
-  float NFA = this->ComputeRectNFA(*rec);
+  double NFA = this->ComputeRectNFA(*rec);
   
   if( NFA > 0. ) return NFA;
  
@@ -457,7 +462,7 @@ LineSegmentDetector<TInputImage, TPrecision>
   this->SetPixelToUsed(index);
   
   /** Neighborhooding */
-  typedef itk::ConstNeighborhoodIterator<InputImageType>  NeighborhoodIteratorType;
+  typedef itk::ConstNeighborhoodIterator<OutputImageType>  NeighborhoodIteratorType;
   typename NeighborhoodIteratorType::SizeType             radius;
   radius.Fill(1);
   NeighborhoodIteratorType                                itNeigh(radius,m_MagnitudeFilter->GetOutput(),
@@ -470,12 +475,12 @@ LineSegmentDetector<TInputImage, TPrecision>
   IndexVectorType                                 reg;
 
   /** Angle of the region*/
-  float regionAngle = 0;
+  double regionAngle = 0;
   
   /** Add the first point to the region */
   reg.push_back(index);
-  float sumX = 0.;
-  float sumY = 0.;
+  double sumX = 0.;
+  double sumY = 0.;
   
   /** 
    * Loop for searching regions 
@@ -492,7 +497,7 @@ LineSegmentDetector<TInputImage, TPrecision>
       while(s < neighSize )
 	{
 	  InputIndexType NeighIndex = itNeigh.GetIndex(s);
-	  float angleComp =   itNeighDir.GetPixel(s);
+	  double angleComp =   itNeighDir.GetPixel(s);
 
 	  if( !this->IsUsed(NeighIndex) && this->IsAligned(angleComp, regionAngle, m_Prec) )
 	    {
@@ -521,9 +526,9 @@ LineSegmentDetector<TInputImage, TPrecision>
 template <class TInputImage, class TPrecision  >
 bool
 LineSegmentDetector<TInputImage, TPrecision>
-::IsAligned(float Angle, float regionAngle, float prec)
+::IsAligned(double Angle, double regionAngle, double prec)
 {
-  float diff = Angle - regionAngle;
+  double diff = Angle - regionAngle;
   
   if( diff < 0.0 ) diff = -diff;
   if( diff > 1.5*M_PI )
@@ -542,15 +547,15 @@ LineSegmentDetector<TInputImage, TPrecision>
 template <class TInputImage, class TPrecision  >
 void
 LineSegmentDetector<TInputImage, TPrecision>
-::Region2Rect(IndexVectorType  region , float angleRegion)
+::Region2Rect(IndexVectorType  region , double angleRegion)
 {
   /** Local Variables*/
-  float weight = 0.,sumWeight = 0.; 
-  float  x = 0., y = 0.;
-  float l_min = 0., l_max = 0.,l =0., w=0. , w_min = 0. , w_max =0.;
+  double weight = 0.,sumWeight = 0.; 
+  double  x = 0., y = 0.;
+  double l_min = 0., l_max = 0.,l =0., w=0. , w_min = 0. , w_max =0.;
   
   /** Neighborhooding again*/
-  typedef itk::ConstNeighborhoodIterator<InputImageType>  NeighborhoodIteratorType;
+  typedef itk::ConstNeighborhoodIterator<OutputImageType>  NeighborhoodIteratorType;
   typename NeighborhoodIteratorType::SizeType             radius;
   radius.Fill(0);
   NeighborhoodIteratorType                                itNeigh(radius,m_MagnitudeFilter->GetOutput(),
@@ -562,8 +567,8 @@ LineSegmentDetector<TInputImage, TPrecision>
     {
       itNeigh.SetLocation(*it);
       weight = *itNeigh.GetCenterValue(); 
-      x += static_cast<float>((*it)[0])* weight;
-      y += static_cast<float>((*it)[1])* weight;
+      x += static_cast<double>((*it)[0])* weight;
+      y += static_cast<double>((*it)[1])* weight;
       sumWeight +=  weight;
       ++it;
     }
@@ -578,7 +583,7 @@ LineSegmentDetector<TInputImage, TPrecision>
     }
   
   /** Compute the orientation of the region*/
-  float theta = this->ComputeRegionOrientation(region , x , y , angleRegion); 
+  double theta = this->ComputeRegionOrientation(region , x , y , angleRegion); 
     
   /* Lenght & Width of the rectangle **/
   
@@ -588,50 +593,50 @@ LineSegmentDetector<TInputImage, TPrecision>
   MagnitudeVector sum_l( 2*Diagonal, itk::NumericTraits<MagnitudePixelType>::Zero);
   MagnitudeVector sum_w( 2*Diagonal, itk::NumericTraits<MagnitudePixelType>::Zero);
   
-  float dx = vcl_cos(theta);
-  float dy = vcl_sin(theta);
+  double dx = vcl_cos(theta);
+  double dy = vcl_sin(theta);
   
   it = region.begin();
   while(it != region.end())
     {
       itNeigh.SetLocation(*it);
       weight = *itNeigh.GetCenterValue(); 
-      l =  ( static_cast<float>((*it)[0]) - x )*dx + ( static_cast<float>((*it)[1]) - y )*dy;
-      w = -( static_cast<float>((*it)[0]) - x )*dy + ( static_cast<float>((*it)[1]) - y )*dx;
+      l =  ( static_cast<double>((*it)[0]) - x )*dx + ( static_cast<double>((*it)[1]) - y )*dy;
+      w = -( static_cast<double>((*it)[0]) - x )*dy + ( static_cast<double>((*it)[1]) - y )*dx;
       
       if(l<l_min) l_min = l;   if(l>l_max)  l_max = l;
       if(w<w_min) w_min = w;   if(w>w_max)  w_max = w;
 
-      sum_l[static_cast< int>(vcl_floor(l)+0.5)+ Diagonal ] +=  weight;
-      sum_w[static_cast< int>(vcl_floor(w)+0.5)+ Diagonal ] +=  weight;
+      sum_l[static_cast< int>(vcl_floor(l)+0.5)+ Diagonal ] +=  static_cast<MagnitudePixelType>(weight);
+      sum_w[static_cast< int>(vcl_floor(w)+0.5)+ Diagonal ] +=  static_cast<MagnitudePixelType>(weight);
       
       ++it;
     }
   
   /** Thresholdinq the width and the length*/
   
-  float sum_th = 0.01 * sumWeight; /* weight threshold for selecting region */
-  float s = 0.;
+  double sum_th = 0.01 * sumWeight; /* weight threshold for selecting region */
+  double s = 0.;
   int i = 0;
   
   for( s=0.0,i = static_cast<int>(l_min); s<sum_th && i<=static_cast<int>(l_max) ; i++) 
     s += sum_l[ Diagonal + i];
   
-  float lb = (static_cast<float>(i-1) - 0.5 );
+  double lb = (static_cast<double>(i-1) - 0.5 );
   
   for(s=0.0,i=static_cast<int>(l_max); s<sum_th && i>=static_cast<int>(l_min); i--) 
     s += sum_l[ Diagonal  + i];
-  float lf =  (static_cast<float>(i+1) + 0.5 ); 
+  double lf =  (static_cast<double>(i+1) + 0.5 ); 
 
   for(s=0.0,i=static_cast<int>(w_min); s<sum_th && i<=static_cast<int>(w_max); i++) 
     s += sum_w[ Diagonal + i];
   
-  float wr = (static_cast<float>(i-1) - 0.5);  
+  double wr = (static_cast<double>(i-1) - 0.5);  
 
   for(s=0.0,i=static_cast<int>(w_max); s<sum_th && i>=static_cast<int>(w_min); i--) 
     s += sum_w[Diagonal  + i];
   
-  float wl = (static_cast<float>(i+1) + 0.5 );  
+  double wl = (static_cast<double>(i+1) + 0.5 );  
 
 
   /** Finally store the rectangle in vector this way : 
@@ -665,19 +670,19 @@ LineSegmentDetector<TInputImage, TPrecision>
  */
 
 template <class TInputImage, class TPrecision  >
-float
+double
 LineSegmentDetector<TInputImage, TPrecision>
-::ComputeRegionOrientation(IndexVectorType  region , float x,float y  , float angleRegion) 
+::ComputeRegionOrientation(IndexVectorType  region , double x,double y  , double angleRegion) 
 {
 
-  float Ixx = 0.0;
-  float Iyy = 0.0;
-  float Ixy = 0.0;
-  float theta = 0.;
-  float weight = 0.,sum = 0.;
+  double Ixx = 0.0;
+  double Iyy = 0.0;
+  double Ixy = 0.0;
+  double theta = 0.;
+  double weight = 0.,sum = 0.;
   
   /** Neighborhooding again*/
-  typedef itk::ConstNeighborhoodIterator<InputImageType>  NeighborhoodIteratorType;
+  typedef itk::ConstNeighborhoodIterator<OutputImageType>  NeighborhoodIteratorType;
   typename NeighborhoodIteratorType::SizeType             radius;
   radius.Fill(0);
   NeighborhoodIteratorType                                itNeigh(radius,m_MagnitudeFilter->GetOutput(),
@@ -689,8 +694,8 @@ LineSegmentDetector<TInputImage, TPrecision>
     {
       itNeigh.SetLocation(*it);
       weight = *itNeigh.GetCenterValue(); 
-      float Iyy2 = static_cast<float>((*it)[0]) - x;
-      float Ixx2 = static_cast<float>((*it)[1]) - y;
+      double Iyy2 = static_cast<double>((*it)[0]) - x;
+      double Ixx2 = static_cast<double>((*it)[1]) - y;
       
       Ixx += Ixx2*Ixx2*weight;
       Iyy += Iyy2*Iyy2*weight;
@@ -700,8 +705,8 @@ LineSegmentDetector<TInputImage, TPrecision>
     }
 
   /** using te itk Eigen analysis*/
-  typedef itk::Matrix<float , 2, 2>    MatrixType;
-  typedef std::vector<float >    MatrixEigenType;
+  typedef itk::Matrix<double , 2, 2>    MatrixType;
+  typedef std::vector<double >    MatrixEigenType;
   MatrixType Inertie, eigenVector;
   MatrixEigenType  eigenMatrix(2,0.);
   Inertie[1][1] =Ixx;
@@ -728,9 +733,9 @@ LineSegmentDetector<TInputImage, TPrecision>
  * Compute the difference betwenn 2 angles modulo 2_PI
  */
 template <class TInputImage, class TPrecision  >
-float 
+double 
 LineSegmentDetector<TInputImage, TPrecision>
-::angle_diff(float a, float b)
+::angle_diff(double a, double b)
 {
   a -= b;
   while( a <= -M_PI ) a += 2*M_PI;
@@ -744,16 +749,16 @@ LineSegmentDetector<TInputImage, TPrecision>
  * compute the number of false alarm of the rectangle 
  */
 template <class TInputImage, class TPrecision  >
-float 
+double 
 LineSegmentDetector<TInputImage, TPrecision>
 ::ComputeRectNFA(RectangleType  rec)
 {
   int NbAligned = 0;
-  float nfa_val = 0.;
+  double nfa_val = 0.;
   
-  float dx = vcl_cos(rec[5]);
-  float dy = vcl_sin(rec[5]);
-  float halfWidth = rec[4]/2;
+  double dx = vcl_cos(rec[5]);
+  double dy = vcl_sin(rec[5]);
+  double halfWidth = rec[4]/2;
   
   /** Determine the four corners of the rectangle*/
   /** Remember : 
@@ -785,7 +790,7 @@ LineSegmentDetector<TInputImage, TPrecision>
    */
   
   /**  Compute the number of points aligned */
-  typedef otb::Polygon<float>       PolygonType;
+  typedef otb::Polygon<double>       PolygonType;
   PolygonType::Pointer  rectangle = PolygonType::New();
 
   /** Fill the rectangle with the points*/
@@ -818,7 +823,7 @@ LineSegmentDetector<TInputImage, TPrecision>
     }
   
   /** Compute the NFA from the rectangle computed below*/
-  float logNT = 5.*(vcl_log10(static_cast<float>(m_Length)) + vcl_log10(static_cast<float>(m_Width)))/2.;
+  double logNT = 5.*(vcl_log10(static_cast<double>(m_Length)) + vcl_log10(static_cast<double>(m_Width)))/2.;
   
   nfa_val = NFA(pts,NbAligned ,m_DirectionsAllowed,logNT);
   
@@ -835,11 +840,11 @@ LineSegmentDetector<TInputImage, TPrecision>
    logNT - logarithm of Number of Tests
  */
 template <class TInputImage, class TPrecision  >
-float 
+double 
 LineSegmentDetector<TInputImage, TPrecision>
 ::NFA(int n, int k, double p, double logNT)
 {
-  float val;
+  double val;
 
   if(k==0) return -logNT;
 
@@ -876,12 +881,12 @@ LineSegmentDetector<TInputImage, TPrecision>
 /**************************************************************************************************************/
 /**************************************************************************************************************/
 template <class TInputImage, class TPrecision  >
-float 
+double 
 LineSegmentDetector<TInputImage, TPrecision>
-::betacf(float a, float b, float x)
+::betacf(double a, double b, double x)
 {
   int m,m2;
-  float aa,c,d,del,h,qab,qam,qap;
+  double aa,c,d,del,h,qab,qam,qap;
 
   qab=a+b;
   qap=a+1.0;
@@ -915,9 +920,9 @@ LineSegmentDetector<TInputImage, TPrecision>
 }
 
 template <class TInputImage, class TPrecision  >
-float 
+double 
 LineSegmentDetector<TInputImage, TPrecision>
-::gammln(float xx)
+::gammln(double xx)
 {
   double x,y,tmp,ser;
   static double cof[6]={76.18009172947146,-86.50532032941677,
@@ -934,13 +939,13 @@ LineSegmentDetector<TInputImage, TPrecision>
 }
 
 template <class TInputImage, class TPrecision  >
-float 
+double 
 LineSegmentDetector<TInputImage, TPrecision>
-::betai(float a, float b, float x)
+::betai(double a, double b, double x)
 {
-  float betacf(float a, float b, float x);
-  float gammln(float xx);
-  float bt;
+  double betacf(double a, double b, double x);
+  double gammln(double xx);
+  double bt;
 
   if (x == 0.0 || x == 1.0) bt=0.0;
   else
@@ -953,12 +958,12 @@ LineSegmentDetector<TInputImage, TPrecision>
 
 
 template <class TInputImage, class TPrecision  >
-float 
+double 
 LineSegmentDetector<TInputImage, TPrecision>
 ::factln(int n)
 {
-  float gammln(float xx);
-  static float a[101];
+  double gammln(double xx);
+  static double a[101];
 
   if (n <= 1) return 0.0;
   if (n <= 100) return a[n] ? a[n] : (a[n]=this->gammln(n+1.0));
diff --git a/Code/FeatureExtraction/otbLineSpatialObjectListToRightAnglePointSetFilter.txx b/Code/FeatureExtraction/otbLineSpatialObjectListToRightAnglePointSetFilter.txx
index f09f975807aa12abdaa17c2bf9d4b06227bce7a0..6e78364cc11ce7929c4067de93a86059c4371ff6 100644
--- a/Code/FeatureExtraction/otbLineSpatialObjectListToRightAnglePointSetFilter.txx
+++ b/Code/FeatureExtraction/otbLineSpatialObjectListToRightAnglePointSetFilter.txx
@@ -159,20 +159,20 @@ LineSpatialObjectListToRightAnglePointSetFilter<TImage,TLinesList ,TPointSet>
   
   /** Verify if the indexes of the line are inside the region*/
   typename InputImageType::SizeType size = this->GetInputImage()->GetRequestedRegion().GetSize();
-  unsigned int width = size[0];
-  unsigned int length = size[1];
+  //unsigned int width = size[0];
+  //unsigned int length = size[1];
   
-  if(IndexBeginSrc[0] > width)  IndexBeginSrc[0] = width-1;
-  if(IndexBeginSrc[1] > length) IndexBeginSrc[1] = length-1;
+//   if(IndexBeginSrc[0] > width)  IndexBeginSrc[0] = width-1;
+//   if(IndexBeginSrc[1] > length) IndexBeginSrc[1] = length-1;
 
-  if(IndexEndSrc[0] > width)   IndexEndSrc[0] = width-1;
-  if(IndexEndSrc[1] > length)  IndexEndSrc[1] = length-1;
+//   if(IndexEndSrc[0] > width)   IndexEndSrc[0] = width-1;
+//   if(IndexEndSrc[1] > length)  IndexEndSrc[1] = length-1;
 
-  if(IndexBeginSrc[0] <0) IndexBeginSrc[0] = 0;
-  if(IndexBeginSrc[1] <0) IndexBeginSrc[1] = 0;
+//   if(IndexBeginSrc[0] <0) IndexBeginSrc[0] = 0;
+//   if(IndexBeginSrc[1] <0) IndexBeginSrc[1] = 0;
 
-  if(IndexEndSrc[0]<0)   IndexEndSrc[0] = 0;
-  if(IndexEndSrc[1]<0)   IndexEndSrc[1] = 0;
+//   if(IndexEndSrc[0]<0)   IndexEndSrc[0] = 0;
+//   if(IndexEndSrc[1]<0)   IndexEndSrc[1] = 0;
 
  
   /** Extract Indexes from the Dst line to instantiate the line iterator*/
@@ -288,7 +288,7 @@ LineSpatialObjectListToRightAnglePointSetFilter<TImage,TLinesList ,TPointSet>
   /** Extract Indexes from the Dst line to instantiate the line iterator*/
   typename LineType::PointListType &pointsListDst = lineDst->GetPoints();
   typename LineType::PointListType::const_iterator itPointsDst = pointsListDst.begin();
-
+  
   float Xq1 = (*itPointsDst).GetPosition()[0];  //xq1
   float Yq1 = (*itPointsDst).GetPosition()[1];  //yq1
   ++itPointsDst;
diff --git a/Code/FeatureExtraction/otbPanTexTextureImageFilter.txx b/Code/FeatureExtraction/otbPanTexTextureImageFilter.txx
index 5a8f6d4efac91c8b8f4f75e0cd6ddc1161a0f053..774d7e7ef4268416009008ff840510adcebb0bc6 100644
--- a/Code/FeatureExtraction/otbPanTexTextureImageFilter.txx
+++ b/Code/FeatureExtraction/otbPanTexTextureImageFilter.txx
@@ -36,75 +36,6 @@ PanTexTextureImageFilter<TInputImage,TOutputImage>
 }
 
 
-// template <class TInputImage, class TOutputImage, class TFunction  >
-// void
-// UnaryFunctorNeighborhoodWithOffsetImageFilter<TInputImage,TOutputImage,TFunction>
-// ::BeforeThreadedGenerateData()
-// {
-//   Superclass::BeforeThreadedGenerateData();
-
-//   for (int i =0; i<this->GetNumberOfThreads(); i++)
-//   {
-//     m_FunctorList.push_back(m_Functor);
-//   }
-// }
-
-
-// template <class TInputImage, class TOutputImage, class TFunction  >
-// void
-// UnaryFunctorNeighborhoodWithOffsetImageFilter<TInputImage,TOutputImage,TFunction>
-// ::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();
-
-//   // pad the input requested region by the operator radius
-//   InputImageSizeType maxRad;
-//   maxRad[0] = m_Radius + vcl_abs(m_Offset[0]);
-//   maxRad[1] = m_Radius + vcl_abs(m_Offset[1]);;
-//   inputRequestedRegion.PadByRadius( maxRad );
-
-//   // 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__);
-//     itk::OStringStream msg;
-//     msg << this->GetNameOfClass()
-//     << "::GenerateInputRequestedRegion()";
-//     e.SetLocation(msg.str().c_str());
-//     e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
-//     e.SetDataObject(inputPtr);
-//     throw e;
-//   }
-// }
-
-
 
 /**
  * ThreadedGenerateData Performs the neighborhood-wise operation
diff --git a/Code/FeatureExtraction/otbTextureFunctorBase.h b/Code/FeatureExtraction/otbTextureFunctorBase.h
index 35b037e9dc800a6727d32e60bb794caf901d1edc..50152c7ecba4367b6c80481bf417c291724e15d9 100755
--- a/Code/FeatureExtraction/otbTextureFunctorBase.h
+++ b/Code/FeatureExtraction/otbTextureFunctorBase.h
@@ -22,6 +22,9 @@
 #include "itkMacro.h"
 #include "itkNumericTraits.h"
 #include "itkNeighborhood.h"
+#include "itkOffset.h"
+#include "itkSize.h"
+#include "itkVariableLengthVector.h"
 
 
 namespace otb
@@ -41,7 +44,7 @@ namespace Functor
  *  \ingroup Statistics
  */
 
-template <class TIterInput, class TOutput>
+template <class TScalarInputPixelType, class TScalarOutputPixelType>
 class TextureFunctorBase
 {
 public:
@@ -62,17 +65,17 @@ public:
   };
   virtual ~TextureFunctorBase() {};
 
-  typedef TIterInput                           IterType;
-  typedef TOutput                              OutputType;
-  typedef typename IterType::OffsetType        OffsetType;
-  typedef typename IterType::RadiusType        RadiusType;
-  typedef typename OutputType::ValueType       OutputPixelType;
-  typedef typename IterType::InternalPixelType InternalPixelType;
-  typedef typename IterType::ImageType         ImageType;
-  typedef itk::Neighborhood<InternalPixelType, ::itk::GetImageDimension<ImageType>::ImageDimension>    NeighborhoodType;
-  typedef std::vector<double>                   DoubleVectorType;
-  typedef std::vector<int>                      IntVectorType;
-  typedef std::vector<IntVectorType>            IntVectorVectorType;
+  typedef TScalarInputPixelType                       InputScalarType;
+  typedef TScalarOutputPixelType                      OutputScalarType;
+  typedef itk::VariableLengthVector<InputScalarType>  InputVectorType;
+  typedef itk::VariableLengthVector<OutputScalarType> OutputVectorType;
+  typedef itk::Offset<>                               OffsetType;
+  typedef itk::Size<>                                 RadiusType;
+  typedef itk::Neighborhood<InputScalarType, 2>       NeighborhoodType;
+  typedef itk::Neighborhood<InputVectorType, 2>       NeighborhoodVectorType;
+  typedef std::vector<double>                         DoubleVectorType;
+  typedef std::vector<int>                            IntVectorType;
+  typedef std::vector<IntVectorType>                  IntVectorVectorType;
 
 
   void SetOffset(OffsetType off){ m_Offset=off; };
@@ -164,21 +167,64 @@ public:
       m_OffsetBinLength = scottCoef*binLengthOff;     
     }
 
-  inline TOutput operator()(const IterType &itOff)
+inline OutputScalarType operator()(const NeighborhoodType &neigh)
+  {
+    RadiusType radiusOff = neigh.GetRadius();
+    //OutputScalarType outPix;
+    //outPix.SetSize( neigh.GetCenterPixel().GetSize() );
+    //outPix.Fill(0);
+    OffsetType offset;
+    offset.Fill(0);
+    // Compute the neighborhood radius from the neigh+offset iterator to extract the neighborhood area from the iterator
+    RadiusType radius;
+    radius[0] = static_cast<unsigned int>( 0.5*( static_cast<double>(neigh.GetSize()[0] - 1) ) - static_cast<double>( vcl_abs(m_Offset[0])) );
+    radius[1] = static_cast<unsigned int>( 0.5*( static_cast<double>(neigh.GetSize()[1] - 1) ) - static_cast<double>( vcl_abs(m_Offset[1])) );
+    
+    NeighborhoodType inNeigh;
+    inNeigh.SetRadius(radius);
+    NeighborhoodType offNeigh;
+    offNeigh.SetRadius(radiusOff);
+    // Extract the neighborhood area 
+    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] =  neigh[offset];//neigh.GetPixel(offset);
+	  }
+      }
+    // Extract the offset area
+    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] =  neigh[offset];//neigh.GetPixel(offset);
+	  }
+      }
+    OutputScalarType outPix = static_cast<OutputScalarType>( this->ComputeOverSingleChannel(inNeigh, offNeigh) );
+    
+    return outPix;
+  }
+
+  inline OutputVectorType operator()(const NeighborhoodVectorType &neigh)
     { 
-      RadiusType radiusOff = itOff.GetRadius();
-      OutputType outPix;
-      outPix.SetSize( itOff.GetCenterPixel().GetSize() );
+      RadiusType radiusOff = neigh.GetRadius();
+      OutputVectorType outPix;
+      outPix.SetSize( neigh.GetCenterValue/*Pixel*/().GetSize() );
       outPix.Fill(0);
       OffsetType offset;
       offset.Fill(0);
       // Compute the neighborhood radius from the neigh+offset iterator to extract the neighborhood area from the iterator
       RadiusType radius;
-      radius[0] = static_cast<unsigned int>( 0.5*( static_cast<double>(itOff.GetSize()[0] - 1) ) - static_cast<double>( vcl_abs(m_Offset[0])) );
-      radius[1] = static_cast<unsigned int>( 0.5*( static_cast<double>(itOff.GetSize()[1] - 1) ) - static_cast<double>( vcl_abs(m_Offset[1])) );
+      radius[0] = static_cast<unsigned int>( 0.5*( static_cast<double>(neigh.GetSize()[0] - 1) ) - static_cast<double>( vcl_abs(m_Offset[0])) );
+      radius[1] = static_cast<unsigned int>( 0.5*( static_cast<double>(neigh.GetSize()[1] - 1) ) - static_cast<double>( vcl_abs(m_Offset[1])) );
   
       // For each channel
-      for ( unsigned int i=0; i<itOff.GetCenterPixel().GetSize(); i++ )
+      for ( unsigned int i=0; i<neigh.GetCenterValue/*Pixel*/().GetSize(); i++ )
 	{
 	  NeighborhoodType inNeigh;
 	  inNeigh.SetRadius(radius);
@@ -191,7 +237,7 @@ public:
 	      for ( int k = -static_cast<int>(radius[1]); k <= static_cast<int>(radius[1]); k++)
 		{
 		  offset[1] = k;
-		  inNeigh[offset] = itOff.GetPixel(offset)[i];
+		  inNeigh[offset] = neigh[offset][i];//neigh.GetPixel(offset)[i];
 		}
 	    }
 	  // Extract the offset area
@@ -202,10 +248,10 @@ public:
 	      for ( int k = -static_cast<int>(radiusOff[1]); k <= static_cast<int>(radiusOff[1]); k++)
 		{
 		  offset[1] = k;
-		  offNeigh[offset] = itOff.GetPixel(offset)[i];
+		  offNeigh[offset] = neigh[offset][i];
 		}
 	    }
-	  outPix[i] = static_cast<OutputPixelType>( this->ComputeOverSingleChannel(inNeigh, offNeigh) );
+	  outPix[i] = static_cast<OutputScalarType>( this->ComputeOverSingleChannel(inNeigh, offNeigh) );
 	}
       return outPix;
     }
diff --git a/Code/FeatureExtraction/otbTextureImageFunction.txx b/Code/FeatureExtraction/otbTextureImageFunction.txx
index b8c6a6364890ec7cf94cdd48f19ea42a896f0c75..e80d82e3c1735e7e7fef022dab4deac124ccfdeb 100644
--- a/Code/FeatureExtraction/otbTextureImageFunction.txx
+++ b/Code/FeatureExtraction/otbTextureImageFunction.txx
@@ -69,8 +69,6 @@ TextureImageFunction<TInputImage, TFunctor, TCoordRep>
     return ( itk::NumericTraits<RealType>::max() );
   }
 
-  IteratorType 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]);
@@ -80,7 +78,7 @@ TextureImageFunction<TInputImage, TFunctor, TCoordRep>
   FunctorType funct;
   funct.SetOffset(m_Offset);
   
-  return static_cast<RealType>(funct.ComputeOverSingleChannel( it.GetNeighborhood(), itOff.GetNeighborhood()) );
+  return static_cast<RealType>(funct( itOff.GetNeighborhood() ));
 }
 
 
diff --git a/Code/IO/otbGDALImageIO.cxx b/Code/IO/otbGDALImageIO.cxx
index 1bf5a872b655609ed06e7629c1abfd911ca9db9e..254d92b56135e92d6d6865ae5e29df3adcbb27ea 100644
--- a/Code/IO/otbGDALImageIO.cxx
+++ b/Code/IO/otbGDALImageIO.cxx
@@ -34,7 +34,7 @@
 #include "otbGDALImageIO.h"
 #include "otbMacro.h"
 #include "otbSystem.h"
-#include "otbImageBase.h"
+// #include "otbImageMetadata.h"
 #include "otbImage.h"
 
 #include "itkMetaDataObject.h"
diff --git a/Code/IO/otbImage.h b/Code/IO/otbImage.h
index 82f015a2a065acadbd09589962c3510f23aade4b..16de1b2f32a9db62755d8bf8a76bc73aa21146ec 100644
--- a/Code/IO/otbImage.h
+++ b/Code/IO/otbImage.h
@@ -23,7 +23,7 @@
 #endif
 
 #include "itkImage.h"
-#include "otbImageBase.h"
+#include "otbImageMetadataInterface.h"
 
 #include <string.h>
 
@@ -34,11 +34,8 @@ namespace otb
  *
  */
 
-// Le 3ieme parametre template est bidon MAIS necessaire pour compiler avec Microsoft Visual C++ 6.0
 template <class TPixel,unsigned int VImageDimension>
-class ITK_EXPORT Image : public itk::Image<TPixel,VImageDimension>,
-      public ImageBase
-
+class ITK_EXPORT Image : public itk::Image<TPixel,VImageDimension>
 {
 public:
   /** Standard class typedefs. */
@@ -48,8 +45,8 @@ public:
   typedef itk::SmartPointer<const Self>  ConstPointer;
   typedef itk::WeakPointer<const Self>  ConstWeakPointer;
 
-  typedef ImageBase::VectorType  VectorType;
-  typedef ImageBase::ImageKeywordlistType  ImageKeywordlistType;
+  typedef ImageMetadataInterface::VectorType  VectorType;
+  typedef ImageMetadataInterface::ImageKeywordlistType  ImageKeywordlistType;
 
   /** Method for creation through the object factory. */
   itkNewMacro(Self);
@@ -192,7 +189,7 @@ protected:
 private:
   Image(const Self&); //purposely not implemented
   void operator=(const Self&); //purposely not implemented
-
+  typename ImageMetadataInterface::Pointer m_ImageMetadataInterface;
 };
 
 
diff --git a/Code/IO/otbImage.txx b/Code/IO/otbImage.txx
index 3bcd3d8bbb5c4bdea9000d0f606303fbfc177cdb..f18ebeccaadff4b1fbc833838c3646f6b7657b99 100644
--- a/Code/IO/otbImage.txx
+++ b/Code/IO/otbImage.txx
@@ -30,114 +30,115 @@ namespace otb
 template <class TPixel, unsigned int VImageDimension>
 Image<TPixel,VImageDimension>::Image()
 {
+  m_ImageMetadataInterface = ImageMetadataInterface::New();
 }
 
 template <class TPixel, unsigned int VImageDimension>
 std::string Image<TPixel, VImageDimension>::GetProjectionRef( void ) const
 {
-  return ( this->ImageBase::GetProjectionRef( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetProjectionRef( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 std::string Image<TPixel, VImageDimension>::GetGCPProjection( void ) const
 {
-  return ( this->ImageBase::GetGCPProjection( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetGCPProjection( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 unsigned int Image<TPixel, VImageDimension>::GetGCPCount( void  ) const
 {
-  return ( this->ImageBase::GetGCPCount( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetGCPCount( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 OTB_GCP & Image<TPixel, VImageDimension>::GetGCPs ( unsigned int GCPnum )
 {
-  return (this->ImageBase::GetGCPs( this->GetMetaDataDictionary(), GCPnum ) );
+  return (m_ImageMetadataInterface->GetGCPs( this->GetMetaDataDictionary(), GCPnum ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 std::string Image<TPixel, VImageDimension>::GetGCPId( unsigned int GCPnum  ) const
 {
-  return ( this->ImageBase::GetGCPId( this->GetMetaDataDictionary(), GCPnum ) );
+  return ( m_ImageMetadataInterface->GetGCPId( this->GetMetaDataDictionary(), GCPnum ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 std::string Image<TPixel, VImageDimension>::GetGCPInfo( unsigned int GCPnum ) const
 {
-  return ( this->ImageBase::GetGCPInfo( this->GetMetaDataDictionary(), GCPnum ) );
+  return ( m_ImageMetadataInterface->GetGCPInfo( this->GetMetaDataDictionary(), GCPnum ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 double Image<TPixel, VImageDimension>::GetGCPRow( unsigned int GCPnum ) const
 {
-  return ( this->ImageBase::GetGCPRow( this->GetMetaDataDictionary(), GCPnum ) );
+  return ( m_ImageMetadataInterface->GetGCPRow( this->GetMetaDataDictionary(), GCPnum ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 double Image<TPixel, VImageDimension>::GetGCPCol( unsigned int GCPnum ) const
 {
-  return ( this->ImageBase::GetGCPCol( this->GetMetaDataDictionary(), GCPnum ) );
+  return ( m_ImageMetadataInterface->GetGCPCol( this->GetMetaDataDictionary(), GCPnum ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 double Image<TPixel, VImageDimension>::GetGCPX( unsigned int GCPnum ) const
 {
-  return ( this->ImageBase::GetGCPX( this->GetMetaDataDictionary(), GCPnum ) );
+  return ( m_ImageMetadataInterface->GetGCPX( this->GetMetaDataDictionary(), GCPnum ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 double Image<TPixel, VImageDimension>::GetGCPY( unsigned int GCPnum ) const
 {
-  return ( this->ImageBase::GetGCPY( this->GetMetaDataDictionary(), GCPnum ) );
+  return ( m_ImageMetadataInterface->GetGCPY( this->GetMetaDataDictionary(), GCPnum ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 double Image<TPixel, VImageDimension>::GetGCPZ( unsigned int GCPnum ) const
 {
-  return ( this->ImageBase::GetGCPZ( this->GetMetaDataDictionary(), GCPnum ) );
+  return ( m_ImageMetadataInterface->GetGCPZ( this->GetMetaDataDictionary(), GCPnum ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
-ImageBase::VectorType Image<TPixel, VImageDimension>::GetGeoTransform( void ) const
+ImageMetadataInterface::VectorType Image<TPixel, VImageDimension>::GetGeoTransform( void ) const
 {
-  return ( this->ImageBase::GetGeoTransform( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetGeoTransform( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
-ImageBase::VectorType Image<TPixel, VImageDimension>::GetUpperLeftCorner( void ) const
+ImageMetadataInterface::VectorType Image<TPixel, VImageDimension>::GetUpperLeftCorner( void ) const
 {
-  return ( this->ImageBase::GetUpperLeftCorner( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetUpperLeftCorner( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
-ImageBase::VectorType Image<TPixel, VImageDimension>::GetUpperRightCorner( void ) const
+ImageMetadataInterface::VectorType Image<TPixel, VImageDimension>::GetUpperRightCorner( void ) const
 {
-  return ( this->ImageBase::GetUpperRightCorner( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetUpperRightCorner( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
-ImageBase::VectorType Image<TPixel, VImageDimension>::GetLowerLeftCorner( void ) const
+ImageMetadataInterface::VectorType Image<TPixel, VImageDimension>::GetLowerLeftCorner( void ) const
 {
-  return ( this->ImageBase::GetLowerLeftCorner( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetLowerLeftCorner( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
-ImageBase::VectorType Image<TPixel, VImageDimension>::GetLowerRightCorner( void ) const
+ImageMetadataInterface::VectorType Image<TPixel, VImageDimension>::GetLowerRightCorner( void ) const
 {
-  return ( this->ImageBase::GetLowerRightCorner( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetLowerRightCorner( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
-ImageBase::ImageKeywordlistType Image<TPixel, VImageDimension>::GetImageKeywordlist( void )
+ImageMetadataInterface::ImageKeywordlistType Image<TPixel, VImageDimension>::GetImageKeywordlist( void )
 {
-  return ( this->ImageBase::GetImageKeywordlist( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetImageKeywordlist( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
-const ImageBase::ImageKeywordlistType Image<TPixel, VImageDimension>::GetImageKeywordlist( void ) const
+const ImageMetadataInterface::ImageKeywordlistType Image<TPixel, VImageDimension>::GetImageKeywordlist( void ) const
 {
-  return ( this->ImageBase::GetImageKeywordlist( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetImageKeywordlist( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
@@ -155,7 +156,7 @@ Image<TPixel, VImageDimension>::PrintSelf(std::ostream& os, itk::Indent indent)
 {
   Superclass::PrintSelf(os,indent);
 
-  this->ImageBase::PrintSelf( os, indent, this->GetMetaDataDictionary() );
+  m_ImageMetadataInterface->PrintSelf( os, indent, this->GetMetaDataDictionary() );
 
 }
 
diff --git a/Code/IO/otbImageFileWriter.txx b/Code/IO/otbImageFileWriter.txx
old mode 100755
new mode 100644
diff --git a/Code/IO/otbImageBase.cxx b/Code/IO/otbImageMetadataInterface.cxx
similarity index 76%
rename from Code/IO/otbImageBase.cxx
rename to Code/IO/otbImageMetadataInterface.cxx
index 46c0f52c0c17c538302937805e57b679bab9d21f..f003f0193ddb4c29c0cc1cd27b6af429e9e7a4d4 100644
--- a/Code/IO/otbImageBase.cxx
+++ b/Code/IO/otbImageMetadataInterface.cxx
@@ -25,19 +25,19 @@
 
 #include "itkMetaDataObject.h"
 
-#include "otbImageBase.h"
+#include "otbImageMetadataInterface.h"
 #include "otbImageKeywordlist.h"
 
 namespace otb
 {
 
 
-ImageBase::ImageBase()
+ImageMetadataInterface::ImageMetadataInterface()
 {
 
 }
 
-std::string ImageBase::GetProjectionRef( const MetaDataDictionaryType & dict ) const
+std::string ImageMetadataInterface::GetProjectionRef( const MetaDataDictionaryType & dict ) const
 {
   std::string metadata;
 
@@ -50,7 +50,7 @@ std::string ImageBase::GetProjectionRef( const MetaDataDictionaryType & dict ) c
     return ("");
 }
 
-std::string ImageBase::GetGCPProjection( const MetaDataDictionaryType & dict ) const
+std::string ImageMetadataInterface::GetGCPProjection( const MetaDataDictionaryType & dict ) const
 {
   std::string metadata;
 
@@ -63,7 +63,7 @@ std::string ImageBase::GetGCPProjection( const MetaDataDictionaryType & dict ) c
     return ("");
 }
 
-unsigned int ImageBase::GetGCPCount( const MetaDataDictionaryType & dict) const
+unsigned int ImageMetadataInterface::GetGCPCount( const MetaDataDictionaryType & dict) const
 {
   unsigned int GCPCount = 0;
 
@@ -76,7 +76,7 @@ unsigned int ImageBase::GetGCPCount( const MetaDataDictionaryType & dict) const
 }
 
 
-OTB_GCP & ImageBase::GetGCPs( MetaDataDictionaryType & dict, unsigned int GCPnum )
+OTB_GCP & ImageMetadataInterface::GetGCPs( MetaDataDictionaryType & dict, unsigned int GCPnum )
 {
   std::string key;
 
@@ -93,7 +93,7 @@ OTB_GCP & ImageBase::GetGCPs( MetaDataDictionaryType & dict, unsigned int GCPnum
 
 }
 
-std::string ImageBase::GetGCPId( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const
+std::string ImageMetadataInterface::GetGCPId( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const
 {
   std::string key;
 
@@ -111,7 +111,7 @@ std::string ImageBase::GetGCPId( const MetaDataDictionaryType & dict, unsigned i
     return ("");
 }
 
-std::string ImageBase::GetGCPInfo( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const
+std::string ImageMetadataInterface::GetGCPInfo( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const
 {
   std::string key;
 
@@ -129,7 +129,7 @@ std::string ImageBase::GetGCPInfo( const MetaDataDictionaryType & dict, unsigned
     return ("");
 }
 
-double ImageBase::GetGCPRow( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const
+double ImageMetadataInterface::GetGCPRow( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const
 {
   std::string key;
 
@@ -148,7 +148,7 @@ double ImageBase::GetGCPRow( const MetaDataDictionaryType & dict, unsigned int G
     return (0);
 }
 
-double ImageBase::GetGCPCol( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const
+double ImageMetadataInterface::GetGCPCol( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const
 {
   std::string key;
 
@@ -167,7 +167,7 @@ double ImageBase::GetGCPCol( const MetaDataDictionaryType & dict, unsigned int G
     return (0);
 }
 
-double ImageBase::GetGCPX( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const
+double ImageMetadataInterface::GetGCPX( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const
 {
   std::string key;
 
@@ -185,7 +185,7 @@ double ImageBase::GetGCPX( const MetaDataDictionaryType & dict, unsigned int GCP
     return (0);
 }
 
-double ImageBase::GetGCPY( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const
+double ImageMetadataInterface::GetGCPY( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const
 {
   std::string key;
 
@@ -203,7 +203,7 @@ double ImageBase::GetGCPY( const MetaDataDictionaryType & dict, unsigned int GCP
     return (0);
 }
 
-double ImageBase::GetGCPZ( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const
+double ImageMetadataInterface::GetGCPZ( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const
 {
   std::string key;
 
@@ -221,7 +221,7 @@ double ImageBase::GetGCPZ( const MetaDataDictionaryType & dict, unsigned int GCP
     return (0);
 }
 
-ImageBase::VectorType ImageBase::GetGeoTransform( const MetaDataDictionaryType & dict ) const
+ImageMetadataInterface::VectorType ImageMetadataInterface::GetGeoTransform( const MetaDataDictionaryType & dict ) const
 {
   VectorType adfGeoTransform;
 
@@ -232,7 +232,7 @@ ImageBase::VectorType ImageBase::GetGeoTransform( const MetaDataDictionaryType &
   return ( adfGeoTransform );
 }
 
-ImageBase::VectorType ImageBase::GetUpperLeftCorner( const MetaDataDictionaryType & dict ) const
+ImageMetadataInterface::VectorType ImageMetadataInterface::GetUpperLeftCorner( const MetaDataDictionaryType & dict ) const
 {
   VectorType UpperLeftCorner;
 
@@ -243,7 +243,7 @@ ImageBase::VectorType ImageBase::GetUpperLeftCorner( const MetaDataDictionaryTyp
   return ( UpperLeftCorner );
 }
 
-ImageBase::VectorType ImageBase::GetUpperRightCorner( const MetaDataDictionaryType & dict ) const
+ImageMetadataInterface::VectorType ImageMetadataInterface::GetUpperRightCorner( const MetaDataDictionaryType & dict ) const
 {
   VectorType UpperRightCorner;
 
@@ -254,7 +254,7 @@ ImageBase::VectorType ImageBase::GetUpperRightCorner( const MetaDataDictionaryTy
   return ( UpperRightCorner );
 }
 
-ImageBase::VectorType ImageBase::GetLowerLeftCorner( const MetaDataDictionaryType & dict ) const
+ImageMetadataInterface::VectorType ImageMetadataInterface::GetLowerLeftCorner( const MetaDataDictionaryType & dict ) const
 {
   VectorType LowerLeftCorner;
 
@@ -265,7 +265,7 @@ ImageBase::VectorType ImageBase::GetLowerLeftCorner( const MetaDataDictionaryTyp
   return ( LowerLeftCorner );
 }
 
-ImageBase::VectorType ImageBase::GetLowerRightCorner( const MetaDataDictionaryType & dict ) const
+ImageMetadataInterface::VectorType ImageMetadataInterface::GetLowerRightCorner( const MetaDataDictionaryType & dict ) const
 {
   VectorType LowerRightCorner;
 
@@ -276,7 +276,7 @@ ImageBase::VectorType ImageBase::GetLowerRightCorner( const MetaDataDictionaryTy
   return ( LowerRightCorner );
 }
 
-ImageBase::ImageKeywordlistType ImageBase::GetImageKeywordlist( MetaDataDictionaryType & dict )
+ImageMetadataInterface::ImageKeywordlistType ImageMetadataInterface::GetImageKeywordlist( MetaDataDictionaryType & dict )
 {
   ImageKeywordlistType ImageKeywordlist;
 
@@ -287,7 +287,7 @@ ImageBase::ImageKeywordlistType ImageBase::GetImageKeywordlist( MetaDataDictiona
   return ( ImageKeywordlist );
 }
 
-const ImageBase::ImageKeywordlistType ImageBase::GetImageKeywordlist(const MetaDataDictionaryType & dict ) const
+const ImageMetadataInterface::ImageKeywordlistType ImageMetadataInterface::GetImageKeywordlist(const MetaDataDictionaryType & dict ) const
 {
   ImageKeywordlistType ImageKeywordlist;
 
@@ -299,7 +299,7 @@ const ImageBase::ImageKeywordlistType ImageBase::GetImageKeywordlist(const MetaD
 }
 
 void
-ImageBase::PrintSelf(std::ostream& os, itk::Indent indent, const MetaDataDictionaryType & dict) const
+ImageMetadataInterface::PrintSelf(std::ostream& os, itk::Indent indent, const MetaDataDictionaryType & dict) const
 {
 
   std::vector<std::string> keys = dict.GetKeys();
diff --git a/Code/IO/otbImageBase.h b/Code/IO/otbImageMetadataInterface.h
similarity index 68%
rename from Code/IO/otbImageBase.h
rename to Code/IO/otbImageMetadataInterface.h
index 8c22b530798a60134a93e48d6b099a0099fcdb00..6aef34ee94f1afb35ece7d6b8cea623f2d2dd328 100644
--- a/Code/IO/otbImageBase.h
+++ b/Code/IO/otbImageMetadataInterface.h
@@ -15,8 +15,8 @@
      PURPOSE.  See the above copyright notices for more information.
 
 =========================================================================*/
-#ifndef __otbImageBase_h
-#define __otbImageBase_h
+#ifndef __otbImageMetadataInterface_h
+#define __otbImageMetadataInterface_h
 
 #if defined(_MSC_VER)
 #pragma warning ( disable : 4786 )
@@ -31,16 +31,26 @@
 
 namespace otb
 {
-/** \class ImageBase
+/** \class ImageMetadataInterface
  *
- * \brief Creation of an "otb" ImageBase that gets metadata.
+ * \brief Creation of an "otb" ImageMetadataInterface that gets metadata.
  *
  */
-class ITK_EXPORT ImageBase
+class ITK_EXPORT ImageMetadataInterface : public itk::Object
 {
 public:
 
-  typedef ImageBase Self;
+  typedef ImageMetadataInterface Self;
+  typedef itk::Object 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(ImageMetadataInterface, itk::Object);
+
 
   typedef itk::MetaDataDictionary   MetaDataDictionaryType;
 
@@ -50,77 +60,56 @@ public:
 
 
   /** Get the projection coordinate system of the image. */
-
   std::string GetProjectionRef( const MetaDataDictionaryType & dict ) const;
-  virtual std::string GetProjectionRef( void ) const = 0;
 
   /** Get the GCP projection coordinates of the image. */
-
   std::string GetGCPProjection( const MetaDataDictionaryType & dict ) const;
-  virtual std::string GetGCPProjection( void ) const = 0;
 
   unsigned int GetGCPCount( const MetaDataDictionaryType & dict ) const;
-  virtual unsigned int GetGCPCount( void ) const = 0;
 
   OTB_GCP & GetGCPs( MetaDataDictionaryType & dict, unsigned int GCPnum );
-  virtual OTB_GCP & GetGCPs( unsigned int GCPnum ) = 0;
 
   std::string GetGCPId( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const;
-  virtual std::string GetGCPId( unsigned int GCPnum ) const = 0;
 
   std::string GetGCPInfo( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const;
-  virtual std::string GetGCPInfo( unsigned int GCPnum ) const = 0;
 
   double GetGCPRow( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const;
-  virtual double GetGCPRow( unsigned int GCPnum ) const = 0;
 
   double GetGCPCol( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const;
-  virtual double GetGCPCol( unsigned int GCPnum ) const = 0;
 
   double GetGCPX( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const;
-  virtual double GetGCPX( unsigned int GCPnum ) const = 0;
 
   double GetGCPY( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const;
-  virtual double GetGCPY( unsigned int GCPnum ) const = 0;
 
   double GetGCPZ( const MetaDataDictionaryType & dict, unsigned int GCPnum ) const;
-  virtual double GetGCPZ( unsigned int GCPnum ) const = 0;
 
   /** Get the six coefficients of affine geoTtransform. */
 
   VectorType GetGeoTransform( const MetaDataDictionaryType & dict ) const;
-  virtual VectorType GetGeoTransform( void ) const = 0;
 
   /** Get image corners. */
 
   VectorType GetUpperLeftCorner( const MetaDataDictionaryType & dict ) const;
-  virtual VectorType GetUpperLeftCorner() const = 0;
 
   VectorType GetUpperRightCorner( const MetaDataDictionaryType & dict ) const;
-  virtual VectorType GetUpperRightCorner() const = 0;
 
   VectorType GetLowerLeftCorner( const MetaDataDictionaryType & dict ) const;
-  virtual VectorType GetLowerLeftCorner() const = 0;
 
   VectorType GetLowerRightCorner( const MetaDataDictionaryType & dict ) const;
-  virtual VectorType GetLowerRightCorner() const = 0;
 
   /** Get the ImageKeywordlist */
   ImageKeywordlistType GetImageKeywordlist( MetaDataDictionaryType & dict );
   const ImageKeywordlistType GetImageKeywordlist(const MetaDataDictionaryType & dict ) const;
 
-  virtual ImageKeywordlistType GetImageKeywordlist() = 0;
-  virtual const ImageKeywordlistType GetImageKeywordlist(void) const = 0;
-
   void PrintSelf(std::ostream& os, itk::Indent indent, const MetaDataDictionaryType & dict) const;
 
 
 protected:
-  ImageBase();
-  virtual ~ImageBase() {};
+  ImageMetadataInterface();
+  virtual ~ImageMetadataInterface() {};
 
 private:
-  ImageBase(const Self&); //purposely not implemented
+  ImageMetadataInterface(const Self&); //purposely not implemented
   void operator=(const Self&); //purposely not implemented
 
   OTB_GCP m_GCP;
diff --git a/Code/IO/otbTileMapImageIO.cxx b/Code/IO/otbTileMapImageIO.cxx
index 37bf0a69c7477d7e30444fab587a3c1d7afcc105..cef7187aa0c4e6a517ed6399a607afc9c2d995eb 100644
--- a/Code/IO/otbTileMapImageIO.cxx
+++ b/Code/IO/otbTileMapImageIO.cxx
@@ -34,7 +34,6 @@
 #include "otbTileMapImageIO.h"
 #include "otbMacro.h"
 #include "otbSystem.h"
-#include "otbImageBase.h"
 
 #include "itkMetaDataObject.h"
 #include "itkPNGImageIO.h"
diff --git a/Code/IO/otbVectorImage.h b/Code/IO/otbVectorImage.h
index e0cb701997001a562a06607577af2bcae3a51d82..d004eac12ae8786b475517f0339dd2aff0dee98a 100644
--- a/Code/IO/otbVectorImage.h
+++ b/Code/IO/otbVectorImage.h
@@ -23,7 +23,7 @@
 #endif
 
 #include "itkVectorImage.h"
-#include "otbImageBase.h"
+#include "otbImageMetadataInterface.h"
 
 #include <string.h>
 
@@ -34,8 +34,7 @@ namespace otb
  *
  */
 template <class TPixel, unsigned int VImageDimension=2>
-class ITK_EXPORT VectorImage : public itk::VectorImage<TPixel, VImageDimension>,
-      public ImageBase
+class ITK_EXPORT VectorImage : public itk::VectorImage<TPixel, VImageDimension>
 {
 public:
 
@@ -46,8 +45,8 @@ public:
   typedef itk::SmartPointer<const Self>  ConstPointer;
   typedef itk::WeakPointer<const Self> ConstWeakPointer;
 
-  typedef ImageBase::VectorType  VectorType;
-  typedef ImageBase::ImageKeywordlistType  ImageKeywordlistType;
+  typedef ImageMetadataInterface::VectorType  VectorType;
+  typedef ImageMetadataInterface::ImageKeywordlistType  ImageKeywordlistType;
 
   /** Method for creation through the object factory. */
   itkNewMacro(Self);
@@ -182,7 +181,7 @@ protected:
 private:
   VectorImage(const Self&); //purposely not implemented
   void operator=(const Self&); //purposely not implemented
-
+  typename ImageMetadataInterface::Pointer m_ImageMetadataInterface;
 };
 
 
diff --git a/Code/IO/otbVectorImage.txx b/Code/IO/otbVectorImage.txx
index 6c1ee3837758d7a698dcd57b1e573cfa5794bf5c..a7ffe306210b80c146933309a0f9ceeed99b44e1 100644
--- a/Code/IO/otbVectorImage.txx
+++ b/Code/IO/otbVectorImage.txx
@@ -31,114 +31,115 @@ namespace otb
 template <class TPixel, unsigned int VImageDimension>
 VectorImage<TPixel,VImageDimension>::VectorImage()
 {
+  m_ImageMetadataInterface = ImageMetadataInterface::New();
 }
 
 template <class TPixel, unsigned int VImageDimension>
 std::string VectorImage<TPixel, VImageDimension>::GetProjectionRef( void ) const
 {
-  return ( ImageBase::GetProjectionRef( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetProjectionRef( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 std::string VectorImage<TPixel, VImageDimension>::GetGCPProjection( void ) const
 {
-  return ( ImageBase::GetGCPProjection( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetGCPProjection( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 unsigned int VectorImage<TPixel, VImageDimension>::GetGCPCount( void  ) const
 {
-  return ( ImageBase::GetGCPCount( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetGCPCount( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 OTB_GCP & VectorImage<TPixel, VImageDimension>::GetGCPs ( unsigned int GCPnum )
 {
-  return ( this->ImageBase::GetGCPs( this->GetMetaDataDictionary(), GCPnum ) );
+  return ( m_ImageMetadataInterface->GetGCPs( this->GetMetaDataDictionary(), GCPnum ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 std::string VectorImage<TPixel, VImageDimension>::GetGCPId( unsigned int GCPnum  ) const
 {
-  return ( ImageBase::GetGCPId( this->GetMetaDataDictionary(), GCPnum ) );
+  return ( m_ImageMetadataInterface->GetGCPId( this->GetMetaDataDictionary(), GCPnum ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 std::string VectorImage<TPixel, VImageDimension>::GetGCPInfo( unsigned int GCPnum ) const
 {
-  return ( ImageBase::GetGCPInfo( this->GetMetaDataDictionary(), GCPnum ) );
+  return ( m_ImageMetadataInterface->GetGCPInfo( this->GetMetaDataDictionary(), GCPnum ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 double VectorImage<TPixel, VImageDimension>::GetGCPRow( unsigned int GCPnum ) const
 {
-  return ( ImageBase::GetGCPRow( this->GetMetaDataDictionary(), GCPnum ) );
+  return ( m_ImageMetadataInterface->GetGCPRow( this->GetMetaDataDictionary(), GCPnum ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 double VectorImage<TPixel, VImageDimension>::GetGCPCol( unsigned int GCPnum ) const
 {
-  return ( ImageBase::GetGCPCol( this->GetMetaDataDictionary(), GCPnum ) );
+  return ( m_ImageMetadataInterface->GetGCPCol( this->GetMetaDataDictionary(), GCPnum ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 double VectorImage<TPixel, VImageDimension>::GetGCPX( unsigned int GCPnum ) const
 {
-  return ( ImageBase::GetGCPX( this->GetMetaDataDictionary(), GCPnum ) );
+  return ( m_ImageMetadataInterface->GetGCPX( this->GetMetaDataDictionary(), GCPnum ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 double VectorImage<TPixel, VImageDimension>::GetGCPY( unsigned int GCPnum ) const
 {
-  return ( ImageBase::GetGCPY( this->GetMetaDataDictionary(), GCPnum ) );
+  return ( m_ImageMetadataInterface->GetGCPY( this->GetMetaDataDictionary(), GCPnum ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
 double VectorImage<TPixel, VImageDimension>::GetGCPZ( unsigned int GCPnum ) const
 {
-  return ( ImageBase::GetGCPZ( this->GetMetaDataDictionary(), GCPnum ) );
+  return ( m_ImageMetadataInterface->GetGCPZ( this->GetMetaDataDictionary(), GCPnum ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
-ImageBase::VectorType VectorImage<TPixel, VImageDimension>::GetGeoTransform( void ) const
+ImageMetadataInterface::VectorType VectorImage<TPixel, VImageDimension>::GetGeoTransform( void ) const
 {
-  return ( ImageBase::GetGeoTransform( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetGeoTransform( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
-ImageBase::VectorType VectorImage<TPixel, VImageDimension>::GetUpperLeftCorner( void ) const
+ImageMetadataInterface::VectorType VectorImage<TPixel, VImageDimension>::GetUpperLeftCorner( void ) const
 {
-  return ( ImageBase::GetUpperLeftCorner( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetUpperLeftCorner( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
-ImageBase::VectorType VectorImage<TPixel, VImageDimension>::GetUpperRightCorner( void ) const
+ImageMetadataInterface::VectorType VectorImage<TPixel, VImageDimension>::GetUpperRightCorner( void ) const
 {
-  return ( ImageBase::GetUpperRightCorner( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetUpperRightCorner( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
-ImageBase::VectorType VectorImage<TPixel, VImageDimension>::GetLowerLeftCorner( void ) const
+ImageMetadataInterface::VectorType VectorImage<TPixel, VImageDimension>::GetLowerLeftCorner( void ) const
 {
-  return ( ImageBase::GetLowerLeftCorner( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetLowerLeftCorner( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
-ImageBase::VectorType VectorImage<TPixel, VImageDimension>::GetLowerRightCorner( void ) const
+ImageMetadataInterface::VectorType VectorImage<TPixel, VImageDimension>::GetLowerRightCorner( void ) const
 {
-  return ( ImageBase::GetLowerRightCorner( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetLowerRightCorner( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
-ImageBase::ImageKeywordlistType VectorImage<TPixel, VImageDimension>::GetImageKeywordlist( void )
+ImageMetadataInterface::ImageKeywordlistType VectorImage<TPixel, VImageDimension>::GetImageKeywordlist( void )
 {
-  return ( ImageBase::GetImageKeywordlist( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetImageKeywordlist( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
-const ImageBase::ImageKeywordlistType VectorImage<TPixel, VImageDimension>::GetImageKeywordlist( void ) const
+const ImageMetadataInterface::ImageKeywordlistType VectorImage<TPixel, VImageDimension>::GetImageKeywordlist( void ) const
 {
-  return ( ImageBase::GetImageKeywordlist( this->GetMetaDataDictionary() ) );
+  return ( m_ImageMetadataInterface->GetImageKeywordlist( this->GetMetaDataDictionary() ) );
 }
 
 template <class TPixel, unsigned int VImageDimension>
@@ -156,7 +157,7 @@ VectorImage<TPixel, VImageDimension>::PrintSelf(std::ostream& os, itk::Indent in
 {
   Superclass::PrintSelf(os,indent);
 
-  this->ImageBase::PrintSelf( os, indent, this->GetMetaDataDictionary() );
+  m_ImageMetadataInterface->PrintSelf( os, indent, this->GetMetaDataDictionary() );
 
 }
 
diff --git a/Code/MultiScale/otbWPCost.h b/Code/MultiScale/otbWPCost.h
new file mode 100644
index 0000000000000000000000000000000000000000..d60982cbc1371822fee7b6e6333034730392fdf2
--- /dev/null
+++ b/Code/MultiScale/otbWPCost.h
@@ -0,0 +1,84 @@
+/*=========================================================================
+
+  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) Institut Telecom / Telecom Bretagne. All rights reserved. 
+  See ITCopyright.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 __otbWPCost__h
+#define __otbWPCost__h
+
+#include "itkObject.h"
+#include "itkMacro.h"
+
+namespace otb {
+
+/** \class FullyDecomposedWaveletPacketCost
+ * \brief Cost evaluation to be used into the Wavelet Packet decomposition class.
+ *
+ * This class implements the criteria to perform fully decomposed wavelet packet. 
+ * It is based on the depth of the decomposition only...
+ *
+ * \sa WaveletPacketForwardTransform
+ */
+template < class TImage >
+class ITK_EXPORT FullyDecomposedWaveletPacketCost
+  : public itk::Object
+{
+public:
+  /** Standard typedefs */
+  typedef FullyDecomposedWaveletPacketCost Self;
+  typedef itk::Object Superclass;
+  typedef itk::SmartPointer<Self> Pointer;
+	typedef itk::SmartPointer<const Self> ConstPointer;
+
+	/** Type macro */
+	itkNewMacro(Self);
+
+	/** Creation through object factory macro */
+	itkTypeMacro(FullyDecomposedWaveletPacketCost,Object);
+
+  typedef TImage ImageType;
+
+  /** Acces to the data */
+  itkGetMacro(NumberOfAllowedDecompositions,unsigned int);
+  itkSetMacro(NumberOfAllowedDecompositions,unsigned int);
+
+  /** Evaluate the cost */
+  virtual bool Evaluate ( unsigned int decomposition, ImageType * image )
+  {
+    return ( decomposition < m_NumberOfAllowedDecompositions );
+  }
+
+protected:
+  FullyDecomposedWaveletPacketCost ()
+  {
+    m_NumberOfAllowedDecompositions = 1;
+  }
+  virtual ~FullyDecomposedWaveletPacketCost() {}
+
+  unsigned int m_NumberOfAllowedDecompositions;
+private:
+  FullyDecomposedWaveletPacketCost ( const Self & ); // not implemented
+  void operator= ( const Self & );
+
+}; // end of class 
+
+} // end of namespace otb
+
+#endif
+
+
diff --git a/Code/MultiScale/otbWaveletForwardTransform.h b/Code/MultiScale/otbWaveletForwardTransform.h
index 1feb5ab3499e722cc61823beb9f316d432f0706e..50011c717987b7517ed90ce68445869814cd6059 100644
--- a/Code/MultiScale/otbWaveletForwardTransform.h
+++ b/Code/MultiScale/otbWaveletForwardTransform.h
@@ -68,7 +68,7 @@ public:
 
   typedef TFilter FilterType;
   typedef typename FilterType::Pointer FilterPointerType;
-  typedef List< FilterType > FilterListType;
+  typedef ObjectList< FilterType > FilterListType;
   typedef typename FilterListType::Pointer FilterListPointerType;
 
   itkGetObjectMacro(FilterList,FilterListType);
diff --git a/Code/MultiScale/otbWaveletPacketForwardTransform.h b/Code/MultiScale/otbWaveletPacketForwardTransform.h
new file mode 100644
index 0000000000000000000000000000000000000000..e73818d8fbf98d8a784d889deb7fbb855251b38d
--- /dev/null
+++ b/Code/MultiScale/otbWaveletPacketForwardTransform.h
@@ -0,0 +1,123 @@
+/*=========================================================================
+
+  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) Institut Telecom / Telecom Bretagne. All rights reserved. 
+  See ITCopyright.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 __otbWaveletPacketForwardTransform__h
+#define __otbWaveletPacketForwardTransform__h
+
+#include "otbImageToImageListFilter.h"
+
+namespace otb {
+
+/** \class WaveletPacketForwardTransform
+ * \brief Wavelet packet transformation framework
+ *
+ * This class defines the Wavelet Packet transformation of an image
+ * by using a (templated) elementary wavelet transform and (templated)
+ * cost criteria to stop the decomposition.
+ *
+ * It yields a list of images...
+ *
+ * TODO: this class implements the decomposition only. But there is no storage 
+ * of the way of decomposition. It should be added before implementing backward
+ * transformation.
+ *
+ * the user is supposed to initialize Cost properly (through GetCost() macro) 
+ * depending on its type before calling an Update(). The Cost class has to contain 
+ * a New() and Evaluate() function.
+ *
+ * \sa FullyDecomposedWaveletPacketCost
+ * \sa StationaryFilterBank
+ * \sa WaveletForwardTransform
+ *
+ */
+template < class TInputImage, class TOutputImage, class TFilter, class TCost >
+class ITK_EXPORT WaveletPacketForwardTransform
+  : public ImageToImageListFilter< TInputImage, TOutputImage >
+{
+public:
+  /** Standard typedefs */
+  typedef WaveletPacketForwardTransform Self;
+  typedef ImageToImageListFilter< TInputImage, TOutputImage > Superclass;
+  typedef itk::SmartPointer<Self> Pointer;
+	typedef itk::SmartPointer<const Self> ConstPointer;
+
+	/** Type macro */
+	itkNewMacro(Self);
+
+	/** Creation through object factory macro */
+	itkTypeMacro(WaveletPacketForwardTransform,ImageToImageListFilter);
+
+  typedef          TInputImage                  InputImageType;
+  typedef typename InputImageType::Pointer      InputImagePointerType;
+  typedef typename InputImageType::ConstPointer InputImageConstPointerType;
+  typedef typename InputImageType::RegionType   InputImageRegionType;
+  typedef typename InputImageType::PixelType    InputImagePixelType;
+  typedef typename InputImageType::SizeType     SizeType;
+  typedef typename InputImageType::ValueType    ValueType;
+
+  typedef typename Superclass::OutputImageType        OutputImageType;
+  typedef typename Superclass::OutputImagePointerType OutputImagePointerType;
+  typedef typename Superclass::OutputImageListType    OutputImageListType;
+  typedef typename OutputImageListType::Pointer       OutputImageListPointerType;
+  typedef typename OutputImageListType::Iterator      OutputImageIterator;
+
+  typedef TFilter                           FilterType;
+  typedef typename FilterType::Pointer      FilterPointerType;
+  typedef ObjectList< FilterType >          FilterListType;
+  typedef typename FilterListType::Pointer  FilterListPointerType;
+  typedef typename FilterListType::Iterator FilterListIterator;
+
+  itkGetObjectMacro(FilterList,FilterListType);
+
+  typedef TCost CostType;
+  typedef typename CostType::Pointer CostPointerType;
+
+  itkGetObjectMacro(Cost,CostType);
+ 
+  itkStaticConstMacro(InputImageDimension, unsigned int,TInputImage::ImageDimension);
+
+protected:
+  WaveletPacketForwardTransform();
+  virtual ~WaveletPacketForwardTransform() {}
+
+  /** Generate data redefinition */
+  virtual void GenerateData ();
+
+  /** Performs (if any) the local decomposition (called recursively) */
+  virtual void PerformDecomposition ( unsigned int depth, OutputImageIterator & outputIt );
+
+private:
+  WaveletPacketForwardTransform ( const Self &);
+  void operator= ( const Self & );
+
+  FilterListPointerType m_FilterList;
+  CostPointerType m_Cost;
+
+}; // end of class
+
+} // end of namespace
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbWaveletPacketForwardTransform.txx"
+#endif
+
+#endif
+
+
diff --git a/Code/MultiScale/otbWaveletPacketForwardTransform.txx b/Code/MultiScale/otbWaveletPacketForwardTransform.txx
new file mode 100644
index 0000000000000000000000000000000000000000..062372d9e21b1b691ebcff20fe3a511011970d53
--- /dev/null
+++ b/Code/MultiScale/otbWaveletPacketForwardTransform.txx
@@ -0,0 +1,95 @@
+/*=========================================================================
+
+  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) Institut Telecom / Telecom Bretagne. All rights reserved. 
+  See ITCopyright.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 __otbWaveletPacketForwardTransform__txx
+#define __otbWaveletPacketForwardTransform__txx
+
+#include "otbWaveletPacketForwardTransform.h"
+
+namespace otb {
+
+template < class TInputImage, class TOutputImage, class TFilter, class TCost >
+WaveletPacketForwardTransform< TInputImage, TOutputImage, TFilter, TCost >
+::WaveletPacketForwardTransform ()
+{
+  this->SetNumberOfRequiredInputs(1);
+  this->SetNumberOfOutputs(1);
+
+  m_FilterList = FilterListType::New();
+  m_Cost = CostType::New();
+}
+
+template < class TInputImage, class TOutputImage, class TFilter, class TCost >
+void
+WaveletPacketForwardTransform< TInputImage, TOutputImage, TFilter, TCost >
+::GenerateData ()
+{
+  //this->AllocateOutputs();
+  
+  /*
+   * Start with a decomposition
+   */
+  this->GetFilterList()->PushBack( FilterType::New() );
+  FilterType * filter = this->GetFilterList()->GetNthElement( 0 );
+  filter->SetInput( this->GetInput() );
+  filter->Update();
+  for ( unsigned int idx = 0; idx < filter->GetNumberOfOutputs(); idx++ )
+    this->GetOutput()->PushBack( filter->GetOutput( idx ) );
+
+  /**
+   * Loop decompositions
+   */
+  for ( OutputImageIterator outputIt = this->GetOutput()->Begin();
+      outputIt != this->GetOutput()->End();
+      ++outputIt ) 
+  {
+    PerformDecomposition( 1, outputIt );
+  }
+}
+
+template < class TInputImage, class TOutputImage, class TFilter, class TCost >
+void
+WaveletPacketForwardTransform< TInputImage, TOutputImage, TFilter, TCost >
+::PerformDecomposition ( unsigned int depth, OutputImageIterator & outputIt )
+{
+  if ( this->GetCost()->Evaluate( depth, outputIt.Get() ) )
+  {
+    this->GetFilterList()->PushBack( FilterType::New() );
+    FilterType * filter = this->GetFilterList()->GetNthElement( this->GetFilterList()->Size()-1 );
+    filter->SetInput( outputIt.Get() );
+    filter->Update();
+
+    depth++;
+    outputIt.Set( filter->GetOutput(0) );
+    PerformDecomposition( depth, outputIt );
+
+    for ( unsigned int idx = 1; idx < filter->GetNumberOfOutputs(); idx++ )
+    {
+      ++outputIt;
+      outputIt = this->GetOutput()->Insert( outputIt, filter->GetOutput( idx ) );
+      PerformDecomposition( depth, outputIt );
+    }
+  }
+}
+
+} // end of namespace otb
+
+#endif
+
diff --git a/Code/Radiometry/otbMultiChannelGAndRIndexImageFilter.h b/Code/Radiometry/otbMultiChannelGAndRIndexImageFilter.h
index 5ec7f433680e96c391df1f9eafe7e4597109307c..adc1734278cd19575c1d75d218df63e8110cd5be 100644
--- a/Code/Radiometry/otbMultiChannelGAndRIndexImageFilter.h
+++ b/Code/Radiometry/otbMultiChannelGAndRIndexImageFilter.h
@@ -73,7 +73,9 @@ protected:
   /// Before generating data, set functor parameters
   virtual void BeforeThreadedGenerateData()
   {
-    if(m_GreenIndex < 1 || m_RedIndex < 1)
+    unsigned int lNbChan = this->GetInput()->GetNumberOfComponentsPerPixel();
+    if(m_GreenIndex < 1 || m_RedIndex < 1 ||
+       m_GreenIndex > lNbChan || m_RedIndex > lNbChan )
       {
       itkExceptionMacro(<<"Channel indices must belong to range [1, ...[");
       }
diff --git a/Code/Radiometry/otbMultiChannelRAndBAndNIRIndexImageFilter.h b/Code/Radiometry/otbMultiChannelRAndBAndNIRIndexImageFilter.h
index 0eb437ecc2a6a74f1549889a4cfbcc1cba07c4e4..800809ebe4178e5e359c2e3f7c909ecfeea4e059 100644
--- a/Code/Radiometry/otbMultiChannelRAndBAndNIRIndexImageFilter.h
+++ b/Code/Radiometry/otbMultiChannelRAndBAndNIRIndexImageFilter.h
@@ -77,7 +77,9 @@ protected:
   /// Before generating data, set functor parameters
   virtual void BeforeThreadedGenerateData()
   {
-    if(m_RedIndex < 1 || m_BlueIndex < 1 || m_NIRIndex < 1)
+    unsigned int lNbChan = this->GetInput()->GetNumberOfComponentsPerPixel();
+    if(m_RedIndex < 1 || m_BlueIndex < 1 || m_NIRIndex < 1 ||
+       m_RedIndex > lNbChan || m_BlueIndex > lNbChan || m_NIRIndex > lNbChan)
       {
       itkExceptionMacro(<<"Channel indices must belong to range [1, ...[");
       }
diff --git a/Code/Radiometry/otbMultiChannelRAndBAndNIRVegetationIndexImageFilter.h b/Code/Radiometry/otbMultiChannelRAndBAndNIRVegetationIndexImageFilter.h
index 15a776e337dfedbfdda5edd7410ae12f870618c5..63e9659a921e0d13feaa14a847ae3dbb742df3fa 100644
--- a/Code/Radiometry/otbMultiChannelRAndBAndNIRVegetationIndexImageFilter.h
+++ b/Code/Radiometry/otbMultiChannelRAndBAndNIRVegetationIndexImageFilter.h
@@ -77,7 +77,9 @@ protected:
   /// Before generating data, set functor parameters
   virtual void BeforeThreadedGenerateData()
   {
-    if(m_RedIndex < 1 || m_BlueIndex < 1 || m_NIRIndex < 1)
+    unsigned int lNbChan = this->GetInput()->GetNumberOfComponentsPerPixel();
+    if(m_RedIndex < 1 || m_BlueIndex < 1 || m_NIRIndex < 1 ||
+       m_RedIndex > lNbChan || m_BlueIndex > lNbChan || m_NIRIndex > lNbChan)
       {
       itkExceptionMacro(<<"Channel indices must belong to range [1, ...[");
       }
diff --git a/Code/Radiometry/otbMultiChannelRAndGAndNIRIndexImageFilter.h b/Code/Radiometry/otbMultiChannelRAndGAndNIRIndexImageFilter.h
index b20d5201cec3b6fdaf8bf221e2b47f2b4ae96af8..6569e5f1dfe9ace786ae6134d3ae56e6a3ae5909 100644
--- a/Code/Radiometry/otbMultiChannelRAndGAndNIRIndexImageFilter.h
+++ b/Code/Radiometry/otbMultiChannelRAndGAndNIRIndexImageFilter.h
@@ -77,7 +77,9 @@ protected:
   /// Before generating data, set functor parameters
   virtual void BeforeThreadedGenerateData()
   {
-    if(m_RedIndex < 1 || m_GreenIndex < 1 || m_NIRIndex < 1)
+    unsigned int lNbChan = this->GetInput()->GetNumberOfComponentsPerPixel();
+    if(m_RedIndex < 1 || m_GreenIndex < 1 || m_NIRIndex < 1 ||
+       m_RedIndex > lNbChan || m_GreenIndex > lNbChan || m_NIRIndex > lNbChan)
       {
       itkExceptionMacro(<<"Channel indices must belong to range [1, ...[");
       }
diff --git a/Code/Radiometry/otbMultiChannelRAndGAndNIRVegetationIndexImageFilter.h b/Code/Radiometry/otbMultiChannelRAndGAndNIRVegetationIndexImageFilter.h
index 890436e829d9588e10d15055b9e4674ec20b871f..c2a2e499619a31f6cc7821f9acf11808b91cc2ea 100644
--- a/Code/Radiometry/otbMultiChannelRAndGAndNIRVegetationIndexImageFilter.h
+++ b/Code/Radiometry/otbMultiChannelRAndGAndNIRVegetationIndexImageFilter.h
@@ -77,7 +77,9 @@ protected:
   /// Before generating data, set functor parameters
   virtual void BeforeThreadedGenerateData()
   {
-    if(m_RedIndex < 1 || m_GreenIndex < 1 || m_NIRIndex < 1)
+    unsigned int lNbChan = this->GetInput()->GetNumberOfComponentsPerPixel();
+    if(m_RedIndex < 1 || m_GreenIndex < 1 || m_NIRIndex < 1 ||
+       m_RedIndex > lNbChan || m_GreenIndex > lNbChan || m_NIRIndex > lNbChan)
       {
       itkExceptionMacro(<<"Channel indices must belong to range [1, ...[");
       }
diff --git a/Code/Radiometry/otbMultiChannelRAndNIRIndexImageFilter.h b/Code/Radiometry/otbMultiChannelRAndNIRIndexImageFilter.h
index c58a23c5a453c08abdbdfd29ac00ce05cc3776f5..e681331466c077bdcd22aa30cd76018203f8bbb4 100644
--- a/Code/Radiometry/otbMultiChannelRAndNIRIndexImageFilter.h
+++ b/Code/Radiometry/otbMultiChannelRAndNIRIndexImageFilter.h
@@ -73,7 +73,9 @@ protected:
   /// Before generating data, set functor parameters
   virtual void BeforeThreadedGenerateData()
   {
-    if(m_RedIndex < 1 || m_NIRIndex < 1)
+    unsigned int lNbChan = this->GetInput()->GetNumberOfComponentsPerPixel();
+    if(m_RedIndex < 1 || m_NIRIndex < 1 ||
+       m_RedIndex > lNbChan || m_NIRIndex > lNbChan)
       {
       itkExceptionMacro(<<"Channel indices must belong to range [1, ...[");
       }
diff --git a/Code/Radiometry/otbMultiChannelRAndNIRVegetationIndexImageFilter.h b/Code/Radiometry/otbMultiChannelRAndNIRVegetationIndexImageFilter.h
index 093e0290d8852146e6e37dc4fd87ef559258ecba..07c7a76dd98e99430c53ed2044487e3dbed6f808 100644
--- a/Code/Radiometry/otbMultiChannelRAndNIRVegetationIndexImageFilter.h
+++ b/Code/Radiometry/otbMultiChannelRAndNIRVegetationIndexImageFilter.h
@@ -74,7 +74,9 @@ protected:
   /// Before generating data, set functor parameters
   virtual void BeforeThreadedGenerateData()
   {
-    if(m_RedIndex < 1 || m_NIRIndex < 1)
+    unsigned int lNbChan = this->GetInput()->GetNumberOfComponentsPerPixel();
+    if(m_RedIndex < 1 || m_NIRIndex < 1 ||
+       m_RedIndex > lNbChan || m_NIRIndex > lNbChan)
       {
       itkExceptionMacro(<<"Channel indices must belong to range [1, ...[");
       }
diff --git a/Code/Radiometry/otbWaterIndicesFunctor.h b/Code/Radiometry/otbWaterIndicesFunctor.h
new file mode 100644
index 0000000000000000000000000000000000000000..b3e9431c3c643267fadc4e1ace4646b33a9403ee
--- /dev/null
+++ b/Code/Radiometry/otbWaterIndicesFunctor.h
@@ -0,0 +1,429 @@
+/*=========================================================================
+
+  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 __otbWaterIndicesFunctor_h
+#define __otbWaterIndicesFunctor_h
+
+#include "otbMath.h"
+#include "itkVariableLengthVector.h"
+
+namespace otb
+{
+namespace Functor
+{
+/**
+   * \class WaterIndexBase
+   * \brief Base class
+   *
+   *  Implement operators for UnaryFunctorImageFilter templated with a
+   *  VectorImage and BinaryFunctorImageFilter templated with single
+   *  images.
+   *  Subclasses should NOT overload operators, they must  re-implement
+   *  the Evaluate() method.
+   *
+   * \ingroup Radiometry
+ */
+template<class TInput1, class TInput2, class TOutput>
+class WaterIndexBase
+{
+public:
+  /// Vector pixel type used to support both vector images and multiple
+  /// input images
+  typedef itk::VariableLengthVector<TInput1> InputVectorType;
+
+  // Operator on vector pixel type
+  inline TOutput operator()(const InputVectorType & inputVector)
+  {
+    return this->Evaluate(inputVector[m_Index1-1],static_cast<TInput2>(inputVector[m_Index2-1]));
+  }
+
+  // Binary operator
+  inline TOutput operator()(const TInput1 &id1, const TInput2 &id2)
+  {
+    return this->Evaluate(id1,id2);
+  };
+  /// Constructor
+  WaterIndexBase() {};
+  /// Desctructor
+  virtual ~WaterIndexBase() {};
+
+  /// Set Index 1
+  void SetIndex1(unsigned int channel)
+  {
+    m_Index1 = channel;
+  }
+  /// Get Index 1
+  unsigned int GetIndex1()
+  {
+    return m_Index1;
+  }
+  /// Set Index 2
+  void SetIndex2(unsigned int channel)
+  {
+    m_Index2 = channel;
+  }
+  /// Get Index 2
+  unsigned int GetIndex2()
+  {
+    return m_Index2;
+  }
+protected:
+  // This method must be reimplemented in subclasses to actually
+  // compute the index value
+  virtual TOutput Evaluate(const TInput1 & id1, const TInput2 & id2) const = 0;
+
+private:
+  unsigned int m_Index1;
+  unsigned int m_Index2;
+};
+
+
+/** \class WaterIndexFunctor
+ *  \brief This functor will be used for most of water index functors.
+ *
+ *  \ingroup Functor
+ * \ingroup Radiometry
+ */
+template <class TInput1, class TInput2, class TOutput>
+class WaterIndexFunctor : public WaterIndexBase<TInput1, TInput2, TOutput>
+{
+public:
+  WaterIndexFunctor() {};
+  ~WaterIndexFunctor() {};
+protected:
+  inline TOutput Evaluate(const TInput1 &id1, const TInput2 &id2) const
+  {
+    double dindex1 = static_cast<double>(id1);
+    double dindex2 = static_cast<double>(id2);
+    double ddenom = dindex1 + dindex2;
+    if ( ddenom == 0 )
+      {
+      return static_cast<TOutput>(0.);
+      }
+    return ( static_cast<TOutput>((dindex1- dindex2)/ddenom));
+  }
+};
+
+
+
+
+
+
+
+
+
+/** \class SRWI
+ *  \brief This functor computes the Simple Ratio Water Index (SRWI)
+ *  \brief For MODIS bands 860 & 1240
+ *
+ *  \ingroup Functor
+ * \ingroup Radiometry
+ */
+template <class TInput1, class TInput2, class TOutput>
+class SRWI : public WaterIndexBase<TInput1, TInput2, TOutput>
+{
+public:
+  SRWI() {};
+  ~SRWI() {};
+protected:
+  inline TOutput Evaluate(const TInput1 &rho860, const TInput2 &rho1240) const
+  {
+    double drho860 = static_cast<double>(rho860);
+    double drho1240 = static_cast<double>(rho1240);
+    if ( drho1240 == 0 )
+      {
+      return static_cast<TOutput>(0.);
+      }
+    return ( static_cast<TOutput>(drho860/drho1240) );
+  }
+};
+
+
+/** \class NDWI
+ *  \brief This functor computes the Normalized Difference Water Index (NDWI)
+ *  \brief Also called :
+ *  \brief NDII (Normalized Difference Infrared Index)
+ *  \brief LSWI (Land Surface Water Index)
+ *  \brief NDWI (Normalized Difference Moisture Index)
+ *
+ *  [Gao, 1996 ]
+ *
+ *  \ingroup Functor
+ * \ingroup Radiometry
+ */
+template <class TInput1, class TInput2, class TOutput>
+class NDWI : public WaterIndexBase<TInput1,TInput2,TOutput>
+{
+public:
+  typedef WaterIndexFunctor<TInput1, TInput2, TOutput> WIFunctorType;
+  /// Constructor
+  NDWI() {};
+  /// Desctructor
+  ~NDWI() {};
+  WIFunctorType GetWIFunctor(void)const
+  {
+    return (m_WIFunctor);
+  }
+  /// Set Index NIR
+  void SetNIRIndex(unsigned int channel)
+  {
+    this->SetIndex1(channel);
+  }
+  /// Get Index NIR
+  unsigned int GetNIRIndex()
+  {
+    return this->GetIndex1;
+  }
+  /// Set Index MIR
+  void SetMIRIndex(unsigned int channel)
+  {
+    this->SetIndex2(channel);
+  }
+  /// Get Index MIR
+  unsigned int GetMIRIndex()
+  {
+    return this->GetIndex2;
+  }
+
+protected:
+  inline TOutput Evaluate(const TInput1 &nir, const TInput2 &mir) const
+  {
+    return ( static_cast<TOutput>(GetWIFunctor()(nir,mir)) );
+  }
+private:
+  // Water Index Classic Functor
+  const WIFunctorType m_WIFunctor;
+};
+
+
+/** \class NDWI2
+ *  \brief This functor computes the Normalized Difference Water Index (NDWI2)
+ *
+ *  [McFeeters, 1996 ]
+ *
+ *  \ingroup Functor
+ * \ingroup Radiometry
+ */
+template <class TInput1, class TInput2, class TOutput>
+class NDWI2 : public WaterIndexBase<TInput1,TInput2,TOutput>
+{
+public:
+  typedef WaterIndexFunctor<TInput1, TInput2, TOutput> WIFunctorType;
+  /// Constructor
+  NDWI2() {};
+  /// Desctructor
+  ~NDWI2() {};
+  WIFunctorType GetWIFunctor(void)const
+  {
+    return (m_WIFunctor);
+  }
+  /// Set Index G
+  void SetGIndex(unsigned int channel)
+  {
+    this->SetIndex1(channel);
+  }
+  /// Get Index G
+  unsigned int GetGIndex()
+  {
+    return this->GetIndex1;
+  }
+  /// Set Index NIR
+  void SetNIRIndex(unsigned int channel)
+  {
+    this->SetIndex2(channel);
+  }
+  /// Get Index NIR
+  unsigned int GetNIRIndex()
+  {
+    return this->GetIndex2;
+  }
+
+protected:
+  inline TOutput Evaluate(const TInput1 &g, const TInput2 &nir) const
+  {
+    return ( static_cast<TOutput>(GetWIFunctor()(g,nir)) );
+  }
+private:
+  // Water Index Classic Functor
+  const WIFunctorType m_WIFunctor;
+};
+
+/** \class MNDWI
+ *  \brief This functor computes the Modified Normalized Difference Water Index (MNDWI)
+ *
+ *  [Xu & al., 2006 ]
+ *
+ *  \ingroup Functor
+ * \ingroup Radiometry
+ */
+template <class TInput1, class TInput2, class TOutput>
+class MNDWI : public WaterIndexBase<TInput1,TInput2,TOutput>
+{
+public:
+  typedef WaterIndexFunctor<TInput1, TInput2, TOutput> WIFunctorType;
+  /// Constructor
+  MNDWI() {};
+  /// Desctructor
+  ~MNDWI() {};
+  WIFunctorType GetWIFunctor(void)const
+  {
+    return (m_WIFunctor);
+  }
+  /// Set Index G
+  void SetGIndex(unsigned int channel)
+  {
+    this->SetIndex1(channel);
+  }
+  /// Get Index G
+  unsigned int GetGIndex()
+  {
+    return this->GetIndex1;
+  }
+  /// Set Index MIR
+  void SetMIRIndex(unsigned int channel)
+  {
+    this->SetIndex2(channel);
+  }
+  /// Get Index MIR
+  unsigned int GetMIRIndex()
+  {
+    return this->GetIndex2;
+  }
+
+protected:
+  inline TOutput Evaluate(const TInput1 &g, const TInput2 &mir) const
+  {
+    return ( static_cast<TOutput>(GetWIFunctor()(g,mir)) );
+  }
+private:
+  // Water Index Classic Functor
+  const WIFunctorType m_WIFunctor;
+};
+
+
+/** \class NDPI
+ *  \brief This functor computes the Normalized Difference Pond Index (NDPI)
+ *
+ *  [J.P Lacaux & al., 2006 ]
+ *
+ *  \ingroup Functor
+ * \ingroup Radiometry
+ */
+template <class TInput1, class TInput2, class TOutput>
+class NDPI : public WaterIndexBase<TInput1,TInput2,TOutput>
+{
+public:
+  typedef WaterIndexFunctor<TInput1, TInput2, TOutput> WIFunctorType;
+  /// Constructor
+  NDPI() {};
+  /// Desctructor
+  ~NDPI() {};
+  WIFunctorType GetWIFunctor(void)const
+  {
+    return (m_WIFunctor);
+  }
+  /// Set Index MIR
+  void SetMIRIndex(unsigned int channel)
+  {
+    this->SetIndex1(channel);
+  }
+  /// Get Index MIR
+  unsigned int GetMIRIndex()
+  {
+    return this->GetIndex1;
+  }
+ /// Set Index G
+  void SetGIndex(unsigned int channel)
+  {
+    this->SetIndex2(channel);
+  }
+  /// Get Index G
+  unsigned int GetGIndex()
+  {
+    return this->GetIndex2;
+  }
+
+protected:
+  inline TOutput Evaluate(const TInput1 &mir, const TInput2 &g) const
+  {
+    return ( static_cast<TOutput>(GetWIFunctor()(mir,g)) );
+  }
+private:
+  // Water Index Classic Functor
+  const WIFunctorType m_WIFunctor;
+};
+
+
+/** \class NDTI
+ *  \brief This functor computes the Normalized Difference Turbidity Index (NDTI)
+ *
+ *  [J.P Lacaux & al., 2006 ]
+ *
+ *  \ingroup Functor
+ * \ingroup Radiometry
+ */
+template <class TInput1, class TInput2, class TOutput>
+class NDTI : public WaterIndexBase<TInput1,TInput2,TOutput>
+{
+public:
+  typedef WaterIndexFunctor<TInput1, TInput2, TOutput> WIFunctorType;
+  /// Constructor
+  NDTI() {};
+  /// Desctructor
+  ~NDTI() {};
+  WIFunctorType GetWIFunctor(void)const
+  {
+    return (m_WIFunctor);
+  }
+  /// Set Index R
+  void SetRIndex(unsigned int channel)
+  {
+    this->SetIndex1(channel);
+  }
+  /// Get Index R
+  unsigned int GetRIndex()
+  {
+    return this->GetIndex1;
+  }
+ /// Set Index G
+  void SetGIndex(unsigned int channel)
+  {
+    this->SetIndex2(channel);
+  }
+  /// Get Index G
+  unsigned int GetGIndex()
+  {
+    return this->GetIndex2;
+  }
+
+protected:
+  inline TOutput Evaluate(const TInput1 &r, const TInput2 &g) const
+  {
+    return ( static_cast<TOutput>(GetWIFunctor()(r,g)) );
+  }
+private:
+  // Water Index Classic Functor
+  const WIFunctorType m_WIFunctor;
+};
+
+
+
+} // namespace Functor
+} // namespace otb
+
+#endif
+
diff --git a/Examples/BasicFilters/CMakeLists.txt b/Examples/BasicFilters/CMakeLists.txt
index 0177c1ba5fcf2b14e621e8cd990119e0899f6b48..76d308c7cb30ee1be22c21f55c82cfbeee1f44f8 100644
--- a/Examples/BasicFilters/CMakeLists.txt
+++ b/Examples/BasicFilters/CMakeLists.txt
@@ -4,6 +4,9 @@ INCLUDE_REGULAR_EXPRESSION("^.*$")
 ADD_EXECUTABLE(LeeImageFilter LeeImageFilter.cxx )
 TARGET_LINK_LIBRARIES(LeeImageFilter OTBCommon OTBIO)
 
+ADD_EXECUTABLE(FrostImageFilter FrostImageFilter.cxx )
+TARGET_LINK_LIBRARIES(FrostImageFilter OTBCommon OTBIO)
+
 ADD_EXECUTABLE(DEMToRainbowExample DEMToRainbowExample.cxx )
 TARGET_LINK_LIBRARIES(DEMToRainbowExample OTBCommon OTBIO)
 
@@ -44,6 +47,18 @@ ADD_TEST(LeeImageFilterTest ${EXE_TESTS}
     3 1
 )
 
+# ------- FrostImageFilterTest----------
+
+ADD_TEST(FrostImageFilterTest ${EXE_TESTS}
+	--compare-n-images ${TOL} 1
+	${BASELINE}/GomaSmallFrostFiltered.png
+	${TEMP}/GomaSmallFrostFiltered.png
+	FrostImageFilterTest
+	${INPUTDATA}/GomaSmall.png
+	${TEMP}/GomaSmallFrostFiltered.png
+    5 0.1
+)
+
 # ------- DEMToRainbowExampleTest ----------
 
 ADD_TEST(DEMToRainbowExampleTest ${EXE_TESTS}
diff --git a/Examples/BasicFilters/FrostImageFilter.cxx b/Examples/BasicFilters/FrostImageFilter.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..d442ca3a10420f91fe2df773056f51282e01071f
--- /dev/null
+++ b/Examples/BasicFilters/FrostImageFilter.cxx
@@ -0,0 +1,189 @@
+/*=========================================================================
+
+  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
+
+#ifdef __BORLANDC__
+#define ITK_LEAN_AND_MEAN
+#endif
+
+//  Software Guide : BeginCommandLineArgs
+//    INPUTS: {GomaSmall.png}
+//    OUTPUTS: {GomaSmallFrostFiltered.png}
+//    5 0.1
+//  Software Guide : EndCommandLineArgs
+
+// Software Guide : BeginLatex
+//
+// This example illustrates the use of the \doxygen{otb}{FrostImageFilter}.
+// This filter belongs to the family of the edge-preserving smoothing
+// filters which are usually used for speckle reduction in radar
+// images.
+//
+// This filter uses a negative exponential convolution kernel.
+// The output of the filter for pixel p is:
+//      $ \hat I_{s}=\sum_{p\in\eta_{p}} m_{p}I_{p} $
+//
+// where :   $ m_{p}=\frac{KC_{s}^{2}\exp(-KC_{s}^{2}d_{s,p})}{\sum_{p\in\eta_{p}} KC_{s}^{2}\exp(-KC_{s}^{2}d_{s,p})} $
+//    and  $ d_{s,p}=\sqrt{(i-i_{p})^2+(j-j_{p})^2} $
+//
+// \begin{itemize}
+// \item $ K $     : the decrease coefficient
+// \item $ (i,j)$ : the coordinates of the pixel inside the region
+// defined by $ \eta_{s} $
+// \item $ (i_{p},j_{p})$ : the coordinates of the pixels belonging to $ \eta_{p} \subset \eta_{s} $
+// \item $ C_{s}$ : the variation coefficient computed over $ \eta_{p}$
+// \end{itemize}
+//
+//
+//
+// Most of this example is similar to the previous one and only the differences
+// will be highlighted.
+//
+// First, we need to include the header:
+//
+// Software Guide : EndLatex
+
+// Software Guide : BeginCodeSnippet
+#include "otbFrostImageFilter.h"
+// Software Guide : EndCodeSnippet
+
+#include "otbImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+
+int main( int argc, char * argv[] )
+{
+
+  if ( argc != 5 )
+  {
+    std::cerr << "Usage: " << argv[0] << " inputImageFile ";
+    std::cerr << " outputImageFile radius deramp" << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  typedef  unsigned char  PixelType;
+
+  typedef otb::Image< PixelType,  2 >   InputImageType;
+  typedef otb::Image< PixelType,  2 >   OutputImageType;
+
+
+  //  Software Guide : BeginLatex
+  //
+  //  The filter can be instantiated using the image types defined previously.
+  //
+  //  Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  typedef otb::FrostImageFilter< InputImageType, OutputImageType >  FilterType;
+  // Software Guide : EndCodeSnippet
+
+
+
+  typedef otb::ImageFileReader< InputImageType >  ReaderType;
+
+  typedef otb::ImageFileWriter< OutputImageType >  WriterType;
+
+  ReaderType::Pointer reader = ReaderType::New();
+  FilterType::Pointer filter = FilterType::New();
+
+
+  WriterType::Pointer writer = WriterType::New();
+  writer->SetInput( filter->GetOutput() );
+  reader->SetFileName( argv[1] );
+
+
+  //  Software Guide : BeginLatex
+  //
+  //  The image obtained with the reader is passed as input to the
+  //  \doxygen{otb}{FrostImageFilter}.
+  //
+  //  \index{otb::FrostImageFilter!SetInput()}
+  //  \index{otb::FileImageReader!GetOutput()}
+  //
+  //  Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  filter->SetInput( reader->GetOutput() );
+  // Software Guide : EndCodeSnippet
+
+
+  //  Software Guide : BeginLatex
+  //
+  //  The method \code{SetRadius()} defines the size of the window to
+  //  be used for the computation of the local statistics. The method
+  //  \code{SetDeramp()} sets the $K$ coefficient.
+  //
+  //  \index{otb::FrostImageFilter!SetRadius()}
+  //  \index{otb::FrostImageFilter!SetDeramp()}
+  //  \index{SetDeramp()!otb::FrostImageFilter}
+  //
+  //  Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  FilterType::SizeType Radius;
+  Radius[0]= atoi(argv[3]);
+  Radius[1]= atoi(argv[3]);
+
+  filter->SetRadius( Radius );
+  filter->SetDeramp( atof(argv[4]) );
+  // Software Guide : EndCodeSnippet
+
+
+  //  Software Guide : BeginLatex
+  //
+  //  The filter is executed by invoking the \code{Update()} method. If the
+  //  filter is part of a larger image processing pipeline, calling
+  //  \code{Update()} on a downstream filter will also trigger update of this
+  //  filter.
+  //
+  //  Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  filter->Update();
+  // Software Guide : EndCodeSnippet
+
+
+  writer->SetFileName( argv[2] );
+  writer->Update();
+
+
+  //  Software Guide : BeginLatex
+  // Figure~\ref{fig:FROST_FILTER} shows the result of applying the Frost
+  // filter to a SAR image.
+  // \begin{figure}
+  // \center
+  // \includegraphics[width=0.44\textwidth]{GomaSmall.eps}
+  // \includegraphics[width=0.44\textwidth]{GomaSmallFrostFiltered.eps}
+  // \itkcaption[Frost Filter Application]{Result of applying the
+  // \doxygen{otb}{FrostImageFilter} to a SAR image.}
+  // \label{fig:FROST_FILTER}
+  // \end{figure}
+  //
+  //  \relatedClasses
+  //  \begin{itemize}
+  //  \item \doxygen{otb}{LeeImageFilter}
+  //  \end{itemize}
+  //
+  //  Software Guide : EndLatex
+
+
+  return EXIT_SUCCESS;
+}
+
diff --git a/Examples/BasicFilters/LeeImageFilter.cxx b/Examples/BasicFilters/LeeImageFilter.cxx
index f31d848ef9eeec42825568a453aaa9a7775b4da6..0fb91ca37ed9277c18178374d68c81e8c7b5190b 100644
--- a/Examples/BasicFilters/LeeImageFilter.cxx
+++ b/Examples/BasicFilters/LeeImageFilter.cxx
@@ -158,7 +158,7 @@ int main( int argc, char * argv[] )
   //  image.
   //
   //  \index{otb::LeeImageFilter!SetRadius()}
-  //  \index{otb::LeeImageFilter!NbLooks()}
+  //  \index{otb::LeeImageFilter!SetNbLooks()}
   //  \index{SetNbLooks()!otb::LeeImageFilter}
   //
   //  Software Guide : EndLatex
diff --git a/Examples/BasicFilters/otbBasicFiltersExamplesTests.cxx b/Examples/BasicFilters/otbBasicFiltersExamplesTests.cxx
index 27e3e623f9c7790bb2a0b6357cf6b6b19e282181..858c43c32735c37ba7cb28d59c3d51e733ee7785 100644
--- a/Examples/BasicFilters/otbBasicFiltersExamplesTests.cxx
+++ b/Examples/BasicFilters/otbBasicFiltersExamplesTests.cxx
@@ -26,6 +26,7 @@
 void RegisterTests()
 {
   REGISTER_TEST(LeeImageFilterTest);
+  REGISTER_TEST(FrostImageFilterTest);
   REGISTER_TEST(DEMToRainbowExampleTest);
   REGISTER_TEST(ScalingFilterExampleTest);
   REGISTER_TEST(PrintableImageFilterExample1Test);
@@ -38,6 +39,10 @@ void RegisterTests()
 #define main LeeImageFilterTest
 #include "LeeImageFilter.cxx"
 
+#undef main
+#define main FrostImageFilterTest
+#include "FrostImageFilter.cxx"
+
 #undef main
 #define main DEMToRainbowExampleTest
 #include "DEMToRainbowExample.cxx"
diff --git a/Examples/FeatureExtraction/AlignmentsExample.cxx b/Examples/FeatureExtraction/AlignmentsExample.cxx
index fe908f463a491b80dba7aca07c91930b67dca677..bc441d79b1c38809169f3f3696b2cf320467c044 100644
--- a/Examples/FeatureExtraction/AlignmentsExample.cxx
+++ b/Examples/FeatureExtraction/AlignmentsExample.cxx
@@ -223,20 +223,13 @@ int main( int argc, char *argv[] )
   // \center
   // \includegraphics[width=0.35\textwidth]{QB_Suburb.eps}
   // \includegraphics[width=0.35\textwidth]{QB_SuburbAlign.eps}
-  // \itkcaption[Lee Filter Application]{Result of applying the
+  // \itkcaption[Alignment Detection Application]{Result of applying the
   // \doxygen{otb}{ImageToPathListAlignFilter} to a VHR image of a suburb.}
   // \label{fig:Align}
   // \end{figure}
   //
   //  Software Guide : EndLatex
 
-  //  \relatedClasses
-  //  \begin{itemize}
-  //  \item \doxygen{otb}{FrostImageFilter}
-  //  \end{itemize}
-  //
-  //  Software Guide : EndLatex
-
 
 
   return EXIT_SUCCESS;
diff --git a/Examples/FeatureExtraction/CMakeLists.txt b/Examples/FeatureExtraction/CMakeLists.txt
index a6fa86f323565e4ca500e0fba904b7d17b7adf20..2d9150a9cad302f3279e284476fe147cccb2b8b2 100644
--- a/Examples/FeatureExtraction/CMakeLists.txt
+++ b/Examples/FeatureExtraction/CMakeLists.txt
@@ -1,6 +1,10 @@
 PROJECT(FeatureExtractionExamples)
 INCLUDE_REGULAR_EXPRESSION("^.*$")
 
+ADD_EXECUTABLE(InnerProductPCAExample InnerProductPCAExample.cxx )
+TARGET_LINK_LIBRARIES(InnerProductPCAExample OTBCommon OTBIO OTBBasicFilters)
+
+
 ADD_EXECUTABLE(AlignmentsExample AlignmentsExample.cxx )
 TARGET_LINK_LIBRARIES(AlignmentsExample OTBFeatureExtraction OTBCommon OTBIO
 ITKIO)
@@ -96,6 +100,9 @@ TARGET_LINK_LIBRARIES(TextureExample OTBIO OTBCommon OTBFeatureExtraction ITKCom
 ADD_EXECUTABLE(EdgeDensityExample EdgeDensityExample.cxx)
 TARGET_LINK_LIBRARIES(EdgeDensityExample OTBIO OTBCommon OTBBasicFilters ITKCommon ITKBasicFilters)
 
+ADD_EXECUTABLE(SIFTDensityExample SIFTDensityExample.cxx)
+TARGET_LINK_LIBRARIES(SIFTDensityExample OTBIO OTBCommon OTBBasicFilters ITKCommon ITKBasicFilters)
+
 
 IF( NOT OTB_DISABLE_CXX_TESTING AND BUILD_TESTING )
 
diff --git a/Examples/FeatureExtraction/EdgeDensityExample.cxx b/Examples/FeatureExtraction/EdgeDensityExample.cxx
index 90e29f6fc8fbd279c7ce0b7ad9685f14e4a14f1a..c7f3a91f7688e565a83fcb9a267524b6b919d6fd 100644
--- a/Examples/FeatureExtraction/EdgeDensityExample.cxx
+++ b/Examples/FeatureExtraction/EdgeDensityExample.cxx
@@ -15,65 +15,222 @@ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 PURPOSE.  See the above copyright notices for more information.
 
 =========================================================================*/
-#include "otbImage.h"
-#include "otbEdgeDensityImageFilter.h"
-#include "otbBinaryImageDensityFunction.h"
-#include "itkCannyEdgeDetectionImageFilter.h"
 
+#include "otbImage.h"
 #include "otbImageFileReader.h"
 #include "otbImageFileWriter.h"
+#include "itkRescaleIntensityImageFilter.h"
+
+//  Software Guide : BeginCommandLineArgs
+//    INPUTS: {suburb2.jpeg}
+//    OUTPUTS: {EdgeDensityOutput.tif}, {PrettyEdgeDensityOutput.png}
+//    7 50 10 1.0 0.01
+//  Software Guide : EndCommandLineArgs
+
+// Software Guide : BeginLatex
+//
+// This example illustrates the use of the
+// \doxygen{otb}{EdgeDensityImageFilter}. 
+// This filter computes a local density of edges on an image and can
+// be useful to detect man made objects or urban areas, for
+// instance. The filter has been implemented in a generic way, so that
+// the way the edges are detected and the way their density is
+// computed can be chosen by the user.
+//
+// The first step required to use this filter is to include its header file.
+//
+// Software Guide : EndLatex
+
+// Software Guide : BeginCodeSnippet
+#include "otbEdgeDensityImageFilter.h"
+// Software Guide : EndCodeSnippet
+// Software Guide : BeginLatex
+//
+// We will also include the header files for the edge detector (a
+// Canny filter) and the density estimation (a simple count on a
+// binary image).
+//
+// The first step required to use this filter is to include its header file.
+//
+// Software Guide : EndLatex
+
+// Software Guide : BeginCodeSnippet
 
-
+#include "itkCannyEdgeDetectionImageFilter.h"
+#include "otbBinaryImageDensityFunction.h"
+// Software Guide : EndCodeSnippet
 
 int main(int argc, char* argv[] )
 {
 
   const  char *     infname       = argv[1];
   const  char *     outfname      = argv[2];
+  const  char *     prettyfilename      = argv[3];
+  
+  const unsigned int radius       = atoi(argv[4]);
+
 
-  /** Variables for the canny detector*/
-  const double    upperThreshold   = atof(argv[3]);
-  const double    lowerThreshold   = atof(argv[4]);
-  const double    variance         = atof(argv[5]);
-  const double    maximumError     = atof(argv[6]);
   /*--*/
     
-  const   unsigned int                                      Dimension = 2;
-  typedef float                                             PixelType;
+  const   unsigned int        Dimension = 2;
+  typedef float               PixelType;
+
+  /** Variables for the canny detector*/
+  const PixelType    upperThreshold   = static_cast<PixelType>(atof(argv[5]));
+  const PixelType    lowerThreshold   = static_cast<PixelType>(atof(argv[6]));
+  const double    variance         = atof(argv[7]);
+  const double    maximumError     = atof(argv[8]);
+
+
+  // Software Guide : BeginLatex
+  //
+  // As usual, we start by defining the types for the images, the reader
+  // and the writer.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  typedef otb::Image< PixelType, Dimension >    ImageType;
+  typedef otb::ImageFileReader<ImageType>       ReaderType;
+  typedef otb::ImageFileWriter<ImageType>       WriterType;
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // We define now the type for the function which will be used by the
+  // edge density filter to estimate this density. Here we choose a
+  // function which counts the number of non null pixels per area. The
+  // fucntion takes as template the type of the image to be processed.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  typedef otb::BinaryImageDensityFunction<ImageType> CountFunctionType;
+    // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // These {\em non null pixels} will be the result of an edge
+  // detector. We use here the classical Canny edge detector, which is
+  // templated over the input and output image types.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  typedef itk::CannyEdgeDetectionImageFilter<ImageType, ImageType>
+                                                     CannyDetectorType;
+
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // Finally, we can define the type for the edge density filter which
+  // takes as template the input and output image types, the edge
+  // detector type, and the count fucntion type..
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  typedef otb::EdgeDensityImageFilter<ImageType, ImageType, CannyDetectorType,
+                     CountFunctionType> EdgeDensityFilterType;
+
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // We can now instantiate the different processing objects of the
+  // pipeline using the \code{New()} method.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  ReaderType::Pointer               reader = ReaderType::New();
+  EdgeDensityFilterType::Pointer    filter = EdgeDensityFilterType::New();
+  CannyDetectorType::Pointer        cannyFilter = CannyDetectorType::New();
+  WriterType::Pointer               writer = WriterType::New();
+
+
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // The edge detection filter needs to be instantiated because we
+  // need to set its parameters. This is what we do here for the Canny
+  // filter.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  cannyFilter->SetUpperThreshold(upperThreshold); 
+  cannyFilter->SetLowerThreshold(lowerThreshold);
+  cannyFilter->SetVariance(variance); 
+  cannyFilter->SetMaximumError(maximumError); 
+
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // After that, we can pass the edge detector to the filter which
+  // will be used it internally.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  filter->SetDetector(cannyFilter);
+  filter->SetNeighborhoodRadius( radius );
+
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // Finally, we set the file names for the input and the output
+  // images and we plug the pipeline. The \code{Update()} method of
+  // the writer will trigger the processing.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  
+  reader->SetFileName(infname);
+  writer->SetFileName(outfname);
+  
+  
+  filter->SetInput(reader->GetOutput());
+  writer->SetInput(filter->GetOutput());
+  writer->Update();
+  // Software Guide : EndCodeSnippet
 
-  typedef otb::Image< PixelType, Dimension >                                      ImageType;
-  typedef ImageType::IndexType                                                    IndexType;
-  typedef otb::ImageFileReader<ImageType>                                         ReaderType;
-  typedef otb::ImageFileWriter<ImageType>                                         WriterType;
+  //  Software Guide : BeginLatex
+  // Figure~\ref{fig:EDGEDENSITY_FILTER} shows the result of applying
+  // the edge density filter to an image.
+  // \begin{figure}
+  // \center
+  // \includegraphics[width=0.25\textwidth]{suburb2.eps}
+  // \includegraphics[width=0.25\textwidth]{PrettyEdgeDensityOutput.eps}
+  // \itkcaption[Edge Density Filter]{Result of applying the
+  // \doxygen{otb}{EdgeDensityImageFilter} to an image. From left to right :
+  // original image, edge density.}
+  // \label{fig:EDGEDENSITY_FILTER}
+  // \end{figure}
+  //
+  //  Software Guide : EndLatex
 
-  typedef otb::BinaryImageDensityFunction<ImageType>                              CountFunctionType;
-  typedef itk::CannyEdgeDetectionImageFilter<ImageType , ImageType>               CannyDetectorType;
+  /************* Image for printing **************/
 
-  typedef otb::EdgeDensityImageFilter<ImageType,CannyDetectorType  ,CountFunctionType , ImageType> EdgeDensityFilterType;
+  typedef otb::Image< unsigned char, 2 > OutputImageType;
 
-  /**Instancitation of an object*/
-  EdgeDensityFilterType::Pointer    filter =      EdgeDensityFilterType::New();
-  ReaderType::Pointer               reader =      ReaderType::New();
-  CannyDetectorType::Pointer        CannyFilter = CannyDetectorType::New();
+  typedef itk::RescaleIntensityImageFilter< ImageType, OutputImageType >
+    RescalerType;
 
-  /** Set The input*/
-  reader->SetFileName(infname);
-  filter->SetInput(reader->GetOutput());
-  
-  /** Update the Canny Filter Information*/
-  CannyFilter->SetUpperThreshold(static_cast<ImageType::PixelType>(upperThreshold)); /** 30*/
-  CannyFilter->SetLowerThreshold(static_cast<ImageType::PixelType>(lowerThreshold));/** 10*/
-  CannyFilter->SetVariance(variance); //1.
-  CannyFilter->SetMaximumError(maximumError); ///0.01f
+  RescalerType::Pointer rescaler = RescalerType::New();
 
-  filter->SetDetector(CannyFilter);
+  rescaler->SetOutputMinimum(0);
+  rescaler->SetOutputMaximum(255);
 
-  /** Write the output*/
-  WriterType::Pointer          writer = WriterType::New();
-  writer->SetFileName(outfname);
-  writer->SetInput(filter->GetOutput());
-  writer->Update();
+  rescaler->SetInput( filter->GetOutput() );
 
+  typedef otb::ImageFileWriter< OutputImageType > OutputWriterType;
+  OutputWriterType::Pointer outwriter = OutputWriterType::New();
+
+  outwriter->SetFileName( prettyfilename );
+  outwriter->SetInput( rescaler->GetOutput() );
+  outwriter->Update();
+  
+  
+  
   return EXIT_SUCCESS;
 }
 
diff --git a/Examples/FeatureExtraction/InnerProductPCAExample.cxx b/Examples/FeatureExtraction/InnerProductPCAExample.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..008e4bf643740a82f326cbec412d7901c62492aa
--- /dev/null
+++ b/Examples/FeatureExtraction/InnerProductPCAExample.cxx
@@ -0,0 +1,184 @@
+/*=========================================================================
+
+  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 "otbImage.h"
+#include "otbVectorImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+#include "otbMultiToMonoChannelExtractROI.h"
+#include "itkRescaleIntensityImageFilter.h"
+
+//  Software Guide : BeginCommandLineArgs
+//    INPUTS: {ROI_QB_MUL_1.png}
+//    OUTPUTS: {InnerProductPCAOutput.tif}, {PrettyInnerProductPCAOutput1.png}, {PrettyInnerProductPCAOutput2.png}, {PrettyInnerProductPCAOutput3.png}
+//    3
+//  Software Guide : EndCommandLineArgs
+
+// Software Guide : BeginLatex
+//
+// This example illustrates the use of the
+// \doxygen{otb}{InnerProductPCAImageFilter}. 
+// This filter computes a Principal Component Analysis using an
+// efficient method based on the inner product in order to compute the
+// covariance matrix.
+//
+// The first step required to use this filter is to include its header file.
+//
+// Software Guide : EndLatex
+
+// Software Guide : BeginCodeSnippet
+#include "otbInnerProductPCAImageFilter.h"
+// Software Guide : EndCodeSnippet
+
+int main( int argc, char* argv[] )
+{
+  typedef double PixelType;
+  const unsigned int Dimension = 2;
+  const char * inputFileName = argv[1];
+  const char * outputFilename = argv[2];
+  const unsigned int numberOfPrincipalComponentsRequired(atoi(argv[6]));
+
+
+  // Software Guide : BeginLatex
+  //
+  // We start by defining the types for the images and the reader and
+  // the writer. We choose to work with a \doxygen{otb}{VectorImage},
+  // since we will produce a multi-channel image (the principal
+  // components) from a multi-channel input image.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  typedef otb::VectorImage<PixelType,Dimension> ImageType;
+  typedef otb::ImageFileReader< ImageType >     ReaderType;
+  typedef otb::ImageFileWriter< ImageType >     WriterType;
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // We instantiate now the image reader and we set the image file name.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  ReaderType::Pointer     reader     = ReaderType::New();
+  reader->SetFileName(inputFileName);
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // We define the type for the filter. It is templated over the input
+  // and the output image types. We the instantiate the filter.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  typedef otb::InnerProductPCAImageFilter<ImageType,ImageType> PCAFilterType;
+  PCAFilterType::Pointer     pcafilter     = PCAFilterType::New();
+
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // The only parameter needed for the PCA is the number of principal
+  // components required as output. We can choose to get less PCs than
+  // the number of input bands.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+
+  pcafilter->SetNumberOfPrincipalComponentsRequired(
+                                    numberOfPrincipalComponentsRequired);
+
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // We now instantiate the writer and set the file name for the
+  // output image.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+
+  WriterType::Pointer     writer     = WriterType::New();
+  writer->SetFileName(outputFilename);
+
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // We finally plug the pipeline and trigger the PCA computation with
+  // the method \code{Update()} of the writer.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+
+  pcafilter->SetInput(reader->GetOutput());
+  writer->SetInput(pcafilter->GetOutput());
+  
+  writer->Update();
+
+  // Software Guide : EndCodeSnippet
+
+  //  Software Guide : BeginLatex
+  // Figure~\ref{fig:INNERPRODUCTPCA_FILTER} shows the result of applying
+  // the PCA to a 3 band RGB image.
+  // \begin{figure}
+  // \center
+  // \includegraphics[width=0.25\textwidth]{ROI_QB_MUL_1.eps}
+  // \includegraphics[width=0.25\textwidth]{PrettyInnerProductPCAOutput1.eps}
+  // \includegraphics[width=0.25\textwidth]{PrettyInnerProductPCAOutput2.eps}
+  // \includegraphics[width=0.25\textwidth]{PrettyInnerProductPCAOutput3.eps}
+  // \itkcaption[Inner Product PCA Filter]{Result of applying the
+  // \doxygen{otb}{InnerProductPCAImageFilter} to an image. From left
+  // to right and top to bottom: 
+  // original image, first PC, second PC, third PC.}
+  // \label{fig:INNERPRODUCTPCA_FILTER}
+  // \end{figure}
+  //
+  //  Software Guide : EndLatex
+
+
+
+  typedef otb::Image<PixelType,Dimension> MonoImageType;
+  
+  typedef otb::MultiToMonoChannelExtractROI< PixelType, PixelType >
+                                                    ExtractROIFilterType;
+
+  typedef otb::Image<unsigned char, 2> OutputImageType;
+  typedef otb::ImageFileWriter< OutputImageType >  WriterType2;
+  typedef itk::RescaleIntensityImageFilter< MonoImageType, OutputImageType> RescalerType;
+  
+  for(unsigned int cpt=0 ; cpt < numberOfPrincipalComponentsRequired ; cpt++)
+  {
+    ExtractROIFilterType::Pointer extractROIFilter = ExtractROIFilterType::New();
+    RescalerType::Pointer rescaler = RescalerType::New();
+    
+    WriterType2::Pointer     writer2     = WriterType2::New();
+    extractROIFilter->SetInput(pcafilter->GetOutput());
+    extractROIFilter->SetChannel(cpt+1);
+
+    rescaler->SetInput( extractROIFilter->GetOutput());
+    rescaler->SetOutputMinimum(0);
+    rescaler->SetOutputMaximum(255);
+    
+    writer2->SetInput(rescaler->GetOutput());
+    writer2->SetFileName(argv[cpt+3]);
+    writer2->Update();
+  }
+
+
+  return EXIT_SUCCESS;
+}
diff --git a/Examples/FeatureExtraction/LineSegmentDetectorExample.cxx b/Examples/FeatureExtraction/LineSegmentDetectorExample.cxx
index 79efad01aa885f35baff2ac617d0b5537b597f59..67fb52c99ff000f8433165a0c5e14d2d63a64548 100644
--- a/Examples/FeatureExtraction/LineSegmentDetectorExample.cxx
+++ b/Examples/FeatureExtraction/LineSegmentDetectorExample.cxx
@@ -17,55 +17,161 @@
 =========================================================================*/
 
 #include "otbImage.h"
-#include "otbLineSegmentDetector.h"
 #include "otbDrawLineSpatialObjectListFilter.h"
 #include "otbLineSpatialObjectList.h"
 
 #include "otbImageFileReader.h"
 #include "otbImageFileWriter.h"
 
+//  Software Guide : BeginCommandLineArgs
+//    INPUTS: {QB_Suburb.png}
+//    OUTPUTS: {QB_SuburbLSD.png}
+//  Software Guide : EndCommandLineArgs
+
+// Software Guide : BeginLatex
+//
+// This example illustrates the use of the
+// \doxygen{otb}{LineSegmentDetector}\cite{LSD}, also known as {\em Lucy in the
+// Sky with Diamonds}.
+//
+// The first step required to use this filter is to include its header file.
+//
+// Software Guide : EndLatex
+
+// Software Guide : BeginCodeSnippet
+#include "otbLineSegmentDetector.h"
+// Software Guide : EndCodeSnippet
+
 int main( int argc, char * argv[] )
 {
   const   char * infname  = argv[1];
   const   char * outfname  = argv[2];
   
   typedef unsigned char   InputPixelType;
+  typedef double   PrecisionType;
   const   unsigned int    Dimension = 2;
   
-  /** Typedefs */
-  typedef otb::Image< InputPixelType,  Dimension >    InputImageType;
-  typedef otb::ImageFileReader<InputImageType>        ReaderType;
+
+// Software Guide : BeginLatex
+//
+// As usual, we start by defining the types for the input image and
+// the image file reader.
+//
+// Software Guide : EndLatex
+
+// Software Guide : BeginCodeSnippet  
+  typedef otb::Image< InputPixelType,  Dimension >    ImageType;
+  typedef otb::ImageFileReader<ImageType>        ReaderType;
+// Software Guide : EndCodeSnippet    
+// Software Guide : BeginLatex
+//
+// We instantiate the reader and set the file name for the input image.
+//
+// Software Guide : EndLatex
+
+// Software Guide : BeginCodeSnippet  
+
   ReaderType::Pointer reader = ReaderType::New();
   reader->SetFileName(infname);
-  
-  
-  
+// Software Guide : EndCodeSnippet    
+// Software Guide : BeginLatex
+//
+// We define now the type for the segment detector filter. It is
+// templated over the input image type and the precision with which
+// the coordinates of the detected segments will be given. It is
+// recommended to set this precision to a real type. The output of the
+// filter will be a list of \doxygen{itk}{LineSpatialObject}s.
+//
+// Software Guide : EndLatex
+
+// Software Guide : BeginCodeSnippet  
+    
+  typedef otb::LineSegmentDetector<ImageType,
+                                                 PrecisionType> LsdFilterType;
+
 
-  typedef otb::LineSegmentDetector<InputImageType,
-                                                 InputPixelType> lsdFilterType;
+  
+  LsdFilterType::Pointer  lsdFilter = LsdFilterType::New();
 
+  // Software Guide : EndCodeSnippet    
+// Software Guide : BeginLatex
+//
+// In order to be able to display the results, we will draw the
+// detected segments on top of the input image. For this matter, we
+// will use the \doxygen{otb}{DrawLineSpatialObjectListFilter} which
+// is templated over the input and output image types.
+//
+// Software Guide : EndLatex
 
-  /** Instanciation of smart pointer*/
-  lsdFilterType::Pointer  lsdFilter = lsdFilterType::New();
+// Software Guide : BeginCodeSnippet  
 
 
-  typedef otb::DrawLineSpatialObjectListFilter< InputImageType,
-                                       InputImageType > DrawLineListType;
+  typedef otb::DrawLineSpatialObjectListFilter< ImageType,
+                                       ImageType > DrawLineListType;
   DrawLineListType::Pointer drawLineFilter =   DrawLineListType::New();
 
 
-  reader->GenerateOutputInformation();
-  lsdFilter->SetInput(reader->GetOutput());
-  drawLineFilter->SetInput(reader->GetOutput());
-  drawLineFilter->SetInputLineSpatialObjectList(lsdFilter->GetOutput());
+  // Software Guide : EndCodeSnippet    
+// Software Guide : BeginLatex
+//
+// We can now define the type for the writer, instantiate it and set
+// the file name for the output image.
+//
+// Software Guide : EndLatex
 
- 
-  typedef otb::ImageFileWriter<InputImageType>        WriterType;
+// Software Guide : BeginCodeSnippet  
+
+  typedef otb::ImageFileWriter<ImageType>        WriterType;
   WriterType::Pointer writer = WriterType::New();
   writer->SetFileName(outfname);
+
+  // Software Guide : EndCodeSnippet    
+// Software Guide : BeginLatex
+//
+// We plug the pipeline.
+//
+// Software Guide : EndLatex
+
+// Software Guide : BeginCodeSnippet  
+
+  lsdFilter->SetInput(reader->GetOutput());
   writer->SetInput(drawLineFilter->GetOutput());
+
+  drawLineFilter->SetInput(reader->GetOutput());
+  drawLineFilter->SetInputLineSpatialObjectList(lsdFilter->GetOutput());
+
+  // Software Guide : EndCodeSnippet    
+// Software Guide : BeginLatex
+//
+// Before calling the \code{Update()} method of the writer in order to
+// trigger the pipeline execution, we call the
+// \doxygen{GenerateOutputInformation()} of the reader, so the LSD
+// filter gets the information about image size and spacing.
+//
+// Software Guide : EndLatex
+
+// Software Guide : BeginCodeSnippet  
+
+
+  reader->GenerateOutputInformation();
   writer->Update();
 
+  // Software Guide : EndCodeSnippet      
+
+  //  Software Guide : BeginLatex
+  // Figure~\ref{fig:LSD} shows the result of applying the line segment
+  // detection to a small patch extracted from a VHR image.
+  // \begin{figure}
+  // \center
+  // \includegraphics[width=0.35\textwidth]{QB_Suburb.eps}
+  // \includegraphics[width=0.35\textwidth]{QB_SuburbLSD.eps}
+  // \itkcaption[LSD Application]{Result of applying the
+  // \doxygen{otb}{LineSegmentDetector} to a VHR image of a suburb.}
+  // \label{fig:LSD}
+  // \end{figure}
+  //
+  //  Software Guide : EndLatex
+
   
   return EXIT_SUCCESS;
 }
diff --git a/Examples/FeatureExtraction/SIFTDensityExample.cxx b/Examples/FeatureExtraction/SIFTDensityExample.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..345f209965448aa577362424aa9e2aa186b05038
--- /dev/null
+++ b/Examples/FeatureExtraction/SIFTDensityExample.cxx
@@ -0,0 +1,235 @@
+/*=========================================================================
+
+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 "itkPointSet.h"
+#include "itkVariableLengthVector.h"
+#include "otbImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+#include "itkRescaleIntensityImageFilter.h"
+
+//  Software Guide : BeginCommandLineArgs
+//    INPUTS: {suburb2.jpeg}
+//    OUTPUTS: {SIFTDensityOutput.tif}, {PrettySIFTDensityOutput.png}
+//    3 3 7
+//  Software Guide : EndCommandLineArgs
+
+// Software Guide : BeginLatex
+//
+// This example illustrates the use of the
+// \doxygen{otb}{KeyPointDensityImageFilter}. 
+// This filter computes a local density of keypoints (SIFT or SURF,
+// for instance) on an image and can
+// be useful to detect man made objects or urban areas, for
+// instance. The filter has been implemented in a generic way, so that
+// the way the keypoints are detected can be chosen by the user.
+//
+// The first step required to use this filter is to include its header file.
+//
+// Software Guide : EndLatex
+
+// Software Guide : BeginCodeSnippet
+#include "otbKeyPointDensityImageFilter.h"
+// Software Guide : EndCodeSnippet
+#include "otbImageToSIFTKeyPointSetFilter.h"
+
+int main(int argc, char* argv[] )
+{
+  const char * infname = argv[1];
+  const char * outfname = argv[2];
+  const char * prettyfname = argv[3];
+  const unsigned int scales = atoi(argv[4]);
+  const unsigned int octaves = atoi(argv[5]);
+  const unsigned int radius = atoi(argv[6]);
+
+  const   unsigned int         Dimension = 2;
+  typedef float                PixelType;
+
+  // Software Guide : BeginLatex
+  //
+  // As usual, we start by defining the types for the images, the reader
+  // and the writer.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  typedef otb::Image< PixelType, Dimension >         ImageType;
+  typedef otb::ImageFileReader<ImageType>            ReaderType;
+  typedef otb::ImageFileWriter<ImageType>            WriterType;
+
+    // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // We define now the type for the keypoint detection. The keypoints
+  // will be stored in vector form (they may contain many descriptors)
+  // into a point set. The filter for detecting the SIFT is templated
+  // over the input image type and the output pointset type.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  
+  typedef itk::VariableLengthVector<PixelType>       RealVectorType;
+  typedef itk::PointSet<RealVectorType,Dimension>   PointSetType;
+  typedef otb::ImageToSIFTKeyPointSetFilter<ImageType, PointSetType>
+                                                        DetectorType;
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // We can now define the filter which will compute the SIFT
+  // density. It will be templated over the input and output image
+  // types and the SIFT detector.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+  typedef otb::KeyPointDensityImageFilter<ImageType, ImageType, DetectorType>
+                                                                   FilterType;
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // We can instantiate the reader and the writer as wella s the
+  // filter and the detector. The detector needs to be instantiated in
+  // order to set its parameters.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+
+
+  ReaderType::Pointer      reader = ReaderType::New();
+  WriterType::Pointer writer = WriterType::New();
+  FilterType::Pointer    filter =     FilterType::New();
+  DetectorType::Pointer  detector = DetectorType::New();
+
+
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // We set the file names for the input and the output images.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+
+  reader->SetFileName(infname);
+  writer->SetFileName(outfname);
+
+
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // We set the parameters for the SIFT detector (the number of
+  // octaves and the number of scales per octave).
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+
+
+  detector ->SetOctavesNumber(octaves);
+  detector->SetScalesNumber(scales);
+
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // And we pass the detector to the filter and we set the radius for
+  // the density estimation.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+
+
+  filter->SetDetector(detector);
+  filter->SetNeighborhoodRadius(radius);
+
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // We plug the pipeline.
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+
+
+  filter->SetInput(reader->GetOutput());
+  writer->SetInput(filter->GetOutput());
+
+  // Software Guide : EndCodeSnippet
+  // Software Guide : BeginLatex
+  //
+  // We trigger the execution by calling th \code{Update()} method on
+  // the writer, but before that we run the
+  // \code{GenerateOutputInformation()} on the reader so the filter
+  // gets the information about the image size (needed for the SIFT
+  // computation). 
+  //
+  // Software Guide : EndLatex
+
+  // Software Guide : BeginCodeSnippet
+
+  reader->GenerateOutputInformation();
+  writer->Update();
+  // Software Guide : EndCodeSnippet
+
+  //  Software Guide : BeginLatex
+  // Figure~\ref{fig:SIFTDENSITY_FILTER} shows the result of applying
+  // the key point density filter to an image using the SIFT
+  // detector. The main difference with respect to figure
+  // \ref{fig:EDGEDENSITY_FILTER} is that for SIFTS, individual trees
+  // contribute to the density.
+  // \begin{figure}
+  // \center
+  // \includegraphics[width=0.25\textwidth]{suburb2.eps}
+  // \includegraphics[width=0.25\textwidth]{PrettySIFTDensityOutput.eps}
+  // \itkcaption[SIFT Density Filter]{Result of applying the
+  // \doxygen{otb}{KeypointDensityImageFilter} to an image. From left
+  // to right : 
+  // original image, SIF density.}
+  // \label{fig:SIFTDENSITY_FILTER}
+  // \end{figure}
+  //
+  //  Software Guide : EndLatex
+
+  /************* Image for printing **************/
+
+  typedef otb::Image< unsigned char, 2 > OutputImageType;
+
+  typedef itk::RescaleIntensityImageFilter< ImageType, OutputImageType >
+    RescalerType;
+
+  RescalerType::Pointer rescaler = RescalerType::New();
+
+  rescaler->SetOutputMinimum(0);
+  rescaler->SetOutputMaximum(255);
+
+  rescaler->SetInput( filter->GetOutput() );
+
+  typedef otb::ImageFileWriter< OutputImageType > OutputWriterType;
+  OutputWriterType::Pointer outwriter = OutputWriterType::New();
+
+  outwriter->SetFileName( prettyfname );
+  outwriter->SetInput( rescaler->GetOutput() );
+  outwriter->Update();
+  
+ 
+  return EXIT_SUCCESS;
+}
+
diff --git a/Testing/Code/BasicFilters/CMakeLists.txt b/Testing/Code/BasicFilters/CMakeLists.txt
index 9e50512b6c60317544a82ca231cf101c18be2356..cbb5b9a68c1251d81a2df8fd12dad9fe8f58a4e6 100644
--- a/Testing/Code/BasicFilters/CMakeLists.txt
+++ b/Testing/Code/BasicFilters/CMakeLists.txt
@@ -1025,7 +1025,7 @@ ADD_TEST(bfTvEdgeDensityImageFilter ${BASICFILTERS_TESTS9}
         otbEdgeDensityImageFilter
 	    ${INPUTDATA}/scene.png
 	    ${TEMP}/bfTvEdgeDensityImageFilterOutputImage.tif
-	    2 # radius
+	    1 # radius
 	    15 3  1. 0.01  #Canny Parameters
 	    
 )
diff --git a/Testing/Code/BasicFilters/otbEdgeDensityImageFilter.cxx b/Testing/Code/BasicFilters/otbEdgeDensityImageFilter.cxx
index b7a000d2089e1f65d6ec3b41c17b2c5367a4484b..05189e887ed4c7e5dfc6c1f827e43467d0e85270 100644
--- a/Testing/Code/BasicFilters/otbEdgeDensityImageFilter.cxx
+++ b/Testing/Code/BasicFilters/otbEdgeDensityImageFilter.cxx
@@ -50,7 +50,7 @@ int otbEdgeDensityImageFilter(int argc, char* argv[] )
   typedef otb::BinaryImageDensityFunction<ImageType>                              CountFunctionType;
   typedef itk::CannyEdgeDetectionImageFilter<ImageType , ImageType>               CannyDetectorType;
 
-  typedef otb::EdgeDensityImageFilter<ImageType,CannyDetectorType  ,CountFunctionType , ImageType> EdgeDensityFilterType;
+  typedef otb::EdgeDensityImageFilter<ImageType, ImageType, CannyDetectorType  ,CountFunctionType> EdgeDensityFilterType;
 
   /**Instancitation of an object*/
   EdgeDensityFilterType::Pointer    filter =      EdgeDensityFilterType::New();
@@ -68,6 +68,7 @@ int otbEdgeDensityImageFilter(int argc, char* argv[] )
   CannyFilter->SetMaximumError(maximumError); ///0.01f
 
   filter->SetDetector(CannyFilter);
+  filter->SetNeighborhoodRadius(radius);
 
   /** Write the output*/
   WriterType::Pointer          writer = WriterType::New();
diff --git a/Testing/Code/BasicFilters/otbEdgeDensityImageFilterNew.cxx b/Testing/Code/BasicFilters/otbEdgeDensityImageFilterNew.cxx
index c324cd6d861145a39f102d10c9e977b85aa7c729..35318804e87487281bd129a7a81d8d04d0c78bac 100644
--- a/Testing/Code/BasicFilters/otbEdgeDensityImageFilterNew.cxx
+++ b/Testing/Code/BasicFilters/otbEdgeDensityImageFilterNew.cxx
@@ -35,7 +35,7 @@ int otbEdgeDensityImageFilterNew(int, char* [] )
   typedef otb::BinaryImageDensityFunction<ImageType>                              CountFunctionType;
   typedef itk::CannyEdgeDetectionImageFilter<ImageType , ImageType>               CannyDetectorType;
 
-  typedef otb::EdgeDensityImageFilter<ImageType,CannyDetectorType  ,CountFunctionType , ImageType> EdgeDensityFilterType;
+  typedef otb::EdgeDensityImageFilter<ImageType , ImageType , CannyDetectorType  ,CountFunctionType> EdgeDensityFilterType;
 
   /**Instancitation of an object*/
   EdgeDensityFilterType::Pointer    filter =      EdgeDensityFilterType::New();
diff --git a/Testing/Code/BasicFilters/otbKeyPointDensityImageFilterNew.cxx b/Testing/Code/BasicFilters/otbKeyPointDensityImageFilterNew.cxx
index d2742f30c3f77ae849853e330aefab65200ca71b..a70c2bdb725c7e4bf29cc966e4a76b8ce5f7c99c 100644
--- a/Testing/Code/BasicFilters/otbKeyPointDensityImageFilterNew.cxx
+++ b/Testing/Code/BasicFilters/otbKeyPointDensityImageFilterNew.cxx
@@ -36,7 +36,7 @@ int otbKeyPointDensityImageFilterNew(int, char* [] )
   typedef itk::PointSet<RealVectorType,Dimension>                               PointSetType;
   typedef otb::ImageToSIFTKeyPointSetFilter<ImageType,PointSetType>             DetectorType;
 
-  typedef otb::KeyPointDensityImageFilter< ImageType,DetectorType, ImageType>   FilterType;
+  typedef otb::KeyPointDensityImageFilter< ImageType, ImageType, DetectorType>   FilterType;
 
   /**Instancitation of an object*/
   FilterType::Pointer    filter =     FilterType::New();
diff --git a/Testing/Code/BasicFilters/otbKeyPointDensityImageFilterTest.cxx b/Testing/Code/BasicFilters/otbKeyPointDensityImageFilterTest.cxx
index 9c036dc1469ddab2f091f274780ef4bcd782e09b..bbd2ba825cd995eed882502b29e9647ad0837a0c 100644
--- a/Testing/Code/BasicFilters/otbKeyPointDensityImageFilterTest.cxx
+++ b/Testing/Code/BasicFilters/otbKeyPointDensityImageFilterTest.cxx
@@ -46,7 +46,7 @@ int otbKeyPointDensityImageFilterTest(int argc, char* argv[] )
   typedef itk::PointSet<RealVectorType,Dimension>                               PointSetType;
   typedef otb::ImageToSIFTKeyPointSetFilter<ImageType,PointSetType>             DetectorType;
 
-  typedef otb::KeyPointDensityImageFilter< ImageType,DetectorType, ImageType>   FilterType;
+  typedef otb::KeyPointDensityImageFilter< ImageType, ImageType, DetectorType>   FilterType;
 
 
   /**Instancitation of an object*/
diff --git a/Testing/Code/Common/CMakeLists.txt b/Testing/Code/Common/CMakeLists.txt
index f03e2cc3fdc058a22aaf26241da783b93e07c3df..fb56092e5fba57230b9d5a542bf1abe399d1c133 100644
--- a/Testing/Code/Common/CMakeLists.txt
+++ b/Testing/Code/Common/CMakeLists.txt
@@ -1,4 +1,3 @@
-
 IF( NOT OTB_DISABLE_CXX_TESTING AND BUILD_TESTING )
 
 SET(BASELINE ${OTB_DATA_ROOT}/Baseline/OTB/Images)
@@ -6,6 +5,11 @@ SET(BASELINE_FILES ${OTB_DATA_ROOT}/Baseline/OTB/Files)
 SET(INPUTDATA ${OTB_DATA_ROOT}/Input)
 SET(TEMP ${OTBTesting_BINARY_DIR}/Temporary)
 
+#Remote sensing images (large images )
+IF(OTB_DATA_USE_LARGEINPUT)
+  SET(LARGEDATA ${OTB_DATA_LARGEINPUT_ROOT} )
+ENDIF(OTB_DATA_USE_LARGEINPUT)
+
 #Tolerance sur diff pixel image
 SET(TOL 0.0)
 
@@ -123,6 +127,25 @@ ADD_TEST(coTvExtractROI_RGB2 ${COMMON_TESTS2}
         ${TEMP}/coExtractROI_RGB_poupees_302_1_134_330.jpg
         303 2 134 330 )
 
+# -------            otb::VectorDataExtractROI   ------------------------------
+ADD_TEST(coTuVectorDataExtractROINew ${COMMON_TESTS2}
+	otbVectorDataExtractROINew)
+
+IF(OTB_DATA_USE_LARGEINPUT)
+
+ADD_TEST(coTvVectorDataExtractROI ${COMMON_TESTS2}
+#  --compare-image ${TOL}   
+#              ${BASELINE}/coVectorDataExtractROIOutput.shp 
+#              ${TEMP}/coVectorDataExtractROIOutput.shp 
+	 otbVectorDataExtractROI	
+        ${LARGEDATA}/TOULOUSE/QuickBird/GIS_FILES/000000128955_01_ORDER_SHAPE.shp
+	${TEMP}/coVectorDataExtractROIOutput.shp
+	10 10 50 50 
+ )
+
+ENDIF(OTB_DATA_USE_LARGEINPUT)
+
+
 
 # -------            otb::MultiChannelExtractROI   ------------------------------
 
@@ -592,6 +615,8 @@ otbExtractROI_RGB.cxx
 otbMultiChannelExtractROI.cxx
 otbMultiChannelExtractROINew.cxx
 otbTestMultiExtractMultiUpdate.cxx
+otbVectorDataExtractROINew.cxx
+otbVectorDataExtractROI.cxx	
 )
 SET(BasicCommon_SRCS3
 otbMultiToMonoChannelExtractROI.cxx
diff --git a/Testing/Code/Common/otbCommonTests2.cxx b/Testing/Code/Common/otbCommonTests2.cxx
index 88349a63bbabd391b0deb7b0a9d58013b2b295ab..808364dae8bf8a0c2613ab4aef4ef14ce0ae5b01 100644
--- a/Testing/Code/Common/otbCommonTests2.cxx
+++ b/Testing/Code/Common/otbCommonTests2.cxx
@@ -34,4 +34,6 @@ void RegisterTests()
   REGISTER_TEST(otbMultiChannelExtractROI );
   REGISTER_TEST(otbMultiChannelExtractROINew );
   REGISTER_TEST(otbTestMultiExtractMultiUpdate);
+  REGISTER_TEST(otbVectorDataExtractROINew);
+  REGISTER_TEST(otbVectorDataExtractROI);
 }
diff --git a/Testing/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.cxx b/Testing/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.cxx
index ab7bc454c8a3a76b3c488984abb11f80afb24b9a..2f811e0cba52b89a179b576dedb5c0fe9720ce9f 100644
--- a/Testing/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.cxx
+++ b/Testing/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilter.cxx
@@ -25,15 +25,20 @@
 
 namespace Functor
 {
-template <class TIter, class TOutput>
+template <class TInput, class TOutput>
 class UnaryFunctorNeighborhoodWithOffseImageFilterTest
 {
 public:
   UnaryFunctorNeighborhoodWithOffseImageFilterTest() {};
   ~UnaryFunctorNeighborhoodWithOffseImageFilterTest() {};
 
-  typedef TIter IterType;
-  typedef typename IterType::OffsetType OffsetType;
+  typedef TInput                                      InputScalarType;
+  typedef TOutput                                     OutputScalarType;
+  typedef itk::VariableLengthVector<InputScalarType>  InputVectorType;
+  typedef itk::VariableLengthVector<OutputScalarType> OutputVectorType;
+  typedef itk::Offset<>                               OffsetType;
+  typedef itk::Neighborhood<InputScalarType, 2>       NeighborhoodType;
+  typedef itk::Neighborhood<InputVectorType, 2>       NeighborhoodVectorType;
 
   void SetOffset(OffsetType off)
   {
@@ -44,9 +49,13 @@ public:
     return m_Offset;
   };
 
-  inline TOutput operator() (const TIter & it)
+  inline OutputScalarType operator() (const NeighborhoodType & neigh)
   {
-    return(static_cast<TOutput>(it.GetCenterPixel()));
+    return(static_cast<OutputScalarType>(neigh.GetCenterValue()));
+  }
+  inline OutputVectorType operator() (const NeighborhoodVectorType & neigh)
+  {
+    return(static_cast<OutputVectorType>(neigh.GetCenterValue()));
 
   }
 private:
@@ -68,7 +77,7 @@ int otbUnaryFunctorNeighborhoodWithOffsetImageFilter(int argc, char * argv[])
   typedef otb::ImageFileReader<ImageType>            ReaderType;
   typedef otb::ImageFileWriter<ImageType>            WriterType;
   typedef itk::ConstNeighborhoodIterator<ImageType>  IterType;;
-  typedef Functor::UnaryFunctorNeighborhoodWithOffseImageFilterTest<IterType, PixelType>  FunctorType;
+  typedef Functor::UnaryFunctorNeighborhoodWithOffseImageFilterTest<InputPixelType, InputPixelType>  FunctorType;
   typedef otb::UnaryFunctorNeighborhoodWithOffsetImageFilter<ImageType, ImageType, FunctorType>     UnaryFunctorNeighborhoodImageFilterType;
 
   // Instantiating object
diff --git a/Testing/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilterNew.cxx b/Testing/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilterNew.cxx
index 930bd2cf0ef7203a91ea880bced15aedea1b2288..52e3a337dcaefc602d58dc65a5a9af623a283214 100644
--- a/Testing/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilterNew.cxx
+++ b/Testing/Code/Common/otbUnaryFunctorNeighborhoodWithOffsetImageFilterNew.cxx
@@ -20,18 +20,24 @@
 #include "otbUnaryFunctorNeighborhoodWithOffsetImageFilter.h"
 #include "otbVectorImage.h"
 #include "itkConstNeighborhoodIterator.h"
+#include "itkOffset.h"
 
 namespace Functor
 {
-template <class TIter, class TOutput>
+template <class TInput, class TOutput>
 class UnaryFunctorNeighborhoodWithOffsetImageFilterFunctorNewTest
 {
 public:
   UnaryFunctorNeighborhoodWithOffsetImageFilterFunctorNewTest() {};
   ~UnaryFunctorNeighborhoodWithOffsetImageFilterFunctorNewTest() {};
 
-  typedef TIter IterType;
-  typedef typename IterType::OffsetType OffsetType;
+  typedef TInput                                      InputScalarType;
+  typedef TOutput                                     OutputScalarType;
+  typedef itk::VariableLengthVector<InputScalarType>  InputVectorType;
+  typedef itk::VariableLengthVector<OutputScalarType> OutputVectorType;
+  typedef itk::Offset<>                               OffsetType;
+  typedef itk::Neighborhood<InputScalarType, 2>       NeighborhoodType;
+  typedef itk::Neighborhood<InputVectorType, 2>       NeighborhoodVectorType;
 
   void SetOffset(OffsetType off)
   {
@@ -42,11 +48,16 @@ public:
     return m_Offset;
   };
 
-  inline TOutput operator() (const TIter & it)
+  inline OutputScalarType operator() (const NeighborhoodType & neigh)
   {
-    return(static_cast<TOutput>(it.GetCenterPixel()));
+    return(static_cast<OutputScalarType>(neigh.GetCenterValue()));
+  }
+  inline OutputVectorType operator() (const NeighborhoodVectorType & neigh)
+  {
+    return(static_cast<OutputVectorType>(neigh.GetCenterValue()));
 
   }
+
 private:
   OffsetType m_Offset;
 };
@@ -60,7 +71,7 @@ int otbUnaryFunctorNeighborhoodWithOffsetImageFilterNew(int argc, char * argv[])
   typedef otb::VectorImage<InputPixelType,Dimension> ImageType;
   typedef ImageType::PixelType PixelType;
   typedef itk::ConstNeighborhoodIterator<ImageType>   IterType;;
-  typedef Functor::UnaryFunctorNeighborhoodWithOffsetImageFilterFunctorNewTest<IterType, PixelType>  FunctorType;
+  typedef Functor::UnaryFunctorNeighborhoodWithOffsetImageFilterFunctorNewTest<InputPixelType, InputPixelType>  FunctorType;
   typedef otb::UnaryFunctorNeighborhoodWithOffsetImageFilter<ImageType, ImageType, FunctorType> UnaryFunctorNeighborhoodWithOffsetImageFilterType;
 
   // Instantiating object
diff --git a/Testing/Code/Common/otbVectorDataExtractROI.cxx b/Testing/Code/Common/otbVectorDataExtractROI.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b646a04193745d163ac1f3192ca49a8539ed760e
--- /dev/null
+++ b/Testing/Code/Common/otbVectorDataExtractROI.cxx
@@ -0,0 +1,54 @@
+/*=========================================================================
+
+  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 "otbVectorData.h"
+#include "otbVectorDataExtractROI.h"
+
+#include "otbVectorDataFileReader.h"
+
+int otbVectorDataExtractROI( int argc, char * argv[] )
+{
+  const char * infname = argv[1];
+
+  const unsigned int  startX = atoi(argv[3]);
+  const unsigned int  startY = atoi(argv[4]);
+  const unsigned int  sizeX = atoi(argv[5]);
+  const unsigned int  sizeY = atoi(argv[6]);
+  
+  typedef otb::VectorData<>                                  VectorDataType;
+  typedef otb::VectorDataExtractROI< VectorDataType >        FilterType;
+  typedef otb::VectorDataFileReader<VectorDataType>          VectorDataFileReaderType;
+
+  /** */
+  FilterType::Pointer filter = FilterType::New();
+  VectorDataFileReaderType::Pointer reader = VectorDataFileReaderType::New();
+  
+  /** */
+  reader->SetFileName(infname);
+  
+  filter->SetInput(reader->GetOutput());
+  filter->SetSizeX(sizeX);
+  filter->SetSizeY(sizeY);
+  filter->SetStartX(startX);
+  filter->SetStartY(startY);
+
+  filter->Update();
+  
+  return EXIT_SUCCESS;
+}
+
+
diff --git a/Testing/Code/Common/otbVectorDataExtractROINew.cxx b/Testing/Code/Common/otbVectorDataExtractROINew.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..ee2f998747267d85da991483305a4f4678e78f51
--- /dev/null
+++ b/Testing/Code/Common/otbVectorDataExtractROINew.cxx
@@ -0,0 +1,33 @@
+/*=========================================================================
+
+  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 "otbVectorData.h"
+#include "otbVectorDataExtractROI.h"
+
+int otbVectorDataExtractROINew( int argc, char * argv[] )
+{
+  typedef otb::VectorData<> VectorDataType;
+  typedef otb::VectorDataExtractROI< VectorDataType >  FilterType;
+  
+
+  FilterType::Pointer filter = FilterType::New();
+
+
+  return EXIT_SUCCESS;
+}
+
+
diff --git a/Testing/Code/FeatureExtraction/CMakeLists.txt b/Testing/Code/FeatureExtraction/CMakeLists.txt
index 7ccc32d480df5249975346f1e43bc31f50a74dff..6cc7bbb0987f2a2d9d05887183236f5796820c64 100644
--- a/Testing/Code/FeatureExtraction/CMakeLists.txt
+++ b/Testing/Code/FeatureExtraction/CMakeLists.txt
@@ -395,7 +395,7 @@ ADD_TEST(feTvLocalHoughDrawBis ${FEATUREEXTRACTION_TESTS5}
 ADD_TEST(feTuFillGapsFilterNew ${FEATUREEXTRACTION_TESTS5}
         otbFillGapsFilterNew)
 
-ADD_TEST(feTuFillGapsFilter ${FEATUREEXTRACTION_TESTS5}
+ADD_TEST(feTvFillGapsFilter ${FEATUREEXTRACTION_TESTS5}
 	  --compare-image ${EPSILON}
 	       ${BASELINE}/feTvFillGapfilter.png
                ${TEMP}/feTvFillGapfilter.png
@@ -1089,6 +1089,17 @@ ADD_TEST(feTvLineSpatialObjectListToRightAnglePointSetFilterBySteps ${FEATUREEXT
     ${TEMP}/feTvLineSpatialObjectListToRightAnglePointSetFilterByStepsOutputAscii.txt
 )
 
+#ADD_TEST(feTvLineSpatialObjectListToRightAnglePointSetFilterOutputImage ${FEATUREEXTRACTION_TESTS11} 
+# --compare-ascii ${EPS}
+#	${BASELINE_FILES}/feTvLineSpatialObjectListToRightAnglePointSetFilterByStepsOutputAscii.txt
+#	${TEMP}/feTvLineSpatialObjectListToRightAnglePointSetFilterByStepsOutputAscii.txt
+#    otbLineSpatialObjectListToRightAnglePointSetFilterOutputImage
+#    ${INPUTDATA}/prison_toulouse.tif
+#    ${TEMP}/feTvLineSpatialObjectListToRightAnglePointSetFilterByStepsOutputAscii.txt
+#    ${TEMP}/feTvLineSpatialObjectListToRightAnglePointSetFilterOutputImage.tif
+#)
+
+
 
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ otbFeatureExtractionTests12 ~~~~~~~~~~~~~~~~~~~~~
@@ -1634,17 +1645,31 @@ ADD_TEST(feTuLineDirectionImageFilterNew ${FEATUREEXTRACTION_TESTS14}
     otbLineDirectionImageFilterNew)
 
 ADD_TEST(feTvLineDirectionImageFilterTest ${FEATUREEXTRACTION_TESTS14} 
- #--compare-image ${EPS}
- #               ${BASELINE}/feTvLineDirectionImageFilterTest.tif
- #               ${TEMP}/feTvLineDirectionImageFilterTest.tif
+--compare-n-images 0 #${EPS}
+		   5
+                   ${BASELINE}/feTvLineDirectionLengthImageFilterTest.tif
+		   ${TEMP}/feTvLineDirectionLengthImageFilterTest.tif
+		   ${BASELINE}/feTvLineDirectionWidthImageFilterTest.tif
+		   ${TEMP}/feTvLineDirectionWidthImageFilterTest.tif
+		   ${BASELINE}/feTvLineDirectionWMeanImageFilterTest.tif
+		   ${TEMP}/feTvLineDirectionWMeanImageFilterTest.tif
+		   ${BASELINE}/feTvLineDirectionRatioImageFilterTest.tif
+		   ${TEMP}/feTvLineDirectionRatioImageFilterTest.tif
+		   ${BASELINE}/feTvLineDirectionSDImageFilterTest.tif
+		   ${TEMP}/feTvLineDirectionSDImageFilterTest.tif
    otbLineDirectionImageFilterTest
-    ${INPUTDATA}/poupees_sub.png
-    ${TEMP}/feTvLineDirectionImageFilterTest.tif
-    8 # spectral threshold
-    3 # spatial threshold
+    ${INPUTDATA}/poupees_sub_c1.png
+    ${TEMP}/feTvLineDirectionLengthImageFilterTest.tif
+    ${TEMP}/feTvLineDirectionWidthImageFilterTest.tif
+    ${TEMP}/feTvLineDirectionWMeanImageFilterTest.tif
+    ${TEMP}/feTvLineDirectionRatioImageFilterTest.tif
+    ${TEMP}/feTvLineDirectionSDImageFilterTest.tif
+    50 # spectral threshold
+    20 # spatial threshold
     15 # direction
-    5 # max min/max number takes in care for ratio
+    4 # max min/max number takes in care for ratio
     0.6 # alpha value
+    # PSI is disable for texture selection test
 )
 
 
@@ -1782,6 +1807,7 @@ otbLineSegmentDetector.cxx
 otbLineSpatialObjectListToRightAnglePointSetFilterNew.cxx
 otbLineSpatialObjectListToRightAnglePointSetFilter.cxx
 otbLineSpatialObjectListToRightAnglePointSetFilterByStepsOutputAscii.cxx
+#otbLineSpatialObjectListToRightAnglePointSetFilterOutputImage.cxx
 )
 SET(BasicFeatureExtraction_SRCS12
 otbTextureFunctorBase.cxx
diff --git a/Testing/Code/FeatureExtraction/otbFeatureExtractionTests11.cxx b/Testing/Code/FeatureExtraction/otbFeatureExtractionTests11.cxx
index 0c6d94a4c37a113a006186f2f3a74d2da4448534..8cfeaf90f11696cbb2bad192c626426f52fcc754 100644
--- a/Testing/Code/FeatureExtraction/otbFeatureExtractionTests11.cxx
+++ b/Testing/Code/FeatureExtraction/otbFeatureExtractionTests11.cxx
@@ -38,4 +38,5 @@ REGISTER_TEST(otbLineSegmentDetector);
 REGISTER_TEST(otbLineSpatialObjectListToRightAnglePointSetFilterNew);
 REGISTER_TEST(otbLineSpatialObjectListToRightAnglePointSetFilter);
 REGISTER_TEST(otbLineSpatialObjectListToRightAnglePointSetFilterByStepsOutputAscii);
+//REGISTER_TEST(otbLineSpatialObjectListToRightAnglePointSetFilterOutputImage);
 }
diff --git a/Testing/Code/FeatureExtraction/otbLineDirectionImageFilterNew.cxx b/Testing/Code/FeatureExtraction/otbLineDirectionImageFilterNew.cxx
index 7c9010a0fed549446799fbdbf8f48ab99a8cbeb3..d92e877f50e027f395b1ec01dcacbedf709b5ee6 100644
--- a/Testing/Code/FeatureExtraction/otbLineDirectionImageFilterNew.cxx
+++ b/Testing/Code/FeatureExtraction/otbLineDirectionImageFilterNew.cxx
@@ -26,31 +26,12 @@ int otbLineDirectionImageFilterNew(int argc, char * argv[])
 {
   const unsigned int Dimension =2;
   typedef double PixelType;
-  typedef otb::VectorImage<PixelType,Dimension>                ImageType;
+  typedef otb::Image<PixelType,Dimension>                ImageType;
   typedef otb::VectorImage<PixelType,Dimension> VectorImageType;
   //typedef otb::StreamingImageFileWriter<VectorImageType> WriterType;
-  typedef otb::LineDirectionImageFilter<ImageType, VectorImageType> FilterType;
+  typedef otb::LineDirectionImageFilter<ImageType, ImageType> FilterType;
   FilterType::Pointer filter = FilterType::New();
 
- /*
-  typedef otb::ImageFileReader<ImageType> ReaderType;
 
-  typedef otb::LineDirectionFunctor<VectorImageTypeImageType, ImageType> FunctorType;
-  typedef otb::UnaryFunctorNeighborhoodImageFilter<ImageType, FunctorType> FilterType;
-
-
-  // Instantiating object
-  FilterType::Pointer filter = FilterType::New();
-  
-  ReaderType::Pointer reader = ReaderType::New();
-  WriterType::Pointer writer = WriterType::New();
-
-  reader->SetFileName(argv[1]);
-  writer->SetFileName(argv[2]);
-  filter->SetInput( reader->GetOutput() );
-  writer->SetInput( panTex->GetOutput() );
-
-  writer->Update();
-  */
   return EXIT_SUCCESS;
 }
diff --git a/Testing/Code/FeatureExtraction/otbLineDirectionImageFilterTest.cxx b/Testing/Code/FeatureExtraction/otbLineDirectionImageFilterTest.cxx
index 7f68a7b2beceeadbe09dc75dfd212e6de36e9709..cb90dbc0c6e0ec2687295088d6a40b5a53b9f21c 100644
--- a/Testing/Code/FeatureExtraction/otbLineDirectionImageFilterTest.cxx
+++ b/Testing/Code/FeatureExtraction/otbLineDirectionImageFilterTest.cxx
@@ -22,7 +22,7 @@
 #include "otbImage.h"
 #include "otbVectorImage.h"
 #include "otbImageFileReader.h"
-//#include "otbStreamingImageFileWriter.h"
+#include "otbStreamingImageFileWriter.h"
 #include "otbImageFileWriter.h"
 
 
@@ -32,47 +32,67 @@ int otbLineDirectionImageFilterTest(int argc, char * argv[])
   const unsigned int Dimension =2;
 
   std::string inName            = argv[1];
-  std::string outName           = argv[2];
-  PixelType spectThresh         = atof(argv[3]);
-  unsigned int spatialThresh    = atoi(argv[4]);
-  unsigned int dirNb            = atoi(argv[5]);
-  unsigned int maxConsideration = atoi(argv[6]);
-  double alpha                  = atof(argv[7]);  
-
-
-  typedef otb::VectorImage<PixelType,Dimension>                           ImageType;
+  std::string outNameLength     = argv[2];
+  std::string outNameWidth      = argv[3];
+  std::string outNameWMean      = argv[4];
+  std::string outNameRatio      = argv[5];
+  std::string outNameSD         = argv[6];
+  PixelType spectThresh         = atof(argv[7]);
+  unsigned int spatialThresh    = atoi(argv[8]);
+  unsigned int dirNb            = atoi(argv[9]);
+  unsigned int maxConsideration = atoi(argv[10]);
+  double alpha                  = atof(argv[11]);  
+
+
+  typedef otb::Image<PixelType,Dimension>                           ImageType;
   typedef ImageType::PixelType                                InputPixelType;
   typedef otb::VectorImage<PixelType,Dimension>                     VectorImageType;
   typedef otb::ImageFileReader<ImageType>                           ReaderType;
   //typedef otb::StreamingImageFileWriter<VectorImageType>            WriterType;
-  typedef otb::ImageFileWriter<VectorImageType>            WriterType;
-  typedef otb::LineDirectionImageFilter<ImageType, VectorImageType> FilterType;
+  typedef otb::ImageFileWriter<ImageType>            WriterType;
+  typedef otb::LineDirectionImageFilter<ImageType, ImageType> FilterType;
 
   FilterType::Pointer filter = FilterType::New(); 
   ReaderType::Pointer reader = ReaderType::New();
-  WriterType::Pointer writer = WriterType::New();
-
-
-
+  WriterType::Pointer writerLength = WriterType::New();
+  WriterType::Pointer writerWidth = WriterType::New();
+  WriterType::Pointer writerWMean = WriterType::New();
+  WriterType::Pointer writerRatio = WriterType::New();
+  WriterType::Pointer writerSD = WriterType::New();
 
   reader->SetFileName(inName);
   reader->GenerateOutputInformation();
-  writer->SetFileName(outName);
-
-  InputPixelType spect;
-  // TO MODIFY
-  //spect.SetSize(reader->GetOutput()->GetNumberOfComponentsPerPixel());
-  //spect.Fill(spectThresh);
+   
   filter->SetSpectralThreshold(spectThresh);
   filter->SetSpatialThreshold(spatialThresh);
   filter->SetNumberOfDirections(dirNb);
   filter->SetRatioMaxConsiderationNumber(maxConsideration);
   filter->SetAlpha(alpha);
+  // disable PSI texture
+  filter->SetTextureStatus(3, false);
   filter->SetInput( reader->GetOutput() );
-  writer->SetInput( filter->GetOutput() );
 
-  writer->Update();
- 
+
+  writerLength->SetFileName(outNameLength);
+  writerLength->SetInput( filter->GetLengthOutput() );
+  writerLength->Update();
+
+  writerWidth->SetFileName(outNameWidth);
+  writerWidth->SetInput( filter->GetWidthOutput() );
+  writerWidth->Update();
+
+  writerWMean->SetFileName(outNameWMean);
+  writerWMean->SetInput( filter->GetWMeanOutput() );
+  writerWMean->Update();
+
+  writerRatio->SetFileName(outNameRatio);
+  writerRatio->SetInput( filter->GetRatioOutput() );
+  writerRatio->Update();
+
+  writerSD->SetFileName(outNameSD);
+  writerSD->SetInput( filter->GetSDOutput() );
+  writerSD->Update();
+  
 
   return EXIT_SUCCESS;
 }
diff --git a/Testing/Code/FeatureExtraction/otbPanTexTextureImageFilter.cxx b/Testing/Code/FeatureExtraction/otbPanTexTextureImageFilter.cxx
index 50b5c649c0151e8b171ad2a8a0cd57e29935d592..5978bc794f781adc030790755faebdc1b3cf6db2 100644
--- a/Testing/Code/FeatureExtraction/otbPanTexTextureImageFilter.cxx
+++ b/Testing/Code/FeatureExtraction/otbPanTexTextureImageFilter.cxx
@@ -24,7 +24,7 @@
 
 
 int otbPanTexTextureImageFilter(int argc, char * argv[])
-{
+{/*
   const unsigned int Dimension =2;
   typedef double PixelType;
   typedef otb::VectorImage<PixelType,Dimension> ImageType;
@@ -43,6 +43,6 @@ int otbPanTexTextureImageFilter(int argc, char * argv[])
   writer->SetInput( panTex->GetOutput() );
 
   writer->Update();
-
+ */
   return EXIT_SUCCESS;
 }
diff --git a/Testing/Code/FeatureExtraction/otbPanTexTextureImageFilterNew.cxx b/Testing/Code/FeatureExtraction/otbPanTexTextureImageFilterNew.cxx
index 9511c798fa2713bba8607dc183afa278d5e13574..3fa6a922ad18a9708cdede6d12248cb2feada52a 100644
--- a/Testing/Code/FeatureExtraction/otbPanTexTextureImageFilterNew.cxx
+++ b/Testing/Code/FeatureExtraction/otbPanTexTextureImageFilterNew.cxx
@@ -22,6 +22,7 @@
 
 int otbPanTexTextureImageFilterNew(int argc, char * argv[])
 {
+  /*
   const unsigned int Dimension =2;
   typedef double PixelType;
   typedef otb::VectorImage<PixelType,Dimension> ImageType;
@@ -29,7 +30,7 @@ int otbPanTexTextureImageFilterNew(int argc, char * argv[])
 
   // Instantiating object
   PanTexType::Pointer object = PanTexType::New();
-
+  */
 
   return EXIT_SUCCESS;
 }
diff --git a/Testing/Code/FeatureExtraction/otbPanTexTextureImageFunctionFilter.cxx b/Testing/Code/FeatureExtraction/otbPanTexTextureImageFunctionFilter.cxx
index 78e1bf5f6e35992c678dda0040ed837b6861890d..ec7a933e46f465a409d123947d72ea044809cf15 100644
--- a/Testing/Code/FeatureExtraction/otbPanTexTextureImageFunctionFilter.cxx
+++ b/Testing/Code/FeatureExtraction/otbPanTexTextureImageFunctionFilter.cxx
@@ -25,6 +25,7 @@
 
 int otbPanTexTextureImageFunctionFilter(int argc, char * argv[])
 {
+  /*
   const unsigned int Dimension =2;
   typedef double PixelType;
   typedef otb::Image<PixelType,Dimension> ImageType;
@@ -43,6 +44,6 @@ int otbPanTexTextureImageFunctionFilter(int argc, char * argv[])
   writer->SetInput( panTex->GetOutput() );
 
   writer->Update();
-
+  */
   return EXIT_SUCCESS;
 }
diff --git a/Testing/Code/FeatureExtraction/otbPanTexTextureImageFunctionFilterNew.cxx b/Testing/Code/FeatureExtraction/otbPanTexTextureImageFunctionFilterNew.cxx
index 0a99f2d9039acd0e693d2a7b36b0e339939afee8..a2543a86a60c2050da4941a4ec4c6aa588257ea4 100644
--- a/Testing/Code/FeatureExtraction/otbPanTexTextureImageFunctionFilterNew.cxx
+++ b/Testing/Code/FeatureExtraction/otbPanTexTextureImageFunctionFilterNew.cxx
@@ -22,6 +22,7 @@
 
 int otbPanTexTextureImageFunctionFilterNew(int argc, char * argv[])
 {
+  /*
   const unsigned int Dimension =2;
   typedef double PixelType;
   typedef otb::Image<PixelType,Dimension> ImageType;
@@ -29,7 +30,7 @@ int otbPanTexTextureImageFunctionFilterNew(int argc, char * argv[])
 
   // Instantiating object
   PanTexType::Pointer object = PanTexType::New();
-
+  */
 
   return EXIT_SUCCESS;
 }
diff --git a/Testing/Code/FeatureExtraction/otbTextureFunctor.cxx b/Testing/Code/FeatureExtraction/otbTextureFunctor.cxx
index 9fe9feeed77776db9995a906675f17a37ca26980..9aaf2ad01b566802bf49685714e6ab7e84676b4f 100644
--- a/Testing/Code/FeatureExtraction/otbTextureFunctor.cxx
+++ b/Testing/Code/FeatureExtraction/otbTextureFunctor.cxx
@@ -70,15 +70,15 @@ int otbTextureFunctor(int argc, char * argv[])
   typedef otb::VectorImage<InputPixelType,Dimension> ImageType;
   typedef ImageType::PixelType                       PixelType;
   typedef itk::ConstNeighborhoodIterator<ImageType>  IteratorType;
-
+  
   if(strArgv == "ENJ")
     {
-      typedef otb::Functor::EnergyTextureFunctor<IteratorType, PixelType> FunctorType;
+      typedef otb::Functor::EnergyTextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureFunctor<ImageType, ImageType, FunctorType>(argc,argv) );
     }
   else if ( strArgv == "ENT" )
     {
-      typedef otb::Functor::EntropyTextureFunctor<IteratorType, PixelType> FunctorType;
+      typedef otb::Functor::EntropyTextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureFunctor<ImageType, ImageType, FunctorType>(argc,argv) );
     }
   else if ( strArgv == "IMD" )
@@ -86,6 +86,7 @@ int otbTextureFunctor(int argc, char * argv[])
       typedef otb::Functor::InverseDifferenceMomentTextureFunctor<IteratorType, PixelType> FunctorType;
       return( generic_TextureFunctor<ImageType, ImageType, FunctorType>(argc,argv) );
     }
+  /*
  else if ( strArgv == "ASM" )
     {
       typedef otb::Functor::AngularSecondMomentumTextureFunctor<IteratorType, PixelType> FunctorType;
@@ -156,10 +157,12 @@ int otbTextureFunctor(int argc, char * argv[])
       typedef otb::Functor::MeanTextureFunctor<IteratorType, PixelType> FunctorType;
       return( generic_TextureFunctor<ImageType, ImageType, FunctorType>(argc,argv) );
     }
+  */
   else
     {
       return EXIT_FAILURE;
     }
   
+
   return EXIT_SUCCESS;
 }
diff --git a/Testing/Code/FeatureExtraction/otbTextureFunctorBase.cxx b/Testing/Code/FeatureExtraction/otbTextureFunctorBase.cxx
index 44e139e3da2c63fa9c27b5131e5c7f8edfd7c642..ff6e103ff4851d6f73dc268741e2e9b0c0fa6422 100644
--- a/Testing/Code/FeatureExtraction/otbTextureFunctorBase.cxx
+++ b/Testing/Code/FeatureExtraction/otbTextureFunctorBase.cxx
@@ -25,28 +25,22 @@
 #include "otbTextureFunctorBase.h"
 
 
-template <class TIterInput, class TOutput>
+template <class TScalarInput, class TScalarOutput>
 class ITK_EXPORT TextureFunctorTest : 
-public otb::Functor::TextureFunctorBase<TIterInput, TOutput>
+public otb::Functor::TextureFunctorBase<TScalarInput, TScalarOutput>
 {
 public:
   TextureFunctorTest()
     {};
   ~TextureFunctorTest(){};
 
-  typedef TIterInput                           IterType;
-  typedef TOutput                              OutputType;
-  typedef typename IterType::OffsetType        OffsetType;
-  typedef typename IterType::RadiusType        RadiusType;
-  typedef typename IterType::InternalPixelType InternalPixelType;
-  typedef typename IterType::ImageType         ImageType;
-  typedef itk::Neighborhood<InternalPixelType,::itk::GetImageDimension<ImageType>::ImageDimension>    NeighborhoodType;
+  typedef itk::Neighborhood<TScalarInput, 2>    NeighborhoodType;
   
   virtual double ComputeOverSingleChannel(const NeighborhoodType &neigh, const NeighborhoodType &neighOff)
   {
     double out = neigh.GetCenterValue(); 
     
-	return out;
+    return out;
   }
 };
 
@@ -65,7 +59,7 @@ int otbTextureFunctorBase(int argc, char * argv[])
   typedef otb::ImageFileWriter<ImageType> WriterType;
 
   typedef itk::ConstNeighborhoodIterator<ImageType>   IterType;;
-  typedef TextureFunctorTest<IterType, PixelType>  FunctorType;
+  typedef TextureFunctorTest<InputPixelType, InputPixelType>  FunctorType;
   typedef otb::UnaryFunctorNeighborhoodWithOffsetImageFilter<ImageType, ImageType, FunctorType> UnaryFunctorNeighborhoodImageFilterType;
 
   // Instantiating object
diff --git a/Testing/Code/FeatureExtraction/otbTextureImageFunction.cxx b/Testing/Code/FeatureExtraction/otbTextureImageFunction.cxx
index 3525d085f80b4d74f571256b7ad0bd812d823d0c..1e19ea35870d83dc72e61ce817ab84670165c780 100644
--- a/Testing/Code/FeatureExtraction/otbTextureImageFunction.cxx
+++ b/Testing/Code/FeatureExtraction/otbTextureImageFunction.cxx
@@ -41,8 +41,6 @@ int generic_TextureImageFunction(int argc, char * argv[])
 
   typedef otb::TextureImageFunction<TInputImage, TFunctor> FunctionType;
   typedef otb::FunctionWithNeighborhoodToImageFilter<TInputImage, TOutputImage, FunctionType> FilterType;
-  //typename FunctionType::Pointer energyFunction = FunctionType::New();
-  //typename FunctionType energyFunction;
   typename FilterType::Pointer filter = FilterType::New();
 
   // Instantiating object
@@ -52,7 +50,6 @@ int generic_TextureImageFunction(int argc, char * argv[])
   writer->SetFileName(outputFileName);
 
   filter->SetInput(reader->GetOutput());
-  //filter->SetFunction(energyFunction);
   SizeType radius;
   radius[0] = atoi(argv[3]);
   radius[1] = atoi(argv[4]);
@@ -83,92 +80,95 @@ int otbTextureImageFunction(int argc, char * argv[])
   typedef itk::VariableLengthVector<double> VectorType;
   typedef itk::ConstNeighborhoodIterator<ImageType> IteratorType;
  
-
+  
   if(strArgv == "ENJ")
     {
-      typedef otb::Functor::EnergyTextureFunctor<IteratorType, VectorType> FunctorType;
+      typedef otb::Functor::EnergyTextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
+  
   else if ( strArgv == "ENT" )
     {
-      typedef otb::Functor::EntropyTextureFunctor<IteratorType, VectorType> FunctorType;
+      typedef otb::Functor::EntropyTextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
   else if ( strArgv == "IMD" )
     {
-      typedef otb::Functor::InverseDifferenceMomentTextureFunctor<IteratorType, VectorType> FunctorType;
+      typedef otb::Functor::InverseDifferenceMomentTextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
+  /*
  else if ( strArgv == "ASM" )
     {
-      typedef otb::Functor::AngularSecondMomentumTextureFunctor<IteratorType, VectorType> FunctorType;
+      typedef otb::Functor::AngularSecondMomentumTextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
  else if ( strArgv == "VAR" )
     {
-      typedef otb::Functor::VarianceTextureFunctor<IteratorType, VectorType> FunctorType;
+      typedef otb::Functor::VarianceTextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
  else if ( strArgv == "COR" )
     {
-      typedef otb::Functor::CorrelationTextureFunctor<IteratorType, VectorType> FunctorType;
+      typedef otb::Functor::CorrelationTextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
   else if ( strArgv == "CON" )
     {
-      typedef otb::Functor::ContrastTextureFunctor<IteratorType, VectorType> FunctorType;
+      typedef otb::Functor::ContrastTextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
   else if ( strArgv == "SAV" )
     {
-      typedef otb::Functor::SumAverageTextureFunctor<IteratorType, VectorType> FunctorType;
+      typedef otb::Functor::SumAverageTextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
   else if ( strArgv == "DEN" )
     {
-      typedef otb::Functor::DifferenceEntropyTextureFunctor<IteratorType, VectorType> FunctorType;
+      typedef otb::Functor::DifferenceEntropyTextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
   else if ( strArgv == "SEN" )
     {
-      typedef otb::Functor::SumEntropyTextureFunctor<IteratorType, VectorType> FunctorType;
+      typedef otb::Functor::SumEntropyTextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
   else if ( strArgv == "SVA" )
     {
-      typedef otb::Functor::SumVarianceTextureFunctor<IteratorType, VectorType> FunctorType;
+      typedef otb::Functor::SumVarianceTextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
   else if ( strArgv == "DVA" )
     {
-      typedef otb::Functor::DifferenceVarianceTextureFunctor<IteratorType, VectorType> FunctorType;
+      typedef otb::Functor::DifferenceVarianceTextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
   else if ( strArgv == "IC1" )
     {
-      typedef otb::Functor::InformationMeasureOfCorrelation1TextureFunctor<IteratorType, VectorType> FunctorType;
+      typedef otb::Functor::InformationMeasureOfCorrelation1TextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
   else if ( strArgv == "IC2" )
     {
-      typedef otb::Functor::InformationMeasureOfCorrelation2TextureFunctor<IteratorType, VectorType> FunctorType;
+      typedef otb::Functor::InformationMeasureOfCorrelation2TextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
   else if ( strArgv == "CSH" )
     {
-      typedef otb::Functor::ClusterShadeTextureFunctor<IteratorType, VectorType> FunctorType;
+      typedef otb::Functor::ClusterShadeTextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
   else if ( strArgv == "CPR" )
     {
-      typedef otb::Functor::ClusterProminenceTextureFunctor<IteratorType, VectorType> FunctorType;
+      typedef otb::Functor::ClusterProminenceTextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
   else if ( strArgv == "MEA" )
     {
-      typedef otb::Functor::MeanTextureFunctor<IteratorType, VectorType> FunctorType;
+      typedef otb::Functor::MeanTextureFunctor<InputPixelType, InputPixelType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
+  */
   else
     {
       return EXIT_FAILURE;
diff --git a/Testing/Code/Radiometry/CMakeLists.txt b/Testing/Code/Radiometry/CMakeLists.txt
index 51c975d14eb7e3f91e8123cf92620d2ce8fc53bc..3049644a02a85e1a184fb3b458a64497cc2db15c 100644
--- a/Testing/Code/Radiometry/CMakeLists.txt
+++ b/Testing/Code/Radiometry/CMakeLists.txt
@@ -870,6 +870,29 @@ ADD_TEST(raTvISU_MultiChannelRAndNIRBuiltUpIndexImageFilter ${RADIOMETRY_TESTS8}
 )
 
 
+# -------            NDWI MutliChannel WaterIndexFilter ------------------------------
+ADD_TEST(raTvNDWI_MultiChannelWaterIndexImageFilter ${RADIOMETRY_TESTS8}  
+ --compare-image ${EPSILON}   ${BASELINE}/raMultiChannelWaterIndex_NDWI_qb_RoadExtract.tif
+                   ${TEMP}/raMultiChannelWaterIndex_NDWI_qb_RoadExtract.tif
+        otbNDWIMultiChannelWaterIndexImageFilter
+        ${INPUTDATA}/qb_RoadExtract2sub200x200.tif
+        ${TEMP}/raMultiChannelWaterIndex_NDWI_qb_RoadExtract.tif
+        3 4 # mir nir channels
+)
+
+
+# -------            NDWI SingleChannel WaterIndexFilter ------------------------------
+ADD_TEST(raTvNDWI_WaterIndexImageFilter ${RADIOMETRY_TESTS8}  
+ --compare-image ${EPSILON}   ${BASELINE}/raWaterIndex_NDWI_verySmallFSATSW.tif
+                   ${TEMP}/raWaterIndex_NDWI_verySmallFSATSW.tif
+        otbNDWIWaterIndexImageFilter
+        ${INPUTDATA}/verySmallFSATSW_r.tif
+        ${INPUTDATA}/verySmallFSATSW_nir.tif
+        ${TEMP}/raWaterIndex_NDWI_verySmallFSATSW.tif
+)
+
+
+
 # A enrichir
 SET(Radiometry_SRCS1
 otbMultiChannelRAndNIRVegetationIndexImageFilterNew.cxx
@@ -940,6 +963,8 @@ otbNDBITM4AndTM5IndexImageFilter.cxx
 otbNDBIMultiChannelTM4AndTM5IndexImageFilter.cxx
 otbISURAndNIRIndexImageFilter.cxx
 otbISUMultiChannelRAndNIRIndexImageFilter.cxx
+otbNDWIMultiChannelWaterIndexImageFilter
+otbNDWIWaterIndexImageFilter
 )
 
 
diff --git a/Testing/Code/Radiometry/otbNDWIMultiChannelWaterIndexImageFilter.cxx b/Testing/Code/Radiometry/otbNDWIMultiChannelWaterIndexImageFilter.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..7eb2307a73d7befde0f3bdbbaaf000a125a7ea7d
--- /dev/null
+++ b/Testing/Code/Radiometry/otbNDWIMultiChannelWaterIndexImageFilter.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 "itkUnaryFunctorImageFilter.h"
+#include "otbImage.h"
+#include "otbVectorImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+#include "otbWaterIndicesFunctor.h"
+
+
+int otbNDWIMultiChannelWaterIndexImageFilter(int argc, char * argv[])
+{
+  const unsigned int Dimension = 2;
+  typedef otb::VectorImage<double ,Dimension>    InputImageType;
+  typedef otb::Image<double,Dimension>           OutputImageType;
+  typedef otb::ImageFileReader<InputImageType>   ReaderType;
+  typedef otb::ImageFileWriter<OutputImageType>  WriterType;
+  typedef otb::Functor::NDWI < InputImageType::InternalPixelType,
+                              InputImageType::InternalPixelType,
+                              OutputImageType::PixelType > FunctorType;
+
+  typedef itk::UnaryFunctorImageFilter<InputImageType,OutputImageType,FunctorType>
+                                                             UnaryFunctorImageFilterType;
+
+  // Instantiating object
+  UnaryFunctorImageFilterType::Pointer filter = UnaryFunctorImageFilterType::New();
+  ReaderType::Pointer reader = ReaderType::New();
+  WriterType::Pointer writer = WriterType::New();
+
+  const char * inputFilename  = argv[1];
+  const char * outputFilename = argv[2];
+
+  unsigned int nirChannel(::atoi(argv[3]));
+  unsigned int mirChannel(::atoi(argv[4]));
+
+  reader->SetFileName( inputFilename );
+  writer->SetFileName( outputFilename  );
+  filter->GetFunctor().SetNIRIndex(nirChannel);
+  filter->GetFunctor().SetMIRIndex(mirChannel);
+  filter->SetInput( reader->GetOutput() );
+
+  writer->SetInput( filter->GetOutput() );
+  writer->Update();
+
+  return EXIT_SUCCESS;
+}
+
+
+
+
diff --git a/Testing/Code/Radiometry/otbNDWIWaterIndexImageFilter.cxx b/Testing/Code/Radiometry/otbNDWIWaterIndexImageFilter.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..cd416f8d28845e2c3052ba12097b4c0b16e0e525
--- /dev/null
+++ b/Testing/Code/Radiometry/otbNDWIWaterIndexImageFilter.cxx
@@ -0,0 +1,71 @@
+/*=========================================================================
+
+  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 "itkBinaryFunctorImageFilter.h"
+#include "otbImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+#include "otbWaterIndicesFunctor.h"
+
+
+int otbNDWIWaterIndexImageFilter(int argc, char * argv[])
+{
+  const unsigned int                      Dimension = 2;
+  typedef double                          PixelType;
+  typedef otb::Image<PixelType,Dimension> InputNIRImageType;
+  typedef otb::Image<PixelType,Dimension> InputMIRImageType;
+  typedef otb::Image<double,Dimension>    OutputImageType;
+
+  typedef otb::ImageFileReader<InputNIRImageType>  NIRReaderType;
+  typedef otb::ImageFileReader<InputMIRImageType>  MIRReaderType;
+  typedef otb::ImageFileWriter<OutputImageType>    WriterType;
+
+  typedef otb::Functor::NDWI  < InputNIRImageType::PixelType,
+                                InputMIRImageType::PixelType,
+                                OutputImageType::PixelType > FunctorType;
+
+  typedef itk::BinaryFunctorImageFilter< InputNIRImageType,
+                                         InputMIRImageType,
+                                         OutputImageType,
+                                         FunctorType > BinaryFunctorImageFilterType;
+
+  // Instantiating object
+  BinaryFunctorImageFilterType::Pointer filter = BinaryFunctorImageFilterType::New();
+  NIRReaderType::Pointer readerNIR = NIRReaderType::New();
+  MIRReaderType::Pointer readerMIR = MIRReaderType::New();
+  WriterType::Pointer writer = WriterType::New();
+
+
+  const char * inputFilenameNIR  = argv[1];
+  const char * inputFilenameMIR  = argv[2];
+  const char * outputFilename = argv[3];
+
+  readerNIR->SetFileName( inputFilenameNIR );
+  readerMIR->SetFileName( inputFilenameMIR );
+  writer->SetFileName( outputFilename  );
+  filter->SetInput1( readerNIR->GetOutput() );
+  filter->SetInput2( readerMIR->GetOutput() );
+
+  writer->SetInput( filter->GetOutput() );
+  writer->Update();
+
+
+  return EXIT_SUCCESS;
+}
+
diff --git a/Testing/Code/Radiometry/otbRadiometryTests8.cxx b/Testing/Code/Radiometry/otbRadiometryTests8.cxx
index 84105e0bf4fdb4c3634c52f8837ec13aafc01977..bc83872c661f701d637df4ec85b94226772f7e11 100644
--- a/Testing/Code/Radiometry/otbRadiometryTests8.cxx
+++ b/Testing/Code/Radiometry/otbRadiometryTests8.cxx
@@ -31,6 +31,8 @@ void RegisterTests()
   REGISTER_TEST(otbNDBIMultiChannelTM4AndTM5IndexImageFilter);
   REGISTER_TEST(otbISURAndNIRIndexImageFilter);
   REGISTER_TEST(otbISUMultiChannelRAndNIRIndexImageFilter);
+  REGISTER_TEST(otbNDWIMultiChannelWaterIndexImageFilter);
+  REGISTER_TEST(otbNDWIWaterIndexImageFilter);
 }