text:
%SerialLink.ikine Numerical inverse kinematics
%
% Q = R.ikine(T) are the joint coordinates (1xN) corresponding to the robot
% end-effector pose T (4x4) which is a homogenenous transform.
%
% Q = R.ikine(T, Q0, OPTIONS) specifies the initial estimate of the joint
% coordinates.
%
% This method can be used for robots with 6 or more degrees of freedom.
%
% Underactuated robots::
%
% For the case where the manipulator has fewer than 6 DOF the solution
% space has more dimensions than can be spanned by the manipulator joint
% coordinates.
%
% Q = R.ikine(T, Q0, M, OPTIONS) similar to above but where M is a mask
% vector (1x6) which specifies the Cartesian DOF (in the wrist coordinate
% frame) that will be ignored in reaching a solution.  The mask vector
% has six elements that correspond to translation in X, Y and Z, and rotation
% about X, Y and Z respectively.  The value should be 0 (for ignore) or 1.
% The number of non-zero elements should equal the number of manipulator DOF.
%
% For example when using a 3 DOF manipulator rotation orientation might be
% unimportant in which case  M = [1 1 1 0 0 0].
%
% For robots with 4 or 5 DOF this method is very difficult to use since
% orientation is specified by T in world coordinates and the achievable
% orientations are a function of the tool position.
%
% Trajectory operation::
%
% In all cases if T is 4x4xM it is taken as a homogeneous transform sequence
% and R.ikine() returns the joint coordinates corresponding to each of the
% transforms in the sequence.  Q is MxN where N is the number of robot joints.
% The initial estimate of Q for each time step is taken as the solution
% from the previous time step.
%
% Options::
% 'pinv'         use pseudo-inverse instead of Jacobian transpose (default)
% 'ilimit',L     set the maximum iteration count (default 1000)
% 'tol',T        set the tolerance on error norm (default 1e-6)
% 'alpha',A      set step size gain (default 1)
% 'varstep'      enable variable step size if pinv is false
% 'verbose'      show number of iterations for each point
% 'verbose=2'    show state at each iteration
% 'plot'         plot iteration state versus time
%
% References::
% - Robotics, Vision & Control, Section 8.4,
%   P. Corke, Springer 2011.
%
% Notes::
% - Solution is computed iteratively.
% - Solution is sensitive to choice of initial gain.  The variable
%   step size logic (enabled by default) does its best to find a balance
%   between speed of convergence and divergence.
% - Some experimentation might be required to find the right values of
%   tol, ilimit and alpha.
% - The pinv option leads to much faster convergence (default)
% - The tolerance is computed on the norm of the error between current
%   and desired tool pose.  This norm is computed from distances
%   and angles without any kind of weighting.
% - The inverse kinematic solution is generally not unique, and
%   depends on the initial guess Q0 (defaults to 0).
% - The default value of Q0 is zero which is a poor choice for most
%   manipulators (eg. puma560, twolink) since it corresponds to a kinematic
%   singularity.
% - Such a solution is completely general, though much less efficient
%   than specific inverse kinematic solutions derived symbolically, like
%   ikine6s or ikine3.
% - This approach allows a solution to be obtained at a singularity, but
%   the joint angles within the null space are arbitrarily assigned.
% - Joint offsets, if defined, are added to the inverse kinematics to
%   generate Q.
% - Joint limits are not considered in this solution.
%
% See also SerialLink.ikcon, SerialLink.ikunc, SerialLink.fkine, SerialLink.ikine6s.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function [qt,histout] = ikine(robot, tr, varargin)

    %  set default parameters for solution
    opt.ilimit = 1000;
    opt.tol = 1e-6;
    opt.alpha = 0.9;
    opt.plot = false;
    opt.pinv = true;
    opt.varstep = false;

    [opt,args] = tb_optparse(opt, varargin);

    n = robot.n;

    % robot.ikine(tr, q)
    if ~isempty(args)
        q = args{1};
        if isempty(q)
            q = zeros(1, n);
        else
            q = q(:)';
        end
    else
        q = zeros(1, n);
    end

    % robot.ikine(tr, q, m)
    if length(args) > 1
        m = args{2};
        m = m(:);
        if numel(m) ~= 6
            error('RTB:ikine:badarg', 'Mask matrix should have 6 elements');
        end
        if numel(find(m)) > robot.n
            error('RTB:ikine:badarg', 'Number of robot DOF must be >= the same number of 1s in the mask matrix')
        end
    else
        if n < 6
            error('RTB:ikine:badarg', 'For a manipulator with fewer than 6DOF a mask matrix argument must be specified');
        end
        m = ones(6, 1);
    end
    % make this a logical array so we can index with it
    m = logical(m);

    npoints = size(tr,3);    % number of points
    qt = zeros(npoints, n);  % preallocate space for results
    tcount = 0;              % total iteration count

    if ~ishomog(tr)
        error('RTB:ikine:badarg', 'T is not a homog xform');
    end

    J0 = jacob0(robot, q);
    J0 = J0(m, :);
    if cond(J0) > 100
        warning('RTB:ikine:singular', 'Initial joint configuration results in a (near-)singular configuration, this may slow convergence');
    end

    history = [];
    failed = false;
    e = zeros(6,1);
    revolutes = [robot.links.sigma] == 0;

    for i=1:npoints
        T = tr(:,:,i);

        nm = Inf;
        % initialize state for the ikine loop
        eprev = [Inf Inf Inf Inf Inf Inf];
        save.e = [Inf Inf Inf Inf Inf Inf];
        save.q = [];
        count = 0;

        while true
            % update the count and test against iteration limit
            count = count + 1;
            if count > opt.ilimit
                warning('ikine: iteration limit %d exceeded (row %d), final err %f', ...
                    opt.ilimit, i, nm);
                failed = true;
                %q = NaN*ones(1,n);
                break
            end

            % compute the error
            Tq = robot.fkine(q');

            e(1:3) = transl(T - Tq);
            Rq = t2r(Tq);
            [th,n] = tr2angvec(Rq'*t2r(T));
            e(4:6) = th*n;

            % optionally adjust the step size
            if opt.varstep
                % test against last best error, only consider the DOF of
                % interest
                if norm(e(m)) < norm(save.e(m))
                    % error reduced,
                    % let's save current state of solution and rack up the step size
                    save.q = q;
                    save.e = e;
                    opt.alpha = opt.alpha * (2.0^(1.0/8));
                    if opt.verbose > 1
                        fprintf('step %d: raise alpha to %f\n', count, opt.alpha);
                    end
                else
                    % rats!  error got worse,
                    % restore to last good solution and reduce step size
                    q = save.q;
                    e = save.e;
                    opt.alpha = opt.alpha * 0.5;
                    if opt.verbose > 1
                        fprintf('step %d: drop alpha to %f\n', count, opt.alpha);
                    end
                end
            end

            % compute the Jacobian
            J = jacob0(robot, q);

            % compute change in joint angles to reduce the error,
            % based on the square sub-Jacobian
            if opt.pinv
                %pinv( J(m,:) )
                dq = opt.alpha * pinv( J(m,:) ) * e(m);
            else
                dq = J(m,:)' * e(m);
                dq = opt.alpha * dq;
            end

            % diagnostic stuff
            if opt.verbose > 1
                fprintf('%d/%d: |e| = %f\n', i, count, nm);
                fprintf('       e  = '); disp(e');
                fprintf('       dq = '); disp(dq');
            end
            if opt.plot
                h.q = q';
                h.dq = dq;
                h.e = e;
                h.ne = nm;
                h.alpha = opt.alpha;
                history = [history; h]; %#ok<*AGROW>
            end

            % update the estimated solution
            q = q + dq';

            % wrap angles for revolute joints
            k = (q > pi) & revolutes;
            q(k) = q(k) - 2*pi;

            k = (q < -pi) & revolutes;
            q(k) = q(k) + 2*pi;

            nm = norm(e(m));

            if norm(e(m)) > 2*norm(eprev(m))
                warning('RTB:ikine:diverged', 'solution diverging at step %d, try reducing alpha', count);
            end
            eprev = e;

            if nm <= opt.tol
                break
            end

        end  % end ikine solution for tr(:,:,i)
        qt(i,:) = q';
        tcount = tcount + count;
        if opt.verbose
            fprintf('%d iterations\n', count);
        end
    end

    if opt.verbose && npoints > 1
        fprintf('TOTAL %d iterations\n', tcount);
    end

    % plot evolution of variables
    if opt.plot
        figure(1);
        subplot(511)
        plot([history.q]');
        ylabel('q');
        grid

        subplot(512)
        plot([history.dq]');
        ylabel('dq');
        grid

        subplot(513)
        plot([history.e]');
        ylabel('e');
        grid

        subplot(514)
        semilogy([history.ne]);
        ylabel('|e|');
        grid

        subplot(515)
        plot([history.alpha]);
        xlabel('iteration');
        ylabel('\alpha');
        grid

        if nargout > 1
            histout = history;
        end
    end
end

author: Allopart

repo: rbpf-gmapping


text:
% MATLAB script that add the PhysioNet Cardiovascular Signal Toolbox to
% Matlab path, aslo check Matlab version

% --- Check Matlab version
if verLessThan('matlab', '9.3.0')
    warning('Matlab:Version', 'Your version of Matlab is too old. R2017b or higher is requested.\n');
    return;
end

fprintf('Adding the PhysioNet Cardiovascular Signal Toolbox to Matlab path\n')
try
    ss = ['..' filesep 'PhysioNet-Cardiovascular-Signal-Toolbox'];
    addpath(genpath(ss));
    fprintf('PhysioNet Cardiovascular Signal Toolbox successfully added to Matlab path\n')
catch
    try
        ss = ['..' filesep 'PhysioNet-Cardiovascular-Signal-Toolbox-master'];
        addpath(genpath(ss));
        fprintf('PhysioNet Cardiovascular Signal Toolbox successfully added to Matlab path\n')
    catch
    end
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
% [MIN, MAX] = range2(MTX)
%
% Compute minimum and maximum values of MTX, returning them as a 2-vector.

% Eero Simoncelli, 3/97.

function [mn, mx] = range2(mtx)

%% NOTE: THIS CODE IS NOT ACTUALLY USED! (MEX FILE IS CALLED INSTEAD)

% fprintf(1,'WARNING: You should compile the MEX code for "range2", found in the MEX subdirectory.  It is MUCH faster.\n');

if (~isreal(mtx))
  error('MTX must be real-valued');
end

% mn = min(min(mtx));
% mx = max(max(mtx));
mn = min(mtx(:));
mx = max(mtx(:));

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ellipse, labels] = imInertiaEllipse(img, varargin)
%IMINERTIAELLIPSE Inertia ellipse of a binary or label image
%
%   Deprecated, use 'imEquivalentEllipse' instead.
%
%   ELLI = imInertiaEllipse(IMG)
%   Compute the inertia ellipses of the particles in labeled image IMG. If
%   the image is binary, one ellipse, corresponding to the foreground (i.e.
%   the pixels with value 1) will be computed.
%
%   The result is a N-by-5 array ELLI = [XC YC A B THETA], containing
%   coordinates of ellipse center, lengths of semi major and minor axes,
%   and the orientation (given in degrees, in the direction of the greatest
%   axis).
%
%   The same result could be obtained with the regionprops function. The
%   advantage of using imInertiaEllipse is that equivalent ellipses can be
%   obtained in one call. Orientation of both functions are not consistent.
%
%   ELLI = imInertiaEllipse(IMG, SPACING);
%   ELLI = imInertiaEllipse(IMG, SPACING, ORIGIN);
%   Specifies the spatial calibration of image. Both SPACING and ORIGIN are
%   1-by-2 row vectors. SPACING = [SX SY] contains the size of a pixel.
%   ORIGIN = [OX OY] contains the center position of the top-left pixel of
%   image.
%   If no calibration is specified, spacing = [1 1] and origin = [1 1] are
%   used. If only the sapcing is specified, the origin is set to [0 0].
%
%   ELLI = imInertiaEllipse(..., LABELS)
%   Specify the labels for which the inertia ellipse needs to be computed.
%   The result is a N-by-5 array with as many rows as the number of labels.
%
%
%   Example
%   % Draw a commplex particle together with its equivalent ellipse
%     img = imread('circles.png');
%     imshow(img); hold on;
%     elli = imInertiaEllipse(img);
%     drawEllipse(elli)
%
%   % Compute and display the equivalent ellipses of several particles
%     img = imread('rice.png');
%     img2 = img - imopen(img, ones(30, 30));
%     lbl = bwlabel(img2 > 50, 4);
%     ellipses = imInertiaEllipse(lbl);
%     imshow(img); hold on;
%     drawEllipse(ellipses, 'linewidth', 2, 'color', 'g');
%
%   See also
%     imEquivalentEllipse
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2011-03-30,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

warning('MatImage:deprecated', ...
    'function imInertiaEllipse is obsolete, use imEquivalentEllipse instead');

%% extract spatial calibration

% default values
spacing = [1 1];
origin  = [1 1];
calib   = false;

% extract spacing
if ~isempty(varargin) && sum(size(varargin{1}) == [1 2]) == 2
    spacing = varargin{1};
    varargin(1) = [];
    calib = true;
    origin = [0 0];
end

% extract origin
if ~isempty(varargin) && sum(size(varargin{1}) == [1 2]) == 2
    origin = varargin{1};
end

%% Initialisations

% check if labels are specified
labels = [];
if ~isempty(varargin) && size(varargin{1}, 2) == 1
    labels = varargin{1};
end

% extract the set of labels, without the background
if isempty(labels)
    labels = imFindLabels(img);
end
nLabels = length(labels);

% allocate memory for result
ellipse = zeros(nLabels, 5);

%% Extract ellipse corresponding to each label

for i = 1:nLabels
    % extract points of the current particle
    [y, x] = find(img==labels(i));

    % transform to physical space if needed
    if calib
        x = (x-1) * spacing(1) + origin(1);
        y = (y-1) * spacing(2) + origin(2);
    end

    % compute centroid, used as center of inertia ellipse
    xc = mean(x);
    yc = mean(y);

    % recenter points (should be better for numerical accuracy)
    x = x - xc;
    y = y - yc;

    % number of points
    n = length(x);

    % compute inertia parameters. 1/12 is the contribution of a single
    % pixel, then for regions with only one pixel the resulting ellipse has
    % positive radii.
    Ixx = sum(x.^2) / n + spacing(1)^2/12;
    Iyy = sum(y.^2) / n + spacing(2)^2/12;
    Ixy = sum(x.*y) / n;

    % compute ellipse semi-axis lengths
    common = sqrt( (Ixx - Iyy)^2 + 4 * Ixy^2);
    ra = sqrt(2) * sqrt(Ixx + Iyy + common);
    rb = sqrt(2) * sqrt(Ixx + Iyy - common);

    % compute ellipse angle and convert into degrees
    theta = atan2(2 * Ixy, Ixx - Iyy) / 2;
    theta = theta * 180 / pi;

    % create the resulting inertia ellipse
    ellipse(i,:) = [xc yc ra rb theta];
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [ wav_files ] = find_dat( path )
%FIND_WAV, find all wav file recursively
wav_files = [];
if(isdir(path) == 0)

    return;
end
path_files = dir(path);
fileNum = length(path_files);
for k= 3:fileNum
    file = [path,'\', path_files(k).name];
    if (path_files(k).isdir == 1)
        ret = find_wav(file);
        if(isempty(ret) ~= 1)
            if(isempty(wav_files))
                wav_files = char(ret);
            else
                wav_files = char(wav_files, ret);
            end
        end
    elseif strfind(path_files(k).name, '.dat')
        if(isempty(wav_files))
            wav_files = char(file);
        else
            wav_files = char(wav_files, file);
        end
    end
end
end

author: vipchengrui

repo: traditional-speech-enhancement


text:
function s = report_this_filefun(level)
    % this script replaced all the individual "disp('This is ...') messages
    % it can be used to track which files are used, or do some function/file specific
    % setup
    %  REPORT_THIS_FILEFUN() display message if ANY level has been set
    %  REPORT_THIS_FILEFUN( LEVEL ) print message if

    %called by script or function?
    % modified from: http://blogs.mathworks.com/loren/2013/08/26/what-kind-of-matlab-file-is-this/

    ZG = ZmapGlobal.Data;

    %% change debugging behavior
    if nargin==0
        level=1;
    end

    if ~ZG.debug || level < ZG.debugLevel
        return
    end

    %% ok. Now do the reporting

    dbk=dbstack('-completenames');

    thecaller = 'user interaction';
    thefun = '- none -'; %#ok<NASGU>
    thefilename = '<anon. funct>'; %#ok<NASGU>
    line = nan;

    switch numel(dbk)
        case 1
            msg.dbdisp('nothing to report. [no stack to speak of]')
            return
        case 2
            thefun = dbk(2).name;
            thefilename=dbk(2).file;
        otherwise
            thefun = dbk(2).name;
            thefilename=dbk(2).file;
            thecaller = dbk(3).name;
            line = dbk(3).line;
    end
    % test to see if it is a script

    try
        [~] = nargin(thefun); % errors if not a function
        fn_id_tag = 'function';
    catch ME
        if strcmp(ME.identifier, 'MATLAB:nargin:isScript')
            fn_id_tag = 'script';
        else
            fn_id_tag =  '';
        end
    end
    if nargin > 0
        s=sprintf('%s %s [%s]\n  called by %s line %d\n\n', fn_id_tag, thefun, thefilename, thecaller, line);
    else
        s=sprintf('[%s]\n  called by %s line %d\n', thefilename, thecaller, line);
    end
    msg.dbdisp(s,['in ' fn_id_tag ' [' thefun ']']);
end

author: CelsoReyes


text:
function [Out] = icaDim(Origdata,DEMDT,VN,Iterate,NDist)
% Matthew F. Glasser, Chad Donahue, Steve Smith, Christian Beckmann

%%%%%%%%%%%%%%%%%%%%
%Variables
Out.VNDIM=VN; %Variance Normalization Dimensionality initially set to 1
lnb = 0.5; %Lower noise bound for Golden Section Search
stabThresh = Iterate; %How many iterations meeting dimThresh criterion (1 for VNDIM=1 results, 2 for converged results)

%%%%%%%%%%%%%%%%%%%%

%Remove Constant Timeseries and Detrend Data
if DEMDT==1
    data = Origdata(std(Origdata,[],2)>0,:);
    data_detrend = detrend(data')';
    data_trend = data - data_detrend;
elseif DEMDT==0
    data = Origdata(std(Origdata,[],2)>0,:);
    data_detrend = demean(data')';
    data_trend = single(zeros(size(data,1),size(data,2)));
elseif DEMDT==-1
    data = Origdata(std(Origdata,[],2)>0,:);
    data_detrend = data;
    data_trend = single(zeros(size(data,1),size(data,2)));
end

%Precompute PCA reconstruction
%octave doesn't have "rng", but it does accept "rand('twister', 0);" for matlab compatibility
rand('twister', 0);
[u,EigS,v]=nets_svds(data_detrend',0); %Eigenvalues of detrended data
Out.DOF=sum(diag(EigS)>(std(diag(EigS))*0.1)); %Compute degrees of freedom

c=1;
stabCount = 0;
while stabCount < stabThresh %Loop until dim output is stable

    c
    rand('twister', 0);
    %Variance normalization via PCA reconstruction: Isolate unstructured noise
    if VN~=0
      noise_unst = (u(:,Out.VNDIM(c):Out.DOF)*EigS(Out.VNDIM(c):Out.DOF,Out.VNDIM(c):Out.DOF)*v(:,Out.VNDIM(c):Out.DOF)')';
      Out.noise_unst_std = max(std(noise_unst,[],2),0.001);
      data_detrend_vn = data_detrend ./ repmat(Out.noise_unst_std,1,size(data_detrend,2));
    elseif VN==0
      data_detrend_vn = data_detrend;
      Out.noise_unst_std = single(ones(size(data_detrend,1),1));
    end
    if size(data_detrend_vn,1)<size(data_detrend_vn,2)
        if DEMDT~=-1
            d_pcaNorm=eig(cov(data_detrend_vn'));
        else
            d_pcaNorm=eig(data_detrend_vn*data_detrend_vn');
        end
    else
        if DEMDT~=-1
            d_pcaNorm=eig(cov(data_detrend_vn));
        else
            d_pcaNorm=eig(data_detrend_vn'*data_detrend_vn);
        end
    end
    Out.lambda_pcaNorm=flipud(d_pcaNorm);

    %Enable fitting multiple noise distributions
    DOF=Out.DOF;
    %DOF=sum(Out.lambda_pcaNorm>std(Out.lambda_pcaNorm)*0.1); %Recompute DOF because of demeaning in cov
    MaxX=length(data_detrend_vn);
    if DEMDT==-1
        MaxX=1000000;
    end
    lambda=Out.lambda_pcaNorm;
    S=0;
    clear EN
    for i=1:NDist
        [x(i),en] = FitWishart(lnb,S,DOF,MaxX,lambda(1:DOF));
        lambda=lambda(1:DOF)-en(1:DOF);
        en=[en;single(zeros(Out.DOF-length(en)+1,1))];
        EN(:,i)=en(1:Out.DOF);
        MaxX=x(i);
        DOF=min(find(lambda <= 0));
        lambda=[lambda(1:DOF);single(zeros(Out.DOF-DOF+1,1))];
        lambda=lambda(1:Out.DOF);
    end

    Out.EN=sum(EN,2);
    %Out.x(c)=x(1);
    Out.x(c)=round(x(end));
    %DOF=round(Out.DOF);

    %MaxS=5;
    %DOF=Out.DOF;
    %MaxX=round(Out.x(c));
    %MaxX=length(data_detrend_vn);
    %lambda=Out.lambda_pcaNorm;
    %[Out.x(c),Out.EN,Out.s(c)] = SmoothEst(lnb,MaxS,DOF,MaxX,lambda);

    %[Out.x(c),Out.EN] = FitWishart(lnb,Out.s(c),DOF,MaxX,lambda(1:DOF));

    %Divide eigenvalues by null for adjusted output values
    Out.lambdaAdj = abs(Out.lambda_pcaNorm(1:DOF)./(Out.EN(1:DOF)));
    %Out.lambdaAdj(DOF:end)=1;

    %Normalize adjusted eigenvalues to 1 for input to laplacian
    %Out.lambdaAdj_norm = Out.lambdaAdj ./ max(Out.lambdaAdj);
    Out.pcaDim_lambdaAdj = pca_dim(Out.lambdaAdj',Out.x(c));

    %Maximum of laplacian is estimate of optimal dimensionality
    [~, Out.calcDim]=max(Out.pcaDim_lambdaAdj.lap(1:DOF));

    if VN~=1
      stabCount = 2;
    end

    %Next loop will use calculated dimensionality
    c=c+1;

    Out.VNDIM(c) = Out.calcDim;
    %manual display of Out fields, to make octave less verbose
    disp('Out =');
    disp(['   VNDIM: ' mat2str(Out.VNDIM)]);
    disp(['   DOF: ' mat2str(Out.DOF)]);
    disp(['   noise_unst_std: array of size ' mat2str(size(Out.noise_unst_std))]);
    disp(['   lambda_pcaNorm: array of size ' mat2str(size(Out.lambda_pcaNorm))]);
    disp(['   EN: array of size ' mat2str(size(Out.EN))]);
    disp(['   x: ' mat2str(Out.x)]);
    disp(['   lambdaAdj: array of size ' mat2str(size(Out.lambdaAdj))]);
    disp(['   pcaDim_lambdaAdj: struct']);
    disp(['   calcDim: ' mat2str(Out.calcDim)]);

    %Store dims in array, check number of occurances to prevent dim loops
    stabCount = sum(Out.VNDIM==Out.calcDim);

end %End while loop for dim calcs
rand('twister', 0);
if DEMDT~=-1
    [u,EigS,v]=nets_svds(demean(data_detrend_vn)',0);
else
    [u,EigS,v]=nets_svds(data_detrend_vn',0);
end
u(isnan(u))=0; v(isnan(v))=0;
Out.EigSAdj=single(zeros(length(EigS),1));
Out.grot_one=diag(EigS(1:length(Out.EN),1:length(Out.EN)));
Out.grot_two=Out.grot_one.^2;
Out.grot_three=(Out.grot_two./max(Out.grot_two)).*max(Out.lambda_pcaNorm);
Out.grot_four=Out.grot_three-(Out.EN/median((Out.lambda_pcaNorm(1:Out.DOF)./max(Out.lambda_pcaNorm(1:Out.DOF)))./(diag(EigS(1:Out.DOF,1:Out.DOF).^2)./max(diag(EigS(1:Out.DOF,1:Out.DOF).^2)))));
Out.grot_five=(Out.grot_four./max(Out.lambda_pcaNorm)).*max(diag(EigS.^2));
firstOne = min(find(Out.grot_five <= 0));
Out.grot_six=Out.grot_five;
if ~length(firstOne)==0
    Out.grot_six(firstOne:end) = 0;
    Out.NewDOF=firstOne-1;
else
    Out.NewDOF=length(Out.grot_six);
end
%manual display of Out fields, to make octave less verbose
disp('Out =');
disp(['   VNDIM: ' mat2str(Out.VNDIM)]);
disp(['   DOF: ' mat2str(Out.DOF)]);
disp(['   noise_unst_std: array of size ' mat2str(size(Out.noise_unst_std))]);
disp(['   lambda_pcaNorm: array of size ' mat2str(size(Out.lambda_pcaNorm))]);
disp(['   EN: array of size ' mat2str(size(Out.EN))]);
disp(['   x: ' mat2str(Out.x)]);
disp(['   lambdaAdj: array of size ' mat2str(size(Out.lambdaAdj))]);
disp(['   pcaDim_lambdaAdj: struct']);
disp(['   calcDim: ' mat2str(Out.calcDim)]);
disp(['   EigSAdj: array of size ' mat2str(size(Out.EigSAdj))]);
disp(['   grot_one: array of size ' mat2str(size(Out.grot_one))]);
disp(['   grot_two: array of size ' mat2str(size(Out.grot_two))]);
disp(['   grot_three: array of size ' mat2str(size(Out.grot_three))]);
disp(['   grot_four: array of size ' mat2str(size(Out.grot_four))]);
disp(['   grot_five: array of size ' mat2str(size(Out.grot_five))]);
disp(['   grot_six: array of size ' mat2str(size(Out.grot_six))]);
disp(['   NewDOF: ' mat2str(Out.NewDOF)]);

Out.EigSAdj(1:length(Out.EN))=sqrt(Out.grot_six);

Out.data=single(zeros(size(Origdata,1),size(Origdata,2)));
if DEMDT~=-1
    Out.data(std(Origdata,[],2)>0,:)=data_trend + (((u*diag(Out.EigSAdj)*v')'+repmat(mean(data_detrend_vn),size(data_detrend_vn,1),1)) .* repmat(Out.noise_unst_std,1,size(data_detrend_vn,2)));
else
    Out.data(std(Origdata,[],2)>0,:)=data_trend + (((u*diag(Out.EigSAdj)*v')') .* repmat(Out.noise_unst_std,1,size(data_detrend_vn,2)));
end
temp=single(zeros(size(Origdata,1),1)); temp(std(Origdata,[],2)>0,:)=Out.noise_unst_std; Out.noise_unst_std=max(temp,0.001); clear temp;

end %End function

function [out] = lpdist(in)
    out=pdist(log(in));
end

function [x,EN] = FitWishart(lnb,S,DOF,MaxX,lambda) %FitWishart(lnb,step,DOF,MaxX,lambda)
    rand('twister', 0);
    EigDn1=round(DOF*lnb); %Isolate search to noise
    EigDn2=round(DOF-1); %Reqd for post MR+FIX deconcatinated tcs
    %EigDn2=round(DOF*0.75); %Reqd for post MR+FIX deconcatinated tcs

    a = DOF; %Lower bound for search range
    b = MaxX; %Upper bound for search range
    epsilon = 1; %Accuracy/stopping criterion
    iter = 500; %# iterations/secondary stopping criterion
    tau = double((sqrt(5)-1)/2); %Golden ratio (constant), 0.618...
    k = 1; %Iteration count

    %Initial section ranges to instantiate optimization
    x1 = a+(1-tau)*(b-a);
    x2 = a+tau*(b-a);
    %x1=a;
    %x2=b;

    %Calculate initial null spectra
    EN_x1 = iFeta([0:0.001:5],DOF,x1)'; %Call feta to calc null spectrum
    %EN_x1=flipud(eig(cov(Smooth(randn(round(x1),DOF),S))));
    EN_x1=EN_x1*median(lambda(EigDn1:EigDn2)./EN_x1(EigDn1:EigDn2)); %Remove offset between null & data
    f_x1 = lpdist([EN_x1(EigDn1:EigDn2)'; lambda(EigDn1:EigDn2)']); %Compute pairwise distance b/w null & data

    EN_x2 = iFeta([0:0.001:5],DOF,x2)';
    %EN_x2=flipud(eig(cov(Smooth(randn(round(x2),DOF),S))));
    EN_x2=EN_x2*median(lambda(EigDn1:EigDn2)./EN_x2(EigDn1:EigDn2));
    f_x2 = lpdist([EN_x2(EigDn1:EigDn2)'; lambda(EigDn1:EigDn2)']);

    while (abs(b-a)>epsilon) && (k<iter) %Loop until low error OR max iter met
        %k=k+1;
        rand('twister', 0);
        disp([num2str(a) ' ' num2str(x1) ' ' num2str(x2) ' ' num2str(b)]);
        %Check both terms for minimal pairwise distance and continue toward minimum
        if (f_x1<f_x2)
            b=x2;
            x2=x1;
            f_x2=f_x1;%don't recompute, this is the entire point of golden search
            EN_x2=EN_x1;
            x1=a+(1-tau)*(b-a);

            EN_x1 = iFeta([0:0.001:5],DOF,x1)';
            %EN_x1=flipud(eig(cov(Smooth(randn(round(x1),DOF),S))));
            EN_x1=EN_x1*median(lambda(EigDn1:EigDn2)./EN_x1(EigDn1:EigDn2));
            f_x1 = lpdist([EN_x1(EigDn1:EigDn2)'; lambda(EigDn1:EigDn2)']);

        else
            a=x1;
            x1=x2;
            f_x1=f_x2;
            EN_x1=EN_x2;
            x2=a+tau*(b-a);

            EN_x2 = iFeta([0:0.001:5],DOF,x2)';
            %EN_x2=flipud(eig(cov(Smooth(randn(round(x2),DOF),S))));
            EN_x2=EN_x2*median(lambda(EigDn1:EigDn2)./EN_x2(EigDn1:EigDn2));
            f_x2 = lpdist([EN_x2(EigDn1:EigDn2)'; lambda(EigDn1:EigDn2)']);
        end
        k = k+1;
    end

    %Compute null spectrum w/ optimal pairwise distance for output
    if (f_x1<f_x2)
        x = x1;
    else
        x = x2;
    end
    rand('twister', 0);
    %disp(randn(1,1));
    EN=flipud(eig(cov(Smooth(randn(round(x),DOF),S))));

    x=x/1+S;

    %Remove offset from null spectrum
    EN = EN*median(lambda(EigDn1:EigDn2) ./ EN(EigDn1:EigDn2));
end

function [x,EN,s] = SmoothEst(lnb,MaxS,DOF,MaxX,lambda)

    EigDn1=round(DOF*lnb); %Isolate search to noise
    EigDn2=round(DOF-1);

    M=randn(MaxX,DOF);

    a = 0; %Lower bound for search range
    b = MaxS; %Upper bound for search range
    epsilon = 0.01; %Accuracy/stopping criterion
    iter = 500; %# iterations/secondary stopping criterion
    tau = double((sqrt(5)-1)/2); %Golden ratio (constant), 0.618...
    k = 1; %Iteration count

    %Initial section ranges to instantiate optimization
    s1 = a+(1-tau)*(b-a);
    s2 = a+tau*(b-a);
    %x1=a;
    %x2=b;

    %Calculate initial null spectra
    %EN_x1 = iFeta([0:step:5],DOF,x1)'; %Call feta to calc null spectrum
    EN_s1=flipud(eig(cov(Smooth(M,s1/(2*sqrt(2*log(2)))))));
    EN_s1=EN_s1*median(lambda(EigDn1:EigDn2)./EN_s1(EigDn1:EigDn2)); %Remove offset between null & data
    f_s1 = lpdist([EN_s1(EigDn1:EigDn2)'; lambda(EigDn1:EigDn2)']); %Compute pairwise distance b/w null & data

    %EN_x2 = iFeta([0:step:5],DOF,x2)';
    EN_s2=flipud(eig(cov(Smooth(M,s2/(2*sqrt(2*log(2)))))));
    EN_s2=EN_s2*median(lambda(EigDn1:EigDn2)./EN_s2(EigDn1:EigDn2));
    f_s2 = lpdist([EN_s2(EigDn1:EigDn2)'; lambda(EigDn1:EigDn2)']);

    while (abs(b-a)>epsilon) && (k<iter) %Loop until low error OR max iter met
        %k=k+1;
        disp([num2str(a) ' ' num2str(s1) ' ' num2str(s2) ' ' num2str(b)]);
        %Check both terms for minimal pairwise distance and continue toward minimum
        if (f_s1<f_s2)
            b=s2;
            s2=s1;
            f_s2=f_s1;%don't recompute, this is the entire point of golden search
            EN_s2=EN_s1;
            s1=a+(1-tau)*(b-a);

            %EN_x1 = iFeta([0:step:5],DOF,x1)';
            EN_s1=flipud(eig(cov(Smooth(M,s1/(2*sqrt(2*log(2)))))));
            EN_s1=EN_s1*median(lambda(EigDn1:EigDn2)./EN_s1(EigDn1:EigDn2));
            f_s1 = lpdist([EN_s1(EigDn1:EigDn2)'; lambda(EigDn1:EigDn2)']);

        else
            a=s1;
            s1=s2;
            f_s1=f_s2;
            EN_s1=EN_s2;
            s2=a+tau*(b-a);

            %EN_x2 = iFeta([0:step:5],DOF,x2)';
            EN_s2=flipud(eig(cov(Smooth(M,s2))));
            EN_s2=EN_s2*median(lambda(EigDn1:EigDn2)./EN_s2(EigDn1:EigDn2));
            f_s2 = lpdist([EN_s2(EigDn1:EigDn2)'; lambda(EigDn1:EigDn2)']);
        end
        k = k+1;
    end

    %Compute null spectrum w/ optimal pairwise distance for output
    if (f_s1<f_s2)
        s = s1;
    else
        s = s2;
    end

    %EN = iFeta([0:step:5],DOF,x)';
    %EN=flipud(eig(cov(randn(round(x),DOF))));
    EN=flipud(eig(cov(Smooth(M,s/(2*sqrt(2*log(2)))))));
    x=MaxX/(1+s); %Number of resels
    %Remove offset from null spectrum
    EN = EN*median(lambda(EigDn1:EigDn2) ./ EN(EigDn1:EigDn2));
end

function [O] = Smooth(M,S)
if S==0
    O=M;
else
sigma = S/(2*sqrt(2*log(2)));
%sz = S*6;    % length of gaussFilter vector
%x = linspace(-sz / 2, sz / 2, sz);
%gaussFilter = exp(-x .^ 2 / (2 * sigma ^ 2));
%gaussFilter = gaussFilter / sum (gaussFilter); % normalize

width = round((6*sigma - 1)/2);
support = (-width:width);
gaussFilter = exp( -(support).^2 ./ (2*sigma^2) );
gaussFilter = gaussFilter/ sum(gaussFilter);

O=single(zeros(size(M,1),size(M,2)));

for i=1:size(M,2)
    op
    O(:,i) = conv (M(:,i), gaussFilter, 'same');
    %O(:,i) = filter (gaussFilter,1,M(:,i));
end
end
end

author: yetianmed

repo: subcortex

save_path: github-repos/MATLAB/yetianmed-subcortex


text:
%DELTA2TR Convert differential motion  to SE(3) homogeneous transform
%
% T = DELTA2TR(D) is a homogeneous transform (4x4) representing differential
% motion D (6x1).
%
% The vector D=(dx, dy, dz, dRx, dRy, dRz) represents infinitessimal translation
% and rotation, and is an approximation to the instantaneous spatial velocity
% multiplied by time step.
%
% Reference::
% - Robotics, Vision & Control: Second Edition, P. Corke, Springer 2016; p67.
%
% See also tr2delta, SE3.delta.

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math

function delta = delta2tr(d)
    d = d(:);
    delta = eye(4,4) + [skew(d(4:6)) d(1:3); 0 0 0 0];

author: petercorke

repo: spatialmath-matlab


text:
function [mach] = kmph2mach(kmph)
% Convert speed from miles per hour to mach number (at STP!)
% Chad A. Greene 2012
mach = kmph*0.0008098477486233;

author: Sable

repo: mcbench-benchmarks


text:
close all;
clear all;
clc;

data_file_path = 'bin/ra_mmv_phase_transition_snr_30db_s_4.mat';
options.export = true;
options.export_dir = 'bin';
options.export_name = 'ra_mmv_snr_30_db_s_4';
options.chosen_ks = [2, 4, 8, 16, 32, 64];
options.subtitle = 'Rank Aware, SNR=30dB s=4';
spx.pursuit.PhaseTransitionAnalysis.print_results(data_file_path, ...
    'CoSaMP', options);

author: indigits

repo: sparse-plex


text:
function A = VBA_spm_mesh_adjacency(F)
% Compute the adjacency matrix of a triangle mesh
% FORMAT A = spm_mesh_adjacency(F)
% F        - a [fx3] faces array or a patch structure
%
% A        - adjacency matrix as a sparse [vxv] array
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_mesh_adjacency.m 4035 2010-08-05 18:54:32Z guillaume $

if ~isnumeric(F) && isfield(F,'vertices')
    N = size(F.vertices,1);
    F = double(F.faces);
    A = sparse([F(:,1); F(:,1); F(:,2); F(:,2); F(:,3); F(:,3)], ...
           [F(:,2); F(:,3); F(:,1); F(:,3); F(:,1); F(:,2)], 1, N, N);
else
    if isstruct(F), F = F.faces; end
    F = double(F);
    A = sparse([F(:,1); F(:,1); F(:,2); F(:,2); F(:,3); F(:,3)], ...
           [F(:,2); F(:,3); F(:,1); F(:,3); F(:,1); F(:,2)], 1);
end

A = double(A > 0);

author: MBB-team

repo: VBA-toolbox


text:
% ==============================================================================
%
%    Software License Agreement (BSD License)
%    Copyright (c) 2019
%    (www.aimlab.wpi.edu)
%
%    All rights reserved.
%
%    Redistribution and use in source and binary forms, with or without
%    modification, are permitted provided that the following conditions
%    are met:
%
%    * Redistributions of source code must retain the above copyright
%    notice, this list of conditions and the following disclaimer.
%
%    * Redistributions in binary form must reproduce the above
%    copyright notice, this list of conditions and the following
%    disclaimer in the documentation and/or other materials provided
%    with the distribution.
%
%    * Neither the name of authors nor the names of its contributors may
%    be used to endorse or promote products derived from this software
%    without specific prior written permission.
%
%    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
%    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
%    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
%    FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
%    COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
%    INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
%    BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
%    LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
%    CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
%    LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
%    ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%    POSSIBILITY OF SUCH DAMAGE.
%
%    \author:    <http://www.aimlab.wpi.edu>
%    \author:    <amunawar@wpi.edu>
%    \author:    Adnan Munawar
%    \version:   0.1$
% ==============================================================================

function benchmarkRRT

clc;
close all;
clear all;

num_of_runs =1;
run_RRTconnect =0;
run_RRTextend = 0;
run_LazyRRT = 0;
run_RRTstar = 1;

dim = 2;
stepsize = [5];

random_world = 0;
radius = 10;
samples = 4000;

show_output = 1;
show_benchmark_results = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t_lazy = [];
t_extend = [];
t_connect = [];
t_star = [];

l_lazy = [];
l_extend = [];
l_connect = [];
l_star = [];

p_lazy = [];
p_extend = [];
p_connect = [];
p_star = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for sits = 1:size(stepsize,2)
    segmentLength = stepsize(sits);
    if run_LazyRRT == 1
        time = 0;
        avg_its = 0;
        avg_path = 0;

        for i = 1:num_of_runs

            [n_its path_n run_time] =  LazyRRT3D(dim,segmentLength,random_world,show_output);
            time = time + run_time;
            avg_its = avg_its + n_its;
            avg_path = avg_path + path_n;
        end

        str1 = ['The time taken by Lazy RRT for ', num2str(num_of_runs), ' runs is ', num2str(time)];
        str2 = ['The averagae time taken by Lazy RRT for each run is ', num2str(time/num_of_runs)];
        str3 = ['The averagae number of states explored by Lazy RRT for each run is ', num2str(avg_its/num_of_runs)];
        str4 = ['The averagae number of state in Path by Lazy RRT for each run is ', num2str(avg_path/num_of_runs)];

        disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
        disp(str1);
        disp(str2);
        disp(str3);
        disp(str4);
        disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');

        t_lazy = [t_lazy time/num_of_runs];
        l_lazy = [l_lazy avg_its/num_of_runs];
        p_lazy = [p_lazy avg_path/num_of_runs];

    end

    if run_RRTstar == 1
        time = 0;
        avg_its = 0;
        avg_path = 0;
        for i = 1:num_of_runs
            [n_its path_n,run_time] = RRTstar3D(dim,segmentLength,radius,random_world,show_output,samples);
            time = time + run_time;
            avg_its = avg_its + n_its;
            avg_path = avg_path + path_n;
        end

        str1 = ['The time taken by RRT-Star for ', num2str(num_of_runs), ' runs is ', num2str(time)];
        str2 = ['The averagae time taken by RRT_Star for each run is ', num2str(time/num_of_runs)];
        str3 = ['The averagae number of states explored by RRT_Star for each run is ', num2str(avg_its/num_of_runs)];
        str4 = ['The averagae number of state in Path by RRT-Star for each run is ', num2str(avg_path/num_of_runs)];

        disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
        disp(str1);
        disp(str2);
        disp(str3);
        disp(str4);
        disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');

        t_star = [t_star time/num_of_runs];
        l_star = [l_star avg_its/num_of_runs];
        p_star = [p_star avg_path/num_of_runs];

    end

    if run_RRTextend == 1
        time = 0;
        avg_its = 0;
        avg_path = 0;
        for i = 1:num_of_runs
            [n_its path_n,run_time] = RRTextend3D(dim,segmentLength,random_world,show_output);
            time = time + run_time;
            avg_its = avg_its + n_its;
            avg_path = avg_path + path_n;
        end

        str1 = ['The time taken by RRT-Extend for ', num2str(num_of_runs), ' runs is ', num2str(time)];
        str2 = ['The averagae time taken by RRT_Extend for each run is ', num2str(time/num_of_runs)];
        str3 = ['The averagae number of states explored by RRT_Extend for each run is ', num2str(avg_its/num_of_runs)];
        str4 = ['The averagae number of state in Path by RRT-Extend for each run is ', num2str(avg_path/num_of_runs)];

        disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
        disp(str1);
        disp(str2);
        disp(str3);
        disp(str4);
        disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');

        t_extend = [t_extend time/num_of_runs];
        l_extend = [l_extend avg_its/num_of_runs];
        p_extend = [p_extend avg_path/num_of_runs];

    end

    if run_RRTconnect == 1

        time = 0;
        avg_its = 0;
        avg_path = 0;

        for i = 1:num_of_runs
            [n_its path_n,run_time] =  RRTconnect3D(dim,segmentLength,random_world,show_output);
            time = time + run_time;
            avg_its = avg_its + n_its;
            avg_path = avg_path + path_n;
        end

        str1 = ['The time taken by RRT-Connect for ', num2str(num_of_runs), ' runs is ', num2str(time)];
        str2 = ['The averagae time taken by RRT-Connect for each run is ', num2str(time/num_of_runs)];
        str3 = ['The averagae number of states explored by RRT-Connect for each run is ', num2str(avg_its/num_of_runs)];
        str4 = ['The averagae number of state in Path by RRT-Connect for each run is ', num2str(avg_path/num_of_runs)];

        disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
        disp(str1);
        disp(str2);
        disp(str3);
        disp(str4);
        disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');

        t_connect = [t_connect time/num_of_runs];
        l_connect = [l_connect avg_its/num_of_runs];
        p_connect = [p_connect avg_path/num_of_runs];
    end
end

if show_benchmark_results == 1

    figure;
    hold on;
    plot(stepsize,t_lazy,'r','LineWidth',2);
    plot(stepsize,t_extend,'g','LineWidth',2);
    plot(stepsize,t_connect,'b','LineWidth',2);
    ylabel('Computational Time');
    xlabel('Step Size');
    dim_str = sprintf('Comparison of computational time for %d Dimensional C-Space',dim);
    title(dim_str)
    hold off;

    figure;
    hold on;
    plot(stepsize,l_lazy,'r','LineWidth',2);
    plot(stepsize,l_extend,'g','LineWidth',2);
    plot(stepsize,l_connect,'b','LineWidth',2);
    ylabel('Number of States Explored');
    xlabel('Step Size');
    dim_str = sprintf(' Comparison of number of states explored for %d Dimensional C-Space',dim);
    title(dim_str)
    hold off;

    figure;
    hold on;
    plot(stepsize,p_lazy,'r','LineWidth',2);
    plot(stepsize,p_extend,'g','LineWidth',2);
    plot(stepsize,p_connect,'b','LineWidth',2);
    ylabel('Number of States in Path');
    xlabel('Step Size');
    dim_str = sprintf('Comparison for number of states in path for %d Dimensional C-Space',dim);
    title(dim_str)
    hold off;

end

%     t_lazy
%     l_lazy
%     p_lazy
%
%
%
%     t_extend
%     l_extend
%     p_extend
%
%     t_connect
%     l_connect
%     p_connect

end

author: adnanmunawar

repo: matlab-rrt-variants


text:
function Mi = get_block_row(M, C, row_range)
	%% ================== File info ==========================
	% Author		: Tiep Vu (http://www.personal.psu.edu/thv102/)
	% Time created	: 1/27/2016 2:44:26 AM
	% Last modified	: 1/27/2016 2:44:30 AM
	% Description	: Get a row block of a big matrix X.
	% 	INPUT
	%		M        : the big matrix. M = [M1 ; M2 ; ... ; MC]
	%		C        : block index
	%		row_range: a vector store the last index of each block. row_range(1) = 0.
	%					i-th block is indexed by col_range(i)+1: col_range(i+1).
	% 	OUTPUT
	%		Mi: output block matrix
	%
	%% ================== end File info ==========================
	id_sel = [];
	for i = 1: numel(C)
		c = C(i);
		id_sel = [id_sel, row_range(c) + 1: row_range(c+1)];
	end
	Mi = M(id_sel, : , :);
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function x = daub2_transform_inverse ( n, y )

%*****************************************************************************80
%
%% DAUB2_TRANSFORM_INVERSE inverts the DAUB2 transform of a vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the vector.
%    N must be a power of 2.
%
%    Input, real Y(N), the transformed vector.
%
%    Output, real X(N), the original vector.
%
  c = [ 7.071067811865475E-01, ...
        7.071067811865475E-01 ];

  x(1:n,1) = y(1:n);
  z(1:n,1) = 0.0;

  m = 1;

  while ( m * 2 <= n )

    for i = 1 : m
      z(2*i-1,1) = c(1) * ( x(i) + x(i+m) );
      z(2*i,1)   = c(2) * ( x(i) - x(i+m) );
    end

    x(1:2*m,1) = z(1:2*m);

    m = m * 2;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ yval, ypval ] = spline_quadratic_val ( ndata, tdata, ydata, tval )

%*****************************************************************************80
%
%% SPLINE_QUADRATIC_VAL evaluates a piecewise quadratic spline at a point.
%
%  Discussion:
%
%    Because of the simple form of a piecewise quadratic spline,
%    the raw data points ( TDATA(I), YDATA(I)) can be used directly to
%    evaluate the spline at any point.  No processing of the data
%    is required.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, int NDATA, the number of data points defining the spline.
%    NDATA should be odd.
%
%    Input, real TDATA(NDATA), YDATA(NDATA), the values of the independent
%    and dependent variables at the data points.  The values of TDATA should
%    be distinct and increasing.
%
%    Input, real TVAL, the point at which the spline is to be evaluated.
%
%    Output, real YVAL, YPVAL, the value of the spline and its first
%    derivative dYdT at TVAL.  YPVAL is not reliable if TVAL is exactly
%    equal to TDATA(I) for some I.
%
  if ( mod ( ndata, 2 ) == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SPLINE_QUADRATIC_VAL - Fatal error!\n' );
    fprintf ( 1, '  NDATA must be odd.\n' );
    error ( 'SPLINE_QUADRATIC_VAL - Fatal error!' );
  end
%
%  Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] that contains, or is
%  nearest to, TVAL.
%
  [ left, right ] = r8vec_bracket ( ndata, tdata, tval );
%
%  Force LEFT to be odd.
%
  if ( mod ( left, 2 ) == 0 )
    left = left - 1;
  end
%
%  Copy out the three abscissas.
%
  t1 = tdata(left);
  t2 = tdata(left+1);
  t3 = tdata(left+2);

  if ( t2 <= t1 | t3 <= t2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SPLINE_QUADRATIC_VAL - Fatal error!\n' );
    fprintf ( 1, '  T2 <= T1 or T3 <= T2.\n' );
    error ( 'SPLINE_QUADRATIC_VAL - Fatal error!' );
  end
%
%  Construct and evaluate a parabolic interpolant for the data
%  in each dimension.
%
  y1 = ydata(left);
  y2 = ydata(left+1);
  y3 = ydata(left+2);

  dif1 = ( y2 - y1 ) / ( t2 - t1 );

  dif2 = ( ( y3 - y1 ) / ( t3 - t1 ) ...
       - ( y2 - y1 ) / ( t2 - t1 ) ) / ( t3 - t2 );

  yval = y1 + ( tval - t1 ) * ( dif1 + ( tval - t2 ) * dif2 );
  ypval = dif1 + dif2 * ( 2.0E+00 * tval - t1 - t2 );

  return
end


text:
function sys = double(X)
%double           Overloaded.

X = flatten(X);
nlmi = length(X.LMIid);

if (nlmi == 0)
 sys = NaN;
end

if nlmi>1
 error('Double not applicable on list of constraints')
end

sys = double(X.clauses{1}.data);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [Y X Z n m p T k1 k3 q1 q2 q3]=maprelim(data_endo,data_exo,const,lags,regimeperiods,names)

% function [Y X Z n m p T k1 k3 q1 q2 q3]=maprelim(data_endo,data_exo,const,lags)
% creates the preliminary values necessary for all the subsequent computations
% inputs:  - matrix 'data_endo': the matrix storing the endogenous time series data used to estimate the model
%          - matrix 'data_exo': the matrix storing the exogenous time series data used to estimate the model
%          - integer 'const': 0-1 value determining whether a constant term should be included in the model
%          - integer 'lags': the number of lags to include in the model
% outputs: - matrix 'Y': the matrix of endogenous variables, defined in (3.5.10)
%          - matrix 'X': the matrix of endogenous regressors, defined in (3.5.10)
%          - matrix 'Z': the matrix of exogenous regressors, defined in (3.5.10)
%          - integer 'n': the number of endogenous variables in the model
%          - integer 'm': the number of exogenous variables in the model
%          - integer 'p': the number of lags in the model
%          - integer 'T': the sample size, i.e. the number of time periods used to estimate the model
%          - integer 'k1': the number of coefficients related to the endogenous variables for each equation in the model
%          - integer 'k3': the number of coefficients related to the exogenous variables for each equation, in the reformulated model (3.5.5)
%          - integer 'q1': the total number of VAR coefficients related to the endogenous variables
%          - integer 'q2': the total number of VAR coefficients related to the exogenous variables
%          - integer 'q3': the total number of VAR coefficients related to the exogenous variables, in the reformulated model (3.5.5)

% first compute p, the number of lags in the model, defined p77
p=lags;

% then compute n, the number of endogenous variables in the model; it is simply the number of columns in the matrix 'data_endo'
n=size(data_endo,2);

% augment the matrix of exogenous with a column of ones to account for the constant
%data_exo=[ones(size(data_endo,1),1) data_exo];
%alternative data exo

if isempty(regimeperiods);
    data_exo=[ones(size(data_endo,1),1)];
else
data_exo=[ones(size(data_endo,1),1) zeros(size(data_endo,1),1)];
data_exo(find(strcmp(names(2:end,1),regimeperiods(1))):find(strcmp(names(2:end,1),regimeperiods(2))),2)=1;
data_exo(find(strcmp(names(2:end,1),regimeperiods(1))):find(strcmp(names(2:end,1),regimeperiods(2))),1)=0;
end

if isempty(regimeperiods);
    data_exo=[ones(size(data_endo,1),1)];
else
data_exo=[ones(size(data_endo,1),1) zeros(size(data_endo,1),1)];
data_exo(find(strcmp(names(2:end,1),regimeperiods(1))):find(strcmp(names(2:end,1),regimeperiods(2))),2)=1;
end
% then compute m, the number of exogenous variables in the model, defined p77
% if data_exo is empty, set m=0
m=size(data_exo,2);

% estimate k1, the number of parameters related to endogenous variables in each equation, defined p77
k1=n*p;

% estimate q1, the total number of parameters related to endogenous variables, defined p77
q1=n*k1;

% estimate q2, the total number of parameters related to exogenous variables, defined p77
q2=n*m;

% estimate k3, the number of parameters related to exogenous variables in each equation of the modified system, defined p78
k3=m*(p+1);

% estimate q3, the total number of parameters related to endogenous variables in the modified system, defined p78
q3=n*k3;

% obtain the matrices Y and X, defined in (3.6.10)
% to do so, use the lagx function on the data matrix
temp=bear.lagx(data_endo,lags);

% to build X, take off the n initial columns of current data
X=temp(:,n+1:end);

% save the n first columns of temp as Y
Y=temp(:,1:n);

% obtain the matrix Z, defined in (3.6.10)
temp=bear.lagx(data_exo,lags);
temp(:,m+1:end)=-temp(:,m+1:end);
Z=temp;

% Define T, the number of periods of the model, as the number of rows of Y
T=size(Y,1);

% obtain eventually y, Xbar and Zbar, as defined in (XXX)
y=Y(:);
Xbar=kron(eye(n),X);
Zbar=kron(eye(n),Z);

author: european-central-bank

repo: BEAR-toolbox


text:
%% clear the workspace and select data
clear; clc; close all;

%% choose data
neuron = Sources2D();
nam = get_fullname('./data_1p.tif');          % this demo data is very small, here we just use it as an example
nam = neuron.select_data(nam);  %if nam is [], then select data interactively

%% parameters
% -------------------------    COMPUTATION    -------------------------  %
pars_envs = struct('memory_size_to_use', 8, ...   % GB, memory space you allow to use in MATLAB
    'memory_size_per_patch', 0.6, ...   % GB, space for loading data within one patch
    'patch_dims', [64, 64]);  %GB, patch size

% -------------------------      SPATIAL      -------------------------  %
gSig = 3;           % pixel, gaussian width of a gaussian kernel for filtering the data. 0 means no filtering
gSiz = 13;          % pixel, neuron diameter
ssub = 1;           % spatial downsampling factor
with_dendrites = true;   % with dendrites or not
if with_dendrites
    % determine the search locations by dilating the current neuron shapes
    updateA_search_method = 'dilate';  %#ok<UNRCH>
    updateA_bSiz = 5;
    updateA_dist = neuron.options.dist;
else
    % determine the search locations by selecting a round area
    updateA_search_method = 'ellipse'; %#ok<UNRCH>
    updateA_dist = 5;
    updateA_bSiz = neuron.options.dist;
end
spatial_constraints = struct('connected', true, 'circular', false);  % you can include following constraints: 'circular'
spatial_algorithm = 'hals_thresh';

% -------------------------      TEMPORAL     -------------------------  %
Fs = 10;             % frame rate
tsub = 1;           % temporal downsampling factor
deconv_options = struct('type', 'ar1', ... % model of the calcium traces. {'ar1', 'ar2'}
    'method', 'foopsi', ... % method for running deconvolution {'foopsi', 'constrained', 'thresholded'}
    'smin', -5, ...         % minimum spike size. When the value is negative, the actual threshold is abs(smin)*noise level
    'optimize_pars', true, ...  % optimize AR coefficients
    'optimize_b', true, ...% optimize the baseline);
    'max_tau', 100);    % maximum decay time (unit: frame);

nk = 3;             % detrending the slow fluctuation. usually 1 is fine (no detrending)
% when changed, try some integers smaller than total_frame/(Fs*30)
detrend_method = 'spline';  % compute the local minimum as an estimation of trend.

% -------------------------     BACKGROUND    -------------------------  %
bg_model = 'ring';  % model of the background {'ring', 'svd'(default), 'nmf'}
nb = 1;             % number of background sources for each patch (only be used in SVD and NMF model)
ring_radius = 18;  % when the ring model used, it is the radius of the ring used in the background model.
%otherwise, it's just the width of the overlapping area
num_neighbors = []; % number of neighbors for each neuron
bg_ssub = 2;        % downsample background for a faster speed

% -------------------------      MERGING      -------------------------  %
show_merge = false;  % if true, manually verify the merging step
merge_thr = 0.65;     % thresholds for merging neurons; [spatial overlap ratio, temporal correlation of calcium traces, spike correlation]
method_dist = 'max';   % method for computing neuron distances {'mean', 'max'}
dmin = 5;       % minimum distances between two neurons. it is used together with merge_thr
dmin_only = 2;  % merge neurons if their distances are smaller than dmin_only.
merge_thr_spatial = [0.8, 0.4, -inf];  % merge components with highly correlated spatial shapes (corr=0.8) and small temporal correlations (corr=0.1)

% -------------------------  INITIALIZATION   -------------------------  %
K = [];             % maximum number of neurons per patch. when K=[], take as many as possible.
min_corr = 0.8;     % minimum local correlation for a seeding pixel
min_pnr = 8;       % minimum peak-to-noise ratio for a seeding pixel
min_pixel = gSig^2;      % minimum number of nonzero pixels for each neuron
bd = 0;             % number of rows/columns to be ignored in the boundary (mainly for motion corrected data)
frame_range = [];   % when [], uses all frames
save_initialization = false;    % save the initialization procedure as a video.
use_parallel = true;    % use parallel computation for parallel computing
show_init = true;   % show initialization results
choose_params = true; % manually choose parameters
center_psf = true;  % set the value as true when the background fluctuation is large (usually 1p data)
% set the value as false when the background fluctuation is small (2p)

% -------------------------  Residual   -------------------------  %
min_corr_res = 0.7;
min_pnr_res = 6;
seed_method_res = 'auto';  % method for initializing neurons from the residual
update_sn = true;

% ----------------------  WITH MANUAL INTERVENTION  --------------------  %
with_manual_intervention = true;

% -------------------------  FINAL RESULTS   -------------------------  %
save_demixed = true;    % save the demixed file or not
kt = 3;                 % frame intervals

% -------------------------    UPDATE ALL    -------------------------  %
neuron.updateParams('gSig', gSig, ...       % -------- spatial --------
    'gSiz', gSiz, ...
    'ring_radius', ring_radius, ...
    'ssub', ssub, ...
    'search_method', updateA_search_method, ...
    'bSiz', updateA_bSiz, ...
    'dist', updateA_bSiz, ...
    'spatial_constraints', spatial_constraints, ...
    'spatial_algorithm', spatial_algorithm, ...
    'tsub', tsub, ...                       % -------- temporal --------
    'deconv_options', deconv_options, ...
    'nk', nk, ...
    'detrend_method', detrend_method, ...
    'background_model', bg_model, ...       % -------- background --------
    'nb', nb, ...
    'ring_radius', ring_radius, ...
    'num_neighbors', num_neighbors, ...
    'bg_ssub', bg_ssub, ...
    'merge_thr', merge_thr, ...             % -------- merging ---------
    'dmin', dmin, ...
    'method_dist', method_dist, ...
    'min_corr', min_corr, ...               % ----- initialization -----
    'min_pnr', min_pnr, ...
    'min_pixel', min_pixel, ...
    'bd', bd, ...
    'center_psf', center_psf);
neuron.Fs = Fs;

%% distribute data and be ready to run source extraction
neuron.getReady(pars_envs);

%% initialize neurons from the video data within a selected temporal range
if choose_params
    % change parameters for optimized initialization
    [gSig, gSiz, ring_radius, min_corr, min_pnr] = neuron.set_parameters();
end

[center, Cn, PNR] = neuron.initComponents_parallel(K, frame_range, save_initialization, use_parallel);
neuron.compactSpatial();
if show_init
    figure();
    ax_init= axes();
    imagesc(Cn, [0, 1]); colormap gray;
    hold on;
    plot(center(:, 2), center(:, 1), '.r', 'markersize', 10);
end

%% estimate the background components
neuron.update_background_parallel(use_parallel);

%% udpate spatial&temporal components, delete false positives and merge neurons
for m=1:2
    % update spatial
    neuron.update_spatial_parallel(use_parallel);

    % update temporal
    neuron.update_temporal_parallel(use_parallel);
end

%% post-process the results automatically
% delete bad neurons
neuron.remove_false_positives();

% merge neurons based on temporal correlation + distances
neuron.merge_neurons_dist_corr(show_merge);

% merge neurons with high correlations
neuron.merge_high_corr(show_merge, merge_thr_spatial);

%% pick neurons from the residual
[center_res, Cn_res, PNR_res] =neuron.initComponents_residual_parallel([],...
    save_initialization, use_parallel, min_corr_res, min_pnr_res, seed_method_res);
% if show_init
%     axes(ax_init);
%     plot(center_res(:, 2), center_res(:, 1), '.g', 'markersize', 10);
% end

% estimate the background components
neuron.update_background_parallel(use_parallel);

% udpate spatial&temporal components, delete false positives and merge neurons
for m=1:2
    % update spatial
    neuron.update_spatial_parallel(use_parallel);

    % update temporal
    neuron.update_temporal_parallel(use_parallel);
end
% delete bad neurons
neuron.remove_false_positives();

% merge neurons based on temporal correlation + distances
neuron.merge_neurons_dist_corr(show_merge);

% merge neurons with high correlations
neuron.merge_high_corr(show_merge, merge_thr_spatial);

%% add a manual intervention and run the whole procedure for a second time
if with_manual_intervention
    neuron.orderROIs('snr');

    % interactively merge neurons
    show_merge = true;
    neuron.merge_close_neighbors(true, dmin_only);

    % interactively delete neurons
    tags = neuron.tag_neurons_parallel();  % find neurons with fewer nonzero pixels than min_pixel and silent calcium transients
    ids = find(tags>0);
    if ~isempty(ids)
        neuron.viewNeurons(ids, neuron.C_raw);
    end

    % go through all neurons

    % update temporal
    neuron.update_temporal_parallel(use_parallel);

    % update spatial
    neuron.update_spatial_parallel(use_parallel);
end

%% save the workspace for future analysis
neuron.orderROIs('snr');
cnmfe_path = neuron.save_workspace();

%% show neuron contours
Coor = neuron.show_contours(0.6);

%% create a video for displaying the
amp_ac = 140;
range_ac = 5+[0, amp_ac];
multi_factor = 10;
range_Y = 1300+[0, amp_ac*multi_factor];

avi_filename = neuron.show_demixed_video(save_demixed, kt, [], amp_ac, ...
    range_ac, range_Y, multi_factor);

%% save neurons shapes
neuron.save_neurons();

author: zhoupc


text:
function [m] = nm2m(nm)
% Convert length from nanometers to meters.
% Chad A. Greene 2012
m = nm*1e-9;


text:
function [xc,yc,radius,h] = fitcircle_manual(hfig,color,ti)

xc = [];
yc = [];
radius = [];
h = [];

if nargin < 1,
  hfig = gcf;
end
hold on;

if nargin < 2,
  color = 'm';
end

if nargin < 3,
  ti = 'Click points on circle';
end
title(ti);
try
  [x,y] = getline(hfig,'closed');
catch %#ok<CTCH>
  return;
end

if isempty(x),
  return;
end

% solve for parameters a, b, and c in the least-squares sense by
% using the backslash operator
abc = [x y ones(length(x),1)] \ -(x.^2+y.^2);
a = abc(1); b = abc(2); c = abc(3);

% calculate the location of the center and the radius
xc = -a/2;
yc = -b/2;
radius  =  sqrt((xc^2+yc^2)-c);

% display the calculated center
h(1) = plot(xc,yc,[color 'x'],'LineWidth',2);
h(2) = text(xc,yc,sprintf('  (%.1f, %.1f)',xc,yc),'Color',color,'FontWeight','bold');

% plot the entire circle
theta = 0:0.01:2*pi;

% use parametric representation of the circle to obtain coordinates
% of points on the circle
Xfit = radius*cos(theta) + xc;
Yfit = radius*sin(theta) + yc;

h(3) = plot(Xfit, Yfit,[color '-'],'LineWidth',2);

message = sprintf('The estimated radius is %2.3f pixels', radius);
h(4) = text(15,15,message,'Color',color,'FontWeight','bold','BackgroundColor','k');

author: kristinbranson

repo: JAABA


text:
function rho45=rho45(h)
sldata;
limitpoints;
slopes;
rho2=rho1*(T2/T1)^-(1+(g/(m12*R)));
rho3=rho2*exp(-g*(h3-h2)/(R*T2));
rho4=rho3*(T4/T3)^-(1+(g/(m34*R)));
rho45=rho4*(T45(h)/T4)^-(1+(g/(m45*R)));


text:
function str = time2human(t)
% time2human(t)
% t = seconds
% str = output

if t >= 2
  strs = {'sec', 'min', 'hr', 'days', 'weeks', 'months', 'years'};
  divs = [1, 60, 60*60, 60*60*24, 60*60*24*7, 60*60*24*30, 60*60*24*365];

  n = t ./ divs;
  idx = find(n >= 1.5, 1, 'last');
  if isempty(idx)
    idx = 1;
  end
  str = [num2str(n(idx), '%0.3f'), ' ', strs{idx}];
else
  str = [num2str(t*1000, '%0.3f'), ' ms'];
end

repo: CategoryShapes


text:
function tsss = tsss_config
% configuration file for cropping
%__________________________________________________________________________
% Copyright (C) 2014 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: tsss_config.m 7703 2019-11-22 12:06:29Z guillaume $

%--------------------------------------------------------------------------
% D
%--------------------------------------------------------------------------
D        = cfg_files;
D.tag    = 'D';
D.name   = 'File Name';
D.filter = 'mat';
D.num    = [1 1];
D.help   = {'Select the M/EEG mat file.'};

%--------------------------------------------------------------------------
% temporal
%--------------------------------------------------------------------------

temporal = cfg_menu;
temporal.tag = 'temporal';
temporal.name = 'Do temporal denoising';
temporal.labels = {'Yes', 'No'};
temporal.val = {1};
temporal.values = {1,0};
temporal.help = {'Determines whether to use temporal denoising (tSSS)',...
    'if not, just SSS is done.'};

%--------------------------------------------------------------------------
% realign
%--------------------------------------------------------------------------
Dref        = cfg_files;
Dref.tag    = 'Dref';
Dref.name   = 'Reference dataset';
Dref.filter = 'mat';
Dref.num    = [0 1];
Dref.val    = {[]};
Dref.help   = {'Select the M/EEG mat file.',...
    'Leave empty to realign within dataset'};

refind         = cfg_entry;
refind.tag     = 'refind';
refind.name    = 'Reference index';
refind.help    = {'Index of the reference sensors within dataset',...
    'Normally should be 1. Set to 0 for not realigning a composite dataset'};
refind.strtype = 'w';
refind.num     = [1 1];
refind.val     = {1};

realign         = cfg_branch;
realign.tag      = 'realign';
realign.name     = 'Realign head location';
realign.val      = {Dref, refind};
realign.help     = {'Realign head location to another dataset'}';

%--------------------------------------------------------------------------
% timewin
%--------------------------------------------------------------------------
timewin         = cfg_entry;
timewin.tag     = 'timewin';
timewin.name    = 'Time window';
timewin.help    = {'Time window (in sec) for temporal correlation',...
    'Ignored for epoched data'};
timewin.strtype = 'r';
timewin.num     = [1 1];
timewin.val     = {1};

%--------------------------------------------------------------------------
% corrlimit
%--------------------------------------------------------------------------
corrlimit         = cfg_entry;
corrlimit.tag     = 'corrlimit';
corrlimit.name    = 'Correlation limit';
corrlimit.help    = {'Correlation limit parameter for tSSS'};
corrlimit.strtype = 'r';
corrlimit.num     = [1 1];
corrlimit.val     = {0.98};

%--------------------------------------------------------------------------
% Lin
%--------------------------------------------------------------------------
Lin         = cfg_entry;
Lin.tag     = 'Lin';
Lin.name    = 'Inner dimension';
Lin.help    = {'Order if the inner SSS basis'};
Lin.strtype = 'n';
Lin.num     = [1 1];
Lin.val     = {8};

%--------------------------------------------------------------------------
% Lout
%--------------------------------------------------------------------------
Lout         = cfg_entry;
Lout.tag     = 'Lout';
Lout.name    = 'Outer dimension';
Lout.help    = {'Order if the outer SSS basis'};
Lout.strtype = 'n';
Lout.num     = [1 1];
Lout.val     = {3};

%--------------------------------------------------------------------------
% condthresh
%--------------------------------------------------------------------------
condthresh         = cfg_entry;
condthresh.tag     = 'condthresh';
condthresh.name    = 'Condition number threshold';
condthresh.help    = {'Threshold on condition number applied for basis regularisation'};
condthresh.strtype = 'r';
condthresh.num     = [1 1];
condthresh.val     = {80};

%--------------------------------------------------------------------------
% ospace
%--------------------------------------------------------------------------

ospace = cfg_menu;
ospace.tag = 'ospace';
ospace.name = 'Output space';
ospace.labels = {'sensor', 'SSS'};
ospace.val = {0};
ospace.values = {0,1};
ospace.help = {'Determines whether the output file is in sensor space',...
    'or has virtual montage trasforming to SSS space.'};

%--------------------------------------------------------------------------
% prefix
%--------------------------------------------------------------------------
prefix         = cfg_entry;
prefix.tag     = 'prefix';
prefix.name    = 'Filename Prefix';
prefix.help    = {'Specify the string to be prepended to the filenames of the output dataset. Default prefix is ''tsss_''.'};
prefix.strtype = 's';
prefix.num     = [1 Inf];
prefix.val     = {'tsss_'};

%--------------------------------------------------------------------------
% tsss
%--------------------------------------------------------------------------
tsss          = cfg_exbranch;
tsss.tag      = 'tsss';
tsss.name     = 'TSSS denoising';
tsss.val      = {D, temporal, realign, timewin, corrlimit, Lin, Lout, condthresh, ospace, prefix};
tsss.help     = {'TSSS clean-up for Neuromag data'}';
tsss.prog     = @eeg_tsss;
tsss.vout     = @vout_eeg_tsss;
tsss.modality = {'EEG'};

%==========================================================================
function out = eeg_tsss(job)
% construct the S struct
S = [];
S.D = char(job.D);
S.tsss       = job.temporal;
S.Dref       = char(job.realign.Dref);
S.refind     = job.realign.refind;
S.t_window   = job.timewin;
S.corr_limit = job.corrlimit;
S.Lin        = job.Lin;
S.Lout       = job.Lout;
S.cond_threshold = job.condthresh;
S.xspace     = job.ospace;
S.prefix     = job.prefix;
out.D        = tsss_spm_enm(S);
out.Dfname   = {fullfile(out.D)};

%==========================================================================
function dep = vout_eeg_tsss(job)
% return dependencies
dep(1)            = cfg_dep;
dep(1).sname      = 'TSSS-ed MEG data';
dep(1).src_output = substruct('.','D');
dep(1).tgt_spec   = cfg_findspec({{'strtype','e'}});

dep(2)            = cfg_dep;
dep(2).sname      = 'TSSS-ed MEG datafile';
dep(2).src_output = substruct('.','Dfname');
dep(2).tgt_spec   = cfg_findspec({{'filter','mat'}});

author: spm


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

%%
clear; close all; clc;

%% Syntax
% |axisGeom(h,fontSize);|

%% Description
% This function sets axis properties commonly used for viewing geometry in
% 3D:
%
% |axis equal; axis vis3d; axis tight;|
%
% The optional inputs include the axis handle and the desired font size.

%% Examples
%

%% Setting axis parameters for 3D geometry viewing

%%
% Create something to visualize
[F,V]=graphicsModels(9);

%%
% Visualize

cFigure; hold on;
gpatch(F,V,'w','k');

axisGeom; %Set axis properties

camlight headlight;
gdrawnow;

%%
%
% <<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 [VecTot, TOT, time] = rmt(P, Vec, Mat)
%rmt returns a matrix containing the best models for all the paths of the Replacement Method.
%TOT contains all the relative results showing the evolution of the method.
%
%
%	   Input:
%             P             Property vector
%             Vec           Initial descriptors vector
%             Mat           Descriptors matrix with descriptors pool
%
%
%     Returns:
%
%            VecTot           vector containing the best model for all the
%                               paths of the Replacement Method
%            TOT               contains all the relative results
%                               showing the evolution of the method.
%
% Andrew G. Mercader, Pablo R. Duchowicz
% INIFTA, La Plata, Argentina
% Created: 5 March 2007

TOT=[];
VecTot=[];
time=cputime;
warning off

[k_v,n_v]=size(Vec);

for k=1:n_v

Sr=rms(P,Vec,Mat);
TOT(k).A(1,:)=[Sr,Vec];

VecA=rmsr(P, Vec, Mat, k);
Po(1)=k;
VecI=VecA;
if n_v==1
    VecTot=[1,VecI];
    TOT=VecI;
    time=cputime-time
    return
end
VecI(1)=[];
COEF=rmder(P,VecI,Mat);
COER=COEF;
COER(Po)=[];
pos=find(COEF==max(COER));
Po(2)=pos;
TOT(k).A(2,:)=VecA;
for i=2:n_v;
    VecA=rmsr(P, VecI, Mat, pos);
    VecI=VecA;
    TOT(k).A(i+1,:)=VecA;
    VecI(1)=[];
    COEF=rmder(P,VecI,Mat);
    if i==n_v
        Po=[];
        break
    end
    COER=COEF;
    COER(Po)=[];
    pos=find(COEF==max(COER));
    Po(i+1)=pos;
end

for j=1:2;
COER=COEF;
COER(Po)=[];
pos=find(COEF==max(COER));
Po(1)=pos;
for i=1:n_v;
    VecA=rmsr(P, VecI, Mat, pos);
    VecI=VecA;
    TOT(k).A(i+(j*n_v),:)=VecA;
    VecI(1)=[];
    COEF=rmder(P,VecI,Mat);
    if i==n_v
        Po=[];
        break
    end
    COER=COEF;
    COER(Po)=[];
    pos=find(COEF==max(COER));
    Po(i+1)=pos;
end
end

for j=3:100;
COER=COEF;
COER(Po)=[];
pos=find(COEF==max(COER));
Po(1)=pos;
    for i=1:n_v;
    VecA=rmsr(P, VecI, Mat, pos);
    VecI=VecA;
    TOT(k).A(i+(j*n_v),:)=VecA;
    VecI(1)=[];
    COEF=rmder(P,VecI,Mat);
    if i==n_v
        Po=[];
        break
    end
    COER=COEF;
    COER(Po)=[];
    pos=find(COEF==max(COER));
    Po(i+1)=pos;
    end
   if TOT(k).A(i+(j*n_v),:)==TOT(k).A(i+(j*n_v)-(2*n_v),:)
        Po=[];
        break
   end

end
VecP=find(TOT(k).A==min(TOT(k).A(:,1)));
VecTot(k,:)=[k,TOT(k).A(VecP(1),:)];
VecTot=sortrows(VecTot,2);
end
time=cputime-time
warning on

% % End of rmt
%


text:
classdef LvdOptimizerAlgoEnum < matlab.mixin.SetGet
    %LvdOptimizerAlgoEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        Fmincon('FMINCON', ...
                'The standard gradient-based nonlinear programming solver in MATLAB.  Supports a variety of subroutines.  Try this algorithm first.', ...
                true);

        PatternSearch('Pattern Search', ...
                      'Uses a "direct search" algorithm to solve NLP problems.  Direct search is a method for solving optimization problems that does not require any information about the gradient of the objective function.  Use this function to find a reasonable guess for a solution and then switch to a gradient-based method.', ...
                      false);

        Nomad('NOMAD', ...
              'NOMAD is a C++ implementation of the Mesh Adaptive Direct Search (MADS) algorithm designed for constrained optimization of blackbox functions.  At the core of NOMAD resides the Mesh Adaptive Direct Search (MADS) algorithm.', ...
              false);

        Ipopt('IPOPT', ...
              'IPOPT is a software library for large scale nonlinear optimization of continuous systems.', ...
              true);

        Surrogate('Surrogate Optimizer', ...
              'Surrogate optimization for global minimization of time-consuming objective functions.  This optimizer does not necessarily handle equality constraints well.', ...
              false);
    end

    properties
        name(1,:) char
        desc(1,:) char
        reqGrad(1,1) logical = true
    end

    methods
        function obj = LvdOptimizerAlgoEnum(name, desc, reqGrad)
            obj.name = name;
            obj.desc = desc;
            obj.reqGrad = reqGrad;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('LvdOptimizerAlgoEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('LvdOptimizerAlgoEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('LvdOptimizerAlgoEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function reversed_sequence = sequence_reverse(sequence)
% sequence_reverse Returns reverse sequence
%
% This sequence converter returns a sequence with all its frames reversed.
%
% Input:
% - sequence (structure): A valid sequence structure.
%
% Output:
% - reversed_sequence (structure): A sequence descriptor of a converted sequence.

reversed_sequence = sequence;
reversed_sequence.groundtruth = sequence.groundtruth(end:-1:1);
reversed_sequence.indices = sequence.indices(end:-1:1);
reversed_sequence.images = sequence.images(end:-1:1);
reversed_sequence.tags.data = sequence.tags.data(end:-1:1, :);
reversed_sequence.values.data = sequence.values.data(end:-1:1, :);

author: votchallenge

repo: toolkit-legacy


text:
function [x_new,P_new,w_new] = gauss_merge(x,P,w,threshold)
% gauss_merge Merge mixture components that fall in close proximity
%
% TODO:
% ----
% * Generate documentation
%
% Credit:
% -------
% * This code has been taken from http://ba-tuong.vo-au.com/codes.html.
% * Minor modifications may have been applied to the original version.

L= length(w); x_dim= size(x,1);
I= 1:L;
el= 1;

if all(w==0)
    x_new = [];
    P_new = [];
    w_new = [];
    return;
end

while ~isempty(I)
    [notused,j]= max(w); j= j(1);
    Ij= []; iPt= inv(P(:,:,j));
    x_new(:,el)= zeros(x_dim,1);
    P_new(:,:,el)= zeros(x_dim,x_dim);
    w_new(1,el)= 0;
    for i= I
        val= (x(:,i)-x(:,j))'*iPt*(x(:,i)-x(:,j));
        if val <= threshold
            Ij= [ Ij i ];
        end
    end

    x_new(:,el)= wsumvec(w(Ij),x(:,Ij),x_dim);
    P_new(:,:,el)= wsummat(w(Ij)',P(:,:,Ij),x_dim);
    w_new(1,el)= sum(w(Ij));

    x_new(:,el)= x_new(:,el)/w_new(el);
    P_new(:,:,el)= P_new(:,:,el)/w_new(el);
    I= setdiff(I,Ij);
    w(Ij)= -1;
    el= el+1;
end

function out = wsumvec(w,vecstack,xdim)
    wmat = repmat(w,[xdim,1]);
    out  = sum(wmat.*vecstack,2);

function out = wsummat(w,matstack,xdim)
    w = reshape(w,[1,1,size(w)]);
    wmat = repmat(w,[xdim,xdim,1]);
    out = sum(wmat.*matstack,3);

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function createBasePanel(handles)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC stateS;
indexS = planC{end};

    %create a tree view and initialize filter parameters;
    w = 0.3; h = 0.6;
    createTreeView(handles, w, h);
    handles = guidata(handles.mainframe);

    fsize = 9;

    %create output list panel
    handles.outputPanel = uipanel('Parent',handles.mainframe,'Units','characters','FontSize',fsize,'Title','Output','Tag','OutputPanel',...
                    'Clipping','on','Units','normalized','Position',[0 0 w 1-h]);

    handles.OutputList = uicontrol('Parent',handles.outputPanel,'Units','normalized',...
                    'FontSize',fsize,'Position',[0 0 1 1],'Style','listbox', 'Tag','OutputList');

    %add generalPanel,toolPanel and optionPanel
    w = 0.3; h1 = 0.30; h2 = 0.6; h3 = 1-h1-h2;
    handles.generalPanel = uipanel('parent', handles.mainframe, 'units', 'normalize', 'FontSize',fsize, 'position', [w h2+h3 1-w h1]);

    handles.toolPanel = uipanel('parent', handles.mainframe, 'FontSize',fsize, 'units', 'normalize', 'position', [w h3 1-w h2]);

    handles.optionPanel = uipanel('parent', handles.mainframe, 'FontSize',fsize, 'Title','Options', 'units', 'normalize', ...
                                  'position', [w 0 1-w h3]);

    %do button
    handles.doButton = uicontrol(handles.optionPanel, 'style', 'pushbutton', 'units', 'pixel', ...
                            'position', [310 10 100 28], 'string', 'Continue', ...
                            'tag', 'continueButton', 'callback', 'CERRRegistrationRigidSetup(''auto_registration'')');

% Options
    %downsample option
    dy = 20;
    handles.dsampleCheck = uicontrol(handles.optionPanel, 'style', 'checkbox', 'units', 'pixel', 'value', 0, ...
                                'position', [2 24 136 20], 'string', 'downSample(2x2x2)', 'tag', 'dsamplecheckbox');
    handles.saveCheck = uicontrol(handles.optionPanel, 'style', 'checkbox', 'units', 'pixel', 'Value', 0, ...
                                'position', [2 2 136 20], 'string', 'save resampled data', 'tag', 'savecheckbox');

    InitTranstext = uicontrol(handles.optionPanel, 'style', 'text', 'units', 'pixel', ...
                                'position', [138 1 50 16], 'string', 'InitTrans:');
    handles.InitTrans = uicontrol(handles.optionPanel, 'style', 'popup', 'units', 'pixel', 'value', 1, ...
                                'position', [188 4 86 18], 'string', 'MomentsOn|GeometryOn|InitTransMOn','tag', 'InitTrans');
%
    handles.flipX = uicontrol(handles.optionPanel, 'style', 'checkbox', 'units', 'pixel', 'Value', 0, ...
                                'position', [138 24 50 20], 'string', 'flipX', 'tag', 'flipX', 'enable', 'off');
    handles.flipY = uicontrol(handles.optionPanel, 'style', 'checkbox', 'units', 'pixel', 'Value', 0, ...
                                'position', [188 24 50 20], 'string', 'flipY', 'tag', 'flipY', 'enable', 'off');
    handles.flipZ = uicontrol(handles.optionPanel, 'style', 'checkbox', 'units', 'pixel', 'Value', 0, ...
                                'position', [238 24 50 20], 'string', 'flipZ', 'tag', 'flipZ', 'enable', 'off');

%     2D Match
%     handles.match2DCheck = uicontrol(handles.optionPanel, 'style', 'checkbox', 'units', 'pixel', 'Value', 0, ...
%                                 'position', [138 15 80 28], 'string', '2D
%                                 Match', 'tag', 'match2Dcheckbox');

%registration method panel
    handles.registrationMethodPanel = uipanel('Parent',handles.generalPanel,'Units','normalize', 'FontSize',fsize,...
                                'Title','Image Registration Method','Tag','regMethodPanel','Clipping','off',...
                                'Position',[0 0 1 1], 'Visible','on');

    dy = -1; dx = -5;
    h21 = uicontrol(...
    'Parent',handles.registrationMethodPanel,...
    'Units','characters',...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[6+dx 10.2307692307692+dy 20.6 2],...
    'String','Transform:',...
    'Style','text',...
    'Tag','text6');

    handles.transform = uicontrol(...
    'Parent',handles.registrationMethodPanel,...
    'Units','characters',...
    'FontSize',10,...
    'Position',[28.2+dx 10.1538461538462+dy 51 2.07692307692308],...
    'String',{ 'Similarity Transform'; 'Versor3D Transform'; 'Euler Transform';  },... %'Affine Transform'
    'Style','popupmenu',...
    'CallBack', 'createParaPanel(guidata(gcbf), ''transform_advance'')', ...
    'Value',1);

%     try
%         [x, map] = imread('help-book-open.gif', 'GIF');
%         a = ind2rgb(x, map);
%     catch
%         a = [];
%     end
%     handles.helpButton = uicontrol(handles.registrationMethodPanel, 'style', 'pushbutton', 'units', 'characters', ...
%                             'cdata', a, 'position', [80+dx 10.1538461538462+dy 6 2.07692307692308], 'string', '', ...
%                             'tag', 'continueButton', 'callback', 'winopen(''irqat.chm'')');

    h23 = uicontrol(...
    'Parent',handles.registrationMethodPanel,...
    'Units','characters',...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[6+dx 7.38461538461539+dy 20.6 2],...
    'String','Interpolator:',...
    'Style','text',...
    'Tag','text7');

    handles.interpolater = uicontrol(...
    'Parent',handles.registrationMethodPanel,...
    'Units','characters',...
    'FontSize',10,...
    'Position',[28.2+dx 7.30769230769231+dy 51 2.07692307692308],...
    'String',{  'Linear Interpolate'; 'Nearest Interpolate'},...
    'Style','popupmenu',...
    'Value',1);

    h25 = uicontrol(...
    'Parent',handles.registrationMethodPanel,...
    'Units','characters',...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[6+dx 4.61538461538462+dy 20.6 2],...
    'String','Similarity Metric:',...
    'Style','text',...
    'Tag','text8');

    handles.metric = uicontrol(...
    'Parent',handles.registrationMethodPanel,...
    'Units','characters',...
    'FontSize',10,...
    'Position',[28.2+dx 4.53846153846154+dy 51 2.07692307692308],...
    'String',{  'Mean Squares'; 'Normalized Correlation' },...
    'Style','popupmenu',...
    'Value',1);

    h27 = uicontrol(...
    'Parent',handles.registrationMethodPanel,...
    'Units','characters',...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[6+dx 1.92307692307692+dy 20.6 2],...
    'String','Optimizer:',...
    'Style','text',...
    'Tag','text9');

    handles.optimizer = uicontrol(...
    'Parent',handles.registrationMethodPanel,...
    'Units','characters',...
    'FontSize',10,...
    'Position',[28.2+dx 1.84615384615385+dy 51 2.07692307692308],...
    'String','RegularStepGradientDescent',...
    'Style','popupmenu',...
    'Value',1);

    %save handles
    guidata(handles.mainframe, handles);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function fig = guipar()
% This is the machine-generated representation of a Handle Graphics object
% and its children.  Note that handle values may change when these objects
% are re-created. This may cause problems with any callbacks written to
% depend on the value of the handle at the time the object was saved.
% This problem is solved by saving the output as a FIG-file.
%
% To reopen this object, just type the name of the M-file at the MATLAB
% prompt. The M-file and its associated MAT-file must be on your path.
%
% NOTE: certain newer features in MATLAB may not have been saved in this
% M-file due to limitations of this format, which has been superseded by
% FIG-files.  Figures which have been annotated using the plot editor tools
% are incompatible with the M-file/MAT-file format, and should be saved as
% FIG-files.

load guipar

h0 = figure('BackingStore','off', ...
	'Color',[0.8 0.8 0.8], ...
	'Colormap',mat0, ...
	'FileName','C:\Program\Matlab53\m-filer\lab\guipar.m', ...
	'MenuBar','none', ...
	'Name','Simulator Parameters', ...
	'NumberTitle','off', ...
	'PaperPosition',[18 180 576 432], ...
	'PaperUnits','points', ...
	'Position',[530 162 369 201], ...
	'Tag','Fig2', ...
	'ToolBar','none');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'Callback',mat1, ...
	'ListboxTop',0, ...
	'Position',[288 244.5 32.25 13.5], ...
	'String','2000', ...
	'Style','edit', ...
	'Tag','k');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'ListboxTop',0, ...
	'Position',[231 239.25 54 21.75], ...
	'String','Heat Transfer Coefficient', ...
	'Style','text', ...
	'Tag','StaticText6');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'ListboxTop',0, ...
	'Position',[10 10 38.25 19.5], ...
	'String','Holdup volume', ...
	'Style','text', ...
	'Tag','StaticText7');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'ListboxTop',0, ...
	'Position',[10 10 34.5 21], ...
	'String','No of nodes', ...
	'Style','text', ...
	'Tag','StaticText8');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'ListboxTop',0, ...
	'Position',[10 10 32.25 13.5], ...
	'String','20', ...
	'Style','edit', ...
	'Tag','n');
if nargout > 0, fig = h0; end


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni

%problem 1

%x(t) =t*exp(-t), 0<t<T , T=6

%x(t),0<t<4T
t=0:.1:6;
x=t.*exp(-t);
xp=repmat(x,1,4);
tp=linspace(-6,18,length(xp));
plot(tp,xp)
legend('x(t) in 4 periods')

%'Approximations
figure %exponential
t0=0;
T=6;
w=2*pi/T;
syms t
x=t.*exp(-t);
k=-40:40;
a=(1/T)*int(x*exp(-j*k*w*t), t,t0,t0+T);
xx=sum(a.*exp(j*k*w*t));
ezplot(xx,[-6 18]);
legend('approximate signal in 4 periods using a_k')

figure %trigonometric
a0=(1/T)*int(x,t0,t0+T);
n=1:81;
b=(2/T)*int(x*cos(n*w*t),t,t0,t0+T);
c=(2/T)*int(x*sin(n*w*t),t,t0,t0+T);
xx=sum(a.*exp(j*k*w*t));
ezplot(xx,[-6 18]);
legend('approximate signal in 4 periods using b_k ,c_k' )


text:
% Do the example from Satnam Alag's PhD thesis, UCB ME dept 1996 p46

% Make the following polytree, where all arcs point down

% 1   2
%  \ /
%   3
%  / \
% 4   5

N = 5;
dag = zeros(N,N);
dag(1,3) = 1;
dag(2,3) = 1;
dag(3, [4 5]) = 1;

ns = [2 1 2 1 2];

bnet = mk_bnet(dag, ns, 'discrete', []);

bnet.CPD{1} = gaussian_CPD(bnet, 1, 'mean', [1 0]', 'cov', [4 1; 1 4]);
bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', 1, 'cov', 1);
B1 = [1 2; 1 0]; B2 = [2 1]';
bnet.CPD{3} = gaussian_CPD(bnet, 3, 'mean', [0 0]', 'cov', [2 1; 1 1], ...
			   'weights', [B1 B2]);
H1 = [1 1];
bnet.CPD{4} = gaussian_CPD(bnet, 4, 'mean', 0, 'cov', 1, 'weights', H1);
H2 = [1 0; 1 1];
bnet.CPD{5} = gaussian_CPD(bnet, 5, 'mean', [0 0]', 'cov', eye(2), 'weights', H2);

engine = {};
engine{end+1} = jtree_inf_engine(bnet);
engine{end+1} = pearl_inf_engine(bnet, 'protocol', 'tree');
engine{end+1} = pearl_inf_engine(bnet, 'protocol', 'parallel');
E = length(engine);

if 1
% no evidence
evidence = cell(1,N);
ll = zeros(1,E);
for e=1:E
  [engine{e}, ll(e)] = enter_evidence(engine{e}, evidence);
  add_ev = 1;
  m = marginal_nodes(engine{e}, 3, add_ev);
  assert(approxeq(m.mu, [3 2]'))
  assert(approxeq(m.Sigma, [30 9; 9 6]))

  m = marginal_nodes(engine{e}, 4, add_ev);
  assert(approxeq(m.mu, 5))
  assert(approxeq(m.Sigma, 55))

  m = marginal_nodes(engine{e}, 5, add_ev);
  assert(approxeq(m.mu, [3 5]'))
  assert(approxeq(m.Sigma, [31 39; 39 55]))
end
end

if 1
% evidence on leaf 5
evidence = cell(1,N);
evidence{5} = [5 5]';
for e=1:E
  [engine{e}, ll(e)] = enter_evidence(engine{e}, evidence);
  add_ev = 1;
  m = marginal_nodes(engine{e}, 3, add_ev);
  assert(approxeq(m.mu, [4.4022 1.0217]'))
  assert(approxeq(m.Sigma, [0.7011 -0.4891; -0.4891 1.1087]))

  m = marginal_nodes(engine{e}, 4, add_ev);
  assert(approxeq(m.mu, 5.4239))
  assert(approxeq(m.Sigma, 1.8315))

  m = marginal_nodes(engine{e}, 1, add_ev);
  assert(approxeq(m.mu, [0.3478 1.1413]'))
  assert(approxeq(m.Sigma, [1.8261 -0.1957; -0.1957 1.0924]))

  m = marginal_nodes(engine{e}, 2, add_ev);
  assert(approxeq(m.mu, 0.9239))
  assert(approxeq(m.Sigma, 0.8315))

  m = marginal_nodes(engine{e}, 5, add_ev);
  assert(approxeq(m.mu, evidence{5}))
  assert(approxeq(m.Sigma, zeros(2)))
end
end

if 1
% evidence on leaf 4 (non-info-state version is uninvertible)
evidence = cell(1,N);
evidence{4} = 10;
for e=1:E
  [engine{e}, ll(e)] = enter_evidence(engine{e}, evidence);
  add_ev = 1;
  m = marginal_nodes(engine{e}, 3, add_ev);
  assert(approxeq(m.mu, [6.5455 3.3636]'))
  assert(approxeq(m.Sigma, [2.3455 -1.6364; -1.6364 1.9091]))

  m = marginal_nodes(engine{e}, 5, add_ev);
  assert(approxeq(m.mu, [6.5455 9.9091]'))
  assert(approxeq(m.Sigma, [3.3455 0.7091; 0.7091 1.9818]))

  m = marginal_nodes(engine{e}, 1, add_ev);
  assert(approxeq(m.mu, [1.9091 0.9091]'))
  assert(approxeq(m.Sigma, [2.1818 -0.8182; -0.8182 2.1818]))

  m = marginal_nodes(engine{e}, 2, add_ev);
  assert(approxeq(m.mu, 1.2727))
  assert(approxeq(m.Sigma, 0.8364))
end
end

if 1
% evidence on leaves 4,5 and root 2
evidence = cell(1,N);
evidence{2} = 0;
evidence{4} = 10;
evidence{5} = [5 5]';
for e=1:E
  [engine{e}, ll(e)] = enter_evidence(engine{e}, evidence);
  add_ev = 1;
  m = marginal_nodes(engine{e}, 3, add_ev);
  assert(approxeq(m.mu, [4.9964 2.4444]'));
  assert(approxeq(m.Sigma, [0.6738 -0.5556; -0.5556 0.8889]));

  m = marginal_nodes(engine{e}, 1, add_ev);
  assert(approxeq(m.mu, [2.2043 1.2151]'));
  assert(approxeq(m.Sigma, [1.2903 -0.4839; -0.4839 0.8065]));
end
end

if 1
  [time, engine] = cmp_inference_static(bnet, engine, 'maximize', 0, 'check_ll', 0, ...
				     'singletons_only', 0, 'observed', [1 3 5]);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [angleRate] = PTrc2deg(X,rcRate,rcExpo,superrate, C)
% raw RCcommand data to RCrate in deg/s, i.e. "set point"
%  with help from: https://github.com/betaflight/betaflight-configurator/blob/master/src/js/RateCurve.js
%   X is a vector containing a single axis of RCcommand data scaled from -500 to 500,
%   rcRate(0-255),rcExpo(0-100),superrate(0-100)

    expoPower=3;
    rcRateConstant=C;
    angleRate=[];

    rcRate=rcRate/100;
    rcExpo=(rcExpo/100);
    superrate=superrate/100;

    maxRC=500;
    rcCommandf = X / maxRC;% scales from -1 to 1 for the math
    rcCommandfAbs = abs(rcCommandf) / 1;%max(abs(rcCommandf));

    if (rcRate > 2)
        rcRate = rcRate + (rcRate - 2) * 14.54;
    end
    if (rcExpo > 0)
       %disp('rcExpo > 0')
        rcCommandf =  rcCommandf .* power(rcCommandfAbs, expoPower) * rcExpo + rcCommandf * (1-rcExpo);
    end
    if (superrate > 0)
        %disp('superrate > 0')
        rcFactor = 1 ./ (1 - rcCommandfAbs * superrate); % this creates the super expo curve needed to convert RCcommand to RCrate
        angleRate = (rcRateConstant * rcRate * rcCommandf);
        angleRate = angleRate .* rcFactor;
       % disp(['angleRate:' num2str(angleRate) ' rcFactor:' num2str(rcFactor)])
    else
        angleRate = (rcRateConstant * rcRate * rcCommandf);
    end

end

repo: PIDtoolbox