text:
function [sourcemodel, cfg] = ft_prepare_leadfield(cfg, data)

% FT_PREPARE_LEADFIELD computes the forward model for many dipole locations
% on a regular 2D or 3D sourcemodel and stores it for efficient inverse modelling
%
% Use as
%   [sourcemodel] = ft_prepare_leadfield(cfg, data)
%
% It is necessary to input the data on which you want to perform the inverse
% computations, since that data generally contain the gradiometer information and
% information about the channels that should be included in the forward model
% computation. The data structure can be either obtained from FT_PREPROCESSING,
% FT_FREQANALYSIS or FT_TIMELOCKANALYSIS. If the data is empty, all channels will be
% included in the forward model.
%
% The configuration should contain
%   cfg.channel            = Nx1 cell-array with selection of channels (default = 'all'),
%                            see FT_CHANNELSELECTION for details
%
% The positions of the sources can be specified as a regular 3-D
% sourcemodel that is aligned with the axes of the head coordinate system
%   cfg.xgrid      = vector (e.g. -20:1:20) or 'auto' (default = 'auto')
%   cfg.ygrid      = vector (e.g. -20:1:20) or 'auto' (default = 'auto')
%   cfg.zgrid      = vector (e.g.   0:1:20) or 'auto' (default = 'auto')
%   cfg.resolution = number (e.g. 1 cm) for automatic sourcemodel generation
%
% Alternatively the position of a few sources at locations of interest can
% be specified, for example obtained from an anatomical or functional MRI
%   cfg.sourcemodel.pos        = N*3 matrix with position of each source
%   cfg.sourcemodel.inside     = N*1 vector with boolean value whether sourcemodel point is inside brain (optional)
%   cfg.sourcemodel.dim        = [Nx Ny Nz] vector with dimensions in case of 3-D sourcemodel (optional)
%
% The volume conduction model of the head should be specified as
%   cfg.headmodel     = structure with volume conduction model, see FT_PREPARE_HEADMODEL
%
% The EEG or MEG sensor positions can be present in the data or can be specified as
%   cfg.elec          = structure with electrode positions or filename, see FT_READ_SENS
%   cfg.grad          = structure with gradiometer definition or filename, see FT_READ_SENS
%
% Optionally, you can modify the leadfields by reducing the rank (i.e. remove the
% weakest orientation), or by normalizing each column.
%   cfg.reducerank      = 'no', or number (default = 3 for EEG, 2 for MEG)
%   cfg.backproject     = 'yes' or 'no',  determines when reducerank is applied whether the
%                         lower rank leadfield is projected back onto the original linear
%                         subspace, or not (default = 'yes')
%   cfg.normalize       = 'yes' or 'no' (default = 'no')
%   cfg.normalizeparam  = depth normalization parameter (default = 0.5)
%   cfg.weight          = number or Nx1 vector, weight for each dipole position to compensate
%                         for the size of the corresponding patch (default = 1)
%
% Depending on the type of headmodel, some additional options may be
% specified.
%
% For OPENMEEG based headmodels:
%   cfg.openmeeg.batchsize    = scalar (default 1e4), number of dipoles
%                               for which the leadfield is computed in a
%                               single call to the low-level code. Trades off
%                               memory efficiency for speed.
%   cfg.openmeeg.dsm          = 'no'/'yes', reuse existing DSM if provided
%   cfg.openmeeg.keepdsm      = 'no'/'yes', option to retain DSM (no by default)
%   cfg.openmeeg.nonadaptive  = 'no'/'yes'
%
% For SINGLESHELL based headmodels:
%   cfg.singleshell.batchsize = scalar or 'all' (default 1), number of dipoles
%                               for which the leadfield is computed in a
%                               single call to the low-level code. Trades off
%                               memory efficiency for speed.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% If you specify this option the input data will be read from a *.mat
% file on disk. This mat files should contain only a single variable named 'data',
% corresponding to the input structure.
%
% See also FT_SOURCEANALYSIS, FT_DIPOLEFITTING, FT_PREPARE_HEADMODEL, FT_PREPARE_SOURCEMODEL

% Undocumented local options:
% cfg.feedback
% cfg.sel50p      = 'no' (default) or 'yes'
% cfg.lbex        = 'no' (default) or a number that corresponds with the radius
% cfg.mollify     = 'no' (default) or a number that corresponds with the FWHM

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar data
ft_preamble provenance data

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% the data can be passed as input arguments or can be read from disk
hasdata = exist('data', 'var');

if ~hasdata
  % the data variable will be passed to the prepare_headmodel function below
  % where it would be used for channel selection
  data = [];
