矩阵最小二乘法求解仿射变换矩阵
一个矩形三个顶点(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{aligned} &\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{aligned} $$
奇异矩阵没有逆矩阵,$(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.]])
上就是最小二乘的一个应用,也给了一篇链接介绍推导,后来我翻阅学习线代时的笔记,其实有从投影方面给的解释,直观易懂,于是另写了篇博文来介绍这个推导。