function [montage, cfg] = ft_nirs_prepare_ODtransformation(cfg, data)

% FT_NIRS_PREPARE_ODTRANSFORMATION returns the transformation matrix from
% optical densities (OD) to chromophore concentrations such as (de-)
% oxygenated hemoglobin.
%
% Use as
%   [montage] = ft_prepare_ODtransformation(cfg, data)
%
% It is necessary to input the data on which you want to perform the
% inverse computations, since that data generally contain the optode
% information and information about the channels that should be included in
% the transformation. The data structure can be either obtained
% from FT_PREPROCESSING, FT_FREQANALYSIS, FT_TIMELOCKANALYSIS or
% FT_COMPONENTANALYSIS. If the data is empty, all channels will be included
% in transformation.
%
% The configuration should contain
%   cfg.channel            = Nx1 cell-array with selection of channels
%                            (default = 'nirs'), see FT_CHANNELSELECTION for
%                            more details
%
% Optional configuration settings are
%   cfg.age                = scalar, age of the subject (necessary to
%                            automatically select the appropriate DPF, or
%   cfg.dpf                = scalar, differential path length factor
%   cfg.dpffile            = string, location to a lookup table for the
%                            relation between participant age and DPF
%
% Note that the DPF might be different across channels, and is usually
% contained in the optode structure contained in the data.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% If you specify this option the input data will be read from a *.mat
% file on disk. This mat files should contain only a single variable named 'data',
% corresponding to the input structure.
%
% See also FT_NIRS_TRANSFORM_ODS, FT_APPLY_MONTAGE

% options to be implemented:
%
% The NIRS positions can be present in the data or can be specified as
% cfg.opto      = structure with optode positions, see FT_DATATYPE_SENS
% cfg.siunits   = ft_getopt(cfg, 'siunits', 'no');   % yes/no, ensure that SI units are used consistently
% cfg.logarithm = string, can be 'natural' or 'base10' (default =

% You are using the FieldTrip NIRS toolbox developed and maintained by
% Artinis Medical Systems (http://www.artinis.com). For more information
% on FieldTrip, see http://www.fieldtriptoolbox.org
%
% This work is licensed under a Creative Commons Attribution-ShareAlike 4.0
% International License. To view a copy of this license, visit
% http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to
% Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
%
% Creative Commons Attribution-ShareAlike 4.0 International License:
% -----------------------------------
% You are free to:
%
%     Share - copy and redistribute the material in any medium or format
%     Adapt - remix, transform, and build upon the material
%     for any purpose, even commercially.
%
%     The licensor cannot revoke these freedoms as long as you follow the
%     license terms.
%
% Under the following terms:
%
%     Attribution - You must give appropriate credit, provide a link to
%                    the license, and indicate if changes were made. You
%                    may do so in any reasonable manner, but not in any way
%                    that suggests the licensor endorses you or your use.
%
%     ShareAlike - If you remix, transform, or build upon the material,
%                   you must distribute your contributions under the same
%                   license as the original.
%
%     No additional restrictions - You may not apply legal terms or
%                                   technological measures that legally
%                                   restrict others from doing anything the
%                                   license permits.
%
% -----------------------------------
%
% This toolbox is not to be used for medical or clinical purposes.
%
% Copyright (c) 2015-2016 by Artinis Medical Systems.
% Contact: askforinfo@artinis.com
%
% Main programmer: Jörn M. Horschig
% $Id$

revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function

% the ft_preamble function works by calling a number of scripts from
% fieldtrip/utility/private that are able to modify the local workspace

ft_defaults                   % this ensures that the path is correct and that the ft_defaults global variable is available
ft_preamble init              % this will reset ft_warning and show the function help if nargin==0 and return an error
ft_preamble debug             % this allows for displaying or saving the function name and input arguments upon an error
ft_preamble loadvar    data   % this reads the input data in case the user specified the cfg.inputfile option
ft_preamble provenance data   % this records the time and memory usage at the beginning of the function

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  % do not continue function execution in case the outputfile is present and the user indicated to keep it
  return
end

% check if the input data is valid for this function
data = ft_checkdata(data);

% set the defaults
cfg.channel = ft_getopt(cfg, 'channel', 'nirs');
cfg.age     = ft_getopt(cfg, 'age', []);
cfg.dpf     = ft_getopt(cfg, 'dpf', []);

% get the optode definition
opto = ft_fetch_sens(cfg, data);

% select the appropriate channels
if isfield(data, 'topolabel')
  % the data reflects a componentanalysis, where the topographic and the
  % timecourse labels are different
  cfg.channel = ft_channelselection(cfg.channel, data.topolabel);
elseif isfield(data, 'label')
  % In the subsequent code, the matching channels in the sensor array and
  % in the configuration will be selected. To ensure that these channels
  % are also present in the data, update the configuration to match the data.
  cfg.channel = ft_channelselection(cfg.channel, data.label);
else
  % update the selected channels based on the optode definition
  cfg.channel = ft_channelselection(cfg.channel, opto.label);
end
cfg.channel = ft_channelselection('nirs', cfg.channel); % only NIRS valid

% check current data type by looking what's inside of []
typeidx = regexp(cfg.channel, sprintf('%s%s', regexptranslate('wildcard', '[*]'), '$'));

% remove channels that do not obey to labelling guidelines
illegalidx = cellfun(@isempty, typeidx);
cfg.channel(illegalidx) = [];

if isempty(cfg.channel)
  ft_error('no valid NIRS channels found')
end

% channel indices wrt optode structure
chanidx     = match_str(opto.label, cfg.channel);

% check on DPF values
if ~isempty(cfg.age) && ~isempty(cfg.dpf)
  error('unsure if age or dpf in configuration should take precedence');
elseif ~isempty(cfg.age)
  error('the use of cfg.age is not implemented, yet');
elseif ~isempty(cfg.dpf)
  % this is where they should be
  dpfs = repmat(cfg.dpf, size(cfg.channel));
elseif isfield(opto, 'DPF')
  % this is for backward compatibility
  dpfs = opto.DPF(chanidx);
elseif isfield(data, 'opto') && isfield(data.opto, 'DPF')
  % this is for backward compatibility
  dpfs = data.opto.DPF(chanidx);
elseif isfield(data, 'hdr') && isfield(data.hdr, 'opto') && isfield(data.hdr.opto, 'DPF')
  % this is for backward compatibility
  dpfs = data.hdr.opto.DPF(chanidx);
end

% which chromophores are desired
chromophoreIdx = [3 2];
chromophoreName = {'O2Hb' 'HHb'};

%% transform to concentrations or to OD
% read in the coefficient table
% FIXME put this into an own function to read this out
fid = fopen(fullfile(fileparts(mfilename('fullpath')), 'private', 'Cope_ext_coeff_table.txt'));
coefs = cell2mat(textscan(fid, '%f %f %f %f %f'));

% extract all optode combinations that are relevant here
tratra         = opto.tra(chanidx, :)';
transmitteridx = tratra>0;
receiveridx    = tratra<0;
optodeidx      = (transmitteridx | receiveridx)'; % transpose to get back to tra order

% extract the wavelengths
wavelengths = opto.wavelength(tratra(transmitteridx));
wlidx = bsxfun(@minus, coefs(:, 1), wavelengths);

% find the relevant channel combinations
chancmb = logical([]); % we need to loop here to get size of chancmb
chanUsed = zeros(numel(chanidx), 1);
for c=1:numel(chanidx)
  % skip channel if it already was in another channel combination
  if chanUsed(c)
    continue;
  end
  % compute the channel combinations
  tupletidx = sum(bsxfun(@minus, optodeidx, optodeidx(c, :))~=0, 2)==0;
  chanUsed = chanUsed|tupletidx;
  chancmb(:, end+1) = tupletidx;
end

% do the transformation
tra      = zeros(size(chancmb, 2)*numel(chromophoreIdx), size(cfg.channel, 1));
labelnew = cell(size(chancmb, 2)*numel(chromophoreIdx), 1);

% transformation has to be done per channel combination
chanidx = []; % we will use this variable here for indexing channels (Rx-Tx cmbs)
optoidx = []; % we will use this variable here for indexing optodes (individual Rx and Tx)
for c=1:size(chancmb, 2)

  % find all other channels of the exact same transmitter / receiver pair
  chanidx = chancmb(:, c);

  % extract coefficient idx of these channels
  [coefidx, colidx] = find(wlidx(:, chanidx)==0);

  % compute the transmitter/receiver distance in cm
  optoidx = find(chanidx, 1, 'first'); % we can take 'first' because the transmitter-optodes have to be physically in the exact same spot to form "one" channel
  dist = sqrt(sum(diff(opto.optopos(optodeidx(optoidx, :), :)).^2));
    
  % select dpf
  dpf = mean(dpfs(chanidx));

  % compute distance factor for correct scaling
  distfactor = dist/10*dpf/10; % in mm

  % select the respective coefficients
  eps=coefs(coefidx,chromophoreIdx)/10; % in mm

  % add to transformation matrix
  traIdx = 1+(c-1)*numel(chromophoreIdx):(c)*numel(chromophoreIdx);
  tra(traIdx, chanidx) = pinv(eps)/(distfactor);

  % create the new labels for the channels
  chanorg = cfg.channel(find(chanidx, 1, 'first'));
  for n=1:numel(chromophoreIdx)
    labelnew(traIdx(n), 1) = regexprep(chanorg, sprintf('%s%s', regexptranslate('wildcard', '[*]'), '$'), sprintf('[%s]', chromophoreName{n}));
  end

end

%% create output
montage = [];
montage.labelorg = cfg.channel(:)';
montage.labelnew = labelnew;
montage.tra      = tra;
