API Reference

TotalVariationImageFiltering.jl provides total-variation denoising primitives for N-dimensional arrays.

Main entry points:

• TVProblem(...) to describe a denoising problem
• solve(problem, config) / solve!(u, problem, config)
• solve_batch(...) / solve_batch!(...)

Abstract supertype for boundary conventions used by differential operators.

Forward-difference gradient with homogeneous Neumann boundary treatment.

Abstract supertype for data-fidelity terms.

L2 data fidelity, 0.5 * ||u - f||^2.

Poisson negative log-likelihood fidelity, sum(u - f * log(u)) up to constants.

Abstract supertype for pointwise convex constraints on the primal variable u.

No additional primal constraint.

Pointwise non-negativity constraint, u[i] >= 0.

Pointwise box constraint, lower <= u[i] <= upper.

Abstract supertype for TV penalty modes.

Isotropic TV, sum(sqrt(sum_d (grad_d u)^2)).

Anisotropic TV, sum(sum_d |grad_d u|).

Abstract supertype for all total-variation solver configurations.

TV-based reconstruction/denoising min_u data_fidelity(u, f) + lambda * TV(u)

spacing[d] is the physical grid spacing along axis d used by differential operators.

TVProblem is immutable, but f is stored by reference. For repeated solves on same shape/eltype data, mutate f in place (for example, copyto!(problem.f, new_f)) and reuse the same TVProblem. If you need a different f array object or different metadata (lambda, spacing, data_fidelity, tv_mode, boundary, constraint), construct a new TVProblem.

Convergence summary for one solver run.

Configuration for the ROF dual-projection solver.

This solver targets the unconstrained Rudin-Osher-Fatemi denoising model with L2Fidelity:

minimize_u 0.5 * ||u - f||_2^2 + lambda * TV(u)

tau is the dual ascent step and must satisfy the stability bound used by this implementation:

tau < 1 / (2 * sum_{d: n_d > 1} h_d^(-2))

where h_d = spacing[d] and n_d = size(f, d).

Fields:

• maxiter: maximum number of iterations.
• tau: dual step size.
• tol: stopping tolerance on relative primal change.
• check_every: evaluate convergence every check_every iterations.

References:

• L. I. Rudin, S. Osher, E. Fatemi, "Nonlinear total variation based noise removal algorithms," Physica D 60(1-4):259-268, 1992. DOI: 10.1016/0167-2789(92)90242-F
• A. Chambolle, "An algorithm for total variation minimization and applications," Journal of Mathematical Imaging and Vision 20:89-97, 2004. DOI: 10.1023/B:JMIV.0000011325.36760.1E

Reusable workspace for the ROF dual projection solver.

ROFState serves two roles:

• scratch buffers (u, u_prev, divp, g, grad_g) reused every call,
• warm-start storage for the dual variable p.

When the same state object is passed to repeated solve! calls, p is not reset between calls.

Configuration for the PDHG / Chambolle-Pock primal-dual solver.

This solver targets TV-regularized models of the form

minimize_u D(u, f) + lambda * TV(u) + I_C(u)

where:

• D is either L2Fidelity or PoissonFidelity,
• C is a pointwise convex set from NoConstraint(), NonnegativeConstraint(), or BoxConstraint(lower, upper).

The primal-dual step sizes must satisfy:

tau * sigma * ||grad||^2 < 1

This implementation uses the conservative bound

||grad||^2 <= 4 * sum_{d: n_d > 1} h_d^(-2)

with h_d = spacing[d] and n_d = size(f, d).

Fields:

• maxiter: maximum number of iterations.
• tau: primal step size.
• sigma: dual step size.
• theta: over-relaxation in [0, 1].
• tol: stopping tolerance on max(relative_primal_change, primal_dual_residual).
• check_every: evaluate convergence every check_every iterations.

References:

• A. Chambolle and T. Pock, "A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging," JMIV 40:120-145, 2011. DOI: 10.1007/s10851-010-0251-1

Reusable workspace for the PDHG / Chambolle-Pock solver.

PDHGState serves two roles:

• scratch buffers (u, u_prev, u_bar, divp, primal_tmp, grad_u_bar) reused every call,
• warm-start storage for the dual variable p.

When the same state object is passed to repeated solve! calls, p is not reset between calls.

Summary returned by select_lambda_discrepancy.

Fields:

• lambda: selected TV weight.
• u: denoised image/volume at lambda.
• solve_stats: solver statistics for u.
• residual_norm2: ||u - f||_2^2.
• target_norm2: discrepancy target (target_scale * length(f) * sigma^2).
• relative_mismatch: relative residual-target mismatch.
• evaluations: number of ROF solves used for parameter search.
• bracket: last (lambda_lo, lambda_hi) search interval.

Summary returned by select_lambda_sure.

Fields:

• lambda: selected TV weight.
• u: denoised image/volume at lambda.
• solve_stats: solver statistics for u.
• sure: selected SURE value.
• residual_norm2: ||u - f||_2^2.
• divergence: MC estimate of divergence at selected lambda.
• epsilon: finite-difference perturbation scale used for MC-SURE.
• evaluations: number of ROF solves used for parameter search.
• lambda_grid: tested lambda grid (sorted ascending).
• sure_values: SURE values aligned with lambda_grid.
• residual_values: residual values aligned with lambda_grid.
• divergence_values: divergence estimates aligned with lambda_grid.

Solve a TV problem and return (u, stats).

config selects the solver backend. Implemented: ROFConfig, PDHGConfig.

Run ROF denoising in place.

