function multiplotTFR(cfg, data)

% multiplotTFR plots time-frequency representations in a toporgaphical layout.
% The plots of the indivual sensors are arranged according to their location
% specified in the layout file.
%
% Use as:
%   multiplotTFR(cfg, data)
%
% where freq is the result from FREQANALYSIS or FREQDESCRIPTIVES and the cfg structure
% contains all configurations details for the plotting:
%
% cfg.layout           = specification of the layout, see below
% cfg.baseline         = 'yes','no' or [time1 time2]      (default = 'no')
% cfg.baselinetype     = 'absolute' or 'relative'         (default = 'absolute')
% cfg.xlim             = 'maxmin' or [xmin xmax]          (default = 'maxmin')
% cfg.ylim             = 'maxmin' or [ymin ymax]          (default = 'maxmin')
% cfg.zlim             = 'maxmin','absmax' or [zmin zmax] (default = 'maxmin')
% cfg.colorbar         = 'yes', 'no'                      (default = 'no')
% cfg.comment          = string of text                   (default = date + zlimits)
%                        Add 'comment' to graph (according to COMNT in the layout)
% cfg.showlabels       = 'yes', 'no', 'numbers'           (default = 'no')
% cfg.fontsize         = fontsize of text                 (default = 8)
% cfg.interactive      = 'yes', 'no'                      (default = 'no')
% cfg.yparam           = 'powspctrm' (default) or 'cohspctrm'
%
% In case cfg.yparam   = 'cohspctrm', you should specify the reference-channel in
% cfg.cohrefchannel    = string
%
% The layout defines how the channels are arranged and what the size of each
% subplot is. You can specify the layout in a variety of ways:
%  - you can give the name of an ascii layout file with extension *.lay
%  - you can give the name of an electrode file
%  - you can give an electrode definition, i.e. "elec" structure
%  - you can give a gradiometer definition, i.e. "grad" structure
% If you do not specify any of these and the data structure contains an
% electrode or gradiometer structure (common for MEG data, since the header
% of the MEG datafile contains the gradiometer information), that will be
% used for creating a layout. If you want to have more fine-grained control
% over the layout of the subplots, you should create your own layout file.
%
% The old way to plot the coherence with respect to a reference channel is
% still supported. Then specify
% cfg.yparam           = 'powspctrm';
% cfg.cohtargetchannel = string with the channel relative to which coherence
%                        is plotted (default is to plot power)
%
% See also:
%   multiplotER, singleplotER, singleplotTFR, topoplotER, topoplotTFR.
%
% Copyright (C) 2003, Ole Jensen
%
% $Log: multiplotTFR.m,v $
% Revision 1.14  2006/03/06 11:47:55  denpas
% Fixed zmin/zmax bug.
%
% Revision 1.13  2006/02/27 15:03:03  denpas
% many changes, most important is added interactive functionality
% made data selection consistent between different plot functions
% changed dimord for consistency
%
% Revision 1.9  2005/05/17 17:50:37  roboos
% changed all "if" occurences of & and | into && and ||
% this makes the code more compatible with Octave and also seems to be in closer correspondence with Matlab documentation on shortcircuited evaluation of sequential boolean constructs
%
% Revision 1.8  2005/04/08 17:26:54  olejen
% TFR was wrongly assigned before calling freqbaseline.m - this is now reversed
%
% Revision 1.7  2005/04/06 14:26:55  olejen
% clf before plotting
% now make use of freqbaseline.m
%
% Revision 1.6  2005/04/06 07:40:56  jansch
% included option cfg.cohrefchannel. updated help.
%
% Revision 1.5  2005/02/07 17:12:00  roboos
% changed handling of layout files (using new function createlayout), now also supports automatic layout creation based on gradiometer/electrode definition in data, updated help, cleaned up indentation
%
% Revision 1.4  2004/09/24 15:54:54  roboos
% included the suggested improvements by Doug Davidson: added option cfg.cohtargetchannel
% and updated the help
%
% Revision 1.3  2004/09/01 18:02:23  roboos
% added copyright statements, removed cfg as output argument
%

clf

% for backward compatibility with old data structures
data = fixdimord(data);

% Set the defaults:
if ~isfield(cfg,'baseline'),        cfg.baseline = 'no';               end
if ~isfield(cfg,'baselinetype'),    cfg.baselinetype = 'absolute';     end
if ~isfield(cfg,'xlim'),            cfg.xlim = 'maxmin';               end
if ~isfield(cfg,'ylim'),            cfg.ylim = 'maxmin';               end
if ~isfield(cfg,'zlim'),            cfg.zlim = 'maxmin';               end
if ~isfield(cfg,'colorbar'),        cfg.colorbar = 'no';               end
if ~isfield(cfg,'comment'),         cfg.comment = date;                end
if ~isfield(cfg,'showlabels'),      cfg.showlabels = 'no';             end
if ~isfield(cfg,'channel'),         cfg.channel = 'all';               end
if ~isfield(cfg,'fontsize'),        cfg.fontsize = 8;                  end
if ~isfield(cfg,'interactive'),     cfg.interactive = 'no';            end

