Fix left of test for linestrings

nyalldawson · nyalldawson · commit e34d7fb63605 · 2017-11-29T06:45:33.000+11:00
The test was returning invalid results for certain geometries
diff --git a/python/core/geometry/qgsgeometryutils.sip b/python/core/geometry/qgsgeometryutils.sip
@@ -125,7 +125,11 @@ class QgsGeometryUtils
 
     static double leftOfLine( double x, double y, double x1, double y1, double x2, double y2 );
 %Docstring
-Returns < 0 if point(x/y) is left of the line x1,y1 -> x2,y2
+ Returns a value < 0 if the point (``x``, ``y``) is left of the line from (``x1``, ``y1``) -> ( ``x2``, ``y2``).
+ A positive return value indicates the point is to the right of the line.
+
+ If the return value is 0, then the test was unsuccessful (e.g. due to testing a point exactly
+ on the line, or exactly in line with the segment) and the result is undefined.
  :rtype: float
 %End
 
diff --git a/src/core/geometry/qgsgeometryutils.h b/src/core/geometry/qgsgeometryutils.h
@@ -147,7 +147,13 @@ class CORE_EXPORT QgsGeometryUtils
      */
     static QVector<SelfIntersection> getSelfIntersections( const QgsAbstractGeometry *geom, int part, int ring, double tolerance ) SIP_SKIP;
 
