text:
%XYZLABEL Label X, Y and Z axes
%
% XYZLABEL() label the x-, y- and z-axes with 'X', 'Y', and 'Z'
% respectiveley.
%
% XYZLABEL(FMT) as above but pass in a format string where %s is substituted
% for the axis label, eg.
%
%          xyzlabel('This is the %s axis')
%
% See also xlabel, ylabel, zlabel, sprintf.

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

author: petercorke

repo: spatialmath-matlab


text:
function shim = cvx_glpk( shim )

% CVX_SOLVER_SHIM	GLPK interface for CVX.
%   This procedure returns a 'shim': a structure containing the necessary
%   information CVX needs to use this solver in its modeling framework.

if ~isempty( shim.solve ),
    return
end
if isempty( shim.name ),
    fname = 'glpk.m';
    ps = pathsep;
    shim.name = 'GLPK';
    shim.dualize = true;
    flen = length(fname);
    fpaths = which( fname, '-all' );
    if ~iscell(fpaths),
      fpaths = { fpaths };
    end
    old_dir = pwd;
    oshim = shim;
    shim = [];
    for k = 1 : length(fpaths),
        fpath = fpaths{k};
        if ~exist( fpath, 'file' ) || any( strcmp( fpath, fpaths(1:k-1) ) ),
            continue
        end
        new_dir = fpath(1:end-flen-1);
        cd( new_dir );
        tshim = oshim;
        tshim.fullpath = fpath;
        tshim.version = 'unknown';
        tshim.location = new_dir;
        if isempty( tshim.error ),
            tshim.check = @check;
            tshim.solve = @solve;
            tshim.eargs = {};
            if k ~= 1,
                tshim.path = [ new_dir, ps ];
            end
        end
        shim = [ shim, tshim ]; %#ok
    end
    cd( old_dir );
    if isempty( shim ),
        shim = oshim;
        shim.error = 'Could not find a GLPK installation.';
    end
else
    shim.check = @check;
    shim.solve = @solve;
end

function found_bad = check( nonls ) %#ok
found_bad = false;

function [ x, status, tol, iters, y, z ] = solve( At, b, c, nonls, quiet, prec, settings )

n  = length( c );
m  = length( b );
lb = -Inf(n,1);
ub = +Inf(n,1);
vtype = 'C';
vtype = vtype(ones(n,1));
ctype = 'S';
ctype = ctype(ones(m,1));
rr = zeros(0,1);
cc = rr;
vv = rr;
zinv = rr;
is_ip = false;
for k = 1 : length( nonls ),
    temp = nonls( k ).indices;
    nn = size( temp, 1 );
    nv = size( temp, 2 );
    tt = nonls( k ).type;
    if strncmp( tt, 'i_', 2 ),
      is_ip = true;
      vartype(temp) = 'I';
      if strcmp(tt,'i_binary'),
        lb(temp) = 0;
        ub(temp) = 1;
      end
    elseif nn == 1 || isequal( tt, 'nonnegative' ),
        lb(temp) = 0;
    elseif isequal( tt, 'lorentz' ),
        if nn == 2,
            rr2  = [ temp ; temp ];
            cc2  = reshape( floor( 1 : 0.5 : 2 * nv + 0.5 ), 4, nv );
            vv2  = [1;1;-1;1]; vv = vv(:,ones(1,nv));
            rr   = [ rr ; rr(:) ];
            cc   = [ cc ; cc(:) ];
            vv   = [ vv ; vv(:) ];
            zinv = [ zinv ; temp(:) ];
        else
            error('GLPK does not support nonlinear constraints.' );
        end
    else
      error('GLPK does not support nonlinear constraints.' );
    end
end
if ~isempty(rr),
  znorm = [1:n]';
  znorm(zinv) = [];
  rr = [ rr ; znorm ];
  cc = [ cc ; znorm ];
  vv = [ vv ; ones(size(znorm)) ];
  reord = sparse( rr, cc, vv, n, n );
  At = reord' * At;
  c  = reord' * c;
end
if quiet,
  param.msglev = 0;
else
  param.msglev = 2;
