diff --git a/Code/FeatureExtraction/otbLineCorrelationDetector.h b/Code/FeatureExtraction/otbLineCorrelationDetector.h
new file mode 100755
index 0000000000000000000000000000000000000000..7189964d49163a6149f562b58b81ed271b7497c7
--- /dev/null
+++ b/Code/FeatureExtraction/otbLineCorrelationDetector.h
@@ -0,0 +1,176 @@
+/*=========================================================================
+
+ Programme : OTB (ORFEO ToolBox)
+ Auteurs : CS - C.Ruffel
+ Language : C++
+ Date : 14 mars 2006
+ Role : Filter of detection of linear features
+ $Id$
+
+=========================================================================*/
+#ifndef __otbLineCorrelationDetector_h
+#define __otbLineCorrelationDetector_h
+
+#include "itkBSplineInterpolateImageFunction.h"
+#include "itkImageToImageFilter.h"
+#include "itkImage.h"
+#include "itkNumericTraits.h"
+
+#define MIN(_A,_B) ((_A) < (_B) ? (_A) : (_B))
+#define MAX(_A,_B) ((_A) > (_B) ? (_A) : (_B))
+#define ROTATION(_x,_y,_theta,_xout,_yout) \
+ (_xout) = (_x)*cos(_theta) - (_y)*sin(_theta); \
+ (_yout) = (_x)*sin(_theta) + (_y)*cos(_theta)
+
+
+namespace otb
+{
+
+/** \class LineCorrelationDetector
+ * \brief Application of the filter of detection of linear features
+ *
+ * This class implements the Tupin's detector D2 used to detect
+ * two parallel lines, based on a new edge detector.
+ *
+ * The region is devided in three zones to delimite two parallel lines.
+ * The size of one zone is defined by the product of the width
+ * of the linear feature by its length.
+ *
+ * For each vertical line, we calculate the cross correlation coefficient
+ * \rho_{ij} between regions i and j.
+ * \[\rho_{ij}=\fract{1}{1+(n_{i}+n_{j})\fract{n_{i}\gamma_{i}^{2}\bar{c}_{ij}^{2}}{n_{i}n_{j}(\bar{c}_{ij}-1)^{2}} \]
+ *
+ * where n_{i} is the pixel number in region i, \bar{c}_{ij} = \fract{\mu_{i}}{\mu_{j}}
+ * is the empiraical contrast between regions i and j, and \gamma_{i} is the varaition
+ * coefficient (ratio of standard deviation and mean) that measures homogeneity in
+ * radar imagery scenes.
+ *
+ * The intensity of detection in the three other directions \rho(\theta_{i})
+ * is determined by rotation of the pixels of each zone around the
+ * pixel central of the region considered. By default, the pixel location after
+ * rotation is determined by the Spline interpolator.
+ *
+ * Finally, the intensity of detection formed by the two parallel lines
+ * is determined by the minimum response of a ration edge detector on both sides
+ * of the linear structure:
+ * \[\rho = \min(\rho_{12};\rho_{13}) \]
+ * where \rho_{12} and \rho_{13} are the maximum response of the ratio edge
+ * detector of \rho(\theta_{i}).
+ *
+ * The exit is an image of intensity of detection.
+ *
+ *
+ */
+
+template <class TInputImage,
+ class TOutputImage,
+ class InterpolatorType = itk::BSplineInterpolateImageFunction<TInputImage> >
+class LineCorrelationDetector : public itk::ImageToImageFilter< TInputImage, TOutputImage >
+{
+public:
+ /** Extract dimensions as well of the images of entry of exit. */
+ itkStaticConstMacro( InputImageDimension,
+ unsigned int,
+ TInputImage::ImageDimension);
+ itkStaticConstMacro( OutputImageDimension,
+ unsigned int,
+ TOutputImage::ImageDimension);
+
+ typedef TInputImage InputImageType;
+ typedef TOutputImage OutputImageType;
+
+ /** typedef for the classes standards. */
+ typedef LineCorrelationDetector Self;
+ typedef itk::ImageToImageFilter< InputImageType, OutputImageType> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Method for management of the "object factory". */
+ itkNewMacro(Self);
+
+ /** Return the nale of the class. */
+ itkTypeMacro(LineCorrelationDetector, ImageToImageFilter);
+
+ /** Typedefs to describe and access Interpolator */
+ typedef typename InterpolatorType::Pointer InterpolatorPointer;
+ typedef typename InterpolatorType::CoordRepType CoordRepType;
+
+ typedef typename InputImageType::PointType TPoint;
+
+
+ /** Definition of the input and output images */
+ typedef typename InputImageType::PixelType InputPixelType;
+ typedef typename OutputImageType::PixelType OutputPixelType;
+
+
+ typedef typename InputImageType::RegionType InputImageRegionType;
+ typedef typename OutputImageType::RegionType OutputImageRegionType;
+
+ /** Definition of the size of the images. */
+ typedef typename InputImageType::SizeType SizeType;
+
+ /** Set the length of the linear feature. */
+ itkSetMacro(LengthLine, unsigned int);
+
+ /** Get the length of the linear feature. */
+ itkGetConstReferenceMacro(LengthLine, unsigned int);
+
+ /** Set the width of the linear feature. */
+ itkSetMacro(WidthLine, unsigned int);
+
+ /** Get the length of the linear feature. */
+ itkGetConstReferenceMacro(WidthLine, unsigned int);
+
+ /** Set the radius of one zone. */
+ itkSetMacro(Radius, SizeType);
+
+ /** Get the radius of one zone. */
+ itkGetConstReferenceMacro(Radius, SizeType);
+
+
+ virtual void GenerateInputRequestedRegion() throw(itk::InvalidRequestedRegionError);
+
+protected:
+ LineCorrelationDetector();
+ virtual ~LineCorrelationDetector() {};
+ void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+
+ /** LineCorrelationDetector can be implemented for a treatment of filter multithreaded.
+ * Thus, the ThreadedGenerateData() method is called for each thread process.
+ * The data image are allocated automatically by the mother class by calling the
+ * ThreadedGenerateData() method. ThreadedGenerateData can only write the portion
+ * of the image specified by the parameter "outputRegionForThread"
+ *
+ * \sa ImageToImageFilter::ThreadedGenerateData(),
+ * ImageToImageFilter::GenerateData() */
+ void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
+ int threadId );
+
+private:
+ LineCorrelationDetector(const Self&); //purposely not implemented
+ void operator=(const Self&); //purposely not implemented
+
+ /** Length of the linear feature = 2*m_LengthLine+1 */
+ unsigned int m_LengthLine;
+
+ /** Width of the linear feature = 2*m_WidthLine+1 */
+ unsigned int m_WidthLine;
+
+ /** Radius of the region*/
+ SizeType m_Radius;
+
+ /** Size of the facelist*/
+ SizeType m_FaceList;
+
+ InterpolatorPointer m_Interpolator;
+
+};
+} // end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbLineCorrelationDetector.txx"
+#endif
+
+
+#endif
diff --git a/Code/FeatureExtraction/otbLineCorrelationDetector.txx b/Code/FeatureExtraction/otbLineCorrelationDetector.txx
new file mode 100755
index 0000000000000000000000000000000000000000..a7cd1f10fcdb1b5f38f1c9b16f98452fd0513bec
--- /dev/null
+++ b/Code/FeatureExtraction/otbLineCorrelationDetector.txx
@@ -0,0 +1,367 @@
+/*=========================================================================
+
+ Programme : OTB (ORFEO ToolBox)
+ Auteurs : CS - C.Ruffel
+ Language : C++
+ Date : 14 mars 2006
+ Role : Filter of detection of linear features
+ $Id$
+
+=========================================================================*/
+#ifndef __otbLineCorrelationDetector_txx
+#define __otbLineCorrelationDetector_txx
+
+#include "otbLineCorrelationDetector.h"
+
+#include "itkDataObject.h"
+#include "itkExceptionObject.h"
+#include "itkConstNeighborhoodIterator.h"
+#include "itkNeighborhoodInnerProduct.h"
+#include "itkImageRegionIterator.h"
+#include "itkNeighborhoodAlgorithm.h"
+#include "itkZeroFluxNeumannBoundaryCondition.h"
+#include "itkProgressReporter.h"
+#include <math.h>
+
+#define M_PI 3.14159265358979323846
+
+namespace otb
+{
+
+/**
+ *
+ */
+template <class TInputImage, class TOutputImage, class InterpolatorType >
+LineCorrelationDetector<TInputImage, TOutputImage, InterpolatorType>::LineCorrelationDetector()
+{
+ m_Radius.Fill(1);
+ m_LengthLine = 1;
+ m_WidthLine = 0;
+ m_FaceList.Fill(0);
+ m_Interpolator = InterpolatorType::New();
+}
+
+template <class TInputImage, class TOutputImage, class InterpolatorType>
+void LineCorrelationDetector<TInputImage, TOutputImage, InterpolatorType>::GenerateInputRequestedRegion() throw (itk::InvalidRequestedRegionError)
+{
+ // 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();
+
+ // Define the size of the region by the radius
+ m_Radius[0] = static_cast<unsigned int>(3*m_WidthLine + 2);
+ m_Radius[1] = m_LengthLine ;
+
+ // Define the size of the facelist by taking into account the rotation of the region
+ m_FaceList[0] = static_cast<unsigned int>( sqrt((m_Radius[0]*m_Radius[0])+ (m_Radius[1]*m_Radius[1]))+1 );
+ m_FaceList[1] = m_FaceList[0];
+
+
+ // pad the input requested region by the operator radius
+ inputRequestedRegion.PadByRadius( m_Radius );
+
+ // 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 << static_cast<const char *>(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, class InterpolatorType>
+void LineCorrelationDetector< TInputImage, TOutputImage, InterpolatorType>
+::ThreadedGenerateData(
+ const OutputImageRegionType& outputRegionForThread,
+ int threadId
+ )
+{
+ unsigned int i;
+ itk::ZeroFluxNeumannBoundaryCondition<InputImageType> nbc;
+ itk::ConstNeighborhoodIterator<InputImageType> bit;
+ typename itk::ConstNeighborhoodIterator<InputImageType>::OffsetType off;
+ itk::ImageRegionIterator<OutputImageType> it;
+
+ // Allocate output
+ typename OutputImageType::Pointer output = this->GetOutput();
+ typename InputImageType::ConstPointer input = this->GetInput();
+
+ m_Interpolator->SetInputImage(input);
+
+ // Find the data-set boundary "faces"
+ typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>::FaceListType faceList;
+ typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator fit;
+
+
+ itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType> bC;
+ faceList = bC(input, outputRegionForThread, m_FaceList);
+
+
+ // support progress methods/callbacks
+ itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
+
+ typename TInputImage::IndexType bitIndex;
+ typename InterpolatorType::ContinuousIndexType Index;
+
+
+ // --------------------------------------------------------------------------
+
+ // Number of direction
+ const int NB_DIR = 4;
+ // Number of zone
+ const int NB_ZONE = 3;
+ // Definition of the 4 directions
+ double Theta[NB_DIR];
+
+ Theta[0] = 0.;
+ Theta[1] = M_PI / 4. ;
+ Theta[2] = M_PI / 2. ;
+ Theta[3] = 3*M_PI / 4. ;
+
+ // Number of the zone
+ unsigned int zone;
+
+ // Pixel numbers in each zone
+ const int NbPixelZone = (2*m_WidthLine+1)*(2*m_LengthLine+1);
+
+ double InvNbPixel = static_cast<double>(1. / NbPixelZone);
+
+ // Contains for the 4 directions the sum of the pixels belonging to each zone
+ double Sum[NB_DIR][NB_ZONE];
+ double Sum2[NB_DIR][NB_ZONE];
+
+ // Mean of zone 1, 2 and 3
+ double M1, M2, M3;
+
+ // Standard deviation
+ double sigma1, sigma2, sigma3;
+
+ // Result of the filter for each direction
+ double rho12_theta[NB_DIR];
+ double rho13_theta[NB_DIR];
+ // Result in one direction
+ double rho12, rho13;
+ // Intensity of the linear feature
+ double rho;
+
+ // Pixel location in the input image
+ int X, Y;
+ // Value of the pixel
+ double ValXY = 0.;
+
+ // Pixel location after rotation in the system axis of the region
+ double xout, yout;
+
+ // Pixel location in the input image after rotation
+ TPoint point;
+
+ // location of the pixel central in the input image
+ int Xc, Yc;
+
+ // location of the pixel central between zone 1 and 2 and between zone 1 and 3
+ int Xc12, Xc13;
+
+ //---------------------------------------------------------------------------
+
+ // Process each of the boundary faces. These are N-d regions which border
+ // the edge of the buffer.
+ for (fit=faceList.begin(); fit != faceList.end(); ++fit)
+ {
+ bit = itk::ConstNeighborhoodIterator<InputImageType>(m_Radius, input, *fit);
+ unsigned int neighborhoodSize = bit.Size();
+
+ it = itk::ImageRegionIterator<OutputImageType>(output, *fit);
+
+ bit.OverrideBoundaryCondition(&nbc);
+ bit.GoToBegin();
+
+
+ while ( ! bit.IsAtEnd() )
+ {
+//std::cout << "Xc,Yc " << bit.GetIndex() << std::endl;
+
+ // Initialisations
+ for (int dir=0; dir<NB_DIR; dir++)
+ {
+ for (int z=0; z<NB_ZONE; z++)
+ {
+ Sum[dir][z] = 0.;
+ Sum2[dir][z] = 0.;
+ }
+ }
+
+
+ // Location of the central pixel of the region
+ bitIndex = bit.GetIndex();
+
+ Xc = bitIndex[0];
+ Yc = bitIndex[1];
+
+ // Location of the central pixel between zone 1 and zone 2
+ Xc12 = Xc - m_WidthLine - 1;
+
+ // Location of the central pixel between zone 1 and zone 3
+ Xc13 = Xc + m_WidthLine + 1;
+
+ // Loop on the region
+ for (i = 0; i < neighborhoodSize; ++i)
+ {
+//std::cout << "---"<< i <<" "<< bit.GetIndex(i)<< std::endl;
+//std::cout << "val(X,Y) "<< static_cast<double>(bit.GetPixel(i)) << " " << std::endl;
+
+ bitIndex = bit.GetIndex(i);
+ X = bitIndex[0];
+ Y = bitIndex[1];
+
+
+ // We determine in the vertical direction with which zone the pixel belongs.
+
+ if ( X < Xc12 )
+ zone = 1;
+ else if ( ( Xc12 < X ) && ( X < Xc13 ) )
+ zone = 0;
+ else if ( X < Xc13 )
+ zone = 2;
+
+ ValXY = static_cast<double>(bit.GetPixel(i));
+
+ Sum[0][zone] += ValXY;
+ Sum2[0][zone] += ValXY*ValXY;
+
+ // Loop on the 3 other directions
+ for ( int dir=1; dir<NB_DIR; dir++ )
+ {
+ ROTATION( (X-Xc), (Y-Yc), Theta[dir], xout, yout);
+
+ Index[0] = static_cast<CoordRepType>(xout + Xc);
+ Index[1] = static_cast<CoordRepType>(yout + Yc);
+
+//std::cout << "X' Y' "<< (xout + Xc) << " " << (yout + Yc) << std::endl;
+//std::cout << "val(X',Y') "<< static_cast<double>(m_Interpolator->EvaluateAtContinuousIndex( Index )) << std::endl;
+
+ ValXY = static_cast<double>(m_Interpolator->EvaluateAtContinuousIndex( Index ));
+
+ Sum[dir][zone] += ValXY;
+ Sum2[dir][zone] += ValXY*ValXY;
+
+ }
+
+ } // end of the loop on the pixels of the region
+
+ rho12 = -1.;
+ rho13 = -1.;
+
+ double NbPixel = static_cast<double>(NbPixelZone);
+
+ // Loop on the 4 directions
+ for ( int dir=0; dir<NB_DIR; dir++ )
+ {
+
+ // Calculation of the mean for the 3 zones
+ M1 = Sum[dir][0] * InvNbPixel;
+ M2 = Sum[dir][1] * InvNbPixel;
+ M3 = Sum[dir][2] * InvNbPixel;
+
+ // Calculation of the standard deviation for the 3 zones
+ sigma1 = sqrt( (Sum2[dir][0] - NbPixel*M1*M1) * InvNbPixel );
+ sigma2 = sqrt( (Sum2[dir][1] - NbPixel*M2*M2) * InvNbPixel );
+ sigma3 = sqrt( (Sum2[dir][2] - NbPixel*M3*M3) * InvNbPixel );
+
+ // Calculation of the cross correlation coefficient
+
+ double d1 = 0.;
+ double d2 = 0.;
+ double d3 = 0.;
+
+ // rho12
+ if ( M2 != 0. )
+ d1 = NbPixel*pow(sigma1,2)*pow(M1/M2,2);
+ d2 = NbPixel*pow(sigma2,2);
+ if ( M2 != 0. )
+ d3 = NbPixel*NbPixel*pow(((M1/M2)-1.),2);
+
+ if ( ( M1 != M2 ) )
+ rho12_theta[dir] = static_cast<double>( 1. / ( 1. + ( 2.*NbPixel*(d1+d2)/d3 ) ) );
+ else
+ rho12_theta[dir] = 0.;
+
+ // rho13
+ if ( M3 != 0. )
+ d1 = NbPixel*pow(sigma1,2)*pow(M1/M3,2);
+ d2 = NbPixel*pow(sigma3,2);
+ if ( M3 != 0. )
+ d3 = NbPixel*NbPixel*pow(((M1/M3)-1.),2);
+
+ if ( ( M1 != M3 ) )
+ rho13_theta[dir] = static_cast<double>( 1. / ( 1. + ( 2.*NbPixel*(d1+d2)/d3 ) ) );
+ else
+ rho13_theta[dir] = 0.;
+
+
+ // Determination of the maximum cross correlation coefficients
+ rho12 = static_cast<double>( MAX( rho12, rho12_theta[dir] ) );
+ rho13 = static_cast<double>( MAX( rho13, rho13_theta[dir] ) );
+
+ } // end of the loop on the directions
+
+ // Intensity of the linear feature
+ rho = MIN ( rho12, rho13 );
+
+ // Assignment of this value to the output pixel
+ it.Set( static_cast<OutputPixelType>(rho) );
+
+ ++bit;
+ ++it;
+ progress.CompletedPixel();
+
+ }
+
+ }
+}
+
+/**
+ * Standard "PrintSelf" method
+ */
+template <class TInputImage, class TOutput, class InterpolatorType>
+void
+LineCorrelationDetector<TInputImage, TOutput, InterpolatorType>::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+ Superclass::PrintSelf( os, indent );
+ os << indent << "Radius: " << m_Radius << std::endl;
+}
+
+
+} // end namespace otb
+
+
+#endif