function  [p_vol,tstat_vol, effect_vol, sdeffect_vol] = cba_linmodel(rbrain, Xin, cor_vol,isbrain,contrast)
%rbrain = nxnxnumslicesxtc
%Xin = tcxnumhrfbasesxnumslices
%cor_vol = nxnxnumslices
%contrast = numcontrasts x numpredictors

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


%output:
tstat_vol = zeros(size(isbrain));
effect_vol = zeros(size(isbrain));
sdeffect_vol = zeros(size(isbrain));
p_vol  = zeros(size(isbrain));
for slice=1:size(rbrain ,3)
    slice
    for i=1:size(rbrain,1)
        for j=1:size(rbrain,2)
            if numresponses>1
    			X = zeros(size(Xin,1),size(Xin,2));
			    for k=1:size(Xin,2)
	               X(:,k)=whiten(Xin(:,k,slice),cor_vol(i,j,slice));%whitened design matrix (no whitening if cor is 0). 
               end
           else
                X=whiten(Xin(:,slice),cor_vol(i,j,slice));%whitened design matrix (no whitening if cor is 0). 
            end
            Df=size(rbrain,4)-rank(X);
            if (isbrain(i,j,slice)>0)
                Y =whiten((squeeze(rbrain(i,j,slice,:))),cor_vol(i,j,slice));%whitened data
                pinvX=pinv(X);
                betahat=pinvX*Y;
                resid=Y-X*betahat; 
                SSE=sum(resid.^2,1);
                sd=sqrt(SSE/Df);
                V=pinvX*pinvX';
                sdbetahat=sqrt(diag(V))*sd;
                mag_ef=contrast*betahat;%contr_mag*betahat;
                mag_sd=sqrt(diag(contrast*V*contrast'))*sd;
                effect_vol(i,j,slice)=mag_ef';
                sdeffect_vol(i,j,slice)=mag_sd';
                tstat_vol(i,j,slice)= (mag_ef./(mag_sd+(mag_sd<=0)).*(mag_sd>0))'; %UNDERSTAND THIS!!!
                p_vol(i,j,slice) = 1-tcdf(tstat_vol(i,j,slice),Df);
            end%if
        end %for j
    end%for i
end %for slice
                
            

