diff --git a/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H b/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H
index ce208cb2591ea72c4c7944c8def26bbd69949ee7..8a365e9ae81b9033383314d97b0f53bc02a8d58c 100644
--- a/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H
+++ b/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -59,10 +59,10 @@ Usage
         bounds   | Limit with bounding box                  | no  |
         surfaceType | Type of surface                       | yes |
         surfaceName | Name of surface in \c triSurface/     | no  | dict name
-        topology    | Topology filter (none/largestRegion/nearestPoints/proximity) | no  | none
+        topology    | Topology filter name                  | no  | none
         nearestPoints | Points for point-based segmentation | no  |
         maxDistance | Max search distance for nearestPoints | no  | GREAT
-        relProximity | Max limit of absolute vs normal distance | no  | 1e3
+        absProximity | Max proximity of face centres        | no  | 1e-5
     \endtable
 
 SourceFiles
diff --git a/src/sampling/surface/distanceSurface/distanceSurface.C b/src/sampling/surface/distanceSurface/distanceSurface.C
index be57d3b122b6aa294f6f73a6bc65c580d2f07265..ed9d6ec631031fea52e0c8443143a9cb4805a3a9 100644
--- a/src/sampling/surface/distanceSurface/distanceSurface.C
+++ b/src/sampling/surface/distanceSurface/distanceSurface.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -52,7 +52,9 @@ Foam::distanceSurface::topoFilterNames_
     { topologyFilterType::NONE, "none" },
     { topologyFilterType::LARGEST_REGION, "largestRegion" },
     { topologyFilterType::NEAREST_POINTS, "nearestPoints" },
-    { topologyFilterType::PROXIMITY, "proximity" },
+    { topologyFilterType::PROXIMITY_REGIONS, "proximityRegions" },
+    { topologyFilterType::PROXIMITY_FACES, "proximityFaces" },
+    { topologyFilterType::PROXIMITY_FACES, "proximity" },
 });
 
 
@@ -184,15 +186,15 @@ static inline void calcNormalDistance_filtered
 // the normal at the edge, which creates a zero-crossing
 // extending to infinity.
 //
-// Ad hoc treatment: require that the surface hit
-// point is within a somewhat generous bounding box
-// for the cell
+// Ad hoc treatment: require that the surface hit point is within a
+// somewhat generous bounding box for the cell (+10%)
 template<bool WantPointFilter = false>
 static bitSet simpleGeometricFilter
 (
     bitSet& ignoreCells,
     const List<pointIndexHit>& nearest,
-    const polyMesh& mesh
+    const polyMesh& mesh,
+    const scalar boundBoxInflate = 0.1  // 10% to catch corners
 )
 {
     // A deny filter. Initially false (accept everything)
@@ -215,9 +217,7 @@ static bitSet simpleGeometricFilter
 
         cellBb.clear();
         cellBb.add(mesh.points(), cPoints);
-
-        // Expand slightly to catch corners
-        cellBb.inflate(0.1);
+        cellBb.inflate(boundBoxInflate);
 
         if (!cellBb.contains(pt))
         {
@@ -455,10 +455,30 @@ void Foam::distanceSurface::createGeometry()
         checkAllHits(nearest);
 
         // Geometric pre-filtering when distance == 0
+
+        // NOTE (2021-05-31)
+        // Can skip the prefilter if we use proximity-regions filter anyhow
+        // but it makes the iso algorithm more expensive and doesn't help
+        // unless we start relying on area-based weighting for rejecting regions.
+
         if (filterCells)
         {
-            ignoreCellPoints =
-                simpleGeometricFilter<false>(ignoreCells, nearest, fvmesh);
+            ignoreCellPoints = simpleGeometricFilter<false>
+            (
+                ignoreCells,
+                nearest,
+                fvmesh,
+
+                // Inflate bound box.
+                // - To catch corners: approx. 10%
+                // - Extra generous for PROXIMITY_REGIONS
+                //   (extra weighting for 'bad' faces)
+                (
+                    topologyFilterType::PROXIMITY_REGIONS == topoFilter_
+                  ? 1
+                  : 0.1
+                )
+            );
         }
 
         if (withSignDistance_)
@@ -617,6 +637,8 @@ void Foam::distanceSurface::createGeometry()
             refineBlockedCells(ignoreCells, isoCutter);
             filterKeepNearestRegions(ignoreCells);
         }
+
+        // Note: apply similar filtering for PROXIMITY_REGIONS later instead
     }
 
     // Don't need point filter beyond this point
