b_projection.m

平行束CT重建算法

function B=b_projection()
tic
Y_Length=40;
L=64;
N=128;
d=Y_Length/L;

F=zeros(N,N,10);

load h_projection;
m=0;
for k=1:10
m=m+1   
%反投影过程
p=0;
for theta=0:pi/L:pi-pi/L
    p=p+1;
  for i=1:N     %    图象离散为256*256
    for j=1:N
        
         x=[i-(N+1)/2]*d;
        
         y=[j-(N+1)/2]*d;
       
         r=sqrt(x^2+y^2);                                                 
         
      if x==0
           fi=pi/2||3*pi/2;
      elseif (x>0&&y>=0)
           fi=atan(y/x);
      elseif x<0
           fi=pi+atan(y/x);
      else fi=2*pi+atan(y/x);
          
      end
      
      
           R=(d*N/2+r*cos(fi-theta))/d;
           h(i,j)=R;
           n_0=floor(R);
           z(i,j)=n_0;
           delta=R-n_0;
             
       if n_0==0
           In_f=(1-delta/d)*Q(p,1,k)+delta/d*Q(p,2,k);
       elseif (n_0>0)&&(n_0<N)
           In_f=(1-delta/d)*Q(p,n_0,k)+delta/d*Q(p,n_0+1,k);
       else
           In_f=0;
       end  
           Add(i,j,k)=In_f;
        
    end
  end
  
F(:,:,k)=F(:,:,k)+Add(:,:,k);

end

f(:,:,k)=F(:,:,k);
end

filename='b_projection';
save(filename,'f');
toc