function cor_vol = cba_cor( z, rbrain, keep,Xin,isbrain)
%z = nxnxnumslices, rbrain=nxnxnumslicesxtc
%keep: the time points used in the analysis
%output:
%rho_vol = nxnxnumslices
numslices = size(z,3);
npix = size(z,1)*size(z,2);

if size(Xin,3)>1
    numresponses=size(Xin,2); %the number of stimuli in the experiment
else
    numresponses=1; %the number of stimuli in the experiment
end
arbrain =cba_sig2avg(rbrain(:,:,:,keep), z,isbrain);%each voxel receives the cluster average value
clear rbrain;

%reshape elements
Y=reshape(arbrain,npix, numslices, size(arbrain,4));
isbrain = reshape(isbrain, npix,numslices);

%output 
cor_vol = zeros(size(isbrain));
for slice=1:numslices
    slice
    if numresponses>1
        X=[squeeze(Xin(keep,:,slice))];
    else
        X=[squeeze(Xin(keep,slice))];
    end
    Df=length(keep)-rank(X);
    ipix = 1:npix;
    pix = ipix(isbrain(:,slice)==1);
    if ~isempty(pix)
        Yslice =(squeeze(Y(pix,slice,:)))';
        pinvX=pinv(X);
        betahat=pinvX*Yslice;
        resid=Yslice-X*betahat;    
temp = sum(resid.^2,1);
if sum(sum(temp==0))>0
    a=5
end
        cor_vol(pix,slice)  = (sum(resid(2:end,:).*resid(1:end-1,:),1)./sum(resid.^2,1))';
    end%if    
 end%slice
cor_vol = reshape(cor_vol,sqrt(npix),  sqrt(npix), numslices);



