function [onsetTimeVector activeMat activeTP] = nic_TCAcore_mod(data,Z0,...
    percentThreshold,Case)
% The core function of an altered Temporal Clustering Analysis (TCA) method
% which could be used to calculate the spontaneous activition time point
% (onesetTimeVector) in resting state fMRI data(data).
%
% FORMAT [onsetTimeVector activeMat activeTP]=nic_TCAcore(data,Z0,...
%        percentThreshold,Case)
%
% Input:
% 1. data - M-by-N matrix of 4D fMRI data [M:the number of time points;
%                                          N:the number of voxels]
% 2. Z0   - the threshold of activation, with Z>Z0
% 3. percentThreshold
%         - set the percentage of the voxels with Z>Z0 at t0
%           which appeared in the set of voxels with Z>Z0 at higher ranked
%           time points defined as activation, e.g,0.3 (30%)
% 4. Case - selecting the method to get spontaneous activation time course
%           [1,2,or 3]
%        1- Find only one reference time point, which has the greatest
%           total voxels of Z>Z0, and exceeds gaussian assumption, Other
%           time points had to satisfy 2 criteria in order to be considered
%           as Spon-TPs:
%           (a) The number of voxels with Z > 2 at this time point exceeded
%               2.5% of the total voxel number in the ROI.
%           (b) A large proportion (more than 30%) of the voxels with Z > 2
%               at this time point must have appeared in the set of voxels 
%               with Z > 2 at higher ranked time points that have been
%               defined as Spon-TPs. We added this criterion to avoid
%               misinterpreting noise as spontaneous activations.
%           NOTE: if a time point is excluded, all lower-ranked time points
%           are excluded, too.
%       2 - like the case 1, except that if a time point is excluded, the
%           lower-ranked time points still have chance to be selected as
%           spontaneous time points, as long as they comply with the rules
%           a) and b)
%       3 - Find two separate reference time points, the first one has the
%           greatest total voxels of Z>Z0, and exceeds gaussian
%           assumption;the second one was the first time point whose voxels
%           with Z>2 appeared in the set of voxels with Z>2 at the first one
%           is less than 30%.
%           other time points was defined as Spon-TPs of reference 1 or 2
%           follows the rule in Case 2.
%
% Output:
% 1. onsetTimeVector
%         - the vector which shows activation time point (i.e. onset time)
%           when the element = 1, active; when = 0, non-active
% 2. activeMat
%         - matrix contains active voxel at the the active time.
%           NOTE: not the final result, just for showing purpose
% 3. activeTP
%         - list the # of activation voxel in onsetTimeVector
% _________________________________________________________________________
% Copyright (C) 2007 Neuroimage Computering Group, State Key Laboratory of
% Cognitive Neuroscience and Learning
%
% Shuangye Ma, mashuangye@gmail.com
% @(#)nic_TCAcore.m  ver 1.0, 07/11/15
% Modified by Han Zhang, napoleon1982@gmail.com
% @(#)nic_TCAcore.m  ver 1.1, 07/11/19
%
% reference:
% 1. Hunter et al, PNAS, 2005
% 2. Wang et al, Cereb Cortex, 2007


A=data;
S=size(A);
B=[];
T=[];
% select the active voxels of every time point
Lm=mean(A);
Lsd=std(A);
for i=1:S(2)
    B(:,i)=A(:,i)-(Lm(i)+Z0*Lsd(i));
end
A=(B>0);
Atemp=A;

% calculate the number of active voxels of every time point
for i=1:S(1)
    T(i,1)=sum(A(i,:));
end

%Gaussian assumptions: the number of active voxels at an active time point
%should exceed 2.5% of N (i.e.total amount of voxels in consideration)
for i=1:S(1)
    if T(i)<=round(0.025*S(2));
        T(i)=0;
        Atemp(i,:)=0;
    end
end

%time points are ranked in order of the number of active voxels
[Torder,IX]=sort(T,'descend');

if Torder(1)==0
    error('No time point meet the Gaussian assumptions! ');
    return; % the # of active voxels does not exeed gaussian assumption
end

% the time point with greatest total voxels of Z>2 was defined as a spontaneous
% activation,other time points were defined as spontaneous activations if they
%contained voxels with Z>2 that are common to all higher-ranked time points


onsetTimeVector(1,:) =zeros(size(A,1),1);
onsetTimeVector(1,IX(1))=1;
activeTP(1,:)=zeros(size(A,1),1);
activeTP(1,IX(1))=T(IX(1));
SetHigherTP(1,:)=A(IX(1),:);
CF(1,:)=A(IX(2),:);


switch Case
    case 1
        i=2;
        while sum(CF(1,:)&SetHigherTP(1,:))>sum(CF(1,:))*percentThreshold      
            onsetTimeVector(1,IX(i))=1;
            activeTP(1,IX(i))=T(IX(i));
            SetHigherTP(1,:)=CF(1,:)|SetHigherTP(1,:);
            CF(1,:)=A(IX(i+1),:);
            i=i+1;
            if i+1>=S(1)
                break;
            end
        end
    case 2
        for i=2:S(1)-1
            if sum(CF(1,:)&SetHigherTP(1,:))>sum(CF(1,:))*percentThreshold
                onsetTimeVector(1,IX(i))=1;
                activeTP(1,IX(i))=T(IX(i));
                SetHigherTP(1,:)=CF(1,:)|SetHigherTP(1,:);
                CF(1,:)=A(IX(i+1),:);
            end
        end
    case 3
        C=CF(1,:);
        Set=SetHigherTP(1,:);
        flag=0;
        for i=2:S(1)-1
            if sum(C&Set)<=sum(C)*percentThreshold
                if Torder(i)==0
                    break;
                end
                SetHigherTP(2,:)=C;
                CF(2,:)=A(IX(i+1),:);
                activeTP(2,IX(i))=T(IX(i));
                onsetTimeVector(2,IX(i))=1;
                j=i+1;
                flag=1;
                break;
            else
                C=A(IX(i+1),:);
            end
        end
        for i=2:S(1)-1
            if sum(CF(1,:)&SetHigherTP(1,:))>sum(CF(1,:))*percentThreshold
                onsetTimeVector(1,IX(i))=1;
                activeTP(1,IX(i))=T(IX(i));
                SetHigherTP(1,:)=CF(1,:)|SetHigherTP(1,:);
                CF(1,:)=A(IX(i+1),:);
            end
        end
        if flag==1
            onsetTimeVector(2,:)=zeros(size(A,1),1);
            activeTP(2,:)=zeros(size(A,1),1);
            for i=j:S(1)-1
                if sum(CF(2,:)&SetHigherTP(2,:))>sum(CF(2,:))*percentThreshold
                    onsetTimeVector(2,IX(i))=1;
                    activeTP(2,IX(i))=T(IX(i));
                    SetHigherTP(2,:)=CF(2,:)|SetHigherTP(2,:);
                    CF(2,:)=A(IX(i+1),:);
                end
            end
        end
        
end
activeMat=Atemp;
return;









