【问题标题】:Code for a multiple quadratic (or polynomial) least squares (surface fit)?多重二次(或多项式)最小二乘法(表面拟合)的代码?
【发布时间】:2018-01-15 10:46:47
【问题描述】:

对于机器视觉项目,我正在尝试搜索二次曲面 (f(x,y) = Ax^2+Bx+Cy^2+Dy+Exy+F) 的图像数据。我的计划是遍历数据区域并执行表面拟合,查看错误,看看它是否是连续表面(这可能表示图像中的特征)。

我之前能够通过使用本网站上的方程式,通过采样线在图像数据中找到二次曲线 (f(x) = Ax^2+Bx+C)

http://mathforum.org/library/drmath/view/72047.html

这很好用,很有希望,但对于我的任务来说,寻找形成连续表面的二维区域会更有用。

我看到很多文章表明最小二乘回归可以扩展到多个维度,但我无法找到此代码希望有一个“封闭形式”(非迭代,只需从您的数据点计算)解决方案,就像上面描述的一维数据。有没有人知道一些源代码或伪代码来完成这个?谢谢。

(对不起,如果我的术语有点不对。)

【问题讨论】:

  • 你写的是f(x,y) = Ax^2+Bx+Cy^2+D+E,但你的意思是f(x,y) = Ax+Bx^2+Cy+Dy^2+Exy+F?
  • 糟糕,谢谢。

标签: image-processing curve-fitting least-squares


【解决方案1】:

我不确定你的背景是什么,但如果你了解一些线性代数,你会发现linear least squares on wikipedia 很有用。

让我们举个例子。假设我们有以下图像

我们想知道这与最小二乘法意义上的 2D 二次函数的拟合程度。

可能最直接的解决问题的方法是在最小二乘意义上计算最优系数,然后检查错误。

首先我们需要描述矩阵。

让X 是一个矩阵,包含图像中的每个 x,y 坐标,采用以下形式

X = [x1  x1^2  y1  y1^2  x1*y1  1;
     x2  x2^2  y2  y2^2  x2*y2  1;
     ...
     xN  xN^2  yN  yN^2  xN*yN  1];

对于上面的示例图像,X 将是一个100x6 矩阵。

设y 为形式向量中的图像强度值

y = [img(x1,y1);
     img(x2,y2);
     ...
     img(xN,yN)]

在这种情况下,y 是一个 100 元素的列向量。

我们希望最小二乘目标函数S 相对于系数向量b

S(b) = |y - X*b|^2

其中|.| 是 L2 范数,b 是所需系数

b = [A;
     B;
     C;
     D;
     E;
     F]

取S(b)对b的向量导数,置零,求解bleads to the standard least squares solution。

b = inv(X'X)*X'*y

其中inv 是逆矩阵,' 是转置,* 是矩阵乘法。

MATLAB 示例。

% Generate an image

% define x,y coordinates for each location in the image
[x,y] = meshgrid(1:10,1:10);

% true coefficients
b_true = [0.1 0.5 0.3 -0.4 0.4 124];
% magnitude of noise
P = 2;
% create image
img = b_true(1).*x + b_true(2).*x.^2 + b_true(3).*y + b_true(4).*y.^2 + b_true(5).*x.*y + b_true(6);
noise = P*randn(10,10);
img = img + noise;

% Begin least squares optimization

% create matrices
X = [x(:) x(:).^2 y(:) y(:).^2 x(:).*y(:) ones(size(x(:)))];
y = img(:);

% estimated coefficients
b = (X.'*X)\(X.')*y

% mean square error (expected to be near P^2)
E = 1/numel(y) * sum((y - X*b).^2)

输出

b =

    0.0906
    0.5093
    0.1245
   -0.3733
    0.3776
  124.5412


E =

    3.4699

在您的应用程序中,您可能希望定义一些阈值,这样当E < threshold 您接受图像(或图像区域)作为二次多项式时。

【讨论】:

  • 天哪……我怎么没看到这个??非常感谢您的详尽回答!现在消化它......
猜你喜欢
  • 2019-07-31
  • 2012-07-14
  • 1970-01-01
  • 1970-01-01
  • 2013-01-21
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多