text:
% Author: Yong Yuan
% Homepage: yongyuan.name

clear all;close all;clc;

% version: matconvnet-1.0-beta17
%run ./matconvnet-1.0-beta17/matlab/vl_compilenn
run ./matconvnet-1.0-beta17/matlab/vl_setupnn

%% Step 1 lOADING PATHS
path_imgDB = './facesDataset/';
addpath(path_imgDB);
addpath tools;

% viesion: matconvnet-1.0-beta17
net = load('vgg-face.mat') ;

%% Step 2 LOADING IMAGE AND EXTRACTING FEATURE
imgFiles = dir(path_imgDB);
imgNamList = {imgFiles(~[imgFiles.isdir]).name};
clear imgFiles;
imgNamList = imgNamList';

numImg = length(imgNamList);
feat = [];
rgbImgList = {};

%parpool;

%parfor i = 1:numImg
for i = 1:numImg
   oriImg = imread(imgNamList{i, 1});
   featVec = extractCNN(oriImg, net);
   feat = [feat; featVec'];
   fprintf('extract %d image\n\n', i);
end

feat_norm = normalize1(feat);
save('feat4096Norml.mat','feat_norm', 'imgNamList', '-v7.3');

author: willard-yuan

repo: CNN-for-Face-Image-Retrieval


text:
function [sAllAtlas, Messages] = import_label(SurfaceFile, LabelFiles, isNewAtlas, GridLoc, FileFormat)
% IMPORT_LABEL: Import an atlas segmentation for a given surface
%
% USAGE: import_label(SurfaceFile, LabelFiles, isNewAtlas=1, GridLoc=[], FileFormat=[]) : Add label information to SurfaceFile
%        import_label(SurfaceFile)                                                      : Ask the user for the label file to import

% @=============================================================================
% 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, 2012-2022

import sun.misc.BASE64Decoder;

%% ===== GET FILES =====
sAllAtlas = repmat(db_template('Atlas'), 0);
Messages = [];
% Parse inputs
if (nargin < 5) || isempty(FileFormat)
    FileFormat = [];
end
if (nargin < 4) || isempty(GridLoc)
    GridLoc = [];
end
if (nargin < 3) || isempty(isNewAtlas)
    isNewAtlas = 1;
end
if (nargin < 2) || isempty(LabelFiles)
    LabelFiles = [];
end

% CALL: import_label(SurfaceFile)
if isempty(LabelFiles)
    % Get last used folder
    LastUsedDirs = bst_get('LastUsedDirs');
    DefaultFormats = bst_get('DefaultFormats');
    % Get label files
    [LabelFiles, FileFormat] = java_getfile( 'open', ...
       'Import labels...', ...        % Window title
       LastUsedDirs.ImportAnat, ...   % Default directory
       'multiple', 'files', ...       % Selection mode
       bst_get('FileFilters', 'labelin'), ...
       DefaultFormats.LabelIn);
    % If no file was selected: exit
    if isempty(LabelFiles)
        return
    end
	% Save last used dir
    LastUsedDirs.ImportAnat = bst_fileparts(LabelFiles{1});
    bst_set('LastUsedDirs',  LastUsedDirs);
    % Save default export format
    DefaultFormats.LabelIn = FileFormat;
    bst_set('DefaultFormats',  DefaultFormats);
    % Warning if trying to import volume atlases on surfaces
    if isempty(GridLoc) && ismember(FileFormat, {'MRI-MASK', 'MRI-MASK-MNI', 'MRI-MASK-NOOVERLAP', 'MRI-MASK-NOOVERLAP-MNI'})
        isConfirm = java_dialog('confirm', ['You are trying to import a volume parcellation on a cortex surface.' 10 ...
            'This usually results in very poor results, with incomplete cortex parcels.' 10 ...
            'Using volume parcellations is only recommended with volume source models.' 10 10 ...
            'Proceed anyways?'], 'Warning');
        if ~isConfirm
            return;
        end
    end
% CALL: import_label(SurfaceFile, LabelFiles, ...)
else
    % Force cell input
    if ~iscell(LabelFiles)
        LabelFiles = {LabelFiles};
    end
    % Detect file format based on file extension
    if isempty(FileFormat)
        [fPath, fBase, fExt] = bst_fileparts(LabelFiles{1});
        switch (fExt)
            case '.annot',  FileFormat = 'FS-ANNOT';
            case '.label',  FileFormat = 'FS-LABEL-SINGLE';
            case '.gii',    FileFormat = 'GII-TEX';
            case '.dfs',    FileFormat = 'DFS';
            case '.dset',   FileFormat = 'DSET';
            case '.mat'
                switch (file_gettype(LabelFiles{1}))
                    case 'scout',        FileFormat = 'BST';
                    case 'subjectimage', FileFormat = 'MRI-MASK-NOOVERLAP';
                    otherwise,           Messages = 'Unsupported Brainstorm file type.'; return;
                end
            otherwise,  Messages = 'Unknown file extension.'; return;
        end
    end
end

%% ===== READ FILES =====
% Read destination surface
sSurf = bst_memory('GetSurface', file_short(SurfaceFile));
if isempty(sSurf)
    isLoadedHere = 1;
    sSurf = bst_memory('LoadSurface', file_short(SurfaceFile));
    panel_scout('SetCurrentSurface', sSurf.FileName);
else
    isLoadedHere = 0;
end
% Process one after the other
for iFile = 1:length(LabelFiles)
    LabelsTable = [];
    % Get updated atlases after the first iteration
    if (iFile > 1)
        sSurf = bst_memory('LoadSurface', file_short(SurfaceFile));
    end
    % Use the grid of points provided in input if available
    if ~isempty(GridLoc)
        Vertices = GridLoc;
        isVolumeAtlas = 1;
    % Otherwise, use the vertices of the cortex surface
    else
        Vertices = sSurf.Vertices;
        isVolumeAtlas = 0;
    end
    % Get filename
    [fPath, fBase, fExt] = bst_fileparts(LabelFiles{iFile});
    % Is it a Brainstorm volume atlas from the subject anatomy?
    isVolatlas = strcmpi(fExt, '.mat') && ~isempty(strfind(fBase, '_volatlas'));
    % New atlas structure: use filename as the atlas name
    if isNewAtlas || isempty(sSurf.Atlas) || isempty(sSurf.iAtlas)
        sAtlas = db_template('Atlas');
        iAtlas = 'Add';
        % Volume sources file
        if ~isempty(GridLoc)
            if isVolatlas
                sAtlas.Name = sprintf('Volume %d: ', length(GridLoc));
            else
                sAtlas.Name = sprintf('Volume %d: %s', length(GridLoc), fBase);
            end
            sAtlas.Name = file_unique(sAtlas.Name, {sSurf.Atlas.Name});
        % Surface sources file
        else
            % Get atlas name for standard FreeSurfer and MarsAtlas files
            if isVolatlas
                sAtlas.Name = 'From volume: ';
            else
                sAtlas.Name = GetAtlasName(fBase);
            end
        end
    % Existing atlas structure
    else
        iAtlas = sSurf.iAtlas;
        sAtlas = sSurf.Atlas(iAtlas);
        % For volume source files: Can only import volume scouts
        if ~isempty(GridLoc)
            % Can only work with volume scouts
            [isVolumeAtlas, nAtlasGrid] = panel_scout('ParseVolumeAtlas', sAtlas.Name);
            if ~isVolumeAtlas
                Messages = [Messages, 'Error: You can only load volume scouts for this sources file.'];
                return;
            end
            % Check the number of sources
            if (length(GridLoc) ~= nAtlasGrid)
                Messages = [Messages, sprintf('Error: The number of grid points in this sources file (%d) does not match the selected atlas (%d).', length(GridLoc), nAtlasGrid)];
                return;
            end
        end
    end
    % Check that atlas have the correct structure
    if isempty(sAtlas.Scouts)
        sAtlas.Scouts = repmat(db_template('scout'), 0);
    end
    % Switch based on file format
    switch (FileFormat)
        % ===== FREESURFER ANNOT =====
        case 'FS-ANNOT'
            % === READ FILE ===
            % Read label file
            try
                [vertices, labels, colortable] = read_annotation(LabelFiles{iFile}, 0);
            catch
                Messages = [Messages, sprintf('%s: read_annotation crashed: %s\n', [fBase, fExt], lasterr)];
                continue
            end
            % Check sizes
            if (length(labels) ~= length(Vertices))
                Messages = [Messages, sprintf('%s: Number of vertices in the surface (%d) and the label file (%d) do not match.\n', [fBase, fExt], length(Vertices), length(labels))];
                continue
            end

            % === CONVERT TO SCOUTS ===
            % Convert to scouts structures
            lablist = unique(labels);
            % Loop on each label
            for i = 1:length(lablist)
                % Find entry in the colortable
                iTable = find(colortable.table(:,5) == lablist(i));
                % If correspondence not defined: ignore label
                if (length(iTable) ~= 1)
                    continue;
                end
                % New scout index
                iScout = length(sAtlas.Scouts) + 1;
                sAtlas.Scouts(iScout).Vertices = find(labels == lablist(i))';
                if ~isempty(colortable.struct_names{iTable})
                    % Strip uselss parts of Schaeffer labels
                    Label = colortable.struct_names{iTable};
                    Label = strrep(Label, '17Networks_LH_', '');
                    Label = strrep(Label, '17Networks_RH_', '');
                    Label = strrep(Label, '7Networks_LH_', '');
                    Label = strrep(Label, '7Networks_RH_', '');
                    sAtlas.Scouts(iScout).Label = file_unique(Label, {sAtlas.Scouts.Label});
                else
                    sAtlas.Scouts(iScout).Label = file_unique('Unknown', {sAtlas.Scouts.Label});
                end
                sAtlas.Scouts(iScout).Color    = colortable.table(iTable,1:3) ./ 255;
                sAtlas.Scouts(iScout).Function = 'Mean';
                sAtlas.Scouts(iScout).Region   = 'UU';
            end
            if isempty(sAtlas.Scouts)
                Messages = [Messages, fBase, ': Could not match labels and color table.' 10];
                continue;
            end

        % ==== FREESURFER LABEL ====
        case {'FS-LABEL', 'FS-LABEL-SINGLE'}
            % === READ FILE ===
            % Read label file
            LabelMat = mne_read_label_file(LabelFiles{iFile});
            % Convert indices from 0-based to 1-based
            LabelMat.vertices = LabelMat.vertices + 1;
            % Check sizes
            if (max(LabelMat.vertices) > length(Vertices))
                Messages = [Messages, sprintf('%s: Numbers of vertices in the label file (%d) exceeds the number of vertices in the surface (%d)\n', fBase, max(LabelMat.vertices), length(Vertices))];
                continue
            end
            % === CONVERT TO SCOUTS ===
            % Number of ROIs
            if strcmpi(FileFormat, 'FS-LABEL-SINGLE')
                uniqueValues = 1;
            else
                uniqueValues = unique(LabelMat.values);
                minmax = [min(uniqueValues), max(uniqueValues)];
            end
            % Loop on each label
            for i = 1:length(uniqueValues)
                % New scout index
                iScout = length(sAtlas.Scouts) + 1;
                % Single ROI
                if strcmpi(FileFormat, 'FS-LABEL-SINGLE')
                    ScoutVert = sort(double(LabelMat.vertices));
                    Label = GetAtlasName(fBase);
                    Color = [];
                % Probability map
                else
                    % Calculate intensity [0,1]
                    if (minmax(1) == minmax(2))
                        c = 0;
                    else
                        c = (uniqueValues(i) - minmax(1)) ./ (minmax(2) - minmax(1));
                    end
                    ScoutVert = sort(double(LabelMat.vertices(LabelMat.values == uniqueValues(i))));
                    Label = file_unique(num2str(uniqueValues(i)), {sAtlas.Scouts.Label});
                    Color = [1 c 0];
                end
                % Create structure
                sAtlas.Scouts(iScout).Vertices = ScoutVert(:)';
                sAtlas.Scouts(iScout).Seed     = [];
                sAtlas.Scouts(iScout).Label    = Label;
                sAtlas.Scouts(iScout).Color    = Color;
                sAtlas.Scouts(iScout).Function = 'Mean';
                sAtlas.Scouts(iScout).Region   = 'UU';
            end
            if isempty(sAtlas.Scouts)
                Messages = [Messages, fBase, ': Could not match labels and color table.' 10];
                continue;
            end

        % ==== BRAINVISA GIFTI =====
        case 'GII-TEX'
            % Remove the "L" and "R" strings from the name
            AtlasName = sAtlas.Name;
            AtlasName = strrep(AtlasName, 'R', '');
            AtlasName = strrep(AtlasName, 'L', '');
            % Get labels
            LabelsTable = mri_getlabels(LabelFiles{iFile});
            % Read .gii file
            [sXml, Values] = in_gii(LabelFiles{iFile});
            % If there is more than one entry: force adding
            if (length(Values) > 1)
                iAtlas = 'Add';
            end
            % Process all the entries
            for ia = 1:length(Values)
                % Atlas name
                if (length(Values) > 1) && isNewAtlas
                    sAtlas(ia).Name = sprintf('%s #%d', AtlasName, ia);
                    sAtlas(ia).Scouts = repmat(db_template('scout'), 0);
                end
                % Check sizes
                if (length(Values{ia}) ~= length(Vertices))
                    Messages = [Messages, sprintf('%s: Entry #%d: Numbers of vertices in the surface (%d) and the label file (%d) do not match\n', fBase, ia, length(Vertices), length(Values{ia}))];
                    continue;
                elseif (all(size(Values{ia}) > 1))
                    Messages = [Messages, sprintf('%s: Entry #%d: Not a list of labels\n', fBase, ia)];
                    continue;
                end
                % Round the label values
                Values{ia} = round(Values{ia} * 1e3) / 1e3;
                % Convert to scouts structures
                lablist = unique(Values{ia});
                % Loop on each label
                for i = 1:length(lablist)
                    % Scout label: atlas or simply value converted to strinf
                    if ~isempty(LabelsTable) && ismember(lablist(i), [LabelsTable{:,1}])
                        scoutLabel = LabelsTable{lablist(i) == [LabelsTable{:,1}], 2};
                        scoutColor = LabelsTable{lablist(i) == [LabelsTable{:,1}], 3} ./ 255;
                    else
                        scoutLabel = num2str(lablist(i));
                        scoutColor = [];
                    end
                    % New scout index
                    iScout = length(sAtlas(ia).Scouts) + 1;
                    % Get the vertices for this annotation
                    sAtlas(ia).Scouts(iScout).Vertices = reshape(find(Values{ia} == lablist(i)), 1, []);
                    sAtlas(ia).Scouts(iScout).Seed     = [];
                    sAtlas(ia).Scouts(iScout).Label    = file_unique(scoutLabel, {sAtlas(ia).Scouts.Label});
                    sAtlas(ia).Scouts(iScout).Color    = scoutColor;
                    sAtlas(ia).Scouts(iScout).Function = 'Mean';
                    sAtlas(ia).Scouts(iScout).Region   = 'UU';
                end
            end

        % ===== SUMA DSET ROIs =====
        case 'DSET'
            % Read file
            sAtlas.Scouts = in_label_dset(LabelFiles{iFile});
            % Force adding new atlas
            iAtlas = 'Add';

        % ===== MRI VOLUMES =====
        case {'MRI-MASK', 'MRI-MASK-MNI', 'MRI-MASK-NOOVERLAP', 'MRI-MASK-NOOVERLAP-MNI'}
            bst_progress('text', 'Reading atlas...');
            % If the file that is loaded has to be interpreted in MNI space
            isMni = strcmpi(FileFormat, 'MRI-MASK-MNI') || strcmpi(FileFormat, 'MRI-MASK-NOOVERLAP-MNI');
            isOverlap = strcmpi(FileFormat, 'MRI-MASK') || strcmpi(FileFormat, 'MRI-MASK-MNI');
            % Read MRI volume  (do not normalize values when reading an atlas)
            if isVolatlas
                sMriMask = in_mri_bst(LabelFiles{iFile});
                sAtlas.Name = [sAtlas.Name, str_remove_parenth(sMriMask.Comment)];
            elseif isMni
                sMriMask = in_mri(LabelFiles{iFile}, 'ALL-MNI', [], 0);
            else
                sMriMask = in_mri(LabelFiles{iFile}, 'ALL', [], 0);
            end
            if isempty(sMriMask)
                return;
            end
            % Select only the first volume (if more than one)
            sMriMask.Cube = double(sMriMask.Cube(:,:,:,1));
            % Display warning when no MNI transformation available
            if isMni && (~isfield(sMriMask, 'NCS') || ...
                ((~isfield(sMriMask.NCS, 'R') || ~isfield(sMriMask.NCS, 'T') || isempty(sMriMask.NCS.R) || isempty(sMriMask.NCS.T)) && ...
                 (~isfield(sMriMask.NCS, 'iy') || isempty(sMriMask.NCS.iy))))
                isMni = 0;
                disp('Error: No MNI normalization available in this file.');
            end
            % Get all the values in the MRI
            bst_progress('text', 'Extract regions...');
            allValues = unique(sMriMask.Cube);
            % If values are not integers, it is not a mask or an atlas: it has to be binarized first
            if any(allValues ~= round(allValues))
                % Warning: not a binary mask
                isConfirm = java_dialog('confirm', ['Warning: This is not a binary mask.' 10 'Try to import this MRI as a surface anyway?'], 'Import binary mask');
                if ~isConfirm
                    return;
                end
                % Analyze MRI histogram
                Histogram = mri_histogram(sMriMask.Cube);
                % Binarize based on background level
                sMriMask.Cube = (sMriMask.Cube > Histogram.bgLevel);
                allValues = [0,1];
            end
            % Skip the first value (background)
            allValues(1) = [];
            % Load the subject SCS coordinates
            sSubject = bst_get('SurfaceFile', SurfaceFile);
            sMriSubj = in_mri_bst(sSubject.Anatomy(sSubject.iAnatomy).FileName);
            % Converting the sSurf vertices to the loaded MRI volume
            if isMni
                % The original volume is in subject space and and the atlas volume is in MNI space
                vertMni = cs_convert(sMriSubj, 'scs', 'mni', Vertices);
                if isempty(vertMni)
                    Messages = [Messages, 'Error: Compute the MNI transformation first.'];
                    return
                end
                vertMri = round(cs_convert(sMriMask, 'mni', 'voxel', vertMni));
            else
                % Check the compatibility of MRI sizes
                if ~isequal(size(sMriSubj.Cube), size(sMriMask.Cube))
                    Messages = [Messages, 'Error: The selected MRI file does not match the size of the subject''s MRI.' 10 10 ...
                                          'If the MRI is an atlas in a normalized MNI space, use the file format ' 10 ...
                                          '"Volume mask or atlas (MNI space)".' 10 10];
                    return;
                end
                % Both the original volume and the atlas volume and in the same subject space
                vertMri = round(cs_convert(sMriSubj, 'scs', 'voxel', Vertices));
            end
            % Get the corresponding MRI labels
            indMri = sub2ind(size(sMriMask.Cube), ...
                        bst_saturate(vertMri(:,1), [1, size(sMriMask.Cube,1)]), ...
                        bst_saturate(vertMri(:,2), [1, size(sMriMask.Cube,2)]), ...
                        bst_saturate(vertMri(:,3), [1, size(sMriMask.Cube,3)]));
            % Try to get volume labels for this atlas
            if isVolatlas
                VolumeLabels = sMriMask.Labels;
            else
                VolumeLabels = mri_getlabels(LabelFiles{iFile});
            end
            % Create one scout for each value in the volume
            for i = 1:length(allValues)
                bst_progress('text', sprintf('Creating scouts... [%d/%d]', i, length(allValues)));
                % Get the binary mask of the current region
                mask = (sMriMask.Cube == allValues(i));
                % Dilate mask
                if isOverlap
                    mask = mri_dilate(mask);
                    mask = mri_dilate(mask);
                end
                % Get the vertices in this mask
                iScoutVert = find(mask(indMri));
                if isempty(iScoutVert)
                    continue;
                end

                % New scout index
                iScout = length(sAtlas.Scouts) + 1;
                sAtlas.Scouts(iScout).Seed     = [];
                sAtlas.Scouts(iScout).Color    = [];
                sAtlas.Scouts(iScout).Function = 'Mean';
                sAtlas.Scouts(iScout).Region   = 'UU';
                % Use an existing label for the loaded atlas
                if ~isempty(VolumeLabels)
                    iLabel = find([VolumeLabels{:,1}] == allValues(i));
                    if ~isempty(iLabel)
                        sAtlas.Scouts(iScout).Label = file_unique(VolumeLabels{iLabel,2}, {sAtlas.Scouts.Label});
                        if (length(VolumeLabels{iLabel,2}) > 2)
                            if (VolumeLabels{iLabel,2}(end-1:end) == ' L')
                                sAtlas.Scouts(iScout).Region = 'LU';
                            elseif (VolumeLabels{iLabel,2}(end-1:end) == ' R')
                                sAtlas.Scouts(iScout).Region = 'RU';
                            end
                        end
                        if (size(VolumeLabels,2) >= 3) && (length(VolumeLabels{iLabel,3}) == 3)
                            sAtlas.Scouts(iScout).Color = VolumeLabels{iLabel,3} ./ 255;
                        end
                    else
                        sAtlas.Scouts(iScout).Label = file_unique(num2str(allValues(i)), {sAtlas.Scouts.Label});
                    end
                else
                    sAtlas.Scouts(iScout).Label = file_unique(num2str(allValues(i)), {sAtlas.Scouts.Label});
                end

                % If importing an atlas in MNI coordinates
                if isMni
                    % Get the coordinates of all the points in the mask
                    Pvox = [];
                    [Pvox(:,1),Pvox(:,2),Pvox(:,3)] = ind2sub(size(mask), indMri(iScoutVert));
                    % Convert to MNI coordinates
                    Pmni = cs_convert(sMriMask, 'voxel', 'mni', Pvox);
                    % Find left and right areas
                    iL = find(Pmni(:,1) < 0);
                    iR = find(Pmni(:,1) >= 0);
                    % If there is about the same number of points in the two sides: consider it's a bilateral region
                    if isempty(iL) || isempty(iR)
                        isSplit = 0;
                    elseif (length(iL) / length(iR) < 1.6) && (length(iL) / length(iR) > 0.4)
                        isSplit = 1;
                    % If the name starts with "S_", "G_" or "N_", it's also a bilateral region
                    elseif (length(sAtlas.Scouts(iScout).Label) >= 3) && ismember(sAtlas.Scouts(iScout).Label(1:2), {'S_', 'G_', 'N_'})
                        isSplit = 1;
                    else
                        isSplit = 0;
                    end
                    % If this is a bilateral region: split in two
                    if isSplit
                        % Duplicate scout
                        sAtlas.Scouts(iScout+1) = sAtlas.Scouts(iScout);
                        % Left scout
                        sAtlas.Scouts(iScout).Vertices = reshape(iScoutVert(iL), 1, []);
                        sAtlas.Scouts(iScout).Label    = [sAtlas.Scouts(iScout).Label, ' L'];
                        sAtlas.Scouts(iScout).Region   = 'LU';
                        % Right scout
                        sAtlas.Scouts(iScout+1).Vertices = reshape(iScoutVert(iR), 1, []);
                        sAtlas.Scouts(iScout+1).Label    = [sAtlas.Scouts(iScout+1).Label, ' R'];
                        sAtlas.Scouts(iScout+1).Region   = 'RU';
                    % Do not split
                    else
                        sAtlas.Scouts(iScout).Vertices = iScoutVert(:)';
                    end
                % Importing regular subject volume: Select all the vertices
                else
                    sAtlas.Scouts(iScout).Vertices = iScoutVert(:)';
                end
            end
            bst_progress('text', 'Saving atlas...');

        % ===== BRAINSTORM SCOUTS =====
        case 'BST'
            % Load file
            ScoutMat = load(LabelFiles{iFile});
            % Convert old scouts structure to new one
            if isfield(ScoutMat, 'Scout')
                ScoutMat.Scouts = ScoutMat.Scout;
            elseif isfield(ScoutMat, 'Scouts')
                % Ok
            else
                Messages = [Messages, fBase, ': Invalid scouts file.' 10];
                continue;
            end
            % Check the number of vertices
            if ~isVolumeAtlas && (length(Vertices) ~= ScoutMat.TessNbVertices)
                Messages = [Messages, sprintf('%s: Numbers of vertices in the surface (%d) and the scout file (%d) do not match\n', fBase, length(Vertices), ScoutMat.TessNbVertices)];
                continue;
            end
            % If name is not defined: use the filename
            if isNewAtlas
                if isfield(ScoutMat, 'Name') && ~isempty(ScoutMat.Name)
                    sAtlas.Name = ScoutMat.Name;
                else
                    [fPath,fBase] = bst_fileparts(LabelFiles{iFile});
                    sAtlas.Name = strrep(fBase, 'scout_', '');
                end
            end
            % Copy the new scouts
            for i = 1:length(ScoutMat.Scouts)
                iScout = length(sAtlas.Scouts) + 1;
                sAtlas.Scouts(iScout).Vertices = ScoutMat.Scouts(i).Vertices(:)';
                sAtlas.Scouts(iScout).Seed     = ScoutMat.Scouts(i).Seed;
                sAtlas.Scouts(iScout).Color    = ScoutMat.Scouts(i).Color;
                sAtlas.Scouts(iScout).Label    = file_unique(ScoutMat.Scouts(i).Label, {sAtlas.Scouts.Label});
                sAtlas.Scouts(iScout).Function = ScoutMat.Scouts(i).Function;
                if isfield(ScoutMat.Scouts(i), 'Region')
                    sAtlas.Scouts(iScout).Region = ScoutMat.Scouts(i).Region;
                else
                    sAtlas.Scouts(iScout).Region = 'UU';
                end
            end

        % ===== BrainSuite/SVReg surface file =====
        case 'DFS'
            % === READ FILE ===
            [VertexLabelIds, labelMap, AtlasName] = in_label_bs(LabelFiles{iFile});

            % === CONVERT TO SCOUTS ===
            % Convert to scouts structures
            lablist = unique(VertexLabelIds);
            if isNewAtlas
                if isempty(AtlasName)
                    sAtlas.Name = 'SVReg';
                else
                    sAtlas.Name = AtlasName;
                end
            end

            % Loop on each label
            if isempty(labelMap)
                new_colors = distinguishable_colors(length(lablist), [0,0,0]);
                new_colors(1,:) = [0.5,0.5,0.5];
            end

            for i = 1:length(lablist)
                % Find label ID
                id = lablist(i);
                % Skip if label id is not in labelMap
                if isempty(labelMap)
                    labelInfo.Name = num2str(id);
                    labelInfo.Color = new_colors(i,:)';
                    %continue;
                elseif ~labelMap.containsKey(num2str(id))
                    labelInfo.Name = num2str(id);
                    labelInfo.Color = [0,0,0]';
                else
                    entry = labelMap.get(num2str(id));
                    labelInfo.Name = entry(1);
                    labelInfo.Color = entry(2);
                end
                % Transpose color vector
                labelInfo.Color = labelInfo.Color(:)';
                % Skip the "background" scout
                if strcmpi(labelInfo.Name, 'background')
                    labelInfo.Color = [0.5,0.5,0.5]';
                    continue;
                end
                % New scout index
                iScout = length(sAtlas.Scouts) + 1;
                sAtlas.Scouts(iScout).Vertices = reshape(find(VertexLabelIds == id), 1, []);
                sAtlas.Scouts(iScout).Label    = file_unique(labelInfo.Name, {sAtlas.Scouts.Label});
                sAtlas.Scouts(iScout).Color    = labelInfo.Color;
                sAtlas.Scouts(iScout).Function = 'Mean';
                sAtlas.Scouts(iScout).Region   = 'UU';
            end
            if isempty(sAtlas.Scouts)
                Messages = [Messages, fBase, ': Could not match vertex labels and label description file.' 10];
                continue;
            end

        % ===== Unknown file =====
        otherwise
            Messages = [Messages, fBase, ': Unknown file extension.' 10];
            continue;
    end

    % Get the scouts colortable
    ColorTable = panel_scout('GetScoutsColorTable');
    % Loop on all the loaded atlases
    for ia = 1:length(sAtlas)
        if isempty(sAtlas(ia).Scouts)
            continue;
        end
        % Brodmann atlas: remove the "Unknown" scout
        iUnknown = find(strcmpi({sAtlas(ia).Scouts.Label}, 'unknown') | strcmpi({sAtlas(ia).Scouts.Label}, 'medial.wall') | strcmpi({sAtlas(ia).Scouts.Label}, 'freesurfer_defined_medial_wall'));
        if ~isempty(iUnknown)
            sAtlas(ia).Scouts(iUnknown) = [];
        end
        % Fix all the scouts seeds and identify regions
        for i = 1:length(sAtlas(ia).Scouts)
            % Fix seed
            if isempty(sAtlas(ia).Scouts(i).Seed) || ~ismember(sAtlas(ia).Scouts(i).Seed, sAtlas(ia).Scouts(i).Vertices)
                sAtlas(ia).Scouts(i) = panel_scout('SetScoutsSeed', sAtlas(ia).Scouts(i), Vertices);
            end
            % Fix regions
            if isempty(sAtlas(ia).Scouts(i).Region) || strcmpi(sAtlas(ia).Scouts(i).Region, 'UU')
                sAtlas(ia).Scouts(i) = tess_detect_region(sAtlas(ia).Scouts(i));
            end
            % Fix color
            if isempty(sAtlas(ia).Scouts(i).Color)
                iColor = mod(i-1, length(ColorTable)) + 1;
                sAtlas(ia).Scouts(i).Color = ColorTable(iColor,:);
            end
        end
        % Sort scouts by alphabetical order
        if ~strcmpi(FileFormat, 'BST')
            [tmp,I] = sort(lower({sAtlas(ia).Scouts.Label}));
            sAtlas(ia).Scouts = sAtlas(ia).Scouts(I);
        end
        % Return new atlas
        sAllAtlas(end+1) = sAtlas(ia);
        % Add atlas to the surface
        panel_scout('SetAtlas', SurfaceFile, iAtlas, sAtlas(ia));
    end
end

%% ===== SAVE IN SURFACE =====
% Unload surface to save it
if isLoadedHere
    bst_memory('UnloadSurface', SurfaceFile);
end

end

%% ===== GET ATLAS NAME =====
function AtlasName = GetAtlasName(fBase)
    fBase = lower(fBase);
    switch (fBase)
        case {'lh.aparc.a2009s', 'rh.aparc.a2009s'}
            AtlasName = 'Destrieux';
        case {'lh.aparc', 'rh.aparc'}
            AtlasName = 'Desikan-Killiany';
        case {'lh.ba', 'rh.ba', 'lh.ba_exvivo', 'rh.ba_exvivo'}
            AtlasName = 'Brodmann';
        case {'lh.ba.thresh', 'rh.ba.thresh', 'lh.ba_exvivo.thresh', 'rh.ba_exvivo.thresh'}
            AtlasName = 'Brodmann-thresh';
        case {'lh.aparc.dktatlas40', 'rh.aparc.dktatlas40'}
            AtlasName = 'DKT40';
        case {'lh.aparc.dktatlas', 'rh.aparc.dktatlas', 'lh.aparc.mapped', 'rh.aparc.mapped'}  % FreeSurfer, FastSurfer
            AtlasName = 'DKT';
        case {'lh.pals_b12_brodmann', 'rh.pals_b12_brodmann'}
            AtlasName = 'PALS-B12 Brodmann';
        case {'lh.pals_b12_lobes', 'rh.pals_b12_lobes'}
            AtlasName = 'PALS-B12 Lobes';
        case {'lh.pals_b12_orbitofrontal', 'rh.pals_b12_orbitofrontal'}
            AtlasName = 'PALS-B12 Orbito-frontal';
        case {'lh.pals_b12_visuotopic', 'rh.pals_b12_visuotopic'}
            AtlasName = 'PALS-B12 Visuotopic';
        case {'lh.yeo2011_7networks_n1000', 'rh.yeo2011_7networks_n1000'}
            AtlasName = 'Yeo 7 Networks';
        case {'lh.yeo2011_17networks_n1000', 'rh.yeo2011_17networks_n1000'}
            AtlasName = 'Yeo 17 Networks';
        case {'lh.prf', 'rh.prf'}
            AtlasName = 'Retinotopy';
        case {'lh.myaparc_36', 'rh.myaparc_36'}
            AtlasName = 'Lausanne-S33';
        case {'lh.myaparc_60', 'rh.myaparc_60'}
            AtlasName = 'Lausanne-S60';
        case {'lh.myaparc_125', 'rh.myaparc_125'}
            AtlasName = 'Lausanne-S125';
        case {'lh.myaparc_250', 'rh.myaparc_250'}
            AtlasName = 'Lausanne-S250';
        case {'lh.bn_atlas', 'rh.bn_atlas'}
            AtlasName = 'Brainnetome';
        case {'lh.oasis.chubs', 'rh.oasis.chubs'}
            AtlasName = 'OASIS cortical hubs';
        case {'lh.mpm.vpnl', 'rh.mpm.vpnl'}
            AtlasName = 'vcAtlas';
        otherwise
            % CAT12 / FreeSurfer
            if ~isempty(strfind(fBase, 'aparc_a2009s'))
                AtlasName = 'Destrieux';
            elseif ~isempty(strfind(fBase, 'dk40'))
                AtlasName = 'Desikan-Killiany';
            elseif ~isempty(strfind(fBase, 'hcp_mmp')) || ~isempty(strfind(fBase, 'hcp-mmp'))
                AtlasName = 'HCP_MMP1';
            elseif ~isempty(strfind(fBase, 'schaefer2018_100parcels_17'))
                AtlasName = 'Schaefer_100_17net';
            elseif ~isempty(strfind(fBase, 'schaefer2018_200parcels_17'))
                AtlasName = 'Schaefer_200_17net';
            elseif ~isempty(strfind(fBase, 'schaefer2018_400parcels_17'))
                AtlasName = 'Schaefer_400_17net';
            elseif ~isempty(strfind(fBase, 'schaefer2018_600parcels_17'))
                AtlasName = 'Schaefer_600_17net';
            elseif ~isempty(strfind(fBase, 'schaefer2018_800parcels_17'))
                AtlasName = 'Schaefer_800_17net';
            elseif ~isempty(strfind(fBase, 'schaefer2018_1000parcels_17'))
                AtlasName = 'Schaefer_1000_17net';
            elseif ~isempty(strfind(fBase, 'schaefer2018_100parcels_7'))
                AtlasName = 'Schaefer_100_7net';
            elseif ~isempty(strfind(fBase, 'schaefer2018_200parcels_7'))
                AtlasName = 'Schaefer_200_7net';
            elseif ~isempty(strfind(fBase, 'schaefer2018_400parcels_7'))
                AtlasName = 'Schaefer_400_7net';
            elseif ~isempty(strfind(fBase, 'schaefer2018_600parcels_7'))
                AtlasName = 'Schaefer_600_7net';
            elseif ~isempty(strfind(fBase, 'schaefer2018_800parcels_7'))
                AtlasName = 'Schaefer_800_7net';
            elseif ~isempty(strfind(fBase, 'schaefer2018_1000parcels_7'))
                AtlasName = 'Schaefer_1000_7net';
            % FreeSurfer left/right
            elseif (length(fBase) > 3) && (strcmpi(fBase(1:3), 'lh.') || strcmpi(fBase(1:3), 'rh.'))
                AtlasName = fBase(4:end);
            % BrainVISA/MarsAtlas
            elseif (~isempty(strfind(fBase, '_lwhite_parcels_marsatlas')) || ~isempty(strfind(fBase, '_rwhite_parcels_marsatlas')))
                AtlasName = 'MarsAtlas';
            elseif (~isempty(strfind(fBase, '_lwhite_parcels_model')) || ~isempty(strfind(fBase, '_rwhite_parcels_model')))
                AtlasName = 'MarsAtlas model';
            elseif (~isempty(strfind(fBase, '_lwhite_pole_cingular')) || ~isempty(strfind(fBase, '_rwhite_pole_cingular')))
                AtlasName = 'MarsAtlas pole cingular';
            elseif (~isempty(strfind(fBase, '_lwhite_pole_insula')) || ~isempty(strfind(fBase, '_rwhite_pole_insula')))
                AtlasName = 'MarsAtlas pole insula';
            elseif (~isempty(strfind(fBase, '_lwhite_sulcalines')) || ~isempty(strfind(fBase, '_rwhite_sulcalines')))
                AtlasName = 'MarsAtlas sulcal lines';
            else
                AtlasName = fBase;
            end
    end
end

author: brainstorm-tools


text:
%EVALSFUN Shape function evaluation driver subroutine.
%
%   [ VBASE, NLDOF, XLDOF, SFUN ] = EVALSFUN( SFUN, I_EVAL, N_SDIM, N_VERT, I_DOF, XI, AINVJAC, VBASE )
%   Evaluates local basis function I_DOF in point XI for shape function SFUN on all
%   cells corresponding to AINVJAC. Returns the evaluated values in VBASE and
%   the local number of degrees of freedom on vertices, edges, faces, and cell
%   interiors in NLDOF and XLDOF as well as the low level function name in SFUN.
%
%   See also EVALDVAR, EVALEXPR

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function msfun_realtime_pacer(block)
% Help for Writing Level-2 M-File S-Functions:
%   web([docroot '/toolbox/simulink/sfg/f7-67622.html']
%   http://www.mathworks.com/access/helpdesk/help/toolbox/simulink/sfg/f7-67622.html

% Copyright 2009, The MathWorks, Inc.

% instance variables
mySimTimePerRealTime = 1;
myRealTimeBaseline = 0;
mySimulationTimeBaseline = 0;
myResetBaseline = true;
myTotalBurnedTime = 0;
myNumUpdates = 0;

setup(block);

%% ---------------------------------------------------
    function setup(block)
        % Register the number of ports.
        block.NumInputPorts  = 0;
        block.NumOutputPorts = 0;

        % Set up the states
        block.NumContStates = 0;
        block.NumDworks = 0;

        % Register the parameters.
        block.NumDialogPrms     = 1; % scale factor
        block.DialogPrmsTunable = {'Nontunable'};

        % Block is fixed in minor time step, i.e., it is only executed on major
        % time steps. With a fixed-step solver, the block runs at the fastest
        % discrete rate.
        block.SampleTimes = [0 1];

        block.SetAccelRunOnTLC(false); % run block in interpreted mode even w/ Acceleration

        % methods called during update diagram/compilation.
        block.RegBlockMethod('CheckParameters', @CheckPrms);

        % methods called at run-time
        block.RegBlockMethod('Start', @Start);
        block.RegBlockMethod('Update', @Update);
        block.RegBlockMethod('SimStatusChange', @SimStatusChange);
        block.RegBlockMethod('Terminate', @Terminate);
    end

%%
    function CheckPrms(block)
        try
            validateattributes(block.DialogPrm(1).Data, {'double'},{'real', 'scalar', '>', 0});
        catch %#ok<CTCH>
            throw(MSLException(block.BlockHandle, ...
                'Simulink:Parameters:BlkParamUndefined', ...
                'Enter a number greater than 0'));
        end
    end

%%
    function Start(block)
        mySimTimePerRealTime = block.DialogPrm(1).Data;
        myTotalBurnedTime = 0;
        myNumUpdates = 0;
        myResetBaseline = true;
        if strcmp(pause('query'),'off')
            fprintf('%s: Enabling MATLAB PAUSE command\n', getfullname(block.BlockHandle));
            pause('on');
        end
    end

%%
    function Update(block)
        if  myResetBaseline
            myRealTimeBaseline = tic;
            mySimulationTimeBaseline = block.CurrentTime;
            myResetBaseline = false;
        else
            if isinf(mySimTimePerRealTime)
                return;
            end
            elapsedRealTime = toc(myRealTimeBaseline);
            differenceInSeconds = ((block.CurrentTime - mySimulationTimeBaseline) / mySimTimePerRealTime) - elapsedRealTime;
            if differenceInSeconds >= 0
                pause(differenceInSeconds);
                myTotalBurnedTime = myTotalBurnedTime + differenceInSeconds;
                myNumUpdates = myNumUpdates + 1;
            end
        end
    end

%%
    function SimStatusChange(block, status)
        if status == 0,
            % simulation paused
            fprintf('%s: Pausing real time execution of the model (simulation time = %g sec)\n', ...
                getfullname(block.BlockHandle), block.CurrentTime);
        elseif status == 1
            % Simulation resumed
            fprintf('%s: Continuing real time execution of the model\n', ...
                getfullname(block.BlockHandle));
            myResetBaseline = true;
        end
    end

%%
    function Terminate(block)
        if myNumUpdates > 0
            fprintf('%s: Average idle real time per major time step = %g sec\n', ...
                getfullname(block.BlockHandle),  myTotalBurnedTime / myNumUpdates);
        end
    end

end

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function [imgs, TR, scanspersess, names, onsets, durations] = spm_mat2batchinput(SPM)
% Extract info from SPM structure and format in batch input mode
%
% :Usage:
% ::
%
%     [imgs, TR, scanspersess, names, onsets, durations] = spm_mat2batchinput(SPM)
%
% :Input into SPM8::
%
%   tmod = time modulator (number per cell, order, linear = 1)
%
% For parametric modulation include a structure array, which is up to 1 x n in size, called pmod. n must be less
% than  or  equal  to  the  number  of  cells in the names/onsets/durations cell arrays. The structure array pmod
% must  have the fields: name, param and poly.  Each of these fields is in turn a cell array to allow the inclusion
% R for regressors of no interest
%
% NOTES:
% This is a preliminary version of this function and does not extract
% covariates, parametric modulators, etc.
%
% Returns output in the format you would use to enter into the SPM GUI to
% set up a model. Same format is used by canlab_spm_fmri_model_job.m
%
% ..
%    Tor Wager, Feb 2011
%    Could be extended to include nuisance covs, parametric modulators, etc.
%    Right now this is very basic.
% ..
%
% CANlab SPM 1st-level model helper tools:
% scn_spm_design_check, spm2canlab_dataset, spm_mat2batchinput, scn_spm_get_events_of_interest

if iscell(SPM.xY.P) % check input image names
    imgs = char(SPM.xY.P{:});
else
    imgs = SPM.xY.P;
end

TR = SPM.xY.RT;
scanspersess = SPM.nscan;

[names, onsets, durations, covs] = deal({});

% list regressor names
for sess = 1:length(SPM.Sess)

    names = [names cat(1, SPM.Sess(sess).U(:).name)'];
    %disp(sprintf('Session %3.0f', sess));

end

% check onsets
% 0 should be the start of the first image acquisition

for sess = 1:length(SPM.Sess)

    %disp(sprintf('Session %3.0f', sess));

    for evt = 1:length(SPM.Sess(sess).U)
        %         disp(['Session ' num2str(sess) ' event type: ' SPM.Sess(sess).U(evt).name])
        %         [SPM.Sess(sess).U(evt).ons  SPM.Sess(sess).U(evt).durations]

        onsets = [onsets {SPM.Sess(sess).U(evt).ons}];
        durations = [durations {SPM.Sess(sess).U(evt).dur}];

    end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function a = rutis1_eigen_right ( )

%*****************************************************************************80
%
%% RUTIS1_EIGEN_RIGHT returns the right eigenvectors of the RUTIS1 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A(4,4), the right eigenvector matrix.
%
  a(1:4,1:4) = [ ...
     1.0,  1.0,  1.0,  1.0; ...
     1.0,  0.0,  0.0, -1.0; ...
     0.0,  1.0, -1.0,  0.0; ...
     1.0, -1.0, -1.0,  1.0 ]';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%LINK Robot manipulator Link class
%
% A Link object holds all information related to a robot link such as
% kinematics parameters, rigid-body inertial parameters, motor and
% transmission parameters.
%
% Methods::
%  A             link transform matrix
%  RP            joint type: 'R' or 'P'
%  friction      friction force
%  nofriction    Link object with friction parameters set to zero
%  dyn           display link dynamic parameters
%  islimit       test if joint exceeds soft limit
%  isrevolute    test if joint is revolute
%  isprismatic   test if joint is prismatic
%  display       print the link parameters in human readable form
%  char          convert to string
%
% Properties (read/write)::
%
%  theta    kinematic: joint angle
%  d        kinematic: link offset
%  a        kinematic: link length
%  alpha    kinematic: link twist
%  sigma    kinematic: 0 if revolute, 1 if prismatic
%  mdh      kinematic: 0 if standard D&H, else 1
%  offset   kinematic: joint variable offset
%  qlim     kinematic: joint variable limits [min max]
%-
%  m        dynamic: link mass
%  r        dynamic: link COG wrt link coordinate frame 3x1
%  I        dynamic: link inertia matrix, symmetric 3x3, about link COG.
%  B        dynamic: link viscous friction (motor referred)
%  Tc       dynamic: link Coulomb friction
%-
%  G        actuator: gear ratio
%  Jm       actuator: motor inertia (motor referred)
%
% Examples::
%
%         L = Link([0 1.2 0.3 pi/2]);
%         L = Link('revolute', 'd', 1.2, 'a', 0.3, 'alpha', pi/2);
%         L = Revolute('d', 1.2, 'a', 0.3, 'alpha', pi/2);
%
% Notes::
% - This is a reference class object.
% - Link objects can be used in vectors and arrays.
%
% References::
% - Robotics, Vision & Control, Chap 7,
%   P. Corke, Springer 2011.
%
% See also Link, Revolute, Prismatic, SerialLink, RevoluteMDH, PrismaticMDH.

% Copyright (C) 1993-2015, 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

classdef Link < handle
    properties
        % kinematic parameters
        theta % kinematic: link angle
        d     % kinematic: link offset
        alpha % kinematic: link twist
        a     % kinematic: link length
        sigma % revolute=0, prismatic=1
        mdh   % standard DH=0, MDH=1
        offset % joint coordinate offset
        name   % joint coordinate name
        flip   % joint moves in opposite direction

        % dynamic parameters
        m  % dynamic: link mass
        r  % dynamic: position of COM with respect to link frame (3x1)
        I  % dynamic: inertia of link with respect to COM (3x3)
        Jm % dynamic: motor inertia
        B  % dynamic: motor viscous friction (1x1 or 2x1)

        Tc % dynamic: motor Coulomb friction (1x2 or 2x1)
        G  % dynamic: gear ratio
        qlim % joint coordinate limits (2x1)
    end

    methods
        function l = Link(varargin)
            %LINK Create robot link object
            %
            % This the class constructor which has several call signatures.
            %
            % L = Link() is a Link object with default parameters.
            %
            % L = Link(LNK) is a Link object that is a deep copy of the link
            % object LNK.
            %
            % L = Link(OPTIONS) is a link object with the kinematic and dynamic
            % parameters specified by the key/value pairs.
            %
            % Options::
            % 'theta',TH    joint angle, if not specified joint is revolute
            % 'd',D         joint extension, if not specified joint is prismatic
            % 'a',A         joint offset (default 0)
            % 'alpha',A     joint twist (default 0)
            % 'standard'    defined using standard D&H parameters (default).
            % 'modified'    defined using modified D&H parameters.
            % 'offset',O    joint variable offset (default 0)
            % 'qlim',L      joint limit (default [])
            % 'I',I         link inertia matrix (3x1, 6x1 or 3x3)
            % 'r',R         link centre of gravity (3x1)
            % 'm',M         link mass (1x1)
            % 'G',G         motor gear ratio (default 1)
            % 'B',B         joint friction, motor referenced (default 0)
            % 'Jm',J        motor inertia, motor referenced (default 0)
            % 'Tc',T        Coulomb friction, motor referenced (1x1 or 2x1), (default [0 0])
            % 'revolute'    for a revolute joint (default)
            % 'prismatic'   for a prismatic joint 'p'
            % 'standard'    for standard D&H parameters (default).
            % 'modified'    for modified D&H parameters.
            % 'sym'         consider all parameter values as symbolic not numeric
            %
            % - It is an error to specify both 'theta' and 'd'
            % - The link inertia matrix (3x3) is symmetric and can be specified by giving
            %   a 3x3 matrix, the diagonal elements [Ixx Iyy Izz], or the moments and products
            %   of inertia [Ixx Iyy Izz Ixy Iyz Ixz].
            % - All friction quantities are referenced to the motor not the load.
            % - Gear ratio is used only to convert motor referenced quantities such as
            %   friction and interia to the link frame.
            %
            % Old syntax::
            % L = Link(DH, OPTIONS) is a link object using the specified kinematic
            % convention  and with parameters:
            %  - DH = [THETA D A ALPHA SIGMA OFFSET] where OFFSET is a constant displacement
            %    between the user joint angle vector and the true kinematic solution.
            %  - DH = [THETA D A ALPHA SIGMA] where SIGMA=0 for a revolute and 1 for a
            %    prismatic joint, OFFSET is zero.
            %  - DH = [THETA D A ALPHA], joint is assumed revolute and OFFSET is zero.
            %
            % Options::
            %
            % 'standard'    for standard D&H parameters (default).
            % 'modified'    for modified D&H parameters.
            % 'revolute'    for a revolute joint, can be abbreviated to 'r' (default)
            % 'prismatic'   for a prismatic joint, can be abbreviated to 'p'
            %
            % Examples::
            % A standard Denavit-Hartenberg link
            %        L3 = Link('d', 0.15005, 'a', 0.0203, 'alpha', -pi/2);
            % since 'theta' is not specified the joint is assumed to be revolute, and
            % since the kinematic convention is not specified it is assumed 'standard'.
            %
            % Using the old syntax
            %        L3 = Link([ 0, 0.15005, 0.0203, -pi/2], 'standard');
            % the flag 'standard' is not strictly necessary but adds clarity.  Only 4 parameters
            % are specified so sigma is assumed to be zero, ie. the joint is revolute.
            %
            %        L3 = Link([ 0, 0.15005, 0.0203, -pi/2, 0], 'standard');
            % the flag 'standard' is not strictly necessary but adds clarity.  5 parameters
            % are specified and sigma is set to zero, ie. the joint is revolute.
            %
            %        L3 = Link([ 0, 0.15005, 0.0203, -pi/2, 1], 'standard');
            % the flag 'standard' is not strictly necessary but adds clarity.  5 parameters
            % are specified and sigma is set to one, ie. the joint is prismatic.
            %
            % For a modified Denavit-Hartenberg revolute joint
            %        L3 = Link([ 0, 0.15005, 0.0203, -pi/2, 0], 'modified');
            %
            % Notes::
            % - Link object is a reference object, a subclass of Handle object.
            % - Link objects can be used in vectors and arrays.
            % - The parameter D is unused in a revolute joint, it is simply a placeholder
            %   in the vector and the value given is ignored.
            % - The parameter THETA is unused in a prismatic joint, it is simply a placeholder
            %   in the vector and the value given is ignored.
            % - The joint offset is a constant added to the joint angle variable before
            %   forward kinematics and subtracted after inverse kinematics.  It is useful
            %   if you  want the robot to adopt a 'sensible' pose for zero joint angle
            %   configuration.
            % - The link dynamic (inertial and motor) parameters are all set to
            %   zero.  These must be set by explicitly assigning the object
            %   properties: m, r, I, Jm, B, Tc.
            % - The gear ratio is set to 1 by default, meaning that motor friction and
            %   inertia will be considered if they are non-zero.
            %
            % See also Revolute, Prismatic.

            %TODO eliminate legacy dyn matrix

            if nargin == 0
                % create an 'empty' Link object
                % this call signature is needed to support arrays of Links

                %% kinematic parameters
                l.alpha = 0;
                l.a = 0;
                l.theta = 0;
                l.d = 0;
                l.sigma = 0;
                l.mdh = 0;
                l.offset = 0;
                l.flip = false;
                l.qlim = [];

                %% dynamic parameters
                % these parameters must be set by the user if dynamics is used
                l.m = 0;
                l.r = [0 0 0];
                l.I = zeros(3,3);

                % dynamic params with default (zero friction)
                l.Jm = 0;
                l.G = 1;
                l.B = 0;
                l.Tc = [0 0];
            elseif nargin == 1 && isa(varargin{1}, 'Link')
                % clone the passed Link object
                l = copy(varargin{1});

            else
                % Create a new Link based on parameters

                % parse all possible options
                opt.theta = [];
                opt.a = 0;
                opt.d = [];
                opt.alpha = 0;
                opt.G = 0;
                opt.B = 0;
                opt.Tc = [0 0];
                opt.Jm = 0;
                opt.I = zeros(3,3);
                opt.m = 0;
                opt.r = [0 0 0];
                opt.offset = 0;
                opt.qlim = [];
                opt.type = {'revolute', 'prismatic', 'fixed'};
                opt.convention = {'standard', 'modified'};
                opt.sym = false;
                opt.flip = false;

                [opt,args] = tb_optparse(opt, varargin);

                % return a parameter as a number of symbol depending on
                % the 'sym' option

                if isempty(args)
                    % This is the new call format, where all parameters are
                    % given by key/value pairs
                    %
                    % eg. L3 = Link('d', 0.15005, 'a', 0.0203, 'alpha', -pi/2);
                    if ~isempty(opt.theta)
                        % constant value of theta means it must be prismatic
                        l.theta = value( opt.theta, opt);
                        l.sigma = 1;
                    end

                    if ~isempty(opt.d)
                        % constant value of d means it must be revolute
                        l.d = value( opt.d, opt);
                        l.sigma = 0;
                    end
                    if ~isempty(opt.d) && ~isempty(opt.theta)
                        error('RTB:Link:badarg', 'specify only one of ''d'' or ''theta''');
                    end

                    l.a =     value( opt.a, opt);
                    l.alpha = value( opt.alpha, opt);

                    l.offset = value( opt.offset, opt);
                    l.flip = value( opt.flip, opt);

                    l.qlim =   value( opt.qlim, opt);

                    l.m = value( opt.m, opt);
                    l.r = value( opt.r, opt);
                    l.I = value( opt.I, opt);
                    l.Jm = value( opt.Jm, opt);
                    l.G = value( opt.G, opt);
                    l.B = value( opt.B, opt);
                    l.Tc = value( opt.Tc, opt);
                else
                    % This is the old call format, where all parameters are
                    % given by a vector containing kinematic-only, or
                    % kinematic plus dynamic parameters
                    %
                    % eg. L3 = Link([ 0, 0.15005, 0.0203, -pi/2, 0], 'standard');
                    dh = args{1};
                    if length(dh) < 4
                        error('must provide params (theta d a alpha)');
                    end

                    % set the kinematic parameters
                    l.theta = dh(1);
                    l.d = dh(2);
                    l.a = dh(3);
                    l.alpha = dh(4);

                    l.sigma = 0;
                    l.offset = 0;
                    l.flip = false;
                    l.mdh = 0;  % default to standard D&H

                    % optionally set sigma and offset
                    if length(dh) >= 5
                        l.sigma = dh(5);
                    end
                    if length(dh) == 6
                        l.offset = dh(6);
                    end

                    if length(dh) > 6
                        % legacy DYN matrix

                        l.sigma = dh(5);
                        l.mdh = 0;  % default to standard D&H
                        l.offset = 0;

                        % it's a legacy DYN matrix
                        l.m = dh(6);
                        l.r = dh(7:9).';     % a column vector
                        v = dh(10:15);
                        l.I = [ v(1) v(4) v(6)
                            v(4) v(2) v(5)
                            v(6) v(5) v(3)];
                        if length(dh) > 15
                            l.Jm = dh(16);
                        end
                        if length(dh) > 16
                            l.G = dh(17);
                        else
                            l.G = 1;
                        end
                        if length(dh) > 17
                            l.B = dh(18);
                        else
                            l.B = 0.0;
                        end
                        if length(dh) > 18
                            l.Tc = dh(19:20);
                        else
                            l.Tc = [0 0];
                        end
                        l.qlim = [];
                    else
                        % we know nothing about the dynamics
                        l.m = [];
                        l.r = [];
                        l.I = [];
                        l.Jm = [];
                        l.G = 0;
                        l.B = 0;
                        l.Tc = [0 0];
                        l.qlim = [];
                    end
                end

                % set the kinematic convention to be used
                if strcmp(opt.convention, 'modified')
                    l.mdh = 1;
                else
                    l.mdh = 0;
                end

            end

            function out = value(v, opt)
                if opt.sym
                    out = sym(v);
                else
                    out = v;
                end
            end

        end % link()

        function  tau = friction(l, qd)
            %Link.friction Joint friction force
            %
            % F = L.friction(QD) is the joint friction force/torque for link velocity QD.
            %
            % Notes::
            % - The returned friction value is referred to the output of the gearbox.
            % - The friction parameters in the Link object are referred to the motor.
            % - Motor viscous friction is scaled up by G^2.
            % - Motor Coulomb friction is scaled up by G.
            % - The appropriate Coulomb friction value to use in the non-symmetric case
            %   depends on the sign of the joint velocity, not the motor velocity.

            tau = l.B * abs(l.G) * qd;

            if issym(l)
                tau = tau + l.Tc;
            elseif qd > 0
                tau = tau + l.Tc(1);
            elseif qd < 0
                tau = tau + l.Tc(2);
            end
            tau = -abs(l.G) * tau;     % friction opposes motion
        end % friction()

        function tau = friction2(l, qd)

                    % experimental code
            qdm = qd / l.G;
            taum = -l.B * qdm;
            if qdm  > 0
                taum = taum - l.Tc(1);
            elseif qdm < 0
                taum = taum - l.Tc(2);
            end
            tau = taum * l.G;
        end

        function  l2 = nofriction(l, only)
            %Link.nofriction Remove friction
            %
            % LN = L.nofriction() is a link object with the same parameters as L except
            % nonlinear (Coulomb) friction parameter is zero.
            %
            % LN = L.nofriction('all') as above except that viscous and Coulomb friction
            % are set to zero.
            %
            % LN = L.nofriction('coulomb') as above except that Coulomb friction is set to zero.
            %
            % LN = L.nofriction('viscous') as above except that viscous friction is set to zero.
            %
            % Notes::
            % - Forward dynamic simulation can be very slow with finite Coulomb friction.
            %
            % See also SerialLink.nofriction, SerialLink.fdyn.

            l2 = Link(l);

            if nargin == 1
                only = 'coulomb';
            end

            switch only
                case 'all'
                    l2.B = 0;
                    l2.Tc = [0 0];
                case 'viscous'
                    l2.B = 0;
                case 'coulomb'
                    l2.Tc = [0 0];
            end
        end

        function v = RP(l)
            %Link.RP Joint type
            %
            % c = L.RP() is a character 'R' or 'P' depending on whether joint is revolute or
            % prismatic respectively.
            % If L is a vector of Link objects return a string of characters in joint order.
            v = '';
            for ll=l
                if ll.sigma == 0
                    v = strcat(v, 'R');
                else
                    v = strcat(v, 'P');
                end
            end
        end % RP()

        function set.r(l, v)
            %Link.r Set centre of gravity
            %
            % L.r = R sets the link centre of gravity (COG) to R (3-vector).
            %
            if isempty(v)
                return;
            end
            if length(v) ~= 3
                error('RTB:Link:badarg', 'COG must be a 3-vector');
            end
            l.r = v(:).';
        end % set.r()

        function set.Tc(l, v)
            %Link.Tc Set Coulomb friction
            %
            % L.Tc = F sets Coulomb friction parameters to [F -F], for a symmetric
            % Coulomb friction model.
            %
            % L.Tc = [FP FM] sets Coulomb friction to [FP FM], for an asymmetric
            % Coulomb friction model. FP>0 and FM<0.  FP is applied for a positive
            % joint velocity and FM for a negative joint velocity.
            %
            % Notes::
            % - The friction parameters are defined as being positive for a positive
            %   joint velocity, the friction force computed by Link.friction uses the
            %   negative of the friction parameter, that is, the force opposing motion of
            %   the joint.
            %
            % See also Link.friction.
            if isempty(v)
                return;
            end

            if isa(v,'sym') && ~isempty(findsym(v))
                l.Tc = sym('Tc');
            elseif isa(v,'sym') && isempty(findsym(v))
                v = double(v);
            end

            if length(v) == 1  ~isa(v,'sym')
                l.Tc = [v -v];
            elseif length(v) == 2 && ~isa(v,'sym')
                if v(1) < v(2)
                    error('RTB:Link:badarg', 'Coulomb friction is [Tc+ Tc-]');
                end
                l.Tc = v;
            else
                error('Coulomb friction vector can have 1 (symmetric) or 2 (asymmetric) elements only')
            end
        end % set.Tc()

        function set.I(l, v)
            %Link.I Set link inertia
            %
            % L.I = [Ixx Iyy Izz] sets link inertia to a diagonal matrix.
            %
            % L.I = [Ixx Iyy Izz Ixy Iyz Ixz] sets link inertia to a symmetric matrix with
            % specified inertia and product of intertia elements.
            %
            % L.I = M set Link inertia matrix to M (3x3) which must be symmetric.
            if isempty(v)
                return;
            end
            if all(size(v) == [3 3])
                if isa(v, 'double') && norm(v-v') > eps
                    error('inertia matrix must be symmetric');
                end
                l.I = v;
            elseif length(v) == 3
                l.I = diag(v);
            elseif length(v) == 6
                l.I = [ v(1) v(4) v(6)
                    v(4) v(2) v(5)
                    v(6) v(5) v(3) ];
            else
                error('RTB:Link:badarg', 'must set I to 3-vector, 6-vector or symmetric 3x3');
            end
        end % set.I()

        function v = islimit(l, q)
            %Link.islimit  Test joint limits
            %
            % L.islimit(Q) is true (1) if Q is outside the soft limits set for this joint.
            %
            % Note::
            % - The limits are not currently used by any Toolbox functions.
            if isempty(l.qlim)
                error('no limits assigned to link')
            end
            v = (q > l.qlim(2)) - (q < l.qlim(1));
        end % islimit()

        function v = isrevolute(L)
            %Link.isrevolute  Test if joint is revolute
            %
            % L.isrevolute() is true (1) if joint is revolute.
            %
            % See also Link.isprismatic.

            v = [L.sigma] == 0;
        end

        function v = isprismatic(L)
            %Link.isprismatic  Test if joint is prismatic
            %
            % L.isprismatic() is true (1) if joint is prismatic.
            %
            % See also Link.isrevolute.
            v = [L.sigma] == 1;
        end

        function T = A(L, q)
            %Link.A Link transform matrix
            %
            % T = L.A(Q) is the link homogeneous transformation matrix (4x4) corresponding
            % to the link variable Q which is either the Denavit-Hartenberg parameter THETA (revolute)
            % or D (prismatic).
            %
            % Notes::
            % - For a revolute joint the THETA parameter of the link is ignored, and Q used instead.
            % - For a prismatic joint the D parameter of the link is ignored, and Q used instead.
            % - The link offset parameter is added to Q before computation of the transformation matrix.
            sa = sin(L.alpha); ca = cos(L.alpha);
            if L.flip
                q = -q + L.offset;
            else
                q = q + L.offset;
            end
            if L.sigma == 0
                % revolute
                st = sin(q); ct = cos(q);
                d = L.d;
            else
                % prismatic
                st = sin(L.theta); ct = cos(L.theta);
                d = q;
            end

            if L.mdh == 0
                % standard DH

                T =    [    ct  -st*ca  st*sa   L.a*ct
                    st  ct*ca   -ct*sa  L.a*st
                    0   sa      ca      d
                    0   0       0       1];
            else
                % modified DH

                T =    [    ct      -st     0   L.a
                    st*ca   ct*ca   -sa -sa*d
                    st*sa   ct*sa   ca  ca*d
                    0       0       0   1];
            end

        end % A()

        function display(l)
            %Link.display Display parameters
            %
            % L.display() displays the link parameters in compact single line format.  If L is a
            % vector of Link objects displays one line per element.
            %
            % Notes::
            % - This method is invoked implicitly at the command line when the result
            %   of an expression is a Link object and the command has no trailing
            %   semicolon.
            %
            % See also Link.char, Link.dyn, SerialLink.showlink.
            loose = strcmp( get(0, 'FormatSpacing'), 'loose');
            if loose
                disp(' ');
            end
            disp([inputname(1), ' = '])
            disp( char(l) );
        end % display()

        function s = char(links, from_robot)
            %Link.char Convert to string
            %
            % s = L.char() is a string showing link parameters in a compact single line format.
            % If L is a vector of Link objects return a string with one line per Link.
            %
            % See also Link.display.

            % display in the order theta d a alpha
            if nargin < 2
                from_robot = false;
            end

            s = '';

            for j=1:length(links)
                l = links(j);

                conv = l.RP;
                if l.mdh == 0
                    conv = strcat(conv, ',stdDH');
                else
                    conv = strcat(conv, ',modDH');
                end
                if length(links) == 1
                    qname = 'q';
                else
                    qname = sprintf('q%d', j);
                end

                if from_robot
                    % invoked from SerialLink.char method, format for table
                    if l.sigma == 1
                        % prismatic joint
                        js = sprintf('|%3d|%11s|%11s|%11s|%11s|%11s|', ...
                            j, ...
                            render(l.theta), ...
                            qname, ...
                            render(l.a), ...
                            render(l.alpha), ...
                            render(l.offset));
                    else
                        % revolute joint
                        js = sprintf('|%3d|%11s|%11s|%11s|%11s|%11s|', ...
                            j, ...
                            qname, ...
                            render(l.d), ...
                            render(l.a), ...
                            render(l.alpha), ...
                            render(l.offset));
                    end
                else

                    if l.sigma == 1
                        % prismatic joint
                        js = sprintf(' theta=%s, d=%s, a=%s, alpha=%s, offset=%s (%s)', ...
                            render(l.theta), ...
                            qname, ...
                            render(l.a), ...
                            render(l.alpha), ...
                            render(l.offset), ...
                            conv);
                    else
                        % revolute
                        js = sprintf(' theta=%s, d=%s, a=%s, alpha=%s, offset=%s (%s)', ...
                            qname, ...
                            render(l.d), ...
                            render(l.a), ...
                            render(l.alpha), ...
                            render(l.offset), ...
                            conv);
                    end
                end
                if isempty(s)
                    s = js;
                else
                    s = char(s, js);
                end
            end

        end % char()

        function dyn(links)
            %Link.dyn Show inertial properties of link
            %
            % L.dyn() displays the inertial properties of the link object in a multi-line
            % format. The properties shown are mass, centre of mass, inertia, friction,
            % gear ratio and motor properties.
            %
            % If L is a vector of Link objects show properties for each link.
            %
            % See also SerialLink.dyn.

            for j=1:numel(links)
                l = links(j);
                if numel(links) > 1
                    fprintf('\nLink %d::', j);
                end
                fprintf('%s\n', l.char());
                if ~isempty(l.m)
                    fprintf('  m    = %s\n', render(l.m))
                end
                if ~isempty(l.r)
                    s = render(l.r);
                    fprintf('  r    = %s %s %s\n', s{:});
                end
                if ~isempty(l.I)
                    s = render(l.I(1,:));
                    fprintf('  I    = | %s %s %s |\n', s{:});
                    s = render(l.I(2,:));
                    fprintf('         | %s %s %s |\n', s{:});
                    s = render(l.I(3,:));
                    fprintf('         | %s %s %s |\n', s{:});
                end
                if ~isempty(l.Jm)
                    fprintf('  Jm   = %s\n', render(l.Jm));
                end
                if ~isempty(l.B)
                    fprintf('  Bm   = %s\n', render(l.B));
                end
                if ~isempty(l.Tc)
                    fprintf('  Tc   = %s(+) %s(-)\n', ...
                        render(l.Tc(1)), render(l.Tc(2)));
                end
                if ~isempty(l.G)
                    fprintf('  G    = %s\n', render(l.G));
                end
                if ~isempty(l.qlim)
                    fprintf('  qlim = %f to %f\n', l.qlim(1), l.qlim(2));
                end
            end
        end % dyn()

        % Make a copy of a handle object.
        % http://www.mathworks.com/matlabcentral/newsreader/view_thread/257925
        function new = copy(this)

            for j=1:length(this)
                % Instantiate new object of the same class.
                %new(j) = feval(class(this(j)));
                new(j) = Link();
                % Copy all non-hidden properties.
                p = properties(this(j));
                for i = 1:length(p)
                    new(j).(p{i}) = this(j).(p{i});
                end
            end
        end

        function links = horzcat(this, varargin)
            links = this.toLink();

            for i=1:length(varargin)
                links = cat(2, links, varargin{i}.toLink());
            end
        end

        function links = vertcat(this, varargin)
            links = this.horzcat(varargin{:});
        end

        function new = toLink(this)
            for j=1:length(this)
                new(j) = Link();

                % Copy all non-hidden properties.
                p = properties(this(j));
                for i = 1:length(p)
                    new(j).(p{i}) = this(j).(p{i});
                end
            end

        end

        function res = issym(l)
            %Link.issym Check if link is a symbolic model
            %
            % res = L.issym() is true if the Link L has symbolic parameters.

            res = isa(l.alpha,'sym');
        end

        function sl = sym(l)

            sl = Link(l);   % clone the link

            if ~isempty(sl.theta)
                sl.theta = sym(sl.theta);
            end
            if ~isempty(sl.d)
                sl.d = sym(sl.d);
            end
            sl.alpha = sym(sl.alpha);
            sl.a = sym(sl.a);
            sl.offset = sym(sl.offset);

            sl.I = sym(sl.I);
            sl.r = sym(sl.r);
            sl.m = sym(sl.m);

            sl.Jm = sym(sl.Jm);
            sl.G = sym(sl.G);
            sl.B = sym(sl.B);
            sl.Tc = sym(sl.Tc);

        end
    end % methods

end % class

function s = render(v, fmt)
    if nargin < 2
        fmt = '%11.4g';
    end
    if length(v) == 1
                if isa(v, 'double')
            s = sprintf(fmt, v);
        elseif isa(v, 'sym')
            s = char(v);
        else
            error('RTB:Link:badarg', 'Link parameter must be numeric or symbolic');
        end
    else

    for i=1:length(v)
        if isa(v, 'double')
            s{i} = sprintf(fmt, v(i));
        elseif isa(v, 'sym')
            s{i} = char(v(i));
        else
            error('RTB:Link:badarg', 'Link parameter must be numeric or symbolic');
        end
    end
    end
end

author: Allopart

repo: rbpf-gmapping


text:
function gX = rbfperiodic2KernDiagGradX(kern, X)

% RBFPERIODIC2KERNDIAGGRADX Gradient of RBFPERIODIC2 kernel's diagonal with respect to X.
% FORMAT
% DESC computes the gradient of the diagonal of the RBF periodic covariance with variying period kernel matrix with
% respect to the elements of the design matrix given in X.
% ARG kern : the kernel structure for which gradients are being computed.
% ARG X : the input data in the form of a design matrix.
% RETURN gX : the gradients of the diagonal with respect to each element
% of X. The returned matrix has the same dimensions as X.
%
% SEEALSO : rbfperiodic2KernParamInit, kernDiagGradX, rbfperiodic2kernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2007, 2009
%
% MODIFICATIONS : Andreas C. Damianou, 2011
%
% MODIFICATIONS : Michalis K. Titsias, 2011

% KERN

gX = zeros(size(X));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function txt = infoTextRoi(roi, format);
%
% Print out info about an ROI file / struct
% in a human-readable set of strings.
%
% txt = infoTextRoi(roi, [format]);
%
% roi: roi data struct, or path to an roi file to load.
%
% format: a flag; if 0 [default], returns as a cell-of-strings
%         for display in a listbox uicontrol;
%         if 1, returns as a single char vector w/ control
%         characters (/n) inserted.
%
% This function could be applied to other structs besides the
% roi.info struct it was first written for. It simply takes
% all fields in the struct that are numeric or char, and prints the
% field name next to the field value.
%
% ras, 10/2005.
if notDefined('format'), format = 0; end

% (1) initialize an empty text cell
txt = {};

% (2) add fields from the struct
fields = setdiff(fieldnames(roi), {'coords' 'lineHandles'});
for f = fields(:)'
    val = roi.(f{1});
    if ischar(val)
        txt{end+1} = sprintf('%s: %s', f{1}, val);

    elseif isnumeric(val) & (size(val,1)==1 | size(val,2)==1)
        txt{end+1} = sprintf('%s: %s', f{1}, num2str(val));

    end
end

% (3) Add list of first <=100 coords:
nVoxels = size(roi.coords,2);
txt{end+1} = sprintf('# of Voxels: %i', nVoxels);
txt{end+1} = sprintf('Total Volume: %3.2f %s^3', ...
                     nVoxels * prod(roi.voxelSize), roi.dimUnits{1});
txt{end+1} = '';
txt{end+1} = sprintf('Coordinates in %s:', roi.reference);
nCoords = min(100, nVoxels);
if nCoords < nVoxels
    txt{end+1} = '(First 100 only)';
end

for n = 1:nCoords
    txt{end+1} = num2str(roi.coords(:,n)');
end

if format==1    % format as a single string
    tmp = '';
    for i = 1:length(txt)
        tmp = sprintf('%s \n%s',tmp,txt{i});
    end
    txt = tmp;
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%
% Written by M. Harper Langston - 5/10/00
% harper@cims.nyu.edu
%
% Using Lemma 10.5, form the orthogonal matrix, Q, the block matrix of
% matrices of eigenvectors of a TST block matrix. Then use this to form Gamma
% See Iselres page 246 or Lemma 10.5 for more info for Q.
% Gamma will be transformed into a plain tridiagonal matrix, using the method
% in Iserles, page 247 and 249.
function [Gam] = Form_Gamma(m,N)
for j = 1:m
   for k = 1:m
      Q(j,k) = sqrt(2/(m+1))*sin(pi*j*k/(m+1));
   end
end
% Now, form Gamma.  As per Iserles, page 247, we can rearrange columns into
% rows of the vectors, so Gamma will be tridiagonal.  The only difference is
% that in the problem, Gamma*t = c, solving for t, c must be converted and
% then t will have to be converted back.  See Iselres, pages 247 and 249
% for more informtion.
if N==9 | N==10
	S = (diag((-10/3)*ones(m,1)) + diag((2/3)*ones(m-1,1),1) + diag((2/3)*ones(m-1,1),-1));
   T = (diag((2/3)*ones(m,1)) + diag((1/6)*ones(m-1,1),1) + diag((1/6)*ones(m-1,1),-1));
elseif N==5
   S = (diag((-4)*ones(m,1)) + diag((1)*ones(m-1,1),1) + diag((1)*ones(m-1,1),-1));
   T = diag((1)*ones(m,1));
elseif N==101 % N = 101 is called by Modified_Right.m only, to  speed up
              % the modified 9=point scheme
   S = (diag((2/3)*ones(m,1)) + diag((1/12)*ones(m-1,1),1) + diag((1/12)*ones(m-1,1),-1));
   T = (diag((1/12)*ones(m,1)));
else
   error('Must either use 5 or 9-point scheme (N=5 or N=9 in poisson.m)');
end
% inv(Q) = Q, so don't need to invert Q to diagonalize S and T
GS = diag(Q*S*Q);
GT = diag(Q*T*Q);
Gam = sparse(zeros(size(Q)));
for j = 1:m
   temp = (j-1)*m +1;
   Gam(temp:j*m,temp:j*m) = diag(GS(j)*ones(m,1)) ...
	        + diag((GT(j))*ones(m-1,1),1) + diag((GT(j))*ones(m-1,1),-1);
end
Gam = sparse(Gam);
%
% Written by M. Harper Langston - 5/10/00
% harper@cims.nyu.edu
%

author: Sable

repo: mcbench-benchmarks


text:
function fh = mycursors(varargin)
% MYCURSORS - Create as many cursors as you want on an axes
%   FH = MYCURSORS(VARARGIN) add cursor on an axes.
%
%   VARARGIN can contain 0, 1 or 2 elements:
%     1) First input is the handle of the axes on which you want to create the
%        cursor.
%        Default value is the current axes (GCA)
%     2) Second input is the color of the cursor.
%        Default value is red ('r' or [1 0 0])
%
%   FH is a structure containing function handles that enable you to manipulate
%   the cursor. FH structure looks like:
%     FH.add(XPOS) : Adds a cursor on the same axes at the position x = XPOS.
%                    XPOS is facultative. Default value is the middle of the
%                    axes.
%     FH.off(NUM)  : Removes the NUMth cursor created
%                    If NUM is a character string then the function removes only
%                    the last cursor.
%     FH.val(NUM)  : Returns a vector containing the positions of the NUMth
%                    cursor.
%                    If you omit NUM then the function returns the position of
%                    all cursors.
%     FH.hdl()     : Returns a matrix containing the handles of the cursors in
%                    the first column and the handles of the associated labels
%                    in the second column.

% Default inputs
inArgs = {...
  gca, ... % Current axes is used by default
  'r'  ... % Default color for cursors is red
  };

% Replace default input arguments by input values
inArgs(1:nargin) = varargin;

% Initializations
axeH = inArgs{1};
figH = ancestor(axeH,'figure');
curH = [];
labH = [];
curC = [];
curL = [];
yLim = ylim(axeH);
curColor = inArgs{2};

% Add a cursor
CursorAdd();

% Assign Output Argument
fh = struct('add',@CursorAdd,'off',@CursorOff,'val',@CursorValue,'hdl',@CursorHandles);

% Drag Cursor shape
curCData = [...
  NaN NaN NaN NaN NaN   2   2   2   2   2 NaN NaN NaN NaN NaN NaN
  NaN NaN NaN NaN NaN   2   1   2   1   2 NaN NaN NaN NaN NaN NaN
  NaN NaN NaN NaN NaN   2   1   2   1   2 NaN NaN NaN NaN NaN NaN
  NaN NaN NaN NaN   2   2   1   2   1   2   2 NaN NaN NaN NaN NaN
  NaN NaN NaN   2   1   2   1   2   1   2   1   2 NaN NaN NaN NaN
  NaN NaN   2   1   1   2   1   2   1   2   1   1   2 NaN NaN NaN
  NaN   2   1   1   1   1   1   2   1   1   1   1   1   2 NaN NaN
    2   1   1   1   1   1   1   2   1   1   1   1   1   1   2 NaN
  NaN   2   1   1   1   1   1   2   1   1   1   1   1   2 NaN NaN
  NaN NaN   2   1   1   2   1   2   1   2   1   1   2 NaN NaN NaN
  NaN NaN NaN   2   1   2   1   2   1   2   1   2 NaN NaN NaN NaN
  NaN NaN NaN NaN   2   2   1   2   1   2   2 NaN NaN NaN NaN NaN
  NaN NaN NaN NaN NaN   2   1   2   1   2 NaN NaN NaN NaN NaN NaN
  NaN NaN NaN NaN NaN   2   1   2   1   2 NaN NaN NaN NaN NaN NaN
  NaN NaN NaN NaN NaN   2   2   2   2   2 NaN NaN NaN NaN NaN NaN
  NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
  ];
curHotSpot = [8 8];

%----------------------------------------------------------

  function CursorAdd(varargin)
    % Adds a new cursor on the axes

    % Check input to define default cursor position
    if nargin
      defPos = varargin{1};
    else
      defPos = diff(xlim(axeH)) / 2;
    end

    % Get axes limits
    yLim2 = ylim(axeH);

    % Create Cursor Line
    curH(end+1) = line([defPos defPos],yLim2,'Color',curColor,'LineWidth',2,'Parent',axeH);

    % Create associated label
    labH(end+1) = text(defPos,yLim2(2),...
      sprintf('%u',length(curH)),...
      'BackgroundColor','w',...
      'EdgeColor',curColor,...
      'LineWidth',2,...
      'Color',curColor,...
      'FontSize',8,...
      'FontWeight','bold');

    % Set the WindowButtonDownFcn
    set(curH(end),'ButtonDownFcn',@ButtonDownFcn)

  end

%----------------------------------------------------------

  function CursorOff(varargin)
    % Remove one or all cursors

    % Execute only if there is a cursor
    if ~isempty(curH)

      if nargin

        if ischar(varargin{1})
          ind2del = length(curH); % Remove last cursor if nothing is specified
        else
          ind2del = varargin{1};
          ind2del(ind2del>length(curH)) = []; % Remove indices outside the range
          ind2del(ind2del<1) = []; % Remove indices outside the range
        end

        % Delete specified cursor and associated label
        delete(curH(ind2del));
        delete(labH(ind2del));

        % Remove cursor and label handles from vectors
        curH(ind2del) = [];
        labH(ind2del) = [];

        % Update labels
        for indL=1:length(labH)
          set(labH(indL),'String',indL)
        end

      else

        % Remove all cursors and all associated labels
        delete([curH labH]);

      end

    end

  end

%----------------------------------------------------------

  function res = CursorValue(varargin)
    % Returns the positions of each cursor in a vector

    % Executes only if there is a cursor
    if ~isempty(curH)

      if nargin

        % Ensure that number specified in input argument is correct
        ind = varargin{1};
        ind(ind>length(curH)) = []; % Remove indices outside the range
        ind(ind<1) = []; % Remove indices outside the range

        % Retrieve data
        res = get(curH(ind),'XData');

      else

        % If no input specified then returns all positions
        if ~iscell(get(curH,'XData'))
          res = cell2mat({get(curH,'XData')});
        else
          res = cell2mat(get(curH,'XData'));
        end

      end

      % Remove the second column
      res(:,2) = [];

    else

      % If there is no cursor then returns an empty vector
      res = [];

    end

  end

%----------------------------------------------------------

  function res = CursorHandles()
    % Returns a matrix with cursor handles in the first column and labels
    % handles in the second

    % Execute only if there is a cursor
    if ~isempty(curH)

      res = [curH(:) labH(:)];

    end

  end

%----------------------------------------------------------

  function ButtonDownFcn(varargin)

    % Update current cursor handle
    curC = varargin{1};

    if ~isempty(curC)

      % Get label handle associated
      curL = labH(curC==curH);

      % Change cursor pointer
      set(figH,'Pointer','custom','PointerShapeCData',curCData,'PointerShapeHotSpot',curHotSpot)

      % Update EraseMode property for the current Cursor
      set(curC,'EraseMode','xor')

      % Update figure callbacks
      set(figH,'WindowButtonMotionFcn',@ButtonMotionFcn)
      set(figH,'WindowButtonUpFcn',@ButtonUpFcn)

    end

  end

%----------------------------------------------------------

  function ButtonUpFcn(varargin)

    % Change cursor pointer
    set(figH,'Pointer','arrow')

    % Update EraseMode property for the current Cursor
    set(curC,'EraseMode','normal')

    % Update Figure Callbacks
    set(figH,'WindowButtonMotionFcn','')

    % Unset current Cursor and label
    curC = [];
    curL = [];

  end

%----------------------------------------------------------

  function ButtonMotionFcn(varargin)

    % Get current point
    curP = get(axeH,'CurrentPoint');

    % Get axes limits
    yLim2 = ylim(axeH);

    % Update cursor position
    set(curC,'XData',[curP(1) curP(1)],'YData',[min(yLim(1),yLim2(1)) max(yLim(2),yLim2(2))])

    % Move label
    set(curL,'Position',[curP(1) yLim2(2)]);

    % Refresh Display
    drawnow;

  end

end


text:
%% numReplace
% Below is a demonstration of the features of the |numReplace| function

%%
clear; close all; clc;

%% REPLACING NUMBERS IN ARRAYS
%%
% An example array
A=[0,-6,3,0,0;-1,-4,-7,-9,-9;-4,-7,11,-5,12;10,-7,5,-7,13;0,11,-2,-6,2;12,4,2,-5,NaN]

%%
% Defining the input array for entries (NaN's allows) that need to be replaced
a=[2 -5 nan 0]
%%
% Defining the numbers (NaN's allows) to take their place
b=[991 992 993 994]; %Numbers to take their place

%%
% Replacing the numbers using |numReplace|
B=numReplace(A,a,b)

%% NOTES ON PERFORMANCE FOR NON-INTEGERS
% The |numReplace| function employs the |ismember| function. Hence it is
% suitable for all number cases where ismember is able to detect
% membership. Numerical precission difficulties may arise for non-integer
% entires. Consider the below:

logicMember=ismember(pi,pi+eps(pi))
logicMember=ismember(pi,pi+eps(pi)/10)

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*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 tapas_sem_plot_fits(fits, dt)
%% Plot the fits of the model.
%
% Input
%       data        -- Structure of the data of a single subject.
%       fits        -- Fits struct array
%       dt          -- Factor to normalize the area
% Output
%

% aponteeduardo@gmail.com
% copyright (C) 2019
%

n = 1;

n = n + 1;
if nargin < n
    dt = 1;
end

nconds = numel(fits);

for i = 1:nconds

    ax = subplot(nconds, 2, (i - 1) * 2 + 1);
    hold on;

    plot(fits(i).t, fits(i).pro * dt, 'k', 'linewidth', 2)

    ax = subplot(nconds, 2, (i - 1) * 2 + 2);
    hold on;
    plot(fits(i).t, fits(i).anti * dt, 'k', 'linewidth', 2)
end

end

author: translationalneuromodeling

repo: tapas


text:
function [Dimension,NodeCoord,NodeWeight,Name]=GetTSPData(infile)
if exist(infile,'file')
    fid=fopen(infile,'r');
else
    disp('Input file no exist!');
    return;
end
if fid<0
    disp('Error while open file!');
    return;
end
NodeWeight = [];
while feof(fid)==0
    temps=fgetl(fid);
    if strcmp(temps,'')
        continue;
    elseif strncmpi('NAME',temps,4)
        k=findstr(temps,':');
        Name=temps(k+1:length(temps));
    elseif strncmpi('DIMENSION',temps,9)
        k=findstr(temps,':');
        d=temps(k+1:length(temps));
        Dimension=str2double(d);
    elseif strncmpi('EDGE_WEIGHT_SECTION',temps,19)
        formatstr = [];
        for i=1:Dimension
            formatstr = [formatstr,'%g '];
        end
        NodeWeight=fscanf(fid,formatstr,[Dimension,Dimension]);
        NodeWeight=NodeWeight';
    elseif strncmpi('NODE_COORD_SECTION',temps,18) || strncmpi('DISPLAY_DATA_SECTION',temps,20)
        NodeCoord=fscanf(fid,'%g %g %g',[3 Dimension]);
        NodeCoord=NodeCoord';
    end
end
fclose(fid);


text:
% PLUGIN_ASKINSTALL - install EEGLAB plugin from a GUI or command line
%                       call.
% Usage:
%  plugin_askinstall(pluginName, pluginFunc, force);
%
% Inputs:
%  pluginName - [string] name of the plugin
%  pluginFunc - [string] function belonging to the plugin
%  force       - [boolean] force install (even if already installed)
%
% Example:
%   Force install BVA-IO
%   plugin_askinstall('bva-io', [], true);

% Copyright (C) 2012- Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function installRes = plugin_askinstall(pluginName, pluginFunc, forceInstall)

if nargin < 3, forceInstall = false; end

if nargin < 2 || ~exist(char(pluginFunc))

    if ~forceInstall
        db = dbstack;
        if length(db) > 2 && ~strcmpi(db(end).file, 'checkouteeglab.m')
            error([ 'Cannot find ' pluginName ' extension, use EEGLAB Extension Manager to install it' ]);
        end

        installRes = 0;

        % check is deactivated
        try, PLUGINLIST = evalin('base', 'PLUGINLIST'); catch, PLUGINLIST = []; end
        if ~isempty(PLUGINLIST) && isfield(PLUGINLIST, 'plugin')
            indPlugin = strmatch(lower(pluginName), lower({ PLUGINLIST.plugin }), 'exact');
            if ~isempty(indPlugin) && strcmpi(PLUGINLIST(indPlugin(1)).status, 'deactivated')
                res = questdlg2( [ pluginName ' extension is de-activated. Do you want to reactivate it now?' ], [ pluginName ' extension installation' ], 'No', 'Yes', 'Yes' );
                if strcmpi(res, 'no'), return, end
                plugin_reactivate(PLUGINLIST(indPlugin(1)).foldername);
                evalin('base', 'eeglab rebuild');
                installRes = 1;
                return;
            end
        end

        % check for installing
        res = questdlg2( [ pluginName ' extension is not installed. Do you want to download it now?' ], [ pluginName ' extension installation' ], 'No', 'Yes', 'Yes' );
    else
        res = 'yes';
    end

    if strcmpi(res, 'no'), return, end
    try
        plugins = plugin_getweb('plugin_install', []);
    catch
        error('Issue with retrieving statistics for extensions, maybe check your connection');
    end
    if isempty(plugins)
        error('Cannot download the extension. Please check your internet connection');
    end
    indPlugin = strmatch(lower(pluginName), lower({ plugins.name }), 'exact');
    if isempty(indPlugin)
        error([ pluginName ' extension not found' ]);
    end
    result = plugin_install(plugins(indPlugin(1)).zip, plugins(indPlugin(1)).name, plugins(indPlugin(1)).version, forceInstall);
    if result == 1, installRes = 1; end

    evalin('base', 'eeglab rebuild');
else
    installRes = 1;
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab