Fdtxjr

I am watching some wonderful videos on Variational Methods in Image Processing, and the presenter talks about how variational methods are used to de-blur or denoise images, as well as other applications--including "active contours." Variational image processing has been an active research program for some time. One example model the presenter discusses is the Rudin, Osher, Fatemi model for image denoising. The energy function for this minimization problem is given below.

$$
E(u) = int_Omega(u - f)^2 + lambda |nabla u|dx
$$

In the equation, $u$ is a function which generates some data and exhibits some properties as per the first and second terms. The presenter indicated that the Variational problem is minimized using gradient descent methods for optimization.

My question was, does anyone know of a good tutorial on numerically implementing gradient descent on this type of optimization problem? I want to understand exactly how the computer will compute the solution to this problem. Gradient Descent is easy enough to compute, but I was not sure about optimizing over the discretized infinite dimensional space of functions. The functional analysis is not the problem, just understand how the computer algorithm is coded to make this work.

I found an article by Getreuer which talks about the optimization algorithm using the Split Bregman method, which also has some C code. In looking at the code though, it is a bit hard to understand the underlying algorithm versus all of the different function calls and typedefs, and all.

I was hoping someone might know a good online tutorial or video, etc., that discusses how to implement this type of optimization routine. Thanks.

The basic idea is to right the problem in a discrete manner.

Let's start with a 1D example.

Imagine the discrete realization of $ f $ is given by the vector $ boldsymbolf in mathbbR^n $.

The problem becomes:

$$ hat boldsymbolu = arg min_ boldsymbolu frac12 boldsymbolu - boldsymbolf right_2^2 + lambda left| D boldsymbolu right|_1 $$

Where $ D $ is finite differences operator in Matrix form such that:

$$ left| D boldsymbolu right|_1 = sum_i = 2^n left| boldsymbolu_i - boldsymbolu_i - 1 right| $$

This is basically Linear Least Squares with Total Variation Regularization.

There are may methods to solve this, for instance, see The Dual Norm of Total Variation Norm (Form of ⟨⋅,⋅⟩) By Smoothing.

Modern way of solving this is using the Alternating Direction Method of Multipliers (ADMM) with great resource on this subject would be Professor Stephen P. Boyd page on ADMM.

Assume the solution to that is given by a Black Box - 1D TV Denoising.

I will point out more than that, let's assume our solution can solve the problem above for any matrix $ D $ (Not specifically Forward / Backward Differences operator).

We'll get to a specific solution and code for it later, but let's assume we have it.

We have in image in 2D space but we could vectorize it into a column stack form using the Vectorization Operator.

It is easy to see the Fidelity Term (The $ L_2 $ term) in our optimization works well. But we have to deal with the 2nd term - The TV Operator.

The TV operator could be defined in (At least?) 2 ways in 2D.

The 2 form of the TV Operators in 2D is $ mathbfTV left( X right) = mathbbR^m times n to mathbbR $

1. 2D Anistoropic Total Variation
   $$ mathbfTV_I left( X right) = sum_i = 1^m - 1 sum_j = 1^n - 1 sqrt left( X_i, j - X_i + 1, j right)^2 + left( X_i, j - X_i, j + 1 right)^2 + sum_i = 1^m - 1 left| X_i, n - X_i + 1, n right| + sum_j = 1^n - 1 left| X_m, j - X_m, j + 1 right| $$
   


2. 2D Isotropic Total Variation

$$ mathbfTV_A left( X right) = sum_i = 1^m - 1 sum_j = 1^m left| X_i, j - X_i + 1, j right| + sum_i = 1^m sum_j = 1^n - 1 left| X_i, j - X_i, j + 1 right| $$

Remark

One could come up with many more variations. For instance, include the diagonals as well.

I will ignore the Isotropic version for this case (Though it is not hard to solve it as well).

The Anisotropic Case is very similar to the 1D case with on difference, we need to create 2 subtractions for each pixel.

So with some tweaking of the matrix we can write is in the same model of the 1D case.

We just need to update $ D $ to be $ D : mathbbR^m n to mathbbR^2 m n - m - n $.

So we need a code to build the matrix $ D $ in the 2D case and the solver to the 1D case and we're done.

2D Gradient Operator

Work in Progress

Results

So, the implemented code yields as following (The reference being CVX based solution):

As can be seen above, the numerical solver can achieve the optimal solution.

The code is available at my StackExchange Mathematics Q3164164 GitHub Repository.

Black Box Solver

Relatively simple yet effective method is An Algorithm for Total Variation Minimization and Applications.

The idea is the use of Dual Norm / Support Function which suggests that:

$$ _1 = max_boldsymbolp left boldsymbolp^T D boldsymbolx mid boldsymbolp right_infty leq 1 right $$

Hence our problem becomes:

$$ arg min_ boldsymbolu frac12 boldsymbolu - boldsymbolf right_2^2 + lambda left| D boldsymbolu right|_1 = arg min_ boldsymbolu max_ left_infty leq 1 frac12 boldsymbolu - boldsymbolf right_2^2 + lambda boldsymbolp^T D boldsymbolu $$

Using the Min Max Theorem (The objective is Convex with respect to $ boldsymbolu $ and Concave with respect to $ boldsymbolp $) one could switch the order of the Maximum and Minimum which yields:

$$ arg max_ left_infty leq 1 min_ boldsymbolu frac12 boldsymbolu - boldsymbolf right_2^2 + lambda boldsymbolp^T D boldsymbolu $$

Which can be solved pretty easily with Projected Gradient Descent.

Remarks

• The solver in its FISTA (A Fast Iterative Shrinkage Thresholding Algorithm for Linear Inverse Problems) form is much faster with a simple tweak.


• The code use matrix form for simplicity. In practice all operation should be done using Convolution in the case of the TV.