text:
function [c,cl] = cluster_nmdsfig(cl,fldname,colm,varargin)
% function [c,cl] = cluster_nmdsfig(cl,fldname,colm,[behavior],[ancova var],[covs of no interest],[overlay],[niter])
%
% Makes a multidim scaling figure from clusters
% Plus some other plots
%
% Example:
% c = cluster_nmdsfig(cl,'BARPLOT.data',2);
%
% c is a structure with lots of info
% cl has names appended, if not entered previously
%
% behavioral input is of three types:
% outcome       covariate of interest (1 only)
% ancova_var        categorical ANCOVA coding var (2-levels, one only);
%                   generally not of interest
%                   IF ancova_var is entered, "PPI" of group x brain
%                   covariance interaction will be computed
% covs of no int    of no interest; will be removed from brain and behavior
%                   before performing ANCOVA (for simplicity); 1 or more OK
%
% Tor Wager
%
% Examples:
%
% Use average data from cl.timeseries and correlate with EXPT.behavior.
% Covary out effects of order before running models; no ANCOVA
% [c,cl] = cluster_nmdsfig(cl,'timeseries',1,EXPT.behavior',[],EXPT.order2);
%
% More examples:
% c =
% cluster_nmdsfig(cl,'CONTRAST.yadj',1,repmat(EXPT.cov(:,1),4,1),[],[],EXPT.overlay,1000);
% [c,cl] = cluster_nmdsfig(cl,'CONTRAST.data',1,R.X(:,1),[],R.X(:,2),EXPT.overlay,200);
% [c,cl] = cluster_nmdsfig(cl,'timeseries',1,SETUP.Y,[],SETUP.X,EXPT.overlay,1000);

% subdirectory: prompt
dosubdir = input('Create a subdirectory for png images? (type name or return to use current dir) ','s');
if ~isempty(dosubdir)
    mkdir(dosubdir)
    cd(dosubdir)
end

diary cluster_nmdsfig_output.txt
doranks = 1;
addbtwn = 0;

if length(varargin) > 0, c.outcome = varargin{1}; else  c.outcome = [];end
if length(varargin) > 1, c.ancova_codes = varargin{2}; else  c.ancova_codes = [];end
if length(varargin) > 2, c.covs_nointerest = varargin{3}; else  c.covs_nointerest = [];end
if length(varargin) > 3, overlay = varargin{4}; else  overlay = [];end
if length(varargin) > 4, niter = varargin{5}; else  niter = 1000;end

cl(1).outcome = c.outcome; cl(1).ancova_codes = c.ancova_codes; cl(1).covs_nointerest=c.covs_nointerest;

c.X = [c.outcome, c.covs_nointerest];
contrasts = [];

c.colors = {'ro' 'go' 'bo' 'yo' 'co' 'mo' 'ko' 'r^' 'g^' 'b^' 'y^' 'c^' 'm^' 'k^'};

if ~isfield(c,'btwn_names') || isempty(c.btwn_names)
    c.btwn_names = {'Behav'};
end

% -----------------------------------------------------------------------
% get data from selected field
% remove covariates of no interest from predictors and data
% -----------------------------------------------------------------------

c = nmdsfig_tools('get_data_matrix',c,cl,fldname,colm,contrasts,doranks);

if addbtwn
    % add between-subjects behavioral scores to dat; just to get stim coords; remove later...
    c.dat = [c.outcome c.dat];
end

% -----------------------------------------------------------------------
% get names
% -----------------------------------------------------------------------

cl = cluster_names(cl,1);
for i = 1:length(cl), c.names{i} = cl(i).shorttitle; end

% -----------------------------------------------------------------------
% get correlations, significance, and number of dimensions
% If groups are entered and ANCOVA is specified, handle that
% Also get c.D, dissimilarity matrix for MDS
% -----------------------------------------------------------------------

c = nmdsfig_tools('get_correlations',c);

% MDS decomposition and plotting
% -----------------------------------------------------------------------
% get number of dimensions that carry most of the variance.
% This is to simplify the space and make it less sparse --
% expected to give better clustering results.
%
% Reduce to distances in ndims space, check this against actual distances
% -----------------------------------------------------------------------

% graphic check on number of dimensions, and MDS
[c.GroupSpace,c.obs,c.implied_dissim] = shepardplot(c.D,[]);
drawnow
c.ndims = size(c.GroupSpace,2); % set here; graphical decision within shepardplot

nclust = 2:(size(c.GroupSpace,1) / 2)+1;     % choose number of clusters to test

save_figure('cluster_nmdsfig_shepardplot');

if addbtwn
    % get coordinates for behavior along with the others
    fprintf(1,'Using between-subjects covariates in MDS model.\n');

    % save coordinates in MDS space for btwn data
    c.btwn_coords = c.GroupSpace(1,:);

    % remove from GroupSpace
    c.GroupSpace = c.GroupSpace(2:end,:);

    % remove from r
    c.r = c.r(2:end,2:end);
end

% -----------------------------------------------------------------------
% testcluster:
%
% Permute objects in space, get null hypothesis cluster quality, and test
% observed cluster solution against this.
% -----------------------------------------------------------------------

% pval is p-values for each number of clusters tested
% classes = class ("network") assignments for best clustering solution
[c.ClusterSolution.pvals,c.ClusterSolution.classes, c.ClusterSolution.classnames ...
    c.ClusterSolution.X,c.ClusterSolution.include,c.ClusterSolution.names]= ...
    testclustnew(c.GroupSpace,nclust,c.ndims,niter,c.names,'keep','average');

if ~any(c.ClusterSolution.pvals < .05)
    % nothing significant
    c.ClusterSolution.classes = ones(size(c.ClusterSolution.classes));
end

%[bestpval,bestXc,bestnames,bestX,where,clustnames]=testcluster(X,clust,[r],[nperm],[names],[remove],[linkagetype]);

save_figure('cluster_nmdsfig_testcluster');

% -----------------------------------------------------------------------
% factor scores with reduced dimensions
% -----------------------------------------------------------------------
r = c.r;
[tmp c.eigenvalues] = cmdscale(c.D);
sqL = sqrt(diag(c.eigenvalues(1:c.ndims)));
B = inv(r) * c.GroupSpace(:,1:c.ndims) * sqL;   % inv(r) * factor loading matrix A = VsqrtL
c.factor_scores = c.dat * B;

% -----------------------------------------------------------------------
% Stepwise regression of dimensions on behavior
% -----------------------------------------------------------------------

if ~isempty(c.outcome)
    fprintf(1,'\n');
    fprintf(1,'Prediction of behavior with component scores (from classic MDS)\n');

    DATA.INDSCAL.STEPWISE = stepwise_tor(c.factor_scores,c.outcome);
    fprintf(1,'\n');
end

% -----------------------------------------------------------------------
% % plot nmds figures of this
% -----------------------------------------------------------------------
dosideplots = 0;
figname = nmdsfig_tools('nmdsfig_plot',c, addbtwn, dosideplots, 'nofill');

save_figure(figname);

% -----------------------------------------------------------------------
% separate plots of high/low behavior
% -----------------------------------------------------------------------
if ~isempty(c.outcome)  % plot interactions between beh and regional covariance
    create_wide_figure;
    subplot(1,2,1);
    [r,p]=corrcoef(c.dat(c.outcome>median(c.outcome),:)); r(r>.9999) = 1; sig = sign(r) .* (p < .05);
    nmdsfig(c.GroupSpace,'classes',c.ClusterSolution.classes,'names',c.names,'sig',sig);
    title('High behavior')

    subplot(1,2,2)
    [r,p]=corrcoef(c.dat(c.outcome<=median(c.outcome),:)); r(r>.9999) = 1; sig = sign(r) .* (p < .05);
    nmdsfig(c.GroupSpace,'classes',c.ClusterSolution.classes,'names',c.names,'sig',sig);
    title('Low behavior')
    drawnow
end

if ~isempty(c.ancova_codes)  % plot interactions between beh and regional covariance
    tor_fig(1,3); subplot(1,3,1);
    [r,p]=corrcoef(c.dat(c.ancova_codes>median(c.ancova_codes),:)); r(r>.9999) = 1; sig = sign(r) .* (p < .05);
    nmdsfig(c.GroupSpace,'classes',c.ClusterSolution.classes,'names',c.names,'sig',sig);
    title('High ANCOVA group')

    subplot(1,3,2)
    [r,p]=corrcoef(c.dat(c.ancova_codes<median(c.ancova_codes),:)); r(r>.9999) = 1; sig = sign(r) .* (p < .05);
    nmdsfig(c.GroupSpace,'classes',c.ClusterSolution.classes,'names',c.names,'sig',sig);
    title('Low ANCOVA group')
    drawnow

    subplot(1,3,3)
    [r,p]=corrcoef(c.dat(c.ancova_codes<median(c.ancova_codes),:)); r(r>.9999) = 1; sig = sign(r) .* (p < .05);
    nmdsfig(c.GroupSpace,'classes',c.ClusterSolution.classes,'names',c.names,'sig',c.STATS.siginteract);
    title('Interactions')

end

% -----------------------------------------------------------------------
% Optional plots of interactions and correlations
% -----------------------------------------------------------------------
if c.doancova && ~isempty(c.ancova_codes)
    [i,j] = find(triu(c.STATS.siginteract));
    if ~isempty(i),
        go = input([num2str(length(i)) ' significant behavior x covariance interactions.  Plot? ']);
        if go
            for ind = 1:length(i)
                cluster_orthviews(cl(i(ind)),{[1 0 0]});
                cluster_orthviews(cl(j(ind)),{[0 0 1]},'add');
                h = axes('Position',[.55 .05 .42 .42]);
                [b,t,p] = ancova(c.ancova_codes,c.dat(:,i(ind)),c.dat(:,j(ind)),1);
                xlabel(cl(i(ind)).shorttitle); ylabel(cl(j(ind)).shorttitle);
                title([]);
                input('Save and press a key for the next one')
            end
        end
    end
end

if c.doancova, c.doancova = 'Yes'; else c.doancova = 'No'; end

if ~any(c.ClusterSolution.pvals < .05)
    go  = input('No significant clustering.  Proceed with analyses anyway?');
    if ~go, return, end
end

% -----------------------------------------------------------------------
% Get average correlations within and between clusters
%
% Apply this clustering solution to the clusters
% Get individual 'network' scores and use them as predictors of behavior
% in multiple regression
% -----------------------------------------------------------------------

c = nmdsfig_tools('apply_clusters',c);

c = nmdsfig_tools('predict_behavior',c,'regions');
c = nmdsfig_tools('predict_behavior',c,'classes');

% -----------------------------------------------------------------------
% brain plots
% -----------------------------------------------------------------------
if isempty(overlay)
    cont = input('Show brain slices? ');
    if ~cont
        diary off

        save nmds_output c cl
        cd ..
        return
    end

    overlay = spm_get(1,'*img','Select overlay image');
else
    % we have overlay, implies we want slices

    cluster_orthviews_classes(cl,c.ClusterSolution.classes,overlay,'axial',0);
    save_figure('cluster_slices_axial');

    %     cluster_orthviews_showcenters(cl,'coronal',overlay);
    %     save_figure('cluster_slices_coronal');
    %
    %     cluster_orthviews_showcenters(cl,'saggital',overlay);
    %     save_figure('cluster_slices_saggital');
end

diary off

save nmds_output c cl
cd ..

return

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

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

% Sub-functions

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

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

% -----------------------------------------------------------------------
% Utility functions
% -----------------------------------------------------------------------
function save_figure(myname)
f1 = gcf;
scn_export_papersetup;
saveas(f1,myname,'png');
return

function f1 = create_wide_figure
scnsize = get(0,'ScreenSize');
f1 = figure('position',[50 50 scnsize(3)-100 scnsize(4)/2],'color','white');
return

function [pthr,sig] = fdr_correct_pvals(p,r)

psq = p; psq(find(eye(size(p,1)))) = 0;
psq = squareform(psq);
pthr = FDR(p,.05);
if isempty(pthr), pthr = 0; end

sig = sign(r) .* (p < pthr);

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%% This function will display the ODE function.

% Last Updated: 2019/05/08
% Coded By: K

function ODEs=Print_ODEs(ODE,n_state,n_control,disp,actual)

z_vars=sym('z',[n_state;1]);
u_vars=sym('u',[n_control;1]);
d_vars=sym('dz',[n_state;1]);

if isa(ODE,'function_handle')
    if n_control~=0
        ODEs=ODE(0,z_vars,u_vars);
    else
        ODEs=ODE(0,z_vars);
    end
end

if disp==1
    if actual==1
        fprintf('\v The actual ODE of the system is/are :\n')
    else
        fprintf('\v The discovered ODE of the system is/are :\n')
    end
     for i=1:n_state
         digits(4)
          fprintf(strcat('\t',char(d_vars(i,1)),'=',char(vpa(ODEs(i,1))),'\n'));
     end
end

author: dynamicslab

repo: SINDy-PI


text:
function [Xi W] = SigmaPoints(xm, P, kappa)
%
%
n  = numel(xm);
Xi = zeros(n, 2*n+1);              % sigma points = col of Xi
W  = zeros(n, 1);

Xi(:, 1) = xm;
W(1)     = kappa / (n + kappa);

U = chol((n+kappa)*P);             % U'*U = (n+kappa)*P

for k=1:n
  Xi(:, k+1) = xm + U(k, :)';      % row of U
  W(k+1)     = 1 / (2*(n+kappa));
end

for k=1:n
  Xi(:, n+k+1) = xm - U(k, :)';
  W(n+k+1)      = 1 / (2*(n+kappa));
end

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
% overlay_conc
%
% callback on mouse down event on spm_image figure
% should display concentration maps on mri
%
% INPUTS:
% in: CSI object
% met: metabolite array
% X: (1,2) vector with the first and last index of x coordinates
% Y: (1,2) vector with the first and last of y coordinates

function overlayConcMapOverMRI
    %call global variable from spm
    global st;
    %default value of 1 for coilNum
    spmFigure = st.fig;

    [metabolitePlot, xIndexBounds, yIndexBounds, voxels] = getCurrentUserDataValues(spmFigure);

    if(isempty(metabolitePlot))
        return
    end
    % cursorPosition in the dimensions (l-r,a-p,s-i)
    cursorPosition = st.centre;
    %dimension size of the MRI
    dimensionSize = st.vols{1}.dim;
    %resolution of MRI (in mm)(not used, might be useful)
    resolution = (st.bb(2,:) - st.bb(1,:) + 1)/dimensionSize(:)';

    [sagital_voxels, coronal_voxels, transverse_voxels] = findIntersectingVoxels(voxels, cursorPosition);

    %3D bounding box of the MRI scan in mm. ie the coordinates where the MRI is
    %plotted onto.
    mriBoundingBox = st.bb;

    deleteCurrentConcentrationMap(st);
    %(needs to be modified if 3D MRSI is to be done)
    voxelIntersectionPositions = {};
    counter = 1;
    metabolitesToPlotVectorized = zeros(1, numel(transverse_voxels));
    for i = 1:numel(transverse_voxels)
        xIndex = transverse_voxels(i).fid_aIndex(1) - xIndexBounds(1) + 1;
        yIndex = transverse_voxels(i).fid_aIndex(2) - yIndexBounds(1) + 1;
        if (xIndex > 0 && yIndex > 0 && xIndex <= size(metabolitePlot, 1) && yIndex <= size(metabolitePlot, 2))
            transverseIntersection = transverse_voxels(i).findIntersection('axial', cursorPosition(3));
            voxelIntersectionPositions{counter} = transverseIntersection(1:2, :) - mriBoundingBox(1,1:2)';
            metabolitesToPlotVectorized(counter) = metabolitePlot(xIndex, yIndex);
            counter = counter + 1;
        end
    end

    for i = 1:length(voxelIntersectionPositions)
        currentCoordiante = voxelIntersectionPositions{i};
        patch(st.vols{1}.ax{4}.ax, currentCoordiante(1, :), ...
            currentCoordiante(2, :), metabolitesToPlotVectorized(i), ...
            'Tag', 'trans_plot');
    end
    %PLOTTING IN THE SAGITAL PLANE
    plotVoxelsOnMRIPlane(st.vols{1}.ax{3}.ax, sagital_voxels, 'sagital', cursorPosition);

    %PLOTTING ON THE CORONAL PLANE
    plotVoxelsOnMRIPlane(st.vols{1}.ax{2}.ax, coronal_voxels, 'coronal', cursorPosition);
end

function [metabolitePlot, xIndexBounds, yIndexBounds, voxels] = getCurrentUserDataValues(spmFigure)
    figureUserData = spmFigure.UserData;
    metabolitePlot = figureUserData.metabolitePlot;
    xIndexBounds = figureUserData.xIndexBounds;
    yIndexBounds = figureUserData.yIndexBounds;
    voxels = figureUserData.voxels;
end

function deleteCurrentConcentrationMap(st)
    concentrationMap = findobj(st.vols{1}.ax{4}.ax,'Tag','trans_plot');
    delete(concentrationMap);
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function simulate_2spk_mix(data_type, wsj0root, output_dir, fs8k, min_max)
% Simulate 2-speaker mixture data for speaker extraction.
% Call:
%     simulate_2spk_mix(data_type, wsj0root, output_dir, fs8k, min_max)
%     e.g., simulate_2spk_mix('tt', '/media/clx214/data/wsj/', '/media/clx214/data/wsj0_2mix_extr_tmp/wav8k', 8000, 'max')
% Paras:
%     data_type: data set to generate, (tr|cv|tt), e.g., 'tt'
%     wsj0root: YOUR_PATH/, the folder containing converted wsj0/, e.g., '/media/clx214/data/wsj/'
%     output_dir: the folder to save simulated data for extraction, e.g., '/media/clx214/data/wsj0_2mix_extr_tmp/wav8k'
%     fs8k: sampling rate of the simulated data, e.g., 8000
%     min_max: get the mininium or maximum wav length, when simulating mixture data, e.g, 'max'
%
% The code is based on "create_wav_2speakers_extr.m" from "http://www.merl.com/demos/deep-clustering"
%
% 1. Assume that WSJ0's wv1 sphere files is converted to wav files. The folder
%    structure and file name are kept same under wsj0/, e.g.,
%    ORG_PATH/wsj0/si_tr_s/01t/01to030v.wv1 is converted to wav and
%    stored in YOUR_PATH/wsj0/si_tr_s/01t/01to030v.wv1.
%    Relevant data ('si_tr_s', 'si_dt_05' and 'si_et_05') are under YOUR_PATH/wsj0/
% 2. Put 'voicebox' toolbox in current folder. (http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html)
% 3. Set your 'YOUR_PATH' and 'OUTPUT_PATH' properly, then run this script in Matlab.
%    (The max lenght of the wav will be kept when generate the mixture. The sampling rate will be 8kHz.)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Copyright (C) 2016 Mitsubishi Electric Research Labs
%                          (Jonathan Le Roux, John R. Hershey, Zhuo Chen)
%   Apache 2.0  (http://www.apache.org/licenses/LICENSE-2.0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Copyright 2018 Chenglin Xu, Nanyang Technological University, Singapore
%
%Licensed under the Apache License, Version 2.0 (the "License");
%you may not use this file except in compliance with the License.
%You may obtain a copy of the License at
%
%    http://www.apache.org/licenses/LICENSE-2.0
%
%Unless required by applicable law or agreed to in writing, software
%distributed under the License is distributed on an "AS IS" BASIS,
%WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
%See the License for the specific language governing permissions and
%limitations under the License.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist([output_dir '/' min_max '/' data_type],'dir')
    mkdir([output_dir '/' min_max '/' data_type]);
end
mkdir([output_dir  '/' min_max '/' data_type '/s1/']);
mkdir([output_dir  '/' min_max '/' data_type '/aux/']);
mkdir([output_dir  '/' min_max '/' data_type '/mix/']);

TaskFile = ['mix_2_spk_' data_type '_extr.txt'];
fid = fopen(TaskFile,'r');
C = textscan(fid,'%s %f %s %f %s');
num_files = length(C{1});

fprintf(1,'Start to generate data for %s\n', [min_max '_' data_type]);
for i = 1:num_files
    [inwav1_dir,invwav1_name,inwav1_ext] = fileparts(C{1}{i});
    [inwav2_dir,invwav2_name,inwav2_ext] = fileparts(C{3}{i});
    [inwav_aux_dir,invwav_aux_name,inwav_aux_ext] = fileparts(C{5}{i});

    inwav1_snr = C{2}(i);
    inwav2_snr = C{4}(i);
    mix_name = [invwav1_name,'_',num2str(inwav1_snr),'_',invwav2_name,'_',num2str(inwav2_snr),'_',invwav_aux_name];

    % get input wavs
    [s1, fs] = audioread([wsj0root C{1}{i}]);
    s2       = audioread([wsj0root C{3}{i}]);
    s_aux    = audioread([wsj0root C{5}{i}]);

    % resample, normalize to 8 kHz file
    s1_8k = resample(s1,fs8k,fs);
    [s1_8k,lev1] = activlev(s1_8k,fs8k,'n'); % y_norm = y /sqrt(lev);
    s2_8k = resample(s2,fs8k,fs);
    [s2_8k,lev2] = activlev(s2_8k,fs8k,'n');
    s_aux_8k = resample(s_aux,fs8k,fs);
    [s_aux_8k,lev_aux] = activlev(s_aux_8k,fs8k,'n');

    weight_1 = 10^(inwav1_snr/20);
    weight_2 = 10^(inwav2_snr/20);

    s1_8k = weight_1 * s1_8k;
    s2_8k = weight_2 * s2_8k;

    switch min_max
        case 'max'
            mix_8k_length = max(length(s1_8k),length(s2_8k));
            s1_8k = cat(1,s1_8k,zeros(mix_8k_length - length(s1_8k),1));
            s2_8k = cat(1,s2_8k,zeros(mix_8k_length - length(s2_8k),1));
        case 'min'
            mix_8k_length = min(length(s1_8k),length(s2_8k));
            s1_8k = s1_8k(1:mix_8k_length);
            s2_8k = s2_8k(1:mix_8k_length);
    end
    mix_8k = s1_8k + s2_8k;

    max_amp_8k = max(cat(1,abs(mix_8k(:)),abs(s1_8k(:)),abs(s2_8k(:)),abs(s_aux_8k(:))));
    mix_scaling_8k = 1/max_amp_8k*0.9;
    s1_8k = mix_scaling_8k * s1_8k;
    mix_8k = mix_scaling_8k * mix_8k;
    s_aux_8k = mix_scaling_8k * s_aux_8k;

    audiowrite([output_dir '/' min_max '/' data_type '/s1/' mix_name '.wav'],s1_8k,fs8k);
    audiowrite([output_dir '/' min_max '/' data_type '/aux/' mix_name '.wav'],s_aux_8k,fs8k);
    audiowrite([output_dir '/' min_max '/' data_type '/mix/' mix_name '.wav'],mix_8k,fs8k);
end
fclose(fid);
fprintf(1,'End of generating data for %s\n', [min_max '_' data_type]);
end

author: gemengtju


text:
function plotSteeringComparison(files, laps, labels)
% Author:       Alexander Wischnewski
% Description:
%   function used to compare different controllers
% Inputs/parameters:
%   files:      Cell array with file names of logs
%   laps:       laps to evaluate (matrix with start and end lap for each datafile)
%   labels:     Cell array with abels for the datasets used for legends
%               (file names are used if not given)

% plotting parameters
LineWidth = 1;
DyLim = 1.0;
axLim = 25;
ayLim = 25;

% load all the relevant data files
for i = 1:1:length(files)
    data{i} = load(files{i});
end

% check if labels are given, if not use file names
if(nargin <= 2)
    labels = files;
end

% find start and end indices
for i = 1:1:length(files)
    idx_start{i} = find((data{i}.debug.debug_mvdc_path_matching_debug_ActualTrajPoint_LapCnt.Data == laps(i,1)), 1, 'first') + 50;
    idx_end{i} = find((data{i}.debug.debug_mvdc_path_matching_debug_ActualTrajPoint_LapCnt.Data == laps(i,2)), 1, 'last') - 50;
end

figure;
hold on; grid on;
colors = {'r', 'b'};
for i = 1:1:length(files)
    plot(data{i}.debug.debug_mvdc_path_matching_debug_ActualTrajPoint_s_glob_m.Data(idx_start{i}:idx_end{i}), ...
        data{i}.debug.debug_mvdc_actuator_debug_RequestSteeringAngle_rad.Data(idx_start{i}:idx_end{i}), ...
        'LineWidth', LineWidth, 'Color', colors{i}, 'DisplayName', [labels{i} ' Target']);
    plot(data{i}.debug.debug_mvdc_path_matching_debug_ActualTrajPoint_s_glob_m.Data(idx_start{i}:idx_end{i}), ...
        data{i}.debug.debug_VehicleSensorData_Delta_Wheel_rad.Data(idx_start{i}:idx_end{i}), '--',...
        'LineWidth', LineWidth, 'Color', colors{i}, 'DisplayName', [labels{i} ' Actual']);
end
xlabel('s in m'); ylabel({'Steering', 'in rad'}); ylim([-0.02, 0.02]);
legend();

author: TUMFTM


text:
function a = band(a,p,q)
%BAND         Extract band from matrix a, lower bandwidth p, upper bandwidth q
%   if parameter q is omitted, q:=p
%
%   res = band(a,p,q)
%

% written  12/02/96     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    improved performance
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if nargin<3
    q = p;
  end

  a = tril(triu(a,-p),q);

repo: OpenMAS


text:
function [reg_c,rho_c,eta_c] = l_corner(rho,eta,reg_param,U,s,b,method,M)
%L_CORNER Locate the "corner" of the L-curve.
%
% [reg_c,rho_c,eta_c] =
%        l_corner(rho,eta,reg_param)
%        l_corner(rho,eta,reg_param,U,s,b,method,M)
%        l_corner(rho,eta,reg_param,U,sm,b,method,M) ,  sm = [sigma,mu]
%
% Locates the "corner" of the L-curve in log-log scale.
%
% It is assumed that corresponding values of || A x - b ||, || L x ||,
% and the regularization parameter are stored in the arrays rho, eta,
% and reg_param, respectively (such as the output from routine l_curve).
%
% If nargin = 3, then no particular method is assumed, and if
% nargin = 2 then it is issumed that reg_param = 1:length(rho).
%
% If nargin >= 6, then the following methods are allowed:
%    method = 'Tikh'  : Tikhonov regularization
%    method = 'tsvd'  : truncated SVD or GSVD
%    method = 'dsvd'  : damped SVD or GSVD
%    method = 'mtsvd' : modified TSVD,
% and if no method is specified, 'Tikh' is default.  If the Spline Toolbox
% is not available, then only 'Tikh' and 'dsvd' can be used.
%
% An eighth argument M specifies an upper bound for eta, below which
% the corner should be found.

% Per Christian Hansen, IMM, July 26, 2007.
%
%
% This file is part of regtools [1] and covered by the BSD License
%
% Copyright (c) 2008, Per Christian Hansen
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
%
% References:
% [1] http://www.mathworks.com/matlabcentral/fileexchange/file_infos/52-regtools

% Set default regularization method.
if (nargin <= 3)
  method = 'none';
  if (nargin==2), reg_param = (1:length(rho))'; end
else
  if (nargin==6), method = 'Tikh'; end
end

% Set this logical variable to 1 (true) if the corner algorithm
% should always be used, even if the Spline Toolbox is available.
alwayscorner = 0;

% Set threshold for skipping very small singular values in the
% analysis of a discrete L-curve.
s_thr = eps;  % Neglect singular values less than s_thr.

% Set default parameters for treatment of discrete L-curve.
deg   = 2;  % Degree of local smooting polynomial.
q     = 2;  % Half-width of local smoothing interval.
order = 4;  % Order of fitting 2-D spline curve.

% Initialization.
if (length(rho) < order)
  error('Too few data points for L-curve analysis')
end
if (nargin > 3)
  [p,ps] = size(s); [m,n] = size(U);
  beta = U'*b;
  if (m>n), b0 = b - U*beta; end
  if (ps==2)
    s = s(p:-1:1,1)./s(p:-1:1,2);
    beta = beta(p:-1:1);
  end
  xi = beta./s;
end

% Restrict the analysis of the L-curve according to M (if specified).
if (nargin==8)
  index = find(eta < M);
  rho = rho(index); eta = eta(index); reg_param = reg_param(index);
end

if (strncmp(method,'Tikh',4) | strncmp(method,'tikh',4))

  % The L-curve is differentiable; computation of curvature in
  % log-log scale is easy.

  % Compute g = - curvature of L-curve.
  g = lcfun(reg_param,s,beta,xi);

  % Locate the corner.  If the curvature is negative everywhere,
  % then define the leftmost point of the L-curve as the corner.
  [gmin,gi] = min(g);
  reg_c = fminbnd('lcfun',...
    reg_param(min(gi+1,length(g))),reg_param(max(gi-1,1)),...
    optimset('Display','off'),s,beta,xi); % Minimizer.
  kappa_max = - lcfun(reg_c,s,beta,xi); % Maximum curvature.

  if (kappa_max < 0)
    lr = length(rho);
    reg_c = reg_param(lr); rho_c = rho(lr); eta_c = eta(lr);
  else
    f = (s.^2)./(s.^2 + reg_c^2);
    eta_c = norm(f.*xi);
    rho_c = norm((1-f).*beta);
    if (m>n), rho_c = sqrt(rho_c^2 + norm(b0)^2); end
  end

elseif (strncmp(method,'tsvd',4) | strncmp(method,'tgsv',4) | ...
        strncmp(method,'mtsv',4) | strncmp(method,'none',4))

  % Use the adaptive pruning algorithm to find the corner, if the
  % Spline Toolbox is not available.
  if ~exist('splines','dir') | alwayscorner
    %error('The Spline Toolbox in not available so l_corner cannot be used')
    reg_c = corner(rho,eta);
    rho_c = rho(reg_c);
    eta_c = eta(reg_c);
    return
  end

  % Othersise use local smoothing followed by fitting a 2-D spline curve
  % to the smoothed discrete L-curve. Restrict the analysis of the L-curve
  % according to s_thr.
  if (nargin > 3)
    if (nargin==8)       % In case the bound M is in action.
      s = s(index,:);
    end
    index = find(s > s_thr);
    rho = rho(index); eta = eta(index); reg_param = reg_param(index);
  end

  % Convert to logarithms.
  lr = length(rho);
  lrho = log(rho); leta = log(eta); slrho = lrho; sleta = leta;

  % For all interior points k = q+1:length(rho)-q-1 on the discrete
  % L-curve, perform local smoothing with a polynomial of degree deg
  % to the points k-q:k+q.
  v = (-q:q)'; A = zeros(2*q+1,deg+1); A(:,1) = ones(length(v),1);
  for j = 2:deg+1, A(:,j) = A(:,j-1).*v; end
  for k = q+1:lr-q-1
    cr = A\lrho(k+v); slrho(k) = cr(1);
    ce = A\leta(k+v); sleta(k) = ce(1);
  end

  % Fit a 2-D spline curve to the smoothed discrete L-curve.
  sp = spmak((1:lr+order),[slrho';sleta']);
  pp = ppbrk(sp2pp(sp),[4,lr+1]);

  % Extract abscissa and ordinate splines and differentiate them.
  % Compute as many function values as default in spleval.
  P     = spleval(pp);  dpp   = fnder(pp);
  D     = spleval(dpp); ddpp  = fnder(pp,2);
  DD    = spleval(ddpp);
  ppx   = P(1,:);       ppy   = P(2,:);
  dppx  = D(1,:);       dppy  = D(2,:);
  ddppx = DD(1,:);      ddppy = DD(2,:);

  % Compute the corner of the discretized .spline curve via max. curvature.
  % No need to refine this corner, since the final regularization
  % parameter is discrete anyway.
  % Define curvature = 0 where both dppx and dppy are zero.
  k1    = dppx.*ddppy - ddppx.*dppy;
  k2    = (dppx.^2 + dppy.^2).^(1.5);
  I_nz  = find(k2 ~= 0);
  kappa = zeros(1,length(dppx));
  kappa(I_nz) = -k1(I_nz)./k2(I_nz);
  [kmax,ikmax] = max(kappa);
  x_corner = ppx(ikmax); y_corner = ppy(ikmax);

  % Locate the point on the discrete L-curve which is closest to the
  % corner of the spline curve.  Prefer a point below and to the
  % left of the corner.  If the curvature is negative everywhere,
  % then define the leftmost point of the L-curve as the corner.
  if (kmax < 0)
    reg_c = reg_param(lr); rho_c = rho(lr); eta_c = eta(lr);
  else
    index = find(lrho < x_corner & leta < y_corner);
    if ~isempty(index)
      [dummy,rpi] = min((lrho(index)-x_corner).^2 + (leta(index)-y_corner).^2);
      rpi = index(rpi);
    else
      [dummy,rpi] = min((lrho-x_corner).^2 + (leta-y_corner).^2);
    end
    reg_c = reg_param(rpi); rho_c = rho(rpi); eta_c = eta(rpi);
  end

elseif (strncmp(method,'dsvd',4) | strncmp(method,'dgsv',4))

  % The L-curve is differentiable; computation of curvature in
  % log-log scale is easy.

  % Compute g = - curvature of L-curve.
  g = lcfun(reg_param,s,beta,xi,1);

  % Locate the corner.  If the curvature is negative everywhere,
  % then define the leftmost point of the L-curve as the corner.
  [gmin,gi] = min(g);
  reg_c = fminbnd('lcfun',...
    reg_param(min(gi+1,length(g))),reg_param(max(gi-1,1)),...
    optimset('Display','off'),s,beta,xi,1); % Minimizer.
  kappa_max = - lcfun(reg_c,s,beta,xi,1); % Maximum curvature.

  if (kappa_max < 0)
    lr = length(rho);
    reg_c = reg_param(lr); rho_c = rho(lr); eta_c = eta(lr);
  else
    f = s./(s + reg_c);
    eta_c = norm(f.*xi);
    rho_c = norm((1-f).*beta);
    if (m>n), rho_c = sqrt(rho_c^2 + norm(b0)^2); end
  end

else
  error('Illegal method')
end

repo: REST


text:
function [res, diams] = imOrientedGranulo(img, angleList, granuloType, strelSizes, varargin)
% Gray level granulometry mean size for various orientations.
%
%   Compute typical size of bright or dark structures within images by
%   performing gray-level granulometries with horizontal structuring
%   elements and for various rotations of the input image.
%   The function "imDirectionalGranulo" is an alternative approach based on
%   line structuring elements with various orientations.
%
%   Usage
%   RES = imOrientedGranulo(IMG, ANGLES, TYPE, SIZES)
%   IMG should be a 2D image (binary, grayscale or color)
%   ANGLES is the list of angles to consider, in degrees, as a 1-by-N row
%     vector
%   GRTYPE can be one of {'opening', 'closing', 'erosion', 'dilation'}.
%   SIZES are given as radius in pixels. Diameters of strels are obtained
%     as 2*R+1.
%   The result GR is a 1-by-N array with as many columns as the number of
%     elements provided in ANGLES array.
%
%   [RES, DIAMS] = imOrientedGranulo(...)
%   Also returns the diameters used for each morphological filtering step.
%
%
%   Example
%   imOrientedGranulo
%
%   See also
%     imGranulometry, imDirectionalGranulo, imGranulo, imGranuloByRegion,
%     granuloMeanSize
%
%   Reference
%   The methodology is described in the following article:
%   "Exploring the microstructure of natural fibre composites by confocal
%   Raman imaging and image analysis", by Antoine Gallos, Gabriel Paës,
%   David Legland, Florent Allais, Johnny Beaugrand (2017).
%   Composites Part A: Applied Science and Manufacturing 94, p. 32-40.
%   doi: https://doi.org/10.1016/j.compositesa.2016.12.005
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2015-12-02,    using Matlab 8.6.0.267246 (R2015b)
% Copyright 2015 INRAE - Cepia Software Platform.

% how to interpolate images
if islogical(img)
    interp = 'nearest';
else
    interp = 'linear';
end

% convert sizes from radius to diameter
diams = 2 * strelSizes + 1;

% allocate result array
nAngles = length(angleList);
res = zeros(1, nAngles);

% iterate over the orientations
for iAngle = 1:nAngles
    angle = angleList(iAngle);
    imgr = imrotate(img, angle, interp);
    grCurve = imGranulo(imgr, granuloType, 'lineh', strelSizes);
    res(iAngle) = granuloMeanSize(grCurve, diams);
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [ radiance ] = dehaze_fast( image, omega, win_size )
%DEHZE Summary of this function goes here
%   Detailed explanation goes here

if ~exist('omega', 'var')
    omega = 0.95;
end

if ~exist('win_size', 'var')
    win_size = 15;
end

r = 15;
res = 0.001;

[m, n, ~] = size(image);

dark_channel = get_dark_channel(image, win_size);

atmosphere = get_atmosphere(image, dark_channel);

trans_est = get_transmission_estimate(image, atmosphere, omega, win_size);

x = guided_filter(rgb2gray(image), trans_est, r, res);

transmission = reshape(x, m, n);

radiance = get_radiance(image, transmission, atmosphere);

end

author: sjtrny

repo: Dark-Channel-Haze-Removal


text:
function [SNR] = bvpsnr(BVP, FS, HR, PlotTF)
%BVPSNR Estimates the signal-to-noise ratio of the blood volume pulse
% signal. Adapted from the method by G. de Haan, TBME, 2013.
% SNR calculated as the ratio (in dB) of power contained within +/- 0.1 Hz
% of the reference heart rate frequency and +/- 0.2 of its first
% harmonic and sum of all other power between 0.5 and 4 Hz.
% Adapted from the method by G. de Haan, TBME, 2013
%
%   Inputs:
%       BVP                     = A BVP timeseries.
%       FS                      = The sample rate of the BVP time series (Hz/fps).
%       HR                      = The reference heart rate (Hz/fps).
%       PlotTF                  = Boolean to turn plotting results on or off.
%
%   Outputs:
%       SNR                     = Blood Volume Pulse Signal-to-Noise Ratio.
%
% Daniel McDuff, Ethan Blackford, January 2019
% Copyright (c)
% Licensed under the MIT License and the RAIL AI License.

HR_F=HR/60;

NyquistF = FS/2;
FResBPM = 0.5; %resolution (bpm) of bins in power spectrum used to determine PR and SNR
N = (60*2*NyquistF)/FResBPM; %number of bins in power spectrum

%% Construct Periodogram
[Pxx,F] = periodogram(BVP,hamming(length(BVP)),N,FS);
GTMask1 = (F >= HR_F-0.1)&(F <= HR_F+0.1);
GTMask2 = (F >= HR_F*2-0.2)&(F <= HR_F*2+0.2);
SPower = sum(Pxx(GTMask1|GTMask2));
FMask2 = (F >= 0.5)&(F <= 4);
AllPower = sum(Pxx(FMask2));
SNR = pow2db(SPower/(AllPower-SPower));

%% Optionally plot the power spectrum and regions used to calculate the SNR
if(PlotTF)
    % Plot spower spectrum and SNR regions
    figure
    plot(F,pow2db(Pxx))
    title('Power Spectrum and SNR Regions')
    xlabel('Frequency (Hz)')
    ylabel('Power (dB)')
    xlim([0.5 4])
    ylimreg=ylim;
    hold on

    % HR peak
    %line([HR_F,HR_F],ylim,'Color','magenta','LineStyle','-.')

    % HR peak region
    line([HR_F-0.1,HR_F-0.1],ylim,'Color','red','LineStyle','--')
    line([HR_F+0.1,HR_F+0.1],ylim,'Color','red','LineStyle','--')

    % First harmonic
    line([HR_F*2-0.2,HR_F*2-0.2],ylim,'Color','red','LineStyle','--')
    line([HR_F*2+0.2,HR_F*2+0.2],ylim,'Color','red','LineStyle','--')

    % Overall power region
    line([0.5,0.5],ylim,'Color','black','LineStyle','-')
    line([4,4],ylim,'Color','black','LineStyle','-')

    % Adjust axes
    xlim([0 4.5])
    ylim(ylimreg);
end

end

author: danmcduff

repo: iphys-toolbox


text:
% TD | HoRPCA-S | HoRPCA with Singleton model solved by ADAL (Goldfarb and Qin, 2013)
% process_video('TD', 'HoRPCA-S', 'dataset/demo.avi', 'output/demo_HoRPCA-S.avi');

alg_path_aux = fullfile(lrs_conf.td_path,'RLRT');
addpath(genpath(alg_path_aux));

pdata.T = T;
pdata.X = T;
N = ndims(pdata.T);
r = 1/sqrt(max(size(pdata.T)));
params.E0 = tenzeros(size(pdata.T));
params.X0 = tenzeros(size(pdata.T));
params.V0 = cell(1, N);
for i = 1:N
  params.V0{i} = tenzeros(size(pdata.T));
end
params.mu0 = 1/(N+1);
params.mode = N;
params.IsTC = false; % is tensor completion
params.rRatio = 1/4;
params.opt_tol = 1e-3;
params.eta = 1/(N+1);
params.max_iter = 1000;
params.mu1fac = 10;
params.mu1 = params.mu1fac*std(T(:));
params.mu2 = params.mu1;
params.mu_min = 1e-4;
params.mu_max = 1e2;
params.lambdaS = 1;
params.lambda = params.lambdaS*r*params.rRatio;
params.verbose = 1;
params.use_cont = true;
params.k = [size(T,1) size(T,2) 1];
%%%%%%%%%% for PROPACK %%%%%%%%%%%%
% declare global var 'sv'
global sv;
global tmode;
global use_propack;
global curr_mu;
sv =  ceil(min(size(pdata.T)) * 0.1) * ones( 1, N );
use_propack = true;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
results = tensor_rpca_adal2(pdata, params); % tensor_rpca_adal
L = double(results.X);
S = double(results.E);
clear sv tmode use_propack curr_mu;

rmpath(genpath(alg_path_aux));

author: andrewssobral

repo: lrslibrary


text:
function T = myreshape(T, sizes)
% MYRESHAPE Like the built-in reshape, except myreshape(T,n) == reshape(T,[n 1])
% T = myreshape(T, sizes)

if length(sizes)==0
  return;
elseif length(sizes)==1
  T = reshape(T, [sizes 1]);
else
  T = reshape(T, sizes(:)');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function feaMap = rescaleFeaMap(feaMap)
feaMap = feaMap - min(feaMap(:));
feaMap = feaMap ./ max(feaMap(:));

author: aimerykong


text:
% NTF | NTD-HALS | Non-negative Tucker Decomposition solved by Hierarchical ALS  (Zhou et al. 2012)
% process_video('NTF', 'NTD-HALS', 'dataset/demo.avi', 'output/demo_NTD-HALS.avi');

alg_path_aux = fullfile(lrs_conf.ntf_path,'lraNTD');
addpath(genpath(alg_path_aux));

R = [size(T,1) size(T,2) 2];
opts = struct('NumOfComp',R,'nlssolver','hals','maxiter',100,...
  'maxiniter',20,'tdalgFile','call_tucker_als_opts.mat');

T_hat = lraNTD_ANLS(T, opts);

L = double(tensor(T_hat));
S = double(T) - L;

rmpath(genpath(alg_path_aux));


text:
classdef MomentumParameterFactory < handle

    methods (Access = public, Static)

        function obj = create(cParams)
            switch cParams.type
                case 'CONSTANT'
                    obj = MomentumParameterConstant(cParams);
                case 'CONVEX'
                    obj = MomentumParameterConvexCase(cParams);
                case 'STRICTLY CONVEX'

            end

        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%Help file for INTLAB Version 2
%
%Changes:
%
%  Input/output:
%
%  - INTLAB output is rigorous [ Note that INTLAB V1 input is already rigorous
%      when using character constants (e.g. intval('.1')) ]
%  - rigorous standard functions
%  - Interval display by uncertainty, e.g. 3.14159_
%  - Interval input by string, also with tolerances (e.g. '3.14_',
%      '[3,4]' or '<-3.1e2,0.01_>' or '<3-4i,.1> etc.
%  - Permanent switch of display mode of intervals thru intvalinit
%  - Function initvar for gradients replaced by gradientinit (replace
%      by global change, exactly the same functionality); was necessary
%      because of ambiguouity with initialization of slope variables
%
%
%  Standard functions:
%
%  - Rigorous standard functions incorporated
%
%
%  Slopes:
%
%  - Slope toolbox added
%
%
%  Long precision:
%
%  - Long toolbox added
%
%
%  Others:
%
%  - Simplified call of verifynlss for one-dimensional nonlinear functions
%  - Extended arithmetic including +/-infinity, e.g. division by
%      zero intervals etc.
%  - Some changes in subdirectory structure
%
%
%
%Further changes:
%
%  - Linear system solver in 2 stages with sharp interval hull for
%      preconditioned system
%  - Components of arrays of intervals, gradients... may be erased by
%      assignment of an empty set, e.g. A=intval(ones(5)); A(2,:)=[];
%  - Functions SetRoundDown, SetRoundUp and SetRoundNear replaced by
%      one routine setround(rnd)
%  - Functions "in", "in0", "isnan" deliver array of 0/1 (like Matlab)
%  - Functions pred/succ with second parameter for k-th predecessor/successor
%  - Function power10tobinary with increased exponent range, changes also
%      power10tobinary.mat
%  - Fix of one or the other bug, thanks to many users
%
%New functions despite the above:
%
%  - mig          Mignitude for points and intervals
%  - compmat      Ostrowski's comparison matrix
%  - band         Extract band out of matrix
%  - bandwidth    Get lower and upper bandwidth
%  - midradcmplx  Always produces complex interval
%  - binom        Binomial coefficient
%  - distbits     Distance in bits of two reals
%  - getround     Get current rounding mode
%  - gregk...     Some sparse test matrices from Gregory/Karney test set
%  - perms_       Fast version of Matlab perms
%  - randint      Random integers in specified range
%  - random       Random numbers uniformly distrubuted in [-1,1] or
%                   within specified range
%  - randomc      Complex random numbers with real and imaginary part as random
%  - relerr       Relative error for scalars, vectors and matrices
%


text:
classdef VonMisesD
%%VONMISESD Functions to handle the von Mises distribution, which is a
%    common circular distribution. The distribution is given in Chapter
%    3.5.4 of [1].
%Implemented methods are: circVar, PDF, CDF,trigMoment, normProdDist,
%                         convDistApprox, params4TrigMoment, rand, entropy
%
%REFERENCES:
%[1] K. V. Mardia and P. E. Jupp, Directional Statistics. Chichester: John
%    Wiley and Sons, 2000.
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

methods(Static)

function val=circVar(kappa)
%%CIRCVAR Obtain the circular variance of a particular von Mises
%         distribution.
%
%INPUTS: kappa The concentration of the distribution 0<=kappa<inf.
%
%OUTPUTS: val The circular variance of the distribution.
%
%The circular variance is defined as 1-abs(r), where r is the mean
%resultant value of the second moment. The mean resultant values for all
%moment of the von Mises distribution are given in Chapter 2.2.4 of [1].
%
%EXAMPLE:
%Here, we verify that the result agrees with the value obtained by
%sampling.
% numRuns=1e5;
% kappa=100;
% mu=0.25;
% x=VonMisesD.rand([numRuns,1],mu,kappa);
% circVarEst=findCircVar(x);
% abs(VonMisesD.circVar(kappa)-circVarEst)./abs(VonMisesD.circVar(kappa))
%One will see that the relative error is on the order of 1e-3.
%
%REFERENCES:
%[1] S. R. Jammalamadaka and A. SenGupta, Topics in Circular Statistics.
%    Singapore: World Scientific, 2001.
%
%April 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

val=1-BesseliRatio(1,kappa);

end

function val=PDF(x,mu,kappa)
%%PDF Evaluate a von Mises PDF at a given point with the specified
%     parameters.
%
%INPUTS: x The point(s) at which the von Mises distribution is to be
%          evaluated.
%       mu The mean direction of the distribution. This is normally between
%          -pi and pi.
%    kappa The concentration of the distribution 0<=kappa<inf.
%
%OUTPUTS: val The value(s) of the von Mises PDF with the specified
%             parameters evaluated at x.
%
%The von Mises distribution is used to model errors on angular quantities.
%It is described in [1] and [2].
%
%REFERENCES:
%[1] N. Fisher, Statistical Analysis of Circular Data. Cambridge, UK:
%    Cambridge University Press, 1993.
%[2] K. V. Mardia and P. E. Jupp, Directional Statistics. Chichester,
%    England: John Wiley & Sons, 2000.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=exp(kappa*cos(x-mu))/(2*pi*besseli(0,kappa));
end

function vals=CDF(x,mu,kappa,startAng)
%%CDF Evaluate  the cumulative distribution function (CDF) of the von Mises
%     distribution. Since this is a circular distribution, the CDF has to
%     be defined in terms of a specific starting angle. Increading density
%     goes in the direction of increasing angle from the starting point
%     (around the circle).
%
%INPUTS: z The point(s) at which the CDF should be evaluated.
%       mu The mean direction of the distribution. This is normally between
%          -pi and pi.
%    kappa The concentration of the distribution 0<=kappa<inf.
%  startAng The starting angle. This is where the CDF is zero. The default
%           if this parameter is omitted or an empty matrix is passed is
%           zero.
%
%OUTPUTS: vals The scalar value(s) of the von Mises CDF with the given PDF
%              evaluated at the points in z.
%
%The function Besseli0Inc gives the CDF of the von Mises distribution
%starting from -pi. Thus, by taking the difference of two evaluations of
%Besseli0Inc (and determining whether we are going the long or the short
%way around the circle) one can evaluate the CDF.
%
%August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

if(nargin<4||isempty(startAng))
    startAng=-pi;
end

vals=Besseli0Inc(x-mu,kappa)-Besseli0Inc(startAng-mu,kappa);

sel=vals<0;
vals(sel)=vals(sel)+1;
end

function [rho,theta,R]=trigMoment(n,mu,kappa)
%%TRIGMOMENT Compute the nth trigonometric moment of a von Mises
%            distribution having the given parameters. A direction theta on
%            the unit circle can be modeled as a complex quantity having
%            unit magnitude as exp(1j*theta). The nth raw complex
%            trigonometric moment is the expected value of exp(1j*n*theta),
%            where the integral for the expected value is taken over any
%            interval of length 2*pi.
%
%INPUTS: n The order of the moment desired. This is >=1.
%       mu The mean direction of the distribution. This is normally between
%          -pi and pi.
%    kappa The concentration of the distribution 0<=kappa<inf. For
%          kappa>500, the results might be inaccurate.
%
%OUTPUTS: rho The (complex) mean resultant value for the nth moment. Note
%             that rho=R*exp(1j*theta).
%       theta The (real) trigonometric mean angle in radians for the nth
%             moment. This is between -pi and pi.
%           R The (real) mean resultant length for the nth moment.
%
%The moments of a von Mises distribution are given in Chapter 2.2.4 of [1]
%and Chapter 3.5.4 of [2].
%
%EXAMPLE:
%Here, we verify the first trigonometric moment based on samples.
% numRuns=1e5;
% kappa=100;
% mu=0.25;
% x=VonMisesD.rand([numRuns,1],mu,kappa);
% [rhoS,thetaS,RS]=findTrigMomentFromSamp(1,x);
% [rho,theta,R]=VonMisesD.trigMoment(1,mu,kappa);
% abs(rhoS-rho)
% abs(thetaS-theta)
% abs(RS-R)
%Here, we will see that the errors tend to be on the order of 1e-4,
%dominated by the angular error.
%
%REFERENCES:
%[1] S. R. Jammalamadaka and A. SenGupta, Topics in Circular Statistics.
%    Singapore: World Scientific, 2001.
%[2] K. V. Mardia and P. E. Jupp, Directional Statistics. Chichester: John
%    Wiley and Sons, 2000.
%
%August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

%For n=1 or 2, there should be no need to use anything but the default
%maximum number of iterations in BesseliRatio.
if(n==1)
    rho=BesseliRatio(1,kappa)*exp(1j*n*mu);
elseif(n==2)
    rho=BesseliRatio(2,kappa)*BesseliRatio(1,kappa)*exp(1j*n*mu);
else
    rho=besseli(n,kappa,1)/besseli(0,kappa,1)*exp(1j*n*mu);
end

if(nargout>1)
    theta=angle(rho);
    R=abs(rho);
end
end

function [mu,kappa]=normProdDist(mu1,kappa1,mu2,kappa2)
%%NORMPRODDIST The product of two von Mises distributions is an
%              unnormalized von Mises distribution. This finds the
%              parameters of the product distribution.
%
%INPUTS: mu1, kappa1 The mean direction and concentration of the first
%                    distribution. Note that 0<=kappa<inf.
%        mu2, kappa2 The mean direction and concentration of the second
%                    distribution.
%
%OUTPUTS: mu,kappa The mean direction and concentration of the product
%                  distribution -pi<=mu<pi.
%
%The product of two von Mises distributions is described in [1].
%
%REFERENCES:
%[1] M. Azmani, S. Reboul, J.-B. Choquel, and M. Benkelloun, "A recursive
%    fusion filter for angular data," in Proceedings of the 2009 IEEE
%    International Conference on Robotics and Biomimetics, Guilin, China,
%    19-23 Dec. 2009, pp. 882-887.
%
%August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

C=kappa1*cos(mu1)+kappa2*cos(mu2);
S=kappa1*sin(mu1)+kappa2*sin(mu2);
mu=atan2(S,C);
kappa=sqrt(C^2+S^2);
end

function [mu,kappa]=convDistApprox(mu1,kappa1,mu2,kappa2)
%%CONVDISTAPPROX The convolution of two von Mises distributions is not a
%          von Mises distribution. However, it can be approximated by a von
%          Mises distribution that matches the moments, as is done in
%          Equation 7 of [1]. This function returns the paramters of the
%          approximation.
%
%INPUTS: mu1, kappa1 The mean direction and concentration of the first
%                    distribution. Note that 0<=kappa<inf.
%        mu2, kappa2 The mean direction and concentration of the second
%                    distribution.
%
%OUTPUTS: mu, kappa The mean direction and concentration of a von Mises
%                   approximation to the product distribution. -pi<=mu<pi.
%
%REFERENCES:
%[1] I. Markovic and I. Petrovic, "Bearings-only tracking with a mixture of
%    von Mises distributions," in IEEE/RSJ International Conference on
%    Intelligent Robots and Systems, Vilamoura, Algarve, Portugal, 7-12
%    Oct. 2012, pp. 707-712.
%
%April 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

    mu=wrapRange(mu1+mu2,-pi,pi);

    AkProd=BesseliRatio(1,kappa1)*BesseliRatio(1,kappa2);
    kappa=BesseliRatioInv(1,AkProd);
end

function [mu,kappa]=params4TrigMoment(rho)
%%PARAMS4TRIGMOMENT Obtain the parameters (mu and kappa) of the von Mises
%           distribution such that the von Mises distribution matches the
%           given complex mean resultant length for the first moment.
%
%INPUTS: rho The (complex) mean resultant value for the first moment. Note
%            that abs(rho)<=1.
%
%OUTPUTS: mu, kappa The mean direction and the concentration parameter of
%                   the von Mises distribution with the specified mean
%                   resultant value.
%
%The moments of a von Mises distribution are given in Chapter 2.2.4 of [1].
%This function just implements an inverse.
%
%REFERENCES:
%[1] S. R. Jammalamadaka and A. SenGupta, Topics in Circular Statistics.
%    Singapore: World Scientific, 2001.
%
%August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

    mu=atan2(imag(rho),real(rho));
    kappa= BesseliRatioInv(1,abs(rho));
end

function vals=rand(N,mu,kappa)
%%RAND Generate von Mises random variables with the given parameters.
%
%INPUTS: N If N is a scalar, then rand returns an NXN matrix of random
%          variables. If N=[M,N1] is a two-element row vector, then rand
%          returns an MXN1 matrix of random variables.
%       mu The mean direction of the distribution. This is normally
%          between -pi and pi.
%    kappa The concentration of the distribution 0<=kappa<inf.
%
%OUTPUTS: vals A matrix whose dimensions are determined by N of the
%              generated von Mises random variables.
%
%The algorithm if that of Section 4 of [1], which is a rejection sampling
%method based on using a wrapped Cauchy distribution as a proposal density.
%However, the range of the random variable u1 in [1] does not appear to be
%correct. The correct formulation is based on Chapter 7.3 of [2] using the
%errata available at http://luc.devroye.org/errors.pdf so u1 is randomly
%generated between -1 and 1 as opposed to between 0 and 1.
%
%REFERENCES:
%[1] D. J. Best and N. I. Fisher, "Efficient simulation of the von Mises
%    distribution," Journal of the Royal Statistical Society. Seriec C
%    (Applied Statistics), vol. 28, no. 2, pp. 152-157, 1979.
%[2] L. Devroye, Non-Uniform random Variate Generation. New York:
%    Springer-Verlag, 1986.
%
%August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(isscalar(N))
        dims=[N N];
    else
        dims=N;
    end

    vals=zeros(dims);%Allocate space
    numSamp=prod(dims);

    %Step 0
    tau=1+sqrt(1+4*kappa^2);
    rho=(tau-sqrt(2*tau))/(2*kappa);
    r=(1+rho^2)/(2*rho);

    for curSamp=1:numSamp
        while(1)
            %Step 1
            u1=2*rand(1)-1;%In range [-1,1]
            z=cos(pi*u1);
            f=(1+r*z)/(r+z);
            c=kappa*(r-f);

            %Step 2
            u2=rand(1);%in range [0,1]
            if(c*(2-c)-u2>0)
               break;
            end

            %Step 3
            if(log(c/u2)+1-c>=0)
                break;
            end
        end
        %Step 4
        u3=2*rand(1)-1;%In range [-1,1]
        vals(curSamp)=wrapRange(sign(u3)*acos(f)+mu,-pi,pi);
    end
end

function entropyVal=entropy(kappa)
%%ENTROPY Obtain the differential entropy of the von Mises distribution
%         given in nats. The differential entropy of a continuous
%         distribution is entropy=-int_x p(x)*log(p(x)) dx where the
%         integral is over all values of x. Units of nats mean that the
%         natural logarithm is used in the definition. Unlike the Shannon
%         entropy for discrete variables, the differential entropy of
%         continuous variables can be both positive and negative.
%
%INPUTS: kappa The concentration of the distribution 0<=kappa<inf.
%
%OUTPUTS: entropyVal The value of the differential entropy in nats.
%
%Differential entropy is defined in Chapter 8 of [1].
%
%If kappa is very large, besseli(0,kappa) will overflow and a warning will
%be issued. In such an instance, an approximation based on the first term
%of the series on page 289 of [2] is used. The approximation becomes
%increasingly bad as kappa becomes larger.
%
%REFERENCES:
%[1] T. M. Cover and J. A. Thomas, Elements of Information Theory, 2nd ed.
%    Hoboken, NJ: Wiley-Interscience, 2006.
%[2] S. R. Jammalamadaka and A. SenGupta, Topics in Circular Statistics.
%    Singapore: World Scientific, 2001.
%
%April 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

    besseliVal=besseli(0,kappa);

    if(isfinite(besseliVal))
        entropyVal=-kappa*BesseliRatio(1,kappa)+log(2*pi*besseliVal);
    else
        %Use the approximation exp(kappa)/sqrt(2*pi*kappa) for values of
        %kappa so large that besseli overflows
        warning('besseli(0,kappa) overflowed; using an approximation.')
        entropyVal=-kappa*BesseliRatio(1,kappa)+(1/2)*(4000+log(pi/1000));
    end
end

end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%VGG_SOLVELIN_BLKSYM  Solves M*x==y where M is (typically huge sparse) symmetric 4-block matrix.
%   It solves the system much more efficiently than a general (sparse) linear system solver.
%   Typical usage to solve normal equations in Levenberg-Marquardt in bundle adjustment.
%
%   X = VGG_SOLVELIN_BLKSYM(A,B,C,p,q [,'nocheck']) solves M*x==Y where
%   M=[A B; B' C] and y=[p;q].
%
%   X = VGG_SOLVELIN_BLKSYM(M,y,sideA [,'nocheck']) does the same but takes M and splits it,
%   M=[A B; B' C] and size(A)==[sideA sideA].
%
%   The function checks whether C is really near to diagonal; if not, a warning is printed.
%   Parmeter 'nocheck' switches off this test.

% (c) {awf,werner}@robots.ox.ac.uk, March 2002

function x = vgg_solvelin_blksym(varargin)

switch nargin
 case 6
  [A,B,C,p,q,check] = deal(varargin{:});
 case 5
  [A,B,C,p,q] = deal(varargin{:});
  check = '';
 case 4
  [M,y,sideA,check] = deal(varargin{:});
 case 3
  [M,y,sideA] = deal(varargin{:});
  check = '';
 otherwise
  error('Bad number of parameters');
end

if nargin<5
  Rtop = 1:sideA;
  Rbot = sideA+1:size(M,1);
  C = M(Rbot, Rbot);
  B = M(Rtop, Rbot);
  A = M(Rtop, Rtop);
  p = y(Rtop);
  q = y(Rbot);
end

if ~strcmp(check,'nocheck') & nnz(C(1,:))/size(C,1) > .1
  warning('It is likely that M was splitted incorrectly. Check diagonality of C and/or value of sideA.');
end

%tic

% Surprisingly, branch 1 is 2x slower than branch 2. We don't know why.
switch 2
 case 1
  invC_Btq = C \ [B' q];
 case 2
  invC_Btq = inv(C) * [B' q];
end

invC_Bt = invC_Btq(:,1:end-1);
invC_q = invC_Btq(:,end);

u = (A - B * invC_Bt) \ (p - B * invC_q);
v = invC_q - invC_Bt*u;

x = [u;v];

%toc

return

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
fprintf('\nHere we train a Convolutional RBM on the MNIST Dataset.\n');
dataset = 'mnistSmall';

load('mnistSmall.mat','trainData');
trainData = reshape(trainData',28,28,10000);

dataSize = [28,28,1];  % [nY x nX x nChannels]

% DEFINE AN ARCHITECTURE
arch = struct('dataSize', dataSize, ...
		'nFM', 9, ...
        'filterSize', [7 7], ...
        'stride', [2 2], ...
        'inputType', 'binary');

% GLOBAL OPTIONS
arch.opts = {'nEpoch', 1, ...
			 'lRate', .05, ...
			 'displayEvery',500, ...
			 'sparsity', .02, ...
			 'sparseGain', 5};%, ...
%  			 'visFun', @visBinaryCRBMLearning}; % UNCOMMENT TO VIEW LEARNING

% INITIALIZE AND TRAIN
cr = crbm(arch);
cr = cr.train(trainData);

% INFER HIDDEN AND POOLING LAYER EXPECTATIONS
% CONDITIONED ON SOME INPUT
[cr,ep] = cr.poolGivVis(trainData(:,:,1));
cr = cr.hidGivVis(trainData(:,:,1));

% DISPLAY NETWORK FEATURES
figure;
[nRows,nCols,nFM]=size(cr.W);
W = reshape(cr.W,nRows*nCols,nFM);
subplot(141);
visWeights(W,1);
title('Learned Filters');

subplot(142);
imagesc(trainData(:,:,1)); colormap gray; axis image; axis off
title(sprintf('Sample Data\nPoint'));

subplot(143);
[nCols,nRows,k]=size(cr.eHid);
eh = reshape(cr.eHid,nRows*nCols,k);
visWeights(eh);
title('Feature Maps')

subplot(144);
[eY,eX,k] = size(ep);
visWeights(reshape(ep,eY*eX,k)); colormap gray
title(sprintf('Pooling Layer\nExpectations'))
drawnow

author: dustinstansbury

repo: medal


text:
function t=collision_time(s,t,x,T);
% finds the collsion time between x(1) at t(1) with speed s(s) and x(2)
% at t(2) with speed s(2).
%
if nargin<4
  T=1e6;
end;
sl=s(1); sr=s(2);
tl=t(1); tr=t(2);
xl=x(1); xr=x(2);
if sr > sl,
  t=T+1;
elseif sr < sl,
  h1=sl*tl;
  h2=sr*tr;
  h3=xl-xr;
  h4=sl-sr;
  h1=(h1-h2-h3);
  t=h1/h4;
elseif xl==xr&tl==tr,
  t=tl;
else
  t=T+1;
end;

repo: Aero-matlab


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% SPHARM-MAT is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function SpharmMatUtilDisplayObjs(confs, objs, cpath)

numSbj = size(objs,2);

tplt = [];
if exist(confs.Template,'file')
    if strcmp(confs.Template(end-3:end),'.mat')==1
        tplt = load(confs.Template);
        if ~isfield(tplt,'fvec')
            disp('Template file does not contain fvec');
            tplt = [];
        end
    else
        disp('Template should be a .mat file');
    end
end

for i = 1:numSbj
    file = objs{i};
    [pa,na,ex]=fileparts(file);
    nastr = strrep(na,'_','-');

    clear adc; clear cdata; rmsd = [];

    suffix = '';

    load(file);
    if exist('bim','var')
        roi = bim;
        DIM = size(roi);

        % Make binary image
        ind = find(roi>0); roi(ind) = 1;
        ind = find(roi<1); roi(ind) = 0;

        % Generate voxel surface for the binary image
        [vertices, faces] =  gen_surf_data(roi,origin,vxsize);
    end

    h = figure('Name', ['Object: ' na],'NumberTitle', 'off');
    cameratoolbar(h, 'Show');

    % calculate SPHARM reconstruction
    if exist('fvec','var')
        % get regular mesh in the parameter space: sph_verts, faces
        switch strtrim(char(confs.Mesh))
            % to-do-item: add an option for reconstruction using original mesh
            case 'quad32'
                load(fullfile(cpath,'ParamMeshes/R32_quad.mat'));
            case 'quad64'
                load(fullfile(cpath,'ParamMeshes/R64_quad.mat'));
            case 'quad128'
                load(fullfile(cpath,'ParamMeshes/R128_quad.mat'));
            case 'quad256'
                load(fullfile(cpath,'ParamMeshes/R256_quad.mat'));
            case 'quad512'
                load(fullfile(cpath,'ParamMeshes/R512_quad.mat'));
            case 'icosa1'
                load(fullfile(cpath,'ParamMeshes/L1_icosa.mat'));
            case 'icosa2'
                load(fullfile(cpath,'ParamMeshes/L2_icosa.mat'));
            case 'icosa3'
                load(fullfile(cpath,'ParamMeshes/L3_icosa.mat'));
            case 'icosa4'
                load(fullfile(cpath,'ParamMeshes/L4_icosa.mat'));
            case 'icosa5'
                load(fullfile(cpath,'ParamMeshes/L5_icosa.mat'));
            case 'icosa6'
                load(fullfile(cpath,'ParamMeshes/L6_icosa.mat'));
        end

        % calculate SPHARM reconstruction using the regular mesh: vertices
        if ~strcmp(deblank(char(confs.Mesh)), 'orig')
            % get the user-specified degree
            if isempty(confs.Degree)
                degree = Inf;
            else
                degree = confs.Degree;
            end
            degree = min(degree,sqrt(size(fvec,1))-1);
            suffix = sprintf('%s_d%d',strtrim(char(confs.Mesh)),degree);
            % do the calculation
            Z = calculate_SPHARM_basis(sph_verts, degree);
            lb = 1;
            ub = (degree+1)^2;
            vertices = real(Z(:,lb:ub)*fvec(lb:ub,:));
            % calculate RMSD if there is a template
            if ~isempty(tplt)
                tfvec = tplt.fvec;
                tcnt = size(tfvec,1);
                if tcnt>=ub
                    tfvec = tfvec(1:ub,:);
                else
                    tfvec = zeros(1:tcnt,3);
                    tfvec(1:tcnt,:) = tplt.fvec;
                end
                rmsd = SPHARM_rmsd(fvec,tfvec);
            end
        end
    end

    % option for display in object space, parameter space, or both
    switch strtrim(char(confs.Space))
        case 'object'
            if ~exist('vertices','var') | ~exist('faces','var')
                disp(sprintf('Cannot display %s in the object space',nastr));
                break;
            end
            subp = 0;
            titleStr = sprintf('Object: %s',nastr);
        case 'param'
            if ~exist('sph_verts','var') | ~exist('faces','var')
                disp(sprintf('Cannot display %s in the parameter space',nastr));
                break;
            end
            subp = 1;
            titleStr = sprintf('Parameterization: %s',nastr);
        case 'both'
            if ~exist('faces','var') | (~exist('vertices','var') & ~exist('sph_verts','var'))
                disp(sprintf('Cannot display %s in object and parameter spaces',nastr));
                break;
            end
            subp = 2;
            titleStrLeft = sprintf('%s',nastr);
            titleStrRight = sprintf('Parameterization');
    end

    % option to visualize parameterization in color
    % ADC is shown together with this option
    over1 = 0;
    switch strtrim(char(confs.Overlay))
        case 'none'
            overl = 0;
        case 'adc_paramap'
            if exist('vertices','var') & exist('sph_verts','var') & exist('faces','var')
                overl = 1;
                % calculate the color coding
                [THETA,PHI,R] = cart2sph(sph_verts(:,1),sph_verts(:,2),sph_verts(:,3));
                ix = find(THETA<0); THETA(ix) = THETA(ix)+pi*2;
                cdata = THETA.*sign(PHI);
                % calculate the area distortion cost: adc
                asr = calc_asr(vertices, sph_verts, faces); adc = asr.prmstc;
            else
                disp(sprintf('Color coded parameterization cannot be calculated for %s due to missing variables',nastr));
            end
    end

    % option for solid shading, mesh view, or both
    switch strtrim(char(confs.Shade))
        case 'solid'
            if overl == 1
                if subp == 0
                    patch_color(vertices, faces, cdata); title(titleStr);
                elseif subp == 1
                    patch_color(sph_verts, faces, cdata); title(titleStr);
                elseif subp == 2
                    subplot(1,2,1);
                    patch_color(vertices, faces, cdata); title(titleStrLeft);
                    subplot(1,2,2);
                    patch_color(sph_verts, faces, cdata); title(sprintf('%s: ADC=%f',titleStrRight,adc));
                end
            else
                if subp == 0
                    patch_light(vertices, faces); title(titleStr);
                elseif subp == 1
                    patch_light(sph_verts, faces); title(titleStr);
                elseif subp == 2
                    subplot(1,2,1);
                    patch_light(vertices, faces); title(titleStrLeft);
                    subplot(1,2,2);
                    patch_light(sph_verts, faces); title(titleStrRight);
                end
            end
        case 'mesh'
            if overl == 1
                disp('Color coded parameterization is not shown in mesh view');
                if subp == 0
                    patch_mesh(vertices, faces); title(titleStr);
                elseif subp == 1
                    patch_mesh(sph_verts, faces); title(titleStr);
                elseif subp == 2
                    subplot(1,2,1);
                    patch_mesh(vertices, faces); title(titleStrLeft);
                    subplot(1,2,2);
                    patch_mesh(sph_verts, faces); title(sprintf('%s: ADC=%f',titleStrRight,adc));
                end
            else
                if subp == 0
                    patch_mesh(vertices, faces); title(titleStr);
                elseif subp == 1
                    patch_mesh(sph_verts, faces); title(titleStr);
                elseif subp == 2
                    subplot(1,2,1);
                    patch_mesh(vertices, faces); title(titleStrLeft);
                    subplot(1,2,2);
                    patch_mesh(sph_verts, faces); title(titleStrRight);
                end
            end
        case 'both'
            if overl == 1
                if subp == 0
                    patch_colormesh(vertices, faces, cdata); title(titleStr);
                elseif subp == 1
                    patch_colormesh(sph_verts, faces, cdata); title(titleStr);
                elseif subp == 2
                    subplot(1,2,1);
                    patch_colormesh(vertices, faces, cdata); title(titleStrLeft);
                    subplot(1,2,2);
                    patch_colormesh(sph_verts, faces, cdata); title(sprintf('%s: ADC=%f',titleStrRight,adc));
                end
            else
                if subp == 0
                    patch_lightmesh(vertices, faces); title(titleStr);
                elseif subp == 1
                    patch_lightmesh(sph_verts, faces); title(titleStr);
                elseif subp == 2
                    subplot(1,2,1);
                    patch_lightmesh(vertices, faces); title(titleStrLeft);
                    subplot(1,2,2);
                    patch_lightmesh(sph_verts, faces); title(titleStrRight);
                end
            end
    end

    if ~isempty(suffix)
        suffix = sprintf('%s_%s',confs.Space(1:3),suffix);
    else
        suffix = confs.Space(1:3);
    end
    if ~isempty(rmsd)
        zstr = sprintf('%s, rmsd=%0.3f',suffix,rmsd);
        zlabel(strrep(zstr,'_','-'));
    else
        zlabel(strrep(suffix,'_','-'));
    end

    % display and/or export
    switch strtrim(char(confs.Export))
        case 'screen'
            % display the figure
        case 'png'
            % save as png
            pa_png = fullfile(pa,'PNG');
            if ~exist(pa_png,'dir')
                mkdir(pa_png);
            end
            saveas(h,fullfile(pa_png,[na '_' suffix]),'png'); close(h);
        case 'both'
            % do both
            pa_png = fullfile(pa,'PNG');
            if ~exist(pa_png,'dir')
                mkdir(pa_png);
            end
            saveas(h,fullfile(pa_png,[na '_' suffix]),'png');
    end

end

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
classdef prtActionWrapperMil < prtAction

    properties (SetAccess = private)
        % Descriptive name of prtAction object.
        name = 'MIL Wrapper';

        % Shortened name for the prtAction object.
        nameAbbreviation = 'MILWrap';
    end

    properties (SetAccess = protected)
        % Specifies if the prtAction requires a labeled dataSet
        isSupervised

        % Indicates whether or not cross-validation is a valid operation
        isCrossValidateValid
    end

    properties (Dependent)
        theAction; % a prtActionObject
    end
    properties (SetAccess = private,Hidden)
        theAction_ = prtPreProcPca;
    end

    methods
        function self = prtActionWrapperMil(varargin)
            self = prtUtilAssignStringValuePairs(self,varargin{:});
            self.isSupervised = self.theAction.isSupervised;
            self.isCrossValidateValid = self.theAction.isCrossValidateValid;
        end

        function val = get.theAction(self)
            val = self.theAction_;
        end

        function self = set.theAction(self,action)

            if isa(action,'prtAction');
                self.theAction_ = action;

                self.isSupervised = self.theAction.isSupervised;
                self.isCrossValidateValid = self.theAction.isCrossValidateValid;

            else
                error('prt:prtActionWrapperMil','theAction must be a subclass of prtAction, provided action was a %s',class(action));
            end
        end
    end

    methods (Access = protected, Hidden = true)
        % prtAction.trainAction() - Primary method for training a prtAction
        %   self = prtAction.trainAction(self,DataSet)
        function self = trainAction(self, dsMil)

            ds = prtDataSetClass(dsMil.expandedData,dsMil.expandedTargets);
            self.theAction = self.theAction.train(ds);

        end
        % prtAction.runAction() - Primary method for evaluating a prtAction
        %   DataSet = runAction(self, DataSet)
        function dsMil = runAction(self, dsMil)

            ds = prtDataSetClass(dsMil.expandedData,dsMil.expandedTargets);
            dsOut = self.theAction.run(ds);
            dsMil = dsMil.fromPrtDataSetClassData(dsOut);

        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text: %% BinghamS2

author: mtex-toolbox

repo: mtex

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


text:
function [ cluster, cluster_center, cluster_energy ] = ...
  cluster_energy_compute ( dim_num, point_num, cluster_num, point, ...
  cluster_center )

%*****************************************************************************80
%
%% CLUSTER_ENERGY_COMPUTE computes the energy of the clusters.
%
%  Discussion:
%
%    The cluster energy is defined as the sum of the distance
%    squared from each point to its cluster center.  It is the goal
%    of the H-means and K-means algorithms to find, for a fixed number
%    of clusters, a clustering that minimizes this energy
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the number of spatial dimensions.
%
%    Input, integer POINT_NUM, the number of data points.
%
%    Input, integer CLUSTER_NUM, the number of clusters.
%
%    Input, real  POINT(DIM_NUM,POINT_NUM), the data points.
%
%    Input, integer CLUSTER(POINT_NUM), the cluster to which each
%    data point belongs.
%
%    Input, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM), the centers.
%
%    Output, real CLUSTER_ENERGY(CLUSTER_NUM), the energy
%    associated with each cluster.
%
  cluster_energy(1:cluster_num) = 0.0;

  for i = 1 : point_num

    j = cluster(i);

    point_energy = sum ( ...
      ( point(1:dim_num,i) - cluster_center(1:dim_num,j) ).^2 );

    cluster_energy(j) = cluster_energy(j) + point_energy;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


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

author: petercorke

repo: robotics-toolbox-matlab


text:
%% Toolox Image - A Toolbox for General Purpose Image Processing
%
% Copyright (c) 2008 Gabriel Peyre
%

%%
% The toolbox can be downloaded from Matlab Central
% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=16201&objectType=FILE

%%
% This includes in the path some additional useful scripts.
path(path, 'toolbox/');

%% Image loading and displaying

%%
% The |load_image| function contains several synthetic examples and can
% also load arbitrary images. You can provide the size of the image, and
% you can load batch of images using cell arrays. Optional parameters for
% the synthetic images are given using the |options| field.

options.radius = .3;
namelist = {'lena' 'disk' 'chessboard' 'phantom'};
Mlist = load_image(namelist, 256, options);

%%
% The |imageplot| function allows to display sets of images together and
% link the axes. It handles color images.

clf;
imageplot(Mlist,namelist, 2,2);

%% Image Manipulation

%%
%  You can rotate an image.

M = load_image('lena', 256);
Mlist = {};
for i=1:6
    Mlist{end+1} = perform_image_rotation(M,i*pi/8);
end
imageplot(Mlist);

%%
% You can (gaussian) blur an image.

M = load_image('lena', 256);
Mlist = {};
for i=1:6
    Mlist{end+1} = perform_blurring(M,i*4);
end
imageplot(Mlist);

%%
% You can crop and resize images.

M = load_image('lena');
M = crop(M,32);
Mi = perform_image_resize(M,[256 256]);
clf; imageplot({M Mi}, {'Original' 'Interpolated'});

%% Image Registration

%%
% You can perform image registration by first selecting two points on the image
% (here using clicking).

M = load_image('barb');
% extract the eyes
clf; imagesc([0,1],[0,1],M);
title('Click on each eye');
colormap gray(256);
axis image; axis off;
[x,y,b] = ginput(2);

%%
% Then a similitude is applied to register the eye position.

% position of the eye
u = [x(1) y(1)];
v = [x(2) y(2)];
% target position
u1 = [0.3,0.3];
v1 = [0.7,0.3];
% registration
M1 = perform_image_similitude(M,u,u1,v,v1);
% display
clf; imageplot({M M1}, {'Original' 'Registered'});

%% Spacially varying filter
% You can provide a set of filters, and for each pixel decide which filter
% to use (for instance according to local image information). The code is implemented with
% a mex C code, so it is quite fast. Here we show
% a simple foveated filtering, where the blurring is stronger in the center
% of the image.

n = 256;
M = load_image('lena', n);

%%
% We first create the set of filters
m = 41; % width of filers
p = 20; % number of filters
sigma = linspace(0.05,10,p);
H = zeros(m,m,p);
for i=1:p
    H(:,:,i) = compute_gaussian_filter([m m],sigma(i)/n,[n n]);
end

%%
% Then we compute an index map that tells which filter to use (foveation effect).
x = linspace(-1,1,n);
[Y,X] = meshgrid(x,x);
R = sqrt(X.^2 + Y.^2);
I = round(rescale(R,1,p));

%%
% At last, we launch the filtering, and display the result.

M1 = perform_adaptive_filtering(M,H,I);
imageplot({M M1},{'Original' 'Foveated'});

%% Median Filtering
% Median filter is effective to remove salt and pepper impulse noise.

n = 256;
M = load_image('lena');
M = rescale(crop(M, n));
% add noise at random locations
r = .3;  % 30% location corrupted
I = randperm(n^2); I = I(1:round(r*end));
M(I) = rand(size(I));
% do the filtering
M1 = perform_median_filtering(M,2);
% display
imageplot({M M1}, {'Original' 'Filtered'});

%% Line integral convolution
% The LIC algorithm filters an image along a vector field.

%%
% create a vector field (this is done using the functions of the diffc
% toolbox).

n = 256;    % size of the image
sigma = 60; % regularity of the vector field
options.bound = 'sym'; % boundary handling
v = perform_blurring(randn(n,n,2), sigma, options);
v = perform_vf_normalization(v);

%%
% Now we perform the LIC of an initial noise image along the flow.

M = randn(n);
% parameters for the LIC
options.histogram = 'linear'; % keep contrast fixed
options.verb = 0;
options.dt = 1.5; % time steping
% size of the features
options.flow_correction = 1;
options.niter_lic = 2; % several iterations gives better results
% iterated lic
Mlist = {};
for i=1:4
    options.M0 = M;
    Mlist{end+1} = perform_lic(v, i*4, options);
end
% display
clf; imageplot(Mlist,'',2,2);

%% TV denoising
% The Sparsity Toolbox implements TV denoising, that minimizes the total variation of
% the image.

%%
% |argmin_f |F-f|^2 + 2*lambda*TV(f)|

%%
% Where |F| is the noisy image and |f| is the optimized denoised image.
% |TV(f)| is the discrete total variation of the image. |lambda| controls
% how much denoising you want to perform.

%%
% It corresponds to the Rudin/Osher/Fatemi image denoiser
% (lagrangian formulation), and works well for images with strong edges
% (otherwise, wavelets performs better). The TV minimization is performed
% using the iterative fixed point algorihtm of Chambolle (published in
% JMIV 2002). A nice feature of the algorithm is that you can give him
% the target residual value (norm of the noise) in |options.etgt| and it will track it for
% you during the iterations, and will automatically select the correct |lambda| value.

%%
% First we load and image and make some noise.

n = 128;
M0 = load_image('lena',256);
M0 = rescale(crop(M0,n));
sigma = .1;
M = M0 + randn(size(M0))*sigma;
% display
clf; imageplot({clamp(M0) clamp(M)},{'Original' 'Noisy'});

%%
% We set the target residual error |options.etgt| just little larger than the noise level
% |n*sigma|, and run the algorithm.

% some parameter for the algorithm
options.verb = 0;
options.display = 0;
options.niter = 300;    % number of iterations
options.etgt = 1.1*sigma*n;
options.lambda = .3; % initial regularization
% now we perform the denoising
[Mtv,err,tv,lambda] = perform_tv_denoising(M,options);

%%
% We display the result.

clf; imageplot({clamp(M) clamp(Mtv)},{'Noisy' 'TV Denoised'});

%%
% We can keep track of the decay of |lambda| during the iteration and the
% Lagrangian decay.

clf;
subplot(2,1,1);
h = plot(lambda); axis tight;
set(h, 'LineWidth', 2);
set(gca, 'FontSize', 20);
title('Evolution of \lambda');
subplot(2,1,2);
lagr = .5*err.^2 + lambda(end)*tv;
h = plot(lagr); axis tight;
set(h, 'LineWidth', 2);
set(gca, 'FontSize', 20);
title('Evolution of |y-U x|^2+\lambda |x|_1');
axis([1 options.niter min(lagr) min(lagr)*1.3]);

author: gpeyre

repo: matlab-toolboxes


text:
function f = bohach1 ( m, x )

%*****************************************************************************80
%
%% BOHACH1 evaluates the Bohachevsky function #1.
%
%  Discussion:
%
%    The minimizer is
%
%      X* = [ 0.0, 0.0 ]
%      F(X*) = 0.0
%
%    Suggested starting point:
%
%      X = [ 0.5, 1.0 ];
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Zbigniew Michalewicz,
%    Genetic Algorithms + Data Structures = Evolution Programs,
%    Third Edition,
%    Springer Verlag, 1996,
%    ISBN: 3-540-60676-9,
%    LC: QA76.618.M53.
%
%  Parameters:
%
%    Input, integer M, the number of variables.
%
%    Input, real X(M), the argument of the function.
%
%    Output, real F, the value of the function at X.
%
  f =       x(1) * x(1) - 0.3 * cos ( 3.0 * pi * x(1) ) ...
    + 2.0 * x(2) * x(2) - 0.4 * cos ( 4.0 * pi * x(2) ) ...
    + 0.7;

  return
end


text:
function monomial_test035 ( )

%*****************************************************************************80
%
%% MONOMIAL_TEST035 tests MONO_NEXT_GRLEX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MONOMIAL_TEST035\n' );
  fprintf ( 1, '  MONO_NEXT_GRLEX computes the next monomial\n' );
  fprintf ( 1, '  in M variables, in graded lexicographic order.\n' );

  m = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Let M =  %d\n', m );

  a = 0;
  b = 3;
  seed = 123456789;

  for i = 1 : 10

    [ x, seed ] = i4vec_uniform_ab ( m, a, b, seed );

    fprintf ( 1, '\n' );
    for j = 1 : m
      fprintf ( 1, '  %1d', x(j) );
    end
    fprintf ( 1, '\n' );

    for j = 1 : 5
      x = mono_next_grlex ( m, x );
      for j = 1 : m
        fprintf ( 1, '  %1d', x(j) );
      end
      fprintf ( 1, '\n' );
    end

  end

  return
end


text:
%% Colour Segmentation - VIBGYOR Colour Segmentation
%% This function can be used for VIBGYOR Colour segmentation from the RGB
%% Color images.
%% Function C = VIBGYORsegmentation(img)
%% input    img = Color image (The input image should be a color image)
%%  Example: C = VIBGYORsegmentation(img);
%%      Posted date : 14 - 07 - 2008
%%
%% Developed By : K.Kannan & Jeny Rajan
%%                  Medical Imaging Research Group (MIRG), NeST, Trivandrum.
%%
function C = VIBGYORsegmentation(img)
%% Electro-Magnetic wavelengths ranges in Nano Meter
% red 760-720
% orange 620-590
% yellow 590-545
% green 525-490
% blue 490-450
% indigo 450-420
% violet 420-380
%%
[row col plane] = size(img);
img = double(img);
C = zeros(row,col,plane);
if plane ~= 3
    disp('Input should be a color image');
    return;
end
GL = 255;
%% Color Choice
clc;
disp('   ');
disp('           r / R - for Red Color');
disp('           o / O - for Orange Color');
disp('           y / Y - for Yellow Color');
disp('           g / G - for Green Color');
disp('           b / B - for Blue Color');
disp('           i / I - for Indigo Color');
disp('           v / V - for Violet Color');
color = input('\nYour Color Choice = ','s');
%%
switch color
    case {'R','r'}
        f1 = (GL * 1);f2 = (GL * 0.4);f3 = (GL * 0.5);
        f4 = (GL * 0);f5 = (GL * 0.5);f6 = (GL * 0);
        C = cfilter(img,f1,f2,f3,f4,f5,f6,1,0);
    case {'O','o'}
        f1 = (GL * 1);f2 = (GL * 0.6);f3 = (GL * 0.68);
        f4 = (GL * 0.3);f5 = (GL * 0.3);f6 = (GL * 0);
        C = cfilter(img,f1,f2,f3,f4,f5,f6,1,0);
    case {'Y','y'}
        f1 = (GL * 1);f2 = (GL * 0.78);f3 = (GL * 1);
        f4 = (GL * 0.72);f5 = (GL * 0.52);f6 = (GL * 0);
        C = cfilter(img,f1,f2,f3,f4,f5,f6,1,1);
    case {'G','g'}
        f1 = (GL * 0.68);f2 = (GL * 0);f3 = (GL * 1);
        f4 = (GL * 0.4);f5 = (GL * 0.68);f6 = (GL * 0);
        C = cfilter(img,f1,f2,f3,f4,f5,f6,2,0);
    case {'B','b'}
        f1 = (GL * 0.5);f2 = (GL * 0);f3 = (GL * 0.68);
        f4 = (GL * 0);f5 = (GL * 1);f6 = (GL * 0.4);
        C = cfilter(img,f1,f2,f3,f4,f5,f6,3,0);
    case {'I','i'}
        f1 = (GL * 0.68);f2 = (GL * 0);f3 = (GL * 1);
        f4 = (GL * 0.6);f5 = (GL * 1);f6 = (GL * 0.6);
        C = cfilter(img,f1,f2,f3,f4,f5,f6,1,1);
    case {'V','v'}
        f1 = (GL * 1);f2 = (GL * 0.6);f3 = (GL * 0.68);
        f4 = (GL * 0);f5 = (GL * 1);f6 = (GL * 0.6);
        C = cfilter(img,f1,f2,f3,f4,f5,f6,3,1);
    otherwise
        disp('unknown method');
end
%% Display
C = uint8(C);
figure, subplot(1,2,1),imshow(uint8(img),[]);title('Original Image');
subplot(1,2,2),imshow(uint8(C),[]);title('Color Segmented Image');

%% Function for Color Filter
function C = cfilter(img,f1,f2,f3,f4,f5,f6,m,flg)
[row col plane] = size(img);
C = zeros(row,col,plane);
for i = 1:row
    for j = 1:col
        if flg == 0
            if (img(i,j,1) <= f1 && img(i,j,1) >= f2 && ...
                img(i,j,2) <= f3 && img(i,j,2) >= f4 && ...
                img(i,j,3) <= f5 && img(i,j,3) >= f6 ...
                && img(i,j,m) == max([img(i,j,1) img(i,j,2) img(i,j,3)]))
                C(i,j,1:3) = img(i,j,1:3);
            else
                C(i,j,1:3) = (img(i,j,1) * 0.3) + (img(i,j,2) * 0.59) + (img(i,j,3) * 0.11);
            end
        else
            if (img(i,j,1) <= f1 && img(i,j,1) >= f2 && ...
                img(i,j,2) <= f3 && img(i,j,2) >= f4 && ...
                img(i,j,3) <= f5 && img(i,j,3) >= f6)
                C(i,j,1:3) = img(i,j,1:3);
            else
                C(i,j,1:3) = (img(i,j,1) * 0.3) + (img(i,j,2) * 0.59) + (img(i,j,3) * 0.11);
            end
        end
    end
end
%%

author: Sable

repo: mcbench-benchmarks


text:
%vgg_F_from_7pts_2img  Computes fundamental matrix from 7 points across 2 images.
%
%   [P,X] = vgg_F_from_7pts_2img(x), where
%      x ... double(3,7,2) or cell{2} of double(3,7), 7 homogeneous points in 2 images
%      F ... double(3,3), fundamental matrix
%   There are 0 to 3 solutions for F. Solutions are pruned by requirement that
%   scalars s in all equations s*cross(e1,x1)==F*x2 are positive.
%   In case of multiple solutions, F has one dimension
%   more such that F(:,:,n) is the n-th solution.
%
%   Also the form F = vgg_F_from_7pts_2img(x1,x2) is accepted.

function F = vgg_F_from_7pts_2img(x1,x2)

if nargin==1
  if iscell(x1)
    x2 = x1{2};
    x1 = x1{1};
  else
    x2 = x1(:,:,2);
    x1 = x1(:,:,1);
  end
end
if any(size(x1)~=[3 7]) | any(size(x2)~=[3 7])
  error('Wrong size of input points.');
end

% Linear step
A = vgg_vec_swap(x1,x2)';
[u,s,v] = svd(A,0);
FF{1} = reshape(v(:,end-1),[3 3]);
FF{2} = reshape(v(:,end  ),[3 3]);

% Solving cubic equation and getting 1 or 3 solutions for F
a = vgg_singF_from_FF(FF);
F = [];
for i = 1:length(a)
  Fi = a(i)*FF{1} + (1-a(i))*FF{2};
  %for n = 1:7, disp(norm(x(:,n,1)'*Fi*x(:,n,2))), end  % test code
  if signs_OK(Fi,x1,x2)
    F = cat(3, F, Fi);
  end
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%

% Checks sign consistence of F and x
function OK = signs_OK(F,x1,x2)
[u,s,v] = svd(F');
e1 = v(:,3);
l1 = vgg_contreps(e1)*x1;
s = sum( (F*x2) .* l1 );
OK = all(s>0) | all(s<0);
return


text:
function [energy, kinetic, potential] = energy_3_link(z,P)
%[ENERGY, KINETIC, POTENTIAL] = ENERGY_3_LINK(Z,P)
%
%FUNCTION:  This function computes the energy of a double
%    pendulum.
%INPUTS:
%
%
%OUTPUTS:
%    energy = [1 X nTime] vector of total energy
%    kinetic = [1 X nTime] vector of kinetic energy
%    potential = [1 X nTime] vector of potential energy
%
%NOTES:
%    This file was automatically generated by writeEnergy.m

g  = P.g ; %gravity
m1 = P.m(1); % Link 1 mass
m2 = P.m(2); % Link 2 mass
m3 = P.m(3); % Link 3 mass
l1 = P.l(1); % Link 1 length
l2 = P.l(2); % Link 2 length
l3 = P.l(3); % Link 3 length
I1 = P.I(1); % Link 1 moment of inertia about its center of mass
I2 = P.I(2); % Link 2 moment of inertia about its center of mass
I3 = P.I(3); % Link 3 moment of inertia about its center of mass
d1 = P.d(1); % Link 1 distance between center of mass and parent joint
d2 = P.d(2); % Link 2 distance between center of mass and parent joint
d3 = P.d(3); % Link 3 distance between center of mass and parent joint

th1 = z(1,:);
th2 = z(2,:);
th3 = z(3,:);

dth1 = z(4,:);
dth2 = z(5,:);
dth3 = z(6,:);

kinetic = (I1.*dth1.^2)./2 + (I2.*dth2.^2)./2 + (I3.*dth3.^2)./2 + (m1.*(d1.^2.*dth1.^2.*cos(th1).^2 + d1.^2.*dth1.^2.*sin(th1).^2))./2 + (m2.*((dth1.*l1.*cos(th1) + d2.*dth2.*cos(th2)).^2 + (d2.*dth2.*sin(th2) + dth1.*l1.*sin(th1)).^2))./2 + (m3.*((dth1.*l1.*cos(th1) + dth2.*l2.*cos(th2) + d3.*dth3.*cos(th3)).^2 + (d3.*dth3.*sin(th3) + dth1.*l1.*sin(th1) + dth2.*l2.*sin(th2)).^2))./2;
potential = g.*m3.*(d3.*sin(th3) + l1.*sin(th1) + l2.*sin(th2)) + g.*m2.*(d2.*sin(th2) + l1.*sin(th1)) + d1.*g.*m1.*sin(th1);
energy = potential + kinetic;

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function p = mtexExamplePath
% returns the default path to EBSD data files
%
% Syntax
%   mtexEBSDPath
%   getMTEXpref('EBSDPath') % long form of mtexEBSDPath
%   setMTEXpref('EBSDPath',path) % change the path to EBSD files
%
% See also
% mtexDataPath mtexCifPath mtexEBSDPath mtexPoleFigurePath mtexODFPath
% mtexTensorPath

examplePath = getMTEXpref('ExamplePath');

if exist(examplePath,'dir')
  p = examplePath;
else
  error('Examples not found!');
end


text:
outputVideo = VideoWriter(fullfile('background-video.avi'));
outputVideo.FrameRate=12;
open(outputVideo);

A=reshape(A,[z1 z2 z3]); ov=imnorm(A);
% B=reshape(B,[z1 z2 z3]); ov=imnorm(B);
% Z=reshape(Z,[z1 z2 z3]); ov=imnorm(Z);

for ii = 1:size(ov,3)
    img = ov(:,:,ii);
    writeVideo(outputVideo,img);
end

close(outputVideo);

outputVideo = VideoWriter(fullfile('foreground-video.avi'));
outputVideo.FrameRate=12;
open(outputVideo);

% A=reshape(A,[z1 z2 z3]); ov=imnorm(A);
B=reshape(B,[z1 z2 z3]); ov=imnorm(B);
% Z=reshape(Z,[z1 z2 z3]); ov=imnorm(Z);

for ii = 1:size(ov,3)
    img = ov(:,:,ii);
    writeVideo(outputVideo,img);
end

close(outputVideo);

outputVideo = VideoWriter(fullfile('raw-video.avi'));
outputVideo.FrameRate=12;
open(outputVideo);

% A=reshape(A,[z1 z2 z3]); ov=imnorm(A);
% B=reshape(B,[z1 z2 z3]); ov=imnorm(B);
Z=reshape(Z,[z1 z2 z3]); ov=imnorm(Z);

for ii = 1:size(ov,3)
    img = ov(:,:,ii);
    writeVideo(outputVideo,img);
end

close(outputVideo);

disp('videos written.')
clear outputVideo img ov


text:
function [C]=skinColor()
C=[255 218*0.9 180*0.9]./255;
end

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function a = cheby_van2 ( n )

%*****************************************************************************80
%
%% CHEBY_VAN2 returns the CHEBY_VAN2 matrix.
%
%  Discussion:
%
%    CHEBY_VAN2 is the Chebyshev Vandermonde-like matrix.
%
%    The formula for this matrix has been slightly modified, by a scaling
%    factor, in order to make it closer to its inverse.
%
%  Formula:
%
%    A(I,J) = ( 1 / sqrt ( N - 1 ) ) * cos ( (I-1) * (J-1) * PI / (N-1) )
%
%  Example:
%
%    N = 4
%
%                 1      1           1           1
%    1/sqrt(3) *  1  COS(PI/3)   COS(2*PI/3) COS(3*PI/3)
%                 1  COS(2*PI/3) COS(4*PI/3) COS(6*PI/3)
%                 1  COS(3*PI/3) COS(6*PI/3) COS(9*PI/3)
%
%  Properties:
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    The entries of A are based on the extrema of the Chebyshev
%    polynomial T(n-1).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  if ( n == 1 )
    a(1,1) = 1.0;
    return
  end

  for i = 1 : n
    for j = 1 : n

      angle = ( ( i - 1 ) * ( j - 1 ) ) * pi / ( n - 1 );

      a(i,j) = cos ( angle );

    end
  end

  a(1:n,1:n) = a(1:n,1:n) / sqrt ( n - 1 );

  return
end


text:
function exact = p48_exact ( )

%*****************************************************************************80
%
%% P48_EXACT returns the exact integral for problem 48.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  exact = -4.0;

  return
end


text:
function value = p23_i4vec ( action, name, dim_num, value )

%*****************************************************************************80
%
%% P23_I4VEC sets or gets I4VEC parameters for problem 23.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION, the action.
%    'D' sets the object to a default value.
%    'G' means the current value of the object should be returned.
%    'S' means the input values of the object and its dimension should
%    be stored.
%
%    Input, string NAME, the name of the parameter.
%    'E' is the exponent vector.
%
%    Input, integer DIM_NUM, the dimension of the object.
%
%    Input/output, integer VALUE(DIM_NUM), the value of the object.
%
  persistent e
  persistent dim_num_save

  if ( size ( dim_num_save ) == 0 )
    dim_num_save = 0;
  end

  if ( dim_num_save ~= dim_num )
    dim_num_save = 0;
    e = [];
  end

  if ( dim_num_save == 0 )
    dim_num_save = dim_num;
  end

  if ( action == 'D' | action == 'd' )

    if ( name == 'E' | name == 'e' | name == '*' )
      e(1:dim_num) = 2;
    end

  elseif ( action == 'G' | action == 'g' )

    if ( name == 'E' | name == 'e' )
      value(1:dim_num) = e(1:dim_num);
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P23_I4VEC - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P23_I4VEC - Fatal error!' );
    end

  elseif ( action == 'S' | action == 's' )

    if ( name == 'E' | name == 'e' )
      e(1:dim_num) = value(1:dim_num);
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P23_I4VEC - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P23_I4VEC - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P23_I4VEC - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action = "%s".\n', action );
    error ( 'P23_I4VEC - Fatal error!' );

  end

  return
end


text:
function res = ge(a,b)
%GE           Implements  a >= b  for slopes, compares only a.r and b.r
%

% written  12/06/98     S.M. Rump
% modified 09/28/01     S.M. Rump  matrices and multi-dimensional arrays
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if ~isa(a,'slope')
    asize = size(a);
    bsize = b.size;
    res = ( a(:)>=b.r(:,1) );
  elseif ~isa(b,'slope')
    asize = a.size;
    bsize = size(b);
    res = ( a.r(:,1)>=b(:) );
  else
    asize = a.size;
    bsize = b.size;
    res = ( a.r(:,1)>=b.r(:,1) );
  end

  if ~isequal(asize,bsize) & ( prod(asize)~=1 ) & ( prod(bsize)~=1 )
    error('incompatible size for slope >= ')
  end

  if prod(asize)==1
    res = reshape(res,bsize);
  else
    res = reshape(res,asize);
  end

repo: OpenMAS


text:
function value = r8vec_eq ( n, a1, a2 )

%*****************************************************************************80
%
%% R8VEC_EQ is true if two R8VEC's are equal.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vectors.
%
%    Input, real A1(N), A2(N), two vectors to compare.
%
%    Output, logical VALUE, is TRUE if every pair of elements A1(I)
%    and A2(I) are equal, and FALSE otherwise.
%
  value = all ( a1(1:n) == a2(1:n) );

  return
end


text:
%compute dyheadsm_mm

function [data,units]=compute_dyheadsm_mm(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
dyheadsm_mm=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    dyheadsm_mm{1,i}=(trx(larva).yheadsm_mm(2:end)-trx(larva).yheadsm_mm(1:end-1))./trx(larva).dt;
end

units=parseunits('mm/s');
data=dyheadsm_mm;

author: kristinbranson

repo: JAABA