end
param.scale = 128;
param.tolbnd = prec(1);
param.toldj = prec(1);
param.tolobj = prec(1);
[ xx, fmin, errnum, extra ] = cvx_run_solver( @glpk, c, At', b, lb, ub, ctype, vtype, 1, param, 'xx', 'fmin', 'errnum', 'extra', settings, 9 );
tol   = [];
iters = [];
x = full( xx );
y = full( extra.lambda );
z = full( extra.redcosts );
if ~isempty( rr ),
  x = reord * x;
  z = reord * z;
  z(zinv) = z(zinv) * 0.5;
end
status = 'Failed';
switch errnum,
case 0,
  switch extra.status,
  case 2,
    if is_ip,
      status = 'Suboptimal';
    elseif errnumn == 0,
      status = 'Solved';
    else
      status = 'Inaccurate/Solved';
    end
  case {3,4},
    status = 'Infeasible';
  case 5,
    status = 'Solved';
  case 6,
    status = 'Unbounded';
  end
case 10,
  status = 'Infeasible';
case 11,
  status = 'Unbounded';
case {5,17,15,19}
  status = 'Failed';
case {6,7,8,9,13,14},
  switch extra.status,
  case {2,5}
    if is_ip,
      status = 'Suboptimal';
    else
      status = 'Inaccurate/Solved';
    end
  case { 3,4 },
    status = 'Inaccurate/Infeasible';
  case 6,
    status = 'Inaccurate/Unbounded';
  end
end
if strcmp(status,'Failed'),
  tol = Inf;
elseif strncmp(status,'Inaccurate/',11),
  tol = prec(3);
else
  tol = prec(2);
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [a,f,B,W]=v_ldatrace(b,w,n,c)
%V_LDATRACE Calculates an LDA transform to maximize trace discriminant [a,f,B,W]=(b,w,n,c)
% If a feature vector X can come from one of several class and W and B are respectively
% the within-class and between-class covariance matrices, then the generalized Fisher discriminant
% F=trace(W\B) is a measure of how well the feature vector discriminates between the classes.
% If we choose a rectangular (tall, skinny) transformation matrix, we can define a smaller
% feature vector Y=A'*X. The aim of this routine is to choose A to maximize the Fisher
% discriminant. We assume that W is positive definite and B is positive semi-definite.
% The input argument C allows the uset to pre-specify some of the columns of A.
%
% Inputs:
%     w[m,m] = within class covariance matrix of x
%     b[m,m] = between class covariance matrix of x [default = I]
%     n is the number of columns in output matrix A [default = M]
%     c[m,r] specifies the first few columns of A to be predefined values [default = null)
%
% Outputs:
%     a[m,n] is the transformation matrix: y=a'*x
%     f[n,1] gives the incremental gain in f value for successive columns of A
%     B(n,n) gives the between-class covariance matrix of y
%     W[n,n] gives the within-class covariance matrix of y

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: v_ldatrace.m 10865 2018-09-21 17:22:45Z 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

m=size(b,1);    % dimension of data vectors
if nargin<4
    r=0;
    if nargin<3
        n=m;
        if nargin<2
            w=eye(m);
        end
    end
else
    r=size(c,2);    % number of columns that are pre-specified
end
if r
    if n>r          % need to find additional vectors
        g=chol(w);
        v=g\null(c'*g');
        [p,l,q]=svd(v'*b*v);
        a(:,r+1:n)=v*p(:,1:n-r);
        a(:,1:r)=c;
    else
        a=c;        % no new vectors to find
    end
    if nargout>1
        ari=a/triu(qr(chol(a'*w*a))); % matrix a must be of full rank
        f=diag(ari'*b*ari);
    end
else
    [g,d]=eig(b,w,'qz');
    [ds,is]=sort(-diag(d));
    a=g(:,is(1:n));
    if nargout>1
        f=-ds(1:n);
    end
end
if nargout > 2
    B=a'*b*a;
    if nargout > 3
        W=a'*w*a;
    end

end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function nlp = IdentityMapConstraint(obj, nlp, src, tar, bounds, varargin)
    % Impose system specific constraints as NLP constraints for the
    % trajectory optimization problem
    %
    % Parameters:
    % nlp: an trajectory optimization NLP object of the system
    % @type TrajectoryOptimization
    % src: the NLP of the source domain system
    % @type TrajectoryOptimization
    % tar: the NLP of the target domain system
    % @type TrajectoryOptimization
    % bounds: a struct data contains the boundary values @type struct
    % varargin: extra arguments @type varargin

    % by default, identity map
    x = obj.States.x;
    xn = obj.States.xn;
    x_map = SymFunction(['xDiscreteMap' obj.Name],x-xn,{x,xn});

    addNodeConstraint(nlp, x_map, {'x','xn'}, 'first', 0, 0, 'Linear');

    if strcmp(obj.Plant.Type,'SecondOrder')
        dx = obj.States.dx;
        dxn = obj.States.dxn;
        dx_map = SymFunction(['dxDiscreteMap' obj.Name],dx-dxn,{dx,dxn});

        addNodeConstraint(nlp, dx_map, {'dx','dxn'}, 'first', 0, 0, 'Linear');
    end

end

author: ayonga

repo: frost-dev

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


text:
function D = spm_eeg_inv_forward(varargin)
% Compute M/EEG leadfield
% FORMAT D = spm_eeg_inv_forward(D,val)
%
% D                - input struct
% (optional) fields of S:
% D                - filename of EEG/MEG mat-file
%
% Output:
% D                - EEG/MEG struct with filenames of Gain matrices)
%__________________________________________________________________________
% Copyright (C) 2008-2018 Wellcome Trust Centre for Neuroimaging

% Jeremie Mattout & Christophe Phillips
% $Id: spm_eeg_inv_forward.m 7702 2019-11-22 11:32:26Z guillaume $

SVNrev = '$Rev: 7702 $';

%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename, SVNrev);

%-Initialisation
%--------------------------------------------------------------------------
[D, val] = spm_eeg_inv_check(varargin{:});

if numel(D.inv{val}.datareg) ~= numel(D.inv{val}.forward)
    error('Separate coregistration is required for every modality.');
end

Fgraph = spm_figure('FindWin','Graphics');
spm_figure('Clear',Fgraph);
spm('Pointer', 'Watch');
if isempty(Fgraph) || spm('CmdLine'), graph = 'no'; else, graph = 'yes'; end

for i = 1:numel(D.inv{val}.forward)
    M    = D.inv{val}.datareg(i).fromMNI*D.inv{val}.mesh.Affine;

    M    = diag([1e-3 1e-3 1e-3 1])*M; % convert to m

    mesh = spm_eeg_inv_transform_mesh(M, D.inv{val}.mesh);

    mesh_correction = [];

    sens = D.inv{val}.datareg(i).sensors;

    if isequal(D.inv{val}.datareg(i).modality, 'MEG')
        sens = ft_datatype_sens(sens, 'amplitude', 'T', 'distance', 'm');
    else
        sens = ft_datatype_sens(sens, 'amplitude', 'V', 'distance', 'm');
    end

    switch D.inv{val}.forward(i).voltype
        case 'EEG interpolated'
            vol = D.inv{val}.forward(i).vol;
            modality = 'EEG';
        case '3-Shell Sphere'
            cfg              = [];
            cfg.feedback     = graph;
            cfg.siunits      = 'yes';
            cfg.showcallinfo = 'no';

            headshape(1) = export(gifti(mesh.tess_scalp),  'ft');
            headshape(2) = export(gifti(mesh.tess_oskull), 'ft');
            headshape(3) = export(gifti(mesh.tess_iskull), 'ft');

            % determine the convex hull of the brain, to determine the support points
            pnt  = mesh.tess_ctx.vert;
            tric = convhulln(pnt);
            sel  = unique(tric(:));

            % create a triangulation for only the support points
            headshape(4).pnt = pnt(sel, :);
            headshape(4).tri = convhulln(pnt(sel, :));

            cfg.method = 'concentricspheres';

            vol  = ft_prepare_headmodel(cfg, headshape);

            cfg = [];
            cfg.headmodel = vol;
            cfg.grid.pos  = mesh.tess_ctx.vert;
            cfg.spherify  = 'yes';
            gridsphere    = ft_prepare_sourcemodel(cfg);

            mesh_correction    = rmfield(cfg, {'headmodel', 'grid'});

            mesh.tess_ctx.vert = gridsphere.pos;
            modality = 'EEG';

        case 'EEG BEM'
            volfile = spm_file(mesh.sMRI, 'suffix','_EEG_BEM', 'ext','mat');
            vol = [];

            if exist(volfile, 'file')
                vol = ft_read_headmodel(volfile);
                if ~isfield(vol, 'unit') || ~isequal(vol.unit, 'm')
                    vol = [];
                end
            end

            if isempty(vol)

                vol.cond   = [0.3300 0.0041 0.3300];
                vol.source = 1; % index of source compartment
                vol.skin   = 3; % index of skin surface
                % brain
                vol.bnd(1) = export(gifti(mesh.tess_iskull), 'ft');
                % skull
                vol.bnd(2) = export(gifti(mesh.tess_oskull), 'ft');
                % skin
                vol.bnd(3) = export(gifti(mesh.tess_scalp),  'ft');

                % create the BEM system matrix
                cfg        = [];
                cfg.method = 'bemcp';
                cfg.showcallinfo = 'no';
                cfg.siunits      = 'yes';
                vol = ft_prepare_headmodel(cfg, vol);

                spm_progress_bar('Set', 1);

                save(volfile, 'vol', spm_get_defaults('mat.format'));

                spm_progress_bar('Clear');
                spm('Pointer', 'Arrow');
            end

            cfg = [];
            cfg.headmodel = vol;
            cfg.grid.pos = mesh.tess_ctx.vert;
            cfg.moveinward = 6e-3; %move to empirically determined BEM safe zone
            gridcorrect = ft_prepare_sourcemodel(cfg);

            mesh_correction    = rmfield(cfg, {'headmodel', 'grid'});

            mesh.tess_ctx.vert = gridcorrect.pos;

            vol = volfile;
            modality = 'EEG';

        case 'OpenMEEG BEM'
            vol        = [];
            vol.cond   = [0.3300 0.0041 0.3300];
            vol.source = 1; % index of source compartment
            vol.skin   = 3; % index of skin surface
            % brain
            vol.bnd(1) = export(gifti(mesh.tess_iskull), 'ft');
            % skull
            vol.bnd(2) = export(gifti(mesh.tess_oskull), 'ft');
            % skin
            vol.bnd(3) = export(gifti(mesh.tess_scalp),  'ft');

            cfg                         = [];
            cfg.method                 = 'openmeeg';
            cfg.siunits                = 'yes';
            cfg.showcallinfo           = 'no';
            vol                        = ft_prepare_headmodel(cfg, vol);

            cfg                        = [];
            cfg.vol                    = vol;
            cfg.grid.pos               = mesh.tess_ctx.vert;
            cfg.moveinward             = 6e-3; % smaller shift might suffice for OpenMEEG
            gridcorrect                = ft_prepare_sourcemodel(cfg);

            mesh_correction            = rmfield(cfg, {'vol', 'grid'});

            mesh.tess_ctx.vert         = gridcorrect.pos;

            modality = 'EEG';
        case 'Single Sphere'
            cfg                        = [];
            cfg.feedback               = 'yes';
            cfg.showcallinfo           = 'no';
            cfg.grad                   = D.inv{val}.datareg(i).sensors;
            cfg.method                 = 'singlesphere';
            cfg.siunits                = 'yes';

            headshape                  = export(gifti(mesh.tess_scalp), 'ft');

            vol                        = ft_prepare_headmodel(cfg, headshape);
            modality                   = 'MEG';
        case 'MEG Local Spheres'
            cfg                        = [];
            cfg.feedback               = 'yes';
            cfg.showcallinfo           = 'no';
            cfg.grad                   = sens;
            cfg.method                 = 'localspheres';
            cfg.siunits                = 'yes';

            headshape                  = export(gifti(mesh.tess_scalp), 'ft');
            vol                        = ft_prepare_headmodel(cfg, headshape);
            modality                   = 'MEG';
        case  'Single Shell'
            cfg                        = [];
            cfg.feedback               = 'yes';
            cfg.showcallinfo           = 'no';
            cfg.grad                   = sens;
            cfg.method                 = 'singleshell';
            cfg.siunits                = 'yes';

            headshape                  = export(gifti(mesh.tess_iskull), 'ft');

            vol                        = ft_prepare_headmodel(cfg, headshape);
            modality                   = 'MEG';

        otherwise
            error('Unsupported volume model type.');
    end

    D.inv{val}.forward(i).vol             = vol;
    D.inv{val}.forward(i).mesh            = mesh.tess_ctx;
    D.inv{val}.forward(i).mesh_correction = mesh_correction;
    D.inv{val}.forward(i).modality        = modality;
    D.inv{val}.forward(i).siunits         = 1;

    D.inv{val}.forward(i).sensors  = sens;

    D.inv{val}.forward(i).toMNI    = D.inv{val}.datareg(i).toMNI*diag([1e3 1e3 1e3 1]);
    D.inv{val}.forward(i).fromMNI  = diag([1e-3 1e-3 1e-3 1])*D.inv{val}.datareg(i).fromMNI;

    spm_figure('Clear',Fgraph);
end

% This is to force recomputing the lead fields
try, D.inv{val} = rmfield(D.inv{val}, 'gainmat'); end

fprintf('%-40s: %30s\n','Completed',spm('time'));                       %-#
spm('Pointer', 'Arrow');

author: spm


text:
function plotData(x, y)
    %% PLOTDATA Plots the data points x and y into a new figure
    %   PLOTDATA(x,y) plots the data points and gives the figure axes labels of
    %   population and profit.

    % ====================== YOUR CODE HERE ======================
    % Instructions: Plot the training data into a figure using the
    %               "figure" and "plot" commands. Set the axes labels using
    %               the "xlabel" and "ylabel" commands. Assume the
    %               population and revenue data have been passed in
    %               as the x and y arguments of this function.
    %
    % Hint: You can use the 'rx' option with plot to have the markers
    %       appear as red crosses. Furthermore, you can make the
    %       markers larger by using plot(..., 'rx', 'MarkerSize', 10);

    figure; % open a new figure window
    plot(x, y, 'rx', 'MarkerSize', 10);
    ylabel('Profit in $10,000s');
    xlabel('Population of City in 10,000s');
end

author: worldveil

repo: coursera-ml


text:
function [W,BC,DV,Q,r,B] = voxelize(V,F,side,varargin)
  % VOXELIZE Given a mesh compute a voxelization of that mesh on a regular grid
  % fit to the bounding box.
  %
  % W = voxelize(V,F,side)
  % [W,BC,DV,Q] = voxelize(V,F,side,'ParameterName',ParameterValue, ...)
  %
  % Inputs:
  %   V  #V by 3 list of vertex positions
  %   F  #F by 3 list of triangle indices into V, expects counterclockwise
  %     orientation to produce outward pointing normal
  %   side  either the number of cells in the x direction or a triple
  %     containing the number of cells in each dimension
  %   Optional:
  %     'Pad'  followed by number of padding to add to each side. 1 would mean
  %       add one extra cell on the top, bottom, left, right, front, and back
  %       of bounding box lattice {0}.
  %     'Boundary'  followed by wether to consider any cell intersecting the
  %       mesh as inside {true}.
  %     'Interior'  followed by whether to consider any cell full inside the
  %        mesh as inside {true}
  %     'Closed'  followed by whether to assume mesh is closed when determining
  %       interior (avoid winding number computation, use flood filling)
  %     'Centers'  followed by prod(side) centers (e.g., already produced by
  %       call to previous call to voxel_grid.m)
  % Outputs:
  %   W  side(1) by side(2) by side(3) matrix with W(i,j,k) ~= if location
  %     cell centered at BC(i,j,k) overlaps with the volume of (V,F)
  %   BC  prod(side) by dim matrix of cell barycenters
  %   DV  side(1)+1 by side(2)+1 by side(3)+1 matrix of cell corner locations
  %   Q  #Q by 4 list of quads indexing DV
  %   r  dim-long list of voxel widths in each direction
  %
  % Known issues: the ouput surface mesh will contain non-manifold edges and
  % vertices at voxels that meet at such edges or vertices.
  %
  % Example:
  %   [W,BC,DV,Q] = voxelize(V,F,50);
  %   % Get triangles from quads
  %   DF = [Q(:,[1 2 3]);Q(:,[1 3 4])];
  %   % Remove unreferenced corners
  %   [SV,IM] = remove_unreferenced(DV,DF);
  %   % re-index
  %   SF = IM(DF);
  %   SQ = IM(Q);
  %
  %   [W,BC] = voxelize(V,F,50,'Pad',1);
  %   % Use matlab's isosurface to extract surface as triangle mesh
  %   BC = reshape(BC,[size(W) 3]);
  %   surf = isosurface(BC(:,:,:,1),BC(:,:,:,2),BC(:,:,:,3),W,0.5);
  %   DV = bsxfun(@plus,bsxfun(@times,bsxfun(@rdivide,surf.vertices-1, ...
  %     [size(W,2) size(W,1) size(W,3)]-1),max(BC)-min(BC)),min(BC));
  %   DF = surf.faces;
  %
  %
  % See also: bwboundaries, isosurface

  with_boundary = true;
  with_interior = true;
  pad_count = 0;
  closed = [];
  BC = [];
  % default values
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Boundary','Interior','Pad','Centers','Closed'}, ...
    {'with_boundary','with_interior','pad_count','BC','closed'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end
  if isempty(closed)
    closed = isempty(boundary_faces(F));
  end

  assert(any(size(side) == 1),'side should be a vector');
  % Make sure we have a row vector
  side = side(:)';

  dim = size(V,2);

  if isempty(BC)
    [BC,side,r] = voxel_grid(V,side,'Pad',pad_count);
  else
    assert(prod(side) == size(BC,1),'side dims must match BC');
    r = (max(BC)-min(BC))./(side-1);
  end
  NV = min(BC);
  XV = max(BC);

  switch dim
  case 2
    W  = zeros([side(2) side(1)]);
    [DV,~,QT,QL] = voxel_surface(W,'Centers',BC);

    if with_interior && ~closed
      WDV = reshape(winding_number(V,F,DV),[side(2) side(1)]+1);
      WDV = abs(WDV) >= 0.5;
      W = ( ...
        WDV(1:end-1,1:end-1) | ...
        WDV(1:end-1,  2:end) | ...
        WDV(  2:end,1:end-1) | ...
        WDV(  2:end,  2:end));
    end
    Q = [QT;QL];
    if with_boundary
      for ei = 1:size(F,1)

        sqrD = segment_segment_squared_distance( ...
          V(F(ei,1),:), V(F(ei,2),:), ...
          DV(Q(:,1),:), DV(Q(:,2),:));
        IF = mod(find(sqrD<eps)-1,size(Q,1))+1;
        iQT = IF(IF<=size(QT,1));
        iQL = IF(IF>size(QT,1))-size(QT,1);

        [II,JJ] = ind2sub([side(2) side(1)]+1,max(QT(iQT,:),[],2));
        W(sub2ind(side([2 1]),II-1,JJ-1)) = 1;
        W(sub2ind(side([2 1]),II  ,JJ-1)) = 1;

        [II,JJ] = ind2sub([side(2) side(1)]+1,max(QL(iQL,:),[],2));
        W(sub2ind(side([2 1]),II-1,JJ-1)) = 1;
        W(sub2ind(side([2 1]),II-1,JJ  )) = 1;

        %clf;
        %hold on;
        %surf(reshape(BC(:,1),size(W)),reshape(BC(:,2),size(W)),reshape(BC(:,1)*0,size(W)),'CData',W*1,fphong,'EdgeColor','none');
        %plot_edges(V,F(ei,:),'y','LineWidth',8);
        %plot_edges(V,F,'g','LineWidth',4);
        %text(BC(:,1),BC(:,2),num2str((1:size(BC,1))'),'BackgroundCOlor',[0.8 0.8 0.8]);
        %text(DV(:,1),DV(:,2),num2str((1:size(DV,1))'),'BackgroundCOlor',[0.8 0.1 0.1]);
        %plot_edges(DV,QT,'--r','LineWidth',1);
        %plot_edges(DV,QL,'--b','LineWidth',1);
        %plot_edges(DV,QT(iQT,:),'r','LineWidth',3);
        %plot_edges(DV,QL(iQL,:),'b','LineWidth',3);
        %hold off;
        %axis equal;
        %colormap([zeros(20,3);1 1 1]);
        %drawnow;

      end
    end
    if with_interior && closed
      W = imfill(W,'holes');
    end
    % Pad winding number with 0s and take differences in all 6 directions
    Wp = padarray(W,[1 1],0);
    Dy = diff(Wp,1,1);
    Dx = diff(Wp,1,2);
    Q = [ ...
      QL(Dx(2:end-1,1:end,2:end-1)>0.5,:); ...
      QT(Dy(1:end,2:end-1,2:end-1)>0.5,:); ...
      fliplr(QL(Dx(2:end-1,1:end,2:end-1)<-0.5,:)); ...
      fliplr(QT(Dy(1:end,2:end-1,2:end-1)<-0.5,:)); ...
      ];
  case 3
    W = zeros([side(2) side(1) side(3)]);
    [DV,~,QT,QF,QL] = voxel_surface(W,'Centers',BC);

    if with_interior && ~closed
      % Winding number is a heavy handed way of determining inside/outside for a
      % simple closed polyhedron.  If the boundary has already been detected then
      % this should be floodfilling instead.

      %W = winding_number(V,F,BC);
      WDV = reshape(winding_number(V,F,DV,'Fast',true),[side(2) side(1) side(3)]+1);
      WDV = abs(WDV) >= 0.5;
      W = ( ...
        WDV(1:end-1,1:end-1,1:end-1) | ...
        WDV(1:end-1,1:end-1,2:end) | ...
        WDV(1:end-1,  2:end,1:end-1) | ...
        WDV(1:end-1,  2:end,2:end) | ...
        WDV(  2:end,1:end-1,1:end-1) | ...
        WDV(  2:end,1:end-1,2:end) | ...
        WDV(  2:end,  2:end,1:end-1) | ...
        WDV(  2:end,  2:end,2:end));
    end

    Q = [QT;QF;QL];

    if with_boundary
      % TODO: Could at least ignore already-inside cells.
      DF = [Q(:,[3 2 1]);Q(:,[4 3 1])];
      IF = intersect_other(V,F,DV,DF);
      JF = IF(:,1);
      IF = mod(IF(:,2)-1,size(Q,1))+1;
      B = zeros(side(2),side(1),side(3));
      % Force winding number to mark voxels intersecting boundary as inside
      % (This should just be replaced with rasterizing the surface)
      iQ = cell(3,1);
      iQ{1} = IF(IF<=size(QT,1));
      iQ{3} = IF(IF>size(QT,1) & IF<=size(QT,1)+size(QF,1))-size(QT,1);
      iQ{2} = IF(IF>size(QT,1)+size(QF,1))-size(QT,1)-size(QF,1);
      jQ{1} = JF(IF<=size(QT,1));
      jQ{3} = JF(IF>size(QT,1) & IF<=size(QT,1)+size(QF,1));
      jQ{2} = JF(IF>size(QT,1)+size(QF,1));

      side123 = side([2 1 3]);
      for d = 1:3
        sz = side123;
        sz(d) = sz(d)+1;
        II = cell(3,1);
        [II{1},II{2},II{3}] = ind2sub(sz,iQ{d});
        for pass = 0:1
          II{d} = II{d}-pass;
          keep = II{d}<=side123(d);
          B(sub2ind(side123,II{1}(keep),II{2}(keep),II{3}(keep))) = jQ{d}(keep);
        end
      end

      W = W | B;
    end

    if with_interior && closed
      W = imfill(W,'holes');
    end

    %trisurf(Q(IF,:),DV(:,1),DV(:,2),DV(:,3),'FaceColor','r');
    %hold on;
    %tsurf(F,V,'FaceColor','g','FaceAlpha',0.2,'EdgeAlpha',0.2);
    %hold off;

    %% From cells to faces:
    %[II,JJ,KK] = ind2sub([side(2) side(1) side(3)],find(abs(W)>0.5));
    %Q = [ ...
    %  QT(sub2ind([side(2)+1 side(1)   side(3)  ],II+1,JJ,KK),:); ...
    %  QT(sub2ind([side(2)+1 side(1)   side(3)  ],II+0,JJ,KK),:); ...
    %  QF(sub2ind([side(2)   side(1)   side(3)+1],II,JJ,KK+1),:); ...
    %  QF(sub2ind([side(2)   side(1)   side(3)+1],II,JJ,KK+0),:); ...
    %  QL(sub2ind([side(2)   side(1)+1 side(3)  ],II,JJ+1,KK),:); ...
    %  QL(sub2ind([side(2)   side(1)+1 side(3)  ],II,JJ+0,KK),:); ...
    %  ];
    %trisurf(Q,DV(:,1),DV(:,2),DV(:,3),'FaceAlpha',0.5,'EdgeAlpha',0.5);
    %hold on;
    %tsurf(F,V,'FaceColor','g','FaceAlpha',0.2,'EdgeAlpha',0.2);
    %hold off;
    %axis equal
    %
    %error

    % Pad winding number with 0s and take differences in all 6 directions
    Wp = padarray(W,[1 1 1],0);
    Dy = diff(Wp,1,1);
    Dx = diff(Wp,1,2);
    Dz = diff(Wp,1,3);
    Q = [ ...
      QF(Dz(2:end-1,2:end-1,1:end)>0.5,:); ...
      QL(Dx(2:end-1,1:end,2:end-1)>0.5,:); ...
      QT(Dy(1:end,2:end-1,2:end-1)>0.5,:); ...
      fliplr(QF(Dz(2:end-1,2:end-1,1:end)<-0.5,:)); ...
      fliplr(QL(Dx(2:end-1,1:end,2:end-1)<-0.5,:)); ...
      fliplr(QT(Dy(1:end,2:end-1,2:end-1)<-0.5,:)); ...
      ];
  end
end

author: alecjacobson

repo: gptoolbox


text:
function [tKuka,flag]=connectToKuka( jPos)
% Copy right: Mohammad SAFEEA
% 16th-Oct-2017

cDir = pwd;
cDir=getTheKSTDirectory(cDir);
addpath(cDir);

warning('off');
ip='172.31.1.147'; % The IP of the controller
% start a connection with the server
tKuka=net_establishConnection( ip );

if ~exist('tKuka','var') || isempty(tKuka)
  warning('Connection could not be establised, script aborted');
  flag=false;
  return;
end
 % Move robot to home position
      setBlueOff(tKuka); % turn Off blue light

      relVel=0.15;
      movePTPJointSpace( tKuka , jPos, relVel); % move to initial configuration

       flag=true;
end

repo: KST-Kuka-Sunrise-Toolbox


text:
% pop_chansel() - pop up a graphic interface to select channels
%
% Usage:
%   >> [chanlist] = pop_chansel(chanstruct); % a window pops up
%   >> [chanlist strchannames cellchannames] = ...
%                        pop_chansel(chanstruct, 'key', 'val', ...);
%
% Inputs:
%   chanstruct     - channel structure. See readlocs()
%
% Optional input:
%   'withindex'      - ['on'|'off'] add index to each entry. May also a be
%                      an array of indices
%   'select'         - selection of channel. Can take as input all the
%                      outputs of this function.
%   'selectionmode' - selection mode 'multiple' or 'single'. See listdlg2().
%
% Output:
%   chanlist      - indices of selected channels
%   strchannames  - names of selected channel names in a concatenated string
%                   (channel names are separated by space characters)
%   cellchannames - names of selected channel names in a cell array
%
% Author: Arnaud Delorme, CNL / Salk Institute, 3 March 2003

% Copyright (C) 3 March 2003 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% 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

function [chanlist,chanliststr, allchanstr] = pop_chansel(chans, varargin);

    if nargin < 1
        help pop_chansel;
        return;
    end;
    if isempty(chans), disp('Empty input'); return; end;
    if isnumeric(chans),
        for c = 1:length(chans)
            newchans{c} = num2str(chans(c));
        end;
        chans = newchans;
    end;
    chanlist    = [];
    chanliststr = {};
    allchanstr  = '';

    g = finputcheck(varargin, { 'withindex'     {  'integer' 'string' } { [] {'on' 'off'} }   'off';
                                'select'        { 'cell' 'string' 'integer' } [] [];
                                'selectionmode' 'string' { 'single' 'multiple' } 'multiple'});
    if isstr(g), error(g); end;
    if ~isstr(g.withindex), chan_indices = g.withindex; g.withindex = 'on';
    else                    chan_indices = 1:length(chans);
    end;

    % convert selection to integer
    % ----------------------------
    if isstr(g.select) & ~isempty(g.select)
        g.select = parsetxt(g.select);
    end;
    if iscell(g.select) & ~isempty(g.select)
        if isstr(g.select{1})
            tmplower = lower( chans );
            for index = 1:length(g.select)
                matchind = strmatch(lower(g.select{index}), tmplower, 'exact');
                if ~isempty(matchind), g.select{index} = matchind;
                else error( [ 'Cannot find ''' g.select{index} '''' ] );
                end;
            end;
        end;
        g.select = [ g.select{:} ];
    end;
    if ~isnumeric( g.select ), g.select = []; end;

    % add index to channel name
    % -------------------------
	tmpstr = {chans};
    if isnumeric(chans{1})
        tmpstr = [ chans{:} ];
        tmpfieldnames = cell(1, length(tmpstr));
        for index=1:length(tmpstr),
            if strcmpi(g.withindex, 'on')
                tmpfieldnames{index} = [ num2str(chan_indices(index)) '  -  ' num2str(tmpstr(index)) ];
            else
                tmpfieldnames{index} = num2str(tmpstr(index));
            end;
        end;
    else
        tmpfieldnames = chans;
        if strcmpi(g.withindex, 'on')
            for index=1:length(tmpfieldnames),
                tmpfieldnames{index} = [ num2str(chan_indices(index)) '  -  ' tmpfieldnames{index} ];
            end;
        end;
    end;
    [chanlist,tmp,chanliststr] = listdlg2('PromptString',strvcat('(use shift|Ctrl to', 'select several)'), ...
                'ListString', tmpfieldnames, 'initialvalue', g.select, 'selectionmode', g.selectionmode);
    if tmp == 0
        chanlist = [];
        chanliststr = '';
        return;
    else
        allchanstr = chans(chanlist);
    end;

    % test for spaces
    % ---------------
    spacepresent = 0;
    if ~isnumeric(chans{1})
        tmpstrs = [ allchanstr{:} ];
        if ~isempty( find(tmpstrs == ' ')) | ~isempty( find(tmpstrs == 9))
            spacepresent = 1;
        end;
    end;

    % get concatenated string (if index)
    % -----------------------
    if strcmpi(g.withindex, 'on') | spacepresent
        if isnumeric(chans{1})
            chanliststr = num2str(celltomat(allchanstr));
        else
            chanliststr = '';
            for index = 1:length(allchanstr)
                if spacepresent
                    chanliststr = [ chanliststr '''' allchanstr{index} ''' ' ];
                else
                    chanliststr = [ chanliststr allchanstr{index} ' ' ];
                end;
            end;
            chanliststr = chanliststr(1:end-1);
        end;
    end;

    return;

author: PatternRecognition

repo: OpenBMI


text:
function fem1d_bvp_quadratic_test01 ( )

%*****************************************************************************80
%
%% FEM1D_BVP_QUADRATIC_TEST01 carries out test case #1.
%
%  Discussion:
%
%    Use A1, C1, F1, EXACT1, EXACT_UX1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dianne O'Leary,
%    Scientific Computing with Case Studies,
%    SIAM, 2008,
%    ISBN13: 978-0-898716-66-5,
%    LC: QA401.O44.
%
  n = 11;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_BVP_QUADRATIC_TEST01\n' );
  fprintf ( 1, '  Solve -( A(x) U''(x) )'' + C(x) U(x) = F(x)\n' );
  fprintf ( 1, '  for 0 < x < 1, with U(0) = U(1) = 0.\n' );
  fprintf ( 1, '  A1(X)  = 1.0\n' );
  fprintf ( 1, '  C1(X)  = 0.0\n' );
  fprintf ( 1, '  F1(X)  = X * ( X + 3 ) * exp ( X )\n' );
  fprintf ( 1, '  U1(X)  = X * ( 1 - X ) * exp ( X )\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes = %d\n', n );
%
%  Geometry definitions.
%
  x_first = 0.0;
  x_last = 1.0;
  x = linspace ( x_first, x_last, n );
  x = x(:);

  u = fem1d_bvp_quadratic ( n, @a1, @c1, @f1, x );

  uexact = exact1 ( x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     I    X         U         Uexact    Error\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %4d  %8f  %8f  %8f  %8e\n', ...
      i, x(i), u(i), uexact(i), abs ( u(i) - uexact(i) ) );
  end

  e1 = l1_error ( n, x, u, @exact1 );
  e2 = l2_error_quadratic ( n, x, u, @exact1 );
  h1s = h1s_error_quadratic ( n, x, u, @exact_ux1 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  l1 norm of error  = %g\n', e1 );
  fprintf ( 1, '  L2 norm of error  = %g\n', e2 );
  fprintf ( 1, '  Seminorm of error = %g\n', h1s );

  return
end
function value = a1 ( x )

%*****************************************************************************80
%
%% A1 evaluates A function #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of A(X).
%
  value = 1.0;

  return
end
function value = c1 ( x )

%*****************************************************************************80
%
%% C1 evaluates C function #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 March 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of C(X).
%
  value = 0.0;

  return
end
function value = exact1 ( x )

%*****************************************************************************80
%
%% EXACT1 evaluates exact solution #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of U(X).
%
  value = x .* ( 1.0 - x ) .* exp ( x );

  return
end
function value = exact_ux1 ( x )

%*****************************************************************************80
%
%% EXACT_UX1 evaluates the derivative of exact solution #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of dUdX(X).
%
  value = ( 1.0 - x - x .* x ) .* exp ( x );

  return
end
function value = f1 ( x )

%*****************************************************************************80
%
%% F1 evaluates right hand side function #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 March 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of F(X).
%
  value = x .* ( x + 3.0 ) .* exp ( x );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% simple wrapper for softmax function
function [y] = softmax(x)

logZ = logsum(x, 2);
y = exp(bsxfun(@minus, x, logZ));

author: kyunghyuncho

repo: deepmat


text:
function msh = meshFV2msh(fv,mmPerVox,windowID, actor)
% Convert Matlab faces/vertices to vistasoft mesh
%
% Example:  See t_meshSimpleShapes, t_meshCurvature
%
% BW (c) VISTASOFT Team, Stanford, 2013

%% Check and initialize parameters
if notDefined('mmPerVox'), mmPerVox = [1 1 1]; end
if notDefined('windowID'), windowID = 1000; end
if notDefined('actor'), actor = 33; end

% Create a vista mesh structure
msh = meshCreate;
msh = meshSet(msh,'window id',windowID);
msh = meshSet(msh,'actor',actor);

%% Copy the Matlab vertices faces and vertices into msh
msh = meshSet(msh,'vertices',double(fv.vertices'));

%% Set the scene origin
vertices = meshGet(msh,'vertices');
msh = meshSet(msh,'origin',-mean(vertices,2)');

% Set the scale
msh = meshSet(msh,'mmPerVox',mmPerVox);

%% Triangles and normals

% Permute the triangle order, for some reason.  Otherwise we have apparent
% holes in the meshes - one side is transparent and the other side opaque.
% We haven't understood this yet.
msh = meshSet(msh,'triangles',double(fv.faces(:,[3 2 1])' - 1));

% Compute and set vertex normals
msh = meshSet(msh,'normals',patchnormals(fv)');

%% Calculate mean curvature

% try
%     % [Umin,Umax,Cmin,Cmax,Cmean,Cgauss,Normal]
%     % Cmax is no good
%     % Cmin is possible
%     % Cmean is good
%     % Cgauss is not so good
%     % n is useless
%     [~,~,~,~,Cmean,~,~] = compute_curvature(fv.vertices,fv.faces);
%     msh = meshSet(msh,'curvature',Cmean);
% catch err
% Different calculation, slower, but it seems to run.
fprintf('Using tricurv_v01\n');
curvatures = tricurv_v01(fv.faces,fv.vertices);
msh = meshSet(msh,'curvature',curvatures.km);
% end

%% Build  color values from curvature

% This should become a separate routine

% Negative curvature is brighter because we multiply by -1.
% We compress because we are mostly interested in the curvature near 0.
% s = 0.1;
% c = ((2*c > 0) - 1) * .25 * 128 + 127.5;  % Formula from meshColor

% If curvature is negative, make the color index bigger
curvatureSaturation = 0.1;
curvatureContrast = 0.45;
curvatureCutpoint = -0.15;
c = -1 * (msh.curvature(:)) ./ (abs(msh.curvature(:)) + curvatureSaturation);
c = (c > curvatureCutpoint) * curvatureContrast * 128 + 96;  % Formula from meshColor

% Place the color map in a 4xN
c = c*ones(1,3);
c(:,4) = 255;  % Make opaque (not working well)
c = round(c);
msh = meshSet(msh,'colors',c');

%% meshVisualize(msh)
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function cov = pc2covv(pc,R0,ncov)

%function cov = pc2covv(pc,R0,ncov)
%   Transforms partial correlations pc into a covariancefunction cov
%   ncov is the number of covariances
%   For a zero-mean process, the covariance function has been defined as
%   R(s)(.) = E{X[n+s]<X[n],.>},
%   or, in old-fashioned notation:
%   R(s) = E{X[n+s]X[n]'} where X[n] is a column vector.

%S. de Waele, March 2003.

if ~isstatv(pc), error('Partial correlations non-stationairy!'), end

s = kingsize(pc);
order = s(3)-1;
dim = s(1); I = eye(dim);

if nargin <  3, ncov = order; end

m = min(order,ncov);

par = zeros(dim,dim,m+1);
parb = zeros(dim,dim,m+1);
cov = zeros(dim,dim,ncov+1);

[rc,rcb,Pf,Pb] = pc2rcv(pc(:,:,1:m+1),R0);

cov(:,:,1) = R0;
par(:,:,1) = I;
parb(:,:,1)= I;
if m,
	cov(:,:,2) = -rc(:,:,2)*R0;
	par(:,:,2) = rc(:,:,2);
	parb(:,:,2)= rcb(:,:,2);
	par_o  = par;
	parb_o = parb;
end

for p = 2:m,

   par(:,:,2:p) = par_o(:,:,2:p) +flipdim(timesv(rc(:,:,p+1),parb_o(:,:,2:p)),3);
   par(:,:,p+1) = rc(:,:,p+1);

   parb(:,:,2:p)= parb_o(:,:,2:p)+flipdim(timesv(rcb(:,:,p+1) ,par_o(:,:,2:p)),3);
   parb(:,:,p+1)= rcb(:,:,p+1);

   cov(:,:,p+1) = -prodsumv(par(:,:,2:p+1),cov(:,:,p:-1:1));

   par_o  = par;
   parb_o = parb;
end %for p = 2:m,
if ncov > order,
	cov(:,:,order+2:end) = armafilterv(zeros(dim,dim,ncov-order),par,I,cov(:,:,2:order+1));
end

author: Sable

repo: mcbench-benchmarks


text:
function plotAsProbe(data, ann_percents, yc, cm, cm2, sqSizeY, active_site_start, rpl)
% data is vector, same lenght as xc and yc
% cm is colormap, size [n x 3]
%
% other usage: if data is empty and cm's first dimension is the same length
% as xc and yc, then use those colors literally
sqCoordsY = [sqSizeY/2 sqSizeY/2 -sqSizeY/2 -sqSizeY/2];

for q = 1:size(ann_percents,1)

    % get x-boundaries corresponding to percent certainty of region
    percent_region = [0 ann_percents(q,1)*100 ann_percents(q,1)*100 0];
    percent_second_region = [100-ann_percents(q,2)*100 100 ...
                                        100 100-ann_percents(q,2)*100];
    y = yc(q)+sqCoordsY;

    if yc(q) < active_site_start  || yc(q) > rpl*10
        cur_alpha = .4;
    else
        cur_alpha = .8;
    end

    fill(percent_region,y,cm(q,:), 'EdgeAlpha', 0, 'FaceAlpha', cur_alpha);
    fill(percent_second_region,y,cm2(q,:), 'EdgeAlpha', 0, 'FaceAlpha', cur_alpha);

    if ann_percents(q,1) + ann_percents(q,2) < .995
        percent_in_between = [ann_percents(q,1)*100 100-ann_percents(q,2)*100 ...
                                        100-ann_percents(q,2)*100 ann_percents(q,1)*100];
        fill(percent_in_between,y,[.8 .8 .8], 'EdgeAlpha', 0, 'FaceAlpha', cur_alpha);
    end

    hold on;
end

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF


text:
function varargout = vertcat(varargin)

% function varargout = vertcat(varargin)
% Prevent vertcat for cfg_item objects.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: vertcat.m 1862 2008-06-30 14:12:49Z volkmar $

rev = '$Rev: 1862 $'; %#ok

cfg_message('matlabbatch:cfg_item:cat', ['Concatenation of cfg_item objects is ' ...
                    'not allowed.']);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%  clears the page
%
%  ::
%
%     document.clearpage();
%
%  Args:
%     none
%
%

author: jmaih


text:
function i4mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% I4MAT_PRINT_SOME prints out a portion of an I4MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, integer A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, a title.
%
  incx = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );

  if ( m <= 0 || n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (None)\n' );
    return
  end

  for j2lo = max ( jlo, 1 ): incx : min ( jhi, n )

    j2hi = j2lo + incx - 1;
    j2hi = min ( j2hi, n );
    j2hi = min ( j2hi, jhi );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col: ' );

    for j = j2lo : j2hi
      fprintf ( 1, '%7d  ', j );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row\n' );
    fprintf ( 1, '\n' );

    i2lo = max ( ilo, 1 );
    i2hi = min ( ihi, m );

    for i = i2lo : i2hi

      fprintf ( 1, '%5d: ', i );
      for j = j2lo : j2hi
        fprintf ( 1, '%7d  ', a(i,j) );
      end

      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
% [threshx,cost,classorder] = onedim2kclass(x,idx,...)
%
% finds the k-1 thresholds that best separate the input one-dimensional
% data vector x into k classes
%
% input:
%
% x: n x 1 is the n one-dimensional data points
% idx: n x 1 is the integer labels for each data point
%
% optional:
%
% 'issorted': true/false, whether the list is sorted already
% 'classorder': 1 x k vector where classorder(1) comes first and
% classorder(k) comes last. otherwise, classorder is set based on means of each
% class
%
% output:
% threshx: k+1x1, class classorder(i) is >= threshx(i) and < thresh(i+1)
% cost: number of misclassified samples
% classorder: 1 x k
%
function [threshx,cost,classorder] = onedim2kclass(x,idx,varargin)

x = x(:);
idx = idx(:);
n = length(x);

[issorted,classorder] = myparse(varargin,'issorted',false,'classorder',nan);

%k = max(idx); % right now, assuming labels are 1,...,k

% make labels 1,...,k
oldidx = idx;
[newidx2oldidx0,tmp,idx] = unique(idx);
%[tmp,oldidx2newidx0] = sort(newidx2oldidx0);
k = length(newidx2oldidx0);
if ~any(isnan(classorder)),
  classorder0 = classorder;
  for i = 1:k,
    classorder(classorder0==newidx2oldidx0(i)) = i;
  end
end

if k == 1,
  threshx = [min(x),max(x)+eps];
  cost = 0;
  if isnan(classorder),
    classorder = newidx2oldidx0;
  end
  return;
end

%if k == 2,
%  % use onedim2class
%  % idx should be -1,1
%  idx = 2*idx-3;
%  [thresh,lower,cost] = onedim2class(x,idx);
%  if lower == -1,
%    order = [1,2];
%  else
%    order = [2,1];
%  end
%  classorder = order;
%  %classorder = newidx2oldidx(order);
%  return;
%end

% choose order for classes based on order of means
if any(isnan(classorder)),

  mu = zeros(1,k);
  for k1 = 1:k,
    mu(k1) = mean(x(idx==k1));
  end
  [sortedmu,order] = sort(mu);
  [tmp,oldidx2newidx] = sort(order);
  classorder = newidx2oldidx0(order);
  % relabel
  idx0 = idx;
  for k1 = 1:k,
    idx(idx0==k1) = oldidx2newidx(k1);
  end

else
  [tmp,oldidx2newidx] = sort(classorder);
end

if ~issorted,
  [x,xorder] = sort(x);
  idx = idx(xorder);
  [tmp,xreorder] = sort(xorder);
end
% deal with duplicates: cannot have a class start at i if x(i) == x(i-1)
% isdup(i) = x(i) == x(i-1)
isdup = [x(1:end-1) == x(2:end);false];
dupidx = find(isdup);

% costs(n1+1,k1) is the min cost of classifying samples 1:n1 as classes 1:k1
costs = inf(n+1,k-1);
% prev(n1+1,k1)+1 is the start of class k1 when classifying samples 1:n1
prev = zeros(n+1,k-1);

% base cases

% for i = 1,...,n
% cost(i+1,1) is cost of classifying 1:i as 1 is
% sum(idx(1:i)~=k) = c(i+1)
% for i = 0, cost(i+1,1) = 0
costs(2:n+1,1) = cumsum(double(idx~=1));
% cannot end at dupidx - 1
costs(dupidx+1) = inf;

% for k1 = 1,...,k cost of classifying 1:0 as 1:k1 is 0
costs(1,:) = 0;

% compute for increasing k, up to k-1
for k1 = 2:k-1,

  c = [0;cumsum(double(idx~=k1))];

  for n1 = 1:n,
    if isdup(n1), continue; end
    % compute one-class cost for i = 1:n1+1, end = n1
    oneclasscost = c(n1+1)-c(1:n1+1);
    % prev(n1+1,k1) is the start of the last class
    [costs(n1+1,k1),prev(n1+1,k1)] = min( costs(1:n1+1,k1-1) + oneclasscost );
  end

end

% compute for n,k
c = [0;cumsum(double(idx~=k))];
oneclasscost = c(n+1)-c(1:n+1);
[cost,prevlast] = min( costs(1:n+1,k-1) + oneclasscost );

% get thresholds
% class k1 goes from thresh(k1) to thresh(k1+1)-1
thresh = ones(k+1,1);
thresh(end) = n+1;
for k1 = k:-1:1,
  if prevlast < 1,
    thresh(1:k1) = 1;
    break;
  end
  thresh(k1) = prevlast;
  if k1 > 1,
    prevlast = prev(prevlast,k1-1);
  end
end

% convert from indices to data points
threshx = zeros(k+1,1);
threshx(thresh<=1) = x(1);
threshx(thresh==n+1) = x(n)+1;
other = thresh > 1 & thresh <= n;
threshx(other) = (x(thresh(other)-1)+x(thresh(other)))/2;
%threshx(other) = x(thresh(other));

if exist('classorder0','var')
  classorder = classorder0;
end

% here is some code for testing
if 0,
n = 100;
k = 3;
ntests = 100;
sig = .1;
for test = 1:ntests,
  %x = sort(rand(n,1));
  idxtrue = randsample(k,n,true);
  mutrue = sort(rand(k,1));
  x = randn(n,1)*sig + mutrue(idxtrue);
  [x,order] = sort(x);
  idxtrue = idxtrue(order);
  tic;
  [mu,cost,idxdp,threshx] = onedimkmeans(x,k);
  t = toc;
  fprintf('DP score = %f, diff = %e, time = %f, centers = ',cost,cost-sum(sumD),t);
  fprintf('%f, ',mu);
  fprintf('\n');
  if cost > sum(sumD),
    if all(idxkm==idxdp),
      fprintf('Bigger cost, diff = %e, but all idx are the same\n',cost-sum(sumD));
      pause(1);
    else
      fprintf('*****Bigger cost, diff = %e, differing indices\n',cost-sum(sumD));
      fprintf('Differing indices: ');
      fprintf('%d, ',find(idxkm~=idxdp));
      fprintf('\n');
      break;
    end
  end
  diffcost(test) = cost-sum(sumD);
  fprintf('\n');
end
end

author: kristinbranson

repo: JAABA


text:
global iiwa;
global isConnected;
global read_state_var;
iiwa=[];
isConnected=false;
read_state_var=false;
%% Feedback values
global feedback_eef_pos;
global feedback_jpos;
feedback_eef_pos=[];
feedback_jpos=[];
%% Interface container
main_frame_width=1200;
main_frame_height=600;
main_frame_x0=10;
main_frame_y0=100;
daCaption='LBR_GUI_interface';
main_frame_h=figure('Name',daCaption,'rend','painters','pos',[main_frame_x0 main_frame_y0 main_frame_width main_frame_height]);
fcn=@closeMainScript;
set(main_frame_h,'CloseRequestFcn',fcn);  % add confirmation dialog to closing the GUI
pause(1);
%% Hand guiding home interface
interface_handguiding_home(main_frame_h)
%% Connection interface
connection_panel_width=400;
connection_panel_height=250;
connection_panel_margin=20;
interface_connection(main_frame_h,main_frame_width,main_frame_height,connection_panel_width,...
    connection_panel_height,connection_panel_margin);
%% Feedback interface
feedback_panel_width=500;
feedback_panel_height=350;
feedback_panel_margin=20;
interface_feedback_panel(main_frame_h, main_frame_width,main_frame_height,feedback_panel_width,...
    feedback_panel_height,feedback_panel_margin);
%% Buttons for adding points to command line
motionCommand_panel_width=220;
motionCommand_panel_height=230;

xPos=main_frame_width-motionCommand_panel_width-feedback_panel_width-2*feedback_panel_margin;
yPos=main_frame_height-motionCommand_panel_height-feedback_panel_margin;

interface_motionCommand_panel(main_frame_h, main_frame_width,main_frame_height,...
    xPos,yPos,motionCommand_panel_width,...
    motionCommand_panel_height);
%% Command line interface
commandLine_width=500;
commandLine_height=200;
commandLine_margin=20;
interface_instructions_commandLine(main_frame_h, main_frame_width,main_frame_height,commandLine_width,...
    commandLine_height,commandLine_margin );
%% Buttons to execute commandLine program
h=main_frame_h;
executeCommandLine_panel_width=motionCommand_panel_width;
executeCommandLine_panel_height=commandLine_height;
xPos=main_frame_width-commandLine_margin-commandLine_width-commandLine_margin-executeCommandLine_panel_width;
yPos=commandLine_margin;
interface_executeCommandLine_panel(h,xPos,yPos,executeCommandLine_panel_width,...
    executeCommandLine_panel_height)
%% Buttons to save/load a program from commandLine
h=main_frame_h;
loadSave_commandLine_Panel_width=motionCommand_panel_width/2;
loadSave_commandLine_Panel_height=commandLine_height;
xPos=main_frame_width...
    -commandLine_margin-commandLine_width...
    -commandLine_margin-executeCommandLine_panel_width...
    -loadSave_commandLine_Panel_width-commandLine_margin;
yPos=commandLine_margin;
interface_save_load_CommandLine_program_panel(h,xPos,yPos,loadSave_commandLine_Panel_width,...
    loadSave_commandLine_Panel_height)


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcImTool(Temp,Event,handles)
handles = guidata(handles.MU_matrix_display);
try
    BMatrix=handles.BMatrix;
    imtool(BMatrix);
    colormap(handles.V.Color_map);
    set(gca,'Clim',[handles.V.C_lower handles.V.C_upper]);
catch me
    error_msg{1,1}='ERROR!!! imtool is not working in this Matlab version.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function y = huber_circ( x, DIM, varargin )

%HUBER_CIRC   Huber penalty function with circular symmetry.
%   For a vector X, HUBER_CIRC(X) computes the Huber penalty function
%
%       HUBER_CIRC(X) =   NORM(X,2)^2 if NORM(X,2)<=1,
%                       2*NORM(X,2)-1 if NORM(X,2)>=1.
%
%   For matrices and N-D arrays, the penalty function is applied to the
%   first dimens
%
%   HUBER_CIRC(X,[],M) computes the penalty function with halfwidth M,
%   M.^2.*HUBER_CIRC(X./M). M must be real and positive.
%
%   HUBER_CIRC(X,[],M,T) computes the penalty function with halfwidth M
%   and concomitant scale T:
%
%       HUBER_CIRC(X,[],M,T) = T.*HUBER_CIRC(X./T,[],M) if T > 0
%                              +Inf                     if T <= 0
%
%   See the help file for HUBER for information about this usage.
%
%   If X is a matrix, the penalty function is applied to the columns of the
%   matrix X, and a row vector is returned. If X is an N-D matrix, the
%   penalties are computed along the first non-singleton dimension.
%
%   HUBER_CIRC(X,DIM), HUBER_CIRC(X,DIM,M), and HUBER_CIRC(X,DIM,M,T)
%   computes the penalty along the dimension DIM.
%
%   Disciplined convex programming information:
%       HUBER_CIRC is jointly convex in X and T. It is nonomonotonic in X
%       and nonincreasing in T. Therefore, when used in CVX specifications,
%       X must be affine and T must be concave (or affine). T must be real.
%       X, on the other hand, may be real or complex.

if ~cvx_isaffine( x ),
    error( 'Disciplined convex programming error:\n    HUBER_CIRC is nonmonotonic in X, so X must be affine.', 1 ); %#ok
end
if nargin < 2, DIM = []; end
y = huber_pos( norms( x, DIM ), varargin{:} );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.


text:
function [xc,xp,nx,ny] = ComputeStripes(fname,graphout);

%[xc,xp] = ComputeStripes(fname,graphout)
%
% This function computes the projector stripe coordinate (at subpixel
% accuracy) at every pixel in the image. The algorithm used in temporal
% processing (with a translationnal pattern of period 32 pixels). A coarse
% to fine pattern projection helpresolve for the period ambiguity (in a
% Gray-Code Scheme).
%
% The naming convention is,
%
% 	fname_pat##p.ras    --> for the positive image from pass ##
% 	fname_pat##n.ras	--> for the negative image from pass ##
%
% 	##=00: Black and White images
%
%   ##=[01 - 10] : Coarse to fine projection (Gary-Scale projection)
%   ##=[11 - 42 ]: 32 translational patterns of period 32 pixels (for
%                   temporal processing)
% INPUT:
%	fname   --> base name of the images
%   graphout --> Set to 1 to show graphical figures
%
% OUTPUT:
%	xc  is a 2 by N matrix of the point in the image plane
%       (convention: (0,0) -> center of the upper left pixel in the camera image)
%	xp	is a N vector of the corresponding projector stripe numbers (at subpixel accuracy).
%       (convention: (0,0) -> center of the upper left pixel in the projector image)
%   (nx,ny): Size of the camera image
%
%	(c) in 1996 by Jean-Yves Bouguet - Updated 11/26/2003

if nargin < 2,
    graphout = 0;
end;

N = 10;

%-- Load the white and black images:

blackIm = imread([fname '_pat00p.bmp']);
whiteIm = imread([fname '_pat00n.bmp']);

if size(blackIm,3) > 1,
    blackIm = 0.299 * double(blackIm(:,:,1)) + 0.5870 * double(blackIm(:,:,2)) + 0.114 * double(blackIm(:,:,3));
    whiteIm = 0.299 * double(whiteIm(:,:,1)) +  0.5870 * double(whiteIm(:,:,2)) + 0.114 * double(whiteIm(:,:,3));
else
    blackIm = double(blackIm);
    whiteIm = double(whiteIm);
end;

[ny,nx] = size(blackIm);

%%% Contrast Thresholding

%% good value for cthresh: 20

totContr = whiteIm - blackIm;
totContr(totContr < 0) = 0;

cthresh = 3; %--> totContr is larger than cthresh for valid pixels

% In order to remove the highlights (reject image regions where whiteIm >254)
hightlight_reject = 1;

%-- Enable processing of a small region of the image instead of the whole image:

xs = 1;
xe = nx;
ys = 1;
ye = ny;

totContr = totContr(ys:ye,xs:xe);
whiteIm = whiteIm(ys:ye,xs:xe);
blackIm = blackIm(ys:ye,xs:xe);

[yPixels xPixels] = size(totContr);

good1 = find((totContr > cthresh)&(whiteIm <= 254)); % the first mask!!! (no need to compute anything outside of this mask)
Ng1 = length(good1);

%--------------------------------------------------------------------------
%-- STEP 1: TEMPORAL PROCESSING -> Finding the edge time at every pixel in the image
%--------------------------------------------------------------------------

period = 32; % Total Period size in pixels of the projected pattern

index_list2 = (0:period-1) + N + 1;

% Read all temporal images (and compute max and min images):

for kk=0:period-1,

    tmp = imread([fname '_pat' sprintf('%.2d',index_list2(kk+1)) 'p.bmp']);

    if size(tmp,3) > 1,
        tmp = 0.299 * double(tmp(ys:ye,xs:xe,1)) + 0.5870 * double(tmp(ys:ye,xs:xe,2)) + 0.114 * double(tmp(ys:ye,xs:xe,3));
    else
        tmp = double(tmp(ys:ye,xs:xe));
    end;

    eval(['I_' num2str(kk) '= tmp;']);

    if kk == 0,
        Imin = tmp;
        Imax = tmp;
    else
        Imin = min(Imin,tmp);
        Imax = max(Imax,tmp);
    end;

end;

% Substract opposite images (to compute a zero crossing):
for kk = 0:period-1,

    eval(['I_kk = I_' num2str(kk) ';']);
    eval(['I_kk2 = I_' num2str(mod(kk+period/2,period)) ';']);

    eval(['J_' num2str(kk) ' = I_kk - I_kk2;']);

end;

% Start computing the edge points:

not_computed = ones(yPixels,xPixels);
xp_crossings = -ones(yPixels,xPixels);

for kk = 1:period,

    eval(['Ja = J_' num2str(mod(kk-3,period)) ';']);
    eval(['Jb = J_' num2str(mod(kk-2,period)) ';']);
    eval(['Jc = J_' num2str(mod(kk-1,period)) ';']);
    eval(['Jd = J_' num2str(mod(kk,period)) ';']);
    eval(['Je = J_' num2str(mod(kk+1,period)) ';']);
    eval(['Jf = J_' num2str(mod(kk+2,period)) ';']);

    % Temporal Smoothing: (before zero crossing computation)
    J_current = (Jb + 4*Jc + 6*Jd + 4*Je + Jf)/16;
    J_prev = (Ja + 4*Jb + 6*Jc + 4*Jd + Je)/16;

    ind_found = find( (J_current >= 0) & (J_prev < 0) & (not_computed) );

    J_current = J_current(ind_found);
    J_prev = J_prev(ind_found);

    xp_crossings(ind_found) = (kk - (J_current ./ (J_current - J_prev))) - 0.5;

    not_computed(ind_found) = zeros(length(ind_found),1);

end;

% Final temporal solution:
xp_crossings = mod(xp_crossings,period);

if graphout,
    figure(3);
    image(xp_crossings*8);
    colormap(gray(256));
    title('STEP1: Subpixel projector coordinate with a 32 pixel ambiguity');
    drawnow;
end;

%--------------------------------------------------------------------------
%-- STEP 2: SPATIAL PROCESSING -> Coarse to fine processing to resolve ambiguity
%--------------------------------------------------------------------------

bin_current = zeros(yPixels,xPixels);

num_period = zeros(yPixels,xPixels);

for i = 1:N-log2(period),

    tmpN = imread([fname '_pat' sprintf('%.2d',i) 'p.bmp']);
    tmpI = imread([fname '_pat' sprintf('%.2d',i) 'n.bmp']);

    if size(tmpN,3) > 1,
        tmpN = 0.299 * double(tmpN(ys:ye,xs:xe,1)) + 0.5870 * double(tmpN(ys:ye,xs:xe,2)) + 0.114 * double(tmpN(ys:ye,xs:xe,3));
        tmpI = 0.299 * double(tmpI(ys:ye,xs:xe,1)) + 0.5870 * double(tmpI(ys:ye,xs:xe,2)) + 0.114 * double(tmpI(ys:ye,xs:xe,3));
    else
        tmpN = double(tmpN(ys:ye,xs:xe));
        tmpI = double(tmpI(ys:ye,xs:xe));
    end;

    diffI = (tmpN-tmpI)>0;

    bin_current = xor(bin_current,diffI);

    num_period = num_period + (2^(N-i))*bin_current;

end;

if graphout,
    figure(4);
    image(num_period/4);
    colormap(gray(256));
    title('STEP2: Period number (for removing the periodic ambiguity)');
    drawnow;
end;

% Finish off the spatial processing to higher resolution:

xp_spatial = num_period;

finer_image = 2;

for i = N-log2(period)+1:N-finer_image+1,

    tmpN = imread([fname '_pat' sprintf('%.2d',i) 'p.bmp']);
    tmpI = imread([fname '_pat' sprintf('%.2d',i) 'n.bmp']);

    if size(tmpN,3) > 1,
        tmpN = 0.299 * double(tmpN(ys:ye,xs:xe,1)) + 0.5870 * double(tmpN(ys:ye,xs:xe,2)) + 0.114 * double(tmpN(ys:ye,xs:xe,3));
        tmpI = 0.299 * double(tmpI(ys:ye,xs:xe,1)) + 0.5870 * double(tmpI(ys:ye,xs:xe,2)) + 0.114 * double(tmpI(ys:ye,xs:xe,3));
    else
        tmpN = double(tmpN(ys:ye,xs:xe));
        tmpI = double(tmpI(ys:ye,xs:xe));
    end;

    diffI = (tmpN-tmpI)>0;

    bin_current = xor(bin_current,diffI);

    xp_spatial = xp_spatial + (2^(N-i))*bin_current;

end;

% Final spatial solution:
xp_spatial = xp_spatial + (2^(finer_image-1)-1);

%--------------------------------------------------------------------------
%-- STEP 3: Solve for periodic ambiguity, and fixing gliches of the temporal processing (due to noise)
% In order to compare  xp_spatial and  xp_crossings; Not discussed in class
%--------------------------------------------------------------------------

% Fix glitches at the stripe boundaries (of width 4 pixels):

for kkk = 1:10,
    pos_cand = ((num_period == ([1e10*ones(yPixels,1) num_period(:,1:end-1)]+period)) | (num_period == ([1e10*ones(yPixels,2) num_period(:,1:end-2)]+period)))&(xp_crossings > 5*period /6);
    neg_cand = ((num_period == ([num_period(:,2:end) zeros(yPixels,1)]-period)) | (num_period == ([num_period(:,3:end) zeros(yPixels,2)]-period)))&(xp_crossings < period /6);
    num_period = num_period - period*pos_cand + period*neg_cand;
end;

xp_crossings2 = xp_crossings + num_period;

period3 = period / 2;

% Fix the little glitch at the stripe boundaries:

% Find single glitches:

for kkk = 1:5,
    delta_x =  xp_crossings2(:,2:end)-xp_crossings2(:,1:end-1);

    pos_glitch = (delta_x > 3*period3/4)&(delta_x < 3*period3);
    neg_glitch = (delta_x < -period3/4)&(delta_x > -3*period3);
    no_glitch = ~pos_glitch & ~neg_glitch;

    % Place to subtract a period:
    sub_places = [ (neg_glitch & [zeros(yPixels,1) pos_glitch(:,1:end-1)]) zeros(yPixels,1)] ;
    add_places = [ (pos_glitch & [zeros(yPixels,1) neg_glitch(:,1:end-1)]) zeros(yPixels,1)] ;

    xp_crossings3 = xp_crossings2 - period3 * sub_places + period3 * add_places;

    xp_crossings2 = xp_crossings3;

end;

% Find double glitches:

for kkk = 1:5,
    delta_x =  xp_crossings2(:,2:end)-xp_crossings2(:,1:end-1);
    pos_glitch = (delta_x > 3*period3/4)&(delta_x < 3*period3);
    neg_glitch = (delta_x < -period3/4)&(delta_x > -3*period3);
    no_glitch = ~pos_glitch & ~neg_glitch;

    sub2_places = [((no_glitch)& [zeros(yPixels,1) pos_glitch(:,1:end-1)] & [neg_glitch(:,2:end) zeros(yPixels,1)])|(neg_glitch & [zeros(yPixels,1) no_glitch(:,1:end-1)] & [zeros(yPixels,2) pos_glitch(:,1:end-2)]) zeros(yPixels,1)];
    add2_places = [((no_glitch)& [zeros(yPixels,1) neg_glitch(:,1:end-1)] & [pos_glitch(:,2:end) zeros(yPixels,1)])|(pos_glitch & [zeros(yPixels,1) no_glitch(:,1:end-1)] & [zeros(yPixels,2) neg_glitch(:,1:end-2)]) zeros(yPixels,1)];

    xp_crossings3 = xp_crossings2 - period3 * sub2_places + period3 * add2_places;

    xp_crossings2 = xp_crossings3;
end;

% End fix

if graphout,
    figure(5);
    image(xp_crossings2/4);
    colormap(gray(256));
    title('STEP3: Subpixel projector coordinate without the 32 pixel ambiguity');
    drawnow;
end;

%--------------------------------------------------------------------------
%-- STEP 4: Compare xp_spatial and xp_crossings2 and retains the
% xp_crossings that are valid
%--------------------------------------------------------------------------

%comparison of spatial and temporal xp for rejecting the bad pixels:
err_xp = xp_spatial - xp_crossings2;

spatial_temporal_agree = (err_xp <= 2^(finer_image-1)) & (err_xp > -1);

mask_temporal = zeros(yPixels,xPixels);
mask_temporal(2:(yPixels-1),2:(xPixels-1)) = ones(yPixels-2,xPixels-2);

if hightlight_reject,
    highlight = (whiteIm > 254);
    mask_good = (totContr > cthresh) & (not_computed >= 0) & mask_temporal & spatial_temporal_agree & ~highlight;
else
    mask_good = (totContr > cthresh) & (not_computed >= 0) & mask_temporal & spatial_temporal_agree;
end;

good1 = find(mask_good);

xp_crossings3 = xp_crossings2;

xp_crossings3(~mask_good) = NaN;

if graphout,
    figure(6);
    image(xp_crossings3/4);
    colormap(gray(256));
    title('STEP4: Final clean subpixel projector coordinates');
    drawnow;
end;

%--------------------------------------------------------------------------
%-- STEP 5: Produce the camera coordinates and the projector coordinates xc, xp
%--------------------------------------------------------------------------

% extract the good pixels only:
xp = xp_crossings2(good1)';

[X,Y] = meshgrid(0:xPixels-1,0:yPixels-1);

xc = [X(good1)';Y(good1)'];

%%% Express the coordinates of the points in the original
%%% image coordinates:

xc(1,:) = xc(1,:) + xs - 1;
xc(2,:) = xc(2,:) + ys - 1;

repo: ekfmonoslam


text:
function [alpha, c] = dlsqrat(t,y,p,q,alpha)

%  [alpha] = dlsqrat(t,y,p,q,alpha)
%
%  A Full Newton non-linear least-squares code for discrete
%  least-squares rational approximation. This code implements the algorithm
%  described in the paper:
%
% Carlos F. Borges, A Full-Newton Approach to Separable Nonlinear Least
% Squares Problems and its Application to Discrete Least Squares Rational
% Approximation, Electronic Transactions on Numerical Analysis, Volume 35,
% pp.57-68, 2009.
%
% All are welcome to use this code as they wish. I only ask that you cite
% the paper above if you do.
%
%Inputs:
% - t,y are the data points.
% - p,q are the degrees of the numerator and denominator.
% - alpha (optional) is the starting guess
%
%Outputs:
% - alpha contains the denominator coefficients starting with alpha_1
% - c contains the numerator coefficients starting with c_0
%
%  Please note that the polynomial coefficients are generated in ascending
%  order so if you want to use Matlab's polyval routine to evaluate things
%  you need to flip the c vector, and you need to flip the alpha vector and
%  then append a 1. Here is a code fragment you can use to view the results
%  of the fit:
%
%    cla;
%    plot(t,y,'b.'); hold on
%    tt = linspace(min(t),max(t),1000)';
%    yy = polyval(flipud(c),tt)./polyval([flipud(alpha); 1],tt);
%    plot(tt,yy); hold off;
%
%
% Copyright (c) 2008 by Carlos F. Borges.
% All rights reserved.
% begin dlsqrat

% Set the convergence tolerance.
TOLERANCE = 10^(-12);

% N is the Vandermonde that will be used to evaluate the numerator.
N = zeros(length(t),p+1);
N(:,1) = ones(length(t),1);
for k=2:p+1
    N(:,k) = N(:,k-1).*t;
end
% M is the Vandermonde that will be used to evaluate the denominator.
M = zeros(length(t),q);
M(:,1) = t;
for k=2:q
    M(:,k) = M(:,k-1).*t;
end

% If we are not given an initial guess then generate one.
if nargin < 5
    tmp_pade = [N -diag(y)*M]\y;
    alpha = tmp_pade(p+2:end);
end

% Construct the model matrix and compute ancillary quantities.
update(alpha);

for iter=1:100

    % Update the error.
    old_err = err;

    % Compute the Jacobian and the Hessian.
    Tmp1 = diag(Py.*D)*M;
    Tmp2 = Q'*diag((Py-r).*D)*M;
    J = Tmp1 - Q*Tmp2;
    H = M'*diag((Py-2*r).*D)*Tmp1 - Tmp2'*Tmp2;

    % Compute the gradient.
    gradient = J'*r;

    % Compute the Cholesky factorization of H.
    [R, not_PD] = chol(H);
    % If H is not positive definite then regularize and factor
    if not_PD
        R = chol(H - 1.2*min(eig(H))*eye(q));
    end

    %Compute the Newton step.
    delta = -R\(R'\gradient);

    % Use stepsize control to take a step.
    step_control;

    % Convergence testing
    if err > old_err
        disp('Failed to find descending step length.');
        break;
    else
        alpha = new_alpha;
        rel_err = abs(old_err - err)/old_err;
        if rel_err <= TOLERANCE
            break;
        end
    end
    % End convergence testing.

end  %End of main loop.

% Compute the coefficients of the numerator.
c = (diag(D)*N)\y;

% Generate an error message if the algorithm failed to converge.
if rel_err > TOLERANCE
    disp('Algorithm did not converge.');
end

%XXXXXXXXXXXXXXXXX Subroutines XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    function update(alpha)
        % Updates the model matrix and computes ancillary quantities.
        D = 1./(1+M*alpha);         % Compute the denominator.
        [Q R] = qr(diag(D)*N,0);    % Compute the QR factorization of A = D*N
        Py = Q*(Q'*y);              % Compute the projection of y onto the range of A.
        r = y - Py;                 % Compute the residual.
        err = r'*r;                 % Compute the current squared error.
    end

    function step_control
        % This function implements stepsize control using a simple
        % backtracking scheme from Dennis & Schnabel.

        % Try taking a full step.
        new_alpha = alpha + delta;

        % Update the model.
        update(new_alpha);

        % If a full step does not sufficiently reduce the error then we
        % use a backtracking line-search method for step-size control.
        % This involves minimizing a function f(lambda) that interpolates the
        % computed error (and its derivatives) at different values of lambda.
        f0 = old_err;
        fprime = gradient'*delta;
        steptol = f0 + .0001*fprime;
        if err > steptol

            errs(1) = err; lams(1) = 1;   % We'll need this if further refinement is necessary.

            % We start with a quadratic model at f(0), f'(0), and f(1)
            % and will take the larger of the computed step or 1/10.
            lambda = max([-fprime/(2*(err - f0 - fprime)) .1]);

            new_alpha = alpha + lambda*delta;

            % Update the model matrix and compute ancillary quantities.
            update(new_alpha);

            % If this doesn't work then we loop with a cubic model at f(0),
            % f'(0), f(lambda), and f(lam2) where the last two are errors at
            % the last two lambda that were tried.
            steptol = f0 + .0001*fprime*lambda;
            while err > steptol

                % Push the current lambda and error to the top of the lams and errs
                % stacks.
                lams = [lambda; lams(1)]; errs = [err; errs(1)];
                rhs = (errs - fprime*lams - [f0 ; f0])./(lams.*lams);
                ab = [lams [1 ; 1]]\rhs;

                lambda = (-ab(2)+sqrt(ab(2)*ab(2) - 3*ab(1)*fprime))/(3*ab(1));

                % It is still important to make certain that the new lambda
                % progresses quickly but not too quickly. So if lambda is less
                % than lam2/10 we just use lam2/10, and if it is larger than
                % lam2/2 then we use lam2/2.
                if lambda < lams(1)/10
                    lambda = lams(1)/10;
                end
                if lambda > lams(1)/2
                    lambda = lams(1)/2;
                end

                new_alpha = alpha + lambda*delta;

                % Update the model matrix and compute ancillary quantities.
                update(new_alpha);

                steptol = f0 + .0001*fprime*lambda;

            end
        end
    end

%XXXXXXXXXXXXXXXXX Subroutines End XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

end
% End of function.


text:
% File: VisualizeDataset.m
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

function VisualizeDataset(Dataset)

figure
for i=1:size(Dataset,1)
    img = ShowPose(reshape(Dataset(i,:,:), [10 3]));
    imshow(img);
    pause(0.3);
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
% takes a 3D matrix and returns a cell with length equal to the length of
% the third dimension of the matrix.
function X = lastdim2cell(Y)
X = arrayfun(@(t) Y(:,:,t), 1:size(Y,3), 'uniformoutput', false);
end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function X = Inv_DetailCurveCoeff(C,IsImageReal);
% Inv_DetailCurveCoeff:  Reconstruct the jth dyadic frequency band
%                       from curvelet coefficients at that scale
%  Usage:
%    X = Inv_DetailCurveCoeff(C);
%  Inputs:
%    C    matrix of curvelet coefficients at scale 2^j
%  Outputs:
%    X    jth dyadic subband
%  See Also
%   DetailCurveCoeff, Inv_SeparateAngles, Inv_Curvelet02Xform
%
% By Emmanuel Candes, 2003-2004

  C =  ClockwisetoWENS(C);
  nn = size(C);
  R = zeros(nn);
  deep = log2(nn(2));

  for j = 1:size(R,1),
    for m = 1:size(R,2),
      W = squeeze(C(j,m,:,:));
      W = fft2_mid0(W)/sqrt(prod(size(W)));
      R(j,m,:,:) = W;
    end
  end

  for w=1:size(R,2)
    tmp = squeeze(R(2,w,:,:));
    R(2,w,:,:) = tmp([2:end,1], [2:end,1]);
  end
  for w=1:size(R,2)
    tmp = squeeze(R(4,w,:,:));
    R(4,w,:,:) = tmp([2:end,1], [2:end,1]);
  end

  MaxIts = 25;
  X = Inv_SeparateAngles(Adj_SqueezeAngularFT(R),deep, ...
                         IsImageReal,MaxIts,1e-6,[]);

author: thomaskuestner


text:
classdef WaveBuffer
    properties
        nCh = 1;            % number of channels
        BufferSize = 400;   % this is the buffer size, and also the data size we get from the buffer everytime we call Get().
        BlockSize = 160;    % this is the amount of data we push into the buffer everytime we call Push.
        FrameOverlap = 240; % this is the overlap between two consecutive frames.
        % Note that, the BlockSize = BufferSize - FrameOverlap. In this way,
        % we can gauarantee that the blockwise processing generate the same
        % Short-time Fourier Transform representation of the signal as in
        % the offline processing.
        precision = 'double';
        data = [];
    end

    methods
        function obj = WaveBuffer(nCh, BufferSize, FrameOverlap, precision)
            obj.nCh = nCh;
            obj.BufferSize = BufferSize;
            obj.FrameOverlap = FrameOverlap;
            obj.BlockSize = BufferSize - FrameOverlap;
            if nargin>=4
                obj.precision = precision;
            end

            % initialize the buffer with 0's
            obj.data = zeros(obj.nCh, obj.BufferSize, obj.precision);
        end

        function obj = Push(obj, data)  % push data into the buffer
            [nCh, nSample] = size(data);

            assert(obj.nCh == nCh, sprintf('Error: number of channel in data (%d) is different from number of channel in buffer (%d)', nCh, obj.nCh));
            assert(obj.BlockSize == nSample, sprintf('Error: number of samples in data (%d) is different from block size(%d)', nSample, obj.BlockSize));

            obj.data = [obj.data(:,obj.BlockSize+1:end) data];
        end

        function data = Get(obj)        % get the data from the buffer
            data = obj.data;
        end

        function obj = Reset(obj)       % reset the content of the buffer to 0
            obj.data(:) = 0;
        end
    end

end

author: singaxiong

repo: SignalGraph


text:
%--- help for abstvar/irf ---
%
%  Compute impulse response function from the given parameter values
%
%  ::
%
%     myirfs = irf(self, shock_names, irf_periods, params, Rfunc);
%
%  Args:
%
%     self (var object): var object
%
%     shock_names (cellstr): shocks to compute IRFs (has to be consistent with the names given in :func:`identification <var.identification>`)
%
%     irf_periods (int): number of periods to compute IRFs (default: 40)
%
%     params (vector): parameter values of the model. If empty MLE/posterior-mode values are used.
%
%     Rfunc (function handle): identification function. This is an output of :func:`identification <var.identification>`. (default: choleski identification)
%
%     girf_setup (struct|{empty}): structure containing information relevant
%        for the computation of generalized impulse response functions. If
%        empty, simple regime-specific impulse responses are computed, else
%        girfs are computed. In that case the relevant information to
%        provide in girf_setup is:
%        - nsims : (default=300) number of simulations for the integration.
%        Note that even setting girf_setup=struct() will trigger the
%        computation of girfs. But in that case only the default options
%        will apply.
%
%  Returns:
%
%     : struct containing IRFs
%
%  Note:
%
%     Only successful IRFs are returned. If the structure does not return some IRFs make sure that IRFs properly identified.
%
%
%    Other functions named irf
%
%       dsge/irf    generic/irf
%


text:
function drawArcEllipse(x,y,a,b,theta0,theta1,color)
%(x,y):= coordinates of the center
% r: = radius of the circle
angle = theta0+0.01:0.01:theta1;
xp = x+ a*cos(angle);
yp = y+ b*sin(angle);
plot(xp,yp,color);

% plot the symetric arc
angle = theta0+pi+0.01:0.01:theta1+pi;
xp = x+ a*cos(angle);
yp = y+ b*sin(angle);
plot(xp,yp,color);


text:
function ctf_plot(ctf,CHAN,TIME,TRIALS,Xhair)

% ctf_plot - plot ctf.data
%
% ctf_plot(ctf,CHAN,TIME,TRIALS,Xhair)
%
% CHAN - see ctf_channel_select for options
% TIME - see ctf_read for options (given in msec)
% TRIALS - select 1 trial to plot (the default is trial = 1)
%
% Xhair - 1 to use the crosshair GUI features (default), 0 otherwise
%
%      <>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%      <                                                      >
%      <                    DISCLAIMER:                       >
%      <                                                      >
%      < THIS PROGRAM IS INTENDED FOR RESEARCH PURPOSES ONLY. >
%      < THIS PROGRAM IS IN NO WAY INTENDED FOR CLINICAL OR   >
%      <                    OFFICIAL USE.                     >
%      <                                                      >
%      <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>
%

% $Revision: 1.1 $ $Date: 2009-01-30 03:49:27 $

% Copyright (C) 2004  Darren L. Weber
%
% 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.

% Modified: 02/2004, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('CHAN','var'), CHAN = 'all'; end
if ~exist('TIME','var'), TIME = 'all'; end
if ~exist('TRIALS','var'), TRIALS = 1; end
if ~exist('Xhair','var'), Xhair = 1; end

if isempty(CHAN), CHAN = 'all'; end
if isempty(TIME), TIME = 'all'; end
if isempty(TRIALS), TRIALS = 1; end
if isempty(Xhair), Xhair = 1; end

% This function calls ctf_channel_select
[CHAN,type] = ctf_channel_select(ctf,CHAN);

switch num2str(TIME),
  case 'all',
    TIME = ctf.setup.time_msec;
    TIME_index = 1:ctf.setup.number_samples;
otherwise
    fprintf('...sorry, time restrictions not implemented correctly (03/2004)\n');
    TIME = ctf.setup.time_msec;
    TIME_index = 1:ctf.setup.number_samples;

    %     % assume the input is a range of times in sec
    %     % check the range
    %     if TIME(1) > ctf.setup.time_msec(1),
    %       fprintf('...setting TIME(1) = ctf.setup.time_msec(1)\n');
    %       TIME(1) = ctf.setup.time_msec(1);
    %     end
    %     if TIME(end) > ctf.setup.time_msec(end),
    %       fprintf('...setting TIME(end) = ctf.setup.time_msec(end)\n');
    %       TIME(end) = ctf.setup.time_msec(end);
    %     end
    %     % now find the nearest indices into the samples matrix
    %     TIME_index = interp1(ctf.setup.time_msec,1:ctf.setup.number_samples,TIME,'nearest');
    %     % now ensure that the TIME array is consistent with ctf.setup.time_sec
    %     TIME = ctf.setup.time_msec(TIME_index);
end
TIME = sort(TIME);

switch num2str(TRIALS),
  case 'all',
    TRIALS = 1:ctf.setup.number_trials;
  otherwise
    % assume the input is an array of trials
end
TRIALS = unique(sort(TRIALS));

% check the input data for more than 1 trial
trials = size(ctf.data,3);
if trials > 1,
    fprintf('...plotting trial %d of %d trials in ctf.data\n',TRIALS,trials);
    data = ctf.data(:,:,TRIALS);
else
    data = ctf.data;
end

% now plot the data
switch type,

    case 'all',

        if isempty(ctf.sensor.index.meg_sens),
            fprintf('...no meg sensors\n');
        else
            figure('Name','MEG sensors');
            MEGCHAN = ctf_channel_select(ctf,'meg');
            plot(ctf.setup.time_msec(TIME_index),data(TIME_index,MEGCHAN))
            axis tight
            if exist('crosshair') & Xhair,
                crosshair;
            end
        end

        if isempty(ctf.sensor.index.eeg_sens),
            fprintf('...no eeg sensors\n');
        else
            figure('Name','EEG sensors');
            EEGCHAN = ctf_channel_select(ctf,'eeg');
            plot(ctf.setup.time_msec(TIME_index),data(TIME_index,EEGCHAN))
            axis tight
        end

    case 'meg',

        if isempty(CHAN),
            fprintf('...no meg sensors\n');
        else
            figure('Name','MEG sensors');
            plot(ctf.setup.time_msec(TIME_index),data(TIME_index,CHAN))
            axis tight
        end

    case 'eeg',

        if isempty(CHAN),
            fprintf('...no eeg sensors\n');
        else
            figure('Name','EEG sensors');
            plot(ctf.setup.time_msec(TIME_index),data(TIME_index,CHAN))
            axis tight
        end

    case 'ref',

        if isempty(CHAN),
            fprintf('...no ref sensors\n');
        else
            figure('Name','MEG REF sensors');
            plot(ctf.setup.time_msec(TIME_index),data(TIME_index,CHAN))
            axis tight
        end

    case 'others',

        if isempty(CHAN),
            fprintf('...no other sensors\n');
        else
            figure('Name','Other sensors');
            plot(ctf.setup.time_msec(TIME_index),data(TIME_index,CHAN))
            axis tight
        end

    otherwise

        figure;
        plot(ctf.setup.time_msec(TIME_index),data(TIME_index,CHAN))
        axis tight

end

if exist('crosshair') & Xhair,
    crosshair;
end

return


text:
function imgOut = DurandTMO(img, target_contrast, filter_type)
%
%       imgOut = DurandTMO(img, target_contrast, filter_type)
%
%
%        Input:
%           -img: input HDR image
%           -target_contrast: how to reduce the dynamic range
%
%        Output:
%           -imgOut: tone mapped image
%
%     Copyright (C) 2010-15  Francesco Banterle
%
%     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/>.
%
%     The paper describing this technique is:
%     "Fast Bilateral Filtering for the Display of High-Dynamic-Range Images"
% 	  by Fredo Durand and Julie Dorsey
%     in Proceedings of SIGGRAPH 2002
%

%is it a three color channels image?
check13Color(img);

checkNegative(img);

%default parameters
if(~exist('target_contrast', 'var'))
    target_contrast = 5; %as in the original paper
end

if(~exist('filter_type', 'var'))
    filter_type = 'approx_importance';
end

%compute luminance channel
L = lum(img);

sigma_s = max(size(L)) * 0.02;

%separate detail and base
[Lbase, Ldetail] = bilateralSeparation(L, sigma_s, 0.4, 'log10', filter_type);

eps = 1e-6;
log_base = log10(Lbase + eps);
max_log_base = max(log_base(:));
log_detail = log10(Ldetail);
c_factor = log10(target_contrast) / (max_log_base - min(log_base(:)));
log_absolute = c_factor * max_log_base;

log_Ld = log_base * c_factor + log_detail  - log_absolute;
Ld = 10.^(log_Ld) - eps;
Ld(Ld < 0.0) = 0.0;

%change luminance
imgOut = ChangeLuminance(img, L, Ld);

end

author: banterle


text:
function y = wprctile(x, p, w)
%WPRCTILE  Percentiles of a weighted sample.
%
%  Description
%    Y = PRCTILE(X, P, W) returns percentiles of the values in X
%    (along the first dimension). P is a scalar or a vector of percent
%    values. W is a vector of unnormalized weights for samples. Length
%    of W has to be same as length of X. X need to be a a vector. Y is
%    the same size as P, and Y(i) contains the P(i)-th percentile.
%
%    Example
%      y = prctile(x,50,w); % the median of x given sample weights w
%
%  See also wmean, prctile
%
% Copyright (c) 2000-2013 Aki Vehtari

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

x=sort(x,1);
p=p./100;
y=zeros(length(p),size(x,2));
ww=cumsum(w);ww=ww./ww(end);
for j=1:length(p)
  wi=min(find(ww>=p(j)));
  if wi==1
    y(j,:)=x(1,:);
  else
    w1=ww(wi-1);x1=x(wi-1,:);
    y(j,:)=x1+(x(wi,:)-x1).*(p(j)-w1)./(ww(wi)-w1);
  end
end

author: gpstuff-dev

repo: gpstuff

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


text:
% Matlab script for RAPID equivalent commands
% please note that variables should be
% declared as global as in RAPID syntaxis
% this script demonstrates the use of functions in the matlab script. This
% script can be translated to RAPID by means of the matlab2RAPID script.
%
% IN ORDER TO SIMULATE THE PROGRAM:
%   A) FIRST, LOAD A ROBOT
%       >> robot = load_robot('abb','irb140');
%       OR
%       >> robot = load_robot('abb','irb52');
%   B) NEXT, LOAD SOME EQUIPMENT.
%       >> robot.equipment{1} = load_robot('equipment','tables/table_two_areas');
%
%   C) NOW, LOAD AN END TOOL
%       >> robot.tool= load_robot('equipment','end_tools/parallel_gripper_0');
%   D) FINALLY, LOAD A PIECE TO GRAB BY THE ROBOT
%       >> robot.piece{1}=load_robot('equipment','cylinders/cylinder_tiny');
%
%   E) IF NECESSARY, CHANGE THE POSITION AND ORIENTATION OF THE piece,
%   relative to the robot's base reference system.
%
%       >> robot.piece{1}.T0= [1 0 0 -0.1; 0 1 0 -0.5;0 0 1 0.2;0 0 0 1];
%
%   The A, B, C, D and E steps can be performed usint the 'teach' application
%   by clicking on the Load equipment, load end tool, and load piece
%   buttons.
%
% during the simulation, call simulation_open_tool; to open the tool and
% simulation_close_tool; to close it.
%
% To grip the piece, call simulation_grip_piece; and
% simulation_release_piece to release it.
%
% The call to each function must be in a correct, otherwise, the result may be
% undesired, thus, typically the correct sequence is:
% simulation_open_tool;
% (approach the piece to grab with MoveL or MoveJ)
% simulation_close_tool;
% simulation_grip_piece; --> the piece will be drawn with the robot
%
% move to a different place with MoveJ or MoveL
% simulation_open_tool;
% simulation_release_piece
%
% In RAPID, consider the use of ConfJ on, or ConfJ off, or ConfL on, or ConfL off
% in the case the controller avoids the change of configurations between target points
function practice_2_bolts
global RT_pos_ini RJ_ini RT_aprox_rec RT_pos_rec
global RT_aprox_dej RT_pos_dej TD_gripper VAR_pieza
global robot

%If desired, move the robot with respect to the world reference system
%by modifying robot.T0. Example: robot.T0(2,4)=1;

TD_gripper=[1,[[0,0,0.125],[1,0,0,0]],[0.1,[0,0,0.100],[1,0,0,0],0,0,0]];

RJ_ini=[0,0,0,0,0,0];
RT_pos_ini=[[0.044,-0.501,0.432],[0.262525,0.622008,0.673812,-0.300276],[0,0,0,0],[9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];

RT_aprox_rec=[[0.03,-0.6,0.45],[0.08,0.60337,0.79283,-0.03126],[0,0,0,0],[9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
RT_pos_rec=[[0.0345,-0.629,0.337],[0.261833,0.652321,0.67413,-0.226869],[-1,0,-1,0],[9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];

RT_aprox_dej=[[0.470,-0.460,0.500],[0.07981,0.603204,0.792983,-0.030902],[0,0,0,0],[9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
RT_pos_dej=[[0.470,-0.460,0.450],[0.07981,0.603204,0.792983,-0.030902],[0,0,0,0],[9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];

%please: comment the following lines to avoid loading the robot and equipment in every simulation.
%Feel free to load other robots/equipments... etc.
robot = load_robot('ABB','IRB140');
%load static equipment, such as tables or
robot.equipment{1} = load_robot('equipment','tables/table_two_areas');
robot.tool= load_robot('equipment','end_tools/parallel_gripper_0');
robot.piece{1}=load_robot('equipment','objects/bolt');

%local function to init simulation variables
init_simulation;

VAR_pieza=0;
%!Reset --> abrir pinza
%!Set --> Cerrar pinza

main
end

function main()

global RJ_ini TD_gripper VAR_pieza

%inicializar pieza
VAR_pieza=0;

MoveAbsJ(RJ_ini,'vmax','z100',TD_gripper, 'wobj0');

%MoveJ(RT_pos_ini,'vmax','z100',TD_gripper, 'wobj0');

for i=1:4,
    ASIR_DEJAR();
    VAR_pieza=VAR_pieza +1;
    %local function to init simulation variables
    init_simulation;
end

end

function ASIR_DEJAR()

global  RT_aprox_rec RT_aprox_dej TD_gripper

%!Moverse a posici�n de aproximacion de recogida
MoveJ(RT_aprox_rec,'vmax','z100',TD_gripper, 'wobj0');
COGER_PIEZA;
%!Moverse a posici�n base de dejada
MoveJ(RT_aprox_dej,'vmax','z100',TD_gripper, 'wobj0');
METER_EN_CAJA;

end

function COGER_PIEZA()
global RT_aprox_rec RT_pos_rec TD_gripper

%!Ahora abrir pinza
simulation_open_tool; %Reset do1;
WaitTime(0.1);!esperar apertura
%!Coger
MoveL(RT_pos_rec,'vmax','fine',TD_gripper, 'wobj0');

%!Ahora cerrar pinza
simulation_close_tool; %Set do1;
simulation_grip_piece;

WaitTime(0.1); !esperar cierre
%!Subir por seguridad antes de ir al siguiente punto
MoveL(RT_aprox_rec,'vmax','z50',TD_gripper, 'wobj0');

end

function  METER_EN_CAJA()
global RT_pos_dej TD_gripper VAR_pieza

if VAR_pieza==0
    %Bajar para meter primera pieza
    MoveL(Offs(RT_pos_dej,0,0,-0.105),'v1000','fine',TD_gripper, 'wobj0');
elseif VAR_pieza==1
    MoveL(Offs(RT_pos_dej,-0.055,0,0),'v1000','z100',TD_gripper, 'wobj0');
    MoveL(Offs(RT_pos_dej,-0.055,0,-0.105),'v1000','fine',TD_gripper, 'wobj0');
elseif VAR_pieza==2
    MoveL(Offs(RT_pos_dej,0,-0.055,0),'v1000','fine',TD_gripper, 'wobj0');
    MoveL(Offs(RT_pos_dej,0,-0.055,-0.105),'v1000','fine',TD_gripper, 'wobj0');
else
    MoveL(Offs(RT_pos_dej,-0.055,-0.055,0),'v1000','fine',TD_gripper, 'wobj0');
    MoveL(Offs(RT_pos_dej,-0.055,-0.055,-0.105),'v1000','fine',TD_gripper, 'wobj0');
end

%Ahora abrir pinza
%yes, release the piece
simulation_open_tool; %Reset do1;
simulation_release_piece;

plot_points;

WaitTime(0.5);
%Subir por seguridad antes de ir al siguiente punto
MoveL(Offs(RT_pos_dej,0,0,0.105),'v1000','z10',TD_gripper, 'wobj0');

end

function init_simulation
global robot
robot.piece{1}.T0=eye(4);
robot.tool.Trel=eye(4);
% Now open the tool
simulation_open_tool; %Reset do1;
simulation_release_piece;
%init the position of the piece at the beginning of the simulation
robot.piece{1}.T0(1:3,4)=[0.03 -0.68 0.26]';
th=-30*pi/180;
u=[1 0 0;
   0 cos(th) -sin(th);
   0 sin(th) cos(th)];
robot.piece{1}.T0(1:3,1:3)=u;%eye(3);
%robot.tool.piece_gripped=0;
drawrobot3d(robot, robot.q);

end

function plot_points
x0=0.47;
y0=-0.46;
z0=0.45-0.155;
delta=0.055;
X=[x0 x0-delta x0 x0-delta];
Y=[y0 y0 y0-delta y0-delta];
Z=[z0 z0 z0 z0];

plot3(X', Y', Z', 'LineWidth', 5)
end

repo: ARTE


text:
function c = butter(c,varargin)
   %butter   create and apply a butterworth filter to the traces
   % C = BUTTER(C,[TYPE],CUTOFF,[POLES])
   % This function creates and applies a butterworth filter to each waveform
   % in the correlation object. The function returns a correlation object.
   % Filtering is performed by an underlying call to the filterobject/filtfilt
   % routine. Note that filtfilt applies the filter once in each direction to
   % minimize phase distortion (effectively a zero-phase filter). This also
   % has the effect of doubling the order of the filter. 2 or 4 poles should
   % be sufficient for most applications.
   %
   % EXAMPLES:
   %  c = BUTTER(c,[1 5])            band pass filter on 1-5 Hz (4 poles)
   %  c = BUTTER(c,'B',[1 5])        same as previous example
   %  c = BUTTER(c,'L',5)            low  pass filter below 5 Hz (4 poles)
   %  c = BUTTER(c,'H',1)            high pass filter above 5 Hz (4 poles)
   %  c = BUTTER(c,...,2)            use 2 poles

   % Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$
   % TODO: should check to ensure low cutoff period is >> trace length

   % GET INPUTS
   if ischar(varargin{1})       % set filter type
      type = varargin{1};
      varargin = varargin(2:end);
   else
      type = 'B';
   end

   cutoff = varargin{1};            % set cutoff frequencies

   if length(varargin)>1
      poles = varargin{2};
   else
      poles = 4;                  % set number of poles
   end

   % CHECK FIRST LETTER OF FILTER TYPE
   type = upper(type(1));
   if (type~='B') && (type~='H') && (type~='L')
      error('Filter type not recognized')
   end;

   % CHECK NUMBER OF CUTOFF FREQUENCIES
   if type=='B'
      if length(cutoff)~=2
         error('Two cutoff frequecies needed for bandpass filter')
      end;
   else
      if length(cutoff)~=1
         error('High and lowpass filters require one cutoff frequency')
      end;
   end;

   % CHECK FREQUENCY RANGE
   if cutoff(end) >= c.traces(1).nyquist()
      error('Frequency cutoffs exceed the Nyquist frequency')
   end;

   % duration = get(c,'DURATION_EPOCH');
   % if 1/cutoff(1) >= duration(1)
   %     warning('Frequency cutoff is very low relative to trace length');
   % end;
   % uration(1)
   % 1/cutoff(1)

   % APPLY FILTER
   f = TraceFilter(type,cutoff,poles);
   c.traces = f.filtfilt(c.traces);
end

author: geoscience-community-codes

repo: GISMO


text:
%GETSTRUCTURINGELEMENT  Returns a structuring element of the specified size and shape for morphological operations
%
%     elem = cv.getStructuringElement('OptionName', optionValue, ...)
%
% ## Output
% * __elem__ Output structuring element of specified shape and size.
%
% ## Options
% * __Shape__ Element shape, default 'Rect'. Could be one of:
%   * __Rect__ a rectangular structuring element: `E(i,j)=1`
%   * __Cross__ a cross-shaped structuring element: `E(i,j)=1` if
%     `i=Anchor(2)` or `j=Anchor(1)`, `E(i,j)=0` otherwise.
%   * __Ellipse__ an elliptic structuring element, that is, a filled ellipse
%     inscribed into the rectangle `[0, 0, KSize(1), KSize(2)]`.
% * __KSize__ Size of the structuring element `[w,h]`. default [3,3].
% * __Anchor__ Anchor position within the element. The default value (-1,-1)
%   means that the anchor is at the center. Note that only the shape of a
%   cross-shaped element depends on the anchor position. In other cases the
%   anchor just regulates how much the result of the morphological operation
%   is shifted.
%
% The function constructs and returns the structuring element that can be
% further passed to cv.erode, cv.dilate or cv.morphologyEx. But you can also
% construct an arbitrary binary mask yourself and use it as the structuring
% element.
%
% See also: cv.sepFilter2D, cv.getDerivKernels, cv.getGaussianKernel
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv