% dip08_3.m                gr 2016
%   Adaptive & frequency-selective linear filters
clear all;

f = double(imread('.\img\cameraman.tif'));
% f = double(imread('.\img\lenar.tif'));
[r,c] = size(f);

return

%% %%%%%%%%%%%%%%%% adaptive linear filter
sigN = 15;
fn = f + randn(r,c) * sigN;       %  additive Gaussian noise

w = [1 1 1; 1 1 1; 1 1 1] / 9;
g1 = filter2(w,fn);        % averaging filter (arithmetic mean) (3*3)

 sigqN = sigN.^2;          % adaptive 5*5 linear filter (dip08-17)
 g2 = fn;          
 for i=3:r-2;   for j=3:c-2;
   dum = fn(i-2:i+2,j-2:j+2);
   sigqL = std(dum(:)).^2;
   g2(i,j) = fn(i,j) - sigqN./sigqL .* (fn(i,j) - mean(dum(:)));
   end; end;

   % same, with automatic noise variance estimation (and fast realization)
locMean = filter2(ones(3,3)/9, fn);       % note this is a 3x3 filter
locVar = filter2(ones(3,3)/9, fn.^2) - locMean.^2;
nsVar = mean2(locVar);  % (over?)estimated noise variance
g3 = fn - min(nsVar./locVar,1).*(fn - locMean);
% [g3,noise] = wiener2(fn,[3 3]);     % equivalent Matlab function

 figure(1);  imshow([fn, g1; g2, g3],[0 255]);   %%%%%%%%%%%%%%%%%%%%%%%%%%%%
 PSNR_noisy = 10 * log10( 256^2 / (mean(mean((fn-f).^2))) )
 PSNR_g1 = 10 * log10( 256^2 / (mean(mean((g1-f).^2))) )
 PSNR_g2 = 10 * log10( 256^2 / (mean(mean((g2-f).^2))) )
 PSNR_g3 = 10 * log10( 256^2 / (mean(mean((g3-f).^2))) )

  return;

 %%  %%%%%%%%%%%%%%%%%%%%%%%    band-reject linear filter
        % only for square images, r=c 
[x,y] = meshgrid(1:r,1:c);    
fn = f + sin(x+.3*y)*10;     % add sinusoidal noise
Fn = fftshift(fft2(fn));
figure(1);  imshow(log(abs(Fn)),[]); pause();

d0 = 43; W = 7;  % sqrt( (1/pi)^2 + (.3/pi)^2 ) = 43/128
[u,v] = meshgrid(1:r,1:c);      %  build filter
z = sqrt((u-round(r/2)-1).^2 + (v-round(c/2)-1).^2); 
 H = double(z<d0-W/2 | z>d0+W/2);         %  ideal band-reject filter
% H = 1 ./ (1 + ((z .* W) ./ (z.^2 - d0.^2)).^2);  %  Butterworth filter, n=1
 figure(3);  mesh(abs(H));     % filter response
 figure(1);  imshow(log(abs(Fn))+abs(H),[]);     % filter superimposed on image transform
G = Fn .* H;
g = real(ifft2(ifftshift(G)));
figure(2); imshow([fn g],[]); 
 PSNR_noisy = 10 * log10( 256^2 / (mean(mean((fn-f).^2))) )
 PSNR_g = 10 * log10( 256^2 / (mean(mean((g-f).^2))) )


%% fft-domain periodic noise removal
clear all; close all;
figure(1);

a = double(imread('.\img\portrait1940.jpg'));
imshow(a,[]); pause;
r0 = 410;  c0 = 1;
a = a(r0:r0+399,c0:c0+749); % crop
[r,c] = size(a);

A = fftshift(fft2(a));
Am = abs(A);
imshow(log(Am),[]); pause;

map = ones(r,c);
map(Am>2e4) = 0;  % clear all large DFT components...
imshow(map,[]); pause;

lpr = 17; lpc = c/2-1;   % ...keeping those deemed useful
map(r/2-lpr:r/2+lpr,c/2-lpc:c/2+lpc) = 1;
lpr = 36; lpc = 18;
map(r/2-lpr:r/2+lpr,c/2-lpc:c/2+lpc) = 1;
imshow(map,[]); pause;

Bm = Am .* map;
imshow(log(Bm+1),[]); pause;

B = Bm .* exp(j*angle(A));
b = real(ifft2(ifftshift(B)));
imshow([a;b],[0 255]);

return;

