text:
function [der,errest,finaldelta] = derivest(fun,x0,varargin)
% DERIVEST: estimate the n'th derivative of fun at x0, provide an error estimate
% usage: [der,errest] = DERIVEST(fun,x0)  % first derivative
% usage: [der,errest] = DERIVEST(fun,x0,prop1,val1,prop2,val2,...)
%
% Derivest will perform numerical differentiation of an
% analytical function provided in fun. It will not
% differentiate a function provided as data. Use gradient
% for that purpose, or differentiate a spline model.
%
% The methods used by DERIVEST are finite difference
% approximations of various orders, coupled with a generalized
% (multiple term) Romberg extrapolation. This also yields
% the error estimate provided. DERIVEST uses a semi-adaptive
% scheme to provide the best estimate that it can by its
% automatic choice of a differencing interval.
%
% Finally, While I have not written this function for the
% absolute maximum speed, speed was a major consideration
% in the algorithmic design. Maximum accuracy was my main goal.
%
%
% Arguments (input)
%  fun - function to differentiate. May be an inline function,
%        anonymous, or an m-file. fun will be sampled at a set
%        of distinct points for each element of x0. If there are
%        additional parameters to be passed into fun, then use of
%        an anonymous function is recommended.
%
%        fun should be vectorized to allow evaluation at multiple
%        locations at once. This will provide the best possible
%        speed. IF fun is not so vectorized, then you MUST set
%        'vectorized' property to 'no', so that derivest will
%        then call your function sequentially instead.
%
%        Fun is assumed to return a result of the same
%        shape as its input x0.
%
%  x0  - scalar, vector, or array of points at which to
%        differentiate fun.
%
% Additional inputs must be in the form of property/value pairs.
%  Properties are character strings. They may be shortened
%  to the extent that they are unambiguous. Properties are
%  not case sensitive. Valid property names are:
%
%  'DerivativeOrder', 'MethodOrder', 'Style', 'RombergTerms'
%  'FixedStep', 'MaxStep'
%
%  All properties have default values, chosen as intelligently
%  as I could manage. Values that are character strings may
%  also be unambiguously shortened. The legal values for each
%  property are:
%
%  'DerivativeOrder' - specifies the derivative order estimated.
%        Must be a positive integer from the set [1,2,3,4].
%
%        DEFAULT: 1 (first derivative of fun)
%
%  'MethodOrder' - specifies the order of the basic method
%        used for the estimation.
%
%        For 'central' methods, must be a positive integer
%        from the set [2,4].
%
%        For 'forward' or 'backward' difference methods,
%        must be a positive integer from the set [1,2,3,4].
%
%        DEFAULT: 4 (a second order method)
%
%        Note: higher order methods will generally be more
%        accurate, but may also suffere more from numerical
%        problems.
%
%        Note: First order methods would usually not be
%        recommended.
%
%  'Style' - specifies the style of the basic method
%        used for the estimation. 'central', 'forward',
%        or 'backwards' difference methods are used.
%
%        Must be one of 'Central', 'forward', 'backward'.
%
%        DEFAULT: 'Central'
%
%        Note: Central difference methods are usually the
%        most accurate, but sometiems one must not allow
%        evaluation in one direction or the other.
%
%  'RombergTerms' - Allows the user to specify the generalized
%        Romberg extrapolation method used, or turn it off
%        completely.
%
%        Must be a positive integer from the set [0,1,2,3].
%
%        DEFAULT: 2 (Two Romberg terms)
%
%        Note: 0 disables the Romberg step completely.
%
%  'FixedStep' - Allows the specification of a fixed step
%        size, preventing the adaptive logic from working.
%        This will be considerably faster, but not necessarily
%        as accurate as allowing the adaptive logic to run.
%
%        DEFAULT: []
%
%        Note: If specified, 'FixedStep' will define the
%        maximum excursion from x0 that will be used.
%
%  'Vectorized' - Derivest will normally assume that your
%        function can be safely evaluated at multiple locations
%        in a single call. This would minimize the overhead of
%        a loop and additional function call overhead. Some
%        functions are not easily vectorizable, but you may
%        (if your matlab release is new enough) be able to use
%        arrayfun to accomplish the vectorization.
%
%        When all else fails, set the 'vectorized' property
%        to 'no'. This will cause derivest to loop over the
%        successive function calls.
%
%        DEFAULT: 'yes'
%
%
%  'MaxStep' - Specifies the maximum excursion from x0 that
%        will be allowed, as a multiple of x0.
%
%        DEFAULT: 100
%
%  'StepRatio' - Derivest uses a proportionally cascaded
%        series of function evaluations, moving away from your
%        point of evaluation. The StepRatio is the ratio used
%        between sequential steps.
%
%        DEFAULT: 2.0000001
%
%        Note: use of a non-integer stepratio is intentional,
%        to avoid integer multiples of the period of a periodic
%        function under some circumstances.
%
%
% See the document DERIVEST.pdf for more explanation of the
% algorithms behind the parameters of DERIVEST. In most cases,
% I have chosen good values for these parameters, so the user
% should never need to specify anything other than possibly
% the DerivativeOrder. I've also tried to make my code robust
% enough that it will not need much. But complete flexibility
% is in there for your use.
%
%
% Arguments: (output)
%  der - derivative estimate for each element of x0
%        der will have the same shape as x0.
%
%  errest - 95% uncertainty estimate of the derivative, such that
%
%        abs(der(j) - f'(x0(j))) < erest(j)
%
%  finaldelta - The final overall stepsize chosen by DERIVEST
%
%
% Example usage:
%  First derivative of exp(x), at x == 1
%   [d,e]=derivest(@(x) exp(x),1)
%   d =
%       2.71828182845904
%
%   e =
%       1.02015503167879e-14
%
%  True derivative
%   exp(1)
%   ans =
%       2.71828182845905
%
% Example usage:
%  Third derivative of x.^3+x.^4, at x = [0,1]
%   derivest(@(x) x.^3 + x.^4,[0 1],'deriv',3)
%   ans =
%       6       30
%
%  True derivatives: [6,30]
%
%
% See also: gradient
%
%
% Author: John D'Errico
% e-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 12/27/2006

par.DerivativeOrder = 1;
par.MethodOrder = 4;
par.Style = 'central';
par.RombergTerms = 2;
par.FixedStep = [];
par.MaxStep = 100;
% setting a default stepratio as a non-integer prevents
% integer multiples of the initial point from being used.
% In turn that avoids some problems for periodic functions.
par.StepRatio = 2.0000001;
par.NominalStep = [];
par.Vectorized = 'yes';

na = length(varargin);
if (rem(na,2)==1)
  error 'Property/value pairs must come as PAIRS of arguments.'
elseif na>0
  par = parse_pv_pairs(par,varargin);
end
par = check_params(par);

% Was fun a string, or an inline/anonymous function?
if (nargin<1)
  help derivest
  return
elseif isempty(fun)
  error 'fun was not supplied.'
elseif ischar(fun)
  % a character function name
  fun = str2func(fun);
end

% no default for x0
if (nargin<2) || isempty(x0)
  error 'x0 was not supplied'
end
par.NominalStep = max(x0,0.02);

% was a single point supplied?
nx0 = size(x0);
n = prod(nx0);

% Set the steps to use.
if isempty(par.FixedStep)
  % Basic sequence of steps, relative to a stepsize of 1.
  delta = par.MaxStep*par.StepRatio .^(0:-1:-25)';
  ndel = length(delta);
else
  % Fixed, user supplied absolute sequence of steps.
  ndel = 3 + ceil(par.DerivativeOrder/2) + ...
     par.MethodOrder + par.RombergTerms;
  if par.Style(1) == 'c'
    ndel = ndel - 2;
  end
  delta = par.FixedStep*par.StepRatio .^(-(0:(ndel-1)))';
end

% generate finite differencing rule in advance.
% The rule is for a nominal unit step size, and will
% be scaled later to reflect the local step size.
fdarule = 1;
switch par.Style
  case 'central'
    % for central rules, we will reduce the load by an
    % even or odd transformation as appropriate.
    if par.MethodOrder==2
      switch par.DerivativeOrder
        case 1
          % the odd transformation did all the work
          fdarule = 1;
        case 2
          % the even transformation did all the work
          fdarule = 2;
        case 3
          % the odd transformation did most of the work, but
          % we need to kill off the linear term
          fdarule = [0 1]/fdamat(par.StepRatio,1,2);
        case 4
          % the even transformation did most of the work, but
          % we need to kill off the quadratic term
          fdarule = [0 1]/fdamat(par.StepRatio,2,2);
      end
    else
      % a 4th order method. We've already ruled out the 1st
      % order methods since these are central rules.
      switch par.DerivativeOrder
        case 1
          % the odd transformation did most of the work, but
          % we need to kill off the cubic term
          fdarule = [1 0]/fdamat(par.StepRatio,1,2);
        case 2
          % the even transformation did most of the work, but
          % we need to kill off the quartic term
          fdarule = [1 0]/fdamat(par.StepRatio,2,2);
        case 3
          % the odd transformation did much of the work, but
          % we need to kill off the linear & quintic terms
          fdarule = [0 1 0]/fdamat(par.StepRatio,1,3);
        case 4
          % the even transformation did much of the work, but
          % we need to kill off the quadratic and 6th order terms
          fdarule = [0 1 0]/fdamat(par.StepRatio,2,3);
      end
    end
  case {'forward' 'backward'}
    % These two cases are identical, except at the very end,
    % where a sign will be introduced.

    % No odd/even trans, but we already dropped
    % off the constant term
    if par.MethodOrder==1
      if par.DerivativeOrder==1
        % an easy one
        fdarule = 1;
      else
        % 2:4
        v = zeros(1,par.DerivativeOrder);
        v(par.DerivativeOrder) = 1;
        fdarule = v/fdamat(par.StepRatio,0,par.DerivativeOrder);
      end
    else
      % par.MethodOrder methods drop off the lower order terms,
      % plus terms directly above DerivativeOrder
      v = zeros(1,par.DerivativeOrder + par.MethodOrder - 1);
      v(par.DerivativeOrder) = 1;
      fdarule = v/fdamat(par.StepRatio,0,par.DerivativeOrder+par.MethodOrder-1);
    end

    % correct sign for the 'backward' rule
    if par.Style(1) == 'b'
      fdarule = -fdarule;
    end

end % switch on par.style (generating fdarule)
nfda = length(fdarule);

% will we need fun(x0)?
if (rem(par.DerivativeOrder,2) == 0) || ~strncmpi(par.Style,'central',7)
  if strcmpi(par.Vectorized,'yes')
    f_x0 = fun(x0);
  else
    % not vectorized, so loop
    f_x0 = zeros(size(x0));
    for j = 1:numel(x0)
      f_x0(j) = fun(x0(j));
    end
  end
else
  f_x0 = [];
end

% Loop over the elements of x0, reducing it to
% a scalar problem. Sorry, vectorization is not
% complete here, but this IS only a single loop.
der = zeros(nx0);
errest = der;
finaldelta = der;
for i = 1:n
  x0i = x0(i);
  h = par.NominalStep(i);

  % a central, forward or backwards differencing rule?
  % f_del is the set of all the function evaluations we
  % will generate. For a central rule, it will have the
  % even or odd transformation built in.
  if par.Style(1) == 'c'
    % A central rule, so we will need to evaluate
    % symmetrically around x0i.
    if strcmpi(par.Vectorized,'yes')
      f_plusdel = fun(x0i+h*delta);
      f_minusdel = fun(x0i-h*delta);
    else
      % not vectorized, so loop
      f_minusdel = zeros(size(delta));
      f_plusdel = zeros(size(delta));
      for j = 1:numel(delta)
        f_plusdel(j) = fun(x0i+h*delta(j));
        f_minusdel(j) = fun(x0i-h*delta(j));
      end
    end

    if ismember(par.DerivativeOrder,[1 3])
      % odd transformation
      f_del = (f_plusdel - f_minusdel)/2;
    else
      f_del = (f_plusdel + f_minusdel)/2 - f_x0(i);
    end
  elseif par.Style(1) == 'f'
    % forward rule
    % drop off the constant only
    if strcmpi(par.Vectorized,'yes')
      f_del = fun(x0i+h*delta) - f_x0(i);
    else
      % not vectorized, so loop
      f_del = zeros(size(delta));
      for j = 1:numel(delta)
        f_del(j) = fun(x0i+h*delta(j)) - f_x0(i);
      end
    end
  else
    % backward rule
    % drop off the constant only
    if strcmpi(par.Vectorized,'yes')
      f_del = fun(x0i-h*delta) - f_x0(i);
    else
      % not vectorized, so loop
      f_del = zeros(size(delta));
      for j = 1:numel(delta)
        f_del(j) = fun(x0i-h*delta(j)) - f_x0(i);
      end
    end
  end

  % check the size of f_del to ensure it was properly vectorized.
  f_del = f_del(:);
  if length(f_del)~=ndel
    error 'fun did not return the correct size result (fun must be vectorized)'
  end

  % Apply the finite difference rule at each delta, scaling
  % as appropriate for delta and the requested DerivativeOrder.
  % First, decide how many of these estimates we will end up with.
  ne = ndel + 1 - nfda - par.RombergTerms;

  % Form the initial derivative estimates from the chosen
  % finite difference method.
  der_init = vec2mat(f_del,ne,nfda)*fdarule.';

  % scale to reflect the local delta
  der_init = der_init(:)./(h*delta(1:ne)).^par.DerivativeOrder;

  % Each approximation that results is an approximation
  % of order par.DerivativeOrder to the desired derivative.
  % Additional (higher order, even or odd) terms in the
  % Taylor series also remain. Use a generalized (multi-term)
  % Romberg extrapolation to improve these estimates.
  switch par.Style
    case 'central'
      rombexpon = 2*(1:par.RombergTerms) + par.MethodOrder - 2;
    otherwise
      rombexpon = (1:par.RombergTerms) + par.MethodOrder - 1;
  end
  [der_romb,errors] = rombextrap(par.StepRatio,der_init,rombexpon);

  % Choose which result to return

  % first, trim off the
  if isempty(par.FixedStep)
    % trim off the estimates at each end of the scale
    nest = length(der_romb);
    switch par.DerivativeOrder
      case {1 2}
        trim = [1 2 nest-1 nest];
      case 3
        trim = [1:4 nest+(-3:0)];
      case 4
        trim = [1:6 nest+(-5:0)];
    end

    [der_romb,tags] = sort(der_romb);

    der_romb(trim) = [];
    tags(trim) = [];
    errors = errors(tags);
    trimdelta = delta(tags);

    [errest(i),ind] = min(errors);

    finaldelta(i) = h*trimdelta(ind);
    der(i) = der_romb(ind);
  else
    [errest(i),ind] = min(errors);
    finaldelta(i) = h*delta(ind);
    der(i) = der_romb(ind);
  end
end

end % mainline end

% ============================================
% subfunction - romberg extrapolation
% ============================================
function [der_romb,errest] = rombextrap(StepRatio,der_init,rombexpon)
% do romberg extrapolation for each estimate
%
%  StepRatio - Ratio decrease in step
%  der_init - initial derivative estimates
%  rombexpon - higher order terms to cancel using the romberg step
%
%  der_romb - derivative estimates returned
%  errest - error estimates
%  amp - noise amplification factor due to the romberg step

srinv = 1/StepRatio;

% do nothing if no romberg terms
nexpon = length(rombexpon);
rmat = ones(nexpon+2,nexpon+1);
switch nexpon
  case 0
    % rmat is simple: ones(2,1)
  case 1
    % only one romberg term
    rmat(2,2) = srinv^rombexpon;
    rmat(3,2) = srinv^(2*rombexpon);
  case 2
    % two romberg terms
    rmat(2,2:3) = srinv.^rombexpon;
    rmat(3,2:3) = srinv.^(2*rombexpon);
    rmat(4,2:3) = srinv.^(3*rombexpon);
  case 3
    % three romberg terms
    rmat(2,2:4) = srinv.^rombexpon;
    rmat(3,2:4) = srinv.^(2*rombexpon);
    rmat(4,2:4) = srinv.^(3*rombexpon);
    rmat(5,2:4) = srinv.^(4*rombexpon);
end

% qr factorization used for the extrapolation as well
% as the uncertainty estimates
[qromb,rromb] = qr(rmat,0);

% the noise amplification is further amplified by the Romberg step.
% amp = cond(rromb);

% this does the extrapolation to a zero step size.
ne = length(der_init);
rhs = vec2mat(der_init,nexpon+2,max(1,ne - (nexpon+2)));
rombcoefs = rromb\(qromb.'*rhs);
der_romb = rombcoefs(1,:).';

% uncertainty estimate of derivative prediction
s = sqrt(sum((rhs - rmat*rombcoefs).^2,1));
rinv = rromb\eye(nexpon+1);
cov1 = sum(rinv.^2,2); % 1 spare dof
errest = s.'*12.7062047361747*sqrt(cov1(1));

end % rombextrap

% ============================================
% subfunction - vec2mat
% ============================================
function mat = vec2mat(vec,n,m)
% forms the matrix M, such that M(i,j) = vec(i+j-1)
[i,j] = ndgrid(1:n,0:m-1);
ind = i+j;
mat = vec(ind);
if n==1
  mat = mat.';
end

end % vec2mat

% ============================================
% subfunction - fdamat
% ============================================
function mat = fdamat(sr,parity,nterms)
% Compute matrix for fda derivation.
% parity can be
%   0 (one sided, all terms included but zeroth order)
%   1 (only odd terms included)
%   2 (only even terms included)
% nterms - number of terms

% sr is the ratio between successive steps
srinv = 1./sr;

switch parity
  case 0
    % single sided rule
    [i,j] = ndgrid(1:nterms);
    c = 1./factorial(1:nterms);
    mat = c(j).*srinv.^((i-1).*j);
  case 1
    % odd order derivative
    [i,j] = ndgrid(1:nterms);
    c = 1./factorial(1:2:(2*nterms));
    mat = c(j).*srinv.^((i-1).*(2*j-1));
  case 2
    % even order derivative
    [i,j] = ndgrid(1:nterms);
    c = 1./factorial(2:2:(2*nterms));
    mat = c(j).*srinv.^((i-1).*(2*j));
end

end % fdamat

% ============================================
% subfunction - check_params
% ============================================
function par = check_params(par)
% check the parameters for acceptability
%
% Defaults
% par.DerivativeOrder = 1;
% par.MethodOrder = 2;
% par.Style = 'central';
% par.RombergTerms = 2;
% par.FixedStep = [];

% DerivativeOrder == 1 by default
if isempty(par.DerivativeOrder)
  par.DerivativeOrder = 1;
else
  if (length(par.DerivativeOrder)>1) || ~ismember(par.DerivativeOrder,1:4)
    error 'DerivativeOrder must be scalar, one of [1 2 3 4].'
  end
end

% MethodOrder == 2 by default
if isempty(par.MethodOrder)
  par.MethodOrder = 2;
else
  if (length(par.MethodOrder)>1) || ~ismember(par.MethodOrder,[1 2 3 4])
    error 'MethodOrder must be scalar, one of [1 2 3 4].'
  elseif ismember(par.MethodOrder,[1 3]) && (par.Style(1)=='c')
    error 'MethodOrder==1 or 3 is not possible with central difference methods'
  end
end

% style is char
valid = {'central', 'forward', 'backward'};
if isempty(par.Style)
  par.Style = 'central';
elseif ~ischar(par.Style)
  error 'Invalid Style: Must be character'
end
ind = find(strncmpi(par.Style,valid,length(par.Style)));
if (length(ind)==1)
  par.Style = valid{ind};
else
  error(['Invalid Style: ',par.Style])
end

% vectorized is char
valid = {'yes', 'no'};
if isempty(par.Vectorized)
  par.Vectorized = 'yes';
elseif ~ischar(par.Vectorized)
  error 'Invalid Vectorized: Must be character'
end
ind = find(strncmpi(par.Vectorized,valid,length(par.Vectorized)));
if (length(ind)==1)
  par.Vectorized = valid{ind};
else
  error(['Invalid Vectorized: ',par.Vectorized])
end

% RombergTerms == 2 by default
if isempty(par.RombergTerms)
  par.RombergTerms = 2;
else
  if (length(par.RombergTerms)>1) || ~ismember(par.RombergTerms,0:3)
    error 'RombergTerms must be scalar, one of [0 1 2 3].'
  end
end

% FixedStep == [] by default
if (length(par.FixedStep)>1) || (~isempty(par.FixedStep) && (par.FixedStep<=0))
  error 'FixedStep must be empty or a scalar, >0.'
end

% MaxStep == 10 by default
if isempty(par.MaxStep)
  par.MaxStep = 10;
elseif (length(par.MaxStep)>1) || (par.MaxStep<=0)
  error 'MaxStep must be empty or a scalar, >0.'
end

end % check_params

% ============================================
% Included subfunction - parse_pv_pairs
% ============================================
function params=parse_pv_pairs(params,pv_pairs)
% parse_pv_pairs: parses sets of property value pairs, allows defaults
% usage: params=parse_pv_pairs(default_params,pv_pairs)
%
% arguments: (input)
%  default_params - structure, with one field for every potential
%             property/value pair. Each field will contain the default
%             value for that property. If no default is supplied for a
%             given property, then that field must be empty.
%
%  pv_array - cell array of property/value pairs.
%             Case is ignored when comparing properties to the list
%             of field names. Also, any unambiguous shortening of a
%             field/property name is allowed.
%
% arguments: (output)
%  params   - parameter struct that reflects any updated property/value
%             pairs in the pv_array.
%
% Example usage:
% First, set default values for the parameters. Assume we
% have four parameters that we wish to use optionally in
% the function examplefun.
%
%  - 'viscosity', which will have a default value of 1
%  - 'volume', which will default to 1
%  - 'pie' - which will have default value 3.141592653589793
%  - 'description' - a text field, left empty by default
%
% The first argument to examplefun is one which will always be
% supplied.
%
%   function examplefun(dummyarg1,varargin)
%   params.Viscosity = 1;
%   params.Volume = 1;
%   params.Pie = 3.141592653589793
%
%   params.Description = '';
%   params=parse_pv_pairs(params,varargin);
%   params
%
% Use examplefun, overriding the defaults for 'pie', 'viscosity'
% and 'description'. The 'volume' parameter is left at its default.
%
%   examplefun(rand(10),'vis',10,'pie',3,'Description','Hello world')
%
% params =
%     Viscosity: 10
%        Volume: 1
%           Pie: 3
%   Description: 'Hello world'
%
% Note that capitalization was ignored, and the property 'viscosity'
% was truncated as supplied. Also note that the order the pairs were
% supplied was arbitrary.

npv = length(pv_pairs);
n = npv/2;

if n~=floor(n)
  error 'Property/value pairs must come in PAIRS.'
end
if n<=0
  % just return the defaults
  return
end

if ~isstruct(params)
  error 'No structure for defaults was supplied'
end

% there was at least one pv pair. process any supplied
propnames = fieldnames(params);
lpropnames = lower(propnames);
for i=1:n
  p_i = lower(pv_pairs{2*i-1});
  v_i = pv_pairs{2*i};

  ind = strmatch(p_i,lpropnames,'exact');
  if isempty(ind)
    ind = find(strncmp(p_i,lpropnames,length(p_i)));
    if isempty(ind)
      error(['No matching property found for: ',pv_pairs{2*i-1}])
    elseif length(ind)>1
      error(['Ambiguous property name: ',pv_pairs{2*i-1}])
    end
  end
  p_i = propnames{ind};

  % override the corresponding default in params
  params = setfield(params,p_i,v_i); %#ok

end

end % parse_pv_pairs

author: yorgoon

repo: minimum-snap-geometric-control


text:
function callback_slider_prob(sliderHandle,~,sRate,epoched,epoch_length,ax1,latency,events2plot,events,axLikelihood,windowInS,ylim1,ylim2)

figh = gcf;
currentFrame = sliderHandle.Value;
data = get(figh,'userdata');
windowlength = data{2};
data(2) = [];
data = data{1};
if gca ~= ax1
    axes(ax1);
end

%children = get(ax,'children');
if ~epoched

    T = ceil(windowlength);
    if sRate*(currentFrame + T)+1 > size(data,2)
        plot(ax1,data(end,(sRate*currentFrame+1:end)),data(1:end-2,(sRate*currentFrame+1:end)))
        set(ax1,'xlim',[currentFrame ceil(size(data,2)/sRate)]);
        xlim = get(ax1,'xlim');
        set(ax1,'xtick',linspace(xlim(1),xlim(2),11))

    else
        plot(ax1,data(end,(sRate*currentFrame:sRate*(currentFrame + T))+1),data(1:end-2,(sRate*currentFrame:sRate*(currentFrame + T))+1))
        set(ax1,'xlim',[currentFrame currentFrame+T]);
        set(ax1,'xtick',linspace(currentFrame,currentFrame+T,11));
    end

    xlim = get(ax1,'xlim');
    set(ax1,'ylim',[-0.1 1.1]);

    ylabel(ax1,'Probability of Model Being Active')
    xlabel(ax1,fastif(epoched,'Trials','Time(sec)'))
    hold(ax1,'on');
    %draw events
    x_coord = [];
    if ~isempty(events2plot)
        count = 1;

        for i = 1:length(latency)
            if latency(i)>sRate*xlim(2)
                break;
            else
                if latency(i)>=currentFrame*sRate && ismember({events(i).type},events2plot)
                    x_coord(count) = latency(i)/sRate;
                    plot(ax1,[x_coord(count) x_coord(count)],[-0.1 1.1],'linestyle','- -')
                    ht = text(x_coord(count),1.12,num2str(events(i).type));
                    set(ht,'rotation',90);
                    count = count + 1;
                end
            end
        end
    end
    hold(ax1,'off')

    if sRate*(currentFrame + T)+1 > size(data,2)
        plot(axLikelihood,data(end,(sRate*currentFrame+1:end)),data(end-1,(sRate*currentFrame+1:end)))
        set(axLikelihood,'xlim',[currentFrame ceil(size(data,2)/sRate)]);

        set(axLikelihood,'xtick',linspace(xlim(1),xlim(2),11))
    else
        plot(axLikelihood,data(end,(sRate*currentFrame:sRate*(currentFrame + T))+1),data(end-1,(sRate*currentFrame:sRate*(currentFrame + T))+1));
        set(axLikelihood,'xlim',[currentFrame currentFrame+T]);
        set(axLikelihood,'xtick',linspace(currentFrame,currentFrame+T,11));
    end

    set(axLikelihood,'ylim',[ylim1 ylim2]);

    xlabel(axLikelihood,fastif(epoched,'Trials','Time(sec)'))
    hold(axLikelihood,'on');

    %draw events
    for i = 1:length(x_coord)
        plot(axLikelihood,[x_coord(i) x_coord(i)],[ylim1 ylim2],'linestyle','- -')
    end

    hold(axLikelihood,'off')
    if T ~= windowlength
        set(findobj('parent',gcf,'tag','WindowSizeTag'),'string',num2str(T))
    end

    sliderHandle.UnitIncrement = ceil(T/10);
    sliderHandle.Maximum = size(data,2)/sRate-T + ceil(T/10);
    sliderHandle.VisibleAmount = sliderHandle.UnitIncrement;

else

    T = ceil(windowlength);

    if epoch_length*(currentFrame + T) > size(data,2)
        currentFrame = size(data,2)/epoch_length - T;
        sliderHandle.Value = currentFrame;
    end

    plot(ax1,data(end,(currentFrame*epoch_length+1:epoch_length*(currentFrame + T))),data(1:end-2,(currentFrame*epoch_length+1:epoch_length*(currentFrame + T))))
    hold(ax1,'on');
    set(ax1,'ylim',[-0.1 1.1]);
    set(ax1,'xlim',[currentFrame+0.5 currentFrame+T+0.5]);
    lim = get(ax1,'xlim');
    set(ax1,'xtick',floor(lim(1):lim(2)))
    ylabel(ax1,'Probability of Model Being Active')
    xlabel(ax1,fastif(epoched,'Trials','Time(sec)'))

    %draw epoch seperating lines
    for i = currentFrame+1:currentFrame+T
        plot(ax1,[i+0.5 i+0.5],[-0.1 1.1],'color',[1 0.5 0.25],'linestyle','- -');
    end
    %-------------------------------
    windowInMs = 1000*windowInS;
    x_coord = [];
    if ~isempty(events2plot)
    % draw event lines

    C = 1/(windowInMs(2)-windowInMs(1));
    count = 1;

    for i = currentFrame+1:currentFrame+T
        lengthevents = length(events(i).event);
        if lengthevents == 1
            try
                temp = events(i).eventtype;
            catch
                temp = events(i).eventtype{1};
            end
        else
            temp = '';
        end

        for j = 1:lengthevents

            if lengthevents == 1
                str = temp{1};
            else
                str = events(i).eventtype{j};
            end

            if ismember(num2str(str),events2plot)
            x_coord(count) = C*([events(i).eventlatency{j}]-windowInMs(1)) + (i - 0.5);
            plot(ax1,[x_coord(count) x_coord(count)],[-0.1 1.1],'linestyle','- -')
            ht = text(x_coord(count),1.12,num2str(events(i).eventtype{j}));
            set(ht,'rotation',90);
            count = count + 1;
            end
        end
    end

    end
    hold(ax1,'off');

    plot(axLikelihood,data(end,(currentFrame*epoch_length+1:epoch_length*(currentFrame + T))),data(end-1,(currentFrame*epoch_length+1:epoch_length*(currentFrame + T))))
    set(axLikelihood,'ylim',[ylim1 ylim2]);
    set(axLikelihood,'xlim',[currentFrame+0.5 currentFrame+T+0.5]);
    set(axLikelihood,'xtick',floor(lim(1):lim(2)))
    xlabel(axLikelihood,fastif(epoched,'Trials','Time(sec)'))
    hold(axLikelihood,'on');

    for i = currentFrame+1:currentFrame+T
        plot(axLikelihood,[i+0.5 i+0.5],[ylim1 ylim2],'color',[1 0.5 0.25],'linestyle','- -');
    end

    % draw event lines

    for i = 1:length(x_coord)
        plot(axLikelihood,[x_coord(i) x_coord(i)],[ylim1 ylim2],'linestyle','- -');
    end

    hold(axLikelihood,'off');
    if T ~= windowlength
        set(findobj('parent',gcf,'tag','WindowSizeTag'),'string',num2str(T))
    end

    sliderHandle.UnitIncrement = ceil(T/5);
    sliderHandle.Maximum = size(data,2)/epoch_length-T + ceil(T/5);
    sliderHandle.VisibleAmount = sliderHandle.UnitIncrement;

end
ylabel(axLikelihood,'Log-likelihood of data under most probable model')

end

author: PatternRecognition

repo: OpenBMI


text:
function Cube = InitCube
% initialize the cube
Cube = [];
SmCube = ...
  [[0 0 0]; [1 0 0]; [0 1 0]; [1 1 0]; ...
   [0 0 1]; [1 0 1]; [0 1 1]; [1 1 1]];
Cube.Nodes = [];
for x=1:3,
  for y=1:3,
    for z=1:3,
      Cube.Nodes(:,1,(z-1)*9+(y-1)*3+x) = SmCube(:,1)*2+2*x-5;
      Cube.Nodes(:,2,(z-1)*9+(y-1)*3+x) = SmCube(:,2)*2+2*y-5;
      Cube.Nodes(:,3,(z-1)*9+(y-1)*3+x) = SmCube(:,3)*2+2*z-5;
    end
  end
end
Cube.Order = [[1 5 7 3];[2 4 8 6];[1 2 6 5];...
              [3 7 8 4];[1 3 4 2];[5 6 8 7]];
NumColor = [[1:3:25];[3:3:27];[1 2 3 10 11 12 19 20 21];...
  [7 8 9 16 17 18 25 26 27];[1:9];[19:27]];
ValColor = [[1 0 0];[1 0 1];[0 0 1];[0 1 0];[1 1 0];[1 1 1]];
Cube.Color = zeros(6,3,27); % black
for x=1:6,
  for y=1:9,
    Cube.Color(x,:,NumColor(x,y)) = ValColor(x,:);
  end
end
Cube.RotCubeCube = ... % for rotate cube
  [[ 7  8  9 16 17 18 25 26 27  4  5  6 13 14 15 22 23 24  1  2  3 10 11 12 19 20 21];...
   [19 20 21 10 11 12  1  2  3 22 23 24 13 14 15  4  5  6 25 26 27 16 17 18  7  8  9];...
   [19 10  1 22 13  4 25 16  7 20 11  2 23 14  5 26 17  8 21 12  3 24 15  6 27 18  9];...
   [ 3 12 21  6 15 24  9 18 27  2 11 20  5 14 23  8 17 26  1 10 19  4 13 22  7 16 25];...
   [ 3  6  9  2  5  8  1  4  7 12 15 18 11 14 17 10 13 16 21 24 27 20 23 26 19 22 25];...
   [ 7  4  1  8  5  2  9  6  3 16 13 10 17 14 11 18 15 12 25 22 19 26 23 20 27 24 21]];
Cube.RotCubeFlat = ... % for rotate cube
  [[1 2 5 6 4 3];...
   [1 2 6 5 3 4];...
   [6 5 3 4 1 2];...
   [5 6 3 4 2 1];...
   [3 4 2 1 5 6];...
   [4 3 1 2 5 6]];
Cube.RotLayerCube = ... % for rotate layer
  [[1:3:25];[3:3:27];[1 2 3 10 11 12 19 20 21];[7 8 9 16 17 18 25 26 27];[1:9];[19:27]];
Cube.NumColor = NumColor;
% Mouse polygons for any rotations
Cube.PgRotCubeX = ... % for rotate the cube around OX, OY, OZ
  [[168;114;171;225;168],[292;292;346;346;292],...
   [232;292;346;287;232],[292;347;347;292;292],...
   [232;171;119;180;232],[119;180;180;119;119],...
   [340;278;225;287;340],[119;180;180;119;119],...
   [240;240;180;180;240],[400;400;346;346;400],...
   [60;60;119;119;60],[292;292;240;240;292]];
Cube.PgRotCubeY(:,:,1) = ...
  [[373;345;323;350;373],[273;214;239;299;273],...
   [296;273;299;324;296],[89;117;176;148;89],...
   [296;323;294;272;296],[107;85;146;170;107],...
   [352;375;350;324;352],[234;209;272;294;234],...
   [123;187;209;146;123],[203;265;239;176;203],...
   [256;191;170;234;256],[148;214;187;123;148]];
Cube.PgRotLayX = ... % for rotate the layer
  [[240;292;232;240],[240;292;292;240],...
   [60;171;114;60],[240;292;292;240],...
   [220;168;225;220],[400;346;346;400],...
   [400;287;346;400],[400;347;347;400],...
   [400;340;287;400],[240;180;180;240],...
   [240;232;180;240],[240;180;180;240],...
   [60;119;171;60],[60;119;119;60],...
   [220;225;278;220],[60;119;119;60],...
   [60;60;119;60],[240;240;292;240],...
   [240;180;240;240],[400;346;400;400],...
   [240;240;180;240],[400;400;347;400],...
   [60;119;60;60],[240;292;240;240]];
Cube.PgRotLayY = ...
  [[246;273;296;246],[60;89;148;60],...
   [319;323;345;319],[246;214;273;246],...
   [398;373;350;398],[329;299;239;329],...
   [329;324;299;329],[142;176;117;142],...
   [329;352;324;329],[246;272;209;246],...
   [246;296;272;246],[60;146;85;60],...
   [319;294;323;319],[133;107;170;133],...
   [398;350;375;398],[319;234;294;319],...
   [319;256;234;319],[246;187;214;246],...
   [246;209;187;246],[329;239;265;329],...
   [60;123;146;60],[142;203;176;142],...
   [133;170;191;133],[60;148;123;60]];
Cube.hand = ...
  [[NaN,NaN,NaN,NaN,NaN,NaN,NaN,  1,  1,NaN,NaN,NaN,NaN,NaN,NaN,NaN];...
   [NaN,NaN,NaN,NaN,  1,  1,  1,  2,  2,  1,  1,  1,NaN,NaN,NaN,NaN];...
   [NaN,NaN,NaN,  1,  2,  2,  1,  2,  2,  1,  2,  2,  1,NaN,NaN,NaN];...
   [NaN,NaN,NaN,  1,  2,  2,  1,  2,  2,  1,  2,  2,  1,NaN,NaN,NaN];...
   [NaN,NaN,NaN,  1,  2,  2,  1,  2,  2,  1,  2,  2,  1,  1,  1,NaN];...
   [NaN,  1,  1,  1,  2,  2,  1,  2,  2,  1,  2,  2,  1,  2,  2,  1];...
   [  1,  2,  2,  1,  2,  2,  1,  2,  2,  1,  2,  2,  1,  2,  2,  1];...
   [  1,  2,  2,  1,  2,  2,  1,  2,  2,  1,  2,  2,  1,  2,  2,  1];...
   [  1,  2,  2,  1,  2,  2,  1,  2,  2,  1,  2,  2,  1,  2,  2,  1];...
   [  1,  2,  2,  1,  2,  2,  1,  2,  2,  1,  2,  2,  1,  2,  2,  1];...
   [  1,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  1,  2,  2,  1];...
   [  1,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  1];...
   [NaN,  1,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  1];...
   [NaN,  1,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  1,NaN];...
   [NaN,NaN,  1,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  1,NaN];...
   [NaN,NaN,  1,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  1,NaN]];
return

author: Sable

repo: mcbench-benchmarks


text:
function [nmps2] = cmps22nmps2(cmps2)
% Convert acceleration from centimeers per second squared to nanometers
% per second-squared
% Chad A. Greene 2012
nmps2 = cmps2*1e+7;
end


text:
% EDGELINK - Link edge points in an image into lists
%
% Usage: [edgelist edgeim] = edgelink3(im, minlength, location)
%
% Arguments:  im         - Binary edge image, it is assumed that edges
%                          have been thinned.
%             minlength  - Minimum edge length of interest
%             location   - Optional complex valued image holding subpixel
%                          locations of edge points. For any pixel the
%                          real part holds the subpixel row coordinate of
%                          that edge point and the imaginary part holds
%                          the column coordinate.  See NONMAXSUP.  If
%                          this argument is supplied the edgelists will
%                          be formed from the subpixel coordinates,
%                          otherwise the the integer pixel coordinates of
%                          points in 'im' are used.
%
% Returns:  edgelist - a cell array of edge lists in row,column coords in
%                      the form
%                     { [r1 c1   [r1 c1   etc }
%                        r2 c2    ...
%                        ...
%                        rN cN]   ....]
%
%           edgeim   - Image with pixels labeled with edge number.
%
%
% This function links edge points together into chains.  Where an edge
% diverges at a junction the function simply tracks one of the branches.
% The other branch is eventually processed as another edge.  These `broken
% branches' can be remerged by MERGESEG after a call to LINESEG.
%
% See also:  DRAWEDGELIST, LINESEG, MAXLINEDEV, MERGESEG, DRAWSEG, NONMAXSUP

% Acknowledgement:
% This code is inspired by David Lowe's Link.c function from the Vista image
% processing library developed at the University of British Columbia
%    http://www.cs.ubc.ca/nest/lci/vista/vista.html

% Copyright (c) 2001-2005 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.

% February  2001 - Original version
% September 2004 - Revised to allow subpixel edge data to be used
% January   2006 - Edited for TF Peaks by Luke Rankine
% March     2006 - Modified to handle images with no components
%                - Corrected the nearest neighbour set for IF estimation

function [edgelist, edgeim] = edgelink3(im, minlength, location)

global EDGEIM;      % Some global variables to avoid passing (and
% copying) of arguments, this improves speed.
global ROWS;
global COLS;

elist = {};

EDGEIM = im ~= 0;                     % make sure image is binary.

%   EDGEIM = bwmorph(EDGEIM,'thin',Inf);  % make sure edges are thinned.
%   show(EDGEIM,1)

EDGEIM = double(EDGEIM);
[ROWS, COLS] = size(EDGEIM);
edgeNo = 1;

% Perform raster scan through image looking for edge points.  When a
% point is found trackedge is called to find the rest of the edge
% points.  As it finds the points the edge image pixels are labeled
% with the -ve of their edge No

edgelist = {}; % Modified, Luke. Initializing variable so it can be returned
for r = 1:ROWS
    for c = 1:COLS
        if EDGEIM(r,c) == 1
            edgepoints = trackedge(r,c, edgeNo, minlength);
            if ~isempty(edgepoints)
                edgelist{edgeNo} = edgepoints;
                edgeNo = edgeNo + 1;
            end
        end
    end
end

% Modified, Luke. Checking to see if any components were linked
if isempty(edgelist) == 1
    edgeim = zeros(ROWS,COLS);
    return
else

    edgeim = -EDGEIM;  % Finally negate image to make edge encodings +ve.

    % If subpixel edge locations are supplied upgrade the integer precision
    % edgelists that were constructed with data from 'location'.
    if nargin == 3
        for I = 1:length(edgelist)
            ind = sub2ind(size(im),edgelist{I}(:,1),edgelist{I}(:,2));
            edgelist{I}(:,1) = real(location(ind))';
            edgelist{I}(:,2) = imag(location(ind))';
        end
    end

    %   show(edgeim,2), colormap(flag);
end

%----------------------------------------------------------------------
% TRACKEDGE
%
% Function to track all the edge points associated with a start point.  From
% a given starting point it tracks in one direction, storing the coords of
% the edge points in an array and labelling the pixels in the edge image
% with the -ve of their edge number.  When no more connected points are
% found the function returns to the start point and tracks in the opposite
% direction.  Finally a check for the overall number of edge points found is
% made and the edge is ignored if it is too short.
%
% Note that when a junction is encountered along an edge the function
% simply tracks one of the branches.  The other branch is eventually
% processed as another edge.  These `broken branches' can be remerged by
% MERGESEG.
%
% Usage:   edgepoints = trackedge(rstart, cstart, edgeNo, minlength)
%
% Arguments:   rstart, cstart   - row and column No of starting point
%              edgeNo           - the current edge number
%              minlength        - minimum length of edge to accept
%
% Returns:     edgepoints       - Nx2 array of row and col values for
%                                 each edge point.
%                                 An empty array is returned if the edge
%                                 is less than minlength long.

function edgepoints = trackedge(rstart, cstart, edgeNo, minlength)

global EDGEIM;
global ROWS;
global COLS;

edgepoints = [rstart cstart];      % Start a new list for this edge.
EDGEIM(rstart,cstart) = -edgeNo;   % Edge points in the image are
% encoded by -ve of their edgeNo.

[thereIsAPoint, r, c] = nextpoint(rstart,cstart); % Find next connected
% edge point.

while thereIsAPoint
    edgepoints = [edgepoints             % Add point to point list
        r    c   ];
    EDGEIM(r,c) = -edgeNo;               % Update edge image
    [thereIsAPoint, r, c] = nextpoint(r,c);
end

edgepoints = flipud(edgepoints);  % Reverse order of points

% Now track from original point in the opposite direction

[thereIsAPoint, r, c] = nextpoint(rstart,cstart);

while thereIsAPoint
    edgepoints = [edgepoints
        r    c   ];
    EDGEIM(r,c) = -edgeNo;
    [thereIsAPoint, r, c] = nextpoint(r,c);
end

% Reject short edges
Npoints = size(edgepoints,1);
if Npoints < minlength
    for i = 1:Npoints   % Clear pixels in the edge image
        EDGEIM(edgepoints(i,1), edgepoints(i,2)) = 0;
    end
    edgepoints = [];    % Return empty array
end

%----------------------------------------------------------------------
%
% NEXTPOINT
%
% Function finds a point that is 8 connected to an existing edge point
%

function [thereIsAPoint, nextr, nextc] = nextpoint(rp,cp)

global EDGEIM;
global ROWS;
global COLS;

% Modified by Luke Rankine, Jan. 2006
% row and column offsets for the 10(originally eight) neighbours of a
% point starting with those that are semi 4-connected.
%     roff = [1  0 -1  0  1  1 -1 -1  2 -2 -2  2];
%     coff = [0  1  0 -1  1 -1 -1  1  1  1 -1 -1];
roff = [0  0 1  1 2  2 -1 -1 -2 -2];
coff = [1 -1 1 -1 1 -1  1 -1  1 -1];

r = rp+roff;
c = cp+coff;

% Search through neighbours and return first connected edge point
for i = 1:length(roff)
    if r(i) >= 1 & r(i) <= ROWS & c(i) >= 1 & c(i) <= COLS
        if EDGEIM(r(i),c(i)) == 1
            nextr = r(i);
            nextc = c(i);
            thereIsAPoint = 1;
            return;             % break out and return with the data
        end
    end
end

% If we get here there was no connecting next point.
nextr = 0;
nextc = 0;
thereIsAPoint = 0;

author: HeLiangHIT


text:
function [origin,site,arrival,ray] = dbload(dbName)

%DBLOAD extracts travel times from a database
% [ORIGIN,SITE,ARRIVAL,RAY] = DBLOAD(DBNAME) reads an Antelope database
% DBNAME and creates structures with information needed to explore the
% traveltime and ray coverage of an earthquake catalog. DBNAME must contain
% the following tables: origin, event, assoc, arrival and site Antelope
% database tables.
%
%
% -- OUTPUT ARGUMENTS --
%
% ORIGIN
% This is a minimally populated catalog object (see help catalog) that
% describes the earthquake hypocenters. It includes the fields lat, lon,
% depth (km, positive down), origin time (Matlab datenum format), magnitude
% (Ml) and origin id (orid). One element for each earthquake.
%
% ARRIVAL
% This structure describes the phase arrivals. Fields include station name
% (sta), phase name (iphase), travel time (seconds), timeRes (modeled
% travel time residual. specifically it is the observed travel time minus
% the model predicted travel time) and orid. One element for each pahse
% arrival.
%
% RAY
% This structure crudely describes the ray path coverage. It's elements are
% the same length as the elements of arrival. They include originLat,
% otiginLon, originDepth (km, positive down), siteLat, siteLon, siteElev
% (km, positive down) and flatDist (the horizontal distance in km between
% the origin and and site. RAY contains information that is redundant to
% other structures. it exists to allow faster plotting of ray paths.
%
% STATION
% This structure descibes the stations. Fields include sta (station name),
% lat, lon, elevation. STATION has one element for each station. Only
% stations that have at least one phase arrival are included.
%
%
% CAVEATS
% Only P and S phase arrivals are loaded from the database since the goal
% of the ttimes package is to prepare local/regional scale tomography data.
% Because there are so many ways users might wish to cull the dataset, this
% is best done directly on the database before final application of this
% function. No subsetting functions are included here.
%
% See also ttimes.map ttimes.depth_section ttimes.tt_curve
% ttimes.arrival_histogram ttimes.write_lotos ttimes.do_all

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

% CHECK FOR VALID DATABASE
if ~admin.antelope_exists
   error('This function requires the Antelope toolbox for Matlab');
end

try
    db = dbopen(dbName,'r');
    dbclose(db);
catch
    error(['Could not open database: ' dbName]);
end

% FOR DEBUGGING
%  dbName = '/Users/west/FIELDDATA/doc/database/origin/2009april-2010april/reviewed';
%  [origin,site,arrival,ray] = ttimes.dbload(dbName);

%% ACCESS DATABASE
% CREATE DATABASE VIEW
try
    % Entire database
    db = dbopen(dbName,'r');

    % Event view
    db = dblookup_table(db,'origin');
    db1 = dblookup_table(db,'event');
    db = dbjoin(db,db1);
    db = dbsubset(db,'orid==prefor');
    nEvents = dbquery(db,'dbRECORD_COUNT');
    oridList = dbgetv(db,'orid');
    %disp(['Number of origins: ' num2str(nEvents) ]);
    [tmp.lat,tmp.lon,tmp.depth,tmp.time,tmp.ml,tmp.orid,tmp.nass] = dbgetv(db,'origin.lat','origin.lon','origin.depth','origin.time','origin.ml','origin.orid','origin.nass');
    tmp.time = epoch2datenum(tmp.time);
    f = find(tmp.ml < -2);
    tmp.ml(f)=-2;
    f = find(tmp.nass < 0);
    tmp.nass(f)=0;
    origin = catalog;
    origin = set(origin,'LAT',tmp.lat,'LON',tmp.lon,'DEPTH',tmp.depth,'DNUM',tmp.time,'mag',tmp.ml,'ORID',tmp.orid,'NASS',tmp.nass);

    % Arrival view
    db1 = dblookup_table(db,'assoc');
    db = dbjoin(db,db1);
    db1 = dblookup_table(db,'arrival');
    db = dbjoin(db,db1);
    db = dbsubset(db,'iphase=="S" || iphase=="P"');
    db1 = dblookup_table(db,'site');
    db = dbjoin(db,db1);
    nArrivals = dbquery(db,'dbRECORD_COUNT');
    %disp(['Number of arrivals: ' num2str(nArrivals) ]);
    [arrival.sta,arrival.iphase,arrival.time,arrival.originTime,arrival.timeres,arrival.orid] = dbgetv(db,'arrival.sta','arrival.iphase','arrival.time','origin.time','timeres','orid');
    arrival.time = epoch2datenum(arrival.time);
    arrival.originTime = epoch2datenum(arrival.originTime);
    arrival.travelTime = 86400*(arrival.time-arrival.originTime);
    f = find(arrival.timeres==-999);
    arrival.timeres(f) = NaN;
    arrival.predTravelTime = arrival.travelTime + arrival.timeres;

    % Ray view
    [ray.originLat,ray.originLon,ray.originDepth,ray.siteLat,ray.siteLon,ray.siteElev] = dbgetv(db,'origin.lat','origin.lon','origin.depth','site.lat','site.lon','site.elev');
    ray.siteElev = -1*ray.siteElev;
    ray.flatDist = deg2km(distance(ray.originLat,ray.originLon,ray.siteLat,ray.siteLon));

    % Site view
    [tmp.sta,tmp.lon,tmp.lat,tmp.elev] = dbgetv(db,'site.sta','site.lon','site.lat','site.elev');
    tmp.elev = -1*tmp.elev;
    site = [];
    [site.sta,m,~] = unique(tmp.sta);
    site.lat = tmp.lat(m);
    site.lon = tmp.lon(m);
    site.elev = tmp.elev(m);

catch
    error(['Problems with one or more database tables (origin, event, assoc, arrival, site)']);
end
dbclose(db);

% % CREATE TRAVELTIME
% for n = 1:numel(origin.orid)
%     f = find(arrival.orid==origin.orid(n));
%     arrival.travelTime(f) = 86400*(arrival.time(f)-origin.dnum(n));
% end
% arrival.travelTime = arrival.travelTime';
% arrival.predTravelTime = arrival.travelTime + arrival.timeres;

author: geoscience-community-codes

repo: GISMO


text:
function [ img_ihs] = rgb2ihs( img_rgb )
r=img_rgb(:,:,1);
g=img_rgb(:,:,2);
b=img_rgb(:,:,3);

% rgbתihs
I=1/sqrt(3).*r+1/sqrt(3).*g+1/sqrt(3).*b;
v1=1/sqrt(6).*r+1/sqrt(6).*g-2/sqrt(6).*b;
v2=1/sqrt(2).*r-1/sqrt(2).*g;
img_ihs=cat(3,I,v1,v2);
end

author: thfylsty


text:
function [lambda_vec, error_train, error_val] = ...
    validationCurve(X, y, Xval, yval)
%VALIDATIONCURVE Generate the train and validation errors needed to
%plot a validation curve that we can use to select lambda
%   [lambda_vec, error_train, error_val] = ...
%       VALIDATIONCURVE(X, y, Xval, yval) returns the train
%       and validation errors (in error_train, error_val)
%       for different values of lambda. You are given the training set (X,
%       y) and validation set (Xval, yval).
%

% Selected values of lambda (you should not change this)
lambda_vec = [0 0.001 0.003 0.01 0.03 0.1 0.3 1 3 10]';

% You need to return these variables correctly.
error_train = zeros(length(lambda_vec), 1);
error_val = zeros(length(lambda_vec), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return training errors in
%               error_train and the validation errors in error_val. The
%               vector lambda_vec contains the different lambda parameters
%               to use for each calculation of the errors, i.e,
%               error_train(i), and error_val(i) should give
%               you the errors obtained after training with
%               lambda = lambda_vec(i)
%
% Note: You can loop over lambda_vec with the following:
%
%       for i = 1:length(lambda_vec)
%           lambda = lambda_vec(i);
%           % Compute train / val errors when training linear
%           % regression with regularization parameter lambda
%           % You should store the result in error_train(i)
%           % and error_val(i)
%           ....
%
%       end
%
%

for i = 1:length(lambda_vec)
  lambda = lambda_vec(i);

  % Compute train / val errors when training linear
  theta = trainLinearReg(X, y, lambda);

  error_train(i) = linearRegCostFunction(X, y, theta, 0);
  error_val(i) = linearRegCostFunction(Xval, yval, theta, 0);
end

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

end

author: benoitvallon

repo: coursera-machine-learning


text:
% MAKE_TIMEWARP - Select a subset of epochs containing a given event sequence, and return
%                   a matrix of latencies for time warping the selected epochs to a common
%                   timebase in NEWTIMEF. Events in the given sequence may be further
%                   restricted to those with specified event field values.
% Usage:
%       >> timeWarp = make_timewarp(EEG, eventSequence, 'key1',value1,
%       'key2',value2,...);
%
% Inputs:
%   EEG                 - dataset structure
%   eventSequence       - cell array containing a sequence of event type strings.
%                         For example, to select epochs containing a 'movement Onset'
%                         event followed by a 'movement peak', use
%                            {'movement Onset' 'movement peak'}
%                         The output timeWarp matrix will contain the epoch latencies
%                         of the two events for each selected epoch.
%
% Optional inputs (in 'key', value format):
%   'baselineLatency'     - (ms) the minimum acceptable epoch latency for the first
%                         event in the sequence {default: 0}
%   'eventConditions'     - cell array specifying conditions on event fields.
%                         For example, for a sequence consisting of two
%                         events with (velocity) fields vx and vy, use
%                            {vx>0' 'vx>0 && vy<1'}
%                         To accept events unconditionally, use empty strings,
%                         for example {'vx>0','',''} {default: no conditions}
%  'maxSTDForAbsolute'    - (positive number of std. devs.) Remove epochs containing events
%                         whose latencies are more than this number of standard deviations
%                         from the mean in the selected epochs. {default: Inf -> no removal}
%  'maxSTDForRelative'    - (positive number of std. devs.) Remove epochs containing inter-event
%                         latency differences larger or smaller than this number of standard deviations
%                         from the mean in the selected epochs. {default: Inf -> no removal}
% Outputs:
%   timeWarp            - a structure with latencies (time-warp matrix with fields
%                          timeWarp.latencies - an (N, M) timewarp matrix for use in NEWTIMEF
%                           where    N = number of selected epochs containing the specified sequence,
%                                    M = number of events in specified event sequence.
%                          timeWarp.epochs - a (1, M) vector giving the index of epochs with the
%                              specified sequence. Only these epochs should be passed to NEWTIMEF.
%                          timeWarp.eventSequence - same as the 'eventSequence' input variable.
% Example:
%   % Create a timewarp matrix for a sequence of events, first an event of type 'movement Onset' followed by
%   % a 'movement peak' event, with event fields vx and vy:
%
%   >> timeWarp = make_timewarp(EEG, {'movement Onset' 'movement%   peak'},'baselineLatency', 0 ...
%                  ,'eventConditions', {'vx>0' 'vx>0 && vy<1'});
%
%   % To remove events with latencies more than 3 standard deviations from the mean OR 2 std. devs.
%   % from the mean inter-event difference:
%
%   >> timeWarp = make_timewarp(EEG, {'movement Onset' 'movement peak'},'baselineLatency', 0, ...
%                  'eventConditions', {'vx>0' 'vx>0 && vy<1'},'maxSTDForAbsolute', 3,'maxSTDForRelative', 2);
%
% Author: Nima Bigdely Shamlo, SCCN/INC/UCSD, 2008
% See also: SHOW_EVENTS, NEWTIMEF

% Copyright (C) Nima Bigdely Shamlo, SCCN/INC/UCSD, 2008
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function timeWarpStructure = make_timewarp(EEG, eventSequence, varargin)

inputKeyValues = finputcheck(varargin, ...
   {'baselineLatency'     'real'    []          0; ...
    'eventConditions'     'cell'    {}        {}  ; ...
    'maxSTDForAbsolute'        'real'    [0 Inf]          Inf; ...
    'maxSTDForRelative'        'real'    [0 Inf]         Inf...
    });

baselineLatency = [];
maxSTDForAbsolute = 0;
maxSTDForRelative = 0;
eventConditions = [];

% place key values into function workspace variables
inputKeyValuesFields = fieldnames(inputKeyValues);
for i=1:length(inputKeyValuesFields)
    eval([inputKeyValuesFields{i} '= inputKeyValues.' inputKeyValuesFields{i} ';']);
end

if length(eventConditions) < length(eventSequence)
    for i = (length(eventConditions)+1):length(eventSequence)
        eventConditions{i} = '';
    end
end

epochsIsAcceptable = ones(1, length(EEG.epoch));
for epochNumber = 1:length(EEG.epoch)
    eventNameID = 1;
    minimumLatency = baselineLatency;
    timeWarp{epochNumber} = [];
    while eventNameID <= length(eventSequence) % go tthrought event names and find the event that comes after a certain event with correct type and higher latency

        firstLatency = eventsOfCertainTypeAfterCertainLatencyInEpoch(EEG.epoch(epochNumber), eventSequence{eventNameID}, minimumLatency, eventConditions{eventNameID});
        if isempty(firstLatency) % means there were no events, so the epoch is not acceptable
            break;
        else
            timeWarp{epochNumber} = [timeWarp{epochNumber}; firstLatency];
            minimumLatency = firstLatency;
            eventNameID = eventNameID + 1;
        end
    end

    if length(timeWarp{epochNumber}) < length(eventSequence)
        epochsIsAcceptable(epochNumber) = false;
    end
end

acceptableEpochs = find(epochsIsAcceptable);

if isempty(acceptableEpochs)
    timeWarp = {}; % no epoch meet the criteria
else
    timeWarp = cell2mat(timeWarp(acceptableEpochs));

    rejectedEpochesBasedOnLateny = union_bc(rejectEventsBasedOnAbsoluteLatency(timeWarp), rejectEventsBasedOnRelativeLatency(timeWarp));
    timeWarp(:,rejectedEpochesBasedOnLateny) = [];
    acceptableEpochs(rejectedEpochesBasedOnLateny) = [];

end

% since latencies and accepted epochs always have to be together, we put them in one structure
timeWarpStructure.latencies = timeWarp'; % make it suitable for newtimef()

if isempty(timeWarpStructure.latencies) % when empty, it becomes a {} instead of [], so we change it to []
    timeWarpStructure.latencies = [];
end

timeWarpStructure.epochs = acceptableEpochs;
timeWarpStructure.eventSequence = eventSequence;

    function rejectedBasedOnLateny =  rejectEventsBasedOnRelativeLatency(timeWarp)
        % remeve epochs in which the time warped event is further than n
        % standard deviations to mean of latency distance between events
        timeWarpDiff =  diff(timeWarp);
        rejectedBasedOnLateny = [];
        for eventNumber = 1:size(timeWarpDiff, 1)
            rejectedBasedOnLateny = [rejectedBasedOnLateny find(abs(timeWarpDiff(eventNumber, :)-  mean(timeWarpDiff(eventNumber, :))) > maxSTDForRelative * std(timeWarpDiff(eventNumber, :)))];
        end
        rejectedEpochesBasedOnLateny = unique_bc(rejectedBasedOnLateny);
    end

    function rejectedBasedOnLateny =  rejectEventsBasedOnAbsoluteLatency(timeWarp)
        % remeve instances in which the time warped event is further than n
        % standard deviations to mean
        rejectedBasedOnLateny = [];
        for eventNumber = 1:size(timeWarp, 1)
            rejectedBasedOnLateny = [rejectedBasedOnLateny find(abs(timeWarp(eventNumber, :)-  mean(timeWarp(eventNumber, :))) >  maxSTDForAbsolute* std(timeWarp(eventNumber, :)))];
        end
        rejectedEpochesBasedOnLateny = unique_bc(rejectedBasedOnLateny);
    end

    function [firstLatency resultEventNumbers] = eventsOfCertainTypeAfterCertainLatencyInEpoch(epoch, certainEventType, certainLatency, certainCondition)
        resultEventNumbers = [];
        firstLatency = []; % first event latency that meets the critria
        for eventNumber = 1:length(epoch.eventtype)
%            if strcmp(certainEventType, epoch.eventtype(eventNumber)) && epoch.eventlatency{eventNumber} >= certainLatency && eventMeetsCondition(epoch, eventNumber, certainCondition)
             if  eventIsOfType(epoch.eventtype(eventNumber), certainEventType)  && epoch.eventlatency{eventNumber} >= certainLatency && eventMeetsCondition(epoch, eventNumber, certainCondition)
                resultEventNumbers = [resultEventNumbers eventNumber];

                if isempty(firstLatency)
                    firstLatency = epoch.eventlatency{eventNumber};
                end
            end
        end
    end

    function result = eventMeetsCondition(epoch, eventNumber, condition)
        if strcmp(condition,'') || strcmp(condition,'true')
            result = true;
        else

            % get the name thatis before themfor field in the epoch, then remove 'event' name
            epochField = fieldnames(epoch);
            for i=1:length(epochField)
                epochField{i} = strrep(epochField{i},'event',''); % remove event from the beginning of field names
                condition = strrep(condition, epochField{i}, ['cell2mat(epoch.event' epochField{i} '(' num2str(eventNumber) '))' ]);
            end

            result = eval(condition);
        end

    end

    function result = eventIsOfType(eventStr, types)
        % if events are numbers, turn them into strings before comparison
        if ischar(types)
            if iscell(eventStr) && isnumeric(cell2mat(eventStr))
                eventStr = num2str(cell2mat(eventStr));
            end

            result = strcmp(eventStr, types);
        else % it must be a cell of strs
            result = false;
            for i=1:length(types)
                result = result ||  strcmp(eventStr, types{i});
            end
        end
    end

end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function score = HV(Population,optimum)
% <max> <multi/many> <real/integer/label/binary/permutation> <large/none> <constrained/none> <expensive/none> <multimodal/none> <sparse/none> <dynamic/none> <robust/none>
% Hypervolume

%------------------------------- Reference --------------------------------
% E. Zitzler and L. Thiele, Multiobjective evolutionary algorithms: A
% comparative case study and the strength Pareto approach, IEEE
% Transactions on Evolutionary Computation, 1999, 3(4): 257-271.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    PopObj = Population.best.objs;
    if size(PopObj,2) ~= size(optimum,2)
        score = nan;
    else
        [N,M]  = size(PopObj);
        fmin   = min(min(PopObj,[],1),zeros(1,M));
        fmax   = max(optimum,[],1);
        PopObj = (PopObj-repmat(fmin,N,1))./repmat((fmax-fmin)*1.1,N,1);
        PopObj(any(PopObj>1,2),:) = [];
        RefPoint = ones(1,M);
        if isempty(PopObj)
            score = 0;
        elseif M < 4
            % Calculate the exact HV value
            pl = sortrows(PopObj);
            S  = {1,pl};
            for k = 1 : M-1
                S_ = {};
                for i = 1 : size(S,1)
                    Stemp = Slice(cell2mat(S(i,2)),k,RefPoint);
                    for j = 1 : size(Stemp,1)
                        temp(1) = {cell2mat(Stemp(j,1))*cell2mat(S(i,1))};
                        temp(2) = Stemp(j,2);
                        S_      = Add(temp,S_);
                    end
                end
                S = S_;
            end
            score = 0;
            for i = 1 : size(S,1)
                p     = Head(cell2mat(S(i,2)));
                score = score + cell2mat(S(i,1))*abs(p(M)-RefPoint(M));
            end
        else
            % Estimate the HV value by Monte Carlo estimation
            SampleNum = 1e6;
            MaxValue  = RefPoint;
            MinValue  = min(PopObj,[],1);
            Samples   = unifrnd(repmat(MinValue,SampleNum,1),repmat(MaxValue,SampleNum,1));
            for i = 1 : size(PopObj,1)
                drawnow('limitrate');
                domi = true(size(Samples,1),1);
                m    = 1;
                while m <= M && any(domi)
                    domi = domi & PopObj(i,m) <= Samples(:,m);
                    m    = m + 1;
                end
                Samples(domi,:) = [];
            end
            score = prod(MaxValue-MinValue)*(1-size(Samples,1)/SampleNum);
        end
    end
end

function S = Slice(pl,k,RefPoint)
    p  = Head(pl);
    pl = Tail(pl);
    ql = [];
    S  = {};
    while ~isempty(pl)
        ql  = Insert(p,k+1,ql);
        p_  = Head(pl);
        cell_(1,1) = {abs(p(k)-p_(k))};
        cell_(1,2) = {ql};
        S   = Add(cell_,S);
        p   = p_;
        pl  = Tail(pl);
    end
    ql = Insert(p,k+1,ql);
    cell_(1,1) = {abs(p(k)-RefPoint(k))};
    cell_(1,2) = {ql};
    S  = Add(cell_,S);
end

function ql = Insert(p,k,pl)
    flag1 = 0;
    flag2 = 0;
    ql    = [];
    hp    = Head(pl);
    while ~isempty(pl) && hp(k) < p(k)
        ql = [ql;hp];
        pl = Tail(pl);
        hp = Head(pl);
    end
    ql = [ql;p];
    m  = length(p);
    while ~isempty(pl)
        q = Head(pl);
        for i = k : m
            if p(i) < q(i)
                flag1 = 1;
            else
                if p(i) > q(i)
                    flag2 = 1;
                end
            end
        end
        if ~(flag1 == 1 && flag2 == 0)
            ql = [ql;Head(pl)];
        end
        pl = Tail(pl);
    end
end

function p = Head(pl)
    if isempty(pl)
        p = [];
    else
        p = pl(1,:);
    end
end

function ql = Tail(pl)
    if size(pl,1) < 2
        ql = [];
    else
        ql = pl(2:end,:);
    end
end

function S_ = Add(cell_,S)
    n = size(S,1);
    m = 0;
    for k = 1 : n
        if isequal(cell_(1,2),S(k,2))
            S(k,1) = {cell2mat(S(k,1))+cell2mat(cell_(1,1))};
            m = 1;
            break;
        end
    end
    if m == 0
        S(n+1,:) = cell_(1,:);
    end
    S_ = S;
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


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

% problem 3 - causality of y(t)=x(t/4)

%not causal
t1=-5:.1:-1;
x1=zeros(size(t1));
t2=-1:.1:2;
x2=ones(size(t2));
t3=2:.1:10;
x3=zeros(size(t3));
t=[t1 t2 t3];
x=[x1 x2 x3];
plot(t,x);
ylim([-0.1 1.1]);
legend('x(t)')

figure
plot(4*t,x);
ylim([-0.1 1.1]);
legend('y(t)')


text:
function [qn,QNq] = normquat(q)

% NORMQUAT Normalize quaternion to unit length
%   NORMQUAT(Q) returns a unit length quaternion Q/norm(Q)
%
%   [qn,QNq] = NORMQUAT(Q) returns also the Jacobian wrt Q. Note that this
%   Jacobian is a symmetric 4x4 matrix.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

nq = sqrt(q(:)'*q(:));
qn = q/nq;

if nargout > 1

    a = q(1);
    b = q(2);
    c = q(3);
    d = q(4);

    nq3 = nq^3;

    QNq = [...
        [ (b^2+c^2+d^2)/nq3,          -a/nq3*b,          -a/nq3*c,          -a/nq3*d]
        [          -a/nq3*b, (a^2+c^2+d^2)/nq3,          -b/nq3*c,          -b/nq3*d]
        [          -a/nq3*c,          -b/nq3*c, (a^2+b^2+d^2)/nq3,          -c/nq3*d]
        [          -a/nq3*d,          -b/nq3*d,          -c/nq3*d, (a^2+b^2+c^2)/nq3]];
end
return

%%

syms a b c d real
q = [a;b;c;d];
[qn,QNq] = normquat(q);

QNq - simple(jacobian(qn,q))

QNq - QNq'

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function imdb = cnn_setup_data_voc07_ssw(varargin)
%CNN_SETUP_DATE_VOC07_SSW  Setup PASCAL VOC 2007 data with precomputed SSW proposals
%
% Copyright (C) 2016 Hakan Bilen.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

opts.sswDir = fullfile('data','SSW') ;
opts.dataDir = 'data' ;
opts.addFlipped = true ;
opts.useDifficult = true ;
opts.fgThresh = 0.5 ;
opts.bgThreshHi = 0.5 ;
opts.bgThreshLo = 0.1 ;
opts = vl_argparse(opts, varargin) ;

%% Get selective search windows
files = {...
  'SelectiveSearchVOC2007trainval.mat', ...
  'SelectiveSearchVOC2007test.mat'} ;

if ~exist(opts.sswDir, 'dir')
  mkdir(opts.sswDir) ;
end

for i=1:numel(files)
  outPath = fullfile(opts.sswDir, files{i}) ;
  if ~exist(outPath, 'file')
    url = sprintf('http://koen.me/research/downloads/%s',files{i}) ;
    fprintf('Downloading %s to %s\n', url, outPath) ;
    urlwrite(url,outPath) ;
  end
end

%% Get image names and gt boxes.
imdb = cnn_setup_data_voc07(...
  'dataDir', opts.dataDir, ...
  'useDifficult', opts.useDifficult, ...
  'addFlipped', opts.addFlipped) ;

%% Load precomputed object proposals
train = load(fullfile(opts.sswDir, files{1})) ;
test = load(fullfile(opts.sswDir, files{2})) ;

% Check that all proposals are correctly matched to an image in imdb
assert(isempty(...
  setxor(vertcat(train.images,test.images), cellfun(@(x)x(1:6),imdb.images.name,'uniformoutput',0)))) ;

if opts.addFlipped
  boxes = horzcat(train.boxes,train.boxes,test.boxes) ;
  imnames = vertcat(train.images,train.images,test.images) ;
else
  boxes = horzcat(train.boxes,test.boxes) ;
  imnames = vertcat(train.images,test.images) ;
end
assert(numel(boxes)==numel(imdb.images.name)) ;

% Sort proposals by ascending image name, as this is also the order in imdb
[imnames,si] = sort(imnames);
boxes = boxes(si);

% Transpose XY in proposals
for i=1:numel(boxes)
  boxes{i} = boxes{i}(:,[2 1 4 3]);
end

% Check boxes and find corresponding image
for i=1:numel(imdb.images.name)

  bbox = boxes{i};
  assert(strcmp(imnames{i},imdb.images.name{i}(1:end-4)));
  assert(all(all(bbox(:,[1 2]) <= bbox(:,[3 4])))) ;
  assert(all(all(1 <= bbox & bsxfun(@le, bbox, imdb.images.size(i,[1 2 1 2]))))) ;

  % Flip the proposals if needed
  if imdb.boxes.flip(i)
    imageWidth = imdb.images.size(i,1) ;
    bbox(:,[3 1]) = imageWidth - bbox(:,[1 3]) + 1 ;
  end

  boxes{i} = bbox ;
end

imdb = attach_proposals(imdb,boxes,opts.fgThresh,opts.bgThreshHi,opts.bgThreshLo);
imdb = add_bboxreg_targets(imdb);

if 0
  for i = 1:100
    for c = 1:21
      sel = find(imdb.boxes.plabel{i} == c) ;
      if isempty(sel), continue ; end
      im = imread(fullfile(imdb.imageDir, imdb.images.name{i})) ;
      if imdb.boxes.flip(i), im = fliplr(im) ; end
      pbox = imdb.boxes.pbox{i}(sel,:) ;
      clf ; bbox_draw(im, pbox) ;
      title(sprintf('Class %d', c)) ;
      keyboard
    end
  end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
% Dynamical model function for the random sine signal demo

% Copyright (C) 2007 Jouni Hartikainen
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.

function x_n = ekf_sine_f(x,param)
    dt = param(1);
    A = [1 dt 0;0 1 0;0 0 1];
    x_n = A*x(1:3,:);
    if size(x,1) == 7 || size(x,1) == 6
        x_n(1:3,:) = x_n(1:3,:) + x(4:6,:);
    end

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


text:
function [vehicle_pose]=apply_kinematics(vehicle_pose,vehicle_m,rad_speed,deltaT)

vehicle_pose=vehicle_pose+[vehicle_m*cos(vehicle_pose(3))*deltaT; ...
    vehicle_m*sin(vehicle_pose(3))*deltaT; rad_speed*deltaT];
end

author: mubowen

repo: Baidu-Apollo-control-algorithm


text:
function [ n_data, tc, p, pr ] = prandtl_values ( n_data )

%*****************************************************************************80
%
%% PRANDTL_VALUES returns some values of the Prandtl number.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Lester Haar, John Gallagher and George Kell,
%    NBS/NRC Steam Tables:
%    Thermodynamic and Transport Properties and Computer Programs
%      for Vapor and Liquid States of Water in SI Units,
%    Hemisphere Publishing Corporation, Washington, 1984,
%    TJ270.H3, page 265.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real TC, the temperature, in degrees Celsius.
%
%    Output, real P, the pressure, in bar.
%
%    Output, real PR, the Prandtl number, dimensionless.
%
  n_max = 35;

  pr_vec = [ ...
     13.50E+00, ...
     13.48E+00, ...
     13.46E+00, ...
     13.39E+00, ...
     13.27E+00, ...
     13.15E+00, ...
     13.04E+00, ...
     12.93E+00, ...
     12.83E+00, ...
     12.73E+00, ...
     12.63E+00, ...
     12.53E+00, ...
     12.43E+00, ...
     12.34E+00, ...
     12.25E+00, ...
     12.08E+00, ...
     11.92E+00, ...
     11.77E+00, ...
     11.62E+00, ...
     11.48E+00, ...
     11.36E+00, ...
     11.23E+00, ...
     11.12E+00, ...
     10.91E+00, ...
     10.72E+00, ...
     10.55E+00, ...
      6.137E+00, ...
      3.555E+00, ...
      2.378E+00, ...
      1.000E+00, ...
      0.974E+00, ...
      0.960E+00, ...
      0.924E+00, ...
      0.899E+00, ...
      0.882E+00 ];

  p_vec = [ ...
        1.0E+00, ...
        5.0E+00, ...
       10.0E+00, ...
       25.0E+00, ...
       50.0E+00, ...
       75.0E+00, ...
      100.0E+00, ...
      125.0E+00, ...
      150.0E+00, ...
      175.0E+00, ...
      200.0E+00, ...
      225.0E+00, ...
      250.0E+00, ...
      275.0E+00, ...
      300.0E+00, ...
      350.0E+00, ...
      400.0E+00, ...
      450.0E+00, ...
      500.0E+00, ...
      550.0E+00, ...
      600.0E+00, ...
      650.0E+00, ...
      700.0E+00, ...
      800.0E+00, ...
      900.0E+00, ...
     1000.0E+00, ...
        1.0E+00, ...
        1.0E+00, ...
        1.0E+00, ...
        1.0E+00, ...
        1.0E+00, ...
        1.0E+00, ...
        1.0E+00, ...
        1.0E+00, ...
        1.0E+00 ];

  tc_vec = [ ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
       0.0E+00, ...
      25.0E+00, ...
      50.0E+00, ...
      75.0E+00, ...
     100.0E+00, ...
     150.0E+00, ...
     200.0E+00, ...
     400.0E+00, ...
     600.0E+00, ...
     800.0E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    tc = 0.0;
    p = 0.0;
    pr = 0.0;
  else
    tc = tc_vec(n_data);
    p = p_vec(n_data);
    pr = pr_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef (Abstract) net_model_ac < mp.net_model% & mp.form_ac
%MP.NET_MODEL_AC Abstract class, explicitly a subclass of MP.NET_MODEL and
%              implicitly assumed to be a subclass of MP.FORM_AC as well

%   MATPOWER
%   Copyright (c) 2019-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

    properties
        zr = [];
        zi = [];
    end

    properties (Access=protected)
        inln_list = {};         %% list of indexes of nme's w/inln
        snln_list = {};         %% list of indexes of nme's w/snln
        inln_hess_list = {};    %% list of indexes of nme's w/inln_hess
        snln_hess_list = {};    %% list of indexes of nme's w/snln_hess
    end

    methods
        function obj = def_set_types(obj)
            def_set_types@mp.net_model(obj);    %% call parent first
            obj.set_types.zr = 'NON-VOLTAGE VARS REAL (zr)';
            obj.set_types.zi = 'NON-VOLTAGE VARS IMAG (zi)';
        end

        function obj = build_params(obj, nm, dm)
            %% call parent to build individual element parameters
            build_params@mp.net_model(obj, nm, dm);

            %% aggregate parameters from individual elements
            obj.Y = obj.stack_matrix_params('Y', 1);
            obj.L = obj.stack_matrix_params('L', 0);
            obj.M = obj.stack_matrix_params('M', 1);
            obj.N = obj.stack_matrix_params('N', 0);
            obj.i = obj.stack_vector_params('i');
            obj.s = obj.stack_vector_params('s');

            %% add general nonlinear function if any element has one defined
            for k = 1:length(obj.elements)
                if ~isempty(obj.elements{k}.inln)
                    obj.inln_list{end+1} = k;
                    if ~isempty(obj.elements{k}.inln_hess)
                        obj.inln_hess_list{end+1} = k;
                    end
                end
                if ~isempty(obj.elements{k}.snln)
                    obj.snln_list{end+1} = k;
                    if ~isempty(obj.elements{k}.snln_hess)
                        obj.snln_hess_list{end+1} = k;
                    end
                end
            end
            if ~isempty(obj.inln_list)
                obj.inln = @(x_, sysx, idx)port_inj_nln(obj, 'i', x_, sysx, idx);
                if ~isempty(obj.inln_hess_list)
                    obj.inln_hess = @(x_, lam, sysx, idx)port_inj_nln_hess(obj, 'i', x_, lam, sysx, idx);
                end
            end
            if ~isempty(obj.snln_list)
                obj.snln = @(x_, sysx, idx)port_inj_nln(obj, 's', x_, sysx, idx);
                if ~isempty(obj.snln_hess_list)
                    obj.snln_hess = @(x_, lam, sysx, idx)port_inj_nln_hess(obj, 's', x_, lam, sysx, idx);
                end
            end
        end

        function [g, gv1, gv2, gzr, gzi] = port_inj_nln(obj, si, x_, sysx, idx)
            if nargin < 5
                idx = [];
                if nargin < 4
                    sysx = 1;
                end
            end

            %% current or power
            fcn = [si 'nln'];
            fcn_list = [fcn '_list'];

            %% initialize
            if isempty(idx)
                sel = 0;        %% all ports
                np = obj.np;
            else
                sel = 1;        %% selected ports only
                np = length(idx);
            end
            nv = obj.get_nv_(sysx);
            nz = obj.nz;
            nc = size(x_, 2);   %% num of cols in x_, for evaluating multiple x_
            g = zeros(np, nc);
            gv1 = sparse(np, nv);
            gv2 = sparse(np, nv);
            gzr = sparse(np, nz);
            gzi = sparse(np, nz);

            %% loop through elements w/gen nonlin fcns, evaluate them
            pp = obj.get_idx('port');
            if ~sysx
                ss = obj.get_idx('state');
            end
            for kk = obj.(fcn_list)
                k = kk{1};      %% index into obj.elements
                nme = obj.elements{k};
                i1 = pp.i1.(nme.name)(1);
                iN = pp.iN.(nme.name)(end);

                %% set up port index vector for nme
                if sel
                    apidx = find(idx >= i1 & idx <= iN);   %% aggregate port indices in range
                    if isempty(apidx)  %% skip if selected ports, but none in range
                        continue;
                    end
                    nme_idx = idx(apidx) - i1 + 1;  %% port index vector for nme
                else
                    nme_idx = [];                   %% all ports for nme
                end

                %% set up proper x_ for nme
                if sysx
                    nme_x_ = x_;
                else
                    if isfield(ss.i1, nme.name)
                        j1 = ss.i1.(nme.name)(1);
                        jN = ss.iN.(nme.name)(end);
                    else
                        j1 = 1;
                        jN = 0;
                    end
                    nme_x_ = [  x_(i1:iN, :);
                                x_(nv+j1:nv+jN, :)  ];
                end

                %% call nonlinear function
                gg = cell(1, nargout);
                [gg{:}] = nme.(fcn)(nme_x_, sysx, nme_idx);

                %% insert the results in aggregate output args
                if sel
                    g(apidx, :) = gg{1};
                    if nargout > 1
                        if sysx
                            gv1(apidx, :) = gg{2};
                            gv2(apidx, :) = gg{3};
                            if nargout > 3 && nme.nz
                                gzr(apidx, :) = gg{4};
                                gzi(apidx, :) = gg{5};
                            end
                        else
                            gv1(apidx, i1:iN) = gg{2};
                            gv2(apidx, i1:iN) = gg{3};
                            if nargout > 3 && nme.nz
                                gzr(apidx, j1:jN) = gg{4};
                                gzi(apidx, j1:jN) = gg{5};
                            end
                        end
                    end
                else
                    g(i1:iN, :) = gg{1};
                    if nargout > 1
                        if sysx
                            gv1(i1:iN, :) = gg{2};
                            gv2(i1:iN, :) = gg{3};
                            if nargout > 3 && nme.nz
                                gzr(i1:iN, :) = gg{4};
                                gzi(i1:iN, :) = gg{5};
                            end
                        else
                            gv1(i1:iN, i1:iN) = gg{2};
                            gv2(i1:iN, i1:iN) = gg{3};
                            if nargout > 3 && nme.nz
                                gzr(i1:iN, j1:jN) = gg{4};
                                gzi(i1:iN, j1:jN) = gg{5};
                            end
                        end
                    end
                end
            end     %% for loop
        end

        function H = port_inj_nln_hess(obj, si, x_, lam, sysx, idx)
            % H = obj.port_inj_power_hess(x_, lam)
            % H = obj.port_inj_power_hess(x_, lam, sysx)
            % H = obj.port_inj_power_hess(x_, lam, sysx, idx)
            if nargin < 6
                idx = [];
                if nargin < 5
                    sysx = 1;
                end
            end

            %% current or power
            fcn = [si 'nln_hess'];
            fcn_list = [fcn '_list'];

            %% initialize
            n = 2 * length(x_);
            H = sparse(n, n);

            %% loop through elements w/gen nonlin Hessians, evaluate them
            pp = obj.get_idx('port');
            if ~sysx
                ss = obj.get_idx('state');
            end
            for kk = obj.(fcn_list)
                k = kk{1};      %% index into obj.elements
                nme = obj.elements{k};
                i1 = pp.i1.(nme.name)(1);
                iN = pp.iN.(nme.name)(end);

                %% set up x_ for nme & corresp row/col indices for nme
                if sysx
                    nme_x_ = x_;
                else
                    nv = obj.get_nv_(sysx);
                    nz = obj.nz;
                    if isfield(ss.i1, nme.name)
                        j1 = ss.i1.(nme.name)(1);
                        jN = ss.iN.(nme.name)(end);
                    else
                        j1 = 1;
                        jN = 0;
                    end
                    nme_x_ = [  x_(i1:iN, :);
                                x_(nv+j1:nv+jN, :)  ];

                    %% indices of rows/cols of H corresponding to nme x_
                    h = [(i1:iN) nv+(i1:iN) 2*nv+(j1:jN) 2*nv+nz+(j1:jN)].';
                end

                %% set up port index and lambda vectors for nme
                if ~isempty(idx)    %% selected ports only
                    apidx = find(idx >= i1 & idx <= iN);    %% aggregate port indices in range
                    if isempty(apidx)  %% skip if selected ports, but none in range
                        continue;
                    end
                    nme_idx = idx(apidx) - i1 + 1;  %% port index vector for nme
                    nme_lam = lam(apidx);           %% corresponding lam
                else                %% all ports
                    nme_idx = [];
                    nme_lam = lam(i1:iN);
                end

                %% call nonlinear function
                nme_H = nme.(fcn)(nme_x_, nme_lam, sysx, nme_idx);

                %% accumulate output
                if sysx
                    H = H + nme_H;
                else
                    H(h,h) = H(h,h) + nme_H;
                end
            end
        end

        function [G, Gv1, Gv2, Gzr, Gzi] = nodal_complex_current_balance(obj, x_)
            %% node incidence matrix
            C = obj.C;

            %% get port current injections with derivatives
            if nargout > 1
                [I, Iv1, Iv2, Izr, Izi] = obj.port_inj_current(x_, 1);
                Gv1 = C * Iv1;      %% Gva or Gvr
                Gv2 = C * Iv2;      %% Gvm or Gvi
                Gzr = C * Izr;
                Gzi = C * Izi;
            else
                I = obj.port_inj_current(x_, 1);
            end

            %% nodal current balance
            G = C * I;
        end

        function [G, Gv1, Gv2, Gzr, Gzi] = nodal_complex_power_balance(obj, x_)
            %% node incidence matrix
            C = obj.C;

            %% get port power injections with derivatives
            if nargout > 1
                [S, Sv1, Sv2, Szr, Szi] = obj.port_inj_power(x_, 1);
                Gv1 = C * Sv1;      %% Gva or Gvr
                Gv2 = C * Sv2;      %% Gvm or Gvi
                Gzr = C * Szr;
                Gzi = C * Szi;
            else
                S = obj.port_inj_power(x_, 1);
            end

            %% nodal power balance
            G = C * S;
        end

        function d2G = nodal_complex_current_balance_hess(obj, x_, lam)
            %% get port power injection hessians
            d2G = obj.port_inj_current_hess(x_, obj.C' * lam);
        end

        function d2G = nodal_complex_power_balance_hess(obj, x_, lam)
            %% get port power injection hessians
            d2G = obj.port_inj_power_hess(x_, obj.C' * lam);
        end

        function obj = port_inj_soln(obj)
            %% compute port injections
%             obj.soln.gi_ = obj.port_inj_current(obj.soln.x);
            obj.soln.gs_ = obj.port_inj_power(obj.soln.x);
        end

        function va = get_va(obj, idx)
            if isfield(obj.soln, 'v')           %% solved value
                if nargin < 2 || isempty(idx)
                    va = angle(obj.soln.v);
                else
                    va = angle(obj.soln.v(idx));
                end
            else                                %% initial value
                if nargin < 2 || isempty(idx)
                    va = obj.initial_voltage_angle();
                else
                    va = obj.initial_voltage_angle(idx);
                end
            end
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function varargout = spm_bsplins(varargin)
% Sample a volume using B-spline interpolation - a compiled routine
% FORMAT [f,dfx,dfy,dfz] = spm_bsplins(c,x,y,z,d)
%   c - volume of B-spline coefficients (from spm_bsplinc)
%   x,y,z - co-ordinates of sampled points
%       d(1:3) - degree of B-spline (from 0 to 7) along different dimensions
%                - these must be same as used by spm_bsplinc
%       d(4:6) - 1/0 to indicate wrapping along the dimensions
%   f - sampled data
%   dfx,dfy,dfz - sampled first derivatives
%
% This function takes B-spline basis coefficients from spm_bsplinc,
% and re-convolves them with B-splines centred at the new sample points.
%
% Note that nearest neighbour interpolation is used instead of 0th
% degree B-splines, and the derivatives of trilinear interpolation are
% returned insted of those of 1st degree B-splines.  The difference is
% extremely subtle.
%
%_______________________________________________________________________
%
% References:
%   M. Unser, A. Aldroubi and M. Eden.
%   "B-Spline Signal Processing: Part I-Theory,"
%   IEEE Transactions on Signal Processing 41(2):821-832 (1993).
%
%   M. Unser, A. Aldroubi and M. Eden.
%   "B-Spline Signal Processing: Part II-Efficient Design and
%   Applications,"
%   IEEE Transactions on Signal Processing 41(2):834-848 (1993).
%
%   M. Unser.
%   "Splines: A Perfect Fit for Signal and Image Processing,"
%   IEEE Signal Processing Magazine, 16(6):22-38 (1999)
%
%   P. Thevenaz and T. Blu and M. Unser.
%   "Interpolation Revisited"
%   IEEE Transactions on Medical Imaging 19(7):739-758 (2000).
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_bsplins.m 2696 2009-02-05 20:29:48Z guillaume $

%-This is merely the help file for the compiled routine
error('spm_bsplins.c not compiled - see Makefile');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function exp_endfun(varargin)
% End the definition of an expression-aware function.
% exp_endfun(Arguments...)
%
% See exp_beginfun. Invokes plugin functions that shall run after the body of the expression-aware
% function (here called poststeps) and handles memoization.
%
% In:
%   Arguments... : optional name-value pairs; some options allowed in exp_beginfun can also be
%                  overridden here: 'set_online', 'append_online', 'add_impure'
%
% Notes:
%   The exp_endfun must be called in every code path that exits the function being defined (except
%   for error conditions), i.e. a premature return must be preceded by and exp_endfun;
%
% See also:
%   exp_beginfun, exp_eval
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-15
dp;

% retrieve the context structure
% If you get an "Undefined function or variable" error here then you forgot to call exp_beginfun
% at the beginning of the same function in which exp_endfun was called
context = evalin('caller','exp_internal_context');
if ~context.enabled
    return; end

% update the new function input & output arguments
for i=1:length(context.outargs)
    try context.ws_output_post.(context.outargs{i}) = evalin('caller',context.outargs{i}); catch, end
end

% if the arg system is used by the function, we have an easier-to-handle version of the posteval expression
if evalin('caller','exist(''arg_nvps'',''var'')')
    context.expression_posteval.parts = evalin('caller','arg_nvps'); end

% assign additional exp_endfun options to the context's options...
if ~iscellstr(varargin(1:2:end))
    error('The given arguments must be name-value pairs.'); end
for k=1:2:length(varargin)
    context.opts.(varargin{k}) = varargin{k+1}; end

if context.opts.add_impure
    % if impure, incorporate the original unevaluated expression into the function result (if it's a struct)
    if isfield(context.ws_output_post,context.outargs{1}) && isstruct(context.ws_output_post.(context.outargs{1}))
        context.ws_output_post.(context.outargs{1}).tracking.expression = context.expression_preeval; end
end

% run the poststep(s) on the workspaces (add online expression, etc.)
for s = 1:length(context.opts.poststeps)
    context = context.opts.poststeps{s}(context); end

% handle fingerprinting of the result
if isempty(context.opts.fingerprint_create)
    context.opts.fingerprint_create = hlp_resolve('fingerprint_create',true,context.exec_ctx); end
if strcmp(context.opts.fingerprint_create,'passthrough')
    % no changes are being made to the fingerprint of the output
elseif isfield(context.ws_output_post,context.outargs{1}) && is_impure_expression(context.ws_output_post.(context.outargs{1}))
    % we have an impure expression as output
    if ~isequal(false,context.opts.fingerprint_create) && ~isequal(false,exp_eval(utl_replacerepeated(context.opts.fingerprint_create,{exp_rule(@expression,context.expression_preeval)}),inf)) ...
        % fingerprint creation turned on: fingerprint the data, so we can later check the consistency of impure expression
        context.ws_output_post.(context.outargs{1}).tracking.fingerprint = hlp_fingerprint(rmfield(context.ws_output_post.(context.outargs{1}),'tracking'));
    else
        % fingerprint creation turned off: make sure that the fingerprint gets removed...
        if isfield(context.ws_output_post.(context.outargs{1}).tracking,'fingerprint')
            context.ws_output_post.(context.outargs{1}).tracking = rmfield(context.ws_output_post.(context.outargs{1}).tracking,'fingerprint'); end
    end
end

% remember the computation time of the object
if isfield(context.ws_output_post,context.outargs{1}) && is_impure_expression(context.ws_output_post.(context.outargs{1}))
    context.ws_output_post.(context.outargs{1}).tracking.computation_time = toc(context.eval_time); end

% assign the changed values in the workspace
for fn=fieldnames(context.ws_output_post)'
    assignin('caller',fn{1},context.ws_output_post.(fn{1})); end

% memoize the results (for the original expression)
if ~isempty(context.memoize_id) && isfield(context.ws_output_post,context.outargs{1})
    obj = {context.ws_output_post.(context.outargs{1})};
    for k=2:length(context.outargs)
        % and also append the other output arguments as far as they are set
        try obj = [obj {context.ws_output_post.(context.outargs{k})}]; catch break; end
    end

    % commit object to cache
    utl_memoize_commit(obj,context.memoize_id,context.input_size);
end

repo: REST


text:
% edgeProcessing processes edges from given image following steps below.
%
%   * Step 1. Edge detection
% Extracts Edge Map from given image and preprocess it to reduce complexity.
%   * Step 2. Edge Contours extraction
% Extracts Edge Contours from given Edge Map.
%   * Step 3. Line segment fitting
% Extracts Line Segments fitting on the Edge Contour.
%   * Step 4. Detecting edge portions with Smooth Curvatures.
% Classify edges using Curvatures so that every edge doesn't have a sudden
% change which is called sharp turn and a inflexion point.
%
% Input     - im : 'jpg' image which taken by normal smart phone.
% Output    - ResizedIm : (possibly) resized image.
%             SmoothCurvatures = The list of edges with Smooth Curvature.
%             Each of them doesn't have any sharp turn and inflexion point.

function [ResizedIm, SmoothCurvatures] = edgeProcessing(im, showFigures)

%% Initialize

% Set default parameters
if nargin < 2
    showFigures = true;
end

% For edge detection.
T_low = 0.2;
T_high = 0.4;
sigma = 1;

% For Line segment fitting.
DeviationThreshold = 2;

% For Smooth Curvatures.
SharpTurnThreshold = 90;

%% Preprocess image

% Resize if input image is too high resolution.
larger_side = max(size(im,1), size(im,2));
max_side = 800;
if larger_side > max_side
    im = imresize(im, max_side/larger_side);
end
ResizedIm = im;

% Convert RGB to gray.
im = rgb2gray(im);
im = histeq(im);

%% Step 1. Edge detection

% Implement Canny edge detection.
EdgeMap = edge(im, 'canny', [T_low, T_high], sigma);

% Apply 'thin' operation to ignore unnecessary parts.
EdgeMap = bwmorph(EdgeMap,'thin',Inf);

% Remove junctions.
[ROW, COL] = size(EdgeMap);
for i = 2:ROW-1
    for j = 2:COL-1
        if EdgeMap(i,j) == 1 && sum(sum(EdgeMap(i-1:i+1,j-1:j+1))) >= 4
            EdgeMap(i-1:i+1,j-1:j+1) = 0;
        end
    end
end

% Remove isolated pixels.
EdgeMap = bwmorph(EdgeMap,'clean');

%% Step 2. Edge Contours extraction

% Extract Edge Contours using preprocessed EdgeMap.
EdgeContours = extractEdgeContours(EdgeMap);

% Draw Edge Contours
if showFigures
    drawLineSegments(size(EdgeMap), EdgeContours);
end

%% Step 3. Line segment fitting

% Extract Line Segments fitting on the Edge Contour.
SegmentList = lineSegmentFitting(EdgeContours, DeviationThreshold);

% Draw Segment Lines
if showFigures
    drawLineSegments(size(EdgeMap), SegmentList);
end

%% Step 4. Detecting edge portions with Smooth Curvatures.

% Classify Smooth Curvatures.
SmoothCurvatures = classifySmoothCurvatures(SegmentList, SharpTurnThreshold);

% Draw Smooth Curvatures.
if showFigures
    drawLineSegments(size(EdgeMap), SmoothCurvatures);
end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm