text:
function hfig = tightfig(hfig)
% tightfig: Alters a figure so that it has the minimum size necessary to
% enclose all axes in the figure without excess space around them.
%
% Note that tightfig will expand the figure to completely encompass all
% axes if necessary. If any 3D axes are present which have been zoomed,
% tightfig will produce an error, as these cannot easily be dealt with.
%
% Input
%
% hfig - handle to figure, if not supplied, the current figure will be used
%   instead.
%
%

    if nargin == 0
        hfig = gcf;
    end

    % There can be an issue with tightfig when the user has been modifying
    % the contnts manually, the code below is an attempt to resolve this,
    % but it has not yet been satisfactorily fixed
%     origwindowstyle = get(hfig, 'WindowStyle');
    set(hfig, 'WindowStyle', 'normal');

    % 1 point is 0.3528 mm for future use

    % get all the axes handles note this will also fetch legends and
    % colorbars as well
    hax = findall(hfig, 'type', 'axes');
    % TODO: fix for modern matlab, colorbars and legends are no longer axes
    hcbar = findall(hfig, 'type', 'colorbar');
    hleg = findall(hfig, 'type', 'legend');

    % get the original axes units, so we can change and reset these again
    % later
    origaxunits = get(hax, 'Units');

    % change the axes units to cm
    set(hax, 'Units', 'centimeters');

    pos = [];
    ti = [];

    % get various position parameters of the axes
    if numel(hax) > 1
%         fsize = cell2mat(get(hax, 'FontSize'));
        ti = cell2mat(get(hax,'TightInset'));
        pos = [pos; cell2mat(get(hax, 'Position')) ];
    else
%         fsize = get(hax, 'FontSize');
        ti = get(hax,'TightInset');
        pos = [pos; get(hax, 'Position') ];
    end

    if ~isempty (hcbar)

        set(hcbar, 'Units', 'centimeters');

        % colorbars do not have tightinset property
        for cbind = 1:numel(hcbar)
            %         fsize = cell2mat(get(hax, 'FontSize'));
            [cbarpos, cbarti] = colorbarpos (hcbar);

            pos = [pos; cbarpos];
            ti = [ti; cbarti];
        end
    end

    if ~isempty (hleg)

        set(hleg, 'Units', 'centimeters');

        % legends do not have tightinset property
        if numel(hleg) > 1
            %         fsize = cell2mat(get(hax, 'FontSize'));
            pos = [pos; cell2mat(get(hleg, 'Position')) ];
        else
            %         fsize = get(hax, 'FontSize');
            pos = [pos; get(hleg, 'Position') ];
        end
        ti = [ti; repmat([0,0,0,0], numel(hleg), 1); ];
    end

    % ensure very tiny border so outer box always appears
    ti(ti < 0.1) = 0.15;

    % we will check if any 3d axes are zoomed, to do this we will check if
    % they are not being viewed in any of the 2d directions
    views2d = [0,90; 0,0; 90,0];

    for i = 1:numel(hax)

        set(hax(i), 'LooseInset', ti(i,:));
%         set(hax(i), 'LooseInset', [0,0,0,0]);

        % get the current viewing angle of the axes
        [az,el] = view(hax(i));

        % determine if the axes are zoomed
        iszoomed = strcmp(get(hax(i), 'CameraViewAngleMode'), 'manual');

        % test if we are viewing in 2d mode or a 3d view
        is2d = all(bsxfun(@eq, [az,el], views2d), 2);

        if iszoomed && ~any(is2d)
           error('TIGHTFIG:haszoomed3d', 'Cannot make figures containing zoomed 3D axes tight.')
        end

    end

    % we will move all the axes down and to the left by the amount
    % necessary to just show the bottom and leftmost axes and labels etc.
    moveleft = min(pos(:,1) - ti(:,1));

    movedown = min(pos(:,2) - ti(:,2));

    % we will also alter the height and width of the figure to just
    % encompass the topmost and rightmost axes and lables
    figwidth = max(pos(:,1) + pos(:,3) + ti(:,3) - moveleft);

    figheight = max(pos(:,2) + pos(:,4) + ti(:,4) - movedown);

    % move all the axes
    for i = 1:numel(hax)

        set(hax(i), 'Position', [pos(i,1:2) - [moveleft,movedown], pos(i,3:4)]);

    end

    for i = 1:numel(hcbar)

        set(hcbar(i), 'Position', [pos(i+numel(hax),1:2) - [moveleft,movedown], pos(i+numel(hax),3:4)]);

    end

    for i = 1:numel(hleg)

        set(hleg(i), 'Position', [pos(i+numel(hax)+numel(hcbar),1:2) - [moveleft,movedown], pos(i+numel(hax)+numel(hcbar),3:4)]);

    end

    origfigunits = get(hfig, 'Units');

    set(hfig, 'Units', 'centimeters');

    % change the size of the figure
    figpos = get(hfig, 'Position');

    set(hfig, 'Position', [figpos(1), figpos(2), figwidth, figheight]);

    % change the size of the paper
    set(hfig, 'PaperUnits','centimeters');
    set(hfig, 'PaperSize', [figwidth, figheight]);
    set(hfig, 'PaperPositionMode', 'manual');
    set(hfig, 'PaperPosition',[0 0 figwidth figheight]);

    % reset to original units for axes and figure
    if ~iscell(origaxunits)
        origaxunits = {origaxunits};
    end

    for i = 1:numel(hax)
        set(hax(i), 'Units', origaxunits{i});
    end

    set(hfig, 'Units', origfigunits);

%      set(hfig, 'WindowStyle', origwindowstyle);

end

function [pos, ti] = colorbarpos (hcbar)

    % 1 point is 0.3528 mm

    pos = hcbar.Position;
    ti = [0,0,0,0];

    if ~isempty (strfind (hcbar.Location, 'outside'))

        if strcmp (hcbar.AxisLocation, 'out')

            tlabels = hcbar.TickLabels;

            fsize = hcbar.FontSize;

            switch hcbar.Location

                case 'northoutside'

                    % make exta space a little more than the font size/height
                    ticklablespace_cm = 1.1 * (0.3528/10) * fsize;

                    ti(4) = ti(4) + ticklablespace_cm;

                case 'eastoutside'

                    maxlabellen = max ( cellfun (@numel, tlabels, 'UniformOutput', true) );

                    % 0.62 factor is arbitrary and added because we don't
                    % know the width of every character in the label, the
                    % fsize refers to the height of the font
                    ticklablespace_cm = (0.3528/10) * fsize * maxlabellen * 0.62;

                    ti(3) = ti(3) + ticklablespace_cm;

                case 'southoutside'

                    % make exta space a little more than the font size/height
                    ticklablespace_cm = 1.1 * (0.3528/10) * fsize;

                    ti(2) = ti(2) + ticklablespace_cm;

                case 'westoutside'

                    maxlabellen = max ( cellfun (@numel, tlabels, 'UniformOutput', true) );

                    % 0.62 factor is arbitrary and added because we don't
                    % know the width of every character in the label, the
                    % fsize refers to the height of the font
                    ticklablespace_cm = (0.3528/10) * fsize * maxlabellen * 0.62;

                    ti(1) = ti(1) + ticklablespace_cm;

            end

        end

    end

end

author: emtiyaz

repo: vadam

save_path: github-repos/MATLAB/emtiyaz-vadam


text:
% Test some operations involving chebfuns with delta functions.

function pass = test_deltaOps(pref)

if ( nargin < 1 )
    pref = chebfunpref();
end

tol = pref.deltaPrefs.deltaTol;

% These tests used to live in the scripts in tests/deltafun, but since they
% create chebfuns, they are better placed here.

% Some tests for cumsum() based on examples provided by LNT.
n = 6;
x = chebfun('x',[0 n]);
f = 0.5*sin(x);
for j = 1:n-1
  f = f + randn*dirac(x-j);
end
err = norm(diff(cumsum(f)) - f);
pass(1) = err < tol;

x = chebfun('x');
f = dirac(x-.5) + dirac(x) + dirac(x+.5) + heaviside(x);
err = norm(diff(cumsum(f)) - f);
pass(2) = err < tol;

x = chebfun('x');
f = sign(x)+sign(x-.5);
f2 = f(-1) + cumsum(diff(f));
err = norm(f - f2);
pass(3) = err < tol;

% A test for diff() based on an example from LNT.
savedPrefs = chebfunpref();
chebfunpref.setDefaults('enableDeltaFunctions', true);
try
    x = chebfun('x', [0 5]);
    f = 0.5*sin(x);
    A = randn(4, 1);
    for j = 1:4
      f = f + A(j)*dirac(x-j);
    end
    F = cumsum(.5*sin(x));
    for j = 1:4
      F = F + A(j)*heaviside(x-j);
    end
    err = norm(diff(F) - (f - f(0)));
    pass(4) = err < tol;
catch ME
    chebfunpref.setDefaults(savedPrefs);
    rethrow(ME)
end

chebfunpref.setDefaults(savedPrefs);

% A test involving innerProduct:
x = chebfun('x');
pass(5) = innerProduct(dirac(x), x) < eps;

% Test quasi matrix construction
x = chebfun('x');
A = [];
A = [ A, dirac(x-1) ];
A = [ A, dirac(x-0.5) ];
A = [ A, dirac(x-0.0) ];
pass(6) = (norm(A(:, 2) - dirac(x-.5)) < eps) && (size(A,2) == 3);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef stressTensor < tensor

  methods
    function sT = stressTensor(varargin)

      sT = sT@tensor(varargin{:},'rank',2);
      if ~sT.isSymmetric, warning('Tensor is not symmetric!'); end

    end
  end

  methods (Static = true)

    function sigma = load(varargin)
      T = load@tensor(varargin{:});
      sigma = stressTensor(T);
    end

    function sigma = uniaxial(v)
      % define uniaxial stress tensor
      %
      % Syntax
      %   sigma = stressTensor.uniaxial(v)
      %
      % Input
      %  v - @vector3d loading direction
      %
      % Output
      %  sigma - @stressTensor
      %

      sigma = stressTensor(dyad(v,2));

    end

    function sigma = rand(varargin)
      t = tensor.rand(varargin{:},'rank',2);
      sigma = stressTensor(t.sym);
    end

  end
end

author: mtex-toolbox

repo: mtex

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


text:
function [leftmap,rightmap] = mne_read_morph_map(from,to,subjects_dir)
%
% [leftmap,rightmap] = mne_read_morph_map(from,to,subjects_dir)
%
% Read the morphing map from subject 'from' to subject 'to'.
% If subjects_dir is not specified, the SUBJECTS_DIR environment
% variable is used
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%
%   Revision 1.1  2008/03/25 21:06:08  msh
%   Added mne_read_morph_map function
%
%
me='MNE:mne_read_morph_map';

global FIFF;
if isempty(FIFF)
    FIFF = fiff_define_constants();
end

if nargin < 3
    subjects_dir=getenv('SUBJECTS_DIR');
    if isempty(subjects_dir)
        error(me,'SUBJECTS_DIR not set');
    end
end
if nargin < 2
    error(me,'Not enough input arguments');
end
%
%  Does the file exist
%
name = sprintf('%s/morph-maps/%s-%s-morph.fif',subjects_dir,from,to);
if ~exist(name,'file')
    name = sprintf('%s/morph-maps/%s-%s-morph.fif',subjects_dir,to,from);
    if ~exist(name,'file')
        error(me,'The requested morph map does not exist');
    end
end
%
%    Open it
%
[fid,tree] = fiff_open(name);
%
%   Locate all maps
%
maps = fiff_dir_tree_find(tree,FIFF.FIFFB_MNE_MORPH_MAP);
if isempty(maps)
    fclose(fid);
    error(me,'Morphing map data not found');
end
%
%   Find the correct ones
%
for k = 1:length(maps)
    tag = find_tag(maps(k),FIFF.FIFF_MNE_MORPH_MAP_FROM);
    if strcmp(tag.data,from)
        tag = find_tag(maps(k),FIFF.FIFF_MNE_MORPH_MAP_TO);
        if strcmp(tag.data,to)
            %
            %  Names match: which hemishere is this?
            %
            tag = find_tag(maps(k),FIFF.FIFF_MNE_HEMI);
            if tag.data == FIFF.FIFFV_MNE_SURF_LEFT_HEMI
                tag     = find_tag(maps(k),FIFF.FIFF_MNE_MORPH_MAP);
                leftmap = tag.data;
                fprintf(1,'\tLeft-hemisphere map read.\n');
            elseif tag.data == FIFF.FIFFV_MNE_SURF_RIGHT_HEMI
                tag      = find_tag(maps(k),FIFF.FIFF_MNE_MORPH_MAP);
                rightmap = tag.data;
                fprintf(1,'\tRight-hemisphere map read.\n');
            end
        end
    end
end
fclose(fid);
if ~exist('leftmap')
    error(me,'Left hemisphere map not found in %s',name);
end
if ~exist('rightmap')
    error(me,'Left hemisphere map not found in %s',name);
end

return;

    function [tag] = find_tag(node,findkind)

        for p = 1:node.nent
            if node.dir(p).kind == findkind
                tag = fiff_read_tag(fid,node.dir(p).pos);
                return;
            end
        end
        tag = [];
        return;
    end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function F = times( F , G )
