Total Variation Denoising

Total variation denoising removes noise from an image by solving an optimization problem: among all images, it looks for the one that stays close to the noisy observation while having the smallest total variation, the summed magnitude of its gradient. The model is due to Rudin, Osher and Fatemi (1992) and is known as the ROF model; Getreuer (2012) calls the minimization the Rudin–Osher–Fatemi or ROF problem. Because total variation charges a jump only by its height, not by how abruptly it happens, the result keeps sharp edges while flattening noise, and denoised images take on a cartoon-like, piecewise-flat look. Image denoising places the method among the other families of denoisers; this article covers the model, its exact solutions, the algorithms that solve it, and the two side effects Getreuer describes, loss of contrast and staircasing.

Total variation is also a building block well beyond denoising: it regularizes deblurring, the TV-L1 optical flow of Zach, Pock and Bischof, and compressed-sensing MRI reconstruction. The optimization article explains why convexity makes such problems reliable to solve.

Problem

An observed image $f$ is a clean image corrupted by additive noise, $f = u + \eta$, with $\eta$ white and Gaussian in the standard setting (image noise describes real camera noise). In our reading, linear smoothing cannot separate noise from edges, because both contain high frequencies; the Wiener filter, the best linear estimate for a Gaussian image model, is reported there not to capture edges. Total variation denoising instead uses a nonlinear prior: clean images are assumed to consist of regions of slowly varying intensity separated by sharp boundaries, and the penalty is chosen so that a boundary costs no more than its height times its length.

Inputs and Outputs

• Input: a noisy grayscale image $f$ (colour images are handled channel by channel or with a vectorial total variation, described under Variants), and a regularization weight that sets the balance between smoothing and fidelity.
• Output: the denoised image $u$, the unique minimizer of the ROF energy. The method is non-blind: the weight encodes the noise level, and choosing it is part of using the method.

Intuition

Consider a single row of pixels. Its total variation is $\sum_n |u_{n+1} - u_n|$. For a monotone row it equals the difference between the last and the first value, whatever happens in between: a sharp step from 0 to 1 and a gentle ramp from 0 to 1 cost the same. Noise adds many small ups and downs, and each one adds to the sum. Minimizing the total variation therefore removes oscillations but has no reason to blur an edge, unlike a quadratic penalty $\sum_n (u_{n+1} - u_n)^2$, which charges a step of height 1 a full 1 but the same rise spread over ten pixels only $10 \times 0.1^2 = 0.1$, and so prefers to smear it. Boyd and Vandenberghe’s comparison of the two penalties, summarized in optimization, makes the same point.

The fidelity term pulls the result back toward the data. Where the noise is small compared with the weight, the result goes flat; where a jump is large and the region on either side is wide, the jump survives. The same trade-off explains both side effects. A small bright region has a short boundary but also a small area to defend its brightness with, so lowering its contrast is cheap: it loses contrast in proportion to its perimeter over its area. And on a noisy ramp, flat segments joined by small jumps cost little more total variation than the ramp itself while absorbing the noise, so the ramp breaks into stairs.

Algorithm

The ROF problem is convex but not differentiable where the gradient vanishes, which is exactly where its solutions spend most of their area. Getreuer (2012, §2) surveys the algorithms: regularizing total variation with a small $\epsilon$, evolving the gradient-descent PDE, solving for its steady state directly, duality-based methods, Newton-based methods, graph cuts, frame shrinkage, and operator splitting, of which split Bregman is his subject. Three are worth knowing.

1. Gradient-descent PDE. Rudin, Osher and Fatemi discretized the gradient-descent flow $\partial_t u = \operatorname{div}(\nabla u / |\nabla u|) + \lambda_G (f - u)$ with Neumann boundary conditions, as Getreuer describes (his Eq. 17 and Algorithm 1, with $0/0$ taken as 0, and a small positive time step). Because the problem is convex, the steady state of the flow is the minimizer.

2. Chambolle’s dual algorithm. It works on a dual vector field $p$, one component per image axis, and recovers the image from it. scikit-image’s denoise_tv_chambolle implements it (scikit-image 0.25.2 source); in the library’s scaling, each iteration:

   • computes $u = f - \operatorname{div} p$ and its forward-difference gradient $g = \nabla u$;
   • updates $p \leftarrow (p - \tau g) / (1 + (\tau/\lambda)|g|)$, with step $\tau = 1/(2d)$ for a $d$-dimensional image;
   • stops after max_num_iter iterations (default 200) or when a cost changes by less than eps (default $2 \times 10^{-4}$) times its first value.

3. Split Bregman. Goldstein and Osher’s method, as Getreuer (2012, §5) presents it, introduces an auxiliary field $d = \nabla u$ and enforces the constraint by Bregman iteration, alternating three steps:

   • a d-step, a vectorial shrinkage $d = \frac{\nabla u + b}{|\nabla u + b|} \max(|\nabla u + b| - 1/\gamma, 0)$;
   • a u-step, the screened Poisson equation $\lambda_G u - \gamma \Delta u = \lambda_G f - \gamma \operatorname{div}(d - b)$;
   • a Bregman update $b \leftarrow b + \nabla u - d$.

   Getreuer solves the u-step approximately with one Gauss–Seidel sweep per iteration, following Goldstein and Osher, notes that it may also be solved in the Fourier or DCT domain, and uses $\gamma = 5$, finding the algorithm fairly insensitive to its value. The example below solves the u-step exactly with the discrete cosine transform, which diagonalizes the Neumann Laplacian.

