NHP-pRF/Analysis_code/MRI/pRF_prepdata_avg.m

function pRF_prepdata_avg(SessionList,doUpsample)
% collects data, concatenates, downsamples stimulus and resaves
% NB: won't run on macbook unless there's a USB-drive with all data 

% Extra regressor are not possible because we're averaging the BOLD here
% z-score per voxel and average over runs in a session

%% WHICH DATA =============================================================
%clear all; clc;
if nargin <2
    fprintf('Not enough arguments specified, will use defaults:\n');
    fprintf('SessionList: pRF_PrepDatalist_M01\n');
    pRF_PrepDatalist_M01;
    %pRF_PrepDatalist_M02;
    doUpsample = true;
else
    eval(SessionList);
end

%% Sweep to volume mapping ------------------------------------------------
TR = 2.5;

SwVolMap_230 = { ... % 15 blanks / 20 steps
    1 , 6:25    ;...
    2 , 41:60   ;...
    3 , 61:80   ;...
    4 , 96:115  ;...
    5 , 116:135 ;...
    6 , 151:170 ;...
    7 , 171:190 ;...
    8 , 206:225 };

SwVolMap_210 = { ... % also for 215 vols
    1 , 6:25    ;... % 10 blanks / 20 steps
    2 , 36:55   ;...
    3 , 56:75   ;...
    4 , 86:105  ;...
    5 , 106:125 ;...
    6 , 136:155 ;...
    7 , 156:175 ;...
    8 , 186:205 };

SwVolMap_218 = { ... % 12 blanks / 20 steps
    1 , 6:25    ;...
    2 , 38:57   ;...
    3 , 58:77   ;...
    4 , 90:109  ;...
    5 , 110:129 ;...
    6 , 142:161 ;...
    7 , 162:181 ;...
    8 , 194:213 };

SwVolMap_436 = { ...
    1 , 6:25    ;...
    2 , 38:57   ;...
    3 , 58:77   ;...
    4 , 90:109  ;...
    5 , 110:129 ;...
    6 , 142:161 ;...
    7 , 162:182 ;...
    8 , 194:213 ;...
    1 , 224:243 ;...
    2 , 256:275 ;...
    3 , 276:295 ;...
    4 , 308:327 ;...
    5 , 328:347 ;...
    6 , 360:379 ;...
    7 , 380:399 ;...
    8 , 412:431 };

Is436 = false;

%% INITIALIZE =============================================================

homefld = pwd;
cd ../../
SHARED_ROOT_FLD = pwd;
cd(homefld)

% Platform specific basepath
tool_basepath = fullfile(SHARED_ROOT_FLD,'Toolboxes');

% Link to the brain mask
BrainMask_file = fullfile(SHARED_ROOT_FLD,'Preprocessed_data','manual-masks',...
    ['sub-' MONKEY],'func',['sub-' MONKEY '_ref_func_mask_res-1x1x1.nii']);

% create a folder to save outputs in
if doUpsample
    out_folder = fullfile(SHARED_ROOT_FLD,'Raw_data','MRI','Raw_matlab',...
        'us-padded', MONKEY);
else
    out_folder = fullfile(SHARED_ROOT_FLD,'Raw_data','MRI','Raw_matlab',...
        'padded', MONKEY);
end

warning off %#ok<*WNOFF>
[~,~,~] = mkdir(out_folder);
warning on %#ok<*WNON>

