%% This is a tutorial script for the generation of innovation-driven co-
% activation patterns (iCAPs) from total activation outputs
%
% Created for BrainHack Zurich (March 3rd 2017)
%
% Written by Younes Farouj and Thomas Bolton, Medical Image Processing
% Laboratory, EPFL, Switzerland
addpath(genpath('Functions'));
addpath(genpath('StartData'));
addpath(genpath('Utilities'));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% 1. Data loading and parameters %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%% Data loading
% Cell array: each cell is n_GM_voxels x n_timepoints. For each location and
% time point, the mask contains +1 for a positive innovation, -1 for a
% negative innovation, and 0 for no innovation
load('mask1.mat')
% Innovation signal: each cell is n_GM_voxels x n_timepoints.
load('Innovation.mat')
% Logical vector of size n_voxels x 1 indexing the elements belonging to
% gray matter with +1. Used for plotting the iCAPs.
load('mask.mat');
% Header containing information about the considered data (dimensions,
% etc.). Used for plotting the iCAPs
load('fHeader.mat');
%%% Parameters not to change
% Number of subjects in the dataset
n_subjects = length(Innovation);
% Number of time points per subject
n_timepoints = size(Innovation{1},2);
% Number of gray matter voxels
n_GM_voxels = size(mask1{1},1);
% Type of distance to use for k-means clustering
DistType = 'cosine';
% Number of folds for which to run k-means clustering
n_folds = 3;
% Number of slices to display in each row of the iCAP plots
n_slices_x = 6;
% Ranges for the positive and negative values to display in iCAPs plotting
Range_Pos = [1.5,4];
Range_Neg = [-1.5,-4];
%%% Parameters that can be changed
% Number of voxels that must show innovation together to retain a given
% frame for the clustering
n_voxels = 500;
% Number of clusters to separate the data into
K = 20;
% Type of slices to display for the iCAP plots: choose between 'axial',
% 'sagittal' and 'coronal'
typeslice = 'axial';
% Associated range of (MNI) values at which to draw the slices. Suggested
% values: for axial slices, -30:10:80; for coronal slices, -100:10:70; for
% sagittal slices, -70:10:70
slice_range = -30:10:80;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% 2. Frame selection process %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% Data will contain our inputs for the k-means clustering process
Data = [];
% Looping through subjects...
for i=1:n_subjects
% I has size n_timepoints x n_voxels and contains the innovations for
% the considered subject
I = Innovation{i}';
%%%%%%% TO FILL %%%%%%%%%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% 3. Generation of iCAPs %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% K-means clustering to separate the data into K different iCAPs
[IDX,iCAPs] = kmeans(Data,K,'Distance',DistType,...
'Replicates',n_folds);
% Z-scoring of the iCAPs
iCAPs_z = ZScore_iCAPs(iCAPs,Data,IDX);
% Rearranges the data and centroids and indices so that '1' denotes
% the cluster with most occurrences in the data
[Data_rear,IDX_rear,iCAPs_rear] = RearrangeTimeCoursesClassic(Data,...
IDX,iCAPs_z);
% Plotting of all K iCAPs
PlotiCAPs(which('MNI152_T1_2mm_brain.nii'),iCAPs_rear(1:K,:),fHeader,...
mask,Range_Pos,Range_Neg,n_slices_x,slice_range,typeslice);