Skip to content

Commit

Permalink
Make isClockwise() available at the OGRCurve level
Browse files Browse the repository at this point in the history
  • Loading branch information
@rouault
rouault committed Jun 10, 2022
1 parent 36df9b5 commit c7203e8
Show file tree
Hide file tree
Showing 4 changed files with 281 additions and 118 deletions.
6 changes: 4 additions & 2 deletions ogr/ogr_geometry.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1092,7 +1092,8 @@ class CPL_DLL OGRCurve : public OGRGeometry
std::unique_ptr<Private> m_poPrivate;
public:
ConstIterator(const OGRCurve* poSelf, bool bStart);
ConstIterator(ConstIterator&& oOther) noexcept; // declared but not defined. Needed for gcc 5.4 at least
ConstIterator(ConstIterator&& oOther) noexcept;
ConstIterator& operator=(ConstIterator&& oOther);
~ConstIterator();
const OGRPoint& operator*() const;
ConstIterator& operator++();
Expand Down Expand Up @@ -1143,6 +1144,7 @@ class CPL_DLL OGRCurve : public OGRGeometry
virtual OGRPointIterator* getPointIterator() const = 0;
virtual OGRBoolean IsConvex() const;
virtual double get_Area() const = 0;
virtual int isClockwise() const;

/** Down-cast to OGRSimpleCurve*.
* Implies prior checking that wkbFlatten(getGeometryType()) == wkbLineString or wkbCircularString. */
Expand Down Expand Up @@ -1443,6 +1445,7 @@ class CPL_DLL OGRLineString : public OGRSimpleCurve
// Non-standard from OGRGeometry.
virtual OGRwkbGeometryType getGeometryType() const override;
virtual const char *getGeometryName() const override;
virtual int isClockwise() const override;

/** Return pointer of this in upper class */
inline OGRSimpleCurve* toUpperClass() { return this; }
Expand Down Expand Up @@ -1528,7 +1531,6 @@ class CPL_DLL OGRLinearRing : public OGRLineString
// Non standard.
virtual const char *getGeometryName() const override;
virtual OGRLinearRing *clone() const override;
virtual int isClockwise() const;
virtual void reverseWindingOrder();
virtual void closeRings() override;
OGRBoolean isPointInRing( const OGRPoint* pt,
Expand Down
158 changes: 158 additions & 0 deletions ogr/ogrcurve.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -580,6 +580,8 @@ struct OGRCurve::ConstIterator::Private
{
CPL_DISALLOW_COPY_ASSIGN(Private)
Private() = default;
Private(Private&&) = delete;
Private& operator=(Private&&) = default;

OGRPoint m_oPoint{};
std::unique_ptr<OGRPointIterator> m_poIterator{};
Expand All @@ -596,6 +598,17 @@ OGRCurve::ConstIterator::ConstIterator(const OGRCurve* poSelf, bool bStart):
}
}

OGRCurve::ConstIterator::ConstIterator(ConstIterator&& oOther) noexcept:
m_poPrivate(std::move(oOther.m_poPrivate))
{
}

OGRCurve::ConstIterator& OGRCurve::ConstIterator::operator=(ConstIterator&& oOther)
{
m_poPrivate = std::move(oOther.m_poPrivate);
return *this;
}

OGRCurve::ConstIterator::~ConstIterator() = default;

const OGRPoint& OGRCurve::ConstIterator::operator*() const
Expand Down Expand Up @@ -624,3 +637,148 @@ OGRCurve::ConstIterator OGRCurve::end() const
{
return {this, false};
}

/************************************************************************/
/* epsilonEqual() */
/************************************************************************/

constexpr double EPSILON = 1.0E-5;

static inline bool epsilonEqual(double a, double b, double eps)
{
return ::fabs(a - b) < eps;
}

/************************************************************************/
/* isClockwise() */
/************************************************************************/

/**
* \brief Returns TRUE if the ring has clockwise winding (or less than 2 points)
*
* Assumes that the line is closed.
*
* @return TRUE if clockwise otherwise FALSE.
*/

int OGRCurve::isClockwise() const

