Nonlinear Total Variation based noise removal algorithms

Nonlinear Total Variation based noise

removal algorithms

最近仔细读了下提出全变差（Total Variation，TV）最小化去噪的那篇文章《Nonlinear Total

Variation based noise removal algorithms》[1]（ROF92）。

在了解TV去噪的过程中，我发现，不管是期刊论文、博硕士论文还是博客，很少有进一步去将文献[1]讲得更加清楚一些的。再加上文献[1]本身公式大多都缺少中间过程、推导过程不易理解，给出结论都非常突兀。还有便是，Guy Gilboa副教授的开源代码中最关键的u 和λ 的更新公式都是不同于文献[1]的公式(2.8a)和(2.9c)的。这些都给水平不高的我造成了很多理解上的困惑。所幸，在大师兄的帮助下，我算是看懂文献[1]。另外，Guy Gilboa副教授的开源代码关键变量的更新公式虽然不同于文献[1]，但也提供了另外一个视角，帮助我对算法的理解更上一层楼。

针对上面提到的问题，我打算把我自己看论文和写程序时想明白的写一下，预计是3篇博文，分别解决以下3个问题：

1.

2.

3.

如何将TV去噪问题从范数约束形式通过Euler-Lagrange方程变为偏微分形式？

如何推导出论文中直接给出的那些公式？

Guy Gilboa副教授的开源代码更新公式表达与论文不同，为何还是正确的？

我们先来解决第1个问题。实际上，第1个还能够进一步细分为3个子问题：

1.

2.

3.

如何得到公式(2.5a)？

如何得到公式(2.6)?

公式(2.8a)是不是有错？

1. 如何将TV去噪问题从范数约束形式通过Euler-Lagrange方程变为偏微分形式？

文献[1]中，公式(2.5a)的提出很突兀，这一小节主要讲述如何导出公式(2.5a)：

ut=∂∂x⎛⎝⎜uxu2x+u2y−−−−−−√⎞⎠⎟+∂∂y⎛⎝⎜uyu2x+u2y−−−−−−√⎞⎠⎟−λ(u−u0)(2.5a)

众所周知，TV去噪问题可以表示成如下形式：

minu∫Ωu2x+u2y−−−−−−√+(u−u0)2dxdy,(my−1)

这里，u 是去噪后也就是待恢复的图像，u0 是噪声图像，一般默认为高斯白噪声。文献[1]中的公式仍然采用原标号，出于解释得更加清楚而我自己加的公式，为了以示区别，会在前面加 my- 。事实上，(my-1)所代表的是各向同性的TV（isotropic TV），还存在各向异性的TV（anisotropic TV），公式如下：

minu∫Ω|ux|+|uy|+(u−u0)2dxdy,(my−2)

不过文献[1]只考虑各向同性的TV，我们这里也只针对(my-1)。要从公式(my-1)得到文献[1]中的公式(2.5a)我们需要借助Euler-Lagrange公式，如下所示（一维情况）：

∂F∂f−ddx∂F∂f′=0,(my−3)

这里

Fffxfy=u2x+u2y−−−−−−√+(u−u0)2=u=ux=uy

把上面的变量代入Euler-Lagrange公式，可以得到

∂F∂f∂F∂fx∂F∂fy∂∂x∂F∂fx∂∂y∂F∂fy=d(λ2(u−u0)2)du=λ(u−u0)=du2x+u2y−−−−−−√dux=uxu2x+u2y−−−−−−√=du2x+u2y−−−−−−√duy=uyu2x+u2y−−−−−−√=∂∂x⎛⎝⎜uxu2x+u2y−−−−−−√⎞⎠⎟=∂∂y⎛⎝⎜uyu2x+u2y−−−−−−√⎞⎠⎟

所以就能得到：

0=∂∂x⎛⎝⎜uxu2x+u2y−−−−−−√⎞⎠⎟+∂∂y⎛⎝⎜uyu2x+u2y−−−−−−√⎞⎠⎟−λ(u−u0)

