- 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()
使用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)\)点到中心的距离,\(D_0\)由自己设置
白点就是所允许通过的频率范围
3d图像如下
我们先把苹果转化成频率域看下效果
- #取绝对值:将复数变化成实数
- #取对数的目的为了将数据变化到0-255
- s1 = np.log(np.abs(fshift))
- plt.subplot(121),plt.imshow(s1,‘gray‘)
- plt.title(‘Frequency Domain‘)
- plt.show()
matplotlib对于不是uint8的图像会自动把图像的数值缩放到0-255上,更多可以查看对该问题的讨论
我们在频率域上试着取不同的\(d_0\)再将其反变换到空间域看下效果
- 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()
- 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()
讲上面过程整理得到频率域低通滤波器的代码如下
- 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")
高通滤波器同低通滤波器非常类似,只不过二者通过的波正好是相反的
\[
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}
\]
- 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()
显然当\(D_0=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()
空间域上的平均滤波和低通滤波一样,只要起去掉无关信息,平滑图像的作用。
Roberts,Laplace等滤波则起的提取边缘的作用。
频率域高斯高通滤波器的公式如下
\[
H(u,v) = 1-e^{\dfrac{-D^2(u,v)}{2D_0^2}}
\]
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()
频率域高斯低通滤波器的公式如下
\[
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()
通常空间域使用高斯滤波来平滑图像,在上一篇已经写过,直接使用上篇文章的代码。
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()
无论是低通滤波器,高通滤波器都是粗暴的一刀切,正如之前那么多空间域的滤波器一样,我们希望它通过的频率和与中心线性相关。
\[
h(u,v) = \frac{1} {{1+(D_0 / D(u,v))}^{2n}}
\]
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()
可以明显的观察出过大的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()
来源: http://www.bubuko.com/infodetail-2412923.html