Re: fitting circle to data

[Primoz PETERLIN]

Hello,

Actually, both solutions

G = @(x,a) (x(:,1) - a(1)).^2 + (x(:,2) - a(2)).^2 - a(3).^2;
pin = [1;1;1];
[f , p] = leasqr( [xt yt], zeros(size(xt)), pin, G);

and 

A = [ 2*xt 2*yt ones(size(xt))];
b = xt.^2 + yt.^2;
c = A\b;
r = sqrt(c(3) + c(1)^2 + c(2)^2)

boil down to Kasa's algorithm (I. Kasa, "A curve fitting procedure and its error analysis", IEEE Trans. Inst. Meas., 25(1976): 8-14). The small difference between the two methods (243.62 vs 243.63 for radius for your set of data) is probably due to some numerical error.

However, with a circle arc smaller than 120˚ (and your data fits well below this margin), Kasa's algorithm is known to systematically give  too small estimates for radius. Any better method for fitting circles (Pratt, Taubin, geometric fit...) will give you an estimate around 247.42 for radius.

Let me advertise Chernov's page once again: http://www.math.uab.edu/~chernov/cl/

All the best, Primož

2009/12/1 Gastón Araguás <address@hidden>

Sorry for the noise, here ist

G = @(x,a) (x(:,1) - a(1)).^2 + (x(:,2) - a(2)).^2 - a(3).^2

[f , p] = leasqr( [xt yt], zeros(size(xt)), pin, G);
p =

  8.2224e-02
  2.4964e+02
  2.3825e+02

thanks again

2009/12/1 Gastón Araguás <address@hidden>:

> Somebody in the list told me that the function i'm using to minimize
> doesn't depends on x
> xt = %data
> yt = %data
>
> F = @(x,a) (xt-a(1)).^2+(yt-a(2)).^2-a(3).^2
>
> so i try to put it in this way
>
> F = @(x,a) (x(1)-a(1)).^2+(x(2)-a(2)).^2-a(3).^2
> pin = [1; 1; 1];
> [f , p] = leasqr([xt yt],zeros(size(xt)),pin,F);
>
> but i get:
> warning: inverse: matrix singular to machine precision, rcond = 7.27849e-18
> error: number of rows must match (1 != 388) near line 324, column 18
> error: evaluating argument list element number 1
> error: evaluating assignment _expression_ near line 324, column 2
> error: called from `leasqr' in file
> `/usr/share/octave/packages/3.0/optim-1.0.5/leasqr.m'
>
> i can't find my mistake, i think the data sizes are ok
>
> size([xt yt])
> ans =
>
>   388     2
>
> size(zeros(size(xt)))
> ans =
>
>   388     1
>
> thank you very much
>
>
>
>
> ---------- Forwarded message ----------
> From: Gastón Araguás <address@hidden>
> Date: 2009/12/1
> Subject: Re: fitting circle to data
> To: address@hidden
>
>
> 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
>
>
>
> --
> 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
>



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

-- 
Primož Peterlin,   Inštitut za biofiziko, Med. fakulteta, Univerza v Ljubljani
Lipičeva 2, SI-1000 Ljubljana, Slovenija.        address@hidden
Tel +386-1-5437612, fax +386-1-4315127,  http://biofiz.mf.uni-lj.si/~peterlin/
F8021D69 OpenPGP fingerprint: CB 6F F1 EE D9 67 E0 2F  0B 59 AF 0D 79 56 19 0F