text:
% topoplotUpdate() - get 2D color data for a topographic map of a scalp data field in a 2-D circular view
%              (looking down at the top of the head) using interpolation on a fine
%              cartesian grid. Can also show specified channnel location(s), or return
%              an interpolated value at an arbitrary scalp location (see 'noplot').
%              By default, channel locations below head center (arc_length 0.5) are
%              shown in a 'skirt' outside the cartoon head (see 'plotrad' and 'headrad'
%              options below). Nose is at top of plot; left is left; right is right.
%              Using option 'plotgrid', the plot may be one or more rectangular grids.
% Usage:
%        >>  topoplot(datavector, EEG.chanlocs);   % plot a map using an EEG chanlocs structure
%        >>  topoplot(datavector, 'my_chan.locs'); % read a channel locations file and plot a map
%        >>  topoplot('example');                  % give an example of an electrode location file
%        >>  [h grid_or_val plotrad_or_grid, xmesh, ymesh]= ...
%                           topoplot(datavector, chan_locs, 'Input1','Value1', ...);
% Required Inputs:
%   datavector        - single vector of channel values. Else, if a vector of selected subset
%                       (int) channel numbers -> mark their location(s) using 'style' 'blank'.
%   chan_locs         - name of an EEG electrode position file (>> topoplot example).
%                       Else, an EEG.chanlocs structure (>> help readlocs or >> topoplot example)
% Optional inputs:
%   'maplimits'       - 'absmax'   -> scale map colors to +/- the absolute-max (makes green 0);
%                       'maxmin'   -> scale colors to the data range (makes green mid-range);
%                       [lo.hi]    -> use user-definined lo/hi limits
%                       {default: 'absmax'}
%   'style'           - 'map'      -> plot colored map only
%                       'contour'  -> plot contour lines only
%                       'both'     -> plot both colored map and contour lines
%                       'fill'     -> plot constant color between contour lines
%                       'blank'    -> plot electrode locations only {default: 'both'}
%   'electrodes'      - 'on','off','labels','numbers','ptslabels','ptsnumbers'. To set the 'pts'
%                       marker,,see 'Plot detail options' below. {default: 'on' -> mark electrode
%                       locations with points ('.') unless more than 64 channels, then 'off'}.
%   'plotchans'       - [vector] channel numbers (indices) to use in making the head plot.
%                       {default: [] -> plot all chans}
%   'chantype'        - cell array of channel type(s) to plot. Will also accept a single quoted
%                       string type. Channel type for channel k is field EEG.chanlocs(k).type.
%                       If present, overrides 'plotchans' and also 'chaninfo' with field
%                       'chantype'. Ex. 'EEG' or {'EEG','EOG'} {default: all, or 'plotchans' arg}
%   'plotgrid'        - [channels] Plot channel data in one or more rectangular grids, as
%                       specified by [channels],  a position matrix of channel numbers defining
%                       the topographic locations of the channels in the grid. Zero values are
%                       given the figure background color; negative integers, the color of the
%                       polarity-reversed channel values.  Ex: >> figure; ...
%                        >> topoplot(values,'chanlocs','plotgrid',[11 12 0; 13 14 15]);
%                       % Plot a (2,3) grid of data values from channels 11-15 with one empty
%                       grid cell (top right) {default: no grid plot}
%   'nosedir'         - ['+X'|'-X'|'+Y'|'-Y'] direction of nose {default: '+X'}
%   'chaninfo'        - [struct] optional structure containing fields 'nosedir', 'plotrad'
%                       and/or 'chantype'. See these (separate) field definitions above, below.
%                       {default: nosedir +X, plotrad 0.5, all channels}
%   'plotrad'         - [0.15<=float<=1.0] plotting radius = max channel arc_length to plot.
%                       See >> topoplot example. If plotrad > 0.5, chans with arc_length > 0.5
%                       (i.e. below ears-eyes) are plotted in a circular 'skirt' outside the
%                       cartoon head. See 'intrad' below. {default: max(max(chanlocs.radius),0.5);
%                       If the chanlocs structure includes a field chanlocs.plotrad, its value
%                       is used by default}.
%   'headrad'         - [0.15<=float<=1.0] drawing radius (arc_length) for the cartoon head.
%                       NOTE: Only headrad = 0.5 is anatomically correct! 0 -> don't draw head;
%                       'rim' -> show cartoon head at outer edge of the plot {default: 0.5}
%   'intrad'          - [0.15<=float<=1.0] radius of the scalp map interpolation area (square or
%                       disk, see 'intsquare' below). Interpolate electrodes in this area and use
%                       this limit to define boundaries of the scalp map interpolated data matrix
%                       {default: max channel location radius}
%   'intsquare'       - ['on'|'off'] 'on' -> Interpolate values at electrodes located in the whole
%                       square containing the (radius intrad) interpolation disk; 'off' -> Interpolate
%                       values from electrodes shown in the interpolation disk only {default: 'on'}.
%   'conv'            - ['on'|'off'] Show map interpolation only out to the convext hull of
%                       the electrode locations to minimize extrapolation.  {default: 'off'}
%   'noplot'          - ['on'|'off'|[rad theta]] do not plot (but return interpolated data).
%                       Else, if [rad theta] are coordinates of a (possibly missing) channel,
%                       returns interpolated value for channel location.  For more info,
%                       see >> topoplot 'example' {default: 'off'}
%   'verbose'         - ['on'|'off'] comment on operations on command line {default: 'on'}.
%
% Plot detail options:
%   'drawaxis'        - ['on'|'off'] draw axis on the top left corner.
%   'emarker'         - Matlab marker char | {markerchar color size linewidth} char, else cell array
%                       specifying the electrode 'pts' marker. Ex: {'s','r',32,1} -> 32-point solid
%                       red square. {default: {'.','k',[],1} where marker size ([]) depends on the number
%                       of channels plotted}.
%   'emarker2'        - {markchans}|{markchans marker color size linewidth} cell array specifying
%                       an alternate marker for specified 'plotchans'. Ex: {[3 17],'s','g'}
%                       {default: none, or if {markchans} only are specified, then {markchans,'o','r',10,1}}
%   'hcolor'          - color of the cartoon head. Use 'hcolor','none' to plot no head. {default: 'k' = black}
%   'shading'         - 'flat','interp'  {default: 'flat'}
%   'numcontour'      - number of contour lines {default: 6}
%   'contourvals'     - values for contour {default: same as input values}
%   'pmask'           - values for masking topoplot. Array of zeros and 1 of the same size as the input
%                       value array {default: []}
%   'color'           - color of the contours {default: dark grey}
%   'whitebk '        -  ('on'|'off') make the background color white (e.g., to print empty plotgrid channels)
%                       {default: 'off'}
%   'gridscale'       - [int > 32] size (nrows) of interpolated scalp map data matrix {default: 67}
%   'colormap'        -  (n,3) any size colormap {default: existing colormap}
%   'circgrid'        - [int > 100] number of elements (angles) in head and border circles {201}
%   'emarkercolor'    - cell array of colors for 'blank' option.
%   'plotdisk'        - ['on'|'off'] plot disk instead of dots for electrodefor 'blank' option. Size of disk
%                       is controled by input values at each electrode. If an imaginary value is provided,
%                       plot partial circle with red for the real value and blue for the imaginary one.
%
% Dipole plotting options:
%   'dipole'          - [xi yi xe ye ze] plot dipole on the top of the scalp map
%                       from coordinate (xi,yi) to coordinates (xe,ye,ze) (dipole head
%                       model has radius 1). If several rows, plot one dipole per row.
%                       Coordinates returned by dipplot() may be used. Can accept
%                       an EEG.dipfit.model structure (See >> help dipplot).
%                       Ex: ,'dipole',EEG.dipfit.model(17) % Plot dipole(s) for comp. 17.
%   'dipnorm'         - ['on'|'off'] normalize dipole length {default: 'on'}.
%   'diporient'       - [-1|1] invert dipole orientation {default: 1}.
%   'diplen'          - [real] scale dipole length {default: 1}.
%   'dipscale'        - [real] scale dipole size {default: 1}.
%   'dipsphere'       - [real] size of the dipole sphere. {default: 85 mm}.
%   'dipcolor'        - [color] dipole color as Matlab code code or [r g b] vector
%                       {default: 'k' = black}.
% Outputs:
%              handle - handle of the colored surface.If
%                       contour only is plotted, then is the handle of
%                       the countourgroup. (If no surface or contour is plotted,
%                       return "gca", the handle of the current plot)
%         grid_or_val - [matrix] the interpolated data image (with off-head points = NaN).
%                       Else, single interpolated value at the specified 'noplot' arg channel
%                       location ([rad theta]), if any.
%     plotrad_or_grid - IF grid image returned above, then the 'plotrad' radius of the grid.
%                       Else, the grid image
%     xmesh, ymesh    - x and y values of the returned grid (above)
%
% Chan_locs format:
%    See >> topoplot 'example'
%
% Examples:
%
%    To plot channel locations only:
%    >> figure; topoplot([],EEG.chanlocs,'style','blank','electrodes','labelpoint','chaninfo',EEG.chaninfo);
%
% Notes: - To change the plot map masking ring to a new figure background color,
%            >> set(findobj(gca,'type','patch'),'facecolor',get(gcf,'color'))
%        - Topoplots may be rotated. From the commandline >> view([deg 90]) {default: [0 90])
%
% Authors: Andy Spydell, Colin Humphries, Arnaud Delorme & Scott Makeig
%          CNL / Salk Institute, 8/1996-/10/2001; SCCN/INC/UCSD, Nov. 2001 -
%
% See also: timtopo(), envtopo()

% Deprecated options:
%           'shrink' - ['on'|'off'|'force'|factor] Deprecated. 'on' -> If max channel arc_length
%                       > 0.5, shrink electrode coordinates towards vertex to plot all channels
%                       by making max arc_length 0.5. 'force' -> Normalize arc_length
%                       so the channel max is 0.5. factor -> Apply a specified shrink
%                       factor (range (0,1) = shrink fraction). {default: 'off'}
%   'electcolor' {'k'}  ... electrode marking details and their {defaults}.
%   'emarker' {'.'}|'emarkersize' {14}|'emarkersizemark' {40}|'efontsize' {var} -
%                       electrode marking details and their {defaults}.
%   'ecolor'          - color of the electrode markers {default: 'k' = black}
%   'interplimits'    - ['electrodes'|'head'] 'electrodes'-> interpolate the electrode grid;
%                       'head'-> interpolate the whole disk {default: 'head'}.

% Unimplemented future options:

% Copyright (C) Colin Humphries & Scott Makeig, CNL / Salk Institute, Aug, 1996
%
% This program 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 2 of the License, or
% (at your option) any later version.
%
% This program 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 this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% Topoplot Version 2.1
% Early development history:
% Begun by Andy Spydell and Scott Makeig, NHRC,  7-23-96
% 8-96 Revised by Colin Humphries, CNL / Salk Institute, La Jolla CA
%   -changed surf command to imagesc (faster)
%   -can now handle arbitrary scaling of electrode distances
%   -can now handle non integer angles in chan_locs
% 4-4-97 Revised again by Colin Humphries, reformatted by SM
%   -added parameters
%   -changed chan_locs format
% 2-26-98 Revised by Colin
%   -changed image back to surface command
%   -added fill and blank styles
%   -removed extra background colormap entry (now use any colormap)
%   -added parameters for electrode colors and labels
%   -now each topoplot axes use the caxis command again.
%   -removed OUTPUT parameter
% 3-11-98 changed default emarkersize, improve help msg -sm
% 5-24-01 made default emarkersize vary with number of channels -sm
% 01-25-02 reformated help & license, added link -ad
% 03-15-02 added readlocs and the use of eloc input structure -ad
% 03-25-02 added 'labelpoint' options and allow Values=[] -ad &sm
% 03-25-02 added details to "Unknown parameter" warning -sm & ad

function [Zi,amin,amax] = topoplotUpdate(Values,loc_file,varargin)

%
%%%%%%%%%%%%%%%%%%%%%%%% Set defaults %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
icadefs                 % read defaults MAXTOPOPLOTCHANS and DEFAULT_ELOC and BACKCOLOR
if ~exist('BACKCOLOR')  % if icadefs.m does not define BACKCOLOR
   BACKCOLOR = [.93 .96 1];  % EEGLAB standard
end
whitebk = 'off';  % by default, make gridplot background color = EEGLAB screen background color

persistent warningInterp;

plotgrid = 'off';
plotchans = [];
noplot  = 'off';
handle = [];
Zi = [];
chanval = NaN;
rmax = 0.5;             % actual head radius - Don't change this!
INTERPLIMITS = 'head';  % head, electrodes
INTSQUARE = 'on';       % default, interpolate electrodes located though the whole square containing
                        % the plotting disk
