Repository: PMEAL/porespy

Issue:
enhance trim_nearby_peaks to accept a tolerance or threshold
At the moment it trims peaks that are closer to another peak than to solid, but it'd be nice if we could put some sort of absolute AND relative tolerance, so it only trims peaks if they are too close to a nearby peak within some acceptable amount.

Discussion:

Patch:
diff --git a/porespy/filters/_snows.py b/porespy/filters/_snows.py
index 346716570..aa17cb1a8 100644
--- a/porespy/filters/_snows.py
+++ b/porespy/filters/_snows.py
@@ -350,7 +350,7 @@ def trim_saddle_points(peaks, dt, max_iters=10):
     return peaks
 
 
-def trim_nearby_peaks(peaks, dt):
+def trim_nearby_peaks(peaks, dt, f=1.0):
     r"""
     Finds pairs of peaks that are nearer to each other than to the solid
     phase, and removes the peak that is closer to the solid.
@@ -363,6 +363,9 @@ def trim_nearby_peaks(peaks, dt):
     dt : ND-array
         The distance transform of the pore space for which the true peaks
         are sought.
+    f : scalar
+        Controls how close peaks must be before they are considered near to each
+        other.  Sets of peaks are tagged as near if ``d_neighbor < f * d_solid``.
 
     Returns
     -------
@@ -399,7 +402,7 @@ def trim_nearby_peaks(peaks, dt):
     dist_to_neighbor = temp[0][:, 1]
     del temp, tree  # Free-up memory
     dist_to_solid = dt[tuple(crds.T)]  # Get distance to solid for each peak
-    hits = np.where(dist_to_neighbor < dist_to_solid)[0]
+    hits = np.where(dist_to_neighbor < f * dist_to_solid)[0]
     # Drop peak that is closer to the solid than it's neighbor
     drop_peaks = []
     for peak in hits:
@@ -415,7 +418,6 @@ def trim_nearby_peaks(peaks, dt):
     return peaks > 0
 
 
-
 def _estimate_overlap(im, mode='dt', zoom=0.25):
     logger.trace('Calculating overlap thickness')
     if mode == 'watershed':


Test patch:
diff --git a/test/unit/test_filters.py b/test/unit/test_filters.py
index d3d44a028..9f3c2b50d 100644
--- a/test/unit/test_filters.py
+++ b/test/unit/test_filters.py
@@ -5,6 +5,7 @@ import porespy as ps
 import scipy.ndimage as spim
 from skimage.morphology import disk, ball, skeletonize_3d
 from skimage.util import random_noise
+from scipy.stats import norm
 
 
 class FilterTest():
@@ -476,6 +477,21 @@ class FilterTest():
         p2 = (im[0, ...]).sum()
         np.testing.assert_approx_equal(np.around(p1 / p2, decimals=1), 1)
 
+    def test_trim_nearby_peaks_threshold(self):
+        np.random.seed(10)
+        dist = norm(loc=7, scale=5)
+        im = ps.generators.polydisperse_spheres([100, 100, 100],
+                                                porosity=0.8, dist=dist)
+        im_dt = edt(im)
+        im_dt = im_dt
+        dt = spim.gaussian_filter(input=im_dt, sigma=0.4)
+        peaks = ps.filters.find_peaks(dt=dt)
+        peaks_far = ps.filters.trim_nearby_peaks(peaks=peaks, dt=dt)
+        peaks_close = ps.filters.trim_nearby_peaks(peaks=peaks, dt=dt, f=0.3)
+        num_peaks_after_far_trim = spim.label(peaks_far)[1]
+        num_peaks_after_close_trim = spim.label(peaks_close)[1]
+        assert num_peaks_after_far_trim <= num_peaks_after_close_trim
+
 
 if __name__ == '__main__':
     t = FilterTest()