mars_img2rois

 creates ROIs from cluster image or image containing ROIs defined by unique nos
 FORMAT mars_img2rois(P, roipath, rootn, flags)

 P       - image (string or struct)
 roipath - path to directory to store ROIs
 rootn   - string to prefix to ROI filenames
 flags   - none or more of: [default = 'i']
             'i' - id image, voxel values identify ROIs
             'c' - cluster image, clusters identified by location
             'x' - label clusters by location of maximum 
                   (default is location of centre of mass)

 $Id$

0001 function mars_img2rois(P, roipath, rootn, flags)
0002 % creates ROIs from cluster image or image containing ROIs defined by unique nos
0003 % FORMAT mars_img2rois(P, roipath, rootn, flags)
0004 %
0005 % P       - image (string or struct)
0006 % roipath - path to directory to store ROIs
0007 % rootn   - string to prefix to ROI filenames
0008 % flags   - none or more of: [default = 'i']
0009 %             'i' - id image, voxel values identify ROIs
0010 %             'c' - cluster image, clusters identified by location
0011 %             'x' - label clusters by location of maximum
0012 %                   (default is location of centre of mass)
0013 %
0014 % $Id$
0015   
0016 if nargin < 1
0017   P = '';
0018 end
0019 if nargin < 2
0020   roipath = '';
0021 end
0022 if nargin < 3
0023   rootn = '';
0024 end
0025 if nargin < 4
0026   flags = ' ';
0027 end
0028 
0029 % Process input arguments
0030 if any(flags == 'i')
0031   Pprompt = 'Image containing ROI ids';
0032 else
0033   Pprompt = 'Image containing clusters';
0034 end
0035 if isempty(P)
0036   P = spm_get(1, mars_veropts('get_img_ext', Pprompt));
0037 end
0038 if isempty(P)
0039   return
0040 end
0041 if ischar(P)
0042   P = spm_vol(P);
0043 end
0044 if isempty(roipath)
0045   roipath = spm_get([-1 0], '', 'Directory to save ROIs');
0046 end
0047 if isempty(roipath)
0048   return
0049 end
0050 if isempty(rootn)
0051   [pn rootn ext] = fileparts(P.fname);
0052   rootn = spm_input('Prefix for ROI filenames', '+1', 's', rootn);
0053 end
0054 if isempty(rootn)
0055   return
0056 end
0057 
0058 if isempty(flags)
0059   flags = 'i';  % id image is default
0060 end
0061 
0062 % read img, get non-zero voxels
0063 img = spm_read_vols(P);
0064 img = img(:)';
0065 dim = P.dim(1:3);
0066 pts = find(img~=0);
0067 
0068 % e2xyz
0069 nz = pts-1;
0070 pl_sz = dim(1)*dim(2);
0071 Z = floor(nz / pl_sz);
0072 nz = nz - Z*pl_sz;
0073 Y = floor(nz / dim(1));
0074 X = nz - Y*dim(1);
0075 XYZ = [X; Y;Z] +1;
0076 
0077 % collect clusters
0078 vals = img(pts);
0079 
0080 % select cluster or id
0081 if any(flags == 'i')
0082   cl_vals = vals;
0083 else
0084   cl_vals = spm_clusters(XYZ);
0085 end
0086 
0087 for c = unique(cl_vals)
0088   % points for this region/cluster
0089   t_cl_is = find(cl_vals == c);
0090 
0091   % corresponding XYZ
0092   cXYZ = XYZ(:, t_cl_is);
0093 
0094   if ~isempty(cXYZ)
0095       % location label for cluster images
0096       if any(flags == 'c')
0097           if any(flags == 'x') % maximum
0098               [mx maxi] = max(vals(t_cl_is));
0099               mi = t_cl_is(maxi);
0100               % voxel coordinate of max
0101               vco = XYZ(:, mi);
0102       else % centre of mass
0103           vco = mean(cXYZ, 2);
0104       end
0105 
0106       % pt coordinates in mm
0107       pt_lab = P.mat * [vco; 1];
0108       pt_lab = pt_lab(1:3);
0109 
0110       % file name and labels
0111       d = sprintf('%s cluster at [%0.1f %0.1f %0.1f]', rootn, pt_lab);
0112       l = sprintf('%s_%0.0f_%0.0f_%0.0f', rootn, pt_lab);
0113 
0114   else % id image labels from voxel values
0115       % file name and labels
0116       d = sprintf('%s: id: %d', rootn, c);
0117       l = sprintf('%s_%d', rootn, c);
0118   end
0119 
0120   fname = maroi('filename', fullfile(roipath, l));
0121   o = maroi_pointlist(struct('XYZ',cXYZ,...
0122   'mat',P.mat,...
0123   'descrip',d,...
0124   'label', l), ...
0125   'vox');
0126   fprintf('\nSaving %s as %s...', d, fname);
0127   saveroi(o, fname);
0128   end
0129 end
0130 fprintf('\nDone...\n');