TV去噪的拙劣理解以及python实现

TV去噪的原理描述为：

受噪声污染的图像总变分比无噪图像总变分大

啥是总变分？根据以下公式：

$J(u)={\iint }_{x,y\in image}^{}|▽u|dxdy={\iint }_{x,y\in image}^{}\sqrt{u_x^2+u_y^2}dxdy$

其中$u_x$,$u_y$是像素在$x$,$y$方向的梯度

可知总变分是梯度幅值的积分，为了消除噪声，需要使总变分（$J$）最小化。

公式推导看着晕，这到底什么意思？

根据公式的字面理解，是要让全图梯度值之和最小？那纯色图的梯度值之和是最小的，那确实也没有噪声了。

难道只是说噪声图像的梯度变化更大，那不是废话。。。

实际上为了使去噪结果更接近原图，避免失真，推导时会加入保真项，最后的推导结果是：

$\lambda (u-u_0)-▽(\frac{▽u}{|▽u|})=0$

这就是$J$取极值的解，以下的思考与其它地方查到的不同，不知是否更准确：

• 当在平滑区域，梯度的梯度同样是0，即$▽(\frac{▽u}{|▽u|})=0$，此时$u=u_0$，$u_0$就是原图。
• 当在图像边缘时，梯度$▽u$变大，梯度的梯度具有一定的值，$▽(\frac{▽u}{|▽u|})$值较小，此时$u\approx u_0$。
• 当在噪声区域中时，突变的噪点会使梯度的梯度增加，此时$▽(\frac{▽u}{|▽u|})$较大，就起到了去噪的效果。

再举个例子，假如离散序列为【100，100，500，500，500】
梯度是【0，400，0，0，？】（$x_{n+1}-x_n$）
梯度的梯度是【400，-400，0，？，？】

假如有噪点，那么离散序列是：【100，100，500，100，100】
那么梯度是【0，400，-400，0，？】
梯度的梯度是【400，-800，400，?，？】

明显噪点的$▽(\frac{▽u}{|▽u|})$更可观

（以上猜想概不负责）---------------------------------------------------------------------

最后根据论文【变分图像去噪】，们得到了一个迭代公式：

$u_{x,y}^{n+1}=u_{x,y}^n-△t\lambda (u_{x,y}^n-u_{x,y}^0)+△t(divp)$

其中$n$代表迭代次数，$n=0$时，$u^0$代表输入图，$△t$是时间步长参数，类似权重，$\lambda$是正则化参数，$divp$是扩散项

$divp=▽(\frac{▽u_{x,y}^n}{|▽u_{x,y}^n|})=\frac{u_y^2{u}_{xx}-2u_x{u}_y{u}_{xy}+u_x^2{u}_{xy}}{u_x^2+u_y^2}$

$n=0$时，$divp=0$

每次迭代，都要用上一次的输出图像$u^n$减去$△t(\lambda (u_{x,y}^n-u_{x,y}^0)-divp)$，理论上当它跟之前的解一样等于0时，每次迭代的结果相同，意味着噪声被去除了。

• 如果某个点的divp较大，那么迭代输出与前一次的差距就更大，起到了去噪的作用。

• 如果divp=0，那么公式变为：$u_{x,y}^{n+1}=(1-△t\lambda )u_{x,y}^n-△t\lambda u_{x,y}^0$，也就是$u_{x,y}^n$与$u_{x,y}^0$的加权和。

然后$u$的一系列导数计算公式为（我看到有一篇文章搞错了加减号）：

$(u_x{)}_{i,j}^n=\frac{u_{i+1,j}^n-u_{i-1,j}^n}{2}$

$(u_y{)}_{i,j}^n=\frac{u_{i,j+1}^n-u_{i,j-1}^n}{2}$

$(u_{xx}{)}_{i,j}^n=u_{i+1,j}^n+u_{i-1,j}^n-2u_{i,j}^n$

$(u_{yy}{)}_{i,j}^n=u_{i,j+1}^n+u_{i,j-1}^n-2u_{i,j}^n$

$(u_{xy}{)}_{i,j}^n=\frac{u_{i+1,j+1}^n-u_{i-1,j-1}^n-u_{i-1,j+1}^n-u_{i+1,j-1}}{4}$

最后贴出论文中的步骤：

其中h是空间步长，△t是时间步长

具体求解过程，论文中更详细，数学好的朋友可以去看看。变分图像去噪

python代码我会贴在文末，注意迭代次数不能太少，$\lambda$=0。

# -*- coding: utf-8 -*-
"""
Created on Sun Mar  7 14:24:12 2021

@author: SEELE
"""

import cv2
import numpy as np
import math

def gasuss_noise(image, mean=0, var=0.001):
  image = np.array(image/255, dtype=float)
  noise = np.random.normal(mean, var ** 0.5, image.shape)
  out = image + noise
  if out.min() < 0:
    low_clip = -1.
  else:
    low_clip = 0.
  out = np.clip(out, low_clip, 1.0)
  out = out*255
  #cv.imshow("gasuss", out)
  return out

img = cv2.imread('F:/hh/2.jpg')
img = cv2.resize(img,(500,500))
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
gray = gray.astype('float32')
size = img.shape  # h w c

#加高斯噪声
gray = gasuss_noise(gray)
cv2.imshow('addnoise',gray/255)

maxloop = 100
n = 0
#h = 1
t = 0.1
lamda = 0
ep = 1

divp = 0
imtmp = gray.copy()#曾经用等号赋值，吃了大亏
imnew = np.zeros([size[0],size[1]])
for loop in range(0,maxloop):
    print("processed %d times" % (loop))
    
    #这段也可以直接用矩阵计算来写
    for y in range(1,size[0] - 1):
        for x in range(1,size[1] - 1):
            
            
            ux = (imtmp[y,x+1] - imtmp[y,x-1]) / 2
            uy = (imtmp[y+1,x] - imtmp[y-1,x]) / 2
            uxx = imtmp[y,x+1] + imtmp[y,x-1] - imtmp[y,x] - imtmp[y,x]
            uyy = imtmp[y+1,x] + imtmp[y-1,x] - imtmp[y,x] - imtmp[y,x]
            
            uxy = (imtmp[y+1,x+1] + imtmp[y-1,x-1] - imtmp[y+1,x-1] - imtmp[y-1,x+1]) / 4
            
            divp = ((uy * uy + ep) * uxx - 2* ux * uy * uxy + (ux * ux + ep) * uyy) / ((ux * ux + uy * uy + ep))#防止除0

            imnew[y,x] = imtmp[y,x] + t *( lamda * (gray[y,x] - imtmp[y,x]) + divp)
            dv = 1* t * lamda * (gray[y,x] - imtmp[y,x])
            
            if x == 25 and y == 15:
                print('dp=%d,div=%d,im=%d' % (divp,dv,imnew[y,x]))
          
    if loop != maxloop-1:
        imtmp = imnew.copy()
    if loop % 20 == 0:
        string = 'loop='+str(loop)
        cv2.imshow(string,imnew.astype('uint8'))
    

#可以限制一下范围
#for y in range(1,size[0] - 1):
#    for x in range(1,size[1] - 1):
#        imnew[y,x] = min(max(imnew[y,x],0),255)
        
imnew = imnew.astype('uint8')
cv2.imshow('result',imnew)
            
            
            

图中分别是循环0、20、40、60、80、100次的结果，可以看出这个算法其实还是不能直接拿来用的