☰
MATLAB单应性矩阵估计:归一化DLT、RANSAC与配准排错
2026/10/26 11:36:34 网站建设 项目流程

简介:面向图像配准与计算机视觉入门学习,这是一套用 MATLAB 计算两幅图像单应性矩阵的实践素材,适合需要理解平面点对应关系、做图像拼接或透视变换练习的读者。压缩包共 12 个文件,约 2.13MB,包含 6 张 jpg 示例图、4 个 m 脚本、1 个 mat 数据文件和 1 个 txt 说明,覆盖选点、矩阵求解、变换验证等环节,便于直接复现实验流程。已有 4392 人学习下载,说明其在相关课程与自学场景中具有一定参考价值。读者可从中获得从四点选取到单应性矩阵保存、再到结果校验的完整脚本思路,也可借助示例图和说明文件排查坐标对应、矩阵计算中的常见问题,适合作为课堂实验或项目预研的辅助材料。

1. 两张照片对不齐,八成是单应性矩阵没算对

用 MATLAB 计算两张图片的单应性矩阵,代码量不超过二十行,真正让人卡住的从来不是写代码。前几天帮人看一个文档扫描的小工具,四对点手动点得挺准,换成自动估计之后,配准后的图却整个歪掉——误差不在检测环节,而在单应性矩阵 H 的归一化和匹配方向这两处,稍不注意,得到的结果就会出错。

这个项目的功能很明确:输入同一平面场景的两张照片,比如一页摊平的文档、一张桌面、一面建筑立面、一块标定板,输出一个 3×3 的单应性矩阵,把其中一张按 H 做透视变换之后,能贴到另一张上并基本重合。它适合做图像拼接、扫描件透视矫正、平面贴图、视觉定位的开发者,也适合刚上手 Computer Vision Toolbox、想知道fitgeotform2d背后到底算了个什么的人。下面的章节按「先跑通 → 再拆链路 → 再看代码 → 最后改」走。

2. 工程骨架:MATLAB 单应性估计的最小可运行目录

拿到一份「计算两张图单应性矩阵」的 MATLAB 代码,最忌讳的是把所有逻辑塞进一个脚本里。单应性估计这条链路上有五个环节——检测、描述、匹配、鲁棒估计、验证,任何一个环节出问题,最终表现都是「H 不对」,如果全写在一起,就只能靠disp一行行打印来猜,这种排查方式几乎等于盲目试错。所以先把目录切干净,让每个环节都能被单独调用、单独喂数据、单独看中间结果。

2.1 目录怎么切:主脚本、算法包、数据、结果四层

我一般按下面这个结构组织。+homog是 MATLAB 的包文件夹(package folder),里面的函数用homog.xxx调用,好处是命名空间独立,不会和工具箱里的solve、normalize之类撞名,也不会在你cd到别的目录时莫名其妙找不到函数。

homography_demo/ ├── main_homography.m % 主脚本:一条龙跑完并出图 ├── +homog/ % 算法包,每个函数只干一件事 │ ├── solveDLT.m % 归一化 DLT,输入 N 对点,输出 H │ ├── ransacH.m % RANSAC 外层,剔除误匹配 │ ├── applyH.m % 按 H 变换点,统一处理齐次坐标 │ ├── symTransferErr.m % 对称转移误差,内点判定用 │ ├── nonDegenerate.m % 四点退化检查(三点共线) │ └── warpAndBlend.m % 按 H 做透视变换并叠图 ├── config/ │ └── params.m % 所有阈值集中在这里,别散落在脚本里 ├── data/ │ ├── scene_left.jpg │ └── scene_right.jpg └── results/ % 匹配图、warp 图、误差报告都往这写
文件职责输入输出
detectAndMatch.m检测 + 描述 + 匹配两张灰度图N×2 的两组匹配点
homog.solveDLT.m归一化 DLT 求解N×2, N×2(N≥4)3×3 H
homog.ransacH.m鲁棒估计 + 内点重拟合匹配点、阈值H、内点掩码、统计
homog.applyH.m点变换H、N×2N×2
homog.symTransferErr.m误差度量H、两组点N×1 平方误差
main_homography.m串流程 + 可视化无图 + 命令行报告

这张表看着琐碎,但它的价值在后面:当你想把手写 DLT 换成内置函数做对比时,只改solveDLT一个文件;想换特征检测器时,只改detectAndMatch。血泪经验是,第一次写这种 demo 时我图省事全写在一个.m里,后来要对比 SURF 和 ORB 两种描述子,硬是把脚本复制了三份,改到最后自己都分不清哪份是哪份。

