function [output] = ft_volumelookup(cfg, volume)

% FT_VOLUMELOOKUP can be used in to combine an anatomical or functional
% atlas with the source reconstruction results. You can use it for forward
% and reverse lookup.
%
% Given the region of interest (ROI) as anatomical or functional label, it
% looks up the locations and creates a mask (as a binary volume) based on
% the label. Given the ROI as point in the brain, it creates a sphere or
% box around that point. In these two case the function is to be used as:
%   mask = ft_volumelookup(cfg, volume)
%
% Given a binary volume that indicates a ROI or a point of interest (POI),
% it looks up the corresponding anatomical or functional labels from the
% atlas. In this case the function is to be used as:
%    labels = ft_volumelookup(cfg, volume)
%
% In both cases the input volume can be:
%   mri    is the output of FT_READ_MRI source is the output of FT_SOURCEANALYSIS
%   stat   is the output of FT_SOURCESTATISTICS
%
% The configuration options for a mask according to an atlas:
%   cfg.atlas               = string, filename of atlas to use, see FT_READ_ATLAS
%   cfg.roi                 = string or cell-array of strings, ROI from anatomical atlas
%
% The configuration options for a spherical/box mask around a POI:
%   cfg.roi                 = Nx3 vector, coordinates of the POI
%   cfg.sphere              = radius of each sphere in cm/mm dep on unit of input
%   cfg.box                 = Nx3 vector, size of each box in cm/mm dep on unit of input
%   cfg.round2nearestvoxel  = 'yes' or 'no' (default = 'no'), voxel closest to point of interest is calculated
%                             and box/sphere is centered around coordinates of that voxel
%
% The configuration options for labels from a mask:
%   cfg.atlas               = string, filename of atlas to use, see FT_READ_ATLAS
%   cfg.maskparameter       = string, field in volume to be looked up, data in field should be logical
%   cfg.minqueryrange       = number, should be odd and <= to maxqueryrange (default = 1)
%   cfg.maxqueryrange       = number, should be odd and >= to minqueryrange (default = 1)
%
% The configuration options for labels around POI:
%   cfg.output              = 'single' always outputs one label; if several POI are provided, they are considered together as describing a ROI (default)
%                             'multiple' outputs one label per POI (e.g., choose to get labels for different electrodes)
%   cfg.roi                 = Nx3 vector, coordinates of the POI
%   cfg.atlas               = string, filename of atlas to use, see FT_READ_ATLAS
%   cfg.minqueryrange       = number, should be odd and <= to maxqueryrange (default = 1)
%   cfg.maxqueryrange       = number, should be odd and >= to minqueryrange (default = 1)
%   cfg.querymethod         = 'sphere' searches voxels around the ROI in a sphere (default)
%                           = 'cube' searches voxels around the ROI in a cube
%   cfg.round2nearestvoxel  = 'yes' or 'no', voxel closest to POI is calculated (default = 'yes')
%
% The label output has a field "names", a field "count" and a field "usedqueryrange".
% To get a list of areas of the given mask you can do for instance:
%      [tmp ind] = sort(labels.count,1,'descend');
%      sel = find(tmp);
%      for j = 1:length(sel)
%        found_areas{j,1} = [num2str(labels.count(ind(j))) ': ' labels.name{ind(j)}];
%      end
% In the "found_areas" variable you can then see how many times which labels are
% found. Note that in the AFNI brick one location can have 2 labels.
%
% Dependent on the input coordinates and the coordinates of the atlas, the
% input MRI is transformed betweem MNI and Talairach-Tournoux coordinates
% See http://www.mrc-cbu.cam.ac.uk/Imaging/Common/mnispace.shtml for more details.
%
% See http://www.fieldtriptoolbox.org/template/atlas for a list of templates and
% atlasses that are included in the FieldTrip release.
%
% See also FT_READ_ATLAS, FT_SOURCEPLOT

% Copyright (C) 2008-2017, Robert Oostenveld, Ingrid Nieuwenhuis
% Copyright (C) 2013, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar volume
ft_preamble provenance volume

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input data is valid for this function
% additional checking of the input data is done further down
volume = ft_checkdata(volume, 'datatype', {'volume', 'source'}, 'feedback', 'yes', 'hasunit', 'yes', 'hascoordsys', 'yes');

% the handling of the default cfg options is done further down
cfg.minqueryrange      = ft_getopt(cfg,'minqueryrange', 1);
cfg.maxqueryrange      = ft_getopt(cfg,'maxqueryrange', 1);
cfg.output             = ft_getopt(cfg,'output', []); % in future, cfg.output could be extended to support both 'label' and 'mask'

