text:
function [elec_realigned] = ft_electroderealign(cfg, elec_original)

% FT_ELECTRODEREALIGN rotates, translates, scales and warps electrode positions. The
% default is to only rotate and translate, i.e. to do a rigid body transformation in
% which only the coordinate system is changed. With the right settings if can apply
% additional deformations to the input sensors (e.g. scale them to better fit the
% skin surface). The different methods are described in detail below.
%
% INTERACTIVE - You can display the skin surface together with the electrode or
% gradiometer positions, and manually (using the graphical user interface) adjust the
% rotation, translation and scaling parameters, so that the electrodes correspond
% with the skin.
%
% FIDUCIAL - You can apply a rigid body realignment based on three fiducial
% locations. After realigning, the fiducials in the input electrode set (typically
% nose, left and right ear) are along the same axes as the fiducials in the template
% electrode set.
%
% TEMPLATE - You can apply a spatial transformation/deformation that automatically
% minimizes the distance between the electrodes or gradiometers and a template or
% sensor array. The warping methods use a non-linear search to minimize the distance
% between the input sensor positions and the corresponding template sensors.
%
% HEADSHAPE - You can apply a spatial transformation/deformation that automatically
% minimizes the distance between the electrodes and the head surface. The warping
% methods use a non-linear search to minimize the distance between the input sensor
% positions and the projection of the electrodes on the head surface.
%
% PROJECT - This projects all electrodes to the nearest point on the
% head surface mesh.
%
% MOVEINWARD - This moves all electrodes inward according to their normals
%
% Use as
%   [elec_realigned] = ft_sensorrealign(cfg)
% with the electrode or gradiometer details in the configuration, or as
%   [elec_realigned] = ft_sensorrealign(cfg, elec_orig)
% with the electrode or gradiometer definition as 2nd input argument.
%
% The configuration can contain the following options
%   cfg.method         = string representing the method for aligning or placing the electrodes
%                        'interactive'     realign manually using a graphical user interface
%                        'fiducial'        realign using three fiducials (e.g. NAS, LPA and RPA)
%                        'template'        realign the electrodes to match a template set
%                        'headshape'       realign the electrodes to fit the head surface
%                        'project'         projects electrodes onto the head surface
%                        'moveinward'      moves electrodes inward along their normals
%   cfg.warp          = string describing the spatial transformation for the template and headshape methods
%                        'rigidbody'       apply a rigid-body warp (default)
%                        'globalrescale'   apply a rigid-body warp with global rescaling
%                        'traditional'     apply a rigid-body warp with individual axes rescaling
%                        'nonlin1'         apply a 1st order non-linear warp
%                        'nonlin2'         apply a 2nd order non-linear warp
%                        'nonlin3'         apply a 3rd order non-linear warp
%                        'nonlin4'         apply a 4th order non-linear warp
%                        'nonlin5'         apply a 5th order non-linear warp
%                        'dykstra2012'     back-project ECoG onto the cortex using energy minimzation
%                        'hermes2010'      back-project ECoG onto the cortex along the local norm vector
%                        'fsaverage'       surface-based realignment with FreeSurfer fsaverage brain (left->left or right->right)
%                        'fsaverage_sym'   surface-based realignment with FreeSurfer fsaverage_sym left hemisphere (left->left or right->left)
%                        'fsinflated'      surface-based realignment with FreeSurfer individual subject inflated brain (left->left or right->right)
%   cfg.channel        = Nx1 cell-array with selection of channels (default = 'all'),
%                        see  FT_CHANNELSELECTION for details
%   cfg.keepchannel    = string, 'yes' or 'no' (default = 'no')
%   cfg.fiducial       = cell-array with the name of three fiducials used for
%                        realigning (default = {'nasion', 'lpa', 'rpa'})
%   cfg.casesensitive  = 'yes' or 'no', determines whether string comparisons
%                        between electrode labels are case sensitive (default = 'yes')
%   cfg.feedback       = 'yes' or 'no' (default = 'no')
%
% The electrode positions can be present in the 2nd input argument or can be specified as
%   cfg.elec          = structure with electrode positions or filename, see FT_READ_SENS
%
% If you want to realign the EEG electrodes using anatomical fiducials, you should
% specify the target location of the three fiducials, e.g.
%   cfg.target.pos(1,:) = [110 0 0]     % location of the nose
%   cfg.target.pos(2,:) = [0  90 0]     % location of the left ear
%   cfg.target.pos(3,:) = [0 -90 0]     % location of the right ear
%   cfg.target.label    = {'NAS', 'LPA', 'RPA'}
%
% If you want to align EEG electrodes to a single or multiple template electrode sets
% (which will be averaged), you should specify the template electrode sets either as
% electrode structures (i.e. when they are already read in memory) or their file
% names using
%   cfg.target          = single electrode set that serves as standard
% or
%   cfg.target{1..N}    = list of electrode sets that will be averaged
%
% If you want to align EEG electrodes to the head surface, you should specify the head surface as
%   cfg.headshape      = a filename containing headshape, a structure containing a
%                        single triangulated boundary, or a Nx3 matrix with surface
%                        points
%
% If you want to align ECoG electrodes to the pial surface, you first need to compute
% the cortex hull with FT_PREPARE_MESH. Then use either the algorithm described in
% Dykstra et al. (2012, Neuroimage) or in Hermes et al. (2010, J Neurosci methods) to
% snap the electrodes back to the cortical hull, e.g.
%   cfg.method         = 'headshape'
%   cfg.warp           = 'dykstra2012', or 'hermes2010'
%   cfg.headshape      = a filename containing headshape, a structure containing a
%                        single triangulated boundary, or a Nx3 matrix with surface
%                        points
%   cfg.feedback       = 'yes' or 'no' (default), feedback of the iteration procedure
%
% Additional configuration options for cfg.warp = 'dykstra2012'
%   cfg.maxiter        = number (default: 50), maximum number of optimization iterations
%   cfg.pairmethod     = 'pos' (default) or 'label', the method for electrode
%                        pairing on which the deformation energy is based
%   cfg.isodistance    = 'yes', 'no' (default) or number, to enforce isotropic
%                        inter-electrode distances (pairmethod 'label' only)
%   cfg.deformweight   = number (default: 1), weight of deformation relative
%                        to shift energy cost (lower increases grid flexibility)
%
% If you want to move the electrodes inward, you should specify
%   cfg.moveinward     = number, the distance that the electrode should be moved
%                        inward (negative numbers result in an outward move)
%
% If you want to align ECoG electrodes to the freesurfer average brain, you should
% specify the path to your headshape (e.g., lh.pial), and ensure you have the
% corresponding registration file (e.g., lh.sphere.reg) in the same directory.
% Moreover, the path to the local freesurfer home is required. Note that, because the
% electrodes are being aligned to the fsaverage brain, the corresponding brain should
% be also used when plotting the data, i.e. use freesurfer/subjects/fsaverage/surf/lh.pial
% rather than surface_pial_left.mat
%   cfg.method         = 'headshape'
%   cfg.warp           = 'fsaverage'
%   cfg.headshape      = string, filename containing subject headshape (e.g. <path to freesurfer/surf/lh.pial>)
%   cfg.fshome         = string, path to freesurfer
%
% See also FT_READ_SENS, FT_VOLUMEREALIGN, FT_INTERACTIVEREALIGN,
% FT_DETERMINE_COORDSYS, FT_PREPARE_MESH

% Copyright (C) 2005-2019, Robert Oostenveld, Arjen Stolk
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% the interactive method uses a global variable to get the data from the figure when it is closed
global norm

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar    elec_original
ft_preamble provenance elec_original
ft_preamble trackconfig

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed',    {'template', 'target'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'realignfiducials', 'fiducial'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'realignfiducial',  'fiducial'});
cfg = ft_checkconfig(cfg, 'renamedval', {'warp', 'homogenous', 'rigidbody'});
cfg = ft_checkconfig(cfg, 'renamedval', {'warp', 'homogeneous', 'rigidbody'});
cfg = ft_checkconfig(cfg, 'forbidden', 'outline');

% set the defaults
cfg.warp          = ft_getopt(cfg, 'warp', 'rigidbody');
cfg.channel       = ft_getopt(cfg, 'channel',  'all');
cfg.keepchannel   = ft_getopt(cfg, 'keepchannel', 'no');
cfg.feedback      = ft_getopt(cfg, 'feedback', 'no');
cfg.casesensitive = ft_getopt(cfg, 'casesensitive', 'no');
cfg.headshape     = ft_getopt(cfg, 'headshape', []);     % for triangulated head surface, without labels
cfg.target        = ft_getopt(cfg, 'target',  []);       % for electrodes or fiducials, always with labels
cfg.coordsys      = ft_getopt(cfg, 'coordsys');          % this allows for automatic template fiducial placement

if isempty(cfg.target)
  % remove the field, otherwise ft_checkconfig will complain
  cfg = rmfield(cfg, 'target');
end

if ~isempty(cfg.coordsys) && isempty(cfg.target)
  % set the template fiducial locations according to the coordinate system
  switch lower(cfg.coordsys)
    case 'ctf'
      cfg.target = [];
      cfg.target.coordsys = 'ctf';
      cfg.target.pos(1,:) = [100  0 0];
      cfg.target.pos(2,:) = [0   80 0];
      cfg.target.pos(3,:) = [0  -80 0];
      cfg.target.label{1} = 'NAS';
      cfg.target.label{2} = 'LPA';
      cfg.target.label{3} = 'RPA';
    otherwise
      ft_error('the %s coordinate system is not automatically supported, please specify fiducial details in cfg.target')
  end
end

% ensure that the right cfg options have been set corresponding to the method
switch cfg.method
  case 'template'        % realign the sensors to match a template set
    cfg = ft_checkconfig(cfg, 'required', 'target', 'forbidden', 'headshape');
  case 'headshape'     % realign the sensors to fit the head surface
    cfg = ft_checkconfig(cfg, 'required', 'headshape', 'forbidden', 'target');
  case 'fiducial'        % realign using the NAS, LPA and RPA fiducials
    cfg = ft_checkconfig(cfg, 'required', 'target', 'forbidden', 'headshape');
  case 'moveinward'      % moves eletrodes inward
    cfg = ft_checkconfig(cfg, 'required', 'moveinward');
end % switch cfg.method

if strcmp(cfg.method, 'fiducial') && isfield(cfg, 'warp') && ~isequal(cfg.warp, 'rigidbody')
  ft_warning('The method ''fiducial'' implies a rigid body tramsformation. See also http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1722');
  cfg.warp = 'rigidbody';
end

if isfield(cfg, 'headshape') && isa(cfg.headshape, 'config')
  % convert the nested config-object back into a normal structure
  cfg.headshape = struct(cfg.headshape);
end

if isfield(cfg, 'target') && isa(cfg.target, 'config')
  % convert the nested config-object back into a normal structure
  cfg.target = struct(cfg.target);
end

% the data can be passed as input arguments or can be read from disk
hasdata = exist('elec_original', 'var');

% get the electrode definition that should be warped
if ~hasdata
  elec_original = ft_fetch_sens(cfg);
else
  % the input electrodes were specified as second input argument
  % or read from cfg.inputfile
end

% ensure that the units are specified
elec_original = ft_determine_units(elec_original);

% ensure up-to-date sensor description (Oct 2011)
elec_original = ft_datatype_sens(elec_original);

% ensure that channel and electrode positions are the same
assert(isequaln(elec_original.elecpos, elec_original.chanpos), 'this function requires same electrode and channel positions');

% remember the original electrode locations and labels and do all the work with a
% temporary copy, this involves channel selection and changing to lower case
elec = elec_original;

% instead of working with all sensors, only work with the fiducials
% this is useful for gradiometer structures
if strcmp(cfg.method, 'fiducial') && isfield(elec, 'fid')
  fprintf('using the fiducials instead of the sensor positions\n');
  elec.fid.unit = elec.unit;
  elec          = elec.fid;
end

