Commit d1e26c75 authored by Cyrille Valladeau's avatar Cyrille Valladeau
Browse files

Ajout d'une classe pour la selection de ROI polygonale.

+ Passage en classe template de otbPolygone et correction des tests l'utilisant.
parent 889b04d9
......@@ -54,7 +54,7 @@ class ITK_EXPORT PolyLineParametricPathWithValue
/** Derived typedefs */
typedef typename Superclass::VertexType VertexType;
typedef typename Superclass::VertexListType VertexListType;
typedef typename Superclass::ContinuousIndexType ContinuousIndexType;
itkGetMacro(Key,std::string);
......
/*=========================================================================
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];
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
......@@ -33,16 +33,17 @@ namespace otb
*
* \sa otb::PolyLineParametricPathWithValue
*/
template<class TValue=double>
class ITK_EXPORT Polygon
: public PolyLineParametricPathWithValue<double,2>
: public PolyLineParametricPathWithValue<TValue,2>
{
public:
/** Standard typedefs */
typedef Polygon Self;
typedef PolyLineParametricPathWithValue<double,2> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Standard typedefs */
typedef Polygon Self;
typedef PolyLineParametricPathWithValue<TValue,2> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
typedef TValue ValueType;
/** Type macro */
itkNewMacro(Self);
......@@ -50,9 +51,10 @@ class ITK_EXPORT Polygon
itkTypeMacro(Polygon, PolyLineParametricPathWithValue);
/** Derived typedefs */
typedef Superclass::VertexType VertexType;
typedef Superclass::VertexListType VertexListType;
typedef VertexListType::ConstIterator VertexListIteratorType;
typedef typename Superclass::VertexType VertexType;
typedef typename Superclass::VertexListType VertexListType;
typedef typename Superclass::ContinuousIndexType ContinuousIndexType;
typedef typename VertexListType::ConstIterator VertexListIteratorType;
itkSetMacro(Epsilon,double);
itkGetMacro(Epsilon,double);
......@@ -125,4 +127,8 @@ private:
double m_Epsilon;
};
}// End namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbPolygon.txx"
#endif
#endif
/*=========================================================================
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.
*/
template<class TValue>
bool
Polygon<TValue>
::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];
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.
*/
template<class TValue>
bool
Polygon<TValue>
::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))