当前位置:   article > 正文

python低通滤波函数_python数字图像处理(四) 频率域滤波

python numpy 滤波函数

import matplotlib.pyplot as plt

import numpy as np

import cv2

%matplotlib inline

首先读入这次需要使用的图像

img = cv2.imread('apple.jpg',0) #直接读为灰度图像

plt.imshow(img,cmap="gray")

plt.axis("off")

plt.show()

0069Secoly1fm0gna2wbhj307a070ab7.jpg

使用numpy带的fft库完成从频率域到空间域的转换。

f = np.fft.fft2(img)

fshift = np.fft.fftshift(f)

低通滤波器

低通滤波器的公式如下

\[H(u,v)=

\begin{cases}

1, & \text{if $D(u,v)$ } \leq D_{0}\\

0, & \text{if $D(u,v)$} \geq D_{0}

\end{cases}

\]

其中D(u,v)为频率域上(u,v)点到中心的距离,D0由自己设置

SouthEast

白点就是所允许通过的频率范围

3d图像如下

SouthEast

我们先把苹果转化成频率域看下效果

#取绝对值:将复数变化成实数

#取对数的目的为了将数据变化到0-255

s1 = np.log(np.abs(fshift))

plt.subplot(121),plt.imshow(s1,'gray')

plt.title('Frequency Domain')

plt.show()

0069Secoly1fm0gnw1oomj305g05jt9f.jpg

matplotlib对于不是uint8的图像会自动把图像的数值缩放到0-255上,更多可以查看对该问题的讨论

我们在频率域上试着取不同的d0再将其反变换到空间域看下效果

def make_transform_matrix(d,image):

transfor_matrix = np.zeros(image.shape)

center_point = tuple(map(lambda x:(x-1)/2,s1.shape))

for i in range(transfor_matrix.shape[0]):

for j in range(transfor_matrix.shape[1]):

def cal_distance(pa,pb):

from math import sqrt

dis = sqrt((pa[0]-pb[0])**2+(pa[1]-pb[1])**2)

return dis

dis = cal_distance(center_point,(i,j))

if dis <= d:

transfor_matrix[i,j]=1

else:

transfor_matrix[i,j]=0

return transfor_matrix

d_1 = make_transform_matrix(10,fshift)

d_2 = make_transform_matrix(30,fshift)

d_3 = make_transform_matrix(50,fshift)

设定距离分别为10,30,50其通过的频率的范围如图

plt.subplot(131)

plt.axis("off")

plt.imshow(d_1,cmap="gray")

plt.title('D_1 10')

plt.subplot(132)

plt.axis("off")

plt.title('D_2 30')

plt.imshow(d_2,cmap="gray")

plt.subplot(133)

plt.axis("off")

plt.title("D_3 50")

plt.imshow(d_3,cmap="gray")

plt.show()

0069Secoly1fm0gnw0tdzj30ai041q2p.jpg

img_d1 = np.abs(np.fft.ifft2(np.fft.ifftshift(fshift*d_1)))

img_d2 = np.abs(np.fft.ifft2(np.fft.ifftshift(fshift*d_2)))

img_d3 = np.abs(np.fft.ifft2(np.fft.ifftshift(fshift*d_3)))

plt.subplot(131)

plt.axis("off")

plt.imshow(img_d1,cmap="gray")

plt.title('D_1 10')

plt.subplot(132)

plt.axis("off")

plt.title('D_2 30')

plt.imshow(img_d2,cmap="gray")

plt.subplot(133)

plt.axis("off")

plt.title("D_3 50")

plt.imshow(img_d3,cmap="gray")

plt.show()

0069Secoly1fm0gnw1n54j30ai041mxm.jpg

讲上面过程整理得到频率域低通滤波器的代码如下

def lowPassFilter(image,d):

f = np.fft.fft2(image)

fshift = np.fft.fftshift(f)

def make_transform_matrix(d):

transfor_matrix = np.zeros(image.shape)

center_point = tuple(map(lambda x:(x-1)/2,s1.shape))

for i in range(transfor_matrix.shape[0]):

for j in range(transfor_matrix.shape[1]):

def cal_distance(pa,pb):

from math import sqrt

dis = sqrt((pa[0]-pb[0])**2+(pa[1]-pb[1])**2)

return dis

dis = cal_distance(center_point,(i,j))

if dis <= d:

transfor_matrix[i,j]=1

else:

transfor_matrix[i,j]=0

return transfor_matrix

d_matrix = make_transform_matrix(d)

new_img = np.abs(np.fft.ifft2(np.fft.ifftshift(fshift*d_matrix)))

return new_img

plt.imshow(lowPassFilter(img,60),cmap="gray")

0069Secoly1fm0gnw1w86j307a070q3o.jpg

高通滤波器

高通滤波器同低通滤波器非常类似,只不过二者通过的波正好是相反的

\[H(u,v)=

\begin{cases}

0, & \text{if $D(u,v)$ } \leq D_{0}\\

1, & \text{if $D(u,v)$} \geq D_{0}

\end{cases}

\]

