function rasterpic(thispic,rastersize,xclip,yclip) 
 
% RASTERPIC(thispic xclip,yclip,rastersize) imports the image 
% thispic which should be a grayscale 8-bit .gif image in the 
% search path.
% RasterPIC clips it to a central (xclip x yclip) region ...
% use for discussion of linear map, matrix representation
%
% transformed images are shown side-by-side... 
% 
% Author: Matthias Kawski, July 2001 
%         revised Dec 2003
%         http://math.la.asu.edu/~kawski   
% 


% Import an image -- the standard format is .tif 
% The following is a photoshop conversion of me.jpg 
%
if nargin < 1
    thispic='einstein.gif'
end
A=imread(thispic);
n=size(thispic,2)
BWnegative = min( thispic(1,n-2:n)=='gif')
pause
[M,N,R]=size(A);  % default: work w/ whole picture 

if nargin < 2    
   rastersize=4;    
end

% managing the clip region for large images
if nargin < 3    
   m = M;
else
   m = xclip;
end
if nargin < 4    
   n = N;
else
   n = yclip;
end

A=double(A); % convert from uint8 to integer 
% 
% To get a grey image take average of RGB channels 
% The following code is awkward -- due to my lack 
% of experience three dimensional m x n x r arrays 
% 
MAX=max(max(max(A)));
MIN=min(min(min(A)));
A=(A-MIN)/(MAX-MIN);

if BWnegative
    A = 1-A;
end
if R > 1
   size(A) 
   disp('flattening original RGB image to greyscale image')
   pp=(A(:,:,1)+A(:,:,2)+A(:,:,3))/3; 
else
   disp('original image is a greyscale image')
   pp=A(:,:,1);
end;
ppp(:,:,1)=pp; 
ppp(:,:,2)=pp; 
ppp(:,:,3)=pp;


image(ppp) 
% 
% clipping the image 
if m > M     
   m = M; 
end
if n > N     
   n = N;    
end

m=rastersize*floor(m/rastersize);
n=rastersize*floor(n/rastersize);
r=rastersize;

m0=floor((M-m)/2)+1; 
n0=floor((N-n)/2)+1; 
pp=pp(m0:m0+m-1,n0:n0+n-1,:); 
ppp=zeros(m,n,3); 
ppp(:,:,1)=pp; 
ppp(:,:,2)=pp; 
ppp(:,:,3)=pp; 
image(ppp) 
% 
disp('Hit any key to continue') 
pause 
% 
qq=zeros(size(pp));
for i=[0:m/r-1]
   for j=[0:n/r-1]
      for k=[1:r]
         for l=[1:r]
            qq(r*i+k,r*j+l)=sum(sum(pp(r*i+1:r*i+r,r*j+1:r*j+r)));
         end
      end
   end
end

qq = qq/r/r;
% 
ppp(:,:,1)=pp; 
ppp(:,:,2)=pp; 
ppp(:,:,3)=pp; 

qqq(:,:,1)=qq; 
qqq(:,:,2)=qq; 
qqq(:,:,3)=qq; 

vstrip=zeros(m,1,3); 
hstrip=zeros(1,2*n+3,3); 
image([hstrip;vstrip,ppp,vstrip,qqq,vstrip;hstrip]) 
%