%% GET THE FILE-PATHS OF THE IMAGING  & STIM-MASK FILES ===================
% All functional runs that are preprocessed with the BIDS pipeline are
% resampled to 1x1x1 mm isotropic voxels, reoriented from sphinx,
% motion corrected, (potentially smoothed with 2 mm FWHM), and
% registered to an example functional volume
% (so they're already in a common space)
% do the analysis in this functional space than we can register to hi-res
% anatomical data and/or the NMT template later
sessions = unique(DATA(:,1)); %#ok<*NODEF>

% These BIDS-basepath files are not included because they are too large
BIDS_basepath = fullfile(SHARED_ROOT_FLD,'Raw_data','MRI','Raw_nifti');

monkey_path_nii = fullfile(BIDS_basepath, 'derivatives',...
    'featpreproc','highpassed_files',['sub-' MONKEY]);
monkey_path_stim = fullfile(BIDS_basepath,['sub-' MONKEY]);
monkey_path_motion.regress = fullfile(BIDS_basepath, 'derivatives',...
    'featpreproc','motion_corrected',['sub-' MONKEY]);
monkey_path_motion.outlier = fullfile(BIDS_basepath, 'derivatives',...
    'featpreproc','motion_outliers',['sub-' MONKEY]);

for s=1:length(sessions)    
    sess_path_nii{s} = fullfile(monkey_path_nii, ['ses-' sessions{s}(1:8)], 'func'); %#ok<*SAGROW>
    sess_path_stim{s} = fullfile(monkey_path_stim, ['ses-' sessions{s}(1:8)], 'func');
    sess_path_motreg{s} = fullfile(monkey_path_motion.regress, ['ses-' sessions{s}(1:8)], 'func');
    sess_path_motout{s} = fullfile(monkey_path_motion.outlier, ['ses-' sessions{s}(1:8)], 'func');
    runs = unique(DATA(strcmp(DATA(:,1),sessions{s}),2));
    for r=1:length(runs)
        if ispc % the ls command works differently in windows
            a=ls( fullfile(sess_path_nii{s},['*run-' runs{r} '*.nii.gz']));
            run_path_nii{s,r} = fullfile(sess_path_nii{s},a(1:end-3));
            b = ls( fullfile(sess_path_stim{s}, ...
                ['*run-' runs{r} '*model*']));
            run_path_stim{s,r}= fullfile(sess_path_stim{s},b,'StimMask.mat');
            c=ls( fullfile(sess_path_motreg{s},['*run-' runs{r} '*.param.1D']));
            run_path_motreg{s,r} = fullfile(sess_path_motreg{s},c);
            d=ls( fullfile(sess_path_motout{s},['*run-' runs{r} '*.outliers.txt']));
            run_path_motout{s,r} = fullfile(sess_path_motout{s},d);
            e = ls( fullfile(sess_path_stim{s}, ...
                ['*run-' runs{r} '*model*']));
            run_path_rew{s,r}= fullfile(sess_path_stim{s},e,'RewardEvents.txt');
        else
            a = ls( fullfile(sess_path_nii{s},['*run-' runs{r} '*.nii.gz']));
            run_path_nii{s,r} = a(1:end-3);
            run_path_nii{s,r} = run_path_nii{s,r}(1:end-1);
            run_path_stim{s,r} = ls( fullfile(sess_path_stim{s}, ...
                ['*run-' runs{r} '*model*'],'StimMask.mat'));
            run_path_stim{s,r} = run_path_stim{s,r}(1:end-1);
            run_path_motreg{s,r} = ls( fullfile(sess_path_motreg{s}, ...
                ['*run-' runs{r} '*.param.1D']));
            run_path_motreg{s,r}=run_path_motreg{s,r}(1:end-1);
            run_path_motout{s,r} = ls( fullfile(sess_path_motout{s}, ...
                ['*run-' runs{r} '*_outliers.txt']));
            run_path_motout{s,r}=run_path_motout{s,r}(1:end-1);
            run_path_rew{s,r} = ls( fullfile(sess_path_stim{s}, ...
                ['*run-' runs{r} '*model*'],'RewardEvents.txt'));
            run_path_rew{s,r}=run_path_rew{s,r}(1:end-1);
        end
        sweepinc{s,r} = DATA( ...
            (strcmp(DATA(:,1),sessions{s}) & strcmp(DATA(:,2),runs{r})),3);
    end
end

%% LOAD & RE-SAVE STIMULUS MASKS & NIFTI ==================================
for s=1:size(run_path_stim,1) % sessions
    fprintf(['Processing session ' sessions{s} '\n']);
    rps = [];
    if strcmp(sessions{s},'20160721')
        TR=TR_3;
        fprintf('!!!! NB: This TR is 3s do not use in averaging !!!!\n');
        % should be excluded in the prepdata list
    end

    for i=1:size(run_path_stim,2)
        if ~isempty(run_path_stim{s,i})
            rps=[rps i];
        end
    end
    
    for r=rps % runs
        % stimulus mask -----
        % loads variable called stimulus (x,y,t) in volumes
        load(run_path_stim{s,r}(1:end-4));
        
        sinc = cell2mat(sweepinc{s,r});
        if size(stimulus,3) == 210 || size(stimulus,3) == 215
            SwVolMap = SwVolMap_210;
            Is436 = false;
        elseif size(stimulus,3) == 218
            SwVolMap = SwVolMap_218;
            Is436 = false;
        elseif size(stimulus,3) == 230
            SwVolMap = SwVolMap_230;
            Is436 = false;
        elseif size(stimulus,3) == 436
            SwVolMap = SwVolMap_436;  
            Is436 = true;
        else
            error('weird number of stimulus frames');
        end
        firstvol = SwVolMap{min(sinc),2}(1) - 5;
        lastvol = SwVolMap{max(sinc),2}(end) + 5;
        v_inc=firstvol:lastvol;
        v_all=1:size(stimulus,3);
        
        bin_vinc = zeros(1,size(stimulus,3));
        bin_vinc(v_inc) = 1;
                
        % volumes ------
        %fprintf('Unpacking nii.gz');
        %uz_nii=gunzip(run_path_nii{s,r});
        % >>> Unpacking is slow from within matlab. Do this in the system
        
        temp_nii=load_nii(run_path_nii{s,r});%load_nii(uz_nii{1});
        %delete(uz_nii{1});
        fprintf(' ...done\n');
              
        % save the session-based stims & vols -----
        for v=1:size(stimulus,3)
            % resample image (160x160 pix gives 10 pix/deg)
            s_run(r).stim{v} = imresize(stimulus(:,:,v_all(v)),[160 160]);
            s_run(r).vol{v} = temp_nii.img(:,:,:,v_all(v));
            if v==1
                s_run(r).hdr = temp_nii.hdr;
            end
        end
        
        % create nan volume/stim to pad with
        nanVol=nan(size(temp_nii.img(:,:,:,1)));
        nanStim=nan(160);
        
        % pad everything to 230 volumes
        p_run(r).stim = cell(1,230);
        p_run(r).vol = cell(1,230);
        p_run(r).hdr = s_run(r).hdr;
        p_run(r).inc = zeros(1,230);
        if Is436
            shiftruns=max(rps);
            p_run(r+shiftruns).stim = cell(1,230);
            p_run(r+shiftruns).vol = cell(1,230);
            p_run(r+shiftruns).hdr = s_run(r).hdr;
            p_run(r+shiftruns).inc = zeros(1,230);
        end
       
        if size(stimulus,3) <= 215 % 210/215
            p_run(r).stim = [ ...
                s_run(r).stim(1:SwVolMap{2,2}(1)-1) ...
                nanStim nanStim nanStim nanStim nanStim ....
                s_run(r).stim(SwVolMap{2,2}(1):SwVolMap{4,2}(1)-1) ...
                nanStim nanStim nanStim nanStim nanStim ....
                s_run(r).stim(SwVolMap{4,2}(1):SwVolMap{6,2}(1)-1) ...
                nanStim nanStim nanStim nanStim nanStim ....
                s_run(r).stim(SwVolMap{6,2}(1):SwVolMap{8,2}(1)-1) ...
                nanStim nanStim nanStim nanStim nanStim ....
                s_run(r).stim(SwVolMap{8,2}(1):SwVolMap{8,2}(end)+5) ...
                ];
            p_run(r).inc = [ ...
                bin_vinc(1:SwVolMap{2,2}(1)-1) ...
                0 0 0 0 0 ....
                bin_vinc(SwVolMap{2,2}(1):SwVolMap{4,2}(1)-1) ...
                0 0 0 0 0 ....
                bin_vinc(SwVolMap{4,2}(1):SwVolMap{6,2}(1)-1) ...
                0 0 0 0 0 ....
                bin_vinc(SwVolMap{6,2}(1):SwVolMap{8,2}(1)-1) ...
                0 0 0 0 0 ....
                bin_vinc(SwVolMap{8,2}(1):SwVolMap{8,2}(end)+5) ...
                ];
            p_run(r).vol = [ ...
                s_run(r).vol(1:SwVolMap{2,2}(1)-1) ...
                nanVol nanVol nanVol nanVol nanVol ....
                s_run(r).vol(SwVolMap{2,2}(1):SwVolMap{4,2}(1)-1) ...
                nanVol nanVol nanVol nanVol nanVol ....
                s_run(r).vol(SwVolMap{4,2}(1):SwVolMap{6,2}(1)-1) ...
                nanVol nanVol nanVol nanVol nanVol ....
                s_run(r).vol(SwVolMap{6,2}(1):SwVolMap{8,2}(1)-1) ...
                nanVol nanVol nanVol nanVol nanVol ....
                s_run(r).vol(SwVolMap{8,2}(1):SwVolMap{8,2}(end)+5) ...
                ];
            % add 5 blanks
        elseif size(stimulus,3) == 218
           p_run(r).stim = [ ...
                s_run(r).stim(1:SwVolMap{2,2}(1)-1) ...
                nanStim nanStim nanStim ....
                s_run(r).stim(SwVolMap{2,2}(1):SwVolMap{4,2}(1)-1) ...
                nanStim nanStim nanStim ...
                s_run(r).stim(SwVolMap{4,2}(1):SwVolMap{6,2}(1)-1) ...
                nanStim nanStim nanStim ....
                s_run(r).stim(SwVolMap{6,2}(1):SwVolMap{8,2}(1)-1) ...
                nanStim nanStim nanStim ....
                s_run(r).stim(SwVolMap{8,2}(1):SwVolMap{8,2}(end)+5) ...
                ];
            p_run(r).inc = [ ...
                bin_vinc(1:SwVolMap{2,2}(1)-1) ...
                0 0 0 ....
                bin_vinc(SwVolMap{2,2}(1):SwVolMap{4,2}(1)-1) ...
                0 0 0 ....
                bin_vinc(SwVolMap{4,2}(1):SwVolMap{6,2}(1)-1) ...
                0 0 0 ....
                bin_vinc(SwVolMap{6,2}(1):SwVolMap{8,2}(1)-1) ...
                0 0 0 ....
                bin_vinc(SwVolMap{8,2}(1):SwVolMap{8,2}(end)+5) ...
                ];
            p_run(r).vol = [ ...
                s_run(r).vol(1:SwVolMap{2,2}(1)-1) ...
                nanVol nanVol nanVol ....
                s_run(r).vol(SwVolMap{2,2}(1):SwVolMap{4,2}(1)-1) ...
                nanVol nanVol nanVol ....
                s_run(r).vol(SwVolMap{4,2}(1):SwVolMap{6,2}(1)-1) ...
                nanVol nanVol nanVol ....
                s_run(r).vol(SwVolMap{6,2}(1):SwVolMap{8,2}(1)-1) ...
                nanVol nanVol nanVol ....
                s_run(r).vol(SwVolMap{8,2}(1):SwVolMap{8,2}(end)+5) ...
                ];
            % add 3 blanks
        elseif size(stimulus,3) == 230
            p_run(r).stim = s_run(r).stim;
            p_run(r).inc = bin_vinc;
            p_run(r).vol = s_run(r).vol;
            % as is
        elseif size(stimulus,3) == 436
            p_run(r).stim = [ ...
                s_run(r).stim(1:SwVolMap{2,2}(1)-1) ...
                nanStim nanStim nanStim ....
                s_run(r).stim(SwVolMap{2,2}(1):SwVolMap{4,2}(1)-1) ...
                nanStim nanStim nanStim ...
                s_run(r).stim(SwVolMap{4,2}(1):SwVolMap{6,2}(1)-1) ...
                nanStim nanStim nanStim ....
                s_run(r).stim(SwVolMap{6,2}(1):SwVolMap{8,2}(1)-1) ...
                nanStim nanStim nanStim ....
                s_run(r).stim(SwVolMap{8,2}(1):SwVolMap{8,2}(end)+5) ...
                ];
            p_run(r).inc = [ ...
                bin_vinc(1:SwVolMap{2,2}(1)-1) ...
                0 0 0 ....
                bin_vinc(SwVolMap{2,2}(1):SwVolMap{4,2}(1)-1) ...
                0 0 0 ....
                bin_vinc(SwVolMap{4,2}(1):SwVolMap{6,2}(1)-1) ...
                0 0 0 ....
                bin_vinc(SwVolMap{6,2}(1):SwVolMap{8,2}(1)-1) ...
                0 0 0 ....
                bin_vinc(SwVolMap{8,2}(1):SwVolMap{8,2}(end)+5) ...
                ];
            p_run(r).vol = [ ...
                s_run(r).vol(1:SwVolMap{2,2}(1)-1) ...
                nanVol nanVol nanVol ....
                s_run(r).vol(SwVolMap{2,2}(1):SwVolMap{4,2}(1)-1) ...
                nanVol nanVol nanVol ....
                s_run(r).vol(SwVolMap{4,2}(1):SwVolMap{6,2}(1)-1) ...
                nanVol nanVol nanVol ....
                s_run(r).vol(SwVolMap{6,2}(1):SwVolMap{8,2}(1)-1) ...
                nanVol nanVol nanVol ....
                s_run(r).vol(SwVolMap{8,2}(1):SwVolMap{8,2}(end)+5) ...
                ];
            
            p_run(r+shiftruns).stim = [ ...
                s_run(r).stim(SwVolMap{9,2}(1)-5:SwVolMap{10,2}(1)-1) ...
                nanStim nanStim nanStim ....
                s_run(r).stim(SwVolMap{10,2}(1):SwVolMap{12,2}(1)-1) ...
                nanStim nanStim nanStim ...
                s_run(r).stim(SwVolMap{12,2}(1):SwVolMap{14,2}(1)-1) ...
                nanStim nanStim nanStim ....
                s_run(r).stim(SwVolMap{14,2}(1):SwVolMap{16,2}(1)-1) ...
                nanStim nanStim nanStim ....
                s_run(r).stim(SwVolMap{16,2}(1):SwVolMap{16,2}(end)+5) ...
                ];
            p_run(r+shiftruns).inc = [ ...
                bin_vinc(1:SwVolMap{2,2}(1)-1) ...
                0 0 0 ....
                bin_vinc(SwVolMap{2,2}(1):SwVolMap{4,2}(1)-1) ...
                0 0 0 ....
                bin_vinc(SwVolMap{4,2}(1):SwVolMap{6,2}(1)-1) ...
                0 0 0 ....
                bin_vinc(SwVolMap{6,2}(1):SwVolMap{8,2}(1)-1) ...
                0 0 0 ....
                bin_vinc(SwVolMap{8,2}(1):SwVolMap{8,2}(end)+5) ...
                ];
            p_run(r+shiftruns).vol = [ ...
                s_run(r).vol(SwVolMap{9,2}(1)-5:SwVolMap{10,2}(1)-1) ...
                nanVol nanVol nanVol ....
                s_run(r).vol(SwVolMap{10,2}(1):SwVolMap{12,2}(1)-1) ...
                nanVol nanVol nanVol ....
                s_run(r).vol(SwVolMap{12,2}(1):SwVolMap{14,2}(1)-1) ...
                nanVol nanVol nanVol ....
                s_run(r).vol(SwVolMap{14,2}(1):SwVolMap{16,2}(1)-1) ...
                nanVol nanVol nanVol ....
                s_run(r).vol(SwVolMap{16,2}(1):SwVolMap{16,2}(end)+5) ...
                ];
            % split and add 3 blanks
        end        

        clear stimulus temp_nii
        
        % if requested, upsample temporal resolution
        if doUpsample
            % stim ---
            tempstim = p_run(r).stim;
            ups_stim = cell(1,2*length(tempstim));
            ups_stim(1:2:end) = tempstim;
            ups_stim(2:2:end) = tempstim;
            p_run(r).stim = ups_stim;
            clear tempstim ups_stim
            
            % selected timepoints
            tempinc = p_run(r).inc;
            ups_inc = nan(1,2*length(tempinc));
            ups_inc(1:2:end) = tempinc;
            ups_inc(2:2:end) = tempinc;
            p_run(r).inc = ups_inc;
            
            % bold ---
            us_nii=[];
            for v=1:length(p_run(r).vol)
                us_nii=cat(4,us_nii,p_run(r).vol{v});
            end
            fprintf('Upsampling BOLD data...\n');
            us_nii = tseriesinterp(us_nii,TR,TR/2,4);
            for v=1:size(us_nii,4)
                p_run(r).vol{v} = us_nii(:,:,:,v);
            end
            clear us_nii
            
            if Is436
                % stim ---
                tempstim = p_run(r+shiftruns).stim;
                ups_stim = cell(1,2*length(tempstim));
                ups_stim(1:2:end) = tempstim;
                ups_stim(2:2:end) = tempstim;
                p_run(r+shiftruns).stim = ups_stim;
                clear tempstim ups_stim
                
                % selected timepoints
                tempinc = p_run(r+shiftruns).inc;
                ups_inc = nan(1,2*length(tempinc));
                ups_inc(1:2:end) = tempinc;
                ups_inc(2:2:end) = tempinc;
                p_run(r+shiftruns).inc = ups_inc;
                
                % bold ---
                us_nii=[];
                for v=1:length(p_run(r+shiftruns).vol)
                    us_nii=cat(4,us_nii,p_run(r+shiftruns).vol{v});
                end
                fprintf('Upsampling BOLD data...\n');
                us_nii = tseriesinterp(us_nii,TR,TR/2,4);
                for v=1:size(us_nii,4)
                    p_run(r+shiftruns).vol{v} = us_nii(:,:,:,v);
                end
                clear us_nii
            end
            
        end
    end
    fprintf(['Saving ses-' sessions{s} '\n']);
    save(fullfile(out_folder, ['ses-' sessions{s} '-230vols']),'p_run','-v7.3');
    fprintf('Saved result in code folder. Please move it manually...\n')
    clear s_run p_run
end