text:
function kern = cmpndKernExpandParamTransformSettings(kern, paramtransformsettings)

% CMPNDKERNEXPANDPARAMTRANSFORMSETTINGS Create kernel structure from CMPND kernel's parameter transformation settings.
%
% FORMAT
%
% DESC returns a compound kernel structure filled with the
% parameter transformation settings in the given cell array. This
% is used as a helper function to enable parameters to be optimised
% in, for example, the NETLAB optimisation functions.
%
% ARG kern : the kernel structure in which the parameters are to be
% placed.
%
% ARG param : cell array of parameter transformation settings which
% are to be placed in the kernel structure. Each setting needs to
% correspond to the transformations used by the kernel. For
% example, some transformation might need knowledge of the desired
% output range, so the transformation setting could be e.g. [0 100],
% whereas some other transformations might need more complicated
% setting information.
%
% RETURN kern : kernel structure with the given parameter
% transformation settings in the relevant locations.
%
% SEEALSO : cmpndKernParamInit, cmpndKernExtractParam, kernExpandParam
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006
%
% COPYRIGHT : Jaakko Peltonen, 2011

% KERN

%fprintf(1,'cmpndKernExpandParamTransformSettings step1\n');

%fprintf(1,'cmpndKernExpandParam, params, kern.paramGroups\n')
%params
%full(kern.paramGroups)

% kern.paramGroups
% paramtransformsettings
% length(paramtransformsettings)
% size(kern.paramGroups,2)

% There should be as many transformation settings provided as there
% are shared parameters in the compound kernel.
if length(paramtransformsettings)~=size(kern.paramGroups,2),
  error(sprintf('Problem in cmpndKernExpandParamTransformSettings: expected %d transformation-settings, received %d\n',size(kern.paramGroups,2),length(paramtransformsettings)));
end;

%fprintf(1,'cmpndKernExpandParamTransformSettings step1b\n');
%pause

% In a compound (cmpnd) kernel, some parameters may be shared
% between groups. Each provided transformation-setting is
% assumed to be specific to a single shared parameter, so
% the transformation-setting has to be duplicated for each group
% using that shared parameter.

expandedsettings=cell(size(kern.paramGroups,1),1);

%fprintf(1,'cmpndKernExpandParamTransformSettings step1b2\n');

for k=1:length(expandedsettings),
  %fprintf(1,'cmpndKernExpandParamTransformSettings step1b3\n');
  % Find the shared parameter that correspond to the k:th non-shared parameter
  k2=find(kern.paramGroups(k,:)==1);

  expandedsettings{k}=paramtransformsettings{k2};
end;

%fprintf(1,'expanding CMPND kernel param-transforms settings:\n');
%paramtransformsettings
%expandedsettings
%pause

%fprintf(1,'cmpndKernExpandParamTransformSettings step2\n');

% Distribute each group of transformation-settings into the
% corresponding component-kernel of the compound kernel.
startVal = 1;
endVal = 0;
for i = 1:length(kern.comp)
  % Each component-kernel is assumed to require as many transformation
  % settings as it has parameters.
  endVal = endVal + kern.comp{i}.nParams;
  %fprintf(1,'cmpndKernExpandParamTransformSettings step3\n');
  kern.comp{i} = kernExpandParamTransformSettings(kern.comp{i}, expandedsettings(startVal:endVal));
  %fprintf(1,'cmpndKernExpandParamTransformSettings step4\n');
  startVal = endVal + 1;
end

%fprintf(1,'cmpndKernExpandParamTransformSettings done\n');

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function f=ftorque(tt)
t1max=34;t2max=12;
t1=tt(1,:);t2=tt(2,:);
f=0;
for i=1:size(tt,2)
    if abs(t1(i)) < t1max
        tc1=0;
    else
        tc1=abs(t1(i))-t1max;
    end
    if abs(t2(i)) < t2max
        tc2=0;
    else
        tc2=abs(t2(i))-t2max;
    end
    f=f+tc1+tc2;
end

author: Sable

repo: mcbench-benchmarks


text:
function qfsum = eiqf ( nt, t, mlt, wts, nwts, ndx, key, f )

%*****************************************************************************80
%
%% EIQF evaluates an interpolatory quadrature formula.
%
%  Discussion:
%
%   The knots, weights and integrand are supplied.
%
%   All knots with nonzero NDX are used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, real T(NT), the knots.
%
%    Input, integer MLT(NT), the multiplicity of the knots.
%
%    Input, real WTS(NWTS), the weights.
%
%    Input, integer NWTS, the number of weights.
%
%    Input, integer NDX(NT), used to index the array WTS.
%    If KEY = 1, then NDX need not be preset.  For more details see the
%    comments in CAWIQ.
%
%    Input, integer KEY, indicates the structure of the WTS
%    array.  It will normally be set to 1.  This will cause the weights to be
%    packed sequentially in array WTS.  For more details see the comments
%    in CAWIQ.
%
%    Input, function F, the name of a routine which
%    evaluates the function and some of its derivatives.  The routine
%    must have the form
%      function value = f ( x, i )
%    and return in VALUE the value of the I-th derivative of the function
%    at X.  The highest value of I will be the maximum value in MLT minus
%    one.  The value X will always be a knot.
%
%    Output, real QFSUM, the value of the quadrature formula
%    applied to F.
%
  l = abs ( key );

  if ( l < 1 || 4 < l )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'EIQF - Fatal error!\n' );
    fprintf ( 1, '  Magnitude of KEY must be between 1 and 4.\n' );
    error ( 'EIQF - Fatal error!' );
  end

  qfsum = 0.0;
  for j = 1 : nt
    l = abs ( ndx(j) );
    if ( l ~= 0 )
      p = 1.0;
      for i = 1 : mlt(j)
        qfsum = qfsum + wts(l+i-1) * f ( t(j), i - 1 ) / p;
        if ( key <= 0 )
          p = p * i;
        end
      end
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = min(varargin)
%MIN   Minimum value of a DISKFUN in one direction.
%   MIN(f) returns a chebfun representing the minimum of the DISKFUN along
%   the latitude direction, i.e, MIN(f) = @(lambda) max( f ( lambda, : ) ).
%
%   MIN(f, [], dim) returns a CHEBFUN representing the minimum of f along
%   direction DIM. DIM = 1 is used to compute max along the radial-direction
%   and DIM = 2 along the angular-direction.
%
%   WARNING: This function is not always accurate to the expected precision.
%
%   For the global minimum use MIN2.
%
% See also DISKFUN/MAX, DISKFUN/MAX2, DISKFUN/MIN2, DISKFUN/MINANDMAX2.

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

[varargout{1:nargout}] = min@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [C] = spm_mc_loss_C(x,P)
% cost function for the mountain car problem
% problem
% FORMAT [C] = spm_mc_loss_C(x,P)
%
% x     - hidden states
% v     - exogenous inputs
% P.x,k - parameters for gradient function:     G(x(1),P.p)
% P.q,p - parameters for cost or loss-function: C(x(1),P.q)
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_mc_loss_C.m 3757 2010-03-08 11:41:53Z guillaume $

% gradient (G) (quadratic potential = (x(1) - P.x)^2*P.k/2)
%--------------------------------------------------------------------------
C   = abs(x(1,:) - 1) > 1/4;

author: spm


text:
classdef prtBrvDpHmm < prtBrvHmm

    methods
        function self = prtBrvDpHmm(varargin)
            self.transitionProbabilities = prtBrvDiscreteStickBreaking;
            self.initialProbabilities = prtBrvDiscreteStickBreaking;

            self = prtUtilAssignStringValuePairs(self, varargin{:});
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function opt = mpoption_info_most(selector)
%MPOPTION_INFO_MOST  Returns MATPOWER option info for MOST.
%
%   DEFAULT_OPTS = MPOPTION_INFO_MOST('D')
%   VALID_OPTS   = MPOPTION_INFO_MOST('V')
%   EXCEPTIONS   = MPOPTION_INFO_MOST('E')
%
%   Returns a structure for MOST options for MATPOWER containing ...
%   (1) default options,
%   (2) valid options, or
%   (3) NESTED_STRUCT_COPY exceptions for setting options
%   ... depending on the value of the input argument.
%
%   This function is used by MPOPTION to set default options, check validity
%   of option names or modify option setting/copying behavior for this
%   subset of optional MATPOWER options.
%
%   See also MPOPTION.

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

if nargin < 1
    selector = 'D';
end
if have_feature('most')
    switch upper(selector)
        case {'D', 'V'}     %% default and valid options
            opt = struct(...
                'most',   struct(...
                    'build_model',                  1, ...          %% was md.CreateQP
                    'solve_model',                  1, ...          %% was md.Solve
                    'resolve_new_cost',             0, ...          %% was md.ReSolveNewCoordCost
                    'dc_model',                     1, ...          %% was md.DCMODEL = []
                    'fixed_res',                    -1, ...         %% was md.IncludeFixedReserves = []
                    'q_coordination',               0, ...          %% was md.QCoordination
                    'security_constraints',         -1, ...         %% was md.SecurityConstrained = []
                    'storage',                      struct(...
                        'terminal_target',              -1, ...     %% was md.Storage.ForceExpectedTerminalStorage
                        'cyclic',                       0), ...     %% was md.Storage.ForceCyclicStorage
                    'uc',                           struct(...
                        'run',                          -1, ...     %% was missing
                        'cyclic',                       0), ...     %% was md.UC.CyclicCommitment
                    'alpha',                        0, ...          %% was md.alpha
                    'solver',                       'DEFAULT', ...  %% was md.QP.opt.alg
                    'skip_prices',                  0, ...          %% was md.QP.opt.skip_prices
                    'price_stage_warn_tol',         1e-7 ...
                ) ...
            );
        case 'E'            %% exceptions used by nested_struct_copy() for applying
            opt = struct([]);   %% no exceptions
        otherwise
            error('mpoption_info_most: ''%s'' is not a valid input argument', selector);
    end
else
    opt = struct([]);       %% MOST is not available
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [p,dp] = spm_LAP_eval(M,qu,qh)
% Evaluate precisions for a LAP model
% FORMAT [p dp] = spm_LAP_eval(M,qu,qh)
%
% p.h     - vector of precisions for causal states (v)
% p.g     - vector of precisions for hidden states (x)
%
% dp.h.dx - dp.h/dx
% dp.h.dv - dp.h/dv
% dp.h.dh - dp.h/dh
%
% dp.g.dx - dp.g/dx
% dp.g.dv - dp.g/dv
% dp.g.dg - dp.g/dg
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_LAP_eval.m 6290 2014-12-20 22:11:50Z karl $

% Get states {qu.v{1},qu.x{1}} in hierarchical form (v{i},x{i})
%--------------------------------------------------------------------------
N          = length(M);
v          = cell(N,1);
x          = cell(N,1);
v(1:N - 1) = spm_unvec(qu.v{1},{M(1 + 1:N).v});
x(1:N - 1) = spm_unvec(qu.x{1},{M(1:N - 1).x});

% precisions
%==========================================================================
for i = 1:N

    % precision of causal and hidden states
    %----------------------------------------------------------------------
    try
        h{i,1} = spm_vec(feval(M(i).ph,x{i},v{i},qh.h{i},M(i)));
    catch
        h{i,1} = sparse(M(i).l,1);
    end
    try
        g{i,1} = spm_vec(feval(M(i).pg,x{i},v{i},qh.g{i},M(i)));
    catch
        g{i,1} = sparse(M(i).n,1);
    end

end

% Concatenate over hierarchical levels
%--------------------------------------------------------------------------
p.h  = spm_cat(h);
p.g  = spm_cat(g);

if nargout < 2, return, end

% gradients
%==========================================================================

% assume precisions can be functions of hyper-parameters and states
%--------------------------------------------------------------------------
try method.h = M(1).E.method.h; catch, method.h = 1; end
try method.g = M(1).E.method.g; catch, method.g = 1; end
try method.x = M(1).E.method.x; catch, method.x = 1; end
try method.v = M(1).E.method.v; catch, method.v = 1; end

% number of variables
%--------------------------------------------------------------------------
nx      = numel(spm_vec(x));
nv      = numel(spm_vec(v));
hn      = numel(spm_vec(qh.h));
gn      = numel(spm_vec(qh.g));
nh      = size(p.h,1);
ng      = size(p.g,1);

dp.h.dh = sparse(nh,hn);
dp.g.dg = sparse(ng,gn);
dp.h.dx = sparse(nh,nx);
dp.h.dv = sparse(nh,nv);
dp.g.dx = sparse(ng,nx);
dp.g.dv = sparse(ng,nv);

% gradients w.r.t. h only (no state-dependent noise)
%----------------------------------------------------------------------
if method.h || method.g

    for i = 1:N

        % precision of causal and hidden states
        %--------------------------------------------------------------
        dhdh{i,i} = spm_diff(M(i).ph,x{i},v{i},qh.h{i},M(i),3);
        dgdg{i,i} = spm_diff(M(i).pg,x{i},v{i},qh.g{i},M(i),3);

    end

    % Concatenate over hierarchical levels
    %------------------------------------------------------------------
    dp.h.dh = spm_cat(dhdh);
    dp.g.dg = spm_cat(dgdg);

end

% gradients w.r.t. causal states
%----------------------------------------------------------------------
if method.v

    for i = 1:N

        % precision of causal states
        %--------------------------------------------------------------
        dhdv{i,i} = spm_diff(M(i).ph,x{i},v{i},qh.h{i},M(i),2);

        % precision of hidden states
        %--------------------------------------------------------------
        dgdv{i,i} = spm_diff(M(i).pg,x{i},v{i},qh.g{i},M(i),2);

    end

    % Concatenate over hierarchical levels
    %------------------------------------------------------------------
    dp.h.dv = spm_cat(dhdv);
    dp.g.dv = spm_cat(dgdv);

end

% gradients w.r.t. hidden states
%----------------------------------------------------------------------
if method.x

    for i = 1:N

        % precision of causal states
        %--------------------------------------------------------------
        dhdx{i,i} = spm_diff(M(i).ph,x{i},v{i},qh.h{i},M(i),1);

        % precision of hidden states
        %--------------------------------------------------------------
        dgdx{i,i} = spm_diff(M(i).pg,x{i},v{i},qh.g{i},M(i),1);

    end

    % Concatenate over hierarchical levels
    %------------------------------------------------------------------
    dp.h.dx = spm_cat(dhdx);
    dp.g.dx = spm_cat(dgdx);

end


text:
function [k,v]=findpeaks(x,m,w)
%FINDPEAKS finds peaks with optional quadratic interpolation [K,V]=(X,M,W)
%
%  Inputs:  X        is the input signal (does not work with UInt datatype)
%           M        is mode:
%                       'q' performs quadratic interpolation
%                       'v' finds valleys instead of peaks
%           W        is the width tolerance; a peak will be eliminated if there is
%                    a higher peak within +-W samples
%
% Outputs:  K        are the peak locations in X (fractional if M='q')
%           V        are the peak amplitudes: if M='q' the amplitudes will be interpolated
%                    whereas if M~='q' then V=X(K).

% Outputs are column vectors regardless of whether X is row or column.
% If there is a plateau rather than a sharp peak, the routine will place the
% peak in the centre of the plateau. When the W input argument is specified,
% the routine will eliminate the lower of any pair of peaks whose separation
% is <=W; if the peaks have exactly the same height, the second one will be eliminated.
% All peak locations satisfy 1<K<length(X).
%
% If no output arguments are specified, the results will be plotted.
%

%	   Copyright (C) Mike Brookes 2005
%      Version: $Id: findpeaks.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<2
    m=' ';
end
nx=length(x);
if any(m=='v')
    x=-x(:);        % invert x if searching for valleys
else
    x=x(:);        % force to be a column vector
end
dx=x(2:end)-x(1:end-1);
r=find(dx>0);
f=find(dx<0);

if length(r)>0 & length(f)>0    % we must have at least one rise and one fall
    dr=r;
    dr(2:end)=r(2:end)-r(1:end-1);
    rc=repmat(1,nx,1);
    rc(r+1)=1-dr;
    rc(1)=0;
    rs=cumsum(rc); % = time since the last rise

    df=f;
    df(2:end)=f(2:end)-f(1:end-1);
    fc=repmat(1,nx,1);
    fc(f+1)=1-df;
    fc(1)=0;
    fs=cumsum(fc); % = time since the last fall

    rp=repmat(-1,nx,1);
    rp([1; r+1])=[dr-1; nx-r(end)-1];
    rq=cumsum(rp);  % = time to the next rise

    fp=repmat(-1,nx,1);
    fp([1; f+1])=[df-1; nx-f(end)-1];
    fq=cumsum(fp); % = time to the next fall

    k=find((rs<fs) & (fq<rq) & (floor((fq-rs)/2)==0));   % the final term centres peaks within a plateau
    v=x(k);

    if any(m=='q')         % do quadratic interpolation
        b=0.5*(x(k+1)-x(k-1));
        a=x(k)-b-x(k-1);
        j=(a>0);            % j=0 on a plateau
        v(j)=x(k(j))+0.25*b(j).^2./a(j);
        k(j)=k(j)+0.5*b(j)./a(j);
        k(~j)=k(~j)+(fq(k(~j))-rs(k(~j)))/2;    % add 0.5 to k if plateau has an even width
    end

    % now purge nearby peaks

    if nargin>2
        j=find(k(2:end)-k(1:end-1)<=w);
        while any(j)
            j=j+(v(j)>=v(j+1));
            k(j)=[];
            v(j)=[];
            j=find(k(2:end)-k(1:end-1)<=w);
        end
    end
else
    k=[];
    v=[];
end
if any(m=='v')
    v=-v;    % invert peaks if searching for valleys
end
if ~nargout
    if any(m=='v')
        x=-x;    % re-invert x if searching for valleys
        ch='v';
    else
        ch='^';
    end
    plot(1:nx,x,'-',k,v,ch);
end

author: zhenghuatan

repo: rVAD

save_path: github-repos/MATLAB/zhenghuatan-rVAD


text:
function [wholepath] = ComputeGeodesicProjection(mesh_total, A, start_vertex, verticeslist)

% --- Step 1: iterative reorganization of the vertices' list
% (each vertex must appear in the right order so that they can be
% subsequently linked together one by one).
% This uses mesh_vertex_nearest (D. Weber) and is based on the assumption
% that verticeslist are relatively regularly situated on a circle.

% NOTE: sometimes this function can return a small loop instead of the full
% perimeter and it's mandatory to check that with e.g. the size (number of
% vertices of the perimeter. Just starting the reorglist
% with another start_vertex usually overcomes that problem.

reorglist=verticeslist(start_vertex);
remaininglist=verticeslist;
remaininglist(start_vertex)=[];

% search the nearest vertex for the first vertex of the list (start_vertex)
[nextindex,nextvalue]=mesh_vertex_nearest(mesh_total.vertices(remaininglist,:),mesh_total.vertices(verticeslist(start_vertex),:));
clear nextvalue;

% create the reorglist (reorganized list)
reorglist=[reorglist remaininglist(nextindex)];

% delete this vertex from the list of vertices (i.e. from the pool of
% vertices to be reorganized)
remaininglist(nextindex)=[];

% repeat that one more time so that we are sufficiently far from the
% start_vertex to add it to the pool of vertices to be reorganized.
% Then continue iteratively until we close the loop by finding the start_vertex.
[nextindex,nextvalue]=mesh_vertex_nearest(mesh_total.vertices(remaininglist,:),mesh_total.vertices(reorglist(2),:));
clear nextvalue;
reorglist=[reorglist remaininglist(nextindex)];
remaininglist(nextindex)=[];

% add the start_vertex
remaininglist= [remaininglist verticeslist(start_vertex)];

% continue to reorganize until start_vertex
for z= 3: size(verticeslist,2)
    [nextindex,nextvalue]=mesh_vertex_nearest(mesh_total.vertices(remaininglist,:),mesh_total.vertices(reorglist(z),:));
    clear nextvalue;
    reorglist=[reorglist remaininglist(nextindex)];
    if remaininglist(nextindex) == verticeslist(start_vertex), break, end
    remaininglist(nextindex)=[];
end

clear remaininglist;
clear nextindex;
clear nextvalue;

% ---- Step 2: find the shortest geodesic path to link each vertex of the
% reorglist. This uses dijk, a function by Michael G. Kay (Matlog toolbox);
% which itself requires the adjacency matrix (A).
% The whole perimeter of the pial ROI is stored in the vector "wholepath".

k=size(reorglist,2)-1;
wholepath=reorglist(1);

for q=1:k
    [D,P]=dijk(A,reorglist(q),reorglist(q+1));
    wholepath=[wholepath P(2:end)];
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = logsumexp(varargin)
%LOGSUMEXP
%
% y = LOGSUMEXP(x)
%
% Computes/declares log of sum of exponentials log(sum(exp(x)))
%
% Implemented as evalutation based nonlinear operator. Hence, the convexity
% of this function is exploited to perform convexity analysis and rigorous
% modelling.

switch class(varargin{1})

    case 'double'
        x = varargin{1};
        varargout{1} = log(sum(exp(x)));

    case 'sdpvar'

        if min(size(varargin{1}))>1
            x = varargin{1};
            y = [];
            for i = 1:size(x,2)
                y = [y yalmip('define',mfilename,x(:,i))];
            end
            varargout{1} = y;
        elseif max(size(varargin{1}))==1
            varargout{1} = varargin{1};
        else
            varargout{1} = yalmip('define',mfilename,varargin{1});
        end

    case 'char'

        X = varargin{3};

        operator = struct('convexity','convex','monotonicity','none','definiteness','none','model','callback');
        operator.bounds = @bounds;
        operator.derivative = @(x) exp(x)./(sum(exp(x)));
        operator.convexhull = @convexhull;

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = X;

    otherwise
        error('SDPVAR/LOG called with CHAR argument?');
end

function [L, U] = bounds(xL,xU)

L = log(sum(exp(xL)));
U = log(sum(exp(xU)));

function [Ax, Ay, b] = convexhull(xL,xU)
Ax = [];
Ay = [];
b = [];

author: shamilmamedov


text:
% An example of nonnegative CP decomposition from partial observations
%% Generate synthetic 4-order tensor
rand('seed',0); randn('seed',0);
N1 = 30; N2 = 30; N3 = 30; N4 = 30; % tensor dimension
R = 10; % tensor rank

% randomly generate factor matrices
A1 = max(0,randn(N1,R));
A2 = max(0,randn(N2,R));
A3 = max(0,randn(N3,R));
A4 = max(0,randn(N4,R));
% generate tensor M using the above factor matrices
M = ktensor({A1,A2,A3,A4});
M = full((arrange(M)));

sr = 0.4; % percentage of samples
% randomly choose samples
known = randsample(N1*N2*N3*N4,round(sr*N1*N2*N3*N4));
data = M.data(known);

sn = 60; % signal to noise
% -- add noise --
Noise = max(0,randn(size(data)));
data = data + 10^(-sn/20)*norm(data)/norm(Noise)*Noise;
ndim = [N1,N2,N3,N4];

%% Solve problem
% exact rank works but a rough rank overestimate is also fine
esr = round(1.25*R); % overestimate of rank
opts.tol = 1e-4; opts.maxit = 1000;
t0 = tic;
[A,Out] = ncpc(data,known,ndim,esr,opts);
time = toc(t0);

%% Reporting
relerr = norm(full(A)-M)/norm(M);
fprintf('time = %4.2e, ',time);
fprintf('solution relative error = %4.2e\n\n',relerr);

figure;
semilogy(1:Out.iter, Out.hist_obj,'k-','linewidth',2);
xlabel('iteration','fontsize',12);
ylabel('objective value','fontsize',12)

figure;
semilogy(1:Out.iter, Out.hist_rel(2,:),'k-','linewidth',2);
xlabel('iteration','fontsize',12);
ylabel('relative residual','fontsize',12)

%% <http://www.caam.rice.edu/~optimization/bcu/ncpc/example_ncpc.m Download this m-file>

author: andrewssobral