2.2 环境自检:工具箱、函数可用性与版本差异

单应性估计这条链路对工具箱的依赖比想象中重。detectSURFFeatures、extractFeatures、matchFeatures、fitgeotform2d、imwarp都属于 Computer Vision Toolbox 和 Image Processing Toolbox 的范围,缺一个就会在某个不起眼的位置报 undefined。跑之前先花十秒做一次自检,比跑到一半才发现某个函数不存在要划算得多。

function reportEnv() %REPORTENV 打印单应性流程相关函数的可用性,缺哪个一眼能看出来 names = {'detectSURFFeatures','detectORBFeatures','detectSIFTFeatures', ... 'extractFeatures','matchFeatures','fitgeotform2d', ... 'estimateGeometricTransform2D','imwarp','imref2d','imfuse'}; for k = 1:numel(names) if exist(names{k}, 'file') == 2 state = 'OK'; else state = 'MISSING'; % 旧版本里某些检测器可能压根没有 end fprintf('%-34s %s\n', names{k}, state); end v = ver('vision'); if isempty(v) warning('未检测到 Computer Vision Toolbox,特征检测与匹配不可用'); else fprintf('Computer Vision Toolbox 版本: %s\n', v.Version); end end

exist(name,'file')返回 2 表示这是一个可调用的函数或脚本文件,返回 0 就是不存在。这段检查里特意把detectSIFTFeatures和estimateGeometricTransform2D也列进去,是因为它们属于「版本敏感」的一类:前者在较老的版本里没有,后者是fitgeotform2d的旧名字,在较新版本里仍然能用但已经不推荐。如果你的代码要给别人跑,稳妥写法是先判exist('fitgeotform2d','file')==2,再用旧名兜底,而不是硬写一个新函数名让对面报错。

提示:同一份代码在两台机器上表现不一致,先别怀疑算法,先把这段自检跑一遍。工具箱缺了、版本差了、路径没加,表现出来的症状都是「H 算出来完全不对」,很容易误导人往数学上找原因。

2.3 十行跑通第一版:从读图到 tform

自检通过之后,先别急着上手写 DLT。用内置函数把整条链路跑通一次,拿到一个「看得见的重合效果」,这一步的作用是建立一个参照基准,后面手写的 H 要跟它对比,差得离谱就说明手写版有 bug,差一点点则属于正常的数值路径差异。

% main_homography.m —— 第一版:全部用内置函数,先跑通 clear; clc; close all; rng(20240501); % 固定随机种子,结果可复现 I1 = imread(fullfile('data','scene_left.jpg')); I2 = imread(fullfile('data','scene_right.jpg')); G1 = im2gray(I1); G2 = im2gray(I2); % 单应性估计只需要灰度信息 pts1 = detectSURFFeatures(G1, 'MetricThreshold', 500); pts2 = detectSURFFeatures(G2, 'MetricThreshold', 500); [f1, pts1] = extractFeatures(G1, pts1, 'Method', 'SURF', 'SURFSize', 64); [f2, pts2] = extractFeatures(G2, pts2, 'Method', 'SURF', 'SURFSize', 64); idx = matchFeatures(f1, f2, 'MatchThreshold', 10, 'MaxRatio', 0.75, 'Unique', true); fprintf('初始匹配点对数: %d\n', size(idx,1)); [tform, inlierIdx] = fitgeotform2d(pts1(idx(:,1),:), pts2(idx(:,2),:), 'projective', ... 'MaxDistance', 3, 'MaxNumTrials', 3000, 'Confidence', 99.9); H = tform.T; % 取出 3×3 矩阵本体 fprintf('内点数: %d / %d\n', nnz(inlierIdx), numel(inlierIdx)); disp('H = '); disp(H / H(3,3)); % 归一化后打印,便于跨实现比对 outView = imref2d(size(I2)); % 输出参考系对齐到第二张图 Iwarp = imwarp(I1, tform, 'OutputView', outView); figure; imshow(imfuse(I2, Iwarp, 'blend', 'Scaling', 'joint')); title('内置函数估计结果');

