网站建设公司业务,望京做网站,高新公司网站建设电话,在深圳注册公司有什么好处目录频率域滤波基础频率域的其他特性频率域滤波基础知识频率域滤波步骤小结空间域和频率域滤波之间的对应关系频率域滤波基础 频率域的其他特性 频率域中的滤波过程如下#xff1a; 首先修改傅里叶变换以在到特定目的然后计算IDFT#xff0c;返回到空间域 # 频率域中的其… 目录频率域滤波基础频率域的其他特性频率域滤波基础知识频率域滤波步骤小结空间域和频率域滤波之间的对应关系频率域滤波基础 频率域的其他特性 频率域中的滤波过程如下 首先修改傅里叶变换以在到特定目的然后计算IDFT返回到空间域 # 频率域中的其他特性 img cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0429(a)(blown_ic).tif, -1)# FFT img_fft np.fft.fft2(img.astype(np.float32))# 非中心化的频谱 spectrum spectrum_fft(img_fft) spectrum np.uint8(normalize(spectrum) * 255)# 中心化 fshift np.fft.fftshift(img_fft) # 将变换的频率图像四角移动到中心# 中心化后的频谱 spectrum_fshift spectrum_fft(fshift) spectrum_fshift_n np.uint8(normalize(spectrum_fshift) * 255)# 对频谱做对数变换 spectrum_log np.log(1 spectrum_fshift)plt.figure(figsize(15, 8)) plt.subplot(1, 2, 1), plt.imshow(img, cmapgray), plt.title(Original), plt.xticks([]), plt.yticks([]) plt.subplot(1, 2, 2), plt.imshow(spectrum_log, cmapgray), plt.title(FFT Shift log), plt.xticks([]), plt.yticks([]) plt.tight_layout() plt.show()频率域滤波基础知识 已知大小为P×QP \times QP×Q像素的一幅数字图像f(x,y)f(x, y)f(x,y)则我们感兴趣取的基本滤波公式是 g(x,y)Real{J−1[H(u,v)F(u,v)]}(4.104)g(x, y) Real\{ \mathfrak{J}^{-1}[H(u, v) F(u, v)] \} \tag{4.104}g(x,y)Real{J−1[H(u,v)F(u,v)]}(4.104) J−1\mathfrak{J}^{-1}J−1是IDFTH(u,v)H(u, v)H(u,v)是滤波器传递函数经常称为滤波器或滤波器函数F(u,v)F(u, v)F(u,v)是输入图像f(x,y)f(x, y)f(x,y)的DFTg(x,y)g(x, y)g(x,y)是滤波后的图像。 中以化用(−1)xy(-1)^{x y}(−1)xy乘以输入图像来完成的 低通滤波器 会衰减高频而通过低频的函数会模糊图像 高通滤波器 将使图像的细节更清晰但会降低图像的对比度 def centralized_2d_loop(arr):centralized 2d array $f(x, y) (-1)^{x y}$ #中心化height, width arr.shape[:2]dst np.zeros([height, width])for h in range(height):for w in range(width):dst[h, w] arr[h, w] * ((-1) ** (h w))return dstdef centralized_2d_index(arr):centralized 2d array $f(x, y) (-1)^{x y}$, about 35 times faster than loopdef get_1_minus1(img):get 1 when image index is even, -1 while index is odd, same shape as input image, need this array to multiply with input imageto get centralized image for DFTParameter: img: input, here we only need img shape to create the arrayreturn such a [[1, -1, 1], [-1, 1, -1]] array, example is 3x3M, N img.shapem, n img.shape# 这里判断一个m, n是否是偶数如果是偶数的话结果不正确需要变换为奇数if m % 2 0 :m 1if n % 2 0 :n 1temp np.arange(m * n).reshape(m, n)dst np.piecewise(temp, [temp % 2 0, temp % 2 1], [1, -1])dst dst[:M, :N]return dst#中心化mask get_1_minus1(arr)dst arr * maskreturn dstdef centralized_2d(arr):centralized 2d array $f(x, y) (-1)^{x y}$, about 4.5 times faster than index, 160 times faster than loop,def get_1_minus1(img):get 1 when image index is even, -1 while index is odd, same shape as input image, need this array to multiply with input imageto get centralized image for DFTParameter: img: input, here we only need img shape to create the arrayreturn such a [[1, -1, 1], [-1, 1, -1]] array, example is 3x3dst np.ones(img.shape)dst[1::2, ::2] -1dst[::2, 1::2] -1return dst#中心化mask get_1_minus1(arr)dst arr * maskreturn dsta np.arange(25).reshape(5, 5) aarray([[ 0, 1, 2, 3, 4],[ 5, 6, 7, 8, 9],[10, 11, 12, 13, 14],[15, 16, 17, 18, 19],[20, 21, 22, 23, 24]])centralized_2d(a)array([[ 0., -1., 2., -3., 4.],[ -5., 6., -7., 8., -9.],[ 10., -11., 12., -13., 14.],[-15., 16., -17., 18., -19.],[ 20., -21., 22., -23., 24.]])def idea_high_pass_filter(source, radius5):M, N source.shape[1], source.shape[0]u np.arange(M)v np.arange(N)u, v np.meshgrid(u, v)D np.sqrt((u - M//2)**2 (v - N//2)**2)D0 radiuskernel D.copy()kernel[D D0] 1kernel[D D0] 0kernel_normal (kernel - kernel.min()) / (kernel.max() - kernel.min())return kernel_normal# 频率域滤波基础知识 img cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0429(a)(blown_ic).tif, -1)# FFT img_fft np.fft.fft2(img.astype(np.float32))# 中心化 fshift np.fft.fftshift(img_fft) # 将变换的频率图像四角移动到中心# 滤波器 h idea_high_pass_filter(img, radius1)# 滤波 res fshift * h# 先反中以化再反变换 f1shift np.fft.ifftshift(res) ifft np.fft.ifft2(f1shift)# 取实数部分 recon np.real(ifft) # 小于0都被置为0 # recon np.where(recon 0, recon, 0) recon np.clip(recon, 0, recon.max()) recon np.uint8(normalize(recon) * 255)plt.figure(figsize(15, 8)) plt.subplot(1, 2, 1), plt.imshow(img, cmapgray), plt.title(Original), plt.xticks([]), plt.yticks([]) plt.subplot(1, 2, 2), plt.imshow(recon, cmapgray), plt.title(FFT Shift log), plt.xticks([]), plt.yticks([]) plt.tight_layout() plt.show()def gauss_high_pass_filter(source, center, radius5):create gaussian high pass filter param: source: input, source imageparam: center: input, the center of the filter, where is the lowest value, (0, 0) is top left corner, source.shape[:2] is center of the source imageparam: radius: input, the radius of the lowest value, greater value, bigger blocker out range, if the radius is 0, then allvalue is 1return a [0, 1] value filterM, N source.shape[1], source.shape[0]u np.arange(M)v np.arange(N)u, v np.meshgrid(u, v)D np.sqrt((u - center[1]//2)**2 (v - center[0]//2)**2)D0 radiusif D0 0:kernel np.ones(source.shape[:2], dtypenp.float)else:kernel 1 - np.exp(- (D**2)/(D0**2)) kernel (kernel - kernel.min()) / (kernel.max() - kernel.min())return kerneldef plot_3d(ax, x, y, z):ax.plot_surface(x, y, z, antialiasedTrue, shadeTrue)ax.view_init(30, 60), ax.grid(bFalse), ax.set_xticks([]), ax.set_yticks([]), ax.set_zticks([])def imshow_img(img, ax, cmapgray):ax.imshow(img, cmapcmap), ax.set_xticks([]), ax.set_yticks([])def frequency_filter(fshift, h):Frequency domain filtering using FFTparam: fshift: input, FFT shift to centerparam: h: input, filter, value range [0, 1]return a uint8 [0, 255] reconstruct image# 滤波res fshift * h# 先反中以化再反变换f1shift np.fft.ifftshift(res)ifft np.fft.ifft2(f1shift)# 取实数部分recon np.real(ifft)# 小于0都被置为0# recon np.where(recon 0, recon, 0)recon np.clip(recon, 0, recon.max())recon np.uint8(normalize(recon) * 255)return recon# 高斯核滤波器与对应的结果 from mpl_toolkits.mplot3d import Axes3D import numpy as np from matplotlib import pyplot as plt from matplotlib import cm# 显示多个3D高斯核函数 k 3 sigma 0.8 X np.linspace(-k, k, 200) Y np.linspace(-k, k, 200) x, y np.meshgrid(X, Y) z np.exp(- ((x)**2 (y)**2)/ (2 * sigma**2))fig plt.figure(figsize(15, 10)) ax_1 fig.add_subplot(2, 3, 1, projection3d) plot_3d(ax_1, x, y, z)ax_2 fig.add_subplot(2, 3, 2, projection3d) plot_3d(ax_2, x, y, 1 - z)# 这里只显示不明显 a 1 x a x y a y ax_3 fig.add_subplot(2, 3, 3, projection3d) plot_3d(ax_3, x, y, 1 - z)img_ori cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0431(d)(blown_ic_crop).tif, -1)# FFT img_fft np.fft.fft2(img_ori.astype(np.float32))# 中心化 fshift np.fft.fftshift(img_fft) # 将变换的频率图像四角移动到中心# 滤波器 radius 20 center img_ori.shape[:2] h_high_1 gauss_high_pass_filter(img_ori, center, radiusradius) h_low 1 - h_high_1 center (1026, 872) h_high_2 gauss_high_pass_filter(img_ori, center, radiusradius)img_low frequency_filter(fshift, h_low) img_high_1 frequency_filter(fshift, h_high_1) img_high_2 frequency_filter(fshift, h_high_2)ax_4 fig.add_subplot(2, 3, 4) imshow_img(img_low, ax_4)ax_5 fig.add_subplot(2, 3, 5) imshow_img(img_high_1, ax_5)ax_6 fig.add_subplot(2, 3, 6) imshow_img(img_high_2, ax_6)plt.tight_layout() plt.show()下面展示的是函数未填充时会产生交叠错误而函数填充0后产生我们期望的结果但是此时的滤波器应该是原来的2倍大小如果是还是原来的大小的话那么图像会比之前更模糊的。 这里采用的方法是对图像进行零填充再直接在频率域中创建期望的滤波器传递函数该函数的到底会给你个填充零后的图像的大小相同。这个明显降低交叠错误只是会出现振铃效应。 # 未填充与填充0的结果 img_ori cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0432(a)(square_original).tif, -1)#未填充 # FFT img_fft np.fft.fft2(img_ori.astype(np.float32))# 中心化 fshift np.fft.fftshift(img_fft) # 将变换的频率图像四角移动到中心# 滤波器 radius 50 center img_ori.shape[:2] h_high_1 gauss_high_pass_filter(img_ori, center, radiusradius) h_low 1 - h_high_1img_low frequency_filter(fshift, h_low)fig plt.figure(figsize(15, 5)) ax_1 fig.add_subplot(1, 3, 1) imshow_img(img_ori, ax_1)ax_2 fig.add_subplot(1, 3, 2) imshow_img(img_low, ax_2)#填充0 fp np.zeros([img_ori.shape[0]*2, img_ori.shape[1]*2]) fp[:img_ori.shape[0], :img_ori.shape[1]] img_oriimg_fft np.fft.fft2(fp.astype(np.float32)) fshift np.fft.fftshift(img_fft) hp 1 - gauss_high_pass_filter(fp, fp.shape, radius100)img_low frequency_filter(fshift, hp) img_dst img_low[:img_ori.shape[0], :img_ori.shape[1]]ax_3 fig.add_subplot(1, 3, 3) imshow_img(img_dst, ax_3)plt.tight_layout() plt.show()另一个合理的做法是构建一个与未填充图像大小相同的函数计算该函数的IDFT得到相应的空间表示在空间域中对这一表示填充然后计算其DFT返回到频率域。这也会出现振铃效应 def centralized(arr):centralized 1d array $f(x) (-1)^{n}$u arr.shape[0]dst np.zeros(arr.shape)# 中心化for n in range(u):dst[n] arr[n] * ((-1) ** n)return dst# 频率域反变换到空间域再填充再DFT会出现振铃效应 h np.zeros([256]) h[:9] 1 h_c np.fft.fftshift(h)fig plt.figure(figsize(15, 10)) grid plt.GridSpec(2, 3, wspace0.2, hspace0.1) ax_1 fig.add_subplot(grid[0, 0]) ax_1.plot(h_c), ax_1.set_xticks([0, 128, 255]), ax_1.set_yticks(np.arange(-0.2, 1.4, 0.2)), ax_1.set_xlim([0, 255])# 用公式法的中以化才能得到理想的结果 h_pad np.pad(h, (0, 256), modeconstant, constant_values0) ifshift centralized_2d(h_pad.reshape(1, -1)).flatten() # 也可以用centralized_2d这个函数 ifft np.fft.ifft(ifshift) ifft ifft.real * 2 ifft[:128] 0 ifft[384:] 0 ax_2 fig.add_subplot(grid[0, 1:3]) ax_2.plot(ifft), ax_2.set_xticks([0, 128, 256, 384, 511]), ax_2.set_yticks(np.arange(-0.01, 0.05, 0.01)), ax_2.set_xlim([0, 511])f ifft[128:384] ax_3 fig.add_subplot(grid[1, 0]) ax_3.plot(f), ax_3.set_xticks([0, 128, 255]), ax_3.set_yticks(np.arange(-0.01, 0.05, 0.01)), ax_3.set_xlim([0, 255])# 已经中心化截取部分构成原函数但出现的效果不是太好但可以看出振铃效应 d 120 f_ np.zeros([h.shape[0]]) f_[:256-d] f[d:] f_pad np.pad(f_, (0, 256), modeconstant, constant_values0) f_centralized centralized(f_pad) fft np.fft.fft(f_centralized) ax_4 fig.add_subplot(grid[1, 1:3]) ax_4.plot(fft.real), ax_4.set_xticks([0, 128, 256, 384, 511]), ax_4.set_yticks(np.arange(-0.2, 1.4, 0.2)), ax_4.set_xlim([0, 511])plt.show()零相移滤波器 相角计算为复数的实部与虚部之比的反正切。因为H(u,v)H(u, v)H(u,v)同是乘以R,IR,IR,I所以当这个比率形成后它会被抵消。对实部和虚部的影响相同且对相角无影响的滤波器。 # 相角对反变换的影响 img_ori cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0431(d)(blown_ic_crop).tif, -1)#未填充 # FFT img_fft np.fft.fft2(img_ori.astype(np.float32))# 中心化 fshift np.fft.fftshift(img_fft) # 频谱与相角 spectrum spectrum_fft(fshift) phase phase_fft(fshift)# 相角乘以-1 phase_1 phase * -1 new_fshift spectrum * np.exp(phase_1 * -1j) new_fshift np.fft.ifftshift(new_fshift) recon_1 np.fft.fft2(new_fshift)# 相角乘以0.5需要乘以一个倍数以便可视化 phase_2 phase * 0.25 new_fshift spectrum * np.exp(phase_2 * -1j) new_fshift np.fft.ifftshift(new_fshift) recon_2 np.fft.fft2(new_fshift) recon_2 np.uint8(normalize(recon_2.real) * 4000)fig plt.figure(figsize(15, 5)) ax_1 fig.add_subplot(1, 3, 1) imshow_img(img_ori, ax_1)ax_2 fig.add_subplot(1, 3, 2) imshow_img(recon_1.real, ax_2)ax_3 fig.add_subplot(1, 3, 3) ax_3.imshow(recon_2, cmapgray)plt.tight_layout() plt.show()频率域滤波步骤小结 已知一幅大小为M×NM\times NM×N的输入图像f(x,y)f(x, y)f(x,y)分别得到填充零后的尺寸PPP和QQQ即P2MP2MP2M和Q2NQ2NQ2N。使用零填充、镜像填充或复制填充形成大小为P×QP\times QP×Q的填充后的图像fP(x,y)f_P(x, y)fP​(x,y)。将fP(x,y)f_P(x, y)fP​(x,y)乘以(−1)xy(-1)^{xy}(−1)xy使用傅里叶变换位于P×QP\times QP×Q大小的频率矩形的中心。计算步骤3得到的图像的DFT即F(u,v)F(u, v)F(u,v)。构建一个实对称滤波器传递函数H(u,v)H(u, v)H(u,v)其大小为P×QP\times QP×Q中心中(P/2,Q/2)(P/2, Q/2)(P/2,Q/2)处。采用对应像素相乘得到G(u,v)H(u,v)F(u,v)G(u, v) H(u, v)F(u, v)G(u,v)H(u,v)F(u,v)。计算G(u,v)G(u, v)G(u,v)的IDFT得到滤波后的图像大小为P×QP\times QP×QgP(x,y){Real[J−1[G(u,v)]]}(−1)xyg_P(x, y) \Big\{\text{Real} \big[\mathfrak{J}^{-1}[G(u, v)] \big]\Big\}(-1)^{xy}gP​(x,y){Real[J−1[G(u,v)]]}(−1)xy。最后从gP(x,y)g_P(x, y)gP​(x,y)的左上象限提取一个大小为M×NM\times NM×N的区域得到与输入图像大小相同的滤波后的结果g(x,y)g(x, y)g(x,y)。 def pad_image(img, modeconstant):pad image into PxQ shape, orginal is in the top left corner.param: img: input imageparam: mode: input, str, is numpy pad parameter, default is constant. for more detail please refere to Numpy padreturn PxQ shape image padded with zeros or other valuesdst np.pad(img, ((0, img.shape[0]), (0, img.shape[1])), modemode)return dst # 频率域滤波过程 img_ori cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0431(d)(blown_ic_crop).tif, -1) M, N img_ori.shape[:2] fig plt.figure(figsize(15, 15))# 这里对原图像进行pad以得到更好的显示 img_ori_show np.ones([img_ori.shape[0]*2, img_ori.shape[1]*2]) * 255 img_ori_show[:img_ori.shape[0], img_ori.shape[1]:] img_ori ax_1 fig.add_subplot(3, 3, 1) imshow_img(img_ori_show, ax_1)fp pad_image(img_ori, modeconstant) ax_2 fig.add_subplot(3, 3, 2) imshow_img(fp, ax_2)fp_cen centralized_2d(fp) fp_cen_show np.clip(fp_cen, 0, 255) # 负数的都截断为0 ax_3 fig.add_subplot(3, 3, 3) ax_3.imshow(fp_cen_show, cmapgray), ax_3.set_xticks([]), ax_3.set_yticks([])fft np.fft.fft2(fp_cen) spectrum spectrum_fft(fft) spectrum np.log(1 spectrum) ax_4 fig.add_subplot(3, 3, 4) imshow_img(spectrum, ax_4)H 1 - gauss_high_pass_filter(fp, centerfp.shape, radius70) ax_5 fig.add_subplot(3, 3, 5) imshow_img(H, ax_5)HF fft * H HF_spectrum np.log(1 spectrum_fft(HF)) ax_6 fig.add_subplot(3, 3, 6) imshow_img(HF_spectrum, ax_6)ifft np.fft.ifft2(HF) gp centralized_2d(ifft.real) ax_7 fig.add_subplot(3, 3, 7) imshow_img(gp, ax_7)# 最终结果有黑边 g gp[:M, :N] ax_8 fig.add_subplot(3, 3, 8) imshow_img(g, ax_8)plt.tight_layout() plt.show()def frequency_filter(fft, h):Frequency domain filtering using FFTparam: fshift: input, FFT shift to centerparam: h: input, filter, value range [0, 1]return a uint8 [0, 255] reconstruct image# 滤波res fft * h# 先反中以化再反变换ifft np.fft.ifft2(res)# 取实数部分recon centralized_2d(ifft.real)# 小于0都被置为0recon np.clip(recon, 0, recon.max())recon np.uint8(normalize(recon) * 255)return recon# 频率域滤波 img_ori cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0431(d)(blown_ic_crop).tif, -1) M, N img_ori.shape[:2]# 填充 fp pad_image(img_ori, modereflect) # 中心化 fp_cen centralized_2d(fp) # 正变换 fft np.fft.fft2(fp_cen) # 滤波器 H 1 - gauss_high_pass_filter(fp, centerfp.shape, radius100) # 滤波 HF fft * H # 反变换 ifft np.fft.ifft2(HF) # 去中心化 gp centralized_2d(ifft.real) # 还回来与原图像的大小 g gp[:M, :N] fig plt.figure(figsize(15, 15)) ax_8 fig.add_subplot(3, 3, 8) imshow_img(g, ax_8) plt.tight_layout() plt.show()空间域和频率域滤波之间的对应关系 空间域中的滤波器核与频率域傅里叶变换对 h(x,y)⇔H(u,v)(4.106)h(x, y) \Leftrightarrow H(u, v) \tag{4.106}h(x,y)⇔H(u,v)(4.106) h(x,y)h(x, y)h(x,y)是空间核。这个核 可以由一个频率域滤波器对一个冲激的响应得到因此有时我们称h(x,y)h(x, y)h(x,y)是H(u,v)H(u, v)H(u,v)的冲激响应。 一维频率域高斯传递函数 H(u)Ae−u2/2σ2(4.107)H(u) A e^{-u^2/2\sigma^2} \tag{4.107}H(u)Ae−u2/2σ2(4.107) 空间域的核 由H(u)H(u)H(u)的傅里叶反变换得到 h(x)2πσAe−2π2σ2x2(4.108)h(x) \sqrt{2\pi}\sigma A e^{-2\pi^2\sigma^2 x^2} \tag{4.108}h(x)2π​σAe−2π2σ2x2(4.108) 这是两个很重要的公式 它们是一个傅里叶变换对两者都是高斯函数和实函数两个函数的表现相反 高斯差 用一个常量减去低通函数可以构建一个高通滤波器。因为使用所谓的高斯差涉及两个低通函数就可更好地控制滤波函数的形状 频率域 A≥B,σ1σ2A \ge B, \sigma_1 \sigma_2A≥B,σ1​σ2​ H(u)Ae−u2/2σ12−Be−u2/2σ22(4.109)H(u) A e^{-u^2/2\sigma_1^2} - B e^{-u^2/2\sigma_2^2}\tag{4.109}H(u)Ae−u2/2σ12​−Be−u2/2σ22​(4.109) 空间域 h(x)2πσ1Ae−2π2σ12x2−2πσ2Be−2π2σ22x2(4.110)h(x) \sqrt{2\pi}\sigma_1 A e^{-2\pi^2\sigma_1^2 x^2} - \sqrt{2\pi}\sigma_2 B e^{-2\pi^2\sigma_2^2 x^2} \tag{4.110}h(x)2π​σ1​Ae−2π2σ12​x2−2π​σ2​Be−2π2σ22​x2(4.110) # 空间域与频率域的滤波器对应关系 u np.linspace(-3, 3, 200)A 1 sigma 0.2 Hu A * np.exp(-np.power(u, 2) / (2 * np.power(sigma, 2)))x u hx np.sqrt(2 * np.pi) * sigma * A * np.exp(-2 * np.pi**2 * sigma**2 * x**2)fig plt.figure(figsize(10, 10)) ax_1 setup_axes(fig, 221) ax_1.plot(u, Hu), ax_1.set_title(H(u), loccenter, y1.05), ax_1.set_ylim([0, 1]), ax_1.set_yticks([])B 1 sigma_2 4 Hu_2 B * np.exp(-np.power(u, 2) / (2 * np.power(sigma_2, 2))) Hu_diff Hu_2 - Huax_2 setup_axes(fig, 222) ax_2.plot(x, Hu_diff), ax_2.set_title(H(u), loccenter, y1.05), ax_2.set_ylim([0, 1.1]), ax_2.set_yticks([])ax_3 setup_axes(fig, 223) ax_3.plot(x, hx), ax_3.set_title(h(x), loccenter, y1.05), ax_3.set_ylim([0, 0.8]), ax_3.set_yticks([])hx_2 np.sqrt(2 * np.pi) * sigma_2 * B * np.exp(-2 * np.pi**2 * sigma_2**2 * x**2) hx_diff hx_2 - hxax_4 setup_axes(fig, 224) ax_4.plot(x, hx_diff), ax_4.set_title(h(x), loccenter, y1.05), ax_4.set_ylim([-1, 10]), ax_4.set_yticks([])plt.tight_layout() plt.show()# 空间核得到频率域 img_ori cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0438(a)(bld_600by600).tif, -1)fp pad_zero(img_ori) fp_cen centralized_2d(fp) fft np.fft.fft2(fp_cen) spectrum spectrum_fft(fft) spectrum np.log(1 spectrum)plt.figure(figsize(18, 15)) plt.subplot(1, 2, 1), plt.imshow(img, cmapgray), plt.title(No Log) plt.subplot(1, 2, 2), plt.imshow(spectrum, cmapgray), plt.title(FFT) plt.tight_layout() plt.show()注频率域的Sobel滤波器还没有解决 # 频域 center (fp.shape[0], fp.shape[1] 300) GHPF gauss_high_pass_filter(fp, center, radius200) GHPF (1 - np.where(GHPF 1e-6, GHPF, 1)) * -2center (fp.shape[0], fp.shape[1] - 300) GHPF_1 gauss_high_pass_filter(fp, center, radius200) GHPF_1 (1 - np.where(GHPF_1 1e-6, GHPF_1, 1)) * 2sobel_f GHPF GHPF_1plt.figure(figsize(15, 5))plt.subplot(1, 3, 1), plt.imshow(GHPF, gray) plt.subplot(1, 3, 2), plt.imshow(GHPF_1, gray) plt.subplot(1, 3, 3), plt.imshow(sobel_f, gray)plt.tight_layout() plt.show()# 频域 fig plt.figure(figsize(15, 5)) HF fft * sobel_f HF_spectrum np.log(1 spectrum_fft(HF)) ax_6 fig.add_subplot(1, 3, 1) imshow_img(HF_spectrum, ax_6)ifft np.fft.ifft2(HF) gp centralized_2d(ifft.real) ax_7 fig.add_subplot(1, 3, 2) imshow_img(gp, ax_7)# 最终结果有黑边 g gp[:M, :N] ax_8 fig.add_subplot(1, 3, 3) imshow_img(g, ax_8)# plt.figure(figsize(15, 5))# plt.subplot(1, 3, 1), plt.imshow(GHPF, gray) # plt.subplot(1, 3, 2), plt.imshow(GHPF_1, gray) # plt.subplot(1, 3, 3), plt.imshow(sobel_f, gray)plt.tight_layout() plt.show()空间域的Sobel滤波器 def kernel_seperate(kernel):这里的分离核跟之前不太一样get kernel seperate vector if separableparam: kernel: input kerel of the filterreturn two separable vector c, rrank np.linalg.matrix_rank(kernel)if rank 1:#-----------------Numpy------------------# here for sobel kernel seperatec kernel[:, -1]r kernel[-1, :]else:print(Not separable)# c np.array(c).reshape(-1, 1) # r np.array(r).reshape(-1, 1)return c, rdef separate_kernel_conv2D(img, kernel):separable kernel convolution 2D, if 2D kernel not separable, then canot use this fuction. in the fuction. first need toseperate the 2D kernel into two 1D vectors.param: img: input image you want to implement 2d convolution with 2d kernelparam: kernel: 2D separable kernel, such as Gaussian kernel, Box Kernelreturn image after 2D convolutionm, n kernel.shape[:2]pad_h int((m - 1) / 2)pad_w int((n - 1) / 2)w_1, w_2 kernel_seperate(kernel) # w_1 np.array([-1, 0, 1]) # w_2 np.array([1, 2, 1])convovle np.vectorize(np.convolve, signature(x),(y)-(k))r_2 convovle(img, w_2)r_r np.rot90(r_2)r_1 convovle(r_r, w_1)r_1 r_1.Tr_1 r_1[:, ::-1]img_dst img.copy().astype(np.float)img_dst r_1[pad_h:-pad_h, pad_w:-pad_w]return img_dst# Sobel梯度增强边缘 img_ori cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0438(a)(bld_600by600).tif, 0)sobel_x np.zeros([3, 3], np.int) sobel_x[0, :] np.array([-1, -2, -1]) sobel_x[2, :] np.array([1, 2, 1])sobel_y np.zeros([3, 3], np.int) sobel_y[:, 0] np.array([-1, -2, -1]) sobel_y[:, 2] np.array([1, 2, 1])gx separate_kernel_conv2D(img_ori, kernelsobel_x) gy separate_kernel_conv2D(img_ori, kernelsobel_y)# gx cv2.filter2D(img_ori, ddepth-1, kernelsobel_x) # gy cv2.filter2D(img_ori, ddepth-1, kernelsobel_y)# 先对gx gy做二值化处理再应用下面的公式 img_sobel np.sqrt(gx**2 gy**2) # 二值化后平方根的效果与绝对值很接近 img_sobel abs(gx) abs(gy) img_sobel np.uint8(normalize(img_sobel) * 255)plt.figure(figsize(15, 12)) plt.subplot(1, 2, 1), plt.imshow(img_ori, gray, vmax255), plt.title(Original), plt.xticks([]), plt.yticks([]) plt.subplot(1, 2, 2), plt.imshow(img_sobel, gray, vmax255), plt.title(Sobel), plt.xticks([]), plt.yticks([]) plt.tight_layout() plt.show()