u is both the initial guess and output buffer.

Constraint support:

• ROFConfig supports only problem.constraint = NoConstraint().
• For constrained problems, use PDHGConfig.

State and buffer reuse:

• Reuse u and state across calls to avoid allocations.
• u is copied into solver state as the initial output buffer.
• For ROF iterations, the effective warm start is the retained dual variable state.p.
• state.p is reused as dual warm start and is not reset by solve!.
• For new image data with the same array storage, update in place with copyto!(problem.f, new_f) and call solve! again.
• If image array object, shape/eltype, or problem metadata changes, create a new TVProblem (and a new ROFState only when shape/eltype changes).

The implementation uses a Chambolle-style dual projected-gradient method for the ROF model. With dual variable p and image f, each iteration computes:

g^k = div(p^k) - f / lambda

p^(k+1) = Proj_B(p^k + tau * grad(g^k))

u^(k+1) = f - lambda * div(p^(k+1))

where Proj_B is projection onto the TV dual unit ball:

• isotropic TV: ||p[i]||_2 <= 1 per pixel/voxel,
• anisotropic TV: |p[d][i]| <= 1 per component.

gradient! and divergence! use forward/backward finite differences with homogeneous Neumann boundary handling and account for problem.spacing.

References:

• L. I. Rudin, S. Osher, E. Fatemi, "Nonlinear total variation based noise removal algorithms," Physica D 60(1-4):259-268, 1992. DOI: 10.1016/0167-2789(92)90242-F
• A. Chambolle, "An algorithm for total variation minimization and applications," Journal of Mathematical Imaging and Vision 20:89-97, 2004. DOI: 10.1023/B:JMIV.0000011325.36760.1E

Run TV denoising/reconstruction in place using PDHG / Chambolle-Pock.

u is both the initial guess and output buffer.

State and buffer reuse:

• Reuse u and state across calls to avoid allocations.
• u is copied into solver state and acts as primal warm start.
• state.p is reused as dual warm start and is not reset by solve!.
• For new image data with the same array storage, update in place with copyto!(problem.f, new_f) and call solve! again.
• If image array object, shape/eltype, or problem metadata changes, create a new TVProblem (and a new PDHGState only when shape/eltype changes).

Stopping criterion:

• solve! checks both relative primal change and a PDHG primal-dual residual.
• The primal-dual residual is normalized by sqrt(length(u)).
• Convergence uses max(relative_primal_change, primal_dual_residual) <= tol.

Constraints:

• If problem.constraint is not NoConstraint(), the primal update applies the exact pointwise constrained proximal map for the supported fidelities.

Solve a batch of same-size images stored in one array.

The last axis is interpreted as batch index, and TV operators act only on the leading spatial axes.

Return values:

• Default (return_per_item_stats = false): (u_batch, summary_stats).
• With return_per_item_stats = true: (u_batch, summary_stats, per_item_stats), where per_item_stats is a Vector{SolverStats} for each batch item.

Select lambda for ROF denoising using Morozov's discrepancy principle.

This function finds lambda >= 0 such that the ROF residual satisfies

||u_lambda - f||_2^2 ≈ target_scale * length(f) * sigma^2,

with u_lambda solving

min_u 0.5 * ||u - f||_2^2 + lambda * TV(u).

The search strategy is:

1. bracket the target residual by expanding lambda_max,
2. run bisection in lambda up to max_bisect iterations.

Works for any dimensionality (N-D arrays), including 2D images and 3D volumes.

References:

• V. A. Morozov, Methods for Solving Incorrectly Posed Problems, 1984. DOI: 10.1007/978-1-4612-5280-1
• Y. Wen and R. H. Chan, "Parameter selection for total-variation based image restoration using discrepancy principle," IEEE TIP 21(4):1770-1781, 2012. DOI: 10.1109/TIP.2011.2181401

Select lambda for ROF denoising by minimizing Monte-Carlo SURE on a grid.

For each lambda in lambda_grid, this evaluates

SURE(lambda) = -N*sigma^2 + ||u_lambda - f||_2^2 + 2*sigma^2*div(u_lambda(f)),

where divergence is approximated with Monte-Carlo finite differences:

div(u_lambda(f)) ≈ (1/epsilon) * b' * (u_lambda(f + epsilon*b) - u_lambda(f)), with b ~ N(0, I) and averaged over mc_samples.

Works for any dimensionality (N-D arrays), including 2D images and 3D volumes.

References:

• C.-A. Deledalle et al., "Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection," 2014. HAL: hal-00987295
• S. Ramani, T. Blu, M. Unser, "Monte-Carlo SURE: A black-box optimization of regularization parameters for general denoising algorithms," IEEE TIP 17(9):1540-1554, 2008. DOI: 10.1109/TIP.2008.2001404
• Y. Lin, B. Wohlberg, H. Guo, "UPRE method for total variation parameter selection," Signal Processing 90(8):2546-2551, 2010. DOI: 10.1016/j.sigpro.2010.02.025

Allocate a tuple of N arrays with the same shape/type as u to store vector-valued fields.

In-place forward finite-difference gradient with Neumann boundary treatment.

For each dimension d, g[d][i] = (u[i + e_d] - u[i]) / spacing[d] on interior points and 0 on the upper boundary.

In-place divergence operator that is negative-adjoint consistent with gradient! under the same boundary convention.

Project dual field p onto the isotropic TV dual ball:

||p[i]||_2 <= radius for every pixel/voxel i.

Project dual field p onto the anisotropic TV dual ball:

|p[d][i]| <= radius for every component and pixel/voxel i.