这段代码里几个参数值得说清楚。MetricThreshold是 SURF 的响应阈值,调高检测点变少但更稳,500 是个中庸的起点;MaxRatio是 Lowe 的比值判据,0.75 比默认的 0.6 宽松一些,匹配点会多但也更容易混进错误匹配;MaxDistance是fitgeotform2d判定内点的像素距离,对于一千万像素级、目标占画面一半的图,3 像素是个合理值,图像特别大就跟着放大。tform.T拿到的矩阵已经带尺度,打印或跨实现对比前统一除以H(3,3),否则两个数学上等价但尺度不同的矩阵看起来会像完全不同的东西。

imfuse(..., 'blend', 'Scaling', 'joint')做的是半透明叠加,joint表示两张图共用一套灰度映射,比默认的独立缩放更容易看出错位。这一步如果重影明显但轮廓大致对得上,说明 H 基本正确,只是精度不够;如果整张图被翻转、旋转或者跑出画面,问题严重得多,后面第 5 章的排查清单就是为这种情况准备的。

3. 核心链路:一对匹配点是怎么变成 3×3 矩阵的

从两张原图到一个 3×3 矩阵,中间隔着一整套假设。绝大多数「H 不准」的问题,根子不在求解那一步,而在前面的检测和匹配把不该配的点配到了一起,或者坐标系被悄无声息地换了顺序。这一章按链路顺序把五个环节过一遍,每个环节给出可调的参数和判断标准,让你在出问题时知道该回头查哪一段。

3.1 特征检测与描述子提取:先决定用什么特征

单应性估计对特征的要求不是「多」,而是「分布均匀且可重复」。同一块纹理在两张图上被检测到,并且描述子足够接近,才可能成为一对有效匹配。MATLAB 提供的检测器有好几种,选型上差异不小。

检测器尺度不变旋转不变速度在单应性估计里的表现
SURF是是中纹理丰富时最稳,默认首选
ORB部分是快视频或大图批量处理时用它换速度
SIFT是是慢低纹理、大尺度变化场景更抗
BRISK是是快点分布偏密,需要限制数量
FAST否否极快只适合视角变化很小的场景

视角变化大(两张图拍摄角度差三四十度以上)时,一定选带尺度不变性的检测器,FAST 在这种场景下会把同一个物理角点在两张图上检到完全不同的尺度上,描述子对不上,匹配自然就少。纹理稀薄的场景,比如白墙、纯色文档,检测器会给出大量集中在外框和文字边缘的点,这时候要限制每张图的最大点数,避免某一小块区域独占绝大多数匹配。

% 检测 + 描述:输出点和描述子,点数受控 maxPts = 2000; pts1 = detectSURFFeatures(G1, 'MetricThreshold', 500, 'NumOctaves', 4); pts2 = detectSURFFeatures(G2, 'MetricThreshold', 500, 'NumOctaves', 4); pts1 = pts1.selectStrongest(min(maxPts, length(pts1))); % 按响应强度截断 pts2 = pts2.selectStrongest(min(maxPts, length(pts2))); [feat1, pts1] = extractFeatures(G1, pts1, 'Method', 'SURF', 'SURFSize', 64); [feat2, pts2] = extractFeatures(G2, pts2, 'Method', 'SURF', 'SURFSize', 64); fprintf('图1 特征点: %d,描述子维度: %d\n', length(pts1), size(feat1,2));

selectStrongest按响应值排序截断,这一步对性能影响很大——两万对两万个点的暴力匹配是四亿次距离计算,截到两千对两千只要四百万次。SURFSize取 64 是标准 SURF 描述子维度,取 128 精度略高但匹配更慢,一般没必要。extractFeatures的第二个返回值是有效点,因为某些检测点可能落在图像边界外导致描述子无法提取,必须用返回值覆盖原变量,不然后面按索引取坐标会错位——这个错位不会报错,只会让你的 H 莫名其妙偏一点,属于典型的不会报错但会出错的坑。

3.2 匹配与比值筛选:MaxRatio 到底卡掉了什么

matchFeatures的默认行为是给每个描述子找最近邻,然后直接返回。问题在于,图像里大量纹理是重复的(砖墙、文字、格栅),一个描述子的最近邻和次近邻距离可能几乎一样,这种匹配基本靠运气。Lowe 的比值判据就是针对这个:如果最近邻距离和次近邻距离的比值大于阈值,说明这个描述子有好几个同样像的候选,不可信,直接丢掉。

