path(path,'F:\fMRI_Data_Analysis\fMRIstatCBA') %path to the fMRIstatCBA software

cd 'F:\fMRI_Data_Analysis\Data' %data directory


%ERRATA: THE RESULTS REPORTED IN THE PAPER ARE ONLY APPROXIMATE DUE TO A BUG THAT HAS BEEN NOW CORRECTED.

%-------------preprocessing step: DO ONLY ONCE! takes few hours-------%

%The processed files are stored in analyze format in 'rbrains_mcf_stc' sub
%directory



% clear;
% rbrain = brain_preprocess('dat+10+cbi_bold_localizer_mcf_stc.img', 'rbrains_mcf_stc\rbrain10');
% clear;
% rbrain = brain_preprocess('dat+09+cbi_bold_localizer_mcf_stc.img', 'rbrains_mcf_stc\rbrain09');
% clear;
% rbrain = brain_preprocess('dat+12+cbi_bold_localizer_mcf_stc.img', 'rbrains_mcf_stc\rbrain12');
% clear;
% rbrain = brain_preprocess('dat+11+cbi_bold_localizer_mcf_stc.img', 'rbrains_mcf_stc\rbrain11');
% clear;

%-----------------brain mask-------------------------------------%
mask_files = { 'dat+10+cbi_bold_localizer_mcf_stc.img';
    'dat+09+cbi_bold_localizer_mcf_stc.img'; 'dat+12+cbi_bold_localizer_mcf_stc.img';
    'dat+11+cbi_bold_localizer_mcf_stc.img'; 'dat+13+cbi_bold_localizer_mcf_stc.img'}

isbrain = brain_thresh(mask_files,16); %3D matrix with 16 slices, 1 for brain, 0 for background. 

%--------------brain cluster------------------------------%
input_files ={'rbrains_mcf_stc\rbrain09.img'; 'rbrains_mcf_stc\rbrain12.img'}
 keep = [7:166 183:342]; %the time points to include in clustering. 
z=brain_cluster(input_files, isbrain,keep); % 5.3862    2.0000   47.0000    4.0000    4.7859


%--------------LOCALIZER design matrix: using FMRISTAT function fmridesign-------------%
%read ER activation times
fid = fopen('dat_LOlocER_r3.txt', 'r');
ind=1;
while 1
    tmp = fgetl(fid);
    if ~isnumeric(tmp) & ~isempty(tmp)
        x = str2num(tmp);
        data(ind,:) = x;
        ind = ind+1;
    else
        break;
    end
