最小二乘与仿射变换(附js代码)

最新推荐文章于 2024-03-11 14:07:37 发布

snoopy_21

最新推荐文章于 2024-03-11 14:07:37 发布

阅读量1.2k

点赞数

分类专栏：数学 # js

本文链接：https://blog.csdn.net/qq_29598161/article/details/109523218

版权

数学同时被 2 个专栏收录

6 篇文章 0 订阅

订阅专栏

4 篇文章 0 订阅

订阅专栏

仿射变换(affine transform)

在这里插入图片描述
仿射变换的变换矩阵统一用 $\left[\begin{array}{lll}a & b & c \\ d & e & f \\ 0 & 0 & 1\end{array}\right]$ 来描述，不同基础变换的 $a, b, c, d, e, f$ 约束不同。所以有以下关系：
$\left[\begin{array}{l} x^{\prime} \\ y^{\prime} \\ 1 \end{array}\right]=\left[\begin{array}{lll} m_{11} & m_{12} & t_{x}\\ m_{21} & m_{22} & t_{x} \\ 0 & 0 & 1 \end{array}\right]\left[\begin{array}{l} x \\ y \\ 1 \end{array}\right]$
注意：当位移值 $t_{x}$ 和 $t_{y}$ 位于第3列时，上面变换矩阵是左乘。
也可以表示为：
$\left[\begin{array}{lll} x & y & 1 \end{array}\right] \times\left[\begin{array}{ccc} m_{11} & m_{21} & 0 \\ m_{12} & m_{22} & 0 \\ t_{x} & t_{y} & 1 \end{array}\right]=\left[\begin{array}{lll} x^{\prime} & y^{\prime} & 1 \end{array}\right]$

位移值 $t_{x}$ 和 $t_{y}$ 位于第3行时，上面变换矩阵是右乘。可以看到变换矩阵其实是进行了转置。在不同库中，例如threejs：2维向量应用3 x 3矩阵，其源码使用的是左乘，所以我们表示时使用前者。

在这里插入图片描述
注意：上面旋转矩阵的 $sin\theta$ 前的负号情况。在常规坐标系下，顺时针旋转角度为正才这么表示。我们更常用的表示是逆时针为正，所以有：

$\left[\begin{array}{ccc} \cos \theta & -\sin \theta & 0 \\ \sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{array}\right]$
更多参考：仿射变换及其变换矩阵的理解

在本文，将仿射变换表示为：
$\left[\begin{array}{ccc} M_{11} & M_{12} & M_{13} \\ M_{21} & M_{22} & M_{23} \\ 0 & 0& 1 \\ \end{array}\right]$
在求解时表示为：
$\left[\begin{array}{ccc} M_{11} \\ M_{12} \\ M_{13} \\ M_{21} \\ M_{22} \\ M_{23} \\ \end{array}\right]$
在给出3对数据点个时(不共线)，我们可以列出6个方程组
$M_{11}x_{0} + M_{11}y_{0} + M_{13} = x_{p0} \\ M_{21}x_{0} + M_{21}y_{0} + M_{23} = y_{p0} \\ M_{11}x_{1} + M_{11}y_{1} + M_{13} = x_{p1} \\ M_{21}x_{1} + M_{21}y_{1} + M_{23} = y_{p1} \\ M_{11}x_{2} + M_{11}y_{2} + M_{13} = x_{p2} \\ M_{21}x_{2} + M_{21}y_{2} + M_{23} = y_{p2} \\$

$Mx = x_p$
$x$ 表示原始坐标， $M$ 表示仿射变换矩阵， $x_p$ 表示变换后坐标

编程中把变换矩阵M描述为待求量X，表示：

[x0 y0 1  0  0  0 ]     [M11]   [xp0]
[0  0  0  x0 y0 1 ]     [M12]   [yp0] 
[x1 y1 1  0  0  0 ]     [M13]   [xp1]
[0  0  0  x1 y1 1 ]  *  [M21] = [yp1] 
[x2 y2 1  0  0  0 ]     [M22]   [xp2] 
[0  0  0  x2 y2 1 ]     [M23]   [yp2]
        A            *    X   =   B

3个点不共线时，A满秩，所|A|不为0，直接对A求逆，便可解出X，得到仿射变换的6个自由度值。

下面使用最小二乘分析数据点大于和小于3时，求解方法。

最小二乘

线性最小二乘法是将数据与未知参数之间的关系以线性形式表达的一种拟合数据模型的方法。

$\begin{array}{l} A x=b \\ J = \underset{x}{\operatorname{argmin}}\|A x-b\|^{2}\\ \end{array}$

取极值，求导要为0

$\frac{\partial J}{\partial a}=A^{T}(A x-b)=0\\ A^{T}A x = A^{T}b\\ x = (A^{T} A)^{-1} A^{T} b$

A的秩>3时， $A^{T} A|$ 不为0，所以可以求得 $A^{T} A)^{-1}$ ，也就有 $x$
A的秩<2时， $A^{T} A|$ 可能为0，需要加正则项来确定唯一。

以下示例展示的matlab代码表明：当数据点刚好为3个时，且不共线，秩为3，可以直接求逆：

function M= affine_least_square(x0,y0, x1,y1, x2,y2, xp0,yp0, xp1,yp1, xp2,yp2)
% This function finds the affine transform between three points

% an affine transformation with a 3x3 affine transformation matrix: 
% 
% [M11 M12 M13]
% [M21 M22 M23]
% [M31 M32 M33]

%Since the third row is always [0 0 1] we can skip that.