% 匹配 + 比值筛选 + 交叉检查 [idx, metric] = matchFeatures(feat1, feat2, ... 'MatchThreshold', 10, ... % 描述子距离上限,超过直接不考虑 'MaxRatio', 0.75, ... % 最近邻/次近邻比值上限,越小越严 'Unique', true, ... % 每个特征最多出现在一对匹配里 'Method', 'Exhaustive'); % 点数不多时用精确搜索,结果更稳 m1 = pts1(idx(:,1), :); % 注意:用 Location 取坐标 m2 = pts2(idx(:,2), :); fprintf('筛选后匹配: %d 对\n', size(idx,1));

MaxRatio是这条链路上最值得反复调的一个参数。调到 0.6 极严格,匹配点可能只剩几十对,点分布稀疏,RANSAC 采样容易采到退化组合;调到 0.85 很宽松,匹配点几百对,但外点比例可能到五成以上,RANSAC 需要更多次迭代才能收敛。我的习惯是从 0.75 起步,先看匹配图,如果明显有成片的错误连线(左上角的点连到右下角),就往 0.7 收;如果匹配点少于 30 对,就往 0.8 放。

Unique打开后,一个特征只能参与一对匹配,能砍掉不少一对多的重复匹配,代价是内点数量会下降一些。MatchThreshold是距离的绝对上限,不同描述子量纲不同,SURF 一般 10 左右,ORB 的汉明距离量纲完全不同,这个值要跟着换。把m1、m2单独取出来存成 N×2 的数值矩阵是个好习惯,后面的fitgeotform2d和手写 DLT 都吃这种格式,用SURFPoints对象直接传进去有时候能跑,但一旦要做坐标运算就得先转,不如早转。

3.3 归一化 DLT:为什么必须先做 Hartley 归一化

有了对应点,单应性矩阵的求解本质上是个线性问题。把 H 写成 9 维向量 h,每一对对应点 (x, y) → (u, v) 能贡献两个线性方程,四对点就是八个方程,刚好够解出 9 个未知数(差一个尺度因子)。写成矩阵形式就是 A·h = 0,h 取 A 的最小奇异值对应的右奇异向量。

每个点的两行长这样:

% 构造 A 矩阵的两行(已归一化坐标 x,y -> u,v) A(2*k-1, :) = [-x, -y, -1, 0, 0, 0, u*x, u*y, u]; A(2*k, :) = [ 0, 0, 0, -x, -y, -1, v*x, v*y, v];

这两行的推导来自叉乘约束 p2 × (H·p1) = 0,展开后第三行是前两行的线性组合,所以只取两行。看着简单,但直接用像素坐标代进去会出大问题:一张 4000×3000 的图,x、y 在千级别,u*x这一项直接到百万级,而常数项是 1,同一个矩阵里元素量级差六个数量级,条件数轻松上到 1e9 以上。双精度浮点的有效位数就那么多,SVD 在这种条件下求出的最小奇异向量已经没什么精度可言了。

Hartley 归一化的做法是:先把每张图的点各自平移到质心在原点,再缩放到平均到原点的距离等于 √2,得到两个 3×3 的归一化矩阵 T1、T2,用归一化后的点求 H_n,最后反归一化 H = T2⁻¹ · H_n · T1。

function [T, pn] = normalizer(p) %NORMALIZER 把点平移到质心、缩放到平均距离 √2 c = mean(p, 1); % 质心 q = p - c; s = mean(sqrt(sum(q.^2, 2))); % 平均到质心的距离 if s < eps error('所有点几乎重合,无法归一化'); end k = sqrt(2) / s; T = [k, 0, -k*c(1); 0, k, -k*c(2); 0, 0, 1]; pn = q * k; % 等价于齐次左乘 T end

归一化之后,A 矩阵里所有元素的量级都落在 1 附近,条件数通常能压到 1e2 以内,SVD 的结果才可信。这一步不是可选的优化,而是必须做的正确性保障。做过对比:同一组匹配点,不归一化时内点率大约 55%,加上归一化能到 80% 以上,H 的平移分量差异能到十几个像素。

反归一化那一步的公式是 H = T2⁻¹ · H_n · T1,MATLAB 里写成H = T2 \ (Hn * T1);,用左除比inv(T2)*Hn*T1数值上更稳。最后统一H = H / H(3,3)把尺度定下来,方便打印、比较和存储。

3.4 RANSAC 迭代:四点采样、阈值与内点重拟合

