Search code examples
image-processingcurve-fittingleast-squares

Code for a multiple quadratic (or polynomial) least squares (surface fit)?

for a machine vision project I am trying to search image data for quadratic surfaces (f(x,y) = Ax^2+Bx+Cy^2+Dy+Exy+F). My plan is to iterate through regions of data and perform a surface-fit, look at the error, see if it's a continuous surface (which would probably indicate a feature in the image).

I was previously able to find quadratic curves (f(x) = Ax^2+Bx+C) in the image data by sampling lines, by using the equations on this site

Link

this worked well, was promising, but it would be much more useful for my task to find 2-D regions that form continuous surfaces.

I see lots of articles indicating that least squares regressions scales up to multiple dimensions, but I'm not able to find code for this Hopefully there is a "closed form" (non-iterative, just compute from your data points) solution, like described above for 1D data. Does anybody know of some source or pseudocode that accomplishes this? Thanks.

(Sorry if my terminology is a bit off.)

Solution

  • I'm not sure what your background is, but if you know some linear algebra you will find linear least squares on wikipedia useful.

    Lets take the following example. Say we have the following image

    enter image description here

    and we want to know how well this fits to a 2D quadratic function in a least squares sense.

    Probably the most straightforward way to solve the problem is to compute the optimal coefficients in a least squares sense, then check the error.

    First we need to describe the matrices.

    Let X be a matrix containing every x,y coordinate in the image, taking the form

    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];

    For the example image above, X would be a 100x6 matrix.

    Let y be the image intensity values in a vector of the form

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

    In this case y is a 100 element column vector.

    We want to minimize the least squares objective function S with respect to the vector of coefficients b

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

    where |.| is the L2 norm and b is the desired coefficients

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

    Taking the vector derivative of S(b) with respect to b, setting to zero, and solving for b leads to the standard least squares solution.

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

    where inv is the matrix inverse, ' is transpose, and * is matrix multiplication.

    MATLAB example.

    % 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)

    Output

    b =

    0.0906
    0.5093
    0.1245
   -0.3733
    0.3776
  124.5412


E =

    3.4699

    In your application you would probably want to define some threshold such that when E < threshold you accept the image (or image region) as a quadratic polynomial.