text:
function A = nnls_spatial(Y, A, C, active_pixel, maxN)
%% run HALS by fixating all spatial components
% input:
%   Y:  d*T, fluorescence data
%   A:  d*K, spatial components
%   C:  K*T, temporal components
%   active_pixel: d*T binary matrix, indicating the nonzero elements of A
%   maxN: scalar, maximum number of neurons overlapping at one pixel
% output:
%   A: d*K, updated spatial components

% Author: Pengcheng Zhou, Carnegie Mellon University, adapted from Johannes

%% options for HALS
d = size(Y, 1);
K = size(C, 1);
if nargin<5;    maxN = 5;    end;    %maximum iteration number
if nargin<4;    active_pixel=true(d, K);
elseif isempty(active_pixel)
    active_pixel = true(d, K);
else
    active_pixel = logical(active_pixel);
end;     %determine nonzero pixels

%% initialization
Ymean = mean(Y,2);
Y = bsxfun(@minus, Y, Ymean);
C = bsxfun(@minus, C, mean(C,2));
CC = C*C';
YC = C*Y';
ind_fit = find(sum(active_pixel,2)>1e-9);
A = zeros(size(A));

%% updating
for m=1:length(ind_fit)
    ind = active_pixel(ind_fit(m), :);
    A(ind_fit(m), ind) = nnls(CC(ind, ind), YC(ind, ind_fit(m)), [], 1e-4, maxN);
end

function s = nnls(A, b, s, tol, maxIter)
%% fast algorithm for solving nonnegativity constrained least squared
% problem minize norm(y-K*s, 2), s.t. s>=0.

%% inputs:
%   A: n x p matrix, K'*K
%   b: n x 1 vector, K'*y
%   s: p x 1 vector, warm started s
%   tol: scalar, smallest nonzero values
%   maxIter: scalar, maximum nonzero values

%% outputs:
%   s: p x 1 vector, solution

%% Authors: Pengcheng Zhou, Carnegie Mellon University, 2016
% ported from the Python implementation from Johannes Friedrich

%% References
% Friedrich J et.al., NIPS 2016, Fast Active Set Method for Online Spike Inference from Calcium Imaging
% Bro R & Jong S, Journal of Chemometrics 1997, A FAST NON-NEGATIVITY-CONSTRAINED LEAST SQUARES ALGORITHM

%% input arguments
p = size(A,2);      % dimension of s
if ~exist('s', 'var') || isempty(s)
    s = zeros(p, 1);
end
if ~exist('tol', 'var') || isempty(tol)
    tol = 1e-9;
end
if ~exist('maxIter', 'var') || isempty(maxIter)
    maxIter = p;
end
if sum(s>0)>maxIter
    s = zeros(p,1);
end
for miter=1:maxIter
    l = b - A*s;            % negative gradient
    Pset = (s>0);       % passive set

    if max(l) < tol         % no more passive set
        break;
    end

    [~, temp] = max(l);         % choose the one with the largest gradient
    Pset(temp) = true;         % add it to the passive set
    if sum(Pset)>maxIter
        break;
    end
    % correct nonnegativity violations
    while any(Pset)
        % run unconstrained least squares for variables in passive sets
        try
            mu = A(Pset, Pset) \ b(Pset);
        catch
            mu = (A(Pset, Pset) + tol*eye(sum(Pset))) \ b(Pset);
        end

        if all(mu>tol)
            break;
        end

        temp = s(Pset) ./ (s(Pset)-mu);
        temp(mu>tol) = [];
        a = min(temp);
        s(Pset) = s(Pset) + a*(mu-s(Pset));
        Pset(s<tol) = false;
    end

    s(Pset) = mu;
end

author: zhoupc


text:
function w = compute(w);

%COMPUTE calculates the fast fourier transform of a waveform object
% WAVEFORM = COMPUTE(WAVEFORM) calculate frequency spectrum of a waveform.
% The results of the fast fourier transform is added as new fields
% in the waveform:
%   FFT_FREQ is a vector of frequencies with N samples
%   FFT_AMP is a vector of spetral amplitudes
%   FFT_PHASE is a vector of phases
%   FFT_DOM is the scalar frequency of the maximum amplitude
%           peak (or dominant frequency. This could change)

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

% CHECK ARGUMENTS
if ~strcmpi(class(w),'waveform')
    error('First input must be a waveform object');
end

% STEP THROUGH WAVFORMS ADDING NEW FIELDS
[N,M] = size(w);
for i = 1:N*M
    Fn = get(w(i),'NYQ');
    x = get(w(i),'DATA');
    NFFT=2.^(ceil(log(length(x))/log(2)));  % Next highest power of 2
    FFTX=fft(x,NFFT);                       % Take fft, padding with zeros.
    NumUniquePts = ceil((NFFT+1)/2);
    FFTX=FFTX(1:NumUniquePts);              % throw out neg frequencies
    MX=abs(FFTX);                           % Take magnitude of X
    MX=MX*2;                                % Multiply by 2
    MX=MX/length(x);
    PX=phase(FFTX);                           % Take magnitude of X
    f=(0:NumUniquePts-1)*2/NFFT;
    f=f*Fn;
    w(i) = addfield(w(i),'FFT_FREQ',f');
    w(i) = addfield(w(i),'FFT_AMP',MX);
    w(i) = addfield(w(i),'FFT_PHASE',PX);
    a = find(MX == max(MX));
    w(i) = addfield(w(i),'FFT_DOM',f(a));
end;

author: geoscience-community-codes

repo: GISMO


text:
%              Apart = extractA(At,Ajc,blk0,blk1,blkstart[,blkstart2])
% EXTRACTA  Fast alternative to
%  Apart = At(blkstart(1):blkstart(2)-1,:).
%  Instead of blkstart(2), it takes "blkstart2" (if supplied) or
%  size(At,1)+1 (if neither blkstart(2) nor blkstart2) are available.
%
%  Extract submatrix of
%  A with subscripts "blkstart(1):blkstart(2)-1" and has its nonzeros per column
%  bounded bij Ajc1:Ajc2. If Ajc1 = [] (Ajc2=[]) then start at column start
%  (end at column end) of A.
%
% ******************** INTERNAL FUNCTION OF SEDUMI ********************
%
%  See also partitA

function Apart = extractA(At,Ajc,blk0,blk1,blkstart,blkstart2) %#ok
%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

disp('The SeDuMi binaries are not installed.')
disp('In Matlab, launch "install_sedumi" in the folder you put the SeDuMi files.')
disp('For more information see the file Install.txt.')
error(' ')

repo: REST


text:
% Test file for ADCHEBFUN feval and jump functions

function pass = test_fevalJump

% List of airy functions to test, with different number of inputs

fevalFun = @(u) feval(u, .3);
jumpFun = @(u) jump(u, .5, -1);
% List of trigonometric functions to test.
funcList = {fevalFun, jumpFun};

% Tolerance for Taylor testing
tolOrder = 1e-2;
tolDiff = 1e-12;
% Initialise vector with pass information
pass = zeros(2, numel(funcList));

% Do the tests.
for k = 1:numel(funcList)
    % Call the valueTesting method, which also returns linearity information
    [err, lin] = adchebfun.valueTesting(funcList{k});

    % First, check that the computed function values match what we expect
    pass(1, k) = ( err == 0 );

    % Call the taylorTesting method
    [order1, order2, nDiff2] = adchebfun.taylorTesting(funcList{k});

    % We expect all elements of ORDER1 to be close to 1. Since the methods being
    % tested in this case are all linear, ORDER2 will be noise. However, since
    % the methods are indeed linear, we should expect nDiff2 to have values all
    % close to machine epsilon, which we can use to check for the correctness of
    % the derivative computed.
    %
    % However, the jump method will run into issues of dividing 0 by 0. This
    % will result in NaNs, which can only occur if it gets the derivative
    % correct. So check instead that nDiff2 has all zero entries.
    if ( k ~= 2)
        pass(2,k) = ( (max(abs(order1 - 1)) < tolOrder) && ...
        (max(abs(nDiff2)) < tolDiff) );
    else
        pass(2,2) = ~any(nDiff2);
    end

    % Check that we received the correct linearity information
    pass(3, k) = ( lin == 1 );
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% MSH_EVAL_BOUNDARY_SIDE: evaluate the parameterization in one boundary side of the domain.
%
%     [msh_side, msh_side_from_interior] = msh_eval_boundary_side (msh, iside, [element_list]);
%
% INPUTS:
%
%    msh:   mesh object (see msh_cartesian)
%    iside: number of the boundary side to compute, from 1 to 2*msh.ndim (see the file geo_specs for documentation about face numbering)
%    element_list: elements on which to compute the parametrization
%
% OUTPUT:
%
%     msh_side: structure that contains the following fields
%
%     FIELD_NAME    (SIZE)                  DESCRIPTION
%     side_number   (scalar)                  number of the side
%     nel           (scalar)                  number of elements of the boundary side
%     nel_dir       (1 x ndim vector)         number of elements in each parametric direction
%     nqn           (scalar)                  number of quadrature nodes per element
%     nqn_dir       (1 x ndim vector)         number of quadrature nodes per element in each parametric direction
%     quad_nodes    (ndim x nqn x nel vector) coordinates of the quadrature nodes in parametric domain
%     quad_weights  (nqn x nel vector)        weights associated to the quadrature nodes
%     geo_map       (rdim x nqn x nel vector) physical coordinates of the quadrature nodes
%     geo_map_jac   (rdim x ndim x nqn x nel) Jacobian matrix of the map evaluated at the quadrature nodes
%     jacdet        (nqn x nel)               element of length, area, volume (if rdim = ndim, determinant of the Jacobian)
%
%     msh_side_from_interior: mesh structure that contains quadrature
%       points on the boundary, but which computes information from the
%       volumetric parametrization, like the derivative in the normal direction
%
% Copyright (C) 2009, 2010 Carlo de Falco
% Copyright (C) 2011, 2014, 2015, 2017 Rafael Vazquez
% Copyright (C) 2017 Luca Coradello
%
%    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/>.

function [msh_side, msh_side_from_interior] = msh_eval_boundary_side (msh, iside, element_list)

if (nargin == 2)
  msh_side = struct (msh_precompute (msh.boundary(iside)));
  msh_side.side_number = iside;

  if (msh.ndim > 1)

    %    ind  = [2 3; 2 3; 1 3; 1 3; 1 2; 1 2] in 3D, %ind  = [2 2 1 1] in 2D;
    %    ind2 = [1 1 2 2 3 3] in 3D,                  %ind2 = [1 1 2 2] in 2D
    ind2 = ceil (iside/2);
    ind = setdiff (1:msh.ndim, ind2);

    % Compute the normal vector. This requires the derivative also in the
    % normal direction (the boundary manifold is not enough).
    qn(ind) = msh_side.qn;
    if (mod (iside, 2) == 0)
      qn{ind2} = msh.breaks{ind2}(end);
    else
      qn{ind2} = msh.breaks{ind2}(1);
    end

    % Auxiliary vector sizes, to use with reshape and permute
    reorder = @(x) x(:)';
    psize = reorder ([msh_side.nqn_dir; msh_side.nel_dir]);
    vorder = [1:2:msh_side.ndim*2, 2:2:msh_side.ndim*2]; % [1 3 5 2 4 6], for ndim = 3

    jac = feval (msh.map_der, cellfun (reorder, qn, 'UniformOutput', false));
    jac = reshape (jac, [msh.rdim, msh.ndim, psize]);
    jac = permute (jac, [1 2 vorder+2]);

    geo_map_jac = reshape(jac, msh.rdim, msh.ndim, msh_side.nqn, msh_side.nel);
    JinvT = geopdes_invT__ (geo_map_jac);
    JinvT = reshape (JinvT, [msh.rdim, msh.ndim, msh_side.nqn, msh_side.nel]);

    normal = zeros (msh.ndim, msh_side.nqn, 1, msh_side.nel);
    normal(ind2,:) = (-1)^iside;
    normal = geopdes_prod__ (JinvT, normal);
    normal = reshape (normal, [msh.rdim, msh_side.nqn, msh_side.nel]);
    % Now normalize
    norms = reshape (geopdes_norm__ (normal), [1, msh_side.nqn, msh_side.nel]);
    msh_side.normal = bsxfun (@rdivide, normal, norms);

    % Compute the characteristic length in the normal direction
    % This has to be computed using information from the interior
    qn = msh.qn; qw = msh.qw; nel_dir = msh.nel_dir;
    if (mod (iside, 2) == 0)
      qn{ind2} = qn{ind2}(:,end);
      if (~isempty(msh.qw))
        qw = qw{ind2}(:,end);
      else
        qw = [];
      end
    else
      qn{ind2} = qn{ind2}(:,1);
      if (~isempty(msh.qw))
        qw = qw{ind2}(:,1);
      else
        qw = [];
      end
    end
    nel_dir(ind2) = 1;

    reorder = @(x) x(:)';
    psize = reorder ([msh.nqn_dir; nel_dir]);
    vorder = [1:2:msh.ndim*2, 2:2:msh.ndim*2]; % [1 3 5 2 4 6], for ndim = 3

    jac = feval (msh.map_der, cellfun (reorder, qn, 'UniformOutput', false));
    jac = reshape (jac, [msh.rdim, msh.ndim, psize]);
    jac = permute (jac, [1 2 vorder+2]);
    jac = reshape (jac, [msh.rdim, msh.ndim, msh.nqn_dir, msh_side.nel]);

    jac = reshape (sqrt (sum (jac(:,ind2,:,:,:,:).^2, 1)), [msh.nqn_dir, msh_side.nel]); % Module

    qsize = ones (1, msh.ndim);
    qsize(ind2) = numel (qw);
    jac_times_qw = bsxfun(@times, jac, reshape (qw, qsize));
    msh_side.charlen = reshape (sum (jac_times_qw, ind2), msh_side.nqn, msh_side.nel);
  end

elseif (nargin == 3)

  if (isempty (element_list))
    msh_side.quad_weights = [];
    msh_side.geo_map = [];
    msh_side.geo_map_jac = [];
    msh_side.geo_map_der2 = [];
    msh_side.jacdet = [];
    msh_side.element_size = [];
    msh_side.normal = [];
    return
  end

  msh_side = msh_evaluate_element_list (msh.boundary(iside), element_list);
  msh_side.side_number = iside;

  if (msh.ndim > 1)

    %    ind  = [2 3; 2 3; 1 3; 1 3; 1 2; 1 2] in 3D, %ind  = [2 2 1 1] in 2D;
    %    ind2 = [1 1 2 2 3 3] in 3D,                  %ind2 = [1 1 2 2] in 2D
    ind2 = ceil (iside/2);
    ind = setdiff (1:msh.ndim, ind2);

    % Compute the normal vector. This requires the derivative also in the
    % normal direction (the boundary manifold is not enough).
    element_list = element_list(:)';
    indices = cell (msh_side.ndim, 1);
    [indices{:}] = ind2sub (msh_side.nel_dir, element_list);
    indices = cell2mat (indices);

    qn_elems = arrayfun(@(ii) {msh.boundary(iside).qn{ii}(:,indices(ii,:))}, 1:msh_side.ndim);
    qqn = cell (1,msh_side.nel);
    for iel = 1:numel(element_list)
      for idim = 1:numel(ind)
        qqn{iel}{ind(idim)} = qn_elems{idim}(:,iel)';
      end
      if (mod (iside, 2) == 0)
        qqn{iel}{ind2} = msh.breaks{ind2}(end);
      else
        qqn{iel}{ind2} = msh.breaks{ind2}(1);
      end
    end

    jac = cellfun (@(x) feval (msh.map_der, x), qqn, 'UniformOutput', false);
    geo_map_jac  = zeros (msh.rdim, msh.ndim, msh_side.nqn, msh_side.nel);
    for iel = 1:numel(element_list)
      geo_map_jac(:,:,:,iel) = jac{iel};
    end
    JinvT = geopdes_invT__ (geo_map_jac);
    JinvT = reshape (JinvT, [msh.rdim, msh.ndim, msh_side.nqn, msh_side.nel]);
    normal = zeros (msh.ndim, msh_side.nqn, 1, msh_side.nel);
    normal(ind2,:) = (-1)^iside;
    normal = geopdes_prod__ (JinvT, normal);
    normal = reshape (normal, [msh.rdim, msh_side.nqn, msh_side.nel]);
    % Now normalize
    norms = reshape (geopdes_norm__ (normal), [1, msh_side.nqn, msh_side.nel]);
    msh_side.normal = bsxfun (@rdivide, normal, norms);

  end

end

if (nargout == 2)
  brk_bnd = msh.breaks; qn_bnd = msh.qn; qw_bnd = msh.qw;
  ind2 = ceil (iside/2);
  if (mod (iside, 2) == 1)
    brk_bnd{ind2} = brk_bnd{ind2}(1:2);
    qn_bnd{ind2} = brk_bnd{ind2}(1);
    qw_bnd{ind2} = 1;
  else
    brk_bnd{ind2} = brk_bnd{ind2}(end-1:end);
    qn_bnd{ind2} = brk_bnd{ind2}(end);
    qw_bnd{ind2} = 1;
  end

  geo.map = msh.map; geo.map_der = msh.map_der; geo.map_der2 = msh.map_der2;
  geo.rdim = msh.rdim;
  msh_side_from_interior = msh_cartesian (brk_bnd, qn_bnd, qw_bnd, geo, 'boundary', false);
  if (nargin == 2)
    msh_side_from_interior = msh_precompute (msh_side_from_interior);
  elseif (nargin == 3)
    msh_side_from_interior = msh_evaluate_element_list (msh_side_from_interior, element_list);
  end
end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [f,t,w]=enframe(x,win,inc,m)
%ENFRAME split signal up into (overlapping) frames: one per row. [F,T]=(X,WIN,INC)
%
% Usage:  (1) f=enframe(x,n)     % split into frames of length n
%         (2) f=enframe(x,hamming(n,'periodic'),n/4)     % use a 75% overlapped Hamming window of length n
%         (3) frequency domain frame-based processing:
%
%             S=...;                              % input signal
%             OV=2;                               % overlap factor of 2 (4 is also often used)
%             INC=20;                             % set frame increment in samples
%             NW=INC*OV;                          % DFT window length
%             W=sqrt(hamming(NW,'periodic'));     % omit sqrt if OV=4
%             W=W/sqrt(sum(W(1:INC:NW).^2));      % normalize window
%             F=rfft(enframe(S,W,INC),NW,2);      % do STFT: one row per time frame, +ve frequencies only
%             ... process frames ...
%             X=overlapadd(irfft(F,NW,2),W,INC);  % reconstitute the time waveform (omit "X=" to plot waveform)
%
%  Inputs:   x    input signal
%          win    window or window length in samples
%          inc    frame increment in samples
%            m    mode input:
%                  'z'  zero pad to fill up final frame
%                  'r'  reflect last few samples for final frame
%                  'A'  calculate the t output as the centre of mass
%                  'E'  calculate the t output as the centre of energy
%
% Outputs:   f    enframed data - one frame per row
%            t    fractional time in samples at the centre of each frame
%                 with the first sample being 1.
%            w    window function used
%
% By default, the number of frames will be rounded down to the nearest
% integer and the last few samples of x() will be ignored unless its length
% is lw more than a multiple of inc. If the 'z' or 'r' options are given,
% the number of frame will instead be rounded up and no samples will be ignored.
%
% Example of frame-based processing:
%          INC=20       						% set frame increment in samples
%          NW=INC*2     						% oversample by a factor of 2 (4 is also often used)
%          S=cos((0:NW*7)*6*pi/NW);				% example input signal
%          W=sqrt(hamming(NW),'periodic'));  	% sqrt hamming window of period NW
%          F=enframe(S,W,INC);               	% split into frames
%          ... process frames ...
%          X=overlapadd(F,W,INC);               % reconstitute the time waveform (omit "X=" to plot waveform)

