function Find_CirclesCenter

%%Find Center of circles from picture
%%PSL@CSU
%%QQ:547423688
%%Email:anyuezhiji@qq.com
%%Edit @ 2010.4.14

path='D:\Program Files\MATLAB71\work\CurrentWork\test\';
outfile='CResult.txt';
clear CResult;
n=0;
maxcm=0;
flag=1;
for i=1:720
fname=['8-',num2str(i)];
filename=[fname,'.spe'];  
if exist([path,filename],'file')
   disp(['正在处理',filename]);
   n=n+1;
   [img,dtype]=ReadSpe([path,filename]);
       %figure;imshow(img);
       minp=min(min(img));
       mask=img<3*minp;
       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
end

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<=CResult(i).num
fprintf(fidout,'%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<=CResult(i).num
fprintf(fidout,'%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;

%读取.spe文件
function [data,data_type] = ReadSpe(filename,startX,endX,startY,endY)
   fid=fopen(filename,'r','ieee-le');

   [header,count]=fread(fid,4100,'uint8');
   dataType=header(109);
   xdim=header(43) + header(44)*2^8;
   ydim=header(657) + header(658)*2^8;
   frames=header(1447) + header(1448)*2^8 +header(1449)*2^16 + header(1450)

*2^24;

    if ~exist('startX','var')
          startX = 1;
    end
    if ~exist('endX','var')
          endX=xdim;
    end
    if ~exist('startY','var')
          startY=1;
    end
    if ~exist('endY','var')
          endY=ydim;
    end
   strDataType={'float32','int32','int16','uint16'};
   data_type=strDataType{dataType+1};
   if (dataType>3)
      fclose(fid);
      fprintf(strcat('Unknown WinSpec data type: ',num2str(dataType)));
      clear data;
      data=[];
      return;
   end;
   data=zeros([xdim,endY-startY+1,frames]);
   dataSize=max((dataType<2)+1)*2;
   for z=1:frames
         fseek(fid,4100+dataSize*(xdim*((startY-1)+ydim*(z-1))),'bof');
         [data(1:xdim,1:(endY-startY+1),z),count]=fread(fid,[xdim (endY-startY+1)],strcat(strDataType{dataType+1},'=>',strDataType{dataType+1}));
   end;
   data=permute(data,[2 1 3]);
   data=data(:,startX:endX,:);
   fclose(fid);