Merge pull request #4771 from m-kuhn/closest_point

Closest point

Author: m-kuhn

python/core/geometry/qgsgeometryutils.sip

@@ -37,10 +37,12 @@ class QgsGeometryUtils
  :rtype: QgsPoint
 %End
 
-    static double closestPointMeasure( const QgsAbstractGeometry &geom, const QgsPoint &pt );
+    static QgsPoint closestPoint( const QgsAbstractGeometry &geometry, const QgsPoint &point );
 %Docstring
- Returns measure of nearest point on a geometry for a specified point or NaN if geometry does not have measures
- :rtype: float
+ Returns the nearest point on a segment of a ``geometry``
+ for the specified ``point``. The z and m values will be linearly interpolated between
+ the two neighbouring vertices.
+ :rtype: QgsPoint
 %End
 
     static double distanceToVertex( const QgsAbstractGeometry &geom, QgsVertexId id );

python/core/geometry/qgslinestring.sip

@@ -29,6 +29,8 @@ class QgsLineString: QgsCurve
     QgsLineString( const QVector<QgsPoint> &points );
 %Docstring
  Construct a linestring from a vector of points.
+ Z and M type will be set based on the type of the first point
+ in the vector.
 .. versionadded:: 3.0
 %End
 
@@ -73,7 +75,7 @@ class QgsLineString: QgsCurve
 %Docstring
  Returns the z-coordinate of the specified node in the line string.
  \param index index of node, where the first node in the line is 0
- :return: z-coordinate of node, or 0.0 if index is out of bounds or the line
+ :return: z-coordinate of node, or ``nan`` if index is out of bounds or the line
  does not have a z dimension
 .. seealso:: setZAt()
  :rtype: float
@@ -83,7 +85,7 @@ class QgsLineString: QgsCurve
 %Docstring
  Returns the m value of the specified node in the line string.
  \param index index of node, where the first node in the line is 0
- :return: m value of node, or 0.0 if index is out of bounds or the line
+ :return: m value of node, or ``nan`` if index is out of bounds or the line
  does not have m values
 .. seealso:: setMAt()
  :rtype: float

src/core/geometry/qgscurve.cpp

@@ -31,9 +31,12 @@ bool QgsCurve::isClosed() const
   //don't consider M-coordinates when testing closedness
   QgsPoint start = startPoint();
   QgsPoint end = endPoint();
-  return ( qgsDoubleNear( start.x(), end.x(), 1E-8 ) &&
-           qgsDoubleNear( start.y(), end.y(), 1E-8 ) &&
-           qgsDoubleNear( start.z(), end.z(), 1E-8 ) );
+
+  bool closed = qgsDoubleNear( start.x(), end.x(), 1E-8 ) &&
+                qgsDoubleNear( start.y(), end.y(), 1E-8 );
+  if ( is3D() && closed )
+    closed &= qgsDoubleNear( start.z(), end.z(), 1E-8 ) || ( qIsNaN( start.z() ) && qIsNaN( end.z() ) );
+  return closed;
 }
 
 bool QgsCurve::isRing() const

src/core/geometry/qgsgeometryutils.cpp

@@ -93,33 +93,38 @@ QgsPoint QgsGeometryUtils::closestVertex( const QgsAbstractGeometry &geom, const
   return minDistPoint;
 }
 