% Bugs/Suggestions:
%  (1) Possible additional mode options:
%        'u'  modify window for first and last few frames to ensure WOLA
%        'a'  normalize window to give a mean of unity after overlaps
%        'e'  normalize window to give an energy of unity after overlaps
%        'wm' use Hamming window
%        'wn' use Hanning window
%        'x'  include all frames that include any of the x samples

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

nx=length(x(:));
if nargin<2 || isempty(win)
    win=nx;
end
if nargin<4 || isempty(m)
    m='';
end
nwin=length(win);
if nwin == 1
    lw = win;
    w = ones(1,lw);
else
    lw = nwin;
    w = win(:).';
end
if (nargin < 3) || isempty(inc)
    inc = lw;
end
nli=nx-lw+inc;
nf = max(fix(nli/inc),0);   % number of full frames
na=nli-inc*nf+(nf==0)*(lw-inc);       % number of samples left over
fx=nargin>3 && (any(m=='z') || any(m=='r')) && na>0; % need an extra row
f=zeros(nf+fx,lw);
indf= inc*(0:(nf-1)).';
inds = (1:lw);
if fx
    f(1:nf,:) = x(indf(:,ones(1,lw))+inds(ones(nf,1),:));
    if any(m=='r')
        ix=1+mod(nf*inc:nf*inc+lw-1,2*nx);
        f(nf+1,:)=x(ix+(ix>nx).*(2*nx+1-2*ix));
    else
        f(nf+1,1:nx-nf*inc)=x(1+nf*inc:nx);
    end
    nf=size(f,1);
else
    f(:) = x(indf(:,ones(1,lw))+inds(ones(nf,1),:));
end
if (nwin > 1)   % if we have a non-unity window
    f = f .* w(ones(nf,1),:);
end
if nargout>1
    if any(m=='E')
        t0=sum((1:lw).*w.^2)/sum(w.^2);
    elseif any(m=='A')
        t0=sum((1:lw).*w)/sum(w);
    else
        t0=(1+lw)/2;
    end
    t=t0+inc*(0:(nf-1)).';
end

repo: ekfmonoslam


text:
classdef Display_Plot < Display_Abstract

    methods (Access = protected)

        function setChartType(obj)
            obj.handle = plot(0,0);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
classdef(Abstract) AbstractEventAction < matlab.mixin.SetGet & matlab.mixin.Heterogeneous
    %AbstractEventAction Summary of this class goes here
    %   Detailed explanation goes here

    properties(Abstract=false)
        id(1,1) double
        event LaunchVehicleEvent
    end

    methods
        newStateLogEntry = executeAction(obj, stateLogEntry)

        initAction(obj, initialStateLogEntry)

        name = getName(obj)

        tf = usesStage(obj, stage)

        tf = usesEngine(obj, engine)

        tf = usesTank(obj, tank)

        tf = usesEngineToTankConn(obj, engineToTank)

        tf = usesStopwatch(obj, stopwatch)

        tf = usesExtremum(obj, extremum)

        tf = usesTankToTankConn(obj, tankToTank)

        function tf = usesCalculusCalc(obj, calculusCalc)
            tf = false;
        end

        function tf = usesEvent(obj, event)
            tf = false;
        end

        function tf = usesPwrSink(obj, powerSink)
            tf = false;
        end

        function tf = usesPwrSrc(obj, powerSrc)
            tf = false;
        end

        function tf = usesPwrStorage(obj, powerStorage)
            tf = false;
        end

        function tf = usesSensor(obj, sensor)
            tf = false;
        end

        [tf,vars] = hasActiveOptimVar(obj)
    end

    methods(Sealed)
        function tf = eq(A,B)
            tf = [A.id] == [B.id];
        end
    end

    methods(Static)
        addActionTf = openEditActionUI(action, lv);
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function mean = lorentz_mean ( )

%*****************************************************************************80
%
%% LORENTZ_MEAN returns the mean of the Lorentz PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 1999
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real MEAN, the mean of the PDF.
%
  mean = 0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Example
% used to test the timing characteristic of the connection between IIWA and PC

% First run the following script in Matlab
% Then start the client on the KUKA iiwa controller

% Mohammad SAFEEA, 24th of April 2018

close all,clear all;clc;

ip='172.31.1.147'; % The IP of the controller
% start a connection with the server
global t_Kuka;
t_Kuka=net_establishConnection( ip );

