text:
function model = destroyDummyModel(model, dummyMetBool, dummyRxnBool, metsOrig, rxnsOrig)
% Remove all traces of dummy metabolites and reactions
%
% INPUT
% model     COBRA model
%
% OPTIONAL INPUT
% dummyMetBool      m x 1 boolean vector indicating dummy metabolites
% dummyRxnBool      n x 1 boolean vector indicating dummy reactions
% metsOrig          Array indicating the metabolites before creating a dummyModel
% rxnsOrig          Array indicating the reactions before creating dummyModel
%
% OUTPUT
% model     COBRA model without dummy metabolites or reactions
%

if nargin < 2 || isempty(dummyMetBool)
    % remove dummy metabolites
    dummyMetBool = contains(model.mets, 'dummy_Met_');
end
if nargin < 3 || isempty(dummyRxnBool)
    % Remove dummy reactions
    dummyRxnBool = contains(model.rxns, 'dummy_Rxn_');
end
% if nargin < 4 || isempty(metsOrig)
%     % Original metabolites
%     metsOrig = model.mets;
% end
% if nargin < 5 || isempty(rxnsOrig)
%     % Original reactions
%     rxnsOrig = model.rxns;
% end

if any(dummyRxnBool) || any(dummyMetBool)
    model.mets = model.mets(~dummyMetBool);
    model.S = model.S(~dummyMetBool, ~dummyRxnBool);
    model.b = model.b(~dummyMetBool);
    model.rxns = model.rxns(~dummyRxnBool);
    model.rxnNames = model.rxnNames(~dummyRxnBool);
    model.lb = model.lb(~dummyRxnBool);
    model.ub = model.ub(~dummyRxnBool);
    model.c = model.c(~dummyRxnBool);

    if isfield(model, 'csense')
        model.csense = model.csense(~dummyMetBool);
    end
    if isfield(model, 'ctrs')
        model.C = model.C(:, ~dummyRxnBool);
    end
    if isfield(model, 'rxnGeneMat')
        model.rxnGeneMat = model.rxnGeneMat(~dummyRxnBool, :);
    end
    if isfield(model, 'rules')
        model.rules = model.rules(~dummyRxnBool);
    end
    if isfield(model, 'subSystems')
        model.subSystems = model.subSystems(~dummyRxnBool);
    end
    if isfield(model, 'SIntMetBool')
        model.SIntMetBool = model.SIntMetBool(~dummyMetBool);
    end
    if isfield(model, 'SIntRxnBool')
        model.SIntRxnBool = model.SIntRxnBool(~dummyRxnBool);
    end
    if isfield(model, 'SConsistentMetBool')
        model.SConsistentMetBool = model.SConsistentMetBool(~dummyMetBool);
    end
    if isfield(model, 'SConsistentRxnBool')
        model.SConsistentRxnBool = model.SConsistentRxnBool(~dummyRxnBool);
    end
    if isfield(model, 'fluxConsistentMetBool')
        model.fluxConsistentMetBool = model.fluxConsistentMetBool(~dummyMetBool);
    end
    if isfield(model, 'fluxConsistentRxnBool')
        model.fluxConsistentRxnBool = model.fluxConsistentRxnBool(~dummyRxnBool);
    end
    if isfield(model, 'thermoFluxConsistentMetBool')
        model.thermoFluxConsistentMetBool = model.thermoFluxConsistentMetBool(~dummyMetBool);
    end
    if isfield(model, 'thermoFluxConsistentRxnBool')
        model.thermoFluxConsistentRxnBool = model.thermoFluxConsistentRxnBool(~dummyRxnBool);
    end
end

if isfield(model,'dummyMetBool')
    model = rmfield(model, 'dummyMetBool');
end
if isfield(model,'dummyRxnBool')
    model = rmfield(model, 'dummyRxnBool');
end

% model.S = model.S(~model.dummyMetBool,~model.dummyRxnBool);
% model.b = model.b(~model.dummyMetBool,1);
% model.csense = model.csense(~model.dummyMetBool,1);
% model.mets = model.mets(~model.dummyMetBool);
% model.rxns = model.rxns(~model.dummyRxnBool);
% model.rxnNames = model.rxnNames(~model.dummyRxnBool);
% model.lb = model.lb(~model.dummyRxnBool);
% model.ub = model.ub(~model.dummyRxnBool);
% model.c = model.c(~model.dummyRxnBool);
% if isfield(model, 'ctrs')
%     model.C = model.C(:,~model.dummyRxnBool);
% end
% if isfield(model, 'rxnGeneMat')
%     model.rxnGeneMat = model.rxnGeneMat(~model.dummyRxnBool,:);
% end
% if isfield(model, 'rules')
%     model.rules = model.rules(~model.dummyRxnBool);
% end
% if isfield(model,'subSystems')
%     model.subSystems = model.subSystems(~model.dummyRxnBool);
% end
% if isfield(model,'SIntRxnBool')
%     model.SIntRxnBool = model.SIntRxnBool(~model.dummyRxnBool);
% end
% if isfield(model,'SConsistentRxnBool')
%     model.SConsistentRxnBool = model.SConsistentRxnBool(~model.dummyRxnBool);
% end
% if isfield(model,'fluxConsistentRxnBool')
%     model.fluxConsistentRxnBool = model.fluxConsistentRxnBool(~model.dummyRxnBool);
% end
% if isfield(model,'thermoFluxConsistentRxnBool')
%     model.thermoFluxConsistentRxnBool = model.thermoFluxConsistentRxnBool(~model.dummyRxnBool);
% end

if exist('metsOrig','var')
    % Remove deleted mets/rxns in all the fields
    if ~isequal(model.mets, metsOrig) || ~isequal(model.rxns, rxnsOrig)
        warning('mets/rxns were not consistently deleted in all fields!')
        cbFields = fieldnames(model);
        mets2Delete = ~ismember(metsOrig, model.mets);
        rxns2Delete = ~ismember(rxnsOrig, model.rxns);
        for i = 1:length(cbFields)
            if length(metsOrig) == length(model.(cbFields{i}))
                model.(cbFields{i})(mets2Delete) = [];
            elseif length(rxnsOrig) == length(model.(cbFields{i}))
                model.(cbFields{i})(rxns2Delete) = [];
            end
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function [TorF, vstr, rdate] = have_feature_opti_clp()
%HAVE_FEATURE_OPTI_CLP  Detect availability/version info for OPTI_CLP
%
%   Feature detection function implementing 'opti_clp' tag for HAVE_FEATURE
%   to detect availability/version of the version of CLP (COIN-OR Linear
%   Programming solver) distributed with OPTI Toolbox
%   (https://www.inverseproblem.co.nz/OPTI/).
%
%   See also HAVE_FEATURE, HAVE_FEATURE_CLP, QPS_MASTER, CLP.

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

TorF = exist('opti_clp', 'file') == 2 && exist('clp', 'file') == 3;
vstr = '';
rdate = '';
if TorF
    str = evalc('clp');
    pat = 'CLP: COIN-OR Linear Programming \[v([^\s,]+), Built ([^\],])+(,[^\]]*)*\]';  %% OPTI, Giorgetti/Currie
    [s,e,tE,m,t] = regexp(str, pat);
    if ~isempty(t)
        vstr = t{1}{1};
        rdate = datestr(t{1}{2}, 'dd-mmm-yyyy');
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
classdef ( Hidden, Sealed ) PointerManager < handle
    %uix.PointerManager  Pointer manager

    %  Copyright 2016 The MathWorks, Inc.
    %  $Revision: 1435 $ $Date: 2016-11-17 17:50:34 +0000 (Thu, 17 Nov 2016) $

    properties( SetAccess = private )
        Figure % figure
    end

    properties( Access = private )
        Tokens % tokens
        Pointers % pointers
        NextToken % next token
        PointerListener % listener
    end

    methods( Access = private )

        function obj = PointerManager( figure )
            %uix.PointerManager  Create pointer manager
            %
            %  m = uix.PointerManager(f) creates a pointer manager for the
            %  figure f.

            obj.Figure = figure;
            obj.Tokens = 0;
            obj.Pointers = {figure.Pointer};
            obj.NextToken = 1;
            obj.PointerListener = event.proplistener( figure, ...
                findprop( figure, 'Pointer' ), 'PostSet', ...
                @obj.onPointerChanged );

        end % constructor

    end % structors

    methods( Access = private )

        function doSetPointer( obj, token, pointer )
            %doSetPointer  Set pointer
            %
            %  m.doSetPointer(t,p) sets the pointer to p with the token t.

            % Remove old entry
            tf = obj.Tokens == token;
            obj.Tokens(tf) = [];
            obj.Pointers(tf) = [];

            % Add new entry
            obj.Tokens(end+1) = token;
            obj.Pointers{end+1} = pointer;

            % Set pointer
            obj.PointerListener.Enabled = false;
            obj.Figure.Pointer = pointer;
            obj.PointerListener.Enabled = true;

        end % doSetPointer

        function doUnsetPointer( obj, token )
            %doUnsetPointer  Unset pointer
            %
            %  m.doUnsetPointer(s) unsets the pointer with the token t.

            % Remove old entry
            tf = obj.Tokens == token;
            obj.Tokens(tf) = [];
            obj.Pointers(tf) = [];

            % Update pointer
            obj.PointerListener.Enabled = false;
            obj.Figure.Pointer = obj.Pointers{end};
            obj.PointerListener.Enabled = true;

        end % doUnsetPointer

    end % private methods

    methods

        function onPointerChanged( obj, ~, ~ )
            %onPointerChanged  Event handler

            % Log as unknown setter
            obj.doSetPointer( 0, obj.Figure.Pointer )

        end % onPointerChanged

    end % event handlers

    methods( Static )

        function token = setPointer( figure, pointer )
            %setPointer  Set pointer
            %
            %  t = uix.PointerManager.setPointer(f,p) sets the pointer of
            %  the figure f to p.  The returned token t can be used
            %  subsequently to unset the pointer.

            % Get pointer manager
            obj = uix.PointerManager.getInstance( figure );

            % Retrieve token
            token = obj.NextToken;

            % Set
            obj.doSetPointer( token, pointer )

            % Increment token
            obj.NextToken = token + 1;

        end % setPointer

        function unsetPointer( figure, token )
            %unsetPointer  Unset pointer
            %
            %  uix.PointerManager.unsetPointer(f,t) unsets the pointer of
            %  the figure f using the token t.

            % Check ID
            validateattributes( token, {'numeric'}, {'scalar','integer','>',0} )

            % Get pointer manager
            obj = uix.PointerManager.getInstance( figure );

            % Unset
            obj.doUnsetPointer( token )

        end % unsetPointer

        function obj = getInstance( figure )
            %getInstance  Get pointer manager
            %
            %  m = uix.PointerManager.getInstance(f) gets the pointer
            %  manager for the figure f.

            % Get pointer manager
            name = 'UIxPointerManager';
            if isprop( figure, name ) % existing, retrieve
                obj = figure.( name );
            else % new, create and store
                obj = uix.PointerManager( figure );
                p = addprop( figure, name );
                p.Hidden = true;
                figure.( name ) = obj;
            end

        end % getInstance

    end % static methods

end % classdef

author: goGPS-Project


text:
%% EVENT ENUMERATION CLASSES (eventType.m) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file contains the enumeration definition of the available events
% that can happen inside the simulation environment.

% Author: James A. Douthwaite 05/10/16

classdef eventType < uint8
   % Assign uint8 type to enumeration values
   enumeration
      event(0)          % Generic event
      detection(1)      % Detection
      warning(2)        % Proximity warning
      collision(3)      % Collision
      waypoint(4)       % Waypoint achieve
      null_detection(5) % Detection loss
      null_warning(6)   % Warning-clear
      null_collision(7) % Collision-clear
      null_waypoint(8)  % Waypoint drop
   end
end

repo: OpenMAS


text:
classdef vfa_t1 < AbstractModel
% vfa_t1: Compute a T1 map using Variable Flip Angle
%
% Assumptions:
%
% Inputs:
%   VFAData         Spoiled Gradient echo data, 4D volume with different flip angles in time dimension
%   (B1map)         Normalized transmit excitation field map (B1+). B1+ is defined
%                   as a  normalized multiplicative factor such that:
%                   FA_actual = B1+ * FA_nominal. (OPTIONAL).
%   (Mask)          Binary mask to accelerate the fitting. (OPTIONAL)
%
% Outputs:
%   T1              Longitudinal relaxation time [s]
%   M0              Equilibrium magnetization
%
% Protocol:
%   VFAData Array [nbFA x 2]:
%       [FA1 TR1; FA2 TR2;...]      flip angle [degrees] TR [s]
%
% Options:
%   None
%
% Example of command line usage:
%   Model = vfa_t1;  % Create class from model
%   Model.Prot.VFAData.Mat=[3 0.015; 20 0.015]; %Protocol: 2 different FAs
%   data = struct;  % Create data structure
%   data.VFAData = load_nii_data('VFAData.nii.gz');
%   data.B1map = load_nii_data('B1map.nii.gz');
%   FitResults = FitData(data,Model); %fit data
%   FitResultsSave_mat(FitResults);
%
%   For more examples: qMRusage(vfa_t1)
%
%
% Author: Ian Gagnon, 2017
%
% References:
%   Please cite the following if you use this module:
%     Fram, E.K., Herfkens, R.J., Johnson, G.A., Glover, G.H., Karis, J.P.,
%     Shimakawa, A., Perkins, T.G., Pelc, N.J., 1987. Rapid calculation of
%     T1 using variable flip angle gradient refocused imaging. Magn. Reson.
%     Imaging 5, 201?208
%   In addition to citing the package:
%     Karakuzu A., Boudreau M., Duval T.,Boshkovski T., Leppert I.R., Cabana J.F.,
%     Gagnon I., Beliveau P., Pike G.B., Cohen-Adad J., Stikov N. (2020), qMRLab:
%     Quantitative MRI analysis, under one umbrella doi: 10.21105/joss.02343

properties (Hidden=true)
 onlineData_url = 'https://osf.io/7wcvh/download?version=3';
end

    properties
        MRIinputs = {'VFAData','B1map','Mask'};
        xnames = {'M0','T1'};
        voxelwise = 0;

        % Protocol
        Prot  = struct('VFAData',struct('Format',{{'FlipAngle' 'TR'}},...
                                         'Mat', [3 0.015; 20 0.015])); % You can define a default protocol here.

        % fitting options
        st           = [2000 0.7]; % starting point
        lb           = [0   0.00001]; % lower bound
        ub           = [6000   5]; % upper bound
        fx           = [0     0]; % fix parameters

        % Model options
        buttons = {};
        options= struct(); % structure filled by the buttons. Leave empty in the code

        % Simulation Options
        Sim_Single_Voxel_Curve_buttons = {'SNR',50};
        Sim_Optimize_Protocol_buttons = {'# of volumes',5,'Population size',100,'# of migrations',100};
    end

methods (Hidden=true)
% Hidden methods goes here.
end

    methods

        function obj = vfa_t1()
            obj.options = button2opts(obj.buttons);
        end

        function Smodel = equation(obj,x)
            % Generates a VFA signal based on input parameters
            x = mat2struct(x,obj.xnames); % if x is a structure, convert to vector

            % Equation: S=M0sin(a)*(1-E)/(1-E)cos(a); E=exp(-TR/T1)
            flipAngles = (obj.Prot.VFAData.Mat(:,1))';
            TR = obj.Prot.VFAData.Mat(1,2);
            E = exp(-TR/x.T1);
            Smodel = x.M0*sin(flipAngles/180*pi)*(1-E)./(1-E*cos(flipAngles/180*pi));
        end

       function FitResult = fit(obj,data)
            % T1 and M0
            flipAngles = (obj.Prot.VFAData.Mat(:,1))';
            TR = obj.Prot.VFAData.Mat(:,2);
            if obj.voxelwise == 0
                if (length(unique(TR))~=1), error('VFA data must have same TR'); end
                if ~isfield(data, 'B1map'), data.B1map = []; end
                if ~isfield(data, 'Mask'), data.Mask = []; end
                [FitResult.T1, FitResult.M0] = Compute_M0_T1_OnSPGR(double(data.VFAData), flipAngles, TR(1), data.B1map, data.Mask);
            elseif obj.voxelwise == 1
                if ~isfield(data,'B1map'), data.B1map=1; end
                [m0, t1] = mtv_compute_m0_t1(double(data.VFAData), flipAngles, TR(1), data.B1map);
                FitResult.T1 = t1;
                FitResult.M0 = m0;
            end
       end

       function plotModel(obj,x,data)
           if nargin<2 || isempty(x), x = obj.st; end
           x = mat2struct(x,obj.xnames);
           disp(x)
           flipAngles = obj.Prot.VFAData.Mat(:,1)';
           TR = obj.Prot.VFAData.Mat(1,2)';
           subplot(2,1,1)
           if exist('data','var')
               if isfield(data,'B1map')
                  if ~isempty(data.B1map)
                   B1map=data.B1map;
                  end
               else
                   B1map=1;
               end
                % Plot data and fitted signal
                plot(flipAngles,data.VFAData,'.','MarkerSize',16)
            else
                B1map=1;
            end
            E = exp(-TR/x.T1);
            Smodel = x.M0*sin(flipAngles/180*pi*B1map)*(1-E)./(1-E*cos(flipAngles/180*pi*B1map));
            hold on
            plot(flipAngles,Smodel,'x','MarkerSize',16)
            hold off
            title('Data points','FontSize',14);
            xlabel('Flip Angle [deg]','FontSize',12);
            ylabel('Signal','FontSize',12);
            legend('data', 'fitted','Location','best')
            set(gca,'FontSize',12)

            % Plot linear fit
            subplot(2,1,2)
            if exist('data','var')
                ydata = data.VFAData./sin(flipAngles/180*pi*B1map)';
                xdata = data.VFAData./tan(flipAngles/180*pi*B1map)';
                plot(xdata,ydata,'xb','MarkerSize',16)
                hold on
            end
            slope = exp(-TR/x.T1);
            intercept = x.M0*(1-slope);
            X=Smodel./tan(flipAngles/180*pi*B1map);
            mval = min(X);
            Mval = max(X);
            plot([mval Mval],intercept+slope*[mval Mval],'-r');
            hold off
            title(sprintf('Linear Fit: T1=%0.4f s; M0=%0.0f;',x.T1,x.M0),'FontSize',14);
            xlabel('[au]','FontSize',12);
            ylabel('[au]','FontSize',12);
            legend('linearized data', 'linear fit','Location','best')
            %txt=strcat('T1=',num2str(x.T1),'s M0=',num2str(x.M0));
            set(gca,'FontSize',12)

%             h = plot( fitresult, xData, yData,'+');
%             set(h,'MarkerSize',30)
%             legend( h, 'y vs. x', 'untitled fit 1', 'Location', 'NorthEast' );
%             p11 = predint(fitresult,x,0.95,'observation','off');
%             hold on
%             plot(x,p11,'m--'); drawnow;
%             hold off
%             % Label axes
%             xlabel( 'x' );
%             ylabel( 'y' );
%             grid on
%             saveas(gcf,['temp.jpg']);
       end
       function [FitResults, data] = Sim_Single_Voxel_Curve(obj, x, Opt,display)
           % Simulates Single Voxel
           %
           % :param x: [struct] fit parameters
           % :param Opt.SNR: [struct] signal to noise ratio to use
           % :param display: 1=display, 0=nodisplay
           % :returns: [struct] FitResults, data (noisy dataset)

           if ~exist('display','var'), display = 1; end
           Smodel = equation(obj, x);
           sigma = max(abs(Smodel))/Opt.SNR;
           data.VFAData = ricernd(Smodel,sigma)';
           data.B1map = 1;

           FitResults = fit(obj,data);
           if display
               plotModel(obj, FitResults, data);
           end
       end

       function SimVaryResults = Sim_Sensitivity_Analysis(obj, OptTable, Opt)
           % SimVaryGUI
           SimVaryResults = SimVary(obj, Opt.Nofrun, OptTable, Opt);
       end

       function SimRndResults = Sim_Multi_Voxel_Distribution(obj, RndParam, Opt)
           % SimVaryGUI
           SimRndResults = SimRnd(obj, RndParam, Opt);
       end

    end

    % CLI-only implemented static methods. Can be called directly from
    % class - no object needed.
    methods(Static)
        function Mz = analytical_solution(params)
            %ANALYTICAL_SOLUTION  Analytical equations for the longitudinal magnetization of
            %steady-state gradient echo experiments.
            %
            %   Reference: Stikov, N. , Boudreau, M. , Levesque, I. R.,
            %   Tardif, C. L., Barral, J. K. and Pike, G. B. (2015), On the
            %   accuracy of T1 mapping: Searching for common ground. Magn.
            %   Reson. Med., 73: 514-522. doi:10.1002/mrm.25135
            %
            %   params: Struct.
            %           Properties: T1, TR, EXC_FA, constant (optional)
            %

            Mz = vfa_equation(params);

        end

        function signMaxAngle = ernst_angle(params)
            %ERNST_ANGLE  Analytical equations for the longitudinal magnetization of
            %steady-state gradient echo experiments.
            %
            %   Reference: Ernst, R. R. (1966). "Application of Fourier
            %   transform spectroscopy to magnetic resonance". Review of
            %   Scientific Instruments. 37: 93. doi:10.1063/1.171996
            %
            %   params: Struct.
            %           Properties: T1, TR
            %

            try
                T1 = params.T1;
                TR = params.TR;
            catch
                error('vfa_t1.ernst_equation: Incorrect parameters.  Run `help vfa_t1.ernst_angle` for more info.')
            end

            signMaxAngle = acosd(exp(-TR./T1));

        end

        function [Mz, Msig] = bloch_sim(params)
            %BLOCH_SIM Bloch simulations of the GRE-IR pulse sequence.
            % Simulates 100 spins params. Nex repetitions of the IR pulse
            % sequences.
            %
            % params: Struct with the following fields:
            %   EXC_FA: Excitation pulse flip angle in degrees.
            %   TI: Inversion time (ms).
            %   TR: Repetition time (ms).
            %   TE: Echo time (ms).
            %   T1: Longitudinal relaxation time (ms).
            %   T2: Transverse relaxation time (ms).
            %   Nex: Number of excitations
            %
            %   (optional)
            %       df: Off-resonance frequency of spins relative to excitation pulse (in Hz)
            %       crushFlag: Numeric flag for perfect spoiling (1) or partial spoiling (2).
            %       partialDephasingFlag: do partialDephasing (see below)
            %       partialDephasing: Partial dephasing fraction (between [0, 1]). 1 = no dephasing, 0 = complete dephasing (sele
            %       inc: Phase spoiling increment in degrees.
            %
            % Outputs:
            %   Mz: Longitudinal magnetization at just prior to excitation pulse.
            %   Msig: Complex signal produced by the transverse magnetization at time TE after excitation.
            %

            %% Setup parameters
            %

            alpha = deg2rad(params.EXC_FA);
            TR = params.TR;
            T1 = params.T1;

            TE = params.TE;
            T2 = params.T2;

            Nex = params.Nex;

            %% Optional parameers

            if isfield(params, 'df')
                df = params.df;
            else
                df = 0;
            end

            if isfield(params, 'crushFlag')
                crushFlag = params.crushFlag;
            else
                crushFlag = 1;
            end

            if isfield(params, 'partialDephasingFlag')
                partialDephasingFlag = params.partialDephasingFlag;
            else
                partialDephasingFlag = 0;
            end

            if isfield(params, 'partialDephasing')
                partialDephasing = params.partialDephasing;
            else
                partialDephasing = 1;
            end

            if isfield(params, 'inc')
                inc = deg2rad(params.inc);
            else
                inc = 0;
            end

            %% Simulate for every flip angless
            %

            for ii = 1:length(alpha)

                [Msig(ii), Mz(ii)] = vfa_blochsim(                  ...
                    alpha(ii),            ...
                    T1,                   ...
                    T2,                   ...
                    TE,                   ...
                    TR,                   ...
                    crushFlag,            ...
                    partialDephasingFlag, ...
                    partialDephasing,     ...
                    df,                   ...
                    Nex,                  ...
                    inc                   ...
                    );

            end
        end

        function EXC_FA = find_two_optimal_flip_angles(params, sigFigs)
            %FIND_TWO_OPTIMAL_FLIP_ANGLES Calculate the two optimal flip
            %angles (having signals 71% of the signal at the Ernst angle).
            %
            % Simulates 100 spins params. Nex repetitions of the IR pulse
            % sequences.
            %
            % References:
            %
            % Deoni, S. C., Rutt, B. K. and Peters, T. M. (2003), Rapid
            % combined T1 and T2 mapping using gradient recalled
            % acquisition in the steady state. Magn. Reson. Med., 49:
            % 515-526.
            %
            % Schabel, M.C. & Morrell, G.R., 2009. Uncertainty in T(1)
            % mapping using the variable flip angle method with two flip
            % angles. Physics in medicine and biology, 54(1), pp.N1?8.

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

            % Set up search space of flip angle and signal values
            flipAngleSearchSpace = 0:1*10^(-sigFigs):90;

            params.EXC_FA = flipAngleSearchSpace;
            signals = vfa_t1.analytical_solution(params);

            % Get the Angle, find index in search space
            maxAngle = vfa_t1.ernst_angle(params);
            [~,maxRangeIndex] = min(abs(flipAngleSearchSpace-maxAngle));

            % Calculate signal at optimal flip angle values
            optAngleSignal= 0.71*signals(maxRangeIndex);

            % Calculate indices of optimal flip angles (smaller and larger
            % than the Ernst angle)
            [~,optRangeIndex_small] = min(abs(signals(1:maxRangeIndex)-optAngleSignal));
            [~,optRangeIndex_large_rel] = min(abs(signals(maxRangeIndex+1:end)-optAngleSignal));
            optRangeIndex_large = maxRangeIndex+optRangeIndex_large_rel;

            % Output optimal flip angle values
            EXC_FA = [flipAngleSearchSpace(optRangeIndex_small), flipAngleSearchSpace(optRangeIndex_large)];

        end
    end
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
%% Multi-target planar tracking
%
% Example of using features2d framework for multiple planar targets tracking
% in a video using homography matching. This sample uses the |PlaneTracker|
% class.
%
% Sample video: <http://www.youtube.com/watch?v=pzVbhxx6aog>.
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.3.1/samples/python/plane_tracker.py>
%

% video file, and multiple targets to track [x,y,w,h]
vid = fullfile(mexopencv.root(), 'test', 'blais.mp4');
assert(exist(vid, 'file') == 2, 'Missing video file');
win = [
    136 0 366 433;   % book
    135 165 285 175  % face
];
N = size(win,1);

% open video feed, and get firs frame
cap = cv.VideoCapture(vid);
pause(1);
assert(cap.isOpened(), 'Failed to open video');
frame = cap.read();
assert(~isempty(frame), 'Failed to read frames');

% create and initialize tracker
tracker = PlaneTracker();
for i=1:N
    tracker.addTarget(frame, win(i,:));
end

% prepare plot
hImg = imshow(frame);
clr = lines(N);
for i=1:N
    % keypoints and bounding box for each tracked target
    hPts(i) = line(NaN, NaN, 'Color',clr(i,:), 'LineStyle','none', 'Marker','o');
    hLin(i) = line(NaN, NaN, 'Color',clr(i,:), 'LineWidth',3);
end

% main loop
while ishghandle(hImg)
    % read new frame
    frame = cap.read();
    if isempty(frame), break; end

    % track and update keypoints/boundaries of matched targets
    tracked = tracker.track(frame);
    for i=1:numel(tracked)
        tr = tracked(i);
        set(hPts(tr.index), 'XData',tr.pt1(:,1), 'YData',tr.pt1(:,2));
        set(hLin(tr.index), ...
            'XData',tr.quad([1:end 1],1), 'YData',tr.quad([1:end 1],2));
    end

    % update tracked objects which were not found in current frame
    idx = setdiff(1:N, [tracked.index]);
    if ~isempty(idx)
        set(hPts(idx), 'XData',NaN, 'YData',NaN);
        set(hLin(idx), 'XData',NaN, 'YData',NaN);
    end

    % display result
    set(hImg, 'CData',frame);
    drawnow;
end
cap.release();

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Fig. 6.49   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%

clf
clear all;
%close all;
clf

num=0.01*[20 1];
den=[1 0 0]+[0 num];
w=logspace(-3,1,100);
[m,p]=bode(num,den,w);
loglog(w,m);
axis([.01 1 .1 10])
text(.12,1.2,'T(j\omega)')
text(.031,.15,'S(j\omega)')
bodegrid;
hold on
% add line at mag = 0.707
w2=[.01 1];
mcl=[.707 .707];
loglog(w2,mcl,'r')

% now error FR
numE=[1 0 0];
denE=den;
[me,pe]=bode(numE,denE,w);
loglog(w,me)

xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Fig. 6.49 Closed-loop frequency response.');
bodegrid;
hold off

author: Sable

repo: mcbench-benchmarks


text:
function [ic,ik] = bspdegelev(d,c,k,t)
% BSPDEGELEV  Degree elevate a univariate B-Spline.
% -------------------------------------------------------------------------
% ADAPTATION of BSPDEGELEV from C Routine
% -------------------------------------------------------------------------
%
% Calling Sequence:
%
%   [ic,ik] = bspdegelev(d,c,k,t)
%
% Parameters:
%
%   d - Degree of the B-Spline.
%   c - Control points, matrix of size (dim,nc).
%   k - Knot sequence, row vector of size nk.
%   t - Raise the B-Spline degree t times.
%
%   ic - Control points of the new B-Spline.
%   ik - Knot vector of the new B-Spline.
%

[mc,nc] = size(c);
                                                          %
                                                          % int bspdegelev(int d, double *c, int mc, int nc, double *k, int nk,
                                                          %                int t, int *nh, double *ic, double *ik)
                                                          % {
                                                          %   int row,col;
                                                          %
                                                          %   int ierr = 0;
                                                          %   int i, j, q, s, m, ph, ph2, mpi, mh, r, a, b, cind, oldr, mul;
                                                          %   int n, lbz, rbz, save, tr, kj, first, kind, last, bet, ii;
                                                          %   double inv, ua, ub, numer, den, alf, gam;
                                                          %   double **bezalfs, **bpts, **ebpts, **Nextbpts, *alfs;
                                                          %
                                                          %   double **ctrl  = vec2mat(c, mc, nc);
% ic = zeros(mc,nc*(t));                                  %   double **ictrl = vec2mat(ic, mc, nc*(t+1));
                                                          %
n = nc - 1;                                               %   n = nc - 1;
                                                          %
bezalfs =  zeros(d+1,d+t+1);                              %   bezalfs = matrix(d+1,d+t+1);
bpts = zeros(mc,d+1);                                     %   bpts = matrix(mc,d+1);
ebpts = zeros(mc,d+t+1);                                  %   ebpts = matrix(mc,d+t+1);
Nextbpts = zeros(mc,d+1);                                 %   Nextbpts = matrix(mc,d+1);
alfs = zeros(d,1);                                        %   alfs = (double *) mxMalloc(d*sizeof(double));
                                                          %
m = n + d + 1;                                            %   m = n + d + 1;
ph = d + t;                                               %   ph = d + t;
ph2 = floor(ph / 2);                                      %   ph2 = ph / 2;
                                                          %
                                                          %   // compute bezier degree elevation coefficeients
bezalfs(1,1) = 1;                                         %   bezalfs[0][0] = bezalfs[ph][d] = 1.0;
bezalfs(d+1,ph+1) = 1;                                    %

for i=1:ph2                                               %   for (i = 1; i <= ph2; i++) {
   inv = 1/bincoeff(ph,i);                                %     inv = 1.0 / bincoeff(ph,i);
   mpi = min(d,i);                                        %     mpi = min(d,i);
                                                          %
   for j=max(0,i-t):mpi                                   %     for (j = max(0,i-t); j <= mpi; j++)
       bezalfs(j+1,i+1) = inv*bincoeff(d,j)*bincoeff(t,i-j);  %       bezalfs[i][j] = inv * bincoeff(d,j) * bincoeff(t,i-j);
   end
end                                                       %   }
                                                          %
for i=ph2+1:ph-1                                          %   for (i = ph2+1; i <= ph-1; i++) {
   mpi = min(d,i);                                        %     mpi = min(d, i);
   for j=max(0,i-t):mpi                                   %     for (j = max(0,i-t); j <= mpi; j++)
       bezalfs(j+1,i+1) = bezalfs(d-j+1,ph-i+1);          %       bezalfs[i][j] = bezalfs[ph-i][d-j];
   end
end                                                       %   }
                                                          %
mh = ph;                                                  %   mh = ph;
kind = ph+1;                                              %   kind = ph+1;
r = -1;                                                   %   r = -1;
a = d;                                                    %   a = d;
b = d+1;                                                  %   b = d+1;
cind = 1;                                                 %   cind = 1;
ua = k(1);                                                %   ua = k[0];
                                                          %
for ii=0:mc-1                                             %   for (ii = 0; ii < mc; ii++)
   ic(ii+1,1) = c(ii+1,1);                                %     ictrl[0][ii] = ctrl[0][ii];
end                                                       %
for i=0:ph                                                %   for (i = 0; i <= ph; i++)
   ik(i+1) = ua;                                          %     ik[i] = ua;
end                                                       %
                                                          %   // initialise first bezier seg
for i=0:d                                                 %   for (i = 0; i <= d; i++)
   for ii=0:mc-1                                          %     for (ii = 0; ii < mc; ii++)
      bpts(ii+1,i+1) = c(ii+1,i+1);                       %       bpts[i][ii] = ctrl[i][ii];
   end
end                                                       %
                                                          %   // big loop thru knot vector
while b < m                                               %   while (b < m)  {
   i = b;                                                 %     i = b;
   while b < m && k(b+1) == k(b+2)                        %     while (b < m && k[b] == k[b+1])
      b = b + 1;                                          %       b++;
   end                                                    %
   mul = b - i + 1;                                       %     mul = b - i + 1;
   mh = mh + mul + t;                                     %     mh += mul + t;
   ub = k(b+1);                                           %     ub = k[b];
   oldr = r;                                              %     oldr = r;
   r = d - mul;                                           %     r = d - mul;
                                                          %
                                                          %     // insert knot u(b) r times
   if oldr > 0                                            %     if (oldr > 0)
      lbz = floor((oldr+2)/2);                            %       lbz = (oldr+2) / 2;
   else                                                   %     else
      lbz = 1;                                            %       lbz = 1;
   end                                                    %

   if r > 0                                               %     if (r > 0)
      rbz = ph - floor((r+1)/2);                          %       rbz = ph - (r+1)/2;
   else                                                   %     else
      rbz = ph;                                           %       rbz = ph;
   end                                                    %

   if r > 0                                               %     if (r > 0) {
                                                          %       // insert knot to get bezier segment
      numer = ub - ua;                                    %       numer = ub - ua;
      for q=d:-1:mul+1                                    %       for (q = d; q > mul; q--)
         alfs(q-mul) = numer / (k(a+q+1)-ua);             %         alfs[q-mul-1] = numer / (k[a+q]-ua);
      end

      for j=1:r                                           %       for (j = 1; j <= r; j++)  {
         save = r - j;                                    %         save = r - j;
         s = mul + j;                                     %         s = mul + j;
                                                          %
         for q=d:-1:s                                     %         for (q = d; q >= s; q--)
            for ii=0:mc-1                                 %           for (ii = 0; ii < mc; ii++)
               tmp1 = alfs(q-s+1)*bpts(ii+1,q+1);
               tmp2 = (1-alfs(q-s+1))*bpts(ii+1,q);
               bpts(ii+1,q+1) = tmp1 + tmp2;              %             bpts[q][ii] = alfs[q-s]*bpts[q][ii]+(1.0-alfs[q-s])*bpts[q-1][ii];
            end
         end                                              %

         for ii=0:mc-1                                    %         for (ii = 0; ii < mc; ii++)
            Nextbpts(ii+1,save+1) = bpts(ii+1,d+1);       %           Nextbpts[save][ii] = bpts[d][ii];
         end
      end                                                 %       }
   end                                                    %     }
                                                          %     // end of insert knot
                                                          %
                                                          %     // degree elevate bezier
   for i=lbz:ph                                           %     for (i = lbz; i <= ph; i++)  {
      for ii=0:mc-1                                       %       for (ii = 0; ii < mc; ii++)
         ebpts(ii+1,i+1) = 0;                             %         ebpts[i][ii] = 0.0;
      end
      mpi = min(d, i);                                    %       mpi = min(d, i);
      for j=max(0,i-t):mpi                                %       for (j = max(0,i-t); j <= mpi; j++)
         for ii=0:mc-1                                    %         for (ii = 0; ii < mc; ii++)
            tmp1 = ebpts(ii+1,i+1);
            tmp2 = bezalfs(j+1,i+1)*bpts(ii+1,j+1);
            ebpts(ii+1,i+1) = tmp1 + tmp2;                %           ebpts[i][ii] = ebpts[i][ii] + bezalfs[i][j]*bpts[j][ii];
         end
      end
   end                                                    %     }
                                                          %     // end of degree elevating bezier
                                                          %
   if oldr > 1                                            %     if (oldr > 1)  {
                                                          %       // must remove knot u=k[a] oldr times
      first = kind - 2;                                                    %       first = kind - 2;
      last = kind;                                        %       last = kind;
      den = ub - ua;                                      %       den = ub - ua;
      bet = floor((ub-ik(kind)) / den);                   %       bet = (ub-ik[kind-1]) / den;
                                                          %
                                                          %       // knot removal loop
      for tr=1:oldr-1                                     %       for (tr = 1; tr < oldr; tr++)  {
         i = first;                                       %         i = first;
         j = last;                                        %         j = last;
         kj = j - kind + 1;                               %         kj = j - kind + 1;
         while j-i > tr                                   %         while (j - i > tr)  {
                                                          %           // loop and compute the new control points
                                                          %           // for one removal step
            if i < cind                                   %           if (i < cind)  {
               alf = (ub-ik(i+1))/(ua-ik(i+1));           %             alf = (ub-ik[i])/(ua-ik[i]);
               for ii=0:mc-1                              %             for (ii = 0; ii < mc; ii++)
                  tmp1 = alf*ic(ii+1,i+1);
                  tmp2 = (1-alf)*ic(ii+1,i);
                  ic(ii+1,i+1) = tmp1 + tmp2;             %               ictrl[i][ii] = alf * ictrl[i][ii] + (1.0-alf) * ictrl[i-1][ii];
               end
            end                                           %           }
            if j >= lbz                                   %           if (j >= lbz)  {
               if j-tr <= kind-ph+oldr                    %             if (j-tr <= kind-ph+oldr) {
                  gam = (ub-ik(j-tr+1)) / den;            %               gam = (ub-ik[j-tr]) / den;
                  for ii=0:mc-1                           %               for (ii = 0; ii < mc; ii++)
                     tmp1 = gam*ebpts(ii+1,kj+1);
                     tmp2 = (1-gam)*ebpts(ii+1,kj+2);
                     ebpts(ii+1,kj+1) = tmp1 + tmp2;      %                 ebpts[kj][ii] = gam*ebpts[kj][ii] + (1.0-gam)*ebpts[kj+1][ii];
                  end                                     %             }
               else                                       %             else  {
                  for ii=0:mc-1                           %               for (ii = 0; ii < mc; ii++)
                     tmp1 = bet*ebpts(ii+1,kj+1);
                     tmp2 = (1-bet)*ebpts(ii+1,kj+2);
                     ebpts(ii+1,kj+1) = tmp1 + tmp2;      %                 ebpts[kj][ii] = bet*ebpts[kj][ii] + (1.0-bet)*ebpts[kj+1][ii];
                  end
               end                                        %             }
            end                                           %           }
            i = i + 1;                                    %           i++;
            j = j - 1;                                    %           j--;
            kj = kj - 1;                                  %           kj--;
         end                                              %         }
                                                          %
         first = first - 1;                               %         first--;
         last = last + 1;                                 %         last++;
      end                                                 %       }
   end                                                    %     }
                                                          %     // end of removing knot n=k[a]
                                                          %
                                                          %     // load the knot ua
   if a ~= d                                              %     if (a != d)
      for i=0:ph-oldr-1                                   %       for (i = 0; i < ph-oldr; i++)  {
         ik(kind+1) = ua;                                 %         ik[kind] = ua;
         kind = kind + 1;                                 %         kind++;
      end
   end                                                    %       }
                                                          %
                                                          %     // load ctrl pts into ic
      for j=lbz:rbz                                       %     for (j = lbz; j <= rbz; j++)  {
         for ii=0:mc-1                                    %       for (ii = 0; ii < mc; ii++)
            ic(ii+1,cind+1) = ebpts(ii+1,j+1);            %         ictrl[cind][ii] = ebpts[j][ii];
         end
         cind = cind + 1;                                 %       cind++;
      end                                                 %     }
                                                          %
      if b < m                                            %     if (b < m)  {
                                                          %       // setup for next pass thru loop
         for j=0:r-1                                      %       for (j = 0; j < r; j++)
            for ii=0:mc-1                                 %         for (ii = 0; ii < mc; ii++)
               bpts(ii+1,j+1) = Nextbpts(ii+1,j+1);       %           bpts[j][ii] = Nextbpts[j][ii];
            end
         end
         for j=r:d                                        %       for (j = r; j <= d; j++)
            for ii=0:mc-1                                 %         for (ii = 0; ii < mc; ii++)
               bpts(ii+1,j+1) = c(ii+1,b-d+j+1);          %           bpts[j][ii] = ctrl[b-d+j][ii];
            end
         end
         a = b;                                           %       a = b;
         b = b+1;                                         %       b++;
         ua = ub;                                         %       ua = ub;
                                                          %     }
      else                                                %     else
                                                          %       // end knot
         for i=0:ph                                       %       for (i = 0; i <= ph; i++)
            ik(kind+i+1) = ub;                            %         ik[kind+i] = ub;
         end
      end
end                                                       %   }
% End big while loop                                      %   // end while loop
                                                          %
                                                          %   *nh = mh - ph - 1;
                                                          %
                                                          %   freevec2mat(ctrl);
                                                          %   freevec2mat(ictrl);
                                                          %   freematrix(bezalfs);
                                                          %   freematrix(bpts);
                                                          %   freematrix(ebpts);
                                                          %   freematrix(Nextbpts);
                                                          %   mxFree(alfs);
                                                          %
                                                          %   return(ierr);
                                                          % }

function b = bincoeff(n,k)
%  Computes the binomial coefficient.
%
%      ( n )      n!
%      (   ) = --------
%      ( k )   k!(n-k)!
%
%  b = bincoeff(n,k)
%
%  Algorithm from 'Numerical Recipes in C, 2nd Edition' pg215.

                                                          % double bincoeff(int n, int k)
                                                          % {
b = floor(0.5+exp(factln(n)-factln(k)-factln(n-k)));      %   return floor(0.5+exp(factln(n)-factln(k)-factln(n-k)));
                                                          % }

function f = factln(n)
% computes ln(n!)
if n <= 1, f = 0; return, end
f = gammaln(n+1); %log(factorial(n));


text:
function varargout = airy(varargin)

switch class(varargin{1})

    case 'sdpvar'
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'

        operator = CreateBasicOperator('callback');
        operator.derivative = @(x)(airy(1,x));
        operator.range = [-4.1901547803e-01 5.35656656014e-01];

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

    otherwise
        error(['SDPVAR/' upper(mfilename) ' called with weird argument']);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [p,q] = detectAndMatchFeatures(image1,image1ROI,image2,image2ROI,featureOpts)
%detectAndMatchFeatures: A wrapper for detecting and matching features
%between two images using various feature detection methods.

%Inputs
%
%image1: The first image (2D unit8, grayscale or rgb).
%
%image1ROI: A region-of-interest specifying which region(s) in image1 we
%should detect features in. If image1ROI=[] then the whole image is used.
%
%image2: The second image (2D unit8, grayscale or rgb).
%
%image2ROI: A region-of-interest specifying which region(s) in image2 we
%should detect features in. If image2ROI=[] then the whole image is used.
%
%featureOpts: Options field specifying the feature detection method.
%Currently only SURF (matlab's built-in) and ASIFT
%(http://www.ipol.im/pub/art/2011/my-asift/) are supported, but it is very
%easy to introduce others.
%for ASIFT, featureOpts should have the structures:
%
%   featureOpts.featureMethod = 'ASIFT';
%   featureOpts.asiftPath (the path to the asift code & compiled executable).
%
%for SURF, featureOpts should have the structures:
%
%   featureOpts.featureMethod = 'SURF'; %careful, unlike ASIFT this only works well
%when the plane's viewpoint is not too tilted!!
%
%   featureOpts.loweRatioThreshold = 1.2 (default) You need to
%set a confidence ratio (see Lowe's SIFT paper for the explanation of
%this). Basically, a high value means only using feature matches that are
%likely to be correct (but at the cost of fewer feature matches). A default
%of 1.2 is usually fine.
%
% outputs:
%p : 2XN matrix holding the matched points in image1
%q : 2XN matrix holding the matched points in image2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of the IPPE package for fast plane-based pose
% estimation from the paper "Infinitesimal Plane-based Pose Estimation" by Toby Collins and Adrien Bartoli,
% published in the International Journal of Computer Vision, September
% 2014. A copy of the author's pre-print version can be found here:
%
% http://isit.u-clermont1.fr/~ab/Publications/Collins_Bartoli_IJCV14.pdf
%
% This package is free and covered by the BSD licence without any warranty. We hope you find this code useful and please cite our paper in your work:
% (c) Toby Collins 2015
%
%
%@article{
%year={2014},
%issn={0920-5691},
%journal={International Journal of Computer Vision},
%volume={109},
%number={3},
%doi={10.1007/s11263-014-0725-5},
%title={Infinitesimal Plane-Based Pose Estimation},
%url={http://dx.doi.org/10.1007/s11263-014-0725-5},
%publisher={Springer US},
%keywords={Plane; Pose; SfM; PnP; Homography},
%author={Collins, Toby and Bartoli, Adrien},
%pages={252-286},
%language={English}
%}
%
%
% This is free software covered by the FreeBSD License (see IPPE_license.txt) with Copyright (c) 2015 Toby Collins

%basic argument checking
if nargin~=5
   error('detectAndMatchFeatures has 5 input arguments.');
end
if nargout>2
   error('detectAndMatchFeatures has 2 output arguments.');
end

assert(nargin ==5);
assert(size(image1,3)==1|size(image1,3)==3);
assert(size(image2,3)==1|size(image2,3)==3);
assert(isa(image1,'uint8'));
assert(isa(image2,'uint8'));

%check whether regions-of-interest are used:
if isempty(image1ROI)
    image1ROI = true(size(image1,1),size(image1,2));
end
if isempty(image2ROI)
    image2ROI = true(size(image2,1),size(image2,2));
end

assert(isa(image1ROI,'logical'));
assert(isa(image1ROI,'logical'));

%convert image to grayscale:
if size(image1,3)==3
    image1 = rgb2gray(image1);
end

if size(image2,3)==3
    image2 = rgb2gray(image2);
end

%perform detection and matching:
switch featureOpts.featureMethod
    case 'SURF'
        %detection:
        pointsTemplate = detectSURFFeatures(image1,'MetricThreshold',200);
        featuresTemplate = extractFeatures(image1,pointsTemplate);
        p = pointsTemplate.Location';

        pointsInput = detectSURFFeatures(image2,'MetricThreshold',200);
        featuresInput = extractFeatures(image2,pointsInput);

        %matching using Lowe's ratio test for rejecting bad matches (see his SIFT paper for details)
        [IDX, D]= knnsearch(featuresInput,featuresTemplate,'K',2);
        vlds = D(:,2)./D(:,1)> featureOpts.loweRatioThreshold;
        q = pointsInput.Location(IDX(:,1),:)';
        q = q(:,vlds);
        p = p(:,vlds);
    case 'ASIFT'
        [p,q] = asiftWrapper(featureOpts.asiftPath,image1,image2);
    otherwise
        error('unknown feature detection method is specified');
end

%keep only features located in the rois:
vlds = interp2(double(image1ROI),p(1,:),p(2,:))==1;
vlds = vlds & interp2(double(image2ROI),q(1,:),q(2,:))==1;
p = p(:,vlds);
q = q(:,vlds);

author: tobycollins

repo: IPPE

save_path: github-repos/MATLAB/tobycollins-IPPE


text:
function tet_mesh_test005 ( )

%*****************************************************************************80
%
%% TET_MESH_TEST005 tests TETRAHEDRON_BARYCENTRIC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2009
%
%  Author:
%
%    John Burkardt
%
  test1_num = 3;
  test2_num = 5;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST005\n' );
  fprintf ( 1, '  TETRAHEDRON_BARYCENTRIC computes the barycentric\n' );
  fprintf ( 1, '  coordinates of a point.\n' );
%
%  Choose a random tetrahedron.
%
  for test1 = 1 : test1_num

    [ tet_xyz, seed ] = r8mat_uniform_01 ( 3, 4, seed );

    r8mat_transpose_print ( 3, 4, tet_xyz, '  Random tetrahedron:' );
%
%  Choose barycentric coordinates C1 at random.
%
%  Define a point P.
%
%  Have TETRAHEDRON_BARYCENTRIC compute C2, the barycentric coordinates of P.
%
    for test2 = 1 : test2_num

      [ c1, seed ] = r8vec_uniform_01 ( 4, seed );
      c1_sum = sum ( c1(1:4) );
      c1(1:4) = c1(1:4) / c1_sum;

      p(1:3) = tet_xyz(1:3,1:4) * c1(1:4)';

      c2 = tetrahedron_barycentric ( tet_xyz, p );

      fprintf ( 1, '\n' );
      fprintf ( 1, '  C1 =  %14.6f  %14.6f  %14.6f  %14.6f\n', c1 );
      fprintf ( 1, '  C2 =  %14.6f  %14.6f  %14.6f  %14.6f\n', c2 );

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function modelRev = convertToReversible(model)
% Converts a model structure from irreversible format to
% reversible format
%
% USAGE:
%
%    modelRev = convertToReversible(model)
%
% INPUT:
%    model:       COBRA model in irreversible format (forward/backward reactions separated)
%                 forward reactions have to be marked with an _f in the id,
%                 backward with an _b. reverse reactions (i.e. reactions
%                 which were originally only allowed to carry negative flux
%                 must have a _r tag at the end of the reaction id.
%                 It is further assumed, that there is an exact 1:1 match
%                 between _f and _b reactions.
%
% OUTPUT:
%    modelRev:    Model in reversible format
%
% .. Author: - Greg Hannum 7/22/05
%
%
%

nRxns = length(model.rxns);
%indicator for reactions already handled
rxnProcessed = false(length(model.rxns),1);
%get the basic names of the reactions.
reactionbasenames = regexprep(model.rxns,'(_f|_r|_b)$','');
%Get the reaction positions which need to be reversed
reversePos = ~cellfun(@isempty ,regexp(model.rxns,'(_r|_b)$'));
%and those, which will be deleted (all backwards. Since backward only
%reactions are marked with an _r, we should always have a _f which is the
%other matching reaction.
reactionsToRemove = ~cellfun(@isempty ,regexp(model.rxns,'(_b)$'));

%and replace the names
model.rxns = reactionbasenames;

%We can directly update the lower bounds and stoichiometries of reverse reactions
%The assumption is, that all reverse reaction bounds and forward reaction
%bounds are consistent.
templbs = model.lb(reversePos);
model.lb(reversePos) = -model.ub(reversePos);
model.ub(reversePos) = -templbs;
model.S(:,reversePos) = -model.S(:,reversePos);

%update all reactions

for i = 1:numel(reactionbasenames)
    if(rxnProcessed(i))
        continue;
    end
    matches = ismember(reactionbasenames,reactionbasenames{i});
    rxnProcessed = rxnProcessed | matches;
    if(sum(matches) > 1)
        %We only need to do something if we have a forward/backward pair.
        %otherwise we don't do anything (since the input format is not
        %matched.
        revpos = reversePos & matches;
        fwpos = ~reversePos & matches;
        model.lb(fwpos) = min(model.lb(revpos), model.lb(fwpos));
        model.ub(fwpos) = max(model.ub(revpos), model.ub(fwpos));
    end
end

%Remove the surplus model fields (they should be the same anyways)
modelRev = removeFieldEntriesForType(model,reactionsToRemove,'rxns',nRxns);

%and remove a potential "match field
if isfield(modelRev, 'match')
    modelRev = rmfield(modelRev,'match');
end
modelRev.reversibleModel = true;


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: derivative check (use data from E3_splineInterpolation2D)
%
%==============================================================================

% setup test data
dataT = flipud([1,2,3,4;1,2,3,4;4,4,4,4])';
m     = size(dataT);
omega = [0,m(1),0,m(2)];
B     = @(i) spdiags(ones(m(i),1)*[1,4,1],[-1:1],m(i),m(i));
T     = B(1)\dataT/B(2);
xf    = reshape(getCellCenteredGrid(omega,10*m),[],2);

fctn = @(x) splineInter(T,omega,x);
figure(1); clf;
[fig,ph,th] = checkDerivative(fctn,xf(:),'fig',1);
%==============================================================================


text:
%% exportGif
% Below is a demonstration of the features of the |exportGif| function

%%
clear; close all; clc;

%% Syntax
% |exportGif(optStruct);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*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 [ n_data, a, b, lambda, x, fx ] = beta_noncentral_cdf_values ( n_data )

%*****************************************************************************80
%
%% BETA_NONCENTRAL_CDF_VALUES returns some values of the noncentral Beta CDF.
%
%  Discussion:
%
%    The values presented here are taken from the reference, where they
%    were given to a limited number of decimal places.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    R Chattamvelli, R Shanmugam,
%    Algorithm AS 310:
%    Computing the Non-central Beta Distribution Function,
%    Applied Statistics,
%    Volume 46, Number 1, 1997, pages 146-156.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0
%    before the first call.  On each call, the routine increments N_DATA by 1,
%    and returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, rea A, B, the shape parameters.
%
%    Output, real LAMBDA, the noncentrality parameter.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 25;

  a_vec = [ ...
        5.0, ...
        5.0, ...
        5.0, ...
       10.0, ...
       10.0, ...
       10.0, ...
       20.0, ...
       20.0, ...
       20.0, ...
       10.0, ...
       10.0, ...
       15.0, ...
       20.0, ...
       20.0, ...
       20.0, ...
       30.0, ...
       30.0, ...
       10.0, ...
       10.0, ...
       10.0, ...
       15.0, ...
       10.0, ...
       12.0, ...
       30.0, ...
       35.0 ];
  b_vec = [ ...
        5.0, ...
        5.0, ...
        5.0, ...
       10.0, ...
       10.0, ...
       10.0, ...
       20.0, ...
       20.0, ...
       20.0, ...
       20.0, ...
       10.0, ...
        5.0, ...
       10.0, ...
       30.0, ...
       50.0, ...
       20.0, ...
       40.0, ...
        5.0, ...
       10.0, ...
       30.0, ...
       20.0, ...
        5.0, ...
       17.0, ...
       30.0, ...
       30.0 ];
  fx_vec = [ ...
       0.4563021, ...
       0.1041337, ...
       0.6022353, ...
       0.9187770, ...
       0.6008106, ...
       0.0902850, ...
       0.9998655, ...
       0.9925997, ...
       0.9641112, ...
       0.9376626573, ...
       0.7306817858, ...
       0.1604256918, ...
       0.1867485313, ...
       0.6559386874, ...
       0.9796881486, ...
       0.1162386423, ...
       0.9930430054, ...
       0.0506899273, ...
       0.1030959706, ...
       0.9978417832, ...
       0.2555552369, ...
       0.0668307064, ...
       0.0113601067, ...
       0.7813366615, ...
       0.8867126477 ];
  lambda_vec = [ ...
        54.0, ...
       140.0, ...
       170.0, ...
        54.0, ...
       140.0, ...
       250.0, ...
        54.0, ...
       140.0, ...
       250.0, ...
       150.0, ...
       120.0, ...
        80.0, ...
       110.0, ...
        65.0, ...
       130.0, ...
        80.0, ...
       130.0, ...
        20.0, ...
        54.0, ...
        80.0, ...
       120.0, ...
        55.0, ...
        64.0, ...
       140.0, ...
        20.0 ];
  x_vec = [ ...
       0.8640, ...
       0.9000, ...
       0.9560, ...
       0.8686, ...
       0.9000, ...
       0.9000, ...
       0.8787, ...
       0.9000, ...
       0.9220, ...
       0.868, ...
       0.900, ...
       0.880, ...
       0.850, ...
       0.660, ...
       0.720, ...
       0.720, ...
       0.800, ...
       0.644, ...
       0.700, ...
       0.780, ...
       0.760, ...
       0.795, ...
       0.560, ...
       0.800, ...
       0.670 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    a = 0.0;
    b = 0.0;
    lambda = 0.0;
    x = 0.0;
    fx = 0.0;
  else
    a = a_vec(n_data);
    b = b_vec(n_data);
    lambda = lambda_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function d = testing(a,d)

  d=set_x(d,feval('get_kernel',a,d,a.dat));

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [HDR] = save2mm(fn,MM,montage);
% SAVE2MM  saves Matrix into MatrixMarket format
%
%       HDR = save2mm(filename,M,comment);
%
% filename	destination file
% M		Matrix
%
% see also: SLOAD, getMontage, regress_eog, save2gdf
%

%
%    BioSig 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.
%
%    BioSig 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 BioSig.  If not, see <http://www.gnu.org/licenses/>.

% 	$Id$
%	Copyright (C) 2010 by Alois Schloegl <a.schloegl@ieee.org>
%       This file is part of the biosig project http://biosig.sf.net/

        [I,J,V] = find(MM);
        fid = fopen(fn,'w+');
        fprintf(fid,'%%%%MatrixMarket matrix coordinate real general\n');
        fprintf(fid,'%% generated on %04i-%02i-%02i %02i:%02i:%02.0f\n',clock);

        if ischar(montage) m = montage; else m = '? (user specified)'; end;
        fprintf(fid,'%% Spatial Filter for %s \n',m);
        fprintf(fid,'%i %i %i\n',size(MM),length(V));

        for k = 1:length(V),
                fprintf(fid,'%2i %2i %f\n',I(k),J(k),V(k));
        end;
        fclose(fid);

        HDR.Calib = MM;
        HDR.FileName = fn;
        HDR.TYPE = 'MatrixMarket';

author: PatternRecognition

repo: OpenBMI


text:
function output =EH_AHE(input,conf)

output=input; %output is the input after AHE
[m,n,~]=size(input);

%ECR
grid_cols=conf.grid;
grid_rows=conf.grid;
limit = conf.limit;

grid_width=int32(fix(m/grid_cols));
grid_height=int32(fix(n/grid_rows));

map=zeros(grid_cols,grid_rows,256);

%for each grid,we create their mapping function
for i=1:grid_cols
    for j=1:grid_rows
        map(i,j,:)=MakeHistogram(input,1+(i-1)*grid_width,1+(j-1)*grid_height,grid_width,grid_height,limit);
    end
end

%interpolate
%boundary cases I followed the Karel Zuiderveld's implement(C version)
xi = 1;
for i = 1:grid_cols+1
    if i == 1
        subx = grid_width/2;
        xu = 1;
        xd = 1;
    elseif i == grid_cols+1
        subx = grid_width/2;
        xu =  grid_cols;
        xd =  grid_cols;
    else
        subx = grid_width;
        xu = i - 1;
        xd = i;
    end
    yi = 1;
    for j = 1:grid_rows+1
        if j == 1
            suby = grid_height/2;
            yl = 1;
            yr = 1;
        elseif j == grid_rows+1
            suby = grid_height/2;
            yl = grid_rows;
            yr = grid_rows;
        else
            suby = grid_height;
            yl = j - 1;
            yr = j;
        end
        UL = map(xu,yl,:);
        UR = map(xu,yr,:);
        DL = map(xd,yl,:);
        DR = map(xd,yr,:);

        subinput = input(xi:xi+subx-1,yi:yi+suby-1);
        subinput = Interpolate(subinput,UL,UR,DL,DR,subx,suby);
        output(xi:xi+subx-1,yi:yi+suby-1) = subinput;
        yi = yi + suby;
    end
    xi = xi + subx;
end

end

author: AomanHao

repo: Matlab-Image-Dehaze-Enhance


text:
function varargout= proc_laplacian(dat, varargin)
%PROC_LAPLACIAN - Apply spatial Laplacian filter to signals
%
%Synopsis:
% [DAT_LAP, LAP_W]= proc_laplacian(DAT, <OPT>)
%
%Arguments:
% DAT: data structure of continuous or epoched data
% OPT: struct or proerty/value list of optional properties:
%  .filterType    - {small, large, horizontal, vertical, diagonal, eight},
%                     default 'small'
%  .clab - channels that are to be obtained, default '*'.
%  .ignoreClab: labels of channels to be ignored, default {'E*'}.
%  .requireCompleteNeighborhood
%  .requireCompleteOutput
%  .verbose
%
%Returns:
%  DAT_LAP: updated data structure
%  LAP_W:   filter matrix that can be used, e.g. in proc_linearDerivation
%
%See also:
%  procutil_getClabForLaplacian, proc_linearDerivation

%        Benjamin Blankertz
% 07-2012 Johannes Hoehne - Updated documentation and parameter naming

props= {'clab'                         '*'          '!CHAR';
        'IgnoreCLab'                   {'E*'}       'CHAR|CELL{CHAR}'
        'CopyCLab'                     {'E*'}       'CHAR|CELL{CHAR}'
        'requireCompleteNeighborhood'  1            '!BOOL'
        'requireCompleteOutput'        0            '!BOOL'
        'appendix'                     ' lap'       'CHAR'
        'verbose'                      0            'BOOL'
        'GridFcn'                      @util_gridForLaplacian 'FUNC'
        'FilterType'                   'small'      'CHAR(small large horizontal vertical bip_to_anterior bip_to_posterior bip_to_left bip_to_right diagonal diagonal_small)'
        };

if nargin==0,
  varargout{1} = props;
  return
end

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

misc_checkType(dat, 'STRUCT(clab)');

if ~iscell(opt.IgnoreCLab),
  opt.IgnoreCLab= {opt.IgnoreCLab};
end

laplace= [];
laplace.grid= opt.GridFcn();
if isnumeric(opt.FilterType),
  laplace.filter= opt.FilterType;
else
  laplace.filter= procutil_lapGetLaplacianFilter(opt.FilterType);
end

if ~isstruct(dat),
  dat= struct('clab',{dat});
end

rc= util_chanind(dat, {'not', opt.IgnoreCLab{:}});
nOrigChans= length(dat.clab);
pos= zeros(2, nOrigChans);
for ic= 1:nOrigChans,
  pos(:,ic)= procutil_lapGetCoordinates(dat.clab{ic}, laplace.grid);
end
pos(:,setdiff(1:nOrigChans,rc,'legacy'))= inf;

idx_tbf= util_chanind(dat, opt.clab);
W= zeros(length(dat.clab), length(idx_tbf));
clab = [];
lc= 0;
for ci= 1:length(idx_tbf),
  cc= idx_tbf(ci);
  refChans= [];
  nRefs= size(laplace.filter,2);
  for ir= 1:nRefs,
    ri= find( pos(1,:)==pos(1,cc)+laplace.filter(1,ir) & ...
      pos(2,:)==pos(2,cc)+laplace.filter(2,ir) );
    refChans= [refChans ri];
  end
  if length(refChans)==nRefs | ~opt.requireCompleteNeighborhood,
    lc= lc+1;
    W(cc,lc)= 1;
    if ~isempty(refChans),
      W(refChans,lc)= -1/length(refChans);
    end
    clab= [clab, dat.clab(cc)];
    if opt.verbose,
      fprintf('%s: ref''ed to: %s\n', ...
        dat.clab{cc}, str_vec2str(dat.clab(refChans)));
    end
  elseif opt.requireCompleteOutput,
    error('channel %s has incomplete neighborhood', dat.clab{cc});
  end
end
W= W(:,1:lc);

if isfield(dat, 'x'),
  out= proc_linearDerivation(dat, W, 'clab', strcat(clab, opt.appendix));
  if ~isempty(opt.CopyCLab),
    idx= util_chanind(dat, opt.CopyCLab);
    if ~isempty(idx),
      out.x= cat(2, out.x, dat.x(:,idx,:));
      out.clab= cat(2, out.clab, dat.clab(idx));
    end
  end
  varargout= {out, W};
else
  varargout= {W};
end

author: bbci


text:
function outState = FDSUpdateMap(state, newMap)
%         enums
        NEW = -1;
        OPEN = 1;
        CLOSED = 0;

        startPos = state.startPos;
        endPos =  state.endPos;
        pattern = state.pattern;
        ucc = state.ucc;
        graph = state.graph;
        kM = state.kM;
        SQRT2 = sqrt(2)-1;
        stack = state.stack;

        comparator = pcmp;
        stack2= java.util.PriorityQueue(180247, comparator);
        comparator = pcmp;
        stack3= java.util.PriorityQueue(180247, comparator);

        map = newMap;
        difference = map - state.map;
        state.map = map;
        %% ---------------addObstacle
        [x, y] = find(difference == -1);
        indices = [x,y];
        indices(:,3:4) = 0;

        Ind = [indices'];

        for n=indices'
            for s=state.pattern
                Ind = [Ind, n+s];
            end
        end

        added = unique(Ind', 'rows')';
        %% ---------------rmvObstacle
        [x, y] = find(difference == 1);
        indices = [x,y];
        indices(:,3:4) = 0;

        Ind = [indices'];

        for n=indices'
            for s=state.pattern
                Ind = [Ind, n+s];
            end
        end

        removed = unique(Ind', 'rows')';
        %% ----------------- process
        for u = removed
             tmp = inf;
                for n = ucc
                       s = u+n;
                            insert(s, h(s))
                end
        end
        for u = added
             tmp = inf;
                for n = ucc
                       s = u+n;
                        if t(s) == CLOSED && h(s) ~= inf
                            seth(s, inf);
                            s(3:4)  = [-k(s)+kM; k(s)+g(s)];
                            add(stack, s);
                            sett(s, OPEN);
                        end
                end
        end
        outState = state;
        outState.graph = graph;
        outState.map = map;
        outState.stack = stack;
%%
function out = inQ(s)
    out = graph(s(1), s(2), 3);
end
function out = t(s)
    out = graph(s(1), s(2), 3);
end
function sett(s, val)
    graph(s(1), s(2), 3) = val;
end
function setQ(s)
     graph(s(1), s(2), 3) = 1;
end
function rsetQ(s)
     graph(s(1), s(2), 3) = 0;
end
function incr(s)
    graph(s(1), s(2), 4) = graph(s(1), s(2), 4)+1;
end
%-----------------------------------------------------------
function out = g(s)
        s = abs(startPos - s);
%         s = abs(endPos - s);
%         out = SQRT2*min(s(1:2)) + max(s(1:2));
        out = sqrt(s(1)^2 + s(2)^2);
end
function out = h(s)
    out = graph(s(1), s(2),1);
end
function seth(s, val)
    graph(s(1), s(2),1) = val;
end
%-----------------------------------------------------------
function out = k(s)
    out = graph(s(1), s(2), 2);
end
function setk(s, val)
    graph(s(1), s(2),2) = val;
end
%-----------------------------------------------------------
function out = calculateKey(x)
    out =  [
        h(x)+g(x)+kM;
%           min( h(x), k(x)+g(x));
          min(h(x), k(x))
        ];
end
function out = testNode(s)
    out = true;
    if template(s(1), s(2)) == 1
        out = false;
        return
    end
    for n = pattern
        pos = s + n;
        if map(pos(1), pos(2)) == 0
            out = false;
            template(s(1), s(2)) = 1;
%             graph(s(1), s(1), 5) = 0;
%             setg(s, -1)
            return
        end
    end
end
function out = cmp(s1, s2)
    out = s1(1) < s2(1) || ( s1(1) ==s2(1) && s1(2) < s2(2));
end
function insert(s, h_new)
        t = inQ(s);
        if t == NEW
            setk(s, h_new);
            seth(s, h_new);
        else
            if t == OPEN
                setk(s, min(k(s), h_new) );
            else
                setk(s, min(h(s), h_new) );
                seth(s, h_new);
            end
            s(3:4)  = [k(s)+kM+g(s); k(s)+g(s)];
            add(stack, s);
            sett(s, OPEN);
        end

end

end

author: LazyFalcon


text:
function cc=lpcar2cc(ar,np)
%LPCAR2CC LPC: Convert ar filter to complex cepstrum CC=(AR,NP)
% the "real" cepstrum is half the complex cepstrum
% cc() does not include c0 whose value can be calculated
% from the prediction residual energy, e, as ln(e)/2
% for both real and complex cepstrum.

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

[nf,p1]=size(ar);
p=p1-1;
if (nargin<2) np=p; end
cc=zeros(nf,np);
cm=(1:np).^(-1);
if np>p
  xm=-(1:p);
  nz=np-p;
  for k=1:nf
    cc(k,:)=filter(1,ar(k,:),[ar(k,2:p1).*xm zeros(1,nz)]).*cm;
  end
else
  p1=np+1;
  xm=-(1:np);
  for k=1:nf
    cc(k,:)=filter(1,ar(k,:),ar(k,2:p1).*xm).*cm;
  end
end

author: decouples


text:
function mm_to_st ( input_filename, output_filename )

%*****************************************************************************80
%
%% MSM_TO_ST converts sparse matrix information from MM to ST format.
%
%  Discussion:
%
%    MM format is the "Matrix Market" format.
%    ST format is the "Sparse Triplet" format.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the MM file.
%
%    Input, string OUTPUT_FILENAME, the name of the ST file.
%
  [ a, rows, cols, entries, rep, field, symm ] = mm_to_msm ( input_filename );

  msm_to_st ( a, output_filename );

  return
end


text:
function D=sim2dis(S)
%function D=sim2dis(S);
%
%PURPOSE
%
%To transform elements of a similarity matrix S into
%dissimilarities D by D=1-S.
%
%INPUT
%
%S (matrix) MxM similarity matrix whose elements must be in [-1,1],
%    for example, absolute values of linear correlation coefficients.
%
%OUTPUT
%
% D  (matrix) MxM dissimilarity matrix
%
%SEE ALSO
%  sim2dis2
%
%USED IN
%  icassoCluster
%  icassoProjection

%COPYRIGHT NOTICE
%This function is a part of Icasso software library
%Copyright (C) 2003-2005 Johan Himberg
%
%This program is free software; you can redistribute it and/or
%modify it under the terms of the GNU General Public License
%as published by the Free Software Foundation; either version 2
%of the License, or 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, write to the Free Software
%Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

% ver 1.2 johan 100105

if any(S(:)>1) | any(S(:)<-1),
  error('Values of similarity must be in [-1,1]');
end

D=1-S;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function nv=nodesurfnorm(node,elem)
%  nv=nodesurfnorm(node,elem)
%
%  calculate a nodal norm for each vertix on a surface mesh (surface
%   can only be triangular or cubic)
%
%  author: Qianqian Fang <q.fang at neu.edu>
%  date: 12/12/2008
%
% parameters:
%      node: node coordinate of the surface mesh (nn x 3)
%      elem: element list of the surface mesh (3 columns for
%            triangular mesh, 4 columns for cubic surface mesh)
%      pt: points to be projected, 3 columns for x,y and z respectively
%
% outputs:
%      nv: nodal norms (vector) calculated from nodesurfnorm.m
%          with dimensions of (size(v,1),3)
%
% Please find more information at http://iso2mesh.sf.net/cgi-bin/index.cgi?metch
%
% this function is part of "metch" toobox, see COPYING for license

nn=size(node,1);
ne=size(elem,1);
nedim=size(elem,2);

ev=surfacenorm(node,elem);

nv=zeros(nn,3);
ev2=repmat(ev,1,3);
for i=1:ne
  nv(elem(i,:),:)=nv(elem(i,:),:)+reshape(ev2(i,:),3,3)';
end
nvnorm=sqrt(sum(nv.*nv,2));
idx=find(nvnorm>0);
if(length(idx)<nn)
	warning(['found interior nodes, their norms will be set to zeros; to remove ',...
                'them, please use removeisolatednodes.m from iso2mesh toolbox']);

	nv(idx,:)=nv(idx,:)./repmat(nvnorm(idx),1,3);
else
	nv=nv./repmat(nvnorm,1,3);
end

author: fangq


text:
function varargout = ismember(varargin)
% (1) To make Octave's `ismember` behaviour match recent Matlab's: the second
% output gives the index of the first occurrence, not the last
% (2) Octave gives error when the first two inputs are char and
% numeric, whereas Matlab works. So I include here conversion from char to
% double in that case.
% (3) Octave's ismember treats complex inputs as their absolute value,
% which is incorrect (see https://savannah.gnu.org/bugs/index.php?52437).
% This seems to have been caused by a "reimplementation" of ismember. So if
% any input is complex we use unique to replace the inputs by unique labels
% (3):
if (isnumeric(varargin{1})&&~isreal(varargin{1})) || (isnumeric(varargin{2})&&~isreal(varargin{2}))
    a = varargin{1};
    b = varargin{2};
    [~, ~, u] = unique([a(:); b(:)]);
    a(:) = real(u(1:numel(a)));
    b(:) = real(u(numel(a)+1:end));
    varargin{1} = a;
    varargin{2} = b;
end
% (2):
if ischar(varargin{1})&&isnumeric(varargin{2})
   varargin{1} = double(varargin{1});
end
if isnumeric(varargin{1})&&ischar(varargin{2})
   varargin{2} = double(varargin{2});
end
% (1):
if nargout==1 % single output: no interception needed
    varargout{1} = builtin('ismember', varargin{:});
elseif nargin>=3 && strcmp(varargin{3},'rows') % two outputs, 'rows' case
    [aim, bim] = builtin('ismember', varargin{1}, varargin{2}(end:-1:1,:), 'rows');
    bim(aim) = size(varargin{2},1)+1-bim(aim);
    varargout = {aim, bim};
else % two outputs, no 'rows'
    [aim, bim] = builtin('ismember', varargin{1}, varargin{2}(end:-1:1));
    bim(aim) = numel(varargin{2})+1-bim(aim);
    varargout = {aim, bim};
end
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
function res = im2row(im, winSize)
%res = im2row(im, winSize)
[sx,sy,sz] = size(im);

res = zeros((sx-winSize(1)+1)*(sy-winSize(2)+1),prod(winSize),sz);
count=0;
for y=1:winSize(2)
    for x=1:winSize(1)
        count = count+1;
        res(:,count,:) = reshape(im(x:sx-winSize(1)+x,y:sy-winSize(2)+y,:),...
            (sx-winSize(1)+1)*(sy-winSize(2)+1),1,sz);
    end
end

author: thomaskuestner


text:
function [err,XYZdic] = AFNI_Index2XYZcontinuous (Indx, Info, CoordCode)
%
%   [err,XYZdic] = AFNI_Index2XYZcontinuous (Indx, Info, [CoordCode])
%
%Purpose:
%   Change from voxel XYZindex (called Voxel Coords in AFNI) to XYZ in mm
%   The mm and voxel coordinates refer to the values displayed
%   on the top left corner of AFNI controller.
%   CoordCode is the one you'd set from the Coord Order plugin
%
%
%Input Parameters:
%   Indx an Mx3 matrix or an  Mx1 vector containing the voxel indices to be
%        transformed to voxel coordinates.  (indices start at 0)
%   Info is the output of BrikInfo
%   CoordCode is an optional parameter used to specify the coordinates system of the output
%      if empty or not specified, the default is 'RAI'. The code can be either a string or a vector
%      of numbers (see AFNI_CoordChange for more on that)
%
%Output Parameters:
%   err : 0 No Problem
%       : 1 Mucho Problems
%   XYZdic : The continuous coordinates corresponding to Indx
%       The coordnate system output is in RAI (DICOM)
%       unless otherwise specified by CoordCode
%
%
%Key Terms:
%
%More Info :
%   BrikInfo
%   Test_AFNI_Index2XYZcontinuous
%   AFNI_XYZcontinuous2Index
%   Test_AFNI_XYZcontinuous2Index
%
% You can also go from index to XYZ using the header field IJK_TO_DICOM_REAL
%  For instance, say you have voxel indices 12, 2, 4 (matlab indexing 13, 3, 5)
%  to go from AFNI index to AFNI DICOM RAI you can do:
%     M = [reshape(Info.IJK_TO_DICOM_REAL, 4, 3)' ; 0 0 0 1];
%     I = [12 2 4 1]';
%     X = M*I;
%     To go from AFNI DICOM RAI to AFNI indices:
%     I = inv(M)*X;
%
%     Author : Ziad Saad
%     Date : Tue Sep 5 21:48:06 PDT 2000           Latest Modification: Feb 18 04
%     LBC/NIMH/ National Institutes of Health, Bethesda Maryland

%Define the function name for easy referencing
FuncName = 'AFNI_Index2XYZcontinuous';

%Debug Flag
DBG = 1;

ChangeCoord = 0;
if (nargin == 3)
	if (~isempty(CoordCode)),
		ChangeCoord = 1;
	end
end

%initailize return variables
err = 1;
XYZmm = [];

%make sure Indx is the right size
switch size(Indx,2),
	case 1, %change 1D index to XYZ index
		[err, Indx] = AfniIndex2AfniXYZ (Indx, Info.DATASET_DIMENSIONS(1), Info.DATASET_DIMENSIONS(2))
	case 3, %OK
	otherwise,
		err = ErrEval(FuncName,'Err_Bad dimension for Indx');
		return
end

XYZmm = Indx;

	%The equations that would change the indices to coordinate system result in a coordinate system that
	% may be any permutation of RAI (like IRA or AIR or IAR or RIA or ARI) so one only needs to find the
	%dimension permutation needed to bring the final result to RAI.

	%determine the ordering map to go from any permutation of RAI to RAI
		%[maploc(1),jnk] = find(Info.Orientation == 'R');
		%[maploc(2),jnk] = find(Info.Orientation == 'A');
		%[maploc(3),jnk] = find(Info.Orientation == 'I');

	%pre - Wed May 23 18:20:56 PDT 2001 - WRONG !
		%XYZmm(:, maploc(1)) = Info.ORIGIN(1) + Indx(:,1) .* Info.DELTA(1);
		%XYZmm(:, maploc(2)) = Info.ORIGIN(2) + Indx(:,2) .* Info.DELTA(2);
		%XYZmm(:, maploc(3)) = Info.ORIGIN(3) + Indx(:,3) .* Info.DELTA(3);

	%post - Wed May 23 18:20:56 PDT 2001 - WRONG!
		%XYZmm(:, 1) = Info.ORIGIN(maploc(1)) + Indx(:,maploc(1)) .* Info.DELTA(maploc(1));
		%XYZmm(:, 2) = Info.ORIGIN(maploc(2)) + Indx(:,maploc(2)) .* Info.DELTA(maploc(2));
		%XYZmm(:, 3) = Info.ORIGIN(maploc(3)) + Indx(:,maploc(3)) .* Info.DELTA(maploc(3));

   %Feb 18 04, back to the original
      XYZmm(:, 1) = Info.ORIGIN(1) + Indx(:,1) .* Info.DELTA(1);
		XYZmm(:, 2) = Info.ORIGIN(2) + Indx(:,2) .* Info.DELTA(2);
		XYZmm(:, 3) = Info.ORIGIN(3) + Indx(:,3) .* Info.DELTA(3);
      %Now this is in the axis orientation which is Info.Orientation(:,1)' called 3dmm in thd_coords.c
      [err,XYZdic, map] = THD_3dmm_to_dicomm (Info, XYZmm);

if (ChangeCoord),
	[err, maplocation, mapsign, XYZdic] = AFNI_CoordChange ('RAI', CoordCode, XYZdic);
end

err = 0;
return;


text:
		% Generarea retelei de noduri pe domeniului considerat
[X,Y]=meshgrid(-3.75:0.05:3.5);
		% Evaluarea functiei de doua variabile de reprezentat
Z=peaks(X,Y);
		% Incarcarea nivelelor negative de reprezentat
v_neg=-3:0.5:-1;
		% Incarcarea nivelelor pozitive de reprezentat
v_poz=1:0.5:7;
		% Desenarea liniilor de contur corespunzatoare
		%  - nivelelor negative cu linii intrerupte albastre
contour(X,Y,Z,v_neg,'b--');
hold on
		%  - nivelelor pozitive cu linii continue rosii
contour(X,Y,Z,v_poz,'r-');
hold off


text:
%% "Spirograph" Code
%  Jonathan Jamieson - September 2013
%  unigamer@gmail.com
%  www.jonathanjamieson.com

function [P_x_global,P_y_global ] = CalculatePenPosition(A_x_global,A_y_global,C_x_global,C_y_global,extension  )

dx = C_x_global-A_x_global;
dy = C_y_global-A_y_global;

theta = atan2(dy,dx);

P_x_global = extension*cos(theta)+C_x_global;
P_y_global = extension*sin(theta)+C_y_global;

end


text:
% Question No: 7

% Given an image. Implement the split and merge procedure for segmenting
% the image with different values for minimum dimensions of the quadtree
% regions.

function quadtree(x)
f=imread(x);
q=2^nextpow2(max(size(f)));
[m n]=size(f);
f=padarray(f,[q-m,q-n],'post');
mindim=2;
s=qtdecomp(f,@split,mindim,@predicate);
lmax=full(max(s(:)));
g=zeros(size(f));
marker=zeros(size(f));
for k=1:lmax
    [vals,r,c]=qtgetblk(f,s,k);
    if ~isempty(vals)
        for i=1:length(r)
            xlow=r(i);ylow=c(i);
            xhigh=xlow+k-1;
            yhigh=ylow+k-1;
            region=f(xlow:xhigh,ylow:yhigh);
            flag=feval(@predicate,region);
            if flag
                 g(xlow:xhigh,ylow:yhigh)=1;
                 marker(xlow,ylow)=1;
            end
        end
    end
end
g=bwlabel(imreconstruct(marker,g));
g=g(1:m,1:n);
f=f(1:m,1:n);
figure, imshow(f),title('Original Image');
figure, imshow(g),title('Segmented Image');
end

function v=split(b,mindim,fun)
k=size(b,3);
v(1:k)=false;
for i=1:k
    quadrgn=b(:,:,i);
    if size(quadrgn,1)<=mindim
        v(i)=false;
        continue;
    end
    flag=feval(fun,quadrgn);
    if flag
        v(i)=true;
    end
end
end

function flag=predicate(region)
sd=std2(region);
m=mean2(region);
flag=(sd>5)&(m>0)&(m<200);
end


text:
function [nconds cforecast_record cforecast_estimates]=panel3cf(N,n,m,p,k,q,cfconds,cfshocks,cfblocks,data_endo_a,data_exo_a,data_exo_p,It,Bu,Fperiods,const,beta_gibbs,D_record,gamma_record,CFt,Fband)

% initiate the cell recording the Gibbs sampler draws
cforecast_record={};
cforecast_estimates={};

% because conditional forecasts can be computed for many (potentially all) units, loop over units
for ii=1:N
% check wether there are any conditions on unit ii
temp=cfconds(:,:,ii);
nconds(ii,1)=numel(temp(cellfun(@(x) any(~isempty(x)),temp)));

   % if there are conditions
   if nconds(ii,1)~=0
   % prepare the elements for conditional forecast estimation, depending on the type of conditional forecasts
   temp1=cfconds(:,:,ii);
      if CFt==1
      temp2={};
      temp3=[];
      elseif CFt==2
      temp2=cfshocks(:,:,ii);
      temp3=cfblocks(:,:,ii);
      end
   % run the Gibbs sampler for unit ii
   cforecast_record(:,:,ii)=bear.cforecast(data_endo_a(:,:,ii),data_exo_a,data_exo_p,It,Bu,Fperiods,temp1,temp2,temp3,CFt,const,beta_gibbs(:,:,ii),D_record(:,:,ii),gamma_record(:,:,ii),n,m,p,k,q);
   % then obtain point estimates and credibility intervals
   cforecast_estimates(:,:,ii)=bear.festimates(cforecast_record(:,:,ii),n,Fperiods,Fband);

   % if there are no conditions, return empty elements
   elseif nconds(ii,1)==0
   cforecast_record(:,:,ii)=cell(n,1);
   cforecast_estimates(:,:,ii)=cell(n,1);
   end
end

author: european-central-bank

repo: BEAR-toolbox


text:
function y=jinc(r)
% y=jinc(r) returns the jinc function [J1(2*pi*r)/(2*pi*r)] evaluated
% at r. Per this definition, the first zero of jinc function is at 0.61.

y=real(2*besselj(1,2*pi*r)./ (2*pi*r));
%Imaginary parts due to numerical error are of order 1E-16.

%Division by zero leads to NaN.
y(isnan(y))=1;

end


text:
function beale_test ( )

%*****************************************************************************80
%
%% BEALE_TEST works with the Beale function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BEALE_TEST:\n' );
  fprintf ( 1, '  Test COMPASS_SEARCH with the Beale function.\n' );
  m = 2;
  delta_tol = 0.00001;
  delta = 0.1;
  k_max = 20000;

  x = [ 1.0, 1.0 ];
  r8vec_print ( m, x, '  Initial point X0:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X0) = %g\n', beale ( m, x ) );

  [ x, fx, k ] = compass_search ( @beale, m, x, delta_tol, delta, k_max );
  r8vec_print ( m, x, '  Estimated minimizer X1:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X1) = %g, number of steps = %d\n', fx, k );
%
%  Repeat with more difficult start.
%
  x = [ 1.0, 4.0 ];
  r8vec_print ( m, x, '  Initial point X0:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X0) = %g\n', beale ( m, x ) );

  [ x, fx, k ] = compass_search ( @beale, m, x, delta_tol, delta, k_max );
  r8vec_print ( m, x, '  Estimated minimizer X1:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X1) = %g, number of steps = %d\n', fx, k );
%
%  Demonstrate correct minimizer.
%
  x = [ 3.0, 0.5 ];
  r8vec_print ( m, x, '  Correct minimizer X*:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X*) = %g\n', beale ( m, x ) );

  return
end


text:
function spm_dem_search_movie(DEM)
% creates a movie of visual search in extrinsic and intrinsic coordinates
% FORMAT spm_dem_search_movie(DEM)
%
% DEM - {DEM} structures from visual search simulations
%
% hidden causes and states
%==========================================================================
% x    - hidden states:
%   o(1) - oculomotor angle
%   o(2) - oculomotor angle
%   x(1) - relative amplitude of visual hypothesis 1
%   x(2) - relative amplitude of visual hypothesis 2
%   x(3) - ...
%
% v    - hidden causes
%
% g    - sensations:
%   g(1) - oculomotor angle (proprioception - x)
%   g(2) - oculomotor angle (proprioception - y)
%   g(3) - retinal input - channel 1
%   g(4) - retinal input - channel 2
%   g(5) - ...
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dem_search_movie.m 6901 2016-10-08 13:21:41Z karl $

% Preliminaries
%--------------------------------------------------------------------------
clf, global STIM
if ~iscell(DEM), DEM = {DEM};           end
if ~isfield(STIM,'W'), STIM.W = 1/6;    end
if ~isfield(STIM,'P'), STIM.P = [0;0];  end
if ~isfield(STIM,'U'), STIM.U = STIM.V; end
if ~isfield(STIM,'S'), STIM.S = STIM.H; end

N   = length(DEM);
S   = spm_read_vols(STIM.U);

% Stimulus
%==========================================================================
Dx  = STIM.P(1)*16 + STIM.U.dim(1)/2;
Dy  = STIM.P(2)*16 + STIM.U.dim(2)/2;

dim = size(STIM.R);
vox = STIM.U.dim(1);
dx  = vox/dim(1)*STIM.W;
di  = dx*([1 dim(1)] - dim(1)/2) + Dx;
dj  = dx*([1 dim(2)] - dim(2)/2) + Dy;
di  = [di;di]; di = di(:);
dj  = [dj;dj]';dj = dj(:);
ax  = [(Dy - vox) (Dy + vox) (Dx - vox) (Dx + vox)]/2;

a     = [];
for i = 1:N

    % i-th saccade - position
    %----------------------------------------------------------------------
    pU = DEM{i}.pU.x{1}(1:2,:)*16;
    qU = DEM{i}.qU.x{1}(1:2,:)*16;
    T  = length(pU);
    a  = [a DEM{i}.qU.a{2}];

    % eye movements in extrinsic coordinates
    %======================================================================
    subplot(2,2,1)
    for t = 1:T
        image((S + 1)*32), axis image, hold on
        plot(qU(2,t) + Dy,qU(1,t) + Dx,'.g','Markersize',8)
        plot(pU(2,t) + Dy,pU(1,t) + Dx,'.r','Markersize',16)
        plot(pU(2,t) + dj,pU(1,t) + di,'+' ,'Markersize',16)

        % show location of image ccentre
        %------------------------------------------------------------------
        plot(Dy,Dx,'+r','Markersize',32);  % axis(ax);

        drawnow, hold off

        % save
        %------------------------------------------------------------------
        Me(t + T*(i - 1)) = getframe(gca);

    end

    % sensory input
    %======================================================================
    subplot(2,2,2)
    for t = 1:T

        o   = DEM{i}.pU.x{1}(:,t);
        s   = ADEM_sample_image(STIM.U,o,STIM.R);
        image(s*64), axis image, drawnow

        % save
        %----------------------------------------------------------------------
        Mi(t + T*(i - 1)) = getframe(gca);

    end

    % i-th saccade - percept
    %----------------------------------------------------------------------
    qU = DEM{i}.qU.x{1}(3:end,:);

    % percept
    %======================================================================
    subplot(2,2,4)
    for t = 1:T

        % hypotheses
        %------------------------------------------------------------------
        h     = spm_softmax(qU(:,t),2);
        H     = 1 + h'*log(h)/log(length(h));

        % retinotopic predictions
        %------------------------------------------------------------------
        s     = 0;
        for j = 1:length(h)
            s = s + h(j)*spm_read_vols(STIM.S{j});
        end
        image(s*H*64), axis image, drawnow

        % save
        %------------------------------------------------------------------
        Mq(t + T*(i - 1)) = getframe(gca);

    end

end

% set ButtonDownFcn
%--------------------------------------------------------------------------
subplot(2,2,3)
plot(a')
title('Action (EOG)','FontSize',16)
xlabel('time')
axis square

subplot(2,2,1)
set(gca,'Userdata',{Me,16})
set(gca,'ButtonDownFcn','spm_DEM_ButtonDownFcn')
title('saccades (click axis for movie)','FontSize',16)

subplot(2,2,2)
set(gca,'Userdata',{Mi,16})
set(gca,'ButtonDownFcn','spm_DEM_ButtonDownFcn')
title('samples (click axis for movie)','FontSize',16)

subplot(2,2,4)
set(gca,'Userdata',{Mq,16})
set(gca,'ButtonDownFcn','spm_DEM_ButtonDownFcn')
title('percept (click axis for movie)','FontSize',16)

author: spm


text:
function h = int_hist(x, n)
% INT_HIST(x, n) is a histogram of all integer values 1..n in x.
% If n is not given, max(x) is used.

% Hans Olsson's one-liner from matlab faq
h = full(sum(sparse(1:length(x(:)),x(:),1)));
if nargin == 2
  if n > length(h)
    % pad with zeros
    h = [h zeros(1,n-length(h))];
  end
end

author: andrewssobral

repo: lrslibrary