图像理想低通滤波_图像处理之滤波（下）

[toc]目录

一、常规滤波

低通
高通
带通
带阻

二、非局部均值滤波

三、维纳滤波

四、卡尔曼滤波

前言

所谓滤波，其实就是从混合在一起的诸多信号中提取出所需要的信号。

信号的分类：

确定型信号，可以表示为确定的时间函数，可确定其在任何时刻的量值。（具有确定的频谱）；一般可通过低通、高通、带通、带阻等模拟滤波器或其他常规滤波算法实现。
随机信号，不能用确定的数学关系式描述，不能预测其未来任何瞬时值，其值的变化服从统计规律。（频谱不确定，功率谱确定）；根据有用信号和干扰信号的功率谱设计滤波器——维纳滤波（Wiener Filtering）或卡尔曼滤波（Kalman Filter）。

一、常规滤波

在图像处理或者计算机应用中，在正式对图像爱那个进行分析处理前一般需要一个预处理的过程。预处理是对图像作一些诸如降维、降噪的操作，主要是为后续处理提供一个体积合适、只包含所需信息的图像，通常会用到一些滤波处理手法。滤波，实际上就是信号处理，而图像本身可以看作是一个二维信号，其中像素点灰度的高低代表信号的强弱。对应的高低频的意义：

高频：图像中灰度变化强烈的点，一般是轮廓或者是噪声。

低频：图像中平坦的，灰度变化不大的点，图像中的大部各区域。

而根据图像的高频与低频的特征，可以设计相应的高通和低通滤波器，高通滤波可以检测图像中尖锐、变化明显的地方，而低通滤波可以让图像变得光滑，滤除图像中的噪声、OpenCv中提供的低通滤波有线性的均值滤波器、高斯滤波器，非线性的双边滤波器、中值滤波器；高通滤波有基于Canny，Sobel算子的各种滤波。其实很多时候低通滤波和高通滤波其实是相互矛盾的，很多时候在边缘检测前需要通过低通滤波降噪，这里就需要调节参数在保证高频的边缘不丢失的前提下尽可能多的去处图像的噪点。

这里使用频域的高通和低通滤波。

低通

理想的低通滤波器的模版为：

$H(u,v) = \begin{equation} \left\{ \begin{array}{lr} 1, \quad D(u, v) \le D_0 \\ 0, \quad D(u, v) \ge D_0 \end{array} \right. \end{equation}$

其中，

$D_0$ 表示通带半径，

$D(u, v)$ 是到频谱中心的距离（欧式距离），计算公式如下：

$D(u, v) = \sqrt {(u - M/ 2)^2 + (v-N/2)^2}$ ，M和N表示频谱图像的大小，

$(M/2, N/2)$ 即为频谱中心。

def low_pass_filter(img, radius=100):r = radiusrows, cols = img.shapecenter = int(rows / 2), int(cols / 2)mask = np.zeros((rows, cols, 2), np.uint8)x, y = np.ogrid[:rows, :cols]mask_area = (x - center[0]) ** 2 + (y - center[1]) ** 2 <= r * rmask[mask_area] = 1return mask

Butterworth低通滤波器为：

$H(u,v) = \frac{1}{1+[\frac{D(u,v)}{D_0}]^{2n}}$

def Butterworth(src, d0, n, ftype):template = np.zeros(src.shape, dtype=np.float32)  # 构建滤波器r, c = src.shapefor i in np.arange(r):for j in np.arange(c):distance = np.sqrt((i - r/2)**2 + (j - c/2)**2)template[i, j] = 1/(1 + (distance/d0)**(2*n))  # Butterworth 滤波函数template[i, j] = np.e ** (-1 * (distance**2 / (2 * d0**2)))  # Gaussian滤波函数if ftype == 'high':template = 1 - templatereturn template

高斯低通滤波器：

$H(u,v) = e^{\frac{-D^2(u, v)}{2D_0^2}}$ ,1减去低通滤波器模板即可得到高通滤波器。

# 定义函数，高斯高/低通滤波模板
def Gaussian(src, d0, ftype):template = np.zeros(src.shape, dtype=np.float32)  # 构建滤波器r, c = src.shapefor i in np.arange(r):for j in np.arange(c):distance = np.sqrt((i - r / 2) ** 2 + (j - c / 2) ** 2)template[i, j] = np.e ** (-1 * (distance ** 2 / (2 * d0 ** 2)))  # Gaussian滤波函数if ftype == 'high':template = 1 - templatereturn template

带通
带阻

def bandreject_filters(img, r_out=300, r_in=35):rows, cols = img.shapecrow, ccol = int(rows / 2), int(cols / 2)radius_out = r_outradius_in = r_inmask = np.zeros((rows, cols, 2), np.uint8)center = [crow, ccol]x, y = np.ogrid[:rows, :cols]mask_area = np.logical_and(((x - center[0]) ** 2 + (y - center[1]) ** 2 >= r_in ** 2),((x - center[0]) ** 2 + (y - center[1]) ** 2 <= r_out ** 2))mask[mask_area] = 1mask = 1 - maskreturn mask

二、非局部均值滤波

图像中的像素点之间不是孤立存在的，某一点的像素与别处的像素点一定存在着某种关联，可以概括为灰度相关性和几何结构相似性。

#coding:utf8
import cv2
import numpy as np
def psnr(A, B):return 10*np.log(255*255.0/(((A.astype(np.float)-B)**2).mean()))/np.log(10)def double2uint8(I, ratio=1.0):return np.clip(np.round(I*ratio), 0, 255).astype(np.uint8)def make_kernel(f):kernel = np.zeros((2*f+1, 2*f+1))for d in range(1, f+1):kernel[f-d:f+d+1, f-d:f+d+1] += (1.0/((2*d+1)**2))return kernel/kernel.sum()def NLmeansfilter(I, h_=10, templateWindowSize=5,  searchWindowSize=11):f = templateWindowSize/2t = searchWindowSize/2height, width = I.shape[:2]padLength = t+fI2 = np.pad(I, padLength, 'symmetric')kernel = make_kernel(f)h = (h_**2)I_ = I2[padLength-f:padLength+f+height, padLength-f:padLength+f+width]average = np.zeros(I.shape)sweight = np.zeros(I.shape)wmax =  np.zeros(I.shape)for i in range(-t, t+1):for j in range(-t, t+1):if i==0 and j==0:continueI2_ = I2[padLength+i-f:padLength+i+f+height, padLength+j-f:padLength+j+f+width]w = np.exp(-cv2.filter2D((I2_ - I_)**2, -1, kernel)/h)[f:f+height, f:f+width]sweight += wwmax = np.maximum(wmax, w)average += (w*I2_[f:f+height, f:f+width])return (average+wmax*I)/(sweight+wmax)if __name__ == '__main__':I = cv2.imread('lena.jpg', 0)sigma = 20.0I1 = double2uint8(I + np.random.randn(*I.shape) *sigma)print u'噪声图像PSNR',psnr(I, I1)R1  = cv2.medianBlur(I1, 5)print u'中值滤波PSNR',psnr(I, R1)R2 = cv2.fastNlMeansDenoising(I1, None, sigma, 5, 11)print u'opencv的NLM算法',psnr(I, R2)R3 = double2uint8(NLmeansfilter(I1.astype(np.float), sigma, 5, 11))print u'NLM PSNR',psnr(I, R3)

尽管opencv中已经有实现，对于彩色图像，首先要先转换到CIELAB颜色空间，然后对L和AB成分分别去噪。而且据说上面的实现会比opencv自带的实现要好一些。

cv2.fastNlMeansDenoising() - 使用单个灰度图像
cv2.fastNlMeansDenoisingColored() - 使用彩色图像。
cv2.fastNlMeansDenoisingMulti() - 用于在短时间内捕获的图像序列（灰度图像）
cv2.fastNlMeansDenoisingColoredMulti() - 与上面相同，但用于彩色图像。fastNlMeansDenoisingColored( InputArray src, OutputArray dst,float h = 3, float hColor = 3,int templateWindowSize = 7, int searchWindowSize = 21)
参数：
• h : 决定过滤器强度。h 值高可以很好的去除噪声,但也会把图像的细节抹去。(取 10 的效果不错)
• hForColorComponents : 与 h 相同,但使用与彩色图像。(与 h 相同,10)
• templateWindowSize : 奇数。(推荐值为 7)
• searchWindowSize : 奇数。(推荐值为 21)

三、维纳滤波

对于运动引起的模糊，最简单的方法就是直接作逆滤波，但是逆滤波对于加性噪声特别敏感，使得恢复的图像几乎不可用。最小均方差（维纳）滤波用来去处含有噪声的模糊图像，其目标是找到未污染图像的一个估计，使得他们之间的均方误差最小，可以去除噪声，同时清晰化模糊图像。

主要参考：写的不错，有公式又证明！里面还有约束最小二乘方滤波！

https://blog.csdn.net/wsp_1138886114/article/details/95024180blog.csdn.net

import matplotlib.pyplot as plt
import numpy as np
from numpy import fft
import math
import cv2# 仿真运动模糊
def motion_process(image_size, motion_angle):PSF = np.zeros(image_size)print(image_size)center_position = (image_size[0] - 1) / 2print(center_position)slope_tan = math.tan(motion_angle * math.pi / 180)slope_cot = 1 / slope_tanif slope_tan <= 1:for i in range(15):offset = round(i * slope_tan)  # ((center_position-i)*slope_tan)PSF[int(center_position + offset), int(center_position - offset)] = 1return PSF / PSF.sum()             # 对点扩散函数进行归一化亮度else:for i in range(15):offset = round(i * slope_cot)PSF[int(center_position - offset), int(center_position + offset)] = 1return PSF / PSF.sum()# 对图片进行运动模糊
def make_blurred(input, PSF, eps):input_fft = fft.fft2(input)             # 进行二维数组的傅里叶变换PSF_fft = fft.fft2(PSF) + epsblurred = fft.ifft2(input_fft * PSF_fft)blurred = np.abs(fft.fftshift(blurred))return blurreddef inverse(input, PSF, eps):                # 逆滤波input_fft = fft.fft2(input)PSF_fft = fft.fft2(PSF) + eps            # 噪声功率，这是已知的，考虑epsilonresult = fft.ifft2(input_fft / PSF_fft)  # 计算F(u,v)的傅里叶反变换result = np.abs(fft.fftshift(result))return resultdef wiener(input, PSF, eps, K=0.01):        # 维纳滤波，K=0.01input_fft = fft.fft2(input)PSF_fft = fft.fft2(PSF) + epsPSF_fft_1 = np.conj(PSF_fft) / (np.abs(PSF_fft) ** 2 + K)result = fft.ifft2(input_fft * PSF_fft_1)result = np.abs(fft.fftshift(result))return resultdef normal(array):array = np.where(array < 0,  0, array)array = np.where(array > 255, 255, array)array = array.astype(np.int16)return arraydef main(gray):channel = []img_h, img_w = gray.shape[:2]PSF = motion_process((img_h, img_w), 60)      # 进行运动模糊处理blurred = np.abs(make_blurred(gray, PSF, 1e-3))result_blurred = inverse(blurred, PSF, 1e-3)  # 逆滤波result_wiener = wiener(blurred, PSF, 1e-3)    # 维纳滤波blurred_noisy = blurred + 0.1 * blurred.std() * np.random.standard_normal(blurred.shape)  # 添加噪声,standard_normal产生随机的函数inverse_mo2no = inverse(blurred_noisy, PSF, 0.1 + 1e-3)   # 对添加噪声的图像进行逆滤波wiener_mo2no = wiener(blurred_noisy, PSF, 0.1 + 1e-3)     # 对添加噪声的图像进行维纳滤波channel.append((normal(blurred),normal(result_blurred),normal(result_wiener),normal(blurred_noisy),normal(inverse_mo2no),normal(wiener_mo2no)))return channelif __name__ == '__main__':image = cv2.imread('./gggg/001.png')b_gray, g_gray, r_gray = cv2.split(image.copy())Result = []for gray in [b_gray, g_gray, r_gray]:channel = main(gray)Result.append(channel)blurred = cv2.merge([Result[0][0][0], Result[1][0][0], Result[2][0][0]])result_blurred = cv2.merge([Result[0][0][1], Result[1][0][1], Result[2][0][1]])result_wiener = cv2.merge([Result[0][0][2], Result[1][0][2], Result[2][0][2]])blurred_noisy = cv2.merge([Result[0][0][3], Result[1][0][3], Result[2][0][3]])inverse_mo2no = cv2.merge([Result[0][0][4], Result[1][0][4], Result[2][0][4]])wiener_mo2no = cv2.merge([Result[0][0][5],  Result[1][0][5], Result[2][0][5]])#========= 可视化 ==========plt.figure(1)plt.xlabel("Original Image")plt.imshow(np.flip(image, axis=2))                         # 显示原图像plt.figure(2)plt.figure(figsize=(8, 6.5))imgNames = {"Motion blurred":blurred,"inverse deblurred":result_blurred,"wiener deblurred(k=0.01)":result_wiener,"motion & noisy blurred":blurred_noisy,"inverse_mo2no":inverse_mo2no,'wiener_mo2no':wiener_mo2no}for i,(key,imgName) in enumerate(imgNames.items()):plt.subplot(231+i)plt.xlabel(key)plt.imshow(np.flip(imgName, axis=2))plt.show()
OpenCV-Python 图像去模糊（维纳滤波，约束最小二乘方滤波）import matplotlib.pyplot as plt
import numpy as np
from numpy import fft
import math
import cv2# 仿真运动模糊
def motion_process(image_size, motion_angle):PSF = np.zeros(image_size)print(image_size)center_position = (image_size[0] - 1) / 2print(center_position)slope_tan = math.tan(motion_angle * math.pi / 180)slope_cot = 1 / slope_tanif slope_tan <= 1:for i in range(15):offset = round(i * slope_tan)  # ((center_position-i)*slope_tan)PSF[int(center_position + offset), int(center_position - offset)] = 1return PSF / PSF.sum()             # 对点扩散函数进行归一化亮度else:for i in range(15):offset = round(i * slope_cot)PSF[int(center_position - offset), int(center_position + offset)] = 1return PSF / PSF.sum()# 对图片进行运动模糊
def make_blurred(input, PSF, eps):input_fft = fft.fft2(input)             # 进行二维数组的傅里叶变换PSF_fft = fft.fft2(PSF) + epsblurred = fft.ifft2(input_fft * PSF_fft)blurred = np.abs(fft.fftshift(blurred))return blurreddef inverse(input, PSF, eps):                # 逆滤波input_fft = fft.fft2(input)PSF_fft = fft.fft2(PSF) + eps            # 噪声功率，这是已知的，考虑epsilonresult = fft.ifft2(input_fft / PSF_fft)  # 计算F(u,v)的傅里叶反变换result = np.abs(fft.fftshift(result))return resultdef wiener(input, PSF, eps, K=0.01):        # 维纳滤波，K=0.01input_fft = fft.fft2(input)PSF_fft = fft.fft2(PSF) + epsPSF_fft_1 = np.conj(PSF_fft) / (np.abs(PSF_fft) ** 2 + K)result = fft.ifft2(input_fft * PSF_fft_1)result = np.abs(fft.fftshift(result))return resultdef normal(array):array = np.where(array < 0,  0, array)array = np.where(array > 255, 255, array)array = array.astype(np.int16)return arraydef main(gray):channel = []img_h, img_w = gray.shape[:2]PSF = motion_process((img_h, img_w), 60)      # 进行运动模糊处理blurred = np.abs(make_blurred(gray, PSF, 1e-3))result_blurred = inverse(blurred, PSF, 1e-3)  # 逆滤波result_wiener = wiener(blurred, PSF, 1e-3)    # 维纳滤波blurred_noisy = blurred + 0.1 * blurred.std() * np.random.standard_normal(blurred.shape)  # 添加噪声,standard_normal产生随机的函数inverse_mo2no = inverse(blurred_noisy, PSF, 0.1 + 1e-3)   # 对添加噪声的图像进行逆滤波wiener_mo2no = wiener(blurred_noisy, PSF, 0.1 + 1e-3)     # 对添加噪声的图像进行维纳滤波channel.append((normal(blurred),normal(result_blurred),normal(result_wiener),normal(blurred_noisy),normal(inverse_mo2no),normal(wiener_mo2no)))return channelif __name__ == '__main__':image = cv2.imread('./gggg/001.png')b_gray, g_gray, r_gray = cv2.split(image.copy())Result = []for gray in [b_gray, g_gray, r_gray]:channel = main(gray)Result.append(channel)blurred = cv2.merge([Result[0][0][0], Result[1][0][0], Result[2][0][0]])result_blurred = cv2.merge([Result[0][0][1], Result[1][0][1], Result[2][0][1]])result_wiener = cv2.merge([Result[0][0][2], Result[1][0][2], Result[2][0][2]])blurred_noisy = cv2.merge([Result[0][0][3], Result[1][0][3], Result[2][0][3]])inverse_mo2no = cv2.merge([Result[0][0][4], Result[1][0][4], Result[2][0][4]])wiener_mo2no = cv2.merge([Result[0][0][5],  Result[1][0][5], Result[2][0][5]])#========= 可视化 ==========plt.figure(1)plt.xlabel("Original Image")plt.imshow(np.flip(image, axis=2))                         # 显示原图像plt.figure(2)plt.figure(figsize=(8, 6.5))imgNames = {"Motion blurred":blurred,"inverse deblurred":result_blurred,"wiener deblurred(k=0.01)":result_wiener,"motion & noisy blurred":blurred_noisy,"inverse_mo2no":inverse_mo2no,'wiener_mo2no':wiener_mo2no}for i,(key,imgName) in enumerate(imgNames.items()):plt.subplot(231+i)plt.xlabel(key)plt.imshow(np.flip(imgName, axis=2))plt.show()

四、卡尔曼滤波

Kalman filtering是一种利用线性系统状态方程，通过系统输入输出观测数据，对系统状态进行最优估计的算法。Kalman滤波在测量方差已知的情况下能够从一系列存在测量噪声的数据中，估计动态系统的状态。A Kalman filter is an optimal estimation algorithm. （最优化自回归数据数据算法）

这都是百度百科的东西，搞得跟数学书上公式定理一样，说的都不是人话！关键是字我都认识，但我没懂什么意思。

1.1 是什么？

举个例子：在海图作业中，航海张通常以前一时刻的船位为基准，根据航向、船速和海流等一系列因素推算下一个船位，这个称之为观测船位；另外还会选择适当的方法，通过仪器得到另一个推算船位，这个称之为推算船位；观测船位和推算船位一般都不重合，航海长需要通过分析和判断选择一个可靠的船位，作为船舰当前的位置。

由此得出卡尔曼滤波思想：以

$K-1$ 时刻的最有估计

$x_{k-1}$ 为准，预测

$K$ 时刻的状态变量

$\hat x_{k/k-1}$ ，同时又对该状态进行观测，得到观测变量

$z_k$ ，再在预测与预测之间进行分析，或者说是以观测量对预测量进行修正，从而得到

$K$ 时刻的最有状态估计

$x_k$ 。

卡尔曼滤波算法在控制领域有着极广泛的应用，比如自动驾驶汽车，想象一下，一个雷达传感器告诉你另一辆车距离15米，一个激光传感器说车辆距离20米，那么如何协调这些传感器测量。再比如，在发动机燃油喷射控制中，可以应用扩展的卡尔曼滤波理论研究瞬态工况下发动机循环进气量的最有估计算法；在雷达中，人们最感兴趣的是跟踪目标，但目标的位置、速度、加速度的测量往往在任何时候都有噪声，而卡尔曼滤波则是利用目标的动态信息，设法去除噪声，得到一个关于目标位置的最好估计。

1.2 怎么回事？

再举个例子：假设要研究一个房间的温度，以一分钟位时间单位，根据经验判断，这个房间的温度是恒定的，但是对于我们的经验不是完全相信，可能存在上下几度的偏差，我们把该偏差看作是高斯白噪声。另外，在房间里放一个温度计，而温度计也不准确，测量值会与实际值存在偏差，我们也把这偏差看作是高斯白噪声。那么现在，我们要根据经验温度和温度计的测量值以及他们各自的噪声来估算出放房间的实际温度。

那么接下来该如何解决呢？假如我们要估算

$k$ 时刻的实际温度，首先要根据

$k-1$ 时刻的温度值，来预测

$k$ 时刻的温度（

$K$ 时刻的经验温度）。因为我们的经验认为温度是恒定的，所以会得到

$k$ 时刻的温度和

$k-1$ 时刻是一样的，假设是23度，同时该值（预测值）的高斯白噪声为5度，（5是这样得到的。如果k-1时刻估算出的最优温度值的偏差是3，对自己预测的不确定度是4度，他们平方相加再开方，就是5），而温度计得到的温度值为25度，同时该值的偏差为4，此时，对于K时刻房间的温度值有两个：估计值23度和测量值25度，那么究竟相信谁？用均方根误差判断，

$H^2 = \frac{5^2}{5^2 + 4^2} \Rightarrow H = 0.78$

所以估计K时刻最优温度值为

$23+0.78*(25-23)=24.56$ 度，得到了K时刻的最优温度，下一步就是对K+1时刻的温度值进行最优估算，需要K时刻的最优温度的偏差，

$\sqrt {(1-H) *5^2} = 2.35$

无控制离散型卡尔曼滤波基本方程

带有控制的离散型卡尔曼滤波基本方程

详情参照：[卡尔曼滤波器分类及其基本公式](卡尔曼滤波器分类及基本公式 - 百度文库)

就这样，卡尔曼滤波就能不断吧均方误差递归，从而估算出最优的温度值，运行速度快，且只保留上一时刻的协方差。

总而言之，Kalman滤波用在当测量值与模型预测值均不准确的情况下，用来计算预测真值的一种滤波算法，在目标识别与追踪任务中经常用到。

1.3 python代码实现

def KalmanFilter(z, n_iter=20):# 这里是假设A=1，H=1的情况# intial parameterssz = (n_iter,)  # size of array# Q = 1e-5 # process varianceQ = 1e-6  # process variance# allocate space for arraysxhat = numpy.zeros(sz)  # a posteri estimate of xP = numpy.zeros(sz)  # a posteri error estimatexhatminus = numpy.zeros(sz)  # a priori estimate of xPminus = numpy.zeros(sz)  # a priori error estimateK = numpy.zeros(sz)  # gain or blending factorR = 0.1 ** 2  # estimate of measurement variance, change to see effect# intial guessesxhat[0] = 0.0P[0] = 1.0A = 1H = 1for k in range(1, n_iter):# time updatexhatminus[k] = A * xhat[k - 1]  # X(k|k-1) = AX(k-1|k-1) + BU(k) + W(k),A=1,BU(k) = 0Pminus[k] = A * P[k - 1] + Q  # P(k|k-1) = AP(k-1|k-1)A' + Q(k) ,A=1# measurement updateK[k] = Pminus[k] / (Pminus[k] + R)  # Kg(k)=P(k|k-1)H'/[HP(k|k-1)H' + R],H=1xhat[k] = xhatminus[k] + K[k] * (z[k] - H * xhatminus[k])  # X(k|k) = X(k|k-1) + Kg(k)[Z(k) - HX(k|k-1)], H=1P[k] = (1 - K[k] * H) * Pminus[k]  # P(k|k) = (1 - Kg(k)H)P(k|k-1), H=1return xhat

每日一句毒鸡汤：

牛逼的算法往往都是来源一个很简单的思想所演化出来的！

继续坚持！加油~