Least Squares Fitting Method with a Circle

Hi everyone,

Forgive me, I am no expert at MATLAB. I would appreciate it greatly if someone could explain to me the method of nonlinear least squares and how to fit it with a circle of random points. Here's what I have so far:

% radius information-------------------------------------------------------
r = 10;                 % constant radius
Rc = @(x,y,z) (x.^2 + y.^2 + z.^2).^(1/2);      % radius as a function
n = rand(1000,3);       % generate random points on nominal circle
xc = n(:,1);
yc = n(:,2);            % separate into coordinates
zc = n(:,3);
%--------------------------------------------------------------------------

% angular information------------------------------------------------------
theta = linspace(0,2*pi,1000);      % plot bounds
theta = theta';
%--------------------------------------------------------------------------

% generate random radii----------------------------------------------------
R = Rc(xc,yc,zc);
R(:,2) = R(:,1);       % same radii in all directions
R(:,3) = R(:,1);
p = n + (r/100)*R;     % generate plot points around the circle
figure(1)
polar(theta,Rc(xc,yc,zc))
%--------------------------------------------------------------------------

% Least Squares fitting----------------------------------------------------
e = dist(p,R');
perf = sse(e);
lsqrfit = lsqr(e,Rc(xc,yc,zc),[],1000);
%--------------------------------------------------------------------------

So from what I understand, the residuals (errors) need to be found and then the sum of squares of the residuals need to be found. After that I am not sure where to go. Thanks for all the help.

Ian

[Andrew Newell]

I'm having trouble figuring out what this code is supposed to be doing. First, the data: are they random points on a circle, or scattered near a circle? Can the center and the radius of the circle vary? What is your fitting criterion? (I would think sum of squares of the distance from the circle would be best).

[Ian Wood]

They are supposed to be scattered near a circle, and then fitted with the least squares method. The location of the center of the circle doesn't matter.

If you look at the last section of my code entitled "Least Squares Fitting", that's what I'm attempting to do (I think..) The "dist" function gives me the minimum distances, "sse" is the sum squared of the errors (residuals) and "lsqr" is supposed to solve for the least squares parameters, though I am unsure if the parameters in that function are correct to do what I intend. Note that they are all pre-defined functions within the library.

Answer 1

[bym]

Maybe this will help, it is adapted from "measuring the radius of a roll of tape" demo in the image processing toolbox:

r = 5+.5*rand(100,1);         %100 random radii
ang = linspace(0,2*pi,100)';  % angles
x = r.*cos(ang)+2; y = r.*sin(ang)+3; % x,y coordinates centered [2,3]

plot(x,y,'+')
axis equal; hold on;

c = [x y ones(length(x),1)]\-(x.^2+y.^2); %least squares fit

xhat = -c(1)/2;
yhat = -c(2)/2;
rhat = sqrt(xhat^2+yhat^2-c(3));

plot(rhat*cos(ang)+xhat,rhat*sin(ang)+yhat,...
    'g','linewidth',2) %best fit circle

[Ian Wood]

Thanks! I was able to modify this to fit what I wanted in my code.

[Rajdeep Chowdhury]

Hey Bym,

Can I use same type of codes to fit a circle with a fixed radius? Thanks.

[Image Analyst]

For a fixed radius, where would you want the center? At the mean of the point?

xCenter = mean(x);
yCenter = mean(y);

[Rajdeep Chowdhury]

I want a least square type fitting to obtain the center (x,y) when the radius is predefined.

Answer 2

[Daniel Shub]

Have you looked at the FEX: http://www.mathworks.com/matlabcentral/fileexchange/?term=fit+chernov

[Ian Wood]

No I haven't, but that should help me somewhat. Thanks!

[Ian Wood]

It's the right idea, but I want to write my own code for the fit. I noticed that none of the files used the polar plot. Is there a reason for this, or am I able to do the fit on polar coordinates?

Answer 3

[Ian Wood]

Is there a problem with my code? Like maybe I shouldn't use a polar plot in this case, or something that could be done better.