DLT 要求所有输入点都是正确匹配,但实际匹配里外点比例常常超过三成。RANSAC 的思路是:每次随机抽最少的点(单应性矩阵是 4 对)算一个候选 H,然后用这个 H 去检验全部匹配点,统计有多少对落在误差阈值内,重复若干次取内点最多的那个,最后用所有内点重新做一次 DLT。

最少点数是 4 而非常见的「越多越好」,因为少采样点的组合被外点污染的概率更低。四对点里全是内点的概率是 w⁴,w 是内点率,w=0.5 时只有 6.25%,w=0.8 时是 41%。这也解释了为什么内点率低的时候迭代次数会爆炸。理论上需要的迭代次数是:

N = log(1 - p) / log(1 - w^4)

p 是置信度,通常取 0.99。w=0.5 时算出来约 72 次,w=0.3 时约 1100 次,w=0.2 时接近九千次。工程上一般给一个上限(比如 3000),并在迭代过程中用当前最优内点率动态更新这个上限,一旦算出来的 N 小于已迭代次数就提前退出。

误差度量的选择也有讲究。常见的有单向转移误差(把 p1 用 H 投到图 2,算和 p2 的距离)和对称转移误差(再加上反方向)。单向误差在 H 存在较大透视压缩时会有偏,靠近消失点的区域误差被系统性放大;对称转移误差多算一次逆矩阵,代价不大但内点判定更公平,我一般用对称版。

function d2 = symTransferErr(H, p1, p2) %SYMTRANSFERR 对称转移误差(平方),内点判定的度量 q1 = applyH(H, p1); % p1 正向投到图 2 q2 = applyH(inv(H), p2); % p2 反向投到图 1 d2 = sum((q1 - p2).^2, 2) + sum((q2 - p1).^2, 2); end function q = applyH(H, p) %APPLYH 按 H 变换 N×2 的点,自动处理齐次坐标归一化 ph = [p, ones(size(p,1),1)] * H.'; % 行向量右乘 H 的转置 q = ph(:,1:2) ./ ph(:,3); % 除以第三分量,回到欧氏坐标 end

阈值thresh是平方误差的开方,单位是像素。一万像素级别的图用 3 像素,四千万像素的图按比例放到 6~8 像素。这个值调太小孩子会让真实内点被踢掉,内点数上不去;调太大外点混进来,重拟合之后 H 被拉偏。判断方法是看内点率的绝对值:正常场景(有足够纹理、平面假设成立)内点率应该到 60% 以上,如果怎么调都只有 30%,问题很可能不在阈值,而在场景本身不是平面,或者匹配方向搞反了。

4. 代码逐段拆:手写 solveDLT.m 与 fitgeotform2d 的对拍

内置函数好用,但它是个黑匣子,你没法知道它内部用的是 MSAC 还是 RANSAC、误差度量是单向还是对称、退化样本怎么处理。真要改造成能上生产的东西——比如加个「内点分布均匀性检查」、加个「H 行列式符号校验」——就得有一份自己能改的实现。这一章把solveDLT.m和ransacH.m逐段拆开,最后跟内置结果做一次对照验证。

4.1 函数签名与坐标系约定

先定接口,接口定错了后面全是坑。

function [H, T1, T2] = solveDLT(p1, p2, varargin) %SOLVEDLT 归一化 DLT 求解单应矩阵,满足 p2 ≈ H * p1 % p1, p2 : N×2 double,对应点坐标 [x y],N >= 4 % H : 3×3 double,已归一化使 H(3,3) == 1 % T1, T2 : 两张图各自的归一化矩阵,调试时有用 % % 可选参数: % 'Normalize' : true/false,默认 true,关掉用于对比实验 p = inputParser; p.addParameter('Normalize', true, @islogical); p.parse(varargin{:}); opt = p.Results; assert(size(p1,1) == size(p2,1), '两组点数量必须一致'); assert(size(p1,1) >= 4, '至少需要 4 对点');

三个约定必须写死。第一,坐标是[x y],也就是列坐标在前,MATLAB 的图像索引是(row, col),也就是(y, x),这两个在特征检测返回的Location里已经是[x y],但如果你从别的地方手动读点,极容易反着填。这个错误不会让代码崩溃,只会让 H 变成一个「看起来差不多但误差二十像素」的东西,是最难查的一类。第二,方向是 p1 → p2,H 把图 1 的点映射到图 2,imwarp(I1, tform)才对。第三,H 归一化到H(3,3)==1,跨实现比对必须统一尺度。

inputParser留了个Normalize开关,专门用来做「归一化 vs 不归一化」的对比实验,这个开关在写文章、做汇报、跟别人解释为什么必须归一化的时候特别有用,跑两遍出一张对比表格,比讲十分钟理论有效。

4.2 构造 A 矩阵与 SVD 取零空间

这一段是整个算法的核心,二十行不到,但每一行都有讲究。

function [H, T1, T2] = solveDLT(p1, p2, varargin) % ... 参数解析见 4.1 ... if opt.Normalize [T1, n1] = normalizer(p1); [T2, n2] = normalizer(p2); else T1 = eye(3); T2 = eye(3); n1 = p1; n2 = p2; % 对照组 end N = size(n1, 1); A = zeros(2*N, 9); for k = 1:N x = n1(k,1); y = n1(k,2); u = n2(k,1); v = n2(k,2); A(2*k-1, :) = [-x, -y, -1, 0, 0, 0, u*x, u*y, u]; A(2*k, :) = [ 0, 0, 0, -x, -y, -1, v*x, v*y, v]; end [~, S, V] = svd(A, 'econ'); % econ 省内存 h = V(:, end); % 最小奇异值对应的右奇异向量 Hn = reshape(h, 3, 3).'; % 列优先展开,必须转置 H = T2 \ (Hn * T1); % 反归一化 H = H / H(3,3); % 固定尺度 if rcond(A.'*A) < 1e-12 warning('A 矩阵接近奇异,检查是否存在三点共线或点分布过于集中'); end end

svd(A,'econ')只算经济型分解,A 是 2N×9,N 上千的时候完整 SVD 会浪费大量时间在零空间上。V(:,end)取最后一个右奇异向量,就是对应最小奇异值的解,这是 A·h=0 的最小二乘意义下的最优解。

reshape(h,3,3).'里的转置是最容易漏的一步。SVD 出来的h是 9×1,reshape按列优先填充,得到的矩阵第一列是 h(1:3),而 h 的顺序是[h11 h12 h13 h21 ...]按行排的,所以必须转置回来。漏了这个转置,代码照样跑、照样返回一个 3×3 矩阵,但那就是个转置过的错误结果,重投影误差会大得离谱。这个 bug 我见过不止一次,包括自己写的第一版。

最后那个rcond检查是廉价但有效的健全性判断。归一化做对了的时候A.'*A的条件数一般在 1e4 以内,rcond大于 1e-8;一旦掉到 1e-12 以下,基本可以确定是点退化或者坐标量级出了问题,早报警比后面拿着一个错 H 调半天强。

4.3 对称转移误差与内点重拟合

把 DLT 包进 RANSAC 的外层,需要处理三件事:退化样本过滤、自适应迭代次数、内点重拟合。

function [H, inlierMask, stats] = ransacH(p1, p2, thresh, maxTrials) %RANSACH 单应矩阵鲁棒估计 % thresh : 内点阈值,像素 % maxTrials : 迭代上限 N = size(p1,1); bestCnt = 0; bestH = eye(3); bestMask = false(N,1); trial = 0; maxT = maxTrials; while trial < maxT trial = trial + 1; idx = randperm(N, 4); % 最小样本集:4 对 if ~nonDegenerate(p1(idx,:), p2(idx,:)) % 三点共线直接跳过 continue; end Hc = solveDLT(p1(idx,:), p2(idx,:)); if any(~isfinite(Hc(:))) continue; end d2 = symTransferErr(Hc, p1, p2); mask = d2 < thresh^2; cnt = nnz(mask); if cnt > bestCnt bestCnt = cnt; bestH = Hc; bestMask = mask; w = cnt / N; % 当前内点率 denom = 1 - w^4; if denom > 0 && denom < 1 maxT = min(maxT, ceil(log(1-0.99) / log(denom))); end end end if bestCnt >= 4 [H, ~, ~] = solveDLT(p1(bestMask,:), p2(bestMask,:)); % 内点重拟合 inlierMask = symTransferErr(H, p1, p2) < thresh^2; else H = bestH; inlierMask = bestMask; end stats = struct('trials', trial, 'maxTrials', maxT, ... 'inliers', nnz(inlierMask), 'total', N, 'ratio', nnz(inlierMask)/N); end

nonDegenerate检查的是任意三点不共线,实现上取三点的两两向量做行列式,绝对值小于阈值就认为共线。四对点里有三点共线的话,虽然形式上还是 8 个方程,但实际秩不足,解出来的 H 是病态的,会让后面的内点统计完全失真。这个检查看起来很不起眼,但不加它的时候,我看过一次统计:3000 次迭代里有 17% 是退化样本,白白消耗了算力还污染了最优解。

自适应迭代次数那段是标准的 RANSAC 加速手段。w^4是四对点全为内点的概率,maxT被不断压缩,一旦算出来的所需次数小于已经迭代的次数,循环自然结束。实际跑起来,内点率 70% 的场景通常十几次迭代就收敛了,比固定跑 3000 次快两个数量级。

重拟合那一步用的是全部内点做一次普通 DLT,这一步能显著改善精度,因为初次的四对点只用了最小样本,对内点噪声没有任何平滑作用。但重拟合也可能把 H 拉偏——如果内点里混进了几个误差刚好卡在阈值边缘的点,它们对最小二乘的影响不小。稳妥做法是重拟合之后重新算一遍内点掩码,再迭代一次,通常两次就稳定了。

4.4 与内置函数对拍:差异应该落在哪里

手写版跑通之后,必须和fitgeotform2d做一次对拍。对拍不是为了证明谁更准,而是为了确认两件事:你的实现没有结构性错误,你的误差量级和内置版在同一档。

% 手写版 [Hmine, maskMine, st] = ransacH(m1, m2, 3, 3000); % 内置版 [tform, maskBuiltin] = fitgeotform2d(m1, m2, 'projective', ... 'MaxDistance', 3, 'MaxNumTrials', 3000, 'Confidence', 99.9); Hbuiltin = tform.T; Hbuiltin = Hbuiltin / Hbuiltin(3,3); % 两个矩阵的直接差异(注意尺度已统一) fprintf('H 矩阵最大元素差: %.4f\n', max(abs(Hmine(:) - Hbuiltin(:)))); % 更有意义的比较:在同一组点上算重投影误差 eMine = sqrt(mean(symTransferErr(Hmine, m1, m2))); eBuiltin = sqrt(mean(symTransferErr(Hbuiltin, m1, m2))); fprintf('手写版 RMSE: %.3f px,内置版 RMSE: %.3f px\n', eMine, eBuiltin); fprintf('内点数 手写 %d / 内置 %d\n', nnz(maskMine), nnz(maskBuiltin)); % 行列式符号必须为正,为负说明匹配方向或坐标顺序反了 fprintf('det(H) = %.3e,符号: %s\n', det(Hmine), ... string(sign(det(Hmine))));

判断标准是这样:如果两个 H 的元素差在 1e-2 量级以内、RMSE 相差不超过 20%、内点数差异在 10% 以内,说明两条路径是一致的,差异来自 RANSAC 的随机采样,属于正常范围。如果元素差到 0.5 以上,或者某个 H 的 det 为负,那就不是随机性了,一定有结构性 bug,回头查 4.2 的转置和 4.1 的坐标约定。

我一般还会把内点画出来看一眼分布。两个原因:内点全挤在图像某个角落,说明这一小块纹理占了主导,H 在外推区域的表现会很差;内点沿一条直线排列,说明场景里有大片直线边缘(比如楼宇的窗框),这种分布下 H 在垂直方向上的约束很弱。这两种情况即使 RMSE 很小也不能用,后面第 5 章会回到这一点。

5. 排查清单:H 算出来不对时先查这五处

单应性估计的问题有个共同特点:算法不报错,只给你一个错的结果。所以排查靠的是「现象 → 定位」,而不是看堆栈。下面五条是实际遇到的问题里频率最高的,每条都按现象、原因、解决三段写,遇到对不上的情况可以按顺序试。

5.1 warp 之后整张图被镜像或者转了 180 度

现象:imwarp出来的图和目标图完全对不上,观感上像是翻转了或者被甩到了画面外,但内点数统计还挺高(几十上百)。

原因:最常见的是匹配方向反了。matchFeatures(f1,f2)返回的idx(:,1)是图 1 的索引,idx(:,2)是图 2 的索引,如果你在调用fitgeotform2d的时候把两个参数写反了,得到的 H 是图 2 → 图 1 的映射,再用imwarp(I1, tform)去变换图 1,方向完全相反。第二种原因是坐标顺序,把[x y]当成了 `[row

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

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询