default_intrad = 1;     % indicator for (no) specified intrad
MAPLIMITS = 'absmax';   % absmax, maxmin, [values]
GRID_SCALE = 67;        % plot map on a 67X67 grid
CIRCGRID   = 201;       % number of angles to use in drawing circles
AXHEADFAC = 1.3;        % head to axes scaling factor
CONTOURNUM = 6;         % number of contour levels to plot
STYLE = 'both';         % default 'style': both,straight,fill,contour,blank
HEADCOLOR = [0 0 0];    % default head color (black)
CCOLOR = [0.2 0.2 0.2]; % default contour color
ELECTRODES = [];        % default 'electrodes': on|off|label - set below
MAXDEFAULTSHOWLOCS = 64;% if more channels than this, don't show electrode locations by default
EMARKER = '.';          % mark electrode locations with small disks
ECOLOR = [0 0 0];       % default electrode color = black
EMARKERSIZE = [];       % default depends on number of electrodes, set in code
EMARKERLINEWIDTH = 1;   % default edge linewidth for emarkers
EMARKERSIZE1CHAN = 20;  % default selected channel location marker size
EMARKERCOLOR1CHAN = 'red'; % selected channel location marker color
EMARKER2CHANS = [];      % mark subset of electrode locations with small disks
EMARKER2 = 'o';          % mark subset of electrode locations with small disks
EMARKER2COLOR = 'r';     % mark subset of electrode locations with small disks
EMARKERSIZE2 = 10;      % default selected channel location marker size
EMARKER2LINEWIDTH = 1;
EFSIZE = get(0,'DefaultAxesFontSize'); % use current default fontsize for electrode labels
HLINEWIDTH = 2;         % default linewidth for head, nose, ears
BLANKINGRINGWIDTH = .035;% width of the blanking ring
HEADRINGWIDTH    = .007;% width of the cartoon head ring
SHADING = 'flat';       % default 'shading': flat|interp
shrinkfactor = [];      % shrink mode (dprecated)
intrad       = [];      % default interpolation square is to outermost electrode (<=1.0)
plotrad      = [];      % plotting radius ([] = auto, based on outermost channel location)
headrad      = [];      % default plotting radius for cartoon head is 0.5
squeezefac = 1.0;
MINPLOTRAD = 0.15;      % can't make a topoplot with smaller plotrad (contours fail)
VERBOSE = 'off';
MASKSURF = 'off';
CONVHULL = 'off';       % dont mask outside the electrodes convex hull
DRAWAXIS = 'off';
PLOTDISK = 'off';
CHOOSECHANTYPE = 0;
ContourVals = Values;
PMASKFLAG   = 0;
COLORARRAY  = { [1 0 0] [0.5 0 0] [0 0 0] };
%COLORARRAY2 = { [1 0 0] [0.5 0 0] [0 0 0] };
gb = [0 0];
COLORARRAY2 = { [gb 0] [gb 1/4] [gb 2/4] [gb 3/4] [gb 1] };

%%%%%% Dipole defaults %%%%%%%%%%%%
DIPOLE  = [];
DIPNORM   = 'on';
DIPSPHERE = 85;
DIPLEN    = 1;
DIPSCALE  = 1;
DIPORIENT  = 1;
DIPCOLOR  = [0 0 0];
NOSEDIR   = '+X';
CHANINFO  = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
%%%%%%%%%%%%%%%%%%%%%%% Handle arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if nargin< 1
   help topoplot;
   return
end

% calling topoplot from Fieldtrip
% -------------------------------
fieldtrip = 0;
if nargin < 2, loc_file = []; end;
if isstruct(Values) | ~isstruct(loc_file), fieldtrip == 1; end;
if isstr(loc_file), if exist(loc_file) ~= 2, fieldtrip == 1; end; end;
if fieldtrip
    error('Wrong calling format, are you trying to use the topoplot Fieldtrip function?');
end;

nargs = nargin;
if nargs == 1
  if isstr(Values)
    if any(strcmp(lower(Values),{'example','demo'}))
      fprintf(['This is an example of an electrode location file,\n',...
               'an ascii file consisting of the following four columns:\n',...
               ' channel_number degrees arc_length channel_name\n\n',...
               'Example:\n',...
               ' 1               -18    .352       Fp1 \n',...
               ' 2                18    .352       Fp2 \n',...
               ' 5               -90    .181       C3  \n',...
               ' 6                90    .181       C4  \n',...
               ' 7               -90    .500       A1  \n',...
               ' 8                90    .500       A2  \n',...
               ' 9              -142    .231       P3  \n',...
               '10               142    .231       P4  \n',...
               '11                 0    .181       Fz  \n',...
               '12                 0    0          Cz  \n',...
               '13               180    .181       Pz  \n\n',...
                                                             ...
               'In topoplot() coordinates, 0 deg. points to the nose, positive\n',...
               'angles point to the right hemisphere, and negative to the left.\n',...
               'The model head sphere has a circumference of 2; the vertex\n',...
               '(Cz) has arc_length 0. Locations with arc_length > 0.5 are below\n',...
               'head center and are plotted outside the head cartoon.\n',...
               'Option plotrad controls how much of this lower-head "skirt" is shown.\n',...
               'Option headrad controls if and where the cartoon head will be drawn.\n',...
               'Option intrad controls how many channels will be included in the interpolation.\n',...
               ])
      return
    end
  end
end
if nargs < 2
  loc_file = DEFAULT_ELOC;
  if ~exist(loc_file)
      fprintf('default locations file "%s" not found - specify chan_locs in topoplot() call.\n',loc_file)
      error(' ')
  end
end
if isempty(loc_file)
  loc_file = 0;
end
if isnumeric(loc_file) & loc_file == 0
  loc_file = DEFAULT_ELOC;
end

