text:
function out = getfisx(fis, arg1, arg2, arg3, arg4, arg5)
%GETFISX Get fuzzy inference system properties.
%
%   See getfis for syntax and explanation.
%
%   It is almost clone of MATALAB's getfis
%   with some modifications, so it's compatible
%   with extendent fuzzy rule structure.
%
%   Compare also lines 189-224 of this function
%   with lines 209-239 of the original getfis.

%   Per Konstantin A. Sidelnikov, 2009.

if isfield(fis, 'input')
    numInputs=length(fis.input);
else
    numInputs=0;
end

if isfield(fis, 'output')
    numOutputs=length(fis.output);
else
    numOutputs=0;
end

switch nargin
    case 1,
        % ===============================================
        % Handle generic inquiries related to the whole fis
        % ===============================================

        fprintf('      Name      = %s\n',fis.name);
        fprintf('      Type      = %s\n',fis.type);
        fprintf('      NumInputs = %s\n',num2str(numInputs));
        fprintf('      InLabels  = \n');
        if numInputs~=0,
            for i=1:length(fis.input)
                fprintf('            %s\n',fis.input(i).name);
            end
        end
        fprintf('      NumOutputs = %s\n',num2str(numOutputs));
        fprintf('      OutLabels = \n');
        if numOutputs~=0,
            for i=1:length(fis.output)
                fprintf('            %s\n',fis.output(i).name);
            end
        end
        fprintf('      NumRules = %s\n',num2str(length(fis.rule)));
        fprintf('      AndMethod = %s\n',fis.andMethod);
        fprintf('      OrMethod = %s\n',fis.orMethod);
        fprintf('      ImpMethod = %s\n',fis.impMethod);
        fprintf('      AggMethod = %s\n',fis.aggMethod);
        fprintf('      DefuzzMethod = %s\n',fis.defuzzMethod);
        out=fis.name;

    case 2,
        % ===============================================
        propName=lower(arg1);
        switch propName
            case 'name'
                out=fis.name;
            case 'type';
                out=fis.type;
            case 'numinputs'
                out=numInputs;
            case 'numoutputs'
                out=numOutputs;
            case 'numinputmfs'
                numInputMFs=[];
                for i=1:length(fis.input)
                    numInputMFs(i)=length(fis.input(i).mf);
                end
                out=numInputMFs;
            case 'numoutputmfs'
                numOutputMFs=[];
                for i=1:length(fis.output)
                    numOutputMFs(i)=length(fis.output(i).mf);
                end
                out=numOutputMFs;

            case 'numrules'
                out=length(fis.rule);
            case 'andmethod'
                out=fis.andMethod;
            case 'ormethod'
                out=fis.orMethod;
            case 'impmethod'
                out=fis.impMethod;
            case 'aggmethod'
                out=fis.aggMethod;
            case 'defuzzmethod'
                out=fis.defuzzMethod;

            case 'inlabels'
                out=[];
                for i=1:numInputs
                    out=strvcat(out, fis.input(i).name);
                end

            case 'outlabels'
                out=[];
                for i=1:numOutputs
                    out=strvcat(out, fis.output(i).name);
                end

            case 'inrange'
                for i=1:numInputs
                    out(i, 1:2)=fis.input(i).range;
                end

            case 'outrange'
                for i=1:numOutputs
                    out(i,1:2)=fis.output(i).range;
                end

            case 'inmfs'
                for i=1:numInputs
                    out(i)=length(fis.input(i).mf);
                end

            case 'outmfs'
                for i=1:numOutputs
                    out(i)=length(fis.output(i).mf);
                end

            case 'inmflabels'
                out=[];
                for i=1:numInputs
                    for j=1:length(fis.input(i).mf)
                        out=strvcat(out, fis.input(i).mf(j).name);
                    end
                end

            case 'outmflabels'
                out=[];
                for i=1:numOutputs
                    for j=1:length(fis.output(i).mf)
                        out=strvcat(out, fis.output(i).mf(j).name);
                    end
                end

            case 'inmftypes'
                out=[];
                for i=1:numInputs
                    for j=1:length(fis.input(i).mf)
                        out=strvcat(out, fis.input(i).mf(j).type);
                    end
                end

            case 'outmftypes'
                out=[];
                for i=1:numOutputs
                    for j=1:length(fis.output(i).mf)
                        out=strvcat(out, fis.output(i).mf(j).type);
                    end
                end

            case 'inmfparams'
                numInputMFs=[];
                for i=1:length(fis.input)
                    numInputMFs(i)=length(fis.input(i).mf);
                end
                totalInputMFs=sum(numInputMFs);
                k=1;
                out=zeros(totalInputMFs, 4);
                for i=1:numInputs
                    for j=1:length(fis.input(i).mf)
                        temp=fis.input(i).mf(j).params;
                        out(k,1:length(temp))=temp;
                        k=k+1;
                    end
                end

            case 'outmfparams'
                numOutputMFs=[];
                for i=1:length(fis.output)
                    numOutputMFs(i)=length(fis.output(i).mf);
                end
                totalOutputMFs=sum(numOutputMFs);
                k=1;
                out=zeros(totalOutputMFs, 4);
                for i=1:numOutputs
                    for j=1:length(fis.output(i).mf)
                        temp=fis.output(i).mf(j).params;
                        out(k, 1:length(temp))=temp;
                        k=k+1;
                    end
                end

            case 'rulelist'
                out = {};

                if ~isempty(fis.rule)
                    ifempty = @(x) isempty(x);

                    in_rules = {fis.rule.antecedent}';
                    i = find(cellfun(ifempty, in_rules), 1, 'first');
                    if ~isempty(i)
                        error('FuzzyLogic:ruleError', ...
                            'antecedent of rule %d is empty.', i);
                    end

                    out_rules = {fis.rule.consequent}';
                    i = find(cellfun(ifempty, out_rules), 1, 'first');
                    if ~isempty(i)
                        error('FuzzyLogic:ruleError', ...
                            'consequent of rule %d is empty.', i);
                    end

                    weights = {fis.rule.weight}';
                    i = find(cellfun(ifempty, weights), 1, 'first');
                    if ~isempty(i)
                        error('FuzzyLogic:ruleError', ...
                            'weight of rule %d is empty.', i);
                    end

                    connections = {fis.rule.connection}';
                    i = find(cellfun(ifempty, connections), 1, 'first');
                    if ~isempty(i)
                        error('FuzzyLogic:ruleError', ...
                            'connection of rule %d is empty.', i);
                    end

                    out = [in_rules, out_rules, weights, connections];
                end

            case 'inputs'
                fprintf('      Name =       %s\n',fis.name);
                out.Name = fis.name;
                fprintf('      NumInputs =  %s\n',num2str(numInputs));
                out.NumInputs = numInputs;
                fprintf('      InLabels  = \n');
                if numInputs~=0,
                    for i=1:length(fis.input)
                        fprintf('            %s\n',fis.input(i).name);
                        out = setfield(out, ['input' num2str(i)], fis.input(i).name);
                    end
                end

            case 'outputs'
                fprintf('      Name =       %s\n',fis.name);
                out.Name = fis.name;
                fprintf('      NumOutputs = %s\n',num2str(numOutputs));
                out.NumOutputs = numOutputs;
                fprintf('      OutLabels = \n');
                if numOutputs~=0,
                    for i=1:length(fis.output)
                        fprintf('            %s\n',fis.output(i).name);
                        out = setfield(out, ['output' num2str(i)], fis.output(i).name);
                    end
                end

            otherwise
                error('FuzzyLogic:FISPropertyError', ...
                    sprintf('There is no FIS system property called ''%s''', propName));
        end

    case 3,
        % ===============================================
        % Handle generic inquiries related to VARIABLES
        % ===============================================
        if strcmp(arg1,'input') | strcmp(arg1,'output'),
            varType=lower(arg1);
            varIndex=arg2;

            numMFs=getfis(fis,varType,varIndex,'NumMFs');
            out.Name = getfis(fis,varType,varIndex,'Name');
            fprintf('      Name =     %s\n',out.Name);
            fprintf('      NumMFs =   %s\n',num2str(numMFs));
            out.NumMFs = numMFs;
            fprintf('      MFLabels = \n');
            if numMFs~=0,
                mfLabels=getfis(fis,varType,varIndex,'MFLabels');
                for n=1:numMFs,
                    fprintf('            %s\n',mfLabels(n,:));
                    out = setfield(out, ['mf' num2str(n)], deblank(mfLabels(n,:)));
                end
            end
            range=getfis(fis,varType,varIndex,'Range');
            fprintf('      Range =    %s\n',mat2str(range));
            out.range = range;

        end

    case 4,
        % ===============================================
        % Handle specific inquiries related to VARIABLES
        % ===============================================
        varType=lower(arg1);
        varIndex=arg2;
        varProp=lower(arg3);
        switch varType
            case 'input',
                if varIndex>numInputs,
                    error('FuzzyLogic:parameterMismatch', ...
                        sprintf('%i is not a valid input index.', varIndex));
                end

                switch varProp
                    case 'name'
                        out=fis.input(varIndex).name;
                    case 'range'
                        out=fis.input(varIndex).range;
                    case 'nummfs'
                        out=length(fis.input(varIndex).mf);
                    case 'mflabels'
                        numMFs=length(fis.input(varIndex).mf);
                        MFList=[];
                        for n=1:numMFs,
                            MFList=strvcat(MFList,fis.input(varIndex).mf(n).name);
                        end
                        out=MFList;
                    otherwise
                        error('FuzzyLogic:FISVariablePropertiesError', ...
                            sprintf(['Invalid variable properties : ''%s'' \n' ...
                            'Valid entries are: \n'...
                            '\tname \n' ...
                            '\trange \n' ...
                            '\tnummfs \n'...
                            '\tmflabels '], varProp));
                end

            case 'output',
                if varIndex>numOutputs,
                    error('FuzzyLogic:parameterMismatch', ...
                        sprintf('%i is not a valid output index.', varIndex));
                end

                switch varProp
                    case 'name'
                        out=fis.output(varIndex).name;
                    case 'range',
                        out=fis.output(varIndex).range;
                    case 'nummfs',
                        out=length(fis.output(varIndex).mf);
                    case 'mflabels',
                        numMFs=length(fis.output(varIndex).mf);
                        MFList=[];
                        for n=1:numMFs,
                            MFList=strvcat(MFList,fis.output(varIndex).mf(n).name);
                        end
                        out=MFList;
                    otherwise
                        error('FuzzyLogic:FISVariablePropertiesError', ...
                            sprintf(['Invalid variable property : ''%s'' \n' ...
                            'Valid properties are: \n'...
                            '\tname \n' ...
                            '\trange \n' ...
                            '\tnummfs \n'...
                            '\tmflabels '], varProp));
                end

            otherwise
                        error('FuzzyLogic:FISVariableError', ...
                            'Variable type must be either "input" or "output"');

        end

    case 5,
        % ===============================================
        % Handle generic inquiries related to MEMBERSHIP FUNCTIONS
        % ===============================================
        if strcmp(arg1,'input') | strcmp(arg1,'output'),
            varType=lower(arg1);
            varIndex=arg2;
            MFIndex=arg4;

            MFLabels=getfis(fis,varType,varIndex,'MFLabels');
            out.Name = getfis(fis,varType,varIndex,'MF',MFIndex,'Name');
            fprintf('      Name = %s\n',out.Name);
            out.Type = getfis(fis,varType,varIndex,'MF',MFIndex,'Type');
            fprintf('      Type = %s\n',out.Type);
            params=getfis(fis,varType,varIndex,'MF',MFIndex,'Params');
            out.params = params;
            fprintf('      Params = %s\n',mat2str(params))
        end

    case 6,
        % ===============================================
        % Handle specific inquiries related to MEMBERSHIP FUNCTIONS
        % ===============================================
        varType=lower(arg1);
        varIndex=arg2;
        MFIndex=arg4;
        MFProp=lower(arg5);

        switch varType
            case 'input'
                if varIndex>numInputs,
                    errStr=['There are only ',int2str(length(fis.input)), ...
                        ' input variables'];
                    error('FuzzyLogic:FISVariablesError', ...
                        errStr)
                end

                if MFIndex>length(fis.input(varIndex).mf),
                    errStr=['There are only ',int2str(length(fis.input(varIndex).mf)), ...
                        ' MFs associated with that variable'];
                    error('FuzzyLogic:FISVariablesError', ...
                        errStr)
                end

                switch MFProp
                    case 'name'
                        out=fis.input(varIndex).mf(MFIndex).name;
                    case 'type'
                        out=fis.input(varIndex).mf(MFIndex).type;
                    case 'params'
                        out=fis.input(varIndex).mf(MFIndex).params;
                    otherwise
                    error('FuzzyLogic:InvalidMembershipFunctionPropertyError', ...
                        sprintf(['Invalid Membership Function property : ''%s'' \n' ...
                            'Valid properties are : \n' ...
                            '\tname \n' ...
                            '\ttype \n' ...
                            '\tparams \n'], MFProp));
                end

            case 'output'
                if varIndex>numOutputs,
                    errStr=['There are only ',int2str(length(fis.output)), ...
                        ' output variables'];
                    error('FuzzyLogic:InvalidMembershipFunctionPropertyError', errStr)
                end

                if MFIndex>length(fis.output(varIndex).mf),
                    errStr=['There are only ',int2str(length(fis.output(varIndex).mf)), ...
                        ' MFs associated with that variable'];
                    error('FuzzyLogic:InvalidMembershipFunctionPropertyError', errStr)
                end

                switch MFProp
                    case 'name'
                        out=fis.output(varIndex).mf(MFIndex).name;
                    case 'type'
                        out=fis.output(varIndex).mf(MFIndex).type;
                    case 'params'
                        out=fis.output(varIndex).mf(MFIndex).params;
                    otherwise
                        error('FuzzyLogic:InvalidMembershipFunctionPropertyError', ...
                            sprintf(['Invalid Membership Function property : ''%s'' \n' ...
                            'Valid properties are : \n' ...
                            '\tname \n' ...
                            '\ttype \n' ...
                            '\tparams \n'], MFProp));
                end

            otherwise
                error('FuzzyLogic:FISVariableError', ...
                    'Variable type must be either "input" or "output"');

        end