usetarget    = isfield(cfg, 'target')    && ~isempty(cfg.target);
useheadshape = isfield(cfg, 'headshape') && ~isempty(cfg.headshape);

if usetarget
  % get the template electrode definitions
  if ~iscell(cfg.target)
    cfg.target = {cfg.target};
  end
  Ntemplate = length(cfg.target);
  for i=1:Ntemplate
    if isstruct(cfg.target{i})
      target(i) = cfg.target{i};
    else
      target(i) = ft_read_sens(cfg.target{i}, 'senstype', 'eeg');
    end
  end
  clear tmp
  for i=1:Ntemplate
    % ensure up-to-date sensor description
    % ensure that the units are consistent with the electrodes
    tmp(i) = ft_convert_units(ft_datatype_sens(target(i)), elec.unit);
  end
  target = tmp;
end

if useheadshape
  % get the surface describing the head shape
  [headshape.pos, headshape.tri] = headsurface([], [], 'headshape', cfg.headshape);

  % ensure that the units are consistent with the electrodes
  headshape = ft_convert_units(headshape, elec.unit);
end

% convert all labels to lower case for string comparisons
cfg.channel = ft_channelselection(cfg.channel, elec.label);
if strcmp(cfg.casesensitive, 'no')
  elec.label  = lower(elec.label);
  cfg.channel = lower(cfg.channel);
  if usetarget
    for j=1:length(target)
      for i=1:length(target(j).label)
        target(j).label{i} = lower(target(j).label{i});
      end
    end
  end
end
[cfgsel, datsel] = match_str(cfg.channel, elec.label);
% keep the original channel labels
label_original = elec_original.label(datsel);

% start with an empty structure, this will be returned at the end
norm = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(cfg.method, 'template')
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % determine electrode selection and overlapping subset for warping
  cfg.channel = ft_channelselection(cfg.channel, elec.label);
  for i=1:Ntemplate
    cfg.channel = ft_channelselection(cfg.channel, target(i).label);
  end

  % make consistent subselection of electrodes
  [cfgsel, datsel] = match_str(cfg.channel, elec.label);
  elec.label = elec.label(datsel);
  elec.elecpos   = elec.elecpos(datsel,:);
  for i=1:Ntemplate
    [cfgsel, datsel] = match_str(cfg.channel, target(i).label);
    target(i).label   = target(i).label(datsel);
    target(i).elecpos = target(i).elecpos(datsel,:);
  end

  % compute the average of the target electrode positions
  average = ft_average_sens(target);

  fprintf('warping electrodes to average template... '); % the newline comes later
  [norm.elecpos, norm.m] = ft_warp_optim(elec.elecpos, average.elecpos, cfg.warp);
  norm.label = elec.label;

  dpre  = mean(sqrt(sum((average.elecpos - elec.elecpos).^2, 2)));
  dpost = mean(sqrt(sum((average.elecpos - norm.elecpos).^2, 2)));
  fprintf('mean distance prior to warping %f, after warping %f\n', dpre, dpost);

  if strcmp(cfg.feedback, 'yes')
    % create an empty figure, continued below...
    figure
    axis equal
    axis vis3d
    hold on
    xlabel('x')
    ylabel('y')
    zlabel('z')

    % plot all electrodes before warping
    ft_plot_sens(elec, 'r*');

    % plot all electrodes after warping
    ft_plot_sens(norm, 'm.', 'label', 'label');

    % plot the template electrode locations
    ft_plot_sens(average, 'b.');

    % plot lines connecting the input and the realigned electrode locations with the template locations
    my_line3(elec.elecpos, average.elecpos, 'color', 'r');
    my_line3(norm.elecpos, average.elecpos, 'color', 'm');
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(cfg.method, 'headshape')
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % determine electrode selection and overlapping subset for warping
  cfg.channel = ft_channelselection(cfg.channel, elec.label);
  [cfgsel, datsel] = match_str(cfg.channel, elec.label);
  elec.label   = elec.label(datsel);
  elec.elecpos = elec.elecpos(datsel,:);

  norm.label = elec.label;
  if strcmp(cfg.warp, 'dykstra2012')
    norm.elecpos = warp_dykstra2012(cfg, elec, headshape);
  elseif strcmp(cfg.warp, 'hermes2010')
    norm.elecpos = warp_hermes2010(cfg, elec, headshape);
  elseif strcmp(cfg.warp, 'fsaverage')
    norm.elecpos = warp_fsaverage(cfg, elec);
  elseif strcmp(cfg.warp, 'fsaverage_sym')
    norm.elecpos = warp_fsaverage_sym(cfg, elec);
  elseif strcmp(cfg.warp, 'fsinflated')
    norm.elecpos = warp_fsinflated(cfg, elec);
  else
    fprintf('warping electrodes to skin surface... '); % the newline comes later
    [norm.elecpos, norm.m] = ft_warp_optim(elec.elecpos, headshape, cfg.warp);

    dpre  = ft_warp_error([],     elec.elecpos, headshape, cfg.warp);
    dpost = ft_warp_error(norm.m, elec.elecpos, headshape, cfg.warp);
    fprintf('mean distance prior to warping %f, after warping %f\n', dpre, dpost);
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(cfg.method, 'fiducial')
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % the fiducials have to be present in the electrodes and in the template set
  label = intersect(lower(elec.label), lower(target.label));

  if ~isfield(cfg, 'fiducial') || isempty(cfg.fiducial)
    % try to determine the names of the fiducials automatically
    option1 = {'nasion' 'left' 'right'};
    option2 = {'nasion' 'lpa' 'rpa'};
    option3 = {'nz' 'left' 'right'};
    option4 = {'nz' 'lpa' 'rpa'};
    option5 = {'nas' 'left' 'right'};
    option6 = {'nas' 'lpa' 'rpa'};
    if length(match_str(label, option1))==3
      cfg.fiducial = option1;
    elseif length(match_str(label, option2))==3
      cfg.fiducial = option2;
    elseif length(match_str(label, option3))==3
      cfg.fiducial = option3;
    elseif length(match_str(label, option4))==3
      cfg.fiducial = option4;
    elseif length(match_str(label, option5))==3
      cfg.fiducial = option5;
    elseif length(match_str(label, option6))==3
      cfg.fiducial = option6;
    else
      ft_error('could not determine consistent fiducials in the input and the target, please specify cfg.fiducial or cfg.coordsys')
    end
  end
  fprintf('matching fiducials {''%s'', ''%s'', ''%s''}\n', cfg.fiducial{1}, cfg.fiducial{2}, cfg.fiducial{3});

  % determine electrode selection
  cfg.channel = ft_channelselection(cfg.channel, elec.label);
  [cfgsel, datsel] = match_str(cfg.channel, elec.label);
  elec.label     = elec.label(datsel);
  elec.elecpos   = elec.elecpos(datsel,:);

  if length(cfg.fiducial)~=3
    ft_error('you must specify exactly three fiducials');
  end

  % do case-insensitive search for fiducial locations
  nas_indx = match_str(lower(elec.label), lower(cfg.fiducial{1}));
  lpa_indx = match_str(lower(elec.label), lower(cfg.fiducial{2}));
  rpa_indx = match_str(lower(elec.label), lower(cfg.fiducial{3}));
  if length(nas_indx)~=1 || length(lpa_indx)~=1 || length(rpa_indx)~=1
    ft_error('not all fiducials were found in the electrode set');
  end
  elec_nas = elec.elecpos(nas_indx,:);
  elec_lpa = elec.elecpos(lpa_indx,:);
  elec_rpa = elec.elecpos(rpa_indx,:);

  % FIXME change the flow in the remainder
  % if one or more template electrode sets are specified, then align to the average of those
  % if no template is specified, then align so that the fiducials are along the axis

  % find the matching fiducials in the template and average them
  tmpl_nas = nan(Ntemplate,3);
  tmpl_lpa = nan(Ntemplate,3);
  tmpl_rpa = nan(Ntemplate,3);
  for i=1:Ntemplate
    nas_indx = match_str(lower(target(i).label), lower(cfg.fiducial{1}));
    lpa_indx = match_str(lower(target(i).label), lower(cfg.fiducial{2}));
    rpa_indx = match_str(lower(target(i).label), lower(cfg.fiducial{3}));
    if length(nas_indx)~=1 || length(lpa_indx)~=1 || length(rpa_indx)~=1
      ft_error('not all fiducials were found in template %d', i);
    end
    tmpl_nas(i,:) = target(i).elecpos(nas_indx,:);
    tmpl_lpa(i,:) = target(i).elecpos(lpa_indx,:);
    tmpl_rpa(i,:) = target(i).elecpos(rpa_indx,:);
  end
  tmpl_nas = mean(tmpl_nas,1);
  tmpl_lpa = mean(tmpl_lpa,1);
  tmpl_rpa = mean(tmpl_rpa,1);

  % realign both to a common coordinate system
  elec2common  = ft_headcoordinates(elec_nas, elec_lpa, elec_rpa);
  templ2common = ft_headcoordinates(tmpl_nas, tmpl_lpa, tmpl_rpa);

  % compute the combined transform
  norm         = [];
  norm.m       = templ2common \ elec2common;

  % apply the transformation to the fiducials as sanity check
  norm.elecpos(1,:) = ft_warp_apply(norm.m, elec_nas, 'homogeneous');
  norm.elecpos(2,:) = ft_warp_apply(norm.m, elec_lpa, 'homogeneous');
  norm.elecpos(3,:) = ft_warp_apply(norm.m, elec_rpa, 'homogeneous');
  norm.label        = cfg.fiducial;

  nas_indx = match_str(lower(elec.label), lower(cfg.fiducial{1}));
  lpa_indx = match_str(lower(elec.label), lower(cfg.fiducial{2}));
  rpa_indx = match_str(lower(elec.label), lower(cfg.fiducial{3}));
  dpre  = mean(sqrt(sum((elec.elecpos([nas_indx lpa_indx rpa_indx],:) - [tmpl_nas; tmpl_lpa; tmpl_rpa]).^2, 2)));
  nas_indx = match_str(lower(norm.label), lower(cfg.fiducial{1}));
  lpa_indx = match_str(lower(norm.label), lower(cfg.fiducial{2}));
  rpa_indx = match_str(lower(norm.label), lower(cfg.fiducial{3}));
  dpost = mean(sqrt(sum((norm.elecpos([nas_indx lpa_indx rpa_indx],:) - [tmpl_nas; tmpl_lpa; tmpl_rpa]).^2, 2)));
  fprintf('mean distance between fiducials prior to realignment %f, after realignment %f\n', dpre, dpost);

  if strcmp(cfg.feedback, 'yes')
    % create an empty figure, continued below...
    figure
    axis equal
    axis vis3d
    hold on
    xlabel('x')
    ylabel('y')
    zlabel('z')

    % plot the first three electrodes before transformation
    my_plot3(elec.elecpos(1,:), 'r*');
    my_plot3(elec.elecpos(2,:), 'r*');
    my_plot3(elec.elecpos(3,:), 'r*');
    my_text3(elec.elecpos(1,:), elec.label{1}, 'color', 'r');
    my_text3(elec.elecpos(2,:), elec.label{2}, 'color', 'r');
    my_text3(elec.elecpos(3,:), elec.label{3}, 'color', 'r');

    % plot the template fiducials
    my_plot3(tmpl_nas, 'b*');
    my_plot3(tmpl_lpa, 'b*');
    my_plot3(tmpl_rpa, 'b*');
    my_text3(tmpl_nas, ' nas', 'color', 'b');
    my_text3(tmpl_lpa, ' lpa', 'color', 'b');
    my_text3(tmpl_rpa, ' rpa', 'color', 'b');

    % plot all electrodes after transformation
    my_plot3(norm.elecpos, 'm.');
    my_plot3(norm.elecpos(1,:), 'm*');
    my_plot3(norm.elecpos(2,:), 'm*');
    my_plot3(norm.elecpos(3,:), 'm*');
    my_text3(norm.elecpos(1,:), norm.label{1}, 'color', 'm');
    my_text3(norm.elecpos(2,:), norm.label{2}, 'color', 'm');
    my_text3(norm.elecpos(3,:), norm.label{3}, 'color', 'm');
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(cfg.method, 'interactive')
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  tmpcfg = [];
  tmpcfg.individual.elec = elec;
  if isfield(cfg, 'headshape') && ~isempty(cfg.headshape)
    tmpcfg.template.headshape = cfg.headshape;
  end
  if isfield(cfg, 'target') && ~isempty(cfg.target)
    if iscell(cfg.target)
      if numel(cfg.target)>1
        ft_notice('computing the average electrode positions');
        tmpcfg.template.elec = ft_average_sens(cfg.target);
      else
        tmpcfg.template.elec = cfg.target{1};
      end
    elseif isstruct(cfg.target)
      tmpcfg.template.elec = cfg.target;
    end
    tmpcfg.template.elecstyle = {'facecolor', 'blue'};
    ft_info('plotting the target electrodes in blue');
  end

  % use the more generic ft_interactiverealign for the actual work
  tmpcfg = ft_interactiverealign(tmpcfg);
  % only keep the transformation, it will be applied to the electrodes further down
  norm.m = tmpcfg.m;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(cfg.method, 'project')
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % determine electrode selection
  cfg.channel = ft_channelselection(cfg.channel, elec.label);
  [cfgsel, datsel] = match_str(cfg.channel, elec.label);
  elec.label     = elec.label(datsel);
  elec.elecpos   = elec.elecpos(datsel,:);

  norm.label = elec.label;
  [dum, norm.elecpos] = project_elec(elec.elecpos, headshape.pos, headshape.tri);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(cfg.method, 'moveinward')
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % determine electrode selection
  cfg.channel = ft_channelselection(cfg.channel, elec.label);
  [cfgsel, datsel] = match_str(cfg.channel, elec.label);
  elec.label     = elec.label(datsel);
  elec.elecpos   = elec.elecpos(datsel,:);

  norm.label = elec.label;
  norm.elecpos = moveinward(elec.elecpos, cfg.moveinward);

else
  ft_error('unknown method');
end % if method

% apply the spatial transformation to all electrodes, and replace the
% electrode labels by their case-sensitive original values
switch cfg.method
  case {'template', 'headshape'}
    if strcmpi(cfg.warp, 'dykstra2012') || strcmpi(cfg.warp, 'hermes2010') || ...
        strcmpi(cfg.warp, 'fsaverage') || strcmpi(cfg.warp, 'fsaverage_sym') || strcmpi(cfg.warp, 'fsinflated')
      elec_realigned = norm;
      elec_realigned.label = label_original;
    else
      % the transformation is a linear or non-linear warp, i.e. a vector
      try
        % convert the vector with fitted parameters into a 4x4 homogenous transformation
        % apply the transformation to the original complete set of sensors
        elec_realigned = ft_transform_geometry(feval(cfg.warp, norm.m), elec_original);
      catch
        % the previous section will fail for nonlinear transformations
        elec_realigned.label = label_original;
        try
          elec_realigned.elecpos = ft_warp_apply(norm.m, elec_original.elecpos, cfg.warp);
        end % FIXME why is an error here not dealt with?
      end
      % remember the transformation
      elec_realigned.(cfg.warp) = norm.m;
    end

  case  {'fiducial' 'interactive'}
    % the transformation is a 4x4 homogenous matrix
    % apply the transformation to the original complete set of sensors
    elec_realigned = ft_transform_geometry(norm.m, elec_original);
    % remember the transformation
    elec_realigned.homogeneous = norm.m;

  case {'project', 'moveinward'}
    % nothing to be done
    elec_realigned = norm;
    elec_realigned.label = label_original;

  otherwise
    ft_error('unknown method');
end

% the coordinate system is in general not defined after transformation
if isfield(elec_realigned, 'coordsys')
  elec_realigned = rmfield(elec_realigned, 'coordsys');
end

% in some cases the coordinate system matches that of the input target or headshape
switch cfg.method
  case 'template'
    if isfield(target, 'coordsys')
      elec_realigned.coordsys = target.coordsys;
    end
  case 'headshape'
    if isfield(headshape, 'coordsys')
      elec_realigned.coordsys = headshape.coordsys;
    end
    if isfield(elec_original, 'coordsys')
      if strcmp(cfg.warp, 'dykstra2012') || strcmp(cfg.warp, 'hermes2010')  % this warp simply moves the electrodes in the same coordinate space
        elec_realigned.coordsys = elec_original.coordsys;
      elseif strcmp(cfg.warp, 'fsaverage')
        elec_realigned.coordsys = 'fsaverage';
      elseif strcmp(cfg.warp, 'fsaverage_sym')
        elec_realigned.coordsys = 'fsaverage_sym';
      end
    end
  case 'fiducial'
    if isfield(target, 'coordsys')
      elec_realigned.coordsys = target.coordsys;
    end
  case 'interactive'
    % the coordinate system is not known
  case {'project', 'moveinward'}
    % the coordinate system remains the same
    if isfield(elec_original, 'coordsys')
      elec_realigned.coordsys = elec_original.coordsys;
    end
  otherwise
    ft_error('unknown method');
end

if istrue(cfg.keepchannel)
  % append the channels that are not realigned
  [dum, idx] = setdiff(elec_original.label, elec_realigned.label);
  idx = sort(idx);
  elec_realigned.label = [elec_realigned.label; elec_original.label(idx)];
  elec_realigned.elecpos = [elec_realigned.elecpos; elec_original.elecpos(idx,:)];
end

% channel positions are identical to the electrode positions (this was checked at the start)
elec_realigned.chanpos = elec_realigned.elecpos;
elec_realigned.tra = eye(numel(elec_realigned.label));

% copy over unit, chantype, chanunit, and tra information in case this was not already done
if ~isfield(elec_realigned, 'unit') && isfield(elec_original, 'unit')
  elec_realigned.unit = elec_original.unit;
end
if ~isfield(elec_realigned, 'chantype') && isfield(elec_original, 'chantype')
  idx = match_str(elec_original.label, elec_realigned.label);
  elec_realigned.chantype = elec_original.chantype(idx);
end
if ~isfield(elec_realigned, 'chanunit') && isfield(elec_original, 'chanunit')
  elec_realigned.chanunit = elec_original.chanunit;
  idx = match_str(elec_original.label, elec_realigned.label);
  elec_realigned.chanunit = elec_original.chanunit(idx);
end

% update it to the latest version
elec_realigned = ft_datatype_sens(elec_realigned);

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble previous   elec_original
ft_postamble provenance elec_realigned
ft_postamble history    elec_realigned

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% some simple SUBFUNCTIONs that facilitate 3D plotting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = my_plot3(xyz, varargin)
h = plot3(xyz(:,1), xyz(:,2), xyz(:,3), varargin{:});
function h = my_text3(xyz, varargin)
h = text(xyz(:,1), xyz(:,2), xyz(:,3), varargin{:});
function my_line3(xyzB, xyzE, varargin)
for i=1:size(xyzB,1)
  line([xyzB(i,1) xyzE(i,1)], [xyzB(i,2) xyzE(i,2)], [xyzB(i,3) xyzE(i,3)], varargin{:})
end

author: spm


text:
% VL_DEMO_SIFT_MATCH  Demo: SIFT: basic matching

pfx = fullfile(vl_root,'figures','demo') ;
randn('state',0) ;
rand('state',0) ;
figure(1) ; clf ;

% --------------------------------------------------------------------
%                                                    Create image pair
% --------------------------------------------------------------------

Ia = imread(fullfile(vl_root,'data','roofs1.jpg')) ;
Ib = imread(fullfile(vl_root,'data','roofs2.jpg')) ;

% --------------------------------------------------------------------
%                                           Extract features and match
% --------------------------------------------------------------------

[fa,da] = vl_sift(im2single(rgb2gray(Ia))) ;
[fb,db] = vl_sift(im2single(rgb2gray(Ib))) ;

[matches, scores] = vl_ubcmatch(da,db) ;

[drop, perm] = sort(scores, 'descend') ;
matches = matches(:, perm) ;
scores  = scores(perm) ;

figure(1) ; clf ;
imagesc(cat(2, Ia, Ib)) ;
axis image off ;
vl_demo_print('sift_match_1', 1) ;

figure(2) ; clf ;
imagesc(cat(2, Ia, Ib)) ;

xa = fa(1,matches(1,:)) ;
xb = fb(1,matches(2,:)) + size(Ia,2) ;
ya = fa(2,matches(1,:)) ;
yb = fb(2,matches(2,:)) ;

hold on ;
h = line([xa ; xb], [ya ; yb]) ;
set(h,'linewidth', 1, 'color', 'b') ;

vl_plotframe(fa(:,matches(1,:))) ;
fb(1,:) = fb(1,:) + size(Ia,2) ;
vl_plotframe(fb(:,matches(2,:))) ;
axis image off ;

vl_demo_print('sift_match_2', 1) ;

author: Cloud-CV

repo: object-proposals


text:
% LTFAT - Signals
%
%  Signal generators
%    CTESTFUN       -  Complex valued test function.
%    NOISE          -  Stochastic noise generator.
%    PINKNOISE      -  Pink noise.
%    EXPCHIRP       -  Exponential chirp.
%
%  Sound signals.
%    BAT            -  Bat chirp.
%    BATMASK        -  Mask for bat signal.
%    GREASY         -  Woman speaking the word 'greasy'
%    COCKTAILPARTY  -  Sentence by male, native English speaker.
%    GSPI           -  Glockenspiel test signal.
%    LINUS          -  Linus pronouncing Linux.
%    LTFATLOGO      -  Synthetic sound from spectrogram reconstruction.
%    OTOCLICK       -  Click-evoked otoacoustic emmision.
%    TRAINDOPPLER   -  Sound of passing train.
%
%  Images.
%    CAMERAMAN      -  Greyscale image of the cameraman.
%    LICHTENSTEIN   -  Color image of the Lichtenstein castle.
%    LTFATTEXT      -  Black and white word: 'LTFAT'.
%
%  For help, bug reports, suggestions etc. please visit
%  http://github.com/ltfat/ltfat/issues
%

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
classdef CelestialBodyIntegration
    %CelestialBodyIntegration Summary of this class goes here
    %   Detailed explanation goes here

    properties
        celBodyData CelestialBodyData
    end

    methods
        function obj = CelestialBodyIntegration(celBodyData)
            obj.celBodyData = celBodyData;
        end

        function integrateCelestialBodies(obj, minUT, maxUT, hWaitbar)
            arguments
                obj(1,1) CelestialBodyIntegration
                minUT(1,1) double
                maxUT(1,1) double
                hWaitbar = [];
            end

%             allBodies = obj.celBodyData.getAllBodyInfo();
            [~, allBodiesCell] = ma_getSortedBodyNames(obj.celBodyData);

            for(i=1:length(allBodiesCell)) %#ok<*NO4LP>
                allBodies(i) = allBodiesCell{i}; %#ok<AGROW>
            end

            bool = [allBodies.propTypeEnum] == BodyPropagationTypeEnum.Numerical;
            if(any(bool))
