text:
function out = MD_hrv_classic(y)
% MD_hrv_classic    Classic heart rate variability (HRV) statistics.
%
% Typically assumes an NN/RR time series in units of seconds.
%
%---INPUTS:
% y, the input time series.
%
% Includes:
%  (i) pNNx
%  cf. "The pNNx files: re-examining a widely used heart rate variability
%           measure", J.E. Mietus et al., Heart 88(4) 378 (2002)
%
%  (ii) Power spectral density ratios in different frequency ranges
%   cf. "Heart rate variability: Standards of measurement, physiological
%       interpretation, and clinical use",
%       M. Malik et al., Eur. Heart J. 17(3) 354 (1996)
%
%  (iii) Triangular histogram index, and
%
%  (iv) Poincare plot measures
%  cf. "Do existing measures of Poincare plot geometry reflect nonlinear
%       features of heart rate variability?"
%       M. Brennan, et al., IEEE T. Bio.-Med. Eng. 48(11) 1342 (2001)
%
% Code is heavily derived from that provided by Max A. Little:
% http://www.maxlittle.net/

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

% Standard defaults
diffy = diff(y);
N = length(y); % time-series length

% ------------------------------------------------------------------------------
% Calculate pNNx percentage
% ------------------------------------------------------------------------------
% pNNx: recommendation as per Mietus et. al. 2002, "The pNNx files: ...", Heart
% strange to do this for a z-scored time series...

Dy = abs(diffy);

% Anonymous function to do the PNNx calcualtion:
% proportion of difference magnitudes greater than X*sigma
PNNxfn = @(x) mean(Dy > x/1000);

out.pnn5  = PNNxfn(5); % 0.005*sigma
out.pnn10 = PNNxfn(10); % 0.01*sigma
out.pnn20 = PNNxfn(20); % 0.02*sigma
out.pnn30 = PNNxfn(30); % 0.03*sigma
out.pnn40 = PNNxfn(40); % 0.04*sigma

% ------------------------------------------------------------------------------
% Calculate PSD
% ------------------------------------------------------------------------------
% [Pxx, F] = psd(series,1024,1,hanning(1024),512);
[Pxx, F] = periodogram(y,hann(N)); % periodogram with hanning window

% ------------------------------------------------------------------------------
% Calculate spectral measures such as subband spectral power percentage, LF/HF ratio etc.
% ------------------------------------------------------------------------------
% LF/HF: as per Malik et. al. 1996, "Heart Rate Variability"
LF_lo = 0.04; % /pi -- fraction of total power (max F is pi)
LF_hi = 0.15;
HF_lo = 0.15;
HF_hi = 0.4;

fbinsize = F(2) - F(1);
indl  = ((F >= LF_lo) & (F <= LF_hi));
indh  = ((F >= HF_lo) & (F <= HF_hi));
indv  = (F <= LF_lo);
lfp   = fbinsize * sum(Pxx(indl));
hfp   = fbinsize * sum(Pxx(indh));
vlfp  = fbinsize * sum(Pxx(indv));
out.lfhf  = lfp / hfp;
total     = fbinsize * sum(Pxx);
out.vlf   = vlfp/total * 100;
out.lf    = lfp/total * 100;
out.hf    = hfp/total * 100;

% ------------------------------------------------------------------------------
% Triangular histogram index
% ------------------------------------------------------------------------------
numBins = 10;
out.tri = length(y)/max(histcounts(y,numBins));

% ------------------------------------------------------------------------------
% Poincare plot measures:
% ------------------------------------------------------------------------------
% cf. "Do Existing Measures ... ", Brennan et. al. (2001), IEEE Trans Biomed Eng 48(11)
rmssd = std(diffy); % std of differenced series
sigma = std(y); % should be 1 for zscored time series
out.SD1 = 1/sqrt(2) * rmssd * 1000;
out.SD2 = sqrt(2 * sigma^2 - (1/2) * rmssd^2) * 1000;

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [meg] = read_ctf_meg4(fname, hdr, begsample, endsample, chanindx)

% READ_CTF_MEG4 reads specified samples from a CTF continous datafile
% It neglects all trial boundaries as if the data was acquired in
% non-continous mode.
%
% [meg] = read_ctf_meg4(filename, hdr, begsample, endsample)
%
% Where
%   filename	name of the datafile, including the .meg4 extension
%   header      with all data information (from read_ctf_meg4)
%   begsample   index of the first sample to read
%   endsample   index of the last sample to read
%   chanindx	index of channels to read (optional, default is all)
%
% Due to a non-dislosure agreement between CTF and the F.C. Donders Centre,
% the allowed use of this function is limited. Do not distribute this function.

% This program is provided to users of CTF MEG systems as a courtesy only. Please
% do not redistribute it without permission from CTF Systems Inc.
% This program has no warranty whatsoever.

% Copyright (C) 2003, Robert Oostenveld
%
% use global flag for feedback
global fb
if isempty(fb)
  fb = 0;
end

if begsample<1
  error('cannot read before the start of the data');
end

if endsample>hdr.nSamples*hdr.nChans*hdr.nTrials
  error('cannot read beyond the end of the data');
end

if begsample>endsample
  error('cannot read a negative number of samples');
end

if nargin<5
  % select all channels
  chanindx = 1:hdr.nChans;
end

if isempty(chanindx)
  error('no channels were specified for reading CTF data')
end

fid = fopen(fname,'r','ieee-be');

if fid == -1
  error('could not open datafile');
end

CTFformat = setstr(fread(fid,8,'char'))';
if (strcmp(CTFformat(1,1:7),'MEG41CP')==0),
    error('datafile is not in CTF MEG4 format')
end

% the data is not channel multiplexed, but stored in trials
% in each trial first all samples for channel 1 are given, then all samples of channel 2 ...
% this means that we have to read each channel for each trial

begtrial = ceil(begsample/hdr.nSamples);
endtrial = ceil(endsample/hdr.nSamples);

