%% Lab 5: Convolution in 1 and 2-D
load('lab5_files.mat'); 
% ** SOLUTIONS ** 

%% Exercise 1: The mechanics of convolution

% Question 1: When you use the conv or conv2 function in matlab, you are 
% implicitly using zero padding by default. Fill in the function 
% manualPadding at the end of this file that replicates the native conv
% function, but by manually applying the appropriate zero padding and using
% the 'valid' mode only.
%
% Bonus: Try replicating conv by building the convolution matrix from
% scratch for each setting, instead of using 'valid'.


% Question 2: Fill in the skeleton code below, then interpret the resulting
% plots. Specifically, think about what it means for the kernel size to
% increse (in this specific case, and in general). 

x = 0:0.01:2*pi;
y = sin(x) + 0.1*randn(size(x));

figure;
set(gcf, 'Position', [0, 0, 1500, 400]);

subplot(1, 4, 1); 
plot(x, y);
title('Raw Signal');
counter = 2;
for kernel_sz = [5, 10, 100]
    % define a kernel for the convlution that takes a rolling average of
    % the signal
    
    % ** YOUR CODE HERE **
    
    % ** ANSWER **
    k = ones(1, kernel_sz) .* (1/kernel_sz);
    
    % convolve the noisy sine wave y with the kernel using 'same'
    % padding
    
    % ** YOUR CODE HERE **
    
    % ** ANSWER **
    smoothed_y = conv(y, k, 'same'); 
    
    % visualize the results
    subplot(1, 4, counter);
    plot(x, smoothed_y); 
    title(['Averager Kernel Size = ', num2str(kernel_sz)]);
    counter = counter + 1;
    
end
    

% Question 3: Interpret the meaning of larger vs. smaller kernel size for
% convolutions with time series data. (In the general case, the above 
% is one possible example). 
%, 
% ** ANSWER **
% Increasing the kernel size increases the number of elements in the
% input that have an effect on any given output. In the case of time series
% data this means that events that input events from further in the past
% effect the output (input has longer temporal persistence on output).
% Or,if the system is not causal it could mean the output at a given time 
% depnds on inputs that are further in the future
%
% For example a kernel that outputs a time centered average, 
% so output[t] = input[t - 1] + input[t + 1] would not be causal.

%% Example 1: Feature Extraction 
% Convolving a 2-D signal with a kernel compares the similarity of a pieces
% of the input with the kernel (in the dot product sense of similarity).
% The image that the convolution results in can be thought of as a heat-map
% of how similar a part of the input is with the kernel. In computer vision
% this operation is called "feature extraction," with the kernel being the
% "feature."
%
% Fill in the code below to implement the procedure described above for
% features called 2-D "gabor functions," which are often used to model V1
% receptive fields.

clear all; 
load('lab5_files.mat');
% Gabor parameters 
wavelength = 20; % scales size of kernel
theta = pi/4; % changes the orientation of the kernel
x = gabor_fn(theta, wavelength);

% visualize some gabors by messing with the value of theta above
% (note: you can try highlighting and running this section of the code by 
% itself to compare figures)
figure;
set(gcf, 'Position', [100, 100, 500, 500]);
imagesc(x);

% display input image (loaded at top of script)
figure; 
imagesc(feature_input_image); 
title('Feature Extraction Input Image');


% compute model convolutions with kernels at different angles/sizes
% (note the difference between 3*pi/4 and pi/4, at least)
theta = 3*pi/4; % controls orientation
wavelength = 30; % scales the size of the kernel
kernel = gabor_fn(theta, wavelength);
responses = 0; % REPLACE THIS LINE
responses = conv2(feature_input_image, kernel); % ** ANSWER ** 
          
% visualize the results
figure; 
imagesc(abs(responses)); 
title(['Convolution with Theta= ', num2str(theta),...
    ', Wavelength= ', num2str(wavelength)]);
  
