传统的图像去噪方法（四）之变换域去噪算法_变换域方法去噪

我们知道图像的经典去噪方法主要有两大类：一类是基于空间域的处理方法，一类是基于频域的处理方法。前面三节讲到的欧式基于空间域的处理方法，基于频域的处理方法主要是用滤波器，把有用的信号和干扰信号分开，它在有用信号和干扰信号的频谱没有重叠的前提下，才能把有用信号和干扰信号完全分开，但实际上很难分开。

因此图像变换域去噪算法的基本思想其实就是首先进行某种变换，将图像从空间域转换到变换域，然后从频率上把噪声分为高中低频噪声，用这种变换域的方法就可以把不同频率的噪声分离，之后进行反变换将图像从变换域转换到原始空间域，最终达到去除图像噪声的目的。

（一）傅里叶变换

傅立叶变换用于分析各种滤波器的频率特性，
对于一幅图像来说在分析其频率特性时，它的边缘，突出部分以及颗粒噪声往往代表图像信号的高频分量，而大面积的图像背景区则代表图像信号的低频分量。因此，使用滤波的方法滤除其高频部分也就可以去除噪声，使得图像得到一定的平滑。
图像是二维离散的，连续与离散都可以用傅里叶进行变换，那么二维信号无非就是在x方向与y方向都进行一次一维的傅里叶变换得到。所有图像的频率构成都认为是这样的，那么不同的就是一幅图的振幅与相位了，也就是说在opencv或者matlab下对图像进行傅里叶变换后其实是可以得到图像的振幅图与相位图的，而想把图像从频域空间恢复到时域空间，必须要同时有图像的振幅图与相位图才可以，缺少一个就恢复的不完整。

其中，F(u,v)是含噪声图像的傅里叶变换，G(u,v)是平滑后图像的傅里叶变换，H(u,v)是低通滤波器传递函数。处理过程如下图所示：

利用numoy包进行傅里叶变换

np.fft.fft2()
第一个参数是输入图像，它是灰度图像
第二个参数是可选的，它决定了输出数组的大小,如果它大于输入图像的大小，则输入图像在计算FFT之前填充了0.如果它小于输入图像，输入图像将被裁剪,如果没有参数传递，输出数组的大小将与输入相同.
一旦得到结果，零频率分量（DC分量）将位于左上角。 如果要将其置于中心位置，则需要在两个方向上将结果移动N2.np.fft.fftshift(),一旦你找到频率变换，你就能找到大小谱.

import cv2
import numpy as np
from matplotlib import pyplot as plt

img=cv2.imread('E:\opencv\yg.jpg',0)  # 直接读取灰度图像
f=np.fft.fft2(img)
fshift=np.fft.fftshift(f)
magnitude_spectrum=20*np.log(np.abs(fshift))  #取绝对值是为了将复数变换成实数，
                                              #取对数的目的是将数据变化到较小的范围（比如0--255）

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

但是上图并没有什么含义，显示出来的可以看成是频域中图像的振幅信息，但是没有相位信息，图像的相位ϕ=atan(实部/虚部），numpy包中自带一个angle函数可以直接根据复数的实部和虚部求出角度（默认是弧度）。

import cv2
import numpy as np
from matplotlib import pyplot as plt

img=cv2.imread('E:\opencv\yg.jpg',0)  # 直接读取灰度图像
f=np.fft.fft2(img)
fshift=np.fft.fftshift(f)

ph_f=np.angle(f)
ph_fshift=np.angle(fshift)

plt.subplot(121),plt.imshow(ph_f, cmap = 'gray')
plt.title('original'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(ph_fshift, cmap = 'gray')
plt.title('center'), plt.xticks([]), plt.yticks([])
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

上图就是图像上每个像素点对应的相位图，理解就是偏移的角度。

下面介绍一下频域下的低通滤波器，高通滤波器，带通带阻滤波器。

1.高通滤波器

图像在变换加移动中心后，从中间到外面，频率上依次是从低频到高频，因此我们如果把中间一小部分去掉，就相当于高通滤波器：
黑色为0，白色为1，把这个模板和图像的傅里叶变换的频域矩阵相乘就得到了高通滤波。

import cv2
import numpy as np
from matplotlib import pyplot as plt
img=cv2.imread('E:\opencv\yg.jpg',0)  # 直接读取灰度图像
f=np.fft.fft2(img)
fshift=np.fft.fftshift(f)
magnitude_spectrum=20*np.log(np.abs(fshift))  #取绝对值是为了将复数变换成实数，
                                              #取对数的目的是将数据变化到较小的范围（比如0--255）
rows,cols=img.shape
mask = np.ones(img.shape,img.dtype)
crow,ccol=int(rows/2),int(cols/2)
mask[crow-20:crow+20,ccol-20:ccol+20]=0  # 作高通滤波模板

fshift=mask*fshift # 滤波
f_ishift=np.fft.ifftshift(fshift)
img_back=np.fft.ifft2(f_ishift) # 傅里叶逆变换
img_back=np.abs(img_back) # 得到的是复数会无法显示

plt.subplot(221),plt.imshow(img, cmap = 'gray')
plt.title('original'), plt.xticks([]), plt.yticks([])
plt.subplot(222),plt.imshow(img_back, cmap = 'gray')
plt.title('Image after HPF'), plt.xticks([]), plt.yticks([])
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

可以看出，高通滤波器有利于提取图像的轮廓。这是因为图像的轮廓或者边缘或者噪声处，灰度变化强烈，那么在它们经过傅里叶变换后，就会变成高频信号。

2.低通滤波器

低通滤波器就是把上述模板的1变成0，0变成1：

import cv2
import numpy as np
from matplotlib import pyplot as plt

img=cv2.imread('E:\opencv\yg.jpg',0)  # 直接读取灰度图像
f=np.fft.fft2(img)
fshift=np.fft.fftshift(f)
magnitude_spectrum=20*np.log(np.abs(fshift))  #取绝对值是为了将复数变换成实数，
                                              #取对数的目的是将数据变化到较小的范围（比如0--255）

rows,cols=img.shape
mask = np.zeros(img.shape,img.dtype)
crow,ccol=int(rows/2),int(cols/2)
mask[crow-20:crow+20,ccol-20:ccol+20]=1  # 作低通滤波模板

fshift=mask*fshift # 滤波
f_ishift=np.fft.ifftshift(fshift)
img_back=np.fft.ifft2(f_ishift) # 傅里叶逆变换
img_back=np.abs(img_back) # 得到的是复数会无法显示

plt.subplot(221),plt.imshow(img, cmap = 'gray')
plt.title('original'), plt.xticks([]), plt.yticks([])
plt.subplot(222),plt.imshow(img_back, cmap = 'gray')
plt.title('Image after LPF'), plt.xticks([]), plt.yticks([])
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

可以看到低通滤波后除了轮廓模糊外，基本没变化，这是因为图像的主要信息都集中到了低频上。

3.带通滤波器

import cv2
import numpy as np
from matplotlib import pyplot as plt

img=cv2.imread('E:\opencv\yg.jpg',0)  # 直接读取灰度图像
f=np.fft.fft2(img)
fshift=np.fft.fftshift(f)

rows,cols=img.shape
mask1 = np.ones(img.shape,img.dtype)
mask2 = np.zeros(img.shape,img.dtype)
crow,ccol=int(rows/2),int(cols/2)
mask1[crow-8:crow+8,ccol-8:ccol+8]=0  # 作高通滤波
mask2[crow-80:crow+80,ccol-80:ccol+80]=1  # 作低通滤波模板
mask=mask1*mask2  # 带通滤波

fshift=mask*fshift # 滤波
f_ishift=np.fft.ifftshift(fshift)
img_back=np.fft.ifft2(f_ishift) # 傅里叶逆变换
img_back=np.abs(img_back) # 得到的是复数会无法显示

plt.subplot(221),plt.imshow(img, cmap = 'gray')
plt.title('original'), plt.xticks([]), plt.yticks([])
plt.subplot(222),plt.imshow(img_back, cmap = 'gray')
plt.title('Image after BPF'), plt.xticks([]), plt.yticks([])
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

带通滤波的效果，既保留了一部分低频，也保留了一部分高频。

（二）基于小波变换的图像去噪技术

傅里叶变换对于非平稳过程有局限性：它只能获取一段信号总体上包含哪些频率的成分，但是对各个成分出现的时刻不清楚。因此时域相差很大的两个信号，可能频谱图一样。
但是对于一个非平稳信号，只知道包含哪些频率成分时不够的，我们还想知道各个成分出现的时间，知道信号频率随时间变化的情况，各个时刻的顺时频率及其幅值----时频分析。
小波把傅里叶变换的基换了，将无限长的三角函数基换成了有限长的会衰减的小波基。这样不仅能够获取频率，还可以定位到时间了。

1.小波阈值去噪介绍

主要思想是经过小波变换后图像和噪声的统计特性不同，其中图像本身的小波系数具有较大幅值，主要集中在高频，噪声小波系数幅值较小，并且存在于小波变换后的所有系数中。因此设置一个阈值门限来进行分离。
通常在去噪时，不处理含有大量图像能量的低通系数，只是就单个高通部分进行处理。因此不能只进行一次阈值去噪，还需要对低频部分进行阈值去噪和小波分解，直到估计噪声和实际图像的偏差值最小。一般进行3-4层小波分解和去噪就可以了。

2.小波阈值去噪方法
算法的基本过程为：
①对原始信号进行小波分解
②对变换后的小波系数进行阈值处理，得到估计小波系数
③根据估计小波系数进行小波重构

3.小波阈值去噪基本问题

1. 小波基的选择：通常我们希望所选取的小波满足以下条件：正交性、高消失矩、紧支性、对称性或反对称性。但事实上具有上述性质的小波是不可能存在的，因为小波是对称或反对称的只有Haar小波，并且高消失矩与紧支性是一对矛盾，所以在应用的时候一般选取具有紧支的小波以及根据信号的特征来选取较为合适的小波。
2. 阈值的选择：直接影响去噪效果的一个重要因素就是阀值的选取，不同的阈值选取将有不同的去噪效果。目前主要有通用阈值（VisuShrink）、SureShrink阈值、Minimax阈值、BayesShrink阈值等。
3. 阈值函数的选择：阈值函数是修正小波系数的规则，不同的反之函数体现了不同的处理小波系数的策略。最常用的阈值函数有两种：一种是硬阈值函数，另一种是软阈值函数。还有一种介于软、硬阈值函数之间的Garrote函数。

利用阈值函数pywt.threshold:

import pywt
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math

# np.linspace()函数用于在指定的间隔内返回均匀间隔的数字。
data=np.linspace(1,10,10)
#[1. 2. 3. 4. 5. 6. 7. 8. 9. 10.]

# pywt.threshold(data,value,mode,substitute) mode 模式有四种，soft,hard,greater,less;substitute是替换值

data_soft=pywt.threshold(data=data,value=6,mode='soft',substitute=12)
# soft模式把小于6的值设置为12，大于等于6的值全部减6
#[12. 12. 12. 12. 12. 0. 1. 2. 3. 4.]

data_hard=pywt.threshold(data=data,value=6,mode='hard',substitute=12)
# hard模式把小于6的值设置为12，其余的值不变
#[12. 12. 12. 12. 12. 6. 7. 8. 9. 10.]

data_greater=pywt.threshold(data,6,'greater',12)
# 把小于6的值设为12，大于等于阈值的值不变
#[12. 12. 12. 12. 12. 6. 7. 8. 9. 10.]

data_less=pywt.threshold(data,6,'less',12)
# 把大于6的值设为12，小于等于阈值的值不变
#[1. 2. 3. 4. 5. 6. 12. 12. 12. 12. 12. ]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

4.小波去噪实验

pywt.dwt_max_level(data_len,filter_len):
pywt.wavedec(data, wavelet, mode=‘symmetric’, level=None, axis=-1)
data:输入数据，wavelet：使用的小波，level：分解级别，默认为dwt_max_len
返回：[cA_n, cD_n, cD_n-1, …, cD2, cD1]
n代表分解的级别，第一个cA_n代表的是近似分解系数，接下来的元素cD_n-cD1代表的是系数数组的详细信息

import pywt
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math

data=np.linspace(1,10,10)

# Get Data
ecg=pywt.data.ecg() # 生成心电信号
index=[]
data=[]
for i in range(len(ecg)-1):
    X=float(i)
    Y=float(ecg[i])
    index.append(X)
    data.append(Y)
    
# Create wavelet object and define parameters 
w=pywt.Wavelet('db8') # 创建一个小波对象 选用Daubechies8小波
#  Family name:    Daubechies
#  Short name:     db
#  Filters length: 16  #滤波器长度为16
#  Orthogonal:     True  #正交
#  Biorthogonal:   True  #双正交
#  Symmetry:       asymmetric  # 对称性：不对称
#  DWT:            True  #离散小波变换
#  CWT:            False
maxlev=pywt.dwt_max_level(len(data),w.dec_len)  # 计算最大有用的分解级别（data_len,filter_len)
print("maximum level is " + str(maxlev))
threshold=0.04 # Threshold for filtering

# Decompose into wavelet components,to the level selected
coeffs=pywt.wavedec(data,'db8',level=maxlev)

plt.figure()
for i in range(1,len(coeffs)):
    coeffs[i]=pywt.threshold(coeffs[i],threshold*max(coeffs[i]))  # 将噪声滤波

datarec=pywt.waverec(coeffs,'db8')  # 将信号进行小波重构

mintime=0
maxtime=mintime+len(data)+1

plt.figure()
plt.subplot(2,1,1)
plt.plot(index[mintime:maxtime],data[mintime:maxtime])
plt.xlabel('time(s)')
plt.ylabel('microvolts (uV)')
plt.title("Raw signal")
plt.subplot(2,1,2)
plt.plot(index[mintime:maxtime],datarec[mintime:maxtime-1])
plt.xlabel('time(s)')
plt.ylabel('microvolts (uV)')
plt.title("De-noised signal using wavelet techniques")

plt.tight_layout()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

参考链接：
Python小波变换去噪
图像滤波和傅里叶变换