傅里叶变换可以得到图像的频率信息,频率信息代表了图像像素变化的快慢,通过对图像频率信息的处理,可以实现图像滤波的效果。
频域滤波的步骤如下
(1)给定一幅M×N大小的图像$f(x,y)$,选择填充参数$P=2M$,$Q=2N$。
(2)对$f(x,y)$填充必要的0,形成大小为M×N的填充图像$f_p(x,y)$。
(3)用$(-1)^{x+y}$乘以$f_p(x,y)$移到其变换的中心。
(4)计算来自步骤3图像的DFT,得到$F(u,v)。
(5)生成一个实的对称的滤波函数$H(u,v)$,其大小为$P*Q$,中心在$(\frac{M}{2},\frac{N}{2})$。用阵列相乘得到$G(u,v)=F(u,v)*H(u,v)$。
(6)得到处理后的图像
(7)通过从$f_p(x,y)$左上象限提取c处理结果$g(x,y)$。
进行图像填充和频率中心平移的自定义函数如下
# 将图像填充到size大小,默认为填充到原图像的两倍 def my_fill(image, size=(-1,-1)): raw = 2 * image.shape[0] column = 2 * image.shape[1] if size[0] > 0: raw, column = size # 这里的数据结构可能更改 filled_image = np.zeros((raw, column), dtype=np.int64) for i in range(image.shape[0]): for j in range(image.shape[1]): filled_image[i][j] = image[i][j] return filled_image # 图像从offset截取size大小,默认从左上角截取一半 def my_ifill(image, size=(-1, -1), offset=(0, 0)): raw = int(image.shape[0] / 2) column = int(image.shape[1] / 2) if size[0] > 0: raw, column = size ifilled_image = np.zeros((raw, column), dtype=np.int64) for i in range(raw): for j in range(column): ifilled_image[i][j] = round(image[i+offset[0]][j+offset[1]]) return ifilled_image
# 对像素乘以(-1)^(x+y) def my_get_fp(image): image_fp = image.copy() raw, column = image.shape for i in range(raw): for j in range(column): image_fp[i][j] = image[i][j] * pow(-1, i + j) return image_fp
按照上述处理步骤自定义的频域滤波函数如下
# 自定义模板频域滤波函数 def myfunc_seldifine(image, model, output_offset=(0, 0)): """ 用于得到频域滤波的结果,滤波模板需要自己设计 :param image: 输入图像,不用填充 :param model: 频域滤波模板,需要将自己设计好的模板输入,尺寸大于图像尺寸。在函数中图像会根据模板大小进行填充 :param output_offset: 在进行反变换后的大尺寸图像截取时的偏移量 :return: 滤波后的输出图像,尺寸与原图像相同 """ # 填充图像 image_fill = my_fill(image, size=model.shape) # 将频率中心移到中点 image_fp = my_get_fp(image_fill) # 求解fft image_fft = np.fft.fft2(image_fp) fft_cal = model * image_fft # 进行傅里叶反变换 image_ifft = np.fft.ifft2(fft_cal) # 将图像移到中心点 image_ifft_real = image_ifft.real image_ifp = my_get_fp(image_ifft_real) result = my_ifill(image_ifp, size=image.shape, offset=output_offset) return result # 典型频域滤波函数 def myfunc_filter(image, d0=60, filetype = filter_type.butterworth_low_pass, n=1): """ 用于得到频域滤波的结果,可以选择典型的滤波类型。 :param image: 输入源图像 :param d0: 滤波范围 :param filetype: 滤波类型,有以下几种 filter_type.butterworth_low_pass:巴特沃斯低通 filter_type.butterworth_high_pass:巴特沃斯高通 filter_type.gaussion_low_pass:高斯低通 filter_type.gaussion_high_pass:高斯高通 filter_type.laplace_pass:Laplace滤波 :param n: 滤波器的阶数 :return: 滤波后的输出图像,和原图像尺寸相同。 """ # 填充图像 image_fill = my_fill(image) # 将频率中心移到中点 image_fp = my_get_fp(image_fill) # 求解fft image_fft = np.fft.fft2(image_fp) # 产生巴特沃斯滤波模板 model = my_get_butterworth_model(image_fft.shape, d0, n) if filetype == filter_type.butterworth_high_pass: model = my_get_butterworth_highpass_model(image_fft.shape, d0, n) elif filetype == filter_type.gaussion_low_pass: model = my_get_gaussion_model(image_fft.shape, d0) elif filetype == filter_type.gaussion_high_pass: model = my_get_gaussion_highpass_model(image_fft.shape, d0) elif filetype == filter_type.laplace_pass: model = my_get_laplace_model(image_fft.shape) # 进行频域滤波 fft_cal = model * image_fft # 进行傅里叶反变换 image_ifft = np.fft.ifft2(fft_cal) # 将图像移到中心点 image_ifft_real = image_ifft.real image_ifp = my_get_fp(image_ifft_real) result = my_ifill(image_ifp) if filetype == filter_type.laplace_pass: result = result / np.max(result) * 255 # 这里有一个问题,就是计算出的图像有没有负的像素 result = cv.convertScaleAbs(result) return result
进行频域滤波的一个关键步骤是确当恰当的截止频率,可以通过规定图像总功率$P_T$的圆来建立一组标准的截止频率轨迹。
低通滤波器
理想低通滤波器
在以原点为中心的一个圆内,无衰减地通过所有频率,而圆外的频率全部滤除掉的滤波器为理想低通滤波器。理想低通滤波器公式如下。
理想低通滤波器虽然在离散域容易实现,但是滤波后的图像会有明显的振铃现象,所以实际中不常使用。理想高通滤波器类似。
布特沃斯低通滤波器
n阶布特沃斯滤波器定义为
由于布特沃斯滤波器在截止频率处为平滑过渡,所以基本上没有振铃现象。
根据定义计算布特沃斯滤模板函数为
# 得到巴特沃斯滤波模板 def my_get_butterworth_model(size, d0, n=1): model = np.zeros(size, dtype=np.float64) wide, height = size for i in range(wide): for j in range(height): length = math.sqrt((i - wide / 2) ** 2 + (j - height / 2) ** 2) model[i][j] = 1 / (1 + pow(length / d0, 2 * n)) return model
高斯低通滤波器
二维高斯低通滤波器定义为
其中$D_0$是截止频率。高斯低通滤波器在截止频率出同样是平滑过渡,所以也没有振铃现象。
高斯低通滤波器模板函数为
# 得到高斯滤波模板 def my_get_gaussion_model(size, d0): model = np.zeros(size, dtype=np.float64) wide, height = size for i in range(wide): for j in range(height): length2 = (i - wide / 2) ** 2 + (j - height / 2) ** 2 model[i][j] = math.exp(-1 * length2 / (2 * d0 ** 2)) return model
布特沃斯高通滤波器
截止频率为$D_0$的n阶布特沃斯高通滤波器定义为
与布特沃斯低通滤波器相比$D_0$和$D(u,v)$的位置相反。模板函数为
# 得到巴特沃斯高通滤波模板 def my_get_butterworth_highpass_model(size, d0, n=1): model = np.zeros(size, dtype=np.float64) wide, height = size for i in range(wide): for j in range(height): length = math.sqrt((i - wide / 2) ** 2 + (j - height / 2) ** 2) + 0.00000001 model[i][j] = 1 / (1 + pow(d0 / length, 2 * n)) return model
高斯高通滤波器
高斯高通滤波器定义为
模板函数
# 得到高斯高通滤波模板 def my_get_gaussion_highpass_model(size, d0): model = np.zeros(size, dtype=np.float64) wide, height = size for i in range(wide): for j in range(height): length2 = (i - wide / 2) ** 2 + (j - height / 2) ** 2 model[i][j] = 1 - math.exp(-1 * length2 / (2 * d0 ** 2)) return model
频域Laplace滤波器
通过变换可以将空域的Laplace算子变换到频域,得到频域Laplace算子如下
在对将低频平移到中心的图像进行处理时,用以下公式计算
计算Laplace滤波模板的函数为
# 得到拉普拉斯滤波模板 def my_get_laplace_model(size): model = np.zeros(size, dtype=np.float64) wide, height = size for i in range(wide): for j in range(height): length2 = (i - wide / 2) ** 2 + (j - height / 2) ** 2 model[i][j] = -4*pow(math.pi, 2) * length2 return model
注意到图像经过Laplace滤波之后,得到的是图像的边缘信息,而且信息像素值有正有负。如果想要显示边缘信息,需要将得到的像素值取绝对值。如果想要通过边缘信息对图像进行锐化,则不需要取绝对值,直接用原图像减去边缘信息的 一定倍数即可。
进行Laplace锐化的函数为
# 拉普拉斯锐化 def myfunc_laplace_sharpen(image): laplace_result = myfunc_seldifine(image, my_get_laplace_model((2 * image.shape[0], 2 * image.shape[1]))) laplace_result = laplace_result / np.max(laplace_result) * 255 result = np.zeros(image.shape, dtype=np.int64) for i in range(image.shape[0]): for j in range(image.shape[1]): result[i][j] = int(image[i][j]) - int(laplace_result[i][j]) sharpen_image = cv.convertScaleAbs(result) return sharpen_image