2017年数学建模A组练习代码

clear,clc
%%表5处理图形
D=xlsread('D:\sparename\Documents\2017\A\A.xls',5)
% imagesc(D)
R=iradon(D,[0:179]+30,512)
imagesc(R)

clear,clc
%%表1处理图形
A=xlsread('D:\sparename\Documents\2017\A\A.xls',1)
imshow(A)
%%附件2
figure
B=xlsread('D:\sparename\Documents\2017\A\A.xls',2)
imagesc(B)

题目1

%% 求探测器距离d
B=xlsread('A题附件.xls',2);
imagesc(B)
W=sum(B>0);
d=80/max(W)
%% 求旋转中心偏移量
imin=min(W),jmin=find(W(:)==imin);%得投影最小长度imin探测器
i=jmin(1);
x=find(B(:,i)>0)%此时对应角度的所有接受的探测器
% 512/2得到近似中心探测器编号，而(max(x)+min(x))/2)得到椭圆中心对应的探测器编号。编号之差乘d得到距离。
nx=(max(x)+169)
idx=(512/2-(max(x)+169)/2)*d;
jmax=150-90;
y=find(B(:,jmax)>0);%此时对应角度的所有接受的探测器
idy=(512/2-(max(y)+min(y))/2)*d;
%%求偏移角度
alpha=180-jmin
j=jmin(2)
x=find(B(:,j)>0);%此时对应角度的所有接受的探测器
jdx=(512/2-(max(x)+169)/2)*d;
jmax=150-90
y=find(B(:,jmax)>0)%此时对应角度的所有接受的探测器
jdy=(512/2-(max(y)+min(y))/2)*d;
dx=-(idx+jdx)/2
dy=(idy+jdy)/2

%%求偏移角度
alpha=180-jmin
%% 验证
B=[zeros(200,180);B;zeros(200,180)]

theta=[0:1:179]+30
% 基于R-L滤波器的滤波反投影法
R= iradon(B,theta,'Hann')
figure
imshow(R)

题目2

clear,clc
%%表三处理
C=xlsread('D:\sparename\Documents\2017\A\A.xls',3)
C=[zeros(200,180);C;zeros(200,180)];
R=iradon(C,[0:179]+30,'Hann');
imagesc(R)
d=0.2768
dx=-9.2042
dy=5.8131
x=ceil(-dx/d)+size(R,1)/2+[ceil(-50/d):ceil(50/d)];
y=ceil(dy/d)+size(R,1)/2+[ceil(-50/d):ceil(50/d)];
% 调整大小
R=imrescale(R(y,x),[256,256])
T=flipud(R)
D=xlsread('D:\sparename\Documents\2017\A\A.xls',4)
d=100/256
c=ceil(D(:,1)/d)
r=ceil(D(:,2)/d)
ix=sub2ind(size(T),r,c);
rate=[D T(ix)];
imagesc([0,100],[0,100],T);
colorbar;
hold on
plot(D(T(ix)==0,1),D(T(ix)==0,2),'wx');
plot(D(T(ix)~=0,1),D(T(ix)~=0,2),'r+');
set(gca,'ydir','normal')
% 
% figure
% imagesc([0,179]+30,[1,512],R)
%函数
function new=imrescale(old,newsize)
scale=numel(old)/prod(newsize);
new=imresize(old,newsize)*scale;
new(new<1e-1)=0;
end