% Isaac D. Gerg / Pushkar U. Durve
% CSE 485 - Digital Image Processing I 
% Project 3
% Group Number # 4

clc;
close all;

imgTemp = imread('lenna.gif', 'gif');
% Normalize Image
imgLenna = double(imgTemp) ./ 255;

%-------------------------------------------------------------------------------
% Part 1 - 2D FFT
	mxLennaSize = size(imgLenna);
	intRows = mxLennaSize(1);
	intCols = mxLennaSize(2);
	
	imgLennaSpectrum = fftshift(FFT2(imgLenna));
	
	% a)    log magnitude and phase
	imgLennaLogMag = log10(1+abs(imgLennaSpectrum));
	imgLennaPhase = angle(imgLennaSpectrum);
	
	% b)    Inverse FFT
	imgLennaInverseFFT = abs(IFFT2(imgLennaSpectrum));
	
	% c)    Plots
	figure('Name', 'Part 1', 'NumberTitle', 'off', 'MenuBar', 'none');
	colormap('gray');
	subplot(2,2,1);
	imagesc(imgLenna);
	title('Original Image (lenna.gif)');
	
	subplot(2,2,2);
	imagesc(imgLennaLogMag);
	title('Log Maganitude of FFT of lenna.gif');
	
	subplot(2,2,3);
	imagesc(imgLennaPhase);
	title('Phase of FFT of lenna.gif');
	
	subplot(2,2,4);
	imagesc(imgLennaInverseFFT);
	title('Inverse FFT');

%-------------------------------------------------------------------------------
% Part 2 - Lowpass Filter Design
    % a)
    LowpassFilter = Fspecial('gaussian', [11 11], 1.3);
    
    % b)  fft of lowpass filter
    LowpassFilterSpectrum = fftshift(fft2(LowpassFilter, intRows, intCols));
    
    % c)  plots of lowpass filter function in frequency domain
    figure('Name', 'Lowpass Filter', 'NumberTitle', 'off', 'MenuBar', 'none');
    colormap('gray');
    subplot(1,2,1);
    imagesc(log10(1+abs(LowpassFilterSpectrum)));
    title('Log Magnitude');
    subplot(1,2,2);
    imagesc(angle(LowpassFilterSpectrum));
    title('Phase');
    
    % d)  Filtering
    imgLennaFiltered = LowpassFilterSpectrum .* imgLennaSpectrum;
    figure('Name', 'Lowpass Filtered Image', 'NumberTitle', 'off', 'MenuBar', 'none');
    colormap('gray');
    subplot(2,2,1);
    imagesc(log10(1+abs(imgLennaFiltered)));
    title('Log Magnitude');
    subplot(2,2,2);
    imagesc(angle(imgLennaFiltered));
    title('Phase');

    % e) Inverse FFT
    subplot(2,2,3); 
    imgLennaLowpassFiltered = abs(ifft2(imgLennaFiltered));
    imgLennaLowpassFiltered = circshift(imgLennaLowpassFiltered,[-1.*floor(length(LowpassFilter)/2) -1.*floor(length(LowpassFilter)/2)]);
    imagesc(imgLennaLowpassFiltered);
    title('Inverse FFT of Lowpass Filtered Image');
    
%-------------------------------------------------------------------------------
% Part 3 - High Pass Filter Design
    % a) 
    HighpassFilter = Fspecial('laplacian'); % Using default alpha value.
    
    % b) fft of highpass filter
    HighpassFilterSpectrum = fftshift(fft2(HighpassFilter, intRows, intCols));
    
    % c) plots of highpass filter function in frequency domain
    figure('Name', 'Highpass Filter', 'NumberTitle', 'off', 'MenuBar', 'none');
    colormap(gray);
    subplot(1,2,1);
    imagesc(log10(1+abs(HighpassFilterSpectrum)));
    title('Log Magnitude');
    subplot(1,2,2);
    imagesc(angle(HighpassFilterSpectrum));
    title('Phase');
    
    % d) Filtering
    imgLennaFiltered = HighpassFilterSpectrum .* imgLennaSpectrum;
    figure('Name', 'Highpass Filtered Image', 'NumberTitle', 'off', 'MenuBar', 'none');
    colormap(gray);
    subplot(2,2,1);
    imagesc(log10(1+abs(imgLennaFiltered)));
    title('Log Magnitude');
    subplot(2,2,2);
    imagesc(angle(imgLennaFiltered));
    title('Phase');
    
    % e) Inverse FFT
    subplot(2,2,3);
    imgLennaHighpassFiltered = abs(ifft2(imgLennaFiltered));
    imgLennaHighpassFiltered = circshift(imgLennaHighpassFiltered,[-1.*floor(length(HighpassFilter)/2) -1.*floor(length(HighpassFilter)/2)]);
    imagesc(imgLennaHighpassFiltered);
    title('Inverse FFT of Highpass Filtered Image');

%-------------------------------------------------------------------------------
% Part 4 - 2D Filter Design
    % a)
        % 1. 
        [f1,f2] = freqspace(21,'meshgrid');
        r = sqrt(f1.^2 + f2.^2);
        Hd = ones(size(f1));
        Bandpass = Hd;
        Lowpass = Hd;
        Highpass = Hd;
        Bandpass((r<0.1) | (r>0.6)) = 0;
        Lowpass(r > 0.6) = 0;
        Highpass(r < 0.6) = 0;
        
    % b)
        % Fsamp2 Filter Design ------------------------------------------------
            BandpassFilter = fsamp2(Bandpass);
            LowpassFilter = fsamp2(Lowpass);
            HighpassFilter = fsamp2(Highpass);
            
            BandpassFilterSpectrum = fftshift(fft2(BandpassFilter, intRows, intCols));
            LowpassFilterSpectrum = fftshift(fft2(LowpassFilter, intRows, intCols));
            HighpassFilterSpectrum = fftshift(fft2(HighpassFilter, intRows, intCols));
            
            figure('Name', 'Bandpass Filter Spectrum (Fsamp2)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(1,2,1);
            imagesc(log10(1+abs(BandpassFilterSpectrum)));
            title('Log Magnitude');
            subplot(1,2,2);
            imagesc(angle(BandpassFilterSpectrum));
            title('Phase');
            
            figure('Name', 'Lowpass Filter Spectrum (Fsamp2)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(1,2,1);
            imagesc(log10(1+abs(LowpassFilterSpectrum)));
            title('Log Magnitude');
            subplot(1,2,2);
            imagesc(angle(LowpassFilterSpectrum));
            title('Phase');
	
            figure('Name', 'Highpass Filter Spectrum (Fsamp2)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(1,2,1);
            imagesc(log10(1+abs(HighpassFilterSpectrum)));
            title('Log Magnitude');
            subplot(1,2,2);
            imagesc(angle(HighpassFilterSpectrum));
            title('Phase');
            
            % Application of filters to image.
            imgLennaFilteredSpectrum = imgLennaSpectrum .* BandpassFilterSpectrum;
            figure('Name', 'Bandpass Filtered Image (Fsamp2)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(2,2,1);
            imagesc(log10(1+abs(imgLennaFilteredSpectrum)));
            title('Log Magnitude');
            subplot(2,2,2);
            imagesc(angle(imgLennaFilteredSpectrum));
            title('Phase');
            subplot(2,2,3);
            imgLennaFiltered = abs(ifft2(imgLennaFilteredSpectrum));
            imgLennaFiltered = circshift(imgLennaFiltered, [-1.*floor(length(BandpassFilter)/2) -1.*floor(length(BandpassFilter)/2)]);
            imagesc(imgLennaFiltered);
            title('Inverse FFT');
	
            imgLennaFilteredSpectrum = imgLennaSpectrum .* LowpassFilterSpectrum;
            figure('Name', 'Lowpass Filtered Image (Fsamp2)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(2,2,1);
            imagesc(log10(1+abs(imgLennaFilteredSpectrum)));
            title('Log Magnitude');
            subplot(2,2,2);
            imagesc(angle(imgLennaFilteredSpectrum));
            title('Phase');
            subplot(2,2,3);
            imgLennaFiltered = abs(ifft2(imgLennaFilteredSpectrum));
            imgLennaFiltered = circshift(imgLennaFiltered, [-1.*floor(length(LowpassFilter)/2) -1.*floor(length(LowpassFilter)/2)]);
            imagesc(imgLennaFiltered);
            title('Inverse FFT');
	
            imgLennaFilteredSpectrum = imgLennaSpectrum .* HighpassFilterSpectrum;
            figure('Name', 'Highpass Filtered Image (Fsamp2)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(2,2,1);
            imagesc(log10(1+abs(imgLennaFilteredSpectrum)));
            title('Log Magnitude');
            subplot(2,2,2);
            imagesc(angle(imgLennaFilteredSpectrum));
            title('Phase');
            subplot(2,2,3);
            imgLennaFiltered = abs(ifft2(imgLennaFilteredSpectrum));
            imgLennaFiltered = circshift(imgLennaFiltered, [-1.*floor(length(HighpassFilter)/2) -1.*floor(length(HighpassFilter)/2)]);
            imagesc(imgLennaFiltered);
            title('Inverse FFT');

        % Fwind1 Filter Design ------------------------------------------------
            BandpassFilter = fwind1(Bandpass, hamming(21));
            HighpassFilter = fwind1(Highpass, hamming(21));
            LowpassFilter = fwind1(Lowpass, hamming(21));
            
            BandpassFilterSpectrum = fftshift(fft2(BandpassFilter, intRows, intCols));
            LowpassFilterSpectrum = fftshift(fft2(LowpassFilter, intRows, intCols));
            HighpassFilterSpectrum = fftshift(fft2(HighpassFilter, intRows, intCols));
            
            figure('Name', 'Bandpass Filtered Image (Fwind1)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(1,2,1);
            imagesc(log10(1+abs(BandpassFilterSpectrum)));
            title('Log Magnitude');
            subplot(1,2,2);
            imagesc(angle(BandpassFilterSpectrum));
            title('Phase'); 
            
            figure('Name', 'Highpass Filtered Image (Fwind1)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(1,2,1);
            imagesc(log10(1+abs(HighpassFilterSpectrum)));
            title('Log Magnitude');
            subplot(1,2,2);
            imagesc(angle(HighpassFilterSpectrum));
            title('Phase'); 
            
            figure('Name', 'Lowpass Filtered Image (Fwind1)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(1,2,1);
            imagesc(log10(1+abs(LowpassFilterSpectrum)));
            title('Log Magnitude');
            subplot(1,2,2);
            imagesc(angle(LowpassFilterSpectrum));
            title('Phase'); 
            
            % Application of filters to image.
            imgLennaFilteredSpectrum = imgLennaSpectrum .* BandpassFilterSpectrum;
            figure('Name', 'Bandpass Filtered Image (Fwind1)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(2,2,1);
            imagesc(log10(1+abs(imgLennaFilteredSpectrum)));
            title('Log Magnitude');
            subplot(2,2,2);
            imagesc(angle(imgLennaFilteredSpectrum));
            title('Phase');
            subplot(2,2,3);
            imgLennaFiltered = abs(ifft2(imgLennaFilteredSpectrum));
            imgLennaFiltered = circshift(imgLennaFiltered, [-1.*floor(length(BandpassFilter)/2) -1.*floor(length(BandpassFilter)/2)]);
            imagesc(imgLennaFiltered);
            title('Inverse FFT');
	
            imgLennaFilteredSpectrum = imgLennaSpectrum .* LowpassFilterSpectrum;
            figure('Name', 'Lowpass Filtered Image (Fwind1)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(2,2,1);
            imagesc(log10(1+abs(imgLennaFilteredSpectrum)));
            title('Log Magnitude');
            subplot(2,2,2);
            imagesc(angle(imgLennaFilteredSpectrum));
            title('Phase');
            subplot(2,2,3);
            imgLennaFiltered = abs(ifft2(imgLennaFilteredSpectrum));
            imgLennaFiltered = circshift(imgLennaFiltered, [-1.*floor(length(LowpassFilter)/2) -1.*floor(length(LowpassFilter)/2)]);
            imagesc(imgLennaFiltered);
            title('Inverse FFT');
	
            imgLennaFilteredSpectrum = imgLennaSpectrum .* HighpassFilterSpectrum;
            figure('Name', 'Highpass Filtered Image (Fwind1)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(2,2,1);
            imagesc(log10(1+abs(imgLennaFilteredSpectrum)));
            title('Log Magnitude');
            subplot(2,2,2);
            imagesc(angle(imgLennaFilteredSpectrum));
            title('Phase');
            subplot(2,2,3);
            imgLennaFiltered = abs(ifft2(imgLennaFilteredSpectrum));
            imgLennaFiltered = circshift(imgLennaFiltered, [-1.*floor(length(HighpassFilter)/2) -1.*floor(length(HighpassFilter)/2)]);
            imagesc(imgLennaFiltered);
            title('Inverse FFT');
            
        % Fwind2 Filter Design ------------------------------------------------       
            window = fspecial('gaussian', 21, 2);
            window = window ./ max(max(window));
        
            BandpassFilter = fwind2(Bandpass, window);
            HighpassFilter = fwind2(Highpass, window);
            LowpassFilter = fwind2(Lowpass, window);
            
            BandpassFilterSpectrum = fftshift(fft2(BandpassFilter, intRows, intCols));
            LowpassFilterSpectrum = fftshift(fft2(LowpassFilter, intRows, intCols));
            HighpassFilterSpectrum = fftshift(fft2(HighpassFilter, intRows, intCols));
            
            figure('Name', 'Bandpass Filtered Image (Fwind2)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(1,2,1);
            imagesc(log10(1+abs(BandpassFilterSpectrum)));
            title('Log Magnitude');
            subplot(1,2,2);
            imagesc(angle(BandpassFilterSpectrum));
            title('Phase'); 
            
            figure('Name', 'Highpass Filtered Image (Fwind2)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(1,2,1);
            imagesc(log10(1+abs(HighpassFilterSpectrum)));
            title('Log Magnitude');
            subplot(1,2,2);
            imagesc(angle(HighpassFilterSpectrum));
            title('Phase'); 
            
            figure('Name', 'Lowpass Filtered Image (Fwind2)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(1,2,1);
            imagesc(log10(1+abs(LowpassFilterSpectrum)));
            title('Log Magnitude');
            subplot(1,2,2);
            imagesc(angle(LowpassFilterSpectrum));
            title('Phase'); 
            
            % Application of filters to image.
            imgLennaFilteredSpectrum = imgLennaSpectrum .* BandpassFilterSpectrum;
            figure('Name', 'Bandpass Filtered Image (Fwind2)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(2,2,1);
            imagesc(log10(1+abs(imgLennaFilteredSpectrum)));
            title('Log Magnitude');
            subplot(2,2,2);
            imagesc(angle(imgLennaFilteredSpectrum));
            title('Phase');
            subplot(2,2,3);
            imgLennaFiltered = abs(ifft2(imgLennaFilteredSpectrum));
            imgLennaFiltered = circshift(imgLennaFiltered, [-1.*floor(length(BandpassFilter)/2) -1.*floor(length(BandpassFilter)/2)]);
            imagesc(imgLennaFiltered);
            title('Inverse FFT');
	
            imgLennaFilteredSpectrum = imgLennaSpectrum .* LowpassFilterSpectrum;
            figure('Name', 'Lowpass Filtered Image (Fwind2)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(2,2,1);
            imagesc(log10(1+abs(imgLennaFilteredSpectrum)));
            title('Log Magnitude');
            subplot(2,2,2);
            imagesc(angle(imgLennaFilteredSpectrum));
            title('Phase');
            subplot(2,2,3);
            imgLennaFiltered = abs(ifft2(imgLennaFilteredSpectrum));
            imgLennaFiltered = circshift(imgLennaFiltered, [-1.*floor(length(LowpassFilter)/2) -1.*floor(length(LowpassFilter)/2)]);
            imagesc(imgLennaFiltered);
            title('Inverse FFT');
	
            imgLennaFilteredSpectrum = imgLennaSpectrum .* HighpassFilterSpectrum;
            figure('Name', 'Highpass Filtered Image (Fwind2)', 'NumberTitle', 'off', 'MenuBar', 'none');
            colormap(gray);
            subplot(2,2,1);
            imagesc(log10(1+abs(imgLennaFilteredSpectrum)));
            title('Log Magnitude');
            subplot(2,2,2);
            imagesc(angle(imgLennaFilteredSpectrum));
            title('Phase');
            subplot(2,2,3);
            imgLennaFiltered = abs(ifft2(imgLennaFilteredSpectrum));
            imgLennaFiltered = circshift(imgLennaFiltered, [-1.*floor(length(HighpassFilter)/2) -1.*floor(length(HighpassFilter)/2)]);
            imagesc(imgLennaFiltered);
            title('Inverse FFT');

%-------------------------------------------------------------------------------