Getreuer notes that split Bregman with this alternating scheme is equivalent to the alternating direction method of multipliers (ADMM), after Esser, and that Bregman iteration with linear constraints is equivalent to the augmented Lagrangian method.

Mathematical Formulation

The ROF energy

Getreuer (2012, Eq. 1) writes the ROF problem as

$$\hat u = \arg\min_{u \in BV(\Omega)} \; \lVert u \rVert_{TV(\Omega)} + \frac{\lambda_G}{2} \int_\Omega \big(f(\mathbf{x}) - u(\mathbf{x})\big)^2 \, d\mathbf{x},$$

where $BV(\Omega)$ is the space of functions of bounded variation and, for smooth $u$, $\lVert u \rVert_{TV} = \int_\Omega |\nabla u| \, d\mathbf{x}$. The general definition, a supremum of $\int u \operatorname{div} \mathbf{g}$ over vector fields with $|\mathbf{g}| \le 1$ (Getreuer’s Eq. 3), gives a step edge the total variation height times length, which is what lets the minimizer jump. Getreuer states that for $f \in L^2$ the minimizer exists, is unique and is stable with respect to perturbations of $f$.

This article and scikit-image’s denoise_tv_chambolle use the equivalent scaling with the weight on the other side,

$$E(u) = \sum_{\text{pixels}} |\nabla u| + \frac{1}{2\lambda} \sum_{\text{pixels}} (u - f)^2, \qquad \lambda = 1/\lambda_G ,$$

so that a larger $\lambda$ means more smoothing. The library’s docstring states the ROF problem with Getreuer’s $\lambda$ and says that weight equals $1/\lambda$; the update in the source confirms it: in our derivation, a fixed point of the update has $p = -\lambda \nabla u / |\nabla u|$ wherever the gradient is nonzero, so with $u = f - \operatorname{div} p$ it satisfies the optimality condition below for exactly this energy. Conventions differ across sources: image denoising writes the unconstrained form with $\lambda$ multiplying $\int (u - z)^2$ without the factor $\tfrac12$, so its $\lambda$ is $1/(2\lambda)$ in the notation here, and scikit-image’s split Bregman function denoise_tv_bregman documents yet another scaling of its weight. Translate the weight before comparing settings across tools.

Rudin, Osher and Fatemi originally posed a constrained problem, minimizing total variation subject to the residual having zero mean and the noise variance, as image denoising summarizes. Getreuer (2012, §6) states the link: minimizing total variation subject to $\int (f - u)^2 = \sigma^2 |\Omega|$ is equivalent to the unconstrained problem for one value of $\lambda_G$, the Lagrange multiplier of the constraint, but no closed form gives that value from $\sigma$.

Bayesian view. Getreuer (2012, Eqs. 5–8) derives the model as a maximum a posteriori estimate: with Gaussian noise of variance $\sigma^2$ and a prior $-\log p(u) = \mu \lVert u \rVert_{TV}$, the MAP estimate minimizes $\mu \lVert u \rVert_{TV} + \frac{1}{2\sigma^2} \int (f - u)^2$, the ROF problem with $1/\lambda_G = \mu \sigma^2$. Other noise models change the fidelity term (Variants below).

Discretization

The discrete energy needs a discrete gradient. Getreuer (2012, Eqs. 12 and 29–40) uses forward differences, $\partial u_n = u_{n+1} - u_n$, set to zero at the last sample, which is the forward difference of the half-sample symmetric extension; the divergence is minus the adjoint, a backward difference, and the Laplacian $\operatorname{div} \nabla$ is the 5-point Laplacian in the interior. The isotropic total variation sums $\sqrt{(\partial_x u)^2 + (\partial_y u)^2}$ over pixels. scikit-image’s Chambolle solver uses the same differences (scikit-image 0.25.2 source), so the two algorithms in the example minimize the same discrete energy, which, being strictly convex, has one minimizer. Getreuer notes that central differences are unsuitable because they ignore one-sample-wide structures, and that a discretization should also be consistent, its total variation converging to the continuous one as the grid is refined, citing Wang and Lucier.

Optimality

Because the energy is convex, $u$ is the minimizer exactly when $0$ lies in its subdifferential. Written with a vector field $\mathbf{z}$ that has $|\mathbf{z}| \le 1$ everywhere and equals $\nabla u / |\nabla u|$ wherever the gradient is nonzero, the condition is

$$u - f = \lambda \operatorname{div} \mathbf{z}, \qquad \mathbf{z} \cdot \boldsymbol{\nu} = 0 \text{ on } \partial\Omega .$$

Exhibiting such a field proves that a candidate is the exact minimizer. The dual variable $p$ of scikit-image’s Chambolle solver is this field up to sign and the factor $\lambda$, which is why the algorithm returns the image as $f$ minus a divergence. Two consequences follow directly. Summing the discrete divergence over the image gives zero, so the mean is preserved: $\operatorname{mean}(u) = \operatorname{mean}(f)$. And where $u$ is flat, $\mathbf{z}$ is only constrained by $|\mathbf{z}| \le 1$, which leaves room for whole regions to be exactly constant.