%.*   Times of two SPHEREFUNV objects.
%   F.*G if F is a SPHEREFUNV and G is double returns the SPHEREFUNV after
%   componentwise multiplication.
%
%   F.*G if F is a SPHEREFUN and G is a SPHEREFUNV returns the SPHEREFUNV
%   after multiplication of F by each component of G.
%
%   F.*G if F is a double and G is a SPHEREFUNV returns the SPHEREFUNV
%   after multiplication of F by each component of G.

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

% Empty check:
if ( isempty(F) || isempty(G) )
    F = spherefunv;
    return
end

if ( ~isa(F, 'spherefunv') )
    F = times(G, F);
    return
end

nF = 3;

if ( isa(G, 'double') )             % SPHEREFUNV.*double
    if ( numel(G) == 1 )            % SPHEREFUNV.*scalar
        scalar = G;
        for jj = 1 : nF
            F.components{jj} = times(F.components{jj}, scalar);
        end
    elseif ( ( size(G, 1) == nF ) || ( F.isTransposed && (size(G, 2) == nF) ) )
        for jj = 1 : nF
            F.components{jj} = times(F.components{jj}, G(jj));
        end
    else
        error('SPHEREFUN:SPHEREFUNV:times:double', ...
            'SPHEREFUNV and double size mismatch.');
    end

elseif ( isa(G, 'spherefunv') )      % SPHEREFUNV . * SPHEREFUNV
    for jj = 1:3
        F.components{jj} = times(F.components{jj}, G.components{jj});
    end

elseif ( isa(G, 'spherefun') )       % SPHEREFUN * SPHEREFUNV
    for jj = 1 : nF
        F.components{jj} = times(F.components{jj}, G);
    end

else  % error
    error( 'SPHEREFUN:SPHEREFUNV:times:inputs', ...
        'Unrecognized input arguments.');
end

end


text:
function L = nodesInLen(nodeNo,inLen,doNoExt,wt)
%NODESINLEN Length of the node input signal
%   Usage:  L = nodesInLen(nodeNo,inLen,doExt,treeStruct);
%
%   Input parameters:
%         nodeNo     : Node index.
%         inLen      : Filter thee input signal length.
%         doNoExt    : Expansive representation indicator.
%         wt         : Structure containing description of the filter tree.
%
%   Output parameters:
%         Lin        : Length of the node input signal
%
%   `nodesInLen(nodeNo,inLen,doExt,treeStruct)` return length of the input
%   signal of the node `nodeNo`. For definition of the structure see `wfbinit`.
%
%   See also: wfbtinit
%

L = zeros(numel(nodeNo),1);
for nn=1:length(nodeNo)
    subPat = [];
    filtLenPat = [];
    tmpNodeNo = nodeNo(nn);

    while(wt.parents(tmpNodeNo))
       parentNo = wt.parents(tmpNodeNo);
       tmpIdx = find(wt.children{parentNo}==tmpNodeNo);
       subPat(end+1) = wt.nodes{parentNo}.a(tmpIdx);
       filtLenPat(end+1) = length(wt.nodes{parentNo}.g{tmpIdx}.h);
       tmpNodeNo=parentNo;
    end

    subPat = subPat(end:-1:1);
    filtLenPat = filtLenPat(end:-1:1);

    L(nn) = inLen;
    if(~doNoExt)
        for ii=1:length(subPat)
            L(nn) = floor((L(nn)+filtLenPat(ii)-1)/subPat(ii));
        end
    else
        for ii=1:length(subPat)
            L(nn) = ceil(L(nn)/subPat(ii));
        end
    end
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function s = fastSampen(y,m,r)
% INPUTS
%  y     num_var x num_samples
%  m     template length
%  r     radius
% OUTPUTS
%  s     Sample Entropy
%
% Written by Shamim Nemati <shamim.nemati@alum.mit.edu>
%
% 01-31-2018 : modified by Giulia Da Poian <giulia.dap@gmail.com>, see
%              comments in the code
%              NOTE : Sample entropy quantifies the likelihood that a
%              sequence of m consecutive data points that matches another
%              sequence of the same length (match within a tolerance of r)
%              will still match the other sequence when their length is
%              increased of one sample (sequences of length m + 1);
%              References :
%              1. Richman JS, Moorman JR. Physiological time-series
%              analysis using approximate entropy and sample entropy.
%              American Journal of Physiology-Heart and Circulatory
%              Physiology. 2000 Jun 1;278(6):H2039-49.
%              2. Humeau-Heurtier A. The multiscale entropy
%              algorithm and its variants: A review. Entropy. 2015 May
%              12;17(5):3110-23.
% LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information

% Ensure y is a row vector
if size(y, 1) > size(y,2)
    y = y';
end

xx = convert_to_lagged_form(y, m)';   % Giulia: replaced m-1 with m to match SampEn definition
Dxx = pdist(xx,'chebychev');

yy = convert_to_lagged_form(y, m+1)'; % Giulia: replaced m with m+1 to match SampEn definition
Dyy = pdist(yy,'chebychev');

A = mean( Dxx < r ) ;
B = mean( Dyy < r );

s = -log(B/A);

function yy = convert_to_lagged_form(y, k)
% Create an observation vector yy(:,t) containing the last k values of y, newest first
% e.g., k=2, y = (a1 a2 a3)     yy  = a2 a3
%                (b1 b2 b3)           b2 b2
%                                     a1 a2
%                                     b1 b2
[s, T] = size(y);
bs = s*ones(1,k);
yy = zeros(k*s, T-k+1);
for i=1:k,     yy(block(i,bs), :) = y(:, k-i+1:end-i+1); end

function sub = block(blocks, block_sizes)
% BLOCK Return a vector of subscripts corresponding to the specified blocks.
% sub = block(blocks, block_sizes)
%
% e.g., block([2 5], [2 1 2 1 2]) = [3 7 8].
blocks = blocks(:)';
block_sizes = block_sizes(:)';
skip = [0 cumsum(block_sizes)];
start = skip(blocks)+1;
fin = start + block_sizes(blocks) - 1;
sub = [];
for j=1:length(blocks)
    sub = [sub start(j):fin(j)];
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function im = rgb_to_nrgb(im)
% Transforms image into NRGB color space
% Code by Ross Girshick

if ~isa(im, 'double')
  im = double(im);
end
im = im./repmat(sum(im,3), [1 1 3]);

im = 255*im;

author: Cloud-CV

repo: object-proposals


text:
function [strshocks_estimates]=strsestimates(strshocks_record,n,T,IRFband)

% create first the cell that will contain the estimates
strshocks_estimates=cell(n,1);

% for each variable and each sample period, compute the median, lower and upper bound from the Gibbs sampler records
% consider variables in turn
for ii=1:n
   % consider sample periods in turn
   for jj=1:T
   % compute first the lower bound
   strshocks_estimates{ii,1}(1,jj)=quantile(strshocks_record{ii,1}(:,jj),(1-IRFband)/2);
   % then compute the median
   strshocks_estimates{ii,1}(2,jj)=quantile(strshocks_record{ii,1}(:,jj),0.5);
   % finally compute the upper bound
   strshocks_estimates{ii,1}(3,jj)=quantile(strshocks_record{ii,1}(:,jj),(1-(1-IRFband)/2));
   end
end

author: european-central-bank

repo: BEAR-toolbox


text:
% angle to wall
function [data,units] = compute_angle2wall_rect(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);
  x = trx(fly).x_mm;
  y = trx(fly).y_mm;

  walla = zeros(1,4);
  walla(1) = atan2(trx.landmark_params{n}.tr_y(fly)-trx.landmark_params{n}.tl_y(fly),trx.landmark_params{n}.tr_x(fly)-trx.landmark_params{n}.tl_x(fly));
  walla(2) = atan2(trx.landmark_params{n}.br_y(fly)-trx.landmark_params{n}.tr_y(fly),trx.landmark_params{n}.br_x(fly)-trx.landmark_params{n}.tr_x(fly));
  walla(3) = atan2(trx.landmark_params{n}.bl_y(fly)-trx.landmark_params{n}.br_y(fly),trx.landmark_params{n}.bl_x(fly)-trx.landmark_params{n}.br_x(fly));
  walla(4) = atan2(trx.landmark_params{n}.tl_y(fly)-trx.landmark_params{n}.bl_y(fly),trx.landmark_params{n}.tl_x(fly)-trx.landmark_params{n}.bl_x(fly));

  [dtop] = getDist(trx.landmark_params{n}.tl_x(fly),trx.landmark_params{n}.tl_y(fly), trx.landmark_params{n}.tr_x(fly),trx.landmark_params{n}.tr_y(fly),...
    x,y);
  [dright]= getDist(trx.landmark_params{n}.tr_x(fly),trx.landmark_params{n}.tr_y(fly), trx.landmark_params{n}.br_x(fly),trx.landmark_params{n}.br_y(fly),...
    x,y);
  [dbottom] = getDist(trx.landmark_params{n}.bl_x(fly),trx.landmark_params{n}.bl_y(fly), trx.landmark_params{n}.br_x(fly),trx.landmark_params{n}.br_y(fly),...
    x,y);
  [dleft] = getDist(trx.landmark_params{n}.tl_x(fly),trx.landmark_params{n}.tl_y(fly), trx.landmark_params{n}.bl_x(fly),trx.landmark_params{n}.bl_y(fly),...
    x,y);

  [~,closest] = min([dtop;dright; dbottom; dleft],[],1);

  closest_walla = walla(closest);
  targeta = trx(fly).theta;

  a = mod(targeta-closest_walla+pi/2+pi,2*pi)-pi;

  data{i} = a;
end
units = parseunits('rad');

function [d,a] = getDist(p1_x,p1_y,p2_x,p2_y,p_x,p_y)

dp1p2 = (p1_x-p2_x).^2 + (p1_y-p2_y).^2;

dotpr_x = (p1_x - p_x).*(p1_x-p2_x);
dotpr_y = (p1_y - p_y).*(p1_y-p2_y);

t = (dotpr_x+dotpr_y)/dp1p2;
proj_x = p1_x + t.*(p2_x-p1_x);
proj_y = p1_y + t.*(p2_y-p1_y);
p2proj_x = p_x-proj_x;
p2proj_y = p_y-proj_y;

d = sqrt(  (p2proj_x).^2 + (p2proj_y).^2);

author: kristinbranson

repo: JAABA


text:
function [ni] = niftiApplyAndCreateXform(nii,xformType)
%First create then apply a specified transform onto the supplied nifti
%struct.
%
% USAGE
%  nii = readNifti(niftiFullPath);
%  xformType = 'Inplane';
%  niftiApplyAndCreateXform(nii,xformType);
%
% INPUTS
%  Nifti struct
%  String specifying the transform to apply
%
% RETURNS
%  Nifti Struct
%
%
% Copyright Stanford VistaLab 2013

xformType = mrvParamFormat(xformType);

xform = niftiCreateXform(nii,xformType);
%Create transform
if (xform == 0)
   %We shouldn't do anything
   ni = nii;
else
   ni = niftiApplyXform(nii,xform);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function vo = getOrthogonalVector(v)

vn = v ./ norm(v);

vo = [0,0,0];
while (isequal(vo,[0,0,0]))
    vo = rand(1,3)-0.5;
    vc = vo*vn';
    vo = vo-(vn*vc);
end

vo = vo ./ norm(vo);


text:
function [total_error, max_error] = eq_area_error(dim,N)
%EQ_AREA_ERROR Total area error and max area error per region of an EQ partition
%
%Syntax
% [total_error, max_error] = eq_area_error(dim,N)
%
%Description
% [TOTAL_ERROR, MAX_ERROR] = EQ_AREA_ERROR(dim,N) does the following:
% 1) uses the recursive zonal equal area sphere partitioning algorithm to
%    partition the unit sphere S^dim into N regions,
% 2) sets TOTAL_ERROR to be the absolute difference between the total area of
%    all regions of the partition, and the area of S^dim, and
% 3) sets MAX_ERROR to be the maximum absolute difference between the area of
%    any region of the partition, and the ideal area of a region as given by
%    AREA_OF_IDEAL_REGION(dim,N), which is 1/N times the area of S^dim.
%
% The argument dim must be a positive integer.
% The argument N must be a positive integer or an array of positive integers.
% The results TOTAL_ERROR and MAX_ERROR will be arrays of the same size as N.
%
%Examples
% > [total_error, max_error] = eq_area_error(2,10)
% total_error =
%    1.7764e-15
%
% max_error =
%    4.4409e-16
%
% > [total_error, max_error] = eq_area_error(3,1:6)
% total_error =
%    1.0e-12 *
%     0.0036    0.0036    0.1847    0.0142    0.0142    0.2132
%
% max_error =
%    1.0e-12 *
%     0.0036    0.0018    0.1954    0.0284    0.0440    0.0777
%
%See also
% EQ_REGIONS, AREA_OF_SPHERE, AREA_OF_IDEAL_REGION

% Copyright 2004-2005 Paul Leopardi for the University of New South Wales.
% $Revision 1.10 $ $Date 2005-06-01 $
% Documentation files renamed
% $Revision 1.00 $ $Date 2005-02-12 $
%
% For licensing, see COPYING.
% For references, see AUTHORS.
% For revision history, see CHANGELOG.

