Filtre produit scalaire. A terminer.

Code/FeatureExtraction/otbNeighborhoodScalarProductFilter.h

+/*=========================================================================
+
+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 _otbNeighborhoodScalarProductFilter_h
+#define _otbNeighborhoodScalarProductFilter_h
+
+#include "otbImageToModulusAndDirectionImageFilter.h"
+
+namespace otb
+{
+/** \class NeighborhoodScalarProductFilter
+ *  \brief 
+ *
+ *  \ingroup Streamed
+ *  \ingroup Threaded
+ */
+template <class TInputImage, class TOutputModulus, class TOutputDirection>
+class ITK_EXPORT NeighborhoodScalarProductFilter
+  : public ImageToModulusAndDirectionImageFilter<TInputImage,TOutputModulus,TOutputDirection>
+{
+ public:
+  /** Standard typedefs */
+  typedef NeighborhoodScalarProductFilter                                                    Self;
+  typedef ImageToModulusAndDirectionImageFilter<TInputImage,TOutputModulus,TOutputDirection> Superclass;
+  typedef itk::SmartPointer<Self>                                                            Pointer;
+  typedef itk::SmartPointer<const Self>                                                      ConstPointer;
+  
+  /** Type macro */
+  itkNewMacro(Self);
+  
+  /** Creation through object factory macro */
+  itkTypeMacro(NeighborhoodScalarProductFilter,ImageToModulusAndDirectionImageFilter);
+  
+  /** Template parameters typedefs */
+  typedef TInputImage InputImageType;
+  typedef typename InputImageType::ConstPointer InputImageConstPointerType;
+  typedef typename InputImageType::PixelType InputPixelType;
+  typedef TOutputModulus OutputModulusType;
+  typedef typename OutputModulusType::Pointer OutputModulusPointerType;
+  typedef typename OutputModulusType::RegionType RegionType;
+  typedef typename OutputModulusType::SizeType SizeType;
+  typedef typename OutputModulusType::IndexType IndexType;
+  typedef TOutputDirection OutputDirectionType;
+  typedef typename OutputDirectionType::Pointer OutputDirectionPointerType;
+  typedef typename OutputDirectionType::RegionType OutputImageRegionType;
+ 
+protected:
+  /** Constructor */
+  NeighborhoodScalarProductFilter();
+  /** Destructor */
+  virtual ~NeighborhoodScalarProductFilter() {};
+ /**PrintSelf method */
+  virtual void PrintSelf(std::ostream& os, itk::Indent indent) const;
+  /** NeighborhoodScalarProductImageFilter can be implemented as a multithreaded filter.
+   * Therefore, this implementation provides a ThreadedGenerateData() routine
+   * which is called for each processing thread. The output image data is
+   * allocated automatically by the superclass prior to calling
+   * ThreadedGenerateData().  ThreadedGenerateData can only write to the
+   * portion of the output image specified by the parameter
+   * "outputRegionForThread"
+   *
+   * \sa ImageToImageFilter::ThreadedGenerateData(),
+   *     ImageToImageFilter::GenerateData()  */
+  void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
+                            int threadId );
+private:
+  NeighborhoodScalarProductFilter(const Self&); //purposely not implemented
+  void operator=(const Self&); //purposely not implemented
+};
+}// End namespace otb
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbNeighborhoodScalarProductFilter.txx"
+#endif
+
+#endif

Code/FeatureExtraction/otbNeighborhoodScalarProductFilter.txx