%                 hFig = uifigure('Position',[0,0,400,125], 'WindowStyle','modal', 'Icon','logoSquare_48px_transparentBg.png');
%                 centerUIFigure(hFig);
                drawnow;

                numericBodies = allBodies(bool);
                numericInds = 1:sum(bool);

                analyticBodies = allBodies(not(bool));
                analyticInds = 1:sum(not(bool));

                epochs = [numericBodies.epoch];
                if(numel(unique(epochs)) == 1)
                    t0 = unique(epochs);

                    [y0, numericMasses] = CelestialBodyIntegration.getInitialStates(numericBodies);
                    if(t0 > minUT)
                        [T,Y,~] = obj.callIntegrator(numericBodies, numericMasses, numericInds, analyticBodies, analyticInds, t0, minUT, y0, hWaitbar);

                        t0 = T(end);
                        y0 = Y(end,:);
                    end

                    [T,~,YY] = obj.callIntegrator(numericBodies, numericMasses, numericInds, analyticBodies, analyticInds, t0, maxUT, y0, hWaitbar);

                    if(not(isempty(hWaitbar)))
                        hWaitbar.Message = 'Performing frame rotations on computed states.  Please wait.';
                        hWaitbar.Value = 0;
                        hWaitbar.Indeterminate = 'off';
                    end

                    for(i=1:size(YY,3))
                        bodyInfo = numericBodies(i);
                        states = YY(:,:,i);

                        times = T;
                        rVects = states(:,1:3)';
                        vVects = states(:,4:6)';
                        cartElem = CartesianElementSet(times, rVects, vVects, obj.celBodyData.globalBaseFrame, true);

                        if(not(isempty(bodyInfo.getParBodyInfo(obj.celBodyData))))
                            parentBody = bodyInfo.getParBodyInfo(obj.celBodyData);
                            convertFrame = parentBody.getBodyCenteredInertialFrame();

                        else
                            convertFrame = obj.celBodyData.globalBaseFrame;
                        end

                        cartElem = convertToFrame(cartElem, convertFrame);

                        times = [cartElem.time]';
                        rVects = [cartElem.rVect];
                        vVects = [cartElem.vVect];
                        states = vertcat(rVects, vVects)';
                        bodyInfo.setStateCacheData(times, states, convertFrame);

                        if(not(isempty(hWaitbar)))
                            hWaitbar.Value = i/size(YY,3);
                        end
                    end
                else
                    error('All bodies must be defined with the same orbit epoch.');
                end

%                 if(isvalid(hFig))
%                     close(hFig);
%                 end
            end
        end
    end

    methods(Access=private)
        function [T,Y,YY] = callIntegrator(obj, numericBodies, numericMasses, numericInds, analyticBodies, analyticInds, t0, tf, y0, hWaitbar)
            numPts = length(numericBodies);

            ptCombsNum = getPtsCombs(numericInds, numericInds, true);
            ptCombsAnalytic = getPtsCombs(numericInds, analyticInds, false);

            massesAnalytic = [analyticBodies.gm];

            hWaitbar.Indeterminate = false;

            odefcn = @(t,y) odefun(t,y, numericMasses, numPts, ptCombsNum, massesAnalytic, ptCombsAnalytic, analyticBodies, obj.celBodyData);
            options = odeset('AbsTol',1E-10, 'RelTol',1E-6, ...
                             'OutputFcn',@(t,y,flag) outputFcn(t,y,flag, t0, tf, hWaitbar)); %'Events',@(t,y) myEventsFcn(t,y, numPts, masses, density, ptCombs),
                                                                             %'OutputFcn',@(t,y,flag) myOutputFcn(t,y,flag), ...

            tspan = [t0 tf];
            [T,Y] = ode113(odefcn, tspan, y0, options);

            YY = reshape(Y, size(Y,1), size(Y,2)/numPts, numPts);
        end
    end

    methods(Static, Access=private)
        function [y0, masses] = getInitialStates(bodies)
            rVects = NaN(3, length(bodies));
            vVects = rVects;
            masses = NaN(1, length(bodies));
            for(i=1:length(bodies))
                bodyInfo = bodies(i);

                [rVectB, vVectB] = getPositOfBodyWRTSun(bodyInfo.epoch, bodyInfo, bodyInfo.celBodyData);
                rVects(:,i) = rVectB(:);
                vVects(:,i) = vVectB(:);
                masses(i) = bodyInfo.gm;
            end

            states = [rVects;vVects];
            y0 = states(:);
        end
    end
end

function ydot = odefun(t,y, masses, numPtsNumeric, ptCombsNumeric, massesAnalytic, ptCombsAnalytic, analyticBodies, celBodyData)
    %Gravity from numerically propagated bodies
    states = reshape(y, 6, numPtsNumeric);
    rVects = states(1:3,:);
    vVects = states(4:6,:);

    ydot = NaN(size(states));
    ydot(1:3,:) = vVects;

    rDiffs = rVects(:, ptCombsNumeric(2,:)) - rVects(:, ptCombsNumeric(1,:));
    rDiffMags = sqrt(sum((rDiffs.^2), 1));

    massProducts = masses(ptCombsNumeric(2,:));

    fVectMag = (massProducts ./ (rDiffMags.^3));
    accelVects = bsxfun(@times, fVectMag, rDiffs);

    accelVects = reshape(accelVects, [3, size(accelVects,2)/numPtsNumeric, numPtsNumeric]);
    accelVectsTotalNumeric = sum(accelVects,2);
    accelVectsTotalNumeric = reshape(accelVectsTotalNumeric,[3,numPtsNumeric]);

    %Gravity from analytic propagated bodies
    if(not(isempty(massesAnalytic)))
        rVectAnalytic = NaN(3,numel(analyticBodies));
        for(i=1:length(analyticBodies))
            rVectAnalytic(:,i) = getPositOfBodyWRTSun(t, analyticBodies(i), celBodyData);
        end

        rDiffs = -1*(rVects(:, ptCombsAnalytic(2,:)) - rVectAnalytic(:, ptCombsAnalytic(1,:)));
        rDiffMags = sqrt(sum((rDiffs.^2), 1));

        massProducts = massesAnalytic(ptCombsAnalytic(1,:));

        fVectMag = (massProducts ./ (rDiffMags.^3));
        accelVects = bsxfun(@times, fVectMag, rDiffs);

        accelVectsTotalAnalytic = NaN([3, numPtsNumeric]);
        for(i=1:numPtsNumeric)
            accelVectsTotalAnalytic(:,i) = sum(accelVects(:, ptCombsAnalytic(2,:) == i), 2);
        end
    else
        accelVectsTotalAnalytic = zeros(size(accelVectsTotalNumeric));
    end

    %Sum up inputs
    accelVectsTotal = accelVectsTotalNumeric + accelVectsTotalAnalytic;
    ydot(4:6,:) = accelVectsTotal; %a = F/m
    ydot = ydot(:);
end

%% Output function (for fun)
function status = outputFcn(t,~,flag, t0, tf, hWaitbar)
    if(isempty(flag))
        tPercent = (t-t0)/(tf-t0);
        hWaitbar.Value = tPercent;
    end
    status = 0;
end

%% Utility function(s)
function ptCombs = getPtsCombs(ptInds1, ptInd2, removeDups)
    ptCombs = combvec(ptInds1,ptInd2);

    if(removeDups)
        ptCombs(:, ptCombs(1,:) == ptCombs(2,:)) = []; %eliminate duplicates
    end

    A = ptCombs(1,:);
    B = ptCombs(2,:);
    ptCombs = [B;A];
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function val = peval(x,coef,deg)
% function val = peval(x,coef,deg)
%
% DESCRIPTION
%   Evaluate a polynomial at x.
%
%   Note: The SUBS function should be used for polynomial evaluations.
%   The PEVAL function has no error checking for speed and it assumes
%   that the polynomial is a row or column vector.  It also assumes that
%   the rows of x are ordered to correspond to the columns in the degree
%   matrix deg.
%
% INPUTS
%   x: nvars-by-npts matrix each column of which specifies a point
%        at which to evaluate the polynomial.
%   coef: nterms-by-lp coefficient matrix where lp is the length
%        of the polynomial
%   deg: nterms-by-nvars degree matrix
%
% OUTPUTS
%   val: lp-by-npts matrix where each column specifies the value of
%        the polynomial evaluated at the corresponding column of x.
%
% SYNTAX
%   val = peval(x,coef,deg)

% 6/15/06  PJS  Initial Coding

% Compute dimensions
[nterms,nvars]=size(deg);
npts = size(x,2);
lp = size(coef,2);
coef = full(coef);
deg = full(deg);

% Compute monomials
xs = shiftdim(x,-1);            % xs is 1-by-nvars-by-npts
xrep = repmat(xs,[nterms 1 1]); % xrep is nterms-by-nvars-by-npts
drep = repmat(deg,[1 1 npts]);  % drep is nterms-by-nvars-by-npts

% npow = xrep;
% idx = find(drep~=1);
% npow(idx) = xrep(idx).^drep(idx);

npow = xrep.^drep;
nmonom = prod(npow,2);          % nmonom is nterms-by-1-by-npts

% Evalute polynomial--mulitply by coefs and sum monomials
nmonomrep = repmat(nmonom,[1 lp 1]); % nmonomrep is nterms-by-lp-by-npts
coefrep = repmat(coef,[1 1 npts]);   % coefrep is nterms-by-lp-by-npts
val = sum(coefrep.*nmonomrep,1);     % val is 1-by-lp-by-npts

% Reshape val from 1-by-lp-by-npts to lp-by-npts
if npts==1
    val = val(:);
elseif lp==1
    val = val(:)';
else
    val = squeeze(val);
end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function set1001d(m)

% sets parameters of DT3157 according to
% http://cmp.felk.cvut.cz/cmp/hardware/pulnix.html#TM-1001digital
% for single camera directly connected to frame grabber

% m=openfg('DT3157');
   setfg(m,'FrameType','NonInterlaced')
   setfg(m,'CameraType',4);				% 8 bit mode
   setfg(m,'ClockSource','External')	%(switch camera to NSP)
   setfg(m,'SyncSentinel','Off')

   setfg(m,'ExtOnLoToHi',0)   % pixel clock sampling edge
   setfg(m,'LineOnLoToHi',0)  % line data valid after edge
   setfg(m,'FieldOnLoToHi',0) % frame data valid after edge

   setfg(m,'TotalPixPerLine',1280)
   setfg(m,'FirstActivePixel',142)
   setfg(m,'ActivePixelCount',1008)

   setfg(m,'TotalLinesperFld',1100)
   setfg(m,'FirstActiveLine',24)
   setfg(m,'ActiveLineCount',1017)

   setfg(m,'FrameLeft',0)
   setfg(m,'HorFrameInc',1)
   setfg(m,'FrameWidth',1008)
   setfg(m,'FrameTop',0)
   setfg(m,'VerFrameInc',1)
   setfg(m,'FrameHeight',1017)

   setfg(m,'digioconfig',255);
   setfg(m,'digitalio',255)

author: Sable

repo: mcbench-benchmarks


text:
function sgmga_unique_index_tests ( )

%****************************************************************************80
%
%% SGMGA_UNIQUE_INDEX_TESTS calls SGMGA_UNIQUE_INDEX_TEST with various arguments.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 June 2010
%
%  Author:
%
%    John Burkardt
%
%  Local Parameters:
%
%    Local, real TOL, a tolerance for point equality.
%    A value of sqrt ( eps ) is reasonable, and will allow the code to
%    consolidate points which are equal, or very nearly so.  A value of
%    -1.0, on the other hand, will force the code to use every point, regardless
%    of duplication.
%
  addpath ( '../sandia_rules' );

  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SGMGA_UNIQUE_INDEX_TESTS\n' );
  fprintf ( 1, '  Call SGMGA_UNIQUE_INDEX_TEST with various arguments\n' );
