单应性矩阵

单应性矩阵

1.1 核心概念

单应性矩阵处理的是平面到平面的变换问题。想象一张倾斜拍摄的文档照片，我们想要将其校正为正视图，这就是典型的单应性矩阵应用。

1. 变换层级关系
   平面变换的自由度从低到高、覆盖范围从窄到宽的层级为：

刚性变换（Rigid）⊂ 相似变换（Similarity）⊂ 仿射变换（Affine）⊂ 单应变换（Homography）⊂ 投影变换

2. 矩阵形式的差异与联系
   两者均用齐次坐标描述平面点的变换（平面点 \((x,y)\)表示为齐次坐标 \((x,y,1)\)），但矩阵的约束和自由度不同。

1.2 仿射变换矩阵（Affine）

维度：\(3 \times 3\) 矩阵，最后一行固定为 \([0, 0, 1]\)（无透视项）
数学形式：

\[\mathbf{H_{affine}} = \begin{bmatrix} a_1 & a_2 & t_x \\ a_3 & a_4 & t_y \\ 0 & 0 & 1 \end{bmatrix} \]

变换公式：齐次坐标下，目标点\(p_1 =H_{affine} \cdot p\),展开为：

\[ u = a_1x + a_2y + t_x \\ y = a_3x + a_4y + t_y \]

1.3 单应性矩阵（Homography）

维度：\(3 \times 3\) 矩阵，无固定行约束（包含透视项），且矩阵是齐次的（整体缩放不影响变换结果，通常归一化使最后一个元素为 1）

数学表达式：

\[\mathbf{x} \sim \mathbf{H} \cdot \mathbf{X} \]

\[\begin{bmatrix} u \\ v \\ 1 \end{bmatrix} \sim \begin{bmatrix} h_{11} & h_{12} & h_{13} \\ h_{21} & h_{22} & h_{23} \\ h_{31} & h_{32} & h_{33} \end{bmatrix} \begin{bmatrix} X \\ Y \\ 1 \end{bmatrix} \]

其中：

• \(\mathbf{X} = (X, Y, 1)^T\) 是世界平面上的齐次坐标
• \(\mathbf{x} = (u, v, 1)^T\) 是图像平面上的齐次坐标
• "\(\sim\)" 表示比例相等（齐次坐标的性质）

展开形式：

\[u = \frac{h_{11}X + h_{12}Y + h_{13}}{h_{31}X + h_{32}Y + h_{33}}, \quad v = \frac{h_{21}X + h_{22}Y + h_{23}}{h_{31}X + h_{32}Y + h_{33}} \]

其中单应矩阵H有几个自由度？
其实只有8个，因为这里使用的是齐次坐标系（可以进行任意尺度的缩放）；比如把\(h_{ij}\)乘以任意一个非零常数k,并不会改变等式结果。

\[u = \frac{kh_{11}X + kh_{12}Y + kh_{13}}{kh_{31}X + kh_{32}Y + kh_{33}}, \quad v = \frac{kh_{21}X + kh_{22}Y + kh_{23}}{kh_{31}X + kh_{32}Y + kh_{33}} \]

与上述等式是一样的。

1.4 求解方法

求解条件

• 至少需要4对不共线的对应点
• 每对点提供2个独立方程
• 单应性矩阵有8个自由度（尺度无关）

8自由度下H计算过程有两种方法：

1. 直接设置\(h_{33} = 1\),那么上述等式变为：

\[u = \frac{h_{11}X + h_{12}Y + h_{13}}{h_{31}X + h_{32}Y + 1}, \quad v = \frac{h_{21}X + h_{22}Y + h_{23}}{h_{31}X + h_{32}Y + 1} \]

2. 将H添加约束条件，将H矩阵模变成1，如下：

\[h_{11}^2 + h_{12}^2 + h_{13}^2+ h_{21}^2 + h_{22}^2 + h_{23}^2 + h_{31}^2 + h_{32}^2 + h_{33}^2 = 1 \]

将上述等式分别乘以分母展开，得到：

\[(h_{31}X + h_{32}Y + h_{33})u = h_{11}X + h_{12}Y + h_{13} \]

\[(h_{31}X + h_{32}Y + h_{33})v = h_{21}X + h_{22}Y + h_{23} \]

整理，得到：

\[h_{11}X + h_{12}Y + h_{13} - h_{31}Xu -h_{32}Yu - h_{33}u = 0 \\ h_{21}X + h_{22}Y + h_{23} - h_{31}Xu -h_{32}Yu - h_{33}u = 0 \]

那么可以得到如下线性方程组：

\[\begin{bmatrix} x & y & 1 & 0 & 0 & 0 & -xu & -yu & -u \\ 0 & 0 & 0 & x & y & 1 & -xv & -yv & -v \end{bmatrix} \begin{bmatrix} h_{11} \\ h_{12} \\ h_{13} \\ h_{21} \\ h_{22} \\ h_{23} \\ h_{31} \\ h_{32} \\h_{33} \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} \]

1.5 实际应用案例

场景：文档扫描校正

输入数据：

世界坐标系（A4纸角点，单位：毫米）：
(0, 0), (210, 0), (210, 297), (0, 297)

图像坐标系（像素坐标）：
(120, 150), (380, 160), (390, 520), (110, 510)

求解过程详解

第一步：构建完整系数矩阵

对于4个点，构建8×9的系数矩阵：

\[\mathbf{A} = \begin{bmatrix} 0 & 0 & 1 & 0 & 0 & 0 & -0*120 & -0*120 & -120 \\ 0 & 0 & 0 & 0 & 0 & 1 & -0*150 & -0*150 & -150 \\ 210 & 0 & 1 & 0 & 0 & 0 & -380×210 & -380×0 & -380 \\ 0 & 0 & 0 & 210 & 0 & 1 & -160×210 & -160×0 & -160 \\ 210 & 297 & 1 & 0 & 0 & 0 & -390×210 & -390×297 & -390 \\ 0 & 0 & 1 & 210 & 297 & 1 & -520×210 & -520×297 & -520 \\ 0 & 297 & 1 & 0 & 0 & 0 & -110×0 & -110×297 & -110 \\ 0 & 0 & 1 & 0 & 297 & 1 & -510×0 & -510×297 & -510 \end{bmatrix} \]

第二步：SVD分解求解

对矩阵 \(\mathbf{A}\) 进行奇异值分解：\(\mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T\)

最小二乘解对应于 \(\mathbf{V}\) 的最后一列（最小奇异值对应的右奇异向量）。

第三步：获得单应性矩阵

通过SVD求解得到的向量（归一化后）：

\[\mathbf{h} = [0.85, -0.02, 120.5, 0.03, 1.20, 148.2, 0.0001, 0.0002, 1.0]^T \]

第四步：构建单应性矩阵

\[\mathbf{H} = \begin{bmatrix} 0.85 & -0.02 & 120.5 \\ 0.03 & 1.20 & 148.2 \\ 0.0001 & 0.0002 & 1.0 \end{bmatrix} \]

验证求解结果

测试点：世界坐标 (105, 148.5)（A4纸中心）

使用单应性矩阵变换：

\[\begin{bmatrix} 0.85 & -0.02 & 120.5 \\ 0.03 & 1.20 & 148.2 \\ 0.0001 & 0.0002 & 1.0 \end{bmatrix} \begin{bmatrix} 105 \\ 148.5 \\ 1 \end{bmatrix} = \begin{bmatrix} 0.85×105 + (-0.02)×148.5 + 120.5 \\ 0.03×105 + 1.20×148.5 + 148.2 \\ 0.0001×105 + 0.0002×148.5 + 1.0 \end{bmatrix} = \begin{bmatrix} 206.78 \\ 329.55 \\ 1.0402 \end{bmatrix} = \begin{bmatrix} 198.79 \\ 316.78 \\ 1 \end{bmatrix} \]

重投影误差检查：

• 理论图像坐标：(198.79, 316.78)
• 实际检测坐标：(约200, 约320)
• 误差：\(\sqrt{(200-198.79)^2 + (320-316.78)^2} = \sqrt{1.46 + 10.37} = 3.44\)

1.6 编程调用

OpenCV 计算单应矩阵
OpenCV 中核心函数是 findHomography，支持最小二乘、RANSAC 等鲁棒算法，适配不同场景（如外点较多的情况）。

Mat findHomography(
    InputArray srcPoints,  // 源点集（至少4对不共线点，如vector<Point2f>）
    InputArray dstPoints,  // 目标点集（数量/顺序与源点一一对应）
    int method = 0,        // 计算方法：
                           // 0：默认，最小二乘（无鲁棒性，易受外点影响）
                           // RANSAC：鲁棒估计算法（推荐，可过滤外点）
                           // LMEDS：中值鲁棒估计
    double ransacReprojThreshold = 3,  // RANSAC重投影阈值（像素），超过则视为外点
    OutputArray mask = noArray()       // 输出掩码（内点=1，外点=0）
);

示例：

#include <opencv2/opencv.hpp>
using namespace cv;
using namespace std;

// 1. 准备源点和目标点（示例：4对不共线点）
vector<Point2f> srcPoints = {
   Point2f(10, 10), Point2f(100, 10), 
   Point2f(100, 100), Point2f(10, 100)
};
vector<Point2f> dstPoints = {
   Point2f(20, 20), Point2f(110, 15), 
   Point2f(105, 110), Point2f(15, 105)
};

// 2. 计算单应矩阵（RANSAC鲁棒估计）
Mat H, mask;
H = findHomography(srcPoints, dstPoints, RANSAC, 3.0, mask);

Halcon 计算单应矩阵
Halcon 中计算单应矩阵分两种场景：

• 无鲁棒性（最小二乘）：vector_to_hom_mat2d；
• 有鲁棒性（RANSAC 过滤外点）：ransac_homography_2d（推荐）

// 最小二乘，无鲁棒性
vector_to_hom_mat2d( : : Row1, Column1, Row2, Column2 : HomMat2D)
// Row1/Column1：源点的行 / 列坐标（Halcon 中Row对应 y，Column对应 x，）
// Row2/Column2：目标点的行 / 列坐标；
// HomMat2D：输出 3x3 单应矩阵。