function clusterstats( tempath, voxdata, real_idx, real_clus, stats)

% get (and unpack) useful variables
nvox            =   size( voxdata, 2);
nsamples        =   stats.nsamples;
sided           =   stats.sided;
cutoff          =   stats.cutoff;
null_methods    =   stats.null_methods;
h0_mean         =   stats.h0_mean;
ttype           =   stats.ttype;

% iterate over methods for generating a null distribution
for nmeth = 1 : length( null_methods)
    
    % preallocate matrices for clusters in null data
    cl_prob     =   NaN(   1, nvox);
    clus_idx    =   zeros( 1, nvox);
    all_clus    =   cell( nsamples, 1);
    
    % iterate over samples
    for ns = 1:nsamples
        
        % create null data for method (current iteration)
        switch null_methods{nmeth}
            case 'perm'
                % randomly multiply each participant by 1 or -1
                this_it = voxdata(:,:,1) .* repmat( (2*round( rand( size( voxdata(:,:,1), 1), 1))-1), 1, nvox);
            case 'labs'
                % retrieve data from shuffled test labels
                this_it = voxdata(:,:,1+ns);
        end
        
        % compute p values based on requested test
        switch ttype
            
            case 'ttest'
                
                [~,pvals]   =   ttest( this_it, h0_mean);
                pvals       =   pvals/(3-sided);
                
            case 'bootstrap'
                
                mus = NaN( nsamples, size( this_it, 2));
                for si = 1 : nsamples
                    mus(si, :)     =   mean( datasample( this_it, size( this_it, 1)));
                end
                
                pvals   =   min( [mean( mus<=h0_mean); 1-mean( mus<=h0_mean)], [], 1)*sided;
                
            case 'rank'
                
                pvals = NaN( 1, nvox);
                for ni = 1:nvox
                    pvals(1, ni)    =   signrank( this_it(:,ni), h0_mean);
                end
                pvals  =   pvals/(3-sided);
        end
        
        % retrieve sign (to separate clusters by sign)
        it_sign         =       (mean( this_it, 1) - h0_mean) > 0;
        ss              =       zeros( size( pvals));
        ss(strfind( [0,it_sign(:)'],[0 1])) = 1;
        ss(strfind( [0,it_sign(:)'],[1 0])) = 1;
        
        % compute maximum cluster size in current iteration
        ii              =       zeros( size( pvals));
        jj              =       pvals < cutoff;             % all sub threshold rois
        ii(strfind( [0,jj(:)'],[0 1])) = 1;                 % onset of sub threshold rois
        ii              =        ii | (ss & jj);            % taking into account sign changes
        idx             =       cumsum(ii).*jj;             % and number the cluster incrementally
        
        % stack cluster sizes of each cluster
        if any( idx)
            all_clus{ns}     =       unique( cellfun( @length, accumarray( idx(jj)', pvals(jj)', [], @(x){x'})));
        else
            all_clus{ns}     =       0;
        end
    end % loop over samples
    
    % find probability per cluster size
    for ncsize = 1:nvox
        cl_prob( ncsize)    =       mean( cellfun( @(x) any(x>=ncsize), all_clus));
    end
    
    % find significant cluster size, given alpha level
    min_clus    =   find( cl_prob < cutoff, 1, 'first');
    
    % list cluster indices of clusters with significant size
    if ~isempty( real_clus)
        sig_clus                =   find( real_clus >= min_clus);
        for ncl = 1:length( sig_clus)
            clus_idx = clus_idx + (real_idx == sig_clus(ncl)).*ncl;
        end
    end
    
    % compute average classifiation data of significant clusters
    k_mean   =   mean( voxdata(:, clus_idx > 0, 1), 2); 
    
    % save data
    save( tempath{nmeth}, 'clus_idx', 'cl_prob', 'min_clus', 'k_mean');
    
end % end of loop over null distribution methods


end