%
% Check number of arguments
%
error(nargchk(2,2,nargin));
error(nargoutchk(2,2,nargout));

%
% Flatten N into a row vector.
%
shape = size(N);
n_partitions = prod(shape);
N = reshape(N,1,n_partitions);

total_error = zeros(size(N));
max_error   = zeros(size(N));
sphere_area = area_of_sphere(dim);

for partition_n = 1:n_partitions
    n = N(partition_n);
    regions = eq_regions(dim,n);
    ideal_area = area_of_ideal_region(dim,n);
    total_area = 0;
    for region_n = 1:size(regions,3)
        area = area_of_region(regions(:,:,region_n));
        total_area = total_area + area;
        region_error = abs(area - ideal_area);
        if region_error > max_error(partition_n)
            max_error(partition_n) = region_error;
        end
    end
    total_error(partition_n) = abs(sphere_area - total_area);
end
%
% Reshape output to same array size as original N.
%
total_error = reshape(total_error,shape);
max_error = reshape(max_error,shape);
%
% end function

function area = area_of_region(region)
%AREA_OF_REGION Area of given region
%
% area = area_of_region(region);

dim = size(region,1);
s_top = region(dim,1);
s_bot = region(dim,2);
if dim > 1
    area = area_of_collar(dim, s_top, s_bot)*area_of_region(region(1:dim-1,:))/area_of_sphere(dim-1);
else
    if s_bot == 0
        s_bot = 2*pi;
    end
    if s_top == s_bot
        s_bot = s_top + 2*pi;
    end
    area = s_bot - s_top;
end
%
% end function

repo: gspbox


text:
function wedge_integrals_test ( )

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

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

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function plotRelativeProjectedAmps(view)
%
% plotRelativeProjectedAmplitudes(view)
%
% Bar plot of the amplitudes for each scan, averaging across
% all pixels (in all slices) in the current ROI.  The bar heights
% and a coarse SEM can be obtained from get(gca,'UserData').
%
% gmb  5/25/98
% bw   2/19/99  Added seY field to the UserData field.
%	    seY is an estimate of the variability in the
%      amplitudes.  It is the SEM of the in the complex
%      (amp*exp(-i*ph)) representation.  The values are
%      computed in vectorMean.m
% fwc   11/07/02 plots data relative to current view

% Compute means across scans, for all pixels in the
% currently selected ROI.  The seZ value is the mean
% distance from the mean.
[meanAmps,meanPhs,seZ] = vectorMeans(view);

%Reference scan is the current scan
refScan = getCurScan(view);

% Compute the amplitude projected onto the reference phase
meanProjectedAmps = meanAmps.*cos(meanPhs-meanPhs(refScan));

meanRelProjAmps=meanProjectedAmps/meanProjectedAmps(refScan);

selectGraphWin

% Header
ROIname = view.ROIs(view.selectedROI).name;
headerStr = ['Mean of projected amplitudes relative to reference scan, ROI ',ROIname];
set(gcf,'Name',headerStr);

%plot the bar graph
clf
fontSize = 14;
h=mybar(meanRelProjAmps);
xlabel('Scan','FontSize',fontSize);
ylabel('Relative Mean Projected Amplitude','FontSize',fontSize);
ylim =get(gca,'YLim');
set(gca,'YLim',ylim*1.1);
set(gca,'FontSize',fontSize);

foo=cell2struct(h,'bar');
hbar=foo(refScan).bar;
set(hbar,'FaceColor','r')

%Save the data in gca('UserData')
data.y = meanRelProjAmps;
data.refScan = refScan;
data.seY = seZ/meanProjectedAmps(refScan); % this should probably be adapted

set(gca,'UserData',data);

return;


min_cost_0_4000 = min_cost_0;
min_cost_1_4000 = min_cost_1;
min_cost_2_4000 = min_cost_2;
min_cost_3_4000 = min_cost_3;
min_cost_4_4000 = min_cost_4;

totalNum_0_4000 = totalNum_0;
totalNum_1_4000 = totalNum_1;
totalNum_2_4000 = totalNum_2;
totalNum_3_4000 = totalNum_3;
totalNum_4_4000 = totalNum_4;


text:
function varargout = EarthGradient(f,varargin)
%GRADIENT Approximate gradient.
%   [FX,FY] = GRADIENT(F) returns the numerical gradient of the
%   matrix F. FX corresponds to dF/dx, the differences in x (horizontal)
%   direction. FY corresponds to dF/dy, the differences in y (vertical)
%   direction. The spacing between points in each direction is assumed to
%   be one. When F is a vector, DF = GRADIENT(F) is the 1-D gradient.
%
%   [FX,FY] = GRADIENT(F,H), where H is a scalar, uses H as the
%   spacing between points in each direction.
%
%   [FX,FY] = GRADIENT(F,HX,HY), when F is 2-D, uses the spacing
%   specified by HX and HY. HX and HY can either be scalars to specify
%   the spacing between coordinates or vectors to specify the
%   coordinates of the points.  If HX and HY are vectors, their length
%   must match the corresponding dimension of F.
%
%   [FX,FY,FZ] = GRADIENT(F), when F is a 3-D array, returns the
%   numerical gradient of F. FZ corresponds to dF/dz, the differences
%   in the z direction. GRADIENT(F,H), where H is a scalar,
%   uses H as the spacing between points in each direction.
%
%   [FX,FY,FZ] = GRADIENT(F,HX,HY,HZ) uses the spacing given by
%   HX, HY, HZ.
%
%   [FX,FY,FZ,...] = GRADIENT(F,...) extends similarly when F is N-D
%   and must be invoked with N outputs and either 2 or N+1 inputs.
%
%   Note: The first output FX is always the gradient along the 2nd
%   dimension of F, going across columns.  The second output FY is always
%   the gradient along the 1st dimension of F, going across rows.  For the
%   third output FZ and the outputs that follow, the Nth output is the
%   gradient along the Nth dimension of F.
%
%   Examples:
%       [x,y] = meshgrid(-2:.2:2, -2:.2:2);
%       z = x .* exp(-x.^2 - y.^2);
%       [px,py] = gradient(z,.2,.2);
%       contour(z), hold on, quiver(px,py), hold off
%
%   Class support for input F:
%      float: double, single
%
%   See also DIFF, DEL2.

%   Copyright 1984-2015 The MathWorks, Inc.

[f,ndim,loc,rflag] = parse_inputs(f,varargin);
nargoutchk(0,ndim);

% Loop over each dimension.

varargout = cell(1,ndim);
siz = size(f);
% first dimension
g  = zeros(size(f),class(f)); % case of singleton dimension
h = loc{1};
n = siz(1);
% Take forward differences on left and right edges
if n > 1
   g(1,:) = (f(2,:) - f(1,:))./h;
   g(n,:) = (f(n,:) - f(n-1,:))./h;
end

% Take centered differences on interior points
if n > 2
   h = 2*h;
   for nn = 2:n-1
       g(nn,:) = (f(nn+1,:)-f(nn-1,:)) ./ h;
   end
end

varargout{1} = g;

% second dimensions and beyond
if ndim == 2
    % special case 2-D matrices to support sparse matrices,
    % which lack support for N-D operations including reshape
    % and indexing
    n = siz(2);
    h = loc{2};
    g = zeros(size(f),class(f));

    % Take forward differences on left and right edges
    if n > 1
        g(:,1) = (f(:,2) - f(:,1))./h;
        g(:,n) = (f(:,n) - f(:,n-1))./h;
    end

    % Take centered differences on interior points
    if n > 2
        h = 2*h;
        g(:,2:n-1) = (f(:,3:n) - f(:,1:n-2)) ./ h;
    end
    varargout{2} = g;

elseif ndim > 2
    % N-D case
    for k = 2:ndim
        n = siz(k);
        newsiz = [prod(siz(1:k-1)) siz(k) prod(siz(k+1:end))];
        nf = reshape(f,newsiz);
        h = reshape(loc{k},1,[]);
        g  = zeros(size(nf),class(nf)); % case of singleton dimension

        % Take forward differences on left and right edges
        if n > 1
            g(:,1,:) = (nf(:,2,:) - nf(:,1,:))/(h(2)-h(1));
            g(:,n,:) = (nf(:,n,:) - nf(:,n-1,:))/(h(end)-h(end-1));
        end

        % Take centered differences on interior points
        if n > 2
            h = h(3:n) - h(1:n-2);
            g(:,2:n-1,:) = (nf(:,3:n,:) - nf(:,1:n-2,:)) ./ h;
        end

        varargout{k} = reshape(g,siz);
    end
end

% Swap 1 and 2 since x is the second dimension and y is the first.
if ndim > 1
    varargout(2:-1:1) = varargout(1:2);
elseif rflag
    varargout{1} = varargout{1}.';
end

%-------------------------------------------------------
function [f,ndim,loc,rflag] = parse_inputs(f,v)
%PARSE_INPUTS
%   [ERR,F,LOC,RFLAG] = PARSE_INPUTS(F,V) returns the spacing
%   LOC along the x,y,z,... directions and a row vector
%   flag RFLAG.

loc = {};

% Flag vector case and row vector case.
ndim = ndims(f);
vflag = false;
rflag = false;
if isvector(f)
    ndim = 1;
    vflag = true;
    if isrow(f) % Treat row vector as a column vector
        rflag = true;
        f = f.';
    end
end

indx = size(f);

% Default step sizes: hx = hy = hz = 1
if isempty(v)
    % gradient(f)
    loc = cell(1, ndims(f));
    for k = 1:ndims(f)
        loc(k) = {1:indx(k)};
    end
elseif isscalar(v) % gradient(f,h)
    % Expand scalar step size
    if isscalar(v{1})
        loc = cell(1, ndims(f));
        for k = 1:ndims(f)
            h = v{1};
            loc(k) = {h*(1:indx(k))};
        end
        % Check for vector case
    elseif vflag
        loc(1) = v(1);
    else
        error(message('MATLAB:gradient:InvalidInputs'));
    end
elseif ndims(f) == numel(v)  % gradient(f,hx,hy,hz,...)
    % Swap 1 and 2 since x is the second dimension and y is the first.
    loc = v;
    if ndim > 1
        loc(2:-1:1) = loc(1:2);
    end
else
    error(message('MATLAB:gradient:InvalidInputs'));
end

author: CHLNDDEV


text:
function mu = bellmf(z,a,b)
%BELLMF Bell-shaped membership function.
%   MU = BELLMF(Z,A,B) computes the bell-shaped fuzzy membership
%   function. Z is the input variable and can be a vector of any length.
%   A and B are scalar shape parameters, ordered such that A <= B.
%
%     MU = SMF(Z,A,B)         Z < B
%     MU = SMF(2*B - Z,A,B)   B <= Z
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

mu = zeros(size(z));

left_side = z < b;
mu(left_side) = smf(z(left_side),a,b);

right_side = z >= b;
mu(right_side) = smf(2*b - z(right_side),a,b);

author: dipum

repo: dipum-toolbox

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


text:
function [data] = ft_denoise_ssp(cfg, data)

% take ssp projectors from data.hdr and construct a montage from them

author: spm


text:
function [params,names] = lmcKernExtractParam(kern)

% LMCKERNEXTRACTPARAM Extract parameters from the LMC kernel struc.
% FORMAT
% DESC Extract parameters from the linear model of coregionalization kernel
% structure into a vector of parameters for optimisation.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the
% kernel. The vector of transforms is assumed to be empty here, any
% transormation of parameters is assumed to be done in the
% component kernels.
%
% FORMAT
% DESC the same that before but also returns the names of the parameters.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the
% kernel. The vector of transforms is assumed to be empty here, any
% transormation of parameters is assumed to be done in the
% component kernels.
% RETURN names : cell array of strings containing parameter names.
%
% SEEALSO multiKernExtractParam
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

% First extract the parameters of the basic kernel

fhandle = str2func([kern.basicKernelType 'KernExtractParam']);

if nargout > 1
  [paramsTemp, namesTemp] = fhandle(kern);
  namesB = cell(kern.nout, kern.rankCorregMatrix);
  for i = 1:kern.nout
      for j =1:kern.rankCorregMatrix
          namesB{i,j} = ['A(' num2str(i) ',' num2str(j) ')'];
      end
  end
  names = {namesTemp{1:kern.nParamsBK}, namesB{:}};
else
  paramsTemp = fhandle(kern);
end

% Add the parameters of the corregionalization matrix