-    //! Returns < 0 if point(x/y) is left of the line x1,y1 -> x2,y2
+    /**
+     * Returns a value < 0 if the point (\a x, \a y) is left of the line from (\a x1, \a y1) -> ( \a x2, \a y2).
+     * A positive return value indicates the point is to the right of the line.
+     *
+     * If the return value is 0, then the test was unsuccessful (e.g. due to testing a point exactly
+     * on the line, or exactly in line with the segment) and the result is undefined.
+     */
     static double leftOfLine( double x, double y, double x1, double y1, double x2, double y2 );
 
     /**
diff --git a/src/core/geometry/qgslinestring.cpp b/src/core/geometry/qgslinestring.cpp
@@ -903,9 +903,16 @@ void QgsLineString::addVertex( const QgsPoint &pt )
 double QgsLineString::closestSegment( const QgsPoint &pt, QgsPoint &segmentPt,  QgsVertexId &vertexAfter, bool *leftOf, double epsilon ) const
 {
   double sqrDist = std::numeric_limits<double>::max();
+  double leftOfDist = std::numeric_limits<double>::max();
+  bool prevLeftOf = false;
+  double prevLeftOfX;
+  double prevLeftOfY;
   double testDist = 0;
   double segmentPtX, segmentPtY;
 
+  if ( leftOf )
+    *leftOf = false;
+
   int size = mX.size();
   if ( size == 0 || size == 1 )
   {
@@ -924,14 +931,36 @@ double QgsLineString::closestSegment( const QgsPoint &pt, QgsPoint &segmentPt,
       sqrDist = testDist;
       segmentPt.setX( segmentPtX );
       segmentPt.setY( segmentPtY );
-      if ( leftOf )
-      {
-        *leftOf = ( QgsGeometryUtils::leftOfLine( pt.x(), pt.y(), prevX, prevY, currentX, currentY ) < 0 );
-      }
       vertexAfter.part = 0;
       vertexAfter.ring = 0;
       vertexAfter.vertex = i;
     }
+    if ( leftOf && qgsDoubleNear( testDist, sqrDist ) )
+    {
+      double left = QgsGeometryUtils::leftOfLine( pt.x(), pt.y(), prevX, prevY, currentX, currentY );
+      // if left equals 0, the test could not be performed (e.g. point in line with segment or on segment)
+      // so don't set leftOf in this case, and hope that there's another segment that's the same distance
+      // where we can perform the check
+      if ( !qgsDoubleNear( left, 0 ) )
+      {
+        if ( qgsDoubleNear( testDist, leftOfDist ) && ( left < 0 ) != prevLeftOf )
+        {
+          // we have two possible segments each with equal distance to point, but they disagree
+          // on whether or not the point is to the left of them.
+          // so we test the segments themselves and flip the result.
+          // see https://stackoverflow.com/questions/10583212/elegant-left-of-test-for-polyline
+          *leftOf = QgsGeometryUtils::leftOfLine( currentX, currentY, prevLeftOfX, prevLeftOfY, prevX, prevY ) > 0;
+        }
+        else
+        {
+          *leftOf = left < 0;
+        }
+        prevLeftOf = *leftOf;
+        leftOfDist = testDist;
+        prevLeftOfX = prevX;
+        prevLeftOfY = prevY;
+      }
+    }
   }
   return sqrDist;
 }
diff --git a/tests/src/core/testqgsgeometry.cpp b/tests/src/core/testqgsgeometry.cpp
@@ -3880,6 +3880,60 @@ void TestQgsGeometry::lineString()
   QCOMPARE( v, QgsVertexId( 0, 0, 2 ) );
   QCOMPARE( leftOf, false );
 
+  l35.setPoints( QgsPointSequence() << QgsPoint( 5, 5 )
+                 << QgsPoint( 6, 4 )
+                 << QgsPoint( 4, 4 )
+                 << QgsPoint( 5, 5 ) );
+  QGSCOMPARENEAR( l35.closestSegment( QgsPoint( 2.35, 4 ), p, v, &leftOf ), 2.7225, 4 * DBL_EPSILON );
+  QCOMPARE( p, QgsPoint( 4, 4 ) );
+  QCOMPARE( v, QgsVertexId( 0, 0, 2 ) );
+  QCOMPARE( leftOf, true );
+
+  l35.setPoints( QgsPointSequence() << QgsPoint( 5, 5 )
+                 << QgsPoint( 4, 4 )
+                 << QgsPoint( 6, 4 )
+                 << QgsPoint( 5, 5 ) );
+  QGSCOMPARENEAR( l35.closestSegment( QgsPoint( 2.35, 4 ), p, v, &leftOf ), 2.7225, 4 * DBL_EPSILON );
+  QCOMPARE( p, QgsPoint( 4, 4 ) );
+  QCOMPARE( v, QgsVertexId( 0, 0, 1 ) );
+  QCOMPARE( leftOf, false );
+
+  l35.setPoints( QgsPointSequence() << QgsPoint( 5, 5 )
+                 << QgsPoint( 6, 4 )
+                 << QgsPoint( 4, 4 )
+                 << QgsPoint( 5, 5 ) );
+  QGSCOMPARENEAR( l35.closestSegment( QgsPoint( 3.5, 2 ), p, v, &leftOf ), 4.250000, 4 * DBL_EPSILON );
+  QCOMPARE( p, QgsPoint( 4, 4 ) );
+  QCOMPARE( v, QgsVertexId( 0, 0, 2 ) );
+  QCOMPARE( leftOf, true );
+
+  l35.setPoints( QgsPointSequence() << QgsPoint( 5, 5 )
+                 << QgsPoint( 4, 4 )
+                 << QgsPoint( 6, 4 )
+                 << QgsPoint( 5, 5 ) );
+  QGSCOMPARENEAR( l35.closestSegment( QgsPoint( 3.5, 2 ), p, v, &leftOf ), 4.250000, 4 * DBL_EPSILON );
+  QCOMPARE( p, QgsPoint( 4, 4 ) );
+  QCOMPARE( v, QgsVertexId( 0, 0, 1 ) );
+  QCOMPARE( leftOf, false );
+
+  l35.setPoints( QgsPointSequence() << QgsPoint( 1, 1 )
+                 << QgsPoint( 1, 4 )
+                 << QgsPoint( 2, 2 )
+                 << QgsPoint( 1, 1 ) );
+  QGSCOMPARENEAR( l35.closestSegment( QgsPoint( 1, 0 ), p, v, &leftOf ), 1, 4 * DBL_EPSILON );
+  QCOMPARE( p, QgsPoint( 1, 1 ) );
+  QCOMPARE( v, QgsVertexId( 0, 0, 1 ) );
+  QCOMPARE( leftOf, true );
+
+  l35.setPoints( QgsPointSequence() << QgsPoint( 1, 1 )
+                 << QgsPoint( 2, 2 )
+                 << QgsPoint( 1, 4 )
+                 << QgsPoint( 1, 1 ) );
+  QGSCOMPARENEAR( l35.closestSegment( QgsPoint( 1, 0 ), p, v, &leftOf ), 1, 4 * DBL_EPSILON );
+  QCOMPARE( p, QgsPoint( 1, 1 ) );
+  QCOMPARE( v, QgsVertexId( 0, 0, 1 ) );
+  QCOMPARE( leftOf, false );
+
   //sumUpArea
   QgsLineString l36;
   double area = 1.0; //sumUpArea adds to area, so start with non-zero value
