基于MATLAB的CT直接反投影算法的实验研究

推导参考我的上一篇博客基于MATLAB的CT平行束投影数据的仿真
实验任务
1、利用“iradon”函数进行直接反投影重建。

I=phantom(256);
theta=0:179;
P=radon(I,theta);
rec=iradon(P,theta,'None');
figure;
imshow(I,[]),title('原始图像');
figure;
imshow(I,[]),title('直接反投影图像');

2、采用本节推导出来的直接反投影算法进行重建。

clear;
N=256;
I=phantom(N);
delta=pi/180;
theta=0:1:179;
theta_num=length(theta);
%%产生投影数据
P=radon(I,theta);
[mm,nn]=size(P);
e=floor((mm-N-1)/2+1)+1;
P=P(e:N+e-1,:);
P1=reshape(P,N,theta_num);
%%投影重建
rec=zeros(N);
for m=1:theta_numpm=P1(:,m);Cm=(N/2)*(1-cos((m-1)*delta)-sin((m-1)*delta));for k1=1:Nfor k2=1:NXrm=Cm+(k2-1)*cos((m-1)*delta)+(k1-1)*sin((m-1)*delta);n=floor(Xrm);t=Xrm-floor(Xrm);n=max(1,n);n=min(n,N-1);p=(1-t)*pm(n)+t*pm(n+1);rec(N+1-k1,k2)=rec(N+1-k1,k2)+p;endend
end
figure;
imshow(I,[]),title('原始图像');
figure;
imshow(rec,[]),title('解析法直接反投影图像');