if ~exist('t_Kuka','var') || isempty(t_Kuka) || strcmp(t_Kuka.Status,'closed')
  disp('Connection could not be establised, script aborted');
  return;
else
    disp('Connection established')
end
pause(1);
disp('This will take up a minute or so to finish, please wait')

[time_stamps,time_delay]=net_updateDelay(t_Kuka);

net_turnOffServer( t_Kuka );

repo: KST-Kuka-Sunrise-Toolbox


text:
%% DEMO_visualization_von_mises_plasticity_01
% Below is a demonstration for:
%
% * Visualization of a Von Mises yield surface in 3D

%%

clear; close all; clc

%% Plot settings.

fontSizeAxis=25;
fontSizeText=35;
colorSet=viridis(3);
lineWidth1=1;
lineWidth2=3;
lineWidthAxis=3;
f=2;

%% Control parameters

sy=1; %Yield stress
tau_octahedral=sy*sqrt(2/3); %Octahedral shear stress = cylinder radius

%%
% Create ellipsoids
t=linspace(0,2*pi,1000); %Angular coordinates
x=tau_octahedral*cos(t);
y=tau_octahedral*sqrt(3).*sin(t);
z=zeros(size(x));
v=[x(:) y(:) z(:)]; %Ellipse coordinates

%Make cylinder "slices" ellipses by rotating above ellipse
R12=euler2DCM([0 0 (45/180)*pi]);
v12=v*R12;

R13_1=euler2DCM([0 (90/180)*pi 0]);
R13_2=euler2DCM([(135/180)*pi 0 0]);
v13=v*R13_1*R13_2;

R23_1=euler2DCM([(90/180)*pi 0 0]);
R23_2=euler2DCM([0 (-45/180)*pi 0 ]);
v23=v*R23_1*R23_2;

%%
% Creating cylinder data

% Creating input structure
inputStruct.cylRadius=tau_octahedral;
inputStruct.numRadial=250;
inputStruct.cylHeight=2*sqrt(2)*tau_octahedral;
inputStruct.numHeight=10;
inputStruct.meshType='quad';
[F,V,C]=patchcylinder(inputStruct);
Q=euler2DCM([0 pi/2 0]);
V=V*Q;

R1=euler2DCM([0 asin(1/sqrt(3)) 0]);
R2=euler2DCM([0 0 -(45/180)*pi]);
V=V*R1*R2;

[x1,y1]=meshgrid(-f*sy:1:f*sy);
z1=zeros(size(x1));

[x2,z2]=meshgrid(-f*sy:1:f*sy);
y2=zeros(size(x2));

[y3,z3]=meshgrid(-f*sy:1:f*sy);
x3=zeros(size(y3));

%%
%Visualize
hf=cFigure; hold on;

h4=gpatch(F,V,'rw','none',0.5);
surf(x1,y1,z1,'EdgeColor','k','faceColor',1*ones(1,3),'EdgeAlpha',0.5,'FaceAlpha',0.1,'LineWidth',lineWidth1);
surf(x2,y2,z2,'EdgeColor','k','faceColor',1*ones(1,3),'EdgeAlpha',0.5,'FaceAlpha',0.1,'LineWidth',lineWidth1);
surf(x3,y3,z3,'EdgeColor','k','faceColor',1*ones(1,3),'EdgeAlpha',0.5,'FaceAlpha',0.1,'LineWidth',lineWidth1);

h1=quiverVec([0 0 0],[1 0 0],sy,colorSet(1,:));
h2=quiverVec([0 0 0],[0 1 0],sy,colorSet(2,:));
h3=quiverVec([0 0 0],[0 0 1],sy,colorSet(3,:));

h5=plotV([-1 -1 -1; 1 1 1],'k--','LineWidth',lineWidth2); %Hydrostatic axis

h6=plotV(v12,'r-','LineWidth',lineWidth2);
h6.Color=(colorSet(1,:)+colorSet(2,:))/2;

h7=plotV(v13,'r-','LineWidth',lineWidth2);
h7.Color=(colorSet(1,:)+colorSet(3,:))/2;

h8=plotV(v23,'r-','LineWidth',lineWidth2);
h8.Color=(colorSet(2,:)+colorSet(3,:))/2;

hAxis=gca;
hAxis.XRuler.FirstCrossoverValue  = 0; % X crossover with Y axis
hAxis.YRuler.FirstCrossoverValue  = 0; % Y crossover with X axis
hAxis.ZRuler.FirstCrossoverValue  = 0; % Z crossover with X axis
hAxis.XRuler.SecondCrossoverValue = 0; % X crossover with Z axis
hAxis.YRuler.SecondCrossoverValue = 0; % Y crossover with Z axis
hAxis.ZRuler.SecondCrossoverValue = 0; % Z crossover with Y axis

text(0.5+f*sy,0,0,'$\sigma_1$','Interpreter','Latex','FontSize',fontSizeText);
text(0,0.5+f*sy,0,'$\sigma_2$','Interpreter','Latex','FontSize',fontSizeText);
text(0,0,0.5+f*sy,'$\sigma_3$','Interpreter','Latex','FontSize',fontSizeText);

xticks(-f*sy:0.5:f*sy);
yticks(-f*sy:0.5:f*sy);
zticks(-f*sy:0.5:f*sy);

legend([h1 h2 h3 h4 h5 h6 h7 h8],{'$\sigma_1$ axis','$\sigma_2$ axis','$\sigma_3$ axis',...
                           'Von Mises yield surface',...
                           'Hydrostatic line',...
                           'yield ellipse $\sigma_3=0$','yield ellipse $\sigma_2=0$','yield ellipse $\sigma_1=0$'...
                           },'Interpreter','Latex');

axis tight; axis equal; axis vis3d; view(3); %box on;
camlight headlight;
axis(f*sy*[-1 1 -1 1 -1 1])
set(gca,'FontSize',fontSizeAxis,'LineWidth',lineWidthAxis);
gdrawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [state,options] = psoplotscorediversity(options,state,flag)
% Plots a histogram containing the best and mean scores of particle swarm.

if strcmp(flag,'init')
    set(gca,'NextPlot','replacechildren',...
        'XLabel',xlabel('Scores'),...
        'YLabel',ylabel('Number of inidividuals'))
end

[n,bins] = hist(state.Score) ;
bar(bins,n,'Tag','scorehistogram','FaceColor',[0.1 0.1 0.5])


text:
function [AA, aa, pp, rankA, p] = rowReduce(A, a, mode, printLevel)
% Eliminates dependent rows from `A` & `a` where :math:`A x = a`
%
% USAGE:
%
%    [AA, aa, pp, rankA, p] = rowReduce(A, a)
%
% INPUT:
%    A:        from :math:`A x = a`
%
% OPTIONAL INPUT:
%    a:        from :math:`A x = a`
%
%    mode:     If mode=1, LUSOL operates on A itself.
%              If mode=2, LUSOL operates on A'.
%
%    printLevel
%
% OUTPUT:
%    AA:       row reduced `A`
%    aa:       row reduced `a` i.e. `aa = a(pp)`
%    pp:       1:rankA indices of independent rows
%    rankA:    rank of `A`
%    p:        row permutation which leaves first `1:rankA` rows independent and
%              last rows dependent
%
% .. Author: - Ronan Fleming, with linear algebra advice from Michael Saunders
%            Dept of Management Science and Engineering (MS&E) Stanford University

if ~exist('a','var') %create a if not provided
    a=sparse(size(A,1),1);
else
    if size(A,1)~=length(a)
        error('Dimensions of A and a are inconsistent');
    end
end

if ~exist('mode','var')
    mode=1;
end

if ~exist('printLevel','var')
    printLevel=1;
end

[mlt,nlt]=size(A);
archstr = computer('arch');
archstr = lower(archstr);
%archstr='';%bypass until issue with lusol tolerance sorted.
switch archstr
    case {'glnx86','glnxa64','maci64'}
        %Eliminate dependent rows
        [AA,aa,p,rankA] = lusolCondense(A,a,mode,printLevel-1);
        if ~(nnz(A)>0 && rankA==0) && printLevel>0
            fprintf('%s',['Eliminated ' int2str(mlt-rankA) ' dependent rows, using lusol.']);
        end
        pp=p(1:rankA);
        %case {'PCWIN','PCWIN64'}
    otherwise
        [AA, aa, pp, rankA, p] = qrRowReduce(A,a, printLevel);
end

if nnz(A)>0 && rankA==0
    %backup in case something has gone wrong with lusolCondense
    [AA, aa, pp, rankA, p] = qrRowReduce(A, a, printLevel);
end

fprintf('\n')
end

