text:
classdef AnnotationState < uint8
    enumeration
        kZoomInState(1)
        kZoomOutState (2)
        kPanState (3)
        kSetTimelineState (4)
        kAddEventState (5)
        kAddRangeState (6)
        kModifyAnnotationState (7)
        kDeleteAnnotationState (8)
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function writeTopple(IC)

fid = fopen('ToppleFromRest.m','w');

fprintf(fid,'function Critical = ToppleFromRest(P)\n');
fprintf(fid,'\n');
fprintf(fid,'AUTOMATICALLY GENERATED FILE  --  DO ONT EDIT\n');
fprintf(fid,'\n');
fprintf(fid,'%%This function uses conservation of energy to compute the dynamcis as the\n');
fprintf(fid,'%%stick topples from rest. This assumes that toppling will be in the\n');
fprintf(fid,'%%positive theta direction.\n');
fprintf(fid,'\n');
fprintf(fid,'%%%%%%%% Find the critical angle at which slipping occurs:\n');
fprintf(fid,'\n');
fprintf(fid,'options = optimset( ''Display'',''off'',...\n');
fprintf(fid,'                    ''TolX'',1e-14,...\n');
fprintf(fid,'                    ''TolFun'',1e-12);\n');
fprintf(fid,'[th_left,~,flag_left] = fzero(@(th)FrictionConeLeft(th,P),pi/4,options);\n');
fprintf(fid,'[th_right,~,flag_right] = fzero(@(th)FrictionConeRight(th,P),pi/4,options);\n');
fprintf(fid,'\n');
fprintf(fid,'if th_left > th_right\n');
fprintf(fid,'    %Then slips forwards first\n');
fprintf(fid,'    th_crit = th_left;\n');
fprintf(fid,'    flag = flag_left;\n');
fprintf(fid,'     Critical.exit = ''SlipForwards'';\n');
fprintf(fid,'else\n');
fprintf(fid,'    %Then slips backwards first\n');
fprintf(fid,'    th_crit = th_right;\n');
fprintf(fid,'    flag = flag_right;\n');
fprintf(fid,'    Critical.exit = ''SlipBackwards'';\n');
fprintf(fid,'end\n');
fprintf(fid,'\n');
fprintf(fid,'if flag~=1 %Then something went wrong\n');
fprintf(fid,'    Critical.exit = ''Fail'';\n');
fprintf(fid,'    th_crit = [];\n');
fprintf(fid,'    dth_crit = [];\n');
fprintf(fid,'    H_crit = [];\n');
fprintf(fid,'elseif th_crit < 0 || th_crit > pi/2\n');
fprintf(fid,'    Critical.exit = ''OutOfBounds'';\n');
fprintf(fid,'    th_crit = [];\n');
fprintf(fid,'    dth_crit = [];\n');
fprintf(fid,'    H_crit = [];   \n');
fprintf(fid,'    \n');
fprintf(fid,'else\n');
fprintf(fid,'    dth_crit = AngularRate(th_crit,P);\n');
fprintf(fid,'    H_crit = ContactHorizontal(th_crit,P);\n');
fprintf(fid,'end\n');
fprintf(fid,'\n');
fprintf(fid,'Critical.th = th_crit;\n');
fprintf(fid,'Critical.dth = dth_crit;\n');
fprintf(fid,'Critical.H = H_crit;\n');
fprintf(fid,'\n');
fprintf(fid,'end\n');
fprintf(fid,'\n');
fprintf(fid,'\n');
fprintf(fid,'function FCL = FrictionConeLeft(th,P)\n');
fprintf(fid,'\n');
fprintf(fid,'u = P.u;\n');
fprintf(fid,'H = ContactHorizontal(th,P);\n');
fprintf(fid,'V = ContactVertical(th,P);\n');
fprintf(fid,'FCL = H+u*V;\n');
fprintf(fid,'\n');
fprintf(fid,'end\n');
fprintf(fid,'\n');
fprintf(fid,'function FCR = FrictionConeRight(th,P)\n');
fprintf(fid,'\n');
fprintf(fid,'u = P.u;\n');
fprintf(fid,'H = ContactHorizontal(th,P);\n');
fprintf(fid,'V = ContactVertical(th,P);\n');
fprintf(fid,'FCR = H-u*V;\n');
fprintf(fid,'end\n');
fprintf(fid,'\n');
fprintf(fid,'function V = ContactVertical(th,P)\n');
fprintf(fid,'\n');
fprintf(fid,'m = P.m;\n');
fprintf(fid,'g = P.g;\n');
fprintf(fid,'L = P.L;\n');
fprintf(fid,'I = P.I;\n');
fprintf(fid,'\n');
fprintf(fid,['V = ' vectorize(IC.V(1)) ';\n']);
fprintf(fid,'\n');
fprintf(fid,'end\n');
fprintf(fid,'\n');
fprintf(fid,'function H = ContactHorizontal(th,P)\n');
fprintf(fid,'\n');
fprintf(fid,'m = P.m;\n');
fprintf(fid,'g = P.g;\n');
fprintf(fid,'L = P.L;\n');
fprintf(fid,'I = P.I;\n');
fprintf(fid,'\n');
fprintf(fid,['H = ' vectorize(IC.H(1)) ';\n']);
fprintf(fid,'\n');
fprintf(fid,'end\n');
fprintf(fid,'\n');
fprintf(fid,'function dth = AngularRate(th,P)\n');
fprintf(fid,'\n');
fprintf(fid,'m = P.m;\n');
fprintf(fid,'g = P.g;\n');
fprintf(fid,'L = P.L;\n');
fprintf(fid,'I = P.I;\n');
fprintf(fid,'\n');
fprintf(fid,['dth = ' vectorize(IC.dth(1)) ';\n']);
fprintf(fid,'\n');
fprintf(fid,'end\n');

fclose(fid);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
A=rand(1000,1000);
A=A'*A;

tic
B=mexInvSym(A);
t=toc;
fprintf('mex-file time: %fs\n',t);

tic
B2=inv(A);
t=toc;
fprintf('matlab-file time: %fs\n',t);

sum((B(:)-B2(:)).^2)

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function artinis_license

% ARTINIS_LICENSE prints the license only once upon the first call to
% this function. If the user does a "clear all", the license will
% again be shown.  This function should be included in every artinis
% function to ensure that the license is displayed at least once.

% Copyright (C) 2009-2011, Robert Oostenveld
% Copyright (C) 2014-2019, Jörn M. Horschig, Artinis Medical Systems, www.artinis.com

persistent status

if isempty(status)
  % this is for the first time
  status = false;
end

