text:
%SerialLink.fkine Forward kinematics
%
% T = R.fkine(Q, OPTIONS) is the pose of the robot end-effector as an SE(3)
% homogeneous transformation (4x4) for the joint configuration Q (1xN).
%
% If Q is a matrix (KxN) the rows are interpreted as the generalized joint
% coordinates for a sequence of points along a trajectory.  Q(i,j) is the
% j'th joint parameter for the i'th trajectory point.  In this case T is a
% 3d matrix (4x4xK) where the last subscript is the index along the path.
%
% [T,ALL] = R.fkine(Q) as above but ALL (4x4xN) is the pose of the link
% frames 1 to N, such that ALL(:,:,k) is the pose of link frame k.
%
% Options::
%  'deg'    Assume that revolute joint coordinates are in degrees not
%           radians
%
% Note::
% - The robot's base or tool transform, if present, are incorporated into the
%   result.
% - Joint offsets, if defined, are added to Q before the forward kinematics are
%   computed.
%
% See also SerialLink.ikine, SerialLink.ikine6s.

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

function [t allt] = fkine(robot, q, varargin)

%
% evaluate fkine for each point on a trajectory of
% theta_i or q_i data
%

n = robot.n;

opt.deg = false;

opt = tb_optparse(opt, varargin);

if opt.deg
    % in degrees mode, scale the columns corresponding to revolute axes
    k = ~robot.links.isprismatic;
    q(:,k) = q(:,k) * pi/180;
end

if nargout > 1
    allt = zeros(4,4,n);
    if isa(q,'sym')
        allt = sym(allt);
    end
end

L = robot.links;
if numel(q) == n
    t = robot.base;
    for i=1:n
        t = t * L(i).A(q(i));

        if nargout > 1
            allt(:,:,i) = t; % intermediate transformations
        end
    end
    t = t * robot.tool;
else
    if numcols(q) ~= n
        error('q must have %d columns', n)
    end
    t = zeros(4,4,0);
    for qv=q',		% for each trajectory point
        tt = robot.base;
        for i=1:n
            tt = tt * L(i).A(qv(i));
        end
        t = cat(3, t, tt * robot.tool);
    end
end

if isa(t, 'sym')
    t = simplify(t);
end

%robot.T = t;
%robot.notify('Moved');

author: Allopart

repo: rbpf-gmapping


text:
function r = lt(p,q)
% this function compares the less than condition of x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

if (isa(p, 'FaceVariable')&&isa(q, 'FaceVariable'))
    r=p;
    r.xvalue = p.xvalue<q.xvalue;
    r.yvalue = p.yvalue<q.yvalue;
    r.zvalue = p.zvalue<q.zvalue;
elseif isa(p, 'FaceVariable')
    r=p;
    r.xvalue = p.xvalue<q;
    r.yvalue = p.yvalue<q;
    r.zvalue = p.zvalue<q;
else
    r=q;
    r.xvalue = p<q.xvalue;
    r.yvalue = p<q.yvalue;
    r.zvalue = p<q.zvalue;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
classdef ImageRatio < Algorithms.CDAlg
    methods
        function obj = ImageRatio()
            obj@Algorithms.CDAlg('Image Ratioing', '')
        end
        DI = detectChange(obj, ~, ~);
    end
end

author: Bobholamovic

repo: ChangeDetectionToolbox


text:
function vars = depends(F)

vars = depends(lmi(F));

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function D = mahalDist(x, m, C, use_log)
% p=gaussian_prob(x, m, C, use_log)
%
% Evaluate the multi-variate density with mean vector m and covariance
% matrix C for the input vector x.
% Vectorized version: Here X is a matrix of column vectors, and p is
% a vector of probabilities for each vector.

if nargin<4, use_log = 0; end

d   = length(m);

if size(x,1)~=d
   x=x';
end
N       = size(x,2);

m       = m(:);
M       = m*ones(1,N);
denom   = (2*pi)^(d/2)*sqrt(abs(det(C)));
invC    = inv(C);
mahal   = sum(((x-M)'*invC).*(x-M)',2);   % Chris Bregler's trick

switch use_log,
case 2,
  D     = mahal;
case 1,
  D     = -0.5*mahal - log(denom);
case 0,
  numer = exp(-0.5*mahal);
  D     = numer/denom;
otherwise
    error('Unsupported log type')
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function yp = p25_fun ( neqn, t, y )

%*****************************************************************************80
%
%% P25_FUN evaluates the function for problem P25.
%
%  Discussion:
%
%    2 equations
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Wayne Enright, John Pryce,
%    Algorithm 648,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 1, pages 28-34.
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T, Y(NEQN), the arguments of the derivative
%    function.
%
%    Output, real YP(NEQN), the value of the derivative function.
%
  yp = zeros ( neqn, 1 );

  yp(1) = y(2);
  yp(2) = sqrt ( 1.0 + y(2).^2 ) / ( 25.0 - t );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [cost,costs] = crossvalidate(model, L, estfct,combinefct)

% Estimate the model performance of a model with [$ l$] -fold crossvalidation
%
% CAUTION!! Use this function only to obtain the value of the crossvalidation score
% function given the tuning parameters. Do not use this function together with
% 'tunelssvm', but use 'crossvalidatelssvm' instead. The latter is a faster
% implementation which uses previously computed results.
%
% >> cost = crossvalidate({Xtrain,Ytrain,type,gam,sig2})
% >> cost = crossvalidate( model)
%
% The data is once permutated randomly, then it is divided into L (by default 10)
% disjoint sets. In the i-th (i=1,...,l) iteration, the i-th set is used to estimate
% the performance ('validation set') of the model trained on the other l-1 sets ('training set').
% Finally, the l (denoted by L) different estimates of the performance are combined (by default by the 'mean').
% The assumption is made that the input data are distributed independent and identically over the
% input space. As additional output, the costs in the different folds ('costs') of the data are returned:
%
% >> [cost, costs] = crossvalidate(model)
%
% Some commonly used criteria are:
%
% >> cost = crossvalidate(model, 10, 'misclass', 'mean')
% >> cost = crossvalidate(model, 10, 'mse', 'mean')
% >> cost = crossvalidate(model, 10, 'mae', 'median')
%
% Full syntax
%
%     1. Using LS-SVMlab with the functional interface:
%
% >> [cost, costs] = crossvalidate({X,Y,type,gam,sig2,kernel,preprocess}, L, estfct, combinefct)
%
%       Outputs
%         cost          : Cost estimation of the L-fold cross validation
%         costs(*)      : L x 1 vector with costs estimated on the L different folds
%
%       Inputs
%         X             : Training input data used for defining the LS-SVM and the preprocessing
%         Y             : Training output data used for defining the LS-SVM and the preprocessing
%         type          : 'function estimation' ('f') or 'classifier' ('c')
%         gam           : Regularization parameter
%         sig2          : Kernel parameter (squared bandwidth in the case of the 'RBF_kernel')
%         kernel(*)     : Kernel type (by default 'RBF_kernel')
%         preprocess(*) : 'preprocess'(*) or 'original'
%         L(*)          : Number of folds (by default 10)
%         estfct(*)     : Function estimating the cost based on the residuals (by default mse)
%         combinefct(*) : Function combining the estimated costs on the different folds (by default mean)
%
%
%     2. Using the object oriented interface:
%
% >> [cost, costs] = crossvalidate(model, L, estfct, combinefct)
%
%       Outputs
%         cost          : Cost estimation of the L-fold cross validation
%         costs(*)      : L x 1 vector with costs estimated on the L different folds
%
%       Inputs
%         model         : Object oriented representation of the LS-SVM model
%         L(*)          : Number of folds (by default 10)
%         estfct(*)     : Function estimating the cost based on the residuals (by default mse)
%         combinefct(*) : Function combining the estimated costs on the different folds (by default mean)
%
%
%     3. Using other modeling techniques::
%
%
% See also:
% leaveoneout, gcrossvalidate, trainlssvm, simlssvm

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

% LS-SVMlab
eval('model = initlssvm(model{:});',' ');
eval('L;','L=min(ceil(sqrt(model.nb_data)),10);');
eval('estfct;','estfct=''mse'';');
eval('combinefct;','combinefct=''mean'';');

%
% initialisation and defaults
%
nb_data = size(model.ytrain,1);
d = size(model.xtrain,2);

if L==nb_data, p = 1:nb_data; else p = randperm(nb_data); end
px = model.xtrain(p,:);
py = model.ytrain(p,:);

[~,Y] = postlssvm(model,[],py);

%initialize: no incremental  memory allocation
costs = zeros(L,length(model.gam));
block_size = floor(nb_data/L);

% calculate matrix for LS-SVM once for the entire data
S = ones(nb_data,1);
Inb = eye(nb_data);
K = kernel_matrix(px,model.kernel_type,model.kernel_pars);
Atot = K+Inb./model.gam;

% Cholesky factor
try R = chol(Atot);
    % Solve full system
    q = R\(R'\[py S]);
    p = q(:,2); q = q(:,1);
    s  = 1/sum(p);
    bias  = s*sum(q);
    alpha  = q - p*bias;

    % Two expensive steps yet more efficient that using LINSOLVE on each fold
    Ri = R\Inb;
    C = Ri*Ri' - s*(p)*p';

catch %R = cholinc(sparse(Atot),1e-5);
    A = [K+Inb./model.gam S; S' 0];
    C = pinv(A);
    alpha = C*[py;0];
    %bias = alpha(nb_data+1);
    alpha = alpha(1:nb_data);
end

% Solve full system
q = R\(R'\[py S]);
p = q(:,2); q = q(:,1);
s  = 1/sum(p);
bias  = s*sum(q);
alpha  = q - p*bias;

% Two expensive steps yet more efficient that using LINSOLVE on each fold
Ri = R\Inb;
C = Ri*Ri' - s*(p)*p';

% start loop over l validations
for l = 1:L,
    % divide data in validation set and trainings data set
    if l==L,
        validation = block_size*(l-1)+1:nb_data;
    else
        validation = block_size*(l-1)+1:block_size*l;
    end
    % Submatrix of C to compute residuals for the l-th fold left out
    Ckk = C(validation,validation);
    % Solution of small linear system (block_size x block_size)
    try % faster
        Rkk = chol(Ckk+eps);
        betak = Rkk\(Rkk'\alpha(validation));
    catch
        betak = Ckk\alpha(validation);
    end
    % latent outputs for validation
    yh = py(validation) - betak;
    [~,yh] = postlssvm(model,[],yh);
    if ~(model.type(1)=='c')
        costs(l,1) = feval(estfct,yh - Y(validation,:));
    else
        costs(l,1) = feval(estfct,Y(validation,:),sign(yh));
    end
end
cost = feval(combinefct,costs);

author: peterhcharlton

repo: RRest


text:
function [params] = gui_CalcOmori(params,nNodeStart, nNodeEnd)
    % gui_CalcOmori Calculate Omori parameters, here setting up the data
    %
    % [params] = gui_CalcOmori(params,nNodeStart, nNodeEnd);
    % ----------------------------------------------------------------
    %
    % Input parameters:
    %   params.mCatalog           Earthquake catalog
    %   params.mPolygon           Polygon (defined by ex_selectgrid)
    %   params.vX                 X-vector (defined by ex_selectgrid)
    %   params.vY                 Y-vector (defined by ex_selectgrid)
    %   params.vUsedNodes         Used nodes vX * vY defining the mPolygon (defined by ex_selectgrid)
    %   params.nCalculation       Number of random simulations
    %   params.bMap               Calculate a map (true) or a cross-section (false)
    %   params.bNumber            Use constant number (true) or constant radius (false)
    %   params.nNumberEvents      Number of earthquakes if bNumber is true
    %   params.fMaxRadius         Maximum Radius using a constant number of events; works only with bNumber is true
    %   params.fRadius            Radius of gridnode if bNumber is false
    %   params.nMinimumNumber     Minimum number of earthquakes per node for determining a b-value
    %   params.fMinMag            Lower limit of magnitude range for testing
    %   params.fMaxMag            Upper limit of magnitude range for testing
    %   params.bTimePeriod        Calculate seismicity difference for 2 periods (0) until start and end of catalog or
    %                             a specific time period before and after fSplitTime (1)
    %   params.fTimePeriodDays        Length of time periods
    %   params.bTstart            Check for starting time of temporal mapping
    %   params.fTstart            Starting time for temporal mapping
    %   params.bBstnum            Check for boostrap sampling
    %   params.fBstnum            Number of bootstrap samples
    %   params.fBinning           Bin size for magnitude binning
    %   params.sComment           Comment on calculation

    % Output parameters:
    %   Same as input parameters including
    %   params.mValueGrid         Matrix of calculated values
    %   params.vcsGridNames       Names of parameters calculated
    %   Check gui_NodeCalcOmori.m for a list of variables!!
    %
    % J. Woessner; j.woessner@sed.ethz.ch
    % updated: 16.02.2006

    report_this_filefun();

    % Check calculation for splitting of nodes
    if nargin < 2
        nNodeStart = 1;
        nNodeEnd = length(params.mPolygon(:,1));
    end

    % Initialize
    vResults = [];
    params.sComment = [];
    if isempty(params.fBinning)
        params.fBinning = 0.1;
    end

    % Determine time period of catalog
    params.fTminCat = min(params.mCatalog.Date);
    params.fTmaxCat = max(params.mCatalog.Date);

    % Init result matrix
    mValueGrid_ = [];

    % Temporary saving the original catalog
    mCatalog = params.mCatalog;

    % Force saving all 100 nodes
    fDivide = length(params.mPolygon(:,1))/100
    fForceSave = length(params.mPolygon(:,1))/fDivide;

    % % Loop over time
    fTstart = params.fTstart;
    % while fTstart < params.fTmaxCat
    % mValueGrid_ = [];
    % params.mCatalog = mCatalog;

    % Create Indices to catalog
    [params.caNodeIndices] = ex_CreateIndexCatalog(params.mCatalog, params.mPolygon, params.bMap, params.nGriddingMode, ...
        params.nNumberEvents, params.fRadius, params.fSizeRectHorizontal, params.fSizeRectDepth);
    % Loop over all grid nodes
    hWaitbar1 = waitbar(0,'Calculating nodes...');
    set(hWaitbar1,'Numbertitle','off','Name','Node percentage');
    for nNode_ = nNodeStart:nNodeEnd
        % Create node catalog
        mNodeCatalog_ = params.mCatalog(params.caNodeIndices{nNode_}, :);
        % Check for constant number of events calculations
        if (params.nGriddingMode == 0)
            [mNodeCatalog_] = ex_CheckMaxRadius(mNodeCatalog_, params.mPolygon, nNode_, params.caNodeIndices, params.fMaxRadius, params.nNumberEvents, params.bMap);
        end
        % Check for number of events in aftershock sequence
        % Create catalog after split time (aftershock sequence)
        vSelAf = (params.fTstart <= mNodeCatalog_(:,3) & mNodeCatalog_(:,3) < params.fTstart+params.fTimePeriodDays);
        [nX,nY] = size(mNodeCatalog_(vSelAf,:));
        if (nX < params.nMinimumNumber)
            [nNumBgr,nY2] = size(mNodeCatalog_(~vSelAf,:));
            mValueGrid_= [mValueGrid_; NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN nNumBgr nNumBgr/params.fTimeBgr NaN NaN];
        else
            [rCalcNodeResult_] = gui_NodeCalcOmori(params,mNodeCatalog_);
            % Store the results
            mValueGrid_= [mValueGrid_;rCalcNodeResult_.pval1 rCalcNodeResult_.cval1 rCalcNodeResult_.kval1 rCalcNodeResult_.nNumEvents rCalcNodeResult_.fLogEqdens...
                rCalcNodeResult_.H rCalcNodeResult_.pmean1 rCalcNodeResult_.cmean1 rCalcNodeResult_.kmean1...
                rCalcNodeResult_.pmeanStd1 rCalcNodeResult_.cmeanStd1 rCalcNodeResult_.kmeanStd1...
                rCalcNodeResult_.fTafseq rCalcNodeResult_.fTafseqmean rCalcNodeResult_.nNumAf rCalcNodeResult_.nNumBgr...
                rCalcNodeResult_.fBgrate rCalcNodeResult_.fLog10Bgrate rCalcNodeResult_.fLog10Tafseq];
        end % End of if on nNode_
        if rem(nNode_,floor(fForceSave)) == 0
            waitbar(nNode_/length(params.mPolygon(:,1)))
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Temporary saving
            params.vcsGridNames = cellstr(char('p-value','c-value','k-value','Number of events','log(EQ_density)',...
                'H','p-mean','c-mean','k-mean','\sigma (p)','\sigma (c)','\sigma (k)','t_a [years]','t_a(Bst) [years]',...
                'Number of aftershock','Num. events background','Background rate [year]','Log10(Bgr rate) [year]',...
                'log10(t_a) [years]'));
            params.mValueGrid = mValueGrid_;
            % Add parameter to params.sComment
            if  params.nGriddingMode == 0;   % Constant number
                params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Constant number: ' num2str(params.nNumberEvents) ', MaxRadius: '...
                    num2str(params.fMaxRadius) ' km'];
                vResults = params;
                save(['tmp_result_Time' num2str(fTstart) '_Constnum_' num2str(params.nNumberEvents) '_MaxRad_' num2str(params.fMaxRadius)...
                    '_Nmin_' num2str(params.nMinimumNumber) '_Node' num2str(nNode_) '.mat'], 'vResults');
            elseif params.nGriddingMode == 1;   % Constant radius
                params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Radius: ' num2str(params.fRadius) ' km, Nmin: ' num2str(params.nMinimumNumber)];
                vResults = params;
                save(['tmp_result_Time' num2str(fTstart) '_Rad_' num2str(params.fRadius) '_Nmin_' num2str(params.nMinimumNumber)  '_Node' num2str(nNode_) '.mat'], 'vResults');
            else  % Rectangle mode
                params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Rect. X: ' num2str(params.fSizeRectHorizontal) ' km, Rect. Y: ' num2str(params.fSizeRectDepth)...
                    ' km, Nmin: ' num2str(params.nMinimumNumber)];
                vResults = params;
                save(['tmp_result_Time' num2str(fTstart) '_RectX_' num2str(params.fSizeRectHorizontal) '_RectY_' num2str(params.fSizeRectDepth)...
                    '_Nmin_' num2str(params.nMinimumNumber)  '_Node' num2str(nNode_) '.mat'], 'vResults');
            end % END of params.nGriddingmode
            vResults =[];
        end % End updating waitbar
    end % for nNode
    close(hWaitbar1);
    % Parameter description
    params.vcsGridNames = cellstr(char('p-value','c-value','k-value','Number of events','log(EQ_density)',...
        'H','p-mean','c-mean','k-mean','\sigma (p)','\sigma (c)','\sigma (k)','t_a [years]','t_a(Bst) [years]',...
        'Number of aftershock','Num. events background','Background rate [year]','Log10(Bgr rate) [year]',...
        'log10(t_a) [years]'));
    params.mValueGrid = mValueGrid_;
    % Add parameter to params.sComment
    if  params.nGriddingMode == 0;   % Constant number
        params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
            ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Constant number: ' num2str(params.nNumberEvents) ', MaxRadius: '...
            num2str(params.fMaxRadius) ' km'];
        vResults = params;
        save(['result_Time' num2str(fTstart) '_Constnum_' num2str(params.nNumberEvents) '_MaxRad_' num2str(params.fMaxRadius)...
            '_Nmin_' num2str(params.nMinimumNumber) '_Nodes_' num2str(nNodeStart) '_' num2str(nNodeEnd) '.mat'], 'vResults');
    elseif params.nGriddingMode == 1;   % Constant radius
        params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
            ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Radius: ' num2str(params.fRadius) ' km, Nmin: ' num2str(params.nMinimumNumber)];
        vResults = params;
        save(['result_Time' num2str(fTstart) '_Rad_' num2str(params.fRadius) '_Nmin_' num2str(params.nMinimumNumber)...
            '_Nodes_' num2str(nNodeStart) '_' num2str(nNodeEnd) '.mat'], 'vResults');
    else  % Rectangle mode
        params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
            ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Rect. X: ' num2str(params.fSizeRectHorizontal) ' km, Rect. Y: ' num2str(params.fSizeRectDepth)...
            ' km, Nmin: ' num2str(params.nMinimumNumber)];
        vResults = params;
        save(['result_Time' num2str(fTstart) '_RectX_' num2str(params.fSizeRectHorizontal) '_RectY_' num2str(params.fSizeRectDepth)...
            '_Nmin_' num2str(params.nMinimumNumber) '_Nodes_' num2str(nNodeStart) '_' num2str(nNodeEnd) '.mat'], 'vResults');
    end

author: CelsoReyes


text:
%ASSEMBLEPROB Assembly of system matrix and right hand side/load vector.
%
%   [ M, A, F, T_M, T_A, T_F, T_SP ] = ASSEMBLEPROB( PROB, VARARGIN ) Calls the
%   assembly routines to compute and assemble a monolithic sparse system
%   matrix and a right hand side/load vector as specified in the PROB struct.
%   Accepts optional propery value pairs in VARARGIN.
%
%       Input       Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       prob        struct                 Problem definition struct
%       icub        scalar/{auto}          Numnerical integration rule used in assembly
%                                                Default 1+max(shape function order)
%       imass       scalar {1}             Mass matrix lumping:  1 = Full mass matrix
%                                          2 = row sum lumping,  3 = diagonal lumping
%                                          4 = HRZ diagonal lumping
%       n_cmax      scalar {50000}         Max number of cells to assemble
%                                          for at once (to limit memory consumption)
%       f_m         logical {false}        Assembly flag for mass matrix
%       f_a         logical {false}        Assembly flag for system matrix
%       f_f         logical {false}        Assembly flag for load vector
%       f_c         logical {true}         Assembly flag for integral constraints
%       f_sparse    logical {false}        Return sparse/struct matrix format
%       solcomp     {all dvars/subd}       Dependent variables/subdomains to assemble for
%                                                                                .
%       Output      Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       M           sparse [n_M]           Assembled mass matrix
%       A           sparse [n_A]           Assembled system matrix
%       f           [n_A,1]                Assembled rhs/load vector
%       t_m         scalar                 Time spent assembling mass matrix
%       t_a         scalar                 Time spent assembling system matrix
%       t_f         scalar                 Time spent assembling rhs/load vector
%       t_sp        scalar                 Time for sparse matrix conversion
%
%   See also ASSEMBLEA, ASSEMBLEF, ASSEMMAT

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function showposlat(model, start, pos, fp_count, overlap)

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2009-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

conf = voc_config();

% get training data
if nargin < 3
  pos = pascal_data(model.class, model.year);
end

numpos = length(pos);
model.interval = conf.training.interval_fg;
pixels = model.minsize * model.sbin / 2;
minsize = prod(pixels);
if nargin < 2
  start = 1;
end

if nargin < 5
  overlap = 0.7;
end
show_slabs = true;

% compute latent filter locations and record target bounding boxes
k = 0;
for i = start:numpos
  fprintf('%s %s: show pos lat: %d/%d\n', procid(), model.class, i, numpos);
  if max(pos(i).sizes) < minsize
    fprintf('  all too small\n');
    continue;
  end

  % do image level operations
  im = color(imreadx(pos(i)));
  [im, boxes] = croppos(im, pos(i).boxes);
  [pyra, model_dp] = gdetect_pos_prepare(im, model, boxes, overlap);
  %fprintf('%s\n', pos(i).im);

  % process each box in the image
  for b = 1:size(boxes,1)
    fg_ibox = b;
    bg_iboxes = 1:size(boxes,1);
    bg_iboxes(b) = [];

    fg_box = boxes(b,:);
    bg_boxes = boxes;
    bg_boxes(b,:) = [];

    if show_slabs
      boxesc = [fg_box zeros(1,model.numfilters*4-4+2) 2];
    else
      boxesc = [fg_box 2];
    end
    if ~isempty(bg_boxes)
      if show_slabs
        boxesc = cat(1, boxesc, padarray(bg_boxes, [0 model.numfilters*4-4+2+1], 3, 'post'));
      else
        boxesc = cat(1, boxesc, padarray(bg_boxes, [0 1], 3, 'post'));
      end
    end

    % skip small examples
    if pos(i).sizes(b) < minsize
      fprintf('  %d: too small\n', b);
      continue;
    end
    [det, bs, trees] = gdetect_pos(pyra, model_dp, 1+fp_count, fg_ibox, ...
                                  overlap, bg_iboxes, 0.5);
    if ~isempty(det)
      o = boxoverlap(clipboxes(im, det), fg_box);
      fprintf('  %d: comp %d  score %.4f  overlap=%0.4f\n', ...
              b, bs(1,end-1), bs(1,end), o(1));
      if size(det,1) > 1
        %losses = lossfunc(o);
        losses = [];
        for ii = 1:length(trees)
          losses = [losses; trees{ii}(13,1)];
        end
        for j = 2:length(o)
          score = det(j,end);
          fprintf('  > %d: score=%0.4f (%0.4f + %0.4f)  comp %d  overlap=%0.4f\n', ...
                  j, score, score-losses(j), losses(j), det(j,end-1), o(j));
        end
        if show_slabs
          boxesc = cat(1, boxesc, padarray(bs(2:end,:), [0 1], 1, 'post'));
        else
          boxesc = cat(1, boxesc, padarray(det(2:end,1:4), [0 1], 1, 'post'));
        end
      end
      if show_slabs
        boxesc = cat(1, boxesc, padarray(bs(1,:), [0 1], 0, 'post'));
        boxesc = cat(1, boxesc, padarray(det(1,1:4), [0 model.numfilters*4-4+2+1], 4, 'post'));
      else
        boxesc = cat(1, boxesc, padarray(det(1,1:4), [0 1], 0, 'post'));
      end
      %showboxes(im, bs(1,:));
      tree = trees{1};
    else
      fprintf('  %d: no overlap\n', b);
      tree = [];
    end

    subplot(1,2,1);
    showboxesc(im, boxesc);
    title('green = fg box;  blue = bg boxes;  red = loss adjusted boxes;  cyan = label completed detection');
    subplot(1,2,2);
    if ~isempty(tree)
      vis_derived_filter(model, tree);
    end

    pause;
  end
end

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function ry = vecroty(angle)
%
% Function Name:
%
%   vecroty - Transformation matrix for a rotation around the y axis.
%
% Calling Sequence:
%
%   ry = vecroty(angle);
%
% Parameters:
%
%   angle		: rotation angle defined in radians
%
%   ry		: (4x4) Transformation matrix.
%
%
% Description:
%
%   Return the (4x4) Transformation matrix for a rotation about the y axis
%   by the defined angle.
%
%   The matrix is:
%
%         [  cos(angle)       0        sin(angle)       0]
%         [      0            1            0            0]
%         [ -sin(angle)       0        cos(angle)       0]
%         [      0            0            0            1]
%
% Examples:
%
%    Rotate the NURBS line (0.0 0.0 0.0) - (3.0 3.0 3.0) by 45 degrees
%    around the y-axis
%
%    line = nrbline([0.0 0.0 0.0],[3.0 3.0 3.0]);
%    trans = vecroty(%pi/4);
%    rline = nrbtform(line, trans);
%
% See:
%
%    nrbtform

%  Dr D.M. Spink
%  Copyright (c) 2000.

sn = sin(angle);
cn = cos(angle);
ry = [cn 0 sn 0; 0 1 0 0; -sn 0 cn 0; 0 0 0 1];

author: Sable

repo: mcbench-benchmarks


text:
% This file serves as a template of defining the topology of the
% beamforming network with cross entropy training.
% You should create a copy for each of your experiments and name them
% differently.
%
% Created by Xiong Xiao, Temasek Laboratories, NTU, Singapore.
% Last Modified: 29 Nov 2016
%
function para = ConfigDereverbNet_Regression(para)
para.topology = SetDefaultValue(para.topology, 'fs', 16000);
if isfield(para.topology, 'useChannel')
    para.topology.nCh = length(para.topology.useChannel);
else
    para.topology.nCh = 1;      % by default use 2 channels
end
para.topology = SetDefaultValue(para.topology, 'useFileName', 0);   % by default load waveforms into memory. If data is too big, we can also use wave file names

para.topology = SetDefaultValue(para.topology, 'fft_len', 512);
para.topology.freqBin = (0:1/para.topology.fft_len:0.5)*2*pi;
para.topology.nFreqBin = length(para.topology.freqBin);
% define the parameters for extracting Fourier coefficients
para.topology.frame_len = para.topology.fs * 0.025;
para.topology.frame_shift = para.topology.fs * 0.01;
para.topology.removeDC = 0;     % do not remove DC for faster speed.
para.topology.win_type = 'hamming';

% define the regression network type
para.topology = SetDefaultValue(para.topology, 'RegressionNetType', 'LSTM');
switch para.topology.RegressionNetType
    case 'DNN'
        para.topology = SetDefaultValue(para.topology, 'hiddenLayerSize', [1024]);
        para.topology = SetDefaultValue(para.topology, 'contextSize', 11);  % for DNN, we use 11 frames of consecutive frames
    case 'LSTM'
        para.topology = SetDefaultValue(para.topology, 'hiddenLayerSize', [512]);
        para.topology = SetDefaultValue(para.topology, 'useDelta', 1);    % for LSTM, we use delta features without splicing
end

end

author: singaxiong

repo: SignalGraph


text:
%% Load Model
data_path = '../';

%% Load Sample
sample_name = 'image00002';
img = imread([data_path sample_name '.jpg']);
load([data_path sample_name '.mat']);
[height, width, nChannels] = size(img);

imshow(img);
hold on;
plot(pt3d_68(1,:), pt3d_68(2,:), 'b.');
hold off;

author: XgTu


text:
function out=diag(in,K)
if isempty(in),out=in;return;end
precision=in(1).precision;
zero=mp(0,precision);
if nargin==1,K=0;end
s=size(in);
switch min(s)
 case 1 %a vector, either row or column
  N=length(in);
  out=zeros(mp(N),mp(N));
  SS.type='()';SS.subs={1:(N+1):(N*N)};
  out=subsasgn(out,SS,in);
 otherwise %a matrix
  out=zero*[];
  for ii=1:size(in,1)
   if (ii+K)<=size(in,2)
    %Force column vector
    out(ii,1)=in(ii,(ii+K));
   end
  end
end

author: opencobra

repo: cobratoolbox


text:
function [estimate] = ft_inverse_music(sourcemodel, sens, headmodel, dat, varargin)

% FT_INVERSE_MUSIC source localization using MUltiple SIgnal Classification.
% This is a signal subspace method, which covers the techniques for
% multiple source localization by using the eigen-structure of the
% measured data matrix.
%
% Use as
%   [estimate] = ft_inverse_music(sourcemodel, sens, headmodel, dat, ...)
% where
%   sourcemodel is the input source model, see FT_PREPARE_SOURCEMODEL
%   sens        is the gradiometer or electrode definition, see FT_DATATYPE_SENS
%   headmodel   is the volume conductor definition, see FT_PREPARE_HEADMODEL
%   dat         is the data matrix with the ERP or ERF
% and
%   estimate    contains the estimated source parameters
%
% Additional input arguments should be specified as key-value pairs and can include
%   'cov'              = data covariance matrix
%   'numcomponent'     = integer number
%   'feedback'         = can be 'none', 'gui', 'dial', 'textbar', 'text', 'textcr', 'textnl' (default = 'text')
%
% These options influence the forward computation of the leadfield
%   'reducerank'      = 'no' or number  (default = 3 for EEG, 2 for MEG)
%   'backproject'     = 'yes' or 'no', in the case of a rank reduction this parameter determines whether the result will be backprojected onto the original subspace (default = 'yes')
%   'normalize'       = 'no', 'yes' or 'column' (default = 'no')
%   'normalizeparam'  = parameter for depth normalization (default = 0.5)
%   'weight'          = number or Nx1 vector, weight for each dipole position to compensate for the size of the corresponding patch (default = 1)
%
% This implements
% - J.C. Mosher, P.S. Lewis and R.M. Leahy, "Multiple dipole modeling and
%   localization from spatiotemporal MEG data", IEEE Trans. Biomed. Eng.,
%   pp 541-557, June, 1992.
%
% See also FT_SOURCEANALYSIS, FT_PREPARE_HEADMODEL, FT_PREPARE_SOURCEMODEL

% Copyright (C) 2004-2008, 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 mod(nargin-4,2)
  % the first 4 arguments are fixed, the other arguments should come in pairs
  ft_error('invalid number of optional arguments');
end

% get the optional input arguments, or use defaults
cov            = ft_getopt(varargin, 'cov');
numcomponent   = ft_getopt(varargin, 'numcomponent');     % this is required, see below
feedback       = ft_getopt(varargin, 'feedback', 'text');

% construct the low-level options for the leadfield computation as key-value pairs, these are passed to FT_COMPUTE_LEADFIELD
leadfieldopt = {};
leadfieldopt = ft_setopt(leadfieldopt, 'reducerank',     ft_getopt(varargin, 'reducerank'));
leadfieldopt = ft_setopt(leadfieldopt, 'backproject',    ft_getopt(varargin, 'backproject'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalize',      ft_getopt(varargin, 'normalize'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalizeparam', ft_getopt(varargin, 'normalizeparam'));
leadfieldopt = ft_setopt(leadfieldopt, 'weight',         ft_getopt(varargin, 'weight'));

if isempty(numcomponent)
  ft_error('you must specify the number of signal components');
end

% flags to avoid calling isfield repeatedly in the loop over grid positions (saves a lot of time)
hasmom        = isfield(sourcemodel, 'mom');
hasleadfield  = isfield(sourcemodel, 'leadfield');
hasfilter     = false; % not used here

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find the dipole positions that are inside/outside the brain
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isfield(sourcemodel, 'inside')
  if hasfilter
    sourcemodel.inside = ~cellfun(@isempty, sourcemodel.filter);
  elseif hasleadfield
    sourcemodel.inside = ~cellfun(@isempty, sourcemodel.leadfield);
  else
    sourcemodel.inside = ft_inside_headmodel(sourcemodel.pos, headmodel);
  end
end

% convert to logical representation
sourcemodel = fixinside(sourcemodel);

% keep the original details on inside and outside positions
originside = sourcemodel.inside;
origpos    = sourcemodel.pos;

% select only the dipole positions inside the brain for scanning
sourcemodel.pos    = sourcemodel.pos(originside,:);
sourcemodel.inside = true(size(sourcemodel.pos,1),1);

if hasmom
  sourcemodel.mom = sourcemodel.mom(:,originside);
end

if hasleadfield
  ft_info('using precomputed leadfields\n');
  sourcemodel.leadfield = sourcemodel.leadfield(originside);
else
  ft_info('computing forward model on the fly\n');
end

if ~isempty(cov)
  % compute signal and noise subspace from covariance matrix
  [u, s, v] = svd(cov);
else
  % compute signal and noise subspace from average data matrix
  [u, s, v] = svd(dat);
end
% select the noise subspace, c.f. equation 25
us = u(:,(numcomponent+1):end);
ps = us * us';

% allocate space to hold the result
estimate = [];
estimate.jr = nan(size(sourcemodel.pos,1),1);

ft_progress('init', feedback, 'scanning grid');
for i=1:size(sourcemodel.pos,1)
  ft_progress(i/size(sourcemodel.pos,1), 'scanning grid %d/%d\n', i, size(sourcemodel.pos,1));

  if hasleadfield && hasmom && size(sourcemodel.mom, 1)==size(sourcemodel.leadfield{i}, 2)
    % reuse the leadfield that was previously computed and project
    lf = sourcemodel.leadfield{i} * sourcemodel.mom(:,i);
  elseif  hasleadfield &&  hasmom
    % reuse the leadfield that was previously computed but don't project
    lf = sourcemodel.leadfield{i};
  elseif  hasleadfield && ~hasmom
    % reuse the leadfield that was previously computed
    lf = sourcemodel.leadfield{i};
  elseif ~hasleadfield &&  hasmom
    % compute the leadfield for a fixed dipole orientation
    lf = ft_compute_leadfield(sourcemodel.pos(i,:), sens, headmodel, leadfieldopt{:}) * sourcemodel.mom(:,i);
  else
    % compute the leadfield
    lf = ft_compute_leadfield(sourcemodel.pos(i,:), sens, headmodel, leadfieldopt{:});
  end

  % compute the music metric, c.f. equation 26
  estimate.jr(i) = (norm(ps * lf)./norm(lf)).^2;
  % as described in the Mosher 1992 paper on page 550, "...the general approach is to
  % evaluare Jr(i) over a fine three-dimensional grid, plot its inverse,
  % and look for p sharp spikes..."

end % for each dipole position
ft_progress('close');

% reassign the estimated values over the inside and outside grid positions
estimate.inside   = originside;
estimate.pos      = origpos;
if isfield(estimate, 'jr')
  estimate.jr( originside) = estimate.jr;
  estimate.jr(~originside) = nan;
end
if isfield(estimate, 'leadfield')
  estimate.leadfield( originside) = estimate.leadfield;
  estimate.leadfield(~originside) = {[]};
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
dataset_path = '/home/zhixuan/data1/HUMBI/HUMBI_uploaded/Gaze_81_140';
subject = 82;
frame = 1;
showAxis = false; % whether show 3D axis when plotting
useFitParams = false; % if true, read fit params and reconstruct face vertices/kps
                     % if false, read reconstructed face vertices/kps directly
vis3D = true; % if true, visualize 3D mesh and kps
               % if false, visualize their reprojections on image
%%%%%%%%%%%% modify above lines accordingly %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
config;
addpath('../face');
load('../face/model/model.mat', 'tri_list', 'shape_basis_3DMM', 'expression_basis_3DMM');
left_eye_tri_list = load_tri_list('./model/mesh_tri_left_eye.txt');
right_eye_tri_list = load_tri_list('./model/mesh_tri_right_eye.txt');
left_eye_vertex_indices = unique(left_eye_tri_list(:));
right_eye_vertex_indices = unique(right_eye_tri_list(:));

% extract camera params
intrinsics_path = sprintf(intrinsics_path_format, subject);
extrinsics_path = sprintf(extrinsics_path_format, subject);
[M, C, ~, ~] = extractCameraParameters(intrinsics_path, extrinsics_path, num_cam);
[ground_plane_normal, ground_plane_center] = getGroundPlaneFromCameras(C(:, 2:2:66), cameras_to_plane_offset);

% read data
reconstruction_dir = sprintf(recon_dir_format, subject, frame);
if useFitParams
    [~, ~, fitted_params] = read_face_recon(reconstruction_dir);
    [scale, R, t, shape_coefficients, expression_coefficients] = parse_fitted_params(fitted_params);
    vertices = reconstruct_face_fitted(meanface_3DMM, shape_basis_3DMM(:, 1:10), expression_basis_3DMM, ...
        shape_coefficients, expression_coefficients, scale, R, t);
else
    [vertices, ~, ~] = read_face_recon(reconstruction_dir);
end

% visualization
set(gcf, 'OuterPosition', get(0, 'Screensize')); % maximize the figure
if vis3D  % visualize 3D mesh and kps
    setup_vis; campos(C(:, 33)); camva(7); camlight('left');
    h = trimesh(tri_list+1, vertices(:,1), vertices(:,2), vertices(:,3), ...
        'FaceColor', [0 150 230]/255, 'EdgeColor', 'none', 'LineWidth', 0.1);
    h.FaceLighting = 'flat'; % or 'gouraud'
    hl = trimesh(left_eye_tri_list, vertices(:,1), vertices(:,2), vertices(:,3), ...
        'FaceColor', 'r', 'EdgeColor', 'none', 'LineWidth', 0.1);
    hl.FaceLighting = 'flat'; % or 'gouraud'
    hr = trimesh(right_eye_tri_list, vertices(:,1), vertices(:,2), vertices(:,3), ...
        'FaceColor', 'r', 'EdgeColor', 'none', 'LineWidth', 0.1);
    hr.FaceLighting = 'flat'; % or 'gouraud'
    title(sprintf('subject %d, frame %d', subject, frame));
else
    bboxes = readmatrix([sprintf(img_dir_format, subject, frame), '/list.txt']); % n x 5, bboxes on original image
    for i = 1:size(bboxes, 1)
        % compute reprojections
        cam = bboxes(i, 1);
        bbox = bboxes(i, 2:end); % 1 x 4, [xmin, xmax, ymin, ymax]
        scale_x = (bbox(2) - bbox(1) + 1) / img_size(1);
        scale_y = (bbox(4) - bbox(3) + 1) / img_size(2);
        % Notice: scale_x and scale_y may not be exactly the same but should be very close
        rvertices = reproject(vertices, M(:, :, cam+1)); % 3448 x 2
        rvertices = (rvertices - [bbox(1), bbox(3)]) ./ [scale_x, scale_y];
        rvertices = rvertices([left_eye_vertex_indices, right_eye_vertex_indices], :);

        % draw
        img_path = sprintf([img_dir_format, '/image%07d.png'], subject, frame, cam);
        img = imread(img_path);
        imshow(img); hold on
        scatter(rvertices(:,1), rvertices(:,2), 100, '.', 'MarkerEdgeColor', 'r', 'MarkerEdgeAlpha', 0.5);
        title(sprintf('subject %d, frame %d, cam %d (press "Enter" to go to next view)', subject, frame, cam));
        pause;
    end
end

function tri_list = load_tri_list(filename, startRow, endRow)
if nargin<=2
    startRow = 1;
    endRow = inf;
end
formatSpec = '%16f%16f%f%[^\n\r]';
fileID = fopen(filename,'r');
dataArray = textscan(fileID, formatSpec, endRow(1)-startRow(1)+1, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string', 'EmptyValue', NaN, 'HeaderLines', startRow(1)-1, 'ReturnOnError', false, 'EndOfLine', '\r\n');
for block=2:length(startRow)
    frewind(fileID);
    dataArrayBlock = textscan(fileID, formatSpec, endRow(block)-startRow(block)+1, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string', 'EmptyValue', NaN, 'HeaderLines', startRow(block)-1, 'ReturnOnError', false, 'EndOfLine', '\r\n');
    for col=1:length(dataArray)
        dataArray{col} = [dataArray{col};dataArrayBlock{col}];
    end
end
fclose(fileID);
tri_list = [dataArray{1:end-1}];
end

author: zhixuany

repo: HUMBI

save_path: github-repos/MATLAB/zhixuany-HUMBI


text:
%% OPT14_RUN
%
%  Modified:
%
%    08 January 2008
%
   %---------------------------------------------------------------------
   %  Nonlinear Equation Example
   %---------------------------------------------------------------------
   fprintf('---------------------------------------------------------\n')
   fprintf('Running testcase_14:  exact solutions (, )\n')
   fprintf('---------------------------------------------------------\n')
   fname = 'opt14_fgh';
   options = [];
   options.verbose            = 0;
   options.max_iterations     = 180;
   options.max_fevals         = 180;
   options.method             = 'newton';

   fprintf('Newton:\n')
   options.globalization      = 'trust_region';

   x0 = [ -.8; .8; -.8 ];
   x = entrust(fname, x0, options);
   fprintf('Newtons method produced  (%10.7e,%10.7e,%10.7e)\n\n',x(1),x(2),x(3))
   f = opt14_fgh ( x, 'f' );
   fprintf('Value of F(X) = %f\n', f );

   x0 = [ 2; -1; 1 ];
   x = entrust(fname, x0, options);
   fprintf('Newtons method produced  (%10.7e,%10.7e,%10.7e)\n\n',x(1),x(2),x(3))
   f = opt14_fgh ( x, 'f' );
   fprintf('Value of F(X) = %f\n', f );

   %---------------------------------------------------------------------
   %  Test Gauss-Newton strategies.
   %---------------------------------------------------------------------
   fprintf('---------------------------------------------------------\n')
   fprintf('Running testcase_14 as least squares problem: \n')
   fprintf('Exact solution (,)\n')
   fprintf('---------------------------------------------------------\n')
   fname = 'opt14_rj';
   options = [];
   options.verbose            = 0;
   options.method             = 'gauss_newton';
   options.step_tolerance     = 1.e-15;
   options.globalization      = 'none';
   options.gradient_tolerance = 1.e-10;
   options.max_iterations     = 450;

   x0 = [ -.8; .8; -.8 ];
   x = entrust(fname, x0, options);
   fprintf('Gauss-Newton produced  (%10.7e, %10.7e, %10.7e)\n\n',x(1),x(2),x(3))
   [ res, jac ] = opt14_rj ( x, 'f' );
   fprintf('Norm of RES(X) = %f\n', norm ( res ) );

   x0 = [ 2; -1; 1 ];
   x = entrust(fname, x0, options);
   fprintf('Gauss-Newton produced  (%10.7e, %10.7e, %10.7e)\n\n',x(1),x(2),x(3))
   [ res, jac ] = opt14_rj ( x, 'f' );
   fprintf('Norm of RES(X) = %f\n', norm ( res ) );


text: %MLDIVIDE Dataset overload

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
% function [up,n1,n2] = p2up(n)
%
% Calculates next power of 2, and left/right padding to center the
% original n locations.
%
% Input:
%   n: non-dyadic integer
% Output:
%   up: next power of 2
%   n1: length on left
%   n2: length on right
%
%---------------------------------------------------------------------------------
%    Synchrosqueezing Toolbox
%    Authors: Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/)
%---------------------------------------------------------------------------------
function [up,n1,n2] = p2up(n)
    up = 2^(1+round(log2(n+eps)));
    n1 = floor((up-n)/2); n2 = n1;
    if (mod(2*n1+n,2)==1), n2 = n1 + 1; end
end

author: HeLiangHIT


text:
function [parm_success] = get_parm_xml(xmlfile,parm,parm_key)
% function which will read xml file and updates the parameter_field matching
% the requested parameter. Will return "FAIL" or "SUCCESS" based on
% whether it found the variable
%
% Bekaert David - Jet Propulsion Laboratory

if nargin<3
    error('Require three inputs')
end

% initialize the field
temp = 0;
count_changes = 0;

% loop over the xml file
fid = fopen(xmlfile,'r');
fid_new = fopen([xmlfile '_temp'],'w');

while temp~=-1
    temp = fgetl(fid);
    if temp==-1
        break
    end
    ix =  strfind(temp,['"' parm '"']);
    if isempty(ix)
        ix = strfind(temp,['"' lower(parm) '"']);
    end
    if ~isempty(ix)
        % see if the value is at this string or not
        ix_value = strfind(temp,'<value>');
        while isempty(ix_value);
            % write the output in case its not a value line
            if isempty(ix_value);
                fprintf(fid_new,[temp '\n']);
            end

            temp = fgetl(fid);
            ix_value = strfind(temp,'<value>');

            % avoid inf loop
            if temp==-1
                break
            end
        end

        if temp~=-1
            % update with new value or string
            if isnumeric(parm_key)
                parm_key = strtrim(num2str(parm_key));
            end
            parm_key_str = [repmat(' ',1,ix_value-1) '<value>' parm_key '</value>'];
            fprintf(fid_new,[parm_key_str '\n']);
            count_changes = count_changes+1;
        end
    else
        fprintf(fid_new,[temp '\n']);
    end
end
fclose(fid);
fclose(fid_new);

% rename the file
delete(xmlfile)
movefile([xmlfile '_temp'],xmlfile)

% see if the search was succesful
if count_changes==0
    fprintf(['Could not find ' parm '... \n'])
else
    parm_success = 'SUCCESS';
end

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
%{
    Copyright (c) 2014, Philipp Krähenbühl
    All rights reserved.

    Redistribution and use in source and binary forms, with or without
    modification, are permitted provided that the following conditions are met:
        * Redistributions of source code must retain the above copyright
        notice, this list of conditions and the following disclaimer.
        * Redistributions in binary form must reproduce the above copyright
        notice, this list of conditions and the following disclaimer in the
        documentation and/or other materials provided with the distribution.
        * Neither the name of the Stanford University nor the
        names of its contributors may be used to endorse or promote products
        derived from this software without specific prior written permission.

    THIS SOFTWARE IS PROVIDED BY Philipp Krähenbühl ''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 Philipp Krähenbühl 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.
%}
addpath( genpath('../build') )

author: Cloud-CV

repo: object-proposals


text:
function cdf = lorentz_cdf ( x )

%*****************************************************************************80
%
%% LORENTZ_CDF evaluates the Lorentz CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 1999
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the CDF.
%
%    Output, real CDF, the value of the CDF.
%
  cdf = 0.5 + atan ( x ) / pi;

  return
end


text:
classdef TabPanel < uix.TabPanel
    %TabPanel  Show one element inside a tabbed panel
    %
    %   obj = uiextras.TabPanel() creates a panel with tabs along one edge
    %   to allow selection between the different child objects contained.
    %
    %   obj = uiextras.TabPanel(param,value,...) also sets one or more
    %   property values.
    %
    %   See the documentation for more detail and the list of properties.
    %
    %   Examples:
    %   >> f = figure();
    %   >> p = uiextras.TabPanel( 'Parent', f, 'Padding', 5 );
    %   >> uicontrol( 'Style', 'frame', 'Parent', p, 'Background', 'r' );
    %   >> uicontrol( 'Style', 'frame', 'Parent', p, 'Background', 'b' );
    %   >> uicontrol( 'Style', 'frame', 'Parent', p, 'Background', 'g' );
    %   >> p.TabNames = {'Red', 'Blue', 'Green'};
    %   >> p.SelectedChild = 2;
    %
    %   See also: uiextras.Panel
    %             uiextras.BoxPanel

    %  Copyright 2009-2014 The MathWorks, Inc.
    %  $Revision: 979 $ $Date: 2014-09-28 14:26:12 -0400 (Sun, 28 Sep 2014) $

    properties( Hidden, Access = public, Dependent )
        Callback
    end

    properties( Access = private )
        Callback_ = '' % backing for Callback
    end

    properties( Hidden, Access = public, Dependent )
        Enable % deprecated
        SelectedChild
        TabEnable
        TabNames
        TabPosition
        TabSize
    end

    properties( Access = private )
        SelectionChangedListener % listener
    end

    methods

        function obj = TabPanel( varargin )

            % Call uix constructor
            obj@uix.TabPanel( varargin{:} )

            % Auto-parent
            if ~ismember( 'Parent', varargin(1:2:end) )
                obj.Parent = gcf();
            end

            % Create listeners
            selectionChangedListener = event.listener( obj, ...
                'SelectionChanged', @obj.onSelectionChanged );

            % Store properties
            obj.SelectionChangedListener = selectionChangedListener;

        end % constructor

    end % structors

    methods

        function value = get.Enable( ~ )

            % Warn
            % warning( 'uiextras:Deprecated', ...
            %     'Property ''Enable'' will be removed in a future release.' )

            % Return
            value = 'on';

        end % get.Enable

        function set.Enable( ~, value )

            % Check
            assert( ischar( value ) && any( strcmp( value, {'on','off'} ) ), ...
                'uiextras:InvalidPropertyValue', ...
                'Property ''Enable'' must be ''on'' or ''off''.' )

            % Warn
            % warning( 'uiextras:Deprecated', ...
            %     'Property ''Enable'' will be removed in a future release.' )

        end % set.Enable

        function value = get.Callback( obj )

            % Get
            value = obj.Callback_;

        end % get.Callback

        function set.Callback( obj, value )

            % Check
            if ischar( value ) % string
                % OK
            elseif isa( value, 'function_handle' ) && ...
                    isequal( size( value ), [1 1] ) % function handle
                % OK
            elseif iscell( value ) && ndims( value ) == 2 && ...
                    size( value, 1 ) == 1 && size( value, 2 ) > 0 && ...
                    isa( value{1}, 'function_handle' ) && ...
                    isequal( size( value{1} ), [1 1] ) %#ok<ISMAT> % cell callback
                % OK
            else
                error( 'uiextras:InvalidPropertyValue', ...
                    'Property ''Callback'' must be a valid callback.' )
            end

            % Set
            obj.Callback_ = value;

        end % set.Callback

        function value = get.SelectedChild( obj )

            % Get
            value = obj.Selection;

        end % get.SelectedChild

        function set.SelectedChild( obj, value )

            % Set
            obj.Selection = value;

        end % set.SelectedChild

        function value = get.TabEnable( obj )

            % Get
            value = transpose( obj.TabEnables );

        end % get.TabEnable

        function set.TabEnable( obj, value )

            % Set
            obj.TabEnables = value;

        end % set.TabEnable

        function value = get.TabNames( obj )

            % Get
            value = transpose( obj.TabTitles );

        end % get.TabNames

        function set.TabNames( obj, value )

            % Set
            obj.TabTitles = value;

        end % set.TabNames

        function value = get.TabPosition( obj )

            % Get
            value = obj.TabLocation;

        end % get.TabPosition

        function set.TabPosition( obj, value )

            % Set
            obj.TabLocation = value;

        end % set.TabPosition

        function value = get.TabSize( obj )

            % Get
            value = obj.TabWidth;

        end % get.TabSize

        function set.TabSize( obj, value )

            % Set
            obj.TabWidth = value;

        end % set.TabSize

    end % accessors

    methods( Access = private )

        function onSelectionChanged( obj, source, eventData )

            % Create legacy event data structure
            oldEventData = struct( 'Source', eventData.Source, ...
                'PreviousChild', eventData.OldValue, ...
                'SelectedChild', eventData.NewValue );

            % Call callback
            callback = obj.Callback_;
            if ischar( callback ) && isequal( callback, '' )
                % do nothing
            elseif ischar( callback )
                feval( callback, source, oldEventData )
            elseif isa( callback, 'function_handle' )
                callback( source, oldEventData )
            elseif iscell( callback )
                feval( callback{1}, source, oldEventData, callback{2:end} )
            end

        end % onSelectionChanged

    end % event handlers

end % classdef

author: goGPS-Project


text:
function [ s ] = xml2struct( file )
%Convert xml file into a MATLAB structure
% [ s ] = xml2struct( file )
%
% A file containing:
% <XMLname attrib1="Some value">
%   <Element>Some text</Element>
%   <DifferentElement attrib2="2">Some more text</DifferentElement>
%   <DifferentElement attrib3="2" attrib4="1">Even more text</DifferentElement>
% </XMLname>
%
% Used to produce:
% s.XMLname.Attributes.attrib1 = "Some value";
% s.XMLname.Element.Text = "Some text";
% s.XMLname.DifferentElement{1}.Attributes.attrib2 = "2";
% s.XMLname.DifferentElement{1}.Text = "Some more text";
% s.XMLname.DifferentElement{2}.Attributes.attrib3 = "2";
% s.XMLname.DifferentElement{2}.Attributes.attrib4 = "1";
% s.XMLname.DifferentElement{2}.Text = "Even more text";
%
% Will produce (gp: to matche the output of xml2struct in XML4MAT, but note that Element(2) is empty):
% Element: Some text
% DifferentElement:
%    attrib2: 2
%    DifferentElement: Some more text
% attrib1: Some value
%
% Element:
% DifferentElement:
%    attrib3: 2
%    attrib4: 1
%    DifferentElement: Even more text
% attrib1:
%
% Note the characters : - and . are not supported in structure fieldnames and
% are replaced by _
%
% Written by W. Falkena, ASTI, TUDelft, 21-08-2010
% Attribute parsing speed increased by 40% by A. Wanner, 14-6-2011
% 2011/12/14 giopia: changes in the main function to make more similar to xml2struct of the XML4MAT toolbox, bc it's used by fieldtrip
% 2012/04/04 roboos: added the original license clause, see also http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=645#c11
% 2012/04/04 roboos: don't print the filename that is being read

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

    if (nargin < 1)
        clc;
        help xml2struct
        return
    end

    %check for existance
    if (exist(file,'file') == 0)
        %Perhaps the xml extension was omitted from the file name. Add the
        %extension and try again.
        if (isempty(strfind(file,'.xml')))
            file = [file '.xml'];
        end

        if (exist(file,'file') == 0)
            ft_error(['The file ' file ' could not be found']);
        end
    end
    %fprintf('xml2struct reading %s\n', file); % gp 11/12/15
    %read the xml file
    xDoc = xmlread(file);

    %parse xDoc into a MATLAB structure
    s = parseChildNodes(xDoc);
    fn = fieldnames(s); s = s.(fn{1}); % gp 11/12/15: output is compatible with xml2struct of xml2mat
end

% ----- Subfunction parseChildNodes -----
function [children,ptext] = parseChildNodes(theNode)
    % Recurse over node children.
    children = struct;
    ptext = [];
    if theNode.hasChildNodes
        childNodes = theNode.getChildNodes;
        numChildNodes = childNodes.getLength;

        for count = 1:numChildNodes
            theChild = childNodes.item(count-1);
            [text,name,attr,childs] = getNodeData(theChild);

            if (~strcmp(name,'#text') && ~strcmp(name,'#comment'))
                %XML allows the same elements to be defined multiple times,
                %put each in a different cell
                if (isfield(children,name))
%                     if 0 % numel(children) > 1 % gp 11/12/15: (~iscell(children.(name)))
%                         %put existsing element into cell format
%                         children.(name) = {children.(name)};
%                     end
                    index = length(children)+1; % gp 11/12/15: index = length(children.(name))+1;
                else
                  index = 1; % gp 11/12/15: new field
                end
                %add new element
                children(index).(name) = childs;
                if isempty(attr)
                  if(~isempty(text))
                    children(index).(name) = text;
                  end
                else
                  fn = fieldnames(attr);
                  for f = 1:numel(fn)
                    children(index).(name)(1).(fn{f}) = attr.(fn{f}); % gp 11/12/15: children.(name){index}.('Attributes') = attr;
                  end
                  if(~isempty(text))
                    children(index).(name).(name) = text; % gp 11/12/15: children.(name){index}.('Text') = text;
                  end
                end

%                 else % gp 11/12/15: cleaner code, don't reuse the same code
%                     %add previously unknown new element to the structure
%                     children.(name) = childs;
%                     if(~isempty(text))
%                         children.(name) = text; % gp 11/12/15: children.(name).('Text') = text;
%                     end
%                     if(~isempty(attr))
%                         children.('Attributes') = attr; % gp 11/12/15 children.(name).('Attributes') = attr;
%                     end
%                 end
            elseif (strcmp(name,'#text'))
                %this is the text in an element (i.e. the parentNode)
                if (~isempty(regexprep(text,'[\s]*','')))
                    if (isempty(ptext))
                        ptext = text;
                    else
                        %what to do when element data is as follows:
                        %<element>Text <!--Comment--> More text</element>

                        %put the text in different cells:
                        % if (~iscell(ptext)) ptext = {ptext}; end
                        % ptext{length(ptext)+1} = text;

                        %just append the text
                        ptext = [ptext text];
                    end
                end
            end
        end
    end
end

% ----- Subfunction getNodeData -----
function [text,name,attr,childs] = getNodeData(theNode)
    % Create structure of node info.

    %make sure name is allowed as structure name
    name = regexprep(char(theNode.getNodeName),'[-:.]','_');

    attr = parseAttributes(theNode);
    if (isempty(fieldnames(attr)))
        attr = [];
    end

    %parse child nodes
    [childs,text] = parseChildNodes(theNode);

    if (isempty(fieldnames(childs)))
        %get the data of any childless nodes
        try
            %faster then if any(strcmp(methods(theNode), 'getData'))
            text = char(theNode.getData);
        catch
            %no data
        end
    end

end

% ----- Subfunction parseAttributes -----
function attributes = parseAttributes(theNode)
    % Create attributes structure.

    attributes = struct;
    if theNode.hasAttributes
       theAttributes = theNode.getAttributes;
       numAttributes = theAttributes.getLength;

       for count = 1:numAttributes
            %attrib = theAttributes.item(count-1);
            %attr_name = regexprep(char(attrib.getName),'[-:.]','_');
            %attributes.(attr_name) = char(attrib.getValue);

            %Suggestion of Adrian Wanner
            str = theAttributes.item(count-1).toString.toCharArray()';
            k = strfind(str,'=');
            attr_name = regexprep(str(1:(k(1)-1)),'[-:.]','_');
            attributes.(attr_name) = str((k(1)+2):(end-1));
       end
    end
end


text:
function [dStates, A, B, Positions] = Derive_EoM()

%This function derives the equations of motion for a tractor trailer truck.

%The outputs assume that the states and inputs are being listed in the
%order in which they are defined

%% Parameters

Lt = sym('Lt');
Lc = sym('Lc');

%% Inputs

v = sym('v');
psi = sym('psi');

%% States

x = sym('x');
y = sym('y');
th = sym('th');
phi = sym('phi');

%% State Derivatives

Dx = sym('Dx');
Dy = sym('Dy');
Dth = sym('Dth');
Dphi = sym('Dphi');

%% Equations of motion

%Unit Vectors
    i = [1;0];
    j = [0;1];

    An = cos(th)*i + sin(th)*j;
    At = -sin(th)*i + cos(th)*j;
    Bn = cos(phi)*An + sin(phi)*At;
    Bt = -sin(phi)*An + cos(phi)*At;
    Cn = cos(psi)*Bn + sin(psi)*Bt;
    Ct = -sin(psi)*Bn + cos(psi)*Bt;

%Derivatives of Unit Vectors
    dAn = Dth*At;
    dAt = -Dth*An;
    dBn = (Dphi + Dth)*Bt;
    dBt = -(Dphi + Dth)*Bn;

%Position Vectors
    Pa = x*i + y*j;
    Pb = Pa + Lt*At;
    Pc = Pb + Lc*Bt;
    Positions = [Pa,Pb,Pc];

%Derivatives of Position Vectors
    dPa = Dx*i + Dy*j;
    dPb = dPa + Lt*dAt;
    dPc = dPb + Lc*dBt;

%Trailer Back Wheel Constraint
    RHS_1 = 0;
    LHS_1 = Dot(dPa,An);

%Cab Back Wheel Constraint
    RHS_2 = 0;
    LHS_2 = Dot(dPb,Bn);

%Cab Front Wheel Constraints
    RHS_3 = dPc(1);
    LHS_3 = v*Ct(1);

    RHS_4 = dPc(2);
    LHS_4 = v*Ct(2);

%Solve for State Derivatives
    Solution = solve(...
        RHS_1-LHS_1, RHS_2-LHS_2, RHS_3-LHS_3, RHS_4-LHS_4,...
        Dx, Dy, Dth, Dphi);

%Express as Vector [Dx, Dy, Dth, Dphi]

    dStates = [...
        Solution.Dx;
        Solution.Dy;
        Solution.Dth;
        Solution.Dphi];

%% Jacobian Calculations

    J_x = diff(dStates,x);
    J_y = diff(dStates,y);
    J_th = diff(dStates,th);
    J_phi = diff(dStates,phi);
    J_v = diff(dStates,v);
    J_psi = diff(dStates,psi);

    A = [J_x, J_y, J_th, J_phi];
    B = [J_v, J_psi];

end

author: MatthewPeterKelly

repo: dscTutorials


text:
classdef MimImageCoordinateUtilities
    % MimImageCoordinateUtilities. Utility functions related to processing 3D
    % image coordinates
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    methods (Static)

        % In Matlab, matrices can be represented as a linear vector. So every
        % point in a 3D matrix has a linear index as well as an i-j-k
        % coordinate. This function returns the 'offset' values which can be
        % added to the linear index of any point to return the linear indices of
        % its nearest neighbours.
        function [linear_offsets, linear_offsets27] = GetLinearOffsets(image_size)
            % Compute linear index offsets for diretion vectors
            dirs = [5, 23, 11, 17, 13, 15];
            linear_offsets = MimImageCoordinateUtilities.GetLinearOffsetsForDirections(dirs, image_size);

            dirs = 1:27;
            linear_offsets27 = MimImageCoordinateUtilities.GetLinearOffsetsForDirections(dirs, image_size);
        end

        function linear_offsets = GetLinearOffsetsForDirections(dirs, image_size)
            direction_vectors = MimImageCoordinateUtilities.CalculateDirectionVectors;
            linear_offsets = zeros(1, numel(dirs));
            for n = 1 : length(dirs)
                direction = dirs(n);
                direction_vector = direction_vectors(direction, :);
                start_point = [2 2 2];
                i_end_point = start_point + direction_vector;
                i = [start_point(1); i_end_point(1)];
                j = [start_point(2); i_end_point(2)];
                k = [start_point(3); i_end_point(3)];
                linear_indices = sub2ind(image_size, i, j, k);
                linear_offsets(n) = linear_indices(2) - linear_indices(1);
            end
        end

        % Returns the coordinates of each point in a 3x3x3 matrix relative to
        % its centre
        function direction_vectors = CalculateDirectionVectors
            [i, j, k] = ind2sub([3 3 3], 1:27);
            direction_vectors = [i' - 2, j' - 2, k' - 2];
        end

        % This function alters matrix indices to transform from a smaller matrix to
        % a bigger one
        function new_indices = OffsetIndices(indices, offset, size_small, size_big)
            indices_i = (indices - 1);
            div1 = (size_small(1));
            div2 = (size_small(1)*size_small(2));

            k_mod = rem(indices_i, div2);

            k = (indices_i - k_mod)/div2; % Equivalent to idivide but quicker

            i = mod(k_mod, div1);

            j = (k_mod - i)/div1; % Equivalent to idivide but quicker

            i = i + offset(1);
            j = j + offset(2);
            k = k + offset(3);

            % Note that i,j,k can be negative, in the case where a border has
            % been added to the image so its coordinates extend beyond the
            % boundaries of the original image

            new_indices = 1 + (i) + (j)*size_big(1) + (k)*size_big(1)*size_big(2);
        end

        % Creates an image cropped to the smallest box size that encloses all
        % the points specified by their linear indices.
        function [offset reduced_image reduced_image_size] = GetMinimalImageForIndices(indices, image_size)
            if size(indices, 1) > 1
                error('GetMinimalImageForIndices requires indices to be in a row vector');
            end
            indices = int32(indices);
            [i, j, k] = MimImageCoordinateUtilities.FastInd2sub(image_size, indices);

            voxel_coordinates = [i' j' k'];
            mins = min(voxel_coordinates, [], 1);
            maxs = max(voxel_coordinates, [], 1);
            reduced_image_size = maxs - mins + int32([1 1 1]);
            reduced_image = false(reduced_image_size);
            offset = mins - 1;
            i = MimImageCoordinateUtilities.FastSub2ind(reduced_image_size, voxel_coordinates(:,1)-offset(1), voxel_coordinates(:,2)-offset(2), voxel_coordinates(:,3)-offset(3));

            reduced_image(i) = true;
        end

        % A faster alternative to Ind2sub
        function [i, j, k] = FastInd2sub(im_size, indices)
            indices = indices - 1;
            div1 = (im_size(1));
            div2 = ((im_size(1)*im_size(2)));

            k_mod = rem(indices, div2);

            k = 1 + (indices - k_mod)/div2; % Equivalent to idivide but quicker
            i = 1 + mod(k_mod, im_size(1));
            j = 1 + (k_mod - i + 1)/div1; % Equivalent to idivide but quicker
        end

        function indices = FastSub2ind(im_size, i, j, k)
            indices = i + (j-1)*im_size(1) + (k-1)*im_size(1)*im_size(2);
        end

        function rot_matrix = GetEulerRotationMatrix(phi, theta, psi)
            rot_matrix = zeros(3,3);

            rot_matrix(1, 1) = cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi);
            rot_matrix(1, 2) = cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi);
            rot_matrix(1, 3) = sin(psi)*sin(theta);

            rot_matrix(2, 1) = -sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi);
            rot_matrix(2, 2) = -sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi);
            rot_matrix(2, 3) = cos(psi)*sin(theta);

            rot_matrix(3, 1) =  sin(theta)*sin(phi);
            rot_matrix(3, 2) = -sin(theta)*cos(psi);
            rot_matrix(3, 3) = cos(theta);
        end

        function affine_matrix = CreateAffineMatrix(x)
            affine_matrix = zeros(3, 4, 'single');
            affine_matrix(:) = x(:);
            affine_matrix = [affine_matrix; [0 0 0 1]];
        end

        function affine_matrix = CreateAffineTranslationMatrix(x)
            affine_matrix = zeros(3, 4, 'single');
            affine_matrix(1, 1) = 1;
            affine_matrix(2, 2) = 1;
            affine_matrix(3, 3) = 1;
            affine_matrix(4, 4) = 1;
            affine_matrix(1:3, 4) = x;
        end

        function affine_matrix = CreateRigidAffineMatrix(x)
            affine_matrix = zeros(3, 4, 'single');

            euler_rot_matrix = MimImageCoordinateUtilities.GetEulerRotationMatrix(x(1), x(2), x(3));
            affine_matrix(1:3, 1:3) = euler_rot_matrix;
            affine_matrix(1:3, 4) = x(4:6);

            affine_matrix = [affine_matrix; [0 0 0 1]];
        end

        function [i, j, k] = TransformCoordsAffine(i, j, k, augmented_matrix)
            [j, i, k] = MimImageCoordinateUtilities.TranslateAndRotateMeshGrid(j, i, k, augmented_matrix(1:3,1:3), augmented_matrix(1:3,4));
        end

        function [i, j, k] = TransformCoordsFluid(i, j, k, deformation_field)
            i = i - deformation_field.RawImage(:,:,:,1);
            j = j - deformation_field.RawImage(:,:,:,2);
            k = k - deformation_field.RawImage(:,:,:,3);
        end

        function [X, Y, Z] = TranslateAndRotateMeshGrid(X, Y, Z, rot_matrix, trans_matrix)
            % Rotates and translates meshgrid generated coordinates in 3D
            % Note coordinates are [XYZ] NOT [IJK]
            [X, Y, Z] = MimImageCoordinateUtilities.RotateMeshGrid(X + trans_matrix(1), Y + trans_matrix(2), Z + trans_matrix(3), rot_matrix);
        end

        function [X, Y, Z] = RotateMeshGrid(X, Y, Z, rot_matrix)
            % Rotates coordinates that are given in 3D meshgrid matrices
            coords = rot_matrix * [ ...
                reshape(X, 1, numel(X)); ...
                reshape(Y, 1, numel(Y)); ...
                reshape(Z, 1, numel(Z)) ...
                ];

            X = reshape(coords(1, :), size(X));
            Y = reshape(coords(2, :), size(Y));
            Z = reshape(coords(3, :), size(Z));
        end

        function affine_matrix = GetAffineTranslationFromPatientPosition(image_1, image_2)

            % Get the coordinates of the centre of the first voxel in image1,
            % relative to the centre of image1
            [i1, j1, k1] = image_1.GlobalCoordinatesToCoordinatesMm([1, 1, image_1.OriginalImageSize(3)]);
            [i1, j1, k1] = image_1.GlobalCoordinatesMmToCentredGlobalCoordinatesMm(i1, j1, k1);
            image_1_origin_coordinates = [i1, j1, k1];
            image_1_centre_coordinates = image_1.GlobalOrigin - image_1_origin_coordinates;

            % Get the coordinates of the centre of the first voxel in image2,
            % relative to centre of image2
            [i2, j2, k2] = image_2.GlobalCoordinatesToCoordinatesMm([1, 1, image_2.OriginalImageSize(3)]);
            [i2, j2, k2] = image_2.GlobalCoordinatesMmToCentredGlobalCoordinatesMm(i2, j2, k2);
            image_2_origin_coordinates = [i2, j2, k2];
            image_2_centre_coordinates = image_2.GlobalOrigin - image_2_origin_coordinates;

            translation = image_1_centre_coordinates - image_2_centre_coordinates;

            translation = translation([2 1 3]);
            translation(3) = - translation(3);
            affine_matrix = MimImageCoordinateUtilities.CreateAffineTranslationMatrix(translation);
        end

        function deformation_field = AdjustDeformationFieldForInitialAffineTransformation(deformation_field, affine_initial_matrix)
            % To combine a rigid transformation with a nonrigid deformation field,
            % compute the change in image coordinates after applying the
            % deformation field and then the rigid affine transformation.

            [df_i, df_j, df_k] = deformation_field.GetCentredGlobalCoordinatesMm;
            [df_i, df_j, df_k] = ndgrid(df_i, df_j, df_k);
            [df_i_t, df_j_t, df_k_t] = MimImageCoordinateUtilities.TransformCoordsFluid(df_i, df_j, df_k, deformation_field);
            [df_i_t, df_j_t, df_k_t] = MimImageCoordinateUtilities.TransformCoordsAffine(df_i_t, df_j_t, df_k_t, affine_initial_matrix);

            deformation_field_raw = zeros(deformation_field.ImageSize);
            deformation_field_raw(:,:,:,1) = df_i - df_i_t;
            deformation_field_raw(:,:,:,2) = df_j - df_j_t;
            deformation_field_raw(:,:,:,3) = df_k - df_k_t;
            deformation_field2 = deformation_field.BlankCopy;
            deformation_field2.ChangeRawImage(deformation_field_raw);
            deformation_field = deformation_field2;
        end

        function [permutation_vector, flip] = GetDimensionPermutationVectorFromDicomOrientation(orientation, reporting)
            % Returns a vector which defines the order in which the dimensions of an
            % DICOM image volume should be permuted in order to align it with the
            % PTK coordinate system

            % DICOM coordinates are XYZ but Matlab's coordinates are YXZ. We
            % convert the orientation to Matlab coordinates, which we refer to as IJK.
            orientation_1 = orientation([5, 4, 6])'; % Direction of first image axis in ijk (=yxz) coordinates
            orientation_2 = orientation([2, 1, 3])'; % Direction of second image axis in ijk (=yxz) coordinates

            % By switching the i and j axes we have inverted the coordinate
            % system, so we need to flip the k dimension
            orientation_1(3) = - orientation_1(3);
            orientation_2(3) = - orientation_2(3);

            % Determine the PTK dimensions to which each of these vectors correspond
            [permutation_vector, dimension_1, dimension_2, dimension_3] = MimImageCoordinateUtilities.GetPermutationFromOrientations(orientation_1, orientation_2, reporting);

            % Calculate flip for each dimension, based on whether the
            % dimension axis lies in the same or opposite direction to the
            % image axis
            orientation_3 = cross(orientation_2, orientation_1);
            flip = MimImageCoordinateUtilities.GetFlipFromOrientations(orientation_1, orientation_2, orientation_3, dimension_1, dimension_2, dimension_3);

            % Check the resulting vector is valid
            if (sum(permutation_vector == 1) ~= 1) || (sum(permutation_vector == 2) ~= 1) || (sum(permutation_vector == 3) ~= 1) || ...
                    ~isempty(setdiff(permutation_vector, [1,2,3]))
                reporting.Error('MimImageCoordinateUtilities:InvalidPermutationVector', 'GetDimensionPermutationVectorFromDicomOrientation() resulted in an invalid permutation vector');
            end
        end

        function [permutation_vector, flip] = GetDimensionPermutationVectorFromNiiOrientation(header, reporting)
            % Returns a vector which defines the order in which the dimensions of a
            % nii image volume should be permuted in order to align it with the
            % PTK coordinate system

            if header.SformCode > 0
                d1 = [-1; -1; 1].*header.SrowX(1:3);
                d2 = [-1; -1; 1].*header.SrowY(1:3);
                [permutation_vector, flip] = MimImageCoordinateUtilities.GetDimensionPermutationVectorFromDicomOrientation([d1; d2], reporting);
            elseif header.QformCode > 0
                B = header.QuaternB;
                C = header.QuaternC;
                D = header.QuaternD;
                A = sqrt(1 - B^2 - C^2 - D^2);
                d1 = [-1; -1; 1].*[A^2+B^2-C^2-D^2; 2*(B*C - A*D); 2*(B*D+A*C)];
                d2 = [-1; -1; 1].*[2*(B*C + A*D); A^2 + C^2 - B^2 - D^2; 2*(C*D - A*B)];
                d3 = [-1; -1; 1].*[2*(B*D - A*C); 2*(C*D + A*B); A^2 + D^2 - B^2 - C^2];
                [permutation_vector, flip] = MimImageCoordinateUtilities.GetDimensionPermutationVectorFromDicomOrientation([d1; d2], reporting);
            else
                d1 = [-1; -1; 1].*[1; 0; 0];
                d2 = [-1; -1; 1].*[0; 1; 0];
                [permutation_vector, flip] = MimImageCoordinateUtilities.GetDimensionPermutationVectorFromDicomOrientation([d1; d2], reporting);
            end
        end

        function [permutation_vector, flip] = GetDimensionPermutationVectorForAnalyze(orientation, reporting)
            permutation_vector = [2, 1, 3];
            flip = [false, false, true];
        end

        function [permutation_vector, dimension_1, dimension_2, dimension_3] = GetPermutationFromOrientations(orientation_1, orientation_2, reporting)
            [dimension_1, dimension_2, dimension_3] = MimImageCoordinateUtilities.GetDimensionIndicesFromOrientations(orientation_1, orientation_2, reporting);

            permutation_vector = [3, 3, 3];
            permutation_vector(dimension_1) = 1;
            permutation_vector(dimension_2) = 2;
        end

        function flip = GetFlipFromOrientations(orientation_1, orientation_2, orientation_3, dimension_1, dimension_2, dimension_3)
            flip = [false, false, false];
            flip(dimension_1) = orientation_1(dimension_1) < 0;
            flip(dimension_2) = orientation_2(dimension_2) < 0;
            flip(dimension_3) = orientation_3(dimension_3) < 0;
        end

        function dicom_cosine = AnatomicalOrientationToDicomCosine(anatomical_orientation_char, reporting)
            switch anatomical_orientation_char
                case 'R'
                    dicom_cosine = [1, 0, 0];
                case 'L'
                    dicom_cosine = [-1, 0, 0];
                case 'A'
                    dicom_cosine = [0, 1, 0];
                case 'P'
                    dicom_cosine = [0, -1, 0];
                case 'I'
                    dicom_cosine = [0, 0, 1];
                case 'S'
                    dicom_cosine = [0, 0, -1];
                otherwise
                    reporting.Error('MimImageCoordinateUtilities:UnknownAnatomicalOrientation', ['MimImageCoordinateUtilities: WARNING: no implementation yet for anatomical orientation ' anatomical_orientation_char '.']);
            end
        end

        function [permutation_vector, flip_orientation] = GetDimensionPermutationVectorFromAnatomicalOrientation(anatomical_orientation_string, reporting)
            direction_cosine_1 = MimImageCoordinateUtilities.AnatomicalOrientationToDicomCosine(anatomical_orientation_string(1), reporting);
            direction_cosine_2 = MimImageCoordinateUtilities.AnatomicalOrientationToDicomCosine(anatomical_orientation_string(2), reporting);
            direction_cosine_3 = MimImageCoordinateUtilities.AnatomicalOrientationToDicomCosine(anatomical_orientation_string(3), reporting);

            [permutation_vector, flip_orientation] = MimImageCoordinateUtilities.GetDimensionPermutationVectorFromMhdCosines(direction_cosine_1, direction_cosine_2, direction_cosine_3, reporting);
        end

        function [permutation_vector, flip_orientation] = GetDimensionPermutationVectorFromMhdCosines(direction_cosine_1, direction_cosine_2, direction_cosine_3, reporting)
            % Determines the permutation and flip transformations required
            % to convert MetaIO data to PTK data

            % DICOM coordinates are XYZ but Matlab's coordinates are YXZ. We
            % convert the orientation to Matlab coordinates, which we refer to as IJK.
            orientation_1 = direction_cosine_2([2, 1, 3])'; % Direction of first image axis in ijk (=yxz) coordinates
            orientation_2 = direction_cosine_1([2, 1, 3])'; % Direction of second image axis in ijk (=yxz) coordinates
            orientation_3 = direction_cosine_3([2, 1, 3])'; % Direction of third image axis in ijk (=yxz) coordinates

            % By switching the i and j axes we have inverted the coordinate
            % system, so we need to flip the k dimension
            orientation_1(3) = - orientation_1(3);
            orientation_2(3) = - orientation_2(3);
            orientation_3(3) = - orientation_3(3);

            % Determine the PTK dimensions to which each of these vectors correspond
            [permutation_vector, dimension_1, dimension_2, dimension_3] = MimImageCoordinateUtilities.GetPermutationFromOrientations(orientation_1, orientation_2, reporting);

            % Calculate flip for each dimension, based on whether the
            % dimension axis lies in the same or opposite direction to the
            % image axis
            flip_orientation = MimImageCoordinateUtilities.GetFlipFromOrientations(orientation_1, orientation_2, orientation_3, dimension_1, dimension_2, dimension_3);

            % Check the resulting vector is valid
            if (sum(permutation_vector == 1) ~= 1) || (sum(permutation_vector == 2) ~= 1) || (sum(permutation_vector == 3) ~= 1) || ...
                    ~isempty(setdiff(permutation_vector, [1,2,3]))
                reporting.Error('MimImageCoordinateUtilities:InvalidPermutationVector', 'GetDimensionPermutationVectorFromAnatomicalOrientation() resulted in an invalid permutation vector');
            end
        end

        function [dimension_number_1, dimension_number_2, dimension_number_3] = GetDimensionIndicesFromOrientations(orientation_vector_1, orientation_vector_2, reporting)
            % The orientation vector is formed of cosines. Typically these will
            % be 1s and 0s but we allow for small variations in the angles.

            [~, dimension_number_1] = max(abs(orientation_vector_1(:)));
            remaining_dimensions = setdiff([1,2,3], dimension_number_1);
            reduced_orientation_vector_2 = orientation_vector_2(remaining_dimensions);
            [~, dimension_number_2_from_reduced_set] = max(abs(reduced_orientation_vector_2(:)));
            dimension_number_2 = remaining_dimensions(dimension_number_2_from_reduced_set);
            dimension_number_3 = setdiff([1, 2, 3], [dimension_number_1, dimension_number_2]);
        end

        function spline_points = CreateSplineCurve(points, num_points)
            number_of_points = size(points, 1);
            number_of_coordinates = size(points, 2);
            extended_points = zeros(number_of_points + 2, number_of_coordinates);
            extended_points(2 : number_of_points + 1, :) = points;
            extended_points(1, :) = extended_points(2, :) - (extended_points(3, :) - extended_points(2, :));
            extended_points(end, :) = extended_points(end - 1, :) - (extended_points(end - 2, : ) - extended_points(end - 1, :));

            extended_number_of_points = size(extended_points, 1);
            interval_values = linspace(0, 1, num_points + 1);
            interval_values2 = interval_values.^2;
            interval_values3 = interval_values.^3;

            for index = 2 : extended_number_of_points - 2
                coeffs = (1/6).*[...
                        extended_points(index - 1,:)  + 4*extended_points(index, :) + extended_points(index + 1, :); ...
                    - 3*extended_points(index - 1, :) + 3*extended_points(index + 1, :); ...
                      3*extended_points(index - 1, :) - 6*extended_points(index, :) + 3*extended_points(index + 1, :); ...
                    -   extended_points(index - 1, :) + 3*extended_points(index, :) - 3*extended_points(index + 1, :) + extended_points(index+2, :) ...
                ]';

                interval = [ones(size(interval_values)); interval_values; interval_values2; interval_values3];
                spline_points(:, (index - 2)*num_points + 1 : (index - 1)*num_points + 1) = coeffs*interval;
            end
        end

        function dicom_coordinates = PTKToDicomCoordinates(ptk_coordinates, template_image)
            offset = template_image.GetDicomOffset;
            dicom_coordinates = ptk_coordinates + repmat(offset, size(ptk_coordinates, 1), 1);
        end

        function [d_x, d_y, d_z] = PTKToDicomCoordinatesCoordwise(p_x, p_y, p_z, template_image)
            offset = template_image.GetDicomOffset;
            d_x = p_x + offset(1);
            d_y = p_y + offset(2);
            d_z = p_z + offset(3);
        end

        function ptk_coordinates = DicomToPTKCoordinates(dicom_coordinates, template_image)
            offset = template_image.GetDicomOffset;
            ptk_coordinates = dicom_coordinates - repmat(offset, size(dicom_coordinates, 1), 1);
        end

        function dicom_coordinates = PTKToCornerCoordinates(ptk_coordinates, template_image)
            offset = template_image.GetCornerOffset;
            dicom_coordinates = ptk_coordinates + repmat(offset, size(ptk_coordinates, 1), 1);
        end

        function [d_x, d_y, d_z] = PTKToCornerCoordinatesCoordwise(p_x, p_y, p_z, template_image)
            offset = template_image.GetCornerOffset;
            d_x = p_x + offset(1);
            d_y = p_y + offset(2);
            d_z = p_z + offset(3);
        end

        function ptk_coordinates = CornerToPTKCoordinates(corner_coordinates, template_image)
            offset = template_image.GetCornerOffset;
            ptk_coordinates = corner_coordinates - repmat(offset, size(corner_coordinates, 1), 1);
        end

        function [ptk_x, ptk_y, ptk_z] = CoordinatesMmToPTKCoordinates(ic, jc, kc)
            ptk_x = jc;
            ptk_y = ic;
            ptk_z = -kc;
        end

        function coordinates_mm = PTKCoordinatesToCoordinatesMm(ptk_coordinates)
            coordinates_mm = [ptk_coordinates(:, 2), ptk_coordinates(:, 1), - ptk_coordinates(:, 3)];
        end

        function dicom_coordinates = PTKCoordinatesToCornerCoordinates(ptk_coordinates, template_image)
            offset = template_image.GetCornerOffset;
            dicom_coordinates = ptk_coordinates + repmat(offset, size(ptk_coordinates, 1), 1);
        end

        function dicom_coordinates = CoordinatesMmToCornerCoordinates(ptk_coordinates, template_image)
            voxel_size = template_image.VoxelSize;

            % Adjust to coordinates at centre of first voxel
            offset = -voxel_size/2;
            offset = [offset(2), offset(1), -offset(3)];

            % Shift the global origin to the first slice of the image
            global_origin = [0, 0, 0];

            % Adjust to Dicom origin
            offset = offset + global_origin;

            dicom_coordinates = [ptk_coordinates(:, 2), ptk_coordinates(:, 1), - ptk_coordinates(:, 3)];
            dicom_coordinates = dicom_coordinates + repmat(offset, size(ptk_coordinates, 1), 1);
        end

        function ptk_coordinates = CornerToCoordinatesMm(dicom_coordinates, template_image)
            voxel_size = template_image.VoxelSize;

            offset = -voxel_size/2;
            offset = [offset(2), offset(1), -offset(3)];

            dicom_coordinates = dicom_coordinates - repmat(offset, size(dicom_coordinates, 1), 1);

            ptk_coordinates = [dicom_coordinates(:, 2), dicom_coordinates(:, 1), - dicom_coordinates(:, 3)];
        end

        function ptk_coordinates = ConvertToPTKCoordinates(coordinates, coordinate_system, template_image)
            switch coordinate_system
                case MimCoordinateSystem.PTK
                    ptk_coordinates = coordinates;
                case MimCoordinateSystem.Dicom
                    ptk_coordinates = MimImageCoordinateUtilities.DicomToPTKCoordinates(coordinates, template_image);
                case MimCoordinateSystem.DicomUntranslated
                    ptk_coordinates = MimImageCoordinateUtilities.CornerToPTKCoordinates(coordinates, template_image);
                otherwise
                    reporting.Error('MimImageCoordinateUtilities:UnsupportedCoordinateSystem', 'The coordinate system specified by parameter coordinate_system is not supported');
            end
        end

        function coordinates = ConvertFromPTKCoordinates(ptk_coordinates, coordinate_system, template_image)
            switch coordinate_system
                case MimCoordinateSystem.PTK
                    coordinates = ptk_coordinates;
                case MimCoordinateSystem.Dicom
                    coordinates = MimImageCoordinateUtilities.PTKToDicomCoordinates(ptk_coordinates, template_image);
                case MimCoordinateSystem.DicomUntranslated
                    coordinates = MimImageCoordinateUtilities.PTKToCornerCoordinates(ptk_coordinates, template_image);
                otherwise
                    reporting.Error('MimImageCoordinateUtilities:UnsupportedCoordinateSystem', 'The coordinate system specified by parameter coordinate_system is not supported');
            end
        end

        function [c_x, c_y, c_z] = ConvertFromPTKCoordinatesCoordwise(p_x, p_y, p_z, coordinate_system, template_image)
            switch coordinate_system
                case MimCoordinateSystem.PTK
                    c_x = p_x;
                    c_y = p_y;
                    c_z = p_z;
                case MimCoordinateSystem.Dicom
                    [c_x, c_y, c_z] = MimImageCoordinateUtilities.PTKToDicomCoordinatesCoordwise(p_x, p_y, p_z, template_image);
                case MimCoordinateSystem.DicomUntranslated
                    [c_x, c_y, c_z] = MimImageCoordinateUtilities.PTKToCornerCoordinatesCoordwise(p_x, p_y, p_z, template_image);
                otherwise
                    reporting.Error('MimImageCoordinateUtilities:UnsupportedCoordinateSystem', 'The coordinate system specified by parameter coordinate_system is not supported');
            end
        end

        function voxel_indices = AddNearestNeighbours(voxel_indices, template_image)
            if isempty(voxel_indices)
                return;
            end
            [~, linear_offsets27] = MimImageCoordinateUtilities.GetLinearOffsets(template_image.ImageSize);
            voxel_indices = repmat(int32(voxel_indices), 27, 1) + repmat(int32(linear_offsets27'), 1, length(voxel_indices));
            voxel_indices = unique(voxel_indices(:));
        end

        function global_coordinates = GetGlobalCoordinatesForPoints(point_list, template_image)
            xc = [point_list.CoordX];
            yc = [point_list.CoordY];
            zc = [point_list.CoordZ];
            ptk_coords = [xc', yc', zc'];
            coordinates_mm = MimImageCoordinateUtilities.PTKCoordinatesToCoordinatesMm(ptk_coords);
            global_coordinates = round(template_image.CoordinatesMmToGlobalCoordinates(coordinates_mm));
        end

        function global_indices = GetGlobalIndicesForPoints(point_list, template_image)
            global_coordinates = MimImageCoordinateUtilities.GetGlobalCoordinatesForPoints(point_list, template_image);
            global_indices = template_image.GlobalCoordinatesToGlobalIndices(global_coordinates);
        end

        function dist = DistanceBetweenPoints(point_1, point_2)
            dist = norm([point_1.CoordX - point_2.CoordX, point_1.CoordY - point_2.CoordY, point_1.CoordZ - point_2.CoordZ]);
        end

        % Select an appropriate image orientation for exporting an image
        function orientation = ChooseOrientation(voxel_size)
            orientation = GemImageOrientation.XY;
            [sorted_voxel_size, sorted_voxel_size_index] = sort(voxel_size, 'descend');
            if abs(sorted_voxel_size(1) - sorted_voxel_size(2)) > abs(sorted_voxel_size(2) - sorted_voxel_size(3))
                if sorted_voxel_size_index(1) == 1
                    orientation = GemImageOrientation.XZ;
                end
            end
        end

        function offset_voxels = GetOriginOffsetVoxels(from_image, to_image)
            if ~isempty(to_image.GlobalOrigin) && ~isempty(from_image.GlobalOrigin)
                offset_mm = to_image.GlobalOrigin - from_image.GlobalOrigin;
                offset_voxels = offset_mm./from_image.VoxelSize;
            else
                offset_voxels = [0, 0, 0];
            end
        end

        function [Y, X] = GetSliceCoordinates(coordinates_3d, dimension)
            % Returns the 2D [Y, X] coordinates for a 2D slice in the
            % specified orientation given a set of 3D [i,j,k] coordinates
           switch dimension
               case GemImageOrientation.XZ
                   Y = coordinates_3d(:, 3);
                   X = coordinates_3d(:, 2);
               case GemImageOrientation.YZ
                   Y = coordinates_3d(:, 3);
                   X = coordinates_3d(:, 1);
               case GemImageOrientation.XY
                   Y = coordinates_3d(:, 1);
                   X = coordinates_3d(:, 2);
               otherwise
                   error('Unsupported dimension');
           end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function cellInfo = StatesBus(varargin)
% STATESBUS returns a cell array containing bus object information
%
% Optional Input: 'false' will suppress a call to Simulink.Bus.cellToObject
%                 when the MATLAB file is executed.
% The order of bus element attributes is as follows:
%   ElementName, Dimensions, DataType, SampleTime, Complexity, SamplingMode, DimensionsMode, Min, Max, DocUnits, Description

suppressObject = false;
if nargin == 1 && islogical(varargin{1}) && varargin{1} == false
    suppressObject = true;
elseif nargin > 1
    error('Invalid input argument(s) encountered');
end

cellInfo = { ...
  { ...
    'StatesBus', ...
    '', ...
    '', ...
    'Auto', ...
    '-1', {...
{'V_body', [3 1], 'double', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'Omega_body', [3 1], 'double', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'Euler', [3 1], 'double', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'Accel_body', [3 1], 'double', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'dOmega_body', [3 1], 'double', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'V_ned', [3 1], 'double', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'X_ned', [3 1], 'double', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'LLA', [1 3], 'double', -1, 'real', 'Sample', 'Fixed', [], [], '', sprintf('Latitude, Longitude, Altitude')}; ...
{'DCM_be', [3 3], 'double', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
    } ...
  } ...
}';

if ~suppressObject
    % Create bus objects in the MATLAB base workspace
    Simulink.Bus.cellToObject(cellInfo)
end

author: Firmament-Autopilot

repo: FMT-Model