From 71c70ccd7129fe0409a308c2d2f0b85044b7de0f Mon Sep 17 00:00:00 2001
From: Ferdinand Deger <ferdinand.deger@gmail.com>
Date: Sat, 13 Jul 2013 09:52:09 +0200
Subject: [PATCH] Added isotropic TV denoising + comment

---
 skimage/filter/_denoise_cy.pyx | 28 ++++++++++++++++++++++++----
 1 file changed, 24 insertions(+), 4 deletions(-)

diff --git a/skimage/filter/_denoise_cy.pyx b/skimage/filter/_denoise_cy.pyx
index 6fd95f03..6e803062 100644
--- a/skimage/filter/_denoise_cy.pyx
+++ b/skimage/filter/_denoise_cy.pyx
@@ -200,7 +200,7 @@ def denoise_bilateral(image, Py_ssize_t win_size=5, sigma_range=None,
     return np.squeeze(out)
 
 
-def denoise_tv_bregman(image, double weight, int max_iter=100, double eps=1e-3):
+def denoise_tv_bregman(image, double weight, int max_iter=100, double eps=1e-3, anisotropic=False):
     """Perform total-variation denoising using split-Bregman optimization.
 
     Total-variation denoising (also know as total-variation regularization)
@@ -225,6 +225,8 @@ def denoise_tv_bregman(image, double weight, int max_iter=100, double eps=1e-3):
     max_iter: int, optional
         Maximal number of iterations used for the optimization.
 
+    anisotropic: boolean, optimal - switch between isotropic and anisotropic tv denoising
+
     Returns
     -------
     u : ndarray
@@ -239,6 +241,7 @@ def denoise_tv_bregman(image, double weight, int max_iter=100, double eps=1e-3):
     .. [3] Pascal Getreuer, "Rudin–Osher–Fatemi Total Variation Denoising
            using Split Bregman" in Image Processing On Line on 2012–05–19,
            http://www.ipol.im/pub/art/2012/g-tvd/article_lr.pdf
+       [4] http://www.math.ucsb.edu/~cgarcia/UGProjects/BregmanAlgorithms_JacquelineBush.pdf 
 
     """
 
@@ -325,9 +328,26 @@ def denoise_tv_bregman(image, double weight, int max_iter=100, double eps=1e-3):
                     bxx = bx[r, c, k]
                     byy = by[r, c, k]
 
-                    s = sqrt((ux + bxx)**2 + (uy + byy)**2)
-                    dxx = s * lam * (ux + bxx) / (s * lam + 1)
-                    dyy = s * lam * (uy + byy) / (s * lam + 1)
+                    # d_subproblem after reference [4]
+                    if anisotropic:
+                        s = ux + bxx
+                        if s > 1/lam:
+                            dxx = s - 1/lam
+                        elif s < -1/lam:
+                            dxx = s + 1/lam
+                        else:
+                            dxx = 0
+                        s = uy + byy
+                        if s > 1/lam:
+                            dyy = s - 1/lam
+                        elif s < -1/lam:
+                            dyy = s + 1/lam
+                        else:
+                            dyy = 0
+                    else:
+                        s = sqrt((ux + bxx)**2 + (uy + byy)**2)
+                        dxx = s * lam * (ux + bxx) / (s * lam + 1)
+                        dyy = s * lam * (uy + byy) / (s * lam + 1)
 
                     dx[r, c, k] = dxx
                     dy[r, c, k] = dyy