end

author: Sable

repo: mcbench-benchmarks


text:
classdef CoreTestException < MException
    % CoreTestException. Part of the MTest unit test framework
    %
    % This class is used by the CoreMockReporting class when generating exceptions
    % during testing. This allows tests to verify that an exception has occurred
    % as expected, and allows easy distinguishing of exceptions generated by the
    % framework compared to other errors.
    %
    %     Licence
    %     -------
    %     Part of CoreMat. https://github.com/tomdoel/coremat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
    end

    methods
        function obj = CoreTestException
            obj = obj@MException('CoreTestException:FakeException', 'Fake exception');
        end
    end

end

author: tomdoel

repo: pulmonarytoolkit


text:
function [B] = KB2B(KB)
% Convert computery things from kilobytes to bytes.
% Chad A. Greene 2012
B = KB*1024;


text:
classdef BernoulliD
%%BERNOULLID Functions to handle the Bernoulli distribution. This is just a
%            binary distribution having a certain probability p of being 1.
%Implemented methods are: mean, var, PMF, CDF, rand, momentGenFun, entropy
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

methods(Static)
function val=mean(p)
%%MEAN Obtain the mean of the Bernoulli distribution.
%
%INPUTS: p The probability of success of a Bernoulli trial.
%
%OUTPUTS: val The mean of the Bernoulli distribution.
%
%As noted on the page opposite the inside front cover of [1], the mean is
%just p.
%
%EXAMPLE:
%We show that the mean matches the sample mean.
% p=0.75;
% meanVal=BernoulliD.mean(p)
% numRuns=1e6;
% meanSamp=mean(BernoulliD.rand([numRuns,1],p))
%One will see that both values are about 0.75.
%
%REFERENCES:
%[1] A. Papoulis and S. U. Pillai, Probability, Random Variables and
%    Stochastic Processes, 4th ed. Boston: McGraw Hill, 2002.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=p;
end

function val=var(p)
%%VAR Obtain the variance of the Bernoulli distribution.
%
%INPUTS: p The probability of success of a Bernoulli trial.
%
%OUTPUTS: val The variance of the Bernoulli distribution.
%
%As noted on the page opposite the inside front cover of [1], the mean is
%just p*(1-p).
%
%EXAMPLE:
%We show that the mean matches the sample mean.
% p=0.75;
% varVal=BernoulliD.var(p)
% numRuns=1e6;
% varSamp=var(BernoulliD.rand([numRuns,1],p))
%One will see that both values are about 0.1875.
%
%REFERENCES:
%[1] A. Papoulis and S. U. Pillai, Probability, Random Variables and
%    Stochastic Processes, 4th ed. Boston: McGraw Hill, 2002.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

   val=p*(1-p);
end

function val=PMF(x,p)
%%PMF Evaluate the Bernoulli probability mass function (PMF) at given
%     points.
%
%INPUTS: x The point(s) at which the Bernoulli PMF is to be evaluated. For
%          nonzero PMF values, x=0 or x=1.
%        p The probability of success of a Bernoulli trial.
%
%OUTPUTS: val The value(s) of the Bernoulli PMF.
%
%The Bernoulli distribution is given on the page opposite the inside front
%cover of [1].
%
%REFERENCES:
%[1] A. Papoulis and S. U. Pillai, Probability, Random Variables and
%    Stochastic Processes, 4th ed. Boston: McGraw Hill, 2002.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=zeros(size(x));
    val(x==0)=1-p;
    val(x==1)=p;
end

function val=CDF(x,p)
%%CDF Evaluate the cumulative distribution function of the Bernoulli
%     distribution at desired points.
%
%INPUTS: x The point(s) at which the Bernoulli PMF is to be evaluated. For
%          nonzero PMF values, x=0 or x=1.
%        p The probability of success of a Bernoulli trial.
%
%OUTPUTS: val The CDF of the Bernoulli distribution evaluated at the
%             desired point(s).
%
%The Bernoulli distribution is given on the page opposite the inside front
%cover of [1].
%
%REFERENCES:
%[1] A. Papoulis and S. U. Pillai, Probability, Random Variables and
%    Stochastic Processes, 4th ed. Boston: McGraw Hill, 2002.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=zeros(size(x));
    sel=(0<=x|x<1);
    val(sel)=1-p;
    sel=(x>=1);
    val(sel)=p;
end

function x=rand(N,p)
%%RAND Generate Bernoulli random variables with a given mean.
%
%INPUTS: N If N is a scalar, then rand returns an NXN matrix of random
%          variables. If N=[M,N1] is a two-element row vector, then rand
%          returns an MXN1 matrix of random variables.
%        p The probability of success of a Bernoulli trial.
%
%OUTPUTS: vals A matrix whose dimensions are determined by N of the
%              generated Bernoulli random variables.
%
%The Bernoulli random variables are generated simply using the comparison
%mentioned in Chapter 4 of [1]. [1].
%
%REFERENCES:
%[1] S. M. Ross, Simulation, 4th ed. Amsterdam: Academic Press, 2006.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(isscalar(N))
        dims=[N, N];
    else
        dims=N;
    end

    x=rand(dims)<=p;
end

function momentVal=momentGenFun(p,numDerivs,t)
%%MOMENTGENFUN Evaluate the moment generating function (or one of its
%              derivatives) of the Bernoulli distribution. Taking the kth
%              derivative of the moment generating function and evaluating
%              it at t=0 provides the kth noncentral moment of the
%              distribution.
%
%INPUTS: p The probability of success of a Bernoulli trial.
% numDerivs The number of derivatives to take with respect to the
%          argument of the moment generating function. numDerivs>=0.
%        t The matrix of points where the moment generating function
%          should be evaluated. If this parameter is omitted or an
%          empty matrix is passed, the default value of 0 is used.
%
%OUTPUTS: momentVal A numPointsX1 vector of the values of the derivatives
%                   of the moment generating function given at the points
%                   in t or at t=0 if t is omitted.
%
%The moment generating function of a random variable is defined to be
%E(exp(t*x)) where E is the expected value operator, x is the random
%variable and t is a real parameter. The moment generating function for the
%Bernoulli distribution is just (1-p)+p*exp(t). Derivatives are
%straightforward.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<3||isempty(t))
        t=0;
    end

    momentVal=exp(t)*p;
    sel=numDerivs==0;
    momentVal(sel)=momentVal(sel)+1-p;
end

function entropyVal=entropy(p)
%%ENTROPY Obtain the Shannon entropy of the Bernoulli distribution given in
%         nats. The Shannon entropy of a discrete distribution is
%         entropy=-sum_x Pr(x)*log(Pr(x)) where the sum is over all
%         discrete values of x. Units of nats mean that the natural
%         logarithm is used in the definition.
%
%INPUTS: p The probability of success of a Bernoulli trial.
%
%OUTPUTS: entropyVal The value of the Shannon entropy in nats.
%
%Shannon entropy is defined in Chapter 2 of [1].
%
%REFERENCES:
%[1] T. M. Cover and J. A. Thomas, Elements of Information Theory, 2nd ed.
%    Hoboken, NJ: Wiley-Interscience, 2006.
%
%April 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(p==0||p==1)
        entropyVal=0;
        return;
    end

    q=1-p;

    entropyVal=-q*log(q)-p*log(p);
end
end
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 [varargout]=image_numeric(varargin)

% function image_numeric(M)
% ------------------------------------------------------------------------
% Plots the intensities of an image (rounded to fit 4 digits) in the image
% at the pixel coordinates.
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2008/08/28: Created
%
%
% ------------------------------------------------------------------------

%%

switch nargin
    case 1
        M=varargin{1};
        hf=gcf;
        numDigits=2;
        fontSize=10;
        textColor='k';
        interpreter='tex';
    case 2
        M=varargin{1};
        hf=varargin{2};
        numDigits=2;
        fontSize=10;
        textColor='k';
        interpreter='tex';
    case 3
        M=varargin{1};
        hf=varargin{2};
        numDigits=varargin{3};
        fontSize=10;
        textColor='k';
        interpreter='tex';
    case 4
        M=varargin{1};
        hf=varargin{2};
        numDigits=varargin{3};
        fontSize=varargin{4};
        textColor='k';
        interpreter='tex';
    case 5
        M=varargin{1};
        hf=varargin{2};
        numDigits=varargin{3};
        fontSize=varargin{4};
        textColor=varargin{5};
        interpreter='tex';
    case 6
        M=varargin{1};
        hf=varargin{2};
        numDigits=varargin{3};
        fontSize=varargin{4};
        textColor=varargin{5};
        interpreter=varargin{6};
end

if isempty(hf)
    hf=gcf;
end

% textFormat=['%6.',num2str(numDigits),'e'];
textFormat=['%0.',num2str(numDigits),'f';];

%%
figure(hf);
H=zeros(1,numel(M));
qh=1;
for qi=1:size(M,1)
    for qj=1:size(M,2)
        m=M(qi,qj);
        switch class(m)
            case'sym'
                try %to convert to double
                    m=double(m);
                    image_text=sprintf(textFormat,m);
                catch %Contains symbolic expressions
                    switch interpreter
                        case 'none'
                            image_text=char(m);
                        case 'latex'
                            image_text=latex(m);
                        case 'tex'
                            image_text=texlabel(m);
                    end
                    %image_text = strsplit(image_text,' '); %Places space
                    %separated elements on new line
                end
            otherwise
                image_text=sprintf(textFormat,m);
        end
        H(qh)=text(qj,qi, image_text,'horizontalAlignment','center','color',textColor,'FontWeight','demi','FontSize',fontSize,'interpreter',interpreter);
        qh=qh+1;
    end
end
if nargout==1
    varargout{1}=H;
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [c,h]=tricontour(tri,x,y,z,nv)
%TRICONTOUR Triangular Contour Plot.
% TRICONTOUR(TRI,X,Y,Z,N) draws scalar N contour lines treating the values
% in Z as heights above a plane. TRI,X,Y,and Z define a triangulation where
% the triangles are defined by the M-by-3 face matrix TRI, such as that
% returned by DELAUNAY. Each row of TRI contains indices into the X,Y, and
% Z vertex vectors to define a single triangular face. Contours are
% computed directly from the triangulation rather than interpolating back
% to a cartesian grid using GRIDDATA.
% TRICONTOUR(TRI,X,Y,Z,V) draws length(V) contour lines at the values
% specified in vector V.
% TRICONTOUR(TRI,X,Y,Z,[v v]) draws a single contour line at the level v.
%
% [C,H] = TRICONTOUR(...) returns contour matrix C as described in CONTOURC
% and a vector of handles H to the created patch objects.
% H can be used to set patch properties.
% CLABEL(C) or CLABEL(C,H) labels the contour levels.
%
% Example:
%           x=linspace(-3,3,39);
%           y=linspace(-2.5,2.5,49);
%           [xx,yy]=meshgrid(x,y);
%           zz=peaks(xx,yy);
%           v=-3:1:5; % contour levels
%           subplot(1,2,1)
%           [C,h]=contour(xx,yy,zz,v);   % standard contour for comparison
%           clabel(C)
%           title Contour
%
%           idx=randperm(numel(zz));     % grab some scattered indices
%           n=idx(1:ceil(numel(zz)/2))'; % one half of them
%           x=xx(n);                     % get scattered data
%           y=yy(n);
%           z=zz(n);
%           tri=delaunay(x,y);           % triangulate scattered data
%           subplot(1,2,2)
%           [C,h]=tricontour(tri,x,y,z,v);
%           clabel(C,h)
%           title TriContour
%
% view(3) displays the contour in 3-D.
%
% See also DELAUNAY, CONTOUR, TRIMESH, TRISURF, TRIPLOT, PATCH.

% D.C. Hanselman, University of Maine, Orono, ME 04469
% MasteringMatlab@yahoo.com
% Mastering MATLAB 7
% 2006-05-07, 2006-05-16, 2006-07-25

if nargin<5
	error('Not Enough Input Arguments.')
end
x=x(:);	% convert input data into column vectors
y=y(:);
z=z(:);
xlen=length(x);
if ~isequal(xlen,length(y),length(z))
   error('X, Y, and Z Must Have the Same Number of Elements.')
end
if size(tri,2)~=3 || any(tri(:)<0) || any(tri(:)>xlen)
   error('TRI Must Be a Valid Triangulation of the Data in X, Y, Z.')
end

zs=z(tri);
zmax=max(max(zs));              % find max and min in z data that is in tri
zmin=min(min(zs));

if length(nv)==1                                 % nv is number of contours
   zlev=linspace(zmax,zmin,nv+2);
elseif length(nv)==2 && nv(1)==nv(2)              % nv is one contour level
   zlev=nv(1);
else                                       % nv is vector of contour levels
   zlev=sort(nv,'descend');
end
zlev(zlev>=zmax | zlev<=zmin)=[];  % eliminate contours outside data limits
nlev=length(zlev);

if nlev==0
   error('No Contours to Plot. Chosen Contours Outside Limits of Data.')
end

% precondition the input data
[zs,zidx]=sort(zs,2);         % sort vertices by z value ascending
for k=1:size(zs,1)            % shuffle triangles to match
   tri(k,:)=tri(k,zidx(k,:));
end

hax=newplot;                  % create new axis if needed
h=[];                         % patch handle storage
C=zeros(2,0);                 % Clabel data storage
cs=[2 1];                     % column swap vector cs(1)=2, cs(2)=1;

% Main Loop ---------------------------------------------------------------
for v=1:nlev                  % one contour level at a time
   zc=zlev(v);                % chosen level
   above=zs>=zc;              % true for vertices above given contour
   numabove=sum(above,2);     % number of triangle vertices above contour
   tri1=tri(numabove==1,:);   % triangles with one vertex above contour
   tri2=tri(numabove==2,:);   % triangles with two vertices above contour
   n1=size(tri1,1);           % number with one vertex above
   n2=size(tri2,1);           % number with two vertices above

   edge=[tri1(:,[1 3])        % first column is indices below contour level
         tri1(:,[2 3])        % second column is indices above contour level
         tri2(:,[1 2])
         tri2(:,[1 3])];
   if n1==0                   % assign edges to triangle number
      n=[1:n2 1:n2]';
   elseif n2==0
      n=[1:n1 1:n1]';
   else
      n=[1:n1 1:n1 n1+(1:n2) n1+(1:n2)]';
   end

   [edge,idx]=sortrows(edge);    % put shared edges next to each other
   n=n(idx);                     % shuffle triangle numbers to match

   idx=all(diff(edge)==0,2);     % find shared edges
   idx=[idx;false]|[false;idx];  % True for all shared edges

   % eliminate redundant edges, two triangles per interior edge
   edgeh=edge(~idx,:);           % hull edges
   nh=n(~idx);                   % hull triangle numbers
   if ~isempty(nh)
      nh(end,2)=0;               % zero second column for hull edges
   end
   edges=edge(idx,:);            % shared edges
   edges=edges(1:2:end-1,:);     % take only unique edges
   ns=n(idx);                    % interior triangle numbers
   ns=[ns(1:2:end) ns(2:2:end)]; % second column is second triangle
   edge=[edgeh;edges];           % unique edges
   nn=[nh;ns];                   % two columns of triangle numbers
   ne=size(edge,1);              % number of edges

   flag=true(ne,2);              % true for each unused edge per triangle
   tmp=zeros(ne+1,1);            % contour data temporary storage

   xe=x(edge);                   % x values at vertices of edges
   ye=y(edge);                   % y values at  vertices of edges
   ze=z(edge);                   % z data at  vertices of edges

   alpha=(zc-ze(:,1))./(ze(:,2)-ze(:,1)); % interpolate all edges
   xc=alpha.*(xe(:,2)-xe(:,1)) + xe(:,1); % x values on this contour
   yc=alpha.*(ye(:,2)-ye(:,1)) + ye(:,1); % y values on this contour

   while any(flag)	% while there are still unused edges -----------------

      xtmp=tmp;
      ytmp=tmp;
      [ir,ic]=find(flag,1);            % find next unused edge
      flag(ir,ic)=false;               % mark this edge used

      k=1;                             % first data point in subcontour
      xtmp(k)=xc(ir);                  % store data from this edge
      ytmp(k)=yc(ir);

      while true     % complete this subcontour ---------------------------

         [ir,ic]=find(flag&nn(ir,ic)==nn,1);% find other edge of triangle
         flag(ir,ic)=false;            % mark this edge used
         k=k+1;
         xtmp(k)=xc(ir);               % store data from this edge
         ytmp(k)=yc(ir);

         ic=cs(ic);                    % other triangle that shares edge

         if nn(ir,ic)==0               % reached hull, subcontour complete
            k=k+1;
            xtmp(k)=nan;               % don't let subcontour close
            ytmp(k)=nan;
            break
         elseif ~flag(ir,ic)           % complete closed subcontour
            break
         else                          % more points remain on subcontour
            flag(ir,ic)=false;         % mark this edge used
         end
      end % while true ----------------------------------------------------
      xtmp(k+1:end)=[];                % throw away unused storage
      ytmp(k+1:end)=[];                % xtmp,ytmp contain subcontour

      if nargout<2                     % plot the subcontour
         patch('XData',xtmp,'YData',ytmp,'CData',repmat(zc,k,1),...
               'Parent',hax,'FaceColor','none','EdgeColor','flat',...
               'UserData',zc)
         C=horzcat(C,[zc xtmp';k ytmp']); % contour label data
      else                             % plot subcontour and create output
         h=[h;patch('XData',xtmp,'YData',ytmp,'CData',repmat(zc,k,1),...
         'Parent',hax,'FaceColor','none','EdgeColor','flat',...
         'UserData',zc)]; %#ok
         C=horzcat(C,[zc xtmp';k ytmp']); % contour label data
      end
   end % while any(flag) --------------------------------------------------
end % for v=1:nlev
if nargout
   c=C;
end


text:
function [ end_index, min_index ] = getTrough( time_slice, begin_index, threshold )
%GETTROUGH Get trough for cyclonic eddy
%   Detailed explanation goes here

end_index = begin_index;
n = length(time_slice);
while end_index < n && time_slice(end_index) >= time_slice(end_index+1) + threshold
    end_index = end_index + 1;
end

min_index = end_index;

while end_index < n && (time_slice(end_index) - time_slice(end_index+1)) <= threshold
    end_index = end_index + 1;
end

end

author: jfaghm

repo: OceanEddies

save_path: github-repos/MATLAB/jfaghm-OceanEddies


text:
function [bodyX, bodyY, bodyZ, R_total] = computeBodyAxesFromEuler(ut, rVect, vVect, bodyInfo, rollAng, pitchAng, yawAng)
    R_body_2_ned = eul2rotmARH([yawAng,pitchAng,rollAng],'zyx');
    [R_ned_2_inert, ~, ~, ~] = computeNedFrame(ut, rVect, bodyInfo);
	R_total = R_ned_2_inert * R_body_2_ned;

    bodyX = R_total(:,1);
    bodyY = R_total(:,2);
    bodyZ = R_total(:,3);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function x = fac_to_vec(A)
%FAC_TO_VEC Converts a set of factor matrices to a vector.
%
%   X = FAC_TO_VEC(A) converts a cell array of factor matrices A to a
%   vector by vectorizing each matrix and stacking them.
%
%   See also CP_VEC_TO_FAC, CP_FUN, CP_OPT.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Set-up
N = length(A);

%% Get sizes
sz = zeros(N,1);
for n = 1:N
    sz(n) = size(A{n},1);
end
R = size(A{1},2);
P = sum(sz)*R;

%% Create x
x = zeros(P,1);
for n = 1:N
    idx1 = sum(sz(1:n-1))*R + 1;
    idx2 = sum(sz(1:n))*R;
    x(idx1:idx2) = reshape(A{n},sz(n)*R,1);
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function hFig = JTattooDemo()
% JTattooDemo - demonstrates effects of Look-and-Feel changes (including JTattoo L&F)
%
% JTattooDemo displays a figure with Matlab uicontrols and Java components,
% and allows the user to modify the figure's look-and-feel (L&F) to any of
% the standard L&Fs installed on the system, as well to the JTattoo L&Fs.
%
% Syntax:
%    hFig = JTattooDemo()
%
% Input Parameters:
%    None.
%
% Output parameters:
%    hFig - handle of the created figure
%
% Technical description:
%    http://UndocumentedMatlab.com/blog/jtattoo-look-and-feel-demo/
%    http://www.jtattoo.net
%
% Bugs and suggestions:
%    Please send to Yair Altman (altmany at gmail dot com)
%
% Release history:
%    1.1 2013-03-19: Fixed JTattoo's name; restored original L&F for other GUI; removed slider Java control; added current L&F title
%    1.0 2013-03-19: First version posted on MathWorks File Exchange

% License to use and modify this code is granted freely to all interested, as long as the original author is
% referenced and attributed as such. The original author maintains the right to be solely associated with this work.

% Programmed and Copyright by Yair M. Altman: altmany(at)gmail.com
% $Revision: 1.0 $  $Date: 2013/03/19 23:22:38 $

    % Add the JTatto java library
    folder = fileparts(mfilename('fullpath'));
    javaclasspath( {folder, [folder '\JTattoo.jar']} );

    % The main figure can't be moved or resized with setting the default Windows L&F
    %originalLnF = javax.swing.UIManager.getLookAndFeel;
    %javax.swing.UIManager.setLookAndFeel(originalLnF); %'com.sun.java.swing.plaf.windows.WindowsLookAndFeel');

    % Create MATLAB figure
    hFig = figure( 'Name', 'JTattoo Look-and-Feel Demo', ...
                   'NumberTitle', 'off', ...
                   'MenuBar', 'none', ...
                   'Toolbar', 'none', ...
                   'CloseRequestFcn',@onExit);

    % Get the MATLAB figure's underlying Java frame
    javaWindow = getJavaFrame( hFig );
    setappdata( hFig, 'javaWindow', javaWindow );

    % Add the figure's main menus
    addMainMenus( hFig );

    % Add the figure's controls
    addControls( hFig );

end  % JTattooDemo

% Add the main menus to the figure
function addMainMenus( hFig )
    % File main menu
    FileMenu = uimenu( hFig, 'Label', 'File' );
    uimenu( FileMenu, 'Label', 'Exit', 'Callback', @onExit );

    % System main menu
    SystemMenu = uimenu( hFig, 'Label', 'System' );
    lnfs = javax.swing.UIManager.getInstalledLookAndFeels();
    for idx = 1 : length( lnfs )
        uimenu( SystemMenu, 'Label',char(lnfs(idx).getName), 'Callback',@(h,e)updateInterface(char(lnfs(idx).getClassName)));
    end

    % JTattoo look and feel object names
    jt_lnf =   {
        'Default'     ''
        'Smart'       'com.jtattoo.plaf.smart.SmartLookAndFeel'
        'Aluminium'   'com.jtattoo.plaf.aluminium.AluminiumLookAndFeel'
        'Acryl'       'com.jtattoo.plaf.acryl.AcrylLookAndFeel'
        'Aero'        'com.jtattoo.plaf.aero.AeroLookAndFeel'
        'Bernstein'   'com.jtattoo.plaf.bernstein.BernsteinLookAndFeel'
        'Graphite'    'com.jtattoo.plaf.graphite.GraphiteLookAndFeel'
        'Fast'        'com.jtattoo.plaf.fast.FastLookAndFeel'
        'Hifi'        'com.jtattoo.plaf.hifi.HiFiLookAndFeel'
        'Luna'        'com.jtattoo.plaf.luna.LunaLookAndFeel'
        'McWin'       'com.jtattoo.plaf.mcwin.McWinLookAndFeel'
        'Mint'        'com.jtattoo.plaf.mint.MintLookAndFeel'
        'Noire'       'com.jtattoo.plaf.noire.NoireLookAndFeel'
        };

    % JTattoo main menu
    JTattooMenu = uimenu( hFig, 'Label', 'JTattoo' );
    for idx = 1 : length( jt_lnf )
        uimenu( JTattooMenu, 'Label',jt_lnf{idx,1}, 'Callback',@(h,e)updateInterface(jt_lnf{idx,2})) ;
    end
end  % addMainMenus

% Add the figure's controls
function addControls( hFig )
    % Example MATLAB controls
    matlabPanel = uipanel(hFig, 'units','norm', 'pos',[0.05,0.05,0.4,0.9], 'title','Matlab controls');
    controlLayout = uiflowcontainer( 'v0', matlabPanel );
    uicontrol(controlLayout, 'Style','list', 'BackgroundColor', 'w', 'String',{'Listbox item #1','Listbox item #2','Listbox item #3'}, 'Parent', controlLayout);
    uicontrol(controlLayout, 'Style','Pushbutton',  'String', 'Button ');
    uicontrol(controlLayout, 'Style','edit',        'String', 'Edit Box');
    uicontrol(controlLayout, 'Style','text',        'String', 'Static Text');
    uicontrol(controlLayout, 'Style','slider',      'String', 'Slider');
    uicontrol(controlLayout, 'Style','radiobutton', 'String', 'Radio Button');
    uicontrol(controlLayout, 'Style','checkbox',    'String', 'Check Box');
    uicontrol(controlLayout, 'Style','popupmenu',   'String', {'Option 1','Option 2', 'Option 3'});

    % Example Java components
    javaPanel = uipanel(hFig, 'units','norm', 'pos',[0.5,0.05,0.4,0.9], 'title','Java controls');
    jspinner     = javax.swing.JSpinner();
    jspinner.setValue(25);
    jbutton      = javax.swing.JButton('Java Button');
    jtextfield   = javax.swing.JTextField('Java TextField');
    jprogressbar = javax.swing.JProgressBar(0,100);
    jprogressbar.setStringPainted(true);
    jprogressbar.setValue(65);
    jslider      = javax.swing.JSlider(0,100,35);
    jslider.setPaintLabels(true);
    jslider.setMajorTickSpacing(20);
    jscrollbar = javax.swing.JScrollBar();
    jscrollbar.setValue(40);
    jscrollbar.setOrientation(jscrollbar.HORIZONTAL);
    jradiobutton = javax.swing.JRadioButton('Java RadioButton');
    jcheckbox    = javax.swing.JCheckBox('Java CheckBox');
    jcombobox    = javax.swing.JComboBox({'Option #1','Option #2','Option #3'});
    placeJavaComponent(jspinner,     [0.1 0.9  0.8 0.07], javaPanel);
    placeJavaComponent(jbutton,      [0.1 0.8  0.8 0.07], javaPanel);
    placeJavaComponent(jtextfield,   [0.1 0.7  0.8 0.07], javaPanel);
    placeJavaComponent(jprogressbar, [0.1 0.60 0.8 0.07], javaPanel);
   %placeJavaComponent(jslider,      [0.1 0.45 0.8 0.10], javaPanel);  % buggy in JTattoo, so leave out
    placeJavaComponent(jscrollbar,   [0.1 0.35 0.8 0.07], javaPanel);
    placeJavaComponent(jradiobutton, [0.1 0.25 0.8 0.07], javaPanel);
    placeJavaComponent(jcheckbox,    [0.1 0.15 0.8 0.07], javaPanel);
    placeJavaComponent(jcombobox,    [0.1 0.05 0.8 0.07], javaPanel);

    % Display the current L&F name
    bgcolor = get(hFig,'color');
    lnfName = char(javax.swing.UIManager.getLookAndFeel.getName);
    uicontrol(hFig, 'Style','text', 'FontWeight','bold', 'units','norm', 'pos',[0.02,0.95,0.9,0.04], 'string',lnfName, 'Background',bgcolor, 'tag','L&F title');
end  % addControls

% Place a Java component within a panel in a specified normalized position
function placeJavaComponent ( jcomponent, position, parent )
    jcomponent = javaObjectEDT( jcomponent );  % ensure component is auto-delegated to EDT
    jcomponent.setOpaque(false);  % useful to demonstrate L&F backgrounds
    [jc,hContainer] = javacomponent( jcomponent, [], parent ); %#ok<ASGLU>
    set(hContainer, 'Units','Normalized', 'Position',position);
end  % placeJavaComponent

% User wants to exit the application
function onExit( h, e ) %#ok<INUSD>
    hFig = gcbf;
    try
        javaWindow = getappdata( hFig, 'javaWindow' );
        SetDefaultLookAndFeel( javaWindow );
    catch
        % ignore...
    end
    delete( hFig );
    drawnow;
end % onExit

% User requested to change the L&F
function updateInterface( lookandfeel )
    originalLnF = javax.swing.UIManager.getLookAndFeel;
    javaWindow = getappdata( gcbf, 'javaWindow' );
    if isempty( lookandfeel )
        javax.swing.UIManager.setLookAndFeel(...
            'com.sun.java.swing.plaf.windows.WindowsLookAndFeel');
    else
        jTattoo = javaObjectEDT( 'UpdateJTattooInterface' );
        javaMethodEDT('setLookAndFeel', jTattoo, lookandfeel, 'Default', javaWindow);
    end
    container = javaMethodEDT('getContentPane', javaWindow);
    javaMethodEDT('updateComponentTreeUI', 'javax.swing.SwingUtilities', container);

    % Update the title
    hTitle = findall(gcbf, 'tag','L&F title');
    lnfName = char(javax.swing.UIManager.getLookAndFeel.getName);
    set(hTitle, 'String',lnfName);
    drawnow;

    % Restore original L&F for other GUI
    javax.swing.UIManager.setLookAndFeel(originalLnF)
end % updateInterface

% Get a Matlab figure's underlying Java Frame (RootPane) reference handle
function JavaFrame = getJavaFrame( hFig )
    try
        %  contentSize = [0,0];  % initialize
        JavaFrame = hFig;
        figName = get(hFig,'name');
        if strcmpi(get(hFig,'number'),'on')
            figName = regexprep(['Figure ' num2str(hFig) ': ' figName],': $','');
        end
        mde = com.mathworks.mde.desk.MLDesktop.getInstance;
        jFigPanel = mde.getClient(figName);
        JavaFrame = jFigPanel;
        JavaFrame = jFigPanel.getRootPane;
    catch
        try
            jFrame = get(handle(hFig),'JavaFrame');
            jFigPanel = get(jFrame,'FigurePanelContainer');
            JavaFrame = jFigPanel;
            JavaFrame = jFigPanel.getComponent(0).getRootPane;
        catch
            % Never mind
        end

    end
    try
        % If invalid RootPane, retry up to N times
        tries = 10;
        while isempty(JavaFrame) && tries>0  % might happen if figure is still undergoing rendering...
            drawnow; pause(0.001);
            tries = tries - 1;
            JavaFrame = jFigPanel.getComponent(0).getRootPane;
        end

        % If still invalid, use FigurePanelContainer which is good enough in 99% of cases... (menu/tool bars won't be accessible, though)
        if isempty(JavaFrame)
            JavaFrame = jFigPanel;
        end
        % contentSize = [JavaFrame.getWidth, JavaFrame.getHeight];

        % Try to get the ancestor FigureFrame
        JavaFrame = JavaFrame.getTopLevelAncestor;
    catch
        % Never mind - FigurePanelContainer is good enough in 99% of cases... (menu/tool bars won't be accessible, though)
    end
end  % getJavaFrame


text:
function [h,s,v] = rgb2hsv_fast(r,g,b,opt_typ,opt_sel)
%RGB2HSV_FAST does the same as RGB2HSV
%but faster and less memory exhaustive,
%especially on images which are stored as integers.
%
%RGB2HSV_FAST(IMG,'single')
%You can specify the optional argument 'single'
%to use single precision.
%This saves memory and is sufficiently exact, anyhow.
%
%RGB2HSV_FAST(IMG,'','H')
%Optionally, you can specify which channels to calculate:
%'HSV' is default.
%'H' calculates hue only.
%'S' calculates saturation only.
%'V' calculates value only.
%'HS', 'HV', or 'SV' are also possible.

% How to make it efficient?
%  -  Don't cast integer images to floating point too soon.
%     Operate on integers as far as possible and
%     take advantage of functions like IMLINCOMB().
%  -  Avoid FIND(). Directly index via binary arrays instead.
%  -  Avoid RESHAPE().
%     http://blogs.mathworks.com/loren/?p=28

%
% Alexander Ihlow, Germany, Oct 2006
% Version 1.1
% 2006-11-02 added support for individual channel selection
%            version 1.0 saved to rgb2hsv_fast_10.m
%
%
% RGB2HSV_FAST has a high cyclomatic complexity
% due to lots of case checking.
%
% Please report bugs or suggestions to ml-user@web.de.

threeD = ndims(r) >= 3;

if (nargin < 4), opt_typ = ''; end
if (nargin < 5), opt_sel = ''; end
if (nargin > 1) && ~isnumeric(g), opt_typ = g; end
if (nargin > 2) && ~isnumeric(b), opt_sel = b; end

%disp(['threeD = ' num2str(threeD)])
%disp(['opt_typ = ' opt_typ])
%disp(['opt_sel = ' opt_sel])

vv = version; vv=str2double(vv(1));
% Check input class
switch (class(r))
  case 'uint8'
    immax = 255;
  case 'uint16'
    immax = 65535;
  case 'double'
    immax = 1;
  case 'single'
    if vv < 7
      error('Matlab version < 7 cannot perform calculations with single precision. Please cast to double first!')
    else
      immax = 1;
      % If input is single, output single per default.
      if isempty(opt_typ), opt_typ = 'single'; end
    end
  otherwise
    if vv < 7
      error('Only images of type uint8, uint16, and double are supported!')
    else
      error('Only images of type uint8, uint16, single, and double are supported!')
    end
end

% Default output is double.
if isempty(opt_typ), opt_typ = 'double'; end

% Use single precision if optionally specified.
% This saves 1/2 of memory and is exact enough
% for those simple colorspace transformations.
% single : out_typ = 1
% double : out_typ = 2
switch lower(opt_typ(1))
  case 's'
    out_typ = 1;
  otherwise
    out_typ = 2;
end

% For Matlab < version 7, no single precision is supported.
if (out_typ == 1) && (vv < 7)
  disp('Single precision is not available in this Matlab version. Using double.')
  out_typ = 2;
end

% Helper functions to map float either to single or double.
if out_typ == 1
  im2float = inline('im2single(x)');
  float = inline('single(x)');
  floatstr = 'single';
else
  im2float = inline('im2double(x)');
  float = inline('double(x)');
  floatstr = 'double';
end

immax = float(immax);

if threeD
  % Split RGB image into R G B channels and
  % thereby support arbitrary dimensional images.
  % R G B is assumed to be contained in the third dimension.
  sizc = num2cell(size(r));
  for k=1:numel(sizc), sizc{k} = ':'; end
  sizc{3} = 2;
  g = r(sizc{:});
  sizc{3} = 3;
  b = r(sizc{:});
  sizc{3} = 1;
  r = r(sizc{:});
elseif nargin==1
  % split N x 3 matrix into three vectors
  g = r(:,2); b = r(:,3); r = r(:,1);
end
siz = size(r);

if isempty(opt_sel), opt_sel = 'HSV'; else opt_sel = upper(opt_sel); end

if any( (opt_sel ~= 'H') & (opt_sel ~= 'S') & (opt_sel ~= 'V') )
  error('Invalid optional argument! Please use ''H'', ''S'', or ''V''')
end

if exist('imsubtract','file')
  % Image Processing Toolbox available
  % --> The calculation of HUE and SAT
  %     is performed without precasting R, G, and B to floats.
  v = max(max(r,g),b);
  if any(opt_sel == 'H') || any(opt_sel == 'S')
    s = imsubtract(v, min(min(r,g),b));
    z = ~s;
    s = im2float(s);
    s(z) = 1;

    if any(opt_sel == 'H')
      % Calculating HUE via IMLINCOMB() is highly efficient.
      sizc = num2cell(siz);
      h(sizc{:}) = float(0);
      k = (r == v);
      h(k) =     imlincomb(1,g(k),-1,b(k),floatstr)./(immax*s(k));
      k = (g == v);
      h(k) = 2 + imlincomb(1,b(k),-1,r(k),floatstr)./(immax*s(k));
      k = (b == v);
      h(k) = 4 + imlincomb(1,r(k),-1,g(k),floatstr)./(immax*s(k));

      h = (1/6) * h;
      k = (h < 0);
      h(k) = h(k) + 1;

      h(z) = 0;
    end

    if any(opt_sel == 'S') || any(opt_sel == 'V')
      k = (v ~= 0);
      v = im2float(v);
      if any(opt_sel == 'S')
        s(k) = (~z(k)).*s(k)./v(k);
        s(~k) = 0;
      end
    end
  else
    v = im2float(v);
  end

else
  % Image Processing Toolbox not available
  % --> Use suboptimal procedure.
  v = im2float(max(max(r,g),b));
  if any(opt_sel == 'H') || any(opt_sel == 'S')
    s = v - im2float(min(min(r,g),b));
    z = ~s;
    s(z) = 1;

    sizc = num2cell(siz);
    h(sizc{:}) = float(0);

    if any(opt_sel == 'H')
      k = (im2float(r) == v);
      h(k) =     (im2float(g(k))-im2float(b(k)))./s(k);
      k = (im2float(g) == v);
      h(k) = 2 + (im2float(b(k))-im2float(r(k)))./s(k);
      k = (im2float(b) == v);
      h(k) = 4 + (im2float(r(k))-im2float(g(k)))./s(k);

      h = (1/6) * h;
      k = (h < 0);
      h(k) = h(k) + 1;

      h(z) = 0;
    end

    if any(opt_sel == 'S')
      k = (v ~= 0);
      s(k) = (~z(k)).*s(k)./v(k);
      s(~k) = 0;
    end
  end

end

if ~any(opt_sel == 'H'), h = float([]); end
if ~any(opt_sel == 'S'), s = float([]); end
if ~any(opt_sel == 'V'), v = float([]); end

if nargout <= 1,
  if (threeD || nargin>=3),
    h = cat(3,h,s,v);
  else
    h = [h s v];
  end
else
  switch opt_sel
    case 'HSV'
      %
    case 'H'
      disp('There is more than one output variable given, but only HUE is assigned to the first one!')
    case 'S'
      h = s;
      disp('There is more than one output variable given, but only SAT is assigned to the first one!')
    case 'V'
      h = v;
      disp('There is more than one output variable given, but only VAL is assigned to the first one!')
    case 'HS'
      if nargout > 2
        disp('There are more than two output variables given, but only HUE and SAT are assigned to the first two!')
      end
    case 'HV'
      s = v;
      if nargout > 2
        disp('There are more than two output variables given, but only HUE and VAL are assigned to the first two!')
      end
    case 'SV'
      h = s;
      s = v;
      if nargout > 2
        disp('There are more than two output variables given, but only SAT and VAL are assigned to the first two!')
      end
    otherwise
      disp('Unrecognized output - please check what happened!')
  end
end


text:
% -*- INTERNAL UNDOCUMENTED FUNCTION -*-
function [spv, spp] = sp_bspline_fluid_3d (element_name, ...
                   knots, nsub_p, degree_p, regularity_p, msh)

warning ('geopdes:obsolete', 'Function SP_BSPLINE_FLUID_3D is obsolete. Using SP_BSPLINE_FLUID instead')

[spv, spp] = sp_bspline_fluid (element_name, knots, nsub_p, degree_p, regularity_p, msh);

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [A, B] = kafbox_covNoise(logtheta, x, z);

% Independent covariance function, ie "white noise", with specified variance.
% The covariance function is specified as:
%
% k(x^p,x^q) = s2 * \delta(p,q)
%
% where s2 is the noise variance and \delta(p,q) is a Kronecker delta function
% which is 1 iff p=q and zero otherwise. The hyperparameter is
%
% logtheta = [ log(sqrt(s2)) ]
%
% For more help on design of covariance functions, try "help covFunctions".
%
% (C) Copyright 2006 by Carl Edward Rasmussen, 2006-03-24.

if nargin == 0, A = '1'; return; end              % report number of parameters

s2 = exp(2*logtheta);                                          % noise variance

if nargin == 2                                      % compute covariance matrix
  A = s2*eye(size(x,1));
elseif nargout == 2                              % compute test set covariances
  A = s2;
  B = 0;                               % zeros cross covariance by independence
else                                                % compute derivative matrix
  A = 2*s2*eye(size(x,1));
end

repo: kafbox


text:
%-----------------------------------------------------------------------
% Job configuration created by cfg_util (rev $Rev: 3355 $)
%-----------------------------------------------------------------------
matlabbatch{1}.cfg_toy{1}.add2{1}.cfg_example_add1.a = '<UNDEFINED>';
matlabbatch{1}.cfg_toy{1}.add2{1}.cfg_example_add1.b = '<UNDEFINED>';
matlabbatch{2}.cfg_toy{1}.cfg_example_div.a = '<UNDEFINED>';
matlabbatch{2}.cfg_toy{1}.cfg_example_div.b(1) = cfg_dep;
matlabbatch{2}.cfg_toy{1}.cfg_example_div.b(1).tname = 'Input b';
matlabbatch{2}.cfg_toy{1}.cfg_example_div.b(1).tgt_spec = {};
matlabbatch{2}.cfg_toy{1}.cfg_example_div.b(1).sname = 'Add1: Add1: a + b';
matlabbatch{2}.cfg_toy{1}.cfg_example_div.b(1).src_exbranch = substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{2}.cfg_toy{1}.cfg_example_div.b(1).src_output = substruct('()',{1});
matlabbatch{3}.cfg_basicio.cfg_assignin.name = 'div_sum_ab_c';
matlabbatch{3}.cfg_basicio.cfg_assignin.output(1) = cfg_dep;
matlabbatch{3}.cfg_basicio.cfg_assignin.output(1).tname = 'Output Item';
matlabbatch{3}.cfg_basicio.cfg_assignin.output(1).tgt_spec = {};
matlabbatch{3}.cfg_basicio.cfg_assignin.output(1).sname = 'div: a div b: mod';
matlabbatch{3}.cfg_basicio.cfg_assignin.output(1).src_exbranch = substruct('.','val', '{}',{2}, '.','val', '{}',{1});
matlabbatch{3}.cfg_basicio.cfg_assignin.output(1).src_output = substruct('.','mod');

author: spm


text:
function [F,totalCost,exitCode]=minCostFlow(AMat,CMat,b,maxIter)
%%MINCOSTFLOW Solve the minimum cost flow problem using a strong polynomial
%             time cycle cancelling algorithm that works with non-integer
%             costs. Minimum cost flow problems include transportation
%             problems.
%
%INPUTS: AMat An NXN matrix of costs in a directed graph such that
%             AMat(i,j) is the cost of an edge going from vertex i to
%             vertex j. If no vertex goes from node i to node j, then a
%             cost of 0 should be inserted.
%        CMat An NXN matrix of capacities of the edges in the graph.
%             CMat(i,j) is the capacity of an edge from node i to node j.
%             If no edge exists, then a capacity of zero should be used. It
%             is assumed that all capacities are finite.
%           b An NX1 vector of the supply provided by each node. It is
%             required that sum(b)=0.
%     maxIter An optional parameter specifying the maximum number of
%             iterations that should be performed by the algorithm. If
%             omitted, a maximum of 100+10N iterations is used.
%
%OUTPUTS: F The flow matrix. F(i,j) is the amount of flow going from node i
%           to node j. Note that F(i,j)=-F(j,i). If the algorithm could not
%           obtain an initial feasible solution, then an empty matrix is
%           returned. If the maximum number of iterations is exceeded, then
%           the algorithm will return a feasible solution that is not
%           optimal.
% totalCost The total cost of the flow. This is the quantity being
%           minimized. If the algorithm could not obtain an initial
%           feasible solution, then an empty matrix is returned.
%  exitCode A parameter indicating whether an error occurred or whether the
%           algorithm terminated successfully. Possible values are:
%           0 The algorithm was successful.
%           1 No feasible solution could be found.
%           2 The maximum number of iterations was reached.
%           3 Inputs are invalid. This means that either C contains NaN
%             values, sum(b)~=0 or AMat contains nonfinite values.
%
%The minimum cost flow problems seeks to find a flow F to minimize
%sum(sum(AMat.*F)) subject to the constraints:
%4) Total Vertex Flow: sum(F(i,:))=b(i) for all i.
%2) Capacity Constraint: 0<=F(u,v)<=CMat(u,v)
%3) Skew Symmetry: F(u,v)=-F(v,u)
%Whereas the maximum flow problem, which is solved in the function
%solveMaxFlowEdmondsKarp seeks to maximize the flow between two nodes the
%minimum cost flow problem seeks to minimize the total cost of a fixed
%amount of flow.
%
%The strong polynomial cycle cancelling algorithm of [1] is used. However,
%the epsilon-based convergence acceleration method described in [1] is not
%used. Rather, the algorithm is more similar to the simple cycle
%cancellation algorithm of Chapter 7.2 of [2], except the minimum mean
%cycle is always augmented rather than any negative cost cycle.
%
%The complexity of the non-epsilon accelerated algorithm in [1] is bounded
%as O(N^2*m^3*log(N)), where N is the number of vertices in the graph, and
%m is the number of edges. The true complexity depends on the efficiency of
%the function computeResidualCapacity, which is used to obtain the negative
%cycles.
%
%EXAMPLE: This is example 7.1 in Chapter 7.2 of [2].
% AMat=[0, 4, 1, 0;
%       0, 0, 2, 5;
%       0, 3, 0, 2;
%       0, 0, 0, 0];
% CMat=[0, 2, 2, 0;
%       0, 0, 1, 1;
%       0, 1, 0, 1;
%       0, 0, 0, 0];%Capacity matrix
% b=[2;0;0;-2];
% [F,totalCost,exitCode]=minCostFlow(AMat,CMat,b)
%One will get a flow matrix of
% F =[ 0     0     2     0;
%      0     0    -1     1;
%     -2     1     0     1;
%      0    -1    -1     0];
%having a total cost of 12.
%
%OTHER EXAMPLES: Transportation problems can be transformed into minimum
%cost flow problems and vice versa, as shown in Chapters 7.4 and 7.5 of
%[2]. Thus, the examples in the comments to the function
%solveTransportationProblem can be taken as additional examples of the
%minimum cost flow problem.
%
%REFERENCES:
%[1] A. V. Goldberg and R. E. Tarjan, "Finding minimum-cost circulations
%    by canceling negative cycles," Journal of the Association for
%    Computing Machinery, vol. 36, no. 4, pp. 873-886, Oct. 1989.
%[2] C. H. Papadimitriou and K. Steiglitz, Combinatorial Optimization:
%    Algorithms and Complexity. Englewood Cliffs, NJ: Prentice-Hall Inc.,
%    1982.
%
%July 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numVertex=length(b);

if(nargin<5||isempty(maxIter))
    maxIter=100+10*numVertex;
end

if(any(isnan(CMat(:)))||sum(b)~=0||any(~isfinite(AMat(:))))
    exitCode=3;%The problem is infeasible.
    F=[];
    totalCost=[];
    return;
end

%Step 1: Find a feasible flow.
%In this step, we must find a feasible flow F. This is a flow that
%satisfies constraints 1-3 but does not necessarily minimize the cost
%sum(sum(AMat.*F)). This can be done by ignoring costs and solving a
%maximum flow problem based only on the constraints. The maximum flow
%problem does not have demand constraints on all of the nodes, rather it
%only has them on a sink and source which are added.
%
%To enforce constraint 1, we have to create a network with an artificial
%sink and source. If b(i) is positive, then we need an arc from the source
%to that node with capacity b(i). If b(i) is negative, then we need an arc
%from the sink to that node with capacity -b(i). If the maximum flow
%problem solved on the modified graph saturates all of the edges coming out
%of the source, then the problem is feasible.

%Augment the capacity matrix.
CMatAugment=[CMat,zeros(numVertex,2);
             zeros(2,numVertex+2)];

source=numVertex+1;
sink=numVertex+2;
totalFlowIn=0;
for curVertex=1:numVertex
    if(b(curVertex)>0)
        CMatAugment(source,curVertex)=b(curVertex);
        totalFlowIn=totalFlowIn+b(curVertex);
    else
        CMatAugment(curVertex,sink)=-b(curVertex);
    end
end

[maxFlow,F]=solveMaxFlowEdmondsKarp(CMatAugment,source,sink);

%Check whether all of the source nodes have been saturated. This is true if
%the total flow in equals the maximum flow found. A comparison to an
%epsilon value is used as it is not clear whether it can be guaranteed
%within finite precision bounds that the two quantities will always be
%equal for feasible problems, though they seem to generally be equal for
%feasible problems.
if(abs(maxFlow-totalFlowIn)>eps(totalFlowIn))
    exitCode=1;%The problem is infeasible.
    F=[];
    totalCost=[];
    return;
end

%Get rid of the artificial source and sink nodes from the flow matrix.
F=F(1:numVertex,1:numVertex);

%To make the computation of the residual matrix simpler, we will transform
%the cost matrix. This simplifies the computation of costs in the residual
%flow network.
AMatOrig=AMat;
AMat=AMat-AMat';

%Compute the residual flow network, also known as the incremental flow
%network. Definition 7.3 from Chapter 7.2 of [2] explains what the
%incremental flow network is and how the costs are assigned. The residual
%capacity is updated during augmentation of the flow matrix F, so that it
%need not be recomputed each loop.

%Find the capacity of the residual flow network.
residualCapacity=CMat-F;
for curIter=1:maxIter
    %Step 2: Find the minimum mean cost cycle with respect to the costs on
    %        the residual flow network, also known as the incremental flow
    %        network.

    %The costs associated with the residual flow network are AMat(i,j) on
    %forward nodes and -AMat(i,j) on backward nodes. Hence the reason why
    %we replaced AMat with AMat-AMat'.
    costMat=Inf(numVertex,numVertex);
    costMat(residualCapacity~=0)=AMat(residualCapacity~=0);

    [minCycleMean,cycleVertices]=findMinMaxCycleMean(costMat,true);

    if(minCycleMean>=0)
        exitCode=0;
        FP=F;
        FP(F<0)=0;
        totalCost=sum(sum(FP.*AMatOrig));
        return;
    end

    %Step 3: The minimum mean cycle is negative. We must cancel the minimum
    %        mean cycle. That is, we augment the graph along the negative
    %        cycle. We adjust the flow around the cycle by as much
    %        as possible without violating capacity constraints, so that
    %        the negative cycle no longer exists.

    %Determine the minimum residual capacity of the cycle.
    numVertexInCycle=length(cycleVertices);
    minCapacity=Inf;
    startVertex=cycleVertices(1);
    for curEndVertex=2:numVertexInCycle
        endVertex=cycleVertices(curEndVertex);

        curCapacity=residualCapacity(startVertex,endVertex);
        minCapacity=min(curCapacity,minCapacity);

        startVertex=endVertex;
    end

    %Push the amount of flow around the negative cycle that will saturate
    %the minimum capacity arc on the minimum mean cycle. Note that the skew
    %symmetry constraint dictates how the flow is to be added (how
    %augmentation is performed).
    startVertex=cycleVertices(1);
    for curEndVertex=2:numVertexInCycle
        endVertex=cycleVertices(curEndVertex);

        %Update the flow
        F(startVertex,endVertex)=F(startVertex,endVertex)+minCapacity;
        F(endVertex,startVertex)=F(endVertex,startVertex)-minCapacity;

        %Update the residual flow network. The update is opposite that of
        %F, because the residual flow is CMat-F.
        residualCapacity(startVertex,endVertex)=residualCapacity(startVertex,endVertex)-minCapacity;
        residualCapacity(endVertex,startVertex)=residualCapacity(endVertex,startVertex)+minCapacity;

        startVertex=endVertex;
    end
end

%Find the cost of the flow.
FP=F;
FP(F<0)=0;
totalCost=sum(sum(FP.*AMatOrig));

%If we are here, the maximum number of iterations has elapsed. We will
%check whether convergence occurred on the final iteration, and it not,
%then we will return accordingly.
costMat=Inf(numVertex,numVertex);
costMat(residualCapacity~=0)=AMat(residualCapacity~=0);

minCycleMean=findMinMaxCycleMean(costMat,true);
if(minCycleMean>=0)
    exitCode=0;
    return;
end

%If we get here, then the algorithm terminated without cancelling all
%negative cycles. Thus, the algorithm did not converge.
exitCode=2;
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.


text:
function y = solveLinearSystem(mode,A,b,para);

switch mode,
  case {'matlab','MATLAB'},
    y = A\b;
%     file = createFileName('prefix','JL');
%     save(file,'A','b');

  case 'MG',
    para.MGlevel      = log2(para.m(1))+1;
    para.MGcycle      = 1;
    para.MGomega      = 0.5;   %% !!! 0.5 should be better
    para.MGsmoother   = 'mfJacobi';
    para.MGpresmooth  = 3;
    para.MGpostsmooth = 1;
    para.dim = length(para.Omega);
    u = zeros(size(b));
    [y,res,r] = mfvcycle(para,u,b,1e-12,para.MGlevel,(max(para.m)>32));
end;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function test_bug1298

% MEM 3gb
% WALLTIME 00:10:00
% DEPENDENCY ft_timelockanalysis ft_prepare_leadfield ft_sourceanalysis

megraw = load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/beamformer/data_all.mat'));

cfg = [];
cfg.covariance = 'yes';
cfg.keeptrials = 'yes';
megtlock = ft_timelockanalysis(cfg,megraw.data_all);

load(dccnpath('/home/common/matlab/fieldtrip/data/test/latest/vol/Subject01vol_localspheres.mat'))

cfg = [];
cfg.headmodel = vol;
grid = ft_prepare_leadfield(cfg,megtlock);

cfg = [];
cfg.headmodel = vol;
cfg.method = 'lcmv';
cfg.sourcemodel = grid;
cfg.keepleadfield = 'yes';
cfg.lcmv.keepfilter = 'yes';
megsource1 = ft_sourceanalysis(cfg,megtlock);

cfg = [];
cfg.headmodel = vol;
cfg.method = 'lcmv';
cfg.sourcemodel = grid;
cfg.sourcemodel.leadfield = megsource1.leadfield;
cfg.sourcemodel.filter = megsource1.avg.filter;
% cfg.keeptrials = 'yes';
cfg.rawtrial = 'yes';
megsource11 = ft_sourceanalysis(cfg,megtlock);

% and then single trial estimates are in megsource11.trial

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% X = NGROUPSK(N, K)
%
% The number of ways N*K objects can be divided into N groups of equal size
% (K objects per group).
%
% N = number of groups
% K = size of the groups
%
% X = number of different group division combinations

% Last modified 2011-01-28
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function x = groupsk(n, k)

x = factorial(n.*k) ./ ( factorial(n) .* factorial(k).^n );

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function a = sec(a)
%SEC          Gradient secant sec(a)
%

% written  10/16/98     S.M. Rump
% modified 10/14/00     S.M. Rump  use Tony's trick
% modified 03/22/04     S.M. Rump  improved performance
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    accelaration for sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 10/08/08     S.M. Rump  improved sparse multiplication: not using intval data type
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  N = getappdata(0,'INTLAB_GRADIENT_NUMVAR');

  ax = tan(a.x(:));
  a.x = sec(a.x);
  ax = a.x(:) .* ax;
  if issparse(a.dx)
    sizeax = size(a.dx,1);
    [ia,ja,sa] = find(a.dx);
    % take care of scalar a.x: cures Matlab V6.0 bug
    % a=7; i=[1 1]; x=a(i), b=sparse(a); y=b(i)  yields row vector x but column vector y
    ax = ax(ia);
    if isa(a.x,'intval')
      adx = times(ax(:),sa(:),0);
      if adx.complex
        a.dx = intval( sparse(ia,ja,adx.mid,sizeax,N) , sparse(ia,ja,adx.rad,sizeax,N) , 'midrad' );
       else
        a.dx = intval( sparse(ia,ja,adx.inf,sizeax,N) , sparse(ia,ja,adx.sup,sizeax,N) , 'infsup' );
      end
    else
      a.dx = sparse(ia,ja,ax(:).*sa(:),sizeax,N);
    end
  else
    a.dx = a.dx .* ax(:,ones(1,N));
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function [params, runConfig] = readConfig(yaml_file)

configs = ReadYaml(yaml_file);
names = fieldnames(configs);
num = length(names);

for i=1:num
    if strcmp(names{i}, 'runConfigs')
        runConfig = getfield(configs, names{i});
        names = [names(1:i-1), names(i+1:end)];
        break;
    end
end

num = length(names);
param_bounds = zeros(num, 2);
discrete_indices = zeros(0, 0);
log_indices = zeros(0, 0);
categorical_indices = zeros(0, 0);
for i=1:num
    field = getfield(configs, names{i});
    param_bounds(i, 1) = field.lowerbound;
    param_bounds(i, 2) = field.upperbound;

    if isfield(field, 'integer') && field.integer
        discrete_indices(length(discrete_indices)+1) = i;
    end

    if isfield(field, 'categorical') && field.categorical
        categorical_indices(length(categorical_indices)+1) = i;
    end

    if isfield(field, 'logscale') && field.logscale
        log_indices(length(log_indices)+1) = i;
    end
end

param_bounds(log_indices, :) = log(param_bounds(log_indices, :));

params.high_dim = size(param_bounds, 1);
params.param_bounds = param_bounds';
params.log_indices = log_indices;
params.discrete_indices = discrete_indices;
params.categorical_indices = categorical_indices;
params.param_range = params.param_bounds(2, :) - params.param_bounds(1, :);
params.param_range(discrete_indices) = params.param_range(discrete_indices) + 1;
params.param_bounds(1, discrete_indices) = ...
    params.param_bounds(1, discrete_indices) - 0.5;

end

author: ziyuw

repo: rembo

save_path: github-repos/MATLAB/ziyuw-rembo


text:
function y=PSNR(noisyImage,restoredImage)

% Compute the PSNR of two gray scale image
% Traditional progarmming using loops
% Class input : [0,1]
% july, 25 , 2012
% KHMOU Youssef

N=size(noisyImage);
if length(N)> 2
    error('Input must be grayscale image');
end
if size(noisyImage)~=size(restoredImage)
    error('The images must have the same size');
end

%if ~isa(noisyImage,'double')
%   noisyImage=double(noisyImage)./255.00;
%end
%if  ~isa(restoredImage,'double')
%    restoredImage=double(restoredImage)./255.00;
%end

% begin

d1=max(noisyImage(:));
d2=max(restoredImage(:));
d=max(d1,d2);

MSE=0;
for i=1:N(1)
    for j=1:N(2)
        if isnan(noisyImage(i,j)) || isinf(restoredImage(i,j))...
                || isnan(restoredImage(i,j)) || isinf(noisyImage(i,j))
            continue;
        end
        MSE=MSE+((abs(noisyImage(i,j)-restoredImage(i,j))).^2);
    end
end

MSE=MSE./(N(1)*N(2));

y=10*log10((d.^2) /MSE)


text:
% 1) You need to check that your FlightGear installation path is configured correctly below.
% 2) You need to install the 'ExperimentalCarrier' aircraft in %FG_ROOT%/data/Aircraft, e.g. C:\Program Files\FlightGear 3.4.0\data\Aircraft\ExperimentalCarrier

setenv('FG_ROOT', 'C:\Program Files\FlightGear 3.4.0\data');
currPath = pwd;
cd 'C:\Program Files\FlightGear 3.4.0'
system('bin\fgfs --aircraft=ExperimentalCarrier --fdm=network,localhost,5501,5502,5503 --fog-fastest --disable-clouds --start-date-lat=2004:06:01:09:00:00 --disable-sound --in-air --enable-freeze --airport=KSFO --runway=10L --altitude=328 --heading=0 --offset-distance=0 --offset-azimuth=0 &')
cd(currPath);

author: jrgenerative

repo: fixed-wing-sim


text:
function qualityMeasure = callQualityMeasure(sr, gt, qualityMeasureName)

    % Assemble name of the quality measure wrapper function.
    funName = sprintf('imageQuality_%s', qualityMeasureName);

    % Call the desired quality measure.
    qualityMeasure = feval(funName, sr, gt);

author: thomas-koehler

repo: SupER


text:
function stsa_mis(filename,outfile)

%
%  Implements the Bayesian estimator based on the modified Itakura-Saito
%  distortion measure [1, Eq. 43].
%
%  Usage:  stsa_mis(noisyFile, outputFile)
%
%         infile - noisy speech file in .wav format
%         outputFile - enhanced output file in .wav format
%
%
%  Example call:  stsa_mis('sp04_babble_sn10.wav','out_mis.wav');
%
%  References:
%   [1] Loizou, P. (2005). Speech enhancement based on perceptually motivated
%       Bayesian estimators of the speech magnitude spectrum. IEEE Trans. on Speech
%       and Audio Processing, 13(5), 857-869.
%
% Author: Philipos C. Loizou
%
% Copyright (c) 2006 by Philipos C. Loizou
% $Revision: 0.0 $  $Date: 10/09/2006 $
%-------------------------------------------------------------------------

if nargin<2
   fprintf('Usage: stsa_mis inFile outFile.wav \n\n');
   return;
end

[x, Srate, bits]= wavread( filename);

% =============== Initialize variables ===============
%

len=floor(20*Srate/1000); % Frame size in samples
if rem(len,2)==1, len=len+1; end;
PERC=50; % window overlap in percent of frame size
len1=floor(len*PERC/100);
len2=len-len1;

win=hanning(len);  % define window
win = win*len2/sum(win);  % normalize window for equal level output

% Noise magnitude calculations - assuming that the first 6 frames is noise/silence
%
nFFT=len;
nFFT2=len/2;
noise_mean=zeros(nFFT,1);
j=1;
for k=1:5
   noise_mean=noise_mean+abs(fft(win.*x(j:j+len-1),nFFT));
   j=j+len;
end
noise_mu=noise_mean/5;
noise_mu2=noise_mu.^2;

%--- allocate memory and initialize various variables

img=sqrt(-1);
x_old=zeros(len1,1);
Nframes=floor(length(x)/len2)-1;
xfinal=zeros(Nframes*len2,1);

%===============================  Start Processing =======================================================
%
k=1;
aa=0.98;
fprintf('\nThis might take some time ...\n');
for n=1:Nframes

   insign=win.*x(k:k+len-1);

   %--- Take fourier transform of  frame ----

   spec=fft(insign,nFFT);
   sig=abs(spec); % compute the magnitude
   sig2=sig.^2;

       gammak=min(sig2./noise_mu2,40);  % post SNR. Limit it to avoid overflows
       if n==1
           ksi=aa+(1-aa)*max(gammak-1,0);
       else
           ksi=aa*Xk_prev./noise_mu2 + (1-aa)*max(gammak-1,0);     % a priori SNR
       end

       vk=ksi.*gammak./(1+ksi);

       sig_hat=log(comp_int(vk,gammak,sig)); % Eq. 41

       Xk_prev=sig_hat.^2;

       xi_w= ifft( sig_hat.* exp(img*angle(spec)));
	   xi_w= real( xi_w);

	% --- Overlap and add ---------------
    %
    xfinal(k:k+ len2-1)= x_old+ xi_w(1:len1);
	x_old= xi_w(len1+ 1: len);

    if rem(n,20)==0, fprintf('Frame: %d Percent completed:%4.2f\n',n,n*100/Nframes); end;

 k=k+len2;
end
%========================================================================================

wavwrite(xfinal,Srate,16,outfile);

%------------------------------E N D  -----------------------------------
function xhat=comp_int(vk,gammak,Yk)

% -- Evaluates Eq. 43 in [1]
%

Yk2=Yk.*Yk;
G2=gammak.^2;
EV=exp(-vk);

N=40; % number of terms to keep in infinite sum (Eq. 43)
L=length(vk)/2+1;
J1=zeros(L,1);
J2=zeros(L,1);

for j=1:L
  sum=0;  sum_b=0;
  for m=0:N
     F=factorial(m);
     d1=(vk(j))^m;
     d2=hyperg(-m,-m,0.5,Yk2(j)/(4*G2(j)),10);
     d2_b=hyperg(-m,-m,1.5,Yk2(j)/(4*G2(j)),10);
     sum=sum+d1*d2/F;
     sum_b=sum_b+gamma(m+1.5)*d1*d2_b/(F*gamma(m+1));
end
 J1(j)=sum;
 J2(j)=sum_b;
end

J1=J1.*EV(1:L);
J2=J2.*EV(1:L).*sqrt(vk(1:L)).*Yk(1:L)./gammak(1:L);

xhat2=max(real(J1+J2),0.00001);
xhat = [xhat2; flipud(xhat2(2:L-1))];

author: jtkim-kaist

repo: Speech-enhancement


text:
function [R,t,PfromTransformed]=rigidTransformation(Pfrom,Pto,varargin)
%% Find the optimal rigid body transformation between two point clouds %%%%%%%%%
% The optimal solution is the least square minimization solution between
% the two point clouds. See (Rubin & Solav 2016)
% The two point clouds must be corresponded (indeces refer to the same point)

%% INPUT:
% Pfrom=  point cloud (to be transformed into TruP)
% Pto= true 3D point cloud
% * optional input:
% estimation method: 'rigid' (default) or 'affine'

%% OUTPUT:
% R = estimated rotation matrix from Pfrom to Pto (3-by-3)
% t = estimated translation vector (3-by-1)
% PfromTransformed = Pfrom transformed using R and t

%%
nArg=numel(varargin);

switch nArg
    case 0
        estimationMethod='rigid';
    case 1
        estimationMethod=varargin{1};
end

%%
a=Pfrom;
b=Pto;

% find NaNs
aNanInd=find(isnan(a(:,1)));
bNanInd=find(isnan(b(:,1)));
abNanInd=[aNanInd; bNanInd];
% delete Nans
aNoNan=a;
aNoNan(abNanInd,:)=[];
bNoNan=b;
bNoNan(abNanInd,:)=[];

ac=(1/length(aNoNan))*sum(aNoNan); % centroid of  RecP.
bc=(1/length(bNoNan))*sum(bNoNan); % centroid of  TruP.

da=aNoNan-ac;
db=bNoNan-bc;

