%   BHPF
%   2/09
%   HE

close all; clear;
im=imread('C:\Hector\EE_552\Images\scene2\VIS_0064.bmp');
[M,N,c]=size(im);

D0=100;
D01=2*[D0,0];
D02=2*[0,D0];
DW=D0/4;
[x,y]=meshgrid(1:M,1:N);
n=sin(2*pi*x*D0/M)+sin(2*pi*y*D0/N);
imR=double(im(:,:,1))+50*n';
figure; imagesc(imR); colormap(gray);
save 'Test1_data2' imR

%Method Section 4.7.3
%Padd image and shift
P=2*M; Q=2*N; n=2;
fp=zeros(P,Q); fp(1:M,1:N)=imR;
fpp=zeros(P,Q);
fpp(P/2-M/2:P/2+M/2-1,Q/2-N/2:Q/2+N/2-1)=imR;
[x,y]=meshgrid(1:Q,1:P); onemtx=repmat(-1,P,Q);
fp=fp.*(onemtx.^(x+y));

%Compute Notch filters at
%D01=[100,0] & D02=[0,100]
Dk1=sqrt((x-Q/2-D01(1)).^2+(y-P/2-D01(2)).^2);
D_k1=sqrt((x-Q/2+D01(1)).^2+(y-P/2+D01(2)).^2);
Dk2=sqrt((x-Q/2-D02(1)).^2+(y-P/2-D02(2)).^2);
D_k2=sqrt((x-Q/2+D02(1)).^2+(y-P/2+D02(2)).^2);
HNR1=(1./(1+(DW./Dk1).^4)).*(1./(1+(DW./D_k1).^4));
HNR2=(1./(1+(DW./Dk2).^4)).*(1./(1+(DW./D_k2).^4));
HNR=HNR1.*HNR2;

F=fft2(fp);
G=F.*HNR;
figure; imagesc(log(abs(G))); colormap(gray);
g=real(ifft2(G)).*(onemtx.^(x+y)); g=g(1:M,1:N);
figure; imagesc(g); colormap(gray); title('High pass filtering with a BHPF');