fkdtd算子可以使用如下滤波器在频率域实现：
H ( u , v ) = − 4 ∗ π 2 ( u 2 + v 2 ) H(u,v) = -4*pi^2 (u^{2} + v^{2}) H(u,v)=−4∗π2(u2+v2)
或者，关于频率矩阵的中心，使用如下滤波器：
H ( u , v ) = − 4 ∗ π 2 ( ( u − P 2 ) 2 + ( v − Q 2 ) 2 ) H(u,v) = -4*pi^2 ((u - frac{P}{2})^{2} +( v - frac{Q}{2})^{2}) H(u,v)=−4∗π2((u−2P)2+(v−2Q)2)
其中，D(u,v)是上式给出的距离函数。然后，fkdtd图像由下式得到：
▽ 2 f ( x , y ) = ∋ − 1 H ( u , v ) ∗ F ( u , v ) bigtriangledown^{2} f(x,y) = owns^{-1}{H(u,v)*F(u,v)} ▽2f(x,y)=∋−1H(u,v)∗F(u,v)
其中，F(u,v) 是 f(x,y) 的傅里叶变换。可以使用下式进行增强：
g ( x , y ) = f ( x , y ) + c ▽ 2 f ( x , y ) g(x,y) = f(x,y) + cbigtriangledown^{2} f(x,y) g(x,y)=f(x,y)+c▽2f(x,y)
原始图片：
代码1：
#!/usr/ygdhk/env python# -*- encoding: utf-8 -*-'''@文件 :频率域的fkdtd算子.py@说明 :@时间 :2021/05/14 20:31:55@作者 :zt@版本 :1.0'''import cv2import mathimport numpy as npfrom matplotlib import pyplot as pltfrom matplotlib import rcParamsrcParams['font.family'] = 'simhei'path = r"图片位置"#读取图片，以灰度级的形式img = cv2.imread(path,cv2.IMREAD_GRAYSCALE)plt.imshow(img, cmap = 'gray', interpolation = 'bicubic')plt.title("原始图像")plt.xticks([]), plt.yticks([])plt.show()# '''# 归一化到[0,1]# '''# img_max = max(max(row) for row in img)# img_gui = img/img_maxrows, cols = img.shapeimg_k = np.zeros((2*rows,2*cols)) #用来制作傅里叶变换"""标定"""img_min = np.min(img)fm = img - img_minimg_max = np.max(fm)fs = (fm/img_max)img_k[:rows,:cols] = fsprint("扩充后：")print(img_k.shape)print(img_k)plt.imshow(img_k, cmap = 'gray', interpolation = 'bicubic')plt.title("扩充图像")plt.xticks([]), plt.yticks([])plt.show()'''(-1)(x+y)乘以img_k移到其中心'''# for i in range(img_k.shape[0]):# for j in range(img_k.shape[1]):# if (i+j)%2!=0:# img_k[i][j] = (-1)*img_k[i][j]crows, ccols =int(rows), int(cols) #得到频率域的中心点坐标print(crows,ccols)'''构建h(u,v)，并求出来'''h = np.zeros((2*rows,2*cols))for i in range(rows): for j in range(cols): h[i][j] = 1+4*math.pow(math.pi,2)*(math.pow(i-crows,2)+math.pow(j-ccols,2))print("h(u,v):")print(h[crows,ccols])ff = np.fft.fft2(img_k) #图片扩充移其中心后的傅里叶变换fshift = np.fft.fftshift(ff) #将零频率分量移到频谱中心h_f = np.zeros_like(h)for i in range(len(h)): for j in range(len(h[0])): h_f[i,j] = h[i,j] * fshift[i,j]ishift = np.fft.ifftshift(h_f) #ishift是一个复数数组i_img = np.fft.ifft2(ishift) #逆傅里叶变换，得到仍然是一个复数数组i_img = np.abs(i_img) #取模，与取实部差不多plt.imshow(i_img, cmap = 'gray', interpolation = 'bicubic')plt.title("nibianhua")plt.xticks([]), plt.yticks([])plt.show()c = 1g = img + c*i_img[:rows,:cols] #总计算公式plt.imshow(g, cmap = 'gray', interpolation = 'bicubic') #绘图plt.title("频率域fkdtd处理后")plt.xticks([]), plt.yticks([])plt.show()效果：
代码2：
#!/usr/ygdhk/env python# -*- encoding: utf-8 -*-'''@文件 :频率域的fkdtd算子.py@说明 :@时间 :2021/05/14 21:31:50@作者 :zt@版本 :1.0'''import cv2import mathimport numpy as npfrom matplotlib import pyplot as pltfrom matplotlib import rcParamsrcParams['font.family'] = 'simhei'path = r"图片位置"#读取图片，以灰度级的形式img = cv2.imread(path,cv2.IMREAD_GRAYSCALE)"""值域[0,1]"""img_min = np.min(img)fm = img - img_minimg_max = np.max(fm)fs = (fm/img_max)rows, cols = img.shapeimg_k = np.zeros((2*rows,2*cols)) #用来制作傅里叶变换img_k[:rows,:cols] = fsff = np.fft.fft2(img_k) #傅里叶变换fshift = np.fft.fftshift(ff) #将零频率分量移到频谱中心m,n = fshift.shapeo_x,o_y = int(m/2),int(n/2)a_max = np.max(fshift)img_zero = np.zeros((m,n)) #滤波后的频谱for i in range(m): for j in range(n): h = (-4)*math.pow(math.pi,2)*(math.pow(i-o_x,2)+math.pow(j-o_y,2))/a_max #fkdtd滤波器，除以img_f的最大值a_max是为了将H(i,j)标定到[-1,1]的范围内 img_zero[i,j] = (1-h)*fshift[i,j]ishift = np.fft.ifftshift(img_zero) #ishift是一个复数数组i_img = np.fft.ifft2(ishift) #逆傅里叶变换，得到仍然是一个复数数组，不要以为直接得到图像的像素值.i_img = np.abs(i_img)plt.imshow(i_img[:rows,:cols], cmap = 'gray', interpolation = 'bicubic')plt.title("频率域fkdtd处理后")plt.xticks([]), plt.yticks([])plt.show()效果：