switch estimationMethod

    case 'rigid'
        M=db'*da;
        [U,~,V] = svd(M);
        S=[1 0 0; 0 1 0; 0 0 det(U*V')];
        R=U*S*V'; %rotation matrix from reconstructed coordinates to true
        t=bc'-R*ac'; %translation vector from centroid of reconstructed coordinates to true- this is the translation of the origin of the coordinate system!

    case 'affine'
        M=db'*da;
        A=da'*da;
        F=M/A;
%         R=F*(F'*F)^-.5; % this gives the same solution as the svd
        [U,~,V] = svd(F);
        S=[1 0 0; 0 1 0; 0 0 det(U*V')];
        R=U*S*V'; %rotation matrix from reconstructed coordinates to true
        t=bc'-R*ac'; %translation vector from centroid of reconstructed coordinates to true- this is the translation of the origin of the coordinate system!

end

PfromTransformed=(R*a'+t)'; %transformed constructed coordinates

% residuals=PfromTransformed-Pto;
% residualsMgn=sqrt(sum(residuals.^2,2));
% residulasMgnRMS=sqrt(nansum(residualsMgn.^2)/sum(~isnan(residualsMgn)));

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function [xhk, pf] = particle_filter(sys, yk, pf, resampling_strategy)
%% Generic particle filter
%
% Note: when resampling is performed on each step this algorithm is called
% the Bootstrap particle filter
%
% Usage:
% [xhk, pf] = particle_filter(sys, yk, pf, resamping_strategy)
%
% Inputs:
% sys  = function handle to process equation
% yk   = observation vector at time k (column vector)
% pf   = structure with the following fields
%   .k                = iteration number
%   .Ns               = number of particles
%   .w                = weights   (Ns x T)
%   .particles        = particles (nx x Ns x T)
%   .gen_x0           = function handle of a procedure that samples from the initial pdf p_x0
%   .p_yk_given_xk    = function handle of the observation likelihood PDF p(y[k] | x[k])
%   .gen_sys_noise    = function handle of a procedure that generates system noise
% resampling_strategy = resampling strategy. Set it either to
%                       'multinomial_resampling' or 'systematic_resampling'
%
% Outputs:
% xhk   = estimated state
% pf    = the same structure as in the input but updated at iteration k
%
% Reference:
% [1] Arulampalam et. al. (2002).  A tutorial on particle filters for
%     online nonlinear/non-gaussian bayesian tracking. IEEE Transactions on
%     Signal Processing. 50 (2). p 174--188

%% Programmed by:
% Diego Andres Alvarez Marin (diegotorquemada@gmail.com)
% Universidad Nacional de Colombia at Manizales, February 29, 2012

%%
k = pf.k;
if k == 1
   error('error: k must be an integer greater or equal than 2');
end

%% Initialize variables
Ns = pf.Ns;                              % number of particles
nx = size(pf.particles,1);               % number of states

wkm1 = pf.w(:, k-1);                     % weights of last iteration
if k == 2
   for i = 1:Ns                          % simulate initial particles
      pf.particles(:,i,1) = pf.gen_x0(); % at time k=1
   end
   wkm1 = repmat(1/Ns, Ns, 1);           % all particles have the same weight
end

%%
% The importance sampling function:
% PRIOR: (this method is sensitive to outliers)   THIS IS THE ONE USED HERE
% q_xk_given_xkm1_yk = pf.p_xk_given_xkm1;

% OPTIMAL:
% q_xk_given_xkm1_yk = q_xk_given_xkm1^i_yk;
% Note this PDF can be approximated by MCMC methods: they are expensive but
% they may be useful when non-iterative schemes fail

%% Separate memory
xkm1 = pf.particles(:,:,k-1); % extract particles from last iteration;
xk   = zeros(size(xkm1));     % = zeros(nx,Ns);
wk   = zeros(size(wkm1));     % = zeros(Ns,1);

%% Algorithm 3 of Ref [1]
for i = 1:Ns
   % xk(:,i) = sample_vector_from q_xk_given_xkm1_yk given xkm1(:,i) and yk
   % Using the PRIOR PDF: pf.p_xk_given_xkm1: eq 62, Ref 1.
   xk(:,i) = sys(k, xkm1(:,i), pf.gen_sys_noise());

   % Equation 48, Ref 1.
   % wk(i) = wkm1(i) * p_yk_given_xk(yk, xk(:,i))*p_xk_given_xkm1(xk(:,i), xkm1(:,i))/q_xk_given_xkm1_yk(xk(:,i), xkm1(:,i), yk);

   % weights (when using the PRIOR pdf): eq 63, Ref 1
   wk(i) = wkm1(i) * pf.p_yk_given_xk(k, yk, xk(:,i));

   % weights (when using the OPTIMAL pdf): eq 53, Ref 1
   % wk(i) = wkm1(i) * p_yk_given_xkm1(yk, xkm1(:,i)); % we do not know this PDF
end;

%% Normalize weight vector
wk = wk./sum(wk);

%% Calculate effective sample size: eq 48, Ref 1
Neff = 1/sum(wk.^2);

%% Resampling
% remove this condition and sample on each iteration:
% [xk, wk] = resample(xk, wk, resampling_strategy);
%if you want to implement the bootstrap particle filter
resample_percentaje = 0.50;
Nt = resample_percentaje*Ns;
if Neff < Nt
   disp('Resampling ...')
   [xk, wk] = resample(xk, wk, resampling_strategy);
   % {xk, wk} is an approximate discrete representation of p(x_k | y_{1:k})
end

%% Compute estimated state
xhk = zeros(nx,1);
for i = 1:Ns;
   xhk = xhk + wk(i)*xk(:,i);
end

%% Store new weights and particles
pf.w(:,k) = wk;
pf.particles(:,:,k) = xk;

return; % bye, bye!!!

%% Resampling function
function [xk, wk, idx] = resample(xk, wk, resampling_strategy)

Ns = length(wk);  % Ns = number of particles

% wk = wk./sum(wk); % normalize weight vector (already done)

switch resampling_strategy
   case 'multinomial_resampling'
      with_replacement = true;
      idx = randsample(1:Ns, Ns, with_replacement, wk);
%{
      THIS IS EQUIVALENT TO:
      edges = min([0 cumsum(wk)'],1); % protect against accumulated round-off
      edges(end) = 1;                 % get the upper edge exact
      % this works like the inverse of the empirical distribution and returns
      % the interval where the sample is to be found
      [~, idx] = histc(sort(rand(Ns,1)), edges);
%}
   case 'systematic_resampling'
      % this is performing latin hypercube sampling on wk
      edges = min([0 cumsum(wk)'],1); % protect against accumulated round-off
      edges(end) = 1;                 % get the upper edge exact
      u1 = rand/Ns;
      % this works like the inverse of the empirical distribution and returns
      % the interval where the sample is to be found
      [~, idx] = histc(u1:1/Ns:1, edges);
   % case 'regularized_pf'      TO BE IMPLEMENTED
   % case 'stratified_sampling' TO BE IMPLEMENTED
   % case 'residual_sampling'   TO BE IMPLEMENTED
   otherwise
      error('Resampling strategy not implemented')
end;

xk = xk(:,idx);                    % extract new particles
wk = repmat(1/Ns, 1, Ns);          % now all particles have the same weight

return;  % bye, bye!!!

author: Sable

repo: mcbench-benchmarks


text:
function value = year_length_months_republican ( y )

%*****************************************************************************80
%
%% YEAR_LENGTH_MONTHS_REPUBLICAN returns the number of months in a Republican year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, integer VALUE, the number of months in the year.
%
  value = 13;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Reads ECG recording in HES format
% Reads ECG recordings in HES (Biosigna) format. Implements the documentation
% available in the help of the application provided with the database (not
% available with the ECHkit).
%
% Arguments:
%   + filename: recording to be read.
%   + start_sample: (opt) start sample to read. Default 1.
%   + end_sample: (opt) end sample to read. Default min(All recording, ECG block of 200 Mbytes)
%
% Output:
%   + ECG: the ECG block
%   + heasig: header with the ECG properties.
%   + ann: annotations for the ECG recordings.
%
% Limits:
% This routine is limited to read blocks smaller than 200 Mbytes for
% performance reasons. You can disable this limit by doing:
% MaxIOread = Inf; %megabytes
%
% See also read_HES_ann, read_HES_header, read_ECG, ECGwrapper
%
% Author: Mariano Llamedo Soria
% <matlab:web('mailto:llamedom@electron.frba.utn.edu.ar','-browser') (email)>
% Version: 0.1 beta
% Birthdate: 17/12/2010
% Last update: 19/11/2014
% Copyright 2008-2015
%
function [ECG heasig ann last_sample] = read_HES_format( filename, start_sample, end_sample )

ann = [];
heasig = [];
ECG = [];
last_sample = [];

tablas_y_constantes;

%No leer bloques mas grandes de 200 megabytes
MaxIOread = 200; %megabytes

if( nargin < 2 || isempty( start_sample ) )
    start_sample = 1;
else
    start_sample = max(1,start_sample);
end

fidECG = fopen( filename, 'r');

if( fidECG > 0 )

    fseek(fidECG, 900, 'bof');

    num_of_leads = fread(fidECG, 1, 'uint8');
    num_of_leads_simult = fread(fidECG, 1, 'uint8');

    if( num_of_leads ~= num_of_leads_simult)
        error('Sampleo no simult�neo, revisar.')
    end

    heasig.nsig = num_of_leads_simult;

    lead_description_idx = fread(fidECG, heasig.nsig, 'uint8');
    [dummy lead_description_table_idx] = intersect(Lead_description_idx, lead_description_idx);
    heasig.desc = char(cLead_description_table(lead_description_table_idx,1));

    fseek(fidECG, 992, 'bof');

    aux = fread(fidECG, 5, 'uint16');

    sample_interval_usec = aux(1);
    heasig.freq = 1/sample_interval_usec*1e6;
    heasig.gain = repmat(1/aux(3),1,heasig.nsig);
    heasig.units = repmat('nV',heasig.nsig,1);
    heasig.adcres = repmat(aux(4),1,heasig.nsig);
    heasig.adczero = repmat(2^(aux(4)-1),1,heasig.nsig);
    heasig.nsamp = fread(fidECG, 1, 'uint32');

    if( nargin < 3 || isempty( end_sample ) )
        %Intento la lectura total por defecto
        samples2read = heasig.nsamp - (start_sample-1);
    else
        samples2read = min(heasig.nsamp, end_sample) - (start_sample-1);
    end

    if( (samples2read*heasig.nsig*2) > (MaxIOread * 1024^2) )
        samples2read = (MaxIOread * 1024^2) / heasig.nsig / 2;
        warning(['No es recomendable leer mas de ' num2str(MaxIOread) ' Mb. Realice varias lecturas.'])
    end

    ECG = nan(samples2read,heasig.nsig);

    try

        fseek(fidECG, 1024+((start_sample-1)*heasig.nsig)*2, 'bof');

        ECG = fread(fidECG, [heasig.nsig samples2read], '*int16')';

        fclose(fidECG);

    catch ME
        fclose(fidECG);
        rethrow(ME)
    end

    last_sample = size(ECG,1) + start_sample - 1;

    [~, heasig.recname] = fileparts(filename);
    heasig.btime = '00:00:00';
    heasig.bdate = '01/01/2000';

    if( nargout > 2 )
        ann = read_HES_ann([ filename(1:end-4) '.lst' ]);
        ann.time = round(ann.time * heasig.freq);
    end

end

author: marianux

repo: ecg-kit

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


text:
function [traces, rel, cSize] = calcTraces(grains, clusterId, varargin)
% traces of subsets of grains
%
% Syntax
%
%   clusterId = [grainId,variantId];
%   [traces, rel, cSize] = calcTraces(grains, clusterId,'shape')
%
% Input
%  grains - @grain2d
%  cId  - cluster Id
%
% Output
%  traces - @vector3d, size max(clusterId(:,1)) x max(clusterId(:,2)) ....
%  rel    - relyability index, same size as traces
%  cSize  - cluster size, , same size as traces
%
% Options
%  minClusterSize - minimum grainSize required for trace computation (default: 100)
%  shape - characteristic shape based algorithm
%  calliper - use shortest calliper instead of eigenvectors
%  hist  - circular histogram based algorithm
%
% References
%
% * <https://arxiv.org/abs/2201.02103 Determination of child phase habit
% planes from two-dimensional reconstructed parent phase orientation maps>,
% arXiv, 2022
%

% generate clusters
if nargin == 1 || ~isnumeric(clusterId), clusterId = ones(length(grains),1); end
notZero = all(clusterId>0,2);
cSize = accumarray(clusterId(notZero,:),grains.grainSize(notZero));
ic = find(cSize > get_option(varargin,'minClusterSize',100));

% prepare output
sz = [max(clusterId),1];
omega = nan(sz);
rel = zeros(sz);

% extract grain geometry
V = grains.boundary.V;
F = grains.boundary.F;
I_BG = grains.boundary.I_FG;
grainId = grains.id;

% get method to use
useHist = check_option(varargin,'hist');
useCalliper = check_option(varargin,'calliper');

% loop over all relevant clusters
for i = 1:length(ic)

  % combine grains according to clusterIndex
  [sub{1:size(clusterId,2)}]  = ind2sub(sz,ic(i));
  ind = all(clusterId == [sub{:}],2);

  % the boundary segments
  lF = F(any(I_BG(:,grainId(ind)),2),:);

  % shifted into the origin
  lS = V(lF(:,1),:) - V(lF(:,2),:);

  if isempty(lF), continue; end

  if useHist % density function method

    rho = atan2(lS(:,2),lS(:,1));

    %fun = calcDensity(rho, 'weights',sqrt(sum(lS.^2,2)),'periodic','sigma',10*degree);
    fun = calcDensity(rho, 'weights',vecnorm(lS,2,2),'periodic','sigma',10*degree);
    fun.antipodal = true;

    phi = linspace(0,2*pi,360);
    [m,ind] = max(real(fun.eval(phi)));

    omega(ic(i)) = phi(ind);
    rel(ic(i)) = (m-1)/m;

  else % characteristic shape method

    %cS = shape2d.byFV(F(any(I_BG(:,ind),2),:),V,'noSimplify');
    %[omega(ic(i)),a,b] = principalComponents(cS);
    %traces(ic(i)) = cS.caliper('shortestPerp'); % this is a bit more precise but slower

    % the following lines are from shape2d.byFV
    % just consider one direction
    fcond = lS(:,2)<0;
    lS(fcond,:)=lS(fcond,:).*-1;
    dxy = [lS; -lS];

    % sort segments according to angle
    [~,id]= sort(atan2(dxy(:,2),dxy(:,1)));
    dxy = dxy(id,:);

    % sum up
    xyn = cumsum(dxy);

    % shift again
    xyn = [xyn(:,1) - mean(xyn(:,1)) xyn(:,2) - mean(xyn(:,2))];

    if useCalliper

      mid = round(size(xyn,1)/2);
      dxyn = xyn - [xyn(1+mid:end,:);xyn(1:mid,:)];
      delta = vecnorm(dxyn,2,2);

      % minimum and maximum Ferret diameter
      a = max(delta); [b,ib] = min(delta);

      % use vector perpendicular to short axis
      omega(ic(i)) = pi/2 + atan2(dxyn(ib,2),dxyn(ib,1));

    else

      % the following lines are taken from grain2d/principleComponent
      % compute length of line segments
      dist = sqrt(sum((xyn(1:end-1,:) - xyn(2:end,:)).^2,2));
      dist = 0.5*(dist(1:end) + [dist(end);dist(1:end-1)]);

      % weight vertices according to half the length of the adjacent faces
      xyn = xyn(1:end-1,:) .* [dist,dist] .* sum(xyn(1:end-1,:).^2,2).^(0.25);

      % compute eigen values and vectors
      [ew, omega(ic(i))] = eig2(xyn' * xyn);

      % halfaxes are square roots of the eigenvalues
      b = sqrt(ew(1)); a = sqrt(ew(2));
    end
    rel(ic(i)) = (a-b)./a;
  end

end

traces = vector3d.byPolar(pi/2,omega,'antipodal');

author: mtex-toolbox

repo: mtex

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


text:
function [W, TE, R, R2, SDS] = ROLLINGSTYLE(F, S, N)
% Rolling style analysis
%
% INPUTS:
%   F... vector of fund returns
%   S... matrix of style factor returns
%   N... number of rolling periods
%
% OUTPUTS:
% W... vector of optimal style index weights
% TE... tracking error between calculated and actual fund
% Fcalc... vector of calculated fund time series
% R2... coefficient of determination between fund and calculated fund
% SDS... Style Drift Score
%
% Andreas Steiner
% performanceanalysis@andreassteiner.net,
% http://www.andreassteiner.net/performanceanalysis

W = []; TE = []; R = []; R2 = [];

for t = N:rows(F)

     [w, te, r, r2] = TEMINQP(F(t-(N-1):t), S(t-(N-1):t,:));

     W = [W; w];
     TE = [TE; te];
     R = [R; r(end)];
     R2 = [R2; r2];

end;

SDS = sqrt(sum(var(W)));

figure
bar(W, 'stacked')
title('Rolling Style Analysis: Style Weights Over Time')
axis tight


text:
function [points2d,z3] = project3dPtsTo2d(points3d, K, Rtilt, trans, crop)
% project 3D point cloud in world coordinate to x-y 2D plane. Provided by
% Shuran Song
% points3d: point cloud in camera coordinate
% K: camera intrinsic
% Rtilt, trans: camera extrinsic
% crop: assume it be 1

    %% inverse of get_aligned_point_cloud
    points3d = bsxfun(@plus, Rtilt * points3d', trans)';

    %% inverse rgb_plane2rgb_world
    if isempty(K)
        camera_params;
    else
        cx_rgb = K(1,3); cy_rgb = K(2,3);
        fx_rgb = K(1,1); fy_rgb = K(2,2);
    end
    % Make the original consistent with the camera location:
    x3 = points3d(:,1);
    y3 = -points3d(:,2); % when doing projection or depth->3d, always flip y dimension.
    z3 = points3d(:,3);

    xx = x3 * fx_rgb ./ z3 + cx_rgb;
    yy = y3 * fy_rgb ./ z3 + cy_rgb;

    if ~exist('crop','var')||isempty(crop)
        xx = xx - 41 + 1;
        yy = yy - 45 + 1;
    else
        xx = xx - crop(2) + 1;
        yy = yy - crop(1) + 1;
    end
    points2d = [xx yy];
end

author: zhirongw