简介傅里叶梅林变换参考MATLAB代码是一份面向图像处理与计算机视觉初学者的示例程序聚焦利用傅里叶与梅林变换的结合实现旋转、平移不变性分析适合学习图像配准、OCR或目标识别时参考。压缩包共6个文件包括1个Register.m源码脚本和5幅bmp格式的Lena测试图像整体仅152KB轻量易用。测试图像覆盖裁剪、旋转、平移及其组合场景可直接运行查看配准效果。已有3982人学习下载说明该主题关注度较高。代码中未实现缩放部分读者可在此基础上自行扩展尺度估计并通过修改Register.m中的参数加深对频域变换流程的理解是入门傅里叶梅林变换的一份实用范例。 做图像配准的时候最让人头疼的往往不是算法本身的取舍而是你手里那张图根本提不出特征点。我之前处理过一批晶圆表面光学照片纹理少得可怜整张图就是灰底上几道圆弧SIFT跑完一个关键点都没有。后来换用傅里梅林变换也就是大家常说的Fourier-Mellin Transform不需要提特征直接在频域把旋转、缩放、平移三个参数一次估计出来整个流程几个小时就写完。这篇我把整理好的参考MATLAB代码完整放出来加上我实际跑过的参数细节和踩坑记录代码可以直接复制改改就用适合图像配准、目标对齐、工业视觉定位这几类场景的读者参考。Fourier-Mellin变换这个名字听着唬人但它的核心并不是某个像FFT一样的新变换而是一套经典组合拳对两幅图像做FFT取幅度谱把幅度谱映射到对数极坐标空间再在这个空间里做一次相位相关最后把估计出的旋转缩放参数代回去做平移精配准。整个过程里傅里叶变换负责处理平移与旋转缩放梅林变换对应的是对数映射这一步两两组合配合相位相关就可以在完全不依赖特征点的情况下把刚性变换参数全部解出来。1. 为什么需要傅里叶梅林变换旋转缩放平移配准的频域解法先说一个反直觉的结论对于某些图像频域配准比空域匹配更稳。很多人在做图像配准时第一反应是特征点匹配SIFT、ORB、SuperPoint提完特征再做特征匹配和单应估计这套流程在纹理丰富的自然图像上确实好用。但一旦遇到低纹理工业件、重复纹理的板材、医学显微镜图片或者是两张只有灰度差异的遥感切片特征点的数量和质量会断崖式下降后续RANSAC都救不回来。傅里叶梅林变换解决的正是这种场景它不检测任何角点或斑点而是把整幅图的全局频率结构作为匹配依据。也就是说只要图像整体上有可辨识的纹理周期性或频谱分布哪怕像素级特征很少也能把旋转、缩放、平移参数估计出来。我当年就是被晶圆图虐完之后彻底转向这套频域思路。整套方法适合谁用概括下来有三类做模板匹配或目标对齐但目标图像纹理单一、没有稳定特征点已知两幅图之间只存在旋转、缩放、平移的刚性变换不需要处理仿射或透视需要快速获得初始配准参数后续再配合精配准或迭代优化。它的设计逻辑也非常清晰不追杀特征的绝对位置而是研究两张图在频率结构上的差异。如果图像A是图像B旋转了某个角度得到的那么它们的频谱幅度图也相差同样的旋转角度如果图像A是图像B缩放得到的那么频谱幅度图则是反向缩放。只要把这两个关系量化旋转和缩放参数就出来了。平移参数则在旋转缩放宽校正之后直接在空间域再做一次普通相位相关一次FFT就搞定。2. 核心原理拆解从FFT幅度谱到对数极坐标2.1 傅里叶变换的三个关键性质傅里叶变换在配准领域立得住靠的是三条经典性质空域平移只改变频谱的相位不改变幅度。图像整体平移(tx, ty)频谱F(u,v)变为F(u,v)乘以相位因子e^{-j2π(utx/Mvty/N)}幅度谱不变。这意味着旋转、缩放估计完全可以只看幅度谱平移不影响它们。空域旋转等价于频谱同角度旋转。图像旋转θ度频谱幅度图也旋转θ度。可以类比成你拿一个印花圆盘整个圆盘转多少盘上的纹理方向就转多少。空域缩放等价于频谱反向缩放。图像放大a倍频谱在频率轴上压缩为原来的1/a幅度整体还要乘一个1/a²因子。换句话说图像变大频谱变小图像变小频谱变大。这三条合在一起旋转和缩放的估计就可以从幅度谱入手了。但是幅度谱的坐标系是直角坐标旋转和缩放作用在上面仍然很复杂旋转让幅度谱的坐标旋转缩放让幅度谱的坐标伸缩两者纠缠在一起没法直接解。2.2 对数极坐标把缩放问题变成平移问题接下来就是对数极坐标变换发挥作用的地方。把幅度谱从直角坐标(u,v)映射到(logρ, θ)坐标横坐标是频率半径的对数纵坐标是0到360度的角度。这一步的效果非常巧妙图像旋转θ₀度幅度谱的每一个点沿着角度方向移动θ₀对应到对数极坐标空间就是上下平移图像缩放a倍幅度谱的半径坐标变化为原来的1/a由于横坐标取了log这个“缩小”的倍数就变成了一个固定偏移log(1/a)对应到对数极坐标空间就是左右平移。所以原来旋转变成了角度轴的平移缩缩放变成了半径对数轴的平移。一个二维旋转加缩放估计问题被转化成对数极坐标空间里的纯平移估计问题。这背后其实就是梅林变换的经典性质它天然把尺度变化转变成平移量。这里要注意一个细节幅度谱是中心对称的也就是|F(u,v)||F(-u,-v)|对一整圈360度来说间隔180度的频谱是重复的。这带来一个非常实际的坑我放在后面专门说。正因为有这个对称性单纯靠幅度谱只能估计出0到180度的角度180度方向存在原理性歧义。2.3 相位相关整数像素到亚像素能力平移量怎么测相位相关是标准答案。它做的事情是计算两幅图频谱的归一化互功率谱再做逆FFT变换。数学形式就是P(u,v) F₁(u,v) · conj(F₂(u,v)) / |F₁(u,v) · conj(F₂(u,v))|逆变换后的空域结果会在平移量处出现一个尖锐的峰。由于分母做了归一化只保留相位信息、丢掉幅度信息所以相位相关对图像的整体亮度变化、对比度变化都不敏感。相位相关的峰越尖说明匹配越可靠。峰的高度可以直接当置信度用接近1表示两图相关性很强低于0.2基本可以判定配准失败。亚像素精度可以通过对峰邻域做抛物线拟合获得但如果不做整数像素精度对大多数工业预配准场景已经够用。3. MATLAB实现可直接运行的参考代码Fourier-Mellin的MATLAB实现网上碎片很多但完整可跑的并不多。我这个版本把主流程、对数极坐标采样、相位相关全部封装成了函数主入口是fm_register输入参考图和待配准图输出旋转角度、缩放因子、平移量和置信度。3.1 主函数与整体流程function [theta, scale, delta, peak] fm_register(im1, im2, opts) % FM_REGISTER 傅里叶梅林变换图像配准 % 输入 % im1 参考图像灰度图double矩阵 % im2 待配准图像灰度图double矩阵尺寸需与im1一致 % opts 可选参数详见代码 % 输出 % theta 旋转角度单位度im2相对im1的旋转角 % scale 缩放因子im2相对im1的缩放1表示放大 % delta [dy, dx]im2相对im1的平移量 % peak 相位相关置信度越接近1越可靠 arguments im1 double im2 double opts.numR (1,1) double 128 opts.numTheta (1,1) double 256 opts.useWindow (1,1) logical true opts.useHPF (1,1) logical true opts.fineMode (1,1) logical true end % 1. 预处理灰度化、去均值、归一化 im1 preprocess(im1); im2 preprocess(im2); % 2. 计算频谱幅度谱 spec1 spectral_amplitude(im1, opts); spec2 spectral_amplitude(im2, opts); % 3. 对数极坐标变换 lp1 log_polar_transform(spec1, opts.numR, opts.numTheta); lp2 log_polar_transform(spec2, opts.numR, opts.numTheta); % 4. 对数极坐标空间相位相关 [shift_lp, peak] phase_correlation(lp1, lp2); dy shift_lp(1); dx shift_lp(2); % 5. 换算旋转角度与缩放因子 theta dy / opts.numTheta * 360; if theta 180 theta theta - 360; end maxR min(floor((size(im1,1)-1)/2), floor((size(im1,2)-1)/2)); logR_step log(maxR) / (opts.numR - 1); scale exp(-dx * logR_step); % 6. 消除180度歧义并估计平移量 if opts.fineMode bestPeak -inf; for k 1:2 candTheta theta 180 * (k - 1); im2c correct_rs(im2, candTheta, scale); [delta, pk] phase_correlation(im1, im2c); if pk bestPeak bestPeak pk; theta candTheta; delta delta; peak pk; end end else im2c correct_rs(im2, theta, scale); [delta, ~] phase_correlation(im1, im2c); end end主流程就是这样。有几个点值得说明一下。第一步的预处理里去均值和归一化是必须的。如果不去均值直流分量会在频谱中心形成一个极大亮点对数极坐标采样时会严重压制周围刻度信息。归一化保证不同亮度的图像处于同一量纲否则频谱幅度差异过大后续高通滤波参数会失准。第六步的fineMode是我强烈建议保留的。前面说过幅度谱的180度对称性导致角度估计存在歧义直接用单一角度校正可能在一个完全错误的方向上完成平移估计。加了双候选之后代码会用theta和theta180两个角度分别校正图像再做一次空间域相位相关取置信度更高的那个。这个思路用一次FFT成本换回角度方向上的确定性非常划算。3.2 对数极坐标变换的采样细节对数极坐标变换是整套代码里最影响精度的部分。我这里的实现是对幅度谱做半径和角度的网格采样然后双三次插值取值。function lp log_polar_transform(amp, numR, numTheta) % 对数极坐标变换 % amp: fftshift后的频谱幅度谱 % numR: 半径采样数 % numTheta: 角度采样数 [M, N] size(amp); cx (N 1) / 2; cy (M 1) / 2; maxR min(floor((M - 1) / 2), floor((N - 1) / 2)); % 半径从1到maxR对数均匀采样 radius exp(linspace(log(1), log(maxR), numR)); % 角度从0到2*pi去掉最后一个重复采样点 theta linspace(0, 2*pi*(1 - 1/numTheta), numTheta); [Rgrid, Tgrid] ndgrid(radius, theta); Xq cx Rgrid .* cos(Tgrid); Yq cy Rgrid .* sin(Tgrid); lp interp2(amp, Xq, Yq, cubic, 0); end半径为什么从1开始而不是0因为log0没有定义。半径0这一个点正好是频谱直流分量的位置本来在去均值之后它就已经被削弱没必要强行采进来。而且半径过小的区域离散化严重一两个像素的变化就会让对数半径抖动很大采进来只会增加噪声。插值方式我默认用cubic实测比linear精度高。尤其是半径较小时线性插值会让频谱的圆环纹理变糊导致相位相关峰不够尖。如果你的MATLAB版本没有cubic选项可以用spline效果相当但速度会慢一些。还有一个细节lp最后做了一个转置。这是因为ndgrid生成的矩阵形状是numR行numTheta列而我习惯把输出整理成“行对应角度、列对应半径”的格式。这样后面的相位相关在行方向的位移直接对应角度变化列方向的位移直接对应log半径变化调试的时候人脑不容易绕晕。3.3 相位相关与坐标环绕处理相位相关的实现不长但有一个坐标处理很容易写错。function [delta, peak] phase_correlation(im1, im2) F1 fft2(im1); F2 fft2(im2); C F1 .* conj(F2); C C ./ (abs(C) eps); r ifft2(C); % 将零频移到中心使平移量以图像中心为原点 r ifftshift(r); [peak, idx] max(r(:)); [dy, dx] ind2sub(size(r), idx); % 处理环绕平移超过半幅的按负数处理 if dy size(r, 1) / 2 dy dy - size(r, 1); end if dx size(r, 2) / 2 dx dx - size(r, 2); end delta [dy, dx]; end这里最容易踩的坑是ifftshift和fftshift的区别。fft2之后零频在矩阵的(1,1)位置逆变换得到的结果也需要做同样的移位才方便直接当坐标读。偶数尺寸下fftshift和ifftshift等价但奇数尺寸下两者差一个像素所以统一用ifftshift更稳。分母加eps是防止除零。这个细节不加的话频谱中某些位置幅度恰好为0时会出现NaN相位相关结果直接废掉。我用的是eps而不是0既能防除零又不会像加一个固定大数那样污染归一化结果。translation部分的环绕处理也很重要。相位相关的峰可能出现在矩阵边界附近例如真实平移量是[-5, -3]但峰的原始位置可能记录为[size-5, size-3]。如果不把超过半幅的坐标折算成负数后面的校正环节会偏移一大截。4. 实验验证与参数调优4.1 构造已知参数的测试图像写算法最怕没有对照真值。我习惯先合成一个已知参数的测试集把图像做一次确定的旋转、缩放、平移然后用算法去反推。这个验证通过了再上真实图像。%% 生成测试图像 rng(42); N 256; [X, Y] meshgrid(1:N, 1:N); im1 0.6 0.3*sin(X/8) 0.25*cos(Y/6) 0.2*sin((XY)/4 1); im1 im1 0.1*randn(N, N); im1 imgaussfilt(im1, 1); im1 mat2gray(im1); %% 设定真实变换参数 theta_true 25; % 逆时针25度 scale_true 1.25; % 放大1.25倍 delta_true [12, -9]; % 平移量 %% 扩边生成待配准图 pad 160; im1p padarray(im1, [pad pad], mean(im1(:)), both); im2p imrotate(im1p, theta_true, bilinear, crop); im2p imresize(im2p, scale_true, bicubic); im2 crop_center(im2p, N, N); im2 circshift(im2, delta_true); %% 配准 [theta_est, scale_est, delta_est, peak] fm_register(im1, im2); fprintf(theta truth/est: %.2f / %.2f\n, theta_true, theta_est); fprintf(scale truth/est: %.4f / %.4f\n, scale_true, scale_est); fprintf(delta truth/est: %d %d / %d %d\n, ... delta_true(1), delta_true(2), delta_est(1), delta_est(2)); fprintf(peak: %.3f\n, peak);这里有一个重要技巧生成待配准图时先padarray扩边再做旋转和缩放最后中心裁剪回原始尺寸。如果不扩边imrotate的crop选项会直接裁掉图像旋转后露出的四角等于是人为制造了信息丢失会让算法性能看起来比实际更差。真实项目中如果已知图像内容大概率超出画面范围也建议先扩边再估计。4.2 典型输出与180度歧义处理在我机器上跑一组典型输出大概是这样theta truth/est: 25.00 / 24.84 scale truth/est: 1.2500 / 1.2493 delta truth/est: 12 -9 / 12 -9 peak: 0.912角度误差在1度以内缩放误差千分之一量级平移量整数像素完全命中。这个精度对预配准来说是完全够用的。不同版本MATLAB和不同插值参数下数字会有小幅波动但整体数量级应该一致。测试时建议专门试一下角度为0度到180度之间的多组值特别注意接近0度和接近180度的边界情况。如果角度为0缩放估计仍然正确如果角度为180度最后双候选逻辑会稳定选择峰值更大的180度方向。我曾经在只做单候选角度估计时遇到过把实际旋转170度的图像配成-10度的情况原因就是180度歧义没有消除所以fineMode一定要开。4.3 参数调整速查表代码里开放了几个参数实际使用时可以按需调整参数默认值作用调整建议numTheta256角度采样数决定角度分辨率需要更高角度精度时加到512耗时线性增加numR128半径采样数决定缩放分辨率缩放范围大时加到192或256useWindowtrue是否加汉宁窗抑制频谱泄漏图像边缘本来就是平滑过渡时关闭useHPFtrue是否启用高通斜坡滤波低纹理图像建议保持开启fineModetrue是否做180度双候选消歧强烈建议保持开启numTheta和numR不是越大越好。采样数翻倍对数极坐标变换的时间也翻倍但精度提升在numTheta达到512之后会大幅放缓。对多数场景256个角度采样已经把分辨率压到1.4度以内配合亚像素抛物线拟合足够。5. 一些容易踩的坑这套代码我前前后后跑了挺多真实图像有几个坑属于那种你不踩一次根本不会注意的地方。第一个坑是角度符号的方向问题。图像坐标系的y轴向下MATLAB的imrotate正角度对应视觉上的逆时针但不同数据形态下“偏移方向”和“频谱旋转方向”的相对关系可能让你估计出的角度恒为负数。我自己的习惯是先跑一遍已知角度的合成测试如果发现角度符号反了在theta输出处加一个负号就行。如果你用OpenCV或其他库生成测试数据更要在实验里先验证方向不要直接照搬。第二个坑是180度歧义。这个在很多中文资料里都没有明说但原理上非常关键。幅度谱满足共轭对称性意味着转180度之后幅度谱不变。只靠幅度谱做角度估计天然无法区分theta和theta180。解决思路就是我代码里的双候选候选角度分别校正后做空间域相位相关取置信度更高的方向。这个方案简单有效但要求你最终一定做空间域平移精配准不能只输出LP域的估计结果就结束。第三个坑是高通滤波的强度。自然图像的频谱能量高度集中在低频如果不对幅度谱做任何处理对数极坐标里几乎只有一两个低频采样点有能量相位相关峰会钝得没法看。我用的斜坡滤波是在0.15到1.0之间随频率线性增强目的就是给高频频谱一个相对公平的参与机会。但这个强度不是越大越好如果图像本身高频信息较弱强行增强会让噪声主导频谱peak值会掉得很快。第四个坑是图像尺寸。Fourier-Mellin的输入图像最好保持正方形或者长宽比不要太极端而且尺寸最好是偶数。奇数尺寸下频谱中心和图像中心的对应关系会差半个像素虽然对粗配准影响不大但在做亚像素精度需求的项目里会莫名其妙地出现系统性偏差。第五个坑是缩放范围。我的代码里半径采样上限是min(M,N)/2如果待配准图相对参考图缩得太小或放得太大频谱的主要能量会落到采样范围之外估计结果就会失效。这种场景的常规解法是先做一层金字塔粗略估计放大或者缩小回到合理范围再用傅里叶梅林精估计。这个方法我后面有空会单独写一篇专门讲。最后想说的是这套方法不怕光照变化不怕低纹理但怕两幅图不是同一尺度内容。比如一张图是汽车全景另一张图是车身局部这种不是全局缩放关系的情况Fourier-Mellin本身就不适用。使用之前先确认你的问题确实是刚性变换这是省时间的最大前提。本文还有配套的精品资源点击获取 傅里叶梅林变换图像配准实战：从原理到MATLAB代码实现
简介傅里叶梅林变换参考MATLAB代码是一份面向图像处理与计算机视觉初学者的示例程序聚焦利用傅里叶与梅林变换的结合实现旋转、平移不变性分析适合学习图像配准、OCR或目标识别时参考。压缩包共6个文件包括1个Register.m源码脚本和5幅bmp格式的Lena测试图像整体仅152KB轻量易用。测试图像覆盖裁剪、旋转、平移及其组合场景可直接运行查看配准效果。已有3982人学习下载说明该主题关注度较高。代码中未实现缩放部分读者可在此基础上自行扩展尺度估计并通过修改Register.m中的参数加深对频域变换流程的理解是入门傅里叶梅林变换的一份实用范例。 做图像配准的时候最让人头疼的往往不是算法本身的取舍而是你手里那张图根本提不出特征点。我之前处理过一批晶圆表面光学照片纹理少得可怜整张图就是灰底上几道圆弧SIFT跑完一个关键点都没有。后来换用傅里梅林变换也就是大家常说的Fourier-Mellin Transform不需要提特征直接在频域把旋转、缩放、平移三个参数一次估计出来整个流程几个小时就写完。这篇我把整理好的参考MATLAB代码完整放出来加上我实际跑过的参数细节和踩坑记录代码可以直接复制改改就用适合图像配准、目标对齐、工业视觉定位这几类场景的读者参考。Fourier-Mellin变换这个名字听着唬人但它的核心并不是某个像FFT一样的新变换而是一套经典组合拳对两幅图像做FFT取幅度谱把幅度谱映射到对数极坐标空间再在这个空间里做一次相位相关最后把估计出的旋转缩放参数代回去做平移精配准。整个过程里傅里叶变换负责处理平移与旋转缩放梅林变换对应的是对数映射这一步两两组合配合相位相关就可以在完全不依赖特征点的情况下把刚性变换参数全部解出来。1. 为什么需要傅里叶梅林变换旋转缩放平移配准的频域解法先说一个反直觉的结论对于某些图像频域配准比空域匹配更稳。很多人在做图像配准时第一反应是特征点匹配SIFT、ORB、SuperPoint提完特征再做特征匹配和单应估计这套流程在纹理丰富的自然图像上确实好用。但一旦遇到低纹理工业件、重复纹理的板材、医学显微镜图片或者是两张只有灰度差异的遥感切片特征点的数量和质量会断崖式下降后续RANSAC都救不回来。傅里叶梅林变换解决的正是这种场景它不检测任何角点或斑点而是把整幅图的全局频率结构作为匹配依据。也就是说只要图像整体上有可辨识的纹理周期性或频谱分布哪怕像素级特征很少也能把旋转、缩放、平移参数估计出来。我当年就是被晶圆图虐完之后彻底转向这套频域思路。整套方法适合谁用概括下来有三类做模板匹配或目标对齐但目标图像纹理单一、没有稳定特征点已知两幅图之间只存在旋转、缩放、平移的刚性变换不需要处理仿射或透视需要快速获得初始配准参数后续再配合精配准或迭代优化。它的设计逻辑也非常清晰不追杀特征的绝对位置而是研究两张图在频率结构上的差异。如果图像A是图像B旋转了某个角度得到的那么它们的频谱幅度图也相差同样的旋转角度如果图像A是图像B缩放得到的那么频谱幅度图则是反向缩放。只要把这两个关系量化旋转和缩放参数就出来了。平移参数则在旋转缩放宽校正之后直接在空间域再做一次普通相位相关一次FFT就搞定。2. 核心原理拆解从FFT幅度谱到对数极坐标2.1 傅里叶变换的三个关键性质傅里叶变换在配准领域立得住靠的是三条经典性质空域平移只改变频谱的相位不改变幅度。图像整体平移(tx, ty)频谱F(u,v)变为F(u,v)乘以相位因子e^{-j2π(utx/Mvty/N)}幅度谱不变。这意味着旋转、缩放估计完全可以只看幅度谱平移不影响它们。空域旋转等价于频谱同角度旋转。图像旋转θ度频谱幅度图也旋转θ度。可以类比成你拿一个印花圆盘整个圆盘转多少盘上的纹理方向就转多少。空域缩放等价于频谱反向缩放。图像放大a倍频谱在频率轴上压缩为原来的1/a幅度整体还要乘一个1/a²因子。换句话说图像变大频谱变小图像变小频谱变大。这三条合在一起旋转和缩放的估计就可以从幅度谱入手了。但是幅度谱的坐标系是直角坐标旋转和缩放作用在上面仍然很复杂旋转让幅度谱的坐标旋转缩放让幅度谱的坐标伸缩两者纠缠在一起没法直接解。2.2 对数极坐标把缩放问题变成平移问题接下来就是对数极坐标变换发挥作用的地方。把幅度谱从直角坐标(u,v)映射到(logρ, θ)坐标横坐标是频率半径的对数纵坐标是0到360度的角度。这一步的效果非常巧妙图像旋转θ₀度幅度谱的每一个点沿着角度方向移动θ₀对应到对数极坐标空间就是上下平移图像缩放a倍幅度谱的半径坐标变化为原来的1/a由于横坐标取了log这个“缩小”的倍数就变成了一个固定偏移log(1/a)对应到对数极坐标空间就是左右平移。所以原来旋转变成了角度轴的平移缩缩放变成了半径对数轴的平移。一个二维旋转加缩放估计问题被转化成对数极坐标空间里的纯平移估计问题。这背后其实就是梅林变换的经典性质它天然把尺度变化转变成平移量。这里要注意一个细节幅度谱是中心对称的也就是|F(u,v)||F(-u,-v)|对一整圈360度来说间隔180度的频谱是重复的。这带来一个非常实际的坑我放在后面专门说。正因为有这个对称性单纯靠幅度谱只能估计出0到180度的角度180度方向存在原理性歧义。2.3 相位相关整数像素到亚像素能力平移量怎么测相位相关是标准答案。它做的事情是计算两幅图频谱的归一化互功率谱再做逆FFT变换。数学形式就是P(u,v) F₁(u,v) · conj(F₂(u,v)) / |F₁(u,v) · conj(F₂(u,v))|逆变换后的空域结果会在平移量处出现一个尖锐的峰。由于分母做了归一化只保留相位信息、丢掉幅度信息所以相位相关对图像的整体亮度变化、对比度变化都不敏感。相位相关的峰越尖说明匹配越可靠。峰的高度可以直接当置信度用接近1表示两图相关性很强低于0.2基本可以判定配准失败。亚像素精度可以通过对峰邻域做抛物线拟合获得但如果不做整数像素精度对大多数工业预配准场景已经够用。3. MATLAB实现可直接运行的参考代码Fourier-Mellin的MATLAB实现网上碎片很多但完整可跑的并不多。我这个版本把主流程、对数极坐标采样、相位相关全部封装成了函数主入口是fm_register输入参考图和待配准图输出旋转角度、缩放因子、平移量和置信度。3.1 主函数与整体流程function [theta, scale, delta, peak] fm_register(im1, im2, opts) % FM_REGISTER 傅里叶梅林变换图像配准 % 输入 % im1 参考图像灰度图double矩阵 % im2 待配准图像灰度图double矩阵尺寸需与im1一致 % opts 可选参数详见代码 % 输出 % theta 旋转角度单位度im2相对im1的旋转角 % scale 缩放因子im2相对im1的缩放1表示放大 % delta [dy, dx]im2相对im1的平移量 % peak 相位相关置信度越接近1越可靠 arguments im1 double im2 double opts.numR (1,1) double 128 opts.numTheta (1,1) double 256 opts.useWindow (1,1) logical true opts.useHPF (1,1) logical true opts.fineMode (1,1) logical true end % 1. 预处理灰度化、去均值、归一化 im1 preprocess(im1); im2 preprocess(im2); % 2. 计算频谱幅度谱 spec1 spectral_amplitude(im1, opts); spec2 spectral_amplitude(im2, opts); % 3. 对数极坐标变换 lp1 log_polar_transform(spec1, opts.numR, opts.numTheta); lp2 log_polar_transform(spec2, opts.numR, opts.numTheta); % 4. 对数极坐标空间相位相关 [shift_lp, peak] phase_correlation(lp1, lp2); dy shift_lp(1); dx shift_lp(2); % 5. 换算旋转角度与缩放因子 theta dy / opts.numTheta * 360; if theta 180 theta theta - 360; end maxR min(floor((size(im1,1)-1)/2), floor((size(im1,2)-1)/2)); logR_step log(maxR) / (opts.numR - 1); scale exp(-dx * logR_step); % 6. 消除180度歧义并估计平移量 if opts.fineMode bestPeak -inf; for k 1:2 candTheta theta 180 * (k - 1); im2c correct_rs(im2, candTheta, scale); [delta, pk] phase_correlation(im1, im2c); if pk bestPeak bestPeak pk; theta candTheta; delta delta; peak pk; end end else im2c correct_rs(im2, theta, scale); [delta, ~] phase_correlation(im1, im2c); end end主流程就是这样。有几个点值得说明一下。第一步的预处理里去均值和归一化是必须的。如果不去均值直流分量会在频谱中心形成一个极大亮点对数极坐标采样时会严重压制周围刻度信息。归一化保证不同亮度的图像处于同一量纲否则频谱幅度差异过大后续高通滤波参数会失准。第六步的fineMode是我强烈建议保留的。前面说过幅度谱的180度对称性导致角度估计存在歧义直接用单一角度校正可能在一个完全错误的方向上完成平移估计。加了双候选之后代码会用theta和theta180两个角度分别校正图像再做一次空间域相位相关取置信度更高的那个。这个思路用一次FFT成本换回角度方向上的确定性非常划算。3.2 对数极坐标变换的采样细节对数极坐标变换是整套代码里最影响精度的部分。我这里的实现是对幅度谱做半径和角度的网格采样然后双三次插值取值。function lp log_polar_transform(amp, numR, numTheta) % 对数极坐标变换 % amp: fftshift后的频谱幅度谱 % numR: 半径采样数 % numTheta: 角度采样数 [M, N] size(amp); cx (N 1) / 2; cy (M 1) / 2; maxR min(floor((M - 1) / 2), floor((N - 1) / 2)); % 半径从1到maxR对数均匀采样 radius exp(linspace(log(1), log(maxR), numR)); % 角度从0到2*pi去掉最后一个重复采样点 theta linspace(0, 2*pi*(1 - 1/numTheta), numTheta); [Rgrid, Tgrid] ndgrid(radius, theta); Xq cx Rgrid .* cos(Tgrid); Yq cy Rgrid .* sin(Tgrid); lp interp2(amp, Xq, Yq, cubic, 0); end半径为什么从1开始而不是0因为log0没有定义。半径0这一个点正好是频谱直流分量的位置本来在去均值之后它就已经被削弱没必要强行采进来。而且半径过小的区域离散化严重一两个像素的变化就会让对数半径抖动很大采进来只会增加噪声。插值方式我默认用cubic实测比linear精度高。尤其是半径较小时线性插值会让频谱的圆环纹理变糊导致相位相关峰不够尖。如果你的MATLAB版本没有cubic选项可以用spline效果相当但速度会慢一些。还有一个细节lp最后做了一个转置。这是因为ndgrid生成的矩阵形状是numR行numTheta列而我习惯把输出整理成“行对应角度、列对应半径”的格式。这样后面的相位相关在行方向的位移直接对应角度变化列方向的位移直接对应log半径变化调试的时候人脑不容易绕晕。3.3 相位相关与坐标环绕处理相位相关的实现不长但有一个坐标处理很容易写错。function [delta, peak] phase_correlation(im1, im2) F1 fft2(im1); F2 fft2(im2); C F1 .* conj(F2); C C ./ (abs(C) eps); r ifft2(C); % 将零频移到中心使平移量以图像中心为原点 r ifftshift(r); [peak, idx] max(r(:)); [dy, dx] ind2sub(size(r), idx); % 处理环绕平移超过半幅的按负数处理 if dy size(r, 1) / 2 dy dy - size(r, 1); end if dx size(r, 2) / 2 dx dx - size(r, 2); end delta [dy, dx]; end这里最容易踩的坑是ifftshift和fftshift的区别。fft2之后零频在矩阵的(1,1)位置逆变换得到的结果也需要做同样的移位才方便直接当坐标读。偶数尺寸下fftshift和ifftshift等价但奇数尺寸下两者差一个像素所以统一用ifftshift更稳。分母加eps是防止除零。这个细节不加的话频谱中某些位置幅度恰好为0时会出现NaN相位相关结果直接废掉。我用的是eps而不是0既能防除零又不会像加一个固定大数那样污染归一化结果。translation部分的环绕处理也很重要。相位相关的峰可能出现在矩阵边界附近例如真实平移量是[-5, -3]但峰的原始位置可能记录为[size-5, size-3]。如果不把超过半幅的坐标折算成负数后面的校正环节会偏移一大截。4. 实验验证与参数调优4.1 构造已知参数的测试图像写算法最怕没有对照真值。我习惯先合成一个已知参数的测试集把图像做一次确定的旋转、缩放、平移然后用算法去反推。这个验证通过了再上真实图像。%% 生成测试图像 rng(42); N 256; [X, Y] meshgrid(1:N, 1:N); im1 0.6 0.3*sin(X/8) 0.25*cos(Y/6) 0.2*sin((XY)/4 1); im1 im1 0.1*randn(N, N); im1 imgaussfilt(im1, 1); im1 mat2gray(im1); %% 设定真实变换参数 theta_true 25; % 逆时针25度 scale_true 1.25; % 放大1.25倍 delta_true [12, -9]; % 平移量 %% 扩边生成待配准图 pad 160; im1p padarray(im1, [pad pad], mean(im1(:)), both); im2p imrotate(im1p, theta_true, bilinear, crop); im2p imresize(im2p, scale_true, bicubic); im2 crop_center(im2p, N, N); im2 circshift(im2, delta_true); %% 配准 [theta_est, scale_est, delta_est, peak] fm_register(im1, im2); fprintf(theta truth/est: %.2f / %.2f\n, theta_true, theta_est); fprintf(scale truth/est: %.4f / %.4f\n, scale_true, scale_est); fprintf(delta truth/est: %d %d / %d %d\n, ... delta_true(1), delta_true(2), delta_est(1), delta_est(2)); fprintf(peak: %.3f\n, peak);这里有一个重要技巧生成待配准图时先padarray扩边再做旋转和缩放最后中心裁剪回原始尺寸。如果不扩边imrotate的crop选项会直接裁掉图像旋转后露出的四角等于是人为制造了信息丢失会让算法性能看起来比实际更差。真实项目中如果已知图像内容大概率超出画面范围也建议先扩边再估计。4.2 典型输出与180度歧义处理在我机器上跑一组典型输出大概是这样theta truth/est: 25.00 / 24.84 scale truth/est: 1.2500 / 1.2493 delta truth/est: 12 -9 / 12 -9 peak: 0.912角度误差在1度以内缩放误差千分之一量级平移量整数像素完全命中。这个精度对预配准来说是完全够用的。不同版本MATLAB和不同插值参数下数字会有小幅波动但整体数量级应该一致。测试时建议专门试一下角度为0度到180度之间的多组值特别注意接近0度和接近180度的边界情况。如果角度为0缩放估计仍然正确如果角度为180度最后双候选逻辑会稳定选择峰值更大的180度方向。我曾经在只做单候选角度估计时遇到过把实际旋转170度的图像配成-10度的情况原因就是180度歧义没有消除所以fineMode一定要开。4.3 参数调整速查表代码里开放了几个参数实际使用时可以按需调整参数默认值作用调整建议numTheta256角度采样数决定角度分辨率需要更高角度精度时加到512耗时线性增加numR128半径采样数决定缩放分辨率缩放范围大时加到192或256useWindowtrue是否加汉宁窗抑制频谱泄漏图像边缘本来就是平滑过渡时关闭useHPFtrue是否启用高通斜坡滤波低纹理图像建议保持开启fineModetrue是否做180度双候选消歧强烈建议保持开启numTheta和numR不是越大越好。采样数翻倍对数极坐标变换的时间也翻倍但精度提升在numTheta达到512之后会大幅放缓。对多数场景256个角度采样已经把分辨率压到1.4度以内配合亚像素抛物线拟合足够。5. 一些容易踩的坑这套代码我前前后后跑了挺多真实图像有几个坑属于那种你不踩一次根本不会注意的地方。第一个坑是角度符号的方向问题。图像坐标系的y轴向下MATLAB的imrotate正角度对应视觉上的逆时针但不同数据形态下“偏移方向”和“频谱旋转方向”的相对关系可能让你估计出的角度恒为负数。我自己的习惯是先跑一遍已知角度的合成测试如果发现角度符号反了在theta输出处加一个负号就行。如果你用OpenCV或其他库生成测试数据更要在实验里先验证方向不要直接照搬。第二个坑是180度歧义。这个在很多中文资料里都没有明说但原理上非常关键。幅度谱满足共轭对称性意味着转180度之后幅度谱不变。只靠幅度谱做角度估计天然无法区分theta和theta180。解决思路就是我代码里的双候选候选角度分别校正后做空间域相位相关取置信度更高的方向。这个方案简单有效但要求你最终一定做空间域平移精配准不能只输出LP域的估计结果就结束。第三个坑是高通滤波的强度。自然图像的频谱能量高度集中在低频如果不对幅度谱做任何处理对数极坐标里几乎只有一两个低频采样点有能量相位相关峰会钝得没法看。我用的斜坡滤波是在0.15到1.0之间随频率线性增强目的就是给高频频谱一个相对公平的参与机会。但这个强度不是越大越好如果图像本身高频信息较弱强行增强会让噪声主导频谱peak值会掉得很快。第四个坑是图像尺寸。Fourier-Mellin的输入图像最好保持正方形或者长宽比不要太极端而且尺寸最好是偶数。奇数尺寸下频谱中心和图像中心的对应关系会差半个像素虽然对粗配准影响不大但在做亚像素精度需求的项目里会莫名其妙地出现系统性偏差。第五个坑是缩放范围。我的代码里半径采样上限是min(M,N)/2如果待配准图相对参考图缩得太小或放得太大频谱的主要能量会落到采样范围之外估计结果就会失效。这种场景的常规解法是先做一层金字塔粗略估计放大或者缩小回到合理范围再用傅里叶梅林精估计。这个方法我后面有空会单独写一篇专门讲。最后想说的是这套方法不怕光照变化不怕低纹理但怕两幅图不是同一尺度内容。比如一张图是汽车全景另一张图是车身局部这种不是全局缩放关系的情况Fourier-Mellin本身就不适用。使用之前先确认你的问题确实是刚性变换这是省时间的最大前提。本文还有配套的精品资源点击获取