function ROI_pipeline_subfunc( ptcs, roi_names, voxstep, voxmax, wdir, delall)


%% ROIsettings

% write settings
if ~exist( 'wdir', 'var')
    wdir = '1';
end
clear_old   =   0;        % remove all existing rois in directory (i.e., clean up)
if exist( 'delall', 'var')
   clear_old = delall; 
end

% Mask names & directories, contrast directories & numbers, and threshold
mask_names  =   {  'LOC_LR',       'PPA_LR',       'EVC_(BA1718)_mask',  'EVL',            'EVR',            'V1_(BA17)_mask', 'V1L',            'V1R',          'V2_(BA18)_mask', 'V2L',          'V2R',            'TOS_LR',       'RSC_LR',       'FFA_LR',       'LOL',          'LOR',          'PPL',           'PPR',          'TOL',          'TOR',          'RSL',          'RSR',          'FFL',          'FFR'          }; 
mask_dirs   =   {  'Func',         'Func',         'AAL',                'AAL',            'AAL',            'AAL',            'AAL',            'AAL',          'AAL',            'AAL',          'AAL',            'Func',         'Func',         'Func',         'Func',         'Func',         'Func',          'Func',         'Func',         'Func',         'Func',         'Func',         'Func',         'Func'         };
con_dirs    =   {  'localizer_s3', 'localizer_s3', 'training_s3',        'training_s3',    'training_s3',    'training_s3',    'training_s3',    'training_s3',  'training_s3',    'training_s3',  'training_s3',    'localizer_s3', 'localizer_s3', 'localizer_s3', 'localizer_s3', 'localizer_s3', 'localizer_s3',  'localizer_s3', 'localizer_s3', 'localizer_s3', 'localizer_s3', 'localizer_s3', 'localizer_s3', 'localizer_s3' };         
con_nums    =   [  1,              3,              1,                    1,                1,                1,                1,                1,              1,                1,              1,                3,              3,              2,              1,              1,              3,               3,              3,              3,              3               3,              2,              2              ];                
vox_alpha   =   [  0.05,           0.05,           0.05,                 0.05,             0.05,             0.05,             0.05,             0.05,           0.05,             0.05,           0.05,             0.05,           0.05,           0.05,           0.05,           0.05,           0.05,            0.05,           0.05,           0.05,           0.05,           0.05,           0.05,           0.05           ];

% retrieve requested ROIs
if isempty( roi_names)
   roi_names = mask_names; 
end
wregs       =   find( contains( mask_names, roi_names));
nrois       =   length( roi_names);

% set defaults
if nrois > 1 && length( voxstep) == 1
    voxstep = repmat( voxstep, 1, nrois);
end
if nrois > 1 && length( voxmax) == 1
    voxmax = repmat( voxmax, 1, nrois);
end


%% file info 

% directories
mask_root   =   fileparts( mfilename('fullpath'));
proj_root   =   fullfile( mask_root, '..','..');

% set filename to store ROI info
rep_dir     =   fullfile( mask_root, 'reports');
roi_info_f  =   fullfile( rep_dir, 'roi_info.mat');

% coregistration job
cjob_templ  =   'coreg_mni.mat';


%% Create ROIs

% run through participants
for nptc = 1 : length( ptcs)

    % retreive participant and report progress
    ptc = sprintf( '%02g', ptcs(nptc));
    fprintf( '\n\nParticipant %g/%g (%s)\n', nptc, length(ptcs), ptc);
    
    % load participant info
    load( fullfile( proj_root, 'SUBJECT level', ptc, 'exp_info.mat'), 'exp_info');
        
    % set participant contrast directory
    ptc_con_dir     =   fullfile( proj_root, 'SUBJECT level', ptc, 'univariate');
    
    % iterate through ROIs
    for nroi = 1:nrois
        
        % retrieve ROI from roi list, and index of voxel counts
        wroi    =   wregs(nroi);
        wvox    =   find( contains( roi_names, mask_names{wroi}(1:3)));
        
        % report progress
        fprintf( ' - ROI %g/%g, started at %s \n', nroi, nrois, datestr(now, 'HH:MM:SS'));
        
        % set participant destination directory
        if contains( con_dirs{wroi}, 'nonorm')
            mask_type = 'native';
        else
            mask_type = 'mni';
        end
        ptc_des_dir     =   fullfile( proj_root, 'SUBJECT level', ptc, 'roi', sprintf( '%s%s', mask_type, wdir));
        
        % delete output directory if requested
        if clear_old && nroi==1
            rmdir( ptc_des_dir, 's');
        end
        
        % create output directory if it does not exist
        if ~exist( ptc_des_dir, 'dir')
            mkdir( ptc_des_dir)
        end
        
        % retrieve suffix for smoothing, and run type for run count
        mysm = strfind( con_dirs{wroi}, '_s');
        if isempty( mysm)
            suffix = '';
            run_type = con_dirs{wroi};
        else
            suffix = con_dirs{wroi}(mysm:mysm+2);
            run_type = con_dirs{wroi}(1:mysm-1);
        end
        
        % determine number of runs, to calculate minimal t-statistic
        nruns   =   height( exp_info( contains( exp_info{ :, 'run_type'}, run_type),:));
        min_t   =   icdf('t', 1-vox_alpha(wroi), nruns-1);   % t cut-off
        
        % retrieve path and file of mask
        p = fullfile( mask_root, mask_dirs{wroi});
        f = mask_names{wroi};
            
        % choose pipeline depending on eventual space
        if strcmp( mask_type, 'native')
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % RESLICE MASKS TO NATIVE SPACE %
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            % define reference for reslicing
            reference       =   fullfile( ptc_con_dir, con_dirs{1}, 'mask.nii');

            % locate inversion parameters
            filelist        =   dir( fullfile( ptc_con_dir, '..', 'anatomical', 'anatomical', 'iy_*.nii'));
            sn_matname      =   fullfile( filelist.folder, filelist.name);
            
            % inverse normalize
            myspm_invnormalise_spm12( fullfile( p, sprintf( '%s.nii', f)), sn_matname, reference)
            
            % smooth resliced mask
            hdr_r       =   spm_vol( fullfile( p, sprintf( 'rw%s.nii', f)));
            newname     =   fullfile( p, sprintf( 's2rw%s.nii', f));
            
            spm_smooth( hdr_r, newname, 2);
            hdr_rs       =   spm_vol( fullfile( p, sprintf( 's2rw%s.nii', f)));
            
            vol = spm_read_vols( hdr_rs);
            vol = vol > 0.05;
            spm_write_vol( hdr_rs, vol);
            
            % organize resliced files
            resliced_mask   =   fullfile( ptc_des_dir, sprintf( 'rw_%s.nii', f));
            movefile( fullfile( p, sprintf( 's2rw%s.nii', f)), resliced_mask);
            delete(   fullfile( p, sprintf( 'rw%s.nii', f)));
            delete(   fullfile( p, sprintf( 'w%s.nii',  f)));

        else
            
            % check if template is already resliced to correct dimensions
            if ~exist( fullfile( p, sprintf( 'r%s.nii', f)), 'file')
                
                % load coregistration job (WARNING: method is deprecated)
                warning( 'Coregistration method is deprecated. Use "ROI_from_AAL" to reslice instead');
                load( cjob_templ);
                
                % set source (anatomical mask for intersecting roi)
                matlabbatch{1}.spm.spatial.coreg.estwrite.source{1}         =   fullfile( p, sprintf( '%s.nii', f));
                
                % set reference (subject/mni whole brain mask)
                my_ref                                                      =   con_dirs{find(~contains(con_dirs,'nonorm'),1)};
                matlabbatch{1}.spm.spatial.coreg.estwrite.ref{1}            =   fullfile( ptc_con_dir, my_ref, 'mask.nii');
                
                % set prefix r
                matlabbatch{1}.spm.spatial.coreg.estwrite.roptions.prefix   =   'r';
                
                % Save & run file:
                save( cjob_templ, 'matlabbatch');
                spm_jobman( 'run', cjob_templ);
            end
            
            % Retrieve resliced mask
            resliced_mask   =   fullfile( p, sprintf( 'r%s.nii', f));
        end
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % INTERSECT (RESLICED) MASKS WITH top voxels from T-MAPS %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        % load resliced mask 
        hdr_mask    =   spm_vol( resliced_mask);
        vol_mask    =   logical( spm_read_vols( hdr_mask));
        
        % load contrast file (t-map)
        hdr_cont    =   spm_vol( fullfile( ptc_con_dir, con_dirs{wroi}, sprintf( 'spmT_%04g.nii', con_nums(wroi))));
        vol_cont    =   spm_read_vols( hdr_cont);
    
        % get all voxels in con-image that are also in mask and t > t_min
        vox_ind             =   find( vol_mask & logical( vol_cont) & vol_cont > min_t);
        vox_all             =   vol_cont(vox_ind);
        
        % sort all voxels according to t-value
        [~,vox_ind_sorted]  =   sort( vox_all,'descend');
        
        % all voxels meeting criteria:
        vox_count   =   numel(vox_ind);
        if ~ isempty( voxstep)
            nVoxels     =   [ voxstep(wvox):voxstep(wvox):voxmax(wvox), vox_count];
        else
            nVoxels     =   vox_count;
        end

        % intersect
        for nvox = 1 : length( nVoxels)
        
            % number of included voxels (maximum or vox count)                
            numvox = min( [nVoxels(nvox), vox_count]);
              
            % name and save ROI
            if nvox == length( nVoxels)
                voxname = 'all';
            else
                voxname = sprintf( '%g', nVoxels(nvox));
            end
            name_roi        = sprintf( '%s_%s_%s%s.nii', mask_names{wroi}(1:3), con_dirs{wroi}(1:4), voxname, suffix);
            hdr_mask.fname  = fullfile( ptc_des_dir, name_roi);
                
                
            % if there are voxels in mask, intersect with t-map
            if vox_count
                
                % select top voxels
                vox_ind_nVoxels                     =   vox_ind_sorted(1:numvox);
                vol_roi                             =   false( size( vol_mask));
                vol_roi(vox_ind(vox_ind_nVoxels))   =   true;
                
                % write volume
                spm_write_vol( hdr_mask, vol_roi);
                                
            else
                % issue warning if insufficient voxels are included
                warning( ' Participant %s has only %g voxels in roi [%s], which is less than %g.\n No ROI was created.', ...
                    ptc, vox_count, name_roi, min_vox);
                hdr_mask.fname = 'NO_ROI_CREATED!';
            end
            
            % store roi info
            roi_info.(con_dirs{wroi})(nvox).voxels(str2double(ptc),1)   =   numvox;
            roi_info.(con_dirs{wroi})(nvox).mask                        =   mask_names{wroi};
            roi_info.(con_dirs{wroi})(nvox).smoothing                   =   suffix(2:end);
            roi_info.(con_dirs{wroi})(nvox).space                       =   mask_type;
            roi_info.(con_dirs{wroi})(nvox).location{str2double(ptc),1} =   hdr_mask.fname;
            roi_info.(con_dirs{wroi})(nvox).voxel_alpha                 =   vox_alpha(wroi);
            
        end % end of voxel loop
    end % end of ROI loop
end % end of subject loop

% save roi info, and report progress
if ~exist( rep_dir, 'dir')
    mkdir( rep_dir);
end
save( roi_info_f, 'roi_info');
fprintf( '\nDone\n\n');
