function sphere=spm_roi_sphere_mask(loc,radius,img)
% creates a mask containing a single sphere
% usage: mk_sphere_mask([-21 -24 0],8) creates a sphere with 8 mm radius
% centered at MNI coordinates [-21 -24 0]

  if ~exist('img'),
    img=spm_get(1,'*.img','choose any image for this subject');
  end;

  if exist(img,'file'),
    v=spm_vol(img);
  else,
    fprintf('%s not found\n',img);
    return;
  end;

  % now find the voxel that is nearest on the grid
  
  x=[v.mat(1,4):v.mat(1,1):v.dim(1)*v.mat(1,1)];
  y=[v.mat(2,4):v.mat(2,2):v.dim(2)*v.mat(2,2)];
  z=[v.mat(3,4):v.mat(3,3):v.dim(2)*v.mat(3,3)];
  if isempty(find(x==loc(1))) | isempty(find(y==loc(2))) | isempty(find(z==loc(3))),
    fprintf('the specified location is not on the grid!\n');
    return;
  end;

  if size(loc,1)<size(loc,2),loc=loc';end;
  
  % create a list of possible coordinates
  xfudge=round(2*radius/v.mat(1,1))*v.mat(1,1);
  yfudge=round(2*radius/v.mat(2,2))*v.mat(2,2);
  zfudge=round(2*radius/v.mat(3,3))*v.mat(3,3);
  
  XYZmm=[];
  c=1;
  for x=loc(1)-xfudge:v.mat(1,1):loc(1)+xfudge,
    for y=loc(2)-yfudge:v.mat(2,2):loc(2)+yfudge,
      for z=loc(3)-zfudge:v.mat(3,3):loc(3)+zfudge,
	fprintf('%d\t%d\t%d\n',x,y,z);
	XYZmm(:,c)=[x;y;z];
	c=c+1;
      end;
    end;
  end;

  j  = find(sum((XYZmm - loc*ones(1,size(XYZmm,2))).^2) <= radius^2);
whos  

  maskimg=zeros(v.dim(1:3));
  for i=j,
    coord=[];
    coord(1)=(XYZmm(1,i)-v.mat(1,4))/v.mat(1,1);
    coord(2)=(XYZmm(2,i)-v.mat(2,4))/v.mat(2,2);
    coord(3)=(XYZmm(3,i)-v.mat(3,4))/v.mat(3,3);
    maskimg(coord(1),coord(2),coord(3))=1;
  end;
  
  sphere=struct('j',j,'v',v,'maskimg',maskimg);
  
  v.fname=sprintf('sphmsk_%dmm_%d_%d_%d.img',radius,loc);
  spm_write_vol(v,maskimg);

