text:
function atlas = ft_read_atlas(filename, varargin)

% FT_READ_ATLAS reads an template/individual segmentation or parcellation from disk.
% The volumetric segmentation or the surface-based parcellation can either represent
% a template atlas (e.g. AAL or the Talairach Daemon), it can represent an
% individualized atlas (e.g. obtained from FreeSurfer) or it can represent an
% unlabeled parcellation/segmentation obtained from an individual's DTi, anatomical,
% or resting state fMRI scan.
%
% Use as
%   atlas = ft_read_atlas(filename, ...)
% or
%   atlas = ft_read_atlas({filenamelabels, filenamemesh}, ...)
%
% Additional options should be specified in key-value pairs and can include
%   'format'      = string, see below
%   'unit'        = string, e.g. 'mm' (default is to keep it in the native units of the file)
%   'map'         = string, 'maxprob' (default), or 'prob', for FSL-based atlases, providing
%                   either a probabilistic segmentation or a maximum a posterior probability map
%   'labelfile'   = string, point to a (generic) text or xml file for interpretation of the values in the atlas
%
% For individual surface-based atlases from FreeSurfer you should specify two
% filenames as a cell-array: the first points to the file that contains information
% with respect to the parcels' labels, the second points to the file that defines the
% mesh on which the parcellation is defined.
%
% The 'format' variable, if not specified, will be determined automatically. In general
% it will not be needed to specify it. The following formats are supported:
%
% Volumetric atlases based on a (gzipped) nifti-file with an companion txt-file for interpretation
%   'aal'               assumes filename starting with 'ROI_MNI'
%   'brainnetome'       assumes companion lookuptable txt-file starting with 'Brainnetome Atlas'
%   'simnibs_v4'        assumes filename starting with 'final_tissues', with companion freesurfer-style lookuptable txt-file
%   'wfu'               assumes specific formatting of companion lookuptable txt-file
%
% Volumetric atlases based on a (gzipped) nifti-file with hard coded assumption on the labels
%   'yeo7'
%   'yeo17'
%
% Volumetric atlases based on a folder with (gzipped) nifti-files with a companion xml-file for interpretation
%   'fsl'               assumes path to folder with data mentioned in the xml-file. Use xml-file as filename
%
% Volumetric atlases based on the freesurfer mgz format with standard lookuptable txt-file for interpretation
%   'freesurfer_volume' assumes the freesurfer LUT file for interpretation, and assumes aparc or aseg in the
%                       filename, used for subject-specific parcellations
%
% Volumetric atlases based on the afni software
%   'afni'              assumes filename containing BRIK or HEAD, assumes generic interpretation of the labels
%                       for the TTatlas+tlrc, or otherwise the interpretation should be in the file
%
% Volumetric atlas based on the spm_anatomy toolbox
%   'spm_anatomy'       pair of .hdr/.img files, and an associated mat-file for the interpretation
%                       Specify the associated mat-file with MPM in filename
%
% Surface based atlases, requiring a pair of files, containing the labels, and the associated geometry
%   'caret_label'       hcp-workbench/caret style .gii, with .label. in filename, requires additional file describing the geometry
%   'freesurfer_surface' freesurfer style annotation file, requires additional file describing the geometry
%
% Miscellaneous formats
%   'mat'               mat-file, with FieldTrip style struct, other matlab data that FieldTrip knows to handle, can also be
%                       Brainstorm derived surfaces
%   'vtpm'
%
% For volume data for whicth the format cannot be automatically detected, or if the volume data does not have a companion file
% for the interpretation of the labels, a list of 'fake' labels will be generated.
%
% The output atlas will be represented as structure according to FT_DATATYPE_SEGMENTATION or
% FT_DATATYPE_PARCELLATION.
%
% The 'lines' and the 'colorcube' colormaps may be useful for plotting the different
% patches, for example using FT_PLOT_MESH, or FT_SOURCEPLOT.
%
% See also FT_READ_MRI, FT_READ_HEADSHAPE, FT_PREPARE_SOURCEMODEL, FT_SOURCEPARCELLATE, FT_PLOT_MESH

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

% deal with multiple filenames, if the geometry and labels are stored
% in different files, as with freesurfer/caret
if isa(filename, 'cell')
  if numel(filename)==2
    filenamemesh = filename{2};
    filename     = filename{1};
  else
    ft_error('with multiple filenames, only 2 files are allowed');
  end % if precisely two input files
end % iscell

% optionally get the data from an URL and make a temporary local copy
filename = fetch_url(filename);

[p, f, x] = fileparts(filename);

labelfile = ft_getopt(varargin, 'labelfile', []);
ftype     = ft_filetype(filename);

% if the original file was a .gz
if isequal(x,'.gz')
  [p, f, x] = fileparts(filename(1:end-3));
end

% do an educated guess of the format based on the input, ensure that the
% required additional functions are available, and specify the companion labelfile
if contains(filename, 'BRIK') || contains(filename, 'HEAD')
  % the above is needed to correctly detect zipped files
  ft_hastoolbox('afni', 1);

  % this is robust for both compressed or uncompressed afni atlases.
  format = 'afni';
elseif startsWith(ftype, 'nifti')

  % first handle the ones for which the labels are implicit in the filename
  if contains(f, 'Yeo2011_7Networks')
    % assume to be conform the shared atlas, with hardcoded index-to-label mapping
    format = 'yeo7';
  elseif contains(f, 'Yeo2011_17Networks')
    % assume to be conform the shared atlas, with hardcoded index-to-label mapping
    format = 'yeo17';
  end

  % try for different conventions of the naming of a potential companion labelfile
  if isempty(labelfile)
    labelfile = fullfile(p, sprintf('%s.txt', f));
    if ~exist(labelfile, 'file')
      % in case the volume was zipped, and the txt only replaced the gz
      labelfile = fullfile(p, sprintf('%s.nii.txt', f));
    end
    if ~exist(labelfile, 'file')
      labelfile = ''; % revert back to empty
    end
  end

  if ~isempty(labelfile)
    % do a quick check on the companion labelfile to specify the format
    fid = fopen_or_error(labelfile, 'rt');
    l1  = fgetl(fid);
    if strcmp(l1(1),'[') && strcmp(l1(end),']')
      % specific to this one is that some of the wakeforest atlases contain non-human data, so no fixed coordsys can be assumed
      format   = 'wfu';
      coordsys = 'unknown';
    elseif strcmp(l1,'Brainnetome Atlas')
      % specific to this one is the coordinate system it seems
      format   = 'brainnetome';
      coordsys = 'mni'; % apparently the image is in radiological convention,
      % so the voxel-axes are left-handed. The even valued parcels should end
      % up in the right hemisphere
    elseif contains(filename, 'ROI_MNI') || contains(filename, 'AAL')
      % newer versions of AAL also exist with a companion .xml file, this is not yet supported
      format   = 'aal';
      coordsys = 'mni';
    end
    fclose(fid);
  end
  if contains(filename, 'final_tissues')
    ft_hastoolbox('freesurfer', 1); % required to read lookuptable

    % assume to be from SimNIBS version 4
    format = 'simnibs_v4';

    labelfile = fullfile(p, sprintf('%s_LUT.txt', f));
  end

  if isempty(labelfile)
    % just a nifti file without interpretation of labels
    format = 'nifti_no_label';
  end

elseif isequal(ftype, 'freesurfer_mgz') && contains(f, 'aparc') || contains(f, 'aseg')
  ft_hastoolbox('freesurfer', 1);

  % individual volume based segmentation from freesurfer
  format = 'freesurfer_volume';

  if isempty(labelfile)
    % use the version for freesurfer that is in fieldtrip/external/freesurfer
    [ftver, ftpath] = ft_version;
    labelfile  = fullfile(ftpath, 'external/freesurfer', 'FreeSurferColorLUT.txt');
  end

elseif isequal(ftype, 'freesurfer_annot')
  ft_hastoolbox('freesurfer', 1);

  % individual volume based segmentation from freesurfer
  format = 'freesurfer_surface';

elseif isequal(ftype, 'caret_label')
  ft_hastoolbox('gifti', 1);

  % this is a gifti file that contains both the values for a set of vertices as well as the labels.
  format = 'caret_label';

elseif contains(filename, 'MPM.mat')
  ft_hastoolbox('spm8up', 1);

  % assume to be from the spm_anatomy toolbox
  format = 'spm_anatomy';

elseif exist(fullfile(p, [f '_MPM.mat']), 'file')
    ft_warning('please specify the corresponding MPM.mat file as the filename');
    filename = fullfile(p, [f '_MPM.mat']); % update the filename
    [p, f, x] = fileparts(filename);

    ft_hastoolbox('spm8up', 1);

    % assume to be from the spm_anatomy toolbox
    format = 'spm_anatomy';

elseif strcmp(x, '.xml')
  ft_hastoolbox('gifti', 1); % required to read the xml

  % fsl-style atlas, this is assumed to consist of an .xml file that specifies the labels, as well as a pointer
  % to the file/folder with the volume data.
  format = 'fsl';

  labelfile = fullfile(p, sprintf('%s.xml',f));

elseif strcmp(x, '.mat')
  % mat-file to contain a well-defined structure
  format = 'mat';

else

  format = 'unknown';
end

% get the optional input arguments
fileformat = ft_getopt(varargin, 'format', format);
unit       = ft_getopt(varargin, 'unit');

