ctfbp.m

CT图像重建的变换类算法

%--卷积反投影重建
%--p投影值，K--投影方向数，N--重建点数N^2,S--每方向投影数

function F=ctfbp(p,K,N,S)
tic
I=K*S;
M=N;
J=M*N;
d=1/(N-1);d1=1/(S-1);
%--滤波函数h(n)
h=[];
for n=1:2*S-1
   h(n)=-2/(pi^2*(1/(2*5))*(4*(n-S)^2-1));
end
%--滤波后投影c1
for n=1:K
    c(:,n)=conv(p(:,n),h');
end
c1=c(S:2*S-1,:);
%--重建图像F
p2=[];
for m=0:K-1
    xr=0;
for i=1:N
for j=1:N
    xr=xr+1;
xrm(i,j)=((j-1)*d-(N-1)/2*d)*cos(m*pi/K)-((i-1)*d-(N-1)/2*d)*sin(m*pi/K)+(N-1)/2*d;
n0=xrm(i,j)/d1-mod(xrm(i,j)/d1,1)+1;

if n0>0&&n0<(S-1)
p2(xr,m+1)=(1-mod(xrm(i,j)/d1,1))*c1(n0,m+1)+mod(xrm(i,j)/d1,1)*c1(n0+1,m+1);
else p2(xr,m+1)=0;
end

end
end
end
F=sum(p2');

for j=1:J
    if F(j)<0
        F(j)=0;
    end
end

F=reshape(F,N,N);
x=-0.5:1/(N-1):0.5;
y=-0.5:1/(N-1):0.5;
[X,Y]=meshgrid(x,y);
figure
mesh(X,Y,F)
toc