%% 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);
% Number of folds for which to run k-means clustering
n_folds = 1;
% 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;
% Type of distance to use for k-means clustering
DistType = 'cosine';
% 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
disp(['Selecting frames for subject ',num2str(i),'...']);
% I has size n_timepoints x n_voxels and contains the innovations for
% the considered subject
I = Innovation{i}';
%%% Q14. Check how many voxels show positive innovation at each time
%%% point (saved in POS, size n_timepoints x n_voxels)
POS = mask1{i}';
POS(POS < 0) = 0;
%%% Q15. Creates the positive innovation matrix, of size n_timepoints x
%%% n_voxels. It contains only the positive values that survive
%%% thresholding
I_pos = I;
I_pos(I < 0) = 0;
I_pos = (I_pos .* POS);
%%% Q16. Removing isolated positive elements from the data
I_pos = check_interconnectedness(I_pos,fHeader,mask);
%%% Q17. Creates a mask denoting what frames must be retained as moments
%%% of positive innovation
mask_POS = zeros(n_timepoints,1);
mask_POS(sum(POS,2) > n_voxels) = 1;
mask_POS = logical(mask_POS);
%%% Q18. Concatenates the frames in question for clustering
Data = [Data;I_pos(mask_POS,:)];
%%% Q19. The same process is run for negative frames, which are inverted
NEG = mask1{i}';
NEG(NEG > 0) = 0;
NEG = abs(NEG);
I_neg = I;
I_neg(I > 0) = 0;
I_neg = -1*(I_neg .* NEG);
I_neg = check_interconnectedness(I_neg,fHeader,mask);
mask_NEG = zeros(n_timepoints,1);
mask_NEG(sum(NEG,2) > n_voxels) = 1;
mask_NEG = logical(mask_NEG);
Data = [Data;I_neg(mask_NEG,:)];
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);