clc
clear all
N = 256;
I = phantom(N);
imshow(I,[])
SO = 980;
SD = 1250;
objD = 200/11*16/N; %%单位mm 标志男人头的直径为200mm左右 占Sheep-Logan图像的11/16 可以计算头部每个体素的长为 200/11*16/512 = 0.5682
detD = objD/SO*SD;
D = round(SO/objD); %%光源到旋转中心的距离包含多少个体素
NumAngle = 360;
theta = (0:360/NumAngle:360-360/NumAngle)*pi/180; %%扫描角度，从X正轴逆时针旋转
cosd = cos(theta);
sind = sin(theta);
% M = ceil(N*sqrt(2)*SO/sqrt(SO^2-2*(N/2*objD)^2))+1; %%求探测的长度
objX = ((1:N)-(N+1)/2)*objD;
objY = ((N+1)/2-(1:N))*objD;
projdata = fanbeam(I,D,'FanSensorGeometry','line','FanSensorSpacing',1,'FanRotationIncrement',1);
PrjD = flipud(projdata');
M = size(PrjD,2);
% if len < M
% rz = M - len; % how many rows of zeros
% PrjD = [zeros(NumAngle,ceil(rz/2)),PrjD,zeros(NumAngle,floor(rz/2))];
% else
% M = len;
% end
%%投影数据
% for i = 1:NumAngle
% for j = 1:N
% for k = 1:N
% XX = objX(j)*cosd(i)+objY(k)*sind(i); %%逆时针选择delta后，体素在新坐标系中的坐标
% YY = -objX(j)*sind(i)+objY(k)*cosd(i);
% dd = (M+1)/2-YY/(SO-XX)*SD/detD; %%在探测器上的位置
% rd = floor(dd); %%双线性插值系数
% eff = dd-rd; %%双线性插值系数
% PrjD(i,rd) = PrjD(i,rd)+I(k,j)*(1-eff);
% PrjD(i,rd+1) = PrjD(i,rd+1)+I(k,j)*eff;
% end
% end
% end
t1 = clock;
%%S-L滤波器，R-L滤波器，新滤波器
hRL = zeros(2*M-1,1);
for i = 1:M-1
if mod(i,2) == 0
hRL(i) = -1/(pi^2*(M-i)^2*detD^2);
hRL(2*M-i) = hRL(i);
end
end
hRL(M) = 1/(4*detD^2);
% hSL = zeros(2*M-1,1);
% for i = 1:M-1
% hSL(i) = -2/(pi^2*(4*(M-i)^2-1)*detD^2);
% hSL(2*M-i) = hSL(i);
% end
% hSL(M) = 2/(pi^2*detD^2);
% hd = zeros(2*M-1,1); %%NEW滤波器
% for i = 1:M-1
% hd(i) = -1/(2*pi^2*(M-i)^2*detD^2);
% hd(2*M-i) = hd(i);
% end
% hd(M) = 1/(6*detD^2);
%%投影数据加权%%%%卷积%%
for p = 1:NumAngle
for q = 1:M
PrjD(p,q) = PrjD(p,q)*SD/sqrt(SD^2+((M/2-q+0.5)*detD)^2);
end
A = conv(PrjD(p,:),hRL);
PrjD(p,:) = A(M:2*M-1);
end
%%反投影%%
ObjPic = zeros(N,N);
for i = 1:NumAngle
for j = 1:N
for k = 1:N
XX = objX(j)*cosd(i)+objY(k)*sind(i); %%逆时针选择delta后，体素在新坐标系中的坐标
YY = -objX(j)*sind(i)+objY(k)*cosd(i);
dd = M/2-YY/(SO-XX)*SD/detD+0.5; %%在探测器上的位置
rd = floor(dd); %%取整
eff = dd-rd; %%双线性插值系数
if rd>0 && rd<M
ObjPic(k,j) = ObjPic(k,j)+(SO/(SO-XX))^2*(PrjD(i,rd)*(1-eff)+PrjD(i,rd+1)*eff);
end
end
end
end
t2 = clock;
etime(t2,t1);
ObjPic1 = ObjPic'/max(ObjPic(:));
figure(2); imshow(ObjPic1);