text:
function cbar = cbarLoadPreset(name, viewer, o);
% Load a preset colorbar, and set it as the colorbar for a mrViewer
% overlay.
%
%  cbar = cbarLoadPreset([name=prompt], [viewer=gcf], [o=1]);
%
% name is the name (filename) of a saved colorbar. viewer is a handle to
% a mrViewer UI. o is the overlay index for the cbar.
%
% ras, 02/24/07.
if notDefined('name')
	name = inputdlg({'Colorbar Name?'}, 'Load a Preset Cbar', 1, {''});
	name = name{1};
end

if notDefined('viewer'),	viewer = gcf;		end
if notDefined('o'),			o = 1;				end

presetsDir = fileparts(which(mfilename));

cbarPath = fullfile(presetsDir, [name '.mat']);
if ~exist(cbarPath, 'file')
	what(presetsDir)
	error(sprintf('Preset %s not found.', name))
end

cbar = cbarLoad(cbarPath);

mrViewSetOverlay(viewer, 'cbar', cbar, o);

return


text:
function c = bitshift(a,b)
	checkinputs(a,b);

	if numel(b)==1
		if b<0
			c = uint64rightshift(uint64(a),uint64(-b));
		else
			c = uint64leftshift(uint64(a),uint64(b));
		end
	else
		c = zeros(size(b),'uint64');
		neg = b < 0;
		if numel(a)==1
			c(neg)  = uint64rightshift(uint64(a),uint64(-b(neg)));
			c(~neg) = uint64leftshift(uint64(a),uint64(b(~neg)));
		else
			c(neg)  = uint64rightshift(uint64(a(neg)),uint64(-b(neg)));
			c(~neg) = uint64leftshift(uint64(a(~neg)),uint64(b(~neg)));
		end
	end
end


text:
function [fragments, fragment_offset] = sequence_fragment(sequence, fragment_length)
% sequence_fragment Returns an array of subsequences
%
% This function returns a cell array of sequence objects that are fragments of the source sequence.
%
% Input:
% - sequence (structure): A valid sequence structure.
% - fragment_length (integer): The length of a single fragment.
%
% Output:
% - fragments (cell): A cell array of sequence descriptors.
% - fragment_offset(vector): A vector of offsets in the original sequence.

fragment_offset = [1:fragment_length:sequence.length, sequence.length+1];

fragments = cell(numel(fragment_offset)-1, 1);

for f = 1:numel(fragment_offset)-1
    indices = fragment_offset(f):fragment_offset(f+1)-1;
    fragments{f} = sequence;

    fragments{f}.groundtruth = sequence.groundtruth(indices, :);
    fragments{f}.indices = sequence.indices(indices);
    fragments{f}.tags.names = sequence.tags.names;
    fragments{f}.tags.data = sequence.tags.data(indices, :);
    fragments{f}.values.names = sequence.values.names;
    fragments{f}.values.data = sequence.values.data(indices, :);
    fragments{f}.images = sequence.images(indices);
    fragments{f}.length = length(indices);

end;

fragment_offset = fragment_offset(1:end-1);

author: votchallenge

repo: toolkit-legacy


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nickabattista@gmail.com
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs, non-invariant beams*)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "Hill+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (nickabattista@gmail.com) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: updates the spring attributes!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function springs_info = update_Springs(dt,current_time,xLag,yLag,springs_info)

%springs_info: col 1: starting spring pt (by lag. discretization)
%              col 2: ending spring pt. (by lag. discretization)
%              col 3: spring stiffness
%              col 4: spring resting lengths

%RL = springs_info(:,4); % resting-length vector

% Contraction Frequency
freq = 1.0;

% CHANGE RESTING LENGTH BTWN SIDES OF JELLYFISH BELL
springs_info(19:end,4) = abs( cos(freq*current_time*pi) );

%NOTE: not 2*pi*ft b/c of abs()

author: nickabattista


text:
function view = measureCortMag(view, roiList, paramsFile)
%
% USAGE: view = measureCortMag(view, [roiList], [paramsFile])
%
%   view: can be either flat or gray.
%
%   roiList: a list of ROI indices specifying which ROIs to
%            use for the cortMag. Defaults to all ROIs.
%
%   paramsFile: a .mat file into which the parms will be saved.
%               If the file exists, it will be read in and those
%               params used as a starting point.
%
% PURPOSE:
%   Compute cortical magnification measurements.
%
% NOTE:
%   This code assumes that the ROI coords have been sorted such that the fovea
%   is at the beginning of the list! If you generate your ROIs with mrFindAreas,
%   you will be fine (it sorts the coords for you). If you draw them with the
%   'create line ROI' tool, you should be OK as long as you start all your lines
%   in the fovea. Otherwise, you are on your own.
%
% FIX THESE ISSUES:
%   * The inter-bin distance should be computed via paths only through layer 1.
%   * The inter-bin distance is currently computed from the bin edge, not the bin center.
%     This won't affect the distance measures, but it does bias the mean phase estimate
%     for each bin. It is tricky to adjust for this- the bins are not regularly spaced.
%
% HISTORY:
%   2002.01.22 RFD (bob@white.stanford.edu): wrote it, porting some
%   code from mrLoadRet-2.5. However, I've taken a fresh start and attempted
%   to streamline things and make them more modular.
%
% 7/16/02 djh, replaced mrSESSION.vAnatomyPath with global vANATOMYPATH

% We need vANATOMYPATH and mrSESSION.subject
global vANATOMYPATH;
global mrSESSION;
global dataTYPES;

if(~exist('roiList','var') | isempty(roiList))
    % default to all
    roiList = [1:length(view.ROIs)];
end

% We currently only support Flat view. All analyses are performed on
% flat view ROIs that _should_ be created by an atlas fit.
%
switch(view.viewType)
case {'Inplane','Volume','Gray'}
    error([mfilename,' doesn''t work for ',view.viewType,'.']);
case 'Flat'
    % Get a gray structure because we need the gray nodes.
    grayView = getSelectedGray;
    if isempty(grayView)
        grayView = initHiddenGray;
    end
    cortMag.ROIs = view.ROIs(roiList);
otherwise
    error([view.viewType,' is unknown!']);
end

if(~exist('paramsFile','var') | ~exist(paramsFile,'file'))
    cortMag.paramsFileName = [];
    nowNum = now;
    cortMag.subject = mrSESSION.subject;
    cortMag.timestamp = datestr(nowNum);
    cortMag.mmPerPix = readVolAnatHeader(vANATOMYPATH);
    cortMag.coThresh = getCothresh(view);
    cortMag.expNumber = getCurScan(view);
    cortMag.nodeIndices = {};
    cortMag.bins = {};
    cortMag.dataType = getDataTypeName(view);
    cortMag.subdir = view.subdir;
    % This is a crude method for determining the hemisphere- but it usually works.
    if(cortMag.ROIs(1).coords(3,1)==1)
        cortMag.hemisphere = 'left';
        cortMag.slice = 1;
    else
        cortMag.hemisphere = 'right';
        cortMag.slice = 2;
    end

else
    % this file _should_ contain a cortmag struct.
    % Note that this file, if it exists, will form the starting point.
    % We will overwrite many of it's fields below.
    tmp = load(paramsFile);
    if(isfield(tmp,'cortMag'))
        cortMag = tmp.cortMag;
    else
        warning(['Params file "',paramsFile,'" exists, but it does '...
                'not contain a cortMag struct! It will be ignored.']);
    end
    % Just in case the paramsFileName is not correct:
    cortMag.paramsFileName = paramsFile;
    % If we were really cool, we'd check the validity of the loaded cortmag struct.
end

% Make sure that the correct data is loaded in the grayView
grayView.curDataType = existDataType(cortMag.dataType);
grayView = loadCorAnal(grayView);

if(isempty(cortMag.paramsFileName))
    cortMag.paramsFileName = ['CortMag_',cortMag.subdir,'_',cortMag.dataType,'_',cortMag.hemisphere,'_',...
            datestr(nowNum,11),datestr(nowNum,5),datestr(nowNum,7),'.mat'];
end

% At this point, we have a gray struct and we have a set of ROIs in flat coordinates.
% Note that the ROI coords are assumed to be sorted already (with the fovea at the
% beginning of the list).

% Now, we should call a function that computes cortmag. But, I'd like to keep all UI code
% here, so we'll do it like this:

satisfied = 0;
while(~satisfied)
    resp = getCortMagParams(cortMag);
    if(~isstruct(resp) & resp==0)
        % a 0 means the user cancelled- they are aborting or are satisfied
        satisfied = 1;
    else
        % if it's not 0, then it's an updated cortMag struct
        oldCortMag = cortMag;
        cortMag = resp;

        % Finding the nodes for each coordinate in all the ROIs is slow. Once we do it, we
        % needn't do it again unless we get a new ROI added to the list. That is, none of the
        % cortMag parameters will affect the nodes that are found. So, we don't need to redo it
        % if a parameter changes.
        if(isempty(cortMag.nodeIndices))
            cortMag = buildCortMagNodes(cortMag, view, grayView);
            cortMag.bins = [];
        end
        % Finding the bins is somewhat slow. We only really need to redo it if the binSize changes.
        if( isempty(cortMag.bins) | cortMag.binDist ~= oldCortMag.binDist ...
            | cortMag.expNumber ~= oldCortMag.expNumber | cortMag.flatDataFlag ~= oldCortMag.flatDataFlag)
            cortMag = buildCortMagBins(cortMag, view, grayView);
            if(strmatch('Atlases', dataTYPES(view.curDataType).name))
                phaseScale = dataTYPES(view.curDataType).atlasParams(cortMag.expNumber).phaseScale(cortMag.slice);
                phaseShift = dataTYPES(view.curDataType).atlasParams(cortMag.expNumber).phaseShift(cortMag.slice);
                fprintf('Fixing atlas phases with shift (%.2f) and scale (%.2f)...', phaseShift, phaseScale);
                for(ii=1:length(cortMag.data))
                    cortMag.data{ii}.ph = mod((cortMag.data{ii}.ph - phaseShift) / phaseScale, 2*pi);
                end
            end
        end
        cortMag = computeCortMagFunction(cortMag, view, grayView);
        plotCortmagResults(cortMag);
    end
end

view.cortMag = cortMag;
if(~isempty(cortMag.paramsFileName))
    %if(~exist(cortMag.paramsFileName,'file') | strcmp(questdlg('Overwrite existing cortMag file?', 'Overwrite?', 'Yes', 'No', 'Yes'),'Yes'))
    if(exist(cortMag.paramsFileName,'file'))
        disp([cortMag.paramsFileName,' exists- launching gui to save it somewhere else.']);
        uisave({'cortMag'}, cortMag.paramsFileName);
    else
        save(cortMag.paramsFileName, 'cortMag');
        disp([cortMag.paramsFileName,' saved.']);
    end
end
return;


text:
function [boxes, blobIndIm, blobIndBoxes, hierarchy] = selective_search_boxes(im, fast_mode, im_width)
% Girshick Wrapper
%
% Based on the demo.m file included in the Selective Search
% IJCV code.
%
% Requires selective search code.

if ~exist('fast_mode', 'var') || isempty(fast_mode)
  fast_mode = true;
end

originalImSize = size(im);

if ~exist('im_width', 'var') || isempty(im_width)
  im_width = [];
  scale = 1;
else
  scale = size(im, 2) / im_width;
end

if scale ~= 1
  im = imresize(im, [NaN im_width]);
end

% Parameters. Note that this controls the number of hierarchical
% segmentations which are combined.
colorTypes = {'Hsv', 'Lab', 'RGI', 'H', 'Intensity'};

% Here you specify which similarity functions to use in merging
simFunctionHandles = {@SSSimColourTextureSizeFillOrig, ...
                      @SSSimTextureSizeFill, ...
                      @SSSimBoxFillOrig, ...
                      @SSSimSize};

% Thresholds for the Felzenszwalb and Huttenlocher segmentation algorithm.
% Note that by default, we set minSize = k, and sigma = 0.8.
% controls size of segments of initial segmentation.
ks = [50 100 150 300];
sigma = 0.8;

% After segmentation, filter out boxes which have a width/height smaller
% than minBoxWidth (default = 20 pixels).
minBoxWidth = 20;

% Comment the following three lines for the 'quality' version
if fast_mode
  colorTypes = colorTypes(1:2); % 'Fast' uses HSV and Lab
  simFunctionHandles = simFunctionHandles(1:2); % Two different merging strategies
  ks = ks(1:2);
end

idx = 1;
for j = 1:length(ks)
  k = ks(j); % Segmentation threshold k
  minSize = k; % We set minSize = k
  for n = 1:length(colorTypes)
    colorType = colorTypes{n};
    [boxesT{idx} blobIndIm{idx} blobIndBoxes{idx} hierarchy{idx} priorityT{idx}] = ...
      Image2HierarchicalGrouping(im, sigma, k, minSize, colorType, simFunctionHandles);
    idx = idx + 1;
  end
end
boxes = cat(1, boxesT{:}); % Concatenate boxes from all hierarchies
priority = cat(1, priorityT{:}); % Concatenate priorities

% Do pseudo random sorting as in paper
priority = priority .* rand(size(priority));
[priority sortIds] = sort(priority, 'ascend');
boxes = boxes(sortIds,:);

boxes = BoxRemoveDuplicates(boxes);

if scale ~= 1
  boxes = floor((boxes - 1) * scale + 1);
  for iii = 1:length(blobIndBoxes)
    blobIndBoxes{iii} = floor((blobIndBoxes{iii} - 1) * scale + 1);
    blobIndIm{iii} = imresize(blobIndIm{iii}, [originalImSize(1) originalImSize(2)], 'nearest');
  end
end

% Filter width of boxes
[nr, nc] = BoxSize(boxes);
idsGood = (nr >= minBoxWidth) & (nc >= minBoxWidth);
boxes = boxes(idsGood,:);

author: nightrome

repo: matconvnet-calvin


