function resFileClean = preprocessResult(resFile, seqName, dataDir)
% reads submitted (raw) MOT16 result from .txt
% and removes all boxes that are associated with ambiguous annotations
% such as sitting people, cyclists or mannequins.

% resFile='/home/amilan/research/projects/bmtt-dev/code/bae/res/0001/ADL-Rundle-6.txt';
% seqName = 'MOT16-09';
assert(~isempty(strfind(seqName, 'MOT16')),'preproccessing should only be done for MOT16')

fprintf('Preprocessing (cleaning) %s...\n',seqName);

%  if nargin<3, dataDir = getDataDir; end

% [seqName, seqFolder, imgFolder, imgExt, F, dirImages] ...
%     = getSeqInfo(seq, dataDir);

[seqName, seqFolder, imgFolder, frameRate, F, imWidth, imHeight, imgExt] ...
    = getSeqInfoFromFile(seqName, dataDir);

% read in result
resRaw = dlmread(resFile);

%%

% read ground truth
gtFolder = [dataDir,seqName,filesep,'gt',filesep];
gtFile = [gtFolder,'gt.txt'];
gtRaw = dlmread(gtFile);

% make sure we have MOT16 ground truth (= 9 columns)
assert(size(gtRaw,2)==9, 'unknown GT format')

% define which classes should be ignored
distractors = {'person_on_vhcl','static_person','distractor','reflection'};

keepBoxes = true(size(resRaw,1),1);

showVis = 0;

td=0.5; % overlap threshold
for t=1:F
    if ~mod(t,100), fprintf('.'); end

    % find all result boxes in this frame
    resInFrame = find(resRaw(:,1)==t); N = length(resInFrame);
    resInFrame = reshape(resInFrame,1,N);

    % find all GT boxes in frame
    GTInFrame = find(gtRaw(:,1)==t); Ngt = length(GTInFrame);
    GTInFrame = reshape(GTInFrame,1,Ngt);

    % compute all overlaps for current frame
    allisects=zeros(Ngt,N);
    g=0;
    for gg=GTInFrame
        g=g+1; r=0;
        bxgt=gtRaw(gg,3); bygt=gtRaw(gg,4); bwgt=gtRaw(gg,5); bhgt=gtRaw(gg,6);
        for rr=resInFrame
            r=r+1;
            bxres=resRaw(rr,3); byres=resRaw(rr,4); bwres=resRaw(rr,5); bhres=resRaw(rr,6);

            if bxgt+bwgt<bxres, continue; end % ignore if no horizontal overlap
            if bxgt>bxres+bwres, continue; end

            if bygt+bhgt<byres, continue; end % ignore if no vertical overlap
            if bygt>byres+bhres, continue; end

            allisects(g,r)=boxiou(bxgt,bygt,bwgt,bhgt,bxres,byres,bwres,bhres);
        end
    end
%     t

    tmpai=allisects;
    tmpai=1-tmpai;
    tmpai(tmpai>td)=Inf;
    [Mtch,Cst]=Hungarian(tmpai);
    [mGT,mRes]=find(Mtch);
%     pause
    nMtch = length(mGT);
    % go through all matches
    for m=1:nMtch
        g=GTInFrame(mGT(m)); % gt box
        gtClassID = gtRaw(g,8);
        gtClassString = classIDToString(gtClassID);

        % if we encounter a distractor, mark to remove box
        if ismember(gtClassString, distractors)
            r = resInFrame(mRes(m)); % result box
            keepBoxes(r) = false;

            if showVis
                bxgt=gtRaw(g,3); bygt=gtRaw(g,4); bwgt=gtRaw(g,5); bhgt=gtRaw(g,6); idgt=gtRaw(g,2);
                bxres=resRaw(r,3); byres=resRaw(r,4); bwres=resRaw(r,5); bhres=resRaw(r,6); idres=resRaw(r,2);

                clf
                im = imread(fullfile(dataDir,seqName,'img1',sprintf('%06d.jpg',t)));
                imshow(im); hold on
                text(50,50,sprintf('%d',t),'color','w')

                % show GT box
%                 text(bxgt,bygt-20,sprintf('%d',idgt),'color','w')
                classString = insertEscapeChars(classIDToString(gtClassID));
                text(bxgt+50,bygt-20,sprintf('%s',classString),'color','w')
                rectangle('Position',[bxgt,bygt,bwgt,bhgt],'EdgeColor','w');

                % show Res box
%                 text(bxres,byres-20,sprintf('%d',idres),'color','y')
                rectangle('Position',[bxres,byres,bwres,bhres],'EdgeColor','y');

                pause(.01)
            end
        end
    end
end

%%
fprintf('\nRemoving %d boxes from solution...\n',numel(find(~keepBoxes)));
resNew = resRaw;
resNew=resNew(keepBoxes,:);

%% leave original copy in a separate dir (preclean) if not already there
[p,f,e]=fileparts(resFile);
cleanDir = [p,filesep,'clean'];
if ~exist(cleanDir, 'dir'), mkdir(cleanDir); end
resFileClean = [cleanDir,filesep,f,e];
dlmwrite(resFileClean, resNew);

% end


text:
function x = stack(X)

% STACK Return column stacked vector of given matrix.
% FORMAT
% DESC returns a column stacked vector of a given matrix. This is
% useful if you wish to stack a column vector from a matrix
% returned by another function (i.e. you can't apply the colon
% operator directly).
% ARG X : the matrix to be stacked.
% RETURN x : stacked column vector from the matrix.
%
% SEEALSO : colon
%
% COPYRIGHT : Neil D. Lawrence, 2004

% NDLUTIL

x = X(:);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% DISPLAY - display an EEG data class underlying structure
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, Nov. 2008

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

function b = display(a)

    i.type = '()';
    i.subs = { ':' ':' ':' };
    b = subsref(a, i); % note that subsref cannot be called directly
    return;

    %struct(a)
    %return;

    if ~strcmpi(a.fileformat, 'transposed')
        a.data.data.x;
    else
        permute(a, [3 1 2]);
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function damped_cosine_test ( )

%*****************************************************************************80
%
%% DAMPED_COSINE_TEST examines the damped_cosine correlation.
%
%  Discussion:
%
%    This code is based substantially on a document by Toby Driscoll.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Toby Driscoll,
%    Mercer's theorem and the Karhunen-Loeve expansion,
%    Oxford University Mathematical Institute,
%    http://www2.maths.ox.ac.uk/chebfun/examples/stats/pdf/MercerKarhunenLoeve.pdf
%
  rmpath ( '/usr/local/matlab/toolbox/datafeed/datafeed' )
  addpath ( '../chebfun' )

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'DAMPED_COSINE_TEST:\n' );
  fprintf ( 1, '  Demonstrate Mercer''s theorem and the KL expansion\n' );
  fprintf ( 1, '  for the damped_cosine kernel.\n' );
%
%  Set the interval.
%
  a = 0.0;
  b = +10.0;
  fprintf ( 1, '  Using interval [%g,%g]\n', a, b );
  s_num = 21;
  s_vec = linspace ( a, b, s_num );
%
%  FRED is a function from the CHEBFUN library.
%
%  It constructs a "chebop" representing the Fredholm integral operator
%  with kernel K for functions in domain [A,B].
%
  F = fred ( @damped_cosine_correlation, domain ( [ a, b ] ) );
%
%  EIGS has been extended to be able to compute the eigenvalues Lambda
%  and eigenfunctions Psi of the Fredholm integral operator represented by F.
%
%  The "LM" switch requests that we return the eigenvalues of largest magnitude.
%
%  Each Psi is a CHEBFUN, that is, it takes a real number argument.
%
  eigen_num = 20;
  [ Psi, Lambda ] = eigs ( F, eigen_num, 'lm' );

  eigen_found = length ( Lambda );
  fprintf ( 1, '  Requested %d eigenmodes, computed %d\n', eigen_num, eigen_found );
  eigen_num = min ( eigen_num, eigen_found );
%
%  Print the eigenvalues.
%
  lambda_vec = diag ( Lambda );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   I       Lambda(I)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : eigen_num
    fprintf ( 1, '  %2d  %10g\n', i, lambda_vec(i) );
  end
%
%  Plot the eigenvalues.
%
  figure ( 1 )
  clf
  plot ( lambda_vec, 'Linewidth', 2 );
  hold on
  plot ( lambda_vec, 'b.', 'Markersize', 20 );
  title ( 'damped_cosine: Mercer eigenvalues' );
  xlabel ( '<--- N --->')
  grid on
  print -dpng 'damped_cosine_figure1.png'
%
%  Plot selected eigenfunctions.
%
  figure ( 2 )
  subplot ( 4, 1, 1 )
  plot ( Psi(:,1), 'Linewidth', 2 );
  title ( 'damped_cosine: Mercer eigenfunction PSI(1)' )
  grid on
  subplot ( 4, 1, 2 )
  plot ( Psi(:,2), 'Linewidth', 2 );
  title ( 'damped_cosine: Mercer eigenfunction PSI(2)' )
  grid on
  subplot ( 4, 1, 3 )
  plot ( Psi(:,5), 'Linewidth', 2 );
  title ( 'damped_cosine: Mercer eigenfunction PSI(5)' )
  grid on
  subplot ( 4, 1, 4 )
  plot ( Psi(:,10), 'Linewidth', 2 );
  title ( 'damped_cosine: Mercer eigenfunction PSI(10)' )
  grid on
  print -dpng 'damped_cosine_figure2.png'
%
%  Orthonormality check.
%
  ptp = Psi' * Psi;
  error_frobenius = r8mat_is_identity ( eigen_num, ptp );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Frobenius norm of I - Psi'' * Psi = %g\n', error_frobenius );
%
%  K(S,S) should be exactly 1.
%  Because we are using a truncated representation of K, our estimate of K(S,S)
%  will be smaller than 1.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Truncated estimate of K(s,s) = 1 for S in the interval.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '        S        K(s,s) estimate\n' );
  fprintf ( 1, '\n' );
  for i = 1 : 21
    s = s_vec(i);
    ptlp = Psi(s,:) * Lambda * Psi(s,:)';
    fprintf ( 1, '  %10g  %14g\n', s, ptlp );
  end
%
%  Look at eigenvalue decay.
%
  x = diff ( log ( ( 1:eigen_num ) ) );
  y = diff ( log ( lambda_vec ) )';
  c = y ./ x;
  figure ( 3 );
  clf
  plot ( c, 'Linewidth', 2 );
  grid on
  hold on
  plot ( c, 'b.', 'Markersize', 20 );
  xlabel ( '<-- N -->' )
  title ( 'damped_cosine: Eigenvalue decay rate' );
  print -dpng 'damped_cosine_figure3.png'
%
%  Look at eigenvalue sum.
%
%  The trace of K(s,s) over [a,b] should be the integral of 1 over [a,b],
%  that is, b - a.
%
%  We compare the trace to the partial sums of the lambda's to see how much
%  of the variance of the process we have captured.
%
  trace_K = b - a;
  lambda_cum = cumsum ( lambda_vec ) / trace_K;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Index    Cumulative Eigenvalue sum\n' );
  fprintf ( 1, '\n' );
  for i = 1 : eigen_num
    fprintf ( 1, '  %5d  %10g\n', i, lambda_cum(i) );
  end
%
%  We decide to use the first 10 eigenfunctions.
%
  eigen_use = 10;
%
%  Find 400 realizations of the process by selecting, for each realization,
%  10 random parameters Z in the truncated KL expansion.
%
  Z = randn ( eigen_use, 400 );
  X = Psi(:,1:eigen_use) * ( sqrt ( Lambda(1:eigen_use,1:eigen_use) ) * Z );
