input_file='g:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc';
frametimes=(0:119)*3;
slicetimes=[0.14 0.98 0.26 1.10 0.38 1.22 0.50 1.34 0.62 1.46 0.74 1.58 0.86];
eventid=kron(ones(10,1),[1; 2]);
eventimes=(0:19)'*18+9;
duration=ones(20,1)*9;
height=ones(20,1);
events=[eventid eventimes duration height] 
X_cache=fmridesign(frametimes,slicetimes,events);
exclude=[1 2];
n_poly=3;
fwhm_rho=15;
contrast=[1 -1 0 0 0 0];

d=fmris_read_image(input_file,0,0);
d.dim
numframes=d.dim(4);
numslices=d.dim(3);
numys=d.dim(2);
numxs=d.dim(1);
numpix=numxs*numys;
Steps=d.vox

allpts = 1:numframes;
allpts(exclude) = zeros(1,length(exclude));
keep = allpts( find( allpts ) );
n=length(keep)

    trend=((2*keep-(max(keep)+min(keep)))./(max(keep)-min(keep)))';
    Trend=(trend*ones(1,n_poly+1)).^(ones(n,1)*(0:n_poly));

numtrends=size(Trend,2)
numresponses=size(X_cache,2)


indk1=((keep(2:n)-keep(1:n-1))==1);
k1=find(indk1)+1;
fwhm_x=fwhm_rho/abs(Steps(1));
ker_x=exp(-(-ceil(fwhm_x):ceil(fwhm_x)).^2*4*log(2)/fwhm_x^2);
ker_x=ker_x/sum(ker_x);
fwhm_y=fwhm_rho/abs(Steps(2));
ker_y=exp(-(-ceil(fwhm_y):ceil(fwhm_y)).^2*4*log(2)/fwhm_y^2);
ker_y=ker_y/sum(ker_y);
ker_xy=ker_x'*ker_y;

slice=4;
X=[squeeze(X_cache(keep,:,1,slice)) Trend];
pinvX=pinv(X);

d=fmris_read_image(input_file,slice,1:numframes);
img=reshape(d.data,numpix,numframes);
Y=img(:,keep)';
betahat_ls=pinvX*Y;
resid=Y-X*betahat_ls;
i=64;
j=64;
p=rank(X);
df=n-1-p;
effect=zeros(numxs,numys);
sdeffect=zeros(numxs,numys);
for j=2:(numys-1)
    for i=2:(numxs-1)
        pix=i+(j-1)*numxs;
        nbr=[pix pix+1 pix-1 pix+numxs pix-numxs];
        %nbr=[pix];
        coef=pinv(resid(k1-1,nbr))*resid(k1,pix);
        Ystar=Y(k1,pix)-Y(k1-1,nbr)*coef;
        Xstar=X(k1,:)-X(k1-1,:)*sum(coef);
        pinvXstar=pinv(Xstar);
        betahat=pinvXstar*Ystar;
        sd=norm(Ystar-Xstar*betahat)/sqrt(df);
        effect(i,j)=contrast*betahat;
        V=pinvXstar*pinvXstar';
        sdeffect(i,j)=sqrt(diag(contrast*V*contrast'))*sd;
    end
end
tstat1=effect./(sdeffect+(sdeffect<=0)).*(sdeffect>0);
imagesc(tstat1'); axis xy; colorbar; colormap(spectral);

figure(1);
imagesc([tstat1 tstat0]'); axis xy; colorbar; colormap(spectral);