% ensure that old and unsupported options are not being relied on by the end-user's script
% instead of specifying cfg.coordsys, the user should specify the coordsys in the data
cfg = ft_checkconfig(cfg, 'forbidden', {'units', 'coordsys', 'inputcoord', 'inputcoordsys', 'coordinates'});
cfg = ft_checkconfig(cfg, 'renamedval', {'output', 'label', 'single'});

roi2mask   = 0;
mask2label = 0;
roi2label  = 0;
if isfield(cfg, 'roi') && ~isempty(cfg.output)
  roi2label = 1;
elseif isfield(cfg, 'roi')
  roi2mask = 1;
elseif isfield(cfg, 'maskparameter')
  mask2label = 1;
else
  ft_error('you should either specify cfg.roi, or cfg.maskparameter')
end

if roi2mask
  % only for volume data
  volume = ft_checkdata(volume, 'datatype', 'volume');

  cfg.round2nearestvoxel = ft_getopt(cfg, 'round2nearestvoxel', 'no');

  isatlas = iscell(cfg.roi) || ischar(cfg.roi);
  ispoi   = isnumeric(cfg.roi);
  if isatlas+ispoi ~= 1
    ft_error('do not understand cfg.roi')
  end

  if isatlas
    ft_checkconfig(cfg, 'forbidden', {'sphere', 'box'}, 'required',  {'atlas'});
  elseif ispoi
    ft_checkconfig(cfg, 'forbidden', {'atlas'});
    if isempty(ft_getopt(cfg, 'sphere')) && isempty(ft_getopt(cfg, 'box'))
      ft_error('you should either specify cfg.sphere or cfg.box')
    end
  end

elseif mask2label || roi2label
  % convert to source representation, this is easier to work with
  volume = ft_checkdata(volume, 'datatype', 'source');
  ft_checkconfig(cfg, 'required', {'atlas'});

  if cfg.minqueryrange > cfg.maxqueryrange
    ft_error('maxqueryrange should be superior or equal to minqueryrange');
  end

  if rem(cfg.minqueryrange, 2) == 0 || rem(cfg.maxqueryrange, 2) == 0
    ft_error('incorrect query range, should be odd numbers');
  end

  cfg.round2nearestvoxel = ft_getopt(cfg, 'round2nearestvoxel', 'yes');
  cfg.querymethod        = ft_getopt(cfg, 'querymethod', 'sphere');
end


