Tuesday, October 20, 2015

Week 2: Image Inpainting

Goal


The main goal of Week 2 is to restore images through interpolation. We apply the inpainting algorithm as a re-formulation of the ROF problem that we saw previously in Week 1.

im1-lambda0.01-iROF300-iInp100.jpg

Mandatory Tasks


Effect of Parameters


Lambda

With a greater lambda, we have a faster restoration than with a smaller lambda, but lower fidelity.
Using a greater lambda gives the smoothness term greater importance, which blurs the image, thus causing a loss of contours.



im2-lambda0.01-iROF300-iInp100.jpg
λ = 0.01


im2-lambda0.05-iROF300-iInp100.jpg
λ = 0.05

paramROF.iterMax

This parameter determines the accuracy with which the inpainting is done.
Fewer iterations with ROF minimizes the u image less precisely (iterates too little to give the best solution) so it will give us a suboptimal u. Then, when we use this u in the Inpainting algorithm, it does not give the best solution to the problem, no matter how many times we “paint” the image.

paramInp.iterMax

The number of maximum iterations for the Inpainting function is directly related to the re-calculation of the u and v values. So, the fewer iterations we make, the fewer times we re-calculate v and u, and therefore, the more pixels we are going to restore with the same u and v values.
At the same time, the number of iterations determines the quantity of pixels we are able to restore; with a higher number of iterations we will be able to restore a greater amount of pixels.
This can be seen as a less precise, albeit faster way of restoring the image.

paramROF.dt

This is the time step used when performing derivatives. Therefore, the smaller it is, the closer the results obtained discretely will be to the continuous results.

Justification of Boundary Conditions


We apply separate boundary conditions to three objects:

F and term4grad
For these terms we use Neumann boundary conditions because, as we observed by applying both Dirichlet and Neumann BCs in Block 1, homogeneous Dirichlet boundary conditions would make the borders entirely black by setting the value of the pixels to 0. Neumann boundary conditions allow the borders to be more faithful by constraining the gradient of the border pixels to as opposed to their actual value.

p
The dual variable p is defined such that p =  u/ |u| where u = (ux,uy)
When u = 0, this means that ux = 0 and uy = 0.
Therefore, p1ux / |u| = 0 and p2uy / |u| = 0, which means that the boundary conditions for p must be Dirichlet, as required by the Euler-Lagrange equation.


ROF Seed


The seed for the numerical resolution of the dual ROF model at the first iteration is the noisy image f. At other iterations, the seed is v.

Gradient Descent of Dual Formulation of ROF


To solve the restoration problem in color images, we have to extend the TV calculations to a vectorial space. We consider 3 methods to relate the RGB channels:


Method 1



Method 2



Method 3


If we restore the channels independently (as method 1 does, applying independently TV to the 3 channels and then summing up the results), we are considering the channels as independent problems. We will restore the image, but there will be a lot of noise in the restoration. To perform a good restoration, we should fix two of the channels while we minimize the factorial for the other one, and repeat this operation for all the channels.

The best method to relate the channels is method 3. However, it is computationally difficult and expensive. We have implemented method 1, and with some iterations we could see how the image started to restore. But we need many more iterations to completely solve the problem, and we don’t have enough time nor computational power to perform them.




Restoration of an RGB image using method 1, λ = 0.01, 300 ROF iterations and 800 Inpainting iterations for each channel.

As we can see in this example, the restoration works, but it needs a lot of iterations. Even so, the result still wouldn’t be optimal, because we should use method 3. Big restored areas will be really smooth.


Difference Between Primal & Dual Implementation

tv(u).PNG
The main issue with this expression for total variation used in the ROF model is its non-differentiability at u = 0 because of the modulus and for functions u that are discontinuous.

For the primal implementation of the ROF model for an energy functional J(u), we used the parameter ϵ  (a small constant) to solve this problem.

For the dual implementation, we instead introduce the vector p = (p1,p2).
We can then rewrite the modulus of the gradient of u in the following way:

norm(grad(u)).PNG
This can be further generalized to

sup1.PNG,
so finally,
sup2.PNG

The dual formulation allows for a representation of total variation such that u can be a non-differentiable function because in this form, the gradient does not operate on it.

As for the behavioral difference, the key is in the relationship between the primal and the dual implementations of the same optimization model.

In the primal implementation, for each suboptimal solution that satisfies the regularity conditions, there exists a direction (or several) for it to move in that minimizes the difference between the candidate solution and the constraints.

