From 1460d7f1b1f1a66b90f3c35e9db52c957eef42b1 Mon Sep 17 00:00:00 2001
From: odebeir <odebeir@ulb.ac.be>
Date: Tue, 30 Oct 2012 16:55:05 +0100
Subject: [PATCH] compare bilateral

---
 doc/examples/plot_16bitbilateral.py    | 34 +++++++++-----
 doc/examples/plot_compare_bilateral.py | 63 ++++++++++++++++++++++++++
 2 files changed, 86 insertions(+), 11 deletions(-)
 create mode 100644 doc/examples/plot_compare_bilateral.py

diff --git a/doc/examples/plot_16bitbilateral.py b/doc/examples/plot_16bitbilateral.py
index 076d03c4..fc30aa6b 100644
--- a/doc/examples/plot_16bitbilateral.py
+++ b/doc/examples/plot_16bitbilateral.py
@@ -1,9 +1,23 @@
 """
 ==============================
-Simplified bilateral filtering
+Bilateral mean
 ==============================
+This example compares
 
-to complete
+* local mean
+* percentile mean
+* bilateral mean
+
+build on the local histogram distribution
+local mean uses all pixels belonging to the structuring element to compute average gray level,
+percentile mean uses only values between percentiles p0 and p1 (here 10% and 90%),
+whereas bilateral mean uses only pixels of the structuring element having a gray level situated inside
+g-s0 and g+s1 (here g-500 and g+500).
+The filters are applied on a 16 bit image (actual bitdepth is 12bit).
+
+Percentile and usual mean give here similar results, these filters smooth the complete image (background and details).
+Bilateral mean exhibits a high filtering rate for continuous area (i.e. background) while image higher frequencies
+remains untouched.
 
 """
 import numpy as np
@@ -13,23 +27,21 @@ from skimage import data
 from skimage.morphology import disk
 import skimage.filter.rank as rank
 
-a8 = (data.coins()).astype('uint8')
-
 a16 = (data.coins()).astype('uint16')*16
-selem = np.ones((20,20),dtype='uint8')
-f1 = rank.percentile_mean(a8,selem = selem,p0=.1,p1=.9)
+selem = disk(20)
+
+f1 = rank.percentile_mean(a16,selem = selem,p0=.1,p1=.9)
 f2 = rank.bilateral_mean(a16,selem = selem,s0=500,s1=500)
-selem = disk(50)
-f3 = rank.equalize(a16,selem = selem)
+f3 = rank.mean(a16,selem = selem)
 
 # display results
-fig, axes = plt.subplots(nrows=3, figsize=(15,15))
+fig, axes = plt.subplots(nrows=3, figsize=(15,10))
 ax0, ax1, ax2 = axes
 
-ax0.imshow(np.hstack((a8,f1)))
+ax0.imshow(np.hstack((a16,f1)))
 ax0.set_title('percentile mean')
 ax1.imshow(np.hstack((a16,f2)))
 ax1.set_title('bilateral mean')
 ax2.imshow(np.hstack((a16,f3)))
-ax2.set_title('local equalization')
+ax2.set_title('local mean')
 plt.show()
diff --git a/doc/examples/plot_compare_bilateral.py b/doc/examples/plot_compare_bilateral.py
new file mode 100644
index 00000000..977b7704
--- /dev/null
+++ b/doc/examples/plot_compare_bilateral.py
@@ -0,0 +1,63 @@
+"""
+====================================================
+Bilateral comparison
+====================================================
+
+In this example, we compare both bilateral implementation
+
+* filter.denoise_bilateral
+* filter.rank.bilateral_mean
+
+The first filter implements a spatial-gaussian and spectral-gaussian kernel bilateral filter whereas the latter implements
+a cylindrical kernel bilateral filter i.e. spatial-flat and spectral-flat kernel.
+
+The timing comparison is just for information since the kernel are not the same.
+
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from skimage import data
+from skimage.filter._denoise import denoise_bilateral
+from skimage.filter.rank import bilateral_mean
+from skimage.morphology import disk
+from skimage.filter import denoise_bilateral
+import time
+
+def exec_and_timeit(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
+        return (res,ms)
+    return wrapper
+
+
+@exec_and_timeit
+def den_bil(image):
+    return denoise_bilateral(a8,win_size=20,sigma_range=255,sigma_spatial=1)[:,:,0]*255
+
+@exec_and_timeit
+def rank_bil(image):
+    return bilateral_mean(a8.astype(np.uint16),disk(20),s0=10,s1=10)
+
+a8 = data.camera()
+selem = disk(10)
+
+f1,t1 = den_bil(a8)
+f2,t2 = rank_bil(a8)
+
+# display results
+fig, axes = plt.subplots(nrows=2, figsize=(15,10))
+ax0, ax1= axes
+
+ax0.imshow(np.hstack((f1,a8-f1)))
+ax0.set_title('denoise bilateral (%f ms)'%t1)
+ax1.imshow(np.hstack((f2,a8-f1)))
+ax1.set_title('bilateral mean (%f ms)'%t2)
+plt.show()