if roi2mask
  % determine location of each anatomical voxel in its own voxel coordinates
  dim = volume.dim;
  i = 1:dim(1);
  j = 1:dim(2);
  k = 1:dim(3);
  [I, J, K] = ndgrid(i, j, k);
  ijk = [I(:) J(:) K(:) ones(prod(dim),1)]';
  % determine location of each anatomical voxel in head coordinates
  xyz = volume.transform * ijk; % note that this is 4xN

  if isatlas
    if ischar(cfg.atlas)
      % assume it to represent a filename
      atlas = ft_read_atlas(cfg.atlas);
    else
      % assume cfg.atlas to be a struct
      atlas = cfg.atlas;
    end

    % determine which field(s) to use to look up the labels,
    % and whether these are boolean or indexed
    fn = fieldnames(atlas);
    isboolean = false(numel(fn),1);
    isindexed = false(numel(fn),1);
    for i=1:length(fn)
      if islogical(atlas.(fn{i})) && isequal(size(atlas.(fn{i})), atlas.dim)
        isboolean(i) = true;
      elseif isnumeric(atlas.(fn{i})) && isequal(size(atlas.(fn{i})), atlas.dim)
        isindexed(i) = true;
      end
    end
    if any(isindexed)
      % let the indexed prevail
      fn = fn(isindexed);
      isindexed = 1;
    elseif any(isboolean)
      % use the boolean
      fn = fn(isboolean);
      isindexed = 0;
    end

    if ischar(cfg.roi)
      cfg.roi = {cfg.roi};
    end

    if isindexed
      sel = zeros(0,2);
      for m = 1:length(fn)
        for i = 1:length(cfg.roi)
          tmp = find(strcmp(cfg.roi{i}, atlas.([fn{m},'label'])));
          sel = [sel; tmp m*ones(numel(tmp),1)];
        end
      end
      fprintf('found %d matching anatomical labels\n', size(sel,1));

      % this is to accommodate for multiple parcellations:
      % the brick refers to the parcellationname
      % the value refers to the value within the given parcellation
      brick = sel(:,2);
      value = sel(:,1);

      % convert between MNI head coordinates and TAL head coordinates
      % coordinates should be expressed compatible with the atlas
      if     strcmp(volume.coordsys, 'mni') && strcmp(atlas.coordsys, 'tal')
        xyz(1:3,:) = mni2tal(xyz(1:3,:));
      elseif strcmp(volume.coordsys, 'mni') && strcmp(atlas.coordsys, 'mni')
        % nothing to do
      elseif strcmp(volume.coordsys, 'tal') && strcmp(atlas.coordsys, 'tal')
        % nothing to do
      elseif strcmp(volume.coordsys, 'tal') && strcmp(atlas.coordsys, 'mni')
        xyz(1:3,:) = tal2mni(xyz(1:3,:));
      elseif ~strcmp(volume.coordsys, atlas.coordsys)
        ft_error('the mismatch between the coordinate system in the atlas and the coordinate system in the data cannot be resolved');
      end

      % determine location of each anatomical voxel in atlas voxel coordinates
      ijk = atlas.transform \ xyz;
      ijk = round(ijk(1:3,:))';

      inside_vol = ijk(:,1)>=1 & ijk(:,1)<=atlas.dim(1) & ...
        ijk(:,2)>=1 & ijk(:,2)<=atlas.dim(2) & ...
        ijk(:,3)>=1 & ijk(:,3)<=atlas.dim(3);
      inside_vol = find(inside_vol);

      % convert the selection inside the atlas volume into linear indices
      ind = sub2ind(atlas.dim, ijk(inside_vol,1), ijk(inside_vol,2), ijk(inside_vol,3));

      brick_val = cell(1,numel(brick));
      % search the bricks for the value of each voxel
      for i=1:numel(brick_val)
        brick_val{i} = zeros(prod(dim),1);
        brick_val{i}(inside_vol) = atlas.(fn{brick(i)})(ind);
      end

      mask = zeros(prod(dim),1);
      for i=1:numel(brick_val)
        %fprintf('constructing mask for %s\n', atlas.descr.name{sel(i)});
        mask = mask | (brick_val{i}==value(i));
      end
    else
      ft_error('support for atlases that have a probabilistic segmentationstyle is not supported yet');
      % NOTE: this may be very straightforward indeed: the mask is just the logical "or" of the specified ROIs.
    end

  elseif ispoi

    if istrue(cfg.round2nearestvoxel)
      for i=1:size(cfg.roi,1)
        cfg.roi(i,:) = poi2voi(cfg.roi(i,:), xyz);
      end
    end

    % sphere(s)
    if isfield(cfg, 'sphere')
      mask = zeros(1,prod(dim));
      for i=1:size(cfg.roi,1)
        dist = sqrt( (xyz(1,:) - cfg.roi(i,1)).^2 + (xyz(2,:) - cfg.roi(i,2)).^2 + (xyz(3,:) - cfg.roi(i,3)).^2 );
        mask = mask | (dist <= cfg.sphere(i));
      end
      % box(es)
    elseif isfield(cfg, 'box')
      mask = zeros(1, prod(dim));
      for i=1:size(cfg.roi,1)
        mask = mask | ...
          (xyz(1,:) <= (cfg.roi(i,1) + cfg.box(i,1)./2) & xyz(1,:) >= (cfg.roi(i,1) - cfg.box(i,1)./2)) & ...
          (xyz(2,:) <= (cfg.roi(i,2) + cfg.box(i,2)./2) & xyz(2,:) >= (cfg.roi(i,2) - cfg.box(i,2)./2)) & ...
          (xyz(3,:) <= (cfg.roi(i,3) + cfg.box(i,3)./2) & xyz(3,:) >= (cfg.roi(i,3) - cfg.box(i,3)./2));
      end
    end
  end

  mask = reshape(mask, dim);
  fprintf('%i voxels in mask, which is %.3f %% of total volume\n', sum(mask(:)), 100*mean(mask(:)));
  output = mask;