In the dual implementation, the dual vector p multiplies the constants that determine the positions of the constraints in the primal. A varying p varies the upper bounds in the primal problem, and we are looking for the lowest possible upper bound. So in the dual implementation, the dual vector is minimized in order to get the candidate positions of the constraints closer to the optimum value.

Implementation of Inpainting Algorithm Based on TV Minimization


The code works in the following way:

The start function initializes all the needed parameters, loads the image to restore and the mask, and calls the G4_DualTV_Inpainting_GD function once per channel. The G4_DualTV_Inpainting_GD function is the one that fixes u and v values and calls G4_DualROF_Den_GD function to do the minimization. The G4_DualROF_Den_GD function is the one that minimizes the Jinp(u, v) energy functional.

G4_DualTV_Inpainting_GD

First of all it initializes v.
                

 v=f;

Then, it does the following process paramInp.iterMax times:

At each iteration calls the G4_DualROF_Den_GD fixing v.

u = G4_DualROF_Den_GD(v,paramROF);

Then it fixes u and minimizes again the energy functional. This time there is no need to call the dual ROF algorithm, because  DU is fixed. So the minimization is reduced to equal v to u.


 v=u;

After the minimizations, it restores the original points where we had values in the image to restore.

 v(~mask)=vOld(~mask);

Finally it computes the error between consecutive iterations to apply the stopping criterium.

err=max(abs(vOld(:)-v(:)));

It does this process paramInp.iterMax while the error is over the introduced paramInp.tol parameter.

while (i<paramInp.iterMax) && (err > (paramInp.tol))

G4_DualROF_Den_GD

It does the minimization of the ROF energy functional using Chambolle’s method. Each time it is called, it uses paramROF.iterMax iterations to do the minimization.

Justification PDE Used to Minimize w.r.t v


Once we have a fixed the value of u, the next step is to obtain v. To do so, we operate in the same way as we do to obtain u; minimizing the energy functional.


Since we have already a fixed value for u, the operation is greatly simplified. The first term of the equation is cancelled and therefore, to minimize the functional, it has to be satisfied the next equation:

v = u

Test With Given Images


Image 1

im1-lambda0.01-iROF300-iInp100.jpg
im1-lambda0.05-iROF300-iInp100.jpg

               λ = 0.01         λ = 0.05



ROF iterations 300, Inpainting iterations = 100.

We can see that when we use λ= 0.01, the restoration is more faithful to the original.

Image 2

im2-lambda0.01-iROF300-iInp100.jpg
 λ = 0.01
im2-lambda0.05-iROF300-iInp100.jpg
 λ = 0.05

ROF iterations 300, Inpainting iterations = 100.

Image 2, for which the deterioration is similar to image 1, the result is also better with   λ = 0.005.

Image 3

im3-lambda0.01-iROF200-iInp300.jpg

 λ = 0.01 ROF iterations 200, Inpainting iterations 300

im3-lambda0.05-iROF200-iInp300.jpg

 λ = 0.05 ROF iterations 200, Inpainting iterations 300

Image 3 is more damaged than images 1 and 2, so we need more Inpainting iterations to restore it, and the final result won’t be as good anyway.
In this example, due to computational time limitations, we have decreased ROF iterations to be able to complete 300 Inpainting iterations. The final result is not faithful enough, but a restoration is afforded.

Image 4

im4-lambda0.01-iROF200-iInp300.jpg
 λ = 0.01,  ROF iterations 200, Inpainting iterations 300

im4-lambda0.05-iROF200-iInp300.jpg
 λ = 0.05,  ROF iterations 200, Inpainting iterations 300

Image 4 is as damaged as image 3, so the restoration result is similar.

Thursday, October 15, 2015

Week 1: Image Denoising

Goal

The main goal of Week 1 is to implement variational calculus techniques learned in M2 in a real image denoising problem. This will be achieved by applying the Rudin-Osher-Fatemi Denoising Model, with particular attention to the impact the variation of different parameters have on the final image.

Mandatory Tasks

Effect of Model and Numerical Parameters


The denoising is afforded minimizing an energy functional. We can vary some parameters in that functional in order to control how the minimization is done, and therefore the final result. In this section we study the influence of those parameters.