params = [paramsTemp(1:kern.nParamsBK)  kern.A(:)'];

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function download_train(data_dir)
% DOWNLOAD_TRAIN Checks, and, if required, downloads the necessary data and networks for the training.
%
%   download_train(DATA_ROOT) checks if the data and networks necessary for running the training script exist.
%   If not it downloads it in the folder structure:
%     DATA_ROOT/train/dbs/         : folder with training database mat files
%     DATA_ROOT/train/ims/         : folder with original images used for training
%     DATA_ROOT/networks/imagenet/ : CNN models pretrained for classification using imagenet data

    % Create data folder if it does not exist
    if ~exist(data_dir, 'dir')
        mkdir(data_dir);
    end

    % Create train folder if it does not exist
    train_dir = fullfile(data_dir, 'train');
    if ~exist(train_dir, 'dir')
        mkdir(train_dir);
    end

    % Download folder train/db/
    src_dir = fullfile('http://cmp.felk.cvut.cz/cnnimageretrieval/data', 'train', 'dbs');
    dst_dir = fullfile(data_dir, 'train', 'dbs');
    dl_files = {'retrieval-SfM-30k.mat',  'retrieval-SfM-30k-whiten.mat', ...
                'retrieval-SfM-120k.mat', 'retrieval-SfM-120k-whiten.mat'};
    if ~exist(dst_dir, 'dir')
        fprintf('>> Database directory does not exist. Creating: %s\n', dst_dir);
        mkdir(dst_dir);
        fprintf('>> Downloading database files from cmp.felk.cvut.cz/cnnimageretrieval\n');
    end
    for i = 1:numel(dl_files)
        src_file = fullfile(src_dir, dl_files{i});
        dst_file = fullfile(dst_dir, dl_files{i});
        if ~exist(dst_file, 'file')
            fprintf('>> DB file %s does not exist. Downloading...\n', dl_files{i});
            system(sprintf('wget %s -O %s', src_file, dst_file));
        end
    end

    % Download folder train/ims/
    src_dir = fullfile('http://cmp.felk.cvut.cz/cnnimageretrieval/data', 'train', 'ims');
    dst_dir = fullfile(data_dir, 'train', 'ims');
    dl_file = 'ims.tar.gz';
    if ~exist(dst_dir, 'dir')
        src_file = fullfile(src_dir, dl_file);
        dst_file = fullfile(dst_dir, dl_file);
        fprintf('>> Image directory does not exist. Creating: %s\n', dst_dir);
        mkdir(dst_dir);
        fprintf('>> Downloading ims.tar.gz...\n');
        system(sprintf('wget %s -O %s', src_file, dst_file));
        fprintf('>> Extracting %s...\n', dst_file);
        system(sprintf('tar -zxf %s -C %s', dst_file, dst_dir));
        fprintf('>> Extracted, deleting %s...\n', dst_file);
        system(sprintf('rm %s', dst_file));
    end
    dl_files = {'retrieval-SfM-30k-imagenames-clusterids.mat', 'retrieval-SfM-120k-imagenames-clusterids.mat'};
    for i = 1:numel(dl_files)
        src_file = fullfile(src_dir, dl_files{i});
        dst_file = fullfile(dst_dir, dl_files{i});
        if ~exist(dst_file, 'file')
            fprintf('>> Dataset image split %s does not exist. Downloading...\n', dl_files{i});
            system(sprintf('wget %s -O %s', src_file, dst_file));
        end
    end

    % Download folder networks/imagenet/
    src_dir = fullfile('http://www.vlfeat.org/matconvnet/', 'models');
    dst_dir = fullfile(data_dir, 'networks', 'imagenet');
    dl_files = {'imagenet-caffe-alex.mat', 'imagenet-vgg-verydeep-16.mat', 'imagenet-googlenet-dag.mat', ...
                'imagenet-resnet-50-dag.mat', 'imagenet-resnet-101-dag.mat', 'imagenet-resnet-152-dag.mat'};
    if ~exist(dst_dir, 'dir')
        fprintf('>> Imagenet networks directory does not exist. Creating: %s\n', dst_dir);
        mkdir(dst_dir);
        fprintf('>> Downloading imagenet networks from http://www.vlfeat.org/matconvnet\n');
    end
    for i = 1:numel(dl_files)
        src_file = fullfile(src_dir, dl_files{i});
        dst_file = fullfile(dst_dir, dl_files{i});
        if ~exist(dst_file, 'file')
            fprintf('>> Network %s does not exist. Downloading...\n', dl_files{i});
            system(sprintf('wget %s -O %s', src_file, dst_file));
        end
    end

author: filipradenovic

repo: cnnimageretrieval


text:
function [Z, H, dnorm] = seminmf( X, k, varargin )
% Matrix sizes
% X: m x n
% Z: m x num_of_components
% H: num_of_components x num_of_components

% Process optional arguments
pnames = {'z0' 'h0' 'bUpdateH' 'maxiter' 'TolFun' 'bUpdateZ' 'verbose'};

% Do SVD initialisation of the init components

[z0, h0] = NNDSVD(abs(X), k, 0);

dflts  = {z0, h0, 1, 300,  1e-5, 1, 1};

[Z, H, bUpdateH, max_iter, tolfun, bUpdateZ, verbose] = ...
        internal.stats.parseArgs(pnames,dflts,varargin{:});

if exist(path, 'file') ~= 0
    load(path);
    return;
end

for i = 1:max_iter

    if bUpdateZ
        Z = X * pinv(H);
    end

    A = Z' * X;
    Ap = (abs(A)+A)./2;
    An = (abs(A)-A)./2;

    B = Z' * Z;
    Bp = (abs(B)+B)./2;
    Bn = (abs(B)-B)./2;

    if bUpdateH
        H = H .* sqrt((Ap + Bn * H) ./ (An + Bp * H + eps));
    end

    if mod(i, 10) == 0 || mod(i+1, 10) == 0

        s = X - Z * H;
        dnorm = sqrt(sum(s(:).^2));
        % dnorm = norm(gX - Z * H, 'fro');

        if mod(i+1, 10) == 0
            dnorm0 = dnorm;
            continue
        end

        if mod(i, 100) == 0 && verbose
            display(sprintf('...Semi-NMF iteration #%d out of %d, error: %f\n', i, max_iter, dnorm));
        end

        if 1 && exist('dnorm0')
            assert(dnorm <= dnorm0, sprintf('Rec. error increasing! From %f to %f. (%d)', dnorm0, dnorm, k));
        end

        % Check for convergence
        if exist('dnorm0') && dnorm0-dnorm <= tolfun*max(1,dnorm0)
            if verbose
                display(sprintf('Stopped at %d: dnorm: %f, dnorm0: %f', i, dnorm, dnorm0));
            end
            break;
        end

    end
end

author: jwyang


text:
function pdf = reciprocal_pdf ( x, a, b )

%*****************************************************************************80
%
%% RECIPROCAL_PDF evaluates the Reciprocal PDF.
%
%  Formula:
%
%    PDF(X)(A,B) = 1.0D+00 / ( X * LOG ( B / A ) )
%    for 0.0D+00 <= X
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < A <= B.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x <= 0.0 )
    pdf = 0.0;
  elseif ( 0.0 < x )
    pdf = 1.0 / ( x * log ( b / a ) );
  end

  return
end


text:
% setup.m
% run this file to set up matlab path etc.
% you may need to modify this depending on how you installed the toolbox
% so this should be considered simply a "guide" not a robust script.

if ~exist('irtdir', 'var')
	disp('The variable "irtdir" is not set, so trying default, assuming')
	disp('that you launched matlab from the irt install directory.')
	disp('You may need to edit setup.m or adjust your path otherwise.')

%	irtdir = pwd; % the default is to assume launch from irt directory

	% default is to look for directory where this setup.m is installed!
	irtdir = which('setup'); % find setup.m
	[irtdir, ~] = fileparts(irtdir);

	disp(['Assuming you installed irt in directory "' irtdir '".'])

%	irtdir = '~fessler/l/src/matlab/alg/'; % where you install this package
%	irtdir = '~fessler/l/web/irt/'; % where you install this package
end

if ~exist(irtdir, 'dir')
	disp(sprintf('The directory "%s" does not exist', irtdir))
	error(sprintf('you need to edit %s to change default path', mfilename))
end

if irtdir(end) ~= filesep % make sure there is a '/' at end of directory
	irtdir = [irtdir filesep];
end

list = {...
'align', ...		% image registration
'align/mex', ...	% image registration mex files
'blob', ...		% blob (KB) basis
'ct', ...		% x-ray CT (polyenergetic) recon
'data', ...		% example data
'data/downloads', ...		% example data
'emission', ...		% emission image reconstruction
'example', ...		% example applications
'fbp', ...		% FBP (filtered backprojection) code
'general', ...		% generic image reconstruction
'graph', ...		% graphics routines
'mri', ...		% MRI reconstruction
'mri/fieldmap', ...	% MRI B0 field map estimation
'mri-rf/yip-spsp', ...	% MRI RF pulse design
'mri/arch',... % MRI examples
'nufft', ...		% nonuniform FFT (for a fast projector)
'nufft/table', ...	% mex files for NUFFT
'penalty', ...		% regularizing penalty functions
'systems', ...		% system "matrices"
'systems/tests', ...	% tests of systems
'transmission', ...	% transmission image reconstruction
'utilities', ...	% various utility functions
'wls', ...		% weighted least-squares (WLS) estimates
'um', ...		% extra directories for UM users only
'um/lustig-poisson-disk'
};
%'mri/recon', ...	% MRI reconstruction - old

for ii=1:numel(list)
	tmp = [irtdir list{ii}];
	if exist(tmp, 'dir'), addpath(tmp), end
end

% Set up path to mex files, possibly depending on matlab version.
% Fortunately it seems that the v6 files also run on v7.
% Unfortunately, it seems that sometimes compiling on one version
% e.g., 7.3 and running on earlier 7.0 won't work.
% If you have mex problems, comment out the following line.
% Much of the toolbox will work without mex, just slower.
addpath([irtdir 'mex/v7']);

% do not add the paths below if you are using Matlab!
if ir_is_octave
	addpath([irtdir 'octave']); % extra stuff for octave only!
elseif isempty(which('dbstack')) % for freemat only!
	addpath([irtdir 'freemat']); % extra stuff for freemat only!
end

% check to see if path setup worked by looking for im() routine.
if strcmp([irtdir 'graph' filesep 'im.m'], which('im'))
	disp('Path setup for irt appears to have succeeded.')
	clear list ii irtdir tmp
else
	disp('Path setup for irt may have failed.')
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [CIJkcore,kn,peelorder,peellevel] = kcore_bu(CIJ,k)
%KCORE_BU       K-core
%
%   [CIJkcore,kn,peelorder,peellevel] = kcore_bu(CIJ,k);
%
%   The k-core is the largest subnetwork comprising nodes of degree at
%   least k. This function computes the k-core for a given binary
%   undirected connection matrix by recursively peeling off nodes with
%   degree lower than k, until no such nodes remain.
%
%   input:          CIJ,        connection/adjacency matrix (binary, undirected)
%                     k,        level of k-core
%
%   output:    CIJkcore,        connection matrix of the k-core.  This matrix
%                               only contains nodes of degree at least k.
%                    kn,        size of k-core
%                    peelorder, indices in the order in which they were
%                               peeled away during k-core decomposition
%                    peellevel, corresponding level - nodes at the same
%                               level were peeled away at the same time
%
%   'peelorder' and 'peellevel' are similar the the k-core sub-shells
%   described in Modha and Singh (2010).
%
%   References: e.g. Hagmann et al. (2008) PLoS Biology
%
%   Olaf Sporns, Indiana University, 2007/2008/2010/2012

%#ok<*AGROW>

peelorder = [];
peellevel = [];
iter = 0;