% [x0 y0 1  0  0  0 ]     [M11]   [xp0]
% [0  0  0  x0 y0 1 ]     [M12]   [yp0] 
% [x1 y1 1  0  0  0 ]     [M13]   [xp1]
% [0  0  0  x1 y1 1 ]  *  [M21] = [yp1] 
% [x2 y2 1  0  0  0 ]     [M22]   [xp2] 
% [0  0  0  x2 y2 1 ]     [M23]   [yp2]
                          
%         A            *    X   =  B
%A -> 6x6
%X -> 6x1 in fact 
%X=A\B


A=[x0 y0 1  0  0  0 ; 0  0  0  x0 y0 1 ; x1 y1 1  0  0  0 ; 0  0  0  x1 y1 1 ; x2 y2 1  0  0  0;0  0  0  x2 y2 1];
B=[xp0; yp0; xp1;  yp1; xp2 ; yp2 ];


% X=A\B; 
% X=inv(A)*B;
X=pinv(A)*B;

M11 =X(1,1);
M12 =X(2,1);  
M13 =X(3,1);
M21 =X(4,1);
M22 =X(5,1);
M23 =X(6,1);
M=[ M11 M12 M13; M21 M22 M23; 0 0 1];
end

当3个数据点共线时，其秩变为2，行列式值为0，即不可逆。但是可以求A的伪逆(广义逆矩阵)，得到其中的一个解。当数据点多余3个且秩大于3时，便不能求得A的逆矩阵，此时也需要求伪逆来解。

线性最小二乘法是将数据与未知参数之间的关系以线性形式表达的一种拟合数据模型的方法。

$\begin{array}{l} A x=b \\ X^{*}=\underset{x}{\operatorname{argmin}}\|A x-b\|^{2}=\left(A^{T} A\right)^{-1} A^{T} b \end{array}$

import * as math from "mathjs";

export default class MathUtils{
    constructor(){

    }
    /**
     * Get generalized inverse use QR decomposition
     * @param A arbitrary matrix
     * @return matrix:PInvA of APInvAA=A | APInvA=I
     */
    static pInvQR(A: math.Matrix | math.MathArray): math.Matrix | math.MathArray{
        // more reference:
        // https://blog.csdn.net/xingozd/article/details/50417233
        const QR = math.qr(A);
        const Q = QR.Q;
        const R = QR.R;
        const RT = math.transpose(R);
        const QT = math.transpose(Q);
        const InvR = math.multiply(math.inv(math.multiply(RT, R)), RT);
        const PInvA = math.multiply(InvR, QT);
        return PInvA;
    }

    /**
     * This function finds the affine transform between three points
     * an affine transformation with a 3x3 affine transformation matrix:
     * [M11 M12 M13]
     * [M21 M22 M23]
     * [M31 M32 M33]
     * Since the third row is always [0 0 1] we can skip that.
     * [x0 y0 1  0  0  0 ]     [M11]   [xp0]
     * [0  0  0  x0 y0 1 ]     [M12]   [yp0]
     * [x1 y1 1  0  0  0 ]     [M13]   [xp1]
     * [0  0  0  x1 y1 1 ]  *  [M21] = [yp1]
     * [x2 y2 1  0  0  0 ]     [M22]   [xp2]
     * [0  0  0  x2 y2 1 ]     [M23]   [yp2]
     *         A            *    X   =  BapplyMatrix3
     * @param pointPairs [A B]
     * @return matrix:x of Ax = B
     */
    static affineLeastSquare(pointPairs: [[number, number], [number, number]][]): math.Matrix{
        // more reference:
        // http://ros-developer.com/2017/12/26/finding-affine-transform-with-linear-least-squares/
        // https://www.cnblogs.com/AndyJee/p/5053354.html

        let A: number[][] = [];
        let B: number[] = [];

        pointPairs.forEach(pointPair => {
            let p1 = pointPair[0];
            let p2 = pointPair[1];
            A.push([p1[0], p1[1], 1, 0, 0, 0]);
            A.push([0, 0, 0, p1[0], p1[1], 1]);
            B.push(p2[0]);
            B.push(p2[1]);
        });


        // (A^TA)^(-1)
        let AT = math.transpose(A);
        let ATA = math.multiply(AT, A);

        let ATAInv: math.Matrix | math.MathArray;
        try {
            ATAInv = math.inv(ATA) as math.MathArray
        }catch (e) {
            ATAInv = this.pInvQR(ATA)//这里应该使用正则项求解，要么直接在最开始就是用pinv求解
        }

        // X = (A^TA)^(-1)A^TB
        let X = math.multiply(math.multiply(ATAInv, AT), B) as number[];

        let x = math.matrix([
            [X[0], X[1], X[2]],
            [X[3], X[4], X[5]],
            [0, 0, 1]
        ]);
        return x;
    }

}

伪逆

在这里插入图片描述

伪逆矩阵和最小二乘估计

在这里插入图片描述

求伪逆

求伪逆的三种方法：直接，SVD，QR及具体的应用
https://blog.csdn.net/xingozd/article/details/50417233?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-2.channel_param&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-2.channel_param

C++SVD分解求伪逆（Eigen库）（附C++代码）
https://blog.csdn.net/weixin_44344462/article/details/88879356

js计算

https://mathnotepad.com/

snoopy_21

关注

0
点赞
踩
7

收藏

觉得还不错? 一键收藏
1
评论
最小二乘与仿射变换(附js代码)

伪逆矩阵和最小二乘估计相关性质矩阵的逆、伪逆、左右逆，最小二乘，投影矩阵https://blog.csdn.net/baidu_38172402/article/details/82931879?utm_medium=distribute.pc_relevant.none-task-blog-title-2&spm=1001.2101.3001.4242伪逆矩阵（广义逆矩阵https://www.jianshu.com/p/609fa0cce409含义：伪逆矩阵的意义及求..
复制链接

扫一扫