text:
function [P] = spm_dcm_fmri_graph_gen(x,v,P)
% Generates adjacency graph for spectral DCM for fMRI
% FORMAT [g] = spm_dcm_fmri_graph_gen(x,v,P)
%
% This routine computes the adjacency matrix (A) for spm_fx_fmri
%
% see also: spm_fx_fmri
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_fmri_graph_gen.m 5821 2013-12-31 14:26:41Z karl $

% compute bias for log connectivity using functional space
%==========================================================================

% spectral deomposition
%--------------------------------------------------------------------------
P.A   = full(P.A);

if isfield(P,'modes')

    % outer product
    %======================================================================
    n     = length(P.A);                          % number of nodes
    m     = length(v);                            % number of modes
    s     = [-exp(-v); (zeros(n - m,1) - 1)];
    P.A   = P.modes*diag(s)*P.modes';
    P.A   = full(P.A + diag(log(-2*diag(P.A)) - diag(P.A)));
    P     = rmfield(P,'modes');

    return

end

if isnumeric(v)

    % static modes (explicit negativity constraints on self excitation)
    %======================================================================
    % P.A   = v'*v;

    % dynamical modes (implicit negative definite constraints)
    %======================================================================
    P.A   = logm(v'*v + eye(size(v,2),size(v,2))*exp(-16))/16;
    P.A   = P.A + diag(log(-2*diag(P.A)) - diag(P.A));

    return

end

% Distance-based bias on (empirical) prior mean of log connectivity
%--------------------------------------------------------------------------
[n m]          = size(v.x);
if size(P.A,3) == 1 && numel(v.a) == 1

    % one-state model of (MoG) connectivity
    %======================================================================
    for i = 1:m
        for j = (i + 1):m

            % Euclidean distance
            %--------------------------------------------------------------
            P.A(i,j) =  ...
                exp(v.a - sum((v.x(:,i) - v.x(:,j)).^2)/2)/4 - ... % excitatory
                exp(v.a - sum((v.x(:,i) - v.x(:,j)).^2)/8)/8;      % inhibitory
            P.A(j,i) = P.A(i,j);

        end
    end

elseif size(P.A,3) == 1 && numel(v.a) == 2

    % one-state model of centre-surround connectivity
    %======================================================================
    for i = 1:m
        for j = (i + 1):m

            % Euclidean distance
            %--------------------------------------------------------------
            D        = exp(-sum((v.x(:,i) - v.x(:,j)).^2)/2);
            P.A(i,j) = exp(v.a(1))*D/16 + v.a(2);
            P.A(j,i) = P.A(i,j);

        end
    end

elseif size(P.A,3) == 2

    % assume two-state model of log connectivity
    %======================================================================
    for i = 1:m
        for j = (i + 1):m

            % Euclidean distance
            %--------------------------------------------------------------
            P.A(i,j,1) = v.a - sum((v.x(:,i) - v.x(:,j)).^2)/2;
            P.A(j,i,1) = P.A(i,j,1);

            % hierarchical distance
            %--------------------------------------------------------------
            P.A(i,j,2) = (sqrt(sum(v.x(:,i).^2)) - sqrt(sum(v.x(:,j)).^2))/2;
            P.A(j,i,2) = -P.A(i,j,2);

        end
    end

end

author: spm


text:
function dq = linevel2dq(v,r,vp,rp)

% LINEVEL2DQ     Transforms a line velocity expressed in vector notation
%                into its dual quaternion representation.
%
%      DQ = LINE2DQ(V,R,VP,RP) transforms the line position, specified by:
%          - the line orientation V
%          - the position of any point of the line, R.
%          - the line orientation rate of change, VP
%          - the velocity component (orthogonal to the line orientation V)
%          of point P, RP
%       V does not need to be unitary, but VP must be orthogonal to V. If
%       RP has a component in the V orientation, it does not matter, since
%       it does not change the expression of the resulting dual quaternion
%       DQ.
%       V,R,VP and RP must have the same size. The inputs (V,R,VP,RP) are
%       either a vector of size 3 or an array of size 3*N (column i
%       represents the input component of Line velocity i) where N is the
%       number of lines. DQ is either a vector of size 8, either an array
%       of size (8*N) depending on the input format. Each column of DQ
%       represents the dual quaternion representation of the corresponding
%       line position velocity.
%
% See also POS2DQ, VEL2DQ, LINE2DQ

sv = size(v);
sr = size(r);
svp = size(vp);
srp = size(rp);
if sv == [1 3],  v = v.'; sv = size(v); end
if sr == [1 3],  r = r.'; sr = size(r); end
if svp == [1 3],  vp = vp.'; svp = size(vp); end
if srp == [1 3],  rp = rp.'; srp = size(rp); end

% check that all inputs have the same size
tab_s = [sv; sr; svp; srp];
if max(tab_s) ~= min(tab_s)
      error('DualQuaternion:linevel2dquat:sizesDoNotMatch',...
        'Arrays v, r, vp and r should be the same size. Size of \n - v is [%d %d] \n - r is [%d %d] \n - vp is [%d %d] \n - rp is [%d %d]',...
           sv(1),sv(2),sr(1),sr(2),svp(1),svp(2),srp(1),srp(2));
end

% if the format is wrong
if sv(1) ~= 3
    error('DualQuaternion:linevel2dquat:wrongsize',...
        '%d rows in the V,R,VP and RP  arrays. It should be 3. ',sv(1));
end

% normalization of the axis vector (if necessary)
n = length(v(1,:));
n2 = sum(v.^2).^0.5;
n2 = repmat(n2,3,1);
v=v./n2;
vp=vp./n2;

% construction of the line velocity dual quaternion
dq = sym(zeros(8,n));
dq(2:4,:) = vp;
dq(6:8,:) = cross(r,vp)+cross(rp,v);

author: Sable

repo: mcbench-benchmarks


text:
function [normalized_speed, actual_speed] = normalize_speed(speed, failures, skipping, tracker, sequence)
% normalize_speed Normalizes tracker speed estimate
%
% This function normalizes speed estimates based on performance profile and some information about
% the way the measurement was obtained (sequence, number of failures, frame skipping).
%
% Input:
% - speed (double): The initial speed estimate.
% - failures (double): Number of failures of the tracker.
% - skipping (integer): Number of skipped frames after each failure.
% - tracker (structure): A valid tracker descriptor.
% - sequence (structure): A valid sequence descriptor.
%
% Output:
% - normalized_speed (double): Normalized speed estimate.
% - actual_speed (double): Corrected raw speed based on supplied information,
%

if ~isfield(tracker, 'performance')
    error('Tracker %s has no performance profile, unable to normalize speed.', tracker.identifier);
end;

performance = tracker.performance;

factor = performance.nonlinear_native;
startup = 0;

if strcmpi(tracker.interpreter, 'matlab')
    if isfield(performance, 'matlab_startup')
        startup = performance.matlab_startup;
    else
        model = get_global_variable('matlab_startup_model', []);
		if ~isempty(model)
			startup = model(1) * performance.reading + model(2);
		end;
    end;
end

failure_count = cellfun(@(x) numel(x), failures, 'UniformOutput', true);

if tracker.trax
	actual_length = sequence.length - (skipping - 1) * failure_count;
	full_length = sequence.length;
	startup_time = startup * (1 + failure_count);
else
	full_length = cellfun(@(x) sum(sequence.length - x - (skipping - 1)), failures, 'UniformOutput', true) + sequence.length;
	actual_length = full_length;
	startup_time = startup * (1 + failure_count);
end;

actual_speed = (((speed .* full_length) - startup_time) ./ actual_length);
normalized_speed = actual_speed / factor;

author: votchallenge

repo: toolkit-legacy


text:
function [D] = biharmonic_distance(V,F,i,dim,p)
  % BIHARMONIC_DISTANCE Takes a mesh (V,F) and returns a distance field D from
  % all points in V to the ith vertex, according to the biharmonic embedding
  %
  % [D] = biharmonic_distance(V,F,i,dim)
  % [D] = biharmonic_distance(V,F,i,dim,p)
  %
  % Input:
  %   V  #V by dim list of vertex positions
  %   F  #F by 3 list of triangle indices
  %   i  index of vertex from which to calculate distances
  %   dim  requested dimension of the embedding
  %   Optional:
  %     p  exponent above eigen values
  %       0.5  "semi-harmonic" embedding
  %       1  commute time embedding, "harmonic"
  %       2  biharmonic {default}
  %       3  "triharmonic" embedding
  % Output:
  %   D  biharmonic distance field
  %

  % if index not given then use n/2th point
  if(~exist('i','var'))
    i = ceil(size(V,1)/2 + sqrt(size(V,1))/2);
  end

  % if dimension is not specfied use 4
  if(~exist('dim','var'))
    dim = 4;
  end

  % if power is not specfied use 2
  if(~exist('p','var'))
    p = 2;
  end

  B = biharmonic_embedding(V,F,dim,p);
  %B = biharmonic_embedding_yaron(V,F);
  D = sqrt(sum((repmat(B(i,:),size(B,1),1)-B(:,:)).^2,2));

  %tsurf(F,[V(:,1) V(:,2) D]);

end

author: alecjacobson

repo: gptoolbox


text:
%  Start a subsubsection
%
%  ::
%
%     document.subsubsection('option_name',option_value);
%
%  Args:
%
%     varargin: arguments need to come in pairs
%
%        - **title** : title of the subsubsection
%        - **numbering** : [{true}|false] whether to include the subsubsection in numbering of the table of contents
%
%

author: jmaih


text:
function P = principalComponents(X,q)
%principalComponents Principal components of a vector population.
%   P = principalComponents(X,Q) computes the principal-component
%   vectors of the vector population contained in the rows of X, a
%   matrix of size K-by-n where K (assumed to be > 1) is the number of
%   vectors and n is their dimensionality. Q, with values in the range
%   [0,n], is the number of eigenvectors used in constructing the
%   principal-components transformation matrix. P is a structure with
%   the following fields:
%
%     P.Y      K-by-Q matrix whose columns are the principal-
%              component vectors.
%     P.A      Q-by-n principal components transformation matrix
%              whose rows are the Q eigenvectors of Cx corresponding
%              to the Q largest eigenvalues.
%     P.X      K-by-n matrix whose rows are the vectors
%              reconstructed from the principal-component vectors.
%              P.X and P.Y are identical if Q = n.
%     P.ems    The mean square error incurred in using only the Q
%              eigenvectors corresponding to the largest
%              eigenvalues. P.ems is 0 if Q = n.
%     P.Cx     The n-by-n covariance matrix of the population in X.
%     P.mx     The n-by-1 mean vector of the population in X.
%     P.Cy     The Q-by-Q covariance matrix of the population in
%              Y. The main diagonal contains the eigenvalues (in
%              descending order) corresponding to the Q
%              eigenvectors.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

K = size(X,1);
X = double(X);

% Obtain the mean vector and covariance matrix of the vectors in X.
[P.Cx,P.mx] = covmatrix(X);
% Convert mean vector to a row vector.
P.mx = P.mx';

% Obtain the eigenvectors and corresponding eigenvalues of Cx.  The
% eigenvectors are the columns of n-by-n matrix V.  D is an n-by-n
% diagonal matrix whose elements along the main diagonal are the
% eigenvalues corresponding to the eigenvectors in V, so that X*V =
% D*V.
[V,D] = eig(P.Cx);

% Sort the eigenvalues in decreasing order. Rearrange the eigenvectors
% to match.
d = diag(D);
[d,idx] = sort(d);
d = flipud(d);
idx = flipud(idx);
D = diag(d);
V = V(:,idx);

% Now form the q rows of A from the first q columns of V.
P.A = V(:,1:q)';

% Compute the principal component vectors.
P.Y = P.A*(X - P.mx)'; % q-by-K matrix.

% Obtain the reconstructed vectors.
P.X = (P.A'*P.Y)' + P.mx;

% Convert P.Y to a K-by-q array and P.mx to an n-by-1 vector.
P.Y = P.Y';
P.mx = P.mx';

% The mean square error is the sum of all the eigenvalues minus the sum
% of the q largest eigenvalues.
d = diag(D);
P.ems = sum(d(q + 1:end));

% Covariance matrix of the Y's.
P.Cy = P.A*P.Cx*P.A';

author: dipum

repo: dipum-toolbox

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


text:
%  internal function

%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)

function save_nii_hdr(hdr, fid)

   if ~isequal(hdr.hk.sizeof_hdr,348),
      error('hdr.hk.sizeof_hdr must be 348.');
   end

   write_header(hdr, fid);

   return;					% save_nii_hdr

%---------------------------------------------------------------------
function write_header(hdr, fid)

        %  Original header structures
	%  struct dsr				/* dsr = hdr */
	%       {
	%       struct header_key hk;            /*   0 +  40       */
	%       struct image_dimension dime;     /*  40 + 108       */
	%       struct data_history hist;        /* 148 + 200       */
	%       };                               /* total= 348 bytes*/

   header_key(fid, hdr.hk);
   image_dimension(fid, hdr.dime);
   data_history(fid, hdr.hist);

   %  check the file size is 348 bytes
   %
   fbytes = ftell(fid);

   if ~isequal(fbytes,348),
      msg = sprintf('Header size is not 348 bytes.');
      warning(msg);
   end

   return;					% write_header

%---------------------------------------------------------------------
function header_key(fid, hk)

   fseek(fid,0,'bof');

	%  Original header structures
	%  struct header_key                      /* header key      */
	%       {                                /* off + size      */
	%       int sizeof_hdr                   /*  0 +  4         */
	%       char data_type[10];              /*  4 + 10         */
	%       char db_name[18];                /* 14 + 18         */
	%       int extents;                     /* 32 +  4         */
	%       short int session_error;         /* 36 +  2         */
	%       char regular;                    /* 38 +  1         */
	%       char dim_info;   % char hkey_un0;        /* 39 +  1 */
	%       };                               /* total=40 bytes  */

   fwrite(fid, hk.sizeof_hdr(1),    'int32');	% must be 348.

   % data_type = sprintf('%-10s',hk.data_type);	% ensure it is 10 chars from left
   % fwrite(fid, data_type(1:10), 'uchar');
   pad = zeros(1, 10-length(hk.data_type));
   hk.data_type = [hk.data_type  char(pad)];
   fwrite(fid, hk.data_type(1:10), 'uchar');

   % db_name   = sprintf('%-18s', hk.db_name);	% ensure it is 18 chars from left
   % fwrite(fid, db_name(1:18), 'uchar');
   pad = zeros(1, 18-length(hk.db_name));
   hk.db_name = [hk.db_name  char(pad)];
   fwrite(fid, hk.db_name(1:18), 'uchar');

   fwrite(fid, hk.extents(1),       'int32');
   fwrite(fid, hk.session_error(1), 'int16');
   fwrite(fid, hk.regular(1),       'uchar');	% might be uint8

   % fwrite(fid, hk.hkey_un0(1),    'uchar');
   % fwrite(fid, hk.hkey_un0(1),    'uint8');
   fwrite(fid, hk.dim_info(1),      'uchar');

   return;					% header_key

%---------------------------------------------------------------------
function image_dimension(fid, dime)

	%  Original header structures
	%  struct image_dimension
	%       {                                /* off + size      */
	%       short int dim[8];                /* 0 + 16          */
	%       float intent_p1;   % char vox_units[4];   /* 16 + 4       */
	%       float intent_p2;   % char cal_units[8];   /* 20 + 4       */
	%       float intent_p3;   % char cal_units[8];   /* 24 + 4       */
	%       short int intent_code;   % short int unused1;   /* 28 + 2 */
	%       short int datatype;              /* 30 + 2          */
	%       short int bitpix;                /* 32 + 2          */
	%       short int slice_start;   % short int dim_un0;   /* 34 + 2 */
	%       float pixdim[8];                 /* 36 + 32         */
	%			/*
	%				pixdim[] specifies the voxel dimensions:
	%				pixdim[1] - voxel width
	%				pixdim[2] - voxel height
	%				pixdim[3] - interslice distance
	%				pixdim[4] - volume timing, in msec
	%					..etc
	%			*/
	%       float vox_offset;                /* 68 + 4          */
	%       float scl_slope;   % float roi_scale;     /* 72 + 4 */
	%       float scl_inter;   % float funused1;      /* 76 + 4 */
	%       short slice_end;   % float funused2;      /* 80 + 2 */
	%       char slice_code;   % float funused2;      /* 82 + 1 */
	%       char xyzt_units;   % float funused2;      /* 83 + 1 */
	%       float cal_max;                   /* 84 + 4          */
	%       float cal_min;                   /* 88 + 4          */
	%       float slice_duration;   % int compressed; /* 92 + 4 */
	%       float toffset;   % int verified;          /* 96 + 4 */
	%       int glmax;                       /* 100 + 4         */
	%       int glmin;                       /* 104 + 4         */
	%       };                               /* total=108 bytes */

   fwrite(fid, dime.dim(1:8),        'int16');
   fwrite(fid, dime.intent_p1(1),  'float32');
   fwrite(fid, dime.intent_p2(1),  'float32');
   fwrite(fid, dime.intent_p3(1),  'float32');
   fwrite(fid, dime.intent_code(1),  'int16');
   fwrite(fid, dime.datatype(1),     'int16');
   fwrite(fid, dime.bitpix(1),       'int16');
   fwrite(fid, dime.slice_start(1),  'int16');
   fwrite(fid, dime.pixdim(1:8),   'float32');
   fwrite(fid, dime.vox_offset(1), 'float32');
   fwrite(fid, dime.scl_slope(1),  'float32');
   fwrite(fid, dime.scl_inter(1),  'float32');
   fwrite(fid, dime.slice_end(1),    'int16');
   fwrite(fid, dime.slice_code(1),   'uchar');
   fwrite(fid, dime.xyzt_units(1),   'uchar');
   fwrite(fid, dime.cal_max(1),    'float32');
   fwrite(fid, dime.cal_min(1),    'float32');
   fwrite(fid, dime.slice_duration(1), 'float32');
   fwrite(fid, dime.toffset(1),    'float32');
   fwrite(fid, dime.glmax(1),        'int32');
   fwrite(fid, dime.glmin(1),        'int32');

   return;					% image_dimension

%---------------------------------------------------------------------
function data_history(fid, hist)

	% Original header structures
	%struct data_history
	%       {                                /* off + size      */
	%       char descrip[80];                /* 0 + 80          */
	%       char aux_file[24];               /* 80 + 24         */
	%       short int qform_code;            /* 104 + 2         */
	%       short int sform_code;            /* 106 + 2         */
	%       float quatern_b;                 /* 108 + 4         */
	%       float quatern_c;                 /* 112 + 4         */
	%       float quatern_d;                 /* 116 + 4         */
	%       float qoffset_x;                 /* 120 + 4         */
	%       float qoffset_y;                 /* 124 + 4         */
	%       float qoffset_z;                 /* 128 + 4         */
	%       float srow_x[4];                 /* 132 + 16        */
	%       float srow_y[4];                 /* 148 + 16        */
	%       float srow_z[4];                 /* 164 + 16        */
	%       char intent_name[16];            /* 180 + 16        */
	%       char magic[4];   % int smin;     /* 196 + 4         */
	%       };                               /* total=200 bytes */

   % descrip     = sprintf('%-80s', hist.descrip);     % 80 chars from left
   % fwrite(fid, descrip(1:80),    'uchar');
   pad = zeros(1, 80-length(hist.descrip));
   hist.descrip = [hist.descrip  char(pad)];
   fwrite(fid, hist.descrip(1:80), 'uchar');

   % aux_file    = sprintf('%-24s', hist.aux_file);    % 24 chars from left
   % fwrite(fid, aux_file(1:24),   'uchar');
   pad = zeros(1, 24-length(hist.aux_file));
   hist.aux_file = [hist.aux_file  char(pad)];
   fwrite(fid, hist.aux_file(1:24), 'uchar');

   fwrite(fid, hist.qform_code,    'int16');
   fwrite(fid, hist.sform_code,    'int16');
   fwrite(fid, hist.quatern_b,   'float32');
   fwrite(fid, hist.quatern_c,   'float32');
   fwrite(fid, hist.quatern_d,   'float32');
   fwrite(fid, hist.qoffset_x,   'float32');
   fwrite(fid, hist.qoffset_y,   'float32');
   fwrite(fid, hist.qoffset_z,   'float32');
   fwrite(fid, hist.srow_x(1:4), 'float32');
   fwrite(fid, hist.srow_y(1:4), 'float32');
   fwrite(fid, hist.srow_z(1:4), 'float32');

   % intent_name = sprintf('%-16s', hist.intent_name);	% 16 chars from left
   % fwrite(fid, intent_name(1:16),    'uchar');
   pad = zeros(1, 16-length(hist.intent_name));
   hist.intent_name = [hist.intent_name  char(pad)];
   fwrite(fid, hist.intent_name(1:16), 'uchar');

   % magic	= sprintf('%-4s', hist.magic);		% 4 chars from left
   % fwrite(fid, magic(1:4),           'uchar');
   pad = zeros(1, 4-length(hist.magic));
   hist.magic = [hist.magic  char(pad)];
   fwrite(fid, hist.magic(1:4),        'uchar');

   return;					% data_history

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [mri] = align_ctf2acpc(mri, opt, template)

% ALIGN_CTF2ACPC performs an approximate rigid body alignment of the anatomical
% volume from CTF towards ACPC coordinates. Only the homogeneous transformation
% matrix is modified and the coordsys-field is updated.
%
% Use as
%   mri = align_ctf2acpc(mri)
%   mri = align_ctf2acpc(mri, opt)
%   mri = align_ctf2acpc(mri, opt, template)
%
% The first input argument is a FieldTrip MRI-structure, and the second optional
% argument specifies how the registration is to be done:
%   method = 0: only an approximate coregistration
%   method = 1: an approximate coregistration, followed by spm_affreg
%   method = 2: an approximate coregistration, followed by spm_normalise (default)
%
% When opt = 1 or 2, an optional template filename can be specified, which denotes
% the filename of the target volume. This is required when running in deployed
% mode.
%
% See also ALIGN_NEUROMAG2ACPC, ALIGN_FSAVERAGE2MNI

if nargin<2
  opt = 2;
end

%--------------------------------------------------------------------------
% do a first round of approximate coregistration using the predefined voxel
% locations of the fiducials and landmarks in the template T1 image from
% SPM

acpcvox2acpchead = [
  2     0     0   -92
  0     2     0  -128
  0     0     2   -74
  0     0     0     1
  ];

% these are the voxel indices of some points in the SPM canonical T1
acpcvox_Ac           = [46 64  37  1]';     % the anterior commissure
acpcvox_Ori          = [46 48  10  1]';     % approximately between the ears in T1.mnc
acpcvox_Nas          = [46 106 13  1]';     % approximately the nasion in T1.mnc
acpcvox_Lpa_canal    = [ 5 48  10  1]';     % Left ear canal
acpcvox_Rpa_canal    = [87 48  10  1]';     % Right ear canal

acpchead_Ac           = acpcvox2acpchead * acpcvox_Ac  ;
acpchead_Ori          = acpcvox2acpchead * acpcvox_Ori ;
acpchead_Nas          = acpcvox2acpchead * acpcvox_Nas ;
acpchead_Lpa_canal    = acpcvox2acpchead * acpcvox_Lpa_canal ;
acpchead_Rpa_canal    = acpcvox2acpchead * acpcvox_Rpa_canal ;

ctfvox2ctfhead   = mri.transform;
acpchead2ctfhead = ft_headcoordinates(acpchead_Nas(1:3), acpchead_Lpa_canal(1:3), acpchead_Rpa_canal(1:3), 'ctf');

%ctfvox2acpchead =  inv(acpchead2ctfhead) *  ctfvox2ctfhead;
ctfvox2acpchead =  acpchead2ctfhead \ ctfvox2ctfhead;

% change the transformation matrix, such that it returns approximate SPM head coordinates
mri.transform     = ctfvox2acpchead;
mri.vox2headOrig  = ctfvox2ctfhead;
mri.vox2head      = ctfvox2acpchead;
mri.head2headOrig = acpchead2ctfhead;
mri.coordsys      = 'acpc';

%--------------------------------------------------------------------------
% Do a second round of affine registration (rigid body) to get improved
% alignment with ACPC coordinate system. this is needed because there may be
% different conventions defining LPA and RPA. The affine registration may
% fail however, e.g. if the initial alignment is not close enough. In that
% case SPM will throw an error

if opt==1
  % use spm_affreg

  switch lower(spm('ver'))
    case 'spm2'
      if isdeployed
        if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using opt==1'); end
      else
        template = fullfile(spm('Dir'),'templates','T1.mnc');
      end

    case 'spm8'
      if isdeployed
        if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using opt==1'); end
      else
        template = fullfile(spm('Dir'),'templates','T1.nii');
      end

    case 'spm12'
      if isdeployed
        if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using opt==1'); end
      else
        template = fullfile(spm('Dir'),'toolbox','OldNorm','T1.nii');
        if ~exist('spm_affreg', 'file')
          addpath(fullfile(spm('Dir'),'toolbox','OldNorm'));
        end
      end
      fprintf('using ''OldNorm'' affine registration\n');

    otherwise
      ft_error('unsupported SPM version');
  end
  mri2 = ft_read_mri(template);

  tname1 = [tempname, '.img'];
  tname2 = [tempname, '.img'];
  V1 = ft_write_mri(tname1, mri.anatomy,  'transform', mri.transform,  'spmversion', spm('ver'), 'dataformat', 'nifti_spm');
  V2 = ft_write_mri(tname2, mri2.anatomy, 'transform', mri2.transform, 'spmversion', spm('ver'), 'dataformat', 'nifti_spm');

  % the below, using just spm_affreg does not work robustly enough in some cases
  flags.regtype = 'rigid';
  [M, scale]    = spm_affreg(V1,V2,flags);

  % some juggling around with the transformation matrices
  ctfvox2acpchead2  = M \ V1.mat;
  acpchead2ctfhead2 = ctfvox2ctfhead / ctfvox2acpchead2;

  % update the transformation matrix
  mri.transform     = ctfvox2acpchead2;

  % this one is unchanged
  mri.vox2headOrig  = ctfvox2ctfhead;

  % these are new
  mri.vox2head      = ctfvox2acpchead2;
  mri.head2headOrig = acpchead2ctfhead2;

  % delete the temporary files
  delete(tname1); delete(strrep(tname1, 'img', 'hdr'));
  delete(tname2); delete(strrep(tname2, 'img', 'hdr'));

elseif opt==2
  % use spm_normalise

  switch lower(spm('ver'))
    case 'spm2'
      if isdeployed
        if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using opt==2'); end
      else
        template = fullfile(spm('Dir'),'templates','T1.mnc');
      end

    case 'spm8'
      if isdeployed
        if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using opt==2'); end
      else
        template = fullfile(spm('Dir'),'templates','T1.nii');
      end

    case 'spm12'
      % this uses the 'OldNorm' functionality, so the path needs to be
      % added, can only be done if non-deployed.
      if isdeployed
        if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using opt==2'); end
      else
        template = fullfile(spm('Dir'),'toolbox','OldNorm','T1.nii');
        if ~exist('spm_normalise', 'file')
          addpath(fullfile(spm('Dir'),'toolbox','OldNorm'));
        end
      end
      fprintf('using ''OldNorm'' normalisation\n');

    otherwise
      ft_error('unsupported SPM version');
  end
  mri2 = ft_read_mri(template);

  tname1 = [tempname, '.img'];
  tname2 = [tempname, '.img'];
  V1 = ft_write_mri(tname1, mri.anatomy,  'transform', mri.transform,  'spmversion', spm('ver'), 'dataformat', 'nifti_spm');
  V2 = ft_write_mri(tname2, mri2.anatomy, 'transform', mri2.transform, 'spmversion', spm('ver'), 'dataformat', 'nifti_spm');

  flags.nits       = 0; %set number of non-linear iterations to zero
  flags.regtype    = 'rigid';
  params           = spm_normalise(V2,V1,[],[],[],flags);
  acpchead2ctfhead2 = acpchead2ctfhead*V1.mat*params.Affine/V2.mat;
  ctfvox2acpchead2  = acpchead2ctfhead2\ctfvox2ctfhead;

  % update the transformation matrix
  mri.transform     = ctfvox2acpchead2;

  % this one is unchanged
  mri.vox2headOrig  = ctfvox2ctfhead;

  % these are new
  mri.vox2head      = ctfvox2acpchead2;
  mri.head2headOrig = acpchead2ctfhead2;

  % delete the temporary files
  delete(tname1); delete(strrep(tname1, 'img', 'hdr'));
  delete(tname2); delete(strrep(tname2, 'img', 'hdr'));
end


text: %MEDIAN Datafile overload

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [mCatalog] = import_start(sFilterDir)

    % Default return value
    mCatalog = nan;

    % Check if filter directory ends with the appropriate directory separator
    sLastChar = sFilterDir(length(sFilterDir));
    if ~(sLastChar == filesep)
        sFilterDir = [sFilterDir filesep];
    end

    % Get directory in filter directory
    vrDir = dir(sFilterDir);
    if isempty(vrDir)
        errordlg('Directory does not exist');
        return;
    elseif length(vrDir) == 2
        errordlg('Directory is empty. No Filters found.');
        return;
    end

    % Initialize filter list
    cFilter = {};
    % Counter for detected filters
    nFilter = 1;

    for nCnt = 1:length(vrDir)
        % Check for directories
        if ~((strcmp(vrDir(nCnt).name, '.')) || (strcmp(vrDir(nCnt).name, '..')))
            % Get name of file
            sName = vrDir(nCnt).name;
            nPos = strfind(sName, '.m');
            % Is it a .m-file
            if ~isempty(nPos)
                sName = sName(1:(nPos - 1));
                try
                    % Is it an import filter
                    sDescription = feval(sName, 0);
                catch
                    sDescription = [];
                end
                % Is it really an import filter
                if ischar(sDescription)
                    % Insert it into filter list
                    cFilter{nFilter, 1} = sDescription;
                    cFilter{nFilter, 2} = sName;
                    nFilter = nFilter + 1;
                end
            end
        end
    end

    % Charmatrix with names of filters
    mFilterNames = char(cFilter{:,1});
    mFilterFiles = char(cFilter{:,2});

    % Invoke the user interface
    hDialog = import_dialog(mFilterNames, mFilterFiles, sFilterDir);

    % Analyze Output
    if ~ishandle(hDialog)
        answer = 0;
    else
        handles = guidata(hDialog);
        answer = handles.answer;
        % OK pressed
        if answer == 1
            % Get the values from figure
            nFilter = get(handles.lstFilter, 'Value');
            % Remove figure from memory
            delete(hDialog);
            % Get the catalog file
            [newfile, newpath] = uigetfile('*.*', 'Choose catalog datafile');
            % Cancel pressed?
            if ~(isequal(newfile, 0) || isequal(newpath, 0))
                newfile = [newpath newfile];
                % Everything ok?
                if length(newfile) > 1
                    mCatalog = feval(cFilter{nFilter,2}, 1, newfile);
                end
            end
        else
            % Remove figure from memory
            delete(hDialog);
        end
    end

author: CelsoReyes


%% REW_gym
% %{
imfolder = 'images\REW_gym';
im_n = 2;
imfile = cell(im_n,1);
imfile{1} = [imfolder '\' 'gym_01.jpg'];
imfile{2} = [imfolder '\' 'gym_02.jpg'];

im = cell(im_n,1);
for ii = 1:im_n
    im{ii} = imread(imfile{ii});
end

edge_list = [1,2];

% imsize = zeros(im_n,3);
%
% for ii = 1:im_n
%     imsize(ii,:) = size(im{ii});
%     if imsize(ii,1) > 720
%         scale = 720/size(im{ii}, 1);
%         im{ii} = imresize(im{ii}, scale);
%
%         imsize(ii,:) = size(im{ii});
%     end
% end

mosaic = REW_mosaic( im, edge_list, 1, 'persp', 0, imfolder );
%}


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

%%begin

% Forward dynamics is the computation of joint accelerations given position and
% velocity state, and actuator torques.  It is useful in simulation of a robot
% control system.
%
% Consider a Puma 560 at rest in the zero angle pose, with zero applied joint
% torques. The joint acceleration would be given by

mdl_puma560
p560.accel(qz, zeros(1,6), zeros(1,6))

% To be useful for simulation this function must be integrated.  fdyn() uses the
% MATLAB function ode45() to integrate the joint acceleration.  It also allows
% for a user written function to compute the joint torque as a function of
% manipulator state.
%
% To simulate the motion of the Puma 560 from rest in the zero angle pose
% with zero applied joint torques

tic
[t q qd] = p560.nofriction().fdyn(10, [], qz);
toc

% and the resulting motion can be plotted versus time

subplot(3,1,1)
plot(t,q(:,1))
xlabel('Time (s)');
ylabel('Joint 1 (rad)')
subplot(3,1,2)
plot(t,q(:,2))
xlabel('Time (s)');
ylabel('Joint 2 (rad)')
subplot(3,1,3)
plot(t,q(:,3))
xlabel('Time (s)');
ylabel('Joint 3 (rad)')

% Clearly the robot is collapsing under gravity, but it is interesting to
% note that rotational velocity of the upper and lower arm are exerting
% centripetal and Coriolis torques on the waist joint, causing it to rotate.

% This can be shown in animation also
clf
p560.plot(q)

author: Allopart

repo: rbpf-gmapping


text:
% Imports a ROS bag cotaining event data recorded with the rpg_dvs_ros
% driver into the same workspace structure as ImportAedat uses.
%
% Code taken from ConvertRosbagToAedat, written by Guillermo Gallego
%
% RPG DVS/DAVIS driver for ROS:
% https://github.com/uzh-rpg/rpg_dvs_ros
% Currently, the code assumes that the ROS data is published/stored in the
% dafault topic names: /dvs/events, /dvs/image_raw and /dvs/imu
%
% To read ROS bags, you need to have in the matlab path the library
% https://github.com/bcharrow/matlab_rosbag/releases
%
% AEDAT 2.0 format is explained in:
% https://inilabs.com/support/software/fileformat/#h.4ydb2xpu03ik
% Timestamps tick is 1 microsecond
%
% Input:
% aedat struct, containing at least field 'importParams',
% containing at least field 'filePath': full name and path of the ROS bag (dataset) to read
%
% Output:
% aedat: data structure used in AedatTools with information and
% content of the input dataset.
%
% Guillermo Gallego
% e-mail: guillermo.gallego@ifi.uzh.ch

function aedat = ImportRosbag(aedat)

if ~isfield(aedat, 'importParams') || ~isfield(aedat.importParams, 'filePath')
    error('No filename')
end

% Load a bag and get information about it
bag = ros.Bag.load(aedat.importParams.filePath);
bag.info()

% Start to fill in aedat structure, for conversion
aedat.data = [];

%% Default ROS topics with the published data
topic_events = '/dvs/events';
topic_images = '/dvs/image_raw';
topic_imu = '/dvs/imu';

%% Get first timestamp of the bag so that it is subtracted from all timestamps
% and hence the 32 bits of the timestamps are better used
bag.resetView({topic_events, topic_images, topic_imu});
msg = bag.read();
% IMU or image messages have the timestamp in the header
% Event messages have header timestamp empty
if isfield(msg, 'events')
    % first message is a packet of events
    time_first_sec_uint32 = uint32( msg.events(3,1) );
else
    % first message is image or IMU
    time_first_sec_uint32 = msg.header.stamp.sec;
end

margin = 10; % some margin. Maybe some events happen before this time
if time_first_sec_uint32 > margin
    time_first_sec_uint32 = time_first_sec_uint32 - margin;
end
bag.resetView({topic_events, topic_images, topic_imu});

%% Frames (i.e., images)
disp('Reading image topic from ROS bag');
msgs = bag.readAll(topic_images);
num_msgs = length(msgs);
disp(['Found ' num2str(num_msgs) ' frames']);
if (num_msgs > 0)
    aedat.data.frame = []; % create images struct
    aedat.data.frame.numEvents = num_msgs; % number of frames
    disp('Converting APS-frame messages...');
    img_width = double(msgs{1}.width);
    img_height = double(msgs{1}.height);
    v_zero = zeros(num_msgs,1);
    v_one = ones(num_msgs,1);
    aedat.data.frame.xPosition = uint16(v_zero);
    aedat.data.frame.yPosition = uint16(v_zero);
    aedat.data.frame.xLength = uint16(img_width * v_one); % 240
    aedat.data.frame.yLength = uint16(img_height* v_one); % 180
    for ii = 1:num_msgs
        % timestamp: remove offset, convert to microseconds and cast to uint32
        time_us = double(msgs{ii}.header.stamp.sec - time_first_sec_uint32) * 1e6 ...
            + double(msgs{ii}.header.stamp.nsec) / 1e3; % microseconds
        aedat.data.frame.timeStampStart(ii) = uint32(time_us);
        aedat.data.frame.timeStampEnd(ii) = uint32(time_us);

        img = reshape(msgs{ii}.data,img_width,img_height)';
        img = double(img) * 4; % convert from 8 bits to 10 bits: [0,255] to [0,1023]
        img = flipud(img);
        if (img_width == 346)
            % miniDAVIS346 needs to flip left-right events and frames
            img = fliplr(img);
        end
        aedat.data.frame.samples{ii,1} = img; % cell array with the grayscale values
    end
    % column vectors
    aedat.data.frame.timeStampStart = aedat.data.frame.timeStampStart';
    aedat.data.frame.timeStampEnd = aedat.data.frame.timeStampEnd';
    disp('...done!');
end

%% IMU
disp('Reading IMU topic from ROS bag');
msgs = bag.readAll(topic_imu);
num_msgs = length(msgs);
disp(['Found ' num2str(num_msgs) ' IMU samples']);
if (num_msgs > 0)
    aedat.data.imu6 = []; % create IMU struct
    aedat.data.imu6.numEvents = num_msgs; % number of measurements
    disp('Converting IMU messages...');
    % Column vectors
    v_zero = zeros(num_msgs,1);
    aedat.data.imu6.timeStamp = uint32(v_zero);
    aedat.data.imu6.accelX = v_zero;
    aedat.data.imu6.accelY = v_zero;
    aedat.data.imu6.accelZ = v_zero;
    aedat.data.imu6.gyroX = v_zero;
    aedat.data.imu6.gyroY = v_zero;
    aedat.data.imu6.gyroZ = v_zero;
    aedat.data.imu6.temperature = 28 + v_zero; % centigrades
    for ii = 1:num_msgs
        % timestamp: remove offset, convert to microseconds and cast to uint32
        time_us = double(msgs{ii}.header.stamp.sec - time_first_sec_uint32) * 1e6 ...
            + double(msgs{ii}.header.stamp.nsec) / 1e3; % microseconds
        aedat.data.imu6.timeStamp(ii) = uint32(time_us);

        aedat.data.imu6.accelX(ii) = msgs{ii}.linear_acceleration(1);
        aedat.data.imu6.accelY(ii) = msgs{ii}.linear_acceleration(2);
        aedat.data.imu6.accelZ(ii) = msgs{ii}.linear_acceleration(3);
        aedat.data.imu6.gyroX(ii) = msgs{ii}.angular_velocity(1);
        aedat.data.imu6.gyroY(ii) = msgs{ii}.angular_velocity(2);
        aedat.data.imu6.gyroZ(ii) = msgs{ii}.angular_velocity(3);
        % aedat.data.imu6.temperature(ii) = 0;
    end

    % Conversion of units.
    % Acceleration is imported as m/s^2, and aedat wants it in multiples of g
    gravity = 9.8; % m/s^2
    aedat.data.imu6.accelX = aedat.data.imu6.accelX / gravity;
    aedat.data.imu6.accelY = aedat.data.imu6.accelY / gravity;
    aedat.data.imu6.accelZ = aedat.data.imu6.accelZ / gravity;
    % Angular velocity is imported as rad/s, and aedat expects deg/s
    rad_in_deg = 180 / pi; % degrees in one radian
    aedat.data.imu6.gyroX = aedat.data.imu6.gyroX * rad_in_deg;
    aedat.data.imu6.gyroY = aedat.data.imu6.gyroY * rad_in_deg;
    aedat.data.imu6.gyroZ = aedat.data.imu6.gyroZ * rad_in_deg;

    disp('...done!');
end

%% Events
disp('Reading event topic from ROS bag');
msgs = bag.readAll(topic_events);
num_msgs = length(msgs);
disp(['Found ' num2str(num_msgs) ' event messages']);
if (num_msgs > 0)
    aedat.data.polarity = []; % create events struct
    disp('Converting event messages...');

    % Get size of the event sensor
    ev_sensor_width = double(msgs{1}.width);
    ev_sensor_height = double(msgs{1}.height);

    % Compute number of events
    numEvents = 0;
    for ii = 1:num_msgs
        numEvents = numEvents + size(msgs{ii}.events,2);
    end
    disp(['Number of events = ' num2str(numEvents)]);
    aedat.data.polarity.numEvents = numEvents;

    % Allocate memory for all events. Read events
    idx_firt_ev_msg = 1;
    for ii = 1:num_msgs
        if mod(ii,100) == 0
            disp(['Message ' num2str(ii) ' of ' num2str(num_msgs)]);
        end
        num_ev_msg = size(msgs{ii}.events,2);
        idx = idx_firt_ev_msg + (0:num_ev_msg-1)';
        idx_firt_ev_msg = idx_firt_ev_msg + num_ev_msg;

        % Read (x, y, t, polarity) of the events

        % timestamp: remove offset, convert to microseconds and cast to uint32
        time_us = (msgs{ii}.events(3,:) - double(time_first_sec_uint32)) * 1e6; % microseconds
        aedat.data.polarity.timeStamp(idx) = uint32(time_us);

        if (ev_sensor_width == 128) || (ev_sensor_width == 346)
            % DVS and miniDAVIS346 need to flip events. Don't ask why
            aedat.data.polarity.x(idx) = uint16(ev_sensor_width-1 -msgs{ii}.events(1,:));
        else
            aedat.data.polarity.x(idx) = uint16(msgs{ii}.events(1,:));
        end
        aedat.data.polarity.y(idx) = uint16(ev_sensor_height-1 -msgs{ii}.events(2,:));
        aedat.data.polarity.polarity(idx) = logical(msgs{ii}.events(4,:));
    end
    % column vectors
    aedat.data.polarity.timeStamp = aedat.data.polarity.timeStamp';
    aedat.data.polarity.x = aedat.data.polarity.x';
    aedat.data.polarity.y = aedat.data.polarity.y';
    aedat.data.polarity.polarity = aedat.data.polarity.polarity';
    disp('...done!');
end

%% Fill in additional information for aedat file
aedat.info = [];

% Set sensor size
if exist('ev_sensor_width','var') && exist('ev_sensor_height','var')
    % Get the sensor size from the event topic
    aedat.info.deviceAddressSpace = [ev_sensor_width, ev_sensor_height];
elseif exist('img_width','var') && exist('img_height','var')
    % Get the sensor size from the image topic
    aedat.info.deviceAddressSpace = [img_width, img_height];
else
    error('Could not determine the sensor size from the events or images');
end

% Set the type of sensor used, according to resolution
% See file BasicSourceName.m
sensor_size = aedat.info.deviceAddressSpace;
if (sensor_size(1) == 240) && (sensor_size(2) == 180)
    aedat.info.source = 'Davis240C';
elseif (sensor_size(1) == 346) && (sensor_size(2) == 260)
    aedat.info.source = 'Davis346BMono';
elseif (sensor_size(1) == 128) && (sensor_size(2) == 128)
    aedat.info.source = 'Dvs128';
else
    error(['Unknown sensor type for resolution: ' ...
        num2str(sensor_size(1)) ' x ' num2str(sensor_size(2)) ' pixels']);
end

% Fill in standard info
aedat = NumEventsByType(aedat);
aedat = FindFirstAndLastTimeStamps(aedat);

author: panpanfei


text:
%UIDrawSubCodes UIDrawSubCodes enumeration for the brick
%
% Notes::
% - UIDrawSubCodes can be found in the EV3 documentation and source code
% (bytecodes.h)
% - USE SCREEN FOR GRAPHICS
%   -----------------------
%
%     opUI_DRAW,LC0(FILLWINDOW),LC0(0),LC0(0),LC0(0)
%
%       (clear entire screen by filling it with zero - background color)
%
%
%     opUI_DRAW, - - - - - -
%     opUI_DRAW, - - - - - -
%     opUI_DRAW, - - - - - -
%
%       (draw graphical stuff on the screen)
%
%
%     opUI_DRAW,LC0(UPDATE)
%
%       (show the stuff by updating the screen)
%
%
%
%   RESTORE RUN SCREEN
%   ------------------
%
%     opUI_WRITE,LC0(INIT_RUN)
%
%       (enable the animated run screen)

classdef UIDrawSubCodes < uint8
    enumeration
         Update(0)
         Clean (1)
         Pixel (2)
         Line (3)
         Circle (4)
         Text (5)
         Icon (6)
         Picture (7)
         Value (8)
         Fillrect(9)
         Rect (10)
         Notification (11)
         Question (12)
         Keyboard (13)
         Browse (14)
         Vertbar (15)
         Inverserect (16)
         SelectFont (17)
         Topline (18)
         Fillwindow (19)
         Scroll (20)
         Dotline (21)
         ViewValue (22)
         ViewUnit (23)
         Fillcircle (24)
         Store (25)
         Restore (26)
         IconQuestion (27)
         Bmppfile(28)
         Popup (29)
         GraphSetup (30)
         GraphDraw (31)
         Textbox (32)
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function minmax(ob)

minmax(struct(ob))

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%  Sorts population pf mcmc draws
%
%  ::
%
%    [pop,tags]=sort_population(pop)
%
%  Args:
%
%     - **pop** [struct]: vector of individuals, each with fields
%
%       - **f** [numeric]: fitness level
%       - **x** [vector]: parameters
%
%  Returns:
%     :
%
%     - **pop** [struct]: sorted vector of individuals
%     - **tags** [vector]: reordering index
%
%


text:
function [tuples,fStar,qStar,u,exitCode]=assign3D(C,maximize,subgradMethod,subgradParams,maxIter,AbsTol,RelTol)
%%ASSIGN3D Approximate the solution to the operations research axial 3D
%         assignment problem using a dual-primal Lagrangian relaxation
%         technique. Such problems are NP-hard. The optimization problem
%         being solved is
%         minimize (or maximize)
%         \sum_{i=1}^{n1}\sum_{j=1}^{n2}\sum_{k=1}^{n3}C_{i,j,k}*\rho_{i,j,k}
%         subject to
%         \sum_{i=1}^{n1}\sum_{j=1}^{n2}\rho_{i,j,k}<=1 for all k
%         \sum_{i=1}^{n1}\sum_{k=1}^{n3}\rho_{i,j,k}<=1 for all j
%         \sum_{j=1}^{n2}\sum_{k=1}^{n3}\rho_{i,j,k} =1 for all i
%         \rho_{i,j,k} = 0 or 1
%         assuming that n1<=n2<=n3, and C is an n1Xn2Xn3 cost matrix.
%         This is equivalent to the optimization problem
%         min (or max) sum_{i} C(i,phi_2(i),phi_3(i))
%         where phi_2 and phi_3 are length n1 arrangements of n2 and n3
%         items over which the minimization is performed.
%
%INPUTS: C An n1Xn2Xn3 cost hypermatrix with n1<=n2<=n3. C cannot contain
%          any NaNs and the largest finite element minus the smallest
%          element is a finite quantity (does not overflow) when performing
%          minimization and where the smallest finite element minus the
%          largest element is finite when performing maximization.
%          Forbidden assignments can be given costs of +Inf for
%          minimization and -Inf for maximization. During minimization,
%          there should be no elements with -Inf cost and no elements with
%          +Inf cost during maximization.
% maximize If true, the minimization problem is transformed into a
%          maximization problem. The default if this parameter is omitted
%          or an empty matrix is passed is false.
% subgradMethod A parameter indicating the subgradient optimization
%          algorithm to use. Possible values are:
%          0 Polyak's Method from [2]. (The default if omitted or an empty
%            matrix is passed) In this method, for minimization, the dual
%            variables are updated as uNew=u+gamma*(fStar-q)/norm(g)^2*g,
%            where g is the subgradient, q is the current dual cost, and g
%            is the subgradient vector. gamma is a design parameter.
%          1 Bragin's Method from [3]. The dual update for minimization is
%            uNew=u+alpha*g
%            where during the first step, alpha=(fStar-q)/norm(g)^2 and for
%            subsequent steps it is
%            alpha=(1-1/(M*k^(1-1/k^r)))*alphaPrev*norm(gPrev)/norm(g)
%            where alphaPrev is alpha for the previous step, k is the step
%            number (starting at 0), gPrev is the subgradient before the
%            current step and M and r are design parameters.
%          2 Bertsekas' Heuristic Method from [4]. The update is
%            uNew=u+alpha*g
%            where alpha=((1+a/beta^b)*qMax-q)/norm(g)^2 where qMax is the
%            highest dual cost value encountered thus far (when
%            minimizing), a and b are design parameters, and beta is a
%            value that increases or decreases in the algorithm depending
%            on whether or not there was an improvement in the best dual
%            cost found after the last step.
%          3 Shor's Space Dilation Algorithm from [5].
%            The update is
%            uNew=u+alpha*H*g;
%            where alpha=(2*M/(M+1))*(fStar-q)/norm(d), d=g and H is a
%            matrix that also depends on M and d, which is a design
%            parameter. Also, after NR iterations, the recursive turms that
%            go into H reset.
%          4 The r Space Dilation Algorithm from [5]. This is the same as 3
%            except d=g-gPrev except for the first iteration and if
%            g=gPrev, in which case d=g.
% subgradParams This is a structure that takes the design parameters for
%          the selected cubgradient algorithm. Possible fields of the
%          structure as well as default values depend on the selected
%          subgradient method and are:
%          Method 0: 'gamma' with default 1.
%          Method 1: 'M' and 'r' with defaults 2.8 and 0.06.
%          Method 2: 'a and 'b' with defaults 0.3 and 1.5.
%          Method 3,4: 'M', 'NR', and 'normBound' with defaults 2, n3, and
%                     eps(). normBound is such that if norm(B.'*d) in the
%                     computation of the matrix H is <= normBound, then B
%                     is reset to the identity matrix.
%  maxIter The maximum number of iterations to perform. The default value
%          if this parameter is omitted or an empty matrix is passed is 20.
%   AbsTol The absolute duality gap to use for convergence determiniation.
%          Convergence is declared if (fStar-qStar)<=AbsTol, where if
%          minimizing, fStar is the smallest value of the primal function
%          found and qStar is the highest value of the dual cost function
%          found. The default if omitted or an empty matrix is passed is
%          1e-10.
%   RelTol The relative duality gap to use for convergence determiniation.
%          Convergence is declared if (fStar-qStar)<=RelTol*abs(qStar). The
%          default if omitted or an empty matrix is passed is 0.05.
%
%OUTPUTS: tuples A 3Xn1 matrix where tuples(:,i) is the ith assigned tuple
%                in C as described above. An empty matrix is returned if
%                the problem is infeasible or the heuristic to obtain a
%                feasible solution failed.
%          fStar The cost of the assignment in tuples. This is the sum of
%                the values of the assigned elements in C.
%          qStar The value of the best dual solution found. The absolute
%                duality gap is fStar-qStar and the relative duality gap is
%                abs(fStar-qStar)/abs(qStar)
%              u The n3X1 vector of dual variables.
%       exitCode A parameter indicating how the algorithm terminated.
%                Possible values are:
%                -3 A non-finite number arose in the dual variables, so the
%                   algorithm stopped.
%                -2 The algorithm did not converge within the alotted
%                   number of iterations. However, a valid feasible
%                   assignment was found.
%                -1 A subproblem in computing the dual or in obtaining a
%                   feasible primal solution was infeasible.The output
%                   tuples will be an empty matrix.
%                >=0 Values that are zero or positive indicate the
%                   convergence was obtained and the returned value is the
%                   number of iterations.
%
%This function implements the algorithm of [1], but modified so that it
%does not have any unconstrained indices and offering different subgradient
%methods. We are solving the "operations research" 3D assignment problem
%rather than the "data fusion" 3D assignment problem of [1].
%
%EXAMPLE:
%Here is an example using a random problem:
% n1=10;
% n2=12;
% n3=15;
% C=randn(n1,n2,n3);
% maximize=false;
% [tuples,fStar,qStar,u,exitCode]=assign3D(C,maximize)
%
%REFERENCES:
%[1] K. Pattipati, S. Deb, Y. Bar-Shalom, and R. B. Washburn Jr., "A
%    new relaxation algorithm and passive sensor data association," IEEE
%    Transactions on Automatic Control, vol. 37, no. 2, pp. 198-213, Feb.
%    1992.
%[2] B. T. Polyak, "Minimization of unsmooth functionals," USSR
%    Computational Mathematics and Mathematical Physics, vol. 9, no. 3, pp.
%    14-29, 1969.
%[3] M. A. Bragin, P. B. Luh, J. H. Yan, N. Yu, and G. A. Stern,
%    "Convergence of the surrogate Lagrangian relaxation method," Journal
%    of Optimization Theory and Applications, vol. 164, no. 1, pp. 173-201,
%    Jan. 2015.
%[4] D. P. Bertsekas, Nonlinear Programming, 3rd ed. Belmont, MA: Athena
%    Scientific, 2016, Chapter 7.5.
%[5] N. Z. Shor, "Utilization of the operation of space dilation in the
%    minimization of convex functions," Cybernetics, vol. 6, no. 1, pp. 7-
%    15, Dec. 1972.
%
%February 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

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

if(nargin<3||isempty(subgradMethod))
    subgradMethod=0;
end

if(nargin<4)
    subgradParams=[];
end

if(nargin<5||isempty(maxIter))
    maxIter=20;
end

if(nargin<6||isempty(AbsTol))
    AbsTol=1e-10;
end

if(nargin<7||isempty(RelTol))
    RelTol=0.05;
end

nVals=size(C);

%Deal with the special case of C being scalar.
if(isscalar(C))
    tuples=[1;1;1];
    fStar=C(1);
    qStar=fStar;
    u=0;
    exitCode=0;
    return;
end

if(isempty(C))%The empty matrix special case.
    tuples=[];
    fStar=0;
    qStar=0;
    u=[];
    exitCode=0;
    return;
end

if(maximize)
    C=-C;
end

if(length(nVals)<3||~(nVals(1)<=nVals(2)&&nVals(2)<=nVals(3)))
    error('It is required that size(C,1)<=size(C,2)<=size(C,3)')
end

n1=nVals(1);
n2=nVals(2);
n3=nVals(3);

if((n1==n2)&&(n2==n3))
    unconstU=true;
else
    unconstU=false;
end

%Get the parameters and initial values for the chosen subgradient method.
switch(subgradMethod)
    case 0%Polyak's method
        gammaParam=1;

        if(nargin>3&&~isempty(subgradParams))
            if(isfield(subgradParams,'gamma'))
                gammaParam=subgradParams.gamma;
            end
        end
    case 1%Bragin's Method
        M=2.8;
        r=0.06;
        gNormPrev=[];
        alphaPrev=[];

        if(nargin>3&&~isempty(subgradParams))
            if(isfield(subgradParams,'M'))
                M=subgradParams.M;
            end

            if(isfield(subgradParams,'r'))
                r=subgradParams.r;
            end
        end
    case 2%Bertsekas' Heuristic Method
        a=0.3;
        b=1.5;
        beta=1;%Initialize

        if(nargin>3&&~isempty(subgradParams))
            if(isfield(subgradParams,'a'))
                a=subgradParams.a;
            end

            if(isfield(subgradParams,'b'))
                b=subgradParams.b;
            end
        end
    case 3%Shor's Space Dilation Algorithm.
        NR=n3;
        M=2;
        rho=(M-1)/(M+1);
        normBound=eps();

        B=eye(n3,n3);%Allocate and initialize
        dPrev=[];

        if(nargin>3&&~isempty(subgradParams))
            if(isfield(subgradParams,'NR'))
                NR=subgradParams.NR;
            end

            if(isfield(subgradParams,'M'))
                M=subgradParams.M;
            end

            if(isfield(subgradParams,'normBound'))
                normBound=subgradParams.normBound;
            end
        end
    case 4%The r Space Dilation Algorithm.
        NR=n3;
        M=2;
        rho=(M-1)/(M+1);
        normBound=eps();

        B=eye(n3,n3);%Allocate and initialize
        gPrev=[];
        dPrev=[];

        if(nargin>3&&~isempty(subgradParams))
            if(isfield(subgradParams,'NR'))
                NR=subgradParams.NR;
            end

            if(isfield(subgradParams,'M'))
                M=subgradParams.M;
            end

            if(isfield(subgradParams,'normBound'))
                normBound=subgradParams.normBound;
            end
        end
    otherwise
        error('Invalid subgradient method specified.')
end

u=zeros(n3,1);%Allocate and initialize the dual variables.
%Allocate space for the return tuples.
tuples=zeros(3,n1);

%qStar is the best (maximum) dual cost encountered thus far.
qStar=-Inf;

%fStar is the best (minimum) feasible primal cost encountered
%thus far.
fStar=Inf;
%If the loop exits without convergence, then exitCode is not modified.
exitCode=-2;
for k=0:(maxIter-1)
%%%%%%
%DUAL COST AND SUBGRADIENT UPDATE. SEE FIG. 3 IN [1].
%%%%%%
    %Note that d3=C;
    [d2,gamma2]=min(bsxfun(@plus,C,reshape(u,[1,1,n3])),[],3);

    %gamma1 is the columns for rows (nonzero elements in omega).
    [gamma1, ~, minVal]=assign2D(d2,false);

    %If the subproblem is infeasible.
    if(isempty(gamma1))
        tuples=[];
        fStar=[];
        qStar=[];
        u=[];
        exitCode=-1;
        return
    end

    q=minVal-sum(u);%The dual cost.

    %Keep track of the maximum q value. This is used for testing
    %convergence.
    if(q>qStar)
        qStar=q;

        %The duality gap.
        costGap=fStar-qStar;

        %The correctness of this on the first iterations requires proper
        %handling of NaNs.
        if(costGap<=AbsTol||(costGap<abs(qStar)*RelTol))
            exitCode=k+1;%The algorithm converged.
            break;
        end
    end

    %Compute the subgradient
    g=-1*ones(n3,1);
    for i1=1:n1
        i2=gamma1(i1);
        i3=gamma2(i1,i2);

        g(i3)=g(i3)+1;
    end

%%%%%%
%TEST THE DUAL VARIABLES AND FINISH THE LOOP. See Fig. 2 in [1].
%%%%%%
    gNorm2=dot(g,g);

    %If no constraints are violated, then we have a feasible solution and
    %the algorithm has converged to a local or global minimum point.
    if(gNorm2==0)
        %If the dual cost is less than the best primal solution found with
        %a heuristic, then we use the tuples for that solution. Otherwise,
        %the best primal solution to this point is returned.
        if(q<fStar)
            fStar=q;
            %Record the tuples.
            for i1=1:n1
                i2=gamma1(i1);
                i3=gamma2(i1,i2);
                tuples(:,i1)=[i1;i2;i3];
            end
        end

        exitCode=k+1;%The algorithm converged.
        break;
    end

    [tuples,fStar,retVal]=updateBestFeasSol3D(C,gamma1,tuples,fStar,qStar,AbsTol,RelTol);

    %If the primal converged.
    if(retVal==0)
        exitCode=k+1;
        break;
    end

    %If an error occurred obtaining a feasible solution.
    if(retVal<0)
        tuples=[];
        fStar=[];
        qStar=[];
        u=[];
        exitCode=retVal;
        return;
    end

    %Perform the subgradient update.
    switch(subgradMethod)
        case 0%Polyak's method
            alpha=gammaParam*((fStar-q)/gNorm2);
            u=u+alpha*g;
        case 1%Bragin's Method
            gNorm=sqrt(gNorm2);
            if(k==0)
                alpha=((fStar-q)/gNorm2);
            else
                alpha=(1-1/(M*k^(1-k^(-r))))*alphaPrev*(gNormPrev/gNorm);
            end
            u=u+alpha*g;

            alphaPrev=alpha;
            gNormPrev=gNorm;
        case 2%Bertsekas' Heuristic Method
            %Note that we do not have to keep track of qStar value from the
            %previous iteration. Though in [1], the comparison is with the
            %previous qStar, if q>qStarPrev, then qStar will be set to q,
            %so the comparison q<qStarPrev has the same result as q<qStar
            %and if q<=qStarPrev, then qStarPrev=qStar, so the comparison
            %is still correct.
            if(q<qStar)
                beta=beta+1;
            else
                beta=max(beta-1,1);
            end

            alpha=abs(((1+a)/(beta^b))*qStar-q)/gNorm2;

            u=u+alpha*g;
        case 3%Shor's Space Dilation Algorithm.
            d=g;

            %The H matrix is initialized on the first iteration and reset
            %to the identity matrix every NR iterations.
            if(mod(k,NR)==0)
                B=eye(n3,n3);
                normVal=norm(d);%=norm(B'*d)
            else
                zeta=B.'*dPrev/norm(B.'*dPrev);
                R=eye(n3,n3)+(rho-1)*(zeta*zeta.');
                B=B*R;
                normVal=norm(B.'*d);
                if(normVal<=normBound)
                    B=eye(n3,n3);
                    normVal=norm(d);
                end
            end
            H=(B*B.')/normVal;

            alpha=(2*M/(M+1))*((fStar-q)/norm(d));

            u=u+alpha*H*d;

            dPrev=d;
        case 4%The r Space Dilation Algorithm.
            if(k==0||all(g==gPrev))
                d=g;
            else
                d=g-gPrev;
            end

            %The H matrix is initialized on the first iteration and reset
            %to the identity matrix every NR iterations.
            if(mod(k,NR)==0)
                B=eye(n3,n3);
                normVal=norm(d);%=norm(B'*d)
            else
                zeta=B.'*dPrev/norm(B.'*dPrev);
                R=eye(n3,n3)+(rho-1)*(zeta*zeta.');
                B=B*R;
                normVal=norm(B.'*d);
                if(norm(B.'*d)<=normBound)
                    B=eye(n3,n3);
                    normVal=norm(d);
                end
            end
            H=(B*B.')/normVal;

            alpha=(2*M/(M+1))*((fStar-q)/norm(d));

            u=u+alpha*H*d;

            gPrev=g;
            dPrev=d;
        otherwise
            error('Invalid subgradient method specified.')
    end

    if(any(~isfinite(u)))
        %This can sometimes occur with big problems and poor stepsizes.
        exitCode=-3;
        break;
    end

    %Unless the constraints are equality constraints, the dual variables
    %have to be clipped.
    if(unconstU==false)
        u(u<0)=0;
    end
end

%Adjust the gains for the case where the initial cost matrix is transformed
%so that maximization can be performed.
if(maximize)
    fStar=-fStar;
    qStar=-qStar;
    u=-u;
end
end

function [optTuples,fStar,exitCode]=updateBestFeasSol3D(C,gamma1,optTuples,fStar,qStar,AbsTol,RelTol)
%%UPDATEBESTFEASSOL3D This function implements a subroutine to
%         heuristically obtain a feasible solution given a partial set of
%         assignments. Here, it is just for the 3D assignment problem.
%         Additionally, this function updates the best feasible solution
%         cost fStar and its associated set of tuples optTuples.
%
%INPUTS: C The original n1Xn2Xn3 3D cost matrix.
%   gamma1 gamma1(i1) gives the value of n2 for each value of i1.
% optTuples The set of tuples associated with the current best primal cost.
%    fStar The current best (lowest) primal cost function value found.
%    qStar The current best (highest) dual cost function value found.
%   AbsTol The absolute threshold for determining convergence of the
%          algorithm according to the duality gap.
%   RelTol The threshold for determining convergence of the algorithm
%          according to the relative duality gap.
%
%OUTPUTS: optTuples The updated set of tuples associated with fStar.
%             fStar The updated best cost function found thus far.
%          exitCode A value indicating the terminating state of the
%                   function. Possible values are:
%                  -1 The subproblem posted is not feasible.
%                   0 Convergence achieved.
%                   1 No errors occurred, but convegrence was not achieved.
%
%February 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

    n1=size(C,1);
    n3=size(C,3);

    CFeas=zeros(n1,n3);
    for i1=1:n1
        CFeas(i1,:)=reshape(C(i1,gamma1(i1),:),[1,n3]);
    end

    [gammaTilde3,~,f]=assign2D(CFeas,false);

    %If the subproblem is not feasible.
    if(isempty(gammaTilde3))
        exitCode=-1;
        return;
    end

    if(f<fStar)
        fStar=f;

        for i1=1:n1
           optTuples(:,i1)=[i1;gamma1(i1);gammaTilde3(i1)];
        end

        %The duality gap
        costGap=fStar-qStar;

        if(costGap<=AbsTol||(costGap<abs(qStar)*RelTol))
            exitCode=0;%The algorithm converged.
        else
            exitCode=1;
        end
    else
        exitCode=1;
    end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function filtS = processImageUsingPyradiomics(planC,strName,filtType,paramS)
%processImageUsingPyradiomics
%filtType : 'LoG', 'wavlelet'
% AI 06/12/2020

%% Get scan & mask
indexS = planC{end};

if ~isempty(strName)
    strC = {planC{indexS.structures}.structureName};
    strNum = getMatchingIndex(strName,strC,'exact');
    mask3M = getStrMask(strNum,planC);

    scanNum = getStructureAssociatedScan(strNum,planC);
else
    %Use entire scan
    scanNum = 1;
    mask3M = false(size(getScanArray(scanNum,planC)));
end
scan3M = double(getScanArray(scanNum,planC));
CToffset = planC{indexS.scan}(scanNum).scanInfo(1).CTOffset;
scan3M = scan3M - CToffset;

%% Get voxel size
scanS = planC{indexS.scan}(scanNum);
[xV,yV,zV] = getScanXYZVals(scanS);
dx = median(abs(diff(xV)));
dy = median(abs(diff(yV)));
dz = median(diff(zV));
voxelSizeV = [dx,dy,dz]*10; %convert to mm

%% Apply filters

%Add python module to system path & iImport
pyModule = 'pyProcessImage';

try
    py.importlib.import_module(pyModule);
catch
    disp('Python module could not be imported, check the pyradiomics path');
end

%Write scan & mask to NRRD format
fprintf('\nWriting scan and mask to NRRD format...\n');

originV = [0,0,0];
encoding = 'raw';

mask3M = uint16(mask3M);
mask3M = permute(mask3M, [2 1 3]);
mask3M = flip(mask3M,3);

scan3M = permute(scan3M, [2 1 3]);
scan3M = flip(scan3M,3);

scanFilename = strcat(tempdir,'scan.nrrd');
scanRes = nrrdWriter(scanFilename, scan3M, voxelSizeV, originV, encoding);

maskFilename = strcat(tempdir, 'mask.nrrd');
maskRes = nrrdWriter(maskFilename, mask3M, voxelSizeV, originV, encoding);

%Call image processing fn
try

    outPyList = py.pyProcessImage.filtImg(scanFilename, maskFilename,...
        filtType, paramS);
    filtImgC = outPyList{1};
    filtTypeC = outPyList{2};

    filtS = struct();
    paramC = fieldnames(paramS);

    for n = 1:length(filtImgC)
        %Convert python dictionary to matlab struct
        for m = 1:length(paramC)
            filtType = char(filtTypeC{n});
            filtType = strrep(filtType,'-','_');
            outFieldname = filtType;
        end
        pyFiltScan3M = double(filtImgC{n});
        pyFiltScan3M = permute(pyFiltScan3M,[2,3,1]);
        pyFiltScan3M = flip(pyFiltScan3M,3);
        filtS.(outFieldname) = pyFiltScan3M;
    end

catch e
    error('Feature extraction failed with message %s',e.message)
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function pathStr = dataDir(vw, ~)
% Directory containing data for a given dataTYPE
%
%  pathStr = dataDir(view)
%
% If you're nice and neat e.g., if you're in an inplane view and looking at
% average tSeries, the dataDir would be HOMEDIR/Inplane/Averages/
%
% If this directory doesn't exist, make it!
%
% djh, 2/2001
global HOMEDIR
if isempty(HOMEDIR), HOMEDIR = pwd; end
dir = viewGet(vw,'subdir');
dt = viewGet(vw,'dtStruct');
subDir = dtGet(dt,'Name');
pathStr = fullfile(HOMEDIR,dir,subDir);
if ~exist(pathStr,'dir')
	pathStr;
	subDir;
	mkdir(dir,subDir);
end
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ ret ] = setBlueOff( t_Kuka )
%% About:
% This function is used to set off the blue light of the KUKA iiwa 7 R 800.

%% Syntax:
% [ ret ] = setBlueOff( t_Kuka )

%% Arreguments:
% t_Kuka: is the TCP/IP connection

%% Return value:
% ret: a boolean variable
% true: if the message is treated successfully
% false: if an error has happened

% Copy right, Mohammad SAFEEA, 3rd of May 2017

theCommand='blueOff';
fprintf(t_Kuka, theCommand);
message=fgets(t_Kuka);
ret=checkAcknowledgment(message);
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function [params, transform] = sqexpKernExtractParam(kern)

% SQEXPKERNEXTRACTPARAM Extract parameters from the SQEXP kernel structure.
% FORMAT
% DESC Extract parameters from the pre-built compound squared
% exponential kernel structure into a vector of parameters for
% optimisation.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel matrix, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
%
% DESC Extract parameters and parameter names from the pre-built
% compound squared exponential kernel structure.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel matrix, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
% RETURN names : celly array of strings containing parameter names.
%
% SEEALSO sqexpKernParamInit, sqexpKernExpandParam, kernExtractParam, scg, conjgrad
%
% COPYRIGHT : Neil D. Lawrence, 2004
%
% KERN

params = [kern.inverseWidth kern.rbfVariance kern.biasVariance ...
          kern.whiteVariance];
if nargout > 1
  names{1} = 'sqexp inverse width';
  names{2} = 'sqexp rbf variance';
  names{3} = 'sqexp bias variance';
  names{4} = 'sqexp white variance';
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function visNNLearning(net,fig,transWeights);
%  visNNLearning(net);

layers = net.layers;

if notDefined('fig')
	fig = 99;
end

figure(fig);
clf

nLayers = numel(layers);

if notDefined('transWeights')
	transWeights = zeros(nLayers,1);
end

for iL = 1:nLayers
	subplot(3,nLayers,iL);
	visWeights(layers{iL}.W,transWeights(iL));
	title(sprintf('Layer %i Weights',iL));
	axis square
end

for iL = nLayers+1:2*nLayers
	subplot(3,nLayers,iL);
	visWeights(layers{iL-nLayers}.dW,transWeights(iL-nLayers));
	title(sprintf('Layer %i Gradients',iL-nLayers));
	axis square
end

for iL = 2*nLayers+1:3*nLayers
	subplot(3,nLayers,iL);
	hist(layers{iL-2*nLayers}.output(:),20);
	title(sprintf('Layer %i -- %s output',iL-2*nLayers,layers{iL-2*nLayers}.actFun));
	axis square
end
drawnow

author: dustinstansbury

repo: medal


text:
function x=usasi(n,fs)
%USASI generates N samples of USASI noise at sample frequency FS X=(N,FS)

% This routine is based on the USASI noise defined in [1] which was later
% reissued as [2]. USASI noise is intended to simulate the long-term average
% of typical audio program material. The routine does not currently implement
% the pulsation at 2.5Hz 12.5% duty cycle that is recommended by the standard.
% Also it should probably be scaled to a well-defined power.
%
%  [1] NRSC AM Reemphasis, Deemphasize, and Broadcast Audio Transmission Bandwidth Specifications,
%      EIA-549 Standard, Electronics Industries Association , July 1988.
%  [2] NRSC AM Reemphasis, Deemphasize, and Broadcast Audio Transmission Bandwidth Specifications,
%      NRSC-1-A Standard, Sept 2007, Online: http://www.nrscstandards.org/SG/NRSC-1-A.pdf

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: usasi.m,v 1.6 2008/05/27 08:38:05 dmb Exp $
%
%   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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<2 fs=8000; end
b=[1 0 -1];
a=poly(exp(-[100 320]*2*pi/fs));

x=randfilt(b,a,n);

author: decouples


text:
function [ZI,f] = spm_eeg_plotScalpData(Z,pos,ChanLabel,in)
% Display M/EEG interpolated sensor data on a scalp image
% FORMAT [ZI,f] = spm_eeg_plotScalpData(Z,pos,ChanLabel,in)
%
% INPUT:
%   Z          - the data matrix at the sensors
%   pos        - the positions of the sensors
%   ChanLabel  - the names of the sensors
%   in         - a structure containing some informations related to the
%                main PRESELECTDATA window. This entry is not necessary
% OUTPUT:
%   ZI         - an image of interpolated data onto the scalp
%   f          - the handle of the figure which displays the interpolated
%                data
%__________________________________________________________________________
%
% This function creates a figure whose purpose is to display an
% interpolation of the sensor data on the scalp (as an image).
%__________________________________________________________________________
% Copyright (C) 2008-2014 Wellcome Trust Centre for Neuroimaging

% Jean Daunizeau
% $Id: spm_eeg_plotScalpData.m 7221 2017-11-16 14:25:37Z vladimir $

ChanLabel  = char(ChanLabel);
ParentAxes = [];
f          = [];
clim       = [min(Z(:))-( max(Z(:))-min(Z(:)) )/63 , max(Z(:))];
figName    = 'Image Scalp data';
noButtons  = 0;
if nargin < 4 || isempty(in)
    in     = [];
else
    if isfield(in,'min') && ...
            isfield(in,'max') && ...
            isfield(in,'type')
        clim    = [in.min, in.max];
        dc      = abs(diff(clim))./63;
        clim(1) = clim(1) - dc;
        figName = ['Image Scalp data: ',in.type,' sensors'];
        if isfield(in,'trN')
            figName = [figName ', trial #',num2str(in.trN),'.'];
        end
    end
    if isfield(in,'f')
        f = in.f;
    else
        f = figure;
    end
    if isfield(in,'ParentAxes')
        ParentAxes = in.ParentAxes;
    else
        ParentAxes = axes('parent',f);
    end
    if isfield(in,'noButtons')
        noButtons = ~~in.noButtons;
    end
end

if ~isfield(in,'cbar')
    in.cbar = 1;
end

if ~isfield(in,'plotpos')
    in.plotpos = 1;
end

if size(pos,2) ~= size(ChanLabel, 1)
    pos = pos';
end

nD = size(pos,1);
if nD ~= 2
    % get 2D positions from 3D positions
   xyz   = pos;
   [pos] = get2Dfrom3D(xyz);
   pos   = pos';
end

% exclude channels ?
goodChannels = find(~isnan(pos(1,:)));
pos          = pos(:,goodChannels);
Z            = Z(goodChannels,:);
ChanLabel    = ChanLabel(goodChannels, :);

if ~isempty(in) && isfield(in,'type') && strcmp(in.type, 'MEGPLANAR')
    [cZ, cpos, cChanLabel] = combineplanar(Z, pos, ChanLabel);
else
    cZ         = Z;
    cpos       = pos;
    cChanLabel = ChanLabel;
end

xmin    = min(cpos(1,:));
xmax    = max(cpos(1,:));
dx      = (xmax-xmin)./100;
ymin    = min(cpos(2,:));
ymax    = max(cpos(2,:));
dy      = (ymax-ymin)./100;
x       = xmin:dx:xmax;
y       = ymin:dy:ymax;
[XI,YI] = meshgrid(x,y);
ZI      = griddata(cpos(1,:)',cpos(2,:)',full(double(cZ')),XI,YI);

try
    figure(f)
catch
    f   = figure(...
        'name',      figName,...
        'color',     [1 1 1],...
        'deleteFcn', @dFcn);
    ParentAxes = axes('parent',f);
end

COLOR = get(f,'color');
d.hi = image(flipud(ZI),...
    'CDataMapping','scaled',...
    'Parent',ParentAxes);
set(ParentAxes,'nextPlot','add',...
    'tag','spm_eeg_plotScalpData')
try
    if length(unique(ZI)) ~= 1
        [C,d.hc] = contour(ParentAxes,flipud(ZI),...
            'linecolor',0.5.*ones(3,1));
    end
end
caxis(ParentAxes,clim);
col = jet;
col(1,:) = COLOR;
colormap(ParentAxes,col)

if in.cbar
    d.cbar = colorbar('peer',ParentAxes);
end

axis(ParentAxes,'off')
axis(ParentAxes,'equal')
axis(ParentAxes,'tight')

fpos = cpos;
fpos(1,:) = fpos(1,:) - xmin;
fpos(2,:) = fpos(2,:) - ymin;
fpos(1,:) = fpos(1,:)./(dx);
fpos(2,:) = fpos(2,:)./(dy);
fpos(2,:) = 100-fpos(2,:);  % for display purposes (flipud imagesc)

figure(f);
if in.plotpos
    d.hp = plot(ParentAxes,...
        fpos(1,:),fpos(2,:),...
        'ko');
end

d.ht = text(fpos(1,:),fpos(2,:),cChanLabel,...
    'Parent',ParentAxes,...
    'visible','off');
axis(ParentAxes,'image')

d.interp.XI = XI;
d.interp.YI = YI;
d.interp.pos = cpos;
d.f = f;
d.pos = fpos;
d.goodChannels = goodChannels;
d.ChanLabel = cChanLabel;
d.origChanLabel = ChanLabel;
d.origpos = pos;
d.ParentAxes = ParentAxes;
d.in = in;

if ~noButtons
    d.hsp = uicontrol(f,...
        'style','pushbutton',...
        'callback',{@dosp},...
        'BusyAction','cancel',...
        'Interruptible','off',...
        'position',[10    50    80    20],...
        'string','channel pos');
    d.hsn = uicontrol(f,...
        'style','pushbutton',...
        'callback',{@dosn},...
        'BusyAction','cancel',...
        'Interruptible','off',...
        'position',[10    80    80    20],...
        'string','channel names');
end
if ~isempty(in) && isfield(in,'handles')
    nT = length(in.gridTime);
    d.hti = uicontrol(f,...
        'style','text',...
        'BackgroundColor',COLOR,...
        'string',[num2str(in.gridTime(in.x)),' (',in.unit,')'],...
        'position',[10    10    120    20]);
    d.hts = uicontrol(f,...
        'style','slider',...
        'Position',[130 10 250 20],...
        'min',1,'max',nT,...
        'value',in.x,'sliderstep',[1./(nT-1) 1./(nT-1)],...
        'callback',{@doChangeTime},...
        'BusyAction','cancel',...
        'Interruptible','off');
    set(d.hti,'userdata',d);
    set(d.hts,'userdata',d);
end
if ~noButtons
    set(d.hsp,'userdata',d);
    set(d.hsn,'userdata',d);
end
set(d.ParentAxes,'userdata',d);

%==========================================================================
% dFcn
%==========================================================================
function dFcn(btn,evd)
hf = findobj('tag','Graphics');
D = get(hf,'userdata');
try delete(D.PSD.handles.hli); end

%==========================================================================
% dosp
%==========================================================================
function dosp(btn,evd)
d = get(btn,'userdata');
switch get(d.hp,'visible');
    case 'on'
        set(d.hp,'visible','off');
    case 'off'
        set(d.hp,'visible','on');
end

%==========================================================================
% dosn
%==========================================================================
function dosn(btn,evd)
d = get(btn,'userdata');
switch get(d.ht(1),'visible')
    case 'on'
        set(d.ht,'visible','off');
    case 'off'
        set(d.ht,'visible','on');
end

%==========================================================================
% doChangeTime
%==========================================================================
function doChangeTime(btn,evd)
d = get(btn,'userdata');
v = get(btn,'value');
% get data
if ishandle(d.in.handles.hfig)
    D = get(d.in.handles.hfig,'userdata');
    if ~isfield(d.in,'trN')
        trN = 1;
    else
        trN = d.in.trN;
    end
    try
        Z = D(d.in.ind,v,trN);
        Z = Z(d.goodChannels);

        if strcmp(d.in.type, 'MEGPLANAR')
            Z = combineplanar(Z, d.origpos, d.origChanLabel);
        end

        clear ud;
        % interpolate data
        ZI = griddata(d.interp.pos(1,:),d.interp.pos(2,:),full(double(Z)),d.interp.XI,d.interp.YI);
        % update data display
        set(d.hi,'Cdata',flipud(ZI));
        % update time index display
        v = round(v);
        set(d.hti,'string',[num2str(d.in.gridTime(v)), ' (', d.in.unit, ')']);
        % update display marker position
        try;set(d.in.hl,'xdata',[v;v]);end
        set(d.ParentAxes,'nextPlot','add')
        try
            % delete current contour plot
            delete(findobj(d.ParentAxes,'type','hggroup'));
            delete(findobj(d.ParentAxes,'type','contour')); % R2014b
            % create new one
            [C,hc] = contour(d.ParentAxes,flipud(ZI),...
                'linecolor',[0.5.*ones(3,1)]);
        end
        axis(d.ParentAxes,'image')
        drawnow
    catch
%     else
        error('Did not find the data!')
    end
else
    error('SPM Graphics Figure has been deleted!')
end

%==========================================================================
% get2Dfrom3D
%==========================================================================
function [xy] = get2Dfrom3D(xyz)
% function [xy] = get2Dfrom3D(xyz)
% This function is used to flatten 3D sensor positions onto the 2D plane
% using a modified spherical projection operation.
% It is used to visualize channel data.
% IN:
%   - xyz: the cartesian sensor position in 3D space
% OUT:
%   - xy: the (x,y) cartesian coordinates of the sensors after projection
%   onto the best-fitting sphere

if size(xyz,2) ~= 3
    xyz = xyz';
end
% exclude channels ?
badChannels  = find(isnan(xyz(:,1)));
goodChannels = find(isnan(xyz(:,1))~=1);
xyz          = xyz(goodChannels,:);
% Fit sphere to 3d sensors and center frame
C            = fitSphere(xyz(:,1),xyz(:,2),xyz(:,3));
xyz          = xyz - repmat(C,size(xyz,1),1);
% apply transformation using spherical coordinates
[TH,PHI,RAD] = cart2sph(xyz(:,1),xyz(:,2),xyz(:,3));
TH           = TH - mean(TH);
[X,Y,Z]      = sph2cart(TH,zeros(size(TH)),RAD.*(cos(PHI+pi./2)+1));
xy           = [X(:),Y(:)];

%==========================================================================
% combineplanar
%==========================================================================
function [Z, pos, ChanLabel] = combineplanar(Z, pos, ChanLabel)

if ~iscell(ChanLabel)
    ChanLabel = cellstr(ChanLabel);
end

chanind = zeros(1, numel(ChanLabel));
for i = 1:numel(ChanLabel)
    chanind(i) = sscanf(ChanLabel{i}, 'MEG%d');
end

pairs = [];
unpaired = [];
paired = zeros(length(chanind));
for i = 1:length(chanind)
    if ~paired(i)

        cpair = find(abs(chanind - chanind(i))<2);

        if length(cpair) == 1
            unpaired = [unpaired cpair];
        else
            pairs = [pairs; cpair(:)'];
        end
        paired(cpair) = 1;
    end
end

if ~isempty(unpaired)
    warning(['Could not pair all channels. Ignoring ' num2str(length(unpaired)) ' unpaired channels.']);
end

Z = sqrt(Z(pairs(:, 1)).^2 + Z(pairs(:, 2)).^2);
pos = (pos(:, pairs(:, 1)) + pos(:, pairs(:, 2)))./2;
ChanLabel = {};
for i = 1:size(pairs,1)
    ChanLabel{i} = ['MEG' num2str(min(pairs(i,:))) '+' num2str(max(pairs(i,:)))];
end

%==========================================================================
% fitSphere
%==========================================================================
function [C,R,out] = fitSphere(x,y,z)
% fitSphere  Fit sphere.
%       A = fitSphere(x,y,z) returns the parameters of the best-fit
%       [C,R,out] = fitSphere(x,y,z) returns the center and radius
%       sphere to data points in vectors (x,y,z) using Taubin's method.
% IN:
%   - x/y/z: 3D carthesian ccordinates
% OUT:
%   - C: the center of sphere coordinates
%   - R: the radius of the sphere
%   - out: an output structure devoted to graphical display of the best fit
%   sphere

% Make sugary one and zero vectors
l = ones(length(x),1);
O = zeros(length(x),1);

% Make design mx
D = [(x.*x + y.*y + z.*z) x y z l];

Dx = [2*x l O O O];
Dy = [2*y O l O O];
Dz = [2*z O O l O];

% Create scatter matrices
M = D'*D;
N = Dx'*Dx + Dy'*Dy + Dz'*Dz;

% Extract eigensystem
[v, evalues] = eig(M);
evalues = diag(evalues);
Mrank = sum(evalues > eps*5*norm(M));

if (Mrank == 5)
    % Full rank -- min ev corresponds to solution
    Minverse = v'*diag(1./evalues)*v;
    [v,evalues] = eig(inv(M)*N);
    [dmin,dminindex] = max(diag(evalues));
    pvec = v(:,dminindex(1))';
else
    % Rank deficient -- just extract nullspace of M
    pvec = null(M)';
    [m,n] = size(pvec);
    if m > 1
        pvec = pvec(1,:);
    end
end

% Convert to (R,C)
if nargout == 1
    if pvec(1) < 0
        pvec = -pvec;
    end
    C = pvec;
else
    C = -0.5*pvec(2:4) / pvec(1);
    R = sqrt(sum(C*C') - pvec(5)/pvec(1));
end

[X,Y,Z] = sphere;
[TH,PHI,R0] = cart2sph(X,Y,Z);
[X,Y,Z] = sph2cart(TH,PHI,R);
X = X + C(1);
Y = Y + C(2);
Z = Z + C(3);

out.X = X;
out.Y = Y;
out.Z = Z;

author: spm


text: mex libsvm_classifier_spider.cpp svm.cpp -I..\compiletools ..\compiletools\mexarg.cpp

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function D = spm_eeg_contrast(S)
% Compute contrasts over trials or trial types
% FORMAT D = spm_eeg_contrast(S)
%
% S         - optional input struct
% fields of S:
% D         - filename of EEG mat-file with epoched data
% c         - contrast matrix, each row computes a contrast of the data
% label     - cell array of labels for the contrasts, the same size as
%             number of rows in c
% weighted  - flag whether average should be weighted by number of
%             replications (yes (1), no (0))
% prefix    - prefix for the output file [default: 'w']
%
% Output:
% D         - EEG data struct (also written to disk)
%__________________________________________________________________________
%
% spm_eeg_contrast computes contrasts of data, over epochs of data. The
% input is a single MEEG file. The argument c must have dimensions
% Ncontrasts X Nepochs, where Ncontrasts is the number of contrasts and
% Nepochs the number of epochs, i.e. each row of c contains one contrast
% vector. The output is a M/EEG file with Ncontrasts epochs. The typical
% use is to compute, for display purposes, contrasts like the difference or
% interaction between trial types in channel space.
%__________________________________________________________________________
% Copyright (C) 2008-2017 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel, Rik Henson
% $Id: spm_eeg_contrast.m 7132 2017-07-10 16:22:58Z guillaume $

SVNrev = '$Rev: 7132 $';

%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename, SVNrev);
spm('FigName','M/EEG Contrast'); spm('Pointer','Watch');

%-Get MEEG object
%--------------------------------------------------------------------------
if ~isfield(S, 'prefix'),       S.prefix     = 'w';          end
if ~isfield(S, 'weighted'),     S.weighted   =  0;           end

if ~(isfield(S, 'c') && isfield(S, 'label') && numel(S.label)==size(S.c, 1))
    error('Invalid contrast specification.');
end

D = spm_eeg_load(S.D);

%-Compute contrasts
%--------------------------------------------------------------------------

c          = S.c;
Ncontrasts = size(c, 1);

% Pad with zeros as in the contrast manager
if size(c, 2) ~= D.ntrials
    error('The number of columns in the contrast matrix does not match the number of trials.');
end

if ~isempty(D.repl)
    weighted = S.weighted;
else
    weighted = 0;
end

if strncmp(D.transformtype, 'TF', 2)
    Dnew = clone(D, [S.prefix fname(D)], [D.nchannels D.nfrequencies D.nsamples Ncontrasts]);
else
    Dnew = clone(D,  [S.prefix fname(D)],[D.nchannels D.nsamples Ncontrasts]);
end

spm_progress_bar('Init', Ncontrasts, 'Contrasts computed');
if Ncontrasts > 100, Ibar = floor(linspace(1, Ncontrasts, 100));
else Ibar = 1:Ncontrasts; end

for i = 1:Ncontrasts

    if weighted
        p = find(c(i,:) == 1);
        if ~isempty(p)
            r = D.repl(p);
            c(i,p) = r/sum(r);
        end

        p = find(c(i,:) == -1);
        if ~isempty(p)
            r = D.repl(p);
            c(i,p) = -r/sum(r);
        end
    end

    disp(['Contrast ', mat2str(i),': ', mat2str(c(i,:),3)])

    if strncmp(D.transformtype, 'TF', 2)
        d = zeros(D.nchannels, D.nfrequencies, D.nsamples);

        for j = 1:D.nchannels
            for f = 1:D.nfrequencies
                d(j, f, :) = c(i,:) * spm_squeeze(D(j, f, :, :), [1 2 4])';
            end
        end

        Dnew(1:Dnew.nchannels, 1:D.nfrequencies, 1:Dnew.nsamples, i) = d;

    else

        d = zeros(D.nchannels, D.nsamples);

        for j = 1:D.nchannels
            d(j, :) = c(i,:) * squeeze(D(j, :, :))';
        end

        Dnew(1:Dnew.nchannels, 1:Dnew.nsamples, i) = d;
    end

    newrepl(i) = sum(D.repl(find(c(i,:)~=0)));

    if any(Ibar == i), spm_progress_bar('Set', i); end

end

spm_progress_bar('Clear');

%-
%--------------------------------------------------------------------------
Dnew = conditions(Dnew, ':', S.label);
Dnew = trialonset(Dnew, ':', []);
Dnew = trialtag(Dnew, ':', []);
Dnew = badtrials(Dnew, ':', 0);
Dnew = repl(Dnew, ':', newrepl);
Dnew = type(Dnew, 'evoked');

% remove previous source reconsructions
if isfield(Dnew,'inv')
    Dnew = rmfield(Dnew,'inv');
end

%-Save new M/EEG data
%--------------------------------------------------------------------------
D = Dnew;
D = D.history('spm_eeg_contrast', S);
save(D);

%-Cleanup
%--------------------------------------------------------------------------
spm('FigName','M/EEG Contrast: done'); spm('Pointer','Arrow');


text:
function [psi] = Pa2psi(Pa)
% Convert pressure from pascals to pounds-per-square-inch.
% Chad Greene 2012
psi = Pa/6894.757;

author: Sable


text:
%PNMDEMO Demonstration of portable bitmap utilities.

%   Author:      Peter J. Acklam
%   Time-stamp:  1998-05-11 16:39:20
%   E-mail:      jacklam@math.uio.no (Internet)
%   URL:         http://www.math.uio.no/~jacklam

more off;
echo on;
home; clc;
%
% First, read the portable pixelmap with the Matlab logo.
%
% The image may be read as an RGB image.
%

[ r, g, b ] = ppmread( 'logo5' );

pause

%
% The image may also be read as an indexed image.
%

[ x, map ] = ppmread( 'logo5' );

pause

%
% Since PPMREAD does not simplify the color map automatically, the number
% of colors in the map equals the number of pixels in the image.
%

pixels_in_image = prod( size(x) )
colors_in_map  = size( map, 1 )

pause

%
% Use IM2IND to convert intensity images, indexed images and RGB images
% to indexed images with the smallest possible color map.
%

[ x, map ] = im2ind( x, map );

pause

%
% Now we see the number of unique colors in the image. The number of
% colors in the map has been reduced significantly.
%

colors_in_map = size( map, 1 )

pause

%
% Display the image.
%

figure( gcf ); clf;
colormap( map ); image( x ); axis image;

pause

%
% Write the image as a portable graymap.
% PGMWRITE automatically converts color images to grayscale images.
% For speed, we use binary (raw) encoding.
%

pgmwrite( x, map, 'logo5', 'binary' );

pause

%
% Read the image. We can use both PGMREAD and PPMREAD for this.
% PGMREAD can read the image as an intensity image or an indexed image.
% We choose the latter.
%

[ x, map ] = pgmread( 'logo5' );

pause

%
% Display the image.
%

figure( gcf ); clf;
colormap( map ); image( x ); axis image;

pause

%
% As with PPMREAD, PGMREAD does not simplify the color map
% automatically, so the number of colors in the map may be larger than
% the number of colors used in the image.
%

colors_in_map = size( map, 1 )

pause

%
% Again, we use IM2IND to get the indexed image with the smallest
% possible color map.
%

[ x, map ] = im2ind( x, map );

pause

%
% Now we see the number of unique colors in the image.
%

colors_in_map = size( map, 1 )

pause

%
% Now write the image as a portable bitmap.
% PBMWRITE automatically converts color images and grayscale images to
% black and white images. For speed, we use binary (raw) encoding.
%

pbmwrite( x, map, 'logo5', 'binary' );

pause

%
% Read the image. We can use both PBMREAD, PGMREAD and PPMREAD to this.
% PBMREAD can read the image as a black and white image or as an indexed
% image. We choose the latter.
%

[ x, map ] = pbmread( 'logo5' );

pause

%
% Display the image.
%

figure( gcf ); clf;
colormap( map ); image( x ); axis image;

pause

%
% That should be all for now, so clean up.
%

delete logo5.pgm logo5.pbm


text:
function calpak_test51 ( )

%*****************************************************************************80
%
%% CALPAK_TEST51 tests WEEKDAY_TO_NAME_REPUBLICAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST51\n' );
  fprintf ( 1, '  For the Republican calendar,\n' );
  fprintf ( 1, '  WEEKDAY_TO_NAME_REPUBLICAN names the days of the week.\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10
    s = weekday_to_name_republican ( i );
    fprintf ( 1, '  %2d  %s\n', i, s );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function r = intrand(p, m, n);
%INTRAND Random matrices from uniform integer distribution.
%
%   R = INTRAND(P) returns a matrix of random numbers chosen
%   from the uniform integer distribution with parameter P. The
%   minimum and maximum output integer numbers are P(1) and P(2).
%   R = INTRAND(P,M,N) returns an M by N matrix.
%
%
% Copyright (c) 2000 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.

if nargin < 1 | any(size([p(:)])~=[2 1])
  error('Incorrect range assignment in INTRAND');
end
if p(1)<p(2)
  rmin=ceil(p(1));
  rrange=floor(p(2))-rmin+1;
else
  rmin=ceil(p(2));
  rrange=floor(p(1))-rmin+1;
end

if nargin < 2
  r=floor(rand.*rrange+rmin);
else
  if nargin < 3
    n=m;
  end
  r=floor(rand([m n]).*rrange+rmin);
end

author: gpstuff-dev

repo: gpstuff

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


text:
function a = sweet3 (  )

%*****************************************************************************80
%
%% SWEET3 returns the SWEET3 matrix.
%
%  Example:
%
%      8    4    1    6    2    3
%      4    8    4    1    6    2
%    -34    4    8    4    1    6
%      5  -34    4    8    4    1
%      3    5  -34    4    8    4
%      1    3    5  -34    4    8
%
%  Properties:
%
%    A is Toeplitz: constant along diagonals.
%
%    A is generally not symmetric: A' /= A.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per Hansen, Tony Chan,
%    FORTRAN Subroutines for General Toeplitz Systems,
%    ACM Transactions on Mathematical Software,
%    Volume 18, Number 3, September 1992, pages 256-273.
%
%    Douglas Sweet,
%    The use of pivoting to improve the numerical performance of
%    Toeplitz solvers,
%    In "Advanced Algorithms and Architectures for Signal Processing",
%    Edited by J M Speiser,
%    Proceedings SPIE 696, 1986, pages 8-18.
%
%  Parameters:
%
%    Output, real A(6,6), the matrix.
%
  n = 6;

  v = [ 1.0, 3.0, 5.0, -34.0, 4.0, 8.0, 4.0, ...
    1.0, 6.0, 2.0, 3.0 ];

  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n

      a(i,j) = v ( j - i + 6 );

    end
  end

  return
end