function nma_EECS_203A_HW4_prob_a

M=480;  % rows
N=640;  % columns
fileName='triangle.raw';
%img=imread(imageFileName,'ras');

 
fid=fopen(fileName,'r');
f=fread(fid,[N,M],'uint8');
fclose(fid);

f=f';

subplot(1,3,1);
imagesc(f,[0 255]);
colormap gray;
colorbar;

fprintf('generating fft2 for original image ....\n');
title(sprintf('Original raster file %s',fileName));
%center it
for(i=1:M)
    for(j=1:N)
        f(i,j)=f(i,j)*(-1)^(i+j);
    end
end

F=fft2(f);

%F = fftshift(F); % Center FFT

F2 = abs(F);  % Takes the magnitude of F
F2 = round(F2);
subplot(1,3,2);
imagesc(20*log(1+F2),[0 255]); colormap gray; 
title(sprintf('centerd FFT of original image %s',fileName));

subplot(1,3,3);
mesh(20*log(1+F2));
title(sprintf('centerd FFT of original, 3D'));


fprintf('creating Gaussian filter....\n');

%create filter
std=5; %pixels

h=zeros(M,N);
ii=0;
jj=0;
for(i=-M/2:M/2 -1)
    jj=0;
    ii=ii+1;
    for(j=-N/2:N/2 -1)
        jj=jj+1;
        h(ii,jj)=gaussian2D(i,j,std);
    end
end

fprintf('finished creating Gaussian filter....now FFT it...\n');

figure;
subplot(1,3,1);
mesh(h);
title(sprintf('Gaussian filter spatial, 3D'));
%center it
for(i=1:M)
    for(j=1:N)
        h(i,j)=h(i,j)*(-1)^(i+j);
    end
end
H=fft2(h);

subplot(1,3,2);
%H2 = fftshift(H); % Center FFT
H2 = abs(H);  % Takes the magnitude of F
H2 = round(H2);
imagesc(40*log(1+H2)); colormap gray; 
title(sprintf('centerd FFT of H'));

subplot(1,3,3);
mesh(20*log(1+H2)); 
title(sprintf('centerd FFT of Guasian filter'));


% width=6*std;
% hight=width;
% gaussianFilter=zeros(width,hight);
% ii=0;
% jj=0;
% for(i=-width/2:width/2 -1)
%     jj=0;
%     ii=ii+1;
%     for(j=-hight/2:hight/2 -1)
%         jj=jj+1;
%         gaussianFilter(ii,jj)=gaussian2D(i,j,std);
%     end
% end

fprintf('Now convolve h with f ...\n');

g=conv2(h,f);

size(g)

fprintf('Now generating FFT for degraded image ...\n');

figure;
subplot(1,3,1);
imagesc(g,[0 255]); colormap gray;
colormap gray;
colorbar;
title('degraded image, covolution of Gaussian with original image');

% NO !!
%center it
%for(i=1:M)
%    for(j=1:N)
%        g(i,j)=g(i,j)*(-1)^(i+j);
%    end
%end
G=fft2(g);

%G2 = fftshift(G); % Center FFT
G2 = abs(G);  % Takes the magnitude of F
G2 = round(G2);
subplot(1,3,2);
imagesc(4*log(1+G2)); colormap gray; 
title(sprintf('G: centerd FFT of degraded image'));

subplot(1,3,3);
mesh(4*log(1+G2));
title(sprintf('Mesh centerd FFT of degraded image'));

fprintf('Now generating restored image ...\n');

FF=zeros(M,N);

ii=0;
jj=0;
for(i=1:M)
    ii=i+M/2;
    jj=0;
    for(j=1:N)
        jj=j+N/2;
        if(abs(real(H(i,j)))<=realmin)
            realPart=0;                
        else
            realPart=real(G(ii,jj))/real(H(i,j));
        end
        
        if(abs(imag(H(i,j)))<=realmin)
            imagPart=0;                
        else
            imagPart=imag(G(ii,jj))/imag(H(i,j));
        end
            
        FF(i,j)=realPart+sqrt(-1)*imagPart;
    end
end

%cutoff
FF_cut=FF(200:end-200,240:end-240);

figure;
subplot(1,4,1);
imagesc(5*log(1+(round(abs(FF_cut))))); colormap gray; 
title(sprintf('F_hat: centerd FFT of restored image'));

subplot(1,4,2);
mesh(4*log(1+FF_cut)); colormap gray; 
title(sprintf('centerd mesh FFT'));


%do a cut off
ff=ifft2(FF_cut);
ff=real(ff);
%decenter it
[M,N]=size(ff);
for(i=1:M)
    for(j=1:N)
        ff(i,j)=ff(i,j)*(-1)^(i+j);
    end
end

 
subplot(1,4,3);
imagesc(ff,[0 255]);colormap gray;
title('restored image, centered');

subplot(1,4,4);
imagesc(ff);colormap gray;
title('restored image, no scale');




c=22;
imagesc(c*log(1+F),[0 255]); colormap gray; 
title(sprintf('|F| for triangle. Scaled with c log(1+F), c=%d',c));
xlabel('v'); ylabel('u');

figure;
mesh(1:N,1:M,.6*log(1+F)); colormap jet;
xlabel('v'); ylabel('u');
title(sprintf('3D view of |F| for triangle. Scaled with c log(1+F), c=%d',c));


figure;
c=1;
for (i=1:25)
     subplot(5,5,c);
     imagesc(i*log(1+F),[0 255]); colormap gray; title(sprintf('c=%d',c));
     c=c+1;
end

 %%FF=zeros(M,N); FF=G(240:end-240,320:end-320) ./ H;
 
%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%
function v=gaussian2D(x,y,std)
v=1/(sqrt(2*pi)*std)*exp(-(x^2+y^2)/(2*std^2));
