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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 【求助】再次求助
31/1返回列表
查看: 669  |  回复: 2
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
本帖产生 1 个 仿真EPI ，点击这里进行查看

飞扬2282

荣誉版主 (著名写手)

[交流] 【求助】再次求助 已有1人参与

http://muchong.com/bbs/viewthread.php?tid=1952310

6楼的程序。

现在要求略有变化，数据存储在C:\TDdownload\data.dat 文件中，数据为512x512x
90,即90幅数据，请问如何读取？
应该在6楼的程序上修改即可？
回复此楼

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

anyuezhiji

银虫 (正式写手)

星空行者
★ ★ ★
飞扬2282(金币+30):高手啊，再次感谢，牛人！ 2010-04-23 22:29
adu886886(金币+3, 仿真EPI+1):AN版造诣很好啊 2010-04-24 11:11
代码如下：

可以和http://muchong.com/bbs/viewthread.php?tid=1952310
6楼做比较
引用回帖:
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;
赞一下(1人) 回复此楼
暗月下没有留下风的痕迹，但它已经寂然飘逝。。By<暗月之寂>:tiger38:
2楼2010-04-23 21:53:49
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

anyuezhiji

银虫 (正式写手)

星空行者
★ ★
小木虫(金币+0.5):给个红包，谢谢回帖交流
飞扬2282(金币+20):追加金币，感谢帮助 2010-04-24 07:56
adu886886(金币+1):尽心尽力。赞一个 2010-04-24 11:11
修改版 读数据 计算坐标
包含生成jpg gif avi 等附加功能
引用回帖:
function Find_Center_of_circles2
%%Find Center of circles from picture
%%PSL@CSU
%%QQ:547423688
%%Email:anyuezhiji@qq.com
%%Edit @ 2010.4.23
%%Modify @ 2010.4.24
clear all
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;%是否执行计算
savemov=0;%是否生成avi文件，1为是，0为否
savejpg=0;%是否生成jpg文件，1为是，0为否
savegif=1;%是否生成gif文件，1为是，0为否
if(savemov==1)
mov=avifile([path,'CResult.avi']);
end
%[header,count]=fread(fid,4100,'uint8');
for i=1:90
   fname=['Rusult',num2str(i)];
   disp(['正在处理第',num2str(i),'张图。。。']);
   n=n+1;
      img = fread(fid,[512, 512],'float32');
      img=img';
      img=img(512:-1:1,;
      %img=img(:,512:-1:1);
   fig=figure(1);
   %image(img);
   imshow(img, []);
   %movie
   pause(0.01);
   %drawnow;%刷新
   if(savegif==1)
   %f=getframe(fig);
   %imind=frame2im(f);
   %[imind,cm] = rgb2ind(imind,256);
   if n==1
       %imwrite(imind,cm,[path 'CResult.gif'],'gif', 'Loopcount',inf,'DelayTime',0.1);
       imwrite(im2uint8(img),[path 'CResult.gif'],'gif', 'Loopcount',inf,'DelayTime',0.1);
   else
       %imwrite(imind,cm,[path 'CResult.gif'],'gif','WriteMode','append','DelayTime',0.1);
       imwrite(im2uint8(img),[path 'CResult.gif'],'gif','WriteMode','append','DelayTime',0.1);
   end
   end
   if(savejpg==1)
   print(1,'-djpeg',[path 'Result_' fname '.jpg']);
   end
   if(savemov==1)
   f=getframe(fig);
   mov=addframe(mov,f);
   end
       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
if(savemov==1)
mov=close(mov);
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;

运行结果
引用回帖:
正在处理第1张图。。。
正在处理第2张图。。。
正在处理第3张图。。。
。。。。。。

结果写入CResult.txt
引用回帖:
Rusult1                 圆心1：(  257，   99)                圆心2：(  257，  457)
Rusult1                 圆心1：(  257，   99)                圆心2：(  257，  457)
Rusult2                 圆心1：(  268，   99)                圆心2：(  243，  457)
。。。。。。

%圆心1坐标               
x1=[  257  268  279  290 。。。。。。]';
y1=[   99   99   99   99 。。。。。。]';
%圆心2坐标               
x2=[  257  243  229  215 。。。。。。]';
y2=[  457  457  457  457 。。。。。。]';

静态图
引用回帖:

。。。。。。

gif动画
引用回帖:

运动轨迹
引用回帖:

avi。。 略。。。。
赞一下(2人) 回复此楼
暗月下没有留下风的痕迹，但它已经寂然飘逝。。By<暗月之寂>:tiger38:
3楼2010-04-24 01:48:55
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 飞扬2282 的主题更新
31/1返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 面上再次挂了，太难了，躺也躺不了，倦也卷不过，小学校之殇！ +16 低垂的野花 2026-07-31 22/1100 2026-08-01 16:39 by Cathyau
[博后之家] 售SCI一区T0P文章，我:8.O55.1.O.54,科目全,可十急 +3 jzg3xk1po8 2026-07-31 9/450 2026-08-01 15:50 by y4XuShGyllvT
[公派出国] 售SCI-T0P文章，我:8O.5.5.1.O.54,科目齐全,可+急 +3 jzg3xk1po8 2026-07-31 6/300 2026-08-01 15:46 by y4XuShGyllvT
[考博] 售SCI-T0P文章，我:8O.5.5.1.O.54,科目齐全,可+急 +3 jzg3xk1po8 2026-07-31 11/550 2026-08-01 15:43 by y4XuShGyllvT
[基金申请] 售SCI一区文章，我:8.O.55.1.O.54,科目齐全,可伽急 +3 jzg3xk1po8 2026-07-31 11/550 2026-08-01 15:43 by y4XuShGyllvT
[教师之家] 售SCI一区T0P文章，我:8.O55.1.O.54,科目全,可十急 +3 jzg3xk1po8 2026-07-31 8/400 2026-08-01 15:42 by y4XuShGyllvT
[考研] 售SCI一区T0P文章，我:8.O.55.1.O.5.4,科目全,可+急 +3 jzg3xk1po8 2026-07-31 11/550 2026-08-01 15:40 by y4XuShGyllvT
[基金申请] 2026年国自然面上资助率 +15 布布和一二 2026-07-30 19/950 2026-08-01 15:18 by 40862026
[基金申请] 微信指数没变化，科研之友没阅读 +14 wangze12014 2026-07-28 18/900 2026-07-31 18:35 by lzchaihuo
[教师之家] 基础研究怎么拉横向，学校到款任务越来越多，难以完成 拉横向，都有哪些途径啊 +8 锦衣卫寒战 2026-07-28 8/400 2026-07-31 18:13 by beefly
[高分子] HXDI做水性聚氨酯乳液，是不是特别容易出渣 15+3 yuyusuv 2026-07-29 3/150 2026-07-31 09:21 by huizingga
[基金申请] 系统今天又提示维护了，估计离放榜不远了 +11 winnerche 2026-07-29 15/750 2026-07-30 23:08 by jnhyjjm
[基金申请] 你们的时间戳变了吗 +3 archvillain 2026-07-30 4/200 2026-07-30 18:53 by levinzhwen
[基金申请] 今年的WR进展到哪一步了？ +5 wsgjhwz 2026-07-25 9/450 2026-07-30 12:59 by 可淡不可忘
[有机交流] 产物和副产物价值比较 30+4 小汤02 2026-07-27 5/250 2026-07-30 12:44 by czyzsu
[基金申请] 时间戳他又来了 +14 晓晓爱翠翠 2026-07-26 17/850 2026-07-29 16:30 by 呵呵二哥
[有机交流] 这个自发加氧反应的机理是什么？ 50+4 YUAN2273 2026-07-25 7/350 2026-07-29 12:27 by xiaobenpu
[基金申请] 同事接到电话了，我却没有 +4 1234567wang 2026-07-27 4/200 2026-07-28 10:11 by GOODLUCKER
[基金申请] 准备明年的基金了 +5 Tide man 2026-07-27 5/250 2026-07-28 09:25 by newfuzzy1
[分析] 水蒸气蒸馏法乳化现象 +3 月光娇娇 2026-07-26 3/150 2026-07-26 23:25 by caiyun
信息提示
关闭
请填处理意见
关闭