-double QgsGeometryUtils::closestPointMeasure( const QgsAbstractGeometry &geom, const QgsPoint &pt )
+QgsPoint QgsGeometryUtils::closestPoint( const QgsAbstractGeometry &geometry, const QgsPoint &point )
 {
-  if ( QgsWkbTypes::hasM( geom.wkbType() ) )
-  {
-    QgsPoint closestPoint;
-    QgsVertexId vertexAfter;
-    bool leftOf;
-    geom.closestSegment( pt, closestPoint, vertexAfter, &leftOf, DEFAULT_SEGMENT_EPSILON );
-    if ( vertexAfter.isValid() )
+  QgsPoint closestPoint;
+  QgsVertexId vertexAfter;
+  bool leftOf;
+  geometry.closestSegment( point, closestPoint, vertexAfter, &leftOf, DEFAULT_SEGMENT_EPSILON );
+  if ( vertexAfter.isValid() )
+  {
+    QgsPoint pointAfter = geometry.vertexAt( vertexAfter );
+    if ( vertexAfter.vertex > 0 )
     {
-      QgsPoint pointAfter = geom.vertexAt( vertexAfter );
-      if ( vertexAfter.vertex > 0 )
-      {
-        QgsVertexId vertexBefore = vertexAfter;
-        vertexBefore.vertex--;
-        QgsPoint pointBefore = geom.vertexAt( vertexBefore );
-        double length = pointBefore.distance( pointAfter );
-        double distance = pointBefore.distance( closestPoint );
-        return pointBefore.m() + ( pointAfter.m() - pointBefore.m() ) * distance / length;
-      }
+      QgsVertexId vertexBefore = vertexAfter;
+      vertexBefore.vertex--;
+      QgsPoint pointBefore = geometry.vertexAt( vertexBefore );
+      double length = pointBefore.distance( pointAfter );
+      double distance = pointBefore.distance( closestPoint );
+
+      if ( qgsDoubleNear( distance, 0.0 ) )
+        closestPoint = pointBefore;
+      else if ( qgsDoubleNear( distance, length ) )
+        closestPoint = pointAfter;
       else
       {
-        return pointAfter.m();
+        if ( QgsWkbTypes::hasZ( geometry.wkbType() ) && length )
+          closestPoint.addZValue( pointBefore.z() + ( pointAfter.z() - pointBefore.z() ) * distance / length );
+        if ( QgsWkbTypes::hasM( geometry.wkbType() ) )
+          closestPoint.addMValue( pointBefore.m() + ( pointAfter.m() - pointBefore.m() ) * distance / length );
       }
     }
   }
-  return std::numeric_limits<double>::quiet_NaN();
+
+  return closestPoint;
 }
 
 double QgsGeometryUtils::distanceToVertex( const QgsAbstractGeometry &geom, QgsVertexId id )

src/core/geometry/qgsgeometryutils.h

@@ -44,9 +44,12 @@ class CORE_EXPORT QgsGeometryUtils
      */
     static QgsPoint closestVertex( const QgsAbstractGeometry &geom, const QgsPoint &pt, QgsVertexId &id SIP_OUT );
 
-    /** Returns measure of nearest point on a geometry for a specified point or NaN if geometry does not have measures
+    /**
+     * Returns the nearest point on a segment of a \a geometry
+     * for the specified \a point. The z and m values will be linearly interpolated between
+     * the two neighbouring vertices.
      */
-    static double closestPointMeasure( const QgsAbstractGeometry &geom, const QgsPoint &pt );
+    static QgsPoint closestPoint( const QgsAbstractGeometry &geometry, const QgsPoint &point );
 
     /** Returns the distance along a geometry from its first vertex to the specified vertex.
      * \param geom geometry

src/core/geometry/qgslinestring.cpp

@@ -380,8 +380,8 @@ QgsPoint QgsLineString::pointN( int i ) const
 
   double x = mX.at( i );
   double y = mY.at( i );
-  double z = 0;
-  double m = 0;
+  double z = std::numeric_limits<double>::quiet_NaN();
+  double m = std::numeric_limits<double>::quiet_NaN();
 
   bool hasZ = is3D();
   if ( hasZ )
@@ -441,15 +441,15 @@ double QgsLineString::zAt( int index ) const
   if ( index >= 0 && index < mZ.size() )
     return mZ.at( index );
   else
-    return 0.0;
+    return std::numeric_limits<double>::quiet_NaN();
 }
 
 double QgsLineString::mAt( int index ) const
 {
   if ( index >= 0 && index < mM.size() )
     return mM.at( index );
   else
-    return 0.0;
+    return std::numeric_limits<double>::quiet_NaN();
 }
 
 void QgsLineString::setXAt( int index, double x )

src/core/geometry/qgslinestring.h

@@ -44,6 +44,8 @@ class CORE_EXPORT QgsLineString: public QgsCurve
 
     /**
      * Construct a linestring from a vector of points.
+     * Z and M type will be set based on the type of the first point
+     * in the vector.
      * \since QGIS 3.0
      */
     QgsLineString( const QVector<QgsPoint> &points );
@@ -81,15 +83,15 @@ class CORE_EXPORT QgsLineString: public QgsCurve
 
     /** Returns the z-coordinate of the specified node in the line string.
      * \param index index of node, where the first node in the line is 0
-     * \returns z-coordinate of node, or 0.0 if index is out of bounds or the line
+     * \returns z-coordinate of node, or ``nan`` if index is out of bounds or the line
      * does not have a z dimension
      * \see setZAt()
      */
     double zAt( int index ) const;
 
     /** Returns the m value of the specified node in the line string.
      * \param index index of node, where the first node in the line is 0
-     * \returns m value of node, or 0.0 if index is out of bounds or the line
+     * \returns m value of node, or ``nan`` if index is out of bounds or the line
      * does not have m values
      * \see setMAt()
      */

src/gui/qgsmaptoolidentify.cpp

@@ -332,8 +332,8 @@ void QgsMapToolIdentify::closestPointAttributes( const QgsAbstractGeometry &geom
   // measure
   if ( QgsWkbTypes::hasM( geometry.wkbType() ) )
   {
-    double measure = QgsGeometryUtils::closestPointMeasure( geometry, QgsPoint( layerPoint.x(), layerPoint.y() ) );
-    QString str = QLocale::system().toString( measure, 'g', 10 );
+    QgsPoint closestPoint = QgsGeometryUtils::closestPoint( geometry, QgsPoint( layerPoint.x(), layerPoint.y() ) );
+    QString str = QLocale::system().toString( closestPoint.m(), 'g', 10 );
     derivedAttributes.insert( QStringLiteral( "Closest point M" ), str );
   }
 }

tests/src/core/testqgsgeometry.cpp

@@ -1317,8 +1317,8 @@ void TestQgsGeometry::lineString()
   QCOMPARE( l9.zAt( 0 ), 3.0 );
   QCOMPARE( l9.zAt( 1 ), 13.0 );
   QCOMPARE( l9.zAt( 2 ), 23.0 );
-  QCOMPARE( l9.zAt( -1 ), 0.0 ); //out of range
-  QCOMPARE( l9.zAt( 11 ), 0.0 ); //out of range
+  QVERIFY( qIsNaN( l9.zAt( -1 ) ) ); //out of range
+  QVERIFY( qIsNaN( l9.zAt( 11 ) ) ); //out of range
 
   l9.setZAt( 0, 53.0 );
   QCOMPARE( l9.zAt( 0 ), 53.0 );
@@ -1330,8 +1330,8 @@ void TestQgsGeometry::lineString()
   QCOMPARE( l9.mAt( 0 ), 4.0 );
   QCOMPARE( l9.mAt( 1 ), 14.0 );
   QCOMPARE( l9.mAt( 2 ), 24.0 );
-  QCOMPARE( l9.mAt( -1 ), 0.0 ); //out of range
-  QCOMPARE( l9.mAt( 11 ), 0.0 ); //out of range
+  QVERIFY( qIsNaN( l9.mAt( -1 ) ) ); //out of range
+  QVERIFY( qIsNaN( l9.mAt( 11 ) ) ); //out of range
 
   l9.setMAt( 0, 54.0 );
   QCOMPARE( l9.mAt( 0 ), 54.0 );
@@ -1346,8 +1346,8 @@ void TestQgsGeometry::lineString()
                 << QgsPoint( QgsWkbTypes::PointM, 21, 22, 0, 24 ) );
 
   //basically we just don't want these to crash