switch fileformat
  case {'aal' 'brainnetome' 'freesurfer_volume' 'nifti_no_label' 'simnibs_v4' 'wfu'}
    atlas = ft_read_mri(filename, 'outputfield', 'tissue');

    % interpret the format specific labelfile
    if isequal(fileformat, 'aal')
      % The labelfile is a combination of nii+txt file, where the txt file may contain three columns like this
      %   FAG	Precentral_L	2001
      %   FAD	Precentral_R	2002
      %   ...

      fid = fopen_or_error(labelfile, 'rt');
      C = textscan(fid, '%s%s%d');
      lab = C{2};
      idx = C{3};
      fclose(fid);
    elseif isequal(fileformat, 'wfu')
      % the download from http://fmri.wfubmc.edu comes with pairs of nii and txt files
      % the text file looks like this, with tabs between the colums
      % the section at the end with three times 191 does not always exist
      %
      % [ TD Labels]
      % 53	Angular Gyrus         191 191 191
      % 39	Anterior Cingulate		191 191 191
      % ...

      fid = fopen_or_error(labelfile, 'rt');
      C = textscan(fid, '%d%s%*[^\n]', 'HeaderLines', 1, 'Delimiter', '\t');
      lab = C{2};
      idx = C{1};
      fclose(fid);
    elseif any(strcmp(fileformat, {'simnibs_v4' 'freesurfer_volume'}))
      [idx, lab, rgba] = read_fscolorlut(labelfile);
      lab = cellstr(lab);
    elseif isequal(fileformat, 'brainnetome')
      % Brainnetome Atlas: L. Fan, et al.The Human Brainnetome Atlas: A New Brain Atlas Based on
      % Connectional Architecture. Cereb Cortex 2016; 26 (8): 3508-3526. doi: 10.1093/cercor/bhw157
      fid = fopen_or_error(labelfile, 'rt');
      fgetl(fid); % this reads: 'Brainnetome Atlas'
      lab = cell(246,1);
      for i=1:246
        lab{i,1}=fgetl(fid);
      end
      idx = (1:246)';
      fclose(fid);

    elseif isequal(fileformat, 'nifti_no_label')
      % the file does not exist
      ft_warning('cannot locate a labelfile, making default tissue labels');

      idx = (1:max(atlas.tissue(:)))';
      lab = cell(size(idx));
      for i = 1:numel(lab)
        % this is consistent with FIXSEGMENTATION
        lab{i} = sprintf('tissue %d', i);
      end
    end

    if ~isfield(atlas, 'coordsys') && exist('coordsys', 'var')
      atlas.coordsys = coordsys;
    end

    uval = unique(atlas.tissue(:));
    sel  = find(ismember(idx, uval));
    fprintf('subselecting %d labels from the total list of %d\n', numel(sel), numel(lab));
    idx  = idx(sel);
    lab  = lab(sel);
    if exist('rgba', 'var')
      tmprgba = rgba(sel,:);
      rgba = zeros(0,4);
    end

    % remap the values in the data, if needed
    if ~isequal(idx(:)', 1:numel(idx))
      dat = zeros(atlas.dim);
      cnt = 0;
      for k = 1:numel(idx)
        sel = atlas.tissue==idx(k);
        if sum(sel(:))
          cnt = cnt+1;
          fprintf('re-indexing label %s to a value of %d (was %d)\n', lab{k}, cnt, idx(k));
          dat(sel)           = cnt;
          tissuelabel{cnt,1} = lab{k};
          if exist('rgba', 'var')
            rgba(cnt,:) = tmprgba(k,:);
          end
        end
      end
      atlas.tissue = dat;
    else
      tissuelabel = lab;
    end

    atlas.tissuelabel = tissuelabel;
    if exist('rgba', 'var'), atlas.rgba = rgba; end

    % reduce memory footprint
    if numel(tissuelabel)<=intmax('uint8')
      atlas.tissue = uint8(atlas.tissue);
    elseif numel(tissuelabel)<=intmax('uint16')
      atlas.tissue = uint16(atlas.tissue);
    elseif numel(tissuelabel)<=intmax('uint32')
      atlas.tissue = uint32(atlas.tissue);
    end

  case 'afni'

    tmp     = ft_read_mri(filename);
    if isfield(tmp, 'coordsys') && ~strcmp(tmp.coordsys, 'unknown')
      coordsys = tmp.coordsys;
    elseif isfield(tmp.hdr, 'TEMPLATE_SPACE') && ~isempty(tmp.hdr.TEMPLATE_SPACE)
      coordsys = lower(tmp.hdr.TEMPLATE_SPACE); % FIXME this is based on AFNI conventions, not easily decodable by FT
    else
      coordsys = 'tal'; % FIXME could be different in other atlases
    end

    if isfield(tmp.hdr, 'ATLAS_LABEL_TABLE') && ~isempty(tmp.hdr.ATLAS_LABEL_TABLE)
      if isfield(tmp.hdr.ATLAS_LABEL_TABLE(1), 'sb_label') && ~all(tmp.anatomy(:)==round(tmp.anatomy(:)))
        % probabilistic atlas
        isprobabilistic = true;
      else
        % indexed atlas
        isprobabilistic = false;
      end
      labels  = {tmp.hdr.ATLAS_LABEL_TABLE.struct}';
      values  = [tmp.hdr.ATLAS_LABEL_TABLE.val]';

    elseif contains(filename, 'TTatlas+tlrc')
      isprobabilistic = false;

      [labels, values] = TTatlas_labels;

    else
      ft_error('no information about the atlas labels is available');
    end

    atlas     = [];
    atlas.dim = tmp.dim(1:3);
    atlas.transform = tmp.transform;
    atlas.hdr = tmp.hdr;
    atlas.coordsys = coordsys;

    nbrick  = size(tmp.anatomy,4);
    for k = 1:nbrick

      if ~isprobabilistic
        brickname = sprintf('brick%d',k-1);
        brick     = tmp.anatomy(:,:,:,k);
        ulabel    = setdiff(unique(brick(:)), 0);
        label     = cell(size(ulabel));
        nlabel    = numel(label);

        % renumber the brick from 1:N and keep track of the label
        newbrick  = zeros(size(brick));
        for i = 1:nlabel
          sel = find(values==ulabel(i));
          if ~isempty(sel)
            label(i) = labels(sel);
            newbrick(brick==ulabel(i)) = i;
          else
            ft_warning('the value %d does not have a label according to the ATLAS_LABEL_TABLE and will be discarded', ulabel(i));
          end
        end
        atlas.(brickname) = newbrick;
        atlas.([brickname 'label']) = label;
      else
        atlas.(labels{k}) = tmp.anatomy(:,:,:,values(k)+1); % indexing is 0-based in this case
      end

    end

  case {'freesurfer_surface'}

    if contains(filename, 'a2009s')
      parcelfield = 'a2009s';
    elseif contains(filename, 'aparc')
      parcelfield = 'aparc';
    elseif contains(filename, 'ba')
      parcelfield = 'BA';
    else
      ft_error('unknown freesurfer parcellation type requested');
    end

    % read the labels
    [v, p, c] = read_annotation(filename);

    label = c.struct_names;
    rgba  = c.table(:,1:4);
    rgb   = c.table(:,5); % compound value that is used for the indexing in vector p
    index = ((1:c.numEntries)-1)';

    switch ft_filetype(filenamemesh)
      case 'freesurfer_triangle_binary'
        [pos, tri] = read_surf(filenamemesh);

        % ensure the triangles to be 1-indexed
        if min(tri(:))==0 && max(tri(:))==size(pos,1)-1
          tri = tri+1;
        end

        bnd.pos    = pos;
        bnd.tri    = tri;
        reindex    = true;
      otherwise
        ft_error('unsupported fileformat for surface mesh');
    end

    % check the number of vertices
    if size(bnd.pos,1) ~= numel(p)
      ft_error('the number of vertices in the mesh does not match the number of elements in the parcellation');
    end

    % reindex the parcels, if needed: I am not fully sure about this, but the caret
    % label files seem to have the stuff numbered with normal numbers, with unknown
    % being -1. assuming them to be in order;
    if reindex
      % this is then freesurfer convention, coding in rgb
      newp = zeros(size(p));
      for k = 1:numel(label)
        newp(p==rgb(k)) = index(k)+1;
      end
    else
      uniquep = unique(p);
      if uniquep(1)<0
        p(p<0) = 0;
      end
      newp   = p;
    end
    atlas       = [];
    atlas.pos   = bnd.pos;
    atlas.tri   = bnd.tri;
    atlas.(parcelfield)            = newp;
    atlas.([parcelfield, 'label']) = label;
    atlas.rgba  = rgba;
    atlas       = ft_determine_units(atlas);

  case 'caret_label'

    g = gifti(filename);

    rgba = [];
    if isfield(g, 'labels')
      label = g.labels.name(:);
      key   = g.labels.key(:);
      if isfield(g.labels, 'rgba')
        rgba = g.labels.rgba; % I'm not sure whether this always exists
      end
    else
      label = g.private.label.name(:);
      key   = g.private.label.key(:);
      if isfield(g.private.label, 'rgba')
        rgba = g.private.label.rgba; % I'm not sure whether this always exists
      end
    end

    %label = g.private.label.name; % provides the name of the parcel
    %key   = g.private.label.key;  % maps value to name

    % Store each column in cdata as an independent parcellation, because
    % each vertex can have multiple values in principle

    atlas = [];
    for k = 1:size(g.cdata,2)
      tmporig  = g.cdata(:,k);
      tmpnew   = nan(size(tmporig));
      tmplabel = cell(0,1);
      tmprgba  = zeros(0,4);
      cnt = 0;
      for m = 1:numel(label)
        sel = find(tmporig==key(m));
        if ~isempty(sel)
          cnt = cnt+1;
          if any(strcmp(tmplabel,deblank(label{m})))
            % one feature of a gifti can be that the same labels can exist (and
            % are treated as a different parcel) while they should be the same
            % parcel, i.e. when there's an empty space at the end of the label
            val = find(strcmp(tmplabel, deblank(label{m})));
          else
            % add as a new label
            tmplabel{end+1,1} = label{m};
            if ~isempty(rgba), tmprgba(end+1,:)  = rgba(m,:); end
            val = cnt;
          end
          tmpnew(tmporig==key(m)) = val;
        end
      end

      % there is some additional meta data that may be useful, but for now
      % stick to the rather uninformative parcellation1/2/3 etc.
      %
      % if strcmp(g.private.data{k}.metadata(1).name, 'Name')
      %   parcelfield = fixname(g.private.data{k}.metadata(1).value);
      % else
      %   ft_error('could not determine parcellation name');
      % end

      if size(g.cdata,2)>1
        parcelfield = ['parcellation' num2str(k)];
      else
        parcelfield = 'parcellation';
      end

      atlas.(parcelfield)           = tmpnew;
      atlas.([parcelfield 'label']) = tmplabel;
      if ~isempty(tmprgba), atlas.rgba = tmprgba; end
    end

    if exist('filenamemesh', 'var')
      tmp       = ft_read_headshape(filenamemesh);
      atlas.pos = tmp.pos;
      atlas.tri = tmp.tri;
      atlas     = ft_determine_units(atlas);
    elseif ~isfield(atlas, 'coordsys')
      atlas.coordsys = 'unknown';
    end

  case 'spm_anatomy'

    % load the map, this is assumed to be the struct-array MAP
    load(filename);
    [p,f,e]      = fileparts(filename);
    mrifilename  = fullfile(p,[strrep(f, '_MPM',''),'.img']);
    atlas        = ft_read_mri(mrifilename, 'dataformat', 'analyze_img', 'outputfield', 'tissue');
    tissue       = round(atlas.tissue); % I don't know why the values are non-integer
    label        = {MAP.name}';
    idx          = [MAP.GV]';

    % check whether all labels are present
    if numel(intersect(idx,unique(tissue(:))))<numel(idx)
      fprintf('there are fewer labels in the volume than in the list\n');
    end

    % remap the values of the labels to run from 1-numel(idx)
    newtissue = zeros(size(tissue));
    for k = 1:numel(idx)
      newtissue(tissue==idx(k)) = k;
    end
    atlas.tissue      = newtissue;
    atlas.tissuelabel = label;
    atlas.coordsys    = 'spm'; % I think this is safe to assume

    clear tissue newtissue;

  case 'fsl'
    map = ft_getopt(varargin, 'map', 'maxprob');
    switch map
      case 'prob'
        imagefile = 'imagefile';
      case 'maxprob'
        imagefile = 'summaryimagefile';
      otherwise
        error('unknown map requested');
    end

    ft_hastoolbox('gifti', 1);
    hdr = xmltree(filename);
    hdr = convert(hdr);

    % get the full path
    [p, f , x]  = fileparts(filename);

    mrifilename = sprintf('%s.nii.gz', hdr.header.images{1}.(imagefile));
    if isequal(mrifilename(1), '/')
      mrifilename = mrifilename(2:end); % remove first backslash to avoid error if no full path was given
    end

    % this uses the thresholded image
    switch map
      case 'maxprob'
        atlas        = ft_read_mri(fullfile(p, mrifilename), 'outputfield', 'tissue');
        atlas.tissuelabel = hdr.data.label(:);
        atlas.coordsys    = 'mni';
      case 'prob'
        tmp   = ft_read_mri(fullfile(p, mrifilename));
        atlas = removefields(tmp, 'anatomy');
        for m = 1:numel(hdr.data.label)
          fn = strrep(hdr.data.label{m}, ' ' , '_');
          fn = strrep(fn, '-', '_');
          fn = strrep(fn, '(', '');
          fn = strrep(fn, ')', '');

          atlas.(fn) = tmp.anatomy(:,:,:,m);
          if ~isa(atlas.(fn), 'double') && ~isa(atlas.(fn), 'single')
            % ensure that the probabilistic values are either double or
            % single precision, do single precision to save memory
            atlas.(fn) = single(atlas.(fn));
          end
          if any(atlas.(fn)(:)>1) && all(atlas.(fn)(:)<=100)
            % convert to probability values, assuming 100 to be max
            atlas.(fn) = atlas.(fn)./100;
          end
        end
        atlas.coordsys = 'mni';
        atlas.dim      = atlas.dim(1:3);
    end

  case 'mat'
    tmp = load(filename);
    if isfield(tmp, 'Vertices') && isfield(tmp, 'Atlas')
      % this applies to BrainStorm *.mat files containing cortical meshes
      atlas.pos = tmp.Vertices;
      % copy some optional fields over with a new name
      atlas = copyfields(tmp, atlas, {'Faces', 'Curvature', 'SulciMap'});
      atlas = renamefields(atlas, {'Faces', 'Curvature', 'SulciMap'}, {'tri', 'curv', 'sulc'});
      for i=1:numel(tmp.Atlas)
        name   = fixname(tmp.Atlas(i).Name);
        nlabel = numel(tmp.Atlas(i).Scouts);
        label  = cell(1, nlabel);
        index  = zeros(size(atlas.pos,1), 1);
        for j=1:nlabel
          index(tmp.Atlas(i).Scouts(j).Vertices) = j;
          label{j} = tmp.Atlas(i).Scouts(j).Label;
        end
        atlas.( name         ) = index;
        atlas.([name 'label']) = label;
      end
    elseif isfield(tmp, 'atlas')
      % this applies to most FieldTrip *.mat files
      atlas = tmp.atlas;
    elseif numel(fieldnames(tmp))==1
      % just take whatever variable is contained in the file
      fn = fieldnames(tmp);
      atlas = tmp.(fn{1});
      if isstruct(atlas)
        ft_warning('assuming that the variable "%s" in "%s" represents the atlas', fn{1}, filename);
      else
        ft_error('cannot read atlas structure from "%s"', filename);
      end
    else
      ft_error('the mat-file %s does not contain a variable called ''atlas''',filename);
    end

  case 'yeo7'
    % this uses Yeo2011_7Networks_MNI152_FreeSurferConformed1mm_LiberalMask_colin27.nii, which is
    % the 7 network parcellation from https://surfer.nmr.mgh.harvard.edu/fswiki/CorticalParcellation_Yeo2011
    % aligned to the colin27 template (skull-stripped version of single_subj_T1_1mm.nii)
    % using AFNI's 3dQwarp and 3dNwarpApply
    atlas = ft_read_mri(filename, 'outputfield', 'tissue');

    atlas.coordsys    = 'mni';
    atlas.tissuelabel = {
      '7Networks_1'
      '7Networks_2'
      '7Networks_3'
      '7Networks_4'
      '7Networks_5'
      '7Networks_6'
      '7Networks_7'
      };

    colors = [
      120 18 134;
      70 130 180;
      0 118 14;
      196 58 250;
      220 248 164;
      230 148 34;
      205 62 78
      ]; % not used

  case 'yeo17'
    % this uses Yeo2011_17Networks_MNI152_FreeSurferConformed1mm_LiberalMask_colin27.nii, which is
    % the 17 network parcellation from https://surfer.nmr.mgh.harvard.edu/fswiki/CorticalParcellation_Yeo2011
    % aligned to the colin27 template (skull-stripped version of single_subj_T1_1mm.nii) using AFNI's 3dQwarp and 3dNwarpApply
    atlas = ft_read_mri(filename, 'outputfield', 'tissue');

    atlas.coordsys    = 'mni';
    atlas.tissuelabel = {
      '17Networks_1'
      '17Networks_2'
      '17Networks_3'
      '17Networks_4'
      '17Networks_5'
      '17Networks_6'
      '17Networks_7'
      '17Networks_8'
      '17Networks_9'
      '17Networks_10'
      '17Networks_11'
      '17Networks_12'
      '17Networks_13'
      '17Networks_14'
      '17Networks_15'
      '17Networks_16'
      '17Networks_17'
      };
    colors = [
      120 18 134;
      255 0 0;
      70 130 180;
      42 204 164;
      74 155 60;
      0 118 14;
      196 58 250;
      255 152 213;
      220 248 164;
      122 135 50;
      119 140 176;
      230 148 34;
      135 50 74;
      12 48 255;
      0 0 130;
      255 255 0;
      205 62 78
      ]; % not used

  otherwise
    ft_error('unsupported format "%s"', fileformat);

end % switch fileformat

if ~isempty(unit)
  % ensure the atlas is in the desired units
  atlas = ft_convert_units(atlas, unit);
else
  % ensure the units of the atlas are specified
  try
    atlas = ft_determine_units(atlas);
  catch
    % ft_determine_units will fail for triangle-only gifties.
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TTatlas labels moved to a subfunction for readability of the above code
function [labels, values] = TTatlas_labels

% the following information is from https://sscc.nimh.nih.gov/afni/doc/misc/afni_ttatlas/index_html
values = [
  68
  71
  20
  21
  22
  24
  25
  26
  27
  28
  29
  30
  31
  32
  33
  34
  35
  36
  37
  39
  40
  41
  42
  43
  44
  45
  46
  47
  48
  49
  50
  51
  52
  70
  72
  73
  74
  75
  76
  77
  124
  125
  126
  128
  129
  130
  131
  132
  133
  134
  135
  136
  137
  138
  144
  145
  151
  146
  147
  148
  149
  81
  82
  83
  84
  85
  86
  87
  88
  89
  90
  91
  93
  94
  95
  96
  97
  98
  99
  100
  101
  102
  103
  104
  105
  106
  107
  108
  109
  110
  111
  112
  113
  114
  115
  116
  117
  118
  119
  120
  121
  122
  123
  53
  54
  55
  56
  57
  58
  59
  60
  61
  62
  63
  66
  65
  127
  64
  67
  ];

labels = {
  'Hippocampus'
  'Amygdala'
  'Posterior Cingulate'
  'Anterior Cingulate'
  'Subcallosal Gyrus'
  'Transverse Temporal Gyrus'
  'Uncus'
  'Rectal Gyrus'
  'Fusiform Gyrus'
  'Inferior Occipital Gyrus'
  'Inferior Temporal Gyrus'
  'Insula'
  'Parahippocampal Gyrus'
  'Lingual Gyrus'
  'Middle Occipital Gyrus'
  'Orbital Gyrus'
  'Middle Temporal Gyrus'
  'Superior Temporal Gyrus'
  'Superior Occipital Gyrus'
  'Inferior Frontal Gyrus'
  'Cuneus'
  'Angular Gyrus'
  'Supramarginal Gyrus'
  'Cingulate Gyrus'
  'Inferior Parietal Lobule'
  'Precuneus'
  'Superior Parietal Lobule'
  'Middle Frontal Gyrus'
  'Paracentral Lobule'
  'Postcentral Gyrus'
  'Precentral Gyrus'
  'Superior Frontal Gyrus'
  'Medial Frontal Gyrus'
  'Lentiform Nucleus'
  'Hypothalamus'
  'Red Nucleus'
  'Substantia Nigra'
  'Claustrum'
  'Thalamus'
  'Caudate'
  'Caudate Tail'
  'Caudate Body'
  'Caudate Head'
  'Ventral Anterior Nucleus'
  'Ventral Posterior Medial Nucleus'
  'Ventral Posterior Lateral Nucleus'
  'Medial Dorsal Nucleus'
  'Lateral Dorsal Nucleus'
  'Pulvinar'
  'Lateral Posterior Nucleus'
  'Ventral Lateral Nucleus'
  'Midline Nucleus'
  'Anterior Nucleus'
  'Mammillary Body'
  'Medial Globus Pallidus'
  'Lateral Globus Pallidus'
  'Putamen'
  'Nucleus Accumbens'
  'Medial Geniculum Body'
  'Lateral Geniculum Body'
  'Subthalamic Nucleus'
  'Brodmann area 1'
  'Brodmann area 2'
  'Brodmann area 3'
  'Brodmann area 4'
  'Brodmann area 5'
  'Brodmann area 6'
  'Brodmann area 7'
  'Brodmann area 8'
  'Brodmann area 9'
  'Brodmann area 10'
  'Brodmann area 11'
  'Brodmann area 13'
  'Brodmann area 17'
  'Brodmann area 18'
  'Brodmann area 19'
  'Brodmann area 20'
  'Brodmann area 21'
  'Brodmann area 22'
  'Brodmann area 23'
  'Brodmann area 24'
  'Brodmann area 25'
  'Brodmann area 27'
  'Brodmann area 28'
  'Brodmann area 29'
  'Brodmann area 30'
  'Brodmann area 31'
  'Brodmann area 32'
  'Brodmann area 33'
  'Brodmann area 34'
  'Brodmann area 35'
  'Brodmann area 36'
  'Brodmann area 37'
  'Brodmann area 38'
  'Brodmann area 39'
  'Brodmann area 40'
  'Brodmann area 41'
  'Brodmann area 42'
  'Brodmann area 43'
  'Brodmann area 44'
  'Brodmann area 45'
  'Brodmann area 46'
  'Brodmann area 47'
  'Uvula of Vermis'
  'Pyramis of Vermis'
  'Tuber of Vermis'
  'Declive of Vermis'
  'Culmen of Vermis'
  'Cerebellar Tonsil'
  'Inferior Semi-Lunar Lobule'
  'Fastigium'
  'Nodule'
  'Uvula'
  'Pyramis'
  'Culmen'
  'Declive'
  'Dentate'
  'Tuber'
  'Cerebellar Lingual'
  };

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%PREX_EIGENFACES   PRTools example on the use of images and eigenfaces

help prex_eigenfaces

echo on
              % Load all faces (may take a while)

  faces = prdataset(orl);
	faces = setprior(faces,0);     % give them equal priors
  a = gendat(faces,ones(1,40));  % select one image per class

              % Compute the eigenfaces
  w = pcam(a);

              % Display them
  newfig(1,3); show(w); drawnow;

              % Project all faces onto the eigenface space

  b = [];
  for j = 1:40
    a = seldat(faces,j);
    b = [b;a*w];
    % Don't echo loops
    echo off
  end
  echo on

              % Show a scatterplot of the first two eigenfaces
  newfig(2,3)
  scatterd(b)
  title('Scatterplot of the first two eigenfaces')

              % Compute leave-one-out error curve
  featsizes = [1 2 3 5 7 10 15 20 30 39];
  e = zeros(1,length(featsizes));
  for j = 1:length(featsizes)
    k = featsizes(j);
     e(j) = testk(b(:,1:k),1);
     echo off
  end
  echo on
              % Plot error curve
  newfig(3,3)
  plot(featsizes,e)
  xlabel('Number of eigenfaces')
  ylabel('Error')
echo off

author: marianux

repo: ecg-kit

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


text:
function x=SeamPlot(x,SeamVector)
% SEAMPLOT takes as input an image and the SeamVector array and produces
% an image with the seam line superimposed upon the input image, x, for
% display purposes.
%
% Author: Danny Luong
%         http://danluong.com
%
% Last updated: 12/20/07

value=1.5*max(x(:));
for i=1:size(SeamVector,1)
    x(i,SeamVector(i))=value;
end

author: Sable

repo: mcbench-benchmarks


text:
function id = initmumps()
%
% id = initmumps
% it returns a default matlab MUMPS structure
%
% Use help mumps_help for detailed information
%
errmsg = nargoutchk(1,1,nargout);
if(~isempty(errmsg))
     disp(errmsg);
     return;
end
id = struct('SYM',0,'JOB',-1,'ICNTL',zeros(1,40)-9998,'CNTL',zeros(1,15)-9998,'PERM_IN',-9999,'COLSCA',-9999,'ROWSCA',-9999,'RHS',-9999,'INFOG',zeros(1,40)-9998,'RINFOG',zeros(1,40)-9998,'VAR_SCHUR',-9999,'SCHUR',-9999,'INST',-9999,'SOL',-9999,'REDRHS',-9999,'PIVNUL_LIST',-9999,'MAPPING',-9999,'SYM_PERM',-9999,'UNS_PERM',-9999,'TYPE',0);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function x = cvx_basis( x )
x = sparse( reshape( x, 1, numel( x ) ) );

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

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function saveSession(queryFlag)
%Save the mrSESSION file in HOMEDIR.
%
%  saveSession([queryFlag=0])
%
% 3/26/2001, djh, updated to version 3.0 to also save dataTYPES structure.
% 06/06, ras -- doesn't save any over other info (like vANATOMYPATH)
% that may already be present in the file.
% 10/06, ras -- auto saves a backup file, to help prevent saving over data.

global mrSESSION
global HOMEDIR
global dataTYPES
global vANATOMYPATH

if ~exist('queryFlag','var')
    queryFlag=0;
end

if isempty(HOMEDIR)
    warning(['Writing mrSESSION to: ',pwd]);
    pathStr = fullfile(pwd,'mrSESSION.mat');
else
    pathStr = fullfile(HOMEDIR,'mrSESSION.mat');
end

% save a backup if there's an existing file
if exist(pathStr, 'file')
    tmp = load(pathStr);
    save('mrSESSION_backup.mat', '-struct', 'tmp');
end

% The query flag says ask if the mrSESSION file already exists.
if exist(pathStr, 'file') && queryFlag
    but = questdlg('Over-write existing mrSESSION file?');
    switch but
    case 'Yes',
        save(pathStr,'mrSESSION','dataTYPES', 'vANATOMYPATH');
    otherwise
    end
else
    if exist(pathStr,'file')
        save(pathStr,'mrSESSION','dataTYPES', 'vANATOMYPATH', '-append');
    else
        save(pathStr,'mrSESSION','dataTYPES', 'vANATOMYPATH');
    end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function p = obj2pnts(obj)

% OBJ2PNTS Get points matrix from 3D object.
%   OBJ2PNTS(OBJ) returns the 3-by-N matrix of points in the structure OBJ,
%   stored in the N-by-3 matrix OBJ.vert.

%   Copyright 2009 Joan Sola @ LAAS-CNRS.

p = obj.vert';

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
% ?Rahul Kala, IIIT Allahabad, Creative Commons Attribution-ShareAlike 4.0 International License.
% The use of this code, its parts and all the materials in the text; creation of derivatives and their publication; and sharing the code publically is permitted without permission.
% Please cite the work in all materials as: R. Kala (2014) Code for Robot Path Planning using Bidirectional Rapidly-exploring Random Trees, Indian Institute of Information Technology Allahabad, Available at: http://rkala.in/codes.html

function feasible=feasiblePoint(point,map)
feasible=true;
% check if collission-free spot and inside maps
if ~(point(1)>=1 && point(1)<=size(map,1) && point(2)>=1 && point(2)<=size(map,2) && map(point(1),point(2))==1)
    feasible=false;
end

author: xuedidi


text:
function out = SD_SurrogateTest(x,surrMeth,numSurrs,extrap,theTestStat,randomSeed)
% SD_SurrogateTest   Analyzes test statistics obtained from surrogate time series
%
% This function is based on information found in:
% "Surrogate data test for nonlinearity including nonmonotonic transforms"
% D. Kugiumtzis Phys. Rev. E 62(1) R25 (2000)
%
% The generation of surrogates is done by the periphery function,
% SD_MakeSurrogates
%
%---INPUTS:
% x, the input time series
%
% surrMeth, the method for generating surrogate time series:
%       (i) 'RP': random phase surrogates that maintain linear correlations in
%                 the data but destroy any nonlinear structure through phase
%                 randomization
%       (ii) 'AAFT': the amplitude-adjusted Fourier transform method maintains
%                    linear correlations but destroys nonlinear structure
%                    through phase randomization, yet preserves the approximate
%                    amplitude distribution,
%       (iii) 'TFT': preserves low-frequency phases but randomizes high-frequency phases (as a way of dealing
%                    with non-stationarity, cf.:
%               "A new surrogate data method for nonstationary time series",
%                   D. L. Guarin Lopez et al., arXiv 1008.1804 (2010)
%
% numSurrs, the number of surrogates to compute (default is 99 for a 0.01
%         significance level 1-sided test)
%
% extrap, extra parameter, the cut-off frequency for 'TFT'
%
% theTestStat, the test statistic to evalute on all surrogates and the original
%           time series. Can specify multiple options in a cell and will return
%           output for each specified test statistic:
%           (i) 'ami': the automutual information at lag 1, cf.
%                 "Testing for nonlinearity in irregular fluctuations with
%                 long-term trends" T. Nakamura and M. Small and Y. Hirata,
%                 Phys. Rev. E 74(2) 026205 (2006)
%           (ii) 'fmmi': the first minimum of the automutual information
%                       function
%           (iii) 'o3': a third-order statistic used in: "Surrogate time
%                 series", T. Schreiber and A. Schmitz, Physica D 142(3-4) 346
%                 (2000)
%           (iv) 'tc3': a time-reversal asymmetry measure. Outputs of the
%                 function include a z-test between the two distributions, and
%                 some comparative rank-based statistics.
%
% randomSeed, whether (and how) to reset the random seed, using BF_ResetSeed

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function 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/>.
% ------------------------------------------------------------------------------

doPlot = 0; % plot outputs to a figure

% ------------------------------------------------------------------------------
%% Check inputs
% ------------------------------------------------------------------------------
if nargin < 2 || isempty(surrMeth)
    surrMeth = 'RP'; % randomize phases
end

if nargin < 3 || isempty(numSurrs)
    numSurrs = 99; % create 99 surrogates for a 0.01 significance level 1-sided test
end

if nargin < 4
    extrap = [];
end

if nargin < 5 || isempty(theTestStat)
    theTestStat = 'AMI'; % automutual information
end

if ischar(theTestStat)
    theTestStat = {theTestStat};
end

% randomSeed: how to treat the randomization
if nargin < 6
    randomSeed = [];
end

N = length(x); % time-series length

% ------------------------------------------------------------------------------
%% Generate surrogate time series using SD_MakeSurrogates:
% ------------------------------------------------------------------------------
z = SD_MakeSurrogates(x,surrMeth,numSurrs,extrap,randomSeed);

% z is matrix where each column is a surrogate time series

% ------------------------------------------------------------------------------
% Evaluate test statistic on each surrogate
% ------------------------------------------------------------------------------
if ismember('ami1',theTestStat)
    % Investigate AMI(1) of surrogates compared to that of signal itself
    % This statistic is used by Nakamura et al. (2006), PRE
    % Could use CO_HistogramAMI or TSTL, but I'll use IN_AutoMutualInfo with
    % Gaussian kernel…
    % FYI: for a histogram method, there are upper and lower bounds on the number of
    % bins to use: [1+log_2(N)], [sqrt(N)].

    % Use the gaussian approximation to estimate automutual information using the
    % Information Dynamics Toolkit:
    ami_fn = @(timeSeries,timeDelay) IN_AutoMutualInfo(timeSeries,timeDelay,'gaussian');

    AMIx = ami_fn(x,1);
    AMIsurr = zeros(numSurrs,1);
    for i = 1:numSurrs
        AMIsurr(i) = ami_fn(z(:,i),1);
    end
    % So we have a value AMIx, and a distribution for the surrogates
    % AMIsurr -- we must compare and return something meaningful
    % surrogates should always have lower AMI than original signal
    someStats = SDgivemestats(AMIx,AMIsurr,'right');
    fnames = fieldnames(someStats);
    for i = 1:length(fnames)
        out.(sprintf('ami_%s',fnames{i})) = someStats.(fnames{i});
    end
end

if ismember('fmmi',theTestStat)
    % Investigate the first minimum of mutual information of surrogates compared to
    % that of signal itself
    fmmix = CO_FirstMin(x,'mi');
    fmmiSurr = zeros(numSurrs,1);
    for i = 1:numSurrs
        try
            fmmiSurr(i) = CO_FirstMin(z(:,i),'mi');
        catch
            fmmiSurr(i) = NaN;
        end
    end
    if any(isnan(fmmiSurr))
        error('fmmi failed');
    end

    % FMMI should be higher for signal than surrogates
    someStats = SDgivemestats(fmmix,fmmiSurr,'right');
    fnames = fieldnames(someStats);
    for i = 1:length(fnames)
        out.(sprintf('fmmi_%s',fnames{i})) = someStats.(fnames{i});
    end
end

if ismember('o3',theTestStat)
    % Third-order statistic in Schreiber, Schmitz (Physica D)
    tau = 1;
    o3x = 1/(N-tau)*sum((x(1+tau:end) - x(1:end-tau)).^3);
    o3surr = zeros(numSurrs,1);
    for i = 1:numSurrs
        o3surr(i) = 1/(N-tau)*sum((z(1+tau:end,i) - z(1:end-tau,i)).^3);
    end
    someStats = SDgivemestats(o3x,o3surr,'both');
    fnames = fieldnames(someStats);
    for i = 1:length(fnames)
        out.(sprintf('o3_%s',fnames{i})) = someStats.(fnames{i});
    end
end

if ismember('tc3',theTestStat)
    % TC3 statistic -- another time-reversal asymmetry measure
    tau = 1;
    tmp = CO_tc3(x,tau);
    tc3x = tmp.raw;
    tc3surr = zeros(numSurrs,1);
    for i = 1:numSurrs
        tmp = CO_tc3(z(:,i),tau);
        tc3surr(i) = tmp.raw;
    end

    someStats = SDgivemestats(tc3x,tc3surr,'both');
    fnames = fieldnames(someStats);
    for i = 1:length(fnames)
        out.(sprintf('tc3_%s',fnames{i})) = someStats.(fnames{i});
    end
end

if ismember('nlpe',theTestStat)
    % Locally constant phase space prediction error
    warning('''nlpe'' can be very time consuming...')
    de = 3; tau = 1; % embedding parameters: fixed like a dummy!
    tmp = NL_MS_nlpe(x,de,tau);
    nlpex = tmp.msqerr;
    nlpesurr = zeros(numSurrs,1);
    for i = 1:numSurrs
        res = MS_nlpe(z(:,i),de,tau);
        msqerr = sum(res.^2);
        nlpesurr(i) = msqerr;
    end

    someStats = SDgivemestats(nlpex,nlpesurr,'right'); % NLPE should be higher than surrogates
    fnames = fieldnames(someStats);
    for i = 1:length(fnames)
        out.(sprintf('nlpe_%s',fnames{i})) = someStats.(fnames{i});
    end
end

if ismember('fnn',theTestStat)
    warning('fnn takes like *literally forever*...')

    % false nearest neighbours at d=2;
    tmp = NL_MS_fnn(x,2,1,5,1);
    fnnx = tmp.pfnn_2;
    fnnsurr = zeros(numSurrs,1);
    for i = 1:numSurrs
        tmp = NL_MS_fnn(z(:,i),2,1,5,1);
        fnnsurr(i) = tmp.pfnn_2;
    end

    someStats = SDgivemestats(fnnx,fnnsurr,'right'); % FNN(2) should be higher than surrogates?
    fnames = fieldnames(someStats);
    for i = 1:length(fnames)
        out.(sprintf('fnn_%s',fnames{i})) = someStats.(fnames{i});
    end
end

% ------------------------------------------------------------------------------
function someStats = SDgivemestats(statx,statsurr,leftrightboth)
    numSurrs = length(statsurr);
    if any(isnan(statsurr))
        error('SDgivemestats failed');
    end
%         leftrightboth = {'left','right','both'}
    % have a distribution of some statistic with samples statsurr
    % and we have the value of the statistic for a given process in
    % statx. Want to return measures of how consistant the measured
    % statistic is with the sample statsurr.

    % ASSUME GAUSSIAN DISTRIBUTION:
    % so can use 1/2-sided z-statistic
    [~, p, ~, zStat] = ztest(statx, mean(statsurr), std(statsurr), 0.05, leftrightboth);
    someStats.p = p; % pvalue
    someStats.zscore = zStat; % z-statistic

    % fit a kernel distribution to zscored distribution:
    if std(statsurr) == 0
        % all surrogates have same value of this statistic
        % cannot do a meaningful zscore -- do it raw
        [f, xi] = ksdensity(statsurr);
        % find where the statx value would be:
        if statx < min(xi) || statx > max(xi)
            someStats.f = 0; % out of range -- assume p=0 here
        else
            [~, minhere] = min(abs(statx-xi));
            someStats.f = f(minhere); % return probability density where the point is
        end
    else
        zscstatsurr = (statsurr-mean(statsurr))/std(statsurr);
        zscstatx = (statx-mean(statsurr))/std(statsurr);
        [f, xi] = ksdensity(zscstatsurr);

        % find where the statx value would be:
        if (zscstatx < min(xi)) || (zscstatx > max(xi))
            someStats.f = 0; % out of range -- assume p=0 here
        else
            [~, minhere] = min(abs(zscstatx-xi));
            someStats.f = f(minhere); % return probability density where the point is
        end
    end

    % What fraction of the range is the sample in?
    medsurr = median(statsurr);
    iqrsurr = iqr(statsurr);
    if iqrsurr == 0
        someStats.mediqr = NaN;
    else
        someStats.mediqr = abs(statx-medsurr)/iqrsurr;
    end

    % rank statistic
    [~, ix] = sort([statx; statsurr]);
    xfitshere = find(ix == 1) - 1;
    if strcmp(leftrightboth,'right') % x statistic smaller than distribution
        xfitshere = numSurrs + 1 - xfitshere; % how far from highest ranked value
    elseif strcmp(leftrightboth,'both')
        xfitshere = min(xfitshere,numSurrs+1-xfitshere);
    end

    if isempty(xfitshere)
        someStats.prank = 1/(numSurrs+1); % rank-based p-value
    else
        someStats.prank = (1+xfitshere)/(numSurrs+1); % rank-based p-value
    end
    if strcmp(leftrightboth,'both')
        someStats.prank = someStats.prank*2; % I think this factor should be in here
    end

    % DO PLOTTING:
    if doPlot
        figure('color','w')
        plot(statsurr,ones(numSurrs,1),'.k');
        hold('on');
        plot(statx*ones(2,1),[0,2],'r')
    end
end
% ------------------------------------------------------------------------------

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function rs = spec(s)

%tstoolbox/@signal/spec
%   Syntax:
%     * rs = spec(s)
%
%   compute power spectrum for real valued scalar signals. Multivariate
%   signals are accepted but may produce unwanted results as only the
%   spectrum of the first column is returned.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(1,1);

c = spec(s.core); 	% call real working routine for parent core object
rs = signal(c, s);	% special constructor calling syntax for working routines

a = getaxis(s, 1);
rs = setaxis(rs, 1, achse(unit(a)^(-1),0, samplerate(a)/dlens(s,1)));
rs = addhistory(rs, 'Calculated spectrum (spec)');
rs = setyunit(rs, yunit(s)^2);
rs = addcommandlines(rs, 's = spec(s');


text:
function [U, S] = pca(X)
%PCA Run principal component analysis on the dataset X
%   [U, S, X] = pca(X) computes eigenvectors of the covariance matrix of X
%   Returns the eigenvectors U, the eigenvalues (on diagonal) in S
%

% Useful values
[m, n] = size(X);

% You need to return the following variables correctly.
U = zeros(n);
S = zeros(n);

% ====================== YOUR CODE HERE ======================
% Instructions: You should first compute the covariance matrix. Then, you
%               should use the "svd" function to compute the eigenvectors
%               and eigenvalues of the covariance matrix.
%
% Note: When computing the covariance matrix, remember to divide by m (the
%       number of examples).
%

Sigma = 1.0/m .* X' * X;

[U, S, V] = svd(Sigma);

% =========================================================================

end

author: Borye

path: github-repos/MATLAB/Borye-machine-learning-coursera-1/machine-learning-coursera-1-033fdc2e6da393eeb1179a09aafe92362021effb/Week 8 Assignments/K-Means Clustering and PCA/mlclass-ex7/pca.m


text:
function [fevd_estimates]=olsfevd(irf_estimates,IRFperiods,gamma,n)

% function [fevd_estimates]=olsfevd(irf_estimates,IRFperiods,gamma,n,endo,datapath)
% computes and displays fevd values for the OLS VAR model
% inputs:  - cell 'irf_estimates': lower bound, point estimates, and upper bound for the IRFs
%          - integer 'IRFperiods': number of periods for IRFs
%          - matrix 'gamma': structural disturbance variance-covariance matrix (defined p 48 of technical guide)
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - cell 'endo': list of endogenous variables of the model
%          - string 'datapath': user-supplied path to excel data spreadsheet
% outputs: - cell 'fevd_estimates': lower bound, point estimates, and upper bound for the FEVD

% this function implements the procedure described p55-56

% preliminary tasks
% create the first cell
temp=cell(n,n+1);

% start by filling the first column of every Tij matrix in the cell
% loop over rows of temp
for jj=1:n
   % loop over columns of temp
   for ii=1:n
   % square each element
   temp{jj,ii}(:,1)=irf_estimates{jj,ii}(:,1).^2;
   end
end
% fill all the other entries of the Tij matrices
% loop over rows of temp
for jj=1:n
   % loop over columns of temp
   for ii=1:n
      % loop over remaining columns
      for kk=2:IRFperiods
      % define the column as the square of the corresponding column in orthogonalised_irf_record
      % additioned to the value of the preceeding columns, which creates the cumulation
      temp{jj,ii}(:,kk)=irf_estimates{jj,ii}(:,kk).^2+temp{jj,ii}(:,kk-1);
      end
   end
end
% multiply each matrix in the cell by the variance of the structural shocks
% loop over rows of temp
for ii=1:n
% loop over columns of temp
   for jj=1:n
   % multiply column jj of the matrix by the variance of the structural shock
   temp{ii,jj}(1,:)=temp{ii,jj}(1,:)*gamma(jj,jj);
   end
end

% obtain now the values for Ti, the (n+1)th matrix of each row
% loop over rows of temp
for ii=1:n
% start the summation over Tij matrices
temp{ii,n+1}=temp{ii,1};
   % sum over remaining columns
   for jj=2:n
   temp{ii,n+1}=temp{ii,n+1}+temp{ii,jj};
   end
end

% create the output cell fevd_record
fevd_estimates=cell(n,n);
% fill the cell
% loop over rows of fevd_estimates
for ii=1:n
   % loop over columns of fevd_estimates
   for jj=1:n
   % define the matrix Vfij as the division (pairwise entry) of Tfij by Tfj
   fevd_estimates{ii,jj}=temp{ii,jj}./temp{ii,n+1};
   end
end

author: european-central-bank

repo: BEAR-toolbox


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

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

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

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d  ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );
    fprintf ( 1, '\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d: ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%7d  ', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ p, seed ] = triangle_reference_sample ( n, seed )

%*****************************************************************************80
%
%% TRIANGLE_REFERENCE_SAMPLE returns random points in the reference triangle.
%
%  Diagram:
%
%       3
%    s  |\
%    i  | \
%    d  |  \
%    e  |   \  side 2
%       |    \
%    3  |     \
%       |      \
%       1-------2
%
%         side 1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points to generate.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real P(2,N), random points in the triangle.
%
%    Output, integer SEED, a seed for the random number generator.
%
  dim_num = 2;

  for j = 1 : n

    [ r, seed ] = r8_uniform_01 ( seed );
%
%  Interpret R as a percentage of the triangle's area.
%
%  Imagine a line L, parallel to side 1, so that the area between
%  vertex 1 and line L is R percent of the full triangle's area.
%
%  The line L will intersect sides 2 and 3 at a fraction
%  ALPHA = SQRT ( R ) of the distance from vertex 1 to vertices 2 and 3.
%
    alpha = sqrt ( r );
%
%  Now choose, uniformly at random, a point on the line L.
%
    [ beta, seed ] = r8_uniform_01 ( seed );

    p(1,j) = ( 1.0 - beta ) * alpha;
    p(2,j) =         beta   * alpha;

  end

  return
end


text:
function u = bsv_upwind ( a, b, alpha, beta, nu, n, output )

%*****************************************************************************80
%
%% BSV_UPWIND applies Newton's method to a discretized steady viscous Burgers equation.
%
%  Discussion:
%
%    An upwind scheme is used to combat the numerical oscillations that
%    can occur otherwise.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the left and right endpoints.
%
%    Input, real ALPHA, BETA, the Dirichlet boundary values at the left
%    and right.
%
%    Input, real NU, the viscosity.  Normally, 0 < NU.
%
%    Input, integer N, the number of nodes to use between A and B.
%
%    Input, logical OUTPUT, is TRUE if printout is desired.
%
%    Output, real U(N), the computed discretized solution.
%
  if ( n < 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BSV_UPWIND - Fatal error!\n' );
    fprintf ( 1, '  N < 2.\n' );
    error ( 'BSV_UPWIND - Fatal error!\n' );
  end
%
%  Set some iteration parameters.
%
  newton_step_max = 50;
  newton_resid_tol = sqrt ( eps );
  newton_step_tol = sqrt ( eps );
%
%  Use equally spaced nodes from A to B.
%
  dx = ( b - a ) / ( n - 1 );
%
%  The initial guess will be the linear interpolant to the boundary conditions.
%
  u = ( linspace ( alpha, beta, n ) )';
%
%  Prepare for the Newton iteration.
%
  J = sparse ( [], [], [], n, n, 3 * n );
  newton_step = 0;

  if ( output )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Step  ||F(U)||\n' );
    fprintf ( 1, '\n' );
  end

  while ( newton_step <= newton_step_max )

    f(1,1) = u(1) - alpha;

    f(2:n-1,1) = - nu * ( u(3:n) - 2.0 * u(2:n-1) + u(1:n-2) ) / ( dx^2 );
%
%  Massaging the results of FIND can be tedious.
%
    i = find ( 0.0 <= u(2:n-1) );
    i = i + 1;
    j = find ( u(2:n-1) < 0.0 );
    j = j + 1;

    f(i,1) = f(i,1) + 0.5 * ( u(i).^2 - u(i-1).^2 ) / dx;
    f(j,1) = f(j,1) + 0.5 * ( u(j+1).^2 - u(j).^2 ) / dx;

    f(n,1) = u(n) - beta;

    f_norm = norm ( f, inf );

    if ( output )
      fprintf ( 1, '  %4d  %g\n', newton_step, f_norm );
    end

    if ( f_norm < newton_resid_tol )
      break
    end
%
%  Define the Jacobian matrix.
%
    J(1,1) = 1.0;
    for i = 2 : n - 1
      if ( 0.0 <= u(i) )
        J(i,i-1) = - u(i-1) / dx        - nu / dx^2;
        J(i,i)   =   u(i)   / dx  + 2.0 * nu / dx^2;
        J(i,i+1) =                      - nu / dx^2;
      else
        J(i,i-1) =                      - nu / dx^2;
        J(i,i)   = - u(i)   / dx  + 2.0 * nu / dx^2;
        J(i,i+1) =   u(i+1) / dx        - nu / dx^2;
      end
    end
    J(n,n) = 1.0;
%
%  Solve the linear system J * du = -f to get the Newton update.
%
    du = J \ f;
    du_norm = norm ( du, inf );

    u_norm = norm ( u, inf );
    if ( du_norm < newton_step_tol * ( u_norm + 1.0 ) )
      break
    end

    u = u - du;
    newton_step = newton_step + 1;

  end

  if ( newton_step_max < newton_step )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BSV_UPWIND - Warning!\n' );
    fprintf ( 1, '  The Newton iteration did not converge.\n' );
  end

  return
end


text:
function line = medianLine(varargin)
%MEDIANLINE Create a median line between two points
%
%   L = medianLine(P1, P2);
%   Create the median line of points P1 and P2, that is the line containing
%   all points located at equal distance of P1 and P2.
%
%   L = medianLine(PTS);
%   Creates the median line of 2 points, given as a 2*2 array. Array has
%   the form:
%   [ [ x1 y1 ] ; [ x2 y2 ] ]
%
%   L = medianLine(EDGE);
%   Creates the median of the edge. Edge is a 1*4 array, containing [X1 Y1]
%   coordinates of first point, and [X2 Y2], the coordinates of the second
%   point.
%
%   Example
%   % Draw the median line of two points
%     P1 = [10 20];
%     P2 = [30 50];
%     med = medianLine(P1, P2);
%     figure; axis square; axis([0 100 0 100]);
%     drawEdge([P1 P2], 'linewidth', 2, 'color', 'k');
%     drawLine(med)
%
%   % Draw the median line of an edge
%     P1 = [50 60];
%     P2 = [80 30];
%     edge = createEdge(P1, P2);
%     figure; axis square; axis([0 100 0 100]);
%     drawEdge(edge, 'linewidth', 2)
%     med = medianLine(edge);
%     drawLine(med)
%
%
%   See also:
%   lines2d, createLine, orthogonalLine
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 31/10/2003.
%

% history
% 2010-08-06 vectorize and change behaviour for N-by-4 inputs

nargs = length(varargin);

if nargs == 1
    tab = varargin{1};
    if size(tab, 2)==2
        % input is an array of two points
        x0 = tab(1,1);
        y0 = tab(1,2);
        dx = tab(2,1)-x0;
        dy = tab(2,2)-y0;
    else
        % input is an edge
        x0 = tab(:, 1);
        y0 = tab(:, 2);
        dx = tab(:, 3) - tab(:, 1);
        dy = tab(:, 4) - tab(:, 2);
    end

elseif nargs==2
    % input is given as two points, or two point arrays
    p1 = varargin{1};
    p2 = varargin{2};
    x0 = p1(:, 1);
    y0 = p1(:, 2);
    dx = bsxfun(@minus, p2(:, 1), x0);
    dy = bsxfun(@minus, p2(:, 2), y0);

else
    error('Too many input arguments');
end

% compute median using middle point of the edge, and the direction vector
% rotated by 90 degrees counter-clockwise
line = [bsxfun(@plus, x0, dx/2), bsxfun(@plus, y0, dy/2), -dy, dx];

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function r = foldl(list, identity, func)
% VAL = FOLDL(ARRAY, ID, FUNC)
%   Perform a left fold over ARRAY with FUNC using ID as the identity.
%   FOLDL applies FUNC to the current identity (initially ID) and each item
%   of ARRAY successively.  Thus, FOLDL iterates over ARRAY with FUNC.  ID
%   is updated to be the result of this function.  The final result of ID
%   is returned from the function.
%
%   FOLDL is most simply viewed as the haskell code
%     foldl []   id f = id
%     foldl x:xs id f = foldl xs (f id x) f
%
%   FOLDL and FOLDR are implemented as a loop for optimization purposes in
%   these libraries.  FOLDL is tail recursive and thus more efficient on
%   most platforms that use the above functional definition.
%
%   Examples of Common FOLDL uses:
%     sum(array) == foldl(array, 0, @(a,b) a + b)
%     all(array) == foldl(array, true, @(a,b) a && b)
%     any(array) == foldl(array, false, @(a,b) a || b)
%

  if isempty(list)
    r = identity;
    return;
  end

  if iscell(list)
    for i = 1:numel(list)
      identity = func(identity, list{i});
    end
  else
    for i = 1:numel(list)
      identity = func(identity, list(i));
    end
  end

  r = identity;
end


text:
function r = erf(a)
%ERF          Hessian (elementwise) error function
%

% written  31/05/13     S.M. Rump
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  a.x = full(a.x);
  K = prod(size(a.x));
  % factorLB <= 2/sqrt(pi) <= factorUB
  INTLAB_STDFCTS_ERF = getappdata(0,'INTLAB_STDFCTS_ERF');
  factorLB = INTLAB_STDFCTS_ERF.TWO_SQRTPIINF;  % round to nearest
  factorUB = INTLAB_STDFCTS_ERF.TWO_SQRTPISUP;  % ~ 1.12

  if K==1                   % scalar hessian

    ax = exp(-a.x(:).^2);
    if isa(a.x,'intval')
        ax = intval(factorLB,factorUB,'infsup') * ax;
    else
      ax = factorLB * ax;
    end
    r.x = erf(a.x);
    r.dx = ax * a.dx;
    r.hx = ax * ( a.hx - reshape( (a.x.*a.dx) * a.dx.' , size(a.hx) ) );

  else                      % matrix hessian

    N = getappdata(0,'INTLAB_HESSIAN_NUMVAR');
    N2 = N^2;

    r.x = erf(a.x);
    a.x = a.x(:);
    if issparse(a.hx)               % input sparse

      ax = exp(-a.x.^2);
      if isa(a.x,'intval')
        ax = intval(factorLB,factorUB,'infsup') * ax;
      else
        ax = factorLB * ax;
      end
      sizeax = length(ax);
      [ia,ja,sa] = find(a.dx);
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if isempty(ia)
        r.dx = sparse([],[],[],N,sizeax);
        r.hx = a.hx;
      else
        if isa(a.x,'intval')          % sparse intval
          rdx = times(ax(ja),sa(:),0);
          adx1 = times(a.x(ja),sa(:),0);
          r.dx = intval( sparse(ia,ja,rdx.inf,N,sizeax) , sparse(ia,ja,rdx.sup,N,sizeax) , 'infsup' );
          rdx = intval(sparse(ia,ja,adx1.inf,N,sizeax),sparse(ia,ja,adx1.sup,N,sizeax),'infsup');
        else                          % sparse point
          r.dx = sparse(ia,ja,ax(ja).*sa(:),N,sizeax);
          rdx = sparse(ia,ja,a.x(ja).*sa(:),N,sizeax);
        end
        r.hx = a.hx - adx2rhx(N,sizeax,rdx,a.dx);
      end
      [ia,ja,sa] = find(r.hx);        % sparse point or intval
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if ~isempty(ia)
        if isa(ax,'intval')
          rhx = times(ax(ja),sa(:),0);
          r.hx = sparse(ia,ja,intval(rhx.inf,rhx.sup,'infsup'),N2,sizeax);
        else
          r.hx = sparse(ia,ja,ax(ja).*sa(:),N2,sizeax);
        end
      end

    else                            % input full

      if isa(a.x,'intval')
        ax = intval(factorLB,factorUB,'infsup') * exp((-((a.x).').^2));
      else
        ax = factorLB * exp((-((a.x).').^2));
      end
      ax = ax(ones(N*N,1),:);
      r.dx = a.dx .* ax(1:N,:);
      adx = repmat(a.x(:).',N,1) .* a.dx;
      r.hx = ( a.hx - adx(repmat(1:N,N,1),:) .* a.dx(repmat(1:N,1,N),:) ) .* ax;

    end

  end

  r = class(r,'hessian');

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function [Population,FrontNo,CrowdDis] = LCSA_NSGAIIEnvironmentalSelection(Population,N)
% -----------------------------------------------------------------------
%  Copyright (C) 2020 Heiner Zille
%
%  This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0
%  International License. (CC BY-NC-SA 4.0). To view a copy of this license,
%  visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or see the
%  pdf-file "License-CC-BY-NC-SA-4.0.pdf" that came with this code.
%
%  You are free to:
%  * Share ? copy and redistribute the material in any medium or format
%  * Adapt ? remix, transform, and build upon the material
%  Under the following terms:
%  * Attribution ? You must give appropriate credit, provide a link to the
%     license, and indicate if changes were made. You may do so in any reasonable
%     manner, but not in any way that suggests the licensor endorses you or your use.
%  * NonCommercial ? You may not use the material for commercial purposes.
%  * ShareAlike ? If you remix, transform, or build upon the material, you must
%    distribute your contributions under the same license as the original.
%  * No additional restrictions ? You may not apply legal terms or technological
%    measures that legally restrict others from doing anything the license permits.
%
%  Author of this Code:
%   Heiner Zille <heiner.zille@ovgu.de> or <heiner.zille@gmail.com>
%
%  This code is based on the following publications:
%
%  1) Heiner Zille
%     "Large-scale Multi-objective Optimisation: New Approaches and a Classification of the State-of-the-Art"
%     PhD Thesis, Otto von Guericke University Magdeburg, 2019
%     http://dx.doi.org/10.25673/32063
%
%  2) Heiner Zille and Sanaz Mostaghim
%     "Linear Search Mechanism for Multi- and Many-Objective Optimisation"
%     10th International Conference on Evolutionary Multi-Criterion Optimization (EMO 2019),
%        Lecture Notes in Computer Science, vol 11411.
%        Deb K. et al. (eds), Springer, Cham, East Lansing, Michigan, USA, March 2019
%     https://doi.org/10.1007/978-3-030-12598-1_32.
%
%  This file is intended to work with the PlatEMO framework version 2.5.
%  Date of publication of this code: 06.04.2020
%  Last Update of this code: 06.04.2020
%  A newer version of this algorithm may be available. Please contact the author
%  or see http://www.ci.ovgu.de/Research/Codes.html.
%
% The files may have been modified in Feb 2021 by the authors of the Platemo framework to work with the Platemo 3.0 release.
% -----------------------------------------------------------------------
% This file is derived from its original version containied in the PlatEMO
% framework.
% -----------------------------------------------------------------------

    %% Non-dominated sorting
    [FrontNo,MaxFNo] = NDSort(Population.objs,Population.cons,N);
    Next = FrontNo < MaxFNo;

    %% Calculate the crowding distance of each solution
    CrowdDis = CrowdingDistance(Population.objs,FrontNo);

    %% Select the solutions in the last front based on their crowding distances
    Last     = find(FrontNo==MaxFNo);
    [~,Rank] = sort(CrowdDis(Last),'descend');
    Next(Last(Rank(1:N-sum(Next)))) = true;

    %% Population for next generation
    Population = Population(Next);
    FrontNo    = FrontNo(Next);
    CrowdDis   = CrowdDis(Next);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
classdef TirePolynomial < VehicleDynamicsLateral.Tire
    % TirePolynomial Polynomial tire model
    %
    % It inherits methods from Tire.

    methods
        % Constructor
        function self = TirePolynomial()
            self.k1 = 115000;
            self.k2 = 560000;
        end

        function p = PlotTire(self)
            % Returns the handle of the curve
            alpha = (0:0.1:15)*pi/180;
            Fy = - self.Characteristic(alpha);
            p = plot(alpha*180/pi,Fy);
            grid on; box on;
            xlabel('Slip angle [deg]')
            ylabel('Lateral force [N]')
        end

        function Fy = Characteristic(self, alpha, varargin)
            % Lateral force
            Fy = - (self.k1 * alpha - self.k2 * alpha.^3);
        end
    end

    properties
        k1 % 1st polynomial coefficient, cornering stiffness [N/rad]
        k2 % 2nd polynomial coefficient [N/rad^3]
    end
end

%% See Also
%
% <../../index.html Home>
%

author: andresmendes

repo: Vehicle-Dynamics-Lateral


text:
function a = r8col_sortr_a ( m, n, a, key )

%*****************************************************************************80
%
%% R8COL_SORTR_A ascending sorts one column of an R8COL, adjusting all entries.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), an unsorted M by N array.
%
%    Input, integer KEY, the column in which the "key" value
%    is stored.  On output, column KEY of the array will be
%    in nondecreasing order.
%
%    Output, real A(M,N), rows of the array have been shifted in such
%    a way that column KEY of the array is in nondecreasing order.
%
  if ( m <= 0 )
    return
  end

  if ( key < 1 || n < key )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8COL_SORTR_A - Fatal error!\n' );
    fprintf ( 1, '  The value of KEY is not a legal column index.\n' );
    fprintf ( 1, '  KEY = %d\n', key );
    fprintf ( 1, '  N = %d\n', n );
    error ( 'R8COL_SORTR_A - Fatal error!' );
  end
%
%  Initialize.
%
  i = 0;
  indx = 0;
  isgn = 0;
  j = 0;
%
%  Call the external heap sorter.
%
  while ( 1 )

    [ indx, i, j ] = sort_heap_external ( m, indx, isgn );
%
%  Interchange the I and J objects.
%
    if ( 0 < indx )

      a = r8row_swap ( m, n, a, i, j );
%
%  Compare the I and J objects.
%
    elseif ( indx < 0 )

      if ( a(i,key) < a(j,key) )
        isgn = -1;
      else
        isgn = +1;
      end

    elseif ( indx == 0 )

      break

    end

  end

  return
end


text:
% -------------------------------------------------------------------------------------------------
function [imgs, pos, target_sz] = load_video_info(base_path, video)
%LOAD_VOT_VIDEO_INFO
%   Loads all the relevant information for the video in the given path:
%   the list of image files (cell array of strings), initial position
%   (1x2), target size (1x2), the ground truth information for precision
%   calculations (Nx4, for N frames), and the path where the images are
%   located. The ordering of coordinates and sizes is always [y, x].
%
%   Joao F. Henriques, 2014
%   http://www.isr.uc.pt/~henriques/
% -------------------------------------------------------------------------------------------------
	%full path to the video's files
	if base_path(end) ~= '/' && base_path(end) ~= '\',
		base_path(end+1) = '/';
	end
	video_path = [base_path video '/img/'];

	%load ground truth from text file
	ground_truth = csvread([base_path '/' video '/' 'groundtruth_rect.txt']);
	region = ground_truth(1, :);
	[cx, cy, w, h] = get_axis_aligned_BB(region);
    pos = [cy cx]; % centre of the bounding box
    target_sz = [h w];

	%load all jpg files in the folder
	img_files = dir([video_path '*.jpg']);
	assert(~isempty(img_files), 'No image files to load.')
	img_files = sort({img_files.name});

	%eliminate frame 0 if it exists, since frames should only start at 1
	img_files(strcmp('00000000.jpg', img_files)) = [];
    img_files = strcat(video_path, img_files);
    % read all frames at once
    imgs = vl_imreadjpeg(img_files,'numThreads', 12);
end

author: shenjianbing

repo: TripletTracking


text:
% Gps
p_drone.sd_gps_pn = 0.21; % standard deviation (sd)
p_drone.sd_gps_pe = 0.21;
p_drone.sd_gps_ph = 0.4;

p_drone.sd_gps_vn = 0.05; % standard deviation (sd)
p_drone.sd_gps_ve = 0.05;
p_drone.sd_gps_vh = 0.01;
p_drone.Ts_gps = 1; % [s]
p_drone.k_gps = 1/1100;

% Gyroscope
p_drone.sd_gyro = deg2rad(0.13); % standard deviation (sd) [deg/s]

% Acceleration
p_drone.sd_accel = 0.0025*p_physics.gravity; % standard deviation (sd) [m/s^2]

% Pressure
p_drone.pres0 = 101325; % static pressure at sea level    [Pa]
p_drone.sd_static_pres = 10; % standard deviation (sd)    [Pa]
p_drone.sd_diff_pres = 2; % standard deviation (sd)       [Pa]

p_drone.bias_gyro_x = 0.01;
p_drone.bias_gyro_y = 0.01;
p_drone.bias_gyro_z = 0.01;

p_drone.a_gyro            = 0.1;
p_drone.a_static_pres     = 0.1;
p_drone.a_diff_pres       = 0.5;
p_drone.a_accel           = 0.1;
p_drone.a_gps_speed       = 0;
p_drone.a_gps_pos         = 0;

p_drone.Q = [ p_drone.sd_gyro^2, 0;...
        0,           p_drone.sd_gyro^2];

p_drone.R = ([ p_drone.sd_accel, 0,          0;...
        0,           p_drone.sd_accel, 0;...
        0,           0,          p_drone.sd_accel]).^2;

p_drone.b_chidot  = 1.2;
p_drone.b_chi     = 1;
p_drone.b_hdot    = 1.7;
p_drone.b_h       = 1.7;
p_drone.b_Va      = 8;
p_drone.gamma_max = 10;

p_drone.bias_gyro_x = 0;
p_drone.bias_gyro_y = 0;
p_drone.bias_gyro_z = 0;

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
% Pick an index from the vector v,
% where each element is 0 or 1.
% Only pick from the elements denoted by 1.
%
% Output
%  item: scalar indx
function item = randset(v)
    assert(isvector(v));
    assert(all(v==0|v==1));
    indx = find(v);
    z = randint(1,1,[1 length(indx)]);
    item = indx(z);
end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
%% set up path

[~,computername] = system('hostname');
computername = strtrim(computername);

switch computername,

  case 'bransonk-desktop',

    JCtrax_path = '/groups/branson/home/bransonk/tracking/code/JCtrax';
    configfilename = 'params/JLabelParams_bransonk-desktop.xml';
    rootdatadir = '/groups/branson/home/bransonk/behavioranalysis/data/larvae_vani';

  case 'bransonk-lw2',

    JCtrax_path = 'C:\Code\JCtrax';
    configfilename = 'params\JLabelParamsKristinChase.xml';
    rootdatadir = 'C:\Data\larvae_vani';

  otherwise

    warning('Unknown computer %s. Paths may not be setup correctly',computername);
    JCtrax_path = '/groups/branson/home/bransonk/tracking/code/JCtrax';
    configfilename = 'params/JLabelParams_bransonk-desktop.xml';
    rootdatadir = '/groups/branson/home/bransonk/behavioranalysis/data/larvae_vani';

end

addpath(genpath(fullfile(JCtrax_path,'pdollar_toolbox')));
addpath(fullfile(JCtrax_path,'misc'));
addpath(fullfile(JCtrax_path,'filehandling'));
jlabelpath = fileparts(which('JLabel'));
addpath(fullfile(jlabelpath,'compute_perframe_features'));

%% parameters

genotype_name = 'genotypetest';
experiment_name = '20111213-144536';

%%

expdir = fullfile(rootdatadir,genotype_name,experiment_name);

ConvertLarvaeVaniTrxToTrx(expdir);

author: kristinbranson

repo: JAABA


text:
function H = per(sequence,isplot)
%
% 'per' estimate the hurst parameter of a given sequence with periodogram
%     method.
%
% Inputs:
%     sequence: the input sequence for estimate
%     isplot: whether display the plot. without a plot if isplot equal to 0
% Outputs:
%     H: the estimated hurst coeffeient of the input sequence

%  Author: Chu Chen
%  Version 1.0,  03/10/2008
%  chen-chu@163.com
%

if nargin == 1
    isplot = 0;
end

n = length(sequence);
Xk = fft(sequence);
P_origin = abs(Xk).^2/(2*pi*n);
P = P_origin(1:floor(n/2)+1);

x = log10((pi/n)*[2:floor(0.5*n)]);
y = log10(P(2:floor(0.5*n)));

% Use the lowest 20% part of periodogram to estimate the similarity.
X = x(1:floor(length(x)/5));
Y = y(1:floor(length(y)/5));
p1 = polyfit(X,Y,1);
Yfit = polyval(p1,X);
H = (1-(Yfit(end)-Yfit(1))/(X(end)-X(1)))/2;

if isplot ~= 0
    figure,clf,hold on;
    plot(x,y,'b.');
    plot(X,Yfit,'r-','LineWidth',3);
    xlabel('log10(Frequency)'),ylabel('log10(Periodogram)'),title('Periodogram Method');
end

author: Sable


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

% Append Current Temporary Matrix
function MU_uploadTmp(h)
handles = guidata(h.MU_matrix_display);
MatrixName=get(h.Matrix_name_edit,'String');
try
    if ~isempty(h.V.ROI)
        switch h.V.ROI.ROI_flag
            case {7,8,9}
                Flag = MU_load_matrix([MatrixName '_tmp'], h.Mask, 0);
            otherwise
                TMatrix=MU_Crop(handles);
                Flag = MU_load_matrix([MatrixName '_tmp'], TMatrix, 0);
        end
    else
        TMatrix=MU_Crop(handles);
        Flag = MU_load_matrix([MatrixName '_tmp'], TMatrix, 0);
    end
catch me
    error_msg{1,1}='ERROR!!! Saving current matrix aborted.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);
    return;
end

if ~Flag
    errordlg('Uploading temporary matrix failed!');
end

end

function TMatrix=MU_Crop(handles)

if numel(handles.V.DimSize)>=3
    dim = inputdlg(['Please specify up to which dimension to save (2 ~ ' num2str(numel(handles.V.DimSize)) ').'],'Specify Dimension',1,{num2str(numel(handles.V.DimSize))});
    if isempty(dim)
        error('Saving matrix was cancelled.');
    end
    try
        if str2double(dim{1})<2 | str2double(dim{1})>numel(handles.V.DimSize)
            error('The input dimension is invalid.');
        end
        if str2double(dim{1}) == 2
            TMatrix=handles.TMatrix(:,:,handles.V.Slice);
        else
            if str2double(dim{1})+1>numel(handles.V.DimSize)
                dimFlag =[];
            else
                dimFlag=num2str(handles.V.DimPointer(str2double(dim{1})+1:end),',%u');
            end
            eval(['TMatrix=handles.TMatrix(:,:' repmat(',:', [1 str2double(dim{1})-2]) dimFlag ');']);
        end
    catch me
        error('The input dimension is invalid.');
    end
else
    TMatrix=handles.TMatrix(:,:);
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function hostip = hlp_hostip()
% Find the host IP address in a platform-independent fashion.

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

import java.net.InetAddress;
hostip = char(InetAddress.getLocalHost());
hostip = hostip(find(hostip=='/',1)+1:end);

repo: REST


text:
%
% A variational approach to SPCP (Aravkin et al. 2014)
%
% RPCA | Lag-SPCP-SPG | Lagrangian SPCP solved by Spectral Projected Gradient (Aravkin et al. 2014)
% process_video('RPCA', 'Lag-SPCP-SPG', 'dataset/demo.avi', 'output/demo_Lag-SPCP-SPG.avi');

alg_path_aux = fullfile(lrs_conf.rpca_path,'SPGL1');
addpath(genpath(alg_path_aux));

nFrames     = size(M,2);
lambda      = 1/sqrt(max(size(M,1),size(M,2)));
L0          = repmat(median(M,2), 1, nFrames);
S0          = M - L0;
epsilon     = 5e-3*norm(M,'fro'); % tolerance for fidelity to data

% Lagrangian SPCP solved by Spectral Projected Gradient
opts    = struct('sum',false,'L0',L0,'S0',S0,'max',false,'tol',1e-3);
lambdaL = 0.25; lambdaS = 0.01; % (Aravkin et al. 2014)
[L,S] = solver_RPCA_Lagrangian(M,lambdaL,lambdaS,[],opts);

rmpath(genpath(alg_path_aux));

author: andrewssobral

repo: lrslibrary


text:
function a = setname(a, newname)

%tstoolbox/@achse/setname
%   Syntax:
%     * a = setname(a, newname)
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

a.name = newname;


text:
% DEMO_MATLAB: Examples of MATLAB-Python integration.
%
% FUNCTIONS:
%    - demo_matlab.py : Original Python example
%    - demo_matlab.m  : Call demo_matlab.py from Matlab
%    - demo_matlab.m  : Reproduce demo_matlab.py with Matlab calls
%
% REFERENCES
%    - Installation: https://neuroimage.usc.edu/brainstorm/MnePython
%    - Dataset: https://neuroimage.usc.edu/brainstorm/Tutorials/MedianNerveCtf
%
% AUTHORS:
%    - Python version: Mainak Jas <mainak.jas@telecom-paristech.fr>
%    - Matlab version: Francois Tadel, 2018-2020

% Initialize Brainstorm-Python
bst_python_init('Initialize', 1);

% Get test dataset
try
    data_path = py.mne.datasets.brainstorm.bst_raw.data_path();
catch
    disp([10 'MATLAB does not support the use of Python function input().' 10 ...
        'You must download the dataset manually from a Python terminal:' 10 ...
        '>>> import mne' 10 ...
        '>>> mne.datasets.brainstorm.bst_raw.data_path()' 10 10]);
end

% Read CTF dataset
raw_fname = fullfile(char(data_path), 'MEG', 'bst_raw', 'subj001_somatosensory_20111109_01_AUX-f.ds');
raw = py.mne.io.read_raw_ctf(raw_fname, pyargs('preload', true));
% Plot raw data
pyFig = raw.plot();
% Change the color of the first axes of the figure
pyAxes = pyFig.get_axes();
pyAxes{1}.set_facecolor({0,1,0});

% Set EOG channel
raw.set_channel_types(py.dict(pyargs('EEG058', 'eog')));
% Alternative syntax using the automatic conversion struct<=>dict
raw.set_channel_types(struct('EEG058', 'eog'));

% Show power line interference
raw.plot_psd();
% Notch filter: 60Hz
% raw.notch_filter(py.numpy.arange(60, 181, 60));

% Read events
events = py.mne.find_events(raw, pyargs('stim_channel', 'UPPT001'));
event_id = uint8(2);    % Right-hand somato

% Pick MEG channels
picks = py.mne.pick_types(raw.info, pyargs('meg', true, 'eeg', false, 'stim', false, 'eog', true, 'exclude', 'bads'));

% Compute epochs: [-100, +300]ms
epochs = py.mne.Epochs(raw, events, event_id, -0.1, 0.3, pyargs(...
    'picks',     picks, ...
    'baseline',  py.tuple({py.None, 0}), ...
    'reject',    py.dict(pyargs('mag', 4e-12, 'eog', 250e-6)), ...
    'preload',   false));

% Compute evoked
evoked = epochs.average();

% Fix stim artifact
py.mne.preprocessing.fix_stim_artifact(evoked);

% Correct delays due to hardware (stim artifact is at 4 ms)
evoked.shift_time(-0.004);

% Plot the result
evoked.plot();
evoked.plot_topomap(pyargs('times', py.numpy.array({0.016, 0.030, 0.060, 0.070})));

author: brainstorm-tools


text:
function moment = truncated_normal_b_moment ( order, mu, sigma, b )

%*****************************************************************************80
%
%% TRUNCATED_NORMAL_B_MOMENT returns a moment of the upper truncated Normal PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Phoebus Dhrymes,
%    Moments of Truncated Normal Distributions,
%    May 2005.
%
%  Parameters:
%
%    Input, integer ORDER, the order of the moment.
%    0 <= ORDER.
%
%    Input, real MU, SIGMA, the mean and standard deviation of the
%    parent Normal distribution.
%
%    Input, real B, the upper truncation limit.
%
%    Output, real MOMENT, the moment of the PDF.
%
  if ( order < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRUNCATED_NORMAL_B_MOMENT - Fatal error!\n' );
    fprintf ( 1, '  ORDER < 0.\n' );
    error ( 'TRUNCATED_NORMAL_B_MOMENT - Fatal error!\n' );
  end

  h = ( b - mu ) / sigma;
  pdf = normal_01_pdf ( h );
  cdf = normal_01_cdf ( h );

  if ( cdf == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRUNCATED_NORMAL_B_MOMENT - Fatal error!\n' );
    fprintf ( 1, '  PDF/CDF ratio fails, because CDF too small.\n', cdf );
    fprintf ( 1, '  PDF = %g\n', pdf );
    fprintf ( 1, '  CDF = %g\n', cdf );
    error ( 'TRUNCATED_NORMAL_B_MOMENT - Fatal error!\n' );
  end

  f = pdf / cdf;

  moment = 0.0;
  irm2 = 0.0;
  irm1 = 0.0;

  for r = 0 : order

    if ( r == 0 )
      ir = 1.0;
    elseif ( r == 1 )
      ir = - f;
    else
      ir = - h ^ ( r - 1 ) * f + ( r - 1 ) * irm2;
    end

    moment = moment + r8_choose ( order, r ) ...
      * mu ^ ( order - r ) ...
      * sigma ^ r * ir;

    irm2 = irm1;
    irm1 = ir;

  end

  return
end


text:
function stroud_test045 ( )

%*****************************************************************************80
%
%% TEST045 tests BALL_UNIT_VOLUME_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST045\n' );
  fprintf ( 1, '  In 3 dimensions:\n' );
  fprintf ( 1, '  BALL_UNIT_VOLUME_3D gets the volume of the unit ball.\n' );
  fprintf ( 1, '  BALL_UNIT_VOLUME_ND will be called for comparison.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  N    Volume    Method\n' );
  fprintf ( 1, '\n' );

  fprintf ( 1, '  %1d  %12f  %s\n', n, ball_unit_volume_3d ( ), ...
    'BALL_UNIT_VOLUME_3D' );

  fprintf ( 1, '  %1d  %12f  %s\n', n, ball_unit_volume_nd ( n ), ...
    'BALL_UNIT_VOLUME_ND' );

  return
end


text:
function [ x,y ] = scd_scheme_deltaDeltaTE_constraints( TE, Treadout, T_RF180, deltaMin )
% tri_Dd( TE, Treadout, T_RF180, deltaMin )
% display 2D space of feasible DELTA delta for the set of parameters in input
%
% Constraints :
% DELTA + delta + Tro/2 <= TE
% DELTA >= delta + T_RF180

DELTAMin = deltaMin + T_RF180;

x1 = deltaMin;
y1 = DELTAMin;

x2 = 0.5*(TE - Treadout - T_RF180);
y2 = -x2 + TE - Treadout;

x3 = deltaMin;
y3 = -deltaMin + TE - Treadout;

x = [x1 x2 x3];
y = [y1 y2 y3];

% fill(x, y, 'r')
% xlabel('\delta (in ms)')
% ylabel('\Delta (in ms)')
% title({'2D space of feasible (\Delta,\delta) combinations' ; ['TE : ', num2str(TE),' ms', '   ','Treadout : ', num2str(Treadout),' ms', '   ']});

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function loadDirs()
	% Forces all directories to be loaded, including external programs normally not loaded.
	% Biafra Ahanonu
	% started: 2020.05.09 [18:39:52]
	% inputs
		%
	% outputs
		%

	% changelog
		%
	% TODO
		%

	%========================
	% DESCRIPTION
	% options.exampleOption = '';
	% get options
	% options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	try
		ciapkg.loadBatchFxns;
	catch err
		disp(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		disp(repmat('@',1,7))
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
%
%      This HDR Toolbox demo creates an HDR radiance map:
%	   1) Read a stack of RAW images
%	   2) Read exposure values from the EXIF
%	   4) Build the radiance map using the stack and stack_exposure
%	   5) Save the radiance map in .hdr format
%	   6) Show the tone mapped version of the radiance map
%
%       Author: Francesco Banterle
%       Copyright 2015 (c)
%

clear all;

disp('Note: this demo does not have RAW files due to room limits in the repository');
disp('please add a folder with RAW files in the variable ''name_folder''.');
name_folder = ''; %Insert your raw stack folder
format = ''; %Insert your raw image format

disp('1) Read a stack of LDR images');
% negative number for saturation -> will be computed from the stack.
% With raw images this usually leads to artefacts in the highlights,
% so it is recommended to check the correct white level from image metadata
stack = ReadRAWStack(name_folder, format, -1);

disp('2) Read exposure values from the exif');
stack_exposure = ReadRAWStackInfo(name_folder, format);

disp('4) Build the radiance map using the stack and stack_exposure');
% Raw images are inherently linear, no need for CRF
imgHDR = BuildHDR(stack, stack_exposure, 'linear', [], 'Deb97', 'log');

disp('5) Save the radiance map in the .hdr format');
hdrimwrite(imgHDR, 'demos/output/stack_hdr_image_raw.exr');

disp('6) Show the tone mapped version of the radiance map with gamma encoding');
h = figure(2);
set(h, 'Name', 'Tone mapped version of the built HDR image');
imgTMO = GammaTMO(ReinhardTMO(imgHDR, 0.18), 2.2, 0, 1);
imwrite(imgTMO, 'demos/output/stack_hdr_image_tmo.png');

author: banterle


text:
function disp(x)
%DISP         Display function for pop-up windows in debugger
%
%Displays x without header
%

% written  04/26/13     S.M. Rump
%

  display(x,[],1);