% Apply realign parameters (rp*.txt generated by SPM) to grey matter, white matter and CSF images. Buy doing so, one could evaluate the partial voluming effects on BOLD signal. These script is written for a porject of Yan, C.G., Cheung, B., Kelly, C., Colcombe, S., Craddock, R.C., Di Martino, A., Li, Q., Zuo, X.N., Castellanos, F.X., Milham, M.P., 2013. A comprehensive assessment of regional variation in the impact of head micromovements on functional connectomics. Neuroimage 76, 183-201. However, we didn't include this part into the final paper.
%__________________________________________________________________________
% Written by YAN Chao-Gan 120229.
% The Nathan Kline Institute for Psychiatric Research, 140 Old Orangeburg Road, Orangeburg, NY 10962, USA
% The Phyllis Green and Randolph Cowen Institute for Pediatric Neuroscience, New York University Child Study Center, New York, NY 10016, USA
% ycg.yan@gmail.com
% http://rfmri.org/Scripts_YAN

%%%%%%%%%%%%%%
%Reslice the summed GM into Head Motion functional space
hld = 1;

AutoDataProcessParameter.DataProcessDir='/Processing';
OutputDir='/Processing';

load('/Processing/SubID.mat');
AutoDataProcessParameter.SubjectID=SubID;
AutoDataProcessParameter.SubjectNum=length(AutoDataProcessParameter.SubjectID);

parfor i=1:AutoDataProcessParameter.SubjectNum
    
    if 7==exist([AutoDataProcessParameter.DataProcessDir,filesep,'RealignParameter',filesep,AutoDataProcessParameter.SubjectID{i}],'dir')
        DirMean=dir([AutoDataProcessParameter.DataProcessDir,filesep,'RealignParameter',filesep,AutoDataProcessParameter.SubjectID{i},filesep,'mean*.img']);
        if isempty(DirMean)  %YAN Chao-Gan, 111114. Also support .nii files.
            DirMean=dir([AutoDataProcessParameter.DataProcessDir,filesep,'RealignParameter',filesep,AutoDataProcessParameter.SubjectID{i},filesep,'mean*.nii']);
        end
        if ~isempty(DirMean)
            RefFile=[AutoDataProcessParameter.DataProcessDir,filesep,'RealignParameter',filesep,AutoDataProcessParameter.SubjectID{i},filesep,DirMean(1).name];
        end
        
    end
    
    DirRP=dir([AutoDataProcessParameter.DataProcessDir,filesep,'RealignParameter',filesep,AutoDataProcessParameter.SubjectID{i},filesep,'rp*']);
    RPFile=[AutoDataProcessParameter.DataProcessDir,filesep,'RealignParameter',filesep,AutoDataProcessParameter.SubjectID{i},filesep,DirRP.name];

    SourceFileList=[];
    SourceFileList=[SourceFileList;{[OutputDir,filesep,'GMWMCSF_BlockSummed',filesep,AutoDataProcessParameter.SubjectID{i},filesep,'GM_BlockSummed.nii']}];
    SourceFileList=[SourceFileList;{[OutputDir,filesep,'GMWMCSF_BlockSummed',filesep,AutoDataProcessParameter.SubjectID{i},filesep,'WM_BlockSummed.nii']}];
    SourceFileList=[SourceFileList;{[OutputDir,filesep,'GMWMCSF_BlockSummed',filesep,AutoDataProcessParameter.SubjectID{i},filesep,'CSF_BlockSummed.nii']}];
    
    OutputFileList=[];
    OutputFileList=[OutputFileList;{[OutputDir,filesep,'GM_HM',filesep,AutoDataProcessParameter.SubjectID{i}]}];
    OutputFileList=[OutputFileList;{[OutputDir,filesep,'WM_HM',filesep,AutoDataProcessParameter.SubjectID{i}]}];
    OutputFileList=[OutputFileList;{[OutputDir,filesep,'CSF_HM',filesep,AutoDataProcessParameter.SubjectID{i}]}];
    

    RP=load(RPFile);
    
    sinq1=sin(RP(:,4));
    sinq2=sin(RP(:,5));
    sinq3=sin(RP(:,6));
    cosq1=cos(RP(:,4));
    cosq2=cos(RP(:,5));
    cosq3=cos(RP(:,6));
    
    [RefData RefHead]   = y_Read(RefFile);
    TargetMat=RefHead.mat;
    TargetDim=RefHead.dim;
    
    % In case there exist reorient matrix (interactive reorient after head motion correction and before T1-Fun coregistration) 
    ReorientMat=eye(4);
    if exist([AutoDataProcessParameter.DataProcessDir,filesep,'ReorientMats'])==7
        if exist([AutoDataProcessParameter.DataProcessDir,filesep,'ReorientMats',filesep,AutoDataProcessParameter.SubjectID{i},'_ReorientFunImgMat.mat'])==2
            ReorientMat_Interactively = load([AutoDataProcessParameter.DataProcessDir,filesep,'ReorientMats',filesep,AutoDataProcessParameter.SubjectID{i},'_ReorientFunImgMat.mat']);
            ReorientMat=ReorientMat_Interactively.mat*ReorientMat;
        end
    end
        
    for iFile=1:length(OutputFileList)
        mkdir(OutputFileList{iFile});
        cd(OutputFileList{iFile});

        [SourceData SourceHead]=y_Read(SourceFileList{iFile});
        
        [x1,x2,x3] = ndgrid(1:TargetDim(1),1:TargetDim(2),1:TargetDim(3));
        d     = [hld*[1 1 1]' [1 1 0]'];
        C = spm_bsplinc(SourceHead, d);
        v = zeros(TargetDim);
        
        for t=1:size(RP,1)

            M1=[1       0        0     0;...
                0    cosq1(t)  sinq1(t)  0;...
                0    -sinq1(t) cosq1(t)  0;...
                0       0        0     1;];
            
            M2=[cosq2(t)  0    sinq2(t)     0;...
                0        1       0        0;...
                -sinq2(t) 0    cosq2(t)     0;...
                0       0        0        1;];
            
            M3=[cosq3(t)   sinq3(t)   0     0;...
                -sinq3(t)  cosq3(t)   0     0;...
                0           0       1     0;...
                0           0       0     1;];
            
            MT=[1    0     0     RP(t,1);...
                0    1     0     RP(t,2);...
                0    0     1     RP(t,3);...
                0    0     0     1;];
            
            M_RigidBodyTransform=MT*M1*M2*M3;

            M = inv(ReorientMat*inv(M_RigidBodyTransform)*inv(ReorientMat)*SourceHead.mat)*TargetMat; % M = inv(TargetMat\SourceHead.mat) in spm_reslice.m
            y1   = M(1,1)*x1+M(1,2)*x2+(M(1,3)*x3+M(1,4));
            y2   = M(2,1)*x1+M(2,2)*x2+(M(2,3)*x3+M(2,4));
            y3   = M(3,1)*x1+M(3,2)*x2+(M(3,3)*x3+M(3,4));
            
            OutVolume    = spm_bsplins(C, y1,y2,y3, d);
            
            OutHead = SourceHead;
            OutHead.mat      = TargetMat;
            OutHead.dim(1:3) = TargetDim;
            OutHead.pinfo = [1;0;0];
            OutHead.dt    =[16,0];
            y_Write(OutVolume,OutHead,[OutputFileList{iFile},filesep,sprintf('%.8d', t),'.nii']);
            
        end
    end
    
end