text:
classdef MimViewerPanelToolbar < GemPanel
    % MimViewerPanelToolbar. Part of the gui for the TD MIM Toolkit.
    %
    %     This class is used internally within the TD MIM Toolkit to help
    %     build the user interface.
    %
    %     MimViewerPanelToolbar is used to build the controls for the MimViewingPanel
    %
    %
    %     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.
    %

    properties (Access = private)
        ViewerPanel

        Tools

        % Handler to listener for when the mouse status changes
        MouseCursorStatusListener

        % User interface controls
        OrientationPanel
        MouseControlPanel
        ImageOverlayPanel
        WindowLevelPanel
        OrientationButtons
        OpacitySlider
        ImageCheckbox
        OverlayCheckbox
        StatusText
        WindowText
        WindowEditbox
        WindowSlider
        LevelText
        LevelEditbox
        LevelSlider
        MouseControlButtons
    end

    methods
        function obj = MimViewerPanelToolbar(viewer_panel)
            obj = obj@GemPanel(viewer_panel);
            obj.ViewerPanel = viewer_panel;
            obj.BackgroundColour = 'black';
            obj.Tools = viewer_panel.GetToolList;

            obj.MouseCursorStatusListener = addlistener(obj.ViewerPanel, 'MouseCursorStatusChanged', @obj.MouseCursorStatusChangedCallback);

            % Add listeners to respond to changes in the viewer panel
            obj.AddPostSetListener(obj.ViewerPanel, 'WindowLimits', @obj.WindowLimitsChangedCallback);
            obj.AddPostSetListener(obj.ViewerPanel, 'LevelLimits', @obj.LevelLimitsChangedCallback);
            obj.AddPostSetListener(obj.ViewerPanel, 'SelectedControl', @obj.SelectedControlChangedCallback);
            obj.AddPostSetListener(obj.ViewerPanel.GetImageSliceParameters, 'Orientation', @obj.GuiChangeCallback);
            obj.AddPostSetListener(obj.ViewerPanel.GetImageSliceParameters, 'SliceNumber', @obj.GuiChangeCallback);
            obj.AddPostSetListener(obj.ViewerPanel.GetBackgroundImageDisplayParameters, 'Level', @obj.GuiChangeCallback);
            obj.AddPostSetListener(obj.ViewerPanel.GetBackgroundImageDisplayParameters, 'Window', @obj.GuiChangeCallback);
            obj.AddPostSetListener(obj.ViewerPanel.GetBackgroundImageDisplayParameters, 'ShowImage', @obj.GuiChangeCallback);
            obj.AddPostSetListener(obj.ViewerPanel.GetOverlayImageDisplayParameters, 'ShowImage', @obj.GuiChangeCallback);
            obj.AddPostSetListener(obj.ViewerPanel.GetOverlayImageDisplayParameters, 'Opacity', @obj.GuiChangeCallback);
            obj.AddPostSetListener(obj.ViewerPanel.GetOverlayImageDisplayParameters, 'BlackIsTransparent', @obj.GuiChangeCallback);
            obj.AddPostSetListener(obj.ViewerPanel.GetOverlayImageDisplayParameters, 'OpaqueColour', @obj.GuiChangeCallback);
            obj.AddPostSetListener(obj.ViewerPanel.GetMarkerImageDisplayParameters, 'ShowLabels', @obj.GuiChangeCallback);

            % Listen for new image events
            obj.AddEventListener(obj.ViewerPanel.GetBackgroundImageSource, 'NewImage', @obj.NewBackgroundImageCallback);
            obj.AddEventListener(obj.ViewerPanel.GetOverlayImageSource, 'NewImage', @obj.NewOverlayImageCallback);
            obj.AddEventListener(obj.ViewerPanel.GetQuiverImageSource, 'NewImage', @obj.NewQuiverImageCallback);
        end

        function CreateGuiComponent(obj, position)
            CreateGuiComponent@GemPanel(obj, position);

            parent = obj.GetParentFigure;
            keypress_function = @parent.CustomKeyPressedFunction;

            font_size = 9;

            % Buttons for coronal/sagittal/axial views
            obj.OrientationPanel = uibuttongroup('Parent', obj.Parent.GetContainerHandle, 'Units', 'pixels', 'BorderType', 'none', 'SelectionChangeFcn', @obj.OrientationCallback, 'BackgroundColor', 'black', 'ForegroundColor', 'white');
            orientation_buttons = [0 0 0];
            orientation_buttons(1) = uicontrol('Style', 'togglebutton', 'Units', 'pixels', 'Parent', obj.OrientationPanel, 'String', 'Cor', 'Units', 'pixels', 'FontSize', font_size, 'Tag', 'Coronal', 'TooltipString', 'View coronal slices (Y-Z)', 'KeyPressFcn', keypress_function);
            orientation_buttons(2) = uicontrol('Style', 'togglebutton', 'Units', 'pixels', 'Parent', obj.OrientationPanel, 'String', 'Sag', 'Units', 'pixels', 'FontSize', font_size, 'Tag', 'Sagittal', 'TooltipString', 'View sagittal slices (X-Z)', 'KeyPressFcn', keypress_function);
            orientation_buttons(3) = uicontrol('Style', 'togglebutton', 'Units', 'pixels', 'Parent', obj.OrientationPanel, 'String', 'Ax', 'Units', 'pixels', 'FontSize', font_size, 'Tag', 'Axial', 'TooltipString', 'View transverse slices (X-Y)', 'KeyPressFcn', keypress_function);
            obj.OrientationButtons = orientation_buttons;
            set(obj.OrientationButtons(obj.ViewerPanel.Orientation), 'Value', 1);

            % Buttons for each mouse tool
            obj.MouseControlPanel = uibuttongroup('Parent', obj.Parent.GetContainerHandle, 'Units', 'pixels', 'BorderType', 'none', 'SelectionChangeFcn', @obj.ControlsCallback, 'BackgroundColor', 'black', 'ForegroundColor', 'white');
            obj.MouseControlButtons = containers.Map;
            for tool_set = obj.Tools.Tools.values
                tool = tool_set{1};
                tag = tool.Tag;
                obj.MouseControlButtons(tag) = uicontrol('Style', 'togglebutton', 'Units', 'pixels', 'Parent', obj.MouseControlPanel, 'String', tool.ButtonText, 'Units', 'pixels', 'FontSize', font_size, 'Tag', tool.Tag, 'TooltipString', tool.ToolTip, 'KeyPressFcn', keypress_function);
            end
            set(obj.MouseControlButtons(obj.ViewerPanel.SelectedControl), 'Value', 1);

            obj.WindowLevelPanel = uipanel('Parent', obj.Parent.GetContainerHandle, 'Units', 'pixels', 'BorderType', 'none', 'BackgroundColor', 'black', 'ForegroundColor', 'white');
            obj.ImageOverlayPanel = uipanel('Parent', obj.Parent.GetContainerHandle, 'Units', 'pixels', 'BorderType', 'none', 'BackgroundColor', 'black', 'ForegroundColor', 'white');

            obj.StatusText = uicontrol('Style', 'text', 'Parent', obj.ImageOverlayPanel, 'Units', 'pixels', 'FontSize', font_size, 'BackgroundColor', 'black', 'ForegroundColor', 'white');

            obj.OpacitySlider = uicontrol('Style', 'slider', 'Parent', obj.ImageOverlayPanel, 'Units', 'pixels', 'Callback', @obj.OpacitySliderCallback, 'TooltipString', 'Change opacity of overlay', 'Min', 0, 'Max', 100, 'KeyPressFcn', keypress_function, 'Value', obj.ViewerPanel.OverlayOpacity);
            obj.ImageCheckbox = uicontrol('Style', 'checkbox', 'Parent', obj.ImageOverlayPanel, 'Units', 'pixels', 'FontSize', font_size, 'Callback', @obj.ImageCheckboxCallback, 'String', 'Image', 'BackgroundColor', 'black', 'ForegroundColor', 'white', 'TooltipString', 'Show image', 'KeyPressFcn', keypress_function, 'Value', obj.ViewerPanel.ShowImage);
            obj.OverlayCheckbox = uicontrol('Style', 'checkbox', 'Parent', obj.ImageOverlayPanel, 'Units', 'pixels', 'FontSize', font_size, 'Callback', @obj.OverlayCheckboxCallback, 'String', 'Overlay', 'BackgroundColor', 'black', 'ForegroundColor', 'white', 'TooltipString', 'Show overlay over image', 'KeyPressFcn', keypress_function, 'Value', obj.ViewerPanel.ShowOverlay);

            obj.WindowText = uicontrol('Style', 'text', 'Parent', obj.WindowLevelPanel, 'Units', 'pixels', 'FontSize', font_size, 'String', 'Window:', 'BackgroundColor', 'black', 'ForegroundColor', 'white', 'HorizontalAlignment', 'right');
            obj.LevelText = uicontrol('Style', 'text', 'Parent', obj.WindowLevelPanel, 'Units', 'pixels', 'FontSize', font_size, 'String', 'Level:', 'BackgroundColor', 'black', 'ForegroundColor', 'white', 'HorizontalAlignment', 'right');
            obj.WindowEditbox = uicontrol('Style', 'edit', 'Parent', obj.WindowLevelPanel, 'Units', 'pixels', 'FontSize', font_size, 'Callback', @obj.WindowTextCallback, 'BackgroundColor', 'black', 'ForegroundColor', 'white', 'TooltipString', 'Change window (contrast)', 'String', num2str(obj.ViewerPanel.Window));
            obj.LevelEditbox = uicontrol('Style', 'edit', 'Parent', obj.WindowLevelPanel, 'Units', 'pixels', 'FontSize', font_size, 'Callback', @obj.LevelTextCallback, 'BackgroundColor', 'black', 'ForegroundColor', 'white', 'TooltipString', 'Change level (brightness)', 'String', num2str(obj.ViewerPanel.Level));
            window_slider_min = min(obj.ViewerPanel.Window, 0);
            window_slider_max = max(obj.ViewerPanel.Window, 1);
            obj.WindowSlider = uicontrol('Style', 'slider', 'Units', 'pixels', 'Min', window_slider_min, 'Max', window_slider_max, 'Value', 1, 'Parent', obj.WindowLevelPanel, 'Callback', @obj.WindowSliderCallback, 'TooltipString', 'Change window (contrast)', 'KeyPressFcn', keypress_function, 'Value', obj.ViewerPanel.Window);
            level_slider_min = min(obj.ViewerPanel.Level, 0);
            level_slider_max = max(obj.ViewerPanel.Level, 1);
            obj.LevelSlider = uicontrol('Style', 'slider', 'Units', 'pixels', 'Min', level_slider_min, 'Max', level_slider_max, 'Value', 0, 'Parent', obj.WindowLevelPanel, 'Callback', @obj.LevelSliderCallback, 'TooltipString', 'Change level (brightness)', 'KeyPressFcn', keypress_function, 'Value', obj.ViewerPanel.Level);

            % Add context menu
            obj.Tools.UpdateTools;

            % Add custom listeners to allow continuous callbacks from the
            % sliders
            obj.AddEventListener(obj.OpacitySlider, 'ContinuousValueChange', @obj.OpacitySliderCallback);
            obj.AddEventListener(obj.WindowSlider, 'ContinuousValueChange', @obj.WindowSliderCallback);
            obj.AddEventListener(obj.LevelSlider, 'ContinuousValueChange', @obj.LevelSliderCallback);

            obj.ResizePanel(position);
        end

        function delete(obj)
            delete(obj.MouseCursorStatusListener);
        end

        function Resize(obj, position)
            Resize@GemPanel(obj, position);

            if obj.ComponentHasBeenCreated
                obj.ResizePanel(position);
            end
        end

        function UpdateWindowLimits(obj)
            % Sets the minimum and maximum values for the window slider

            window_limits = obj.ViewerPanel.WindowLimits;
            if ~isempty(window_limits)
                set(obj.WindowSlider, 'Min', window_limits(1), 'Max', window_limits(2));
            end

        end

        function UpdateLevelLimits(obj)
            % Sets the minimum and maximum values for the level slider

            level_limits = obj.ViewerPanel.LevelLimits;
            if ~isempty(level_limits)
                set(obj.LevelSlider, 'Min', level_limits(1), 'Max', level_limits(2));
            end
        end

        function SetControl(obj, tag_value)
            % Enable the button for this tool
            set(obj.MouseControlButtons(tag_value), 'Value', 1);
        end

        function UpdateGui(obj)
            main_image = obj.ViewerPanel.BackgroundImage;
            if obj.ComponentHasBeenCreated
                set(obj.OverlayCheckbox, 'Value', obj.ViewerPanel.ShowOverlay);
                set(obj.ImageCheckbox, 'Value', obj.ViewerPanel.ShowImage);
                set(obj.OrientationButtons(obj.ViewerPanel.Orientation), 'Value', 1);
                set(obj.MouseControlButtons(obj.ViewerPanel.SelectedControl), 'Value', 1);

                if ~isempty(main_image) && main_image.ImageExists
                    set(obj.WindowEditbox, 'String', num2str(obj.ViewerPanel.Window));
                    set(obj.WindowSlider, 'Value', obj.ViewerPanel.Window);
                    set(obj.LevelEditbox, 'String', num2str(obj.ViewerPanel.Level));
                    set(obj.LevelSlider, 'Value', obj.ViewerPanel.Level);
                    set(obj.OpacitySlider, 'Value', obj.ViewerPanel.OverlayOpacity);
                end
            end
        end
    end

    methods (Access = private)

        function MouseCursorStatusChangedCallback(obj, ~, ~)
            mouse_cursor_status = obj.ViewerPanel.MouseCursorStatus;
            orientation = obj.ViewerPanel.Orientation;

            if ~mouse_cursor_status.ImageExists
                status_text = 'No image';
            else
                rescale_text = '';

                i_text = int2str(mouse_cursor_status.GlobalCoordX);
                j_text = int2str(mouse_cursor_status.GlobalCoordY);
                k_text = int2str(mouse_cursor_status.GlobalCoordZ);

                if ~isempty(mouse_cursor_status.ImageValue)
                    voxel_value = mouse_cursor_status.ImageValue;
                    if isinteger(voxel_value)
                        value_text = int2str(voxel_value);
                    else
                        value_text = num2str(voxel_value, 3);
                    end

                    rescaled_value = mouse_cursor_status.RescaledValue;
                    rescale_units = mouse_cursor_status.RescaleUnits;
                    if ~isempty(rescale_units) && ~isempty(rescaled_value)
                        rescale_text = [rescale_units ':' int2str(rescaled_value)];
                    end

                    if isempty(mouse_cursor_status.OverlayValue)
                        overlay_text = [];
                    else
                        overlay_value = mouse_cursor_status.OverlayValue;
                        if isinteger(overlay_value)
                            overlay_value_text = int2str(overlay_value);
                        else
                            overlay_value_text = num2str(overlay_value, 3);
                        end
                        overlay_text = [' O:' overlay_value_text];
                    end
                else
                    overlay_text = '';
                    value_text = '-';
                    switch orientation
                        case GemImageOrientation.XZ
                            j_text = '--';
                            k_text = '--';
                        case GemImageOrientation.YZ
                            i_text = '--';
                            k_text = '--';
                        case GemImageOrientation.XY
                            i_text = '--';
                            j_text = '--';
                    end

                end

                status_text = ['X:' j_text ' Y:' i_text ' Z:' k_text ' I:' value_text ' ' rescale_text overlay_text];
            end

            obj.SetStatus(status_text);
        end

        function ResizePanel(obj, position)
            panel_width_pixels = position(3);
            panel_height_pixels = position(4);

            button_width = 32;

            number_of_enabled_control_buttons = 0;

            for tool_tag = obj.MouseControlButtons.keys
                button = obj.MouseControlButtons(tool_tag{1});
                tool = obj.Tools.GetTool(tool_tag{1});
                if tool.IsEnabled(obj.ViewerPanel.Mode, obj.ViewerPanel.SubMode)
                    number_of_enabled_control_buttons = number_of_enabled_control_buttons + 1;
                    set(button, 'Units', 'Pixels', 'Position', [1+button_width*(number_of_enabled_control_buttons - 1), 1, button_width, button_width], 'Visible', 'on');
                else
                    set(button, 'Visible', 'off');
                end
            end

            orientation_width = numel(obj.OrientationButtons)*button_width;
            mouse_controls_width = (number_of_enabled_control_buttons)*button_width;
            mouse_controls_position = panel_width_pixels - mouse_controls_width + 1;
            central_panels_width = max(1, (panel_width_pixels - mouse_controls_width - orientation_width)/2);
            windowlevel_panel_position = orientation_width + central_panels_width + 1;
            set(obj.MouseControlPanel, 'Units', 'Pixels', 'Position', [mouse_controls_position 1 mouse_controls_width panel_height_pixels]);
            set(obj.ImageOverlayPanel, 'Units', 'Pixels', 'Position', [orientation_width+1 1 central_panels_width panel_height_pixels]);
            set(obj.WindowLevelPanel, 'Units', 'Pixels', 'Position', [windowlevel_panel_position 1 central_panels_width panel_height_pixels]);

            % Setup orientation panel
            set(obj.OrientationPanel, 'Units', 'Pixels', 'Position', [1 1 orientation_width panel_height_pixels]);
            set(obj.OrientationButtons(1), 'Units', 'Pixels', 'Position', [1 1 button_width  button_width]);
            set(obj.OrientationButtons(2), 'Units', 'Pixels', 'Position', [1+button_width 1 button_width button_width]);
            set(obj.OrientationButtons(3), 'Units', 'Pixels', 'Position', [1+button_width*2 1 button_width button_width]);

            % Setup image/overlay panel
            halfpanel_height = 16;
            checkbox_width = 70;
            status_width = max(1, central_panels_width - checkbox_width);
            set(obj.ImageCheckbox, 'Units', 'Pixels', 'Position', [1 1+halfpanel_height checkbox_width halfpanel_height]);
            set(obj.StatusText, 'Units', 'Pixels', 'Position', [1+checkbox_width halfpanel_height status_width halfpanel_height]);
            set(obj.OverlayCheckbox, 'Units', 'Pixels', 'Position', [1 1 checkbox_width halfpanel_height]);
            set(obj.OpacitySlider, 'Units', 'Pixels', 'Position', [checkbox_width 1 status_width halfpanel_height]);

            % Setup window/level panel
            windowlevel_text_width = 60;
            windowlevel_editbox_width = 60;
            windowlevel_editbox_height = 19;
            windowlevel_slider_width = max(1, central_panels_width - windowlevel_text_width - windowlevel_editbox_width);
            windowlevel_editbox_position =  1+windowlevel_text_width;
            windowlevel_slider_position = 1+windowlevel_text_width + windowlevel_editbox_width;
            set(obj.LevelText, 'Units', 'Pixels', 'Position', [1 0 windowlevel_text_width halfpanel_height]);
            set(obj.LevelEditbox, 'Units', 'Pixels', 'Position', [windowlevel_editbox_position 0 windowlevel_editbox_width windowlevel_editbox_height]);
            set(obj.LevelSlider, 'Units', 'Pixels', 'Position', [windowlevel_slider_position 1 windowlevel_slider_width halfpanel_height]);
            set(obj.WindowText, 'Units', 'Pixels', 'Position', [1 halfpanel_height windowlevel_text_width halfpanel_height]);
            set(obj.WindowEditbox, 'Units', 'Pixels', 'Position', [windowlevel_editbox_position halfpanel_height-1 windowlevel_editbox_width windowlevel_editbox_height]);
            set(obj.WindowSlider, 'Units', 'Pixels', 'Position', [windowlevel_slider_position 1+halfpanel_height windowlevel_slider_width halfpanel_height]);
        end

        function SetStatus(obj, status_text)
            set(obj.StatusText, 'String', status_text);
        end

        function ImageCheckboxCallback(obj, hObject, ~, ~)
            % Called when show image checkbox is changed
            obj.ViewerPanel.ShowImage = get(hObject, 'Value');
        end

        function OverlayCheckboxCallback(obj, hObject, ~, ~)
            % Called when show overlay checkbox is changed
            obj.ViewerPanel.ShowOverlay = get(hObject, 'Value');
        end

        function WindowSliderCallback(obj, hObject, ~, ~)
            % Window slider
            obj.ViewerPanel.Window = round(get(hObject,'Value'));
        end

        function LevelSliderCallback(obj, hObject, ~, ~)
            % Level slider
            obj.ViewerPanel.Level = round(get(hObject,'Value'));
        end

        function WindowTextCallback(obj, hObject, ~, ~)
            % Window edit box
            new_value = round(str2double(get(hObject, 'String')));

            if isnan(new_value)
                new_value = obj.ViewerPanel.Window;
                set(hObject, 'String', num2str(new_value, '%.6g'));
            else
                obj.ViewerPanel.Window = new_value;
            end
        end

        function LevelTextCallback(obj, hObject, ~, ~)
            % Level edit box
            new_value = round(str2double(get(hObject, 'String')));

            if isnan(new_value)
                new_value = obj.ViewerPanel.Level;
                set(hObject, 'String', num2str(new_value, '%.6g'));
            else
                obj.ViewerPanel.Level = new_value;
            end
        end

        function OrientationCallback(obj, ~, eventdata, ~)
            % Coronal/sagittal/axial orientation
            switch get(eventdata.NewValue, 'Tag')
                case 'Coronal'
                    obj.ViewerPanel.Orientation = GemImageOrientation.XZ;
                case 'Sagittal'
                    obj.ViewerPanel.Orientation = GemImageOrientation.YZ;
                case 'Axial'
                    obj.ViewerPanel.Orientation = GemImageOrientation.XY;
            end
        end

        function ControlsCallback(obj, ~, eventdata, ~)
            % New tool selected
            obj.ViewerPanel.SetControl(get(eventdata.NewValue, 'Tag'));
        end

        function OpacitySliderCallback(obj, hObject, ~, ~)
            % The slider controlling the opacity of the transparent overlay

            obj.ViewerPanel.OverlayOpacity = get(hObject,'Value');
        end

        function WindowLimitsChangedCallback(obj, ~, ~, ~)
            % This methods is called when the window limits have changed
            obj.UpdateWindowLimits;
        end

        function LevelLimitsChangedCallback(obj, ~, ~, ~)
            % This methods is called when the window limits have changed

            obj.UpdateLevelLimits;
        end

        function SelectedControlChangedCallback(obj, ~, ~, ~)
            % Need to resize the control panel as the number of tools may have changed
            obj.SetControl(obj.ViewerPanel.SelectedControl);
            obj.ResizeControlPanel;
        end

        function GuiChangeCallback(obj, ~, ~)
            obj.UpdateGui;
        end

        function ResizeControlPanel(obj)
            control_panel_position = obj.Position;
            control_panel_position(4) = obj.ViewerPanel.ControlPanelHeight;
            obj.Resize(control_panel_position);
        end

        function NewBackgroundImageCallback(obj, ~, ~)
            obj.UpdateGui;
        end

        function NewOverlayImageCallback(obj, ~, ~)
            obj.UpdateGui;
        end

        function NewQuiverImageCallback(obj, ~, ~)
            obj.UpdateGui;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function cy = zcopy ( n, cx, incx, cy, incy )