上式和(2.5a)相比，就是左边一个是0 ，一个是ut 。其实吧，这个0 只有在达到稳态条件，也就是图像不再变换时才能取到，如果还没到稳态条件，这就是ut ，由此我们得到了公式(2.5a)：

ut=∂∂x⎛⎝⎜uxu2x+u2y−−−−−−√⎞⎠⎟+∂∂y⎛⎝⎜uyu2x+u2y−−−−−−√⎞⎠⎟−λ(u−u0)(2.5a)

不同的论文里，公式(2.5a)也可以有其他两种本质一模一样的形式：

ut=∇⋅(∇u|∇u|)−λ(u−u0)(2.5a)

或者

ut=div(∇u|∇u|)−λ(u−u0)(2.5a)

至此，我们就将TV去噪问题从范数约束形式通过Euler-Lagrange方程变为偏微分形势。

去噪的迭代公式为：

un+1=un+Δt⋅ut

2. 如何得到公式(2.6)?

文献[1]提到了公式(2.6)的得到，只需要在公式(2.5a)左右乘上u−u0 并对整幅图像积分就可以了，但过程还是不容易理解。这一小节就是要把这个过程讲清楚。先放一下我们最终的目标，公式(2.6)：

λ=−12σ2∫⎡⎣⎢u2x+u2y−−−−−−√−⎛⎝⎜(u0)xuxu2x+u2y−−−−−−√+(u0)yuyu2x+u2y−−−−−−√⎞⎠⎟⎤⎦⎥dxdy(2.6)

在公式(2.5a)左右乘上u−u0 并对整幅图像积分我们得到：

∫Ωut(u−u0)dxdy=∫Ω⎡⎣⎢∂∂x⎛⎝⎜uxu2x+u2y−−−−−−√⎞⎠⎟+∂∂y⎛⎝⎜uyu2x+u2y−−−−−−√⎞⎠⎟⎤⎦⎥(u−u0)−λ(u−u0)2dxdy(2.5a−1)

在达到稳态条件时，ut 应该为零（u 不在变化了），所以上式的左边为0 。

又因为根据(2.3c)公式

∫Ω12(u−u0)2dxdy=σ2(2.3c)

公式(2.5a-1)可以变为

λ=12σ2∫Ω⎡⎣⎢∂∂x⎛⎝⎜uxu2x+u2y−−−−−−√⎞⎠⎟+∂∂y⎛⎝⎜uyu2x+u2y−−−−−−√⎞⎠⎟⎤⎦⎥(u−u0)dxdy(2.5a−2)

进一步的，我们把公式(2.5a-2)写得更加简洁一点：

λ=12σ2∫Ω∇⋅(∇u|∇u|)(u−u0)dxdy(2.5a−3)

其实，(2.5a-3)可以看作是下面这个分部积分中的一部分：

∫Ω∇⋅(∇u|∇u|(u−u0))dxdy=∫Ω∇⋅(∇u|∇u|)(u−u0)dxdy+∫Ω∇u|∇u|∇(u−u0)dxdy

通过，高斯公式，我们可以得到：

∫Ω∇⋅(∇u|∇u|(u−u0))dxdy=∫∂Ω∇u|∇u|(u−u0)dl⃗

利用文献[1]公式(2.5c)的边界条件，

∂u∂n=0

我们能够得到：

∫Ω∇⋅(∇u|∇u|(u−u0))dxdy=∫∂Ω∇u|∇u|(u−u0)dl⃗ =0

所以

∫Ω∇⋅(∇u|∇u|)(u−u0)dxdy=−∫Ω∇u|∇u|∇(u−u0)dxdy

由此，公式(2.5a-3)就可以变成

λ=12σ2∫Ω∇⋅(∇u|∇u|)(u−u0)dxdy=−12σ2∫Ω∇u|∇u|∇(u−u0)dxdy