20131208162032343

def highPassFilter(image,d):

f = np.fft.fft2(image)

fshift = np.fft.fftshift(f)

def make_transform_matrix(d):

transfor_matrix = np.zeros(image.shape)

center_point = tuple(map(lambda x:(x-1)/2,s1.shape))

for i in range(transfor_matrix.shape[0]):

for j in range(transfor_matrix.shape[1]):

def cal_distance(pa,pb):

from math import sqrt

dis = sqrt((pa[0]-pb[0])**2+(pa[1]-pb[1])**2)

return dis

dis = cal_distance(center_point,(i,j))

if dis <= d:

transfor_matrix[i,j]=0

else:

transfor_matrix[i,j]=1

return transfor_matrix

d_matrix = make_transform_matrix(d)

new_img = np.abs(np.fft.ifft2(np.fft.ifftshift(fshift*d_matrix)))

return new_img

img_d1 = highPassFilter(img,10)

img_d2 = highPassFilter(img,30)

img_d3 = highPassFilter(img,50)

plt.subplot(131)

plt.axis("off")

plt.imshow(img_d1,cmap="gray")

plt.title('D_1 10')

plt.subplot(132)

plt.axis("off")

plt.title('D_2 30')

plt.imshow(img_d2,cmap="gray")

plt.subplot(133)

plt.axis("off")

plt.title("D_3 50")

plt.imshow(img_d3,cmap="gray")

plt.show()

0069Secoly1fm0gnw22igj30ai041t9y.jpg

显然当D0=10时,苹果的边缘最清楚

不同滤波的比较

import imagefilter

thread_img = imagefilter.RobertsAlogrithm(img)

laplace_img = imagefilter.LaplaceAlogrithm(img,"fourfields")

mean_img = cv2.blur(img,(3,3))

plt.subplot(131)

plt.imshow(thread_img,cmap="gray")

plt.title("ThreadImage")

plt.axis("off")

plt.subplot(132)

plt.imshow(laplace_img,cmap="gray")

plt.axis("off")

plt.title("LaplaceImage")

plt.subplot(133)

plt.imshow(mean_img,cmap="gray")

plt.title("meanImage")

plt.axis("off")

plt.show()

0069Secoly1fm0gnw2h96j30ai041wfo.jpg

空间域上的平均滤波和低通滤波一样,只要起去掉无关信息,平滑图像的作用。

Roberts,Laplace等滤波则起的提取边缘的作用。

频率域高通滤波器

高斯高通滤波器

频率域高斯高通滤波器的公式如下

\[H(u,v) = 1-e^{\dfrac{-D^2(u,v)}{2D_0^2}}

\]

20131208165324703

def GaussianHighFilter(image,d):

f = np.fft.fft2(image)

fshift = np.fft.fftshift(f)

def make_transform_matrix(d):

transfor_matrix = np.zeros(image.shape)

center_point = tuple(map(lambda x:(x-1)/2,s1.shape))

for i in range(transfor_matrix.shape[0]):

for j in range(transfor_matrix.shape[1]):

def cal_distance(pa,pb):

from math import sqrt

dis = sqrt((pa[0]-pb[0])**2+(pa[1]-pb[1])**2)

return dis

dis = cal_distance(center_point,(i,j))

transfor_matrix[i,j] = 1-np.exp(-(dis**2)/(2*(d**2)))

return transfor_matrix

d_matrix = make_transform_matrix(d)

new_img = np.abs(np.fft.ifft2(np.fft.ifftshift(fshift*d_matrix)))

return new_img

使用高斯滤波器d分别为10,30,50实现的效果

img_d1 = GaussianHighFilter(img,10)

img_d2 = GaussianHighFilter(img,30)

img_d3 = GaussianHighFilter(img,50)

plt.subplot(131)

plt.axis("off")

plt.imshow(img_d1,cmap="gray")

plt.title('D_1 10')

plt.subplot(132)

plt.axis("off")

plt.title('D_2 30')

plt.imshow(img_d2,cmap="gray")

plt.subplot(133)

plt.axis("off")

plt.title("D_3 50")

plt.imshow(img_d3,cmap="gray")

plt.show()

0069Secoly1fm0gnw4pg1j30ai041aba.jpg

高斯低通滤波器

频率域高斯低通滤波器的公式如下

\[H(u,v) = e^{\dfrac{-D^2(u,v)}{2D_0^2}}

\]

def GaussianLowFilter(image,d):

f = np.fft.fft2(image)

fshift = np.fft.fftshift(f)

def make_transform_matrix(d):

transfor_matrix = np.zeros(image.shape)

center_point = tuple(map(lambda x:(x-1)/2,s1.shape))

for i in range(transfor_matrix.shape[0]):

for j in range(transfor_matrix.shape[1]):

def cal_distance(pa,pb):

from math import sqrt

dis = sqrt((pa[0]-pb[0])**2+(pa[1]-pb[1])**2)

return dis

