text:
% Init Script
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     08/11/2016  Royi Avital
%   *   First release.
%

%% General Parameters
close('all');
clear();
clc();

% set(0, 'DefaultFigureWindowStyle', 'docked');
% defaultLoosInset = get(0, 'DefaultAxesLooseInset');
% set(0, 'DefaultAxesLooseInset', [0.05, 0.05, 0.05, 0.05]);

figPosSmall     = [100, 100, 400, 300];
figPosMedium    = [100, 100, 800, 600];
figPosLarge     = [100, 100, 960, 720];
figPosDefault   = [100, 100, 560, 420];

fontSizeTitle   = 14;
fontSizeAxis    = 12;
fontSizeString  = 12;

lineWidthThin   = 1;
lineWidthNormal = 3;
lineWidthThick  = 4;

markerSizeSmall     = 4;
markerSizeNormal    = 8;
markerSizeLarge     = 10;

randomNumberStream  = RandStream('mlfg6331_64', 'NormalTransform', 'Ziggurat');
subStreamNumber     = round(sum(clock()));
% subStreamNumber     = 162;
% subStreamNumber     = 2122;
set(randomNumberStream, 'Substream', subStreamNumber);
RandStream.setGlobalStream(randomNumberStream);

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

author: RoyiAvital


text:
function val=soft_threshoda(w,alpha,mu)
%soft_threshoda
%
val=sign(w).*max(abs(w)-alpha*mu,0)/(1+alpha*mu);
end

author: xiubooth


text:
function sF = acos(sF, varargin)
% cost of a function
% Syntax
%   sF = acos(sF)
%   sF = acos(sF, 'bandwidth', bandwidth)
%
% Input
%  sF - @S2FunHarmonic
%
% Output
%  sF - @S2FunHarmonic
%
% Options
%  bandwidth - minimal degree of the spherical harmonic
%

sF = sF.quadrature(@(v) acos(max(-1,min(1,sF.eval(v)))),varargin{:});

end

author: mtex-toolbox

repo: mtex

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


text:
function a = atan(a)
%ATAN         Gradient inverse tangent atan(a)
%

% written  10/16/98     S.M. Rump
% modified 10/14/00     S.M. Rump  use Tony's trick
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    accelaration for sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 10/08/08     S.M. Rump  improved sparse multiplication: not using intval data type
% modified 08/26/12     S.M. Rump  global variables removed
%

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

  N = getappdata(0,'INTLAB_GRADIENT_NUMVAR');

  % use full(a.x(:)): cures Matlab V6.0 bug
  % a=7; i=[1 1]; x=a(i), b=sparse(a); y=b(i)  yields row vector x but column vector y
  % ax is full anyway
  ax = 1 ./ ( 1 + sqr(full(a.x(:))) );
  a.x = atan(a.x);
  if issparse(a.dx)
    sizeax = size(a.dx,1);
    [ia,ja,sa] = find(a.dx);
    if isa(a.x,'intval')
      adx = times(ax(ia),sa,0);
      if adx.complex
        a.dx = intval( sparse(ia,ja,adx.mid,sizeax,N) , sparse(ia,ja,adx.rad,sizeax,N) , 'midrad' );
      else
        a.dx = intval( sparse(ia,ja,adx.inf,sizeax,N) , sparse(ia,ja,adx.sup,sizeax,N) , 'infsup' );
      end
    else
      a.dx = sparse(ia,ja,ax(ia).*sa,sizeax,N);
    end
  else
    a.dx = a.dx .* ax(:,ones(1,N));
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function varargout = calcChildVariant(varargin)

warning('calcChildVariant is obsolete. Please use "calcVariantId" instead.')

[varargout{1:nargout}] = calcVariantId(varargin{:});


text:
function [AGC_atten_dB] = load_AGC_block( file_pre, first_rg_line, ...
                                        Nrg_lines_blk, block, UseMATfiles )
%
% -------------------------------------------------------------------------
% This program returns the RADARSAT-1 receiver attenuation values for
% one block.  They are integers from 2 to 17, in 1 dB steps.
% -------------------------------------------------------------------------
% file_pre        : File prefix ( path + prefix + block number)
% first_rg_line   : First range line for the extracted data
% Nrg_lines_blk   : Number of range lines in a block
% block           : the current block number
% UseMATfiles     : (0) using extract_AGC and the '_aux.dat' file
%                   (1) load 'AGC_attenuation_values.mat',
%                          which has the AGC values for every range line.
% -------------------------------------------------------------------------
% Created:  Nov 01, 2004  by Kaan Ersahin
% Modified: Nov 29, 2004  by Kaan Ersahin
% Modified: Dec 13, 2004  by Ian Cumming
%            - added UseMATfiles option
% -------------------------------------------------------------------------

if     UseMATfiles == 0  % extract AGC values from binary aux data file

    file_in = strcat(file_pre, '_aux.dat');
    tic
    AGC_atten_dB = extract_AGC( file_in, Nrg_lines_blk );
    toc

elseif UseMATfiles == 1   % load all AGC values from .mat file
    tic
    file_in = 'AGC_attenuation_values.mat'
    load( file_in )

    first = first_rg_line + (Nrg_lines_blk * (block-1));
    last  = first + Nrg_lines_blk - 1;
    AGC_atten_dB = nom_attenuation( first : last );
    toc
end

% Use check sum to verify read operation
% lenAGC = length(AGC_atten_dB);
% multfac = ([1:lenAGC] - lenAGC/2)';
% checksumAGC = sum( multfac .* AGC_atten_dB );
% fprintf('\nCheck sum of AGC values:%10.0f\n', checksumAGC )

%{
figure(201),   clf      %  Plot AGC values
plot( AGC_atten_dB ),   grid
axis([-9 length(AGC_atten_dB)+10  0.5 17.5])
xlabel('Azimuth (sample no. within this block)  \rightarrow', 'FontS', 12 )
ylabel('Magnitude  (dB)  \rightarrow', 'FontS', 12 )
title('AGC attenuation values for this block of data', 'FontS', 13 )
% text( lenAGC/2, 1.6, sprintf(...
%     '\nCheck sum of AGC values:%10.0f\n', checksumAGC ), 'Hor', 'c' )
%}
pause(0.1)

author: denkywu

repo: SAR-Synthetic-Aperture-Radar


text:
function c = horzcat(varargin)
% Implements the horzcat ([a b]) operator on image_vector objects across voxels.
% Requires that each object has an equal number of columns and voxels
%
% :Usage:
% ::
%
%    function s = horzcat(varargin)
%
% :Examples:
% ::
%
%    c = [dat1 dat2];
%
% ..
%    Programmer Notes
%    Created 3/14/14 by Luke Chang
% ..

error('This method is deprecated.  Use image_math instead.');

dat = [];
for i = 1:nargin
    %Check if image_vector object
    if ~isa(varargin{i}, 'image_vector')
        error('Input Data is not an image_vector object')
    end

    dat = [dat, varargin{i}.dat];
end

c = varargin{1};
c.dat = dat;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [GW] = MW2GW(MW)
% Convert power from megawatts to gigawatts.
% Chad A. Greene 2012
GW = MW/1000;

author: Sable

repo: mcbench-benchmarks


text:
%%
f = imread('Fig0208(a).tif');
imshow(f);
figure,imhist(f);  %  显示直方图
ylim('auto')
%%
g = histeq(f,256); %  直方图均衡
figure,imshow(g);
figure,imhist(g);
ylim('auto')

author: chengchengzi

repo: Digital-Image-Processing


text:
function tabledata = highlight_best_rows(tabledata, columns, varargin)
% highlight_best_rows Adds highlight to the best three cells in a given
% row.
%
% The function modifies the input matrix by converting it to cell array
% and augmenting the cells in individual columns with class labels for
% first three places acording to specified ordering. The resulting array
% can be then handed over to matrix2html function.
%
% Input:
% - tabledata (matrix): <description>
% - columns (cell): An array of strings that denote type of ordering for
% individual columns.
%     - Decreasing: Highlight in decreasing order.
%     - Increasing: Highlight in increasing order.
%     - None: Do not highlight.
%
% Output:
% - tabledata (cell): Modified array.
%

if size(tabledata, 2) ~= numel(columns)
    error('Number of columns does not match the sorting instructuions.');
end;

labels = {'first', 'second', 'third'};

for i = 1:numel(columns)

    switch lower(columns{i})
        case {'descending', 'descend', 'high'}
            values = cell2mat(tabledata(:, i));
            usable = find(~isnan(values));
            levels = sort(unique(values(usable)), 'descend');

        case {'ascending', 'ascend', 'low'}
            values = cell2mat(tabledata(:, i));
            usable = find(~isnan(values));
            levels = sort(unique(values(usable)), 'ascend');

        otherwise
            levels = [];
    end;

    if isempty(levels)
        continue;
    end;

    for j = 1:min(numel(labels), numel(levels))
        tabledata(values == levels(j), i) = cellfun(...
            @(x) struct('text', x, 'class', labels{j}), tabledata(values == levels(j), i), 'UniformOutput', false);

    end
end;

author: votchallenge

repo: toolkit-legacy


text:
function result = hlp_config(configname, operation, varargin)
% helper function to process human-readable config scripts.
% Result = hlp_config(FileName,Operation,VariableName,Value,NVPs...)
%
% Config scripts consist of assignments of the form name = value; to set configuration options. In
% addition, there may be any type of comments, conditional control flow, etc - e.g., setting certain
% values on some platforms and others on others. This function allows to get or set the value
% assigned to a variable in the place of the script where it is actually assigned on the current
% platform. Note that the respective variable has to be already in the config file for this function
% to work.
%
% In:
%   FileName : name of the configuration file to process
%
%   Operation : operation to perform on the config file
%               'get' : get the currently defined value of a given variable
%               'set' : replace the current defintion of a given variable
%
%   VariableName : name of the variable to be affected (must be a MATLAB identifier)
%
%   Value : the new value to be assigned, if the operation is 'set', as a string
%           note that most data structures can be converted into a string via hlp_tostring
%
%   NVPs... : list of further name-value pairs, where each name denotes a config variables and the subsequent
%             value is the string expression that should be written into the config file. It is
%             generally a good idea to use hlp_tostring() to turn a data structure into such a string
%             representation.
%
% Out:
%   Result : the current value of the variable of interest, when using the 'get'
%            operation
%
% Notes:
%   There can be multiple successive variable name / value pairs for the set mode.
%   If an error occurs during a set operation, any changes will be rolled back.
%
% Examples:
%   % read out the value of the 'data' config variable from a config file
%   data = hlp_config('/home/christian/myconfig.m','get','data')
%
%   % override the values of the 'files' and 'capacity' config variables in the given config script
%   hlp_config('/home/christian/myconfig.m', 'set', 'files',myfiles, 'capacity',1000)
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-11-19

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

if ~exist(configname,'file')
    error('hlp_config:file_not_found','The specified config file was not found.'); end

switch operation
    case 'get'
        varname = varargin{1};
        if ~isvarname(varname)
            error('hlp_config:bad_varname','The variable name must be a valid MATLAB identifier.'); end
        % get the currently defined value of a variable...
        result = get_value(configname,varname);
    case 'set'
        backupfile = [];
        try
            % apply first assignment
            backupfile = set_value(configname,varargin{1},varargin{2},true);
            for k = 4:2:length(varargin)
                % apply all other assignments
                set_value(configname,varargin{k-1},varargin{k},false); end
        catch e
            % got an error; roll back changes if necessary
            if ~isempty(backupfile)
                try
                    movefile(backupfile,configname);
                catch
                    disp(['Could not roll back changes. You can manually revert changes by replacing ' configname ' by ' backupfile '.']);
                end
            end
            rethrow(e);
        end
    otherwise
        error('hlp_config:unsupported_option','Unsupported config operation.');
end

% run the given config script and obtain the current value of the given variable...
function res = get_value(filename__,varname__)
try
    run_script(filename__);
catch e
    error('hlp_config:erroneous_file',['The config file is erroneous; Error message: ' e.message]);
end
if ~exist(varname__,'var')
    error('hlp_config:var_not_found','The variable is not being defined in the config file.'); end
res = eval(varname__);

function backup_name = set_value(filename,varname,newvalue,makebackup)
backup_name = [];
if ~exist(filename,'file')
    error('hlp_config:file_not_found','The config file was not found.'); end
if ~isvarname(varname)
    error('hlp_config:incorrect_value','The variable name must be a valid MATLAB identifier.'); end
if ~ischar(newvalue)
    error('hlp_config:incorrect_value','The value to be assigned must be given as a string.'); end
try
    % read the config file contents
    contents = {};
    f = fopen(filename,'r');
    while 1
        l = fgetl(f);
        if ~ischar(l)
            break; end
        contents{end+1} = [l 10];
    end
    fclose(f);
    % turn it into one str
    contents = [contents{:}];