% Run the above code fora  variety of choices of theta. 
% What effect does changing the size and orientation of the 
% kernel have on the output. Can you link this to the analogy 
% between V1 cells receptive fields and convolution that was hinted
% at/discussed?


% ** ANSWER ** 
% The convolution is bright where the image is similar to the kernel, so
% changing the orientation changes the pattern in the output. For example,
% there is a long edge in the image with orientation ~3pi/4 that is
% highligted in the output with a filter at this orientation (which has 
% an edge at that orientation by definition), but dimmed in the output when
% the kernel orientationis set to pi/4, which is approximately orthogonal
% to the orientation of the long edge in the image. 
%
% Increasing the wavelength of the gabor/size of the kernel changes the
% size of the features the output is sensitive to. One effect of this is
% that more distant elements of the input are connected to each other in
% the convolution matrix, so the output features appear "blurred," or
% pooled together. 


%% Example 2: Recovering the PSF

% You notice that the images coming out of your imaging system seem like
% they might be distorted. You decide to test this by using the system to
% measure/image an input that you already know the value of (this test 
% input image is loaded and named X below), and the output you get is the
% matrix/image Y.
% 
% In the example below, we know the point spread function that is
% associated with some imaging system, and so given any input signal the
% output of the imaging system is the convolution of the input with the
% point spread function. Below an example input output pair are visualized
% along with the point spread function. 
%
% Note: you can manually change the figure size by adjusting the size of
% the window that pops up.
figure;
subplot(1, 3, 2)
imagesc(PSF_input_image); 
title('Imaging System Input Image');

% psf is gaussian (largest in middle and falls off in all directions)
subplot(1, 3, 1)
psf = fspecial('gaussian', 6, 3);  
imagesc(psf); colorbar;  
title('Imaging System PSF')

% output is input convolved with psf
subplot(1, 3, 3)
output_image = conv2(PSF_input_image, psf); 
imagesc(output_image); 
title('Imaging System Output Image');

% In the real world you may not know a-priori a measurement 
% system's point spread function. Explain a strategy to estimate the
% PSF/covoluition kernel of an imaging system using linear regression,
% given access to known inputs (X) and their outputs (Y). You can explain
% in words, or pseudocode, or actual code, or in any other way. 

% ** ANSWER ** 
% Each element of the observed output is the dot product of the kernel
% vector (will call the kernel vector "r") with a piece of the input
% signal. Writing this statement out as a matrix equation we have 
% y = Xr  with r being the kernel vector, and with  the i_th row of X 
% being the piece of the input the kernel is dotted with to get the i_th
% output.
% 
% Then, since r is an unknown, we can find its least squares estimate given
% our data by multiplying each side by the pseudoinverse of X (on the
% left), so r_est = pseudoInverse(X) * y.

%% Exercise 1 Function
% NOTE: For the purposes of this functin we will assume both signal and
% kernel are column vectors
function output = manualPadConv(signal, kernel, pad_mode)
    if strcmp(pad_mode, 'valid')
        output = conv(signal, kernel, 'valid');
    elseif strcmp(pad_mode, 'full')
        % pad the input signal with the correct amount of zeros for full
        num_pad = 0; % ** replace this value ** 
        % ** ANSWER **x
        num_pad = size(kernel, 1) - 1;
        
        padded_signal = [zeros(num_pad, 1); signal; zeros(num_pad, 1)]; 
        output = conv(padded_signal,kernel, 'valid');
    elseif strcmp(pad_mode, 'same')
        % NOTE: You only need to replicate the built-in function for 'same'
        % in the case that the length of the kernel is odd. (In the case 
        % of even kernel size, we need to pad asymetrically in order to 
        % replicate, for those curious).
        % ** YOUR CODE HERE ** 
        
        % ** ANSWER ** 
        num_pad = floor((size(kernel, 1) - 1)/2);
        
        padded_signal = [zeros(num_pad, 1); signal; zeros(num_pad, 1)]; 
        output = conv(padded_signal,kernel, 'valid');    
    else
        error('unsupported pad mode');
    end
     
end


