text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   [T, vx, ax, phi, phid]=directkinematic_5R(robot, q)
%
%   Direct kinematics for the 5R planar parallel robot.
%
%   This script comprises the computation of the position of the end
%   effector T as a function of q. As well, the speed and acceleration
%   of the end effector can be also computed. The speed of the passive
%   joints phid are also computed.
%
%
%   In general, given a value of active joint coordinates q1, q2,
%   there exists two possible solutions. The position of the end effector P(x,y) is returned as two different
%   homogeneous matrices T= [T1, T2] with two possible positions of P in
%   each. If the robot is placed at a singular point P (direct kinematic
%   singularity), then both solutions are coincident. If this is the case,
%   a differential movement
%
%   In addition, if qd (q dot, speed), or qdd (q dot dot, acceleration)
%   of the joint variables, is given, the vx (speed) and ax (acceleration)
%   in cartesian coordinates are also computed.
%
%   Given q1 and q2, in this 5R robot, two possible solutions are feasible.
%   In consequence, vx=[vx1, vx3] where each column stores the speed [vx, vy]
%   for each of the possible solutions
%
%   Please, note that, given a value of active joint coordinates q1, q2,
%   there exists different possible solutions for P()
%
%   Author: Arturo Gil Aparicio, arturo.gil@umh.es
%   Date: 13/09/2013
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE 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.
%
% ARTE 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 ARTE.  If not, see <http://www.gnu.org/licenses/>.

function [T, vx, ax, phi, phid]=directkinematic_5R(robot, q)
%close all

a1=eval(robot.robot1.DH.a);
a2=eval(robot.robot2.DH.a);

%Link lengths
l1=abs(a1(1));
l2=abs(a1(2));
l3=abs(a2(1));
l4=abs(a2(2));
L=robot.L;%2.5;

xA=cos(q(1));
yA=sin(q(1));

xB=cos(q(2))+L;
yB=sin(q(2));

%figure,
plot(0,0,'r.'), %hold
plot(xA, yA,'k.')
plot(xB, yB,'k.')
plot_line([0 0 0], [xA yA 0], 'r', 2)
plot_line([L 0 0], [xB yB 0], 'r', 2)
%plot(L, 0,'*g')

%draw circles around xA, yA and xB yB as a first approximation
x=-0.47:0.001:2;
rr=[];
for i=1:length(x),
    r=solve_poly(x(i), xA, yA, l2);
    rr=[rr r];
end
plot(x, real(rr(1,:)))
plot(x, real(rr(2,:)))

x=0.5:0.001:3.085;
rr=[];
for i=1:length(x),
    r=solve_poly(x(i), xB, yB, l3);
    rr=[rr r];
end
plot(x, real(rr(1,:)), 'b')
plot(x, real(rr(2,:)), 'b')

%compute polynomial in terms of alpha, beta and delta
af=xA-xB;
be=l2^2-l1^2-l3^2+xB^2+yB^2;
de=(yB-yA);

a=(de^2+af^2);
b=af*be-2*xA*de^2-2*yA*af*de;
c=de^2*(xA^2+yA^2)+be^2/4-yA*be*de-l2^2*de^2;

%two different solutions for x
x1=(-b+sqrt(b^2-4*a*c))/(2*a);
x2=(-b-sqrt(b^2-4*a*c))/(2*a);

%find y from polynomial, given x1 and x2
r1=solve_poly(x1, xA, yA, l2);
r2=solve_poly(x2, xA, yA, l2);
r3=solve_poly(x1, xB, yB, l3);
r4=solve_poly(x2, xB, yB, l3);

%find those solutions that comply with both equations
R1=[r1; r2];
R2=[r3; r4];
%to do this, obtain only those roots that are repeated
y=[];
for i=1:length(R1),
   val=is_double_root(R1(i),R2);
   if val==1
       y=[y R1(i)];
   end
end

%Plot error if an unfeasible solution is found
if (~isreal(y))
    disp('ERROR: directkinematic_5R: unfeasible solution');
end

y1=y(1);
y2=y(2);

%return solutions in homogeneous matrices
T1=eye(4);
T2=eye(4);

T1(1,4)=x1;
T1(2,4)=y1;
T2(1,4)=x2;
T2(2,4)=y2;
T=[];
T=[T1 T2];

% avoid warnings for imaginary parts when plottin
% the above message should be triggered when an unfeasible solution is
% found
y1=real(y(1));
y2=real(y(2));
x1=real(x1);
x2=real(x2);

%plot solutions
plot_line([xA yA 0], [x1 y1 0], 'y', 2)
plot_line([xB yB 0], [x1 y1 0], 'y', 2)

plot_line([xA yA 0], [x2 y2 0], 'g', 2)
plot_line([xB yB 0], [x2 y2 0], 'g', 2)

% in this case, q1, q2, qd1, qd2
if length(q)==4
    fprintf('directkinematic_5R:: Computing speeds');
    q1=q(1);
    q2=q(2);
    qd1=q(3);
    qd2=q(4);
    %compute direct kinematic speeds,
    %given joint values q1 and q2 and joint speeds qd1 and qd2,
    %compute the speed
    [vx1, phi1, phid1]=compute_direct_speeds(robot, [q1 q2], [qd1 qd2], T1);
    [vx2, phi2, phid2]=compute_direct_speeds(robot, [q1 q2], [qd1 qd2], T2);
    vx=[vx1, vx2];
    phi=[phi1, phi2]; %position of passive joints
    phid=[phid1, phid2];%speed of passive joints
end

ax=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Build de coefficients of a polynomial as a function of the
%   points xA, yA and the length of the link L.
%   Solve for second order polynomial:
%   ax^2+bx+c=0;
%
%   as x=(-b +- sqrt(b^2-4*a*c))/(2*a)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function r=solve_poly(x, xAB, yAB, L)
a=1;
b=-2*yAB;
c=x^2-2*x*xAB-L^2+xAB^2+yAB^2;

r =[(-b+sqrt(b^2-4*a*c))/(2*a); (-b-sqrt(b^2-4*a*c))/(2*a)];

function yes_no=is_double_root(R1_i,R2)
thres=0.0001;
[j, k, v]=find(abs(R2-R1_i)<thres);

%in case a match has been found
if length(j)>0
        yes_no=1;
else
        yes_no=0;
end

function plot_line(p0, p1, color, w)
x0 = p0(1);
y0 = p0(2);
z0 = p0(3);
x1 = p1(1);
y1 = p1(2);
z1 = p1(3);
% Draw a line between p0 and p1
plot3([x0;x1],[y0;y1],[z0;z1], color, 'LineWidth',w);

function [vx, phi, phid]=compute_direct_speeds(robot, q, qd, T)
% caution, q1 and q2 are the active coordinates,
% the compute jacobians requires every of the 4 joint values
% in this order q={q1, phi1, q2, phi2}

q1=q(1);
q2=q(2);
qd1=qd(1);
qd2=qd(2);

%yes, compute the inverse kinematic to find the passive joint values that
%correspond to the position specified by T. We choose the only solution
% in which q1 and q2 are equal to the given values
qinv=inversekinematic(robot, T);
for i=1:4,
    qq=qinv(:,i);
    %simple threshold are used to find the most similar solution.
    if(abs(qq(1)-q1)<0.05) && (abs(qq(3)-q2)<0.05)
        % yes, we have used the inversekinematic function to compute
        % the phi1 and phi2 angles
        phi1=qq(2);
        phi2=qq(4);
        [JX, Jq, Jphi]=compute_jacobians_5R(robot, [q1 phi1 q2 phi2]);

        if det([JX Jphi])==0
            fprintf('\n WARNING: directkinematic_5R: direct kinematic singularity detected');
        end
        %Now that the Jacobians are known, compute vx and Phid
        [vxphid]=-inv([JX Jphi])*Jq*[qd1; qd2];
        vx=vxphid(1:2);
        phid=vxphid(3:4);
        phi=[phi1, phi2]';
        return;
    end
end

repo: ARTE


text:
function computeAllFinalRF_HN(pathRF,nBoot,seed)

startpath = pwd;

cd(pathRF), load('training')
nameOutcomes = fieldnames(training.outcomes); nOutcomes = numel(nameOutcomes);
fSetNames = fieldnames(training.textures.(nameOutcomes{1})); nFset = numel(fSetNames);

for o = 1:nOutcomes
    if strcmp(nameOutcomes{o},'DeathSign')
        outcome = training.outcomes.Death;
    else
        outcome = training.outcomes.(nameOutcomes{o});
    end
    for f = 1:nFset
        text = training.textures.(nameOutcomes{o}).(fSetNames{f}); nText = size(text,2);
        indClinic = training.clinical.bestAdd.(nameOutcomes{o}).(fSetNames{f});
        cost = training.cost.(nameOutcomes{o}).(fSetNames{f});
        tableTrain = [text,training.clinical.table(:,indClinic)];
        cat = logical([zeros(1,nText),training.clinical.categories(indClinic)]);
        rng(seed), [RF] = trainRF_table(tableTrain,outcome,cat,nBoot,cost);
        RF = compact(RF); % Compact version
        save(['RF_',[fSetNames{f},'clinic'],'_',nameOutcomes{o}],'RF')
    end
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
% File: lognormpdf.m
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

function [log_prob] = lognormpdf(x,mu,sigma)

% LOGNORMPDF Natural logarithm of the normal probability density function (pdf)
% Y = lognormpdf(X,MU,SIGMA) returns the log of the pdf of the normal
% distribution parameterized by mean MU and standard deviation SIGMA evaluated
% at each value in the vector X. Thus, the size of the return
% vector Y is the size of X.
%
% MU and X should have the same dimensions.

log_prob = -log(sigma*sqrt(2*pi))-(x-mu).^2 ./ (2*sigma.^2);

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function [A, B] = gsp_filterbank_bounds(G,W,param)
%GSP_FILTERBANK_BOUNDS Compute approximated frame bounds for a filterbank
%   Usage: [A, B] = gsp_filterbank_bounds(G,W);
%          [A, B] = gsp_filterbank_bounds(G,W,param);
%          [A, B] = gsp_filterbank_bounds([xmin, xmax],W);
%          [A, B] = gsp_filterbank_bounds([xmin, xmax],W,param);
%
%   Input parameters
%       G   : Graph structure or interval to compute the bound
%       W   : Filterbank (cell array of inline function)
%       param: optional parameter
%   Output parameters
%       A   : Filterbank lower bound
%       B   : Filterbank Upper bound
%
%   *param* is a Matlab structure containing the following fields:
%
%   * *param.N* : Number of points for the line search default (default 999)
%   * *param.use_eigenvalues* : Use eigenvalues if possible (default 1). To
%     be used, the eigenvalues have to be computed first using
%     |gsp_compute_fourier_basis|.
%

% Author: Nathanael Perraudin
% Date : 26 March 2014

if nargin < 3
    param = struct;
end

if ~isfield(param,'N'), param.N = 999; end
if ~isfield(param,'use_eigenvalues'),  param.use_eigenvalues = 1; end

if iscell(G)
    NG = numel(G);
    A = cell(NG,1);
    B = cell(NG,1);
    for ii = 1:NG
        [A{ii}, B{ii}] = gsp_filterbank_bounds(G{ii},W{ii},param);
    end
    return
end

if isstruct(G)
    if ~isfield(G,'lmax')
        G = gsp_estimate_lmax(G);
    warning(['GSP_FILTERBANK_BOUNDS: To be more efficient you should run: ',...
        'G = gsp_estimate_lmax(G); before using this proximal operator.']);
    end
    xmax = G.lmax;
    xmin = 0;
else
    xmin = G(1);
    xmax = G(2);
end

if param.use_eigenvalues && isstruct(G) && isfield(G,'e')
    lambda = G.e;
else
    lambda = linspace(xmin,xmax , param.N);
end

Nf = numel(W);

sum_filters = sum(abs(gsp_filter_evaluate(W,lambda).^2),2);
% for ii=1:Nf
%     sum_filters = sum_filters + (W{ii}(lambda)).^2;
% end

A = min(sum_filters);
B = max(sum_filters);

end

repo: gspbox


text:
function value = year_length_months_julian ( y )

%*****************************************************************************80
%
%% YEAR_LENGTH_MONTHS_JULIAN returns the number of months in a Julian year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, integer VALUE, the number of
%    months in the year.
%
  value = 12;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function L= ShRec(shCoeff)
% Generates bandpass data as per the passed shearlet coefficient
%Input:
%        shCoeff    : Shearlet Coefficient as obtained in Rec phase
%Output:
%        L          : Reconstructed BandPass Data . Reconstuction
%                   : Can be Partial as per passed shearlet coeifficient
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function L= ShRec(shCoeff)
level = length(shCoeff);
L=cell(1,level);

for l= 1:level
   [sizel2,sizel1] =size(shCoeff{l});
   D=zeros(size(shCoeff{l}{1,1}));
   for  l1= 1:sizel1
     for l2 =1:sizel2
      D=D+shCoeff{l}{l2,l1};
     end
   end
   L{l}=D;
end

author: thomaskuestner


text:
%  FILE:   cnn_widerface.m
%
%    This function serves as the main function for evaluating a model. It loads
%    detection and ground truth and evaluates at all three difficulty levels at
%    once.
%
%    Note that this evaluation script consistently produces slightly lower
%    numbers (mAPs) comparing to the official evaluation script. I implement my
%    own version because it runs faster and it is more customizable, in terms of
%    adding visualization etc.
%
%  INPUT:  configuration (see code for details)
%
%  OUTPUT: none

function cnn_widerface_eval(varargin)

%
opts.noUpsampling = false;
opts.noOrgres = false;
opts.noDownsampling = false;
opts.testMultires = false;
opts.testTag = '';
opts.inputSize = [500, 500];
[opts, varargin] = vl_argparse(opts, varargin);

%
opts.evalWithReg = true;
opts.tag = '';
opts.bboxReg = true;
opts.skipLRMult = [1, 0.001, 0.0001, 0.00001];
opts.penalizeDuplicate = true;
opts.nmsThresh = 0.3;
opts.probThresh = 0.8;
opts.draw = false;
opts.minOverlap = [0.3:0.2:0.7];
opts.testEpoch = 0;
opts.testIter = 0;
opts.testSet = 'val';
opts.resDir = 'results/';
opts.testName = '';
opts.clusterNum = 25;
opts.clusterName = '';
opts.sampleSize = 256;
opts.posFraction = 0.5;
opts.posThresh = 0.7;
opts.negThresh = 0.3;
opts.pretrainModelPath = 'matconvnet/pascal-fcn8s-tvg-dag.mat';
opts.dataDir = fullfile('data','widerface') ;
opts.modelType = 'pascal-fcn8s-tvg-dag' ;
opts.networkType = 'dagnn' ;
opts.batchNormalization = true ;
opts.weightInitMethod = 'gaussian' ;
[opts, varargin] = vl_argparse(opts, varargin) ;

%
opts.minClusterSize = [2, 2];
opts.maxClusterSize = opts.inputSize;
[opts, varargin] = vl_argparse(opts, varargin) ;

%
sfx = opts.modelType ;
sfx = [sfx '-' 'sample' num2str(opts.sampleSize)] ;
sfx = [sfx '-' 'posfrac' num2str(opts.posFraction)] ;
sfx = [sfx '-' 'N' num2str(opts.clusterNum)];

if opts.bboxReg, sfx = [sfx '-' 'bboxreg']; end

if any(opts.inputSize~=500),
    sz = opts.inputSize;
    sfx = [sfx '-input' num2str(sz(1)) 'x' num2str(sz(2))];
end
if ~isempty(opts.clusterName)
    sfx = [sfx '-' 'cluster-' opts.clusterName];
end

%
opts.expDir = fullfile('models', ['widerface-' sfx]);
if ~isempty(opts.tag)
    opts.expDir = [opts.expDir '-' opts.tag];
end
[opts, varargin] = vl_argparse(opts, varargin) ;

% imdb
opts.imdbPath = fullfile(opts.dataDir, 'imdb.mat');
[opts, varargin] = vl_argparse(opts, varargin) ;

%

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

%
if opts.testEpoch == 0
    testEpoch = findLastCheckpoint(opts.expDir);
    if testEpoch==0
        error('No available model for evaluation.');
    end
else
    testEpoch = opts.testEpoch;
end
testSet = opts.testSet;

if isempty(opts.testName)
    testName = strrep(opts.expDir, 'models/', '');
    testName = sprintf('%s-epoch%d-%s', testName, testEpoch, testSet);
    testName = [testName '-probthresh' num2str(opts.probThresh)];
    testName = [testName '-nmsthresh' num2str(opts.nmsThresh)];
    if opts.testMultires, testName = [testName '-multires']; end
    if opts.noOrgres, testName = [testName '-noorgres']; end
    if opts.noUpsampling, testName = [testName '-noupsampling']; end
    if opts.noDownsampling, testName = [testName '-nodownsampling']; end
    if opts.evalWithReg, testName = [testName '-evalWithReg']; end
else
    testName = opts.testName;
end

if ~isempty(opts.testTag),
    testName = [testName '-' opts.testTag];
end

if strcmp(opts.resDir, 'results/')
    resDir = fullfile(opts.resDir, testName);
    if numel(resDir) > 255, resDir = strrep(resDir, '-nmsthresh0.3', ''); end
    if numel(resDir) > 255, resDir = strrep(resDir, 'widerface-resnet-50-',''); end
    if numel(resDir) > 255, resDir = strrep(resDir, 'bboxreg-logistic-cluster-',''); end
else
    resDir = opts.resDir;
end

% load ground truth info
%fprintf('Loading ground truth annotation for %s.\n', opts.testSet);
switch opts.testSet
  case 'val'
    test = find(imdb.images.set == 2);
  case 'test'
    error('Evaluting on test set - you must be kidding me...\n')
end

hard_info = load('eval_tools/ground_truth/wider_hard_val.mat');
med_info = load('eval_tools/ground_truth/wider_medium_val.mat');
easy_info = load('eval_tools/ground_truth/wider_easy_val.mat');

% thresholds
minOverlap = opts.minOverlap;
nt = numel(opts.minOverlap);

%
img_idx_map = containers.Map();
for i = 1:numel(hard_info.file_list)
    for j = 1:numel(hard_info.file_list{i})
        img_idx_map(hard_info.file_list{i}{j}) = j;
    end
end

% number of files
easy_npos = 0 ;
med_npos = 0 ;
hard_npos = 0 ;

gt(numel(test)) = struct('BB', [], 'det', []);
easy_gt(numel(test)) = struct('BB', [], 'det', []);
med_gt(numel(test)) = struct('BB', [], 'det', []);
hard_gt(numel(test)) = struct('BB', [], 'det', []);
for i = 1:numel(test)
    eventId = imdb.labels.eventid(test(i));

    gtrects = imdb.labels.rects{test(i)};
    gt(i).BB = gtrects';
    easy_gt(i).det = false(size(gtrects,1), nt);
    med_gt(i).det = false(size(gtrects,1), nt);
    hard_gt(i).det = false(size(gtrects,1), nt);

    gt_h = gtrects(:,4)-gtrects(:,2)+1;
    gt_w = gtrects(:,3)-gtrects(:,1)+1;

    % NOTE: using released evaluation tools
    img_size = imdb.images.size(test(i),:);
    [~, img_name, ~] = fileparts(imdb.images.name{test(i)});
    img_idx = img_idx_map(img_name);

    easy_gt(i).ign = ones(size(gtrects, 1), 1);
    med_gt(i).ign = ones(size(gtrects, 1), 1);
    hard_gt(i).ign = ones(size(gtrects, 1), 1);

    % easy_gt(i).ign(easy_info.ignore_list{eventId}{img_idx}) = 0;
    % med_gt(i).ign(med_info.ignore_list{eventId}{img_idx}) = 0;
    % hard_gt(i).ign(hard_info.ignore_list{eventId}{img_idx}) = 0;
    easy_gt(i).ign(easy_info.gt_list{eventId}{img_idx}) = 0;
    med_gt(i).ign(med_info.gt_list{eventId}{img_idx}) = 0;
    hard_gt(i).ign(hard_info.gt_list{eventId}{img_idx}) = 0;

    easy_npos = easy_npos + sum(easy_gt(i).ign == 0);
    med_npos = med_npos + sum(med_gt(i).ign == 0);
    hard_npos = hard_npos + sum(hard_gt(i).ign == 0);
end

% collect all predictions
detrespath = fullfile(resDir, sprintf('detections.txt'));
if ~exist(detrespath) || 1
    fout = fopen(detrespath, 'w');
    for i = 1:numel(test)
        eventId = imdb.labels.eventid(test(i));

        event = imdb.events.name{eventId};
        [~,imname,~] = fileparts(imdb.images.name{test(i)});
        resfile = fullfile(resDir, event, [imname '.txt']);
        if ~exist(resfile), continue; end
        o = fopen(resfile);
        C = textscan(o,'%d %d %d %d %f','HeaderLines', 2);
        % if this is an older result
        if any(cellfun(@(x)(any(isnan(x))), C))
            fclose(o);
            o = fopen(resfile);
            C = textscan(o,'%d %d %d %d %d %f','HeaderLines', 2);
        end
        x1 = C{1}; y1 = C{2}; w = C{3}; h = C{4};
        x2 = x1 + (w-1); y2 = y1 + (h-1);
        confidence = C{end};

        fclose(o);
        for j = 1:numel(confidence)
            fprintf(fout, '%d %f %d %d %d %d\n',i,confidence(j), ...
                    x1(j),y1(j),x2(j),y2(j));
        end
    end
    fclose(fout);
end

%
[ids,confidence,b1,b2,b3,b4] = textread(detrespath, '%d %f %d %d %d %d');
BB = [b1 b2 b3 b4]';

%
bh = b4 - b2 + 1;
bw = b3 - b1 + 1;

% sort detections by decreasing confidence
[sc,si]=sort(-confidence);
ids=ids(si);
BB=BB(:,si);
%tids = tids(si);

% assign detections to ground truth objects
nd=length(confidence);

easy_tp=zeros(nd,nt);
easy_fp=zeros(nd,nt);
med_tp=zeros(nd,nt);
med_fp=zeros(nd,nt);
hard_tp=zeros(nd,nt);
hard_fp=zeros(nd,nt);

ovmaxs=zeros(nd,1);
tic;
for d=1:nd
    % display progress
    if toc>1
        %fprintf('face: pr: compute: %d/%d\n',d,nd);
        drawnow;
        tic;
   end

    % find ground truth image
    i=ids(d);

    % assign detection to ground truth object if any
    bb=BB(:,d);
    ovmax=-inf;
    jmax = 0 ;

    for j=1:size(gt(i).BB,2)
        bbgt=gt(i).BB(:,j);
        bi=[max(bb(1),bbgt(1)) ; max(bb(2),bbgt(2)) ; min(bb(3),bbgt(3)) ; min(bb(4),bbgt(4))];
        iw=bi(3)-bi(1)+1;
        ih=bi(4)-bi(2)+1;
        if iw>0 & ih>0
            % compute overlap as area of intersection / area of union
            ua=(bb(3)-bb(1)+1)*(bb(4)-bb(2)+1)+...
               (bbgt(3)-bbgt(1)+1)*(bbgt(4)-bbgt(2)+1)-...
               iw*ih;
            ov=iw*ih/ua;
            if ov>ovmax
                ovmax=ov;
                jmax=j;
            end
        end
    end

    % assign detection as true positive/don't care/false positive
    for k=1:nt
        if ovmax>=opts.minOverlap(k)
            %if ~easy_gt(i).ign(j)
            if ~easy_gt(i).ign(jmax)
                if ~easy_gt(i).det(jmax,k) || ~opts.penalizeDuplicate
                    easy_tp(d,k)=1;            % true positive
                    easy_gt(i).det(jmax,k)=true;
                else
                    easy_fp(d,k)=1;            % false positive (multiple detection)
                end
            end
            %if ~med_gt(i).ign(j)
            if ~med_gt(i).ign(jmax)
                if ~med_gt(i).det(jmax,k) || ~opts.penalizeDuplicate
                    med_tp(d,k)=1;            % true positive
                    med_gt(i).det(jmax,k)=true;
                else
                    med_fp(d,k)=1;            % false positive (multiple detection)
                end
            end
            %if ~hard_gt(i).ign(j)
            if ~hard_gt(i).ign(jmax)
                if ~hard_gt(i).det(jmax,k) || ~opts.penalizeDuplicate
                    hard_tp(d,k)=1;            % true positive
                    hard_gt(i).det(jmax,k)=true;
                else
                    hard_fp(d,k)=1;            % false positive (multiple detection)
                end
            end
        else
            easy_fp(d,k)=1;                    % false positive
            med_fp(d,k)=1;                    % false positive
            hard_fp(d,k)=1;                    % false positive
        end
    end
    ovmaxs(d)=ovmax;
end

% compute precision/recall
easy_fp_ = easy_fp;
easy_tp_ = easy_tp;

easy_fp=cumsum(easy_fp,1);
easy_tp=cumsum(easy_tp,1);
easy_rec=easy_tp/easy_npos;
easy_prec=easy_tp./(easy_fp+easy_tp);

easy_aps = zeros(1, nt);
for k = 1:nt
    easy_aps(k)=VOCap(easy_rec(:,k),easy_prec(:,k));
    fprintf('class: face, epoch: %d, subset: %s (easy), minOverlap=%.2f, AP = %.3f\n', ...
            testEpoch, testSet, opts.minOverlap(k), easy_aps(k))
end

med_fp_ = med_fp;
med_tp_ = med_tp;

med_fp=cumsum(med_fp,1);
med_tp=cumsum(med_tp,1);
med_rec=med_tp/med_npos;
med_prec=med_tp./(med_fp+med_tp);

med_aps = zeros(1, nt);
for k = 1:nt
    med_aps(k)=VOCap(med_rec(:,k),med_prec(:,k));
    fprintf('class: face, epoch: %d, subset: %s (medium), minOverlap=%.2f, AP = %.3f\n', ...
            testEpoch, testSet, opts.minOverlap(k), med_aps(k))
end

hard_fp_ = hard_fp;
hard_tp_ = hard_tp;

hard_fp=cumsum(hard_fp,1);
hard_tp=cumsum(hard_tp,1);
hard_rec=hard_tp/hard_npos;
hard_prec=hard_tp./(hard_fp+hard_tp);

hard_aps = zeros(1, nt);
for k = 1:nt
    hard_aps(k)=VOCap(hard_rec(:,k),hard_prec(:,k));
    fprintf('class: face, epoch: %d, subset: %s (hard), minOverlap=%.2f, AP = %.3f\n', ...
            testEpoch, testSet, opts.minOverlap(k), hard_aps(k))
end

if opts.draw
    close all;
    for i = 1:3
        switch i
          case 1
            rec = easy_rec;
            prec = easy_prec;
            aps = easy_aps;
            diffLevel = 'easy';
          case 2
            rec = med_rec;
            prec = med_prec;
            aps = med_aps;
            diffLevel = 'medium';
          case 3
            rec = hard_rec;
            prec = hard_prec;
            aps = hard_aps;
            diffLevel = 'hard';
        end

        figure(i);
        %% plot our precision/recall
        hold on;
        for k = 1:nt
            plot(rec(:,k),prec(:,k),'-','linewidth',3);
            [~,ii] = max(rec(:,k).*prec(:,k));
        end

        hold off;
        grid;
        xlabel 'recall'
        ylabel 'precision'
        %title(sprintf('class: face, subset: %s, AP = %.3f',testSet,ap));
        title(sprintf('class face, epoch %d, subset %s(%s), AP(0.5) %.3f', ...
                      testEpoch, testSet, diffLevel, aps(opts.minOverlap==0.5)));
        xlim([0,1]);
        ylim([0,1]);
        print('-dpng', fullfile(resDir, ['wider_' diffLevel '.png']));
    end
end

% -------------------------------------------------------------------------
function [epoch, iter] = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;

epoch = 0;

tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
% find latest epoch
for i = 1:numel(tokens)
  token = tokens{i}{1};
  ep = str2num(token{1});
  if ep >= epoch
    epoch = ep;
  end
end

author: peiyunh

repo: tiny

save_path: github-repos/MATLAB/peiyunh-tiny


text:
function [updatedStates, rho_i_pj_i_est] = updateStateStruct( currentStates, lmIds, rho_i_pj_i_est, dx)
%UPDATESTATESTRUCT Updates states by applying a step dx

updatedStates = currentStates;
numStates = length(currentStates);
for stIdx = 1:numStates-1
    dr = dx(1+(stIdx-1)*6:3+(stIdx-1)*6);
    phi = dx(4+(stIdx-1)*6:6+(stIdx-1)*6);

    updatedStates{stIdx+1}.r_vi_i = currentStates{stIdx+1}.r_vi_i + dr;
    updatedStates{stIdx+1}.C_vi = Cfrompsi(phi)*currentStates{stIdx+1}.C_vi;
end

initialIdx = (numStates-1)*6 + 1;
for lm_i = 1:length(lmIds)
    idx = initialIdx + (lm_i-1)*3;
    dp = dx(idx:idx+2);
    rho_i_pj_i_est(:, lmIds(lm_i)) = rho_i_pj_i_est(:, lmIds(lm_i)) + dp;
end

end

repo: MSCKF


text:
function niederreiter2_test01 ( )