elseif mask2label || roi2label
  if ischar(cfg.atlas)
    % assume it to represent a filename
    atlas = ft_read_atlas(cfg.atlas);
  else
    % assume cfg.atlas to be a struct
    atlas = cfg.atlas;
  end

  % determine which field(s) to use to look up the labels,
  % and whether these are boolean or indexed
  fn = fieldnames(atlas);
  isboolean = false(numel(fn),1);
  isindexed = false(numel(fn),1);
  for i=1:length(fn)
    if islogical(atlas.(fn{i})) && isequal(size(atlas.(fn{i})), atlas.dim)
      isboolean(i) = true;
    elseif isnumeric(atlas.(fn{i})) && isequal(size(atlas.(fn{i})), atlas.dim)
      isindexed(i) = true;
    end
  end
  if any(isindexed)
    % let the indexed prevail
    fn = fn(isindexed);
    isindexed = 1;
  elseif any(isboolean)
    % use the boolean
    fn = fn(isboolean);
    isindexed = 0;
  end

  labels.name = cell(0,1);
  for k = 1:numel(fn)
    % ensure that they are concatenated as column
    tmp = atlas.([fn{k},'label']);
    labels.name = cat(1, labels.name(:), tmp(:));
  end
  labels.name{end+1} = 'no_label_found';
  labels.count = zeros(length(labels.name),1);
  for iLab = 1:length(labels.name)
    labels.usedqueryrange{iLab} = [];
  end

  if mask2label
    sel = find(volume.(cfg.maskparameter)(:));
  elseif roi2label
    if istrue(cfg.round2nearestvoxel)
      % determine location of each anatomical voxel in head coordinates
      xyz = [volume.pos ones(size(volume.pos,1),1)]'; % note that this is 4xN
      nSel = size(cfg.roi, 1);
      for i=1:nSel
        cfg.roi(i,:) = poi2voi(cfg.roi(i,:), xyz);
      end
    end % round2nearestvoxel
    sel = zeros(size(cfg.roi,1), 1);
    for i = 1:size(cfg.roi,1)
        sel(i) = find(volume.pos(:, 1) == cfg.roi(i, 1) & volume.pos(:, 2) == cfg.roi(i, 2) & volume.pos(:, 3) == cfg.roi(i, 3));
    end
  end
  if strcmp(cfg.output, 'multiple')
      labels = repmat(labels, length(sel), 1);
  end
  for iVox = 1:length(sel)
      label = {};
      for qr = cfg.minqueryrange:2:cfg.maxqueryrange
          if isempty(label)
              label = atlas_lookup(atlas, volume.pos(sel(iVox), :), 'coordsys', volume.coordsys, 'queryrange', qr, 'method', cfg.querymethod);
              usedQR = qr;
          end
      end

      if isempty(label)
          label = {'no_label_found'};
      elseif length(label) == 1
          label = {label};
      end

      if strcmp(cfg.output, 'multiple')
          iLabOut = iVox;
      else
          iLabOut = 1;
      end
      ind_lab = [];
      for iLab = 1:length(label)
          ind_lab = find(strcmp(label{iLab}, labels(iLabOut).name));
          labels(iLabOut).count(ind_lab) = labels(iLabOut).count(ind_lab)+1; % labels.count should give the number of times a label was found within a query range
      end

      %     labels.count(ind_lab) = labels.count(ind_lab) + (1/length(ind_lab));
      %     ^this gives each label a weight depending on the number of
      %     labels found within the query range. Using this method, all labels
      %     that were found will have the same number listed for labels.count,
      %     which defeats the point of the labels.count field as I understand it
      for iFoundLab = 1:length(ind_lab)
          if isempty(labels(iLabOut).usedqueryrange{ind_lab(iFoundLab)})
              labels(iLabOut).usedqueryrange{ind_lab(iFoundLab)} = usedQR;
          else
              labels(iLabOut).usedqueryrange{ind_lab(iFoundLab)} = [labels(iLabOut).usedqueryrange{ind_lab(iFoundLab)} usedQR];
          end
      end
  end %iVox

  output = labels;

end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble provenance


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION point of interest to voxel of interest
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function voi = poi2voi(poi, xyz)
xmin = min(abs(xyz(1,:) - poi(1))); xcl = round(abs(xyz(1,:) - poi(1))) == round(xmin);
ymin = min(abs(xyz(2,:) - poi(2))); ycl = round(abs(xyz(2,:) - poi(2))) == round(ymin);
zmin = min(abs(xyz(3,:) - poi(3))); zcl = round(abs(xyz(3,:) - poi(3))) == round(zmin);
xyzcls = xcl + ycl + zcl; ind_voi = xyzcls == 3;
if sum(ind_voi) > 1
  fprintf('%i voxels at same distance of poi, taking first voxel\n', sum(ind_voi))
  ind_voi_temp = find(ind_voi); ind_voi_temp = ind_voi_temp(1);
  ind_voi = zeros(size(ind_voi));
  ind_voi(ind_voi_temp) = 1;
  ind_voi = logical(ind_voi);
end
voi = xyz(1:3,ind_voi);
fprintf('coordinates of voi: %.1f  %.1f %.1f\n', voi(1), voi(2), voi(3));
