New class used in vector RCC8

Code/Common/otbPolygon.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.
+
+=========================================================================*/
+#ifndef _otbPolygon_txx
+#define _otbPolygon_txx
+
+#include "otbPolygon.h"
+
+namespace otb
+{
+  /**
+   * Check wether point is strictly inside the polygon.
+   * \param point The point to check.
+   * \return True if the point is inside the polygon.
+   */
+  bool 
+  Polygon
+  ::IsInside(VertexType point)
+  {
+    double x = point[0];
+    double y = point[1];
+    unsigned int crossingCount = 0;
+    VertexListIteratorType it = this->GetVertexList()->Begin();
+    double xa = it.Value()[0];
+    double ya = it.Value()[1];
+    ++it; 
+    while(it != this->GetVertexList()->End())
+      {      
+	double xb = it.Value()[0];
+	double yb = it.Value()[1];
+	if(vcl_abs(xb-xa)<m_Epsilon)
+	{
+	  if(ya>yb && xa>x && y>=yb && y<ya)
+	    {
+	      ++crossingCount;
+	    }
+	  else if (ya<yb && xa>x && y>=ya && y<yb) 
+	    {
+	      ++crossingCount;
+	    }
+	}
+	else if(vcl_abs(yb-ya)>=m_Epsilon)
+	  {
+	    double xcross = xa + (xb - xa)*(y - ya)/(yb-ya);
+	    
+	    if(ya>yb && xcross>x && y >= yb && y < ya)
+	      {
+		++crossingCount;
+	      }
+	    else if(ya<yb && xcross>x && y >= ya && y < yb)
+	      {
+		++crossingCount;
+	      }
+	  }
+	xa = xb;
+	ya = yb;
+	++it;
+      }
+    double xb = this->GetVertexList()->Begin().Value()[0];
+    double yb = this->GetVertexList()->Begin().Value()[1];
+    double xmin = std::min(xa,xb);
+    double xmax = std::max(xa,xb);
+    if(vcl_abs(xb-xa)< m_Epsilon)
+      {
+	if(ya>yb && xa>x && y>=yb && y<ya)
+	    {
+	      ++crossingCount;
+	    }
+	  else if (ya<yb && xa>x && y>=ya && y<yb)
+	    {
+	      ++crossingCount;
+	    }
+      }
+    else if(vcl_abs(yb-ya)>=m_Epsilon)
+      {
+	double xcross = xa + (xb - xa)*(y - ya)/(yb-ya);
+	
+	if(ya>yb && xcross>x && y >= yb && y < ya)
+	  {
+	    ++crossingCount;
+	  }
+	else if(ya<yb && xcross>x && y >= ya && y < yb)
+	  {
+	    ++crossingCount;
+	  }
+      }
+    //std::cout<<"Crossing count: "<<crossingCount<<std::endl;
+    return (crossingCount%2 == 1);
+  }
+  /**
+   * Check wether point is strictly on the edge of the polygon.
+   * \param point The point to check.
+   * \return True if the point is on the edge of the polygon.
+   */
+  bool 
+  Polygon
+  ::IsOnEdge(VertexType point)
+  {
+    //std::cout<<"Checking point "<<point<<std::endl;
+    bool resp = false;
+    double x = point[0];
+    double y = point[1];
+    double xb,yb; 
+    VertexListIteratorType it = this->GetVertexList()->Begin();
+    double xa = it.Value()[0];
+    double ya = it.Value()[1];
+    double xbegin = xa;
+    double ybegin = ya;
+    ++it;
+    while(!resp && it != this->GetVertexList()->End())
+      {
+	xb = it.Value()[0];
+	yb = it.Value()[1];
+	if(vcl_abs(xb-xa)>=m_Epsilon)
+	  {
+	    double cd = (yb - ya)/(xb - xa);
+	    double oo = (ya - cd * xa);
+	    double xmin = std::min(xa,xb);
+	    double xmax = std::max(xa,xb);
+	    if ((vcl_abs(y - cd * x - oo)<m_Epsilon)
+		&& ( x <= xmax)
+		&& ( x >= xmin))
+	      {
+		//std::cout<<"Case 1: Point is on segment "<<xa<<", "<<ya<<" <-> "<<xb<<", "<<yb<<std::endl;
+		resp = true;
+	      }
+	  }
+	else
+	  {
+	    double ymin = std::min(ya,yb);
+	    double ymax = std::max(ya,yb);
+	    if((vcl_abs(x-xa)<m_Epsilon)
+	       &&(y<=ymax)
+	       &&(y>=ymin))
+	      {
+		resp = true;
+		//std::cout<<"Case 2: Point is on segment "<<xa<<", "<<ya<<" <-> "<<xb<<", "<<yb<<std::endl;
+	    }
+	  }
+	xa = xb;
+	ya = yb;
+	++it;
+      }
+    xb = xbegin;
+    yb = ybegin;
+    if(vcl_abs(xb-xa)>=m_Epsilon)
+      {
+	double cd = (yb - ya)/(xb - xa);
+	double oo = (ya - cd * xa);
+	double xmin = std::min(xa,xb);
+	double xmax = std::max(xa,xb);
+	
+	if ((vcl_abs(y - cd * x - oo)<m_Epsilon)
+	    && ( x <= xmax)
+	    && ( x >= xmin))
+	  {
+	    resp = true;
+	    //std::cout<<"Case 1: Point is on segment "<<xa<<", "<<ya<<" <-> "<<xb<<", "<<yb<<std::endl;	
+	  }
+      }
+    else 
+      {
+	double ymin = std::min(ya,yb);
+	double ymax = std::max(ya,yb);	
+	if((vcl_abs(x-xa)<=m_Epsilon)
+	   &&(y<=ymax)
+	   &&(y>=ymin))
+	  {
+	    resp = true;
+	    //std::cout<<"Case 2: Point is on segment "<<xa<<", "<<ya<<" <-> "<<xb<<", "<<yb<<std::endl;
+	  }
+      }
+    //std::cout<<"Result: "<<resp<<std::endl;
+    return resp;
+  }
+  /**
+   * Returns the number of crossings of the polygon with a given segment.
+   * \param a First point of the segment,
+   * \param b Second point of the segment,
+   * \return the number of strict crossings of segment [ab] with the polygon.
+   */
+  unsigned int 
+  Polygon
+  ::NbCrossing(VertexType a, VertexType b)
+  {
+    unsigned int resp = 0;
+    VertexListIteratorType it = this->GetVertexList()->Begin();
+    VertexListIteratorType it_end = this->GetVertexList()->End();
+    VertexType current = it.Value();
+    VertexType first = current;
+    ++it;
+    while ( it != it_end)
+      {
+	//std::cout<<"Testing Is crossing "<<a<<" "<<b<<current<<it.Value()<<" = "<<IsCrossing(a,b,current,it.Value())<<std::endl;
+	if(IsCrossing(a,b,current,it.Value()))
+	  {
+	    ++resp;
+	  }
+	current = it.Value();
+	++it;
+      }
+    //std::cout<<"Testing Is crossing "<<a<<" "<<b<<current<<first<<" = "<<IsCrossing(a,b,current,first)<<std::endl;
+    if(IsCrossing(a,b,current,first))
+      {
+	++resp;
+      }
+    return resp;
+  }
+  /**
+   * Returns the number of touchings without crossing of the polygon with a given segment.
+   * \param a First point of the segment,
+   * \param b Second point of the segment,
+   * \return the number of touchings without crossing of segment [ab] with the polygon.
+   */
+  unsigned int 
+  Polygon
+  ::NbTouching(VertexType a, VertexType b)
+  {
+    unsigned int resp = 0;
+    VertexListIteratorType it = this->GetVertexList()->Begin();
+    VertexListIteratorType it_end = this->GetVertexList()->End();
+    VertexType current = it.Value();
+    VertexType first = current;
+    ++it;
+    while ( it != it_end)
+      {
+	//std::cout<<"IsTouching "<<a<<" "<<b<<", "<<current<<" "<<it.Value()<<" -> "<<IsTouching(a,b,current,it.Value())<<std::endl;
+	if (IsTouching(a,b,current,it.Value()))
+	  {
+	    ++resp;
+	  }
+	current = it.Value();
+	++it;
+      }
+    //std::cout<<"IsTouching "<<a<<" "<<b<<", "<<current<<" "<<first<<" -> "<<IsTouching(a,b,current,first)<<std::endl;
+    if(IsTouching(a,b,current,first))
+      {
+	++resp;
+      }
+    return resp;
+  }
+  /**
+   * Check wether two segments [a1a2] and [b1b2] are strictly crossing.
+   * \param a1 First point of the first segment,
+   * \param a1 Second point of the first segment,
+   * \param a1 First point of the second segment,
+   * \param a1 Second point of the second segment.
+   * \return True if the two segments are strictly crossing.
+   */
+  bool 
+  Polygon
+  ::IsCrossing(VertexType a1, VertexType a2, VertexType b1, VertexType b2)
+  {
+    bool resp = false;
+    double xbmin = std::min(b1[0],b2[0]);
+    double xbmax = std::max(b1[0],b2[0]);
+    double ybmin = std::min(b1[1],b2[1]);
+    double ybmax = std::max(b1[1],b2[1]);
+    double xamin = std::min(a1[0],a2[0]);
+    double xamax = std::max(a1[0],a2[0]);
+    double yamin = std::min(a1[1],a2[1]);
+    double yamax = std::max(a1[1],a2[1]);    
+  if(vcl_abs(a1[0]-a2[0])<m_Epsilon && vcl_abs(b1[0]-b2[0])<m_Epsilon)
+    {
+      resp = false;
+    }
+  else if(vcl_abs(a1[0]-a2[0])<m_Epsilon)
+    {
+      double cd2 = (b2[1] - b1[1]) / (b2[0] - b1[0]);
+      double oo2 = b1[1] - cd2 * b1[0];      
+      double ycross = cd2 * a1[0] + oo2;
+      resp = ( xbmin < a1[0] && xbmax > a1[0]
+	   && yamin < ycross && yamax > ycross);
+    }
+  else if(vcl_abs(b1[0]-b2[0])<m_Epsilon)
+    {
+      double cd1 = (a2[1] - a1[1]) / (a2[0] - a1[0]);
+      double oo1 = a1[1] - cd1 * a1[0];
+      double ycross = cd1 * b1[0] + oo1;
+      resp = ( xamin < b1[0] && xamax > b1[0]
+	   && ybmin < ycross && ybmax > ycross);
+    }
+  else
+    {
+      double cd1 = (a2[1] - a1[1]) / (a2[0] - a1[0]);
+      double oo1 = a1[1] - cd1 * a1[0];
+      double cd2 = (b2[1] - b1[1]) / (b2[0] - b1[0]);
+      double oo2 = b1[1] - cd2 * b1[0];
+      if(cd1 != cd2)
+	{
+	  double xcross = (oo2 - oo1) / (cd1 - cd2);	  
+	  resp = (xamin < xcross  && xbmin < xcross
+	     &&   xamax > xcross  && xbmax > xcross);
+	}
+    }
+  return resp;
+  }
+  /**
+   * Check wether two segments[a1a2] and [b1b2] are touching without crossing.
+   * \param a1 First point of the first segment,
+   * \param a1 Second point of the first segment,
+   * \param a1 First point of the second segment,
+   * \param a1 Second point of the second segment.
+   * \return True if the two segments are touching without crossing.
+   */
+  bool 
+  Polygon
+  ::IsTouching(VertexType a1, VertexType a2, VertexType b1, VertexType  b2)
+  {
+    bool resp = false;
+    double xbmin = std::min(b1[0],b2[0]);
+    double xbmax = std::max(b1[0],b2[0]);
+    double ybmin = std::min(b1[1],b2[1]);
+    double ybmax = std::max(b1[1],b2[1]);
+    double xamin = std::min(a1[0],a2[0]);
+    double xamax = std::max(a1[0],a2[0]); 
+    double yamin = std::min(a1[1],a2[1]);
+    double yamax = std::max(a1[1],a2[1]);
+    if(vcl_abs(a1[0]-a2[0])<m_Epsilon && vcl_abs(b1[0]-b2[0])<m_Epsilon)
+      {
+	resp = (vcl_abs(a1[0]-b1[0])<m_Epsilon)
+	      && ((a1[1]<=ybmax && a1[1]>=ybmin)
+	      ||  (a2[1]<=ybmax && a2[1]>=ybmin)
+	      ||  (b1[1]<=yamax && b1[1]>=yamin)
+	      ||  (b2[1]<=yamax && b2[1]>=yamin) );
+      }
+    else if(vcl_abs(a1[0]-a2[0])<m_Epsilon)
+      {
+	double cd2 = (b2[1] - b1[1]) / (b2[0] - b1[0]);
+	double oo2 = b1[1] - cd2 * b1[0];
+	
+	if(vcl_abs(a1[1]-cd2*a1[0]-oo2)<m_Epsilon)
+	  {
+	    resp = ( a1[0]>=xbmin && a1[0]<=xbmax);
+	  }
+	else if(vcl_abs(a2[1]-cd2*a2[0]-oo2)<m_Epsilon)
+	  {
+	    resp = ( a2[0]>=xbmin && a2[0]<=xbmax);
+	  } 
+      }
+    else if(vcl_abs(b1[0]-b2[0])<m_Epsilon)
+      {
+	double cd1 = (a2[1] - a1[1]) / (a2[0] - a1[0]);
+	double oo1 = a1[1] - cd1 * a1[0];
+	
+	if(vcl_abs(b1[1]-cd1*b1[0]-oo1)<m_Epsilon)
+	  {
+	    resp = ( b1[0]>=xamin && b1[0]<=xamax);
+	  }
+	else if(vcl_abs(b2[1]-cd1*b2[0]-oo1)<m_Epsilon)
+	  {
+	    resp = ( b2[0]>=xamin && b2[0]<=xamax);
+	  } 
+      }
+    else
+      {
+	double cd1 = (a2[1] - a1[1]) / (a2[0] - a1[0]);
+	double oo1 = a1[1] - cd1 * a1[0];
+	double cd2 = (b2[1] - b1[1]) / (b2[0] - b1[0]);
+	double oo2 = b1[1] - cd2 * b1[0];
+	if(vcl_abs(cd1-cd2)<m_Epsilon && vcl_abs(oo1-oo2)<m_Epsilon)
+	  {
+	    resp =((xamin <= xbmax && xamin >= xbmin)
+	      ||   (xamax <= xbmax && xamax >= xbmin)
+	      ||   (xbmin <= xamax && xbmin >= xamin)
+	      ||   (xbmax <= xamax && xbmax >= xamin) );
+	  }
+	else
+	  {
+	    if(vcl_abs(a1[1]-cd2*a1[0]-oo2)<m_Epsilon)
+	      {
+		resp = ( a1[0]>=xbmin && a1[0]<=xbmax);
+	      }
+	    else if(vcl_abs(a2[1]-cd2*a2[0]-oo2)<m_Epsilon)
+	      {
+		resp = ( a2[0]>=xbmin && a2[0]<=xbmax);
+	      } 
+	    if(vcl_abs(b1[1]-cd1*b1[0]-oo1)<m_Epsilon)
+	      {
+		resp = ( b1[0]>=xamin && b1[0]<=xamax);
+	      }
+	    else if(vcl_abs(b2[1]-cd1*b2[0]-oo1)<m_Epsilon)
+	      {
+		resp = ( b2[0]>=xamin && b2[0]<=xamax);
+	      } 
+	  }
+      }
+    return resp;
+  }
+/**
+ * PrintSelf Method
+ */
+void
+Polygon
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+  Superclass::PrintSelf(os, indent);
+}
+} // End namespace otb
+#endif

