text:
function [data_train, labels_train, data_devel, labels_devel, raw_devel, PC, means_norm, stds_norm, vid_ids_devel_string] = ...
    Prepare_HOG_AU_data_generic_dynamic(train_users, devel_users, au_train, rest_aus, semaine_dir, feature_dir)

%%
addpath(genpath('../data extraction/'));

% First extracting the labels
[ labels_train, valid_ids_train, vid_ids_train ] = extract_SEMAINE_labels(semaine_dir, train_users, au_train);

[ labels_other, ~, ~ ] = extract_SEMAINE_labels(semaine_dir, train_users, rest_aus);
labels_other = cat(1, labels_other{:});

% Reading in the HOG data (of only relevant frames)
[train_appearance_data, valid_ids_train_hog, vid_ids_train_string] = Read_HOG_files_dynamic(train_users, vid_ids_train, feature_dir);

[train_geom_data] = Read_geom_files_dynamic(train_users, vid_ids_train, feature_dir);

% Subsample the data to make training quicker
labels_train = cat(1, labels_train{:});
valid_ids_train = logical(cat(1, valid_ids_train{:}));

reduced_inds = false(size(labels_train,1),1);
reduced_inds(labels_train == 1) = true;

% make sure the same number of positive and negative samples is taken
pos_count = sum(labels_train == 1);
neg_count = sum(labels_train == 0);

num_other = floor(pos_count / (size(labels_other, 2)));

inds_all = 1:size(labels_train,1);

for i=1:size(labels_other, 2)+1

    if(i > size(labels_other, 2))
        % fill the rest with a proportion of neutral
        inds_other = inds_all(sum(labels_other,2)==0 & ~labels_train );
        num_other_i = min(numel(inds_other), pos_count - sum(labels_train(reduced_inds,:)==0));
    else
        % take a proportion of each other AU
        inds_other = inds_all(labels_other(:, i) & ~labels_train );
        num_other_i = min(numel(inds_other), num_other);
    end
    inds_other_to_keep = inds_other(round(linspace(1, numel(inds_other), num_other_i)));
    reduced_inds(inds_other_to_keep) = true;

end

% Remove invalid ids based on CLM failing or AU not being labelled
reduced_inds(~valid_ids_train) = false;
reduced_inds(~valid_ids_train_hog) = false;

labels_other = labels_other(reduced_inds, :);
labels_train = labels_train(reduced_inds,:);
train_appearance_data = train_appearance_data(reduced_inds,:);
train_geom_data = train_geom_data(reduced_inds,:);
vid_ids_train_string = vid_ids_train_string(reduced_inds,:);

%% Extract devel data

% First extracting the labels
[ labels_devel, valid_ids_devel, vid_ids_devel ] = extract_SEMAINE_labels(semaine_dir, devel_users, au_train);

% Reading in the HOG data (of only relevant frames)
[devel_appearance_data, valid_ids_devel_hog, vid_ids_devel_string] = Read_HOG_files_dynamic(devel_users, vid_ids_devel, feature_dir);

[devel_geom_data] = Read_geom_files_dynamic(devel_users, vid_ids_devel, feature_dir);

labels_devel = cat(1, labels_devel{:});

% Peforming zone specific masking
if(au_train < 8 || au_train == 43 || au_train == 45) % upper face AUs ignore bottom face
    % normalise the data
    pca_file = '../../pca_generation/generic_face_upper.mat';
    load(pca_file);
elseif(au_train > 9) % lower face AUs ignore upper face and the sides
    % normalise the data
    pca_file = '../../pca_generation/generic_face_lower.mat';
    load(pca_file);
elseif(au_train == 9) % Central face model
    % normalise the data
    pca_file = '../../pca_generation/generic_face_rigid.mat';
    load(pca_file);
end

% Grab all data for validation as want good params for all the data
raw_devel = cat(2, devel_appearance_data, devel_geom_data);

devel_appearance_data = bsxfun(@times, bsxfun(@plus, devel_appearance_data, -means_norm), 1./stds_norm);
train_appearance_data = bsxfun(@times, bsxfun(@plus, train_appearance_data, -means_norm), 1./stds_norm);

data_train = train_appearance_data * PC;
data_devel = devel_appearance_data * PC;

data_train = cat(2, data_train, train_geom_data);
data_devel = cat(2, data_devel, devel_geom_data);

PC_n = zeros(size(PC)+size(train_geom_data, 2));
PC_n(1:size(PC,1), 1:size(PC,2)) = PC;
PC_n(size(PC,1)+1:end, size(PC,2)+1:end) = eye(size(train_geom_data, 2));
PC = PC_n;

means_norm = cat(2, means_norm, zeros(1, size(train_geom_data,2)));
stds_norm = cat(2, stds_norm, ones(1, size(train_geom_data,2)));

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function varargout = lvd_EditActionSetClampStateGUI(varargin)
% LVD_EDITACTIONSETCLAMPSTATEGUI MATLAB code for lvd_EditActionSetClampStateGUI.fig
%      LVD_EDITACTIONSETCLAMPSTATEGUI, by itself, creates a new LVD_EDITACTIONSETCLAMPSTATEGUI or raises the existing
%      singleton*.
%
%      H = LVD_EDITACTIONSETCLAMPSTATEGUI returns the handle to a new LVD_EDITACTIONSETCLAMPSTATEGUI or the handle to
%      the existing singleton*.
%
%      LVD_EDITACTIONSETCLAMPSTATEGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in LVD_EDITACTIONSETCLAMPSTATEGUI.M with the given input arguments.
%
%      LVD_EDITACTIONSETCLAMPSTATEGUI('Property','Value',...) creates a new LVD_EDITACTIONSETCLAMPSTATEGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before lvd_EditActionSetClampStateGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to lvd_EditActionSetClampStateGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help lvd_EditActionSetClampStateGUI

% Last Modified by GUIDE v2.5 03-Dec-2018 16:58:14

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @lvd_EditActionSetClampStateGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @lvd_EditActionSetClampStateGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before lvd_EditActionSetClampStateGUI is made visible.
function lvd_EditActionSetClampStateGUI_OpeningFcn(hObject, eventdata, handles, varargin)
    % This function has no output args, see OutputFcn.
    % hObject    handle to figure
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    % varargin   command line arguments to lvd_EditActionSetClampStateGUI (see VARARGIN)

    % Choose default command line output for lvd_EditActionSetClampStateGUI
    handles.output = hObject;

    centerUIFigure(hObject);

    action = varargin{1};
    setappdata(hObject,'action',action);

    populateGUI(handles, action);

    % Update handles structure
    guidata(hObject, handles);

    % UIWAIT makes lvd_EditActionSetClampStateGUI wait for user response (see UIRESUME)
    uiwait(handles.lvd_EditActionSetClampStateGUI);

function populateGUI(handles, action)

    if(action.activeStateToSet)
        set(handles.stateCombo,'Value',1);
    else
        set(handles.stateCombo,'Value',2);
    end

% --- Outputs from this function are returned to the command line.
function varargout = lvd_EditActionSetClampStateGUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
    if(isempty(handles))
        varargout{1} = false;
    else
        action = getappdata(hObject,'action');

        contents = cellstr(get(handles.stateCombo,'String'));
        stateStr = contents{get(handles.stateCombo,'Value')};

        switch stateStr
            case 'Active'
                state = true;
            case 'Inactive'
                state = false;
            otherwise
                error('Invalid state string: %s', stateStr);
        end

        action.activeStateToSet = state; %TODO Create this field

        varargout{1} = true;
        close(handles.lvd_EditActionSetClampStateGUI);
    end

% --- Executes on button press in saveAndCloseButton.
function saveAndCloseButton_Callback(hObject, eventdata, handles)
% hObject    handle to saveAndCloseButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    uiresume(handles.lvd_EditActionSetClampStateGUI);

% --- Executes on button press in cancelButton.
function cancelButton_Callback(hObject, eventdata, handles)
% hObject    handle to cancelButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    close(handles.lvd_EditActionSetClampStateGUI);

% --- Executes on selection change in stateCombo.
function stateCombo_Callback(hObject, eventdata, handles)
% hObject    handle to stateCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns stateCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from stateCombo

