ENH: improve to handle neighborhood

Code/BasicFilters/otbSpectralAngleDataNodeFeatureFunction.h

@@ -23,6 +23,7 @@
 #include "otbBinarySpectralAngleFunctor.h"
 #include "otbDataNode.h"
 #include "otbPolyLineImageConstIterator.h"
+#include "itkLineConstIterator.h"
 #include "itkVariableLengthVector.h"
 
 namespace otb
@@ -31,10 +32,20 @@ namespace otb
   * \brief Compute a spectral angle based feature alongside a
   * datanode.
   *
-  * This function compute a feature alongside a datanode.
+  * This function compute a spectral angle alongside a datanode.
   * The feature is the mean spectral angle regarding a
   * reference pixel alongside the tested datanode.
   *
+  * Furethemore, it compute the spectral in a neighborhood region located between
+  * +/- m_StartNeighborhoodRadius and StopNeighborhoodRadius.
+  *
+  * The output has five elements:
+  * - #0: \$f \frac{mean spectral angle alongside the datanode}{mean spectral angle of the neighborhood} \$f
+  * - #1: cumulate spectral angle alongside the datanode
+  * - #2: number of pixel alongside the datanode
+  * - #3: cumulate spectral angle of the neighborhood
+  * - #4: number of pixel in the neighborhood
+  *
   * \ingroup Functions
   * \sa DataNodeImageFunction
   * \sa NDVIDataNodeFeatureFunction
@@ -65,6 +76,8 @@ public:
   /** Some typedefs. */
   typedef typename Superclass::DataNodeType           DataNodeType;
   typedef typename DataNodeType::LineType             LineType;
+  typedef typename LineType::Pointer                  LinePointer;
+  typedef typename LineType::ContinuousIndexType      ContinuousIndexType;
 
   typedef TImage                                      InputImageType;
   typedef typename InputImageType::ConstPointer       InputImageConstPointer;
@@ -79,11 +92,13 @@ public:
 
   typedef itk::VariableLengthVector<PrecisionType>    ReferencePixelType;
 
-  typedef PolyLineImageConstIterator<InputImageType, LineType>
-                                                      ImageLineIteratorType;
+  typedef PolyLineImageConstIterator<InputImageType, LineType> ImageLineIteratorType;
+  typedef itk::LineConstIterator<InputImageType>      LineIteratorType;
+
   typedef Functor::BinarySpectralAngleFunctor<PixelType, ReferencePixelType, PrecisionType>
                                                       SpectralAngleFunctorType;
 
+  typedef std::pair<IndexType, IndexType>             IndexPairType; 
   typedef std::vector<PrecisionType>                  OutputType;
 
   virtual OutputType Evaluate( const DataNodeType& node ) const;
@@ -92,6 +107,12 @@ public:
   itkGetConstMacro(RefPixel, PixelType);
   itkSetMacro(RefPixel, PixelType);
 
+  itkGetConstMacro(StartNeighborhoodRadius, unsigned int);
+  itkSetMacro(StartNeighborhoodRadius, unsigned int);
+  itkGetConstMacro(StopNeighborhoodRadius, unsigned int);
+  itkSetMacro(StopNeighborhoodRadius, unsigned int);
+
+
 protected:
   SpectralAngleDataNodeFeatureFunction();
   ~SpectralAngleDataNodeFeatureFunction() {}
@@ -104,6 +125,10 @@ private:
   /** SpectralAngle Functor & ReferencePixel*/
   ReferencePixelType          m_RefPixel;
   SpectralAngleFunctorType    m_SpectralAngleFunctor;
+  /** Start neighborhod radius */
+  unsigned int m_StartNeighborhoodRadius;
+ /** Stop neighborhod radius */
+  unsigned int m_StopNeighborhoodRadius;
 };
 
 }

Code/BasicFilters/otbSpectralAngleDataNodeFeatureFunction.txx

@@ -20,6 +20,8 @@
 
 #include "otbSpectralAngleDataNodeFeatureFunction.h"
 
