% Filter Lenna to demonstrate band-limiting

% need to be in directory 
% cd c:\jl\4F8_imageprocessing\imageprocessing3\imageproc3mat\imageproc3demo
% imagefilter3_demo.m


clear,clc,close all,home

disp('Use the following for demonstration');
disp('Vertical cut-off frequencies FL=0.05, FU=.1');
disp('Horizontal cut-off frequencies FL=0.1, FU=.5');
disp('Single-sided filter order N=16');
disp('Note that zoom is turned on if required'); 
FL1=input('Lower normalised cut-off frequency in vertical direction = ');

% program to show the effects of filtering Lenna with an ideal rectangular
% bandpass filter

% need to be in directory c:\jl\I8_imageprocessing\imageprocessing3\imageproc3mat\imageproc3demo



FU1=input('Upper normalised cut-off frequency in vertical direction = ');

FL2=input('Lower normalised cut-off frequency in horizontal direction = ');
FU2=input('Upper normalised cut-off frequency in horizontal direction = ');

N1=input('Number of samples in single-sided impulse response = ');

n=(-N1:N1);
h=FU2*FU1*sinc(FU2*n)'*sinc(FU1*n);
h=(1/(pi^2))*(h-FL2*FL1*sinc(FL2*n)'*sinc(FL1*n));


figure(1)
load lenna256.mat
N=256;
LPx=filter2(h,x);
subplot(1,2,1),grayimage(x);
title('Original 256x256');
zoom on
subplot(1,2,2),grayimage(LPx);
%string=['Horizontal FL=',num2str(FL2),'  FU=',num2str(FU2),',  Vertical FL=',num2str(FL1),'  FU=',num2str(FU1)];
%title(string);
title('Filtered');
zoom on
pause;

% now do the alternative ideal rectangular filter
n=(-N1:N1);
h1=FU1*sinc(FU1*n) - FL1*sinc(FL1*n);
h2=FU2*sinc(FU2*n) - FL2*sinc(FL2*n);
h3=(1/(pi^2))*h1'*h2;


figure(2);
LPy=filter2(h3,x);
subplot(1,2,1),grayimage(x);
title('Original 256x256');
zoom on
subplot(1,2,2),grayimage(LPy);
%string=['Horizontal FL=',num2str(FL2),'  FU=',num2str(FU2),',  Vertical FL=',num2str(FL1),'  FU=',num2str(FU1)];
%title(string);
title('Filtered');
zoom on
pause;

% now do the ideal circularly symmetric bandpass filter
n=(-N1:N1);
for i=1:2*N1+1
 i1 = -pi + (i-1)*(2*pi/100); 
 for j=1:2*N1+1
   j1 = -pi + (j-1)*(2*pi/100); 
   k = sqrt(i1^2+j1^2);
   h4(i,j) = (1/(2*pi*k))*(FU2*besselj(1,FU2*k)-FL2*besselj(1,FL2*k));
 end
end


figure(3);
LPz=filter2(h4,x);
subplot(1,2,1),grayimage(x);
title('Original 256x256');
zoom on
subplot(1,2,2),grayimage(LPz);
%string=['Horizontal FL=',num2str(FL2),'  FU=',num2str(FU2),',  Vertical FL=',num2str(FL1),'  FU=',num2str(FU1)];
%title(string);
title('Filtered');
zoom on


