
% Note you also need to put Atrium2.jpg in your current work directory

function lab4()
%% exercise 1
im=zeros(128,128);
im(40:89,63:66)=1;
subplot(3,5,1);imshow(im);

f1=fft2(im);
f1=fftshift(f1);

fabs1=FH_abs(f1);
subplot(3,5,2);imshow(fabs1);title('magnitude');

fangle1=angle(f1);
fangle1=(fangle1-min(fangle1(:)))/(max(fangle1(:))-min(fangle1(:)));
subplot(3,5,3);imshow(fangle1);title('phase angles');

%% exercise 2
im2=zeros(128,128);
im2(20:69,33:36)=1;
subplot(3,5,6);imshow(im2);

f2=fft2(im2);
f2=fftshift(f2);

fabs2=FH_abs(f2);
subplot(3,5,7);imshow(fabs2);title('magnitude');

fangle2=angle(f2);
fangle2=(fangle2-min(fangle2(:)))/(max(fangle2(:))-min(fangle2(:)));
subplot(3,5,8);imshow(fangle2);title('phase angles');


im3=imrotate(im,45,'crop');
subplot(3,5,11);imshow(im3);

f3=fft2(im3);
f3=fftshift(f3);

fabs3=FH_abs(f3);
subplot(3,5,12);imshow(fabs3);title('magnitude');

fangle3=angle(f3);
fangle3=(fangle3-min(fangle3(:)))/(max(fangle3(:))-min(fangle3(:)));
subplot(3,5,13);imshow(fangle3);title('phase angles');

%% exercise 3
im4=imread('Atrium2.jpg');
im_grey=sum(im4,3);
im_grey=im_grey/max(im_grey(:));
subplot(3,5,4);imshow(im_grey);title('grey image');

f_im=fft2(im_grey);
f_im=fftshift(f_im);

fabs_im=FH_abs(f_im);
subplot(3,5,5);imshow(fabs_im);title('magnitude');

fangle_im=angle(f_im);
fangle_im=(fangle_im-min(fangle_im(:)))/(max(fangle_im(:))-min(fangle_im(:)));
% imshow(fangle_im);title('phase angles');

f_im=fft2(im_grey);
f_im=fftshift(f_im);
fangle_im=angle(f_im);
fabs_im=abs(f_im);
radius=40;
[width,height]=size(im_grey);
for i=1:width
    for j=1:height
        if sqrt((i-width/2)^2+(j-height/2)^2)>radius
            fabs_im(i,j)=0;  % low pass filter
        end
    end
end
subplot(3,5,10);imshow(fabs_im);title('low pass filter');

f_recover=FH_complex(fabs_im,fangle_im);
im_recover=ifft2(fftshift(f_recover));
subplot(3,5,9);imshow(im_recover); title('image recovered');

fabs_im=abs(f_im);
radius=40;
[width,height]=size(im_grey);
for i=1:width
    for j=1:height
        if sqrt((i-width/2)^2+(j-height/2)^2)<radius
            fabs_im(i,j)=0;  % high pass filter
        end
    end
end
subplot(3,5,15);imshow(fabs_im/max(fabs_im(:)));title('high pass filter');

f_recover=FH_complex(fabs_im,fangle_im);
im_recover=ifft2(fftshift(f_recover));
subplot(3,5,14);imshow(im_recover);title('image recovered');

end







function fabs=FH_abs(f)
% Output the magnitude for complex matrix f. The magnitude goes throught
% the logarithm enhancement and the normalization procedure, in order to
% display with a higher contrast.
fabs=abs(f);
fabs=log(1+fabs);   % logarithm enhancement
fabs=(fabs-min(fabs(:)))/(max(fabs(:))-min(fabs(:)));
end





function c=FH_complex(magnitude,angle)
%FH_complex Recover the complex matrix from its magnitude part and phase part 
% 
% [ Syntax ]
%   - c = FH_complex(magnitude, angle);
%
% [ Description ]
%
%	angle: ranges from [-pi,pi], it is usually generated from atan2
%           function which is used in angle function
%
%   The relation between magnitude, angle, the real part and imaginary part are as follows:
%       
%   magnitude = sqrt(real(X).^2 + imag(X).^2)
%   angle = atan2(imag(X),real(X))
%
% [ Demo usage ]
%
%   im=rand(100,100);
%   f_im=fft2(im);
%   mag=abs(f_im);
%   ang=angle(f_im);
%   f_recover=FH_complex(mag,ang);
%   difference1=max(max(real(f_recover)-real(f_im)));
%   difference2=max(max(imag(f_recover)-imag(f_im)));
%   fprintf('The maximum difference between the original complex matrix and the recovered matrix is %f and %f.\n',difference1,difference2);
%

assert(isreal(magnitude));
assert(isreal(angle));

coefficent = tan(angle);
coefficent = abs(coefficent);

real=sqrt(magnitude.^2./(1+coefficent.^2));
imag=real.*coefficent;

index=find(angle<0);
imag(index)=-imag(index);

index=find(angle < -pi/2 | angle > pi/2);
real(index)=-real(index);

c=complex(real,imag);

end