+/*=========================================================================
+
+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 _otbNeighborhoodScalarProduct_txx
+#define _otbNeighborhoodScalarProduct_txx
+
+#include "otbNeighborhoodScalarProductFilter.h"
+#include "itkImageRegionIterator.h"
+#include "itkConstNeighborhoodIterator.h"
+#include "itkNeighborhoodAlgorithm.h"
+#include "itkProgressReporter.h"
+
+namespace otb
+{
+/**
+ * Constructor
+ */
+template <class TInputImage, class TOutputModulus, class TOutputDirection>
+NeighborhoodScalarProductFilter<TInputImage,TOutputModulus,TOutputDirection>
+::NeighborhoodScalarProductFilter()
+{}
+
+template <class TInputImage, class TOutputModulus, class TOutputDirection>
+void
+NeighborhoodScalarProductFilter<TInputImage,TOutputModulus,TOutputDirection>
+::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
+{
+  // some typedefs
+  typedef itk::ConstNeighborhoodIterator<InputImageType> NeighborhoodIteratorType;
+  typedef typename NeighborhoodIteratorType::RadiusType RadiusType;
+  typedef typename NeighborhoodIteratorType::OffsetType OffsetType;
+  typedef itk::ZeroFluxNeumannBoundaryCondition<InputImageType> BoundaryConditionType ;
+  typedef itk::ImageRegionIterator<OutputModulusType> OutputIteratorType;
+  typedef itk::ImageRegionIterator<OutputDirectionType> OutputDirectionIteratorType;
+  typedef itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType> BoundaryFacesCalculatorType;
+  typedef typename BoundaryFacesCalculatorType::FaceListType FaceListType;
+  typedef typename FaceListType::iterator FaceListIteratorType;
+  
+  // Pointers on inputs/outputs 
+  InputImageType * inputPtr = const_cast<InputImageType *>(this->GetInput());
+  OutputModulusPointerType outputPtr = this->GetOutput();
+  OutputDirectionPointerType outputDirPtr = this->GetOutputDirection();
+  
+  // Neighborhood radius
+  RadiusType r;
+  r.Fill(1);
+
+  // Find the data-set boundary "faces"
+  BoundaryFacesCalculatorType bC;
+  FaceListType faceList = bC(inputPtr, outputRegionForThread, r);
+  FaceListIteratorType 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.
+  for (fit=faceList.begin(); fit != faceList.end(); ++fit)
+    { 
+      NeighborhoodIteratorType neighInputIt(r, inputPtr, *fit);
+      OutputIteratorType outputIt(outputPtr,outputRegionForThread);
+      OutputDirectionIteratorType outputDirIt(outputDirPtr,outputRegionForThread);
+      neighInputIt.GoToBegin();
+      outputIt.GoToBegin();
+      outputDirIt.GoToBegin();
+    
+      while ((!neighInputIt.IsAtEnd()) && (!outputIt.IsAtEnd()) && (!outputDirIt.IsAtEnd()) )
+	{
+	  // local variable intialisation
+	  int neighborhoodNumberMax = 0;
+	  double scalarMaxValue= 0.0;
+	  int flagPosNegDirection = 0;
+	 
+	  // walk through each case
+	  for (int neighborhoodNumber = 0;  neighborhoodNumber<4; ++neighborhoodNumber)
+	    {
+	    double scalarCurrentValue = 0.0;
+	    OffsetType offset1;
+	    OffsetType offset2;
+	    switch(neighborhoodNumber){
+	    case 0:
+	      offset1[0]=1;
+	      offset1[1]=-1;
+	      offset2[0]=-1;
+	      offset2[1]=1;
+	      break;
+	    case 1:
+	      offset1[0]=1;
+	      offset1[1]=0;
+	      offset2[0]=-1;
+	      offset2[1]=0;
+	      break;
+	    case 2:
+	      offset1[0]=1;
+	      offset1[1]=1;
+	      offset2[0]=-1;
+	      offset2[1]=-1;
+	      break;
+	    case 3:
+	      offset1[0]=0;
+	      offset1[1]=1;
+	      offset2[0]=0;
+	      offset2[1]=-1;
+	      break;                
+	    }  
+	    // Get the gradient values
+	    InputPixelType pixel1 = neighInputIt.GetPixel(offset1);
+	    InputPixelType pixel2 = neighInputIt.GetPixel(offset2);
+	    
+	    // Compute the scalar product
+	    scalarCurrentValue = - (pixel1[0]*pixel2[0]+pixel1[1]*pixel2[1]);
+	    
+	    // If the value is upper than the current max value
+	    if (scalarCurrentValue > scalarMaxValue)
+	      {
+		// keep this configuration
+		scalarMaxValue = scalarCurrentValue;
+		neighborhoodNumberMax = neighborhoodNumber;
+		
+		// Also keep the direction
+		if (pixel1[0] <0) 
+		  {
+		    flagPosNegDirection = 1;
+		  } 
+		else 
+		  {
+		  flagPosNegDirection = 0;
+		  }
+		
+	      }
+	    }
+	  // Compute the direction 
+	  double angle = (1+neighborhoodNumberMax) * M_PI/4;
+	  if (flagPosNegDirection) 
+	    {
+	      angle -= M_PI;
+	    }
+
+	  // Set the ouptut values	  
+	  outputIt.Set(scalarMaxValue);
+	  outputDirIt.Set(angle);
+	  ++neighInputIt;
+	  ++outputIt;
+	  ++outputDirIt;
+	  progress.CompletedPixel();
+	}
+    }
+}
+/**
+ * PrintSelf Method
+ */
+template <class TInputImage, class TOutputModulus, class TOutputDirection>
+void
+NeighborhoodScalarProductFilter<TInputImage,TOutputModulus,TOutputDirection>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+  Superclass::PrintSelf(os, indent);
+}
+} // End namespace otb
+#endif

Testing/Code/FeatureExtraction/otbNeighborhoodScalarProductFilter.cxx

+/*=========================================================================
+
+  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 "otbNeighborhoodScalarProductFilter.h"
+#include "otbImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+#include "itkGradientRecursiveGaussianImageFilter.h"
+
+int otbNeighborhoodScalarProductFilter(int argc, char * argv[])
+{
+  try
+    {
+      const char * infname = argv[1];
+      const char * outfname = argv[2];
+      const char * diroutfname = argv[3];
+      const double sigma = atof(argv[4]);
+
+      const unsigned int Dimension = 2;
+      typedef double PixelType;
+      typedef itk::CovariantVector<PixelType,Dimension> VectorPixelType;
+      typedef otb::Image<VectorPixelType,Dimension> VectorImageType;
+      typedef otb::Image<PixelType,Dimension> ImageType;
+      typedef otb::NeighborhoodScalarProductFilter<VectorImageType,ImageType,ImageType> FilterType;
+      typedef otb::ImageFileReader<ImageType> ReaderType;
+      typedef otb::ImageFileWriter<ImageType> WriterType;
+      typedef itk::GradientRecursiveGaussianImageFilter<ImageType,VectorImageType> GradientFilterType;
+      
+
+      // Instantiating object
+      FilterType::Pointer filter = FilterType::New();
+      ReaderType::Pointer reader = ReaderType::New();
+      WriterType::Pointer writer = WriterType::New();
+      GradientFilterType::Pointer gradient = GradientFilterType::New();
+      
+      reader->SetFileName(infname);
+      gradient->SetInput(reader->GetOutput());
+      gradient->SetSigma(sigma);
+      filter->SetInput(gradient->GetOutput());
+      writer->SetInput(filter->GetOutput());
+      writer->SetFileName(outfname);
+      writer->Update();
+      writer = WriterType::New();
+      writer->SetFileName(diroutfname);
+      writer->SetInput(filter->GetOutputDirection());
+      writer->Update();
+    }
+
+  catch( itk::ExceptionObject & err ) 
+    { 
+    std::cout << "Exception itk::ExceptionObject thrown !" << std::endl; 
+    std::cout << err << std::endl; 
+    return EXIT_FAILURE;
+    } 
+
+  catch( ... ) 
+    { 
+    std::cout << "Unknown exception thrown !" << std::endl; 
+    return EXIT_FAILURE;
+    } 
+  return EXIT_SUCCESS;
+}

Testing/Code/FeatureExtraction/otbNeighborhoodScalarProductFilterNew.cxx

+/*=========================================================================
+
+  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 "otbNeighborhoodScalarProductFilter.h"
+#include "otbImage.h"
+
+int otbNeighborhoodScalarProductFilterNew(int argc, char * argv[])
+{
+  try
+    {
+      const unsigned int Dimension = 2;
+      typedef double PixelType;
+      typedef itk::CovariantVector<PixelType,Dimension> VectorPixelType;
+      typedef otb::Image<VectorPixelType,Dimension> InputImageType;
+      typedef otb::Image<PixelType,Dimension> OutputImageType;
+      typedef otb::NeighborhoodScalarProductFilter<InputImageType,OutputImageType,OutputImageType> FilterType;
+
+      // Instantiating object
+      FilterType::Pointer object = FilterType::New();
+    }
+
+  catch( itk::ExceptionObject & err ) 
+    { 
+    std::cout << "Exception itk::ExceptionObject thrown !" << std::endl; 
+    std::cout << err << std::endl; 
+    return EXIT_FAILURE;
+    } 
+
+  catch( ... ) 
+    { 
+    std::cout << "Unknown exception thrown !" << std::endl; 
+    return EXIT_FAILURE;
+    } 
+  return EXIT_SUCCESS;
+}