function roi_anat_generate(varargin)
% Generate anatomically defined ROIs, based on Talairach Daemon (TD)
% fields.
  
% Input arguments must be strings coming from the TD abbraviated field
% names, e.g. ('BA10','mFG','sFG'). The extent of the generated ROI is
% determined by, first, taking the union of strings within column and,
% second, taking the intersection bewteen (groups of) strings across columns.

  
if nargin<1, error('Insufficient number of input arguments...'), end
  
  load TDdb;
  
  % parse input
  % ------------------------------------------------------------------
  for i = 1:nargin
    
    argi{i} = zeros(5,1);
    tmp = strmatch(varargin{i},hem_l,  'exact');if tmp, argi{i}(1)=tmp; end
    tmp = strmatch(varargin{i},lobe_l, 'exact');if tmp, argi{i}(2)=tmp; end
    tmp = strmatch(varargin{i},gyrus_l,'exact');if tmp, argi{i}(3)=tmp; end
    tmp = strmatch(varargin{i},ttype_l,'exact');if tmp, argi{i}(4)=tmp; end
    tmp = strmatch(varargin{i},brod_l, 'exact');if tmp, argi{i}(5)=tmp; end
    clear tmp;

    % check to see that there is no trouble with the TD database
    % or that the input args are spelled wrong
    if length(find(argi{i})) < 1, 
      error(['There was no match for ',sprintf(varargin{i}),...
	     ' in the list of abbreviations. Possible typo.\n\n']);
    elseif length(find(argi{i})) > 1,
      error(['The same abbreviation was used for two entries. \nThis ', ...
	     'may be a bug in the ROI program - please report.\n\n']);
    end
    
  end

  % form union of all labels within field
  % ----------------------------------------------------------------
  
  indeces=[argi{:}]; % these are the consecutive (not the the assigned)
                     % indeces for the field variables
  tmp = find(indeces(1,:)); if tmp, roi_hem_i   = indeces(1,tmp); end
  tmp = find(indeces(2,:)); if tmp, roi_lobe_i  = indeces(2,tmp); end
  tmp = find(indeces(3,:)); if tmp, roi_gyrus_i = indeces(3,tmp); end
  tmp = find(indeces(4,:)); if tmp, roi_ttype_i = indeces(4,tmp); end
  tmp = find(indeces(5,:)); if tmp, roi_brod_i  = indeces(5,tmp); end
  
  
  % indeces from the field variables (in principle consecutive),
  % corresponding to which coordinates from XYZtd comprise the ROI
  % hem
  if exist('roi_hem_i','var')
    hem_idx =[]; hem_name=[];
    for i = 1: length(roi_hem_i)
      hem_idx  =sort([hem_idx; find(hem==hem_i(roi_hem_i(i)))]);
      hem_name =[hem_name,char(hem_l(roi_hem_i(i))),'-'];
    end
  end
  % lobe
  if exist('roi_lobe_i','var')
    lobe_idx =[]; lobe_name =[];
    for i = 1: length(roi_lobe_i)
      lobe_idx  =sort([lobe_idx; find(lobe==lobe_i(roi_lobe_i(i)))]);
      lobe_name =[lobe_name,char(lobe_l(roi_lobe_i(i))),'-'];
     end
  end
  % gyrus
  if exist('roi_gyrus_i','var')
    gyrus_idx =[];gyrus_name =[];
    for i = 1: length(roi_gyrus_i)
      gyrus_idx  =sort([gyrus_idx; find(gyrus==gyrus_i(roi_gyrus_i(i)))]);
      gyrus_name =[gyrus_name,char(gyrus_l(roi_gyrus_i(i))),'-'];
    end
  end
  % ttype
  if exist('roi_ttype_i','var')
    ttype_idx =[];ttype_name =[];
    for i = 1: length(roi_ttype_i)
      ttype_idx  =sort([ttype_idx; find(ttype==ttype_i(roi_ttype_i(i)))]);
      ttype_name =[ttype_name,char(ttype_l(roi_ttype_i(i))),'-'];
    end
  end
  % brod
  if exist('roi_brod_i','var')
    brod_idx =[];brod_name =[];
    for i = 1: length(roi_brod_i)
      brod_idx  =sort([brod_idx; find(brod==brod_i(roi_brod_i(i)))]);
      brod_name =[brod_name,char(brod_l(roi_brod_i(i))),'-'];
    end
  end
  

  % form the intersection across columns
  % -------------------------------------------------------------
  
  % generate a zero index matrix with one column for each field
  % variable and fill in with ones the selected cells in each column
  XYZidx = zeros(size(XYZtd,1),5);
  if exist('hem_idx'  ,'var'), XYZidx(hem_idx,  1)=1; end
  if exist('lobe_idx' ,'var'), XYZidx(lobe_idx, 2)=1; end
  if exist('gyrus_idx','var'), XYZidx(gyrus_idx,3)=1; end
  if exist('ttype_idx','var'), XYZidx(ttype_idx,4)=1; end
  if exist('brod_idx' ,'var'), XYZidx(brod_idx, 5)=1; end
  
  % compress across columns, leaving the  intersection
  XYZidx = sum(XYZidx,2);

  % fount the number of fields inputed
  num_fields = exist('hem_idx','var')+...
	       exist('lobe_idx','var')+...
	       exist('gyrus_idx','var')+...
	       exist('ttype_idx','var')+...
	       exist('brod_idx','var');
  
  % ROIidx is the index defining the roi in the 6 fields, incl.mm
  ROIidx = find(XYZidx==num_fields);
  
  
  % write out an image mask file 
  % ----------------------------------------------------------------
    
  % figure out file name for the roi image file
  % (insert '_' for intersection and '-' for union)
  roi_fname = [];
  if exist('hem_name','var')
    roi_fname = [roi_fname, hem_name(1:end-1), '_'];
  end
  if exist('ttype_name','var')
    roi_fname = [roi_fname, ttype_name(1:end-1), '_'];
  end
  if exist('lobe_name','var')
    roi_fname = [roi_fname, lobe_name(1:end-1), '_'];
  end
  if exist('brod_name','var')
    roi_fname = [roi_fname, brod_name(1:end-1), '_'];
  end
  if exist('gyrus_name','var')
    roi_fname = [roi_fname, gyrus_name(1:end-1), '_'];
  end
  roi_fname = roi_fname(1:end-1);    %to remove the last '_'
    

  % build output image file vol structure
  ROIvol = struct ...
	   ('fname',	roi_fname, ...
	    'dim',	[91 109 91 2], ...
	    'mat',	[2 0 0 -92; 0 2 0 -128; 0 0 2 -74; 0 0 0 1], ...
	    'pinfo',	[1;0;0], ...
	    'descrip',	'Tal-defined ROI mask');
  	    
  % find the XYZmm coordinates for this roi and
  % convert to pixel coordinates
  ROImm = XYZtd(ROIidx,:)';
  iM = inv(ROIvol.mat);
  XYZ = iM(1:3,:)*[ROImm; ones(1,size(ROImm,2))];
  
  % fill in a zeros image and then fill with ones the ROI region
  Y	= zeros(ROIvol.dim(1:3));
  idx	= sub2ind(ROIvol.dim(1:3),XYZ(1,:),XYZ(2,:),XYZ(3,:));
  Y(idx)= ones;
  
  % create the image file
  spm_write_vol(ROIvol,Y);
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  