%*****************************************************************************80
%
%% ZCOPY copies a complex vector X to a vector Y.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer N, the number of elements in CX and CY.
%
%    Input, complex CX(*), the first vector.
%
%    Input, integer INCX, the increment between successive entries of CX.
%
%    Input, complex CY(*), the second vector.
%
%    Input, integer INCY, the increment between successive entries of CY.
%
%    Output, complex CY(*), the second vector, with certain elements
%    copied from CX.
%
  cy(1:incy:1+(n-1)*incy) = cx(1:incx:1+(n-1)*incx);

  return
end


text:
function [I1,wl1,rgb1,g,H1,extra,s,G,S] = estimate_PSF_SRF(I,wl,bbl,rgb,s,extra,options)
%ESTIMATE_PSF_SRF Estimate point spread function (PSF) and spectral
%response function (SRF)
disp('Start estimating PSF and SRF');
show_fig = parse_param(options,'show_fig',0);

I1 = I(:,:,bbl);
wl1 = wl(bbl);
rgb1 = rgb;

Y = reshape_hsi(I1);
X = reshape(rgb1, [size(rgb1,1)*size(rgb1,2), size(rgb1,3)]);
[N,B] = size(Y);
b = size(X,2);
R2 = s + extra*2;
R = prod(R2);

convergence_t = 1e-3;

SRF_range = parse_param(options,'SRF_range','color');
% construct S
[I_sel,bands_sel] = select_relevant_bands(I1,wl1,SRF_range);
wl_sel = wl1(bands_sel)';
S = zeros(B,length(wl_sel));
for i = 1:size(S,2)
    S(bands_sel(i),i) = 1;
end

% construct C
Cs = construct_C(I1,rgb1,s,extra);
C = cat(2,Cs{:});

Y1 = [ones(N,1),Y*S];

% construct D
D = zeros(N*b,R);
for i = 1:N
    D((i-1)*b+1:i*b,:) = X'*Cs{i}';
end

% contruct options
options = [];
options.D = D;
options.Y1 = Y1;
options.D1D = D'*D;
options.D1YI = D'*kron(Y1,eye(b));

g = (1/R) * ones(R,1);
G = g2G(C,g,N);

delta_t0 = 1e-4;
delta_t_g = delta_t0;

errors = [];
for iter = 1:200
    % Update H1
    H1 = solve_for_H1(G*X, Y1, struct('lambda',1e-4));
    options.H1 = H1;

    % Update G
    der_g = calc_der_g(g, options);
    options.der_g = der_g;
    delta_t_g = calc_time_step_adaptive(@eval_obj_fun_g, @update_g, ...
        g, options, delta_t_g, delta_t0);
    g = update_g(g, options, delta_t_g);
    options.g = g;
    G = g2G(C,g,N);

    % Test convergence
    errors(end+1) = eval_obj_fun(options);

    if test_convergence(errors, convergence_t)
        break;
    end

    if mod(iter,100) == 0
        disp(['Process iteration ',num2str(iter)]);
    end
end

if show_fig
    figure('name','PSF'), mesh(reshape(g,R2));
    figure('name','SRF'), plot(H1);
end

function der_g = calc_der_g(g, options)
vecH1 = (options.H1)';
vecH1 = vecH1(:);
der_g = options.D1D * g - options.D1YI * vecH1;

function val = eval_obj_fun_g(g, options)
vecYH = options.Y1 * options.H1;
vecYH = vecYH';
vecYH = vecYH(:);
val = sum((options.D * g - vecYH).^2);

function g_new = update_g(g, options, delta_t)
der_g = options.der_g;