The disc: an exact solution

Let $f$ equal $h$ inside a disc of radius $R$ and 0 outside. On the whole plane, Meyer showed, as Getreuer (2012, Eqs. 65–66) reports, that the solution keeps the edge exactly and lowers the disc:

$$u = \max\!\left(0, \; h - \frac{2\lambda}{R}\right) \text{ inside}, \qquad u = 0 \text{ outside}$$

(Getreuer writes the drop as $2/(\lambda_G r)$). Getreuer adds that Strong and Chan made a similar analysis assuming the edges do not move, and showed similar behaviour on compact domains and for shapes other than discs.

The result can be checked with the optimality condition. For $u = a \mathbf{1}_{\text{disc}}$, $\lambda$ times the energy is $\lambda \, a \, 2\pi R + \tfrac12 (h - a)^2 \pi R^2$, minimized at $a = h - 2\lambda/R$. The field $\mathbf{z} = -\mathbf{x}/R$ inside the disc and $\mathbf{z} = -R\,\mathbf{x}/|\mathbf{x}|^2$ outside has $|\mathbf{z}| \le 1$, divergence $-2/R$ inside and $0$ outside, and equals minus the outward normal on the circle, so it certifies the solution. The drop $2\lambda/R$ is $\lambda$ times the disc’s perimeter divided by its area: it does not depend on the height $h$, it grows as the disc shrinks, and a disc with $R \le 2\lambda/h$ disappears entirely. Total variation denoising is in this sense scale-dependent. The figure in Getreuer (2012, p. 91) shows three discs of equal radius and different heights dropping by almost the same amount (0.23 and 0.24), a larger disc dropping less (0.13) and a small one more (0.68).

A bounded image (our derivation). A real image has edges, and with the Neumann boundary condition the mean is preserved, so the outside cannot stay at 0. In our derivation, if the solution keeps the two-level form, the inside field is unchanged and the inside level is still $a = h - 2\lambda/R$, while preserving the mean lifts the outside to

$$b = (h - a) \frac{A}{|\Omega| - A},$$

with $A$ the disc’s area and $|\Omega|$ the image’s. For a disc-shaped image of radius $L$, the radial field $z(r) = [-R + R(r^2 - R^2)/(L^2 - R^2)]/r$ outside the disc satisfies every condition, so the form is exact there; for a square image we have not proved that such a field exists, and the example tests the prediction instead. On the disc-shaped image, the same radial construction shows that once $a \le b$ the solution is the constant $\operatorname{mean}(f)$; the example tests this on the square too.

The ramp: an exact discrete solution

Let a row be the ramp $f_n = s n$, $n = 0, \ldots, M - 1$. In our derivation, the exact solution of the discrete energy keeps every interior sample unchanged and flattens the $k + 1$ samples at each end, to

$$c = \frac{\lambda}{k + 1} + \frac{s k}{2} \text{ at the start}, \qquad 1 - c \text{ at the end (for } s(M-1) = 1\text{)},$$

where $k$ is the largest integer with $s k (k+1)/2 \le \lambda$. The discrete optimality condition $u_n - f_n = \lambda (z_n - z_{n-1})$, with $z_{-1} = z_{M-1} = 0$ and $z_n$ equal to the sign of $u_{n+1} - u_n$ wherever that difference is nonzero, is met by $z_n = (n + 1)(c - s n / 2)/\lambda$ on the flat end, which rises to exactly 1 at $n = k$ and stays in $[0, 1]$, and by $z_n = 1$ in the interior, where it gives $u_n = f_n$. In the continuum limit the flat end is $\sqrt{2\lambda/s}$ long. A 2D image whose rows are all the same ramp has the same solution in every row, because averaging a candidate over the rows cannot raise the energy and the minimizer is unique. So a clean ramp is not staircased: total variation keeps its slope, apart from its two ends.

Parameters

• The weight ($\lambda$ here, weight in denoise_tv_chambolle, $1/\lambda_G$ in Getreuer) is the only model parameter. It sets the smallest feature that survives. By analogy with the disc solution (our extrapolation; Getreuer reports that Strong and Chan found similar behaviour for other shapes), a feature with perimeter-to-area ratio $P/A$ loses roughly $\lambda P/A$ of its contrast. Getreuer’s examples show the trade-off on a photograph: too strong a weight gives a cartoon, too weak leaves noise, and PSNR peaks in between.
• Choosing it from the noise level. Getreuer (2012, §6) uses the discrepancy principle, choosing the weight so that the residual’s variance matches the noise variance, and solves for it with an iteration he attributes to Chambolle (his Algorithm 6), $\lambda_G \leftarrow \lambda_G \lVert f - u \rVert / \sigma$, started from an empirical fit and run five times, which he finds sufficient for $\lVert f - u \rVert$ to be usually within 10% of $\sigma$. He notes that the principle tends to over-smooth slightly compared with the weight of least mean squared error.
• Solver settings. For Chambolle’s algorithm, max_num_iter and eps; for split Bregman, $\gamma$ and the tolerance. The cost scikit-image uses to stop is, in our reading of the source, not the ROF energy itself: it is computed as E = (d**2).sum() followed by E += weight * norm.sum() and divided by the number of pixels, with d the residual $u - f$ and norm the gradient magnitudes, so it is $\lVert u - f \rVert^2 + \lambda \, TV(u)$, whereas $2\lambda E(u) = \lVert u - f \rVert^2 + 2\lambda \, TV(u)$; it serves only as a stopping rule. The example shows that the defaults can stop far from the minimizer when the weight is large.

