text:
function f = compose(f, op, varargin)
%COMPOSE   Compose command for SPHEREFUN objects.
%   F = COMPOSE(F, OP)  returns the SPHEREFUN that approximates OP(F).
%
%   F = COMPOSE(F, OP, G)  returns the SPHEREFUN that approximates OP(F).
%
%   F = COMPOSE(F, G) with a CHEBFUN G with one column returns a SPHEREFUN that
%   approximates G(F).  If G has 3 columns, the result is a SPHEREFUNV. If G is
%   a CHEBFUN2 or CHEBFUN2V, the composition is interpreted as G(real(F),
%   imag(F)).
%
%   This command is a wrapper for the SPHEREFUN constructor.

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

if ( isempty(op) )
    return
elseif ( isempty(f) )
    f = op;

elseif ( isa(op, 'chebfun') )
    % Composition OP(f) of SPHEREFUN object f and CHEBFUN OP

    if ( length(op.domain) > 2 )
        % If OP has several pieces, OP(SPHEREFUN) might be inaccurate.
        warning('CHEBFUN:SPHEREFUN:compose:pieces', ...
            ['The composition of a CHEBFUN with several pieces and a SPHEREFUN\n', ...
            'might be inaccurate.']);
    end

    % Check that image(f) is contained in domain(OP).
    vals = minandmax2est(f);        % Estimate of image(f).
    tol = 100 * chebfun2eps * max(vscale(f), vscale(op)) * ...
            norm(f.domain, inf);    % Tolerance.
    if ( ~isSubset(vals, op.domain, tol) )
        error('CHEBFUN:SPHEREFUN:COMPOSE:DomainMismatch', ...
            'OP(F) is not defined, since image(F) is not contained in domain(OP).')
    end

    nColumns = size(op, 2);
    if ( nColumns == 1 )
        % Call constructor:
        f = spherefun(@(x,y) op(feval(f, x, y)), f.domain);

    elseif ( nColumns == 3 )
        % Extract columns of the CHEBFUN OP:
        op1 = op(:,1);
        op2 = op(:,2);
        op3 = op(:,3);

        % Call constructor:
        f = spherefunv(@(x,y) op1(feval(f, x, y)), ...
            @(x,y) op2(feval(f, x, y)), @(x,y) op3(feval(f, x, y)));

    else
        % The CHEBFUN object OP has a wrong number of columns.
        error('CHEBFUN:SPHEREFUN:COMPOSE:Columns', ...
            'The CHEBFUN object must have 1 or 3 columns.')

    end

elseif ( isa(op, 'chebfun2') )
    % Composition OP(f) of SPHEREFUN object f and CHEBFUN2 OP, interpreted as
    % OP(real(f), imag(f)).  For now SPHEREFUNS are real, but this might change
    % in the future.

    % Check that image(f) is contained in domain(OP).
    vals = minandmax2est(f);        % Estimate of image(f).
    tol = 100 * chebfun2eps * max(vscale(f), vscale(op)) * ...
            norm(f.domain, inf);    % Tolerance.
    if ( ~isSubset(vals, op.domain(1:2), tol) )
        error('CHEBFUN:SPHEREFUN:COMPOSE:DomainMismatch2', ...
            'OP(F) is not defined, since image(F) is not contained in domain(OP).')
    end

    % Call constructor:
    f = spherefun(@(x,y) op(feval(real(f), x, y), feval(imag(f), x, y)), ...
        f.domain);

elseif ( isa(op, 'chebfun2v') )
    % Composition OP(f) of SPHEREFUN object f and CHEBFUN2V OP with three
    % components, interpreted as OP(real(f), imag(f)).
    % For now SPHEREFUNS are real, but this might change in the future.

    % Check that OP has three components:
    if ( op.nComponents ~= 3 )
        error('CHEBFUN:SPHEREFUN:COMPOSE:C2VnComponents', ...
            'The Chebfun2v objects must have three components.')
    end

    % Get the components:
    op1 = op(1);
    op2 = op(2);
    op3 = op(3);

    % Check that image(f) is contained in domain(OP).
    vals = minandmax2est(f);        % Estimate of image(f).
    tol = 100 * chebfun2eps * max(vscale(f), vscale(op)) * ...
            norm(f.domain, inf);    % Tolerance.
    if ( ~isSubset(vals, op1.domain(1:2), tol) )
        error('CHEBFUN:SPHEREFUN:COMPOSE:DomainMismatch2v', ...
            'OP(F) is not defined, since image(F) is not contained in domain(OP).')
    end

    % Get real and imaginary parts of f:
    realf = real(f);
    imagf = imag(f);

    % Call constructor:
    f = spherefunv(@(x,y) op1(feval(realf, x, y), feval(imagf, x, y)), ...
        @(x,y) op2(feval(realf, x, y), feval(imagf, x, y)), ...
        @(x,y) op3(feval(realf, x, y), feval(imagf, x, y)));

elseif ( nargin == 2 && nargin(op) == 1 )
    % OP has one input variable.

    % Call constructor:
    f = spherefun(@(x,y) op( feval(f, x, y) ), f.domain);

elseif ( nargin == 3 && nargin(op) == 2 )
    % OP has two input variables.

    g = varargin{1};
    if ( isa(g, 'double') )     % promote
        g = spherefun(@(x,y,z) g + 0*x, f.domain);
    end

    if ( isa(f, 'double') )     % promote
        f = spherefun(@(x,y,z) f + 0*x, g.domain);
    end

    % Call constructor:
    f = spherefun(@(x,y) op( feval(f, x, y), feval(g, x, y) ), f.domain);

else
    % Not sure what to do, error:
    error('CHEBFUN:SPHEREFUN:COMPOSE:OP', 'NARGIN(OP) not correct.')
end

end

author: chebfun

repo: chebfun


text:
function planC = postProcStruct(planC,userOptS)
% planC = postProcStruct(planC,userOptS);
% Function to post-process segementations.
%-------------------------------------------------------------------------
%
%
% AI 10/14/19

% Get list of structures to be post-processed
postS = [];
if isfield(userOptS,'postProc')
    postS = userOptS.postProc;
end

if ~isempty(postS)

    strC = fieldnames(postS);

    indexS = planC{end};

    %Loop over structures from postS
    for iStr = 1:length(strC)

        outMask3M = [];

        methodC = postS.(strC{iStr}); %.method;
        if ~iscell(postS.(strC{iStr})) %.method)
            methodCStruct = methodC;
            methodC = {};
            for i = 1:numel(methodCStruct)
              methodC{i} = methodCStruct(i);
            end
        end

        maskC = cell(length(methodC),1);

        %Loop over post-processing methods
        for iMethod = 1:length(methodC)
            strListC = {planC{indexS.structures}.structureName};
            %get updated structure number in case iMethod > 1
            strNum = getMatchingIndex(strC{iStr},strListC,'EXACT');
            strNum = strNum(end); %Use most recent str if >1 matches found
            scanNum = getStructureAssociatedScan(strNum,planC);

            switch methodC{iMethod}.method

                case 'getLargestConnComps'

                    numConnComponents = methodC{iMethod}.params.numCC;
                    [maskC{iMethod}, planC]= getLargestConnComps(strNum,numConnComponents,planC);

                case 'getLargestOverlappingComp'

                    roiName = methodC{iMethod}.params.roiName;
                    roiStrNum = getMatchingIndex(roiName,strListC,'EXACT');
                    [maskC{iMethod}, planC] = getLargestOverlappingComp(strNum,roiStrNum,planC);

                case 'getSegInROI'

                    roiName = methodC{iMethod}.params.roiName;
                    roiStrNum = getMatchingIndex(roiName,strListC,'EXACT');

                    [roiMask3M, planC] = getStrMask(roiStrNum,planC);
                    [strMask3M, planC] = getStrMask(strNum,planC);

                    maskC{iMethod} = roiMask3M & strMask3M;

                case 'removeBackgroundFP'
                    ctOffset = planC{indexS.scan}(scanNum).scanInfo(1).CTOffset;
                    scan3M = double(getScanArray(scanNum,planC)) - ctOffset;
                    threshold = -400; % default for CT
                    if isfield(methodC{iMethod},'params') && isfield(methodC{iMethod}.params,'threshold')
                        threshold = methodC{iMethod}.params.threshold;
                    end
                    connPtMask3M = getPatientOutline(scan3M,[],threshold);
                    [strMask3M, planC] = getStrMask(strNum,planC);

                    maskC{iMethod} = strMask3M & connPtMask3M;

                case 'none'

                    [maskC{iMethod}, planC] = getStrMask(strNum,planC);

                otherwise
                    %Custom post-processing function
                    customMethod = methodC{iMethod};
                    if isfield(customMethod,'params')
                        paramS = customMethod.params;
                    else
                        paramS = [];
                    end
                    [maskC{iMethod},planC] = feval(customMethod.method,strNum,paramS,planC);

            end

            %Combine masks
            if iMethod>1
                switch lower(methodC{iMethod}.operator)
                    case 'union'
                        outMask3M = or(maskC{iMethod-1},maskC{iMethod});
                        maskC{iMethod} = outMask3M;
                    case 'intersection'
                        outMask3M = and(maskC{iMethod-1},maskC{iMethod});
                        maskC{iMethod} = outMask3M;
                end
            end

            %Replace original structure
            outMask3M = maskC{iMethod};
            planC = maskToCERRStructure(outMask3M,0,scanNum,strC{iStr},planC);
            roiGenAlg = planC{indexS.structures}(strNum).roiGenerationAlgorithm;
            roiGenDesc = planC{indexS.structures}(strNum).roiGenerationDescription;
            planC{indexS.structures}(end).roiGenerationAlgorithm = roiGenAlg;
            planC{indexS.structures}(end).roiGenerationDescription = roiGenDesc;
            planC{indexS.structures}(end).structureDescription = roiGenDesc;
            planC = deleteStructure(planC,strNum);
        end

    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [H, WtW, WtV, infos] = nnls_solver(V, W, in_options)
%
% This file is part of NMFLibrary.
%
% Created by H.Kasai on June 20, 2022
%
% Change log:
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.alg_name = 'no_name';
    local_options.algo = 'hals';
    local_options.delta = 1e-6;

    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize
    W = full(W);
    WtW = W' * W;
    WtV = W' * V;
    infos = [];

    % If no initial matrices are provided, H is initialized as follows:
    if ~isfield(options,'init') || isempty(options.init)
        H = nnls_init_nmflibrary(V, W, WtW, WtV);
    else
        H = options.init;
    end

    % switch solvers

    switch options.algo
        case 'hals'

            % perform main routine
            options.V = V;
            options.W = W;
            options.main_routine_mode = true;
            %options.max_epoch = 500;
            options.eps0 = 0;
            options.eps = 1;
            eit1 = cputime;
            alpha = 2;

            % call HALSupdt
            [H, infos] = HALSupdt(H, WtW, WtV, eit1, alpha, options.delta, options);

            if options.verbose > 2
                fprintf('nnls_solver [%d]: H = %.16e, WtW = %.16e, WtV = %.16e\n', length(infos.cost), norm(H), norm(WtW), norm(WtV));
            end

        case 'fpgm'

            options_fpgm.init = H;
            options_fpgm.inner_max_epoch = options.inner_max_epoch;

            % call HALSupdt
            [H, WtW, WtV] = nnls_fpgm(V, W, options_fpgm);

        case 'anls_bpp'

            % call HALSupdt
            H = nnlsm_blockpivot(WtW, WtV, 1, H);

        case 'anls_asgroup'

            % call HALSupdt
            [H, ~, ~] = nnlsm_activeset(WtW, WtV, 0, 1, H);

        otherwise
    end

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function data = loadSUN3Dv2(sequenceName, frameIDs)

if ~exist('sequenceName','var')
    sequenceName = '2014-04-29_14-39-49_094959634447';
end

SUN3Dpath = '/n/fs/sun3d/sun3dv2/';

%{
fileID = fopen(fullfile(SUN3Dpath,sequenceName,'image/time.dat'));
imageTimestamp = fread(fileID,'int64');
fclose(fileID);
data.image = VideoReader(fullfile(SUN3Dpath,sequenceName,'image/image.mp4'));
%}

imageFiles = dirSmart(fullfile(SUN3Dpath,sequenceName,'image/'),'jpg');
imageFrameID = zeros(1,length(imageFiles));
imageTimestamp = zeros(1,length(imageFiles));
for i=1:length(imageFiles)
    id_time = sscanf(imageFiles(i).name, '%d-%ld.jpg');
    imageFrameID(i) = id_time(1);
    imageTimestamp(i) = id_time(2);

    data.imageAll{i}= fullfile(fullfile(SUN3Dpath,sequenceName,'image',imageFiles(i).name));
end

depthFiles = dirSmart(fullfile(SUN3Dpath,sequenceName,'depth/'),'tif');
irFiles = dirSmart(fullfile(SUN3Dpath,sequenceName,'ir/'),'tif');

depthFrameID = zeros(1,length(depthFiles));
depthTimestamp = zeros(1,length(depthFiles));
for i=1:length(depthFiles)
    id_time = sscanf(depthFiles(i).name, '%d-%ld.tif');
    depthFrameID(i) = id_time(1);
    depthTimestamp(i) = id_time(2);

    data.depthAll{i}= fullfile(fullfile(SUN3Dpath,sequenceName,'depth',depthFiles(i).name));
end

irFrameID = zeros(1,length(irFiles));
irTimestamp = zeros(1,length(irFiles));
for i=1:length(irFiles)
    id_time = sscanf(irFiles(i).name, '%d-%ld.tif');
    irFrameID(i) = id_time(1);
    irTimestamp(i) = id_time(2);

    data.irAll{i}= fullfile(fullfile(SUN3Dpath,sequenceName,'ir',irFiles(i).name));
end

data.imageTimestamp = imageTimestamp;
data.depthTimestamp = depthTimestamp;

data.imageTotalFrames = length(data.imageTimestamp);
data.depthTotalFrames = length(data.depthTimestamp);

% synchronize: find a depth for each image
frameCount = length(imageTimestamp);
IDimage2depth = zeros(1,frameCount);
for i=1:frameCount
    [~, IDimage2depth(i)]=min(abs(double(depthTimestamp)-double(imageTimestamp(i))));
end

if ~exist('frameIDs','var') || isempty(frameIDs)
    frameIDs = 1:frameCount;
end

data.sequenceName = sequenceName;

