赞
踩
一、常规滤波
二、非局部均值滤波
三、维纳滤波
四、卡尔曼滤波
所谓滤波,其实就是从混合在一起的诸多信号中提取出所需要的信号。
信号的分类:
在图像处理或者计算机应用中,在正式对图像爱那个进行分析处理前一般需要一个预处理的过程。预处理是对图像作一些诸如降维、降噪的操作,主要是为后续处理提供一个体积合适、只包含所需信息的图像,通常会用到一些滤波处理手法。滤波,实际上就是信号处理,而图像本身可以看作是一个二维信号,其中像素点灰度的高低代表信号的强弱。对应的高低频的意义:
高频:图像中灰度变化强烈的点,一般是轮廓或者是噪声。
低频:图像中平坦的,灰度变化不大的点,图像中的大部各区域。
而根据图像的高频与低频的特征,可以设计相应的高通和低通滤波器,高通滤波可以检测图像中尖锐、变化明显的地方,而低通滤波可以让图像变得光滑,滤除图像中的噪声、OpenCv中提供的低通滤波有线性的均值滤波器、高斯滤波器,非线性的双边滤波器、中值滤波器;高通滤波有基于Canny,Sobel算子的各种滤波。其实很多时候低通滤波和高通滤波其实是相互矛盾的,很多时候在边缘检测前需要通过低通滤波降噪,这里就需要调节参数在保证高频的边缘不丢失的前提下尽可能多的去处图像的噪点。
这里使用频域的高通和低通滤波。
理想的低通滤波器的模版为:
其中,
- def low_pass_filter(img, radius=100):
- r = radius
-
- rows, cols = img.shape
- center = 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 * r
- mask[mask_area] = 1
- return mask
Butterworth低通滤波器为:
- def Butterworth(src, d0, n, ftype):
- template = np.zeros(src.shape, dtype=np.float32) # 构建滤波器
- r, c = src.shape
- for 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 - template
- return template
高斯低通滤波器:
- # 定义函数,高斯高/低通滤波模板
- def Gaussian(src, d0, ftype):
- template = np.zeros(src.shape, dtype=np.float32) # 构建滤波器
- r, c = src.shape
- for 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 - template
- return template
- def bandreject_filters(img, r_out=300, r_in=35):
- rows, cols = img.shape
- crow, ccol = int(rows / 2), int(cols / 2)
-
- radius_out = r_out
- radius_in = r_in
-
- mask = 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] = 1
- mask = 1 - mask
- return 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/2
- t = searchWindowSize/2
- height, width = I.shape[:2]
- padLength = t+f
- I2 = 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:
- continue
- I2_ = 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 += w
- wmax = 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.0
- I1 = 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) / 2
- print(center_position)
-
- slope_tan = math.tan(motion_angle * math.pi / 180)
- slope_cot = 1 / slope_tan
- if 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)] = 1
- return PSF / PSF.sum() # 对点扩散函数进行归一化亮度
- else:
- for i in range(15):
- offset = round(i * slope_cot)
- PSF[int(center_position - offset), int(center_position + offset)] = 1
- return PSF / PSF.sum()
-
- # 对图片进行运动模糊
- def make_blurred(input, PSF, eps):
- input_fft = fft.fft2(input) # 进行二维数组的傅里叶变换
- PSF_fft = fft.fft2(PSF) + eps
- blurred = fft.ifft2(input_fft * PSF_fft)
- blurred = np.abs(fft.fftshift(blurred))
- return blurred
-
- def inverse(input, PSF, eps): # 逆滤波
- input_fft = fft.fft2(input)
- PSF_fft = fft.fft2(PSF) + eps # 噪声功率,这是已知的,考虑epsilon
- result = fft.ifft2(input_fft / PSF_fft) # 计算F(u,v)的傅里叶反变换
- result = np.abs(fft.fftshift(result))
- return result
-
- def wiener(input, PSF, eps, K=0.01): # 维纳滤波,K=0.01
- input_fft = fft.fft2(input)
- PSF_fft = fft.fft2(PSF) + eps
- PSF_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 result
-
- def normal(array):
- array = np.where(array < 0, 0, array)
- array = np.where(array > 255, 255, array)
- array = array.astype(np.int16)
- return array
-
- def 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 channel
-
- if __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) / 2
- print(center_position)
-
- slope_tan = math.tan(motion_angle * math.pi / 180)
- slope_cot = 1 / slope_tan
- if 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)] = 1
- return PSF / PSF.sum() # 对点扩散函数进行归一化亮度
- else:
- for i in range(15):
- offset = round(i * slope_cot)
- PSF[int(center_position - offset), int(center_position + offset)] = 1
- return PSF / PSF.sum()
-
- # 对图片进行运动模糊
- def make_blurred(input, PSF, eps):
- input_fft = fft.fft2(input) # 进行二维数组的傅里叶变换
- PSF_fft = fft.fft2(PSF) + eps
- blurred = fft.ifft2(input_fft * PSF_fft)
- blurred = np.abs(fft.fftshift(blurred))
- return blurred
-
- def inverse(input, PSF, eps): # 逆滤波
- input_fft = fft.fft2(input)
- PSF_fft = fft.fft2(PSF) + eps # 噪声功率,这是已知的,考虑epsilon
- result = fft.ifft2(input_fft / PSF_fft) # 计算F(u,v)的傅里叶反变换
- result = np.abs(fft.fftshift(result))
- return result
-
- def wiener(input, PSF, eps, K=0.01): # 维纳滤波,K=0.01
- input_fft = fft.fft2(input)
- PSF_fft = fft.fft2(PSF) + eps
- PSF_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 result
-
- def normal(array):
- array = np.where(array < 0, 0, array)
- array = np.where(array > 255, 255, array)
- array = array.astype(np.int16)
- return array
-
- def 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 channel
-
- if __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. (最优化自回归数据数据算法)
这都是百度百科的东西,搞得跟数学书上公式定理一样,说的都不是人话!关键是字我都认识,但我没懂什么意思。
举个例子:在海图作业中,航海张通常以前一时刻的船位为基准,根据航向、船速和海流等一系列因素推算下一个船位,这个称之为观测船位;另外还会选择适当的方法,通过仪器得到另一个推算船位,这个称之为推算船位;观测船位和推算船位一般都不重合,航海长需要通过分析和判断选择一个可靠的船位,作为船舰当前的位置。
由此得出卡尔曼滤波思想:以
卡尔曼滤波算法在控制领域有着极广泛的应用,比如自动驾驶汽车,想象一下,一个雷达传感器告诉你另一辆车距离15米,一个激光传感器说车辆距离20米 ,那么如何协调这些传感器测量。再比如,在发动机燃油喷射控制中,可以应用扩展的卡尔曼滤波理论研究瞬态工况下发动机循环进气量的最有估计算法;在雷达中,人们最感兴趣的是跟踪目标,但目标的位置、速度、加速度的测量往往在任何时候都有噪声,而卡尔曼滤波则是利用目标的动态信息,设法去除噪声,得到一个关于目标位置的最好估计。
再举个例子:假设要研究一个房间的温度,以一分钟位时间单位,根据经验判断,这个房间的温度是恒定的,但是对于我们的经验不是完全相信,可能存在上下几度的偏差,我们把该偏差看作是高斯白噪声。另外,在房间里放一个温度计,而温度计也不准确,测量值会与实际值存在偏差,我们也把这偏差看作是高斯白噪声。那么现在,我们要根据经验温度和温度计的测量值以及他们各自的噪声来估算出放房间的实际温度。
那么接下来该如何解决呢? 假如我们要估算
所以估计K时刻最优温度值为
详情参照:[卡尔曼滤波器分类及其基本公式](卡尔曼滤波器分类及基本公式 - 百度文库)
就这样,卡尔曼滤波就能不断吧均方误差递归,从而估算出最优的温度值,运行速度快,且只保留上一时刻的协方差。
总而言之,Kalman滤波用在当测量值与模型预测值均不准确的情况下,用来计算预测真值的一种滤波算法,在目标识别与追踪任务中经常用到。
- def KalmanFilter(z, n_iter=20):
- # 这里是假设A=1,H=1的情况
-
- # intial parameters
-
- sz = (n_iter,) # size of array
-
- # Q = 1e-5 # process variance
- Q = 1e-6 # process variance
- # allocate space for arrays
- xhat = numpy.zeros(sz) # a posteri estimate of x
- P = numpy.zeros(sz) # a posteri error estimate
- xhatminus = numpy.zeros(sz) # a priori estimate of x
- Pminus = numpy.zeros(sz) # a priori error estimate
- K = numpy.zeros(sz) # gain or blending factor
-
- R = 0.1 ** 2 # estimate of measurement variance, change to see effect
-
- # intial guesses
- xhat[0] = 0.0
- P[0] = 1.0
- A = 1
- H = 1
-
- for k in range(1, n_iter):
- # time update
- xhatminus[k] = A * xhat[k - 1] # X(k|k-1) = AX(k-1|k-1) + BU(k) + W(k),A=1,BU(k) = 0
- Pminus[k] = A * P[k - 1] + Q # P(k|k-1) = AP(k-1|k-1)A' + Q(k) ,A=1
- # measurement update
- K[k] = Pminus[k] / (Pminus[k] + R) # Kg(k)=P(k|k-1)H'/[HP(k|k-1)H' + R],H=1
- xhat[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=1
- P[k] = (1 - K[k] * H) * Pminus[k] # P(k|k) = (1 - Kg(k)H)P(k|k-1), H=1
- return xhat
牛逼的算法往往都是来源一个很简单的思想所演化出来的!
继续坚持!加油~
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。