% Set x/y/zparam defaults according to data.dimord value:
if strcmp(data.dimord, 'chan_time')
  if ~isfield(cfg, 'xparam'),      cfg.xparam='time';                  end
  if ~isfield(cfg, 'yparam'),      cfg.yparam='';                      end
  if ~isfield(cfg, 'zparam'),      cfg.zparam='avg';                   end
elseif strcmp(data.dimord, 'chan_freq')
  if ~isfield(cfg, 'xparam'),      cfg.xparam='freq';                  end
  if ~isfield(cfg, 'yparam'),      cfg.yparam='';                      end
  if ~isfield(cfg, 'zparam'),      cfg.zparam='powspctrm';             end
elseif strcmp(data.dimord, 'chan_freq_time')
  if ~isfield(cfg, 'xparam'),      cfg.xparam='time';                  end
  if ~isfield(cfg, 'yparam'),      cfg.yparam='freq';                  end
  if ~isfield(cfg, 'zparam'),      cfg.zparam='powspctrm';             end
end

% Old style coherence plotting with cohtargetchannel is no longer supported:
if isfield(cfg,'cohtargetchannel'), 
  error('cfg.cohtargetchannel is obsolete. Check the fieldtrip documentation for help about coherence plotting.');
end

% Check for unconverted coherence spectrum data:
if (strcmp(cfg.zparam,'cohspctrm')) & (~isfield(data, 'labelcmb'))
  % A reference channel is required:
  if ~isfield(cfg,'cohrefchannel'),
    error('no reference channel specified');
  end
  
  % Convert 2-dimensional channel matrix to a single dimension:
  sel1           = strmatch(cfg.cohrefchannel, data.labelcmb(:,2));
  sel2           = strmatch(cfg.cohrefchannel, data.labelcmb(:,1));
  data.cohspctrm = datavector([sel1;sel2],:,:);
  data.label     = [data.labelcmb(sel1,1);data.labelcmb(sel2,2)];
  data           = rmfield(data, 'labelcmb');
end

% Apply baseline correction:
if ~strcmp(cfg.baseline, 'no')
  data = freqbaseline(cfg, data);
end

% Get physical x-axis range:
if strcmp(cfg.xlim,'maxmin')
  xmin = min(data.(cfg.xparam));
  xmax = max(data.(cfg.xparam));
else
  xmin = cfg.xlim(1);
  xmax = cfg.xlim(2);
end

% Find corresponding x-axis bins:
xidc = find(data.(cfg.xparam) >= xmin & data.(cfg.xparam) <= xmax);

% Align physical x-axis range to the array bins:
xmin = data.(cfg.xparam)(xidc(1));
xmax = data.(cfg.xparam)(xidc(end));

% Get physical y-axis range:
if strcmp(cfg.ylim,'maxmin')
  ymin = min(data.(cfg.yparam));
  ymax = max(data.(cfg.yparam));
else
  ymin = cfg.ylim(1);
  ymax = cfg.ylim(2);
end

% Find corresponding y-axis bins:
yidc = find(data.(cfg.yparam) >= ymin & data.(cfg.yparam) <= ymax);

% Align physical y-axis range to the array bins:
ymin = data.(cfg.yparam)(yidc(1));
ymax = data.(cfg.yparam)(yidc(end));

% Read or create the layout that will be used for plotting:
if ~isfield(cfg, 'layout')
  if isfield(cfg, 'layoutname')
    cfg.layout = cfg.layoutname;       % backward compatible cfg fieldname
  elseif isfield(data, 'grad')
    cfg.layout = data.grad;            % create layout from gradiometer definition
  elseif isfield(data, 'elec')
    cfg.layout = data.elec;            % create layout from electrode definition
  else
    cfg.layout = 'CTF151s.lay';        % revert to the FCDC default
  end
end

% Select the channels in the data that match with the layout:
lay = createlayout(cfg.layout);
[seldat, sellay] = match_str(data.label, lay.label);
datavector = data.(cfg.zparam)(seldat,yidc,xidc);
chanX = lay.prj(sellay, 1);
chanY = lay.prj(sellay, 2);
chanWidth = lay.box(sellay, 1);
chanHeight = lay.box(sellay, 2);
chanLabels = lay.label(sellay);

% Get physical z-axis range (color axis):
if strcmp(cfg.zlim,'maxmin')
  zmin = min(datavector(:));
  zmax = max(datavector(:));
elseif strcmp(cfg.zlim,'absmax')
  zmin = -max(abs(datavector(:)));
  zmax = max(abs(datavector(:)));
else
  zmin = cfg.zlim(1);
  zmax = cfg.zlim(2);
end

hold on;