g_new = g - delta_t * der_g;
g_new = project_to_simplex(g_new');
g_new = g_new';

function val = eval_obj_fun(options)
g = options.g;
val = eval_obj_fun_g(g,options);

author: zhouyuanzxcv

repo: Hyperspectral


text:
function triangulation_test23 ( )

%*****************************************************************************80
%
%% TEST23 tests TRIANGULATION_ORDER6_BOUNDARY_EDGE_COUNT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 August 2006
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST23\n' );
  fprintf ( 1, '  TRIANGULATION_ORDER6_BOUNDARY_EDGE_COUNT counts the\n' );
  fprintf ( 1, '    boundary edges in an order 6 triangulation.\n' );

  [ node_num, triangle_num, hole_num ] = ...
    triangulation_order6_example1_size ( );

  [ node_xy, triangle_node, triangle_neighbor ] = ...
    triangulation_order6_example1 ( node_num, triangle_num );

  boundary_edge_num = triangulation_order6_boundary_edge_count ( ...
    triangle_num, triangle_node );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of boundary edges = %d\n', boundary_edge_num );
  fprintf ( 1, '  Correct number =           %d\n', 16 );

  return
end


text:
function summary = reportLineNoise(fid, noiseDetection, numbersPerRow, indent)
% Outputs a summary of line noise removal to file fid and returns a cell array of important messages
    summary = {};
    if isfield(noiseDetection.errors, 'lineNoise') && ...
            ~isempty(noiseDetection.errors.lineNoise)
        summary{end+1} =  noiseDetection.errors.lineNoise;
        fprintf(fid, '%s\n', summary{end});
    end
    if ~isfield(noiseDetection, 'lineNoise')
        summary{end+1} = 'Signal didn''t have line noise removed';
        fprintf(fid, '%s\n', summary{end});
        return;
    end
    lineNoise = noiseDetection.lineNoise;
    if isfield(lineNoise, 'lineNoiseMethod')
        summary{end + 1} = ['Line noise method: ' ...
                             num2str(lineNoise.lineNoiseMethod)];
    else
        summary{end + 1} = 'Line noise method: clean';
    end
    fprintf(fid, '%s\n', summary{end});
    fprintf(fid, 'Sampling frequency Fs: %g Hz\n', lineNoise.Fs);
    fprintf(fid, 'Line noise frequencies:\n');
    printList(fid, lineNoise.lineFrequencies, numbersPerRow, indent);
    fprintf(fid, 'Maximum iterations: %d\n', lineNoise.maximumIterations);
    fprintf(fid, 'Significant frequency p-value: %g\n', lineNoise.p);
    fprintf(fid, '+/- frequency BW for significant peaks (fScanBandWidth): %g\n', ...
        lineNoise.fScanBandWidth);
    fprintf(fid, 'Taper bandwidth: %d Hz\n', lineNoise.taperBandWidth);
    fprintf(fid, 'Taper window size (seconds): %d\n', lineNoise.taperWindowSize);
    fprintf(fid, 'Taper step size (seconds): %d\n', lineNoise.taperWindowStep);
    fprintf(fid, 'Sigmoidal smoothing factor (tau): %g\n',  ...
        lineNoise.tau);
    fprintf(fid, 'Spectral pad factor: %d\n', lineNoise.pad);
    fprintf(fid, 'Analysis frequency interval(fPassBand): [ %g, %g ] Hz\n', ...
        lineNoise.fPassBand);
    if isfield(lineNoise, 'taperTemplate')
        fprintf(fid, 'Taper template: [ %g, %g, %g ]\n', ...
        lineNoise.taperTemplate);
    end
    fprintf(fid, 'Line noise channels (%d channels):\n', ...
        length(lineNoise.lineNoiseChannels));
    printListCompressed(fid, lineNoise.lineNoiseChannels, numbersPerRow, indent);
end

author: VisLab

repo: EEG-Clean-Tools


text:
function r = gitAdd(repo, pname)
r = -1;
if ~exist('repo','var') || isempty(repo)
    repo = pwd;
end
repoFull = filesepStandard_startup(repo,'full');
currdir = pwd;

if ~ispathvalid(pname)
    return
end

ii = 1;
cmds{ii,1} = sprintf('cd %s', repoFull); ii = ii+1;
cmds{ii,1} = sprintf('git add %s', pname); ii = ii+1;

[errs, msgs] = exeShellCmds(cmds, false, true);

if all(errs == 0)
    r = 0;
else
    k = find(errs ~= 0);
    for ii = 1:length(k)
        fprintf('%s\n', msgs{k(ii)});
    end
end

cd(currdir);

author: BUNPC


text:
classdef ARdrone < quadcopter
    properties
        config = 'OUTDOOR'; % Indoor/outdoor dynamics
        % Performance Parameters
        battery_capacity  = 1000;	% Maximum battery capacity (mAh)
        battery_voltage   = 11.1;	% Maximum battery voltage (V)
        flight_time       = 720; 	% Rated flight time (s)
    end
    methods
        % Constructor
        function [this] = ARdrone(varargin)
            % Call the super class
            this@quadcopter(varargin);

            % Create the dynamics of a quadcopter
            this.SENSORS  = this.CreateSENSORS();     % Get the sensor parameters
            this.DYNAMICS = this.CreateDYNAMICS();
            this.GEOMETRY = this.CreateGEOMETRY();
            this.radius   = this.GEOMETRY.length/2;   % Use the "length" to define the radius

            % //////////////// Check for user overrides ///////////////////
            this = this.ApplyUserOverrides(varargin); % Recursive overrides
            % /////////////////////////////////////////////////////////////
        end
        % Setup
        function [this] = setup(this,localXYZVelocity,localXYZrotations)
            % This function calculates the intial state for a quadrotor
            % object.

            % For a state vector defined as x_k:
            % [x,y,z,phi,theta,psi,dx,dy,dz,wx,wy,wz]

            % THIS MODEL OPERATES IN THE GLOBAL FRAME
            %             p0 = this.GetGLOBAL('position');   % True global position
            %             v0 = this.GetGLOBAL('velocity');   % True global velocity
            %             % GET THE ROTATION MATRIX fixed local axes -> rotated global pose
            %             eta0 = OMAS_geometry.eulersToRotationMatrix(localXYZrotations);
            %             % Build an initial (global) state vector
            %             vecR0   = reshape(R0',9,1);        % Defines local vector to global vector
            %             omega0  = zeros(3,1);
            %             x0      = [p0;v0;vecR0;omega0];

            p0 = this.GetGLOBAL('position');   % True global position
            v0 = this.GetGLOBAL('velocity');   % True global velocity
            eta0   = localXYZrotations;
            omega0 = zeros(3,1);
            x0 = [p0;eta0;v0;omega0];

            % ASSIGN THE LOCAL FRD STATE
            this.SetGLOBAL('priorState',x0);
            this.localState = x0;
        end
        % Main
        function [this] = main(this,ENV,varargin)

            % //////////// CHECK FOR NEW INFORMATION UPDATE ///////////////
            % Update the agent with the environmental data
            [this,~,~] = this.GetAgentUpdate(ENV,varargin{1});
            % /////////////////////////////////////////////////////////////

            % The state reference
            desiredPosition = [1;1;0];

            % Call the controller loop
            this = this.Controller(ENV,desiredPosition);

            % /////////////// RECORD THE AGENT-SIDE DATA //////////////////
            this.DATA.inputNames = {'$\dot{x}$ (m/s)','$\dot{y}$ (m/s)','$\dot{z}$ (m/s)',...
                '$p$ (rad/s)','$q$ (rad/s)','$r$ (rad/s)'};
            % Record the control inputs
            this.DATA.inputs(1:length(this.DATA.inputNames),ENV.currentStep) = this.localState(7:end);
        end
    end
    %% ////////////////////// AUXILLARY METHODS ///////////////////////////
    methods
        % Get the state update (using ode45)
        function [x_k_plus] = UpdateLocalState(this,ENV,x_k,y_desired)
            % This function computes the state update for the current agent
            % using the ode45 function.

            % Integrate across the time delta "dt"
            opts = odeset('RelTol',1e-2,'AbsTol',ENV.dt*1E-1);

            [~,Xset] = ode45(@(t,X) this.ClosedLoopDynamics_position(X,y_desired),[0 ENV.dt],x_k,opts);
            %            [~,Xset] = ode45(@(t,X) this.ClosedLoopDynamics_velocity(X,y_desired),[0 ENV.dt],x_k,opts);

            % Assign the integral state
            x_k_plus = Xset(end,:)';
        end
        % ARdrone Dynamics (Closed-loop)
        function [dxdt] = ClosedLoopDynamics_position(this,x_k,p_desired)

            % Extract the current state properties
            %             p_k   = x_k(1:3,1);
            %             eta_k = x_k(4:6,1);
            %             v_k   = x_k(7:9,1);
            %             w_k   = x_k(10:12,1);
            %             x_k = [p_desired(1:3);zeros(3,1);[0;0;psi_desired];zeros(3,1)];

            psi_desired = 0;
            % Define the desired state
            x_desired = [...
                p_desired(1:3);...
                [0;0;psi_desired];...
                zeros(3,1);...
                zeros(3,1)];

            % Extract DYNAMIC parameters
            e3 = this.DYNAMICS.e3;
            g  = this.DYNAMICS.g;
            m  = this.DYNAMICS.m;
            I  = this.DYNAMICS.I;

            % X = [dx dy dz wx wy wz ddx ddy ddz dwx dwy dwz]
            % State matrix
            A = this.DYNAMICS.A;
            A(1:6,7:12) = eye(6);
            A(7:9,4:6)  = g*[ sin(psi_desired), cos(psi_desired), 0;
                             -cos(psi_desired), sin(psi_desired), 0;
                                            0,                0, 0];
%             A_prev = zeros(12);
%             A_prev(1:3,4:6) = eye(3);
%             A_prev(4:6,7:9) = g*[ sin(psi_desired), cos(psi_desired), 0;
%                                  -cos(psi_desired), sin(psi_desired), 0;
%                                                  0,                0, 0];
%             A_prev(7:9,10:12)=eye(3);

            % Input matrix
            B = this.DYNAMICS.B;
            B(7:9,1) = e3/m;
            B(10:12,2:4) = inv(I);

            % LQR control gain
            [K,~,~] = lqr(A,B,eye(12),eye(4));
            % Calculate the error
            e_k  = x_k - x_desired;
            du_k = -K*e_k;
            % Calculate the inputs
            u_k = [m*g;0;0;0] + du_k;

            % Provide the inputs to the open-loop dynamics
            [dxdt] = this.OpenLoopDynamics(x_k,u_k);
        end
        % Quadcopter Dynamics (Open-loop)
        function [dxdt] = OpenLoopDynamics(this,x_k,u_k)
            % Designed for a state vector defined as:
            % x_k = [x,y,z,phi,theta,psi,dx,dy,dz,p,q,r]

            % Extract the state parameters
            eta_k = x_k(4:6,1);
            v_k   = x_k(7:9,1);
            w_k   = x_k(10:12,1);

            % Define the rotation matrix from the euler angles
%             R_k = reshape(x_k(7:15),3,3);
            R_k = OMAS_geometry.eulersToRotationMatrix(eta_k);

            % Extract the inputs
            f   = u_k(1);
            tau = u_k(2:4);

            % Get the constants
            m  = this.DYNAMICS.m;
            g  = this.DYNAMICS.g;
            I  = this.DYNAMICS.I;
            e3 = this.DYNAMICS.e3;

            % nonlinear dynamics based on rotation dynamics
            p_dot    = v_k;
            v_dot    = f/m*R_k*e3-g*e3;
            % Euler angles
            %R_dot    = R_k*skew(w_k);
            %eta_dot  = OMAS_geometry.rotationMatrixToEulers(R_dot);
            eta_dot = w_k;
            %vecR_dot = reshape(R_dot,9,1);

            w_dot    = inv(I)*(tau-skew(w_k)*I*w_k);

            % THE STATE DIFFERENCE
            dxdt = [p_dot;eta_dot;v_dot;w_dot];
        end
        % Global update from the new state
        function [this] = GlobalUpdate(this,x_k_plus)
            % This function updates the global structure from the new state
            % definition: [p_k;eta_k;v_k;omega_k]

            % Sanity Check
            assert(IsColumn(x_k_plus,12),"Expecting a 12DOF state vector.");

            % Define the rotation matrix from the euler angles
            R_k_plus = OMAS_geometry.eulersToRotationMatrix(x_k_plus(4:6,1));
            % Define the new quaternion from R
            q_k_plus = OMAS_geometry.rotationMatrixToQuaternion(R_k_plus');

            % //////////////// UPDATE GLOBAL PROPERTIES ///////////////////
            [this] = this.GlobalUpdate_direct(...
                x_k_plus(1:3),...   % The global position is within the state
                x_k_plus(7:9),...   % The global velocity is within the state
                q_k_plus);          % Append the global quaternion

            % Ensure the local state is re-assigned
            this.localState = x_k_plus;
        end
    end
    methods (Static)
        % Define the open-loop dynamics
        function [dxdt] = ARdrone_OpenLoopDynamics(x_k,u_k)
            %#codegen

            %    This function was generated by the Symbolic Math Toolbox version 8.1.
            %    17-Jul-2018 15:13:41

            PHI_t = x_k(4,:);
            PHI_dot = x_k(10,:);
            PSI_dot = x_k(12,:);
            THETA_t = x_k(5,:);
            THETA_dot = x_k(11,:);
            omega_1 = u_k(1,:);
            omega_2 = u_k(2,:);
            omega_3 = u_k(3,:);
            omega_4 = u_k(4,:);
            x_dot = x_k(7,:);
            y_dot = x_k(8,:);
            z_dot = x_k(9,:);
            t2 = cos(PHI_t);
            t3 = cos(THETA_t);
            t4 = sin(PHI_t);
            t5 = sin(THETA_t);
            t6 = omega_1.^2;
            t7 = sqrt(2.0);
            t8 = omega_2.^2;
            t9 = omega_3.^2;
            t10 = omega_4.^2;
            t11 = t7.*t9.*8.516167950913242e-5;
            t12 = t7.*t10.*8.516167950913242e-5;
            t13 = PSI_dot.^2;
            dxdt = [x_dot;y_dot;z_dot;PHI_dot;THETA_dot;PSI_dot;t5.*(-9.81e2./1.0e2)+THETA_dot.*t4.*y_dot+THETA_dot.*t2.*z_dot-PSI_dot.*t2.*t3.*y_dot+PSI_dot.*t3.*t4.*z_dot;-PHI_dot.*z_dot+t3.*t4.*(9.81e2./1.0e2)+PSI_dot.*t5.*z_dot-THETA_dot.*t4.*x_dot+PSI_dot.*t2.*t3.*x_dot;t6.*(-3.73475e-5)-t8.*3.73475e-5-t9.*3.73475e-5-t10.*3.73475e-5+PHI_dot.*y_dot+t2.*t3.*(9.81e2./1.0e2)-PSI_dot.*t5.*y_dot-THETA_dot.*t2.*x_dot-PSI_dot.*t3.*t4.*x_dot;t11+t12+THETA_dot.*omega_1.*3.231842180365297e-3-THETA_dot.*omega_2.*3.231842180365297e-3+THETA_dot.*omega_3.*3.231842180365297e-3-THETA_dot.*omega_4.*3.231842180365297e-3-t6.*t7.*8.516167950913242e-5-t7.*t8.*8.516167950913242e-5-THETA_dot.^2.*t4.*9.996789383561644e-1-PSI_dot.*THETA_dot.*t2-t3.*t4.*t13+PSI_dot.*THETA_dot.*t2.*t3.*9.996789383561644e-1;-t11+t12+PHI_dot.*PSI_dot-PHI_dot.*omega_1.*3.231842180365297e-3+PHI_dot.*omega_2.*3.231842180365297e-3-PHI_dot.*omega_3.*3.231842180365297e-3+PHI_dot.*omega_4.*3.231842180365297e-3+t6.*t7.*8.516167950913242e-5-t7.*t8.*8.516167950913242e-5-t5.*t13+PHI_dot.*THETA_dot.*t4.*9.996789383561644e-1-PHI_dot.*PSI_dot.*t2.*t3.*9.996789383561644e-1;t6.*4.175141847767905e-4-t8.*4.175141847767905e-4+t9.*4.175141847767905e-4-t10.*4.175141847767905e-4-PHI_dot.*THETA_dot.*1.000321164757521+PHI_dot.*THETA_dot.*t2.*1.000321164757521+PSI_dot.*THETA_dot.*t5.*1.000321164757521+PHI_dot.*PSI_dot.*t3.*t4.*1.000321164757521];
        end
    end
    % Properties
    methods
        % Get the (ARdrone) SENSOR structure
        function [SENSORS]  = CreateSENSORS(this)
            % This function gets the ARdrone 2.0's sensor-specific data and
            % imports them to OMAS's obj.SENSOR structure.

            % Get the default SENSOR structure
            SENSORS = this.SENSORS;
            SENSORS.ultrasound_freq = 40E3;     % Ultrasonic range-finder frequency (Hz)
            SENSORS.ultrasound_range = 6;       % Ultrasonic range-fidner range (m)
            SENSORS.camera_freq = 30;           % Camera capture frequency (fps)
        end
        % Get the (generic) quadcopter dynamic properties
        function [DYNAMICS] = CreateDYNAMICS(this)

            % Import the default properties of a simple quadcopter
            DYNAMICS = this.CreateDYNAMICS_default();

            % ALTER THE DYNAMICS BASED ON THE CONFIGURATION
            switch upper(this.config)
                case 'NONE'
                    DYNAMICS.m = 0.366;
                case 'OUTDOOR'
                    DYNAMICS.m = 0.400; % Arm length equiv : 0.3196
                case 'INDOOR'
                    DYNAMICS.m = 0.436;
                otherwise
                    error('Configuration not recognised');
            end
            % Assign the ARdrone properties
%             DYNAMICS.I = 1.0E-06.*[...
%                 2.8032E+04,0.0000E+00,0.0000E+00;
%                 0.0000E+00,2.8032E+04,0.0000E+00;
%                 0.0000E+00,0.0000E+00,2.8023E+04];
%
            % Control properties
            % Plant matrix
            DYNAMICS.A = zeros(12);
%             DYNAMICS.A(4:6,7:9)   = eye(3)*DYNAMICS.g;
%             DYNAMICS.A(1:3,4:6)   = eye(3);
%             DYNAMICS.A(7:9,10:12) = eye(3);
            % Input matrix
            DYNAMICS.B = zeros(12,4);
%             DYNAMICS.B(4:6,1)     = DYNAMICS.e3/DYNAMICS.m;
%             DYNAMICS.B(10:12,2:4) = inv(DYNAMICS.I);
            % Observation matrix
            DYNAMICS.C = eye(12);
            % Feed-forward matrix
            DYNAMICS.D = eye(4);
        end
        % Get the (ARdrone) GEOMETRY structure
        function [GEOMETRY] = CreateGEOMETRY(this)

            % Parse the geometry associated with the object
            [GEOMETRY] = this.GetObjectGeometry(this);

            % ALTER THE DYNAMICS BASED ON THE CONFIGURATION
            switch upper(this.config)
                case 'NONE'
                    GEOMETRY.width  = 0.450;
                    GEOMETRY.length = 0.290;
                case 'OUTDOOR'
                    GEOMETRY.width  = 0.452;
                    GEOMETRY.length = 0.452; % Arm length equiv : 0.3196
                case 'INDOOR'
                    GEOMETRY.width  = 0.515;
                    GEOMETRY.length = 0.515;
                otherwise
                    error('Configuration not recognised');
            end
        end
    end
end

repo: OpenMAS


text:
function P = subsref(P,String)
%subsref           Overloaded indexing

if length(String)>1
    R = subsref(P,String(1));
    P = subsref(R,String(2:end));
    return
end
switch String.type
    case '.'
        switch String.subs
            case 'Constraints'
                P = P.Constraints;
            case 'Objective'
                P = P.Objective;
            case 'Options'
                P = P.Options;
            case 'maximize'
                maximize(P);
            case 'minimize'
                minimize(P);
            otherwise
                error('Indexing type not supported');
        end
    case '()'
        if isa(String.subs{1},'struct')
            P.Options = String.subs{1};
        elseif isa(String.subs{1},'lmi') | isa(String.subs{1},'constraint')
            P = [P, String.subs{1}];
        else
            error('Indexing type not supported');
        end
    otherwise
        error('Indexing type not supported');
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function write_input_data(off_path, data_path, classes, volume_size, pad_size, angle_inc)
% Put the mesh object in a volume grid and save the volumetric
% represenation file.
% This is the input volumetric data for 3D ShapeNets.
% off_path: root off data folder
% data_path: destination volumetric data folder

phases = {'train', 'test'};

data_size = pad_size * 2 + volume_size;
for c = 1 : length(classes)
    fprintf('writing the %s category\n', classes{c});
    category_path = [off_path '/' classes{c}];
    dest_path = [data_path '/' classes{c} '/' num2str(data_size)];
    if ~exist(dest_path, 'dir')
        mkdir(dest_path);
    end
    % for train and test phases
    for t = 1 : numel(phases)
        phase = phases{t};
        off_list = [category_path '/' phase];
        dest_tsdf_path = [dest_path '/' phase];
        if ~exist(dest_tsdf_path, 'dir')
            mkdir(dest_tsdf_path);
        end
        files = dir(off_list);
        for i = 1 : length(files)
            if strcmp(files(i).name, '.') || strcmp(files(i).name, '..') || files(i).isdir == 1 || ~strcmp(files(i).name(end-2:end), 'off')
                continue;
            end
            filename = [off_list '/' files(i).name];
            for viewpoint = 1 : 360/angle_inc
                destname = [dest_tsdf_path '/' files(i).name(1:end-4) '_' num2str(viewpoint) '.mat'];
                off_data = off_loader(filename, (viewpoint-1)*angle_inc);
                instance = polygon2voxel(off_data, [volume_size, volume_size, volume_size], 'auto');
                instance = padarray(instance, [pad_size, pad_size, pad_size]);
                instance = int8(instance);
                save(destname, 'instance');
            end
        end
    end
end

function offobj = off_loader(filename, theta, axis, stretch)

offobj = struct();
fid = fopen(filename, 'rb');
OFF_sign = fscanf(fid, '%c', 3);
assert(strcmp(OFF_sign, 'OFF') == 1);

info = fscanf(fid, '%d', 3);
offobj.vertices = reshape(fscanf(fid, '%f', info(1)*3), 3, info(1))';
offobj.faces = reshape(fscanf(fid, '%d', info(2)*4), 4, info(2))';

% do some translation and rotation
center = (max(offobj.vertices) + min(offobj.vertices)) / 2;
offobj.vertices = bsxfun(@minus, offobj.vertices, center);
if exist('axis', 'var')
    switch axis
        case 'x',
            offobj.vertices(:,1) = offobj.vertices(:,1) * stretch;
        case 'y',
            offobj.vertices(:,2) = offobj.vertices(:,2) * stretch;
        case 'z',
            offobj.vertices(:,3) = offobj.vertices(:,3) * stretch;
        otherwise,
            error('off_loader axis set wrong');
    end
end
theta = theta * pi / 180;
R = [cos(theta), -sin(theta), 0;
     sin(theta), cos(theta) , 0;
        0      ,    0       , 1];

offobj.vertices = offobj.vertices * R;

% These vertices to define faces should be offset by one to follow the matlab convention.
offobj.faces = offobj.faces(:,2:end) + 1;

fclose(fid);

author: zhirongw


text:
function y = sigma_max( varargin )

%SIGMA_MAX    Maximum singular value.
%   SIGMA_MAX(X) returns the maximum singular value of X. X must be a 2-D
%   matrix, real or complex. SIGMA_MAX(X) is synonymous with NORM(X).
%
%   Disciplined convex programming information:
%       SIGMA_MAX(X) is convex and nonmontonic in X, so X must be affine.

persistent params
if isempty( params ),
    params.nargs     = 1;
    params.args      = [];
    params.empty     = 0;
	params.constant  = @sigma_max_diag;
	params.diagonal  = @sigma_max_diag;
	params.affine    = @sigma_max_aff;
    params.structure = 'svd';
end

try
    y = cvx_matrix_op( params, varargin );
catch exc
    if strncmp( exc.identifier, 'CVX:', 4 ), throw(exc);
    else rethrow(exc); end
end

function y = sigma_max_diag( D )
y = max( abs( D ) );

function z = sigma_max_aff( X )
[ m, n ] = size( X );
cvx_begin sdp
    epigraph variable z nonnegative_
    z * speye(m+n) >= [zeros(m,m),X;X',zeros(n,n)]; %#ok
cvx_end

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

author: yu-jiang

repo: radpbook

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


text:
function HeadModel = in_bst_headmodel(HeadModelFile, ApplyOrient, varargin)
% IN_BST_HEADMODEL: Read a headmodel file in Brainstorm format.
%
% USAGE:  HeadModel = in_bst_headmodel(HeadModelFile, ApplyOrient, FieldsList) : Read the specified fields
%         HeadModel = in_bst_headmodel(HeadModelFile, ApplyOrient)             : Read all the fields
%         HeadModel = in_bst_headmodel(HeadModelFile)                          : Read all the fields and do NOT constrain orientations
%
% INPUT:
%    - HeadModelFile : Absolute or relative path to the headmodel file to read
%    - ApplyOrient   : If 1, return the Gain matrix for a constrained orientation (orientations from GridOrient field)
%    - FieldsList    : List of fields to read from headmodel file

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2009-2016

%% ===== PARSE INPUTS =====
if (nargin < 2) || isempty(ApplyOrient)
    ApplyOrient = 0;
end
% Get protocol folders
ProtocolInfo = bst_get('ProtocolInfo');
% Filename: Relative / absolute
if ~file_exist(HeadModelFile)
    HeadModelFile = bst_fullfile(ProtocolInfo.STUDIES, HeadModelFile);
    if ~file_exist(HeadModelFile)
        error('Head model file not found.');
    end
end
% Specfied fields
if (nargin < 3)
    % Read all fields
    HeadModel = load(HeadModelFile);
else
    % Get fields to read
    FieldsToRead = varargin;
    % If Gain is required, and ApplyOrient is set to 1, we need also the GridOrient field
    if ismember('Gain', FieldsToRead) && ~ismember('GridOrient', FieldsToRead) && ApplyOrient
        FieldsToRead{end + 1} = 'GridOrient';
    end
    % If SurfaceFile required, also add GridLoc
    if ismember('SurfaceFile', FieldsToRead) && ~ismember('GridLoc', FieldsToRead)
        FieldsToRead{end + 1} = 'GridLoc';
    end
    if ismember('GridLoc', FieldsToRead) && ~ismember('SurfaceFile', FieldsToRead)
        FieldsToRead{end + 1} = 'SurfaceFile';
    end
    % If Comment required, also add HeadModelName
    if ismember('Comment', FieldsToRead) && ~ismember('HeadModelName', FieldsToRead)
        FieldsToRead{end + 1} = 'HeadModelName';
    end
    % Read specified files only
    warning off MATLAB:load:variableNotFound
    HeadModel = load(HeadModelFile, FieldsToRead{:});
    warning on MATLAB:load:variableNotFound
end

%% ===== REMOVE CELLS =====
% In previous verisons of Brainstorm, the values were contained in cell lists
% If some values are enclosed in cells, extract them from the cells
fieldsToProcess = {'Gain', 'GridLoc', 'GridOrient'};
% Loop over all fields
for i = 1:length(fieldsToProcess)
    if isfield(HeadModel, fieldsToProcess{i}) && iscell(HeadModel.(fieldsToProcess{i}))
        HeadModel.(fieldsToProcess{i}) = HeadModel.(fieldsToProcess{i}){1};
    end
end

%% ===== READ SOURCE LOCATIONS =====
SurfaceMat = [];
% If GridLoc is a filename
if isfield(HeadModel, 'GridLoc') && ischar(HeadModel.GridLoc)
    SurfaceFile = HeadModel.GridLoc;
    % Load surface file
    SurfaceMat = in_tess_bst(SurfaceFile, 1);
    HeadModel.GridLoc    = SurfaceMat.Vertices;
    HeadModel.SurfaceFile = SurfaceFile;
elseif ~isfield(HeadModel, 'SurfaceFile')
    HeadModel.SurfaceFile = '';
end

% Clean surface filename
if isfield(HeadModel, 'SurfaceFile') && ~isempty(HeadModel.SurfaceFile)
    HeadModel.SurfaceFile = file_short(HeadModel.SurfaceFile);
end

%% ===== READ SOURCE ORIENTATIONS =====
% If GridOrient is empty
if isfield(HeadModel, 'GridOrient') && isempty(HeadModel.GridOrient) && ~isempty(HeadModel.SurfaceFile)
    if isempty(SurfaceMat)
        SurfaceMat = in_tess_bst(HeadModel.SurfaceFile);
    end
    HeadModel.GridOrient = SurfaceMat.VertNormals;
end

%% ===== READ GAIN MATRIX =====
if isfield(HeadModel, 'Gain')
    if ischar(HeadModel.Gain)
        % Get the proper gain covariance matrix
        gainfile = bst_fullfile(bst_fileparts(HeadModelFile), HeadModel.Gain);
        % If the 3 orientations are required, read the _xyz file
        if ~ApplyOrient
            gainfile = strrep(gainfile, '.bin', '_xyz.bin');
        end
        % Read full gain matrix
        HeadModel.Gain = old_read_gain(gainfile);
    else
        HeadModel.Gain = double(HeadModel.Gain);
    end

    % === APPLY ORIENTATIONS ====
    % If the gain matrix has 3 orientations per source, we need to apply the orientation
    if ApplyOrient && (size(HeadModel.Gain,2) == 3 * length(HeadModel.GridOrient))
        % Apply the fixed orientation to the Gain matrix (normal to the cortex)
        HeadModel.Gain = bst_gain_orient(HeadModel.Gain, HeadModel.GridOrient);
    end
end

%% ===== READ COMMENT =====
if isfield(HeadModel, 'HeadModelName') && ~isfield(HeadModel, 'Comment')
    HeadModel.Comment = HeadModel.HeadModelName;
end

%% ===== READ HEAD MODEL TYPE =====
if isfield(HeadModel, 'HeadModelType')
    if ~ismember(HeadModel.HeadModelType, {'surface', 'volume', 'mixed'})
        HeadModel.HeadModelType = 'surface';
    end
end

%% ===== GRID LOC =====
if isfield(HeadModel, 'GridLoc') && ~isempty(HeadModel.GridLoc)
    % Check matrix orientation
    if (size(HeadModel.GridLoc,2) ~= 3)
        HeadModel.GridLoc = HeadModel.GridLoc';
    end
end

%% ===== GRID LOC =====
if isfield(HeadModel, 'GridOrient') && ~isempty(HeadModel.GridOrient)
    % Check matrix orientation
    if (size(HeadModel.GridOrient,2) ~= 3)
        HeadModel.GridOrient = HeadModel.GridOrient';
    end
end

%% ===== ADD MISSING FIELDS =====
if (nargin >= 3)
    for i = 1:length(FieldsToRead)
        if ~isfield(HeadModel, FieldsToRead{i})
            HeadModel.(FieldsToRead{i}) = [];
        end
    end
end

author: brainstorm-tools