if nargs > 2
    if ~(round(nargs/2) == nargs/2)
        error('Odd number of input arguments??')
    end
    for i = 1:2:length(varargin)
        Param = varargin{i};
        Value = varargin{i+1};
        if ~isstr(Param)
            error('Flag arguments must be strings')
        end
        Param = lower(Param);
        switch Param
            case 'conv'
                CONVHULL = lower(Value);
                if ~strcmp(CONVHULL,'on') & ~strcmp(CONVHULL,'off')
                    error('Value of ''conv'' must be ''on'' or ''off''.');
                end
            case 'colormap'
                if size(Value,2)~=3
                    error('Colormap must be a n x 3 matrix')
                end
                colormap(Value)
            case 'plotdisk'
                PLOTDISK = lower(Value);
                if ~strcmp(PLOTDISK,'on') & ~strcmp(PLOTDISK,'off')
                    error('Value of ''plotdisk'' must be ''on'' or ''off''.');
                end
            case 'intsquare'
                INTSQUARE = lower(Value);
                if ~strcmp(INTSQUARE,'on') & ~strcmp(INTSQUARE,'off')
                    error('Value of ''intsquare'' must be ''on'' or ''off''.');
                end
            case 'emarkercolors'
                COLORARRAY = Value;
            case {'interplimits','headlimits'}
                if ~isstr(Value)
                    error('''interplimits'' value must be a string')
                end
                Value = lower(Value);
                if ~strcmp(Value,'electrodes') & ~strcmp(Value,'head')
                    error('Incorrect value for interplimits')
                end
                INTERPLIMITS = Value;
            case 'verbose'
                VERBOSE = Value;
            case 'nosedir'
                NOSEDIR = Value;
                if isempty(strmatch(lower(NOSEDIR), { '+x', '-x', '+y', '-y' }))
                    error('Invalid nose direction');
                end;
            case 'chaninfo'
                CHANINFO = Value;
                if isfield(CHANINFO, 'nosedir'), NOSEDIR      = CHANINFO.nosedir; end;
                if isfield(CHANINFO, 'shrink' ), shrinkfactor = CHANINFO.shrink;  end;
                if isfield(CHANINFO, 'plotrad') & isempty(plotrad), plotrad = CHANINFO.plotrad; end;
                if isfield(CHANINFO, 'chantype')
                    chantype = CHANINFO.chantype;
                    if ischar(chantype), chantype = cellstr(chantype); end
                    CHOOSECHANTYPE = 1;
                end
            case 'chantype'
                chantype = Value;
                CHOOSECHANTYPE = 1;
                if ischar(chantype), chantype = cellstr(chantype); end
                if ~iscell(chantype), error('chantype must be cell array. e.g. {''EEG'', ''EOG''}'); end
            case 'drawaxis'
                DRAWAXIS = Value;
            case 'maplimits'
                MAPLIMITS = Value;
            case 'masksurf'
                MASKSURF = Value;
            case 'circgrid'
                CIRCGRID = Value;
                if isstr(CIRCGRID) | CIRCGRID<100
                    error('''circgrid'' value must be an int > 100');
                end
            case 'style'
                STYLE = lower(Value);
            case 'numcontour'
                CONTOURNUM = Value;
            case 'electrodes'
                ELECTRODES = lower(Value);
                if strcmpi(ELECTRODES,'pointlabels') | strcmpi(ELECTRODES,'ptslabels') ...
                        | strcmpi(ELECTRODES,'labelspts') | strcmpi(ELECTRODES,'ptlabels') ...
                        | strcmpi(ELECTRODES,'labelpts')
                    ELECTRODES = 'labelpoint'; % backwards compatability
                elseif strcmpi(ELECTRODES,'pointnumbers') | strcmpi(ELECTRODES,'ptsnumbers') ...
                        | strcmpi(ELECTRODES,'numberspts') | strcmpi(ELECTRODES,'ptnumbers') ...
                        | strcmpi(ELECTRODES,'numberpts')  | strcmpi(ELECTRODES,'ptsnums')  ...
                        | strcmpi(ELECTRODES,'numspts')
                    ELECTRODES = 'numpoint'; % backwards compatability
                elseif strcmpi(ELECTRODES,'nums')
                    ELECTRODES = 'numbers'; % backwards compatability
                elseif strcmpi(ELECTRODES,'pts')
                    ELECTRODES = 'on'; % backwards compatability
                elseif ~strcmp(ELECTRODES,'off') ...
                        & ~strcmpi(ELECTRODES,'on') ...
                        & ~strcmp(ELECTRODES,'labels') ...
                        & ~strcmpi(ELECTRODES,'numbers') ...
                        & ~strcmpi(ELECTRODES,'labelpoint') ...
                        & ~strcmpi(ELECTRODES,'numpoint')
                    error('Unknown value for keyword ''electrodes''');
                end
            case 'dipole'
                DIPOLE = Value;
            case 'dipsphere'
                DIPSPHERE = Value;
            case 'dipnorm'
                DIPNORM = Value;
            case 'diplen'
                DIPLEN = Value;
            case 'dipscale'
                DIPSCALE = Value;
            case 'contourvals'
                ContourVals = Value;
            case 'pmask'
                ContourVals = Value;
                PMASKFLAG   = 1;
            case 'diporient'
                DIPORIENT = Value;
            case 'dipcolor'
                DIPCOLOR = Value;
            case 'emarker'
                if ischar(Value)
                    EMARKER = Value;
                elseif ~iscell(Value) | length(Value) > 4
                    error('''emarker'' argument must be a cell array {marker color size linewidth}')
                else
                    EMARKER = Value{1};
                end
                if length(Value) > 1
                    ECOLOR = Value{2};
                end
                if length(Value) > 2
                    EMARKERSIZE = Value{3};
                end
                if length(Value) > 3
                    EMARKERLINEWIDTH = Value{4};
                end
            case 'emarker2'
                if ~iscell(Value) | length(Value) > 5
                    error('''emarker2'' argument must be a cell array {chans marker color size linewidth}')
                end
                EMARKER2CHANS = abs(Value{1}); % ignore channels < 0
                if length(Value) > 1
                    EMARKER2 = Value{2};
                end
                if length(Value) > 2
                    EMARKER2COLOR = Value{3};
                end
                if length(Value) > 3
                    EMARKERSIZE2 = Value{4};
                end
                if length(Value) > 4
                    EMARKER2LINEWIDTH = Value{5};
                end
            case 'shrink'
                shrinkfactor = Value;
            case 'intrad'
                intrad = Value;
                if isstr(intrad) | (intrad < MINPLOTRAD | intrad > 1)
                    error('intrad argument should be a number between 0.15 and 1.0');
                end
            case 'plotrad'
                plotrad = Value;
                if isstr(plotrad) | (plotrad < MINPLOTRAD | plotrad > 1)
                    error('plotrad argument should be a number between 0.15 and 1.0');
                end
            case 'headrad'
                headrad = Value;
                if isstr(headrad) & ( strcmpi(headrad,'off') | strcmpi(headrad,'none') )
                    headrad = 0;       % undocumented 'no head' alternatives
                end
                if isempty(headrad) % [] -> none also
                    headrad = 0;
                end
                if ~isstr(headrad)
                    if ~(headrad==0) & (headrad < MINPLOTRAD | headrad>1)
                        error('bad value for headrad');
                    end
                elseif  ~strcmpi(headrad,'rim')
                    error('bad value for headrad');
                end
            case {'headcolor','hcolor'}
                HEADCOLOR = Value;
            case {'contourcolor','ccolor'}
                CCOLOR = Value;
            case {'electcolor','ecolor'}
                ECOLOR = Value;
            case {'emarkersize','emsize'}
                EMARKERSIZE = Value;
            case {'emarkersize1chan','emarkersizemark'}
                EMARKERSIZE1CHAN= Value;
            case {'efontsize','efsize'}
                EFSIZE = Value;
            case 'shading'
                SHADING = lower(Value);
                if ~any(strcmp(SHADING,{'flat','interp'}))
                    error('Invalid shading parameter')
                end
                if strcmpi(SHADING,'interp') && isempty(warningInterp)
                    warning('Using interpolated shading in scalp topographies prevent to export them as vectorized figures');
                    warningInterp = 1;
                end;
            case 'noplot'
                noplot = Value;
                if ~isstr(noplot)
                    if length(noplot) ~= 2
                        error('''noplot'' location should be [radius, angle]')
                    else
                        chanrad = noplot(1);
                        chantheta = noplot(2);
                        noplot = 'on';
                    end
                end
            case 'gridscale'
                GRID_SCALE = Value;
                if isstr(GRID_SCALE) | GRID_SCALE ~= round(GRID_SCALE) | GRID_SCALE < 32
                    error('''gridscale'' value must be integer > 32.');
                end
            case {'plotgrid','gridplot'}
                plotgrid = 'on';
                gridchans = Value;
            case 'plotchans'
                plotchans = Value(:);
                if find(plotchans<=0)
                    error('''plotchans'' values must be > 0');
                end
                % if max(abs(plotchans))>max(Values) | max(abs(plotchans))>length(Values) -sm ???
            case {'whitebk','whiteback','forprint'}
                whitebk = Value;
%             case 'handle'
%                 cmaphand = Value;
            otherwise
                error(['Unknown input parameter ''' Param ''' ???'])
        end
    end
end
if strcmpi(whitebk, 'on')
    BACKCOLOR = [ 1 1 1 ];
end;

% cmap = get(cmaphand,'Colormap');
% cmaplen = size(cmap,1);

if strcmp(STYLE,'blank')    % else if Values holds numbers of channels to mark
    if length(Values) < length(loc_file)
        ContourVals = zeros(1,length(loc_file));
        ContourVals(Values) = 1;
        Values = ContourVals;
    end;
end;

%
%%%%%%%%%%%%%%%%%%%%%%%%%%% test args for plotting an electrode grid %%%%%%%%%%%%%%%%%%%%%%
%
if strcmp(plotgrid,'on')
   STYLE = 'grid';
   gchans = sort(find(abs(gridchans(:))>0));

   % if setdiff(gchans,unique(gchans))
   %      fprintf('topoplot() warning: ''plotgrid'' channel matrix has duplicate channels\n');
   % end

   if ~isempty(plotchans)
     if intersect(gchans,abs(plotchans))
        fprintf('topoplot() warning: ''plotgrid'' and ''plotchans'' have channels in common\n');
     end
   end
end

%
%%%%%%%%%%%%%%%%%%%%%%%%%%% misc arg tests %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if isempty(ELECTRODES)                     % if electrode labeling not specified
  if length(Values) > MAXDEFAULTSHOWLOCS   % if more channels than default max
    ELECTRODES = 'off';                    % don't show electrodes
  else                                     % else if fewer chans,
    ELECTRODES = 'on';                     % do
  end
end

if isempty(Values)
   STYLE = 'blank';
end
[r,c] = size(Values);
if r>1 & c>1,
  error('input data must be a single vector');
end
Values = Values(:); % make Values a column vector
ContourVals = ContourVals(:); % values for contour

if ~isempty(intrad) & ~isempty(plotrad) & intrad < plotrad
   error('intrad must be >= plotrad');
end

if ~strcmpi(STYLE,'grid')                     % if not plot grid only

%
%%%%%%%%%%%%%%%%%%%% Read the channel location information %%%%%%%%%%%%%%%%%%%%%%%%
%
  if isstr(loc_file)
      [tmpeloc labels Th Rd indices] = readlocs( loc_file);
  elseif isstruct(loc_file) % a locs struct
      [tmpeloc labels Th Rd indices] = readlocs( loc_file );
      % Note: Th and Rd correspond to indices channels-with-coordinates only
  else
       error('loc_file must be a EEG.locs struct or locs filename');
  end
  Th = pi/180*Th;                              % convert degrees to radians
  allchansind = 1:length(Th);

  if ~isempty(plotchans)
      if max(plotchans) > length(Th)
          error('''plotchans'' values must be <= max channel index');
      end
  end

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% channels to plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if ~isempty(plotchans)
    plotchans = intersect(plotchans, indices);
end;
if ~isempty(Values) & ~strcmpi( STYLE, 'blank') & isempty(plotchans)
    plotchans = indices;
end
if isempty(plotchans) & strcmpi( STYLE, 'blank')
    plotchans = indices;
end

%
%%%%%%%%%%%%%%%%%%%%%%%%%%% filter for channel type(s), if specified %%%%%%%%%%%%%%%%%%%%%
%

if CHOOSECHANTYPE,
    newplotchans = eeg_chantype(loc_file,chantype);
    plotchans = intersect(newplotchans, plotchans);
end

%
%%%%%%%%%%%%%%%%%%%%%%%%%%% filter channels used for components %%%%%%%%%%%%%%%%%%%%%
%
if isfield(CHANINFO, 'icachansind') & ~isempty(Values) & length(Values) ~= length(tmpeloc)

    % test if ICA component
    % ---------------------
    if length(CHANINFO.icachansind) == length(Values)

        % if only a subset of channels are to be plotted
        % and ICA components also use a subject of channel
        % we must find the new indices for these channels

        plotchans = intersect(CHANINFO.icachansind, plotchans);
        tmpvals   = zeros(1, length(tmpeloc));
        tmpvals(CHANINFO.icachansind) = Values;
        Values    = tmpvals;
        tmpvals   = zeros(1, length(tmpeloc));
        tmpvals(CHANINFO.icachansind) = ContourVals;
        ContourVals = tmpvals;

    end;
end;

%
%%%%%%%%%%%%%%%%%%% last channel is reference? %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if length(tmpeloc) == length(Values) + 1 % remove last channel if necessary
                                         % (common reference channel)
    if plotchans(end) == length(tmpeloc)
        plotchans(end) = [];
    end;

end;

%
%%%%%%%%%%%%%%%%%%% remove infinite and NaN values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if length(Values) > 1
    inds          = union(find(isnan(Values)), find(isinf(Values))); % NaN and Inf values
    plotchans     = setdiff(plotchans, inds);
end;
if strcmp(plotgrid,'on')
    plotchans = setxor(plotchans,gchans);   % remove grid chans from head plotchans
end

[x,y]     = pol2cart(Th,Rd);  % transform electrode locations from polar to cartesian coordinates
plotchans = abs(plotchans);   % reverse indicated channel polarities
allchansind = allchansind(plotchans);
Th        = Th(plotchans);
Rd        = Rd(plotchans);
x         = x(plotchans);
y         = y(plotchans);
labels    = labels(plotchans); % remove labels for electrodes without locations
labels    = strvcat(labels); % make a label string matrix
if ~isempty(Values) & length(Values) > 1
    Values      = Values(plotchans);
    ContourVals = ContourVals(plotchans);
end;

%
%%%%%%%%%%%%%%%%%% Read plotting radius from chanlocs  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if isempty(plotrad) & isfield(tmpeloc, 'plotrad'),
    plotrad = tmpeloc(1).plotrad;
    if isstr(plotrad)                        % plotrad shouldn't be a string
        plotrad = str2num(plotrad)           % just checking
    end
    if plotrad < MINPLOTRAD | plotrad > 1.0
       fprintf('Bad value (%g) for plotrad.\n',plotrad);
       error(' ');
    end
    if strcmpi(VERBOSE,'on') & ~isempty(plotrad)
       fprintf('Plotting radius plotrad (%g) set from EEG.chanlocs.\n',plotrad);
    end
end;
if isempty(plotrad)
  plotrad = min(1.0,max(Rd)*1.02);            % default: just outside the outermost electrode location
  plotrad = max(plotrad,0.5);                 % default: plot out to the 0.5 head boundary
end                                           % don't plot channels with Rd > 1 (below head)

if isempty(intrad)
  default_intrad = 1;     % indicator for (no) specified intrad
  intrad = min(1.0,max(Rd)*1.02);             % default: just outside the outermost electrode location
else
  default_intrad = 0;                         % indicator for (no) specified intrad
  if plotrad > intrad
     plotrad = intrad;
  end
end                                           % don't interpolate channels with Rd > 1 (below head)
if isstr(plotrad) | plotrad < MINPLOTRAD | plotrad > 1.0
   error('plotrad must be between 0.15 and 1.0');
end

%
%%%%%%%%%%%%%%%%%%%%%%% Set radius of head cartoon %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if isempty(headrad)  % never set -> defaults
  if plotrad >= rmax
     headrad = rmax;  % (anatomically correct)
  else % if plotrad < rmax
     headrad = 0;    % don't plot head
     if strcmpi(VERBOSE, 'on')
       fprintf('topoplot(): not plotting cartoon head since plotrad (%5.4g) < 0.5\n',...
                                                                    plotrad);
     end
  end
elseif strcmpi(headrad,'rim') % force plotting at rim of map
  headrad = plotrad;
end

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Shrink mode %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if ~isempty(shrinkfactor) | isfield(tmpeloc, 'shrink'),
    if isempty(shrinkfactor) & isfield(tmpeloc, 'shrink'),
        shrinkfactor = tmpeloc(1).shrink;
        if strcmpi(VERBOSE,'on')
            if isstr(shrinkfactor)
                fprintf('Automatically shrinking coordinates to lie above the head perimter.\n');
            else
                fprintf('Automatically shrinking coordinates by %3.2f\n', shrinkfactor);
            end;
        end
    end;

    if isstr(shrinkfactor)
        if strcmpi(shrinkfactor, 'on') | strcmpi(shrinkfactor, 'force') | strcmpi(shrinkfactor, 'auto')
            if abs(headrad-rmax) > 1e-2
             fprintf('     NOTE -> the head cartoon will NOT accurately indicate the actual electrode locations\n');
            end
            if strcmpi(VERBOSE,'on')
                fprintf('     Shrink flag -> plotting cartoon head at plotrad\n');
            end
            headrad = plotrad; % plot head around outer electrodes, no matter if 0.5 or not
        end
    else % apply shrinkfactor
        plotrad = rmax/(1-shrinkfactor);
        headrad = plotrad;  % make deprecated 'shrink' mode plot
        if strcmpi(VERBOSE,'on')
            fprintf('    %g%% shrink  applied.');
            if abs(headrad-rmax) > 1e-2
                fprintf(' Warning: With this "shrink" setting, the cartoon head will NOT be anatomically correct.\n');
            else
                fprintf('\n');
            end
        end
    end
end; % if shrink

%
%%%%%%%%%%%%%%%%% Issue warning if headrad ~= rmax  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

if headrad ~= 0.5 & strcmpi(VERBOSE, 'on')
   fprintf('     NB: Plotting map using ''plotrad'' %-4.3g,',plotrad);
   fprintf(    ' ''headrad'' %-4.3g\n',headrad);
   fprintf('Warning: The plotting radius of the cartoon head is NOT anatomically correct (0.5).\n')
end
%
%%%%%%%%%%%%%%%%%%%%% Find plotting channels  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

pltchans = find(Rd <= plotrad); % plot channels inside plotting circle

if strcmpi(INTSQUARE,'on') % interpolate channels in the radius intrad square
  intchans = find(x <= intrad & y <= intrad); % interpolate and plot channels inside interpolation square
else
  intchans = find(Rd <= intrad); % interpolate channels in the radius intrad circle only
end

%
%%%%%%%%%%%%%%%%%%%%% Eliminate channels not plotted  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

allx      = x;
ally      = y;
intchans; % interpolate using only the 'intchans' channels
pltchans; % plot using only indicated 'plotchans' channels

if length(pltchans) < length(Rd) & strcmpi(VERBOSE, 'on')
        fprintf('Interpolating %d and plotting %d of the %d scalp electrodes.\n', ...
                   length(intchans),length(pltchans),length(Rd));
end;

% fprintf('topoplot(): plotting %d channels\n',length(pltchans));
if ~isempty(EMARKER2CHANS)
    if strcmpi(STYLE,'blank')
       error('emarker2 not defined for style ''blank'' - use marking channel numbers in place of data');
    else % mark1chans and mark2chans are subsets of pltchans for markers 1 and 2
       [tmp1 mark1chans tmp2] = setxor(pltchans,EMARKER2CHANS);
       [tmp3 tmp4 mark2chans] = intersect(EMARKER2CHANS,pltchans);
    end
end

if ~isempty(Values)
	if length(Values) == length(Th)  % if as many map Values as channel locs
		intValues      = Values(intchans);
		intContourVals = ContourVals(intchans);
        Values         = Values(pltchans);
		ContourVals    = ContourVals(pltchans);
	end;
end;   % now channel parameters and values all refer to plotting channels only

allchansind = allchansind(pltchans);
intTh = Th(intchans);           % eliminate channels outside the interpolation area
intRd = Rd(intchans);
intx  = x(intchans);
inty  = y(intchans);
Th    = Th(pltchans);              % eliminate channels outside the plotting area
Rd    = Rd(pltchans);
x     = x(pltchans);
y     = y(pltchans);

labels= labels(pltchans,:);
%
%%%%%%%%%%%%%%% Squeeze channel locations to <= rmax %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

squeezefac = rmax/plotrad;
intRd = intRd*squeezefac; % squeeze electrode arc_lengths towards the vertex
Rd = Rd*squeezefac;       % squeeze electrode arc_lengths towards the vertex
                          % to plot all inside the head cartoon
intx = intx*squeezefac;
inty = inty*squeezefac;
x    = x*squeezefac;
y    = y*squeezefac;
allx    = allx*squeezefac;
ally    = ally*squeezefac;
% Note: Now outermost channel will be plotted just inside rmax

else % if strcmpi(STYLE,'grid')
   intx = rmax; inty=rmax;
end % if ~strcmpi(STYLE,'grid')

%
%%%%%%%%%%%%%%%% rotate channels based on chaninfo %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if strcmpi(lower(NOSEDIR), '+x')
     rotate = 0;
else
    if strcmpi(lower(NOSEDIR), '+y')
        rotate = 3*pi/2;
    elseif strcmpi(lower(NOSEDIR), '-x')
        rotate = pi;
    else rotate = pi/2;
    end;
    allcoords = (inty + intx*sqrt(-1))*exp(sqrt(-1)*rotate);
    intx = imag(allcoords);
    inty = real(allcoords);
    allcoords = (ally + allx*sqrt(-1))*exp(sqrt(-1)*rotate);
    allx = imag(allcoords);
    ally = real(allcoords);
    allcoords = (y + x*sqrt(-1))*exp(sqrt(-1)*rotate);
    x = imag(allcoords);
    y = real(allcoords);
end;

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Make the plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if ~strcmpi(STYLE,'blank') % if draw interpolated scalp map
 if ~strcmpi(STYLE,'grid') %  not a rectangular channel grid
  %
  %%%%%%%%%%%%%%%% Find limits for interpolation %%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  if default_intrad % if no specified intrad
   if strcmpi(INTERPLIMITS,'head') % intrad is 'head'
    xmin = min(-rmax,min(intx)); xmax = max(rmax,max(intx));
    ymin = min(-rmax,min(inty)); ymax = max(rmax,max(inty));

   else % INTERPLIMITS = rectangle containing electrodes -- DEPRECATED OPTION!
    xmin = max(-rmax,min(intx)); xmax = min(rmax,max(intx));
    ymin = max(-rmax,min(inty)); ymax = min(rmax,max(inty));
   end
  else % some other intrad specified
    xmin = -intrad*squeezefac; xmax = intrad*squeezefac;   % use the specified intrad value
    ymin = -intrad*squeezefac; ymax = intrad*squeezefac;
  end
  %
  %%%%%%%%%%%%%%%%%%%%%%% Interpolate scalp map data %%%%%%%%%%%%%%%%%%%%%%%%
  %
  xi = linspace(xmin,xmax,GRID_SCALE);   % x-axis description (row vector)
  yi = linspace(ymin,ymax,GRID_SCALE);   % y-axis description (row vector)

  try
      [Xi,Yi,Zi] = griddata(inty,intx,double(intValues),yi',xi,'v4');%,'v4'); % interpolate data
%       [Xi,Yi,ZiC] = griddata(inty,intx,double(intContourVals),yi',xi,'v4'); % interpolate data
  catch,
      [Xi,Yi,Zi] = griddata(inty,intx,intValues',yi,xi'); % interpolate data (Octave)
%       [Xi,Yi,ZiC] = griddata(inty,intx,intContourVals',yi,xi'); % interpolate data
  end;
  %
  %%%%%%%%%%%%%%%%%%%%%%% Mask out data outside the head %%%%%%%%%%%%%%%%%%%%%
  %
  mask = (sqrt(Xi.^2 + Yi.^2) <= rmax); % mask outside the plotting circle
  ii = find(mask == 0);
  Zi(ii)  = NaN;                         % mask non-plotting voxels with NaNs
%   ZiC(ii) = NaN;                         % mask non-plotting voxels with NaNs
  grid = plotrad;                       % unless 'noplot', then 3rd output arg is plotrad
  %
  %%%%%%%%%% Return interpolated value at designated scalp location %%%%%%%%%%
  %
  if exist('chanrad')   % optional first argument to 'noplot'
      chantheta = (chantheta/360)*2*pi;
      chancoords = round(ceil(GRID_SCALE/2)+GRID_SCALE/2*2*chanrad*[cos(-chantheta),...
                                                      -sin(-chantheta)]);
      if chancoords(1)<1 ...
         | chancoords(1) > GRID_SCALE ...
            | chancoords(2)<1 ...
               | chancoords(2)>GRID_SCALE
          error('designated ''noplot'' channel out of bounds')
      else
        chanval = Zi(chancoords(1),chancoords(2));
        grid = Zi;
        Zi = chanval;  % return interpolated value instead of Zi
      end
  end
  %
  %%%%%%%%%%%%%%%%%%%%%%%%%% Return interpolated image only  %%%%%%%%%%%%%%%%%
  %
   if strcmpi(noplot, 'on')
    if strcmpi(VERBOSE,'on')
       fprintf('topoplot(): no plot requested.\n')
    end
    return;
   end
  %
  %%%%%%%%%%%%%%%%%%%%%%% Calculate colormap limits %%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  if isstr(MAPLIMITS)
    if strcmp(MAPLIMITS,'absmax')
      amax = max(max(abs(Zi)));
      amin = -amax;
    elseif strcmp(MAPLIMITS,'maxmin') | strcmp(MAPLIMITS,'minmax')
      amin = min(min(Zi));
      amax = max(max(Zi));
    else
      error('unknown ''maplimits'' value.');
    end
  elseif length(MAPLIMITS) == 2
    amin = MAPLIMITS(1);
    amax = MAPLIMITS(2);
  else
    error('unknown ''maplimits'' value');
  end
  delta = xi(2)-xi(1); % length of grid entry

 end % if ~strcmpi(STYLE,'grid')
  %
  %%%%%%%%%%%%%%%%%%%%%%%%%% Scale the axes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  %cla  % clear current axis
%   hold on
%   h = gca; % uses current axes

                          % instead of default larger AXHEADFAC
  if squeezefac<0.92 & plotrad-headrad > 0.05  % (size of head in axes)
    AXHEADFAC = 1.05;     % do not leave room for external ears if head cartoon
                          % shrunk enough by the 'skirt' option
  end

%   set(gca,'Xlim',[-rmax rmax]*AXHEADFAC,'Ylim',[-rmax rmax]*AXHEADFAC);
                          % specify size of head axes in gca

  unsh = (GRID_SCALE+1)/GRID_SCALE; % un-shrink the effects of 'interp' SHADING

  %
  %%%%%%%%%%%%%%%%%%%%%%%% Plot grid only %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  if strcmpi(STYLE,'grid')                     % plot grid only

    %
    % The goal below is to make the grid cells square - not yet achieved in all cases? -sm
    %
    g1 = size(gridchans,1);
    g2 = size(gridchans,2);
    gmax = max([g1 g2]);
    Xi = linspace(-rmax*g2/gmax,rmax*g2/gmax,g1+1);
    Xi = Xi+rmax/g1; Xi = Xi(1:end-1);
    Yi = linspace(-rmax*g1/gmax,rmax*g1/gmax,g2+1);
    Yi = Yi+rmax/g2; Yi = Yi(1:end-1); Yi = Yi(end:-1:1); % by trial and error!
    %
    %%%%%%%%%%% collect the gridchans values %%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    gridvalues = zeros(size(gridchans));
    for j=1:size(gridchans,1)
      for k=1:size(gridchans,2)
         gc = gridchans(j,k);
         if gc > 0
              gridvalues(j,k) = Values(gc);
         elseif gc < 0
              gridvalues(j,k) = -Values(gc);
         else
              gridvalues(j,k) = nan; % not-a-number = no value
         end
      end
    end
  end
end

    %
    %%%%%%%%%%% reset color limits for grid plot %%%%%%%%%%%%%%%%%%%%%%%%%
    %
%     if isstr(MAPLIMITS)
%       if strcmp(MAPLIMITS,'maxmin') | strcmp(MAPLIMITS,'minmax')
%         amin = min(min(gridvalues(~isnan(gridvalues))));
%         amax = max(max(gridvalues(~isnan(gridvalues))));
%       elseif strcmp(MAPLIMITS,'absmax')
%         % 11/21/2005 Toby edit
%         % This should now work as specified. Before it only crashed (using
%         % "plotgrid" and "maplimits>absmax" options).
%         amax = max(max(abs(gridvalues(~isnan(gridvalues)))));
%         amin = -amax;
%         %amin = -max(max(abs([amin amax])));
%         %amax = max(max(abs([amin amax])));
%       else
%         error('unknown ''maplimits'' value');
%       end
%     elseif length(MAPLIMITS) == 2
%       amin = MAPLIMITS(1);
%       amax = MAPLIMITS(2);
%     else
%       error('unknown ''maplimits'' value');
%     end
%     %
%     %%%%%%%%%% explicitly compute grid colors, allowing BACKCOLOR  %%%%%%
%     %
%     gridvalues = 1+floor(cmaplen*(gridvalues-amin)/(amax-amin));
%     gridvalues(find(gridvalues == cmaplen+1)) = cmaplen;
%     gridcolors = zeros([size(gridvalues),3]);
%     for j=1:size(gridchans,1)
%       for k=1:size(gridchans,2)
%          if ~isnan(gridvalues(j,k))
%              gridcolors(j,k,:) = cmap(gridvalues(j,k),:);
%          else
%             if strcmpi(whitebk,'off')
%                 gridcolors(j,k,:) = BACKCOLOR; % gridchans == 0 -> background color
%                 % This allows the plot to show 'space' between separate sub-grids or strips
%             else % 'on'
%                 gridcolors(j,k,:) = [1 1 1]; BACKCOLOR; % gridchans == 0 -> white for printing
%             end
%          end
%       end
%     end
%
%     %
%     %%%%%%%%%% draw the gridplot image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     %
%     handle=imagesc(Xi,Yi,gridcolors); % plot grid with explicit colors
%     axis square
%   %
%   %%%%%%%%%%%%%%%%%%%%%%%% Plot map contours only %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   %
%   elseif strcmp(STYLE,'contour')                     % plot surface contours only
%     [cls chs] = contour(Xi,Yi,ZiC,CONTOURNUM,'k');
%     handle = chs;                                   % handle to a contourgroup object
%     % for h=chs, set(h,'color',CCOLOR); end
%   %
%   %%%%%%%%%%%%%%%%%%%%%%%% Else plot map and contours %%%%%%%%%%%%%%%%%%%%%%%%%
%   %
%   elseif strcmp(STYLE,'both')  % plot interpolated surface and surface contours
%       if strcmp(SHADING,'interp')
%        tmph = surface(Xi*unsh,Yi*unsh,zeros(size(Zi))-0.1,Zi,...
%                'EdgeColor','none','FaceColor',SHADING);
%     else % SHADING == 'flat'
%        tmph = surface(Xi-delta/2,Yi-delta/2,zeros(size(Zi))-0.1,Zi,...
%                'EdgeColor','none','FaceColor',SHADING);
%     end
%     if strcmpi(MASKSURF, 'on')
%         set(tmph, 'visible', 'off');
%         handle = tmph;
%     end;
%
%     warning off;
%     if ~PMASKFLAG
%         [cls chs] = contour(Xi,Yi,ZiC,CONTOURNUM,'k');
%     else
%         ZiC(find(ZiC > 0.5 )) = NaN;
%         [cls chs] = contourf(Xi,Yi,ZiC,0,'k');
%         subh = get(chs, 'children');
%         for indsubh = 1:length(subh)
%             numfaces = size(get(subh(indsubh), 'XData'),1);
%             set(subh(indsubh), 'FaceVertexCData', ones(numfaces,3), 'Cdatamapping', 'direct', 'facealpha', 0.5, 'linewidth', 2);
%         end;
%     end;
%     handle = tmph;                                   % surface handle
%     for h=chs, set(h,'color',CCOLOR); end
%     warning on;
%   %
%   %%%%%%%%%%%%%%%%%%%%%%%% Else plot map only %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   %
%   elseif strcmp(STYLE,'straight') | strcmp(STYLE,'map') % 'straight' was former arg
%
%       if strcmp(SHADING,'interp') % 'interp' mode is shifted somehow... but how?
%          tmph = surface(Xi*unsh,Yi*unsh,zeros(size(Zi)),Zi,'EdgeColor','none',...
%                   'FaceColor',SHADING);
%       else
%          tmph = surface(Xi-delta/2,Yi-delta/2,zeros(size(Zi)),Zi,'EdgeColor','none',...
%                  'FaceColor',SHADING);
%       end
%     if strcmpi(MASKSURF, 'on')
%         set(tmph, 'visible', 'off');
%         handle = tmph;
%     end;
%     handle = tmph;                                   % surface handle
%   %
%   %%%%%%%%%%%%%%%%%% Else fill contours with uniform colors  %%%%%%%%%%%%%%%%%%
%   %
%   elseif strcmp(STYLE,'fill')
%     [cls chs] = contourf(Xi,Yi,Zi,CONTOURNUM,'k');
%
%     handle = chs;                                   % handle to a contourgroup object
%
%     % for h=chs, set(h,'color',CCOLOR); end
%     %     <- 'not line objects.' Why does 'both' work above???
%
%   else
%     error('Invalid style')
%   end
%   %
%   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Set color axis  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   %
% %   caxis([amin amax]); % set coloraxis
%
% % 7/30/2014 Ramon: +-5% for the color limits were added
% cax_sgn = sign([amin amax]);                                                  % getting sign
% caxis([amin+cax_sgn(1)*(0.05*abs(amin)) amax+cax_sgn(2)*(0.05*abs(amax))]);   % Adding 5% to the color limits
%
% else % if STYLE 'blank'
% %
% %%%%%%%%%%%%%%%%%%%%%%% Draw blank head %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
%   if strcmpi(noplot, 'on')
%    if strcmpi(VERBOSE,'on')
%       fprintf('topoplot(): no plot requested.\n')
%    end
%    return;
%   end
%   %cla
%   hold on
%
%   set(gca,'Xlim',[-rmax rmax]*AXHEADFAC,'Ylim',[-rmax rmax]*AXHEADFAC)
%    % pos = get(gca,'position');
%    % fprintf('Current axes size %g,%g\n',pos(3),pos(4));
%
%   if strcmp(ELECTRODES,'labelpoint') |  strcmp(ELECTRODES,'numpoint')
%     text(-0.6,-0.6, ...
%     [ int2str(length(Rd)) ' of ' int2str(length(tmpeloc)) ' electrode locations shown']);
%     text(-0.6,-0.7, [ 'Click on electrodes to toggle name/number']);
%     tl = title('Channel locations');
%     set(tl, 'fontweight', 'bold');
%   end;
% end % STYLE 'blank'
%
% if exist('handle') ~= 1
%     handle = gca;
% end;
%
% if ~strcmpi(STYLE,'grid')                     % if not plot grid only
%
% %
% %%%%%%%%%%%%%%%%%%% Plot filled ring to mask jagged grid boundary %%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% hwidth = HEADRINGWIDTH;                   % width of head ring
% hin  = squeezefac*headrad*(1- hwidth/2);  % inner head ring radius
%
% if strcmp(SHADING,'interp')
%   rwidth = BLANKINGRINGWIDTH*1.3;             % width of blanking outer ring
% else
%   rwidth = BLANKINGRINGWIDTH;         % width of blanking outer ring
% end
% rin    =  rmax*(1-rwidth/2);              % inner ring radius
% if hin>rin
%   rin = hin;                              % dont blank inside the head ring
% end
%
% if strcmp(CONVHULL,'on') %%%%%%%%% mask outside the convex hull of the electrodes %%%%%%%%%
%   cnv = convhull(allx,ally);
%   cnvfac = round(CIRCGRID/length(cnv)); % spline interpolate the convex hull
%   if cnvfac < 1, cnvfac=1; end;
%   CIRCGRID = cnvfac*length(cnv);
%
%   startangle = atan2(allx(cnv(1)),ally(cnv(1)));
%   circ = linspace(0+startangle,2*pi+startangle,CIRCGRID);
%   rx = sin(circ);
%   ry = cos(circ);
%
%   allx = allx(:)';  % make x (elec locations; + to nose) a row vector
%   ally = ally(:)';  % make y (elec locations, + to r? ear) a row vector
%   erad = sqrt(allx(cnv).^2+ally(cnv).^2);  % convert to polar coordinates
%   eang = atan2(allx(cnv),ally(cnv));
%   eang = unwrap(eang);
%   eradi =spline(linspace(0,1,3*length(cnv)), [erad erad erad], ...
%                                       linspace(0,1,3*length(cnv)*cnvfac));
%   eangi =spline(linspace(0,1,3*length(cnv)), [eang+2*pi eang eang-2*pi], ...
%                                       linspace(0,1,3*length(cnv)*cnvfac));
%   xx = eradi.*sin(eangi);           % convert back to rect coordinates
%   yy = eradi.*cos(eangi);
%   yy = yy(CIRCGRID+1:2*CIRCGRID);
%   xx = xx(CIRCGRID+1:2*CIRCGRID);
%   eangi = eangi(CIRCGRID+1:2*CIRCGRID);
%   eradi = eradi(CIRCGRID+1:2*CIRCGRID);
%   xx = xx*1.02; yy = yy*1.02;           % extend spline outside electrode marks
%
%   splrad = sqrt(xx.^2+yy.^2);           % arc radius of spline points (yy,xx)
%   oob = find(splrad >= rin);            %  enforce an upper bound on xx,yy
%   xx(oob) = rin*xx(oob)./splrad(oob);   % max radius = rin
%   yy(oob) = rin*yy(oob)./splrad(oob);   % max radius = rin
%
%   splrad = sqrt(xx.^2+yy.^2);           % arc radius of spline points (yy,xx)
%   oob = find(splrad < hin);             % don't let splrad be inside the head cartoon
%   xx(oob) = hin*xx(oob)./splrad(oob);   % min radius = hin
%   yy(oob) = hin*yy(oob)./splrad(oob);   % min radius = hin
%
%   ringy = [[ry(:)' ry(1) ]*(rin+rwidth) yy yy(1)];
%   ringx = [[rx(:)' rx(1) ]*(rin+rwidth) xx xx(1)];
%
%   ringh2= patch(ringy,ringx,ones(size(ringy)),BACKCOLOR,'edgecolor','none'); hold on
%
%   % plot(ry*rmax,rx*rmax,'b') % debugging line
%
% else %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% mask the jagged border around rmax %%%%%%%%%%%%%%%5%%%%%%
%
%   circ = linspace(0,2*pi,CIRCGRID);
%   rx = sin(circ);
%   ry = cos(circ);
%   ringx = [[rx(:)' rx(1) ]*(rin+rwidth)  [rx(:)' rx(1)]*rin];
%   ringy = [[ry(:)' ry(1) ]*(rin+rwidth)  [ry(:)' ry(1)]*rin];
%
%   if ~strcmpi(STYLE,'blank')
%     ringh= patch(ringx,ringy,0.01*ones(size(ringx)),BACKCOLOR,'edgecolor','none'); hold on
%   end
%   % plot(ry*rmax,rx*rmax,'b') % debugging line
% end
%
%   %f1= fill(rin*[rx rX],rin*[ry rY],BACKCOLOR,'edgecolor',BACKCOLOR); hold on
%   %f2= fill(rin*[rx rX*(1+rwidth)],rin*[ry rY*(1+rwidth)],BACKCOLOR,'edgecolor',BACKCOLOR);
%
% % Former line-style border smoothing - width did not scale with plot
% %  brdr=plot(1.015*cos(circ).*rmax,1.015*sin(circ).*rmax,...      % old line-based method
% %      'color',HEADCOLOR,'Linestyle','-','LineWidth',HLINEWIDTH);    % plot skirt outline
% %  set(brdr,'color',BACKCOLOR,'linewidth',HLINEWIDTH + 4);        % hide the disk edge jaggies
%
% %
% %%%%%%%%%%%%%%%%%%%%%%%%% Plot cartoon head, ears, nose %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% if headrad > 0                         % if cartoon head to be plotted
% %
% %%%%%%%%%%%%%%%%%%% Plot head outline %%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% headx = [[rx(:)' rx(1) ]*(hin+hwidth)  [rx(:)' rx(1)]*hin];
% heady = [[ry(:)' ry(1) ]*(hin+hwidth)  [ry(:)' ry(1)]*hin];
%
% if ~isstr(HEADCOLOR) | ~strcmpi(HEADCOLOR,'none')
%    %ringh= patch(headx,heady,ones(size(headx)),HEADCOLOR,'edgecolor',HEADCOLOR,'linewidth', HLINEWIDTH); hold on
%    headx = [rx(:)' rx(1)]*hin;
%    heady = [ry(:)' ry(1)]*hin;
%    ringh= plot(headx,heady);
%    set(ringh, 'color',HEADCOLOR,'linewidth', HLINEWIDTH); hold on
% end
%
% % rx = sin(circ); rX = rx(end:-1:1);
% % ry = cos(circ); rY = ry(end:-1:1);
% % for k=2:2:CIRCGRID
% %   rx(k) = rx(k)*(1+hwidth);
% %   ry(k) = ry(k)*(1+hwidth);
% % end
% % f3= fill(hin*[rx rX],hin*[ry rY],HEADCOLOR,'edgecolor',HEADCOLOR); hold on
% % f4= fill(hin*[rx rX*(1+hwidth)],hin*[ry rY*(1+hwidth)],HEADCOLOR,'edgecolor',HEADCOLOR);
%
% % Former line-style head
% %  plot(cos(circ).*squeezefac*headrad,sin(circ).*squeezefac*headrad,...
% %      'color',HEADCOLOR,'Linestyle','-','LineWidth',HLINEWIDTH);    % plot head outline
%
% %
% %%%%%%%%%%%%%%%%%%% Plot ears and nose %%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
%   base  = rmax-.0046;
%   basex = 0.18*rmax;                   % nose width
%   tip   = 1.15*rmax;
%   tiphw = .04*rmax;                    % nose tip half width
%   tipr  = .01*rmax;                    % nose tip rounding
%   q = .04; % ear lengthening
%   EarX  = [.497-.005  .510  .518  .5299 .5419  .54    .547   .532   .510   .489-.005]; % rmax = 0.5
%   EarY  = [q+.0555 q+.0775 q+.0783 q+.0746 q+.0555 -.0055 -.0932 -.1313 -.1384 -.1199];
%   sf    = headrad/plotrad;                                          % squeeze the model ears and nose
%                                                                     % by this factor
%   if ~isstr(HEADCOLOR) | ~strcmpi(HEADCOLOR,'none')
%     plot3([basex;tiphw;0;-tiphw;-basex]*sf,[base;tip-tipr;tip;tip-tipr;base]*sf,...
%          2*ones(size([basex;tiphw;0;-tiphw;-basex])),...
%          'Color',HEADCOLOR,'LineWidth',HLINEWIDTH);                 % plot nose
%     plot3(EarX*sf,EarY*sf,2*ones(size(EarX)),'color',HEADCOLOR,'LineWidth',HLINEWIDTH)    % plot left ear
%     plot3(-EarX*sf,EarY*sf,2*ones(size(EarY)),'color',HEADCOLOR,'LineWidth',HLINEWIDTH)   % plot right ear
%   end
% end
%
% %
% % %%%%%%%%%%%%%%%%%%% Show electrode information %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
%  plotax = gca;
%  axis square                                           % make plotax square
%  axis off
%
%  pos = get(gca,'position');
%  xlm = get(gca,'xlim');
%  ylm = get(gca,'ylim');
%  % textax = axes('position',pos,'xlim',xlm,'ylim',ylm);  % make new axes so clicking numbers <-> labels
%                                                        % will work inside head cartoon patch
%  % axes(textax);
%  axis square                                           % make textax square
%
%  pos = get(gca,'position');
%  set(plotax,'position',pos);
%
%  xlm = get(gca,'xlim');
%  set(plotax,'xlim',xlm);
%
%  ylm = get(gca,'ylim');
%  set(plotax,'ylim',ylm);                               % copy position and axis limits again
%
% axis equal;
% set(gca, 'xlim', [-0.525 0.525]); set(plotax, 'xlim', [-0.525 0.525]);
% set(gca, 'ylim', [-0.525 0.525]); set(plotax, 'ylim', [-0.525 0.525]);
%
% %get(textax,'pos')    % test if equal!
% %get(plotax,'pos')
% %get(textax,'xlim')
% %get(plotax,'xlim')
% %get(textax,'ylim')
% %get(plotax,'ylim')
%
%  if isempty(EMARKERSIZE)
%    EMARKERSIZE = 10;
%    if length(y)>=160
%     EMARKERSIZE = 3;
%    elseif length(y)>=128
%     EMARKERSIZE = 3;
%    elseif length(y)>=100
%     EMARKERSIZE = 3;
%    elseif length(y)>=80
%     EMARKERSIZE = 4;
%    elseif length(y)>=64
%     EMARKERSIZE = 5;
%    elseif length(y)>=48
%     EMARKERSIZE = 6;
%    elseif length(y)>=32
%     EMARKERSIZE = 8;
%    end
%  end
% %
% %%%%%%%%%%%%%%%%%%%%%%%% Mark electrode locations only %%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% ELECTRODE_HEIGHT = 2.1;  % z value for plotting electrode information (above the surf)
%
% if strcmp(ELECTRODES,'on')   % plot electrodes as spots
%   if isempty(EMARKER2CHANS)
%     hp2 = plot3(y,x,ones(size(x))*ELECTRODE_HEIGHT,...
%         EMARKER,'Color',ECOLOR,'markersize',EMARKERSIZE,'linewidth',EMARKERLINEWIDTH);
%   else % plot markers for normal chans and EMARKER2CHANS separately
%     hp2 = plot3(y(mark1chans),x(mark1chans),ones(size((mark1chans)))*ELECTRODE_HEIGHT,...
%         EMARKER,'Color',ECOLOR,'markersize',EMARKERSIZE,'linewidth',EMARKERLINEWIDTH);
%     hp2b = plot3(y(mark2chans),x(mark2chans),ones(size((mark2chans)))*ELECTRODE_HEIGHT,...
%         EMARKER2,'Color',EMARKER2COLOR,'markerfacecolor',EMARKER2COLOR,'linewidth',EMARKER2LINEWIDTH,'markersize',EMARKERSIZE2);
%   end
% %
% %%%%%%%%%%%%%%%%%%%%%%%% Print electrode labels only %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% elseif strcmp(ELECTRODES,'labels')  % print electrode names (labels)
%     for i = 1:size(labels,1)
%     text(double(y(i)),double(x(i)),...
%         ELECTRODE_HEIGHT,labels(i,:),'HorizontalAlignment','center',...
% 	'VerticalAlignment','middle','Color',ECOLOR,...
% 	'FontSize',EFSIZE)
%   end
% %
% %%%%%%%%%%%%%%%%%%%%%%%% Mark electrode locations plus labels %%%%%%%%%%%%%%%%%%%
% %
% elseif strcmp(ELECTRODES,'labelpoint')
%   if isempty(EMARKER2CHANS)
%     hp2 = plot3(y,x,ones(size(x))*ELECTRODE_HEIGHT,...
%         EMARKER,'Color',ECOLOR,'markersize',EMARKERSIZE,'linewidth',EMARKERLINEWIDTH);
%   else
%     hp2 = plot3(y(mark1chans),x(mark1chans),ones(size((mark1chans)))*ELECTRODE_HEIGHT,...
%         EMARKER,'Color',ECOLOR,'markersize',EMARKERSIZE,'linewidth',EMARKERLINEWIDTH);
%     hp2b = plot3(y(mark2chans),x(mark2chans),ones(size((mark2chans)))*ELECTRODE_HEIGHT,...
%         EMARKER2,'Color',EMARKER2COLOR,'markerfacecolor',EMARKER2COLOR,'linewidth',EMARKER2LINEWIDTH,'markersize',EMARKERSIZE2);
%   end
%   for i = 1:size(labels,1)
%     hh(i) = text(double(y(i)+0.01),double(x(i)),...
%         ELECTRODE_HEIGHT,labels(i,:),'HorizontalAlignment','left',...
% 	'VerticalAlignment','middle','Color', ECOLOR,'userdata', num2str(allchansind(i)), ...
% 	'FontSize',EFSIZE, 'buttondownfcn', ...
% 	    ['tmpstr = get(gco, ''userdata'');'...
% 	     'set(gco, ''userdata'', get(gco, ''string''));' ...
% 	     'set(gco, ''string'', tmpstr); clear tmpstr;'] );
%   end
% %
% %%%%%%%%%%%%%%%%%%%%%%% Mark electrode locations plus numbers %%%%%%%%%%%%%%%%%%%
% %
% elseif strcmp(ELECTRODES,'numpoint')
%   if isempty(EMARKER2CHANS)
%     hp2 = plot3(y,x,ones(size(x))*ELECTRODE_HEIGHT,...
%         EMARKER,'Color',ECOLOR,'markersize',EMARKERSIZE,'linewidth',EMARKERLINEWIDTH);
%   else
%     hp2 = plot3(y(mark1chans),x(mark1chans),ones(size((mark1chans)))*ELECTRODE_HEIGHT,...
%         EMARKER,'Color',ECOLOR,'markersize',EMARKERSIZE,'linewidth',EMARKERLINEWIDTH);
%     hp2b = plot3(y(mark2chans),x(mark2chans),ones(size((mark2chans)))*ELECTRODE_HEIGHT,...
%         EMARKER2,'Color',EMARKER2COLOR,'markerfacecolor',EMARKER2COLOR,'linewidth',EMARKER2LINEWIDTH,'markersize',EMARKERSIZE2);
%   end
%   for i = 1:size(labels,1)
%     hh(i) = text(double(y(i)+0.01),double(x(i)),...
%         ELECTRODE_HEIGHT,num2str(allchansind(i)),'HorizontalAlignment','left',...
% 	'VerticalAlignment','middle','Color', ECOLOR,'userdata', labels(i,:) , ...
% 	'FontSize',EFSIZE, 'buttondownfcn', ...
% 	    ['tmpstr = get(gco, ''userdata'');'...
% 	     'set(gco, ''userdata'', get(gco, ''string''));' ...
% 	     'set(gco, ''string'', tmpstr); clear tmpstr;'] );
%   end
% %
% %%%%%%%%%%%%%%%%%%%%%% Print electrode numbers only %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% elseif strcmp(ELECTRODES,'numbers')
%   for i = 1:size(labels,1)
%     text(double(y(i)),double(x(i)),...
%         ELECTRODE_HEIGHT,int2str(allchansind(i)),'HorizontalAlignment','center',...
% 	'VerticalAlignment','middle','Color',ECOLOR,...
% 	'FontSize',EFSIZE)
%   end
% %
% %%%%%%%%%%%%%%%%%%%%%% Mark emarker2 electrodes only  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% elseif strcmp(ELECTRODES,'off') & ~isempty(EMARKER2CHANS)
%     hp2b = plot3(y(mark2chans),x(mark2chans),ones(size((mark2chans)))*ELECTRODE_HEIGHT,...
%         EMARKER2,'Color',EMARKER2COLOR,'markerfacecolor',EMARKER2COLOR,'linewidth',EMARKER2LINEWIDTH,'markersize',EMARKERSIZE2);
% end
% %
% %%%%%%%% Mark specified electrode locations with red filled disks  %%%%%%%%%%%%%%%%%%%%%%
% %
% try,
%     if strcmpi(STYLE,'blank') % if mark-selected-channel-locations mode
%         for kk = 1:length(1:length(x))
%             if abs(Values(kk))
%                 if PLOTDISK
%                     angleRatio = real(Values(kk))/(real(Values(kk))+imag(Values(kk)))*360;
%                     radius     = real(Values(kk))+imag(Values(kk));
%                     allradius  = [0.02 0.03 0.037 0.044 0.05];
%                     radius     = allradius(radius);
%                     hp2 = disk(y(kk),x(kk),radius, [1 0 0], 0 , angleRatio, 16);
%                     if angleRatio ~= 360
%                         hp2 = disk(y(kk),x(kk),radius, [0 0 1], angleRatio, 360, 16);
%                     end;
%                 else
%                     tmpcolor = COLORARRAY{max(1,min(Values(kk), length(COLORARRAY)))};
%                     hp2 = plot3(y(kk),x(kk),ELECTRODE_HEIGHT,EMARKER,'Color', tmpcolor, 'markersize', EMARKERSIZE1CHAN);
%                     hp2 = disk(y(kk),x(kk),0.04, tmpcolor, 0, 360, 10);
%                 end;
%             end;
%         end
%     end
% catch, end;
% %
% %%%%%%%%%%%%%%%%%%%%%%%%%%% Plot dipole(s) on the scalp map  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% if ~isempty(DIPOLE)
%     hold on;
%     tmp = DIPOLE;
%     if isstruct(DIPOLE)
%         if ~isfield(tmp,'posxyz')
%            error('dipole structure is not an EEG.dipfit.model')
%         end
%         DIPOLE = [];  % Note: invert x and y from dipplot usage
%         DIPOLE(:,1) = -tmp.posxyz(:,2)/DIPSPHERE; % -y -> x
%         DIPOLE(:,2) =  tmp.posxyz(:,1)/DIPSPHERE; %  x -> y
%         DIPOLE(:,3) = -tmp.momxyz(:,2);
%         DIPOLE(:,4) =  tmp.momxyz(:,1);
%     else
%         DIPOLE(:,1) = -tmp(:,2);                    % same for vector input
%         DIPOLE(:,2) =  tmp(:,1);
%         DIPOLE(:,3) = -tmp(:,4);
%         DIPOLE(:,4) =  tmp(:,3);
%     end;
%     for index = 1:size(DIPOLE,1)
%         if ~any(DIPOLE(index,:))
%              DIPOLE(index,:) = [];
%         end
%     end;
%     DIPOLE(:,1:4)   = DIPOLE(:,1:4)*rmax*(rmax/plotrad); % scale radius from 1 -> rmax (0.5)
%     DIPOLE(:,3:end) = (DIPOLE(:,3:end))*rmax/100000*(rmax/plotrad);
%     if strcmpi(DIPNORM, 'on')
%         for index = 1:size(DIPOLE,1)
%             DIPOLE(index,3:4) = DIPOLE(index,3:4)/norm(DIPOLE(index,3:end))*0.2;
%         end;
%     end;
%     DIPOLE(:, 3:4) =  DIPORIENT*DIPOLE(:, 3:4)*DIPLEN;
%
%     PLOT_DIPOLE=1;
%     if sum(DIPOLE(1,3:4).^2) <= 0.00001
%       if strcmpi(VERBOSE,'on')
%         fprintf('Note: dipole is length 0 - not plotted\n')
%       end
%       PLOT_DIPOLE = 0;
%     end
%     if 0 % sum(DIPOLE(1,1:2).^2) > plotrad
%       if strcmpi(VERBOSE,'on')
%         fprintf('Note: dipole is outside plotting area - not plotted\n')
%       end
%       PLOT_DIPOLE = 0;
%     end
%     if PLOT_DIPOLE
%       for index = 1:size(DIPOLE,1)
%         hh = plot( DIPOLE(index, 1), DIPOLE(index, 2), '.');
%         set(hh, 'color', DIPCOLOR, 'markersize', DIPSCALE*30);
%         hh = line( [DIPOLE(index, 1) DIPOLE(index, 1)+DIPOLE(index, 3)]', ...
%                    [DIPOLE(index, 2) DIPOLE(index, 2)+DIPOLE(index, 4)]',[10 10]);
%         set(hh, 'color', DIPCOLOR, 'linewidth', DIPSCALE*30/7);
%       end;
%     end;
% end;
%
% end % if ~ 'gridplot'
%
% %
% %%%%%%%%%%%%% Plot axis orientation %%%%%%%%%%%%%%%%%%%%
% %
% if strcmpi(DRAWAXIS, 'on')
%     axes('position', [0 0.85 0.08 0.1]);
%     axis off;
%     coordend1 = sqrt(-1)*3;
%     coordend2 = -3;
%     coordend1 = coordend1*exp(sqrt(-1)*rotate);
%     coordend2 = coordend2*exp(sqrt(-1)*rotate);
%
%     line([5 5+round(real(coordend1))]', [5 5+round(imag(coordend1))]', 'color', 'k');
%     line([5 5+round(real(coordend2))]', [5 5+round(imag(coordend2))]', 'color', 'k');
%     if round(real(coordend2))<0
%          text( 5+round(real(coordend2))*1.2, 5+round(imag(coordend2))*1.2-2, '+Y');
%     else text( 5+round(real(coordend2))*1.2, 5+round(imag(coordend2))*1.2, '+Y');
%     end;
%     if round(real(coordend1))<0
%          text( 5+round(real(coordend1))*1.2, 5+round(imag(coordend1))*1.2+1.5, '+X');
%     else text( 5+round(real(coordend1))*1.2, 5+round(imag(coordend1))*1.2, '+X');
%     end;
%     set(gca, 'xlim', [0 10], 'ylim', [0 10]);
% end;
%
% %
% %%%%%%%%%%%%% Set EEGLAB background color to match head border %%%%%%%%%%%%%%%%%%%%%%%%
% %
% try,
%   icadefs;
%   set(gcf, 'color', BACKCOLOR);
%   catch,
% end;
%
% hold off
% axis off
% return
%
% %
% %%%%%%%%%%%%% Draw circle %%%%%%%%%%%%%%%%%%%%%%%%
% %
% function h2 = disk(X, Y, radius, colorfill, oriangle, endangle, segments)
% 	A = linspace(oriangle/180*pi, endangle/180*pi, segments-1);
%     if endangle-oriangle == 360
%      	 A  = linspace(oriangle/180*pi, endangle/180*pi, segments);
%          h2 = patch( [X   + cos(A)*radius(1)], [Y   + sin(A)*radius(end)], zeros(1,segments)+3, colorfill);
%     else A  = linspace(oriangle/180*pi, endangle/180*pi, segments-1);
%          h2 = patch( [X X + cos(A)*radius(1)], [Y Y + sin(A)*radius(end)], zeros(1,segments)+3, colorfill);
%     end;
%     set(h2, 'FaceColor', colorfill);
% 	set(h2, 'EdgeColor', 'none');


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% Laplace Transform properties

% Initial and Final Value Theorems

%Initial Value Theorem
x=sin(t);
X=laplace(x,s);
limit(s*X,s,inf)

%Final Value Theorem
x=1-exp(-t);
limit(x,t,inf)

X=laplace(x,s);
limit(s*X,s,0)

author: Sable

repo: mcbench-benchmarks


text:
classdef MAC < dagnn.Filter
% MAC layer that computes global Max Pooling of convolution activations.
%
% Authors: F. Radenovic, G. Tolias, O. Chum. 2017.

  methods
    function outputs = forward(self, inputs, params)
      outputs{1} = vl_nnpool(inputs{1}, [size(inputs{1},1), size(inputs{1},2)], 'method', 'max', 'cuDNN');
    end

    function [derInputs, derParams] = backward(self, inputs, params, derOutputs)
      derInputs{1} = vl_nnpool(inputs{1}, [size(inputs{1},1), size(inputs{1},2)], derOutputs{1}, 'method', 'max', 'cuDNN');
      derParams = {};
    end

    function obj = MAC(varargin)
      obj.load(varargin);
    end
  end
end

author: filipradenovic

repo: cnnimageretrieval


text:
%MDL_STANFORD_MDH Create model of Stanford arm using MDH conventions
%
% MDL_STANFORD is a script that creates the workspace variable stanf which
% describes the kinematic and dynamic characteristics of the Stanford
% (Scheinman) arm using modified Denavit-Hartenberg parameters.
%
% Also defines the vectors:
%   qz   zero joint angle configuration.
%
% Notes::
% - SI units are used.
%
% References::
% - Kinematic data from "Modelling, Trajectory calculation and Servoing of
%   a computer controlled arm".  Stanford AIM-177.  Figure 2.3
% - Dynamic data from "Robot manipulators: mathematics, programming and control"
%   Paul 1981, Tables 6.5, 6.6
%
% See also mdl_puma560, mdl_puma560akb, SerialLink.

% MODEL: Stanford, Stanford arm, prismatic, 6DOF, modified_DH

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

clear L

qz = [0 0 0 0 0 0];

stanf = SerialLink([
            RevoluteMDH('d', 0.412)
            RevoluteMDH('d', 0.154, 'alpha', -pi/2)
            PrismaticMDH('alpha', pi/2, 'qlim', [0.2032 0.9144])
            RevoluteMDH()
            RevoluteMDH('alpha', -pi/2)
            RevoluteMDH('d', 0.263, 'alpha', pi/2)
        ], ...
    'name', 'Stanford arm MDH', ...
    'plotopt', {'workspace', [-2 2 -2 2 -2 2]} ...
    );

author: petercorke

repo: robotics-toolbox-matlab


text:
function [decisions, decscores] = slclassify(scores, clabels, op, varargin)
%SLCLASSIFY Classifies a set of samples according to final scores
%
% $ Syntax $
%   - decisions = slclassify(scores, clabels, op, ...)
%   - [decisions, decscores] = slclassify(scores, clabels, op, ...)
%
% $ Arguments $
%   - scores:       the score matrix
%   - clabels:      the class labels of reference samples
%   - op:           the score attribute
%   - decisions:    the classification decisions
%   - decscores:    the scores of the classified targets
%
% $ Description $
%   - decisions = slclassify(scores, clabels, op, ...) classifies
%     a set of query samples to classes. Suppose there are m referenced
%     targets and n query samples, scores should be an m x n matrix, with
%     each column representing the scores of the corresponding sample to
%     all targets.
%     If op is 'high', then the samples are classified to the target of
%     highest score, otherwise, the samples are classified to the target
%     of lowest score.
%     Moreover, you can specify following properties to have more
%     control on the classification process.
%       \*
%       \t   Properties of Classification
%       \h     name     &  description
%             'scheme'  &  The classification scheme
%                          'nn':  using normal nearest sample
%                                 classification (default)
%                          'loo': leave-one-out nearest sample scheme
%                                 (only for the case, when gallery and
%                                  query sets are the same)
%       \*
%
%   - [decisions, decscores] = slclassify(scores, clabels, op, ...)
%     additionally outputs the scores on classified samples.
%
% $ Remarks $
%   - The outputs are 1xn row vectors.
%   - In leave-one-out scheme, it is assumed that the referenced samples
%     and the query samples are actually the same set and in same order.
%
% $ History $
%   - Created by Dahua Lin, on Aug 9th, 2006
%   - Modified by Dahua Lin, on Aug 16th, 2006
%       - eliminate the qlabel parameters, which are essentially not
%         needed.
%       - add functionality to support schemes. In current revision,
%         it supports nearest-neighbor ('nn') and leave-one-out ('loo').
%

%% parse and verify input arguments

if nargin < 3
    raise_lackinput('slclassify', 3);
end
if ~isnumeric(scores) || ndims(scores) ~= 2
    error('sltoolbox:invalidarg', ...
        'scores should be an 2D numeric matrix');
end
m = size(scores, 1);
if length(clabels) ~= m
    error('sltoolbox:sizmismatch', ...
        'The sizes of labels are inconsistent with the score matrix');
end

opts.scheme = 'nn';
opts = slparseprops(opts, varargin{:});

%% Main skeleton

switch opts.scheme
    case 'nn'
        [decisions, decscores] = classify_nn(scores, clabels, op);
    case 'loo'
        [decisions, decscores] = classify_loo(scores, clabels, op);
    otherwise
        error('sltoolbox:invalidarg', ...
            'Invalid scheme for classification: %s', opts.scheme);
end

%% Decision-making core routines

%% NN

function [decisions, decscores] = classify_nn(scores, clabels, op)

switch op
    case 'high'
        [decscores, decinds] = max(scores, [], 1);
    case 'low'
        [decscores, decinds] = min(scores, [], 1);
    otherwise
        error('sltoolbox:invalidarg', ...
            'Invalid score option %s', op);
end

clabels = clabels(:)';
decisions = clabels(decinds);

%% LOO

function [decisions, decscores] = classify_loo(scores, clabels, op)

n = size(scores, 1);
if size(scores, 2) ~= n
    error('sltoolbox:sizmismatch', ...
        'In leave-one-out scheme, the score matrix should be square');
end
if n < 2
    error('sltoolbox:invalidarg', ...
        'In leave-one-out scheme, the set should have at least two elements');
end

% preprocessing to disable the selection of self
inds_diag = (1 : n+1 : n^2);
switch op
    case 'high'
        scores(inds_diag) = -Inf;
    case 'low'
        scores(inds_diag) = Inf;
    otherwise
        error('sltoolbox:invalidarg', ...
            'Invalid score option %s', op);
end

% classify
switch op
    case 'high'
        [decscores, decinds] = max(scores, [], 1);
    case 'low'
        [decscores, decinds] = min(scores, [], 1);
end

clabels = clabels(:)';
decisions = clabels(decinds);

author: lmthang


text:
% sim_spinecho_shaped.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% sim_spinecho_shaped(n,sw,Bfield,linewidth,sys,TE,RF,Tp,grad,pos,ph)
%
% DESCRIPTION:
% This function simulates a localized spin-echo sequence with a shaped
% refocusing pulse.  It enables choice of the echo-time as well as the
% choice of the phase refocusing pulse.  This allows phase cycling of the
% refocusing pulses by repeating simulations with different pulse phases,
% which is necessary to remove unwanted coherences from outside the volume
% of interest.  For the refocusing pulse, a two step phase cycling scheme
% is typically sufficient, where the refocusing pulse is phase cycled by 0
% and 90 degrees the phase are combined by subtraction.
%
% INPUTS:
% n          = number of points in fid/spectrum
% sw         = desired spectral width in [Hz]
% Bfield     = main magnetic field strength in [T]
% linewidth  = linewidth in [Hz]
% sys        = Metabolite spin system definition structure;
% TE         = Echo time in [ms]
% RF         = RF pulse definition structure for refoc pulse (obtain using 'io_loadRFwaveform.m');
% Tp         = duration of refocusing pulse in [ms]
% grad       = gradient strength for the selective refocusing pulse [G/cm]
% pos        = position offset in the direction corresponding to the refocusing pulse [cm]
% ph         = the phase of the refocusing pulse in [degrees];
%
% OUTPUTS:
% out       = simulated spectrum, in FID-A structure format, using spin-echo
%             sequence.

function out = sim_spinecho_shaped(n,sw,Bfield,linewidth,sys,TE,RF,Tp,grad,pos,ph)

%Set water to centre
centreFreq=4.65;
for k=1:length(sys)
    sys(k).shifts=sys(k).shifts-centreFreq;
end

%Calculate Hamiltonian matrices and starting density matrix.
[H,d]=sim_Hamiltonian(sys,Bfield);

%Calculate new delay by subtracting the pulse durations from TE;
delay=TE-Tp;
if delay<0
    error(['ERROR!  TE is too short.']);
end

%BEGIN PULSE SEQUENCE************
d=sim_excite(d,H,'x');                            %EXCITE
d=sim_evolve(d,H,delay/2000);                      %Evolve by delay/2
d=sim_shapedRF(d,H,RF,Tp,180,90+ph,pos,grad);   %shaped 180 degree refocusing pulse about y' axis.
d=sim_evolve(d,H,delay/2000);                      %Evolve by delay/2
[out,dout]=sim_readout(d,H,n,sw,linewidth,90);  %Readout along y (90 degree phase);
%END PULSE SEQUENCE**************

%Correct the ppm scale:
out.ppm=out.ppm-(4.65-centreFreq);

%Fill in structure header fields:
out.seq='spinecho';
out.te=TE;
out.sim='shaped';

%Additional fields for compatibility with FID-A processing tools.
out.sz=size(out.specs);
out.date=date;
out.dims.t=1;
out.dims.coils=0;
out.dims.averages=0;
out.dims.subSpecs=0;
out.dims.extras=0;
out.averages=1;
out.rawAverages=1;
out.subspecs=1;
out.rawSubspecs=1;
out.flags.writtentostruct=1;
out.flags.gotparams=1;
out.flags.leftshifted=0;
out.flags.filtered=0;
out.flags.zeropadded=0;
out.flags.freqcorrected=0;
out.flags.phasecorrected=0;
out.flags.averaged=1;
out.flags.addedrcvrs=1;
out.flags.subtracted=1;
out.flags.writtentotext=0;
out.flags.downsampled=0;
out.flags.isFourSteps=0;

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
% Copyright 2017 Lime Microsystems Ltd.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%    http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%
% Compile Oct functions to speed up Octave
%
% run build.m
%
printf "(c) 2017 Lime Microsystems Ltd.\n"
printf "Compiling Oct functions...\n"
if version >= "4.0.0"
	mkoctfile WCDMADLmake.cpp WCDMAlib.cpp FFTlib.cpp -O3
	mkoctfile WCDMADLdemod.cpp WCDMAlib.cpp FFTlib.cpp -O3
	mkoctfile WCDMAULmake.cpp WCDMAlib.cpp FFTlib.cpp -O3
	mkoctfile WCDMAULdemod.cpp WCDMAlib.cpp FFTlib.cpp -O3
else % 3.8.2 or earlier, no -O3 flag support
	mkoctfile WCDMADLmake.cpp WCDMAlib.cpp FFTlib.cpp
	mkoctfile WCDMADLdemod.cpp WCDMAlib.cpp FFTlib.cpp
	mkoctfile WCDMAULmake.cpp WCDMAlib.cpp FFTlib.cpp
	mkoctfile WCDMAULdemod.cpp WCDMAlib.cpp FFTlib.cpp
end
printf "cleaning up...\n"
delete *.o
delete *~
printf "Done.\n"

author: myriadrf


text:
function [ ccd_out, phase_out ] = ccdmem( reference_image, match_image, corr_window_size )
%CCDMEM Coherent change detection
%    ccd_out = ccdmem(reference_image, match_image, corr_window_size)
%
% Compares two images by performing coherent change detection.  Implements
% CCD equation as described in Jakowatz, et al., "Spotlight-mode Synthetic
% Aperture Radar: A Signal Processing Approach".  Assumes images (and all
% processing) fit into memory and that input images are already registered.
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

conjf_times_g = conv2(conj(reference_image).*match_image, ones(corr_window_size), 'same');
f_squared = conv2(abs(reference_image).^2, ones(corr_window_size), 'same');
g_squared = conv2(abs(match_image).^2, ones(corr_window_size), 'same');
ccd_out = abs(conjf_times_g)./sqrt(f_squared.*g_squared); % Jakowatz, eq. 5.102
ccd_out(~isfinite(ccd_out))=0; % Fix divide by zeros
if nargout>1
    phase_out = angle(conjf_times_g);
    phase_out(~isfinite(phase_out))=0; % Fix divide by zeros
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function W = nullSpaceOperatorApply(nullS, V)
% Computes a sparse matrix `W` from `V` such that :math:`S W = 0`.
% First, `nullS = nullSpaceOperator(S)`
% computes a structure `nullS` from an `m x n` sparse matrix `S` :math:`(m < n)`, with rank `r`.
%
% Second, if `V` is an `(n-r) x k` sparse matrix :math:`(k \geq 1)`,
% `W = nullSpaceOperatorApply(nullS, V)`;
% computes an `n x k` sparse matrix `W` from `V` such that :math:`S W = 0`.
%
% This is an operator form of finding an `n x (n-r)` matrix `Z`
% such that :math:`S Z = 0` and then computing :math:`W = Z V`.
% The aim is to obtain `W` without forming `Z` explicitly.
%
% .. 16 May 2008: (MAS) First version of nullspaceLUSOLapply.m.
%    See nullspaceLUSOLtest.m for testing.
% .. 20 Jan 2015: Ronan Fleming - Updated to use Nick Henderson's 64 bit LUSOL interface

Cinv    = nullS.Cinv; % Column scales
L       = nullS.L;    % Strictly triangular L
p       = nullS.p;    % Column permutation for S
rankS   = nullS.rank; % rank(S)

[n ,n ] = size(L);
[mV,nV] = size(V);

B       = [sparse(rankS,nV)
           V        ];
Z       = (L')\B;
W       = Z;
W(p,:)  = Z;
W       = Cinv*W;
return

author: opencobra

repo: cobratoolbox


text:
function [L,filters] = gaussSmooth( I, sigmas, shape, radius )
% Applies Gaussian smoothing to a (multidimensional) image.
%
% Smooths the n-dimensional array I with a n-dimensional gaussian with
% standard deviations specified by sigmas.  This operation in linearly
% seperable and is implemented as such.
%
% USAGE
%  [L,filters] = gaussSmooth( I, sigmas, [shape], [radius] )
%
% INPUTS
%  I       - input image
%  sigmas  - either n dimensional or 1 dimensional vector of standard devs
%            if sigmas(n)<=.3 then does not smooth along that dimension
%  shape   - ['full'] shape flag 'valid', 'full', 'same', or 'smooth'
%  radius  - [2.25] radius in units of standard deviation
%
% OUTPUTS
%  L       - smoothed image
%  filters - actual filters used, cell array of length n
%
% EXAMPLE
%  load trees; I=ind2gray(X,map);
%  I2 = gaussSmooth( I, 1, 'same' );
%  figure(1); im(I); figure(2); im(I2);
%
% See also FILTERGAUSS
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

nd = ndims(I);  if(length(sigmas)==1); sigmas=repmat(sigmas,[1,nd]); end
if( nd > length(sigmas)); error('Incorrect # of simgas specified'); end
sigmas = sigmas(1:nd);

if( isa( I, 'uint8' ) ); I = double(I); end
if( nargin<3 || isempty(shape) ); shape='full'; end
if( nargin<4 || isempty(radius) ); radius=2.25; end

% create and apply 1D gaussian masks along each dimension
L = I;  filters = cell(1,nd);
for i=1:nd
  if (sigmas(i)>.3)
    r = ceil( sigmas(i)*radius );
    f = filterGauss( 2*r+1, [], sigmas(i)^2 );
    f = permute( f, circshift(1:nd,[1,i-1]) );
    filters{i} = f;
    L = convnFast( L, f, shape );
  else
    filters{i} = 1;
  end
end

author: zouchuhang

repo: LayoutNet


text:
function [pos, zp] = planePosition(point, plane)
%PLANEPOSITION Compute position of a point on a plane.
%
%   COORDS2D = planePosition(POINT, PLANE)
%   POINT has format [X Y Z], and plane has format
%   [X0 Y0 Z0  DX1 DY1 DZ1  DX2 DY2 DZ2], where :
%   - (X0, Y0, Z0) is a point belonging to the plane
%   - (DX1, DY1, DZ1) is a first direction vector
%   - (DX2, DY2, DZ2) is a second direction vector
%   Result COORDS2D has the form [XP YP], with XP and YP the coordinates of
%   the point in the coordinate system of the plane.
%
%   [COORDS2D, Z] = planePosition(POINT, PLANE)
%   Also returns the coordinates along the normal of the plane. When the
%   point is within the plane, this coordinate should be zero.
%
%   Example
%     plane = [10 20 30  1 0 0  0 1 0];
%     point = [13 24 35];
%     pos = planePosition(point, plane)
%     pos =
%         3   4
%
%   Example
%     % plane with non unit direction vectors
%     p0 = [30 20 10]; v1 = [2 1 0]; v2 = [-2 4 0];
%     plane = [p0 v1 v2];
%     pts = [p0 ; p0 + v1 ; p0 + v2 ; p0 + 3 * v1 + 2 * v2];
%     pos = planePosition(pts, plane)
%     pos =
%          0     0
%          1     0
%          0     1
%          3     2
%
%
%   See also
%     geom3d, planes3d, points3d, planePoint
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2005-02-21
% Copyright 2005-2022 INRA - TPV URPOI - BIA IMASTE

% size of input arguments
npl = size(plane, 1);
npt = size(point, 1);

% check inputs have compatible sizes
if npl ~= npt && npl > 1 && npt > 1
    error('geom3d:planePosition:inputSize', ...
        'plane and point should have same size, or one of them must have 1 row');
end

% origin and direction vectors of the plane(s)
p0 = plane(:, 1:3);
v1 = plane(:, 4:6);
v2 = plane(:, 7:9);

% Principle
%
% for each (recentered) point, we need to solve the system:
%   s * v1x + t * v2x + u * v3x  = px
%   s * v1y + t * v2y + u * v3y  = py
%   s * v1z + t * v2z + u * v3z  = pz
% Assuming the point belong to the place, the value of u is 0. The last
% equation is kept only for homogeneity.
% We rewrite in matrix form:
%   [ v1x v2x v3x ]   [ s ]   [ xp ]
%   [ v1y v2y v3y ] * [ t ] = [ yp ]
%   [ v1z v2z v3z ]   [ u ]   [ zp ]
% Or:
%   A * X = B
% We need to solve X = inv(A) * B. In practice, we use the b/A notation,
% obtained after using transpositon on the A\b notation.
%   X' = (inv(mat) * bsxfun(@minus, point, p0)')';
%   X' = bsxfun(@minus, point, p0) * inv(mat');
%   X' = bsxfun(@minus, point, p0) / mat';
%   (and we compute the transpose of the basis transform)

% Compute dot products with direction vectors of the plane
if npl == 1
    % we have npl == 1 and npt > 1
    % build transpose of matrix that changes plane basis to global basis
    mat = [v1 ; v2 ; cross(v1, v2, 2)];
    tmp = bsxfun(@minus, point, p0) / mat;
    pos = tmp(:, 1:2);
    zp = tmp(:,3);

else
    % NPL > 0 -> iterate over planes
    % Number of points can be either 1, or the same number as planes
    pos = zeros(npl, 2);
    zp = zeros(npl, 1);
    for ipl = 1:npl
        mat = [v1(ipl,:) ; v2(ipl,:) ; cross(v1(ipl,:), v2(ipl,:), 2)];
        % choose either point with same index, or the single point
        ind = min(ipl, npt);
        tmp = bsxfun(@minus, point(ind,:), p0(ipl,:)) / mat;
        pos(ipl,:) = tmp(1,1:2);
        zp(ipl) = tmp(1,3);
    end
end

% % old version (:
% s = dot(point-p0, d1, 2) ./ vectorNorm3d(d1);
% t = dot(point-p0, d2, 2) ./ vectorNorm3d(d2);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function walpha(varargin)

	query = reshape(transpose([strvcat(varargin{:}),[repmat('+',size(strvcat(varargin{:}),1)-1,1);' ']]),1,[]);

	url = ['''http://www.wolframalpha.com/input/?i=',query,''''];

	url(url==' ') = [];

	eval(['web ', url ,' -browser'])

end


text:
function [ vec ] = ellipseParametricFunction( a,b,theta )
%% About:
% Calculate the X and Y coordinates of the ellipse point corrsponding to
% angle (theta) from the parametric equation

%% inputs:
% theta: angle of point on the ellipse (from the parametric equation of the ellipse)
% a,b: the big and the small radious of the ellipse

%% Return value:
% The return value is a column vector, vec=[x;y;0] where:
% x: the x coordiante of the point on the ellipse correspodning to angle
% (theta) as described in the frame of the ellipse
% y: the y coordiante of the point on the ellipse  correspodning to angle
% (theta) as described in the frame of the ellipse

% Mohammad SAFEEA 10-Nov-2017

x=a*cos(theta);
y=b*sin(theta);
vec=[x;y;0];

end

repo: KST-Kuka-Sunrise-Toolbox


text:
% read a file containing ROBOTLASER1 in CARMEN logfile format
function laser=read_robotlaser(filename)

fid = fopen(filename, 'r');
laser = cell();

while true
  ln = fgetl(fid);
  if (ln == -1)
    break
  endif
  tokens = strsplit(ln, ' ', true);

  if (strcmp(tokens(1), "ROBOTLASER1") == 0)
    continue;
  endif

  num_tokens = str2double(tokens);

  currentReading = struct (
    "start_angle", 0,
    "angular_resolution", 0,
    "maximum_range", 0,
    "ranges", [],
    "pose", zeros(3,1),
    "laser_offset", zeros(3,1),
    "timestamp", 0
  );

  tk = 3;
  currentReading.start_angle = num_tokens(tk++);
  tk++; % skip FOV
  currentReading.angular_resolution = num_tokens(tk++);
  currentReading.maximum_range = num_tokens(tk++);
  tk += 2; % skip accuracy, remission_mode

  num_readings = int32(num_tokens(tk++));
  currentReading.ranges = num_tokens(tk:tk+num_readings-1);
  tk += num_readings;

  num_remissions = int32(num_tokens(tk++)); % skip reading the remission values

  tk += num_remissions;

  laser_pose = num_tokens(tk:tk+2);
  tk += 3;
  currentReading.pose = num_tokens(tk:tk+2);
  tk += 3;

  currentReading.laser_offset = t2v(inv(v2t(currentReading.pose)) * v2t(laser_pose));

  tk += 5; % skip tv, rv, forward, side, turn
  currentReading.timestamp = num_tokens(tk++);

  laser{end+1} = currentReading;

endwhile

laser = cell2mat(laser);

end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
classdef IpoptOptions < matlab.mixin.SetGet
    %IpoptOptions Summary of this class goes here
    %   Detailed explanation goes here

    properties
        %Convergence
        constr_viol_tol(1,1) double = 0.0001
        dual_inf_tol(1,1) double = 1
        max_cpu_time(1,1) double = 1E6 %seconds
        max_iter(1,1) double = 3000
        tol(1,1) double = 1E-8

        %Line Search
        alpha_for_y(1,1) IpOptAlphaForYEnum = IpOptAlphaForYEnum.Primal
        alpha_red_factor(1,1) double = 0.5
        corrector_type(1,1) IpOptCorrectorTypeEnum = IpOptCorrectorTypeEnum.None;
        nu_inc(1,1) double = 0.0001
        nu_init(1,1) double = 1E-6
        recalc_y(1,1) IpOptRecalcYEnum = IpOptRecalcYEnum.No

        %Parallel
        useParallel(1,1) IpoptUseParallelEnum = IpoptUseParallelEnum.DoNotUseParallel;
        numWorkers(1,1) double = feature('numCores');
    end

    methods
        function obj = IpoptOptions()

        end

        function ipoptDotOptions = getOptionsForOptimizer(obj)
            ipoptDotOptions.print_level = 5;
            ipoptDotOptions.hessian_approximation = 'limited-memory';
            ipoptDotOptions.honor_original_bounds = 'yes';
            ipoptDotOptions.print_timing_statistics = 'yes';
%             ipoptDotOptions.print_user_options = 'yes';
            ipoptDotOptions.mu_strategy = 'adaptive';
            ipoptDotOptions.linear_solver = 'mumps';

            %Convergence
            if(not(isnan(obj.constr_viol_tol)))
                ipoptDotOptions.constr_viol_tol = obj.constr_viol_tol;
            end

            if(not(isnan(obj.dual_inf_tol)))
                ipoptDotOptions.dual_inf_tol = obj.dual_inf_tol;
            end

            if(not(isnan(obj.max_cpu_time)))
                ipoptDotOptions.max_cpu_time = obj.max_cpu_time;
            end

            if(not(isnan(obj.max_iter)))
                ipoptDotOptions.max_iter = obj.max_iter;
            end

            if(not(isnan(obj.tol)))
                ipoptDotOptions.tol = obj.tol;
            end

            %Line Search
            ipoptDotOptions.alpha_for_y = obj.alpha_for_y.optStr;

            if(not(isnan(obj.alpha_red_factor)))
                ipoptDotOptions.alpha_red_factor = obj.alpha_red_factor;
            end

            ipoptDotOptions.corrector_type = obj.corrector_type.optStr;

            if(not(isnan(obj.nu_inc)))
                ipoptDotOptions.nu_inc = obj.nu_inc;
            end

%             if(not(isnan(obj.nu_init)))
%                 ipoptDotOptions.nu_init = obj.nu_init;
%             end
            obj.nu_init = 1000;

            ipoptDotOptions.recalc_y = obj.recalc_y.optStr;
        end

        function tf = usesParallel(obj)
            tf = obj.useParallel;
%             tf = IpoptUseParallelEnum.UseParallel;
        end

        function numWorkers = getNumParaWorkers(obj)
            numWorkers = obj.numWorkers;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot