text:
function moviefunction(cfg, varargin)

% we need cfg.plotfun to plot the data
% data needs to be 3D, N x time x freq (last can be singleton)
%   N needs to correspond to number of vertices (channels, gridpoints, etc)

% new UI artwork
%
% [main window] -------------------------------------\
% | [uipanel: plot]            | [uipanel: colormap] |
% |  <80%> ^70%v               |   <20%>             |
% |  [axes: mainaxes]          |  [axes: coloraxes]  |
% |   [axes: subaxes]          |                     |
% |                            |                     |
% |                            |                     |
% |                            |                     |
% |                            |                     |
% |                            |                     |
% |                            |                     |
% |                            |                     |
% |                            |                     |
% |                            |                     |
% |                            |                     |
% |                            |[Checkbox: automatic]|
% |                            |[Checkbox: symmetric]|
% |--[uipanel: control] <100%>-----------------------|
% | [Button:Go] [Checkbox:Rec]                       |
% | [Slider: frequency]                              |
% | [Slider: time]                                   |
% | [TextField: speed]                               |
% \--------------------------------------------------/
%
% Right click on mainaxes:
%     rotation options (+automatic), zooming options, open subplot (new subaxes added)
%     start, stop, record, time, freq, colorbar
% Left click on subaxes:  dragging it around
% Right click on subaxes: close
%
% Consider following neat extras
%   new figure for zooming and and rotation timing options during playback
%   add an intro a la Jan-Mathijs for source, similar for normal topos
%   make panels foldable (like on mobile devices)

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamedval', {'zlim',         'absmax',  'maxabs'});
cfg = ft_checkconfig(cfg, 'renamed',    {'zparam',       'funparameter'});
cfg = ft_checkconfig(cfg, 'renamed',	  {'parameter',    'funparameter'});
cfg = ft_checkconfig(cfg, 'renamed',	  {'mask',         'maskparameter'});
cfg = ft_checkconfig(cfg, 'renamed',	  {'framespersec', 'framerate'});
cfg = ft_checkconfig(cfg, 'deprecated', {'xparam'});
cfg = ft_checkconfig(cfg, 'forbidden',  {'yparam'});

% set default cfg settings
cfg.xlim            = ft_getopt(cfg, 'xlim',   'maxmin');
cfg.ylim            = ft_getopt(cfg, 'ylim',   'maxmin');
cfg.zlim            = ft_getopt(cfg, 'zlim',   'maxmin');
cfg.xparam          = ft_getopt(cfg, 'xparam', 'time');
cfg.yparam          = ft_getopt(cfg, 'yparam', []);
cfg.maskparameter   = ft_getopt(cfg, 'maskparameter');
cfg.inputfile       = ft_getopt(cfg, 'inputfile',    []);
cfg.moviefreq       = ft_getopt(cfg, 'moviefreq',    []);
cfg.movietime       = ft_getopt(cfg, 'movietime',    []);
cfg.movierpt        = ft_getopt(cfg, 'movierpt',     1);

for i=1:numel(varargin)
  if ~isempty(cfg.yparam) && ~isfield(varargin{i}, 'freq')
    ft_error('data argument %i does not have a .freq field, while all former had', i);
  elseif isempty(cfg.yparam) && isfield(varargin{i}, 'freq')
    cfg.yparam        = 'freq';
  else
    cfg.yparam        = [];
  end
end

% set guidata flags
opt                 = [];
opt.valx            = 1;
opt.valy            = 1;
opt.valz            = 1;
opt.record          = ~istrue(ft_getopt(cfg, 'interactive', 'yes'));
opt.framesfile      = ft_getopt(cfg, 'framesfile',   []);
opt.samperframe     = ft_getopt(cfg, 'samperframe',  1);
opt.framerate       = ft_getopt(cfg, 'framerate',    5);
opt.fixedframesfile = ~isempty(opt.framesfile);
opt.xvalues   = varargin{1}.(cfg.xparam); % below consistency is checked
if ~isempty(cfg.yparam)
  opt.yvalues = varargin{i}.(cfg.yparam); % below consistency is checked
else
  opt.yvalues   = [];
end

% check data options and consistency
Ndata = numel(varargin);
for i = 1:Ndata
  opt.issource(i)   = ft_datatype(varargin{i}, 'source');
  opt.isfreq(i)     = ft_datatype(varargin{i}, 'freq');
  opt.istimelock(i) = ft_datatype(varargin{i}, 'timelock');
  if opt.issource(i)+opt.isfreq(i)+opt.istimelock(i) ~= 1
    ft_error('data argument %i cannot be definitely identified as frequency, timelock or source data', i);
  end
end

if Ndata>1
  if all(opt.issource==1)
    % this is allowed maybe in the future: multiple source input arguments
    ft_error('currently, not more than 1 data argument is allowed, unless one of them is a parcellation, and the other a channel level structure');
  elseif all(opt.isfreq==1)
    % this is allowed maybe in the future: multiple freq input arguments
    ft_error('currently, not more than 1 data argument is allowed, unless one of them is a parcellation, and the other a channel level structure');
  elseif all(opt.istimelock==1)
    % this is allowed maybe in the future: multiple timelock input arguments
    ft_error('currently, not more than 1 data argument is allowed, unless one of them is a parcellation, and the other a channel level structure');
  elseif Ndata==2 && sum(opt.issource==1)
    % this may be allowed, provided the source data is a parcellation and
    % the other argument a parcellated data structure
    if ~ft_datatype(varargin{opt.issource}, 'parcellation')
      ft_error('the source data structure should be a parcellation');
    end
  end
end

% set the funparameter
if any(opt.isfreq)
  opt.funparameter = ft_getopt(cfg, 'funparameter', 'powspctrm'); % use power as default
elseif any(opt.istimelock)
  opt.funparameter = ft_getopt(cfg, 'funparameter', 'avg'); % use power as default
elseif any(opt.issource)
  % FIXME a call to ft_checkdata wouldn't hurt here :-)
  opt.funparameter = ft_getopt(cfg, 'funparameter', 'avg.pow'); % use power as default
end

if any(opt.issource) && (any(opt.isfreq) || any(opt.istimelock))
  opt.anatomy  = varargin{opt.issource};
  opt.ismesh   = isfield(opt.anatomy, 'tri');
  opt.isvolume = isfield(opt.anatomy, 'dim');

  varargin     = varargin(~opt.issource);
  opt.dat{1}   = getsubfield(varargin{1}, opt.funparameter);
  funtok = tokenize(opt.funparameter, '.');
  opt.dimord   = getdimord(varargin{1}, funtok{end});
elseif any(opt.issource)
  opt.anatomy  = varargin{opt.issource};
  opt.ismesh   = isfield(opt.anatomy, 'tri');
  opt.isvolume = isfield(opt.anatomy, 'dim');

  opt.dat{1}   = getsubfield(varargin{opt.issource}, opt.funparameter);
  funtok = tokenize(opt.funparameter, '.');
  opt.dimord   = getdimord(varargin{1}, funtok{end});
else
  opt.layout   = ft_prepare_layout(cfg); % let ft_prepare_layout do the error handling for now

  % ensure that the data in all inputs has the same channels, time-axis, etc.
  tmpcfg = keepfields(cfg, {'frequency', 'latency', 'channel'});
  [varargin{:}] = ft_selectdata(tmpcfg, varargin{:});
  % restore the provenance information
  [cfg, varargin{:}] = rollback_provenance(cfg, varargin{:});

  [opt.seldat, opt.sellay] = match_str(varargin{1}.label, opt.layout.label);
  opt.chanx = opt.layout.pos(opt.sellay,1);
  opt.chany = opt.layout.pos(opt.sellay,2);

  for i = 1:numel(varargin)
    opt.dat{i} = getsubfield(varargin{i}, opt.funparameter);
  end
  opt.dimord = getdimord(varargin{1}, opt.funparameter);
end

dimtok = tokenize(opt.dimord, '_');
if any(strcmp(dimtok, 'rpt') | strcmp(dimtok, 'subj'))
  ft_error('the input data cannot contain trials or subjects, please average first using ft_selectdata');
end

opt.ydim = find(strcmp(cfg.yparam, dimtok));
opt.xdim = find(strcmp(cfg.xparam, dimtok));
opt.zdim = setdiff(1:ndims(opt.dat{1}), [opt.ydim opt.xdim]);
if opt.zdim ~=1
  ft_error('input %i data does not have the correct format of N x time (x freq)', i);
end

% permute the data matrix
for i = 1:numel(opt.dat)
  opt.dat{i} = permute(opt.dat{i}, [opt.zdim(:)' opt.xdim opt.ydim]);
end
opt.xdim = 2;
if ~isempty(cfg.yparam)
  opt.ydim = 3;
else
  opt.ydim = [];
end

% get the mask
for i = 1:numel(varargin)
  if ~isempty(cfg.maskparameter) && ischar(cfg.maskparameter)
    opt.mask{i} = double(getsubfield(varargin{i}, cfg.maskparameter)~=0);
  else
    opt.mask{i} = ones(size(opt.dat{i}));
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% draw the figure with the GUI
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
opt = createGUI(opt);

if isempty(cfg.yparam)
%  set(opt.handles.label.yparam, 'Visible', 'off');
%  set(opt.handles.slider.yparam, 'Visible', 'off');
%  opt.yparam = '';
else
  opt.yparam = [upper(cfg.yparam(1)) lower(cfg.yparam(2:end))];
  set(opt.handles.label.yparam, 'String', opt.yparam);
end
opt.xparam = [upper(cfg.xparam(1)) lower(cfg.xparam(2:end))];
set(opt.handles.label.xparam, 'String', opt.xparam);

% set timer
opt.timer = timer;
set(opt.timer, 'timerfcn', {@cb_timer, opt.handles.figure}, 'period', 0.1, 'executionmode', 'fixedSpacing');

opt = plot_geometry(opt);
opt = plot_other(opt);

opt.speed = 1;
guidata(opt.handles.figure, opt);
% init first screen

cb_slider(opt.handles.figure);
%uicontrol(opt.handles.colorbar);
if strcmp(cfg.zlim, 'maxmin')
  set(opt.handles.checkbox.automatic, 'Value', 1)
  set(opt.handles.checkbox.symmetric, 'Value', 0)
  cb_colorbar(opt.handles.figure);
elseif strcmp(cfg.zlim, 'maxabs')
  set(opt.handles.checkbox.automatic, 'Value', 1)
  set(opt.handles.checkbox.symmetric, 'Value', 1)
  cb_colorbar(opt.handles.figure);
else
  set(opt.handles.checkbox.automatic, 'Value', 0)
  set(opt.handles.checkbox.symmetric, 'Value', 0)
  cb_colorbar(opt.handles.figure, cfg.zlim);
end

if opt.record
  opt.record = false;
  opt.quit = true;
  guidata(opt.handles.figure, opt);
  cb_recordbutton(opt.handles.button.record);
else
  opt.quit = false;
  guidata(opt.handles.figure, opt);
end

end

%%
%  **************************************************************
%  ********************* CREATE GUI *****************************
%  **************************************************************
function opt = createGUI(opt)

%% main figure
opt.handles.figure = figure(    ...
  'Units',       'normalized',  ...
  'Name',        'ft_movieplot',...
  'Menu',        'none',        ...
  'Toolbar',     'figure',      ...
  'NumberTitle', 'off',         ...
  'UserData',    [],            ...
  'Tag',         'mainFigure',  ...
  'Visible',     'on',          ...
  'windowbuttondownfcn', @cb_getposition);
  %'WindowButtonUpFcn', @cb_stopDrag, ...

%%  panels
clf

% visualization panel for the geometrical information
opt.handles.panel.visualization_geometry = uipanel(...
  'tag',      'mainPanels1',      ... % tag according to position
  'parent',   opt.handles.figure, ...
  'units',    'normalized',       ...
  'title',    'Visualization_geometry', ...
  'clipping', 'on',               ...
  'visible',  'on');

% rearrange
ft_uilayout(opt.handles.figure,     ...
  'tag',             'mainPanels1', ...
  'backgroundcolor', [.8 .8 .8],    ...
  'hpos',            'auto',        ...
  'vpos',            .15,           ...
  'halign',          'left',        ...
  'width',           1,             ...
  'height',          0.1);

% visualization panel for the non-geometrical information
opt.handles.panel.visualization = uipanel(...
  'tag',      'mainPanels2',      ... % tag according to position
  'parent',   opt.handles.figure, ...
  'units',    'normalized',       ...
  'title',    'Visualization',    ...
  'clipping', 'on',               ...
  'visible',  'on');

% rearrange
ft_uilayout(opt.handles.figure,     ...
  'tag',             'mainPanels2', ...
  'backgroundcolor', [.8 .8 .8],    ...
  'hpos',            'auto',        ...
  'vpos',            .55,           ...
  'halign',          'left',        ...
  'width',           0.5,           ...
  'height',          0.1);

% settings panel (between different views can be switched)
% opt.handles.panel.view = uipanel(...
%   'tag',    'sidePanels', ... % tag according to position
%   'parent', opt.handles.figure,...
%   'units',  'normalized', ...
%   'title',  'View',...
%   'clipping','on',...
%   'visible', 'on');

% settings panel ()
opt.handles.panel.settings = uipanel(...
  'tag',      'sidePanels',       ... % tag according to position
  'parent',   opt.handles.figure, ...
  'units',    'normalized',       ...
  'title',    'Settings',         ...
  'clipping', 'on',               ...
  'visible',  'on');

% rearrange panel
ft_uilayout(opt.handles.figure,    ...
  'tag',             'sidePanels', ...
  'backgroundcolor', [.8 .8 .8],   ...
  'hpos',            'auto',       ...
  'vpos',            0.15,         ...
  'halign',          'right',      ...
  'width',           0.15,         ...
  'height',          0.55);

% control panel
opt.handles.panel.controls = uipanel(...
  'tag',    'lowerPanels', ... % tag according to position
  'parent', opt.handles.figure,...
  'units',  'normalized', ...
  'title',  'Controls',...
  'clipping','on',...
  'visible', 'on');

% rearrange
ft_uilayout(opt.handles.figure,     ...
  'tag',             'lowerPanels', ...
  'backgroundcolor', [0.8 0.8 0.8], ...
  'hpos',            'auto',        ...
  'vpos',            0,             ...
  'halign',          'right',       ...
  'width',           1,             ...
  'height',          0.15);

ft_uilayout(opt.handles.figure, ...
  'tag',   'mainPanels1',       ...
  'retag', 'sidePanels');

ft_uilayout(opt.handles.figure, ...
  'tag',   'mainPanels2',       ...
  'retag', 'sidePanels');

ft_uilayout(opt.handles.figure,    ...
  'tag',             'sidePanels', ...
  'backgroundcolor', [.8 .8 .8],   ...
  'hpos',            'auto',       ...
  'vpos',            'align',      ...
  'halign',          'left',       ...
  'valign',          'top',        ...
  'height',          .85);

% add axes
% 3 axes for switching to topo viewmode
% opt.handles.axes.A = axes(...
%   'Parent',opt.handles.panel.view,...
%   'Position',[0 0 1 .33], ...
%   'Tag','topoAxes',...
%   'HandleVisibility', 'on', ...
%   'UserData',[]);
%
% opt.handles.axes.B = axes(...
%   'Parent',opt.handles.panel.view,...
%   'Position',[0 0.33 1 .33], ...
%   'Tag','topoAxes',...
%   'HandleVisibility', 'on', ...
%   'UserData',[]);
%
% opt.handles.axes.C = axes(...
%   'Parent',opt.handles.panel.view,...
%   'Position',[0 0.66 1 .33], ...
%   'Tag','topoAxes',...
%   'HandleVisibility', 'on', ...
%   'UserData',[]);

% control panel uicontrols

opt.handles.button.play = uicontrol(...
  'parent',opt.handles.panel.controls,...
  'tag', 'controlInput', ...
  'string','Play',...
  'units', 'normalized', ...
  'style', 'pushbutton', ...
  'callback',@cb_playbutton,...
  'userdata', 'space');

opt.handles.button.record = uicontrol(...
  'parent',opt.handles.panel.controls,...
  'tag', 'controlInput', ...
  'string','Record?',...
  'units', 'normalized', ...
  'style', 'checkbox', ...
  'Callback',@cb_recordbutton,...
  'userdata', 'enter');

ft_uilayout(opt.handles.panel.controls, ...
  'tag', 'controlInput', ...
  'width', 0.15, ...
  'height', 0.25, ...
  'hpos', 'auto', ...
  'vpos', .75, ...
  'halign', 'right', ...
  'valign', 'top');

opt.handles.text.speed = uicontrol(...
  'parent',opt.handles.panel.controls,...
  'tag', 'speedInput', ...
  'string','Speed',...
  'units', 'normalized', ...
  'style', 'text');

opt.handles.edit.speed = uicontrol(...
  'parent',opt.handles.panel.controls,...
  'tag', 'speedInput', ...
  'string','1.0',...
  'units', 'normalized', ...
  'Callback',@cb_speed, ...
  'style', 'edit');

ft_uilayout(opt.handles.panel.controls, ...
  'tag', 'speedInput', ...
  'width', 0.15, ...
  'height', 0.25, ...
  'hpos', 'auto', ...
  'vpos', .25, ...
  'halign', 'right', ...
  'valign', 'top');

% opt.handles.button.faster = uicontrol(...
%   'parent',opt.handles.panel.controls,...
%   'tag', 'controlInput', ...
%   'string','+',...
%   'units', 'normalized', ...
%   'style', 'pushbutton', ...
%   'userdata', '+', ...
%   'Callback',@cb_slider);
%
%
% opt.handles.button.slower = uicontrol(...
%   'parent',opt.handles.panel.controls,...
%   'tag', 'controlInput', ...
%   'string','-',...
%   'units', 'normalized', ...
%   'style', 'pushbutton', ...
%   'userdata', '-', ...
%   'Callback',@cb_slider);

% ft_uilayout(opt.handles.panel.controls, ...
%   'tag', 'controlInput', ...
%   'width', 0.15, ...
%   'height', 0.25, ...
%   'hpos', 'auto', ...
%   'vpos', .75, ...
%   'valign', 'top');
%
% set(opt.handles.button.slower, 'tag', 'speedButtons');
% set(opt.handles.button.faster, 'tag', 'speedButtons');
%
%
% ft_uilayout(opt.handles.panel.controls, ...
%   'tag', 'speedButtons', ...
%   'width', 0.05, ...
%   'height', 0.125, ...
%   'vpos', 'auto');
%
% ft_uilayout(opt.handles.panel.controls, ...
%   'tag', 'speedButtons', ...
%   'hpos', 'align');
%
%
% ft_uilayout(opt.handles.panel.controls, ...
%   'tag', 'speedButtons', ...
%   'retag', 'controlInput');

% speed control

opt.handles.slider.xparam = uicontrol(...
  'Parent',opt.handles.panel.controls,...
  'Units','normalized',...
  'BackgroundColor',[0.9 0.9 0.9],...
  'Callback',@cb_slider,...
  'Position',[0.0 0.75 0.625 0.25],...
  'SliderStep', [1/numel(opt.xvalues) 1/numel(opt.xvalues)],...
  'String',{  'Slider' },...
  'Style','slider',...
  'Tag','xparamslider');

opt.handles.label.xparam = uicontrol(...
  'Parent',opt.handles.panel.controls,...
  'Units','normalized',...
  'Position',[0.0 0.5 0.625 0.25],...
  'String','xparamLabel',...
  'Style','text',...
  'Tag','xparamLabel');

if ~isempty(opt.ydim)
  opt.handles.slider.yparam = uicontrol(...
    'Parent',opt.handles.panel.controls,...
    'Units','normalized',...
    'BackgroundColor',[0.9 0.9 0.9],...
    'Callback',@cb_slider,...
    'Position',[0.0 0.25 0.625 0.25],...
    'SliderStep', [1/numel(opt.yvalues) 1/numel(opt.yvalues)],...
    'String',{  'Slider' },...
    'Style','slider',...
    'Tag','yparamSlider');

  opt.handles.label.yparam = uicontrol(...
    'Parent',opt.handles.panel.controls,...
    'Units','normalized',...
    'Position',[0.0 0.0 0.625 0.25],...
    'String','yparamLabel',...
    'Style','text',...
    'Tag','text3');
end

opt.handles.axes.colorbar = axes(...
  'Parent',opt.handles.panel.settings,...
  'Units','normalized',...
  'Position',[0 0 1.0 1],...
  'ButtonDownFcn',@cb_color,...
  'HandleVisibility', 'on', ...
  'Tag','colorbarAxes', ...
  'YLim', [0 1], ...
  'YLimMode', 'manual', ...
  'XLim', [0 1], ...
  'XLimMode', 'manual');

opt.handles.colorbar = colorbar('ButtonDownFcn', []);
set(opt.handles.colorbar, 'Parent', opt.handles.panel.settings);
set(opt.handles.colorbar, 'Position', [0.4 0.2 0.2 0.65]);
nColors = size(colormap, 1);
set(opt.handles.colorbar, 'YTick', [1 nColors/4 nColors/2 3*nColors/4 nColors]);
if (gcf~=opt.handles.figure)
  close gcf; % sometimes there is a new window that opens up
end
%
% set lines
%YLim = get(opt.handles.colorbar, 'YLim');
opt.handles.lines.upperColor = line([-1 0], [32 32], ...
  'Color', 'black', ...
  'LineWidth', 4, ...
  'ButtonDownFcn', @cb_startDrag, ...
  'Parent', opt.handles.colorbar, ...
  'Visible', 'off', ...
  'Tag', 'upperColor');

opt.handles.lines.lowerColor = line([1 2], [32 32], ...
  'Color', 'black', ...
  'LineWidth', 4, ...
  'ButtonDownFcn', @cb_startDrag, ...
  'Parent', opt.handles.colorbar, ...
  'Visible', 'off', ...
  'Tag', 'lowerColor');

opt.handles.menu.colormap = uicontrol(...
  'Parent',opt.handles.panel.settings,...
  'Units','normalized',...
  'BackgroundColor',[1 1 1],...
  'Callback',@cb_colormap,...
  'Position',[0.25 0.95 0.5 0.05],...
  'String',{  'jet'; 'hot'; 'cool'; 'ikelvin'; 'ikelvinr' },...
  'Style','popupmenu',...
  'Value',1, ...
  'Tag','colormapMenu');

opt.handles.checkbox.automatic = uicontrol(...
  'Parent',opt.handles.panel.settings,...
  'Units','normalized',...
  'Callback',@cb_colorbar,...
  'Position',[0.10 0.05 0.75 0.05],...
  'String','automatic',...
  'Style','checkbox',...
  'Value', 1, ... % TODO make this dependet on cfg.zlim
  'Tag','autoCheck');

opt.handles.checkbox.symmetric = uicontrol(...
  'Parent',opt.handles.panel.settings,...
  'Units','normalized',...
  'Callback',@cb_colorbar,...
  'Position',[0.10 0.0 0.75 0.05],...
  'String','symmetric',...
  'Style','checkbox',...
  'Value', 0, ... % TODO make this dependet on cfg.zlim
  'Tag','symCheck');

% buttons

opt.handles.button.decrLowerColor = uicontrol(...
  'Parent', opt.handles.panel.settings,...
  'Units','normalized',...
  'Enable', 'off', ...
  'Callback',@cb_colorbar,...
  'Position',[0.35 0.125 0.15 0.05],...
  'String','-',...
  'Tag','decrLowerColor');

opt.handles.button.incrLowerColor = uicontrol(...
  'Parent', opt.handles.panel.settings,...
  'Units','normalized',...
  'Enable', 'off', ...
  'Callback',@cb_colorbar,...
  'Position',[0.5 0.125 0.15 0.05],...
  'String','+',...
  'Tag','incrLowerColor');

opt.handles.button.decrUpperColor = uicontrol(...
  'Parent', opt.handles.panel.settings,...
  'Units','normalized',...
  'Enable', 'off', ...
  'Callback',@cb_colorbar,...
  'Position',[0.35 0.875 0.15 0.05],...
  'String','-',...
  'Tag','decrUpperColor');

opt.handles.button.incrUpperColor = uicontrol(...
  'Parent', opt.handles.panel.settings,...
  'Units','normalized',...
  'Enable', 'off', ...
  'Callback',@cb_colorbar,...
  'Position',[0.5 0.875 0.15 0.05],...
  'String','+',...
  'Tag','incrUpperColor');

% Handle to the axes that will contain the geometry
opt.handles.axes.movie = axes(...
  'Parent',             opt.handles.panel.visualization_geometry,...
  'Position',           [0 0 1 1],                     ...
  'CameraPosition',     [0.5 0.5 9.16025403784439],    ...
  'CameraPositionMode', get(0,'defaultaxesCameraPositionMode'),...
  'CLim',               get(0,'defaultaxesCLim'),      ...
  'CLimMode',           'manual',                      ...
  'Color',              [0.9 0.9 0.94],                ...
  'ColorOrder',         get(0,'defaultaxesColorOrder'),...
  'XColor',             get(0,'defaultaxesXColor'),    ...
  'YColor',             get(0,'defaultaxesYColor'),    ...
  'ZColor',             get(0,'defaultaxesZColor'),    ...
  'HandleVisibility',   'on',                          ...
  'ButtonDownFcn',      @cb_view,                      ...
  'Tag',                'geometry');

% Handle to the axes that will contain the non-geometry
opt.handles.axes.other = axes(...
  'Parent',             opt.handles.panel.visualization,...
  'Position',           [0 0 1 1],                     ...
  'CameraPosition',     [0.5 0.5 9.16025403784439],    ...
  'CameraPositionMode', get(0,'defaultaxesCameraPositionMode'),...
  'CLim',               get(0,'defaultaxesCLim'),      ...
  'CLimMode',           'manual',                      ...
  'Color',              [0.9 0.9 0.94],                ...
  'ColorOrder',         get(0,'defaultaxesColorOrder'),...
  'XColor',             get(0,'defaultaxesXColor'),    ...
  'YColor',             get(0,'defaultaxesYColor'),    ...
  'ZColor',             get(0,'defaultaxesZColor'),    ...
  'HandleVisibility',   'on',                          ...
  'ButtonDownFcn',      @cb_view,                      ...
  'Tag',                'other');

% Disable axis labels
axis(opt.handles.axes.movie,    'equal');
axis(opt.handles.axes.colorbar, 'equal');
% axis(opt.handles.axes.A, 'equal');
% axis(opt.handles.axes.B, 'equal');
% axis(opt.handles.axes.C, 'equal');
%
axis(opt.handles.axes.movie,    'off');
axis(opt.handles.axes.colorbar, 'off');
% axis(opt.handles.axes.A, 'off');
% axis(opt.handles.axes.B, 'off');
% axis(opt.handles.axes.C, 'off');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function update_panels(opt)
  for i=1:numel(opt.dat)
    if ~any(opt.issource)
      set(opt.handles.grid{i}, 'cdata', griddata(opt.chanx{i}, opt.chany{i}, opt.mask{i}(:,opt.valx,opt.valy).*opt.dat{i}(:,opt.valx,opt.valy), opt.xdata{i}, opt.nanmask{i}.*opt.ydata{i}, 'v4'));
    else
      set(opt.handles.mesh{i}, 'FaceVertexCData',     squeeze(opt.dat{i}(:,opt.valx,opt.valy)));
      set(opt.handles.mesh{i}, 'FaceVertexAlphaData', squeeze(opt.mask{i}(:,opt.valx,opt.valy)));
    end
  end

  if opt.doplot
    opt.valz
    set(get(opt.handles.axes.other,'children'), 'ydata', opt.dat{1}(opt.valz,:));
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALLBACK FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_speed(h, eventdata)
if ~ishandle(h)
  return
end
opt = guidata(h);

val = get(h, 'String');
% val = get(opt.handles.slider.speed, 'value');
% val = exp(log(opt.MAX_SPEED) * (val-.5)./0.5);
%
% % make sure we can easily get back to normal speed
% if abs(val-opt.AVG_SPEED) < 0.08
%   val = opt.AVG_SPEED;
%   set(opt.handles.slider.speed, 'value', 0.5);
% end

speed = str2num(val);
if isempty(speed)
  speed = opt.speed;
end
opt.speed = speed;

set(h, 'String', opt.speed)

% if val >=100
%   set(opt.handles.label.speed, 'String', ['Speed ' num2str(opt.speed, '%.1f'), 'x'])
% elseif val >= 10
%   set(opt.handles.label.speed, 'String', ['Speed ' num2str(opt.speed, '%.2f'), 'x'])
% else
%   set(opt.handles.label.speed, 'String', ['Speed ' num2str(opt.speed, '%.3f'), 'x'])
% end

guidata(h, opt);
uiresume;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALLBACK FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_timer(obj, info, h)
if ~ishandle(h)
  return
end
opt = guidata(h);
delta = opt.speed/numel(opt.xvalues);
val = get(opt.handles.slider.xparam, 'value');
val = val + delta;

if opt.record
  if val>1
    % stop recording
    stop(opt.timer);
    % reset again
    val = 0;
    set(opt.handles.slider.xparam, 'value', val);
    cb_slider(h);
    cb_recordbutton(opt.handles.button.record);
    % TODO FIXME add some message here
    guidata(h, opt);
    return;
  end
end

while val>1
    val = val-1;
end
set(opt.handles.slider.xparam, 'value', val);
cb_slider(h);

if opt.record
  pause(.1);
  drawnow;
  vs = version('-release');
  vs = vs(1:4);
  % get starting position via parameter panel
  currFrame = getframe(opt.handles.figure);
  for i=1:opt.samperframe
    if vs<2010
      opt.vidObj = addframe(opt.vidObj, currFrame);
    else
      writeVideo(opt.vidObj, currFrame);
    end

  end
  guidata(h, opt);
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALLBACK FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_slider(h, eventdata)
opt = guidata(h);
valx = get(opt.handles.slider.xparam, 'value');
valx = round(valx*(numel(opt.xvalues)-1))+1;
valx = min(valx, numel(opt.xvalues));
valx = max(valx, 1);
set(opt.handles.label.xparam, 'String', [opt.xparam ' ' num2str(opt.xvalues(valx), '%.2f') 's']);

if ~isempty(opt.yvalues)
  valy = get(opt.handles.slider.yparam, 'value');
  valy = round(valy*(numel(opt.yvalues)-1))+1;
  valy = min(valy, numel(opt.yvalues));
  valy = max(valy, 1);

  if valy ~= opt.valy
    cb_colorbar(h);
  end

  set(opt.handles.label.yparam, 'String', [opt.yparam ' ' num2str(opt.yvalues(valy), '%.2f') 'Hz']);
else
  valy = 1;
end

opt.valx = valx;
opt.valy = valy;

guidata(h, opt);

update_panels(opt);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALLBACK FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_playbutton(h, eventdata)
if ~ishandle(h)
  return
end
opt = guidata(h);
switch get(h, 'string')
  case 'Play'
    set(h, 'string', 'Stop');
    start(opt.timer);
  case 'Stop'
    set(h, 'string', 'Play');
    stop(opt.timer);
end
uiresume;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALLBACK FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_recordbutton(h, eventdata)
if ~ishandle(h)
  return
end

opt = guidata(h);

opt.record = ~opt.record; % switch state
guidata(h, opt);

if opt.record

  if ~opt.fixedframesfile
    % open a save-file dialog
    [FileName,PathName,FilterIndex] = uiputfile('*.avi', 'Save AVI-file' , 'ft_movie');

    if (FileName == 0 && PathName == 0) % aborted
      cb_recordbutton(h);
      return;
    end

    opt.framesfile = fullfile(PathName, FileName);

    if (FilterIndex==1) % remove .avi again (4 chars)
      opt.framesfile = opt.framesfile(1:end-4);
    end

  end

  % FIXME open new window to play in there, so that frame getting works
  vs = version('-release');
  vs = vs(1:4);
  if vs<2010
    opt.vidObj = avifile(opt.framesfile, 'FPS', opt.framerate);
  else
    opt.vidObj = VideoWriter(opt.framesfile, 'Uncompressed AVI');
    opt.vidObj.FrameRate = opt.framerate;
    open(opt.vidObj);
  end

  %set(opt.handles.figure,'renderer','opengl')
  %opengl software;
  set(opt.handles.figure,'renderer','zbuffer');
  %opt.vidObj = avifile(opt.framesfile, 'fps', opt.framerate, 'quality', 75);

  set(h, 'string', 'Stop');
  guidata(h, opt);
  start(opt.timer);
else
  % FIXME set handle back to old window
  stop(opt.timer);
  if ~isempty(opt.framesfile)
    vs = version('-release');
    vs = vs(1:4);
    if vs<2010
      opt.vidObj = close(opt.vidObj);
    else
      close(opt.vidObj);
    end
  end
  set(h, 'string', 'Record');
  guidata(h, opt);

  if (opt.quit)
    close(opt.handles.figure);
  end
end

%
% % This function should open a new window, plot in there, extract every
% % frame, store the movie in opt and return again
%
% % this is not needed, no new window is needed
% %scrsz = get(0, 'ScreenSize');
% %f = figure('Position',[1 1 scrsz(3) scrsz(4)]);
%
% % FIXME disable buttons (apart from RECORD) when recording
% % if record is pressed, stop recording and immediately return
%
% % adapted from ft_movieplotTFR
%
% % frequency/time selection
% if ~isempty(opt.yvalues) && any(~isnan(yvalues))
%   if ~isempty(cfg.movietime)
%     indx = cfg.movietime;
%     for iFrame = 1:floor(size(opt.dat, opt.xdim)/cfg.samperframe)
%       indy = ((iFrame-1)*cfg.samperframe+1):iFrame*cfg.samperframe;
%       updateMovie(opt, indx, indy);
%       F(iFrame) = getframe;
%     end
%   elseif ~isempty(cfg.moviefreq)
%     indy = cfg.moviefreq;
%     for iFrame = 1:floor(size(opt.dat, opt.ydim)/cfg.samperframe)
%       indx = ((iFrame-1)*cfg.samperframe+1):iFrame*cfg.samperframe;
%       updateMovie(opt, indx, indy);
%       F(iFrame) = getframe;
%     end
%   else
%     ft_error('Either moviefreq or movietime should contain a bin number')
%   end
% else
%   for iFrame = 1:floor(size(opt.dat, opt.xdim)/cfg.samperframe)
%     indx = ((iFrame-1)*cfg.samperframe+1):iFrame*cfg.samperframe;
%     updateMovie(opt, indx, 1);
%     F(iFrame) = getframe;
%   end
% end
%
% % save movie
% if ~isempty(cfg.framesfile)
%   save(cfg.framesfile, 'F');
% end
% % play movie
% movie(F, cfg.movierpt, cfg.framespersec);

% uiresume;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALLBACK FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_colormap(h, eventdata)
maps = get(h, 'String');
val = get(h, 'Value');

while ~strcmp(get(h, 'Tag'), 'mainFigure')
  h = get(h, 'Parent');
end

opt =  guidata(h);

cmap = feval(maps{val}, size(colormap, 1));
% if strcmp(maps{val}, 'ikelvin')
%   cmap = ikelvin(size(colormap, 1));
% elseif strcmp(maps{val}, 'kelvin')
%   cmap = kelvin(size(colormap, 1));
% else
%   cmap = ft_colormap(opt.handles.axes.movie, maps{val});
% end

if get(opt.handles.checkbox.automatic, 'Value')
  ft_colormap(opt.handles.axes.movie_subplot{1}, cmap);
end

adjust_colorbar(opt);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALLBACK FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_colorbar(h, eventdata)
if strcmp(get(h, 'Tag'), 'mainFigure') % this is the init call
  incr = false;
  decr = false;
else
  incr = strcmp(get(h, 'String'), '+');
  decr = strcmp(get(h, 'String'), '-');
  lower = strfind(get(h, 'Tag'), 'Lower')>0;
  while ~strcmp(get(h, 'Tag'), 'mainFigure')
    h = get(h, 'Parent');
  end
end

opt =  guidata(h);
if (~incr&&~decr&&exist('eventdata', 'var')&&~isempty(eventdata)) % init call
  caxis(opt.handles.axes.movie, eventdata);
end
zmin = inf;
zmax = -inf;
for i=1:numel(opt.dat)
  [tmpmin, tmpmax] = caxis(opt.handles.axes.movie_subplot{i});
  zmin = min(tmpmin, zmin);
  zmax = max(tmpmax, zmax);
end
yLim = get(opt.handles.colorbar, 'YLim');

if incr
  yTick = linspace(zmin, zmax, yLim(end));
  if get(opt.handles.checkbox.symmetric, 'Value')
    zmin = zmin - mean(diff(yTick));
    zmax = zmax + mean(diff(yTick));
  elseif (lower)
    zmin = zmin + mean(diff(yTick));
  else
    zmax = zmax + mean(diff(yTick));
  end
elseif decr
  yTick = linspace(zmin, zmax, yLim(end));
  if get(opt.handles.checkbox.symmetric, 'Value')
    zmin = zmin + mean(diff(yTick));
    zmax = zmax - mean(diff(yTick));
  elseif (lower)
    zmin = zmin - mean(diff(yTick));
  else
    zmax = zmax - mean(diff(yTick));
  end
elseif get(opt.handles.checkbox.automatic, 'Value') % if automatic
  set(opt.handles.lines.upperColor, 'Visible', 'off');
  set(opt.handles.lines.lowerColor, 'Visible', 'off');
  set(opt.handles.lines.upperColor, 'YData', [yLim(end)/2 yLim(end)/2]);
  set(opt.handles.lines.lowerColor, 'YData', [yLim(end)/2 yLim(end)/2]);
  set(opt.handles.button.incrLowerColor, 'Enable', 'off');
  set(opt.handles.button.decrUpperColor, 'Enable', 'off');
  set(opt.handles.button.incrUpperColor, 'Enable', 'off');
  set(opt.handles.button.decrLowerColor, 'Enable', 'off');

  if get(opt.handles.checkbox.symmetric, 'Value') % maxabs
    zmax = -inf;
    for i=1:numel(opt.dat)
      tmpmax = max(max(abs(opt.dat{i}(:,:,opt.valy))));
      zmax = max(tmpmax, zmax);
    end
    zmin = -zmax;
  else   % maxmin
    zmax = -inf;
    zmin = inf;
    for i=1:numel(opt.dat)
      tmpmin = min(min(opt.dat{i}(:,:,opt.valy)));
      tmpmax = max(max(opt.dat{i}(:,:,opt.valy)));
      zmax = max(tmpmax, zmax);
      zmin = min(tmpmin, zmin);
    end

  end
else
  set(opt.handles.lines.upperColor, 'Visible', 'on');
  set(opt.handles.lines.lowerColor, 'Visible', 'on')
  if get(opt.handles.checkbox.symmetric, 'Value') % maxabs
    set(opt.handles.button.incrLowerColor, 'Enable', 'off');
    set(opt.handles.button.decrUpperColor, 'Enable', 'off');
    set(opt.handles.button.incrUpperColor, 'Enable', 'on');
    set(opt.handles.button.decrLowerColor, 'Enable', 'on');
    for i=1:numel(opt.dat)
      [tmpmin, tmpmax] = caxis(opt.handles.axes.movie_subplot{i});
      zmax = max(tmpmax, zmax);
    end
    zmin = -zmax;
  else
    set(opt.handles.button.incrLowerColor, 'Enable', 'on');
    set(opt.handles.button.decrUpperColor, 'Enable', 'on');
    set(opt.handles.button.incrUpperColor, 'Enable', 'on');
    set(opt.handles.button.decrLowerColor, 'Enable', 'on');
    for i=1:numel(opt.dat)
      [tmpmin, tmpmax] = caxis(opt.handles.axes.movie_subplot{i});
      zmin = min(tmpmin, zmin);
      zmax = max(tmpmax, zmax);
    end
  end
end % incr, decr, automatic, else

maps = get(opt.handles.menu.colormap, 'String');
cmap = feval(maps{get(opt.handles.menu.colormap, 'Value')}, size(colormap, 1));
for i=1:numel(opt.dat)
  ft_colormap(opt.handles.axes.movie_subplot{i}, cmap);
end

adjust_colorbar(opt);

if (gcf~=opt.handles.figure)
  close gcf; % sometimes there is a new window that opens up
end
for i=1:numel(opt.dat)
  caxis(opt.handles.axes.movie_subplot{i}, [zmin zmax]);
end

nColors = size(colormap, 1);
%yTick = (zmax-zmin)*(get(opt.handles.colorbar, 'YTick')/nColors)+zmin
yTick = (zmax-zmin)*[0 .25 .5 .75 1]+zmin;

%yTick = linspace(zmin, zmax, yLim(end));
% truncate intelligently/-ish
%yTick = get(opt.handles.colorbar, 'YTick')/nColors;
yTick = num2str(yTick', 5);
set(opt.handles.colorbar, 'YTickLabel', yTick, 'FontSize', 8);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALLBACK FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_startDrag(h, eventdata)
f = get(h, 'Parent');
while ~strcmp(get(f, 'Tag'), 'mainFigure')
  f = get(f, 'Parent');
end
opt =  guidata(f);
opt.handles.current.line = h;

if strfind(get(h, 'Tag'), 'Color')>0
  opt.handles.current.axes = opt.handles.colorbar;
  opt.handles.current.color = true;
else
  disp('Figure out if it works for xparam and yparam');
  keyboard
end

set(f, 'WindowButtonMotionFcn', @cb_dragLine);

guidata(h, opt);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_getposition(h, eventdata)

h   = findobj(h, 'tag', 'mainFigure');
opt = guidata(h);

pos = get(get(h, 'currentaxes'), 'currentpoint');
switch get(get(h, 'currentaxes'), 'tag'),
  case 'geometry'
    if opt.ismesh
      % get the intersection with the mesh
      [ipos, d] = intersect_line(opt.anatomy.pos, opt.anatomy.tri, pos(1,:), pos(2,:));
      [md, ix]  = min(abs(d));

      dpos     = opt.anatomy.pos - ipos(ix*ones(size(opt.anatomy.pos,1),1),:);
      opt.valz = nearest(sum(dpos.^2,2),0);

    elseif opt.isvolume
    else
    end

  case 'other'
  otherwise
end

% if strcmp(get(get(h, 'currentaxes'), 'tag'), 'timecourse')
%   % get the current point
%   pos = get(opt.hy, 'currentpoint');
%   set(opt.sliderx, 'value', nearest(opt.xparam, pos(1,1))./numel(opt.xparam));
%   if isfield(opt, 'hline')
%     set(opt.slidery, 'value', nearest(opt.yparam, pos(1,2))./numel(opt.yparam));
%   end
% elseif strcmp(get(get(h, 'currentaxes'), 'tag'), 'mesh')
%   % get the current point, which is defined as the intersection through the
%   % axis-box (in 3D)
%   pos       = get(opt.hx, 'currentpoint');
%
%   % get the intersection with the mesh
%   [ipos, d] = intersect_line(opt.pos, opt.tri, pos(1,:), pos(2,:));
%   [md, ix]  = min(abs(d));
%
%   dpos      = opt.pos - ipos(ix*ones(size(opt.pos,1),1),:);
%   opt.vindx = nearest(sum(dpos.^2,2),0);
%
%   if isfield(opt, 'parcellation')
%     opt.pindx = find(opt.parcellation(opt.vindx,:));
%     disp(opt.pindx);
%   end
% elseif strcmp(get(get(h, 'currentaxes'), 'tag'), 'mesh2')
%   % get the current point, which is defined as the intersection through the
%   % axis-box (in 3D)
%   pos       = get(opt.hz, 'currentpoint');
%
%   % get the intersection with the mesh
%   [ipos, d] = intersect_line(opt.pos, opt.tri, pos(1,:), pos(2,:));
%   [md, ix]  = min(abs(d));
%
%   dpos      = opt.pos - ipos(ix*ones(size(opt.pos,1),1),:);
%   opt.vindx = nearest(sum(dpos.^2,2),0);
%
%   if isfield(opt, 'parcellation')
%     opt.pindx2 = find(opt.parcellation(opt.vindx,:));
%     disp(opt.pindx2);
%   end
%
% end
guidata(h, opt);
cb_slider(h);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_dragLine(h, eventdata)
opt =  guidata(h);
pt = get(opt.handles.current.axes, 'CurrentPoint');
yLim = get(opt.handles.colorbar, 'YLim');

% upper (lower) bar must not below (above) lower (upper) bar
if ~(opt.handles.current.line == opt.handles.lines.upperColor && ...
    (any(pt(3)*[1 1]<get(opt.handles.lines.lowerColor, 'YData')) || ...
    yLim(end) <= pt(3))) ...
    && ~(opt.handles.current.line == opt.handles.lines.lowerColor && ...
    (any(pt(3)*[1 1]>get(opt.handles.lines.upperColor, 'YData')) || ...
    yLim(1) >= pt(3)))
  set(opt.handles.current.line, 'YData', pt(3)*[1 1]);
end

adjust_colorbar(opt);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_stopDrag(h, eventdata)
while ~strcmp(get(h, 'Tag'), 'mainFigure')
  h = get(h, 'Parent');
end
set(h, 'WindowButtonMotionFcn', '');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function adjust_colorbar(opt)
  % adjust colorbar
  upper = get(opt.handles.lines.upperColor, 'YData');
  lower = get(opt.handles.lines.lowerColor, 'YData');
  if any(round(upper)==0) || any(round(lower)==0)
    return;
  end
  maps = get(opt.handles.menu.colormap, 'String');
  cmap = feval(maps{get(opt.handles.menu.colormap, 'Value')}, size(colormap, 1));
  cmap(round(lower(1)):round(upper(1)), :) = repmat(cmap(round(lower(1)), :), 1+round(upper(1))-round(lower(1)), 1);
  for i=1:numel(opt.dat)
    ft_colormap(opt.handles.axes.movie_subplot{i}, cmap);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function opt = plot_geometry(opt)
  numArgs = numel(opt.dat);
  numRows = floor(sqrt(numArgs));
  numCols = ceil(sqrt(numArgs));
  for i=1:numArgs
    axes(opt.handles.axes.movie);
    opt.handles.axes.movie_subplot{i} = gca;
    if isfield(opt, 'anatomy') && opt.ismesh
      if isfield(opt.anatomy, 'sulc') && ~isempty(opt.anatomy.sulc)
        vdat = opt.anatomy.sulc;
        vdat(vdat>0.5) = 0.5;
        vdat(vdat<-0.5)= -0.5;
        vdat = vdat-min(vdat);
        vdat = 0.35.*(vdat./max(vdat))+0.3;
        vdat = repmat(vdat,[1 3]);
        mesh = ft_plot_mesh(opt.anatomy, 'edgecolor', 'none', 'vertexcolor', vdat);
      else
        mesh = ft_plot_mesh(opt.anatomy, 'edgecolor', 'none', 'facecolor', [0.5 0.5 0.5]);
      end
      lighting gouraud
      % set(mesh, 'Parent', opt.handles.axes.movie);
      % mesh = ft_plot_mesh(source, 'edgecolor', 'none', 'vertexcolor', 0*opt.dat(:,1,1), 'facealpha', 0*opt.mask(:,1,1));
      opt.handles.mesh{i} = ft_plot_mesh(opt.anatomy, 'edgecolor', 'none', 'vertexcolor', opt.dat{i}(:,1,1));
      set(opt.handles.mesh{i}, 'AlphaDataMapping', 'scaled');
      set(opt.handles.mesh{i}, 'FaceVertexAlphaData', opt.mask{i}(:,opt.valx,opt.valy));
      % TODO FIXME below does not work
      %set(opt.handles.mesh, 'FaceAlpha', 'flat');
      %set(opt.handles.mesh, 'EdgeAlpha', 'flat');

      lighting gouraud
      cam1 = camlight('left');
%       set(cam1, 'Parent', opt.handles.axes.movie);
      cam2 = camlight('right');
%       set(cam2, 'Parent', opt.handles.axes.movie);
%       set(opt.handles.mesh, 'Parent', opt.handles.axes.movie);
      %   cameratoolbar(opt.handles.figure, 'Show');
    else
      axes(opt.handles.axes.movie)
      [dum, opt.handles.grid{i}] = ft_plot_topo(opt.layout{i}.pos(opt.sellay,1), opt.layout{i}.pos(opt.sellay,2), zeros(numel(opt.sellay{i}),1), 'mask', opt.layout{i}.mask, 'outline', opt.layout{i}.outline, 'interpmethod', 'v4', 'interplim', 'mask', 'parent', opt.handles.axes.movie);
      %[dum, opt.handles.grid] = ft_plot_topo(layout.pos(sellay,1), layout.pos(sellay,2), zeros(numel(sellay),1), 'mask',layout.mask,  'outline', layout.outline, 'interpmethod', 'v4', 'interplim', 'mask', 'parent', opt.handles.axes.movie);
      % set(opt.handles.grid, 'Parent', opt.handles.axes.movie);
      opt.xdata{i}   = get(opt.handles.grid{i}, 'xdata');
      opt.ydata{i}   = get(opt.handles.grid{i}, 'ydata');
      opt.nanmask{i} = 1-get(opt.handles.grid{i}, 'cdata');
      if (gcf~=opt.handles.figure)
        close gcf; % sometimes there is a new window that opens up
      end
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function opt = plot_other(opt)
  dimord  = opt.dimord;

  switch dimord
    case {'pos_time' 'pos_freq' 'chan_time' 'chan_freq'}
      opt.doplot    = true;
      opt.doimagesc = false;
    case {'pos_freq_time' 'chan_freq_time' 'chan_chan_freq' 'chan_chan_time' 'pos_pos_freq' 'pos_pos_time'}
      opt.doplot    = false;
      opt.doimagesc = false;
    otherwise
  end

  if opt.doplot
    plot(opt.handles.axes.other, opt.xvalues, nanmean(opt.dat{1}(opt.valz,:),1));
  elseif opt.doimagesc
  end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c = ikelvin(m)
%  pos    hue   sat   value
cu = [
  0.0     1/2   0     1.0
  0.125   1/2   0.6   0.95
  0.375   2/3   1.0   0.8
  0.5     2/3   1.0   0.3
  ];

cl = cu;
cl(:, 3:4) = cl(end:-1:1, 3:4);
cl(:, 2)   = cl(:, 2) - 0.5;
cu(:,1)    = cu(:,1)+.5;

x = linspace(0, 1, m)';
l = (x < 0.5); u = ~l;
for i = 1:3
  h(l, i) = interp1(cl(:, 1), cl(:, i+1), x(l));
  h(u, i) = interp1(cu(:, 1), cu(:, i+1), x(u));
end
c = hsv2rgb(h);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c = ikelvinr(m)
%  pos    hue   sat   value
cu = [
  0.0     1/2   0     1.0
  0.125   1/2   0.6   0.95
  0.375   2/3   1.0   0.8
  0.5     2/3   1.0   0.3
  ];

cl = cu;
cl(:, 3:4) = cl(end:-1:1, 3:4);
cl(:, 2)   = cl(:, 2) - 0.5;
cu(:,1)    = cu(:,1)+.5;

x = linspace(0, 1, m)';
l = (x < 0.5); u = ~l;
for i = 1:3
  h(l, i) = interp1(cl(:, 1), cl(:, i+1), x(l));
  h(u, i) = interp1(cu(:, 1), cu(:, i+1), x(u));
end
c = hsv2rgb(h);

c = flipud(c);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef (InferiorClasses = {?chebfun}) domain < double
%DOMAIN   Utility class for CHEBFUN. Mostly for backward compatibility.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DOMAIN Class Description:
%
% DOMAIN inherits from a standard Matlab DOUBLE. A domain object only
% contains vector for the endpoints and breakpoints of the interval it
% represents. This class is lightly documented, since it is mostly intended
% for backward compatibility.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS PROPERTIES:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties ( Access = public )

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS CONSTRUCTOR:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        function obj = domain(varargin)
            %Constructor for the DOMIAN class.

            % Return an empty DOMAIN on null input:
            if ( nargin == 0 )
                data = [];
            else
                data = horzcat(varargin{:});
            end

            % Create the domain:
            obj = obj@double(data);

        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        function out = validate(dom)
            out = true;
            if ( ~isnumeric(dom) || any(isnan(dom)) )
                out = false;
            elseif ( isempty(dom) )
                return
            elseif ( size(dom, 1) > 1 || size(dom, 2) < 2 )
                out = false;
            elseif ( any(diff(dom) <= 0) )
                out = false;
            end
            if ( out == false && nargout == 0 )
                error('CHEBFUN:DOMAIN:domain:invalid', ...
                    'Ends must be a 1xM vector of ordered doubles.');
            end
        end

        function display(dom)
            disp(double(dom))
        end

        function out = subsref(d, s)
            out = subsref@double(d, s);
            out = double(out);
        end

        function varargout = sprintf(varargin)
            % This is required as built-in subsref does not know what to do with
            % a DOMAIN object.

            varargin = domain2double(varargin{:});

            % Call built-in SPRINTF:
            varargout{1:nargout} = sprintf(varargin{:});
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% STATIC METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = true )

        % Merge two domains:
        newDom = merge(varargin)

        function varargin = toDouble(varargin)
            % Cast DOMAIN to DOUBLE:
            for k = 1:nargin
                if ( isa(varargin{k}, 'domain') )
                    varargin{k} = double(varargin{k});
                end
            end
        end

    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function cluster_length = save_declus(catalog)

    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    global cluslength %[OUT]
    %FIXME: this is torn to shreds by the updated catalog format
    report_this_filefun();
    storedcat=a;
    hodis = fullfile(hodi, 'external');
    cd(hodis);

    str = [];

    s = [  ZG.primeCatalog.Date.Year  ZG.primeCatalog.Date.Month  ZG.primeCatalog.Date.Day ZG.primeCatalog.Date.Hour ZG.primeCatalog.Date.Minute ZG.primeCatalog.Magnitude   ZG.primeCatalog.Latitude  ZG.primeCatalog.Longitude  ZG.primeCatalog.Depth   ];
    fid = fopen(['data'],'w') ;
    fprintf(fid,'%4.0f%2.0f%2.0f%2.0f%2.0f  %3.1fmb%7.3f%8.3f%5.1fA\n',s');
    fclose(fid);
    clear s

    s = [taumin*60*24 taumax*60*24 P xk xmeff rfact err derr ];
    fid = fopen(['para.dat'],'w') ;
    fprintf(fid,'%5.0f  %5.0f  %5.3f  %5.3f  %5.3f  %5.3f  %5.3f  %5.3f\n',s');
    fclose(fid);
    clear s

    % This executes the clus.exe FORTRAN code
    unix(['.' filesep 'myclus ']);

    %open datafile
    fid = 'outf.clu';

    try
        format = ['%12c %3f %f %f %f %d'];
        %[dat,mag,lat,lon,dep,clu] = textread(fid,format,'whitespace',' \b\r\t\n mb A ');
        C = textscan(fid, format, 'Whitespace',' \b\r\t\n mb A ');
        dat=C{1}; mag=C{2}; lat = C{3}; lon=C{4}; dep = C{5}; clu = C{6};
    catch ME
        l=ME.message;
        l1 = strfind(l,',');
        anz = str2double(l(53:l1-1));
        %[dat,mag,lat,lon,dep,clu] = textread(fid,format,anz-1,'whitespace',' \b\r\t\n mb A ');
        C = textscan(fid, format, anz-1, 'Whitespace',' \b\r\t\n mb A ');
        dat=C{1}; mag=C{2}; lat = C{3}; lon=C{4}; dep = C{5}; clu = C{6};

        disp(['Error in Line ' num2str(anz) ' read only lines  1 - ' num2str(anz-1) ]);

    end

    %transform data to ZMAP format
    watchon;
    disp('Reloading data ...')

    yr = str2double(dat(:,1:4));
    mo=  str2double(dat(:,5:6));
    da=  str2double(dat(:,7:8));
    hr=  str2double(dat(:,9:10));
    mi=  str2double(dat(:,11:12));

    replaceMainCatalog([lon lat ZG.primeCatalog.Date mo da mag storedcat.Depth hr mi clu]);

    cluslength=[];
    n=0;
    k1=max(clu);
    for j=1:k1                         %for all clusters
        cluslength(j)=length(find(clu==j));  %length of each clusters
    end

    tmp=find(cluslength);      %numbers of clusters that are not empty

    %cluslength,bg,mbg only for events which are not zero
    cluslength=cluslength(tmp);

    clustnumbers=(1:length(tmp));    %stores numbers of clusters
    l = a(:,10) > 0;
    clus = ZG.primeCatalog.subset(l);
    ZG.primeCatalog(l,:) = [];

    % plot the results
    zmap_update_displays();
    set(gca,'NextPlot','add')
    plot(clus(:,1),clus(:,2),'m+');

    st1 = [' The declustering found ' num2str(max(clu)) ' clusters of earthquakes, a total of '...
        ' ' num2str(length(clus(:,1))) ' events (out of ' num2str(storedcat.Count) '). '...
        ' The map window now display the declustered catalog containing ' num2str(ZG.primeCatalog.Count) ' events . The individual clusters are displayed as magenta o in the map. ' ];

    msgbox(st1,'Declustering Information')
    watchoff;
end

author: CelsoReyes


text:
function density = p05_density ( m, n, point )

%*****************************************************************************80
%
%% P05_DENSITY returns the density for problem 05.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real POINT(M,N), the coordinates of the points.
%
%    Output, real DENSITY(N), the mesh density at each point.
%
  density = ones ( 1, n );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function h = splitROIleftright(h)
h.dat.cl.iscell = [h.dat.stat.iscell];

% h.dat.cl.iscell = double(h.dat.cl.iscell);

% overwrite manual selections
% h.dat.cl.iscell(h.dat.cl.manual>1e-3) = 1;
% h.dat.cl.iscell(h.dat.cl.manual<-1e-3) = 0;
%%
IScell = zeros(h.dat.cl.Ly, h.dat.cl.Lx);
ix = h.dat.res.iclust>0;
IScell(ix) = h.dat.cl.iscell(h.dat.res.iclust(ix));
h.dat.cl.k1 = reshape(IScell, h.dat.cl.Ly, h.dat.cl.Lx);

author: cortex-lab


text:
function [g,A] = minusOne_v2(f)
%minusOne Multiplies an input array by (-1)^x+y.
%   [G,A] = minusOne(F) multiplies F by (-1)^x+y to produce G. F can be
%   a 1-D or 2-D array. On the output, G = (-1)^(x+y).*F, and A is the
%   array of values (-1)^(x+y). Output G is floating point with values
%   in the same numerical range as F.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Preliminaries.

% Convert input to floating point for multiplication later.
f = im2double(f);

% Image size.
[M,N] = size(f);

x = (0:M - 1)';
y = 0:N - 1;

A = (-1).^(x + y);

% Multiply input by matrix.
g = A.*f;

author: dipum

repo: dipum-toolbox

save_path: github-repos/MATLAB/dipum-dipum-toolbox


text:
%% Stress Tensors
%
%% PLEASE HELP AND ADD CONTENT HERE
%

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function A = create_2d_image_downsample_matrix(szin, k)
    szout = ceil(szin / k);

    npixin = prod(szin);
    npixout = prod(szout);

    A = sparse(npixout, npixin);
    [m1, m2] = ndgrid(1:szin(1), 1:szin(2));

    mm1 = ceil(m1 / k);
    mm2 = ceil(m2 / k);

    m1 = m1(:); m2 = m2(:);
    mm1 = mm1(:); mm2 = mm2(:);

    rr = (mm2-1)*szout(1) + mm1;
    cc = (m2-1)*szin(1) + m1;
    A = sparse(rr, cc, 1/k^2, npixout, npixin);

%     return;

end

repo: pTVreg


text:
function [Torr] = inH2O2Torr(inH2O)
% Convert pressure from inches of water column at 4 degrees to torr
% Chad Greene 2012
Torr = inH2O*1.86832;

author: Sable

repo: mcbench-benchmarks


text:
function c = abs(a)
	c = a;
end


text:
function [ n_data, x, fx ] = cin_values ( n_data )

%*****************************************************************************80
%
%% CIN_VALUES returns some values of the alternate cosine integral function.
%
%  Discussion:
%
%    The alternate cosine integral is defined by
%
%      CIN(X) = gamma + log(X) + integral ( 0 <= T <= X ) ( cos ( T ) - 1 ) / T  dT
%
%    In Mathematica, the function can be evaluated by:
%
%      EulerGamma + Log[x] - CosIntegral[x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 16;

  fx_vec = [ ...
     0.6185256314820045E-01, ...
     0.8866074809482194E-01, ...
     0.1200260139539026E+00, ...
     0.1557934976348559E+00, ...
     0.1957873187759337E+00, ...
     0.2398117420005647E+00, ...
     0.3390780388012470E+00, ...
     0.4516813164280685E+00, ...
     0.5754867772153906E+00, ...
     0.7081912003853150E+00, ...
     0.8473820166866132E+00, ...
     0.1207635200410304E+01, ...
     0.1556198167561642E+01, ...
     0.1862107181909382E+01, ...
     0.2104491723908354E+01, ...
     0.2274784183779546E+01 ];

  x_vec = [ ...
     0.5E+00, ...
     0.6E+00, ...
     0.7E+00, ...
     0.8E+00, ...
     0.9E+00, ...
     1.0E+00, ...
     1.2E+00, ...
     1.4E+00, ...
     1.6E+00, ...
     1.8E+00, ...
     2.0E+00, ...
     2.5E+00, ...
     3.0E+00, ...
     3.5E+00, ...
     4.0E+00, ...
     4.5E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
% RPCA (De la Torre and Black, 2001)
% process_video('RPCA', 'RPCA', 'dataset/demo.avi', 'output/demo_RPCA.avi');
sizeim = [params.rows params.cols];
[L,S] = RPCA(M,sizeim);

author: andrewssobral

repo: lrslibrary


text:
% Contents of JPRINTF
%
% jprintf.m -- main routine
% jprintf_demo.m -- example script showing how to use JPRINTF
% JavaBoxLibrary.mdl -- Simulink model library for formated text I/O
% test.mdl -- example Simulink model showing use of Simulink text I/O blocks
% jprintsfun.m -- Simulink S-function providing interface to jprintf.m
%


text:
%% BIHARMONICMIXEDFEM solves the biharmonic equation
%
%     laplace^2 u = f;   [0,1]^2
%               u = g_D;
%            Dn u = g_N.
%
% by writing in a mixed form
%
%       w = Lap u
%   Lap w = f
%
% The boundary condition u = g_D is imposed into the space (Dirichlet) and
% the boundary condition Dn u = g_N is imposed as a Newmann boundary
% condition.
%
% The rate of convergence for u is optimal but for w is sub-optimal. For
% linear element, optimal order for w is also observed.
%
% Created by Jie Zhou. Clean up by Long Chen. Further clean up on
% biharmonic functions are needed and multigrid solvers should be included.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

close all
clear variables;

%% Parameters
maxIt = 4;
N = zeros(maxIt,1);
h = zeros(maxIt,1);
erruL2 = zeros(maxIt,1);
erruH1 = zeros(maxIt,1);
errwL2 = zeros(maxIt,1);
errwH1 = zeros(maxIt,1);

%%  Generate an initial mesh
[node,elem] = squaremesh([0 1 0 1], 0.25);
bdFlag = setboundary(node,elem,'Dirichlet');    % Dirichlet boundary condition

for k = 1:1
    [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
%     [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
end

%% Set up PDE data
pde = biharmonicdata;

%% Finite Element Method
for k = 1:maxIt
    % refine mesh
   [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
   [w,u] = biharmonicP1(node,elem,bdFlag,pde);
%    [w,u] = biharmonicP2(node,elem,bdFlag,pde);
%    [w,u] = biharmonicP3(node,elem,bdFlag,pde);
   N(k) = size(w,1)+size(u,1);
   h(k) = 1./(sqrt(size(node,1))-1);
   erruL2(k) = getL2error(node,elem,pde.exactu,u);
   erruH1(k) = getH1error(node,elem,pde.Du,u);
   errwL2(k) = getL2error(node,elem,pde.exactw,w);
   errwH1(k) = getH1error(node,elem,pde.Dw,w);
end

%% Plot convergence rates and display error table
figure;
subplot(1,2,1);
showrateh2(h,erruH1,1,'-*','|| Du - Du_h||',...
           h,erruL2,1,'k-+','|| u - u_h||');
subplot(1,2,2)
showrateh2(h,errwH1,1,'c-*','|| Dw - Dw_h||',...
           h,errwL2,1,'m-+','|| w - w_h||');

fprintf('\n');
disp('Table: Error')
colname = {'#Dof','h','||u-u_h||','||Du-Du_h||','||Dw-Dw_h||','||w - w_h||'};
disptable(colname,N,[],h,'%0.3e',erruL2,'%0.5e',erruH1,'%0.5e',errwH1,'%0.5e',errwL2,'%0.5e');

repo: ifem


text:
%% immesh
% Below is a demonstration of the features of the |immesh| function

%%
clear; close all; clc;

%% Syntax
% |[I,J,K]=immesh(M);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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 3 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, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function svgFlood(s, color, opacity, result)
% Adds a feFlood SVG filter
% PRELIMINARY IMPLEMENTATION (Parameters may change)
%
% svgFlood(s, color, opacity, result)
% Parameters:
%   s : Array of plot object handles
%   color  : A 1x3 color vector [0...1[
%   opacity : Opacity value [0...1]
%   result : String that identifies the filter result for following filter
%            stages.
for i = 1:length(s)
    userdata = get(s(i),'UserData');
    if isfield(userdata, 'svg') && isfield(userdata.svg, 'Filter')
        next = length(userdata.svg.Filter) + 1;
    else
        next = 1;
    end
    userdata.svg.Filter(next).Subfilter.Type = 'feFlood';
    userdata.svg.Filter(next).Subfilter.Color = color;
    userdata.svg.Filter(next).Subfilter.Opacity = opacity;
    userdata.svg.Filter(next).Subfilter.Result = result;
    set(s(i),'UserData', userdata);
end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
% test find root node
clear;clc;close all;
path(path,'toolbox') ;
t3 = 6; % for small cycles;
% sk_filename='../result/cylinder1_contract_t(2)_nn(14)_WL(10.633697)_WH(1.000000)_sl(3.000000)_skeleton.mat';
% sk_filename='../result/simplejoint_v4770_contract_t(3)_nn(30)_WL(15.378798)_WH(1.000000)_sl(3.000000)_skeleton.mat';
sk_filename='../result/horse_v1987_contract_t(3)_nn(24)_WL(7.786614)_WH(1.000000)_sl(3.000000)_skeleton.mat';

load(sk_filename,'M');

%%
[joints, segments] = find_joints(M,false);
[joints, segments] = remove_small_cycles(M, joints, segments,t3, true);

author: taiya

repo: cloudcontr

save_path: github-repos/MATLAB/taiya-cloudcontr


text:
function eps = deformationTensor(sS)
% deformation tensor
%
% Syntax
%   eps = deformationTensor(sS)
%
% Input
%  sS - list of @slipSystem
%
% Output
%  eps - deformationTensor
%

eps = dyad(sS.b.normalize,sS.n.normalize);


text:
function [ w, f ] = jed_to_weekday ( jed )

%*****************************************************************************80
%
%% JED_TO_WEEKDAY computes the day of the week from a JED.
%
%  Discussion:
%
%    BC 4713/01/01 => JED = 0.0 was noon on a Monday.
%
%    jedmod = mod ( 0.0D+00, 7.0D+00 ) = 0.0D+00
%    j = mod ( nint ( 0 ), 7 ) = 0
%    f = ( 0.0D+00 + 0.5D+00 ) - real ( j ) = 0.5D+00
%    w = i4_wrap ( 0 + 2, 1, 7 ) = 2 = MONDAY
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 April 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Richards,
%    Mapping Time, The Calendar and Its History,
%    Oxford, 1999.
%
%  Parameters:
%
%    Input, real JED, the Julian Ephemeris Date.
%
%    Output, integer W, the day of the week of the date.
%    The days are numbered from Sunday through Saturday, 1 through 7.
%
%    Output, real F, the fractional part of the day.
%
  jedmod = mod ( jed, 7.0 );

  j = mod ( floor ( jedmod ), 7 );

  f = ( jedmod + 0.5 ) - j;

  w = i4_wrap ( j + 2, 1, 7 );

  return
end


text:
classdef UniformGraphFuxin < AbstractGraph
%UNIFORMGRAPHFUXIN is a variant of UniformGraph, which generates graphs
% with seeds and unaries which are not dependent on image information (for
% instance unaries are only based on superpixel size). This class is
% derived from AbstractGraph. To undestand how unaries are generated, its
% worth seeing how AbstractGraph.prepare_graphs() calls create_unaries() in
% this class to compute unary values.
%
% @authors:     Fuxin Li
% @contact:     ahumayun@cc.gatech.edu
% @affiliation: Georgia Institute of Technology
% @date:        Fall 2013 - Summer 2014

    properties
    end

    methods
        function obj = UniformGraphFuxin(seg_obj, segm_index)
            obj.init(seg_obj, segm_index);
        end

        function create_unaries(obj, graph_sub_method, seed_sets)
            switch graph_sub_method
                case 'internal'
                    obj.generate_uniform_unaries('internal', seed_sets);
                case 'external'
                    obj.generate_uniform_unaries('external', seed_sets);
                case 'external2'
                    obj.generate_uniform_unaries('external2', seed_sets);
                case 'subframe'
                    error('UniformGraphFuxin:compute_graph', ...
                          '''subframe'' not implemented yet');
                otherwise
                    error('UniformGraphFuxin:compute_graph', ...
                          '''%s'' is an invalid graph method', ...
                          graph_sub_method);
            end
        end

        function generate_uniform_unaries(obj, internal_external, ...
                                          seed_sets)
            inside_frame = ~obj.seg_obj.sp_frame_set;

            funcs = create_unary_aux_funcs(obj, internal_external);

            % find the offset according to the number of unary graphs
            % already added to graph_unaries_all
            ofst = size(obj.graph_unaries_all.nonlambda_s,2);

            % initialize the matrices for holding all the unaries (see the
            % comment for this structure in AbstractGraph.prepare_graphs())
            obj.graph_unaries_all.nonlambda_s = ...
                [obj.graph_unaries_all.nonlambda_s, zeros(size(seed_sets))];
            obj.graph_unaries_all.nonlambda_t = ...
                [obj.graph_unaries_all.nonlambda_t, zeros(size(seed_sets))];
            obj.graph_unaries_all.lambda_s = ...
                [obj.graph_unaries_all.lambda_s, zeros(size(seed_sets))];
            obj.graph_unaries_all.lambda_t = ...
                [obj.graph_unaries_all.lambda_t, zeros(size(seed_sets))];

            % iterate over all seed locations
            for i = 1:size(seed_sets, 2)
                if strcmp(internal_external, 'internal')
                    [nonlambda_s, lambda_s, lambda_t] = ...
                        internal_uniform(obj, seed_sets(:,i), ...
                                         funcs.parametric_func);

                    obj.graph_unaries_all.nonlambda_s(:,ofst+i) = ...
                                                               nonlambda_s;
                    obj.graph_unaries_all.lambda_s(:,ofst+i) = lambda_s;
                    obj.graph_unaries_all.lambda_t(:,ofst+i) = lambda_t;
                elseif strcmp(internal_external, 'external')
                    [nonlambda_s, nonlambda_t, lambda_s, lambda_t] = ...
                        external_uniform(obj, seed_sets(:,i), ...
                                         inside_frame, ...
                                         funcs.parametric_func, ...
                                         funcs.frame_cost_func);

                    obj.graph_unaries_all.nonlambda_s(:,ofst+i) = ...
                                                               nonlambda_s;
                    obj.graph_unaries_all.nonlambda_t(:,ofst+i) = ...
                                                               nonlambda_t;
                    obj.graph_unaries_all.lambda_s(:,ofst+i) = lambda_s;
                    obj.graph_unaries_all.lambda_t(:,ofst+i) = lambda_t;
                else
                    [nonlambda_s, lambda_s, lambda_t] = ...
                        external_uniform2(obj, seed_sets(:,i), ...
                                          inside_frame, ...
                                          funcs.parametric_func, ...
                                          funcs.frame_cost_func);

                    obj.graph_unaries_all.nonlambda_s(:,ofst+i) = ...
                                                               nonlambda_s;
                    obj.graph_unaries_all.lambda_s(:,ofst+i) = lambda_s;
                    obj.graph_unaries_all.lambda_t(:,ofst+i) = lambda_t;
                end
            end
        end

        function [nonlambda_s, lambda_s, lambda_t] = ...
                internal_uniform(obj, the_rectangle, parametric_func, ...
                                 frame_cost_func)
            % see the comment for these output matrices in
            % AbstractGraph.prepare_graphs()
            nonlambda_s = zeros(length(the_rectangle),1);
            lambda_s = zeros(length(the_rectangle),1);
            lambda_t = zeros(length(the_rectangle),1);

            bground_pixels = ~the_rectangle;
%            frame_cost = frame_cost_func(the_rectangle);
% Internal is always without frame_cost
            nonlambda_s(the_rectangle) = inf;
            lambda_t(bground_pixels) = parametric_func(bground_pixels) / 2; %* obj.graph_seed_frame_weight;
            lambda_s(bground_pixels) = parametric_func(bground_pixels); %*obj.graph_seed_frame_weight;
        end

        function [nonlambda_s, nonlambda_t, lambda_s, lambda_t] = ...
                external_uniform(obj, the_rectangle, inside_frame, ...
                                 parametric_func, frame_cost_func)
            % see the comment for these output matrices in
            % AbstractGraph.prepare_graphs()
            nonlambda_s = zeros(length(the_rectangle),1);
            nonlambda_t = zeros(length(the_rectangle),1);
            lambda_s = zeros(length(the_rectangle),1);
            lambda_t = zeros(length(the_rectangle),1);

            inside_frame_but_rect = ~the_rectangle;
            lambda_weights = parametric_func(inside_frame_but_rect); %* obj.graph_seed_frame_weight;
%             bground_pixels = ~(inside_frame | obj.seg_obj.sp_frame_set);

            nonlambda_s(the_rectangle) = inf;
            lambda_s(inside_frame_but_rect) = lambda_weights;

            frame_cost = frame_cost_func(the_rectangle);
            % With frame cost
            nonlambda_t(obj.seg_obj.sp_frame_set) = ...
                frame_cost * 2 .* parametric_func(obj.seg_obj.sp_frame_set);
            lambda_t(inside_frame_but_rect) = lambda_weights / 256;
            % no frame edge if the frame superpixel is seed
            nonlambda_t(the_rectangle) = 0;
        end

        function [nonlambda_s, lambda_s, lambda_t] = ...
                external_uniform2(obj, the_rectangle, inside_frame, ...
                                  parametric_func, frame_cost_func)
            % see the comment for these output matrices in
            % AbstractGraph.prepare_graphs()
            nonlambda_s = zeros(length(the_rectangle),1);
            lambda_s = zeros(length(the_rectangle),1);
            lambda_t = zeros(length(the_rectangle),1);

            inside_frame_but_rect = ~the_rectangle;
            lambda_weights = parametric_func(inside_frame_but_rect); %* obj.graph_seed_frame_weight;

            % W/O frame cost
            nonlambda_s(the_rectangle) = inf;
            lambda_s(inside_frame_but_rect) = lambda_weights;
            lambda_t(inside_frame_but_rect) = lambda_weights / 256;
            % TODO: Want to add horiz/vert, these 2 are useful.
        end
    end
end

author: Cloud-CV

repo: object-proposals


text:
function [ final_dist ] = re_ranking( feat, M, W, query_num, k1, k2, lambda)
% k-reciprocal re-ranking

%% initial ranking list
original_dist = MahDist(M, feat' * W, feat' * W);
[~, initial_rank] = sort(original_dist, 2, 'ascend');
gallery_num = size(original_dist,2);
%% compute k-reciprocal feature
V = zeros(size(original_dist), 'single');
original_dist = original_dist./ repmat(max(original_dist, [], 2), 1, size(original_dist, 2));
for i = 1:size(original_dist, 1)
    %% k-reciprocal neighbors
    forward_k_neigh_index = initial_rank(i, 1: k1 + 1);
    backward_k_neigh_index  = initial_rank(forward_k_neigh_index, 1: k1+1);
    [fi, ~, ~]= find(backward_k_neigh_index == i);
    k_reciprocal_index  = forward_k_neigh_index(fi);
    k_reciprocal_expansion_index = k_reciprocal_index;
    for j = 1: length(k_reciprocal_index)
        candidate = k_reciprocal_index(j);
        candidate_forward_k_neigh_index = initial_rank(candidate, 1: round((k1+1)/2));
        candidate_backward_k_neigh_index = initial_rank(candidate_forward_k_neigh_index, 1: round((k1+1)/2));
        [fi_candidate, ~, ~]= find(candidate_backward_k_neigh_index == candidate);
        candidate_k_reciprocal_index = candidate_forward_k_neigh_index(fi_candidate);
        if length(intersect(k_reciprocal_index, candidate_k_reciprocal_index)) > 2/3*length(candidate_k_reciprocal_index)
            k_reciprocal_expansion_index = [k_reciprocal_expansion_index candidate_k_reciprocal_index];
        end
    end
    k_reciprocal_expansion_index = unique(k_reciprocal_expansion_index);
    %% feature encoding
    weight = exp(-original_dist(i, k_reciprocal_expansion_index));
    V(i, k_reciprocal_expansion_index) = weight/sum(weight);
end
%% local query expansion
if k2 ~=1
    V_qe = zeros(size(V), 'single');
    for i = 1:size(V, 1)
        V_qe(i, :) = single(mean(V(initial_rank(i, 1:k2), :), 1));
    end
    V = V_qe;
    V_qe = [];
end

%% Inverted Index
% Inpsired by Song Bai, Xiang Bai,
% Sparse Contextual Activation for Efficient Visual Re-ranking, TIP, 2016.
% We apply the inverted index to quickly compute the Jaccard distance.

invIndex = cell(gallery_num, 1);
for i = 1:gallery_num
    invIndex{i} = find(V(:, i) ~=0);
end

jaccard_dist = zeros(size(original_dist), 'single');

for i = 1:query_num
    temp_min = zeros(1, gallery_num, 'single');
    indNonZero = find( V( i, : ) ~= 0 );
    indImages = invIndex( indNonZero );
    for j = 1 : length( indNonZero )
        temp_min( 1, indImages{j} ) = temp_min( 1, indImages{j} )...
            + single( min( V(i, indNonZero(j)), V(indImages{j}, indNonZero(j)) ) )';
    end
    jaccard_dist(i, :) = bsxfun(@minus, 1, temp_min./(2 - temp_min));
end

final_dist = jaccard_dist*(1-lambda) + original_dist*lambda;
final_dist = final_dist(1:query_num,query_num+1:end)';

end

repo: Learning-via-Translation


text:
function [i,P] = spm_voice_i(str)
% Get indices, word strings or priors from lexicon
% FORMAT [str] = spm_voice_i(i)
% FORMAT [i  ] = spm_voice_i(str)
% FORMAT [i,P] = spm_voice_i(str)
%
% str  - string or cell array
% i    - index in lexicon (VOX.LEX)
% P    - corresponding array of prior probabilities

% requires the following in the global variable VOX:
% LEX  - lexical structure array
%
%  This routine returns the index or indices of a word if supplied with a
%  string or cell array. Alternatively, it returns the string corresponding
%  to an index or vector of indices. if called with the second argument, it
%  returns a prior probability matrix where the specified words (for
%  strings) have a prior odds ratio of 64.
%
%  NB: If a string is not in the lexicon, 0 is returned.
%__________________________________________________________________________
% Copyright (C) 2019 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_voice_i.m 7750 2019-12-05 17:54:29Z spm $

% get timeseries from audio recorder(or from a path
%--------------------------------------------------------------------------

% words in lexicon
%==========================================================================
global VOX
word  = {VOX.LEX.word};                  % words in lexicon

% return cell array of indexed words
%--------------------------------------------------------------------------
if isnumeric(str)
    i = word(str);
    return
end

% return indices
%--------------------------------------------------------------------------
if iscell(str)
    if nargout < 2
        for w = 1:numel(str)
            i(w) = spm_voice_i(str{w});
        end
        return
    else

        % return indices and prior probabilities
        %------------------------------------------------------------------
        nw    = numel(str);
        P     = zeros(numel(VOX.LEX),nw) + 1/64;
        for w = 1:nw
            i      = spm_voice_i(str{w});
            P(i,w) = 1;
        end

        % sum to one constraint
        %------------------------------------------------------------------
        P     = bsxfun(@rdivide,P,sum(P));
        return
    end
end

% get indices
%--------------------------------------------------------------------------
i    = find(strcmp(str,word));
if isempty(i), i = 0; end

author: spm


text:
function res = imThumbnail(img, siz, varargin)
%IMTHUMBNAIL Resize an image to bound the size in each direction
%
%   RES = imThumbnail(IMG, SIZ)
%   Resizes the image given in IMG such that the dimension is lower or
%   equal to the dimension given by SIZ.
%
%   Example
%     img = imread('cameraman.tif');
%     img2 = imThumbnail(img, [64 64]);
%     imshow(img2)
%
%     img = imread('peppers.png');
%     img2 = imThumbnail(img, [64 64]);
%     imshow(img2)
%
%   See also
%   imresize
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2013-02-14,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2013 INRA - Cepia Software Platform.

k1 = siz(1) / size(img, 1);
k2 = siz(2) / size(img, 2);
k = min(k1, k2);

res = imresize(img, k, varargin{:});

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function output = prefilt(img, fc)
% ima = prefilt(img, fc);
% fc  = 4 (default)
%
% Input images are double in the range [0, 255];
% You can also input a block of images [ncols nrows 3 Nimages]
%
% For color images, normalization is done by dividing by the local
% luminance variance.

if nargin == 1
    fc = 4; % 4 cycles/image
end

w = 5;
s1 = fc/sqrt(log(2));

% Pad images to reduce boundary artifacts
img = log(img+1);
img = padarray(img, [w w], 'symmetric');
[sn, sm, c, N] = size(img);
n = max([sn sm]);
n = n + mod(n,2);
img = padarray(img, [n-sn n-sm], 'symmetric','post');

% Filter
[fx, fy] = meshgrid(-n/2:n/2-1);
gf = fftshift(exp(-(fx.^2+fy.^2)/(s1^2)));
gf = repmat(gf, [1 1 c N]);

% Whitening
output = img - real(ifft2(fft2(img).*gf));
clear img

% Local contrast normalization
localstd = repmat(sqrt(abs(ifft2(fft2(mean(output,3).^2).*gf(:,:,1,:)))), [1 1 c 1]);
output = output./(.2+localstd);

% Crop output to have same size than the input
output = output(w+1:sn-w, w+1:sm-w,:,:);

author: CSAILVision

repo: LabelMeToolbox


text:
%% Slip Systems
%
%%
%


text:
function [SD1, SD2, SD1_SD2_ratio] = EvalPoincareOnWindows(rr, t_rr, HRVparams, tWin, sqi)

%   OVERVIEW:
%       Calculates SD1 SD2 and SD1_SD2_ratio features from Poincare plot
%       for each defined windows
%   INPUT
%       rr           - (seconds) rr intervals
%       t_rr         - (seconds)  time stamp of rr intervals
%       HRVparams    - struct of settings for hrv_toolbox analysis
%       sqi          - Signal Quality Index; Requires a matrix with
%                      at least two columns. Column 1 should be
%                      timestamps of each sqi measure, and Column 2
%                      should be SQI on a scale from 0 to 1.
%       tWin         - Starting time of each windows to analyze
%
%   OUTPUTS:
%       SD1           : (ms) standard  deviation  of  projection  of  the
%                       PP  on  the line perpendicular to the line of
%                       identity (y=-x)
%       SD2           : (ms) standard deviation of the projection of the PP
%                       on the line of identity (y=x)
%       SD1_SD2_ratio : (ms) SD1/SD2 ratio
%
%
%   Written by: Giulia Da Poian <giulia.dap@gmail.com>
%	REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%	COPYRIGHT (C) 2016
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information
%

% Make vector a column
rr = rr(:);

if nargin < 4
    error('no data provided')
end
if nargin <5 || isempty(sqi)
    sqi(:,1) = t_rr;
    sqi(:,2) = ones(length(t_rr),1);
end

windowlength = HRVparams.windowlength;
SQI_th = HRVparams.sqi.LowQualityThreshold;        % SQI threshold
WinQuality_th = HRVparams.RejectionThreshold; % Low quality windows threshold

% Preallocation (all NaN)

SD1 = ones(length(tWin),1)*NaN;
SD2 = ones(length(tWin),1)*NaN;
SD1_SD2_ratio = ones(length(tWin),1)*NaN;

% Run PoincareMetrics by Windows
% Loop through each window of RR data
for i_win = 1:length(tWin)
    if ~isnan(tWin(i_win))
        % Isolate data in this window
        sqi_win = sqi( sqi(:,1) >= tWin(i_win) & sqi(:,1) < tWin(i_win) + windowlength,:);
        nn_win = rr( t_rr >= tWin(i_win) & t_rr < tWin(i_win) + windowlength );
        lowqual_idx = find(sqi_win(:,2) < SQI_th);         % Analysis of SQI for the window
        % If enough data has an adequate SQI, perform the calculations
        if numel(lowqual_idx)/length(sqi_win(:,2)) < WinQuality_th
            [SD1(i_win), SD2(i_win), SD1_SD2_ratio(i_win)] = PoincareMetrics(nn_win);
        end % end of conditional statements run when SQI is adequate
    end % end of check for sufficient data
end % end of loop through windows

end % end of function

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
% TEST_cannonDynamics.m
%
% This script runs a simple simulation to demonstrate the dynamics of the
% cannon ball problem.
%

clc; clear;

% Cannon ball starts at the origin
x0 = 0;
y0 = 0;

% Prescribe initial conditions with velocity and angle
nSim = 20;
v0 = logspace(-2,6,nSim);
th0 = 45*(pi/180);

% Quadratic drag coefficient:
c = 0.5;

% Store the range as a function of initial angle:
xRange = zeros(1,nSim);

for iSim=1:nSim

% Set up initial conditions for ode45
dx0 = v0(iSim)*cos(th0);
dy0 = v0(iSim)*sin(th0);

% Set up arguments for ode45
userFun = @(t,z)cannonDynamics(t,z,c);  %Dynamics function
tSpan = [0,100];  %Never plan on reaching final time
z0 = [x0;y0;dx0;dy0];  %Initial condition
options = odeset('Events',@groundEvent,'Vectorized','on','AbsTol',1e-8, 'RelTol',1e-8);

% Run a simulation
sol = ode45(userFun, tSpan, z0, options);

xRange(iSim) = sol.y(1,end);

fprintf('Initial Speed: %6.6g,  Horizontal Range: %6.6g \n',v0(iSim),xRange(iSim));

end

% Plot the range as a function of the initial speed
plot(v0,xRange,'kx','LineWidth',2,'MarkerSize',10)
set(gca, 'Xscale', 'log')
set(gca, 'Yscale', 'log')
xlabel('Initial Speed')
ylabel('Horizontal Range');
title('Range vs Initial Speed');

author: MatthewPeterKelly

repo: dscTutorials


text:
function VBA_disp (str, options, carriageReturn)
% // VBA toolbox //////////////////////////////////////////////////////////
%
% VBA_disp (str, options)
% display string or cell-arry of strings according to verbose option
%
% IN:
%   - str: string or cell array of strings to display
%   - option structure containing the verbose flag (default: true)
%
% /////////////////////////////////////////////////////////////////////////

% get flag
try
    verbose = options.verbose;
catch
    verbose = true;
end

if nargin < 3
    carriageReturn = true;
end

% conditional display function
if verbose
    if iscell (str)
        for i = 1 : numel(str)
            fprintf(1, str{i});
        end
    else
       fprintf (1, str);
    end
    if carriageReturn
        fprintf ('\n');
    end
end

author: MBB-team

repo: VBA-toolbox


text:
function [xInt,t]=findEllipsLineIntersect(xc,A,x0,a,aType)
%%FINDELLIPSLINEINTERSECT Given the parameters for an ellipsoid (in 2D an
%       ellipse) of the form (x-xc)'*A*(x-xc)=1 and a line either given
%       parameterically as x=x0+a*t or whether x0 and a are two points on
%       the line, find the points of intersection between the ellipsoid and
%       the line.
%
%INPUTS: xc The xDimX1 center of the ellipsoid.
%         A The xDimXxDim symmetric matrix defining the above equation for
%           an ellipsoid.
%     x0, a These xDimX1 values define a line and the type of line is given
%           by aType. If aType=0, then the line is parameteric as x=x0+a*t
%           where t is a scalar value. If aType=1, then x0 and a are both
%           points on the line.
%     aType As mentioned above, this specifies how the line is
%           parameterized. The default if omitted or an empty matrix is
%           passed is 0.
%
%OUTPUTS: xInt Either xDimX2 (or rarely xDimX1) matrix holding the vector
%              points of intersection of the line and the ellipsoid or an
%              empty matrix if there are no real solutions (the line does
%              not intersect the ellipsoid).
%            t If solutions exist, this is the 2X1 (or rarely 1X1) set of
%              parametric parameters. If aType=0, then these are the t
%              values in the parameteric equation. If aType=1, then these
%              are the t values in the parameteric equation x=x0+(a-x0)*t,
%              in which case values between 0 and 1 indicate that the
%              intersection occured between the two specified points.
%
%Start with the equations x=x0+a*t and (x-xc)'*A*(x-xc)=1. Substitite the
%first into the second and one ends up with a quadratic equation
%c1*t^2+c2*t+c3=0 with c1=a'*A*a, c2=2*xTilde'*A*a, c3=xTilde'*A*xTile-1
%where xTilde=x0-xc.
%
%EXAMPLE:
%Find and plot the intersection points of a line an an ellipse. Also
%consider a line that does not intersect and one sees that nothing is
%plotted from that, because no solution exists.
% %The ellipse.
% xc=[1;2];
% A=[1,0.5;
%    0.5,2];
% %Two points defining a line that does intersect.
% x01=[0;3];
% x02=[2;1];
% %Two points defining a line that does not intersect.
% x11=[3;3];
% x12=[3;1];
%
% aType=1;
% xInt0=findEllipsLineIntersect(xc,A,x01,x02,aType);
% xInt1=findEllipsLineIntersect(xc,A,x11,x12,aType);%This is empty.
%
% figure(1)
% clf
% hold on
% drawEllipse(xc,A,1,'linewidth',2)
% plot([x01(1);x02(1)],[x01(2);x02(2)],'-k','linewidth',2)
% plot([x11(1);x12(1)],[x11(2);x12(2)],'-k','linewidth',2)
%
% numSol0=size(xInt0,2);
% numSol1=size(xInt1,2);%Should be 0.
% for k=1:numSol0
%     scatter(xInt0(1,k),xInt0(2,k),400,'.r')
% end
% for k=1:numSol1
%     scatter(xInt1(1,k),xInt1(2,k),400,'.g')
% end
%
%July 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<5||isempty(aType))
    aType=0;
end

if(aType==1)
    %Make it parameteric.
    a=a-x0;
end

xTilde=x0-xc;

c1=a'*A*a;
c2=2*xTilde'*A*a;
c3=xTilde'*A*xTilde-1;

t=roots([c1;c2;c3]);

if(isempty(t)||any(~isreal(t)))
    t=[];
    xInt=[];
    return
end

xDim=size(x0,1);
numSol=length(t);

xInt=zeros(xDim,numSol,1);
for curSol=1:numSol
    xInt(:,curSol)=x0+a*t(curSol);
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function concurrencyPoints = traceConcurrencyPoints( boundaryPoints, cornerPoints )
%TRACECONCURRENCYPOINTS
%
variance = 0;

sizeOfBoundaryPoints = size(boundaryPoints); % 1 element of the vector is the number of row
sizeOfCornerPoints = size(cornerPoints); % 1 element of the vector is the number of row
numberOfConcurrencyPoints = 0;

numberOfToleratedPoints = 0;

    for i=1:sizeOfBoundaryPoints(1)
        for j=1:sizeOfCornerPoints(1)
            if boundaryPoints(i,1) == cornerPoints(j,1) & boundaryPoints(i,2) == cornerPoints(j,2)
                numberOfConcurrencyPoints = numberOfConcurrencyPoints + 1;
                concurrencyPoints (numberOfConcurrencyPoints,1) = cornerPoints(j,1);
                concurrencyPoints (numberOfConcurrencyPoints,2) = cornerPoints(j,2);
            end

            % tolerance implicated
            if (abs(boundaryPoints(i,1) - cornerPoints(j,1)) + abs(boundaryPoints(i,2) - cornerPoints(j,2))) < variance
            %if (abs(boundaryPoints(i,1) - cornerPoints(j,1)) < variance)
            %    if (abs(boundaryPoints(i,2) - cornerPoints(j,2)) < variance)
                numberOfConcurrencyPoints = numberOfConcurrencyPoints + 1;
                                numberOfToleratedPoints = numberOfToleratedPoints +1;

                concurrencyPoints (numberOfConcurrencyPoints,1) = cornerPoints(j,1);
                concurrencyPoints (numberOfConcurrencyPoints,2) = cornerPoints(j,2);
            %    end
            %end
            end

        end
    end

    disp(numberOfToleratedPoints);
end


text:
function [G1,C,impact,eu,F]=schur_solver(g0,g1,c,psi,pi,continuous,check_exist,check_uniq,n_v)
% Solver based on Schur decomposition (This is based on Sims's gensys, but adjusted
%    to take advantage of special case that applies to heterogeneous agent models).
%
% by SeHyoun Ahn, Dec 2016
%
% REFERENCES:
%    * Sims, Christopher A. "Solving linear rational expectations models."
%         Computational economics 20.1 (2002): 1-20.
%    * Ahn, SeHyoun, Greg Kaplan, Benjamin Moll, Thomas Winberry, and
%         Christian Wolf. "When Inequality Matters for Macro and Macro Matters
%         for Inequality."
%
% PARAMETERS:
%    Check Sims's paper for notation or Our Paper for notation
%
%    XXXXX DO NOT USE DISCRETE TIME VERSION. IT PROBABLY IS WRONG XXXXXX
%    The issue is noted in < https://github.com/gregkaplan/phact/issues/1 >
%
%    continuous = 1 for a continuous time problem (default)
%                0 for a discrete time problem
%
%    check_exist = 1 check for existence of solution (default)
%                 0 existence has been checked before
%    check_uniq = 0 do not check for uniqueness (default)
%                1 check for uniqueness
%    n_v = number of unstable roots
%
% OUPUTS:
%    eu(1) = 1 solution exists
%            0 no stable solution
%           -5 Flag for checking existence is off
%    eu(2) = 1 solution is unique
%            0 solution is not unique
%           -5 Flag for checking uniqueness is off
%
% SYNTAX:
% [G1,C,impact,eu,F]=schur_solver(g0,g1,c,psi,pi,continuous,check_exist,check_uniq,varargin)

%% Set Default Values %%
switch nargin
    case 8
        n_v = -1;
    case 7
        check_uniq = 0;
        n_v = -1;
    case 6
        check_uniq = 0;
        check_exist = 1;
        n_v = -1;
    case 5
        check_uniq = 0;
        check_exist = 1;
        continuous = 1;
        n_v = -1;
    case 1:4
        error('Insufficient Number of Input Variables');
end

realsmall = sqrt(eps)*10;
eu = [-5;-5];
n = size(g1,1);

%% Schur Decomposition
[U,T] = schur(full(g1),'real');
if continuous
    g_eigs = real(ordeig(T));
    nunstab = sum(g_eigs>0);
    if n_v > -1
        [aux,ind] = sort(g_eigs,'descend');
        locs = ones(n,1);
        locs(ind(1:n_v)) = 0;
        [U,T] = ordschur(U,T,locs);
        if nunstab > n_v
            warning('<schur_solver>: There are more than n_v number of positive eigenvalues with smallest values:');
            disp(aux(n_v+1:nunstab));
        elseif nunstab < n_v
            warning('<schur_solver>: There are less than n_v number of positive eigenvalues:');
            disp(aux(nunstab+1:n_v));
        end
        nunstab = n_v;
    else
        [U,T] = ordschur(U,T,'lhp');
    end
else
    g_eigs = abs(ordeig(T));
    nunstab = sum(g_eigs>1);
    if n_v > -1
        [aux,ind] = sort(g_eigs,'descend');
        locs = ones(n,1);
        locs(ind(1:n_v)) = 0;
        [U,T] = ordschur(U,T,locs);
        if nunstab > n_v
            warning('<schur_solver>: There are more than n_v number of positive eigenvalues with smallest values:');
            disp(aux(n_v+1:nunstab));
        elseif nunstab < n_v
            warning('<schur_solver>: There are less than n_v number of positive eigenvalues:');
            disp(aux(nunstab+1:n_v));
        end
        nunstab = n_v;
    else
	[U,T] = ordschru(U,T,'udi');
    end
end

u1 = U(:,1:n-nunstab)';
u2 = U(:,n-nunstab+1:n)';

etawt = u2*pi;
[ueta,deta,veta] = svd(etawt);
md = min(size(deta));
bigev = find(diag(deta(1:md,1:md))>realsmall);
ueta = ueta(:,bigev);
veta = veta(:,bigev);
deta = deta(bigev,bigev);

if check_exist
    zwt = u2*psi;
    [uz,dz,vz] = svd(zwt);
    md = min(size(dz));
    bigev = find(diag(dz(1:md,1:md))>realsmall);
    uz = uz(:,bigev);
    vz = vz(:,bigev);
    dz = dz(bigev,bigev);

    if isempty(bigev)
        eu(1) = 1;
    else
        eu(1) = (norm(uz-ueta*ueta'*uz,'fro') < realsmall*n);
    end

    if (~eu(1) && (n_v == -1))
        warning('<schur_solver>: Solution does not exist');
    end
    impact = real(-pi*veta*(deta\ueta')*uz*dz*vz'+psi);
else
    impact = real(-pi*veta*(deta\ueta')*u2*psi+psi);
end

G1 = U*T*spdiags([ones(n-nunstab,1);zeros(nunstab,1)],0,n,n)*U';
G1 = real(G1);

if check_uniq
    [~,deta1,veta1] = svd(u1*pi);
    md = min(size(deta1));
    bigev = find(diag(deta1(1:md,1:md))>realsmall);
    veta1 = veta1(:,bigev);
    if isempty(veta1)
        eu(2) = 1;
    else
        eu(2) = norm(veta1-veta*veta'*veta1,'fro')<realsmall*n;
    end
end
F = u1(:,1:nunstab)'*inv(u1(:,nunstab+1:end)');
impact = [F*psi(nunstab+1:end,:);psi(nunstab+1:end,:)];
C = c;  	% constant term is not coded yet

author: gregkaplan

repo: phact

save_path: github-repos/MATLAB/gregkaplan-phact


text:
function [m, bl] = model_add_block(m, varargin);
% Add a block of parameters to the model.
%   [m, bl] = model_add_block(m, varargin)
%
% Return values
%   m         Updated model
%
% Arguments
%   m         Model to update
%   varargin  (key, value) pairs that can specify the following:
%   key             value
%   ---             -----
%   reg_mult        Block regularization cost
%   learn           0 => parameters are not learned; 1 => they are learned
%   lower_bounds    Lower-bound box contraints for each parameter in this block
%   shape           Dimensions of this block (e.g., [6 6 32])
%   w               Block parameter values
%   type            Block type from the block_type enumeration

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2009-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

valid_opts = {'reg_mult', 'learn', 'lower_bounds', 'shape', 'w', 'type'};
opts = getopts(varargin, valid_opts);

% OPT: w
if opts.isKey('w')
  w = opts('w');
else
  w = [];
end

% OPT: shape
if opts.isKey('shape')
  shape = opts('shape');
  if isempty(w)
    w = zeros(prod(shape), 1);
  end
else
  shape = size(w);
end

% Dimension of the block's parameter vector
dim = prod(shape);

% OPT: reg_mult
if opts.isKey('reg_mult')
  reg_mult = opts('reg_mult');
else
  reg_mult = 1;
end

% OPT: learn
if opts.isKey('learn')
  learn = opts('learn');
else
  learn = 1;
end

% OPT: lower_bounds
if opts.isKey('lower_bounds')
  lower_bounds = opts('lower_bounds');
else
  lower_bounds = -inf*ones(dim, 1);
end

% OPT: type
if opts.isKey('type')
  btype = opts('type');
else
  btype = block_types.Other;
end

% Sanity checks
assert(~isempty(w));
assert(numel(w) == dim);
assert(size(lower_bounds,1) == dim);

% Allocate new block
bl = m.numblocks + 1;
m.numblocks = bl;

% Install block
m.blocks(bl).w        = w(:);
m.blocks(bl).lb       = lower_bounds;
m.blocks(bl).learn    = learn;
m.blocks(bl).reg_mult = reg_mult;
m.blocks(bl).dim      = dim;
m.blocks(bl).shape    = shape;
m.blocks(bl).type     = btype;

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm