text:
function [l c] = path_histogram(G,varargin)
% PATH_HISTOGRAM Compute a histogram of all shortest paths in graph G
%
% [l c] = path_histogram(G) computes all shortest paths in A one at a time
% and forms the histogram of the shortest distances between all vertices.
%
% [l c] = path_histogram(G,struct('sample',N)) uses N random samples instead of an
% exact solution.  (N=0 uses exact solution)
%
% Options:
%   options.weighted: use Dijkstra algorithm for weighted paths [{0} | 1]
%   options.intcount: use integer histogram bins [0 | {1}]
%   options.nbins: number of bins when intcount=0 [{50} | positive integer]
%   options.sample: use N iid samples instead of all sources
%     [{0} | positive integer]
%   options.verbose: output status for long runs [{0} | 1]
%
% Example:
%  load('graphs/cs-stanford.mat');
%  [l c] = path_histogram(A);
%  bar(l,c,[min(l),max(l)],'hist');

%
% David Gleich
% Copyright, Stanford University, 2008
%

% History
% 2008-03-14: Initial coding by David

options=struct('weighted',0,'intcount',1,'nbins',50','sample',0,'istrans',0,...
    'verbose',0);
if ~isempty(varargin), options = merge_structs(varargin{1}, options); end

if ~options.istrans, G = G'; end
bfsoptions=struct('istrans',1);

n=size(G,1); verb=options.verbose;

N = options.sample;
if N==0, vs=1:n; N = n;
else vs=ceil(n*rand(options.sample,1));
end

% sample from all vertices
H=sparse(n,1); t0=clock; p=1; Np=min(N,100);
for vi=1:N
    if verb && vi==(p*floor(N/Np) + max(mod(N,Np)-(Np-p),0))
        dt=etime(clock,t0); p=p+1;
        fprintf('  %5.1f : %9i of %9i; etr=%7f sec\n', ...
            100*(p-1)/Np, vi, N, N*dt/vi-dt);
    end
    di=bfs(G,vi,bfsoptions);
    di=di(di>0);             % only use reachable points
    H=H+accumarray(di,1,[n 1],@sum,0,true);
    bfsoptions.nocheck=1;    % don't repeat options check
end
l = find(H); c = nonzeros(H);

end

