Super-optimised version of geometry bounding box intersects test

Apply some fancy logic to make this test as cheap as possible
to run

python/core/auto_generated/geometry/qgsabstractgeometry.sip.in

@@ -521,6 +521,16 @@ Returns ``True`` if the geometry is empty
     virtual bool hasCurvedSegments() const;
 %Docstring
 Returns ``True`` if the geometry contains curved segments
+%End
+
+    virtual bool boundingBoxIntersects( const QgsRectangle &rectangle ) const /HoldGIL/;
+%Docstring
+Returns ``True`` if the bounding box of this geometry intersects with a ``rectangle``.
+
+Since this test only considers the bounding box of the geometry, is is very fast to
+calculate and handles invalid geometries.
+
+.. versionadded:: 3.20
 %End
 
     virtual QgsAbstractGeometry *segmentize( double tolerance = M_PI / 180., SegmentationToleranceType toleranceType = MaximumAngle ) const /Factory/;

python/core/auto_generated/geometry/qgscompoundcurve.sip.in

@@ -82,6 +82,8 @@ of the curve.
 
     virtual bool removeDuplicateNodes( double epsilon = 4 * DBL_EPSILON, bool useZValues = false );
 
+    virtual bool boundingBoxIntersects( const QgsRectangle &rectangle ) const /HoldGIL/;
+
 
     int nCurves() const /HoldGIL/;
 %Docstring

python/core/auto_generated/geometry/qgscurve.sip.in

@@ -280,6 +280,7 @@ Returns the curve's orientation, e.g. clockwise or counter-clockwise.
 
 
 
+
 };
 
 /************************************************************************

python/core/auto_generated/geometry/qgscurvepolygon.sip.in

@@ -71,6 +71,8 @@ Curve polygon geometry type
 
     virtual bool removeDuplicateNodes( double epsilon = 4 * DBL_EPSILON, bool useZValues = false );
 
+    virtual bool boundingBoxIntersects( const QgsRectangle &rectangle ) const /HoldGIL/;
+
 
 
     int numInteriorRings() const /HoldGIL/;

python/core/auto_generated/geometry/qgsgeometrycollection.sip.in

@@ -102,6 +102,8 @@ Returns a geometry from within the collection.
 
     virtual int vertexNumberFromVertexId( QgsVertexId id ) const;
 
+    virtual bool boundingBoxIntersects( const QgsRectangle &rectangle ) const /HoldGIL/;
+
 
     void reserve( int size ) /HoldGIL/;
 %Docstring

python/core/auto_generated/geometry/qgslinestring.sip.in

@@ -424,6 +424,8 @@ segment in the line.
 
     virtual bool isClosed() const /HoldGIL/;
 
+    virtual bool boundingBoxIntersects( const QgsRectangle &rectangle ) const /HoldGIL/;
+
 
     QVector< QgsVertexId > collectDuplicateNodes( double epsilon = 4 * DBL_EPSILON, bool useZValues = false ) const;
 %Docstring

python/core/auto_generated/geometry/qgspoint.sip.in

@@ -436,6 +436,8 @@ Angle undefined. Always returns 0.0
 
     virtual double segmentLength( QgsVertexId startVertex ) const;
 
+    virtual bool boundingBoxIntersects( const QgsRectangle &rectangle ) const /HoldGIL/;
+
 
     virtual bool addZValue( double zValue = 0 );
 

src/core/geometry/qgsabstractgeometry.cpp

@@ -307,6 +307,11 @@ bool QgsAbstractGeometry::hasCurvedSegments() const
   return false;
 }
 
+bool QgsAbstractGeometry::boundingBoxIntersects( const QgsRectangle &rectangle ) const
+{
+  return boundingBox().intersects( rectangle );
+}
+
 QgsAbstractGeometry *QgsAbstractGeometry::segmentize( double tolerance, SegmentationToleranceType toleranceType ) const
 {
   Q_UNUSED( tolerance )

src/core/geometry/qgsabstractgeometry.h

@@ -537,6 +537,16 @@ class CORE_EXPORT QgsAbstractGeometry
      */
     virtual bool hasCurvedSegments() const;
 
