RE: fitting circle to data

[Richardson, Anthony]

I'm not sure what you mean when you say that this is a linear version:

> A = [ xt yt ones(size(xt))];
> b = xt.^2 + yt.^2;
> c = A\b;

This is an overdetermined system (at least when length(xt) > 3).  
In that case,  the '\' operator gives the least squares solution  
(see the documentation on the '\' operator under the Arithmetic 
Ops section of the documentation).

c = A\b  is equivalent to c = (A'*A)^(-1)*A'*b  (the pseudoinverse
solution).  

The 3-element c vector contains the least squares values of the circle
x and y coordinates (x0, y0) and the third element is the least squares value
of r^2 - x0^2 - y0^2.

I suspect that the leasqr routine uses an iterative method to compute
the least squares solution and that is the reason for the difference
in answers from the two methods.

Tony Richardson

> -----Original Message-----
> From: address@hidden [mailto:address@hidden On Behalf Of
> Gastón Araguás
> Sent: Tuesday, December 01, 2009 7:45 AM
> To: address@hidden
> Subject: Re: fitting circle to data
> 
> Thank you very much to all. This code is working well now:
> 
> xt = %data
> yt = %data
> 
> F = @(x,a) (xt-a(1)).^2+(yt-a(2)).^2-a(3).^2
> pin = [1;1;1];
> [f , p] = leasqr(xt,yt,pin,F);
> 
> p =
> 
>    8.2224e-02
>    2.4914e+02
>    2.3772e+02   <---- r
> 
> on the other hand this linear version gives good results too, without
> initializations (thanks to Richardson)
> 
> A = [ xt yt ones(size(xt))];
> b = xt.^2 + yt.^2;
> c = A\b;
> 
> r = sqrt(c(3) + c(1)^2 + c(2)^2);
> r =  238.25
> 
> I'm not sure wich of them is the best result, I'll try to test it with
> real and artifcial data
> 
> Thank you all again
> 
> --
> Gastón Araguás
> ______________________________________________________
> CIII - Centro de Investigación en Informática para la Ingeniería
> Univ. Tecnológica Nacional Facultad Regional Córdoba. Argentina
> 
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://www-old.cae.wisc.edu/mailman/listinfo/help-octave