应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 【求助】计算N幅图像中圆的中心点坐标【已解决】
61/1返回列表
查看: 1693  |  回复: 5
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前主题已经存档。
本帖产生 1 个 仿真EPI ，点击这里进行查看

飞扬2282

荣誉版主 (著名写手)

[交流] 【求助】计算N幅图像中圆的中心点坐标【已解决】 已有2人参与

已完成！！

[ Last edited by 飞扬2282 on 2010-4-15 at 10:46 ]
回复此楼

» 猜你喜欢
1楼 2010-04-12 11:17:29
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

anyuezhiji

银虫 (正式写手)

星空行者
★ ★
飞扬2282(金币+20):我要读取的是SPE格式的文件，而且是N幅，在一个文件夹里，不知你有什么好办法，我再追加金币！ 2010-04-12 20:09
adu886886(金币+2):辛苦了 2010-04-13 08:26
下面的代码应该可以解决你的问题
引用回帖:
function Find_Center_of_circles


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

%边缘提取
I=imread('picture.TIF');
I=rgb2gray(I);
I=imadjust(I);
Image=I;
figure, imshow(I);

[m,n]=size(Image);
Model=zeros(m,n);
for i=1:m;
for j=1:n;
Model(i,j)=Image(i,j);
end
end

%观察目标的灰度范围
%如果imshow白色的点显示的图不是要查找的两个圆的边缘，调整Model<210
%L=bwlabel(mask)标记对象
mask=Model<210;
figure,imshow(mask);
L=bwlabel(mask);

%show出圆的边缘
for i=1:2
mask1=mask;
mask1(find(L~=i))=0;
BW = edge(uint8(mask1),'sobel');
figure,imshow(BW);
[r,c]=find(BW~=0);
eval(['edge',num2str(i),'=[r c];']);
end
save('edge.mat','edge1', 'edge2');

for i=1:2
eval(['edgei=edge',num2str(i),';']);      
edgei=pixels2Coordinates(edgei);
[r,center]=findcenter(edgei,3)
figure
hold on
plot(edgei(:,1),edgei(:,2),'.');
plot(center(1),center(2),'r.');
end

