algorithm - 联立方程的最小二乘解

Question

我正在尝试适应从一组坐标到另一组坐标的转换。

x' = R + Px + Qy
y' = S - Qx + Py
Where P,Q,R,S are constants, P = scale*cos(rotation). Q=scale*sin(rotation)

有一个众所周知的“手动”公式，用于将 P、Q、R、S 拟合到一组对应点。 但我需要对拟合进行误差估计 - 所以我需要一个最小二乘解。

阅读“数值食谱”，但我无法弄清楚如何对其中包含 x 和 y 的数据集执行此操作。

任何人都可以指出如何执行此操作的示例/教程/代码示例吗？
不太在意语言。
但是 - 仅使用 Matlab/Lapack/numpy/R 的内置功能可能没有帮助!

编辑: 我有一大套 old(x,y) new(x,y) 适合。问题是超定的(数据点多于未知数)，因此简单的矩阵求逆是不够的 - 正如我所说，我确实需要拟合误差。

Answer 1

下面的代码应该可以解决问题。我对残差使用了以下公式:

residual[i] =   (computed_x[i] - actual_x[i])^2
              + (computed_y[i] - actual_y[i])^2

然后根据general procedure推导出最小二乘公式在 Wolfram 的 MathWorld 中描述。

我在 Excel 中测试了这个算法，它按预期执行。我使用了十个随机点的集合，然后通过随机生成的变换矩阵对其进行旋转、平移和缩放。

在没有对输出数据应用随机噪声的情况下，该程序产生四个参数(P、Q、R 和 S) 与输入参数相同，rSquared 值为零。

随着越来越多的随机噪声应用于输出点，常量开始偏离正确值，rSquared 值相应增加。

代码如下:

// test data
const int N = 1000;
float oldPoints_x[N] = { ... };
float oldPoints_y[N] = { ... };
float newPoints_x[N] = { ... };
float newPoints_y[N] = { ... };

// compute various sums and sums of products
// across the entire set of test data
float Ex =  Sum(oldPoints_x, N);
float Ey =  Sum(oldPoints_y, N);
float Exn = Sum(newPoints_x, N);
float Eyn = Sum(newPoints_y, N);
float Ex2 = SumProduct(oldPoints_x, oldPoints_x, N);
float Ey2 = SumProduct(oldPoints_y, oldPoints_y, N);
float Exxn = SumProduct(oldPoints_x, newPoints_x, N);
float Exyn = SumProduct(oldPoints_x, newPoints_y, N);
float Eyxn = SumProduct(oldPoints_y, newPoints_x, N);
float Eyyn = SumProduct(oldPoints_y, newPoints_y, N);

// compute the transformation constants
// using least-squares regression
float divisor = Ex*Ex + Ey*Ey - N*(Ex2 + Ey2);
float P = (Exn*Ex + Eyn*Ey - N*(Exxn + Eyyn))/divisor;
float Q = (Exn*Ey + Eyn*Ex + N*(Exyn - Eyxn))/divisor;
float R = (Exn - P*Ex - Q*Ey)/N;
float S = (Eyn - P*Ey + Q*Ex)/N;

// compute the rSquared error value
// low values represent a good fit
float rSquared = 0;
float x;
float y;
for (int i = 0; i < N; i++)
{
    x = R + P*oldPoints_x[i] + Q*oldPoints_y[i];
    y = S - Q*oldPoints_x[i] + P*oldPoints_y[i];
    rSquared += (x - newPoints_x[i])^2;
    rSquared += (y - newPoints_y[i])^2;
}