fix(geometry): correct area calculation for QgsLineString orientation

qgis · Oct 21, 2024 · fc2a221 · fc2a221
1 parent 372ab66
commit fc2a221
Show file tree

Hide file tree

Showing 3 changed files with 35 additions and 13 deletions.
diff --git a/src/core/geometry/qgslinestring.cpp b/src/core/geometry/qgslinestring.cpp
@@ -2214,6 +2214,11 @@ QgsPoint QgsLineString::centroid() const
  * See details in QEP #17
  ****************************************************************************/
 
+double area2( const double ax, const double ay, const double bx, const double by, const double cx, const double cy )
+{
+  return ( bx - ax ) * ( cy - ay ) - ( cx - ax ) * ( by - ay );
+}
+
 void QgsLineString::sumUpArea( double &sum ) const
 {
   if ( mHasCachedSummedUpArea )
@@ -2224,24 +2229,30 @@ void QgsLineString::sumUpArea( double &sum ) const
 
   mSummedUpArea = 0;
   const int maxIndex = mX.size();
-  if ( maxIndex < 2 )
+  if ( maxIndex < 3 )
   {
     mHasCachedSummedUpArea = true;
     return;
   }
 
   const double *x = mX.constData();
   const double *y = mY.constData();
-  double prevX = *x++;
-  double prevY = *y++;
-  for ( int i = 1; i < maxIndex; ++i )
+  const double x0 = *x++;
+  const double y0 = *y++;
+
+  double x1 = *x++;
+  double y1 = *y++;
+
+  for ( int i = 2; i < maxIndex; ++i )
   {
-    mSummedUpArea += prevX * ( *y - prevY ) - prevY * ( *x - prevX );
-    prevX = *x++;
-    prevY = *y++;
+    const double x2 = *x++;
+    const double y2 = *y++;
+    mSummedUpArea += area2( x0, y0, x1, y1, x2, y2 );
+    x1 = x2;
+    y1 = y2;
   }
-  mSummedUpArea *= 0.5;
 
+  mSummedUpArea *= 0.5;
   mHasCachedSummedUpArea = true;
   sum += mSummedUpArea;
 }

diff --git a/tests/src/core/geometry/testqgslinestring.cpp b/tests/src/core/geometry/testqgslinestring.cpp
@@ -2234,24 +2234,23 @@ void TestQgsLineString::sumUpArea()
 
   ls.setPoints( QgsPointSequence() << QgsPoint( 5, 10 ) << QgsPoint( 10, 10 ) );
   ls.sumUpArea( area );
-  QGSCOMPARENEAR( area, -24, 4 * std::numeric_limits<double>::epsilon() );
+  QGSCOMPARENEAR( area, 1, 4 * std::numeric_limits<double>::epsilon() );
 
   ls.setPoints( QgsPointSequence() << QgsPoint( 0, 0 )
                 << QgsPoint( 2, 0 ) << QgsPoint( 2, 2 ) );
   ls.sumUpArea( area );
-  QGSCOMPARENEAR( area, -22, 4 * std::numeric_limits<double>::epsilon() );
-
+  QGSCOMPARENEAR( area, 3, 4 * std::numeric_limits<double>::epsilon() );
   ls.setPoints( QgsPointSequence() << QgsPoint( 0, 0 ) << QgsPoint( 2, 0 )
                 << QgsPoint( 2, 2 ) << QgsPoint( 0, 2 ) );
   ls.sumUpArea( area );
-  QGSCOMPARENEAR( area, -18, 4 * std::numeric_limits<double>::epsilon() );
+  QGSCOMPARENEAR( area, 7, 4 * std::numeric_limits<double>::epsilon() );
 
   double shift = 10.0;
   ls.setPoints( QgsPointSequence() << QgsPoint( shift + 0, shift + 0 ) << QgsPoint( shift + 2, shift + 0 )
                 << QgsPoint( shift + 2, shift + 2 ) << QgsPoint( shift + 0, shift + 2 )
                 << QgsPoint( shift + 0, shift + 0 ) );
   ls.sumUpArea( area );
-  QGSCOMPARENEAR( area, -14, 4 * std::numeric_limits<double>::epsilon() );
+  QGSCOMPARENEAR( area, 11, 4 * std::numeric_limits<double>::epsilon() );
 
   // the length of the equator ~ 40 075.014 172 304 363 km
   shift = 40075.014172304363;

diff --git a/tests/src/python/test_qgslinestring.py b/tests/src/python/test_qgslinestring.py
@@ -389,6 +389,18 @@ def test_simplify_by_distance(self):
         self.assertEqual(p.simplifyByDistance(0).asWkt(),
                          'LineString (2 0, 2 2, 0 2, 0 0, 2 0)')
 
+    def test_orientation(self):
+        """
+        test orientation. From https://github.com/qgis/QGIS/issues/58333
+        """
+        geom = QgsLineString()
+        geom.fromWkt('LineString (1 1, 2 1, 2 2, 1 2, 1 1)')
+        self.assertEqual(geom.sumUpArea(), 1.0)
+        self.assertEqual(geom.orientation(), Qgis.AngularDirection.CounterClockwise)
+        geom.fromWkt('LineString (1 1, 1 2, 2 2, 2 1, 1 1)')
+        self.assertEqual(geom.sumUpArea(), -1.0)
+        self.assertEqual(geom.orientation(), Qgis.AngularDirection.Clockwise)
+
 
 if __name__ == '__main__':
     unittest.main()