--- Saving session to:
    ECE480.PH437_28-Jan-2002.txt
--- Processed startup.m ---
; ; ; ; ; ; ; 
edit fft2demo
inv3
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
All done!

edit inv3
fschange('c:\personal\class\2001-02\winter\ece480\matlab\inv3.m');
clear inv3
inv3
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
Press a key to continue...
All done!

type inv3

% EC480 Intro to Digital Image Processing (S00)
% (basically same as inv3.m from MA490 Advanced 
% Mathematical Methods in Image Processing, F98

% Interesting parameters
nstd=1;
Wmax=10000;

% Load the image
%x=rgb2gray(imread('btire.jpg'));
%x=rgb2gray(imread('hats256.jpg'));
x=imread('abcross.png');

%x=ones(256)*64;
%x(127:129,127:129)=192;
%x=uint8(x);

% Dimensions of the image
[numrows rowsize]=size(x);
n=numrows*rowsize;

% Display the image
f1=figure;
imshow(x),bigfig
bigtitle('Noise-Free Image: Spatial Domain')
xlabel('x'), ylabel('y')
disp('Press a key to continue...')
pause

% Display a row slice in the frequency domain
xmean=mean(x(:));
X = fftshift(fft2(double(x)-xmean)/n);
plot(abs(X(129,:)))
bigtitle('Noise-Free Image: Frequency Domain')
xlabel('x'), ylabel('y')
disp('Press a key to continue...')
pause


% Make a gaussian PSF
wn=31;

if 0
[xx yy]=meshgrid(-(wn-1)/2:(wn-1)/2,-(wn-1)/2:(wn-1)/2);
r = sqrt(xx.^2 + yy.^2);

% Define a radially symmetric Gaussian lowpass filter:
%  Full width at half maximum (FWHM) in pixels: fwhm  
%  (FWHM in frequency domain is n/fwhm)
fwhm=3;
h=exp(-(r.^2/(2*(fwhm/2.354)^2)));
h=h./sum(h(:));
end

% This is a square FIR filter the same coefficients throughout
if 1
h=ones(wn);
h=h./sum(h(:));
end

% Convert PSF into OTF (center the PSF in the padded image)
hpad=zeros(size(x));
loc=(numrows/2)+1-((wn-1)/2):(numrows/2)+1+((wn-1)/2);
hpad(loc,loc)=h;
hpad=fftshift(hpad);
%imshow(uint8((hpad/max(hpad(:)))*255))
H = fftshift(fft2(hpad));

% Display a row slice of the OTF
plot(abs(H(129,:)))
bigtitle('OTF Image: Frequency Domain')
xlabel('x'), ylabel('y')
disp('Press a key to continue...')
pause

% Use convolution to blur the original image
%g=filter2(h,x,'same');
%imshow(uint8(g))
%bigtitle('Blurred Image')
%xlabel('x'), ylabel('y')
%disp('Press a key to continue...')
%pause

% Convert the blurred image to the frequency domain
%G = fftshift(fft2(g)/n);

% Generate the blurred image in the frequency domain
G = H.*X;

% Display a row slice of the blurred image
plot(abs(G(129,:)))
bigtitle('Blurred Image: Frequency Domain')
xlabel('x'), ylabel('y')
disp('Press a key to continue...')
pause

% Convert blurred image to spatial domain
g = ifft2(fftshift(G))*n;
g = real(g)+xmean;
imshow(uint8(g),'notruesize')
bigtitle('Blurred Image')
xlabel('x'), ylabel('y')
disp('Press a key to continue...')
pause

% Add noise to the blurred image (see parameter
% at top of file)
%nstd=1;
%nstd=0.1;
%nstd=0;
%gn = imnoise(g,'gaussian',0,nvar);
gn = g + nstd*randn(size(g)) + 0;

imshow(uint8(gn),'notruesize')
bigtitle('Blurred + Noise Image')
xlabel('x'), ylabel('y')
disp('Press a key to continue...')
pause

% Convert the image to the frequency domain
gnmean=mean(gn(:));
GN = fftshift(fft2(gn-gnmean)/n);

% Display a row slice of the blurred+noise image
plot(abs(GN(129,:)))
bigtitle('Blurred+Noise Image: Frequency Domain')
xlabel('x'), ylabel('y')
disp('Press a key to continue...')
pause


% Compute the inverse filter function
W = 1./H;
%Wmax=100; %see parameter at top of file
W(find(abs(W)>Wmax))=Wmax;
%W(find(abs(W)>10000))=10000;

% Display a row slice of the inverse filter function
plot(abs(W(129,:)))
bigtitle('Pseudo-Inverse Filter')
xlabel('x'), ylabel('y')
disp('Press a key to continue...')
pause

% Apply the inverse filter
F=W.*GN;

% Display a row slice of the inverse filtered function
plot(abs(F(129,:)))
bigtitle('Restored Image: Frequency Domain')
xlabel('x'), ylabel('y')
disp('Press a key to continue...')
pause


% Convert result back to spatial domain
f = ifft2(fftshift(F))*n;
f = real(f)+gnmean;

% Clip negative values to zero
f(find(f<0))=0;

% Clip values to 255
f(find(f>255))=255;

% Convert to byte image
%f=uint8(f);

% Display the inverse-filtered image
imshow(uint8(f),'notruesize')
bigtitle('Restored Image -- Spatial Domain')
xlabel('x'), ylabel('y')
disp('Press a key to continue...')
pause

% Compare the observed image to the restored image
subplot(2,2,1)
imshow(x)
title('Noise-Free Image')
subplot(2,2,2)
imshow(uint8(f))
title('Restored Image')
subplot(2,2,3)
imshow(uint8(gn))
title('Observed Image')
disp('Press a key to continue...')
pause

figure
row=128;
plot(gn(row,:))
axis([0 255 0 255])
hold on
plot(f(row,:),'r')
%legx=100; legy=100; legoff=20;
%text(legx,legy,['Gaussian PSF, fwhm=' num2str(fwhm) ' pixels'])
%text(legx,legy-legoff,['Noise std=' num2str(255*sqrt(nvar)) ' gray levels'])
%text(legx,legy-2*legoff,['gamma=' num2str(gamma)])
%text(legx,legy-3*legoff,['||gn-Hf||=' num2str(rstd)])
%text(legx,legy-4*legoff,['row=' num2str(row)])
bigtitle('Inv: Observed & Restored Images')
hold off
disp('Press a key to continue...')
pause





%figure(f2);
%imshow(g),bigfig
%figure(f3);
%imshow(f),bigfig
%iblink([f2,f3])
%figure(f1)

if 0
   figure
   subplot(2,1,1)
   plot(u',H(128,:))
   title('MTF')
   subplot(2,1,2)
   plot(u',W(128,:),'g')
   a=axis;
   hold on
   plot(u',1./H(128,:))
   axis(a)
 end
   

disp('All done!')
return


exit