Complexity

Every algorithm above costs a few passes over the image per iteration, $O(n)$ for $n$ pixels, plus $O(n \log n)$ for a transform-domain u-step in split Bregman. The number of iterations decides the running time, and it grows with the weight: in the example, Chambolle’s algorithm with 20,000 iterations had not reached split Bregman’s energy at the disc weights (1 to 8) and missed the exact ramp at $\lambda = 2$ while reaching it at $\lambda = 0.5$. Memory is a few image-sized arrays (the dual field has one array per dimension).

Implementation

The example tests the exact solutions above against two solvers of the same discrete energy: a split Bregman iteration written in NumPy, following Getreuer’s Algorithm 5 with the u-step solved exactly through SciPy’s DCT ($\gamma = 32$, 3,000 iterations, started from $f$), and scikit-image’s denoise_tv_chambolle with eps=0 and max_num_iter=20000, so that it runs exactly 20,000 iterations; the disc part also runs denoise_tv_chambolle with its defaults. Arrays are indexed [row, column], grey levels lie in $[0, 1]$, the discs are centred at $((N-1)/2, (N-1)/2)$, and the boundaries are Neumann, as in the formulation above. w in the code is $\lambda$.

• Part A, the disc. A $256 \times 256$ image with an antialiased disc (8 × 8 supersampling per pixel) of radius $R \in \{8, 16, 32\}$, height $h = 1$, and weight chosen so that $2\lambda/R$ is 0.25 or 0.5; the same at $h = 0.6$ for $R = 16$; and $R = 16$ with $\lambda = 10$ ($2\lambda/R = 1.25$), where the disc should vanish. The inside level is the mean over $r \le R/2$, the outside level the mean over $r \ge R + 8$; they are compared with $a$ and $b$ from the formulas above, with $A$ the summed antialiased disc.
• Part B, the ramp image. $64 \times 128$, each row running from 0 to 1 in steps of 1/127, at $\lambda = 0.5$ and 2.0, against the exact discrete solution.
• Part C, staircasing. A 1D ramp of 256 samples from 0 to 1, with Gaussian noise of standard deviation 0, 0.01 or 0.03 (seeds 0–4), denoised at $\lambda = 0.05$. On samples 32–223, away from the flat ends, a step counts as flat if $|u_{n+1} - u_n|$ is below a tenth of the ramp’s slope, and a plateau is three or more samples joined by flat steps. The code prints, as ranges over seeds, the share of flat steps in the input and the output, the number of plateaus, their mean length, the RMS error against the clean ramp, and the largest difference between the two solvers.

The checks and predictions were written down, with their tolerances, before the example was first run, and the values of $R$, the ratios $2\lambda/R$, the heights, the weights and the noise levels were chosen at planning to cover the regimes of the derivations (a moderate and a large drop, a vanishing disc, a clean and two noisy ramps), not tuned on results; the ratios were set high on purpose so that the predicted drop is large. The disc part is therefore an illustration of the formulas and a test of how far the discrete problem is from the continuous one, not an independent test of Meyer’s result. The checks were: the mean is preserved to $10^{-12}$; the split Bregman solver reproduces the exact ramp to $10^{-8}$ at both weights, and Chambolle’s at $\lambda = 0.5$; split Bregman’s vanishing disc is constant and equal to $\operatorname{mean}(f)$ to $10^{-5}$; the two solvers agree on the 1D signals to $3.9 \times 10^{-5}$, a tenth of the flat-step threshold, and on every disc’s inside level to $10^{-3}$; and the data have the planned properties. The predictions were:

• H1: the gap between the discrete inside level and $a$ shrinks strictly as $R$ grows from 8 to 16 to 32, at both ratios;
• H2: the outside level is within 10% of $b$ in every disc;
• H3: at $R = 16$ the drop $h -$ inside differs between $h = 0.6$ and $h = 1$ by at most 0.005, at both ratios;
• H5: with the library’s defaults the inside level differs from split Bregman’s by more than 0.05 in every disc;
• H6: at these large weights Chambolle’s 20,000 iterations leave a higher energy than split Bregman’s in every disc, a vanishing disc that is not flat to $10^{-3}$, and a ramp error above $10^{-6}$ at $\lambda = 2$;
• H8: with noise, at least 40% of the output’s interior steps are flat on every seed while at most 5% of the input’s are, and without noise none are;
• H9: the largest output error over seeds is below the smallest input error, at both noise levels;
• H10: the shortest mean plateau at $\sigma = 0.03$ is longer than the longest at $\sigma = 0.01$.

The vanishing disc (H4) and the kept ramp (H7) are decided by the checks. The example uses NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2 on one thread, runs in about two minutes, and printed identical output on two runs.

