简介：这份MATLAB参考代码围绕傅里叶梅林变换展开，面向图像处理、计算机视觉方向的学习者与研究者，用于解决图像旋转不变性分析与配准问题。资源包共6个文件，以5个bmp测试图像和1个m脚本为主，压缩包约152KB，其中脚本承载核心算法实现，bmp图像则提供lena系列裁剪、旋转、平移等对照样本，便于直接运行验证。内容涉及二维傅里叶变换、梅林参数估计、频谱构建与逆变换等关键环节，并配有图像配准示例，可帮助读者理解旋转与平移场景下的变换流程。需注意描述中提到缩放部分尚未实现，读者可据此自行补充尺度处理逻辑。目前已有3984人学习下载，适合作为入门傅里叶梅林变换、搭建图像配准实验的参考起点。

1. 傅里叶梅林变换参考 MATLAB 代码：从旋转缩放配准到可复现的工程落地

做图像配准的同行大概率都遇到过这种场景：两幅图之间存在旋转、缩放甚至平移，用普通的相位相关法直接算，峰值散得一塌糊涂，根本锁不住位移量。这时候傅里叶梅林变换（Fourier-Mellin Transform，FMT）就是绕不开的一条路——它把旋转和缩放转换成对数极坐标下的平移问题，再用相位相关一次性解出来。这份参考 MATLAB 代码就是围绕这个思路写的，包含完整的频域变换、对数极坐标重采样、相位相关求解和参数反解流程。适合正在做遥感图像配准、医学影像对齐、工业视觉定位的从业者，也适合想搞懂 FMT 底层实现而不是只会调imregcorr的人。代码结构清晰，变量命名直白，拿来改比从头推公式快得多。

2. 傅里叶梅林变换的数学链路：为什么先取模再取对数

2.1 从傅里叶变换到对数极坐标的推导逻辑

FMT 的核心思路并不复杂，但每一步都有明确的数学动机。假设两幅图像 $f_1(x,y)$ 和 $f_2(x,y)$ 之间只存在旋转 $\theta_0$、缩放 $s$ 和平移 $(x_0, y_0)$，那么空域关系可以写成：

$$f_2(x,y) = f_1\left(s(x\cos\theta_0 + y\sin\theta_0) - x_0,\ s(-x\sin\theta_0 + y\cos\theta_0) - y_0\right)$$

直接对这个关系做相位相关是解不出来的，因为旋转和缩放不是平移。FMT 的做法是分三步走：

第一步，对两幅图分别做二维傅里叶变换，得到频谱 $F_1(u,v)$ 和 $F_2(u,v)$。根据傅里叶变换的旋转和缩放性质，空域的旋转对应频域的旋转，空域的缩放对应频域的缩放加幅度因子。

第二步，取频谱的幅度谱 $|F_1|$ 和 $|F_2|$。取模这一步很关键——它丢掉了相位信息，但也丢掉了平移分量。也就是说，取模之后，平移不再影响幅度谱，只剩下旋转和缩放待解决。

第三步，把幅度谱从笛卡尔坐标映射到对数极坐标 $(\log \rho, \theta)$。在对数极坐标下，旋转表现为 $\theta$ 方向的平移，缩放表现为 $\log \rho$ 方向的平移。到这里，旋转和缩放就都变成了平移问题，可以直接用相位相关求解。

这个推导链路是 FMT 的骨架，代码里每一步都有对应的函数实现。我见过不少人直接跳到调用相位相关，结果旋转角度解出来差个几度，回头查半天才发现是对数极坐标采样的范围没设对。

2.2 相位相关求解平移量的实现细节

在对数极坐标域里拿到两个幅度谱之后，下一步就是经典的相位相关。原理很直接：两幅图如果有平移 $(dx, dy)$，它们的互功率谱相位里就编码了这个平移量。

% 相位相关求解平移量 function [shift, peak_value] = phase_correlation(A, B) % A, B 为同尺寸的频域幅度谱（已转到对数极坐标） FA = fft2(A); FB = fft2(B); % 互功率谱，分母加小量防止除零 cross_power = FA .* conj(FB); cross_power = cross_power ./ (abs(cross_power) + eps); % 逆变换得到相关峰 corr_surface = real(ifft2(cross_power)); % 找峰值位置 [peak_value, idx] = max(corr_surface(:)); [row_shift, col_shift] = ind2sub(size(corr_surface), idx); % 处理循环移位，把超过一半的位移映射到负方向 [M, N] = size(corr_surface); if row_shift > M/2 row_shift = row_shift - M; end if col_shift > N/2 col_shift = col_shift - N; end shift = [row_shift, col_shift]; end

这段代码有几个参数需要留意。eps是防止互功率谱归一化时出现除零，虽然理论上不会为零，但实际图像频谱在某些频点可能非常小，不加这个保护会出 NaN。峰值位置到实际位移的映射要注意循环移位的问题——FFT 的移位是循环的，所以超过一半尺寸的位移要减去全长，否则解出来的角度会差 180 度。peak_value本身也是一个质量指标，峰值越尖锐、越接近 1，说明配准越可信；如果峰值扁平或者低于 0.1，基本可以判定这次配准不可靠。

2.3 对数极坐标重采样的参数选择

对数极坐标变换是 FMT 里最容易翻车的一步。代码里用的是logpolar或者自己写的插值函数，核心参数有三个：径向采样点数、角度采样点数、以及径向范围。

% 对数极坐标变换参数设置 num_angles = 360; % 角度方向采样数，决定旋转分辨率 num_radii = 256; % 径向采样数，决定缩放分辨率 radius_min = 1; % 最小半径，避免 log(0) radius_max = min(size(img))/2; % 最大半径，通常取图像短边一半 % 构建对数极坐标网格 theta = linspace(0, 2*pi, num_angles); log_rho = linspace(log(radius_min), log(radius_max), num_radii); [Theta, LogRho] = meshgrid(theta, log_rho); Rho = exp(LogRho); % 映射回笛卡尔坐标 X = Rho .* cos(Theta) + center_x; Y = Rho .* sin(Theta) + center_y; % 双线性插值采样 lp_image = interp2(double(img), X, Y, 'linear', 0);

num_angles决定了旋转角度的分辨率。360 个采样点对应 1 度分辨率，如果应用场景需要亚度级配准，得加到 720 甚至 1440。num_radii决定缩放分辨率，256 个点通常够用，但如果缩放比例跨度大（比如 0.5 到 2 倍），可以适当增加。radius_min不能取 0，因为 $\log(0)$ 是负无穷，代码里取 1 是常见做法。radius_max取图像短边的一半，超过这个范围采样点会落到图像外面，插值出来的全是零，反而干扰相位相关。

注意：对数极坐标变换前最好加一个汉宁窗或者汉明窗，抑制频谱泄漏。不加窗的话，图像边缘的突变会在频谱里产生十字条纹，直接影响相位相关的峰值质量。

3. 参考代码的模块拆解与运行流程

3.1 主脚本的调用顺序与数据流转

这份参考代码的主脚本结构很清晰，基本按照「读图 → 预处理 → FMT → 相位相关 → 参数反解 → 验证」的顺序组织。我把它拆成几个关键模块来看。

%% 主流程 % 1. 读取图像并转为灰度 img1 = imread('reference.png'); img2 = imread('target.png'); if size(img1, 3) == 3 img1 = rgb2gray(img1); end if size(img2, 3) == 3 img2 = rgb2gray(img2); end % 2. 统一尺寸，取短边对齐 min_size = min([size(img1); size(img2)]); img1 = imresize(img1, [min_size(1), min_size(2)]); img2 = imresize(img2, [min_size(1), min_size(2)]); % 3. 加窗抑制频谱泄漏 win = hanning(min_size(1)) * hanning(min_size(2))'; img1 = double(img1) .* win; img2 = double(img2) .* win; % 4. 傅里叶变换取幅度谱 F1 = abs(fftshift(fft2(img1))); F2 = abs(fftshift(fft2(img2))); % 5. 对数极坐标变换 lp1 = logpolar_transform(F1, num_angles, num_radii); lp2 = logpolar_transform(F2, num_angles, num_radii); % 6. 相位相关求旋转和缩放 [shift_rs, peak_rs] = phase_correlation(lp1, lp2); % 7. 反解旋转角和缩放因子 angle_shift = shift_rs(2) * (360 / num_angles); % 角度方向 scale_shift = shift_rs(1) * (log(radius_max/radius_min) / num_radii); scale_factor = exp(scale_shift); % 8. 对 img2 做反向旋转缩放，再求平移 img2_corrected = imrotate(img2, -angle_shift, 'bilinear', 'crop'); img2_corrected = imresize(img2_corrected, 1/scale_factor); [shift_xy, peak_xy] = phase_correlation(img1, img2_corrected);

