thearn · April 2, 2018 14:21 · Apr 2, 2018 · Apr 2, 2018 · Apr 2, 2018
diff --git a/rkmd.py b/rkmd.py
@@ -1,40 +1,52 @@
 import numpy as np
 import matplotlib.pyplot as plt
-from scipy.optimize import fmin
+from scipy.optimize import fmin, fminbound, fmin_bfgs
 from scipy.misc import ascent
-from skimage.restoration import denoise_tv_chambolle
+import pywt
 from scipy.stats import kurtosis
 true_image = ascent()
 
-# add some Gaussian noise 
-sigma = 22.3
+sigma = 25.0
 test_image = true_image + sigma * np.random.randn(*true_image.shape)
 
-# explore the kurtosis of the residual under different denoising parameter levels
-ks = []
-w_sample = np.linspace(0.001, 45, 500)
-for weight in w_sample:
-
-    # keep it simple: use tv denoising
-    restored_image = denoise_tv_chambolle(test_image, weight)
-    residual = test_image - restored_image
-    k = np.abs(kurtosis(residual.flatten()))
-    ks.append(k)
-
-# minimum seen on plot, use numerical optimization to locate it
-optimal_weight = 23.1
-restored_image = denoise_tv_chambolle(test_image, optimal_weight)
-
-# plot residual kurtosis
-plt.figure()
-plt.plot(w_sample, ks)
-
-# plot noisy vs. restored image
-plt.figure()
+wavelet = 'haar'
+levels = 4
+coeffs = pywt.wavedec2(test_image, wavelet, level=levels)
+
+lowest_kurtosis = 1e99
+for i in range(1, levels + 1)[::-1]:
+    hn, vn, dn = coeffs[i]
+
+    def filter(t):
+        filtered_coeffs = [pywt.threshold(hn, t),
+                           pywt.threshold(vn, t),
+                           pywt.threshold(dn, t)]
+        return filtered_coeffs
+
+    def restore(t):
+        coeffs[i] = filter(t)
+        return pywt.waverec2(coeffs, wavelet)
+
+    def objective(t):
+        restored = restore(np.abs(t))
+        k = kurtosis((test_image - restored).flatten())
+
+        return abs(k) 
+
+    #t_opt = fminbound(objective, 0.0, 5e2, xtol=1e-08)
+    t_opt = fmin_bfgs(objective, 10.0, disp=0)
+    this_k = objective(t_opt)
+    print(i, t_opt, this_k)
+    if t_opt > lowest_kurtosis:
+        break
+    filtered_coeffs = filter(t_opt)
+    coeffs[i] = filtered_coeffs
+
+final_restored = pywt.waverec2(coeffs, wavelet)
+
 plt.gray()
 plt.subplot(121)
 plt.imshow(test_image)
 plt.subplot(122)
-plt.imshow(restored_image)
-plt.show()
-
+plt.imshow(final_restored)
+plt.show()
diff --git a/rkmd.py b/rkmd.py
@@ -6,24 +6,30 @@
 from scipy.stats import kurtosis
 true_image = ascent()
 
+# add some Gaussian noise 
 sigma = 22.3
 test_image = true_image + sigma * np.random.randn(*true_image.shape)
 
+# explore the kurtosis of the residual under different denoising parameter levels
 ks = []
 w_sample = np.linspace(0.001, 45, 500)
 for weight in w_sample:
-
+
+    # keep it simple: use tv denoising
     restored_image = denoise_tv_chambolle(test_image, weight)
     residual = test_image - restored_image
     k = np.abs(kurtosis(residual.flatten()))
     ks.append(k)
 
-optimal_weight = 22.3
+# minimum seen on plot, use numerical optimization to locate it
+optimal_weight = 23.1
 restored_image = denoise_tv_chambolle(test_image, optimal_weight)
 
+# plot residual kurtosis
 plt.figure()
 plt.plot(w_sample, ks)
 
+# plot noisy vs. restored image
 plt.figure()
 plt.gray()
 plt.subplot(121)

diff --git a/rkmd.py b/rkmd.py
@@ -0,0 +1,34 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.optimize import fmin
+from scipy.misc import ascent
+from skimage.restoration import denoise_tv_chambolle
+from scipy.stats import kurtosis
+true_image = ascent()
+
+sigma = 22.3
+test_image = true_image + sigma * np.random.randn(*true_image.shape)
+
+ks = []
+w_sample = np.linspace(0.001, 45, 500)
+for weight in w_sample:
+
+    restored_image = denoise_tv_chambolle(test_image, weight)
+    residual = test_image - restored_image
+    k = np.abs(kurtosis(residual.flatten()))
+    ks.append(k)
+
+optimal_weight = 22.3
+restored_image = denoise_tv_chambolle(test_image, optimal_weight)
+
+plt.figure()
+plt.plot(w_sample, ks)
+
+plt.figure()
+plt.gray()
+plt.subplot(121)
+plt.imshow(test_image)
+plt.subplot(122)
+plt.imshow(restored_image)
+plt.show()
+