%*****************************************************************************80
%
%% NIEDERREITER2_TEST01 tests SETFLD2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NIEDERREITER2_TEST01\n' );
  fprintf ( 1, '  SETFLD2 returns the addition, multiplication, and\n' );
  fprintf ( 1, '  subtraction tables for base 2.\n' );

  [ add, mul, sub ] = setfld2 ( 0 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Addition table:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '+  0  1\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '0  %d  %d\n', add(1,1), add(1,2) );
  fprintf ( 1, '1  %d  %d\n', add(2,1), add(2,2) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Multiplication table:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '*  0  1\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '0  %d  %d\n', mul(1,1), mul(1,2) );
  fprintf ( 1, '1  %d  %d\n', mul(2,1), mul(2,2) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Subtraction table:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '-  0  1\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '0  %d  %d\n', sub(1,1), sub(1,2) );
  fprintf ( 1, '1  %d  %d\n', sub(2,1), sub(2,2) );

  return
end


text:
close all
clear

% system setup and initial condition
timestep = 0.1;
system_param.A = [[1 ,timestep, 0]; [0, 1, timestep]; [0, 0, 1]];
system_param.B = [0; 0; timestep];
system_param.ul = -1;
system_param.uu = 1;
system_param.timestep = timestep;
x0 = [-2;0;1];
t0 = 0.0;
total_time = 1.0;

% problem setup
simulator = CBFDT(system_param, x0, t0);

param_mpcdcbf = ParamMPCDCBF(6, 0.1, 10 * eye(3), 10.0 * eye(3), 1.0);
simulator.setOpt('mpcdcbf', param_mpcdcbf);
simulator.sim(total_time);

param_mpcgcbf = ParamMPCGCBF(6, 0.1, 10 * eye(3), 10.0 * eye(3), 1.0);
simulator.setOpt('mpcgcbf', param_mpcgcbf);
simulator.sim(total_time);

param_nmpcdcbf = ParamNMPCDCBF(6, 6, 0.1, 10 * eye(3), 10.0 * eye(3), 1.0, 10.0);
simulator.setOpt('nmpcdcbf', param_nmpcdcbf);
simulator.sim(total_time);

param_dclfdcbf = ParamDCLFDCBF(0.1, 0.1, 10 * eye(3), 1.0, 10.0);
simulator.setOpt('dclfdcbf', param_dclfdcbf);
simulator.sim(total_time);

param_nmpcdclfdcbf = ParamNMPCDCLFDCBF(6, 6, 6, 0.1, 0.1, 10.0 * eye(3), 10.0 * eye(3), 1.0, 10.0, 10.0);
simulator.setOpt('nmpcdclfdcbf', param_nmpcdclfdcbf);
simulator.sim(total_time);

author: HybridRobotics

repo: NMPC-DCLF-DCBF


text:
function [transformationMatrix] = ...
    phase_difference_linear_registration3d(moving, movingCertainty, ...
    fixed, fixedCertainty, varargin)
% PHASE_DIFFERENCE_LINEAR_REGISTRATION3D Estimates a transformation matrix using phase-difference
%
% [transformationMatrix] = ...
%     phase_difference_linear_registration3d(moving, movingCertainty, ...
%     fixed, fixedCertainty)
%
% INPUT ARGUMENTS
% moving                    - Moving image
% movingCertainty           - Certainty mask of moving image
% fixed                     - Fixed image
% fixedCertainty            - Certainty mask of fixed image
%
% OPTIONAL INPUT ARGUMENTS
% 'transformationModel'     - Transformation model for estimating the
%                             displacement field
%                             'translation', 'affine' (default)
%
% OUTPUT ARGUMENTS
% transformationMatrix      - Estimated transformation matrix

% Copyright (c) 2012 Daniel Forsberg
% danne.forsberg@outlook.com
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

%% Setup default parameters
% translation, affine, non-rigid
transformationModel = 'affine';

% Overwrites default parameter
for k=1:2:length(varargin)
    eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

% Initialize transformation matrix
transformationMatrix = eye(4);

% Load quadrature filters
load quadratureFiltersLinearRegistration3D

% Perform quadrature filtering
q21 = imfilter(moving,f1,'same','conv');
q22 = imfilter(moving,f2,'same','conv');
q23 = imfilter(moving,f3,'same','conv');

q11 = imfilter(fixed,f1,'same','conv');
q12 = imfilter(fixed,f2,'same','conv');
q13 = imfilter(fixed,f3,'same','conv');

filterSize = (size(real(f1),1)-1)/2;
q11 = q11(filterSize+1:end-filterSize,filterSize+1:end-filterSize,filterSize+1:end-filterSize);
q12 = q12(filterSize+1:end-filterSize,filterSize+1:end-filterSize,filterSize+1:end-filterSize);
q13 = q13(filterSize+1:end-filterSize,filterSize+1:end-filterSize,filterSize+1:end-filterSize);
q21 = q21(filterSize+1:end-filterSize,filterSize+1:end-filterSize,filterSize+1:end-filterSize);
q22 = q22(filterSize+1:end-filterSize,filterSize+1:end-filterSize,filterSize+1:end-filterSize);
q23 = q23(filterSize+1:end-filterSize,filterSize+1:end-filterSize,filterSize+1:end-filterSize);

% Compute phase-difference, certainties and phase-gradients
dphi(:,:,:,1) = angle(q11.*conj(q21));
dphi(:,:,:,2) = angle(q12.*conj(q22));
dphi(:,:,:,3) = angle(q13.*conj(q23));

% Estimate certainties
certainty(:,:,:,1) = abs(q11.*q21).*((cos(dphi(:,:,:,1)/2)).^2);
certainty(:,:,:,2) = abs(q12.*q22).*((cos(dphi(:,:,:,2)/2)).^2);
certainty(:,:,:,3) = abs(q13.*q23).*((cos(dphi(:,:,:,3)/2)).^2);

% Add certainty masks
if ~isempty(fixedCertainty)
    fixedCertainty = fixedCertainty(...
        filterSize+1:end-filterSize,filterSize+1:end-filterSize,filterSize+1:end-filterSize);
    certainty = bsxfun(@times,certainty,fixedCertainty);
end
if ~isempty(movingCertainty)
    movingCertainty = movingCertainty(...
        filterSize+1:end-filterSize,filterSize+1:end-filterSize,filterSize+1:end-filterSize);
    certainty = bsxfun(@times,certainty,movingCertainty);
end

% Estimate gradients of phi
grad_x_dphi_n1 = zeros(size(q11));
grad_x_dphi_n1(:,2:end-1,:) = ...
    angle(q11(:,3:end,:).*conj(q11(:,2:end-1,:)) + ...
    q11(:,2:end-1,:).*conj(q11(:,1:end-2,:)) + ...
    q21(:,3:end,:).*conj(q21(:,2:end-1,:)) + ...
    q21(:,2:end-1,:).*conj(q21(:,1:end-2,:)));

grad_y_dphi_n2 = zeros(size(q11));
grad_y_dphi_n2(2:end-1,:,:) = ...
    angle(q12(3:end,:,:).*conj(q12(2:end-1,:,:)) + ...
    q12(2:end-1,:,:).*conj(q12(1:end-2,:,:)) + ...
    q22(3:end,:,:).*conj(q22(2:end-1,:,:)) + ...
    q22(2:end-1,:,:).*conj(q22(1:end-2,:,:)));

grad_z_dphi_n3 = zeros(size(q11));
grad_z_dphi_n3(:,:,2:end-1) = ...
    angle(q13(:,:,3:end).*conj(q13(:,:,2:end-1)) + ...
    q13(:,:,2:end-1).*conj(q13(:,:,1:end-2)) + ...
    q23(:,:,3:end).*conj(q23(:,:,2:end-1)) + ...
    q23(:,:,2:end-1).*conj(q23(:,:,1:end-2)));

% Save all the phase gradients nicely...
phaseGradient(:,:,:,1) = grad_x_dphi_n1;
phaseGradient(:,:,:,2) = grad_y_dphi_n2;
phaseGradient(:,:,:,3) = grad_z_dphi_n3;

% Build A and h for A p = h
[A, h] = build_A_h_3d(dphi, certainty, phaseGradient, transformationModel);

% Solve the equation system
switch transformationModel
    case 'translation'
        d = A \ h;
        transformationMatrix(1:3,4) = d;
    case {'rigid','affine'}
        p = A \ h;
        transformationMatrix(1:3,1:4) = [1+p(1) p(2) p(3) p(10);...
            p(4) 1+p(5) p(6) p(11);...
            p(7) p(8) 1+p(9) p(12)];
end

author: fordanic

repo: image-registration


text:
function Z = replace(X,Y,W)
%REPLACE Substitutes variables
%
%Z = REPLACE(Y,X,W)  Replaces any occurence of variable Y in X with W

Z = ncvar_replace(X,Y,W);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function material_cell = mean_material_node(grid3d, gt, material_node)

chkarg(istypesizeof(grid3d, 'Grid3d'), '"grid3d" should be instance of Grid.');
chkarg(istypesizeof(gt, 'GT'), '"gt" should be instance of GT.');
chkarg(istypesizeof(material_node, 'complexcell', [1 Axis.count], grid3d.N), ...
	'"material_node" should be length-%d cell array, whose each element is %d-by-%d-by-%d array with complex elements.', ...
	Axis.count, grid3d.Ncell{:});

for w = Axis.elems
	material_node{w} = expand_node_array(grid3d, material_node{w});  % (Nx+2) x (Ny+2) x (Nz+2)
end

if gt == GT.prim
	% material parameters for fields on primary grid
	material_cell = arithmetic_mean_material_node(material_node);
else  % gt == GT.dual
	% material parameters for fields on dual grid
	material_cell = harmonic_mean_material_node(material_node);
end

function material_edge_cell = arithmetic_mean_material_node(material_node)
material_edge_cell = cell(1, Axis.count);

material_edge_cell{Axis.x} = (...
	material_node{Axis.x}(2:end-1, 2:end-1, 2:end-1) ...
	+ material_node{Axis.x}(2:end-1, 1:end-2, 2:end-1) ...
	+ material_node{Axis.x}(2:end-1, 2:end-1, 1:end-2) ...
	+ material_node{Axis.x}(2:end-1, 1:end-2, 1:end-2)...
) / 4;

material_edge_cell{Axis.y} = (...
	material_node{Axis.y}(2:end-1, 2:end-1, 2:end-1) ...
	+ material_node{Axis.y}(2:end-1, 2:end-1, 1:end-2) ...
	+ material_node{Axis.y}(1:end-2, 2:end-1, 2:end-1) ...
	+ material_node{Axis.y}(1:end-2, 2:end-1, 1:end-2)...
) / 4;

material_edge_cell{Axis.z} = (...
	material_node{Axis.z}(2:end-1, 2:end-1, 2:end-1) ...
	+ material_node{Axis.z}(1:end-2, 2:end-1, 2:end-1) ...
	+ material_node{Axis.z}(2:end-1, 1:end-2, 2:end-1) ...
	+ material_node{Axis.z}(1:end-2, 1:end-2, 2:end-1)...
) / 4;

function material_face_cell = harmonic_mean_material_node(material_node)
material_face_cell = cell(1, Axis.count);
material_face_cell{Axis.x} = 2./(1./material_node{Axis.x}(1:end-2, 2:end-1, 2:end-1) + 1./material_node{Axis.x}(2:end-1, 2:end-1, 2:end-1));
material_face_cell{Axis.y} = 2./(1./material_node{Axis.y}(2:end-1, 1:end-2, 2:end-1) + 1./material_node{Axis.y}(2:end-1, 2:end-1, 2:end-1));
material_face_cell{Axis.z} = 2./(1./material_node{Axis.z}(2:end-1, 2:end-1, 1:end-2) + 1./material_node{Axis.z}(2:end-1, 2:end-1, 2:end-1));

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function [im, im_scale] = prep_im_for_blob(im, im_means, target_size, max_size)
    im = single(im);
    if size(im,3) < 3
      im = im(:,:,[1 1 1]);
    end
    if ~isa(im, 'gpuArray')
        try
            im = bsxfun(@minus, im, im_means);
        catch
            im_means = imresize(im_means, [size(im, 1), size(im, 2)], 'bilinear', 'antialiasing', false);
            if  size(im,3)>3 && size(im,3) ~= size(im_means,3)
              im(:,:,1:3,:) = bsxfun(@minus, im(:,:,1:3,:), im_means);
            else
              im = bsxfun(@minus, im, im_means);
            end
        end
        im_scale = prep_im_for_blob_size(size(im), target_size, max_size);

        target_size = round([size(im, 1), size(im, 2)] * im_scale);
        im = imresize(im, target_size, 'bilinear', 'antialiasing', false);
    else
        % for im as gpuArray
        try
            im = bsxfun(@minus, im, im_means);
        catch
            im_means_scale = max(double(size(im, 1)) / size(im_means, 1), double(size(im, 2)) / size(im_means, 2));
            im_means = imresize(im_means, im_means_scale);
            y_start = floor((size(im_means, 1) - size(im, 1)) / 2) + 1;
            x_start = floor((size(im_means, 2) - size(im, 2)) / 2) + 1;
            im_means = im_means(y_start:(y_start+size(im, 1)-1), x_start:(x_start+size(im, 2)-1));
            im = bsxfun(@minus, im, im_means);
        end

        im_scale = prep_im_for_blob_size(size(im), target_size, max_size);
        im = imresize(im, im_scale);
    end
end

author: feichtenhofer

repo: Detect-Track


text:
M = 5;
S1 = 4;
S2 = 4;
S = S1 + S2;
K = 2;
Y = zeros(M, S);
Y(1, 1:4) = [1 2 3 4];
Y(2, 1:4) = [3 3 2 6];
Y(3, 5:8) = [3 3 6 6];
Y(4, 5:8) = [4 5 6 7];
Y
Y = spx.norm.normalize_l2(Y)
C = spx.fast.batch_flipped_omp_spr(Y, K, 0)

Y * C
abs(Y - Y * C)

author: indigits

repo: sparse-plex


text:
function [mneObj, DataMat, ChannelMat, iChannels] = out_mne_data(DataFiles, ObjType, ChannelFile, SensorTypes)
% OUT_MNE_DATA: Converts a data file or a set of data files into a MNE-Python object (Raw, Epoched or Evoked)
%
% USAGE:  [mneObj, DataMat, ChannelMat, iChannels] = out_mne_data( DataFiles, ObjType, ChannelFile=[], SensorTypes/iChannels=[]);
%         [mneObj, DataMat, ChannelMat, iChannels] = out_mne_data( DataMat,   ObjType, ChannelMat=[],  SensorTypes/iChannels=[]);
%
% INPUTS:
%    - DataFiles    : String or cell-array of strings, relative path(s) to data file(s) available in the database
%    - DataMat      : Brainstorm data file structure (if multiple files, the F matrix matrix must be [nEpochs, nChans, nTime])
%    - ObjType      : Class of the returned MNE-Python object {'Raw', 'Epoched', 'Evoked'}
%    - ChannelFile  : Relative path to a channel file available in the database (if not provided: look for it based on the DataFile)
%    - ChannelMat   : Brainstorm channel file structure
%    - iChannels    : Vector of selected channel indices
%    - SensorTypes  : Names or types of channels, separated with commas

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

%% ===== PARSE INPUTS =====
if (nargin < 4) || isempty(SensorTypes)
    SensorTypes = [];
end
if (nargin < 3) || isempty(ChannelFile)
    ChannelMat  = [];
    ChannelFile = [];
elseif isstruct(ChannelFile)
    ChannelMat   = ChannelFile;
    ChannelFile = [];
else
    ChannelMat = [];
end
if isempty(DataFiles)
    error('No input data files.');
elseif ischar(DataFiles)
    DataFiles = {DataFiles};
    DataMat = [];
elseif isstruct(DataFiles)
    DataMat = DataFiles;
    DataFiles = [];
end
% Only 'Epoched' objects can have multiple data files in input
if ~ismember(ObjType, {'Raw', 'Epoched', 'Evoked'})
    error('ObjType must be one of the following types of data: ''Raw'', ''Epoched'', ''Evoked''');
elseif (length(DataFiles) > 1) && ~strcmpi(ObjType, 'Epoched')
    error('Only "Epoched" objects accept multiple input files.');
end
% Check that data files are available in the database
MeasDate = [];
if ~isempty(DataFiles) && isappdata(0, 'BrainstormRunning')
    sStudy = bst_get('DataFile', DataFiles{1});
    % Get study date
    if isempty(sStudy)
        warning(['File not found in database: ' DataFiles{1}]);
    else
        MeasDate = sStudy.DateOfStudy;
    end
end

%% ===== LOAD CHANNEL FILE =====
% Get ChannelFile if not provided
if isempty(ChannelFile) && isempty(ChannelMat)
    if ~isempty(DataFiles)
        ChannelFile = bst_get('ChannelFileForStudy', DataFiles{1});
    else
        error('Missing ChannelMat or ChannelFile in input.');
    end
end
% Load channel file
if ~isempty(ChannelFile) && isempty(ChannelMat)
    ChannelMat = in_bst_channel(ChannelFile);
end
% Make sure that the channel file is defined
if isempty(ChannelMat)
    error('No channel file available for the input files.');
end
% Find sensors by names/types
if ~isempty(SensorTypes)
    if ischar(SensorTypes)
        iChannels = channel_find(ChannelMat.Channel, SensorTypes);
        if isempty(iChannels)
            error(['Channels not found: ' SensorTypes]);
        end
    elseif isnumeric(SensorTypes)
        iChannels = SensorTypes;
    else
        error('Invalid input type for parameter "SensorTypes".');
    end
% Default channel selection: all
else
    iChannels = 1:length(ChannelMat.Channel);
end

%% ===== LOAD DATA =====
% Load data file
if ~isempty(DataFiles)
    epochsComment = cell(1,length(DataFiles));
    % Loop on input files
    for iFile = 1:length(DataFiles)
        % Load .mat files
        tmpData = in_bst_data(DataFiles{iFile});
        % Raw files: Read entire recordings
        if isstruct(tmpData.F)
            sFile = tmpData.F;
            RemoveBaseline = 'no';
            UseSsp = 1;
            [tmpData.F, tmpData.Time] = panel_record('ReadRawBlock', sFile, ChannelMat, 1, tmpData.Time([1 end]), 0, 1, RemoveBaseline, UseSsp, iChannels);
            tmpData.Events = sFile.events;
            % Get date
            if isfield(sFile, 'acq_date') && ~isempty(sFile.acq_date)
                MeasDate = sFile.acq_date;
            end
        % Imported data: Keep only selected
        else
            tmpData.F = tmpData.F(iChannels,:);
        end
        tmpData.ChannelFlag = tmpData.ChannelFlag(iChannels);
        % Concatenate all data in the same structure (F = [nEpoch x nChan x nTime])
        if strcmpi(ObjType, 'Epoched')
            tmpData.F = reshape(tmpData.F, 1, size(tmpData.F,1), size(tmpData.F,2));
            if (iFile == 1)
                DataMat = tmpData;
            else
                if ~isequal(size(DataMat.F(1,:,:)), size(tmpData.F))
                    error('All input data must have the same dimensions.');
                end
                DataMat.F = cat(1, DataMat.F, tmpData.F);
                DataMat.ChannelFlag(tmpData.ChannelFlag ~= 1) = -1;
            end
            epochsComment{iFile} = str_remove_parenth(tmpData.Comment, '(');
            epochsComment{iFile} = str_remove_parenth(epochsComment{iFile}, '[');
        else
            DataMat = tmpData;
        end
    end
% Data structures passed in input
else
    % No epochs classification
    epochsComment = [];
    % Select data channels
    DataMat.F = DataMat.F(iChannels,:);
    DataMat.ChannelFlag = DataMat.ChannelFlag(iChannels);
end

%% ===== CREATE INFO OBJECT =====
% Create info object
mneInfo = out_mne_channel(ChannelFile, iChannels);
% Unlocking Info object
mneSetStatus = py.getattr(mneInfo, '__setstate__');
mneSetStatus(py.dict(pyargs('_unlocked', true)));
% Sampling frequency
mneInfo{'sfreq'} = 1 ./ (DataMat.Time(2) - DataMat.Time(1));
% Description
mneInfo{'description'} = DataMat.Comment;
% Bad channels
iBad = find(DataMat.ChannelFlag == -1);
if ~isempty(iBad)
    mneInfo{'bads'} = py.list({ChannelMat.Channel(iBad).Name});
end
% Mark projectors as applied (data loaded with UseSSP=1)
for iProj = 1:length(mneInfo{'projs'})
    mneInfo{'projs'}{iProj}{'active'} = py.bool(true);
end
% Measurement date
if ~isempty(MeasDate)
    try
        % Read date string to a py.datetime object
        dt = py.dateutil.parser.parse(MeasDate);
        % Convert to UTC
        tz = py.datetime.datetime.now().astimezone().tzinfo;
        dt = dt.replace(pyargs('tzinfo', tz));
        dt = dt.astimezone(py.datetime.timezone(py.datetime.timedelta(0)));
        mneInfo{'meas_date'} = dt;
    catch
    end
end
% Locking Info object again
mneSetStatus(py.dict(pyargs('_unlocked', false)));

% Object: Raw
switch ObjType
    % RAW:  https://www.nmr.mgh.harvard.edu/mne/stable/generated/mne.io.Raw.html
    case 'Raw'
        % Create Raw object
        first_samp = round(DataMat.Time(1) .* mneInfo{'sfreq'});
        mneObj = py.mne.io.RawArray(pyargs(...
            'data',       bst_mat2py(DataMat.F), ...
            'info',       mneInfo, ...
            'first_samp', first_samp));

        % Add events
        for iEvt = 1:length(DataMat.Events)
            % No occurrences: skip
            if isempty(DataMat.Events(iEvt).times)
                continue;
            end
            % Get onsets
            annotOnset = DataMat.Events(iEvt).times(1,:);
            % Extended events / simple events
            if (size(DataMat.Events(iEvt).times,1) == 2)
                annotDuration = DataMat.Events(iEvt).times(2,:) - DataMat.Events(iEvt).times(1,:);
            else
                annotDuration = 0 .* annotOnset;
            end
            % Add annotations to MNE objet
            mneObj.annotations.append(annotOnset, annotDuration, repmat({DataMat.Events(iEvt).label}, 1, size(annotOnset,2)));
        end

    case 'Epoched'
        % Sort trials by type, based on the comment of the files
        evts = uint32([(1:size(DataMat.F, 1))', repmat([0, 1], size(DataMat.F, 1), 1)]);
        event_id = py.dict();
        if ~isempty(epochsComment)
            uniqueTypes = unique(epochsComment);
            for iType = 1:length(uniqueTypes)
                evts(strcmpi(epochsComment, uniqueTypes{iType}), 3) = iType;
                event_id{uniqueTypes{iType}} = uint32(iType);
            end
        end
        % Create Epoched object from concatenated trials
        mneObj = py.mne.EpochsArray(pyargs(...
            'data',     DataMat.F, ...
            'info',     mneInfo, ...
            'events',   bst_mat2py(evts, 1), ...
            'tmin',     DataMat.Time(1), ...
            'event_id', event_id));

    case 'Evoked'
        % Create Evoked object
        mneObj = py.mne.EvokedArray(pyargs(...
            'data',    DataMat.F, ...
            'info',    mneInfo, ...
            'tmin',    DataMat.Time(1), ...
            'comment', DataMat.Comment, ...
            'nave',    uint32(DataMat.nAvg)));
end

author: brainstorm-tools


text:
function OUT = multiscaleRetinex(I, method, varargin)

if nargin == 0
%     I = imload; %imread(imgetfile); %imread('office_1.jpg');
%     J = multiscaleRetinex(I, 'MSRCR');
%     ezFig I J
%     return;

    I = imload; %imread(imgetfile); %('office_1.jpg');
    method = Popupmenu({'MSRCR','MSR','SSR'});
    J = ImCtrl(@multiscaleRetinex, I, method);
    ezFig I J
end

if nargin == 1
   method = 'MSRCR';
end

I = im2double(I).*255+1; % +1 to avoid Inf
SSR = @SSRlog; % @SSR; %

% TODO: gray or RGB

switch (method)
    case {'SSR' 'single scale retinex'}, f = SSR;
    case {'MSR' 'multi scale retinex'}, f = @MSR;
    case 'MSRCR', f = @MSRCR;
    %case 'MSRCP', f = @MSRCP;
    otherwise
        return;
        %f = str2func(method);
end

OUT = f(I, varargin{:});

end

function OUT = SSR(I, varargin)
T = imgaussfilt(I, varargin{:});
OUT = I./(T); % avoid NaN
end

function OUT = SSRlog(I, varargin)
T = imgaussfilt(I, varargin{:});
OUT = log(I) - log(T+1) + 0.5;
end

function OUT = MSR(I, varargin)
if numel(varargin) == 0
    varargin = {25 100 240};
end
OUT = 0; N = numel(varargin);
for n = 1:N
    OUT = OUT + (1/N)* multiscaleRetinex(I,'SSR',varargin{n});
end
end

function OUT = MSRCR(I, lowScale, medScale, highScale, leftChop, rightChop)
if ~exist('lowScale', 'var'), lowScale = 15; end
if ~exist('MedScale', 'var'), medScale = 80; end
if ~exist('HighScale', 'var'), highScale = 250; end
if ~exist('s1', 'var'), leftChop = 1; end
if ~exist('s2', 'var'), rightChop = 1; end

MSR = multiscaleRetinex(I, 'MSR', lowScale, medScale, highScale);

for c = 1:3
    CR = (log(125*I(:,:,c))-log(I(:,:,1)+I(:,:,2)+I(:,:,3)));
    OUT(:,:,c) = colorBalance(mat2gray(CR.*MSR(:,:,c)), 'simplest', leftChop, rightChop);
end
%OUT = max(0, min(1, OUT));
end

author: dawnlh

repo: awesome-low-light-image-enhancement


text:
function [g,filtertype] = gsp_jtv_design_meyer(G,Nf)
%GSP_DESIGN_JTV_MEYER Design the jtv Meyer tight filterbank
%   Usage: [g,filtertype] = gsp_design_jtv_meyer(G);
%          [g,filtertype] = gsp_design_jtv_meyer(G,Nf);
%
%   Input parameters:
%       G       : Time-Vertex graph structure
%       Nf      : Number of filters for each domain (total number Nf^2, default Nf = 4)
%   Output parameters:
%       g          : Cell array of time-vertex filters
%       filtertype : Filter domain js
%

% Author :  Francesco Grassi
% Date   : September 2016

if nargin<2
    Nf = 4;
end

%graph meyer
g1 = gsp_design_meyer(G,Nf);

%time meyer
g2 = gsp_design_meyer(0.5/G.jtv.fs,Nf);

%building jtv meyer separable filterbank
n = 0;
g = cell(Nf^2,1);
for ii=1:Nf
    for jj=1:Nf
        n = n+1;
        g{n} = @(x,y) g1{ii}(x).*g2{jj}(abs(y));
    end
end

filtertype = 'js';

end


text:
function varargout = root(F)
%ROOT   Find one common zero of a CHEBFUN3V object.
%   R = ROOT(F) finds one common zero of the three CHEBFUN3 objects F(1),
%   F(2) and F(3) in their domain of definition under the assumption that
%   the solution set is zero-dimensional. R is a row vector storing the
%   x-value, y-value, and z-value of the common root.
%   This function is also called by the syntax ROOT(F, G, H), where F, G and
%   H are CHEBFUN3 objects.
%
%   [x, y, z] = ROOT(F) returns the x-value, y-value and z-value as
%   three separate entries.
%
% See also CHEBFUN3/ROOT and CHEBFUN2/ROOTS.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Empty check:
if ( isempty(F) )
    varargout = {[]};
    return
end

f = F.components{1};
g = F.components{2};
h = F.components{3};

r = root(f, g, h);

if ( nargout <= 1 )
    varargout{1} = r;
else
    varargout = {r(1), r(2), r(3)};
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function EarTraining()

%Richard Moore
%November 18, 2012

%This program is intended to test/train for perfect pitch.  Currently my
%plan is to:
%-Start with just a couple of pitches, probably C and G.
%-Play the tones in different octaves and have the observer guess the
%pitch.
%-As proficiency is demonstrated within a moving average, add the next note
%in the circle of fifths.
%-Also, when a mistake is made, notify the user that they have made a
%mistake and let them know what the correct note name is.

%C, G, D, A, E, B, F#, C#, G#, Eb, Bb, F
NoteNames = {'C'; 'G'; 'D'; 'A'; 'E'; 'B'; 'F#'; 'C#'; 'G#'; 'Eb'; 'Bb'; 'F'};
%Row = presented
%Column = responded
ResponseMat = zeros(12,12);
A = 440;
halfStep = 2^(1/12);
baseFreqs = [A*halfStep^3; A*halfStep^(-2); A*halfStep^5; A; A*halfStep^7; A*halfStep^2; A*halfStep^9; A*halfStep^4; A*halfStep^(-1); A*halfStep^6; A*halfStep^1; A*halfStep^8];
%numPitches = 8;
AccThresh = 0.9;
MinNumTrials = 20;

disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
disp('EAR TRAINING FOR PITCH RECOGNITION');
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
disp(' ');
disp('You will hear notes and be asked to identify them.');
disp(['Once you have correctly identified at least ' num2str(AccThresh*100) '% of ' num2str(MinNumTrials) ' trials, another note will be added.']);
numPitches = input('How many pitches would you like to start out with? ');
disp(['Very weel, then.  You will be starting with ' num2str(numPitches) ' different notes.']);
disp([NoteNames{1:numPitches}]);
disp(' ');
pause(3);
disp('May the force be with you.');
pause(1);
clc;
%So, if the probability exceeds a threshold, add another pitch.
trialNum = 0;
testFlag = 1;
sampRate = 44100;
noteLength = 0.5; %In seconds
while testFlag
    %
    trialNum = trialNum + 1;
    %Present the pitch:
    oct = ceil(rand*4)-2;
    currPitch = ceil(rand*numPitches);
    soundsc(sin(2*pi*(2^oct)*baseFreqs(currPitch)*[0:(1/sampRate):noteLength]),sampRate);
    %Take the response:
    resp = input('Name that note: ','s');
    %This won't work quite yet.  The strings are not all the same length.
    switch resp
        case {'c','C'}
            resp2 = 1;
        case {'g','G'}
            resp2 = 2;
        case {'d','D'}
            resp2 = 3;
        case {'a','A'}
            resp2 = 4;
        case {'e','E'}
            resp2 = 5;
        case {'b','B'}
            resp2 = 6;
        case {'f#','F#','gb','Gb'}
            resp2 = 7;
        case {'c#','C#','db','Db'}
            resp2 = 8;
        case {'g#','G#','ab','Ab'}
            resp2 = 9;
        case {'eb','Eb','d#','D#'}
            resp2 = 10;
        case {'bb','Bb','a#','A#'}
            resp2 = 11;
        case {'f','F'}
            resp2 = 12;
        otherwise
            disp('Ummm, I think that you just made up a new note.');
    end
    %Populate the confusion matrix:
    ResponseMat(currPitch,resp2) = ResponseMat(currPitch,resp2) + 1;
    %Notify the user of a mistake, if applicable:
    if currPitch~=resp2
        disp(['Whoopsie-daisies!  That was actually ' NoteNames{currPitch}]);
    end
    %If skill level is mastered, increase numPitches and notify user:
    currNum = sum(sum(ResponseMat.*eye(size(ResponseMat,1))));
    currDen = sum(ResponseMat(:));
    %Require a minimum score of 90% and a minimum number of 10 trials
    if (currDen>=MinNumTrials)&&((currNum/currDen)>=AccThresh)
        numPitches = numPitches + 1;
        %disp(ResponseMat);
        clc;
        disp('%%%%%%%%%%%%%%%%%%%%');
        disp('Congrats!  You just leveled up!');
        disp(['You now have to identify ' num2str(numPitches) ' different notes: ' NoteNames{1:numPitches}]);
        disp('%%%%%%%%%%%%%%%%%%%%');
        ResponseMat = zeros(12,12);
    end
    if rem(trialNum,10)==0
        %disp(ResponseMat);
        disp(['Since the last level-up your accuracy has been ' num2str(100*currNum/currDen) '%.']);
        testFlag = input('Wanna keep testing? 1=yes, 0=no: ');
    end
end

clc;
disp('Sure, go do something else.');
disp(' ');
pause(1);
disp('I''m sure that Netflix is more important than becoming a ear-training-jedi-master.');
pause(3);

author: Sable

repo: mcbench-benchmarks


text:
function [X_p] = gmm_project(X, mapping)
[N,B] = size(X);
X_p = (X - repmat(mapping.mean,N,1)) * mapping.M;
end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function [Phi,dim1,dim2] = spm_def2sparse(PY,PI)
% Generate a sparse matrix encoding a deformation
% [Phi,dim1,dim2] = spm_def2sparse(PY,PI)
% PY - Filename of deformation field
% PI - Filename of image defining field of view etc
%_______________________________________________________________________
% Copyright (C) Wellcome Trust Centre for Neuroimaging (2009)

% John Ashburner
% $Id: spm_def2sparse.m 4861 2012-08-24 15:56:39Z john $

NY=nifti(PY);
NI=nifti(PI);

dim1 = [size(NY.dat) 1 1];
dim1 = dim1(1:3);

dim2 = [size(NI.dat) 1 1];
dim2 = dim2(1:3);

X1=NY.dat(:,:,:,1,1);
X2=NY.dat(:,:,:,1,2);
X3=NY.dat(:,:,:,1,3);

M = inv(NI.mat);

Y1=M(1,1)*X1+M(1,2)*X2+M(1,3)*X3+M(1,4);
Y2=M(2,1)*X1+M(2,2)*X2+M(2,3)*X3+M(2,4);
Y3=M(3,1)*X1+M(3,2)*X2+M(3,3)*X3+M(3,4);

if false
    % Nearest Neighbour
    fY1 = (round(Y1));
    fY2 = (round(Y2));
    fY3 = (round(Y3));
    I   = find(fY1>=1 & fY1<=dim2(1) & fY2>=1 & fY2<=dim2(2) & fY3>=1 & fY3<=dim2(3));
    J   = fY1(I) + dim2(1)*(fY2(I)-1 + dim2(2)*(fY3(I)-1));
    Phi = sparse(I,J,1,prod(dim1),prod(dim2));
else
    % Trilinear
    fY1 = (floor(Y1));
    fY2 = (floor(Y2));
    fY3 = (floor(Y3));

    I   = find(fY1>=1 & fY1<=dim2(1) & fY2>=1 & fY2<=dim2(2) & fY3>=1 & fY3<=dim2(3));
    J   = fY1(I) + dim2(1)*(fY2(I)-1 + dim2(2)*(fY3(I)-1));
    S   = (1-Y1+fY1).*(1-Y2+fY2).*(1-Y3+fY3);
    S   = S(I);
    Phi = sparse(I,J,S,prod(dim1),prod(dim2));

    I   = find(fY1>=0 & fY1<=dim2(1)-1 & fY2>=1 & fY2<=dim2(2) & fY3>=1 & fY3<=dim2(3));
    J   = fY1(I)+1 + dim2(1)*(fY2(I)-1 + dim2(2)*(fY3(I)-1));
    S   = (Y1-fY1).*(1-Y2+fY2).*(1-Y3+fY3);
    S   = S(I);
    Phi = Phi + sparse(I,J,S,prod(dim1),prod(dim2));

    I   = find(fY1>=1 & fY1<=dim2(1) & fY2>=0 & fY2<=dim2(2)-1 & fY3>=1 & fY3<=dim2(3));
    J   = fY1(I) + dim2(1)*(fY2(I) + dim2(2)*(fY3(I)-1));
    S   = (1-Y1+fY1).*(Y2-fY2).*(1-Y3+fY3);
    S   = S(I);
    Phi = Phi + sparse(I,J,S,prod(dim1),prod(dim2));

    I   = find(fY1>=0 & fY1<=dim2(1)-1 & fY2>=0 & fY2<=dim2(2)-1 & fY3>=1 & fY3<=dim2(3));
    J   = fY1(I)+1 + dim2(1)*(fY2(I) + dim2(2)*(fY3(I)-1));
    S   = (Y1-fY1).*(Y2-fY2).*(1-Y3+fY3);
    S   = S(I);
    Phi = Phi + sparse(I,J,S,prod(dim1),prod(dim2));

    I   = find(fY1>=1 & fY1<=dim2(1) & fY2>=1 & fY2<=dim2(2) & fY3>=0 & fY3<=dim2(3)-1);
    J   = fY1(I) + dim2(1)*(fY2(I)-1 + dim2(2)*fY3(I));
    S   = (1-Y1+fY1).*(1-Y2+fY2).*(Y3-fY3);
    S   = S(I);
    Phi = Phi + sparse(I,J,S,prod(dim1),prod(dim2));

    I   = find(fY1>=0 & fY1<=dim2(1)-1 & fY2>=1 & fY2<=dim2(2) & fY3>=0 & fY3<=dim2(3)-1);
    J   = fY1(I)+1 + dim2(1)*(fY2(I)-1 + dim2(2)*fY3(I));
    S   = (Y1-fY1).*(1-Y2+fY2).*(Y3-fY3);
    S   = S(I);
    Phi = Phi + sparse(I,J,S,prod(dim1),prod(dim2));

    I   = find(fY1>=1 & fY1<=dim2(1) & fY2>=0 & fY2<=dim2(2)-1 & fY3>=0 & fY3<=dim2(3)-1);
    J   = fY1(I) + dim2(1)*(fY2(I) + dim2(2)*fY3(I));
    S   = (1-Y1+fY1).*(Y2-fY2).*(Y3-fY3);
    S   = S(I);
    Phi = Phi + sparse(I,J,S,prod(dim1),prod(dim2));

    I   = find(fY1>=0 & fY1<=dim2(1)-1 & fY2>=0 & fY2<=dim2(2)-1 & fY3>=0 & fY3<=dim2(3)-1);
    J   = fY1(I)+1 + dim2(1)*(fY2(I) + dim2(2)*fY3(I));
    S   = (Y1-fY1).*(Y2-fY2).*(Y3-fY3);
    S   = S(I);
    Phi = Phi + sparse(I,J,S,prod(dim1),prod(dim2));
end

author: spm


text:
function CrossCutter(img, cb_name, cb_format)
%
%        CrossCutter(img, cb_name, cb_format)
%
%
%        Input:
%           -img: input environment map encoded as a cubemap
%           -name: a string representing the prefix name for each face of
%            the cubemap. For example: 'output_cubemap' (default)
%           -format: the output format of each face of the cubemap. For
%           example: 'hdr' (default)
%        Output:
%           -ret: it is set to 1 if the function succeeded
%
%     Copyright (C) 2011  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

if(~exist('cb_name', 'var'))
    cb_name = 'output_cube_map';
end

if(~exist('cb_format', 'var'))
    cb_name = 'hdr';
end

c = size(img, 2);

cubeSize = round(c / 3);

%CUBE_POS_Y
imgPosY = img(1:cubeSize,(cubeSize+1):(2*cubeSize),:);
imgPosY = imrotate(imgPosY, -270);
hdrimwrite(imgPosY,[cb_name,'_POS_Y.',cb_format]);

%CUBE_POS_X
hdrimwrite(img((cubeSize+1)    :(2*cubeSize),(cubeSize+1):(2*cubeSize),:),[cb_name,'_POS_X.',cb_format]);

%CUBE_NEG_Y
imgNegY = img((2*cubeSize+1):(3*cubeSize),(cubeSize+1):(2*cubeSize),:);
hdrimwrite(imrotate(imgNegY,-90),[cb_name,'_NEG_Y.',cb_format]);

%CUBE_NEG_X
imgNegX = img((3*cubeSize+1):(4*cubeSize),(cubeSize+1):(2*cubeSize),:);
for i=1:3
    imgNegX(:,:,i) = flipud(imgNegX(:,:,i));
    imgNegX(:,:,i) = fliplr(imgNegX(:,:,i));
end
hdrimwrite(imgNegX,[cb_name,'_NEG_X.',cb_format]);

%CUBE_POS_Z
hdrimwrite(img((cubeSize+1):(2*cubeSize),1:cubeSize,:),[cb_name,'_POS_Z.',cb_format]);

%CUBE_NEG_Z
hdrimwrite(img((cubeSize+1):(2*cubeSize),(2*cubeSize+1):(3*cubeSize),:),[cb_name,'_NEG_Z.', cb_format]);

end

author: banterle


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

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

% Last Modified by GUIDE v2.5 17-Dec-2018 16:42:03

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

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

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

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

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

    lvdData = varargin{2};
    setappdata(hObject, 'lvdData', lvdData);

	populateGUI(handles, action, lvdData);

    % Update handles structure
    guidata(hObject, handles);

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

function populateGUI(handles, action, lvdData)
    [tanksListStr, tanks] = lvdData.launchVehicle.getTanksListBoxStr();

    if(isempty(action.tank))
        indVal = 1;
    else
        indVal = find(tanks == action.tank, 1);
        if(isempty(indVal))
            indVal = 1;
        end
    end

    handles.tankCombo.String = tanksListStr;
    handles.tankCombo.Value = indVal;

    handles.massToAddText.String = fullAccNum2Str(action.massToAdd);

    optVar = action.optVar;
	if(isempty(optVar))
        optVar = AddMassToTankActionOptimVar(action);
	end

    handles.optCheckbox.Value = double(optVar.getUseTfForVariable());

    [lb, ub] = optVar.getAllBndsForVariable();
    handles.lbText.String = fullAccNum2Str(lb);
    handles.ubText.String = fullAccNum2Str(ub);

    optCheckbox_Callback(handles.optCheckbox, [], handles);

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

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

        [~, tanks] = lvdData.launchVehicle.getTanksListBoxStr();

        action.massToAdd = str2double(handles.massToAddText.String);
        action.tank = tanks(handles.tankCombo.Value);

        optVar = action.optVar;
        if(not(isempty(optVar))) %need to remove existing var if it exists
            lvdData.optimizer.vars.removeVariable(optVar);
        end

        %Set Opt T/F
        useTf = get(handles.optCheckbox,'Value');

        if(isempty(optVar))
            optVar = AddMassToTankActionOptimVar(action);
        end

        optVar.setUseTfForVariable(useTf);

        %Bnds
        lb = str2double(get(handles.lbText,'String'));
        ub = str2double(get(handles.ubText,'String'));

        optVar.setUseTfForVariable(true(size(lb))); %need this to get the full lb/set in there
        optVar.setBndsForVariable(lb, ub);
        optVar.setUseTfForVariable(useTf);

        lvdData.optimizer.vars.addVariable(optVar);

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

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

    if(isempty(errMsg))
        uiresume(handles.lvd_AddMassToTankActionGUI);
    else
        msgbox(errMsg,'Invalid Tank Inputs','error');
    end

function errMsg = validateInputs(handles)
    errMsg = {};

    val = str2double(get(handles.massToAddText,'String'));
    enteredStr = get(handles.massToAddText,'String');
    numberName = 'Mass to Add';
    lb = -Inf;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(val, numberName, lb, ub, isInt, errMsg, enteredStr);

    lwrBnd = str2double(get(handles.lbText,'String'));
    enteredStr = get(handles.lbText,'String');
    numberName = 'Mass to Add (Lower Bound)';
    lb = -Inf;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(lwrBnd, numberName, lb, ub, isInt, errMsg, enteredStr);

    uprBnd = str2double(get(handles.ubText,'String'));
    enteredStr = get(handles.ubText,'String');
    numberName = 'Mass to Add (Upper Bound)';
    lb = -Inf;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(uprBnd, numberName, lb, ub, isInt, errMsg, enteredStr);

    if(isempty(errMsg))
        uprBnd = str2double(get(handles.ubText,'String'));
        enteredStr = get(handles.ubText,'String');
        numberName = 'Mass to Add (Upper Bound)';
        lb = lwrBnd;
        ub = Inf;
        isInt = false;
        errMsg = validateNumber(uprBnd, numberName, lb, ub, isInt, errMsg, enteredStr);
    end

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

function massToAddText_Callback(hObject, eventdata, handles)
% hObject    handle to massToAddText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of massToAddText as text
%        str2double(get(hObject,'String')) returns contents of massToAddText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

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

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

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

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

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

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

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

% Hint: get(hObject,'Value') returns toggle state of optCheckbox
    if(get(hObject,'Value'))
        handles.lbText.Enable = 'on';
        handles.ubText.Enable = 'on';
    else
        handles.lbText.Enable = 'off';
        handles.ubText.Enable = 'off';
    end

function lbText_Callback(hObject, eventdata, handles)
% hObject    handle to lbText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of lbText as text
%        str2double(get(hObject,'String')) returns contents of lbText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

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

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

function ubText_Callback(hObject, eventdata, handles)
% hObject    handle to ubText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of ubText as text
%        str2double(get(hObject,'String')) returns contents of ubText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

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

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

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

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


clear;
%list of images
imglist=importdata('imglist.txt');
%minimum size of face
minsize=20;

%path of toolbox
caffe_path='E:/matlab';
pdollar_toolbox_path='E:/matlab/edges-master/toolbox-master'
caffe_model_path='./model'
addpath(genpath(caffe_path));
addpath(genpath(pdollar_toolbox_path));

%use cpu
%caffe.set_mode_cpu();
gpu_id=0;
caffe.set_mode_gpu();
caffe.set_device(gpu_id);

%three steps's threshold
threshold=[0.6 0.7 0.7]

%scale factor
factor=0.709;

%load caffe models
prototxt_dir =strcat(caffe_model_path,'/det1.prototxt');
model_dir = strcat(caffe_model_path,'/det1.caffemodel');
PNet=caffe.Net(prototxt_dir,model_dir,'test');
prototxt_dir = strcat(caffe_model_path,'/det2.prototxt');
model_dir = strcat(caffe_model_path,'/det2.caffemodel');
RNet=caffe.Net(prototxt_dir,model_dir,'test');
prototxt_dir = strcat(caffe_model_path,'/det3.prototxt');
model_dir = strcat(caffe_model_path,'/det3.caffemodel');
ONet=caffe.Net(prototxt_dir,model_dir,'test');
prototxt_dir =  strcat(caffe_model_path,'/det4.prototxt');
model_dir =  strcat(caffe_model_path,'/det4.caffemodel');
LNet=caffe.Net(prototxt_dir,model_dir,'test');
faces=cell(0);
for i=1:length(imglist)
    i
	img=imread(imglist{i});
	%we recommend you to set minsize as x * short side
	%minl=min([size(img,1) size(img,2)]);
	%minsize=fix(minl*0.1)
    tic
    [boudingboxes points]=detect_face(img,minsize,PNet,RNet,ONet,LNet,threshold,false,factor);
	toc
    faces{i,1}={boudingboxes};
	faces{i,2}={points'};
	%show detection result
	numbox=size(boudingboxes,1);
	imshow(img)
	hold on;
	for j=1:numbox
		plot(points(1:5,j),points(6:10,j),'g.','MarkerSize',10);
		r=rectangle('Position',[boudingboxes(j,1:2) boudingboxes(j,3:4)-boudingboxes(j,1:2)],'Edgecolor','g','LineWidth',3);
    end
    hold off;
	pause
end
save result box landmark


text:
function [Q] = VBA_spm_Q(a,n,q)
% returns an (n x n) (inverse) autocorrelation matrix for an AR(p) process
% FORMAT [Q] = spm_Q(a,n,q)
%
% a  - vector of (p) AR coefficients
% n  - size of Q
% q  - switch to return inverse autocorrelation or precision [default q = 0]
%__________________________________________________________________________
% spm_Q uses a Yule-Walker device to compute K where:
%
% y = K*z
%
% such that y is an AR(p) process generated from an i.i.d innovation
% z.  This means
%
% cov(y) = <K*z*z'*K> = K*K'
%
% If called with q ~= 0, a first order process is assumed when evaluating
% the precision (inverse covariance) matrix; i.e., a = a(1)
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_Q.m 5838 2014-01-18 18:40:37Z karl $

% default
%--------------------------------------------------------------------------
try, q; catch, q = 0; end

if q

    % compute P (precision)
    %----------------------------------------------------------------------
    A    = [-a(1) (1 + a(1)^2) -a(1)];
    Q    = spdiags(ones(n,1)*A,(-1:1),n,n);

else

    % compute Q (covariance)
    %----------------------------------------------------------------------
    p    = length(a);
    A    = [1 -a(:)'];
    P    = spdiags(ones(n,1)*A,-(0:p),n,n);
    K    = inv(P);
    K    = K.*(abs(K) > 1e-4);
    Q    = K*K';
    Q    = toeplitz(Q(:,1));

end

author: MBB-team

repo: VBA-toolbox


text:
function [tag, val, typ, dep, chk, cj] = harvest(item, cj, dflag, rflag)

% function [tag, val, typ, dep, chk, cj] = harvest(item, cj, dflag, rflag)
% Harvest a cfg_repeat object.
% Input arguments:
% item  - item to be harvested
% cj    - configuration tree (passed unmodified)
% dflag - if true, harvest defaults tree, otherwise filled tree
% rflag - if true, resolve dependencies in leaf nodes
% Output arguments:
% tag - tag of harvested item
% val - harvested value
% typ - class of harvested item (currently unused)
% dep - list of unresolved dependencies
% chk - meaningful if ~dflag and all dependencies are resolved. Then it
%       returns success status of this items .check function and its
%       childrens check functions. A job is ready to run if all
%       dependencies are resolved and chk status is true.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: harvest.m 5678 2013-10-11 14:58:04Z volkmar $

rev = '$Rev: 5678 $'; %#ok

typ = class(item);
tag = gettag(item);
dep = []; % Placeholder for dependencies. Will be classified during
          % first call to dep_add.
chk = ~dflag && rflag;

tname = treepart(item, dflag);
ntgt_input = substruct('.', tname, '{}', {});
citems = subsref(item, ntgt_input(1));
if numel(item.values)==1 && isa(item.values{1},'cfg_branch') && ~item.forcestruct,
    if numel(citems) == 0
        % return empty struct
        cargs = cell(2*numel(item.values{1}.val),1);
        for i=1:numel(item.values{1}.val),
            cargs{2*i-1} = gettag(item.values{1}.val{i});
            cargs{2*i}   = {};
        end;
        val = struct(cargs{:});
    end;
    if ~dflag
        njtsubs(1).type = '()';
    end;
else
    if ~dflag
        val = cell(size(citems));
        njtsubs(1).type = '{}';
    else
        val = struct([]);
    end;
end;
for i=1:numel(citems),
    [ctag, cval, unused, cdep, cchk, cj] = harvest(citems{i}, cj, dflag, rflag);
    if numel(item.values)==1 && isa(item.values{1},'cfg_branch') && ~item.forcestruct,
        % FIXME: don't know how to best preinit this without raising
        % warnings about incompatible assignments ...
        val(i) = cval;
    else
        if numel(item.values)>1 || item.forcestruct,
            if iscell(cval)
                cval = struct(ctag,{cval});
            else
                cval = struct(ctag,cval);
            end;
            if ~dflag
                njtsubs(2).type = '.';
                njtsubs(2).subs  = ctag;
            end;
        end;
        if dflag
            % return a struct containing defaults of all child nodes
            % instead of a cell array. This makes defaults easier to
            % read.
            if numel(item.values)>1 || item.forcestruct,
                val(1).(ctag) = cval.(ctag);
            else
                val = {cval};
            end;
        else
            val{i} = cval;
        end;
    end;
    if ~dflag && ~isempty(cdep)
        % augment cdep tsubs references
        ntgt_input(2).subs = {i};
        njtsubs(1).subs = {i};
        dep = dep_add(cdep, dep, ntgt_input, njtsubs);
    end;
    chk = chk && cchk;
end;
if chk
    chk = docheck(item, val);
end;


text:
function P = sammon(D, P, varargin)

%SAMMON Computes Sammon's mapping of a data set.
%
% P = sammon(D, P, [value], [mode], [alpha], [Mdist])
%
%  P = sammon(D,2);            % projection to 2-dim space
%  P = sammon(sMap,3);         % projects the codebook vectors
%  P = sammon(sMap,3,[],[],[],Md) % uses distance matrix Md
%  som_grid(sMap,'Coord',P)    % visualization of map projection
%
%  Input and output arguments ([]'s are optional):
%   D        (matrix) size dlen x dim, data to be projected
%            (struct) data or map struct
%   P        (scalar) output dimension
%            (matrix) size dlen x odim, initial projection matrix
%   [value]  (scalar) all different modes (the next argument) require
%                     a value, default = 100
%   [mode]   (string) 'steps' or 'errlimit' or 'errchange' or 'seconds',
%                     see below, default is 'steps'
%   [alpha]  (scalar) iteration step size, default = 0.2
%   [Dist]   (matrix) pairwise distance matrix, size dlen x dlen.
%                     If the distances in the input space should
%                     be calculated otherwise than as euclidian
%                     distances, the distance from each vector
%                     to each other vector can be given here,
%                     size dlen x dlen. For example PDIST
%                     function can be used to calculate the
%                     distances: Dist = squareform(pdist(D,'mahal'));
%
%   P        (matrix) size dlen x odim, the projections
%
% The output dimension must be 2 or higher but (naturally) lower
% than data set dimension.
%
% The mode argument determines the end condition for iteration. If
% the mode argument is used, also the value argument has to be
% specified. Different mode possibilities are:
% 'steps'      the iteration is terminated when it is run <value>
% 'errlimit'   steps, the iteration is terminated when projection error
%              is lower than <value>,
% 'errchange'  the iteration is terminated when change between
%              projection error on two successive iteration rounds
%	       is less than <value> percent of total error, and
% 'seconds'    the iteration is terminated after <value> seconds
%              of iteration.
%
% See also CCA, PCAPROJ, SOM_GRID.

% Reference: Sammon, J.W. Jr., "A nonlinear mapping for data
%   structure analysis", IEEE Transactions on Computers, vol. C-18,
%   no. 5, 1969, pp. 401-409.

% Contributed to SOM Toolbox vs2, February 2nd, 2000 by Juha Vesanto
% Copyright (c) by Juha Vesanto
% http://www.cis.hut.fi/projects/somtoolbox/

% juuso 040100

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% check arguments

error(nargchk(2, 6, nargin));  % check no. of input arguments is correct

% input data
if isstruct(D),
  if isfield(D, 'data'),         D = D.data; % data struct
  elseif isfield(D, 'codebook'), D = D.codebook; % map struct
  else error('Invalid structure');
  end
end
if any(isnan(D(:))),
  error('Cannot make Sammon''s projection for data with unknown components')
end

% compute data dimensions
orig_si = size(D);
dim = orig_si(end);
noc = prod(orig_si)/dim;
if length(orig_si)>2, D = reshape(D,[noc dim]); end

% output dimension / initial projection matrix
if prod(size(P))==1,
  odim = P;
  P = rand(noc,odim)-0.5;
else
  si = size(P);
  odim = si(end);
  if prod(si) ~= noc*odim,
    error('Initial projection matrix size does not match data size');
  end
  if length(si)>2, P = reshape(P,[noc odim]); end
  inds = find(isnan(P));
  if length(inds), P(inds) = rand(size(inds)); end
end
if odim > dim | odim < 2,
  error('Output dimension must be within [2, dimension of data]');
end

% determine operating mode
if nargin < 3 | isempty(varargin{1}) | isnan(varargin{1}), value=100;
else value = varargin{1};
end

if nargin < 4 | isempty(varargin{2}) | isnan(varargin{2}), mode='steps';
else mode = varargin{2};
end
switch mode,
case 'steps',     runlen = value;
case 'errlimit',  errlimit = value;
case 'errchange', errchange = value; e_prev = 0;
case 'seconds',   endtime = value;
otherwise, error(['Illegal mode: ' mode]);
end

% iteration step size
if nargin > 4, alpha = varargin{3}; else alpha = NaN; end
if isempty(alpha) | isnan(alpha), alpha = 0.2; end

% mutual distances
if nargin > 5, Mdist = varargin{4}; else Mdist = []; end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% initialization

% these are used quite frequently
noc_x_1  = ones(noc, 1);
odim_x_1 = ones(odim,1);

% compute mutual distances between vectors
if isempty(Mdist) | all(isnan(Mdist(:))),
  fprintf(2, 'computing mutual distances\r');
  dim_x_1 = ones(dim,1);
  for i = 1:noc,
    x = D(i,:);
    Diff = D - x(noc_x_1,:);
    N = isnan(Diff);
    Diff(find(N)) = 0;
    Mdist(:,i) = sqrt((Diff.^2)*dim_x_1);
    N = find(sum(N')==dim); %mutual distance unknown
    if ~isempty(N), Mdist(N,i) = NaN; end
  end
else
  % if the distance matrix is output from PDIST function
  if size(Mdist,1)==1, Mdist = squareform(Mdist); end
  if size(Mdist,1)~=noc,
    error('Mutual distance matrix size and data set size do not match');
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% action

if strcmp(mode, 'seconds'), tic; end;
fprintf(2, 'iterating                    \r');

% sammon iteration

x  = P ;
xu = zeros(noc, odim);
xd = zeros(noc, odim);
dq = zeros(noc, 1);
dr = zeros(noc, 1);

i = 0;
ready = 0;
while ~ready
  for j = 1:noc,
    xd      = -x + x(j*noc_x_1,:);
    xd2     = xd.^2;
    dpj     = sqrt(sum(xd2'))';
    dq      = Mdist(:,j) - dpj;
    dr      = Mdist(:,j) .* dpj;
    ind     = find(dr ~= 0);
    term    = dq(ind) ./ dr(ind);
    e1      = sum(xd(ind,:) .* term(:,odim_x_1));
    term2   = ((1.0 + dq(ind) ./ dpj(ind)) ./ dpj(ind)) ./ dr(ind);
    e2      = sum(term) - sum(xd2(ind,:) .* term2(:,odim_x_1));
    xu(j,:) = x(j,:) + alpha * e1 ./ abs(e2);
  end

  % move the center of mass to the center

  c = sum(xu) / noc;
  x = xu - c(noc_x_1, :);

  i = i + 1;

  % compute mapping error
  % doing this adds about 25% to computing time
  if 0,
    e = 0; tot = 0;
    for j = 2:noc,
      d   = Mdist(1:(j - 1), j);
      tot = tot + sum(d);
      ind = find(d ~= 0);
      xd  = -x(1:(j - 1), :) + x(j * ones(j - 1, 1), :);
      ee  = d - sqrt(sum(xd'.^2))';
      e   = e + sum(ee(ind).^2 ./ d(ind));
    end
    e = e/tot;
    fprintf(2, '\r%d iterations, error %f', i, e);
  else
    fprintf(2, '\r%d iterations', i);
  end

  % determine is the iteration ready

  switch mode
    case 'steps',
      if i == runlen, ready = 1; end;
    case 'errlimit',
      if e < errlimit, ready = 1; end;
    case 'errchange',
      if i > 1
	change = 100 * abs(e - e_prev) / e_prev;
	if change < errchange, ready = 1; end;
	fprintf(2, ', change of error %f %%    ', change);
      end
      e_prev = e;
    case 'seconds'
      if toc > endtime, ready = 1; end;
      fprintf(2, ', elapsed time %f seconds  ', toc);
  end
  fprintf(2, '        ');

  % If you want to see the Sammon's projection plotted (in 2-D and 3-D case),
  % execute the code below; it is not in use by default to speed up
  % computation.
  if 0,
    clf
    if odim == 1,     plot(x(:,1), noc_x_1, 'o');
    elseif odim == 2, plot(x(:,1), x(:,2), 'o');
    else              plot3(x(:,1), x(:,2), x(:,3), 'o')
    end
    drawnow
  end
end

fprintf(2, '\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% clean up

% reshape
orig_si(end) = odim;
P = reshape(x, orig_si);

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

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% The COBRAToolbox: testChangeIBMCplexParams.m
%
% Purpose:
%    This script tests whether input parameters for the IBM ILOG Cplex solver
%    are in effect when calling solveCobraLP.m, and also tests directly the function
%    setCplexParam used to change the IBM ILOG Cplex parameters in solveCobraLP.m

% save the current path
currentDir = pwd;

%get the warning settings and turn on all warnings
cwarn = warning;
warning('on');

% initialize the test
fileDir = fileparts(which('testChangeIBMCplexParams'));
cd(fileDir);

try
    ibm_cplex = changeCobraSolver('ibm_cplex', 'LP');
catch
    ibm_cplex = false;
end
if ibm_cplex
    LP = struct();
    LP.A = 1;
    LP.b = 1;
    LP.lb = -1000;
    LP.ub = 1000;
    LP.c = 1;
    LP.osense = -1;
    LP.csense = 'L';
    sol = solveCobraLP(LP);
    assert(sol.full == 1)

    if isunix
        % change the time limit for the solver
        % Note: On windows, the timelimit parameter has no effect
        CplexParam = struct('timelimit', 0);
        sol = solveCobraLP(LP, CplexParam);
        % no solution because of time limit
        assert(isempty(sol.full) & sol.origStat == 11)
    end

    % print parameters not recognized by Cplex
    CplexParam.notexist = 1;
    diary('testChangeIBMCplexParams.txt');
    sol = solveCobraLP(LP, CplexParam, 'printLevel', 1);
    diary off
    fid = fopen('testChangeIBMCplexParams.txt', 'r');
    text1 = '';
    line = fgetl(fid);
    while ~isequal(line, -1)
        text1 = [strtrim(text1), ' ', line];
        line = fgetl(fid);
    end
    fclose(fid);
    % check if the warning is printed
    assert(~isempty(strfind(text1, 'Warning: *.notexist cannot be identified as a valid cplex parameter. Ignore.')))
    delete('testChangeIBMCplexParams.txt');

    % print parameters that cannot be uniquely identified
    CplexParam = struct('timelimit', 0, 'feasibility', 1e-8);
    diary('testChangeIBMCplexParams.txt');
    sol = solveCobraLP(LP, CplexParam, 'printLevel', 1);
    diary off
    fid = fopen('testChangeIBMCplexParams.txt', 'r');
    text1 = '';
    line = fgetl(fid);
    while ~isequal(line, -1)
        text1 = [strtrim(text1), ' ', line];
        line = fgetl(fid);
    end
    fclose(fid);
    % check if the warning is printed
    assert(~isempty(strfind(text1, 'Warning: *.feasibility cannot be uniquely identified as a valid cplex parameter. Ignore.')))
    delete('testChangeIBMCplexParams.txt');

    % test changing other Cplex parameters
    LP = Cplex('test setCplexParam');
    CplexParam = struct();
    [CplexParam.simplex.display, CplexParam.tune.display] = deal(0);
    [CplexParam.simplex.tolerances.optimality, CplexParam.simplex.tolerances.feasibility] = deal(1e-9, 1e-8);
    % change parameters using setCplexParam
    LP = setCplexParam(LP, CplexParam);
    % Cplex internal method for getting all changed parameters
    ChangedParams = LP.getChgParam;

    % test also with .Cur at the end of the input struct
    LP = Cplex('test setCplexParam');
    CplexParam = struct();
    [CplexParam.simplex.display.Cur, CplexParam.tune.display.Cur] = deal(0);
    [CplexParam.simplex.tolerances.optimality.Cur, CplexParam.simplex.tolerances.feasibility.Cur] = deal(1e-9, 1e-8);
    % change parameters using setCplexParam
    LP = setCplexParam(LP, CplexParam);
    % Cplex internal method for getting all changed parameters
    ChangedParams2 = LP.getChgParam;

    % the changed parameters should be the same using either input structures
    assert(isequal(ChangedParams, ChangedParams2))
    % check if the parameters are really changed
    assert(isfield(ChangedParams, 'simplex') && isfield(ChangedParams.simplex, 'display') ...
        && ChangedParams.simplex.display.Cur == 0 && isfield(ChangedParams.simplex, 'tolerances') ...
        && isfield(ChangedParams.simplex.tolerances, 'optimality') && ChangedParams.simplex.tolerances.optimality.Cur == 1e-9 ...
        && isfield(ChangedParams.simplex.tolerances, 'feasibility') && ChangedParams.simplex.tolerances.feasibility.Cur == 1e-8 ...
        && isfield(ChangedParams, 'tune') && isfield(ChangedParams.tune, 'display') && ChangedParams.tune.display.Cur == 0)
end

%Reset the warning settings
warning(cwarn);

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function[res] = prependNotEmpty(str, prefix)
% [res] = prependNotEmpty(str, prefix)
%
% Prepend a prefix to a string, if it is not empty.
%
% Copyrights by Holger Caesar, 2015

if ~isempty(str),
    res = [prefix, str];
else
    res = '';
end;

author: nightrome

repo: matconvnet-calvin


text:
function gui_LoadResult
% function pf_load
%   Opens a open dialog to load result data
%   and invokes the display dialog
%
% Danijel Schorlemmer
% August 22, 2002

global bDebug;
if bDebug
  report_this_filefun(mfilename('fullpath'));
end

% Invoke the open dialog
[newfile, newpath] = uigetfile('*.mat', 'Load calculated data');
% Cancel pressed?
if newfile == 0
  return;
end
% Everything ok?
newfile = [newpath newfile];
if length(newfile) > 1
  % Load the data
  load(newfile, 'vResults');
  % Open dialog
  hResultFig = gui_result(vResults);
end

author: CelsoReyes


text:
function segment_length = p14_boundary_segment_length ( segment_index, h )

%*****************************************************************************80
%
%% P14_BOUNDARY_SEGMENT_LENGTH returns boundary segment lengths in problem 14.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer SEGMENT_INDEX, the index of one of the boundary segments.
%
%    Input, real H, the suggested spacing between points.
%
%    Output, integer SEGMENT_LENGTH, the number of points in the segment.
%
  n1 = 90;
  n2 = 11;

  v1 = [ ...
    316.43027, 404.47559; ...
    291.04946, 400.70917; ...
    265.16504, 409.77890; ...
    241.46794, 402.40310; ...
    216.55145, 396.52064; ...
    163.28492, 411.37102; ...
    142.81752, 391.16355; ...
    111.95404, 346.70264; ...
    100.03538, 325.72710; ...
    103.98723, 302.51587; ...
    128.72978, 285.72802; ...
    147.49111, 266.23345; ...
    196.65261, 242.24055; ...
    213.56835, 221.67192; ...
    226.49969, 198.09326; ...
    248.37126, 183.50473; ...
    262.21952, 165.39102; ...
    278.42330, 149.91715; ...
    300.71846, 145.82601; ...
    311.12698, 166.71094; ...
    326.66315, 184.58335; ...
    359.78574, 225.48049; ...
    357.08892, 252.88958; ...
    358.76685, 285.34403; ...
    361.50834, 303.71287; ...
    371.68926, 314.92452; ...
    380.49890, 324.58632; ...
    396.37634, 328.88990; ...
    412.59116, 327.25238; ...
    425.48394, 315.28623; ...
    435.84305, 302.44664; ...
    458.34025, 297.55121; ...
    479.66439, 288.99238; ...
    493.09812, 270.20636; ...
    518.87309, 264.56427; ...
    547.18014, 268.18846; ...
    600.49708, 240.62570; ...
    625.96183, 238.40347; ...
    633.90530, 260.70629; ...
    621.50451, 285.88914; ...
    576.87224, 322.14121; ...
    570.51915, 348.85423; ...
    567.16400, 378.24075; ...
    558.00668, 406.86552; ...
    565.19008, 435.75599; ...
    567.56437, 465.33407; ...
    550.87626, 490.96358; ...
    532.98174, 515.84491; ...
    500.66817, 551.89078; ...
    478.75120, 562.17222; ...
    430.03371, 583.94286; ...
    401.20454, 587.69910; ...
    368.32214, 581.10110; ...
    354.26303, 585.86085; ...
    346.75200, 601.10367; ...
    332.85137, 628.74602; ...
    308.02188, 645.84180; ...
    295.52344, 647.18525; ...
    286.51519, 651.60328; ...
    285.98846, 662.07339; ...
    298.93455, 665.66316; ...
    301.70226, 682.79570; ...
    278.65857, 689.63850; ...
    266.25737, 712.11005; ...
    287.28701, 732.77147; ...
    318.19548, 736.85151; ...
    343.83067, 753.60957; ...
    375.53164, 758.35231; ...
    405.73444, 768.98687; ...
    406.33873, 785.59001; ...
    378.35436, 789.44240; ...
    350.02151, 795.02238; ...
    338.68030, 788.87325; ...
    325.67930, 786.10177; ...
    319.05995, 798.04657; ...
    301.78158, 795.34254; ...
    280.69272, 773.86634; ...
    254.55844, 758.02898; ...
    234.07759, 737.42090; ...
    218.38337, 711.41500; ...
    220.99086, 682.17833; ...
    224.50640, 651.96297; ...
    240.25971, 631.36117; ...
    259.86174, 612.60253; ...
    291.85381, 556.70385; ...
    315.52139, 537.56387; ...
    341.63663, 520.12519; ...
    351.37130, 458.75372; ...
    349.33183, 431.31454; ...
    328.80465, 412.43055 ]';
  v2 = [ ...
    238.64853, 266.58978; ...
    235.14026, 287.95183; ...
    238.20736, 303.46785; ...
    250.13902, 303.71290; ...
    258.51675, 297.46973; ...
    274.55300, 291.27357; ...
    284.66230, 280.72063; ...
    279.73288, 267.83455; ...
    270.68478, 255.55440; ...
    255.73801, 249.16872; ...
    241.72690, 256.73448 ]';

  if ( h <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P14_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );
    fprintf ( 1, '  Nonpositive H = %f\n', h );
    error ( 'P14_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );
  end

  if ( segment_index == 1 )

    length = polyloop_length_nd ( 2, n1, v1 ) ;

  elseif ( segment_index == 2 )

    length = polyloop_length_nd ( 2, n2, v2 );

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P14_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );
    fprintf ( 1, '  Illegal SEGMENT_INDEX = %d\n', segment_index );
    error ( 'P14_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );

  end

  n = round ( length / h );

  segment_length = n;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [X,P]=subPx(g,x,gMax,N)
% clear, load subPxData.mat, N=100;
xi=linspace(x(1),x(end),N*length(x));
gi = interp1(x,g,xi,'spline');
[pk,loc]=peakDetect(gi,0,0);
% in case multiple peaks are detected, the one closest to gMax is chosen
if length(loc)>1
    disp('WARNING: multiple peaks detected')
    [val,id]=min(abs(xi(loc)-x(find(g==gMax))));
else
    id=1;
end
X=xi(loc(id)); P=pk(id);
% figure(1), plot(x,g,'o',xi,gi,'.',xi(L),gi(L),'o')