import os
os.environ["OPENBLAS_NUM_THREADS"] = "1"          # one BLAS thread, before NumPy
os.environ["OMP_NUM_THREADS"] = "1"
import numpy as np
from scipy.fft import dctn, idctn
from skimage.restoration import denoise_tv_chambolle

# Discrete ROF energy, the one both solvers minimize (weight w = lambda):
#   E(u) = sum |grad u| + sum (u - f)^2 / (2 w),  forward differences, zero at the last index
def grad(u):
    g = np.zeros((u.ndim,) + u.shape)
    for ax in range(u.ndim):
        sl = [slice(None)] * u.ndim
        sl[ax] = slice(0, -1)
        g[(ax,) + tuple(sl)] = np.diff(u, axis=ax)
    return g

def div(v):                                        # minus the adjoint of grad
    out = np.zeros(v.shape[1:])
    for ax in range(v.ndim - 1):
        p = np.moveaxis(v[ax], ax, 0)
        d = p.copy()
        d[1:] -= p[:-1]
        d[-1] = -p[-2]
        out += np.moveaxis(d, 0, ax)
    return out

def energy(u, f, w):
    return np.sqrt((grad(u) ** 2).sum(0)).sum() + ((u - f) ** 2).sum() / (2 * w)

def split_bregman(f, w, gamma=32.0, iters=3000):
    """Getreuer's (2012) split Bregman iteration, u-step solved exactly with the DCT."""
    lam = 1.0 / w
    mu = 0.0
    for ax, n in enumerate(f.shape):               # eigenvalues of -Laplacian (Neumann)
        shape = [1] * f.ndim
        shape[ax] = n
        mu = mu + (2 - 2 * np.cos(np.pi * np.arange(n) / n)).reshape(shape)
    u = f.copy()
    d = np.zeros((f.ndim,) + f.shape)
    b = np.zeros_like(d)
    for _ in range(iters):
        s = grad(u) + b                            # d-step: vectorial shrinkage by 1/gamma
        mag = np.sqrt((s ** 2).sum(0))
        d = s * (np.maximum(mag - 1 / gamma, 0) / np.where(mag > 0, mag, 1))
        rhs = lam * f - gamma * div(d - b)         # u-step: (lam - gamma*Laplacian) u = rhs
        u = idctn(dctn(rhs, norm="ortho") / (lam + gamma * mu), norm="ortho")
        b = b + grad(u) - d                        # Bregman update
    return u

def chambolle(f, w, iters=20000):                  # eps=0: run exactly `iters` iterations
    return denoise_tv_chambolle(f, weight=w, eps=0, max_num_iter=iters)

# ---------------- A. Disc: the exact solution against both solvers ----------------
N = 256
c = (N - 1) / 2
yy, xx = np.mgrid[:N, :N] - c
rr = np.hypot(yy, xx)

def disc(R, h, ss=8):                              # antialiased disc, 8x8 supersampling
    t = (np.arange(N * ss) + 0.5) / ss - 0.5 - c
    Y, X = np.meshgrid(t, t, indexing="ij")
    return h * (X**2 + Y**2 < R * R).reshape(N, ss, N, ss).mean(axis=(1, 3))

print("A. Disc of radius R and height h on a 256 x 256 image, weight w = (2w/R) * R / 2")
print("   a = h - 2w/R (inside, Meyer), b = (h - a) * area / (N^2 - area) (outside, mean kept)")
print("    R    h  2w/R      w |     a       b | SB inside  SB outside | Ch inside  Ch outside"
      " | SB in-a   (SB out-b)/b | (E_Ch-E_SB)/E_SB | Ch default inside")
conds = [(R, 1.0, q) for q in (0.25, 0.5) for R in (8, 16, 32)] + [(16, 0.6, 0.25), (16, 0.6, 0.5)]
drops = {}
mean_err = 0.0
for R, h, q in conds:
    w = q * R / 2
    f = disc(R, h)
    area = f.sum() / h
    a = h - 2 * w / R
    b = (h - a) * area / (N * N - area)
    inside, outside = rr <= R / 2, rr >= R + 8
    us, uc = split_bregman(f, w), chambolle(f, w)
    ud = denoise_tv_chambolle(f, weight=w)         # library defaults: eps=2e-4, 200 iterations
    mean_err = max(mean_err, abs(us.mean() - f.mean()), abs(uc.mean() - f.mean()))
    drops[(R, h, q)] = h - us[inside].mean()
    print(f"   {R:3d}  {h:.1f}  {q:.2f}  {w:5.2f} | {a:.4f}  {b:.6f} | {us[inside].mean():.5f}"
          f"   {us[outside].mean():.5f}   | {uc[inside].mean():.5f}   {uc[outside].mean():.5f}  "
          f" | {us[inside].mean() - a:+.5f}  {(us[outside].mean() - b) / b:+.4f}      "
          f" | {(energy(uc, f, w) - energy(us, f, w)) / energy(us, f, w):+.2e}"
          f"        | {ud[inside].mean():.5f}")
for q in (0.25, 0.5):
    print(f"   drop h - inside at R=16, 2w/R={q}: h=1.0 {drops[(16, 1.0, q)]:.5f}"
          f"  h=0.6 {drops[(16, 0.6, q)]:.5f}  difference {drops[(16, 0.6, q)] - drops[(16, 1.0, q)]:+.5f}")
R, w = 16, 10.0                                    # 2w/R = 1.25: the disc should vanish
f = disc(R, 1.0)
us, uc = split_bregman(f, w), chambolle(f, w)
print(f"   vanish R=16 w=10 (2w/R=1.25): mean(f) {f.mean():.6f}; SB min {us.min():.6f} max {us.max():.6f};"
      f" Ch min {uc.min():.6f} max {uc.max():.6f}")
mean_err = max(mean_err, abs(us.mean() - f.mean()), abs(uc.mean() - f.mean()))
print(f"   largest |mean(u) - mean(f)| over all disc runs: {mean_err:.1e}")

# ---------------- B. Noise-free ramp: kept except for flat ends ----------------
print("B. Ramp image 64 x 128, f = column / 127; exact discrete solution: flat ends of k+1 pixels")
s = 1 / 127
ramp = np.tile(np.arange(128) * s, (64, 1))
for w in (0.5, 2.0):
    k = int(np.floor((np.sqrt(1 + 8 * w / s) - 1) / 2))   # largest k with s k (k+1) / 2 <= w
    cval = w / (k + 1) + s * k / 2
    exact = ramp[0].copy()
    exact[:k + 1] = cval
    exact[-(k + 1):] = 1 - cval
    us, uc = split_bregman(ramp, w), chambolle(ramp, w)
    print(f"   w={w}: flat end {k + 1} px at {cval:.6f}; max|SB - exact| {np.abs(us - exact).max():.1e};"
          f" max|Ch - exact| {np.abs(uc - exact).max():.1e}")

# ---------------- C. Staircasing on a 1-D ramp, with and without noise ----------------
print("C. 1-D ramp of 256 samples (slope 1/255), w = 0.05, seeds 0-4; interior samples 32-223")
print("   flat step: |u[n+1] - u[n]| < slope/10; plateau: 3 or more samples joined by flat steps")
n = np.arange(256)
clean = n / 255
thr = (1 / 255) / 10
def stairs(u):
    flat = np.abs(np.diff(u[32:224])) < thr
    runs, cur = [], 0
    for fl in flat:
        if fl:
            cur += 1
        else:
            if cur >= 2:
                runs.append(cur + 1)
            cur = 0
    if cur >= 2:
        runs.append(cur + 1)
    return flat.mean(), len(runs), (np.mean(runs) if runs else 0.0)
print("   sigma | input flat frac | TV flat frac  | plateaus | mean plateau len | RMSE input | RMSE TV"
      "   | max|Ch-SB|")
for sigma in (0.0, 0.01, 0.03):
    rows = []
    for seed in range(5):
        g = clean + sigma * np.random.default_rng(seed).standard_normal(256)
        uc, us = chambolle(g, 0.05), split_bregman(g, 0.05)
        fi, _, _ = stairs(g)
        ft, npl, ml = stairs(uc)
        rows.append((fi, ft, npl, ml, np.sqrt(((g - clean) ** 2).mean()),
                     np.sqrt(((uc - clean) ** 2).mean()), np.abs(uc - us).max()))
    r = np.array(rows)
    rng = lambda j, fmt: f"{r[:, j].min():{fmt}}-{r[:, j].max():{fmt}}"
    print(f"   {sigma:.2f}  | {rng(0, '.3f')}     | {rng(1, '.3f')}   | {rng(2, '.0f'):>8s} |"
          f" {rng(3, '.2f'):>16s} | {rng(4, '.4f')} | {rng(5, '.4f')} | {r[:, 6].max():.1e}")

Output:

A. Disc of radius R and height h on a 256 x 256 image, weight w = (2w/R) * R / 2
   a = h - 2w/R (inside, Meyer), b = (h - a) * area / (N^2 - area) (outside, mean kept)
    R    h  2w/R      w |     a       b | SB inside  SB outside | Ch inside  Ch outside | SB in-a   (SB out-b)/b | (E_Ch-E_SB)/E_SB | Ch default inside
     8  1.0  0.25   1.00 | 0.7500  0.000771 | 0.73567   0.00081   | 0.73574   0.00080   | -0.01433  +0.0464       | +2.92e-03        | 0.89921
    16  1.0  0.25   2.00 | 0.7500  0.003106 | 0.74224   0.00319   | 0.74236   0.00318   | -0.00776  +0.0285       | +5.52e-03        | 0.92732
    32  1.0  0.25   4.00 | 0.7500  0.012906 | 0.74587   0.01313   | 0.74612   0.01304   | -0.00413  +0.0172       | +8.60e-03        | 0.97243
     8  1.0  0.50   2.00 | 0.5000  0.001542 | 0.46445   0.00163   | 0.46460   0.00163   | -0.03555  +0.0603       | +7.09e-03        | 0.57873
    16  1.0  0.50   4.00 | 0.5000  0.006212 | 0.47877   0.00646   | 0.47903   0.00642   | -0.02123  +0.0393       | +1.35e-02        | 0.75864
    32  1.0  0.50   8.00 | 0.5000  0.025811 | 0.48686   0.02648   | 0.48746   0.02629   | -0.01314  +0.0261       | +2.24e-02        | 1.00000
    16  0.6  0.25   2.00 | 0.3500  0.003106 | 0.33998   0.00322   | 0.34010   0.00320   | -0.01002  +0.0366       | +1.05e-02        | 0.47413
    16  0.6  0.50   4.00 | 0.1000  0.006212 | 0.07660   0.00650   | 0.07736   0.00646   | -0.02340  +0.0459       | +3.33e-02        | 0.56392
   drop h - inside at R=16, 2w/R=0.25: h=1.0 0.25776  h=0.6 0.26002  difference +0.00226
   drop h - inside at R=16, 2w/R=0.5: h=1.0 0.52123  h=0.6 0.52340  difference +0.00217
   vanish R=16 w=10 (2w/R=1.25): mean(f) 0.012271; SB min 0.012271 max 0.012271; Ch min 0.009432 max 0.017001
   largest |mean(u) - mean(f)| over all disc runs: 1.4e-17
