text:
% dipoledensity() - compute and optionally plot a measure of the 3-D spatial
%                   (in)homogeneity of a specified (large) set of 1- or 2-dipole
%                   component models, either as physical dipole density or as
%                   dipole-position entropy across subjects. In either case,
%                   take into account either all the dipoles, or only the nearest
%                   dipole from each of the subjects. If no output arguments,
%                   or if 'plot','on', paints a 3-D density|entropy brain image
%                   on slices of the Montreal Neurological Institute (MNI) mean
%                   MR brain image ('standard_BESA/avg152t1.mat'). Calls
%                   dipplot(),
%                   mri3dplot(), and Fieldtrip function find_inside_vol().
% Usage:
%               >> [dens3d mri] = dipoledensity( dipoles, 'key',val, ... );
%
% Inputs:
%    dipoles - this may be either the same dipole structure given as input to
%              the dipplot() function, a 3 by n array of dipole localization or
%              a cell array containing arguments for the dipplot function. Note that
%              the 'coordformat' option below defines the coordinate space for these
%              dipoles (default is MNI). See help dipplot for more information.
%
% Optional 'key', val input pairs:
%  'mri'        - [string or struct] mri file (matlab format or file format read
%                 by fcdc_read_mri). See dipplot.m help for more information.
%  'method'     - ['alldistance'|'distance'|'entropy'|'relentropy'] method for
%                            computing density:
%                 'alldistance' - {default} take into account the gaussian-weighted
%                            distances from each voxel to all the dipoles. See
%                            'methodparam' (below) to specify a standard deviation
%                            (in mm) for the gaussian weight kernel.
%                 'distance' - take into account only the distances to the nearest
%                              dipole for each subject. See 'methodparam' (below).
%                 'entropy' - taking into account only the nearest dipole to each
%                             voxel for each subject. See 'methodparam' below.
%                 'relentropy' - as in 'entropy,' but take into account all the
%                             dipoles for each subject.
% 'methodparam' - [number] for 'distance'|'alldistance' methods (see above), the
%                 standard deviation (in mm) of the 3-D gaussian smoothing kernel.
%                 For 'entropy'|'relentropy' methods, the number of closest dipoles
%                 to include {defaults: 20 mm | 20 dipoles }
% 'subsample'   - [integer] subsampling of native MNI image {default: 2 -> 2x2x2}
% 'weight'      - [(1,ncomps) array] for 'distance'|'alldistance' methods, the
%                 relative weight of each component dipole {default: ones()}
% 'coordformat' - ['mni'|'spherical'] coordinate format if dipole location or
%                 a structure is given as input. Default is 'mni'.
% 'subjind'     - [(1,ncomps) array] subject index for each dipole model. If two
%                 dipoles are in one component model, give only one subject index.
% 'nsessions'   - [integer] for 'alldistance' method, the number of sessions to
%                 divide the output values by, so that the returned measure is
%                 dipole density per session {default: 1}
% 'plot'        - ['on'|'off'] force plotting dipole density|entropy
%                 {default: 'on' if no output arguments, else 'off'}
% 'dipplot'     - ['on'|'off'] plot the dipplot image (used for converting
%                 coordinates (default is 'off')
% 'plotargs'    - {cell array} plotting arguments for mri3dplot() function.
% 'volmesh_fname' - [string] precomputed mesh volume file name. If not
%                 given as input the function will recompute it (it can take from
%                 five to 20 minutes). By default this function save the volume file
%                 mesh into a file named volmesh_local.mat in the current
%                 folder.
%
% Outputs:
%  dens3d       - [3-D num array] density in dipoles per cubic centimeter. If output
%                 is returned, no plot is produced unless 'plot','on' is specified.
%  mri          - {MRI structure} used in mri3dplot().
%
% Example:
%         >> fakedipoles = (rand(3,10)-0.5)*80;
%         >> [dens3d mri] = dipoledensity( fakedipoles, 'coordformat', 'mni');
%         >> mri3dplot(dens3d,mri); % replot if no output is given above
%                                   % function is called automatically
%
% ------------------------------------
% NOTES:   to do multiple subject causal-weighted density map,
% (1) concatenate dipplot coord matrices for all subject
% (2) make g.subjind vector [ones(1,ncompsS1) 2*ones(1,ncompsS2) ... N*ones(1,ncompssN)]
% (3) concatenate normalized outflows for all subjects to form weight vector
% (4) call dipoledensity function with method = 'entropy' or 'relentropy'
% ------------------------------------
%
% See also:
%           EEGLAB: dipplot(), mri3dplot(), Fieldtrip: find_inside_vol()
%
% Author: Arnaud Delorme & Scott Makeig, SCCN, INC, UCSD

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

function [prob3d, mri] = dipoledensity(dipplotargs, varargin)

    % TO DO: return in dipplot() the real 3-D location of dipoles (in posxyz)
    %        FIX the dimension order here

prob3d = []; mri = [];
if nargin < 1
   help dipoledensity
   return
end

g = finputcheck(varargin, { 'subjind'     'integer'  []               [];
                            'method' 'string' { 'relentropy' 'entropy' 'distance' 'alldistance' } 'alldistance';
                            'methodparam' 'real'     []               20;
                            'weight'      { 'real' 'cell' }  []               [];
                            'smooth'      'real'     []               0;
                            'nsessions'   'integer'  []               1;
                            'subsample'   'integer'  []               2;
                            'plotargs'    'cell'     []               {};
                            'plot'        'string'  { 'on' 'off' }    fastif(nargout == 0, 'on', 'off');
                            'dipplot'     'string'  { 'on' 'off' }   'off';
                            'coordformat' 'string'  { 'mni' 'spherical' }   'mni';
                            'normalization' 'string'  { 'on' 'off' } 'on';
                            'volmesh_fname' 'string'  []  'volmesh_local.mat';
                            'mri'         { 'struct' 'string' } [] ''});
if isstr(g), error(g); end;
if ~strcmpi(g.method, 'alldistance') & isempty(g.subjind)
    error('Subject indices are required for this method');
end;
if ~iscell(g.weight), g.weight = { g.weight }; end;

% plotting dipplot
% ----------------
if ~iscell(dipplotargs) % convert input
    if ~isstruct(dipplotargs)
        if size(dipplotargs,1) == 3, dipplotargs = dipplotargs';
        elseif size(dipplotargs,2) ~= 3
            error('If an array of dipoles is given as entry, there must be 3 columns or 3 rows for x y z');
        end;
        model = [];
        for idip = 1:length(dipplotargs)
            model(idip).posxyz = dipplotargs(idip,:);
            model(idip).momxyz = [1 0 0];
            model(idip).rv = 0.5;
        end;
        dipplotargs = model;
    end;
    dipplotargs = { dipplotargs 'coordformat' g.coordformat };
else
    dipplotargs = { dipplotargs{:} 'coordformat' g.coordformat };
end;
struct = dipplot(dipplotargs{:}, 'plot', g.dipplot);
if nargout == 0
    drawnow;
end;

% retrieve coordinates in MNI space
% ---------------------------------
if 0 % deprecated
     % find dipoles
     % ------------
    hmesh = findobj(gcf, 'tag', 'mesh');
    if isempty(hmesh), error('Current figure must contain dipoles'); end;
    hh = [];
    disp('Finding dipoles...');
    dips = zeros(1,200);
    for index = 1:1000
        hh = [ hh(:); findobj(gcf, 'tag', ['dipole' int2str(index) ]) ];
        dips(index) = length(findobj(gcf, 'tag', ['dipole' int2str(index) ]));
    end;

    disp('Retrieving dipole positions ...');
    count = 1;
    for index = 1:length(hh)
        tmp = get(hh(index), 'userdata');
        if length(tmp) == 1
            allx(count) = tmp.eleccoord(1,1);
            ally(count) = tmp.eleccoord(1,2);
            allz(count) = tmp.eleccoord(1,3);
            alli(count) = index;
            count = count + 1;
        end;
    end;
end;

% check weights
% -------------
if ~isempty(g.weight{1})
    if ~iscell(g.weight)
        if length(g.weight) ~= length(struct)
            error('There must be as many elements in the weight matrix as there are dipoles')
        end;
    else
        if length(g.weight{1}) ~= length(struct) || length(g.weight{1}) ~= length(g.weight{end})
            error('There must be as many elements in the weight matrix as there are dipoles')
        end;
    end;
else
    g.weight = { ones( 1, length(struct)) };
end;
if ~isempty(g.subjind)
    if length(g.subjind) ~= length(struct)
        error('There must be as many element in the subject matrix as there are dipoles')
    end;
else
    g.subjind = ones( 1, length(struct));
end;

% decoding dipole locations
% -------------------------
disp('Retrieving dipole positions ...');
count = 1;
for index = 1:length(struct)
    dips = size(struct(index).eleccoord,1);
    for dip = 1:dips
        allx(count) = struct(index).eleccoord(dip,1);
        ally(count) = struct(index).eleccoord(dip,2);
        allz(count) = struct(index).eleccoord(dip,3);
        alli(count) = index;
        allw1(count) = g.weight{1}(  index)/dips;
        allw2(count) = g.weight{end}(index)/dips;
        alls(count) = g.subjind(index);
        count = count + 1;
    end;
end;
g.weight{1}    = allw1;
g.weight{end}  = allw2;
g.subjind = alls;

% read MRI file
% -------------
if isempty(g.mri) % default MRI file
    dipfitdefs;
    load('-mat', template_models(1).mrifile); % load mri variable
    g.mri = mri;
end
if isstr(g.mri)
    try,
        mri = load('-mat', g.mri);
        mri = mri.mri;
    catch,
        disp('Failed to read Matlab file. Attempt to read MRI file using function read_fcdc_mri');
        try,
            warning off;
            mri = read_fcdc_mri(g.mri);
            mri.anatomy = round(gammacorrection( mri.anatomy, 0.8));
            mri.anatomy = uint8(round(mri.anatomy/max(reshape(mri.anatomy, prod(mri.dim),1))*255));
            % WARNING: if using double instead of int8, the scaling is different
            % [-128 to 128 and 0 is not good]
            % WARNING: the transform matrix is not 1, 1, 1 on the diagonal, some slices may be
            % misplaced
            warning on;
        catch,
            error('Cannot load file using read_fcdc_mri');
        end;
    end;
    g.mri = mri; % output the anatomic mri image
end;

% reserve array for density
% -------------------------
prob3d = {zeros(ceil(g.mri.dim/g.subsample)) };
for i = 2:length(g.weight), prob3d{i} = prob3d{1}; end;

% compute voxel size
% ------------------
point1 = g.mri.transform * [ 1 1 1 1 ]';
point2 = g.mri.transform * [ 2 2 2 1 ]';
voxvol = sum((point1(1:3)-point2(1:3)).^2)*g.subsample^3; % in mm

% compute global subject entropy if necessary
% -------------------------------------------
vals   = unique(g.subjind); % the unique subject indices
if strcmpi(g.method, 'relentropy') | strcmpi(g.method, 'entropy') %%%%% entropy %%%%%%%
    newind = zeros(size(g.subjind));
    for index = 1:length(vals) % foreach subject in the cluster
        tmpind = find(g.subjind == vals(index)); % dipoles for the subject
        totcount(index) = length(tmpind); % store the number of subject dipoles
        newind(tmpind) = index; % put subject index into newind
    end;
    g.subjind = newind;
    gp = totcount/sum(totcount);
    globent = -sum(gp.*log(gp));
end;

% compute volume inside head mesh
% -------------------------------
dipfitdefs; % get the location of standard BEM volume file
tmp = load('-mat',DIPOLEDENSITY_STDBEM); % load MNI mesh

if isempty(g.volmesh_fname)  % default
    filename = [ '/home/arno/matlab/MNI_VoxelTsearch' int2str(g.subsample) '.mat' ];
else
    filename = g.volmesh_fname; %
end
if ~exist(filename)
    disp('Computing volume within head mesh...');
    [X Y Z]           = meshgrid(g.mri.xgrid(1:g.subsample:end)+g.subsample/2, ...
                                 g.mri.ygrid(1:g.subsample:end)+g.subsample/2, ...
                                 g.mri.zgrid(1:g.subsample:end)+g.subsample/2);
    [indX indY indZ ] = meshgrid(1:length(g.mri.xgrid(1:g.subsample:end)), ...
                                 1:length(g.mri.ygrid(1:g.subsample:end)), ...
                                 1:length(g.mri.zgrid(1:g.subsample:end)));
    allpoints = [ X(:)'    ; Y(:)'   ; Z(:)' ];
    allinds   = [ indX(:)' ; indY(:)'; indZ(:)' ];
    allpoints = g.mri.transform * [ allpoints ; ones(1, size(allpoints,2)) ];
    allpoints(4,:) = [];

    olddir = pwd;
    tmppath = which('ft_electroderealign');
    tmppath = fullfile(fileparts(tmppath), 'private');
    cd(tmppath);
    [Inside Outside] = find_inside_vol(allpoints', tmp.vol); % from Fieldtrip
    cd(olddir);
    disp('Done.');

    if 0 % old code using Delaunay %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        P = tmp.vol.bnd(1).pnt;
        T = delaunayn(P); % recompute triangularization (the original one is not compatible
                          % with tsearchn) get coordinates of all points in the volume
        % search for points inside or ouside the volume (takes about 14 minutes!)
        IO = tsearchn(P, T, allpoints');
        Inside        = find(isnan(IO));
        Outside       = find(~isnan(IO));
        disp('Done.');
    end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    try,
        save('-mat', filename, 'allpoints', 'allinds', 'Inside', 'Outside');
        disp('Saving file containing inside/outide voxel indices...');
    catch, end;
else
    disp('Loading file containing inside/outide voxel indices...');
    load('-mat',filename);
end;
InsidePoints  = allpoints(:, Inside);
InsideIndices = allinds(:, Inside);

% scan grid and compute entropy at each voxel
% -------------------------------------------
edges = [0.5:1:length(vals)+0.5];

if ~strcmpi(g.method, 'alldistance')
    fprintf('Computing (of %d):', size(InsideIndices,2));
    % entropy calculation: have to scan voxels
    % ----------------------------------------
    for i = 1:size(InsideIndices,2)

        alldists = (InsidePoints(1,i) - allx).^2 ...
                 + (InsidePoints(2,i) - ally).^2 ...
                 + (InsidePoints(3,i) - allz).^2;
        [tmpsort indsort] = sort(alldists); % sort dipoles by distance
        tmpweights{1}   = g.weight{1}(  indsort);
        tmpweights{end} = g.weight{end}(indsort);

        if strcmpi(g.method, 'relentropy') | strcmpi(g.method, 'entropy') %%%%% entropy %%%%%%%

            subjs  = g.subjind(indsort(1:g.methodparam)); % get subject indices of closest dipoles
            p      = histc(subjs, edges);
            if strcmpi(g.method, 'relentropy')
                p      = p(1:end-1)./totcount;
                % this should be uniform if p conforms to global count for all subjects
            end;
            p      = p/sum(p);
            p(find(p == 0)) = [];
            for tmpi = 1:length(g.weight)
                prob3d{1}(InsideIndices(1,i), InsideIndices(2,i), InsideIndices(3,i)) = -sum(p.*log(p));
            end;
        else
            % distance to each subject
            ordsubjs  = g.subjind(indsort);
            for index = 1:length(vals) % for each subject
                tmpind = find(ordsubjs == vals(index));
                if strcmpi(g.method,'distance')
                    use_dipoles(index) = tmpind(1); % find their nearest dipole
                end
            end;
            for tmpi = 1:length(g.weight)
                prob3d{tmpi}(InsideIndices(1,i), InsideIndices(2,i), InsideIndices(3,i)) = ...
                   sum(tmpweights{tmpi}(use_dipoles).*exp(-tmpsort(use_dipoles)/ ...
                           (2*g.methodparam^2))); % 3-D gaussian smooth
            end;
        end;
        if mod(i,100) == 0, fprintf('%d ', i); end;
    end;
else % 'alldistance'
    % distance calculation: can scan dipoles instead of voxels (since linear)
    % --------------------------------------------------------
    %alldists = allx.^2 + ally.^2 + allz.^2;
    %figure; hist(alldists); return; % look at distribution of distances

    fprintf('Computing (of %d):', size(allx,2));
    for tmpi=1:length(g.weight)
        tmpprob{tmpi} = zeros(1, size(InsidePoints,2));
    end;
    if length(g.weight) > 1, tmpprob2 = tmpprob; end;
    for i = 1:size(allx,2)
        alldists = (InsidePoints(1,:) - allx(i)).^2 + ...
                   (InsidePoints(2,:) - ally(i)).^2 + ...
                   (InsidePoints(3,:) - allz(i)).^2;
%         alldists = 1;  % TM
        for tmpi=1:length(g.weight)
            tmpprob{tmpi} = tmpprob{tmpi} + g.weight{tmpi}(i)*exp(-alldists/(2*g.methodparam^2)); % 3-D gaussian smooth
            if any(isinf(tmpprob{tmpi})), error('Infinite value in probability calculation'); end;
        end;
        if mod(i,50) == 0, fprintf('%d ', i); end;
    end;
    % copy values to 3-D mesh
    % -----------------------
    for i = 1:length(Inside)
        pnts = allinds(:,Inside(i));
        for tmpi = 1:length(g.weight)
            prob3d{tmpi}(pnts(1), pnts(2), pnts(3)) = tmpprob{tmpi}(i);
        end;
    end;

end;
fprintf('\n');

% normalize for points inside and outside the volume
% --------------------------------------------------
if strcmpi(g.method, 'alldistance') && strcmpi(g.normalization,'on')
    for i =1:length(g.weight)
        disp('Normalizing to dipole/mm^3');
        if any(prob3d{i}(:)<0)
            fprintf('WARNING: Some probabilities are negative, this will likely cause problems when normalizing probabilities.\n');
            fprintf('It is highly recommended to turn normaliziation off by using ''normalization'' key to ''off''.\n');
        end;
        totval = sum(prob3d{i}(:));  % total values in the head
        totdip = size(allx,2);   % number of dipoles
        voxvol;                  % volume o af a voxel in mm^3
        prob3d{i} = prob3d{i}/totval*totdip/voxvol*1000; % time 1000 to get cubic centimeters
        prob3d{i} = prob3d{i}/g.nsessions;
    end;
end;

% resample matrix
% ----------------
if g.subsample ~= 1
    for i =1:length(g.weight)
        prob3d{i} = prob3d{i}/g.subsample;
        newprob3d = zeros(g.mri.dim);
        X = ceil(g.mri.xgrid/g.subsample);
        Y = ceil(g.mri.ygrid/g.subsample);
        Z = ceil(g.mri.zgrid/g.subsample);
        for index = 1:size(newprob3d,3)
            newprob3d(:,:,index) = prob3d{i}(X,Y,Z(index));
        end;
        prob3d{i} = newprob3d;
    end;
end;

% 3-D smoothing
% -------------
if g.smooth ~= 0
    disp('Smoothing...');
    for i =1:length(g.weight)
        prob3d{i} = smooth3d(prob3d{i}, g.smooth);
    end;
end;

% plotting
% --------
if strcmpi(g.plot, 'off')
    close gcf;
else
    mri3dplot( prob3d, g.mri, g.plotargs{:}); % plot the density using mri3dplot()
end;
return;

author: PatternRecognition

repo: OpenBMI


text:
function value = trinomial ( i, j, k )

%*****************************************************************************80
%
%% TRINOMIAL computes a trinomial coefficient.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/triangle_integrals/trinomial.m
%
%  Discussion:
%
%    The trinomial coefficient is a generalization of the binomial
%    coefficient.  It may be interpreted as the number of combinations of
%    N objects, where I objects are of type 1, J of type 2, and K of type 3.
%    and N = I + J + K.
%
%    T(I,J,K) = N! / ( I! J! K! )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, J, K, the factors.
%    All should be nonnegative.
%
%    Output, integer VALLUE, the trinomial coefficient.
%

%
%  Each factor must be nonnegative.
%
  if ( i < 0 || j < 0 || k < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRINOMIAL - Fatal error%\n' );
    fprintf ( 1, '  Negative factor encountered.\n' );
    error ( 'TRINOMIAL - Fatal error%' );
  end

  value = 1;

  t = 1;

  for l = 1 : i
%   value = value * t / l;
    t = t + 1;
  end

  for l = 1 : j
    value = value * t / l;
    t = t + 1;
  end

  for l = 1 : k
    value = value * t / l;
    t = t + 1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [V]=ellipseCoord(A,t)

% function [V]=ellipseCoord(A,t)
% ------------------------------------------------------------------------
% Calculates ellipse coordiantes for the angles in t based on the vector A
% which defines the centre coordinates, the radii and the angle
% respectively.
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
% 2013/24/09
%------------------------------------------------------------------------

%%
x0=A(1);
y0=A(2);
x=A(3).*cos(t);
y=A(4).*sin(t);
V=[x(:) y(:) zeros(size(x(:)))];
[R,~]=euler2DCM([0 0 -A(5)]);
V=(R*V')';
V(:,1)=V(:,1)+x0;
V(:,2)=V(:,2)+y0;

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function magcumu() % autogenerated function wrapper
    %   This subroutine selects earthquakes in a magnitude, time
    %   and depth range for plotting cumulative curves
    %   Operates on the last subset of the catalogue (ZG.newcat).
    %   Changes ZG.newt2
    %
    %minma2 =  input('Please input Minimum Magnitude (inclusive):')
    %maxma2 =  input('Please input Maximum Magnitude:')
    %minde =  input('Please input minimum depth (inclusive):')
    %maxde =  input('Please input Maximum depth:')
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun(mfilename('fullpath'));
    ZG=ZmapGlobal.Data;
    % make selection from  catalogue ZG.newcat
    % ZG.newt2 is changed

    ZG.newt2 = ZG.newcat;

    l = ZG.newt2.Magnitude >= minma2 &...
        ZG.newt2.Magnitude <= maxma2 & ...
        ZG.newt2.Date >= mint & ...
        ZG.newt2.Date <= maxt;
    ZG.newt2 = ZG.newt2.subset(l);

    l = ZG.newt2.Depth >= minde &...
        ZG.newt2.Depth <= maxde ;
    ZG.newt2 = ZG.newt2.subset(l);

    %l = ZG.newt2.Date >= minti & ZG.newt2.Date <= maxti ;
    %ZG.newt2 = ZG.newt2.subset(l);

    stri = ['# ' stri1 '#  ' num2str(minma2) ' <= M <= ' num2str(maxma2) ...
        '#  ' num2str(minde) ' <= h(km) < ' num2str(maxde) ];

    ZG.t0b = min(ZG.newt2.Date);
    timeplot(ZG.newt2)
end

author: CelsoReyes


text:
% Chapter 3 - Complex Iterative Maps.
% Program 3b - The Mandelbrot Set.
% Thanks to Steve Lord from The MathWorks for his help.
% Copyright Birkhauser 2013. Stephen Lynch.

% Vectorized program.
% Plot the Mandelbrot set in black and white (Figure 3.2).
Nmax = 50; scale = 0.005;
xmin = -2.4; xmax  = 1.2;
ymin = -1.5; ymax  = 1.5;

% Generate x and y coordinates and z complex values
[x,y]=meshgrid(xmin:scale:xmax,ymin:scale:ymax);
z = x+1i*y;

% Generate w accumulation matrix and k counting matrix
w = zeros(size(z));
k = zeros(size(z));

N = 0;
while N<Nmax && ~all(k(:))
    w = w.^2+z;
    N = N+1;
    k(~k & abs(w)>4) = N;
end
k(k==0) = Nmax;
figure
s = pcolor(x, y, mod(k, 2));
colormap([0 0 0;1 1 1])
set(s,'edgecolor','none')

axis([xmin xmax -ymax ymax])
fsize=15;
set(gca,'XTick',xmin:0.4:xmax,'FontSize',fsize)
set(gca,'YTick',-ymax:0.5:ymax,'FontSize',fsize)
xlabel('Re z','FontSize',fsize)
ylabel('Im z','FontSize',fsize)
% End of Program 3b

author: Sable

repo: mcbench-benchmarks


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%

% even or odd
 t1=0:4
 t2=0:-1:-4

 %even
 x1=t1.^2
 x2=t2.^2

%  odd
 y1=t1.^3
 y2=t2.^3

 %analysis in even and odd parts
 n=-5:5;
 u=(n>=0);
 stem(n,u);

 figure
 u_n=(n<=0)
 ue=1/2*(u+ u_n);
 stem(n,ue);

 figure
 uo=1/2*(u- u_n);
 stem(n,uo);

 figure
 stem(n,ue+uo)


text:
function [Q,H,M4,D,A,G,int] = hessian_squared(V,F,varargin)
  % HESSIAN_SQUARED Construct a matrix to compute the integrated squared
  % Hessian of a function over a triangle mesh using the mixed finite element
  % method.
  %
  % [Q] = hessian_squared(V,F)
  %
  % Inputs:
  %   V  #V by dim list of vertex positions
  %   F  #F by 3 list of triangle mesh indices
  %   extraint #any boundary vertices that are explicitly to be treated as
  %   the interior
  % Outputs:
  %   Q  #V by #V sparse matrix so that X'*Q*X measures the integrated squared
  %     Hessian energy of a scalar function X
  %

  extraints = [];
  if(length(varargin)>0)
      extraints = varargin{1};
  end

  % just special case curves
  if size(F,2) == 2
    E = F;
    n = size(V,1);
    int = find(accumarray(E(:),1,[n 1])==2);
    b = find(accumarray(E(:),1,[n 1])==1);
    M = massmatrix(V,E);
    H = cotmatrix(V,E);
    % Should make this #V by #V with zeros to match other size(F,2) cases
    H(b,:) = 0;
    Q = H'*(M\H);
    M4 = M;
    return;
  end

  % Number of faces
  m = size(F,1);
  % Number of vertices
  n = size(V,1);
  % Number of dimension
  dim = size(V,2);
  m = size(F,1);
  n = size(V,1);
  G = grad(V,F);
  assert(size(G,1) == dim*m,'Gradient should equal dim*m');
  M = massmatrix(V,F);
  M4 = repdiag(M,dim^2);
  % Block transpose of G
  GG = sparse(m,0);
  for d = 1:dim
    GG = [GG G((d-1)*m+(1:m),:)];
  end
  D = repdiag(GG,dim);

  %% This is not necessary. The gradient operator G is guaranteed by
  %%construction to generate vectors in the face-plane. Therefore, even if the
  %%matrix divergence lies off this plane, the dot product with the gradient will
  %%precisely ignore the off-plane component.
  % project_div = true;
  % if project_div
  %   warning('projectin''');
  %   N = normalizerow(normals(V,F));
  %   DN = sparse(repmat(1:m,3,1)',reshape(1:m*3,m,3),N,m,m*3);
  %   D = (speye(m*3)-DN'*DN)*D;
  % end

  switch size(F,2)
  case 3
    A = repdiag(diag(sparse(doublearea(V,F)*0.5)),dim);
  case 4
    A = repdiag(diag(sparse(volume(V,F))),dim);
  end
  H = D'*A*G;

  b = unique(boundary_faces(F));
  int = setdiff(1:n,b);
  int = [int(:); extraints];
  int4 = (0:(dim^2-1))*n + int;
  notint4 = setdiff(1:(n*dim^2), int4);
  %M4 = M4(int4,int4);
  %H = H(int4,:);
  H(notint4,:) = 0;
  Q = H'*(M4\H);
end

author: alecjacobson

repo: gptoolbox


text:
function U_out = updateSubspace(U_in, residual, w, param)
%  update U

    eta = param.eta;
    residualNorm = norm(residual);
    wNorm  = norm(w);

    sigma = param.lambda * residualNorm * wNorm;

    p = residual / residualNorm;

    q = w / wNorm;

  %  t = eta * sigma; % dynamic size

     t =eta;

    U_out = U_in + (cos(t)-1) * U_in * (q * q')  - sin(t) * p * q';

%       normU = norm(U_out)
%       diffI = norm(U_out'*U_out - eye(size(U_out,2)))
end

author: andrewssobral

repo: lrslibrary


text:
function mpe = find_mpe(engine, evidence, varargin)
% FIND_MPE Find the most probable explanation of the data (assignment to the hidden nodes)
% function mpe = find_mpe(engine, evidence,...)
%
% evidence{i} = [] if X(i) is hidden, and otherwise contains its observed value (scalar or column vector).
%
% The following optional arguments can be specified in the form of name/value pairs:
% [default value in brackets]
%
% soft    - a cell array of soft/virtual evidence;
%           soft{i} is a prob. distrib. over i's values, or [] [ cell(1,N) ]
%

bnet = bnet_from_engine(engine);
ns = bnet.node_sizes(:);
N = length(bnet.dag);

engine.evidence = evidence;

% set default params
exclude = [];
soft_evidence = cell(1,N);

% parse optional params
args = varargin;
nargs = length(args);
for i=1:2:nargs
  switch args{i},
   case 'soft',    soft_evidence = args{i+1};
   otherwise,
    error(['invalid argument name ' args{i}]);
  end
end
engine.maximize = 1;

onodes = find(~isemptycell(evidence));
hnodes = find(isemptycell(evidence));
pot_type = determine_pot_type(bnet, onodes);
 if strcmp(pot_type, 'cg')
  check_for_cd_arcs(onodes, bnet.cnodes, bnet.dag);
end

hard_nodes = 1:N;
soft_nodes = find(~isemptycell(soft_evidence));
S = length(soft_nodes);
if S > 0
  assert(pot_type == 'd');
  assert(mysubset(soft_nodes, bnet.dnodes));
end

% Evaluate CPDs with evidence, and convert to potentials
pot = cell(1, N+S);
for n=1:N
  fam = family(bnet.dag, n);
  e = bnet.equiv_class(n);
  if isempty(bnet.CPD{e})
    error(['must define CPD ' num2str(e)])
  else
    pot{n} = convert_to_pot(bnet.CPD{e}, pot_type, fam(:), evidence);
  end
end

for i=1:S
  n = soft_nodes(i);
  pot{N+i} = dpot(n, ns(n), soft_evidence{n});
end
clqs = engine.clq_ass_to_node([hard_nodes soft_nodes]);

[clpot, seppot] = init_pot(engine, clqs, pot, pot_type, onodes);
[clpot, seppot] = collect_evidence(engine, clpot, seppot);
mpe = find_max_config(engine, clpot, seppot, evidence); % instead of distribute evidence

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [solutions, goodRxns]=randomSampling(model,nSamples,replaceBoundsWithInf,supressErrors,showProgress,goodRxns,minFlux)
% randomSampling
%   Returns a number of random solutions
%
%   model                   a model structure
%   nSamples                the number of solutions to return
%                           (opt, default 1000)
%   replaceBoundsWithInf    replace the largest upper bounds with Inf and
%                           the smallest lower bounds with -Inf. This is
%                           needed in order to get solutions without loops
%                           if your model has for example 1000/-1000 as
%                           arbitary large bounds. If your model only has
%                           "biologically relevant" bounds, then set this
%                           to false (opt, default true)
%   supressErrors           the program will halt if it has problems
%                           finding non-zero solutions which are not
%                           involved in loops. This could be because the
%                           constraints on the model are too relaxed (such
%                           as unlimited glucose uptake) or too strict
%                           (such as too many and too narrow constraints)
%                           (opt, default false)
%   showProgress            if true, it will display in the command window
%                           how many iterations have been done (opt, default
%                           false)
%   goodRxns                double vector of indexes of those reactions
%                           that are not involved in loops and can be used
%                           as random objective functions, as generated by
%                           a previous run of randomSampling on the same
%                           model (opt, default empty)
%   minFlux                 determines if a second optimization should be
%                           performed for each random sample, to minimize
%                           the number of fluxes and thereby preventing
%                           loops. Typically, loops are averaged out when a
%                           large number of samples are taken, but this is
%                           not always the case (opt, default false)
%
%   solutions               matrix with the solutions
%   goodRxns                double vector of indexes of those reactions
%                           that are not involved in loops and can be used
%                           as random objective functions
%
%   The solutions are generated by maximizing (with random weights) for a
%   random set of three reactions. For reversible reactions it randomly
%   chooses between maximizing and minimizing.
%
%   Usage: solutions=randomSampling(model,nSamples,replaceBoundsWithInf,supressErrors,showProgress,goodRxns,minFlux)

if nargin<2
    nSamples=1000;
end
if nargin<3
    replaceBoundsWithInf=true;
end
if nargin<4
    supressErrors=false;
end
if nargin<5
    showProgress=false;
end
if nargin<7
    minFlux=false;
end

nRxns=2; %Number of reactions in the objective function in each iteration

%First check that the model is feasible given the constraints
[sol,~]=solveLP(model);
if isempty(sol.x)
    EM='The model has no feasible solution, likely due to incompatible constraints';
    dispEM(EM);
elseif sol.f==0 && showProgress
    warning('The model objective function cannot reach a non-zero value. This might be intended, so randomSampling will continue, but this could indicate problems with your model')
end

%Simplify the model to speed stuff up a little. Keep original mapping
originalRxns=model.rxns;
model=simplifyModel(model,false,false,true,true);

%Then change the bounds to +/- Inf. This is needed in order to not have
%loops in the solutions
if replaceBoundsWithInf==true
    model.ub(model.ub==max(model.ub))=Inf;
    if min(model.lb)<0 % Only negative lower bounds should be set to -Inf
        model.lb(model.lb==min(model.lb))=-Inf;
    end
end

[~,hsSol]=solveLP(model);
%Reactions which can be involved in loops should not be optimized for.
%Check which reactions reach an arbitary high upper bound
if nargin<6 || isempty(goodRxns)
    goodRxns=true(numel(model.rxns),1);
    if showProgress
        fprintf('Prepare goodRxns list of reactions not involved in loops...   0%% complete');
    end
    for i=1:numel(model.rxns)
        if showProgress && rem(i,100) == 0
            progress=num2str(floor(100*i/numel(model.rxns)));
            progress=pad(progress,3,'left');
            fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b%s%% complete',progress);
        end
        if goodRxns(i)==true
            testModel=setParam(model,'eq',model.rxns(i),1000);
            sol=solveLP(testModel,0,[],hsSol);
            if ~isempty(sol.f)
                goodRxns(abs(sol.x)>999)=false;
            else
                %If the reaction is reversible, also check in that direction
                if model.rev(i)
                    testModel=setParam(model,'eq',model.rxns(i),-1000);
                    sol=solveLP(testModel,0,[],hsSol);
                    if ~isempty(sol.f)
                        goodRxns(abs(sol.x)>999)=false;
                    end
                end
            end
        end
    end
    if showProgress
        fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b100%% complete\n');
    end
    goodRxns=find(goodRxns);
end

%Reserve space for a solution matrix
sols=zeros(numel(model.rxns),nSamples);

%Main loop
counter=1;
badSolutions=0;
if showProgress
    itrStr=num2str(nSamples);
    fprintf('Performing random sampling: ready with%s0/%s iterations',repmat(' ',1,length(itrStr)),itrStr);
end
while counter<=nSamples
    if showProgress && rem(counter,100) == 0
        fprintf([repmat('\b',1,length(itrStr)*2) '\b\b\b\b\b\b\b\b\b\b\b\b%s/%s iterations'],pad(num2str(counter),length(itrStr),'left'),itrStr);
    end
    rxns=randsample(numel(goodRxns),nRxns);
    model.c=zeros(numel(model.rxns),1);
    multipliers=randsample([-1 1],nRxns,true);
    multipliers(model.rev(goodRxns(rxns))==0)=1;
    model.c(goodRxns(rxns))=rand(nRxns,1).*multipliers;
    if true(minFlux)
        sol=solveLP(model,1,[],hsSol);
    else
        sol=solveLP(model,0,[],hsSol);
    end
    if any(sol.x) && abs(sol.f)>10^-8
        sols(:,counter)=sol.x;
        counter=counter+1;
        badSolutions=0;
    else
        badSolutions=badSolutions+1;
        %If it only finds bad solutions then throw an error.
        if badSolutions==100 && supressErrors==false
            error('The program is having problems finding non-zero solutions (ignoring reactions that might be involved in loops). Review the constraints on your model. Set supressErrors to true to ignore this error');
        end
    end
end
if showProgress
    fprintf('\n')
end

%Map to original model
[~, I]=ismember(model.rxns,originalRxns);
solutions=zeros(numel(originalRxns),nSamples);
solutions(I,:)=sols;
solutions=sparse(solutions);
end

%To use instead of the normal Matlab randsample function. This is in order
%to not depend on the Matlab statistical toolbox.
function I=randsample(n,k,replacement)
if nargin<3
    replacement=false;
end
%n can be a integer, which leads to I being sampled from 1:n, or it can be
%a population to sample from.
if numel(n)==1 && isnumeric(n)
    n=1:n;
end
%Loop and get random numbers until the list is unique. This is only a good
%option is the number of samples is small compared to the population. There
%are several checks that should be made here, for example regarding size
%and that the number of samples is <=population size if replacement==false.
%This is not the case in randomSampling, so such checks are ignored
while true
    J=randi(numel(n),[k,1]);
    if replacement==true || numel(J)==numel(unique(J))
        I=n(J);
        break;
    end
end
I=I(:);
end

author: SysBioChalmers

repo: RAVEN


text:
% VL_RODR  Rodrigues' formula
%   R = VL_RODR(OM) where OM a 3-dimensional column vector computes the
%   Rodrigues' formula of OM, returning the rotation matrix R =
%   expm(vl_hat(OM)).
%
%   [R,DR] = VL_RODR(OM) computes also the derivative of the Rodrigues
%   formula. In matrix notation this is the expression
%
%           d(vec expm(vl_hat(OM)) )
%     dR = ----------------------.
%                  d om^T
%
%   [R,DR]=VL_RODR(OM) when OM is a 3xK matrix repeats the operation for
%   each column (or equivalently matrix with 3*K elements). In this
%   case R and DR are arrays with K slices, one per rotation.
%
%   See also: VL_IRODR(), VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function self = retrieve(dataformat, varargin)
%CATALOG.RETRIEVE Read seismic events from common file formats & data sources.
% retrieve can read events from many different earthquake catalog file
% formats (e.g. Seisan, Antelope) and data sources (e.g. IRIS DMC) into a
% GISMO Catalog object.
%
% Usage:
%       catalogObject = CATALOG.RETRIEVE(dataformat, 'param1', _value1_, ...
%                                                   'paramN', _valueN_)
%
% dataformat may be:
%
%   * 'iris' (for IRIS DMC, using irisFetch.m),
%   * 'antelope' (for a CSS3.0 Antelope/Datascope database)
%   * 'seisan' (for a Seisan database with a REA/YYYY/MM/ directory structure)
%   * 'zmap' (converts a Zmap data strcture to a Catalog object)
%
% The name-value parameter pairs supported are the same as those supported
% by irisFetch.Events(). Currently these are:
%
%     startTime
%     endTime
%     eventId
%     fetchLimit
%     magnitudeType
%     minimumLongitude
%     maximumLongitude
%     minimumLatitude
%     maximumLatitude
%     minimumMagnitude
%     maximumMagnitude
%     minimumDepth
%     maximumDepth
%
% And the two convenience parameters:
%
% radialcoordinates = [ centerLatitude, centerLongitude, maximumRadius ]
%
% boxcoordinates = [ minimumLatitude maximumLatitude minimumLongitude maximumLongitude ]
%
% For examples, see Catalog_cookbook. Also available at:
% https://geoscience-community-codes.github.io/GISMO/tutorials/html/Catalog_cookbook.html
%
%
% See also CATALOG, IRISFETCH, CATALOG_COOKBOOK

% Author: Glenn Thompson (glennthompson1971@gmail.com)

%% To do:
% Implement name-value parameter pairs for all methods
% Test the Antelope method still works after factoring out db_load_origins
% Test the Seisan method more
% Add in support for 'get_arrivals'

    debug.printfunctionstack('>')

    switch lower(dataformat)
        case 'iris'
            if exist('irisFetch.m','file')
                    ev = irisFetch.Events(varargin{:});
                    self = Catalog.read_catalog.iris(ev);
            else
                warning('Cannot find irisFetch.m')
            end
        case {'css3.0','antelope', 'datascope'}
            if admin.antelope_exists()
                self = Catalog.read_catalog.antelope(varargin{:});
            else
                warning('Sorry, cannot read event Catalog from Antelope database as Antelope toolbox for MATLAB not found')
                self = Catalog();
            end
        case 'seisan'
            self = Catalog.read_catalog.seisan(varargin{:});
        case 'aef'
            self = Catalog.read_catalog.aef(varargin{:});
        case 'sru'
            self = Catalog.read_catalog.sru(varargin{:});
        case 'vdap'
            self = Catalog.read_catalog.vdap(varargin{:});
        case 'zmap'
            self = Catalog.read_catalog.zmap(varargin{:});
        otherwise
            self = NaN;
            fprintf('format %s unknown\n\n',data_source);
    end
    if isempty(self)
        self=Catalog();
    end

    debug.printfunctionstack('<')
end

author: geoscience-community-codes

repo: GISMO


text:
function cl = extract_roi_averages(obj, mask, varargin)
% This image_vector method a extracts and averages data stored in an fmri_data object
% from a set of ROIs defined in a mask.
% It is *slightly* different from the fmri_data method, as fmri_data has
% more fields.
%
% This version resamples the mask_image to be in the same space as the obj if needed.
%
% Regions to average over can be either regions of contiguous voxels
% bounded by voxels with values of 0 or NaN, which are considered non-data
% values, or regions defined by unique integer codes in the mask image
% (i.e., for atlas images with unique codes for each defined region.)
%
% Mask/Atlas image does NOT have to be in the same space as the images to
% extract from.  It will be remapped/resliced.
%
% Extracted data is returned in single data format.
%
% :Usage:
% ::
%
%    cl = extract_roi_averages(image_vector obj, mask, [average_over])
%
% :Inputs:
%
%   - char array of strings containing 4D image file names (data extracted from these)
%   - mask_image to extract from.
%
% :Optional inputs:
%
%   **average_over:**
%      - Default: 'contiguous_regions' to average over contiguous voxels
%        bounded by voxels of 0 or NaN (non-data values)
%      - Alt. option = 'unique_mask_values' to average over unique integer codes in the mask image
%        (i.e., for atlas images with unique codes for each defined region)
%
% :Examples:
% ::
%
%    imgs_to_extract_from = filenames('w*.nii','char');
%    mask_image = which('anat_lbpa_thal.img');
%    [cl, imgdat] = extract_image_data(imgs_to_extract_from, mask_image);
%
% :See also:
%
% For an non-object-oriented alternative, see extract_image_data.m
% fmri_data.extract_roi_averages, apply_parcellation.m

% Programmers' notes:
% This function is different from fmri_data.extract_roi_averages
% Better to have only one function of record in the future...
% Note:
% cl = extract_roi_averages(imgs, atlas_obj{1});
% accomplishes the same task as apply_parcellation, returns slightly different values due to interpolation.

doverbose = 1;
if any(strcmp(varargin, 'noverbose'))
    doverbose = 0;
end

if doverbose
    fprintf('image_vector.extract_roi_averages: ');
end

if ~isa(mask, 'image_vector')
    mask = image_vector('image_names', mask);
end

obj = replace_empty(obj);
mask = replace_empty(mask);

% --------------------------------------------
% Resample space if needed
% --------------------------------------------

isdiff = compare_space(obj, mask);

if isdiff == 0 || isdiff == 3 % Same space, diff voxels
   % do nothing
elseif isdiff == 2
    error('Invalid object: volInfo structure missing for one or more objects');
else
    mask = resample_space(mask, obj);
end

% --------------------------------------------
% Redefine elements of new mask to apply
% to keep only coords in original dataset
% --------------------------------------------

is_inmask = obj.volInfo.image_indx & mask.volInfo.image_indx;
mask.volInfo.image_indx = is_inmask;

orig_wh_inmask = mask.volInfo.wh_inmask; % save for later -

mask.volInfo.wh_inmask = find(mask.volInfo.image_indx);
mask.volInfo.n_inmask = length(mask.volInfo.wh_inmask);

% need to limit coords in new mask to ONLY those in the orig mask as
% well.
wh_to_keep = is_inmask(orig_wh_inmask);
mask.volInfo.xyzlist = mask.volInfo.xyzlist(wh_to_keep, :);
mask.volInfo.cluster = mask.volInfo.cluster(wh_to_keep, :);

% --------------------------------------------
% Redefine data (temporarily; not passed out)
% to keep only coords in new mask
% --------------------------------------------

% need to limit data in obj to ONLY voxels that are also in the new
% mask. Index in space of original fmri_data object for which to keep:
wh_to_keep = is_inmask(obj.volInfo.wh_inmask);

% eliminate out-of-mask voxels before indexing into them with new mask
 % transpose, cols are voxels
obj.dat = obj.dat(wh_to_keep, :)';

if doverbose
    fprintf('\n');
end

% ---------------------------------
% define region object based on choices
% ---------------------------------

average_over = 'unique_mask_values'; %'contiguous_regions'  or 'unique_mask_values';

for varg = 1:length(varargin)
    if ischar(varargin{varg})
        switch varargin{varg}

            % reserved keywords
            case 'contiguous_regions', average_over = 'contiguous_regions';
            case 'unique_mask_values', average_over = 'unique_mask_values';

            case 'pattern_expression', error('pattern expression only defined for fmri_data.extract_roi_averages.');

            case 'noverbose'  % do nothing; already done

            otherwise
                disp('fmri_data.extract_roi_averages: Illegal string value for average_over.');
                fprintf('You entered ''%s''\n Valid values are %s or %s\n', varargin{varg}, '''contiguous_regions''', '''unique_mask_values''');
                error('Exiting');
        end
    end
end

if isa(mask, 'atlas') % If atlas object
    cl = atlas2region(mask);                % handles a few things better
else
    cl = region(mask, average_over);
end

cl(1).source_images = obj.fullpath;

% ---------------------------------
% Now get averages by cluster
% ---------------------------------

if doverbose
    fprintf('Averaging data. ');
end

switch size(mask.dat, 1)
    case length(mask.volInfo.wh_inmask)
        %maskData = mask.dat(logical(mask.volInfo.wh_inmask), :);
        % We already have in-mask data; do not select
        maskData = mask.dat;

    case length(mask.volInfo.image_indx)
        maskData = mask.dat(logical(mask.volInfo.image_indx), :);
    otherwise
        error('mask .dat is wrong size.')
end

switch average_over

    % Define integer codes for sets of voxels to average over.

    case 'unique_mask_values'
        maskData = round(maskData);
        u = unique(maskData)'; u(u == 0) = [];
        nregions = length(u);
        if doverbose, fprintf('Averaging over unique mask values, assuming integer-valued mask: %3.0f regions\n', nregions); end

    case 'contiguous_regions'
        u = unique(mask.volInfo.cluster); u(u == 0) = [];
        maskData = mask.volInfo.cluster;

    case 'none'
        cl = [];
        return

    otherwise
        error('Illegal value for average_over.  See help for this function.');
end

% Now get the average activity in each region for the defined regions

nregions = length(u);

if nregions ~= length(cl)
    disp('Num of regions in mask does not equal num of clusters to extract from.')
    disp('The most likely cause is an ill-formed/outdated obj.volInfo.cluster field');
    disp('in the region-defining object.  Try obj = reparse_contiguous(obj);');
    disp('Stopping in debugger so you can check/debug.');
    keyboard
end

for i = 1:nregions
    imgvec = maskData == u(i);

    regiondat = obj.dat(:, imgvec); % transposed above, cols are voxels. images x voxels

    if ~isempty(regiondat)
        if size(regiondat, 2) == 1              % only one image, mean is image val
            regionmean = double(regiondat);     % row vector
        else
            regionmean = double(nanmean(regiondat, 1)); % mean over images
        end

        cl(i).all_data = single(regiondat);     % images x voxels

    else
        regionmean = NaN .* zeros(size(dat, 1), 1);
    end

    cl(i).dat = regionmean;

end

if doverbose
    fprintf('Done.\n');
end

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function cost = DLSI_term(D, D_range)
% * Syntax: cost = DLSI_term(D, D_range)
% * Calculating the structured incoherence term in DLSI [[5]](#fn_dls).
% * $\sum_{c=1}^C \sum_{i \neq c} \|D_i^TD_c\|_F^2$
% -----------------------------------------------
% Author: Tiep Vu, thv102@psu.edu, 5/11/2016
%         (http://www.personal.psu.edu/thv102/)
% -----------------------------------------------
    if nargin == 0
        d = 100;
        D_range = 10*(0:10);
        D = rand(d, D_range(end));
    end
    %% MAIN
    A = erase_diagonal_blocks(D'*D, D_range, D_range);
    cost = normF2(A);
    %%
    if nargin == 0
        cost = [];
    end
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function varargout = drawLine3d(lin, varargin)
%DRAWLINE3D Draw a 3D line clipped by the current axes
%
%   drawLine3d(LINE) draws the line LINE on the current axis, by clipping
%   with the current axis.
%
%   drawLine3d(LINE, PARAM, VALUE) accepts parameter/value pairs, like
%   for plot function. Color of the line can also be given as a single
%   parameter.
%
%   drawLine3d(AX,...) plots into AX instead of GCA.
%
%   H = drawLine3d(...) returns a handle to the created line object.
%   If the line is not clipped by the axis, function returns -1.
%
%   See also:
%   lines3d, createLine3d, clipLine3d
%
% ---------
% author : David Legland
% INRA - TPV URPOI - BIA IMASTE
% created the 17/02/2005.
%

% Parse and check inputs
isLine3d = @(x) validateattributes(x,{'numeric'},...
    {'nonempty','nonnan','real','finite','size',[nan,6]});
defOpts.Color='b';
[hAx, lin, varargin]=...
    parseDrawInput(lin, isLine3d, 'line', defOpts, varargin{:});

% extract limits of the bounding box
lim = get(hAx, 'xlim');
xmin = lim(1);
xmax = lim(2);
lim = get(hAx, 'ylim');
ymin = lim(1);
ymax = lim(2);
lim = get(hAx, 'zlim');
zmin = lim(1);
zmax = lim(2);

% clip the line with the limits of the current axis
edge = clipLine3d(lin, [xmin xmax ymin ymax zmin zmax]);

% draw the clipped line
if sum(isnan(edge))==0
    h  = drawEdge3d(hAx, edge);
    if ~isempty(varargin)
        set(h, varargin{:});
    end
else
    h  = -1;
end

% process output
if nargout>0
    varargout{1}=h;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
%% _getAllFiles_ demo
% The _getAllFiles_ function recursively collects a list of files from a
% folder tree, allowing you to specify the selection criteria for which
% files are collected and how they are formatted. For these examples, we'll
% be using the main MATLAB toolbox path.
%%

rootPath = 'C:\Program Files\MATLAB\R2016b\toolbox\matlab';
format compact

%% The 'FileFilter' option
% We can specify a
% <https://www.mathworks.com/help/matlab/matlab_prog/regular-expressions.html
% regular-expression> pattern to filter the file names on, collecting a
% list of files that match. Here's an example that uses the 'FileFilter'
% option to recursively find every file with a '.m' extension:

fileList = getAllFiles(rootPath, 'FileFilter', '\.m$');
fprintf('%d files found.\n', size(fileList, 1));
fprintf('%s\n', fileList{1:5}, '...');

%%
% It's a pretty long list, so I've only shown the first 5 files it finds.
% Notice they are listed with the full path prepended by default.

%%
% If you have multiple match expressions to filter on, you can use
% <https://www.mathworks.com/help/matlab/matlab_prog/regular-expressions.html#f0-42884
% grouping operators> to include them all in one expression. For example,
% this will find every '.jpg', '.bmp', and '.tif' file:

fileList = getAllFiles(rootPath, 'FileFilter', '\.(jpg|bmp|tif)$');
fprintf('%d files found.\n', size(fileList, 1));
fprintf('%s\n', fileList{1:5}, '...');

%% The 'Depth' option
% If we don't want to search quite so far down the folder tree, we can
% limit the search depth with the 'Depth' option. Let's see how many '.m'
% files are in the root folder:

fileList = getAllFiles(rootPath, 'FileFilter', '\.m$', 'Depth', 0);
fprintf('%d files found.\n', size(fileList, 1));

%%
% Looks like none are. They are all contained in subfolders. Let's see how
% many are located in just the immediate subfolders of the root folder:

fileList = getAllFiles(rootPath, 'FileFilter', '\.m$', 'Depth', 1);
fprintf('%d files found.\n', size(fileList, 1));

%% The 'PrependPath' option
% Maybe we just want the file names, but don't care about the absolute
% paths. In this case, we just set the 'PrependPath' option to _false_:

fileList = getAllFiles(rootPath, 'FileFilter', '\.m$', ...
                                 'PrependPath', false);
fprintf('%d files found.\n', size(fileList, 1));
fprintf('%s\n', fileList{1:5}, '...');

%% The 'ValidateFcn' option
% Sometimes we might want to select files based on a more complicated
% criteria than just what's in their names. In this case, we can use the
% 'ValidateFcn' option to specify a function that is to be run on each file
% found. This function should accept as input a structure of the form
% returned by the |dir| function and return a logical value (|true| to
% collect it in the list, |false| to ignore it). First, let's find all the
% '.png' files:

fileList = getAllFiles(rootPath, 'FileFilter', '\.png$');
fprintf('%d files found.\n', size(fileList, 1));

%%
% Now, we can specify an anonymous function that gets the byte size of each
% file and returns |true| for only those greater than 250KB:

bigFcn = @(s) (s.bytes > 512^2);
fileList = getAllFiles(rootPath, 'FileFilter', '\.png$', ...
                                 'ValidateFcn', bigFcn);
fprintf('%s\n', fileList{:});

%%
% Just the one.

author: jtkim-kaist

repo: VAD

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


text:
%==============================================================================
% Copyright (C) 2006, Jan Modersitzki and Nils Papenberg, see copyright.m;
% this file is part of the FLIRT Package, all rights reserved,
% http://www.math.uni-luebeck.de/SAFIR/FLIRT-MATLAB.html
%==============================================================================
% this is only for output and does not need any explaination

function varargout = plotGrid(y,Omega,m,varargin);

if length(m) > 2,
  % 3D !!
  varargout = {[]};
  return;
end;

y1 = reshape(y(1:prod(m)),m);
y2 = reshape(y(1+prod(m):end),m);

grid1 = 1;
grid2 = 1;
color = 'r';

for k=1:1:length(varargin)/2,
  %disp([varargin{2*k-1},'=varargin{',int2str(2*k),'};']);
  eval([varargin{2*k-1},'=varargin{',int2str(2*k),'};']);
end;

K1 = 1:grid1:m(1); %K1 = [K1,m(1)];
K2 = 1:grid2:m(2); %K2 = [K2,m(2)];

% and plot the grid
hold on
l = 0;
for j=1:length(K2),  l=l+1; pp(l) = plot(y1(:,K2(j)),y2(:,K2(j)));  end;
for j=1:length(K1),  l=l+1; pp(l) = plot(y1(K1(j),:),y2(K1(j),:));  end;
hold off

if length(color) == 1, color = char(color); end;
set(pp,'color',color);

if nargout == 1,
  varargout = {pp};
end;

return;
%==========================================================================

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function FaceMod(FNameOut)

[FName, FPath]=uigetfile('*.bmp','Choose a file to modify:');

img=imread(FName);
Msk=imread('MaskF.bmp')/255;
MskI=imread('MaskFi.bmp')/255;

imgMsk=img.*Msk;
imgMsk=img.*MskI;

Norm=100/mean(mean(imgMsk));
imgNorm=imgMsk*Norm;
disp(mean(mean(imgNorm)));
imwrite(imgNorm(:,:),sprintf('%s%s%s%s%s%s',FNameOut,'O','1','_F','0','.bmp'),'BMP');
imwrite(flipdim(imgNorm(:,:),1),sprintf('%s%s%s%s%s%s',FNameOut,'O','2','_F','0','.bmp'),'BMP');

FSize=[1 2 4 8 16 32 64 128];
for i=1:8;

    h = fspecial('average', [FSize(i),FSize(i)]);
    imgFiltMsk=imfilter(imgNorm,h).*Msk;
    Norm=100/mean(mean(imgFiltMsk));
    imgFiltNorm=imgFiltMsk*Norm;
    disp(mean(mean(imgFiltNorm)));

    imwrite(imgFiltNorm,sprintf('%s%s%s%s%s%s',FNameOut,'O','1','_F',num2str(i),'.bmp'),'BMP');
    imwrite(flipdim(imgFiltNorm,1),sprintf('%s%s%s%s%s%s',FNameOut,'O','2','_F',num2str(i),'.bmp'),'BMP');
end

author: bids-standard

repo: bids-examples


text:
function newSteps = dtiWarpStep(origSteps,param)
% Warp (shift and scale) the steps into a new coord system
%
%    newSteps = dtiWarpStep(origSteps,param)
%
%    Shifts and scales according to: f(x,p) = ( x+p(2) ) * p(1)
%
%    origSteps - nx1 vector
%    param     - 2x1 vector [scale shift]
%
%    newSteps  - nx1 vector
%
% Example:
%    % original x's
%    x = 1:10;
%    % stretch by 2 times after shifting 5 units to the right
%    nx = dtiWarpStep(x,[2 5]);
%
% History:
%    2007/01/17 shc wrote it.
%

if ieNotDefined('origSteps'), error('Require input data!');         end
if ieNotDefined('param'),     error('Require warping parameters!'); end

newSteps = (origSteps+param(2))*param(1);

return


text:
function gX = whiteblockKernGradX(kern, X, X2)

% WHITEBLOCKKERNGRADX Gradient of WHITEBLOCK kernel wrt input locations.
% FORMAT
% DESC computes the gradident of the white noise block kernel with respect
% to the input positions where both the row positions and column positions
% are provided separately.
% ARG kern : kernel structure for which gradients are being computed.
% ARG x1 : row locations against which gradients are being computed.
% ARG x2 : column locations against which gradients are being computed.
% RETURN g : the returned gradients. The gradients are returned in
% a matrix which is numData2 x numInputs x numData1. Where numData1 is
% the number of data points in X1, numData2 is the number of data
% points in X2 and numInputs is the number of input
% dimensions in X.
%
% SEEALSO whiteblockKernParamInit, kernGradX, whiteblockKernDiagGradX
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if nargin<3
    X2 = X;
end

gX = cell(1, kern.nout);
for i=1:kern.nout
    gX{i} = zeros(size(X2, 1), size(X2, 2), size(X, 1));
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function varargout = clipMeshVertices(v, f, b, varargin)
%CLIPMESHVERTICES Clip vertices of a surfacic mesh and remove outer faces.
%
%   [V2, F2] = clipMeshVertices(V, F, B)
%   Clip a mesh represented by vertex array V and face array F, with the
%   box represented by B. The result is the set of vertices contained in
%   the box, and a new set of faces corresponding to original faces with
%   all vertices within the box.
%
%   [V2, F2] = clipMeshVertices(..., 'shape', 'sphere') Specify the shape.
%   Default is 'box'. But it is also possible to use 'sphere' or 'plane'.
%
%   [V2, F2] = clipMeshVertices(..., 'inside', false) removes the inner
%   faces instead of the outer faces.
%
%   [V2, F2] = clipMeshVertices(..., 'trimMesh', TF)
%   Also specifies if the isolated vertices need to be removed (TF=true) ot
%   not (TF=false). Default is false.
%
%
%   Example
%     [v, f] = createSoccerBall;
%     f = triangulateFaces(f);
%     box = [0 2 -1 2 -.5 2];
%     [v2, f2] = clipMeshVertices(v, f, box, 'inside', false);
%     figure('color','w'); view(3); axis equal
%     drawMesh(v, f, 'faceColor', 'none', 'faceAlpha', .2);
%     drawBox3d(box)
%     drawMesh(v2, f2, 'faceAlpha', .7);
%
%   See also
%   meshes3d, clipPoints3d
%

% ------
% Author: David Legland, oqilipo
% e-mail: david.legland@inra.fr
% Created: 2011-04-07,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% if input is given as a structure, parse fields
if isstruct(v)
    if nargin > 2
        varargin = [b, varargin];
    end
    b = f;
    f = v.faces;
    v = v.vertices;
end

parser = inputParser;
validStrings = {'box', 'sphere', 'plane'};
addParameter(parser, 'shape', 'box', @(x) any(validatestring(x, validStrings)));
addParameter(parser, 'inside', true, @islogical);
addParameter(parser, 'trimMesh', false, @islogical);
parse(parser, varargin{:});

% clip the vertices
[v2, indVertices] = clipPoints3d(v, b,...
    'shape', parser.Results.shape, 'inside', parser.Results.inside);

% create index array for face indices relabeling
refInds = zeros(size(indVertices));
for i = 1:length(indVertices)
    refInds(indVertices(i)) = i;
end

% select the faces with all vertices within the box
if isnumeric(f)
    % Faces given as numeric array
    indFaces = sum(~ismember(f, indVertices), 2) == 0;
    f2 = refInds(f(indFaces, :));

elseif iscell(f)
    % Faces given as cell array
    nFaces = length(f);
    indFaces = false(nFaces, 1);
    for i = 1:nFaces
        indFaces(i) = sum(~ismember(f{i}, indVertices), 2) == 0;
    end
    f2 = f(indFaces, :);

    % re-label indices of face vertices (keeping horizontal index array)
    for i = 1:length(f2)
        f2{i} = refInds(f2{i})';
    end
end

if parser.Results.trimMesh
    [v2, f2] = trimMesh(v2, f2);
end

varargout = formatMeshOutput(nargout, v2, f2);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [map2] = unifyMetabolicPPImapCD(map)
% Unify a metabolic and protein-protein interaction map as a standard.
% Reaction will be grey and Metabolites/Complexes will be White.
%
% USAGE:
%
%   [map2] = unifyMetabolicPPImapCD(map)
%
% INPUT:
%   map:    MATLAB structure of CellDesigner map
%
% OUTPUT:
%   map2:   Map with grey reactions colour, width 1 and white nodes colour.
%           Change colour complex to white
%
% .. Authors:
%       - J.Modamio LCSB, Belval, Luxembourg. 19.08.2017 LCSB. Belval. Luxembourg
%       - N.Sompairac - Institut Curie, Paris, 20/10/2017

    map2 = map;

    % Change reaction colour to ligth grey and and width to 1.0
    for i = 1:length(map2.rxnName)
        map2.rxnColor{i, 1} = 'FFDCDCDC';
        map2.rxnWidth{i, 1} = 1;
    end
    % Change molecules color to white
    for j = 1:length(map2.molColor)
        map2.molColor{j, 1} = 'FFFFFFFF';
    end
    % Change complex colour to white
    for c = 1:length(map2.complexColor)
        map2.complexColor{c, 1} = 'FFFFFFFF';
    end

    % Change reactions secondary links to light grey (products)
    % and width to 1
    index = find(~cellfun(@isempty, map2.rxnProductLineColor));
    for t = 1:length(index)
        if length(map2.rxnProductLineColor{t, 1}) > 1
            for d = 1:length(map2.rxnProductLineColor{t, 1})
                map2.rxnProductLineColor{t, 1}{d, 1} = 'FFDCDCDC';
                map2.rxnProductLineWidth{t, 1}{d, 1} = 1;
            end
        elseif length(map2.rxnProductLineColor{t, 1}) == 1
            map2.rxnProductLineColor{t, 1}{1, 1} = 'FFDCDCDC';
            map2.rxnProductLineWidth{t, 1}{1, 1} = 1;
        end
    end
    clear id index

    % Change reactions secondary links to light grey (reactants)
    % and width to 1
    index = find(~cellfun(@isempty, map2.rxnReactantLineColor));
    for r = 1:length(index)
        if length(map2.rxnReactantLineColor{r, 1}) > 1
            for p = 1:length(map2.rxnReactantLineColor{r, 1})
                map2.rxnReactantLineColor{r, 1}{p, 1} = 'FFDCDCDC';
                map2.rxnReactantLineWidth{r, 1}{p, 1} = 1;
            end
        elseif length(map2.rxnReactantLineColor{r, 1}) == 1
            map2.rxnReactantLineColor{r, 1}{1, 1} = 'FFDCDCDC';
            map2.rxnReactantLineWidth{r, 1}{1, 1} = 1;
        end
    end
    clear id index

end

author: opencobra

repo: cobratoolbox


text:
function subNo = nodesSub(nodeNo,wt)

if(any(nodeNo>numel(wt.nodes)))
   error('%s: Invalid node index range. Number of nodes is %d.\n',upper(mfilename),numel(wt.nodes));
end

nodeNoa = cellfun(@(nEl) nEl.a,wt.nodes(nodeNo),'UniformOutput',0);
nodeNoUps = nodesFiltUps(nodeNo,wt);

nodesCount = numel(nodeNo);
subNo = cell(1,nodesCount);
for ii=1:nodesCount
   subNo{ii} = nodeNoUps(ii)*nodeNoa{ii};
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
img = rgb2gray(im2double(imread('corner2.gif')));
[map r c] = susanCorner(img);
figure,imshow(img),hold on
plot(c,r,'o')

author: Sable

repo: mcbench-benchmarks


text:
function r82row_order_type_test ( )

%*****************************************************************************80
%
%% R82ROW_ORDER_TYPE_TEST tests R82ROW_ORDER_TYPE.
%
%  Discussion:
%
%    An R82ROW is a (2,N) array of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 4;
  test_num = 10;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R82ROW_ORDER_TYPE_TEST\n' );
  fprintf ( 1, '  R82ROW_ORDER_TYPE classifies an R8VEC as\n' );
  fprintf ( 1, '  -1: no order\n' );
  fprintf ( 1, '   0: all equal;\n' );
  fprintf ( 1, '   1: ascending;\n' );
  fprintf ( 1, '   2: strictly ascending;\n' );
  fprintf ( 1, '   3: descending;\n' );
  fprintf ( 1, '   4: strictly descending.\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    [ x, seed ] = r8mat_uniform_01 ( 2, n, seed );

    x(1:2,1:n) = round ( 3.0 * x(1:2,1:n) );

    order = r82row_order_type ( n, x );

    string = sprintf ( '  Order type = %d\n', order );

    r82row_print ( n, x, string );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% OPT09_RUN
%
%  Modified:
%
%    08 January 2008
%
   %---------------------------------------------------------------------
   %  Running the trigonometric function
   %  This is used to test globalization methods.
   %---------------------------------------------------------------------
   fprintf('---------------------------------------------------------\n')
   fprintf('Running testcase_9:  F(X*) = 0\n')
   fprintf('---------------------------------------------------------\n')
   fname = 'opt09_fgh';
   options = [];
   options.verbose            = 0;
   options.method             = 'newton';

   fprintf('Newton:\n')
   options.globalization      = 'line_search';
   x0 = [ .25; .25; .25; .25 ];
   x = entrust(fname, x0, options);
   fprintf('Line search produced  (%8.5e,%8.5e,%8.5e,%8.5e)\n\n',...
     x(1),x(2),x(3),x(4))
   f = opt09_fgh ( x, 'f' );

   fprintf('Value of F(X) = %f\n', f );
   fprintf('Newton:\n')
   options.globalization      = 'trust_region';
   x = entrust(fname, x0, options);
   fprintf('Trust-region produced (%8.5e,%8.5e,%8.5e,%8.5e)\n\n',...
     x(1),x(2),x(3),x(4))
   f = opt09_fgh ( x, 'f' );
   fprintf('Value of F(X) = %f\n', f );

   %---------------------------------------------------------------------
   %  Test Gauss-Newton strategies.
   %---------------------------------------------------------------------
   fprintf('---------------------------------------------------------\n')
   fprintf('Running testcase_9 as least squares problem: \n')
   fprintf(' ||RES(X*)|| = 0\n')
   fprintf('---------------------------------------------------------\n')
   fname = 'opt09_rj';
   options = [];
   options.verbose            = 0;
   options.method             = 'gauss_newton';
   options.step_tolerance     = 1.e-15;
   options.globalization      = 'none';
   options.gradient_tolerance = 1.e-10;
   options.max_iterations     = 800;

   x0 = [ .25; .25; .25; .25 ];
   x = entrust(fname, x0, options);

   fprintf('Gauss-Newton produced (%8.5e,%8.5e,%8.5e,%8.5e)\n\n',...
     x(1),x(2),x(3),x(4))
   [ res, jac ] = opt09_rj ( x, 'f' );
   fprintf('Norm of RES(X) = %f\n', norm ( res ) );


text:
function dX = dynamics(X,P)
% dX = dynamics(X,P)
%
% This function computes the first-order dynamics for the classical
% "three-body" problem in newtonian mechanics: three point masses with
% inverse-square gravity.
%
% The input state contains only information about the first two particles.
% The state of the third particle is computed from the assumption (without
% loss of generality) that the center of mass of the system is at the
% origin and not moving.
%
% INPUTS:
%   X = [12, N] = first order state vector
%       X(1:3,:) = X1 = [x1;y1;z1]
%       X(4:6,:) = X2 = [x2;y2;z2]
%       X(7:9,:) = V1 = [dx1;dy1;dz1]
%       X(10:12,:) = V2 = [dx2;dy2;dz2]
%   P = parameter struct:
%       .G = gravity constant
%       .m1 = mass one
%       .m2 = mass two
%       .m3 = mass three
%
% OUTPUTS:
%   dX = [V1;V2; dV1; dV2];
%
%

% Positions
X1 = X(1:3,:);
X2 = X(4:6,:);

% Velocity
V1 = X(7:9,:);
V2 = X(10:12,:);

% Compute X3 from definition that CoM is at the origin
%   x = (x1*m1 + x2*m2 + x3*m3)/(m1 + m2 + m3);
X3 = -(X1*P.m1 + X2*P.m2)/P.m3;
% V3 = -(V1*P.m1 + V2*P.m2)/P.m3;   %Unused, but good to know

% Force Vectors
F12 = getForce(X1, X2, P.m1, P.m2, P.G);  % Acting on 1 from 2
F23 = getForce(X2, X3, P.m2, P.m3, P.G);  % Acting on 2 from 3
F13 = getForce(X1, X3, P.m1, P.m3, P.G);  % Acting on 1 from 3

% Acceleration Vectors
A1 = (F12 + F13)/P.m1;
A2 = (-F12 + F23)/P.m2;
% A3 = (-F23 - F13)/P.m3;   %Unused, but good to know

% Pack up solution:
dX = [V1; V2; A1; A2];

end

function F = getForce(Xa, Xb, ma, mb, G)
%
% Computes the force acting on particle a due to particle b, given
% parameters in the struct P
%

r = Xa - Xb;    %Vector from b to a
r2 = ones(3,1)*sum(r.^2,1);   % Length squared

lFl = G*ma*mb./r2;  %Force magnitude

F = -lFl .* (r./sqrt(r2));   %Force vector

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function issym = optiIsSymmetric(H)
%Returns true if a matrix is symmetric, using MATLAB routine if available

%Check if MATLAB version is available (>= R2014a)
if(~isempty(which('issymetric')))
    issym = issymmetric(H);
else %jc method
    if(nnz(triu(H,1) - tril(H,-1).') == 0)
        issym = true;
    else
        issym = false;
    end
end


text:
function res = dbcsme(ya,yb)
% derivative of F respect to ya
Ba=zeros(5,5);
Ba(1,1)=1;Ba(2,2)=1;Ba(3,3)=1;Ba(4,4)=1;
res.Ba=Ba;
% derivative of F respect to yb
Bb=zeros(5,5);
Bb(5,3)=-1;Bb(5,5)=1;
res.Bb=Bb;
end


text:
classdef MultiObjectiveEGO < ALGORITHM
% <multi> <real/integer> <constrained/none> <expensive>
% Multi-objective efficient global optimization
% alpha --- 0.7 --- portion of samples for Kriging construction
% num_k ---   5 --- number of infill points per iteration
% H     ---  21 --- number of reference directions

%------------------------------- Reference --------------------------------
% R. Hussein, K. Deb, A Generative Kriging Surrogate Model for Constrained
% and Unconstrained Multi-objective Optimization, in: Proc. Genet. Evol.
% Comput. Conf. 2016, Denver, 2016, 573-580.
%--------------------------------------------------------------------------

% This function is written by Youwei He (email: 1554748356@qq.com)

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            [alpha,num_k,H] = Algorithm.ParameterSet(0.7,5,21);
            % parameter for AASF in equation (10)
            rho = 1e-3;
            %% Generate the initial design points
            % number of design variables
            D = Problem.D;
            % number of initial design points 11*D-1
            N = Problem.N;
            %% step1-1: generate initial design points using Latin Hypercube sampling
            PopDec = repmat(Problem.upper-Problem.lower,N,1).*UniformPoint(N,D,'Latin') ...
                +repmat(Problem.lower,N,1);
            %% step1-2: evaluate initial design points
            Population = Problem.Evaluation(PopDec);

            %% step 2: Generate the reference direction set
            [R,N_R] = UniformPoint(H,Problem.M);
            R = R./sqrt(sum(R.^2,2));
            %% Optimization
            while Algorithm.NotTerminated(Population)
                %% step 2 diversity_preserver procedure: neighborhood approach
                i_direction = linspace(1,N_R,N_R);
                for i = i_direction
                    for j = 1 : num_k
                        Algorithm.NotTerminated(Population);
                        PopDecT = Population.decs;
                        PopObjT = Population.objs;
                        PopConT = Population.cons;% cons>=0: meet constraints
                        num_con = size(PopConT,2);
                        N_Pop=length(PopDecT(:,1));
                        % determine whether the problem is constrained or not
                        if all(all(PopConT==0)) && N_Pop==1
                            constrained=1;
                        else
                            constrained=0;
                        end
                        %% step 3: Points_Selector procedure
                        normW = sqrt(sum(R(i,:).^2,2));
                        normP = sqrt(sum(PopObjT.^2,2));
                        CosineP = sum((PopObjT).*R(i,:),2)./normP./normW;
                        % orthogonal distance of each point to the given reference direction
                        distB = normP.*sqrt(1-CosineP.^2);
                        [~,indx] = sort(distB);
                        N = ceil(alpha*N_Pop);
                        PopObj = PopObjT(indx(1:N),:);
                        PopDec = PopDecT(indx(1:N),:);
                        PopCon = PopConT(indx(1:N),:);
                        %% step 4-1
                        index = sum(PopCon >= 0, 2) == num_con;% feasible solutions
                        PopCon(index,:)=0;
                        PopCon(~index,:)=-PopCon(~index,:);
                        PopObjScaled = (PopObj-repmat(min(PopObj),N,1))./(repmat(max(PopObj),N,1)-repmat(min(PopObj),N,1));
                        if constrained
                            PopConScaled = (PopCon-repmat(min(PopCon),N,1))./(repmat(max(PopCon),N,1)-repmat(min(PopCon),N,1));
                        end
                        %% step 4-2
                        PopSmetric = zeros(N,1);
                        PopSmetricT = PopObjScaled(index,:)./repmat(R(i,:),sum(index),1);
                        PopSmetric(index,:)=max(PopSmetricT,[],2)+rho*sum(PopSmetricT,2);
                        if constrained
                            ASF_max = max(PopSmetric(index,1));
                            PopSmetric(~index,:) = ASF_max.*ones(sum(~index),1) + sum(PopConScaled(~index,:),2);
                        end
                        %% step 4-3
                        kriging_obj= dacefit(PopDec,PopSmetric,'regpoly0','corrgauss',1*ones(1,D),0.001*ones(1,D),1000*ones(1,D));
                        f_min = min(PopSmetric);
                        %% step 5: Optimization
                        infill_criterion = @(x)Infill_Standard_EI(x, kriging_obj, f_min);
                        best_x=rGA(infill_criterion,Problem);
                        % infill point too close (not used in the paper,
                        % but this happens sometimes especially for low dimensional problem)
                        if min(sqrt(sum((PopDecT-best_x).^2,2)))<1E-8
                            best_x = rGA(@(x)Infill_Maximal_Distance(x, PopDecT),Problem);
                        end
                        Population = [Population,Problem.Evaluation(best_x)];
                    end
                end
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [x y]=bresenham(x1,y1,x2,y2)

%Matlab optmized version of Bresenham line algorithm. No loops.
%Format:
%               [x y]=bham(x1,y1,x2,y2)
%
%Input:
%               (x1,y1): Start position
%               (x2,y2): End position
%
%Output:
%               x y: the line coordinates from (x1,y1) to (x2,y2)
%
%Usage example:
%               [x y]=bham(1,1, 10,-5);
%               plot(x,y,'or');
x1=round(x1); x2=round(x2);
y1=round(y1); y2=round(y2);
dx=abs(x2-x1)
dy=abs(y2-y1)
steep=abs(dy)>abs(dx)
if steep t=dx;dx=dy;dy=t; end

%The main algorithm goes here.
if dy==0
    q=zeros(dx+1,1);
else
    q=[0;diff(mod([floor(dx/2):-dy:-dy*dx+floor(dx/2)]',dx))>=0]
end

%and ends here.

if steep
    if y1<=y2 y=[y1:y2]'; else y=[y1:-1:y2]'; end
    if x1<=x2 x=x1+cumsum(q);else x=x1-cumsum(q); end
else
    if x1<=x2 x=[x1:x2]'; else x=[x1:-1:x2]'; end
    if y1<=y2 y=y1+cumsum(q);else y=y1-cumsum(q); end
end


text:
function x = dif1cyclic_null_left ( m, n )

%*****************************************************************************80
%
%% DIF1CYCLIC_NULL_LEFT returns a left null vector of the DIF1CYCLIC matrix.
%
%  Discussion:
%
%    (1,1,1,...,1) is always a null vector.
%
%    If M is even,
%
%    (A,B,A,B,A,B,...,A,B) is also a null vector, for any A and B.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the order of A.
%
%    Output, real X(M,1), the null vector.
%
  x = zeros ( m, 1 );

  if ( mod ( m, 2 ) ~= 0 )
    x(1:m,1) = 1.0;
  else
    a = 1.0;
    b = 2.0;
    x(1:2:m-1,1) = a;
    x(2:2:m,  1) = b;
  end

  return
end


text:
function grains = smooth(grains,iter,varargin)
% constraint laplacian smoothing of grain boundaries
%
% Input
%  grains - @grain2d
%  iter   - number of iterations (default: 1)
%
% Output
%  grains - @grain2d
%
% Options
%  moveTriplePoints  - do not exclude triple/quadruple points from smoothing
%  moveOuterBoundary - do not exclude outer boundary from smoothing
%  second_order, S2  - second order smoothing
%  rate              - default smoothing kernel
%  gauss             - gaussian smoothing kernel
%  exp               - exponential smoothing kernel
%  umbrella          - umbrella smoothing kernel

if nargin < 2 || isempty(iter), iter = 1; end

% compute incidence matrix vertices - faces
I_VF = [grains.boundary.I_VF,grains.innerBoundary.I_VF];

% compute vertice adjacency matrix
A_V = I_VF * I_VF';
t = size(A_V,1);

% do not consider triple points
if check_option(varargin,'moveTriplePoints')
  ignore = false(size(A_V,1),1);
else
  ignore = full(diag(A_V)) > 2;
end

% ignore outer boundary
if ~check_option(varargin,'moveOuterBoundary')
  ignore(grains.boundary.F(any(grains.boundary.grainId==0,2),:)) = true;
end

if check_option(varargin,{'second order','second_order','S','S2'})
  A_V = logical(A_V + A_V*A_V);
  A_V = A_V - diag(diag(A_V));
end

weight = get_flag(varargin,{'gauss','expotential','exp','umbrella','rate'},'rate');
lambda = get_option(varargin,weight,.5);

V = full(grains.V);
isNotZero = ~all(~isfinite(V) | V == 0,2) & ~ignore;

for l=1:iter
  if ~strcmpi(weight,'rate')
    [i,j] = find(A_V);
    d = sqrt(sum((V(i,:)-V(j,:)).^2,2)); % distance
    switch weight
      case 'umbrella'
        w = 1./(d);
        w(d==0) = 1;
      case 'gauss'
        w = exp(-(d./lambda).^2);
      case {'expotential','exp'}
        w = lambda*exp(-lambda*d);
    end

    A_V = sparse(i,j,w,t,t);
  end

  % take the mean over the neigbours
  Vt = A_V*V;

  m = sum(A_V,2);

  dV = V(isNotZero,:)-bsxfun(@rdivide,Vt(isNotZero,:),m(isNotZero,:));

  isZero = any(~isfinite(dV),2);
  dV(isZero,:) = 0;

  V(isNotZero,:) = V(isNotZero,:) - lambda*dV;

end

grains.V = V;

author: mtex-toolbox

repo: mtex

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


text:
function [sMriT1, errMsg] = bst_normalize_mni(T1File, Method, T2File)
% BST_NORMALIZE_MNI: Compute deformation fields to the MNI ICBM152 space.
%
% USAGE:  [sMriT1, errMsg] = bst_normalize_mni(T1File, Method='maff8', T2File=[])
%         [sMriT1, errMsg] = bst_normalize_mni(sMriT1, Method='maff8', sMriT2=[])
%                            bst_normalize_mni('install')               % Only installs default SPM tpm.nii
%
% INPUTS:
%    - T1File : Relative path to a T1 MRI file in the Brainstorm database
%    - sMriT1 : Brainstorm T1 MRI structure
%    - Method : String defining the method to use for the registration
%               'maff8'   : SPM mutual information algorithm (affine transform)
%               'segment' : SPM12 segment
%    - T2File : Relative path to a T2 MRI file in the Brainstorm database
%    - sMriT2 : Brainstorm T2 MRI structure

% @=============================================================================
% 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-2023

%% ===== PARSE INPUTS =====
% Inializations
global GlobalData;
errMsg = [];
% Usage: bst_normalize_mni('install')
if isequal(T1File, 'install')
    isInstall = 1;
    sMriT1 = [];
else
    isInstall = 0;
    % Usage: bst_normalize_mni(sMriT1)
    if ~ischar(T1File)
        sMriT1 = T1File;
        T1File = [];
    % Usage: bst_normalize_mni(T1File)
    else
        sMriT1 = [];
    end
end
% Default method
if (nargin < 2) || isempty(Method)
    Method = 'maff8';
end
% T2 MRI
if (nargin < 3) || isempty(T2File)
    T2File = [];
    sMriT2 = [];
else
    if ~ischar(T2File)
        sMriT2 = T2File;
        T2File = [];
    else
        sMriT2 = [];
    end
end

%% ===== GET SPM TEMPLATE =====
% Open progress bar
isProgress = bst_progress('isVisible');
if ~isProgress
    bst_progress('start', 'MNI normalization', 'Initialization...');
end
bst_plugin('SetProgressLogo', 'spm12');
% Get template file
TpmFile = bst_get('SpmTpmAtlas');
% If it does not exist: download
if isempty(TpmFile) || ~file_exist(TpmFile)
    % Create folder
    if ~file_exist(bst_fileparts(TpmFile))
        mkdir(bst_fileparts(TpmFile));
    end
    % URL to download
    tmpUrl = 'http://neuroimage.usc.edu/bst/getupdate.php?t=SPM_TPM';
    % Path to downloaded file
    tpmZip = bst_fullfile(bst_get('BrainstormUserDir'), 'defaults', 'spm', 'SPM_TPM');
    % Download file
    errMsg = gui_brainstorm('DownloadFile', tmpUrl, tpmZip, 'Download template');
    % Error message
    if ~isempty(errMsg)
        errMsg = ['Impossible to download template:' 10 errMsg];
        return;
    end
    % Progress bar
    bst_progress('text', 'Importing SPM template...');
    % URL: Download zip file
    try
        unzip(tpmZip, bst_fileparts(tpmZip));
    catch
        errMsg = ['Could not unzip anatomy template: ' 10 10 lasterr];
        disp(['BST> Error: ' errMsg]);
        file_delete(tpmZip, 1);
        if ~isProgress
            bst_progress('stop');
        end
        return;
    end
    % Delete zip file
    file_delete(tpmZip, 1);
    % Get template file
    TpmFile = bst_get('SpmTpmAtlas');
    if isempty(TpmFile) || ~file_exist(TpmFile)
        errMsg = 'Missing file TPM.nii';
        return;
    end
end
% If only installing: exit
if isInstall
    return;
end

%% ===== LOAD ANATOMY =====
% T1 MRI
if isempty(sMriT1)
    % Progress bar
    bst_progress('text', 'Loading input MRI...');
    % Check if it is loaded in memory
    [sMriT1, iLoadedMri] = bst_memory('GetMri', T1File);
    % If not: load it from the file
    if isempty(sMriT1)
        sMriT1 = in_mri_bst(T1File);
    end
else
    iLoadedMri = [];
end
% T2 MRI
if isempty(sMriT2) && ~isempty(T2File)
    sMriT2 = in_mri_bst(T2File);
end

%% ===== MNI NORMALIZATION =====
TpmFiles = [];
try
    switch (Method)
        % SPM12 LINEAR MUTUAL INFORMATION
        case 'maff8'
            % Progress bar
            bst_progress('text', 'Resampling MRI...');
            % Resample volume if needed
            if any(abs(sMriT1.Voxsize - [1 1 1]) > 0.001)
                [sMriRes, Tres] = mri_resample(sMriT1, [256 256 256], [1 1 1], 'linear');
            else
                sMriRes = sMriT1;
                Tres = [];
            end
            % Compute affine transformation to MNI space
            Tmni = mri_normalize_maff(sMriRes, TpmFile);
            % Append the resampling transformation matrix
            if ~isempty(Tres)
                Tmni = Tmni * Tres;
            end
            % Save results into the MRI structure
            sMriT1.NCS.R = Tmni(1:3,1:3);
            sMriT1.NCS.T = Tmni(1:3,4);

        % SPM12 SEGMENT
        case 'segment'
            % Initialize SPM
            [isInstalled, errMsg] = bst_plugin('Install', 'spm12');
            if ~isInstalled
                return;
            end
            % Progress bar
            bst_progress('text', 'Running SPM batch... (see command window)');
            % Compute non-linear registration to MNI space
            [sMriT1, TpmFiles] = mri_normalize_segment(sMriT1, TpmFile, sMriT2);
            if isempty(sMriT1)
                errMsg = 'SPM Segment failed.';
                return;
            end
    end
catch
    errMsg = ['bst_normalize_mni/' Method ': ' lasterr()];
    sMriT1 = [];
    return;
end

%% ===== SAVE RESULTS =====
bst_progress('text', 'Saving normalization...');
% Compute default fiducials positions based on MNI coordinates
sMriT1 = mri_set_default_fid(sMriT1, 'maff8');
% Save modifications in the MRI file
if ~isempty(T1File)
    bst_save(file_fullpath(T1File), sMriT1, 'v6');
end

%% ===== UPDATE LOADED FIGURES =====
% If the MRI is currently loaded
if ~isempty(iLoadedMri)
    % Update structures
    GlobalData.Mri(iLoadedMri).NCS.AC = sMriT1.NCS.AC;
    GlobalData.Mri(iLoadedMri).NCS.PC = sMriT1.NCS.PC;
    GlobalData.Mri(iLoadedMri).NCS.IH = sMriT1.NCS.IH;
    GlobalData.Mri(iLoadedMri).NCS.Origin = sMriT1.NCS.Origin;
    % Linear MNI transformation
    if isfield(sMriT1.NCS,'R') && isfield(sMriT1.NCS,'T')
        GlobalData.Mri(iLoadedMri).NCS.R  = sMriT1.NCS.R;
        GlobalData.Mri(iLoadedMri).NCS.T  = sMriT1.NCS.T;
    end
    % Non-linear MNI transformtion
    if isfield(sMriT1.NCS,'y') && isfield(sMriT1.NCS,'iy') && isfield(sMriT1.NCS,'y_vox2ras')
        GlobalData.Mri(iLoadedMri).NCS.y         = sMriT1.NCS.y;
        GlobalData.Mri(iLoadedMri).NCS.iy        = sMriT1.NCS.iy;
        GlobalData.Mri(iLoadedMri).NCS.y_vox2ras = sMriT1.NCS.y_vox2ras;
    end
    % SCS coordinates
    GlobalData.Mri(iLoadedMri).SCS.R   = sMriT1.SCS.R;
    GlobalData.Mri(iLoadedMri).SCS.T   = sMriT1.SCS.T;
    GlobalData.Mri(iLoadedMri).SCS.NAS = sMriT1.SCS.NAS;
    GlobalData.Mri(iLoadedMri).SCS.LPA = sMriT1.SCS.LPA;
    GlobalData.Mri(iLoadedMri).SCS.RPA = sMriT1.SCS.RPA;
    GlobalData.Mri(iLoadedMri).SCS.Origin = sMriT1.SCS.Origin;
end

%% ===== TISSUE CLASSIFICATION =====
% Import tissue classification
if ~isempty(TpmFiles) && ~isempty(T1File)
    bst_progress('text', 'Loading tissue segmentations...');
    % Get subject
    [sSubject, iSubject] = bst_get('MriFile', T1File);
    % Import tissue classification
    import_mri(iSubject, TpmFiles, 'SPM-TPM', 0, 1, 'tissues_segment');
    % Delete the temporary folder
    file_delete(bst_fileparts(TpmFiles{1}), 1, 1);
end

% Close progress bar
bst_plugin('SetProgressLogo', []);
if ~isProgress
    bst_progress('stop');
end

author: brainstorm-tools


text:
function dat = one_class_loss(algo,dat)

  [x y]=get_xy(dat);

  y = ones(size(x));
  lss=sum((x~=y))/length(y);

  dat=data([get_name(dat) ' -> one_class_loss='  num2str(lss,4) ],[],lss);


text:
%% Point Cloud Viewer (OpenCV)
%
% Display point cloud image using OpenCV.
%
% Copyright 2012 The MathWorks, Inc.

%% Input Signals
%
% * Pause (boolean): Pause point cloud image (1: Pause / 0: Not pause).
%
% * Pos (int32): Start position of XYZ and Image input pixel data ([y x]).
% Note that Pos is zero index based.
%
% * XYZ (double): XYZ data of points.
%
% * Image (uint8): RGB data of points.
%
% * View (double): View point ([x y z]). View input is available when "Set view point by input" block parameter was checked.
%

%% Output Signals
%
% * Image (uint8): RGB point cloud image. Image is available when "Image output" block parameter was checked.
%

%% Block Parameters
%
% * Window size: Size of Point Cloud Viewer window
%
% * Frame size: Size of XYZ and Image input data ([rows colums]).
%
% * Full frame size: Size of whole XYZ and Image data ([rows colums]).
%
% * Background color: Background color (R G B).
%
% * Look at point: Look at point of point cloud image ([x y z]).
%
% * Set view point by Input: Configure view point is set by an input signal or keyboard input (ON: Input signal / OFF: Keyboard input).
%
% Note that when this parameter was OFF, view point can be changed by keyboard input while point cloud display window is active.
%
% Note that keyboard input is not case sensitive.
%
% - Press R(r) : Initial view point
%
% - Press K(k) : +y
%
% - Press J(j) : -y
%
% - Press H(h) : +x
%
% - Press L(l) : -x
%
% - Press I(i)  : +z
%
% - Press M(m): -z
%
% * Initial view point: Initial view point when "Set view point by Input" parameter was OFF ([x y z]).
%
% * View point step: View point step for keyboard input. Step value has to be greater than zero.
%
% * Camera distortion coefficients: The input vector of camera distortion coefficients ([k1 k2 p1 p2 k3]) Zero value means no distortion.
%
% * Image output: Configure output of point cloud image (ON: Image output/OFF:No image output).
%


text:
% More complex example

% Written by Benjamin Irving 2013/03/25

% CREATING THE MATLAB EXAMPLE SURFACE
[X, Y, Z] = peaks;
Z=Z/2;
colormap hsv;

% example plot in matlab (just to illustrate the matlab output)
surf(X,Y,Z)
axis equal;

% CONVERT SURFACE TO A PATCH (faces and vertices)
fvc=surf2patch(X, Y, Z, Z, 'triangles');

%CONVERT THE COLORMAP INTO A COLORMAP FOR THE VERTICES (optional if you
%want to include color)
cmap=colormap;
%getting cdata
cdata=fvc.facevertexcdata;
%normalising cdata rangle
cdata=(cdata-min(cdata));cdata=cdata/max(cdata);
%scaling by cmap range and converting to integers
cdata=round(cdata.*(size(cmap,1)-1));
%getting cdata values
cdata=cmap(cdata+1,:);

% Create a second set of vertices with z positions inverted (as an example
% of a deformation)
fvc2=fvc;
fvc2.vertices(:,3)=-fvc2.vertices(:,3);

% EXPORT THE MESH TO HTML USING X3MESH
x3mesh_deform(fvc.faces, fvc.vertices, fvc2.vertices, 'name', 'Example2', 'color', cdata, 'speed', 5)


text:
% This is a demonstration of inpainting using the Potts model.
% Note that it is an intrinsic property of the Potts model to
% create two triple junctions instead of a four-junction in the center.

% load image
img = double(imread('colors.png'))/255;

% set inpainting mask
simg = sum(img, 3);
weights = simg ~= 0;

% add noise
rng(123) % set seed value of random number gen. for reproducibility
imgNoisy = img + 0.3 * randn(size(img)) .* cat(3, weights, weights, weights);

%% Potts restoration
gamma = 2;
clear opts;
opts.weights = weights;
opts.verbose = true;
opts.muStep = 1.1;
opts.muInit = 1e-4;
tic
u = minL2Potts2DADMM(imgNoisy, gamma, opts);
toc

%% Show result
subplot(1,2,1)
imshow(imgNoisy)
title('Noisy data (Black pixels are missing')
subplot(1,2,2)
imshow(u)
title('Potts inpainting')

author: mstorath

repo: Pottslab

save_path: github-repos/MATLAB/mstorath-Pottslab


text:
function r = roots(f, varargin)
%ROOTS   Roots of a CLASSICFUN in the interval [a,b].
%   ROOTS(F) returns the real roots of the CLASSICFUN F in the interval [a,b].
%
%   ROOTS(F, OPTIONS) modifies the default ROOTS properties, by passing the
%   OPTIONS to the rootfinding method of the ONEFUN of F.
%
%   If F is an array-valued CLASSICFUN then there is no reason to expect each column to
%   have the same number of roots. In order to return a useful output, the roots
%   of each column are computed and then padded with NaNs so that a matrix may
%   be returned. The columns of R = ROOTS(F) correspond to the columns of F.
%
% See also ONEFUN. [TODO]: Why may we also want to see ONEFUN?

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

% Deal with empty case:
if ( isempty(f) )
    r = [];
    return
end

% Find the roots of the ONEFUN of F:
onefunRoots = roots(f.onefun, varargin{:});

% Map the roots found on [-1,1] to the interval [a,b]:
r = f.mapping.For(onefunRoots);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


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

%%begin

% Frequently we want to define a smooth sequence of positions (or poses) from
% one point to another.  First consider the 1-dimensional case.
%
% We define the start and end position

p0 = -1;
p1 = 2;

% and a smooth path from p0 to p1 in 50 time steps is given by

p = tpoly(p0, p1, 50);
about p
% which we see has 50 rows.  We can plot this

plot(p)

% and see that it does indeed move smoothly from p0 to p1 and that the initial
% and final derivative (and second derivative) is zero.

% We can also get the velocity and acceleration

[p,pd,pdd] = tpoly(p0, p1, 50);
subplot(3,1,1); plot(p); xlabel('Time'); ylabel('p');
subplot(3,1,2); plot(pd); xlabel('Time'); ylabel('pd');
subplot(3,1,3); plot(pdd); xlabel('Time'); ylabel('pdd');

% This path is a 5th order polynomial and it suffers from the disadvantage that
% the velocity is mostly below the maximum possible value.  An alternative is

[p,pd,pdd] = lspb(p0, p1, 50);
subplot(3,1,1); plot(p); xlabel('Time'); ylabel('p');
subplot(3,1,2); plot(pd); xlabel('Time'); ylabel('pd');
subplot(3,1,3); plot(pdd); xlabel('Time'); ylabel('pdd');
% which we see has a trapezoidal velocity profile.

% Frequently the start and end values are vectors, not scalars, perhaps a 3D
% position or Euler angles.  In this case we apply the scalar trajectory function
% to a vector with

p = mtraj(@tpoly, [0 1 2], [2 1 0], 50);
about p
% and p again has one row per time step, and one column per vector dimension

clf; plot(p)

%---
% Finally, we may wish to interpolate poses.  We will define a start and end pose

T0 = transl(0.4, 0.2, 0) * trotx(pi);
T1 = transl(-0.4, -0.2, 0.3) * troty(pi/2) * trotz(-pi/2);

% and a smooth sequence between them in 50 steps is

T = ctraj(T0, T1, 50);
about T
% which is a 4x4x50 matrix.  The first pose is

T(:,:,1)

% and the 10th pose is

T(:,:,10)

% We can plot the motion of this coordinate frame by

clf; tranimate(T)

author: Allopart

repo: rbpf-gmapping