@@ -627,6 +649,7 @@ void Foam::distanceSurface::createGeometry()
     {
         Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl;
         cellDistance.write();
+
         pointScalarField pDist
         (
             IOobject
@@ -660,6 +683,21 @@ void Foam::distanceSurface::createGeometry()
     );
 
 
+    // Restrict ignored cells to those actually cut
+    if (filterCells && topologyFilterType::PROXIMITY_REGIONS == topoFilter_)
+    {
+        isoSurfaceBase isoCutter
+        (
+            mesh_,
+            cellDistance,
+            pointDistance_,
+            distance_
+        );
+
+        refineBlockedCells(ignoreCells, isoCutter);
+    }
+
+
     // ALGO_POINT still needs cell, point fields (for interpolate)
     // The others can do straight transfer
 
@@ -667,7 +705,8 @@ void Foam::distanceSurface::createGeometry()
     if
     (
         isoParams_.algorithm() != isoSurfaceParams::ALGO_POINT
-     || topologyFilterType::PROXIMITY == topoFilter_
+     || topologyFilterType::PROXIMITY_FACES == topoFilter_
+     || topologyFilterType::PROXIMITY_REGIONS == topoFilter_
     )
     {
         surface_.transfer(static_cast<meshedSurface&>(*isoSurfacePtr_));
@@ -678,9 +717,13 @@ void Foam::distanceSurface::createGeometry()
         pointDistance_.clear();
     }
 
-    if (topologyFilterType::PROXIMITY == topoFilter_)
+    if (topologyFilterType::PROXIMITY_FACES == topoFilter_)
+    {
+        filterFaceProximity();
+    }
+    else if (topologyFilterType::PROXIMITY_REGIONS == topoFilter_)
     {
-        filterByProximity();
+        filterRegionProximity(ignoreCells);
     }
 
     if (debug)
diff --git a/src/sampling/surface/distanceSurface/distanceSurface.H b/src/sampling/surface/distanceSurface/distanceSurface.H
index 22738393291b3783e2b36059adca04144f461a9c..4979257cce29c92ebc1e04d7276874d88a632d78 100644
--- a/src/sampling/surface/distanceSurface/distanceSurface.H
+++ b/src/sampling/surface/distanceSurface/distanceSurface.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -77,13 +77,26 @@ Usage
         bounds   | Limit with bounding box                  | no  |
         surfaceType | Type of surface                       | yes |
         surfaceName | Name of surface in \c triSurface/     | no  | dict name
-        topology    | Topology filter (none/largestRegion/nearestPoints/proximity) | no  | none
+        topology    | Topology filter name                  | no  | none
         nearestPoints | Points for point-based segmentation | no  |
         maxDistance | Max search distance for nearestPoints | no  | GREAT
         absProximity | Max proximity of face centres        | no  | 1e-5
     \endtable
 
-    Filtering (for zero-distance only)
+    Topology/Filtering (for zero-distance only)
+
+    - \c none :
+        No filtering
+
+    - \c proximityFaces or \c proximity (post-filter):
+        Checks the resulting faces against the original search surface
+        and rejects faces with a distance greater than \c absProximity.
+
+    - \c proxityRegions (post-filter):
+        Checks the resulting faces against the original search surface
+        as well as checking the cut cells for topological connectivity.
+        If the area-weighted distance for a region is greater than
+        \c absProximity, the region is discarded.
 
     - \c largestRegion (pre-filter):
         The cut cells are checked for topological connectivity and the
@@ -94,10 +107,6 @@ Usage
         The cut cells split into regions, the regions closest to the
         user-defined points are retained.
         Uses \c maxDistance for additional control.
-
-    - \c proximity (post-filter):
-        Checks the resulting faces against the original search surface
-        and rejects faces with a distance greater than \c absProximity.
     .
 
 Note
@@ -146,7 +155,9 @@ class distanceSurface
             NONE,               //!< No additional filtering
             LARGEST_REGION,     //!< Retain largest region
             NEAREST_POINTS,     //!< Retain regions nearest to the points
-            PROXIMITY           //!< Post-filter for surface proximity
+            PROXIMITY_REGIONS,  //!< Retain regions with good surface proximity
+            PROXIMITY_FACES,    //!< Retain faces with good surface proximity
+            PROXIMITY = PROXIMITY_FACES
         };
 
         //- Names for the topology filter