第 2 步统一尺寸容易被忽略。如果两幅图尺寸不一致，FFT 之后频率轴对不上，对数极坐标采样也会错位。常见做法是取短边对齐，或者都缩放到同一尺寸。第 3 步加窗是血泪经验——不加窗的话，图像边界的突变会在频谱里产生明显的十字亮线，对数极坐标变换后这些亮线变成干扰条纹，相位相关峰值直接掉到 0.05 以下。第 7 步的反解公式里，角度方向的分辨率是 $360/\text{num_angles}$，径向的分辨率是 $\log(\text{radius_max}/\text{radius_min})/\text{num_radii}$，这两个换算关系必须和前面的采样参数严格对应，否则解出来的角度和缩放全是错的。

3.2 参数反解与图像校正的验证方法

解出旋转角、缩放因子和平移量之后，怎么验证结果对不对？代码里给了一个简单的验证方式：把 target 图按解出的参数做反向变换，然后和 reference 图做差值或者叠加显示。

% 验证配准结果 % 构建仿射变换矩阵 tform = affine2d([scale_factor*cosd(angle_shift), -scale_factor*sind(angle_shift), 0; scale_factor*sind(angle_shift), scale_factor*cosd(angle_shift), 0; 0, 0, 1]); % 注意：这里只是示意，实际平移量需要单独叠加 % 或者直接用 imwarp output_view = imwarp(img2, tform, 'OutputView', imref2d(size(img1))); % 叠加显示 figure; imshowpair(img1, output_view, 'blend'); title('配准叠加结果');

验证的时候重点看两个指标。一个是相位相关的峰值，旋转缩放阶段的峰值和最终平移阶段的峰值都应该在 0.3 以上，低于这个值基本说明配准失败。另一个是叠加显示的边缘对齐程度，如果边缘有重影或者错位，说明参数反解有问题。我一般还会算一下配准后的互信息或者归一化互相关，作为定量指标。

提示：如果峰值一直上不去，先检查对数极坐标变换的范围是不是覆盖了图像的有效频谱区域。频谱中心是低频，往外是高频，如果radius_max取太小，高频信息全丢了，配准精度自然差。

4. 避坑与排查：FMT 配准中最容易翻车的五个点

4.1 峰值弥散导致角度解算偏差

现象：相位相关峰值低于 0.1，解出的旋转角和真实值差 5 度以上，缩放因子也偏得离谱。

原因：最常见的原因是对数极坐标变换时没有加窗，或者加的窗类型不对。图像边界的突变在频谱里产生十字条纹，这些条纹在对数极坐标下变成贯穿整个角度方向的亮线，严重干扰相位相关。另一个原因是两幅图的频谱幅度差异太大，比如一幅图整体偏暗，另一幅偏亮，幅度谱的动态范围不一致。

解决：加汉宁窗或汉明窗，窗函数要同时作用于行和列。如果亮度差异大，先做直方图均衡或者归一化。还可以对幅度谱取对数 $\log(1 + |F|)$，压缩动态范围，让高低频信息更均衡。

4.2 角度解算的 180 度歧义

现象：解出的旋转角总是差 180 度，或者在某些图像对上解出两个可能的峰值。

原因：傅里叶幅度谱具有中心对称性，$|F(u,v)| = |F(-u,-v)|$，这意味着旋转 $\theta$ 和 $\theta + 180°$ 在幅度谱上是一样的。相位相关只能解出 $[0, 180°)$ 范围内的角度，超过这个范围就有歧义。

解决：代码里通常只解 $[0, 180°)$ 的角度，然后分别尝试 $\theta$ 和 $\theta + 180°$ 两个候选，用最终的空域配准峰值来判断哪个是对的。这一步在参考代码里是通过两次相位相关实现的——先解旋转缩放，再解平移，平移阶段的峰值能区分 180 度歧义。

4.3 缩放因子反解时的对数域溢出

现象：解出的缩放因子是 NaN 或者无穷大，程序直接报错。

原因：radius_min设成了 0，导致 $\log(0)$ 溢出。或者radius_max设得比图像尺寸还大，采样点落到图像外面，插值结果全是零，相位相关解出的位移没有意义。

