第十单元 第3课 实验 2017年全国大学生 数学建模竞赛A题

问题一的代码： function[x0,y0,sepdist,start_angle,direction180,gain] = problem1(input_args) %模板:template template=xlsread('a题附件.xls','附件1'); figure;imshow(template,[]); title('模板'); %1.计算公切线的斜率 [b,fval]=fmincon(@(b)myopt(b),b0,[],[],[],[],40,43,@(b)nonlcon(b)); k1=(-90*b+sqrt(8100*b^2-8036*(b^2-16)))/4018 k2=(-90*b-sqrt(8100*b^2-8036*(b^2-16)))/4018; theta1=atan(k1)*180/pi; theta2=atan(k2)*180/pi; x1=real(((90-2*k1*b)+sqrt((90-2*k1*b)^2-4*(1+k1^2)*(2009+b^2)))/(2*(1+k1^2))); y1=k1*x1+b; x2=real(((90-2*k2*b)-sqrt((90-2*k2*b)^2-4*(1+k2^2)*(2009+b^2)))/(2*(1+k2^2))); y2=k2*x2+b;

%从上到下4条公切线的端点与斜率,截距 xykb=[[0,x1;b,y1;k1,b],[0,x2;b,y2;k2,b],[0,x2;-b,-y2;-k2,-b],[0,x1;-b,-y1;-k1,-b]]; %ct_radon为该CT系统对template_img进行radon变换之后的投影图 ct_radon=xlsread('a题附件.xls','附件2'); figure;imshow(ct_radon,[]); title('旋转中心有偏移的模板的CT图'); %画模板的CT图中各个投影方向的投影 figure; interval=zeros(512,180); for i=1:180   plot(ct_radon(:,i));   grid on;   hold on;   pos=find(ct_radon(:,i)~=0);   for j=1:length(pos)-1       interval(pos(j),i)=pos(j+1)-pos(j);   end  end

 xlabel('投影地址');  ylabel('投影值');  title('旋转中心有偏移的模板的CT图中各个方向的投影');  saveas(gcf,'有偏移的模板的CT图中各个方向的投影','jpg');  %2.计算公切线投影方向在ct_radon中对应的投影线号  k=0;  for i=1:180      ifsum(interval(:,i))==length(find(interval(:,i)~=0))          k=k+1;          cprojection(k)=i;      end  end  radioe=cprojection(end);  radios=cprojection(1);  figure;  plot([1:512],ct_radon(:,radios),'r',[1:512],ct_radon(:,radioe),'b');  grid on;

xlabel('投影地址');  ylabel('投影值');  legend('第1条内公切线方向的投影','第2条内公切线方向的投影');  title('旋转中心有偏移的模板的CT图中椭圆与圆的两条内公切线方向的投影');  saveas(gcf,'有偏移的模板的CT图中椭圆与圆的两条内公切线方向的投影','jpg');  %3.CT系统相邻投影方向的角度差  delta_angle=2*abs(theta2)/(radioe-radios);  %4.CT系统的起始投影方向角度  start_angle=theta2-(radios-1)*delta_angle;  %第1条投影线方向对应的投影检测器与水平方向的夹角  horizen_start_angle=90+start_angle;  %5.CT系统使用的X射线的180个投影方向

 direction180=start_angle:delta_angle:(start_angle+179*delta_angle);  for i=1:18     direction180paper(2*i-1,:)=(i-1)*10+1:i*10;    direction180paper(2*i,:)=direction180((i-1)*10+1:i*10);  end  xlswrite('CT系统使用的X射线的180个方向paper.xls',direction180paper);  %6.确定探测器单元间距sepdist  projnum=sum(interval(:,find(abs(direction180)<delta_angle)));  sepdist=80/projnum;  %探测器总长度  probe_length=(512-1)*sepdist;  xlswrite('CT系统的探测器单元间距与总长.xls',[sepdist,probe_length]);  %7.模板的理想CT图  theta=0:179;  N=512;

[template_radon,xp] = radon(template,theta,N);  figure,imshow(template_radon,[]);  title('旋转中心没有偏移的模板的CT图');  imwrite(uint8(template_radon),'没有偏移的模板的CT图.jpg','jpg');  xlswrite('没有偏移的模板的CT数据',template_radon);  figure;  plot([1:N],template_radon(:,1),'b',[1:N],template_radon(:,91),'r');  grid on;  xlabel('投影地址');  ylabel('投影值');  legend('无偏移0度投影曲线','无偏移的90度投影曲线');  title('旋转中心无偏移时0度与90度的投影曲线');  axis([0 512 0 max(max(template_radon))+10]) ;  saveas(gcf,'无偏移时0度与90度的投影曲线','jpg'); %8.CT系统旋转中心在正方形托盘中的位置(x0,y0)

template_radon90=template_radon(:,1); ct_col90=find(abs(direction180-90)<delta_angle/2); ct_radon90=flipud(ct_radon(:,ct_col90)); ct_col0=find(abs(direction180)<delta_angle); ct_radon0=ct_radon(:,ct_col0); figure; plot([1:N],template_radon0,'b',[1:N],ct_radon0,'r'); grid on; xlabel('投影地址'); ylabel('投影值'); legend(sprintf('无偏移的%d度水平方向投影',0),sprintf('有偏移的%d度水平方向投影',0)); title(sprintf('旋转中心有与无偏移的%d度水平方向投影',0)); saveas(gcf,'旋转中心有与无偏移的水平方向的投影','jpg'); plot([1:N],template_radon90,'b',[1:N],ct_radon90,'r');

xlabel('投影地址'); ylabel('投影值'); legend(sprintf('无偏移的%d度垂直方向投影',90),sprintf('有偏移的%d度垂直方向投影',90)); title(sprintf('旋转中心有与无偏移的%d度垂直方向的投影',90)); saveas(gcf,'旋转中心有与无偏移的垂直方向的投影','jpg'); wtemplate90=sum(template_radon90); wtemplatex=[1:N]*template_radon90; wct90=sum(ct_radon90); wctx=[1:N]*ct_radon90; x0=(wtemplatex/wtemplate90-wctx/wct90)*sepdist; wtemplate0=sum(template_radon0); wtemplatey=[1:N]*template_radon0; wct0=sum(ct_radon0); wcty=[1:N]*ct_radon0; y0=(wtemplatey/wtemplate0-wcty/wct0)*sepdist; p1=find(template_radon0~=0); p2=find(ct_radon0~=0);

y0=(p1(end)-p2(end)+p1(1)-p2(1))*sepdist/2; xlswrite('CT系统旋转中心在正方形托盘中的位置.xls',[x0,y0]); function fval=myopt(b)    k=(-90*b-sqrt(8100*b^2-8036*(b^2-16)))/4018;     fval=324*k^2*b^2-36*(64+9*k^2)*(b^2-1600); end                           function [c,ceq] =nonlcon(b)    delta1=324*k^2*b^2-36*(64+9*k^2)*(b^2-1600);     c=-delta1;     ceq = []; end functionF=myfun(gain,rxv,y)  F=sum((rxv-gain*y).^2);

问题二的代码： function[reconstruct_attach3_256,xy0,semiaxis,rotate_angle,ellipse_absorb6,attach4_absorb10]=problem2(input_args) %1.比较模板的拉东变换图与附件2中给出的拉东变换投影线的对应关系 template_radon0=template_radon(:,1); template_radon90=template_radon(:,91); attach2_radon_col180=find(abs(direction180-180)<delta_angle/2); attach2_radon180=attach2_radon(:,attach2_radon_col180); attach2_radon_col90=find(abs(direction180-90)<delta_angle/2); attach2_radon90=attach2_radon(:,attach2_radon_col90); %2.有无偏移CT投影配准 alpha=atan(-x0/y0)*180/pi; reconstruct_radon=zeros(size(attach2_radon)); figure; for i=1:180     theta=(i-1)*delta_angle;    attach2_radon_anglei=find(abs(direction180-theta)<delta_angle/2);

if isempty(attach2_radon_anglei)           attach2_radon_anglei=find(abs(direction180-180-theta)<delta_angle/2);     end     if theta<start_angle-delta_angle        reconstruct_radon(:,i)=flipud(attach2_radon(:,attach2_radon_anglei));     else        reconstruct_radon(:,i)=attach2_radon(:,attach2_radon_anglei);     if abs(theta)<delta_angle         d=abs(x0);     elseif  abs(theta-180)<delta_angle          d=y0;         k=tan(theta*pi/180+pi/2);         d=abs(-k*x0+y0)/sqrt(1+k^2);

address_num=round(d/sepdist);     if address_num>0           if theta<=alpha |  theta>180+alpha                reconstruct_radon(:,i)=[reconstruct_radon(address_num+1:end,i);zeros(address_num,1)];           else               reconstruct_radon(:,i)=[zeros(address_num,1);reconstruct_radon(1:end-address_num,i)];           end     end    plot([1:N],template_radon(:,i),'r',[1:N],reconstruct_radon(:,i),'b');     hold on; end  A=reshape(template_radon,512*180,1);  B=reshape(reconstruct_radon,512*180,1);  stdtr=std(A-B);  legend(sprintf('无偏移的投影'),sprintf('有偏移的投影'));

grid on;  xlabel('投影地址');  ylabel('投影值');  title(sprintf('重新配准之后的旋转中心无偏移与有偏移的模板的各个方向的投影'));  saveas(gcf,sprintf('配准后的旋转中心无偏移与有偏移的模板的各个方向的投影'),'jpg');  %3.根据投影配准之后的reconstruct_radon利用iradon反变换重构模板(附件2)原图像  reconstruct_template362=iradon(reconstruct_radon,theta180);  reconstruct_template256=imresize(reconstruct_template362,[256 256]);  figure;imshow(reconstruct_template256,[]);title('重构模板的256X256图像');  imwrite(reconstruct_template256,'重构模板的256X256图像.jpg','jpg');  xlswrite('重构模板的256X256图像数据.xls',reconstruct_template256);

reconstruct_radon=zeros(size(attach3_radon));  for i=1:180     theta=(i-1)*delta_angle;    attach3_radon_anglei=find(abs(direction180-theta)<delta_angle/2);     if isempty(attach3_radon_anglei)          attach3_radon_anglei=find(abs(direction180-180-theta)<delta_angle/2);     end     if theta<start_angle-delta_angle        reconstruct_radon(:,i)=flipud(attach3_radon(:,attach3_radon_anglei));     else        reconstruct_radon(:,i)=attach3_radon(:,attach3_radon_anglei);     if abs(theta)<delta_angle         d=abs(x0);     elseif  abs(theta-180)<delta_angle          d=y0;

  k=tan(theta*pi/180+pi/2);         d=abs(-k*x0+y0)/sqrt(1+k^2);     end     address_num=round(d/sepdist);     if address_num>0           if theta<=alpha |  theta>180+alpha                reconstruct_radon(:,i)=[reconstruct_radon(address_num+1:end,i);zeros(address_num,1)];           elsereconstruct_radon(:,i)=[zeros(address_num,1);reconstruct_radon(1:end-address_num,i)];           end  end  %4.根据投影配准之后的reconstruct_radon利用iradon反变换重构附件3原图像  reconstruct_attach3_362=iradon(reconstruct_radon,theta180);  reconstruct_attach3_256=imresize(reconstruct_attach3_362,[256 256]);  [row,col]=find(reconstruct_attach3_256<0);

for i=1:length(row)      reconstruct_attach3_256(row(i),col(i))=0;  end  figure;imshow(reconstruct_attach3_256,[]);title('重构附件3的256X256图像');  imwrite(reconstruct_attach3_256,'重构附件3的256X256图像.jpg','jpg');  xlswrite('problem2.xls',reconstruct_attach3_256); %5.根据重构的附件3的原图像,确定各个椭圆的位置和形状  attach3_256bw_edge=read_sixellipse_edge; [xy0,semiaxis,rotate_angle,pixelrc0]=compute_sixellipse_parameters(attach3_256bw_edge); %6.重构的附件3的原图像6个椭圆区域的吸收率,附件4中10个给定位置处的吸收率 [ellipse_absorb6,attach4_absorb10]=compute_problem2_absorb(); function fval=myopt(b)     k=(-90*b-sqrt(8100*b^2-8036*(b^2-16)))/4018;     fval=324*k^2*b^2-36*(64+9*k^2)*(b^2-1600); end                         

function [c,ceq] =nonlcon(b)    k=(-90*b-sqrt(8100*b^2-8036*(b^2-16)))/4018;    delta1=324*k^2*b^2-36*(64+9*k^2)*(b^2-1600);     c=-delta1;     ceq = []; end functionF=myfun(gain,rxv,y)  F=sum((rxv-gain*y).^2); functionattach3_256bw_edge=read_sixellipse_edge     for i=1:6       attach3_gray_edge=imread(sprintf('重构附件3的256X256图像中第%d个图形的2值边沿图像.jpg',i));       attach3_256bw_edge(:,:,i)=zeros(size(attach3_gray_edge));       [I,J]=find(attach3_gray_edge>200);       for j=1:length(I)           attach3_256bw_edge(I(j),J(j),i)=1;       end  figure;imshow(attach3_256bw_edge(:,:,i));title(sprintf('重构附件3的256X256图像中第%d个图形的2值边沿图像',i));   end