dir_current=pwd;
nsub=input('type in the number of subjects you were averaging :');


est_x=0.0;
est_y=0.0;
est_z=0.0;

xdim=79;
ydim=95;
zdim=34;

QQ=zeros(xdim,ydim,zdim);

for sub_idx=1:nsub

P = input(['type in the result directory for subject ' num2str(sub_idx) ' ' ],'s');
eval(['cd 'P]);

load SPM;

ss(sub_idx,:)=W;
end

for itr=1:10


for sub_idx=1:nsub


s=ss(sub_idx,:);


x  = round(6*s(1)); x = [-x:x];
y  = round(6*s(2)); y = [-y:y];
z  = round(6*s(3)); z = [-z:z];
x  = exp(-(x).^2/(2*(s(1)).^2)); 
y  = exp(-(y).^2/(2*(s(2)).^2)); 
z  = exp(-(z).^2/(2*(s(3)).^2));
x  = x/sum(x);
y  = y/sum(y);
z  = z/sum(z);

i  = (length(x) - 1)/2;
j  = (length(y) - 1)/2;
k  = (length(z) - 1)/2;

P=randn(xdim,ydim,zdim);
Q=ones(xdim,ydim,zdim);
Vi=[xdim,ydim,zdim,2.0,2.0,4.0,1,3,0,xdim/2,ydim/2,zdim/2]';

spm_conv_vol(Vi,Q,x,y,z,(-[i,j,k]),P);

QQ=QQ+Q;
end

Q=QQ./nsub;
v_vol=std(reshape(Q,xdim*ydim*zdim,1))^2.0;

dx=zeros(xdim,ydim,zdim);
for j=1:ydim
for k=1:zdim

dx(1:xdim-1,j,k)=Q(2:xdim,j,k)-Q(1:xdim-1,j,k);

end
end
;sx=mean(mean(mean(dx~=0)));
tmp=reshape(dx(1:xdim-1,:,:),(xdim-1)*ydim*zdim,1);
vxx=std(tmp)*std(tmp);

s_xx=sqrt(-1.0/(4.0*log(1.0-vxx/(2.0*v_vol))));

dy=zeros(xdim,ydim,zdim);
for i=1:xdim
for k=1:zdim

dy(i,1:ydim-1,k)=Q(i,2:ydim,k)-Q(i,1:ydim-1,k);

end
end
;sy=mean(mean(mean(dy~=0)));
tmp=reshape(dy(:,1:ydim-1,:),xdim*(ydim-1)*zdim,1);
vyy=std(tmp)*std(tmp);

s_yy=sqrt(-1.0/(4.0*log(1.0-vyy/(2.0*v_vol))));

dz=zeros(xdim,ydim,zdim);
for i=1:xdim
for j=1:ydim

dz(i,j,1:zdim-1)=Q(i,j,2:zdim)-Q(i,j,1:zdim-1);

end
end
;sz=mean(mean(mean(dz~=0)));
tmp=reshape(dz(:,:,1:zdim-1),xdim*ydim*(zdim-1),1);
vzz=std(tmp)*std(tmp);

s_zz=sqrt(-1.0/(4.0*log(1.0-vzz/(2.0*v_vol))));

est_x=est_x+s_xx;
est_y=est_y+s_yy;
est_z=est_z+s_zz;



end

est_x=est_x/10.;
est_y=est_y/10.;
est_z=est_z/10.;

W=[est_x est_y est_z]
eval(['cd 'dir_current]);
save SPM B C CONTRAST Dnames Fdf Fnames G H RT S SIGMA TH UF V W df -V4


