图像滤波这话题听起来老生常谈但每次有人拿Python做图像处理/计算机视觉项目最终都会回到滤波器这三个字上来。我最近正好把均值滤波、低通滤波、高通滤波、逆滤波这四类在Python环境里重新捋了一遍边写代码边记笔记趁着热乎把整个思路、原理和可以直接跑的代码都整理出来。这篇内容适合图像处理入门到进阶之间那段最容易被卡住的阶段——你清楚滤波是个什么事但真正在numpy的API里操作时又觉得隔了一层。做图像处理的人应该都有这种感受滤波器的名字都认识均值滤波好像就是窗口平均低通滤掉高频、高通保留边缘逆滤波能做图像去模糊。但真到动手写代码很多人会卡在几个点上核怎么构造、频域怎么乘、傅里叶变换之后要不要fftshift、逆滤波出来的图为什么全是噪点。这次我不打算只丢几个现成API而是从原理到代码把每一类都拆开你照着跑一遍回头再看OpenCV里的各种滤波函数会顺畅很多。1. 项目定位与整体设计思路1.1 为什么一口气做四类滤波器一次做四类滤波器不是贪多而是它们刚好串起了图像处理的三条主线。均值滤波属于空间域模板操作直接在像素邻域上做加权平均最简单的卷积思想从这里起步。低通滤波和高通滤波属于频率域操作借助傅里叶变换把图像从空间域换到频率域然后对频谱的各个分量做乘法再把结果逆变换回来。逆滤波则属于图像复原范畴它处理的是已知退化模型、如何反推原图这类更高一层的问题。四条线放在一起你就能看出来均值滤波是空域卷积的直观案例频域低通/高通是卷积定理的应用逆滤波是频域除法在复原问题里的延伸。它们之间不是孤立的知识点而是从卷积→傅里叶变换→退化模型的完整认知链条。任何一个要做计算机视觉方向的开发者把这四类滤波器吃透后面再看边缘检测、图像增强、去模糊都会轻松很多。1.2 技术选型numpy手写为主OpenCV只做对比这套代码我尽量只用numpy和matplotlib核心逻辑全部手写。有人会问OpenCV一个cv2.blur、cv2.GaussianBlur就解决了何必手写我的想法是滤波器这类基础模块API调用越少越能看清原理。cv2.blur底层做了很多优化你直接调用感知不到核怎么滑过图像、边界如何处理、频域模板怎么构造。而手写一遍代码虽然啰嗦但每个步骤都摊开在面前。这篇内容里numpy负责数组运算和FFTmatplotlib负责可视化OpenCV只在速度对比和工程化建议里提到。1.3 核心概念预备卷积定理与频谱中心化在开始写代码前有两个概念必须拎清楚因为它们直接决定后续代码怎么写。第一个是卷积定理。空间域里的卷积对应频率域里的乘法。换句话说对一张图像做卷积操作等于先把图像通过傅里叶变换转到频率域让频谱乘以一个滤波器模板再逆变换回空间域。对计算机来说如果图像尺寸大、卷积核大直接在空间域滑窗计算非常慢而通过FFT到频域做乘法计算量大为降低。更重要的是许多滤波器在频域里描述非常简洁比如保留中心低频、抑制边缘高频这种操作在空间域很难直观设计在频域却只是一句半径判断的事。第二个是频谱中心化。numpy的fft2得到的结果四个角对应低频分量中心对应高频分量。但我们在设计频域滤波器时通常希望低频在图像中心这样才能用一个半径去判断保留还是滤除。所以标准流程是fft2之后调用fftshift把低频移到中心乘完滤波器模板后再调用ifftshift把频率恢复到原来的排列最后用ifft2回到空间域。这一步做错输出图会是颠倒或者错位的。2. 准备测试图像与通用工具函数2.1 合成一张测试图不依赖外部素材为了让你复制代码就能跑我不用现成的lena或cameraman图片而是自己生成一张合成测试图。合成图的好处是可以精确控制里面有什么——一个亮矩形、一个亮圆、一段渐变纹理、一个高频细节区域这样能直观看到滤波器对边缘和平坦区的处理差异。import numpy as np import matplotlib.pyplot as plt from numpy.fft import fft2, ifft2, fftshift, ifftshift def build_test_image(size256): img np.zeros((size, size), dtypenp.float64) # 背景渐变 yy, xx np.mgrid[0:size, 0:size] img 30 0.3 * (xx yy) / size * 255 # 亮矩形 img[70:110, 50:150] 160 # 亮圆 yy_c, xx_c np.mgrid[0:size, 0:size] circle (xx_c - 170) ** 2 (yy_c - 130) ** 2 35 ** 2 img[circle] 140 # 高频纹理区域棋盘格子 patch np.indices((40, 40)).sum(axis0) % 2 0 img[180:220, 180:220] patch * 80 return np.clip(img, 0, 255)这段代码生成一张256x256的灰度图里面包含明显的平坦区域、规则边缘和棋盘纹理。后面看滤波效果时主要盯三处矩形和圆的边缘是否保持、棋盘纹理是否被抹平、背景渐变是否出现波纹。2.2 频谱可视化先看看图像在频域长什么样做频域滤波最重要的直觉来源就是频谱图。我的习惯是每次跑滤波之前先把原图的频谱打印出来看一眼。频谱的动态范围非常大直接显示会变成一片黑中间一个白点所以要先取幅值再做log变换。def show_spectrum(img, titleSpectrum): spectrum np.log(1 np.abs(fftshift(fft2(img)))) plt.figure(figsize(6, 6)) plt.imshow(spectrum, cmapgray) plt.title(title) plt.axis(off) plt.show()原始图像里的大块平坦区域对应低频集中在频谱中心矩形和圆的边界属于高频表现为中心向外辐射的亮线棋盘格纹理则是与格子周期对应的离散亮点。你如果看到频谱中心有一个明显的十字或放射状亮纹说明原图里有强边缘后面做完低通滤波这些结构会被压暗。2.3 统一的结果处理函数滤波后的结果通常是float数组数值可能不在0到255之间尤其是高通滤波后会有负值。我准备了两个转换函数一个用于保持亮度的clip式转换另一个用于显示细节的归一化转换。def clip_to_uint8(img): return np.clip(np.real(img), 0, 255).astype(np.uint8) def normalize_to_uint8(img): img np.real(img) img img - img.min() if img.max() 0: img img / img.max() * 255 return img.astype(np.uint8)均值滤波和低通滤波的输出与原始图像亮度范围一致用clip_to_uint8。高通滤波和逆滤波的结果往往有负值且动态范围很宽如果不做min-max归一化直接取整会丢失大量细节。这个选择后面会反复用到。3. 均值滤波空域卷积的baseline3.1 均值滤波原理与核的角色均值滤波的思想一句话就能说完把每个像素替换成它邻域内所有像素的平均值。一个3x3均值核长这样1/9 * | 1 1 1 | | 1 1 1 | | 1 1 1 |这个核在图像上滑动每次取9个像素的平均。它本质上是把图像和一个矩形窗做卷积。矩形窗在频域里的响应是sinc函数主瓣让低频通过旁瓣让部分高频漏过去还会产生轻微振铃这是均值滤波比较原生的缺点。不过因为核小空域里直接算这个缺点通常不明显。均值滤波最基本的作用是降噪它把孤立噪声点通过邻域平均抹进周围像素里。代价是边缘也被平均了图像变模糊。核越大模糊越严重。3.2 手写均值滤波代码与边界处理下面这段代码用纯numpy实现重点是边界处理这里我用了reflect镜像填充。为什么不用zero padding如果你用全零填充图像四周会在滤波后出现一大圈黑边因为边缘像素被零拉低了亮度。用reflect模式相当于把边缘外的像素对称翻转进来视觉上自然很多。def mean_filter_manual(img, ksize3): kh kw ksize pad ksize // 2 img_pad np.pad(img, pad, modereflect) kernel np.ones((kh, kw), dtypenp.float64) / (kh * kw) out np.zeros_like(img, dtypenp.float64) for i in range(img.shape[0]): for j in range(img.shape[1]): patch img_pad[i:i kh, j:j kw] out[i, j] np.sum(patch * kernel) return np.clip(out, 0, 255).astype(np.uint8)要说明的是这个双重循环在256x256的图上跑起来还能接受但拿到大图就会很慢。真正做工程时我会用cv2.blur或cv2.filter2D它们底层有SIMD加速。手写循环的意义是看清卷积的滑窗过程以及核的每个位置如何影响输出。3.3 均值滤波的适用场景与频域视角实际项目里均值滤波很少单独作为图像增强手段因为它对边缘的破坏太直接。它更适合在管道最前端做简单的噪声抑制比如在阈值化之前让图像更平滑。如果想在降噪的同时保留边缘通常会上高斯滤波或双边滤波。从频域看3x3均值核对高频的抑制能力有限它压制了高频但并没有彻底切掉。如果你在图像上反复做多次均值滤波图像会变得越来越模糊这种迭代操作在数学上等价于一个更大的等效核但效率远不如直接做一个大核。4. 低通滤波频域平滑与振铃问题4.1 理想低通、巴特沃斯低通、高斯低通低通滤波的目标是保留低频分量、抑制高频分量在视觉上表现为平滑和模糊。但保留多少、抑制多少有不同的策略于是出现了三种常见滤波模板。理想低通最简单以频谱中心为圆心在半径之内的频率全部保留半径之外的频率全部置零。问题在于这个截止太硬频域里产生一个矩形窗逆变换回空间域后会因为sinc的旁瓣产生振铃图像边缘出现一圈圈明暗交替的光晕。巴特沃斯低通有一个可调阶数阶数越低过渡越平缓振铃越弱但高频残留也会多一些。高斯低通则完全没有振铃因为它在频域里是平滑的高斯形状逆变换回来仍然是高斯形状空域里不会产生负值旁瓣。实际项目中高斯低通最常用。4.2 频域低通滤波器的实现代码构造频域滤波器模板时我先用fftshift把低频搬到中心再根据到中心的距离判断每个位置该保留多少。def gaussian_lowpass(shape, cutoff30): rows, cols shape crow, ccol rows // 2, cols // 2 u np.arange(rows) - crow v np.arange(cols) - ccol U, V np.meshgrid(v, u) D np.sqrt(U ** 2 V ** 2) H np.exp(-(D ** 2) / (2 * cutoff ** 2)) return H def butterworth_lowpass(shape, cutoff30, order2): rows, cols shape crow, ccol rows // 2, cols // 2 u np.arange(rows) - crow v np.arange(cols) - ccol U, V np.meshgrid(v, u) D np.sqrt(U ** 2 V ** 2) H 1 / (1 (D / cutoff) ** (2 * order)) return H频域滤波的统一流程如下def apply_frequency_filter(img, H): F fft2(img) F_shift fftshift(F) G_shift F_shift * H F_back ifftshift(G_shift) filtered np.real(ifft2(F_back)) return filtered这段代码我会在低通、高通里共用。注意H必须和图像同尺寸并且H的中心位置要对应fftshift后F_shift的中心。H的最大值通常设计为1这样直流分量不受影响滤波后图像整体亮度不变。4.3 三种低通振铃对比为什么高斯最实用我在同一张测试图上分别用理想低通、二阶巴特沃斯、高斯低通跑了一遍截止半径都取40。理想低通的结果在矩形边缘附近出现了明显的明暗波动棋盘纹理区域甚至能看到环状光晕。巴特沃斯的过渡稍微缓和振铃减弱但并没有完全消失。高斯低通的表现最干净边缘被平滑了但周围没有任何人工波纹。有一个判断技巧如果滤波后的图像边缘出现了水波纹首先怀疑是不是用了理想低通。工程上我几乎不用理想低通它的硬截止在数学上好看视觉上往往得不偿失。滤波器类型过渡特性振铃风险常用场景理想低通陡峭截止严重教学演示巴特沃斯低通可控过渡中等频域设计实验高斯低通平滑过渡无图像平滑、特征预处理5. 高通滤波频域细节增强与边缘提取5.1 从1减去低通得到高通高通滤波和低通滤波是一对你甚至可以简单地把高通看成1减低通。因为频域模板H_lp在中心低频处接近1在边缘高频处接近0用1减去它中心变成0边缘变成1正好保留高频、压制低频。视觉效果就是边缘突出、平坦区域变暗。边缘在图像中对应像素值剧烈变化的位置也就是高频分量集中的地方。高通滤波把平坦区域压暗之后剩下的亮线基本都是轮廓。对测试图来说矩形边界、圆形边界、棋盘格纹理都会被保留。5.2 高通滤波实现代码def gaussian_highpass(shape, cutoff20): lowpass gaussian_lowpass(shape, cutoff) return 1 - lowpass def butterworth_highpass(shape, cutoff20, order2): lowpass butterworth_lowpass(shape, cutoff) return 1 - lowpass然后调用apply_frequency_filter。但这里有一个展示问题高通滤波后的图像平均值接近0直接clip会造成大量像素变黑所以必须用normalize_to_uint8做min-max归一化。这样才能看到暗背景里亮边缘的效果。5.3 高通结果的归一化与图像叠加如果想把高通结果用于锐化正确的做法不是直接显示高通图而是把高通分量叠加回原图。这个操作叫非锐化掩模公式是sharp img alpha * highpassalpha控制锐化程度一般取0.5到2之间。alpha太大噪声会被一起放大边缘周围出现白边。def unsharp_mask(img, cutoff20, alpha1.0): hp apply_frequency_filter(img, gaussian_highpass(img.shape, cutoff)) hp normalize_to_uint8(hp) sharp img alpha * hp return np.clip(sharp, 0, 255).astype(np.uint8)高通滤波结果对噪声非常敏感因为噪声本身也是高频分量。如果原图有噪声直接做高通会得到一堆雪花点。这种情况下要先做低通降噪再做高通提取边缘或者使用带通滤波只保留特定频段。6. 逆滤波退化模型与图像复原6.1 图像退化模型g h * f n逆滤波是图像复原方向的基础它处理的不是增强而是反推。假设原始清晰图像是f退化过程可以建模为g h * f n其中h是点扩散函数PSF也就是退化核可能是运动模糊、失焦模糊或者高斯模糊n是噪声。在频域里卷积变成乘法G(u,v) H(u,v) * F(u,v) N(u,v)于是最直观的复原思路是已知G和H将G除以H就得到F的估计值。psf make_gaussian_psf(21, 5) psf / psf.sum() H fft2(psf, simg.shape) G fft2(blurred) F_hat G / H restored np.abs(ifft2(F_hat))这里make_gaussian_psf生成一个二维高斯核然后归一化让它成为一个合法的模糊核。把原图卷积上这个核再加一点高斯噪声就得到退化图。原理不复杂但直接跑这个代码你会看到恢复出来的图像几乎全是噪点。6.2 为什么直接逆滤波会爆炸噪声放大问题逆滤波最大的坑是H在频域里有大量接近0的值。高斯PSF的频谱H随着频率升高快速衰减到高频段H的值已经接近0。此时噪声N还在G H*F N你除以一个接近0的H相当于把噪声放大了无数倍。高频段被噪声完全淹没恢复出来的图像就变成了一锅粥。解决思路是给逆滤波加一个截止半径只在H值还比较可靠的频段内做除法高频段不恢复直接置零。def inverse_filter_with_radius(blurred, psf, cutoff40): H fft2(psf, sblurred.shape) G fft2(blurred) rows, cols blurred.shape crow, ccol rows // 2, cols // 2 u np.arange(rows) - crow v np.arange(cols) - ccol U, V np.meshgrid(v, u) D np.sqrt(U ** 2 V ** 2) F_hat np.zeros_like(G, dtypenp.complex128) mask (D cutoff) (np.abs(H) 1e-6) F_hat[mask] G[mask] / H[mask] return np.real(ifft2(F_hat))cutoff取多少需要看H的频谱能量分布。一般来说退化核越大、模糊越厉害H的高频衰减越快有效截止半径就要取小一点。实际调试时我会把H的频谱打印出来肉眼看到哪个半径之外幅值已经低于1e-4就以那个为基准再收紧一点。6.3 维纳滤波与更稳的改进思路加半径截断是一种粗放的思路。更成熟的是维纳滤波它用信噪比信息对除法做了正则化F_hat conj(H) / (|H|^2 K) * G其中K是一个小的正常数用来抑制H接近0时的除法放大。K越大恢复结果越平滑但细节也会丢失。K越小越接近直接逆滤波噪声风险越大。def wiener_filter(blurred, psf, K0.01): H fft2(psf, sblurred.shape) conjs np.conj(H) denominator np.abs(H) ** 2 K F_hat conjs / denominator * fft2(blurred) return np.real(ifft2(F_hat))K的取值没有标准答案需要根据图像和噪声水平试。我的习惯是先取0.01看恢复图的噪声水平如果噪声太明显把K增加到0.05、0.1。如果你知道原图的噪声方差K也可以近似取噪声方差与信号方差的比值。7. 完整代码与运行结果串联7.1 一体化可复现脚本前面各节已经把函数拆开了这里给出一个完整的可运行脚本把均值、低通、高通、逆滤波全部串起来。把下面代码保存成filter_demo.py直接跑就行。import numpy as np import matplotlib.pyplot as plt from numpy.fft import fft2, ifft2, fftshift, ifftshift def build_test_image(size256): img np.zeros((size, size), dtypenp.float64) yy, xx np.mgrid[0:size, 0:size] img 30 0.3 * (xx yy) / size * 255 img[70:110, 50:150] 160 yy_c, xx_c np.mgrid[0:size, 0:size] circle (xx_c - 170) ** 2 (yy_c - 130) ** 2 35 ** 2 img[circle] 140 patch np.indices((40, 40)).sum(axis0) % 2 0 img[180:220, 180:220] patch * 80 return np.clip(img, 0, 255) def clip_to_uint8(img): return np.clip(np.real(img), 0, 255).astype(np.uint8) def normalize_to_uint8(img): img np.real(img) img img - img.min() if img.max() 0: img img / img.max() * 255 return img.astype(np.uint8) def mean_filter_manual(img, ksize3): pad ksize // 2 img_pad np.pad(img, pad, modereflect) kernel np.ones((ksize, ksize), dtypenp.float64) / (ksize * ksize) out np.zeros_like(img, dtypenp.float64) for i in range(img.shape[0]): for j in range(img.shape[1]): patch img_pad[i:i ksize, j:j ksize] out[i, j] np.sum(patch * kernel) return clip_to_uint8(out) def gaussian_lowpass(shape, cutoff30): rows, cols shape u np.arange(rows) - rows // 2 v np.arange(cols) - cols // 2 U, V np.meshgrid(v, u) D np.sqrt(U ** 2 V ** 2) return np.exp(-(D ** 2) / (2 * cutoff ** 2)) def gaussian_highpass(shape, cutoff20): return 1 - gaussian_lowpass(shape, cutoff) def apply_frequency_filter(img, H): F_shift fftshift(fft2(img)) F_back ifftshift(F_shift * H) return np.real(ifft2(F_back)) def make_gaussian_psf(size, sigma): ax np.linspace(-(size // 2), size // 2, size) g1d np.exp(-ax ** 2 / (2 * sigma ** 2)) g2d np.outer(g1d, g1d) return g2d / g2d.sum() def inverse_filter_with_radius(blurred, psf, cutoff40): H fft2(psf, sblurred.shape) G fft2(blurred) rows, cols blurred.shape u np.arange(rows) - rows // 2 v np.arange(cols) - cols // 2 U, V np.meshgrid(v, u) D np.sqrt(U ** 2 V ** 2) F_hat np.zeros_like(G, dtypenp.complex128) mask (D cutoff) (np.abs(H) 1e-6) F_hat[mask] G[mask] / H[mask] return normalize_to_uint8(np.real(ifft2(F_hat))) # 主流程 img build_test_image() mean_img mean_filter_manual(img, ksize5) lp_img apply_frequency_filter(img, gaussian_lowpass(img.shape, cutoff30)) hp_img apply_frequency_filter(img, gaussian_highpass(img.shape, cutoff20)) hp_img normalize_to_uint8(hp_img) psf make_gaussian_psf(21, 5) blurred np.real(ifft2(fft2(img) * fft2(psf, simg.shape))) blurred normalize_to_uint8(blurred np.random.normal(0, 8, img.shape)) inv_img inverse_filter_with_radius(blurred, psf, cutoff40) fig, axes plt.subplots(2, 3, figsize(14, 9)) axes[0, 0].imshow(img, cmapgray); axes[0, 0].set_title(Original) axes[0, 1].imshow(mean_img, cmapgray); axes[0, 1].set_title(Mean Filter) axes[0, 2].imshow(clip_to_uint8(lp_img), cmapgray); axes[0, 2].set_title(Gaussian Lowpass) axes[1, 0].imshow(hp_img, cmapgray); axes[1, 0].set_title(Gaussian Highpass) axes[1, 1].imshow(blurred, cmapgray); axes[1, 1].set_title(Blurred Noise) axes[1, 2].imshow(inv_img, cmapgray); axes[1, 2].set_title(Restored) for ax in axes.flat: ax.axis(off) plt.tight_layout() plt.show()7.2 参数速查与调优建议滤波器关键参数调参方向均值滤波ksize核越大噪声越少图越糊高斯低通cutoffcutoff越小平滑越强高斯高通cutoffcutoff越小保留的高频越多逆滤波截断cutoff半径越小越稳但细节丢失越多维纳滤波KK越大越平滑K越小越接近直接逆滤波调参时永远记得看频谱。低通滤波器的cutoff应该在频谱的能量分布图上找一个大概的半径而不是拍脑袋定值。高通阈值同理cutoff如果设得太大会把中频信息也滤掉图像只剩很细的轮廓线线条变碎。8. 高频踩坑与排查技巧实录8.1 滤波结果出现棋盘纹或密集条纹最典型的原因是频域模板和fftshift后的频谱没有对齐。比如你忘了对F做fftshift直接用中心化前后不一致的模板去乘结果就会乱。另一个常见原因是ifftshift用错方向最后的图会上下颠倒或者出现周期性的错位。排查方法很简单先把原图的频谱fftshift后的坐标和模板的值打印成图叠在一张图上肉眼对一下中心点是否重合。另外做完ifft2之后要取实部而不是直接取abs否则负值被强行镜像到正半轴会产生不自然的条纹。8.2 频谱图一片黑只有一个中心点这是log变换没有做好或者图像本身没有归一化。图像的频谱动态范围极大如果不加log直接显示边缘部分在0附近几乎不可见。建议使用np.log(1 np.abs(...))这里的1是防止log0出现。如果只看到一条竖直或水平的亮线说明原图有非常强的方向性结构比如棋盘格。这不是bug而是图像本身的频谱特征。8.3 逆滤波恢复结果全是噪点几乎可以确定是噪声放大问题。直接逆滤波在H接近0的频率处除以了极小值把噪声乘以了一个极大的放大系数。解决方式是加截止半径或者上维纳滤波。K值从0.01开始逐步增大。另外如果PSF方向不对比如运动模糊核没有对齐恢复也会出现放射状条纹。8.4 均值滤波大图跑起来太慢手写双重循环的复杂度是O(HWksize^2)图一大就卡到没法用。工程上建议直接用卷积定理或者OpenCV。真要在numpy里提速可以用np.lib.stride_tricks.sliding_window_view但理解成本比OpenCV更高。我的建议是学习时手写项目里用库函数。8.5 输出图像边界发黑或发白频域滤波本身不会引入边界问题但空域均值滤波会。前面说过zero padding会导致黑边用reflect模式能解决。另外np.clip之后可能让原本正常的像素溢出所以浮点结果在clip之前一定要确认范围。最后再分享一个小技巧跑滤波代码时我几乎每改一个参数都会把中间频谱图打出来。滤波器的效果、参数的选择单看输出图很难判断是该滤的滤掉了还是不该滤的也滤没了。只要把频谱图打开把低通或高通的半径可视化出来所有问题都会变成这个半径合不合适的直观判断。动手跑一遍比背熟任何参数都管用。这也是我这次整理代码最深的体会。