B. Ramp image 64 x 128, f = column / 127; exact discrete solution: flat ends of k+1 pixels
   w=0.5: flat end 11 px at 0.084825; max|SB - exact| 3.9e-15; max|Ch - exact| 2.6e-10
   w=2.0: flat end 23 px at 0.173571; max|SB - exact| 2.9e-10; max|Ch - exact| 1.2e-03
C. 1-D ramp of 256 samples (slope 1/255), w = 0.05, seeds 0-4; interior samples 32-223
   flat step: |u[n+1] - u[n]| < slope/10; plateau: 3 or more samples joined by flat steps
   sigma | input flat frac | TV flat frac  | plateaus | mean plateau len | RMSE input | RMSE TV   | max|Ch-SB|
   0.00  | 0.000-0.000     | 0.000-0.000   |      0-0 |        0.00-0.00 | 0.0000-0.0000 | 0.0023-0.0023 | 5.6e-16
   0.01  | 0.005-0.016     | 0.513-0.560   |    26-32 |        3.44-3.88 | 0.0092-0.0102 | 0.0053-0.0064 | 2.5e-06
   0.03  | 0.000-0.021     | 0.712-0.770   |    31-33 |        4.82-5.58 | 0.0275-0.0305 | 0.0103-0.0132 | 5.3e-08

Every check held. The mean is preserved to $10^{-17}$. Both solvers reproduce the exact ramp at $\lambda = 0.5$, split Bregman to $4 \times 10^{-15}$ and Chambolle’s algorithm to $3 \times 10^{-10}$, and split Bregman also at $\lambda = 2$, to $3 \times 10^{-10}$; the flat ends have the 11 and 23 pixels and the levels 0.084825 and 0.173571 that the formula gives, and the clean 1D input has zero error, as constructed. The prediction that a clean ramp is kept (H7) held through these checks. Split Bregman’s vanishing disc is exactly flat at $\operatorname{mean}(f) = 0.012271$. The two solvers agree on the 1D signals to $2.5 \times 10^{-6}$ or better, and on the inside levels of the discs to within 0.00076, the largest gap being for the low disc at $2\lambda/R = 0.5$.

The disc formula, discretized. The inside levels are close to $a$ but below it in every case: by 0.004 to 0.014 at $2\lambda/R = 0.25$ and by 0.013 to 0.036 at 0.5. H1 held: at each ratio the gap shrinks strictly from $R = 8$ to 16 to 32, by a factor of 1.6 to 1.9 per doubling of the radius, so the discrete minimizer approaches the continuous solution as the disc covers more pixels. In our reading, which the example does not test, a pixelated circle has a larger discrete total variation per unit length than a true circle, so the discrete problem sees a slightly larger perimeter-to-area ratio, and the excess is a boundary effect whose share falls as the radius grows. H2 held too: the outside level is not Meyer’s 0 on the plane but sits within 1.7% to 6.0% of $b$, the level the bounded-domain derivation predicts, so the outside of a $256 \times 256$ image rises by up to 0.026 when the disc loses contrast. The derivation’s two-level form was not proved for the square, and the printed outside levels agree with it to within those percentages. H3 held: the drop at $h = 0.6$ is larger than at $h = 1$ by 0.00226 and 0.00217, well inside 0.005, so the height-independence of the continuous solution survives discretization almost exactly. At $2\lambda/R = 1.25$ the disc vanishes, as derived: split Bregman’s solution is the constant mean, and H4 held.

Solvers and defaults. H5 held: with denoise_tv_chambolle’s defaults (200 iterations, eps $2 \times 10^{-4}$), every inside level is at least 0.114 above split Bregman’s, and at $R = 32$, $\lambda = 8$ the default run leaves the inside level at 1.00000, the input’s height. H6 held in all three parts: with 20,000 iterations Chambolle’s algorithm still ends with an energy 0.3% to 3.3% above split Bregman’s in every disc, leaves the vanishing disc between 0.009432 and 0.017001 instead of flat, and misses the exact ramp by $1.2 \times 10^{-3}$ at $\lambda = 2$, while it reaches the exact ramp at $\lambda = 0.5$. The inside levels in Part A agree with split Bregman’s to within 0.00076 nonetheless, so the energy gap barely shows in the region means the example measures. The larger the weight, the more iterations the dual algorithm needs; with a weight of several units on $[0, 1]$ data its defaults are not a solution of the model, and a library result should be checked for convergence before it is read as the ROF minimizer.