+    /**
+     * Returns TRUE if the bounding box of this geometry intersects with a \a rectangle.
+     *
+     * Since this test only considers the bounding box of the geometry, is is very fast to
+     * calculate and handles invalid geometries.
+     *
+     * \since QGIS 3.20
+     */
+    virtual bool boundingBoxIntersects( const QgsRectangle &rectangle ) const SIP_HOLDGIL;
+
     /**
      * Returns a version of the geometry without curves. Caller takes ownership of
      * the returned geometry.

src/core/geometry/qgscompoundcurve.cpp

@@ -465,6 +465,35 @@ bool QgsCompoundCurve::removeDuplicateNodes( double epsilon, bool useZValues )
   return result;
 }
 
+bool QgsCompoundCurve::boundingBoxIntersects( const QgsRectangle &rectangle ) const
+{
+  if ( mCurves.empty() )
+    return false;
+
+  // if we already have the bounding box calculated, then this check is trivial!
+  if ( !mBoundingBox.isNull() )
+  {
+    return mBoundingBox.intersects( rectangle );
+  }
+
+  // otherwise loop through each member curve and test the bounding box intersection.
+  // This gives us a chance to use optimisations which may be present on the individual
+  // curve subclasses, and at worst it will cause a calculation of the bounding box
+  // of each individual member curve which we would have to do anyway... (and these
+  // bounding boxes are cached, so would be reused without additional expense)
+  for ( const QgsCurve *curve : mCurves )
+  {
+    if ( curve->boundingBoxIntersects( rectangle ) )
+      return true;
+  }
+
+  // even if we don't intersect the bounding box of any member curves, we may still intersect the
+  // bounding box of the overall compound curve.
+  // so here we fall back to the non-optimised base class check which has to first calculate
+  // the overall bounding box of the compound curve..
+  return QgsAbstractGeometry::boundingBoxIntersects( rectangle );
+}
+
 const QgsCurve *QgsCompoundCurve::curveAt( int i ) const
 {
   if ( i < 0 || i >= mCurves.size() )

src/core/geometry/qgscompoundcurve.h

@@ -72,6 +72,7 @@ class CORE_EXPORT QgsCompoundCurve: public QgsCurve
 
     QgsCompoundCurve *snappedToGrid( double hSpacing, double vSpacing, double dSpacing = 0, double mSpacing = 0 ) const override SIP_FACTORY;
     bool removeDuplicateNodes( double epsilon = 4 * std::numeric_limits<double>::epsilon(), bool useZValues = false ) override;
+    bool boundingBoxIntersects( const QgsRectangle &rectangle ) const override SIP_HOLDGIL;
 
     /**
      * Returns the number of curves in the geometry.

src/core/geometry/qgscurve.h

@@ -292,10 +292,13 @@ class CORE_EXPORT QgsCurve: public QgsAbstractGeometry
                             QVector<double> &outX, QVector<double> &outY, QVector<double> &outZ, QVector<double> &outM ) const;
 #endif
 
-  private:
-
+    /**
+     * Cached bounding box.
+     */
     mutable QgsRectangle mBoundingBox;
 
+  private:
+
     mutable bool mHasCachedValidity = false;
     mutable QString mValidityFailureReason;
 };

src/core/geometry/qgscurvepolygon.cpp

@@ -612,6 +612,39 @@ bool QgsCurvePolygon::removeDuplicateNodes( double epsilon, bool useZValues )
   return result;
 }
 
+bool QgsCurvePolygon::boundingBoxIntersects( const QgsRectangle &rectangle ) const
+{
+  if ( !mExteriorRing && mInteriorRings.empty() )
+    return false;
+
+  // if we already have the bounding box calculated, then this check is trivial!
+  if ( !mBoundingBox.isNull() )
+  {
+    return mBoundingBox.intersects( rectangle );
+  }
+
+  // loop through each ring and test the bounding box intersection.
+  // This gives us a chance to use optimisations which may be present on the individual
+  // ring geometry subclasses, and at worst it will cause a calculation of the bounding box
+  // of each individual ring geometry which we would have to do anyway... (and these
+  // bounding boxes are cached, so would be reused without additional expense)
+  if ( mExteriorRing && mExteriorRing->boundingBoxIntersects( rectangle ) )
+    return true;
+
+  for ( const QgsCurve *ring : mInteriorRings )
+  {
+    if ( ring->boundingBoxIntersects( rectangle ) )
+      return true;
+  }
+
+  // even if we don't intersect the bounding box of any rings, we may still intersect the
+  // bounding box of the overall polygon (we are considering worst case scenario here and
+  // the polygon is invalid, with rings outside the exterior ring!)
+  // so here we fall back to the non-optimised base class check which has to first calculate
+  // the overall bounding box of the polygon..
+  return QgsSurface::boundingBoxIntersects( rectangle );
+}
+
 QgsPolygon *QgsCurvePolygon::toPolygon( double tolerance, SegmentationToleranceType toleranceType ) const
 {
   std::unique_ptr< QgsPolygon > poly( new QgsPolygon() );

src/core/geometry/qgscurvepolygon.h

@@ -66,6 +66,7 @@ class CORE_EXPORT QgsCurvePolygon: public QgsSurface
     QgsAbstractGeometry *boundary() const override SIP_FACTORY;
     QgsCurvePolygon *snappedToGrid( double hSpacing, double vSpacing, double dSpacing = 0, double mSpacing = 0 ) const override SIP_FACTORY;
     bool removeDuplicateNodes( double epsilon = 4 * std::numeric_limits<double>::epsilon(), bool useZValues = false ) override;
+    bool boundingBoxIntersects( const QgsRectangle &rectangle ) const override SIP_HOLDGIL;
 
     //curve polygon interface
 

src/core/geometry/qgsgeometry.cpp

@@ -1160,14 +1160,7 @@ bool QgsGeometry::boundingBoxIntersects( const QgsRectangle &rectangle ) const
     return false;
   }
 