+//#include "itkImageIteratorWithIndex.h"
+
 namespace otb
 {
 
@@ -28,7 +30,7 @@ namespace otb
  */
 template <class TImage, class TCoordRep, class TPrecision>
 SpectralAngleDataNodeFeatureFunction<TImage, TCoordRep, TPrecision>
-::SpectralAngleDataNodeFeatureFunction()
+::SpectralAngleDataNodeFeatureFunction() : m_StartNeighborhoodRadius(2), m_StopNeighborhoodRadius(3)
 {
   //Example for QuickBird images (on a specific image)
   m_RefPixel.SetSize(4);
@@ -60,6 +62,7 @@ typename SpectralAngleDataNodeFeatureFunction<TImage, TCoordRep, TPrecision>
 ::Evaluate( const DataNodeType& node ) const
 {
 
+  // TODO faire avce un ikk
   const typename ImageLineIteratorType::PathType* path;
 
   switch (node.GetNodeType())
@@ -79,60 +82,186 @@ typename SpectralAngleDataNodeFeatureFunction<TImage, TCoordRep, TPrecision>
           path = node.GetPolygonExteriorRing();
           break;
           }
-        default:
+          default:
           {
           itkExceptionMacro(<< "This DataNode type is not handle yet");
           break;
           }
         }
 
-  ImageLineIteratorType lineIt(this->GetInputImage(), path);
-  lineIt.GoToBegin();
+  std::vector< std::pair<IndexType, IndexType> > splitedLineIdNeigh;
+  std::vector< std::pair<IndexType, IndexType> > splitedLineIdCentral;
+  // Split line and polygon into segment (ie. line with two vertex
+  typename LineType::VertexListConstIteratorType it1 = path->GetVertexList()->Begin();
+  typename LineType::VertexListConstIteratorType it2 = path->GetVertexList()->Begin();
+  typename LineType::VertexListConstIteratorType itStop = path->GetVertexList()->End();
+   
+  ++it2;
+  if( it2 == itStop )
+    {
+      itkExceptionMacro(<< "Invalid DataNode, must at least contain two points");
+    }
+  
+ 
 
-  double accSpectralAngle=0.;
-  double nbVisitedPixel=0.;
+  //unsigned int count = 0;
+  while ( it1 != itStop && it2 != itStop)
+    {
+      IndexType id1, id2;
+      id1[0] = static_cast<int>(it1.Value()[0]);
+      id1[1] = static_cast<int>(it1.Value()[1]);
+      id2[0] = static_cast<int>(it2.Value()[0]);
+      id2[1] = static_cast<int>(it2.Value()[1]);
+      splitedLineIdCentral.push_back(IndexPairType( id1, id2 ));
+      
 
-  while(!lineIt.IsAtEnd())
+      for(unsigned int j=m_StartNeighborhoodRadius; j<m_StopNeighborhoodRadius; j++)
+        {
+          IndexType shift11, shift12;
+          shift11[0] = id1[0]-j;
+          shift11[1] = id1[1]-j;
+          shift12[0] = id1[0]+j;
+          shift12[1] = id1[1]+j;
+          IndexType shift21, shift22;
+          shift21[0] = id2[0]-j;
+          shift21[1] = id2[1]-j;
+          shift22[0] = id2[0]+j;
+          shift22[1] = id2[1]+j;
+          
+          splitedLineIdNeigh.push_back(IndexPairType( shift11, shift21 ));
+          splitedLineIdNeigh.push_back(IndexPairType( shift12, shift22 ));
+        }
+      ++it1;
+      ++it2;
+    }
+
+  if( node.GetNodeType() == FEATURE_POLYGON )
     {
-    if(this->IsInsideBuffer(lineIt.GetIndex()))
-      {
-      /*
-      std::cout << "m_SpectralAngleFunctor : "
-                << m_SpectralAngleFunctor(lineIt.Get(), this->GetRefPixel())
-                << "  -  currentPixel : "
-                << lineIt.Get()
-                << "  -  RefPixel : "
-                << this->GetRefPixel()
-                << std::endl;
-      */
-      PixelType currPixel;
-      currPixel.SetSize(std::min(this->GetRefPixel().Size(), lineIt.Get().Size()));
-      for (unsigned int i=0; i<std::min(this->GetRefPixel().Size(), lineIt.Get().Size()); i++)
+      it2--;
+      IndexType id1, id2;
+      id1[0] = static_cast<int>(path->GetVertexList()->Begin().Value()[0]);
+      id1[1] = static_cast<int>(path->GetVertexList()->Begin().Value()[1]);
+      id2[0] = static_cast<int>(it2.Value()[0]);
+      id2[1] = static_cast<int>(it2.Value()[1]);
+
+       splitedLineIdCentral.push_back(IndexPairType( id1, id2 ));
+
+      for(unsigned int j=m_StartNeighborhoodRadius; j<m_StopNeighborhoodRadius; j++)
         {
-        currPixel[i] = (lineIt.Get())[i];
+          IndexType shift11, shift12;
+          shift11[0] = id1[0]-j;
+          shift11[1] = id1[1]-j;
+          shift12[0] = id1[0]+j;
+          shift12[1] = id1[1]+j;
+          IndexType shift21, shift22;
+          shift21[0] = id2[0]-j;
+          shift21[1] = id2[1]-j;
+          shift22[0] = id2[0]+j;
+          shift22[1] = id2[1]+j;
+          
+          splitedLineIdNeigh.push_back(IndexPairType( shift11, shift21 ));
+          splitedLineIdNeigh.push_back(IndexPairType( shift12, shift22 ));
         }
-      accSpectralAngle += m_SpectralAngleFunctor(currPixel, this->GetRefPixel());
-      nbVisitedPixel += 1;
-      }
-    ++lineIt;
     }
 
+
+  double centralAccSpectralAngle=0.;
+  double centralNbVisitedPixel=0.;
+
+   for(unsigned int i=0; i<splitedLineIdCentral.size(); i++)
+     {
+       LineIteratorType lineIt( this->GetInputImage(), splitedLineIdCentral[i].first, splitedLineIdCentral[i].second);
+       lineIt.GoToBegin();
+       //itk::ImageIteratorWithIndex<TImage> lol(this->GetInputImage());
+       while(!lineIt.IsAtEnd())
+         {
+           if(this->IsInsideBuffer(lineIt.GetIndex()))
+             {
+               /*
+                 std::cout << "m_SpectralAngleFunctor : "
+                 << m_SpectralAngleFunctor(lineIt.Get(), this->GetRefPixel())
+                 << "  -  currentPixel : "
+                 << lineIt.Get()
+                 << "  -  RefPixel : "
+                 << this->GetRefPixel()
+                 << std::endl;
+               */
+               PixelType currPixel;
+               currPixel.SetSize(std::min(this->GetRefPixel().Size(), lineIt.Get().Size()));
+               for (unsigned int i=0; i<std::min(this->GetRefPixel().Size(), lineIt.Get().Size()); i++)
+                 {
+                   currPixel[i] = (lineIt.Get())[i];
+                 }
+               centralAccSpectralAngle += m_SpectralAngleFunctor(currPixel, this->GetRefPixel());
+               centralNbVisitedPixel += 1;
+             }
+           ++lineIt;
+         }
+     }
+
+   double neighAccSpectralAngle=0.;
+   double neighNbVisitedPixel=0.;
+
+   for(unsigned int i=0; i<splitedLineIdNeigh.size(); i++)
+     {
+       LineIteratorType lineIt( this->GetInputImage(), splitedLineIdCentral[i].first, splitedLineIdCentral[i].second);
+       lineIt.GoToBegin();
+       
+       while(!lineIt.IsAtEnd())
+         {
+           if(this->IsInsideBuffer(lineIt.GetIndex()))
+             {
+               /*
+                 std::cout << "m_SpectralAngleFunctor : "
+                 << m_SpectralAngleFunctor(lineIt.Get(), this->GetRefPixel())
+                 << "  -  currentPixel : "
+                 << lineIt.Get()
+                 << "  -  RefPixel : "
+                 << this->GetRefPixel()
+                 << std::endl;
+               */
+               PixelType currPixel;
+               currPixel.SetSize(std::min(this->GetRefPixel().Size(), lineIt.Get().Size()));
+               for (unsigned int i=0; i<std::min(this->GetRefPixel().Size(), lineIt.Get().Size()); i++)
+                 {
+                   currPixel[i] = (lineIt.Get())[i];
+                 }
+               neighAccSpectralAngle += m_SpectralAngleFunctor(currPixel, this->GetRefPixel());
+               neighNbVisitedPixel += 1;
+             }
+           ++lineIt;
+         }
+     }
+
   OutputType output;
 
-  if(nbVisitedPixel == 0)
+  double meanCurr = 0.;
+  if( centralNbVisitedPixel != 0.)
+    {
+      meanCurr /= static_cast<double>(centralAccSpectralAngle);
+    }
+
+  double meanNeigh = 0.;
+  if( neighNbVisitedPixel != 0.)
+    {
+      meanNeigh /= static_cast<double>(neighAccSpectralAngle);
+    }
+
+  if(meanNeigh == 0.)
     {
-    //std::cout << "nbValidPixel: " << nbValidPixel << "nbVisitedPixel" << nbVisitedPixel << std::endl;
-    //itkExceptionMacro(<< "The DataNode and the Support Image are disjointed");
-    output.push_back(static_cast<PrecisionType>(0.));
+      output.push_back(static_cast<PrecisionType>(0.));
     }
   else
     {
-    output.push_back(static_cast<PrecisionType>(accSpectralAngle/nbVisitedPixel));
+      output.push_back(static_cast<PrecisionType>(meanCurr/meanNeigh));
     }
 
-  output.push_back(static_cast<PrecisionType>(accSpectralAngle));
-  output.push_back(static_cast<PrecisionType>(nbVisitedPixel));
+  output.push_back(static_cast<PrecisionType>(centralAccSpectralAngle));
+  output.push_back(static_cast<PrecisionType>(centralNbVisitedPixel));
 
+  output.push_back(static_cast<PrecisionType>(neighAccSpectralAngle));
+  output.push_back(static_cast<PrecisionType>(neighNbVisitedPixel));
+  
   return output;
 }
 

Code/Common/otbDataNode.h

@@ -267,7 +267,7 @@ protected:
   /** Constructor */
   DataNode();
   /** Destructor */
-  ~DataNode() {}
+  virtual ~DataNode() {}
   /** PrintSelf method */
   void PrintSelf(std::ostream& os, itk::Indent indent) const;