text:
function dist = circle_imp_point_dist_2d ( r, center, p )

%*****************************************************************************80
%
%% CIRCLE_IMP_POINT_DIST_2D: distance ( implicit circle, point ) in 2D.
%
%  Discussion:
%
%    The distance is zero if the point is on the circle.
%
%    An implicit circle in 2D satisfies the equation:
%
%      ( X - CENTER(1) )**2 + ( Y - CENTER(2) )**2 = R**2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%
%    Input, real CENTER(2,1), the center of the circle.
%
%    Input, real P(2,1), the point to be checked.
%
%    Output, real DIST, the distance of the point to the circle.
%
  dim_num = 2;

  r2 = sqrt ( sum ( ( p(1:2,1) - center(1:2,1) ).^2 ) );

  dist = abs ( r2 - r );

  return
end


text:
function tetrahedron_monte_carlo ( )

%*****************************************************************************80
%
%% TETRAHEDRON_MONTE_CARLO_TEST tests the TETRAHEDRON_MONTE_CARLO library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 January 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TETRAHEDRON_MONTE_CARLO_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the TETRAHEDRON_MONTE_CARLO library.\n' );

  tetrahedron_monte_carlo_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TETRAHEDRON_MONTE_CARLO_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Offspring = Operator(Problem,Particle,Pbest,Gbest,Parameter)
