矩阵最小二乘法求解仿射变换矩阵

作者: walkerwzy | 来源:发表于2021-03-04 02:45 被阅读0次

一个矩形三个顶点(0,0), (50, 0), (50, 50), 变换后为(30, 30), (30, 130), (130, 130), 求其仿射矩阵。

我们分别设起始和结束矩阵的坐标为： $(a_x, a_y), (b_x, b_y), (c_x, c_y)$ ，变换后的加一个prime（ $^\prime$ )符号，以此类推。
要知道，一个3X2的矩阵是不可能右乘一个矩阵得到一个3X2的矩阵（只能左乘一个3X3的），
然后，每一个新坐标，都是由原坐标的(x, y)经过变换得到(x', y‘），即使是新坐标的X值，也是需要原坐标的(x, y)值参与过来进行变化的（乘以合适的系数），然后还要加上偏移的系数，以x'为例，应该是这样： $a^\prime_x = a_x m_{00} + a_y m_{01} + m_{02}$
我们根据矩阵特征，补一个1，构造这个矩阵看看效果：
$\begin{bmatrix} \color{red}{a_x} & \color{red}{a_y} & \color{red}1 \\ b_x & b_y & 1 \\ c_x & c_y & 1 \\ \end{bmatrix} \begin{bmatrix} \color{red}{m_{00}} \\ \color{red}{m_{01}} \\ \color{red}{m_{02}} \end{bmatrix} = \begin{bmatrix} \color{red}{a^\prime_x} \\ b^\prime_x \\ c^\prime_x \end{bmatrix} \tag{红色部分即为上面的等式}$
这只是把三个x给变换出来了，其实你也可以认为这是把y给变换出来了（因为原理一样，只是系数不同）。
做到这一步，我们已经知道要如何求y坐标了，即我们只补一列的话，只能得到一个坐标的x值（或y值），要求另一半，根据坐标相乘的原理，看来只能把前三列置零，再把后三列复制进去了（这样仿射矩阵也就变成6X1了），其实就是上面矩阵乘法的重复，只不过交错一下形成x,y交错的排列：
$\begin{bmatrix} a_x & a_y & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & a_x & a_y & 1 \\ b_x & b_y & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & b_x & b_y & 1 \\ c_x & c_y & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & c_x & c_y & 1 \end{bmatrix} \begin{bmatrix} m_{00} \\ m_{01} \\ m_{02} \\ m_{10} \\ m_{11} \\ m_{12} \end{bmatrix} = \begin{bmatrix} a^\prime_x \\ a^\prime_y \\ b^\prime_x \\ b^\prime_y \\ c^\prime_x \\ c^\prime_y \\ \end{bmatrix}$
原理当然就是把第一个公式补全：
$\begin{cases} \; a^\prime_x = a_x m_{00} + a_y m{01} + m_{02} \\ \; a^\prime_y = a_x m_{10} + a_y m{11} + m_{12} \\ \\ \; b^\prime_x = b_x m_{00} + b_y m{01} + m_{02} \\ \; b^\prime_y = b_x m_{10} + b_y m{11} + m_{12} \\ \\ \; c^\prime_x = c_x m_{00} + c_y m{01} + m_{02} \\ \; c^\prime_y = c_x m_{10} + c_y m{11} + m_{12} \\ \end{cases}$
最小二乘的公式如下：
$\begin{equation} \lVert A\beta - Y \rVert{^2_2} \quad A \in \mathbb{R}^{(m\times n+1)}, \beta \in \mathbb{R}^{(n+1)\times 1}, Y \in \mathbb{R}^{m\times 1} \\ \hat \beta = (A^TA)^{-1}A^TY \end{equation}$
推导过程见此

奇异矩阵没有逆矩阵， $(A^TA)^{-1}$ 会出现无法求解的问题，也就是该方法对数据是有约束的，这个有解，另议。

我们把A和Y都做出来了，直接套用公式即可，为了编程方便，我们把前后矩阵设为A和B，仿射矩阵为M，就成了：
$M = (A^TA)^{-1}A^TB$

import numpy as np

A = [[0,0], [50, 0], [50, 50]]
B = [[30, 30], [130, 30], [130, 130]]

# 分别整理成上面分析的6x6和6x1的矩阵
# 先定义变量保留6个坐标的值
(ax, ay), (bx, by), (cx, cy) = A
(ax1, ay1), (bx1, by1), (cx1, cy1) = B

A = np.array([
    [ax, ay, 1, 0, 0, 0],
    [0, 0, 0, ax, ay, 1],
    [bx, by, 1, 0, 0, 0],
    [0, 0, 0, bx, by, 1],
    [cx, cy, 1, 0, 0, 0],
    [0, 0, 0, cx, cy, 1]
])
B = np.array([ax1, ay1, bx1, by1, cx1, cy1]).reshape(6, 1) # 比手写6X1矩阵要省事
M = np.linalg.inv(A.T @ A) @ A.T @ B # 套公式
M.reshape(2, 3)

输出：

array([[ 2.,  0., 30.],
       [ 0.,  2., 30.]])

上就是最小二乘的一个应用，也给了一篇链接介绍推导，后来我翻阅学习线代时的笔记，其实有从投影方面给的解释，直观易懂，于是另写了篇博文来介绍这个推导。

矩阵最小二乘法求解仿射变换矩阵

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读