%In this lab, we will process an LFP trace in the time domain to illustrate
%1D convolution and an image to illustrate 2D convolution 

%% 1 Init - Clear workspace and load the data
clear all
close all
clc 
load('LFPdata.mat'); 
%After loading, the workspace will contain a variable "A" that contains
%10,000 voltage samples per second, but only for 300 milliseconds. They
%represent the "local field potential", voltages in microvolt at the
%electrode tip. 

%% 2a Exploratory data analysis - looking at the raw data
figure
plot(A)
shg
%By visual inspection, we note the presence of a signal, but also of a lot
%of "noise". The "sharp stuff" in particular does not look very
%physiological. There are no known processes in the brain that could
%generate such signals. But there are a lot of known problems with
%electrical noise that could. Solution: Smoothing by convolution.
%Carmen's problem: How big should we make the kernel? 

%% 2b Convolution. If you have domain knowledge, you can pick a suitable kernel size. But we
%are first years, so we don't. So, we are going to try a lot of them, then
%pick the one that looks best. 
for ii = 1:length(A) %Shelby's gambit - go all the way
    kernelSize = ii; %This is the length of any given kernel
    kernel = ones(kernelSize,1); %Symmetric, equal weights, the simplest one
    pause(0.3)
    Aconv = conv(A,kernel,'full')./sum(kernel); %If we don't normalize, the output values scale with kernel length
    plot(Aconv)
    title(['Kernel size = ', num2str(ii), ' samples'])
    ylim([-40 20])
    xlim([-100 3100])
    shg
end

%% Moving on to 2D convolution. In the image domain. 
%Before we move on to actual images, let's do that in a case where we know
%"ground truth". This is the reason some physiologists use artificial
%images. In praise of artifice. So you can see what convolution is doing. 

%% 1. Let's create luminance ramps - introducing meshgrid
numSteps = 7; %We want this many luminance steps
minLum = 0; %Lowest luminance value
maxLum = 255; %Highest luminance value in an 8 bit system
temp = linspace(minLum, maxLum, numSteps); %This creates them, we want to know the stepsize
stepSize = unique(round(diff(temp),4)); %Determine the stepsize from that
%To create an artificial 2d image, it is best to introduce the function
%"meshgrid" at this point. It has 2 outputs, x and y, which are a full
%cross of the inputs, keeping one dimension constant while incrementing the
%other one and vice versa. 
[x,y] = meshgrid(minLum:stepSize:maxLum,minLum:stepSize:maxLum); %This will create luminance ramps
figure
subplot(1,2,1)
imagesc(x)
subplot(1,2,2)
imagesc(y)
colorbar
colormap(bone)
shg

%% 2 Let's use this to create an artificial image
%Matlab represents images as 3D matrices. The first 2 dimensions are x and
%y coordinates, the third dimension are R,G and B sheets. The values
%represent luminances in 8 bit, so they range from 0 to 255.
clear IM %This will be a 3D matrix that represents the image. It will have 3 sheets, we create it one by one
IM(:,:,1) = uint8(x); %Red channel - the x's
IM(:,:,2) = uint8(zeros(length(x)));  %Silence the green channel to get purple
IM(:,:,3) = uint8(y); %Uint8 interprets/renders the values in 8 bit, as an unsigned integer.

figure
subplot(1,2,1)
image(IM)
axis square
title('Unfiltered')

kernel = ones(5); %Small summing kernel 
kernelWeight = sum(kernel(:)); %Sum of the flattened kernel. We need to sum in both dimensions
convIM = uint8(round(convn(IM, kernel, 'valid')) ./kernelWeight); %If we don't normalize by the kernelWeight, the outputs will be - in this case - 25 times too bright
subplot(1,2,2)
image(convIM)
axis square
title('Convolved')


%% 4 Let's use an actual image - writing the dress code
dress = imread('thedress.jpg'); %Putting the dress into a matrix
figure
image(dress)
axis equal %Fix the aspect ratio to unstretch the dress. To make an equal

%unit size

%% 5 Let's convolve the dress with a summing kernel = Blurring
kernelWidth = 50; 
figure
subplot(1,3,1)
image(dress)
title('Original')
subplot(1,3,2)
kernel = ones(kernelWidth);
kernelWeight = sum(kernel(:)); 
lp = uint8(convn(dress,kernel,'same')./kernelWeight); %Same makes sense for image processing
image(lp)
title('Blurry')
%% 6 Contrast enhancement with a difference filter = Sharpening
kernel = [];
kernel(1,:) = [0 -1 0]; 
kernel(2,:) = [-1 5 -1]; 
kernel(3,:) = [0 -1 0]; 
subplot(1,3,3)
kernelWeight = sum(kernel(:)); 
sharper = uint8(convn(dress,kernel,'same')./kernelWeight);
image(sharper)
title('The sharper image')

%% 7 Implementing reflectance as an option for non-zero padding convolution
%Define the pieces
upsideDown = flipud(dress); %This goes on top and bottom
leftRightFlip = fliplr(dress); %This goes left and right
both = flipud(fliplr(dress)); %This goes 

%Piece it together
biggerDress = [both upsideDown both;leftRightFlip dress leftRightFlip; both upsideDown both];

%Check that it worked
figure
image(biggerDress)
pause
kernelWidth = 50; 

%Do the convolution - this will take a while
figure
kernel = ones(kernelWidth);
kernelWeight = sum(kernel(:)); 
convBiggerDress = uint8(convn(biggerDress,kernel,'full')./kernelWeight);
image(convBiggerDress)
pause

%Cut to size (find the right indices) - we might be off by one
dim1 = size(convBiggerDress,1)/3;
dim2 = size(convBiggerDress,2)/3; 
image(convBiggerDress(dim1:dim1*2,dim2:dim2*2,:))
axis equal