Code/Common/otbPolygon.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 _otbPolygon_h
+#define _otbPolygon_h
+
+#include "otbPolyLineParametricPathWithValue.h"
+
+namespace otb
+{
+/** \class Polygon
+ *  \brief This class represent a 2D polygon.
+ *
+ *  It derives from otb::PolyLineParametricPathWithValue. The polygon is considered to be the closed
+ *  path represented by the PolyLineParametricPathWithValue.
+ * 
+ *  It implements some useful methods in the spatial reasoning context, such as IsInside, IsOnEdge, IsTouching,
+ *  IsCrossing ...
+ *
+ * \sa otb::PolyLineParametricPathWithValue
+ */
+class ITK_EXPORT Polygon
+  : public PolyLineParametricPathWithValue<double,2>
+{
+ public:
+  /** Standard typedefs */
+  typedef Polygon                            Self;
+  typedef PolyLineParametricPathWithValue<double,2> Superclass;
+  typedef itk::SmartPointer<Self>            Pointer;
+  typedef itk::SmartPointer<const Self>      ConstPointer;
+  
+  /** Type macro */
+  itkNewMacro(Self);
+  
+  /** Creation through object factory macro */
+  itkTypeMacro(Polygon, PolyLineParametricPathWithValue);
+  
+  /** Derived typedefs */
+  typedef Superclass::VertexType VertexType;
+  typedef Superclass::VertexListType VertexListType;
+  typedef VertexListType::ConstIterator VertexListIteratorType;
+
+  itkSetMacro(Epsilon,double);
+  itkGetMacro(Epsilon,double);
+
+  /**
+   * Check wether point is strictly inside the polygon.
+   * \param point The point to check.
+   * \return True if the point is inside the polygon.
+   */
+  bool IsInside(VertexType point);
+  
+  /**
+   * Check wether point is strictly on the edge of the polygon.
+   * \param point The point to check.
+   * \return True if the point is on the edge of the polygon.
+   */
+  bool IsOnEdge(VertexType point);
+  
+  /**
+   * Returns the number of crossings of the polygon with a given segment.
+   * \param a First point of the segment,
+   * \param b Second point of the segment,
+   * \return the number of strict crossings of segment [ab] with the polygon.
+   */
+  unsigned int NbCrossing(VertexType a, VertexType b);
+
+  /**
+   * Returns the number of touchings without crossing of the polygon with a given segment.
+   * \param a First point of the segment,
+   * \param b Second point of the segment,
+   * \return the number of touchings without crossing of segment [ab] with the polygon.
+   */
+  unsigned int NbTouching(VertexType a, VertexType b);
+
+  /**
+   * Check wether two segments [a1a2] and [b1b2] are strictly crossing.
+   * \param a1 First point of the first segment,
+   * \param a1 Second point of the first segment,
+   * \param a1 First point of the second segment,
+   * \param a1 Second point of the second segment.
+   * \return True if the two segments are strictly crossing.
+   */
+  bool IsCrossing(VertexType a1, VertexType a2, VertexType b1, VertexType b2);
+  
+  /**
+   * Check wether two segments[a1a2] and [b1b2] are touching without crossing.
+   * \param a1 First point of the first segment,
+   * \param a1 Second point of the first segment,
+   * \param a1 First point of the second segment,
+   * \param a1 Second point of the second segment.
+   * \return True if the two segments are touching without crossing.
+   */
+   bool IsTouching(VertexType a1, VertexType a2, VertexType b1, VertexType b2);
+
+protected:
+  /** Constructor */
+  Polygon() 
+    {
+      m_Epsilon = 0.000001;
+    };
+  /** Destructor */
+  virtual ~Polygon() {};
+ /**PrintSelf method */
+  virtual void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+private:
+  Polygon(const Self&); //purposely not implemented
+  void operator=(const Self&); //purposely not implemented
+
+  double m_Epsilon;
+};
+}// End namespace otb
+#endif