Variational Image Denoising is a technique used to remove noise from an image by minimizing a cost function. The cost function is typically chosen to be the sum of two terms: a data fidelity term and a regularization term. The data fidelity term measures the difference between the noisy image and the denoised image, while the regularization term enforces certain properties on the denoised image, such as smoothness.

## Noise Models

When dealing with Gaussian noise, the data fidelity term is typically chosen to be the mean squared error (MSE) between the noisy image and the denoised image. The regularization term is chosen to promote smoothness in the denoised image, and can be implemented using a total variation (TV) regularizer. The TV regularizer promotes piecewise constant solutions, which results in a denoised image that is piecewise constant and preserves edges. Let $$f\in\mathbb{R}^n$$ be a noisy image and $$\mathbb{K}:\mathbb{R}^n\to\mathbb{R}^{m\times 2}$$ a discrete gradient operator, the denoised image $$u^*\in\mathbb{R}^n$$ is the solution of the following optimization problem:

$$u^* = arg\min_u \frac12 \|u-f\|^2 + \alpha\sum_{i=1}^m\|(\mathbb{K} u)_j\|,$$ where $$\alpha\ge0$$ is a balance parameter that controls the contribution of the regularization term in the final reconstruction.

When dealing with Poisson noise, the data fidelity term is typically chosen to be the negative log-likelihood of the Poisson distribution. The regularization term is again chosen to promote smoothness in the denoised image, and can be implemented using a total variation (TV) regularizer.

$$u^* = arg\min_u \sum_{i=1}^n (u_i - f_i \log(u_i)) + \alpha\sum_{i=1}^m\|(\mathbb{K} u)_j\|$$

There are many algorithms that can be used to minimize the cost function for variational image denoising. One popular method is the alternating direction method of multipliers (ADMM), which has been shown to be effective for both Gaussian and Poisson noise. Another method is the primal-dual algorithm, which has also been shown to be effective for both Gaussian and Poisson noise.

## Proximal of the Total Variation Regularizer

The proximal operator of the TV seminorm is a key ingredient in the implementation of a numerical method for solving this optimization problem numerically. Let us recall the definition of the proximal operator of a function $$f$$: $$prox_{\tau,f}(v) = arg\min_v f(v) + \frac{1}{2\tau}\|x-v\|^2$$

## Chambolle-Pock for Gaussian Denoising

Here, to use this method we need to reformulate the optimality conditions for the gaussian denoising problem as a saddle point problem

import numpy as np

def chambolle_pock_denoising(image, lambda_, theta=0.25, max_iter=100, epsilon=1e-4):
"""
Denoises an image using the Chambolle-Pock algorithm.
:param image: numpy array representing the image to be denoised
:param lambda_: regularization parameter for the total variation seminorm
:param theta: overrelaxation parameter
:param max_iter: maximum number of iterations
:param epsilon: stopping criterion for the relative change in the solution
:return: denoised image
"""
# Define the gradient operator
gradient = lambda x: np.dstack((np.diff(x, axis=0), np.diff(x, axis=1)))

# Initialize the primal and dual variables
u = np.zeros_like(image)

# Iterate until the maximum number of iterations is reached
for i in range(max_iter):
u_prev = u
p = p + theta * gradient_u
p = np.clip(p, -lambda_, lambda_)
div_p = np.clip(np.cumsum(p, axis=0), 0, np.inf)
div_p = np.clip(np.cumsum(div_p, axis=1), 0, np.inf)
u = image - div_p
u = np.clip(u, 0, np.inf)

# Check the stopping criterion
relative_change = np.linalg.norm(u - u_prev) / np.linalg.norm(u)
if relative_change < epsilon:
break

return u


## Chambolle-Pock for Poisson Denoising

import numpy as np

def chambolle_pock_denoising(image, lambda_, theta=0.25, max_iter=100, epsilon=1e-4):
"""
Denoises an image using the Chambolle-Pock algorithm.
:param image: numpy array representing the image to be denoised
:param lambda_: regularization parameter for the total variation seminorm
:param theta: overrelaxation parameter
:param max_iter: maximum number of iterations
:param epsilon: stopping criterion for the relative change in the solution
:return: denoised image
"""
# Define the gradient operator
gradient = lambda x: np.dstack((np.diff(x, axis=0), np.diff(x, axis=1)))

# Initialize the primal and dual variables
u = np.zeros_like(image)

# Iterate until the maximum number of iterations is reached
for i in range(max_iter):
u_prev = u
p = p + theta * gradient_u
p = np.clip(p, -lambda_, lambda_)
div_p = np.clip(np.cumsum(p, axis=0), 0, np.inf)
div_p = np.clip(np.cumsum(div_p, axis=1), 0, np.inf)
u = image - div_p
u = np.clip(u, 0, np.inf)

# Check the stopping criterion
relative_change = np.linalg.norm(u - u_prev) / np.linalg.norm(u)
if relative_change < epsilon:
break

return u


In summary, variational image denoising is a powerful technique that can be used to remove noise from an image while preserving important features such as edges. It can be applied to both Gaussian and Poisson noise, and various algorithms such as ADMM and primal-dual algorithm can be used to minimize the cost function.