function reconstruction
clear;
close all;
N=256;%N＊N大小图像
delta=3; %相当于Φ
a = 100;
b = 50;
d = 1;
%生成椭圆
for m =1:180/delta;
L=a^2*cos(m*delta)^2+b^2*sin(m*delta)^2;
for t=-N/2+1:N/2;
if (t^2 <= L)
proj(t+N/2,m)=2*a*b*sqrt(L-t^2)/L;
else
proj(t+N/2,m)=0;
end
end
end
%%%%各种滤波函数
for t = -N+1:N-1;
h1(t+N)=-2/(pi^2*d^2*(4*t^2-1)); % S-L
end;
for t=-N+1:N-1;
if t==0
h2(t+N)=1/(4*d^2); % R-L
elseif mod(t,2) == 1
h2(t+N)=-1/(t^2*pi^2*d^2);
else
h2(t+N)=0;
end;
end;
% Lewitt
esp=0.5;
for t=-N+1:N-1;
if t==0
h3(t+N)=(1-2*esp/3)/(4*d^2);
elseif mod(t,2)==1
h3(t+N)=-(1-esp)/(t^2*pi^2*d^2);
else
h3(t+N)=-esp/(t^2*pi^2*d^2);
end;
end;
%存储各种反投影结果
rProj1=zeros(N,N);
rProj2=zeros(N,N);
rProj3=zeros(N,N);
rProj4=zeros(N,N);
nProj=[];
%%调用重建函数
produce('ori');
produce('sl');
produce('rl');
produce('lew');
function produce(filter)
for m=1:180/delta;
creat_nProj(m);%补零，用于卷积
c_h1=conv(nProj,h1);%求卷积
c_h2=conv(nProj,h2);
c_h3=conv(nProj,h3);
c=0.5*(N-1)*(1-cos(m*delta)-sin(m*delta));
for i=1:N;
for j=1:N;
L(i,j)=c+(i-1)*cos(m*delta)+(j-1)*sin(m*delta); %内插
n = fix(L(i,j));
cL=L(i,j)-n;
if strcmp(filter,'ori')
if (n>0)&(n<255)
rProj1(i,j)=rProj1(i,j)+(1-cL)*nProj(n+N)+cL*nProj(n+1+N);
elseif n==255
rProj1(i,j) =rProj1(i,j)+nProj(n+N);
elseif n==0
rProj1(i,j) =rProj1(i,j)+nProj(n+N+1);
end
elseif strcmp(filter,'sl')
if (n>0)&(n<255)
rProj2(i,j)=rProj2(i,j)+(1 - cL)*c_h1(n+N+N-1)+cL*c_h1(n+N+N);
elseif n==255
rProj2(i,j) =rProj2(i,j)+c_h1(n+N+N-1);
elseif n==0
rProj2(i,j)=rProj2(i,j)+c_h1(n+N+N);
end
elseif strcmp(filter,'rl')
if (n>0)&(n<255)
rProj3(i,j)=rProj3(i,j)+(1 - cL)*c_h2(n+N+N-1)+cL*c_h2(n+N+N);
elseif n==255
rProj3(i,j)=rProj3(i,j)+c_h2(n+N+N-1);
elseif n==0
rProj3(i,j) =rProj3(i,j)+c_h2(n+N+N);
end
elseif strcmp(filter,'lew')
if (n>0)&(n<255)
rProj4(i,j)=rProj4(i,j)+(1 - cL)*c_h3(n+N+N-1)+cL*c_h3(n+N+N);
elseif n==255
rProj4(i,j) =rProj4(i,j)+c_h3(n+N+N-1);
elseif n==0
rProj4(i,j) =rProj4(i,j)+c_h3(n+N+N);
end
end
end
end
end
end
function creat_nProj(m)
for t=1:((N-1) * 3 + 1);
if (t<N)
nProj(t)=(proj(1,m)+proj(2,m))/2;
elseif (t >= 2*N)
nProj(t)=(proj(N-1,m)+proj(N,m))/2;
else
nProj(t)=proj(t-255,m);
end
end
end
%%显示输出结果
subplot(2,2,1);
imshow(rProj1, [min(min(rProj1)), max(max(rProj1))]);
title('直接反投影');
subplot(2,2,2);
imshow(rProj2, [min(min(rProj2)), max(max(rProj2))]);
title('S-L滤波');
subplot(2,2,3);
imshow(rProj3, [min(min(rProj3)), max(max(rProj3))]);
title('R-L滤波');
subplot(2,2,4);
imshow(rProj4, [min(min(rProj4)), max(max(rProj4))]);
title('Lewitt滤波');
end
