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