★ FEATURED ARTICLE 1. 从一张满是噪点的照片说起为什么滤波和空间变换是图像处理的必修课刚接触数字图像处理那会儿我拿到一张在暗光环境下拍的电路板照片满屏的椒盐噪声焊点边缘糊成一团连丝印字符都看不清。当时第一反应是调亮度、拉对比度结果噪声跟着一起被放大画面更糟。后来才明白图像处理里有两类操作是绕不开的基本功一类是滤波处理用来压制噪声、平滑细节或者反过来锐化边缘另一类是空间变换用来改变像素的空间位置关系比如旋转、缩放、平移、镜像。这两类操作几乎出现在所有后续任务的前置环节里——目标检测前要降噪医学影像配准前要几何校正遥感图像拼接前要做投影变换。Matlab在这个领域的位置很特殊。它不像OpenCV那样追求极致的运行效率但它的矩阵运算天然贴合图像的本质——一张灰度图就是一个二维矩阵彩色图就是三维张量。这意味着你可以用非常直观的方式去理解每一个像素级操作而不用被指针、内存对齐这些底层细节分散注意力。对于正在学数字图像处理课程、准备做课程设计或者想从理论过渡到动手实现的人来说Matlab是极佳的试验台。这篇内容会围绕图像的滤波处理和图像空间变换两条主线展开。滤波部分我会讲清楚空间域滤波和频率域滤波各自适合什么场景均值、中值、高斯、拉普拉斯这些算子到底在做什么以及怎么用Matlab把它们跑起来。空间变换部分会覆盖平移、旋转、缩放、镜像、仿射变换和透视变换重点解释插值方法的选择逻辑和变换矩阵的构造方式。每一块都会附上可直接运行的代码并且说明我在实际调试中踩过的坑和总结出的参数调节经验。不管你是刚学完理论课想找代码练手还是已经在做项目但效果不理想想回头补基础这些内容都能直接拿去用。2. 滤波处理从空间域到频率域的完整实现路径2.1 空间域滤波的底层逻辑与Matlab实现空间域滤波的本质很简单输出图像中每个像素的值由输入图像中该像素及其邻域内像素的值共同决定。这个共同决定的方式就是滤波核也叫卷积核、模板、窗口。核是一个小矩阵比如3×3、5×5里面的数值决定了邻域像素如何加权求和。先看均值滤波。它把邻域内所有像素的权重设为相同值然后取平均。3×3的均值核就是每个元素为1/9的矩阵。Matlab里可以用imfilter配合fspecial(average, [3 3])来实现也可以直接构造核。均值滤波的优点是实现简单、计算快缺点是会把边缘也一起模糊掉——因为边缘两侧的像素被同等对待过渡带被拉宽了。% 读取图像并转为灰度 img imread(circuit_board.jpg); if size(img, 3) 3 img_gray rgb2gray(img); else img_gray img; end % 添加椒盐噪声用于测试 img_noisy imnoise(img_gray, salt pepper, 0.05); % 均值滤波 h_mean fspecial(average, [3 3]); img_mean imfilter(img_noisy, h_mean, replicate); % 显示对比 figure; subplot(1,3,1); imshow(img_noisy); title(噪声图像); subplot(1,3,2); imshow(img_mean); title(均值滤波结果); subplot(1,3,3); imshow(img_gray); title(原始图像);这里有个细节值得注意imfilter的边界处理选项。默认是replicate也就是用边界像素的值填充外部区域。还有symmetric镜像填充、circular循环填充等。对于大多数自然图像symmetric通常比replicate更自然因为镜像填充不会在边界处引入突变。我实测下来如果图像边缘有重要内容边界选项的影响会很明显建议根据具体图像试一下不同选项。中值滤波是另一个必须掌握的算子。它不用加权求和而是取邻域内所有像素的中值作为输出。这个取中值的操作让它对椒盐噪声有奇效——因为椒盐噪声的像素值要么极大要么极小排序后会被挤到两端中值自然就把它们剔除了。Matlab里用medfilt2函数指定窗口大小即可。% 中值滤波窗口大小3x3 img_median medfilt2(img_noisy, [3 3]); % 对比不同窗口大小 img_median5 medfilt2(img_noisy, [5 5]); figure; subplot(1,3,1); imshow(img_noisy); title(噪声图像); subplot(1,3,2); imshow(img_median); title(3x3中值滤波); subplot(1,3,3); imshow(img_median5); title(5x5中值滤波);窗口大小的选择是个权衡。3×3能去掉大部分椒盐噪声同时保留较多细节5×5去噪更彻底但细小的纹理和边缘会被抹掉。我的经验是如果图像里有用信息集中在高频区域比如文字、细线优先用3×3如果是大面积平滑区域为主5×5甚至7×7都可以考虑。另外中值滤波对高斯噪声的效果不如均值滤波因为高斯噪声的像素值分布是连续的中值操作没有针对性。高斯滤波则是用高斯函数生成核中心权重最大向外逐渐衰减。它的好处是平滑效果自然不会像均值滤波那样产生明显的块状感。fspecial(gaussian, [h w], sigma)里的sigma控制衰减速度sigma越大核越矮胖平滑范围越广。% 高斯滤波不同sigma对比 h_gauss1 fspecial(gaussian, [5 5], 0.5); h_gauss2 fspecial(gaussian, [5 5], 1.5); h_gauss3 fspecial(gaussian, [5 5], 3.0); img_g1 imfilter(img_noisy, h_gauss1, symmetric); img_g2 imfilter(img_noisy, h_gauss2, symmetric); img_g3 imfilter(img_noisy, h_gauss3, symmetric); figure; subplot(1,3,1); imshow(img_g1); title(sigma0.5); subplot(1,3,2); imshow(img_g2); title(sigma1.5); subplot(1,3,3); imshow(img_g3); title(sigma3.0);sigma的选择直接决定平滑程度。sigma0.5时核几乎只对中心像素加权平滑效果很弱sigma1.5是比较常用的默认值sigma3.0时图像会明显变模糊。我一般会先用sigma1.0试一下如果噪声残留多就往上调如果细节丢失严重就往下调。还有一个技巧核大小一般取6*sigma1再向上取奇数这样能保证高斯函数在边界处已经衰减到接近零不会因为截断引入振铃。2.2 锐化滤波拉普拉斯算子与高频增强平滑滤波是压制高频锐化滤波则是增强高频。拉普拉斯算子是最常用的二阶微分算子它检测的是像素值的变化率。在Matlab里可以用fspecial(laplacian, alpha)生成alpha控制中心像素的权重。% 拉普拉斯锐化 h_lap fspecial(laplacian, 0.2); img_lap imfilter(double(img_gray), h_lap, replicate); % 将拉普拉斯结果叠加回原图实现锐化 img_sharp double(img_gray) - img_lap; img_sharp uint8(min(max(img_sharp, 0), 255)); figure; subplot(1,3,1); imshow(img_gray); title(原始图像); subplot(1,3,2); imshow(img_lap, []); title(拉普拉斯响应); subplot(1,3,3); imshow(img_sharp); title(锐化结果);这里有个关键点拉普拉斯响应有正有负直接imshow会看不到负值部分所以要用[]让Matlab自动缩放显示范围。另外锐化时是原图减去拉普拉斯还是原图加上拉普拉斯取决于核的符号约定。fspecial(laplacian)生成的核中心是负的所以要用减法。如果你自己构造核中心为正、周围为负那就用加法。这个符号问题我见过不少人搞反结果图像越锐化越模糊。锐化带来的副作用是噪声也会被放大。所以实际流程通常是先平滑再锐化或者用非锐化掩蔽unsharp masking先对图像做高斯模糊然后用原图减去模糊图得到高频分量再把高频分量按一定比例加回原图。% 非锐化掩蔽 img_blur imfilter(img_gray, fspecial(gaussian, [5 5], 1.0), symmetric); img_high double(img_gray) - double(img_blur); img_unsharp double(img_gray) 1.5 * img_high; img_unsharp uint8(min(max(img_unsharp, 0), 255)); figure; subplot(1,2,1); imshow(img_gray); title(原始图像); subplot(1,2,2); imshow(img_unsharp); title(非锐化掩蔽结果);那个1.5是增益系数控制锐化强度。系数越大锐化越猛但超过2.0之后容易出现光晕halo效应边缘两侧会出现亮暗条纹。我一般控制在1.2到1.8之间。2.3 频率域滤波傅里叶变换后的操作空间域滤波和频率域滤波是同一件事的两种视角。空间域的卷积对应频率域的乘法所以你可以把图像做傅里叶变换在频域乘以一个滤波器函数再逆变换回来。Matlab的fft2和ifft2让这个过程很直接。% 频率域低通滤波 img_double double(img_gray); F fft2(img_double); F_shifted fftshift(F); % 构造理想低通滤波器 [M, N] size(img_gray); D0 30; % 截止频率 u 0:M-1; v 0:N-1; [V, U] meshgrid(v, u); D sqrt((U - M/2).^2 (V - N/2).^2); H_lowpass double(D D0); % 滤波并逆变换 G F_shifted .* H_lowpass; img_lowpass real(ifft2(ifftshift(G))); img_lowpass uint8(min(max(img_lowpass, 0), 255)); figure; subplot(1,2,1); imshow(img_gray); title(原始图像); subplot(1,2,2); imshow(img_lowpass); title(理想低通滤波);理想低通滤波器在频域是一个硬截断这会导致空间域出现振铃效应ringing。原因是硬截断的逆变换是sinc函数sinc函数有旁瓣。要减轻振铃可以用巴特沃斯低通或高斯低通它们的过渡带是平滑的。% 巴特沃斯低通滤波 n 2; % 阶数 H_butter 1 ./ (1 (D ./ D0).^(2*n)); G_butter F_shifted .* H_butter; img_butter real(ifft2(ifftshift(G_butter))); img_butter uint8(min(max(img_butter, 0), 255)); % 高斯低通滤波 H_gauss exp(-(D.^2) / (2 * D0^2)); G_gauss F_shifted .* H_gauss; img_gauss_freq real(ifft2(ifftshift(G_gauss))); img_gauss_freq uint8(min(max(img_gauss_freq, 0), 255)); figure; subplot(1,3,1); imshow(img_lowpass); title(理想低通); subplot(1,3,2); imshow(img_butter); title(巴特沃斯低通); subplot(1,3,3); imshow(img_gauss_freq); title(高斯低通);巴特沃斯低通的阶数n控制过渡带的陡峭程度。n1时过渡很平缓n越大越接近理想滤波器振铃也越明显。n2是个不错的折中。高斯低通则完全没有振铃因为高斯函数的傅里叶变换还是高斯函数没有旁瓣。代价是过渡带比较宽截止频率附近的频率成分会被逐渐衰减而不是被干净地保留或去除。频率域滤波的优势在于可以精确控制频率成分。比如你要去除周期性噪声像扫描线、摩尔纹在频域找到对应的尖峰用一个陷波滤波器把它抠掉效果比空间域好得多。但频率域滤波的计算量更大而且对局部特征的适应性不如空间域。我的建议是常规去噪用空间域周期性噪声或需要精确频带控制时用频率域。3. 图像空间变换几何操作背后的矩阵与插值3.1 平移、旋转、缩放的矩阵表达与实现空间变换的核心是坐标映射输出图像中每个像素的位置对应输入图像中的哪个位置。这个映射可以用矩阵乘法表示。对于二维图像齐次坐标下的变换矩阵是3×3的。平移矩阵[1 0 tx] [0 1 ty] [0 0 1 ]旋转矩阵绕原点逆时针旋转θ[cosθ -sinθ 0] [sinθ cosθ 0] [0 0 1]缩放矩阵[sx 0 0] [0 sy 0] [0 0 1]Matlab的imtranslate、imrotate、imresize封装了这些操作但理解底层矩阵有助于你组合出更复杂的变换。% 平移 img_translated imtranslate(img_gray, [50, 30]); % 旋转crop表示裁剪到原图大小loose表示保留完整内容 img_rotated imrotate(img_gray, 30, bilinear, crop); % 缩放 img_scaled imresize(img_gray, 1.5, bilinear); figure; subplot(2,2,1); imshow(img_gray); title(原始); subplot(2,2,2); imshow(img_translated); title(平移[50,30]); subplot(2,2,3); imshow(img_rotated); title(旋转30度); subplot(2,2,4); imshow(img_scaled); title(放大1.5倍);旋转时的crop和loose选项很关键。crop保持输出尺寸和输入一致超出边界的部分被裁掉loose会扩大输出画布保证旋转后的图像完整。如果你后续要做图像配准或拼接通常用loose因为裁掉的信息可能正是你需要的。3.2 插值方法的选择最近邻、双线性与双三次空间变换必然涉及插值因为变换后的坐标通常不是整数。比如旋转30度后输出像素(100, 100)对应的输入坐标可能是(98.7, 103.2)你需要根据周围整数位置的像素值估算这个非整数位置的值。最近邻插值直接取最近的整数坐标像素值。速度最快但会产生锯齿边缘呈阶梯状。适合类别标签图、索引图这类不能做加权平均的数据。双线性插值用周围4个像素做两次线性插值。结果比最近邻平滑但会轻微模糊边缘。速度适中是最常用的默认选项。双三次插值用周围16个像素做三次插值。结果最平滑边缘保持最好但计算量最大且在边缘处可能产生轻微过冲。% 对比三种插值方法 img_nn imrotate(img_gray, 30, nearest, crop); img_bl imrotate(img_gray, 30, bilinear, crop); img_bc imrotate(img_gray, 30, bicubic, crop); figure; subplot(1,3,1); imshow(img_nn); title(最近邻); subplot(1,3,2); imshow(img_bl); title(双线性); subplot(1,3,3); imshow(img_bc); title(双三次);放大图像时插值方法的影响更明显。你可以把一张小图放大4倍然后放大局部区域观察。最近邻会看到明显的方块双线性是平滑的过渡双三次在边缘处更锐利。我的经验是如果后续要做边缘检测或分割用双三次如果只是显示或做粗略分析双线性足够如果处理的是掩膜或标签必须用最近邻否则会引入不存在的类别值。3.3 仿射变换与透视变换从矩阵构造到实际应用仿射变换是平移、旋转、缩放、错切的组合保持平行线仍然平行。透视变换则允许平行线在变换后不再平行更符合真实相机的成像模型。Matlab里可以用affine2d和projective2d构造变换对象然后用imwarp应用。% 仿射变换旋转缩放平移 theta 20; scale 1.2; tform_affine affine2d([scale*cosd(theta) -scale*sind(theta) 0; scale*sind(theta) scale*cosd(theta) 0; 0 0 1]); tform_affine.T(3,1:2) [30, -20]; % 平移量 img_affine imwarp(img_gray, tform_affine, OutputView, imref2d(size(img_gray))); % 透视变换用四对对应点估计 movingPoints [10 10; 200 10; 200 150; 10 150]; fixedPoints [30 20; 220 15; 210 170; 20 160]; tform_proj fitgeotrans(movingPoints, fixedPoints, projective); img_proj imwarp(img_gray, tform_proj, OutputView, imref2d(size(img_gray))); figure; subplot(1,3,1); imshow(img_gray); title(原始); subplot(1,3,2); imshow(img_affine); title(仿射变换); subplot(1,3,3); imshow(img_proj); title(透视变换);fitgeotrans的projective模式需要至少4对对应点。对应点的选择直接决定变换质量。如果点选得不准变换后的图像会扭曲。我通常会在图像中找明显的角点或标记点用cpselect工具交互式选取这样精度比手动估坐标高得多。透视变换的一个典型应用是文档矫正。拍一张倾斜的文档照片找到文档四个角在图像中的坐标再指定矫正后矩形的四个角坐标就能把文档拉平。这个流程在OCR预处理里非常常见。% 文档矫正示例 % 假设已找到文档四角坐标按顺时针 srcPoints [120 80; 480 95; 500 420; 100 400]; dstPoints [1 1; 400 1; 400 550; 1 550]; tform_doc fitgeotrans(srcPoints, dstPoints, projective); img_corrected imwarp(img_gray, tform_doc, OutputView, imref2d([550 400])); figure; subplot(1,2,1); imshow(img_gray); title(倾斜文档); subplot(1,2,2); imshow(img_corrected); title(矫正结果);这里imref2d指定了输出图像的空间范围。如果不指定imwarp会自动计算能容纳整个变换后图像的范围但有时候你希望输出固定尺寸就需要手动指定。4. 滤波与变换的组合实战一个完整的预处理流水线单独会用滤波和变换还不够实际项目里它们通常是组合出现的。我拿一个模拟场景来串一遍假设你有一批在传送带上拍摄的零件图像零件位置有偏移和旋转图像有噪声你需要把它们统一到标准姿态并降噪方便后续测量。第一步是去噪。传送带环境通常有椒盐噪声传感器坏点和高斯噪声光照不均。我的做法是先中值滤波去椒盐再高斯滤波去高斯噪声。顺序不能反——如果先高斯椒盐噪声会被扩散成一片中值滤波就不好处理了。% 第一步去噪 img_denoised medfilt2(img_noisy, [3 3]); img_denoised imfilter(img_denoised, fspecial(gaussian, [5 5], 1.0), symmetric);第二步是二值化和边缘提取用来找零件的姿态。可以用Otsu阈值分割然后做形态学操作去掉小噪点。% 第二步二值化与形态学清理 level graythresh(img_denoised); img_bw imbinarize(img_denoised, level); img_bw bwareaopen(img_bw, 50); % 去掉小于50像素的连通域 img_bw imclose(img_bw, strel(disk, 3)); % 闭运算填充小孔第三步是找零件的方向和位置。可以用regionprops获取连通域的主轴方向或者用最小外接矩形。% 第三步获取姿态参数 stats regionprops(img_bw, Orientation, Centroid, BoundingBox); if ~isempty(stats) angle stats(1).Orientation; centroid stats(1).Centroid; % 第四步旋转校正 img_corrected imrotate(img_denoised, -angle, bilinear, crop); % 第五步平移校正将质心移到图像中心 [M, N] size(img_corrected); tx N/2 - centroid(1); ty M/2 - centroid(2); img_corrected imtranslate(img_corrected, [tx, ty]); end这个流水线里每一步的参数都需要根据实际图像调整。比如bwareaopen的面积阈值如果零件本身很小阈值设大了会把零件也去掉imclose的圆盘半径如果零件有细小凹槽半径大了会把凹槽填平。我一般会先把中间结果都显示出来确认每一步的效果再往下走而不是一次性跑完再看最终结果。还有一个容易忽略的点旋转校正后图像边缘会出现黑色填充区域。如果后续要做全局统计比如平均灰度这些黑色区域会拉低均值。解决办法是在旋转后用掩膜标记有效区域统计时只计算掩膜内的像素。% 生成有效区域掩膜 mask true(size(img_denoised)); mask_rotated imrotate(mask, -angle, nearest, crop); mask_translated imtranslate(mask_rotated, [tx, ty]); % 统计时只计算有效区域 valid_pixels img_corrected(mask_translated); mean_gray mean(valid_pixels);这个掩膜思路在图像配准、拼接里也常用。核心思想是几何变换会引入无效区域后续处理必须知道哪些像素是可信的。5. 调试过程中那些文档不会告诉你的细节5.1 数据类型陷阱double与uint8的隐式转换Matlab图像处理里最常见的bug来源就是数据类型。imread读进来通常是uint8范围0-255。但很多滤波操作尤其是拉普拉斯、频域滤波会产生负值或超过255的值。如果你直接在uint8上做减法Matlab会饱和到0负值信息全丢了。% 错误做法uint8直接相减 a uint8(50); b uint8(80); c a - b; % 结果是0不是-30 % 正确做法先转double c double(a) - double(b); % 结果是-30我的习惯是只要涉及加减乘除或频域操作先把图像转成double算完再根据需要用uint8(min(max(x,0),255))转回去。虽然多了一步但能避免大量隐蔽的错误。5.2 滤波核归一化为什么你的滤波结果整体变暗了自定义滤波核时如果核内元素之和不为1滤波后的图像整体亮度会改变。比如一个3×3的核元素全是1和是9卷积后每个像素变成邻域和亮度变成原来的9倍直接饱和成白色。所以平滑类核必须归一化。% 自定义均值核并归一化 h ones(3, 3); h h / sum(h(:)); % 归一化和为1锐化核则通常不归一化因为它的和接近1中心正、周围负正负抵消。但如果和明显偏离1也需要调整。判断方法很简单用一个全白的图像做滤波如果输出还是全白或接近说明核的增益合适如果变灰或变黑就需要调整。5.3 旋转角度的正负与坐标系约定Matlab的imrotate正角度是逆时针旋转但regionprops返回的Orientation是主轴与水平轴的夹角范围-90到90度正负定义和imrotate不完全一致。我在这上面栽过跟头用regionprops得到角度后直接传给imrotate结果旋转方向反了。解决办法是先用几个已知角度的测试图像验证一下。比如把图像旋转30度用regionprops测角度看返回的是30还是-30然后确定校正时应该用angle还是-angle。这个验证花不了几分钟但能省掉后面大量调试时间。5.4 边界填充对滤波结果的影响imfilter的边界选项默认是replicate但不同选项在边界处的表现差异很大。如果图像边界有重要内容建议对比一下replicate、symmetric和circular的效果。% 对比边界选项 img_r imfilter(img_gray, h_gauss2, replicate); img_s imfilter(img_gray, h_gauss2, symmetric); img_c imfilter(img_gray, h_gauss2, circular); % 只看边界区域 border 20; subplot(1,3,1); imshow(img_r(1:border,:)); title(replicate边界); subplot(1,3,2); imshow(img_s(1:border,:)); title(symmetric边界); subplot(1,3,3); imshow(img_c(1:border,:)); title(circular边界);circular会把图像当成周期性的右边界和左边界相连。如果图像本身没有周期性这会在边界处引入不存在的突变。symmetric通常是最安全的选择因为它不会引入新的极值。5.5 频域滤波的频谱中心化与逆变换用fft2得到的频谱零频率在左上角。直接可视化会看到四个角是亮的中间是暗的不直观。所以要用fftshift把零频率移到中心。但滤波完做逆变换前必须用ifftshift移回去否则结果会错位。F fft2(double(img_gray)); F_shifted fftshift(F); % 零频率移到中心 % ... 在F_shifted上做滤波 ... G F_shifted .* H; img_result real(ifft2(ifftshift(G))); % 先移回再逆变换这个fftshift和ifftshift的配对使用是频域滤波的固定套路漏掉任何一个都会得到错误结果。而且ifft2的结果可能有微小虚部浮点误差所以取real是必要的。6. 从这些基础操作延伸出去的方向把滤波和空间变换跑通之后你会发现很多更复杂的任务都是它们的组合。图像金字塔高斯金字塔、拉普拉斯金字塔本质上是缩放加滤波Canny边缘检测里包含了高斯滤波、梯度计算、非极大值抑制SIFT特征点检测里有多尺度空间构造也是高斯滤波加缩放。甚至深度学习里的卷积层本质上就是可学习的滤波核。如果你想继续深入我建议从两个方向选一个。一个是图像配准把不同时间、不同视角拍的同一场景对齐核心是空间变换加相似度度量。另一个是图像增强在滤波的基础上做自适应处理比如直方图均衡化、Retinex增强让图像在视觉上更清晰。Matlab的Image Processing Toolbox里还有很多现成函数可以探索比如imguidedfilter引导滤波、imdiffusefilt各向异性扩散滤波、imregister自动配准。但前提是你对基础滤波和变换的底层逻辑有清晰的认识否则调参时只能瞎试不知道每个参数在控制什么。我在实际项目里最深的体会是不要一上来就追求复杂算法。把均值、中值、高斯、拉普拉斯这几个滤波器的特性摸透把平移、旋转、缩放、仿射、透视的矩阵构造和插值方法搞清楚后面遇到再复杂的需求你都能把它拆解成这些基础操作的组合。基础扎实了调参就是有方向的微调而不是碰运气。
阅读完成 · 觉得有帮助？