dis = cal_distance(center_point,(i,j))

transfor_matrix[i,j] = np.exp(-(dis**2)/(2*(d**2)))

return transfor_matrix

d_matrix = make_transform_matrix(d)

new_img = np.abs(np.fft.ifft2(np.fft.ifftshift(fshift*d_matrix)))

return new_img

img_d1 = GaussianLowFilter(img,10)

img_d2 = GaussianLowFilter(img,30)

img_d3 = GaussianLowFilter(img,50)

plt.subplot(131)

plt.axis("off")

plt.imshow(img_d1,cmap="gray")

plt.title('D_1 10')

plt.subplot(132)

plt.axis("off")

plt.title('D_2 30')

plt.imshow(img_d2,cmap="gray")

plt.subplot(133)

plt.axis("off")

plt.title("D_3 50")

plt.imshow(img_d3,cmap="gray")

plt.show()

0069Secoly1fm0gnw68kgj30ai041jru.jpg

空间域的高斯滤波

通常空间域使用高斯滤波来平滑图像,在上一篇已经写过,直接使用上篇文章的代码。

def GaussianOperator(roi):

GaussianKernel = np.array([[1,2,1],[2,4,2],[1,2,1]])

result = np.sum(roi*GaussianKernel/16)

return result

def GaussianSmooth(image):

new_image = np.zeros(image.shape)

image = cv2.copyMakeBorder(image,1,1,1,1,cv2.BORDER_DEFAULT)

for i in range(1,image.shape[0]-1):

for j in range(1,image.shape[1]-1):

new_image[i-1,j-1] =GaussianOperator(image[i-1:i+2,j-1:j+2])

return new_image.astype(np.uint8)

new_apple = GaussianSmooth(img)

plt.subplot(121)

plt.axis("off")

plt.title("origin image")

plt.imshow(img,cmap="gray")

plt.subplot(122)

plt.axis("off")

plt.title("Gaussian image")

plt.imshow(img,cmap="gray")

plt.subplot(122)

plt.axis("off")

plt.show()

0069Secoly1fm0gnw6gq8j30ai05jjs0.jpg

巴特沃斯滤波器

无论是低通滤波器,高通滤波器都是粗暴的一刀切,正如之前那么多空间域的滤波器一样,我们希望它通过的频率和与中心线性相关。

\[h(u,v) = \frac{1} {{1+(D_0 / D(u,v))}^{2n}}

\]

20131208165014296

def butterworthPassFilter(image,d,n):

f = np.fft.fft2(image)

fshift = np.fft.fftshift(f)

def make_transform_matrix(d):

transfor_matrix = np.zeros(image.shape)

center_point = tuple(map(lambda x:(x-1)/2,s1.shape))

for i in range(transfor_matrix.shape[0]):

for j in range(transfor_matrix.shape[1]):

def cal_distance(pa,pb):

from math import sqrt

dis = sqrt((pa[0]-pb[0])**2+(pa[1]-pb[1])**2)

return dis

dis = cal_distance(center_point,(i,j))

transfor_matrix[i,j] = 1/((1+(d/dis))**n)

return transfor_matrix

d_matrix = make_transform_matrix(d)

new_img = np.abs(np.fft.ifft2(np.fft.ifftshift(fshift*d_matrix)))

return new_img

plt.subplot(231)

butter_100_1 = butterworthPassFilter(img,100,1)

plt.imshow(butter_100_1,cmap="gray")

plt.title("d=100,n=1")

plt.axis("off")

plt.subplot(232)

butter_100_2 = butterworthPassFilter(img,100,2)

plt.imshow(butter_100_2,cmap="gray")

plt.title("d=100,n=2")

plt.axis("off")

plt.subplot(233)

butter_100_3 = butterworthPassFilter(img,100,3)

plt.imshow(butter_100_3,cmap="gray")

plt.title("d=100,n=3")

plt.axis("off")

plt.subplot(234)

butter_100_1 = butterworthPassFilter(img,30,1)

plt.imshow(butter_100_1,cmap="gray")

plt.title("d=30,n=1")

plt.axis("off")

plt.subplot(235)

butter_100_2 = butterworthPassFilter(img,30,2)

plt.imshow(butter_100_2,cmap="gray")

plt.title("d=30,n=2")

plt.axis("off")

plt.subplot(236)

butter_100_3 = butterworthPassFilter(img,30,3)

plt.imshow(butter_100_3,cmap="gray")

plt.title("d=30,n=3")

plt.axis("off")

plt.show()

0069Secoly1fm0gqhaufpj30ai07c771.jpg

可以明显的观察出过大的n造成的振铃现象

butter_5_1 = butterworthPassFilter(img,5,1)

plt.imshow(butter_5_1,cmap="gray")

plt.title("d=5,n=3")

plt.axis("off")

plt.show()

0069Secoly1fm0gqm6ronj307a07cmyn.jpg

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/Monodyee/article/detail/78358
推荐阅读
相关标签
  

闽ICP备14008679号