end
events = [data  (6:2:165)'];
events09 = events(data>0,:);

fid = fopen('dat_LOlocER_r6.txt', 'r');
ind=1;
while 1
    tmp = fgetl(fid);
    if ~isnumeric(tmp) & ~isempty(tmp)
        x = str2num(tmp);
        data(ind,:) = x;
        ind = ind+1;
    else
        break;
    end
end
events = [data  (6:2:165)'];
events12 = events(data>0,:);
events12(:,2) = events12(:,2)+176;
events = [events09; events12]

%design matrix
frametimes=(0:176*2-1);
slicetimes = zeros(1,16);
 X_cache=fmridesign(frametimes,slicetimes,events);
 
 keep = [7:166 183:342]; %the time points to analyze;
 X=[squeeze(X_cache.X(keep,:,1,:))];

 X(:,3,:)=ones(size(X,1),size(X,3));

  %-----------------Find ROI---------------%
 q = 0.0001; %FDR threshold
 contrast = [1 -1 0];
 input_files = {'rbrains_mcf_stc\rbrain09.img'; 'rbrains_mcf_stc\rbrain12.img'}; %localizer experiments. 
 isROI = brain_LM(input_files,  z, z>0,X,keep,contrast,q,0); %232 active voxels
 temp = z;
temp(isROI==0)=0;
for i=1:size(temp,3)
    figure;
    imagesc(squeeze(temp(:,:,i)))
end

 %--------------- IC, noIC eventrelated expt Design matrix--------------
%read ER activation times
fid = fopen('dat_ICSR_ER_1bk_r5_m3.txt', 'r');
ind=1;
while 1
    tmp = fgetl(fid);
    if ~isnumeric(tmp) & ~isempty(tmp)
        x = str2num(tmp);
        data(ind,:) = x;
        ind = ind+1;
    else
        break;
    end
end
events = [data  (6:2:253)'];
events = events(data>0,:);

%design matrix
frametimes=(0:263);
slicetimes = zeros(1,16);
 X_cache=fmridesign(frametimes,slicetimes,events);
 %find the correlation matrix
 keep = 7:254;
 X=[squeeze(X_cache.X(keep,:,1,:))];
X(:,5,:)=ones(size(X,1),size(X,3));
 %-----------------Find Activations in IC, noIC experiment within ROI---------------%
 q = 0.05; %FDR threshold
 contrastwIntercept = [1 -1 0 0 0];
 input_files = {'rbrains_mcf_stc\rbrain11.img'} 
 isActive = brain_LM(input_files,  z, isROI,X,keep,contrastwIntercept,q,1); %44 active  voxels, pA=0.0051

 save('F:\fMRI_Data_Analysis\savedir\isActiveER05','isActive')
 zROI = zeros(size(z));
 zROI(isROI>0) = z(isROI>0);
 u=unique(zROI(zROI>0));
 newindex = 1:20;
 newzROI= zROI;
 for i=1:length(newindex)
     newzROI(zROI==u(i))= newindex(i);
 end
save('F:\fMRI_Data_Analysis\savedir\zROI','newzROI')
 
 temp = z;
temp(isActive==0)=0;
for i=10:11
    figure;
    imagesc(squeeze(temp(:,:,i)))
end

 
 %--For comparison, find activations using VBA--

%every brain voxel gets an index
zv = zeros(size(isbrain));
k=1;
for h=1:size(isbrain,3)
    for i=1:size(isbrain,1)
        for j=1:size(isbrain,2)
            if isbrain(i,j,h)>0
                zv(i,j,h)=k;
                k=k+1;
            end
        end
    end
end

  isActiveV = brain_LM(input_files,  zv, isROI,X,keep,contrastwIntercept,q,1);%0 active voxels
  
   %--------------IC, noIC block design experiment Design matrix-------------%
  events = [[1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 ] ;6:16:(16*16+5); 16*ones(1,16)]';
%design matrix
frametimes=(0:135)*2;
slicetimes = zeros(1,16);
 X_cache=fmridesign(frametimes,slicetimes,events);
 %find the correlation matrix
 keep = 4:131;
 X=[squeeze(X_cache.X(keep,:,1,:))];
 X(:,3,:)=ones(size(X,1),size(X,3));

 
 %-------------Find Activations in IC, noIC block design experiment within ROI---------------%
 q = 0.1; %FDR threshold

 contrastwIntercept = [1 -1 0];
  input_files = {'rbrains_mcf_stc\rbrain10.img'} 
 isActive = brain_LM(input_files,  z, isROI,X,keep,contrastwIntercept,q,1); %174 active voxels
save('F:\fMRI_Data_Analysis\savedir\isActiveBLA10','isActive')
 

  temp = z;
temp(isActive==0)=0;
for i=10:11
    figure;
    imagesc(squeeze(temp(:,:,i)))
end

 isActive = brain_LM(input_files,  z, isROI,X,keep,contrastwIntercept,q,0); %122 active voxels
save('F:\fMRI_Data_Analysis\savedir\isActiveBL10','isActive')

  temp = z;
temp(isActive==0)=0;
for i=10:11
    figure;
    imagesc(squeeze(temp(:,:,i)))
end

 %--For comparison, find activations using VBA--
 isActiveV = brain_LM(input_files,  zv, isROI,X,keep,contrast,q,1);%84 active voxels
save('F:\fMRI_Data_Analysis\savedir\isActiveBLA10','isActive')

 temp = z;
temp(isActiveV==0)=0;
for i=10:11
    figure;
    imagesc(squeeze(temp(:,:,i)))
end
 