function out=pixels2Coordinates(Pixs)
%在像素图片上建立坐标系
%Lx,Ly 长、宽方向像素
%x0 y0 坐标系原点 这里取图像中点
Lx=1024;
Ly=768;
Tx=3.78;
Ty=3.78;
x0=(1+Lx)/2;
y0=(1+Ly)/2;
out=((Pixs(:,2)-y0))/Tx;
out=[out,(Ly-Pixs(:,1)-x0)/Ty];

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 );
       %'PopulationSize',50,'TolFun',1e-7,'TolCon',1e-7,'Display','diagnose');
    [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;

[ Last edited by anyuezhiji on 2010-4-12 at 21:42 ]
赞一下(2人) 回复此楼
暗月下没有留下风的痕迹，但它已经寂然飘逝。。By&amp;amp;lt;暗月之寂&amp;amp;gt;:tiger38:
2楼2010-04-12 18:17:40
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

anyuezhiji

银虫 (正式写手)

星空行者
★ ★
adu886886(金币+2):辛苦了 2010-04-13 08:27
spe格式没碰到过
能存为JPG格式么

JPG格式是没问题的

文件多的话写个循环就可以了
引用回帖:
%%PSL@CSU
%%QQ:547423688
%%Email:anyuezhiji@qq.com
%%Edit @ 2010.4.12

files = dir('*.jpg ')；

k=length(files );

R=zeros(2*k,1);
C=zeros(2*k,2);

for i=1:k

I=imread([files(k).name]);

%读入图片

%开始图片处理
[r,c]=Find_Center_of_Circle(I)
R([2*i-1 2*i])=r;
C([2*i-1 2*i], : )=c;

end

%图片处理代码
引用回帖:
function [R,C]=Find_Center_of_Circle(I)

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

I=rgb2gray(I);
Image=I;
figure, imshow(I);

[py,px]=size(Image);
Model=zeros(py,px);
for i=1:py;
for j=1:px;
Model(i,j)=Image(i,j);
end
end

%观察目标的灰度范围
%如果imshow白色的点显示的图不是要查找的两个圆的边缘，调整Model<50
%L=bwlabel(mask)标记对象
mask=Model<50;
L=bwlabel(mask);

%show出圆的边缘
r=[];c=[];
for i=1:2
mask1=mask;
mask1(find(L~=i))=0;
BW = edge(uint8(mask1),'sobel');
[r,c]=find(BW~=0);
eval(['BW',num2str(i),'=BW;']);
eval(['edge',num2str(i),'=[r c];']);
end
Model=(BW1|BW2)*1;



for i=1:2
eval(['edgei=edge',num2str(i),';']);  
[r,center]=findcenter(edgei,2)
R=[R;r];
C=[C;center]
Model(round(center(1)),round(center(2)))=1;
end
figure,imshow(Model);



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;

读入的图片


处理后输出的圆边缘和圆心


[ Last edited by anyuezhiji on 2010-4-14 at 01:04 ]
赞一下(1人) 回复此楼
暗月下没有留下风的痕迹，但它已经寂然飘逝。。By&amp;amp;lt;暗月之寂&amp;amp;gt;:tiger38:
3楼2010-04-12 21:16:13
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

飞扬2282

荣誉版主 (著名写手)
up
4楼2010-04-13 15:21:41
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

anyuezhiji

银虫 (正式写手)

星空行者
★ ★
飞扬2282(金币+100):谢谢帮助，太感谢你了。 2010-04-13 23:54
adu886886(金币+2):非常感谢指导，欢迎常来仿真模拟板块 2010-04-14 08:11
引用回帖:
function Find_CirclesCenter

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

path='D:\Program Files\MATLAB71\work\CurrentWork\test\';
clear CResult;
n=0;
for i=1:720
filename=['6-',num2str(i),'.spe'];  
if exist([path,filename],'file')
   disp(['正在处理',filename]);
   n=n+1;
   img=ReadSpe([path,filename]);
   mask=img<0.22;
   %figure;imshow(mask);
   CResult(n).fname=['6-',num2str(i)];
   CResult(n).R=[];
   CResult(n).C=[];
   [L,m]=bwlabel(mask);
   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
   CResult(n).num=m;
   [tc,p]=sort(CResult(n).C(:,2));
   CResult(n).C=CResult(n).C(p,: );
end
end

maxcm=0;
fidout=fopen([path,'CResult.txt'],'w');   
for i=1:n
fprintf(fidout,'%-8s',CResult(i).fname);  
cc=CResult(i).C;
[cm,cn]=size(cc);
if cm>maxcm
maxcm=cm;
end
for j=1:cm
fprintf(fidout,'\t\t圆心%d：(%5d，%5d)',j,cc(j,1),cc(j,2));
end
fprintf(fidout,'\n');  
end
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

fclose(fidout);

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文件
%读取.spe文件并显示图像如下面那样调用函数即可
%img=ReadSpe(filename)
%figure;imshow(img);
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);

自动读取spe文件，处理图片，并计算结果

输出文件内容如下：
引用回帖:
6-60                    圆心1：(   66，  144)                圆心2：( 1887， 1955)
6-120                   圆心1：(  627，  119)                圆心2：( 1405， 1937)


%圆心1坐标
x1=[   66  627]';
y1=[  144  119]';
%圆心2坐标
x2=[ 1887 1405]';
y2=[ 1955 1937]';

[ Last edited by anyuezhiji on 2010-4-14 at 01:19 ]
赞一下(1人) 回复此楼
暗月下没有留下风的痕迹，但它已经寂然飘逝。。By&amp;amp;lt;暗月之寂&amp;amp;gt;:tiger38:
5楼2010-04-13 22:27:05
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

anyuezhiji

银虫 (正式写手)

星空行者
★ ★ ★
小木虫(金币+0.5):给个红包，谢谢回帖交流
adu886886(金币+2, 仿真EPI+1):谢谢认真指导 2010-04-15 08:27
最终修改的程序 稍有变化 比如滤掉杂点的功能

有兴趣的虫友话可以把它和前面的代码比较
找找有什么不同
CODE:
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);

[ Last edited by anyuezhiji on 2010-4-14 at 21:07 ]
赞一下(3人) 回复此楼
暗月下没有留下风的痕迹，但它已经寂然飘逝。。By&amp;amp;lt;暗月之寂&amp;amp;gt;:tiger38:
6楼2010-04-14 21:03:47
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 飞扬2282 的主题更新
61/1返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
信息提示
关闭
请填处理意见
关闭