Staircasing is a noise effect. H8 held. On the clean ramp no interior step is flat and there are no plateaus, as the exact solution says. With noise of standard deviation 0.01, between 51% and 56% of the output’s interior steps are flat, against at most 1.6% of the input’s, and the output has 26 to 32 plateaus over 192 samples; at 0.03 the flat share rises to 71–77%, the input’s staying at most 2.1%. H10 held: the plateaus lengthen from a mean of 3.44–3.88 samples to 4.82–5.58. In our reading, the plateaus are the regions where the dual field stays strictly inside $|\mathbf{z}| < 1$, as the optimality condition allows: a run of noise that would cost total variation is cheaper to absorb into a flat segment than to follow. H9 held as well: the output’s RMS error against the clean ramp (at most 0.0064 and 0.0132) is below the input’s (at least 0.0092 and 0.0275) at both noise levels, so the staircase is the price of a real reduction in error, not a failure to denoise. The clean ramp’s output error of 0.0023 is the two flattened ends. This matches the staircase artifact that Getreuer (2012, p. 91) attributes to Nikolova, a tendency to produce small flat regions with artificial edges, which he shows on a noisy image of a smooth dome with one jump; in the example, the clean ramp has no stairs, so they come with the noise.

Properties and Behavior

• Convex, unique, stable. The energy is strictly convex, so there is one minimizer, every correct solver finds it, and Getreuer states it is stable in $L^2$ under perturbations of the data. Different algorithms differ only in speed and in how close they get within a budget, as the example’s Chambolle runs show.
• Edges are kept in place. In the disc solution the edge does not move and stays sharp; only the levels on either side change. The bilateral filter and non-local means also keep edges, but by weighted averaging of similar pixels or patches rather than by solving a global minimization problem.
• Contrast loss depends on scale. A feature loses contrast in proportion to its perimeter over its area, independently of its height, until it disappears. Small, thin and textured structures are lost first, which is why total variation removes fine texture along with noise; the non-local means paper reports that the method noise of total variation minimization is proportional to the curvature of the result, so straight edges are kept and curved ones and texture are degraded.
• Mean preservation. With Neumann boundaries the image mean never changes; the contrast a feature loses reappears as a slight rise of its surroundings.
• Staircasing. Noisy smooth gradients become piecewise constant, while clean ones are kept. The staircases are artificial edges, which, in our reading, edge detectors and segmentation may then treat as real.
• Constant regions are exact. Unlike linear filters, the solution contains regions where the image is exactly flat, which gives total variation results their cartoon-like look.

Limitations

• Texture and fine detail. The piecewise-flat prior treats texture as noise. The non-local means paper assesses total variation minimization, with the other local and transform-domain methods it compared, as degrading fine detail and texture, and reports lower errors for non-local means on its test images; image denoising compares the families, including BM3D and wavelet shrinkage.
• Systematic loss of contrast. Even a perfectly clean disc loses $2\lambda/R$ of its contrast, so measurements of intensity in small structures are biased. Getreuer points to Chan, Esedoğlu, Park and Yip for methods developed to reduce the loss of contrast and staircasing.
• Staircasing on smooth gradients, as shown above.
• One weight for the whole image. A single $\lambda$ cannot suit flat and detailed regions at once; spatially varying weights address this (Variants).
• Gaussian noise model. The quadratic fidelity term assumes Gaussian noise; Getreuer’s impulsive-noise example shows the Gaussian model working poorly where a Laplace (L1) fidelity works well.
• Convergence in practice. The fast first-order solvers need many iterations for large weights, and library defaults may stop early, as the example shows.

Variants

Getreuer (2012, §§1, 2, 5 and 7) describes several, and solves the Laplace and Poisson models with the same split Bregman machinery (§7):

• TV-L1. An absolute-value fidelity, $\lVert u \rVert_{TV} + \lambda_G \int |f - u|$, for Laplace and impulsive noise, which Getreuer attributes to Alliney and to Chan and Esedoğlu. The same pairing of an L1 data term with a total-variation regularizer is the TV-L1 model of optical flow described in Horn–Schunck.
• Poisson fidelity. $\lVert u \rVert_{TV} + \lambda_G \int (u - f \log u)$ for photon-limited data, attributed by Getreuer to Le, Chartrand and Asaki; he also cites extensions to multiplicative and Rician noise.
• Spatially varying weight. $\tfrac12 \int \lambda_G(\mathbf{x}) (f - u)^2$, with TV inpainting as the case where $\lambda_G = 0$ inside the region to be filled.
• Vectorial TV for colour. The gradient magnitudes of the channels are combined under one square root (Getreuer’s Eq. 50), and the grayscale algorithm extends to it directly.
• Anisotropic TV. $\sum |\partial_x u| + |\partial_y u|$, which is not rotation invariant but allows exact combinatorial algorithms; Getreuer names Hochbaum’s algorithm and graph cuts.
• Deblurring and reconstruction. Adding a blur or projection operator to the fidelity gives TV deblurring, used in deconvolution microscopy as image deblurring describes, and the total-variation terms of compressed-sensing MRI. In deep learning, total-variation penalties appear as regularizers in feature visualization and alongside perceptual losses.