-  QCOMPARE( l9.zAt( 0 ), 0.0 );
-  QCOMPARE( l9.zAt( 1 ), 0.0 );
+  QVERIFY( qIsNaN( l9.zAt( 0 ) ) );
+  QVERIFY( qIsNaN( l9.zAt( 1 ) ) );
   l9.setZAt( 0, 53.0 );
   l9.setZAt( 1, 63.0 );
 
@@ -1357,8 +1357,8 @@ void TestQgsGeometry::lineString()
                 << QgsPoint( 21, 22 ) );
 
   //basically we just don't want these to crash
-  QCOMPARE( l9.mAt( 0 ), 0.0 );
-  QCOMPARE( l9.mAt( 1 ), 0.0 );
+  QVERIFY( qIsNaN( l9.mAt( 0 ) ) );
+  QVERIFY( qIsNaN( l9.mAt( 1 ) ) );
   l9.setMAt( 0, 53.0 );
   l9.setMAt( 1, 63.0 );
 

tests/src/core/testqgsgeometryutils.cpp

@@ -55,6 +55,7 @@ class TestQgsGeometryUtils: public QObject
     void testGradient();
     void testCoefficients();
     void testPerpendicularSegment();
+    void testClosestPoint();
 };
 
 
@@ -633,5 +634,41 @@ void TestQgsGeometryUtils::testPerpendicularSegment()
   QCOMPARE( line.pointN( 1 ), line_r.pointN( 1 ) );
 }
 