% this counts the files and offset if the 2GB file size boundary is encountered
multifilecount  = 0;
multifileoffset = 0;
fseek(fid, 0, 'eof');
multifilelength = ftell(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read only the selected data, channel-wise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if endtrial==begtrial
  rawbegsample = begsample - (begtrial-1)*hdr.nSamples;
  rawendsample = endsample - (begtrial-1)*hdr.nSamples;
  for chan=1:length(chanindx)
    % jump to the begin of this channel in this trial
    channeloffset = 8*(multifilecount+1)+(begtrial-1)*hdr.nSamples*4*hdr.nChans+(chanindx(chan)-1)*hdr.nSamples*4-multifileoffset;
    if channeloffset>=multifilelength
      % data goes beyond 2GB file boundary, jump to the next file
      channeloffset   = channeloffset - multifilelength+8;	% change the current offset, keep the header
      multifileoffset = multifileoffset + multifilelength;	% remember for the next trial
      multifilecount  = multifilecount + 1;			% increase the file counter
      nextname = sprintf('%s.%d_meg4', fname(1:(end-5)), multifilecount);
      fprintf('data goes beyond 2GB file boundary, continuing with %s\n', nextname);
      fclose(fid);
      fid = fopen(nextname,'r','ieee-be');
      fseek(fid, 0, 'eof');
      multifilelength = ftell(fid);				% determine the length of the current file
    end
    fseek(fid, channeloffset, 'bof');
    % jump to the first sample of interest
    fseek(fid, 4*(rawbegsample-1), 'cof');
    % read the data from this channel
    [tmp, count] = fread(fid,[(endsample-begsample+1),1],'int32');
    raw(:,chan) = tmp;
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read and concatenate the raw data of all trials
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else
  raw = zeros((endtrial-begtrial+1)*hdr.nSamples, length(chanindx));
  for trial=begtrial:endtrial

    if length(chanindx)==hdr.nChans & all(chanindx(:)'==1:hdr.nChans)
      % read the data from all channels
      rawbegsample = (trial-begtrial)*hdr.nSamples + 1;
      rawendsample = (trial-begtrial)*hdr.nSamples + hdr.nSamples;
      % jump to the begin of this trial
      trialoffset = 8*(multifilecount+1)+(trial-1)*hdr.nSamples*4*hdr.nChans-multifileoffset;
      if trialoffset>=multifilelength
        % data goes beyond 2GB file boundary, jump to the next file
        trialoffset   = trialoffset - multifilelength+8;	% change the current offset, keep the header
        multifileoffset = multifileoffset + multifilelength;	% remember for the next trial
        multifilecount  = multifilecount + 1;			% increase the file counter
        nextname = sprintf('%s.%d_meg4', fname(1:(end-5)), multifilecount);
        fprintf('data goes beyond 2GB file boundary, continuing with %s\n', nextname);
        fclose(fid);
        fid = fopen(nextname,'r','ieee-be');
        fseek(fid, 0, 'eof');
        multifilelength = ftell(fid);				% determine the length of the current file
      end
      fseek(fid, trialoffset, 'bof');
      [tmp, count] = fread(fid,[hdr.nSamples,hdr.nChans],'int32');
      raw(rawbegsample:rawendsample,:) = tmp;

    else
      % read the data from the selected channels
      rawbegsample = (trial-begtrial)*hdr.nSamples + 1;
      rawendsample = (trial-begtrial)*hdr.nSamples + hdr.nSamples;
      for chan=1:length(chanindx)
        % jump to the begin of this channel in this trial
        channeloffset = 8*(multifilecount+1)+(trial-1)*hdr.nSamples*4*hdr.nChans+(chanindx(chan)-1)*hdr.nSamples*4-multifileoffset;
        if channeloffset>=multifilelength
          % data goes beyond 2GB file boundary, jump to the next file
          channeloffset   = channeloffset - multifilelength+8;	% change the current offset, keep the header
          multifileoffset = multifileoffset + multifilelength;	% remember for the next trial
          multifilecount  = multifilecount + 1;			% increase the file counter
          nextname = sprintf('%s.%d_meg4', fname(1:(end-5)), multifilecount);
          fprintf('data goes beyond 2GB file boundary, continuing with %s\n', nextname);
          fclose(fid);
          fid = fopen(nextname,'r','ieee-be');
          fseek(fid, 0, 'eof');
          multifilelength = ftell(fid);				% determine the length of the current file
        end
        fseek(fid, channeloffset, 'bof');
        [tmp, count] = fread(fid,[hdr.nSamples,1],'int32');
        raw(rawbegsample:rawendsample,chan) = tmp;
      end %chan

    end %read all channels
  end %trial

  % select the raw data corresponding to the samples of interest
  rawoffset    = (begtrial-1)*hdr.nSamples;
  rawbegsample = begsample - rawoffset;
  rawendsample = endsample - rawoffset;
  raw = raw(rawbegsample:rawendsample, :);
end
fclose(fid);

% multiply the dimensionless values with the calibration value
gain = hdr.gainV(chanindx);	% only for selected channels
meg = raw';			% transpose the raw data
for i=1:size(meg,1)
  meg(i,:) = gain(i)*meg(i,:);
end

author: PatternRecognition

repo: OpenBMI


text:
% HealpixBorderGetNext for the north polar cap area
% flag = 0 : where dTheta_dPhi >= 0
% flag = 1 : where dTheta_dPhi <= 0
function [phi, theta, dTheta_dPhi] = HealpixBorderGetNextPC(n, k, prev_phi, delta_phi, flag)

phi = prev_phi + delta_phi;
a = 1;
b = - k^2 * pi^2 / (12 * n^2);

if flag == 0
    theta = acos(a + b / phi^2);
    dTheta_dPhi = 2 * b / (phi^3* sin(theta));
else
    theta = acos(a + b / (phi - pi/2)^2);
    dTheta_dPhi = 2 * b / ((phi - pi/2)^3 * sin(theta));
end

author: mtex-toolbox

repo: mtex

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


text:
function d = det(A)
    % Symbolic matrix determinant.
    %   DET(A) is the determinant of the symbolic matrix A.
    %
    %   Examples:
    %       det([a b;c d]) is a*d-b*c.

    % Convert inputs to SymExpression
    % A = SymExpression(A);

    % check if A is a square matrix
    ret = eval_math(['SquareMatrixQ[' A.s ']']);

    if strcmp('False',ret)
        error(message('symbolic:det:SquareMatrix'));
    end

    % construct the operation string
    sstr = ['Det[' A.s ']'];

    % create a new object with the evaluated string
    d = SymExpression(sstr);
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function test_issue1214

% WALLTIME 00:10:00
% MEM 1gb
% DEPENDENCY ft_singleplotER ft_multiplotER ft_plot_vector

%% make some data

timelock1 = ft_timelockanalysis([], ft_freqsimulation([]));

cfg = [];
cfg.parameter = 'avg';
cfg.operation = '-1*x1';
timelock2 = ft_math(cfg, timelock1);

%% this works

cfg = [];
cfg.channel = 'mix';
figure; ft_singleplotER(cfg, timelock1, timelock2);

%% this works

cfg = [];
cfg.channel = 'mix';
cfg.graphcolor = 'my';
figure; ft_singleplotER(cfg, timelock1, timelock2);

%% this failed

cfg = [];
cfg.channel = 'mix';
cfg.graphcolor = [
  1 0 0
  0 1 0
  ];
figure; ft_singleplotER(cfg, timelock1, timelock2);

%% this failed

cfg = [];
cfg.layout = 'vertical';
cfg.graphcolor = [
  1 0 0
  0 1 0
  ];
figure; ft_multiplotER(cfg, timelock1, timelock2);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [x,y,z,info]=csdp(At,b,c,K,pars)
% CSDP Uses CSDP to solve a problem in SeDuMi format
%
% [x,y,z,info]=csdp(At,b,c,K,pars)
%
% Uses CSDP to solve a problem in SeDuMi format.
%
% Input:
%        At, b, c, K      SDP problem in SeDuMi format.
%        pars             CSDP parameters (optional parameter.)
%
% Output:
%
%        x, y, z          solution.
%        info             CSDP return code.
%                         info=100 indicates a failure in the MATLAB
%                         interface, such as inability to write to
%                         a temporary file.
%
% Note: This interface makes use of temporary files with names given by the
% tempname function.  This will fail if there is no working temporary
% directory or there isn't enough space available in this directory.
%
% Note: This code writes its own param.csdp file in the current working
% directory.  Any param.csdp file already in the directory will be deleted.
%
% Note: It is assumed that csdp is the search path made available through
% the "system" or "dos" command.  Typically, having the csdp executable in
% current working directory will work, although some paranoid system
% administrators keep . out of the path.  In that case, you'll need to
% install csdp in one of the directories that is in the search path.
% A simple test is to run csdp from a command line prompt.
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    %
    % First, put a dummy pars in place if no argument was given.  Also
    % set pars.printlevel if not given.
    %
    if (nargin<5)
        pars.printlevel=1;
    else
        if (isfield(pars,'printlevel'))
            pars.printlevel=pars.printlevel;
        else
            pars.printlevel=1;
        end
    end
    %
    % Write out the param.csdp file.
    %
    fid=fopen('param.csdp','w');
    if (fid==-1)
        if (pars.printlevel ~= 0)
            fprintf('Could not open param.csdp\n');
        end
        info=100;
        return
    end
    %
    % Now, go through the parameters.
    %

    if (isfield(pars,'axtol'))
        fprintf(fid,'axtol= %e\n',pars.axtol);
    else
        fprintf(fid,'axtol=%e\n',1.0e-8);
    end

    if (isfield(pars,'atytol'))
        fprintf(fid,'atytol= %e\n',pars.atytol);
    else
        fprintf(fid,'atytol=%e\n',1.0e-8);
    end

    if (isfield(pars,'objtol'))
        fprintf(fid,'objtol= %e\n',pars.objtol);
    else
        fprintf(fid,'objtol=%e\n',1.0e-8);
    end

    if (isfield(pars,'pinftol'))
        fprintf(fid,'pinftol= %e\n',pars.pinftol);
    else
        fprintf(fid,'pinftol=%e\n',1.0e8);
    end

    if (isfield(pars,'dinftol'))
        fprintf(fid,'dinftol= %e\n',pars.dinftol);
    else
        fprintf(fid,'dinftol=%e\n',1.0e8);
    end

    if (isfield(pars,'maxiter'))
        fprintf(fid,'maxiter= %d\n',pars.maxiter);
    else
        fprintf(fid,'maxiter=%d\n',100);
    end
    if (isfield(pars,'minstepfrac'))
        fprintf(fid,'minstepfrac= %e\n',pars.minstepfrac);
    else
        fprintf(fid,'minstepfrac=%e\n',0.90);
    end

    if (isfield(pars,'maxstepfrac'))
        fprintf(fid,'maxstepfrac= %e\n',pars.maxstepfrac);
    else
        fprintf(fid,'maxstepfrac=%e\n',0.97);
    end

    if (isfield(pars,'minstepp'))
        fprintf(fid,'minstepp= %e\n',pars.minstepp);
    else
        fprintf(fid,'minstepp=%e\n',1.0e-8);
    end

    if (isfield(pars,'minstepd'))
        fprintf(fid,'minstepd= %e\n',pars.minstepd);
    else
        fprintf(fid,'minstepd=%e\n',1.0e-8);
    end

    if (isfield(pars,'usexzgap'))
        fprintf(fid,'usexzgap= %d\n',pars.usexzgap);
    else
        fprintf(fid,'usexzgap=%d\n',1);
    end

    if (isfield(pars,'tweakgap'))
        fprintf(fid,'tweakgap= %d\n',pars.tweakgap);
    else
        fprintf(fid,'tweakgap=%d\n',0);
    end

    if (isfield(pars,'affine'))
        fprintf(fid,'affine= %d\n',pars.affine);
    else
        fprintf(fid,'affine=%d\n',0);
    end

    if (isfield(pars,'printlevel'))
        fprintf(fid,'printlevel= %d\n',pars.printlevel);
    else
        fprintf(fid,'printlevel=%d\n',1);
    end

    if (isfield(pars,'perturbobj'))
        fprintf(fid,'printlevel= %d\n',pars.perturbobj);
    else
        fprintf(fid,'printlevel=%d\n',1);
    end

    if (isfield(pars,'fastmode'))
        fprintf(fid,'printlevel= %d\n',pars.fastmode);
    else
        fprintf(fid,'printlevel=%d\n',0);
    end

    %
    % close the parameter file.
    %
    fclose(fid);
    %
    % Write the problem out.
    %
    fname=tempname;
    ret=writesdpa([fname '.dat-s'],At,b,c,K,pars);
    if (ret==1)
        info=100;
        x=[];
        y=[];
        z=[];
        return
    end
    %
    % Solve the problem.
    %
    directory = fileparts(which('csdp'));
    if any(strcmp(computer, {'PCWIN', 'PCWIN64'}))
        % Call for Windows
        info = dos(['"' directory '\csdp.exe" ' fname '.dat-s' ' ' fname '.sol'],'-echo');
    elseif any(strcmp(computer, {'GLNX86', 'GLNXA64'})) || strcmp(computer, 'GLNXA64')
        % Call for Linux
        info = system([directory '/csdplinux ' fname '.dat-s' ' ' fname '.sol']);
    elseif any(strcmp(computer, {'MACI', 'MACI64'}))
        % Call for Intel Mac
        info = system([directory '/csdpmaci ' fname '.dat-s' ' ' fname '.sol']);
    elseif strcmp(computer, 'MAC')
        % Call for Power PC Mac
        info = system([directory '/csdpmac ' fname '.dat-s' ' ' fname '.sol']);
    else
        error('There is no binary of CSDP available for your platform. You should obtain this binary yourself and edit techniques/csdp.m to solve this problem.');
    end
    %
    % Read back the solution.
    %
    [x,y,z]=readsol([fname '.sol'],K,length(b));
    %
    % Delete the temporary files, including param.csdp if we wrote one!
    %
    try
        delete([fname '.dat-s']);
        delete([fname '.sol']);
        delete('param.csdp');
    catch
        warning('Some temporary files may not be deleted.');
    end

repo: cluster


text:
function [trl, event] = ft_trialfun_imotions(cfg)

% FT_TRIALFUN_IMOTIONS makes a trial definition for an iMotions event structure
%
% The trialdef structure can contain the following specifications
% cfg.trialdef.eventtype  = string or cell-array of strings (default = 'StimulusName')
% cfg.trialdef.eventvalue = string or cell-array of strings (default = [])
%
% See also FT_DEFINETRIAL, FT_TRIALFUN_GENERAL

assert(isfield(cfg, 'event'), 'the configuration should contain the event structure');
assert(isfield(cfg, 'fsample'), 'the configuration should contain the sampling frequency');

% set the defaults
cfg.trialdef            = ft_getopt(cfg, 'trialdef', []);
cfg.trialdef.eventtype  = ft_getopt(cfg.trialdef, 'eventtype', 'StimulusName');
cfg.trialdef.eventvalue = ft_getopt(cfg.trialdef, 'eventvalue', []);
cfg.trialdef.offset     = ft_getopt(cfg.trialdef, 'offset', 'absolute'); % absolute or relative

event = cfg.event;

% start by selecting all events
sel = true(size(event));
% make a subselection of events
if ~isempty(cfg.trialdef.eventtype)
  sel = sel & strcmp({event.type}, cfg.trialdef.eventtype);
end
if ~isempty(cfg.trialdef.eventvalue)
  sel = sel & strcmp({event.value}, cfg.trialdef.eventvalue);
end

sample   = [event(sel).sample];
duration = [event(sel).duration];
type     = {event(sel).type};
value    = {event(sel).value};

begsample = sample(:);
endsample = sample(:)+duration(:)-1;

switch cfg.trialdef.offset
  case 'relative'
    % start of each trial/segment is t=0
    offset    = zeros(size(begsample));
  case 'absolute'
    % start of recording is t=0
    offset    = begsample - 1;
end

type  = type(:);
value = value(:);
trl   = table(begsample, endsample, offset, type, value);

author: spm


text:
%GOBJ_EXTRUDE Extrude face/geometry object.
%
%   [ GOBJ, STAT ] = GOBJ_EXTRUDE( GOBJ, D, V, FACE_WPL )
%   Extrude geometry object or face from GOBJ a distance D in the
%   direction of the vector V. FACE_WPL can in 3D indicate a face to
%   extrude or in 2D is a workplane struct defined by (p, n, t).

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
%File Description:
%  Unknown parameters in accelerometer error model are calculated using the
%  lm optimization algorithm.
close all
clc
clear

load AccRaw
%@ g = 9.8;
m = length(AccRaw);

y_dat = g*ones(m,1);  % expected gravitational acceleration
p0 = [1 1 1 0 0 0]';
p_init = [1.0 1.0 1.0 0.1 0.1 0.1]';  % initial value of unknown parameters

y_raw = calFunc(AccRaw, p0);  % gravitational acceleration measured by accelerometer
y_raw = y_raw(:);
r_raw = y_dat - y_raw;
p_fit = lm('calFunc', p_init, AccRaw, y_dat);
y_lm = calFunc(AccRaw,p_fit);  % gravitational acceleration measured by calibrated accelerometer
y_lm = y_lm(:);
r_lm = y_dat - y_lm;

kx = p_fit(1);
ky = p_fit(2);
kz = p_fit(3);
bx = p_fit(4);
by = p_fit(5);
bz = p_fit(6);

Ka1=[kx 0 0;0 ky 0;0 0 kz]
ba1=[bx by bz]'
save('calP1','Ka1','ba1')

figure
bar([r_raw'*r_raw, r_lm'*r_lm])
grid on;
set(gca,'XTickLabel', {'raw','lm'});
ylabel('fa');

t=1:m;
figure
title('Accelerometer Calibration')
plot(t, r_raw, t, r_lm)
legend('Uncalibrated', 'Calibrated-LM')
xlabel('Sampling sequence')
ylabel('Residual')

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
function [caNodeIndices, vResolution_] = ex_CreateIndexCatalog3D(mCatalog, mPolygon, bMap, nGriddingMode, fSmpValue, fSmpBnd, fSizeRectHorizontal, fSizeRectDepth)
    % Creates a cell-array with subcatalogs for every grid node defined by mPolygon.
    %
    % [caNodeIndices] = ex_CreateIndexCatalog(mCatalog, mPolygon, bMap, nGriddingMode,
    %                                                  nNumberEvents, fRadius, fSizeRectHorizontal, fSizeRectDepth)
    % -------------------------------------------------------------------------------------------------------------
    % Creates a cell-array with subcatalogs for every grid node defined by mPolygon. These subcatalogs
    %   contain only indices to the earthquake "rows" in mCatalog.
    %
    % Input parameters:
    %   mCatalog              Earthquake catalog
    %   mPolygon              Polygon (defined by ex_selectgrid)
    %   bMap                  Calculate cell-array for a map (true) or a cross-section (false)
    %   nGriddingMode         Mode of creating grid node subcatalogs
    %                         0: Constant number of events
    %                         1: Constant radius
    %                         2: Rectangular grid node samples
    %                         3: Spherical grid node samples with constant no
    %                            of events
    %                         4: Spherical grid node samples with constant radius
    %   nNumberEvents         Number of events per grid node (nGriddingMode == 0)
    %   fRadius               Radius of grid node sample (nGriddingMode == 1)
    %   fSizeRectHorizontal   Latitude/horizontal size of rectangle (nGriddingMode == 2)
    %   fSizeRectDepth        Longitude/depth size of rectangle (nGriddingMode == 2)
    %
    % Output parameters:
    %   caNodeIndices         Cell-array with index-catalogs per grid node of mPolygon
    %
    % Danijel Schorlemmer
    % June 17, 2002

    report_this_filefun();

    % Create the catalogs for each node with pointers to the overall catalog
    nNumberNodes_ = length(mPolygon(:,1));
    caNodeIndices = cell(nNumberNodes_, 1);

    % If cross-section calculate the length along cross-section
    if ~bMap
        nRow_ = mCatalog.Count;
        vXSecX_ = mCatalog(:,nColumn_);  % length along x-section
        vXSecY_ = (-1) * mCatalog.Depth;  % depth of hypocenters
    end

    % vResolution give the radius (for nGriddingMode = 4) and no. of events
    % (for nGriddingMode = 3)
    vResolution_(:,1)=nan(nNumberNodes_,1)
    if ((nGriddingMode == 3) || (nGriddingMode ==4)) clear vResolution_; end

    % Loop over all points of the polygon
    for nNode_ = 1:nNumberNodes_
        % Get the grid node coordinates
        fX_ = mPolygon(nNode_, 1);
        fY_ = mPolygon(nNode_, 2);
        if size(mPolygon,2)==3
            fZ_ = mPolygon(nNode_, 3);
        end

        if (nGriddingMode == 0) | (nGriddingMode == 1)  % Fixed radius or fixed number
            % Calculate distance from center point
            if bMap
                vDistances_ = sqrt(((mCatalog.Longitude-fX_)*cosd(fY_)*111).^2 + ((mCatalog.Latitude-fY_)*111).^2);
            else
                vDistances_ = sqrt(((vXSecX_ - fX_)).^2 + ((vXSecY_ - fY_)).^2);
            end
            if nGriddingMode == 0 % Fixed number
                if mCatalog.Count == 0
                    caNodeIndices{nNode_} = [];
                    % NaN for no events
                    vResolution_(nNode_) = nan;
                elseif nNumberEvents > mCatalog.Count
                    caNodeIndices{nNode_} = vIndices(1:mCatalog.Count);
                    % take the maximal distance for all eq. in the catalog
                    vResolution_(nNode_) = max(vdistances_);
                else
                    % Use first nNumberEvents events
                    [vTmp, vIndices] = sort(vDistances_);
                    caNodeIndices{nNode_} = vIndices(1:nNumberEvents);
                    % radius of the nNumberEvents-th event in the sorted vDistances_
                    vResolution_(nNode_) = vTmp(nNumberEvents);
                end
            else % Fixed radius
                % Use all events within fRadius
                caNodeIndices{nNode_} = find(vDistances_ <= fRadius);
                vResolution_(nNode_) = length(find(vDistances_ <= fRadius));
            end
        elseif nGriddingMode ==2  % Rectangular gridding (nGriddingMode == 2)
            if bMap
                vSel_ = ((mCatalog.Longitude >= (fX_ - fSizeRectHorizontal/2)) & (mCatalog.Longitude < (fX_ + fSizeRectHorizontal/2)) & ...
                    (mCatalog.Latitude >= (fY_ - fSizeRectDepth/2)) & (mCatalog.Latitude < (fY_ + fSizeRectDepth/2)));
                vResolution_(nNode_) = length(find(vSel_ > 0))
            else
                vSel_ = ((vXSecX_ >= (fX_ - fSizeRectHorizontal/2)) & (vXSecX_ < (fX_ + fSizeRectHorizontal/2)) & ...
                    (vXSecY_ >= (fY_ - fSizRectDepth/2)) & (vXSecY_ < (fY_ + fSizeRectDepth/2)));
                vResolution_(nNode_) = length(find(vSel_ > 0))
            end
            caNodeIndices{nNode_} = find(vSel_ == 1);
        elseif ((nGriddingMode == 3) || (nGriddingMode == 4))     % Spherical grid node samples with constant no of events
            vDistances_ = sqrt(((mCatalog.Longitude-fX_)*cosd(fY_)*111).^2 + ((mCatalog.Latitude-fY_)*111).^2 + (mCatalog.Depth-fZ_).^2);
            if nGriddingMode == 3  % Spherical grid node samples with constant no of events
                nNumberEvents=fSmpValue;
                fRadius=fSmpBnd;
                if mCatalog.Count == 0
                    caNodeIndices{nNode_} = [];
                    % NaN for no events
                    vResolution_{nNode_} = nan;
                    %               vResolution_(nNode_) = nan;
                elseif nNumberEvents > mCatalog.Count
                    caNodeIndices{nNode_} = vIndices(1:mCatalog.Count);
                    % take the maximal distance for all eq. in the catalog
                    vResolution_{nNode_} = vdistances_;
                    %               vResolution_(nNode_) = max(vdistances_);
                else
                    % Use first nNumberEvents events
                    [vTmp, vIndices] = sort(vDistances_);
                    caNodeIndices{nNode_} = vIndices(1:nNumberEvents);
                    % radius of the nNumberEvents-th event in the sorted vDistances_
                    vResolution_{nNode_} = vTmp(nNumberEvents);
                    %               vResolution_(nNode_) = vTmp(nNumberEvents);
                end
            elseif nGriddingMode == 4     % Spherical grid node samples with constant radius
                nNumberEvents=fSmpBnd;
                fRadius=fSmpValue;
                %  muss noch gemacht werden....

            end
        end
    end % of for nNode_
end

author: CelsoReyes


text:
function img = createImage(size, type, varargin)
%CREATEIMAGE Create a new image with given size and type
%
%   Deprecated, use 'imCreate' instead.
%
%   Usage
%   IMG = createImage(SIZE, TYPE)
%   IMG = createImage(SIZE, TYPE, VALUE)
%
%   Description
%   IMG = createImage(SIZE, TYPE)
%   SIZE is a row vector containing the size of the result image, and TYPE
%   is a string representating the class of result image.
%   The effect is similar to the function 'zeros' or 'false', but attends
%   to provide a unified method to create image.
%
%   Example
%   % create uint8 image
%   img = createImage([256 256], 'uint8');
%   imshow(img);
%
%   % create binary image
%   img = createImage([256 256], 'logical');
%   imshow(img);
%
%   % create an image with same size and type as another image
%   baseImage = imread('cameraman.tif');
%   img = createImage(size(baseImage), class(baseImage));
%   imshow(img);
%
%   % create uint8 image filled with dark gray pixels
%   img = createImage([256 256], 'uint8', 100);
%   imshow(img);
%
%   See also
%   zeros, ones, false
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-07-19,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

warning('malimpa:imFilters:createImage:deprecated', ...
    'function createImage is obsolete, use imCreate instead');

% allocate memory depending on type
if strcmp(type, 'logical')
    img = false(size);
else
    img = zeros(size, type);
end

% initialiaze with given value
if ~isempty(varargin)
    img(:) = varargin{1};
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function lam = line_adj_eigenvalues ( n )

%*****************************************************************************80
%
%% LINE_ADJ_EIGENVALUES returns the eigenvalues of the LINE_ADJ matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real LAM(N), the eigenvalues.
%
  lam = zeros ( n, 1 );

  for i = 1 : n
    angle = i * pi / ( n + 1 );
    lam(i) = 2.0 * cos ( angle );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ sym ] = search_symmetryset( candidate_list, obblist, Rmat, Dmat, LABELLEN )

objnum = size(obblist,1);
sym = {};
height_thre = 0.1*7;

%% search for the symmetry without central object
for i = 1:length(candidate_list)
    clist = candidate_list{i};
    clabellist = zeros(size(clist));
    csublists = {};
    for j = 1:length(clist)
        label = obblist(clist(j),end-LABELLEN+1:end);
        label = binary2dec(label);
        clabellist(j) = label;
    end
    ls = unique(clabellist);
    for j = 1:length(ls)
        label = ls(j);
        idx = find(clabellist==label);
        csublists{length(csublists)+1} = clist(idx);
        % "csublists" stores the sublists with same labels (and same size)
    end

    for j = 1:length(csublists)
        cslist = csublists{j};
        if(length(cslist)<=1)
            continue;
        end
        dists = Dmat(cslist, cslist);
        symdist_thre = 0.05*7;

        dl = unique(dists);
        dl = sort(dl);

        ismatched = zeros(1,length(cslist));
        for k = 1:length(dl)
            mind = dl(k);

            cssublists = {};
            for p = 1:length(cslist)
                for q = 1:length(cslist)
                    if(p~=q&&abs(Dmat(cslist(p),cslist(q))-mind)<symdist_thre)
                        isadded = -1;
                        for v = 1:size(cssublists,1)
                            if(ismember(cslist(p),cssublists{v,1})&&isadded==-1)
                                if(ismatched(q)~=1)
                                    if(abs(obblist(cslist(p),2)-obblist(cslist(q),2))<height_thre)
                                        cssublists{v,1} = union(cssublists{v,1},cslist(q));
                                        ismatched(q)=1;
                                        isadded = 1;
                                        break;
                                    end
                                end
                            end
                            if(ismember(cslist(q),cssublists{v,1})&&isadded==-1)
                                if(ismatched(p)~=1)
                                    if(abs(obblist(cslist(p),2)-obblist(cslist(q),2))<height_thre)
                                        cssublists{v,1} = union(cssublists{v,1},cslist(p));
                                        ismatched(p)=1;
                                        isadded = 1;
                                        break;
                                    end
                                end
                            end
                        end
                        if(isadded==-1)
                            if(ismatched(p)~=1&&ismatched(q)~=1)
                                if(abs(obblist(cslist(p),2)-obblist(cslist(q),2))<height_thre)
                                    cssublists{size(cssublists,1)+1,1} = [cslist(p),cslist(q)];
                                    cssublists{size(cssublists,1),2} = mind;
                                    ismatched(p)=1;
                                    ismatched(q)=1;
                                end
                            end
                        end
                    end
                end
            end
            for v = 1:size(cssublists,1)
                csslist = cssublists{v,1};
                % add sym
                sym{size(sym,1)+1,1} = csslist;
                sym{size(sym,1),2} = cssublists{v,2};
            end
        end

    end
end

final_sym = {};
% one object can be involved in one symmetry relation
isinvolved = zeros(1,objnum);
for i = 1:objnum
    if(isinvolved(i)~=1)
        symlist = {};
        symdist = [];
        if(size(sym,1)>0)
            for j = 1:size(sym,1)
                if(ismember(i,sym{j,1}))
                    symlist{length(symlist)+1}= sym{j,1};
                    sdist = Dmat(i,sym{j,1});
                    maxd = max(sdist(:));
                    for k = 1:size(sdist,1)
                        sdist(k,k) = k+maxd;
                    end
                    symdist(length(symlist)) = min(sdist);
                end
            end
            if(length(symdist)>0)
                [md,I] = min(symdist);
                symlist = symlist{I};
            end
        end

        if(length(symlist)>1)
            final_sym{size(final_sym,1)+1,1} = symlist;
            final_sym{size(final_sym,1),2} = md;
            isinvolved(symlist) = 1;
        end
    end
end
sym = final_sym;

end

repo: GRAINS


text:
function ap=aperiodiccomp(apv,dpv,ashift,f0,nshift,imgi)
%	ap=aperiodiccomp(apv,dpv,ashift,f0,nshift,fftl,imgi);
%	Calculate aperiodicity index
%	Input parameters
%		apv, dpv : Upper and lower envelope
%		ashift		: shift step for aperiodicity index calculation (ms)
%		f0			: fundamental frequency (Hz)
%		nshift		: shift step for f0 information (ms)
%		fftl		: FFT size
%		imgi		: display indicator, 1: display on (default) 0: off

%   modified to add the waitbar on 08/Dec./2002
%	modified by Takahashi 10/Aug./2005
%	modified by Kawahara 10/Sept./2005

if nargin==5; imgi=1; end;
%[nn,mm]=size(nsgram);
mm=length(f0);
%%nn=fftl/2+1;
[~,m2]=size(apv);

x=(0:m2-1)'*ashift;
xi=(0:mm-1)'*nshift;
xi=min(max(x),xi);

if imgi==1; hpg=waitbar(0.1,'Interpolating periodicity information'); end;
if imgi==1; drawnow; end;
%ap=interp1q(x,(dpv-apv)',xi)';%,'*linear')';
ap = interp1(x, (dpv-apv)',xi, 'linear', 'extrap')';
if imgi==1; close(hpg); end;

author: HidekiKawahara


text:
%SE2 Create planar translation and rotation transformation
%
% T = SE2(X, Y, THETA) is an SE(2) homogeneous transformation (3x3)
% representing translation X and Y, and rotation THETA in the plane.
%
% T = SE2(XY) as above where XY=[X,Y] and rotation is zero
%
% T = SE2(XY, THETA) as above where XY=[X,Y]
%
% T = SE2(XYT) as above where XYT=[X,Y,THETA]
%
% See also TRANSL2, ROT2, ISHOMOG2, ISROT2, TRPLOT2.

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

function t = se2(a, b, c, varargin)

    opt.deg = false;

    opt = tb_optparse(opt, varargin);

    if length(a) == 3
        x = a(1);
        y = a(2);
        th = a(3);
    elseif length(a) == 2
        x = a(1);
        y = a(2);
        if nargin < 2
            th = 0;
        else
            th = b;
        end
    else
        x = a;
        y = b;
        if nargin < 3
            th = 0;
        else
            th = c;
        end
    end

    if opt.deg
        th = th * pi/180.0;
    end
    cth = cos(th);
    sth = sin(th);
    R = [cth -sth; sth cth];

    t = [R [x; y]; 0 0 1];

author: Allopart

repo: rbpf-gmapping


text:
function c = char(s,varargin)
% object -> string

if check_option(varargin,'verbose')
  c = option2str([{s.pointGroup},s.alignment]);
  if ~isempty(s.mineral), c = [s.mineral ' (' c ')']; end
elseif check_option(varargin,'latex')
  c = ['$' regexprep(s.pointGroup,'-(\w)','\\bar{$1}') '$'];
elseif check_option(varargin,'compact')
  if ~isempty(s.mineral)
    c = s.mineral;
  else
    c = s.pointGroup;
  end

  if ~getMTEXpref('generatingHelpMode')
    id = pushTemp(s);
    c = ['' c ''];
  end
else
  c = ['"',s.pointGroup,'"'];
end

if check_option(varargin,'symmetryType'), c = ['  crystal symmetry : ' c]; end


text:
% op_complexConj.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% out=op_complexConj(in)
%
% DESCRIPTION:
% take the complex conjugate of the data;
%
% INPUTS:
% in	= Input data in matlab structure format.
%
% OUTPUTS:
% out   = Output following conjugation.

function out=op_complexConj(in);

fids=in.fids;
sz=size(fids);

fids=conj(fids);

%re-calculate Specs using fft
specs=fftshift(ifft(fids,[],in.dims.t),in.dims.t);

%Calculate t and ppm arrays using the calculated parameters:
f=[(-in.spectralwidth/2)+(in.spectralwidth/(2*sz(1))):in.spectralwidth/(sz(1)):(in.spectralwidth/2)-(in.spectralwidth/(2*sz(1)))];
%ppm=-f/(in.Bo*42.577);
% ppm=-f/(3*42.577);
ppm = -f/(in.Bo*in.gamma);
if strcmp(in.nucleus,'1H')
    ppm=ppm+4.65;
end

%t=[0:in.dwelltime:(sz(1)-1)*in.dwelltime];
t=[in.dwelltime:in.dwelltime:sz(1)*in.dwelltime];

%FILLING IN DATA STRUCTURE
out=in;
out.fids=fids;
out.specs=specs;
out.sz=sz;
out.ppm=ppm;
out.t=t;

%FILLING IN THE FLAGS
out.flags=in.flags;
out.flags.writtentostruct=1;

author: CIC-methods

repo: FID-A

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


text:
function [dA,twt]=diffrd3(fname)
% diffs an rd3 file... (requires the .rad file to be present)
% syntax: [dA,twt]=loadrd3('prof4')
%
% Aslak Grinsted 2002
%

% load details from .RAD file

fname=strrep(lower(fname),'.rd3','');

tracecount=inf;
starttrace=1;
stack=0;

infos=rd3info(fname,'scans','samples','timewindow');
sz=infos{2};
timewindow=infos{3};

twt=((1:sz)'-1)*timewindow*1e-9/sz;

%loads the rd3 file
fid=fopen([fname '.rd3'],'r');
if (starttrace>1)
   fseek(fid,(starttrace-1)*sz*2,0);
end

blocksize=1+fix(1000000/sz); %number of traces to read at a time...

lastTrace=[];
traceIdx=1;
dA=zeros(1,1);
tracesread=0;
while ((~feof(fid))&(tracesread<tracecount))
    B=fread(fid,[sz blocksize],'short');
    if isempty(lastTrace)
        lastTrace=B(:,1); %ugly hack! (but ok)
    end
    for ii=1:size(B,2)
        dA(traceIdx)=sum(abs(lastTrace-B(:,ii)).*twt); %maybe multiply by twt to do some correction for decay of amps...
        traceIdx=traceIdx+1;
        lastTrace=B(:,ii);
    end
    tracesread=tracesread+size(B,2);
end

dA=dA/mean(dA);

fclose(fid);

author: Sable

repo: mcbench-benchmarks


text:
function fv= proc_rSquareSigned(fv, varargin)
%PROC_RSQUARESIGNED - computes signed r^2 values (measure for discriminance)
%
%Synopsis:
% 	fv = proc_rSquareSigned(fv, <opt>)
%
%Returns:
% FV_RVAL - data structure of signed squared biserial correlation coefficients
%  .x     - signed squared biserial correlation between each featur and the
%           class label
%  .se    - contains the standard error of atanh(r), if opt.Stats==1
%  .p     - contains the p value of null hypothesis that there is zero
%           correlation between feature and class-label, if opt.Stats==1
%  .sgnlogp - contains the signed log10 p-value, if opt.Stats==1
%             if opt.Bonferroni==1, the p-value is multiplied by
%             fv_rval.corrfac
%  .sgnlogp - contains the signed log10 p-value, if opt.Stats==1
%           if opt.Bonferroni==1, the p-value is multiplied by
%           fv_rval.corrfac and then logarithmized
%  .sigmask - binary array indicating significance at alpha level
%             opt.Alphalevel, if opt.Stats==1 and opt.Alphalevel > 0
%  .corrfac - Bonferroni correction factor (number of simultaneous tests),
%             if opt.Bonferroni==1
%
%Properties:
% 'TolerateNans': observations with NaN value are skipped
%    (nanmean/nanstd are used instead of mean/std). Deafult: 0
% 'ValueForConst': constant feauture dimensions are assigned this
%    value. Default: NaN.
% 'MulticlassPolicy': possible options: 'pairwise' (default),
%    'all-against-last', 'each-against-rest', or provide specified
%    pairs as an [nPairs x 2] sized matrix. ('specified_pairs' is obsolete)
% 'Stats' - if true, additional statistics are calculated, including the
%           standard error of atanh(r), the p-value for the null
%           Hypothesis that the correlation is zero,
%           and the "signed log p-value"
% 'Bonferroni' - if true, Bonferroni corrected is used to adjust p-values
%                and their logarithms
% 'Alphalevel' - if provided, a binary indicator of the significance to the
%                alpha level is returned for each feature in fv_rval.sigmask
%
%Description:
% Computes the r^2 value for each feature, multiplied by the sign of
% of r value. The r^2 value is a measure
% of how much variance of the joint distribution can be explained by
% class membership.
%
% Example:
%  [cnt, mrk]= file_readBV(some_file);   %load EEG-data in BV-format
%  mrk= mrk_defineClasses(mrk, {1, 2; 'target','nontarget'});
%  epo= proc_segmentation(cnt, mrk, [-200 800], 'CLab', {'Fz','Cz','Pz'});
%  epo_r = proc_rSquareSigned(epo);
%
% See also proc_classmeanDiff, proc_rValues, proc_rSquare
%
% 03-03 Benjamin Blankertz
% 09-2012 stefan.haufe@tu-berlin.de

if nargin==0,
  fv=proc_rValues; return
end

fv= proc_rValues(fv, varargin{:});
fv.x= fv.x .* abs(fv.x);
for cc= 1:length(fv.className),
  fv.className{cc}= ['sgn r^2' fv.className{cc}(2:end)];
end
fv.yUnit= 'sgn r^2';

author: bbci


text:
function element_size = p11_element_size ( )

%*****************************************************************************80
%
%% P11_ELEMENT_SIZE returns a typical element size for problem 11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, real ELEMENT_SIZE, a typical element size.
%
  element_size = 0.15;

  return
end


text:
function [ tri_num, tri_vert, tri_nabe ] = r8tris2 ( point_num, p )

%*****************************************************************************80
%
%% R8TRIS2 constructs a Delaunay triangulation of 2D vertices.
%
%  Discussion:
%
%    The routine constructs the Delaunay triangulation of a set of 2D vertices
%    using an incremental approach and diagonal edge swaps.  Vertices are
%    first sorted in lexicographically increasing (X,Y) order, and
%    then are inserted one at a time from outside the convex hull.
%
%  Modified:
%
%    07 February 2005
%
%  Author:
%
%    Original FORTRAN77 version by Barry Joe.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Barry Joe,
%    GEOMPACK - a software package for the generation of meshes
%    using geometric algorithms,
%    Advances in Engineering Software,
%    Volume 13, pages 325-331, 1991.
%
%  Parameters:
%
%    Input, integer POINT_NUM, the number of vertices.
%
%    Input, real P(2,POINT_NUM), the vertices.
%
%    Output, integer TRI_NUM, the number of triangles in the triangulation;
%    TRI_NUM is equal to 2*POINT_NUM - NB - 2, where NB is the number
%    of boundary vertices.
%
%    Output, integer TRI_VERT(3,TRI_NUM), the nodes that make up each triangle.
%    The elements are indices of P.  The vertices of the triangles are
%    in counter clockwise order.
%
%    Output, integer TRI_NABE(3,TRI_NUM), the triangle neighbor list.
%    Positive elements are indices of TIL; negative elements are used for links
%    of a counter clockwise linked list of boundary edges; LINK = -(3*I + J-1)
%    where I, J = triangle, edge index; TRI_NABE(J,I) refers to
%    the neighbor along edge from vertex J to J+1 (mod 3).
%
  tri_num = 0;
  tri_vert = [];
  tri_nabe = [];

  tol = 100.0 * r8_epsilon ( );
%
%  Sort the vertices by increasing (x,y).
%
  indx = r82vec_sort_heap_index_a ( point_num, p );

  p = r82vec_permute ( point_num, p, indx );
%
%  Make sure that the data points are "reasonably" distinct.
%
  m1 = 1;

  for i = 2 : point_num

    m = m1;
    m1 = i;

    k = 0;

    for j = 1 : 2

      cmax = max ( abs ( p(j,m) ), abs ( p(j,m1) ) );

      if ( tol * ( cmax + 1.0 ) < abs ( p(j,m) - p(j,m1) ) )
        k = j;
        break
      end

    end

    if ( k == 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8TRIS2 - Fatal error!\n' );
      fprintf ( 1, '  Fails for point number I = %d\n', i );
      fprintf ( 1, '  M = %d\n', m );
      fprintf ( 1, '  M1 = %d\n', m1 );
      fprintf ( 1, '  X,Y(M)  = %f  %f\n', p(1,m), p(2,m) );
      fprintf ( 1, '  X,Y(M1) = %f  %f\n', p(1,m1), p(2,m1) );
      error ( 'R8TRIS2 - Fatal error!' )
      return
    end

  end
%
%  Starting from points M1 and M2, search for a third point M that
%  makes a "healthy" triangle (M1,M2,M)
%
  m1 = 1;
  m2 = 2;
  j = 3;

  while ( 1 )

    if ( point_num < j )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8TRIS2 - Fatal error!\n' );
      error ( 'R8TRIS2 - Fatal error!' )
      return
    end

    m = j;

    lr = lrline ( p(1,m), p(2,m), p(1,m1), p(2,m1), p(1,m2), p(2,m2), 0.0 );

    if ( lr ~= 0 )
      break
    end

    j = j + 1;

  end
%
%  Set up the triangle information for (M1,M2,M), and for any other
%  triangles you created because points were collinear with M1, M2.
%
  tri_num = j - 2;

  if ( lr == -1 )

    tri_vert(1,1) = m1;
    tri_vert(2,1) = m2;
    tri_vert(3,1) = m;
    tri_nabe(3,1) = -3;

    for i = 2 : tri_num

      m1 = m2;
      m2 = i+1;
      tri_vert(1,i) = m1;
      tri_vert(2,i) = m2;
      tri_vert(3,i) = m;
      tri_nabe(1,i-1) = -3 * i;
      tri_nabe(2,i-1) = i;
      tri_nabe(3,i) = i - 1;

    end

    tri_nabe(1,tri_num) = -3 * tri_num - 1;
    tri_nabe(2,tri_num) = -5;
    ledg = 2;
    ltri = tri_num;

  else

    tri_vert(1,1) = m2;
    tri_vert(2,1) = m1;
    tri_vert(3,1) = m;
    tri_nabe(1,1) = -4;

    for i = 2 : tri_num
      m1 = m2;
      m2 = i+1;
      tri_vert(1,i) = m2;
      tri_vert(2,i) = m1;
      tri_vert(3,i) = m;
      tri_nabe(3,i-1) = i;
      tri_nabe(1,i) = -3 * i - 3;
      tri_nabe(2,i) = i - 1;
    end

    tri_nabe(3,tri_num) = -3 * tri_num;
    tri_nabe(2,1) = -3 * tri_num - 2;
    ledg = 2;
    ltri = 1;

  end
%
%  Insert the vertices one at a time from outside the convex hull,
%  determine visible boundary edges, and apply diagonal edge swaps until
%  Delaunay triangulation of vertices (so far) is obtained.
%
  top = 0;

  for i = j+1 : point_num

    m = i;
    m1 = tri_vert(ledg,ltri);

    if ( ledg <= 2 )
      m2 = tri_vert(ledg+1,ltri);
    else
      m2 = tri_vert(1,ltri);
    end

    lr = lrline ( p(1,m), p(2,m), p(1,m1), p(2,m1), p(1,m2), p(2,m2), 0.0 );

    if ( 0 < lr )
      rtri = ltri;
      redg = ledg;
      ltri = 0;
    else
      l = -tri_nabe(ledg,ltri);
      rtri = floor ( l / 3 );
      redg = mod(l,3) + 1;
    end

    [ ltri, ledg, rtri, redg ] = vbedg ( p(1,m), p(2,m), point_num, p, ...
      tri_num, tri_vert, tri_nabe, ltri, ledg, rtri, redg );

    n = tri_num + 1;
    l = -tri_nabe(ledg,ltri);

    while ( 1 )

      t = floor ( l / 3 );
      e = mod ( l, 3 ) + 1;
      l = -tri_nabe(e,t);
      m2 = tri_vert(e,t);

      if ( e <= 2 )
        m1 = tri_vert(e+1,t);
      else
        m1 = tri_vert(1,t);
      end

      tri_num = tri_num + 1;
      tri_nabe(e,t) = tri_num;
      tri_vert(1,tri_num) = m1;
      tri_vert(2,tri_num) = m2;
      tri_vert(3,tri_num) = m;
      tri_nabe(1,tri_num) = t;
      tri_nabe(2,tri_num) = tri_num - 1;
      tri_nabe(3,tri_num) = tri_num + 1;
      top = top + 1;

      if ( point_num < top )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'R8TRIS2 - Fatal error!\n' );
        fprintf ( 1, '  Stack overflow.\n' );
        error ( 'R8TRIS2 - Fatal error!' )
      end

      work(top) = tri_num;

      if ( t == rtri & e == redg )
        break
      end

    end

    tri_nabe(ledg,ltri) = -3 * n - 1;
    tri_nabe(2,n) = -3 * tri_num - 2;
    tri_nabe(3,tri_num) = -l;
    ltri = n;
    ledg = 2;

    [ top, ltri, ledg, tri_vert, tri_nabe ] = swapec ( ...
      m, top, ltri, ledg, point_num, p, tri_num, tri_vert, tri_nabe, work );

  end
%
%  Now account for the sorting that we did.
%
  for i = 1 : 3
    for j = 1 : tri_num
      tri_vert(i,j) = indx ( tri_vert(i,j) );
    end
  end

  indx = perm_inverse ( point_num, indx );

  p = r82vec_permute ( point_num, p, indx );

  return
end


text:
function W = compute_skinning_weight(V,F,b)
%%%%%%%%%%
% TASK: implement (unbounded) biharmonic weights, which solves the
% following optimization problem:
%      min   0.5 * trace( W'*Q*W )
%      s.t.  W(b,:) = bc
%
% black list:
% - min_quad_with_fixed
%
% white list (hints)
% - "cotmatrix" for computing cotangent Laplacian
% - "massmatrix" for computing mass matrix
%
% Hint:
% - give a diagonal matrix A, you can invert A with
%   >> invA = diag(diag(A).^(-1));
% - constraints are 1 for a handle and 0 for other handls, so it is like
%   >> bc = eye(length(b));
% - you shouldn't need any other functions in gptoolbox
%%%%%%%%%%

% naive closest point (replace this one with your solution)
idx = knnsearch(V(b,:),V);
W = full(sparse([1:size(V,1)], idx, 1, size(V,1), length(b)));

end

author: odedstein

repo: sgi-introduction-course


text:
function value = r8vec_dif_f ( x )

%*****************************************************************************80
%
%% R8VEC_DIF_F evaluates the function used in R8VEC_DIF_TEST.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  value = exp ( x );

  return
end


text:
classdef (Sealed) DMDicomLibrary < DMDicomLibraryInterface
    % DMDICOMLIBRARY DicoMat implementation of DMDicomLibraryInterface for
    % reading Dicom data.
    %
    %
    %     Licence
    %     -------
    %     Part of DicoMat. https://github.com/tomdoel/dicomat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the BSD 3-Clause license. Please see the file LICENSE for details.
    %

    methods
        function isDicom = isdicom(obj, fileName)
            % Tests whether a file is in DICOM format

            isDicom = DMisdicom(fileName);
        end

        function metaheader = dicominfo(obj, varargin)
            % Reads the metaheader data from a Dicom file

            metaheader = DMdicominfo(varargin{:});
        end

        function imageData = dicomread(obj, fileName_or_metaHeader)
            % Reads the image data from a Dicom file

            imageData = DMdicomread(fileName_or_metaHeader);
        end

    end

    methods (Access = private)
        function obj = DMDicomLibrary
        end
    end

    methods (Static)
        function singleObj = getLibrary
            persistent singleton
            if isempty(singleton) || ~isvalid(singleton)
                singleton = DMDicomLibrary;
            end
            singleObj = singleton;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function clearStructMeshes(varargin)
%function clearStructMeshes()
%This function clears the surface meshes for passed structure uidC.
%
%APA, 06/21/07
%
% 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/>.

global planC
indexS = planC{end};

if isempty(varargin)
    structNumToClear = 1:length(planC{indexS.structures});
else
    uidC = varargin{1};
    structNumToClear = getAssociatedStr(uidC);
end

%Set Matlab path to directory containing the library
currDir = cd;

if ispc
    meshDir = fileparts(which('libMeshContour.dll'));
    cd(meshDir);
    loadlibrary('libMeshContour','MeshContour.h');
elseif isunix
    meshDir = fileparts(which('libMeshContour.so'));
    cd(meshDir);
    loadlibrary('libMeshContour.so','MeshContour.h');
end
%Clear existing structure meshes
waitbarH = waitbar(0,'Clearing surface meshes for structures...');
try
    for structNum = structNumToClear
        structUID = planC{indexS.structures}(structNum).strUID;
        calllib('libMeshContour','clear',structUID)
        planC{indexS.structures}(structNum).meshRep = 0;
        planC{indexS.structures}(structNum).meshS   = [];
        waitbar(structNum/length(planC{indexS.structures}),waitbarH)
    end
end
close(waitbarH)
%switch back the current irectory
cd(currDir)

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function pSum=polySum(p1,p2)
%%POLYSUM Add together two vectors representing 1D polynomials as power
%         series. The format of the polynomial vectors is the same as that
%         used in Matlab's polyval function: The coefficients are ordered
%         y(z)=p(end)+p(end-1)*z+p(end-2)*z^2+...
%
%INPUTS: p1,p2 Coefficients for the two power series. The number of terms
%              in the different series can be different. Thus p1 and p2 do
%              not have the same length. p1(1) is the coefficient of the
%              highest order term. p(end) is the additive constant term. p1
%              and p2 can be row or column vectors.
%
%OUTPUTS: pSum A column vector holding the sum of p1 and p2. The first
%              element corresponds to the highest term, as used in Matlab's
%              polyval function. Note that the vector is not shortened when
%              the sum of high order terms equals zero (terms cancel).
%
%Polynomials are summed by adding corresponding terms.
%
%As an example, consider adding
%x^2-3*x+2
%to
%x^5+3*x^3-4*x^2+1
%This is done as
% p1=[1;-3;2];
% p2=[1;0;3;-4;0;1];
% pSum=polySum(p1,p2)
%One will see that the result is pSum'=[1,0,3,-3,-3,3], which corresponds
%to x^5+3*x^3-3*x^2-3*x+3.
%
%August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

pSum=[zeros(length(p1)-length(p2),1);p2(:)]+[zeros(length(p2)-length(p1),1);p1(:)];

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function g = sigmoidGradient(z)
%SIGMOIDGRADIENT returns the gradient of the sigmoid function
%evaluated at z
%   g = SIGMOIDGRADIENT(z) computes the gradient of the sigmoid function
%   evaluated at z. This should work regardless if z is a matrix or a
%   vector. In particular, if z is a vector or matrix, you should return
%   the gradient for each element.

g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the gradient of the sigmoid function evaluated at
%               each value of z (z can be a matrix, vector or scalar).

g = exp(-z) ./ ((1.0 + exp(-z)) .^ 2);

% =============================================================

end

author: UtkarshPathrabe


text:
function varargout = sccn_xdf(filename, hdr, begsample, endsample, chanindx)

% This is a wrapper to the reading function from the XDF MATLAB toolbox.
%
% Use as
%   hdr = sccn_xdf(filename);
%   dat = sccn_xdf(filename, hdr, begsample, endsample, chanindx);
%   evt = sccn_xdf(filename, hdr);
%
% See also FT_FILETYPE, FT_READ_HEADER, FT_READ_DATA, FT_READ_EVENT, XDF2FIELDTRIP

% Copyright (C) 2019, Robert Oostenveld
%
% 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$

% ensure this is on the path
ft_hastoolbox('xdf', 1);

needhdr = (nargin==1);
needevt = (nargin==2);
needdat = (nargin==5);

streams = load_xdf(filename);

iscontinuous = false(size(streams));

% figure out which streams contain continuous/regular and discrete/irregular data
for i=1:numel(streams)

    % if the nominal srate is non-zero, the stream is considered continuous
    if ~strcmpi(streams{i}.info.nominal_srate, '0')

       iscontinuous(i) =  true;

       if ~isfield(streams{i}.info, 'effective_srate')

           % in case effective srate field value is missing, add one
           num_samples  = numel(streams{i}.time_stamps);
           t_begin      = streams{i}.time_stamps(1);
           t_end        = streams{i}.time_stamps(end);
           duration     = t_end - t_begin;
           streams{i}.info.effective_srate = (num_samples - 1) / duration;

       elseif isempty(streams{i}.info.effective_srate)

           % in case effective srate field value is missing, add one
           num_samples  = numel(streams{i}.time_stamps);
           t_begin      = streams{i}.time_stamps(1);
           t_end        = streams{i}.time_stamps(end);
           duration     = t_end - t_begin;
           streams{i}.info.effective_srate = (num_samples - 1) / duration;

       end

    end

end

% determine the stream with the highest sampling rate
srate = nan(size(streams));
for i=1:numel(streams)
  if iscontinuous(i)
    srate(i) = streams{i}.info.effective_srate;
  end
end
[dum, indx] = max(srate);

% only keep the stream with the maximum sampling rate
% this is probably the EEG stream
stream = streams{indx};

if needhdr
  % this section of code is shared with xdf2fieldtrip
  hdr = [];
  if isfield(stream.info, 'effective_srate')
    % the stream contains continuously sampled data
    hdr.Fs                  = stream.info.effective_srate;
    hdr.nSamplesPre         = 0;
    hdr.nSamples            = length(stream.time_stamps);
    hdr.nTrials             = 1;
    hdr.FirstTimeStamp      = stream.time_stamps(1);
    hdr.TimeStampPerSample  = (stream.time_stamps(end)-stream.time_stamps(1)) / (length(stream.time_stamps) - 1);
  else
    % the stream does not contain continuously sampled data
    hdr.Fs                  = NaN;
    hdr.nSamplesPre         = NaN;
    hdr.nSamples            = NaN;
    hdr.nTrials             = NaN;
    hdr.FirstTimeStamp      = NaN;
    hdr.TimeStampPerSample  = NaN;
  end
  if isfield(stream.info.desc, 'channels')
    hdr.nChans    = numel(stream.info.desc.channels.channel);
  else
    hdr.nChans    = str2double(stream.info.channel_count);
  end
  hdr.label       = cell(hdr.nChans, 1);
  hdr.chantype    = cell(hdr.nChans, 1);
  hdr.chanunit    = cell(hdr.nChans, 1);

  prefix = stream.info.name;
  for j=1:hdr.nChans
    if isfield(stream.info.desc, 'channels')
      hdr.label{j} = [prefix '_' stream.info.desc.channels.channel{j}.label];
      hdr.chantype{j} = stream.info.desc.channels.channel{j}.type;
      hdr.chanunit{j} = stream.info.desc.channels.channel{j}.unit;
    else
      % the stream does not contain continuously sampled data
      hdr.label{j} = num2str(j);
      hdr.chantype{j} = 'unknown';
      hdr.chanunit{j} = 'unknown';
    end
  end

  % keep the original header details
  hdr.orig = stream.info;

  % return the header
  varargout = {hdr};

elseif needevt
  streams = streams(~iscontinuous);
  event = [];
  for i=1:numel(streams)
    for j=1:numel(streams{i}.time_series)

      % convert the timestamps to the corresponding sample in the selected data stream
      % the first sample in the data stream corresponds to 1
      timestamp = streams{i}.time_stamps(j);
      sample = round((timestamp - hdr.FirstTimeStamp)/hdr.TimeStampPerSample) + 1;

      event(end+1).type      = streams{i}.info.type;
      event(end  ).value     = streams{i}.time_series{j};
      event(end  ).timestamp = timestamp;
      event(end  ).sample    = sample;
    end
  end

  % return the events
  varargout = {event};

elseif needdat
  dat = stream.time_series(chanindx, begsample:endsample);

  % return the data
  varargout = {dat};
end

author: fieldtrip


text:
function solution = optimalConservationVectors(S, lambda, delta)
% DC programming for solving the cardinality optimization problem
%
% .. math::
%
%    min  ~& \lambda ||x||_0 - \delta ||y||_0 \\
%    s.t. ~& x + S^T y = 0 \\
%         ~& 0 \leq y \leq 1e4
%
% USAGE:
%
%    solution = optimalConservationVectors(S, lambda, delta)
%
% INPUT:
%    S:           `m` x `n` stoichiometric matrix
%
% OPTIONAL INPUTS:
%    lambda:      default = 1
%    delta:       default = 1
%
% OUTPUT:
%    solution:    result of the solved cardinality optimization problem

if ~exist('lambda','var')
    lambda=1;
end
if ~exist('delta','var')
    delta=1;
end

[mlt,nlt]=size(S');
prob.p=mlt;
prob.q=nlt;
prob.r=0;
prob.c=zeros(nlt+mlt,1);
prob.A=[speye(mlt,mlt),S'];
prob.b=zeros(mlt,1);
prob.lb=[-inf*ones(mlt,1);zeros(nlt,1)];
prob.ub=[inf*ones(nlt,1);1e4*ones(mlt,1)];
prob.csense(1:mlt,1)='E';
params.lamda=lambda;
params.delta=delta;
solution = optimizeCardinality(prob,params);
% DC programming for solving the cardinality optimization problem
% The l0 norm is approximated by capped-l1 function.
% min       c'(x,y,z) + lambda*||x||_0 - delta*||y||_0
% s.t.      A*(x,y,z) <= b
%           l <= (x,y,z) <=u
%           x in R^p, y in R^q, z in R^r
%
% solution = optimizeCardinality(problem,params)
%
%  problem                  Structure containing the following fields describing the problem
%       p                   size of vector x
%       q                   size of vector y
%       r                   size of vector z
%       c                   (p+q+r) x 1 linear objective function vector
%       lambda              trade-off parameter of ||x||_0
%       delta               trade-off parameter of ||y||_0
%       A                   s x (p+q+r) LHS matrix
%       b                   s x 1 RHS vector
%       csense              s x 1 Constraint senses, a string containting the constraint sense for
%                           each row in A ('E', equality, 'G' greater than, 'L' less than).
%       lb                  (p+q+r) x 1 Lower bound vector
%       ub                  (p+q+r) x 1 Upper bound vector
%
% OPTIONAL INPUTS
% params                    parameters structure
%       nbMaxIteration      stopping criteria - number maximal of iteration (Defaut value = 1000)
%       epsilon             stopping criteria - (Defaut value = 10e-6)
%       theta               parameter of the approximation (Defaut value = 2)
%
% OUTPUT
% solution                  Structure containing the following fields
%       x                   p x 1 solution vector
%       y                   q x 1 solution vector
%       z                   r x 1 solution vector
%       stat                status
%                           1 =  Solution found
%                           2 =  Unbounded
%                           0 =  Infeasible
%                           -1=  Invalid input

author: opencobra

repo: cobratoolbox


text:
%Vehicle Abstract vehicle class
%
% This abstract class models the kinematics of a mobile robot moving on
% a plane and with a pose in SE(2).  For given steering and velocity inputs it
% updates the true vehicle state and returns noise-corrupted odometry
% readings.
%
% Methods::
%   Vehicle          constructor
%   add_driver       attach a driver object to this vehicle
%   control          generate the control inputs for the vehicle
%   f                predict next state based on odometry
%   init             initialize vehicle state
%   run              run for multiple time steps
%   run2             run with control inputs
%   step             move one time step and return noisy odometry
%   update           update the vehicle state
%
% Plotting/display methods::
%   char             convert to string
%   display          display state/parameters in human readable form
%   plot             plot/animate vehicle on current figure
%   plot_xy          plot the true path of the vehicle
%   Vehicle.plotv    plot/animate a pose on current figure
%
% Properties (read/write)::
%   x               true vehicle state: x, y, theta (3x1)
%   V               odometry covariance (2x2)
%   odometry        distance moved in the last interval (2x1)
%   rdim             dimension of the robot (for drawing)
%   L               length of the vehicle (wheelbase)
%   alphalim        steering wheel limit
%   speedmax        maximum vehicle speed
%   T               sample interval
%   verbose         verbosity
%   x_hist          history of true vehicle state (Nx3)
%   driver          reference to the driver object
%   x0              initial state, restored on init()
%
% Examples::
%
% If veh is an instance of a Vehicle class then we can add a driver object
%      veh.add_driver( RandomPath(10) )
% which will move the vehicle within the region -10<x<10, -10<y<10 which we
% can see by
%      veh.run(1000)
% which shows an animation of the vehicle moving for 1000 time steps
% between randomly selected wayoints.
%
% Notes::
% - Subclass of the MATLAB handle class which means that pass by reference semantics
%   apply.
%
% Reference::
%
%   Robotics, Vision & Control, Chap 6
%   Peter Corke,
%   Springer 2011
%
% See also Bicycle, Unicycle, RandomPath, EKF.

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

classdef Vehicle < handle

    properties
        % state
        x           % true state (x,y,theta)
        x_hist      % x history

        % parameters

        speedmax    % maximum speed
        dim         % dimension of the world -dim -> +dim in x and y
        rdim    % dimension of the robot
        dt           % sample interval
        V           % odometry covariance
        odometry    % distance moved in last interval
        verbose
        driver      % driver object
        x0          % initial state
        options

        vhandle     % handle to vehicle graphics object
        vtrail      % vehicle trail
    end

    methods(Abstract)
        f
    end

    methods

        function veh = Vehicle(varargin)
        %Vehicle Vehicle object constructor
        %
        % V = Vehicle(OPTIONS) creates a Vehicle object that implements the
        % kinematic model of a wheeled vehicle.
        %
        % Options::
        % 'covar',C       specify odometry covariance (2x2) (default 0)
        % 'speedmax',S    Maximum speed (default 1m/s)
        % 'L',L           Wheel base (default 1m)
        % 'x0',x0         Initial state (default (0,0,0) )
        % 'dt',T          Time interval (default 0.1)
        % 'rdim',R        Robot size as fraction of plot window (default 0.2)
        % 'verbose'       Be verbose
        %
        % Notes::
        % - The covariance is used by a "hidden" random number generator within the class.
        % - Subclasses the MATLAB handle class which means that pass by reference semantics
        %   apply.

            % vehicle common
            opt.covar = [];
            opt.rdim = 0.2;
            opt.dt = 0.1;
            opt.x0 = zeros(3,1);
            opt.speedmax = 1;
            opt.vhandle = [];

            [opt,args] = tb_optparse(opt, varargin);

            veh.V = opt.covar;
            veh.rdim = opt.rdim;
            veh.dt = opt.dt;
            veh.x0 = opt.x0(:);
            assert(isvec(veh.x0, 3), 'Initial configuration must be a 3-vector');
            veh.speedmax = opt.speedmax;
            veh.options = args;  % unused options go back to the subclass
            veh.vhandle = opt.vhandle;
            veh.x_hist = [];
        end

        function init(veh, x0)
            %Vehicle.init Reset state
            %
            % V.init() sets the state V.x := V.x0, initializes the driver
            % object (if attached) and clears the history.
            %
            % V.init(X0) as above but the state is initialized to X0.

            % TODO: should this be called from run?

            if nargin > 1
                veh.x = x0(:);
            else
                veh.x = veh.x0;
            end
            veh.x_hist = [];

            if ~isempty(veh.driver)
                veh.driver.init();
            end

            veh.vhandle = [];
        end

        function yy = path(veh, t, u, y0)
            %Vehicle.path Compute path for constant inputs
            %
            % XF = V.path(TF, U) is the final state of the vehicle (3x1) from the initial
            % state (0,0,0) with the control inputs U (vehicle specific).  TF is  a scalar to
            % specify the total integration time.
            %
            % XP = V.path(TV, U) is the trajectory of the vehicle (Nx3) from the initial
            % state (0,0,0) with the control inputs U (vehicle specific).  T is a vector (N) of
            % times for which elements of the trajectory will be computed.
            %
            % XP = V.path(T, U, X0) as above but specify the initial state.
            %
            % Notes::
            % - Integration is performed using ODE45.
            % - The ODE being integrated is given by the deriv method of the vehicle object.
            %
            % See also ODE45.

                if length(t) == 1
                    tt = [0 t];
                else
                    tt = t;
                end

                if nargin < 4
                    y0 = [0 0 0];
                end
                out = ode45( @(t,y) veh.deriv(t, y, u), tt, y0);

                y = out.y';
                if nargout == 0
                    plot(y(:,1), y(:,2));
                    grid on
                    xlabel('X'); ylabel('Y')
                else
                    yy = y;
                    if length(t) == 1
                        % if scalar time given, just return final state
                        yy = yy(end,:);
                    end
                end
        end

        function add_driver(veh, driver)
            %Vehicle.add_driver Add a driver for the vehicle
            %
            % V.add_driver(D) connects a driver object D to the vehicle.  The driver
            % object has one public method:
            %        [speed, steer] = D.demand();
            % that returns a speed and steer angle.
            %
            % Notes::
            % - The Vehicle.step() method invokes the driver if one is attached.
            %
            % See also Vehicle.step, RandomPath.
            veh.driver = driver;
            driver.veh = veh;
        end

        function odo = update(veh, u)
            %Vehicle.update Update the vehicle state
            %
            % ODO = V.update(U) is the true odometry value for
            % motion with U=[speed,steer].
            %
            % Notes::
            % - Appends new state to state history property x_hist.
            % - Odometry is also saved as property odometry.

            xp = veh.x; % previous state
            veh.x(1) = veh.x(1) + u(1)*veh.dt*cos(veh.x(3));
            veh.x(2) = veh.x(2) + u(1)*veh.dt*sin(veh.x(3));
            veh.x(3) = veh.x(3) + u(1)*veh.dt/veh.L * u(2);
            odo = [colnorm(veh.x(1:2)-xp(1:2)) veh.x(3)-xp(3)];
            veh.odometry = odo;

            veh.x_hist = [veh.x_hist; veh.x'];   % maintain history
        end

        function odo = step(veh, varargin)
            %Vehicle.step Advance one timestep
            %
            % ODO = V.step(SPEED, STEER) updates the vehicle state for one timestep
            % of motion at specified SPEED and STEER angle, and returns noisy odometry.
            %
            % ODO = V.step() updates the vehicle state for one timestep of motion and
            % returns noisy odometry.  If a "driver" is attached then its DEMAND() method
            % is invoked to compute speed and steer angle.  If no driver is attached
            % then speed and steer angle are assumed to be zero.
            %
            % Notes::
            % - Noise covariance is the property V.
            %
            % See also Vehicle.control, Vehicle.update, Vehicle.add_driver.

            % get the control input to the vehicle from either passed demand or driver
            u = veh.control(varargin{:});

            % compute the true odometry and update the state
            odo = veh.update(u);

            % add noise to the odometry
            if ~isempty(veh.V)
                odo = veh.odometry + randn(1,2)*sqrtm(veh.V);
            end
        end

        function u = control(veh, speed, steer)
            %Vehicle.control Compute the control input to vehicle
            %
            % U = V.control(SPEED, STEER) is a control input (1x2) = [speed,steer]
            % based on provided controls SPEED,STEER to which speed and steering angle
            % limits have been applied.
            %
            % U = V.control() as above but demand originates with a "driver" object if
            % one is attached, the driver's DEMAND() method is invoked. If no driver is
            % attached then speed and steer angle are assumed to be zero.
            %
            % See also Vehicle.step, RandomPath.
            if nargin < 2
                % if no explicit demand, and a driver is attached, use
                % it to provide demand
                if ~isempty(veh.driver)
                    [speed, steer] = veh.driver.demand();
                else
                    % no demand, do something safe
                    speed = 0;
                    steer = 0;
                end
            end

            % clip the speed
            if isempty(veh.speedmax)
                u(1) = speed;
            else
                u(1) = min(veh.speedmax, max(-veh.speedmax, speed));
            end

            % clip the steering angle
            if isprop(veh, 'steermax') && ~isempty(veh.steermax)
                u(2) = max(-veh.steermax, min(veh.steermax, steer));
            else
                u(2) = steer;
            end
        end

        function p = run(veh, nsteps)
            %Vehicle.run Run the vehicle simulation
            %
            % V.run(N) runs the vehicle model for N timesteps and plots
            % the vehicle pose at each step.
            %
            % P = V.run(N) runs the vehicle simulation for N timesteps and
            % return the state history (Nx3) without plotting.  Each row
            % is (x,y,theta).
            %
            % See also Vehicle.step, Vehicle.run2.

            if nargin < 2
                nsteps = 1000;
            end
            if ~isempty(veh.driver)
                veh.driver.init()
            end
            %veh.clear();
            if ~isempty(veh.driver)
                veh.driver.plot();
            end

            veh.plot();
            for i=1:nsteps
                veh.step();
                if nargout == 0
                    % if no output arguments then plot each step
                    veh.plot();
                    drawnow
                end
            end
            p = veh.x_hist;
        end

        % TODO run and run2 should become superclass methods...

        function p = run2(veh, T, x0, speed, steer)
            %Vehicle.run2 Run the vehicle simulation with control inputs
            %
            % P = V.run2(T, X0, SPEED, STEER) runs the vehicle model for a time T with
            % speed SPEED and steering angle STEER.  P (Nx3) is the path followed and
            % each row is (x,y,theta).
            %
            % Notes::
            % - Faster and more specific version of run() method.
            % - Used by the RRT planner.
            %
            % See also Vehicle.run, Vehicle.step, RRT.
            veh.init(x0);

            for i=1:(T/veh.dt)
                veh.update([speed steer]);
            end
            p = veh.x_hist;
        end

        function h = plot(veh, varargin)
        %Vehicle.plot Plot vehicle
        %
        % The vehicle is depicted graphically as a narrow triangle that travels
        % "point first" and has a length V.rdim.
        %
        % V.plot(OPTIONS) plots the vehicle on the current axes at a pose given by
        % the current robot state.  If the vehicle has been previously plotted its
        % pose is updated.
        %
        % V.plot(X, OPTIONS) as above but the robot pose is given by X (1x3).
        %
        % H = V.plotv(X, OPTIONS) draws a representation of a ground robot as an
        % oriented triangle with pose X (1x3) [x,y,theta].  H is a graphics handle.
        %
        % V.plotv(H, X) as above but updates the pose of the graphic represented
        % by the handle H to pose X.
        %
        % Options::
        % 'scale',S    Draw vehicle with length S x maximum axis dimension
        % 'size',S     Draw vehicle with length S
        % 'color',C    Color of vehicle.
        % 'fill'       Filled
        % 'trail',S    Trail with line style S, use line() name-value pairs
        %
        % Example::
        %          veh.plot('trail', {'Color', 'r', 'Marker', 'o', 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r', 'MarkerSize', 3})

        % Notes::
        % - The last two calls are useful if animating multiple robots in the same
        %   figure.
        %
        % See also Vehicle.plotv, plot_vehicle.

            if isempty(veh.vhandle)
                veh.vhandle = Vehicle.plotv(veh.x, varargin{:});
            end

            if ~isempty(varargin) && isnumeric(varargin{1})
                % V.plot(X)
                pos = varargin{1}; % use passed value
            else
                % V.plot()
                pos = veh.x;    % use current state
            end

            % animate it
            Vehicle.plotv(veh.vhandle, pos);

            end

        function out = plot_xy(veh, varargin)
            %Vehicle.plot_xy Plots true path followed by vehicle
            %
            % V.plot_xy() plots the true xy-plane path followed by the vehicle.
            %
            % V.plot_xy(LS) as above but the line style arguments LS are passed
            % to plot.
            %
            % Notes::
            % - The path is extracted from the x_hist property.

            xyt = veh.x_hist;
            if nargout == 0
                plot(xyt(:,1), xyt(:,2), varargin{:});
            else
                out = xyt;
            end
        end

        function verbosity(veh, v)
        %Vehicle.verbosity Set verbosity
        %
        % V.verbosity(A) set verbosity to A.  A=0 means silent.
            veh.verbose = v;
        end

        function display(nav)
        %Vehicle.display Display vehicle parameters and state
        %
        % V.display() displays vehicle parameters and state in compact
        % human readable form.
        %
        % Notes::
        % - This method is invoked implicitly at the command line when the result
        %   of an expression is a Vehicle object and the command has no trailing
        %   semicolon.
        %
        % See also Vehicle.char.

            loose = strcmp( get(0, 'FormatSpacing'), 'loose');
            if loose
                disp(' ');
            end
            disp([inputname(1), ' = '])
            disp( char(nav) );
        end % display()

        function s = char(veh)
        %Vehicle.char Convert to string
        %
        % s = V.char() is a string showing vehicle parameters and state in
        % a compact human readable format.
        %
        % See also Vehicle.display.

            s = '  Superclass: Vehicle';
            s = char(s, sprintf(...
            '    max speed=%g, dT=%g, nhist=%d', ...
                veh.speedmax, veh.dt, ...
                numrows(veh.x_hist)));
            if ~isempty(veh.V)
                s = char(s, sprintf(...
                '    V=(%g, %g)', ...
                    veh.V(1,1), veh.V(2,2)));
            end
            s = char(s, sprintf('    configuration: x=%g, y=%g, theta=%g', veh.x));
            if ~isempty(veh.driver)
                s = char(s, '    driven by::');
                s = char(s, [['      '; '      '] char(veh.driver)]);
            end
        end

    end % method

    methods(Static)

        function h = plotv(varargin)
        %Vehicle.plotv Plot ground vehicle pose
        %
        % H = Vehicle.plotv(X, OPTIONS) draws a representation of a ground robot as an
        % oriented triangle with pose X (1x3) [x,y,theta].  H is a graphics handle.
        % If X (Nx3) is a matrix it is considered to represent a trajectory in which case
        % the vehicle graphic is animated.
        %
        % Vehicle.plotv(H, X) as above but updates the pose of the graphic represented
        % by the handle H to pose X.
        %
        % Options::
        % 'scale',S       Draw vehicle with length S x maximum axis dimension
        % 'size',S        Draw vehicle with length S
        % 'fillcolor',C   Color of vehicle.
        % 'fps',F         Frames per second in animation mode (default 10)
        %
        % Example::
        %
        % Generate some path 3xN
        %         p = PRM.plan(start, goal);
        % Set the axis dimensions to stop them rescaling for every point on the path
        %         axis([-5 5 -5 5]);
        %
        % Now invoke the static method
        %         Vehicle.plotv(p);
        %
        % Notes::
        % - This is a class method.
        %
        % See also Vehicle.plot.

        if isstruct(varargin{1})
            plot_vehicle(varargin{2}, 'handle', varargin{1});
        else
            h = plot_vehicle(varargin{1}, 'fillcolor', 'b', 'alpha', 0.5);
        end

        end
    end % static methods

end % classdef

author: petercorke

repo: robotics-toolbox-matlab


text:
function month_cal_gregorian ( y, m )

%*****************************************************************************80
%
%% MONTH_CAL_GREGORIAN prints a Gregorian month calendar.
%
%  Format:
%
%    GREGORIAN CALENDAR
%        APRIL 1997 AD
%
%    Su  M Tu  W Th  F Sa
%           1  2  3  4  5
%     6  7  8  9 10 11 12
%    13 14 15 16 17 18 19
%    20 21 22 23 24 25 26
%    27 28 29 30
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer  Y, M, the YM date.
%

%
%  Make local copies of the input.
%
  m2 = m;
  y2 = y;
%
%  Check the month and year.  After this call, month is
%  guaranteed to be between 1 and 12.
%
  [ y2, m2, ierror ] = ym_check_gregorian ( y2, m2 );

  if ( ierror ~= 0 )
    return
  end
%
%  Find the day of the week for Y M 1.
%
  d = 1;
  f = 0.0;

  w = ymdf_to_weekday_gregorian ( y2, m2, d, f );
%
%  Find the appropriate label for the first box in the calendar.
%
  iday = 2 - w;
%
%  Print out a heading.
%
  s1 = month_to_month_name_common ( m2 );
  s2 = y_to_s_gregorian ( y2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Gregorian Calendar\n' );
  fprintf ( 1, '%s %s\n', s1, s2 );
  fprintf ( 1, '\n' );
%
%  Get the days of the week.
%
  fprintf ( 1, ' ' );
  for w = 1 : 7
    lab = weekday_to_name_common2 ( w );
    fprintf ( 1, '%3s', lab );
  end
  fprintf ( 1, '\n' );
%
%  Print out a line of day numbers.
%
  while ( iday <= month_length_gregorian ( y2, m2 ) )

    fprintf ( 1, ' ' );

    for w = 1 : 7

      if ( iday < 1 )
        fprintf ( 1, '   ' );
      elseif ( month_length_gregorian ( y2, m2 ) < iday )
        fprintf ( 1, '   ' );
      else
        fprintf ( 1, ' %2d', iday );
      end

      iday = iday + 1;

    end

    fprintf ( 1, '\n' );

  end

  return
end


text:
function out = Rmdm(A, B)

for i = 1:size(A, 3)
    for j = 1:size(B, 3)
        res(i, j) = Rdis(A(:, :, i), B(:, :, j));
    end
    out(i) = find(res(i, :) == min(res(i, :)));
end

author: PatternRecognition

repo: OpenBMI


text:
function gray_3x3 = image_denoise_gray_3x3 ( gray )

%*****************************************************************************80
%
%% IMAGE_DENOISE_GRAY_3X3 filters out noise from a gray scale image.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, uint8 GRAY(:,:), the noisy grayscale data.
%
%    Output, uint8 GRAY_3X3(:,:), the grayscale data for the filtered image.
%
  [ m, n ] = size ( gray );
  gray3d = zeros ( 9, m, n );
%
%  For pixel (I,J), GRAY3D(*,I,J) contains the current pixel value
%  and the 8 neighboring values.
%
  gray3d(1,:,:) = circshift ( gray, [ -1, -1 ] );
  gray3d(2,:,:) = circshift ( gray, [ -1,  0 ] );
  gray3d(3,:,:) = circshift ( gray, [ -1,  1 ] );
  gray3d(4,:,:) = circshift ( gray, [  0, -1 ] );
  gray3d(5,:,:) = circshift ( gray, [  0,  0 ] );
  gray3d(6,:,:) = circshift ( gray, [  0,  1 ] );
  gray3d(7,:,:) = circshift ( gray, [  1, -1 ] );
  gray3d(8,:,:) = circshift ( gray, [  1,  0 ] );
  gray3d(9,:,:) = circshift ( gray, [  1,  1 ] );
%
%  By taking the median of the 9 values, we hope to drop any noisy pixels.
%
  gray3d = median ( gray3d );
%
%  The MEDIAN command returned a 3D array with a first dimension of 1.
%  For convenience, we suppress the first dimension.
%
  gray_3x3(:,:) = gray3d(1,:,:);
%
%  GRAY is a UINT8 array.  The MEDIAN operation returned a DOUBLE value.
%  So we have to convert the array back to UINT8.
%
  gray_3x3 = uint8 ( gray_3x3 );

  return
end


text:
function [dist, pos] = distancePointEdge(point, edge)
%DISTANCEPOINTEDGE Minimum distance between a point and an edge.
%
%   DIST = distancePointEdge(POINT, EDGE);
%   Return the euclidean distance between edge EDGE and point POINT.
%   EDGE has the form: [x1 y1 x2 y2], and POINT is [x y].
%
%   If EDGE is N-by-4 array, result is 1-by-4 array computed for each edge.
%   If POINT is a N-by-2 array, the result is a N-by-1 array.
%   If both POINT and EDGE are array, the result is computed for each
%   point-edge couple, and stored into a NP-by-NE array.
%
%   [DIST POS] = distancePointEdge(POINT, EDGE);
%   Also returns the position of closest point on the edge. POS is
%   comprised between 0 (first point) and 1 (last point).
%
%   Eaxmple
%     % Distance between a point and an edge
%     distancePointEdge([3 4], [0 0 10 0])
%     ans =
%         4
%
%     % Distance between several points and one edge
%     points = [10 15; 15 10; 30 10];
%     edge   = [10 10 20 10];
%     distancePointEdge(points, edge)
%     ans =
%         5
%         0
%        10
%
%     % Distance between a point a several edges
%     point = [14 33];
%     edges  = [10 30 20 30; 20 30 20 40;20 40 10 40;10 40 10 30];
%     distancePointEdge(point, edges)
%     ans =
%         3    6    7    4
%
%
%   See also
%   edges2d, points2d, distancePoints, distancePointLine
%

% ------
% Author: David Legland
% E-mail: david.legland@nantes.inra.fr
% Created: 2004-04-07
% Copyright 2004-2022 INRA - BIA-BIBS

% direction vector of each edge (row vectors)
vx = (edge(:, 3) - edge(:,1))';
vy = (edge(:, 4) - edge(:,2))';

% squared length of edges, with a check of validity
delta = vx .* vx + vy .* vy;
invalidEdges = delta < eps;
delta(invalidEdges) = 1;

% difference of coordinates between point and edge first vertex
% (NP-by-NE arrays)
dx  = bsxfun(@minus, point(:, 1), edge(:, 1)');
dy  = bsxfun(@minus, point(:, 2), edge(:, 2)');

% compute position of points projected on the supporting line, by using
% normalised dot product (NP-by-NE array)
pos = bsxfun(@rdivide, bsxfun(@times, dx, vx) + bsxfun(@times, dy, vy), delta);

% ensure degenerated edges are correclty processed (consider the first
% vertex is the closest)
pos(:, invalidEdges) = 0;

% change position to ensure projected point is located on the edge
pos(pos < 0) = 0;
pos(pos > 1) = 1;

% compute distance between point and its projection on the edge
dist = hypot(bsxfun(@times, pos, vx) - dx, bsxfun(@times, pos, vy) - dy);

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% Copyright (C) 2016, by Arturo Gil Aparicio
%
% 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 direct_kinematics_symbolic_KUKA_LBR
% link lengths

syms q1 q2 q3 q4 q5 q6 Q7
robot = load_robot('KUKA', 'LBR_IIWA_R820_7DOF')

d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

% matrices DH
A01 = dh_sym(q1, d(1), a(1), alpha(1));
A12 = dh_sym(q2, d(2), a(2), alpha(2));
A23 = dh_sym(q3, d(3), a(3), alpha(3));
A34 = dh_sym(q4, d(4), a(4), alpha(4));
A45 = dh_sym(q5, d(5), a(5), alpha(5));
A56 = dh_sym(q6, d(6), a(6), alpha(6));
A67 = dh_sym(q7, d(7), a(7), alpha(7));

A02 = A01*A12;
A03 = A02*A23;
A04 = A03*A34;
A05 = A04*A45;
A06 = A05*A56;
A07 = A06*A67;

A07 = simplify(A07)

A03 = simplify(A03)

function A = dh_sym(theta, d, a, alpha)
syms q1 q2 q3 q4 q5 q6
% avoid almost zero elements in cos(alpha) and sin(alpha)
ca = cos(alpha);
sa = sin(alpha);
if abs(ca) < 1e-6
    ca = 0;
end
if abs(sa) < 1e-6
    sa = 0;
end

A=[cos(theta)  -ca*sin(theta)   sa*sin(theta)   a*cos(theta);
   sin(theta)   ca*cos(theta)  -sa*cos(theta)   a*sin(theta);
            0         sa             ca             d;
            0         0               0             1];

repo: ARTE


text:
% [INPUT]
% data = A float t-by-2 matrix (-Inf,Inf) representing the model input.
% e = A float (0,2] representing the exponent of the euclidean distance used to calculate the Distance Correlation (optional, default=1).
%
% [OUTPUT]
% dcor = A float [0,1] representing the Distance Correlation.
% rmss = A float [0,1] representing the RMS Similarity.

function [dcor,rmss] = similarity_statistics(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('data',@(x)validateattributes(x,{'double'},{'real' 'finite' '2d' 'nonempty' 'size' [NaN 2]}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    data = validate_input(ipr.data);

    nargoutchk(2,2);

    [dcor,rmss] = similarity_statistics_internal(data);

end

function [dcor,rmss] = similarity_statistics_internal(data)

    x1 = data(:,1);
    x2 = data(:,2);

    dcor = calculate_dcor(x1,x2);
    rmss = calculate_rmss(x1,x2);

end

function dcor = calculate_dcor(x1,x2)

    d1 = sqrt(bsxfun(@minus,x1,x1.').^2);
    m1 = mean(d1,1);
    k1 = bsxfun(@minus,bsxfun(@minus,d1,m1.'),m1) + mean(mean(d1));

    d2 = sqrt(bsxfun(@minus,x2,x2.').^2);
    m2 = mean(d2,1);
    k2 = bsxfun(@minus,bsxfun(@minus,d2,m2.'),m2) + mean(mean(d2));

    v1 = sqrt(mean(mean(k1 .* k1)));
    v2 = sqrt(mean(mean(k2 .* k2)));
    v = sqrt(v1 * v2);

    dcov = sqrt(mean(mean(k1 .* k2)));

    if (v > 0)
        dcor = dcov / v;
    else
        dcor = 0;
    end

end

function rmss = calculate_rmss(x1,x2)

    s = 1 - (abs(x1 - x2) ./ (abs(x1) + abs(x2)));
    rmss = sqrt(mean(s.^2,'omitnan'));

end

function data = validate_input(data)

    t = size(data,1);

    if (t < 5)
        error('The value of ''data'' is invalid. Expected input to be a matrix with at least 5 rows.');
    end

end

author: TommasoBelluzzo

repo: SystemicRisk