catch e
    try fclose(f); catch,end
    error('hlp_config:cannot_read_config',['Cannot read the config file; Error message: ' e.message]);
end

% now check if the file is actually writable
try
    f = fopen(filename,'r+');
    if f ~= -1
        fclose(f);
    else
        error('hlp_config:permissions_error','Could not update the config file %s. Please check file permissions and try again.',filename);
    end
catch
    error('hlp_config:permissions_error','Could not update the config file %s. Please check file permissions and try again.',filename);
end

% if the desired config variables does not occur at all in the file (assuming that this is a legacy file)...
% (note that this check will not catch cases where the variable name might plausibly appear in comments or
% some sub-string)
if isempty(strfind(contents,varname))
    % ... append the assignment to the end
    f = fopen(filename,'a');
    fprintf(f,['\n%% the following variable was appended by hlp_config\n' varname ' = ' newvalue ';\n']);
    fclose(f);
    fprintf('Note: the variable %s did not occur in your config file %s; appended it to the end.\n',varname,filename);
else
    % temporarily replace stray semicolons by a special character and contract ellipses,
    % so that the subsequent assignment regex matching will not get derailed)
    evalstr = contents;
    comment_flag = false;
    string_flag = false;
    bracket_level = 0;
    ellipsis_flag = false;
    substitute = false(1,length(evalstr)); % this mask indicates where we have to subsitute reversibly by special characters
    spaceout = false(1,length(evalstr));   % this mask indicates where we can substitute irreversibly by whitespace characters...
    for k=1:length(evalstr)
        % everything that follows an ellipsis will be spaced out (including the subsequent newline that resets it)
        if ellipsis_flag
            spaceout(k) = true; end
        switch evalstr(k)
            case ';' % semicolon
                % in strs, brackets or comments: indicate need for substitution
                if string_flag || bracket_level>0 || comment_flag
                    substitute(k) = true; end
            case '''' % quotes
                % flip str flag, unless in comment
                if ~comment_flag
                    string_flag = ~string_flag; end
            case 10 % newline
                % reset bracket level, unless in ellipsis
                if ~ellipsis_flag
                    bracket_level = 0; end
                % reset comment flag, str flag and ellipsis flag
                comment_flag = false;
                string_flag = false;
                ellipsis_flag = false;
            case {'[','{'} % opening array bracket
                % if not in str nor comment, increase bracket level
                if ~string_flag && ~comment_flag
                    bracket_level = bracket_level+1; end
            case {']','}'} % closing array bracket
                % if not in str nor comment, decrease bracket level
                if ~string_flag && ~comment_flag
                    bracket_level = bracket_level-1; end
            case '%' % comment character
                % if not in str, switch on comment flag
                if ~string_flag
                    comment_flag = true; end
            case '.' % potential ellipsis character
                % if not in comment nor in str, turn on ellipsis and comment
                if ~string_flag && ~comment_flag && k>2 && strcmp(evalstr(k-2:k),'...')
                    ellipsis_flag = true;
                    comment_flag = true;
                    % we want to replace the ellipsis and everything that follows up to and including the next newline
                    spaceout(k-2:k) = true;
                end
        end
    end
    % replace the characters that need to be substituted (by the bell character)
    evalstr(substitute) = 7;
    evalstr(spaceout) = ' ';
    % replace all assignments of the form "varname = *;" by "varname{end+1} = num;"
    [starts,ends] = regexp(evalstr,[varname '\s*=[^;\n]*;']);
    for k=length(starts):-1:1
        evalstr = [evalstr(1:starts(k)-1) varname '{end+1} = struct(''assignment'',' num2str(k) ');' evalstr(ends(k)+1:end)]; end
    % add initial assignment
    evalstr = [sprintf('%s = {};\n',varname) evalstr];
    % back-substitute the special character by semicolons
    evalstr(evalstr==7) = ';';

    % evaluate contents and get the matching assignment id's
    ids = run_protected(evalstr,varname);

    % check validity of the updated value, and of the updated config file
    try
        % check if the value str can in fact be evaluated
        newvalue_eval = eval(newvalue);
    catch
        error('hlp_config:incorrect_value','The value "%s" (to be assigned to variable "%s") cannot be evaluated properly. Note that, for example, string values need to be quoted.',newvalue,varname);
    end
    % evaluate the original config script and record the full variable assignment
    [dummy,wspace_old] = run_protected(contents); %#ok<ASGLU>
    % splice the new value into the config file contents, for the last assignment in ids
    id = ids{end}.assignment;
    contents = [contents(1:starts(id)-1) varname ' = ' newvalue ';' contents(ends(id)+1:end)];
    % evaluate the new config script and record the full variable assignment
    [dummy,wspace_new] = run_protected(contents); %#ok<ASGLU>
    % make sure that the only thing that has changed is the assignment to the variable of interest
    wspace_old.(varname) = newvalue_eval;
    if ~isequalwithequalnans(wspace_old,wspace_new)
        error('hlp_config:update_failed','The config file can not be properly updated.'); end

    % apparently, everything went well, except for the following possibilities
    %  * the newly assigned value makes no sense (--> usage error)
    %  * the settings were changed for unanticipated platforms (--> this needs to be documented properly)
    if makebackup
        try
            % make a backup of the original config file using a fresh name (.bak00X)
            [p,n,x] = fileparts(filename);
            files = dir([p filesep n '*.bak*']);
            backup_numbers = cellfun(@(n)str2num(n(end-2:end)),{files.name},'UniformOutput',false);
            backup_numbers = [backup_numbers{:}];
            if ~isempty(backup_numbers)
                new_number = 1 + max(backup_numbers);
            else
                new_number = 1;
            end
            backup_name = [p filesep n '.bak' sprintf('%03i',new_number)];
            copyfile(filename,backup_name);
            % set read permissions
            warning off MATLAB:FILEATTRIB:SyntaxWarning
            fileattrib(backup_name,'+w','a');
        catch
            error('hlp_config:permissions_error','Could not create a backup of the original config file %s. Please check file permissions and try again.',filename);
        end
    end

    % split the contents into lines again
    contents = strsplit(contents,10);
    try
        % re-create the file, line by line
        f = fopen(filename,'w+');
        for k=1:length(contents)
            fwrite(f,contents{k});
            fprintf(f,'\n');
        end
        fclose(f);
        % set file attributes
        warning off MATLAB:FILEATTRIB:SyntaxWarning
        fileattrib(filename,'+w','a');
    catch
        try fclose(f); catch,end
        error('hlp_config:permissions_error','Could not override the config file %s. Please check file permissions and try again.',filename);
    end
end

% run the given config script and obtain the current value of the given variable...
function [res,wspace] = run_protected(code__,varname__)
try
    eval(code__);
    % collect all variables into a workspace struct
    infos = whos();
    for n = {infos.name}
        if ~any(strcmp(n{1},{'code__','varname__'}))
            wspace.(n{1}) = eval(n{1}); end
    end
    if exist('varname__','var')
        % if a specific variable was to be inspected...
        res = eval(varname__);
        if ~iscell(res) || length(res) < 1 || ~all(cellfun('isclass',res,'struct')) || ~all(cellfun(@(x)isfield(x,'assignment'),res))
            error('Not all assignments to the variable were correctly identified.'); end
    else
        res = [];
    end
catch e
    error('hlp_config:update_error',['The config file could not be parsed (probably it is ill-formed); Debug message: ' e.message]);
end

% split a string without fusing delimiters (unlike hlp_split)
function strs = strsplit(str, delim)
idx = strfind(str, delim);
strs = cell(numel(idx)+1, 1);
idx = [0 idx numel(str)+1];
for k = 2:numel(idx)
    strs{k-1} = str(idx(k-1)+1:idx(k)-1); end

% for old MATLABs that can't properly move files...
function movefile(src,dst)
try
    builtin('movefile',src,dst);
catch e
    if any([src dst]=='$') && hlp_matlab_version <= 705
        if ispc
            [errcode,text] = system(sprintf('move ''%s'' ''%s''',src,dst)); %#ok<NASGU>
        else
            [errcode,text] = system(sprintf('mv ''%s'' ''%s''',src,dst)); %#ok<NASGU>
        end
        if errcode
            error('Failed to move %s to %s.',src,dst); end
    else
        rethrow(e);
    end
end

repo: REST


text:
function blkStruct = slblocks

%
% Name of the subsystem which will show up in the Simulink Blocksets
% and Toolboxes subsystem.
%
blkStruct.Name = ['Self-Tuning Controllers Simulink Library - Standard'];

%
% The function that will be called when the user double-clicks on
% this icon.
%
blkStruct.OpenFcn = '';

%
% The argument to be set as the Mask Display for the subsystem.  You
% may comment this line out if no specific mask is desired.
% Example:  blkStruct.MaskDisplay = 'plot([0:2*pi],sin([0:2*pi]));';
% No display for Simulink Extras.
%
blkStruct.MaskDisplay = '';

%
% Define the Browser structure array, the first element contains the
% information for the Simulink block library and the second for the
% Simulink Extras block library.
%
Browser(1).Library = 'stcsl_std';
Browser(1).Name    = 'STCSL standard';
Browser(1).IsFlat  = 0;% Is this library "flat" (i.e. no subsystems)?

blkStruct.Browser = Browser;

% End of slblocks


text:
function movie = close (movie)
% CLOSE (avi class name) closes the class deletes all values
% NOTE: You must set the class equal to this function

fclose (movie.fid);
clear movie;


function [color, marker, line_style] = color_spec()
% -------------------------------------------------------------------------
%   Description:
%       return color information to plot figures in our paper
%
%   Citation:
%       A Comparative Study for Single Image Blind Deblurring
%       Wei-Sheng Lai, Jia-Bin Huang, Zhe Hu, Narendra Ahuja, and Ming-Hsuan Yang
%       IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2016
%
%   Contact:
%       Wei-Sheng Lai
%       wlai24@ucmerced.edu
%       University of California, Merced
% -------------------------------------------------------------------------

color = {};                             marker = {};            line_style = {};
color{end+1} = [1, 0.5, 39/255];        marker{end+1} = 'o';    line_style{end+1} = '--';
color{end+1} = [0, 1, 0];               marker{end+1} = 'o';    line_style{end+1} = '--';
color{end+1} = [136/255, 0, 21/255];    marker{end+1} = 'o';    line_style{end+1} = '--';
color{end+1} = [1, 0, 1];               marker{end+1} = 'o';    line_style{end+1} = '--';
color{end+1} = [0, 0, 1];               marker{end+1} = 'o';    line_style{end+1} = '--';
color{end+1} = [1, 0, 0];               marker{end+1} = 'o';    line_style{end+1} = '--';
color{end+1} = [0, 162/255, 232/255];   marker{end+1} = 'o';    line_style{end+1} = '--';

color{end+1} = [1, 0.5, 39/255];        marker{end+1} = 's';    line_style{end+1} = '-';
color{end+1} = [0, 1, 0];               marker{end+1} = 's';    line_style{end+1} = '-';
color{end+1} = [136/255, 0, 21/255];    marker{end+1} = 's';    line_style{end+1} = '-';
color{end+1} = [1, 0, 1];               marker{end+1} = 's';    line_style{end+1} = '-';
color{end+1} = [0, 0, 1];               marker{end+1} = 's';    line_style{end+1} = '-';
color{end+1} = [1, 0, 0];               marker{end+1} = 's';    line_style{end+1} = '-';
color{end+1} = [0, 162/255, 232/255];   marker{end+1} = 's';    line_style{end+1} = '-';


text:
%SERIALLINK.RNE_MDH Compute inverse dynamics via recursive Newton-Euler formulation
%
% Recursive Newton-Euler for modified Denavit-Hartenberg notation.  Is invoked by
% R.RNE().
%
% See also SERIALLINK.RNE.

% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, 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

