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


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

%size of image
[im_height,im_width,rgbext] = 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
IMr = fftshift(fft2(im(:,:,1)));
IMg = fftshift(fft2(im(:,:,2)));
IMb = fftshift(fft2(im(:,:,3)));

%filter image with lowpass butterw (matrix dot product)

IMFr= IMr.*HPE;
IMFg= IMg.*HPE;
IMFb= IMb.*HPE;

%image reconstruction from spectra IMF (inverse fft)

imfr = real(ifft2(fftshift(IMFr)));
imfg = real(ifft2(fftshift(IMFg)));
imfb = real(ifft2(fftshift(IMFb)));

%saturation <0,1>
imfrs=imfr;
imfrs(any(imfr>0.999,3)) = 1;
imfrs(any(imfr<0.001,3)) = 0;

imfgs=imfg;
imfgs(any(imfg>0.999,3)) = 1;
imfgs(any(imfg<0.001,3)) = 0;

imfbs=imfb;
imfbs(any(imfb>0.999,3)) = 1;
imfbs(any(imfb<0.001,3)) = 0;



%allocate mem and generate RGB matrix
imf = ones(im_height,im_width,3);
imf(:,:,1) = real(imfrs(:,:));
imf(:,:,2) = real(imfgs(:,:));
imf(:,:,3) = real(imfbs(:,:));

%imagesc(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 HFE butterw. order=%d,D=%d',n,Do));

end

%print butter filter
figure(2);
colormap jet; colorbar;
mesh(HPE);
title(sprintf('HFE 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. order=%d,D=%d,a=%d,b=%d',n,Do,a,b));
figure(4);
imshow(im);
title(sprintf('original'));