else
  % check if the input data is valid for this function
  data = ft_checkdata(data);
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels'}); % prevent accidental typos, see issue 1729
cfg = ft_checkconfig(cfg, 'renamed',    {'hdmfile',   'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'vol',       'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'grid',      'sourcemodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'om',        'openmeeg'});
cfg = ft_checkconfig(cfg, 'renamed',    {'elecfile',  'elec'});
cfg = ft_checkconfig(cfg, 'renamed',    {'gradfile',  'grad'});
cfg = ft_checkconfig(cfg, 'renamed',    {'optofile',  'opto'});
cfg = ft_checkconfig(cfg, 'deprecated', {'patchindx', 'patchsize'});

% set the defaults
cfg.lbex           = ft_getopt(cfg, 'lbex',      'no');
cfg.sel50p         = ft_getopt(cfg, 'sel50p',    'no');
cfg.feedback       = ft_getopt(cfg, 'feedback',  'text');
cfg.mollify        = ft_getopt(cfg, 'mollify',   'no');
cfg.patchsvd       = ft_getopt(cfg, 'patchsvd',  'no');

cfg = ft_checkconfig(cfg, 'renamed', {'tightgrid',   'tight'});  % this is moved to cfg.sourcemodel.tight by the subsequent createsubcfg
cfg = ft_checkconfig(cfg, 'renamed', {'sourceunits', 'unit'});   % this is moved to cfg.sourcemodel.unit by the subsequent createsubcfg

% put the low-level options pertaining to the sourcemodel in their own field
cfg = ft_checkconfig(cfg, 'createsubcfg', {'sourcemodel'});
% move some fields from cfg.sourcemodel back to the top-level configuration
cfg = ft_checkconfig(cfg, 'createtopcfg', {'sourcemodel'});

% this code expects the inside to be represented as a logical array
cfg.sourcemodel = ft_checkconfig(cfg.sourcemodel, 'renamed',  {'pnt' 'pos'});
cfg = ft_checkconfig(cfg, 'inside2logical', 'yes');

if strcmp(cfg.sel50p, 'yes') && strcmp(cfg.lbex, 'yes')
  ft_error('subspace projection with either lbex or sel50p is mutually exclusive');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% collect and preprocess the electrodes/gradiometer and head model
[headmodel, sens, cfg] = prepare_headmodel(cfg, data);

% construct the sourcemodel for which the leadfield will be computed
tmpcfg           = keepfields(cfg, {'sourcemodel', 'mri', 'headshape', 'symmetry', 'smooth', 'threshold', 'spheremesh', 'inwardshift', 'xgrid' 'ygrid', 'zgrid', 'resolution', 'tight', 'warpmni', 'template', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
tmpcfg.headmodel = headmodel;
if ft_senstype(sens, 'eeg')
  tmpcfg.elec = sens;
elseif ft_senstype(sens, 'meg')
  tmpcfg.grad = sens;
end
sourcemodel = ft_prepare_sourcemodel(tmpcfg);

% find the indices of all sourcemodel points that are inside the brain
insideindx = find(sourcemodel.inside);

% check whether units are equal (NOTE: this was previously not required,
% this check can be removed if the underlying bug is resolved. See
% http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2387
if ~isfield(headmodel, 'unit') || ~isfield(sourcemodel, 'unit') || ~isfield(sens, 'unit')
  ft_warning('cannot determine the units of all geometric objects required for leadfield computation (headmodel, sourcemodel, sensor configuration). THIS CAN LEAD TO WRONG RESULTS! (refer to http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2387)');
else
  if ~strcmp(headmodel.unit, sourcemodel.unit) || ~strcmp(sourcemodel.unit, sens.unit)
    ft_error('geometric objects (headmodel, sourcemodel, sensor configuration) are not expressed in the same units (this used to be allowed, and will be again in the future, but for now there is a bug which prevents a correct leadfield from being computed; see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2387)');
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% construct the low-level options for the leadfield computation as key-value pairs, these are passed to FT_COMPUTE_LEADFIELD
leadfieldopt = {};
leadfieldopt = ft_setopt(leadfieldopt, 'reducerank',     ft_getopt(cfg, 'reducerank'));
leadfieldopt = ft_setopt(leadfieldopt, 'backproject',    ft_getopt(cfg, 'backproject'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalize',      ft_getopt(cfg, 'normalize'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalizeparam', ft_getopt(cfg, 'normalizeparam'));
leadfieldopt = ft_setopt(leadfieldopt, 'weight',         ft_getopt(cfg, 'weight'));

if ft_headmodeltype(headmodel, 'openmeeg')

  ft_hastoolbox('openmeeg', 1);  % add to path (if not yet on path)

  cfg.openmeeg = ft_getopt(cfg, 'openmeeg', []);
  batchsize    = ft_getopt(cfg.openmeeg, 'batchsize', 1e4);  % number of voxels per DSM batch; set to e.g. 1000 if not much RAM available
  keepdsm      = ft_getopt(cfg.openmeeg, 'keepdsm', 'no');   % retain DSM

  leadfieldopt = ft_setopt(leadfieldopt, 'dsm',         ft_getopt(cfg.openmeeg, 'dsm')); % reuse existing DSM if provided
  leadfieldopt = ft_setopt(leadfieldopt, 'nonadaptive', ft_getopt(cfg.openmeeg, 'nonadaptive', 'no'));

  % repeated system calls to the openmeeg executable makes it rather slow, calling it once is much more efficient
  fprintf('calculating leadfield for %d positions at a time, this may take a while...\n', batchsize);

  % a dsm in the input cfg currently assumes that the 'content' of the dsm matches
  % exactly the positions that were passed in sourcemodel.pos. This is not guaranteed
  % of course. If anything, it would make sense to represent the dsm in the input
  % sourcemodel, so that the pos, and dsm are bound together. Still no guarantee, but
  % better than nothing. This means that the cfg.openmeeg.dsm option should be
  % deprecated

  ndip       = length(insideindx);
  numchunks  = ceil(ndip/batchsize);
  dippos     = sourcemodel.pos(insideindx,:);

  if(numchunks > 1)
    if istrue(keepdsm)
      ft_warning('Keeping DSM output not supported when the computation is split into batches')
    end
    keepdsm = false;
  end

  % DSM computation is computationally intensive: As it can be reused with same voxel
  % sourcemodel (i.e. if voxels are defined in MRI coordinates rather than MEG
  % coordinates), optionally save result. Dense voxel grids may require several
  % gigabytes of RAM, so optionally split into smaller batches
  dsm = ft_getopt(leadfieldopt, 'dsm');
  if istrue(keepdsm) && ~isempty(dsm)
    % dsm needs to be computed outside ft_compute_leadfield, because it needs to be passed on in the output
    dsm          = ft_sysmat_openmeeg(dippos, headmodel, sens, ft_getopt(leadfieldopt, 'nonadaptive'));
    leadfieldopt = ft_setopt(leadfieldopt, 'dsm', dsm);
  end

  sourcemodel.leadfield = cell(size(sourcemodel.pos,1),1);
  ft_progress('init', cfg.feedback, 'computing leadfield');
  for k = 1:numchunks
    ft_progress(k/numchunks, 'computing leadfield %d/%d\n', k, numchunks);
    diprange = (((k-1)*batchsize + 1):(min(k*batchsize,ndip)));
    tmp      = ft_compute_leadfield(dippos(diprange,:), sens, headmodel, leadfieldopt{:});
    % distribute the columns of the leadfield matrix over the individual dipole positions
    % avoid using the options reducerank and backproject, see https://github.com/fieldtrip/fieldtrip/issues/1410#issuecomment-646994620
    [m, n] = size(tmp);
    sourcemodel.leadfield(insideindx(diprange)) = mat2cell(tmp, m, repmat(n/numel(diprange), 1, numel(diprange)));
  end
  ft_progress('close');

  if istrue(keepdsm)
    % retain DSM in cfg if desired -> FIXME this should not be kept in
    % the cfg. If anything, it is sourcemodel specific, so it should be
    % retained in the output sourcemodel
    cfg.openmeeg.dsm = dsm;
  end

elseif ft_headmodeltype(headmodel, 'singleshell')
  cfg.singleshell = ft_getopt(cfg, 'singleshell', []);
  batchsize       = ft_getopt(cfg.singleshell, 'batchsize', 1);
  if ischar(batchsize) && strcmp(batchsize, 'all')
    batchsize = length(insideindx);
  end

  dippos     = sourcemodel.pos(insideindx,:);
  ndip       = length(insideindx);
  numchunks  = ceil(ndip/batchsize);

  sourcemodel.leadfield = cell(size(sourcemodel.pos,1),1);
  ft_progress('init', cfg.feedback, 'computing leadfield');
  for k = 1:numchunks
    ft_progress(k/numchunks, 'computing leadfield %d/%d\n', k, numchunks);
    diprange = (((k-1)*batchsize + 1):(min(k*batchsize,ndip)));
    tmp      = ft_compute_leadfield(dippos(diprange,:), sens, headmodel, leadfieldopt{:});
    % distribute the columns of the leadfield matrix over the individual dipole positions
    % avoid using the options reducerank and backproject, see https://github.com/fieldtrip/fieldtrip/issues/1410#issuecomment-646994620
    [m, n] = size(tmp);
    sourcemodel.leadfield(insideindx(diprange)) = mat2cell(tmp, m, repmat(n/numel(diprange), 1, numel(diprange)));
  end
  ft_progress('close');

elseif ft_headmodeltype(headmodel, 'duneuro')
%   ft_hastoolbox('duneuro', 1); %does not look necessary here? check
  % repeated system calls to the duneuro executable makes it rather slow
  % calling it once for all dipoles is much more efficient

  % find the indices of all grid points that are inside the brain
  insideindx = find(sourcemodel.inside);

  ft_progress('init', cfg.feedback, 'computing leadfield');
  % compute the leadfield on all grid positions inside the brain
  lf = ft_compute_leadfield(sourcemodel.pos(insideindx,:), sens, headmodel, 'reducerank', cfg.reducerank, 'normalize', cfg.normalize, 'normalizeparam', cfg.normalizeparam, 'backproject', cfg.backproject);
  lf = mat2cell(lf, size(lf,1), repmat(3,1,size(lf,2)/3));
  sourcemodel.leadfield(sourcemodel.inside) = lf;
  for i=1:length(insideindx)
    thisindx = insideindx(i);
    if isfield(cfg, 'grid') && isfield(cfg.grid, 'mom')
      % multiply with the normalized dipole moment to get the leadfield in the desired orientation
      sourcemodel.leadfield{thisindx} = sourcemodel.leadfield{thisindx} * sourcemodel.mom(:,thisindx);
    end
  end % for all grid locations inside the brain

elseif ft_headmodeltype(headmodel, 'interpolate')

  lf = ft_compute_leadfield(sourcemodel.pos(insideindx,:), sens, headmodel, leadfieldopt{:});
  sourcemodel.leadfield(insideindx) = mat2cell(lf, 3, 3.*ones(1,numel(insideindx)));

else
  ft_progress('init', cfg.feedback, 'computing leadfield');
  for i=1:length(insideindx)
    % compute the leadfield on all sourcemodel positions inside the brain
    ft_progress(i/length(insideindx), 'computing leadfield %d/%d\n', i, length(insideindx));
    thisindx = insideindx(i);
    sourcemodel.leadfield{thisindx} = ft_compute_leadfield(sourcemodel.pos(thisindx,:), sens, headmodel, leadfieldopt{:});
  end % for all sourcemodel locations inside the brain
  ft_progress('close');
end

if isfield(cfg, 'sourcemodel') && isfield(cfg.sourcemodel, 'mom')
  for i=1:length(insideindx)
    % multiply with the normalized dipole moment to get the leadfield in the desired orientation
    % FIXME mom and ori seem to be mixed up here, see https://github.com/fieldtrip/fieldtrip/issues/1399
    thisindx = insideindx(i);
    sourcemodel.leadfield{thisindx} = sourcemodel.leadfield{thisindx} * sourcemodel.mom(:,thisindx);
  end
end

% represent the leadfield for positions outside the brain as empty array
sourcemodel.leadfield(~sourcemodel.inside) = {[]};

% add the label of the channels
sourcemodel.label           = sens.label;
sourcemodel.leadfielddimord = '{pos}_chan_ori';

% mollify the leadfields
if ~isequal(cfg.mollify, 'no')
  sourcemodel = mollify(cfg, sourcemodel);
end

% combine leadfields in patches and do an SVD on them
if ~isequal(cfg.patchsvd, 'no')
  sourcemodel = patchsvd(cfg, sourcemodel);
end

% compute the 50 percent channel selection subspace projection
if ~isequal(cfg.sel50p, 'no')
  sourcemodel = sel50p(cfg, sourcemodel, sens);
end

% compute the local basis function expansion (LBEX) subspace projection
if ~isequal(cfg.lbex, 'no')
  sourcemodel = lbex(cfg, sourcemodel);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   data
ft_postamble provenance sourcemodel
ft_postamble history    sourcemodel
ft_postamble savevar    sourcemodel


text:
function plane = createPlane(varargin)
%CREATEPLANE Create a plane in parametrized form.
%
%   PLANE = createPlane(P1, P2, P3)
%   creates a plane containing the 3 points
%
%   PLANE = createPlane(PTS)
%   The 3 points are packed into a single 3x3 array.
%
%   PLANE = createPlane(P0, N);
%   Creates a plane from a point and from a normal to the plane. The
%   parameter N is given either as a 3D vector (1-by-3 row vector), or as
%   [THETA PHI], where THETA is the colatitute (angle with the vertical
%   axis) and PHI is angle with Ox axis, counted counter-clockwise (both
%   given in radians).
%
%   PLANE = createPlane(P0, Dip, DipDir);
%   Creates a plane from a point and from a dip and dip direction angles
%   of the plane. Parameters Dip and DipDir angles are given as numbers.
%   Dip : maximum inclination to the horizontal.
%   DipDir : direction of the horizontal trace of the line of dip,
%            measured clockwise from north.
%
%   The created plane data has the following format:
%   PLANE = [X0 Y0 Z0  DX1 DY1 DZ1  DX2 DY2 DZ2], with
%   - (X0, Y0, Z0) is a point belonging to the plane
%   - (DX1, DY1, DZ1) is a first direction vector
%   - (DX2, DY2, DZ2) is a second direction vector
%   The 2 direction vectors are normalized and orthogonal.
%
%   See also
%   planes3d, medianPlane

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2005-02-18
% Copyright 2005-2022 INRA - TPV URPOI - BIA IMASTE

if length(varargin) == 1
    var = varargin{1};

    if iscell(var)
        plane = zeros([length(var) 9]);
        for i=1:length(var)
            plane(i,:) = createPlane(var{i});
        end
    elseif size(var, 1) >= 3
        % 3 points in a single array
        p1 = var(1,:);
        p2 = var(2,:);
        p3 = var(3,:);

        % create direction vectors
        v1 = p2 - p1;
        v2 = p3 - p1;

        % create plane
        plane = normalizePlane([p1 v1 v2]);
        return;
    end

elseif length(varargin) == 2
    % plane origin
    p0 = varargin{1};

    % second parameter is either a 3D vector or a 3D angle (2 params)
    var = varargin{2};
    if size(var, 2) == 2
        % normal is given in spherical coordinates
        n = sph2cart2([var ones(size(var, 1))]);
    elseif size(var, 2)==3
        % normal is given by a 3D vector
        n = normalizeVector3d(var);
    else
        error ('wrong number of parameters in createPlane');
    end

    % ensure same dimension for parameters
    if size(p0, 1)==1
        p0 = repmat(p0, [size(n, 1) 1]);
    end
    if size(n, 1)==1
        n = repmat(n, [size(p0, 1) 1]);
    end

    % find a vector not colinear to the normal
    v0 = repmat([1 0 0], [size(p0, 1) 1]);
    inds = vectorNorm3d(cross(n, v0, 2))<1e-14;
    v0(inds, :) = repmat([0 1 0], [sum(inds) 1]);
%     if abs(cross(n, v0, 2))<1e-14
%         v0 = repmat([0 1 0], [size(p0, 1) 1]);
%     end

    % create direction vectors
    v1 = normalizeVector3d(cross(n, v0, 2));
    v2 = -normalizeVector3d(cross(v1, n, 2));

    % concatenate result in the array representing the plane
    plane = [p0 v1 v2];
    return;

elseif length(varargin)==3
    var1 = varargin{1};
    var2 = varargin{2};
    var3 = varargin{3};

    if size(var1, 2) == 3 && size(var2, 2) == 3 && size(var3, 2) == 3
        p1 = var1;
        p2 = var2;
        p3 = var3;

        % create direction vectors
        v1 = p2 - p1;
        v2 = p3 - p1;

        plane = normalizePlane([p1 v1 v2]);
        return;
    elseif size(var1, 2) == 3 && size(var2, 2) == 1 && size(var3, 2) == 1
        p0 = var1;
        n = [sin(var2)*sin(var3) sin(var2)*cos(var3) cos(var2)];

        % find a vector not colinear to the normal
        v0 = repmat([1 0 0], [size(p0, 1) 1]);
        inds = vectorNorm3d(cross(n, v0, 2))<1e-14;
        v0(inds, :) = repmat([0 1 0], [sum(inds) 1]);

        % create direction vectors
        v1 = normalizeVector3d(cross(n, v0, 2));
        v2 = -normalizeVector3d(cross(v1, n, 2));

        % concatenate result in the array representing the plane
        plane = [p0 v1 v2];
        return;
    else
        error('Wrong argument in "createPlane".');
    end
else
    error('Wrong number of arguments in "createPlane".');
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% floatwrite() - Write data matrix to float file.
%
% Usage: >> floatwrite(data,filename, 'format')
%
% Inputs:
%   data     - write matrix data to specified file as four-byte floating point numbers.
%   filename - name of the file
%   'format' - The option FORMAT argument specifies the storage format as
%              defined by fopen. Default format is 'native'.
%   'transp|normal' - save the data transposed (.dat files) or not.
%
% Author: Sigurd Enghoff, CNL / Salk Institute, La Jolla, 7/1998
%
% See also: floatread(), fopen()

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

% 07-08-99  FORMAT argument added -se
% 02-08-00  new version included in toolbox -sm
% 01-25-02 reformated help & license, added links -ad

function A = floatwrite(A, fname, fform, transp)

if ~exist('fform')
	fform = 'native';
end
if nargin < 4
    transp = 'normal';
end;

if strcmpi(transp,'normal')
    if strcmpi(class(A), 'mmo')
        A = changefile(A, fname);
        return;
    elseif strcmpi(class(A), 'memmapdata')
        % check file to overwrite
        % -----------------------
        [fpath1 fname1 ext1] = fileparts(fname);
        [fpath2 fname2 ext2] = fileparts(A.data.Filename);
        if isempty(fpath1), fpath1 = pwd; end;

        fname1 = fullfile(fpath1, [fname1 ext1]);
        fname2 = fullfile(fpath2, [fname2 ext2]);
        if ~isempty(findstr(fname1, fname2))
            disp('Warning: raw data already saved in memory mapped file (no need to resave it)');
            return;
        end;

        fid = fopen(fname,'wb',fform);
        if fid == -1, error('Cannot write output file, check permission and space'); end;
        if size(A,3) > 1
            for ind = 1:size(A,3)
                tmpdata = A(:,:,ind);
                fwrite(fid,tmpdata,'float');
            end;
        else
            blocks = [ 1:round(size(A,2)/10):size(A,2)];
            if blocks(end) ~= size(A,2), blocks = [blocks size(A,2)]; end;
            for ind = 1:length(blocks)-1
                tmpdata = A(:, blocks(ind):blocks(ind+1));
                fwrite(fid,tmpdata,'float');
            end;
        end;
    else
        fid = fopen(fname,'wb',fform);
        if fid == -1, error('Cannot write output file, check permission and space'); end;
        fwrite(fid,A,'float');
    end;
else
    % save transposed
    for ind = 1:size(A,1)
        fwrite(fid,A(ind,:),'float');
    end;
end;
fclose(fid);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% OverlayRegCheck: shell for mrAlign3 to compare inplanes for a session
%                  with interpolated inplanes, given the bestrotvol.
%                   This is intended as an update for regCheckRotAll.
%
%
%
% (regCheckRotAll)- Script that interpolates the inplanes corresponding
%                to the estimated rotation (rot) and translation (trans)
%                and displays slice by slice the original inplanes, the
%                interpolated inplanes and a mosaic of both.
%
%  Oscar Nestares - 5/99
%  Rory Sayres - 8/04, updated to use the overlayVol interface.

% Size of the original inplanes
[NyI NxI NzI] = size(INPLANE.anat);
inp = INPLANE.anat;

% interpolating the inplanes
hmsgbox = msgbox('Wait while interpolating the inplanes...'); drawnow
inpMf = regInplanes(reshape(volume,[sagSize, numSlices]),...
                    NxI, NyI, NzI, scaleFac, rot, trans);

% Correct intensity (selecting 'No' doesn't do anything useful)
fprintf('Correcting intensity...\n');
Limit = 4;
IntFunc = 'regEstFilIntGrad'; PbyPflag = 0;
inp = regCorrMeanInt(inp);
% intensity estimation
[Int Noise] = feval(IntFunc, inp, PbyPflag);
% intensity normalization
inp = regCorrIntGradWiener(inp, Int, Noise);
% robust mean and contrast normalization
inp = regCorrContrast(inp,Limit);
% intensity estimation
[IntM NoiseM] = feval(IntFunc, inpMf, PbyPflag);
% intensity normalization
inpMf = regCorrIntGradWiener(inpMf, IntM, NoiseM);
% robust mean and contrast normalization
[inpMf, pM] = regCorrContrast(inpMf,Limit);

close(hmsgbox);

% ensure inplanes and interp inplanes are same size
if size(inp) ~= size(inpMf)
    for i = 1:size(inp,3)
        tmp(:,:,i) = imresize(inp(:,:,i),size(inpMf(:,:,i)));
    end
    inp = tmp;
end

% call external overlay interface
FF = overlayVolumes(inp,inpMf);

% % checking the alignment
% FF = figure;
% SS = get(0,'ScreenSize');
% set(FF,'Position', [1 -40 SS(3)/3  SS(4)-80])
% for k=1:size(inpMf,3)
%    figure(FF)
%    subplot(3,1,1)
%    imagesc(inp(:,:,k));
%    axis('image'); colormap('gray'); axis('off');
%    subplot(3,1,2)
%    imagesc(regMosaic(regNormal(inpMf(:,:,k),0,1),...
%                      regNormal(inp(:,:,k),0,1)));
%    axis('image'); colormap('gray'); axis('off');
%    subplot(3,1,3)
%    imagesc(inpMf(:,:,k));
%    axis('image'); colormap('gray'); axis('off');
%    disp('Press a key to continue...')
%    pause
% end
%
% close(FF)

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


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

%*****************************************************************************80
%
%% LOG_UNIFORM_PDF evaluates the Log Uniform PDF.
%
%  Discussion:
%
%    PDF(A,B;X) = 1 / ( X * ( log ( B ) - log ( A ) ) ) for A <= X <= B
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%
%    Input, real A, B, the parameters of the PDF.
%    1.0 < A < B.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x < a )
    pdf = 0.0;
  elseif ( x <= b )
    pdf = 1.0 / ( x * ( log ( b ) - log ( a ) ) );
  else
    pdf = 0.0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ x, ierror ] = r8mat_solve2 ( n, a, b, x )

%*****************************************************************************80
%
%% R8MAT_SOLVE2 computes the solution of an N by N linear system.
%
%  Discussion:
%
%    The linear system may be represented as
%
%      A*X = B
%
%    If the linear system is singular, but consistent, then the routine will
%    still produce a solution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of equations.
%
%    Input, real A(N,N), the coefficient matrix to be inverted.
%
%    Input, real B(N), the right hand side of the system.
%
%    Output, real X(N), the solution of the linear system.
%
%    Output, integer IERROR.
%    0, no error detected.
%    1, consistent singularity.
%    2, inconsistent singularity.
%
  ierror = 0;

  ipiv(1:n) = 0;
  x(1:n) = 0.0;
%
%  Process the matrix.
%
  for k = 1 : n
%
%  In column K:
%    Seek the row IMAX with the properties that:
%      IMAX has not already been used as a pivot;
%      A(IMAX,K) is larger in magnitude than any other candidate.
%
    amax = 0.0;
    imax = 0;
    for i = 1 : n
      if ( ipiv(i) == 0 )
        if ( amax < abs ( a(i,k) ) )
          imax = i;
          amax = abs ( a(i,k) );
        end
      end
    end
%
%  If you found a pivot row IMAX, then,
%    eliminate the K-th entry in all rows that have not been used for pivoting.
%
    if ( imax ~= 0 )

      ipiv(imax) = k;
      a(imax,k+1:n) = a(imax,k+1:n) / a(imax,k);
      b(imax) = b(imax) / a(imax,k);
      a(imax,k) = 1.0;

      for i = 1 : n

        if ( ipiv(i) == 0 )
          a(i,k+1:n) = a(i,k+1:n) - a(i,k) * a(imax,k+1:n);
          b(i) = b(i) - a(i,k) * b(imax);
          a(i,k) = 0.0;
        end

      end
    end
  end
%
%  Now, every row with nonzero IPIV begins with a 1, and
%  all other rows are all zero.  Begin solution.
%
  for j = n : -1 : 1

    imax = 0;
    for k = 1 : n
      if ( ipiv(k) == j )
        imax = k;
      end
    end

    if ( imax == 0 )

      x(j) = 0.0;

      if ( b(j) == 0.0 )
        ierror = 1;
        fprintf ( 1, '\n' );
        fprintf ( 1, 'R8MAT_SOLVE2 - Warning:\n' );
        fprintf ( 1, '  Consistent singularity, equation = %d\n', j );
      else
        ierror = 2;
        fprintf ( 1, '\n' );
        fprintf ( 1, 'R8MAT_SOLVE2 - Error:\n' );
        fprintf ( 1, '  Inconsistent singularity, equation = %d\n', j );
      end

    else

      x(j) = b(imax);

      for i = 1 : n
        if ( i ~= imax )
          b(i) = b(i) - a(i,j) * x(j);
        end
      end

    end
  end

  return
end


text:
function verify_caffe_crops()

CONTEXT_PAD = 16;
CROP_SIZE = 227;
opts.data_mean = '~/working/caffe-rbg/matlab/caffe/ilsvrc_2012_mean.mat';
d = load(opts.data_mean);
off = floor((256 - 227)/2)+1;
IMAGE_MEAN = d.image_mean(off:off+227-1, off:off+227-1, :);
clear d;
%IMAGE_MEAN(:) = 0;

s = dir('dump/*.txt');
for i = 1:2:length(s)

  info_file = s(i+1).name;
  data_file = s(i).name;
  info = textread(['dump/' info_file], '%s');
  im_path = info{1};
  x1 = str2num(info{2});
  y1 = str2num(info{3});
  x2 = str2num(info{4});
  y2 = str2num(info{5});
  flip = str2num(info{6});
  label = str2num(info{7});
  is_fg = str2num(info{8});
  ov = 0; %str2num(info{9});

  %M1 = textread(['dump/' data_file], '%f');
  fid = fopen(['dump/' data_file], 'rb');
  M1 = fread(fid, 227*227*3, 'single');
  fclose(fid);

  if length(M1) ~= 227*227*3
    warning('skipping wrong size image');
    continue;
  end

  im1 = reshape(M1, [227 227 3]);
  im1 = permute(im1, [2 1 3]);
  im1 = im1(:,:,[3 2 1]);
  im1 = (im1 - min(im1(:))) / (max(im1(:)) - min(im1(:)));
  subplot(1,2,1);
  imagesc(im1);
  axis image;
  title(sprintf('max %.2f min %.2f  label %d  ov %.3f', max(M1), min(M1), label, ov));

  im = imread(im_path);
  im = single(im(:,:,[3 2 1]));
  if CONTEXT_PAD > 0
    scale = CROP_SIZE/(CROP_SIZE - CONTEXT_PAD*2);
    bbox = [x1 y1 x2 y2];
    %figure(1); showboxesc(im/256, bbox, 'b', '-');
    height = bbox(4)-bbox(2)+1;
    width = bbox(3)-bbox(1)+1;
    center = [bbox(1)+width/2 bbox(2)+height/2];
    bbox = round([center center] + [-width -height width height]/2*scale);
    height = bbox(4)-bbox(2)+1;
    width = bbox(3)-bbox(1)+1;
    %figure(1); showboxesc([], bbox, 'r', '-');
    pad_x1 = max(0, 1 - bbox(1));
    pad_y1 = max(0, 1 - bbox(2));
    pad_x2 = max(0, bbox(3) - size(im,2));
    pad_y2 = max(0, bbox(4) - size(im,1));
    bbox(1) = max(1, bbox(1));
    bbox(2) = max(1, bbox(2));
    bbox(3) = min(size(im,2), bbox(3));
    bbox(4) = min(size(im,1), bbox(4));
    window = im(bbox(2):bbox(4), bbox(1):bbox(3), :);

    window_width = round((bbox(3)-bbox(1)+1)*CROP_SIZE/width);
    window_height = round((bbox(4)-bbox(2)+1)*CROP_SIZE/height);
    pad_x1 = round(pad_x1*CROP_SIZE/width);
    pad_y1 = round(pad_y1*CROP_SIZE/height);
    pad_x2 = round(pad_x2*CROP_SIZE/width);
    pad_y2 = round(pad_y2*CROP_SIZE/height);

    pad_h = pad_y1;
    pad_w = pad_x1;
    if flip
      pad_w = pad_x2;
    end

    if pad_h + window_height > CROP_SIZE
      window_height = CROP_SIZE - pad_h;
    end
    if pad_w + window_width > CROP_SIZE
      window_width = CROP_SIZE - pad_w;
    end
    tmp = imresize(window, [window_height window_width], 'bilinear', 'antialiasing', false);
    if flip
      tmp = tmp(:, end:-1:1, :);
    end
    tmp = tmp - IMAGE_MEAN(1+pad_h:window_height+pad_h, 1+pad_w:window_width+pad_w, :);
    %figure(2); window_ = tmp; imagesc((window_-min(window_(:)))/(max(window_(:))-min(window_(:)))); axis image;
    window = zeros(CROP_SIZE, CROP_SIZE, 3, 'single');
    window(1+pad_h:window_height+pad_h, 1+pad_w:window_width+pad_w, :) = tmp;
    %figure(3); imagesc((window-min(window(:)))/(max(window(:))-min(window(:)))); axis image; pause;
    window = permute(window, [2 1 3]);
  else
    window = im(y1:y2, x1:x2, :);
    fprintf('im: %s\n', im_path);
    fprintf('box: %d %d %d %d\n', x1, y1, x2, y2);
    window = imresize(window, [CROP_SIZE CROP_SIZE], 'bilinear', 'Antialiasing', false);
    if flip
      window = window(:, end:-1:1, :);
    end
    window = window - IMAGE_MEAN;
    % permute to make width the fastest dimension (for caffe)
    window = permute(window, [2 1 3]);
  end

  M2 = window(:);

  im2 = reshape(M2, [227 227 3]);
  im2 = permute(im2, [2 1 3]);
  im2 = im2(:,:,[3 2 1]);
  im2 = (im2 - min(im2(:))) / (max(im2(:)) - min(im2(:)));
  subplot(1,2,2);
  imagesc(im2);
  axis image;
  title(sprintf('max %.3f min %.3f', max(M2), min(M2)));

  fprintf('max diff: %f\n', max(abs(M1 - M2)));

  pause;
end

author: rbgirshick

repo: rcnn

save_path: github-repos/MATLAB/rbgirshick-rcnn


text:
classdef FourthOrderVoigtPlaneStressRotator < FourthOrderVoigtRotator

    methods (Access = public)

        function obj = FourthOrderVoigtPlaneStressRotator(angle,dir)
            obj.compute(angle,dir)
        end

    end

    methods (Access = protected,Static)

        function r = createStressRotator(a,d)
            r = StressVoigtPlaneStressRotator(a,d);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%THIS FUMCTION ILLUSTRATES HOW TO USE THE USER-DEFINED ALGORITHM
function [AH,XH] = user_alg4(Y,r,X)
%
% The example of implementing the user-defined algorithm (this is the Lee-Seung algorithm based on the KL divergence)
%
% INPUTS:
% Y - mixed signals (matrix of size [m by T])
% r - number of estimated signals
% X - true source signals
%
% OUTPUTS
% AH - estimated mixing matrix (matrix of size [m by r])
% XH - estimated source signals (matrix of size [r by T])
%
% #########################################################################
% Initialization
[m,T]=size(Y);
Y(Y <=0) = eps; % this enforces the positive value in the data
AH=rand(m,r);
XH=rand(r,T);

IterNo = 1000; % number of alternating steps

% Iterations
for k = 1:IterNo
    XH = XH.*(AH'*(Y./(AH*XH + eps)));
    AH = AH.*((Y./(AH*XH + eps))*XH')./repmat(sum(XH,2)',m,1);
    AH = AH*diag(1./(sum(AH,1) + eps));
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function a = hadamard ( m, n )

%*****************************************************************************80
%
%% HADAMARD returns a HADAMARD matrix.
%
%  Definition:
%
%    A Hadamard matrix is a square matrix A of order N, whose entries are
%    only +1's or -1's, with the property that:
%
%      A * A' = N * I.
%
%  Notes:
%
%    A Hadamard matrix must be of order 1, 2, or else a multiple of 4.
%    It is not known whether a Hadamard matrix exists for every multiple
%    of 4.
%
%    The method used here allows the user to request a Hadamard matrix
%    of any rectangular order, M by N.  The algorithm then essentially
%    finds the largest powers of 2 that are less than or equal to M and
%    N, and produces a Hadamard-like matrix in that space, setting the
%    rest of the matrix to 0.  Thus, the matrix returned by this routine
%    is only a Hadamard matrix if M = N = a power of 2.
%
%  Formula:
%
%    The following recursive formula is used to produce a series of
%    Hadamard matrices of increasing size.
%
%    H(0) = [1]
%
%    H(1) = [ H(0)  H(0) ] = [ 1  1]
%           [ H(0) -H(0) ]   [ 1 -1]
%
%    H(2) = [ H(1)  H(1) ] = [ 1  1  1  1]
%           [ H(1) -H(1) ]   [ 1 -1  1 -1]
%                            [ 1  1 -1 -1]
%                            [ 1 -1 -1  1]
%
%    and so on.
%
%  Properties:
%
%    All entries of a Hadamard matrix are either +1 or -1.  Matrices
%    produced by this routine will be +1 or -1 up to a certain row
%    and column, beyond which the entries will be zero.
%
%    The Hadamard matrices produced by this routine have the property
%    that the first row and column are entirely 1's, although this
%    is not a requirement for a Hadamard matrix.
%
%    The matrices produced by this algorithm are (loosely) symmetric,
%    although that is not required for a Hadamard matrix.
%
%    Hadamard matrices exhibit the maximum possible relative growth of pivot
%    elements during Gaussian elimination with complete pivoting.
%
%    The inverse of a Hadamard matrix of order N is A itself,
%    scaled by 1.0/N.  Thus 1.0/sqrt(N) times a Hadamard matrix
%    yields a symmetric matrix which is its own inverse, or
%    "involutional".
%
%    A is integral: int ( A ) = A.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Gregory and David Karney,
%    Example 3.14,
%    A Collection of Matrices for Testing Computational Algorithms,
%    Wiley, New York, 1969, page 42,
%    LC: QA263.G68.
%
%    William Pratt,
%    Digital Image Processing,
%    John Wiley and Sons, 1978.
%
%    Herbert Ryser,
%    Combinatorial Mathematics,
%    John Wiley and Sons, 1963.
%
%  Parameters:
%
%    Input, integer M, the row order of A.
%
%    Input, integer N, the column order of A.
%
%    Output, real A(M,N), the matrix.
%
  a = zeros ( m, n );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HADAMARD - Fatal error!\n' );
    fprintf ( 1, '  Input value of M = %d\n', m );
    fprintf ( 1, '  but M must be positive.\n' );
    error ( 'HADAMARD - Fatal error!' );
  end

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HADAMARD - Fatal error!\n' );
    fprintf ( 1, '  Input value of N = %d\n', n );
    fprintf ( 1, '  but N must be positive.\n' );
    error ( 'HADAMARD - Fatal error!' );
  end

  a(1,1) = 1.0;

  nn = 1;

  while ( nn < n | nn < m )

    for i = 1 : nn
      for j = 1 : nn

        if ( i <= m & j+nn <= n )
          if ( 2 * nn <= n )
            a(i,j+nn) = a(i,j);
          else
            a(i,j+nn) = 0.0;
          end
        end

        if ( i + nn <= m & j <= n )
          if ( 2 * nn <= m )
            a(i+nn,j) = a(i,j);
          else
            a(i+nn,j) = 0.0;
          end
        end

        if ( i + nn <= m & j + nn <= n )
          if ( 2 * nn <= m & 2 * nn <= n )
            a(i+nn,j+nn) = - a(i,j);
          else
            a(i+nn,j+nn) = 0.0;
          end
        end

      end
    end

    nn = 2 * nn;

  end

  return
end


text:
function y=bvpme(t,x,flag,U)
y =[0.5*x(1)*(x(3)-x(1))/x(2);
    -0.5*(x(3)-x(1));
    (0.9-1000*(x(3)-x(5))-0.5*x(3)*(x(3)-x(1)))/x(4);
    0.5*(x(3)-x(1));
    -100*(x(5)-x(3));];
end

author: Sable

repo: mcbench-benchmarks


text:
function b = bernoulli_number ( n )

%*****************************************************************************80
%
%% BERNOULLI_NUMBER computes the value of the Bernoulli numbers B(0) through B(N).
%
%  Discussion:
%
%    The Bernoulli numbers are rational.
%
%    If we define the sum of the M-th powers of the first N integers as:
%
%      SIGMA(M,N) = sum ( 0 <= I <= N ) I**M
%
%    and let C(I,J) be the combinatorial coefficient:
%
%      C(I,J) = I! / ( ( I - J )! * J! )
%
%    then the Bernoulli numbers B(J) satisfy:
%
%      SIGMA(M,N) = 1/(M+1) * sum ( 0 <= J <= M ) C(M+1,J) B(J) * (N+1)**(M+1-J)
%
%  First values:
%
%   B0  1                   =         1.00000000000
%   B1 -1/2                 =        -0.50000000000
%   B2  1/6                 =         1.66666666666
%   B3  0                   =         0
%   B4 -1/30                =        -0.03333333333
%   B5  0                   =         0
%   B6  1/42                =         0.02380952380
%   B7  0                   =         0
%   B8 -1/30                =        -0.03333333333
%   B9  0                   =         0
%  B10  5/66                =         0.07575757575
%  B11  0                   =         0
%  B12 -691/2730            =        -0.25311355311
%  B13  0                   =         0
%  B14  7/6                 =         1.16666666666
%  B15  0                   =         0
%  B16 -3617/510            =        -7.09215686274
%  B17  0                   =         0
%  B18  43867/798           =        54.97117794486
%  B19  0                   =         0
%  B20 -174611/330          =      -529.12424242424
%  B21  0                   =         0
%  B22  854,513/138         =      6192.123
%  B23  0                   =         0
%  B24 -236364091/2730      =    -86580.257
%  B25  0                   =         0
%  B26  8553103/6           =   1425517.16666
%  B27  0                   =         0
%  B28 -23749461029/870     = -27298231.0678
%  B29  0                   =         0
%  B30  8615841276005/14322 = 601580873.901
%
%  Recursion:
%
%    With C(N+1,K) denoting the standard binomial coefficient,
%
%    B(0) = 1.0
%    B(N) = - ( sum ( 0 <= K < N ) C(N+1,K) * B(K) ) / C(N+1,N)
%
%  Warning:
%
%    This recursion, which is used in this routine, rapidly results
%    in significant errors.
%
%  Special Values:
%
%    Except for B(1), all Bernoulli numbers of odd index are 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 December 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the highest Bernoulli number to compute.
%
%    Output, real B(1:N+1), B(I+1) contains the I-th Bernoulli number.
%
  if ( n < 0 )
    b = [];
    return
  end

  b(1) = 1.0;

  if ( n < 1 )
    return
  end

  b(2) = -0.5;

  c = zeros ( n + 2, 1 );

  c(1) = 1;
  c(2) = 2;
  c(3) = 1;

  for i = 2 : n

    c = comb_row_next ( i + 1, c );

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

      b(i+1) = 0.0;

    else

      b_sum = 0.0;
      for j = 0 : i - 1
        b_sum = b_sum + b(j+1) * c(j+1);
      end

      b(i+1) = -b_sum / c(i+1);

    end

  end

  return
end


text:
% DEMO  --  Chebyshev Polynomial Integral
% UPDATED  --  October 28, 2013
% Written by Matthew Kelly, Cornell University
%
% This demo is designed to show how to compute integrals of chebyshev
% polynomials.
%
% Results - Note that there are relative large errors between the chebyshev
% approximation of the integral and the analytic integral, BUT the last
% value (which is the value of the definite integral) is very accurate.
%
% The last point is so accurate because it is the solution to the definite
% integral over the domain using Clenshaw-Curtis quadrature for an order
% "order" accurate polynomial.
%

%% General Settings
clear; clc;

%What order should the approximation be?
order = 50;

%What domain should we be looking at?
d = [0,1];

%Get chebyshev values at the nodes:
[x,w] = chebyshevPoints(order+1,d);
[~,f] = testFunction(x);
If= chebyshevIntegral(f,d);

%Set the time vector
time = linspace(d(1),d(2),10000);
[Ig, g] = testFunction(time);

%Adjust the constants of integration S.T. they all start at zero:
Ig = Ig-Ig(1);

%Interpolate the chebyshev approximations:
y = chebyshevInterpolate(f,time,d);
Iy = chebyshevInterpolate(If,time,d);

%Show results
figure(407); clf;
subplot(2,2,1); hold on;
    plot(time,g,'b-','LineWidth',2)
    plot(x,f,'ko','MarkerSize',5);
    title(['function approximation - order ' num2str(order)]);
    legend('Analytic','Cheb. nodes','Location','NorthWest')
subplot(2,2,3); hold on;
    plot(time,Ig,'b-','LineWidth',2)
    plot(x,If,'ko','MarkerSize',5);
    title(['integral approximation - order ' num2str(order)]);
        legend('Analytic','Cheb. nodes','Location','NorthWest')
subplot(2,2,2);
    semilogy(time,abs(g-y))
    title('error in function');
subplot(2,2,4); hold on;
    semilogy(time,abs(Ig-Iy))
    semilogy(time(end),abs(Ig(end)-Iy(end)),...
        'bo','MarkerSize',10,'LineWidth',2)
    title(['error in definite integral: ' num2str(Ig(end)-Iy(end))]);
    set(gca,'yscale','log')

author: MatthewPeterKelly

repo: dscTutorials


text:
% GENERATE BILINEAR TIME-FREQUENCY TRANSFORMATIONS USING THE DIRECT METHOD
%
% BILINEAR TRANSFORMATIONS transform the time domain signal (the
% variable INPUT SIGNAL) to an output time-frequency distribution (the
% variable OUTPUT TIME-FREQUENCY ARRAY)
%
% These functions use a direct implementation rather than using the
% quadratic time-frequency implementation. This results in a
% computationally optimised routine.
%
% The different types of time-frequency distributions that can be
% implemented directly are:
%
%   (i)    Wigner-Ville Distribution
%   (ii)   Short-Time Fourier Transform
%   (iii)  Short-Time Fourier Transform(overlap)
%   (iv)   Rihaczek
%   (v)    Windowed-Rihaczek
%
%   The various parameters associated with the distributions are as follows:
%
%   TIME-FREQUENCY ARRAY (tfd)
%
%      The computed time-frequency distribution.  size(tfd) will
%      return [a, b], where a is the next largest power of two above
%      FFT length, and b is floor(length(signal)/time_res) - 1.
%
%   INPUT SIGNAL
%
%      Input one dimensional signal to be analysed. An analytic signal
%      is required for this function, however, if signal is real, a
%      default analytic transformer routine will be called from this
%      function before computing tfd.
%
%   TIME RESOLUTION
%
%      The number of time samples to skip between successive time slices.
%
%   LAG WINDOW LENGTH
%
%      This is the lag window length and controls the size of the
%      signal kernel (or instantaneous autocorrelation function) used
%      for analysis (lag_window_length must be odd). The kernel used
%      will be defined from -(lag_window_length+1)/2 to
%      +(lag_window_length+1)/2 in both time and lag dimensions.
%
%   FFT Length:
%
%      Zero-padding at the FFT stage of the analysis may be specified
%      by giving an FFT length larger than lag window length.  If
%      FFT length is not specified, or is smaller than the
%      lag window length, then the next highest power of two above
%      lag window length is used.  If FFT length is not a power of
%      two, the next highest power of two is used.
%
%
%
%  See Also:  wvd, spec, rihaczek, analyt

author: HeLiangHIT


text:
function [rOUT,sOUT] = FindLocalCoords2D(k, xout, yout)

% function [rOUT,sOUT] = FindLocalCoords2D(k, xout, yout)
% Purpose: find local (r,s) coordinates in the k'th element of given coordinates
%          [only works for straight sided triangles]

Globals2D;

v1  = EToV(k,1); v2  = EToV(k,2); v3  = EToV(k,3);
xy1 = [VX(v1);VY(v1)];  xy2 = [VX(v2);VY(v2)];  xy3 = [VX(v3);VY(v3)];

A(:,1) = xy2-xy1;  A(:,2) = xy3-xy1;

rOUT = zeros(size(xout));  sOUT = zeros(size(xout));

for i=1:length(xout(:))
  rhs = 2.0*[xout(i);yout(i)] -xy2 -xy3;

  tmp = A\rhs;
  rOUT(i) = tmp(1); sOUT(i) = tmp(2);
end
return

repo: Aero-matlab


text:
function B = design_fir(N,F,A,nfft,W)
% B = design_fir(N,F,A,nFFT,W)
% Design an FIR filter using the frequency-sampling method.
%
% The frequency response is interpolated cubically between the specified
% frequency points.
%
% In:
%   N : order of the filter
%
%   F : vector of frequencies at which amplitudes shall be defined
%       (starts with 0 and goes up to 1; try to avoid too
%        sharp transitions)
%
%   A : vector of amplitudes, one value per specified frequency
%
%   nFFT : optionally number of FFT bins to use
%
%   W : optionally the window function to use (default: Hamming)
%
% Out:
%   B : designed filter kernel
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2013-08-14

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

if nargin < 4 || isempty(nfft)
    nfft = max(512,2^ceil(log(N)/log(2))); end
if nargin < 5
    W = 0.54 - 0.46*cos(2*pi*(0:N)/N); end

% calculate interpolated frequency response
F = interp1(round(F*nfft),A,(0:nfft),'pchip');

% set phase & transform into time domain
F = F .* exp(-(0.5*N)*sqrt(-1)*pi*(0:nfft)./nfft);
B = real(ifft([F conj(F(end-1:-1:2))]));

% apply window to kernel
B = B(1:N+1).*W(:)';

author: VisLab

repo: EEG-Clean-Tools


text:
function str=getSummaryStr(obj, format)

    formatters=struct();
    formatters.text=@obj2text;
    formatters.xml=@obj2xml;

    f=formatters.(format);

    str=f(obj);

function str=obj2text(obj)
    content=getSummaryContent(obj);

    % set outcome to 'failure', 'error', 'skipped', or 'passed'
    outcome_verbosity=2;
    outcome=getOutcomeStr(obj, outcome_verbosity);

    if isstruct(content)
        % error or failure
        stack_prefix='  ';
        stack_str=moxunit_util_stack2str(content.stack, stack_prefix);
        str=sprintf('%s: %s\n%s',...
                outcome,content.message,stack_str);

    elseif ischar(content)
        % skipped
        str=sprintf('%s: %s', outcome, content);

    elseif isequal(content,[])
        % passed
        str='';

    else
        assert(false,'this should not happen');
    end

function str=obj2xml(obj)
    test_=getTest(obj);

    prefix=sprintf('<testcase classname="%s" name="%s" time="%.3f"',...
                moxunit_util_escape_xml(get_classname(test_)),...
                moxunit_util_escape_xml(getName(test_)),...
                getDuration(obj));

    if isSuccess(obj)
        infix='';
        suffix=' />';
    else
        content=getSummaryContent(obj);

        % outcome is 'failure', 'error', or 'skipped'
        outcome_verbosity=2;
        outcome=getOutcomeStr(obj, outcome_verbosity);

        if isNonFailure(obj)
            % skipped
            % note: showing the reason is not supported
            infix=sprintf('>\n  <%s />',outcome);
        else
            % XXX currently the entire error message is shown, together
            % with the stack trace. It may be more convenient to show a
            % summary of this information
            raw_message=content.message;
            message=moxunit_util_remove_matlab_anchor_tag(raw_message);
            message=moxunit_util_escape_xml(message);

            raw_stack_str=moxunit_util_stack2str(content.stack);
            stack_str=moxunit_util_escape_xml(raw_stack_str);

            infix=sprintf('>\n  <%s message="%s">%s</%s>',...
                        outcome,message,stack_str,outcome);
        end

        suffix=sprintf('\n</testcase>');
    end

    str=sprintf('%s%s%s',prefix,infix,suffix);

function s=get_classname(test_)
    location=getLocation(test_);
    [pth,file_name]=fileparts(location);

    s=sprintf('%s',file_name);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function s = getStatisticsHeaderSizeWordAligned(this)

% Assume a word is 4 bytes

s = getStatisticsHeaderSize(this);
s = ceil(s/4)*4;


text:
%-----------------------------------------------------------------------
% Job configuration created by cfg_util (rev $Rev: 325M $)
%-----------------------------------------------------------------------
matlabbatch{1}.cfg_basicio.cfg_named_dir.name = 'Subject directory';
matlabbatch{1}.cfg_basicio.cfg_named_dir.dirs = {'<UNDEFINED>'};
matlabbatch{2}.cfg_basicio.cfg_cd.dir = '<UNDEFINED>';
matlabbatch{3}.cfg_basicio.cfg_mkdir.parent = '<UNDEFINED>';
matlabbatch{3}.cfg_basicio.cfg_mkdir.name = 'categorical';
matlabbatch{4}.spm.spatial.realign.estwrite.data = {'<UNDEFINED>'};
matlabbatch{4}.spm.spatial.realign.estwrite.eoptions.quality = 0.9;
matlabbatch{4}.spm.spatial.realign.estwrite.eoptions.sep = 4;
matlabbatch{4}.spm.spatial.realign.estwrite.eoptions.fwhm = 5;
matlabbatch{4}.spm.spatial.realign.estwrite.eoptions.rtm = 0;
matlabbatch{4}.spm.spatial.realign.estwrite.eoptions.interp = 2;
matlabbatch{4}.spm.spatial.realign.estwrite.eoptions.wrap = [0 0 0];
matlabbatch{4}.spm.spatial.realign.estwrite.eoptions.weight = {};
matlabbatch{4}.spm.spatial.realign.estwrite.roptions.which = [2 1];
matlabbatch{4}.spm.spatial.realign.estwrite.roptions.interp = 4;
matlabbatch{4}.spm.spatial.realign.estwrite.roptions.wrap = [0 0 0];
matlabbatch{4}.spm.spatial.realign.estwrite.roptions.mask = 1;
matlabbatch{4}.spm.spatial.realign.estwrite.roptions.prefix = 'r';
matlabbatch{5}.spm.temporal.st.scans = {'<UNDEFINED>'};
matlabbatch{5}.spm.temporal.st.nslices = 24;
matlabbatch{5}.spm.temporal.st.tr = 2;
matlabbatch{5}.spm.temporal.st.ta = 1.91666666666667;
matlabbatch{5}.spm.temporal.st.so = [24 23 22 21 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1];
matlabbatch{5}.spm.temporal.st.refslice = 12;
matlabbatch{5}.spm.temporal.st.prefix = 'a';
matlabbatch{6}.spm.spatial.coreg.estimate.ref = '<UNDEFINED>';
matlabbatch{6}.spm.spatial.coreg.estimate.source = '<UNDEFINED>';
matlabbatch{6}.spm.spatial.coreg.estimate.other = {''};
matlabbatch{6}.spm.spatial.coreg.estimate.eoptions.cost_fun = 'nmi';
matlabbatch{6}.spm.spatial.coreg.estimate.eoptions.sep = [4 2];
matlabbatch{6}.spm.spatial.coreg.estimate.eoptions.tol = [0.02 0.02 0.02 0.001 0.001 0.001 0.01 0.01 0.01 0.001 0.001 0.001];
matlabbatch{6}.spm.spatial.coreg.estimate.eoptions.fwhm = [7 7];
matlabbatch{7}.spm.spatial.preproc.data = '<UNDEFINED>';
matlabbatch{7}.spm.spatial.preproc.output.GM = [0 0 1];
matlabbatch{7}.spm.spatial.preproc.output.WM = [0 0 1];
matlabbatch{7}.spm.spatial.preproc.output.CSF = [0 0 0];
matlabbatch{7}.spm.spatial.preproc.output.biascor = 1;
matlabbatch{7}.spm.spatial.preproc.output.cleanup = 0;
matlabbatch{7}.spm.spatial.preproc.opts.tpm = {
                                               '/export/spm-devel/spm/trunk/tpm/grey.nii'
                                               '/export/spm-devel/spm/trunk/tpm/white.nii'
                                               '/export/spm-devel/spm/trunk/tpm/csf.nii'
                                               };
matlabbatch{7}.spm.spatial.preproc.opts.ngaus = [2
                                                 2
                                                 2
                                                 4];
matlabbatch{7}.spm.spatial.preproc.opts.regtype = 'mni';
matlabbatch{7}.spm.spatial.preproc.opts.warpreg = 1;
matlabbatch{7}.spm.spatial.preproc.opts.warpco = 25;
matlabbatch{7}.spm.spatial.preproc.opts.biasreg = 0.0001;
matlabbatch{7}.spm.spatial.preproc.opts.biasfwhm = 60;
matlabbatch{7}.spm.spatial.preproc.opts.samp = 3;
matlabbatch{7}.spm.spatial.preproc.opts.msk = {''};
matlabbatch{8}.spm.spatial.normalise.write.subj.matname = '<UNDEFINED>';
matlabbatch{8}.spm.spatial.normalise.write.subj.resample = '<UNDEFINED>';
matlabbatch{8}.spm.spatial.normalise.write.roptions.preserve = 0;
matlabbatch{8}.spm.spatial.normalise.write.roptions.bb = [-78 -112 -50
                                                          78 76 85];
matlabbatch{8}.spm.spatial.normalise.write.roptions.vox = [3 3 3];
matlabbatch{8}.spm.spatial.normalise.write.roptions.interp = 1;
matlabbatch{8}.spm.spatial.normalise.write.roptions.wrap = [0 0 0];
matlabbatch{8}.spm.spatial.normalise.write.roptions.prefix = 'w';
matlabbatch{9}.spm.spatial.smooth.data = '<UNDEFINED>';
matlabbatch{9}.spm.spatial.smooth.fwhm = [8 8 8];
matlabbatch{9}.spm.spatial.smooth.dtype = 0;
matlabbatch{9}.spm.spatial.smooth.prefix = 's';
matlabbatch{10}.spm.stats.fmri_spec.dir = '<UNDEFINED>';
matlabbatch{10}.spm.stats.fmri_spec.timing.units = 'scans';
matlabbatch{10}.spm.stats.fmri_spec.timing.RT = 2;
matlabbatch{10}.spm.stats.fmri_spec.timing.fmri_t = 24;
matlabbatch{10}.spm.stats.fmri_spec.timing.fmri_t0 = 12;
matlabbatch{10}.spm.stats.fmri_spec.sess.scans = '<UNDEFINED>';
matlabbatch{10}.spm.stats.fmri_spec.sess.cond = struct('name', {}, 'onset', {}, 'duration', {}, 'tmod', {}, 'pmod', {});
matlabbatch{10}.spm.stats.fmri_spec.sess.multi = '<UNDEFINED>';
matlabbatch{10}.spm.stats.fmri_spec.sess.regress = struct('name', {}, 'val', {});
matlabbatch{10}.spm.stats.fmri_spec.sess.multi_reg = '<UNDEFINED>';
matlabbatch{10}.spm.stats.fmri_spec.sess.hpf = 128;
matlabbatch{10}.spm.stats.fmri_spec.fact(1).name = 'Fam';
matlabbatch{10}.spm.stats.fmri_spec.fact(1).levels = 2;
matlabbatch{10}.spm.stats.fmri_spec.fact(2).name = 'Rep';
matlabbatch{10}.spm.stats.fmri_spec.fact(2).levels = 2;
matlabbatch{10}.spm.stats.fmri_spec.bases.hrf.derivs = [1 1];
matlabbatch{10}.spm.stats.fmri_spec.volt = 1;
matlabbatch{10}.spm.stats.fmri_spec.global = 'None';
matlabbatch{10}.spm.stats.fmri_spec.mask = {''};
matlabbatch{10}.spm.stats.fmri_spec.cvi = 'AR(1)';
matlabbatch{11}.spm.stats.fmri_est.spmmat = '<UNDEFINED>';
matlabbatch{11}.spm.stats.fmri_est.method.Classical = 1;
matlabbatch{12}.spm.stats.con.spmmat = '<UNDEFINED>';
matlabbatch{12}.spm.stats.con.consess{1}.fcon.name = 'Effects of interest';
%%
matlabbatch{12}.spm.stats.con.consess{1}.fcon.convec = {
                                                        [1 0 0 0 0 0 0 0 0 0 0 0
                                                        0 1 0 0 0 0 0 0 0 0 0 0
                                                        0 0 1 0 0 0 0 0 0 0 0 0
                                                        0 0 0 1 0 0 0 0 0 0 0 0
                                                        0 0 0 0 1 0 0 0 0 0 0 0
                                                        0 0 0 0 0 1 0 0 0 0 0 0
                                                        0 0 0 0 0 0 1 0 0 0 0 0
                                                        0 0 0 0 0 0 0 1 0 0 0 0
                                                        0 0 0 0 0 0 0 0 1 0 0 0
                                                        0 0 0 0 0 0 0 0 0 1 0 0
                                                        0 0 0 0 0 0 0 0 0 0 1 0
                                                        0 0 0 0 0 0 0 0 0 0 0 1]
                                                        }';
%%
matlabbatch{12}.spm.stats.con.consess{1}.fcon.sessrep = 'none';
matlabbatch{12}.spm.stats.con.delete = 0;
matlabbatch{13}.spm.stats.results.spmmat = '<UNDEFINED>';
matlabbatch{13}.spm.stats.results.conspec.titlestr = '';
matlabbatch{13}.spm.stats.results.conspec.contrasts = Inf;
matlabbatch{13}.spm.stats.results.conspec.threshdesc = 'FWE';
matlabbatch{13}.spm.stats.results.conspec.thresh = 0.05;
matlabbatch{13}.spm.stats.results.conspec.extent = 0;
matlabbatch{13}.spm.stats.results.conspec.mask = struct('contrasts', {}, 'thresh', {}, 'mtype', {});
matlabbatch{13}.spm.stats.results.print = 1;

author: spm


text:
classdef GFTTDetector < handle
    %GFTTDETECTOR  Wrapping class for feature detection using the goodFeaturesToTrack function
    %
    % The function finds the most prominent corners in the image or in the
    % specified image region, as described in [Shi94]:
    %
    % 1. Function calculates the corner quality measure at every source image
    %    pixel using the cv.cornerMinEigenVal or cv.cornerHarris.
    % 2. Function performs a non-maximum suppression (the local maximums in
    %    `3x3` neighborhood are retained).
    % 3. The corners with the minimal eigenvalue less than
    %    `QualityLevel * max_{x,y}(qualityMeasureMap(x,y))` are rejected.
    % 4. The remaining corners are sorted by the quality measure in the
    %    descending order.
    % 5. Function throws away each corner for which there is a stronger
    %    corner at a distance less than `maxDistance`.
    %
    % ## References
    % [Shi94]:
    % > Jianbo Shi and Carlo Tomasi. "Good features to track".
    % > In Computer Vision and Pattern Recognition, 1994. Proceedings CVPR'94.,
    % > 1994 IEEE Computer Society Conference on, pages 593-600. IEEE, 1994.
    %
    % See also: cv.goodFeaturesToTrack, cv.cornerHarris, cv.cornerMinEigenVal,
    %  cv.FeatureDetector
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % Maximum number of corners to return.
        %
        % If there are more corners than are found, the strongest of them is
        % returned. Default 1000
        MaxFeatures
        % Parameter characterizing the minimal accepted quality of image
        % corners.
        %
        % The parameter value is multiplied by the best corner quality
        % which is the minimal eigenvalue (see cv.cornerMinEigenVal) or the
        % Harris function response (see cv.cornerHarris). The corners with the
        % quality measure less than the product are rejected. For example, if
        % the best corner has the quality measure = 1500, and the
        % `QualityLevel=0.01`, then all the corners with the quality measure
        % less than 15 are rejected. Default 0.01
        QualityLevel
        % Minimum possible Euclidean distance between the returned corners.
        %
        % Default 1.0
        MinDistance
        % Size of an average block for computing a derivative covariation
        % matrix over each pixel neighborhood.
        %
        % See cv.cornerEigenValsAndVecs. Default 3
        BlockSize
        % Parameter indicating whether to use a Harris detector
        % (see cv.cornerHarris) or cv.cornerMinEigenVal.
        %
        % Default false
        HarrisDetector
        % Free parameter of the Harris detector.
        %
        % Default 0.04
        K
    end

    methods
        function this = GFTTDetector(varargin)
            %GFTTDETECTOR  Constructor
            %
            %     obj = cv.GFTTDetector()
            %     obj = cv.GFTTDetector(..., 'OptionName',optionValue, ...)
            %
            % ## Options
            % * __MaxFeatures__ See cv.GFTTDetector.MaxFeatures, default 1000
            % * __QualityLevel__ See cv.GFTTDetector.QualityLevel,
            %   default 0.01
            % * __MinDistance__ See cv.GFTTDetector.MinDistance, default 1
            % * __BlockSize__ See cv.GFTTDetector.BlockSize, default 3
            % * __GradientSize__ default 3
            % * __HarrisDetector__ See cv.GFTTDetector.HarrisDetector,
            %   default false
            % * __K__ See cv.GFTTDetector.K, default 0.04
            %
            % See also: cv.GFTTDetector.detect
            %
            this.id = GFTTDetector_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.GFTTDetector
            %
            if isempty(this.id), return; end
            GFTTDetector_(this.id, 'delete');
        end

        function typename = typeid(this)
            %TYPEID  Name of the C++ type (RTTI)
            %
            %     typename = obj.typeid()
            %
            % ## Output
            % * __typename__ Name of C++ type
            %
            typename = GFTTDetector_(this.id, 'typeid');
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.GFTTDetector.empty, cv.GFTTDetector.load
            %
            GFTTDetector_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Checks if detector object is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the detector object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.GFTTDetector.clear, cv.GFTTDetector.load
            %
            b = GFTTDetector_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.GFTTDetector.load
            %
            GFTTDetector_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.GFTTDetector.save
            %
            GFTTDetector_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.GFTTDetector.save, cv.GFTTDetector.load
            %
            name = GFTTDetector_(this.id, 'getDefaultName');
        end
    end

    %% Features2D: FeatureDetector
    methods
        function keypoints = detect(this, img, varargin)
            %DETECT  Detects keypoints in an image or image set
            %
            %     keypoints = obj.detect(img)
            %     keypoints = obj.detect(imgs)
            %     [...] = obj.detect(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __img__ Image (first variant), 8-bit grayscale image where
            %   keypoints (corners) are detected.
            % * __imgs__ Image set (second variant), cell array of images.
            %
            % ## Output
            % * __keypoints__ The detected keypoints. In the first variant, a
            %   1-by-N structure array. In the second variant of the method,
            %   `keypoints{i}` is a set of keypoints detected in `imgs{i}`.
            %
            % ## Options
            % * __Mask__ A mask specifying where to look for keypoints
            %   (optional). It must be a logical or 8-bit integer matrix with
            %   non-zero values in the region of interest. In the second
            %   variant, it is a cell-array of masks for each input image,
            %   `masks{i}` is a mask for `imgs{i}`. Not set by default.
            %
            % See also: cv.GFTTDetector.GFTTDetector
            %
            keypoints = GFTTDetector_(this.id, 'detect', img, varargin{:});
        end
    end

    %% Getters/Setters
    methods
        function value = get.MaxFeatures(this)
            value = GFTTDetector_(this.id, 'get', 'MaxFeatures');
        end
        function set.MaxFeatures(this, value)
            GFTTDetector_(this.id, 'set', 'MaxFeatures', value);
        end

        function value = get.QualityLevel(this)
            value = GFTTDetector_(this.id, 'get', 'QualityLevel');
        end
        function set.QualityLevel(this, value)
            GFTTDetector_(this.id, 'set', 'QualityLevel', value);
        end

        function value = get.MinDistance(this)
            value = GFTTDetector_(this.id, 'get', 'MinDistance');
        end
        function set.MinDistance(this, value)
            GFTTDetector_(this.id, 'set', 'MinDistance', value);
        end

        function value = get.BlockSize(this)
            value = GFTTDetector_(this.id, 'get', 'BlockSize');
        end
        function set.BlockSize(this, value)
            GFTTDetector_(this.id, 'set', 'BlockSize', value);
        end

        function value = get.HarrisDetector(this)
            value = GFTTDetector_(this.id, 'get', 'HarrisDetector');
        end
        function set.HarrisDetector(this, value)
            GFTTDetector_(this.id, 'set', 'HarrisDetector', value);
        end

        function value = get.K(this)
            value = GFTTDetector_(this.id, 'get', 'K');
        end
        function set.K(this, value)
            GFTTDetector_(this.id, 'set', 'K', value);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% A = riemann_mean(B,epsilon,tol)
%
% Calcul du barycentre des matrice de covariances.
% A : baricentre des K matrices NxN
%
% B : Matrice NxNxK
% epsilon : Pas de la descente de gradient
% tol : arret de la descente si le crit�re < tol

function [A critere niter] = opttransp_mean(B,args)
Imat = size(B,3);
N_itermax = 200;
if (nargin<2)||(isempty(args))
    tol = 10^-3;
    %A = mean(B,3);
    A = eye(size(B,1));
else
    tol = args{1};
    A = args{2};
end

niter = 0;
fc = 0;
K = A^(0.5);

while (niter<N_itermax)
    niter = niter+1;

    Ktmp = 0;
    for i=1:Imat
       Ktmp = Ktmp + (K*B(:,:,i)*K)^(0.5);
    end
    Ktmp = (Ktmp)^(0.5);

    fcn = norm(Ktmp-K,'fro');
    K = Ktmp;
    % improvement
    conv = abs((fcn-fc)/fc);
    if conv<tol % break if the improvement is below the tolerance
       break;
    end
    fc = fcn;
end

A = ((1/Imat)*K)^2;

if niter==N_itermax
    disp('Warning : Nombre d''itérations maximum atteint');
end

critere = fc;

author: alexandrebarachant

repo: covariancetoolbox


text:
% get_MBS_GCI_intervals.m
% Function to detect intervals to be used for searching the LP-residual to
% determine GCI locations using the method described in Drugman et al
% (2012).
%
% Octave compatible
%
% Description
% Function to detect intervals to be used for searching the LP-residual to
% determine GCI locations using the method described in Drugman et al
% (2012).
%
% Inputs
%  MBS             : [samples] [Nx1] Mean-based signal
%  fs              : [Hz]      [1x1] sampling frequency
%  T0mean          : [samples] [1x1] Glottal period
%  F0max           : [Hz]      [1x1] Maximum fundamental frequency
%
% Outputs
%  interval        : [samples] [Mx2] Search intervals for GCI estimation
%
% Example
%  interval = get_MBS_GCI_intervals(MBS,fs,T0mean,F0max)
%
% References
%  [1] Kane, J., Gobl, C., (2013) `Evaluation of glottal closure instant
%       detection in a range of voice qualities', Speech Communication
%       55(2), pp. 295-314.
%
% Copyright (c) 2013 Trinity College Dublin
%
% License
%  This code is a part of the Voice Analysis Toolkit with the following
%  licence:
%  The software product (and any modifications thereof) is distributed under
%  a dual licence: an open source license for individual, non-commercial
%  purposes, and a commercial license. The opensource licence under which
%  the product is distributed is GNU GPL v2. For individual users, this
%  licence suits their use as these are not distributing proprietary
%  modifications, additions to, or derivatives of the product and don't
%  require legal protection of a commercial licence. For commercial users,
%  where open source does not meet their requirements, we offer commercial
%  licensing of the product. A commercial license permits customers to
%  modify, add or produce derivative products without the obligation of
%  making the subsequent code open source. For more information regarding
%  our commercial licence, please contact john.whelan@tcd.ie
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  John Kane kanejo@tcd.ie
%
% $Id <info set by the versioning system> $

function interval = get_MBS_GCI_intervals(MBS,fs,T0mean,F0max)

% Function to detect intervals to be used for searching the LP-residual to
% determine GCI locations using the method described in Drugman et al
% (2012).

%% Initial settings
if nargin < 4
    F0max=500;
end
F0max=F0max*2;
T0max=round(fs/F0max);
[idx,TMP]=v_findpeaks(MBS*-1,'minpeakdistance',T0max); % Find locations of negative peaks
N=length(idx);
search_rate=0.28;
search_left_rate=0.01;
interval=zeros(N,2);

%% Do processing
for n=1:N

    if length(T0mean)>1
        start=idx(n)-round(T0mean(idx(n))*search_left_rate);
        stop=idx(n)+round(T0mean(idx(n))*search_rate);
    else start=idx(n)-round(T0mean*search_left_rate);
        stop=idx(n)+round(T0mean*search_rate);
    end

    if start < 1
        start=1;
    end

    % Check start and end points of detected intervals
    if stop > length(MBS) && start < length(MBS)
        stop=length(MBS);
    elseif stop > length(MBS) && start >= length(MBS)
        break
    end

    interval(n,1)=start;
    interval(n,2)=stop;
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function hogpic = jettify(hogpic)
%Turn a matrix which can be viewed via imagesc into an image using
%the same jet scheme, but now the result can be written to a file
NC = 200;
colorsheet = jet(NC);
dists = hogpic(:);
dists = dists - min(dists);
dists = dists / (max(dists)+eps);
dists = round(dists*(NC-1)+1);
colors = colorsheet(dists,:);
hogpic = reshape(colors,[size(hogpic,1) size(hogpic,2) 3]);

author: quantombone

repo: exemplarsvm


text:
function test_bug2032

% WALLTIME 00:20:00
% MEM 2gb
% DEPENDENCY

% read each of the atlasses and ensure that it specifies the units and coordsys

[ftver, ftpath] = ft_version;

%%

filename = {
  'template/atlas/aal/ROI_MNI_V4.nii'
  'template/atlas/afni/TTatlas+tlrc.BRIK'
  'template/atlas/brainweb/brainweb_discrete.mat'
  'template/atlas/brainweb/brainweb_fuzzy.mat'
  'template/atlas/spm_anatomy/AllAreas_v17.img'
  'template/atlas/spm_anatomy/AllAreas_v18.img'
  'template/atlas/vtpm/vtpm.mat'
  };

for i=1:numel(filename)
  atlas = ft_read_atlas(fullfile(ftpath, filename{i}));
  assert(isfield(atlas, 'unit'), sprintf('unit missing in %s', filename{i}));
  assert(isfield(atlas, 'coordsys'), sprintf('coordsys missing in %s', filename{i}));
end

%%

% the inflated ones do not have a coordsys
%  'template/anatomy/surface_inflated_both.mat'
%  'template/anatomy/surface_inflated_both_caret.mat'
%  'template/anatomy/surface_inflated_left.mat'
%  'template/anatomy/surface_inflated_left_caret.mat'
%  'template/anatomy/surface_inflated_right.mat'
%  'template/anatomy/surface_inflated_right_caret.mat'

filename = {
  'template/anatomy/surface_pial_both.mat'
  'template/anatomy/surface_pial_left.mat'
  'template/anatomy/surface_pial_right.mat'
  'template/anatomy/surface_white_both.mat'
  'template/anatomy/surface_white_left.mat'
  'template/anatomy/surface_white_right.mat'
  };

for i=1:numel(filename)
  surface = ft_read_headshape(fullfile(ftpath, filename{i}));
  assert(isfield(surface, 'unit'), sprintf('unit missing in %s', filename{i}));
  assert(isfield(surface, 'coordsys'), sprintf('coordsys missing in %s', filename{i}));
end


text:
function varargout = lzmadecode(varargin)
%
% output = lzmadecode(input)
%    or
% output = lzmadecode(input,info)
%
% Decompressing an LZMA-compressed byte-stream to recover the original data
% This function depends on the ZMat toolbox (http://github.com/fangq/zmat)
%
% authors:Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%      input: a string, int8/uint8 vector or numerical array to store LZMA-compressed data
%      info (optional): a struct produced by the zmat/lzmaencode function during
%            compression; if not given, the inputs/outputs will be treated as a
%            1-D vector
%
% output:
%      output: the decompressed byte stream stored in a uint8 vector; if info is
%            given, output will restore the original data's type and dimensions
%
% examples:
%      [bytes, info]=lzmaencode(eye(10));
%      orig=lzmadecode(bytes,info);
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

if(nargin==0)
    error('you must provide at least 1 input');
end
if(exist('zmat','file')==2 || exist('zmat','file')==3)
    if(nargin>1)
        [varargout{1:nargout}]=zmat(varargin{1},varargin{2:end});
    else
        [varargout{1:nargout}]=zmat(varargin{1},0,'lzma',varargin{2:end});
    end
else
    error('you must install ZMat toolbox to use this feature: http://github.com/fangq/zmat')
end


text:
function A = hmf(A,n)

%HMF Hybrid median filtering.
%   B = HMF(A,N) performs hybrid median filtering of the matrix A using a
%   NxN box. Hybrid median filtering preserves edges better than a NxN
%   square kernel-based median filter because data from different spatial
%   directions are ranked separately. Three median values are calculated in
%   the NxN box: MR is the median of horizontal and vertical R pixels, and
%   MD is the median of diagonal D pixels. The filtered value is the median
%   of the two median values and the central pixel C: median([MR,MD,C]).
%   For N = 5:
%        |D  *  R  *  D|
%        |*  D  R  D  *|
%        |R  R  C  R  R|
%        |*  D  R  D  *|
%        |D  *  R  *  D|
%
%   B = HMF(A) uses N = 5 (default value).
%
%   A can be a 2-D array or an RGB image. If A is an RGB image, hybrid
%   median filtering is performed in the HSV color space.
%
%   Notes
%   -----
%   1) N must be odd. If N is even then N is incremented by 1.
%   2) The Image Processing Toolbox is required.
%   3) If the function NANMEDIAN exists (Statistics Toolbox), NaN are
%      treated as missing values and are ignored.
%
%   Examples
%   --------
%     % -- original grayscale image --
%     I = imread('eight.tif');
%     % noisy image
%     J = imnoise(I,'salt & pepper',0.03);
%     % hybrid median filtering
%     K = hmf(J);
%     % figures
%     subplot(121),imshow(J),subplot(122),imshow(K)
%
%     % -- original RGB image --
%     [I,map] = imread('trees.tif');
%     I = ind2rgb(I,map);
%     % noisy image
%     J = imnoise(I,'salt & pepper',0.02);
%     % hybrid median filtering (using a 9x9 box)
%     K = hmf(J,9);
%     % figures
%     subplot(121),imshow(J),subplot(122),imshow(K)
%
%   See also MEDFILT2, MEDFILT2RGB, MEDFILT3, COLFILT
%
%   -- Damien Garcia -- 2007/08, revised 2010/02

error(nargchk(1,2,nargin));
if nargin==1, n = 5; end

if ~isscalar(n) || n<0, n = 5; end

% --- n must be odd. If not then n = n+1
n = round(n);
if rem(n+1,2)~=0, n = n+1; end

% --- Do we have an RGB image?
% RGB images can be only be uint8, uint16, single, or double
isRGB = ndims(A)==3 && (isfloat(A) || isa(A,'uint8') || isa(A,'uint16'));
% ---- Adapted from the obsolete function ISRGB ----
if isRGB
   if isfloat(A)
      % At first just test a small chunk to get a possible quick negative
      mm = size(A,1);
      nn = size(A,2);
      chunk = A(1:min(mm,10),1:min(nn,10),:);
      isRGB = (min(chunk(:))>=0 && max(chunk(:))<=1);
      % If the chunk is an RGB image, test the whole image
      if isRGB
         isRGB = (min(A(:))>=0 && max(A(:))<=1);
      end
   end
end
% ---- end of isrgb ----

assert(isRGB | ndims(A)==2,...
    'The input must be a 2-D array or an RGB image.')

classA = class(A);

% --- If the input is an RGB image, HMF is used in the HSV color space
if isRGB
    A = rgb2hsv(A);
    for k = 1:3, A(:,:,k) = hmf(A(:,:,k),n); end
    A = hsv2rgb(A);
    % HSV2RGB returns a double: change the class if necessary
    switch classA
        case 'uint8'
            A = uint8(A*255);
        case 'uint16'
            A = uint16(A*65535);
        case 'single'
            A = single(A);
    end
    return
end

% --- Plus & Cross masks
Plus = false(n,n);
Plus((n+1)/2,:) = true;
Plus(:,(n+1)/2) = true;
Plus = Plus(:);
Cross = false(n,n);
Cross((1:n)+n*(0:n-1)) = true;
Cross((1:n)+n*((n-1):-1:0)) = true;
Cross = Cross(:);

%% --- Hybrid median filtering

% Note: NANMEDIAN is used if this function exists (Statistics Toolbox)
existNaNmedian = exist('nanmedian','file');

% --- the COLFILT function zero-pads! => replicate boundaries
A = padarray(A,[(n-1)/2 (n-1)/2],'replicate');

M1 = colfilt(A,[n n],'sliding',@CrossMedian);
M2 = colfilt(A,[n n],'sliding',@PlusMedian);
if existNaNmedian
    A = nanmedian(cat(3,A,M1,M2),3);
else
    A = median(cat(3,A,M1,M2),3);
end

% remove the borders that were added by PADARRAY
A = A((n+1)/2:end-(n-1)/2,(n+1)/2:end-(n-1)/2);

function CM = CrossMedian(X)
    ncol = size(X,2);
    I = repmat(Cross,[1 ncol]);
    X = reshape(X(I),[2*n-1 ncol]);
    if existNaNmedian
        CM = nanmedian(X);
    else
        CM = median(X);
    end
end

function PM = PlusMedian(X)
    ncol = size(X,2);
    I = repmat(Plus,[1 ncol]);
    X = reshape(X(I),[2*n-1 ncol]);
    if existNaNmedian
        PM = nanmedian(X);
    else
        PM = median(X);
    end
end

end


text:
%[1983]-"Optimization by Simulated Annealing"

% (9/12/2020)

function SA = jSimulatedAnnealing(feat,label,opts)
% Parameters
c  = 0.93;  % cooling rate
T0 = 100;   % initial temperature

if isfield(opts,'T'), max_Iter = opts.T; end
if isfield(opts,'c'), c = opts.c; end
if isfield(opts,'T0'), T0 = opts.T0; end

% Objective function
fun = @jFitnessFunction;
% Number of dimensions
dim = size(feat,2);
% Initial
X   = jInitialization(1,dim);
% Fitness
fit = fun(feat,label,X,opts);
% Initial best
Xgb  = X;
fitG = fit;
% Pre
curve = zeros(1,max_Iter);
t = 2;
% Iterations
while t <= max_Iter
	% Probabilty of swap, insert, flip & eliminate
  prob = randi([1,4]);
  % Swap operation
  if prob == 1
    Xnew  = X;
    % Find index with bit '0' & '1'
    bit0  = find(X == 0);
    bit1  = find(X == 1);
    len_0 = length(bit0);
    len_1 = length(bit1);
    % Solve issue with missing bit '0' or '1'
    if len_0 ~= 0 && len_1 ~= 0
      % Get one random index from x1 & x2
      ind0 = randi([1,len_0]);
      ind1 = randi([1,len_1]);
      % Swap between two index
      Xnew(bit0(ind0)) = 1;
      Xnew(bit1(ind1)) = 0;
    end

  % Insert operation
  elseif prob == 2
    Xnew  = X;
    % Find index with zero
    bit0  = find(X == 0);
    len_0 = length(bit0);
    % Solve problem when all index are '1'
    if len_0 ~= 0
      ind = randi([1,len_0]);
      % Add one feature
      Xnew(bit0(ind)) = 1;
    end

  % Eliminate operation
  elseif prob == 3
    Xnew  = X;
    % Find index with one
    bit1  = find(X == 1);
    len_1 = length(bit1);
    % Solve problem when all index are '0'
    if len_1 ~= 0
      ind = randi([1,len_1]);
      % Remove one feature
      Xnew(bit1(ind)) = 0;
    end

  % Flip operation
  elseif prob == 4
    Xnew = X;
    % Flip all variables
    for d = 1:dim
      Xnew(d) = 1 - Xnew(d);
    end
  end

  % Fitness
  Fnew = fun(feat,label,Xnew,opts);
  % Global best update
  if Fnew <= fitG
    Xgb  = Xnew;
    fitG = Fnew;
    X    = Xnew;
  % Accept worst solution with probability
  else
    % Delta energy
    delta = Fnew - fitG;
    % Boltzmann Probility
    P     = exp(-delta / T0);
    if rand() <= P
      X = Xnew;
    end
  end
  % Temperature update
  T0 = c * T0;
  % Save
  curve(t) = fitG;
  fprintf('\nIteration %d Best (SA)= %f',t,curve(t))
  t = t + 1;
end
% Select features
Pos   = 1:dim;
Sf    = Pos(Xgb == 1);
sFeat = feat(:,Sf);
% Store results
SA.sf = Sf;
SA.ff = sFeat;
SA.nf = length(Sf);
SA.c  = curve;
SA.f  = feat;
SA.l  = label;
end

function X = jInitialization(N,dim)
% Initialize X vectors
X = zeros(N,dim);
for i = 1:N
  for d = 1:dim
    if rand() > 0.5
      X(i,d) = 1;
    end
  end
end
end

author: JingweiToo

repo: Wrapper-Feature-Selection-Toolbox


text:
function bin_seq = getBits(text)
matrix  = dec2bin(uint8(text),8);
bin_seq = reshape(matrix', 1, 8*length(text));
end

author: ktekeli

repo: audio-steganography-algorithms


text:
function [m0, l0, a0, b0] = tapas_rdcm_get_prior(DCM)
% [m0, l0, a0, b0] = tapas_rdcm_get_prior(DCM)
%
% Returns prior parameters on model parameters (theta) and noise precision
% (tau) for the particular connectivity pattern
%
%   Input:
%   	DCM             - model structure
%
%   Output:
%       m0              - prior mean
%       l0              - prior covariance
%       a0              - prior shape parameter
%       b0              - prior rate parameter
%

% ----------------------------------------------------------------------
%
% Authors: Stefan Fraessle (stefanf@biomed.ee.ethz.ch), Ekaterina I. Lomakina
%
% Copyright (C) 2016-2022 Translational Neuromodeling Unit
%                         Institute for Biomedical Engineering
%                         University of Zurich & ETH Zurich
%
% This file is part of the TAPAS rDCM Toolbox, which is released under the
% terms of the GNU General Public License (GPL), version 3.0 or later. You
% can redistribute and/or modify the code under the terms of the GPL. For
% further see COPYING or <http://www.gnu.org/licenses/>.
%
% Please note that this toolbox is in an early stage of development. Changes
% are likely to occur in future releases.
%
% ----------------------------------------------------------------------

% get the number of regions and inputs
[nr, nu] = size(DCM.c);

% get the priors on model parameters
[pE,pC] = tapas_rdcm_spm_dcm_fmri_priors(DCM.a,DCM.b,DCM.c,DCM.d);

% set the prior mean of endogenous parameters to zero
pE.A = zeros(size(pE.A))+diag(diag(pE.A));

% specify the prior mean
m0 = [pE.A reshape(pE.B,nr,nr*nu) pE.C];

% prior precision
pC.A       = 1./pC.A;
pC.B       = 1./pC.B;
pC.C       = 1./pC.C;
pC.D       = 1./pC.D;
pC.transit = 1./pC.transit;
pC.decay   = 1./pC.decay;
pC.epsilon = 1./pC.epsilon;

% prior precision of baseline
if ( any(pC.C(:,end)) )
    pC.C(:,end) = 10^-8;
end

% prior precision
l0 = [pC.A reshape(pC.B,nr,nr*nu) pC.C];

% Setting priors on noise precision
a0 = 2;
b0 = 1;

end

author: translationalneuromodeling

repo: tapas


text:
function mask=getMask_onlineEvaluation(fn_sb)
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% @InProceedings{ZhengICCV09,
%   author = {Yuanjie Zheng and Chandra Kambhamettu},
%   title = {Learning Based Digital Matting},
%   booktitle = {The 20th IEEE International Conference on Computer Vision},
%   year = {2009},
%   month = {September--October}
% }
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% Copyright Yuanjie Zheng @ PICSL @ UPenn on May 9th, 2011
% zheng.vision@gmail.com
% http://sites.google.com/site/zhengvision/
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
I=imread(fn_sb);
mask=zeros(size(I,1),size(I,2));
fore=(I==255);
back=(I==0);

mask(fore)=1;
mask(back)=-1;


text:
function [err, errK] = getHcurlerror3ND(node,elem,curlE,Eh,markedElem)
%% GETHCURLERROR3ND Hcurl norm of approximation error for the lowest order Nedelect element in 3-D.
%
% err = GETHCURLERROR3ND(node,elem,curlE,Eh,markedElem);
%
% Example
%
%     node = [-1,-1,-1; 1,-1,-1; 1,1,-1; -1,1,-1; -1,-1,1; 1,-1,1; 1,1,1; -1,1,1];
%     elem = [1,2,3,7; 1,6,2,7; 1,5,6,7; 1,8,5,7; 1,4,8,7; 1,3,4,7];
%     maxIt = 4;
%     HcurlErr = zeros(maxIt,1);
%     N = zeros(maxIt,1);
%     for k = 1:maxIt
%         [node,elem] = uniformbisect3(node,elem);
%         [elem2dof,edge] = dof3edge(elem);
%         pde = Maxwelldata2;
%         uI = edgeinterpolate(pde.exactu,node,edge);
%         HcurlErr(k) = getHcurlerror3ND(node,elem,pde.curlu,uI);
%         N(k) = length(uI);
%     end
%     r = showrate(N,HcurlErr,1,'b-+');
%     legend('||u-u_I||_{curl}',['N^{' num2str(r) '}'],'LOCATION','Best');
%
% See also getHcurlerror3ND1, getHcurlerror3ND2, getL2error3ND
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Sort elem to ascend ordering
elem = sort(elem,2);

%% Construct Data Structure
elem2dof = dof3edge(elem);
NT = size(elem,1);% Ndof = max(elem2dof(:)); %N = size(node,1);
% compute curl phi
[Dlambda,volume] = gradbasis3(node,elem);
% curl phi = 2*Dlambda_i cross Dlambda_j;
curlPhi(:,:,6) = 2*cross(Dlambda(:,:,3),Dlambda(:,:,4),2);
curlPhi(:,:,1) = 2*cross(Dlambda(:,:,1),Dlambda(:,:,2),2);
curlPhi(:,:,2) = 2*cross(Dlambda(:,:,1),Dlambda(:,:,3),2);
curlPhi(:,:,3) = 2*cross(Dlambda(:,:,1),Dlambda(:,:,4),2);
curlPhi(:,:,4) = 2*cross(Dlambda(:,:,2),Dlambda(:,:,3),2);
curlPhi(:,:,5) = 2*cross(Dlambda(:,:,2),Dlambda(:,:,4),2);
curlEhp = zeros(NT,3);
for k = 1:6
    curlEhp = curlEhp + ...
    repmat(Eh(elem2dof(:,k)),1,3).*curlPhi(:,:,k);
end

%% compute Hcurl error element-wise
[lambda,w] = quadpts3(2);
nQuad = size(lambda,1);
errK = zeros(NT,1);
for p = 1:nQuad
    % quadrature points in the x-y-z coordinate
    pxy = lambda(p,1)*node(elem(:,1),:) ...
        + lambda(p,2)*node(elem(:,2),:) ...
        + lambda(p,3)*node(elem(:,3),:) ...
        + lambda(p,4)*node(elem(:,4),:);
    if isnumeric(curlE) % a constant vector
        curlEp = repmat(curlE,NT,1);
    else % function handel
        curlEp = curlE(pxy);
    end
%     curlEp = curlE(pxy);
    % compute Ehp at quadrature points
    errK = errK + w(p)*sum((curlEp - curlEhp).^2,2);
end
errK = errK.*volume;
% modify the error
errK(isnan(errK)) = 0; % remove the singular part
if (nargin == 5) && ~isempty(markedElem)
    errK = errK(markedElem); % error on some marked region
end
err = sqrt(sum(errK));
%% TODO write more M-lint

repo: ifem


text:
function events = in_events_graph(sFile, EventFile)
% IN_EVENTS_GRAPH: Read Neuromag Graph events files
%
% This file reads different variations of the .evl file in which
% Megin (Elekta) Graph users save their events.
%
% The following example shows the default evl export format from Graph.
% Example 1 of evl file:
%   | (beamformer::saved-event-list
%   |  :source-file "/path/file.fif"
%   |  :events '(
%   |   ((:time  28.8425) (:class :manual) (:length  0.069))
%   |   ((:time  31.194) (:class :manual) (:length  0.0415))
%   | ))
%
% The following example shows a variation of the default evl format which
% is used in some Megin (Elekta) labs.
% Example 2 of evl file:
%   | (beamformer::saved-event-list
%   |  :source-file "/path/file.fif"
%   |  :events '(
%   | ((:time  991.73) (:class "A") (:level  2.47042716145e-11))
%   | ((:time  1019.15) (:class "B") (:level  2.6210856402e-11))
%   | ))

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2015
%          Juan GPC, 2020

% Open file
fid = fopen(EventFile, 'r');
if (fid < 0)
    error(['Cannot open file ', EventFile]);
end

%list of reader functions for different formats
listOfReaders={@readDefaultFormat,...
          @readAlternativeFormat};
errorMsg = 'Unable to read the Graph events-list file.';

for i = 1:length(listOfReaders)
    try
        frewind(fid);
        events = listOfReaders{i}(fid);
    catch
        errorMsg = [errorMsg, char(10), lasterr];
        if (i < length(listOfReaders))
          continue;
        else
          error(errorMsg);
        end
    end
end

% Close file
fclose(fid);

function events = readDefaultFormat(fid)
% Initialize returned structure
events = repmat(db_template('event'), 0);
evtMat = zeros(2,0);
evtLabel = {};

% Read file line by line
while 1
    % Read one line
    newLine = fgetl(fid);
    if ~ischar(newLine)
        break;
    end
    % Strip spaces
    newLine(newLine == ' ') = [];
    % Lines to skip
    if isempty(newLine)
        continue;
    end
    % If the line does not contain ":time": skip
    if isempty(strfind(newLine, ':time'))
        continue;
    end
    % Parse the line: "((:timeXXXXXXX)(:class:manual)(:lengthYYYYYYY))"
    res = sscanf(newLine, '((:time%f)(:class:manual)(:length%f))');
    if (length(res) ~= 2)
        continue;
    end
    evtMat(:,end+1) = res(:);
    evtLabel{end+1} = 'Graph';
end

% If nothing was read: return
if isempty(evtMat)
    error('This file does not contain any Neuromag Graph events.');
end

% Find all the event types
uniqueLabel = unique(evtLabel);
% Convert to a structure matrix
for iEvt = 1:length(uniqueLabel)
    iOcc = find(strcmpi(uniqueLabel{iEvt}, evtLabel));
    events(iEvt).label       = uniqueLabel{iEvt};
    events(iEvt).epochs      = ones(1,length(iOcc));
    events(iEvt).times       = round(([evtMat(1,iOcc); evtMat(1,iOcc) + evtMat(2,iOcc)]) .* sFile.prop.sfreq) ./ sFile.prop.sfreq;
    events(iEvt).reactTimes  = [];
    events(iEvt).select      = 1;
    events(iEvt).channels   = [];
    events(iEvt).notes      = [];
end

end

function events = readAlternativeFormat(fid)

% Initialize returned structure
events = repmat(db_template('event'), 0);
evtMat = zeros(1,0);
evtLabel = {};

% Read file line by line
while 1
    % Read one line
    newLine = fgetl(fid);
    if ~ischar(newLine)
        break;
    end
    % Strip spaces
    newLine(newLine == ' ') = [];
    % Lines to skip
    if isempty(newLine)
        continue;
    end
    % If the line does not contain ":time": skip
    if ( isempty(strfind(newLine, ':time')) || ...
         isempty(strfind(newLine, '"')) ) %#ok<*STREMP>
        continue;
    end
    % Parse the line: "((:timeXXXXXX) (:class "XXX") (:levelXXXXXXXXX))"
    %
    res = sscanf(newLine, '((:time%f)');
    if isempty(res)
        continue;
    end

    posQuotes = strfind(newLine,'"');

    evtMat(:,end+1) = res;
    evtLabel{end+1} = newLine(posQuotes(1)+1:posQuotes(2)-1);
end

% If nothing was read: return
if isempty(evtMat)
    error('This file does not contain any Neuromag Graph events.');
end

% Find all the event types
uniqueLabel = unique(evtLabel);
% Convert to a structure matrix
for iEvt = 1:length(uniqueLabel)
    iOcc = find(strcmpi(uniqueLabel{iEvt}, evtLabel));
    events(iEvt).label       = uniqueLabel{iEvt};
    events(iEvt).epochs      = ones(1,length(iOcc));
    events(iEvt).times       = round(evtMat(1,iOcc)* sFile.prop.sfreq) ./ sFile.prop.sfreq;
    events(iEvt).reactTimes  = [];
    events(iEvt).select      = 1;
    events(iEvt).channels   = [];
    events(iEvt).notes      = [];
end

end

end

author: brainstorm-tools


text:
function d = dot(sF1, sF2)
% inner product
%
% Syntax
%   sF = dot(sF1, sF1)
%
% Input
%  sF1, sF2   - @S2FunHarmonic
%
% Output
%  sF - @S2FunHarmonic
%

bw = min(sF1.bandwidth,sF2.bandwidth);

sF1.bandwidth = bw;
sF2.bandwidth = bw;

d = sum(sF1.fhat .* conj(sF2.fhat));

end

author: mtex-toolbox

repo: mtex

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