function tau = rne_mdh(robot, a1, a2, a3, a4, a5)

	z0 = [0;0;1];
	grav = robot.gravity;	% default gravity from the object
	fext = zeros(6, 1);

	% Set debug to:
	%	0 no messages
	%	1 display results of forward and backward recursions
	%	2 display print R and p*
	debug = 0;

	n = robot.n;
	if numcols(a1) == 3*n,
		Q = a1(:,1:n);
		Qd = a1(:,n+1:2*n);
		Qdd = a1(:,2*n+1:3*n);
		np = numrows(Q);
		if nargin >= 3,
			grav = a2(:);
		end
		if nargin == 4,
			fext = a3;
		end
	else
		np = numrows(a1);
		Q = a1;
		Qd = a2;
		Qdd = a3;
		if numcols(a1) ~= n | numcols(Qd) ~= n | numcols(Qdd) ~= n | ...
			numrows(Qd) ~= np | numrows(Qdd) ~= np,
			error('bad data');
		end
		if nargin >= 5,
			grav = a4(:);
		end
		if nargin == 6,
			fext = a5;
		end
	end

	tau = zeros(np,n);

	for p=1:np,
		q = Q(p,:)';
		qd = Qd(p,:)';
		qdd = Qdd(p,:)';

		Fm = [];
		Nm = [];
		pstarm = [];
		Rm = [];
		w = zeros(3,1);
		wd = zeros(3,1);
		v = zeros(3,1);
		vd = grav(:);

	%
	% init some variables, compute the link rotation matrices
	%
		for j=1:n,
			link = robot.link{j};
			Tj = link(q(j));
			if link.RP == 'R',
				D = link.D;
			else
				D = q(j);
			end
			alpha = link.alpha;
			pm = [link.A; -D*sin(alpha); D*cos(alpha)];	% (i-1) P i
			if j == 1,
				pm = t2r(robot.base) * pm;
				Tj = robot.base * Tj;
			end
			Pm(:,j) = pm;
			Rm{j} = t2r(Tj);
			if debug>1,
				Rm{j}
				Pm(:,j)'
			end
		end

	%
	%  the forward recursion
	%
		for j=1:n,
			link = robot.link{j};

			R = Rm{j}';	% transpose!!
			P = Pm(:,j);
			Pc = link.r;

			%
			% trailing underscore means new value
			%
			if link.RP == 'R',
				% revolute axis
				w_ = R*w + z0*qd(j);
				wd_ = R*wd + cross(R*w,z0*qd(j)) + z0*qdd(j);
				%v = cross(w,P) + R*v;
				vd_ = R * (cross(wd,P) + ...
					cross(w, cross(w,P)) + vd);

			else
				% prismatic axis
				w_ = R*w;
				wd_ = R*wd;
				%v = R*(z0*qd(j) + v) + cross(w,P);
				vd_ = R*(cross(wd,P) + ...
					cross(w, cross(w,P)) + vd ...
				      ) + 2*cross(R*w,z0*qd(j)) + z0*qdd(j);
			end
			% update variables
			w = w_;
			wd = wd_;
			vd = vd_;

			vdC = cross(wd,Pc) + ...
				cross(w,cross(w,Pc)) + vd;
			F = link.m*vdC;
			N = link.I*wd + cross(w,link.I*w);
			Fm = [Fm F];
			Nm = [Nm N];
			if debug,
				fprintf('w: '); fprintf('%.3f ', w)
				fprintf('\nwd: '); fprintf('%.3f ', wd)
				fprintf('\nvd: '); fprintf('%.3f ', vd)
				fprintf('\nvdbar: '); fprintf('%.3f ', vdC)
				fprintf('\n');
			end
		end

	%
	%  the backward recursion
	%

		fext = fext(:);
		f = fext(1:3);		% force/moments on end of arm
		nn = fext(4:6);

		for j=n:-1:1,

			%
			% order of these statements is important, since both
			% nn and f are functions of previous f.
			%
			link = robot.link{j};

			if j == n,
				R = eye(3,3);
				P = [0;0;0];
			else
				R = Rm{j+1};
				P = Pm(:,j+1);		% i/P/(i+1)
			end
			Pc = link.r;

			f_ = R*f + Fm(:,j);
			nn_ = Nm(:,j) + R*nn + cross(Pc,Fm(:,j)) + ...
				cross(P,R*f);

			f = f_;
			nn = nn_;

			if debug,
				fprintf('f: '); fprintf('%.3f ', f)
				fprintf('\nn: '); fprintf('%.3f ', nn)
				fprintf('\n');
			end
			if link.RP == 'R',
				% revolute
				tau(p,j) = nn'*z0 + ...
					link.G^2 * link.Jm*qdd(j) + ...
					abs(link.G) * friction(link, qd(j));
			else
				% prismatic
				tau(p,j) = f'*z0 + ...
					link.G^2 * link.Jm*qdd(j) + ...
					abs(link.G) * friction(link, qd(j));
			end
		end
	end

author: petercorke

repo: robotics-toolbox-matlab


text:
function y = vl_nnbilinearpool(x, varargin)
% VL_NNBILINEARPOOL computes self outer product of a feature x and pool the features across
% all locations
%
% Author: Subhransu Maji, Aruni RoyChowdhury, Tsung-Yu Lin

%
% This file is part of the BCNN and is made available under
% the terms of the BSD license (see the COPYING file).

% input:
% forward pass:
% x: input featre of size [hight, width, channels, batches]
% backward pass:
% x: input featre of size [hight, width, channels, batches]
% dzdy: the gradient with respect to output y.

% output:
% forward pass:
% y: self outer product of x.Features are pooled across locations
%    The output size if [1, 1, channels*channels, batches]
% backward pass:
% y: graident with respect to x. y will have the same size as input x.

backMode = numel(varargin) > 0 && ~isstr(varargin{1}) ;
if backMode
  dzdy = varargin{1};
end

gpuMode = isa(x, 'gpuArray');
[h, w, ch, bs] = size(x);

if backMode
    if gpuMode
        y = gpuArray(zeros(size(x), 'single'));
    else
        y = zeros(size(x), 'single');
    end
    for b=1:bs
        dzdy_b = reshape(dzdy(1,1,:,b), [ch, ch]);
        a = reshape(x(:,:,:,b), [h*w, ch]);
        y(:, :, :, b) = reshape(a*dzdy_b, [h, w, ch])/(h*w);
    end
