%Choosing working directory
cd ~/surfstat

clear all; close all; clc

surf = SurfStatReadSurf( 'av.obj' );

%Adding the midline mask
mask = SurfStatMaskCut( surf );

%Adding the brainstem mask
maskb = (mask & SurfStatROI([0; -16; -8], 20, surf ) == 0);

cd ~/surfstat/glimfiles

%Creating the matrix of data.
[SubjectNum, n3, ThickFileLeft, ThickFileRight]=textread('glimfile.csv', ' %s %s %s %s');

%To read and assign the actual thickness data
T1 = SurfStatReadData( [ThickFileLeft] );
T2 = SurfStatReadData( [ThickFileRight] );
T3 = horzcat(T1, T2);

%Create terms
N3 = term (n3);

%% 1 You are now ready to create your model and to estimate it

Y = 1 + N3 + random(SubjectNum) + I;
figure(1)
image(Y);
slm = SurfStatLinMod( T3, Y, surf );

% MAIN EFFECT OF N3, DIRECTION 1 (N3=50 > N3=100)

%To get your t statistic for group
contrast_direction1 = N3.50 - N3.100;
slm_direction1 = SurfStatT ( slm, contrast_direction1);
figure(2)
SurfStatView ( slm_direction1.t.*maskb, surf, 'tmap N3=50 > N3=100' );
set(gcf,'numbertitle','off','name', 'tmap N3=50 > N3=100')
%SurfStatColLim([-12,12]);

%To get thresholded p values using Random Field Theory
[ pval, peak, clus ] = SurfStatP( slm_direction1, maskb );
figure(3)
SurfStatView( pval, surf, 'RFT map N3=50 > N3=100');
set(gcf,'numbertitle','off','name', 'RFT map N3=50 > N3=100')

% MAIN EFFECT OF N3, DIRECTION 2 (N3=100 > N3=50)

%To get your t statistic for group
contrast_direction2 = N3.100 - N3.50;
slm_direction2 = SurfStatT ( slm, contrast_direction2);
figure(4)
SurfStatView ( slm_direction2.t.*maskb, surf, 'tmap N3=100 > N3=50' );
set(gcf,'numbertitle','off','name', 'tmap N3=100 > N3=50')
%SurfStatColLim([-12,12]);

%To get thresholded p values using Random Field Theory
[ pval, peak, clus ] = SurfStatP( slm_direction2, maskb );
figure(5)
SurfStatView( pval, surf, 'RFT map N3=100 > N3=50');
set(gcf,'numbertitle','off','name', 'RFT map N3=100 > N3=50')

%Mean cortical thickness by N3
figure(6)
SurfStatView( mean(T3(1:5,:)), surf, 'mean cortical thickness N3=50');
set(gcf,'numbertitle','off','name','mean cortical thickness N3=50')
SurfStatColLim([0,5.5]);

figure(7)
SurfStatView( mean(T3(6:10,:)), surf, 'mean cortical thickness N3=100');
set(gcf,'numbertitle','off','name','mean cortical thickness N3=100')
SurfStatColLim([0,5.5]);

%Standard deviation by N3
figure(8)
SurfStatView( std(T3(1:5,:)), surf, 'std cortical thickness N3=50');
set(gcf,'numbertitle','off','name','std cortical thickness N3=50')
SurfStatColLim([0,.7]);

figure(9)
SurfStatView( std(T3(6:10,:)), surf, 'std cortical thickness N3=100');
set(gcf,'numbertitle','off','name','std cortical thickness N3=100')
SurfStatColLim([0,.7]);