function [flag]=nic_Seed2corr(Inpath,Filefilter,seedTC,covariant,poly,Outpath,Prefix,mask)
% Calculate  Pearson correlation coefficients between seed TC and all 
% voxels' TC and output to a correlation image for one subject
% 
% FORMAT [flag]=nic_Seed2corr(Inpath,seedTC,covariant,Outpath)
% 
% Input:
% 1. Inpath - the directory of all subjects. Each sub-directory contains 
%             one subject's fMRI images.  
%             e.g.  ['E:\tempData\schizophrenia_filtered']
% 2. Filefilter - file filter set. e.g., ['w'] or ['sw'] meaning that only
%               process the image involving the 'w' or 'sw' as the prefix
% 3. seedTC - a row or column vector, which is the TC of specific voxel or 
%             mean TC of ROI and used to correlated with whole brain'TCs
% 4. covariant - a matrix, whose column vectors are the confounding TCs
%                e.g. [head motion, TC of WM and CSF,global trend]
% 5. poly -  degrees of Polynomial curve fitting, [] means to remove
%            nothing, [0] means to remove mean value, [1] means to remove
%            linear trend, and so on
% 6. Outpath - the directory for images of correlation 
%              e.g. ['E:\tempData\schizophrenia_filtered_detrend']
% 7. Prefix - the outfile's name can be ['Prefix'_directory name]
% 8. mask - a binary mask [path\*.img] is required so that the value of voxels 
%           outside the mask will be set to 0
% 
% 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.
% 4. refering to REST coded by Song Xiaowei,resting-fmri.sourceforge.net
%_______________________________________________________________________
% Copyright (C) 2007 Neuroimage Computing Group
%
% Chun Meng, charlieofmeng@gmail.com
% @(#)nic_Seed2corr.m  ver 1.0, 07/12/27

t1 = clock;% start time 
if exist(Outpath)~=7
    mkdir(Outpath);
end
currentpath=pwd;
% Add path spm5 and use it
spmPath=nic_getSPMPath('nifti');
addpath(spmPath);

% Make seedTC a column vector
if size(seedTC,1)==1
    seedTC=seedTC';
elseif size(seedTC,2)~=1
    error('seedTC is not a column or a row');
end

cd(Inpath);
Datafolder=dir;% refer to nic_BandPassFilter.m
Datafolder=Datafolder(3:end);
% Determine whether the directory is the deepest catalog
if isdir(Datafolder(1).name)
    nsubjects=size(Datafolder,1);
    Seconddir=1;    
else
    Seconddir=0;
    nsubjects=1;
end
for iFolder=1:nsubjects
    if Seconddir ==1
        Path=[Inpath '\' Datafolder(iFolder).name];
        subjectname=Datafolder(iFolder).name;
    else
        Path=Inpath;
        tmp=find(Path=='\');
        % take the directory name of subject as the name of outfile 
        tmpname=Path(tmp(end)+1:end);
        if isempty(tmpname)
            tmpname='wholebraincorrelation';
        end
        subjectname=tmpname;
    end

    % get file hdr
    FileString=nic_spm_select('List',Path, ['^' Filefilter '.*\.img$']);%
    Vin=nic_spm_vol(FileString);
    
    TClength= size(FileString,1);% number of time points
    %     if size(seedTC,1)~=TClength
    %         error('lengths of seedTC and sliceTC are different!');
    %     end
    nslices= Vin(1).dim(3);
    slicesize=Vin(1).dim(1:2);
    
    % create new header files for output image, one correlation image for 
    % one subject
    Vout = Vin(1);
    [pth,nm,xt,vr] = fileparts(Vout.fname);
    Vout.fname  = fullfile(Outpath,[Prefix '_' subjectname xt vr]);
    if isfield(Vout,'descrip'),
        desc = [Vout.descrip ' '];
    else
        desc = '';
    end
    Vout = nic_spm_create_vol(Vout);    

    % read TCs slice by slice
    slices = zeros([slicesize TClength]);% 3D matrix, 3rd dim is time
    data_r = zeros([slicesize nslices]);
    for k = 1:nslices
        B  = nic_spm_matrix([0 0 k]);
        for m=1:TClength
            slices(:,:,m) = nic_spm_slice_vol(Vin(m),B,slicesize,1);
        end

        % Revise the dataset with covariables at first
        % 去协变量，不用regress，mean等函数，防止溢出
        if ~isempty(poly)
            xx=[1:TClength];
            xx=xx';
            for iDegree=0:poly
                poly_covariant(:,iDegree+1)=xx.^poly;
            end
            covariant=[covariant poly_covariant];
        end
                
        if ~isempty(covariant)
            seedTC =seedTC_RegressOutCovariables(seedTC, covariant);
            slices =sliceTC_RegressOutCovariables(slices, covariant);
        end

        % Compute the Pearson coeffecient, return a 3D brain whose elements
        % are the coefficients
        data_r(:,:,k) =PearsonCorrTC(seedTC, slices);
        fprintf('.');% 动态显示进度
%         Vout = nic_spm_write_plane(Vout,data_r(:,:,k),k);
    end
%     Vout = nic_spm_close_vol(Vout);
   
    if ~isempty(mask)
        % apply the mask
        [Ymask,Head] = nic_ReadNiftiImage(mask);
        Ymask=logical(Ymask);
        if any(size(Ymask)-size(data_r))
            error('mask size must be the same with input images');
        end
        data_r=data_r.*Ymask;
    end
    
    % write out
    Vout = nic_spm_write_vol(Vout,data_r);
    fprintf('Done:%s\n',subjectname);
end

cd(currentpath);
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-------------------------------
% regress out Covariables and calculate correlation coefficients
% ------------------------------------------------------------------------

% To calculate correlation coefficients with 1D and 3D data
% Refer to fc.m coded by Song Xiaowei
function CorrCoefslice=PearsonCorrTC(seedTC, sliceTC);
    if length(seedTC)~=size(sliceTC,3)
        error('lengths of seedTC and sliceTC are different!');
    end
    [nDim1, nDim2, nDim3]=size(sliceTC);% M*N*TClength
    % r=E[(X-EX)(Y-EY)]/[std(X)*std(Y)]

    %Remove the mean
    sliceTC = sliceTC - repmat(sum(sliceTC,3)/nDim3,[1, 1, nDim3]);
    seedTC =seedTC -repmat(sum(seedTC)/nDim3, [size(seedTC,1), size(seedTC,2)]);

    % For denominator  相关系数的分母[std(X)*std(Y)]
    sliceSTD = squeeze(std(sliceTC, 0, 3));% Set flag to 0 to normalize Y by n-1
    tmp1 = sliceSTD * std(seedTC);

    % (1*sampleLength) A matrix * B matrix (sampleLength * VoxelCount)
    % 相关系数的分子E[(X-EX)(Y-EY)]
    tc2D =reshape(sliceTC, nDim1*nDim2, nDim3)';% 将voxel拉成一行，列为时间序列
    seedTC =reshape(seedTC, 1, nDim3);% 不管原来是行还是列，转为一行
    tmp2 = seedTC * tc2D /(nDim3 -1);% 结果为一行
    tmp3 =squeeze(reshape(tmp2', nDim1, nDim2, 1));% 转回一片

    tmp1(find(tmp1==0))=inf;% 0/0=NaN,suppress NaN to zero when denominator is zero
    CorrCoefslice =tmp3 ./tmp1;

function Result =sliceTC_RegressOutCovariables(sliceTC, ABasisFunc)
    % Result =( E - X(X'X)~X')Y
    [nDim1, nDim2, nDim3]=size(sliceTC);

    %Make sure the 1st dim of ABasisFunc is nDim3 long
    if size(ABasisFunc,1)~=nDim3, error('The length of Covariables don''t match with the volume.'); end

    % (1*sampleLength) A matrix * B matrix (sampleLength * VoxelCount)
    tmp =reshape(sliceTC, nDim1*nDim2, nDim3)';
    Result =(eye(nDim3) - ABasisFunc * inv(ABasisFunc' * ABasisFunc)* ABasisFunc') * tmp;
    Result =reshape(Result', nDim1, nDim2, nDim3);

function Result =seedTC_RegressOutCovariables(seedTC, ABasisFunc)
    % Result =( E - X(X'X)~X')Y
    % seedTC would better be a column vector   
    seedTC =reshape(seedTC, prod(size(seedTC)), 1);% 不管原来是行还是列，转为一列
    %Make sure the 1st dim of ABasisFunc is nDim3 long
    if size(ABasisFunc,1)~=length(seedTC), error('The length of Covariables don''t match with the seedTC.'); end
    % (1*sampleLength) A matrix * B matrix (sampleLength * VoxelCount)
    Result =(eye(size(seedTC, 1)) - ABasisFunc * inv(ABasisFunc' * ABasisFunc)* ABasisFunc') * seedTC;
