# Program to fit a circle
# acooding to the minimization of the
# sum of squares of the algebraic distance
#
# Fraga, 19/06/2010
 
# Idea based in papers:
# 
# V. Pratt. Direct least-squares fitting of algebraic surfaces
# ACM SIGGRAPH Computer Graphics, (21) 4, 1987, pp. 145 - 152
# 
# A. Fitzgibbon, M. Pilu, and R.B. Fisher,
# Direct least square fitting of ellipses
# IEEE Trans, of Pattern Anal. & Mach. Intell. (21) 5, 1999,  pp. 476-480
# 

if ( nargin != 1 )
	printf ("Args: file\n");
	exit (1);
endif

file = char( argv() );
# file = char( argv(1) );
# printf ( "# %s\n", file );

# Format of input file
# number of 2d points
# x1 y1
# x2 y2
# .  .
# .  .
# xn yn
[FILE1, msg] = fopen( file, "r" );
if  ( FILE1 < 0 )
	# printf ( "%s\n", msg );
	exit (1);
endif

[n, count] = fscanf (FILE1, "%d", 1 );
# printf ( "%d\n", n );

S = zeros( 4, 4); 
X = zeros( 4, 1);
for i=1:n
	[x, count] = fscanf (FILE1, "%f %f", 2);
	# printf ( "%f %f\n", x(1), x(2) );
	X(1,1) = x(1)*x(1) + x(2)*x(2);
	X(2,1) = x(1);
	X(3,1) = x(2);
	X(4,1) = 1.0;
	S += X * X';
endfor
# S

fclose( FILE1 );


# Constraint matrix
C = [ 0  0  0 -2;
      0  1  0  0;
      0  0  1  0;
     -2  0  0  0];

[P, D2] = eig( S );

D2i = zeros( 4, 4); 
for i=1:4
	D2i(i,i) = 1/sqrt( D2(i,i) );
endfor
# D2i

M = D2i * P' * C * P * D2i;
 
[Y, D0] = eig( M );

X = P*D2i*Y;

k = 1;
for ii=2:4
	if ( D0(ii,ii) > D0(k,k) )
		k = ii;
	end
endfor

# Eigenvector k is the result
x0 = -X(2,k)/(2*X(1,k));
y0 = -X(3,k)/(2*X(1,k));

r = sqrt( (X(2,k)*X(2,k) + X(3,k)*X(3,k))/(4.0*X(1,k)*X(1,k)) - X(4,k)/X(1,k) );

printf("%f %f %f", x0, y0, r );