%
%  Set the point equality tolerance.
%
  tol = sqrt ( eps );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  All tests will use a point equality tolerance of %f\n', tol );

  dim_num = 2;
  for dim = 1 : dim_num
    importance(dim,1) = 1.0;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 1 ]';
  growth = [ 6, 6 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  sgmga_unique_index_test ( dim_num, importance, level_weight, ...
    level_max_min, level_max_max, rule, growth, np, p, tol );

  dim_num = 2;
  for dim = 1 : dim_num
    importance(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 1 ]';
  growth = [ 6, 6 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  sgmga_unique_index_test ( dim_num, importance, level_weight, ...
    level_max_min, level_max_max, rule, growth, np, p, tol );

  dim_num = 3;
  for dim = 1 : dim_num
    importance(dim,1) = 1.0;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 1, 1 ]';
  growth = [ 6, 6, 6 ]';
  np = [ 0, 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  sgmga_unique_index_test ( dim_num, importance, level_weight, ...
    level_max_min, level_max_max, rule, growth, np, p, tol );

  dim_num = 3;
  for dim = 1 : dim_num
    importance(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 1, 1 ]';
  growth = [ 6, 6, 6 ]';
  np = [ 0, 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  sgmga_unique_index_test ( dim_num, importance, level_weight, ...
    level_max_min, level_max_max, rule, growth, np, p, tol );

  dim_num = 2;
  for dim = 1 : dim_num
    importance(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 3 ]';
  growth = [ 6, 6 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  sgmga_unique_index_test ( dim_num, importance, level_weight, ...
    level_max_min, level_max_max, rule, growth, np, p, tol );

  dim_num = 2;
  for dim = 1 : dim_num
    importance(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 4 ]';
  growth = [ 6, 3 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  sgmga_unique_index_test ( dim_num, importance, level_weight, ...
    level_max_min, level_max_max, rule, growth, np, p, tol );

  dim_num = 2;
  for dim = 1 : dim_num
    importance(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 7 ]';
  growth = [ 6, 3 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  sgmga_unique_index_test ( dim_num, importance, level_weight, ...
    level_max_min, level_max_max, rule, growth, np, p, tol );

  dim_num = 2;
  for dim = 1 : dim_num
    importance(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 8 ]';
  growth = [ 6, 3 ]';
  np = [ 0, 1 ]';
  np_sum = sum ( np(1:dim_num) );
  p(1) = 1.5;
  sgmga_unique_index_test ( dim_num, importance, level_weight, ...
    level_max_min, level_max_max, rule, growth, np, p, tol );

  dim_num = 2;
  for dim = 1 : dim_num
    importance(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 2, 9 ]';
  growth = [ 6, 3 ]';
  np = [ 0, 2 ]';
  np_sum = sum ( np(1:dim_num) );
  p(1,1) = 0.5;
  p(2,1) = 1.5;
  sgmga_unique_index_test ( dim_num, importance, level_weight, ...
    level_max_min, level_max_max, rule, growth, np, p, tol );

  dim_num = 2;
  for dim = 1 : dim_num
    importance(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 6, 10 ]';
  growth = [ 3, 4 ]';
  np = [ 1, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p(1) = 2.0;
  sgmga_unique_index_test ( dim_num, importance, level_weight, ...
    level_max_min, level_max_max, rule, growth, np, p, tol );

  dim_num = 3;
  for dim = 1 : dim_num
    importance(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 4, 5 ]';
  growth = [ 6, 3, 3 ]';
  np = [ 0, 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  sgmga_unique_index_test ( dim_num, importance, level_weight, ...
    level_max_min, level_max_max, rule, growth, np, p, tol );
%
%  Try a case with a dimension of "0 importance".
%
  dim_num = 3;
  importance = [ 1.0, 0.0, 1.0 ];
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max_min = 0;
  level_max_max = 3;
  rule = [ 1, 1, 1 ]';
  growth = [ 6, 6, 6 ]';
  np = [ 0, 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  sgmga_unique_index_test ( dim_num, importance, level_weight, ...
    level_max_min, level_max_max, rule, growth, np, p, tol );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SGMGA_UNIQUE_INDEX_TESTS:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../sandia_rules' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function p = addEquality(p,row)
% Equalities on top so easy
p.F_struc = [row;p.F_struc];
p.K.f = p.K.f + size(row,1);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function varargout = isreal(varargin)
%ISREAL   Real-valued CHEBFUN2 test.
%   ISREAL(F) returns logical true if F does not have an imaginary part and
%   false otherwise.
%
%   ~ISREAL(F) detects CHEBFUN2 object that have an imaginary part even if it is
%   all zero.

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

[varargout{1:nargout}] = isreal@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%COPYMAKEBORDER  Forms a border around an image
%
%     dst = cv.copyMakeBorder(src, top, bottom, left, right)
%     dst = cv.copyMakeBorder(src, [top, bottom, left, right])
%     [...] = cv.copyMakeBorder(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __src__ Source image.
% * __top__, __bottom__, __left__, __right__ Parameter specifying how many
%   pixels in each direction from the source image rectangle to extrapolate.
%   For example, `top=1`, `bottom=1`, `left=1`, `right=1` mean that 1
%   pixel-wide border needs to be built.
%
% ## Output
% * __dst__ Destination image of the same type as `src` and the size
%   `[size(src,1)+top+bottom, size(src,2)+left+right, size(src,3)]`.
%
% ## Options
% * __BorderType__ Border type, one of:
%   * __Constant__ `iiiiii|abcdefgh|iiiiiii` with some specified `i`
%   * __Replicate__ `aaaaaa|abcdefgh|hhhhhhh`
%   * __Reflect__ `fedcba|abcdefgh|hgfedcb`
%   * __Reflect101__ `gfedcb|abcdefgh|gfedcba`
%   * __Wrap__ `cdefgh|abcdefgh|abcdefg`
%   * __Default__ same as 'Reflect101' (default)
% * __Value__ Border value when `BorderType` is 'Constant'. default zeros
%
% The function copies the source image into the middle of the destination
% image. The areas to the left, to the right, above and below the copied
% source image will be filled with extrapolated pixels. This is not what
% filtering functions based on it do (they extrapolate pixels on-fly), but
% what other more complex functions, including your own, may do to simplify
% image boundary handling.
%
% See also: cv.borderInterpolate, padarray
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [node,elem] = squaremesh(square,h)
%% SQUAREMESH uniform mesh of a square
%
% [node,elem] = squaremesh([x0,x1,y0,y1],h) generates a uniform mesh of the
% rectangle [x0,x1]*[y0,y1] with mesh size h.
%
% Example
%
%   [node,elem] = squaremesh([0,1,0,1],0.2);
%   showmesh(node,elem);
%   findnode(node);
%   findelem(node,elem);
%
% See also: squarequadmesh, cubehexmesh
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Generate nodes
x0 = square(1); x1 = square(2);
y0 = square(3); y1 = square(4);
[x,y] = meshgrid(x0:h:x1,y0:h:y1);
node = [x(:),y(:)];

%% Generate elements
ni = size(x,1); % number of rows
N = size(node,1);
t2nidxMap = 1:N-ni;
topNode = ni:ni:N-ni;
t2nidxMap(topNode) = [];
k = (t2nidxMap)';
elem = [k+ni k+ni+1 k; k+1 k k+ni+1];
% 4 k+1 --- k+ni+1 3
%    |        |
% 1  k  ---  k+ni  2

repo: ifem


text:
function [y, yvar, gene, times, scale, rawExp, rawVar] = gpsimLoadBarencoTestData

% GPSIMLOADBARENCOTESTDATA Load in Martino Barenco's test data as processed by mmgMOS.
% FORMAT
% DESC loads in from the two Excel spread sheets
% (resultsMartino_exprs.xls and resultsMartino_se.xls) the data
% from the Barenco et al paper as processed by mmgMOS.
% RETURN y : the normalised expression levels.
% RETURN yvar : the variance of the normalised expression levels.
% RETURN gene : the gene names and Affymetrix array tags.
% RETURN times : the times of the expression measurements.
% RETURN scale : the scaling factor applied to normalise.
% RETURN rawExp : the raw gene expresion level.
% RETURN rawVar : the raw variance of the gene expression.
%
% SEEALSO : demBarencoRank
%
% COPYRIGHT : Neil D. Lawrence, 2007

% SHEFFIELDML

if exist('./data/barencoDataTest.mat') == 2
  load('./data/barencoDataTest.mat');
else

  % These excel files include results processed directly from the
  % cel files using the mmgMOS algorithm (Xuejun's code).

  % These are the expression levels.
  [numeric1, txt1] = xlsread('./data/resultsMartino_exprs.xls');
  headTxt1 = txt1(1, 2:end);
  tagTxt1 = txt1(2:end, 1);

  % These are the standard deviations.
  [numeric2, txt2] = xlsread('./data/resultsMartino_se.xls');
  headTxt2 = txt2(1, 2:end);
  tagTxt2 = txt2(2:end, 1);

  if(any(~strcmp(tagTxt2(:), tagTxt1(:))))
    error('Two files are not in same order');
  end
  if(any(~strcmp(headTxt2(:), headTxt1(:))))
    error('Two files are not in same order');
  end

  clear gene, clear ind
  % Gene IDs
  % DDB2
  gene{1, 1} = '203409_at';
  gene{1, 2} = 'DDB2';
  % BIK
  gene{2, 1} = '205780_at';
  gene{2, 2} = 'BIK';
  % TNFRSF10b (other tags include 209294_x_at and 210405_x_at)
  gene{3, 1} = '209295_at';
  gene{3, 2} = 'TNFRSF10b';
  % p21 --- we think this is CIp1/p21
  gene{4, 1} = '202284_s_at';
  gene{4, 2} = 'CIp1/p21';
  % p26 --- named as sesn1 in the platform.
  gene{5, 1} = '218346_s_at';
  gene{5, 2} = 'p26 sesn1';

  for i = 1:length(gene)
    match = find([strcmp(gene{i, 1}, tagTxt1(:))]);
    if length(match)~=1
      error('Too many or too few matches.');
    else
      ind(i) = match(1);
    end
  end
  order = [1 4 5 6 7 2 3 8 11 12 13 14 9 10 15 18 19 20 21 16 17];

  % Perform some normalisation.
  % Make sure that the average for each slide in log space is the
  % same.
  mVal = zeros(size(mean(numeric1)));
  mVal = mVal - mean(mVal);

  rawExp = numeric1(ind, order)';

  for i = 1:size(rawExp, 2)
    rawExp(:, i) = rawExp(:, i) - mVal';
  end

  rawVar = numeric2(ind, order)';
  rawVar = rawVar.*rawVar; % convert standard deviations to variances.

  yFull = exp(rawExp + rawVar/2);  % Logs are normally distributed
                                   % ... recover mean in exp space.
  yFullVar = (exp(rawVar)-1).*exp(2*rawExp + rawVar); % Logs are
                                                      % normally
                                                      % distributed
                                                      % ... recover
                                                      % variance in exp
                                                      % space.

  % Rescale so that average standard deviation of curves is 1.
  scale = mean(sqrt(var(yFull)));
  yFull = yFull/scale;
  yFullVar = yFullVar/(scale*scale);
  y{1} = yFull(1:7, :);
  y{2} = yFull(8:14, :);
  y{3} = yFull(15:21, :);
  yvar{1} = yFullVar(1:7, :);
  yvar{2} = yFullVar(8:14, :);
  yvar{3} = yFullVar(15:21, :);
  times = [0 2 4 6 8 10 12]';
  save('./data/barencoDataTest.mat', 'y', 'yvar', 'gene', 'times', 'scale', 'rawVar', 'rawExp');
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function psodemo(DemoMode)
% Runs the PSO on a few demonstration functions, which should be located
% in the <./testfcns> directory.
%
% For less intensive 3D graphics, run with input argument 'fast', viz.:
% >> psodemo('fast')
%
% S. Chen, Dec 2009.
% Available as part of "Another Particle Swarm Toolbox" at:
% http://www.mathworks.com/matlabcentral/fileexchange/25986
% Distributed under BSD license.

workingdir = pwd ;
testdir = ls('testf*') ;
if ~isempty(testdir), cd(testdir), end

[testfcn,testdir] = uigetfile('*.m','Load demo function for PSO') ;
if ~testfcn
    cd(workingdir)
    return
elseif isempty(regexp(testfcn,'\.m(?!.)','once'))
    error('Test function must be m-file')
else
    cd(testdir)
end

fitnessfcn = str2func(testfcn(1:regexp(testfcn,'\.m(?!.)')-1)) ;
cd(workingdir)

options = fitnessfcn('init') ;

if any(isfield(options,{'options','Aineq','Aeq','LB','nonlcon'}))
    % Then the test function gave us a (partial) problem structure.
    problem = options ;
else
    % Aineq = [1 1] ; bineq = [1.2] ; % Test case for linear constraint
    problem.options = options ;
    problem.Aineq = [] ; problem.bineq = [] ;
    problem.Aeq = [] ; problem.beq = [] ;
    problem.LB = [] ; problem.UB = [] ;
    problem.nonlcon = [] ;
end

problem.fitnessfcn = fitnessfcn ;
problem.nvars = 2 ;

if ~nargin
    problem.options.DemoMode = 'pretty' ;
else
    problem.options.DemoMode = DemoMode ;
end
problem.options.PlotFcns = {@psoplotbestf,@psoplotswarmsurf} ;
% problem.options.VelocityLimit = 0.2 ;
if isfield(problem.options,'PopulationType') && ...
        ~strcmp(problem.options.PopulationType,'bitstring')
    problem.options.HybridFcn = @fmincon ;
end
% problem.options.Display = 'off' ;

if isfield(problem.options,'UseParallel') && ...
        strcmp(problem.options.UseParallel,'always')
    poolopen = false ;
    if ~matlabpool('size')
        matlabpool('open','AttachedFiles',{[pwd '\testfcns']}) ;
    else
        poolopen = true ;
        pctRunOnAll addpath([pwd '\testfcns']) ;
    end
end

pso(problem)

if isfield(problem.options,'UseParallel') && ...
        strcmp(problem.options.UseParallel,'always') && ~poolopen
    matlabpool('close');
end


text:
function [ a, l, r ] = r8vec_part_quick_a ( n, a )

%*****************************************************************************80
%
%% R8VEC_PART_QUICK_A reorders an R8VEC as part of a quick sort.
%
%  Discussion:
%
%    The routine reorders the entries of A.  Using A(1) as a key,
%    all entries of A that are less than or equal to the key will
%    precede the key which precedes all entries that are greater than the key.
%
%  Example:
%
%    Input:
%
%      N = 8
%
%      A = ( 6, 7, 3, 1, 6, 8, 2, 9 )
%
%    Output:
%
%      L = 3, R = 6
%
%      A = ( 3, 1, 2, 6, 6, 8, 9, 7 )
%            -------        -------
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 October 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries of A.
%
%    Input, realA(N), the array to be checked.
%
%    Output, real A(N), has been reordered as described above.
%
%    Output, integer L, R, the indices of A that define the three segments.
%    Let KEY = the input value of A(1).  Then
%    I <= L                 A(I) < KEY;
%         L < I < R         A(I) = KEY;
%                 R <= I    KEY < A(I).
%
  if ( n < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_PART_QUICK_A - Fatal error!\n' );
    fprintf ( 1, '  N < 1.\n' );
    error ( 'R8VEC_PART_QUICK_A - Fatal error!' );
  elseif ( n == 1 )
    l = 0;
    r = 2;
    return;
  end

  key = a(1);
  m = 1;
%
%  The elements of unknown size have indices between L+1 and R-1.
%
  l = 1;
  r = n + 1;

  for i = 2 : n

    if ( key < a(l+1) )
      r = r - 1;
      [ a(r), a(l+1) ] = r8_swap ( a(r), a(l+1) );
    elseif ( a(l+1) == key )
      m = m + 1;
      [ a(m), a(l+1) ] = r8_swap ( a(m), a(l+1) );
      l = l + 1;
    elseif ( a(l+1) < key )
      l = l + 1;
    end

  end
%
%  Now shift small elements to the left, and KEY elements to center.
%
  for i = 1 : l - m
    a(i) = a(i+m);
  end
%
%  Out of bounds here, occasionally.
%
  l = l - m;

  a(l+1:l+m) = key;

  return
end


text:
function [cm, cSq] = DiscreteFrechetDist(P,Q,dfcn)
% Calculates the discrete Frechet distance between curves P and Q
%
% [cm, cSq] = DiscreteFrechetDist(P,Q)
% [cm, cSq] = DiscreteFrechetDist(...,dfcn)
%
% P and Q are two sets of points that define polygonal curves with rows of
% vertices (data points) and columns of dimensionality. The points along
% the curves are taken to be in the order as they appear in P and Q.
%
% Returned in cm is the discrete Frechet distance, aka the coupling
% measure, which is zero when P equals Q and grows positively as the curves
% become more dissimilar.
%
% The optional dfcn argument allows the user to specify a function with
% which to calculate distance between points in P and Q. If not provided,
% the L2 norm is used.
%
% The secondary output, cSq, is the coupling sequence, that is, the
% sequence of steps along each curve that must be followed to achieve the
% minimum coupling distance, cm. The output is returned in the form of a
% matrix with column 1 being the index of each point in P and column 2
% being the index of each point in Q. (NOTE: the coupling sequence is not
% unique in general)
%
% Explanation:
% The Frechet distance is a measure of similarity between to curves, P and
% Q. It is defined as the minimum cord-length sufficient to join a point
% traveling forward along P and one traveling forward along Q, although the
% rate of travel for either point may not necessarily be uniform.
%
% The Frechet distance, FD, is not in general computable for any given
% continuous P and Q. However, the discrete Frechet Distance, also called
% the coupling measure, cm, is a metric that acts on the endpoints of
% curves represented as polygonal chains. The magnitude of the coupling
% measure is bounded by FD plus the length of the longest segment in either
% P or Q,  and approaches FD in the limit of sampling P and Q.
%
% This function implements the algorithm to calculate discrete Frechet
% distance outlined in:
% T. Eiter and H. Mannila. Computing discrete Frechet distance. Technical
% Report 94/64, Christian Doppler Laboratory, Vienna University of
% Technology, 1994.
%
%
%
% EXAMPLE:
% % create data
% t = 0:pi/8:2*pi;
% y = linspace(1,3,6);
% P = [(2:7)' y']+0.3.*randn(6,2);
% Q = [t' sin(t')]+2+0.3.*randn(length(t),2);
% [cm, cSq] = DiscreteFrechetDist(P,Q);
% % plot result
% figure
% plot(Q(:,1),Q(:,2),'o-r','linewidth',3,'markerfacecolor','r')
% hold on
% plot(P(:,1),P(:,2),'o-b','linewidth',3,'markerfacecolor','b')
% title(['Discrete Frechet Distance of curves P and Q: ' num2str(cm)])
% legend('Q','P','location','best')
% line([2 cm+2],[0.5 0.5],'color','m','linewidth',2)
% text(2,0.4,'dFD length')
% for i=1:length(cSq)
%   line([P(cSq(i,1),1) Q(cSq(i,2),1)],...
%        [P(cSq(i,1),2) Q(cSq(i,2),2)],...
%        'color',[0 0 0]+(i/length(cSq)/1.35));
% end
% axis equal
% % display the coupling sequence along with each distance between points
% disp([cSq sqrt(sum((P(cSq(:,1),:) - Q(cSq(:,2),:)).^2,2))])
%
%
%
% %%% ZCD June 2011 %%%
% %%% edits ZCD May 2013: 1) remove excess arguments to internal functions
% and persistence for speed, 2) added example, 3) allowed for user defined
% distance function, 4) added aditional output option for coupling sequence
%

% size of the data curves
sP = size(P);
sQ = size(Q);

% check validity of inputs
if sP(2)~=sQ(2)
    error('Curves P and Q must be of the same dimension')
elseif sP(1)==0
    cm = 0;
    return;
end

% initialize CA to a matrix of -1s
CA = ones(sP(1),sQ(1)).*-1;

% distance function
if nargin==2
    dfcn = @(u,v) sqrt(sum( (u-v).^2 ));
end

% final coupling measure value
cm = c(sP(1),sQ(1));

% obtain coupling measure via backtracking procedure
if nargout==2
    cSq = zeros(sQ(1)+sP(1)+1,2);    % coupling sequence
    CApad = [ones(1,sQ(1)+1)*inf; [ones(sP(1),1)*inf CA]];  % pad CA
    Pi=sP(1)+1; Qi=sQ(1)+1; count=1;  % counting variables
    while Pi~=2 || Qi~=2
        % step down CA gradient
        [v,ix] = min([CApad(Pi-1,Qi) CApad(Pi-1,Qi-1) CApad(Pi,Qi-1)]);
        if ix==1
            cSq(count,:) = [Pi-1 Qi];
            Pi=Pi-1;
        elseif ix==2
            cSq(count,:) = [Pi-1 Qi-1];
            Pi=Pi-1; Qi=Qi-1;
        elseif ix==3
            cSq(count,:) = [Pi Qi-1];
            Qi=Qi-1;
        end
        count=count+1;
    end
    % format output: remove extra zeroes, reverse order, subtract off
    % padding value, and add in the last point
    cSq = [flipud(cSq(1:find(cSq(:,1)==0,1,'first')-1,:))-1; sP(1) sQ(1)];
end

% debug
% assignin('base','CAw',CA)

function CAij = c(i,j)
    % coupling search function
    if CA(i,j)>-1
        % don't update CA in this case
        CAij = CA(i,j);
    elseif i==1 && j==1
        CA(i,j) = dfcn(P(1,:),Q(1,:));     % update the CA permanent
        CAij = CA(i,j);                    % set the current relevant value
    elseif i>1 && j==1
        CA(i,j) = max( c(i-1,1), dfcn(P(i,:),Q(1,:)) );
        CAij = CA(i,j);
    elseif i==1 && j>1
        CA(i,j) = max( c(1,j-1), dfcn(P(1,:),Q(j,:)) );
        CAij = CA(i,j);
    elseif i>1 && j>1
        CA(i,j) = max( min([c(i-1,j), c(i-1,j-1), c(i,j-1)]),...
            dfcn(P(i,:),Q(j,:)) );
        CAij = CA(i,j);
    else
        CA(i,j) = inf;
    end
end     % end function, c

end     % end main function


% Gets the number of seconds since Jan 1 1970 (epoch)
%
% @return seconds Number of seconds since epoch
function [seconds] = get_seconds_since_epoch()
    if (is_octave())
        seconds = uint64(time());
    else
        seconds = uint64(convertTo(datetime('now', 'TimeZone', 'UTC'), 'epochtime'));
    end
end


text:
function [xV,yV,zV] = mtoxyz_Dose(rV,cV,sV,planC,doseNum)
%"mtoxyz_Dose"
%   Convert from rcs coordinates to xyz coordinates, in the nonuniformized
%   dataset.
%
%   if flag = 'uniform', uses the uniformized data.
%
%   JRA 07/15/04
%
%Usage:
%   [xV,yV,zV] = mtoxyz_Dose(rV,cV,sV,planC,doseNum)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};
 [xVals, yVals, zVals] = getDoseXYZVals(planC{indexS.dose}(doseNum));

 xV = interp1(1:length(xVals), xVals, cV, 'linear', 'extrap');
yV = interp1(1:length(yVals), yVals, rV, 'linear', 'extrap');
zV = interp1(1:length(zVals), zVals, sV, 'linear', 'extrap');
return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function sparse_grid_monomial_test ( rule, dim_num, level_max, degree_max )

%*****************************************************************************80
%
%% SPARSE_GRID_MONOMIAL_TEST tests monomial exactness of the sparse grid rules.
%
%  Discussion:
%
%    This test is going to check EVERY monomial of total degree DEGREE_MAX
%    or less.  Even for a moderately high dimension of DIM_NUM = 10, you
%    do NOT want to use a large value of DEGREE_MAX, since there are
%
%      1         monomials of total degree 0,
%      DIM_NUM   monomials of total degree 1,
%      DIM_NUM^2 monomials of total degree 2,
%      DIM_NUM^3 monomials of total degree 3, and so on.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 July 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer RULE, the index of the rule.
%    1, "CC", Clenshaw Curtis Closed Fully Nested rule.
%    2, "F1", Fejer 1 Open Fully Nested rule.
%    3, "F2", Fejer 2 Open Fully Nested rule.
%    4, "GP", Gauss Patterson Open Fully Nested rule.
%    5, "GL", Gauss Legendre Open Weakly Nested rule.
%    6, "GH", Gauss Hermite Open Weakly Nested rule.
%    7, "LG", Gauss Laguerre Open Non Nested rule.
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the level.
%
%    Input, integer DEGREE_MAX, the maximum monomial total
%    degree to check.
%
  level_min = max ( 0, level_max + 1 - dim_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_MONOMIAL_TEST\n' );
  fprintf ( 1, '  Check the exactness of a sparse grid quadrature rule,\n' );
  fprintf ( 1, '  applied to all monomials of orders 0 to DEGREE_MAX.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  For cases where the dimension is greater than 1,\n' );
  fprintf ( 1, '  many sparse grid of this level have accuracy through\n' );
  fprintf ( 1, '  monomials of total degree   %d\n', 2 * level_max + 1 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension DIM_NUM = %d\n', dim_num );
  fprintf ( 1, '  LEVEL_MIN =                 %d\n', level_min );
  fprintf ( 1, '  LEVEL_MAX =                 %d\n', level_max );
  fprintf ( 1, '  1D quadrature index RULE =  %d\n', rule );
  fprintf ( 1, '  Check up to DEGREE_MAX =    %d\n', degree_max );
%
%  Determine the number of points in the rule.
%
  point_num = levels_index_size ( dim_num, level_max, rule );

  fprintf ( 1, '  Unique points in the grid = %d\n', point_num );
%
%  Compute the weights and points.
%
  [ grid_weight, grid_point ] = sparse_grid ( dim_num, level_max, rule, ...
    point_num );
%
%  Compare exact and estimated values of the integrals of various monomials.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '      Error      Total   Monomial\n' );
  fprintf ( 1, '                 Degree  Exponents\n' );
  fprintf ( 1, '\n' );

  for degree = 0 : degree_max

    expon = [];
    more = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ expon, more, h, t ] = comp_next ( degree, dim_num, expon, more, h, t );

      quad_error = monomial_quadrature ( dim_num, expon, point_num, ...
        grid_weight, grid_point, rule );

      fprintf ( 1, '  %14e  %2d  ', quad_error, degree );
      for dim = 1 : dim_num
        fprintf ( 1, '%2d', expon(dim) );
      end
      fprintf ( 1, '\n' );

      if ( ~more )
        break
      end

    end

    if ( 1 < dim_num )
      fprintf ( 1, '\n' );
    end

  end

  return
end


text:
function h = qq_plot(varargin)
%  Quantile-quantile plot with patch option
%
%   QQ_PLOT(Y) displays a quantile-quantile plot of the sample quantiles of
%   Y versus theoretical quantiles from a normal distribution. If the
%   distribution of Y is normal, the plot will be close to linear.
%
%   QQ_PLOT(X,Y) displays a quantile-quantile plot of two samples. If the
%   samples come from the same distribution, the plot will be linear.
%
%   The inputs X and Y should be numeric and have an equal number of
%   elements; every element is treated as a member of the sample.
%
%   The plot displays the sample data with the plot symbol 'x'.
%   Superimposed on the plot is a dashed straight line connecting the first
%   and third quartiles.
%
%   QQ_PLOT(...,MODE) allows the appearance of the plot to be configured.
%   With MODE='line' (default), the plot appears as described above. With
%   MODE='patch', the data are plotted as a patch object, with the area
%   bound by the x-distribution and the linear fit shaded grey. With
%   mode='both' the two appearances are combined.
%
%   QQ_PLOT(...,MODE,METHOD) and QQ_PLOT(...,[],METHOD) allows the method
%   for calculating the quartiles, used for the fit line, to be specified.
%   The default is 'R-8'. Type 'help quantile2' for more information. The
%   latter form of the function call uses the default mode.
%
%   H = QQ_PLOT(...) returns a two- or three-element vector of handles to
%   the plotted object. The nature of the handles depends upon the mode. In
%   all cases, the first handle is to the sample data, the second handle is
%   to the fit line. With MODE='patch' or MODE='both', there is third
%   handle to the patch object.
%
%   Example
%
%      % Display Q-Q plots for the rand and randn functions
%      figure
%      subplot(2,1,1)
%      qq_plot(rand(20),'patch')
%      subplot(2,1,2)
%      h = qq_plot(randn(20),'patch');
%      set(h(3),'FaceColor','r') % change fill color

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LICENSE FILE:
% -----------------------------------------------------------------------
% Copyright (c) 2015, Christopher Hummersone
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ==========================================================
% Last changed:     $Date: 2015-05-20 16:28:56 +0100 (Wed, 20 May 2015) $
% Last committed:   $Revision: 380 $
% Last changed by:  $Author: ch0022 $
% ==========================================================

    %% determine X and Y

    IXn = cellfun(@(x) isnumeric(x) & ~isempty(x),varargin);

    switch sum(IXn)
        case 0
            error('No input data specified')
        case 1
            % compare to normal distrbution
            Y = get_input_sample(varargin,IXn);
            p = (.5:length(Y))/length(Y);
            X = sqrt(2)*erfinv(2*p - 1);
            x_label = 'Standard normal quantiles';
            y_label = 'Sample quantiles';
        case 2
            % compare to input data distribution
            Y = get_input_sample(varargin,find(IXn,1,'last'));
            X = get_input_sample(varargin,find(IXn,1,'first'));
            assert(isequal(size(X),size(Y)),'Input data must be the same size')
            x_label = 'X quantiles';
            y_label = 'Y quantiles';
        otherwise
            error('Unknown input specified')
    end

    %% determine mode and method

    % find inputs
    IXc = cellfun(@(x) ischar(x) | isempty(x),varargin);
    switch sum(IXc)
        case 0
            mode = [];
            method = [];
        case 1
            mode = varargin{IXc};
            method = [];
        case 2
            mode = varargin{find(IXc,1,'first')};
            method = varargin{find(IXc,1,'last')};
        otherwise
            error('Unknown string specified')
    end

    % defaults
    if isempty(mode)
        mode = 'line';
    end
    if isempty(method)
        method = 'R-8';
    end

    %% calculate fit to first and third quartile

    % quartiles
    q1x = quantile2(X,.25,[],method);
    q3x = quantile2(X,.75,[],method);
    q1y = quantile2(Y,.25,[],method);
    q3y = quantile2(Y,.75,[],method);

    % slope
    slope = (q3y-q1y)./(q3x-q1x);
    centerx = (q1x+q3x)/2;
    centery = (q1y+q3y)/2;

    % fit
    maxx = max(X);
    minx = min(X);
    maxy = centery + slope.*(maxx - centerx);
    miny = centery - slope.*(centerx - minx);

    % lines
    X_fit = linspace(minx,maxx,length(X));
    Y_fit = linspace(miny,maxy,length(X));

    %% plot data

    hold on

    if strcmpi(mode,'patch') || strcmpi(mode,'both')
        Hp = patch([X fliplr(X_fit)],[Y fliplr(Y_fit)],[0.5 0.5 0.5]);
        set(Hp,'edgecolor','none')
    else
        Hp = NaN;
    end

    switch lower(mode)
        case 'line'
            linestyle = 'xk';
        case 'patch'
            linestyle = '-k';
        case 'both'
            linestyle = 'xk';
        otherwise
            error('Unknown mode specified')
    end

    H = plot(X,Y,linestyle,X_fit,Y_fit,'--k');

    hold off

    % axis labels
    xlabel(x_label)
    ylabel(y_label)

    box on; axis tight
    set(gca,'layer','top')

    % return handle
    if nargout>0
        h = H;
        if isobject(Hp) || ishandle(Hp)
            h = [h; Hp];
        end
    end

end

function Z = get_input_sample(z,IX)
%GET_INPUT_SAMPLE get sample, order, and convert to vector
    Z = z{IX};
    Z = sort(Z(:))';
end

author: brainstorm-tools


text:
function xyz = plane_normal_qr_to_xyz ( pp, normal, pq, pr, n, qr )

%*****************************************************************************80
%
%% PLANE_NORMAL_QR_TO_XYZ: QR_TO_XYZ coordinates for a normal form plane.
%
%  Discussion:
%
%    The normal form of a plane in 3D is:
%
%      PP is a point on the plane,
%      NORMAL is a normal vector to the plane.
%
%    Two vectors PQ and PR can be computed with the properties that
%    * NORMAL, PQ and PR are pairwise orthogonal;
%    * PQ and PR have unit length;
%    * every point P in the plane has a "QR" representation
%      as P = PP + q * PQ + r * PR.
%
%    This function is given the QR coordinates of a set of points on the
%    plane, and returns the XYZ coordinates.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real PP(3), a point on the plane.
%
%    Input, real NORMAL(3), a normal vector N to the plane.  The
%    vector must not have zero length, but it is not necessary for N
%    to have unit length.
%
%    Input, real PQ(3), a vector of unit length,
%    perpendicular to the vector N and the vector PR.
%
%    Input, real PR(3), a vector of unit length,
%    perpendicular to the vector N and the vector PQ.
%
%    Input, integer N, the number of points on the plane.
%
%    Input, real QR(2,N), the QR coordinates of the points.
%
%    Output, real XYZ(3,N), the XYZ coordinates of the points.
%
  xyz = repmat ( pp, 1, n ) + [ pq'; pr' ]' * qr(1:2,1:n);

  return
end


text:
function x = isemptydir(p,k)
  if(isdir(p))
    f = dir(p);
    x = ~(length(f) > k);
  else
    error('Error: %s is not a directory');
  end
end

author: andrewssobral

repo: lrslibrary


text:
function [rlvm, frvals, frvecs, trnsfrmd, mn, dv] = TSTL_pca(data, mode, maxpercent, sil)

%   [rlvm, frvals, frvecs, trnsfrmd, mn, dv] = pca(data, mode, maxpercent, silent)
%
%   principal component analysis of column orientated data set <data>
%
%   input arguments :
%
%   - each row of data is one 'observation', e.g. the sample values of
%     all channels in a multichannel measurement at one point in time
%
%   - mode can be one of the following : 'normalized' (default), 'mean', 'raw'
%     - in mode 'normalized' each column of data is centered by removing its mean
%       and then normalized by dividing through its standard deviation before
%       the covariance matrix is calculated
%     - in mode 'mean' only the mean of every column of data is removed
%     - in mode 'raw' no preprocessing is applied to data
%
%   - maxpercent gives the limit of the accumulated percentage of the resulting
%     eigenvalues, default is 95 %
%
%   - silent is an optional flag which supresses output of text and plot on the matlab
%     screen. Returned values (see below) are in no way affected
%
%   output arguments :
%
%   - rlvm : number of relevant modes to reach maxpercent
%   - frvals : first rlvm relevant eigenvalues
%   - frvecs : first rlvm relevant eigenvectors (eigenmodes, or principal components)
%   - trnsfrmd : data points transformed to the coordinate system given be the eigenvectors
%   - mn : mean of the original data set (only in mode 'mean' and 'normalized')
%   - dv : standard deviation of the original data set (only in mode 'normalized')
%
%
%   To compute an approximation of data : data = trnsfrmd * frvecs'
%
%   Christian Merkwirth (cmerk) Maerz 1997
%   cmerk Jan.  1998

global silent

if nargin  < 1, help(mfilename); end

if nargin < 2
	mode = 'normalized';
end
if nargin < 3
	maxpercent = 95;
end
if nargin < 4
	silent = 0;
else
	silent = 1;
end

%rang = rank(data);
[n,m] = size(data);

printline('principal component analysis')
printline(['on data set of size ' num2str(n) 'x' num2str(m)]);	% ' with rank ' num2str(rang)])
printline(['eigenvalues are computed up to ' num2str(maxpercent) ' percent']);

maxpercent = maxpercent/100;
mode = lower(mode); 		% no problems with uppercase letters

if strncmp(mode, 'r',1)
	mode = 'raw';
	printline('no data preprocessing');
elseif strncmp(mode, 'm',1)
	mode = 'mean';
	printline('removing mean from data set');
	mn = mean(data);
	data = data - repmat(mn, n, 1);
else
	mode = 'normalized';
	printline('removing mean and normalizing data');
	mn = mean(data);
	dv = std(data);
	data = data - repmat(mn, n, 1);
	data = data ./ repmat(dv, n, 1);
end

%sum(mean(data))/m	% test
%sum(std(data))/m	% test

if n>m
	printline('using direct method to compute covariance matrix');
	K =   data' * data;       % oder K = corrcoef(data)
	[Q,D] = eig(K);
 else
    printline('using indirect method to compute covariance matrix');
 	C = data * data';
 	[Q,D] = eig(C);
end

[evalues,index] = sort(diag(D));
evalues = flipud(evalues);
index = flipud(index);

total = sum(evalues);
rlvm = min(find(cumsum(evalues) >= (total*maxpercent)));

if isempty(rlvm)			% in case percentage was choosen ober 100 %,
	rlvm = length(evalues); % return all eigenvalues
end

frvals = evalues(1:rlvm);

if n>m
	frvecs = Q(:, index(1:rlvm));
	trnsfrmd=data*frvecs;
 else
 	scalefac = 1./ sqrt(evalues(1:rlvm));
 	for i = 1:rlvm
 		P(:,i) = Q(:,index(i)) * scalefac(i);
 	end
 	frvecs = data' * P; % '
	trnsfrmd=C*P;
end

if silent~=1
	bar(100*frvals/total);
	title('Eigenvalues in percent');
end

function printline(string)
global silent
if silent~=1
	disp(string)
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [Set_PSM, Set_TDS, OperatingTime] = ORCD (CT, OLF_Min, OLF_Max, LC, FC, C1, C2, TOP_Desired, IDMT_Saturation, Tolerance_Dn, Discrimination_Time, MinPSM, MaxPSM, StepPSM, MinTDS, MaxTDS, StepTDS)

[PSM_Max, PSM_Min] = PSM(CT, OLF_Max, OLF_Min, LC, MaxPSM, MinPSM, StepPSM);

TOP_Desired_For_Present_Relay=TOP_Desired+Discrimination_Time;
n=round((PSM_Max-PSM_Min)/StepPSM);
PSM_Temp=PSM_Min;
D=zeros(1,4);
    for i=1:n+1
        M=(FC/(PSM_Temp*CT));
        if (M>=IDMT_Saturation)
            M=IDMT_Saturation;
        end

        K=C1/((M^C2)-1);

        TDS_Desired=TOP_Desired_For_Present_Relay/K;

        TDS_Rounded_Up=ceil(TDS_Desired/StepTDS)*StepTDS;
        Difference_Up=abs(TDS_Rounded_Up-TDS_Desired);

        TDS_Rounded_Dn=floor(TDS_Desired/StepTDS)*StepTDS;
        Difference_Dn=abs(TDS_Rounded_Dn-TDS_Desired);

        Difference=min(Difference_Up,Difference_Dn);

        D(i,1)=PSM_Temp;
        D(i,2)=K;
        D(i,3)=TDS_Desired;
        D(i,4)=Difference;

        PSM_Temp=PSM_Temp+StepPSM;
        if (i~=n+1)
            D=[D;zeros(1,4)];
        end
    end
    D;
    [~,I]=min(D(:,4));
    Set_PSM=D(I,1);
    Set_TDS_Up=ceil(D(I,3)/StepTDS)*StepTDS;
    Set_TDS_Dn=floor(D(I,3)/StepTDS)*StepTDS;
    OperatingTime_Up=D(I,2)*Set_TDS_Up;
    OperatingTime_Dn=D(I,2)*Set_TDS_Dn;
    TOP_Error_Up=abs(TOP_Desired_For_Present_Relay-OperatingTime_Up);
    TOP_Error_Dn=abs(TOP_Desired_For_Present_Relay-OperatingTime_Dn);
    TOP_Error=min(TOP_Error_Up,TOP_Error_Dn);
    if (TOP_Error==TOP_Error_Dn && TOP_Error<=Tolerance_Dn)
        OperatingTime=OperatingTime_Dn;
        Set_TDS=Set_TDS_Dn;
    else
        OperatingTime=OperatingTime_Up;
        Set_TDS=Set_TDS_Up;
    end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Q = INVERSEKINEMATIC_tx200l(robot, T)
%   Solves the inverse kinematic problem for the STAUBLI tx200l robot
%   where:
%   robot stores the robot parameters.
%   T is an homogeneous transform that specifies the position/orientation
%   of the end effector.
%
%   A call to Q=INVERSEKINEMATIC_tx200l returns 8 possible solutions, thus,
%   Q is a 6x8 matrix where each column stores 6 feasible joint values.
%
%
%   Example code:
%
%   robot=load_robot('staubli', 'tx200l');
%   q = [0 0 0 0 0 0];
%   T = directkinematic(robot, q);
%   %Call the inversekinematic for this robot
%   qinv = inversekinematic(robot, T);
%   check that all of them are feasible solutions!
%   and every Ti equals T
%   for i=1:8,
%        Ti = directkinematic(robot, qinv(:,i))
%   end
%	See also DIRECTKINEMATIC.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Carlos Pardo Pla
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function q = inversekinematic_tx200l(robot, T)

%initialize q,
%eight possible solutions are generally feasible
q=zeros(6,8);

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%See geometry at the reference for this robot
L6=abs(d(6));

A1 = a(1);

%T= [ nx ox ax Px;
%     ny oy ay Py;
%     nz oz az Pz];
Px=T(1,4);
Py=T(2,4);
Pz=T(3,4);

%Compute the position of the wrist, being W the Z component of the end effector's system
W = T(1:3,3);

% Pm: wrist position
Pm = [Px Py Pz]' - L6*W;

%first joint, two possible solutions admited:
% if q(1) is a solution, then q(1) + pi is also a solution
q1=atan2(Pm(2), Pm(1));

%solve for q2
q2_1=solve_for_theta2(robot, [q1 0 0 0 0 0 0], Pm);

q2_2=solve_for_theta2(robot, [q1+pi 0 0 0 0 0 0], Pm);

%solve for q3
q3_1=solve_for_theta3(robot, [q1 0 0 0 0 0 0], Pm);

q3_2=solve_for_theta3(robot, [q1+pi 0 0 0 0 0 0], Pm);

%Arrange solutions, there are 8 possible solutions so far.
% if q1 is a solution, q1* = q1 + pi is also a solution.
% For each (q1, q1*) there are two possible solutions
% for q2 and q3 (namely, elbow up and elbow up solutions)
% So far, we have 4 possible solutions. Howefer, for each triplet (theta1, theta2, theta3),
% there exist two more possible solutions for the last three joints, generally
% called wrist up and wrist down solutions. For this reason,
%the next matrix doubles each column. For each two columns, two different
%configurations for theta4, theta5 and theta6 will be computed. These
%configurations are generally referred as wrist up and wrist down solution
q = [q1         q1         q1        q1       q1+pi   q1+pi   q1+pi   q1+pi;
     q2_1(1)    q2_1(1)    q2_1(2)   q2_1(2)  q2_2(1) q2_2(1) q2_2(2) q2_2(2);
     q3_1(1)    q3_1(1)    q3_1(2)   q3_1(2)  q3_2(1) q3_2(1) q3_2(2) q3_2(2);
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0];

%leave only the real part of the solutions
q=real(q);

%normalize q to [-pi, pi]
q(1,:) = normalize(q(1,:));
q(2,:) = normalize(q(2,:));

% solve for the last three joints
% for any of the possible combinations (theta1, theta2, theta3)
for i=1:2:size(q,2),
    qtemp = solve_spherical_wrist_tx200l(robot, q(:,i), T, 1); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i)=qtemp;

    qtemp = solve_spherical_wrist_tx200l(robot, q(:,i), T, -1); %wrist down
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i+1)=qtemp;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for second joint theta2, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q2 = solve_for_theta2(robot, q, Pm)

%Evaluate the parameters
d = eval(robot.DH.d);
a = eval(robot.DH.a);

%See geometry
L2=abs(a(2));
L3=abs(d(4));

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

beta = atan2(-p1(2), p1(1));
gamma =(acos((L2^2+r^2-L3^2)/(2*r*L2)));

if ~isreal(gamma)
    disp('WARNING:inversekinematic_tx200l: the point is not reachable for this configuration, imaginary solutions');
    %gamma = real(gamma);
end

%return two possible solutions
%elbow up and elbow down
%the order here is important and is coordinated with the function
%solve_for_theta3
q2(1) = pi/2 - beta - gamma; %elbow up
q2(2) = pi/2 - beta + gamma; %elbow down

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for third joint theta3, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q3 = solve_for_theta3(robot, q, Pm)

%Evaluate the parameters
d = eval(robot.DH.d);
a = eval(robot.DH.a);

%See geometry
L2=abs(a(2));
L3=abs(d(4));

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

eta =(acos((L2^2 + L3^2 - r^2)/(2*L2*L3)));

if ~isreal(eta)
   disp('WARNING:inversekinematic_tx200l: the point is not reachable for this configuration, imaginary solutions');
   %eta = real(eta);
end

%return two possible solutions
%elbow up and elbow down solutions
%the order here is important
q3(1) = pi - eta;
q3(2) = eta - pi;

% Solve the special case of this spherical wrist
% For wrists that whose reference systems have been placed as in the
% ABB IRB 140--> use solve_spherical_wrist2
% For wrists with the same orientation as in the KUKA KR30_jet
%--> use solve_spherical_wrist
function q = solve_spherical_wrist_tx200l(robot, q, T, wrist)

% T is the noa matrix defining the position/orientation of the end
% effector's reference system
vx6=T(1:3,1);
vz5=T(1:3,3); % The vector a z6=T(1:3,3) is coincident with z5

% Obtain the position and orientation of the system 3
% using the already computed joints q1, q2 and q3
T01=dh(robot, q, 1);
T12=dh(robot, q, 2);
T23=dh(robot, q, 3);
T03=T01*T12*T23;

vx3=T03(1:3,1);
vy3=T03(1:3,2);
vz3=T03(1:3,3);

% find z4 normal to the plane formed by z3 and a
vz4=cross(vz3, vz5);	% end effector's vector a: T(1:3,3)

% in case of degenerate solution,
% when vz3 and vz6 are parallel--> then z4=0 0 0, choose q(4)=0 as solution
if norm(vz4) <= 0.00000001
    if wrist == 1 %wrist up
        q(4)=0;
    else
        q(4)=-pi; %wrist down
    end
else
    %this is the normal and most frequent solution
    cosq4=wrist*dot(vy3,vz4);
    sinq4=wrist*dot(-vx3,vz4);
    q(4)=atan2(sinq4, cosq4);
end
%propagate the value of q(4) to compute the system 4
T34=dh(robot, q, 4);
T04=T03*T34;
vx4=T04(1:3,1);
vy4=T04(1:3,2);

% solve for q5
cosq5=dot(-vy4,vz5);
sinq5=dot(vx4,vz5);
q(5)=atan2(sinq5, cosq5);

%propagate now q(5) to compute T05
T45=dh(robot, q, 5);
T05=T04*T45;
vx5=T05(1:3,1);
vy5=T05(1:3,2);

% solve for q6
cosq6=dot(vx6,vx5);
sinq6=dot(vx6,vy5);
q(6)=atan2(sinq6, cosq6);

repo: ARTE


text:
% clear; close all; clc;

fs = 200;
t = 0:1/fs:5-1/fs;
x = sin(2*pi*1*t);

figure;
set(gcf,'position',[500 500 640 480]);
figsize = get(gcf,'position');
h = plot(t,x);

SliderH = uicontrol('style','slider','position',[70 6 512 20],'min',1,'max',5,'value',1);

while(1)
    plot(t, sin(2*pi*SliderH.Value*t));
    pause(0.01);
    cla
end

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols