function out=fourierme(NN,myimage)
% FOURIERME(NN) Naive implementation of Fourier approximation
% of a picture, imported as a .tif image. The parameter NN is
% the degree of highest order desired Fourier approximant.
% myimage is the name of a tif file -- if none is supplied
% the file 'oldme.tif' will be used. 
%
% This implementation is for demo-purposes only, and is
% highly inefficient -- e.g. it stores the function values
% of all cos(mx)*cos(ny) at all grid-points and takes NO
% advantage of FFT ideas. Furthermore no attention is paid
% to correct scaling or to lack of orthogonality of the 
% sinusoids on a discrete grid....
%
% The original version imports a 60 x 80 RGB tiff image, 
% converted it into grey-scale, truncated the whites along
% the edges. Next it tabulates the values of cos(nx), sin(nx)
% and cos(my), sin(my) on the same grid. The Fourier coeffs
% are computed in the naivest way and stored in a 
% (2*NN+1)x(2*NN+1) matrix with negative indices correspon
% ding to sines, postive ones to cosines (and obvious pre-
% view of complex exponentials).
% After calculating all Fourier coeff's up to order NN, 
% images corresponding to the k-th Fourier approximants,
% k=0,1,2,..NN are displayed.
% The spectrum, i.e. Fourier coeff's are returned as a
% (2*NN+1)x(2*NN+1) matrix.
%
% The main purpose is as a demo in a first linear algebra 
% class, to illustrate the idea of an image as a vector,
% and a Fourier approx as a projection onto a subspace.
% Some particularly nice connections are made by displaying
% in side-by-side the surface and its top-view (=image).
%
% All rights reserved: Matthias kawski, July 2000
%                      http://math.la.asu.edu/~kawski
if nargin < 1
   % Self-demonstrating. Degree of highest approximant.
   NN=11
end
if nargin < 2
   % Self-demonstrating. Default image.
   myimage='me.tif';
end
myface=imread(myimage);
SU=figure
set(SU,'Position',[50 300 300 300])
image(myface)
title(strcat('The original image. Size=',...
              num2str(size(myface))))
pause
IM=figure
set(IM,'Position',[350 300 300 300])
% myface=double(myface(24:60,14:50,:));
myface=double(myface);
greyface=(myface(:,:,1)+myface(:,:,2)+myface(:,:,3))/256/3;
colormap(gray)
figure(IM)
image(64*greyface)
title('Converted into grey-scale. Using colormap(grey)')
pause
figure(SU)
surf(64*greyface)
colormap(gray)
title('Rotate the image into near-top view!')
view(-40,30)
pause
greyface=min(0.6,greyface);
surf(64*greyface)
title('Truncated the bright areas in the background!')
pause
figure(IM)
image(64*greyface)
title('The truncated images -- grey-er background!')
pause
Nx=size(greyface,2);
Ny=size(greyface,1);
Nx0=floor(Nx/2);
Ny0=floor(Ny/2);
[xx,yy]=meshgrid(-Nx0:-Nx0+Nx-1,-Ny0:-Ny0+Ny-1);
xx=pi/Nx0*xx;
yy=pi/Ny0*yy;
oo(:,:,1)=ones(size(xx));
for k=[1:NN] 
   cx(:,:,k)=cos(k*xx);
   sx(:,:,k)=sin(k*xx);
   cy(:,:,k)=cos(k*yy);
   sy(:,:,k)=sin(k*yy);
end;
FC=zeros(2*NN+1,2*NN+1);
% the coefficient of the (0,0) frequency, cos(0x)*cos(0y)
FC(NN+1,NN+1)=sum(sum(greyface))/4;
for i=[1:NN]
   % coeff's of the (k,0) frequencies, cos(kx)*cos(0y)
   FC(NN+1+i,NN+1)=sum(sum(cx(:,:,i).*greyface/2));
   % coeff's of the (-k,0) frequencies, sin(kx)*cos(0y)
   FC(NN+1-i,NN+1)=sum(sum(sx(:,:,i).*greyface/2));
   % coeff's of the (0,k) frequencies, cos(0x)*cos(ky)
   FC(NN+1,NN+1+i)=sum(sum(cy(:,:,i).*greyface/2));
   % coeff's of the (0,-k) frequencies, cos(0x)*sin(ky)
   FC(NN+1,NN+1-i)=sum(sum(sy(:,:,i).*greyface/2));
   for j=[1:NN]
      % coeff's of the (i,j) frequencies, cos(ix)*cos(jy)
      FC(NN+1+i,NN+1+j)=sum(sum(cx(:,:,i).*cy(:,:,j).*greyface));
      % coeff's of the (-i,j) frequencies, sin(ix)*cos(jy)
      FC(NN+1-i,NN+1+j)=sum(sum(sx(:,:,i).*cy(:,:,j).*greyface));
      % coeff's of the (i,-j) frequencies, cos(ix)*sin(jy)
      FC(NN+1+i,NN+1-j)=sum(sum(cx(:,:,i).*sy(:,:,j).*greyface));
      % coeff's of the (-i,-j) frequencies, sin(ix)*sin(jy)
      FC(NN+1-i,NN+1-j)=sum(sum(sx(:,:,i).*sy(:,:,j).*greyface));
   end
end
% FC=FC/N/N;
figure(IM)
image(FC)
title('Spectrum')
figure(SU)
surf(FC)
colormap(hot)
title('Spectrum')
pause
Fapp=FC(NN+1,NN+1)/Nx/Ny*4*oo;
figure(IM)
image(64*Fapp)
title('Order of approximation N=0')
figure(SU)
surf(64*Fapp)
colormap(gray)
title('Order of approximation N=0')
pause
for k=[1:NN]
   Fapp=Fapp+FC(NN+1,NN+1+k)*cy(:,:,k);
   Fapp=Fapp+FC(NN+1,NN+1-k)*sy(:,:,k);
   Fapp=Fapp+FC(NN+1+k,NN+1)*cx(:,:,k);
   Fapp=Fapp+FC(NN+1-k,NN+1)*sx(:,:,k);
   for j=[1:k]
       Fapp=Fapp+FC(NN+1+k,NN+1+j)*cx(:,:,k).*cy(:,:,j);
       Fapp=Fapp+FC(NN+1+k,NN+1-j)*cx(:,:,k).*sy(:,:,j);
       Fapp=Fapp+FC(NN+1-k,NN+1+j)*sx(:,:,k).*cy(:,:,j);
       Fapp=Fapp+FC(NN+1-k,NN+1-j)*sx(:,:,k).*sy(:,:,j);
   end
   %avoid double couning the corners
   for i=[1:k-1]
       Fapp=Fapp+FC(NN+1+i,NN+1+k)*cx(:,:,k).*cy(:,:,i);
       Fapp=Fapp+FC(NN+1-i,NN+1+k)*sx(:,:,k).*cy(:,:,i);
       Fapp=Fapp+FC(NN+1+i,NN+1-k)*cx(:,:,k).*sy(:,:,i);
       Fapp=Fapp+FC(NN+1-i,NN+1-k)*sx(:,:,k).*sy(:,:,i);
   end
   mm=min(min(Fapp))
   MM=max(max(Fapp))
   figure(SU)
   surf((Fapp-mm)/(MM+1-mm)*64)
   title(strcat('Order of approximation N=',num2str(k)))
   figure(IM)
   image((Fapp-mm)/(MM+1-mm)*64)
   title(strcat('Order of approximation N=',num2str(k)))
   pause
end
disp('Strike enter to clear the figures and end the demo')
pause
delete(IM)
delete(SU)
out=FC;