diff --git a/python/core/auto_generated/geometry/qgsgeometryengine.sip.in b/python/core/auto_generated/geometry/qgsgeometryengine.sip.in
index bac69b84c67..7d747e3d299 100644
--- a/python/core/auto_generated/geometry/qgsgeometryengine.sip.in
+++ b/python/core/auto_generated/geometry/qgsgeometryengine.sip.in
@@ -166,6 +166,13 @@ Calculate the convex hull of this.
 Calculates the distance between this and ``geom``.
 
 .. versionadded:: 3.0
+%End
+
+    virtual bool distanceWithin( const QgsAbstractGeometry *geom, double maxdistance, QString *errorMsg = 0 ) const = 0;
+%Docstring
+Checks if ``geom`` is within ``maxdistance`` distance from this geometry
+
+.. versionadded:: 3.22
 %End
 
     virtual bool intersects( const QgsAbstractGeometry *geom, QString *errorMsg = 0 ) const = 0;
diff --git a/src/core/geometry/qgsgeometryengine.h b/src/core/geometry/qgsgeometryengine.h
index 23255c7e22c..543fcc26a3a 100644
--- a/src/core/geometry/qgsgeometryengine.h
+++ b/src/core/geometry/qgsgeometryengine.h
@@ -184,6 +184,13 @@ class CORE_EXPORT QgsGeometryEngine
      */
     virtual double distance( const QgsAbstractGeometry *geom, QString *errorMsg = nullptr ) const = 0;
 
+    /**
+     * Checks if \a geom is within \a maxdistance distance from this geometry
+     *
+     * \since QGIS 3.22
+     */
+    virtual bool distanceWithin( const QgsAbstractGeometry *geom, double maxdistance, QString *errorMsg = nullptr ) const = 0;
+
     /**
      * Checks if \a geom intersects this.
      *
diff --git a/src/core/geometry/qgsgeos.cpp b/src/core/geometry/qgsgeos.cpp
index a5913fec753..d6a7ededfa6 100644
--- a/src/core/geometry/qgsgeos.cpp
+++ b/src/core/geometry/qgsgeos.cpp
@@ -475,6 +475,44 @@ double QgsGeos::distance( const QgsAbstractGeometry *geom, QString *errorMsg ) c
   return distance;
 }
 
+bool QgsGeos::distanceWithin( const QgsAbstractGeometry *geom, double maxdist, QString *errorMsg ) const
+{
+  if ( !mGeos )
+  {
+    return false;
+  }
+
+  geos::unique_ptr otherGeosGeom( asGeos( geom, mPrecision ) );
+  if ( !otherGeosGeom )
+  {
+    return false;
+  }
+
+  // TODO: optimize implementation of this function to early-exit if
+  // any part of othergeosGeom is found to be within the given
+  // distance
+
+  double distance;
+  try
+  {
+#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=9 )
+    if ( mGeosPrepared )
+    {
+      GEOSPreparedDistance_r( geosinit()->ctxt, mGeosPrepared.get(), otherGeosGeom.get(), &distance );
+    }
+    else
+    {
+      GEOSDistance_r( geosinit()->ctxt, mGeos.get(), otherGeosGeom.get(), &distance );
+    }
+#else
+    GEOSDistance_r( geosinit()->ctxt, mGeos.get(), otherGeosGeom.get(), &distance );
+#endif
+  }
+  CATCH_GEOS_WITH_ERRMSG( false )
+
+  return distance <= maxdist;
+}
+
 double QgsGeos::hausdorffDistance( const QgsAbstractGeometry *geom, QString *errorMsg ) const
 {
   double distance = -1.0;
diff --git a/src/core/geometry/qgsgeos.h b/src/core/geometry/qgsgeos.h
index ed4c1f2b0ce..84622cdd81c 100644
--- a/src/core/geometry/qgsgeos.h
+++ b/src/core/geometry/qgsgeos.h
@@ -180,6 +180,7 @@ class CORE_EXPORT QgsGeos: public QgsGeometryEngine
     QgsPoint *pointOnSurface( QString *errorMsg = nullptr ) const override;
     QgsAbstractGeometry *convexHull( QString *errorMsg = nullptr ) const override;
     double distance( const QgsAbstractGeometry *geom, QString *errorMsg = nullptr ) const override;
+    bool distanceWithin( const QgsAbstractGeometry *geom, double maxdistance, QString *errorMsg = nullptr ) const override;
 
     /**
      * Returns the Hausdorff distance between this geometry and \a geom. This is basically a measure of how similar or dissimilar 2 geometries are.
diff --git a/src/core/vector/qgsvectorlayerfeatureiterator.cpp b/src/core/vector/qgsvectorlayerfeatureiterator.cpp
index 61685c81e11..35593d9cd26 100644
--- a/src/core/vector/qgsvectorlayerfeatureiterator.cpp
+++ b/src/core/vector/qgsvectorlayerfeatureiterator.cpp
@@ -903,7 +903,7 @@ bool QgsVectorLayerFeatureIterator::postProcessFeature( QgsFeature &feature )
 
   if ( result && mDistanceWithinEngine && feature.hasGeometry() )
   {
-    result = mDistanceWithinEngine->distance( feature.geometry().constGet() ) <= mDistanceWithin;
+    result = mDistanceWithinEngine->distanceWithin( feature.geometry().constGet(), mDistanceWithin );
   }
 
   return result;