else
    if gpuMode
        y = gpuArray(zeros([1, 1, ch*ch, bs], 'single'));
    else
        y = zeros([1, 1, ch*ch, bs], 'single');
    end
    for b = 1:bs,
        a = reshape(x(:,:,:,b), [h*w, ch]);
        y(1,1,:, b) = reshape(a'*a, [1 ch*ch])/(h*w);
    end
end

repo: DBCNN


text:
function [val, status] = ft_findcfg(cfg, var)

% FT_FINDCFG searches for an element in the cfg structure
% or in the nested previous cfgs
%
% Use as
%   val = ft_findcfg(cfg, var)
% where the name of the variable should be specified as string.
%
% e.g.
%   trl   = ft_findcfg(cfg, 'trl')
%   event = ft_findcfg(cfg, 'event')
%
% See also FT_GETOPT, FT_CFG2KEYVAL

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

% if var(1)~='.'
%   var = ['.' var];
% end
val    = [];
depth  = 0;
status = 0;

while ~status
  depth = depth + 1;
  if ~isempty(cfg)
    if issubfield(cfg,  var)
      val = getsubfield(cfg, var);
      status = 1;
    elseif issubfield(cfg, 'previous')
      [val, status] = ft_findcfg(cfg.previous, var);
      if status, break; end
    elseif iscell(cfg)
      for i=1:length(cfg)
        [val, status] = ft_findcfg(cfg{i}, var);
        if status, break; end
      end
    else
      status = -1;
      break
    end
  else
    status = -1;
    break
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function train_spatial_diff_part_neighbour(expidx,cidx)

fprintf('train_spatial_diff_part_app_dx_dy()\n');

if (ischar(expidx))
    expidx = str2num(expidx);
end

if (ischar(cidx))
    cidx = str2num(cidx);
end

p = exp_params(expidx);

pairwiseDir = p.pairwiseDir;
if (~exist(pairwiseDir,'dir')),mkdir(pairwiseDir);end

if isfield(p, 'bidirect') && p.bidirect
    next_joint = neighbour_joint_list();
    if cidx > length(next_joint)
        return;
    end
    start_joint = cidx;
    end_joint = next_joint(cidx);
    cidxs = sort([start_joint end_joint]);
end

if isfield(p, 'allpairs') && p.allpairs
    pwIdxsAllrel = build_joint_pairs(p);
    cidxs = pwIdxsAllrel{cidx};
end

modelFname = [pairwiseDir '/spatial_model_cidx_' num2str(cidxs(1)) '_' num2str(cidxs(2))];
try
    %assert(false);
    load(modelFname,'spatial_model');
    fprintf('cidx: %d - %d\n',cidxs);
    fprintf('spatial model file loaded. quitting.\n');
catch
    [X_pos, keys_pos, boxes_pos, X_neg, keys_neg, boxes_neg] = get_spatial_features_neighbour(expidx,cidx);

    opts.X_pos_mean = mean(X_pos);
    opts.X_pos_std = std(X_pos);
    %idxs1 = X_pos(:,1) >= opts.X_pos_mean(:,1) - 3*opts.X_pos_std(:,1) & X_pos(:,1) <= opts.X_pos_mean(:,1) + 3*opts.X_pos_std(:,1);
    %idxs2 = X_pos(:,2) >= opts.X_pos_mean(:,2) - 3*opts.X_pos_std(:,2) & X_pos(:,2) <= opts.X_pos_mean(:,2) + 3*opts.X_pos_std(:,2);
    %idxs3 = X_pos(:,3) >= opts.X_pos_mean(:,3) - 3*opts.X_pos_std(:,3) & X_pos(:,3) <= opts.X_pos_mean(:,3) + 3*opts.X_pos_std(:,3);
    %X_pos = X_pos(idxs1 & idxs2 & idxs3,:);
    if isfield(p, 'neighbor_locref') && p.neighbor_locref
        X_pos = get_augm_spatial_features_diff_neighbour_locref(X_pos); %opts.X_pos_mean, p
        X_neg = get_augm_spatial_features_diff_neighbour_locref(X_neg); %opts.X_pos_mean, p
    else
        X_pos = get_augm_spatial_features_diff_neighbour(X_pos, p); %opts.X_pos_mean, p
        X_neg = get_augm_spatial_features_diff_neighbour(X_neg, p); %opts.X_pos_mean, p
    end

    if (size(X_neg,1) > size(X_pos,1))
        nFeat = size(X_pos,1);
        idxs_rnd = randperm(size(X_neg,1));
        idxsSamp = idxs_rnd(1:nFeat);
        X_neg = X_neg(idxsSamp,:);
    end

    [X_norm, opts.X_min, opts.X_max] = getFeatNorm([X_pos;X_neg]);
    X_pos_norm = X_norm(1:size(X_pos,1),:);
    X_neg_norm = X_norm(size(X_pos,1)+1:end,:);
    %X_pos_norm = X_pos;
    %X_neg_norm = X_neg;

    reg_type = 0; % L2
    C = 1e-3;

    nPos = size(X_pos_norm,1);
    nNeg = size(X_neg_norm,1);
    lab_pos = ones(nPos,1);
    lab_neg = zeros(nNeg,1);
    lab = [lab_pos; lab_neg];
    ex = sparse(double([X_pos_norm; X_neg_norm]));
    model = train(lab, ex, ['-s ' num2str(reg_type) ' -B 1 -c ' num2str(C)]);

    spatial_model.training_opts = opts;
    spatial_model.log_reg = model;
    save(modelFname, 'spatial_model');

    if (1)
        visDir = [pairwiseDir '/vis/'];
        if (~exist(visDir,'dir')),mkdir(visDir);end
        [~,acc,pred] = predict(lab, ex, model, '-b 1');
        scrsz = get(0,'ScreenSize');
        figure('Position',[1 scrsz(4) scrsz(3)/2 scrsz(4)]);
        vis_logreg(pred,acc,1:nPos,1+nPos:nPos+nNeg);
        print(gcf,'-dpng',[visDir '/logreg_cidx_' num2str(cidxs(1)) '_' num2str(cidxs(2)) '.png']);
        close all;
    end
end
% ------------------------------------------------------------------------

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [grains,grainId,mis2mean] = calcGrains(ebsd,varargin)
% grains reconstruction from 2d EBSD data
%
% Syntax
%
%   [grains, ebsd.grainId] = calcGrains(ebsd,'angle',10*degree)
%
%   % reconstruction low and high angle grain boundaries
%   lagb = 2*degree;
%   hagb = 10*degree;
%   grains = calcGrains(ebsd,'angle',[hagb lagb])
%
%   % allow grains to grow into not indexed regions
%   grains = calcGrains(ebsd('indexed'),'angle',10*degree)
%
%   % do not allow grains to grow into not indexed regions
%   grains = calcGrains(ebsd,'unitCell')
%
%   % follow non convex outer boundary
%   grains = calcGrains(ebsd,'boundary','tight')
%
%   % specify phase dependent thresholds
%   % thresholds follow the same order as ebsd.CSList and should have the same length
%   grains = calcGrains(ebsd,'angle',{angl_1 angle_2 angle_3})
%
%   % markovian clustering algorithm
%   p = 1.5;    % inflation power (default = 1.4)
%   maxIt = 10; % number of iterations (default = 4)
%   delta = 5*degree % variance of the threshold angle
%   grains = calcGrains(ebsd,'method','mcl',[p maxIt],'soft',[angle delta])
%
% Input
%  ebsd   - @EBSD
%
% Output
%  grains       - @grain2d
%  ebsd.grainId - grainId of each pixel
%
% Options
%  threshold, angle - array of threshold angles per phase of mis/disorientation in radians
%  boundary         - bounds the spatial domain ('convexhull', 'tight')
%  maxDist          - maximum distance to for two pixels to be in one grain (default inf)
%  fmc       - fast multiscale clustering method
%  mcl       - markovian clustering algorithm
%  custom    - use a custom property for grain separation
%
% Flags
%  unitCell - omit voronoi decomposition and treat a unitcell lattice
%  qhull    - use qHull for the voronin decomposition
%
% References
%
% * F.Bachmann, R. Hielscher, H. Schaeben, Grain detection from 2d and 3d
% EBSD data - Specification of the MTEX algorithm: Ultramicroscopy, 111,
% 1720-1733, 2011
%
% * C. McMahon, B. Soe, A. Loeb, A. Vemulkar, M. Ferry, L. Bassman,
%   Boundary identification in EBSD data with a generalization of fast
%   multiscale clustering, <https://doi.org/10.1016/j.ultramic.2013.04.009
%   Ultramicroscopy, 2013, 133:16-25>.
%
% See also
% GrainReconstruction GrainReconstructionAdvanced

% subdivide the domain into cells according to the measurement locations,
% i.e. by Voronoi teselation or unit cell
[V,F,I_FD] = spatialDecomposition([ebsd.prop.x(:), ebsd.prop.y(:)],ebsd.unitCell,varargin{:});
% V - list of vertices
% F - list of faces
% D - cell array of cells
% I_FD - incidence matrix faces to vertices

% determine which cells to connect
[A_Db,I_DG] = doSegmentation(I_FD,ebsd,varargin{:});
% A_db - neigbhouring cells with (inner) grain boundary
% I_DG - incidence matrix cells to grains

% compute grain ids
[grainId,~] = find(I_DG.');

% phaseId of each grain
phaseId = full(max(I_DG' * ...
  spdiags(ebsd.phaseId,0,numel(ebsd.phaseId),numel(ebsd.phaseId)),[],2));
phaseId(phaseId==0) = 1; % why this is needed?

% compute boundary this gives
% I_FDext - faces x cells external grain boundaries
% I_FDint - faces x cells internal grain boundaries
[I_FDext, I_FDint, Fext, Fint] = calcBoundary;

if check_option(varargin,'removeQuadruplePoints')
  qAdded = removeQuadruplePoints;
end

% setup grains
grains = grain2d( makeBoundary(Fext,I_FDext), ...
  calcPolygons(I_FDext * I_DG,Fext,V), ...
  [], ebsd.CSList, phaseId, ebsd.phaseMap, varargin{:});

grains.grainSize = full(sum(I_DG,1)).';
grains.innerBoundary = makeBoundary(Fint,I_FDint);
grains.scanUnit = ebsd.scanUnit;

% merge quadruple grains
if check_option(varargin,'removeQuadruplePoints') && qAdded > 0
  mergeQuadrupleGrains;
end

% calc mean orientations, GOS and mis2mean
% ----------------------------------------

[d,g] = find(I_DG);

grainRange    = [0;cumsum(grains.grainSize)];        %
firstD        = d(grainRange(2:end));
q             = quaternion(ebsd.rotations);
meanRotation  = q(firstD);
GOS = zeros(length(grains),1);

% choose between equivalent orientations in one grain such that all are
% close together
for pId = grains.indexedPhasesId
  ndx = ebsd.phaseId(d) == pId;
  if ~any(ndx), continue; end
  q(d(ndx)) = project2FundamentalRegion(q(d(ndx)),ebsd.CSList{pId},meanRotation(g(ndx)));
end

% TODO: this can be done more efficiently using accumarray

% compute mean orientation and GOS
doMeanCalc = find(grains.grainSize>1 & grains.isIndexed);
for k = 1:numel(doMeanCalc)
  qind = subSet(q,d(grainRange(doMeanCalc(k))+1:grainRange(doMeanCalc(k)+1)));
  mq = mean(qind,'robust');
  meanRotation = setSubSet(meanRotation,doMeanCalc(k),mq);
  GOS(doMeanCalc(k)) = mean(angle(mq,qind));
end

% save
grains.prop.GOS = GOS;
grains.prop.meanRotation = reshape(meanRotation,[],1);
mis2mean = inv(rotation(q(:))) .* grains.prop.meanRotation(grainId(:));

% assign variant and parent Ids for variant-based grain computation
if check_option(varargin,'variants')
    variantId = get_option(varargin,'variants');
    grains.prop.variantId = variantId(firstD,1);
    grains.prop.parentId = variantId(firstD,2);
end

  function [A_Db,I_DG] = doSegmentation(I_FD,ebsd,varargin)
    % segmentation
    %
    %
    % Output
    %  A_Db - adjecency matrix of grain boundaries
    %  A_Do - adjecency matrix inside grain connections

    % extract segmentation method
    grainBoundaryCiterions = dir([mtex_path '/EBSDAnalysis/@EBSD/private/gbc*.m']);
    grainBoundaryCiterions = {grainBoundaryCiterions.name};
    gbcFlags = regexprep(grainBoundaryCiterions,'gbc_(\w*)\.m','$1');

    gbc      = get_flag(varargin,gbcFlags,'angle');
    gbcValue = ensurecell(get_option(varargin,{gbc,'threshold','delta'},15*degree,{'double','cell'}));

    if numel(gbcValue) == 1 && length(ebsd.CSList) > 1
      gbcValue = repmat(gbcValue,size(ebsd.CSList));
    end

    % get pairs of neighbouring cells {D_l,D_r} in A_D
    A_D = I_FD'*I_FD==1;
    [Dl,Dr] = find(triu(A_D,1));

    if check_option(varargin,'maxDist')
      xyDist = sqrt((ebsd.prop.x(Dl)-ebsd.prop.x(Dr)).^2 + ...
        (ebsd.prop.y(Dl)-ebsd.prop.y(Dr)).^2);

      dx = sqrt(sum((max(ebsd.unitCell)-min(ebsd.unitCell)).^2));
      maxDist = get_option(varargin,'maxDist',3*dx);
      % maxDist = get_option(varargin,'maxDist',inf);
    else
      maxDist = 0;
    end

    connect = zeros(size(Dl));

    for p = 1:numel(ebsd.phaseMap)

      % neighboured cells Dl and Dr have the same phase
      if maxDist > 0
        ndx = ebsd.phaseId(Dl) == p & ebsd.phaseId(Dr) == p & xyDist < maxDist;
      else
        ndx = ebsd.phaseId(Dl) == p & ebsd.phaseId(Dr) == p;
      end

      connect(ndx) = true;

      % check, whether they are indexed
      ndx = ndx & ebsd.isIndexed(Dl) & ebsd.isIndexed(Dr);

      % now check for the grain boundary criterion
      if any(ndx)

        connect(ndx) = feval(['gbc_' gbc],...
          ebsd.rotations,ebsd.CSList{p},Dl(ndx),Dr(ndx),gbcValue{p},varargin{:});

      end
    end

    % adjacency of cells that have no common boundary
    ind = connect>0;
    A_Do = sparse(double(Dl(ind)),double(Dr(ind)),connect(ind),length(ebsd),length(ebsd));
    if check_option(varargin,'mcl')

      param = get_option(varargin,'mcl');
      if isempty(param), param = 1.4; end
      if length(param) == 1, param = [param,4]; end

      A_Do = mclComponents(A_Do,param(1),param(2));
      A_Db = sparse(double(Dl),double(Dr),true,length(ebsd),length(ebsd));
      A_Db(A_Do~=0) = false;

    else

      A_Db = sparse(double(Dl(connect<1)),double(Dr(connect<1)),true,...
        length(ebsd),length(ebsd));

    end
    A_Do = A_Do | A_Do.';

    % adjacency of cells that have a common boundary
    A_Db = A_Db | A_Db.';

    % compute I_DG connected components of A_Do
    % I_DG - incidence matrix cells to grains
    I_DG = sparse(1:length(ebsd),double(connectedComponents(A_Do)),1);

  end

  function [I_FDext, I_FDint, Fext, Fint] = calcBoundary
    % distinguish between interior and exterior grain boundaries

    % cells that have a subgrain boundary, i.e. a boundary with a cell
    % belonging to the same grain
    sub = ((A_Db * I_DG) & I_DG)';                 % grains x cell
    [i,j] = find( diag(any(sub,1))*double(A_Db) ); % all adjacence to those
    sub = any(sub(:,i) & sub(:,j),1);              % pairs in a grain

    % split grain boundaries A_Db into interior and exterior
    A_Db_int = sparse(i(sub),j(sub),1,size(I_DG,1),size(I_DG,1));
    A_Db_ext = A_Db - A_Db_int;                    % adjacent over grain boundray

    % create incidence graphs
    I_FDbg = diag( sum(I_FD,2)==1 ) * I_FD;
    D_Fbg  = diag(any(I_FDbg,2));

    [ix,iy] = find(A_Db_ext);
    D_Fext  = diag(sum(abs(I_FD(:,ix)) & abs(I_FD(:,iy)),2)>0);

    I_FDext = (D_Fext| D_Fbg)*I_FD;

    [ix,iy] = find(A_Db_int);
    D_Fsub  = diag(sum(abs(I_FD(:,ix)) & abs(I_FD(:,iy)),2)>0);
    I_FDint = D_Fsub*I_FD;

    % remove empty lines from I_FD, F, and V
    isExt = full(any(I_FDext,2));
    I_FDext = I_FDext.'; I_FDext = I_FDext(:,isExt).';

    isInt = full(any(I_FDint,2));
    I_FDint = I_FDint.'; I_FDint = I_FDint(:,isInt).';

    % remove vertices that are not needed anymore
    [inUse,~,F] = unique(F(isExt | isInt,:));
    V = V(inUse,:);
    F = reshape(F,[],2);
    Fext = F(isExt(isExt | isInt),:); % external boundary segments
    Fint = F(isInt(isExt | isInt),:); % internal boundary segments

  end

  function gB = makeBoundary(F,I_FD)

    % compute ebsdInd
    [eId,fId] = find(I_FD.');
    eId = eId(:); fId = fId(:);

    % replace fId that appears a second time by fId + length(F)+1
    % such that it refers to the second column
    d = diff([0;fId]);
    fId = cumsum(d>0) + (d==0)*size(F,1);

    %  ebsdInd - [Id1,Id2] list of adjecent EBSD pixels for each segment
    ebsdInd = zeros(size(F,1),2);
    ebsdInd(fId) = eId;

    % compute misrotations
    mori = rotation.nan(size(F,1),1);
    isNotBoundary = all(ebsdInd,2);
    mori(isNotBoundary) = ...
      inv(ebsd.rotations(ebsdInd(isNotBoundary,2))) ...
      .* ebsd.rotations(ebsdInd(isNotBoundary,1));

    gB = grainBoundary(V,F,ebsdInd,grainId,ebsd.phaseId,mori,ebsd.CSList,ebsd.phaseMap,ebsd.id);

  end

  function qAdded = removeQuadruplePoints

    quadPoints = find(accumarray(reshape(Fext(full(any(I_FDext,2)),:),[],1),1) == 4);
    qAdded = 0;

    if isempty(quadPoints), return; end

    % find the 4 edges connected to the quadpoints
    I_FV = sparse(repmat((1:size(Fext,1)).',1,2),Fext,ones(size(Fext)));

    quadPoints = find(sum(I_FV) == 4).';
    [iqF,~] = find(I_FV(:,quadPoints));

    % this is a length(quadPoints x 4 list of edges
    iqF = reshape(iqF,4,length(quadPoints)).';

    % find the 4 vertices adfacent to each quadruple point
    qV = [Fext(iqF.',1).';Fext(iqF.',2).'];
    qV = qV(qV ~= reshape(repmat(quadPoints.',8,1),2,[]));
    qV = reshape(qV,4,[]).';

    % compute angle with respect to quadruple point
    qOmega = reshape(atan2(V(qV,1) - V(repmat(quadPoints,1,4),1),...
      V(qV,2) - V(repmat(quadPoints,1,4),2)),[],4);

    % sort the angles
    [~,qOrder] = sort(qOmega,2);

    % find common pixels for pairs of edges - first we try 1/4 and 2/3
    s = size(iqF);
    orderSub = @(i) sub2ind(s,(1:s(1)).',qOrder(:,i));

    iqD = I_FDext(iqF(orderSub(1)),:) .* I_FDext(iqF(orderSub(4)),:) + ...
      I_FDext(iqF(orderSub(2)),:) .* I_FDext(iqF(orderSub(3)),:);

    % if not both have one common pixel
    switchOrder = full(sum(iqD,2))~= 2;

    % switch to 3/4 and 1/2
    qOrder(switchOrder,:) = qOrder(switchOrder,[4 1 2 3]);
    orderSub = @(i) sub2ind(s,(1:s(1)).',qOrder(:,i));

    iqD = I_FDext(iqF(orderSub(1)),:) .* I_FDext(iqF(orderSub(4)),:) + ...
      I_FDext(iqF(orderSub(2)),:) .* I_FDext(iqF(orderSub(3)),:);

    % some we will not be able to remove
    ignore = full(sum(iqD,2)) ~= 2;
    iqD(ignore,:) = [];
    quadPoints(ignore) = [];
    iqF(ignore,:) = [];
    qV(ignore,:) = [];
    qOrder(ignore,:) = [];
    s = size(iqF);
    orderSub = @(i) sub2ind(s,(1:s(1)).',qOrder(:,i));

    % add an additional vertex (with the same coordinates) for each quad point
    newVid = (size(V,1) + (1:length(quadPoints))).';
    V = [V;V(quadPoints,:)];

    % include new vertex into face list, i.e. replace quadpoint -> newVid
    Ftmp = Fext(iqF(orderSub(1)),:).';
    Ftmp(Ftmp == quadPoints.') = newVid;
    Fext(iqF(orderSub(1)),:) = Ftmp.';

    Ftmp = Fext(iqF(orderSub(2)),:).';
    Ftmp(Ftmp == quadPoints.') = newVid;
    Fext(iqF(orderSub(2)),:) = Ftmp.';

    %F(iqF(orderSub(1)),:) = [qV(orderSub(1)),newVid];
    %F(iqF(orderSub(2)),:) = [newVid,qV(orderSub(2))];
    sw = Fext(:,1) > Fext(:,2);
    Fext(sw,:) = fliplr(Fext(sw,:));

    [iqD,~] = find(iqD.'); iqD = reshape(iqD,2,[]).';

    % if we have different grains - we need a new boundary
    newBd = full(sum(I_DG(iqD(:,1),:) .* I_DG(iqD(:,2),:),2)) == 0;

    % add new edges
    Fext = [Fext; [quadPoints(newBd),newVid(newBd)]];
    qAdded = sum(newBd);

    % new rows to I_FDext
    I_FDext = [I_FDext; ...
      sparse(repmat((1:qAdded).',1,2), iqD(newBd,:), 1, ...
      qAdded,size(I_FDext,2))];

    % new empty rows to I_FDint
    %I_FDint = [I_FDint; sparse(qAdded,size(I_FDint,2))];

  end

  function mergeQuadrupleGrains

    gB = grains.boundary; gB = gB(length(gB)+1-(1:qAdded));
    toMerge = false(size(gB));

    for iPhase = ebsd.indexedPhasesId

      % restrict to the same phase
      iBd = all(gB.phaseId == iPhase,2);

      if ~any(iBd), continue; end

      % check for misorientation angle % TODO
      toMerge(iBd) = angle(gB(iBd).misorientation) < 5 * degree;

    end

    [grains, parentId] = merge(grains,gB(toMerge));

    % update I_DG
    I_PC = sparse(1:length(parentId),parentId,1);
    I_DG = I_DG * I_PC;

    % update grain ids
    [grainId,~] = find(I_DG.');
  end

end

function poly = calcPolygons(I_FG,F,V)
%
% Input:
%  I_FG - incidence matrix faces to grains
%  F    - list of faces
%  V    - list of vertices

poly = cell(size(I_FG,2),1);

if isempty(I_FG), return; end

% for all grains
for k=1:size(I_FG,2)

  % inner and outer boundaries are circles in the face graph
  EC = EulerCycles(F(I_FG(:,k)>0,:));

  % first cicle should be positive and all others negatively oriented
  for c = 1:numel(EC)
    if xor( c==1 , polySgnArea(V(EC{c},1),V(EC{c},2))>0 )
      EC{c} = fliplr(EC{c});
    end
  end

  % this is needed
  for c=2:numel(EC), EC{c} = [EC{c} EC{1}(1)]; end

  poly{k} = [EC{:}];

end

end


text:
function bti_test03 ( )

%*****************************************************************************80
%
%% BTI_TEST03 tests BURGERS_TIME_INVISCID.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BTI_TEST03\n' );
  fprintf ( 1, '  Test BURGERS_TIME_INVISCID\n' );

  method = 3;
  nx = 81;
  nt = 200;
  t_max = 2.0;
  bc = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Method: 3, Lax Friedrichs.\n' );
  fprintf ( 1, '  Initial condition: expansion\n' );
  fprintf ( 1, '  Number of space nodes = %d\n', nx );
  fprintf ( 1, '  Number of time steps = %d\n', nt );
  fprintf ( 1, '  Final time T_MAX = %g\n', t_max );
  fprintf ( 1, '  Boundary condition = %d\n', bc );

  U = burgers_time_inviscid ( method, @ic_expansion, nx, nt, t_max, bc );

  x = linspace ( -1.0, +1.0, nx );

  figure ( 3 )

  plot ( x, U(1:50:(nt+1),:), 'Linewidth', 3 )
  grid on
  xlabel ( '<-- X -->' )
  ylabel ( '<-- U(X,T) -->' )
  title ( 'Burgers, inviscid, initial expansion, Lax Friedrichs' )

  filename = 'bti_test03.png';
  print ( '-dpng', filename )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saved plot as "%s"\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Gtomosyn_test.m
% Copyright 2010-07-27, Jeff Fessler, University of Michigan

% small system for basic test
if ~isvar('As'), printm 'setup small'
	f.down = 16;
	f.dx = 0.1; % 100 um
	f.dz = 1; % 1 mm
	igs = image_geom('nx', 2304, 'ny', 48, 'nz', 1920, ...
		'dx', f.dx, 'dy', f.dz, 'dz', f.dx, ...
		'down', f.down);

	cgs = ct_geom('fan', 'ns', 2304, 'nt', 1920, 'na', 1, ...
			'orbit', 0, ...
			'ds', f.dx, 'dt', f.dx, ...
			'dfs', inf, ... % flat
			'dsd', 650, ...
			'dod', abs(igs.ny*igs.dy/2), ... % compressed on detector
			'down', f.down);

	As = Gtomosyn(cgs, igs, 'chat', 0); % set to 1 to see geometry!

	tmp = As * igs.unitv;
	im(tmp), cbar

%	Fatrix_test_basic(As, igs.mask)
	tester_tomo2(As, igs.mask, 'halt', 0, 'nblock', 0) % todo: ordering problem

	test_adjoint(As, 'big', 1, 'tol', 2e-6)
prompt
end

% big system for timing test
if ~isvar('A'), printm 'setup big'
	f.down = 4;
	f.dx = 0.1; % 100 um
	f.dz = 1; % 1 mm
	ig = image_geom('nx', 2304, 'ny', 48, 'nz', 1920, ...
		'dx', f.dx, 'dy', f.dz, 'dz', f.dx, ...
		'down', f.down);

	cg = ct_geom('fan', 'ns', 2304, 'nt', 1920, 'na', 1, ...
			'orbit', 0, ...
			'ds', f.dx, 'dt', f.dx, ...
			'dfs', inf, ... % flat
			'dsd', 650, ...
			'dod', abs(igs.ny*igs.dy/2), ... % compressed on detector
			'down', f.down);

	A = Gtomosyn(cg, ig, 'chat', 0);
end

if 1 && jf('ncore') >= 8, printm('time with %d threads', jf('ncore'))
	x = ig.circ;

	tmp1 = sprintf('%s, A*x time:', A.arg.blocks{1}.type);
	cpu etic
	y1 = A * x;
	cpu('etoc', tmp1) % 0.8 s on ir71 (20 cores)

	A2 = Gtomosyn(cg, ig, 'chat', 0, 'type', 'sf2');
	tmp2 = sprintf('%s, A*x time:', A2.arg.blocks{1}.type);
	cpu etic
	y2 = A2 * x;
	cpu('etoc', tmp2) % 1.4 s

	equivs(y1, y2)
end

% todo: compare to analytical - oblique views look funny

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [net, info] = cnn_dicnn_of(varargin)
%CNN_DICNN_OF Fine-tunes a pre-trained CNN with dynamic images on optical
% (DOF in pami journal) flow frames on UCF101 dataset

run(fullfile(fileparts(mfilename('fullpath')), ...
  '..', 'matconvnet', 'matlab', 'vl_setupnn.m')) ;

run(fullfile(fileparts(mfilename('fullpath')), ...
  '..', 'matconvnet', 'contrib', 'mcnExtraLayers', 'setup_mcnExtraLayers.m')) ;

run(fullfile(fileparts(mfilename('fullpath')), ...
  '..', 'matconvnet', 'contrib', 'autonn', 'setup_autonn.m')) ;

addpath Layers Datasets

opts.dataDir = fullfile('data','UCF101') ;
opts.expDir  = fullfile('exp', 'UCF101') ;
opts.modelPath = fullfile('models','resnext_50_32x4d-pt-mcn.mat.mat') ;
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.numFetchThreads = 8 ;

opts.lite = false ;
opts.imdbPath = fullfile(opts.dataDir, 'imdb-of.mat');
opts.pool1Layer = 'conv0'; % before conv1
opts.pool1Type = 'arpool'; % before conv1
opts.pool2Layer = 'fc6'; % before conv1
opts.DropOutRate = 0.85 ;
opts.datasetFn = @cnn_ucf101_of_setup_data ;
opts.networkFn = @cnn_init_resnext ;
opts.network = [] ;

opts.split = 1; % data split
opts.reverseDyn = 0; % reverse video frames e.g.[N:-1:1]
opts.numDynImgs = 10 ;
opts.epochFactor = 5 ;

opts.train = struct() ;
opts.train.gpus = [];
opts.train.batchSize = 128 ;
opts.train.numSubBatches = 32 ;
opts.train.solver = [] ;
opts.train.prefetch = true ;
opts.train.learningRate = 1e-2 ;
opts.train.numEpochs = 30 ;
% opts.train.savePreds = true ;
opts.train.randomSeed = 0 ;

opts = vl_argparse(opts, varargin) ;
if ~isfield(opts.train, 'gpus'), opts.train.gpus = []; end;

% -------------------------------------------------------------------------
%                                                              Prepare data
% -------------------------------------------------------------------------

if exist(opts.imdbPath,'file')
  imdb = load(opts.imdbPath) ;
else
  imdb = opts.datasetFn('dataDir', opts.dataDir, 'lite', opts.lite) ;
  mkdir(opts.expDir) ;
  save(opts.imdbPath, '-struct', 'imdb') ;
end

% UCF101 has 3 data splits
if opts.split>3
  error('split should be <=3');
end
imdb.images.set = imdb.images.sets(opts.split,:);

% reverse frame order
if opts.reverseDyn
  for i=1:numel(imdb.images.names)
    imdb.images.names{i} = imdb.images.names{i}(end:-1:1);
  end
end
% -------------------------------------------------------------------------
%                                                             Prepare model
% -------------------------------------------------------------------------
if isempty(opts.network)
  net = load(opts.modelPath);
  if isfield(net,'net')
    net = net.net;
  end
  opts.nCls = max(imdb.images.label) ;
  % net = dagnn.DagNN.loadobj(net) ;
  net = opts.networkFn(net,opts) ;

  % two channels instead of 3 RGB
  net.params(1).value = net.params(1).value(:,:,1:2,:) ;

  % Set the class names in the network
  net.meta.classes.name = imdb.classes.name ;
  net.meta.classes.description = imdb.classes.name ;
else
  assert(isa(opts.network,'dagnn.DagNN')) ;
  net = opts.network ;
end

% -------------------------------------------------------------------------
%                                                                     Learn
% -------------------------------------------------------------------------
if opts.epochFactor>0
  opts.train.train = repmat(find(imdb.images.set==1),[1 opts.epochFactor]) ;
else
  opts.train.train = NaN ;
  opts.train.numEpochs = 1 ;
end
opts.train.val = find(imdb.images.set==3) ;

[net, info] = cnn_train_dicnn_dag(net, imdb, getBatchFn(opts, net.meta), ...
                      'expDir', opts.expDir, ...
                      opts.train) ;

% -------------------------------------------------------------------------
%                                                          Report accuracy
% -------------------------------------------------------------------------
errlayer = net.getLayerIndex('errMC') ;

if ~isnan(errlayer)
  cats = imdb.classes.name ;
  accs = net.layers(errlayer).block.accuracy ;

  if numel(cats)~=numel(accs)
    error('wrong number of classes\n') ;
  end

  for i=1:numel(cats)
    fprintf('%s acc %.1f\n',cats{i},100*accs(i)) ;
  end
  fprintf('Mean accuracy %.1f\n',100*mean(accs)) ;
end
% -------------------------------------------------------------------------
function fn = getBatchFn(opts, meta)
% -------------------------------------------------------------------------
useGpu = numel(opts.train.gpus) > 0 ;

bopts.numThreads = opts.numFetchThreads ;
bopts.imageSize = meta.normalization.imageSize ;
if isfield(meta.normalization,'border')
  bopts.border = meta.normalization.border ;
else
  bopts.border = meta.normalization.imageSize(1:2) ./ ...
    meta.normalization.cropSize - meta.normalization.imageSize(1:2);
end

bopts.averageImage = 128 * ones([1 1 2],'single') ;
bopts.numDynImgs = opts.numDynImgs ;

fn = @(x,y) getDagNNBatch(bopts,useGpu,x,y) ;

% -------------------------------------------------------------------------
function inputs = getDagNNBatch(opts, useGpu, imdb, batch)
% -------------------------------------------------------------------------

% batch refers to videos (not for frames)
if isempty(batch)
  inputs = {'input', [], 'label', [], 'VideoId1', [], 'VideoId2', []};
  return;
end

isVal = ~isempty(batch) && imdb.images.set(batch(1)) ~= 1 ;

if ~isVal, transformation='multiScaleRegular'; else transformation='none';end

names = imdb.images.names(batch);

% images = strcat([imdb.imageDir filesep], imdb.images.name(batch)) ;

namesM = {};
nVids = numel(batch);

VideoId1 = [];
VideoId2 = [];

% step-size
stepSize = 6;

% pool nFrames into a dynamic image
nFrames = 10;
% number of dynamic images to be max pooled later
nDynImgs = opts.numDynImgs ;
opts = rmfield(opts,'numDynImgs') ;

c1 = 1;
for v=1:nVids

  name = names{v};
  nFrms = numel(name)/2;

  nSample = nFrames;

  if isVal
    startF = 1 ;
  else
    startF = ceil(stepSize/2) ;
  end
  nr = numel(startF:stepSize:nFrms);

  % jitter by removing 50 % and limit a batch to nMaxs * nSamples images
  if nr > 1 && (~isVal && nr>nDynImgs)
    rat = min(nDynImgs,ceil(0.50*nr));
    ri = randperm(nr);
    ri = ri(1:rat);
    r = zeros(1,nr);
    r(ri) = 1;
  else
    r = ones(1,nr);
  end

  c3 = 1;
  c2 = 0;

  for f=startF:stepSize:nFrms
    if r(c3)
      idx = f:min(f+nSample-1,nFrms) ;
      if numel(idx)<nFrames
        idx = [idx idx(end) * ones(1,nFrames-numel(idx))];
      end
      idxu = 2*idx - 1;
      idxv = 2*idx;
      idxuv = zeros(1,2 * numel(idxu)) ;
      idxuv(1:2:end) = idxu ;
      idxuv(2:2:end) = idxv ;

      namesM{end+1} = name(idxuv);
      VideoId1 = [VideoId1 c1 * ones(1,numel(idx))];
      c1 = c1 + 1;
      c2 = c2 + 1;
    end
    c3 = c3 + 1;
  end
  VideoId2 = [VideoId2 v * ones(1,c2) ] ;
end

images = strcat([imdb.imageDir filesep], horzcat(namesM{:}) ) ;

im = cnn_video_of_get_batch(images, VideoId1, opts, ...
  'transformation', transformation, 'prefetch', nargout == 0) ;

if nargout > 0
  if useGpu
    im = gpuArray(im) ;
  end
  inputs = {'input', im, 'label', imdb.images.label(batch), ...
    'VideoId1', VideoId1, 'VideoId2', VideoId2};

end

author: hbilen

repo: dynamic-image-nets


text:
% Create new IM
newIndex = createNewIM();

% Compute dose for a beam
gantryAngle = 70;
isoCenter.x = 0;
isoCenter.y = -55;
isoCenter.z = -115;
newIndex = 2;
planC = batchCalcDose(newIndex,gantryAngle,isoCenter);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function s = dirichlet_initial_s(m, bar_p)

K = length(m);
m = m/sum(m);
s = (K-1)/2/(-sum(m.*bar_p)+sum(m.*log(m)));

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function [ x, seed ] = uniform_01_order_sample ( n, seed )

%*****************************************************************************80
%
%% UNIFORM_01_ORDER_SAMPLE samples the Uniform 01 Order PDF.
%
%  Discussion:
%
%    In effect, this routine simply generates N samples of the
%    Uniform 01 PDF; but it generates them in order.  (Actually,
%    it generates them in descending order, but stores them in
%    the array in ascending order).  This saves the work of
%    sorting the results.  Moreover, if the order statistics
%    for another PDF are desired, and the inverse CDF is available,
%    then the desired values may be generated, presorted, by
%    calling this routine and using the results as input to the
%    inverse CDF routine.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Jerry Banks, editor,
%    Handbook of Simulation,
%    Engineering and Management Press Books, 1998, page 168.
%
%  Parameters:
%
%    Input, integer N, the number of elements in the sample.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X(N), N samples of the Uniform 01 PDF, in
%    ascending order.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  v = 1.0;
  for i = n : -1 : 1
    [ u, seed ] = r8_uniform_01 ( seed );
    v = v * u^( 1.0 / i );
    x(i) = v;
  end

  return
end


text:
function [obs,nobs]=matchobs(rtk,imud,obss)

global gls
imu_sample_rate=rtk.opt.ins.sample_rate;
time0=obss.data(:,1)+obss.data(:,2);

time1=imud.time.time+imud.time.sec;

idx=abs(time1-time0)<(0.501/imu_sample_rate);

if any(idx)
    obs0=obss.data(idx,:);  nobs=size(obs0,1);
else
    obs=NaN; nobs=0; return;
end

obs_tmp=repmat(gls.obs_tmp,nobs,1);
for i=1:nobs
    obs_tmp(i).time.time=obs0(i,1);
    obs_tmp(i).time.sec=obs0(i,2);
    obs_tmp(i).sat=obs0(i,3);
    obs_tmp(i).P=obs0(i,4:6);
    obs_tmp(i).L=obs0(i,7:9);
    obs_tmp(i).D=obs0(i,10:12);
    obs_tmp(i).S=obs0(i,13:15);
    obs_tmp(i).LLI=obs0(i,16:18);
    obs_tmp(i).code=obs0(i,19:21);
end

[obs,nobs]=exclude_sat(obs_tmp,rtk);

return

repo: GINav


text:
% Routine to display the time for a line for ZNLDetect
%
% Matt Allen, July 2005-Aug 2006

global ZNLD

% Find Extrapolation line and add to plot
lineh = gco;
if isempty(lineh) | isempty(get(get(lineh,'Parent'),'Parent'))
    error('No Line on Figure 102 was selected. Select a line and run the script again.');
end
if get(get(lineh,'Parent'),'Parent') ~= 102;
    error('No Line on Figure 102 was selected. Select a line and run the script again.');
end

disp(['Current line is number: ',num2str(get(gco,'UserData')),', for time response starting at: ',...
    num2str(ZNLD.ts_zc(get(gco,'UserData'))),' seconds.']);


text:
function box = scimat_box(scimat, m, a)
% SCIMAT_BOX  Compute tight box around SCIMAT segmentation.
%
% X = scimat_box(SCIMAT, M, A)
%
%   This function finds the vertices of a box tangent to the edges of a
%   segmentation. The box can be vertical or have any other orientation.
%
%   SCIMAT is a struct with the segmentation (see "help scimat" for
%   details).
%
%   X is a (3, 8)-matrix with the coordinates of the box vertices.
%
%   M is a 3-vector with the coordinates of the rotation centre. By
%   default, M=[0 0 0].
%
%   A is a (3, 3)-rotation matrix, in case the data has to be rotated
%   before computing the box. By default, A is the identity matrix and no
%   rotation is performed. To be consistent with the ITK (Insight Toolkit)
%   convention, the A matrix has to be the "backwards transformation", i.e.
%   A is the rotation from output to input voxel coordinates.
%
%   The "backwards transformation" is the transpose of the forward
%   transformation.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2010,2014 University of Oxford
% Version: 0.2.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(1, 3);
nargoutchk(0, 1);

% default
if (nargin < 2 || isempty(m))
    m = [0 0 0];
end
if (nargin < 3 || isempty(a))
    a = eye(3); % identity matrix
end

% remove the dummy dimension
scimat = scimat_squeeze(scimat);

% get hold of all the points in the segmentation mask

% extract linear indices of voxels in the segmentation
idx = find(scimat.data);

% get volume size
sz = size(scimat.data);

% convert linear index to multiple subscripts
[ir, ic, iz] = ind2sub(sz, idx);

% convert indices to real world coordinates and make column vectors
x = scimat_index2world([ir, ic, iz], scimat)';

% to make the interface consistent, we ask the user to input the backwards
% rotation (whenever we input the scimat volume, backwards rotation;
% whenever we input points, forwards rotation). But actually, we are going
% to operate on points, so we need to get the forwards rotation
a = a';

% avoid unnecessary operations if rotation is identity matrix
if all(all(a ~= eye(3)))

    % move points to centre of rotation
    for I = 1:3
        x(I, :) = x(I, :) - m(I);
    end

    % rotate real world coordinates
    x = a * x;

end

% find the limits of the data
minx = min(x, [], 2);
maxx = max(x, [], 2);

% create the box vertices
box = [...
    minx ...
    [minx(1) maxx(2) minx(3)]' ...
    [maxx(1) maxx(2) minx(3)]' ...
    [maxx(1) minx(2) minx(3)]' ...
    [minx(1) minx(2) maxx(3)]' ...
    [minx(1) maxx(2) maxx(3)]' ...
    maxx ...
    [maxx(1) minx(2) maxx(3)]' ];

if all(all(a ~= eye(3)))

    % undo rotation of vertices
    box = a' * box;

    % move vertices back from centre of rotation
    for I = 1:3
        box(I, :) = box(I, :) + m(I);
    end
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [x_world] = Robot2World(x_robot,x)
    %% Change of coordinates from robot fram to world frame

    x_world = zeros(size(x));
    for i = 1:size(x,2)
    x_world(1,i) = cos(x_robot(3))*(x(1,i)) - sin(x_robot(3))*(x(2,i)) + x_robot(1);
    x_world(2,i) = cos(x_robot(3))*(x(2,i)) + sin(x_robot(3))*(x(1,i)) + x_robot(2);
    end
end

author: Allopart

repo: rbpf-gmapping


text:
%% SUBSET_SUM_JOB_LOCAL distributes trials of subset_sum problem to multiple tasks.
%
%  Discussion:
%
%    Each task is given a portion of the range of possible combinations to check.
%
%    Each task returns a result that works or empty values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 May 2012
%
%  Author:
%
%    John Burkardt
%
  clear

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SUBSET_SUM_JOB_LOCAL\n' );
  fprintf ( 1, '  Set up and execute a job of independent tasks.\n' );
  fprintf ( 1, '  We seek a subset of weights that sum to a given target.\n' );
  fprintf ( 1, '  Each task examines a specific range of subsets.\n' );
%
%  Define the problem data.
%
  target = 2463098;

  weights = [ ...
    518533, 1037066, 2074132, 1648264, 796528, ...
   1593056,  686112, 1372224,  244448, 488896, ...
    977792, 1955584, 1411168,  322336, 644672, ...
   1289344,   78688,  157376,  314752, 629504, ...
   1259008 ];
%
%  Create the job.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Define the job:\n' );

  job = createJob ( ...
    'configuration', 'local', ...
    'FileDependencies', { 'subset_sum_task.m' } );
%
%  Compute the number of combinations to check.
%
  combos = fix ( 2^length(weights) );
%
%  Divide the range [0, COMBOS-1] among the tasks.
%
  task_num = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Define %d tasks:\n', task_num );

  endValue = -1;

  for task = 1 : task_num

    startValue = endValue + 1;
    endValue = ( task / task_num ) * combos - 1;

    task_id = createTask ( job, @subset_sum_task, 2, ...
                          { weights, target, [startValue, endValue] } );
  end
%
%  Submit the job.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Submit the job:\n' );

  submit ( job );
%
%  Wait for the job to finish.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Wait for the job:\n' );
  timestamp ( )
  wait ( job );
  timestamp ( )
%
%  The following commands will capture and print error messages
%  from the tasks, if something bad occurs.
%
  errmsgs = get ( job.Tasks, {'ErrorMessage'} );
  nonempty = ~cellfun ( @isempty, errmsgs );
  celldisp ( errmsgs(nonempty) );
%
%  Retrieve the results from all the tasks.
%
  results = getAllOutputArguments ( job );
%
%  Display the results (if any) from the tasks.
%
  fprintf ( 1, '\n' );
  for task = 1 : task_num
    if ~isempty ( results{task,1} )
      disp ( 'weights');       disp(results{task,1});
      disp ( 'weight_values'); disp(results{task,2});
    end
  end
%
%  Clean up by freeing memory associated with the job.
%
  destroy ( job );
%
%  Terminate;
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SUBSET_SUM_JOB_LOCAL:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );


text:
function [group_metrics individual_metrics values gwcsf gwcsfmean gwcsf_l2norm] = qc_metrics_second_level(obj, varargin)
%
% Quality metrics for a 2nd-level analysis (set of images from different subjects)
% The goal is to obtain measures that are simple and scale invariant, and
% so can be compared across image datasets
%
% :Usage:
% ::
%
%     [group_metrics individual_metrics values gwcsf gwcsfmean gwcsf_l2norm] = qc_metrics_second_level(list inputs here, [optional inputs])
%
% For objects: Type methods(object_name) for a list of special commands
%              Type help object_name.method_name for help on specific
%              methods.
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2016 Tor Wager
%
%     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/>.
% ..
%
% :Inputs:
%
%   **obj:**
%        An fmri_data object with a set of images
%        Intended to be individual-participant beta or contrast images from
%        first-level analyses.
%
% :Optional Inputs:
%   **noverbose:**
%        Turn off text display reporting of warnings and metrics
%
%   **param2:**
%        description of param2
%
% :Outputs:
%
%   **group_metrics:**
%        Quality metrics indicating potential problems, as described below
%
%   **individual_metrics:**
%        Individual image-level values for quality metrics
%        These can be used to identify individuals with problematic images
%        or as covariates in analyses.
%
%   **values:**
%        Grand Mean Gray, white, and CSF values - one number per metric per
%        image
%
%   **gwcsf:**
%        Gray, white, and CSF image objects with values for all images
%
%   **gwcsfmean:**
%        Gray, white, and CSF image objects mean values across images
%
%   **gwcsfl2norm:**
%        L2 norm of Gray, white, and CSF components for each image
%
% :References:
%   None.
%
%
% --------------------------------------------------------------------------------
% Coverage:  Percentage of standard MNI gray-matter mask with non-zero/non-nan values
%       NOT IMPLEMENTED YET
%       Range: [0 100],  100 is good (full coverage). Low values indicate
%       dropout due to artifacts or space mis-match
%
% Template match:  Dice coefficient (point-biserial r?) with standard MNI gray-matter mask
%       NOT IMPLEMENTED YET
%       Range: [0 100],  100 is good (full coverage). Low values indicate
%       dropout due to artifacts or space mis-match
%
% Bright-edge bias:  "brightedge"
%       NOT IMPLEMENTED YET
%       High inter-image variability around edge of brain
%       median absolute deviation (MAD) in brain edge-mask / MAD in center
%       of brain.
%       Range: [0 Inf],  1 is good (homogeneous). High values may indicate
%       head movement artifacts
%
% Signal in ventricles: "csf_to_gm_signal_ratio"
%       In most contrast maps, there should be less non-zero signal in the CSF space (e.g., ventricles)
%       than in gray matter.  This measures the ratio of mean absolute signal in CSF / GM.
%       Range: [0 Inf]. Values < 1 indicate low ventricle signal, which is good.
%       Higher values are worse.
%       Apply to: Individual or group. Contrast images (differences across conditions) for a series of individual participants
%
% Effect size in ventricles: "global_d_ventricles"
%       d = Mean / standard error (STD) of whole-ventricle average signal
%       Indicates a shift towards global activation or deactivation of the
%       ventricles in the group.
%       Range: [-Inf Inf], farther from zero is bad.
%       Apply to: Group. Contrast images (differences across conditions) for group of participants
%
% Significant global activation of ventricles: "global_logp_ventricles"
%       -log(p-value) in ventricle global mean across subjects, divided by -log(0.05)
%       Any value > 1 indicates a significant global activation in the
%       ventricles.
%       Range: [0 Inf], higher is worse. Values > 1 are bad.
%       Apply to: Group. Contrast images (differences across conditions) for group of participants
%
% Effect size in white matter: "global_d_wm"
%       d = Mean / standard error (STD) of whole-WM average signal
%       Indicates a shift towards global activation or deactivation of the
%       WM in the group. This is not completely redundant with CSF, because
%       some artifacts can affect one tissue compartment but not another.
%       WM values are, however, more likely to be influenced by true
%       signal.
%       Range: [-Inf Inf], farther from zero is bad.
%
% Significant global activation of WM: "global_logp_wm"
%       -log(p-value) in WM global mean across subjects, divided by -log(0.05)
%       Any value > 1 indicates a significant global activation in WM.
%       Range: [0 Inf], higher is worse. Values > 1 are bad.
%
% Non-spatially specific signal contamination: "r2_explained_by_csf"
%       Variance in individual differences in whole-gray-matter activity
%       Range: [-1 1], usually [0 1].  Farther from zero is bad.
%       Apply to: Group. Contrast images (differences across conditions) for group of participants
%
% Scale inhomogeneity across participants: "csf_scale_inhom"
%       Assumption: Subjects should be on the same scale (e.g., if you look
%       at histograms of contrast images for each subject). We assume that
%       the signal in gray matter may vary for meaningful reasons, but the
%       scale of the signal in CSF (and maybe WM?) should be similar across
%       participants.
%       Measure: The coefficient of variation of the L1-norm across CSF voxels for each participant.
%       The L1-norm is a measure of scale (across voxels, for each observation/image)
%       that is more robust than the variance.
%       High variability in these individual scale parameters is a sign of inhomogeneity.
%       The coefficient of variation is the standard dev. of this whole-image scale value divided by its mean.
%       This is a scale-invariant measure of inhomogeneity across images.
%       It also avoids the instability that occurs by dividing by values that may
%       be near zero, such as the mean.
%       Range: [0 Inf].  Lower is better.
%       Apply to: Group. Condition images (betas) for group of participants
%
% Scale inhomogeneity across participants in gray matter: "gm_scale_inhom"
%       The coefficient of variation of the L1-norm across GM voxels for each participant.
%       The L1-norm is a measure of scale (across voxels, for each observation/image)
%       that is more robust than the variance.
%       High variability in these individual scale parameters is a sign of inhomogeneity.
%       The coefficient of variation is the standard dev. of this whole-image scale value divided by its mean.
%       This is a scale-invariant measure of inhomogeneity across images.
%       It also avoids the instability that occurs by dividing by values that may
%       be near zero, such as the mean.
%       Range: [0 Inf].  Lower is better, but GM values may also be influenced by large differences in true signal.
%
% :Examples:
% ::
%    % Load sample images
%    obj = load_image_set('emotionreg');
%
%    % Run QC metrics
%    [group_metrics individual_metrics] = qc_metrics_second_level(obj);
%
%    % Plot histograms of gray, white, CSF values
%    hist_han = histogram(obj, 'byimage', 'by_tissue_type');
%
% :See also:
%   - extract_gray_white_csf, histogram
%

% ..
%    DEFAULTS AND INPUTS
% ..

doverbose = 1;

% optional inputs with default values
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case 'noverbose', doverbose = 0;

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

[values, components, gwcsf, gwcsf_l2norm] = extract_gray_white_csf(obj);

% Get mean signal in each tissue compartment for visualization, if desired
% -------------------------------------------------------------------------
gwcsfmean = mean(gwcsf{1});
gwcsfmean = replace_empty(gwcsfmean);
for i = 2:3
    m = mean(gwcsf{i});
    m = replace_empty(m);
    gwcsfmean.dat(:, i) = m.dat;
end

% -------------------------------------------------------------------------
% METRICS
% -------------------------------------------------------------------------

% Get coverage of gray-matter a priori mask space (do we have data in all vox?)
% -------------------------------------------------------------------------

objmasked = apply_mask(obj, mean(gwcsf{1}));

wh = isnan(objmasked.dat) | objmasked.dat == 0;
individual_metrics.gray_matter_coverage = sum(~wh) ./ size(objmasked.dat, 1);
group_metrics.mean_gray_matter_coverage = mean(individual_metrics.gray_matter_coverage);

%brightedge
% NOT DONE YET

% Global activation/deactivation in CSF and white matter
% -------------------------------------------------------------------------

individual_metrics.global_ventricle_values = mean(values(:, 3));

group_metrics.global_d_ventricles = mean(values(:, 3)) ./ std(values(:, 3));

n = size(values, 1);
t = group_metrics.global_d_ventricles .* sqrt(n);   % t-value
p = 2 .* min(tcdf(t, n - 1), tcdf(-t, n - 1));      % two-tailed

group_metrics.global_logp_ventricles = -log(p) ./ (-log(.05));

% --- wm ---
individual_metrics.global_white_matter_values = mean(values(:, 2));

group_metrics.global_d_wm = mean(values(:, 2)) ./ std(values(:, 2));

t = group_metrics.global_d_wm .* sqrt(n);   % t-value
p = 2 .* min(tcdf(t, n - 1), tcdf(-t, n - 1));      % two-tailed

group_metrics.global_logp_wm = -log(p) ./ (-log(.05));

% Gray-matter signal explained by CSF
% -------------------------------------------------------------------------

[r, p] = corr(values(:, 1), values(:, 3));

group_metrics.gm_explained_by_csf_pvalue = p;
group_metrics.r2_explained_by_csf = r .^ 2;

% Gray-matter scale (l2 norm) explained by CSF scale
% -------------------------------------------------------------------------

[r, p] = corr(gwcsf_l2norm(:, 1), gwcsf_l2norm(:, 3));

group_metrics.gm_l2norm_explained_by_csf_pvalue = p;
group_metrics.r2_l2norm_explained_by_csf = r .^ 2;

% Signal in ventricles
% -------------------------------------------------------------------------

csfsignal = mean(abs(gwcsf{3}.dat)) ./ mean(abs(gwcsf{1}.dat));
group_metrics.csf_to_gm_signal_ratio = mean(csfsignal);
individual_metrics.csf_to_gm_signal_ratio = csfsignal;

% Scale Inhomogeneity
% -------------------------------------------------------------------------

L1norm = sum(abs(gwcsf{1}.dat))';
group_metrics.gm_scale_inhom = std(L1norm) ./ mean(L1norm);
individual_metrics.gm_L1norm = L1norm;

L1norm = sum(abs(gwcsf{3}.dat))';
group_metrics.csf_scale_inhom = std(L1norm) ./ mean(L1norm);
individual_metrics.csf_L1norm = L1norm;

% Add warnings
% -------------------------------------------------------------------------
group_metrics.warnings = {};

if any(individual_metrics.gray_matter_coverage < .9)
    group_metrics.warnings{end + 1} = sprintf('Warning: Incomplete coverage, missing gray-matter data in some images.\n   - Mean coverage is d = %3.2f%%\n', group_metrics.mean_gray_matter_coverage .* 100);
end

if group_metrics.global_logp_ventricles > 1
    group_metrics.warnings{end + 1} = sprintf('Warning: Significant global activation in CSF space/ventricles.\n   - Effect size is d = %3.2f\n', group_metrics.global_d_ventricles);
end

if group_metrics.global_logp_wm > 1
    group_metrics.warnings{end + 1} = sprintf('Warning: Significant global activation in white matter.\n   - Effect size is d = %3.2f\n', group_metrics.global_d_wm);
end

if group_metrics.gm_explained_by_csf_pvalue < .05
    group_metrics.warnings{end + 1} = sprintf('Warning: Gray-matter individual diffs significantly correlated with mean CSF value.\n   - Var explained (r^2) = %3.2f%%', group_metrics.r2_explained_by_csf .* 100);
end

if group_metrics.gm_l2norm_explained_by_csf_pvalue < .05
    group_metrics.warnings{end + 1} = sprintf('Warning: Gray-matter scale (L2 norm) significantly correlated with mean CSF L2 norm.\n   - Var explained (r^2) = %3.2f%%', group_metrics.r2_l2norm_explained_by_csf .* 100);
end

if group_metrics.csf_to_gm_signal_ratio > 1
    group_metrics.warnings{end + 1} = sprintf('Warning: Strong non-zero signal in CSF relative to gray matter.\n   - Ratio is = %3.2f', group_metrics.csf_to_gm_signal_ratio);
end

if group_metrics.csf_scale_inhom > .3
    group_metrics.warnings{end + 1} = sprintf('Warning: High individual diffs in image scaling estimated from CSF.\n   - CV is = %3.2f', group_metrics.csf_scale_inhom);
end

% Display warnings
% -------------------------------------------------------------------------

if doverbose

    disp(' ');
    disp(group_metrics)
    disp(char(group_metrics.warnings{:}))

end

end % function

repo: CanlabCore


text:
filename='Cantilever_triangle_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'IPOPT';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.4;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function test_suite=test_meeg_senstype_collection
% tests for cosmo_meeg_senstype_collection
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_meeg_senstype_collection_
    if cosmo_skip_test_if_no_external('fieldtrip')
        return;
    end

    % get senstype properties
    % order is:
    % - senstype name
    % - acquisition system
    % - sensor type
    % - number of rows in channels
    % - number of columns in channels
    % - first and last channel label

    senstype_props=get_senstype_properties();

    % test each one
    n=numel(senstype_props);

    clear('cosmo_meeg_senstype_collection');
    stc=cosmo_meeg_senstype_collection();
    stc_cached=cosmo_meeg_senstype_collection();
    assertEqual(stc,stc_cached);
    for k=1:n
        prop=senstype_props{k};
        key=prop{1};
        is_required=numel(prop)>=7 && prop{7};

        if isfield(stc,key)
            sens=prop{2};
            type=prop{3};
            size_=[prop{4:5}];
            labels=prop{6}(:);

            assert(isfield(stc,key));
            st=stc.(key);

            assertEqual(st.sens,sens);
            assertEqual(st.type,type);
            assertEqual(size(st.label),size_);
            assertEqual(st.label([1;end]),labels);
        elseif is_required
            error('Missing key %s',key);
        end
    end

function test_meeg_senstype_collection_ft_regr_2016june
    if cosmo_skip_test_if_no_external('fieldtrip')
            return;
    end

    % since the following commit, ft_senstype supports
    %     ft_senslabel('neuromag306_combined')
    % but its output is different than was expected in
    % cosmo_meeg_senstype_collection. In particular the recent commit
    % includes the magnetometers, whereas the meeg_senstype
    % does not include those. This test is added as check for
    % such regressions

    % FieldTrip (https://github.com/fieldtrip/fieldtrip.git)
    % commit 882dba3426db583f7f4f9ac0cdf4eb3c26aaefc1
    % Author: Robert Oostenveld <r.oostenveld@gmail.com>
    % Date:   Wed Jun 22 12:47:47 2016 +0200
    %
    %     ENH - improve combined planar MEG sensor handling, see
    %     http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=3144

    stc=cosmo_meeg_senstype_collection();
    label=stc.neuromag306_planar_combined.label;
    assertEqual(size(label),[102 1]);

    present_labels={'MEG 2642+2643','MEG 0112+0113'};
    absent_labels={'MEG0111'',MEG0112','MEG0113',...
                    'MEG2641'',MEG2642','MEG2643'};
    assert(all(cosmo_match(present_labels,label(:))));
    assert(~any(cosmo_match(absent_labels,label(:))));

function props=get_senstype_properties()

    props={{ 'biosemi64',...
                    'biosemi64','eeg',...
                    64,1,{'A1','B32'}},...
            {'biosemi128',...
                    'biosemi128','eeg',...
                    128,1,{'A1','D32'}},...
            {'biosemi256',...
                    'biosemi256','eeg',...
                    256,1,{'A1','H32'}},...
            {'bti148',...
                    'bti148','meg_axial',...
                    148,1,{'A1','A148'}},...
            {'bti148_planar',...
                    'bti148','meg_planar',...
                    148,2,{'A1_dH','A148_dV'}},...
            {'bti148_planar_combined',...
                    'bti148','meg_planar_combined',...
                    148,1,{'A1','A148'},false},...  % not required
            {'bti248',...
                    'bti248','meg_axial',...
                    248,1,{'A1','A248'}},...
            {'bti248_planar',...
                    'bti248','meg_planar',...
                    248,2,{'A1_dH','A248_dV'}},...
            {'bti248_planar_combined',...
                    'bti248','meg_planar_combined',...
                    248,1,{'A1','A248'},false},...  % not required
            {'ctf151',...
                    'ctf151','meg_axial',...
                    151,1,{'MLC11','MZP02'}},...
            {'ctf151_planar',...
                    'ctf151','meg_planar',...
                    151,2,{'MLC11_dH','MZP02_dV'}},...
            {'ctf151_planar_combined',...
                    'ctf151','meg_planar_combined',...
                    151,1,{'MLC11','MZP02'}},...
            {'ctf275',...
                    'ctf275','meg_axial',...
                    275,1,{'MLC11','MZP01'}},...
            {'ctf275_planar',...
                    'ctf275','meg_planar',...
                    275,2,{'MLC11_dH','MZP01_dV'}},...
            {'ctf275_planar_combined',...
                    'ctf275','meg_planar_combined',...
                    275,1,{'MLC11','MZP01'}},...
            {'eeg1005',...
                    'ext1020','eeg',...
                    343,1,{'Fp1','T6'}},...
            {'eeg1010',...
                    'ext1020','eeg',...
                    94,1,{'Fp1','T6'}},...
            {'eeg1020',...
                    'ext1020','eeg',...
                    29,1,{'Fp1','T6'}},...
            {'egi32',...
                    'egi32','eeg',...
                    34,1,{'E1','Cz'}},...
            {'egi64',...
                    'egi64','eeg',...
                    66,1,{'E1','Cz'}},...
            {'egi128',...
                    'egi128','eeg',...
                    130,1,{'E1','Cz'}},...
            {'egi256',...
                    'egi256','eeg',...
                    258,1,{'E1','Cz'}},...
            {'neuromag122_planar_combined',...
                    'neuromag122','meg_planar_combined',...
                    61,1,{'MEG 001+002','MEG 121+122'}},...
            {'neuromag306_planar_combined',...
                    'neuromag306','meg_planar_combined',...
                    102,1,{'MEG 0112+0113','MEG 2642+2643'}},...
            {'itab153',...
                    'itab153','meg',...
                    153,1,{'MAG_000','MAG_152'}},...
            {'itab153_planar',...
                    'itab153','meg_planar',...
                    153,2,{'MAG_000_dH','MAG_152_dV'}},...
            {'itab153_planar_combined',...
                    'itab153','meg_planar_combined',...
                    153,1,{'MAG_000','MAG_152'}},...
            {'yokogawa64',...
                    'yokogawa64','meg_axial',...
                    64,1,{'AG001','AG064'}},...
            {'yokogawa64_planar',...
                    'yokogawa64','meg_planar',...
                    64,2,{'AG001_dH','AG064_dV'}},...
            {'yokogawa64_planar_combined',...
                    'yokogawa64','meg_planar_combined',...
                    64,1,{'AG001','AG064'}},...
            {'yokogawa160',...
                    'yokogawa160','meg_axial',...
                    160,1,{'AG001','AG160'}},...
            {'yokogawa160_planar',...
                    'yokogawa160','meg_planar',...
                    160,2,{'AG001_dH','AG160_dV'}},...
            {'yokogawa160_planar_combined',...
                    'yokogawa160','meg_planar_combined',...
                    160,1,{'AG001','AG160'}},...
            {'yokogawa440',...
                    'yokogawa440','meg_axial',...
                    412,1,{'AG001','RM412'}},...
            {'yokogawa440_planar',...
                    'yokogawa440','meg_planar',...
                    210,2,{'AG001_dH','AG392_dV'}},...
            {'yokogawa440_planar_combined',...
                    'yokogawa440','meg_planar_combined',...
                    210,1,{'AG001','AG392'}},...
            {'neuromag122alt_planar_combined',...
                    'neuromag122','meg_planar_combined',...
                    61,1,{'MEG001+002','MEG121+122'}},...
            {'neuromag306alt_planar_combined',...
                    'neuromag306','meg_planar_combined',...
                    102,1,{'MEG0112+0113','MEG2642+2643'}},...
            {'neuromag306_planar',...
                    'neuromag306','meg_planar',...
                    102,2,{'MEG 0112','MEG 2643'}},...
            {'neuromag306_mag',...
                    'neuromag306','meg_axial',...
                    102,1,{'MEG 0111','MEG 2641'}},...
            {'neuromag306alt_planar',...
                    'neuromag306','meg_planar',...
                    102,2,{'MEG0112','MEG2643'}},...
            {'neuromag306alt_mag',...
                    'neuromag306','meg_axial',...
                    102,1,{'MEG0111','MEG2641'}},...
            {'neuromag122_planar',...
                    'neuromag122','meg_planar',...
                    61,2,{'MEG 001','MEG 122'}},...
            {'neuromag122alt_planar',...
                    'neuromag122','meg_planar',...
                    61,2,{'MEG001','MEG122'}}};

    % the data anobe was produced by the code below
    %
    % sc=cosmo_meeg_senstype_collection();
    % names=fieldnames(sc);
    %
    % names={'biosemi64','biosemi128','biosemi256',...
    %         'bti148','bti148_planar','bti148_planar_combined',...
    %         'bti248','bti248_planar','bti248_planar_combined',...
    %         'ctf151','ctf151_planar','ctf151_planar_combined',...
    %         'ctf275','ctf275_planar','ctf275_planar_combined',...
    %         'eeg1005','eeg1010','eeg1020',...
    %         'egi32','egi64','egi128','egi256',...
    %         'neuromag122_planar_combined',...
    %         'neuromag306_planar_combined',...
    %         'itab153','itab153_planar','itab153_planar_combined',...
    %         'yokogawa64','yokogawa64_planar','yokogawa64_planar_combined',...
    %         'yokogawa160','yokogawa160_planar',...
    %                             'yokogawa160_planar_combined',...
    %         'yokogawa440','yokogawa440_planar',...
    %                             'yokogawa440_planar_combined',...
    %         'neuromag122alt_planar_combined',...
    %         'neuromag306alt_planar_combined',...
    %         'neuromag306_planar',...
    %         'neuromag306_mag',...
    %         'neuromag306alt_planar',...
    %         'neuromag306alt_mag',...
    %         'neuromag122_planar',...
    %         'neuromag122alt_planar'};
    % for k=1:numel(names)
    %     name=names{k};
    %     s=sc.(name);
    %     fprintf(['{''%s'',...\n\t\t''%s'',''%s'',...\n\t\t'...
    %                 '%d,%d,{''%s'',''%s''}},...\n'],...
    %             name,s.sens,s.type,size(s.label),s.label{1},s.label{end})
    % end
    %

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function coeff = eq_dist_coeff(dim,N,varargin)
%EQ_DIST_COEFF Coefficient of minimum distance of an EQ point set
%
%Syntax
% coeff = eq_dist_coeff(dim,N,options);
%
%Description
% COEFF = EQ_DIST_COEFF(dim,N) does the following:
% 1) uses the recursive zonal equal area sphere partitioning algorithm to
%    partition the unit sphere S^dim into N regions,
% 2) finds the EQ point set, the set of center points of each region,
% 3) finds the minimum Euclidean distance between points of the EQ point set,
% 4) sets COEFF to be the coefficient in the expression for the lower bound on
%    the minimum distance of a minimum energy point set:
%
%    DIST >= COEFF N^(-1/dim).
%
% The argument dim must be a positive integer.
% The argument N must be a positive integer or an array of positive integers.
% The result COEFF will be an array of the same size as N.
%
% COEFF = EQ_DIST_COEFF(dim,N,'offset','extra'), for dim == 2 or dim == 3, uses
% experimental extra rotation offsets to try to maximize the minimum distance.
% For dim > 3, extra offsets are not used.
%
%Notes
% The expression for the lower bound on minimum distance of a minimum r^(-s)
% energy point set on S^dim was given by [RakSZ95] for s == 0 and dim = 2,
% [Dahl78] for s == dim-1, [KuiSS04 Theorem 8] for dim-1 <= s < dim and
% [KuiS98 (1.12) p. 525] for s > dim.
%
% Ideally eq_dist_coeff(dim,N) should tend to area_of_sphere(dim)^(1/dim) as
% N goes to infinity.
%
%Examples
% > coeff=eq_dist_coeff(2,10)
%  coeff =
%      3.3250
%
% > coeff=eq_dist_coeff(3,1:6)
%  coeff =
%      2.0000    2.5198    2.0396    2.2449    2.4183    2.5698
%
%See also
% PARTITION_OPTIONS, EQ_MIN_DIST

% Copyright 2004-2005 Paul Leopardi for the University of New South Wales.
% $Revision 1.10 $ $Date 2005-06-01 $
% Documentation files renamed
% $Revision 1.00 $ $Date 2005-02-12 $
%
% For licensing, see COPYING.
% For references, see AUTHORS.
% For revision history, see CHANGELOG.

%
% Check number of arguments
%
error(nargchk(2,4,nargin));
%
% dim is the number of dimensions
% N is the number of regions
%
dist = eq_min_dist(dim,N,varargin{:});
coeff =  dist .* N.^(1/dim);
%
% end function

repo: gspbox