%
%  Plot 40 of the realizations;
%  Plot their mean, computed from all 400.
%
  figure ( 4 )
  clf
  plot ( X(:,1:40) )
  mu = sum ( X, 2 ) / 400;
  hold on;
  plot ( mu, 'k', 'linewidth', 3 );
  title ( 'damped_cosine: 40 Random Realizations X(t,omega), and their Mean.' )
  print -dpng 'damped_cosine_figure4.png'
%
%  Estimate the covariance from the data.
%
  figure ( 5 )
  clf
  [ S, T ] = meshgrid ( s_vec, s_vec );
  C = cov ( X(s_vec,:)' );
  mesh ( S, T, C );
  hold on;
  D = damped_cosine_correlation ( S, T );
  plot3 ( S, T, D, 'k.', 'markersize', 10 )
  title ( 'damped_cosine: Covariance K(S,T) (dots), and Estimate from Realizations (mesh)' )
  print -dpng 'damped_cosine_figure5.png'
%
%  Using just 10 functions in the expansion,
%  reduce the correlation length,
%  and examine the sum of the lambda's.
%
  if ( 0 )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Use a fixed number of eigenfunctions = %d\n', eigen_use );
  fprintf ( 1, '  but vary the correlation length RHOBAR.\n' );
  fprintf ( 1, '  (We used RHOBAR = 1 above.)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The sum of the eigenvalues, divided by (B-A),\n' );
  fprintf ( 1, '  discloses the relative amount of the total variation\n' );
  fprintf ( 1, '  that is captured by the truncated expansion.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      RHOBAR      VARSUM\n' );
  fprintf ( 1, '\n' );
  rhobar = 4.0;
  for i = 1 : 10
    K = @ ( s, t ) exp ( - abs ( s - t ) / rhobar );
    F = fred ( K, domain ( [ a, b ] ) );
    lambda_vec = eigs ( F, 20, 'lm' );
    varsum = sum ( lambda_vec(1:eigen_use) ) / ( b - a );
    fprintf ( 1, '  %10g  %10g\n', rhobar, varsum );
    rhobar = rhobar / 2.0;
  end

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'DAMPED_COSINE_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../chebfun' )

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [u,Du,eqn,info] = PoissonWG(node,elem,pde,bdFlag,option)
%% POISSONWG Poisson equation: lowest order weak Galerkin element
%
%   u = POISSONWG(node,elem,pde,bdFlag) produces the linear finite element
%   approximation of the Poisson equation
%
%       -div(d*grad(u))=f  in \Omega, with
%       Dirichlet boundary condition u=g_D on \Gamma_D,
%       Neumann boundary condition   d*grad(u)*n=g_N on \Gamma_N,
%       Robin boundary condition     g_R*u + d*grad(u)*n=g_N on \Gamma _R
%
%   The mesh is given by node and elem and the boundary edge is given by
%   bdFlag. See meshdoc, bddoc for details. The data is given by the
%   structure pde which contains function handles f, g_D, g_N, g_R, or d.
%   For general elliptic equations with convection and reaction
%   coefficients, see ellipticpde.
%
%   The function Poisson assembes the matrix equation AD*u = b and solves
%   it by the direct solver (small size <= 2e3) or the multigrid solver
%   (large size > 2e3). The Dirichlet boundary condition is built into the
%   matrix AD and the Neumann boundary condition is build into b.
%
%   The diffusion coefficient d is a scalar function or a column array with
%   the same length as the elem array.
%
%   u = POISSONWG(node,elem,pde,bdFlag,option) specifies the options.
%    - option.dquadorder: quadrature order for diffusion coefficients
%    - option.fquadorder: quadrature order for computing right hand side f
%    - option.solver
%      'direct': the built in direct solver \ (mldivide)
%      'mg':     multigrid-type solvers mg is used.
%      'amg':    algebraic multigrid method is used.
%      'none':   only assemble the matrix equation but not solve.
%   The default setting is to use the direct solver for small size problems
%   and multigrid solvers for large size problems. For more options on the
%   multigrid solver mg, type help mg.
%
%   When only one type of boundary condition is imposed, the input argument
%   bdFlag can be skipped. The boundary condition is implicitly given in
%   the pde structure by specifying g_D or g_N only. See examples below.
%
%   [u,A] = POISSONWG(node,elem,pde,bdFlag) returns also the
%   non-modified stiffness matrix A, which is semi-definite. The kernel of
%   A consists of constant vectors. The matrix A can be used to evulate the
%   bilinear form A(u,v) = u'*A*v, especially the enery norm of a finite
%   element function u by sqrt(u'*A*u).
%
%   [u,A,eqn] = POISSONWG(node,elem,pde,bdFlag) returns also the equation
%   structure eqn, which includes:
%     - eqn.AD:  modified stiffness matrix AD;
%     - eqn.b:   the right hand side.
%   The solution u = AD\b. The output eqn can be used to test other solvers.
%
%   [u,A,eqn,info] = POISSONWG(node,elem,pde,bdFlag) returns also the
%   information on the assembeling and solver, which includes:
%     - info.assembleTime: time to assemble the matrix equation
%     - info.solverTime:   time to solve the matrix equation
%     - info.itStep:       number of iteration steps for the mg solver
%     - info.error:        l2 norm of the residual b - A*u
%     - info.flag:         flag for the mg solver.
%       flag = 0: converge within max iteration
%       flag = 1: iterated maxIt times but did not converge
%       flag = 2: direct solver
%       flag = 3: no solve
%
%   Example
%     squarePoissonWG;
%
%   See also Poisson, squarePoisson, Lshape, crack, mg
%
%   Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('option','var'), option = []; end
NT = size(elem,1);

tic;  % record assembling time
%% Diffusion coefficient
if ~isfield(pde,'d'), pde.d = []; end
if ~isfield(option,'dquadorder'), option.dquadorder = 1; end
if ~isempty(pde.d) && isnumeric(pde.d)
   K = pde.d;                                 % d is an array
end
if ~isempty(pde.d) && ~isnumeric(pde.d)       % d is a function
    [lambda,weight] = quadpts(option.dquadorder);
    nQuad = size(lambda,1);
    K = zeros(NT,1);
    for p = 1:nQuad
		pxy = lambda(p,1)*node(elem(:,1),:) ...
			+ lambda(p,2)*node(elem(:,2),:) ...
			+ lambda(p,3)*node(elem(:,3),:);
        K = K + weight(p)*pde.d(pxy);
   end
end

%% Construct data structure
[elem2edge,edge] = dofedge(elem);
NE = size(edge,1);
Ndof = NT + NE;
elem2dof = NT + elem2edge;

%% Assemble stiffness matrix
A = sparse(Ndof,Ndof);

% compute ct2 = 1/mean(||x-xc||^2)
center = (node(elem(:,1),:) + node(elem(:,2),:) + node(elem(:,3),:))/3;
mid1 = (node(elem(:,2),:) + node(elem(:,3),:))/2;
mid2 = (node(elem(:,3),:) + node(elem(:,1),:))/2;
mid3 = (node(elem(:,1),:) + node(elem(:,2),:))/2;
ct2 = 3./sum((mid1 - center).^2 + (mid2 - center).^2 + (mid3 - center).^2,2);
[Dphi,area] = gradbasis(node,elem);
clear center mid1 mid2 mid3

% Mbb: edge - edge
for i = 1:3
    for j = i:3
        % local to global index map
        ii = double(elem2dof(:,i));
        jj = double(elem2dof(:,j));
        % local stiffness matrix
        Aij = 4*dot(Dphi(:,:,i),Dphi(:,:,j),2).*area + 4/9*ct2.*area;
        if ~isempty(pde.d)
            Aij = K.*Aij;
        end
        if (j==i)
            A = A + sparse(ii,jj,Aij,Ndof,Ndof);
        else
            A = A + sparse([ii,jj],[jj,ii],[Aij; Aij],Ndof,Ndof);
        end
    end
end

% Mob: interior - edge
Aij = -4/3*ct2.*area;
if ~isempty(pde.d)
    Aij = K.*Aij;
end
Mob = sparse([(1:NT)', (1:NT)', (1:NT)'], ...
             double(elem2dof(:)), [Aij, Aij, Aij], Ndof, Ndof);
A = A + Mob + Mob';

% Moo: diagonal of interor
Aij = 4*ct2.*area;
if ~isempty(pde.d)
    Aij = K.*Aij;
end
A =  A + sparse(1:NT, 1:NT, Aij, Ndof, Ndof);

clear K Aij

%% Assemble the right hand side
b = zeros(Ndof,1);
if ~isfield(option,'fquadorder')
    option.fquadorder = 2;   % default order
end
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if ~isempty(pde.f)
    [lambda,weight] = quadpts(option.fquadorder);
	nQuad = size(lambda,1);
    bt = zeros(NT,1);
    for p = 1:nQuad
		% quadrature points in the x-y coordinate
		pxy = lambda(p,1)*node(elem(:,1),:) ...
			+ lambda(p,2)*node(elem(:,2),:) ...
			+ lambda(p,3)*node(elem(:,3),:);
		fp = pde.f(pxy);
        bt = bt + weight(p)*fp;
    end
    bt = bt.*area;
    b(1:NT) = bt;
end
clear pxy bt

%% Set up boundary conditions
if nargin<=3, bdFlag = []; end
[AD,b,u,freeDof,isPureNeumann] = getbdWG(A,b);

%% Record assembling time
assembleTime = toc;
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 2
    fprintf('Time to assemble matrix equation %4.2g s\n',assembleTime);
end

%% Solve the system of linear equations
if isempty(freeDof), return; end
% Set up solver type
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if NE <= 2e3  % Direct solver for small size systems
        option.solver = 'direct';
    else            % MGCG  solver for large size systems
        option.solver = 'mg';
    end
end
solver = option.solver;
% solve
switch solver
    case 'direct'
        tic;
        u(freeDof) = AD(freeDof,freeDof)\b(freeDof);
        residual = norm(b - AD*u);
        info = struct('solverTime',toc,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'none'
        info = struct('solverTime',[],'itStep',0,'err',[],'flag',3,'stopErr',[]);
    case 'mg'
        % eleminate elementwise dof
        option.solver = 'CG';
        option.x0 = u(NT+1:end);
        Aoinv = spdiags(1./diag(AD(1:NT,1:NT)),0,NT,NT);
        Aob = AD(1:NT,NT+1:end);
        Abo = Aob';
        Abb = AD(NT+1:end,NT+1:end);
        Abbm = Abb - Abo*Aoinv*Aob;
        bm = -Abo*Aoinv*b(1:NT) + b(NT+1:end);
        [ub,info] = mg(Abbm,bm,elem,option,edge);
        u(1:NT) = Aoinv*(b(1:NT) - Aob*ub);
        u(NT+1:end) = ub;
%         option.x0 = u;
%         [u,info] = mg(AD,b,elem,option,edge);
    case 'amg'
        option.solver = 'CG';
        [u(freeDof),info] = amg(AD(freeDof,freeDof),b(freeDof),option);
end
% post-process for pure Neumann problem
if isPureNeumann
    uc = sum(u(1:NT).*area);
    u = u - uc;   % normalization for pure Neumann problem
end

%% Output information
eqn = struct('A',AD,'b',b,'edge',edge,'freeDof',freeDof);
info.assembleTime = assembleTime;

%% Compute Du
dudx =  u(elem2dof(:,1)).*Dphi(:,1,1) + u(elem2dof(:,2)).*Dphi(:,1,2) ...
      + u(elem2dof(:,3)).*Dphi(:,1,3);
dudy =  u(elem2dof(:,1)).*Dphi(:,2,1) + u(elem2dof(:,2)).*Dphi(:,2,2) ...
      + u(elem2dof(:,3)).*Dphi(:,2,3);
Du = -2*[dudx, dudy];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions getbdWG
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [AD,b,u,freeDof,isPureNeumann]= getbdWG(A,b)
    %% GETBDCR Boundary conditions for Poisson equation: WG element.

    u =zeros(Ndof,1);
    %% Initial check
    if ~isfield(pde,'g_D'), pde.g_D = []; end
    if ~isfield(pde,'g_N'), pde.g_N = []; end
    if ~isfield(pde,'g_R'), pde.g_R = []; end

    %% Part 1: Modify the matrix for Dirichlet and Robin condition
    % Robin boundary condition
    Robin = [];
    idxR = (bdFlag(:) == 3);      % index of Robin edges in bdFlag
    if any(idxR)
        isRobin = false(NE,1);
        isRobin(elem2edge(idxR)) = true;
        Robin = edge(isRobin,:);  % Robin edges
    end
    if ~isempty(Robin) && ~isempty(pde.g_R) && ~(isnumeric(pde.g_R) && (pde.g_R == 0))
        ve = node(Robin(:,1),:) - node(Robin(:,2),:);
        edgeLength = sqrt(sum(ve.^2,2));
        mid = (node(Robin(:,1),:) + node(Robin(:,2),:))/2;
        ii = NT + find(isRobin);  % for WG: edge dof is after elem dof
        ss = pde.g_R(mid).*edgeLength; % exact for linear g_R
        A = A + sparse(ii,ii,ss,Ndof,Ndof);
    end

    % Find Dirichlet boundary nodes: fixedEdge
    fixedEdge = []; freeEdge = [];
    if ~isempty(bdFlag)              % find boundary edges
        idxD = (bdFlag(:) == 1);     % all Dirichlet edges in bdFlag
        isFixedEdge = false(NE,1);
        isFixedEdge(elem2edge(idxD)) = true;  % index of fixed boundary edges
        fixedEdge = find(isFixedEdge);
        freeEdge = find(~isFixedEdge);
    end
    if isempty(bdFlag) && ~isempty(pde.g_D) && isempty(pde.g_N) && isempty(pde.g_R)
        % no bdFlag, only pde.g_D is given in the input
        s = accumarray(elem2edge(:), 1, [NE 1]);
        fixedEdge = find(s == 1);
        freeEdge = find(s == 2);
    end
    isPureNeumann = false;
    if isempty(fixedEdge) && isempty(Robin)  % pure Neumann boundary condition
        % pde.g_N could be empty which is homogenous Neumann boundary condition
        isPureNeumann = true;
        fixedEdge = 1;
        freeEdge = (2:NE)';    % eliminate the kernel by enforcing u(1) = 0;
    end
    % Modify the matrix
    % Build Dirichlet boundary condition into the matrix AD by enforcing
    % AD(fixedEdge,fixedEdge)=I, AD(fixedEdge,freeEdge)=0, AD(freeEdge,fixedEdge)=0.
    if ~isempty(fixedEdge)
        bdidx = zeros(Ndof,1);
        bdidx(NT + fixedEdge) = 1;
        Tbd = spdiags(bdidx,0,Ndof,Ndof);
        T = spdiags(1-bdidx,0,Ndof,Ndof);
        AD = T*A*T + Tbd;
    else
        AD = A;
    end

    %% Part 2: Find boundary edges and modify the right hand side b
    % Find boundary edges: Neumann
    Neumann = [];
    if ~isempty(bdFlag)  % bdFlag specifies different bd conditions
        idxN = (bdFlag(:) == 2);      % all Neumann edges in bdFlag
        isNeumann = elem2edge(idxN | idxR); % index of Neumann and Robin edges
        % since boundary integral is also needed for Robin edges
        Neumann = edge(isNeumann,:);      % Neumann edges
    end
    if isempty(bdFlag) && (~isempty(pde.g_N) || ~isempty(pde.g_R))
        % no bdFlag, only pde.g_N or pde.g_R is given in the input
        s = accumarray(elem2edge(:), 1, [NE 1]);
        Neumann = edge(s == 1,:);
    end

    % Neumann boundary condition
    if ~isempty(Neumann) && ~isempty(pde.g_N) && ~(isnumeric(pde.g_N) && (pde.g_N == 0))
        if ~isfield(option,'gNquadorder')
            option.gNquadorder = 3;   % default order exact for linear gN
        end
        [lambdagN,weightgN] = quadpts1(option.gNquadorder);
        nQuadgN = size(lambdagN,1);
        ge = zeros(size(Neumann,1),1);
        el = sqrt(sum((node(Neumann(:,1),:) - node(Neumann(:,2),:)).^2,2));
        for pp = 1:nQuadgN
            % quadrature points in the x-y coordinate
            ppxy = lambdagN(pp,1)*node(Neumann(:,1),:) ...
                 + lambdagN(pp,2)*node(Neumann(:,2),:);
            gNp = pde.g_N(ppxy);
            ge = ge+ weightgN(pp)*gNp;
        end
        ge = ge.*el;
        b(NT+isNeumann) = b(NT+isNeumann) + ge;
    end
    % The case with non-empty Neumann edges but g_N=0 or g_N=[] corresponds to
    % the zero flux boundary condition on Neumann edges and no modification of
    % A,u,b is needed.

    % Dirichlet boundary condition
    if ~isPureNeumann && ~isempty(fixedEdge) && ...
       ~isempty(pde.g_D) && ~(isnumeric(pde.g_D) && all(pde.g_D == 0))    % nonzero g_D
        if isnumeric(pde.g_D)  % pde.g_D could be a numerical array
            u(fixedEdge) = pde.g_D(fixedEdge);
        else % pde.g_D is a function handle
            mid = (node(edge(fixedEdge,1),:) + node(edge(fixedEdge,2),:))/2;
            u(NT+fixedEdge) = pde.g_D(mid);
        end
        b = b - A*u;
        b(NT+fixedEdge) = u(NT+fixedEdge);
    end
    % The case with non-empty Dirichlet nodes but g_D=0 or g_D=[] corresponds
    % to the zero Dirichlet boundary condition and no modification of u,b is
    % needed.

    % Pure Neumann boundary condition
    if isPureNeumann
        b = b - mean(b);   % compatilbe condition: sum(b) = 0
        b(1) = 0;
    end

    freeDof = [(1:NT)'; NT+freeEdge];
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end % end of PoissonWG

repo: Aero-matlab


text:
function demo3_minimum_snap_close_form()
clear,clc;

%% condition
waypts = [0,0;
          1,2;
          2,-1;
          4,8;
          5,2]';
v0 = [0,0];
a0 = [0,0];
v1 = [0,0];
a1 = [0,0];
T = 5;
ts = arrangeT(waypts,T);
n_order = 5;

%% trajectory plan
polys_x = minimum_snap_single_axis_close_form(waypts(1,:),ts,n_order,v0(1),a0(1),v1(1),a1(1));
polys_y = minimum_snap_single_axis_close_form(waypts(2,:),ts,n_order,v0(2),a0(2),v1(2),a1(2));

%% result show
figure(1)
plot(waypts(1,:),waypts(2,:),'*r');hold on;
plot(waypts(1,:),waypts(2,:),'b--');
title('minimum snap trajectory');
color = ['grc'];
for i=1:size(polys_x,2)
    tt = ts(i):0.01:ts(i+1);
    xx = polys_vals(polys_x,ts,tt,0);
    yy = polys_vals(polys_y,ts,tt,0);
    plot(xx,yy,color(mod(i,3)+1));
end

end

function polys = minimum_snap_single_axis_close_form(wayp,ts,n_order,v0,a0,v1,a1)
n_coef = n_order+1;
n_poly = length(wayp)-1;
% compute Q
Q_all = [];
for i=1:n_poly
    Q_all = blkdiag(Q_all,computeQ(n_order,3,ts(i),ts(i+1)));
end

% compute Tk   Tk(i,j) = ts(i)^(j-1)
tk = zeros(n_poly+1,n_coef);
for i = 1:n_coef
    tk(:,i) = ts(:).^(i-1);
end

% compute A (n_continuous*2*n_poly) * (n_coef*n_poly)
n_continuous = 3;  % 1:p  2:pv  3:pva  4:pvaj  5:pvajs
A = zeros(n_continuous*2*n_poly,n_coef*n_poly);
for i = 1:n_poly
    for j = 1:n_continuous
        for k = j:n_coef
            if k==j
                t1 = 1;
                t2 = 1;
            else %k>j
                t1 = tk(i,k-j+1);
                t2 = tk(i+1,k-j+1);
            end
            A(n_continuous*2*(i-1)+j,n_coef*(i-1)+k) = prod(k-j+1:k-1)*t1;
            A(n_continuous*2*(i-1)+n_continuous+j,n_coef*(i-1)+k) = prod(k-j+1:k-1)*t2;
        end
    end
end

% compute M
M = zeros(n_poly*2*n_continuous,n_continuous*(n_poly+1));
for i = 1:n_poly*2
    j = floor(i/2)+1;
    rbeg = n_continuous*(i-1)+1;
    cbeg = n_continuous*(j-1)+1;
    M(rbeg:rbeg+n_continuous-1,cbeg:cbeg+n_continuous-1) = eye(n_continuous);
end

% compute C
num_d = n_continuous*(n_poly+1);
C = eye(num_d);
df = [wayp,v0,a0,v1,a1]';% fix all pos(n_poly+1) + start va(2) + end va(2)
fix_idx = [1:3:num_d,2,3,num_d-1,num_d];
free_idx = setdiff(1:num_d,fix_idx);
C = [C(:,fix_idx) C(:,free_idx)];

AiMC = inv(A)*M*C;
R = AiMC'*Q_all*AiMC;

n_fix = length(fix_idx);
Rff = R(1:n_fix,1:n_fix);
Rpp = R(n_fix+1:end,n_fix+1:end);
Rfp = R(1:n_fix,n_fix+1:end);
Rpf = R(n_fix+1:end,1:n_fix);

dp = -inv(Rpp)*Rfp'*df;

p = AiMC*[df;dp];

polys = reshape(p,n_coef,n_poly);

end

author: symao


text:
function [str] = listBiGGModels()
% Reads from url http://bigg.ucsd.edu/api/v2/models
%
% USAGE:
%
%    [str] = listBiGGModels()
%
% OUTPUT:
%    str:    result of url read
%
% .. http://bigg.ucsd.edu/web_api
% .. Get a list of models
str = urlread('http://bigg.ucsd.edu/api/v2/models');

%{
%str = "results_count": 78, "results": [{"organism": "Escherichia coli str. K-12 substr. MG1655", "metabolite_count": 72, "bigg_id": "e_coli_core", "gene_count": 137, "reaction_count": 95}, {"organism": "Homo sapiens", "metabolite_count": 342, "bigg_id": "iAB_RBC_283", "gene_count": 346, "reaction_count": 469}, {"organism": "Escherichia coli str. K-12 substr. MG1655", "metabolite_count": 1668, "bigg_id": "iAF1260", "gene_count": 1261, "reaction_count": 2382}, {"organism": "Methanosarcina barkeri str. Fusaro", "metabolite_count": 628, "bigg_id": "iAF692", "gene_count": 692, "reaction_count": 690}, {"organism": "Geobacter metallireducens GS-15", "metabolite_count": 1109, "bigg_id": "iAF987", "gene_count": 987, "reaction_count": 1285}, {"organism": "Escherichia coli APEC O1", "metabolite_count": 1944, "bigg_id": "iAPECO1_1312", "gene_count": 1313, "reaction_count": 2736}, {"organism": "Homo sapiens", "metabolite_count": 738, "bigg_id": "iAT_PLT_636", "gene_count": 636, "reaction_count": 1008}, {"organism": "Escherichia coli BL21(DE3)", "metabolite_count": 1945, "bigg_id": "iB21_1397", "gene_count": 1337, "reaction_count": 2742}, {"organism": "Escherichia coli BW2952", "metabolite_count": 1949, "bigg_id": "iBWG_1329", "gene_count": 1328, "reaction_count": 2742}, {"organism": "Escherichia coli CFT073", "metabolite_count": 1938, "bigg_id": "ic_1306", "gene_count": 1307, "reaction_count": 2727}, {"organism": "Escherichia coli O127:H6 str. E2348/69", "metabolite_count": 1919, "bigg_id": "iE2348C_1286", "gene_count": 1284, "reaction_count": 2704}, {"organism": "Escherichia coli 042", "metabolite_count": 1926, "bigg_id": "iEC042_1314", "gene_count": 1314, "reaction_count": 2715}, {"organism": "Escherichia coli 55989", "metabolite_count": 1953, "bigg_id": "iEC55989_1330", "gene_count": 1330, "reaction_count": 2757}, {"organism": "Escherichia coli ABU 83972", "metabolite_count": 1944, "bigg_id": "iECABU_c1320", "gene_count": 1320, "reaction_count": 2732}, {"organism": "Escherichia coli B str. REL606", "metabolite_count": 1953, "bigg_id": "iECB_1328", "gene_count": 1329, "reaction_count": 2749}, {"organism": "Escherichia coli 'BL21-Gold(DE3)pLysS AG'", "metabolite_count": 1954, "bigg_id": "iECBD_1354", "gene_count": 1354, "reaction_count": 2749}, {"organism": "Escherichia coli BL21(DE3)", "metabolite_count": 1945, "bigg_id": "iECD_1391", "gene_count": 1333, "reaction_count": 2742}, {"organism": "Escherichia coli str. K-12 substr. DH10B", "metabolite_count": 1947, "bigg_id": "iECDH10B_1368", "gene_count": 1327, "reaction_count": 2743}, {"organism": "Escherichia coli DH1", "metabolite_count": 1949, "bigg_id": "iEcDH1_1363", "gene_count": 1363, "reaction_count": 2751}, {"organism": "Escherichia coli DH1", "metabolite_count": 1950, "bigg_id": "iECDH1ME8569_1439", "gene_count": 1439, "reaction_count": 2756}, {"organism": "Escherichia coli E24377A", "metabolite_count": 1974, "bigg_id": "iEcE24377_1341", "gene_count": 1341, "reaction_count": 2764}, {"organism": "Escherichia coli ED1a", "metabolite_count": 1929, "bigg_id": "iECED1_1282", "gene_count": 1279, "reaction_count": 2707}, {"organism": "Escherichia coli O157:H7 str. EC4115", "metabolite_count": 1918, "bigg_id": "iECH74115_1262", "gene_count": 1262, "reaction_count": 2695}, {"organism": "Escherichia coli HS", "metabolite_count": 1965, "bigg_id": "iEcHS_1320", "gene_count": 1321, "reaction_count": 2754}, {"organism": "Escherichia coli IAI1", "metabolite_count": 1970, "bigg_id": "iECIAI1_1343", "gene_count": 1343, "reaction_count": 2766}, {"organism": "Escherichia coli IAI39", "metabolite_count": 1957, "bigg_id": "iECIAI39_1322", "gene_count": 1321, "reaction_count": 2722}, {"organism": "Escherichia coli NA114", "metabolite_count": 1927, "bigg_id": "iECNA114_1301", "gene_count": 1301, "reaction_count": 2719}, {"organism": "Escherichia coli O103:H2 str. 12009", "metabolite_count": 1958, "bigg_id": "iECO103_1326", "gene_count": 1327, "reaction_count": 2759}, {"organism": "Escherichia coli O111:H- str. 11128", "metabolite_count": 1959, "bigg_id": "iECO111_1330", "gene_count": 1328, "reaction_count": 2761}, {"organism": "Escherichia coli O26:H11 str. 11368", "metabolite_count": 1965, "bigg_id": "iECO26_1355", "gene_count": 1355, "reaction_count": 2781}, {"organism": "Escherichia coli IHE3034", "metabolite_count": 1943, "bigg_id": "iECOK1_1307", "gene_count": 1304, "reaction_count": 2730}, {"organism": "Escherichia coli ATCC 8739", "metabolite_count": 1971, "bigg_id": "iEcolC_1368", "gene_count": 1368, "reaction_count": 2769}, {"organism": "Escherichia coli 536", "metabolite_count": 1943, "bigg_id": "iECP_1309", "gene_count": 1309, "reaction_count": 2740}, {"organism": "Escherichia coli O157:H7 str. Sakai", "metabolite_count": 1923, "bigg_id": "iECs_1301", "gene_count": 1301, "reaction_count": 2721}, {"organism": "Escherichia coli S88", "metabolite_count": 1944, "bigg_id": "iECS88_1305", "gene_count": 1305, "reaction_count": 2730}, {"organism": "Escherichia coli SE11", "metabolite_count": 1957, "bigg_id": "iECSE_1348", "gene_count": 1348, "reaction_count": 2769}, {"organism": "Escherichia coli SE15", "metabolite_count": 1951, "bigg_id": "iECSF_1327", "gene_count": 1327, "reaction_count": 2743}, {"organism": "Escherichia coli SMS-3-5", "metabolite_count": 1949, "bigg_id": "iEcSMS35_1347", "gene_count": 1347, "reaction_count": 2747}, {"organism": "Escherichia coli O157:H7 str. TW14359", "metabolite_count": 1920, "bigg_id": "iECSP_1301", "gene_count": 1299, "reaction_count": 2713}, {"organism": "Escherichia coli UMN026", "metabolite_count": 1935, "bigg_id": "iECUMN_1333", "gene_count": 1332, "reaction_count": 2741}, {"organism": "Escherichia coli W", "metabolite_count": 1975, "bigg_id": "iECW_1372", "gene_count": 1372, "reaction_count": 2783}, {"organism": "Escherichia coli KO11FL", "metabolite_count": 1974, "bigg_id": "iEKO11_1354", "gene_count": 1354, "reaction_count": 2779}, {"organism": "Escherichia coli ETEC H10407", "metabolite_count": 1964, "bigg_id": "iETEC_1333", "gene_count": 1333, "reaction_count": 2757}, {"organism": "Escherichia coli O55:H7 str. CB9615", "metabolite_count": 1919, "bigg_id": "iG2583_1286", "gene_count": 1283, "reaction_count": 2705}, {"organism": "Clostridium ljungdahlii DSM 13528", "metabolite_count": 698, "bigg_id": "iHN637", "gene_count": 637, "reaction_count": 785}, {"organism": "Helicobacter pylori 26695", "metabolite_count": 485, "bigg_id": "iIT341", "gene_count": 339, "reaction_count": 554}, {"organism": "Synechocystis sp. PCC 6803", "metabolite_count": 795, "bigg_id": "iJN678", "gene_count": 622, "reaction_count": 863}, {"organism": "Pseudomonas putida KT2440", "metabolite_count": 909, "bigg_id": "iJN746", "gene_count": 746, "reaction_count": 1056}, {"organism": "Escherichia coli str. K-12 substr. MG1655", "metabolite_count": 1805, "bigg_id": "iJO1366", "gene_count": 1367, "reaction_count": 2583}, {"organism": "Escherichia coli str. K-12 substr. MG1655", "metabolite_count": 761, "bigg_id": "iJR904", "gene_count": 904, "reaction_count": 1075}, {"organism": "Escherichia coli LF82", "metabolite_count": 1940, "bigg_id": "iLF82_1304", "gene_count": 1302, "reaction_count": 2727}, {"organism": "Thermotoga maritima MSB8", "metabolite_count": 570, "bigg_id": "iLJ478", "gene_count": 482, "reaction_count": 652}, {"organism": "Mus musculus", "metabolite_count": 2775, "bigg_id": "iMM1415", "gene_count": 1375, "reaction_count": 3726}, {"organism": "Saccharomyces cerevisiae S288c", "metabolite_count": 1226, "bigg_id": "iMM904", "gene_count": 905, "reaction_count": 1577}, {"organism": "Saccharomyces cerevisiae S288c", "metabolite_count": 1059, "bigg_id": "iND750", "gene_count": 750, "reaction_count": 1266}, {"organism": "Mycobacterium tuberculosis H37Rv", "metabolite_count": 826, "bigg_id": "iNJ661", "gene_count": 661, "reaction_count": 1025}, {"organism": "Escherichia coli O83:H1 str. NRG 857C", "metabolite_count": 1945, "bigg_id": "iNRG857_1313", "gene_count": 1311, "reaction_count": 2736}, {"organism": "Yersinia pestis CO92", "metabolite_count": 1552, "bigg_id": "iPC815", "gene_count": 815, "reaction_count": 1961}, {"organism": null, "metabolite_count": 1706, "bigg_id": "iRC1080", "gene_count": 1086, "reaction_count": 2191}, {"organism": "Shigella flexneri 2a str. 2457T", "metabolite_count": 1914, "bigg_id": "iS_1188", "gene_count": 1188, "reaction_count": 2620}, {"organism": "Staphylococcus aureus subsp. aureus N315", "metabolite_count": 655, "bigg_id": "iSB619", "gene_count": 619, "reaction_count": 743}, {"organism": "Shigella boydii CDC 3083-94", "metabolite_count": 1912, "bigg_id": "iSbBS512_1146", "gene_count": 1147, "reaction_count": 2592}, {"organism": "Shigella boydii Sb227", "metabolite_count": 1910, "bigg_id": "iSBO_1134", "gene_count": 1134, "reaction_count": 2592}, {"organism": "Shigella dysenteriae Sd197", "metabolite_count": 1890, "bigg_id": "iSDY_1059", "gene_count": 1059, "reaction_count": 2540}, {"organism": "Shigella flexneri 2a str. 301", "metabolite_count": 1917, "bigg_id": "iSF_1195", "gene_count": 1195, "reaction_count": 2631}, {"organism": "Shigella flexneri 5 str. 8401", "metabolite_count": 1917, "bigg_id": "iSFV_1184", "gene_count": 1184, "reaction_count": 2622}, {"organism": "Shigella flexneri 2002017", "metabolite_count": 1918, "bigg_id": "iSFxv_1172", "gene_count": 1169, "reaction_count": 2639}, {"organism": "Shigella sonnei Ss046", "metabolite_count": 1938, "bigg_id": "iSSON_1240", "gene_count": 1240, "reaction_count": 2694}, {"organism": "Escherichia coli UM146", "metabolite_count": 1944, "bigg_id": "iUMN146_1321", "gene_count": 1319, "reaction_count": 2736}, {"organism": "Escherichia coli UMNK88", "metabolite_count": 1971, "bigg_id": "iUMNK88_1353", "gene_count": 1353, "reaction_count": 2778}, {"organism": "Escherichia coli UTI89", "metabolite_count": 1942, "bigg_id": "iUTI89_1310", "gene_count": 1310, "reaction_count": 2726}, {"organism": "Escherichia coli W", "metabolite_count": 1975, "bigg_id": "iWFL_1372", "gene_count": 1372, "reaction_count": 2783}, {"organism": "Escherichia coli str. K-12 substr. W3110", "metabolite_count": 1953, "bigg_id": "iY75_1357", "gene_count": 1358, "reaction_count": 2760}, {"organism": "Klebsiella pneumoniae subsp. pneumoniae MGH 78578", "metabolite_count": 1658, "bigg_id": "iYL1228", "gene_count": 1229, "reaction_count": 2262}, {"organism": "Bacillus subtilis subsp. subtilis str. 168", "metabolite_count": 991, "bigg_id": "iYO844", "gene_count": 844, "reaction_count": 1250}, {"organism": "Escherichia coli O157:H7 str. EDL933", "metabolite_count": 1923, "bigg_id": "iZ_1308", "gene_count": 1308, "reaction_count": 2722}, {"organism": "Homo sapiens", "metabolite_count": 2766, "bigg_id": "RECON1", "gene_count": 1905, "reaction_count": 3742}, {"organism": "Salmonella enterica subsp. enterica serovar Typhimurium str. LT2", "metabolite_count": 1802, "bigg_id": "STM_v1_0", "gene_count": 1271, "reaction_count": 2545}]}

%curl 'http://bigg.ucsd.edu/api/v2/models'
%}

author: opencobra

repo: cobratoolbox


text:
function [ joints, s_mat, w_mat ] = InverseKinematics( pose, nsparam, rconf, robot_type )
% INVERSEKINEMATICS
%
% Calculates the inverse kinematics for a KUKA LBR iiwa manipulator.
%
% Input:  pose    - Homogeneous Matrix size(4,4)
%         nsparam - Arm Angle
%         rconf   - Robot Configuration 8-bit number
% Output: joints  - Joint values size(1,7).
%         s_mat   - Shoulder joint matrices As, Bs and Cs
%         w_mat   - Wrist joint matrices Aw, Bw and Cw
%         robot_type - iiwa7 (0) or iiwa14 (~0)

[ arm, elbow, wrist ] = Configuration(rconf);

%Tolerance
tol = 1e-8;

%Robot parameters
%Link length
l = zeros(1,6);
if robot_type == 0
    l = [0.34 0.4 0.4 0.126]; % iiwa7
else
    l = [0.36 0.42 0.4 0.126]; % iiwa14
end

%Denavit-Hartenberg parameters 7 DoF
%DH: [a, alpha,    d, theta]
dh = [0	 -pi/2	 l(1)    0;
      0	  pi/2      0    0;
      0	  pi/2   l(2)	 0;
      0  -pi/2      0    0;
      0  -pi/2   l(3)    0;
      0   pi/2      0    0;
      0      0   l(4)    0];

%Number of joints
nj = size(dh,1);

% Joint values of virtual manipulator
joints = zeros(1,nj);

% Shoulder rotation matrices
s_mat = zeros(3,3,3);

% Wrist rotation matrices
w_mat = zeros(3,3,3);

xend = pose(1:3,4); % end-effector position from base
xs = [0 0 dh(1,3)]'; % shoulder position from base
xwt = [0 0 dh(end,3)]'; % end-effector position from wrist
xw = xend - pose(1:3,1:3)*xwt; % wrist position from base
xsw = xw - xs; % shoulder to wrist vector
usw = unit(xsw);

lbs = l(1);
lse = l(2); % upper arm length (shoulder to elbow)
lew = l(3); % lower arm length (elbow to wrist)

%Check if pose is within arm+forearm reach
assert(norm(xsw) < lse + lew && norm(xsw) > lse - lew, 'Specified pose outside reacheable workspace');

% -- Joint 4 --
% Elbow joint can be directly calculated since it does only depend on the
% robot configuration and the xsw vector
assert(abs((norm(xsw)^2 - lse^2 - lew^2)-(2*lse*lew)) > tol, 'Elbow singularity. Tip at reach limit.');
% Cosine law - According to our robot, joint 4 rotates backwards
joints(4) = elbow * acos((norm(xsw)^2 - lse^2 - lew^2)/(2*lse*lew));

%Added

T34 = dh_calc(dh(4,1),dh(4,2),dh(4,3),joints(4));
R34 = T34(1:3,1:3);

% Shoulder Joints
% First compute the reference joint angles when the arm angle is zero.
[~, R03_o, ~] = ReferencePlane(pose, elbow, robot_type);

skew_usw = skew(usw);
% Following eq. (15), the auxiliary matrixes As Bs and Cs can be calculated
% by substituting eq. (6) into (9).
% R0psi = I3 + sin(psi)*skew_usw + (1-cos(psi))*skew_usw²    (6)
% R03 = R0psi * R03_o                                         (9)
% Substituting (distributive prop.) we get:
% R03 = R03_o*skew_usw*sin(psi) + R03_o*(-skew_usw²)*cos(psi) + R03_o(I3 + skew_usw²)
% R03 =      As       *sin(psi) +        Bs         *cos(psi) +          Cs
As = skew_usw * R03_o;
Bs = -skew_usw^2 * R03_o;
Cs = (usw*usw') * R03_o;

psi = nsparam;
R03 = As*sin(psi) + Bs*cos(psi) + Cs;

% T03 transformation matrix (DH parameters)
%[ cos(j1)*cos(j2)*cos(j3) - sin(j1)*sin(j3), cos(j1)*sin(j2), cos(j3)*sin(j1) + cos(j1)*cos(j2)*sin(j3), 0.4*cos(j1)*sin(j2)]
%[ cos(j1)*sin(j3) + cos(j2)*cos(j3)*sin(j1), sin(j1)*sin(j2), cos(j2)*sin(j1)*sin(j3) - cos(j1)*cos(j3), 0.4*sin(j1)*sin(j2)]
%[                          -cos(j3)*sin(j2),         cos(j2),                          -sin(j2)*sin(j3),  0.4*cos(j2) + 0.34]
%[                                         0,               0,                                         0,                   1]
joints(1) = atan2(arm * R03(2,2), arm * R03(1,2));
joints(2) = arm * acos(R03(3,2));

joints(3) = atan2(arm * -R03(3,3), arm * -R03(3,1));

Aw = R34' * As' * pose(1:3,1:3);
Bw = R34' * Bs' * pose(1:3,1:3);
Cw = R34' * Cs' * pose(1:3,1:3);

R47 = Aw*sin(psi) + Bw*cos(psi) + Cw;

% T47 transformation matrix (DH parameters)
%[ cos(j5)*cos(j6)*cos(j7) - sin(j5)*sin(j7), - cos(j7)*sin(j5) - cos(j5)*cos(j6)*sin(j7), cos(j5)*sin(j6), (63*cos(j5)*sin(j6))/500]
%[ cos(j5)*sin(j7) + cos(j6)*cos(j7)*sin(j5),   cos(j5)*cos(j7) - cos(j6)*sin(j5)*sin(j7), sin(j5)*sin(j6), (63*sin(j5)*sin(j6))/500]
%[                          -cos(j7)*sin(j6),                             sin(j6)*sin(j7),         cos(j6),   (63*cos(j6))/500 + 2/5]
%[                                         0,                                           0,               0,                        1]
joints(5) = atan2(wrist * R47(2,3), wrist * R47(1,3));
joints(6) = wrist * acos(R47(3,3));
joints(7) = atan2(wrist * R47(3,2), wrist * -R47(3,1));

% Grouping Shoulder and Wrist matrices that will be used by the joint
% limit algorithms
s_mat(:,:,1) = As;
s_mat(:,:,2) = Bs;
s_mat(:,:,3) = Cs;
w_mat(:,:,1) = Aw;
w_mat(:,:,2) = Bw;
w_mat(:,:,3) = Cw;

author: xuhuairuogu

repo: V-REP-Simulation-Projects

path: github-repos/MATLAB/xuhuairuogu-V-REP-Simulation-Projects/V-REP-Simulation-Projects-841b944af4ea3a8fb250578d36434515f577f411/admittance control(adapted from an example in the book--A Systematic Approach to Learning Robot Programming with ROS)/iiwa14_kinematics/InverseKinematics.m


text:
function resultingMatrix = blockDiagonalMatrix(param,varargin)
% blockDiagonalMatrix provides an interface for creating block diagonal matrices

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.
  prev_class  = class(varargin{1});

  % Check if the matrices provided for building the block diagonal are all variables
  % of the same class
  for i=1:nargin-1
    if(isa(varargin{i},prev_class))
      prev_class=class(varargin{i});
    else
      error('All the inputs arguments must be of the same class')
    end
  end

  % When running the model using symbolic variables, the creation of the block diagonal has to be carried out
  % differently if Octave is used
  if(isa(varargin{1},'casadi.SX') || isa(varargin{1},'casadi.MX'))
    if(param.isMatlab)
      resultingMatrix = blkdiag(varargin{:});
    else
      size_tot                                                          = [0,0];
      for i=1:nargin-1
        size_tot = size_tot + size(varargin{i});
      end
      resultingMatrix                                                             = SX.zeros(size_tot);
      start_position = [1,1];
      % Create the block diagonal
      for i=1:nargin-1
        resultingMatrix(start_position(1):start_position(1)+size(varargin{i},1)-1,start_position(2):start_position(2)+size(varargin{i},2)-1) = varargin{i};
        start_position(1) = start_position(1) + size(varargin{i},1);
        start_position(2) = start_position(2) + size(varargin{i},2);
      end
    end
  else
    resultingMatrix = blkdiag(varargin{:});
  end

end

author: lionsimbatoolbox

repo: LIONSIMBA


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

% attempt to find the markers using simple vision methods
i2=imorph(i,ones(3,3),'max');	% eliminate the markers
iff=flatfield(i2);		% and compute the flat field
idisp(iff-i)
pause
[k,s]=kmeans(iff-i,3);		% now k-means segmentation
idisp(s)
pause
[l,nblobs]=ilabel(s==3);	% and assume brightest points are markers
idisp(l)
colormap(hot(nblobs));
pause
f=iblobs(l)			% now extract binary image features
hold on
for i=1:12
	if f(i).area < 100,
		plot(f(i).xc,f(i).yc,'*w')
	end
end

% now enter the refinement phase

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [dat] = read_ricoh_data(filename, hdr, begsample, endsample, chanindx)

%function [dat] = read_ricoh_data(filename, hdr, begsample, endsample, chanindx)
%
%% READ_RICOH_DATA reads continuous or averaged MEG data
%% generated by the RICOH MEG system and software,
%% and allows the data to be used in FieldTrip.
%%
%% Use as
%%   [dat] = read_ricoh_data(filename, hdr, begsample, endsample, chanindx)
%%
%% This is a wrapper function around the function getRData
%%
%% See also READ_RICOH_HEADER, READ_RICOH_EVENT

% ensure that the required toolbox is on the path
ft_hastoolbox('ricoh_meg_reader', 1);

hdr = hdr.orig; % use the original Ricoh header, not the FieldTrip header

% default is to select all channels
if nargin<5
  chanindx = 1:hdr.channel_count;
end

handles = definehandles;

switch hdr.acq_type
  case handles.AcqTypeEvokedAve
    % dat is returned as double
    start_sample  = begsample - 1; % samples start at 0
    sample_length = endsample - begsample + 1;
    epoch_count   = 1;
    start_epoch   = 0;
    dat = getRData(filename, start_sample, sample_length);

  case handles.AcqTypeContinuousRaw
    % dat is returned as double
    start_sample  = begsample - 1; % samples start at 0
    sample_length = endsample - begsample + 1;
    epoch_count   = 1;
    start_epoch   = 0;
    dat = getRData(filename, start_sample, sample_length);

    % Unlike Yokogawa system, "AcqTypeEvokedRaw" is not supported for Ricoh system.

  otherwise
    ft_error('unknown data type');
end

if size(dat,1)~=hdr.channel_count
  ft_error('could not read all channels');
elseif size(dat,2)~=(endsample-begsample+1)
  ft_error('could not read all samples');
end

% select only the desired channels
dat = dat(chanindx,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this defines some usefull constants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function handles = definehandles
handles.output = [];
handles.sqd_load_flag = false;
handles.mri_load_flag = false;
handles.NullChannel         = 0;
handles.MagnetoMeter        = 1;
handles.AxialGradioMeter    = 2;
handles.PlannerGradioMeter  = 3;
handles.RefferenceChannelMark = hex2dec('0100');
handles.RefferenceMagnetoMeter       = bitor( handles.RefferenceChannelMark, handles.MagnetoMeter );
handles.RefferenceAxialGradioMeter   = bitor( handles.RefferenceChannelMark, handles.AxialGradioMeter );
handles.RefferencePlannerGradioMeter = bitor( handles.RefferenceChannelMark, handles.PlannerGradioMeter );
handles.TriggerChannel      = -1;
handles.EegChannel          = -2;
handles.EcgChannel          = -3;
handles.EtcChannel          = -4;
handles.NonMegChannelNameLength = 32;
handles.DefaultMagnetometerSize       = (4.0/1000.0);       % Square of 4.0mm in length
handles.DefaultAxialGradioMeterSize   = (15.5/1000.0);      % Circle of 15.5mm in diameter
handles.DefaultPlannerGradioMeterSize = (12.0/1000.0);      % Square of 12.0mm in length
handles.AcqTypeContinuousRaw = 1;
handles.AcqTypeEvokedAve     = 2;
handles.AcqTypeEvokedRaw     = 3;
handles.sqd = [];
handles.sqd.selected_start  = [];
handles.sqd.selected_end    = [];
handles.sqd.axialgradiometer_ch_no      = [];
handles.sqd.axialgradiometer_ch_info    = [];
handles.sqd.axialgradiometer_data       = [];
handles.sqd.plannergradiometer_ch_no    = [];
handles.sqd.plannergradiometer_ch_info  = [];
handles.sqd.plannergradiometer_data     = [];
handles.sqd.eegchannel_ch_no   = [];
handles.sqd.eegchannel_data    = [];
handles.sqd.nullchannel_ch_no   = [];
handles.sqd.nullchannel_data    = [];
handles.sqd.selected_time       = [];
handles.sqd.sample_rate         = [];
handles.sqd.sample_count        = [];
handles.sqd.pretrigger_length   = [];
handles.sqd.matching_info   = [];
handles.sqd.source_info     = [];
handles.sqd.mri_info        = [];
handles.mri                 = [];

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function CPT2 = mk_named_CPT(family_names, names, dag, CPT1)
% MK_NAMED_CPT Permute the dimensions of a CPT so they agree with the internal numbering convention
% CPT2 = mk_named_CPT(family_names, names, dag, CPT1)
%
% This is best explained by example.
% Consider the following directed acyclic graph
%
%     C
%   /   \
%  R     S
%   \   /
%     W
%
% where all arcs point down.
% When we create the CPT for node W, we consider S as its first parent, and R as its
% second, and hence write
%
%      S R W
% CPT1(1,1,:) = [1.0 0.0];
% CPT1(2,1,:) = [0.2 0.8];  % P(W=1 | R=1, S=2) = 0.2
% CPT1(1,2,:) = [0.1 0.9];
% CPT1(2,2,:) = [0.01 0.99];
%
% However, when we create the dag using mk_adj_mat, the nodes get topologically sorted,
% and by chance, node R preceeds node S in this ordering.
% Hence we should have written
%
%      R S W
% CPT2(1,1,:) = [1.0 0.0];
% CPT2(2,1,:) = [0.1 0.9];
% CPT2(1,2,:) = [0.2 0.8]; % P(W=1 | R=1, S=2) = 0.2
% CPT2(2,2,:) = [0.01 0.99];
%
% Since we do not know the order of the nodes in advance, we can write
%   CPT2 = mk_named_CPT({'S', 'R', 'W'}, names, dag, CPT1)
% where 'S', 'R', 'W' are the order of the dimensions we assumed (the child node must be last in this list),
% and names{i} is the name of the i'th node.

n = length(family_names);
family_nums = zeros(1,n);
for i=1:n
  family_nums(i) = stringmatch(family_names{i}, names); % was strmatch
end

fam = family(dag, family_nums(end));
perm = zeros(1,n);
for i=1:n
  %  perm(i) = find(family_nums(i) == fam);
  perm(i) = find(fam(i) == family_nums);
end

CPT2 = permute(CPT1, perm);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [level,bw] = thresh_tool(im,cmap,defaultLevel) %mainfunction
%THRESH_TOOL  Interactively select intensity level for image thresholding.
%   THRESH_TOOL launches a GUI (graphical user interface) for thresholding
%   an intensity input image, IM. IM is displayed in the top left corner. A
%   colorbar and IM's histogram are displayed on the bottom. A line on the
%   histogram indicates the current threshold level. A binary image is
%   displayed in the top right based on the selected level.  To change the
%   level, click and drag the line. The output image updates automatically.
%
%   There are two ways to use this tool.
%
%Mode 1 - nonblocking behavior:
%   THRESH_TOOL(IM) launches GUI tool.  You can continue using the MATLAB
%   Desktop.  Since no results are needed, the function does not block
%   execution of other commands.
%
%   THRESH_TOOL(IM,CMAP) allows the user to specify the colormap, CMAP.  If
%   not specified, the default colormap is used.
%
%   THRESH_TOOL(IM,CMAP,DEFAULTLEVEL) allows the user to specify the
%   default threshold level. If not specified, DEFAULTLEVEL is determined
%   by GRAYTHRESH. Valid values for DEFAULTLEVEL must be consistent with
%   the data type of IM for integer intensity images: uint8 [0,255], uint16
%   [0,65535], int16 [-32768,32767].
%
%   Example
%       x = imread('rice.png');
%       thresh_tool(x)          %no return value, so MATLAB keeps running
%
%Mode 2 - blocking behavior:
%   LEVEL = THRESH_TOOL(...) returns the user selected level, LEVEL, and
%   MATLAB waits for the result before proceeding.  This blocking behavior
%   mode allows the tool to be inserted into an image processing algorithm
%   to support an automated workflow.
%
%   [LEVEL,BW] = THRESH_TOOL(...) also returns the thresholded binary
%   output image, BW.
%
%   Example
%       x = imread('rice.png');
%       lev = thresh_tool(x)    %MATLAB waits for GUI tool to finish
%
%See also COLORMAP, GRAYTHRESH, IM2BW.

%defensive programming
error(nargchk(1,3,nargin))
error(nargoutchk(0,2,nargout))

%validate defaultLevel within range
if nargin>2 %uer specified DEFAULTLEVEL
  dataType = class(im);
  switch dataType
    case 'uint8','uint16','int16'
      if defaultLevel<intmin(dataType) | defaultLevel>intmax(dataType)
        error(['Specified DEFAULTLEVEL outside class range for ' dataType])
      elseif defaultLevel<min(im(:)) | defaultLevel>max(im(:))
        error('Specified DEFAULTLEVEL outside data range for IM')
      end
    case 'double','single'
      %okay, do nothing
    otherwise
      error(['Unsupport image type ' dataType])
  end %switch
end

max_colors=1000;    %practical limit

%calculate bins centers
color_range = double(limits(im));
if isa(im,'uint8') %special case [0 255]
  color_range = [0 255];
  num_colors = 256;
  di = 1;
elseif isinteger(im)
  %try direct indices first
  num_colors = diff(color_range)+1;
  if num_colors<max_colors %okay
    di = 1;                                 %inherent bins
  else %too many levels
    num_colors = max_colors;                %practical limit
    di = diff(color_range)/(num_colors-1);
  end
else %noninteger
  %try infering discrete resolution first (intensities often quantized)
  di = min(diff(sort(unique(im(:)))));
  num_colors = round(diff(color_range)/di)+1;
  if num_colors>max_colors %too many levels
    num_colors = max_colors;                %practical limit
    di = diff(color_range)/(num_colors-1);
  end
end
bin_ctrs = [color_range(1):di:color_range(2)];
FmtSpec = ['%.' num2str(ceil(-log10(di))) 'f'];

%new figure - interactive GUI tool for level segmenting
h_fig = figure;
set(h_fig,'ToolBar','Figure')
if nargin>1 && isstr(cmap) && strmatch(lower(cmap),'gray')
  full_map = gray(num_colors);
elseif nargin>1 && isnumeric(cmap) && length(size(cmap))==2 && size(cmap,2)==3
  full_map = cmap;
else
  full_map = jet(num_colors);
end
setappdata(h_fig,'im',im)
setappdata(h_fig,'FmtSpec',FmtSpec)

%top left - input image
h_ax1 = axes('unit','norm','pos',[0.05 0.35 0.4 0.60]);
rgb = im2rgb(im,full_map);
  function rgb = im2rgb(im,full_map); %nested
    %coerce intensities into gray range [0,1]
    gray = imadjust(im,[],[0 1]);
    %generate indexed image
    num_colors = size(full_map,1);
    ind = gray2ind(gray,num_colors);
    %convert indexed image to RGB
    rgb = ind2rgb(ind,full_map);
  end %im2rgb

image(rgb), axis image
%subimage(im,full_map)
axis off, title('Input Image')

%top right - segmented (eventually)
h_ax2 = axes('unit','norm','pos',[0.55 0.35 0.4 0.60]);
axis off
setappdata(h_fig,'h_ax2',h_ax2)

%next to bottom - intensity distribution
h_hist = axes('unit','norm','pos',[0.05 0.1 0.9 0.2]);
n = hist(double(im(:)),bin_ctrs);
bar(bin_ctrs,n)
axis([color_range limits(n(2:end-1))]) %ignore saturated end scaling
set(h_hist,'xtick',[],'ytick',[])
title('Intensity Distribution')

%very bottom - colorbar
h_cbar = axes('unit','norm','pos',[0.05 0.05 0.9 0.05],'tag','thresh_tool_cbar');
subimage(color_range,[0.5 1.5],1:num_colors,full_map)
set(h_cbar,'ytick',[],'xlim',color_range)
axis normal

v=version;
if str2num(v(1:3))>=7
  %link top axes (pan & zoom)
  linkaxes([h_ax1 h_ax2])
  %link bottom axes (X only - pan & zoom)
  linkaxes([h_hist h_cbar],'x')
end

%colorbar tick locations
set(h_cbar,'xtick',color_range)

%threshold level - initial guess (graythresh)
if nargin>2 %user specified default level
  my_level = defaultLevel;
else %graythresh default
  lo = double(color_range(1));
  hi = double(color_range(2));
  norm_im = (double(im)-lo)/(hi-lo);
  norm_level = graythresh(norm_im); %GRAYTHRESH assumes DOUBLE range [0,1]
  my_level = norm_level*(hi-lo)+lo;
end

%display level as vertical line
axes(h_hist)
h_lev = vline(my_level,'-');
set(h_lev,'LineWidth',2,'color',0.5*[1 1 1],'UserData',my_level)
setappdata(h_fig,'h_lev',h_lev)

%attach draggable behavior for user to change level
move_vline(h_lev,@update_plot);

axes(h_cbar)
y_lim = get(h_cbar,'ylim');

% PLACE TEXT LOCATION ON COLORBAR (Laurens)
%h_text = text(my_level,mean(y_lim),num2str(round(my_level)));
h_text = text(my_level,mean(y_lim),'dummy','HorizontalAlignment','Center');
if nargin<2
  text_color = 0.5*[1 1 1];
else
  text_color = 'm';
end
set(h_text,'FontWeight','Bold','color',text_color,'Tag','cbar_text')
movex_text(h_text,my_level)
%%%%%%%%%%%%%%%%%%%%%%%%

%segmented image
bw = im>my_level;
axes(h_ax2)
hold on
subimage(bw), axis off, axis ij
hold off
title('Segmented')

update_plot

%add reset button (resort to initial guess)
h_reset = uicontrol('unit','norm','pos',[0.0 0.95 .1 .05]);
set(h_reset,'string','Reset','callback',@ResetOriginalLevel)

if nargout>0 %return result(s)
  h_done = uicontrol('unit','norm','pos',[0.9 0.95 0.1 0.05]);
  set(h_done,'string','Done','callback','delete(gcbo)') %better
  %inspect(h_fig)
  set(h_fig,'WindowStyle','modal')
  waitfor(h_done)
  if ishandle(h_fig)
    h_lev = getappdata(gcf,'h_lev');
    level = mean(get(h_lev,'xdata'));
    if nargout>1
      h_im2 = findobj(h_ax2,'type','image');
      bw = logical(rgb2gray(get(h_im2,'cdata')));
    end
    delete(h_fig)
  else
    warning('THRESHTOOL:UserAborted','User Aborted - no return value')
    level = [];
  end
end

end %thresh_tool (mainfunction)

function ResetOriginalLevel(hObject,varargin) %subfunction
h_lev = getappdata(gcf,'h_lev');
init_level = get(h_lev,'UserData');
set(h_lev,'XData',init_level*[1 1])
text_obj = findobj('Type','Text','Tag','cbar_text');
movex_text(text_obj,init_level)
update_plot
end %ResetOriginalLevel (subfunction)

function update_plot %subfunction
im = getappdata(gcf,'im');
h_lev = getappdata(gcf,'h_lev');
my_level = mean(get(h_lev,'xdata'));
h_ax2 = getappdata(gcf,'h_ax2');
h_im2 = findobj(h_ax2,'type','image');
%segmented image
bw = (im>my_level);
rgb_version = repmat(double(bw),[1 1 3]);
set(h_im2,'cdata',rgb_version)
end %update_plot (subfunction)

%function rgbsubimage(im,map), error('DISABLED')

%----------------------------------------------------------------------
function move_vline(handle,DoneFcn) %subfunction
%MOVE_VLINE implements horizontal movement of line.
%
%  Example:
%    plot(sin(0:0.1:pi))
%    h=vline(1);
%    move_vline(h)
%
%Note: This tools strictly requires MOVEX_TEXT, and isn't much good
%      without VLINE by Brandon Kuczenski, available at MATLAB Central.
%<http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=1039&objectType=file>

% This seems to lock the axes position
set(gcf,'Nextplot','Replace')
set(gcf,'DoubleBuffer','on')

h_ax=get(handle,'parent');
h_fig=get(h_ax,'parent');
setappdata(h_fig,'h_vline',handle)
if nargin<2, DoneFcn=[]; end
setappdata(h_fig,'DoneFcn',DoneFcn)
set(handle,'ButtonDownFcn',@DownFcn)

  function DownFcn(hObject,eventdata,varargin) %Nested--%
    set(gcf,'WindowButtonMotionFcn',@MoveFcn)           %
    set(gcf,'WindowButtonUpFcn',@UpFcn)                 %
  end %DownFcn------------------------------------------%

  function UpFcn(hObject,eventdata,varargin) %Nested----%
    set(gcf,'WindowButtonMotionFcn',[])                 %
    DoneFcn=getappdata(hObject,'DoneFcn');              %
    if isstr(DoneFcn)                                   %
      eval(DoneFcn)                                     %
    elseif isa(DoneFcn,'function_handle')               %
      feval(DoneFcn)                                    %
    end                                                 %
  end %UpFcn--------------------------------------------%

  function MoveFcn(hObject,eventdata,varargin) %Nested------%
    h_vline=getappdata(hObject,'h_vline');                  %
    h_ax=get(h_vline,'parent');                             %
    cp = get(h_ax,'CurrentPoint');                          %
    xpos = cp(1);                                           %
    x_range=get(h_ax,'xlim');                               %
    if xpos<x_range(1), xpos=x_range(1); end                %
    if xpos>x_range(2), xpos=x_range(2); end                %
    XData = get(h_vline,'XData');                           %
    XData(:)=xpos;                                          %
    set(h_vline,'xdata',XData)                              %
    %update text                                            %
    text_obj = findobj('Type','Text','Tag','cbar_text');    %
    movex_text(text_obj,xpos)                               %
  end %MoveFcn----------------------------------------------%

end %move_vline(subfunction)

%----------------------------------------------------------------------
function [x,y] = limits(a) %subfunction
% LIMITS returns min & max values of matrix; else scalar value.
%
%   [lo,hi]=LIMITS(a) returns LOw and HIgh values respectively.
%
%   lim=LIMITS(a) returns 1x2 result, where lim = [lo hi] values

if nargin~=1 | nargout>2 %bogus syntax
  error('usage: [lo,hi]=limits(a)')
end

siz=size(a);

if prod(siz)==1 %scalar
  result=a;                         % value
else %matrix
  result=[min(a(:)) max(a(:))];     % limits
end

if nargout==1 %composite result
  x=result;                         % 1x2 vector
elseif nargout==2 %separate results
  x=result(1);                      % two scalars
  y=result(2);
else %no result
  ans=result                        % display answer
end

end %limits (subfunction)

%----------------------------------------------------------------------
function movex_text(h_txt,x_pos) %subfunction
FmtSpec=getappdata(get(get(h_txt,'parent'),'parent'),'FmtSpec');
msg=sprintf(FmtSpec,x_pos);
pos=get(h_txt,'position');
pos(1)=x_pos;
set(h_txt,'Position',pos,'String',msg)
end %movex_text

%--------------------------------------------------------------------------------------------------------------
function hhh=vline(x,in1,in2) %subfunction
% function h=vline(x, linetype, label)
%
% Draws a vertical line on the current axes at the location specified by 'x'.  Optional arguments are
% 'linetype' (default is 'r:') and 'label', which applies a text label to the graph near the line.  The
% label appears in the same color as the line.
%
% The line is held on the current axes, and after plotting the line, the function returns the axes to
% its prior hold state.
%
% The HandleVisibility property of the line object is set to "off", so not only does it not appear on
% legends, but it is not findable by using findobj.  Specifying an output argument causes the function to
% return a handle to the line, so it can be manipulated or deleted.  Also, the HandleVisibility can be
% overridden by setting the root's ShowHiddenHandles property to on.
%
% h = vline(42,'g','The Answer')
%
% returns a handle to a green vertical line on the current axes at x=42, and creates a text object on
% the current axes, close to the line, which reads "The Answer".
%
% vline also supports vector inputs to draw multiple lines at once.  For example,
%
% vline([4 8 12],{'g','r','b'},{'l1','lab2','LABELC'})
%
% draws three lines with the appropriate labels and colors.
%
% By Brandon Kuczenski for Kensington Labs.
% brandon_kuczenski@kensingtonlabs.com
% 8 November 2001

% Downloaded 8/7/03 from MATLAB Central
% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=1039&objectType=file

if length(x)>1  % vector input
    for I=1:length(x)
        switch nargin
        case 1
            linetype='r:';
            label='';
        case 2
            if ~iscell(in1)
                in1={in1};
            end
            if I>length(in1)
                linetype=in1{end};
            else
                linetype=in1{I};
            end
            label='';
        case 3
            if ~iscell(in1)
                in1={in1};
            end
            if ~iscell(in2)
                in2={in2};
            end
            if I>length(in1)
                linetype=in1{end};
            else
                linetype=in1{I};
            end
            if I>length(in2)
                label=in2{end};
            else
                label=in2{I};
            end
        end
        h(I)=vline(x(I),linetype,label);
    end
else
    switch nargin
    case 1
        linetype='r:';
        label='';
    case 2
        linetype=in1;
        label='';
    case 3
        linetype=in1;
        label=in2;
    end

    g=ishold(gca);
    hold on

    y=get(gca,'ylim');
    h=plot([x x],y,linetype);
    if length(label)
        xx=get(gca,'xlim');
        xrange=xx(2)-xx(1);
        xunit=(x-xx(1))/xrange;
        if xunit<0.8
            text(x+0.01*xrange,y(1)+0.1*(y(2)-y(1)),label,'color',get(h,'color'))
        else
            text(x-.05*xrange,y(1)+0.1*(y(2)-y(1)),label,'color',get(h,'color'))
        end
    end

    if g==0
    hold off
    end
    set(h,'tag','vline','handlevisibility','off')
end % else

if nargout
    hhh=h;
end

end %vline (subfunction)

author: Sable

repo: mcbench-benchmarks


text:
function [k, sk] = whitefixedKernCompute(kern, x, x2)

% WHITEFIXEDKERNCOMPUTE Compute the WHITEFIXED kernel given the parameters and X.
% FORMAT
% DESC computes the kernel parameters for the fixed parameter white noise
% kernel given inputs associated with rows and columns.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : the input matrix associated with the rows of the kernel.
% ARG x2 : the inpute matrix associated with the columns of the kernel.
% RETURN k : the kernel matrix computed at the given points.
%
% FORMAT
% DESC computes the kernel matrix for the fixed parameter white noise
% kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : the kernel matrix computed at the given points.
%
% SEEALSO : whitefixedKernParamInit, kernCompute, kernCreate, whitefixedKernDiagCompute
%
% COPYRIGHT : Nathaniel J. King, 2006
%
% MODIFICATIONS : Neil D. Lawrence, 2009

% KERN

if nargin < 3
  k = whiteKernCompute(kern, x);
else
  k = whiteKernCompute(kern, x, x2);
end
if nargout > 1
  sk = k;
end


text:
function dat = proc_BeerLambert(dat, varargin)
% PROC_BEERLAMBERT - apply the BeerLambert function, which is useful to analyze
% NIRS data:  relative concentration is calculated as a function of total photon path length
%
% Synopsis:
%   DAT = nirs_LB(DAT, 'Property1',Value1, ...)
%
% Arguments:
%   DAT: data struct with NIRS cnt data split in the x field
%
%   OPT - struct or property/value list of optional properties:
%   'Citation'   - if set, the epsilon (extinction values) is taken from
%                  the specified citation number (see GetExtinctions.m for
%                  details). (default 1 = Gratzer et al)
%   'Epsilon'    - sets extinction coefficients manually (wl1 wl2 vs deoxy
%                  and oxy). In this case, citation is ignored. State in
%                  millimol/liter(?).
%   'Opdist'     - optode (source-detector) distance in cm (default 2.5)
%   'Ival'       - either 'all' (default) or a vector [start end] in samples
%                  specifying the baseline for the LB transformation
%   'DPF'        - differential pathlength factor: probabilistic average
%                  distance travelled by photons, default [5.98 7.15]
%
% Returns:
%   DAT - updated data struct with oxy/deoxy fields specifying absorption
%         values in mmol/l.
%
%
% See also: nirs_* nirsfile_* GetExtinctions
%
% Note: Based on the nirX Nilab toolbox functions u_LBG and u_popLBG.

% matthias.treder@tu-berlin.de 2011, mail@aewald.net 2013
% AE: changed channels label from wl1&wl2 to oxa and deoxy.
% Markus Wenzel 2013 (adapted it to the new toolbox)
% Jan Mehnert February 2014 (ready for public BBCI toolbox) (jan@mehnert.org)
% ToDo: possibility to define individual baseline; online functionability

props={'Citation'   1             'INT'
        'Opdist'    2.5           'DOUBLE'
        'Ival'      'all'         'CHAR|DOUBLE'
        'DPF'       [5.98 7.15]   'DOUBLE'
        'Epsilon'   []            'DOUBLE'
        'Verbose'   0             'BOOL'};

if nargin==0,
    dat= props; return
end

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);

opt_checkProplist(opt, props);
misc_checkType(dat, 'STRUCT');

if strcmp(opt.Ival,'all')
  opt.Ival = [1 size(dat.x,1)];
end

s1=size(dat.x,1);
s2=size(dat.x,2)/2;

wl1 = dat.x(:,1:end/2) + eps;
wl2 = dat.x(:,end/2+1:end) + eps;

%% Get epsilon
if isempty(opt.Epsilon)
  if ~isfield(dat,'wavelengths')
    error('Wavelengths should be given in the .wavelengths field.')
  end
  [ext,nfo] = procutil_getExtinctions(dat.wavelengths,opt.Citation);
  if opt.Verbose, fprintf('Citation: %s\n',nfo), end
  epsilon = ext(:,1:2);

  % Divide by 1000 to obtain the required unit
  epsilon = epsilon/1000;

else
  epsilon = opt.Epsilon;
end

%% Arrange epsilon so that higher wavelength is on top
[mw,idx] = max(dat.wavelengths);
if idx==2 % higher wavelength is on bottom
  epsilon = flipud(epsilon);
  if opt.Verbose, fprintf('Epsilon matrix was rearranged so that higher WL is on top\n'), end
end

%% Apply Lambert-Beer law
Att_highWL= real(-log10( wl2 ./ ...
    ( repmat(mean(wl2(opt.Ival(1):opt.Ival(2),:),1), [s1,1]))   ));

Att_lowWL= real(-log10( wl1./ ...
    ( repmat(mean(wl1(opt.Ival(1):opt.Ival(2),:),1), [s1,1]))   ));

A=[];
A(:,1)=reshape(Att_highWL,s1*s2,1);
A(:,2)=reshape(Att_lowWL,s1*s2,1);

%----------------------------------
%       3.cc
%----------------------------------
% e=...looks like this
%               oxy-Hb         deoxy-Hb
% higherWL: 830 | e: 0.974       0.693
% lowerWL : 690 | e: 0.35         2.1

e= epsilon/10;

e2=   e.* [opt.DPF' opt.DPF']  .*  opt.Opdist;
c= ( inv(e2)*A'  )';

dat.x = reshape(c(:,1),s1,s2); %in mmol/l
dat.x = [dat.x reshape(c(:,2),s1,s2)]; %in mmol/l

%% Change Channel labels wl1 & wl2 to 'oxy' and 'deoxy'

dat.clab = strrep(dat.clab, 'highWL', 'oxy');
dat.clab = strrep(dat.clab, 'lowWL', 'deoxy');

dat.signal = 'NIRS (oxy, deoxy)';

dat.yUnit = 'mmol/l';

author: bbci


text:
function [data] = ft_determine_coordsys(data, varargin)

% FT_DETERMINE_COORDSYS plots a geometrical object, allowing you to perform
% a visual check on the coordinatesystem, the units and on the anatomical
% labels for the coordinate system axes.
%
% Use as
%   [dataout] = ft_determine_coordsys(datain, ...)
% where the input data structure can be either
%  - an anatomical MRI
%  - a cortical or head surface mesh
%  - an electrode, gradiometer or optode definition
%  - a volume conduction model of the head
% or most other FieldTrip structures that represent geometrical information.
%
% Additional optional input arguments should be specified as key-value pairs
% and can include
%   interactive  = string, 'yes' or 'no' (default = 'yes')
%   axisscale    = scaling factor for the reference axes and sphere (default = 1)
%   clim         = lower and upper anatomical MRI limits (default = [0 1])
%
% This function wil pop up a figure that allows you to check whether the
% alignment of the object relative to the coordinate system axes is correct
% and what the anatomical labels of the coordinate system axes are. You
% should switch on the 3D rotation option in the figure panel to rotate and
% see the figure from all angles. To change the anatomical labels of the
% coordinate system, you should press the corresponding keyboard button.
%
% Recognized and supported coordinate systems are 'ctf', 'bti', '4d', 'yokogawa',
% 'eeglab', 'neuromag', 'itab', 'acpc', 'spm', 'mni', 'fsaverage', 'tal', 'scanras',
% 'scanlps', 'dicom'.
%
% Furthermore, supported coordinate systems that do not specify the origin are 'ras',
% 'als', 'lps', etc. See https://www.fieldtriptoolbox.org/faq/coordsys for more
% details.
%
% See also FT_CONVERT_COORDSYS, FT_DETERMINE_UNITS, FT_CONVERT_UNITS, FT_PLOT_AXES, FT_PLOT_XXX

% Copyright (C) 2015-2021, Jan-Mathijs Schoffelen
%
% 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$

dointeractive = ft_getopt(varargin, 'interactive', 'yes');
axisscale     = ft_getopt(varargin, 'axisscale', 1);  % this is used to scale the axmax and rbol
clim          = ft_getopt(varargin, 'clim', [0 1]);   % this is used to scale the orthoplot
fontsize      = ft_getopt(varargin, 'fontsize');      % this is passed to ft_plot_axes

data  = ft_checkdata(data, 'hasunit', 'yes');
dtype = ft_datatype(data);

% the high-level data structures are detected with ft_datatype, but there are
% also some low-level data structures that need to be supproted here
if strcmp(dtype, 'unknown')
  if isfield(data, 'fid') || (isfield(data, 'tri') && isfield(data, 'pos'))
    dtype = 'headshape';
  elseif isfield(data, 'hex') && isfield(data, 'pos')
    dtype = 'mesh';
  elseif isfield(data, 'tet') && isfield(data, 'pos')
    dtype = 'mesh';
  elseif ~strcmp(ft_headmodeltype(data), 'unknown')
    dtype = 'headmodel';
  elseif ~strcmp(ft_senstype(data), 'unknown')
    dtype = 'sens';
  end
elseif strcmp(dtype, 'mesh+label')
  % we don't care about the labels here
  dtype = 'mesh';
end

if isfield(data, 'coordsys') && ~isempty(data.coordsys)
  % ensure that it is in lower case
  data.coordsys = lower(data.coordsys);
  % print the interpretation of the coordinate system
  [labelx, labely, labelz] = coordsys2label(data.coordsys, 2, 0);
  fprintf('The positive x-axis is pointing towards %s\n', labelx);
  fprintf('The positive y-axis is pointing towards %s\n', labely);
  fprintf('The positive z-axis is pointing towards %s\n', labelz);
end

% plot the geometrical object
% the plotting style depends on the data content
figure;
switch dtype
  case 'volume'
    funparam = [];
    if isfield(data, 'anatomy')
      funparam = data.anatomy;
    elseif isfield(data, 'gray')
      funparam = data.gray;
    elseif isfield(data, 'white')
      funparam = data.white;
    elseif isfield(data, 'brick0')
      funparam = data.brick0; % used for an atlas
    elseif isfield(data, 'brick1')
      funparam = data.brick1; % used for an atlas
    else
      % try to determine it automatically
      fn = fieldnames(data);
      for i=1:length(fn)
        if isequal(size(data.(fn{i})), data.dim)
          funparam = data.(fn{i});
          break;
        end
      end
    end

    if isempty(funparam)
      ft_error('don''t know which volumetric parameter to plot');
    end

    % the volumetric data needs to be interpolated onto three orthogonal planes
    % determine a resolution that is close to, or identical to the original resolution
    [corner_vox, corner_head] = cornerpoints(data.dim, data.transform);
    diagonal_head = norm(range(corner_head));
    diagonal_vox  = norm(range(corner_vox));
    resolution    = (diagonal_head+eps)/(diagonal_vox+eps); % this is in units of "data.unit"

    % scale funparam between 0 and 1
    if ~isa(funparam, 'double') % avoid integer datatypes to allow for scaling
      funparam = double(funparam);
    end
    dmin = min(funparam(:));
    dmax = max(funparam(:));
    funparam  = (funparam-dmin)./(dmax-dmin);

    clear ft_plot_slice
    ft_plot_ortho(funparam, 'transform', data.transform, 'unit', data.unit, 'resolution', resolution, 'style', 'intersect', 'clim', clim);
    axis vis3d
    view([110 36]);

  case 'source'
    if isfield(data, 'inside') && ~isfield(data, 'tri')
      % only plot the source locations that are inside the volume conduction model
      ft_plot_mesh(data.pos(data.inside, :));
    else
      ft_plot_mesh(data, 'edgecolor','none', 'facecolor', [0.6 0.8 0.6], 'facealpha', 0.6);
    end
    camlight;

  case 'dip'
    ft_plot_mesh(data, 'edgecolor','none', 'facecolor', 'none');
    camlight;

  case 'headshape'
    ft_plot_headshape(data);
    camlight;

  case {'mesh', 'source+mesh'}
    ft_plot_mesh(data);
    camlight;

  case 'headmodel'
    ft_plot_headmodel(data);
    camlight;

  case {'grad' 'elec' 'opto' 'sens'}
    ft_plot_sens(data, 'label', 'label');
    camlight;

  case {'raw', 'timelock', 'freq', 'mvar', 'freqmvar', 'comp'}
    % the data may contain a gradiometer or electrode definition
    if isfield(data, 'grad')
      ft_plot_sens(data.grad);
    elseif isfield(data, 'elec')
      ft_plot_sens(data.elec, 'label', 'label');
    elseif isfield(data, 'opto')
      ft_plot_sens(data.opto, 'label', 'label');
    end

  case 'unknown'
end % switch dtype{k}

if isfield(data, 'tri')
  % this makes the 3-D object easier to understand
  camlight
  lighting gouraud
end

% plot the 3-D axes, labels, and sphere at the origin
ft_plot_axes(data, 'axisscale', axisscale, 'fontsize', fontsize);

if istrue(dointeractive)
  % ensure the figure is updated prior to asking the question
  % this was needed for FT_ELECTRODEPLACEMENT in combination with MATLAB 2022a
  drawnow

  if ~isfield(data, 'coordsys') || isempty(data.coordsys)
    % default is yes
    fprintf('The coordinate system is not specified.\n')
    value = smartinput('Do you want to change the anatomical labels for the axes [Y, n]? ', 'y');
  else
    % default is no
    fprintf('The coordinate system is specified as "%s".\n', data.coordsys)
    value = smartinput('Do you want to change the anatomical labels for the axes [y, N]? ', 'n');
  end

  if strcmpi(value, 'n')
    return
  end

  % interactively determine orientation
  orientation = '   ';
  while ~any(strcmp(orientation(1), {'r', 'l', 'a', 'p', 's', 'i'}))
    orientation(1) = smartinput('What is the anatomical label for the positive X-axis [r, l, a, p, s, i]? ', '');
  end
  while ~any(strcmp(orientation(2), {'r', 'l', 'a', 'p', 's', 'i'}))
    orientation(2) = smartinput('What is the anatomical label for the positive Y-axis [r, l, a, p, s, i]? ', '');
  end
  while ~any(strcmp(orientation(3), {'r', 'l', 'a', 'p', 's', 'i'}))
    orientation(3) = smartinput('What is the anatomical label for the positive Z-axis [r, l, a, p, s, i]? ', '');
  end

  % interactively determine origin
  origin = ' ';
  while ~any(strcmp(origin, {'a', 'i', 's', 'n'}))
    origin = input('Is the origin of the coordinate system at the a(nterior commissure), i(nterauricular), s(scanner origin), n(ot a landmark)? ', 's');
  end

  % some coordinate systems are identical or very similar, see https://www.fieldtriptoolbox.org/faq/coordsys
  if origin=='a' && strcmp(orientation, 'ras')
    coordsys = 'acpc'; % also used for spm, mni, tal
  elseif origin=='i' && strcmp(orientation, 'als')
    coordsys = 'ctf'; % also used for 4d, bti, eeglab
  elseif origin=='i' && strcmp(orientation, 'ras')
    coordsys = 'neuromag'; % also used for itab
  elseif origin=='s' && strcmp(orientation, 'ras')
    coordsys = 'scanras'; % also used for nifti
  elseif origin=='s' && strcmp(orientation, 'lps')
    coordsys = 'scanlps'; % also used for dicom
  else
    % only use the orientation, not the origin
    coordsys = orientation;
  end

  data.coordsys = coordsys;
end % if interactive