function OUTPUT = A_MotionDetect(MSource,Th,Int,MinDur,fs,ThM)
% OUTPUT = A_MotionDetect(MSource,Th,Int,MinDur,fs,ThM) 
% This program detects motion vs resting with trehshold, and export 
%      timepoints (by frame number) that indicate the Moving or silencing timpoint internals. 
% OUPUT={DetPoints, DetPairs};
%  Output.MoveLable:  marked motion as 1 in the corresponding slots of 
%                     "ESource" array. (Same demension and size as ESource)
%  Output.MoveInt: provide an n x 2 array for time slots of detected moving
%                   Interval. [On slot, Off slot; ...] of ESource 
%  Output.RestInt: provide an n x 2 array for time slots of detected
%                  rest interval
% Detailed descrption
%"ESource": a 1D array for threshold detection
% "Th": Threshold
% "Int": in ms, the largest time interval between two detected signal to be
%        view as a continuous signal
%       (if Int is 1 sec, two above Th points are 0.7 sec, they are treated
%       as one continuous signal)
% "MinDur": in sec; minimal duration for detecting continuous motion (or
%           silent) periods
% "F": sampling frequency for the Source
% "ThM": Thresholding method: 1==positive thresholding; -1==negative
%        tresholding
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%Main function%%
%% unify the array orientation %%
if size(MSource,1)>=2;
    MSource=MSource';
end
if Int*fs/1000==1
    TpInt=1;
else
    TpInt(1:Int*fs/1000-1)=1; % defined the Timepoints(Tp)for detecting interval
end
%% thresholding source
if ThM==1;
    Data=(MSource>=Th);
    %EPSource(2,:)=(EPSource>=Th);
else if ThM==-1;
        Data=(MSource<=Th);
        %EPSource(2,:)=(EPSource<=Th);        
    else
        error('Wrong threhold method')
        return
    end
end
    
%% define detection points
Data=double(Data);
DTh=conv(Data,TpInt);
Data(2,:)=DTh(1:length(Data)); 
Data=logical(Data);
%% detect duration 
Changepoints=Data(2,:)-[0,Data(2,1:end-1)];
DetPairs(:,1)=find(Changepoints==1); % find "on" points
DetPairs(:,2)=length(Changepoints)+1;
Offpoints=find(Changepoints==-1);
DetPairs(1:length(Offpoints),2)=Offpoints; % on-off pairs for moving periods
DetPairs(:,3)=DetPairs(:,2)-length(TpInt); % trim the last 0 due to convul of each interval ends
DetPairs(:,4)=DetPairs(:,3)-DetPairs(:,1)+1;

DetPair2(:,1)=[1;DetPairs(:,2)]; % off-on pairs for rest periods
DetPair2(:,2)=[DetPairs(:,1);length(Changepoints)+1];
if DetPair2(1,2)==1; DetPair2(1,:)=[];end;
if DetPair2(end,1)==DetPair2(end,2); DetPair2(end,:)=[];end;
DetPair2(:,3)=DetPair2(:,2)-length(TpInt); % trim to minimize transition effect
DetPair2(:,4)=DetPair2(:,3)-DetPair2(:,1)+1;

MoveInt=DetPairs(find(DetPairs(:,4)>=MinDur*fs),[1,3]);
RestInt=DetPair2(find(DetPair2(:,4)>=MinDur*fs),[1,3]);

%% add lables to datapoints based on MoveInt and RestInt;
LableTimePoints(1:length(Changepoints))=0;
for i = 1:1:size(MoveInt,1);
    LableTimePoints(MoveInt(i,1):MoveInt(i,2))=1;
end
for j=1:1:size(RestInt,1);
    LableTimePoints(RestInt(j,1):RestInt(j,2))=-1;
end
%% Define output %%%%%
OUTPUT.RawMoveLable=Data(2,:);
OUTPUT.MoveInt=MoveInt;
OUTPUT.RestInt=RestInt;
OUTPUT.AssignedMoveLable=LableTimePoints;
end