{
const int nPointCount = getNumPoints();
if( nPointCount < 3 )
return TRUE;

bool bUseFallback = false;

// Find the lowest rightmost vertex.
auto oIter = begin();
const OGRPoint oStartPoint = *oIter;
OGRPoint oPointBefore = oStartPoint;
OGRPoint oPointBeforeSel;
OGRPoint oPointSel = oStartPoint;
OGRPoint oPointNextSel;
bool bNextPointIsNextSel = false;
int v = 0;

for( int i = 1; i < nPointCount - 1; i++ )
{
++oIter;
OGRPoint oPointCur = *oIter;
if( bNextPointIsNextSel )
{
oPointNextSel = oPointCur;
bNextPointIsNextSel = false;
}
if( oPointCur.getY() < oPointSel.getY() ||
( oPointCur.getY() == oPointSel.getY() &&
oPointCur.getX() > oPointSel.getX() ) )
{
v = i;
oPointBeforeSel = oPointBefore;
oPointSel = oPointCur;
bUseFallback = false;
bNextPointIsNextSel = true;
}
else if( oPointCur.getY() == oPointSel.getY() &&
oPointCur.getX() == oPointSel.getX() )
{
// Two vertex with same coordinates are the lowest rightmost
// vertex. Cannot use that point as the pivot (#5342).
bUseFallback = true;
}
oPointBefore = oPointCur;
}
const OGRPoint oPointN_m2 = *oIter;

if( bNextPointIsNextSel )
{
oPointNextSel = oPointN_m2;
}

// Previous.
if( v == 0 )
{
oPointBeforeSel = oPointN_m2;
}

if( epsilonEqual(oPointBeforeSel.getX(), oPointSel.getX(), EPSILON) &&
epsilonEqual(oPointBeforeSel.getY(), oPointSel.getY(), EPSILON) )
{
// Don't try to be too clever by retrying with a next point.
// This can lead to false results as in the case of #3356.
bUseFallback = true;
}

const double dx0 = oPointBeforeSel.getX() - oPointSel.getX();
const double dy0 = oPointBeforeSel.getY() - oPointSel.getY();

// Following.
if( v + 1 >= nPointCount - 1 )
{
oPointNextSel = oStartPoint;
}

if( epsilonEqual(oPointNextSel.getX(), oPointSel.getX(), EPSILON) &&
epsilonEqual(oPointNextSel.getY(), oPointSel.getY(), EPSILON) )
{
// Don't try to be too clever by retrying with a next point.
// This can lead to false results as in the case of #3356.
bUseFallback = true;
}

const double dx1 = oPointNextSel.getX() - oPointSel.getX();
const double dy1 = oPointNextSel.getY() - oPointSel.getY();

const double crossproduct = dx1 * dy0 - dx0 * dy1;

if( !bUseFallback )
{
if( crossproduct > 0 ) // CCW
return FALSE;
else if( crossproduct < 0 ) // CW
return TRUE;
}

// This is a degenerate case: the extent of the polygon is less than EPSILON
// or 2 nearly identical points were found.
// Try with Green Formula as a fallback, but this is not a guarantee
// as we'll probably be affected by numerical instabilities.
oIter = begin();
oPointBefore = oStartPoint;
++oIter;
auto oPointCur = *oIter;
double dfSum = oStartPoint.getX() * (oPointCur.getY() - oStartPoint.getY());

for( int i = 1; i < nPointCount-1; i++ )
{
++oIter;
auto oPointNext = *oIter;
dfSum += oPointCur.getX() * (oPointNext.getY() - oPointBefore.getY());
oPointBefore = oPointCur;
oPointCur = oPointNext;
}

dfSum += oPointCur.getX() * (oStartPoint.getY() - oPointBefore.getY());

return dfSum < 0;
}
116 changes: 0 additions & 116 deletions ogr/ogrlinearring.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -424,122 +424,6 @@ OGRLinearRing *OGRLinearRing::clone() const
return poNewLinearRing;
}

/************************************************************************/
/* epsilonEqual() */
/************************************************************************/

constexpr double EPSILON = 1.0E-5;

static inline bool epsilonEqual(double a, double b, double eps)
{
return ::fabs(a - b) < eps;
}

/************************************************************************/
/* isClockwise() */
/************************************************************************/

/**
* \brief Returns TRUE if the ring has clockwise winding (or less than 2 points)
*
* @return TRUE if clockwise otherwise FALSE.
*/

int OGRLinearRing::isClockwise() const

{
if( nPointCount < 2 )
return TRUE;

bool bUseFallback = false;

// Find the lowest rightmost vertex.
int v = 0; // Used after for.
for( int i = 1; i < nPointCount - 1; i++ )
{
// => v < end.
if( paoPoints[i].y< paoPoints[v].y ||
( paoPoints[i].y== paoPoints[v].y &&
paoPoints[i].x > paoPoints[v].x ) )
{
v = i;
bUseFallback = false;
}
else if( paoPoints[i].y == paoPoints[v].y &&
paoPoints[i].x == paoPoints[v].x )
{
// Two vertex with same coordinates are the lowest rightmost
// vertex. Cannot use that point as the pivot (#5342).
bUseFallback = true;
}
}

// Previous.
int next = v - 1;
if( next < 0 )
{
next = nPointCount - 1 - 1;
}

if( epsilonEqual(paoPoints[next].x, paoPoints[v].x, EPSILON) &&
epsilonEqual(paoPoints[next].y, paoPoints[v].y, EPSILON) )
{
// Don't try to be too clever by retrying with a next point.
// This can lead to false results as in the case of #3356.
bUseFallback = true;
}

const double dx0 = paoPoints[next].x - paoPoints[v].x;
const double dy0 = paoPoints[next].y - paoPoints[v].y;

// Following.
next = v + 1;
if( next >= nPointCount - 1 )
{
next = 0;
}

if( epsilonEqual(paoPoints[next].x, paoPoints[v].x, EPSILON) &&
epsilonEqual(paoPoints[next].y, paoPoints[v].y, EPSILON) )
{
// Don't try to be too clever by retrying with a next point.
// This can lead to false results as in the case of #3356.
bUseFallback = true;
}

const double dx1 = paoPoints[next].x - paoPoints[v].x;
const double dy1 = paoPoints[next].y - paoPoints[v].y;

const double crossproduct = dx1 * dy0 - dx0 * dy1;

if( !bUseFallback )
{
if( crossproduct > 0 ) // CCW
return FALSE;
else if( crossproduct < 0 ) // CW
return TRUE;
}

// This is a degenerate case: the extent of the polygon is less than EPSILON
// or 2 nearly identical points were found.
// Try with Green Formula as a fallback, but this is not a guarantee
// as we'll probably be affected by numerical instabilities.

double dfSum =
paoPoints[0].x * (paoPoints[1].y - paoPoints[nPointCount-1].y);

for( int i = 1; i < nPointCount-1; i++ )
{
dfSum += paoPoints[i].x * (paoPoints[i+1].y - paoPoints[i-1].y);
}

dfSum +=
paoPoints[nPointCount-1].x *
(paoPoints[0].y - paoPoints[nPointCount-2].y);

return dfSum < 0;
}

/************************************************************************/
/* reverseWindingOrder() */
/************************************************************************/
Expand Down
Loading

0 comments on commit c7203e8

Please sign in to comment.