解决：radius_min至少取 1，radius_max取图像短边的一半。如果缩放范围已知，比如只在 0.8 到 1.2 倍之间，可以把radius_max设成radius_min * 1.5，缩小搜索范围，提高精度。

4.4 大角度旋转时的频谱混叠

现象：旋转角超过 45 度时，配准精度明显下降，峰值也变低。

原因：对数极坐标变换在角度方向是均匀采样的，但笛卡尔坐标下的旋转在角度方向并不是均匀的。大角度旋转时，角度方向的采样密度相对不足，产生混叠。

解决：增加num_angles，比如从 360 加到 720 或 1440。同时可以在对数极坐标变换前对频谱做一次低通滤波，抑制高频混叠。

4.5 平移量解算时的边界效应

现象：平移量解出来总是偏大或偏小，尤其是当真实平移量接近图像尺寸的一半时。

原因：FFT 的循环移位特性导致边界处的位移出现环绕。如果真实平移量超过图像尺寸的一半，相位相关解出的位移会映射到相反方向。

解决：代码里已经做了循环移位的处理，把超过一半的位移映射到负方向。但如果真实平移量确实很大，建议先做粗配准，把平移量降到图像尺寸的四分之一以内，再用 FMT 做精配准。

5. 进阶技巧：把 FMT 配准精度再压一压

5.1 亚像素精度插值与峰值拟合

相位相关解出的位移是整数像素级的，如果应用场景需要亚像素精度，得在峰值附近做插值拟合。常见做法是取峰值及其周围 3×3 邻域，用抛物线拟合或者高斯拟合找亚像素峰值位置。

% 亚像素峰值拟合 function [sub_shift, sub_peak] = subpixel_peak(corr_surface, peak_idx) [M, N] = size(corr_surface); [row, col] = ind2sub([M, N], peak_idx); % 取 3x3 邻域，处理边界 rows = mod((row-1:row+1)-1, M) + 1; cols = mod((col-1:col+1)-1, N) + 1; patch = corr_surface(rows, cols); % 行方向抛物线拟合 denom_r = 2 * (2*patch(2,2) - patch(1,2) - patch(3,2)); if abs(denom_r) > eps delta_r = (patch(3,2) - patch(1,2)) / denom_r; else delta_r = 0; end % 列方向抛物线拟合 denom_c = 2 * (2*patch(2,2) - patch(2,1) - patch(2,3)); if abs(denom_c) > eps delta_c = (patch(2,3) - patch(2,1)) / denom_c; else delta_c = 0; end sub_shift = [row + delta_r, col + delta_c]; sub_peak = patch(2,2) - 0.25 * (patch(3,2)-patch(1,2))*delta_r ... - 0.25 * (patch(2,3)-patch(2,1))*delta_c; end

这段拟合能把位移精度从 1 像素压到 0.1 像素左右。delta_r和delta_c是行和列方向的亚像素偏移量，sub_peak是拟合后的峰值。注意边界处理，峰值如果在图像边缘，邻域会环绕，用mod处理循环索引。

5.2 多尺度策略与鲁棒性提升

如果图像之间不仅存在旋转缩放平移，还有局部形变或者噪声，单次 FMT 可能不够稳。我一般会加一个多尺度策略：先降采样做粗配准，再逐级升采样做精配准。

粗配准阶段用低分辨率快速算一遍，把旋转角锁定到 5 度以内，缩放锁定到 10% 以内。中配准阶段在这个范围内细化，精配准阶段做最终优化。这样比直接上高分辨率快得多，而且不容易陷入局部极值。

5.3 从参考代码到工程落地的几个习惯

这份参考代码拿来跑通不难，但要用到实际项目里，有几个习惯我建议从一开始就养成。第一，每次配准都记录相位相关峰值，峰值低于阈值直接标记为失败，不要硬着头皮往下走。第二，对数极坐标变换的参数不要写死，根据图像尺寸和预期缩放范围动态计算。第三，配准结果一定要做可视化验证，imshowpair的blend模式或者montage模式都行，肉眼过一遍比什么都靠谱。

从那以后我每次做 FMT 配准，都强制走一遍「加窗 → 对数极坐标 → 相位相关 → 峰值检查 → 可视化验证」的流程，少一步都不敢往下走。希望帮到你。

本文还有配套的精品资源，点击获取