Lambda
J(u).PNG
The λ is a tradeoff term to control the importance of the regularity and the data fidelity terms in the minimization.
If varying λ we increase the regularity term, the result will be smoother. So we will delete better the noise, but we will smooth object edges.  If we increase the data fidelity term, the result will be more similar to the noisy image. So we will conserve the object edges, but won’t delete as much noise. Consequently, we should give more importance to the data fidelity term when we rely on the image to restore. If the image to restore has a lot of noise, we should give more importance to the regularity term in order to remove it (losing some object edges definition).
When we chose λ we also have to take into account how big is the regularity term chosen. For instance, in the L2-norm regularization method the regularity term is much bigger than in the TV regularization, so we will use λ to compensate that.

λ = 0.05
λ = 0.5
 * In the implementation 1/λ multiplies the fidelity term. So using a higher λ means giving more importance to the regularity term.
Discretization step
From time derivative of image u, used in regularization to determine time step of energy gradient.
Epsilon
The expression for total variation
tv(u).PNG
is not differentiable at the origin or for discontinuous functions u.  To solve this problem, we approximate by a differentiable functional adding the ϵ so the term smooths the peak that the absolute value function has at its minimum in y.
Jrof(u).PNG
Jrofe(u).PNG
If we choose a high ϵ value, the solution will be smoother. But if we choose an ϵ that is too small, the numerical resolution will fail. We don’t want an image that is too smooth (which is something we want to control with the λ parameter instead), but we don’t want the numerical resolution to fail.

Effect of using different seeds


In this section we study the relevance of the chosen seed. The seed we use for the denoising of the image is related with the data fidelity term. For each of the algorithm iterations, we compute the difference between the original image and the one we are looking for (at the beginning is the seed image). Then, the obtained result is used in the next iteration, and so on.
u-f.PNG
At first, the image we are looking for is given a known value (the seed image). Using this image as a starting point, in each of the iterations, we denoise it and compare it with the original image. Thus, the more iterations we make, the difference between the original and the used seed image is less. Therefore, the election of the seed image does not matter as long as we make enough iterations.

Results using 2500 iterations. As we can see in the following images, the chosen quantity of iterations gives us the same results even with different seeds, because they are enough for convergence.

Original
Lenna
Random
Zeros

When we do not run the algorithm enough iterations (200 for example), the obtained image is still similar to the seed image. Therefore, the final result changes for each of the seed images.
Original
Lenna
Random
Zeros

Energy Evolution

To solve the denoising problem we minimize the energy functional:
The lambda is a tradeoff term that should be changed to control the importance of the smoothness and the similarity term in the minimization.
At each iteration, the energy should decrease, which means that we are approaching to the ideal solution of the minimization problem.

L2-norm regularization with random seed and lambda = 15.
In the example above, E(u) decreases until nearly 0. Maybe the denoising problem has been solved correctly, but not necessarily; it’s possible that the chosen lambda is not correct. For instance, if the lambda is too big, maybe the resulting image is too smooth, without edges.

Effect of Boundary Conditions

Neumann boundary conditions specify the values that the gradients of the pixels need to take along the boundary, whereas Dirichlet boundary conditions specify the values that the pixels themselves need to take along the boundary.
Dirichlet boundary conditions imply that the borders are a fixed quantity (0, or black), whereas Neumann conditions imply that the gradient is constant, thus completing the image in a more continuous manner, as it bases the value of the next pixel on the previous.
Tikhonov (λ = 15)


The black border and its expansion due to the blurring can be seen by comparing the energies. With Neumann, the energy approaches zero as the iterations approach infinity. With Dirichlet, the asymptote is at Energy 45.
ROF (λ = 0.5)


We can see that with ROF, the black border when using Dirichlet boundary conditions is more defined and less intrusive in the image than it is with Tikhonov. Mathematically, Dirichlet boundary conditions give the border pixels a constant value. Since the Tikhonov denoising model blurs a lot due to the square value in the smoothness term (p = 2), the edges of the black boundary blur inwards. With ROF, since p = 1, the edges are more defined, so there is less blurring and the constant nature of the Dirichlet boundaries are more obvious.

Solution without Data Fidelity

Without the data fidelity term, the energy functional would be
J(u) without datafidelity.PNG
The minimization would yield the minimum value needed for noise removal, which would blur the image without taking into consideration how similar it might be to the original.
With few iterations, it would simply blur the seed image. After a certain amount of iterations, this would eventually smooth out into a constant color.

Conclusion: Are the results a good solution of the real world problem? 

The results from this denoising method are a good solution of the real world problem. However, energy minimization doesn’t always yield the “best” solution. This is slightly subjective, so we have to look at the solutions ourselves to evaluate them visually, and thus see what parameters were used.