function out=legTrackerErrorCheck(data)
%legTrackerErrorCheck takes the output of legTrackerCleanData and checks
%the leg-position vectors for tracking errors. These are defined as
%instances in which the frame-to-frame position changes by more than 5 x
%the standard deviation of the frame-to-frame motion across the entire
%vector. Subsequent frames are considered to be errors until the
%leg position is within 5-sigma of the previous non-error position, or
%within one sigma of the overall median position.
%Usage:
% out=legTrackerErrorCheck(data)
% data = cleaned data, post- median filtering and interpolation, output
%        from legTrackerCleanData


numFrames=size(data,1);


% perform check on leg-position vectors, not ball motion vectors
for i=2:13
    
    % identify "baseline value"
    medI=median(data(:,i));
    
    % error cutoff in frame-to-frame motion: 5 sigma
    deltaCutoff=5*std(data(1:end-1,i)-data(2:end,i));
    
    % "data definitely back in normal range" cutoff: 1 sigma
    medCutoff=1*std(data(1:end-1,i)-medI);
    
    % are we looking at an error frame? indicator boolean 
    onNaN=0;
    
    % inspect all frames
    for j=2:numFrames
        
        %if we are not looking at error frames
        if onNaN==0;
            
            % if jump is too big, surpassing motion cutoff
            if abs(data(j,i)-data(j-1,i))>deltaCutoff
                
                % start sequence of error frames
                NaNStart=j;
                onNaN=1;
            end
        else   %if we are looking at error frames
            
            % if we within motion cutoff of last non-error frame or within
            % normal range cutoff of the global median value
            if abs(data(j,i)-data(NaNStart-1,i))<deltaCutoff || abs(data(j,i)-medI)<medCutoff
                
                % no longer in an sequence of error frames
                NaNEnd=j-1;
                onNaN=0;
                a1=data(NaNStart-1,i);
                a2=data(NaNEnd+1,i);
                
                % replace sequence of error frames with linear
                % interpolation of flanking frame values
                if a1==a2
                    data(NaNStart-1:NaNEnd+1,i)=ones(size(data(NaNStart-1:NaNEnd+1,i)))*a1;
                else
                data(NaNStart-1:NaNEnd+1,i)=a1:(a2-a1)/(NaNEnd-NaNStart+2):a2;
                end
            end
        end
        
    end
    
end

out=data;