% Motion_Data: size T x 6 x n_subjects x n_sessions
% TemporalMask_All has size T x n_subjects x n_sessions, with 1s if
% scrubbed out
function [Motion_Matrix,Motion_Matrix_Pos,Motion_Matrix_Neg,Space_Labels,...
    Time_Labels,Session_Labels] = MOT_ANA_ComputeParameters_DS(Motion_Data,...
    TemporalMask_All,n_TB,DS_Factor)

    
    
    n_subjects = size(Motion_Data,3);
    n_sessions = size(Motion_Data,4);
    
    % Differential
    DeltaMot(2:1200,:,:,:) = diff(Motion_Data,1);
    DeltaMot(1,:,:,:) = zeros(6,n_subjects,n_sessions);

    
    
    
    % Number of samples available
    samples_session = size(Motion_Data,1);
    
    % Number of samples to include in each temporal bin
    dur = floor(samples_session/n_TB);
    
    
    
    % Retains only the moments of positive/negative motion change; this is to verify
    % that positive and negative-valued motion moments are equally frequent
    DMPos = DeltaMot;
    DMPos(DMPos < 0) = NaN;

    DMNeg = DeltaMot;
    DMNeg(DMNeg > 0) = NaN;

    

    % Time at which we stand (in number of samples) and associated index
    current_time = 1;

    % To fill with the spatiotemporal information
    Motion_Matrix = [];

    Motion_Matrix_Pos = [];
    Motion_Matrix_Neg = [];

    % As long as we still consider a time lying within the session
    % duration...
    while current_time <= n_TB * dur

        % Separately for positive and negative instantaneous changes, samples
        % the data and sums

        idx_sess = 1;
        
        for ses = 1:n_sessions
            
            idx_fill = 1;
            
            for s = 1:n_subjects
                tmp_pos(idx_fill,1+((idx_sess-1)*6):6+((idx_sess-1)*6)) = nanmean(DMPos(~TemporalMask_All(current_time:current_time+dur-1,s,ses),:,s,ses));
                tmp_neg(idx_fill,1+((idx_sess-1)*6):6+((idx_sess-1)*6)) = nanmean(abs(DMNeg(~TemporalMask_All(current_time:current_time+dur-1,s,ses),:,s,ses)));

                tmp_nonscrubbed(idx_fill,1+((idx_sess-1)*6):6+((idx_sess-1)*6)) = mean(abs(DeltaMot(~TemporalMask_All(current_time:current_time+dur-1,s),:,s,ses)));

                idx_fill = idx_fill + 1;
            end
            
            idx_sess = idx_sess + 1;
        end

        Motion_Matrix = [Motion_Matrix,tmp_nonscrubbed];

        % output data matrices are filled in appropriately: each column
        % highlights one spatiotemporal index

        Motion_Matrix_Pos = [Motion_Matrix_Pos,tmp_pos];
        Motion_Matrix_Neg = [Motion_Matrix_Neg,tmp_neg];
        
        current_time = current_time + dur;
        
        clear tmp_nonscrubbed
        clear tmp_neg
        clear tmp_pos
    end
    
    %%%%% Construction of data labels
    
    n_max = size(Motion_Matrix,2);
    

    
    % Space
    Space_Labels(1,:) = [1:6:n_max-5];
    
    for t = 2:6
        Space_Labels(t,:) = 1 + Space_Labels(t-1,:);
    end
    
    % Time
    Time_Labels(1,:) = [1:6*n_sessions];
    
    for t = 2:n_TB
        Time_Labels(t,:) = 6*n_sessions + Time_Labels(t-1,:);
    end
    
    % Sessions
    Slstart = [1:n_sessions*6:n_max];
    
    Session_Labels = Slstart;
    for a = 1:n_TB-1
        Session_Labels = [Session_Labels,Slstart+a];
    end
    
    for ses = 2:n_sessions
        Session_Labels(ses,:) = 6*(ses-1) + Session_Labels(1,:);
    end
    
    
end