Monday 14 October 2013

Total Variation Denoising

Lately I became interested in optical flow estimation and met Thomas Pock who answered nearly every of my questions by "Yes, this can be done using variational methods" :) So I started to look for papers, found out there are many and was a bit lost until I found the paper: A. Chambolle, T. Pock. A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging. Journal of Mathematical Imaging and Vision, 2011. Before implementing the optical flow I implemented all the examples for denoising, which is based on the same algorithm and is easier to grasp at the beginning.

The Matlab code can be downloaded from GitHub:

git clone https://github.com/JanSochman/TVdenoising.git

Feel free to use it.


Implemented total variational models for denoising
  • ROF model (preserves discontinuities due to L1 norm) $$ \min_{u\in X} ||\nabla u||_1 + \frac{\lambda}{2}||u - g||_2^2 $$
  • TV-L1 ROF (adds robustness to noise outliers) $$ \min_{u\in X} ||\nabla u||_1 + \lambda||u - g||_1 $$
  • Huber-ROF (avoids staircasing effect of TV methods) $$ \min_{u\in X} |\nabla u|_\alpha + \frac{\lambda}{2}||u - g||_2^2 \quad \quad |x|_\alpha = \left\{\begin{array}{lr} |x|^2/(2\alpha), \, & |x| \leq \alpha\\ |x| - \alpha/2, \, & |x| > \alpha \end{array} \right. $$

Where $g$ is the input image, $u$ the sought solution, and $X$ is the set of $M\times N$ images.

The optimisation
  • All problems are convex but only Huber-ROF is smooth
  • All problems can be expressed as $$ \min_{x\in X} F(Kx) + G(x) $$ with corresponding primal-dual saddle-point problem being $$ \min_{x\in X} \max_{y\in Y} \langle Kx,y\rangle + G(x) - F^{*}(y) $$
  • This can be solved by an iterative primal-dual algorithm (see the paper)
    • Initialisation: Choose $\tau$, $\sigma \gt 0$, $\theta \in [0, 1]$, $(x^0,y^0)\in X\times Y$ and set $\bar{x}^0 = x^0$.
    • Iterations ($n\ge 0$): Update $x^n$, $y^n$, $\bar{x}^n$ as follows: $$\left\{ \begin{array}{l} y^{n+1} = (I + \sigma\partial F^*)^{-1}(y^n + \sigma K \bar{x}^n)\\ x^{n+1} = (I + \tau \partial G)^{-1}(x^n-\tau K^*y^{n+1})\\ \bar{x}^{n+1} = x^{n+1} + \theta (x^{n+1}-x^n)\end{array}\right.$$
  • Two more variants with faster convergence proposed in the paper in case $G$ or/and $F^*$ are uniformly convex (used when applicable)

Methods comparison: salt & pepper noise

original image salt & pepper noise ROF (non-robust)
TV-L1 ROF
(robust but staircasing)
Huber ROF
(smooth but non-robust)
Huber-L1 ROF
(smooth and robust)


Methods comparison: gaussian noise

original image gaussian noise ROF (non-robust)
TV-L1 ROF
(robust but staircasing)
Huber ROF
(smooth but non-robust)
Huber-L1 ROF
(smooth and robust)

5 comments:

  1. Hi,

    I am trying to run the code on my matlab environment and I am not getting any denoising result. the method leaves the image with noise untouched.

    Following is the sequence of commands I used

    >i1 = imread('D:\projects\libraries\opencv\samples\cpp\lena.jpg')

    >i1g = rgb2gray(i1)

    >i1g_noise = imnoise(i1g,'salt & pepper',0.25)

    >i1gd_noise = im2double( i1g_noise ) //I also tried simply converting the image to double

    >[denoised_ig,criteria]=TVdenoising(i1gd_noise,HuberROFalg3,100,5,i1g,0.05,0)

    >i1gd_denoise = reshape( i1gd_Denoise , 512 , 512)

    >imshow(i1gd_denoise)


    The result I get remains completely untouched. Can you please tell me , what mistake I am making.

    Regards

    Avanindra

    ReplyDelete
  2. I have tried exactly what you are trying and it works for me. The only problem I had was that your code is not syntactically correct. For instance the parameter HuberROFalg3 should be 'HuberROFalg3' and also you save the result into denoised_ig, but then work with i1gd_Denoise. Could this be a problem?

    ReplyDelete
  3. thank you for the response. 'HuberROFalg3' is working now. I was making some mistake before.

    I need to ask one more question though , does 'HubeROFalg3' stands for TVL1 with Huber-ROF model , using the algorithm 3 mentioned in the paper?..

    Thanks

    Avanindra

    ReplyDelete
  4. Also I noticed , the implementation was working only when the color range are normalized to [0 - 1] range. If I pass the gray image with with color range 0-255 , the the denoising algorithm leaves the noise image untouched. This was the mistake I initially was doing.

    ReplyDelete
  5. Yes, you are right. The images are supposed to be normalised to [0, 1].

    HuberROFalg3 is not using TVL1 as you can see in TVdenoising.m. Unfortunately, not all the combinations are possible due to convexity requirements of the individual algorithms. See the original paper for more details.

    Anyway, thank you for trying it out and I hope it will be useful for you :)

    ReplyDelete