联合方程的最小二乘解

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的内置功能可能没有帮助!

编辑:我有一大套旧(x,y)新(x,y)适合.问题是超定的(比未知数更多的数据点)如此简单的矩阵求逆是不够的 - 正如我所说,我真的需要拟合上的误差.

Answer 1

下面的代码应该可以解决问题。我使用以下公式计算残差：

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

然后根据Wolfram 的 MathWorld 中描述的一般过程导出最小二乘公式。

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

由于没有对输出数据应用随机噪声，因此该程序会生成与输入参数相同的四个参数（ P、Q、R和），并且值为零。SrSquared

随着越来越多的随机噪声施加到输出点，常数开始偏离正确值，并且该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;
}