function scan(filein,fileout)
% scan reads in an object density file (such as generated by 
% 'logo.m' and returns a scan data set.
% Parameters: input- and output file names,

% Read in data
   fid = fopen(filein,'r');
   if fid==-1 ; error('Input file could not be opened') ; end;
   mn  = fread(fid,2,'double'); nn = mn(1);
   a   = fread(fid,[nn,nn],'double');  fclose(fid);
   
   nn2  = (nn+1)/2; xv = [-nn2+1:nn2-1];  


n = input('Enter number of rays to sample the object with (odd number): ');
   if mod(n,2)==0 error('Number not odd'); end;
m = input('Enter number of angles to use: ');

% Create xv,xp1,xp2 describing basic (non-rotated) grid   
   n2 = (n+1)/2;  xw=[-n2+1:n2-1]; xz=xw*(nn+1)/(n+1); [xp1,xp2]=meshgrid(xz,xz);
   
% Loop over m angles
   for j=0:m-1
      th = pi*j/m; st = sin(th); ct = cos(th);
      x1  = xp1*ct-xp2*st;    % Generate the x1,x2-coordinates for
      x2  = xp1*st+xp2*ct;    % the grid points in the xp1,xp2 grid.
      s  = interp2(xv,xv,a,x1,x2,'cubic');
      s(isnan(s))=0;          % Get rid of all NaN:s from extrapolation.
      sc(:,j+1) = sum(s,2);   % Sum along each ray
   end
   sc = sc/max(max(sc));      % Normalize so max value = 1.
 
% Display scan data as surface plot plus grey scale image below
   gr = gray; colormap(gr(64:-1:1,:))
   mesh ([0:m-1],-xw,sc+2)
   axis ([0 m-1 -n2 n2 0 4])
   hold on
   pcolor ([0:m-1],-xw,sc*4); shading interp
   hold off
   xlabel('\theta'); ylabel('r')
   
% Write out scan data array   
   fid  = fopen(fileout,'w');
   if fid==-1 ; error('Output file could not be opened') ; end;
   fwrite (fid,[size(sc),sc(:).'],'double');
   fclose (fid);