if ~status
  % show the license
  try
    tmp = which('artinis_license');
  catch
    tmp = mfilename;
  end
  [p, f] = fileparts(tmp);
  clc
  fprintf('==============================================================================\n');
  type(fullfile(p, 'LICENSE.txt'));
  fprintf('==============================================================================\n');
  % remember that it has been shown
  status = true;
  pause(3)
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%compute central speed in the tail-central direction
function [data,units]=compute_velcentraltc(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
velcentraltc=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    % this is just slightly faster
    velcentraltc{i} = trx(larva).dxcentral_mm.*cos(trx(larva).tailcentralang(1:end-1)) + ...
      trx(larva).dycentral_mm.*sin(trx(larva).tailcentralang(1:end-1));
    %velcentraltc{1,i}=trx(larva).velmagcentral.*(cos(trx(larva).velangcentral).*cos(trx(larva).tailcentralang(1,1:end-1))+sin(trx(larva).velangcentral).*sin(trx(larva).tailcentralang(1,1:end-1)));
end

units=parseunits('mm/s');
data=velcentraltc;

author: kristinbranson

repo: JAABA


text:
function [newTree, pdeSign] = splitTreePDE(treeIn)
%SPLITTREEPDE   Split a syntax tree, specific to PDE problems.
%   [NEWTREE, PDESIGN] = SPLITTREPDE(TREEIN) goes through the syntax tree
%   TREEIN, and isolates the part where the PDE variable, e.g. u_t, appears.
%   splits the syntax tree TREEIN into two trees, NEWTREE and LAMBDATREE.
%   LAMBDATREE contains the syntax tree that the eigenvalue parameter LAMBDA
%   appears in, NEWTREE contains the other part of TREEIN. The value of
%   LAMBDASIGN corresponds to the sign in front of LAMBDA in the original syntax
%   tree TREEIN.
%
% See also: STRINGPARSER/SPLITTREE.

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

% Begin by replacing the subtree which contains the pde_variable with a 0
[newTree, ignored, pdeSign] = findPDE(treeIn, 1);

% Do the basic splitting (converting = into -) in newTree
newTree = stringParser.splitTree(newTree);

end

function [newTree, pdeTree, pdeSign] = findPDE(treeIn, pdeSign)
%FINDPDE    Find where the PDE variable (e.g. u_t) appears in TREEIN.

% Initialization.
newTree = treeIn;
% Create some empty trees.
leftEmpty = 1;
rightEmpty = 1;
pdeTree = [];
pdeTreeLeft = [];
pdeTreeRight = [];
treeCenter = treeIn.center;

% Go recursively through the tree on the left.
if ( isfield(treeIn, 'left') )
    [newLeft, pdeTreeLeft, pdeSign] = findPDE(treeIn.left, pdeSign);
    newTree.left = newLeft;
    leftEmpty = 0;
end

% Go recursively through the tree on the right.
if ( isfield(treeIn, 'right') )
    [newRight, pdeTreeRight, pdeSign] = findPDE(treeIn.right, pdeSign);
    newTree.right = newRight;
    rightEmpty = 0;
end

% Return a new pdeTree. If the operator in the center of treeIn is a *,
% we want to return the whole treeIn (e.g. when we see 1*u_t). If not,
% we return the latest pdeTree (e.g. when we see u_t+1).
%
% Start by looking through the left tree.
if ( ~isempty(pdeTreeLeft) )
    if ( strcmp(treeCenter{2}, 'OP*') )
        pdeTree = treeIn;
    else
        pdeTree = pdeTreeLeft;
    end
end

% Now look through the right tree.
if ( ~isempty(pdeTreeRight) )
    if ( strcmp(treeCenter, 'OP*') )
        pdeTree = treeIn;
    % If we have a =, and u_t is on the right, we need to switch signs on
    % the pdeTree.
    elseif ( strcmp(treeCenter{2}, 'OP=') )
        disp('PDE on right')
        pdeSign = -1*pdeSign;
        pdeTree = pdeTreeRight;
    % If we have a -, and we have u_t on the right we need to switch signs
    % on the pdeTree.
    elseif ( strcmp(treeCenter{2}, 'OP-') || strcmp(treeCenter{2}, 'UN-') )
        pdeSign = -1*pdeSign;
        pdeTree = pdeTreeRight;
    else
        pdeTree = pdeTreeRight;
    end
end

% Both left and right trees are empty. We must be at a leaf! Replace the PDEVAR
% with a 0.
if ( leftEmpty && rightEmpty )
    % We encounter a PDE variable. Replace it by a zero.
    if ( strcmp(treeCenter{2}, 'PDEVAR') )
        pdeTree = newTree;
        newTree = struct('center', {{'0', 'NUM'}});
    end
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function so = scaleStream(si,factor)
% SCALESTREAM  Scale a stream (of numerics)
so = {head(si)*factor,delayEval(@scaleStream,{tail(si),factor})};

author: Sable

repo: mcbench-benchmarks


text:
function dx = pend_dx(u, x)

M = 1;   % kg
m = 0.3; % kg
g = 9.81;% m s^-2
L = 0.6; % m

a = x(2);
da = x(4);

if a == 0
	sinca = 1;
else
	sinca = sin(a)/a;
end

E = [...
	eye(2,2) zeros(2,2);
	zeros(2,2) [M/m+1, L*cos(a); cos(a), 4/3*L]];
A = [...
	zeros(2,2) eye(2,2);
	[0 0; 0 g*sinca] [0 L*da*sin(a); 0 0]];

B = [0 0 1/m 0]';

% E dx = A x + B u
piE = pinv(E);
dx = piE*A*x + piE*B*u;


text:
function xyzankurAnim(X, fps)

% XYZANKURANIM Animate point cloud of stick man from Agarwal & Triggs dataset.
% FORMAT
% DESC animates a matrix of x,y,z point clound positions representing the
% motion of the figure used to generate the silhouttes for Agarwal &
% Triggs silhouette data.
% ARG y : the data to animate.
% ARG fps : the number of frames per second to animate (defaults to 24).
%
% SEEALSO : xyzankurVisualise, xyzankurModify
%
% COPYRIGHT : Carl Henrik Ek and Neil Lawrence, 2008

% MOCAP

if(nargin<3)
  fps = 24;
  if(nargin<2)
    fid = 1;
    if(nargin<1)
      error('Too few arguments');
    end
  end
end

for(i = 1:1:size(X,1))
  if(i==1)
    handle = xyzankurVisualise(X(i,:),1);
  else
    xyzankurModify(handle,X(i,:));
  end
  pause(1/fps);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function pfplot(X,Factors,Weights,Option);
%PFPLOT plot parafac model
%
% See also:
% 'parafac'
%
%
% pfplot(X,Factors,Weights,Option);
% Different aspects for evaluation of the solution.
%
% Option # = 1
% 1	NOT ACCESIBLE
% 2	NOT ACCESIBLE
% 3	DIAGONALITY PLOT
% 4	PLOTS OF RESIDUAL VARIANCE
% 5	PLOTS OF LEVERAGE
% 6	RESIDUALS (STANDARD DEVIATION) VERSUS LEVERAGE
% 7	NORMAL PROBABILITY PLOT
% 8	LOADING PLOT
%
% You HAVE to input all four inputs. If you have no weights, just input [].
% The last input must be an 8-vector with ones if you want the plot and
% zeros else. E.g.
%
% pfplot(X,factors,[],[0 0 1 0 0 0 0 1]);
%
% to have the diagonality and the loading plot
%

% $ Version 1.02 $ Date 28. July 1998 $ Not compiled $
% $ Version 1.03 $ Date 6. October 1999 $ Changed to handle missing values correctly$
% $ Version 2.00 $ May 2001 $ Changed to array notation $ RB $ Not compiled $

% Copyright (C) 1995-2006  Rasmus Bro & Claus Andersson
% Copenhagen University, DK-1958 Frederiksberg, Denmark, rb@life.ku.dk
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
% Street, Fifth Floor, Boston, MA  02110-1301, USA.

warning off

DimX = size(X);
X = reshape(X,DimX(1),prod(DimX(2:end)));

% Convert to old format
NewLoad = Factors;
ff = [];
for f=1:length(Factors)
  ff=[ff;Factors{f}(:)];
end
Factors = ff;

factors = Factors;
ord=length(DimX);
Fac=length(factors)/sum(DimX);
lidx(1,:)=[1 DimX(1)*Fac];
for i=2:ord
  lidx=[lidx;[lidx(i-1,2)+1 sum(DimX(1:i))*Fac]];
end
if Option(3)==1
  % ESTIMATE DIAGONALITY OF T3-CORE
  diagonality=corcond(reshape(X,DimX),NewLoad,Weights,1);
end
model=nmodel(NewLoad);
model = reshape(model,DimX(1),prod(DimX(2:end)));
if Option(4)==1
  % PLOTS OF RESIDUAL VARIANCE
  figure,eval(['set(gcf,''Name'',''Residual variance'');']);
  aa=ceil(sqrt(ord));bb=ceil(ord/aa);
  for i=1:ord
    r=nshape(reshape(X-model,DimX),i)';
    varian=stdnan(r).^2;
    subplot(aa,bb,i)
    plot(varian)
    if DimX(i)<30
      hold on
      plot(varian,'r+')
    end
    eval(['xlabel(''Mode ', num2str(i),''');']);
    ylabel('Residual variance');
  end
end
if Option(5)==1
  % PLOTS OF LEVERAGE
  figure
  eval(['set(gcf,''Name'',''Leverage'');']);
  aa=ceil(sqrt(ord));
  bb=ceil(ord/aa);
  for i=1:ord
    A=reshape(factors(lidx(i,1):lidx(i,2)),DimX(i),Fac);
    lev=diag(A*pinv(A'*A)*A');
    subplot(aa,bb,i)
    if std(lev)>eps
      plot(lev+100*eps,'+')
      for j=1:DimX(i)
        text(j,lev(j),num2str(j))
      end
    else
      warning('Leverage is constant')
    end
    eval(['xlabel(''Mode ', num2str(i),''');']);
    ylabel('Leverage');
  end
end
if Option(6)==1
  % RESIDUALS (STANDARD DEVIATION) VERSUS LEVERAGE
  figure
  eval(['set(gcf,''Name'',''Residuals vs. Leverages'');']);
  aa=ceil(sqrt(ord));
  bb=ceil(ord/aa);
  for i=1:ord
    subplot(aa,bb,i)
    A=reshape(factors(lidx(i,1):lidx(i,2)),DimX(i),Fac);
    lev=diag(A*pinv(A'*A)*A')'+100*eps;
    r=nshape(reshape(X-model,DimX),i)';
    stand=stdnan(r);
    if std(lev)>eps
      plot(lev,stand,'+')
      for j=1:DimX(i)
        text(lev(j),stand(j),num2str(j))
      end
      eval(['xlabel(''Leverage in mode ', num2str(i),''');']);
      ylabel('Standard deviation');
    else
      warning('Leverage is constant')
    end
  end
end
if Option(7)==1
  % NORMAL PROBABILITY PLOT
  if exist('normplot')
    disp(' ')
    disp(' Normal probability plots are time-consuming')
    disp(' They are made in the statistics toolbox though, so we can''t change that!')
    figure,
    eval(['set(gcf,''Name'',''Normal probability of residuals'');']);
    aa=ceil(sqrt(ord));
    bb=ceil(ord/aa);
    r=nshape(reshape(X-model,DimX),i)';
    r=r(:);
    normplot(r(find(~isnan(r))))
  end
end
if Option(8)==1
  % LOADING PLOT
  if sum(Option)>1
    figure
  end
  eval(['set(gcf,''Name'',''Loadings'');']);
  aa=ceil(sqrt(ord));
  bb=ceil(ord/aa);
  for i=1:ord
    subplot(aa,bb,i)
    A=reshape(factors(lidx(i,1):lidx(i,2)),DimX(i),Fac);
    plot(A)
    eval(['xlabel(''Mode ', num2str(i),''');']);
    ylabel('Loading');
  end
end
drawnow

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function value = f1sd1 ( x )

%*****************************************************************************80
%
%% F1SD1 evaluates the function 1.0D+00/ sqrt ( 1.1 - x**2 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the function.
%
%    Output, real VALUE, the value of the function.
%
  f1sd1 = 1.0 / sqrt ( 1.1 - x * x );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = m_sxbfbt(x_it,aparams,mparams)
%m_sxbfbt   ideal stabilized least squares commutator preconditioner
%   y = m_sxbfbt(x_it,aparams,mparams);
%   input
%          x_it         operand for preconditioning operator
%          aparams      structure defining coefficient matrix
%          mparams      structure defining preconditioning matrix
%   output
%          y            result of preconditioning operation
%
%   IFISS function: HCE; 15 March 2005.

nv = length(aparams.F);
nu = nv/2;
np = size(aparams.B,1);
Gdiag=spdiags(diag(mparams.G),0,nv,nv);

sigma = mparams.viscosity;

rv=x_it(1:nv); rp=x_it(nv+1:nv+np);

%% pressure solve
xB = (Gdiag\aparams.B')';
BBt = aparams.B*xB';

if mparams.domain==1,
   n_null = mparams.n_null;
   minor = [1:n_null-1,n_null+1:np]';
   rp1 = zeros(np,1);
   rp1(minor) = (BBt(minor,minor)+mparams.Cp1(minor,minor))\rp(minor);
   rp2 =  xB*(aparams.F*(xB'*rp1)) + sigma*(mparams.Cp2*rp1);
   zp = zeros(np,1);
   zp(minor) = - (BBt(minor,minor)+mparams.Cp1(minor,minor))\rp2(minor);
else
   zp = (BBt+mparams.Cp1)\rp;
   zp = -(BBt+mparams.Cp1) \ (xB*(aparams.F*(xB'*zp)) + sigma*(mparams.Cp2*zp) );
end

%% velocity solve
rv = rv-(aparams.B')*zp;
zv = aparams.F \ rv;
y = [zv;zp];


text:
addpath(genpath('../mutils/My/'));
addpath(genpath('../ptv'));
basepth = '../../../data_prj/dir_dataset/DIR_files/';

TIME_e = zeros(2,2,2, 10);
TRE_e = zeros(2,2,2, 10);

for use_refinement = 0 : 1
for resize = 0 : 1
for fast_lcc = 0 : 1
if use_refinement && ~resize
    continue;
end

for idx = 1:10
    pts_struct = DIR_get_all_points_for_the_case(idx, basepth);
    [volmov, spc] = read_DIR_volume_4dCT(idx, 5, basepth);
    volmov = double(volmov);
    pts_mov = pts_struct.extreme.e;
    [volfix, spc] = read_DIR_volume_4dCT(idx, 0, basepth);
    volfix = double(volfix);
    pts_fix = pts_struct.extreme.b;

    volmov = img_thr(volmov, 80, 900, 1);
    volfix = img_thr(volfix, 80, 900, 1);

    % crop images
    init_size = size(volmov);
    min_max1 = [ min(pts_mov, [], 1)', max(pts_mov, [], 1)'];
    min_max2 = [ min(pts_fix, [], 1)', max(pts_fix, [], 1)'];
    min_max = [ min(min_max1(:, 1), min_max2(:, 1)), max(min_max1(:, 2), min_max2(:, 2))];
    d = [10, 10, 5];
    crop_v = [ max(1, min_max(1,1) - d(1)), min(size(volmov, 1), min_max(1,2) + d(1)); ...
               max(1, min_max(2,1) - d(2)), min(size(volmov, 2), min_max(2,2) + d(2)); ...
               max(1, min_max(3,1) - d(3)), min(size(volmov, 3), min_max(3,2) + d(3));];
    volmov = crop_data(volmov, crop_v);
    volfix = crop_data(volfix, crop_v);

    spc_orig = spc;
    if resize
        bszv = size(volmov);
        spc_tmp = [1, 1, 1];
        volfix = volresize(volfix, round(bszv .* spc .* spc_tmp), 1);
        volmov = volresize(volmov, round(bszv .* spc .* spc_tmp), 1);
        spc = [1,1,1] ./ spc_tmp;
    end

    % configure registration
    opts = [];
    opts.loc_cc_approximate = fast_lcc;
    if use_refinement
        if resize
            opts.grid_spacing = [4, 4, 4]*2;
        else
            opts.grid_spacing = [4, 4, 3]*2;  % grid spacing in pixels
        end
        opts.cp_refinements = 1;
    else
        if resize
            opts.grid_spacing = [4, 4, 4];
        else
            opts.grid_spacing = [4, 4, 3];  % grid spacing in pixels
        end
        opts.cp_refinements = 0;
    end
    opts.display = 'off';
    opts.k_down = 0.7;
    opts.interp_type = 0;
    opts.metric = 'loc_cc_fftn_gpu';
    opts.metric_param = [1,1,1] * 2.1;
    opts.scale_metric_param = true;
    opts.isoTV = 0.11;
    opts.csqrt = 5e-3;
    opts.spline_order = 1;
    opts.border_mask = 5;
    opts.max_iters =  80;
    opts.check_gradients = 100*0;
    opts.pix_resolution = spc;

    timer = tic;
    [voldef, Tptv, Kptv] = ptv_register(volmov, volfix, opts);
    TIME_e(use_refinement+1, resize+1, fast_lcc+1, idx) = toc(timer);

    if resize
        Tptv_rsz = cat(4, volresize(Tptv(:,:,:,1), bszv), volresize(Tptv(:,:,:,2), bszv), volresize(Tptv(:,:,:,3), bszv));
        voldef = volresize(voldef, bszv);
    else
        Tptv_rsz = Tptv;
    end
    [~, Tptv_rsz] = uncrop_data(voldef, Tptv_rsz, crop_v, init_size);
    % move points and measure TRE
    [pt_errs_phys, pts_moved_pix, TRE_phys, TREstd_phys] = DIR_movepoints(pts_mov, pts_fix, Tptv_rsz, spc_orig, []);
    TREs(idx) = mean(TRE_phys);

    TRE_e(use_refinement+1, resize+1, fast_lcc+1, idx) = mean(TRE_phys);
    save('results/DIR_TRE.mat', 'TRE_e', 'TIME_e');
end
end
end
end

repo: pTVreg


text:
function [ ohyps ] = rotateHyps( hyps, R )
%ROTATEHYPS Rotate room layout hypothesis
%
numHyps = length(hyps);
ohyps = hyps; %repmat(hyps(1), [numHyps 1]);
for hid = 1:numHyps
    ohyps(hid).extLine = rotateLines(hyps(hid).extLine, R);
    ohyps(hid).srcLine = rotateLines(hyps(hid).srcLine, R);
    ohyps(hid).hCorner = rotatePoint(hyps(hid).hCorner, R);
end

end

author: zouchuhang

repo: LayoutNet


text:
function FindMother(j)
global uLINK

if j ~= 0
    if j == 1
        uLINK(j).mother = 0;
    end
    if uLINK(j).child ~= 0
        uLINK(uLINK(j).child).mother = j;
        FindMother(uLINK(j).child);
    end
    if uLINK(j).sister ~= 0
        uLINK(uLINK(j).sister).mother = uLINK(j).mother;
        FindMother(uLINK(j).sister);
    end
end

author: s-kajita

repo: IntroductionToHumanoidRobotics


text:
% GSPBOX - Filters
%
%  Design - N filter
%    gsp_design_mexican_hat       -  Design a mexican hat filterbank
%    gsp_design_abspline          -  Design a abspline filterbank
%    gsp_design_meyer             -  Design a Meyer filterbank (tight)
%    gsp_design_simple_tf         -  Design a simple tight frame filterbank
%    gsp_design_itersine          -  Design a itersine filterbank (tight)
%    gsp_design_half_cosine       -  Design a half cosine filterbank (tight)
%    gsp_design_warped_translates -  Design a filterbank with a warping function
%
%  Design - 2 filter (LP - HP) tight filterbank
%    gsp_design_regular           -  Design 2 filter with the "regular" construction
%    gsp_design_held              -  Design 2 filter with the "Held" construction
%    gsp_design_simoncelli        -  Design 2 filter with the "Simoncelli" construction
%    gsp_design_papadakis         -  Design 2 filter with the "Papadakis" construction
%
%  Dual filterbank
%    gsp_design_can_dual          -  Design the canonical dual filterbank
%    gsp_evaluate_can_dual        -  Evaluate the canonical dual of a filterbank
%    gsp_test_duality             -  Test if 2 filterbanks are dual
%    gsp_test_duality_coefficient -  Test if 2 discrete filterbanks are dual
%
%  Low pass filters
%    gsp_design_heat              -  Design a heat kernel filter
%    gsp_design_expwin            -  Design a expwin filter
%    gsp_design_smooth_indicator  -  Design a smooth indicator function
%
%  Application
%    gsp_filter                   -  Shortcut to gsp_filter_analysis
%    gsp_filter_analysis          -  Analysis operator for filterbank
%    gsp_filter_synthesis         -  Synthesis operator for filterbank
%    gsp_filter_inverse           -  Inverse operator for filterbank
%
%  Joint Time-Vertex Filter Design
%    gsp_jtv_design_diffusion         -  Design a diffusion filterbank
%    gsp_jtv_design_wave              -  Design a wave filterbank
%    gsp_jtv_design_damped_wave       -  Design a damped wave filterbank
%    gsp_jtv_design_dgw               -  Design a generic dynamic graph wavelet
%
%  Joint Time-Vertex Filter Application
%    gsp_jtv_filter_analysis      -  Analysis operator for time-vertex filterbank
%    gsp_jtv_filter_synthesis     -  Synthesis operator for time-vertex filterbank
%    gsp_jtv_filter_evaluate      -  Evaluate a time-vertex filterbank
%    gsp_jtv_filter_array         -  Convert a ts/js time-vertex filter to a ts/js-array time-vertex filterbank
%    gsp_jtv_compute_frame        -  Return the matrix operator associated to a time-vertex filterbank
%    gsp_jtv_evaluate_can_dual    -  Evaluate the canonical dual of a time-vertex filterbank
%    gsp_jtv_design_can_dual      -  Design the canonical dual of a time-vertex filterbank
%    gsp_filter_inverse           -  Inverse operator for a time-vertex filterbank
%
%  Size Handling
%    gsp_mat2vec                  -  Matrix to vector representation for filterbanks
%    gsp_vec2mat                  -  Vector to matrix representation for filterbanks
%
%  Utils
%    gsp_approx_filter            -  Create an approximation of a filterbank with Chebyshev
%    gsp_wlog_scales              -  Compute log scale vector for wavelets
%    gsp_filter_evaluate          -  Evaluate a filterbank
%    gsp_filterbank_bounds        -  Bound for the filterbank
%    gsp_tighten_filter           -  Create a filter that tighten the filterbank
%    gsp_warp_filter              -  Warp a filter
%    gsp_multiply_filters         -  Multiply two filters
%    gsp_filterbank_matrix        -  Return the matrix operator associated to a filterbank
%
%  For help, bug reports, suggestions etc. please send email to
%  gspbox 'dash' support 'at' groupes 'dot' epfl 'dot' ch
%

% To be done
%   - add function in doc from meyer and simple_tf

repo: gspbox


text:
function  [image_data,vlat,vlon] = mygrid_sand2(region,ssfname)
%  MYGRID_SAND2  Read bathymetry data from Sandwell Database
%    [Z,LAT,LON] = MYGRID_SAND2(REGION) extracts data from
%    the Sandwell and Smith bathymetry, which is now at 1-minute
%    resolution.
%
%
% WARNING: change ssfname and ssversion to the correct one for
%          your machine
%
%                                               Catherine de Groot-Hedlin
%                                               modified Rich Pawlowicz
%
% latitudes must be between -80.738 and 80.738;
%       input:
%               REGION =[west east south north];
%       output:
%               Z - matrix of sandwell bathymetry/topography
%               LAT - vector of latitudes associated with image_data
%               LON - vector of longitudes
%
%  An odd Z value of say -2001m signifies that this pixel was constrained
%  by a real depth sounding while an even depth of say -2000m is
%  a predicted depth.  Plot rem(Z,2) to see this.
%
%  Changes  Jan/18 - modified to work with post v9 database versions
%                    also order of arguments, and various other
%                    things that bothered me.

if nargin<1
    region=[-140 -121 46 52];
end

if nargin<2
  %ssfname='/ocean/rich/more/mmapbase/s_and_s/topo_8.2.img';
  %ssversion=8.2;
  ssfname='/ocean/rich/more/mmapbase/s_and_s/topo_18.1.img';
  ssversion=18.1;
end

%------no changes below here-------------

% determine the requested region
wlon = region(1);
elon = region(2);
blat = region(3);
tlat = region(4);

% Setup the parameters for reading Sandwell data version
if ssversion<9
   db_res         = 2/60;          % 2 minute resolution
   db_loc         = [-72.006 72.006 0.0 360-db_res];
   db_size        = [6336 10800];
else
   db_res         = 1/60;          % 1 minute resolution
   db_loc         = [-80.738  80.738 0.0 360-db_res];
   db_size        = [17280 21600];
end
nbytes_per_lat = db_size(2)*2;  % 2-byte integers

% Check ranges

if blat<db_loc(1)
    blat=db_loc(1);
end
if tlat>db_loc(2)
    tlat=db_loc(2);
end

% Determine if the database needs to be read twice (overlapping prime meridian)
if ((wlon<0)&&(elon>=0))
      wlon      = [wlon           0];
      elon      = [360-db_res  elon];
end

% Calculate number of "records" down to start (latitude) (0 to db_size(1)-1)
% (mercator projection)
rad=pi/180;

arg1=log(tand(45+db_loc(1)/2));
arg2=log(tand(45+blat/2));
iblat = fix(db_size(1) +1 - (arg2-arg1)/(db_res*rad));

arg2=log(tand(45+tlat/2));
itlat = fix(db_size(1) +1 - (arg2-arg1)/(db_res*rad));

if (iblat < 0 ) || (itlat > db_size(1)-1)
        errordlg([' Requested latitude is out of file coverage ']);
end

% Go ahead and read the database
% Open the data file
fid = fopen( ssfname, 'r','b');
if (fid < 0)
        error(['Could not open database: '  ssfname],'Error');
end

image_data     = [];
for i = 1:length(wlon)

        % Make sure the longitude data goes from 0 to 360
        if wlon(i) < 0
                wlon(i) = 360 + wlon(i);
        end

        if elon(i) < 0
                elon(i) = 360 + elon(i);
        end

        % Calculate the longitude indices into the matrix (0 to db_size(1)-1)
        iwlon(i) = fix((wlon(i)-db_loc(3))/db_res);
        ielon(i) = fix((elon(i)-db_loc(3))/db_res);
        if (iwlon(i) < 0 ) || (ielon(i) > db_size(2)-1)
               error([' Requested longitude is out of file coverage ']);
        end

        % allocate memory for the data
        data = zeros(iblat-itlat+1,ielon(i)-iwlon(i)+1);

        % Skip into the appropriate spot in the file, and read in the data
        disp('Reading in bathymetry data');
        for ilat = itlat:iblat
                offset = ilat*nbytes_per_lat + iwlon(i)*2;
                status = fseek(fid, offset, 'bof');
                data(iblat-ilat+1,:)=fread(fid,[1,ielon(i)-iwlon(i)+1],'integer*2');
        end

        % put the two files together if necessary
        if (i>1)
                image_data = [image_data data];
        else
                image_data = data;
        end
end
% close the file
fclose(fid);

% Determine the coordinates of the image_data
vlat=zeros(1,iblat-itlat+1);
arg2 = log(tand(45+db_loc(1)/2.));
for ilat=itlat+1:iblat+1
        arg1 = rad*db_res*(db_size(1)-ilat+0.5);
        term=exp(arg1+arg2);
        vlat(iblat-ilat+2)=2*atand(term)-90;
end
vlon=db_res*((iwlon+1:ielon+1)-0.5);

% to plot it up

if nargout==0
  [xx,yy]=meshgrid(vlon,vlat);
  pcolor(xx,yy,image_data),shading flat,colormap(jet),colorbar('vert')
  xlabel('longitude'),ylabel('latitude'),title('Smith and Sandwell bathymetry Example')
  clear image_data vlon vlat
end

author: goGPS-Project


text:
function y = round(varargin)
% Implements Matlab's `round` with 2 or 3 inputs
% For symbolic data: only 1 input is supported by the original round function, so only that case is
% covered here
if nargin==1
    if ~isa(varargin{1},'sym')
        y = builtin('round',varargin{1});
    else
        y = builtin('@sym/round',varargin{1});
    end
elseif nargin==2 || (nargin==3 && strcmp(varargin{3}, 'decimal'))
    e = 10^varargin{2};
    y = builtin('round',varargin{1}*e)/e;
elseif nargin==3 && strcmp(varargin{3}, 'significant')
    e = 10^(varargin{2} - (floor(log10(abs(varargin{1})))+1));
    y = builtin('round',varargin{1}*e)/e;
else
    error('Incorrect inputs')
end
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
classdef dmce_shunt_mpc2 < mp.dmc_element % & mp.dmce_shunt
%MP.DMCE_SHUNT_MPC2  Data model converter for shunt elements for MATPOWER case v2.

%   MATPOWER
%   Copyright (c) 2021-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
        bus
    end     %% properties

    methods
        function name = name(obj)
            name = 'shunt';
        end

        function df = data_field(obj)
            df = 'bus';
        end

        function [nr, nc, r] = get_import_size(obj, mpc)
            %% define named indices into data matrices
            [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
               VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;

            if isfield(mpc, obj.data_field())
                tab = mpc.(obj.data_field());
            else
                tab = [];
            end
            if isempty(tab)
                nr = 0;
                nc = 0;
                r = [];
            else
                r = find(tab(:, GS) | tab(:, BS));
                obj.bus = r;
                nr = size(r, 1);
                nc = size(tab, 2);          %% use nc of default table
            end
        end

        function [nr, nc, r] = get_export_size(obj, dme)
            [nr, nc] = size(dme.tab);   %% use size of default table
            r = dme.tab.source_uid;     %% rows in bus matrix
        end

        function vmap = table_var_map(obj, dme, mpc)
            vmap = table_var_map@mp.dmc_element(obj, dme, mpc);

            %% define named indices into data matrices
            [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
               VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;

            %% mapping for each name, default is {'col', []}
            vmap.uid        = {'IDs'};      %% consecutive IDs, starting at 1
            vmap.name       = {'cell', ''}; %% empty char
            vmap.status     = {'num', 1};   %% ones
            vmap.source_uid = {'r'};        %% row index in mpc.bus
            vmap.bus{2}     = BUS_I;
            vmap.gs{2}      = GS;
            vmap.bs{2}      = BS;
            vmap.p          = {'num', 0};   %% zeros
            vmap.q          = {'num', 0};   %% zeros
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function str = hyperlink(url, label, msg)
%HYPERLINK  create a string that is displayable as hyperlink in Matlab console
%
% Usage examples:
%   fprintf('Search on %s\n', hyperlink('http://google.com','Google'));
%   fprintf(hyperlink('http://google.com','Google','Search on Google\n'));
%
% HYPERLINK converts the specified URL and text label into a string that is
% displayed as a hyperlink in the Matlab console (Command Window).
% In a deployed (compiled)  program, the URL and text label (if different
% from the URL) are displayed in a non-hyperlinked plain-text manner.
% If the optional 3rd input argument (msg) is specified, then all instances of
% the specified label within msg will be handled as above (hyperlinked etc.)
%
% IN:
%   url   - (mandatory) URL of webpage or Matlab command (e.g., 'matlab:which(...')
%   label - (optional)  text label of the hyperlink. Default: url
%   msg   - (optional)  string in which all label instances should be hyperlinked
%
% OUT:
%   str - string output

% Copyright: Yair Altman 2020-
%{
% 15/01/20 - Initial version
%}

    error(nargchk(1,3,nargin)); %#ok<NCHKN> narginchk is only available in R2011b+
    if nargin > 2 % msg was specified
        % replace all instances of label within msg with corresponding hyperlink
        str = strrep(msg, label, hyperlink(url,label));
        return
    end
    if nargin < 2, label = url; end
    isWebpage = strncmpi(url,'http',4);

    % Only hyperlink in interactive Matlab sessions, not in a deployed program
    if ~isdeployed  % interactive Matlab session
        if isWebpage  % open in a web browser
            str = ['' label ''];
        else  % Matlab command e.g. 'matlab:which(...'
            str = ['' label ''];
        end
    else  % deployed (compiled)
        if isWebpage && ~strcmp(label,url)  % display label next to url
            str = [label ' (' url ')'];
        elseif isWebpage  % text==url - display only the url
            str = url;
        else  % Matlab command (not a webpage) - just display the label
            str = label;
        end
    end
end

author: altmany


text:
filename='Cantileverbeam_Tetrahedra_Linear_Structured';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'holes';
cost = {'compliance'};
weights = [1];
constraint = {'volumeConstraint'};
optimizer = 'SLERP';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';
constraint_case = 'INEQUALITY';

nsteps = 10;
Vfrac_final = 0.2;
Perimeter_target=3.5;
optimality_final =1e-5;
constr_final =1e-5;

BCscale_factor = 0.3;
HJiter0 = 1;
e2 = 1;
N_holes = [6 3 3];
R_holes = 0.1;
phase_holes = [0 0 0];

Vfrac_initial = 0.8;
optimality_initial = 5e-2;
constr_initial = 5e-2;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;

% maxiter = 10;
% maxiter = 1;
TOL.nu_plus = 0.3;
TOL.nu_minus = 0.3;

plotting = 1;
printing = 1;
monitoring = 1;
monitoring_interval = 1;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


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

%%
clear; close all; clc;

%% Syntax
% |S=nowStr(f);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [xvals, yvals, color] = hintmat(w);
%HINTMAT Evaluates the coordinates of the patches for a Hinton diagram.
%
%	Description
%	[xvals, yvals, color] = hintmat(w)
%	  takes a matrix W and returns coordinates XVALS, YVALS for the
%	patches comrising the Hinton diagram, together with a vector COLOR
%	labelling the color (black or white) of the corresponding elements
%	according to their sign.
%
%	See also
%	HINTON
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Set scale to be up to 0.9 of maximum absolute weight value, where scale
% defined so that area of box proportional to weight value.

w = flipud(w);
[nrows, ncols] = size(w);

scale = 0.45*sqrt(abs(w)/max(max(abs(w))));
scale = scale(:);
color = 0.5*(sign(w(:)) + 3);

delx = 1;
dely = 1;
[X, Y] = meshgrid(0.5*delx:delx:(ncols-0.5*delx), 0.5*dely:dely:(nrows-0.5*dely));

% Now convert from matrix format to column vector format, and then duplicate
% columns with appropriate offsets determined by normalized weight magnitudes.

xtemp = X(:);
ytemp = Y(:);

xvals = [xtemp-delx*scale, xtemp+delx*scale, ...
         xtemp+delx*scale, xtemp-delx*scale];
yvals = [ytemp-dely*scale, ytemp-dely*scale, ...
         ytemp+dely*scale, ytemp+dely*scale];

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
moviename = '/home/mayank/Work/FlySpaceTime/walkMovies/SS03500_test/movie.ufmf';
trackfilename = '/home/mayank/Work/FlySpaceTime/walkMovies/SS03500_test/trx.mat';
stationary = true;
method = 'deep-sup';
flowname = 'DS';

ff = fopen('DebugDeepCrash.txt','w');

tracks = load(trackfilename);
tracks = tracks.trx;

[readfcn,nframes,fid,headerinfo] = get_readframe_fcn(moviename);
fclose(fid);

blocksize = 50;

minfirst = min([tracks.firstframe]);
maxlast = max([tracks.endframe]);
nframes = maxlast-minfirst+1;
nblocks = ceil((nframes-1)/blocksize);

allftrs ={};

% compute features in parallel for different intervals of frames.
% parfor ndx = 1:nblocks
for ndx = 41:80
  [readfcn,nframes,fid,headerinfo] = get_readframe_fcn(moviename);
  fstart = minfirst + (ndx-1)*blocksize;
  fend = min(maxlast,ndx*blocksize);
  tic;
  allftrs{ndx} = genFeatures(readfcn,headerinfo,fstart,fend,tracks,stationary,method);
  telapsed = toc;
  fclose(fid);
  fprintf(ff,'%d\n',ndx);

end


text:
% Compare various inference engines on the following network (from Jensen (1996) p84 fig 4.17)
%    1
%  / | \
% 2  3  4
% |  |  |
% 5  6  7
%  \/ \/
%  8   9
% where all arcs point downwards
seed = 0;
rand('state', seed);
randn('state', seed);

N = 9;
dag = zeros(N,N);
dag(1,2)=1; dag(1,3)=1; dag(1,4)=1;
dag(2,5)=1; dag(3,6)=1; dag(4,7)=1;
dag(5,8)=1; dag(6,8)=1; dag(6,9)=1; dag(7,9) = 1;

dnodes = 1:N;
false = 1; true = 2;
ns = 2*ones(1,N); % binary nodes

onodes = [2 4];
bnet = mk_bnet(dag, ns, 'observed', onodes);
% use random params
for i=1:N
  bnet.CPD{i} = tabular_CPD(bnet, i);
end

%USEC = exist('@jtree_C_inf_engine/collect_evidence','file');
query = [3];
engine = {};
engine{end+1} = jtree_inf_engine(bnet);
engine{end+1} = jtree_sparse_inf_engine(bnet);
%engine{end+1} = jtree_ndx_inf_engine(bnet, 'ndx_type', 'SD');
%engine{end+1} = jtree_ndx_inf_engine(bnet, 'ndx_type', 'B');
%engine{end+1} = jtree_ndx_inf_engine(bnet, 'ndx_type', 'D');
%if USEC, engine{end+1} = jtree_C_inf_engine(bnet); end
%engine{end+1} = var_elim_inf_engine(bnet);
%engine{end+1} = enumerative_inf_engine(bnet);
%engine{end+1} = jtree_onepass_inf_engine(bnet, query, onodes);

maximize = 0;  % jtree_ndx crashes on max-prop
[err, time] = cmp_inference_static(bnet, engine, 'maximize', maximize);


text:
function out = IN_AutoMutualInfo(y,timeDelay,estMethod,extraParam)
% IN_AutoMutualInfo     Time-series automutual information
%
%---INPUTS:
%
% y: input time series (column vector)
%
% timeDelay: time lag for automutual information calculation
%
% estMethod: the estimation method used to compute the mutual information:
%           (*) 'gaussian'
%           (*) 'kernel'
%           (*) 'kraskov1'
%           (*) 'kraskov2'
%
% cf. Kraskov, A., Stoegbauer, H., Grassberger, P., Estimating mutual
% information: http://dx.doi.org/10.1103/PhysRevE.69.066138

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

% ------------------------------------------------------------------------------
% Check inputs:
% ------------------------------------------------------------------------------
if nargin < 2 || isempty(timeDelay)
    timeDelay = 1;
end
if ischar(timeDelay) && ismember(timeDelay,{'ac','tau'})
    timeDelay = CO_FirstCrossing(y,'ac',0,'discrete');
end

if nargin < 3 || isempty(estMethod)
    estMethod = 'kernel';
end

if nargin < 4
    extraParam = [];
end

N = length(y);
doPlot = false; % plot outputs to screen
minSamples = 5; % minimum 5 samples to compute a mutual information (could make higher?)

% Ensure y is a column vector
if size(y,2) > size(y,1)
    warning('Please input a column vector for y')
    y = y';
end

% ------------------------------------------------------------------------------
% Loop over time delays if a vector
numTimeDelays = length(timeDelay);
amis = nan(numTimeDelays,1);

if numTimeDelays > 1
    timeDelay = sort(timeDelay);
end

% Initialize miCalc object (needs to be reinitialized within the loop for kraskov):
if ~strcmp(estMethod,'gaussian')
    miCalc = IN_Initialize_MI(estMethod,extraParam,false); % NO ADDED NOISE!
end

for k = 1:numTimeDelays

    % Check enough samples to compute an automutual information
    if timeDelay(k) > N - minSamples
        % Time series is too short -- keep the remaining values as NaNs
        break
    end

    % Form the time-delay vectors y1 and y2
    y1 = y(1:end-timeDelay(k));
    y2 = y(1+timeDelay(k):end);

    if strcmp(estMethod,'gaussian')
        r = corr(y1,y2,'type','Pearson');
        amis(k) = -0.5*log(1 - r^2);
    else
        % Reinitialize for Kraskov:
        miCalc.initialise(1,1);

        % Set observations to time-delayed versions of the time series:
        miCalc.setObservations(y1,y2);

        % Compute:
        amis(k) = miCalc.computeAverageLocalOfObservations();
    end

    % Plot:
    if doPlot
        plot(y1,y2,'.k')
        title(sprintf('ami = %.3f',amis(k)))
        pause(0.1)
    end
end

if any(isnan(amis))
    warning(['Time series (N=%u) is too short for automutual information calculations',...
                ' up to lags of %u'],N,max(timeDelay))
end

if doPlot
    plot(amis,'-k')
end

%-------------------------------------------------------------------------------
% Outputs:
%-------------------------------------------------------------------------------
if numTimeDelays == 1
    out = amis; % a scalar
else
    % A structure
    for k = 1:numTimeDelays
        out.(sprintf('ami%u',timeDelay(k))) = amis(k);
    end
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [structVol, planC] = getStructureVol(structNum,planC)
% getStructureVol
% This function returns the abslute volume of a structure. Can be accessed
% from MATLAB command line or from CERR command.
%
% Created DK
% Usage
% structVol = getStructureVol(structNum)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

if ~exist('planC','var')
    global planC
end
indexS = planC{end};
optS    = planC{indexS.CERROptions};

%Get the scan number associated with the requested structure.
scanSet = getStructureAssociatedScan(structNum, planC);

deltaY = planC{indexS.scan}(scanSet).scanInfo(1).grid1Units;

ROIImageSize = [planC{indexS.scan}(scanSet).scanInfo(1).sizeOfDimension1  planC{indexS.scan}(scanSet).scanInfo(1).sizeOfDimension2];

%Get raster segments for structure.
[segmentsM, planC, isError] = getRasterSegments(structNum, planC);
numSegs = size(segmentsM,1);

%Relative sampling of ROI voxels in this place, compared to CT spacing.
%Set when rasterSegments are generated (usually on import).
sampleRate = optS.ROISampleRate;

%Sample the rows
indFullV =  1 : numSegs;
if sampleRate ~= 1
    rV = 1 : length(indFullV);
    rV([rem(rV+sampleRate-1,sampleRate)~=0]) = [];
    indFullV = rV;
end

%Block process to avoid swamping on large structures
if isfield(optS, 'DVHBlockSize') & ~isempty(optS.DVHBlockSize)
    DVHBlockSize = optS.DVHBlockSize;
else
    DVHBlockSize = 5000;
end

blocks = ceil(length(indFullV)/DVHBlockSize);
volsV  = [];
scansV = [];

start = 1;

for b = 1 : blocks

    %Build the interpolation points matrix

    dummy = zeros(1,DVHBlockSize * ROIImageSize(1));
    x1V = dummy;
    y1V = dummy;
    z1V = dummy;
    volsSectionV =  dummy;

    if start+DVHBlockSize > length(indFullV)
        stop = length(indFullV);
    else
        stop = start + DVHBlockSize - 1;
    end

    indV = indFullV(start:stop);

    mark = 1;

    for i = indV
        tmpV = segmentsM(i,1:10);
        delta = tmpV(5) * sampleRate;
        xV = tmpV(3): delta : tmpV(4);
        len = length(xV);
        rangeV = ones(1,len);
        yV = tmpV(2) * rangeV;
        zV = tmpV(1) * rangeV;
        sliceThickness = tmpV(10);
        %v = delta^2 * sliceThickness;
        v = delta * (deltaY*sampleRate) * sliceThickness;
        x1V(mark : mark + len - 1) = xV;
        y1V(mark : mark + len - 1) = yV;
        z1V(mark : mark + len - 1) = zV;
        volsSectionV(mark : mark + len - 1) = v;
        mark = mark + len;
    end

    %cut unused matrix elements
    x1V = x1V(1:mark-1);
    y1V = y1V(1:mark-1);
    z1V = z1V(1:mark-1);
    volsSectionV = volsSectionV(1:mark-1);

    %Get transformation matrices for both scan and structure.
    transMScan    = getTransM('scan', scanSet, planC);
    transMStruct  = getTransM('struct', structNum, planC);

    %Forward transform the structure's coordinates.
    if ~isempty(transMStruct)
        [x1V, y1V, z1V] = applyTransM(transMStruct, x1V, y1V, z1V);
    end

    %Back transform the coordinates into the scans' coordinate system.
    if ~isempty(transMScan)
        [x1V, y1V, z1V] = applyTransM(inv(transMScan), x1V, y1V, z1V);
    end

    %Interpolate.
    [scansSectionV] = getScanAt(scanSet, x1V, y1V, z1V, planC);

    scansV = [scansV, scansSectionV];
    volsV  = [volsV, volsSectionV];

    start = stop + 1;

end

volsV = volsV * sampleRate^2;  %must account for sampling rate!
structVol = sum(sum(volsV));

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [minM, maxM] = getMinMaxIM(IM,NHOOD)
% function [minM, maxM] = getMinMaxIM(IM,NHOOD)
%
% This function returns min and max for the local neighborhood NHOOD
%
% APA, 04/26/2012

% ---- No reflection required if our NHOOD is symmetric ------

% % NHOOD is reflected across its origin in order for IMDILATE
% % to return the local maxima of I in NHOOD if it is asymmetric. A symmetric NHOOD
% % is naturally unaffected by this reflection.
% reflectH = NHOOD(:);
% reflectH = flipud(reflectH);
% reflectH = reshape(reflectH, size(NHOOD));
% maxMat = imdilate(IM,reflectH);

maxM = imdilate(IM,NHOOD);

% IMERODE returns the local minima of IM in NHOOD.
minM = imerode(IM,NHOOD);

end


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

% problem 6 - convolution of x(t) and h(t)

t1=0:.1:2;
t2=2.1:.1:4;
t3=4.1:.1:10;
x1=t1;
x2=4-t2;
x3=zeros(size(t3));
x=[x1 x2 x3];
t=0:.1:10;
h=t.*exp(-t);
y=conv(x,h)*0.1;
plot(0:.1:20,y);
title('System response y(t)')


text:
function [s,iter]=IAAFT(x,c,maxiter)
%Syntax: [s,iter]=IAAFT(x,c,maxiter)
%___________________________________
%
% Makes c Iterative Amplitude Adjusted Fourier Transformed (AAFT) surrogates
% of a time series x.
%
% s is the IAAFT surrogate time series.
% iter is the number of iterations needed of the i-th surrogate series.
% x is the original time series.
% c is the number of surrogates.
% maxiter is the maximum number of iterations allowed.
%
%
% References:
%
% Schreiber T, Schmitz A (1996): Improved surrogate data for nonlinearity
% tests. Physical Review Letters 77: 635-638
%
% I am grateful to Dr D. Kugiumtzis for providing me a version of his Matlab
% function concerning this algorithm.
%
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110 - Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 20 November 2001.

if nargin<1 | isempty(x)==1
   error('You should provide a time series.');
else
   % x must be a vector
   if min(size(x))>1
      error('Invalid time series.');
   end
   x=x(:);
end

if nargin<2 | isempty(c)==1
   c=1;
else
   % c must be scalar
   if sum(size(c))>2
      error('c must be scalar.');
   end
   % c must be greater or equal than 1
   if c<1
      error('c must be greater or equal than 1.');
   end
end

if nargin<3 | isempty(maxiter)==1
   maxiter=1000;
else
   % maxiter must be scalar
   if sum(size(maxiter))>2
      error('maxiter must be scalar.');
   end
   % maxiter must be greater or equal than 1
   if maxiter<1
      error('maxiter must be greater or equal than 1.');
   end
end

% The magnitudes of x
amp=abs(fft(x));

% Shuffle x
s=shuffle(x,c);

% Sort x
[x,r]=sort(x);

for j=1:c

    % Calculate the phases of the shuffled series
    phase=angle(fft(s(:,j)));

    % Initialize the loop
    k=1;
    indold=r;
    converge = 0;
    while k<=maxiter & converge == 0
        % Make phase-randomized surrogates ...
        s(:,j)=amp.*exp(phase.*i);
        s(:,j)=real(ifft(s(:,j)));
        % ... and give them the distribution of x
        [s(:,j),T]=sort(s(:,j));
        [s(:,j),indnew]=sort(T);
        s(:,j)=x(indnew);
        % Check the convergence
        if indnew==indold
            converge=1;
        else
            indold=indnew;
            k=k+1;
        end
        % Loop again if needed, calculating the phases once more
        phase=angle(fft(s(:,j)));
    end

    % Get the iterations of each surrogate series
    iter(j)=k;

end


text:
function x = double( v ) %#ok
cvx_throw( 'Disciplined convex programming error:\n   Constraints may not appear in if/then statements.' );

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

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


text:
% THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION IS RELEASED "AS IS."  THE U.S. GOVERNMENT MAKES NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, CONCERNING THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION, INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT WILL THE U.S. GOVERNMENT BE LIABLE FOR ANY DAMAGES, INCLUDING ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE, OR INABILITY TO USE, THIS SOFTWARE OR ANY ACCOMPANYING DOCUMENTATION, EVEN IF INFORMED IN ADVANCE OF THE POSSIBILITY OF SUCH DAMAGES.
%
% file: get_ci_ml.m
% computes confidence interval
%  - assumes x is the number of successes in n Bernoulli trials
%  - finds confidence interval around the estimate for p
%  - with confidence level 1 - alpha
%  - method used is based on integrating f(p|x)
%  - finds an interval of minimal length
%
% 010312 tdr created from get_ci_nibp

function ci = get_ci_ml(x,n,alpha)

if nargin < 3, error('Requires three input arguments (x,n,alpha)'); end;

p_hat=x/n;

if ( p_hat == 0 | p_hat == 1)
    ci = interval_est_p(x,n,alpha,'cs1'); % min length is same a "symmetric" when it bumps up against zero or one.
    return;
end;

% the following assumes that 0 < p_hat < 1

y_too_small = 0;
y_too_big = binopdf_01ok(x,n,p_hat);
y_guess = (y_too_small + y_too_big)/2;
continue_search = 1;

%close_enough = alpha/1000;
close_enough = 0.00001;
convergence_limit = 100;
convergence_count = 0;
hit_limit = 0;

p_ll = p_hat/2;
p_ul = (1 + p_hat)/2;

% start search
while (continue_search)
    convergence_count = convergence_count + 1;

    % find p lower limit: p_ll
        continue_p_search = 1;
        p_convergence_count = 0;
        hit_p_limit = 0;
        p_close_enough = 0;
        p_limit_too_small = 0; % 0 for ll, p_hat for ul
        p_limit_too_big = p_hat; % p_hat for ll, 1 for ul
        p_limit = (p_limit_too_big + p_limit_too_small) / 2;
        while (continue_p_search)
            p_convergence_count = p_convergence_count + 1;
            p_limit_y_guess = binopdf_01ok(x,n,p_limit);

            if (p_limit_y_guess > y_guess)
                p_limit_too_big = p_limit; % too_big for ll, too_small for ul
            else
                p_limit_too_small = p_limit; % too_small for ll, too_big for ul
            end;
            % stop or adjust p_ll_guess
            if (convergence_limit == p_convergence_count) hit_p_limit = 1; end;
            if ((p_limit_too_big - p_limit_too_small) < (close_enough / (10*n))) p_close_enough = 1; end; %see DNp2454.
            if (p_close_enough | hit_p_limit)
                continue_p_search = 0;
            else
                p_limit = (p_limit_too_big + p_limit_too_small) / 2;
            end;
        end;
        p_ll = p_limit; % ll or ul
        if hit_p_limit warning('method ''ml'' convergence limit reached for lower interval value.'); end;

    % find p upper limit: p_ul
        continue_p_search = 1;
        p_convergence_count = 0;
        hit_p_limit = 0;
        p_close_enough = 0;
        p_limit_too_small = p_hat; % 0 for ll, p_hat for ul
        p_limit_too_big = 1; % p_hat for ll, 1 for ul
        p_limit = (p_limit_too_big + p_limit_too_small) / 2;
        while (continue_p_search)
            p_convergence_count = p_convergence_count + 1;
            p_limit_y_guess = binopdf_01ok(x,n,p_limit);

            if (p_limit_y_guess > y_guess)
                p_limit_too_small = p_limit; % too_big for ll, too_small for ul
            else
                p_limit_too_big = p_limit; % too_small for ll, too_big for ul
            end;
            % stop or adjust p_ul_guess
            if (convergence_limit == p_convergence_count) hit_p_limit = 1; end;
            if ((p_limit_too_big - p_limit_too_small) < (close_enough / (10*n))) p_close_enough = 1; end;
            if (p_close_enough | hit_p_limit)
                continue_p_search = 0;
            else
                p_limit = (p_limit_too_big + p_limit_too_small) / 2;
            end;
        end;
        p_ul = p_limit; % ll or ul
        if hit_p_limit warning('method ''ml'' convergence limit reached for upper interval value.'); end;

    alpha_guess = (n+1)*(integrate_binopdf(x,n,0,p_ll) + integrate_binopdf(x,n,p_ul,1));

    % adjust y_too_small, y_too_big
    if (alpha_guess > alpha)
        y_too_big = y_guess;
    else
        y_too_small = y_guess;
    end;

    % stop or compute new y_guess
    delta_alpha = abs(alpha - alpha_guess);
    if (convergence_count == convergence_limit) hit_limit = 1; end;
    if ((delta_alpha <= close_enough) | hit_limit)
        continue_search = 0;
    else
        y_guess = (y_too_big + y_too_small)/2;
    end;
end;

if hit_limit warning('method ''ml'' convergence limit reached - results may not be accurate.'); end;

ci=[p_hat p_ll p_ul];

return;


text:
function node_boundary = triangulation_order3_boundary_node ( node_num, ...
  triangle_num, triangle_node )

%*****************************************************************************80
%
%% TRIANGULATION_ORDER3_BOUNDARY_NODE indicates which nodes are on the boundary.
%
%  Discussion:
%
%    This routine is given a triangulation, an abstract list of triples
%    of nodes.  It is assumed that the nodes in each triangle are listed
%    in a counterclockwise order, although the routine should work
%    if the nodes are consistently listed in a clockwise order as well.
%
%    It is assumed that each edge of the triangulation is either
%    * an INTERIOR edge, which is listed twice, once with positive
%      orientation and once with negative orientation, or;
%    * a BOUNDARY edge, which will occur only once.
%
%    This routine should work even if the region has holes - as long
%    as the boundary of the hole comprises more than 3 edges!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer TRIANGLE_NUM, the number of triangles.
%
%    Input, integer TRIANGLE_NODE(3,TRIANGLE_NUM), the nodes that make up the
%    triangles.  These should be listed in counterclockwise order.
%
%    Output, logical NODE_BOUNDARY(NODE_NUM), is TRUE if the node
%    is on a boundary edge.
%
  m = 2;
  n = 3 * triangle_num;
%
%  Set up the edge array.
%
  edge(1:2,               1:  triangle_num) = triangle_node(1:2,1:triangle_num);
  edge(1:2,  triangle_num+1:2*triangle_num) = triangle_node(2:3,1:triangle_num);
  edge(1,  2*triangle_num+1:3*triangle_num) = triangle_node(3,  1:triangle_num);
  edge(2,  2*triangle_num+1:3*triangle_num) = triangle_node(1,  1:triangle_num);
%
%  In each column, force the smaller entry to appear first.
%
  e1(1:n) = min ( edge(1:2,1:n) );
  e2(1:n) = max ( edge(1:2,1:n) );

  edge(1,1:n) = e1(1:n);
  edge(2,1:n) = e2(1:n);
%
%  Ascending sort the column array.
%
  edge = ( sortrows ( edge' ) )';
%
%  Records which appear twice are internal edges and can be ignored.
%
  node_boundary(1:node_num) = 0;

  j = 0;

  while ( j < 3 * triangle_num )

    j = j + 1;

    if ( j == 3 * triangle_num )
      node_boundary(edge(1:m,j)) = 1;
    elseif ( all ( edge(1:m,j) == edge(1:m,j+1) ) )
      j = j + 1;
    else
      node_boundary(edge(1:m,j)) = 1;
    end

  end

  return
end


text:
% Find the number of frames in a tiff file
%
% This function reports number of frames in a tiff file by
% jumping over seekInterval frames until the end of the file.
% This is in general faster than Matlab's imfinfo function.
%
%  USAGE
%   n = nFrames(tiff, seekInterval)
%   tiff            Path to tiff file or a handler of an already opened file.
%   seekInterval    Optional number of frames to jump over.
%                   Default is 1000.
%   n               Number of frames
%
function n = nFrames(tiff, seekInterval)
    if nargin < 2
        seekInterval = 1000;
    end
    %keep guessing until we seek too far
    guess = seekInterval;
    overSeeked = false;
    n = 0;

    if ischar(tiff) || isstring(tiff)
        tiff = Tiff(tiff, 'r');
        closeTiff = onCleanup(@() close(tiff));
    end

    while ~overSeeked
      try
        tiff.setDirectory(guess);
        guess = 2*guess; %double the guess
      catch ex
        overSeeked = true; %we tried to seek past the last directory
      end
    end
    %when overseeking occurs, the current directory/frame will be the last one
    n = tiff.currentDirectory;
end

author: cortex-lab


text:
function lp = lagrange_derivative ( nd, xd, ni, xi )

%*****************************************************************************80
%
%% LAGRANGE_DERIVATIVE evaluates the derivative of the Lagrange basis polynomials.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 November 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ND, the number of data points.
%
%    Input, real XD(ND,1), the data nodes.
%
%    Input, integer NI, the number of evaluation points.
%
%    Input, real XI(NI,1), the evaluation points.
%
%    Output, real LP(NI,ND), the value, at the I-th point XI, of the
%    Jth basis function.
%
  lp = zeros ( ni, nd );

  for i = 1 : ni
    for j = 1 : nd

      for j1 = 1 : nd

        if ( j1 ~= j )
          p = 1.0;
          for j2 = 1 : nd
            if ( j2 == j1 )
              p = p / ( xd(j) - xd(j2) );
            elseif ( j2 ~= j )
              p = p * ( xi(i) - xd(j2) ) / ( xd(j) - xd(j2) );
            end
          end
          lp(i,j) = lp(i,j) + p;
        end

      end

    end
  end

  return
end


text:
% Special wrapper to regular plot
% Works on regularly sampled data
% When there is a gap of data, it insert a nan value
% to avoid linear interpolation of the data
%
% INPUT
%   t           column array of epochs
%   data        columns of data (could be a matrix)
%   varagin     add other useful parameters of the plot
%
% SYNTAX
%   Core_Utils.plotSep(t, data, varagin);
%
% SEE ALSO
%   plot

%--------------------------------------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0 beta 5
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Andrea Gatti
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

function lh = plotSep(t, data, varargin)
    if not(isempty(data))
        if isa(t, 'handle')
            ax = t;
            t = data;
            data = varargin{1};
            varargin = varargin{2:end};
        else
            ax = gca;
        end
        try
            if numel(t) ~= numel(data) && numel(t) ~= size(data, 1)
                % this means that there is only one data parameter and no time
                varargin = [{data} varargin];
                data = t;
                if size(data, 1) == 1
                    % I want the data to be columnwise
                    data = data';
                end
                t = 1 : size(data, 1);
            else
                if size(data, 1) == 1
                    % I want the data to be columnwise
                    data = data';
                end
            end
        catch
            % probably data is undefined
            data = t;
            if size(data, 1) == 1
                % I want the data to be columnwise
                data = data';
            end
            t = (1 : size(data, 1))';
        end
        for c = 1 : size(data, 2)
            if numel(data) == numel(t)
                [t_col, data_col] = insertNan4Plots(t(:,c), data(:,c));
            else
                [t_col, data_col] = insertNan4Plots(t, data(:,c));
            end
            if isempty(varargin)
                lh = plot(ax, t_col, data_col);
            else
                lh = plot(ax, t_col, data_col, varargin{:});
            end
            hold on;
        end
    end
end

function [t, data_set] = insertNan4Plots(t, data_set)
    % Insert a Nan in a regularly sampled dataset to make
    % plots interrupt continuous lines
    %
    % INPUT
    %   t      epoch of the data [matrix of column arrays]
    %   data   epoch of the data [matrix of column arrays]
    %
    % SYNTAX
    %   [t, data] = Core_Utils.insertNan4Plots(t, data)

    t = t(:);
    if size(t, 1) ~= size(data_set, 1)
        % data should be a column array
        data_set = data_set';
    end
    n_set = size(data_set, 2);
    dt = diff(t);
    rate = median(dt);
    id_in = find(dt > 1.5 * rate);
    for x = numel(id_in) : -1 : 1
        t = [t(1 : id_in(x)); (t(id_in(x)) + 1.5 * rate); t((id_in(x)+1) : end)];
        data_set = [data_set(1 : id_in(x), :); nan(1, n_set); data_set((id_in(x)+1) : end, :)];
    end
end


text:
function [data clustPoints idx centers slopes lengths] = generateData( ...
    slope, ...
    slopeStd, ...
    numClusts, ...
    xClustAvgSep, ...
    yClustAvgSep, ...
    lengthAvg, ...
    lengthStd, ...
    lateralStd, ...
    totalPoints ...
    )
% generateData Generates 2D data for clustering; data is created along
%              straight lines, which can be more or less parallel depending
%              on slopeStd argument.
%
% Inputs:
%        slope - Base direction of the lines on which clusters are based.
%     slopeStd - Standard deviation of the slope; used to obtain a random
%                slope variation from the normal distribution, which is
%                added to the base slope in order to obtain the final slope
%                of each cluster.
%    numClusts - Number of clusters (and therefore of lines) to generate.
% xClustAvgSep - Average separation of line centers along the X axis.
% yClustAvgSep - Average separation of line centers along the Y axis.
%                line centers along each dimension.
%    lengthAvg - The base length of lines on which clusters are based.
%    lengthStd - Standard deviation of line length; used to obtain a random
%                length variation from the normal distribution, which is
%                added to the base length in order to obtain the final
%                length of each line.
%   lateralStd - "Cluster fatness", i.e., the standard deviation of the
%                distance from each point to the respective line, in both x
%                and y directions; this distance is obtained from the
%                normal distribution.
%  totalPoints - Total points in generated data (will be
%                randomly divided among clusters).
%
% Output:
%         data - Matrix (totalPoints x 2) with the generated data
%  clustPoints - Vector (numClusts x 1) containing number of points in each
%                cluster
%          idx - Vector (totalPoints x 1) containing the cluster indices of
%                each point
%      centers - Matrix (numClusts x 2) containing centers from where
%                clusters were generated
%       slopes - Vector (numClusts x 1) containing the effective slopes
%                used to generate clusters
%      lengths - Vector (numClusts x 1) containing the effective lengths
%                used to generate clusters
%
% ----------------------------------------------------------
% Usage example:
%
%   [data cp idx] = generateData(1, 0.5, 5, 15, 15, 5, 1, 2, 200);
%
% This creates 5 clusters with a total of 200 points, with a base slope
% of 1 (std=0.5), separated in average by 15 units in both x and y
% directions, with average length of 5 units (std=1) and a "fatness" or
% spread of 2 units.
%
% To take a quick look at the clusters just do:
%
%   scatter(data(:,1), data(:,2), 8, idx);
% ----------------------------------------------------------
%
%  N. Fachada
%  Instituto Superior Técnico
%  Aug 10, 2010
%  Updated: Aug 2, 2012

% Make sure totalPoints >= numClusts
if totalPoints < numClusts
    error('Number of points must be equal or larger than the number of clusters.');
end;

% Determine number of points in each cluster
clustPoints = abs(randn(numClusts, 1));
clustPoints = clustPoints / sum(clustPoints);
clustPoints = round(clustPoints * totalPoints);

% Make sure totalPoints is respected
while sum(clustPoints) < totalPoints
    % If one point is missing add it to the smaller cluster
    [C,I] = min(clustPoints);
    clustPoints(I(1)) = C + 1;
end;
while sum(clustPoints) > totalPoints
    % If there is one extra point, remove it from larger cluster
    [C,I] = max(clustPoints);
    clustPoints(I(1)) = C - 1;
end;

% Make sure there are no empty clusters
emptyClusts = find(clustPoints == 0);
if ~isempty(emptyClusts)
    % If there are empty clusters...
    numEmptyClusts = size(emptyClusts, 1);
    for i=1:numEmptyClusts
        % ...get a point from the largest cluster and assign it to the
        % empty cluster
        [C,I] = max(clustPoints);
        clustPoints(I(1)) = C - 1;
        clustPoints(emptyClusts(i)) = 1;
    end;
end;

% Initialize data matrix
data = zeros(sum(clustPoints), 2);

% Initialize idx (vector containing the cluster indices of each point)
idx = zeros(totalPoints, 1);

% Initialize lengths vector
lengths = zeros(numClusts, 1);

% Determine cluster centers
xCenters = xClustAvgSep * numClusts * (rand(numClusts, 1) - 0.5);
yCenters = yClustAvgSep * numClusts * (rand(numClusts, 1) - 0.5);
centers = [xCenters yCenters];

% Determine cluster slopes
slopes = slope + slopeStd * randn(numClusts, 1);

% Create clusters
for i=1:numClusts
    % Determine length of line where this cluster will be based
    lengths(i) = abs(lengthAvg + lengthStd*randn);
    % Determine how many points have been assigned to previous clusters
    sumClustPoints = 0;
    if i > 1
        sumClustPoints = sum(clustPoints(1:(i - 1)));
    end;
    % Create points for this cluster
    for j=1:clustPoints(i)
        % Determine where in the line the next point will be projected
        position = lengths(i) * rand - lengths(i) / 2;
        % Determine x coordinate of point projection
        delta_x = cos(atan(slopes(i))) * position;
        % Determine y coordinate of point projection
        delta_y = delta_x * slopes(i);
        % Get point distance from line in x coordinate
        delta_x = delta_x + lateralStd * randn;
        % Get point distance from line in y coordinate
        delta_y = delta_y + lateralStd * randn;
        % Determine the actual point
        data(sumClustPoints + j, :) = [(xCenters(i) + delta_x) (yCenters(i) + delta_y)];
    end;
    % Update idx
    idx(sumClustPoints + 1 : sumClustPoints + clustPoints(i)) = i;
end;

author: Sable


text:
% pop_writeeeg - write EEGLAB dataset to disk in EDF/GDF or BDF format
%
% pop_writeeeg( EEG ) % pops up a window
% pop_writeeeg( EEG, filename, 'key', 'val' )
%
% Inputs:
%  EEG        - EEGLAB dataset
%  filename   - [string] filename
%
% Optional keys (same as writeeeg):
%  'TYPE'         - ['GDF'|'EDF'|'BDF'|'CFWB'|'CNT'] file format for writing
%                   default is 'EDF'.
%  See writeeeg for more information
%
% Author: Arnaud Delorme, SCCN, UCSD/CERCO, 2009
%         Based on BIOSIG, sopen and swrite

% Copyright (C) 22 March 2002 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [command] = pop_writeeeg(EEG, filename, varargin)
command = '';

% enforce use the str2double of Biosig
biosigpathfirst

% Check for BIOSIG Toolbox
global PLUGINLIST
 if ~isempty(PLUGINLIST) && ~any(strcmpi({PLUGINLIST.plugin}','biosig'))
     fprintf(2,'pop_writeeeg error: This function requires you to install BIOSIG Toolbox plug-in. \n');
     return;
 elseif isempty(PLUGINLIST)
     warning('pop_writeeeg check for BIOSIG Toolbox could not be performed.');
 end

if nargin < 2
    if EEG.trials > 1
        res = questdlg2( [ 'This dataset contains data epochs.' 10 'Do you want to export the concatenated' 10 'data epochs?' ], '', 'No', 'Yes', 'Yes');
        if strcmpi(res, 'No')
            return;
        end;
    end;

	% ask user
	[filename, filepath] = uiputfile('*.*', 'Enter a file name -- pop_writeeeg()');
	if filename == 0 return; end;
	filename = fullfile(filepath,filename);

    % file format
    % -----------
    fileformats = { 'EDF' 'GDF' 'BDF' };
    uilist = { { 'style' 'text' 'String' 'File format' } ...
               { 'style' 'listbox' 'string' strvcat(fileformats) 'value' 1 } };
    geom = [1 1];
    result = inputgui( 'geometry', geom, 'uilist', uilist, 'helpcom', 'pophelp(''pop_writeeeg'')', ...
                     'title', 'Write data using BIOSIG -- pop_writeeeg()', 'geomvert', [1 2.5]);
    if length(result) == 0 return; end;

    options = { 'TYPE' fileformats{result{1}} };
else
    options = varargin;
end;

warning('off', 'MATLAB:intConvertNonIntVal');
if ~isempty(EEG.chanlocs)
    tmpchanlocs = EEG.chanlocs;
    writeeeg(filename, EEG.data(:,:), EEG.srate, 'label', { tmpchanlocs.labels }, 'EVENT', EEG.event, options{:});
else
    writeeeg(filename, EEG.data(:,:), EEG.srate, 'EVENT', EEG.event, options{:});
end;
warning('on', 'MATLAB:intConvertNonIntVal');

command = sprintf('pop_writeeeg(EEG, ''%s'', %s);', filename, vararg2str(options));
biosigpathlast

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function varargout = plot(sF,varargin)
%
% Syntax
%
%   plot(sF)
%
% Input
%  sF - @S2FunHarmonicSym
%
% See also
% S2Fun/plot

% create a new figure if needed
%[mtexFig,isNew] = newMtexFigure('datacursormode',@tooltip,varargin{:});

% get plotting region
if sF.antipodal, varargin = [varargin,'antipodal']; end
sR = sF.s.fundamentalSector(varargin{:});

% set axes alignment
if isa(sF.s,'crystalSymmetry')
  varargin = [sF.s.plotOptions,varargin];
end

% perform plotting
[varargout{1:nargout}] = sF.plot@S2Fun(sR,sF.s,varargin{:});

function txt = tooltip(varargin)

[h_local,~,value] = getDataCursorPos(mtexFig);

h_local = Miller(h_local,sF.s,'uvw');
h_local = round(h_local,'tolerance',3*degree);
txt = [xnum2str(value) ' at ' char(h_local)];

end

end

author: mtex-toolbox

repo: mtex

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


text:
function [c,ceq] = nl_terminalconstraints(t, x, mpcModel , varargin)

    c(1) = x(2)-mpcModel.battery.range(2);
    c(2) = mpcModel.battery.range(1)-x(2);

    ceq = [];
end

author: juchengquan


text:
function str = printstruct(name, val, varargin)

% PRINTSTRUCT converts a MATLAB structure into a multiple-line string that, when
% evaluated by MATLAB, results in the original structure. It also works for most
% other standard MATLAB classes, such as numbers, vectors, matrices, and cell-arrays.
%
% Use as
%   str = printstruct(val)
% or
%   str = printstruct(name, val)
% where "val" is any MATLAB variable, e.g. a scalar, vector, matrix, structure, or
% cell-array. If you pass the name of the variable, the output is a piece of MATLAB code
% that you can execute, i.e. an ASCII serialized representation of the variable.
%
% Example
%   a.field1 = 1;
%   a.field2 = 2;
%   s = printstruct(a)
%
%   b = rand(3);
%   s = printstruct(b)
%
%   s = printstruct('c', randn(10)>0.5)
%
% See also DISP

% Copyright (C) 2006-2018, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if nargin==1
  val  = name;
  name = inputname(1);
end

% get the options and set defaults
transposed    = ft_getopt(varargin, 'transposed', false);     % print a column as a transposed row, or a row as a transposed column
linebreaks    = ft_getopt(varargin, 'linebreaks', true);      % print with or without linebreaks
lastnewline   = ft_getopt(varargin, 'lastnewline', false);    % end with a newline
lastsemicolon = ft_getopt(varargin, 'lastsemicolon', true);   % end with a semicolon

% when called recursively, it should keep the last newline and semicolon intact
varargin = ft_setopt(varargin, 'lastnewline', true);
varargin = ft_setopt(varargin, 'lastsemicolon', true);

% Note that because we don't know the final size of the string, iteratively appending
% is actually faster than creating a cell-array and subsequently doing a cat(2,
% strings{:}). Note also that sprintf() is slow.
str = '';

% Note further that in the string concatenations I use the numerical value
% of a newline (\n), which is 10.

if numel(val) == 0
  if iscell(val)
    str = [name ' = {};' 10];
  elseif isstruct(val)
    str = [name ' = [];' 10];
  elseif isnumeric(val)
    str = [name ' = [];' 10];
  end

elseif isstruct(val)
  if numel(val)>1
    % print it as a struct-array
    str = cell(size(val));
    for i=1:numel(val)
      str{i} = printstruct(sprintf('%s(%d)', name, i), val(i), varargin{:});
    end
    str = cat(2, str{:});
  elseif numel(val)==1
    % print it as a named structure
    fn = fieldnames(val);
    for i=1:length(fn)
      fv = val.(fn{i});
      switch class(fv)
        case {'char' 'string'}
          line = printstr([name '.' fn{i}], fv, linebreaks, transposed);
        case {'single' 'double' 'int8' 'int16' 'int32' 'int64' 'uint8' 'uint16' 'uint32' 'uint64' 'logical'}
          if ismatrix(fv)
            line = printmat([name '.' fn{i}], fv, linebreaks, transposed);
          else
            line = '''ERROR: multidimensional arrays are not supported''';
          end
        case 'cell'
          line = printcell([name '.' fn{i}], fv, linebreaks, transposed);
        case 'struct'
          line = [printstruct([name '.' fn{i}], fv, varargin{:}) 10];
        case 'function_handle'
          line = printstr([name '.' fn{i}], func2str(fv), linebreaks, transposed);
        otherwise
          ft_error('unsupported');
      end
      str = [str line newline];
    end % for fieldnames
    if isempty(fn)
      % print it as an empty structure
      str = [name ' = struct();'];
    end
  end

elseif ~isstruct(val)
  % print it as a named variable
  switch class(val)
    case {'char' 'string'}
      str = printstr(name, val, linebreaks, transposed);
    case {'double' 'single' 'int8' 'int16' 'int32' 'int64' 'uint8' 'uint16' 'uint32' 'uint64' 'logical'}
      str = printmat(name, val, linebreaks, transposed);
    case 'cell'
      str = printcell(name, val, linebreaks, transposed);
    otherwise
      ft_error('unsupported');
  end
end

% replace double newlines with a single one
str = strrep(str, [newline newline], newline);

if isempty(name)
  % remove the ' = ' part from the string
  str = str(4:end);
end

if str(end)==newline && ~lastnewline
  % remove the last newline
  str = str(1:end-1);
end

if str(end)==';' && ~lastsemicolon
  % remove the last semicolon
  str = str(1:end-1);
end

if ~nargout
  disp(str)
  clear str
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function str = printstr(name, val, linebreaks, transposed)
siz = size(val);
if siz(1)>1
  str = sprintf('%s = \n', name);
  for i=1:siz(1)
    str = [str sprintf('  %s\n', printval(val(i,:), linebreaks, transposed))];
  end
elseif siz(1)==1
  str = sprintf('%s = %s;\n', name, printval(val, linebreaks, transposed));
else
  str = sprintf('%s = '''';\n', name);
end
if ~linebreaks
  str(str==10) = [];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function str = printcell(name, val, linebreaks, transposed)
siz = size(val);
if isempty(val)
  str = sprintf('%s = {};\n', name);
elseif numel(siz) == 2 && siz(2) == 1 && transposed
  % print column vector as (non-conjugate) transposed row
  str = printcell(name, transpose(val), linebreaks, false);
  str = sprintf('%s'';\n', str(1:end-2));
elseif numel(siz) == 2 && siz(1) == 1 && transposed
  % print row vector as (non-conjugate) transposed column
  str = printcell(name, transpose(val), linebreaks, false);
  str = sprintf('%s'';\n', str(1:end-2));
else
  str = sprintf('%s = {\n', name);
  for i=1:siz(1)
    for j=1:(siz(2)-1)
      dum = printval(val{i,j}, linebreaks, transposed);
      if endsWith(dum, ';')
        dum = dum(1:end-1);
      elseif endsWith(dum, [';' newline])
        dum = dum(1:end-2);
      end
      str = [str ' ' dum ',']; % add the element with a comma
    end
    dum = printval(val{i,siz(2)}, linebreaks, transposed);
    if endsWith(dum, ';')
      dum = dum(1:end-1);
    elseif endsWith(dum, [';' newline])
      dum = dum(1:end-2);
    end
    str = [str ' ' dum 10];  % add the last element of each row with a linebreak
  end
  if siz(1)==1
    str(str==10) = [];
  end
  str = sprintf('%s };\n', str);
end
if ~linebreaks
  str(str==10) = [];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function str = printmat(name, val, linebreaks, transposed)
siz = size(val);
if numel(val) == 0
  str = sprintf('%s = [];\n', name);
elseif numel(val) == 1
  % an integer will never get trailing decimals when using %g
  str = sprintf('%s = %g;', name, val);
elseif numel(siz) == 2 && siz(2) == 1 && transposed
  % print column vector as (non-conjugate) transposed row
  str = printmat(name, transpose(val), linebreaks, false);
  str = sprintf('%s''', str);
elseif numel(siz) == 2 && siz(1) == 1 && transposed
  % print row vector as (non-conjugate) transposed column
  str = printmat(name, transpose(val), linebreaks, false);
  str = sprintf('%s''', str);
elseif ismatrix(val)
  if isa(val, 'double')
    str = sprintf('%s = %s;', name, mat2str(val));
  else
    % add class information for non-double numeric matrices
    str = sprintf('%s = %s;', name, mat2str(val, 'class'));
  end
  % add a newline after each semicolon
  str = strrep(str, ';', [';' newline]);
else
  ft_warning('multidimensional arrays are not supported');
  str = '''ERROR: multidimensional arrays are not supported''';
end
if ~linebreaks
  str(str==10) = [];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function str = printval(val, linebreaks, transposed)

switch class(val)
  case 'char'
    str = ['''' val ''''];

  case 'string'
    % these use " in the declaration rather than '
    str = ['"' char(val) '";'];

  case {'single' 'double' 'int8' 'int16' 'int32' 'int64' 'uint8' 'uint16' 'uint32' 'uint64' 'logical'}
    str = printmat('', val, linebreaks, transposed);
    str = str(4:end);
    if endsWith(str, ';')
      str = str(1:end-1);
    end

  case 'function_handle'
    str = ['@' func2str(val)];

  case 'struct'
    % print it as an anonymous structure
    str = 'struct(';
    fn = fieldnames(val);
    for i=1:numel(fn)
      str = [str '''' fn{i} '''' ', ' printval(val.(fn{i}), linebreaks, transposed)];
    end
    str = [str ')'];

  otherwise
    ft_warning('cannot print unknown object at this level');
    str = '''ERROR: cannot print unknown object at this level''';
end
if ~linebreaks
  str(str==10) = [];
end

author: fieldtrip

repo: fieldtrip


text:
function uv = compute_flow(this, init, gt)
%
%COMPUTE_FLOW   Compute flow field
%   UV = COMPUTE_FLOW(THIS[, INIT]) computes the flow field UV with
%   algorithm THIS and the optional initialization INIT.
%
%   This is a member function of the class 'ba_optical_flow'.
%
% Authors: Deqing Sun, Department of Computer Science, Brown University
%          Stefan Roth, Department of Computer Science, TU Darmstadt
% Contact: dqsun@cs.brown.edu, sroth@cs.tu-darmstadt.de
% $Date: $
% $Revision: $
%
% Copyright 2007-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

  % Frame size
  sz = [size(this.images, 1), size(this.images, 2)];

  % If we have no initialization argument, initialize with all zeros
  if (nargin < 2)
    uv = zeros([sz, 2]);
  else
    uv = init;
  end

  % Preprocess input (gray) sequences
  if this.texture
      % Perform ROF structure texture decomposition
      images  = structure_texture_decomposition_rof( this.images, 1/8, 100, this.alp);

  elseif this.fc

      % Laplacian in flowfusion
      f = fspecial('gaussian', [5 5], 1.5);
      % Linear combine Laplacian image with original images
      images  = this.images- this.alp*imfilter(this.images, f, 'symmetric');

      % Gaussian pre-filering
      %         f = fspecial('gaussian', [3 3], 0.5);  % Li & Huttenlocher
      %         images  = imfilter(this.images, f, 'symmetric');

      % Sobel edge magnitude
      %         Dy = fspecial('sobel')/8;
      %         Dx = Dy';
      %         Ix = imfilter(this.images, Dx, 'symmetric');
      %         Iy = imfilter(this.images, Dy, 'symmetric');
      %         images = sqrt(Ix.^2 + Iy.^2);

      images  = scale_image(images, 0, 255);
  else
      images  = scale_image(this.images, 0, 255);
  end;

  if this.auto_level
      % Automatic determine pyramid level
      this.pyramid_levels  =  1 + floor( log(min(size(images, 1), size(images,2))/16) / log(this.pyramid_spacing) );
  end;

  % Construct image pyramid, using filter setting in Bruhn et al in "Lucas/Kanade.." (IJCV2005') page 218

  % For gnc stage 1
  factor            = sqrt(2);  % sqrt(3) worse
  smooth_sigma      = sqrt(this.pyramid_spacing)/factor;   % or sqrt(3) recommended by Manuel Werlberger
  f                 = fspecial('gaussian', 2*round(1.5*smooth_sigma) +1, smooth_sigma);
  pyramid_images    = compute_image_pyramid(images, f, this.pyramid_levels, 1/this.pyramid_spacing);

  % For gnc stage 2 to gnc_iters
  smooth_sigma      = sqrt(this.gnc_pyramid_spacing)/factor;
  f                 = fspecial('gaussian', 2*round(1.5*smooth_sigma) +1, smooth_sigma);
  gnc_pyramid_images= compute_image_pyramid(images, f, this.gnc_pyramid_levels, 1/this.gnc_pyramid_spacing);

  tic

  for ignc = 1:this.gnc_iters

      if this.display
          disp(['GNC stage: ', num2str(ignc)])
      end

      if ignc == 1
          pyramid_levels = this.pyramid_levels;
      else
          pyramid_levels = this.gnc_pyramid_levels;
      end;

      % Iterate through all pyramid levels starting at the top
      for l = pyramid_levels:-1:1

          if this.display
              disp(['-Pyramid level: ', num2str(l)])
          end

          % Generate copy of algorithm with single pyramid level and the
          % appropriate subsampling
          small = this;

          if ignc == 1
              nsz               = [size(pyramid_images{l}, 1) size(pyramid_images{l}, 2)];
              small.images      = pyramid_images{l};
              small.max_linear  = 1;             % number of linearization performed per warping, 1 OK for quadratic formulation

          else
              small.images         = gnc_pyramid_images{l};
              nsz   = [size(gnc_pyramid_images{l}, 1) size(gnc_pyramid_images{l}, 2)];

          end;

          % Rescale the flow field
          uv        = resample_flow(uv, nsz);

          % Run flow method on subsampled images
          uv = compute_flow_base(small, uv);
      end

      % Update GNC parameters (linearly)
      if this.gnc_iters > 1
          new_alpha  = 1 - ignc / (this.gnc_iters-1);
          this.alpha = min(this.alpha, new_alpha);
          this.alpha = max(0, this.alpha);
      end;

      if true % this.display

          fprintf('energy of solution \t%3.3e\n', -evaluate_log_posterior(small, uv));
          if nargin == 3
              [aae stdae aepe] = flowAngErr(gt(:,:,1), gt(:,:,2), uv(:,:,1), uv(:,:,2), 0);        % ignore 0 boundary pixels
              fprintf('AAE %3.3f STD %3.3f average end point error %3.3f \n', aae, stdae, aepe);
          end;

      end;

      fprintf('GNC stage %d finished, %3.2f minutes passed \n', ignc, toc/60);
  end;

  % minaly improves H&S, slightly improves others
%   if ~isempty(this.median_filter_size)
%       uv(:,:,1) = medfilt2(uv(:,:,1), this.median_filter_size, 'symmetric');
%       uv(:,:,2) = medfilt2(uv(:,:,2), this.median_filter_size, 'symmetric');
%   end;

author: kristinbranson

repo: JAABA


text:
function O = imwarp(I, u, v, nopad)
%IMWARP   Warp image with flow field
%   O = IMWARP(I, U, V[, NOPAD]) warps the image I with optical flow field
%   U and V.  The flow field is to be given in the coordinate system of
%   O, i.e. the operation warps I toward O.  If the optional argument
%   NOPAD is given and true, the undefined pixels (source pixel outside
%   the boundary) are given by the nearest boundary pixels, otherwise
%   they are set to NaN.
%
%   Author:  Stefan Roth, Department of Computer Science, TU Darmstadt
%   Contact: sroth@cs.tu-darmstadt.de
%   $Date: 2007-03-27 14:09:11 -0400 (Tue, 27 Mar 2007) $
%   $Revision: 252 $

% Copyright 2004-2007, Brown University, Providence, RI. USA
% Copyright 2007-2010 TU Darmstadt, Darmstadt, Germany.
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

  % Image size
  sx = size(I, 2);
  sy = size(I, 1);

  % Image size w/ padding
  spx = sx + 2;
  spy = sy + 2;

  if (nargin > 3 & nopad)
    % Warped image coordinates
    [X, Y] = meshgrid(1:sx, 1:sy);
    XI = reshape(X + u, 1, sx * sy);
    YI = reshape(Y + v, 1, sx * sy);

    % Bound coordinates to valid region
    XI = max(1, min(sx - 1E-6, XI));
    YI = max(1, min(sy - 1E-6, YI));

    % Perform linear interpolation (faster than interp2)
    fXI = floor(XI);
    cXI = ceil(XI);
    fYI = floor(YI);
    cYI = ceil(YI);

    alpha_x = XI - fXI;
    alpha_y = YI - fYI;

    O = (1 - alpha_x) .* (1 - alpha_y) .* I(fYI + sy * (fXI - 1)) + ...
        alpha_x .* (1 - alpha_y) .* I(fYI + sy * (cXI - 1)) + ...
        (1 - alpha_x) .* alpha_y .* I(cYI + sy * (fXI - 1)) + ...
        alpha_x .* alpha_y .* I(cYI + sy * (cXI - 1));

  else
    % Pad image with NaNs
    Z = [NaN(1, sx+2); NaN(sy, 1), I, NaN(sy, 1); NaN(1, sx+2)];

    % Warped image coordinates in padded image
    [X, Y] = meshgrid(2:sx+1, 2:sy+1);
    XI = reshape(X + u, 1, sx * sy);
    YI = reshape(Y + v, 1, sx * sy);

    % Bound coordinates to valid region
    XI = max(1, min(spx - 1E-6, XI));
    YI = max(1, min(spy - 1E-6, YI));

    % Perform linear interpolation (faster than interp2)
    fXI = floor(XI);
    cXI = ceil(XI);
    fYI = floor(YI);
    cYI = ceil(YI);

    alpha_x = XI - fXI;
    alpha_y = YI - fYI;

    O = (1 - alpha_x) .* (1 - alpha_y) .* Z(fYI + spy * (fXI - 1)) + ...
        alpha_x .* (1 - alpha_y) .* Z(fYI + spy * (cXI - 1)) + ...
        (1 - alpha_x) .* alpha_y .* Z(cYI + spy * (fXI - 1)) + ...
        alpha_x .* alpha_y .* Z(cYI + spy * (cXI - 1));
  end

  O = reshape(O, sy, sx);