function Find_Center_of_circles2
%%Find Center of circles from picture
%%PSL@CSU
%%QQ:547423688
%%Email:anyuezhiji@qq.com
%%Edit @ 2010.4.23
path='F:\Program Files\MATLAB\work\CurrentWork\test\';
filename='zhanli 512 512 90-float.dat';
outfile='CResult.txt';
clear CResult;
fid = fopen([path filename],'rb','ieee-le');
n=0;
maxcm=0;
flag=1;
for i=1:90
   fname=['Rusult',num2str(i)];
   disp(['正在处理第',num2str(i),'张图。。。']);
   n=n+1;
   img = fread(fid,[512, 512],'float32');
   %figure
   %image(img);
   %imshow(img, []);
       minp=max(max(img));
       mask=img>minp/3;
       mask = bwmorph(mask,'clean');
       mask = bwmorph(mask,'majority',inf);
   %figure;imshow(mask);
   CResult(n).fname=fname;
   CResult(n).R=[];
   CResult(n).C=[];
   [L,m]=bwlabel(mask);
   %m
   if flag
   for i=1:m
   mask1=mask;
   mask1(find(L~=i))=0;
   BW = edge(uint8(mask1),'sobel');
   %figure,imshow(BW);
   [r,c]=find(BW~=0);
   %eval(['BW',num2str(i),'=BW;']);
   %eval(['edge',num2str(i),'=[r c];']);
   CircleEdge=[r,c];
   [r,center]=findcenter(CircleEdge,2);
   CResult(n).R=[CResult(n).R;round(r)];
   CResult(n).C=[CResult(n).C;round(center([2 1]))];
   end
   if m>maxcm
   maxcm=m;
   end
   CResult(n).num=m;
   [tc,p]=sort(CResult(n).C(:,2));
   CResult(n).C=CResult(n).C(p,;
   fidout=fopen([path,outfile],'a+');
   fprintf(fidout,'%-8s',CResult(n).fname);
   cc=CResult(n).C;
   for j=1:m
   fprintf(fidout,'\t\t圆心%d：(%5d，%5d)',j,cc(j,1),cc(j,2));
   end
   fprintf(fidout,'\n');
   fclose(fidout);
   end%end if flag
end
status = fclose(fid);
if flag
fidout=fopen([path,outfile],'a+');
fprintf(fidout,'\n');
fprintf(fidout,'\n');
for j=1:maxcm
fprintf(fidout,'%%圆心%d坐标\n',j);
fprintf(fidout,'x%d=[',j);
for i=1:n
if j'%5d',CResult(i).C(j,1));
else
fprintf(fidout,'  nan');
end
end
fprintf(fidout,']'';\n');
fprintf(fidout,'y%d=[',j);
for i=1:n
if j'%5d',CResult(i).C(j,2));
else
fprintf(fidout,'  nan');
end
end
fprintf(fidout,']'';\n');
end
fprintf(fidout,'\n');
fprintf(fidout,'\n');
fclose(fidout);
end%end if flag
function [R,C]=findcenter(edge,mode)
%mode1 相加取中点
%mode2 最小二乘拟合圆x^2+y^2+a(1)*x+a(2)*y+a(3)=0
%mode3 用遗传算法，找到一点和一个半径R，各点到这点的距离与R之差的平方和最小
[m,n]= size(edge);
switch mode
    case 1
    C=sum(edge)/m;
    R=sum(sqrt((edge(:,1)-C(1)).^2+(edge(:,2)-C(2)).^2))/m;
    case 2
    A=[edge(:,1) edge(:,2) ones(m,1)];
    B=-[edge(:,1).*edge(:,1)+edge(:,2).*edge(:,2)];
    a=A\B;
    C=-.5*a([1 2])';
    R = sqrt((a(1)^2+a(2)^2)/4-a(3));
    case 3
    save('edge_m','edge','m');
    %[x, fval, reason]=ga(fitnessfcn,nvars,A,b,Aeq,beq,LB,UB,nonlcon,options)
    C=sum(edge)/m;
    R=sum(sqrt((edge(:,1)-C(1)).^2+(edge(:,2)-C(2)).^2))/m;
    options = gaoptimset('Generations',200,'InitialPopulation',[R,C],...
        'PopulationSize',50,'TolFun',1e-7,'TolCon',1e-7,'MutationFcn',@mutationadaptfeasible );
    [sol, fval, exitflag] = ga(@f,3,[],[],[],[],[0,min(edge)],[2*R,max(edge)],[],options);
    R=sol(1);
    C=sol([2,3]);
    fval
end
%适应度函数的matlab代码
function [eval]=f(sol)
load edge_m
R=sol(1);
C=sol([2,3]);
eval=sum((sqrt((edge(:,1)-C(1)).^2+(edge(:,2)-C(2)).^2)-R).^2)/m;