disp('not meant to be executed! Evaluate commands line by line');
%return

%SFN1 - IIn this exercise we will compress/denoise an image sequence of
%neural activity in a brain slice.

clear all;

%load data time series already in a matlab array
%note - importing movies requires a few lines of code. See loadtiff.m (provided) for
%an example.

load freerun.mat;


%Variables loaded from this array include:
%stack16 - movie frames in 16bit unsigned integer format
%frames - number of images in the time series
%height - number of pixels in the y direction
%width - number of pixels in the x direction
%sampfreq - frame rate of the camera, in hz.

%It is important to keep in mind the format of your data, hence the
%variable name "stack16."

%movie display kludge: 
dispframes=100;
delay=0;
dispsize=3;
dispstack(stack16,dispframes,dispsize,delay);

%make a time series of the change in fluorescence relative to the first
%frame:
dffo_stackDB=dffoSFN(stack16,dispframes,1,0);

% and display it:
dispstack(dffo_stackDB,dispframes,dispsize,delay);

%in order to process the movie we must convert it to vectors, as discussed.
vectordataDB=reshape(double(stack16), width*height, frames);

%calculate the SVD
% [evecs,evals,coeffs]=svd(vectordataDB,0);
% evals=diag(evals);

% OR just load it if your laptop is challenged:
load 'freerunSVD.mat';

%we don't need the original data anymore
clear stack16 vectordataDB;

%Plot the first 4 spatial modes and associated coefficients
doPCAPlot;

%Plot the next 6 spatial modes and associated coefficients
doPCAPlot2;

%let's look more carefully at the singular values:
figure;
semilogy(evals(2:921));
figure;
semilogy(evals(2:100));

%choose the number of singular values for the reconstructed movie:
recon_keep=10;

%reconstruct the movie by adding up these modes, singular values and
%coefficients:
recondataDB = evecs(:,1:recon_keep) * diag(evals(1:recon_keep)) * coeffs(:,1:recon_keep)';

%The data is still in vector format; reshape it back to movie frames
recon_stackDB=reshape(recondataDB,height,width,frames);

%clear the raw reconstruction
clear recondataDB

% look at it if you want to
%dispstack(recon_stackDB,dispframes,dispsize,delay);

%make a time series of the change in fluorescence relative to the first
%frame:
dffo_reconDB=dffoSFN(recon_stackDB,dispframes,1,0);

% and display it:
dispstack(dffo_reconDB,dispframes,dispsize,delay);



