function [a,b,r]=cire(X,Y);
[w,len]=size(X);
a=0;
b=0;
r=0;
S40=sum(X.*X.*X.*X);
S22=sum(X.*X.*Y.*Y);
S04=sum(Y.*Y.*Y.*Y);
S30=sum(X.*X.*X);
S21=sum(X.*X.*Y);
S12=sum(X.*Y.*Y);
S03=sum(Y.*Y.*Y);
S20=sum(X.*X);
S11=sum(X.*Y);
S02=sum(Y.*Y);
S10=sum(X);
S01=sum(Y);
p=len*(S30+S12)-S10*(S20+S02);
q=len*(S21+S03)-S01*(S20+S02);

nom= (S10*S10-S20*len)*(S01*S01-S02*len)-(S10*S01-S11*len)^2;
nom
if (nom > 1)
	a=q*(S10*S01-S11*len) - p*(S01*S01-S02*len);
	a=a/(2*nom);

	b=p*(S10*S01-S11*len) - q*(S10*S10-S20*len);
	b=b/(2*nom);

	c=(S20+S02)-2*S10*a-2*S01*b;
	c=c/len;
	r=sqrt(c+a*a+b*b);
error= S40 +2*S22 +S04 -4*a*S30-4*b*S21-4*a*S12-4*b*S03;
error=error + (4*a*a-2*c)*S20 + 8*a*b*S11 + (4*b*b-2*c)*S02;
error=error + 4*a*c*S10 + 4*b*c*S01 + c*c*len;

error =sqrt(error/len);
error=sqrt(error);
error

end;
end;