@@ -262,15 +273,18 @@ protected:
         //- Prepare blockedFaces for region split
         bitSet filterPrepareRegionSplit(const bitSet& ignoreCells) const;
 
-        //- Adjust extracted iso-surface to remove far faces
-        void filterByProximity();
-
         //- Keep region with the most cuts (after region split)
         void filterKeepLargestRegion(bitSet& ignoreCells) const;
 
         //- Keep region(s) closest to the nearest points
         void filterKeepNearestRegions(bitSet& ignoreCells) const;
 
+        //- Remove region(s) that have far faces
+        void filterRegionProximity(bitSet& ignoreCells) const;
+
+        //- Adjust extracted iso-surface to remove far faces
+        void filterFaceProximity();
+
 
 public:
 
diff --git a/src/sampling/surface/distanceSurface/distanceSurfaceFilter.C b/src/sampling/surface/distanceSurface/distanceSurfaceFilter.C
index 9599058302400ac6dc8ee707ca9b6996880841fd..b007b38820a7e8169113ab5e5bc5776ce0cb9199 100644
--- a/src/sampling/surface/distanceSurface/distanceSurfaceFilter.C
+++ b/src/sampling/surface/distanceSurface/distanceSurfaceFilter.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,6 +28,7 @@ License
 #include "distanceSurface.H"
 #include "regionSplit.H"
 #include "syncTools.H"
+#include "ListOps.H"
 #include "vtkSurfaceWriter.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -306,16 +307,146 @@ void Foam::distanceSurface::filterKeepNearestRegions
 }
 
 
-void Foam::distanceSurface::filterByProximity()
+void Foam::distanceSurface::filterRegionProximity
+(
+    bitSet& ignoreCells
+) const
 {
     const searchableSurface& geom = geometryPtr_();
 
-    // Filtering for faces
+    // For face distances
     const pointField& fc = surface_.faceCentres();
 
-    bitSet faceSelection(fc.size());
-    label nTrimmed = 0;
+    // For region split
+    bitSet blockedFaces(filterPrepareRegionSplit(ignoreCells));
 
+    // Split region
+    regionSplit rs(mesh_, blockedFaces);
+    blockedFaces.clearStorage();
+
+    const labelList& regionColour = rs;
+
+    // For each face
+    scalarField faceDistance(fc.size(), GREAT);
+
+    {
+        List<pointIndexHit> nearest;
+        geom.findNearest
+        (
+            fc,
+            // Use initialized field (GREAT) to limit search too
+            faceDistance,
+            nearest
+        );
+        calcAbsoluteDistance(faceDistance, fc, nearest);
+    }
+
+    // Identical number of regions on all processors
+    scalarField areaRegion(rs.nRegions(), Zero);
+    scalarField distRegion(rs.nRegions(), Zero);
+
+    forAll(meshCells_, facei)
+    {
+        const label celli = meshCells_[facei];
+        const label regioni = regionColour[celli];
+
+        const scalar faceArea = surface_[facei].mag(surface_.points());
+        distRegion[regioni] += (faceDistance[facei] * faceArea);
+        areaRegion[regioni] += (faceArea);
+    }
+
+    Pstream::listCombineGather(distRegion, plusEqOp<scalar>());
+    Pstream::listCombineGather(areaRegion, plusEqOp<scalar>());
+
+    if (Pstream::master())
+    {
+        forAll(distRegion, regioni)
+        {
+            distRegion[regioni] /= (areaRegion[regioni] + VSMALL);
+        }
+    }
+
+    Pstream::listCombineScatter(distRegion);
+
+
+    // Define the per-face acceptance
+    // - region average distance is good
+    // Could also check against ignored cells (if selection refined later)
+
+    bitSet acceptFaces(fc.size());
+    bool prune(false);
+
+    forAll(meshCells_, facei)
+    {
+        const label celli = meshCells_[facei];
+        const label regioni = regionColour[celli];
+
+        if (distRegion[regioni] < absProximity_)
+        {
+            acceptFaces.set(facei);
+        }
+        else
+        {
+            prune = true;
+        }
+    }
+
+    // Heavier debugging
+    if (debug & 4)
+    {
+        const fileName outputName(surfaceName() + "-region-proximity-filter");
+
+        Info<< "Writing debug surface: " << outputName << nl;
+
+        surfaceWriters::vtkWriter writer
+        (
+            surface_.points(),
+            surface_,  // faces
+            outputName
+        );
+
+        writer.write("absolute-distance", faceDistance);
+
+        // Region segmentation
+        labelField faceRegion
+        (
+            ListOps::create<label>
+            (
+                meshCells_,
+                [&](const label celli){ return regionColour[celli]; }
+            )
+        );
+        writer.write("face-region", faceRegion);
+
+        // Region-wise filter state
+        labelField faceFilterState
+        (
+            ListOps::createWithValue<label>(surface_.size(), acceptFaces, 1, 0)
+        );
+        writer.write("filter-state", faceFilterState);
+    }
+
+
+    if (prune)
+    {
+        labelList pointMap, faceMap;
+        meshedSurface filtered
+        (
+            surface_.subsetMesh(acceptFaces, pointMap, faceMap)
+        );
+        surface_.transfer(filtered);
+
+        meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
+    }
+}
+
+
+void Foam::distanceSurface::filterFaceProximity()
+{
+    const searchableSurface& geom = geometryPtr_();
+
+    // For face distances
+    const pointField& fc = surface_.faceCentres();
 
     // For each face
     scalarField faceDistance(fc.size(), GREAT);
@@ -356,18 +487,19 @@ void Foam::distanceSurface::filterByProximity()
 
     // Using the absolute proximity of the face centres is more robust.
 
+    bitSet acceptFaces(fc.size());
+    bool prune(false);
 
     // Consider the absolute proximity of the face centres
     forAll(faceDistance, facei)
     {
         if (faceDistance[facei] <= absProximity_)
         {
-            faceSelection.set(facei);
+            acceptFaces.set(facei);
         }
         else
         {
-            ++nTrimmed;
-
+            prune = true;
             if (debug & 2)
             {
                 Pout<< "trim reject: "
@@ -380,14 +512,7 @@ void Foam::distanceSurface::filterByProximity()
     // Heavier debugging
     if (debug & 4)
     {
-        labelField faceFilterStatus(faceSelection.size(), Zero);
-
-        for (const label facei : faceSelection)
-        {
-            faceFilterStatus[facei] = 1;
-        }
-
-        const fileName outputName(surfaceName() + "-proximity-filter");
+        const fileName outputName(surfaceName() + "-face-proximity-filter");
 
         Info<< "Writing debug surface: " << outputName << nl;
 
@@ -400,16 +525,22 @@ void Foam::distanceSurface::filterByProximity()
 
         writer.write("absolute-distance", faceDistance);
         writer.write("normal-distance", faceNormalDistance);
-        writer.write("filter-state", faceFilterStatus);
+
+        // Region-wise filter state
+        labelField faceFilterState
+        (
+            ListOps::createWithValue<label>(surface_.size(), acceptFaces, 1, 0)
+        );
+        writer.write("filter-state", faceFilterState);
     }
 
 
-    if (returnReduce(nTrimmed, sumOp<label>()) != 0)
+    if (prune)
     {
         labelList pointMap, faceMap;
         meshedSurface filtered
         (
-            surface_.subsetMesh(faceSelection, pointMap, faceMap)
+            surface_.subsetMesh(acceptFaces, pointMap, faceMap)
         );
         surface_.transfer(filtered);