+void TestQgsGeometryUtils::testClosestPoint()
+{
+  QgsLineString linestringZ( QVector<QgsPoint>()
+                             << QgsPoint( 1, 1, 1 )
+                             << QgsPoint( 1, 3, 2 ) );
+
+  QgsPoint pt1 = QgsGeometryUtils::closestPoint( linestringZ, QgsPoint( 1, 0 ) );
+  QGSCOMPARENEAR( pt1.z(), 1, 0.0001 );
+  QVERIFY( qIsNaN( pt1.m() ) );
+
+  QgsLineString linestringM( QVector<QgsPoint>()
+                             << QgsPoint( 1, 1, std::numeric_limits<double>::quiet_NaN(), 1 )
+                             << QgsPoint( 1, 3, std::numeric_limits<double>::quiet_NaN(), 2 ) );
+
+  QgsPoint pt2 = QgsGeometryUtils::closestPoint( linestringM, QgsPoint( 1, 4 ) );
+  QVERIFY( qIsNaN( pt2.z() ) );
+  QGSCOMPARENEAR( pt2.m(), 2, 0.0001 );
+
+  QgsLineString linestringZM( QVector<QgsPoint>()
+                              << QgsPoint( 1, 1, 1, 1 )
+                              << QgsPoint( 1, 3, 2, 2 ) );
+
+  QgsPoint pt3 = QgsGeometryUtils::closestPoint( linestringZM, QgsPoint( 2, 2 ) );
+  QGSCOMPARENEAR( pt3.z(), 1.5, 0.0001 );
+  QGSCOMPARENEAR( pt3.m(), 1.5, 0.0001 );
+
+  QgsLineString linestringDuplicatedPoint( QVector<QgsPoint>()
+      << QgsPoint( 1, 1, 1, 1 )
+      << QgsPoint( 1, 1, 1, 1 )
+      << QgsPoint( 1, 3, 2, 2 ) );
+
+  QgsPoint pt4 = QgsGeometryUtils::closestPoint( linestringDuplicatedPoint, QgsPoint( 1, 0 ) );
+  QGSCOMPARENEAR( pt4.z(), 1, 0.0001 );
+  QGSCOMPARENEAR( pt4.m(), 1, 0.0001 );
+}
+
 QGSTEST_MAIN( TestQgsGeometryUtils )
 #include "testqgsgeometryutils.moc"