% dip08_4.m             gr_nov2016  nov22
%   Rational Filter for edge-preserving noise smoothing
%   Rational filter for speckle noise: equiv. linear filter
%   Bilateral filter
clear all;

 A = imread('.\img\cameraman.tif');
% A = imread('.\img\lena.tif');
A = double(A);
[nr,nc] = size(A);
An = A + randn(nr,nc)*15;
An(An>255) = 255;  An(An<0) = 0;
% figure(1);  imshow(A,[0 255]);
err = (An-A).^2; 
psnr_init = 10 * log10( 256 .^2 / mean(err(:)) )

k = .1;   % edge protection factor  [.1]
W = 1;  %  linear filter tweaking factor [0 ... 20 or so]
itermax = 30;
B = An;
figure(1);   imshow(B,[0 255]);  pause()
psnr = zeros(1,itermax);
for iter=1:itermax
                            	d = [B(2:nr,:);B(nr,:)];
	a = [B(:,2:nc),B(:,nc)];    b = B;     c = [B(:,1),B(:,1:nc-1)];
                            	e = [B(1,:);B(1:nr-1,:)];
    B = b + (a+c-2*b)./(k*(a-c).^2+5+W) + (d+e-2*b)./(k*(d-e).^2+5+W);
    err = (B-A).^2; 
    psnr(iter) = 10 * log10( 256 .^2 / mean(err(:)) );
    figure(1);   imshow(B,[0 255]);  pause(.1);
end;
figure(2); plot(0.2,psnr_init,'x',1:itermax,psnr); pause()

[psnr_max,ind] = max(psnr)		% show max(psnr) result
B = An;
for iter=1:ind
                                d = [B(2:nr,:);B(nr,:)];
	a = [B(:,2:nc),B(:,nc)];    b = B;     c = [B(:,1),B(:,1:nc-1)];
                                e = [B(1,:);B(1:nr-1,:)];
    B = b + (a+c-2*b)./(k*(a-c).^2+5+W) + (d+e-2*b)./(k*(d-e).^2+5+W);
end;
figure(1);   imshow(B,[0 255]);
return


%%   Rational filter for speckle noise
% % filter response when k=0
% for A = 10:-1:6;
%    h = [1 1 1; 1 (A-8) 1; 1 1 1]./A;
%    mesh(fftshift(abs(fft2(h,256,256))));  pause(1)
% end
% % k>0 makes the denominator "locally larger" than A


%% bilateral filter
CC = [];
halfwindowsize = 5;
 sigma_r = .2;   for sigma_d = [1, 2, 4]
% sigma_d = 5;   for sigma_r = [.1, .2, .4]
C = bfilter2(An/255,halfwindowsize,[sigma_d, sigma_r])*255;
  %         ( input,     [5],      [  3,      .1    ]);
CC = [CC C];
 end;
figure(2);  imshow([An CC],[0 255]); 

	 
	 