多重二次(或多项式)最小二乘(表面拟合)的代码?
Code for a multiple quadratic (or polynomial) least squares (surface fit)?
对于机器视觉项目,我正在尝试搜索二次曲面 (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
这很有效,很有前途,但对于我的任务来说,找到形成连续表面的二维区域会更有用。
我看到很多文章表明最小二乘回归可以扩展到多个维度,但我找不到这方面的代码希望有一个 "closed form"(非迭代,只是从你的计算数据点)解决方案,如上文针对一维数据所述。有人知道实现此目的的某些源代码或伪代码吗?谢谢
(对不起,如果我的术语有点不对。)
我不确定你的背景是什么,但如果你了解一些线性代数,你会发现 linear least squares on wikipedia 很有用。
举个例子。假设我们有下面的图像
我们想知道这对最小二乘意义上的二维二次函数的拟合程度如何。
解决这个问题最直接的方法可能是在最小二乘意义上计算最优系数,然后检查错误。
首先我们需要描述矩阵。
设 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 元素的列向量。
我们想要最小化最小二乘 objective 函数 S
关于系数向量 b
S(b) = |y - X*b|^2
其中 |.|
是 L2 范数,b
是所需系数
b = [A;
B;
C;
D;
E;
F]
S(b)
关于b
的矢量导数,设置为零,求解b
leads 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
您接受图像(或图像区域)作为二次多项式时。
对于机器视觉项目,我正在尝试搜索二次曲面 (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
这很有效,很有前途,但对于我的任务来说,找到形成连续表面的二维区域会更有用。
我看到很多文章表明最小二乘回归可以扩展到多个维度,但我找不到这方面的代码希望有一个 "closed form"(非迭代,只是从你的计算数据点)解决方案,如上文针对一维数据所述。有人知道实现此目的的某些源代码或伪代码吗?谢谢
(对不起,如果我的术语有点不对。)
我不确定你的背景是什么,但如果你了解一些线性代数,你会发现 linear least squares on wikipedia 很有用。
举个例子。假设我们有下面的图像
我们想知道这对最小二乘意义上的二维二次函数的拟合程度如何。
解决这个问题最直接的方法可能是在最小二乘意义上计算最优系数,然后检查错误。
首先我们需要描述矩阵。
设 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 元素的列向量。
我们想要最小化最小二乘 objective 函数 S
关于系数向量 b
S(b) = |y - X*b|^2
其中 |.|
是 L2 范数,b
是所需系数
b = [A;
B;
C;
D;
E;
F]
S(b)
关于b
的矢量导数,设置为零,求解b
leads 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
您接受图像(或图像区域)作为二次多项式时。