function [ROI]=spm_roi_define(radius,Finter)
% extracts all of the clusters of activation
% creates a .img file for each cluster, with zeros everywhere outside roi
%=======================================================================
% Russ Poldrack, 5/15/00
% based on roi_define.m by Kalina Christoff / 2000-01-24

verbose=0;
if ~exist('radius'),
   radius=8;
   fprintf('Using default radius of %d mm\n',radius);
else
   fprintf('Using radius of %d mm\n',radius);
end;
extent=radius;
cluster_distance=radius;

%-------------------------
% load the desired statistics file

%[SPM,VOL,xX,xCon,xSDM]=spm_getSPM;
 [hReg,SPM,VOL,xX,xCon,xSDM] = spm_results_ui;

mask_info=strcat(SPM.swd,'/mask.img');
if exist(mask_info),
  [dim vx scale typ offset origin descrip] = spm_hread(mask_info);
else,
  fprintf('Could not find mask.img, please choose one...\n');
  mask_info   = spm_get(1,'mask.img','select mask.img');
  [dim vx scale typ offset origin descrip] = spm_hread(mask_info);
end;


%-Characterize excursion set in terms of maxima
% (sorted on Z values and grouped by regions)
%==================================================================

TabDat= spm_list('List',SPM,VOL,cluster_distance,Inf,SPM.swd,hReg);



% using indexs in AA, get indexes for each roi Q.roi#
%--------------------------------------------------
Y=[];
y=[];
% need to deal with the difference in TabDat structure arising from
% use with FDR toolbox from Tom Nichols

location_column=size(TabDat.dat,2); 

for cluster=1:size(TabDat.dat,1)
     point=TabDat.dat{cluster,location_column};
%     point=vox2tal(point,VOL);
     fprintf('Checking Max #%d: ',cluster);
     sphere(cluster)=find_sphere(point,radius,SPM,VOL);
     numvox=size(sphere(cluster).XYZ,2);

      %-Print cluster info and ask for name of ROI
      roi_name=sprintf('con%d_clust%d',SPM.Ic(1),cluster);
        
      if roi_name
           
            % write a .tal text file with list of coordinates
            name_tal=[roi_name '.tal'];
            fid=fopen(name_tal,'w');
	    
	    % changed 10/05/02 to print out mm rather than voxel coords	    
            fprintf(fid,'%2d %2d %2d\n', sphere.XYZmm);
            fclose(fid);

            % write an image mask (.img and .hdr files) for the roi
            V.dim=[VOL.DIM(1) VOL.DIM(2) VOL.DIM(3) 4]; % write as unsigned int file
            V.fname=sprintf('%s.img',roi_name);
            V.mat=VOL.M;
            V.pinfo=[1;0;0];
            V.descrip=sprintf('funcROI: %d %d %d',point(1),point(2),point(3));
            
            roi_img=zeros(VOL.DIM(1),VOL.DIM(2),VOL.DIM(3));
            numvox=size(sphere(cluster).XYZ,2);
            for vox=1:numvox,
             %  fprintf('setting %d to 1000\n',roi_vox(:,vox));
               if sphere(cluster).Z > 0,
		 sphere(cluster).Z=TabDat.dat{cluster,7};
	       else, 
		 sphere(cluster).Z=100;
	       end;
               sphere(cluster).clust_p_corr=TabDat.dat{cluster,3};
               sphere(cluster).vox_p_corr=TabDat.dat{cluster,6};
               sphere(cluster).full_clust_size=TabDat.dat{cluster,4};
       
               roi_img(sphere(cluster).XYZ(1,vox),sphere(cluster).XYZ(2,vox),sphere(cluster).XYZ(3,vox))=sphere(cluster).Z;
            end
            spm_write_vol(V,roi_img);
 
       end
  
  

end

% print TabDat list to file - this is from spm_roi_list.m

point=zeros(1,3);
info=sprintf('ROI information for %s: Contrast %d\n',pwd,SPM.Ic(1));
info=sprintf('%s\nCluster #\tX\tY\tZ\tMax T\tNumVox\tPvox\tPclust\tClustSize',info);
for cluster=1:size(TabDat.dat,1),
       point=sphere(cluster).center;
       info=sprintf('%s\n%d:\t\t%d\t%d\t%d\t%3.2f\t%d\t%0.3f\t%0.3f\t%d',...
         info,cluster,point,max(sphere(cluster).Z),...
	 size(sphere(cluster).XYZ,2),sphere(cluster).vox_p_corr,...
	 sphere(cluster).clust_p_corr,sphere(cluster).full_clust_size);
end;
info=sprintf('%s\n',info);
roilistfile=sprintf('con%d_clust_list.txt',SPM.Ic(1));

fid=fopen(roilistfile,'w');
fprintf(fid,'%s',info);
fclose(fid);


ROI=struct(...
'SPM',     SPM,...
'VOL',     VOL,...
'xX',      xX,...
'xCon',    xCon,...
'xSDM',    xSDM,... 
 'sphere', sphere,...
'origin',  origin ...   
); 

save_cmd=sprintf('save con%d_roi ROI',SPM.Ic(1));
eval(save_cmd);
fprintf('spm_roi_define - DONE\n\n');
%clear