% --- Executes during object creation, after setting all properties.
function stateCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to stateCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on key press with focus on lvd_EditActionSetClampStateGUI or any of its controls.
function lvd_EditActionSetClampStateGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_EditActionSetClampStateGUI (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'enter'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'escape'
            close(handles.lvd_EditActionSetClampStateGUI);
    end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%diag_mat(A,B,2)=[B 0 0;0 A 0;0 0 A];
function [res]=diagmat(A,B,rep)

if (isempty(A))
    res=B;
    return;
end

if (isempty(rep))
    rep=1;
end

nrow=size(A,1);
ncol=size(A,2);

if (isempty(B))
    res=zeros(nrow*rep,ncol*rep);
    for in=1:rep
        res((in-1)*nrow+1:in*nrow,(in-1)*ncol+1:in*ncol)=A;
    end
else
    nrow1=size(B,1);
    ncol1=size(B,2);
    res=zeros(nrow*rep+nrow1,ncol*rep+ncol1);
    res(1:nrow1,1:ncol1)=B;
    for in=1:rep
        res((in-1)*nrow+1+nrow1:in*nrow+nrow1,(in-1)*ncol+1+ncol1:in*ncol+ncol1)=A;
    end
end

repo: ekfmonoslam


text:
% FDLA and FMMC solutions for a 50-node, 200-edge graph
% S. Boyd, et. al., "Convex Optimization of Graph Laplacian Eigenvalues"
% ICM'06 talk examples (www.stanford.edu/~boyd/cvx_opt_graph_lapl_eigs.html)
% Written for CVX by Almir Mutapcic 08/29/06
% (figures are generated)
%
% In this example we consider a graph described by the incidence matrix A.
% Each edge has a weight W_i, and we optimize various functions of the
% edge weights as described in the referenced paper; in particular,
%
% - the fastest distributed linear averaging (FDLA) problem (fdla.m)
% - the fastest mixing Markov chain (FMMC) problem (fmmc.m)
%
% Then we compare these solutions to the heuristics listed below:
%
% - maximum-degree heuristic (max_deg.m)
% - constant weights that yield fastest averaging (best_const.m)
% - Metropolis-Hastings heuristic (mh.m)

% randomly generate a graph with 50 nodes and 200 edges
% and make it pretty for plotting
n = 50; threshold = 0.2529;
rand('state',209);
xy = rand(n,2);

angle = 10*pi/180;
Rotate = [ cos(angle) sin(angle); -sin(angle) cos(angle) ];
xy = (Rotate*xy')';

Dist = zeros(n,n);
for i=1:(n-1);
  for j=i+1:n;
    Dist(i,j) = norm( xy(i,:) - xy(j,:) );
  end;
end;
Dist = Dist + Dist';
Ad = Dist < threshold;
Ad = Ad - eye(n);
m = sum(sum(Ad))/2;

% find the incidence matrix
A = zeros(n,m);
l = 0;
for i=1:(n-1);
  for j=i+1:n;
    if Ad(i,j)>0.5
      l = l + 1;
      A(i,l) =  1;
      A(j,l) = -1;
    end;
  end;
end;
A = sparse(A);

% Compute edge weights: some optimal, some based on heuristics
[n,m] = size(A);

[ w_fdla, rho_fdla ] = fdla(A);
[ w_fmmc, rho_fmmc ] = fmmc(A);
[ w_md,   rho_md   ] = max_deg(A);
[ w_bc,   rho_bc   ] = best_const(A);
[ w_mh,   rho_mh   ] = mh(A);

tau_fdla = 1/log(1/rho_fdla);
tau_fmmc = 1/log(1/rho_fmmc);
tau_md   = 1/log(1/rho_md);
tau_bc   = 1/log(1/rho_bc);
tau_mh   = 1/log(1/rho_mh);

eig_opt  = sort(eig(eye(n) - A * diag(w_fdla) * A'));
eig_fmmc = sort(eig(eye(n) - A * diag(w_fmmc) * A'));
eig_mh   = sort(eig(eye(n) - A * diag(w_mh)   * A'));
eig_md   = sort(eig(eye(n) - A * diag(w_md)   * A'));
eig_bc   = sort(eig(eye(n) - A * diag(w_bc)   * A'));

fprintf(1,'\nResults:\n');
fprintf(1,'FDLA weights:\t\t rho = %5.4f \t tau = %5.4f\n',rho_fdla,tau_fdla);
fprintf(1,'FMMC weights:\t\t rho = %5.4f \t tau = %5.4f\n',rho_fmmc,tau_fmmc);
fprintf(1,'M-H weights:\t\t rho = %5.4f \t tau = %5.4f\n',rho_mh,tau_mh);
fprintf(1,'MAX_DEG weights:\t rho = %5.4f \t tau = %5.4f\n',rho_md,tau_md);
fprintf(1,'BEST_CONST weights:\t rho = %5.4f \t tau = %5.4f\n',rho_bc,tau_bc);

% plot results
figure(1), clf
gplot(Ad,xy);
hold on;
plot(xy(:,1), xy(:,2), 'ko','LineWidth',4, 'MarkerSize',4);
axis([0.05 1.1 -0.1 0.95]);
title('Graph')
hold off;

figure(2), clf
v_fdla = [w_fdla; diag(eye(n) - A*diag(w_fdla)*A')];
[ifdla, jfdla, neg_fdla] = find( v_fdla.*(v_fdla < -0.001 ) );
v_fdla(ifdla) = [];
wbins = [-0.6:0.012:0.6];
hist(neg_fdla,wbins); hold on,
h = findobj(gca,'Type','patch');
set(h,'FaceColor','r')
hist(v_fdla,wbins); hold off,
axis([-0.6 0.6 0 12]);
xlabel('optimal FDLA weights');
ylabel('histogram');

figure(3), clf
xbins = (-1:0.015:1)';
ymax  = 6;
subplot(3,1,1)
hist(eig_md, xbins); hold on;
max_md = max(abs(eig_md(1:n-1)));
plot([-max_md -max_md],[0 ymax], 'b--');
plot([ max_md  max_md],[0 ymax], 'b--');
axis([-1 1 0 ymax]);
text(0,5,'MAX DEG');
title('Eigenvalue distributions')
subplot(3,1,2)
hist(eig_bc, xbins); hold on;
max_opt = max(abs(eig_bc(1:n-1)));
plot([-max_opt -max_opt],[0 ymax], 'b--');
plot([ max_opt  max_opt],[0 ymax], 'b--');
axis([-1 1 0 ymax]);
text(0,5,'BEST CONST');
subplot(3,1,3)
hist(eig_opt, xbins); hold on;
max_opt = max(abs(eig_opt(1:n-1)));
plot([-max_opt -max_opt],[0 ymax], 'b--');
plot([ max_opt  max_opt],[0 ymax], 'b--');
axis([-1 1 0 ymax]);
text(0,5,'FDLA');

figure(4), clf
xbins = (-1:0.015:1)';
ymax  = 6;
subplot(3,1,1)
hist(eig_md, xbins); hold on;
max_md = max(abs(eig_md(1:n-1)));
plot([-max_md -max_md],[0 ymax], 'b--');
plot([ max_md  max_md],[0 ymax], 'b--');
axis([-1 1 0 ymax]);
text(0,5,'MAX DEG');
title('Eigenvalue distributions')
subplot(3,1,2)
hist(eig_mh, xbins); hold on;
max_opt = max(abs(eig_mh(1:n-1)));
plot([-max_opt -max_opt],[0 ymax], 'b--');
plot([ max_opt  max_opt],[0 ymax], 'b--');
axis([-1 1 0 ymax]);
text(0,5,'MH');
subplot(3,1,3)
hist(eig_fmmc, xbins); hold on;
max_opt = max(abs(eig_fmmc(1:n-1)));
plot([-max_opt -max_opt],[0 ymax], 'b--');
plot([ max_opt  max_opt],[0 ymax], 'b--');
axis([-1 1 0 ymax]);
text(0,5,'FMMC');

figure(5), clf
v_fmmc = [w_fmmc; diag(eye(n) - A*diag(w_fmmc)*A')];
[ifmmc, jfmmc, nonzero_fmmc] = find( v_fmmc.*(v_fmmc > 0.001 ) );
hist(nonzero_fmmc,80);
axis([0 1 0 10]);
xlabel('optimal positive FMMC weights');
ylabel('histogram');

figure(6), clf
An = abs(A*diag(w_fmmc)*A');
An = (An - diag(diag(An))) > 0.0001;
gplot(An,xy,'b-'); hold on;
h = findobj(gca,'Type','line');
set(h,'LineWidth',2.5)
gplot(Ad,xy,'b:');
plot(xy(:,1), xy(:,2), 'ko','LineWidth',4, 'MarkerSize',4);
axis([0.05 1.1 -0.1 0.95]);
title('Subgraph with positive transition prob.')
hold off;

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [ v_num, e_num, t_num ] = ffmsh_2d_size_example ( )

%*****************************************************************************80
%
%% FFMSH_2D_SIZE_EXAMPLE returns sizes for the 2D example.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 December 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer V_NUM, the number of vertices.
%
%    Output, integer E_NUM, the number of boundary edges.
%
%    Output, integer T_NUM, the number of triangles.
%
  e_num = 10;
  t_num = 18;
  v_num = 15;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [iNewSurfaces, OutputSurfacesFiles, nVertices] = import_surfaces(iSubject, SurfaceFiles, FileFormat, isApplyMriOrient, OffsetMri, SelLabels, Comment)
% IMPORT_SURFACES: Import a set of surfaces in a Subject of Brainstorm database.
%
% USAGE: iNewSurfaces = import_surfaces(iSubject, SurfaceFiles, FileFormat, offset=[], SelLabels=[all], Comment=[])
%        iNewSurfaces = import_surfaces(iSubject)   : Ask user the files to import
%
% INPUT:
%    - iSubject     : Indice of the subject where to import the surfaces
%                     If iSubject=0 : import surfaces in default subject
%    - SurfaceFiles : Cell array of full filenames of the surfaces to import (format is autodetected)
%                     => if not specified : files to import are asked to the user
%    - FileFormat   : String representing the file format to import.
%                     Please see in_tess.m to get the list of supported file formats
%    - isApplyMriOrient: {0,1}
%    - OffsetMri    : (x,y,z) values to add to the coordinates of the surface before converting it to SCS
%    - SelLabels    : Cell-array of labels, when importing atlases
%    - Comment      : Comment of the output file
%
% OUTPUT:
%    - iNewSurfaces : Indices of the surfaces added in database

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

%% ===== PARSE INPUTS =====
% Check command line
if ~isnumeric(iSubject) || (iSubject < 0)
    error('Invalid subject indice.');
end
if (nargin < 3) || isempty(SurfaceFiles)
    SurfaceFiles = {};
    FileFormat = [];
else
    if ischar(SurfaceFiles)
        SurfaceFiles = {SurfaceFiles};
    end
    if (nargin == 2) || ((nargin >= 3) && isempty(FileFormat))
        error('When you pass a SurfaceFiles argument, FileFormat must be defined too.');
    end
end
if (nargin < 4) || isempty(isApplyMriOrient)
    isApplyMriOrient = [];
end
if (nargin < 5) || isempty(OffsetMri)
    OffsetMri = [];
end
if (nargin < 6) || isempty(SelLabels)
    SelLabels = [];
end
if (nargin < 7) || isempty(Comment)
    Comment = [];
end
iNewSurfaces = [];
OutputSurfacesFiles = {};
nVertices = [];
Labels = [];

% Get Protocol information
ProtocolInfo = bst_get('ProtocolInfo');
% Get subject directory
sSubject = bst_get('Subject', iSubject);
subjectSubDir = bst_fileparts(sSubject.FileName);
% Check the presence of the MRI: warning if no MRI
if isempty(sSubject.Anatomy)
    res = java_dialog('confirm', ...
        ['WARNING: To import correctly surface files, the subject''s MRI is needed.' 10 10 ...
        'Import subject''s MRI now?' 10 10], 'Import surfaces');
    if res
        import_mri(iSubject, [], [], 1);
        return
    end
end

%% ===== SELECT SURFACE FILES =====
% If surface files to load are not defined : open a dialog box to select it
if isempty(SurfaceFiles)
    % Get last used directories and formats
    LastUsedDirs = bst_get('LastUsedDirs');
    DefaultFormats = bst_get('DefaultFormats');
    if isempty(DefaultFormats.SurfaceIn)
        DefaultFormats.SurfaceIn = 'ALL';
    end
    % Get Surface files
    [SurfaceFiles, FileFormat, FileFilter] = java_getfile( 'open', ...
       'Import surfaces...', ...     % Window title
       LastUsedDirs.ImportAnat, ...   % Default directory
       'multiple', 'files', ...      % Selection mode
       bst_get('FileFilters', 'surface'), ...
       DefaultFormats.SurfaceIn);
    % If no file was selected: exit
    if isempty(SurfaceFiles)
        return
    end
    % Save default import directory
    LastUsedDirs.ImportAnat = bst_fileparts(SurfaceFiles{1});
    bst_set('LastUsedDirs', LastUsedDirs);
    % Save default import format
    DefaultFormats.SurfaceIn = FileFormat;
    bst_set('DefaultFormats',  DefaultFormats);
end

%% ===== APPLY MRI TRANSFORM =====
% Load MRI
if ~isempty(sSubject.Anatomy)
    sMri = bst_memory('LoadMri', sSubject.Anatomy(sSubject.iAnatomy).FileName);
else
    sMri = [];
end
% If user transformation on MRI: ask to apply transformations on surfaces
isMni = isequal(FileFormat, 'MRI-MASK-MNI');
if ~isMni && isempty(isApplyMriOrient) && ~isempty(sMri) && isfield(sMri, 'InitTransf') && ~isempty(sMri.InitTransf)
    isApplyMriOrient = java_dialog('confirm', ['MRI orientation was non-standard and had to be reoriented.' 10 10 ...
                                   'Apply the same transformation to the surfaces ?' 10 ...
                                   'Default answer is: NO', 10 10], 'Import surfaces');
else
    isApplyMriOrient = 0;
end

%% ===== LOAD EACH SURFACE =====
% Process all the selected surfaces
for iFile = 1:length(SurfaceFiles)
    TessFile = SurfaceFiles{iFile};

    % ===== LOAD SURFACE FILE =====
    bst_progress('start', 'Importing tesselation', ['Loading file "' TessFile '"...']);
    % Load surfaces(s)
    [Tess, Labels] = in_tess(TessFile, FileFormat, sMri, OffsetMri, SelLabels);
    if isempty(Tess)
        bst_progress('stop');
        return
    end

    % ===== INITIALIZE NEW SURFACE =====
    % Get imported base name
    if strcmpi(FileFormat, 'FS')
        [tmp__, fBase, fExt] = bst_fileparts(TessFile);
        importedBaseName = [fBase, strrep(fExt, '.', '_')];
    else
        [tmp__, importedBaseName] = bst_fileparts(TessFile);
    end
    importedBaseName = strrep(importedBaseName, 'tess_', '');
    importedBaseName = strrep(importedBaseName, '_tess', '');
    importedBaseName = strrep(importedBaseName, 'subjectimage_', '');
    importedBaseName = strrep(importedBaseName, '_volatlas', '');
    % Only one surface
    if (length(Tess) == 1)
        % Surface mesh
        if isfield(Tess, 'Faces')
            NewTess = db_template('surfacemat');
            NewTess.Comment  = Tess(1).Comment;
            NewTess.Vertices = Tess(1).Vertices;
            if isfield(Tess, 'Faces')   % Volume meshes do not have Faces field
                NewTess.Faces = Tess(1).Faces;
            end
        % Volume FEM mesh
        else
            NewTess = Tess;
        end
    % Multiple surfaces
    else
        [Tess(:).Atlas] = deal(db_template('Atlas'));
        NewTess = tess_concatenate(Tess, [], [], Labels);
        NewTess.iAtlas  = find(strcmpi({NewTess.Atlas.Name}, 'Structures'));
        NewTess.Comment = importedBaseName;
    end
    % Comment
    if ~isempty(Comment)
        NewTess.Comment = Comment;
    elseif isempty(NewTess.Comment)
        NewTess.Comment = importedBaseName;
    end

    % ===== APPLY MRI ORIENTATION =====
    if isApplyMriOrient
        % History: Apply MRI transformation
        NewTess = bst_history('add', NewTess, 'import', 'Apply transformation that was applied to the MRI volume');
        % Apply MRI transformation
        NewTess = ApplyMriTransfToSurf(sMri.InitTransf, NewTess);
    end

    % ===== SAVE BST FILE =====
    % History: File name
    NewTess = bst_history('add', NewTess, 'import', ['Import from: ' TessFile]);
    % Produce a default surface filename (surface of volume mesh)
    if isfield(NewTess, 'Faces')
        BstTessFile = bst_fullfile(ProtocolInfo.SUBJECTS, subjectSubDir, ['tess_' importedBaseName '.mat']);
    else
        BstTessFile = bst_fullfile(ProtocolInfo.SUBJECTS, subjectSubDir, ['tess_fem_' importedBaseName '.mat']);
    end
    % Make this filename unique
    BstTessFile = file_unique(BstTessFile);
    % Save new surface in Brainstorm format
    bst_save(BstTessFile, NewTess, 'v7');

    % ===== UPDATE DATABASE ======
    % Add new surface to database
    BstTessFileShort = file_short(BstTessFile);
    iNewSurfaces(end+1) = db_add_surface(iSubject, BstTessFileShort, NewTess.Comment);
    % Unload surface from memory (if this surface with the same name was previously loaded)
    bst_memory('UnloadSurface', BstTessFile);
    % Save output filename
    OutputSurfacesFiles{end+1} = BstTessFile;
    % Return number of vertices
    nVertices(end+1) = length(NewTess.Vertices);
end

% Save database
db_save();
bst_progress('stop');
end

%% ===== APPLY MRI ORIENTATION =====
function sSurf = ApplyMriTransfToSurf(MriTransf, sSurf)
    % Apply transformation to vertices
    sSurf.Vertices = ApplyMriTransfToPts(MriTransf, sSurf.Vertices);
    % Update faces order: If the surfaces were flipped an odd number of times, invert faces orientation
    if (mod(nnz(strcmpi(MriTransf(:,1), 'flipdim')), 2) == 1)
        sSurf.Faces = sSurf.Faces(:,[1 3 2]);
    end
end

function pts = ApplyMriTransfToPts(MriTransf, pts)
    % Apply step by step all the transformations that have been applied to the MRI
    for i = 1:size(MriTransf,1)
        ttype = MriTransf{i,1};
        val   = MriTransf{i,2};
        switch (ttype)
            case 'flipdim'
                % Detect the dimensions that have constantly negative coordinates
                iDimNeg = find(sum(sign(pts) == -1) == size(pts,1));
                if ~isempty(iDimNeg)
                    pts(:,iDimNeg) = -pts(:,iDimNeg);
                end
                % Flip dimension
                pts(:,val(1)) = val(2)/1000 - pts(:,val(1));
                % Restore initial negative values
                if ~isempty(iDimNeg)
                    pts(:,iDimNeg) = -pts(:,iDimNeg);
                end
            case 'permute'
                pts = pts(:,val);
            case 'vox2ras'
                % Do nothing, applied earlier
        end
    end
end

author: brainstorm-tools


text:
function dt6 = dtiRawFixDt6File(dt6FileName,t1FileName)
% Fix the paths in a dt6 file.
% dt6 = dtiRawFixDt6File(dt6FileName,[t1FileName='t1/t1.nii.gz'])
%
% Note that this function enforces unix-style file seperators (/) because
% these work on all platforms (Windows prefers \, but will take /).
%
% Also note that if you catch the output argument, the fixed dt6 struct
% will NOT be saved.
%
% TODO: Allow t1 file to be outside the subject dir.
%
% 2007.08.01 RFD wrote it based on a script from Sherbondy.

if(~exist('dt6FileName','var')||isempty(dt6FileName))
    [f,p] = uigetfile({'*.mat';'*.*'},'Select a dt6 file to fix...');
    if(isnumeric(f)), disp('User canceled.'); return; end
    dt6FileName = fullfile(p,f);
end
if(~exist('t1FileName','var')||isempty(t1FileName))
    t1FileName='t1/t1.nii.gz';
end

[dt6Dir,fName,ext] = fileparts(dt6FileName);
if(isempty(dt6Dir)) dt6Dir = pwd; end
[subDir,dt6DirName] = fileparts(dt6Dir);
if(isempty(subDir)) subDir = pwd; end

dt6 = load(dt6FileName);
if isfield(dt6.files,'homeDir'); dt6.files = rmfield(dt6.files,'homeDir'); end
if isfield(dt6.files,'binDir'); dt6.files = rmfield(dt6.files,'binDir'); end
fieldList = fieldnames(dt6.files);
for ff=1:length(fieldList)
    [p,f,e] = fileparts(getfield(dt6.files,fieldList{ff}));
    f = [f e];
	dt6.files = setfield(dt6.files,fieldList{ff},[dt6DirName '/bin/',f]);
end
dt6.files.t1 = t1FileName;

if(~exist(fullfile(subDir,dt6.files.t1),'file'))
    [f,p] = uigetfile({'*.mat';'*.*'},'Select the associated t1 file (cancel for none)...');
    if(isnumeric(f))
        dt6.files.t1 = '';
    else
        dt6.files.t1 = [p '/' f];
    end
end
if(nargout==0)
    save(dt6FileName, '-struct', 'dt6');
end

return;

baseDir = '/teal/scr1/dti/probtrack_compare';
subjDirs = {'bg040719' 'md040714'  'mho040625'  'rfd040630'  'ss040804' 'db061209'};
ogDir = pwd;
for dd=1:length(subjDirs)
    disp(['Fixing dt6 for subject: ' subjDirs{dd} ' ...']);
    dt6FileName = fullfile(baseDir,subjDirs{dd},'dti06','dt6.mat');
    dtiRawFixDt6File(dt6FileName);
end
cd(ogDir);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Copyright 2009 - 2010 The MathWorks, Inc.
function y = kalman01(z) %#eml
% Initialize state transition matrix
dt=1;
A=[ 1 0 dt 0 0 0;...
    0 1 0 dt 0 0;...
    0 0 1 0 dt 0;...
    0 0 0 1 0 dt;...
    0 0 0 0 1 0 ;...
    0 0 0 0 0 1 ];

% Measurement matrix
H = [ 1 0 0 0 0 0; 0 1 0 0 0 0 ];
Q = eye(6);
R = 1000 * eye(2);

% Initial conditions
persistent x_est p_est
if isempty(x_est)
    x_est = zeros(6, 1);
    p_est = zeros(6, 6);
end

% Predicted state and covariance
x_prd = A * x_est;
p_prd = A * p_est * A' + Q;

% Estimation
S = H * p_prd' * H' + R;
B = H * p_prd';
klm_gain = (S \ B)';

% Estimated state and covariance
x_est = x_prd + klm_gain * (z - H * x_prd);
p_est = p_prd - klm_gain * H * p_prd;

% Compute the estimated measurements
y = H * x_est;
end

author: Sable

repo: mcbench-benchmarks


text:
function [H, H_cross]= stimutil_cueArrows(dirs, varargin)

props= {'cross'         0       '!DOUBLE[1]'
       };
props_drawArrow= stimutil_drawArrow;
all_props= opt_catProps(props, props_drawArrow);

if nargin==0,
  H= all_props;
  return
end

opt= opt_proplistToStruct(varargin{:});
opt= opt_setDefaults(opt, props);
opt_checkProplist(opt, all_props);

for dd= 1:length(dirs),
  if dd<length(dirs),
    h_tmp= stimutil_drawArrow(dirs{dd}, opt, 'cross',0);
  else
    h_tmp= stimutil_drawArrow(dirs{dd}, opt);
  end
  H(dd)= h_tmp.arrow;
end

if opt.cross,
  H_cross= h_tmp.cross;
else
  H_cross= [];
end

set([H H_cross], 'Visible','off');

author: bbci


text:
function element_node = grid_t10_element ( nelemx, nelemy )

%*****************************************************************************80
%
%% GRID_T10_ELEMENT produces a grid of pairs of 10 node triangles.
%
%  Example:
%
%    Input:
%
%      NELEMX = 2, NELEMY = 2
%
%    Output:
%
%      ELEMENT_NODE =
%         1,  2,  3,  4, 10, 16, 22, 15,  8,  9;
%        25, 24, 23, 22, 16, 10,  4, 11, 18, 17;
%         4,  5,  6,  7, 13, 19, 25, 18, 11, 12;
%        28, 27, 26, 25, 19, 13,  7, 14, 21, 20;
%        22, 23, 24, 25, 31, 37, 43, 36, 29, 30;
%        46, 45, 44, 43, 37, 31, 25, 32, 39, 38;
%        25, 26, 27, 28, 34, 40, 46, 39, 31, 33;
%        49, 48, 47, 46, 40, 34, 28, 35, 42, 41.
%
%  Grid:
%
%   43-44-45-46-47-48-49
%    |\     6 |\     8 |
%    | \      | \      |
%   36 37 38 39 40 41 42
%    |   \    |   \    |
%    |    \   |    \   |
%   29 30 31 32 33 34 35
%    |      \ |      \ |
%    | 5     \| 7     \|
%   22-23-24-25-26-27-28
%    |\     2 |\     4 |
%    | \      | \      |
%   15 16 17 18 19 20 21
%    |   \    |   \    |
%    |    \   |    \   |
%    8  9 10 11 12 13 14
%    |      \ |      \ |
%    | 1     \| 3     \|
%    1--2--3--4--5--6--7
%
%  Reference Element T10:
%
%    |
%    1  10
%    |  |\
%    |  | \
%    |  8  9
%    |  |   \
%    S  |    \
%    |  5  6  7
%    |  |      \
%    |  |       \
%    0  1--2--3--4
%    |
%    +--0----R---1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NELEMX, NELEMY, the number of elements along the
%    X and Y directions.  The number of elements generated will be
%    2 * NELEMX * NELEMY.
%
%    Output, integer ELEMENT_NODE(10,2*NELEMX*NELEMY), the nodes that form
%    each element.
%
  element = 0;

  for j = 1 : nelemy
    for i = 1 : nelemx

      base = ( j - 1 ) * 3 * ( 3 * nelemx + 1 ) + 3 * i - 2;

      element = element + 1;

      element_node( 1,element) = base;
      element_node( 2,element) = base                          + 1;
      element_node( 3,element) = base                          + 2;
      element_node( 4,element) = base                          + 3;
      element_node( 5,element) = base +     ( 3 * nelemx + 1 ) + 2;
      element_node( 6,element) = base + 2 * ( 3 * nelemx + 1 ) + 1;
      element_node( 7,element) = base + 3 * ( 3 * nelemx + 1 );
      element_node( 8,element) = base + 2 * ( 3 * nelemx + 1 );
      element_node( 9,element) = base +     ( 2 * nelemx + 1 ) + 2;
      element_node(10,element) = base +     ( 2 * nelemx + 1 ) + 3;

      element = element + 1;

      element_node( 1,element) = base + 3 * ( 3 * nelemx + 1 ) + 3;
      element_node( 2,element) = base + 3 * ( 3 * nelemx + 1 ) + 2;
      element_node( 3,element) = base + 3 * ( 3 * nelemx + 1 ) + 1;
      element_node( 4,element) = base + 3 * ( 3 * nelemx + 1 );
      element_node( 5,element) = base + 2 * ( 3 * nelemx + 1 ) + 1;
      element_node( 6,element) = base +     ( 3 * nelemx + 1 ) + 2;
      element_node( 7,element) = base                          + 3;
      element_node( 8,element) = base +     ( 3 * nelemx + 1 ) + 3;
      element_node( 9,element) = base + 2 * ( 3 * nelemx + 1 ) + 3;
      element_node(10,element) = base + 2 * ( 3 * nelemx + 1 ) + 2;

    end
  end

  return
end


text:
% testGraphCutsScript2

if 0
    outdir = '../data';
    load([outdir '/rand_indices.mat']);
    load([outdir '/allimsegs2.mat']);
    load([outdir '/mcmcSuperpixelClassifier.mat']);
    load([outdir '/mcmcEdgeData.mat']);
    load([outdir '/mcmcEdgeClassifier.mat']);
end

alpha = [0:0.5:3];

for k = 1:numel(alpha);

    disp(num2str(alpha(k)))

    for cvf = 1:numel(cv_images)

        %disp([num2str(cvf) ': ' imsegs(cv_images(cvf)).imname])

        f = cv_images(cvf);
        c = ceil(cvf/50);

        [labv, labh] = testImageGraphCuts2(pg{cvf}, efeatures{f}, adjlist{f}, eclassifier, ecal{c}, alpha(k));

        nsp = numel(labv);
        pg2{cvf} = zeros(nsp, 7);
        lab = (labv==1) + (labv==2).*(1+labh) + 7*(labv==3);
        % set most likely label
        pg2{cvf}((lab-1)*nsp + [1:nsp]') = 1;
        % set so that h accuracy can be determined
        pg2{cvf}((labh)*nsp + [1:nsp]') = 0.5;
    end

    [vaccgc(k), haccgc(k)] = mcmcProcessResult(imsegs(cv_images), pg2);
    disp(num2str([vaccgc ; haccgc]))
end

author: zouchuhang

repo: LayoutNet


text:
function [newnode,newelem]=mergesurf(node,elem,varargin)
%
% [newnode,newelem]=mergesurf(node1,elem1,node2,elem2,...)
%
% merge two or more triangular meshes and split intersecting elements
%
% author: Qianqian Fang <q.fang at neu.edu>
%
% input:
%      node: node coordinates, dimension (nn,3)
%      elem: tetrahedral element or triangle surface (nn,3)
%
% output:
%      newnode: the node coordinates after merging, dimension (nn,3)
%      newelem: tetrahedral element or surfaces after merging (nn,4) or (nhn,5)
%
% note: you can call meshcheckrepair for the output newnode and
% newelem to remove the duplicated nodes or elements
%
% example:
%
%   [node1,face1,elem1]=meshabox([0 0 0],[10 10 10],1,1);
%   [node2,face2,elem2]=meshasphere([5 5 10],3,0.3,3);
%   [newnode,newface]=mergesurf(node1,face1,node2,face2);
%   plotmesh(newnode,newface,'x>5');
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

len=length(varargin);
newnode=node;
newelem=elem;
if(len>0 && mod(len,2)~=0)
   error('you must give node and element in pairs');
end
for i=1:2:len
   no=varargin{i};
   el=varargin{i+1};
   [newnode,newelem]=surfboolean(newnode,newelem,'all',no,el);
end

author: fangq


text:
function L = noiseExpectationLogLikelihood(noise, mu, varsigma, y);

% NOISEEXPECTATIONLOGLIKELIHOOD Return the expectation of the log likelihood.
% FORMAT
% DESC returns the expectation of the log likelihood for a gven noise model.
% ARG noise : the noise structure for which the expectation of the log
% likelihood is required.
% ARG mu : input mean locations for the likelihood.
% ARG varSigma : input variance locations for the likelihood.
% ARG y : target locations for the likelihood.
%
% SEEALSO : noiseParamInit, noiseLogLikelihood
%
% COPYRIGHT : Neil D. Lawrence, 2007

% NOISE

fhandle = str2func([noise.type 'NoiseExpectationLogLikelihood']);
L = fhandle(noise, mu, varsigma, y);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function embedding = optimize_layout(head_embedding, tail_embedding, ...
    head, tail, n_epochs, n_vertices, epochs_per_sample, a, b, ...
    gamma, initial_alpha, negative_sample_rate, verbose)
%OPTIMIZE_LAYOUT Improve an embedding using stochastic gradient descent to
% minimize the fuzzy set cross entropy between the 1-skeletons of the high
% dimensional and low dimensional fuzzy simplicial sets. In practice this
% is done by sampling edges based on their membership strength (with the
% (1-p) terms coming from negative sampling similar to word2vec). This
% function is only called if the UMAP method is 'MATLAB'.
%
% embedding = OPTIMIZE_LAYOUT(head_embedding, tail_embedding, head, tail,
% n_epochs, n_vertices, epochs_per_sample, a, b)
%
% Parameters
% ----------
% head_embedding: array of size (n_samples, n_components)
%     The initial embedding to be improved by SGD.
%
% tail_embedding: array of size (source_samples, n_components)
%     The reference embedding of embedded points. If not embedding new
%     previously unseen points with respect to an existing embedding this
%     is simply the head_embedding (again); otherwise it provides the
%     existing embedding to embed with respect to.
%
% head: array of size (n_1_simplices, 1)
%     The indices of the heads of 1-simplices with non-zero membership.
%
% tail: array of size (n_1_simplices, 1)
%     The indices of the tails of 1-simplices with non-zero membership.
%
% n_epochs: double
%     The number of training epochs to use in optimization.
%
% n_vertices: double
%     The number of vertices (0-simplices) in the dataset.
%
% epochs_per_samples: array of size (n_1_simplices, 1)
%     A double value of the number of epochs per 1-simplex. 1-simplices with
%     weaker membership strength will have more epochs between being sampled.
%
% a: double
%     Parameter of differentiable approximation of right adjoint functor.
%
% b: double
%     Parameter of differentiable approximation of right adjoint functor.
%
% gamma: double (optional, default 1)
%     Weight to apply to negative samples.
%
% initial_alpha: double (optional, default 1)
%     Initial learning rate for the SGD.
%
% negative_sample_rate: double (optional, default 5)
%     Number of negative samples to use per positive sample.
%
% verbose: boolean (optional, default false)
%     Whether to report information on the current progress of the algorithm.
%
% Returns
% -------
% embedding: array of size (n_samples, n_components)
%     The optimized embedding.
%
% See also: OPTIMIZE_LAYOUT2
%
%   AUTHORSHIP
%   Math Lead & Primary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Secondary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause

    if nargin < 13
        verbose = false;
        if nargin < 12
            negative_sample_rate = 5;
            if nargin < 11
                initial_alpha = 1;
                if nargin < 10
                    gamma = 1;
                end
            end
        end
    end

    dim = size(head_embedding, 2);
    n_1_simplices = size(epochs_per_sample, 1);
    same_embedding = isequal(head_embedding, tail_embedding);
    alpha = initial_alpha;
    ONES=ones(1, dim);
    BG2S=2*gamma* b*ONES;
    FOURS=4*ONES;
    ABNEG2=-2.0*a*b;
    BNEG1=b-1;

    epochs_per_negative_sample = epochs_per_sample / single(negative_sample_rate);
    epoch_of_next_negative_sample = epochs_per_negative_sample;
    epoch_of_next_sample = epochs_per_sample;
    if verbose
        fprintf('\t0/%d epochs done\n', int32(n_epochs));
    end
    for n = 1:n_epochs
        for i = 1:n_1_simplices
            if epoch_of_next_sample(i) <= n
                j = head(i);
                k = tail(i);

                current = head_embedding(j,:);
                other = tail_embedding(k,:);

                dist_squared = norm(current - other).^2;

                grad_coeff = (dist_squared > 0)*(ABNEG2*(dist_squared).^BNEG1)./ (a*dist_squared.^b + 1);
                grad_coeff(isnan(grad_coeff)) = 0;

                grad= max(-4, min(4, grad_coeff .* (current - other)));
                current = current + grad * alpha;

                epoch_of_next_sample(i) = epoch_of_next_sample(i) + epochs_per_sample(i);

                n_neg_samples = floor((single(n) - epoch_of_next_negative_sample(i)) / epochs_per_negative_sample(i));

                if same_embedding
                    other = other - grad * alpha;
                    head_embedding(k,:) = other;
                    tail_embedding(k,:) = other;
                end

                for p = 1:n_neg_samples
                    k = int32(randi(n_vertices));

                    other = tail_embedding(k,:);

                    dist_squared = norm(current - other).^2;

                    if same_embedding && j == k
                        continue
                    end

                    grad_coeff = (dist_squared > 0)*BG2S./(0.001 + dist_squared)./(a*dist_squared.^b + 1);
                    grad_coeff(isnan(grad_coeff)) = 0;
                    grad=(grad_coeff > 0).*max(-4, min(4, grad_coeff .* (current - other))) + ~(grad_coeff > 0).*FOURS;
                    current = current + grad * alpha;
                end

                head_embedding(j,:) = current;
                if same_embedding
                    tail_embedding(j,:) = current;
                end

                epoch_of_next_negative_sample(i) = epoch_of_next_negative_sample(i)+(n_neg_samples * epochs_per_negative_sample(i));
            end

        end
        alpha = initial_alpha * (1 - single(n)/single(n_epochs));

        progress_checkpoint = min(floor(n_epochs / 10), 50);

        if verbose && mod(n, progress_checkpoint) == 0
            fprintf('\t%d/%d epochs done\n', int32(n), int32(n_epochs));
        end
    end

    embedding = head_embedding;
 end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [f0,irms,df,amp]=f0track5(f0v,vrv,dfv,pwt,pwh,aav,shiftm,imgi)

%	F0 trajectory tracker
%	[f0,irms,df,amp]=f0track2(f0v,vrv,dfv,shiftm,imgi)
%		f0	: extracted F0 (Hz)
%		irms	: relative interfering energy in rms
%
%	f0v	: fixed point frequency vector
%	vrv	: relative interfering energy vector
%	dfv	: fixed point slope vector
%	pwt	: total power
%	pwh	: power in higher frequency range
%	aav	: amplitude list for fixed points
%	shiftm	: frame update period (ms)
%	imgi	: display indicator, 1: display on (default), 0: off
%
%	This is a very primitive and conventional algorithm.

%	coded by Hideki Kawahara
%	copyright(c) Wakayama University/CREST/ATR
%	10/April/1999 first version
%	17/May/1999 relative fq jump thresholding
%	01/August/1999 parameter tweeking
%   07/Dec./2002 waitbar was added
%   13/Jan./2005 bug fix on lines 58, 97 (Thanx Ishikasa-san)
%	30/April/2005 modification for Matlab v7.0 compatibility
%	10/Aug./2005 modified  by Takahashi on waitbar
%	10/Sept./2005 modified by Kawahara on waitbar

if nargin==7; imgi=1; end; %10/Sept./2005
vrv=sqrt(vrv);
[~,mm]=size(vrv);
mm=min(mm,length(pwt));

f0=zeros(1,mm);
irms=ones(1,mm);
df=ones(1,mm);
amp=zeros(1,mm);
von=0;
[mxvr,ixx]=min(vrv);
hth=0.12;	% highly confident voiced threshould (updated on 01/August/1999)
lth=0.9;	% threshold to loose confidence
bklm=100;	% back track length for voicing decision
lalm=10;	% look ahead length for silence decision
bkls=bklm/shiftm;
lals=lalm/shiftm;
htr=10*log10(pwh./pwt);

thf0j=0.04*sqrt(shiftm); %  4 % of F0 is the limit of jump
ii=1;
f0ref=0;
htrth=-2.0; % was -3 mod 2002.6.3
if imgi==1; hpg=waitbar(0,'F0 tracking'); end; % 07/Dec./2002 by H.K.%10/Aug./2005
while ii < mm+1
	if (von == 0) && (mxvr(ii)<hth) && (htr(ii)<htrth)
		von = 1;
		f0ref=f0v(ixx(ii),ii); % start value for search
		for jj=ii:-1:max(1,ii-bkls)
			[gomi,jxx]=min(abs((f0v(:,jj)-f0ref)/f0ref));
			gomi=gomi+(f0ref>10000)+(f0v(jxx,jj)>10000);
			if (((gomi>thf0j) || (vrv(jxx,jj)>lth) || (htr(jj)>htrth))&&(f0v(jxx,jj)<1000)) && htr(jj)>-18
%				disp(['break pt1 at ' num2str(jj)])
				break
			end;
			if (gomi>thf0j)
%				disp(['break pt2 at ' num2str(jj)])
				break
			end;

			f0(jj)=f0v(jxx,jj);
			irms(jj)=vrv(jxx,jj);
			df(jj)=dfv(jxx,jj);
			amp(jj)=aav(jxx,jj);
			f0ref=f0(jj);
		end;
		f0ref=f0v(ixx(ii),ii);
	end;
	if (f0ref>0) && (f0ref<10000)
		[gomi,jxx]=min(abs((f0v(:,ii)-f0ref)/f0ref));
	else
		gomi=10;
	end;
	if (von ==1) && (mxvr(ii)>hth)
		for jj=ii:min(mm,ii+lals)
			ii=jj;
			[gomi,jxx]=min(abs((f0v(:,ii)-f0ref)/f0ref));
			gomi=gomi+(f0ref>10000)+(f0v(jxx,ii)>10000);
			if (gomi< thf0j) && ((htr(ii)<htrth)||(f0v(jxx,ii)>=1000))
				f0(ii)=f0v(jxx,ii);
				irms(ii)=vrv(jxx,ii);
				df(ii)=dfv(jxx,ii);
				amp(ii)=aav(jxx,ii);
				f0ref=f0(ii);
			end;
			if (gomi>thf0j) || (vrv(jxx,ii)>lth) || ((htr(ii)>htrth)&&(f0v(jxx,ii)<1000))
				von = 0;f0ref=0;
				break
			end;
		end;
	elseif (von==1) && (gomi < thf0j)  && ((htr(ii)<htrth)||(f0v(jxx,ii)>=1000))
		f0(ii)=f0v(jxx,ii);
		irms(ii)=vrv(jxx,ii);
		df(ii)=dfv(jxx,ii);
		amp(ii)=aav(jxx,ii);
		f0ref=f0(ii);
	else
		von=0;
	end;
    if imgi==1; waitbar(ii/mm); end; %,hpg); % 07/Dec./2002 by H.K.%10/Aug./2005
	ii=ii+1;
end;
if imgi==1; close(hpg); end;%10/Aug./2005

author: HidekiKawahara


text:
function [rc] = calc_ratechange(a,time_as,step,mintime, maxtime,timestep)
    % function [rc] = calc_ratechange(a,time_as,step,mintime,maxtime,timestep);
    % ----------------------------------------------------------------
    % Determines ratechanges within aftershock sequences
    %
    % Input parameters:
    %   a           earthquake catalog
    %   time_as     delay times (days)
    %   step        number of quakes to determine forecast period
    %   mintime     tart time
    %   maxtime     maximal time at which Omori parameters are calculated
    %   timestep    Timesteps between learning periods
    %
    % Output parameters:
    %   rc      Matrix containing: time, absdiff, sigma, numreal, nummod,
    %           pval, pvalstd, cval, cvalstd, kval, kvalstd, t_forecast, fMc
    %
    % Info: Version for ZMAP, original function by ratechange.m by S. Neukomm May 27, 2002
    % J. Woessner
    % last update: 09.07.03

    if min(time_as) >= 1
        rc = [];
        return
    end

    fMc0 = mcwithtime(a,time_as,7);

    rc = [];
    time = mintime;
    while time <= maxtime

        % determination of magnitude of completeness
        if time < 7
            fMc = mcwithtime(a,time_as,time);
        else
            fMc = fMc0;
        end

        % estimation of Omori parameters
        l = time_as <= time & a.Magnitude >= fMc;
        [pval, pvalstd, cval, cvalstd, kval, kvalstd, loopout] = bruteboot(time_as(l));

        if isnan(pval) == 0
            nummod = step; % forecasted number of aftershocks
            if pval == 1
                pv = 1-10^(-6);
            else
                pv = pval;
            end

            t_forecast = (-nummod*(pv-1)/kval+(time+cval)^(1-pv))^(1/(1-pv))-cval; % forecast interval
            if isreal(t_forecast) == 1  &&  t_forecast > time
                if t_forecast > max(time_as)
                    return
                end

                l = time_as <= t_forecast & a.Magnitude >= fMc & time_as >= time;
                numreal = sum(l); % observed number of aftershocks
                absdiff = numreal-nummod;

                % calculate uncertainty sigma in forecasted number of aftershocks
                time1 = t_forecast; mpm1 = 1-pv; t1c = time1+cval; t0c = time+cval;
                sigma = (((-t1c^mpm1+t0c^mpm1)/(pv-1)*kvalstd)^2+...
                    (kval/(pv-1)*(-t1c^mpm1*mpm1/t1c+t0c^mpm1*mpm1/t0c)*cvalstd)^2+...
                    (kval/(pv-1)*(t1c^mpm1*log(t1c)+t1c^mpm1/(pv-1)-t0c^mpm1*log(t0c)-t0c^mpm1/(pv-1))*pvalstd)^2)^0.5;

                rc = [rc; time absdiff sigma numreal nummod pval pvalstd cval cvalstd kval kvalstd t_forecast fMc];
            end
        end
        time = time+timestep;
    end

author: CelsoReyes


text:
% A simple tutorial file to interface with RF
% Options copied from http://cran.r-project.org/web/packages/randomForest/randomForest.pdf

%run plethora of tests
clc
close all

%compile everything
if strcmpi(computer,'PCWIN') |strcmpi(computer,'PCWIN64')
   compile_windows
else
   compile_linux
end

total_train_time=0;
total_test_time=0;

%load the twonorm dataset
load data/twonorm

%modify so that training data is NxD and labels are Nx1, where N=#of
%examples, D=# of features

X = inputs';
Y = outputs;

[N D] =size(X);
%randomly split into 250 examples for training and 50 for testing
randvector = randperm(N);

X_trn = X(randvector(1:250),:);
Y_trn = Y(randvector(1:250));
X_tst = X(randvector(251:end),:);
Y_tst = Y(randvector(251:end));

% example 1:  simply use with the defaults
    model = classRF_train(X_trn,Y_trn);
    Y_hat = classRF_predict(X_tst,model);
    fprintf('\nexample 1: error rate %f\n',   length(find(Y_hat~=Y_tst))/length(Y_tst));

% example 2:  set to 100 trees
    model = classRF_train(X_trn,Y_trn, 100);
    Y_hat = classRF_predict(X_tst,model);
    fprintf('\nexample 2: error rate %f\n',   length(find(Y_hat~=Y_tst))/length(Y_tst));

% example 3:  set to 100 trees, mtry = 2
    model = classRF_train(X_trn,Y_trn, 100,2);
    Y_hat = classRF_predict(X_tst,model);
    fprintf('\nexample 3: error rate %f\n',   length(find(Y_hat~=Y_tst))/length(Y_tst));

% example 4:  set to defaults trees and mtry by specifying values as 0
    model = classRF_train(X_trn,Y_trn, 0, 0);
    Y_hat = classRF_predict(X_tst,model);
    fprintf('\nexample 4: error rate %f\n',   length(find(Y_hat~=Y_tst))/length(Y_tst));

% example 5: set sampling without replacement (default is with replacement)
    extra_options.replace = 0 ;
    model = classRF_train(X_trn,Y_trn, 100, 4, extra_options);
    Y_hat = classRF_predict(X_tst,model);
    fprintf('\nexample 5: error rate %f\n',   length(find(Y_hat~=Y_tst))/length(Y_tst));

% example 6: Using classwt (priors of classes)
    clear extra_options;
    extra_options.classwt = [1 1]; %for the [-1 +1] classses in twonorm
    % if you sort the labels in training and arrange in ascending order then
    % for twonorm you have -1 and +1 classes, with here assigning 1 to
    % both classes
    % As you have specified the classwt above, what happens that the priors are considered
    % also is considered the freq of the labels in the data. If you are
    % confused look into src/rfutils.cpp in normClassWt() function

    model = classRF_train(X_trn,Y_trn, 100, 4, extra_options);
    Y_hat = classRF_predict(X_tst,model);
    fprintf('\nexample 6: error rate %f\n',   length(find(Y_hat~=Y_tst))/length(Y_tst));

% example 7: modify to make class(es) more IMPORTANT than the others
    %  extra_options.cutoff (Classification only) = A vector of length equal to
    %                       number of classes. The 'winning' class for an observation is the one with the maximum ratio of proportion
    %                       of votes to cutoff. Default is 1/k where k is the number of classes (i.e., majority
    %                       vote wins).    clear extra_options;
    extra_options.cutoff = [1/4 3/4]; %for the [-1 +1] classses in twonorm
    % if you sort the labels in training and arrange in ascending order then
    % for twonorm you have -1 and +1 classes, with here assigning 1/4 and
    % 3/4 respectively
    % thus the second class needs a lot less votes to win compared to the first class

    model = classRF_train(X_trn,Y_trn, 100, 4, extra_options);
    Y_hat = classRF_predict(X_tst,model);
    fprintf('\nexample 7: error rate %f\n',   length(find(Y_hat~=Y_tst))/length(Y_tst));
    fprintf('   y_trn is almost 50/50 but y_hat now has %f/%f split\n',length(find(Y_hat~=-1))/length(Y_tst),length(find(Y_hat~=1))/length(Y_tst));

%  extra_options.strata = (not yet stable in code) variable that is used for stratified
%                       sampling. I don't yet know how this works.

% example 8: sampsize example
    %  extra_options.sampsize =  Size(s) of sample to draw. For classification,
    %                   if sampsize is a vector of the length the number of strata, then sampling is stratified by strata,
    %                   and the elements of sampsize indicate the numbers to be drawn from the strata.
    clear extra_options
    extra_options.sampsize = size(X_trn,1)*2/3;

    model = classRF_train(X_trn,Y_trn, 100, 4, extra_options);
    Y_hat = classRF_predict(X_tst,model);
    fprintf('\nexample 8: error rate %f\n',   length(find(Y_hat~=Y_tst))/length(Y_tst));

% example 9: nodesize
    %  extra_options.nodesize = Minimum size of terminal nodes. Setting this number larger causes smaller trees
    %                   to be grown (and thus take less time). Note that the default values are different
    %                   for classification (1) and regression (5).
    clear extra_options
    extra_options.nodesize = 2;

    model = classRF_train(X_trn,Y_trn, 100, 4, extra_options);
    Y_hat = classRF_predict(X_tst,model);
    fprintf('\nexample 9: error rate %f\n',   length(find(Y_hat~=Y_tst))/length(Y_tst));

% example 10: calculating importance
    clear extra_options
    extra_options.importance = 1; %(0 = (Default) Don't, 1=calculate)

    model = classRF_train(X_trn,Y_trn, 100, 4, extra_options);
    Y_hat = classRF_predict(X_tst,model);
    fprintf('\nexample 10: error rate %f\n',   length(find(Y_hat~=Y_tst))/length(Y_tst));

    %model will have 3 variables for importance importanceSD and localImp
    %importance = a matrix with nclass + 2 (for classification) or two (for regression) columns.
    %           For classification, the first nclass columns are the class-specific measures
    %           computed as mean decrease in accuracy. The nclass + 1st column is the
    %           mean decrease in accuracy over all classes. The last column is the mean decrease
    %           in Gini index. For Regression, the first column is the mean decrease in
    %           accuracy and the second the mean decrease in MSE. If importance=FALSE,
    %           the last measure is still returned as a vector.
    figure('Name','Importance Plots')
    subplot(2,1,1);
    bar(model.importance(:,end-1));xlabel('feature');ylabel('magnitude');
    title('Mean decrease in Accuracy');

    subplot(2,1,2);
    bar(model.importance(:,end));xlabel('feature');ylabel('magnitude');
    title('Mean decrease in Gini index');

    %importanceSD = The ?standard errors? of the permutation-based importance measure. For classification,
    %           a D by nclass + 1 matrix corresponding to the first nclass + 1
    %           columns of the importance matrix. For regression, a length p vector.
    model.importanceSD

% example 11: calculating local importance
    %  extra_options.localImp = Should casewise importance measure be computed? (Setting this to TRUE will
    %                   override importance.)
    %localImp  = a D by N matrix containing the casewise importance measures, the [i,j] element
    %           of which is the importance of i-th variable on the j-th case. NULL if
    %          localImp=FALSE.
    clear extra_options
    extra_options.localImp = 1; %(0 = (Default) Don't, 1=calculate)

    model = classRF_train(X_trn,Y_trn, 100, 4, extra_options);
    Y_hat = classRF_predict(X_tst,model);
    fprintf('\nexample 11: error rate %f\n',   length(find(Y_hat~=Y_tst))/length(Y_tst));

    model.localImp

% example 12: calculating proximity
    %  extra_options.proximity = Should proximity measure among the rows be calculated?
    clear extra_options
    extra_options.proximity = 1; %(0 = (Default) Don't, 1=calculate)

    model = classRF_train(X_trn,Y_trn, 100, 4, extra_options);
    Y_hat = classRF_predict(X_tst,model);
    fprintf('\nexample 12: error rate %f\n',   length(find(Y_hat~=Y_tst))/length(Y_tst));

    model.proximity

% example 13: use only OOB for proximity
    %  extra_options.oob_prox = Should proximity be calculated only on 'out-of-bag' data?
    clear extra_options
    extra_options.proximity = 1; %(0 = (Default) Don't, 1=calculate)
    extra_options.oob_prox = 0; %(Default = 1 if proximity is enabled,  Don't 0)

    model = classRF_train(X_trn,Y_trn, 100, 4, extra_options);
    Y_hat = classRF_predict(X_tst,model);
    fprintf('\nexample 13: error rate %f\n',   length(find(Y_hat~=Y_tst))/length(Y_tst));

% example 14: to see what is going on behind the scenes
%  extra_options.do_trace = If set to TRUE, give a more verbose output as randomForest is run. If set to
%                   some integer, then running output is printed for every
%                   do_trace trees.
    clear extra_options
    extra_options.do_trace = 1; %(Default = 0)

    model = classRF_train(X_trn,Y_trn, 100, 4, extra_options);
    Y_hat = classRF_predict(X_tst,model);
    fprintf('\nexample 14: error rate %f\n',   length(find(Y_hat~=Y_tst))/length(Y_tst));

% example 14: to see what is going on behind the scenes
%  extra_options.keep_inbag Should an n by ntree matrix be returned that keeps track of which samples are
%                   'in-bag' in which trees (but not how many times, if sampling with replacement)
%
    clear extra_options
    extra_options.keep_inbag = 1; %(Default = 0)

    model = classRF_train(X_trn,Y_trn, 100, 4, extra_options);
    Y_hat = classRF_predict(X_tst,model);
    fprintf('\nexample 15: error rate %f\n',   length(find(Y_hat~=Y_tst))/length(Y_tst));

    model.inbag

% example 16: getting the OOB rate. model will have errtr whose first
% column is the OOB rate. and the second column is for the 1-st class and
% so on
    model = classRF_train(X_trn,Y_trn);
    Y_hat = classRF_predict(X_tst,model);
    fprintf('\nexample 16: error rate %f\n',   length(find(Y_hat~=Y_tst))/length(Y_tst));

    figure('Name','OOB error rate');
    plot(model.errtr(:,1)); title('OOB error rate');  xlabel('iteration (# trees)'); ylabel('OOB error rate');

% example 17: getting prediction per tree, votes etc for test set
    model = classRF_train(X_trn,Y_trn);

    test_options.predict_all = 1;
    [Y_hat, votes, prediction_pre_tree] = classRF_predict(X_tst,model,test_options);
    fprintf('\nexample 17: error rate %f\n',   length(find(Y_hat~=Y_tst))/length(Y_tst));

repo: sr-metric


text:
%FASTNLMEANSDENOISINGMULTI  Modification of fastNlMeansDenoising function for colored images sequences
%
%     dst = cv.fastNlMeansDenoisingMulti(srcImgs, imgToDenoiseIndex, temporalWindowSize)
%     dst = cv.fastNlMeansDenoisingMulti(..., 'OptionName',optionValue, ...)
%
% ## Input
% * __srcImgs__ Input 8-bit (or 16-bit with 'L1' norm) 1-channel, 2-channel,
%   3-channel or 4-channel images sequence. All images should have the same
%   type and size.
% * __imgToDenoiseIndex__ Target image to denoise index in `srcImgs` sequence
%   (0-based index).
% * __temporalWindowSize__ Number of surrounding images to use for target
%   image denoising. Should be odd. Images from
%   `imgToDenoiseIndex - temporalWindowSize/2` to
%   `imgToDenoiseIndex + temporalWindowSize/2` from `srcImgs` will be used to
%   denoise `srcImgs{imgToDenoiseIndex}` image.
%
% ## Output
% * __dst__ Output image with the same size and type as `srcImgs` images.
%
% ## Options
% * __H__ Array of parameters regulating filter strength, either one parameter
%   applied to all channels or one value per channel in `dst`. Big `H` value
%   perfectly removes noise but also removes image details, smaller `H` value
%   preserves details but also preserves some noise. default [3]
% * __TemplateWindowSize__ Size in pixels of the template patch that is used
%   to compute weights. Should be odd. Recommended value 7 pixels. default 7
% * __SearchWindowSize__ Size in pixels of the window that is used to compute
%   weighted average for given pixel. Should be odd. Affect performance
%   linearly: greater `SearchWindowsSize` - greater denoising time.
%   Recommended value 21 pixels. default 21
% * __NormType__ Type of norm used for weight calculation. Can be either:
%   * __L2__ (default)
%   * __L1__
%
% Modification of cv.fastNlMeansDenoising function for images sequence where
% consecutive images have been captured in small period of time. For example
% video. This version of the function is for grayscale images or for manual
% manipulation with colorspaces. For more details see:
% [CiteSeerX](http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.131.6394)
%
% See also: cv.fastNlMeansDenoising, cv.fastNlMeansDenoisingColored,
%  cv.fastNlMeansDenoisingColoredMulti
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%kvskthick 'Exoskeleton by Thickening'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vskthick.pane file
%
% Parameters:
% InputFile: i 'Input Image ', required: 'input image'
% InputFile: str1 'Struct Element 1', required: 'Structuring Element'
% InputFile: str2 'Struct Element 2', required: 'Structuring Element'
% MultiChoice: r 'Struct El. Rotation Step', default: 1: 'Rotation step'
%    Choices are:
%   1: '1 '
%   2: '2 '
%   3: '4 '
% OutputFile: o 'Output Image ', required: 'resulting output image'
%
% Example: o = kvskthick({i, str1, str2}, {'i','';'str1','';'str2','';'r',1;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vskthick - Exoskeleton by Thickening
%
%  DESCRIPTION
% .I vskthick
% will do the exoskeleton by thickening of an input image by a set of structuring elements.
%
% The output image data type will be the same as that of the highest data type
% of the input image.
%
% Some useful pairs of structuring elements:
%
%                               0 0 0    1 1 1
%   Conditional Segmentation    0 1 0    0 0 0
%                               1 1 1    0 0 0
%
%                               0 0 0    1 0 0
%   Pseudo Convex Hull          0 1 1    1 0 0
%                               0 0 0    1 0 0
%
%                               1 0 0    0 0 0
%   Convex Hull                 1 0 0    0 1 0
%                               1 0 0    0 0 0
%
%                               0 0 0    1 1 1
%   Homotopic Pruning           0 1 0    1 0 1
%                               0 0 0    1 1 1
%
%
%
%  EXAMPLES
% vskthick -i ball.xv -str1 b272.str -str2 b273.str -r 1 -o outimage.xv
%
% Will do the exoskeleton by thickening of image 1 "ball.xv" by the structuring element b272.str and b273.str, rotating the struct. els. 1 step after each iteration, with the resulting image written to "outimage.xv".
%
%  "SEE ALSO"
% vcast1(1)
%
%  RESTRICTIONS
% .I vskthick
% can be defined for all data types supported by Khoros, but at the moment it has been implemented just for the bit and unsigned char types.
% The structuring elements are subsets of the 3x3 matrix and the origin is always at the center of this matrix.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993-1997 Junior Barrera, Roberto Lotufo.  All rights reserved.
%

function varargout = kvskthick(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kvskthick(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'str1', '__input';'str2', '__input';'r', 1;'o', '__output'};
maxval={0,0,0,0,0};
minval={0,0,0,0,0};
istoggle=[0,0,0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','InputFile','MultiChoice','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'vskthick"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%% Axes and Antipodal Symmetry
%
%% Directions vs. Axes
%
% In MTEX it is possible to consider three dimensional vectors either as
% directions or as axes. The key option to distinguish between both
% interpretations is *antipodal*.
%
% Consider a pair of vectors

v1 = vector3d(1,1,2);
v2 = vector3d(1,1,-2);

%%
% and plots them in a spherical projection

plot([v1,v2],'label',{'v_1','v_2'})

%%
% These vectors will appear either on the upper or on the lower hemisphere.
% In order to treat these vectors as axes, i.e. in order to assume
% antipodal symmetry - one has to use the keyword *antipodal*.

plot([v1,v2],'label',{'v_1','v_2'},'antipodal')

%%
% Now the direction *v_2* is identified with the direction *-v_2* which
% plots at the upper hemisphere.

%% The Angle between Directions and Axes
%
% As a consequence the angle between two axes *v1*, *v2* will always be the
% smallest angle between the directions *v1*, *v2* and *v1*, *-v2*, i.e. it
% will always be smaller than 90 degree. In the absence of antipodal
% symmetry we obtain

angle(v1,v2) / degree

%%
% whereas, if antipodal symmetry is assumed we obtain

angle(v1,v2,'antipodal') / degree

%% Antipodal Symmetry in Density Estimation
%
% Another example, where antipodal symmetry matters is
% <VectorsDensityEstimation.html density estimation>. For ordinary
% directions we obtain an arbitrary spherical function

v = vector3d.rand(100)
density = v.calcDensity;
plot(density)

%%
% Whereas, if antipodal symmetry is present the resulting density function
% will have antipodal symmetry as well

density = v.calcDensity('antipodal')
plot(density,'complete')

%% Antipodal Symmetry in Experimental Pole Figures
%
% Due to Friedel's law experimental pole figures always provide antipodal
% symmetry. One consequence of this fact is that MTEX plots pole figure
% data always on the upper hemisphere. Moreover if you annotate a certain
% direction to pole figure data, it is always interpreted as an axis, i.e.
% projected to the upper hemisphere if necessary

mtexdata dubna
CS = pf.CS;

% plot the first pole figure
plot(pf({1}))

% annotate a axis on the souther hemisphere
annotate(vector3d(1,0,-1),'labeled','backgroundColor','w')

%% Antipodal Symmetry in Recalculated Pole Figures
%
% However, in the case of pole figures calculated from an ODF antipodal
% symmetry is in general not present.

% some prefered orientation
o = orientation.byEuler(20*degree,30*degree,0,'ZYZ',CS);

% define an unimodal ODF
odf = unimodalODF(o);

% plot pole figures
plotPDF(odf,[Miller(1,2,2,CS),-Miller(1,2,2,CS)])

%%
% Hence, if one wants to compare calculated pole figures with experimental
% ones, one has to add antipodal symmetry.

plotPDF(odf,Miller(1,2,2,CS),'antipodal')

%% Antipodal Symmetry in Inverse Pole Figures
%
% The same reasoning as above holds true for inverse pole figures. If we
% look at complete, inverse pole figures they do not posses antipodal symmetry
% in general

plotIPDF(odf,[yvector,-yvector],'complete')

%%
% However, if we add the keyword antipodal, antipodal symmetry is enforced.

plotIPDF(odf,yvector,'antipodal','complete')

%%
% Notice how MTEX, automatically reduces the fundamental region of inverse
% pole figures in the case that antipodal symmetry is present.

plotIPDF(odf,yvector)

%%
plotIPDF(odf,yvector,'antipodal')

author: mtex-toolbox

repo: mtex

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


text:
function[C] = localCorrections(x,y,a,G,rq,tol,gradOpt)
% C = localCorrections(x,y,a,G,rq,rMax,tol)
% Local correction matrix for the EBD.
% inputs :
% - x and y : arrays of size N1x2 and N2x2 (clouds of points) inside the
% unit disk.
% - a : the cutoff parameter of the EBD method
% - G : the kernel
% - rq : the radial quadrature that has been computed for G
% - rMax : such that initially, the data was X = rMax*x, Y = rMax*y (before
% rescaling).
% - tol : required tolerance.
% output : the sparse correction matrix C.

N1 = size(x,1);
N2 = size(y,1);
[I,rxyTemp] = rangesearch(y,x,a*1.05);
jdx = cell2mat(I')';
rxy = cell2mat(rxyTemp')';
idx = zeros(size(jdx));
sp_ind = 1;
for x_ind=1:length(I)
    idx(sp_ind:(sp_ind+length(I{x_ind})-1)) = x_ind;
    sp_ind = sp_ind + length(I{x_ind});
end
% Save memory
clear I;
NCI = length(rxy);

if NCI ~= 0

    if gradOpt
        exact_Interactions = G.evalDer(rxy); % Exact local interactions
        C1 = sum(abs(rq.alpha0.*Cp(rq.rho)).*rq.rho.^3); % Bound for the second derivative.
        Ninterp = fix(sqrt(C1)*a/sqrt(8*tol))+10; % Guarantees interpolation error < tol.
        xinterp = linspace(0,a*1.1,Ninterp);
        yinterp = rq.evalDer(xinterp);
        radial_quadratureNear0 = interp1(xinterp,yinterp,rxy); % we remove the radial
    else
        exact_Interactions = G.eval(rxy); % Exact local interactions
        C1 = sum(abs(rq.alpha0.*Cp(rq.rho)).*rq.rho.^2); % Bound for the second derivative.
        Ninterp = fix(sqrt(C1)*a/sqrt(8*tol))+10; % Guarantees interpolation error < tol.
        xinterp = linspace(0,a*1.1,Ninterp);
        yinterp = rq.eval(xinterp);
        radial_quadratureNear0 = interp1(xinterp,yinterp,rxy); % we remove the radial
    end

    % quadrature contribution (using interpolation to avoid evaluating at
    % all points).

    C_val = exact_Interactions - radial_quadratureNear0;
    C = sparse(idx,jdx,C_val,N1,N2);

    if gradOpt
        Tx = sparse(idx,jdx,y(jdx,1) - x(idx,1),N1,N2);
        Ty = sparse(idx,jdx,y(jdx,2) - x(idx,2),N1,N2);
        idx = idx(rxy> 1e-12);
        jdx = jdx(rxy> 1e-12);
        rxy = rxy(rxy> 1e-12);
        R =  sparse(idx,jdx,1./rxy,N1,N2);
        Cx = C.*R.*Tx;
        Cy = C.*R.*Ty;
        C = {Cx,Cy};
    end
else
    % No close interactions
    if gradOpt
        C = {sparse(N1,N2),sparse(N1,N2)};
    else
        C = sparse(N1,N2); % all zeros
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% Local Regression and Likelihood, Figure 9.1.
%
% Hardle's Motorcycle accelaration dataset. Just a scatterplot!
%
% Author: Catherine Loader

load mcyc;
figure('Name','fig9_1: motorcycle scatterplot');
plot(time,accel,'.');
xlabel('Time');
ylabel('Acceleration');

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function newLabel=GetLabel(trnData,tstData,trnLabel)

% Shared by algorithms.
% copyright @ guan naiyang

dist=EuDist2(tstData',trnData',0);
[junk, sortedIdx]=sort(dist,2);
newLabel=trnLabel(sortedIdx);
clear dist;

return;

author: andrewssobral

repo: lrslibrary


text:
function [a] = save_warper(string, proposal,mask)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

 a=0;
 save(string,'proposal','mask');
end

author: xjqicuhk

repo: SIMS

save_path: github-repos/MATLAB/xjqicuhk-SIMS


text:
function f=comp_idgtreal(coef,g,a,M,lt,phasetype)
%COMP_IDGTREAL  Compute IDGTREAL
%   Usage:  f=comp_idgtreal(c,g,a,M,lt,phasetype);
%
%   Input parameters:
%         c     : Array of coefficients.
%         g     : Window function.
%         a     : Length of time shift.
%         M     : Number of modulations.
%         lt    : lattice type
%   Output parameters:
%         f     : Signal.
%

%   AUTHOR : Peter L. Søndergaard.
%   TESTING: TEST_DGT
%   REFERENCE: OK

N=size(coef,2);
L=N*a;

if lt(2)==1
    f = comp_isepdgtreal(coef,g,L,a,M,phasetype);
else
    % Quinqux lattice
    f=comp_inonsepdgtreal_quinqux(coef,g,a,M);
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% function wp = WT_inv(x,L,N,wavename)
%
% Perform a discrete wavelet transform of the data and of the shifted
% version of the data for every decomposition level up to N-L. The shift
% helps the reconstruction
%
%
% INPUTS:
% x:            1D signal on which to perform the wavelet transform
% L:            Lowest wavelet scale used in the analysis
% N:            Number of wavelet levels in which the signal can be decomposed
% wavename:     name of the wavelet used for the decomposition
%
%
% OUTPUTS:
% wp:           Wavelet decomposition coefficient matrix ( #of time
%               points x #of levels+1). The first column contains the last approximation while
%               from the second to the end the details at all levels
%
%
% LOG:
% 10/17/2012 by S. Brigadoi

function wp = WT_inv(x,L,N,wavename)

D = N-L;
n = length(x);
wp = zeros(n,D+1);
dwtmode('per');  % set the wavelet mode to periodization

wp(:,1) = x';
for d=0:(D-1)
    n_blocks = 2^d; % number of blocks in the level
    l_blocks = n/n_blocks; % length of the blocks in the level
    for b=0:(2^d-1)
        s = wp(b*l_blocks+1:b*l_blocks+l_blocks,1)'; % first time take signal, from the second the approximation
        s_shift = [s(end) s(1:end-1)]; % create a shift version of the block

        [cA,cD] = dwt(s,wavename);  % discrete wavelet transform
        [cA_shift,cD_shift] = dwt(s_shift,wavename); % discrete wavelet transform of the shifted version

        wp(b*l_blocks+1:b*l_blocks+l_blocks/2,1) = cA;
        wp(b*l_blocks+l_blocks/2+1:b*l_blocks+l_blocks,1) = cA_shift;

        wp(b*l_blocks+1:b*l_blocks+l_blocks/2,d+2) = cD;
        wp(b*l_blocks+l_blocks/2+1:b*l_blocks+l_blocks,d+2) = cD_shift;
    end
end

author: BUNPC


text:
%kkasc2val 'Import ASCII Data File to Value Segment'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kasc2val.pane file
%
% Parameters:
% InputFile: i1 'ASCII Input File', required: 'ASCII input data file containing new data segment'
% InputFile: i2 'Data Object Input', optional: 'data object input to insert the new segment into'
% Integer: so 'Start Offset ', default: 0: 'starting point in ASCII file to start importing at'
% Integer: sf 'Skip Factor ', default: 0: 'number of points in ASCII file to skip between reads'
% Integer: wsize ' ', default: 1: 'width size must be > 0'
% Integer: hsize ' ', default: 1: 'height size must be > 0'
% Integer: dsize ' ', default: 1: 'depth size must be > 0'
% Integer: tsize ' ', default: 1: 'time size must be > 0'
% Integer: esize ' ', default: 1: 'elements size must be > 0'
% OutputFile: o 'Output Object', required: 'output data object'
%
% Example: o = kkasc2val({i1, i2}, {'i1','';'i2','';'so',0;'sf',0;'wsize',1;'hsize',1;'dsize',1;'tsize',1;'esize',1;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kasc2val - Import ASCII Data File to Value Segment
%
%  DESCRIPTION
% This routine will import an ASCII data file and insert it into a data
% object if it is specified.  If there is no source object to insert the
% new segment into, the data object specified as the output file will be
% created.  The new segment will be inserted into this data object in the
% value segment.  If a source object was specified, then all of its
% data segments and the global object attributes will be copied to the
% destination.
% The index order is specified by setting the -windex, -hindex, -dindex,
% -tindex and -eindex flags.  They can not be set to the same value.
% The start offset is the point to start importing in the
% ASCII data.  All points before this one are ignored.
% The skip factor is the number of points to skip in the ASCII
% data between reads.  If the size of the new segment's lowest
% index is 40 and the skip factor is 5, then this routine will
% read in 40 points skip 5 and read in 40 more, etc. until the
% entire ASCII file read in.
% Complex values in ASCII are represented by number pairs, surrounded
% by parenthesis and separated by a comma.  For example, (8, 3) represents
% the complex number 8+3i.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% kasc2mask, kasc2loc, kasc2time, kasc2map, kimportasc
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkasc2val(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkasc2val(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'i2', '__input';'so', 0;'sf', 0;'wsize', 1;'hsize', 1;'dsize', 1;'tsize', 1;'esize', 1;'o', '__output'};
maxval={0,1,1,1,2,2,2,2,2,0};
minval={0,1,1,1,2,2,2,2,2,0};
istoggle=[0,1,1,1,1,1,1,1,1,0];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','Integer','Integer','Integer','Integer','Integer','Integer','Integer','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kasc2val"  '],Inputs,narglist);


text:
function s = cardinal_sin ( j, m, n, t )

%*****************************************************************************80
%
%% CARDINAL_SIN evaluates the J-th cardinal sine basis function.
%
%  Discussion:
%
%    The base points are T(I) = pi * I / ( M + 1 ), 0 <= I <= M + 1.
%    Basis function J is 1 at T(J), and 0 at T(I) for I /= J
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    John Boyd,
%    Exponentially convergent Fourier-Chebyshev quadrature schemes on
%    bounded and infinite intervals,
%    Journal of Scientific Computing,
%    Volume 2, Number 2, 1987, pages 99-109.
%
%  Parameters:
%
%    Input, integer J, the index of the basis function.
%    0 <= J <= M + 1.
%
%    Input, integer M, indicates the size of the basis set.
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real T(N), one or more points in [0,pi] where the
%    basis function is to be evaluated.
%
%    Output, real S(N), the value of the function at T.
%
  tj = pi * j / ( m + 1 );

  s = ( -1.0 ) ^ mod ( j + 1, 2 ) ...
    * sin ( tj ) ...
    * sin ( ( m + 1 ) * t ) ...
    / ( m + 1 ) ...
    ./ ( cos ( t ) - cos ( tj ) );

  i = find ( abs ( t - tj ) < eps );
  s(i) = 1.0;

  return
end

author: johannesgerer


text:
% Demo script for fitpath.m
%
% Example:
%   demo
%
% Will bring up the demonstration GUI. Click on the button, browse for a file with
% a long path and see it truncated in the GUI. Don't worry it won't do anything
% with the actual file.

%=================================================================================
%
%        AUTHOR: Erik Newton (Newtek Software Ltd)
%
%           WEB: www.newteksoftware.co.uk
%
%	Copyright Newtek Software Ltd 2006
%
% This function is free for all use, but please leave in my credit.
% Check out more tutorials & functions at our website.
%
%==================================================================================

h_figure = figure(...
   'Name','fitpath demo', ...
   'MenuBar','none', ...
   'NumberTitle','Off', ...
   'Position',[50 50 260 30], ...
   'Resize','Off');

LABEL_BACK = get(h_figure,'Color');

% Init filename
uicontrol(h_figure, ...
   'BackgroundColor',LABEL_BACK, ...
   'Position', [5 5 30 18], ...
   'Horizontal', 'Right', ...
   'String','File:', ...
   'Style','text');
h_filename = uicontrol(h_figure, ...
   'BackgroundColor',LABEL_BACK, ...
   'Position',[35 5 200 18], ...
   'Horizontal', 'Left', ...
   'Style','text', ...
   'Tag','InitFilename');
uicontrol(h_figure, ...
   'Position', [235 5 20 18], ...
   'Horizontal', 'Center', ...
   ...'Callback', '[init_file,new_path]=uigetfile; set(h_filename,''string'',fullfile(new_path,init_file));', ...
   'Callback', '[init_file,new_path]=uigetfile; fitpath(h_filename,fullfile(new_path,init_file));', ...
   'String','...', ...
   'TooltipString', 'Select file.', ...
   'Style','PushButton');

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

author: Sable


text:
function plot_distribution(QSM,fig,rela,cumu,dis,dis2,dis3,dis4)

% ---------------------------------------------------------------------
% PLOT_DISTRIBUTION     Plots the specified distribution(s) in the
%                           "treedata" field of the QSM structure array.
%
% Version 1.1.0
% Latest update     3 May 2022
%
% Copyright (C) 2020-2022 Pasi Raumonen
% ---------------------------------------------------------------------
%
% Inputs:
% QSM       The output of treeqsm function, may contain multiple models if
%               only one distribution. If multiple distributions are plotted,
%               then only one model.
% fig       Figure number
% rela      If rela = 1, then plots relative values (%), otherwise plots
%               absolute values
% cumu      If cumu = 1, then plot cumulative distribution
% dis       Distribution to be plotted, string name, e.g. 'VolCylDia'.
%               The name string is the one used in the "treedata"
% dis2      Optional, Second distribution to be plotted. Notice with more
%               than one distribution, only one model.
% dis3      Optional, Third distribution to be plotted
% dis4      Optional, Fourth distribution to be plotted
% ---------------------------------------------------------------------

% Changes from version 1.0.0 to 1.1.0, 3 May 2022:
% 1) Added new input "cum" for plottig the distributions as cumulative.
% 2) Added return if distributions are empty or all zero

% Generate strings for title, xlabel and ylabel:
if strcmp(dis(1:3),'Vol')
  str = 'volume';
  ylab = 'Volume (L)';
elseif strcmp(dis(1:3),'Are')
  str = 'area';
  ylab = 'Area (m^2)';
elseif strcmp(dis(1:3),'Len')
  str = 'length';
  ylab = 'Length (m)';
elseif strcmp(dis(1:3),'Num')
  str = 'number';
  ylab = 'Number';
end

if strcmp(dis(end-2:end),'Dia')
  str2 = 'diameter';
  xlab = 'diameter (cm)';
elseif strcmp(dis(end-2:end),'Hei')
  str2 = 'height';
  xlab = 'height (m)';
elseif strcmp(dis(end-2:end),'Ord')
  str2 = 'order';
  xlab = 'order';
elseif strcmp(dis(end-2:end),'Ang')
  str2 = 'angle';
  xlab = 'angle (deg)';
elseif strcmp(dis(end-2:end),'Azi')
  str2 = 'azimuth direction';
  xlab = 'azimuth direction (deg)';
elseif strcmp(dis(end-2:end),'Zen')
  str2 = 'zenith direction';
  xlab = 'zenith direction (deg)';
end

% Collect the distributions
if nargin == 5
  % Multiple QSMs, one and the same distribution
  m = max(size(QSM));
  D = QSM(1).treedata.(dis);
  n = size(D,2);
  for i = 2:m
    d = QSM(i).treedata.(dis);
    k = size(d,2);
    if k > n
      n = k;
      D(m,n) = 0;
      D(i,1:n) = d;
    elseif k < n
      D(i,1:k) = d;
    else
      D(i,:) = d;
    end
  end
  D = D(:,1:n);
else
  % One QSM, multiple distributions of the same type
  % (e.g. diameter distributions: 'NumCylDia', 'VolCylDia' and 'LenCylDia')
  m = nargin-4;
  D = QSM.treedata.(dis);
  n = size(D,2);
  if n == 0 || all(D == 0)
    return
  end
  for i = 2:m
    if i == 2
      D(m,n) = 0;
      D(i,:) = QSM.treedata.(dis2);
    elseif i == 3
      D(i,:) = QSM.treedata.(dis3);
    else
      D(i,:) = QSM.treedata.(dis4);
    end
  end
end

if rela
  % use relative value
  for i = 1:m
    D(i,:) = D(i,:)/sum(D(i,:))*100;
  end
  ylab = 'Relative value (%)';
end

if cumu
  % use cumulative distribution
  D = cumsum(D,2);
end

% Generate the bar plot
figure(fig)
if strcmp(dis(end-3:end),'hAzi') || strcmp(dis(end-3:end),'1Azi') || strcmp(dis(end-2:end),'Azi')
  bar(-170:10:180,D')
elseif strcmp(dis(end-2:end),'Zen') || strcmp(dis(end-2:end),'Ang')
  bar(10:10:10*n,D')
else
  bar(1:1:n,D')
end

% Generate the title of the plot
if strcmp(dis(end-2:end),'Ord') && ~strcmp(dis(1:3),'Num')
  tit = ['Branch ',str,' per branching order'];
elseif strcmp(dis(end-2:end),'Ord')
  tit = 'Number of branches per branching order';
elseif strcmp(dis(1:3),'Num')
  tit = ['Number of branches per ',str2,' class'];
elseif strcmp(dis(end-3),'h') || strcmp(dis(end-3),'1')
  tit = ['Branch ',str,' per ',str2,' class'];
else
  tit = ['Tree segment ',str,' per ',str2,' class'];
end
n = nargin;
if n > 5
  if ~strcmp(dis(1:3),dis2(1:3))
    if strcmp(dis(4),'C')
      tit = 'Tree segment distribution';
    else
      tit = 'Branch distribution';
    end
  elseif n > 6
    if ~strcmp(dis(1:3),dis3(1:3))
      if strcmp(dis(4),'C')
        tit = 'Tree segment distribution';
      else
        tit = 'Branch distribution';
      end
    elseif n > 7
      if ~strcmp(dis(1:3),dis4(1:3))
        if strcmp(dis(4),'C')
          tit = 'Tree segment distribution';
        else
          tit = 'Branch distribution';
        end
      end
    end
  end
end
title(tit)

% Generate the x-axis label
if strcmp(dis(end-5:end-3),'Cyl')
  xlab = ['Cylinder ',xlab];
else
  xlab = ['Branch ',xlab];
end
xlabel(xlab)

% Generate the y-axis label
ylabel(ylab);

% Tight axes and grid lines
axis tight
grid on

m = max(size(QSM));
% Add legends, if needed
if m > 1
  L = cell(m,1);
  for i = 1:m
    L{i} = ['model',num2str(i)];
  end
  legend(L,'location','best')
elseif nargin > 5
  m = nargin-4;
  L = cell(m,1);
  for i = 1:m
    if i == 1
      L{i} = dis(1:end-3);
    elseif i == 2
      L{i} = dis2(1:end-3);
    elseif i == 3
      L{i} = dis3(1:end-3);
    else
      L{i} = dis4(1:end-3);
    end
  end
  legend(L,'location','best')
end

author: InverseTampere

repo: TreeQSM


text:
clear;
clc;
close all;
format

%% 参考: https://x-io.co.uk/oscillatory-motion-tracking-with-x-imu/

% 利用速度，位置高通滤波来获得短时位移曲线

%load filter_vel.mat
% load example_ins4.mat

data = ch_data_import('data_20210609_163153.csv');
acc = data.imu.acc*9.8;
gyr = deg2rad(data.imu.gyr);
%% 读取数据
acc = acc';
gyr = gyr';

q= [1 0 0 0]';

Fs = 100;
dt = 1 / Fs;
N = length(acc);

fprintf("共%d数据，用时:%.3fs\n", N, N/Fs);

% 惯导解算, 初始化
p = zeros(3, 1);
v = zeros(3, 1);
linVel = zeros(N,3);

for i=1:N
    [p , v , q] = ch_nav_equ_local_tan(p, v, q, acc(i,:)', gyr(i,:)', 1 / Fs, [0, 0, -9.8]');

    % 获得惯性系下加速度并扣除重力
    linAcc(i,:) = ch_qmulv(q, acc(i,:));
    linAcc(i,3) = linAcc(i,3) -9.8;

%   linVel(i,:) = v;
end

% 速度积分

for i = 2:N
    linVel(i,:) = linVel(i-1,:) + linAcc(i,:) * dt;
end

order = 1;
filtCutOff = 0.01;
[b, a] = butter(order, (2*filtCutOff)/(1/dt), 'high');
%linVelHP = filtfilt(b, a, linVel);
linVelHP = filter(b, a, linVel);

% 位置积分
linPos = zeros(size(linVelHP));
for i = 2:length(linVelHP)
    linPos(i,:) = linPos(i-1,:) + linVelHP(i,:) * dt;
end

order = 1;
filtCutOff = 0.5;
[b, a] = butter(order, (2*filtCutOff)/(1/dt), 'high');
%linPosHP = filtfilt(b, a, linPos);
linPosHP = filter(b, a, linPos);

%% Plot
figure('NumberTitle', 'off', 'Name', '速度');
subplot(2,1,1);
plot(linVel);
title('速度');
legend('X', 'Y', 'Z');
subplot(2,1,2);
plot(linVelHP);
title('HP后速度');
legend('X', 'Y', 'Z');

figure('NumberTitle', 'off', 'Name', '位置');
subplot(2,1,1);
plot(linPos);
title('位置');
legend('X', 'Y', 'Z');
subplot(2,1,2);
plot(linPosHP);
title('HP后位置');
legend('X', 'Y', 'Z');

figure('NumberTitle', 'off', 'Name', '原始数据');
subplot(2,1,1);
plot(acc);
legend("X", "Y", "Z");
title("加速度");
subplot(2,1,2);
plot(gyr);
title("陀螺");
legend("X", "Y", "Z");

figure('NumberTitle', 'off', 'Name', '3D位置');
plot3(linPosHP(1,1), linPosHP(1,2), linPosHP(1,3), '-ks');
hold on;
plot3(linPosHP(:,1), linPosHP(:,2), linPosHP(:,3), '.b');
axis equal
xlabel('X(m)');  ylabel('Y(m)');   zlabel('Z(m)');
title('3D位置');
legend('起始', '3D');

linPosHP(end,:)

% SamplePlotFreq = 4;
%
% SixDOFanimation(linPosHP, R, ...
%                 'SamplePlotFreq', SamplePlotFreq, 'Trail', 'Off', ...
%                 'Position', [9 39 400 400], ...
%                 'AxisLength', 0.1, 'ShowArrowHead', false, ...
%                 'Xlabel', 'X (m)', 'Ylabel', 'Y (m)', 'Zlabel', 'Z (m)', 'ShowLegend', false, 'Title', 'Unfiltered',...
%                 'CreateAVI', false, 'AVIfileNameEnum', false, 'AVIfps', ((1/dt) / SamplePlotFreq));

author: yandld