cnt = 0;
for frameID=frameIDs
    cnt = cnt + 1;
    data.depth{cnt} = fullfile(fullfile(SUN3Dpath,sequenceName,'depth',depthFiles(IDimage2depth(frameID)).name));
    data.ir{cnt} = fullfile(fullfile(SUN3Dpath,sequenceName,'ir',irFiles(IDimage2depth(frameID)).name));
end

kinectID = strsplit(sequenceName,'_');
kinectID = kinectID{end};

data.camera = load(fullfile(SUN3Dpath,'intrinsics',[kinectID '.mat']));

end

function files = dirSmart(page, tag)
    [files, status] = urldir(page, tag);
    if status == 0
        files = dir(fullfile(page, ['*.' tag]));
    end
end

function [files, status] = urldir(page, tag)
    if nargin == 1
        tag = '/';
    else
        tag = lower(tag);
        if strcmp(tag, 'dir')
            tag = '/';
        end
        if strcmp(tag, 'img')
            tag = 'jpg';
        end
    end
    nl = length(tag);
    nfiles = 0;
    files = [];

    % Read page
    page = strrep(page, '\', '/');
    [webpage, status] = urlread(page);

    if status
        % Parse page
        j1 = findstr(lower(webpage), '<a href="');
        j2 = findstr(lower(webpage), '');
        Nelements = length(j1);
        if Nelements>0
            for f = 1:Nelements
                % get HREF element
                chain = webpage(j1(f):j2(f));
                jc = findstr(lower(chain), '">');
                chain = deblank(chain(10:jc(1)-1));

                % check if it is the right type
                if length(chain)>length(tag)-1
                    if strcmp(chain(end-nl+1:end), tag)
                        nfiles = nfiles+1;
                        chain = strrep(chain, '%20', ' '); % replace space character
                        files(nfiles).name = chain;
                        files(nfiles).bytes = 1;
                    end
                end
            end
        end
    end
end

function XYZcamera = depth2XYZcamera(K, depth)
    sz = size(depth);
    [x,y] = meshgrid(1:sz(2), 1:sz(1));
    XYZcamera(:,:,1) = (x-K(1,3)).*depth/K(1,1);
    XYZcamera(:,:,2) = (y-K(2,3)).*depth/K(2,2);
    XYZcamera(:,:,3) = depth;
    XYZcamera(:,:,4) = depth~=0;
end

author: jianxiongxiao

repo: ProfXkit


text:
function [S, v, f, ORIGIN] = read_vf(fname, patch)
% function [S, v, f, ORIGIN] = read_vf(fname, patch)
% This function reads a matfile with vertices/faces-information in it.
% if [fname '.mat'] cannot be found, a '.geo'- or 'asc'-file is searched
% converted and saved as a matfile.
% S(1) = nr of surfaces
% S(2) = nr of vertices
% S(3) = nr of faces
% v: n x 3 matrix of vertices
% f: m x 4 matrix of faces
% ORIGIN: coordinates (2D or 3D) of ORIGIN in mm space, where the coordinates
% [0 0 0] in vertices space are mapped to.

% fname: filename of mat-file to be opened
% patch: flag: if 1, try to get patch-information, i.e. indices of vertices and faces
% in generating original surface, if 0, we don't need this information

%
% read_vf.m
%
% Original Author: Bruce Fischl
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:12 $
%    $Revision: 1.3 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

fn = spm_str_manip(fname, 'r');

source = 0;

if ~exist([fn '.mat'])
   % Try to open asc-file
   if patch == 0 | ~exist([fn '.asc'], 'file')
      % Try to open geo-file
      if patch == 1
         % no mat- or asc-file
         error(sprintf('Cannot find patch-infomation, file %s not found', [fn '.asc']))
      elseif ~exist([fn '.geo'], 'file')
         error(sprintf('Cannot find %s or convert a geo- or asc-file to it.', [fn '.mat']))
      else
         [S, v, f] = read_moviebyu([fn '.geo']);
         source = 2;
      end
   else
      [S, v, fi] = read_asc([fn '.asc']);
      source = 1;
   end
else
   load([fn '.mat']);
   if patch == 1 & exist('fi') ~= 1
      % used asked for patch information, but mat-file didn't contain this information
      % so open asc-file (can happen, if asc-file wasn't present when mat-file
      % was generated
      if exist ([fn '.asc'], 'file')
         [S, v, fi] = read_asc([fn '.asc']);
         save([fn '.mat'], 'fi', 'v', '-append');
         source = 1;
      else
         error(sprintf('Cannot find %s', [fn '.asc']));
      end
   elseif patch == 0 & exist('f') ~= 1
      % no patch information needed, try to get plain f
      if exist ([fn '.geo'], 'file')
         [S, v, f] = read_moviebyu([fn '.geo']);
         save([fn '.mat'], 'f', '-append');
         source = 2;
      else
         error(sprintf('Cannot find %s', [fn '.geo']));
      end
   end

   if patch == 0 & size(v, 2) ~= 3
      v = v(:, 2:4);
   end
end

if source == 1 | source == 2
   % Assume that ORIGIN is [128 128 128]
   ORIGIN = [128 128 128];
end

if source == 1
   save([fn '.mat'], 'S', 'v', 'fi', 'ORIGIN');
elseif source == 2
   save([fn '.mat'], 'S', 'v', 'f', 'ORIGIN');
end

if patch == 1
   f = fi;
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function y=gaussmix(m1,sigma1,m2,sigma2,lambda)

n=length(m1);

for i=1:n
    if rand>lambda
y(i)=random('normal',m1(i),sigma1);
    else
y(i)=random('normal',m2(i),sigma2);
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function Y = sec(X)
    % Symbolic secant.

    % Convert inputs to SymExpression
    % X = SymExpression(X);

    % construct the operation string
    sstr = ['Sec[' X.s ']'];

    % create a new object with the evaluated string
    Y = SymExpression(sstr);
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function z=state2RuvRR(xTar,useHalfRange,xTx,xRx,M,includeW)
%%STATE2RUVRR Convert state vectors consisting of at least 3D position and
%             velocity in 3D space into local bistatic r-u-v coordinates
%             with non-relativistic range rate.
%
%INPUTS: xTar An xDimXN matrix of N target states consisting of 3D position
%             and velocity components in the order
%             xTar=[xPosition;xVelocity] and possibly other components,
%             which will be ignored.
% useHalfRange A boolean value specifying whether the bistatic range value
%             should be divided by two. This normally comes up when
%             operating in monostatic mode, so that the range reported is a
%             one-way range. The default if this parameter is not provided
%             is false.
%         xTx An xTxDimXN matrix of the states of the transmitters
%             consisting of stacked 3D position and velocity components.
%             Other components will be ignored. If this parameter is
%             omitted, the transmitters are assumed to be stationary at the
%             origin. If only a single vector is passed, then the
%             transmitter state is assumed the same for all of the target
%             states being converted.
%         xRx An xRxDimXN matrix of the states of the receivers consisting
%             of stacked 3D position and velocity components. Other
%             components will be ignored. If this parameter is omitted, the
%             receivers are assumed to be stationary at the origin. If only
%             a single vector is passed, then the receiver state is assumed
%             the same for all of the target states being converted.
%           M A 3X3XN hypermatrix of the rotation matrices to go from the
%             alignment of the global coordinate system to that at the
%             receiver. The z-axis of the local coordinate system of the
%             receiver is the pointing direction of the receiver. If
%             omitted, then it is assumed that the local coordinate system
%             is aligned with the global and M=eye(3) --the identity matrix
%             is used. If only a single 3X3 matrix is passed, then is is
%             assumed to be the same for all of the N conversions.
%    includeW An optional boolean value indicating whether a third
%             direction cosine component should be included. The u and v
%             direction cosines are two parts of a 3D unit vector.
%             Generally, one might assume that the target is in front of
%             the sensor, so the third component would be positive and is
%             not needed. However, the third component can be included if
%             ambiguity exists. The default if this parameter is omitted or
%             an empty matrix is passed is false.
%
%OUTPUTS: z A 4XN (or 5XN if includeW is true) matrix of the target states
%           in xTar converted into bistatic r-u-v(-w)  and bistatic range
%           rate coordinates. If useHalfRange=true, then the r component is
%           half the bistatic range and the range rate is correspondingly
%           halved.
%
%Details of the conversions are given in [1].
%
%REFERENCES:
%[1] D. F. Crouse, "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

N=size(xTar,2);

if(nargin<6||isempty(includeW))
   includeW=false;
end

if(nargin<5||isempty(M))
    M=repmat(eye(3),[1,1,N]);
elseif(size(M,3)==1&&N>1)
    M=repmat(M,[1,1,N]);
end

if(nargin<4||isempty(xRx))
    xRx=zeros(6,N);
elseif(size(xRx,2)==1&&N>1)
    xRx=repmat(xRx,[1,N]);
end

if(nargin<3||isempty(xTx))
    xTx=zeros(6,N);
elseif(size(xTx,2)==1&&N>1)
    xTx=repmat(xTx,[1,N]);
end

if(nargin<2||isempty(useHalfRange))
   useHalfRange=false;
end

%Allocate space for the converted states.
if(includeW)
    z=zeros(5,N);
else
    z=zeros(4,N);
end

%Convert the positions.
z(1:(end-1),:)=Cart2Ruv(xTar,useHalfRange,xTx,xRx,M,includeW);

%Compute the bistatic range rates.
z(end,:)=getRangeRate(xTar(1:6,:),useHalfRange,xTx,xRx);
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [v,H,R,azel,exc,stat,rtk]=ppp_res(post,x,rtk,obs,nav,sv,dr,exc)

global glc
stat=1; opt=rtk.opt; MAXSAT=glc.MAXSAT; VAR_GLO_IFB=0.6^2;%#ok
nobs=size(obs,1); nf=rtk.NF;
v=zeros(2*nobs*nf,1); H=zeros(2*nobs*nf,rtk.nx); var=zeros(2*nobs*nf,1);
azel=zeros(nobs,2); dants=zeros(3,1); dgrav=0;
obsi=zeros(64,1); frqi=zeros(64,1); ve=zeros(64,1);

for i=1:glc.MAXSAT
    for j=1:opt.nf
        rtk.sat(i).vsat(j)=0;
    end
end

% earth tide correction
rr=x(1:3)+dr;
pos=ecef2pos(rr);

nv=1;ne=0;
for i=1:nobs

    sat=obs(i).sat; lam=nav.lam(sat,:);lam_=[nav.lam(sat,:),nav.lam(sat+1,:)];
    if lam_(fix(j/2)+1)==0||lam(1)==0,continue;end

    % satellite information
    rs=sv(i).pos; dts=sv(i).dts; var_rs=sv(i).vars; vs=sv(i).vel; svh=sv(i).svh;

    % distance/light of sight/azimuth/elevation
    [r,LOS]=geodist(rs,rr); azel(i,:)=satazel(pos,LOS);
    if r<=0||azel(i,2)<opt.elmin,continue;end

    [sys,~]=satsys(sat);
    if sys==0||rtk.sat(sat).vs==0||satexclude(sat,var_rs,svh,opt)==0||exc(i)==1
        exc(i)=1; continue;
    end

    % tropospheric and ionospheric model
    [dtdx,dtrop,vart,stat_t]=model_trop(obs(i).time,pos,azel(i,:),rtk,x,nav);
    [diono,vari,stat_i]    =model_iono(obs(i).time,pos,azel(i,:),rtk,x,nav,sat);
    if stat_t==0||stat_i==0,continue;end

    % satellite and receiver antenna model
    if opt.posopt(1)==1
        dants=satantpcv(rs,rr,nav.pcvs(sat));
    end
    dantr=antmodel(sat,opt.pcvr,opt.antdel(1,:),azel(i,:),opt.posopt(2));

    % phase windup model
    [rtk.sat(sat).phw,stat_tmp]=model_phw(rtk.sol.time,sat,nav.pcvs(sat).type,opt.posopt(3),rs,rr,vs,rtk.sat(sat).phw);
    if stat_tmp==0,continue;end

    % gravitational delay correction
    if opt.posopt(7)==1
        dgrav=model_grav(sys,rr,rs);
    end

    % corrected phase and code measurements
    [L,P,Lc,Pc]=corr_meas(rtk,obs(i),nav,dantr,dants,rtk.sat(sat).phw);

    j=0;
    while j<2*nf
        dcb=0; bias=0;

        if opt.ionoopt==glc.IONOOPT_IFLC
            if rem(j,2)==0,y=Lc;else,y=Pc;end
            if y==0,j=j+1;continue;end
        else
            if rem(j,2)==0,y=L(fix(j/2)+1);else,y=P(fix(j/2)+1);end
            if y==0,j=j+1;continue;end
            if sys==glc.SYS_GLO,mm=2;else,mm=1;end
            if fix(j/2)==1,dcb=-nav.rbias(mm,1);end
        end

        if rem(j,2)==0,C_K1=-1;else,C_K1=1;end
        gama=(lam(fix(j/2)+1)/lam(1))^2;
        C=gama*ionmapf(pos,azel(i,:))*C_K1;

        H(nv,:)=zeros(1,rtk.nx);
        H(nv,1:3)=-LOS;

        % receiver clock
        dtr=0;
        if     sys==glc.SYS_GPS
            dtr=x(rtk.ic+1);
            H(nv,rtk.ic+1)=1;
        elseif sys==glc.SYS_GLO
            dtr=x(rtk.ic+1)+x(rtk.ic+2);
            H(nv,rtk.ic+1)=1; H(nv,rtk.ic+2)=1;
            % for GLONASS icb
            if opt.gloicb==glc.GLOICB_LNF
                if (nf==1&&(fix(j/2)==0&&rem(j,2)==1))||(nf==2&&(fix(j/2)==1&&rem(j,2)==1))
                    frq=get_glo_fcn(sat,nav);
                    dtr=dtr+frq*x(rtk.iicb+1);
                    H(nv,rtk.iicb+1)=frq;
                end
            elseif opt.gloicb==glc.GLOICB_QUAD
                if (nf==1&&(fix(j/2)==0&&rem(j,2)==1))||(nf==2&&(fix(j/2)==1&&rem(j,2)==1))
                    frq=get_glo_fcn(sat,nav);
                    dtr=dtr+frq*x(rtk.iicb+1);
                    H(nv,rtk.iicb+1)=frq;
                    dtr=dtr+frq^2*x(rtk.iicb+2);
                    H(nv,rtk.iicb+2)=frq^2;
                end
            end
        elseif sys==glc.SYS_GAL
            dtr=x(rtk.ic+1)+x(rtk.ic+3);
            H(nv,rtk.ic+1)=1; H(nv,rtk.ic+3)=1;
        elseif sys==glc.SYS_BDS
            dtr=x(rtk.ic+1)+x(rtk.ic+4);
            H(nv,rtk.ic+1)=1; H(nv,rtk.ic+4)=1;
        elseif sys==glc.SYS_QZS
            dtr=x(rtk.ic+1)+x(rtk.ic+5);
            H(nv,rtk.ic+1)=1; H(nv,rtk.ic+5)=1;
        end

        % troposphere
        if opt.tropopt==glc.TROPOPT_EST||opt.tropopt==glc.TROPOPT_ESTG
            H(nv,rtk.it+1)=dtdx(1);
            if opt.tropopt==glc.TROPOPT_ESTG
                H(nv,rtk.it+2)=dtdx(2);H(nv,rtk.it+3)=dtdx(3);
            end
        end

        % ionosphere
        if opt.ionoopt==glc.IONOOPT_EST
            if rtk.x(rtk.ii+sat)==0,j=j+1;continue;end
            H(nv,rtk.ii+sat)=C;
        end

        % L5-receiver-dcb
        if fix(j/2)==2&&rem(j,2)==1
            dcb=dcb+rtk.x(rtk.id+1);
            H(nv,rtk.id+1)=1;
        end

        % ambiguity
        if rem(j,2)==0
            bias=x(rtk.ib+fix(j/2)*MAXSAT+sat);
            if bias==0,j=j+1;continue;end
            H(nv,rtk.ib+fix(j/2)*MAXSAT+sat)=1;
        end

        % residual
        dtS=dts*glc.CLIGHT;
        v(nv,1)=y-(r+dtr-dtS+dtrop+C*diono+dcb+bias-dgrav);

        if rem(j,2)==0,rtk.sat(sat).resc(fix(j/2)+1)=v(nv,1);
        else          ,rtk.sat(sat).resp(fix(j/2)+1)=v(nv,1);
        end

        % variance
        var_rr=varerr_ppp(sys,azel(i,2),fix(j/2),rem(j,2),opt);
        var(nv,1)=var_rr+var_rs+vart+C^2*vari;
        %if sys==glc.SYS_GLO&&rem(j,2)==1,var(nv,1)=var(nv,1)+VAR_GLO_IFB;end

        % reject satellite by pre-fit residuals
        if post==0&&opt.maxinno>0&&abs(v(nv))>opt.maxinno
            exc(i)=1;
            rtk.sat(sat).rejc(rem(j,2)+1)=rtk.sat(sat).rejc(rem(j,2)+1)+1;
            j=j+1; continue;
        end

        % record large post-fit residuals
        if post~=0&&abs(v(nv))>sqrt(var(nv))*4
            obsi(ne+1)=i;frqi(ne+1)=j;ve(ne+1)=v(nv);
            ne=ne+1;
        end

        if rem(j,2)==0,rtk.sat(sat).vsat(fix(j/2)+1)=1;end

        nv=nv+1;
        j=j+1;
    end

end

if nv-1==0
    v=NaN;H=NaN;var=NaN;
elseif nv-1<2*nobs*nf
    v(nv:end)=[]; H(nv:end,:)=[]; var(nv:end)=[];
end
obsi(ne+1:end,:)=[];frqi(ne+1:end,:)=[];ve(ne+1:end,:)=[];

% reject satellite with large and max post-fit residual
if post~=0 && ne>0
    vmax=ve(1);maxobs=obsi(1);maxfrqi=frqi(1);rej=1; %#ok
    j=2;
    while j<=ne
        if abs(vmax)>=abs(ve(j))
            j=j+1; continue;
        end
        vmax=ve(j);maxobs=obsi(j);maxfrqi=frqi(j);rej=j; %#ok
        j=j+1;
    end
    sat=obs(maxobs).sat;
    exc(maxobs)=1;
    rtk.sat(sat).rejc(rem(j,2)+1)=rtk.sat(sat).rejc(rem(j,2)+1)+1;
    stat=0;
    ve(rej)=0; %#ok
end

% measurement noise matrix
if nv-1>0
    R=zeros(nv-1,nv-1);
    for i=1:nv-1
        for j=1:nv-1
            if i==j
                R(i,j)=var(i);
            else
                R(i,j)=0;
            end
        end
    end
else
    R=NaN;
end

if post==0
    stat=nv-1;
end

return

repo: GINav


text:
%%*******************************************************************
%% HSDHKMrhsfun: compute the right-hand side vector of the
%%               Schur complement equation for the  HKM direction.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*******************************************************************

function [rhs,EinvRc,hRd] = HSDHKMrhsfun(blk,At,par,X,Z,rp,Rd,sigmu,hRd,dX,dZ)

m = par.m;
if (nargin > 8)
    corrector = 1;
else
    corrector = 0;
    hRd = zeros(m+2,1);
end
hEinvRc = zeros(m+2,1);
EinvRc = cell(size(blk,1),1);
if length(sigmu)==1; sigmu = sigmu*ones(1,size(blk,1)); end
%%
for p = 1:size(blk,1)
    pblk = blk(p,:);
    n = sum(pblk{2});
    if strcmp(pblk{1},'l')
        if (corrector)
            Rq = dX{p}.*dZ{p};
        else
            Rq = sparse(n,1);
            tmp  = par.dd{p}.*Rd{p};
            tmp2 = mexMatvec(At{p},tmp,1);
            hRd = hRd + tmp2;
        end
        EinvRc{p} = sigmu(p)./Z{p}-X{p} -Rq./Z{p};
        tmp2 = mexMatvec(At{p,1},EinvRc{p},1);
        hEinvRc = hEinvRc + tmp2;
    elseif strcmp(pblk{1},'q')
        if (corrector)
            ff{p} = qops(pblk,1./par.gamz{p},Z{p},3); %#ok
            hdx = qops(pblk,par.gamz{p},ff{p},5,dX{p});
            hdz = qops(pblk,par.gamz{p},ff{p},6,dZ{p});
            hdxdz = Arrow(pblk,hdx,hdz);
            Rq = qops(pblk,par.gamz{p},ff{p},6,hdxdz);
        else
            Rq  = sparse(n,1);
            tmp = par.dd{p}.*Rd{p} ...
                + qops(pblk,qops(pblk,Rd{p},par.Zinv{p},1),X{p},3) ...
                + qops(pblk,qops(pblk,Rd{p},X{p},1),par.Zinv{p},3);
            tmp2 = mexMatvec(At{p,1},tmp,1);
            hRd = hRd + tmp2;
        end
        EinvRc{p} = sigmu(p)*par.Zinv{p}-X{p} -Rq;
        tmp2 = mexMatvec(At{p,1},EinvRc{p},1);
        hEinvRc = hEinvRc + tmp2;
    elseif strcmp(pblk{1},'s')
        if (corrector)
            Rq = Prod3(pblk,dX{p},dZ{p},par.Zinv{p},0);
            Rq = 0.5*(Rq+Rq');
        else
            Rq = sparse(n,n);
            tmp = Prod3(pblk,X{p},Rd{p},par.Zinv{p},0,par.nzlistAy{p});
            EinvRc{p} = 0.5*(tmp+tmp');
            tmp2 = AXfun(pblk,At(p,:),par.permA(p,:),EinvRc(p));
            hRd = hRd + tmp2;
        end
        EinvRc{p} = sigmu(p)*par.Zinv{p}-X{p}- Rq;
        tmp2 = AXfun(pblk,At(p,:),par.permA(p,:),EinvRc(p));
        hEinvRc = hEinvRc + tmp2;
    end
end
%%
rhs = rp + hRd - hEinvRc;
rhs(m+1) = rhs(m+1) + (par.mu/par.tau - par.kap);
if (corrector)
    rhs(m+1) = rhs(m+1) - par.dtau*par.dkap/par.tau;
end
%%*******************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function f = viewfig(obs,pred,flag)
%CREATEFIGURE1(RES_VERF1,YMATRIX1)
%  RES_VERF1:  vector of x data
%  YMATRIX1:  matrix of y data

%  Auto-generated by MATLAB on 10-Jun-2008 18:14:54
%load mytest
YMatrix1 = [obs pred];
x=1:length(obs);
n = length(obs);
% res_verf1=res_verf(:,1);
% YMatrix1 = res_verf(:,2:3);
% Create figure
f = figure('PaperSize',[20.98 29.68],'Position',[50 300 600 500],'visible',flag);
%f = figure('PaperSize',[20.98 29.68]);
% Create axes
axes1 = axes('Parent',f,'Position',[0.08043 0.6038 0.875013 0.3554]);
box('on');
hold('all');

% Create multiple lines using matrix input to plot
plot1 = plot(x,YMatrix1,'Parent',axes1);
set(plot1(1),'DisplayName','observed');
set(plot1(2),'DisplayName','predicted');
set(plot1(1),'color',[0.5 0.5 0.5],'LineWidth',2)
set(plot1(2),'color','b')

xlim(axes1,[0 n])
ylim(axes1,[0 max(max(YMatrix1))])
% Create xlabel
xlabel('Time (hr)');

% Create ylabel
ylabel('Median of discharge');
%legend(axes1,'Observed','Predicted')

%% Zoomed axis
x=4200:4500;

% Create axes
axes2 = axes('Parent',f,'Position',[0.08043 0.09201 0.4089 0.4118]);
% Uncomment the following line to preserve the X-limits of the axes
% xlim([4200 4500]);
box('on');
hold('all');

% Create multiple lines using matrix input to plot
plot2 = plot(x,YMatrix1(x,:),'Parent',axes2);
set(plot2(1),'color',[0.5 0.5 0.5],'LineWidth',2)
set(plot2(2),'color','b')
ylim(axes2,[0 max(max(YMatrix1(x,:)))])
x=4700:5000;
% Create axes
axes3 = axes('Parent',f,'Position',[0.5703 0.09201 0.3835 0.4133]);
% Uncomment the following line to preserve the X-limits of the axes
% xlim([4700 5000]);
box('on');
hold('all');

% Create multiple lines using matrix input to plot
plot3 = plot(x,YMatrix1(x,:),'Parent',axes3);
set(plot3(1),'color',[0.5 0.5 0.5],'LineWidth',2)
set(plot3(2),'color','b')
ylim(axes3,[0 max(max(YMatrix1(x,:)))])
% Create legend
legend(axes1,'show','Location','Best');
%print -dtiff -r900 picture1
% [pathstr, name, ext, versn] = fileparts(figName) ;
% switch  ext
%     case '.fig'
%         saveas(f,figName)
%     case '.tiff'
%       print(f,'-dtiff', '-r900', figName)
%     case '.emf'
%         print(f,'-dmeta', '-r2400', figName)
%     case '.eps'
%         print(f,'-depsc', '-r2400', figName)
%     otherwise
%         error('does not know the figure format')
% end
%delete(f)


text:
%SimuNF_nl.m
%tic
function []=timecourse_VF(x0,i)
para
tf=1.5;
[t1,y1]=ode_yu_VF(0,tf,[-x0(1),-x0(2),0,0,0,0,0,0]',dt);
[t2,y2]=ode_yu_VF(0,tf,[-x0(1),-x0(2),0,0,0,0,0,0]',dt);
[t3,y3]=ode_yu_VF(0,tf,[-x0(1),-x0(2),0,0,0,0,0,0]',dt);
[t4,y4]=ode_yu_VF(0,tf,[-x0(1),-x0(2),0,0,0,0,0,0]',dt);
[t5,y5]=ode_yu_VF(0,tf,[-x0(1),-x0(2),0,0,0,0,0,0]',dt);

%%
figure(1)
subplot(222)
%
% hold on
% draw_bg
% line([-0.03 -0.03],[-0.2 0.01], 'linewidth', 2, 'color', [0,0,0])
% line([ 0.03  0.03],[-0.2 0.01], 'linewidth', 2, 'color', [0,0,0])
plot(y1(:,1)+x0(1),y1(:,2)+x0(2), 'b-', ...
    y2(:,1)+x0(1),y2(:,2)+x0(2), 'b-', ...
    y3(:,1)+x0(1),y3(:,2)+x0(2), 'b-', ...
    y4(:,1)+x0(1),y4(:,2)+x0(2), 'b-', ...
    y5(:,1)+x0(1),y5(:,2)+x0(2), 'b-', ...
    'Linewidth', 1.2)

figure(2)
subplot(8,3,(i-1)*3+2)
plot(t1,sqrt(y1(:,3).^2+y1(:,4).^2), 'b-', ...
     t2,sqrt(y2(:,3).^2+y2(:,4).^2), 'b-', ...
     t3,sqrt(y3(:,3).^2+y3(:,4).^2), 'b-', ...
     t4,sqrt(y4(:,3).^2+y4(:,4).^2), 'b-', ...
     t5,sqrt(y5(:,3).^2+y5(:,4).^2), 'b-', ...
    'Linewidth', 1.2)

figure(1)
%hold on
%hold off
%axis equal
%axis([-0.075 0.075 -0.3 .075])

end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [Yl,Yh,Yscale] = dtwavexfm2(X,nlevels,biort,qshift);

% Function to perform a n-level DTCWT-2D decompostion on a 2D matrix X
%
% [Yl,Yh,Yscale] = dtwavexfm2(X,nlevels,biort,qshift);
%
%     X -> 2D real matrix/Image
%
%     nlevels -> No. of levels of wavelet decomposition
%
%     biort ->  'antonini'   => Antonini 9,7 tap filters.
%               'legall'     => LeGall 5,3 tap filters.
%               'near_sym_a' => Near-Symmetric 5,7 tap filters.
%               'near_sym_b' => Near-Symmetric 13,19 tap filters.
%
%     qshift -> 'qshift_06' => Quarter Sample Shift Orthogonal (Q-Shift) 10,10 tap filters,
%                              (only 6,6 non-zero taps).
%               'qshift_a' =>  Q-shift 10,10 tap filters,
%                              (with 10,10 non-zero taps, unlike qshift_06).
%               'qshift_b' => Q-Shift 14,14 tap filters.
%               'qshift_c' => Q-Shift 16,16 tap filters.
%               'qshift_d' => Q-Shift 18,18 tap filters.
%
%
%     Yl     -> The real lowpass image from the final level
%     Yh     -> A cell array containing the 6 complex highpass subimages for each level.
%     Yscale -> This is an OPTIONAL output argument, that is a cell array containing
%               real lowpass coefficients for every scale.
%
%
% Example: [Yl,Yh] = dtwavexfm2(X,3,'near_sym_b','qshift_b');
% performs a 3-level transform on the real image X using the 13,19-tap filters
% for level 1 and the Q-shift 14-tap filters for levels >= 2.
%
% Nick Kingsbury and Cian Shaffrey
% Cambridge University, Sept 2001

if isstr(biort) & isstr(qshift)		%Check if the inputs are strings
   biort_exist = exist([biort '.mat']);
   qshift_exist = exist([qshift '.mat']);
   if biort_exist == 2 & qshift_exist == 2;        		%Check to see if the inputs exist as .mat files
      load (biort);
      load (qshift);
   else
      error('Please enter the correct names of the Biorthogonal or Q-Shift Filters, see help DTWAVEXFM2 for details.');
   end
else
   error('Please enter the names of the Biorthogonal or Q-Shift Filters as shown in help DTWAVEXFM2.');
end

orginal_size = size(X);

if ndims(X) >= 3;
   error(sprintf('The entered image is %dx%dx%d, please enter each image slice separately.',orginal_size(1),orginal_size(2),orginal_size(3)));
end

% The next few lines of code check to see if the image is odd in size, if so an extra ...
% row/column will be added to the bottom/right of the image
initial_row_extend = 0;  %initialise
initial_col_extend = 0;
if any(rem(orginal_size(1),2)), %if sx(1) is not divisable by 2 then we need to extend X by adding a row at the bottom
   X = [X; X(end,:)];           %Any further extension will be done in due course.
   initial_row_extend = 1;
end
if any(rem(orginal_size(2),2)), 	%if sx(2) is not divisable by 2 then we need to extend X by adding a col to the left
   X = [X X(:,end)];          %Any further extension will be done in due course.
   initial_col_extend = 1;
end
extended_size = size(X);

if nlevels == 0, return; end

%initialise
Yh=cell(nlevels,1);
if nargout == 3
   Yscale=cell(nlevels,1);   %this is only required if the user specifies a third output component.
end

S = [];
sx = size(X);
if nlevels >= 1,

   % Do odd top-level filters on cols.
   Lo = colfilter(X,h0o).';
   Hi = colfilter(X,h1o).';

   % Do odd top-level filters on rows.
   LoLo = colfilter(Lo,h0o).';			% LoLo
   Yh{1} = zeros([size(LoLo)/2  6]);
   Yh{1}(:,:,[1 6]) = q2c(colfilter(Hi,h0o).');			% Horizontal pair
   Yh{1}(:,:,[3 4]) = q2c(colfilter(Lo,h1o).');			% Vertical pair
   Yh{1}(:,:,[2 5]) = q2c(colfilter(Hi,h1o).');	      % Diagonal pair
   S = [ size(LoLo) ;S];
   if nargout == 3
      Yscale{1} = LoLo;
   end
end

if nlevels >= 2;
   for level = 2:nlevels;
      [row_size col_size] = size(LoLo);
      if any(rem(row_size,4)),		% Extend by 2 rows if no. of rows of LoLo are divisable by 4;
         LoLo = [LoLo(1,:); LoLo; LoLo(end,:)];
      end
      if any(rem(col_size,4)),		% Extend by 2 cols if no. of cols of LoLo are divisable by 4;
         LoLo = [LoLo(:,1)  LoLo  LoLo(:,end)];
      end

      % Do even Qshift filters on rows.
      Lo = coldfilt(LoLo,h0b,h0a).';
      Hi = coldfilt(LoLo,h1b,h1a).';

      % Do even Qshift filters on columns.
      LoLo = coldfilt(Lo,h0b,h0a).';	%LoLo
      Yh{level} = zeros([size(LoLo)/2  6]);
      Yh{level}(:,:,[1 6]) = q2c(coldfilt(Hi,h0b,h0a).');	% Horizontal
      Yh{level}(:,:,[3 4]) = q2c(coldfilt(Lo,h1b,h1a).');	% Vertical
      Yh{level}(:,:,[2 5]) = q2c(coldfilt(Hi,h1b,h1a).');	% Diagonal
      S = [ size(LoLo) ;S];
      if nargout == 3
         Yscale{level} = LoLo;
      end
   end
end

Yl = LoLo;

if initial_row_extend == 1 & initial_col_extend == 1;
   warning(sprintf(' \r\r The image entered is now a %dx%d NOT a %dx%d \r The bottom row and rightmost column have been duplicated, prior to decomposition. \r\r ',...
      extended_size(1),extended_size(2),orginal_size(1),orginal_size(2)));
end

if initial_row_extend == 1 ;
   warning(sprintf(' \r\r The image entered is now a %dx%d NOT a %dx%d \r Row number %d has been duplicated, and added to the bottom of the image, prior to decomposition. \r\r',...
      extended_size(1),extended_size(2),orginal_size(1),orginal_size(2),orginal_size(1)));
end

if initial_col_extend == 1;
   warning(sprintf(' \r\r The image entered is now a %dx%d NOT a %dx%d \r Col number %d has been duplicated, and added to the right of the image, prior to decomposition. \r\r',...
      extended_size(1),extended_size(2),orginal_size(1),orginal_size(2),orginal_size(2)));
end
return

%==========================================================================================
%						**********  	INTERNAL FUNCTION    **********
%==========================================================================================

function z = q2c(y)

% function z = q2c(y)
% Convert from quads in y to complex numbers in z.

sy = size(y);
t1 = 1:2:sy(1); t2 = 1:2:sy(2);
j2 = sqrt([0.5 -0.5]);

% Arrange pixels from the corners of the quads into
% 2 subimages of alternate real and imag pixels.
%  a----b
%  |    |
%  |    |
%  c----d

% Combine (a,b) and (d,c) to form two complex subimages.
p = y(t1,t2)*j2(1) + y(t1,t2+1)*j2(2);     % p = (a + jb) / sqrt(2)
q = y(t1+1,t2+1)*j2(1) - y(t1+1,t2)*j2(2); % q = (d - jc) / sqrt(2)

% Form the 2 subbands in z.
z = cat(3,p-q,p+q);

return

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
classdef distance

    methods(Static)

        function result = pairwise_distances(X, distance)
            % Returns the pairwise distance matrix
            if (nargin < 2)
                distance = 'euclidean';
            end
            % Pair wise distance matrix (symmetric)
            result  =  squareform(pdist(X', distance));
        end

        function result = sqrd_l2_distances_cw(X)
            % The vectors are stored column wise.
            % Number of vectors = number of columns
            result = spx.commons.distance.pairwise_x_y(X);
        end

        function result = sqrd_l2_distances_rw(X)
            % The vectors are stored row wise.
            % Number of vectors = number of rows
            result = spx.commons.distance.pairwise_x_y(X');
        end

        function [ distance_matrix ] = pairwise_x_y( X, Y )
            % Computes pairwise distance squared between vectors in X and Y
            % c.f. pdist built-in function
            % We assume that vectors are column vectors
            % Dimension of space to which the vectors belong
            %N = size(X,1);
            % Number of vectors in X
            %nx = size(X,2);
            % Number of vectors in Y
            %ny = size(Y,2);
            % The result
            % distance_matrix(i,j) = distance (x_i, y_j)
            if nargin < 2
                Y = X;
            end
            distance_matrix = bsxfun(@plus,dot(X,X,1)',dot(Y,Y,1))-2*(X'*Y);
        end

    end

end

author: indigits

repo: sparse-plex


text:
[dummy,T] = lu(rand(5));
b = rand(5,1);
%fprintf('The columns should be equal:\n');
r = [solve_triu(T,b) T\b];
assert(all(abs(r(:,1) - r(:,2)) < 1e-10))
r = [solve_tril(T',b) T'\b];
assert(all(abs(r(:,1) - r(:,2)) < 1e-10))
fprintf('Verified that solve_triu and solve_tril results match backslash.\n');

d = 100;
niter = (20000/d)^2;
A = rand(d);
[dummy,T] = lu(A);
b = rand(d,1);
tic; for i = 1:niter T\b; end; t1=toc/niter;
tic; for i = 1:niter solve_triu(T,b); end; t2=toc/niter;
fprintf('backslash: \t%g\nsolve_triu: \t%g (%g times faster)\n',t1,t2,t1/t2);
% backslash is detecting triangularity as a preprocessing step, which doubles
% the time.
%fprintf('by flops, should be %g times faster\n',...
%    flops_solve(T,b)/flops_solve_tri(T,b));

niter = ceil(niter/d);
tic; for i = 1:niter inv(T); end; t1=toc/niter;
%I = eye(size(T));
%tic; for i = 1:niter solve_triu(T,I); end; t2=toc;
tic; for i = 1:niter inv_triu(T); end; t2=toc/niter;
fprintf('inv: \t%g\ninv_triu: \t%g (%g times faster)\n',t1,t2,t1/t2);
fprintf('by flops, should be %g times faster\n',...
    flops_inv(rows(T))/flops_solve_tri(T,eye(size(T))));

author: Cloud-CV

repo: object-proposals


text:
function [M_Ru1, M_Ru2, M_v1, M_v2] = vox2ras_rsolve(Vc_C, inPlaneRotation)
%%
%%	(USE vox2ras_rsolveAA INSTEAD OF THIS METHOD!)
%%
%% NAME
%%
%%     vox2ras_rsolve.m (vox2ras_r{otation matrix}solve)
%%
%% AUTHOR
%%
%%	Rudolph Pienaar
%%
%% SYNOPSIS
%%
%%     [M_Ru1 M_Ru2 M_v1 M_v2 ] = vox2ras_rsolve(Vc_C, inPlaneRotation)
%%
%% ARGUMENTS
%%
%%      Vc_C		in      column vector defining a direction cosine for
%%					a volume
%%	inPlaneRotation	in	scalar float defining the in plane rotation as
%%					read from the data meas.asc file
%%	M_Ru1		out	v1 rotated by f(inPlaneRotation)
%%	M_Ru2		out	v2 rotated by f(inPlaneRotation)
%%      M_v1		out     1st candidate vox2ras matrix (no rotation)
%%	M_v2		out	2nd candidate vox2ras matrix (no rotation)
%%
%%		where f() is a linear function on inPlaneRotation
%%
%% DESCRIPTION
%%
%%	"vox2ras_rsolve" attempts to find a candidate vox2ras rotation matrix that
%%	contains the vector C = [ci cj ck]' such that:
%%
%%				ai	bi	ci 	0
%%		vox2ras	=	-1	bj	cj	0
%%				ak	-1	ck	0
%%				 0	 0	 0	1
%%
%%	Two candidate matrices are returned due to the quadratic nature of the
%%	solutions, although for all practical purposes only the first solution
%%	is important.
%%
%%	This method originally used an empirical linear function to "rotate" the
%%	solution vox2ras matrices to a Siemens-based reference. Soon afterwards,
%%	a method was found that exactly determined the Siemens-based reference
%% 	orientation that depreciated the empirical linear function - and formed
%%	the basis of the vox2ras_rsolveAA function.
%%
%%	For completeness sake, this method was folded back into this function, but
%%	the vox2ras_rsolveAA is simpler and can address different slab orientations.
%%
%%
%% PRECONDITIONS
%%
%%	o The vector C is read from a Siemens meas.asc file such that
%%		ci	= sSliceArray.asSlice[0].sNormal.dSag
%%		cj	= sSliceArray.asSlice[0].sNormal.dCor
%%		ck	= sSliceArray.asSlice[0].sNormal.dTra
%%
%% POSTCONDITIONS
%%
%%	o All returned matrices are 4x4.
%%	o Only the rotations of the vox2ras matrix are determined by this function.
%%		The center of k-space is not determined.
%%	o M_v1 M_v2: 2 vox2ras matrices corresponding to two quadratic solutions.
%%	o M_Ru1 M_Ru2: v1 v2 rotated by f(inPlaneRotation).
%%
%% SEE ALSO
%%
%%	vox2ras_rsolveAA- determine the rotational component of a vox2ras matrix
%%				using Siemens reference orientations directly
%%	vox2ras_ksolve	- determine the k-space col in RAS of a vox2ras matrix
%%	vox2ras_dfmeas	- main function: determines the vox2ras matrix from a
%%			  Siemens meas.asc file.
%%
%% HISTORY
%%
%% 18 May 2004
%% o Initial design and coding.
%%
%% 26 May 2004
%% o Touching up...
%% o 4x4 return sizes fixed.
%%
%% 02 June 2004
%% o Changed return calling parameters
%% o Normalize only first two columns
%% o Incorporated Siemens reference orientation calculations
%%

ci	= Vc_C(1);
cj	= Vc_C(2);
ck	= Vc_C(3);
vox2ras	= zeros(3, 3);
v1	= zeros(3, 3);
v2	= zeros(3, 3);

%% We know that each column vector in vox2ras = [A B C] is orthogonal, i.e. the
%% dot products of each column vector is zero (A.B = B.C = A.C = 0). We also
%% know that C = B X A. Thus, we can expand the cross and dot products and solve
%% for ak first, arriving at the following quadratic equation:
%%
%%	(ci^2 + ck^2)ak^2 - (2cj ck + cicjck)ak + (ci+1)(cj^2 + ci^2) = 0

%
% vox2ras_rsolve.m
%
% Original Author: Rudolph Pienaar
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

a	= ci^2 + ck^2;
b	= -(2*cj*ck + ci*cj*ck);
c	= (ci+1)*(cj^2 + ci^2);

ak1	= (-b + sqrt(b^2-4*a*c))/2/a;
ak2	= (-b - sqrt(b^2-4*a*c))/2/a;

ai1	= (cj - ak1*ck)/ci;
ai2	= (cj - ak2*ck)/ci;

bj1	= (ci + 1) / ak1;
bi1	= (ck - cj*bj1) / ci;
bj2	= (ci + 1) / ak2;
bi2	= (ck - cj*bj2) / ci;

M3_v1 = [
	ai1	bi1	ci
	-1	bj1	cj
	ak1	-1	ck
];

M3_v2 = [
	ai2	bi2	ci
	-1	bj2	cj
	ak2	-1	ck
];

for col=1:2,
	M3_v1(:,col) = M3_v1(:,col) ./ norm(M3_v1(:,col));
end

for col=1:2,
	M3_v2(:,col) = M3_v2(:,col) ./ norm(M3_v2(:,col));
end

%% First candidate
theta_c1	= atan(-ak1);
theta_f1	= theta_c1 - inPlaneRotation;
M3_Mu1		= [	 cos(theta_f1)	 sin(theta_f1)	0
			-sin(theta_f1)	 cos(theta_f1)	0
			 	0		0	1];
M3_Ru1		= M3_v1 * M3_Mu1;

%% Second candidate
theta_c2	= atan(-ak2);
theta_f2	= theta_c2 - inPlaneRotation;
M3_Mu2		= [	 cos(theta_f2)	 sin(theta_f2)	0
			-sin(theta_f2)	 cos(theta_f2)	0
			 	0		0	1];
M3_Ru2		= M3_v2 * M3_Mu2;

%% ***********************************************************
%% DEPRECIATED - Using theta_c1 improves accuracy
%% ***********************************************************
%% The above calculated rotation matrices define an (x,y) plane
%%	given by the first two column vectors. Since we have
%%	"hardcoded" two components of these vectors, we have
%%	fixed them along a default orientation. This default
%%	orientation needs to be rotated by an additional theta_f
%%	radians in order to arrive at the final vox2ras matrix.
%%	By studying several existing vox2ras matrices and their
%%	meas.asc InPlaneRotation values, a simple linear relationship
%%	between this theta_f and the InPlaneRotation value was
%%	found:
%%
%%		theta_f	= m * InPlaneRotation + b
%%
%%	where
%%		m = -1.0025
%%		b = -0.5188 (roughly equal to pi/6 = )
%%

%  m		= -1.0025;
%  b		= -0.5188;
%  theta_f		= m*inPlaneRotation + b;
%  %theta_f		= inPlaneRotation;
%
%  M3_Mu	= [	 cos(theta_f)	 sin(theta_f)	0
%  		-sin(theta_f)	 cos(theta_f)	0
%  		 	0		0	1];
%
%  M3_Ru1	= M3_v1 * M3_Mu;
%  M3_Ru2	= M3_v2 * M3_Mu;
%% ***********************************************************
%%                                                 DEPRECIATED
%% ***********************************************************

M_v1	= eye(4);	M_v1(1:3, 1:3)	= M3_v1;
M_v2	= eye(4);	M_v2(1:3, 1:3)	= M3_v2;
M_Ru1	= eye(4);	M_Ru1(1:3, 1:3)	= M3_Ru1;
M_Ru2	= eye(4);	M_Ru2(1:3, 1:3)	= M3_Ru2;

%% All done!

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function X = gpsample_vbmc(vp,gp,Ns,origflag)
%GPSAMPLE_VBMC Sample from GP obtained through VBMC.

if nargin < 4 || isempty(origflag); origflag = true; end

D = size(gp.X,2);

if isfield(gp,'s2') && ~isempty(gp.s2)
    % Evaluate GP input length scale (use geometric mean)
    Ns_gp = numel(gp.post);
    ln_ell = zeros(D,Ns_gp);
    for s = 1:Ns_gp; ln_ell(:,s) = gp.post(s).hyp(1:D); end
    gplengthscale = exp(mean(ln_ell,2))';
    X_rescaled = bsxfun(@rdivide,gp.X,gplengthscale); % Rescaled GP training inputs

    % Evaluate GP observation noise on training inputs
    sn2new = zeros(size(gp.X,1),Ns_gp);
    for s = 1:Ns_gp
        hyp_noise = gp.post(s).hyp(gp.Ncov+1:gp.Ncov+gp.Nnoise); % Get noise hyperparameters
        if isfield(gp,'s2')
            s2 = gp.s2;
        else
            s2 = [];
        end
        % s2 = noiseshaping_vbmc(s2,gp.y,options);
        sn2new(:,s) = gplite_noisefun(hyp_noise,gp.X,gp.noisefun,gp.y,s2);
    end
    sn2new = mean(sn2new,2);

    % Estimate observation noise variance over variational posterior
    xx = vbmc_rnd(vp,2e4,0,0);
    [~,pos] = min(sq_dist(bsxfun(@rdivide,xx,gplengthscale),X_rescaled),[],2);
    sn2_avg = mean(sn2new(pos));    % Use nearest neighbor approximation
else
    sn2_avg = 0;
end

VarThresh = max(1,sn2_avg);

W = 2*(D+1);
x0 = vbmc_rnd(vp,W,0,0);
X = gplite_sample(gp,Ns,x0,'parallel',[],[],VarThresh);
if origflag
    X = warpvars_vbmc(X,'inv',vp.trinfo);
end

end

%SQ_DIST Compute matrix of all pairwise squared distances between two sets
% of vectors, stored in the columns of the two matrices, a (of size n-by-D)
% and b (of size m-by-D).
function C = sq_dist(a,b)

n = size(a,1);
m = size(b,1);
mu = (m/(n+m))*mean(b,1) + (n/(n+m))*mean(a,1);
a = bsxfun(@minus,a,mu); b = bsxfun(@minus,b,mu);
C = bsxfun(@plus,sum(a.*a,2),bsxfun(@minus,sum(b.*b,2)',2*a*b'));
C = max(C,0);

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
% Test Projection onto L1 Ball
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     27/06/2017  Royi Avital
%   *   First release.

%% General Parameters

run('InitScript.m');

numRows     = 5;
ballRadius  = 3;
stopThr     = 1e-6;

%% Generating Data

vY = 10 * rand([numRows, 1]) - 5;

%% Solution by CVX

cvx_begin('quiet')
    % cvx_precision('best');
    variable vXCvx(numRows)
    minimize( norm(vXCvx - vY) )
    subject to
        norm(vXCvx ,1) <= ballRadius;
cvx_end

disp([' ']);
disp(['CVX Solution Summary']);
disp(['The CVX Solver Status - ', cvx_status]);
disp(['The Optimal Value Is Given By - ', num2str(cvx_optval)]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vXCvx.'), ' ]']);
disp([' ']);

%% Solution by Dual Function and Newton Iteration

% vX = ProjectL1Ball(vY, ballRadius, stopThr);
vX = ProjectL1BallExact(vY, ballRadius);

disp([' ']);
disp(['Dual Function Solution Summary']);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);

%% Display Results

disp([' ']);
disp(['CVX Solution L1 Norm - ', num2str(norm(vXCvx, 1))]);
disp(['Dual Function Solution L1 Norm - ', num2str(norm(vX, 1))]);
disp(['Solutions Difference L1 Norm - ', num2str(norm(vXCvx - vX, 1))]);
disp([' ']);

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
%%%%%%%%%%%%%%%%%%%%%%%
% Compute features for a set of video files from datasets
%
close all;
clear;
warning('off','all');
addpath(genpath('./lib/vqm/TLVQM_naive_release'));

%%
% parameters
algo_name = 'TLVQM_naive';
data_name = 'KONVID_1K';

%% *You need to customize here*
if strcmp(data_name, 'TEST_VIDEOS')
    data_path = 'videos'; % dataset video path
elseif strcmp(data_name, 'KONVID_1K')
    data_path = '/media/ztu/Seagate-ztu-ugc/KONVID_1K/KoNViD_1k_videos';
elseif strcmp(data_name, 'LIVE_VQC')
    data_path = '/media/ztu/Seagate-ztu-ugc/LIVE_VQC/VideoDatabase';
elseif strcmp(data_name, 'YOUTUBE_UGC')
    data_path = '/media/ztu/Seagate-ztu-ugc/YT_UGC/original_videos';
elseif strcmp(data_name, 'LIVE_VQA')
    data_path = '/media/ztu/Seagate-ztu/LIVE_VQA/videos';
end

video_tmp = './video_tmp';
if ~exist(video_tmp, 'dir'), mkdir(video_tmp); end
feat_path = './mos_feat_files';
filelist_csv = fullfile(feat_path, [data_name,'_metadata.csv']);
filelist = readtable(filelist_csv);
num_videos = size(filelist,1);
if ~exist(feat_path, 'dir'), mkdir(feat_path); end
out_feat_name = fullfile(feat_path, [data_name,'_',algo_name,'_feats.mat']);
feats_mat = zeros( num_videos, 75 );
%===================================================

tic
for i = 1:num_videos
% get video full path and decoded video name
	if strcmp(data_name, 'TEST_VIDEOS')
	    video_name = fullfile(data_path,  filelist.video_name{i});
	    yuv_name = fullfile(video_tmp, [filelist.video_name{i}, '.yuv']);
	elseif strcmp(data_name, 'KONVID_1K')
	    video_name = fullfile(data_path,  [num2str(filelist.flickr_id(i)),'.mp4']);
	    yuv_name = fullfile(video_tmp, [num2str(filelist.flickr_id(i)), '.yuv']);
	elseif strcmp(data_name, 'LIVE_VQC')
	    video_name = fullfile(data_path, filelist.File{i});
	    yuv_name = fullfile(video_tmp, [filelist.File{i}, '.yuv']);
	elseif strcmp(data_name, 'YOUTUBE_UGC')
	    video_name = fullfile(data_path, filelist.category{i},...
		[num2str(filelist.resolution(i)),'P'],[filelist.vid{i},'.mkv']);
	    yuv_name = fullfile(video_tmp, [filelist.vid{i}, '.yuv']);
	elseif strcmp(data_name, 'LIVE_VQA')
	    strs = strsplit(filelist.filename{i}, '_');
	    video_name = fullfile(data_path, [strs{1}(1:2), '_Folder'], filelist.filename{i});
	    yuv_name = video_name;
	end
	fprintf('\n---\nComputing features for %d-th sequence: %s\n', i, video_name);
    % decode video and store in video_tmp dir
    if ~strcmp(video_name, yuv_name)
    cmd = ['ffmpeg -loglevel error -y -i ', video_name, ' -pix_fmt yuv420p -vsync 0 ', yuv_name];
    system(cmd);
    end

    % get video meta data
    width = filelist.width(i);
    height = filelist.height(i);
    framerate = round(filelist.framerate(i));
    nb_frames = filelist.nb_frames(i);

    % calculate video features
    tic
    feats_mat(i,:) = compute_nrvqa_features(yuv_name, [width, height], framerate);
    toc
    delete(yuv_name)
end
toc
save(out_feat_name, 'feats_mat');

author: vztu

repo: VIDEVAL

save_path: github-repos/MATLAB/vztu-VIDEVAL


text:
function [keypointsAll, imgidxs_missing]= det2keypointsImgMulticutBBoxMulti(expidx,parts,annolist_bbox,bDisjoint)

fprintf('det2keypointsImgMulticutBBoxMulti()\n');

if (nargin < 4)
    bDisjoint = false;
end

p = exp_params(expidx);

multicutDir = p.multicutDir;
cidxs = p.cidxs;
keypointsAll = [];

if (isfield(p,'nPeopleMax'))
    nPeopleMax = p.nPeopleMax;
else
    nPeopleMax = inf;
end

if (isfield(p,'nPeopleMin'))
    nPeopleMin = p.nPeopleMin;
else
    nPeopleMin = 0;
end

load(p.testGT,'annolist');

assert(length(annolist) == length(annolist_bbox));

nmissing = 0;
bVis = false;
imgidxs_missing = [];

for imgidx = 1:length(annolist_bbox)

    fprintf('.');
    fname = [multicutDir '/imgidx_' padZeros(num2str(imgidx),4) '_cidx_' num2str(cidxs(1)) '_' num2str(cidxs(end))];

    try
        assert(length(annolist(imgidx).annorect) >= nPeopleMin && ...
            length(annolist(imgidx).annorect) <= nPeopleMax);
        load(fname,'unLab','unPos','unProb');
    catch
        for i = 1:length(p.pidxs)
            pidx = p.pidxs(i);
            jidx = parts(pidx+1).pos(1);
            keypointsAll(imgidx).det{1}{jidx+1} = [];
        end
        imgidxs_missing = [imgidxs_missing; imgidx];
        nmissing = nmissing + 1;
        if (~mod(imgidx,100))
            fprintf(' %d/%d\n',imgidx,length(annolist_bbox));
        end
        continue;
    end

    for ridx = 1:length(annolist_bbox(imgidx).annorect)

        bbox = annolist_bbox(imgidx).annorect(ridx).bbox;
        if (~isempty(bbox))
            idxs = find((unPos(:,1) >= bbox(1) & unPos(:,2) >= bbox(2) & ...
                unPos(:,1) <= bbox(3) & unPos(:,2) <= bbox(4)));

            for i = 1:length(p.pidxs)
                pidx = p.pidxs(i);
                assert(parts(pidx+1).pos(1) == parts(pidx+1).pos(2));
                jidx = parts(pidx+1).pos(1);
                det = [];
                id = find(cidxs == i);
                if (~isempty(id))
                    [val,idx] = max(unProb(idxs,id));
                    det = [unPos(idxs(idx),:) val];
                    if (bDisjoint)
                        unProb(idxs(idx),:) = -inf;
                    end
                end
                keypointsAll(imgidx).det{ridx}{jidx+1} = det;
            end
        else
            for i = 1:length(p.pidxs)
                pidx = p.pidxs(i);
                assert(parts(pidx+1).pos(1) == parts(pidx+1).pos(2));
                jidx = parts(pidx+1).pos(1);
                keypointsAll(imgidx).det{ridx}{jidx+1} = [];
            end
        end
    end

    if (bVis)
        colors = {'r','g','b','c','m','y'};
        markers = {'+','o','s'}; %,'x','.','-','s','d'
        img = imread(annolist_bbox(imgidx).image.name);
        figure(100); clf; imagesc(img); axis equal; hold on;
        for j = 1:length(keypointsAll(imgidx).det)
            labels = cell(0);
            n = 0;
            points = zeros(1,2);
            for i = 1:length(p.pidxs)
                pidx = p.pidxs(i);
                jidx = parts(pidx+1).pos(1);
                pp = keypointsAll(imgidx).det{j}{jidx+1};
                if (~isempty(pp))
                    n = n  + 1;
                    labels{n} = num2str(i);
%                     labels{n} = sprintf('%1.2f',pp(3));
                    points(n,:) = pp(1:2);
                end
            end
            lp = j;
            if (lp <= 5)
                m = markers{1};
            elseif (lp > 5 && lp <= 11)
                m = markers{2};
            else
                m = markers{3};
            end
            plot(points(:,1),points(:,2),[colors{mod(lp,6)+1} m],'MarkerSize',5);
            text(points(:,1)+10,points(:,2),labels,'FontSize',6,'BackgroundColor','w');
            axis off;
        end
    end

    if (~mod(imgidx, 100))
        fprintf(' %d/%d\n',imgidx,length(annolist_bbox));
    end
end
fprintf(' done\n');
fprintf('nmissing: %d\n',nmissing);

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [u,p,info,Ai,Bi,invMt,Pro_u,Pro_p] = mgstokesRT0(A,B,f,g,u,p,node,elem,ufreeDof,option,varargin)
%% MGSTOKES
%
% Created by Ming Wang and Jie Zhou based on discussion with Long Chen.
% Improvement by Long Chen: improve efficiency, include smoother and add
% bisection case and BDM1B.
%
% See also: asmgstokes, mgstokes, mg, mgMaxwell
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

tic;

% Assign default values to unspecified parameters
if ~exist('option','var')
    option = [];
end
inputoption = option; % record the input one
option = mgoptions(option,length(f)+length(g));    % parameters
% specific choices for mgstokesRT0
if ~isfield(inputoption,'smoothingstep')  % smoothing steps
    option.smoothingstep = 2;
end
tol = option.tol; maxIt = option.solvermaxit;
smoothingstep = option.smoothingstep; solver = option.solver;
printlevel = option.printlevel; setupflag = option.setupflag;

if setupflag == true
%% Reconstruced hierarchical meshes and transfer operators
N0 = 8;
N = size(node,1);
% NT = size(elem,1);
level = 20;
NL(level+1) = N; % now NL(1:level) = 0;
nodei = cell(level,1);
elemi = cell(level,1);
Pro_u = cell(level,1);
Pro_p = cell(level,1);
nodei{level} = node;
elemi{level} = elem;
for j = level: -1 : 2
    switch option.refType
        case 'red'
            [elemi{j-1},newHB] = uniformcoarsenred(elemi{j}); % coasen red refinement
            if ~isempty(newHB)
                Pro_u{j-1} = transferedgered(elemi{j-1},elemi{j}); % transfer operator of u
                Pro_p{j-1} = repmat(speye(size(elemi{j-1},1)),4,1);   % transfer operator of p
            end
        case 'bisect'
            % merge two coarsen of bisection grids s.t. the ratio is 1/4.
            [tempelem,newHB1,tree] = uniformcoarsen(elemi{j}); % coarse bisection
            if ~isempty(newHB1)
                % first coarsen
                tempPro_u = transferedgecoarsen(tempelem,elemi{j},tree);
                tempPro_p = transferelem(tempelem,elemi{j},tree);
                % second coarsen
                [elemi{j-1},newHB2,tree] = uniformcoarsen(tempelem); % coarse bisection
                Pro_u{j-1} = transferedgecoarsen(elemi{j-1},tempelem,tree);
                Pro_u{j-1} = tempPro_u*Pro_u{j-1};
                Pro_p{j-1} = transferelem(elemi{j-1},tempelem,tree);
                Pro_p{j-1} = tempPro_p*Pro_p{j-1};
                newHB = [newHB1; newHB2];
            end
    end
    if (isempty(newHB)) || (size(elemi{j-1},1)< 2*N0)
    % no nodes are removed or it reaches the coarsest level
        NL = NL(j:end);
        break;
    end
    NL(j) = NL(j+1) - size(newHB,1); % update NL(k)
    nodei{j-1} = nodei{j};           % node in i-th level
    nodei{j-1}(newHB(:,1),:) = [];   % remove fine nodes
end
level = length(NL)-1;    % actual level
nodei = nodei(end-level+1: end);
elemi = elemi(end-level+1: end);
Pro_u = Pro_u(end-level+1: end);
Pro_p = Pro_p(end-level+1: end);

%% No coarsening
if level == 1  % no coarsening. use exact solver
    Np = size(p,1);
    Nu = size(f,1);
    bigA = [A B'; B sparse(Np,Np)];
    bigF = [f; g-mean(g)];
    bigu = zeros(Nu+Np,1);
    bigu(1:end-1) = bigA(1:end-1,1:end-1)\bigF(1:end-1);
    u = bigu(1:Nu);
    p = bigu(Nu+1:end);
    flag = 2; itStep = 0; err = norm(bigF-bigA*bigu)/norm(bigF); time = toc;
    info = struct('solverTime',time,'itStep',itStep,'error',err,'flag',flag,'stopErr',max(err(end,:)));
    fprintf('direct solver time = %4.3g s,\n',time);
    Ai=[]; Bi=[]; invMt =[]; Pro_u =[]; Pro_p =[];
    return
end

%% Matrices in each level
Ai = cell(level,1);
Ai{level} = A;
Bi = cell(level,1);
Bi{level} = B;
invMt = cell(level,1);
area = simplexvolume(node,elem);
NTf = size(elem,1);
invMt{level} = spdiags(1./area,0,NTf,NTf);
NE = size(Pro_u{level-1},1);
ufreeDofj = ufreeDof(ufreeDof <= NE);
for j = level:-1:2
    % get mesh and data in the j-1 level
    elem = sortelem(elemi{j-1});
    node = nodei{j-1};
    [elem2edge,edge] = dofedge(elem);
    [Dlambda,area,elemSign] = gradbasis(node,elem);
    NTj = size(elem,1);

    % assemble matrix using the mesh information.
    % Mv: Lumped mass matrix for vertex: P1 element
    Mv = accumarray([elem(:,1);elem(:,2);elem(:,3)],[area;area;area]/3,[N,1]);
    Mv((Mv == 0)) = [];
    Nj = length(Mv);
    invMv = spdiags(1./Mv,0,Nj,Nj);
    Mv = spdiags(Mv,0,Nj,Nj);
    % Me: Mass matrix for RT0 element
    Me = getmassmatvec(elem2edge,area,Dlambda,'RT0');
    %invMt: the inverse of Mass matrix for P0 element
    invMt{j-1} = spdiags(1./area,0,NTj,NTj);
    % B: - divergence operator
    Bj = -icdmat(double(elem2edge),elemSign*[1 -1 1]);
    % C: curl operator
    C = icdmat(double(edge),[-1 1]);
    % R: weak rot operator
    R = invMv*C'*Me;
    % Vector Laplacian
    Ai{j-1} =  Bj'*invMt{j-1}*Bj + R'*Mv*R;

    % find free dof in each level
    isFixedDof = false(size(Ai{j-1},1),1);
    isDirichlet = false(size(Ai{j-1},1),1);
    % add coarsening into bdFlag
    bdFlag = setboundary(node,elem,'Dirichlet');
    if ~isempty(bdFlag)
        isDirichlet(elem2edge(bdFlag(:)==1)) = true;
        isFixedDof(isDirichlet) = true;
        ufreeDofc = find(~isFixedDof);
    end

    % modify transfer operators
    Pro_u{j-1} = Pro_u{j-1}(ufreeDofj,ufreeDofc);
    Ai{j-1} = Ai{j-1}(ufreeDofc,ufreeDofc);
    Bi{j-1} = Bj(:,ufreeDofc);
    ufreeDofj = ufreeDofc;
end
clear nodei elemi
end % end for setup

%% Only set up the transfer operators
if strcmp(solver,'NO')
    flag = 0; itStep = 0; err = 0; time = toc;
    info = struct('solverTime',time,'itStep',itStep,'error',err,'flag',flag,'stopErr',max(err(end,:)));
    return
end

%% No need of set up
if setupflag == false
    Ai = varargin{1};
    Bi = varargin{2};
   invMt = varargin{3};
   Pro_u = varargin{4};
   Pro_p = varargin{5};
   level = size(Ai,1);
end

%% Construct other necessary operators
Nu = zeros(level,1);
Nu(level) = size(A,1);
Np = zeros(level,1);
Np(level) = size(B,1);
bigAi = cell(level,1);
bigAi{level} = [Ai{level} Bi{level}'; Bi{level} sparse(Np(level),Np(level))];
Res_u = cell(level,1);
Res_p = cell(level,1);
for j = level:-1:2
    Res_u{j} = Pro_u{j-1}';
    Res_p{j} = Pro_p{j-1}';
    Nu(j-1) = size(Ai{j-1},1);
    Np(j-1) = size(Bi{j-1},1);
    bigAi{j-1} = [Ai{j-1} Bi{j-1}'; Bi{j-1} sparse(Np(j-1),Np(j-1))];
end
Ndof = Nu + Np;

%% Smoothers
Su = cell(level,1);
Sp = cell(level,1);
for k = 1:level
    Sp{k} = tril(Bi{k}*(Bi{k})');
    Su{k} = tril(Ai{k});
end

%% Solver
bigF = [f; g-mean(g)];
bigu = [u(ufreeDof); p];
bigr = bigF - bigAi{level}*bigu;
nb = norm(bigF);
err = zeros(maxIt,1);
err(1) = norm(bigr)/nb;
k = 1;
% condest(bigAi{J}(1:end-1,1:end-1))
if strcmp(solver,'GMRES') % GMRES solver
    restart = 50; prefunc = @(r)kcycle(r,level);
    tic;
    [bige, itStep, relRes_u] = PFGMRES(bigAi{level}, bigr, bigu, maxIt, restart, tol, prefunc,0);
    bigu = bigu + bige;
    time = toc;
    fprintf('itStep=%d, relRes_u = %10.6g, time = %4.3g s,\n', ...
             itStep, relRes_u(itStep), time);
else  % MG
    while (max(err(k)) > tol) && (k <= maxIt)
        k = k + 1;
        % one multigrid cycle
        switch (solver)
            case 'VCYCLE'
                biguerr = vcycle(bigr);
            case 'WCYCLE'
                biguerr = wcycle(bigr);
        end
        % update solution
        bigu = bigu + biguerr;
        % compute residual
        bigr = bigr - bigAi{level}*biguerr;
        % compute the relative errror
        err(k) = norm(bigr)/nb;
        if printlevel >= 2
            fprintf('#dof: %8.0u, MG %8s iter: %2.0u, err = %8.4e\n',...
                 Ndof(level), solver, k-1, err(k));
        end
    end
end
err = err(1:k);
itStep = k-1;
u = bigu(1:Nu(level));
p = bigu(Nu(level)+1:end);

%% Output
if k > maxIt
    flag = 1;
else
    flag = 0;
end
time = toc;
if printlevel >= 1
    fprintf('#dof: %6.0u,  #nnz: %6.0u, level: %2.0u  MG %6s iter: %2.0u,  err = %8.4e,  time = %4.2g s\n',...
             Ndof(level), nnz(bigAi{level}), level, solver, itStep, err(end), time)
end
if (flag == 1) && (printlevel>0)
   fprintf('NOTE: the iterative method does not converge! \n');
end
info = struct('solverTime',time,'itStep',itStep,'error',err,'flag',flag,'stopErr',max(err(end,:)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions vcycle, wcycle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Vcycle MG
    function e = vcycle(r,J)
        if nargin < 2, J = level; end
        if J == 1 % exact solver in the coaresest grid
           e = zeros(size(r));
           e(1:end-1) = bigAi{J}(1:end-1,1:end-1)\r(1:end-1);
           e(Nu(J)+1:end) = e(Nu(J)+1:end)-mean(e(Nu(J)+1:end));
           return
        end
        ru = r(1:Nu(J));
        rp = r(Nu(J)+1:end);

        % pre-smoothing in the fine grid
        [eu,ep] = DGSRT0(Ai{J},Bi{J},Bi{J}',zeros(Nu(J),1),zeros(Np(J),1),...
                         ru,rp,smoothingstep,Su{J},Sp{J},invMt{J});

        % form residual and restrict onto coarse grid
        ruc = Res_u{J}*(ru - Ai{J}*eu-(Bi{J})'*ep);
        rpc = Res_p{J}*(rp - Bi{J}*eu);

        % coarse grid correction twice
        rc = [ruc;rpc];
        ec = vcycle(rc,J-1);

        % prolongate coarse grid correction to the fine grid
        tempeu = Pro_u{J-1}*ec(1:Nu(J-1));
        tempep = Pro_p{J-1}*ec(Nu(J-1)+1:end);
        eu = eu + tempeu;
        ep = ep + tempep;

        % post-smoothing in the fine grid
        [eu,ep] = DGSRT0(Ai{J},Bi{J},Bi{J}',eu,ep,ru,rp,smoothingstep,...
                         Su{J},Sp{J},invMt{J});
        e = [eu;ep];
    end

%% Wcycle MG
    function e = wcycle(r,J)
        if nargin < 2, J = level; end
        if J == 1 % exact solver in the coaresest grid
           e = zeros(size(r));
           e(1:end-1) = bigAi{J}(1:end-1,1:end-1)\r(1:end-1);
           e(Nu(J)+1:end) = e(Nu(J)+1:end)-mean(e(Nu(J)+1:end));
           return
        end
        ru = r(1:Nu(J));
        rp = r(Nu(J)+1:end);

        % pre-smoothing in the fine grid
        [eu,ep] = DGSRT0(Ai{J},Bi{J},Bi{J}',zeros(Nu(J),1),zeros(Np(J),1),...
                         ru,rp,smoothingstep,Su{J},Sp{J},invMt{J});

        % form residual and restrict onto coarse grid
        ruc = Res_u{J}*(ru - Ai{J}*eu-(Bi{J})'*ep);
        rpc = Res_p{J}*(rp - Bi{J}*eu);

        % coarse grid correction twice
        rc = [ruc;rpc];
        ec = wcycle(rc,J-1);
        ec = ec + wcycle(rc-bigAi{J-1}*ec,J-1);

        % prolongate coarse grid correction to the fine grid
        tempeu = Pro_u{J-1}*ec(1:Nu(J-1));
        tempep = Pro_p{J-1}*ec(Nu(J-1)+1:end);
        eu = eu + tempeu;
        ep = ep + tempep;

        % post-smoothing in the fine grid
        [eu,ep] = DGSRT0(Ai{J},Bi{J},Bi{J}',eu,ep,ru,rp,smoothingstep,...
                         Su{J},Sp{J},invMt{J});
        e = [eu;ep];
    end

%% DGS for RT0
    function [u,p] = DGSRT0(A,B,Bt,u,p,f,g,itStep,Su,Sp,invMp)
        for s = 1: itStep
            % Step 1: relax Momentum eqns by G-S
            for i = 1:2
                u = u + Su\(f-Bt*p-A*u);
            end
            % Step 2: relax Continuity eqns by G-S
            rp = g - B*u;
            dq = Sp\rp;
            % Step 3: update u and p
            u = u + Bt*dq;
            p = p - invMp*rp;
%             p = p - Sp\(B*(A*(Bt*dq)));
        end
    end
end

repo: ifem


text:
% superline() - dashed (moving) line between two points
%
% Usage:
%   >> [handles] = superline( X, Y, color, linewidth, offset, middle );
%
% Inputs:
%   X          - abscices of the points (2 abscices for the 2 points)
%   Y          - ordinates of the points (2 ordinates for the 2 points)
%   color      - color of the line (default black)
%   linewidth  - thickness of the line (same as line 'linewidth' property, default:1)
%   offset     - change the offset of the dash (in between 0 and 1) (default:0)
%   middle     - put a small line at the middle
%
% Outputs:
%   handles    - handles of all the objects composing the line

% arno@salk.edu, Arnaud Delorme, CNL / Salk Institute, 2001

% This program is free software; you can redistribute it and/or
% modify it.
% 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.

function h = superline(X, Y, color, thickness, offset, middle);

if nargin < 2
	error('Not enough arguments. Type help superline');
	return;
end;
if nargin < 3
	color = 'k';
end;
if nargin < 4
	thickness = 1;
end;
if nargin < 5
	offset = 0;
end;
if nargin < 6
	middle = 0;
end;

% ratio X/Y
% ---------
if abs(Y(1) - Y(2)) == 0
	ratioXY = 1;
else
	ratioXY = abs(X(1) - X(2))/abs(Y(1) - Y(2));
end;

PADSTRIPE  = 0.8;
SIZESTRIPE = 0.2;
%PADSTRIPE  = 0.18;
%SIZESTRIPE = 0.09;

% original line vector
% --------------------
if ratioXY >= 0.05
	slope = (Y(2)-Y(1)) / (X(2)-X(1));
	slopenormalize = [1 slope]/sqrt(slope^2+1); % (1,slope) are the coordinate of the non unitary vector
else
	slopenormalize = sign(Y(2)-Y(1))*[0 1];
end;
% reverse vector if necessary
if X(1) > X(2)
	slopenormalize = - slopenormalize;
end;
vectpad    = slopenormalize*PADSTRIPE*abs(complex(X(1), Y(1)) - complex(X(2), Y(2)));
vectstripe = slopenormalize*SIZESTRIPE*abs(complex(X(1), Y(1)) - complex(X(2), Y(2)));
nblines    = 1/(PADSTRIPE+SIZESTRIPE);

% draw the discontinuous line
% ---------------------------
cont = 1;
compter = 1;
firstpass = 1;
posx1 = X(1);
posy1 = Y(1);
offsetshifted = mod( 0.5 + offset + SIZESTRIPE/2,1);

while cont == 1
	if firstpass == 1
		if offsetshifted/nblines < SIZESTRIPE
			posx1 = posx1 + vectstripe(1)*(offsetshifted/nblines)/SIZESTRIPE;
			posy1 = posy1 + vectstripe(2)*(offsetshifted/nblines)/SIZESTRIPE;
			posx2 = posx1 + vectpad(1);
			posy2 = posy1 + vectpad(2);
		else
			posx2 = posx1 + vectstripe(1)*(offsetshifted/nblines-SIZESTRIPE)/SIZESTRIPE; %consider the offset from 0 to 1
			posy2 = posy1 + vectstripe(2)*(offsetshifted/nblines-SIZESTRIPE)/SIZESTRIPE; %consider the offset from 0 to 1
		end;
		firstpass = 0;
		%posy1
		%posy2
	else
		posx2 = posx1 + vectpad(1);
		posy2 = posy1 + vectpad(2);
		%posy1
		%posy2
	end;
	%fprintf(['Abscice1 ' num2str(X(1)) '  Abscice2 ' num2str(X(1)) '  equal ' num2str(X(1) == X(2)) '\n']);
	if ratioXY < 0.1
		if Y(1) > Y(2)
			if posy2 <= Y(2)
				posx2 = X(2);
				posy2 = Y(2);
				cont = 0;
			end;
		else
			if posy2 >= Y(2)
				posx2 = X(2);
				posy2 = Y(2);
				cont = 0;
			end;
		end;
	else
		if X(1) > X(2)
			if posx2 <= X(2)
				posx2 = X(2);
				posy2 = Y(2);
				cont = 0;
			end;
		else
			if posx2 >= X(2)
				posx2 = X(2);
				posy2 = Y(2);
				cont = 0;
			end;
		end;
	end;
	h(compter) = line([posx1 posx2], [posy1 posy2]);
	set(h(compter), 'color', color);
	set(h(compter), 'linewidth', thickness);
	posx1 = posx2+vectstripe(1);
	posy1 = posy2+vectstripe(2);
	compter = compter +1;
	if ratioXY < 0.1
		if Y(1) > Y(2)
			if posy1 <= Y(2)
				cont = 0;
			end;
		else
			if posy1 >= Y(2)
				cont = 0;
			end;
		end;
	else
		if X(1) > X(2)
			if posx1 <= X(2)
				cont = 0;
			end;
		else
			if posx1 >= X(2)
				cont = 0;
			end;
		end;
	end;
end;

%patch(patchccordinatesX, patchccordinatesY, [0 0 0 0], color);

if middle
	hold on; h = plot((X(1)+X(2))/2, (Y(1)+Y(2))/2, '*k', 'markersize', thickness/2);
end;

return

% tests
figure;
for i=9:31
	fprintf('******************************* %d \n', mod(i/10, 1));
	superline( [ 0+i 1+i ] , [ 10 20 ], 'b', 5, mod(i/10, 1),1);
end;

author: PatternRecognition

repo: OpenBMI


text:
% Test if the A B Matrices are valid

% Coefficients
mdlCoeff

A = [ 0 1 0 0; [-ks -bs ks bs]/mb ; ...
      0 0 0 1; [ks bs -ks-kt -bs]/mw];
B = [0; 10000/mb ; 0;  -10000/mw];

susp_sys_linear(rand(4,1),rand, A, B);

P_init = lyap(A',eye(4));
K1 = 1/r*B'*P_init;

[P,~,K] = care(A,B,eye(4),1)

% Optimal values
%
% P =
%
%     1.8639    0.0703   -1.5871    0.0123
%     0.0703    0.0638   -0.9956    0.0069
%    -1.5871   -0.9956   40.3489   -0.0763
%     0.0123    0.0069   -0.0763    0.0063
%
%
% K =
%
%     0.2868    0.9727  -20.4658   -0.8259


text:
function x = half_normal_cdf_inv ( cdf, a, b )

%*****************************************************************************80
%
%% HALF_NORMAL_CDF_INV inverts the Half Normal CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%    0.0 <= CDF <= 1.0.
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < B.
%
%    Output, real X, the corresponding argument.
%
  if ( cdf < 0.0 | 1.0 < cdf )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HALF_NORMAL_CDF_INV - Fatal error!\n' );
    fprintf ( 1, '  CDF < 0 or 1 < CDF.\n' );
    error ( 'HALF_NORMAL_CDF_INV - Fatal error!' );
  end

  cdf2 = 0.5 * ( cdf + 1.0 );

  x = normal_cdf_inv ( cdf2, a, b );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%--------------------------------------------------------------------------
% This function takes the groups resulted from spectral clutsering and the
% ground truth to compute the misclassification rate.
% groups: [grp1,grp2,grp3] for three different forms of Spectral Clustering
% s: ground truth vector
% Missrate: 3x1 vector with misclassification rates of three forms of
% spectral clustering
%--------------------------------------------------------------------------
% Copyright @ Ehsan Elhamifar, 2012
%--------------------------------------------------------------------------

function Missrate = Misclassification(groups,s)

n = max(s);
for i = 1:size(groups,2)
    Missrate(i,1) = missclassGroups( groups(:,i),s,n ) ./ length(s);
end

repo: BUPT-ICS-Courseware


text:
function node_num = grid_t3_node_num ( nelemx, nelemy )

%*****************************************************************************80
%
%% GRID_T3_NODE_NUM counts the nodes in a grid of T3 elements.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NELEMX, NELEMY, the number of elements along the
%    X and Y directions.  The number of elements generated will be
%    NELEMX * NELEMY.
%
%    Output, integer NODE_NUM, the number of nodes in the grid.
%
  node_num = ( nelemx + 1 ) * ( nelemy + 1 );

  return
end


text:
function [assignment, cost] = assignmentsuboptimal1(distMatrix)

%ASSIGNMENTSUBOPTIMAL1    Compute suboptimal assignment
%		ASSIGNMENTSUBOPTIMAL1(DISTMATRIX) computes a suboptimal assignment for the
%		given rectangular distance (or weight) matrix, for example the assignment
%		of tracks (in rows) to observations (in columns). The result is a column
%		vector containing the assigned column number in each row (or 0 if no
%		assignment could be done).
%
%		[ASSIGNMENT, COST] = ASSIGNMENTSUBOPTIMAL1(DISTMATRIX) returns the
%		assignment vector and the overall cost.
%
%		The algorithm is designed for distance matrices with many forbidden and
%		singly validated assignments, (rows or columns containing only one finite
%		element). The algorithm first searches the matrix for singly validated
%		columns and rejects all assignments with multiply validated row.
%		Afterwards, singly validated rows are searched and assignments to
%		multiply validated columns are rejected. Then, for each row that
%		validates only with singly validated columns (and the other way around),
%		the minimum elements is chosen and the assignment is made. If there are
%		still assignments open, the minimum element in the distMatrix is
%		searched and the corresponding assignment is made.
%
%		In scenarios without any forbidden assignments, the algorithm reduceds
%		to the last step is will give the same result as ASSIGNMENTOPTIMAL2. If
%		there are only some assignments forbidden, the algorithm will perform
%		poorly because singly validated assignments are preferred.
%
%		The last step can still be optimized, see the comments in
%		ASSIGNMENTOPTIMAL2.
%
%		Written by Markus Buehren, www.Lss.uni-stuttgart.de
%		Last modified 14.12.2004

% initialize
[nOfRows, nOfColumns] = size(distMatrix);
nOfValidObservations  = zeros(nOfRows,1);
nOfValidTracks        = zeros(1,nOfColumns);
assignment            = zeros(nOfRows,1);
cost                  = 0;

% compute number of validations for each track
for row=1:nOfRows
	nOfValidObservations(row) = length(find(isfinite(distMatrix(row,:))));
end

if any(nOfValidObservations < nOfColumns)

	if all(nOfValidObservations == 0)
		return
	end

	repeatSteps = 1;
	while repeatSteps

		repeatSteps = 0;

		% step 1: reject assignments of multiply validated tracks to singly validated observations
		for col=1:nOfColumns
			index = isfinite(distMatrix(:,col));
			nOfValidTracks(col) = length(find(index));
			if any(nOfValidObservations(index) == 1)
				index = index & (nOfValidObservations > 1);
				if any(index)
					distMatrix(index, col)      = inf;
					nOfValidObservations(index) = nOfValidObservations(index) - 1;
					nOfValidTracks(col)         = nOfValidTracks(col) - length(find(index));
					repeatSteps = 1;
				end
			end
		end

		% step 2: reject assignments of multiply validated observations to singly validated tracks
		if nOfColumns > 1
			for row=1:nOfRows
				index = isfinite(distMatrix(row,:));
				if any(nOfValidTracks(index) == 1)
					index = index & (nOfValidTracks > 1);
					if any(index)
						distMatrix(row, index)    = inf;
						nOfValidTracks(index)     = nOfValidTracks(index) - 1;
						nOfValidObservations(row) = nOfValidObservations(row) - length(find(index));
						repeatSteps = 1;
					end
				end
			end
		end

	end % while repeatSteps
	%disp(sprintf('xx = %d', xx));

	% for each multiply validated track that validates only with singly validated
	% observations, choose the observation with minimum distance
	for row=1:nOfRows
		if nOfValidObservations(row) > 1
			index = isfinite(distMatrix(row,:));
			if all(nOfValidTracks(index) == 1)
				[minDist, col] = min(distMatrix(row,:));
				assignment(row)    = col;
				cost               = cost + minDist;
				distMatrix(row,:)  = inf;
				distMatrix(:,col)  = inf;
			end
		end
	end

	% for each multiply validated observation that validates only with singly validated
	% track, choose the track with minimum distance
	for col=1:nOfColumns
		if nOfValidTracks(col) > 1
			index = isfinite(distMatrix(:,col));
			if all(nOfValidObservations(index) == 1)
				[minDist, row] = min(distMatrix(:,col));
				assignment(row)    = col;
				cost               = cost + minDist;
				distMatrix(row,:)  = inf;
				distMatrix(:,col)  = inf;
			end
		end
	end

end

% now, recursively search for the minimum element and do the assignment
while 1

	% find minimum distance observation-to-track pair
	[minDist, index1] = min(distMatrix, [], 1);
	[minDist, index2] = min(minDist);
	row = index1(index2);
	col = index2;

	if isfinite(minDist)

		% make the assignment
		assignment(row)    = col;
		cost               = cost + minDist;
		distMatrix(row, :) = inf;
		distMatrix(:, col) = inf;

	else
		break
	end

end

author: kristinbranson

repo: JAABA


text:
function [C, D] = WatsonSHCoeff(k)
% function [C, D] = WatsonSHCoeff(k)
% Computes the spherical harmonic (SH) coefficients of the Watson's
% distribution with the concentration parameter k (kappa) up to the 12th order
% and the derivatives if requested.
%
% Truncating at the 12th order gives good approximation for kappa up to 64.
%
% INPUTS:
%
% k should be an array of positive numbers, specifying a set of
% concentration parameters for the Watson's distribution.
%
% OUTPUTS:
%
% C will be a 2-D array and each row contains the SH coefficients of the
% orders 0, 2, 4, ..., to 2n for the parameter in the corresponding row in
% k.
%
% Note that the SH coefficients of the odd orders are always zero.
%
% D will be the 1st order derivative of C.
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

large = find(k>30);
exact = find(k>0.1);
approx = find(k<=0.1);
% Necessary to make matlab happy when k is a single value
exact = exact(:);
approx = approx(:);
large = large(:);

% The maximum order of SH coefficients (2n)
n = 6;

% Computing the SH coefficients
C = zeros(length(k),n+1);

% 0th order is a constant
C(:,1) = 2*sqrt(pi);

% Precompute the special function values
sk = sqrt(k(exact));
sk2 = sk.*k(exact);
sk3 = sk2.*k(exact);
sk4 = sk3.*k(exact);
sk5 = sk4.*k(exact);
sk6 = sk5.*k(exact);
sk7 = sk6.*k(exact);
k2 = k.^2;
k3 = k2.*k;
k4 = k3.*k;
k5 = k4.*k;
k6 = k5.*k;
k7 = k6.*k;

erfik = NODDI_erfi(sk);
ierfik = 1./erfik;
ek = exp(k(exact));
dawsonk = 0.5*sqrt(pi)*erfik./ek;

% for large enough kappa
C(exact,2) = 3*sk - (3 + 2*k(exact)).*dawsonk;
C(exact,2) = sqrt(5)*C(exact,2).*ek;
C(exact,2) = C(exact,2).*ierfik./k(exact);

C(exact,3) = (105 + 60*k(exact) + 12*k2(exact)).*dawsonk;
C(exact,3) = C(exact,3) -105*sk + 10*sk2;
C(exact,3) = .375*C(exact,3).*ek./k2(exact);
C(exact,3) = C(exact,3).*ierfik;

C(exact,4) = -3465 - 1890*k(exact) - 420*k2(exact) - 40*k3(exact);
C(exact,4) = C(exact,4).*dawsonk;
C(exact,4) = C(exact,4) + 3465*sk - 420*sk2 + 84*sk3;
C(exact,4) = C(exact,4)*sqrt(13*pi)/64./k3(exact);
C(exact,4) = C(exact,4)./dawsonk;

C(exact,5) = 675675 + 360360*k(exact) + 83160*k2(exact) + 10080*k3(exact) + 560*k4(exact);
C(exact,5) = C(exact,5).*dawsonk;
C(exact,5) = C(exact,5) - 675675*sk + 90090*sk2 - 23100*sk3 + 744*sk4;
C(exact,5) = sqrt(17)*C(exact,5).*ek;
C(exact,5) = C(exact,5)/512./k4(exact);
C(exact,5) = C(exact,5).*ierfik;

C(exact,6) = -43648605 - 22972950*k(exact) - 5405400*k2(exact) - 720720*k3(exact) - 55440*k4(exact) - 2016*k5(exact);
C(exact,6) = C(exact,6).*dawsonk;
C(exact,6) = C(exact,6) + 43648605*sk - 6126120*sk2 + 1729728*sk3 - 82368*sk4 + 5104*sk5;
C(exact,6) = sqrt(21*pi)*C(exact,6)/4096./k5(exact);
C(exact,6) = C(exact,6)./dawsonk;

C(exact,7) = 7027425405 + 3666482820*k(exact) + 872972100*k2(exact) + 122522400*k3(exact)  + 10810800*k4(exact) + 576576*k5(exact) + 14784*k6(exact);
C(exact,7) = C(exact,7).*dawsonk;
C(exact,7) = C(exact,7) - 7027425405*sk + 1018467450*sk2 - 302630328*sk3 + 17153136*sk4 - 1553552*sk5 + 25376*sk6;
C(exact,7) = 5*C(exact,7).*ek;
C(exact,7) = C(exact,7)/16384./k6(exact);
C(exact,7) = C(exact,7).*ierfik;

% for very large kappa
if size(large,1) > 0
  lnkd = log(k(large)) - log(30);
  lnkd2 = lnkd.*lnkd;
  lnkd3 = lnkd2.*lnkd;
  lnkd4 = lnkd3.*lnkd;
  lnkd5 = lnkd4.*lnkd;
  lnkd6 = lnkd5.*lnkd;
  C(large,2) = 7.52308 + 0.411538*lnkd - 0.214588*lnkd2 + 0.0784091*lnkd3 - 0.023981*lnkd4 + 0.00731537*lnkd5 - 0.0026467*lnkd6;
  C(large,3) = 8.93718 + 1.62147*lnkd - 0.733421*lnkd2 + 0.191568*lnkd3 - 0.0202906*lnkd4 - 0.00779095*lnkd5 + 0.00574847*lnkd6;
  C(large,4) = 8.87905 + 3.35689*lnkd - 1.15935*lnkd2 + 0.0673053*lnkd3 + 0.121857*lnkd4 - 0.066642*lnkd5 + 0.0180215*lnkd6;
  C(large,5) = 7.84352 + 5.03178*lnkd - 1.0193*lnkd2 - 0.426362*lnkd3 + 0.328816*lnkd4 - 0.0688176*lnkd5 - 0.0229398*lnkd6;
  C(large,6) = 6.30113 + 6.09914*lnkd - 0.16088*lnkd2 - 1.05578*lnkd3 + 0.338069*lnkd4 + 0.0937157*lnkd5 - 0.106935*lnkd6;
  C(large,7) = 4.65678 + 6.30069*lnkd + 1.13754*lnkd2 - 1.38393*lnkd3 - 0.0134758*lnkd4 + 0.331686*lnkd5 - 0.105954*lnkd6;
end

% for small kappa
C(approx,2) = 4/3*k(approx) + 8/63*k2(approx);
C(approx,2) = C(approx,2)*sqrt(pi/5);

C(approx,3) = 8/21*k2(approx) + 32/693*k3(approx);
C(approx,3) = C(approx,3)*(sqrt(pi)*0.2);

C(approx,4) = 16/693*k3(approx) + 32/10395*k4(approx);
C(approx,4) = C(approx,4)*sqrt(pi/13);

C(approx,5) = 32/19305*k4(approx);
C(approx,5) = C(approx,5)*sqrt(pi/17);

C(approx,6) = 64*sqrt(pi/21)*k5(approx)/692835;

C(approx,7) = 128*sqrt(pi)*k6(approx)/152108775;

if nargout == 1
	return;
end

% Computing the derivatives
dawsonk2 = dawsonk.^2;
idawsonk2 = 1./dawsonk2;

D = zeros(length(k),n+1);
D(:,1) = 0.0;

% exact
D(exact,2) = -k(exact) + (2*sk2 -sk).*dawsonk + 2*dawsonk2;
D(exact,2) = (.75*sqrt(5*pi))*D(exact,2)./k2(exact).*idawsonk2;

D(exact,3) = 21*k(exact) - 2*k2(exact);
D(exact,3) = D(exact,3) + (63*sk -44*sk2 + 4*sk3).*dawsonk;
D(exact,3) = D(exact,3) - (84 + 24*k(exact)).*dawsonk2;
D(exact,3) = D(exact,3)*(15*sqrt(pi)/32)./k3(exact).*idawsonk2;

D(exact,4) = -165*k(exact) + 20*k2(exact) - 4*k3(exact);
D(exact,4) = D(exact,4) + (-825*sk + 390*sk2 - 44*sk3 + 8*sk4).*dawsonk;
D(exact,4) = D(exact,4) + (990 + 360*k(exact) + 40*k2(exact)).*dawsonk2;
D(exact,4) = D(exact,4)*(21*sqrt(13*pi)/128)./k4(exact).*idawsonk2;

D(exact,5) = 225225*k(exact) - 30030*k2(exact) + 7700*k3(exact) - 248*k4(exact);
D(exact,5) = D(exact,5) + (1576575*sk - 600600*sk2 + 83160*sk3 - 15648*sk4 + 496*sk5).*dawsonk;
D(exact,5) = D(exact,5) - (1801800 + 720720*k(exact) + 110880*k2(exact) + 6720*k3(exact)).*dawsonk2;
D(exact,5) = D(exact,5)*(3*sqrt(17*pi)/2048)./k5(exact).*idawsonk2;

D(exact,6) = -3968055*k(exact) + 556920*k2(exact) - 157248*k3(exact) + 7488*k4(exact) - 464*k5(exact);
D(exact,6) = D(exact,6) + (-35712495*sk + 11834550*sk2 - 1900090*sk3 + 336960*sk4 - 15440*sk5 + 928*sk6).*dawsonk;
D(exact,6) = D(exact,6) + (39680550 + 16707600*k(exact) + 2948400*k2(exact) + 262080*k3(exact) + 10080*k4(exact)).*dawsonk2;
D(exact,6) = D(exact,6)*(11*sqrt(21*pi)/8192)./k6(exact).*idawsonk2;

D(exact,7) = 540571185*k(exact) - 78343650*k2(exact) + 23279256*k3(exact) - 1319472*k4(exact) + 119504*k5(exact) - 1952*k6(exact);
D(exact,7) = D(exact,7) + (5946283035*sk - 1786235220*sk2 + 319642092*sk3 - 53155872*sk4 + 2997456*sk5 - 240960*sk6 + 3904*sk7).*dawsonk;
D(exact,7) = D(exact,7) - (6486854220 + 2820371400*k(exact) + 537213600*k2(exact) + 56548800*k3(exact) + 3326400*k4(exact) + 88704*k5(exact)).*dawsonk2;
D(exact,7) = D(exact,7)*(65*sqrt(pi)/65536)./k7(exact).*idawsonk2;

% approximation
D(approx,2) = 4/3 + 16/63*k(approx) - 16/315*k2(approx) - 128/6237*k3(approx);
D(approx,2) = D(approx,2)*sqrt(pi/5);

D(approx,3) = 16/105*k(approx) + 32/1155*k2(approx) - 3712/675675*k3(approx) - 5888/2837835*k4(approx);
D(approx,3) = D(approx,3)*sqrt(pi);

D(approx,4) = 16/231*k2(approx) + 128/10395*k3(approx) - 256/106029*k4(approx);
D(approx,4) = D(approx,4)*sqrt(pi/13);

D(approx,5) = 128/19305*k3(approx) + 256/220077*k4(approx);
D(approx,5) = D(approx,5)*sqrt(pi/17);

D(approx,6) = 64/138567*k4(approx);
D(approx,6) = D(approx,6)*sqrt(pi/21);

D(approx,7) = 256/50702925*k5(approx);
D(approx,7) = D(approx,7)*sqrt(pi);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function a = i4mat_elim ( m, n, a )

%*****************************************************************************80
%
%% I4MAT_ELIM carries out exact Gauss elimination on an I4MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in A.
%
%    Input, integer N, the number of columns in A.
%
%    Input, integer A(M,N), the matrix to be Gauss eliminated.
%
%    Output, integer A(M,N), the Gauss-eliminated matrix.
%

%
%  Initialize the swap parity counter.
%
  iswap = 1;
%
%  For each column JCOL...
%
  for jcol = 1 : min ( m, n )
%
%  Find the maximum element in rows JCOL through M.
%
    amax = abs ( a(jcol,jcol) );
    imax = jcol;

    for i = jcol + 1 : m
      if ( amax < abs ( a(i,jcol) ) )
        amax = abs ( a(i,jcol) );
        imax = i;
      end
    end
%
%  If the maximum entry is nonzero, then...
%
    if ( amax ~= 0 )
%
%  If the maximum entry does not occur in row JCOL, then swap rows.
%
      if ( imax ~= jcol )
        iswap = -iswap;
        temp(1:n) = a(jcol,1:n);
        a(jcol,1:n) = a(imax,1:n);
        a(imax,1:n) = temp(1:n);
      end
%
%  Eliminate all nonzero entries in column JCOL, below the diagonal entry.
%
      for i = jcol+1 : m

        if ( a(i,jcol) ~= 0 )

          jmult = a(i,jcol);
          imult = a(jcol,jcol);
          ifact = i4_gcd ( imult, jmult );
          imult = imult / ifact;
          jmult = jmult / ifact;

          for j = jcol : n
            a(i,j) = jmult * a(jcol,j) - imult * a(i,j);
          end

        end

      end
%
%  Remove any row or column factors.
%
      [ a, irow, icol ] = i4mat_red ( m, n, a );

    end

  end

  return
end


text:
%[ P, inls ] = ht_simple_ransac_p3p( u, X, rthr, maxiter )
%u: 3 x n image points
%X: 3 x n 3D points
%rthr: inlier threshold
%maxiter: default=1000

function [ P, inls ] = ht_lo_ransac_p3p( u, X, rthr, max_iter )
if nargin < 4
    max_iter = 1000;
end

%initialization
u = bsxfun(@rdivide, u, sqrt(sum(u.^2, 1)));
Npts = size(u, 2);
rthr = cos(rthr);
max_inlsnum = 3;
no_iter = 0;
P = [];
inls = false(1, Npts);
%ransac
while no_iter < max_iter
    no_iter = no_iter + 1;

    idx = randperm(Npts, 3);
    P_cand = P3PSolver([u(:, idx); X(:, idx)]);
    [inls_cand, inls_cand_num] = calculate_inls_angular(P_cand, u, X, rthr);
    if length(P_cand) == 0
%         no_iter = no_iter - 1;
        continue;
    elseif length(P_cand) > 1
        [inls_cand_num, inls_cand_idx] = max(inls_cand_num);
        inls_cand = inls_cand{inls_cand_idx};
        P_cand = P_cand{inls_cand_idx};
    else
        inls_cand_num = inls_cand_num(1);
        inls_cand = inls_cand{1};
        P_cand = P_cand{1};
    end

    %nonlinuear local optimization
    if inls_cand_num > 3
        lo_cnt = 0;
        while lo_cnt < 10
            lo_cnt = lo_cnt + 1;
            [lo_P, lo_inls, lo_inls_num] = ht_PnPnonlin(P_cand, u, X, inls_cand, rthr);
            if lo_inls_num >= inls_cand_num
                inls_cand_num = lo_inls_num;
                inls_cand = lo_inls;
                P_cand = lo_P;
            else
                break;
            end
        end
    end

    if inls_cand_num >= max_inlsnum
        max_inlsnum = inls_cand_num;
        P = P_cand;
        inls = inls_cand;
        max_iter = min([max_iter, nsamples(max_inlsnum, Npts, 3, 0.95)]);
    end

end

end

function [SampleCnt, q] = nsamples(ni, ptNum, pf, conf)
    q  = prod (((ni-pf+1) : ni) ./ ((ptNum-pf+1) : ptNum));
    if q < eps
      SampleCnt = Inf;
    else
      %   SampleCnt  = log(1 - conf) / log(1 - q);
      if q > conf
        SampleCnt = 1;
      else
        SampleCnt  = log(1 - conf) / log(1 - q);
      end
    end
end

function [inls, inls_num] = calculate_inls_angular(Pcand, u, X, rthr)
    inls = cell(1, length(Pcand));
    inls_num = zeros(1, length(Pcand));
    for ii = 1:1:length(Pcand)
        X_reproj = Pcand{ii} * [X; ones(1, size(X, 2))];
        X_reproj = bsxfun(@rdivide, X_reproj, sqrt(sum(X_reproj.^2, 1)));

        res = sum(u .* X_reproj, 1);
        inls{ii} = res > rthr;
        inls_num(ii) = sum(inls{ii});
    end

end

function [Poptim, inls, inls_num] = ht_PnPnonlin(Pcand, u, X, idx, rthr)
    Poptim = PnP_mex_wrapper( u(:, idx), X(:, idx), Pcand );

    X_reproj = Poptim * [X; ones(1, size(X, 2))];
    X_reproj = bsxfun(@rdivide, X_reproj, sqrt(sum(X_reproj.^2, 1)));
    res = sum(u .* X_reproj, 1);
    inls = res > rthr;
    inls_num = sum(inls);
end

author: HajimeTaira