From 35ebb64c3fa5f119bbce6534d72bd9d216020985 Mon Sep 17 00:00:00 2001
From: Olivier Debeir <odebeir@ulb.ac.be>
Date: Tue, 16 Oct 2012 16:20:12 +0200
Subject: [PATCH] compare ndimage.percentile in demo

---
 skimage/rank/local/demo_benchmark.py | 75 ++++++++++++++++------------
 1 file changed, 44 insertions(+), 31 deletions(-)

diff --git a/skimage/rank/local/demo_benchmark.py b/skimage/rank/local/demo_benchmark.py
index feae3a8b..d3c3d00a 100644
--- a/skimage/rank/local/demo_benchmark.py
+++ b/skimage/rank/local/demo_benchmark.py
@@ -1,44 +1,53 @@
 import numpy as np
 import matplotlib.pyplot as plt
+import time
 
 from skimage import data
-from skimage.morphology import dilation
-import skimage.rank as rank
+from skimage.morphology import dilation,disk
 from skimage.filter import median_filter
+from scipy.ndimage.filters import percentile_filter
+import skimage.rank as rank
 
-from skimage.rank.local.tools import log_timing
+def log_timing(func):
+    """ Decorator that returns both function results and execution time
+    (result, ms)
+    """
+    def wrapper(*arg):
+        t1 = time.time()
+        res = func(*arg)
+        t2 = time.time()
+        ms = (t2-t1)*1000.0
+        print '%s took %0.3f ms' % (func.func_name, ms)
+        return (res,ms)
+    return wrapper
 
-@log_timing
-def cr_max(image,selem):
-    return rank.maximum(image=image,selem = selem)
 
 @log_timing
 def cr_med(image,selem):
     return rank.median(image=image,selem = selem)
 
-@log_timing
-def cm_dil(image,selem):
-    return dilation(image=image,selem = selem)
-
 @log_timing
 def ctmf_med(image,radius):
     return median_filter(image=image,radius=radius)
 
+@log_timing
+def ndi_med(image,n):
+    return percentile_filter(image,50,size=n*2-1)
 
 def compare_dilate():
-    """comparison between
+    """ Comparison between
     - crank.maximum rankfilter implementation
     - cmorph.dilate cython implementation
 
     on increasing structuring element size and increasing image size
     """
-#    a = (np.random.random((500,500))*256).astype('uint8')
     a = data.camera()
 
     rec = []
     e_range = range(1,20,1)
     for r in e_range:
-        elem = np.ones((r,r),dtype='uint8')
+#        elem = np.ones((r,r),dtype='uint8')
+        elem = disk(r+1)
         #        elem = (np.random.random((r,r))>.5).astype('uint8')
         rc,ms_rc = cr_max(a,elem)
         rcm,ms_rcm = cm_dil(a,elem)
@@ -56,7 +65,8 @@ def compare_dilate():
     plt.imshow(np.hstack((rc,rcm)))
 
     r = 9
-    elem = np.ones((r,r),dtype='uint8')
+#    elem = np.ones((r,r),dtype='uint8')
+    elem = disk(r+1)
 
     rec = []
     s_range = range(100,1000,100)
@@ -79,7 +89,7 @@ def compare_dilate():
     plt.show()
 
 def compare_median():
-    """comparison between
+    """ Comparison between
     - crank.median rankfilter implementation
     - ctmf.median_filter filter
 
@@ -88,47 +98,50 @@ def compare_median():
     a = data.camera()
 
     rec = []
-    e_range = range(2,40,4)
+    e_range = range(2,30,4)
     for r in e_range:
-        elem = np.ones((2*r+1,2*r+1),dtype='uint8')
-        #        elem = (np.random.random((r,r))>.5).astype('uint8')
+        elem = disk(r+1)
         rc,ms_rc = cr_med(a,elem)
         rctmf,ms_rctmf = ctmf_med(a,r)
-        rec.append((ms_rc,ms_rctmf))
+        rndi,ms_ndi = ndi_med(a,r)
+        rec.append((ms_rc,ms_rctmf,ms_ndi))
         # check if results are identical
-#        assert  (rc==rctmf).all()
+        # obviously they cannot be identical since structuring element are different (octagon<>disk)
+        # assert  (rc==rctmf).all()
 
     rec = np.asarray(rec)
 
     plt.figure()
     plt.title('increasing element size')
     plt.plot(e_range,rec)
-    plt.legend(['rank.median','ctmf.median_filter'])
-    plt.figure()
-    plt.imshow(np.hstack((rc,rctmf)))
+    plt.legend(['rank.median','ctmf.median_filter','ndimage.percentile'])
     plt.ylabel('time (ms)')
     plt.xlabel('element radius')
+    plt.figure()
+    plt.imshow(np.hstack((rc,rctmf,rndi)))
+    plt.xlabel('rank.median vs ctmf.median_filter vs ndimage.percentile')
 
     r = 9
-    elem = np.ones((r*2+1,r*2+1),dtype='uint8')
+    elem = disk(r+1)
 
     rec = []
-    s_range = range(100,1000,100)
+    s_range = [100,200,500,1000,2000]
     for s in s_range:
         a = (np.random.random((s,s))*256).astype('uint8')
-        (rc,ms_rc) = cr_max(a,elem)
+        (rc,ms_rc) = cr_med(a,elem)
         rctmf,ms_rctmf = ctmf_med(a,r)
-        rec.append((ms_rc,ms_rctmf))
-#        assert  (rc==rcm).all()
+        rndi,ms_ndi = ndi_med(a,r)
+        rec.append((ms_rc,ms_rctmf,ms_ndi))
+        # check if results are identical
+        # obviously they cannot be identical since structuring element are different (octagon<>disk)
+        # assert  (rc==rctmf).all()
 
     rec = np.asarray(rec)
 
     plt.figure()
     plt.title('increasing image size')
     plt.plot(s_range,rec)
-    plt.legend(['rank.median','ctmf.median_filter'])
-    plt.figure()
-    plt.imshow(np.hstack((rc,rctmf)))
+    plt.legend(['rank.median','ctmf.median_filter','ndimage.percentile'])
     plt.ylabel('time (ms)')
     plt.xlabel('image size')