while 1

    % get degrees of matrix
    [deg] = degrees_und(CIJ);

    % find nodes with degree <k
    ff = find((deg<k)&(deg>0));

    % if none found -> stop
    if (isempty(ff)) break; end;            %#ok<SEPEX>

    % peel away found nodes
    iter = iter+1;
    CIJ(ff,:) = 0;
    CIJ(:,ff) = 0;

    peelorder = [peelorder; ff'];
    peellevel = [peellevel; iter.*ones(1,length(ff))'];

end;

CIJkcore = CIJ;
kn = sum(deg>0);


text:
function [V,D]=gabmuleigs(K,c,p3,varargin)
%GABMULEIGS  Eigenpairs of Gabor multiplier
%   Usage:  h=gabmuleigs(K,c,g,a);
%           h=gabmuleigs(K,c,a);
%           h=gabmuleigs(K,c,ga,gs,a);
%
%   Input parameters:
%         K     : Number of eigenvectors to compute.
%         c     : symbol of Gabor multiplier
%         g     : analysis/synthesis window
%         ga    : analysis window
%         gs    : synthesis window
%         a     : Length of time shift.
%   Output parameters:
%         V     : Matrix containing eigenvectors.
%         D     : Eigenvalues.
%
%   `gabmuleigs` has been deprecated. Please use construct a frame multiplier
%   and use |framemuleigs| instead.
%
%   A call to `gabmuleigs(K,c,ga,gs,a)` can be replaced by ::
%
%     [Fa,Fs]=framepair('dgt',ga,gs,a,M);
%     [V,D]=framemuleigs(Fa,Fs,s,K);
%
%   Original help:
%   --------------
%
%   `gabmuleigs(K,c,g,a)` computes the *K* largest eigenvalues and eigen-
%   vectors of the Gabor multiplier with symbol *c* and time shift *a*.  The
%   number of channels is deduced from the size of the symbol *c*.  The
%   window *g* will be used for both analysis and synthesis.
%
%   `gabmuleigs(K,c,ga,gs,a)` does the same using the window the window *ga*
%   for analysis and *gs* for synthesis.
%
%   `gabmuleigs(K,c,a)` does the same using the a tight Gaussian window of
%   for analysis and synthesis.
%
%   If *K* is empty, then all eigenvalues/pairs will be returned.
%
%   `gabmuleigs` takes the following parameters at the end of the line of input
%   arguments:
%
%     'tol',t      Stop if relative residual error is less than the
%                  specified tolerance. Default is 1e-9
%
%     'maxit',n    Do at most n iterations.
%
%     'iter'       Call `eigs` to use an iterative algorithm.
%
%     'full'       Call `eig` to sole the full problem.
%
%     'auto'       Use the full method for small problems and the
%                  iterative method for larger problems. This is the
%                  default.
%
%     'crossover',c
%                  Set the problem size for which the 'auto' method
%                  switches. Default is 200.
%
%     'print'      Display the progress.
%
%     'quiet'      Don't print anything, this is the default.
%
%   See also: gabmul, dgt, idgt, gabdual, gabtight

warning(['LTFAT: GABMULEIGS has been deprecated, please use FRAMEMULEIGS ' ...
         'instead. See the help on FRAMEMULEIGS for more details.']);

% Change this to 1 or 2 to see the iterative method in action.
printopts=0;

if nargin<3
  error('%s: Too few input parameters.',upper(mfilename));
end;

if nargout==2
  doV=1;
else
  doV=0;
end;

M=size(c,1);
N=size(c,2);

istight=1;
if numel(p3)==1
  % Usage: h=gabmuleigs(c,K,a);
  a=p3;
  L=N*a;
  ga=gabtight(a,M,L);
  gs=ga;
  arglist=varargin;
else
  if numel(varargin{1})==1
    % Usage: h=gabmuleigs(c,K,g,a);
    ga=p3;
    gs=p3;
    a=varargin{1};
    L=N*a;
    arglist=varargin(2:end);
  else
    if numel(varargin{2})==1
      % Usage: h=gabmuleigs(c,K,ga,gs,a);
      ga=p3;
      gs=varargin{1};
      a =varargin{2};
      L=N*a;
      istight=0;
      arglist=varargin(3:end);
    end;
  end;
end;

definput.keyvals.maxit=100;
definput.keyvals.tol=1e-9;
definput.keyvals.crossover=200;
definput.flags.print={'quiet','print'};
definput.flags.method={'auto','iter','full'};

[flags,kv]=ltfatarghelper({},definput,arglist);

% Do the computation. For small problems a direct calculation is just as
% fast.

if (flags.do_iter) || (flags.do_auto && L>kv.crossover)

  if flags.do_print
    opts.disp=1;
  else
    opts.disp=0;
  end;
  opts.isreal = false;
  opts.maxit  = kv.maxit;
  opts.tol    = kv.tol;

  % Setup afun
  afun(1,c,ga,gs,a,M,L);

  if doV
    [V,D] = eigs(@afun,L,K,'LM',opts);
  else
    D     = eigs(@afun,L,K,'LM',opts);
  end;

else
  % Compute the transform matrix.
  bigM=tfmat('gabmul',c,ga,gs,a);

  if doV
    [V,D]=eig(bigM);
  else
    D=eig(bigM);
  end;

end;

% The output from eig and eigs is a diagonal matrix, so we must extract the
% diagonal.
D=diag(D);

% Sort them in descending order
[~,idx]=sort(abs(D),1,'descend');
D=D(idx(1:K));

if doV
  V=V(:,idx(1:K));
end;

% Clean the eigenvalues, if we know that they are real-valued
%if isreal(ga) && isreal(gs) && isreal(c)
%  D=real(D);
%end;

% The function has been written in this way, because Octave (at the time
% of writing) does not accept additional parameters at the end of the
% line of input arguments for eigs
function y=afun(x,c_in,ga_in,gs_in,a_in,M_in,L_in)
  persistent c;
  persistent ga;
  persistent gs;
  persistent a;
  persistent M;
  persistent L;

  if nargin>1
    c  = c_in;
    ga = ga_in;
    gs = gs_in;
    a  = a_in;
    M  = M_in;
    L  = L_in;
  else
    y=comp_idgt(c.*comp_dgt(x,ga,a,M,[0 1],0,0,0),gs,a,[0 1],0,0);
  end;


text:
function [ errors ] = test_gsp_solve_l1(  )

% Testin of gsp_solve_l1
gsp_reset_seed(0)
errors = 0;

N = 100;
Nf = 5;
G = gsp_sensor(N);
%G = gsp_2dgrid(sqrt(N));
G = gsp_compute_fourier_basis(G);

g = gsp_design_itersine(G,Nf);

s = gsp_filter_analysis(G,g{1},randn(N,1));
[alpha1, info1]  = gsp_solve_l1(G, g, s, 0 );

sp = gsp_filter_synthesis(G,g,alpha1);

if norm(sp-s)/norm(s)<1e-8
    fprintf('Test SOLVE_L1 OK\n')
else
    errors = errors +1;
    fprintf('Test SOLVE_L1 PAS OK\n')
    norm(sp-s)/norm(s)
end

[alpha2, info2]  = gsp_solve_l1(G, g, s, 0.01 );
sp = gsp_filter_synthesis(G,g,alpha2);
if norm(alpha2,1)<norm(alpha1,1)
    fprintf('Test SOLVE_L1 2 OK\n')
else
    errors = errors +1;
    fprintf('Test SOLVE_L1 2 PAS OK\n')
    norm(alpha2,1)
    norm(alpha1,1)
end

sp2 = gsp_filter_synthesis(G,g,alpha2);
if norm(sp2-s)/norm(s)<0.05
    fprintf('Test SOLVE_L1 3 OK\n')
else
    errors = errors +1;
    fprintf('Test SOLVE_L1 3 PAS OK\n')
    norm(sp2-s)/norm(s)
end

% old code... probably stupid
% % Testin of gsp_solve_l1
% gsp_reset_seed(0)
% errors = 0;
%
% N = 100;
% Nf = 5;
% param = struct;
% %param.distribute = 1;
% G = gsp_sensor(N,param);
% %G = gsp_2dgrid(sqrt(N));
% G = gsp_compute_fourier_basis(G);
% G = gsp_spectrum_cdf_approx(G);
% paramfilter.filter = gsp_design_meyer(G,Nf);
% g = gsp_design_warped_translates(G, Nf, paramfilter);
% %g = gsp_design_itersine(G,Nf);
%
% % figure
% % gsp_plot_filter(G,g)
%
% alphain = zeros(N,Nf);
%
% for ii = 1:4
%     alphain(ii,ii) = 1;
% end
%
% alphain = gsp_mat2vec(alphain);
% signal = gsp_filter_synthesis(G,g,alphain);
%
% % figure
% % gsp_plot_signal(G, signal);
% param.verbose = 0;
% param.maxit = 1000;
% % fig=figure(100);
% % param.do_sol=@(x) plot_objective(x,fig);
%
% param.tol = 1e-6;
% param.guess = alphain;
% lambda = 0.05;
% alphaout = gsp_solve_l1(G,g,signal,lambda,param);
%
% % param.guess = alphaout;
% % lambda2 = 0.1;
% % param.tol = 1e-10;
% % param.maxit = 200;
% % alphaout = gsp_solve_l0(G,g,signal,lambda2,param);
% % norm(alphaout-alphain)
%
% if norm(alphaout-alphain)<2
%     fprintf('Test SOLVE_L1 OK\n')
% else
%     errors = errors +1;
%     fprintf('Test SOLVE_L1 PAS OK\n')
%     norm(alphaout-alphain)
% end

end


text:
function b = cc_mv ( m, n, ncc, icc, ccc, acc, x )

%*****************************************************************************80
%
%% CC_MV multiplies a CC matrix by a vector
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992
%
%  Parameters:
%
%    Input, integer M, the number of rows.
%
%    Input, integer N, the number of columns.
%
%    Input, integer NCC, the number of CC values.
%
%    Input, integer ICC(NCC), the CC rows.
%
%    Input, integer CCC(N+1), the compressed CC columns
%
%    Input, real ACC(NCC), the CC values.
%
%    Input, real X(N), the vector to be multiplied.
%
%    Output, real B(M), the product A*X.
%
  b = zeros(m,1);

  for j = 1 : n
    for k = ccc(j) : ccc(j+1) - 1
      i = icc(k);
      b(i) = b(i) + acc(k) * x(j);
    end
  end

  return
end


text:
function a = sylvester_kac_inverse ( n )

%*****************************************************************************80
%
%% SYLVESTER_KAC_INVERSE returns the inverse of the SYLVESTER_KAC matrix.
%
%  Example:
%
%    N = 6:
%
%      0     1/5    0  -2/15  0   8/15
%      1      0     0    0    0    0
%      0      0     0   1/3   0  -4/3
%    -4/3     0    1/3   0    0    0
%      0      0     0    0    0    1
%     8/15    0   -2/15  0   1/5   0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  if ( mod ( n, 2 ) == 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SYLVESTER_KAC_INVERSE - Fatal error!\n' );
    fprintf ( 1, '  The matrix is singular for odd N.\n' );
    error ( 'SYLVESTER_KAC_INVERSE - Fatal error!' );
  end

  a = zeros ( n, n );

  for i = 1 : n

    if ( mod ( i, 2 ) == 1 )

      for j = i : 2 : n-1

        if ( j == i )
          p1 = 1.0 / ( n - j );
          p2 = 1.0 /       j;
        else
          p1 = - p1 * (     j - 1 ) / ( n - j );
          p2 = - p2 * ( n - j + 1 ) /       j;
        end

        a(i,j+1) = p1;
        a(j+1,i) = p2;

      end

    end

  end

  return
end


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

%%
clear; close all; clc;

%% Syntax
% |[V]=selectBox(varargin);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function ci =  cusum(cat)
    % This function calculates the CUMSUm function (Page 1954).
    %

    report_this_filefun(mfilename('fullpath'));

    [existFlag,figNumber]=figure_exists('CUSUM',1); %find figure with name 'CUSUM', and do not pop to foreground

    if existFlag
        cfig = figNumber;
    else
        cfig=figure_w_normalized_uicontrolunits(...                  %build figure for plot
            'Units','normalized','NumberTitle','off',...
            'Name','CUSUM',...
            'MenuBar','none',...
            'visible','off',...
            'pos',[ 0.300  0.3 0.4 0.6]);
        ho=false;

        matdraw
    end   % if fig exist

    m  = cat(:,6);
    me = mean(m);
    i = (1:1:length(m));
    ci = cumsum(m)' - i.*me;

    figure_w_normalized_uicontrolunits(cfig)
    delete(gca);delete(gca);
    plot(cat(:,3),ci,'o')
    %plot(i,ci,'o')
    set(gca,'visible','on','FontSize',10,'FontWeight','normal',...
        'FontWeight','normal','LineWidth',1.0,...
        'Box','on')
    xlabel('Time [yrs]')
    ylabel('CUSUM [yrs]')
    grid

author: CelsoReyes


text:
function z = makebw(or,ph)
% function z = makebw(or,ph)
%
% Create BWT wavelets
% Version 1.2
%
% Arguments:
%  or: Orientation { -1 (DC), 0, 45, 90, 135 }
%  ph: Phase { 0 (odd), 90 (even) }
%
% Result:
%  z: The selected wavelet, or NaN if the parameters are not recognized
%
% Citation:
%  Willmore B, Prenger RJ, Wu MC and Gallant JL (2008). The Berkeley
%  Wavelet Transform: A biologically-inspired orthogonal wavelet transform.
%  Neural Computation 20:6, 1537-1564
%
% The article is available at:
%  <http://dx.doi.org/10.1162/neco.2007.05-07-513>
%
% Copyright (c) 2008 Ben Willmore
%
% Permission is hereby granted, free of charge, to any person
% obtaining a copy of this software and associated documentation
% files (the "Software"), to deal in the Software without
% restriction, including without limitation the rights to use,
% copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the
% Software is furnished to do so, subject to the following
% conditions:
%
% The above copyright notice and this permission notice shall be
% included in all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
% EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
% OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
% NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
% HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
% WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
% FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
% OTHER DEALINGS IN THE SOFTWARE.

z = nan;

if (or==-1)
    z = ones(3,3)/sqrt(9);

elseif (or==0)
    if (ph==0)
        z = [-1 0 1; -1 0 1; -1 0 1]/sqrt(6);
    elseif (ph==90)
        z = [-1 2 -1; -1 2 -1; -1 2 -1]/sqrt(18);
    end

elseif (or==45)
    if (ph==0)
        z = [-1 1 0; 1 0 -1; 0 -1 1]/sqrt(6);
    elseif (ph==90)
        z = [-1 -1 2; -1 2 -1; 2 -1 -1]/sqrt(18);
    end

elseif (or==90)
    if (ph==0)
        z = [-1 -1 -1; 0 0 0; 1 1 1]/sqrt(6);
    elseif (ph==90)
        z = [-1 -1 -1; 2 2 2; -1 -1 -1]/sqrt(18);
    end

elseif (or==135)
    if (ph==0)
        z = [0 -1 1; 1 0 -1; -1 1 0]/sqrt(6);
    elseif (ph==90)
        z = [2 -1 -1; -1 2 -1; -1 -1 2]/sqrt(18);
    end

end

if (isnan(z))
    disp('Wavelet parameters not recognized');
end

author: Sable

repo: mcbench-benchmarks


text:
% Calculate the number of independent cycles (use L = m-n+c)
%      where L = num cycles, m - num edges, n - num nodes,
%                      c - number of connected components
% This is also known as the "cyclomatic number": the number of edges
% that needs to be removed so that the graph doesn't have cycles.
%
% INPUTS: adjacency matrix, nxn
% OUTPUTs: number of independent cycles (or cyclomatic number)
%
% Other routines used: numNodes.m, numEdges.m, findConnComp.m
% GB: last updated, Oct 5 2012

function L = numCycles(adj)

n=numNodes(adj);
m=numEdges(adj);
comp_mat = findConnComp(adj);

L=m-n+length(comp_mat);

author: aeolianine

repo: octave-networks-toolbox


text:
function [surf, labels] = imSurfaceArea(img, varargin)
% Surface area of the regions within a 3D binary or label image.
%
%   S = imSurfaceArea(IMG)
%   Estimates the surface area of the 3D binary structure represented by
%   IMG.
%
%   S = imSurfaceArea(IMG, NDIRS)
%   Specifies the number of directions used for estimating surface area.
%   NDIRS can be either 3 or 13, default is 13.
%
%   S = imSurfaceArea(..., SPACING)
%   Specifies the spatial calibration of the image. SPACING is a 1-by-3 row
%   vector containing the voxel size in the X, Y and Z directions, in that
%   orders.
%
%   S = imSurfaceArea(LBL)
%   [S, L] = imSurface(LBL)
%   When LBL is a label image, returns the surface area of each label in
%   the 3D array, and eventually returns the indices of processed labels.
%
%   S = imSurfaceArea(..., LABELS)
%   In the case of a label image, specifies the labels of the region to
%   analyse.
%
%
%   Example
%     % Create a binary image of a ball
%     [x y z] = meshgrid(1:100, 1:100, 1:100);
%     img = sqrt( (x-50.12).^2 + (y-50.23).^2 + (z-50.34).^2) < 40;
%     % compute surface area of the ball
%     S = imSurfaceArea(img)
%     S =
%         2.0103e+04
%     % compare with theoretical value
%     Sth = 4*pi*40^2;
%     100 * (S - Sth) / Sth
%     ans =
%         -0.0167
%
%     % compute surface area of several regions in a label image
%     img = uint8(zeros(10, 10, 10));
%     img(2:3, 2:3, 2:3) = 1;
%     img(5:8, 2:3, 2:3) = 2;
%     img(5:8, 5:8, 2:3) = 4;
%     img(2:3, 5:8, 2:3) = 3;
%     img(5:8, 5:8, 5:8) = 8;
%     [surfs, labels] = imSurfaceArea(img)
%     surfs =
%        16.4774
%        29.1661
%        29.1661
%        49.2678
%        76.7824
%     labels =
%          1
%          2
%          3
%          4
%          8
%
%
%   See also
%     imVolume, imMeanBreadth, imSurfaceAreaDensity, imJointSurfaceArea
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2010-07-26,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRAE - Cepia Software Platform.

%% Parse input arguments

% check image dimension
if ndims(img) ~= 3
    error('first argument should be a 3D image');
end

% the labels to compute
labels = [];

% default number of directions
nDirs = 13;

% default image calibration
delta = [1 1 1];

% methods to compute direction weights. Can be {'Voronoi'}, 'isotropic'.
directionWeights = 'voronoi';

% Process user input arguments
while ~isempty(varargin)
    var1 = varargin{1};

    if isnumeric(var1)
        % option is either connectivity or resolution
        if isscalar(var1)
            nDirs = var1;
        elseif all(size(var1) == [1 3])
            delta = var1;
        elseif size(var1, 2) == 1
            labels = var1;
        end
        varargin(1) = [];

    elseif ischar(var1)
        if length(varargin) < 2
            error('optional named argument require a second argument as value');
        end
        if strcmpi(var1, 'directionweights')
            directionWeights = varargin{2};
        end
        varargin(1:2) = [];

    else
        error('option should be numeric');
    end

end

%% Process label image

% in case of a label image, return a vector with a set of results
if ~islogical(img)
    % extract labels if necessary (considers 0 as background)
    if isempty(labels)
        labels = imFindLabels(img);
    end

    % allocate result array
    nLabels = length(labels);
    surf = zeros(nLabels, 1);

    % compute bounding box of each region
    boxes = imBoundingBox(img, labels);

    % Compute surface area of each region considered as binary image
    % The computation is performed on a subset of the image for reducing
    % memory footprint.
    for i = 1:nLabels
        label = labels(i);

        % convert bounding box to image extent, in x, y and z directions
        box = boxes(i,:);
        i0 = ceil(box([3 1 5]));
        i1 = floor(box([4 2 6]));

        % crop image of current label
        bin = img(i0(1):i1(1), i0(2):i1(2), i0(3):i1(3)) == label;
        surf(i) = imSurfaceArea(bin, nDirs, delta, 'directionWeights', directionWeights);
    end

    return;
end

%% Process binary images

% in case of binary image, compute only one label...
labels = 1;

% distances between a pixel and its neighbours.
d1  = delta(1); % x
d2  = delta(2); % y
d3  = delta(3); % z

% volume of a voxel (used for computing line densities)
vol = d1 * d2 * d3;

%% Main processing for 3 directions

% number of voxels
nv = sum(img(:));

% number of connected components along the 3 main directions
% (Use Graph-based formula: chi = nVertices - nEdges)
n1 = nv - sum(sum(sum(img(:,1:end-1,:) & img(:,2:end,:)))); % x
n2 = nv - sum(sum(sum(img(1:end-1,:,:) & img(2:end,:,:)))); % y
n3 = nv - sum(sum(sum(img(:,:,1:end-1) & img(:,:,2:end)))); % z

if nDirs == 3
    % compute surface area by averaging over the 3 main directions
    surf = 4/3 * (n1/d1 + n2/d2 + n3/d3) * vol;
    return;
end

%% Additional processing for 13 directions

% Number of connected components along diagonals contained in the three
% main planes
% XY planes
n4 = nv - sum(sum(sum(img(2:end,1:end-1,:)   & img(1:end-1,2:end,:))));
n5 = nv - sum(sum(sum(img(1:end-1,1:end-1,:) & img(2:end,2:end,:))));
% XZ planes
n6 = nv - sum(sum(sum(img(:,2:end,1:end-1)   & img(:,1:end-1,2:end))));
n7 = nv - sum(sum(sum(img(:,1:end-1,1:end-1) & img(:,2:end,2:end))));
% YZ planes
n8 = nv - sum(sum(sum(img(2:end,:,1:end-1)   & img(1:end-1,:,2:end))));
n9 = nv - sum(sum(sum(img(1:end-1,:,1:end-1) & img(2:end,:,2:end))));

%TODO: add the case of 9 directions ?

% Number of connected components along lines corresponding to diagonals of
% the unit cube
n10 = nv - sum(sum(sum(img(1:end-1,1:end-1,1:end-1) & img(2:end,2:end,2:end))));
n11 = nv - sum(sum(sum(img(2:end,1:end-1,1:end-1) & img(1:end-1,2:end,2:end))));
n12 = nv - sum(sum(sum(img(1:end-1,2:end,1:end-1) & img(2:end,1:end-1,2:end))));
n13 = nv - sum(sum(sum(img(2:end,2:end,1:end-1) & img(1:end-1,1:end-1,2:end))));

% space between 2 voxels in each direction
d12  = hypot(d1, d2);
d13  = hypot(d1, d3);
d23  = hypot(d2, d3);
d123 = sqrt(d1^2 + d2^2 + d3^2);

% Compute weights corresponding to surface fraction of spherical caps
% For isotropic case, weights correspond to:
% c1 = 0.04577789120476 * 2;  % Ox
% c2 = 0.04577789120476 * 2;  % Oy
% c3 = 0.04577789120476 * 2;  % Oz
% c4 = 0.03698062787608 * 2;  % Oxy
% c6 = 0.03698062787608 * 2;  % Oxz
% c8 = 0.03698062787608 * 2;  % Oyz
% c10 = 0.03519563978232 * 2;  % Oxyz
if strcmp(directionWeights, 'isotropic')
    c = zeros(13,1);
    c(1) = 0.04577789120476 * 2;  % Ox
    c(2) = 0.04577789120476 * 2;  % Oy
    c(3) = 0.04577789120476 * 2;  % Oz
    c(4:5)  = 0.03698062787608 * 2;  % Oxy
    c(6:7)  = 0.03698062787608 * 2;  % Oxz
    c(8:9)  = 0.03698062787608 * 2;  % Oyz
    c(10:13) = 0.03519563978232 * 2;  % Oxyz

else
    c = computeDirectionWeights3d13(delta);
end

% compute the weighted sum of each direction
% intersection count * direction weight / line density
surf = 4 * vol * (...
    n1*c(1)/d1 + n2*c(2)/d2 + n3*c(3)/d3 + ...
    (n4+n5)*c(4)/d12 + (n6+n7)*c(6)/d13 + (n8+n9)*c(8)/d23 + ...
    (n10 + n11 + n12 + n13)*c(10)/d123 );

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
 function[M1,M2,M3] = integrate_vf(k,w,f,v)
      M1   = (k*sum(v .* w .* f));    % Number density
      M2   = (k*sum((v .^2 ) .* w .* f));   % Macrospic moment in x
      M3   = (k*sum(1/2*( v .*abs(v).^2 ).* w .* f)); % Energy Density
 end

repo: Aero-matlab


text:
function stroud_test28 ( )

%*****************************************************************************80
%
%% TEST28 tests SIMPLEX_VOLUME_ND.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST28\n' );
  fprintf ( 1, '  SIMPLEX_VOLUME_ND computes the volume of a simplex\n' );
  fprintf ( 1, '    in N dimensions.\n' );
  fprintf ( 1, '\n' );

  for n = 2 : 4

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Spatial dimension N = %d\n', n );
%
%  Set the values of the simplex.
%
    v = setsim ( n );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Simplex vertices:\n' );
    fprintf ( 1, '\n' );

    for i = 1 : n+1
      for j = 1 : n
        fprintf ( 1, '  %4f', v(i,j) );
      end
      fprintf ( 1, '\n' );
    end

    volume = simplex_volume_nd ( n, v );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Volume is %f\n', volume );

  end

  return
end


text:
function dat = read_egis_data(filename, hdr, begtrial, endtrial, chanindx);

% READ_EGIS_DATA reads the data from an EGI EGIS format file
%
% Use as
%   dat = read_egis_data(filename, hdr, begtrial, endtrial, chanindx);
% where
%   filename       name of the input file
%   hdr            header structure, see READ_HEADER
%   begtrial       first trial to read, mutually exclusive with begsample+endsample
%   endtrial       last trial to read,  mutually exclusive with begsample+endsample
%   chanindx       list with channel indices to read
%
% This function returns a 3-D matrix of size Nchans*Nsamples*Ntrials.
% Note that EGIS session files are defined as always being epoched.
% For session files the trials are organized with the members of each cell grouped
% together.  For average files the "trials" (subjects) are organized with the cells
% also grouped together (e.g., "cell1sub1, cell1sub2, ...).
%_______________________________________________________________________
%
%
% Modified from EGI's EGI Toolbox with permission 2007-06-28 Joseph Dien

% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: read_egis_data.m 3063 2011-03-07 05:02:52Z josdie $

fh=fopen([filename],'r');
if fh==-1
    error('wrong filename')
end
fclose(fh);

[fhdr,chdr,ename,cnames,fcom,ftext] = read_egis_header(filename);
fh=fopen([filename],'r');
fhdr(1)=fread(fh,1,'int32'); %BytOrd
[str,maxsize,cEndian]=computer;
if fhdr(1)==16909060
    if cEndian == 'B'
        endian = 'ieee-be';
    elseif cEndian == 'L'
        endian = 'ieee-le';
    end;
elseif fhdr(1)==67305985
    if cEndian == 'B'
        endian = 'ieee-le';
    elseif cEndian == 'L'
        endian = 'ieee-be';
    end;
else
    error('This is not an EGIS average file.');
end;

if fhdr(2) == -1
    fileType = 'ave';
elseif fhdr(2) == 3
    fileType = 'ses';
else
    error('This is not an EGIS file.');
end;

dat=zeros(hdr.nChans,hdr.nSamples,endtrial-begtrial+1);

%read to end of file header
status=fseek(fh,(130+(2*fhdr(18))+(4*fhdr(19))),'bof');
status=fseek(fh,2,'cof');
for loop=1:fhdr(18)
  temp=fread(fh,80,'uchar');
  theName=strtok(temp);
  theName=strtok(theName,char(0));
  cnames{loop}=deblank(char(theName))';
  status=fseek(fh,fhdr(24+(loop-1))-80,'cof');
end
fcom=fread(fh,fhdr(20),'uchar');
ftext=fread(fh,fhdr(21),'uchar');
fpad=fread(fh,fhdr(22),'uchar');
status=fseek(fh,-2,'cof');

%read to start of desired data
fseek(fh, ((begtrial-1)*hdr.nChans*hdr.nSamples*2), 'cof');

for segment=1:(endtrial-begtrial+1)
    dat(:,:,segment) = fread(fh, [hdr.nChans, hdr.nSamples],'int16',endian);
end
dat=dat(chanindx, :,:);

if fileType == 'ave'
    dat=dat/fhdr(12); %convert to microvolts
elseif fileType == 'ses'
    dat=dat/5; %convert to microvolts (EGIS sess files created by NetStation use a 5 bins per microvolt scaling)
end;

fclose(fh);

author: PatternRecognition

repo: OpenBMI


text:
function spm_run_bms_vis(varargin)
% Show results for Bayesian Model Selection Maps
% SPM job execution function
% takes a harvested job data structure (or no input) and calls SPM
% functions to show results from Bayesian Model Selection of
% Log. Evidence Maps
%
% Input:
% Varargin - can be harvested job data structure (see matlabbatch help).
% Output:
% No output.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Maria Joao Rosa
% $Id: spm_run_bms_vis.m 7577 2019-04-24 08:59:56Z guillaume $

% Input
% -------------------------------------------------------------------------
if isempty(varargin)
   job.file{1}=''; job.img{1}=''; job.thres = []; job.scale = []; job.k=[];
else
   job = varargin{1};
end
file_str = length(job.file{1});
image    = length(job.img{1});
thres_b  = length(job.thres);
odds_b   = length(job.scale);
ext_thre = length(job.k);

% Initialise SPM Interactive window
% -------------------------------------------------------------------------
[Finter,F,CmdLine] = spm('FnUIsetup','SPM BMC Toolbox: Visualise options',0);
F = spm_figure('FindWin','Graphics');
spm('Clear',Finter);
WS = spm('WinScale');
FS = spm('FontSizes');

% Select results image from BMS Maps
% -------------------------------------------------------------------------
if file_str
   file_fname = job.file{1};
else
   file_fname = spm_select(1,'mat','Load BMS file (.mat)');
end

% Load BMS
load(file_fname);
wd = fileparts(file_fname);

% Select results image from BMS Maps
% -------------------------------------------------------------------------
if image
   post = job.img{1};
else
   post = spm_select(1,'image','Select BMS Maps results (ex: ppm, alpha or epm image)',[],wd);
end

% Select threshold to apply to image
% -------------------------------------------------------------------------
if thres_b
   threshold = job.thres;
else
   threshold = spm_input('Probability Threshold:', '+1', 'r', '0.5', [1, inf]);
end

% Extent threshold
% -------------------------------------------------------------------------
if ext_thre
   k = job.k;
else
   k = spm_input('& extent threshold {voxels}','+1','r',0,1,[0,Inf]);
end

% Select scale
% -------------------------------------------------------------------------
if odds_b
   odds = job.scale;
else
   odds = spm_input('Change scale to Log-Odds?', '+1', 'y/n', [1,0], 2);
end

% Image dimensions
% -------------------------------------------------------------------------
V             = spm_vol(post);
M             = V.mat;
iM            = inv(M);
vox           = sqrt(sum(M(1:3,1:3).^2));
DIM           = V.dim(1:3)';
xdim          = DIM(1); ydim  = DIM(2); zdim  = DIM(3);
[xords,yords] = ndgrid(1:xdim,1:ydim);
xords         = xords(:)';  yords = yords(:)';
I             = 1:xdim*ydim;
zords_init    = ones(1,xdim*ydim);
svol          = xdim*ydim*zdim;

% Legend of results being plotted
% -------------------------------------------------------------------------
[pathstr,name_image] = fileparts(V.fname);
undersc              = find(name_image=='_');

if ~isempty(undersc) && length(undersc)>1
    res_name             = name_image(undersc(end)+1);
    subset_model         = name_image(1:undersc(end-1)-1);
else
    res_name = '';
end

% Show results being displayed on graphics window
switch res_name
       case 'a'
            xSPM.str = sprintf('Alphas: %s',subset_model);
       case 'p'
            xSPM.str = sprintf('PPM: %s',subset_model);
       case 'x'
            xSPM.str = sprintf('PPM: %s',subset_model);
       case 'e'
            xSPM.str = sprintf('EPM: %s',subset_model);
       otherwise
            xSPM.str = '';
end

% Loop over slices (get voxels above threshold)
% -------------------------------------------------------------------------
xyz_above = [];
z_above   = [];
xyz_total = [];

for z = 1:zdim,

    zords = z*zords_init;
    xyz   = [xords(I); yords(I); zords(I)];
    zvals = spm_get_data(V,xyz);
    above = find(zvals>threshold);
    if ~isempty(above)
        xyz_above = [xyz_above,xyz(:,above)];
        z_above   = [z_above,zvals(above)];
    end
    xyz_total = [xyz_total,xyz];
end

% SPM figure: SPM, MIP and GUI
% -------------------------------------------------------------------------
if ~isempty(z_above)

    z_orig = z_above;
    % Log odds transform
    if odds
       z_odds = log(z_above./(ones(1,length(z_above))-z_above));
       if isempty(find(z_odds == Inf))  % Data don't contain Infs
          z_above = z_odds;
        else
          odds = 0;                     % Don't plot odds
        end
    end

    % Calculate extent threshold filtering
    % ---------------------------------------------------------------------
    A     = spm_clusters(xyz_above);
    Q     = [];
    for i = 1:max(A)
        j = find(A == i);
        if length(j) >= k; Q = [Q j]; end
    end

    % ...eliminate voxels
    %----------------------------------------------------------------------
    Z     = z_above(:,Q);
    XYZ   = xyz_above(:,Q);
    z_ps  = z_orig(:,Q);
    if isempty(Q)
        warning(sprintf('No voxels survive extent threshold k=%0.2g',k))
    end

    % Save data in xSPM structure
    voxels_mm         = M*[XYZ;ones(1,size(XYZ,2))];
    voxels_mm         = voxels_mm(1:3,:);
    xSPM.swd          = file_fname;
    xSPM.STAT         = '';
    xSPM.Z            = Z;
    xSPM.XYZ          = XYZ;
    xSPM.XYZmm        = voxels_mm;
    xSPM.xVol.M       = M;
    xSPM.xVol.DIM     = DIM;
    xSPM.SPM.xVol.XYZ = xyz_total;
    xSPM.SPM.xVol.S   = svol;
    xSPM.SPM.xVol.M   = M;
    xSPM.SPM.xVol.DIM = DIM;
    xSPM.DIM          = DIM;
    xSPM.M            = M;
    xSPM.iM           = iM;
    xSPM.n            = 1;
    xSPM.k            = k;
    xSPM.df           = [];
    xSPM.u            = threshold;
    xSPM.FWHM         = 10;                    % arbritary value
    xSPM.STAT         = '';
    xSPM.VOX          = vox;
    xSPM.thres        = threshold;
    xSPM.vols         = post;
    xSPM.scale        = odds;
    xSPM.units        = {'mm'  'mm'  'mm'};
    xSPM.Ps           = [];
    xSPM.S            = 0;
    xSPM.z_ps         = z_ps;
    xSPM.VRpv         = [];
    xSPM.Vspm         = [];
    BMS.xSPM          = xSPM;

    % Display on workspace
    assignin('base','BMS',BMS);
    assignin('base','xSPM',xSPM);

    % Display results
    hReg = spm_bms_display(BMS,'Init');

    assignin('base','hReg',hReg);

else

    % No voxels for selected threshold
    disp('No voxels above threshold!');
    uicontrol(Finter,'Style','Text','String','No voxels above threshold!',...
            'Position',[020 010 200 020].*WS,...
            'FontAngle','Italic',...
            'FontSize',FS(12),...
            'HorizontalAlignment','Left',...
            'ForegroundColor','w')
    do_results = spm_input('Return to Results?', '+1', 'y/n', [1,0], 2);

    % Go back to BMS Maps (Results)
    if do_results

        % Return to display routine
        spm_run_bms_vis;

    else

        spm('Clear',Finter); % Clear Finter window
        return               % Finish

    end

end


text:
function [x,D_struc,problem,r,res,solvertime,prob] = call_mosek_geometric(model);

x = [];
D_struc = [];
r = [];
res = [];
solvertime = 0;
[prob,problem] = yalmip2geometric(model.options,model.F_struc,model.c,model.Q,model.K,model.ub,model.lb,model.monomtable,model.linear_variables,model.extended_variables);
if problem == 0

    % Mosek does not support equalities
    if ~isempty(prob.G)
        prob.A = [prob.A;prob.G;-prob.G];
        prob.b = [prob.b;prob.h;1./prob.h];
        prob.map = [prob.map;max(prob.map) + (1:2*length(prob.h))'];
    end
    if model.options.savedebug
        save mosekdebug prob
    end

    param = model.options.mosek;

    % Call MOSEK
    showprogress('Calling MOSEK',model.options.showprogress);
    if model.options.verbose == 0
        solvertime = tic;
        res = mskgpopt(prob.b,prob.A,prob.map,param,'minimize echo(0)');
        solvertime = toc(solvertime);
    else
        solvertime = tic;
    	res = mskgpopt(prob.b,prob.A,prob.map,param,'minimize');
        solvertime = toc(solvertime);
    end
    sol = res.sol;

    x = zeros(length(model.c),1);
    x(model.linear_variables)=exp(res.sol.itr.xx);
    D_struc = [];

    % Check, currently not exhaustive...
    switch res.sol.itr.prosta
        case 'PRIMAL_AND_DUAL_FEASIBLE'
            problem = 0;
        case 'PRIMAL_INFEASIBLE'
            problem = 1;
        case 'DUAL_INFEASIBLE'
            problem = 2;
        case 'UNKNOWN'
            problem = 9;
        otherwise
            problem = -1;
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [price, CI]  = blsSobol(S,E,r,T,sigma,nSims)
%
%GetSobolVanillaPrice - Vanilla option pricing using simulation and Sobol
% generator
%
% Return the pice of a vanilla option and the standard deviation of the
% price simulated
%
% Inputs:
%   S   	- Current price of the underlying asset.
%
%   E        - Strike (i.e., exercise) price of the option.
%
%   r        - Annualized continuously compounded risk-free rate of return
%                 over the life of the option, expressed as a positive decimal
%                 number.
%
%   T        - Time to expiration of the option, expressed in years.
%
%   sigma    - Annualized asset price volatility (i.e., annualized standard
%                 deviation of the continuously compounded asset return),
%                 expressed as a positive decimal number.
%
%
%   divYield  - Annualized continuously compounded yield of the underlying
%                 asset over the life of the option, expressed as a decimal
%                 number. If Yield is empty or missing. the default value is
%                 zero.
%
%                 For example, this could represent the dividend yield (annual
%                 dividend rate expressed as a percentage of the price of the
%                 security) or foreign risk-free interest rate for options
%                 written on stock indices and currencies, respectively.
%   nSims      - Number of Simulation used for the pricing
%   nSteps     - Number of time steps used to simulate
%  [SobolPrice,stdSobol] =   GetSobolVanillaPrice(S,E,r,T,sigma,divYield,nsim,nSteps);

Dt = T;

%Generate the random numbers using SOBOL sequences

% Sobol sequences have some zeros
% a common approach in the litterature is to suppress the 64 first points
% the sobol generator has been found on the web

P = sobolset(1);
SobolRandomNumbers = net(P,nSims);

% Sobol numbers are between 0 and  1
% We need to get a normal distribution from this pseudo uniform drawing

RandomNumbers = norminv(SobolRandomNumbers');
mat = exp( (r-sigma^2/2)*Dt + sigma*sqrt(Dt).*RandomNumbers  );
mat = cumprod(mat , 1);
mat = mat.*S;

% Discount and calculate the option price

V = exp(-r*T) * max(mat(end,:)-E , 0);
[price,VarParice,CI]= normfit(V);


text:
function z = quad_pos_over_lin( x, y, dim )

%QUAD_POS_OVER_LIN Sum of squares of positives over linear.
%   Z=QUAD_POS_OVER_LIN(X,Y), where X is a vector and Y is a scalar, is equal to
%   SUM(MAX(X,0).^2)./Y if Y is positive, and +Inf otherwise. Both X and Y must
%   be real.
%
%   For matrices, QUAD_POS_OVER_LIN(X,Y) is a row vector containing the
%   application of QUAD_POS_OVER_LIN to each column. For N-D arrays, the
%   operation is applied to the first non-singleton dimension of X.
%
%   QUAD_POS_OVER_LIN(X,Y,DIM) takes the sum along the dimension DIM of X.
%   A special value of DIM == 0, is accepted here, which is automatically
%   replaced with DIM == NDIMS(X) + 1. This has the effect of eliminating
%   the sum; thus QUAD_POS_OVER_LIN( X, Y, NDIMS(X) + 1 ) = MAX(X,0).^2./Y.
%
%   In all cases, both X and Y must be real, and Y must either be a scalar
%   or a matrix of the same size as SUM(X,DIM).
%
%   Disciplined convex programming information:
%       QUAD_POS_OVER_LIN is convex, nondecreasing in X, and nonincreasing
%       in Y. Thus when used with CVX expressions, X must be convex (or
%       affine) and Y must be concave (or affine).

%
% Check arguments
%

error( nargchk( 2, 3, nargin ) );
if ~isreal( x ),
    error( 'First argument must be real.' );
elseif ~isreal( y ),
    error( 'Second argument must be real.' );
elseif nargin < 3 || isempty( dim ),
    dim = cvx_default_dimension( size( x ) );
elseif ~cvx_check_dimension( dim ),
    error( 'Third argument, if supplied, must be a positive integer.' );
end

%
% Perform calculation
%

z = quad_over_lin( max( x, 0 ), y, dim );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function [res min_cs] = tracking_dp(dres, c_en, c_ex, c_ij, betta, thr_cost, max_it, nms_in_loop)

if ~exist('max_it')
  max_it = 1e5;
end
if ~exist('thr_cost')
  thr_cost = 0;
end

thr_nms = 0.5;

dnum = length(dres.x);

dres.c = betta - dres.r;

dres.dp_c     = [];
dres.dp_link  = [];
dres.orig     = [];

min_c     = -inf;
it        = 0;
k         = 0;
inds_all  = zeros(1,1e5);
id_s      = zeros(1,1e5);
redo_nodes = [1:dnum]';
while (min_c < thr_cost) && (it < max_it)
  it = it+1;

  dres.dp_c(redo_nodes,1) = dres.c(redo_nodes) + c_en;
  dres.dp_link(redo_nodes,1) = 0;
  dres.orig(redo_nodes,1) = redo_nodes;

  for ii=1:length(redo_nodes)
    i = redo_nodes(ii);
    f2 = dres.nei(i).inds;
    if isempty(f2)
      continue
    end

    [min_cost j] = min(c_ij + dres.c(i) + dres.dp_c(f2));
    min_link = f2(j);
    if dres.dp_c(i,1) > min_cost
      dres.dp_c(i,1) = min_cost;
      dres.dp_link(i,1) = min_link;
      dres.orig(i,1) = dres.orig(min_link);
    end
  end

  [min_c ind] = min(dres.dp_c + c_ex);

  inds = zeros(dnum,1);

  k1 = 0;
  while ind~=0
    k1 = k1+1;
    inds(k1) = ind;
    ind = dres.dp_link(ind);
  end
  inds = inds(1:k1);

  inds_all(k+1:k+length(inds)) = inds;
  id_s(k+1:k+length(inds)) = it;
  k = k+length(inds);

  if nms_in_loop
    supp_inds = nms_aggressive(dres, inds, thr_nms);
    origs = unique(dres.orig(supp_inds));
    redo_nodes = find(ismember(dres.orig, origs));
  else
    supp_inds = inds;
    origs = inds(end);
    redo_nodes = find(dres.orig == origs);
  end
  redo_nodes = setdiff(redo_nodes, supp_inds);
  dres.dp_c(supp_inds) = inf;
  dres.c(supp_inds) = inf;

  min_cs(it) = min_c;
end
inds_all = inds_all(1:k);
id_s = id_s(1:k);

res = sub(dres, inds_all);
res.id = id_s';

author: yuxng


text:
function [audio_obj, features] = get_files_and_mfcc_features(input_dir, l, fn_filter)
    current_dir = cd();
    cd(input_dir);
    files = get_files_by_filter(fn_filter);
    % label according to l-th character in file name
    audio_obj = miraudio(files, 'Label', l, 'Normal');
    features = struct();
    mfcc = mirmfcc(audio_obj, 'Frame');
    features.mfcc = mirstat(mfcc);
    mfcc_delta = mirmfcc(audio_obj, 'Frame', 'Delta');
    features.mfcc_delta = mirstat(mfcc_delta);
    centroid = mircentroid(audio_obj, 'Frame');
    features.centroid = mirstat(centroid);
    rms = mirrms(audio_obj, 'Frame');
    features.rms = mirstat(rms);
    cd(current_dir);
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function degrees = degrees(radians)

% DEGREES (RADIANS)
%
% Conversion function from Radians to Degrees.
% Richard Medlock 12-03-2002
%
% Last Updated: 18-09-2009
% - Calculation simplified to make it more efficient
%   based on a suggestion by Joel Parker.

% Original calculation
% radians = radians/((2*pi)/360);

degrees = radians/(pi/180);


text:
% StackExchange Signal Processing Q61273
% https://dsp.stackexchange.com/questions/61273
% Estimate Noise Variance Given Multiple Realization of the Same Image
% References:
%   1.  A
% Remarks:
%   1.  B
% TODO:
% 	1.  C
% Release Notes
% - 1.0.000     15/10/2019
%   *   First release.

%% General Parameters

subStreamNumberDefault = 12;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

inputImageFileName  = 'InputImage.png';

vNumRealizations    = 2:60;
noiseStd            = 0.1;

%% Generate Data

mI = im2double(imread(inputImageFileName));

numRows = size(mI, 1);
numCols = size(mI, 2);

tI              = mI + (noiseStd * randn(numRows, numCols, vNumRealizations(end)));
mMeanImage      = zeros(numRows, numCols);
tN              = zeros(numRows, numCols, vNumRealizations(end));
mEstNoiseStd    = zeros(length(vNumRealizations), 2);
vN              = zeros(numRows * numCols * vNumRealizations(end), 1);

for ii = 1:length(vNumRealizations)
    numRealizations             = vNumRealizations(ii);
    mMeanImage(:)               = mean(tI(:, :, 1:numRealizations), 3);
    tN(:, :, 1:numRealizations) = tI(:, :, 1:numRealizations) - mMeanImage;
    numNoiseSamples             = numRows * numCols * numRealizations;
    vN(1:numNoiseSamples)       = reshape(tN(:, :, 1:numRealizations), numNoiseSamples, 1); %<! In MATLAB R2018b and above can use 'all' to skip this
    mEstNoiseStd(ii, 1)         = sqrt(mean(vN(1:numNoiseSamples) .^ 2)); % std(vN);
    mEstNoiseStd(ii, 2)         = mean(reshape(std(tI(:, :, 1:numRealizations), 0, 3), numRows * numCols, 1));
    mEstNoiseStd(ii, 3)         = sqrt(mean(reshape(var(tI(:, :, 1:numRealizations), 0, 3), numRows * numCols, 1)));
end

%% Analysis

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge);
hAxes       = axes();
hLineObj    = line(vNumRealizations, [mEstNoiseStd, noiseStd * ones(length(vNumRealizations), 1)]);
set(hLineObj, 'LineWidth', lineWidthNormal);
set(hLineObj(end), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', {['Estimation of the Noise STD as a Function of Number of Realizations']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Number of Realizations']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['STD']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['Estimated STD - Method 1'], ['Estimated STD - Method 2'], ['Estimated STD - Method 3'], ['Ground Truth']});

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Restore Defaults

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

author: RoyiAvital

repo: StackExchangeCodes


text:
function generateOpenfieldPB(fieldSizeX, fieldSizeY, beamletSize);

xP = -(fieldSizeX-beamletSize)/2 : beamletSize : (fieldSizeX-beamletSize)/2;
yP = -(fieldSizeY-beamletSize)/2 : beamletSize : (fieldSizeY-beamletSize)/2;
[xPosV, yPosV] = meshgrid(xP, yP);
xPosV = xPosV(:);
yPosV = yPosV(:);
w_field = ones(size(xPosV));
beamlet_delta_x = beamletSize*ones(size(xPosV));
beamlet_delta_y = beamlet_delta_x;

filename = ['PB_', num2str(fieldSizeX), 'x', num2str(fieldSizeY), '_PBsize', num2str(beamletSize), 'cm.mat']
disp ('save PB info to above filename.')
save (filename, 'xPosV', 'yPosV', 'beamlet_delta_x', 'beamlet_delta_y', 'w_field');

% Display generated PB information.
figure;hAxis2 = axes;hold on;
%axis(hAxis2, 'manual');
           %     w_colors = floor((w_field ./ max(w_field))*255)+1;
                %set(gcf, 'doublebuffer', 'on');
                for i=1:length(xPosV)
                    patch([xPosV(i) - beamlet_delta_x(i)/2 xPosV(i) - beamlet_delta_x(i)/2 xPosV(i) + beamlet_delta_x(i)/2 xPosV(i) + beamlet_delta_x(i)/2 xPosV(i) - beamlet_delta_x(i)/2], [yPosV(i) - beamlet_delta_y(i)/2 yPosV(i) + beamlet_delta_y(i)/2 yPosV(i) + beamlet_delta_y(i)/2 yPosV(i) - beamlet_delta_y(i)/2 yPosV(i) - beamlet_delta_y(i)/2], w_field(i));
                end
            %    axis([hAxis1 hAxis2], 'ij');
            %    kids = get(hAxis2, 'children');
            %     set(kids, 'edgecolor', 'none');
            %    cMap = colormap('jet');
            %   set(hAxis2, 'color', cMap(1,:));

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% resolve more uncorrelated sources than the number of array elements using
% co-prime arrays
clear(); close all;

wavelength = 1; % normalized
d = wavelength / 2;
% this co-prime array only have 9 elements
design_cp = design_array_1d('coprime', [3 4], d);
% but we have 10 sources here
doas = linspace(-pi/3, pi/3, 10);
power_source = 1;
power_noise = 1;
snapshot_count = 100;
source_count = length(doas);

% stochastic (unconditional) model
[~, R] = snapshot_gen_sto(design_cp, doas, wavelength, snapshot_count, power_noise, power_source);

% SS-MUSIC and DA-MUSIC
[Rss, dss] = virtual_ula_cov_1d(design_cp, R, 'SS');
[Rda, dda] = virtual_ula_cov_1d(design_cp, R, 'DA');
sp_ss = music_1d(Rss, source_count, dss, wavelength, 1440);
sp_da = music_1d(Rda, source_count, dda, wavelength, 1440);
sp_ss.true_positions = doas;
sp_da.true_positions = doas;
figure;
subplot(2,1,1);
plot_sp(sp_ss, 'Title', ['SS-MUSIC using ' design_cp.name], 'ReuseFigure', true);
subplot(2,1,2);
plot_sp(sp_da, 'Title', ['DA-MUSIC using ' design_cp.name], 'ReuseFigure', true);

% apply MVDR directly
sp_mvdr = mvdr_1d(R, source_count, design_cp, wavelength, 1440, 'RefineEstimates', true);
sp_mvdr.true_positions = doas;
plot_sp(sp_mvdr, 'Title', ['Direct MVDR using ' design_cp.name]);

% use sparse recovery
sp_sparse = sparse_bpdn_1d(R, source_count, design_cp, wavelength, 360, 9, ...
    'Formulation', 'ConstrainedL2');
sp_sparse.true_positions = doas;
plot_sp(sp_sparse, 'Title', 'Sparse Recovery based DOA Estimation');

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
function binv = beta_inv (x, a, b)
%  BETA_INV Inverse of the beta cumulative distribution function (cdf).
%     X = BETA_INV(P,A,B) returns the inverse of the beta cdf with
%     parameters A and B at the values in P.
%
% Copyright (c) 1995-1997, 2005-2007 Kurt Hornik

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if (nargin ~= 3)
  error('must provide 3 parameters')
end

if (~isscalar (a) || ~isscalar(b))
  if ((size(a,1) ~= size(b,1)) || (size(b,1) ~= size(x,1)))
    error ('betainv: x, a and b must be of common size or scalars');
  end
end

[n,nin] = size(x);
binv = zeros(n,nin);

k = find ((x<0) || (x>1) || (a<0) || (b<0) || isnan(x));
if (any(k))
  binv(k) = NaN;
end

k = find ((x==1) && (a>0) && (b>0));
if (any (k))
  binv(k) = 1;
end

k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0));
if (any (k))
  if (~isscalar(a) || ~isscalar(b))
    a = a(k);
    b = b(k);
    y = a./(a+b);
  else
    y = a/(a + b)*ones(size(k));
  end
  x = x(k);
  l = find (y < eps);
  if (any (l))
    y(l) = sqrt (eps) * ones (length (l), 1);
  end
  l = find (y > 1 - eps);
  if (any (l))
    y(l) = 1 - sqrt (eps) * ones (length (l), 1);
  end

  y_old = y;
  for i = 1 : 10000
    h     = (beta_cdf (y_old, a, b) - x) ./ beta_pdf (y_old, a, b);
    y_new = y_old - h;
    ind   = find (y_new <= eps);
    if (any (ind))
      y_new (ind) = y_old (ind) / 10;
    end
    ind = find (y_new >= 1 - eps);
    if (any (ind))
      y_new (ind) = 1 - (1 - y_old (ind)) / 10;
    end
    h = y_old - y_new;
    if (max (abs (h)) < sqrt (eps))
      break;
    end
    y_old = y_new;
  end

  binv(k) = y_new;
end

end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
classdef ipfHKLKey < ipfColorKey
  % defines an orientation mapping based on a certain inverse pole figure
  %   Detailed explanation goes here

  methods
    function oM = ipfHKLKey(varargin)
      oM = oM@ipfColorKey(varargin{:});
      oM.CS1 = oM.CS1.Laue;

      oM.dirMap = HKLDirectionKey(oM.CS1);
    end

  end
end

author: mtex-toolbox

repo: mtex