p_projection.m

平行束CT重建算法

function P=P_projection()
clear all;
close all;

Y_Length=40;
L=32;
N=128;
d=Y_Length/N;

ellip=sub_shepp;
P=zeros(L,N,10);

k=0;
for i=1:size(ellip,1)
    k=k+1;
x_0=ellip(i,1)*Y_Length/2;
y_0=ellip(i,2)*Y_Length/2;
A=ellip(i,4)*Y_Length/2;
B=ellip(i,5)*Y_Length/2;
rf=ellip(i,6)*pi/180;
rol=ellip(i,7);

   p=0;
   for theta=0:pi/L:pi-pi/L
   p=p+1;
   n=-N/2:N/2-1;
   C=2*rol*A*B;
   T=(A^2)*(cos(theta-rf))^2+(B^2)*(sin(theta-rf))^2;
   I=(n.*d-x_0*cos(theta)-y_0*sin(theta)).^2;
   G=(C*sqrt(T-I)./T).*(I<=T)+0.*(I>T);
   subplot(221);   
   plot(G);
   hold on;
   grid on;
   P1(p,:)=G;
   end

   P(:,:,k)=P1;
   hold on
   plot(P)
end

%filename='P_projection';
%save(filename,'P');