function [flag]=nic_Detrend(Inpath,N,RemainMean,Outpath)
% Remove the linear trend for the images (Inpath) and output to specified
% directory (Outpath). Whether remain the mean or not is due to the
% option (RetrendMean). The name of image after detrending will be added
% the prefix 'd'.
% 
% FORMAT [flag]=nic_Detrend(Inpath,N,RemainMean,Outpath)
% 
% Input:
% 1. Inpath - the directory of all images to detrend or the directory of 
%             all subjects. Each sub-directory contains one subject's fMRI
%             images.  i.e., the normalised image 
%             e.g.  ['E:\tempData\schizophrenia_filtered']
% 2. N - degrees of Polynomial curve fitting, usually use 1 for removing 
%        linear trend. e.g. [1]
% 2. RemainMean - 1 or 0. 1 means remaining the mean of time
%             course while 0 remove.
%             e.g. [1]
% 3. Outpath - the directory for detrend images with the same directory 
%             structure of Inpath.
%             e.g.   ['E:\tempData\schizophrenia_filtered_detrend']
% 
% Output:
% 1. flag - a flag for all done, 1: successful, 0: fail
% 
% Note:
% 1. All sub directories under Inpath will be involved in this processing.
% So extra files should not exist.
% 2. It is required that only images or only subdirectories containing
% images exist under Inpath.
% 3. It is suggested to detrend linear trend (N=1). If N>1, it takes long
% time to calculate.
% 
%_______________________________________________________________________
% Copyright (C) 2007 Neuroimage Computing Group
%
% Chun Meng, charlieofmeng@gmail.com
% @(#)nic_Detrend.m  ver 1.0, 07/11/29, checked on 07/12/19

t1 = clock;% start time 

if exist(Outpath)~=7
    mkdir(Outpath);
end

tmp=dir([Inpath '\' '*.img']);
if ~isempty(tmp) % only one subject's images under Inpath
    fprintf('detrend: %s\n',Inpath);
    remove(Inpath,N,RemainMean,Outpath);
else
    % subdirectories exist under Inpath
    sub=dir(Inpath);
    num=length(sub);
    for i=3:num   % loop to process every subdirectory
        direct=fullfile(Inpath,sub(i).name);
        if exist(direct)~=7
            error('There is some file other than the image');
        else
            outdirect=fullfile(Outpath,[sub(i).name '_detrend']);
            if exist(outdirect)~=7
                mkdir(outdirect);
                fprintf('Total: %d\tdetrend: %d - %s\n',num-2,i-2,direct);
                remove(direct,N,RemainMean,outdirect);
            else
                fprintf('----%s has existed\n',outdirect);
            end
        end
    end
end

flag=1;

if flag == 1
    t2 = etime(clock, t1);
    fprintf('--------All done! \tElapsed time: %s s\n',num2str(t2));
else
    error('Please carefully check the whole processing step, do again')
end

% ----------------------------sub function-------------------------------
% Read all images and detrend time course under one directory
% ------------------------------------------------------------------------
function remove(Inpath,N,RemainMean,Outpath)
tmp=nic_spm_get('files',Inpath,['*.img']);
Vin=nic_spm_vol(tmp);
TClength = size(tmp,1);
nslices= Vin(1).dim(3);
slicesize=Vin(1).dim(1:2);

Vout=Vin;% only change outpath and filename
for i=1:TClength
    [pth,nm,xt,vr] = fileparts(deblank(Vin(i).fname));
    pth = Outpath;
    Vout(i).fname  = fullfile(pth,['d' nm xt vr]);%改名字,路径，作为输出文件的头信息
end

Vout = nic_spm_create_vol(Vout,'noopen');

% read TC and detrend slice by slice
for i=1:nslices
    A  = nic_spm_matrix([0 0 i]);
    % read TC of a slice
    for j=1:TClength
        slices(:,:,j) = nic_spm_slice_vol(Vin(j),A,slicesize,1);
    end
    
    meanValue=mean(slices,3);
    
    tmp2=reshape(slices,prod(slicesize),TClength);
    if N==1 % remove linear trend        
        tmp3=detrend(tmp2');        
    else
        x=1:TClength;
        tmp3=[];
        for k=1:prod(slicesize)
            p = polyfit(x,tmp2(k,:),N);
            conf=polyval(p,x);
            tmp3(k,:)=tmp2(k,:)-conf;
        end 
        tmp3=tmp3';
    end
    slices=reshape(tmp3',[slicesize TClength]);
    
    if RemainMean
        slices=slices+repmat(meanValue,[1 1 TClength]);
    end
    
    % write out the slice for all volumes
    for m=1:TClength
        Vout(m) = nic_spm_write_plane(Vout(m),slices(:,:,m),i);
    end
end