% Plot channels:
for k=1:length(seldat)
  % Get cdata:
  cdata = squeeze(datavector(k,:,:));
  
  % Get axes for this panel
  xas = chanX(k) + linspace(0,1,size(cdata,2))*chanWidth(k);
  yas = chanY(k) + linspace(0,1,size(cdata,1))*chanHeight(k);
  
  % Draw plot:
  imagesc(xas, yas, cdata, [zmin zmax]);
  
  % Draw channel labels:
  if strcmp(cfg.showlabels,'yes')
    text(chanX(k), chanY(k)+chanHeight(k), sprintf(' %0s\n ', chanLabels{k}), 'Fontsize', cfg.fontsize);
  end
  
  % Convert channel coordinates to the center of each channel plot:
  chanX(k) = chanX(k) + 0.5 * chanWidth(k);
  chanY(k) = chanY(k) + 0.5 * chanHeight(k);
end

% write comment:
k = cellstrmatch('COMNT',lay.label);
if ~isempty(k)
  comment = cfg.comment;
  comment = sprintf('%0s\nxlim=[%.3g %.3g]', comment, xmin, xmax);
  comment = sprintf('%0s\nylim=[%.3g %.3g]', comment, ymin, ymax);
  comment = sprintf('%0s\nzlim=[%.3g %.3g]', comment, zmin, zmax);
  text(lay.prj(k,1), lay.prj(k,2), sprintf(comment), 'Fontsize', cfg.fontsize);
  axis off;
end

% plot scale:
k = cellstrmatch('SCALE',lay.label);
if ~isempty(k)
  % Get average cdata across channels:
  cdata = squeeze(mean(datavector, 1));
  
  % Get axes for this panel:
  xas = lay.prj(k,1) + linspace(0,1,size(cdata,2))*lay.box(k,1);
  yas = lay.prj(k,2) + linspace(0,1,size(cdata,1))*lay.box(k,2);
  
  % Draw plot:
  imagesc(xas, yas, cdata, [zmin zmax]);
end

% plot colorbar:
if isfield(cfg, 'colorbar') & (strcmp(cfg.colorbar, 'yes'))
  colorbar;
end

% Make the figure interactive:
if strcmp(cfg.interactive, 'yes')
  userData.hFigure = gcf;
  userData.hAxes = gca;
  for i=1:10
    userData.hSelection{i} = plot(mean(chanX), mean(chanY));
    set(userData.hSelection{i}, 'XData', [mean(chanX)]);
    set(userData.hSelection{i}, 'YData', [mean(chanY)]);
    set(userData.hSelection{i}, 'Color', [0 0 0]);
    set(userData.hSelection{i}, 'EraseMode', 'xor');
    set(userData.hSelection{i}, 'LineStyle', '--');
    set(userData.hSelection{i}, 'LineWidth', 1.5);
    set(userData.hSelection{i}, 'Visible', 'on');
    userData.range{i} = [];
  end
  userData.iSelection = 0;
  userData.plotType = 'multiplot';
  userData.selecting = 0;
  userData.selectionType = '';
  userData.selectAxes = 'z';
  userData.lastClick = [];
  userData.cfg = cfg;
  userData.data = data;
  userData.chanX = chanX;
  userData.chanY = chanY;
  userData.chanLabels = chanLabels;
  tag = sprintf('%.5f', 10000 * rand(1));
  set(gcf, 'Renderer', 'opengl');
  set(gcf, 'Tag', tag);
  set(gcf, 'UserData', userData);
  set(gcf, 'WindowButtonMotionFcn', ['plotSelection(get(findobj(''Tag'', ''' tag '''), ''UserData''), 0);']);
  set(gcf, 'WindowButtonDownFcn', ['plotSelection(get(findobj(''Tag'', ''' tag '''), ''UserData''), 1);']);
  set(gcf, 'WindowButtonUpFcn', ['plotSelection(get(findobj(''Tag'', ''' tag '''), ''UserData''), 2);']);
end

axis tight;
orient landscape;
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function l = cellstrmatch(str,strlist)
l = [];
for k=1:length(strlist)
  if strcmp(char(str),char(strlist(k)))
    l = [l k];
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z = pwrspctm2cohspctrm(cfg, freq)
% This creates a copy of freq with new entries in the freq.pwrspctrm. The
% TFR corresponding to cfg.cohtargetchannel is set to ones, but all other
% channels are replaced with the coherence TFRs from freq.cohspctrm.
%
% [z] = pwrspctrm2cohspctrm(cfg, freq);
%
% freq should be organised in a structure as obtained from
% the freqdescriptives function.
%
% doudav; 24.09.04
z = freq;
[a a1] = match_str( cfg.cohtargetchannel, z.label);
[a b1] = match_str( cfg.cohtargetchannel, z.labelcmb(:,1));
[a b2] = match_str( cfg.cohtargetchannel, z.labelcmb(:,2));
[a c1] = match_str( z.labelcmb(b1,2), z.label);
[a c2] = match_str( z.labelcmb(b2,1), z.label);
targets  = [c1; c2];
targets2 = [b1; b2];
z.powspctrm(a1,:,:) = ones(size(z.powspctrm(a1,:,:)));
for i=1:size(targets,1),
  z.powspctrm(targets(i),:,:) = squeeze(z.cohspctrm(targets2(i),:,:));
end
clear a a1 b1 b2 c1 c2;