% Particle swarm optimization in MPSO/D
% c1   ---   2 --- Parameter in updating particle's velocity
% c2   ---   2 --- Parameter in updating particle's velocity
% CR   --- 0.5 --- Parameter CR in differental evolution
% F    --- 0.5 --- Parameter F in differental evolution
% proM ---   1 --- The expectation of number of bits doing mutation
% disM ---  20 --- The distribution index of polynomial mutation

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Parameter setting
    if nargin > 4
        [c1,c2,CR,F,proM,disM] = deal(Parameter{:});
    else
        [c1,c2,CR,F,proM,disM] = deal(2,2,0.5,0.5,1,20);
    end
    ParticleDec = Particle.decs;
    PbestDec    = Pbest.decs;
    GbestDec    = Gbest.decs;
    [N,D]       = size(ParticleDec);
    ParticleVel = Particle.adds(zeros(N,D));

    %% Particle swarm optimization
    Lower = repmat(Problem.lower,N,1);
    Upper = repmat(Problem.upper,N,1);
    DoPSO = repmat(rand(N,1)<0.5,1,D);
    W     = 0.9 - Problem.FE./Problem.maxFE*0.8;
    r1    = repmat(rand(N,1),1,D);
    r2    = repmat(rand(N,1),1,D);
    OffVel        = ParticleVel;
    OffDec        = ParticleDec;
    OffVel(DoPSO) = W.*ParticleVel(DoPSO) + c1.*r1(DoPSO).*(PbestDec(DoPSO)-ParticleDec(DoPSO)) + c2.*r2(DoPSO).*(GbestDec(DoPSO)-ParticleDec(DoPSO));
    OffDec(DoPSO) = ParticleDec(DoPSO) + OffVel(DoPSO);
    % Set the infeasible decision variables to the value of their parents
    Invalid         = OffDec < Lower | OffDec > Upper;
    OffDec(Invalid) = ParticleDec(Invalid);

    %% DE
    Site = ~DoPSO & rand(N,D)<CR;
    OffDec(Site) = ParticleDec(Site) + F.*(GbestDec(Site)-PbestDec(Site));
    % Set the infeasible decision variables to boundary values
    OffDec = max(min(OffDec,Upper),Lower);

    %% Polynomial mutation
    Site  = rand(N,D) < proM/D;
    mu    = rand(N,D);
    temp  = Site & mu<=0.5;
    OffDec(temp) = OffDec(temp)+(Upper(temp)-Lower(temp)).*((2.*mu(temp)+(1-2.*mu(temp)).*...
                   (1-(OffDec(temp)-Lower(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1))-1);
    temp = Site & mu>0.5;
    OffDec(temp) = OffDec(temp)+(Upper(temp)-Lower(temp)).*(1-(2.*(1-mu(temp))+2.*(mu(temp)-0.5).*...
                   (1-(Upper(temp)-OffDec(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1)));
    Offspring = Problem.Evaluation(OffDec,OffVel);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%% data
clc;
clear all;
close all;
snri=(-10:2:30);
qam=[2,4,8,16,32,64];
capacity = zeros(length(qam),length(snri));
for qami=1:length(qam)
    for index=1:length(snri)
        capacity(qami,index) = QAMCapacity(snri(index),1,qam(qami));
    end
end
figure;
hold
for qami=1:length(qam)
    plot(snri,capacity(qami,:),'LineWidth',1.6);
end
GaussianC = zeros(1,length(snri));
for index=1:length(snri)
    GaussianC(index) = log2(1+ 10^(snri(index)/10));
end
plot(snri,GaussianC,'LineWidth',1.6);
set(gca,'FontSize', 12, 'FontName', 'Times New Roman');
grid on;
box on;
xlabel('SNR (Es/No) dB','FontSize', 14, 'FontName', 'Times New Roman');
ylabel('Capacity (bits/Tx)','FontSize',14,'FontName', 'Times New Roman');

save ('DiscretInContinuousOutCapacityAWGN.mat','snri','capacity','GaussianC')

author: Sable

repo: mcbench-benchmarks


text:
classdef BIF < handle
    %BIF  Implementation of bio-inspired features (BIF)
    %
    % An implementation of the bio-inspired features (BIF) approach for
    % computing image descriptors, applicable for human age estimation. For
    % more details we refer to [BIF09] and [BIF2015].
    %
    % ## References
    % [BIF09]:
    % > Guo, Guodong, et al. "Human age estimation using bio-inspired
    % > features". Computer Vision and Pattern Recognition, 2009. CVPR 2009.
    %
    % [BIF2015]:
    % > Spizhevoi, A. S., and A. V. Bovyrin. "Estimating human age using
    % > bio-inspired features and the ranking method." Pattern Recognition and
    % > Image Analysis 25.3 (2015): 547-552.
    %
    % See also: cv.BIF.BIF, cv.BIF.compute
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent, SetAccess = private)
        % The number of filter bands used for computing BIF
        NumBands
        % The number of image rotations
        NumRotations
    end

    %% Constructor/destructor
    methods
        function this = BIF(varargin)
            %BIF  Constructor
            %
            %     obj = cv.BIF()
            %     obj = cv.BIF('OptionName',optionValue, ...)
            %
            % ## Options
            % * __NumBands__ The number of filter bands (`<= 8`) used for
            %   computing BIF. default 8
            % * __NumRotations__ The number of image rotations for computing
            %   BIF. default 12
            %
            % See also: cv.BIF.compute
            %
            this.id = BIF_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.BIF
            %
            if isempty(this.id), return; end
            BIF_(this.id, 'delete');
        end
    end

    %% BIF
    methods
        function features = compute(this, img)
            %COMPUTE  Computes features sby input image
            %
            %     features = model.compute(img)
            %
            % ## Input
            % * __image__ Input image (1-channel `single`).
            %
            % ## Output
            % * __features__ Feature vector (`single` type).
            %
            % See also: cv.BIF.BIF
            %
            features = BIF_(this.id, 'compute', img);
        end
    end

    %% Algorithm
    methods (Hidden)
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.BIF.empty, cv.BIF.load
            %
            BIF_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Checks if detector object is empty.
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the detector object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.BIF.clear, cv.BIF.load
            %
            b = BIF_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.BIF.load
            %
            BIF_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.BIF.save
            %
            BIF_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.BIF.save, cv.BIF.load
            %
            name = BIF_(this.id, 'getDefaultName');
        end
    end

    %% Getters/Setters
    methods
        function value = get.NumBands(this)
            value = BIF_(this.id, 'get', 'NumBands');
        end

        function value = get.NumRotations(this)
            value = BIF_(this.id, 'get', 'NumRotations');
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function ppgOnsets = qppg_fast(data,fs,from,to)

% This function is rewriten from wabp_pleth_new.c and wabp.c
% /* file wabp.c          Wei Zong       23 October 1998
%    			Last revised:   9 April 2010 (by G. Moody)
% -----------------------------------------------------------------------------
% wabp: beat detector for arterial blood presure (ABP) signal
% Copyright (C) 1998-2010 Wei Zong
%
% 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 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.
%
% You may contact the author by e-mail (wzong@mit.edu) or postal mail
% (MIT Room E25-505, Cambridge, MA 02139, USA).  For updates to this software,
% please visit PhysioNet (http://www.physionet.org/).
% ------------------------------------------------------------------------------
%
% This program detects heart beats (pulse waveforms) in a continuous arterial
% blood pressure (ABP) signal.  This version of wabp works best with ABP signals
% sampled at 125 Hz, but it can analyze ABPs sampled at any frequency
% using on-the-fly resampling provided by the WFDB library.  'wabp' has been
% optimized for adult human  ABPs. For other ABPs, it may be necessary to
% experiment with the input sampling frequency and the time constants indicated
% below.
%
% `wabp' can process records containing any number of signals, but it uses only
% one signal for ABP pulse detection (by default, the lowest-numbered signal
% labelled `ABP', `ART', or `BP';  this can be changed using the `-s' option, see
% below).
%
% To compile this program under GNU/Linux, MacOS/X, MS-Windows, or Unix, use gcc:
%         gcc -o wabp wabp.c -lwfdb
% You must have installed the WFDB library, available at
%         http://www.physionet.org/physiotools/wfdb.shtml
% gcc is standard with GNU/Linux and is available for other platforms from:
%         http://www.gnu.org/             (sources and Unix binaries)
%         http://fink.sourceforge.net     (Mac OS/X only)
%         http://www.cygwin.com/          (MS-Windows only)
%
% For a usage summary, see the help text at the end of this file.  The input
% record may be in any of the formats readable by the WFDB library, and it may
% be anywhere in the WFDB path (in a local directory or on a remote web or ftp
% server).  The output of 'wabp' is an annotation file named RECORD.wabp (where
% RECORD is replaced by the name of the input record).  Within the output
% annotation file, the time of each NORMAL annotation marks an ABP pulse wave
% onset.
%
%
%
% input :   data: PPG data
%           fs:   sampling frequency, default 125 Hz
%           from: begin point to analysis
%           to  : end point to analysis
% output:   ppgOnsets: onset position of PPG beats in samples
%
%
%
% CHANGE LOG (should be moved to GitHub as commit messages...)
%
% 03 Aug 2010:  by Qiao Li
%   1) sample() function sometimes get WFDB_INVALID_SAMPLE value even when data is good.
%   2) Changed to mysample(), using isigsettime() and getvec() instead, but slow!!!
%   3) Add a mysetbuf() function, read the data to a databuf first, then mysample()
%      get data from the buf. It's much more fast than before.
%
%
% 04 Jul 2011 by Qiao Li
%   1) Set eye-closing period for PPG 0.34s, slope width 0.17s
%   2) Remove physical units to ADC units transform //Tm = physadu((unsigned)sig, Tm);
%   3) Add maxmin_2_3_threshold to modify the logic of eye-closing in order to minimize
%      the double beats
%	4) Change t += EyeClosing; to t = tpq+EyeClosing;
%
%
% 19 Jul 2011 by Qiao Li
%
%   1) add: (before reading data,at line 480)
%      isigsettime(from);
%      dataL=from;
%      to read data from 'from' setting
% 	2) Changed: (after learning period, at line 544)
%      (void)sample(sig, tpq);
%      if (sample_valid() == 0) break;
%       to:
%      if (dataend) break;
%
% 18 Oct 2016 by Qiao Li
%   1) add: input parameter fs for different sampling frequency data
%	2) add: re-scale data to ~ +/- 2000
%   3) add: find valley from the original data around 0.25s of tpq
%
% 03 Mar 2017 by Adriana Vest
%   Changed name of function to qppg to avoid confusion with wabp
%	Previous name: wabp_pleth_new.m
%
% 12 Sep 2018 by Giulia Da Poian
%   Changed output variable name to ppgOnsets
%
% 22 Dec 2018 by Giulia Da Poian
%   Start from the qppg.m code, replace use of global variables and use
%   function arguments, rename the function qppg_fast

if nargin<3
    from=1;
    to=length(data);
end

if nargin<2
    error('Wrong Number of Input Arguments: Sampling Frequency is required')
end

ppgOnsets=[];
beat_n=1;

sps=fs; % Sampling Frequency

BUFLN = 4096;           % /* must be a power of 2, see slpsamp() */
EYE_CLS = 0.34;         % /* eye-closing period is set to 0.34 sec (340 ms) for PPG */
LPERIOD  = sps*8;       % /* learning period is the first LPERIOD samples */
SLPW = 0.17;            % /* Slope width (170ms) for PPG */
NDP = 2.5;              % /* adjust threshold if no pulse found in NDP seconds */
TmDEF = 5;              % /* minimum threshold value (default) */
Tm = TmDEF;

BUFLN2 = (BUFLN*2);

INVALID_DATA=-32768;
if data(1)<=INVALID_DATA+10
    data(1)=mean(data);
end
inv=find(data<=INVALID_DATA+10);
for i=1:length(inv)
    data(inv(i))=data(inv(i)-1);
end

% re-scale data to ~ +/- 2000
if length(data)<5*60*sps
    data=(data-min(data))./(max(data)-min(data)).*4000-2000;
else
% find max/min every 5 minute for re-scaling data
    n=1;
    for i=1:5*60*sps:length(data)
        max_data(n)=max(data(i:min(i+5*60*sps-1,length(data))));
        min_data(n)=min(data(i:min(i+5*60*sps-1,length(data))));
        n=n+1;
    end
    data=(data-median(min_data))./(median(max_data)-median(min_data)).*4000-2000;
end

samplingInterval = 1000.0/sps;
spm = 60 * sps;
EyeClosing = round(sps * EYE_CLS);   % /* set eye-closing period */
ExpectPeriod = round(sps * NDP);	  % /* maximum expected RR interval */
SLPwindow = round(sps * SLPW);       % /* slope window size */
timer=0;

ebuf(1:BUFLN)=0;
lbuf=ebuf;
if from>BUFLN
    tt_2=from-BUFLN;
else
    tt_2=0;
end
aet=0;

t1=8*sps;
t1 = t1+from;
T0 = 0;
n=0;
for t = from:t1
   [temp,ebuf,lbuf,tt_2, aet] = slpsamp(t,data,BUFLN,ebuf,lbuf,tt_2, aet,SLPwindow);
    if temp > INVALID_DATA+10
        T0 = T0+temp;
        n=n+1;
    end
end
T0 = T0/n; % T0=T0/(t1-from);
Ta = 3 * T0;

learning=1;

%    /* Main loop */
t = from;
while t <= to

    if (learning)
        if (t > from + LPERIOD)
    		learning = 0;
        	T1 = T0;
            t = from;	% /* start over */
        else
            T1 = 2*T0;
        end
    end

	[temp,ebuf,lbuf,tt_2, aet] = slpsamp(t,data,BUFLN,ebuf,lbuf,tt_2, aet,SLPwindow);

    if (temp > T1)    % /* found a possible ABP pulse near t */
	    timer = 0;
            % /* used for counting the time after previous ABP pulse */
	    maxd = temp;
        mind = maxd;
        tmax=t;
        for (tt = t + 1: t + EyeClosing-1)
            [temp2 ,ebuf,lbuf,tt_2, aet] = slpsamp(tt,data,BUFLN,ebuf,lbuf,tt_2, aet,SLPwindow);
            if temp2 > maxd
                maxd=temp2;
                tmax=tt;
            end
        end
        if (maxd == temp)
            t=t+1;
            continue;
        end

        for tt = tmax :-1: (t - EyeClosing / 2 +1)
            [temp2 ,ebuf,lbuf,tt_2, aet] = slpsamp(tt,data,BUFLN,ebuf,lbuf,tt_2, aet,SLPwindow);
            if temp2< mind
                mind=temp2;
            end
        end
        if maxd>mind+10
            onset=(maxd-mind)/100+2;
            tpq=t-round(0.04*fs);
            maxmin_2_3_threshold=(maxd-mind)*2.0/3;
            for tt=tmax:-1:t-EyeClosing/2+1
                [temp2, ebuf,lbuf,tt_2, aet] = slpsamp(tt,data,BUFLN,ebuf,lbuf,tt_2, aet,SLPwindow);
                if temp2<maxmin_2_3_threshold
                    break;
                end
            end
            for tt=tt:-1:t - EyeClosing / 2 + round(0.024*fs)
                [temp2 ,ebuf,lbuf,tt_2, aet] = slpsamp(tt,data,BUFLN,ebuf,lbuf,tt_2, aet,SLPwindow);
                [temp3 ,ebuf,lbuf,tt_2, aet] = slpsamp(tt-round(0.024*fs),data,BUFLN,ebuf,lbuf,tt_2, aet,SLPwindow);
                if temp2-temp3<onset
                    tpq=tt-round(0.016*fs);
                    break;
                end
            end

            % find valley from the original data around 0.25s of tpq
            valley_v = round(tpq);
            for valley_i=round(max(2,tpq-round(0.20*fs))):round(min(tpq+round(0.05*fs),length(data)-1))

                % If vally is too low, it cannot serve as an index, so move to the next time.
                if valley_v <= 0
                    t = t + 1;
                    continue;
                end

                if data(valley_v)>data(valley_i) && data(valley_i)<=data(valley_i-1) && data(valley_i)<=data(valley_i+1)
                    valley_v=valley_i;
                end
            end

            if (~learning)

                % If we are looking for the first peak
                if beat_n == 1

                    % If the proposed peak index > 0
                    if round(valley_v) > 0
                        ppgOnsets(beat_n) = round(valley_v);
                        beat_n = beat_n + 1;
                    end
                else
                    % Check if rounded valley_v is greater than the prior beat index
                    if round(valley_v) > ppgOnsets(beat_n-1)
                        ppgOnsets(beat_n) = round(valley_v);
                        beat_n = beat_n + 1;
                    end
                end
            end

            % /* Adjust thresholds */
            Ta = Ta + (maxd - Ta)/10;
            T1 = Ta / 3;

            % /* Lock out further detections during the eye-closing period */
            t = tpq+EyeClosing;
        end
    else
        if (~learning)
	    % /* Once past the learning period, decrease threshold if no pulse
	    %   was detected recently. */
            timer = timer+1;
            if (timer > ExpectPeriod && Ta > Tm)
                Ta=Ta-1;
                T1 = Ta / 3;
            end
        end
    end

    t=t+1;

end

% Discard first beat because algorithm always finds first minimum value, so trace-back logic
% will find a fir

end

function [beat1,ebuf,lbuf,tt_2, aet] = slpsamp(t,data,BUFLN,ebuf,lbuf,tt_2, aet,SLPwindow)

    while (t > tt_2)
        prevVal=0;

        if (tt_2>0) && (tt_2-1>0) && (tt_2<length(data)) && (tt_2-1<length(data))
            val2=data(tt_2 - 1);
            val1=data(tt_2);
        else
            val2=prevVal;
            val1=val2;
        end
    	prevVal=val2;
        dy =  val1-val2;
        if (dy < 0)
            dy = 0;
        end
        tt_2=tt_2+1;
        M=round(mod(tt_2,(BUFLN-1))+1);
        et=dy;
        ebuf(M)=et;
%         M2=round(mod(tt_2-SLPwindow,(BUFLN-1))+1);
%         aet=aet+et-ebuf(M2);
        aet=0;
        for i=0:SLPwindow-1
            p=M-i;
            if p<=0
                p=p+BUFLN;
            end
            aet=aet+ebuf(p);
        end
        lbuf(M) = aet;

    end
    M3=round(mod(t,(BUFLN-1))+1);
    beat1=lbuf(M3);
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function hrfs = getcanonicalhrflibrary(duration,tr)

% function hrfs = getcanonicalhrflibrary(duration,tr)
%
% <duration> is the duration of the stimulus in seconds.
%   should be a multiple of 0.1 (if not, we round to the nearest 0.1).
%   0 is automatically treated as 0.1.
% <tr> is the TR in seconds.
%
% generate a library of 20 predicted HRFs to a stimulus of
% duration <duration>, with data sampled at a TR of <tr>.
%
% the resulting HRFs are returned as 20 x time. the first point is
% coincident with stimulus onset. each HRF is normalized such
% that the maximum value is one.
%
% example:
% hrfs = getcanonicalhrflibrary(4,1);
% figure; plot(0:size(hrfs,2)-1,hrfs,'o-');

% inputs
if duration == 0
  duration = 0.1;
end

% load the library
file0 = strrep(which('getcanonicalhrflibrary'),'getcanonicalhrflibrary.m','getcanonicalhrflibrary.tsv');
hrfs = load(file0)';  % 20 HRFs x 501 time points

% convolve to get the predicted response to the desired stimulus duration
trold = 0.1;
hrfs = conv2(hrfs,ones(1,max(1,round(duration/trold))));

% resample to desired TR
hrfs = interp1((0:size(hrfs,2)-1)*trold,hrfs',0:tr:(size(hrfs,2)-1)*trold,'pchip')';  % 20 HRFs x time

% make the peak equal to one
hrfs = hrfs ./ repmat(max(hrfs,[],2),[1 size(hrfs,2)]);

%%%%%%%%%%%%%%%%%%% FOR OUR RECORDS BELOW:

% % params taken from the Natural Scenes Dataset
% a1 = load('~/nsd/nsddata/templates/hrfparams.mat');
%
% % obtain canonical response to a 0.1-s stimulus
% hrfs = [];
% for p=1:size(a1.params,1)
%   hrfs(p,:) = spm_hrf(0.1,a1.params(p,:));
% end
%
% dlmwrite('getcanonicalhrflibrary.tsv',hrfs','delimiter','\t','precision',5);
% test=load('getcanonicalhrflibrary.tsv');
% figure; plot(hrfs');
% figure; plot(test);

author: cvnlab

repo: GLMsingle

save_path: github-repos/MATLAB/cvnlab-GLMsingle


text:
%compute sumspinerelangle

function [data,units]=compute_sumspinerelangle(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
spineangles=cell(1,numlarvae);
spinerelangle=cell(1,numlarvae);
sumspinerelangle=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    spineangles{1,i}=bsxfun(@atan2,trx(larva).yspine_mm(1:end-1,:)-trx(larva).yspine_mm(2:end,:),trx(larva).xspine_mm(1:end-1,:)-trx(larva).xspine_mm(2:end,:));
    % spinerelangle{1,i}=bsxfun(@atan2,sin(spineangles{1,i}(2:end,:))-sin(spineangles{1,i}(1:end-1,:)),cos(spineangles{1,i}(2:end,:))-cos(spineangles{1,i}(1:end-1,:)));
    % KB: this method of measuring the difference between angles didn't
    % make sense to me.
    spinerelangle{1,i} = modrange(-diff(spineangles{i},1,1),-pi,pi);
    sumspinerelangle{1,i}=sum(spinerelangle{1,i},1);
end

units=parseunits('rad');
data=sumspinerelangle;

author: kristinbranson

repo: JAABA


text:
function inside = quad_contains_point_2d ( q, p )

%*****************************************************************************80
%
%% QUAD_CONTAINS_POINT_2D finds if a point is inside a convex quadrilateral in 2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real Q(2,4), the vertices of the quadrilateral.
%
%    Input, real P(2,1), the point to be checked.
%
%    Output, logical INSIDE, is TRUE if the point is in the quadrilateral.
%

%
%  This will only handle convex quadrilaterals.
%
  inside = 0;

  angle_1 = angle_rad_2d ( q(1:2,1), q(1:2,2), q(1:2,3) );
  angle_2 = angle_rad_2d ( q(1:2,1), q(1:2,2), p(1:2,1) );

  if ( angle_1 < angle_2 )
    return
  end

  angle_1 = angle_rad_2d ( q(1:2,2), q(1:2,3), q(1:2,4) );
  angle_2 = angle_rad_2d ( q(1:2,2), q(1:2,3), p(1:2,1) );

  if ( angle_1 < angle_2 )
    return
  end

  angle_1 = angle_rad_2d ( q(1:2,3), q(1:2,4), q(1:2,1) );
  angle_2 = angle_rad_2d ( q(1:2,3), q(1:2,4), p(1:2,1) );

  if ( angle_1 < angle_2 )
    return
  end

  angle_1 = angle_rad_2d ( q(1:2,4), q(1:2,1), q(1:2,2) );
  angle_2 = angle_rad_2d ( q(1:2,4), q(1:2,1), p(1:2,1) );

  if ( angle_1 < angle_2 )
    return
  end

  inside = 1;

  return
end


text:
function varargout = quiver3( Z, F, varargin )
%QUIVER3  3-D quiver plot of a CHEBFUN2V at data mapped by a SEPARABLEAPPROX.
%
% QUIVER3(Z, F) plots velocity vectors at the equally spaced surface points
% specified by the SEPARABLEAPPROX Z. We use Z to map a uniform grid. F should be
% a CHEBFUN2V.
%
% QUIVER3(X, Y, Z, F) plots velocity vectors at (x,y,z), where X, Y, Z are
% SEPARABLEAPPROX objects which we use to to map a uniform grid. F should be a
% CHEBFUN2V.
%
% Alternative syntax for this command is:
% QUIVER3(X,Y,Z,[f;g;h]) or QUIVER3(X,Y,Z,f,g,h), where f, g, and h are
% SEPARABLEAPPROX objects.
%
% QUIVER(...,'numpts',N) plots arrows on a N x N uniform grid.
%
% This command is a wrapper to CHEBFUN2V/QUIVER3, and is required because
% SEPARABLEAPPROX methods take priority over CHEBFUN2V methods.
%
% See also CHEBFUN2V/QUIVER3.

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

numpts = 20;

if ( ~isa(Z, 'separableApprox') )
    error('CHEBFUN:SEPARABLEAPPROX:quiver3:inputs1', ...
        'First argument to this command should be SEPARABLEAPPROX.');
end

if ( isempty(varargin) )
    varargin = {};
end

% Number of points to plot
j = 1; argin = {};
while ( ~isempty(varargin) )
    if ( strcmpi(varargin{1}, 'numpts') )
        numpts = varargin{2};
        varargin(1:2) = [];
    else
        argin{j} = varargin{1};
        varargin(1) = [];
        j = j+1;
    end
end
varargin = argin;

if ( isa(F, 'chebfun2v') )              % quiver(Z,F,...)
    zz = new_data_locations(Z, numpts);
    h = quiver3(zz, F, varargin{:});
elseif ( nargin > 3 )
    if ( isa(varargin{2}, 'chebfun2v') ) % quiver(X,Y,Z,F,...)
        xx = new_data_locations(Z, numpts);
        yy = new_data_locations(F, numpts);
        zz = new_data_locations(varargin{1}, numpts);
        h = quiver3( xx, yy, zz, varargin{2}, varargin{3:end} );
    elseif ( isa(Z, 'separableApprox') && isa(F, 'separableApprox') &&...
            isa(varargin{1}, 'separableApprox') &&...
            isa(varargin{2}, 'separableApprox') &&...
            ( nargin ==4 || ~isa(varargin{3}, 'separableApprox')) )
        FF = vertcat( vertcat(F, varargin{1}), varargin{2}) ;
        h = quiver3( Z, FF );     % call quiver3(Z,F,...)
    elseif ( ( nargin > 5 ) && isa(Z,'separableApprox') && isa(F,'separableApprox') &&...
            isa(varargin{1}, 'separableApprox') &&...
            isa(varargin{2}, 'separableApprox') &&...
            isa(varargin{3}, 'separableApprox') &&...
            isa(varargin{4}, 'separableApprox') &&...
            ( nargin == 6 || ~isa(varargin{5}, 'separableApprox')) )
        FF = vertcat( vertcat(varargin{2}, varargin{3}), varargin{4} );
        h = quiver3(Z, F, varargin{1}, FF, varargin{5:end} );  % call quiver3(X,Y,Z,F,...)
    else
        error('CHEBFUN:SEPARABLEAPPROX:quiver3:inputs2', ...
            'Unrecognised input arguments.');
    end
else
    error('CHEBFUN:SEPARABLEAPPROX:quiver3:inputs3', 'Unrecognised input arguments.');
end

if ( nargout > 0 )
    varargout = { h };
end

end

function newloc = new_data_locations(f1, numpts)
% Generate new arrow location if first two inputs are SEPARABLEAPPROX objects.

dom = f1.domain;

% mesh 'em up for the quiver arrows.
x = linspace(dom(1), dom(2), numpts);
y = linspace(dom(3), dom(4), numpts);

[xx, yy] = meshgrid(x, y);
newloc = feval(f1, xx, yy);      % use SEPARABLEAPPROX to generate data locations.

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [W]=downbeam_lambda(A,H,K,M,grt,beta,lambda)
    W=zeros(M,K);
    A=A+lambda.*eye(M);
    A=A^(-1);
    for k0=1:K
        W(:,k0)=sqrt(grt(k0))*beta(k0)*A*H(k0,:)';
    end
end

author: guohuayan


text:
function [ y, next ] = p37_equil ( neqn, next )

%*****************************************************************************80
%
%% P37_EQUIL returns equilibrium solutions of problem p37.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, integer NEXT, the index of the previous
%    equilibrium, which should be 0 on first call.
%
%    Output, real Y(NEQN), the "next" equilibrium solution, if any.
%
%    Output, integer NEXT, the index of the current equilibrium,
%    or 0 if there are no more.
%
  if ( next == 0 )
    next = 1;
    y(1:neqn,1) = [ 0.0; 0.0 ];
  else
    next = 0;
    y = [];
  end

  return
end


text:
function [LR] = imresize_down(im, scale, type, sigma)

if nargin ==3 && strcmp(type,'Gaussian')
    sigma = 1.6;
end

if strcmp(type,'Gaussian') && fix(scale) == scale
    if mod(scale,2)==1
        kernelsize = ceil(sigma*3)*2+1;
        if scale==3 && sigma == 1.6
            kernelsize = 7;
        end
        kernel  = fspecial('gaussian',kernelsize,sigma);
        blur_HR = imfilter(im,kernel,'replicate');

        if isa(blur_HR, 'gpuArray')
            LR = blur_HR(scale-1:scale:end-1,scale-1:scale:end-1,:);
        else
            LR      = imresize(blur_HR, 1/scale, 'nearest');
        end

        % LR      = im2uint8(LR);
    elseif mod(scale,2)==0
        kernelsize = ceil(sigma*3)*2+2;
        kernel     = fspecial('gaussian',kernelsize,sigma);
        blur_HR    = imfilter(im, kernel,'replicate');
        LR= blur_HR(scale/2:scale:end-scale/2,scale/2:scale:end-scale/2,:);
        % LR         = im2uint8(LR);
    end
else
    LR = imresize(im, 1/scale, type);
end

author: cszn

repo: IRCNN

save_path: github-repos/MATLAB/cszn-IRCNN


text:
function p = presolveOneMagicRound(p);

p = presolve_bounds_from_domains(p);
p = presolve_bounds_from_modelbounds(p);
p = propagate_bounds_from_convex_quadratic_ball(p);
p = propagate_bounds_from_evaluations(p);
p = update_monomial_bounds(p);
p = propagate_bounds_from_equalities(p);
p = presolve_bounds_from_inequalities(p);
p = propagate_bounds_from_evaluations(p);
p = update_monomial_bounds(p);
p = propagate_bounds_from_evaluations(p);
p = update_monomial_bounds(p);
p = propagate_bounds_from_arbitrary_quadratics(p);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


% apply the ECO tracking for each target (modified based on the original ECO code)
function tracker = ECO_tracking(frame_id, im, bboxes_det, tracker, opt)

tracker.eco.n_frame = tracker.eco.n_frame + 1;
frame = tracker.eco.n_frame;

% variables from ECO_initialize.m
params = tracker.eco.params;
max_train_samples = tracker.eco.max_train_samples;
features = tracker.eco.features;
global_fparams = tracker.eco.global_fparams;
pos = tracker.eco.pos;
target_sz = tracker.eco.target_sz;
currentScaleFactor = tracker.eco.currentScaleFactor;
base_target_sz = tracker.eco.base_target_sz;
img_support_sz = tracker.eco.img_support_sz;
feature_dim = tracker.eco.feature_dim;
num_feature_blocks = tracker.eco.num_feature_blocks;
feature_reg = tracker.eco.feature_reg;
feature_extract_info = tracker.eco.feature_extract_info;
compressed_dim = tracker.eco.compressed_dim;
compressed_dim_cell = tracker.eco.compressed_dim_cell;
filter_sz = tracker.eco.filter_sz;
filter_sz_cell = tracker.eco.filter_sz_cell;
output_sz = tracker.eco.output_sz;
pad_sz = tracker.eco.pad_sz;
kx = tracker.eco.kx;
ky = tracker.eco.ky;
yf = tracker.eco.yf;
cos_window = tracker.eco.cos_window;
interp1_fs = tracker.eco.interp1_fs;
interp2_fs = tracker.eco.interp2_fs;
reg_filter = tracker.eco.reg_filter;
reg_energy = tracker.eco.reg_energy;
nScales = tracker.eco.nScales;
scaleFactors = tracker.eco.scaleFactors;
scale_filter = tracker.eco.scale_filter;
min_scale_factor = tracker.eco.min_scale_factor;
max_scale_factor = tracker.eco.max_scale_factor;
init_CG_opts = tracker.eco.init_CG_opts;
CG_opts = tracker.eco.CG_opts;
rect_position = tracker.eco.rect_position;
prior_weights = tracker.eco.prior_weights;
sample_weights = tracker.eco.sample_weights;
samplesf = tracker.eco.samplesf;
score_matrix = tracker.eco.score_matrix;
latest_ind = tracker.eco.latest_ind;
frames_since_last_train = tracker.eco.frames_since_last_train;
num_training_samples = tracker.eco.num_training_samples;
minimum_sample_weight = tracker.eco.minimum_sample_weight;
res_norms = tracker.eco.res_norms;
is_color_image = tracker.eco.is_color_image;

% variables which are useful when frame == 1
sample_pos = tracker.eco.sample_pos;
sample_scale = tracker.eco.sample_scale;
xl = tracker.eco.xl;
xlf = tracker.eco.xlf;
projection_matrix = tracker.eco.projection_matrix;
shift_samp = tracker.eco.shift_samp;
xlf_proj = tracker.eco.xlf_proj;
hf = tracker.eco.hf;
lf_ind = tracker.eco.lf_ind;
proj_energy = tracker.eco.proj_energy;
sample_energy = tracker.eco.sample_energy;
rhs_samplef = tracker.eco.rhs_samplef;
diag_M = tracker.eco.diag_M;
p = tracker.eco.p;
rho = tracker.eco.rho;
r_old = tracker.eco.r_old;
init_samplef = tracker.eco.init_samplef;
init_samplef_H = tracker.eco.init_samplef_H;
projection_matrix_init = tracker.eco.projection_matrix_init;
init_samplef_proj = tracker.eco.init_samplef_proj;
init_hf = tracker.eco.init_hf;
fyf = tracker.eco.fyf;
res_norms_temp = tracker.eco.res_norms_temp;

% other variables
hf_full = tracker.eco.hf_full;

% narrow the bbox to avoid tracking drift
target_sz(2) = target_sz(2) * 0.5;

% load image
if size(im,3) > 1 && is_color_image == false
    im = im(:,:,1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Target localization step
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Do not estimate translation and scaling on the first frame, since we
% just want to initialize the tracker there

if frame > 1
    old_pos = inf(size(pos));
    iter = 1;

    %translation search
    while iter <= params.refinement_iterations && any(old_pos ~= pos)
        % Extract features at multiple resolutions
        sample_pos = round(pos);
        det_sample_pos = sample_pos;
        sample_scale = currentScaleFactor*scaleFactors;
        xt = extract_features(im, sample_pos, sample_scale, features, global_fparams, feature_extract_info);

        % Project sample
        xt_proj = project_sample(xt, projection_matrix);

        % Do windowing of features
        xt_proj = cellfun(@(feat_map, cos_window) bsxfun(@times, feat_map, cos_window), xt_proj, cos_window, 'uniformoutput', false);

        % Compute the fourier series
        xtf_proj = cellfun(@cfft2, xt_proj, 'uniformoutput', false);

        % Interpolate features to the continuous domain
        xtf_proj = interpolate_dft(xtf_proj, interp1_fs, interp2_fs);

        % Compute convolution for each feature block in the Fourier domain
        scores_fs_feat = cellfun(@(hf, xf, pad_sz) padarray(sum(bsxfun(@times, hf, xf), 3), pad_sz), hf_full, xtf_proj, pad_sz, 'uniformoutput', false);

        % Also sum over all feature blocks.
        % Gives the fourier coefficients of the convolution response.
        scores_fs = permute(sum(cell2mat(scores_fs_feat), 3), [1 2 4 3]);

        % Optimize the continuous score function with Newton's method.
        [trans_row, trans_col, scale_ind] = optimize_scores(scores_fs, params.newton_iterations);

        % Compute the translation vector in pixel-coordinates and round
        % to the closest integer pixel.
        translation_vec = [trans_row, trans_col] .* (img_support_sz./output_sz) * currentScaleFactor * scaleFactors(scale_ind);
        scale_change_factor = scaleFactors(scale_ind);

        % update position
        old_pos = pos;
        pos = sample_pos + translation_vec;

        if params.clamp_position
            pos = max([1 1], min([size(im,1) size(im,2)], pos));
        end

        % Do scale tracking with the scale filter
        if nScales > 0 && params.use_scale_filter
            scale_change_factor = scale_filter_track(im, pos, base_target_sz, currentScaleFactor, scale_filter, params);
        end

        % Update the scale
        currentScaleFactor = currentScaleFactor * scale_change_factor;

        % Adjust to make sure we are not to large or to small
        if currentScaleFactor < min_scale_factor
            currentScaleFactor = min_scale_factor;
        elseif currentScaleFactor > max_scale_factor
            currentScaleFactor = max_scale_factor;
        end

        iter = iter + 1;
    end
end

pre_pos = tracker.eco.pos;
pre_target_sz = tracker.eco.target_sz;
pos_diff = double(sqrt((pos(1)-pre_pos(1)).^2 + (pos(2)-pre_pos(2)).^2));
if pos_diff > 1.2 * double(pre_target_sz(2))
    pos = pre_pos;
    is_drift = 1;
else
    is_drift = 0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Model update step
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extract sample and init projection matrix
if frame == 1
    % Extract image region for training sample
    sample_pos = round(pos);
    sample_scale = currentScaleFactor;
    xl = extract_features(im, sample_pos, currentScaleFactor, features, global_fparams, feature_extract_info);

    % Do windowing of features
    xlw = cellfun(@(feat_map, cos_window) bsxfun(@times, feat_map, cos_window), xl, cos_window, 'uniformoutput', false);

    % Compute the fourier series
    xlf = cellfun(@cfft2, xlw, 'uniformoutput', false);

    % Interpolate features to the continuous domain
    xlf = interpolate_dft(xlf, interp1_fs, interp2_fs);

    % New sample to be added
    xlf = compact_fourier_coeff(xlf);

    % Initialize projection matrix
    xl1 = cellfun(@(x) reshape(x, [], size(x,3)), xl, 'uniformoutput', false);
    xl1 = cellfun(@(x) bsxfun(@minus, x, mean(x, 1)), xl1, 'uniformoutput', false);

    if strcmpi(params.proj_init_method, 'pca')
        [projection_matrix, ~, ~] = cellfun(@(x) svd(x' * x), xl1, 'uniformoutput', false);
        projection_matrix = cellfun(@(P, dim) single(P(:,1:dim)), projection_matrix, compressed_dim_cell, 'uniformoutput', false);
    elseif strcmpi(params.proj_init_method, 'rand_uni')
        projection_matrix = cellfun(@(x, dim) single(randn(size(x,2), dim)), xl1, compressed_dim_cell, 'uniformoutput', false);
        projection_matrix = cellfun(@(P) bsxfun(@rdivide, P, sqrt(sum(P.^2,1))), projection_matrix, 'uniformoutput', false);
    elseif strcmpi(params.proj_init_method, 'none')
        projection_matrix = [];
    else
        error('Unknown initialization method for the projection matrix: %s', params.proj_init_method);
    end
    clear xl1 xlw

    % Shift sample
    shift_samp = 2*pi * (pos - sample_pos) ./ (sample_scale * img_support_sz);
    xlf = shift_sample(xlf, shift_samp, kx, ky);

    % Project sample
    xlf_proj = project_sample(xlf, projection_matrix);
elseif params.learning_rate > 0
    if ~params.use_detection_sample
        % Extract image region for training sample
        sample_pos = round(pos);
        sample_scale = currentScaleFactor;
        xl = extract_features(im, sample_pos, currentScaleFactor, features, global_fparams, feature_extract_info);

        % Project sample
        xl_proj = project_sample(xl, projection_matrix);

        % Do windowing of features
        xl_proj = cellfun(@(feat_map, cos_window) bsxfun(@times, feat_map, cos_window), xl_proj, cos_window, 'uniformoutput', false);

        % Compute the fourier series
        xlf1_proj = cellfun(@cfft2, xl_proj, 'uniformoutput', false);

        % Interpolate features to the continuous domain
        xlf1_proj = interpolate_dft(xlf1_proj, interp1_fs, interp2_fs);

        % New sample to be added
        xlf_proj = compact_fourier_coeff(xlf1_proj);
    else
        % Use the sample that was used for detection
        sample_scale = sample_scale(scale_ind);
        xlf_proj = cellfun(@(xf) xf(:,1:(size(xf,2)+1)/2,:,scale_ind), xtf_proj, 'uniformoutput', false);
    end

    % Shift the sample so that the target is centered
    shift_samp = 2*pi * (pos - sample_pos) ./ (sample_scale * img_support_sz);
    xlf_proj = shift_sample(xlf_proj, shift_samp, kx, ky);
end

xlf_proj_perm = cellfun(@(xf) permute(xf, [4 3 1 2]), xlf_proj, 'uniformoutput', false);

if params.use_sample_merge
    % Find the distances with existing samples
    dist_vector = find_cluster_distances(samplesf, xlf_proj_perm, num_feature_blocks, num_training_samples, max_train_samples, params);

    [merged_sample, new_cluster, merged_cluster_id, new_cluster_id, score_matrix, prior_weights,num_training_samples] = ...
        merge_clusters(samplesf, xlf_proj_perm, dist_vector, score_matrix, prior_weights,...
                       num_training_samples,num_feature_blocks,max_train_samples,minimum_sample_weight,params);
else
    % Do the traditional adding of a training sample and weight update
    % of C-COT
    [prior_weights, replace_ind] = update_prior_weights(prior_weights, sample_weights, latest_ind, frame, params);
    latest_ind = replace_ind;

    merged_cluster_id = 0;
    new_cluster = xlf_proj_perm;
    new_cluster_id = replace_ind;
end

if frame > 1 && params.learning_rate > 0 || frame == 1 && ~params.update_projection_matrix
    % Insert the new training sample
    for k = 1:num_feature_blocks
        if merged_cluster_id > 0
            samplesf{k}(merged_cluster_id,:,:,:) = merged_sample{k};
        end

        if new_cluster_id > 0
            samplesf{k}(new_cluster_id,:,:,:) = new_cluster{k};
        end
    end
end

sample_weights = prior_weights;

train_tracker = (frame < params.skip_after_frame) || (frames_since_last_train >= params.train_gap);

if train_tracker && is_drift == 0
    % Used for preconditioning
    new_sample_energy = cellfun(@(xlf) abs(xlf .* conj(xlf)), xlf_proj, 'uniformoutput', false);

    if frame == 1
        if params.update_projection_matrix
            hf = cell(2,1,num_feature_blocks);
            lf_ind = cellfun(@(sz) sz(1) * (sz(2)-1)/2 + 1, filter_sz_cell, 'uniformoutput', false);
            proj_energy = cellfun(@(P, yf) 2*sum(abs(yf(:)).^2) / sum(feature_dim) * ones(size(P), 'single'), projection_matrix, yf, 'uniformoutput', false);
        else
            hf = cell(1,1,num_feature_blocks);
        end
        % Initialize the filter
        for k = 1:num_feature_blocks
            hf{1,1,k} = complex(zeros([filter_sz(k,1) (filter_sz(k,2)+1)/2 compressed_dim(k)], 'single'));
        end

        % Initialize Conjugate Gradient parameters
        CG_opts.maxit = params.init_CG_iter; % Number of initial iterations if projection matrix is not updated
        init_CG_opts.maxit = ceil(params.init_CG_iter / params.init_GN_iter);
        sample_energy = new_sample_energy;
        rhs_samplef = cell(size(hf));
        diag_M = cell(size(hf));
        p = []; rho = []; r_old = [];
    else
        CG_opts.maxit = params.CG_iter;

        if params.CG_forgetting_rate == inf || params.learning_rate >= 1
            % CG will be reset
            p = []; rho = []; r_old = [];
        else
            rho = rho / (1-params.learning_rate)^params.CG_forgetting_rate;
        end
        % Update the approximate average sample energy using the learning
        % rate. This is only used to construct the preconditioner.
        sample_energy = cellfun(@(se, nse) (1 - params.learning_rate) * se + params.learning_rate * nse, sample_energy, new_sample_energy, 'uniformoutput', false);
    end

    % Do training
    if frame == 1 && params.update_projection_matrix
        % Initial Gauss-Newton optimization of the filter and
        % projection matrix.

        % Construct stuff for the proj matrix part
        init_samplef = cellfun(@(x) permute(x, [4 3 1 2]), xlf, 'uniformoutput', false);
        init_samplef_H = cellfun(@(X) conj(reshape(X, size(X,2), [])), init_samplef, 'uniformoutput', false);

        % Construct preconditioner
        diag_M(1,1,:) = cellfun(@(m, reg_energy) (1-params.precond_reg_param) * bsxfun(@plus, params.precond_data_param * m, (1-params.precond_data_param) * mean(m,3)) + params.precond_reg_param*reg_energy, sample_energy, reg_energy, 'uniformoutput',false);
        diag_M(2,1,:) = cellfun(@(m) params.precond_proj_param * (m + params.projection_reg), proj_energy, 'uniformoutput',false);

        projection_matrix_init = projection_matrix;

        for iter = 1:params.init_GN_iter
            % Project sample with new matrix
            init_samplef_proj = cellfun(@(x,P) mtimesx(x, P, 'speed'), init_samplef, projection_matrix, 'uniformoutput', false);
            init_hf = cellfun(@(x) permute(x, [3 4 1 2]), hf(1,1,:), 'uniformoutput', false);

            % Construct the right hand side vector for the filter part
            rhs_samplef(1,1,:) = cellfun(@(xf, yf) bsxfun(@times, conj(permute(xf, [3 4 2 1])), yf), init_samplef_proj, yf, 'uniformoutput', false);

            % Construct the right hand side vector for the projection matrix part
            fyf = cellfun(@(f, yf) reshape(bsxfun(@times, conj(f), yf), [], size(f,3)), hf(1,1,:), yf, 'uniformoutput', false);
            rhs_samplef(2,1,:) = cellfun(@(P, XH, fyf, fi) (2*real(XH * fyf - XH(:,fi:end) * fyf(fi:end,:)) - params.projection_reg * P), ...
                projection_matrix, init_samplef_H, fyf, lf_ind, 'uniformoutput', false);

            % Initialize the projection matrix increment to zero
            hf(2,1,:) = cellfun(@(P) zeros(size(P), 'single'), projection_matrix, 'uniformoutput', false);

            % do conjugate gradient
            [hf, ~, ~, ~, res_norms_temp] = pcg_ccot(...
                @(x) lhs_operation_joint(x, init_samplef_proj, reg_filter, feature_reg, init_samplef, init_samplef_H, init_hf, params.projection_reg),...
                rhs_samplef, init_CG_opts, ...
                @(x) diag_precond(x, diag_M), ...
                [], hf);

            % Make the filter symmetric (avoid roundoff errors)
            hf(1,1,:) = symmetrize_filter(hf(1,1,:));

            % Add to the projection matrix
            projection_matrix = cellfun(@plus, projection_matrix, hf(2,1,:), 'uniformoutput', false);

            res_norms = [res_norms; res_norms_temp];
        end

        % Extract filter
        hf = hf(1,1,:);

        % Re-project and insert training sample
        xlf_proj = project_sample(xlf, projection_matrix);
        for k = 1:num_feature_blocks
            samplesf{k}(1,:,:,:) = permute(xlf_proj{k}, [4 3 1 2]);
        end

        if debug
            norm_proj_mat_init = sqrt(sum(cellfun(@(P) norm(P(:))^2, projection_matrix_init)));
            norm_proj_mat = sqrt(sum(cellfun(@(P) norm(P(:))^2, projection_matrix)));
            norm_proj_mat_change = sqrt(sum(cellfun(@(P,P2) norm(P(:) - P2(:))^2, projection_matrix_init, projection_matrix)));
            fprintf('Norm init: %f, Norm final: %f, Matrix change: %f\n', norm_proj_mat_init, norm_proj_mat, norm_proj_mat_change / norm_proj_mat_init);
        end
    else
        % Construct the right hand side vector
        rhs_samplef = cellfun(@(xf) permute(mtimesx(sample_weights, 'T', xf, 'speed'), [3 4 2 1]), samplesf, 'uniformoutput', false);
        rhs_samplef = cellfun(@(xf, yf) bsxfun(@times, conj(xf), yf), rhs_samplef, yf, 'uniformoutput', false);

        % Construct preconditioner
        diag_M = cellfun(@(m, reg_energy) (1-params.precond_reg_param) * bsxfun(@plus, params.precond_data_param * m, (1-params.precond_data_param) * mean(m,3)) + params.precond_reg_param*reg_energy, sample_energy, reg_energy, 'uniformoutput',false);

        % do conjugate gradient
        [hf, ~, ~, ~, res_norms, p, rho, r_old] = pcg_ccot(...
            @(x) lhs_operation(x, samplesf, reg_filter, sample_weights, feature_reg),...
            rhs_samplef, CG_opts, ...
            @(x) diag_precond(x, diag_M), ...
            [], hf, p, rho, r_old);
    end

    % Reconstruct the full Fourier series
    hf_full = full_fourier_coeff(hf);

    frames_since_last_train = 0;
else
    frames_since_last_train = frames_since_last_train+1;
end

% Update the scale filter
if nScales > 0 && params.use_scale_filter
    scale_filter = scale_filter_update(im, pos, base_target_sz, currentScaleFactor, scale_filter, params);
end

% Update the target size (only used for computing output box)
target_sz = base_target_sz * currentScaleFactor;

% restore the size of bbox (narrowed before to avoid tracking drift)
target_sz(2) = target_sz(2) * 2.0;

%save position and calculate FPS
rect_position(frame,:) = round([pos([2,1]) - (target_sz([2,1]) - 1)/2, target_sz([2,1])]);

sampled_scores_display = fftshift(sample_fs(scores_fs(:,:,scale_ind), 10*output_sz));

pre_pos = tracker.eco.pos;
pre_target_sz = tracker.eco.target_sz;
bbox.x = double(pre_pos(2) - (pre_target_sz(2) - 1)/2);
bbox.y = double(pre_pos(1) - (pre_target_sz(1) - 1)/2);
bbox.w = pre_target_sz(2);
bbox.h = pre_target_sz(1);

if isempty(bboxes_det.fr) == 0
    o = calc_overlap(bbox, 1, bboxes_det, 1:numel(bboxes_det.fr));
    [overlap_box, index] = max(o);
    pre_area = bbox.w * bbox.h;
    cur_area = bboxes_det.w(index) * bboxes_det.h(index);
    if overlap_box > opt.det_overlap_thre && max(bbox.h, bboxes_det.h(index))/min(bbox.h, bboxes_det.h(index)) <= 1.3
        pos_det = single([bboxes_det.y(index) + (bboxes_det.h(index) - 1)/2, bboxes_det.x(index) + (bboxes_det.w(index) - 1)/2]);
        target_sz_det = [bboxes_det.h(index), bboxes_det.w(index)];
        pos = single(overlap_box * pos_det + (1 - overlap_box) * pos);
        target_sz = overlap_box * target_sz_det + (1 - overlap_box) * target_sz;
    end
end

% variables from ECO_initialize.m
tracker.eco.params = params;
tracker.eco.max_train_samples = max_train_samples;
tracker.eco.features = features;
tracker.eco.global_fparams = global_fparams;
tracker.eco.pos = pos;
tracker.eco.target_sz = target_sz;
tracker.eco.currentScaleFactor = currentScaleFactor;
tracker.eco.base_target_sz = base_target_sz;
tracker.eco.img_support_sz = img_support_sz;
tracker.eco.feature_dim = feature_dim;
tracker.eco.num_feature_blocks = num_feature_blocks;
tracker.eco.feature_reg = feature_reg;
tracker.eco.feature_extract_info = feature_extract_info;
tracker.eco.compressed_dim = compressed_dim;
tracker.eco.compressed_dim_cell = compressed_dim_cell;
tracker.eco.filter_sz = filter_sz;
tracker.eco.filter_sz_cell = filter_sz_cell;
tracker.eco.output_sz = output_sz;
tracker.eco.pad_sz = pad_sz;
tracker.eco.kx = kx;
tracker.eco.ky = ky;
tracker.eco.yf = yf;
tracker.eco.cos_window = cos_window;
tracker.eco.interp1_fs = interp1_fs;
tracker.eco.interp2_fs = interp2_fs;
tracker.eco.reg_filter = reg_filter;
tracker.eco.reg_energy = reg_energy;
tracker.eco.nScales = nScales;
tracker.eco.scaleFactors = scaleFactors;
tracker.eco.scale_filter = scale_filter;
tracker.eco.min_scale_factor = min_scale_factor;
tracker.eco.max_scale_factor = max_scale_factor;
tracker.eco.init_CG_opts = tracker.eco.init_CG_opts;
tracker.eco.CG_opts = CG_opts;
tracker.eco.rect_position = rect_position;
tracker.eco.prior_weights = prior_weights;
tracker.eco.sample_weights = sample_weights;
tracker.eco.samplesf = samplesf;
tracker.eco.score_matrix = score_matrix;
tracker.eco.latest_ind = latest_ind;
tracker.eco.frames_since_last_train = frames_since_last_train;
tracker.eco.num_training_samples = num_training_samples;
tracker.eco.minimum_sample_weight = minimum_sample_weight;
tracker.eco.res_norms = res_norms;
tracker.eco.is_color_image = is_color_image;

% variables which are useful when frame == 1
tracker.eco.sample_pos = sample_pos;
tracker.eco.sample_scale = sample_scale;
tracker.eco.xl = xl;
tracker.eco.xlf = xlf;
tracker.eco.projection_matrix = projection_matrix;
tracker.eco.shift_samp = shift_samp;
tracker.eco.xlf_proj = xlf_proj;
tracker.eco.hf = hf;
tracker.eco.lf_ind = lf_ind;
tracker.eco.proj_energy = proj_energy;
tracker.eco.sample_energy = sample_energy;
tracker.eco.rhs_samplef = rhs_samplef;
tracker.eco.diag_M = diag_M;
tracker.eco.p = p;
tracker.eco.rho = rho;
tracker.eco.r_old = r_old;
tracker.eco.init_samplef = init_samplef;
tracker.eco.init_samplef_H = init_samplef_H;
tracker.eco.projection_matrix_init = projection_matrix_init;
tracker.eco.init_samplef_proj = init_samplef_proj;
tracker.eco.init_hf = init_hf;
tracker.eco.fyf = fyf;
tracker.eco.res_norms_temp = res_norms_temp;

% other variables
tracker.eco.hf_full = hf_full;

tracker.eco.bb = double([pos([2,1]) - (target_sz([2,1]) - 1)/2, pos([2,1]) + (target_sz([2,1]) - 1)/2]);
tracker.eco.score = max(max(fftshift(sample_fs(scores_fs(:,:,scale_ind), 10*output_sz))));

if tracker.eco.score > opt.tracking_score_thre
    tracker.eco.is_confident = 1;
else
    tracker.eco.is_confident = 0;
end


text:
function about

%Tools for maintaining and programming in the GISMO environment.
% The admin package provides a series of tools to assist with programming
% in GISMO.

help admin;

author: geoscience-community-codes

repo: GISMO


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

%%
clear; close all; clc;

%% Syntax
% |A=tri_area(F,V);|

%% 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 cosmo_plot_slices(data, dim, slice_step, slice_start, slice_stop)
% Plots a set of slices from a dataset, nifti image, or 3D data array
%
% cosmo_plot_slices(data[, dim][, slice_step][, slice_start][, slice_stop])
%
% Inputs:
%  data         an fmri dataset (e.g., from cosmo_fmri_dataset), or a 3D
%               array with data. data should contain samples from a single
%               volume (sample) only.
%  dim          dimension according to which slices are plotted
%               (default: 3). Values between 1 and 3 allows using a
%               saggital, axial or coronal view; but the mapping between
%               the 3 numbers and 3 views depends on the particular
%               orientation of the data.
%  slice_step   step between slices (default: 1). If negative then
%               -slice_step indicates the number of slices
%  slice_start  the index of the first slice to plot (default: 1).
%  slice_stop   the index of the last slice to plot (default: the number of
%               slices in the dim-th dimension).
%
% Examples:
%    % plot an fMRI dataset struct with default options
%    cosmo_plot_slices(ds)
%
%    % plot an fMRI dataset struct along the second spatial dimension
%    cosmo_plot_slices(ds, 2)
%
%    % plot a random gaussian 3D array along the first dimension
%    cosmo_plot_slices(randn([40,50,20]),1)
%
%    % plot an fMRI dataset struct along the default spatial dimension
%    % every 5-th slice
%    cosmo_plot_slices(ds, [], 5)
%
%    % plot an fMRI dataset struct along the third spatial dimension
%    % with about 12 slices
%    cosmo_plot_slices(ds, 3, -12)
%
%    % plot an fMRI dataset struct along the third spatial dimension
%    % with about 12 slices, starting at slice 10 and stopping at slice 25
%    cosmo_plot_slices(ds, 3, -12, 10, 25)
%
% Notes:
%  - Using this function only really makes sense for fMRI-like data.
%  - This function does not provide a consistent orientation for slices,
%    as this depends on the voxel-to-world transformation matrix, which is
%    completely ignored in this function. Thus left-right and top-down
%    swaps can occur. Different datasets may provide different views, for
%    example dim=1 may give a saggital view if the dataset comes from one
%    program and an axial view if it comes from another program.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    if nargin<2 || isempty(dim), dim=3; end
    if nargin<3 || isempty(slice_step), slice_step=-20; end
    if nargin<4 || isempty(slice_start), slice_start=1; end
    if nargin<5 || isempty(slice_stop), slice_stop=[]; end % set later

    data3D=get_data3D(data);

    % get min and max values across the entire volume
    data_lin=data3D(:);
    mn=min(data_lin);
    mx=max(data_lin);

    % shift it so that we can walk over the first dimension
    data_sh=shiftdim(data3D, dim-1);

    if isempty(slice_stop)
        slice_stop=size(data_sh,1);
    end

    if slice_step<0
        nslices=-slice_step;
        slice_step=ceil((slice_stop-slice_start+1)/nslices);
    end

    % determine which slices to show
    slice_idxs=slice_start:slice_step:slice_stop;
    nslices=numel(slice_idxs);

    plot_ratio=.8; % ratio between number of rows and colums
    nrows=ceil(sqrt(nslices)*plot_ratio);
    ncols=ceil(nslices/nrows);

    % use header depending on dim
    header_labels={'i','j','k'};

    % order of slices and whether the slice should be transposed
    xorder=[-1 -1 1];
    yorder=[-1 1 -1];
    do_transpose=[true false true];

    for k=1:nslices
        slice_idx=slice_idxs(k);
        slice=squeeze(data_sh(slice_idx,:,:));

        if xorder(dim)<0
            slice=slice(end:-1:1,:);
        end
        if yorder(dim)<0
            slice=slice(:,end:-1:1);
        end
        if do_transpose(dim)
            slice=slice';
        end

        subplot(nrows, ncols, k);
        imagesc(slice, [mn, mx]);
        title(sprintf('%s = %d', header_labels{dim}, slice_idx));
    end

function data3D=get_data3D(data)
    if isstruct(data)
        cosmo_check_dataset(data,'fmri');
        data4D=cosmo_unflatten(data);
        sz=size(data4D);
        if sz(1)>1
            error(['expected single volume data, but found %d '...
                    'volumes. To select a single volume, use '...
                    'cosmo_slice. For example, to show the %d-th '...
                    'volume from a dataset struct ds, use:\n\n   '...
                    'cosmo_plot_slices(cosmo_slice(ds,%d))\n'],...
                    sz(1),sz(1),sz(1));
        end

        if numel(sz)>4
            error('data must be 4D');
        end

        data3D=reshape(data4D, sz(2:end));
    elseif isnumeric(data) || islogical(data)
        ndim=numel(size(data));
        if ndim>3
            error('expected 3D image - did you select a single volume?');
        end
        data3D=data;
    else
        error('illegal input: expected dataset or 3D array');
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
%% SQUAREPOISSONWG Poisson equation in a square domain solved by Weak Galerkin method.
%
%   squarePoissoWGn computes the lowest order Weak Galerkin approximations
%   of the Poisson equation in the unit square on a sequence of meshes
%   obtained by uniform refinement.
%
% See also PoissonWG,
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

close all; clear variables;

%% Parameters
maxIt = 4;
N = zeros(maxIt,1);
h = zeros(maxIt,1);
erruIuh = zeros(maxIt,1);

%% Generate an initial mesh
[node,elem] = squaremesh([0 1 0 1], 0.25);
%bdFlag = setboundary(node,elem,'Dirichlet','~(x==0)','Neumann','x==0');
bdFlag = setboundary(node,elem,'Dirichlet','all','Neumann','y==1');
%bdFlag = setboundary(node,elem,'Dirichlet');
for k = 1:3
    [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
%     [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
end

%% Get the data of the pde
pde = sincosdata;
% pde = mixBCdata;
% option.solver = 'mg';

%% Finite Element Method
for k = 1:maxIt
    % solve the equation
    [soln,eqn,info] = PoissonWG(node,elem,bdFlag,pde);
    N(k) = size(elem,1);
    h(k) = 1./(sqrt(size(node,1))-1);
    % compute error
    uI = zeros(N(k)+size(eqn.edge,1),1);
    uI(1:N(k)) = pde.exactu((node(elem(:,1),:)+node(elem(:,2),:)+node(elem(:,3),:))/3);
    uI(N(k)+1:end) = pde.exactu((node(eqn.edge(:,1),:)+node(eqn.edge(:,2),:))/2);
    erruIuh(k) = sqrt((soln.u-uI)'*eqn.A*(soln.u-uI));
    % refine mesh
   [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
%     [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
end

%% Plot convergence rates
figure(2);
showrateh(h,erruIuh);

repo: ifem


text:
function [fitresult, gof, x, y] = createFit(num_feats, mean_r,varargin)
%CREATEFIT(NUM_FEATS,MEAN_R)
%  Create a fit.
%
%  Data for 'untitled fit 1' fit:
%      X Input : num_feats
%      Y Output: mean_r
%  Output:
%      fitresult : a fit object representing the fit.
%      gof : structure with goodness-of fit info.
%
%  See also FIT, CFIT, SFIT.

%  Auto-generated by MATLAB on 09-Jan-2018 14:30:55

%% Fit: 'untitled fit 1'.
[xData, yData] = prepareCurveData( num_feats, mean_r );

% Set up fittype and options.
ft = fittype( '(1-b*exp(-a*x))+c', 'independent', 'x', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.Lower = [0 -Inf -Inf];
opts.StartPoint = [0.0358727641861837 0.448641991059645 0.835416899160867];
opts.Upper = [1 Inf Inf];

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );

% % Plot fit with data.
hold on;
h = plot( fitresult, xData, yData);
 set(gca,'XScale','log')
 set(h,varargin{:})
 x=h(2).XData;y=h(2).YData;
% legend( h, 'mean_r vs. num_feats', 'untitled fit 1', 'Location', 'NorthEast' );
% Label axes
% xlabel num_feats
% ylabel mean_r
% grid on
clf;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function notDefined = ieNotDefined( varString )
%
% notDefined = ieNotDefined( varString )
%
% Author: ImagEval
% Purpose:
%    Determine if a variable is defined in the calling function's workspace.
%  A variable is defined if (a) it exists and (b) it is not empty.
%
% notDefined: 1 (true) if the variable is not defined in the calling workspace
%             0 (false) if the variable is defined in the calling workspace
%
%  Defined means the variable exists and is not empty in the function that
%  called this function.
%
%  This routine should replace the many calls of the form
%    if ~exist('varname','var') | isempty(xxx), ...
%
%    with the call
%
%    if ieNotDefined('varname')
%
%

if (~ischar(varString)), error('Varible name must be a string'); end

notDefined = 0;  % Assume the variable is defined

str = sprintf('''%s''',varString);
cmd1 = ['~exist(' str ',''var'') == 1'];
cmd2 = ['isempty(',varString ') == 1'];
cmd = [cmd1, ' | ',cmd2];

% If either of these conditions holds, then not defined is true.
notDefined = evalin('caller',cmd1);     % Check that the variable exists in the caller space
if notDefined, return;                  % If it does not, return with a status of 0
else
    notDefined = evalin('caller',cmd2); % Check if the variable is empty in the caller space
    if notDefined return;
    end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [H, pValue, KSstatistic] = kstest_2s_2d(x1, x2, alpha)

% kstest_2s_2d - FUNCTION Two-sample Two-diensional Kolmogorov-Smirnov Test
%
% Usage:[H, pValue, KSstatistic] = kstest_2s_2d(x1, x2 <, alpha>)
%
% The paired-sample Kolmogorov-Smirnov test is a statistical test used to
% determine whether two sets of data arise from the same or different
% distributions.  The null hypothesis is that both data sets were drawn
% from the same continuous distribution.
%
% The algorithm in this function is taken from Peacock [1].
%
% 'x1' is an [Nx2] matrix, each row containing a two-dimensional sample.
% 'x2' is an [Mx2] matrix, each row likewise containing a two-dimensional
% sample. The optional argument 'alpha' is used to set the desired
% significance level for rejecting the null hypothesis.
%
% 'H' is a logical value: true indicates that the null hypothesis should be
% rejected.  'pValue' is an estimate for the P value of the test statistic.
% 'KSstatistic' is the raw value for the test statistic ('D' in [1]).
%
% In contrast to kstest2, this function can only perform a two-tailed test.
% This is because Peacock does not provide a method for estimating P in the
% one-tailed case [1].  Suggestions for a one-tailed test are welcome.
%
% References: [1] J. A. Peacock, "Two-dimensional goodness-of-fit testing
%  in astronomy", Monthly Notices Royal Astronomy Society 202 (1983)
%  615-627.
%    Available from: http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?1983MNRAS.202..615P&defaultprint=YES&filetype=.pdf
%

% Author: Dylan Muir (From kstest_2s_2d by Qiuyan Peng @ ECE/HKUST) Date:
% 13th October, 2012

%%

if nargin < 2
   error('stats:kstest2:TooFewInputs','At least 2 inputs are required.');
end

%%
%
% x1,x2 are both 2-column matrices
%

if ((size(x1,2)~=2)||(size(x2,2)~=2))
   error('stats:kstest2:TwoColumnMatrixRequired','The samples X1 and X2 must be two-column matrices.');
end
n1 = size(x1,1);
n2 = size(x2,1);

%%
%
% Ensure the significance level, ALPHA, is a scalar
% between 0 and 1 and set default if necessary.
%

if (nargin >= 3) && ~isempty(alpha)
   if ~isscalar(alpha) || (alpha <= 0 || alpha >= 1)
      error('stats:kstest2:BadAlpha',...
         'Significance level ALPHA must be a scalar between 0 and 1.');
   end
else
   alpha  =  0.05;
end

%%
%
% Calculate F1(x) and F2(x), the empirical (i.e., sample) CDFs.
%

% - A function handle to perform comparisons in all possible directions
fhCounts = @(x, edge)([(x(:, 1) >= edge(1)) & (x(:, 2) >= edge(2))...
   (x(:, 1) <= edge(1)) & (x(:, 2) >= edge(2))...
   (x(:, 1) <= edge(1)) & (x(:, 2) <= edge(2))...
   (x(:, 1) >= edge(1)) & (x(:, 2) <= edge(2))]);

KSstatistic = -inf;

for iX = 1:(n1+n2)
   % - Choose a starting point
   if (iX<=n1)
      edge = x1(iX,:);
   else
      edge = x2(iX-n1,:);
   end

   % - Estimate the CDFs for both distributions around this point
   vfCDF1 = sum(fhCounts(x1, edge)) ./ n1;
   vfCDF2 = sum(fhCounts(x2, edge)) ./ n2;

   % - Two-tailed test statistic
   vfThisKSTS = abs(vfCDF1 - vfCDF2);
   fKSTS = max(vfThisKSTS);

   % - Final test statistic is the maximum absolute difference in CDFs
   if (fKSTS > KSstatistic)
      KSstatistic = fKSTS;
   end
end

%% Peacock Z calculation and P estimation

n      =  n1 * n2 /(n1 + n2);
Zn = sqrt(n) * KSstatistic;
Zinf = Zn / (1 - 0.53 * n^(-0.9));
pValue = 2 * exp(-2 * (Zinf - 0.5).^2);

% Clip invalid values for P
if (pValue > 0.2)
   pValue = 0.2;
end

H = (pValue <= alpha);


text:
function Yi = qinterp1(x,Y,xi,methodflag)
% Performs fast interpolation compared to interp1
%
% qinterp1 provides a speedup over interp1 but requires an evenly spaced
% x array.  As x and y increase in length, the run-time for interp1 increases
% linearly, but the run-time for
% qinterp1 stays constant.  For small-length x, y, and xi, qinterp1 runs about
% 6x faster than interp1.
%
%
% Usage:
%   yi = qinterp1(x,Y,xi)  - Same usage as interp1
%   yi = qinterp1(x,Y,xi,flag)
%           flag = 0       - Nearest-neighbor
%           flag = 1       - Linear (default)
%
% Example:
%   x = [-5:0.01:5];   y = exp(-x.^2/2);
%   xi = [-4.23:0.3:4.56];
%   yi = qinterp1(x,y,xi,1);
%
% Usage restrictions
%    x must be monotonically and evenly increasing
%    e.g.,  x=-36:0.02:123;
%
%    Y may be up to two-dimensional
%
% Using with non-evenly spaced arrays:
%   Frequently the user will wish to make interpolations "on the fly" from
%   a fixed pair of library (i.e., x and y) vectors.  In this case, the
%   user can generate an equally-spaced set of library data by calling
%   interp1 once, and then storing this library data in a MAT-file or
%   equivalent.  Because the speed of qinterp1 is independent of the length
%   of the library vectors, the author recommends over-sampling this
%   generated set untill memory considerations start limitting program speed.
%
%   If the user wishes to use two or more spacings (i.e., a closely-spaced
%   library in the region of fine features, and a loosely-spaced library in
%   the region of coarse features), just create multiple libraries, record
%   the switching points, and send the search data to different qinterp1
%   calls depending on its value.
%
%   Example:
%       x1 = [-5:0.01:5];   x2 = [-40:1:-5 5:1:40];
%       y1 = exp(-x1.^2/3); y2 = exp(-x2.^2/3);
%       xi = [-30:0.3:30];
%       in = xi < 5 & xi > -5;
%       yi(in) = qinterp1(x1,y1,xi(in));
%       yi(~in) = qinterp1(x2,y2,xi(~in));

% Author: N. Brahms
% Copyright 2006

% Forces vectors to be columns
x = x(:); xi = xi(:);
sx = size(x); sY = size(Y);
if sx(1)~=sY(1)
    if sx(1)==sY(2)
        Y = Y';
    else
        error('x and Y must have the same number of rows');
    end
end

if nargin>=4
    method=methodflag;
else
    method = 1;    % choose nearest-lower-neighbor, linear, etc.
                   % uses integer over string for speed
end

% Gets the x spacing
ndx = 1/(x(2)-x(1)); % one over to perform divide only once
xi = xi - x(1);      % subtract minimum of x

% Fills Yi with NaNs
s = size(Y);
if length(s)>2
    error('Y may only be one- or two-dimensional');
end
Yi = NaN*ones(length(xi),s(2));

switch method
    case 0 %nearest-neighbor method
        rxi = round(xi*ndx)+1;        % indices of nearest-neighbors
        flag = rxi<1 | rxi>length(x) | isnan(xi);
                                      % finds indices out of bounds
        nflag = ~flag;                % finds indices in bounds
        Yi(nflag,:) = Y(rxi(nflag),:);
    case 1 %linear interpolation method
        fxi = floor(xi*ndx)+1;          % indices of nearest-lower-neighbors
        flag = fxi<1 | fxi>length(x)-1 | isnan(xi);
                                        % finds indices out of bounds
        nflag = ~flag;                  % finds indices in bounds
        Yi(nflag,:) = (fxi(nflag)-xi(nflag)*ndx).*Y(fxi(nflag),:)+...
            (1-fxi(nflag)+xi(nflag)*ndx).*Y(fxi(nflag)+1,:);
end


text:
function bsv_test06 ( )

%*****************************************************************************80
%
%% BSV_TEST06 estimates the expected value of the zero crossing.
%
%  Discussion:
%
%    We assume that the left boundary condition ALPHA can vary like
%    a Gaussian variable with mean 1 and standard deviation 0.05.
%
%    Estimate the integral:
%
%      E(X0(ALPHA)) = integral ( -oo < alpha < +oo ) x0(alpha) e^(-0.5*alpha^2/sigma^2) dalpha
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BSV_TEST06:\n' );
  fprintf ( 1, '  For the Burgers equation on [A,B] with viscosity NU and \n' );
  fprintf ( 1, '  boundary conditions U(A)=ALPHA, U(B) = BETA,\n' );
  fprintf ( 1, '  with ALPHA and BETA of opposite sign,\n' );
  fprintf ( 1, '  let X0 be the point where the solution U changes sign.\n' );
  fprintf ( 1, '  Assume ALPHA is Gaussian with mean 0 and standard deviation 0.05.\n' );
  fprintf ( 1, '  Estimate E(X0(ALPHA)) using M Gaussian samples.\n' );

  a = -1.0;
  b = +1.0;
  beta = -1.0;
  nu = 0.1;
  n = 81;
  output = 0;

  x = linspace ( a, b, n );
%
%  Monte Carlo Estimate.
%  Choose M normal random samples with the correct STD,
%  find X0 for each, and average.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '     M    E(X0(ALPHA)) estimate\n' );
  fprintf ( 1, '\n' );

  for j = 4 : 12

    m = 2^j;
    mu = 1.0;
    sigma = 0.05;
    alpha_vec = mu + 0.05 * randn ( m, 1 );

    x0_bar = 0.0;

    for i = 1 : m
      alpha = alpha_vec(i);
      print = 0;
      u = bsv ( a, b, alpha, beta, nu, n, output );
      x0 = bsv_crossing ( a, b, n, x, u );
      x0_bar = x0_bar + x0;
    end

    x0_bar = x0_bar / m;

    fprintf ( 1, '  %4d  %14.6g\n', m, x0_bar );

  end

  return
end


text:
function [cIX,gIX] = GrowClustersFromSeedsItr(thres_merge,thres_cap,thres_minsize,cIX,gIX,M_0)
%disp('find ROIs')

% Set params

% may be determined by hist of distances between all cells (in a sample)
% in correlation distance, i.e. 1 - corr.coeff
% thres_merge = 0.4;
% thres_cap = 0.5;

% thres_minsize = 10; % cell number in final clusters

%%
gIX = SqueezeGroupIX(gIX);
M = M_0(cIX,:);
C = FindCentroid_Direct(gIX,M); % import functional supervoxels ('foxel')
nFoxels = size(C,1);

%% Calculate correlation distance between all cluster-centroids
temp = pdist(C,'correlation');
Dist = squareform(temp);
for i = 1:nFoxels,
    Dist(i,i) = NaN;
end

%% ITERATION:
% initialize
ROI = [];
ROI(nFoxels).fxlist = []; % not capping
ROI(nFoxels).numcell = [];
ROIcount = 0;

% analogy ~ mask to store new ROI
BW = zeros(nFoxels,1);

%%

while true, % loop through all qualified seeds
    % look for next seed
    [~,ix] = min(Dist(:));
    [I,J] = ind2sub(size(Dist),ix);
    % and add seed position to mask
    BW(I) = 1;
    BW(J) = 1;

    if Dist(I,J)>(1-thres_merge),
        break;
    end

    while true % grow ROI: expand mask and examine next neighbors
        IX_in = find(BW);
        IX_out = find(BW==0);

        % find next closest foxel
        D = 1-corr(C(IX_in,:)',C(IX_out,:)'); % distance matrix between potential foxels to foxels already included in ROI
        if isempty(D),
            break;
        end
        D2 = min(D,[],1);
        [a,ix] = min(D2);
        ix_pretend = IX_out(ix);

        BW_pretend = BW;
        BW_pretend(ix_pretend) = 1;

        % find correlation between potential foxel and new core
        list = find(BW_pretend);
        IX = [];
        for i_gIX = 1:length(list),
            IX = [IX;find(gIX==list(i_gIX))];
        end
        M_core = M(IX,:);
        [~,newCore] = kmeans(M_core,1,'distance','correlation');
        coredist = 1-corr(C(ix_pretend,:)',newCore');

        if a<(1-thres_merge) && coredist<(1-thres_cap),
            BW(ix_pretend) = 1;
        else
            break; % finished expanding this seed
        end
    end

    % save this ROI if bigger than thres
    list = find(BW);
    IX = [];
    for i_gIX = 1:length(list),
        IX = [IX;find(gIX==list(i_gIX))]; %#ok<AGROW>
    end
    numcell = length(IX);
    if numcell >= thres_minsize,
        ROIcount = ROIcount+1;
        ROI(ROIcount).fxlist = find(BW);
        ROI(ROIcount).numcell = numcell;
    end

    % update/reset
    Dist(list,:) = NaN;
    Dist(:,list) = NaN;
    C(list,:) = NaN;
    BW = zeros(nFoxels,1);

    if isempty(find(isnan(Dist)==0)),
        break;
    end
end

ROI(ROIcount+1:end) = [];

%% Merge accordingly
for i = 1:ROIcount,
    fxlist = ROI(i).fxlist;
    for j = 2:length(fxlist),
        gIX(gIX==fxlist(j)) = fxlist(1);
    end
end
[gIX,numU] = SqueezeGroupIX(gIX);
%disp(numU);

end

author: xiuyechen

repo: FishExplorer


text:
function [dietPerformance] = compareDietPerformance(model, diets,rois, varargin)
% Function to check the effect of different diets on the possible flux
% range for specified reactions.
%
% [dietPerformance] = compareDietPerformance(model, diets, rois, varargin)
%
% Example: [dietPerformance] = compareDietPerformance(model, diets, rois,
% 'optPercentage', 0.95, 'microbiotaExcretionBound', 0.5)
%
% INPUTS:
%   model:          COBRA model where the rois are present
%   diets:          Cell array where each cell contains an cell array with
%                   The metabolite composition of the diet {'met', value}
%   rois:           Cell array with reactions to investigate under
%                   different diets
%
% OPTIONAL INPUTS:
%   optPercentage:  Integer, used in fluxVariability to indicate the
%                   percentage of the objective for which FVA is solved.
%                   Defaults to 99
%   microbiotaExretionBound: Integer, used to constrain both bounds of
%                   Excretion_EX_microbiota_LI_biomass. Defaults to 1
%   dietNames:      Cell array with names for diets in that will be used in
%                   the output. Defaults to Diet1, Diet2 etcetc.
% OUTPUT:
%   dietPerformance:    Table for each given diet, the min and max flux for
%                   each roi.
%
% .. Authors: - Bram Nap 06-2022

% Initialise the inputs
parser = inputParser();
parser.addRequired('model', @isstruct);
parser.addRequired('diets', @iscell);
parser.addRequired('rois', @iscell);
parser.addParameter('optPercentage', 99, @isnumeric);
parser.addParameter('microbiotaExcretionBound', 1, @isnumeric);
parser.addParameter('dietNames', cell(0), @iscell);

parser.parse(model, diets, rois, varargin{:});

model = parser.Results.model;
diets = parser.Results.diets;
rois = parser.Results.rois;
optPercentage = parser.Results.optPercentage;
microbiotaExcretionBound = parser.Results.microbiotaExcretionBound;
dietNames = parser.Results.dietNames;

% If the dietNames are not specified, create general names
if isempty(dietNames)
    for i = 1:length(diets)
        dietNames(i) = {['Diet' + string(i)]};
    end
    % Required to prevent errors
    dietNames = string(dietNames);
end

% Initialise array to store results in
dietPerformance = cell(2*length(rois), length(diets)+2);

% Loop over the diets
for i = 1:length(diets)
    diet = diets{i};

    % Constrain the model to the diet and microbiota exchange bounds
    modelDiet = setFoodConstraints(model, diet);
    modelDiet = changeRxnBounds(modelDiet, 'Excretion_EX_microbiota_LI_biomass',microbiotaExcretionBound,'b');

    % Solve with FVA and store results
    [minFlux, maxFlux] = fluxVariability(modelDiet, optPercentage, 'max', rois);
    for j = 1:length(rois)
        dietPerformance(j*2-1:j*2,1) = {rois(j)};
        dietPerformance(j*2-1:j*2,2) = {{'min'};{'max'}};
        dietPerformance(j*2-1:j*2,i+2) = {minFlux(j); maxFlux(j)};
    end
end

% Create the final table
dietPerformance = cell2table(dietPerformance);
dietPerformance.Properties.VariableNames(1:2) = {'Reactions', 'MinMax'};
dietPerformance.Properties.VariableNames(3:end) = dietNames;

end

author: opencobra

repo: cobratoolbox


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Function (ldiv) : calculate inverse Z-transform by long division
% Author : Tamer mohamed samy abdelazim Mellik
% Contact information :
%Department of Electrical & Computer Engineering,
%University of Calgary,
%2500 University Drive N.W. ,
%Calgary, AB T2N 1N4 ,
%Canada .
% email :abdelasi@enel.ucalgary.ca
% email : tabdelaz@ucalgary.ca
% Webpage : http://www.enel.ucalgary.ca/~abdelasi/
% Date    : 2-5-2002
% Version : 1.0.0
%Example
% This function like deconv but it help if the numerator less or equal degree of denominator
% if you have this function (It must arranged in terms of minus power of Z):
%              1
% G(z)= -----------------
%             -1     -2
%      ( 5 - Z  - 3 Z  )
% and you want to calculate long division or inverse Z transform :
% The numerator  is a=[1] and  the denominator is b= [5 -1 -3 ]
% call the function ldiv(a,b) to get the funresult 20 items (default)
% another example :
%                -2     -3
%       ( 5 - 3 Z  + 4 Z  )
%  G(z)-----------------
%              -1     -2
%        (5 - Z  - 3 Z )
% a=[5 0 -3 4] , b= [5 -1 -3 ] and you want the funresult 100 terms !
% ldiv(a,b,100)
% Note : The author doesn't have any responsibility for any harm caused by the use of this file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function funresult=ldiv(a,b,N)
%a numerator
%b denominator
%default order of the filter == 20
funresult=[];
if nargin < 3
    if nargin > 1
            N=20;
    else
            disp('Usage: M = ldiv(a,b,N)')
            disp('a:numerator , b denominator and N is the order of the resultant filter')
            return
    end
end

if size(a) < 1
    disp('Error: numerator must at least have one element not empty')
    return
end
if size(b) < 1
    disp('Error: denominator must at least have one element not empty')
    return
end

if b(1)==0
    disp('Error: The first element of denominator must have nonzero value')
    return
end
if size(b) < 2
    funresult=a./b;
    for i =length(funresult)+1:N
       funresult(i)=0;
    end
    return
end

for i = length(a)+1:N
    a(i)=0;
end
for i = 1 : N
    funresult(i)=a(1)/b(1);
    if length(a)>1
            for k= 2:length(b)
                if k > length(a)
                    a(k)=0;
                end
                a(k)=a(k)-funresult(length(funresult))*b(k);
            end

        for i = 1:length(a)-1
                a(i)=a(i+1);
        end
        a=a(1:length(a)-1);
    end
end


text:
function [p,y,a] = qint(ym1,y0,yp1)
%QINT   Quadratic interpolation of 3 uniformly spaced samples
%
%               [p,y,a] = qint(ym1,y0,yp1)
%
%       returns extremum-location p, height y, and curvature a
%       of a parabolic fit through three points.
%       The parabola is given by y(x) = a*(x-p)^2+b,
%       where y(-1)=ym1, y(0)=y0, y(1)=yp1.

   p = (yp1 - ym1)/(2*(2*y0 - yp1 - ym1));
   if nargout>1
     y = y0 - 0.25*(ym1-yp1)*p;
   end;
   if nargout>2
     a = 0.5*(ym1 - 2*y0 + yp1);
   end;

repo: REST


text:
function [p,f] = oct3bank(x);
% OCT3BANK Simple one-third-octave filter bank.
%    OCT3BANK(X) plots one-third-octave power spectra of signal vector X.
%    Implementation based on ANSI S1.11-1986 Order-3 filters.
%    Sampling frequency Fs = 44100 Hz. Restricted one-third-octave-band
%    range (from 100 Hz to 5000 Hz). RMS power is computed in each band
%    and expressed in dB with 1 as reference level.
%
%    [P,F] = OCT3BANK(X) returns two length-18 row-vectors with
%    the RMS power (in dB) in P and the corresponding preferred labeling
%    frequencies (ANSI S1.6-1984) in F.
%
%    See also OCT3DSGN, OCT3SPEC, OCTDSGN, OCTSPEC.

% Author: Christophe Couvreur, Faculte Polytechnique de Mons (Belgium)
%         couvreur@thor.fpms.ac.be
% Last modification: Aug. 23, 1997, 10:30pm.

% References:
%    [1] ANSI S1.1-1986 (ASA 65-1986): Specifications for
%        Octave-Band and Fractional-Octave-Band Analog and
%        Digital Filters, 1993.
%    [2] S. J. Orfanidis, Introduction to Signal Processing,
%        Prentice Hall, Englewood Cliffs, 1996.

pi = 3.14159265358979;
Fs = 44100; 				% Sampling Frequency
N = 3; 					% Order of analysis filters.
F = [ 100 125 160, 200 250 315, 400 500 630, 800 1000 1250, ...
	1600 2000 2500, 3150 4000 5000 ]; % Preferred labeling freq.
ff = (1000).*((2^(1/3)).^[-10:7]); 	% Exact center freq.
P = zeros(1,18);
m = length(x);

% Design filters and compute RMS powers in 1/3-oct. bands
% 5000 Hz band to 1600 Hz band, direct implementation of filters.
for i = 18:-1:13
   [B,A] = oct3dsgn(ff(i),Fs,N);
   y = filter(B,A,x);
   P(i) = sum(y.^2)/m;
end
% 1250 Hz to 100 Hz, multirate filter implementation (see [2]).
[Bu,Au] = oct3dsgn(ff(15),Fs,N); 	% Upper 1/3-oct. band in last octave.
[Bc,Ac] = oct3dsgn(ff(14),Fs,N); 	% Center 1/3-oct. band in last octave.
[Bl,Al] = oct3dsgn(ff(13),Fs,N); 	% Lower 1/3-oct. band in last octave.
for j = 3:-1:0
   x = decimate(x,2);
   m = length(x);
   y = filter(Bu,Au,x);
   P(j*3+3) = sum(y.^2)/m;
   y = filter(Bc,Ac,x);
   P(j*3+2) = sum(y.^2)/m;
   y = filter(Bl,Al,x);
   P(j*3+1) = sum(y.^2)/m;
end

% Convert to decibels.
Pref = 1; 				% Reference level for dB scale.
idx = (P>0);
P(idx) = 10*log10(P(idx)/Pref);
P(~idx) = NaN*ones(sum(~idx),1);

% Generate the plot
if (nargout == 0)
  bar(P);
  ax = axis;
  axis([0 19 ax(3) ax(4)])
  set(gca,'XTick',[2:3:18]); 		% Label frequency axis on octaves.
  set(gca,'XTickLabels',F(2:3:length(F)));  % MATLAB 4.1c
%  set(gca,'XTickLabel',F(2:3:length(F)));  % MATLAB 5.1
  xlabel('Frequency band [Hz]'); ylabel('Power [dB]');
  title('One-third-octave spectrum')
% Set up output parameters
elseif (nargout == 1)
  p = P;
elseif (nargout == 2)
  p = P;
  f = F;
end

repo: mcbench-benchmarks


text:
function plotRendezvous(hAxis, timeArr, bodyInfo, iniOrbit, finOrbit, xfrOrbit, dv1, dv2, celBodyData)
%plotRendezvous Summary of this function goes here
%   Detailed explanation goes here

    gmuXfr = bodyInfo.gm;

    cla(hAxis, 'reset');

    iniOrbBodyInfo = getBodyInfoStructFromOrbit(iniOrbit);
    finOrbBodyInfo = getBodyInfoStructFromOrbit(finOrbit);

    xfrOrbit2 = xfrOrbit;
    xfrOrbit2(7) = timeArr(1);
    xfrOrbit2(6) = computeMeanFromTrueAnom(xfrOrbit(6), xfrOrbit(2));
    xfrOrbBodyInfo = getBodyInfoStructFromOrbit(xfrOrbit2);

    axes(hAxis);
    hold on;
    axes(hAxis);
    set(hAxis,'XTickLabel',[]);
    set(hAxis,'YTickLabel',[]);
    set(hAxis,'ZTickLabel',[]);
    grid on;
    view(3);

%     dRad = bodyInfo.radius;
%     [X,Y,Z] = sphere(30);
%     surf(dRad*X,dRad*Y,dRad*Z);
%     colormap(bodyInfo.bodycolor);

    ma_initOrbPlot([], hAxis, bodyInfo);

    parentBodyInfo = bodyInfo.getParBodyInfo(celBodyData);

    [iniLB, iniUB, ~] = getOrbitTAPlotBnds(bodyInfo, parentBodyInfo, iniOrbit);
    [finLB, finUB, ~] = getOrbitTAPlotBnds(bodyInfo, parentBodyInfo, finOrbit);

    time1 = timeArr(1);
    time2 = timeArr(2);

    if(iniOrbit(2) >= 1)
        [rVect, vVect] = getStateAtTime(iniOrbBodyInfo, time1, iniOrbit(8));
        [~, ~, ~, ~, ~, iniTA1] = getKeplerFromState(rVect,vVect,iniOrbit(8));
        iniLB = min(iniLB,iniTA1);

        [rVect, vVect] = getStateAtTime(iniOrbBodyInfo, time2, finOrbit(8));
        [~, ~, ~, ~, ~, iniTA2] = getKeplerFromState(rVect,vVect,finOrbit(8));
        iniUB = max(iniUB,iniTA2);
    end

    if(finOrbit(2) >= 1)
        [rVect, vVect] = getStateAtTime(finOrbBodyInfo, time1, finOrbit(8));
        [~, ~, ~, ~, ~, finTA1] = getKeplerFromState(rVect,vVect,finOrbit(8));
        finLB = min(finLB,finTA1);

        [rVect, vVect] = getStateAtTime(finOrbBodyInfo, time2, finOrbit(8));
        [~, ~, ~, ~, ~, finTA2] = getKeplerFromState(rVect,vVect,finOrbit(8));
        finUB = max(finUB,finTA2);
    end

    plotOrbit('r', iniOrbit(1), iniOrbit(2), iniOrbit(3), iniOrbit(4), iniOrbit(5), iniLB, iniUB, gmuXfr);
    plotOrbit('b', finOrbit(1), finOrbit(2), finOrbit(3), finOrbit(4), finOrbit(5), finLB, finUB, gmuXfr);
    plotOrbit('k', xfrOrbit(1), xfrOrbit(2), xfrOrbit(3), xfrOrbit(4), xfrOrbit(5), xfrOrbit(6), xfrOrbit(7), gmuXfr,[],[],[],'--');

    [iniRVectT1, ~] = getStateAtTime(iniOrbBodyInfo, time1, gmuXfr);
    [iniRVectT2, ~] = getStateAtTime(iniOrbBodyInfo, time2, gmuXfr);
    [finRVectT1, ~] = getStateAtTime(finOrbBodyInfo, time1, gmuXfr);
    [finRVectT2, ~] = getStateAtTime(finOrbBodyInfo, time2, gmuXfr);

    hold on;
    plot3(iniRVectT1(1),iniRVectT1(2),iniRVectT1(3),...
                    '-ro',...
                    'LineWidth',1.0,...
                    'MarkerEdgeColor','k',...
                    'MarkerFaceColor','r',...
                    'MarkerSize',10);

    hold on;
    plot3(iniRVectT2(1),iniRVectT2(2),iniRVectT2(3),...
                    '-ro',...
                    'LineWidth',1.0,...
                    'MarkerEdgeColor','k',...
                    'MarkerFaceColor','r',...
                    'MarkerSize',10);

    hold on;
    plot3(finRVectT1(1),finRVectT1(2),finRVectT1(3),...
                    '-bo',...
                    'LineWidth',1.0,...
                    'MarkerEdgeColor','k',...
                    'MarkerFaceColor','b',...
                    'MarkerSize',10);

    hold on;
    plot3(finRVectT2(1),finRVectT2(2),finRVectT2(3),...
                    '-bo',...
                    'LineWidth',1.0,...
                    'MarkerEdgeColor','k',...
                    'MarkerFaceColor','b',...
                    'MarkerSize',10);

    hold off;
    axis equal;

    drawnow;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef PTKThresholdLungFiltered < PTKPlugin
    % PTKThresholdLungFiltered. Plugin to detect airlike voxels using thresholding.
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %     PTKThresholdLungFiltered uses the library routine PTKThresholdAirway to
    %     detect air-like voxels in the lung. Before thresholding, a Gaussian
    %     filter is applied to the image.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Lung threshold <br>(filtered)'
        ToolTip = 'Detect air-like voxels in the lung through thresholding'
        Category = 'Lungs'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Visibility = 'Developer'
        Version = 2
    end

    methods (Static)
        function threshold_image = RunPlugin(dataset, ~)
            if dataset.IsGasMRI
                lung_roi = dataset.GetResult('PTKInvertImage');
                filtered_image = MimGaussianFilter(lung_roi, 2);
            else
                lung_roi = dataset.GetResult('PTKLungROI');
                filtered_image = MimGaussianFilter(lung_roi, 0.5, true);
            end

            if dataset.IsGasMRI
                mri_lung_threshold = dataset.GetResult('PTKMRILungThreshold');
                limits = mri_lung_threshold.Bounds;
                raw_image = filtered_image.RawImage;
                raw_image = (raw_image >= limits(1) & raw_image <= limits(2));
                threshold_image = lung_roi.BlankCopy;
                threshold_image.ChangeRawImage(raw_image);
                threshold_image.ImageType = PTKImageType.Colormap;
            elseif strcmp(dataset.GetImageInfo.Modality, 'MR')
                mri_lung_threshold = dataset.GetResult('PTKMRILungThreshold');
                limits = mri_lung_threshold.Bounds;
                raw_image = filtered_image.RawImage;
                raw_image = (raw_image >= limits(1) & raw_image <= limits(2));
                threshold_image = lung_roi.BlankCopy;
                threshold_image.ChangeRawImage(raw_image);
                threshold_image.ImageType = PTKImageType.Colormap;
            else
                limit_1 = lung_roi.RescaledToGreyscale(-1024);
                limit_2 = lung_roi.RescaledToGreyscale(-775);
                limit_3 = lung_roi.RescaledToGreyscale(-400);

                % Voxles within the wider filtered threshold are given value 2
                filtered_image = uint8(2*(filtered_image.RawImage >= limit_1 & filtered_image.RawImage <= limit_3));

                % Voxles within the narrow unfiltered threshold are given value 1
                filtered_image(lung_roi.RawImage >= limit_1 & lung_roi.RawImage <= limit_2) = 1;

                threshold_image = lung_roi.BlankCopy;
                threshold_image.ChangeRawImage(filtered_image);

                threshold_image.ImageType = PTKImageType.Colormap;
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function varargout=phaseplot(f,varargin)
%PHASEPLOT  Phase plot
%   Usage: phaseplot(f,op1,op2, ... );
%          phaseplot(f,fs,op1,op2, ... );
%
%   `phaseplot(f)` plots the phase of f using a |dgt|.
%
%   `phaseplot(f,fs)` does the same for a signal with sampling rate *fs* Hz.
%
%   `phaseplot` should only be used for short signals (shorter than the
%   resolution of the screen), as there will otherwise be some visual
%   aliasing, such that very fast changing areas will look very smooth.
%   `phaseplot` always calculates the phase of the full time/frequency plane
%   (as opposed to |sgram|), and you therefore risk running out of memory
%   for long signals.
%
%   `phaseplot` takes the following flags at the end of the line of input
%   arguments:
%
%     'tfr',v     Set the ratio of frequency resolution to time resolution.
%                 A value $v=1$ is the default. Setting $v>1$ will give better
%                 frequency resolution at the expense of a worse time
%                 resolution. A value of $0<v<1$ will do the opposite.
%
%     'wlen',s    Window length. Specifies the length of the window
%                 measured in samples. See help of |pgauss| on the exact
%                 details of the window length.
%
%     'nf'        Display negative frequencies, with the zero-frequency
%                 centered in the middle. For real signals, this will just
%                 mirror the upper half plane. This is standard for complex
%                 signals.
%
%     'tc'        Time centering. Move the beginning of the signal to the
%                 middle of the plot. This is usefull for visualizing the
%                 window functions of the toolbox.
%
%     'thr',r     Keep the coefficients with a magnitude larger than r times the
%                 largest magnitude. Set the phase of the rest of the
%                 coefficients to zero. This is useful, because for small
%                 amplitude the phase values can be meaningless.
%
%     'timeinv'   Display the phase as computed by a time-invariant
%                 |dgt|. This is the default.
%
%     'freqinv'   Display the phase as computed by a frequency-invariant
%                 |dgt|.
%
%     'fmax',y    Display y as the highest frequency.
%
%     'colorbar'   Display the colorbar. This is the default.
%
%     'nocolorbar'  Do not display the colorbar.
%
%   For the best result when using `phaseplot`, use a circulant color
%   map, for instance `hsv`.
%
%   Examples:
%   ---------
%
%   The following code shows the phaseplot of a
%   periodic, hyperbolic secant visualized using the `hsv` colormap:::
%
%     phaseplot(psech(200),'tc','nf');
%     colormap(hsv);
%
%   The following phaseplot shows the phase of white, Gaussian noise:::
%
%     phaseplot(randn(200,1));
%     colormap(hsv);
%
%   See also: phaselock
%
%   Demos: demo_phaseplot
%
%   References: Carmona98practical

%   AUTHOR: Peter L. Søndergaard
%   TESTING: NA

if nargin<1
  error('Too few input arguments.');
end;

if sum(size(f)>1)>1
  error('Input must be a vector.');
end;

definput.import={'ltfattranslate','setnorm','tfplot'};
% Override the setting from tfplot, because phaseplot only uses the 'lin'
% plotting.
definput.flags.log={'lin'};
definput.importdefaults={'lin'};

% Define initial value for flags and key/value pairs.
definput.flags.wlen={'nowlen','wlen'};
definput.flags.tc={'notc','tc'};

definput.flags.fmax={'nofmax','fmax'};
definput.flags.phase={'timeinv','freqinv'};

if isreal(f)
  definput.flags.posfreq={'posfreq','nf'};
else
  definput.flags.posfreq={'nf','posfreq'};
end;

definput.keyvals.tfr=1;
definput.keyvals.wlen=0;
definput.keyvals.fmax=[];
definput.keyvals.thr=[];

[flags,kv,fs]=ltfatarghelper({'fs'},definput,varargin);

% Downsample
if ~isempty(kv.fmax)
  if ~isempty(fs)
    resamp=kv.fmax*2/fs;
  else
    resamp=kv.fmax*2/length(f);
  end;

  f=fftresample(f,round(length(f)*resamp));
  kv.fs=2*kv.fmax;
end;

% Always do the full STFT
L=length(f);
a=1;
b=1;
M=L;
N=L;

% Set an explicit window length, if this was specified.
if flags.do_wlen
  kv.tfr=kv.wlen^2/L;
end;

g={'gauss',kv.tfr,flags.norm};

if flags.do_nf
  coef=dgt(f,g,a,M,flags.phase);
else
  coef=dgtreal(f,g,a,M,flags.phase);
end;

if ~isempty(kv.thr)
  % keep only the largest coefficients.
  maxc=max(abs(coef(:)));
  mask=abs(coef)<maxc*kv.thr;
  coef(mask)=0;
end

coef = angle(coef);

if flags.do_nf
  plotdgt(coef,a,'argimport',flags,kv);
else
  plotdgtreal(coef,a,M,'argimport',flags,kv);
end;

if nargout>0
  varargout={coef};
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat