ADD : PointSetDensity Function & related tests

Code/BasicFilters/otbPointSetDensityFunction.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 __otbPointSetDensityFunction_h
+#define __otbPointSetDensityFunction_h
+
+#include "otbPointSetFunction.h"
+#include "itkPoint.h"
+#include "itkProcessObject.h"
+
+namespace otb
+{
+
+/**
+ * \class PointSetDensityFunction
+ * \brief Calculate the density in the neighborhood of a pixel
+ *
+ * \ingroup PointSetFunctions
+ */
+template <class TPointSet, class  TOutput>
+ class ITK_EXPORT PointSetDensityFunction : public PointSetFunction< TPointSet , TOutput >
+{
+public:
+  /** Standard class typedefs. */
+  typedef PointSetDensityFunction                    Self;
+  typedef PointSetFunction< TPointSet ,TOutput >     Superclass;
+  typedef itk::SmartPointer<Self>                    Pointer;
+  typedef itk::SmartPointer<const Self>              ConstPointer;
+
+  
+  /** Run-time type information (and related methods). */
+  itkTypeMacro(PointSetDensityFunction, PointSetFunction);
+  
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self); 
+
+  /** PointSet Type typedef Support*/
+  typedef TPointSet                            PointSetType;
+  typedef typename Superclass::InputType       InputType;
+  typedef typename  PointSetType::Pointer      PointSetPointerType;
+  
+  /** TOutput typedef suppoty*/
+  typedef TOutput                              OutputType;
+
+  /** Set/Get the number of scales*/
+  itkSetMacro(Radius,unsigned int);
+  itkGetMacro(Radius,unsigned int);
+  
+  /** Evaluate Method */
+  virtual OutputType Evaluate(const InputType& input ) const;
+ 
+protected:
+  PointSetDensityFunction();
+  ~PointSetDensityFunction(){};
+
+  void PrintSelf(std::ostream& os, itk::Indent indent) const;
+  
+
+private:
+  PointSetDensityFunction( const Self& ); //purposely not implemented
+  void operator=( const Self& ); //purposely not implemented
+ 
+  unsigned int m_Radius;
+};
+
+} // end namespace otb
+
+
+#ifndef OTB_MANUAL_INSTANTIATION 
+#include "otbPointSetDensityFunction.txx"
+#endif
+
+#endif
+

Code/BasicFilters/otbPointSetDensityFunction.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 __otbPointSetDensityFunction_txx
+#define __otbPointSetDensityFunction_txx
+
+#include "otbPointSetDensityFunction.h"
+#include "otbImage.h"
+#include "otbMath.h"
+
+
+namespace otb
+{
+
+/**
+ * Constructor
+ */
+template <class TPointSet, class  TOutput >
+PointSetDensityFunction< TPointSet,   TOutput>
+::PointSetDensityFunction()
+{
+  m_Radius = 1;
+}
+
+
+/**
+ *
+ */
+template <class TPointSet, class  TOutput > 
+typename PointSetDensityFunction< TPointSet,   TOutput>
+::OutputType
+PointSetDensityFunction< TPointSet,   TOutput>
+::Evaluate(const  InputType& input ) const
+{
+      
+  typename otb::Image<OutputType,2>::IndexType  index;
+  index[0] = input[0];
+  index[1] = input[1];
+  
+  int accu = 0;
+  double surface = M_PI*vcl_pow(2.,static_cast<double>(m_Radius));
+  
+  if(this->GetPointSet()->GetNumberOfPoints() != 0)
+    {
+      typedef typename TPointSet::PointsContainer::ConstIterator     iteratorType;
+      iteratorType it = this->GetPointSet()->GetPoints()->Begin();
+      
+      while( it != this->GetPointSet()->GetPoints()->End())
+	{
+	  float distX2 =( index[0]-it.Value()[0])*( index[0]-it.Value()[0]);
+	  float distY2 =( index[1]-it.Value()[1])*( index[1]-it.Value()[1]);
+	  float dist = vcl_sqrt(distX2 + distY2);
+	  
+	  if(dist <= m_Radius)
+	    accu++;
+	  
+	  ++it;
+	}
+    }
+  else
+    return 0.;
+  
+  return static_cast<float>(accu/surface);
+}
+
+
+/**
+ *
+ */
+template <class TPointSet, class  TOutput > 
+void
+PointSetDensityFunction< TPointSet,   TOutput>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+  this->Superclass::PrintSelf(os,indent);
+}
+
+
+} // end namespace otb
+
+#endif

Code/BasicFilters/otbPointSetFunction.h

@@ -33,12 +33,12 @@ namespace otb
  */
 template <class TPointSet, class  TOutput>
  class ITK_EXPORT PointSetFunction :
-  public itk::SpatialFunction<TOutput >
+  public itk::SpatialFunction< TOutput , 2/* TODO : change 2 by PointType::PointDimension*/, typename TPointSet::PointType >
 {
 public:
   /** Standard class typedefs. */
-  typedef PointSetFunction                    Self;
-  typedef itk::SpatialFunction< TOutput >             Superclass;
+typedef PointSetFunction                                       Self;
+ typedef itk::SpatialFunction< TOutput, 2 ,  typename TPointSet::PointType >       Superclass;
     
   /** Run-time type information (and related methods). */
   itkTypeMacro(PointSetFunction, itk::SpatialFunction);
@@ -51,9 +51,12 @@ public:
   typedef TOutput           OutputType;
   
   /** Set the input image (reimplemented since we need to set the detector input) */
-  virtual void SetPointSet(PointSetType* PointSet);
-  virtual PointSetType * GetPointSet();
+  itkGetConstObjectMacro(PointSet,PointSetType);
 
+  void SetPointSet( PointSetType* PointSet)
+    {
+      m_PointSet = PointSet;
+    }
  
 protected:
   PointSetFunction();

Code/BasicFilters/otbPointSetFunction.txx

@@ -49,28 +49,29 @@ PointSetFunction< TPointSet,   TOutput>
  }
 
 
-/**
- * SetDetector method
- */
-template <class TPointSet, class  TOutput > 
-void
-PointSetFunction< TPointSet,   TOutput>
-::SetPointSet(PointSetType* PointSet)
-{
-  m_PointSet = PointSet;
-}
+// /**
+//  * SetDetector method
+//  */
+// template <class TPointSet, class  TOutput > 
+// void
+// PointSetFunction< TPointSet,   TOutput>
+// ::SetPointSet(PointSetType* PointSet) const
+// {
+//   m_PointSet = PointSet;
+// }
 
-/**
- * GetDetector method
- */
-template <class TPointSet, class  TOutput > 
-typename PointSetFunction< TPointSet,TOutput>
-::PointSetType *
-PointSetFunction< TPointSet,   TOutput>
-::GetPointSet() 
-{
-  return m_PointSet;
-}
+// /**
+//  * GetDetector method
+//  */
+// template <class TPointSet, class  TOutput > 
+// const
+// typename PointSetFunction< TPointSet,TOutput>
+// ::PointSetType *
+// PointSetFunction< TPointSet,   TOutput>
+// ::GetPointSet() const
+// {
+//   return m_PointSet;
+// }
 } // end namespace otb
 
 #endif

Testing/Code/BasicFilters/CMakeLists.txt

@@ -64,6 +64,7 @@ ADD_TEST(bfTvFiltreFrost ${BASICFILTERS_TESTS1}
 ADD_TEST(bfTuImageToPointSetFilterTest ${BASICFILTERS_TESTS1}
          otbImageToPointSetFilterTest)
 
+
 # -------            otb::OpeningClosingMorphologicalFilter   ------------------------------
 
 ADD_TEST(bfTuOpeningClosingFilterNew ${BASICFILTERS_TESTS1}
@@ -1120,6 +1121,30 @@ ADD_TEST(bfTuPointSetFunctionNew ${BASICFILTERS_TESTS11}
 	 otbPointSetFunctionNew
 	 )
 
+
+
+#------------ otbPointSetDensityFunction ---------------------
+
+ADD_TEST(bfTuPointSetDensityFunctionNew ${BASICFILTERS_TESTS11}
+         otbPointSetDensityFunctionNew)
+
+
+ADD_TEST(bfTvPointSetDensityFunctionTest ${BASICFILTERS_TESTS11}
+--compare-ascii ${TOL}
+	    ${BASELINE_FILES}/bfTvPointSetDensityFunctionOutputAscii.txt
+	    ${TEMP}/bfTvPointSetDensityFunctionOutputAscii.txt
+	 otbPointSetDensityFunctionTest
+	 ${TEMP}/bfTvPointSetDensityFunctionOutputAscii.txt
+)
+
+#------------ otbPointSetToDensityImageFilter ---------------------
+
+ADD_TEST(bfTuPointSetToDensityImageFilterNew ${BASICFILTERS_TESTS11}
+         otbPointSetToDensityImageFilterNew)
+
+
+
+
 #
 #ADD_TEST(bfTvCountImageFunction ${BASICFILTERS_TESTS11}
 #--compare-ascii ${TOL}
@@ -1306,6 +1331,9 @@ ENDIF(USE_FFTWD)
 SET(BasicFilters_SRCS11
 otbExtractROIResample.cxx
 otbPointSetFunctionNew.cxx
+otbPointSetDensityFunctionNew.cxx
+otbPointSetDensityFunctionTest.cxx
+otbPointSetToDensityImageFilterNew.cxx
 #otbCountImageFunctionTest.cxx	
 #otbCountImageFilterNew.cxx
 #otbCountImageFilterTest.cxx	

Testing/Code/BasicFilters/otbBasicFiltersTests11.cxx

@@ -29,6 +29,10 @@ void RegisterTests()
 {
 REGISTER_TEST(otbExtractROIResample);
 REGISTER_TEST(otbPointSetFunctionNew); 
+REGISTER_TEST(otbPointSetToDensityImageFilterNew); 
+REGISTER_TEST(otbPointSetDensityFunctionNew); 
+REGISTER_TEST(otbPointSetDensityFunctionTest);
+
 //REGISTER_TEST(otbCountImageFunctionTest);
 //REGISTER_TEST(otbCountImageFilterNew); 
 //REGISTER_TEST(otbCountImageFilterTest);  

Testing/Code/BasicFilters/otbCountImageFilterTest.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 "otbImage.h"
-#include "otbImageFileReader.h"
-#include "otbImageFileWriter.h"
-
-#include <stdio.h>
-#include "otbCountImageFilter.h"
-#include "otbSiftFastImageFilter.h"
-#include "otbSimplePointCountStrategy.h"
-#include "itkPointSet.h"
-#include "itkVariableLengthVector.h"
-
-#include <iostream>
-
-int otbCountImageFilterTest(int argc, char* argv[] )
-{
-  const char *  infname = argv[1];            
-  const char * outfname = argv[2];
-  const unsigned char scales = atoi(argv[3]);
-  const unsigned char radius = atoi(argv[4]);
-  const   unsigned int                                      Dimension = 2;
-  typedef float                                             PixelType; 
-
-  typedef otb::Image< PixelType, Dimension >                ImageType;
-  typedef ImageType::IndexType                              IndexType;
-  
-  typedef otb::ImageFileReader<ImageType>                   ReaderType;
-  typedef otb::ImageFileWriter<ImageType>                   WriterType;
-  
-  typedef itk::VariableLengthVector<PixelType>              RealVectorType;
-  typedef itk::PointSet<RealVectorType,Dimension>           PointSetType;
-  typedef otb::SiftFastImageFilter<ImageType,PointSetType>  DetectorType;
-  
-  typedef otb::Count<PointSetType,unsigned int ,IndexType>  CounterType;
-  
-  typedef otb::CountImageFilter< ImageType,DetectorType,
-                                   CounterType, ImageType>  FilterType;
-
-
-  /** Instanciation of the reader */
-  ReaderType::Pointer reader = ReaderType::New();
-
-  reader->SetFileName(infname);
-  reader->Update();
-  
-  ImageType::SizeType                     size;
-  size = reader->GetOutput()->GetLargestPossibleRegion().GetSize();
-  
-  /**Instancitation of an object*/
-  FilterType::Pointer    filter =     FilterType::New();
-  
-  /** Instanciation of the detector : */
-  DetectorType::Pointer detector = filter->GetDetector();
-  detector->SetNumberOfScales(scales);
-  
-
-  filter->SetInput(reader->GetOutput()); 
-  filter->SetNeighborhoodRadius(radius);
-  filter->SetDetector(detector);  /*Set the number of scales for the detector**/
-
-  
-  /** Output*/
-  WriterType::Pointer writer = WriterType::New();
-  writer->SetFileName(outfname);
-  writer->SetInput(filter->GetOutput());
-  writer->Update();
-
-  return EXIT_SUCCESS;
-}
-

Testing/Code/BasicFilters/otbPointSetDensityFunctionNew.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 <stdio.h>
+
+#include "otbPointSetDensityFunction.h"
+#include "itkPointSet.h"
+#include "itkVariableLengthVector.h"
+
+
+int otbPointSetDensityFunctionNew(int, char* [] )
+{
+
+  const   unsigned int                                             Dimension = 2;
+  typedef float                                                    PixelType; 
+
+  typedef itk::VariableLengthVector<PixelType>                     RealVectorType;
+  typedef itk::PointSet<RealVectorType,Dimension>                  PointSetType;
+  typedef otb::PointSetDensityFunction <PointSetType,PixelType>    FunctionType;
+  
+  /**Instancitation of an object*/
+  FunctionType    filter();
+  //FunctionType::Pointer    filter = FunctionType::New();
+
+  return EXIT_SUCCESS;
+}
+

Testing/Code/BasicFilters/otbPointSetDensityFunctionTest.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 "otbImageFileReader.h"
+#include "otbPointSetDensityFunction.h"
+#include "itkPointSet.h"
+#include "itkVariableLengthVector.h"
+#include "otbSiftFastImageFilter.h"
+
+#include <iostream>
+
+int otbPointSetDensityFunctionTest(int argc, char* argv[] )
+{
+  
+  const char * outfname = argv[1];
+  
+  const   unsigned int                                             Dimension = 2;
+  typedef float                                                    PixelType; 
+
+  typedef itk::VariableLengthVector<PixelType>                     RealVectorType;
+  typedef itk::PointSet<RealVectorType,Dimension>                  PointSetType;
+  typedef otb::PointSetDensityFunction <PointSetType,PixelType>    FunctionType;
+
+  
+  /**Instancitation of an object*/
+  PointSetType::Pointer pointSet = PointSetType::New();
+  FunctionType::Pointer     filter = FunctionType::New();
+  std::ofstream outfile(outfname);
+  
+  /** Construction of the pointSet */
+  PointSetType::PointIdentifier count = 0;
+  
+  PointSetType::PointType  pDst ,pSrc;
+  pDst[0] = 12.78 ; 
+  pDst[1] = 18.76 ;
+  pointSet->SetPoint(count++,pDst);
+
+  pDst[0] = 15.78 ; 
+  pDst[1] = 23.76 ;
+  pointSet->SetPoint(count++,pDst);
+  
+  pDst[0] = 9.78 ; 
+  pDst[1] = 5.76 ;
+  pointSet->SetPoint(count++,pDst);
+  
+  
+  filter->SetPointSet(pointSet);
+  filter->SetRadius(2);
+  
+  /**Point we serach around*/
+  pDst[0] = 14.9; pDst[1] = 24;
+  outfile << "Density computed for the point : " << pDst << " is "<< filter->Evaluate(pDst) << std::endl;
+
+  pDst[0] = 9; pDst[1] = 5;
+  outfile << "Density computed for the point : " << pDst << " is "<< filter->Evaluate(pDst) << std::endl;
+
+  outfile.close();
+  
+  return EXIT_SUCCESS;
+}
+

Testing/Code/BasicFilters/otbPointSetFunctionNew.cxx

@@ -31,7 +31,7 @@ int otbPointSetFunctionNew(int, char* [] )
 
   typedef itk::VariableLengthVector<PixelType>              RealVectorType;
   typedef itk::PointSet<RealVectorType,Dimension>           PointSetType;
-  typedef otb::PointSetFunction <PointSetType,PixelType>   FunctionType;
+  typedef otb::PointSetFunction <PointSetType,PixelType>    FunctionType;
   
   /**Instancitation of an object*/
   FunctionType   filter();