nomsg

Code/FeatureExtraction/otbLineRatioDetectorImageFilter.h

+/*=========================================================================
+
+  Programme :   OTB (ORFEO ToolBox)
+  Auteurs   :   CS - C.Ruffel
+  Language  :   C++
+  Date      :   14 mars 2006
+  Role      :   Filter of detection of linear features 
+  $Id$ 
+
+=========================================================================*/
+#ifndef __otbLineRatioDetectorImageFilter_h
+#define __otbLineRatioDetectorImageFilter_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 LineRatioDetectorImageFilter
+ * \brief Application of the filter of detection of linear features 
+ *
+ * This class implements the Tupin's detector D1 used to detect 
+ * two parallel lines. This detector is derived from the coupling of two
+ * ratio edge detectors (Touzi detector) on both side of a region.
+ *
+ * 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 intensity of detection 
+ * R_{12}(\theta_{0}) between zones 1 and 2 and R_{13}(\theta_{0}) between  
+ * zones 1 and 3 according to the principle of the Touzi's filter.
+ *
+ * The response of the edge detector between two zones i and j is:
+ *	\[R_{ij}=1-\min (\fract{\mu_{i}}{\mu_{j}};\fract{\mu_{j}}{\mu_{i}}) \]
+ *
+ * The intensity of detection in the three other directions R(\theta_{i})
+ * is determined by rotation of the pixels of each zone around the 
+ * central pixel 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:
+ *	\[R = \min (R_{12};R_{13}) \]
+ * where R_{12} and R_{13} are the maximum response of the ratio edge 
+ * detector of R(\theta_{i}). The intensity of detection lies in
+ * the interval [0, 1].
+ * 
+ * The output is an image of intensity of detection. 
+ * 
+ */
+
+template <class TInputImage, 
+	  class TOutputImage,
+	  class InterpolatorType = itk::BSplineInterpolateImageFunction<TInputImage> >
+class ITK_EXPORT LineRatioDetectorImageFilter :  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 LineRatioDetectorImageFilter 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 name of the class. */
+  itkTypeMacro(LineRatioDetectorImageFilter, itk::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);
+  
+  const OutputImageType * GetOutputDirection();
+
+protected:
+  LineRatioDetectorImageFilter();
+  virtual ~LineRatioDetectorImageFilter() {};
+  void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+  void BeforeThreadedGenerateData();
+
+  /** LineRatioDetectorImageFilter 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:
+  LineRatioDetectorImageFilter(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;
+
+  typename OutputImageType::Pointer m_DirectionOuputImage;
+
+};
+} // end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbLineRatioDetectorImageFilter.txx"
+#endif
+
+  
+#endif

Code/FeatureExtraction/otbLineRatioDetectorImageFilter.txx

+/*=========================================================================
+
+  Programme :   OTB (ORFEO ToolBox)
+  Auteurs   :   CS - C.Ruffel
+  Language  :   C++
+  Date      :   14 mars 2006
+  Role      :   Filter of detection of linear features 
+  $Id$ 
+
+=========================================================================*/
+#ifndef __otbLineRatioDetectorImageFilter_txx
+#define __otbLineRatioDetectorImageFilter_txx
+
+#include "otbLineRatioDetectorImageFilter.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 >
+LineRatioDetectorImageFilter<TInputImage, TOutputImage, InterpolatorType>::LineRatioDetectorImageFilter()
+{
+  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 LineRatioDetectorImageFilter<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;
+    }
+}
+
+/**
+ * Set up state of filter before multi-threading.
+ * InterpolatorType::SetInputImage is not thread-safe and hence
+ * has to be set up before ThreadedGenerateData
+ */
+template <class TInputImage, class TOutputImage, class InterpolatorType>
+void 
+LineRatioDetectorImageFilter< TInputImage, TOutputImage, InterpolatorType>
+::BeforeThreadedGenerateData()
+{
+
+  typename OutputImageType::RegionType  region;    
+  typename OutputImageType::Pointer     output = this->GetOutput();
+
+  m_DirectionOuputImage = OutputImageType::New();
+
+  region.SetSize(output->GetLargestPossibleRegion().GetSize());
+  region.SetIndex(output->GetLargestPossibleRegion().GetIndex());
+  m_DirectionOuputImage->SetRegions( region );
+  m_DirectionOuputImage->SetOrigin(output->GetOrigin());
+  m_DirectionOuputImage->SetSpacing(output->GetSpacing());
+  m_DirectionOuputImage->Allocate();
+
+}
+
+template <class TInputImage, class TOutputImage, class InterpolatorType>
+const typename LineRatioDetectorImageFilter< TInputImage, TOutputImage, InterpolatorType>::OutputImageType *
+LineRatioDetectorImageFilter< TInputImage, TOutputImage, InterpolatorType>
+::GetOutputDirection()
+{
+	this->Update();
+	return 	static_cast< const OutputImageType *> (m_DirectionOuputImage);
+}
+
+template< class TInputImage, class TOutputImage, class InterpolatorType>
+void LineRatioDetectorImageFilter< 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;
+  itk::ImageRegionIterator<OutputImageType> 			itdir;
+  
+  // Allocate output
+  typename InputImageType::ConstPointer input  = this->GetInput();
+  typename OutputImageType::Pointer     output = this->GetOutput();  
+  typename OutputImageType::Pointer     outputDir = m_DirectionOuputImage;
+  
+  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 number in each zone
+  const int NbPixelZone = (2*m_WidthLine+1)*(2*m_LengthLine+1);
+  
+  double InvNbPixelZone = 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];
+  // Mean of zone 1, 2 and 3
+  double M1, M2, M3;
+  
+  // Result of the filter for each direction
+  double R12_theta[NB_DIR];
+  double R13_theta[NB_DIR];
+  // Result for each direction
+  double R_theta[NB_DIR];
+  double Sum_R_theta = 0.;
+  // Intensity of the linear feature
+  double R;
+  
+  // Direction of detection
+  double Direction = 0.; 
+
+  // Pixel location in the input image
+  int X, Y;
+  
+  // 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 central pixel in the input image
+  int Xc, Yc;
+
+  // location of the central pixel 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);
+    itdir = itk::ImageRegionIterator<OutputImageType>(outputDir, *fit);
+    
+    bit.OverrideBoundaryCondition(&nbc);
+    bit.GoToBegin();
+
+
+    while ( ! bit.IsAtEnd() )
+      {
+//std::cout << "Xc,Yc " << bit.GetIndex() << std::endl;
+
+      // Initializations
+      for (unsigned int dir=0; dir<NB_DIR; dir++)
+        {
+        for (unsigned int z=0; z<NB_ZONE; z++) 
+          Sum[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;
+        else
+          continue;
+    
+        Sum[0][zone] += static_cast<double>(bit.GetPixel(i));
+         
+        // Loop on the 3 other directions
+        for (unsigned 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;
+           
+            Sum[dir][zone] += static_cast<double>(m_Interpolator->EvaluateAtContinuousIndex( Index ));
+          
+          }      
+
+        } // end of the loop on the pixels of the region 
+          
+      R = -1.;
+      Direction = 0.;
+      Sum_R_theta = 0.;
+           
+      // Loop on the 4 directions
+      for (unsigned int dir=0; dir<NB_DIR; dir++ )
+        {
+        		
+        // Calculation of the mean for the 3 zones	
+        M1 = Sum[dir][0] * InvNbPixelZone;
+        M2 = Sum[dir][1] * InvNbPixelZone;
+        M3 = Sum[dir][2] * InvNbPixelZone;
+     	
+        // Calculation of the intensity of detection
+        if (( M1 != 0. ) && (M2 != 0. )) 
+          R12_theta[dir] = static_cast<double>( 1 - MIN( (M1/M2), (M2/M1) ) );
+        else
+	  R12_theta[dir] = 0.;
+	  
+	if (( M1 != 0. ) && (M3 != 0. )) 
+          R13_theta[dir] = static_cast<double>( 1 - MIN( (M1/M3), (M3/M1) ) );
+        else
+	  R13_theta[dir] = 0.;
+	  
+	// Determination of the minimum intensity of detection between R12 et R13
+	R_theta[dir] = static_cast<double>( MIN( R12_theta[dir], R13_theta[dir] ) );
+	
+	R = static_cast<double>( MAX ( R, R_theta[dir] ) );
+                
+        // Calculation of the directions
+        Direction += Theta[dir] * R_theta[dir];        
+        Sum_R_theta += R_theta[dir];
+	  
+      		
+        } // end of the loop on the directions
+      
+      // Determination of the direction of the contour
+      if ( Sum_R_theta != 0. )	
+        Direction = Direction / Sum_R_theta;
+
+      // Assignment of this value to the output pixel
+      it.Set( static_cast<OutputPixelType>(R) );
+      
+      // Assignment of this value to the "outputdir" pixel
+      itdir.Set( static_cast<OutputPixelType>(Direction) );  
+
+      ++bit;
+      ++it;
+      ++itdir;
+      progress.CompletedPixel();  
+      
+      }
+      
+    }
+}
+
+/**
+ * Standard "PrintSelf" method
+ */
+template <class TInputImage, class TOutput, class InterpolatorType>
+void 
+LineRatioDetectorImageFilter<TInputImage, TOutput, InterpolatorType>::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+  Superclass::PrintSelf( os, indent );
+  os << indent << "Length: " << m_LengthLine << std::endl;
+  os << indent << "Width: " << m_WidthLine << std::endl;
+  
+}
+
+
+} // end namespace otb
+
+
+#endif