-  // optimise trivial case for point intersections
-  if ( QgsWkbTypes::flatType( d->geometry->wkbType() ) == QgsWkbTypes::Point )
-  {
-    const QgsPoint *point = qgsgeometry_cast< const QgsPoint * >( d->geometry.get() );
-    return rectangle.contains( QgsPointXY( point->x(), point->y() ) );
-  }
-
-  return d->geometry->boundingBox().intersects( rectangle );
+  return d->geometry->boundingBoxIntersects( rectangle );
 }
 
 bool QgsGeometry::boundingBoxIntersects( const QgsGeometry &geometry ) const
@@ -1177,7 +1170,7 @@ bool QgsGeometry::boundingBoxIntersects( const QgsGeometry &geometry ) const
     return false;
   }
 
-  return d->geometry->boundingBox().intersects( geometry.constGet()->boundingBox() );
+  return d->geometry->boundingBoxIntersects( geometry.constGet()->boundingBox() );
 }
 
 bool QgsGeometry::contains( const QgsPointXY *p ) const

src/core/geometry/qgsgeometrycollection.cpp

@@ -200,6 +200,35 @@ int QgsGeometryCollection::vertexNumberFromVertexId( QgsVertexId id ) const
   return -1; // should not happen
 }
 
+bool QgsGeometryCollection::boundingBoxIntersects( const QgsRectangle &rectangle ) const
+{
+  if ( mGeometries.empty() )
+    return false;
+
+  // if we already have the bounding box calculated, then this check is trivial!
+  if ( !mBoundingBox.isNull() )
+  {
+    return mBoundingBox.intersects( rectangle );
+  }
+
+  // otherwise loop through each member geometry and test the bounding box intersection.
+  // This gives us a chance to use optimisations which may be present on the individual
+  // geometry subclasses, and at worst it will cause a calculation of the bounding box
+  // of each individual member geometry which we would have to do anyway... (and these
+  // bounding boxes are cached, so would be reused without additional expense)
+  for ( const QgsAbstractGeometry *geometry : mGeometries )
+  {
+    if ( geometry->boundingBoxIntersects( rectangle ) )
+      return true;
+  }
+
+  // even if we don't intersect the bounding box of any member geometries, we may still intersect the
+  // bounding box of the overall collection.
+  // so here we fall back to the non-optimised base class check which has to first calculate
+  // the overall bounding box of the collection..
+  return QgsAbstractGeometry::boundingBoxIntersects( rectangle );
+}
+
 void QgsGeometryCollection::reserve( int size )
 {
   mGeometries.reserve( size );

src/core/geometry/qgsgeometrycollection.h

@@ -125,6 +125,7 @@ class CORE_EXPORT QgsGeometryCollection: public QgsAbstractGeometry
     QgsAbstractGeometry *boundary() const override SIP_FACTORY;
     void adjacentVertices( QgsVertexId vertex, QgsVertexId &previousVertex SIP_OUT, QgsVertexId &nextVertex SIP_OUT ) const override;
     int vertexNumberFromVertexId( QgsVertexId id ) const override;
+    bool boundingBoxIntersects( const QgsRectangle &rectangle ) const override SIP_HOLDGIL;
 
     /**
      * Attempts to allocate memory for at least \a size geometries.

src/core/geometry/qgslinestring.cpp

@@ -387,6 +387,84 @@ bool QgsLineString::isClosed() const
   return closed;
 }
 
+bool QgsLineString::boundingBoxIntersects( const QgsRectangle &rectangle ) const
+{
+  if ( mX.empty() )
+    return false;
+
+  if ( !mBoundingBox.isNull() )
+  {
+    return mBoundingBox.intersects( rectangle );
+  }
+  const int nb = mX.size();
+
+  // We are a little fancy here!
+  if ( nb > 40 )
+  {
+    // if a large number of vertices, take some sample vertices at 1/5th increments through the linestring
+    // and test whether any are inside the rectangle. Maybe we can shortcut a lot of iterations by doing this!
+    // (why 1/5th? it's picked so that it works nicely for polygon rings which are almost rectangles, so the vertex extremities
+    // will fall on approximately these vertex indices)
+    if ( rectangle.contains( mX.at( 0 ), mY.at( 0 ) ) ||
+         rectangle.contains( mX.at( static_cast< int >( nb * 0.2 ) ), mY.at( static_cast< int >( nb * 0.2 ) ) ) ||
+         rectangle.contains( mX.at( static_cast< int >( nb * 0.4 ) ), mY.at( static_cast< int >( nb * 0.4 ) ) ) ||
+         rectangle.contains( mX.at( static_cast< int >( nb * 0.6 ) ), mY.at( static_cast< int >( nb * 0.6 ) ) ) ||
+         rectangle.contains( mX.at( static_cast< int >( nb * 0.8 ) ), mY.at( static_cast< int >( nb * 0.8 ) ) ) ||
+         rectangle.contains( mX.at( nb - 1 ), mY.at( nb - 1 ) ) )
+      return true;
+  }
+
+  // Be even MORE fancy! Given that bounding box calculation is non-free, cached, and we don't
+  // already have it, we start performing the bounding box calculation while we are testing whether
+  // each point falls inside the rectangle. That way if we end up testing the majority of the points
+  // anyway, we can update the cached bounding box with the results we've calculated along the way
+  // and save future calls to calculate the bounding box!
+  double xmin = std::numeric_limits<double>::max();
+  double ymin = std::numeric_limits<double>::max();
+  double xmax = -std::numeric_limits<double>::max();
+  double ymax = -std::numeric_limits<double>::max();
+
+  const double *x = mX.constData();
+  const double *y = mY.constData();
+  bool foundPointInRectangle = false;
+  for ( int i = 0; i < nb; ++i )
+  {
+    const double px = *x++;
+    xmin = std::min( xmin, px );
+    xmax = std::max( xmax, px );
+    const double py = *y++;
+    ymin = std::min( ymin, py );
+    ymax = std::max( ymax, py );
+
+    if ( !foundPointInRectangle && rectangle.contains( px, py ) )
+    {
+      foundPointInRectangle = true;
+
+      // now... we have a choice to make. If we've already looped through the majority of the points
+      // in this linestring then let's just continue to iterate through the remainder so that we can
+      // complete the overall bounding box calculation we've already mostly done. If however we're only
+      // just at the start of iterating the vertices, we shortcut out early and leave the bounding box
+      // uncalculated
+      if ( i < nb * 0.5 )
+        return true;
+    }
+  }
+
+  // at this stage we now know the overall bounding box of the linestring, so let's cache
+  // it so we don't ever have to calculate this again. We've done all the hard work anyway!
+  mBoundingBox = QgsRectangle( xmin, ymin, xmax, ymax, false );
+
+  if ( foundPointInRectangle )
+    return true;
+
+  // NOTE: if none of the points in the line actually fell inside the rectangle, it doesn't
+  // exclude that the OVERALL bounding box of the linestring itself intersects the rectangle!!
+  // So we fall back to the parent class method which compares the overall bounding box against
+  // the rectangle... and this will be very cheap now that we've already calculated and cached
+  // the linestring's bounding box!
+  return QgsCurve::boundingBoxIntersects( rectangle );
+}
+
 QVector< QgsVertexId > QgsLineString::collectDuplicateNodes( double epsilon, bool useZValues ) const
 {
   QVector< QgsVertexId > res;

src/core/geometry/qgslinestring.h

@@ -591,6 +591,7 @@ class CORE_EXPORT QgsLineString: public QgsCurve
     QgsLineString *snappedToGrid( double hSpacing, double vSpacing, double dSpacing = 0, double mSpacing = 0 ) const override SIP_FACTORY;
     bool removeDuplicateNodes( double epsilon = 4 * std::numeric_limits<double>::epsilon(), bool useZValues = false ) override;
     bool isClosed() const override SIP_HOLDGIL;
+    bool boundingBoxIntersects( const QgsRectangle &rectangle ) const override SIP_HOLDGIL;
 
     /**
      * Returns a list of any duplicate nodes contained in the geometry, within the specified tolerance.

src/core/geometry/qgspoint.cpp

@@ -537,6 +537,11 @@ double QgsPoint::segmentLength( QgsVertexId ) const
   return 0.0;
 }
 
+bool QgsPoint::boundingBoxIntersects( const QgsRectangle &rectangle ) const
+{
+  return rectangle.contains( mX, mY );
+}
+
 /***************************************************************************
  * This class is considered CRITICAL and any change MUST be accompanied with
  * full unit tests.

src/core/geometry/qgspoint.h

@@ -529,6 +529,7 @@ class CORE_EXPORT QgsPoint: public QgsAbstractGeometry
     QgsPoint vertexAt( QgsVertexId /*id*/ ) const override;
     QgsPoint *toCurveType() const override SIP_FACTORY;
     double segmentLength( QgsVertexId startVertex ) const override;
+    bool boundingBoxIntersects( const QgsRectangle &rectangle ) const override SIP_HOLDGIL;
 
     bool addZValue( double zValue = 0 ) override;
     bool addMValue( double mValue = 0 ) override;