% internal copy of merge_structs
function S=merge_structs(A,B)
S = A; fn = fieldnames(B);
for ii = 1:length(fn), if (~isfield(A, fn{ii})), S.(fn{ii}) = B.(fn{ii}); end, end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [isConsistent, m, model] = checkStoichiometricConsistency(model, printLevel, method)
% Verification of stoichiometric consistency by checking for at least one
% strictly positive basis in the left nullspace of `S`.
% If `S` is not stoichiometrically consistent detect conserved and unconserved
% metabolites, by returning a maximal conservation vector, which is
% a non-negative basis with as many strictly positive entries as possible.
% This omits rows of `S` that are entirely zero, when any exchange reactions
% are removed.
% The strictly positive and zero entries in `m` correspond to the conserved
% and unconserved metabolites respectively.
%
% Verification of stoichiometric consistency is as initially described in:
% `A. Gevorgyan, M. G. Poolman, and D. A. Fell
% Detection of stoichiometric inconsistencies in biomolecular
% models. Bioinformatics, 24(19):2245–2251, 2008`.
%
% Detection of conserved and unconserved metabolites based on a new
% implementation by Nikos Vlassis & Ronan Fleming.
%
% USAGE:
%
%    [isConsistent, m, model] = checkStoichiometricConsistency(model, printLevel, method)
%
% INPUT:
%    model:         structure with fields:
%
%                     * model.S - `m` x `n` stoichiometric matrix
%                     * model.mets - if exists, but `SIntRxnBool does not exist, then
%                       `findSExRxnBool` will attempt to identify exchange
%                       reactions. (optional)
%                     * model.SIntMetBool - `m` x 1 boolean vector indicating the metabolites that
%                       are thought to be exclusively involved in non-exchange
%                       reactions (optional)
%                     * model.SIntRxnBool  `n` x 1 boolean vector indicating the reactions that are
%                       thought to be non-exchange reactions. If this is not present,
%                       then we will attempt to identify such reactions,
%                       `model.rxns` exists, otherwise, we will assume all
%                       reactions are supposed to be non-exchange reactions. (optional)
%
% OPTIONAL INPUTS:
%    printLevel:    {(0), 1}
%    method:        structure with fields:
%
%                     * method.interface - {('SDCCO'),'LP', 'MILP', 'DCCO'} interface called to do the consistency check
%                     * method.solver - {(default solver as specified by CBT_LP_SOLVER), or any other CBT compatible LP solver}
%                     * method.param - solver specific parameter structure
%
% OUTPUTS:
%    isConsistent:        Solver status in standardized form:
%
%                     * 1 - Optimal solution (Stoichiometrically consistent)
%                     * 2 - Unbounded solution (Should never happen)
%                     * 0 - Infeasible (Stoichiometrically INconsistent)
%                     * -1 - No solution reported (timelimit, numerical problem etc)
%
%    m:             `m` x 1 strictly positive vector in left nullspace
%                   (empty if it does not exist)
%    model:         structure with fields:
%
%                     * .SConsistentMetBool - m x 1 boolean vector indicating metabolites involved
%                       in the maximal consistent vector
%                     * .SConsistentRxnBool - `n` x 1 boolean vector non-exchange reaction involving
%                       a stoichiometrically consistent metabolite
%
% .. Author: - Ronan Fleming   2012 initial coding
%                              2013 update with detection of conserved metabolites based on an algorithm by Nikos Vlassis.
%                              2014 update to omit trivial rows corresponding to S row that are all zero
%                              2022 interface with findStoichConsistentSubset

if ~exist('printLevel','var')
    printLevel=0;
end

resetSolver = 0;
if exist('method','var')
    if ~isfield(method,'interface')
        method.interface = 'SDCCO';
    end
    if isfield(method,'solver')
        global CBT_LP_SOLVER
        oldLPSolver = CBT_LP_SOLVER;
        resetSolver = 1;
        %set the solver and solver parameters
        solverOK = changeCobraSolver(method.solver,'LP');
    else
        global CBT_LP_SOLVER
        method.solver=CBT_LP_SOLVER;
    end
else
    method.interface='SDCCO';
    global CBT_LP_SOLVER
    method.solver=CBT_LP_SOLVER;
end

%set parameters according to feastol
feasTol = getCobraSolverParams('LP', 'feasTol');
epsilon = feasTol*10;

[nMet,nRxn]=size(model.S);
if ~isfield(model,'mets')
    %assume all reactions are internal
    model.SIntRxnBool=true(nRxn,1);
    intR='';
else
    if ~isfield(model,'SIntRxnBool') || ~isfield(model,'SIntMetBool')
        %Requires the openCOBRA toolbox
        model=findSExRxnInd(model);
        intR='internal reaction ';
    else
        intR='';
    end
end

% Check the stoichiometric consistency of the network by
% solving the following linear problem
%       min sum(m_i)
%           s.t     S'*m = 0
%                   m >= 1
% where l  is is a  mx1 vector of the molecular mass of m molecular species
N=model.S(:,model.SIntRxnBool);
LPproblem.A=N';
LPproblem.b=zeros(size(LPproblem.A,1),1);
LPproblem.lb=ones(size(LPproblem.A,2),1);
LPproblem.ub=inf*ones(size(LPproblem.A,2),1);
LPproblem.c=1*ones(size(LPproblem.A,2),1);
LPproblem.osense=1;
LPproblem.csense(1:size(LPproblem.A,1),1)='E';

%Requires the openCOBRA toolbox
solution = solveCobraLP(LPproblem,'printLevel',printLevel);

%OUTPUT
% solution Structure containing the following fields describing a LP
% solution
%  full     Full LP solution vector
%  obj      Objective value
%  rcost    Reduced costs
%  dual     Dual solution
%  solver   Solver used to solve LP problem
%
%  stat     Solver status in standardized form
%            1   Optimal solution
%            2   Unbounded solution
%            0   Infeasible
%           -1   No solution reported (timelimit, numerical problem etc)

isConsistent=solution.stat;

% If the network is not stoichiometrically consistent then one maximizes
% the number of  positive component of the molecular masses vector
if isConsistent~=1
    switch method.interface
        case 'none'
            warning(['Stoichiometrically INconsistent ' intR 'stoichiometry.']);
        case 'SDCCO'
            tic
            massBalanceCheck=0;
            fileName=[];
%             [SConsistentMetBool, SConsistentRxnBool, SInConsistentMetBool, SInConsistentRxnBool, unknownSConsistencyMetBool, unknownSConsistencyRxnBool, model, stoichConsistModel] =...
%                 findStoichConsistentSubset(model, massBalanceCheck, printLevel, fileName, epsilon)
             [~, ~, ~, ~, ~, ~, model, ~] = findStoichConsistentSubset(model, massBalanceCheck, printLevel,fileName, epsilon);
            %TODO
            m = 1*model.SConsistentMetBool;
            solution = [];
            timetaken=toc;
        case 'cvx'
            cvx_solver(method.solver)

            [nMet,~]=size(N);
            %maximal conservation vector
            %cvx code from Nikos Vlassis
            tic
            cvx_begin quiet

            variable m(nMet);
            variable z(nMet);

            maximize( ones(1,nMet) * z );

            z>=0; z<=epsilon;

            m>=z; m<=(1/epsilon);

            N'*m==0;

            cvx_end
            timetaken=toc;
            if any(isnan(m))
                error('NaN in maximal conservation vector')
            end
            %boolean indicating metabolites involved in the maximal consistent vector
            model.SConsistentMetBool=m>epsilon & model.SIntMetBool;
        case 'LP'

            [nMet,~]=size(N);

            % Solve the linear problem
            %   max sum(z_i)
            %       s.t S'*m = 0
            %           z <= m
            %           0 <= m <= 1/epsilon
            %           0 <= z <= epsilon
            nInt=nnz(model.SIntRxnBool);
            LPproblem.A=[N'      , sparse(nInt,nMet);
                         speye(nMet),      -speye(nMet)];

            LPproblem.b=zeros(nInt+nMet,1);

            LPproblem.lb=[zeros(nMet,1);zeros(nMet,1)];
            LPproblem.ub=[ones(nMet,1)*(1/epsilon);ones(nMet,1)*epsilon];

            LPproblem.c=zeros(nMet+nMet,1);
            LPproblem.c(nMet+1:2*nMet,1)=1;
            LPproblem.osense=-1;
            LPproblem.csense(1:nInt,1)='E';
            LPproblem.csense(nInt+1:nInt+nMet,1)='G';

            %Requires the COBRA toolbox
            tic
            if isfield(method,'param')
                solution = solveCobraLP(LPproblem,'printLevel',printLevel-1,method.param);
            else
                solution = solveCobraLP(LPproblem,'printLevel',printLevel-1);
            end
            timetaken=toc;

            if solution.stat==1
                m=solution.full(1:nMet,1);
                z=solution.full(nMet+1:end,1);
                if isfield(model,'SIntMetBool') && 0
                    %boolean indicating metabolites involved in the maximal consistent vector
                    model.SConsistentMetBool=m>epsilon & model.SIntMetBool;
                else
                    %boolean indicating metabolites involved in the maximal consistent vector
                    model.SConsistentMetBool=m>epsilon;
                end
            else
                disp(solution)
                error('solve for maximal conservation vector failed')
            end

            %mosek
            algorithm=solution.algorithm;
            %change back the solver
%            solverOK = changeCobraSolver(oldSolver,'LP');
        case 'DCCO'

            tic
            [~, ~, solution] = maxCardinalityConservationVector(N);
            timetaken=toc;
            method.solver='cappedL1';
            if solution.stat==1
                m = solution.l;
                %dummy z
                z = zeros(nMet,1);
                if isfield(model,'SIntMetBool')  && 0
                    %boolean indicating metabolites involved in the maximal consistent vector
                    model.SConsistentMetBool=m>epsilon & model.SIntMetBool;
                else
                    %boolean indicating metabolites involved in the maximal consistent vector
                    model.SConsistentMetBool=m>epsilon;
                end
            else
                disp(solution)
                error('solve for maximal conservation vector failed')
            end
        case 'MILP'
            [nMet,~]=size(N);

            % Solve the MILP problem
            %   max sum(z_i)
            %       s.t S'*m = 0
            %           z <= m
            %           z binary
            nInt=nnz(model.SIntRxnBool);
            MILPproblem.A=[N'      , sparse(nInt,nMet);
                         speye(nMet),      -speye(nMet)];

            MILPproblem.b=zeros(nInt+nMet,1);

            MILPproblem.lb=[zeros(nMet,1);zeros(nMet,1)];
            MILPproblem.ub=[ones(nMet,1)*(1/epsilon);ones(nMet,1)];

            MILPproblem.c=zeros(nMet+nMet,1);
            MILPproblem.c(nMet+1:2*nMet,1)=1;
            MILPproblem.osense=-1;
            MILPproblem.csense(1:nInt,1)='E';
            MILPproblem.csense(nInt+1:nInt+nMet,1)='G';
            MILPproblem.vartype(1:nMet,1)='C';
            MILPproblem.vartype(nMet+1:nMet+nMet,1)='B';
            MILPproblem.x0 = zeros(nMet+nMet,1);

%             solverOK = changeCobraSolver('gurobi','MILP');
            %Requires the COBRA toolbox
            tic
            if isfield(method,'param')
                solution = solveCobraMILP(MILPproblem,'printLevel',printLevel-1,method.param);
            else
                solution = solveCobraMILP(MILPproblem,'printLevel',printLevel-1);
            end
            timetaken=toc;

            if solution.stat==1
                m=solution.full(1:nMet,1);
                z=solution.full(nMet+1:end,1);
                if isfield(model,'SIntMetBool')  && 0
                    %boolean indicating metabolites involved in the maximal consistent vector
                    model.SConsistentMetBool=m>epsilon & model.SIntMetBool;
                else
                    %boolean indicating metabolites involved in the maximal consistent vector
                    model.SConsistentMetBool=m>epsilon;
                end
            else
                disp(solution)
                error('solve for maximal conservation vector failed')
            end
        case 'maxEnt'
            %does not work very well
            tic
            m=maxEntConsVector(N,printLevel);
            timetaken=toc;
            if isfield(model,'SIntMetBool')  && 0
                %boolean indicating metabolites involved in the maximal consistent vector
                model.SConsistentMetBool=m>epsilon & model.SIntMetBool;
            else
                %boolean indicating metabolites involved in the maximal consistent vector
                model.SConsistentMetBool=m>epsilon;
            end
            z=zeros(nMet,1);
        otherwise
            error(['unregognised method.interface = ' method.interface]);
    end
    if any(m < -feasTol)
        error('m should be greater than or equal to zero')
    end
    m(m<0)=0;
    if printLevel>0
        if isfield(method,'param')
            fprintf('%s%s%s%s%s%s%s%g%s\n','Maximal conservation vector, using ', method.interface, ' ', method.solver,' ',algorithm,', in time ',timetaken,' sec.')
        else
            fprintf('%s%s%s%s%s%g%s\n','Maximal conservation vector, using ', method.interface, ' ', method.solver,', in time ',timetaken,' sec.')
        end
        fprintf('%10f%s\n',ones(1,nMet) * m,' = Optimal objective (i.e. 1''*m)')
        fprintf('%10d%s\n', nnz(model.SConsistentMetBool),' = Number of stoichiometrically consistent rows')
        fprintf('%10g%s\n',norm(m'*N),' = || S''*m ||_inf for non-exchange reactions of S')
    end
else
    m=solution.full;
    %The only consistent rows are those corresponding to non-exchange reactions
    model.SConsistentMetBool=model.SIntMetBool;
    if printLevel>0
        fprintf('%s\n','--- Summary of stoichiometric consistency ----')
        fprintf('%6s\t%6s\n','#mets','#rxns')
        fprintf('%6u\t%6u\t%s\n',nMet,nRxn,' totals.')
        fprintf('%6u\t%6u\t%s\n',nnz(~model.SIntMetBool),nnz(~model.SIntRxnBool),' heuristically external.')
        fprintf('%6u\t%6u\t%s\n',nnz(model.SIntMetBool),nnz(model.SIntRxnBool),' heuristically internal:')
        fprintf('%6u\t%6u\t%s\n',nnz(model.SConsistentMetBool),nnz(model.SIntRxnBool),' ... of which are stoichiometrically consistent.')
    end
end

isConsistent = all(model.SConsistentMetBool==1);

if ~isfield(model,'SConsistentRxnBool')

    %OLD - incorrect way July 14th 2016 - Ronan.
    % %find every non-exchange reaction involving a stoichiometrically consistent metabolite
    % model.SConsistentRxnBool =(sum(model.S(model.SConsistentMetBool,:)~=0,1)~=0)';
    % model.SConsistentRxnBool(~model.SIntRxnBool)=0;

    %corresponding reactions exclusively involving consistent metabolites
    model.SConsistentRxnBool = ~any(model.S(~model.SConsistentMetBool, :), 1)' & model.SIntRxnBool;

    model.SConsistentRxnBool = getCorrespondingCols(model.S,model.SConsistentMetBool,model.SIntRxnBool,'inclusive');

end

if resetSolver
    %reset the solver
    solverOK = changeCobraSolver(oldLPSolver,'LP');
end

author: opencobra

repo: cobratoolbox


text:
function vec=DelaunayVar(TDB1,TDB2)
%%DELAUNAYVAR Get the Delaunay variables that are needed for tide models
%             involving Doodson numbers using the equations in the IERS
%             conventions.
%
%INPUTS: T The time as measured in Julian centuries Barycentric dynamical
%          time (TDB).
%
%OUTPUTS: vec A vector of the Delaunay variables in the order with all
%             units in RADIANS.
%             vec(1) Mean anomaly of the moon: l
%             vec(2) Mean anomaly of the sun: l'
%             vec(3) F=L-Omega, The mean longitude of the moon minus the
%                    mean longitude of the ascending node of the moon.
%             vec(4) Mean Elongation of the moon from the sun: D
%             vec(5) Mean Longitude of the Ascending Node of the Moon:
%                    Omega
%
%This function implements Equation 5.43 in Section 5.7.2 of [1]. The units
%in those equations are degrees (and arcseconds). However, in Equation 6.8
%in Section 6.2, it looks like the angles should be in radians, since one
%is using trigonometric functions. Thus, the return value here is put into
%radians to be convenient.
%
%REFERENCES:
%[1] G. Petit and B. Luzum, IERS Conventions (2010), International Earth
%    Rotation and Reference Systems Service Std. 36, 2010.
%
%March 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    %The time as measured in Julian centuries Barycentric dynamical time
    %(TDB). The precision of the model is low enough that it shouldn't
    %matter that the two parts are added together.
    T=(TDB1+TDB2)/36525;

    vec=zeros(5,1);

    %The conversion from arcseconds to degrees
    sec2Deg=1/3600;

    %Mean anomaly of the moon: l
    vec(1)=134.96340251+sec2Deg*(1717915923.2178*T+31.8792*T^2+0.051635*T^3-0.00024470*T^4);
    %Mean anomaly of the sun: l'
    vec(2)=357.52910918+sec2Deg*(129596581.0481*T-0.5532*T^2+0.000136*T^3-0.00001149*T^4);
    %L-Omega, The mean longitude of the moon minus the mean longitude of the
    %ascending node of the moon.
    vec(3)=93.27209062+sec2Deg*(1739527262.8478*T-12.7512*T^2-0.001037*T^3+0.00000417*T^4);
    %Mean Elongation of the moon from the sun: D
    vec(4)=297.85019547+sec2Deg*(1602961601.2090*T-6.3706*T^2+0.006593*T^3-0.00003169*T^4);
    %Mean Longitude of the Ascending Node of the Moon: Omega
    vec(5)=125.04455501+sec2Deg*(-6962890.5431*T+7.4722*T^2+0.007702*T^3-0.00005939*T^4);

    %Convert from degrees to radians.
    vec=mod(vec*(pi/180),2*pi);
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function node_xy = grid_nodes_01 ( x_num, y_num )

%*****************************************************************************80
%
%% GRID_NODES_01 returns an equally spaced rectangular grid of nodes in the unit square.
%
%  Example:
%
%    X_NUM = 5
%    Y_NUM = 3
%
%    NODE_XY =
%    ( 0, 0.25, 0.5, 0.75, 1, 0,   0.25, 0.5, 0.75, 1,   0, 0.25, 0.5, 0.75, 1;
%      0, 0,    0,   0,    0, 0.5, 0.5,  0.5, 0.5,  0.5, 1, 1.0,  1.0, 1.0,  1 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 May 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X_NUM, Y_NUM, the number of nodes in the X and Y directions.
%
%    Output, real NODE_XY(2,X_NUM*Y_NUM), the coordinates of the nodes.
%
  node_num = x_num * y_num;

  node_xy(1:2,1:node_num) = 0.0;

  if ( x_num == 1 )
    node_xy(1,1:node_num) = 0.5;
  else
    for i = 1 : x_num
      node_xy(1,i:x_num:i+(y_num-1)*x_num) = ( i - 1 ) / ( x_num - 1 );
    end
  end

  if ( y_num == 1 )
    node_xy(2,1:node_num) = 0.5;
  else
    for j = 1 : y_num
      node_xy(2,1+(j-1)*x_num:j*x_num) = ( j - 1 ) / ( y_num - 1 );
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [conds,labels,stimData] = StimCSV_Reader()
    %Open File
%     filename = "StimCSV_Sample_Error.csv";
    [file,path,~] = uigetfile({'*.csv','CSV File (*.csv)'});
    filename = [path file];
    stimFileID = fopen(filename);

    %Initial Variables
    i = 1;
    stimNameIndicator = true;
    stimLabelIndicator = true;
    labels = cell(1,1);
    data = [];
    stimData = cell(1,1);

    %Get Line
    linestr = fgetl(stimFileID);

    while(linestr ~= -1)
        strCell = strsplit(linestr,',');
        strCell = strCell(~cellfun(@isempty,strCell));
        if isempty(strCell) %Update Stim (Finalize)
            try
                stimData{:,i} = data;
            catch
                errordlg('Invalid Stim CSV Syntax. Please check reference file for guidance');
                return;
            end
            i = i+1;
            stimNameIndicator = true;
            stimLabelIndicator = true;
            data = [];
        elseif stimNameIndicator %Update Stim Name
            try
                conds{i} = strCell(1);
            catch
                errordlg('Invalid Stim CSV Syntax. Please check reference file for guidance');
                return;
            end
            stimNameIndicator = false;
        elseif stimLabelIndicator %Update Stim Label
            labels{i} = strCell;
            stimLabelIndicator = false;
        else %Update Stim Data
            try
                data = [data; strCell];
            catch
                errordlg('Invalid Stim CSV Syntax. Please check reference file for guidance');
                return;
            end
        end
        %Read Next Line
        linestr = fgetl(stimFileID);
    end

    %Final Stim Update
    stimData{:,i} = data;

    %Close File
    fclose(stimFileID);
end

author: BUNPC


text:
function view = makePointROI(view, pt, createNew)
% Create a region of interest comprising a single voxel.
%
% view = makePointROI([view], [pt=get from UI], [createNew=1]);
%
% view: mrVista view structure. [Default: current view]
%
% pt: 3-element vector specifying the point. In inplane, this is
%	  [row col slice] of the coordinates, relative to the inplane anatomy;
%	  in gray/volume views, this is [axi cor sag] slice specification; in
%	  flat views, it's [x y hemisphere].
%	  [Default: for 3-view windows, take from the crosshairs position. For
%	  other views, prompts the user to type in the coords.]
%
% createNew: flag to make a new ROI with the point. [default 1] Otherwise,
% will add the point to the view's selected ROI -- much like addROIpoints,
% but allows for a more precise way of specifying the coords.
%
% ras, 09/2008.
if notDefined('view'),			view = getCurView;				end
if notDefined('createNew'),		createNew = 1;					end
if notDefined('pt')
	if checkfields(view, 'loc')
		pt = view.loc;
	else
		pt = inputdlg({'Enter the point coordinates:'}, mfilename, 1);
		pt = str2num(pt{1});
	end
end

% check that the point is valid
if length(pt)~=3
	error('point must be a 3-element vector.')
end

if createNew==1
	view = newROI(view, sprintf('Point %s', num2str(pt)), [], [], pt(:));

else
	if isempty(view.ROIs) || view.selectedROI < 1
		error('No existing ROIs in view.')
	end
	roi = view.ROIs(view.selectedROI);
	if isempty(roi.coords)
		roi.coords = pt(:);
	else
		roi.coords(:,end+1) = pt(:);
	end
	view.ROIs(view.selectedROI) = roi;
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef findquar < ZmapHGridFunction
    % description of this function
    %
    % in the function that generates the figure where this function can be called:
    %
    %     % create some menu items...
    %     h=sample_ZmapFunction.AddMenuItem(hMenu,@catfn) %create subordinate to menu item with handle hMenu
    %     % create the rest of the menu items...
    %
    %  once the menu item is clicked, then sample_ZmapFunction.interative_setup(true,true) is called
    %  meaning that the user will be provided with a dialog to set up the parameters,
    %  and the results will be automatically calculated & plotted once they hit the "GO" button
    %

    properties
        inDaytime         (24,1) logical    = false(24,1) % true for hours that are "daytime" hours
        localNoonEstimate (1,1)  double     = 12; % estimate time where local noon is. used with dayLength
        dayLength         (1,2)  double   {mustBePositive}    = [4, 6] % hours BEFORE noon to hours AFTER noon
    end

    properties(Constant)
        PlotTag         = 'QuarryRatios'
        ReturnDetails   = cell2table({...VariableNames, VariableDescriptions, VariableUnits
            'day_night_ratio',      'Day-Night event ratio', '';
            'day_night_ratio_norm', 'Day-Night event ratio (normalized by hrs in day)', '';
            'n_day',                'Number of events during day','';
            'n_night',              'Number of events during night',''...
            }, 'VariableNames', {'Names','Descriptions','Units'})

        CalcFields      = {'day_night_ratio','day_night_ratio_norm','n_day','n_night'}
        DayColor        = [0.8 0.8 0.2] % for histogram
        NightColor      = [0.1 0.0 0.6] % for histogram

        ParameterableProperties = ["inDaytime" "localNoonEstimate" "dayLength" "NodeMinEventCount"];
        References="";
    end

    methods
        function obj=findquar(zap,varargin) %CONSTRUCTOR
            % create a [...]

            obj@ZmapHGridFunction(zap, 'day_night_ratio_norm');
            report_this_filefun();

            obj.CalcLocalNoon();
            obj.NodeMinEventCount = 20;
            % set the deafult days & nights based on  local noon and the "day" length
            dayStart = mod(obj.localNoonEstimate - obj.dayLength(1),24);
            dayEnd = mod(obj.localNoonEstimate + obj.dayLength(2),24);
            eachHour = 0:23;
            if dayStart < dayEnd
                obj.inDaytime = dayStart <= eachHour  & eachHour < dayEnd;
            else
                obj.inDaytime = dayStart <= eachHour | eachHour < dayEnd;
            end

            obj.parseParameters(varargin);
            obj.StartProcess();
        end

        function InteractiveSetup(obj)
            % allow user to determine grid and selection parameters

            obj.InteractiveSetup_part2();

        end

        function InteractiveSetup_part2(obj)
            % allow user to define which hours are in a "day"

            fifhr = figure(...
                'Name','Daytime (explosion) hours',...
                'NumberTitle','off', ...
                'NextPlot','new', ...
                'units','points',...
                'Visible','on', ...
                'Tag','fifhr',...
                'Position',[ 100 100 500 650]);
            axis off
            ax = gca;

            uicontrol(fifhr,'Style','text','String','Detect Quarry Events',...
                'FontWeight','Bold','FontSize',14,'Units','points','Position',[50 620 300 20]);

            set(ax,'NextPlot','add');
            hax = axes(fifhr,'Units','points','pos', [50 320 300 270]);
            eventHours = obj.ToHourlyCategorical(obj.RawCatalog.Date.Hour);

            dayHist = histogram(hax,eventHours,'DisplayName','day','FaceColor',obj.DayColor);
            set(hax,'NextPlot','add');

            nightHist = histogram(hax,eventHours,'DisplayName','night','FaceColor', obj.NightColor);
            title(' Select the daytime hours and then "GO"')
            [X,N,B] = histcounts(obj.RawCatalog.Date.Hour,-0.5:1:24.5);
            %[X,~,B] = histcounts(eventHours);
            %[X,~] = hist(obj.RawCatalog.Date.Hour,-0.5:1:24.5);

            xlabel(hax,'Hr of the day')
            ylabel(hax,'Number of events per hour')

            evsel=EventSelectionChoice(fifhr,'evsel', [40,100], obj.EventSelector);

            chkpos = @(n)[.80 1-n/28-0.03 .17 1/26];
            for i = 24:-1:1
                hHourly(i)=uicontrol(fifhr,'Style','checkbox',...
                    'string',[num2str(i-1) ' - ' num2str(i) ],...
                    'Units','normalized',...
                    'Position',chkpos(i),'tag',num2str(i),...
                    'Callback', @(~,~)cb_flip(i));
            end

            obj.CalcLocalNoon();
            dayStart = mod(obj.localNoonEstimate - obj.dayLength(1),24);
            dayEnd = mod(obj.localNoonEstimate + obj.dayLength(2),24);
            eachHour = 0:23;
            if dayStart < dayEnd
                obj.inDaytime = dayStart <= eachHour  & eachHour < dayEnd;
            else
                obj.inDaytime = dayStart <= eachHour | eachHour < dayEnd;
            end
            % turn on checkboxes according to their percentile score
            %  idx = X(1:end-1) > prctile2(X,60); % then set values from this

            for hr = 1:24
                hHourly(hr).Value = obj.inDaytime(hr);
            end
            update_histograms();
            %nightHist.Data(obj.inDaytime)=nan;
            legend(hax,'show');

            if isempty(findobj(fifhr, 'Tag', 'quarryinfo'))
                add_menu_divider();
                uimenu(fifhr,'Label','Info',...
                MenuSelectedField(), @(~,~)web(['file:' ZmapGlobal.Data.hodi '/help/quarry.htm'], '-browser'),...
                 'tag', 'quarryinfo');
            end

            uicontrol(fifhr,'style','pushbutton','String','GO','Callback',@cb_go,...
            'Position',[330 10 60 25]);

            uicontrol(fifhr,'style','pushbutton','String','Cancel','Callback',@(~,~)close,...
            'Position',[400 10 60 25]);

            function cb_flip(i)
                obj.inDaytime(i) = ~obj.inDaytime(i);
                update_histograms();
            end
            function update_histograms()
                dayCats = obj.ToHourlyCategorical(find(obj.inDaytime)-1);
                nightCats = obj.ToHourlyCategorical(find(~obj.inDaytime)-1);
                dayHist.Data = eventHours(ismember(eventHours, dayCats));
                nightHist.Data = eventHours(ismember(eventHours, nightCats));
            end
            function cb_go(~,~)

                obj.EventSelector = EventSelectionParameters.fromStruct(evsel.toStruct());
                obj.inDaytime=logical([hHourly.Value]); %same as idx
                close;
                obj.doIt();
            end

        end

        function results=Calculate(obj)
            % once the properties have been set, either by the constructor or by interactive_setup

            close(findobj('Tag','fifhr'));

            ld = sum(obj.inDaytime);

            assert(ld ~= 0, 'No daytime hours chosen. This calculation will have no meaning.');
            assert(ld ~= 24, 'No nighttime hours chosen. This calculation will have no meaning.');

            ln = sum(~obj.inDaytime);
            daynight_hr_ratio=ld/ln;

            % loop over all points in polygon. Evaluated for earthquakes that may extend outside
            % the points.
            %{
            [valueMap,nEv,r]=gridfun(@calculate_day_night_ratio, obj.RawCatalog, mygrid, obj.EvtSel);
            bvg=[valueMap(mygrid.ActivePoints), mygrid.ActiveGrid, r(mygrid.ActivePoints)];
            %}

            obj.gridCalculations(@calculate_day_night_ratio);

            if nargout
                results=obj.Result.values;
            end

            % plot the results
            % obj.oldratios and valueMap (initially ) is the b-value matrix

            %obj.oldratios = valueMap;

            % results of the calculation should be stored in fields belonging to obj.Result

            %obj.Result.bvg=bvg;
            %obj.Result.msg='bvg is [daynightratios x y maxdist_km]';
            %obj.Result.valueMap=valueMap;

            function val = calculate_day_night_ratio(catalog)
                hrofday= hour(catalog.Date);
                nDay= sum(obj.inDaytime(hrofday+1));
                nNight = catalog.Count - nDay;
                myratio = (nDay / nNight);
                normalizedRatioByHoursInDay = myratio * daynight_hr_ratio;
                val=[myratio normalizedRatioByHoursInDay nDay nNight];
            end
        end

        function ModifyGlobals(obj)
            % change the ZmapGlobal variable, if appropriate
            % obj.ZG.SOMETHING = obj.Result.SOMETHING
        end

        function CalcLocalNoon(obj)
            catMedianLatitude = median(obj.RawCatalog.Longitude);
            obj.localNoonEstimate = round( catMedianLatitude  / 15) + 12;
        end

    end %methods

    methods(Static)
        function h = AddMenuItem(parent, zapFcn, varargin)
            % create a menu item that will be used to call this function/class
            label = 'Find Quarry Events';
            h = uimenu(parent, 'Label', label,...
                'MenuSelectedFcn', @(~,~)XYfun.findquar(zapFcn()),...
                varargin{:});
        end

        function ct = ToHourlyCategorical(val)
            ct = categorical(val, 0:23, 'ordinal', true);
        end
    end % static methods

end %classdef

%% Callbacks

% All callbacks should set values within the same field. Leave
% the gathering of values to the SetValuesFromDialog button.

function cb_info(~,~)
    web(['file:' ZmapGlobal.Data.hodi '/help/quarry.htm'], '-browser') ;
end

author: CelsoReyes


text:
% -------------------------------------------------------------------------------------------------------------------------
function [newTargetPosition, bestScale] = tracker_eval(net_x, s_x, scoreId, z_features, x_crops, targetPosition, window, p)
%TRACKER_STEP
%   runs a forward pass of the search-region branch of the pre-trained Fully-Convolutional Siamese,
%   reusing the features of the exemplar z computed at the first frame.
%
%   Luca Bertinetto, Jack Valmadre, Joao F. Henriques, 2016
% -------------------------------------------------------------------------------------------------------------------------
    % forward pass, using the pyramid of scaled crops as a "batch"
    net_x.eval({p.id_feat_z, z_features, 'instance', x_crops});
    responseMaps = reshape(net_x.vars(scoreId).value, [p.scoreSize p.scoreSize p.numScale]);
    responseMaps = gather(responseMaps);
    responseMapsUP = single(zeros(p.scoreSize*p.responseUp, p.scoreSize*p.responseUp, p.numScale));
    % Choose the scale whose response map has the highest peak
    if p.numScale>1
        currentScaleID = ceil(p.numScale/2);
        bestScale = currentScaleID;
        bestPeak = -Inf;
        for s=1:p.numScale
            if p.responseUp > 1
                % upsample to improve accuracy
                responseMapsUP(:,:,s) = imresize(responseMaps(:,:,s), p.responseUp, 'bicubic');
            else
                responseMapsUP(:,:,s) = responseMaps(:,:,s);
            end
            thisResponse = responseMapsUP(:,:,s);
            % penalize change of scale
            if s~=currentScaleID, thisResponse = thisResponse * p.scalePenalty; end
            thisPeak = max(thisResponse(:));
            if thisPeak > bestPeak, bestPeak = thisPeak; bestScale = s; end
        end
        responseMap = responseMapsUP(:,:,bestScale);
    else
        responseMap = responseMapsUP;
        bestScale = 1;
    end
    % make the response map sum to 1
    responseMap = responseMap - min(responseMap(:));
    responseMap = responseMap / sum(responseMap(:));
    % apply windowing
    responseMap = (1-p.wInfluence)*responseMap + p.wInfluence*window;
    [r_max, c_max] = find(responseMap == max(responseMap(:)), 1);
    [r_max, c_max] = avoid_empty_position(r_max, c_max, p);
    p_corr = [r_max, c_max];
    % Convert to crop-relative coordinates to frame coordinates
    % displacement from the center in instance final representation ...
    disp_instanceFinal = p_corr - ceil(p.scoreSize*p.responseUp/2);
    % ... in instance input ...
    disp_instanceInput = disp_instanceFinal * p.totalStride / p.responseUp;
    % ... in instance original crop (in frame coordinates)
    disp_instanceFrame = disp_instanceInput * s_x / p.instanceSize;
    % position within frame in frame coordinates
    newTargetPosition = targetPosition + disp_instanceFrame;
end

function [r_max, c_max] = avoid_empty_position(r_max, c_max, params)
    if isempty(r_max)
        r_max = ceil(params.scoreSize/2);
    end
    if isempty(c_max)
        c_max = ceil(params.scoreSize/2);
    end
end

author: shenjianbing

repo: TripletTracking


text:
%------------------------------------------------------------------
%---------------Load Data and Preprocessing -----------------------
%------------------------------------------------------------------

%   Copyright 2003-2009 Dimitrios Ververidis, AIIA Lab.
%   $Revision: 0.0$ $Date: 09/01/2009$

% Do by your will:
% 1. Remove redudant features: Features with many Nans or unique
%    values
% 2. Normalize feature values to 0-1

function [Patterns, Targets] = DataLoadAndPreprocess(DatasetToUse)

if strfind(DatasetToUse,'.mat')
  DatasetToUse = DatasetToUse(1:(length(DatasetToUse)-4));
end

%--------------------- Load data ----------------------------------
cd('PatTargMatrices')
load([DatasetToUse]);
cd('..');
%------------------- End Load data  -------------------------------

%-----------------Prepare Targets and Patterns --------------------
[NPatterns, KInitialFeatures] = eval(['size(' DatasetToUse ')']);
patterns = eval([DatasetToUse '(1:NPatterns,1:KInitialFeatures)']);

if strcmp(DatasetToUse,'finalvecDES') || ...
       strcmp(DatasetToUse,'finalvecSUSAS') || ...
               strcmp(DatasetToUse,'finalvecKidsVR')
   emotions     = eval([DatasetToUse '(1:NPatterns, 114)']);
   if strcmp(DatasetToUse, 'finalvecSUSAS')
       EmotionsOfInterest = ...
       (emotions == 1 | emotions == 2 | ...
        emotions == 6 | emotions == 8); % | ...
%        emotions == 8 | emotions == 9 | ...
 %       emotions == 10 | emotions == 11);
%   1     2     3     4     5     6     7     8     9     10   11
% 'ARY' 'CAR' 'C50' 'C70' 'FST' 'LRD' 'LUD' 'NAL' 'QON' 'SOW' 'SFT'

       patterns  = patterns(EmotionsOfInterest, :);
       emotions = emotions( EmotionsOfInterest);
       emotions(emotions ==6) = emotions(emotions == 6) -3;
       emotions(emotions ==8) = emotions(emotions == 8) -4;
%        emotions(emotions ==8) = emotions(emotions ==8) -3;
%        emotions(emotions ==9) = emotions(emotions ==9) - 3;
%        emotions(emotions ==10) = emotions(emotions ==10) - 3;
%        emotions(emotions ==11) = emotions(emotions ==11) - 3;
   end

   [NPatterns,  KInitialFeatures]  = size(patterns);
   KFeatures           =     90;   % 90 Remained after Processing

   %------------------- Missing data handle --------
   for FeatureIndex = 1:KInitialFeatures
       FeatureMean(FeatureIndex)=nanmean(patterns(:,FeatureIndex));
       for PatternIndex = 1:NPatterns  ,
           if isnan(patterns(PatternIndex,FeatureIndex)),
               patterns(PatternIndex, FeatureIndex) = ...
                                         FeatureMean(FeatureIndex);
           end
       end
   end

   %--------------Exponential Normalization------------------------
   tabexp   = [13 14 30:32 37 39:40 46 53 64:66 68 ...
       71 86:104  106:113];   % Exponentially distributed features

   for FeaturesIndex=1:length(tabexp),
       lamda=1/mean(patterns(:,tabexp(FeaturesIndex)));
       for PatternsIndex=1:NPatterns
           patterns(PatternsIndex,tabexp(FeaturesIndex)) = (1 - ...
              exp(-lamda* patterns(PatternsIndex,...
              tabexp( FeaturesIndex ))))/(1 - exp(-lamda));
       end
   end

   %-----------------Linear Normalization--------------------------
   for FeaturesIndex=1:KInitialFeatures,
       a = min(patterns(:,FeaturesIndex));
       b = max(patterns(:,FeaturesIndex));
       patterns(:, FeaturesIndex)=(patterns(:, FeaturesIndex) -...
                                                          a)/(b-a);
   end

    %---------------- Remove features ----------------
    tabNans = [23:29 48 57:63];               % Feature with Nans
    tabbias = [8 33:34 41 60 67 75 82 105];   % Features with Bias
    tabnw   = sort([tabNans tabbias]);        % Useless features

    tabfin = 1:KInitialFeatures;
    for FeatureIndex = 1:length(tabnw)
        tabfin(tabnw(FeatureIndex)) = 0;
    end

    tabfin = tabfin(find(tabfin~=0));
    patterns = patterns(:,tabfin);

    Patterns = patterns(1:NPatterns,1:KFeatures);
    Targets = emotions;
elseif strcmp(DatasetToUse,'finalvecCOLONCANCER')
    for FeaturesIndex= 1:KInitialFeatures-1,
        a = min(patterns(:,FeaturesIndex));
        b = max(patterns(:,FeaturesIndex));
        patterns(:, FeaturesIndex)=(patterns(:, FeaturesIndex) -...
            a)/(b-a);
    end
    Patterns = patterns(:,1:(KInitialFeatures-1));
    Targets  = eval([DatasetToUse '(1:NPatterns, end)']);

    NewTargets = zeros(NPatterns,1);
    NewTargets(find(Targets>0)) = 1;
    NewTargets(find(Targets<0)) = 2;
    Targets = NewTargets;
    disp('1 responds to Benign, 2 responds to Malignant');
else   % Your Data is loaded here
       % Linear transformation, Normalization to [0,1]
    for FeaturesIndex= 1:KInitialFeatures-1,
        a = min(patterns(:,FeaturesIndex));
        b = max(patterns(:,FeaturesIndex));
        patterns(:, FeaturesIndex)=(patterns(:, FeaturesIndex) -...
            a)/(b-a);
    end
    Patterns = patterns(:,1:(KInitialFeatures-1));
    Targets  = eval([DatasetToUse '(1:NPatterns, end)']);
    disp(['No Preprocessing Code, however PatternTargets format'...
        'is appropriate']);
end

disp('Size of Patterns');
size(Patterns)
disp('Size of Targets');
size(Targets)
save  'Patterns.mat' Patterns
save  'Targets.mat' Targets
return

author: Sable

repo: mcbench-benchmarks


text:
function adjList = meshFaceAdjacency(vertices, edges, faces)
%MESHFACEADJACENCY Compute adjacency list of face around each face.
%
%
%   Example
%     % Create a sample 3D mesh
%     [v, e, f] = createDodecahedron;
%     adjList = meshFaceAdjacency(v, e, f);
%     figure; hold on; axis equal; view([100 40]);
%     drawMesh(v, f);
%     % draw sample face in a different color
%     drawMesh(v, f(1, :), 'faceColor', 'b');
%     % draw the neighbors of a sample face
%     drawMesh(v, f(adjList{1}, :), 'faceColor', 'g')
%
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-10-04,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

edgeFaceList = meshEdgeFaces(vertices, edges, faces);

% allocate memory for adjacency list
nFaces = max(edgeFaceList(:));
adjList = cell(1, nFaces);

% iterate over edges to populate adjacency list
for iEdge = 1:size(edgeFaceList)
    f1 = edgeFaceList(iEdge, 1);
    f2 = edgeFaceList(iEdge, 2);
    adjList{f1} = [adjList{f1} f2];
    adjList{f2} = [adjList{f2} f1];
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function fem_basis_test02 ( )

%*****************************************************************************80
%
%% FEM_BASIS_TEST02 tests FEM_BASIS_2D
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 October 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM_BASIS_TEST02\n' );
  fprintf ( 1, '  FEM_BASIS_2D evaluates an arbitrary triangular\n' );
  fprintf ( 1, '  basis function.\n' );

  i1 = 1;
  j1 = 0;
  k1 = 2;
  d = i1 + j1 + k1;
  x1 = i1 / d;
  y1 = j1 / d;
  fprintf ( 1, '\n' );
  fprintf ( 1, '   I   J   K        X           Y      L(I,J,K)(X,Y)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  %2d  %2d  %2d  %10.4f  %10.4f  %14.6g\n', ...
    i1, j1, k1, x1, y1, 1.0 );
  fprintf ( 1, '\n' );
  for j2 = 0 : d
    for i2 = 0 : d - j2
      k2 = d - i2 - j2;
      x2 = i2 / d;
      y2 = j2 / d;
      lijk = fem_basis_2d ( i1, j1, k1, x2, y2 );
      fprintf ( 1, '  %2d  %2d  %2d  %10.4f  %10.4f  %14.6g\n', ...
        i2, j2, k2, x2, y2, lijk );
    end
  end

  return
end


text:
function psf = estimate_psf(blurred_x, blurred_y, latent_x, latent_y, weight, psf_size)

    %----------------------------------------------------------------------
    % these values can be pre-computed at the beginning of each level
%     blurred_f = fft2(blurred);
%     dx_f = psf2otf([1 -1 0], size(blurred));
%     dy_f = psf2otf([1;-1;0], size(blurred));
%     blurred_xf = dx_f .* blurred_f; %% FFT (Bx)
%     blurred_yf = dy_f .* blurred_f; %% FFT (By)

    latent_xf = fft2(latent_x);
    latent_yf = fft2(latent_y);
    blurred_xf = fft2(blurred_x);
    blurred_yf = fft2(blurred_y);
    % compute b = sum_i w_i latent_i * blurred_i
    b_f = conj(latent_xf)  .* blurred_xf ...
        + conj(latent_yf)  .* blurred_yf;
    b = real(otf2psf(b_f, psf_size));

    p.m = conj(latent_xf)  .* latent_xf ...
        + conj(latent_yf)  .* latent_yf;
    %p.img_size = size(blurred);
    p.img_size = size(blurred_xf);
    p.psf_size = psf_size;
    p.lambda = weight;

    psf = ones(psf_size) / prod(psf_size);
    psf = conjgrad(psf, b, 20, 1e-5, @compute_Ax, p);

    psf(psf < max(psf(:))*0.05) = 0;
    psf = psf / sum(psf(:));
end

function y = compute_Ax(x, p)
    x_f = psf2otf(x, p.img_size);
    y = otf2psf(p.m .* x_f, p.psf_size);
    y = y + p.lambda * x;
end

author: IVRL

repo: Kernel-Modeling-Super-Resolution


text:
% This demo script runs the STRCF tracker with deep features on the
% included "Human3" video.

% Add paths
setup_paths();

%  Load video information
base_path  =  './sequences';
%video  = choose_video(base_path);
video = 'Human3';

video_path = [base_path '/' video];
[seq, gt_boxes] = load_video_info(video_path);

% Run STRCF
results = run_STRCF(seq);
%results = run_DeepSTRCF(seq);

pd_boxes = results.res;
thresholdSetOverlap = 0: 0.05 : 1;
success_num_overlap = zeros(1, numel(thresholdSetOverlap));
res = calcRectInt(gt_boxes, pd_boxes);
for t = 1: length(thresholdSetOverlap)
    success_num_overlap(1, t) = sum(res > thresholdSetOverlap(t));
end
cur_AUC = mean(success_num_overlap) / size(gt_boxes, 1);
FPS_vid = results.fps;
display([video  '---->' '   FPS:   ' num2str(FPS_vid)   '    op:   '   num2str(cur_AUC)]);

repo: STRCF


text:
chi_095_2 = 5.9915;
% chi_099_2 = 9.2103;
chi_095_3 = 7.8147;
% chi_099_3 = 11.3449;

% plot image stuff
figure(figure_all);
subplot(im_fig);
hold off;
imagesc(tracker.lastImg);
colormap gray;
hold on;
title('Thick red: low innovation inliers. Thin red: high innovation inliers. \newline Magenta: rejected by 1-point RANSAC. Blue: No match found by pyramid KLT.');
half_patch_size_when_initialized = 20;
half_patch_size_when_matching = 6;
features_info=filter.features_info;

for i=1:length(features_info)

    if (~isempty(features_info(i).h)&&~isempty(features_info(i).S))

        %         imagesc( features_info(i).h(1) - half_patch_size_when_matching,...
        %             features_info(i).h(2) - half_patch_size_when_matching, features_info(i).patch_when_matching);

        if features_info(i).low_innovation_inlier
            plotUncertainEllip2D( features_info(i).S,...
                features_info(i).h, chi_095_2, 'r', 4 )
            plot( features_info(i).h(1), features_info(i).h(2),'r+','Markersize',10);
        end
        if features_info(i).high_innovation_inlier
            plotUncertainEllip2D( features_info(i).S,...
                features_info(i).h, chi_095_2, 'r', 2 )
            plot( features_info(i).h(1), features_info(i).h(2),'r+','Markersize',10);
        end
        if (~isempty(features_info(i).z))&&(features_info(i).high_innovation_inlier==0)&&(features_info(i).low_innovation_inlier==0)
            plotUncertainEllip2D( features_info(i).S,...
                features_info(i).h, chi_095_2, 'm', 2 )
            plot( features_info(i).h(1), features_info(i).h(2),'m+','Markersize',10);
        end

        if (isempty(features_info(i).z)&&~isempty(features_info(i).h))
            plotUncertainEllip2D( features_info(i).S,...
                features_info(i).h, chi_095_2, 'b', 2 )
            plot( features_info(i).h(1), features_info(i).h(2),'b+','Markersize',10);
        end

        if (~isempty(features_info(i).z))
            plot( features_info(i).z(1), features_info(i).z(2),'g+','Markersize',10);
            plot([features_info(i).z(1),features_info(i).h(1)],[features_info(i).z(2),features_info(i).h(2)],'Color','r','LineWidth',2);
        end
    end

end

% plot predicted and measured
% which_are_predicted = predicted_measurements(:,1)>0;
% plot( predicted_measurements(which_are_predicted,1), predicted_measurements(which_are_predicted,2),'r+','Markersize',10);
% which_are_measured = measurements(:,1)>0;
% plot( measurements(which_are_measured,1), measurements(which_are_measured,2),'g+');

axes_handler = get(gcf,'CurrentAxes');
set(axes_handler,'XTick',[],'YTick',[]);

% plot 3D stuff
figure(figure_all);
subplot(near3D_fig);
hold off;

imgctr=step-firstGroupFrameId;
qs02c0=filter.camPose(1:4);
% obtain qc2c0 current camera frame c to start camera frame c0, T c in c0
qs02c=quatmult_v001(filter.camPose(1:4),filter.rvqs0(7:10),0);
qc2c0=quatmult_v001(qs02c0,qs02c,2);
Tc2c0=-quatrot_v000(qc2c0, filter.camPose(5:7),0)+...
    quatrot_v000(qs02c0,filter.rvqs0(1:3),0)+filter.camPose(5:7);
trajectory(:,imgctr) =[Tc2c0; qc2c0];
draw_camera( trajectory(:,imgctr), 'k' );
hold on;
title(sprintf('Camera motion and scene [m] in top view at frmId(1 based)%d.',step));

plot3( trajectory(1, 1:imgctr), trajectory(2, 1:imgctr),...
    trajectory(3, 1:imgctr), 'k', 'LineWidth', 2 );

% for each group frame, compute Cci2c0 and Tcj2c0
grpLeg=length(filter.groupPose);
groupTransform=zeros(7, length(filter.groupPose));
qs02c0=filter.camPose(1:4);
for apple=1:grpLeg
    qs02cj=quatmult_v001(filter.camPose(1:4),filter.groupPose(apple).pose(1:4),0);
    groupTransform(1:4,apple)=quatmult_v001(qs02c0,qs02cj,2);
    groupTransform(5:7,apple)=-quatrot_v000(groupTransform(1:4,apple), filter.camPose(5:7),0)+...
        quatrot_v000(qs02c0,filter.groupPose(apple).pose(5:7),0)+...
        filter.camPose(5:7);
end
groupIds=[filter.groupPose.grpId];
covDim=filter.groupFrameSIP-1+length(filter.groupPose)*6+length(filter.features_info);
featSIP=filter.groupFrameSIP+length(filter.groupPose)*6;
%for each feature with lostNo==0, obtain its position in c0 frame
for wasp = 1:length(features_info)
    if (features_info(wasp).lostNo==0)
        % the feature is still being tracked, note some point
        % may be lost. But before a new group comes, it is still in
        % the states
        apple=find(groupIds==filter.features_info(wasp).grpId);

        % because grpId contains gap, this instruction can be sped up by making use of bisection search
        rho=filter.features_info(wasp).invDepth;
        XYZ=quatrot_v000(groupTransform(1:4,apple), [filter.features_info(wasp).xyn; 1]/rho,0)+...
            groupTransform(5:7,apple);
        if rho-3*sqrt(filter.p_k_k(featSIP+wasp-1,featSIP+wasp-1))<0
            if ( rho>0)
                plot3(XYZ(1),XYZ(2),XYZ(3),'k^','Markersize',10)
            end
        else
            if ( rho>0)
                plot3(XYZ(1),XYZ(2),XYZ(3),'r+','Markersize',10)
            end
        end
    end
end

axes_handler = get(gcf,'CurrentAxes');
axis([-100 150 -5 6 -200 210]);
xlabel('x');
ylabel('y');
zlabel('z');
grid on;
hold off
view([0, 0]);% only x and z

repo: ekfmonoslam


text:
function test_pdfs_vbmc()
%TEST_PDFS_VBMC Test pdfs introduced in the VBMC package.

lb = [-1.1,-4.1];
ub = [3.2,-2.8];
a = [-1,-4];
b = [3,-3];
n = 1e6;

tolerr = 1e-3;      % Error tolerance on normalization constant
tolrmse = 0.05;     % Error tolerance on histogram vs pdf

%% Test multivariate uniform box distribution
pdf1 = @(x) munifboxpdf(x,a(1),b(1));
pdf2 = @(x) munifboxpdf(x,a,b);
pdf1log = @(x) exp(munifboxlogpdf(x,a(1),b(1)));
pdf2log = @(x) exp(munifboxlogpdf(x,a,b));
pdfrnd = @(n) munifboxrnd(a,b,n);
name = 'munifbox';

test_pdf1_normalization(pdf1,lb,ub,tolerr,name);
test_pdf2_normalization(pdf2,lb,ub,tolerr,name);
test_pdf1_normalization(pdf1log,lb,ub,tolerr,name);
test_pdf2_normalization(pdf2log,lb,ub,tolerr,name);
test_rnd(pdfrnd,pdf1,a,b,n,tolrmse,name);

%% Test multivariate trapezoidal distribution
u = [-0.5,-3.8];
v = [1.5,-3.4];
pdf1 = @(x) mtrapezpdf(x,a(1),u(1),v(1),b(1));
pdf2 = @(x) mtrapezpdf(x,a,u,v,b);
pdf1log = @(x) exp(mtrapezlogpdf(x,a(1),u(1),v(1),b(1)));
pdf2log = @(x) exp(mtrapezlogpdf(x,a,u,v,b));
pdfrnd = @(n) mtrapezrnd(a,u,v,b,n);
name = 'mtrapez';
test_pdf1_normalization(pdf1,lb,ub,tolerr,name);
test_pdf2_normalization(pdf2,lb,ub,tolerr,name);
test_pdf1_normalization(pdf1log,lb,ub,tolerr,name);
test_pdf2_normalization(pdf2log,lb,ub,tolerr,name);
test_rnd(pdfrnd,pdf1,a,b,n,tolrmse,name);

%% Test multivariate spline trapezoidal distribution
pdf1 = @(x) msplinetrapezpdf(x,a(1),u(1),v(1),b(1));
pdf2 = @(x) msplinetrapezpdf(x,a,u,v,b);
pdf1log = @(x) exp(msplinetrapezlogpdf(x,a(1),u(1),v(1),b(1)));
pdf2log = @(x) exp(msplinetrapezlogpdf(x,a,u,v,b));
pdfrnd = @(n) msplinetrapezrnd(a,u,v,b,n);
name = 'msplinetrapez';
test_pdf1_normalization(pdf1,lb,ub,tolerr,name);
test_pdf2_normalization(pdf2,lb,ub,tolerr,name);
test_pdf1_normalization(pdf1log,lb,ub,tolerr,name);
test_pdf2_normalization(pdf2log,lb,ub,tolerr,name);
test_rnd(pdfrnd,pdf1,a,b,n,tolrmse,name);

%% Test multivariate smoothbox distribution
sigma = [0.7,0.45];
pdf1 = @(x) msmoothboxpdf(x,a(1),b(1),sigma(1));
pdf2 = @(x) msmoothboxpdf(x,a,b,sigma);
pdf1log = @(x) exp(msmoothboxlogpdf(x,a(1),b(1),sigma(1)));
pdf2log = @(x) exp(msmoothboxlogpdf(x,a,b,sigma));
pdfrnd = @(n) msmoothboxrnd(a,b,sigma,n);
name = 'msmoothbox';
lb = [-5,-7];
ub = [5,0];
test_pdf1_normalization(pdf1,lb,ub,tolerr,name);
test_pdf2_normalization(pdf2,lb,ub,tolerr,name);
test_pdf1_normalization(pdf1log,lb,ub,tolerr,name);
test_pdf2_normalization(pdf2log,lb,ub,tolerr,name);
test_rnd(pdfrnd,pdf1,a,b,n,tolrmse,name);

%close all;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function test_pdf1_normalization(pdf1,lb,ub,tol,name)
%TEST_PDF1_NORMALIZATION Test normalization of univariate pdf.

% Check 1D integral
y = integral(@(x) pdf1(x), lb(1), ub(1), 'ArrayValued', true);
fprintf('%s: 1D integral: %.6f\n', name, y);
assert(abs(y - 1) < tol, ['Test error: univariate ' name ' does not integrate to 1.']);

end

function test_pdf2_normalization(pdf2,lb,ub,tol,name)
%TEST_PDF2_NORMALIZATION Test normalization of bivariate pdf.

% Check 2D integral
y = integral2(@(x1,x2) reshape(pdf2([x1(:),x2(:)]),size(x1)), ...
    lb(1), ub(1), lb(2), ub(2));
fprintf('%s: 2D integral: %.6f\n', name, y);
assert(abs(y - 1) < tol, ['Test error: bivariate ' name ' does not integrate to 1.']);

end

function test_rnd(pdfrnd,pdf1,a,b,n,tol,name)
%TEST_RND Test random sample generation (histogram vs pdf).

r = pdfrnd(n);
h = histogram(r(:,1),100,'BinLimits',[a(1),b(1)],'Normalization','pdf');
x = 0.5*(h.BinEdges(1:end-1) + h.BinEdges(2:end))';
y = pdf1(x)';
rmse = sqrt(sum(((y - h.Values)).^2*h.BinWidth));
fprintf('%s: histogram rmse: %.6f\n', name, rmse);
assert(rmse < tol, ['Test error: generated histogram does not match ' name ' pdf.']);

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function svgTile(s, source, result)
% Adds a feTile SVG filter
% PRELIMINARY IMPLEMENTATION (Parameters may change)
%
% svgTile(s, source, result)
% Parameters:
%   s : Array of plot object handles
%   source : Any previous defined filter result string, 'SourceGraphic',
%            or 'SourceAlpha'.
%   result : String that identifies the filter result for following filter
%            stages.
for i = 1:length(s)
    userdata = get(s(i),'UserData');
    if isfield(userdata, 'svg') && isfield(userdata.svg, 'Filter')
        next = length(userdata.svg.Filter) + 1;
    else
        next = 1;
    end
    userdata.svg.Filter(next).Subfilter.Type = 'feTile';
    userdata.svg.Filter(next).Subfilter.Source = source;
    userdata.svg.Filter(next).Subfilter.Result = result;
    set(s(i),'UserData', userdata);
end

author: marianux

repo: ecg-kit

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


text:
% Multi-Frame Analysis based on derivative of Linear Interpolation (dLI-MFA)
%
% Input
%  af             : Array of structures with matrices containing sinusoidal
%                   parameters (as the output of sin_analysis.m).
%                   Each matrix is made of:
%                      The 1st line for the frequency of each sinusoid [Hz]
%                      The 2nd line for their amplitude (linear scale)
%                   The DC has to be included (in the first column)
%  fs             : [Hz] Signal's sampling frequency
%  [extrap_dcny] : If true (default), extrapolate sinusoidal components at DC
%                  and up to Nyquist.
%                  If false, use the sinusoidal components as they are.
%
% Output
%  E              : The amplitude cepstral envelope
%
% Reference
%  [1] G. Degottex, "A Time Regularization Technique for Discrete Spectral
%      Envelopes Through Frequency Derivative", Signal Processing Letters, IEEE,
%      22(7):978-982, July 2015.
%
% Copyright (c) 2013 Foundation for Research and Technology-Hellas - Institute
%                    of Computer Science (FORTH-ICS)
%
% License
%  This file is part of libphoni. libphoni is free software: you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. libphoni is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function [E] = env_dli_mfa(af, fs, dftlen, extrap_dcny)

    if nargin<4; extrap_dcny=true; end

    % If asked, extrapolate components at DC and up to Nyquist
    if extrap_dcny; af = env_extrap_sins_dcny(af, fs); end

    % For each frame, compute the frequency derivative of the linear envelope
    F = fs*(0:dftlen/2)/dftlen;
    dA = ones(numel(af),dftlen/2);
    for fi=1:numel(af)
        A = interp1(af(fi).sins(1,:), log(af(fi).sins(2,:)), F, 'linear', 'extrap');
        dA(fi,:) = diff(A); % derivative approximation
    end

    % Smooth the frequency derivative across time
    win = hamming(size(dA,1));
    win = win./sum(win);
    dAw = dA.*repmat(win, 1, size(dA,2));
    mA = sum(dAw,1);

    % Retrieve the envelope
    E = cumsum(mA);
    E = [E(1), E];

    % Align the envelope on the harmonics of the central frame
    ci = floor((numel(af)-1)/2)+1;
    ek = interp1(F', E, af(ci).sins(1,:));
    ak = log(af(ci).sins(2,:));
    idx = find(af(ci).sins(1,:)>0 & af(ci).sins(1,:)<4000);
    E = E + mean(ak(idx)) - mean(ek(idx));

    E = exp(E);

    % Plot the final solution
    if 0
        subplot(211);
            hold off;
            plot(F(1:end-1), dA, 'r');
            hold on;
            plot(F(1:end-1), mA, 'b');
            xlim([0 fs/2]);
            xlabel('Frequency [Hz]');
        subplot(212);
            hold off;
            plot(F, mag2db(E), 'b');
            hold on;
            plot(af(ci).sins(1,:), mag2db(af(ci).sins(2,:)), 'xk');
            xlim([0 fs/2]);
            xlabel('Frequency [Hz]');
            ylabel('Amplitude [dB]');
            title(['dLI-MFA']);
        pause
%          keyboard
    end

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function batch_fuse_atlas_seg_new(pcDirName,atlasDirName,registeredDirLoc)
%
% batch_fuse_atlas_seg.m
%
% This script fuses atlas segmentations.
%
% INPUTS: pcDirName - directory containing the PC representations in CERR format.
%       : atlasDirName - directory containing the atlas of CERR files
%       : registeredDirLoc - directory to write out the registrations and
%       to generate the atlas segmented results
%
% APA, 8/14/2018

% % directory containing all files
% pcDirName = 'H:\Public\Aditya\mimExtensions\CERR_files_Sandra_contours_PC';
% pcDirName = 'H:\Public\Aditya\mimExtensions\CERR_files_Sandra_contours_CT';
% pcDirName = 'H:\Public\Aditya\mimExtensions\Atlas_Sanne\PC_cerr';
% % pcDirName = 'H:\Public\Aditya\mimExtensions\Atlas_Sanne\CT_cropped_cerr';
% pcDirName = '/lab/deasylab1/Aditya/AtlasSeg/PC_cerr';
% pcDirName = '/lab/deasylab1/Aditya/AtlasSeg/CT_cropped_cerr';
% pcDirName = '/lab/deasylab1/Aditya/AtlasSeg/PC';
% pcDirName = '/lab/deasylab1/Aditya/AtlasSeg/PengAtlas/PC_all';
% pcDirName = '/lab/deasylab1/Aditya/AtlasSeg/SanneAtlas/PC_all';
% pcDirName = '/lab/deasylab1/Aditya/AtlasSeg/hnAtlasDec2017/PC/test';
% pcDirName = '/lab/deasylab1/Aditya/AtlasSeg/ProspectiveEval/PC/test';
%
% atlasDirName = '/lab/deasylab1/Aditya/AtlasSeg/PengAtlas/PC_atlas';
% atlasDirName = '/lab/deasylab1/Aditya/AtlasSeg/hnAtlasDec2017/PC/train';
%
% % directory for writing the registered files (must have \ or / as last character)
% % registeredDir = 'H:\Public\Aditya\mimExtensions\registered_to_ROBINSON^HEATH_35487047\';
% % registeredDir = 'H:\Public\Aditya\mimExtensions\registered_to_ROBINSON^HEATH_35487047_CT\';
% % registeredDir = 'H:\Public\Aditya\mimExtensions\registered_to_MT160_PC\';
% % registeredDir = 'H:\Public\Aditya\mimExtensions\registered_to_MT160_CT\';
%
% registeredDirLoc = 'H:\Public\Aditya\mimExtensions\registered\';
% registeredDirLoc = '/lab/deasylab1/Aditya/AtlasSeg/registered_PC/';
% registeredDirLoc = '/lab/deasylab1/Aditya/AtlasSeg/registered_PC_25pts/';
% % registeredDirLoc = '/lab/deasylab1/Aditya/AtlasSeg/registered_CT_cropped/';
% registeredDirLoc = '/lab/deasylab1/Aditya/AtlasSeg/registered_CT_cropped_25_pts/';
% registeredDirLoc = '/lab/deasylab1/Aditya/AtlasSeg/registered_PC/';
% registeredDirLoc = '/lab/deasylab1/Aditya/AtlasSeg/PengAtlas/registered_PC_all_loocv/';
% registeredDirLoc = '/lab/deasylab1/Aditya/AtlasSeg/SanneAtlas/registered_PC_all_loocv/';
% registeredDirLoc = '/lab/deasylab1/Aditya/AtlasSeg/hnAtlasDec2017/PC/registered_test/';
% registeredDirLoc = '/lab/deasylab1/Aditya/AtlasSeg/ProspectiveEval/PC/registered_test/';

%distcomp.feature( 'LocalUseMpiexec', false )
%pool = parpool(15);

% atlas scans
atlasDirS = dir(atlasDirName);
atlasDirS(1:2) = [];
movScanC = fullfile(atlasDirName,{atlasDirS.name});

dirS = dir(pcDirName);
dirS(1:2) = [];

for indBase = 1:length(dirS)

    %     registeredDir = fullfile(registeredDirLoc,['registered_to_',...
    %         strtok(dirS(indBase).name,'.mat')]);
    [~,fname] = fileparts(dirS(indBase).name);
    registeredDir = fullfile(registeredDirLoc,['registered_to_',fname]);

    % mkdir(registeredDir)

    % base scan file name
    % indBase = baseScanNum; %9 for Sandra's atlas, 3 for Sanne's
    baseScan = fullfile(pcDirName,dirS(indBase).name);

    % Leave-one-out validation file list creation
    % % moving scan file names
    %indV = 1:length(dirS);
    %indV(indBase) = [];
    %movScanC = fullfile(pcDirName,{dirS(indV).name});
    % % Leave-one-out validation file list creation ends

    % Non-LOOCV validation
    dirPatS = dir(registeredDir);
    dirPatS(1:2) = [];
    movScanC = fullfile(registeredDir,{dirPatS(:).name});

    % registration callback
    strNameToWarp = 'Parotid_L_SvD';
    % registerToAtlas(baseScan,movScanC,registeredDir,strNameToWarp)
    %registerToAtlasMultipleScans(baseScan,movScanC,registeredDir,strNameToWarp)
    segStrNamC = {'Parotid_LT_BABS','Parotid_RT_BABS'};

    %%% ---------- Fuse results from multiple atlases
    numScans = 1;
    numStructs = 2;
    structNumV = 1:numScans*numStructs; % structures on all comps
    % structNumV = 1:6; % comment for structures on all comps
    planC = loadPlanC(baseScan,tempdir);
    indexS = planC{end};
    for structNum = 1:numStructs %structNumV
        % combine using the STAPLE and the GRE metric
        regDirS = dir(registeredDir);
        regDirS(1:2) = [];
        regFilesC = fullfile(registeredDir,{regDirS.name});
        %structNum = 2; % used for just one component atlas
        doseNum = 1;
        scanNum = 1;
        numAtlases = length(regFilesC);
        uniScanSiz = getUniformScanSize(planC{indexS.scan}(scanNum));
        strAllM = zeros(prod(uniScanSiz),numAtlases,'single');

        structNumOffset = 3*numStructs; % CT + 2 PC scans
        structNumOffset = numStructs; % Only CT
        parfor i = 1:numAtlases
            for scanNum = 1:numScans
                regPlanC = loadPlanC(regFilesC{i},tempdir);
                str3M = getUniformStr(structNum+numStructs*(scanNum-1)+structNumOffset,regPlanC);
                strAllM(:,i) = str3M(:);
            end
        end

        % APA Commented this out on 7/13/2018
        % Remove outliers
        agreeV = sum(strAllM > 0,2);
        agree50V = agreeV > numAtlases*numScans/2;
        matchM = bsxfun(@xor,strAllM > 0,agree50V);
        numMatchV = sum(~matchM,1);
        p50 = quantile(numMatchV,0.1);
        noOutV = numMatchV > p50;
        %p15 = quantile(numVoxelsV,0.15);
        %p85 = quantile(numVoxelsV,0.85);
        %noOutV = numVoxelsV > p15 & numVoxelsV < p85;
        %noOutV = true(1,size(strAllM,2)); % APA added 7/13/2018

        %confidenceV = linspace(0.3,0.6,4);
        confidenceV = 0.5;
        numConf = length(confidenceV);

        % Average agreement with and without outliers
        meanAgreeV = mean(strAllM > 0, 2);

        meanAgreeAfterOutRemV = mean(strAllM(:,noOutV) > 0, 2);

        % Structure Name
        structName = planC{indexS.structures}(structNum).structureName;
        newStrNam = segStrNamC{structNum};

        for iConf = 1:numConf
            confidence = confidenceV(iConf);
            isUniform = 1;
            scanNum = 1;

            % Majority Vote
            majVote3M = reshape(meanAgreeAfterOutRemV >= confidence, uniScanSiz);
            %newStrNam = [structName,'_MjV_',num2str(confidence*100)]
            planC = maskToCERRStructure(majVote3M,isUniform,scanNum,...
                newStrNam,planC);
            planC = deleteStructureSegments(length(planC{indexS.structures}),...
                0.05,planC);
            planC = smoothContour(length(planC{indexS.structures}),planC);

        end % end of staple-confidence loop

    end % end of loop over structNumv

    % Save planC
    outputFileNam = fullfile(registeredDir,dirS(indBase).name);
    save_planC(planC,[],'passed',outputFileNam);

end

function planC = smoothContour(structNum,planC)
indexS = planC{end};
for slc = 1:length(planC{indexS.structures}(structNum).contour)
    for seg = 1:length(planC{indexS.structures}(structNum).contour(slc).segments)
        ptsM = planC{indexS.structures}(structNum).contour(slc).segments(seg).points;
        if isempty(ptsM)
            continue;
        end
        numPts = size(ptsM,1);
        intrvl = ceil(numPts*0.2/10);
        pts1M = spcrv(ptsM(1:intrvl:end,1:2)',3,100)';
        pts1M(:,3) = ptsM(1,3)*pts1M(:,1).^0;
        pts1M(end+1,:) = pts1M(1,:);
        planC{indexS.structures}(structNum).contour(slc).segments(seg).points = pts1M;
    end
end
planC = getRasterSegs(planC, structNum);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function vscl = vscale( f )
%VSCALE   Estimate the vertical scale of a function.
%   VSCALE(F) estimates the vertical scale (also known as the dynamical range)
%   of a function. This is required because a CHEBTECH does not store its
%   interpolation data. If F is an array-valued CHEBTECH with K columns, then
%   the result is a row vector of length K.

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

% Get the coefficients:
c = f.coeffs;

% Check isempty:
if ( isempty( c ) )
    vscl = 0;
elseif ( size(c, 1) == 1 )
    vscl = abs(c);
else
    % Compute values:
    vals = f.coeffs2vals( c );
    % Take max:
    vscl = max(abs(vals), [], 1);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function varargout=ceil(varargin)
%CEIL (overloaded)

switch class(varargin{1})

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.

        x = varargin{1};
        dim = size(x);
        x = reshape(x,prod(dim),1);
        y = [];
        for i = 1:prod(dim)
            y = [y;yalmip('define',mfilename,extsubsref(x,i))];
        end
        y = reshape(y,dim);
        varargout{1} = y;

    case 'char' % YALMIP send 'graph' when it wants the epigraph or hypograph
        switch varargin{1}
            case {'milp','graph','exact'}
                % Description using epigraphs
                t = varargin{2};
                X = varargin{3};

                F = ([X <= t <= X + 1]) + (integer(t));

                varargout{1} = F;
                varargout{2} = struct('convexity','none','monotonicity','none','definiteness','none','model','integer');
                varargout{3} = X;

            otherwise
                error('SDPVAR/CEIL called with CHAR argument?');
        end
    otherwise
        error('Strange type on first argument in SDPVAR/CEIL');
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function y = cvx_isnonzero( x, full )
error( nargchk( 1, 2, nargin ) );
y = any( x.basis_, 1 );
if nargin < 2,
    y = all( y );
else
    y = cvx_reshape( y, x.size_ );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function b = r8col_normalize_li ( m, n, a )

%*****************************************************************************80
%
%% R8COL_NORMALIZE_LI normalizes an R8COL with the column infinity norm.
%
%  Discussion:
%
%    Each column is scaled so that the entry of maximum norm has the value 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), the array to be examined.
%
%    Output, real B(M,N), the normalize array.
%
  b(1:m,1:n) = a(1:m,1:n);

  for j = 1 : n
    c = a(1,j);
    for i = 2 : m
      if ( abs ( c ) < abs ( a(i,j) ) )
        c = a(i,j);
      end
    end
    if ( c ~= 0.0 )
      b(1:m,j) = b(1:m,j) / c;
    end
  end

  return
end


text:
function roiName = mrmNameROI(pos);
%
%   roiName = mrmNameROI(pos);
%
%Author: Wandell
%Purpose:
%   Generate a standard name for an ROI produced by the mrMesh routines.
%   By separating this routine out, we can assure that there is consistency
%   between the different ROI names for searching.
%
%   pos is the 3D position in the mrMesh window.
%
% Example:
%   pos =mrmGet(msh,'cursor');
%   roiName = mrmNameROI(pos);
%

if ieNotDefined('pos'), pos = [0,0,0]; end

roiName = sprintf('mrm-%.0f-%.0f-%.0f',pos(1),pos(2),pos(3));

return;


text:
function ground_truth = get_ground_truth( database, image_id )
    if strcmp(database,'pascal2012')
        ground_truth.object = imread(fullfile(database_root_dir(database), 'SegmentationObject', [image_id '.png']));
        ground_truth.class  = imread(fullfile(database_root_dir(database), 'SegmentationClass', [image_id '.png']));
    elseif strcmp(database,'bsds500')
        ground_truth = loadvar(fullfile(database_root_dir(database), 'ground_truth', [image_id '.mat']),'gt_seg');
    else
        error(['Unknown database: ' database]);
    end
end

author: jponttuset

repo: mcg

save_path: github-repos/MATLAB/jponttuset-mcg


text:
function dXdt = dynamics_SRB(t,Xt,Ut,Xd,U_ext,p)

%% parameters
mass = p.mass;
J = p.J;       % inertia tensor in body frame {B}
g = 9.81;

%% decompose
% X = [pc dpc vR wb pf]'
pc = reshape(Xt(1:3),[3,1]);
dpc = reshape(Xt(4:6),[3,1]);
R = reshape(Xt(7:15),[3,3]);
wb = reshape(Xt(16:18),[3,1]);
pf34 = reshape(Xt(19:30),[3,4]);
pfd34 = reshape(Xd(19:30),[3,4]);

% r
r34 = pf34 - repmat(pc,[1,4]);

% GRF
f34 = reshape(Ut,[3,4]);

%% dynamics
ddpc = 1/mass * (sum(f34,2) + U_ext) + [0;0;-g];
dR = R * hatMap(wb);

tau_s = zeros(3,1);       % body torque expressed in {S}
for ii = 1:4
    tau_s = tau_s + hatMap(r34(:,ii)) * f34(:,ii);
end
tau_ext = hatMap(R * p.p_ext) * U_ext;
tau_tot = sum(tau_s,2) + tau_ext;
dwb = J \ (R' * tau_tot - hatMap(wb) * J * wb);

dpf = p.Kp_sw * (pfd34(:) - pf34(:));

dXdt = [dpc;ddpc;dR(:);dwb;dpf];

end

author: YanranDing

repo: RF-MPC

save_path: github-repos/MATLAB/YanranDing-RF-MPC


text:
function calc_McGrid_plot(mCatalog, fBinning)
% function calc_McGrid_plot(mCatalog, fBinning);
% --------------------------------------------
% Determine Mc using maximum likelihood score
% Fitting non-cumulative frequency magnitude distribution above and below Mc:
% below Mc with an exponential function: M = c*exp(d*M)
% above: Gutenberg-Richter law
% Normalized version !!!!
% Incoming variables:
% mCatalog   : EQ catalog
% fBinning   : Binning interval, usually 0.1
%
% Outgoing variables:
%
% J. Woessner: woessner@seismo.ifg.ethz.ch
% last update: 15.11.02

% Initialize
vProbability = [];
vMc = [];
vABValue =[];
vX_res = [];
vNCumTmp = [];
mDataPred = [];

% Determine exact time period
fPeriod1 = max(mCatalog(:,3)) - min(mCatalog(:,3));

% Determine max. and min. magnitude
fMinMag = floor(min(mCatalog(:,6)));
if fMinMag > 0
  fMinMag = 0;
end
fMaxMag = ceil(10 * max(mCatalog(:,6))) / 10;

% Calculate FMD for original catalog
%  [vFMD, vNonCFMD] = calc_FMD(mCatalog);
% % Time normalization
% vFMD(2,:) = ceil(vFMD(2,:)./fPeriod1);
% vNonCFMD(2,:) = ceil(vNonCFMD(2,:)./fPeriod1);
%  vNonCFMD = fliplr(vNonCFMD);
for fMc = 0.8:0.1:3.5
    % Calculate FMD for original catalog
    [vFMD, vNonCFMD] = calc_FMD(mCatalog);
    vNonCFMD = fliplr(vNonCFMD);
    [nIndexLo, fMagHi, vSel, vMagnitudes] = fMagToFitBValue(mCatalog, vFMD, fMc);
    [fMeanMag, fBValue, fStdDev, fAValue] =  calc_bmemag(mCatalog(vSel,:), fBinning);
    % Normalize FMDs
    vFMD(2,:) = ceil(vFMD(2,:)./fPeriod1);
    vNonCFMD(2,:) = ceil(vNonCFMD(2,:)./fPeriod1);
    % Compute quantity of earthquakes by power law
    vMstep = [fMinMag:0.1:fMaxMag];
    vNCum = 10.^(fAValue-fBValue.*vMstep); % Cumulative number
    % Compute non-cumulative numbers vN
    fNCumTmp = 10^(fAValue-fBValue*(fMaxMag+0.1));
    vNCumTmp  = [vNCum fNCumTmp];
    vN = abs(diff(vNCumTmp));
    % Normalizea-value
    fAValue = fAValue./fPeriod1;
    % Normlize vN
    vN = vN./fPeriod1;
    % Data selection above and below Mc
    mData = [vN' vNonCFMD']
    vSel = (mData(:,2) >= fMc);
    mDataTest = mData(~vSel,:)
    mDataTmp = mData(vSel,:)
    % Curve fitting: Non cumulative part below Mc
    options = optimset;
    options = optimset('Display','iter','Tolfun',1e-6,'TolX',0.0001);
    [vX, resnorm, resid, exitflag, output, lambda, jacobian]=lsqcurvefit(@calc_expdecay2,[0 1], mDataTest(:,2), mDataTest(:,3));
    mDataTest(:,1) = vX(1).*exp(vX(2).*mDataTest(:,2))-1;

     %% Set predicted data together
    mDataPred = [mDataTest; mDataTmp]
    vProb_ = calc_log10poisspdf(vNonCFMD(2,:)', ceil(mDataPred(:,1)));

    % Sum the probabilities
    fProbability = (-1) * sum(vProb_);
    vProbability = [vProbability; fProbability];
    vMc = [vMc; fMc];
    vABValue = [vABValue; fAValue fBValue];

    %% Confidence interval
    [vPred,delta] = nlpredci(@calc_expdecay2,mDataTest(:,2),vX, resid, jacobian);

    figure_w_normalized_uicontrolunits(200)
    subplot(3,1,1)
    %plot(vFactor, exp(vFactor)-1);
    plot(mDataTest(:,2), mDataTest(:,1),'-r', mDataTest(:,2), mDataTest(:,3), '*')
    hold on;
    plot(vNonCFMD(1,:)', vNonCFMD(2,:),'g^')
    %figure_w_normalized_uicontrolunits(300)
    subplot(3,1,2)
    semilogy(vNonCFMD(1,:)', vNonCFMD(2,:)', '^', vNonCFMD(1,:)', vN, '*', vNonCFMD(1,:)',mDataPred(:,1),'o')
    subplot(3,1,3)
    %figure_w_normalized_uicontrolunits(310)
    plot(mDataTest(:,2),vPred)
    hold on
    plot(mDataTest(:,2), mDataTest(:,3),'+r')
    plot(mDataTest(:,2),vPred + delta,'--g')
    plot(mDataTest(:,2),vPred - delta,'--g')
    hold off;
    if (fProbability == min(vProbability))
        vPredBest = vPred;
        vDeltaBest = delta;
        mDat = mDataTest;
        vNBest = vN;
        fMcBest = fMc;
        mDatPredBest = mDataPred;
        vX_res = [vX resnorm exitflag];
    end
    % Clear variables
    vNCumTmp = [];
    mModelDat = [];
    vNCum = [];
    vSel = [];
    mDataTest = [];
    mDataPred = [];
end
figure_w_normalized_uicontrolunits(400)
plot(vMc, vProbability,'*');
figure_w_normalized_uicontrolunits(410)
plot(mDat(:,2),vPredBest)
hold on;
plot(mDat(:,2), mDat(:,3),'+r')
plot(mDat(:,2),vPredBest + vDeltaBest,'--g')
plot(mDat(:,2),vPredBest - vDeltaBest,'--g')
hold off;
figure_w_normalized_uicontrolunits(420)
semilogy(vNonCFMD(1,:)', vNonCFMD(2,:)', '^', vNonCFMD(1,:)', vNBest, '*', vNonCFMD(1,:)',mDatPredBest(:,1),'o')
sTitlestr = ['N = ' num2str(vX(1)) ' * exp( ' num2str(vX(2)) ' * M) at ' num2str(fMcBest)];
title(sTitlestr)


text:
%DECOMPOSEESSENTIALMAT  Decompose an essential matrix to possible rotations and translation
%
%     S = cv.decomposeEssentialMat(E)
%
% ## Input
% * __E__ The input essential matrix, 3x3.
%
% ## Output
% * __S__ Decomposed `E`. A scalar struct with the following fields:
%   * __R1__ One possible rotation matrix, 3x3.
%   * __R2__ Another possible rotation matrix, 3x3.
%   * __t__ One possible translation, 3x1.
%
% This function decompose an essential matrix `E` using SVD decomposition
% [HartleyZ00]. Generally 4 possible poses exists for a given `E`. They are
% `[R1,t]`, `[R1,-t]`, `[R2,t]`, `[R2,-t]`. By decomposing `E`, you can only
% get the direction of the translation, so the function returns unit `t`.
%
% ## References
% [HartleyZ00]:
% > Richard Hartley and Andrew Zisserman. "Multiple view geometry in computer
% > vision". Cambridge university press, 2003.
%
% See also: cv.findEssentialMat, cv.SVD
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
classdef LevelSetSphere < LevelSetSphereNdim

    methods (Access = public)

        function obj = LevelSetSphere(cParams)
            obj.fracRadius = cParams.fracRadius;
            obj.compute(cParams);
        end
    end

    methods (Access = protected)
        function computeLevelSetValue(obj)
            ls = obj.dist - 1;
            obj.levelSet = ls;
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [X,stat]=fix_state(x_float,x_fb)

nx=size(x_float,1);
X=zeros(nx,1); stat=1;

Cnb = att2Cnb(x_float(1:3));
Cnb = (eye(3)+askew(x_fb(1:3)))*Cnb;
if ~isreal(Cnb)
    stat=0; return;
end

att = Cnb2att(Cnb);
vel = x_float(4:6)-x_fb(4:6);
pos = x_float(7:9)-x_fb(7:9);
bg  = x_float(10:12)+x_fb(10:12);
ba  = x_float(13:15)+x_fb(13:15);

X(1:15,1)  = [att;vel;pos;bg;ba];
X(16:nx,1) = x_float(16:end)+x_fb(16:end);

return

repo: GINav


text:
function [Q] = spm_P_clusterFDR(k,df,STAT,R,n,ui,Ps)
% Return the corrected FDR q-value
% FORMAT [Q] = spm_P_clusterFDR(k,df,STAT,R,n,ui,Ps)
%
% k        - extent {RESELS}
% df       - [df{interest} df{residuals}]
% STAT     - Statistical field
%            'Z' - Gaussian field
%            'T' - T - field
%            'X' - Chi squared field
%            'F' - F - field
% R        - RESEL Count {defining search volume}
% n        - Conjunction number
% ui       - feature-inducing threshold
% Ps       - Vector of sorted (ascending) p-values

% Q        - FDR q-value
%__________________________________________________________________________
%
% References
%
% J.R. Chumbley and K.J. Friston, "False discovery rate revisited: FDR and
% topological inference using Gaussian random fields". NeuroImage,
% 44(1):62-70, 2009.
%
% J.R. Chumbley, K.J. Worsley, G. Flandin and K.J. Friston, "Topological
% FDR for NeuroImaging". Under revision.
%__________________________________________________________________________
% Copyright (C) 2009 Wellcome Trust Centre for Neuroimaging

% Justin Chumbley & Guillaume Flandin
% $Id: spm_P_clusterFDR.m 2764 2009-02-19 15:30:03Z guillaume $

% Compute uncorrected p-values based on k using Random Field Theory
%--------------------------------------------------------------------------
[P, Z] = spm_P_RF(1, k, ui, df, STAT, R, n);

% q value using the  Benjamini & Hochberch False Discovery Rate procedure
%--------------------------------------------------------------------------
Q = spm_P_FDR(Z, df, 'P', n, Ps);

author: spm


text:
function graph_gmm(X,mi,sig,c,coefs,ft)
%
% graph_gmm(X,mi,sig,c,<coefs,ft>)
%
% plots the distribution of coefficients

DEBUG=0;
PRINT=0;
[L,T]=size(X);

if (nargin<5), coefs=1:L; end
if (nargin<6), ft=0; end

LL=length(coefs);

li=fix(sqrt(LL));
co=ceil(LL/li);

figure(1);
clf;

for ll=1:LL
l=coefs(ll);
  xm=min(X(l,:));
xM=max(X(l,:));
x=(-ft*(xM-xm)+xm):((ft+1)*(xM-xm)./100):(xM+ft*(xM-xm));

subplot(li,co,ll);

GMMImpl.histn(X(l,:),300);
hold on;

if DEBUG size(x),end

[laux,lmulti]=GMMImpl.lmultigauss(x,mi(l,:),sig(l,:),c);
aux=exp(laux);
multi=exp(lmulti);

if DEBUG size(x),size(multi'),pause,end

hp=plot(x,multi','r','Linewidth',1);
%xlim([ -xM xM ]);

ha=get(gca,'Children');
%it seem that the bars are children number 4

set(ha(2),'FaceColor',[ 0.8 0.8 0.8 ]);
set(ha(2),'EdgeColor',[ 0.8 0.8 0.8 ]);%*
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone