function [sd,dist]=rbsdmap(cfg,varargin)
%
% [sd,dist]=rbsdmap(cfg,maxdist)
% [sd,dist]=rbsdmap(cfg,'maxdist',10,'excludesrc',[1,2],'excludedet',[1],'wavesrc',...,'wavedet',...)
%
% Create a source-detector mapping table (sd) by setting a maximum separation distance
%
% author: Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%     cfg: the redbird data structure
%     maxdist: the maximum source detector separation distance (in mm)
%
% output:
%     sd: the source detector mapping table
%     dist: the full table of all source detector separations (#row=#sources, #col=#detectors)
%
% license:
%     GPL version 3, see LICENSE_GPLv3.txt files for details 
%
% -- this function is part of Redbird-m toolbox
%

if((~isfield(cfg,'srcpos') && ~isfield(cfg,'widesrc')) || (isempty(cfg.srcpos) && ~isfield(cfg,'widesrc')) || (~isfield(cfg,'detpos') && ~isfield(cfg,'widedet')) || (isempty(cfg.detpos) && ~isfield(cfg,'widedet')))
    error('you must define at least 1 source and 1 detector in order to use this function');
end

if(length(varargin)==1)
    maxdist=varargin{1};
elseif(isempty(varargin))
    maxdist=inf;
elseif(length(varargin)>1)
    opt=varargin2struct(varargin{:});
    maxdist=jsonopt('maxdist',inf,opt);
end

if ~exist('opt','var')
    opt = struct();
end

srcnum=size(cfg.srcpos,1);
detnum=size(cfg.detpos,1);
if isfield(cfg,'widesrc')
    widesrcnum = size(cfg.widesrc,1);
    if isfield(cfg,'excludewfsrc')
        badwfsrc = cfg.excludewfsrc;
    elseif isfield(opt,'excludewfsrc')
        badwfsrc = jsonopt('excludewfsrc',[],opt);
    else
        badwfsrc = [];
    end
    badwfsrc = unique([badwfsrc find((sum(cfg.widesrc,2) == 0) | isnan(sum(cfg.widesrc,2)))']);
else
    widesrcnum = 0;
    badwfsrc = [];
end
if isfield(cfg,'widedet')
    widedetnum = size(cfg.widedet,1);
    if isfield(cfg,'excludewfdet')
        badwfdet = cfg.excludewfdet;
    elseif isfield(opt,'excludewfdet')
        badwfdet = jsonopt('excludewfdet',[],opt);
    else
        badwfdet = [];
    end
    badwfdet = unique([badwfdet find((sum(cfg.widedet,2) == 0) | isnan(sum(cfg.widedet,2)))']);
else
    widedetnum = 0;
    badwfdet = [];
end

if isfield(cfg,'excludesrc')
    badsrc = cfg.excludesrc;
else
    badsrc=jsonopt('excludesrc',[],opt);
end
if isfield(cfg,'excludedet')
    baddet = cfg.excludedet;
else
    baddet=jsonopt('excludedet',[],opt);
end


if(isfield(cfg,'srcpos') && ((size(cfg.srcpos,2) == size(cfg.face,1)) || (size(cfg.srcpos,2) == size(cfg.node,1))) )
    dist=zeros(srcnum,detnum);
else
    if isfield(cfg,'widesrc')
        widesrc = cfg.widesrc;
    else
        widesrc = [];
    end
    if isfield(cfg,'widedet')
        widedet = cfg.widedet;
    else
        widedet = [];
    end
    dist=rbgetdistance(cfg.srcpos,cfg.detpos,badsrc,baddet,widesrc,widedet,cfg)';
end

src = 1:srcnum;widesrc = [1:widesrcnum]+srcnum;
goodsrc=sort(setdiff(1:srcnum,badsrc));
goodwfsrc = sort(setdiff(1:widesrcnum,badwfsrc));
if ~isempty(goodwfsrc)
    goodwfsrc = goodwfsrc+srcnum;
end
det = 1:detnum;widedet = [1:widedetnum]+detnum;
gooddet=sort(setdiff(1:detnum,baddet));
goodwfdet = sort(setdiff(1:widedetnum,badwfdet));
if ~isempty(goodwfdet)
    goodwfdet = goodwfdet+detnum;
end

if(isfield(cfg,'prop') && isa(cfg.prop,'containers.Map'))
    wavelengths=cfg.prop.keys;
    sd=containers.Map();
    [ss,dd]=meshgrid([src widesrc],srcnum+widesrcnum+[det widedet]);

    for wv=wavelengths
        wid=wv{1};
        if (isfield(cfg,'wfsrcmapping'))
            wfsrcmap = cfg.wfsrcmapping(wid);
        end
        if (isfield(cfg,'wfdetmapping'))
            wfdetmap = cfg.wfdetmapping(wid);
        end
        
        if ((isfield(cfg,'wavesrc') && ~isempty(cfg.wavesrc)) || ...
           (isfield(cfg,'wavedet') && ~isempty(cfg.wavedet)))
             if(~isempty(cfg.wavesrc(wid)))
                 wavesrc = cfg.wavesrc(wid);
                 if exist('wfsrcmap','var')
                     wavesrc = rbremapsrc(wavesrc,wfsrcmap,srcnum);
                 end
                 goodwavesrc = intersect([goodsrc goodwfsrc],wavesrc);
             else
                 wavesrc = goodsrc;
             end
             if(~isempty(cfg.wavedet(wid)))
                 wavedet  = cfg.wavedet(wid);
                 if exist('wfdetmap','var')
                     wavedet = rbremapsrc(wavedet,wfdetmap,detnum);
                 end
                 goodwavedet = intersect([gooddet goodwfdet],wavedet);
             else
                wavedet = gooddet;
             end
             [ss,dd]=meshgrid(wavesrc,srcnum+widesrcnum+wavedet);
        else
            goodwavesrc = [goodsrc goodwfsrc];
            goodwavedet = [gooddet goodwfdet];
        end
        sdwv=[ss(:),dd(:)];
%         if(nargin<2 || (size(cfg.srcpos,2) == size(cfg.face,1)))
        if( ~isinf(maxdist) )
            [s2,d2] = meshgrid(1:srcnum,(1:detnum)+srcnum+widesrcnum);
            sd2 = [s2(:) d2(:)];
            [~,idx] = ismember(sdwv,sd2,'rows');idx = idx(find(idx));
            sdwv(:,3) = 1;
            sdwv(:,3)=(reshape(dist(unique(sdwv(:,2)) - (srcnum+widesrcnum),unique(sdwv(:,1))),[],1)<maxdist);
            %sdwv(idx,3)=(dist(:)<maxdist);
        else
            sdwv(:,3) = 0;
            sdwv(ismember(sdwv(:,1),goodwavesrc) & ismember(sdwv(:,2),goodwavedet+srcnum+widesrcnum),3) = 1;
        end
        
        if (isfield(cfg,'rfcw'))
            modes = cfg.rfcw.src.keys;
            sdwv(:,4) = zeros(size(sdwv,1),1);
            for md = modes
                mid = md{1};
                mdSRC = cfg.rfcw.src(mid);
                if (exist('wfsrcmap','var') && any(ismember(wfsrcmap(:,1),mdSRC)))
                    mdSRC = rbremapsrc(mdSRC,wfsrcmap,srcnum);
                end
                mdDET = cfg.rfcw.det(mid);
                if (exist('wfdetmap','var') && any(ismember(wfdetmap(:,1),mdDET)))
                    mdDET = rbremapsrc(mdDET,wfdetmap,detnum);
                end
                mdDET = mdDET + srcnum + widesrcnum;
                
                mdChan = (ismember(sdwv(:,1), mdSRC) & ismember(sdwv(:,2),mdDET));
                if strcmp(mid,'RF')
                    sdwv(mdChan,4) = sdwv(mdChan,4) + 1;
                elseif strcmp(mid,'CW')
                    sdwv(mdChan,4) = sdwv(mdChan,4) + 2;
                end
            end
            sdwv(sdwv(:,4)==0,:) = [];
%             sdwv((sdwv(:,3) == 0 | sdwv(:,4) == 0),:) = [];
        else
            sdwv((sdwv(:,3) == 0),:) = [];
        end
        sd(wid)=sdwv;
    end
else
    [ss,dd]=meshgrid([goodsrc goodwfsrc+srcnum],[gooddet goodwfdet+detnum]+(srcnum+widesrcnum));
    sd=[ss(:),dd(:)];
    if(nargin<2 || (size(cfg.srcpos,2) == size(cfg.face,1)))
        sd(:,3)=1;
    else
        sd(:,3)=(dist(:)<maxdist);
    end
end


function remap = rbremapsrc(sourcelist,wflist,srcnum)
if (isempty(intersect(sourcelist,wflist)) && ~isempty(wflist))
    remap = sourcelist;
else
    [idA,idB] = ismember(wflist(:,1),sourcelist);
    remap = sourcelist;
    remap(idB(find(idB))) = [];
    wf = wflist(idA,:);
    for ii = 1:size(wf,1)
        remap = [remap srcnum+(wf(ii,2):wf(ii,3))];
    end
end
