clear all
% input parametrs:
n = 2; Do = 5; %filter parameters
a = 1; b = 1.5;


%load image
imrgb = imread('Pevnost_na_skaleg.jpg');
im = im2double(imrgb);

%size of image
[im_height,im_width] = size(im);

%euclidean distance matrix
[U,V] = rc2uv(im_height,im_width);
D = sqrt(U.^2+V.^2); 

%generate butterworth filter from euclidean distance matrix
H = 1./(1+(D./Do).^(2*n));

%high pass
H = 1-H;

A = ones(im_height,im_width);

HPE = a.*A + b.*H;

%spectra of image
IM = fftshift(fft2(im));

%filter image with lowpass butterw (matrix dot product)

IMF=IM.*HPE;

%image reconstruction from spectra IMF (inverse fft)

imf = ifft2(fftshift(IMF));

if 0
%print output:
figure(1); clf;


subplot(1,2,1);


%print input figure
imshow(im);
title('original figure');

%print output figure
subplot(1,2,2);
imshow(imf);
title(sprintf('filtered figure with HPE butterw. order=%d,D=%d,a=%d,b=%d',n,Do,a,b));

end

%print butter filter
figure(2);
colormap jet; colorbar;
mesh(HPE);
title(sprintf('HPE filter, param order=%d,diam(D)=%d,a=%d,b=%d',n,Do,a,b));



%whole fullscale
figure(3);
imshow(imf);
title(sprintf('filtered figure with HFE butterw. \n order=%d,D=%d,a=%d,b=%d',n,Do,a,b));
figure(4);
imshow(im);
title(sprintf('original'));

%writefile:
imwrite(imf,'out.png','png','BitDepth',16);

