%  Find the direction for projection in the CT image
%  and show the projections as a film
%
%  Version 1.1, 7/11-97, Jørgen Arendt Jensen

Nx=max(size(img));
Nx2=Nx/2;
N=Nx;
M=Nx;
for theta=0:pi/10:pi

  %  Perform the projection

  x = ((1:N)'-N/2) * ones(1,M);
  y = ones(N,1) * ((1:M)-M/2);

  %  Make the index matrix

  xm = floor( x*cos(theta) + y*sin(theta) + N/2 + 0.5);
  ym = floor(-x*sin(theta) + y*cos(theta) + M/2 + 0.5);

  inside = (xm>0) & (xm<=N) & (ym>1) & (ym<=M);
  xm = xm.*inside + 1*(1-inside);
  ym = ym.*inside + 1*(1-inside);

  %  Find the projection

  ct_projection = sum (img(xm + (ym-1)*N) .* inside)';

  %  Plot the result

  axes(hp1)
  imagesc(img,'XData',[-0.5 0.5],'YData',[0.5 -0.5])
  axis xy
  colormap(gray); 
  axis('image');
  hold on
  r=(-400:400)/700;
  y=r*cos(theta);
  x=r*sin(theta);
  if ((abs(theta) > pi/4) & (theta<3*pi/4))
    for y1=(-Nx2:Nx2/6:Nx2)/Nx
      hl=plot(x,y+y1,'r');
      set(hl,'LineWidth',2);
      end
  else
    for x1=(-Nx2:Nx2/6:Nx2)/Nx
       hl=plot(x+x1,y,'r');
       set(hl,'LineWidth',2);
      end
    end
  hold off
  axis([-Nx2 Nx2 -Nx2 Nx2]/Nx)

  axes(hp2)
  Nc=length(ct_projection);
  plot(-(1:Nc)/Nc+0.5,ct_projection)
  xlabel('Normalized distance')
  ylabel('Projected value')
  drawnow
  pause(1)
  end


