% Filtering residuals
% this implements compcorr (Behzadi et al 2007 NI)
% the results should be similar to those after removing GS because GS is
% significantly correlated to CSF
% This procedure will need the tSTD map.
% Ze Wang Dec 2011
% Get parameters
par
org_pwd=pwd;
spm_defaults;
defaults.stats.maxmem=2^22;
defaults.mask.thresh=-inf;
cutoff=0.4;
for sb =1:length(PAR.subjects) % for each subject 
    clear P;
    clear SPM;
    sprintf('Filtering resting data for subject #%g''s data...%g subjects left...',sb,length(PAR.subjects)-sb)
    for c=1:PAR.ncond
    
       
        tsnrdir=PAR.condirs{sb,c};
        P=[];
   
        fprintf('Reading %s\n', PAR.condirs{sb,c});
%         P=spm_select('FPList', PAR.subs(sb).ses(b).asldir{c},['^r.*\w*.img']);
        P=spm_select('ExtFPList', PAR.condirs{sb,c}, ['^ra.*\.nii$'], 1:1000);
        if isempty(P)
            P=spm_select('ExtFPList', PAR.condirs{sb,c}, ['^ra.*\.img$'], 1:1000);
        end
        if isempty(P)
        		P=spm_select('ExtFPList', PAR.condirs{sb,c}, ['^r.*\.nii$'], 1:1000);
        		if isempty(P)
            		P=spm_select('ExtFPList', PAR.condirs{sb,c}, ['^r.*\.img$'], 1:1000);
            	end
        end
        WM  = spm_select('FPList',PAR.condirs{sb,c}, ['^rbk_c2.*\.nii$']);
        CSF = spm_select('FPList',PAR.condirs{sb,c}, ['^rbk_c3.*\.nii$']);
        
        maskfile=spm_select('FPList', PAR.condirs{sb,c}, '^mask_bet.nii');
       maskfile=maskfile(1,:);
       vm=spm_vol(maskfile);
       mask=spm_read_vols(vm);
       mask=mask(:)>0;
        tSTDfile=spm_select('FPList',tsnrdir,'^tSTD_epi.*.nii');
       if size(tSTDfile,1)>1
            tSTDfile=tSTDfile(2,:);
       end
       vts=spm_vol(tSTDfile);
       tsmask=spm_read_vols(vts);
       tsmask=tsmask(mask>0);
       [sts,loc]=sort(tsmask);  % sorting
       nmskvx=length(find(mask>0));
       csfloc=loc( floor(0.98*nmskvx):end);
       
%         vcsf=spm_vol(CSF);
%         csfdat=spm_read_vols(vcsf);
%         csfloc=find(csfdat>0.9);  % used to be 0.85
        
        vwm=spm_vol(WM);
        wmdat=spm_read_vols(vwm);
        wmloc=find(wmdat>0.95);
        
        vd=spm_vol(P);
        startfrom2=0;
        if mod(size(vd,1),2) vd=vd(2:end); startfrom2=1; end
        sts=spm_check_orientations(vd);
        if sts
            dat=spm_read_vols(vd);
        else
            dat=spm_read_vols(vd(1));
            for im=2: size(P,1)
                tmp=spm_read_vols(vd(im));
                dat=cat(4,dat,tmp);
            end
        end
        Tdim=size(dat,4);
        dat=reshape(dat,vd(1).dim(1)*vd(1).dim(2)*vd(1).dim(3),Tdim);
        dat(find(isnan(dat)))=0;
        csfsig=mean(dat(csfloc,:),1);
        csfsig=csfsig-mean(csfsig);
        csfsig=csfsig';
        
        wmdat=mean(dat(wmloc,:),1);
        wmdat=wmdat-mean(wmdat);
        wmdat=wmdat./std(wmdat);
        %pca on the csfsig
        %cleaning up the zigzag pattern first
        ref=-ones(size(P,1),1);  % note we assume label first
        ref(2:2:end)=1;
        wmdat=wmdat';
        
        
        dat = dat(mask,:);
        mdat=mean(dat,2);
        dat=dat-repmat(mdat, 1, Tdim);   
            movefile = spm_select('FPList', PAR.condirs{sb,c}, ['^rp_a.*\.txt$']);
            
            moves = spm_load(movefile);
            moves = moves(1:Tdim,:);

            move_names={'m1' 'm2' 'm3' 'm4' 'm5' 'm6' 'gs' 'wm' 'csf1' 'csf2' 'csf3' 'csf4' 'csf5' 'csf6'};
            refs=moves(:, 1:6);
            
            refs=cat(2, csfsig, wmdat, refs);
                        
            [u,s,v]=svd(refs);
        % take the first 6 components
        nu=u(:,1:8);
%         dat=dat-(nu*(nu'*dat')/(nu'*nu))';
%         dat=dat-(dat'\nu*nu');   % note: this could fail due to the
%         singularity, and extremely large fitting weights can be found
        dat=dat-(dat*nu*nu');
        nmat=mean(dat,2);
        dat=dat-repmat(nmat,1,Tdim);
        [lb,la]=fltbutter(1,0.4);    % if TR=2, omega=1~1/2TR =0.25 (fMax)
        clear filter
        dat=filter(lb,la,dat,[],2); 
        
        [lb,la]=fltbutter(1,0.04,'high');  % high pass  
        clear filter
        dat=filter(lb,la,dat,[],2); 
        
        dat=repmat(mdat, 1, Tdim)+ dat;
        
        img=zeros(vd(1).dim(1),vd(1).dim(2),vd(1).dim(3));
        vo=vd(1);
        vo.fname=fullfile(spm_str_manip(vo.fname,'dh'),['flt_rimages.nii' ]);
        for i=1:Tdim
            
            vo.n=[i 1];
            img(mask)=dat(:,i);
            vo=spm_write_vol(vo,img);
        end
     
    end  % for each condition
end
    cd(org_pwd);