function [AA, aa, pp, rankA, p] = qrRowReduce(A, a,  printLevel)
    [mlt,nlt]=size(A);
    A=full(A);
    %Eliminate dependent rows
    [Q,R,P] = qr(A');
    [p,q,s] = find(P);
    s      = diag(R);
    tol    = 1e-8;
    rankA  = length(find(abs(s) > tol));
    % nnz(abs(diag(R))>1.e-15)
    AA = sparse(A(p(1:rankA),:));
    pp = p(1:rankA);
    aa = sparse(a(pp));
    if  printLevel>0
             fprintf('%s',['Eliminated ' int2str(mlt-rankA) ' dependent rows, using qr factorisation.']);
    end
end

author: opencobra

repo: cobratoolbox


text:
function [ M, N ] = getPsiMat( data, image )
%GETPSIMAT gets \Psi matrix (LBP features) for all components
%   Detailed explanation goes here
%
% Synopsis:
%  [ M ] = getPsiMat(data, image)
%
% INPUT:
%   data    ...     struct, containes filed lbp (cell array - holds
%                   neccessary parameters for lbppyr_features_sparse
%                   function), imSize (size of image), Images (matrix
%                   [imSize(1)*imSize(2) x nImages]), nImages (count of
%                   images), structure options (with field M - count of
%                   components)
%   image   ...     idx of image on which LBP features are computed
%
% OUTPUT:
%   M       ...     cell array [1 x data.options.M] - \Psi matrices with
%                   LBP features
%   N       ...     cell array [1 x data.options.M] - \Psi matrices with
%                   indices of nonzero elements of LBP features
%
% 05-08-10 Michal Uricar
% 11-07-11 Michal Uricar, corners dataset (checked-only)

    M = cell(1, data.options.M);
    N = cell(1, data.options.M);

    for j = 1 : data.options.M;
        [M{j} N{j}] = lbp2sparse(data.Images(:, image), data.imSize, data.lbp{j}.wins, data.lbp{j}.winSize, data.lbp{j}.hop, 0);
    end;

end

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
 function [cost, fit, reg] = pwls_cost(xs, A, W, yi, R, mask)
%function [cost, fit, reg] = pwls_cost(xs, A, W, yi, R, mask)
%|
%| compute PWLS cost for each column of xs
%|
%| in
%|	xs	[np niter]	iterates
%|	A	[nd np]		system matrix
%|	W	[nd nd]		data weighting matrix, usually diag_sp(wi)
%|	yi	[nd 1]		data
%|	R			penalty object (see Reg1.m or Robject.m)
%|				with penalty cost method: R.penal(R, x)
%|				or just *sparse* C matrix (quadratic only)
%|	mask	[nx ny ...]	optional mask, iff xs is [nx ny ... niter]
%|
%| out
%|	cost	[niter 1]	cost
%|	fit	[niter 1]	(y-A*x)'W(y-Ax)/2
%|	reg	[niter 1]	R(x)
%|
%| Copyright 2002-2-12, Jeff Fessler, University of Michigan

if nargin < 4, ir_usage, end

if ~isvar('R'), R = []; end

if isvar('mask') && ~isempty(mask)
	xs = reshapee(xs, numel(mask), []);	% [(*N) niter]
	xs = xs(mask(:), :);			% [np niter]
end

niter = size(xs,2);
reg = zeros(niter,1);

if isempty(R)
	warning 'empty R means no penalty'

elseif issparse(R) || isa(R, 'Fatrix')
	C = R; % trick!
	if size(C,2) == size(C,1)
		warning 'square C is quite unusual!?'
	end
	for ii=1:niter
		reg(ii) = sum(abs(C * xs(:,ii)).^2)/2;
	end

elseif isstruct(R) || isa(R, 'strum')
	for ii=1:niter
		reg(ii) = R.penal(R, xs(:,ii));
	end

else
	keyboard
	error 'bad R'
end

fit = zeros(niter,1);
for ii=1:niter
	resid = yi - A * xs(:,ii); % predicted measurements
	fit(ii) = resid' * (W * resid) / 2;
end

cost = fit + reg;
cost = reale(cost, 'warn'); % trick: x'*W*x is not always real for complex values

if ~nargout
	pr fit
	pr reg
	pr cost
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [qout] = Filter_Slerp(trajNoisy,hlpf)
%This function will use SLERP  to lowpass filter a noisy trajectory which
%is represented as quaternions.

% The new filter state should be moved toward the current input  by a step
% size that is proportional to the distance between the current input and
% the current filter state. To do this, both the dist and slerp functions
% are used. The dist function returns a measurement in radians of the
% difference in rotation applied by two quaternions. The range of the dist
% function is the half-open interval [0,pi).

%Set parameters - The interpolation parameter to slerp is in the
%closed-interval [0,1] so the output of dist must be re-normalized to this
%range. However, the full range of [0,1] for the interpolation parameter
%gives poor performance so it is limited to a smaller range hrange centered
%at hbias.

% All options have been kept the same as the matlab example...
% Wrtten by Jayishni Maharaj (UQ)
% 15 Feb 2018

hrange = 0.7;%0.4
hbias = 0.7;%0.4
%Limit low and high to the interval [0, 1].
low  = max(min(hbias - (hrange./2), 1), 0);
high = max(min(hbias + (hrange./2), 1), 0);
hrangeLimited = high - low;

%Initialize the filter and preallocate outputs.
y = trajNoisy(1); % initial filter state
qout = zeros(size(y), 'like', y); % preallocate filter output
qout(1) = y;

for ii=2:numel(trajNoisy)
    x = trajNoisy(ii); clear d
    d = dist(y, x);

    % Renormalize dist output to the range [low, high]
    %hlpf = (d./pi).*hrangeLimited + low;
    y = slerp(y,x,hlpf);
    qout(ii,:) = y;
end
end


text:
fprintf('\nHere we train Sparse and Denoising Autoencoders on the MNIST dataset.\n');

% LOAD THE MINIST DATASET
load('mnistSmall.mat','trainData','testData');

[nObs,nIn]= size(trainData);

nHid = 500;

% DEFINE SINGLE-LAYERED ARCHITECTURE
arch.size = [nIn nHid];
arch.actFun = {'sigmoid','sigmoid'}; % ALL SIGMOIDS
arch.lRate = [.1 .1];

% TRAINING OPTIONS
arch.opts = {'nEpoch',10, ...
             'sparsity',.01, ...
             'costFun','xent'};

% INITIALIZE AUTOENCODER
a1 = ae(arch);

% TRAIN & TEST
fprintf('\nTraining a Sparse Autoencoder...\n');
a1 = a1.train(trainData);
[~,xent1] = a1.test(testData);

nVis = 225;
figure; visWeights(a1.layers{1}.W',0,[-1 1],1,nVis);
title(sprintf('Encoding Features of the Sparse Autoencoder\n(sorted by L2 norm)'));
drawnow

% DEFINE DENOISING ARCHITECTURE
arch2 = arch;

% ADD 20% DENOISING RATIO
arch2.opts = {'denoise',.2, ...
              'nEpoch',10, ...
              'costFun','xent'};

% INITIALIZE NETWORK
a2 = ae(arch2);

% TRIAN & TEST
fprintf('\n\nNow training a Denoising Autoencoder...\n');
a2 = a2.train(trainData);
[~,xent2] = a2.test(testData);

figure; visWeights(a2.layers{1}.W',0,[-1 1],1,nVis);
title(sprintf('Encoding Features of the Denoising Autoencoder\n(sorted by L2 norm)'))
drawnow

% DISPLAY
fprintf('\nTesting Mean Squared Error for Autoencoder: %1.2f\n',xent1);
fprintf('\nTesting Mean Squared Error for Denoising Autoencoder: %1.2f\n',xent2);

author: dustinstansbury

repo: medal


text:
function b = ssisl ( a, lda, n, kpvt, b )

%*****************************************************************************80
%
%% SSISL solves a real symmetric system factored by SSIFA.
%
%  Discussion:
%
%    To compute inverse(A) * C where C is a matrix with P columns
%
%      call ssifa ( a, lda, n, kpvt, info )
%
%      if ( info == 0 ) then
%        do j = 1, p
%          call ssisl ( a, lda, n, kpvt, c(1,j) )
%        end do
%      end if
%
%    A division by zero may occur if the inverse is requested
%    and SSICO has set RCOND == 0.0D+00 or SSIFA has set INFO /= 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real A(LDA,N), the output from SSIFA.
%
%    Input, integer LDA, the leading dimension of the array A.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer KPVT(N), the pivot vector from SSIFA.
%
%    Input, real B(N), the right hand side.
%
%    Output, real B(N), the solution.
%

%
%  Loop backward applying the transformations and D inverse to B.
%
  k = n;

  while ( 0 < k )

    if ( 0 <= kpvt(k) )
%
%  1 x 1 pivot block.
%
      if ( k ~= 1 )

        kp = kpvt(k);
%
%  Interchange.
%
        if ( kp ~= k )
          temp = b(k);
          b(k) = b(kp);
          b(kp) = temp;
        end
%
%  Apply the transformation.
%
        b(1:k-1) = saxpy ( k-1, b(k), a(1:k-1,k), 1, b(1:k-1), 1 );

      end
%
%  Apply D inverse.
%
      b(k) = b(k) / a(k,k);
      k = k - 1;

    else
%
%  2 x 2 pivot block.
%
      if ( k ~= 2 )

        kp = abs ( kpvt(k) );
%
%  Interchange.
%
        if ( kp ~= k-1 )
          temp = b(k-1);
          b(k-1) = b(kp);
          b(kp) = temp;
        end
%
%  Apply the transformation.
%
        b(1:k-2) = saxpy ( k-2, b(k), a(1:k-2,k), 1, b(1:k-2), 1 );
        b(1:k-2) = saxpy ( k-2, b(k-1), a(1:k-2,k-1), 1, b(1:k-2), 1 );

      end
%
%  Apply D inverse.
%
      ak = a(k,k) / a(k-1,k);
      akm1 = a(k-1,k-1) / a(k-1,k);
      bk = b(k) / a(k-1,k);
      bkm1 = b(k-1) / a(k-1,k);
      denom = ak * akm1 - 1.0;
      b(k) = ( akm1 * bk - bkm1 ) / denom;
      b(k-1) = ( ak * bkm1 - bk ) / denom;
      k = k - 2;

    end

  end
%
%  Loop forward applying the transformations.
%
  k = 1;

  while ( k <= n )

    if ( 0 <= kpvt(k) )
%
%  1 x 1 pivot block.
%
      if ( k ~= 1 )
%
%  Apply the transformation.
%
        b(k) = b(k) + b(1:k-1)' * a(1:k-1,k);
        kp = kpvt(k);
%
%  Interchange.
%
        if ( kp ~= k )
          temp = b(k);
          b(k) = b(kp);
          b(kp) = temp;
        end

      end

      k = k + 1;

    else
%
%  2 x 2 pivot block.
%
      if ( k ~= 1 )
%
%  Apply the transformation.
%
        b(k) = b(k) + b(1:k-1) * a(1:k-1,k);
        b(k+1) = b(k+1) + b(1:k-1) * a(1:k-1,k+1);
        kp = abs ( kpvt(k) );
%
%  Interchange.
%
        if ( kp ~= k )
          temp = b(k);
          b(k) = b(kp);
          b(kp) = temp;
        end

      end

      k = k + 2;

    end

  end

  return
end


text:
function y = ngaussian(x)

% NGAUSSIAN Compute a Gaussian with mean 0 and variance 1.
% FORMAT
% DESC computes a the likelihood of a normalised Gaussian
% distribution, i.e. with mean 0 and variance 1.
% ARG x : input value(s) for which to compute the distribution.
% RETURN y : probability of the input values under the Gaussian.
%
% SEEALSO : cumGaussian, gaussOverDiffcumGaussian
%
% COPYRIGHT : Neil D. Lawrence, 2004

% NDLUTIL

x2 = x.*x;
y = exp(-.5*x2);
y = y/sqrt(2*pi);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef graphics

methods(Static)

    function result = plot_styles()
        % plot styles for multiple plots
        result = {...
        '--+','-.o',':*',...
        '--.','-.x',':s',...
        '--d','-.^',':v',...
        '-->','-.<',':p',...
        '--h'};
    end

end

end

author: indigits

repo: sparse-plex


text:
% Data File DTX2
% Free damped vibrations
% of a particle
  m    = 'm';    % mass of the particle
  Fx   = '-k*v - c*x'; % projections of forces on x
  x0   = '0.1';  % initial coordinate
  v0   = '10';   % initial velocity
  Tend = 20;     % upper bound of integration
  eps  = 1e-10;  % desirable accuracy
  np   = 3;      % number of parameters
  P{1} = 'm';    % mass of the particle
  P{2} = 'k';    % coefficient of resistance
  P{3} = 'c';    % spring stiffness

author: Sable

repo: mcbench-benchmarks


text:
function test_ft_denoise_prewhiten

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_denoise_prewhiten

% create some data
data = [];
data.label = {'chan01';'chan02';'chan03'};

% full rank data
mix = [1 0 0;0 1 0;0.5 0.5 0.1];
for k = 1:10
  data.trial{k} = mix*randn(3,1000);
  data.time{k}  = (0:999)./1000;
end

cfg = [];
cfg.covariance = 'yes';
tlck = ft_timelockanalysis(cfg, data);

cfg = [];
datawhite = ft_denoise_prewhiten(cfg, data, tlck);

cfg = [];
cfg.covariance = 'yes';
tlckwhite = ft_timelockanalysis(cfg, datawhite);

cfg = [];
tlckwhite2 = ft_denoise_prewhiten(cfg, tlck, tlck);

% rank deficient data
mix = [1 0 0;0 1 0;0.5 0.5 0];
for k = 1:10
  data.trial{k} = mix*randn(3,1000);
  data.time{k}  = (0:999)./1000;
end

cfg = [];
cfg.covariance = 'yes';
tlck = ft_timelockanalysis(cfg, data);

cfg = [];
datawhite = ft_denoise_prewhiten(cfg, data, tlck);

cfg = [];
cfg.covariance = 'yes';
tlckwhite = ft_timelockanalysis(cfg, datawhite);

cfg = [];
tlckwhite2 = ft_denoise_prewhiten(cfg, tlck, tlck);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%IM2GIF Convert a multiframe image to an animated GIF file
%
% Examples:
%   im2gif infile
%   im2gif infile outfile
%   im2gif(A, outfile)
%   im2gif(..., '-nocrop')
%   im2gif(..., '-nodither')
%   im2gif(..., '-ncolors', n)
%   im2gif(..., '-loops', n)
%   im2gif(..., '-delay', n)
%
% This function converts a multiframe image to an animated GIF.
%
% To create an animation from a series of figures, export to a multiframe
% TIFF file using export_fig, then convert to a GIF, as follows:
%
%    for a = 2 .^ (3:6)
%       peaks(a);
%       export_fig test.tif -nocrop -append
%    end
%    im2gif('test.tif', '-delay', 0.5);
%
%IN:
%   infile - string containing the name of the input image.
%   outfile - string containing the name of the output image (must have the
%             .gif extension). Default: infile, with .gif extension.
%   A - HxWxCxN array of input images, stacked along fourth dimension, to
%       be converted to gif.
%   -nocrop - option indicating that the borders of the output are not to
%             be cropped.
%   -nodither - option indicating that dithering is not to be used when
%               converting the image.
%   -ncolors - option pair, the value of which indicates the maximum number
%              of colors the GIF can have. This can also be a quantization
%              tolerance, between 0 and 1. Default/maximum: 256.
%   -loops - option pair, the value of which gives the number of times the
%            animation is to be looped. Default: 65535.
%   -delay - option pair, the value of which gives the time, in seconds,
%            between frames. Default: 1/15.

% Copyright (C) Oliver Woodford 2011

function im2gif(A, varargin)

% Parse the input arguments
[A, options] = parse_args(A, varargin{:});

if options.crop ~= 0
    % Crop
    A = crop_borders(A, A(ceil(end/2),1,:,1));
end

% Convert to indexed image
[h, w, c, n] = size(A);
A = reshape(permute(A, [1 2 4 3]), h, w*n, c);
map = unique(reshape(A, h*w*n, c), 'rows');
if size(map, 1) > 256
    dither_str = {'dither', 'nodither'};
    dither_str = dither_str{1+(options.dither==0)};
    if options.ncolors <= 1
        [B, map] = rgb2ind(A, options.ncolors, dither_str);
        if size(map, 1) > 256
            [B, map] = rgb2ind(A, 256, dither_str);
        end
    else
        [B, map] = rgb2ind(A, min(round(options.ncolors), 256), dither_str);
    end
else
    if max(map(:)) > 1
        map = double(map) / 255;
        A = double(A) / 255;
    end
    B = rgb2ind(im2double(A), map);
end
B = reshape(B, h, w, 1, n);

% Bug fix to rgb2ind
map(B(1)+1,:) = im2double(A(1,1,:));

% Save as a gif
imwrite(B, map, options.outfile, 'LoopCount', round(options.loops(1)), 'DelayTime', options.delay);
end

%% Parse the input arguments
function [A, options] = parse_args(A, varargin)
% Set the defaults
options = struct('outfile', '', ...
                 'dither', true, ...
                 'crop', true, ...
                 'ncolors', 256, ...
                 'loops', 65535, ...
                 'delay', 1/15);

% Go through the arguments
a = 0;
n = numel(varargin);
while a < n
    a = a + 1;
    if ischar(varargin{a}) && ~isempty(varargin{a})
        if varargin{a}(1) == '-'
            opt = lower(varargin{a}(2:end));
            switch opt
                case 'nocrop'
                    options.crop = false;
                case 'nodither'
                    options.dither = false;
                otherwise
                    if ~isfield(options, opt)
                        error('Option %s not recognized', varargin{a});
                    end
                    a = a + 1;
                    if ischar(varargin{a}) && ~ischar(options.(opt))
                        options.(opt) = str2double(varargin{a});
                    else
                        options.(opt) = varargin{a};
                    end
            end
        else
            options.outfile = varargin{a};
        end
    end
end

if isempty(options.outfile)
    if ~ischar(A)
        error('No output filename given.');
    end
    % Generate the output filename from the input filename
    [path, outfile] = fileparts(A);
    options.outfile = fullfile(path, [outfile '.gif']);
end

if ischar(A)
    % Read in the image
    A = imread_rgb(A);
end
end

%% Read image to uint8 rgb array
function [A, alpha] = imread_rgb(name)
% Get file info
info = imfinfo(name);
% Special case formats
switch lower(info(1).Format)
    case 'gif'
        [A, map] = imread(name, 'frames', 'all');
        if ~isempty(map)
            map = uint8(map * 256 - 0.5); % Convert to uint8 for storage
            A = reshape(map(uint32(A)+1,:), [size(A) size(map, 2)]); % Assume indexed from 0
            A = permute(A, [1 2 5 4 3]);
        end
    case {'tif', 'tiff'}
        A = cell(numel(info), 1);
        for a = 1:numel(A)
            [A{a}, map] = imread(name, 'Index', a, 'Info', info);
            if ~isempty(map)
                map = uint8(map * 256 - 0.5); % Convert to uint8 for storage
                A{a} = reshape(map(uint32(A{a})+1,:), [size(A) size(map, 2)]); % Assume indexed from 0
            end
            if size(A{a}, 3) == 4
                % TIFF in CMYK colourspace - convert to RGB
                if isfloat(A{a})
                    A{a} = A{a} * 255;
                else
                    A{a} = single(A{a});
                end
                A{a} = 255 - A{a};
                A{a}(:,:,4) = A{a}(:,:,4) / 255;
                A{a} = uint8(A(:,:,1:3) .* A{a}(:,:,[4 4 4]));
            end
        end
        A = cat(4, A{:});
    otherwise
        [A, map, alpha] = imread(name);
        A = A(:,:,:,1); % Keep only first frame of multi-frame files
        if ~isempty(map)
            map = uint8(map * 256 - 0.5); % Convert to uint8 for storage
            A = reshape(map(uint32(A)+1,:), [size(A) size(map, 2)]); % Assume indexed from 0
        elseif size(A, 3) == 4
            % Assume 4th channel is an alpha matte
            alpha = A(:,:,4);
            A = A(:,:,1:3);
        end
end
end

author: optimaltransport


text:
%% =======================================================================
%  ARDrone Simulation Example: Hovering and Position Control
%  =======================================================================
%
%  The simulation is used to validate the controller and guidance logic of
%  the ARDRone before flight testing. Control blocks are
%  exactly the same for both simulation and real-time Wi-Fi control.
%
%  Authors:
%       David Escobar Sanabria -> descobar@aem.umn.edu
%       Pieter J. Mosterman -> pieter.mosterman@mathworks.com
%  =======================================================================

%%
%  Cleaning workspace
bdclose all;
clear all;
clc

%%
% Adding ARDrone library path
addpath ../lib;
%% Simulation parameters

% Flight management system sample time. This is the sample time at which
% the control law is executed.
FMS.Ts = 0.065;

% Time delay due to communication between drone and host computer
timeDelay = FMS.Ts*4;

%% Vehicle model based on linear dynamics

% Loading state space representation of vehicle dynamics
setupARModel;

%%
% Loading list of waypoints
waypoints = getWaypoints() ;

%%
% Simulation time
simDT = 0.005 ;

%%
% Loading Simulink model of ARDrone
ARDroneHoverSim ;


text:
function [H] = rotate(f)

% ROTATE returns the homogenous coordinate transformation matrix
% corresponding to a rotation around the x, y and z-axis. The direction of
% the rotation is according to the right-hand rule.
%
% Use as
%   [H] = rotate(R)
% where
%   R       [rx, ry, rz] in degrees
%   H       corresponding homogenous transformation matrix
%
% Note that the order in which the rotations are performs matters. The
% rotation is first done around the z-axis, then the y-axis and finally the
% x-axis.

% Copyright (C) 2000-2005, Robert Oostenveld
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% 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$

if numel(f)~=3
  ft_error('incorrect input vector');
end

% convert degrees to radians
f = f*pi/180;

% get the individual angles (in radians)
rx = f(1);
ry = f(2);
rz = f(3);

% precompute the sin/cos values of the angles
cX = cos(rx);
cY = cos(ry);
cZ = cos(rz);
sX = sin(rx);
sY = sin(ry);
sZ = sin(rz);

% according to Roger Woods' http://bishopw.loni.ucla.edu/AIR5/homogenous.html
% it should be this, but I cannot reproduce his rotation matrix
% H = eye(4,4);
% H(1,1) = cZ*cY + sZ*sX*sY;
% H(1,2) = sZ*cY - cZ*sX*sY;
% H(1,3) =            cX*sY;
% H(2,1) = -sZ*cX;
% H(2,2) =  cZ*cX;
% H(2,3) =     sX;
% H(3,1) =  sZ*sX*cY - cZ*sY;
% H(3,2) = -cZ*sX*cY - sZ*sY;
% H(3,3) =             cX*cY;

% instead, the following rotation matrix does work according my
% expectations. It rotates according to the right hand rule and first
% rotates around z, then y and then x axis
H = [
           cZ*cY,          -sZ*cY,              sY,               0
  cZ*sY*sX+sZ*cX, -sZ*sY*sX+cZ*cX,          -cY*sX,               0
 -cZ*sY*cX+sZ*sX,  sZ*sY*cX+cZ*sX,           cY*cX,               0
               0,               0,               0,               1
];

if 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following code can be used to construct the combined rotation matrix
% for either xyz or zyx ordering (using the MATLAB symbolic math toolbox)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  syms sX sY sZ cX cY cZ
  % this is for only rotating around x
  Rx = [
    1   0    0   0
    0   cX  -sX  0
    0   sX   cX  0
    0   0    0   1
    ];
  % this is for only rotating around y
  Ry = [
    cY  0   sY  0
    0   1   0   0
    -sY  0   cY  0
    0   0   0   1
    ];
  % this is for only rotating around z
  Rz = [
    cZ -sZ  0   0
    sZ  cZ  0   0
    0   0   1   0
    0   0   0   1
    ];
  % combine them
  Rzyx = Rz * Ry * Rx  % rotate around x, y, then z
  Rxyz = Rx * Ry * Rz  % rotate around z, y, then x
end

author: spm


text:
function [d,c]=runenc(x);
ind=1;
d(ind)=x(1);
c(ind)=1;

for i=2 :length(x)
    if x(i-1)==x(i)
       c(ind)=c(ind)+1;
    else ind=ind+1;
         d(ind)=x(i);
         c(ind)=1;
    end
end


text:
function full_mem(layer_idx)
    global mem;

end

author: jimmy-ren


text:
function [EEGout cfg] = pop_sim_varmodel(EEG,typeproc,varargin)
%
% Simulate a Vector Autoregressive dynamical system
%
% Input:
% Optional:
%
%   EEG:            existing EEG dataset (configs will be retrived from
%                   here)
%   typeproc:       if 'nogui' don't generate GUI
%
%   <'Name',value> pairs as defined in sim_varmodel()
%
% Output:
%
%   EEGout:         Simulated EEG structure(s).
%                   Optionally this may be an array of
%                   EEG structs with the second struct
%                   being the ground truth model.
%   cfg:            Argument specification structure.
%
%
% See Also: sim_varmodel()
%
% References:
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual.
%   Available at: http://www.sccn.ucsd.edu/wiki/Sift
%
% Author: Tim Mullen 2013, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu
%

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

if nargin<1
    EEG = [];
end
if nargin<2
    typeproc = 0;
end
EEGout = EEG;
cfg    = [];

fcnName     = strrep(mfilename,'pop_','');
fcnHandle   = str2func(fcnName);

if isempty(hlp_checkeegset(EEG,{'cat'})) && isfield(EEG.CAT.configs,fcnName)
    % get default configuration (from prior use) and merge with varargin
    varargin = [hlp_struct2varargin(EEG.CAT.configs.(fcnName)) varargin];
end

if strcmpi(typeproc,'nogui')
    % get the default config from function and overload supplied args
    cfg = arg_tovals(arg_report('rich',fcnHandle,varargin),false);
else
    % render the GUI
    [PGh figh] = feval(['gui_' fcnName],varargin{:});

    if isempty(PGh)
        % user chose to cancel
        return;
    end

    % get the specification of the PropertyGrid
    ps = PGh.GetPropertySpecification;
    cfg = arg_tovals(ps,false);
end

drawnow;

if ~cfg.makeEEGset.arg_selection
    error('SIFT:sim_varmodel',['If using pop_' fcnName '(), you must enable the BuildEEGLABStructure option.\n' ...
                               'Use ' fcnName '() from the command-line to return a raw dataset']);
end

if strcmpi(typeproc,'cfg_only')
    return;
end

% execute the low-level function
[EEGout GroundTruth] = feval(fcnHandle,cfg);

if ~isempty(GroundTruth)
    EEGout = eeg_store(EEGout,GroundTruth,2);
end
if ~isempty(cfg)
    for k=1:length(EEGout)
        % store the configuration structure
        EEGout(k).CAT.configs.(fcnName) = cfg;
    end
end

repo: REST


text:
% Wrapper for vl_nntranslate2D block
% inputs{1} :   X     : 1 x 2 x n x b
% inputs{2} :   T     : 1 x 2 x 1 x b
% outputs{1}:   y     : 1 x 2 x n x b

classdef translate2D < dagnn.Layer
    methods
        function outputs = forward(~, inputs, ~)

            useGPU = isa(inputs{1}, 'gpuArray');
            if useGPU
                outputs{1} = gpuArray ( vl_nntranslate2D(gather(inputs{1}), gather(inputs{2})) );
            else
                outputs{1} = vl_nntranslate2D(inputs{1}, inputs{2});
            end

        end

        function [derInputs, derParams] = backward(~, inputs, ~, derOutputs)

            useGPU = isa(inputs{1}, 'gpuArray');
            if useGPU
                [y,dsdy] = vl_nntranslate2D(gather(inputs{1}),gather(inputs{2}), gather(derOutputs{1}));
                derInputs = {gpuArray(y),gpuArray(dsdy)};
            else
                [y,dsdy] = vl_nntranslate2D(inputs{1},inputs{2}, derOutputs{1});
                derInputs = {y,dsdy};
            end

            derParams = {};
        end

        function outputSizes = getOutputSizes(~, inputSizes)
            outputSizes = inputSizes{1};
        end

        function obj = translate2D(varargin)
            obj.load(varargin);
        end
    end
end

author: anilbas


text:
% distance to closest region of interest
function [data,units,mind] = compute_closestroi2(trx,n,dosave_d)

if nargin < 3,
  dosave_d = true;
end

flies = trx.exp2flies{n};
nflies = numel(flies);
closestfly = cell(1,nflies);
mind = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);
  ROIdata=trx.roi2{n}.data{fly};
  mind{i} = inf(1,trx(fly).nframes);
  closestfly{i} = nan(1,trx(fly).nframes);
  for j = 1:size(ROIdata,1),
    fn = sprintf('dist2roi2_%d',j);
    dcurr = trx(fly).(fn);
    idx = dcurr < mind{i};
    mind{i}(idx) = dcurr(idx);
    closestfly{i}(idx) = j;
  end
end

% so that we don't compute dcenter twice
if dosave_d,
  data = mind; %#ok<NASGU>
  units = parseunits('mm'); %#ok<NASGU>
  filename = trx.GetPerFrameFile('mindist2roi2',n);
  save(filename,'data','units');
end

data = closestfly;
units = parseunits('unit');

author: kristinbranson

repo: JAABA


text:
function sF = unimodal(varargin)
% defines a unimodal spherical function
%
% Syntax
%
%   sF = S2Fun.unimodal
%   sF = S2Fun.unimodal('halfwidth',10*degree)
%
%   v = vector3d(1,1,1)
%   psi = S2DeLaValleePoussinKernel('halfwidth',20*degree)
%   sF = S2Fun.unimodal(v,psi)
%
% Input
%  v - symmetry axis @vector3d
%  psi - @S2Kernel
%
% Output
%  sF - @S2FunHarmonic
%

% extract kernel
psi = S2DeLaValleePoussinKernel('halfwidth',get_option(varargin,'halfwidth',25*degree));
psi = getClass(varargin,'S2Kernel',psi);

% define a radially symmetric function
bw = psi.bandwidth;

% 1 3 5  7  9 11 13
% 1 3 7 13 21 31 43

% the indice of diagonal
l = 0:bw; l = 1+l.^2+l;

f_hat = zeros((2*bw+1)^2,1);
f_hat(l) = psi.A ./ sqrt(2*l+1);

sF = S2FunHarmonic(f_hat);

v = getClass(varargin,'vector3d',vector3d.Z);

if angle(v,vector3d.Z) > 0
  rot = rotation.map(v,vector3d.Z);
  sF = rot * sF;
end

end

author: mtex-toolbox

repo: mtex

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


text:
% eeglab
try
clear all
close all
clc
catch ME
    disp(ME.message)
clear all
close all
end

%%
% eeglabPath = fileparts(which('eeglab'));
mobilabPath = 'C:\DevelCore\mobilab'; %'/Users/timmullen/Documents/WORK/Toolboxes/mobilab'; %[eeglabPath filesep 'plugins' filesep 'mobilab'];
% addpath(genpath(mobilabPath));

templateFile = [mobilabPath filesep 'data' filesep 'headModelColin27_11997.mat'];
standardMontage = [mobilabPath filesep 'data' filesep 'standard_1020.elc'];
template = load([mobilabPath filesep 'data' filesep 'headModelColin27_11997.mat']);
individualMontage = [mobilabPath filesep 'data' filesep 'personA_chanlocs.sfp'];
surfFile = 'demo_Colin27_11997.mat';
surfData = template.surfData;
save(surfFile,'surfData');

%% display template head model: Colin27 + aal atlas + 10/20 montage
[elec,label] = readMontage(standardMontage);
hmObj = headModel('surfaces',surfFile,'atlas',template.atlas,'fiducials',template.fiducials,'channelSpace',elec,'label',label);
plotHeadModel(hmObj);

%% warping individual channel space to template
individualSurfaces = 'demo_Colin27_11997.mat';

hmObj = headModel(individualMontage);
plotMontage(hmObj);

aff = hmObj.warpChannelSpace2Template(templateFile,individualSurfaces,'affine');
plotHeadModel(hmObj);

%% warping template to individual channel space
individualSurfaces = 'warped_demo_Colin27_11997.mat';  % name of the output file

hmObj = headModel(individualMontage);
plotMontage(hmObj);

hmObj.warpTemplate2channelSpace(templateFile,individualSurfaces);
plotHeadModel(hmObj);

%% solving the forward problem with OpenMEEG
conductivity = [0.33 0.022 0.33]; % brain and scalp = 0.33 S/m, skull = 0.022 S/m; these conductivies were taken from
                                  % Valdes-Hernandez et al., 2009, Oostendrop TF, 2000; Wendel and Malmivuo, 2006
normal2surface = true;
hmObj.computeLeadFieldBEM(conductivity,normal2surface);

%% solving the forward problem with NFT
% conductivity = [0.33 0.022 0.33]; % brain and scalp = 0.33 S/m, skull = 0.022 S/m; these conductivies were taken from
%                                   % Valdes-Hernandez et al., 2009, Oostendrop TF, 2000; Wendel and Malmivuo, 2006
% normal2surface = true;
% hmObj.computeLeadFieldBEM_NFT(conductivity,normal2surface);

%% solving the inverse problem with sLORETA
% K: lead field matrix
% L: Laplaciian operator
% rmIndices: indices to be removed (the Thalamus)
% surfData(3).vertices: source space

brainStructsToKeep = {'Occipital_Sup_L','Occipital_Mid_L','Occipital_Inf_L'};
brainStructsToRemove = setdiff_bc(hmObj.atlas.label,brainStructsToKeep); %{'Thalamus_L','Thalamus_R'};
[sourceSpace,K,L,rmIndices] = getSourceSpace4PEB(hmObj,brainStructsToRemove);
load(hmObj.surfaces)
n = size(surfData(3).vertices,1);
J = zeros(n,1);
Jtrue = J;
Jest = J;
ind = setdiff_bc(1:n,rmIndices);

% simulating some Gaussian sources
x0 = [-81.6328 17.9887 93.8088];  % occipital l
x1 = [30.5384 -12.5882 117.1195]; % frontal superior r
x2 = [0.8625 -58.2585 78.9345];   % supra marginal r
x3 = [25.9583 48.3573 35.8547];   % insula l

d = sqrt(sum((surfData(3).vertices(ind,:) - ones(length(ind),1)*x0).^2,2));
Jtrue(ind) = normpdf(d,0,10);
Vtrue = K*Jtrue(ind);

nlambdas = 100;
plotGCV  = true;
threshold = [25 85]; % threshold = [1 99]; threshold = [];
Jt = inverseSolutionLoreta(Vtrue,K,L,nlambdas,plotGCV,threshold);  % Jt contains the solution for the source potentials within the chosen space, the remaining elements are zero
Jest(ind) = Jt;
Vest = K*Jest(ind);
hmObj.plotOnModel(Jtrue,Vtrue,'True source');
hmObj.plotOnModel(Jest,Vest,'Estimated source');

disp('adding noise');
% adding noise
snr = 5;
vn = std(Vtrue)/snr*randn(length(Vtrue),1);
Vtrue_noise = Vtrue + vn;

Jt = inverseSolutionLoreta(Vtrue_noise,K,L,nlambdas,plotGCV,threshold);
Jest(ind) = Jt;
Vest = K*Jest(ind);
disp('plotting');
hmObj.plotOnModel(Jtrue,Vtrue,'True source (noisy)');
hmObj.plotOnModel(Jest,Vest,'Estimated source (noisy)');


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ldldown.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [L,d]=ldldown(L,d,j)
% downdates LDL^T factorization when j-th row and column are replaced
% by j-th unit vector
%
% d contains diag(D) and is assumed positive
%
function [L,d]=ldldown(L,d,j);

n=size(d,1);

test=0;
if test,
  disp('enter ldldown')
  A=L*diag(d)*L';A(:,j)=zeros(n,1);A(j,:)=zeros(1,n);A(j,j)=1;
end;

if j<n,
  I=1:j-1;K=j+1:n;
  [LKK,d(K)]=ldlrk1(L(K,K),d(K),d(j),L(K,j));
  % work around expensive sparse L(K,K)=LKK
  L=[L(I,:);
     sparse(1,n);
     L(K,I),sparse(n-j,1),LKK];
  L(j,j)=1;
else
  L(n,1:n-1)=sparse(1,n-1);
end;
d(j)=1;

if test,
  A1=L*diag(d)*L',A
  quot=norm(A1-A,1)/norm(A,1),
  disp('leave ldldown')
end;

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function labels = context_labels(cls, ds, train_set, train_year)
% Get classification training labels for training the context rescoring
% classifier.
%   labels = context_labels(cls, ds, train_set, train_year)
%
% Return value
%   labels      Binary labels {-1,+1} for each detection in boxes
%
% Arguments
%   cls         Object class
%   ds          Detections
%   train_set   Training dataset
%   train_year  Training dataset year

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
% Copyright (C) 2008, 2009, 2010 Pedro Felzenszwalb, Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

conf = voc_config('pascal.year', train_year);
cachedir = conf.paths.model_dir;
VOCopts  = conf.pascal.VOCopts;

try
  load([cachedir cls '_context_labels_' train_set '_' train_year]);
catch
  fprintf('Constructing training labels (this will take a little while)...\n');
  [gt, npos] = get_ground_truth(cls, train_set, train_year);
  [gtids, t] = textread(sprintf(VOCopts.imgsetpath, train_set),'%s %d');

  labels = cell(length(gtids),1);

  L = 0;
  for i = 1:length(gtids)
    L = L + size(ds{i},1);
  end

  detections = zeros(L,7);
  I = 1;
  for i = 1:length(gtids)
    if ~isempty(ds{i})
      l = size(ds{i},1);
      % Detection scores
      detections(I:I+l-1,1) = ds{i}(:,end);
      % Detection windows
      detections(I:I+l-1,2:5) = ds{i}(:,1:4);
      % The image (i) the detections came from
      detections(I:I+l-1,6) = i;
      % The index in ds{i} for each detection
      detections(I:I+l-1,7) = 1:l;
      labels{i} = zeros(l,1);
      I = I+l;
    else
      labels{i} = [];
    end
  end

  [sc, si] = sort(-detections(:,1));
  ids = detections(si,6);
  idx = detections(si,7);
  BB = detections(si,2:5)';

  % Adapted from the VOCdevkit m-file VOCevaldet.m

  % assign detections to ground truth objects
  nd=length(si);
  for d=1:nd
    % find ground truth image
    i=ids(d);

    % assign detection to ground truth object if any
    bb=BB(:,d);
    ovmax=-inf;
    for j=1:size(gt(i).boxes,2)
      bbgt=gt(i).boxes(:,j);
      bi=[max(bb(1),bbgt(1)) ; max(bb(2),bbgt(2)) ; min(bb(3),bbgt(3)) ; min(bb(4),bbgt(4))];
      iw=bi(3)-bi(1)+1;
      ih=bi(4)-bi(2)+1;
      if iw>0 & ih>0
        % compute overlap as area of intersection / area of union
        ua=(bb(3)-bb(1)+1)*(bb(4)-bb(2)+1)+...
           (bbgt(3)-bbgt(1)+1)*(bbgt(4)-bbgt(2)+1)-...
           iw*ih;
        ov=iw*ih/ua;
        if ov>ovmax
          ovmax=ov;
          jmax=j;
        end
      end
    end
    % assign detection as true positive/don't care/false positive
    if ovmax>=VOCopts.minoverlap
      if ~gt(i).diff(jmax)
        if ~gt(i).det(jmax)
          % True positive
          gt(i).det(jmax)=true;
          labels{i}(idx(d)) = 1;
        else
          % false positive (multiple detection)
          labels{i}(idx(d)) = -1;
        end
      else
        labels{i}(idx(d)) = 1;   % difficult
      end
    else
      % false positive (low overlap)
      labels{i}(idx(d)) = -1;
    end
  end
  save([cachedir cls '_context_labels_' train_set '_' train_year], 'labels');
  fprintf('done!\n');
end

function [gt, npos] = get_ground_truth(cls, dataset, year)
% Load and cache ground-truth annontation data.
% Most of this code is borrowed from the PASCAL devkit.

conf = voc_config('pascal.year', year);
cachedir = conf.paths.model_dir;
VOCopts  = conf.pascal.VOCopts;
VOCyear  = conf.pascal.year;

try
  load([cachedir cls '_gt_anno_' dataset '_' VOCyear]);
catch
  % load ground truth objects
  [gtids, t] = textread(sprintf(VOCopts.imgsetpath,dataset),'%s %d');
  npos = 0;
  for i = 1:length(gtids)
    % display progress
    tic_toc_print('%s: loading ground truth: %d/%d\n',cls,i,length(gtids));

    % read annotation
    rec = PASreadrecord(sprintf(VOCopts.annopath,gtids{i}));

    % extract objects of class
    clsinds = strmatch(cls,{rec.objects(:).class},'exact');
    gt(i).boxes = cat(1,rec.objects(clsinds).bbox)';
    gt(i).diff = [rec.objects(clsinds).difficult];
    gt(i).det = false(length(clsinds),1);
    npos = npos+sum(~gt(i).diff);
  end
  save([cachedir cls '_gt_anno_' dataset '_' VOCyear], 'gt', 'npos');
end

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
clc;clear;
models = {'','useCN10_ECSSD'};

Num = length(models);
method_col = {[1,0,0],[0,1,0],[0.5 0.5 0.5],[136 0 21]/255,[255 127 39]/255,...
    [0 162 232]/255,[163 73 164]/255,[1 0 1]};
for i = 1 : Num
    load(models{i});
    mmFmeasure = (1+0.3).*mPre.*mRecall./(0.3.*mPre+mRecall);
    plot([0:255/20:255],mmFmeasure,'Color',method_col{i},'LineWidth',2);
    axis([0 255 0 0.9]);
    grid on;
    hold on;
end
legend( models );

repo: Saliency-Evaluation-Toolbox


text:
function [center,radius] = minboundcircle(x,y,hullflag)
% minboundcircle: Compute the minimum radius enclosing circle of a set of (x,y) pairs
% usage: [center,radius] = minboundcircle(x,y,hullflag)
%
% arguments: (input)
%  x,y - vectors of points, describing points in the plane as
%        (x,y) pairs. x and y must be the same size. If x and y
%        are arrays, they will be unrolled to vectors.
%
%  hullflag - boolean flag - allows the user to disable the
%        call to convhulln. This will allow older releases of
%        matlab to use this code, with a possible time penalty.
%        It also allows minboundellipse to call this code
%        efficiently.
%
%        hullflag = false --> do not use the convex hull
%        hullflag = true  --> use the convex hull for speed
%
%        default: true
%
%
% arguments: (output)
%  center - 1x2 vector, contains the (x,y) coordinates of the
%        center of the minimum radius enclosing circle
%
%  radius - scalar - denotes the radius of the minimum
%        enclosing circle
%
%
% Example usage:
%   x = randn(50000,1);
%   y = randn(50000,1);
%   tic,[c,r] = minboundcircle(x,y);toc
%
%   Elapsed time is 0.171178 seconds.
%
%   c: [-0.2223 0.070526]
%   r: 4.6358
%
%
% See also: minboundrect
%
%
% Author: John D'Errico
% E-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 1/10/07

% default for hullflag
if (nargin<3) || isempty(hullflag)
  hullflag = true;
elseif ~islogical(hullflag) && ~ismember(hullflag,[0 1])
  error 'hullflag must be true or false if provided'
end

% preprocess data
x=x(:);
y=y(:);

% not many error checks to worry about
n = length(x);
if n~=length(y)
  error 'x and y must be the same sizes'
end

% start out with the convex hull of the points to
% reduce the problem dramatically. Note that any
% points in the interior of the convex hull are
% never needed.
if hullflag && (n>3)
  edges = convhulln([x,y],{'QJ' 'Pp'});

  % list of the unique points on the convex hull itself
  % convhulln returns them as edges
  edges = unique(edges(:));

  % exclude those points inside the hull as not relevant
  x = x(edges);
  y = y(edges);

end

% now we must find the enclosing circle of those that
% remain.
n = length(x);

% special case small numbers of points. If we trip any
% of these cases, then we are done, so return.
switch n
  case 0
    % empty begets empty
    center = [];
    radius = [];
    return
  case 1
    % with one point, the center has radius zero
    center = [x,y];
    radius = 0;
    return
  case 2
    % only two points. center is at the midpoint
    center = [mean(x),mean(y)];
    radius = norm([x(1),y(1)] - center);
    return
  case 3
    % exactly 3 points
    [center,radius] = enc3(x,y);
    return
end

% more than 3 points.

% Use an active set strategy.
aset = 1:3; % arbitrary, but quite adequate
iset = 4:n;

% pick a tolerance
tol = 10*eps*(max(abs(mean(x) - x)) + max(abs(mean(y) - y)));

flag = true;
while flag
  % get the enclosing circle for the current set
  [center,radius] = enc3(x(aset),y(aset));

  % are all the inactive set points inside the circle?
  r = sqrt((x(iset) - center(1)).^2 + (y(iset) - center(2)).^2);
  [rmax,k] = max(r);
  if (rmax - radius) <= tol
    % the active set enclosing circle also enclosed
    % all of the inactive points.
    flag = false;
  else
    % it must be true that we can replace one member of aset
    % with iset(k). Which one?
    s1 = [aset([2 3]),iset(k)];
    [c1,r1] = enc3(x(s1),y(s1));
    if (norm(c1 - [x(aset(1)),y(aset(1))]) <= r1)
      center = c1;
      radius = r1;

      % update the active/inactive sets
      swap = aset(1);
      aset = [iset(k),aset([2 3])];
      iset(k) = swap;

      % bounce out to the while loop
      continue
    end
    s1 = [aset([1 3]),iset(k)];
    [c1,r1] = enc3(x(s1),y(s1));
    if (norm(c1 - [x(aset(2)),y(aset(2))]) <= r1)
      center = c1;
      radius = r1;

      % update the active/inactive sets
      swap = aset(2);
      aset = [iset(k),aset([1 3])];
      iset(k) = swap;

      % bounce out to the while loop
      continue
    end
    s1 = [aset([1 2]),iset(k)];
    [c1,r1] = enc3(x(s1),y(s1));
    if (norm(c1 - [x(aset(3)),y(aset(3))]) <= r1)
      center = c1;
      radius = r1;

      % update the active/inactive sets
      swap = aset(3);
      aset = [iset(k),aset([1 2])];
      iset(k) = swap;

      % bounce out to the while loop
      continue
    end

    % if we get through to this point, then something went wrong.
    % Active set problem. Increase tol, then try again.
    tol = 2*tol;

  end

end

% =======================================
%  begin subfunction
% =======================================
function [center,radius] = enc3(X,Y)
% minimum radius enclosing circle for exactly 3 points
%
% x, y are 3x1 vectors

% convert to complex
xy = X + sqrt(-1)*Y;

% just in case the points are collinear or nearly so, get
% the interpoint distances, and test the farthest pair
% to see if they work.
Dij = @(XY,i,j) abs(XY(i) - XY(j));
D12 = Dij(xy,1,2);
D13 = Dij(xy,1,3);
D23 = Dij(xy,2,3);

% Find the most distant pair. Test if their circumcircle
% also encloses the third point.
if (D12>=D13) && (D12>=D23)
  center = (xy(1) + xy(2))/2;
  radius = D12/2;
  if abs(center - xy(3)) <= radius
    center = [real(center),imag(center)];
    return
  end
elseif (D13>=D12) && (D13>=D23)
  center = (xy(1) + xy(3))/2;
  radius = D13/2;
  if abs(center - xy(2)) <= radius
    center = [real(center),imag(center)];
    return
  end
elseif (D23>=D12) && (D23>=D13)
  center = (xy(2) + xy(3))/2;
  radius = D23/2;
  if abs(center - xy(1)) <= radius
    center = [real(center),imag(center)];
    return
  end
end

% if we drop down to here, then the points cannot
% be collinear, so the resulting 2x2 linear system
% of equations will not be singular.
A = 2*[X(2)-X(1), Y(2)-Y(1); X(3)-X(1), Y(3)-Y(1)];
rhs = [X(2)^2 - X(1)^2 + Y(2)^2 - Y(1)^2; ...
       X(3)^2 - X(1)^2 + Y(3)^2 - Y(1)^2];

center = (A\rhs)';
radius = norm(center - [X(1),Y(1)]);


text:
function [db,f]=lpccc2db(cc,np,nc,c0)
%LPCCC2DB Convert complex cepstrum to dB power spectrum DB=(CC,NP,NC)
%
%  Inputs: cc(nf,n)     Complex ceptral coefficients excluding c(0), one frame per row
%          np           Size of output spectrum is np+1 [n]
%                       Alternatively, np can be a vector of output frequencies in the range 0 to 0.5
%          nc           Highest cepstral coefficient to use [np or, if np is a vector, n]
%                       Set nc=-1 to use n coefficients
%          c0(nf,1)     Cepstral coefficient cc(0) [0]
%
% Outputs: db(nf,np+2)  Power spectrum from DC to Nyquist in dB
%          f(1,np+2)    Normalized frequencies (0 to 0.5)
%
% The "complex cepstral coefficients", cc(n), are the inverse discrete-time Fourier transform
% of the log of the complex-valued spectrum. The cc(n) are real-valued and, for n<0, cc(n)=0.
% The "real cepstral coeffcients", rc(n), are the inverse discrete-time Fourier transform
% of the log of the magnitude spectrum; rc(0)=cc(0) and rc(n)=0.5*cc(n) for n~=0.
% For highest speed, choose np to be a power of 2.

%      Copyright (C) Mike Brookes 1998-2014
%      Version: $Id: lpccc2db.m 5025 2014-08-22 17:07:24Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nf,mc]=size(cc);
if nargin<2 || ~numel(np)
    if nargout
        np=mc;
    else
        np=128;
    end
end
k=10/log(10);
if nargin>=3 && numel(nc)==1 && nc==-1 nc=mc; end
if nargin<4 || ~numel(c0) c0=zeros(nf,1); end
if numel(np)>1 || np(1)<1
    if nargin<3 || ~numel(nc) nc=mc; end
    f=np(:)';
    if nc==mc
        db=k*(2*[c0 cc]*cos(2*pi*(0:mc)'*f));
    else
        db=k*(2*[c0 lpccc2cc(cc,nc)]*cos(2*pi*(0:nc)'*f));
    end
else
    if nargin<3 || ~numel(nc) nc=np; end
    if nc==mc
        db=k*(2*real(rfft([c0 cc].',2*np).'));
    else
        db=k*(2*real(rfft([c0 lpccc2cc(cc,nc)].',2*np).'));
    end
    f=linspace(0,0.5,np+1);
end
if ~nargout
    plot(f,db.');
    xlabel('Normalized frequency f/f_s');
    ylabel('Gain (dB)');
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function tests = test_rotateStack90(varargin)
%TEST_ROTATESTACK90  One-line description here, please.
%   output = test_rotateStack90(input)
%
%   Example
%   test_rotateStack90
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-05-18,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

tests = functiontests(localfunctions);

function testRotateX(testCase)

img = createBasicTestImage;

rotX1th = cat(3, [5 6;1 2], [7 8;3 4]);
imgRX1 = rotateStack90(img, 2, 1);
assertEqual(testCase, rotX1th, imgRX1);
imgRX1 = rotateStack90(img, 'x', 1);
assertEqual(testCase, rotX1th, imgRX1);

rotX2th = cat(3, [7 8;5 6], [3 4;1 2]);
imgRX2 = rotateStack90(img, 2, 2);
assertEqual(testCase, rotX2th, imgRX2);
imgRX2 = rotateStack90(img, 'x', 2);
assertEqual(testCase, rotX2th, imgRX2);

rotX3th = cat(3, [3 4;7 8], [1 2;5 6]);
imgRX3 = rotateStack90(img, 2, 3);
assertEqual(testCase, rotX3th, imgRX3);
imgRX3 = rotateStack90(img, 'x', 3);
assertEqual(testCase, rotX3th, imgRX3);

function testRotateY(testCase)

img = createBasicTestImage;

rotY1th = cat(3, [2 6;4 8], [1 5;3 7]);
imgRY1 = rotateStack90(img, 1, 1);
assertEqual(testCase, rotY1th, imgRY1);
imgRY1 = rotateStack90(img, 'y', 1);
assertEqual(testCase, rotY1th, imgRY1);

rotY2th = cat(3, [6 5;8 7], [2 1;4 3]);
imgRY2 = rotateStack90(img, 1, 2);
assertEqual(testCase, rotY2th, imgRY2);
imgRY2 = rotateStack90(img, 'y', 2);
assertEqual(testCase, rotY2th, imgRY2);

rotY3th = cat(3, [5 1;7 3], [6 2;8 4]);
imgRY3 = rotateStack90(img, 1, 3);
assertEqual(testCase, rotY3th, imgRY3);
imgRY3 = rotateStack90(img, 'y', 3);
assertEqual(testCase, rotY3th, imgRY3);

function testRotateZ(testCase)

img = createBasicTestImage;

rotZ1th = cat(3, [3 1;4 2], [7 5;8 6]);
imgRZ1 = rotateStack90(img, 3, 1);
assertEqual(testCase, rotZ1th, imgRZ1);
imgRZ1 = rotateStack90(img, 'z', 1);
assertEqual(testCase, rotZ1th, imgRZ1);

rotZ2th = cat(3, [4 3;2 1], [8 7;6 5]);
imgRZ2 = rotateStack90(img, 3, 2);
assertEqual(testCase, rotZ2th, imgRZ2);
imgRZ2 = rotateStack90(img, 'z', 2);
assertEqual(testCase, rotZ2th, imgRZ2);

rotZ3th = cat(3, [2 4;1 3], [6 8;5 7]);
imgRZ3 = rotateStack90(img, 3, 3);
assertEqual(testCase, rotZ3th, imgRZ3);
imgRZ3 = rotateStack90(img, 'z', 3);
assertEqual(testCase, rotZ3th, imgRZ3);

function testRotateGrayscaleY(testCase)

img = createTestImageGray;
dim = size(img);

% rotate around Y-axis
img2 = rotateStack90(img, 1, 1);

% check dimension
assertEqual(testCase, dim([1 3 2]), size(img2));

function testRotateGrayscaleX(testCase)

img = createTestImageGray;
dim = size(img);

% rotate around X-axis
img2 = rotateStack90(img, 2, 1);
assertEqual(testCase, dim([3 2 1]), size(img2));

function testRotateGrayscaleZ(testCase)

img = createTestImageGray;
dim = size(img);

% rotate around Z-axis
img2 = rotateStack90(img, 3, 1);
assertEqual(testCase, dim([2 1 3]), size(img2));

function testColorImage(testCase)

% Create rgb image with non equal size in each dimension
lx = 1:50;
ly = 1:52;
lz = 1:54;
r = uint8(discreteBall(lx, ly, lz, [20 30 30], 15)*255);
g = uint8(discreteBall(lx, ly, lz, [30 20 30], 15)*255);
b = uint8(discreteBall(lx, ly, lz, [30 30 20], 15)*255);
img = imMergeChannels(r, g, b);

% basic checks
assertEqual(testCase, 'uint8', class(img));
assertEqual(testCase, [52 50 3 54], size(img));

% dimension of base image, without color
dim = size(r);

% rotate around Y-axis
img2 = rotateStack90(img, 1, 1);
assertEqual(testCase, [dim(1) dim(3) 3 dim(2)], size(img2));
assertEqual(testCase, img(1,1,:,1), img2(1, end,:,1));
assertEqual(testCase, img(1,1,:,end), img2(1, 1,:,1));

% rotate around X-axis
img2 = rotateStack90(img, 2, 1);
assertEqual(testCase, [dim(3) dim(2) 3 dim(1)], size(img2));

% rotate around Z-axis
img2 = rotateStack90(img, 3, 1);
assertEqual(testCase, [dim(2) dim(1) 3 dim(3)], size(img2));

function img = createBasicTestImage
img = cat(3, [1 2;3 4], [5 6;7 8]);

function img = createTestImageGray
% create a 3D gray-scale test image
%  ______
% |1    2|\
% | 5    |6|
% |3____4| |
%  \7____\8|
%
lx = 1:10;
ly = 30:5:60;
lz = 50:10:200;
[x, y, z] = meshgrid(lx, ly, lz);
img = uint8(x+y+z);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function attr = export_module(moduleName, varargin)
%"export_module"
%   Returns a populated Java attr representing a requested module,
%   given its moduleName string and the data required to populate the
%   fields.
%
%   Currently supported module names and their parameters are:
%       'patient'               scanS
%       'general_study'         scanS
%       'general_series'        scanS, doseS, structureS
%       'frame_of_reference'    scanS
%       'general_equipment'     'scan', scanS
%                               'structures', structureS
%       'general_image'         scanInfoS, scanS
%       'image_plane'           scanInfoS, scanS
%       'image_pixel'           scanInfoS, scanS
%       'ct_image'              scanInfoS, scanS
%       'sop_common'            'scan', scanS
%                               'structures', structureS
%       'structure_set'         structureS
%       'roi_contour'           structureS
%       'rt_series'             'structure', structureS
%       'rt_roi_observations'   structureS
%       'rt_dose'
%       'rt_dvh'
%
%JRA 06/19/06
%NAV 07/19/16 updated to dcm4che3
%    Used addAll over copyTo
%
%Usage:
%   attr = export_module(moduleName, varargin)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

%convert to dcm4che3 attribute
attr = javaObject('org.dcm4che3.data.Attributes');

switch lower(moduleName)
    %Get the top level tags used in this module.
    case 'patient'
        tagS = patient_module_tags;
        export_function = @export_patient_module_field;
    case 'general_study'
        tagS = general_study_module_tags;
        export_function = @export_general_study_module_field;
    case 'general_series'
        tagS = general_series_module_tags;
        export_function = @export_general_series_module_field;
    case 'frame_of_reference'
        tagS = frame_of_reference_module_tags;
        export_function = @export_frame_of_reference_module_field;
    case 'general_equipment'
        tagS = general_equipment_module_tags;
        export_function = @export_general_equipment_module_field;
    case 'content_identification'
        tagS = content_identification_module_tags;
        export_function = @export_content_identification_module_field;
    case 'general_image'
        tagS = general_image_module_tags;
        export_function = @export_general_image_module_field;
    case 'image_plane'
        tagS = image_plane_module_tags;
        export_function = @export_image_plane_module_field;
    case 'image_pixel'
        tagS = image_pixel_module_tags;
        export_function = @export_image_pixel_module_field;
    case 'ct_image'
        tagS = CT_image_module_tags;
        export_function = @export_CT_image_module_field;
    case 'pt_image'
        tagS = PT_image_module_tags;
        export_function = @export_PT_image_module_field;
    case 'mr_image'
        tagS = MR_image_module_tags;
        export_function = @export_MR_image_module_field;
    case 'sop_common'
        tagS = SOP_common_module_tags;
        export_function = @export_SOP_common_module_field;
    case 'structure_set'
        tagS = structure_set_module_tags;
        export_function = @export_structure_set_module_field;
    case 'roi_contour'
        tagS = ROI_contour_module_tags;
        export_function = @export_roi_contour_module_field;
    case 'rt_series'
        tagS = RT_series_module_tags;
        export_function = @export_rt_series_module_field;
    case 'rt_roi_observations'
        tagS = RT_ROI_observations_module_tags;
        export_function = @export_rt_roi_observations_module_field;
    case 'rt_dose'
        tagS = RT_dose_module_tags;
        export_function = @export_rt_dose_module_field;
    case 'rt_dvh'
        tagS = RT_DVH_module_tags;
        export_function = @export_RT_DVH_module_field;
    case 'multi_frame'
        tagS = multi_frame_module_tags;
        export_function = @export_multi_frame_module_field;
    case 'gsps'
        tagS = gsps_module_tags;
        export_function = @export_gsps_module_field;
%     case 'graphic_layer'
%         tagS = graphic_layer_module_tags;
%         export_function = @export_graphic_layer_module_field;
    otherwise
        error('Unrecognized or unsupported module export requested.')
end

%Prepare a attr template to be used by the export function for blank
%elements.
template = build_module_template(moduleName);

%Convert top level tags to decimal form.
%allTags = hex2dec({tagS.tag});
allTags = [tagS.tagdec];

%Iterate over all tags, calling the appropriate export function.
for i=1:length(allTags)

   %Prepare argments struct to pass to export function.
   args = [];
   args.tag         = allTags(i);
   args.data        = varargin;
   args.template    = template;

   %Get the Attributes object with data properly filled.
   %In this line "export_function" is a function HANDLE not an actual
   %function.  Need to check backwards compatibility with previous ML vers.

   el = export_function(args);

   if (~isempty(el))
       attr.addAll(el);
   else
       %Data in planC was insufficent to construct this element.
   end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function Image = AddTextToImageWithBorder(Image,String,Position,Color,Font,FontSize,BorderWidth,BorderColor)

% Image = AddTextToImage(Image,String,Position,Color,Font,FontSize,BorderWidth,BorderColor)
%
%   Works as AddTextToImage, except also allows for a border of thickness
%   BorderWidth pixels and color given by BorderColor. BorderColor should
%   be a 1- or 3-element vector giving Y or RGB in the range 0 to 1
%   regardless of the image class.
%
%	If no border parameters are provided, default is 1 pixel black.
%
% Daniel Warren
% Particle Therapy Cancer Research Institute
% University of Oxford

if ~exist('BorderWidth','var')
	BorderWidth = 1;
end
if ~exist('BorderColor','var')
	BorderColor = 0;
end

if BorderWidth == 0
	Image = AddTextToImage(Image,String,Position,Color,Font,FontSize);
	return;
end

% uint8 images go from 0 to 255, whereas double ones go from 0 to 1
if isa(Image, 'uint8')
    ScaleFactor = 255;
else
    ScaleFactor = 1;
end

% monochrome images need monochrome text, colour images need colour text
if ndims(Image) == 2 %#ok<ISMAT>
    BorderColor = mean(BorderColor(:));
end
if ndims(Image) == 3 && numel(BorderColor) == 1
    BorderColor = [BorderColor BorderColor BorderColor];
end

Mask = AddTextToImage(false(size(Image(:,:,1))),String,Position,1,Font,FontSize);
ConvElement = double(CircleMask(1+2*BorderWidth,1+2*BorderWidth,BorderWidth+0.5,BorderWidth+0.5,0.5+BorderWidth));
Outline = xor(conv2(double(Mask),ConvElement,'same'),Mask);
Image = AddTextToImage(Image,String,Position,Color,Font,FontSize);

for i = 1:size(Image,3)
tmp = Image(:,:,i);
tmp(Outline) = ScaleFactor*BorderColor(i);
Image(:,:,i) = tmp;
end

end

function out = CircleMask(w,h,cx,cy,r)
% function out = CircleMask(w,h,cx,cy,r)
%     Outputs a 1D bitmap of a solid circle of radius r pixels
%     centered at pixel (cx,cy) within an matrix of size w x h
%     pixels.

[x y] = meshgrid(1:w,1:h);
out = ((x-cx-0.5).^2+(y-cy-0.5).^2) <= r^2;

end

author: zouchuhang

repo: LayoutNet