text:
%clear all
global planC

directoryname = uigetdir;
dirdata = dir(directoryname);
files={dirdata.name};
ln=length(files);
files=files(3:ln); ln=ln-2;

structNum_roi=1;
structNum_40p=2;
M=[];
Thresholds=[40:10:80];
for i=20:20
    name=files{i}
    load([directoryname,'\',name]);
    %     roiVol = getStructureVol(structNum_40p);
    %    [Dx, Vx(i,:), mean_suv_40p, max_suv_40p, min_suv_40p, slope_ivh_40p(i)]=analyze_ivh1(structNum_40p, 1, 0);
    %slope(i) = calc_slope_grigsby(structNum_roi,Thresholds);
    [energy_40p,contrast_40p,Entropy_40p,Homogeneity_40p,standard_dev_40p,Ph_40p] = getHaralicParams(structNum_40p);
    %     [Eccentricity_40p,EulerNumber_40p,Solidity_40p,Extent_40p]=getShapeParams(structNum_40p);
    %     row=[roiVol, mean_suv_40p, standard_dev_40p, max_suv_40p, min_suv_40p, slope_ivh_40p, Dx, Vx,energy_40p,contrast_40p,Entropy_40p,Homogeneity_40p,...
    %         Eccentricity_40p,EulerNumber_40p,Solidity_40p,Extent_40p];
    %     M=[M;row];
end
%xlswrite('pet_hetro_cervix_grigsby.xls',M)

% clear PlanC
% save pet_hetro_cervix

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function f = laplacian(f)
%LAPLACIAN   Scalar laplacian of a SPHEREFUN.
%   L = LAPLACIAN(F) computes the spherical Laplacian of the SPHEREFUN F.
%
% This is gvien, in intrinsic spherical coordinates as
%
% laplacian(F) = d_{theta} d_{theta} F + tan(theta) d_{theta} F +
%                1/sin(theta) + d_{lambda} d_{lambda} F
%
% where 0 <= theta <= pi is the co-latitude variable and
% -pi <= lambda <= pi  is the longitude variable.
%
% See also GRADIENT, DIVERGENCE.

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

% We implement the computation of the Laplacian using the spherefun/diff
% command.  This function computes tangential derivatives in the x, y, and
% z directions, so that the Laplacian can be implemented as we usually
% think of it in 3D, i.e.
% laplacian(f) = f_xx + f_yy + f_zz.

% Compute the derivatives
fxx = diff(f, 1, 2);
fyy = diff(f, 2, 2);
fzz = diff(f, 3, 2);

% Determine which derivative has the most coefficients in its
% representation for the columns and rows.
[mxx,nxx] = length(fxx);
[myy,nyy] = length(fyy);
[mzz,nzz] = length(fzz);

m = max([mxx myy mzz]);
n = max([nxx nyy nzz]);

% Ensure m and n are even:
m = m + mod(m, 2);
n = n + mod(n, 2);

% Construct a spherefun by sampling the derivatives on the same size grid.
f = spherefun(sample(fxx, m, n/2) + sample(fyy, m, n/2) + ...
               sample(fzz, m, n/2));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% sieveEr: returns a vector with prime numbers from 2 up to N
% assumes: N >= 2
function  y = sieveER(N)
  % precondition
  assert(N >= 2,"N must be >= 2")
  tmp = [false,true(1,N-1)]; % indexed by all numbers from 2 up to N

  % labels all composite number with false
  for i = 2 : 1 : sqrt(N)
    if (tmp(i))
      for j = i^2 : i : N
        tmp(j) = false;
      endfor
    endif
  endfor

  % fills up all prime numbers in vector y
  y = find(tmp);

endfunction

author: TheAlgorithms

repo: MATLAB-Octave


text:
function blas1_c_test05 ( )

%*****************************************************************************80
%
%% TEST05 tests CDOTC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2006
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  x = [ ...
     2.0 - 1.0 * i, ...
    -4.0 - 2.0 * i, ...
     3.0 + 1.0 * i, ...
     2.0 + 2.0 * i, ...
    -1.0 - 1.0 * i ];
  y = [ ...
    -1.0 + 0.0 * i, ...
     0.0 - 3.0 * i, ...
     4.0 + 0.0 * i, ...
    -3.0 + 4.0 * i, ...
    -2.0 + 0.0 * i ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST05\n' );
  fprintf ( 1, '  CDOTC computes the conjugated dot product of\n' );
  fprintf ( 1, '  two complex vectors.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  X =\n' );
  fprintf ( 1, '\n' );
  for j = 1 : n
    fprintf ( 1, '  %6d  %10f  %10f\n', j, real ( x(j) ), imag ( x(j) ) );
  end

  x_norm = cdotc ( n, x, 1, x, 1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The square of the norm of X, computed as\n' );
  fprintf ( 1, '  CDOTC(X,X) = %f  %f\n', real ( x_norm ), imag ( x_norm ) );

  xy_dot = cdotc ( n, x, 1, y, 1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Y = \n' );
  fprintf ( 1, '\n' );
  for j = 1 : n
    fprintf ( 1, '  %6d  %10f  %10f\n', j, real ( y(j) ), imag ( y(j) ) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The dot product X.Y* is %f  %f\n', real ( xy_dot ), imag ( xy_dot ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function output = callsdpnal(interfacedata)

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
x0      = interfacedata.x0;
ub      = interfacedata.ub;
lb      = interfacedata.lb;

% Bounded variables converted to constraints
if ~isempty(ub)
    [F_struc,K] = addStructureBounds(F_struc,K,ub,lb);
end

[blk,A,C,b,oldKs]=sedumi2sdpt3(F_struc(:,1),F_struc(:,2:end),c,K,1);

if options.savedebug
    ops = options.sdpnal;
    save sdpnaldebug blk A C b ops -v6
end

if options.showprogress;showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);end
solvertime = tic;
if options.verbose==0
    evalc('[obj,X,s,y,Z,Z2,y2,v,info,runhist] = sdpnalplus(blk,A,C,b,[],[],[],[],[],options.sdpnal);');
else
    [obj,X,s,y,Z,Z2,y2,v,info,runhist] = sdpnalplus(blk,A,C,b,[],[],[],[],[],options.sdpnal);
end
solvertime = toc(solvertime);

% Create YALMIP dual variable and slack
Dual = [];
Slack = [];
top = 1;
if K.f>0
    Dual = [Dual;X{top}(:)];
    Slack = [Slack;Z{top}(:)];
    top = top+1;
end
if K.l>0
    Dual = [Dual;X{top}(:)];
    Slack = [Slack;Z{top}(:)];
    top = top + 1;
end
if K.q(1)>0
    Dual = [Dual;X{top}(:)];
    Slack = [Slack;Z{top}(:)];
    top = top + 1;
end
if K.s(1)>0
    % Messy format in SDPT3 to block and sort small SDPs
    u = blk(:,1);
    u = find([u{:}]=='s');
    s = 1;
    for top = u
        ns = blk(top,2);ns = ns{1};
        k = 1;
        for i = 1:length(ns)
            Xi{oldKs(s)} = X{top}(k:k+ns(i)-1,k:k+ns(i)-1);
            Zi{oldKs(s)} = Z{top}(k:k+ns(i)-1,k:k+ns(i)-1);
            s = s + 1;
            k = k+ns(i);
        end
    end
    for i = 1:length(Xi)
        Dual = [Dual;Xi{i}(:)];
        Slack = [Slack;Zi{i}(:)];
    end
end
Primal = -y;  % Primal variable in YALMIP

% No error code available
if isfield(info,'termcode')
    switch info.termcode
        case -1
            problem = 4;
        case -2
            problem = 3;
        case 0
            problem = 0;
        case 1
            problem = 5;
        case 2
            problem = 3;
        otherwise
            problem = 11;
    end
else
    if isfield(info,'msg')
        if isequal(info.msg,'maximum iteration reached')
            problem = 3;
        else
            problem = 9;
        end
    end
end

infostr = yalmiperror(problem,interfacedata.solver.tag);

if options.savesolveroutput
    solveroutput.obj = obj;
    solveroutput.X = X;
    solveroutput.y = y;
    solveroutput.Z = Z;
    solveroutput.info = info;
    solveroutput.runhist = runhist;
 else
    solveroutput = [];
end

if options.savesolverinput
    solverinput.blk = blk;
    solverinput.A   = A;
    solverinput.C   = C;
    solverinput.b   = b;
    solverinput.options   = options.sdpnal;
else
    solverinput = [];
end

% Standard interface
output = createOutputStructure(Primal,Dual,[],problem,infostr,solverinput,solveroutput,solvertime);

function [F_struc,K] = deblock(F_struc,K);
X = any(F_struc(end-K.s(end)^2+1:end,:),2);
X = reshape(X,K.s(end),K.s(end));
[v,dummy,r,dummy2]=dmperm(X);
blks = diff(r);

lint = F_struc(1:end-K.s(end)^2,:);
logt = F_struc(end-K.s(end)^2+1:end,:);

newlogt = [];
for i = 1:size(logt,2)
    temp = reshape(logt(:,i),K.s(end),K.s(end));
    temp = temp(v,v);
    newlogt = [newlogt temp(:)];
end
logt = newlogt;

pattern = [];
for i = 1:length(blks)
    pattern = blkdiag(pattern,ones(blks(i)));
end

F_struc = [lint;logt(find(pattern),:)];
K.s(end) = [];
K.s = [K.s blks];
K.m = blks;

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function a = fibonacci2 ( n )

%*****************************************************************************80
%
%% FIBONACCI2 returns the FIBONACCI2 matrix.
%
%  Example:
%
%    N = 5
%
%    0 1 0 0 0
%    1 1 0 0 0
%    0 1 1 0 0
%    0 0 1 1 0
%    0 0 0 1 1
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    A is tridiagonal.
%
%    Because A is tridiagonal, it has property A (bipartite).
%
%    A is banded, with bandwidth 3.
%
%    A is integral: int ( A ) = A.
%
%    A is a zero/one matrix.
%
%    If N = 1 then
%      det ( A ) = 0
%    else
%      det ( A ) = (-1)^(N-1)
%
%    If 1 < N, then A is unimodular.
%
%    For 2 <= N, A has the eigenvalues:
%
%      PHI   (once),
%      1     (N-2) times,
%      1-PHI (once).
%
%    When applied to a Fibonacci1 matrix B, the Fibonacci2 matrix
%    A produces the "next" Fibonacci1 matrix C = A*B.
%
%    Let PHI be the golden ratio (1+sqrt(5))/2.
%
%    For 2 <= N, the eigenvalues and eigenvectors are:
%
%    LAMBDA(1)     = PHI,     vector = (1,PHI,PHI^2,...PHI^(N-1));
%    LAMBDA(2:N-1) = 1        vector = (0,0,0,...,0,1);
%    LAMBDA(N)     = 1 - PHI. vector = ((-PHI)^(N-1),(-PHI)^(N-2),...,1)
%
%    Note that there is only one eigenvector corresponding to 1.
%    Hence, for 3 < N, the matrix is defective.  This fact means,
%    for instance, that the convergence of the eigenvector in the power
%    method will be very slow.
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

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

      if ( i == 1 )

        if ( j == 2 )
          a(i,j) = 1.0;
        else
          a(i,j) = 0.0;
        end

      else

        if ( j == i-1 | j == i )
          a(i,j) = 1.0;
        else
          a(i,j) = 0.0;
        end

      end

    end
  end

  return
end


text:
function[y] = expand(x)

[N,M] = size(x);
N = N*2;
M = M*2;

y = zeros(N,M);
y(1:2:N,1:2:M) = x;
y(2:2:N,2:2:M) = x;
y(1:2:N,2:2:M) = x;
y(2:2:N,1:2:M) = x;

author: thomaskuestner


text:
function o = mean_CS(q,CS)
% fast mean of
%
% Syntax
%   m = mean(o)
%
% Input
%  o        - list of @orientation
%
% Output
%  m      - mean @orientation
%
% See also
% orientation/mean

if length(q) > 1
  q = project2FundamentalRegion(q,CS,q.subSet(1));
  q = mean(q);
end

o = orientation(q,CS);

author: mtex-toolbox

repo: mtex

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


text:
function N = NormDirection(vertex, tri)

% norm of each triangles
pt1 = vertex(:, tri(1, :));
pt2 = vertex(:, tri(2, :));
pt3 = vertex(:, tri(3, :));
n_tri = cross(pt1 - pt2, pt1 - pt3);

% norm of each vertex
N = zeros(3, size(vertex, 2));
% for i = 1 : size(tri, 2)
%     N(:, tri(1, i)) = N(:, tri(1, i)) + n_tri(:, i);
%     N(:, tri(2, i)) = N(:, tri(2, i)) + n_tri(:, i);
%     N(:, tri(3, i)) = N(:, tri(3, i)) + n_tri(:, i);
% end

N = Tnorm_VnormC(double(n_tri), double(tri), double(size(tri,2)), double(size(vertex,2)));

% normalize to unit length
mag = sum(N .* N);
% deal with zero vector
co = find(mag == 0);
mag(co) = 1;
N(1, co) = ones(length(co),1);
N = N ./ sqrt(repmat(mag, 3, 1));
N = -N;

author: XgTu


text:
function [dat] = read_nexstim_nxe(filename, begsample, endsample, chanindx)

% READ_NEXSTIM_NXE reads specified samples from a NXE continuous datafile
%
% Use as
%   [hdr] = read_nexstim_nxe(filename)
% where
%    filename        name of the datafile, including the .bdf extension
% This returns a header structure with the following elements
%   hdr.Fs           sampling frequency
%   hdr.nChans       number of channels
%   hdr.nSamples     number of samples per trial
%   hdr.nSamplesPre  number of pre-trigger samples in each trial
%   hdr.nTrials      number of trials
%   hdr.label        cell-array with labels of each channel
%
% Or use as
%   [dat] = read_nexstim_nxe(filename, begsample, endsample, chanindx)
% where
%    filename        name of the datafile, including the .nxe extension
%    begsample       index of the first sample to read
%    endsample       index of the last sample to read
%    chanindx        index of channels to read (optional, default is all)
% This returns a Nchans X Nsamples data matrix

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

if nargin==1
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % read the header
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    hdr.Fs          = 1450;
    hdr.nChans      = 64;
    hdr.label  = cell(64,1);
    hdr.label(1:4)= {'GATE', 'TRIG1', 'TRIG2','EOG'};
    hdr.label(5:64) = {
        'Fp1'
        'Fpz'
        'Fp2'
        'AF1'
        'AFz'
        'AF2'
        'F7'
        'F3'
        'F1'
        'Fz'
        'F2'
        'F4'
        'F8'
        'FT9'
        'FT7'
        'FC5'
        'FC3'
        'FC1'
        'FCz'
        'FC2'
        'FC4'
        'FC6'
        'FT8'
        'FT10'
        'T7'
        'C5'
        'C3'
        'C1'
        'Cz'
        'C2'
        'C4'
        'C6'
        'T8'
        'TP9'
        'TP7'
        'CP5'
        'CP3'
        'CP1'
        'CPz'
        'CP2'
        'CP4'
        'CP6'
        'TP8'
        'TP10'
        'P9'
        'P7'
        'P3'
        'P1'
        'Pz'
        'P2'
        'P4'
        'P8'
        'P10'
        'PO3'
        'POz'
        'PO4'
        'O1'
        'Oz'
        'O2'
        'Iz'};

    % it is continuous data, therefore append all records in one trial
    hdr.nTrials     = 1;

    fid=fopen_or_error(filename,'r','l');
    fseek(fid,0,'eof');
    numBytes = ftell(fid);
    hdr.nSamples = (numBytes/2)/hdr.nChans;

    hdr.nSamplesPre = 0;

    fclose(fid);

    % return the header
    dat = hdr;

else
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % read the data
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    sfEEG = (1/2000) * (10/65535) * 1000000;
    sfEOG = (1/400) * (10/65535) * 1000000;
    sfTRIG = 10 * (10/65535);

    numChannels = 64;

    fid = fopen_or_error(filename,'r','l');
    fseek(fid, 2*numChannels*(begsample-1),'bof');
    data = fread(fid,[numChannels endsample-begsample+1],'short');
    fclose(fid);

    data(1:3,:) = sfTRIG.*data(1:3,:);
    data(4,:) = sfEOG.*data(4,:);
    data(5:64,:) = sfEEG.*data(5:64,:);

    if nargin<4
        chanindx = 1:numChannels;
    end

    dat = data(chanindx,:);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef math_model_opf_acps_legacy < mp.math_model_opf_acps & mp.mm_shared_opf_legacy
%MP.MATH_MODEL_OPF_ACPS_LEGACY  MATPOWER mathematical model for AC optimal power flow (OPF) problem.
%   ?
%
%   MP.MATH_MODEL_OPF_ACPS_LEGACY ... power flow ...
%
%   Properties
%       ? - ?
%
%   Methods
%       ?

%   MATPOWER
%   Copyright (c) 2021-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end

    methods
        %% constructor
        function obj = math_model_opf_acps_legacy()
            obj@mp.math_model_opf_acps();
            if nargin > 0 && isstruct(mpc)
                obj.mpc = mpc;
            end

            %% Due to a bug related to inheritance in constructors in
            %% Octave 5.2 and earlier (https://savannah.gnu.org/bugs/?52614),
            %% INIT_SET_TYPES() cannot be called directly in the
            %% MP_IDX_MANAGER constructor, as desired.
            %%
            %% WORKAROUND:  INIT_SET_TYPES() is called explicitly as needed
            %%              (if om.var is empty) in ADD_VAR(), DISPLAY() and
            %%              INIT_INDEXED_NAME(), after object construction,
            %%              but before object use.
        end

        function obj = add_named_set(obj, varargin)
            % call parent method (also checks for valid type for named set)
            add_named_set@mp.math_model_opf_acps(obj, varargin{:});
            obj.add_named_set_legacy(varargin{:});
        end

        function obj = def_set_types(obj)
            obj.def_set_types_legacy();
        end

        function obj = init_set_types(obj)
            init_set_types@mp.math_model_opf_acps(obj);
            obj.init_set_types_legacy();
        end

        function obj = build(obj, nm, dm, mpopt)
            obj.mpc = dm.source;
            build@mp.math_model_opf_acps(obj, nm, dm, mpopt);
            obj.build_legacy(nm, dm, mpopt);
        end

        function obj = add_vars(obj, nm, dm, mpopt)
            add_vars@mp.math_model_opf_acps(obj, nm, dm, mpopt);    %% call parent

            %% legacy user-defined variables
            if isfield(dm.userdata, 'legacy_opf_user_mods')
                obj.add_legacy_user_vars(nm, dm, mpopt);
            end
        end

        function add_system_costs(obj, nm, dm, mpopt)
            add_system_costs@mp.math_model_opf_acps(obj, nm, dm, mpopt);    %% call parent

            %% legacy user-defined costs
            if isfield(dm.userdata, 'legacy_opf_user_mods')
                obj.add_legacy_user_costs(nm, dm, 0);
            end
        end

        function obj = add_system_constraints(obj, nm, dm, mpopt)
            %% call parent
            add_system_constraints@mp.math_model_opf_acps(obj, nm, dm, mpopt);

            %% legacy user-defined constraints
            if isfield(dm.userdata, 'legacy_opf_user_mods')
                obj.add_legacy_user_constraints_ac(nm, dm, mpopt);
            end
        end

        function names = legacy_user_var_names(obj)
            names = {'Va', 'Vm', 'Pg', 'Qg'};
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%kvgamth 'Generate a Binary Image by Thresholding a Color Image (K1)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vgamth.pane file
%
% Parameters:
% InputFile: i 'Input Image ', required: 'input image'
% Integer: b 'High Threshold Map Column 1 (RED) ', default: 128: 'high threshold for map column 1'
% Integer: c 'Low Threshold Map Column 1 (RED) ', default: 128: 'low threshold for map column 1'
% Integer: d 'High Threshold Map Column 2 (GREEN)', default: 128: 'high threshold for map column 2'
% Integer: e 'Low Threshold Map Column 2 (GREEN) ', default: 128: 'low threshold for map column 2'
% Integer: f 'High Threshold Map Column 3 (BLUE) ', default: 128: 'high threshold for map column 3'
% Integer: g 'Low Threshold Map Column 3 (BLUE) ', default: 128: 'low threshold for map column 3'
% Integer: v 'Non-Zero Output Value ', default: 128: 'non-zero output value'
% OutputFile: o 'Output Image', required: 'resulting image'
%
% Example: o = kvgamth(i, {'i','';'b',128;'c',128;'d',128;'e',128;'f',128;'g',128;'v',128;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vgamth - Generate a Binary Image by Thresholding a Color Image  (K1)
%
%  DESCRIPTION
% .I vgamth
% generates a binary image by thresholding a pseudo color input image.
% The result image is of data type BYTE, and the input image must also
% be of type BYTE.
%
% The output pixel is given the specified value if:
% 1) only the lower value option is given and the pixel
% value is larger then the lower threshold,
% 2) only the upper value option is given and the pixel
% value is less than the upper threshold,
% 3) the upper value and lower value options are given and the pixel
% value is between the upper and lower threshold,
% or
% 4) same as 3 but both thresholds are equal and the pixel
% value is same as the threshold level.
% Otherwise, the output pixel is given a zero value.
%
% The input image must contain a map with three columns
% The thresholding is actually performed
% on the maps; thus the input image must have a map scheme of
% FORCE, and the map data type must be BYTE.
%
% By choosing the appropriate limits is it possible to identify
% areas in an image that are above or below a certain value, or to
% isolate certain bands of intensities in each color map. There are
% many different combinations that can be used to produce the results
% required.
%
% Note that this type of image is NOT a Khoros BIT image!
% The output image is of the data type BYTE, but the only
% pixel values are 0 and non-zero.
%
%
%
%  EXAMPLES
% vgamth -i input.image -o output.image -b 125 -c 75 -d 125 -e 75 -f 125 -g 75
% -v 100
%
% This command performs a threshold on input.image and stores the
% result in output.image.  All pixel values between 75 and 125 are
% given the value 100 in the output image.  All pixel values that
% do not fall within this range are give the value 0.
%
%  "SEE ALSO"
%
%  RESTRICTIONS
% .I vgamth
% will only operate on 1 band images of data storage type BYTE, with the
% map enable FORCE, and map data type BYTE.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
close all
clear all
path(path,'..\..\..\FUZZCLUST')
%the data
load motorcycle.txt
data.X = motorcycle(:,[1 2]);

%parameters
param.c=4;
param.m=2;
param.e=1e-4;
normalas = 1;
param.val=1;
%normalization
data=clust_normalize(data,'range');
%clustering
result = FCMclust(data,param);
param.c=result.data.f;
result = GGclust(data,param);
plot(data.X(:,1),data.X(:,2),'b.',result.cluster.v(:,1),result.cluster.v(:,2),'ro');
hold on
%draw contour-map
new.X=data.X;
eval=clusteval(new,result,param);
%validation
result = validity(result,data,param);
result.validity

author: Sable

repo: mcbench-benchmarks


text:
% GAMMA_LOGPDF - The log pdf of a Gamma distribution.
%
% Computes the logarithm of the Gamma probability density function
% G(ALPHA|A,B), where ALPHA is the random variable, A is the shape and B is
% the inverse scale (note the difference to MATLAB's built-in functions).
%
% The function is called as
%
% L = GAMMA_LOGPDF(ALPHA, A, B)
% L = GAMMA_LOGPDF(ALPHA, A, B, LOG_ALPHA)
%
% where LOG_ALPHA is LOG(ALPHA). Letting the user to evaluate this logarithm
% allows greater flexibility (e.g., for variational Bayesian cost function
% one can give the expectation of the logarithm). However, A and B are
% assumed to be fixed parameters in order to keep the function simple.

% Last modified 2010-06-07
% Copyright (c) Jaakko Luttinen

function [lp, dlp] = gamma_logpdf(alpha, a, b, log_alpha)

if nargin < 4
  log_alpha = log(alpha);
% $$$   if nargin == 4
% $$$     error(['Check the order of your input variables - the function has been ' ...
% $$$            'changed!'])
% $$$   end
% $$$ else
% $$$   error('Improper number of inputs')
end

lp = (a-1).*log_alpha - b.*alpha - gammaln(a) + a.*log(b);

if nargout >= 2
  inv_alpha = 1./alpha;
  dlp = (a-1).*inv_alpha - b;
end

%lp(alpha<=0) = -inf;
%dlp(alpha<=0) = nan;

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function MatingPool = MatingSelection(PopObj,div)
% The mating selection of GrEA

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    N = size(PopObj,1);

    %% Calculate the grid location of each solution
    fmax = max(PopObj,[],1);
    fmin = min(PopObj,[],1);
    lb   = fmin-(fmax-fmin)/2/div;
    ub   = fmax+(fmax-fmin)/2/div;
    d    = (ub-lb)/div;
    lb   = repmat(lb,N,1);
    d    = repmat(d,N,1);
    GLoc = floor((PopObj-lb)./d);
    GLoc(isnan(GLoc)) = 0;

    %% Calculate the GD value of each solution
    GD = zeros(N)+inf;
    for i = 1 : N-1
        for j = i+1 : N
            GD(i,j) = sum(abs(GLoc(i,:)-GLoc(j,:)));
            GD(j,i) = GD(i,j);
        end
    end

    %% Calculate the GCD value of each solution
    GD  = max(size(PopObj,2)-GD,0);
    GCD = sum(GD,2);

    %% Binary tournament selection
    Parents1   = randi(N,1,N);
    Parents2   = randi(N,1,N);
    Dominate   = any(PopObj(Parents1,:)<PopObj(Parents2,:),2) - any(PopObj(Parents1,:)>PopObj(Parents2,:),2);
    GDominate  = any(GLoc(Parents1,:)<GLoc(Parents2,:),2) - any(GLoc(Parents1,:)>GLoc(Parents2,:),2);
    MatingPool = [Parents1(Dominate==1 | GDominate==1),...
                  Parents2(Dominate==-1 | GDominate==-1),...
                  Parents1(Dominate==0 & GDominate==0 & GCD(Parents1)<=GCD(Parents2)),...
                  Parents2(Dominate==0 & GDominate==0 & GCD(Parents1)>GCD(Parents2))];
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function value = rhs_test01 ( x_num, x, t )

%*****************************************************************************80
%
%% RHS_TEST01 evaluates the right hand side for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X_NUM, the number of nodes.
%
%    Input, real X(X_NUM,1), the node coordinates.
%
%    Input, real T, the current time.
%
%    Output, real VALUE(X_NUM,1), the source term.
%
  value(1:x_num,1) = 0.0;

  return
end


text:
function kappa = CalcCurvature(v, vtheta, indx, dt)

%CALCCURVATURE
%   Usage:  kappa = CalcCurvature(v,theta,indx,dt)
%
% This function calculates the curvature of the path that the fly took
% while walking in the frames indexed by indx.  The inputs are the speed v
% , the angle of the velocity vector vtheta,the index of walking frames,and the time resolution
% (1/sampling rate) dt. I am convinced that this function can be
% simplified, but haven't had the time or energy to find a more efficient way.

%Written by Dan Valente
%September 2007

list = [];
phi = [];
kappa = [];
indx(end+1)=-1;
for i = 1:length(indx)-1
    if (indx(i)+1 == indx(i+1))
        list = [list indx(i)];
    else
        list = [list indx(i)];
        temp = vtheta(list)';
        if (length(temp)~= 1)
            for j = 1:length(temp)-1
                d = temp(j+1)-temp(j);
                if abs(d) >= pi
                    if (temp(j+1) < 0)
                        d = (temp(j+1)+2*pi)-temp(j);
                    elseif (temp(j+1) >= 0)
                        d = (temp(j+1)-2*pi)-temp(j);
                    end
                end
                phi(j) = d/dt;
            end

            phi2 = phi./v(list(1:end-1))'; %transpose for FAnalyze;
            kappa = [kappa phi2];
        end
        list = [];
        phi = [];
    end
end

return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function W = randInitializeWeights(L_in, L_out)
%RANDINITIALIZEWEIGHTS Randomly initialize the weights of a layer with L_in
%incoming connections and L_out outgoing connections
%   W = RANDINITIALIZEWEIGHTS(L_in, L_out) randomly initializes the weights
%   of a layer with L_in incoming connections and L_out outgoing
%   connections.
%
%   Note that W should be set to a matrix of size(L_out, 1 + L_in) as
%   the first row of W handles the "bias" terms
%

% You need to return the following variables correctly
W = zeros(L_out, 1 + L_in);

% ====================== YOUR CODE HERE ======================
% Instructions: Initialize W randomly so that we break the symmetry while
%               training the neural network.
%
% Note: The first row of W corresponds to the parameters for the bias units
%

% Randomly initialize the weights to small values
epsilon_init = 0.12;
W = rand(L_out, 1+L_in)*2*epsilon_init - epsilon_init;

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

end

author: yhyap

repo: machine-learning-coursera


text:
function [elist,nextfront]=advancefront(edges,loop,elen)
%
% [elist,nextfront]=advancefront(edges,loop,elen)
%
% advance an edge-front on an oriented surface to the next separated by
% one-element width
%
% author: Qianqian Fang, <q.fang at neu.edu>
% date: 2012/02/09
%
% input:
%      edges: edge list of an oriented surface mesh, must be in CCW order
%      loop: a 2-column array, specifying a closed loop in CCW order
%      elen: node number inside each element, if ignored, elen is set to 3
%
% output:
%      elist: list of triangles that is enclosed between the two
%             edge-fronts
%      nextfront: a new edge loop list representing the next edge-front
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

elist=[];
nextfront=[];

if(nargin<3) elen=3; end

[hasedge, loc]=ismember(loop,edges,'rows');

if(~all(hasedge))
   error('loop edge is not defined in the mesh');
end

nodenum=size(edges,1)/elen;

elist=unique(mod(loc-1,nodenum))+1;
nextfront=edges(elist,:);
for i=1:elen-1
    nextfront=[nextfront;edges(elist+nodenum*i,:)];
end
nextfront=setdiff(nextfront,loop,'rows');

% remove reversed edge pairs
[flag,loc]=ismember(nextfront,nextfront(:,[2 1]),'rows');
id=find(flag);
if(~isempty(id))
    delmark=flag;
    delmark(loc(find(loc>0)))=1;
    nextfront(find(delmark),:)=[];
end
nextfront=nextfront(:,[2 1]); % reverse this loop, as it is reversed to the input loop


text:
% This file is part of the project NILM-Eval (https://github.com/beckel/nilm-eval).
% Licence: GPL 2.0 (http://www.gnu.org/licenses/gpl-2.0.html)
% Copyright: ETH Zurich, 2014
% Author: Romano Cicchetti

% parameters
dataset = 'thun';
granularity = 60;
houses = [5];
appliance_name = 'microwave';
dates = ['2012-07-11'];
start_date = '2012-07-11';
end_date = '2012-07-11';
consumptionArray = zeros((24*60*60)/granularity,5);

% obtain consumption data
appliance_id =  getApplianceID(appliance_name);
numOfDaysPerHousehold = zeros(length(houses), 1);
counter = 0;
for house = houses
    consumption = read_plug_data(dataset, house, appliance_id, dates, granularity );
    counter = counter + 1;
    consumptionArray(1:length(consumption),counter) = consumption;
end

% plot consumption data
f = figure();
startDate = datenum(strcat(start_date, ' 00:00:00'));
endDate = datenum(strcat(end_date, ' 23:59:00'));
timeTics = linspace(startDate,endDate,86400/granularity);
plot(timeTics, consumptionArray);
datetickzoom('x', 'HH:MM');
ylabel('Power [W]');
xlabel('Time of day');
%legend('Household 1', 'Household 2', 'Household 4', 'Household 5', 'Household 6');
grid on;
%title(sprintf('Power Consumption of %s', appliance_name));
axis tight;

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
%%
% Simple example of gradient flow in metric space.

rep = '../results/gradflow-metric/';
[~,~] = mkdir(rep);

f = @(x,y)x.^2 + y.^2;
Gx = @(x,y)2*x;
Gy = @(x,y)2*y;

% background image
N = 256;
tx = linspace(-.05,1,N);
ty = linspace(-.08,.5,N);
[Y,X] = meshgrid(ty,tx);
F = f(X,Y);
clf; hold on;
imagesc(tx,ty,F');
contour(tx,ty,F', 15, 'k');

% for the implicit stepping
N1 = 2048;
tx1 = linspace(min(tx),max(tx),N1);
ty1 = linspace(min(ty),max(ty),N1);
[Y1,X1] = meshgrid(ty1,tx1);

plist = [1 1];
plist = [1 1.3 1.5 1.8 2 3 5 10 Inf];

gmode = 'explicit';
gmode = 'implicit';

tau = .25; niter = 70;

lgd = {};
for ip=1:length(plist)
    p=plist(ip);

    % conjugate exponent
    q = conjexp(p);
    ri = @(r,p)sign(r).*abs(r).^(1/(p-1));

    grad = @(gx,gy)ri([gx, gy],p);
    if p==1
        grad = @(gx,gy)[gx.*(gx==max(gx,gy)), gy.*(gy==max(gx,gy))];
    elseif p==Inf
        grad = @(gx,gy)[sign(gx),sign(gy)];
    end
    grad = @(gx,gy)grad(gx,gy)*norm([gx,gy],q)/norm(grad(gx,gy),q);

    % initial point
    x = .95; y = .45;
    U = []; V = [];
    % do the descent
    for i=1:niter
        U(end+1) = x; V(end+1) = y;
        switch gmode
            case 'explicit'
                g = grad(Gx(x,y), Gy(x,y));
                x = x - tau*g(1);
                y = y - tau*g(2);
            case 'implicit'
                E = ( abs(X1-x).^p + abs(Y1-y).^p ).^(2/p) + tau * ( X1.^2 + Y1.^2 );
                if p==Inf
                    E = max(abs(X1-x),abs(Y1-y)).^2 + tau * ( X1.^2 + Y1.^2 );
                end
                [tmp,i] = min(E(:));
                x = X1(i); y = Y1(i);
        end
    end

    lgd{end+1} = ['p=' num2str(p)];
    c = (ip-1)/(length(plist)-1);
    plot(U,V, '.-', 'Color', [c 0 1-c], 'LineWidth', 2, 'MarkerSize', 23);
    colormap parula(256);
    axis image; axis off;
    axis([min(tx) max(tx) min(ty) max(ty)]);
    drawnow;

end

saveas(gcf, [rep 'gradflow-' gmode '.eps'], 'epsc');

author: optimaltransport


text:
function x = icholpre(b,A,L,R)

%% Pre-smoothing by Gauss-Seidel
x = tril(A)\b;
r = b - A*x;

%%
e = L\r;
e = R\e;      % solve by achol decomposition L*L'
x = x + e;

%% Post-smoothing
x = x + triu(A)\(b-A*x);

repo: ifem


text:
clear
clc

addpath(genpath(pwd));
addpath('../../edges-master');
deep_model_path = 'model_deep_contour/';
deepModel.model_del_file = [deep_model_path 'deploy_cov4.prototxt'];
deepModel.model_file = [deep_model_path 'deep_contour_model'];
deepModel.mean_proto_file = [deep_model_path 'mean.binaryproto.txt'];
deepModel.use_gpu = 1;
modelDir = [deep_model_path 'models'];
use_gpu = deepModel.use_gpu;
gpu_id = 0;
[nBatch, nC, patchSize, ~, nOutputNum] = parse_del_file(deepModel.model_del_file);

PATCH_MEAN = single(dlmread(deepModel.mean_proto_file));

PATCH_MEAN = single(reshape(PATCH_MEAN, patchSize, patchSize, nC));

if exist('use_gpu', 'var')
  matcaffe_init(deepModel.use_gpu, deepModel.model_del_file, deepModel.model_file, gpu_id);
else
  matcaffe_init();
end
featSelecttion = false;

if featSelecttion
    if exist([modelDir '/feat_selection.mat'], 'file')
        load([modelDir '/feat_selection.mat']);
    else
        [selected_dims, ndim] = deepSparseFeatSelection([modelDir '/feat/'], 1e-3, 0.7);
        if(~isempty(ndim) && ~isempty(selected_dims))
            save([modelDir '/feat_selection.mat'], 'selected_dims', 'ndim');
        end
    end
else
    selected_dims = [];
end

dlPara.selected_dims = selected_dims;
dlPara.patch_mean = PATCH_MEAN;
dlPara.modelDir = modelDir;

edgesDLDemo(dlPara);

repo: DeepContour


text:
classdef (ConstructOnLoad) CoreEventData < event.EventData
    % CoreEventData. Encapsulate data to be passed when an event occurs.
    %
    %     Use this class when you trigger an event which needs to provide data
    %     to its listeners.
    %
    %     Example
    %     -------
    %
    %     This example triggers a MouseClick event and passes the image coordiates
    %
    %         coords = GetImageCoordinates;
    %         notify(obj, 'MouseClick', CoreEventData(coords));
    %
    %
    %     Licence
    %     -------
    %     Part of CoreMat. https://github.com/tomdoel/coremat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        Data
    end

    methods
        function obj = CoreEventData(data)
            obj.Data = data;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%MATCHGMS  GMS (Grid-based Motion Statistics) feature matching strategy
%
%     matchesGMS = cv.matchGMS(size1, keypoints1, size2, keypoints2, matches1to2)
%     matchesGMS = cv.matchGMS(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __size1__ Size of first image `[w,h]`.
% * __keypoints1__ Keypoints from the first source image. A 1-by-N structure
%   array with the following fields:
%   * __pt__ coordinates of the keypoint `[x,y]`
%   * __size__ diameter of the meaningful keypoint neighborhood
%   * __angle__ computed orientation of the keypoint (-1 if not applicable).
%     Its possible values are in a range [0,360) degrees. It is measured
%     relative to image coordinate system (y-axis is directed downward), i.e
%     in clockwise.
%   * __response__ the response by which the most strong keypoints have been
%     selected. Can be used for further sorting or subsampling.
%   * __octave__ octave (pyramid layer) from which the keypoint has been
%     extracted.
%   * **class_id** object id that can be used to clustered keypoints by an
%     object they belong to.
% * __size2__ Size of second image `[w,h]`.
% * __keypoints2__ Keypoints from the second source image. Same format as
%   `keypoints1`.
% * __matches1to2__ Input 1-nearest neighbor matches from the first image to
%   the second one. A 1-by-M structure array with the following fields:
%   * __queryIdx__ query descriptor index (zero-based index)
%   * __trainIdx__ train descriptor index (zero-based index)
%   * __imgIdx__ train image index (zero-based index)
%   * __distance__ distance between descriptors (scalar)
%
% ## Output
% * __matchesGMS__ Matches returned by the GMS matching strategy.
%
% ## Options
% * __WithRotation__ Take rotation transformation into account. default false
% * __WithScale__ Take scale transformation into account. default false
% * __ThresholdFactor__ The higher, the less matches. default 6.0
%
% GMS feature matching strategy by [Bian2017gms].
%
% ## Notes
% * Since GMS works well when the number of features is large, we recommend to
%   use cv.ORB features and set `FastThreshold` to 0 to get as many features
%   as possible quickly.
% * If matching results are not satisfying, please add more features (we use
%   10000 for images with 640x480 size).
% * If your images have big rotation and scale changes, please set
%   `WithRotation` or `WithScale` to true.
%
% ## References
% [Bian2017gms]:
% > JiaWang Bian, Wen-Yan Lin, Yasuyuki Matsushita, Sai-Kit Yeung, Tan Dat
% > Nguyen, and Ming-Ming Cheng. "GMS: Grid-based motion statistics for fast,
% > ultra-robust feature correspondence".
% > In IEEE Conference on Computer Vision and Pattern Recognition, 2017.
%
% See also: cv.FeatureDetector, cv.DescriptorMatcher, cv.drawMatches
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [Oq,Oqval,DEOPA,goodperiods] = OP(simppeak,SIG,dSIG,valid,doublepeaks,COEF,resampC,Ts,method)
% OP: Opening Peaks detection.
% Input files:
% - <simppeak>: contains a list of closing-peak times in ms, and amplitudes
% (<simppeak> is a 2-column matrix)
% - <dSIG>: Derivative of EGG signal
% - <valid> : vector containing information on peaks that have several summits
% - <doublepeaks>: vector containing information on peaks that cluster together
% - <COEF>: contains, among other coefficients, the smoothing step chosen by the
% user. It is not used in analyzing the positive part of the signal, only in the
% negative part, for the calculation of Oq.
% - <resampC>: coefficient for reinterpolation
% - <Ts>: the inverse of the sampling rate
% - <method>: indication on the method which is to be used in the handling of
% double / multiple peaks: selecting either the first or the highest.
% (The fact that there is not one single peak will also be passed on
% to the following functions, so this choice is not of much consequence.)
%
% Output files:
% <Oq>: vector containing open quotient values detected using the maxima
% <Oqval>: open quotient values computed by peak detection (single peaks only)
% <DEOPA>: vector containing DEOPA values
% <SM>, for Success of Method: variable indicating the number of Oq values that
% could be calculated using peak detection.
% <goodperiods>: beginning and end of cycles for which closing peaks are unique.
% <validGOI>: ratio of highest to second highest peak.
%
% No use of autocorrelation or comparison across cycles, otherwise the method
% would not give results up to the end of voicing: border phenomena.
% Passages where voice quality changes quickly are not found only at the offset
% of voicing: they can be word-internal, as in the case of consonantal glottal
% stops, ejectives, implosives, or glottalized tones, e.g. in Vietnamese tone
% C2, which has medial glottal constriction.
%
% Method: detecting maxima; checking the plausibility of results; if results are
% implausible, detecting maxima anew with increased smoothing step

% Setting the output vector Oqval at []; otherwise, in case no single period has
% single closing and opening peaks, this output argument is not assigned during
% call to OP.
Oqval = [];

% creating a matrix containing the beginning and end of periods for which Oq
% will be calculated. Values are converted to indices.
goodperiods = [];

if ~isempty(simppeak)
    if length(simppeak(:,1)) > 1
        for i = 1: length(simppeak(:,1)) - 1
            goodperiods(i,1) = round(simppeak(i,1) * (1 / Ts));
            goodperiods(i,2) = round(simppeak(i + 1,1) * (1 / Ts));
        end

        % Previous calculation method: eliminating periods for which closing peaks are
        % not unique:
        % %%%%%%%%%%%%%%%%%%%%%% Determining periods for which closings are unique
        % % Open quotient measurements are attempted only when closings are unique. For
        % % each glottal cycle, there are two condition:
        % % 1) that the values associated to the peak in <valid> be
        % % above 0.6, i.e. within the peak there is no other summit that is higher than
        % % 60% of the highest summit
        % % 2) that none of the neighbouring peaks be higher than 60% of the candidate peak.
        %
        % goodperiods = [];
        % for i = 1:length(valid) - 1
        %     if (doublepeaks(i) == 0) & (valid(i) < 0.6) & (doublepeaks(i + 1) == 0) & (valid(i) < 0.6)
        %         goodperiodsnb = goodperiodsnb + 1;
        %         goodperiods(goodperiodsnb,1) = Tgci(i,1);
        %         goodperiods(goodperiodsnb,2) = Tgci(i + 1,1);
        %     end
        % end

        %%%%%%%%%%%%%%%%%%%%%% First measurement of openings: done by detecting
        %%%%%%%%%%%%%%%%%%%%%% local minima, without any condition on the values obtained.
        texis = exist('goodperiods');
        if texis == 1
            if length(goodperiods(:,1)) > 0
                for i = 1:length(goodperiods(:,1))
                    % A rounding of the indices is necessary: as there has been
                    % reinterpolation, the values are not integers.
                    [DEOPA(i),GOI(i)] = min(    dSIG( round(goodperiods(i,1)):round(goodperiods(i,2) ))   );
                    % correction: adding the index corresponding to the first point of the
                    % excerpt
                    GOI(i) = GOI(i) + goodperiods(i,1);
                    %% Calculation of Oq, in %
                    Oq(i) = 100 * ((  goodperiods(i,2) - GOI(i)  ) / (  goodperiods(i,2) - goodperiods(i,1)  ));
                end
        %     else
        %         DEOPA = 0;
        %         GOI = 0;
        %         Oq = 0;
            end
        end

        %%%%%%%%%%%%%%%%%%%%%%%% Measurement by peak detection and barycentre calculation
        texis = exist('goodperiods');
        if texis == 1
            if length(goodperiods(:,1)) > 0
                for i = 1:length(goodperiods(:,1))
                    % if you want to see information displayed as analysis is being
                    % conducted, uncomment the two lines below:
                    %% disp(['Treating period number ',num2str(i)])
                    %% disp(goodperiods(i,:))
                    % placing the indices of beginning and end of the relevant portion of
                    % the signal in variables
                    % A rounding of the indices is necessary: as there has been
                    % reinterpolation, the values are not integers. In the present state of
                    % the programme, this rounding is performed earlier, when the
                    % <goodperiods> matrix is built up.
                    BE = goodperiods(i,1);           % BEginning of period
                    EN = goodperiods(i,2);           % ENd of period
                    % placing the extract in a vector, inverting it because CRO operates on
                    % a positive signal. An extra precaution must be taken: excluding the
                    % first and last eight points of the signal (by setting them at zero):
                    % it is physically impossible that a clear opening will take place
                    % at these points, so close to the closing; technically, a small part of
                    % the signal immediately before the closing peak that marks the
                    % beginning of the period (res. after the closing peak that marks its
                    % end) may be included and result in program crash.
                    dSIGextr = - dSIG(BE:EN);
                    L = length(dSIGextr);
                    if L > 16
                        dSIGextr(1:8) = 0;
                        dSIGextr(L-7:L) = 0;
                    end
                    % Also creating an extract for the EGG signal, which is used later on,
                    % by AMPOS
                    SIGextr = - SIG(BE:EN);

                    %%%%%%%%%%%%%%%%%% Calling the function CRO to detect crossings.
                    % The <method> toggle (3rd value in <COEF>) must be set at 0, and the
                    % amplitude threshold (4th value in <COEF>) passed on to the function.
                    COEFOPEN = COEF;
                    COEFOPEN(3) = 0;

                    % setting the threshold for peak detection (empirically). Choice in
                    % Henrich 2001:123 : setting it at 70% of the highest point.
                    COEFOPEN(4) = 0.70 * max(dSIGextr);

                    % finding out the crossings with the DEGG signal. In the best cases, it is
                    % expected that there will be only two crossings, one up and one down.
                    % But this is not frequent: the signal is often dented.
                    [rimsOPEN] = CRO(dSIGextr,COEFOPEN);

                    if isempty(rimsOPEN)
                        % If there is no single detected threshold crossing: the user should
                        % changed the maximum F0 value and ask for the results to be
                        % calculated anew. The user is advised to do so inside the <rim>
                        % function.
                        Oqval(i) = 0;
                        DEOPA(i) = 0;
                    else
                        % finding out the amplitude and position of the peaks. No
                        % reinterpolation of the signal is useful, as the signal-to-noise ratio in
                        % the negative part of the DEGG signal is poor. Quite the opposite: the
                        % signal is smoothed (over 3 samples, i.e. 0.07 ms).

                        % The threshold for exclusion of small peaks is set at 0.7 in the case of
                        % opening peaks, following Henrich 2001:123. This has in fact already
                        % been applied before, in detection of peak rims by RIM; but the
                        % function PEAKSHAPE called by AMPOS needs this input.
                        propthresh = 0.7;
                       % figure(1)
                       % clf
                       % plot(SIGextr)
                        [TGOI,TGOIFo,validGOI,validtimeGOI] = AMPOS(SIGextr,rimsOPEN,1,Ts,method,propthresh);

                        % detection of double / multiple closing peaks ("peak clusters"). The
                        % threshold passed to this function is: (the inverse of) one-fifth of
                        % the period, under the assumption that peaks separated by more than
                        % one-fifth of the period should not be considered as a bundle. In these
                        % cases, no Oq is calculated.
                        maxF = 1 / ( (goodperiods(i,2) - goodperiods(i,1)) / 5);
                        doublepeaks = detectmult(TGOI,maxF);

                        % choice of opening peak, leaving only one, following the method
                        % specified by the user. If there are 2 peaks that are too far apart for
                        % averaging to make sense, there will be 2 lines in simppeak: two "real
                        % peaks" according to the criteria used in the SIMP function.
                        [simppeak,BUND] = simp(TGOI,doublepeaks,method);

                        % calculating the open quotient when possible. The condition is: if
                        % there is one single value in <simppeak>, the matrix where detected
                        % peaks are stored. If there are two (or more) values, this indicates
                        % that there were two (or more) peaks that were too far apart for
                        % averaging to make sense, and it seems best not to give an Fo value;
                        % the Oq value is then set at zero. This is also what happens if it
                        % remains unaffected: it comes out as a zero inside the vector; but if
                        % the last values are zero, then the length of the vector is less than
                        % the size of the Fo vector, resulting in assignment problems. So the
                        % value is set at zero in this programme.
                        % inside the vector).
                        if length(simppeak(:,1)) == 1
                            % transforming the time of peak to an index
                            indexsimppeak = simppeak(1,1) * (1 / Ts);
                            Oqval(i) = 100 * (       ((goodperiods(i,2) - goodperiods(i,1)) - indexsimppeak ) / ...
                                  (goodperiods(i,2) - goodperiods(i,1))       );
                        else
                              Oqval(i) = 0;
                        end
                    % end of "if" loop: condition on existence of rimsOPEN values.
                    end
                end
            end
        end
    else
        goodperiods = []; Oq = []; Oqval = []; DEOPA = [];
    end
else
    goodperiods = []; Oq = []; Oqval = []; DEOPA = [];
end
        %%%%%%%%%%% The method for detecting opening peaks is fairly similar to
        %%%%%%%%%%% that for detecting closing peaks. Concerning the detection
        %%%%%%%%%%% of multiple peaks, however, in the case of closing peaks
        %%%%%%%%%%% the threshold for detection was set rather low, and the
        %%%%%%%%%%% peaks were later sorted according to size and closeness to
        %%%%%%%%%%% one another. In the case of opening peaks, the period is
        %%%%%%%%%%% known already, so the threshold is fixed straight away,
        %%%%%%%%%%% on the basis of the highest point in this part of the
        %%%%%%%%%%% signal: if the second highest peak is 70% as high as the
        %%%%%%%%%%% highest peak, it is decided that the peak is not unique
        %%%%%%%%%%% (following Henrich 2001:123).

        %%%%%%%%%%% It is an issue whether this threshold must apply
        %%%%%%%%%%% within the peak region as well as in the comparison of
        %%%%%%%%%%% neighbouring peaks. The fine detail of the peak is often
        %%%%%%%%%%% slightly dented (see figures); this will result in
        %%%%%%%%%%% "secondary peaks" being detected very close to the maximum
        %%%%%%%%%%% (about 0.1 ms to its left and right). A time condition could
        %%%%%%%%%%% be added: if the secondary peak is very close
        %%%%%%%%%%% (threshold empirically set at a time value, e.g. 0.5 ms, or
        %%%%%%%%%%% as a proportion of the period), it would not count as
        %%%%%%%%%%% a doubled peak. In the present programme, the solution
        %%%%%%%%%%% chosen instead is to use smoothing of the signal.

%         % In an earlier version, the algorithm was:
%         % 1) retrieve highest peak in TGOI; 2) check that it is unique
%         % internally; 3) check that the neighbouring peaks are less than 70% of its
%         % height. GOIC stands for Glottis-Opening-Instant Candidate.
%         % First step:
%         [GOICvalue,GOICindex] = max (TGOI(:,2));
%         % Second step: finding out whether the peak is unique internally. For
%         % this: look up the corresponding slot in vector <validtimeGOI>.
%         if validtimeGOI(GOICindex) < threshtime
%             % retrieving the index of maximum value in TGOI relative to the
%             % beginning of the <dSIGextr> extract: it is the value
%             % in the first column of <TGOI> at the line indicated in
%             % <GOICindex>
%             GOICindex = TGOI(GOICindex,1);
%             % changing the relative index <GOICindex> into an index relative
%             % to the <SIG> signal as a
%             % whole, for computation of the open quotient
%             GOICindex = GOICindex + goodperiods(i,1);
%             % calculation of open quotient
%             Oqval(i) = 100 * (       (goodperiods(i,2) - GOICindex ) / ...
%                 (goodperiods(i,2) - goodperiods(i,1))       );
%         end

            %%%%%%%%%%% Note:
            %% In the first version of the programme, the threshold for detection
            %% of peak regions was first set at 40% of the signal, and the
            %% condition on peak amplitude was added later, by recovering the
            % second highest peak in TGOI and
            % comparing its amplitude to that of the highest peak.
            %% In the present version of the function, only the peaks above 70% of the maximum
            %% within the dSIG extract <dSIGextr> are detected, so that
            %% this part of the function is not necessary anymore.
            %% The calculations were as follows:

%             COEFOPEN(4) = 0.7 * max(dSIGextr);
%             % Setting value for TGOI(GOICindex,2) at zero and re-calculating the
%             % maximum:
%             TGOI(GOICindex,2) = 0;
%             [GOICvalue2,GOICindex2] = max (TGOI(:,2));
%             % condition on whether the proportion of the two peaks is below
%             % threshold
%             if GOICvalue2 / GOICvalue < threshGOI
%                 % retrieving the index of maximum value in TGOI: it is the value
%                 % in the first column of <TGOI> at the line indicated in
%                 % <GOICindex>
%                 GOICindex = TGOI(GOICindex,1);
%                 % changing the relative index <GOICindex> into an index relative
%                 % to the <SIG> signal as a
%                 % whole, for computation of the open quotient
%                 GOICindex = GOICindex + goodperiods(i,1);
%                 % calculation of open quotient
%                 Oqval(i) = 100 * (       (goodperiods(i,2) - GOICindex ) / ...
%                     (goodperiods(i,2) - goodperiods(i,1))       );
%             end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
%%% [2d_matrix, int, int, int] = convert_video3d_to_2d(3d_matrix)
%
function [M,m,n,p] = convert_video3d_to_2d(V)
  [m,n,p] = size(V);

  if(isa(V,'uint8'))
    M = uint8(zeros(m*n,p));
  else
    M = zeros(m*n,p);
  end

  for i = 1:p
    M(:,i) = reshape(V(:,:,i),[],1);
  end
end

author: andrewssobral

repo: lrslibrary


text:
close all; clear all; clc;

audio = audioload();

file = 'text.txt';
fid  = fopen(file, 'r');
text = fread(fid,'*char')';
fclose(fid);

msg = echo_dec_ts(audio.data);
err = BER(text,msg);
nc  = NC(text, msg);

fprintf('Text: %s\n', msg);
fprintf('BER : %d\n', err);
fprintf('NC  : %d\n', nc);

author: ktekeli

repo: audio-steganography-algorithms


text:
function [y, x] = tapas_dcm_euler_gen(DCM, Ep)
% [y, x] = tapas_dcm_euler_gen(DCM, Ep)
%
% Generates synthetic fMRI data under a given signal to noise ratio (SNR)
% with the fixed hemodynamic convolution kernel
%
%   Input:
%   	DCM         - model structure
%       Ep          - data-generating parameters
%
%   Output:
%       y           - generated BOLD signal
%       x           - generated neuronal signal
%

% ----------------------------------------------------------------------
%
% Authors: Stefan Fraessle (stefanf@biomed.ee.ethz.ch), Ekaterina I. Lomakina
%
% Copyright (C) 2016-2018 Translational Neuromodeling Unit
%                         Institute for Biomedical Engineering
%                         University of Zurich & ETH Zurich
%
% This file is part of the TAPAS rDCM Toolbox, which is released under the
% terms of the GNU General Public License (GPL), version 3.0 or later. You
% can redistribute and/or modify the code under the terms of the GPL. For
% further see COPYING or <http://www.gnu.org/licenses/>.
%
% Please note that this toolbox is in an early stage of development. Changes
% are likely to occur in future releases.
%
% ----------------------------------------------------------------------

% number of regions
nr = size(Ep.A,1);
A = Ep.A';

% Transpose each nrxnr matrix
mArrayB = permute(Ep.B, [2 1 3]);

% create an empty D-matrix
if (sum(Ep.D(:))==0)
    Ep.D = zeros(nr,nr,nr);
end

% Transpose each nrxnr matrix
mArrayD = permute(Ep.D, [2 1 3]);
if (isempty(mArrayD))
    mArrayD = zeros(nr,nr,nr);
end

% driving inputs
C = DCM.U.u*(Ep.C'/16);
U = full(DCM.U.u);

% hemodynamic constants
H = [0.64 0.32 2.00 0.32 0.32];

% constants for hemodynamic model
oxygenExtractionFraction = 0.32*ones(nr, 1);
alphainv                 = 1/H(4);
tau                      = H(3)*exp(Ep.transit);
gamma                    = H(2);
kappa                    = H(1)*exp(Ep.decay);
epsilon                  = 1.0*exp(Ep.epsilon)*ones(nr, 1);

% parameter list
paramList = [DCM.U.dt size(U,1) nr size(U,2) 0 1 1];

% neuronal signal and time courses for hemodynamic parameters
[x,~,~,v,q] = dcm_euler_integration(A,C,U,mArrayB,mArrayD,...
                   oxygenExtractionFraction,alphainv,tau,gamma,kappa,paramList);

% constants for BOLD signal equation
relaxationRateSlope      = 25;
frequencyOffset          = 40.3;
oxygenExtractionFraction = 0.4*ones(1,nr);
echoTime                 = 0.04;
restingVenousVolume      = 4;

% coefficients of BOLD signal equation
coefficientK1  = 4.3*frequencyOffset*echoTime*oxygenExtractionFraction;
coefficientK2  = epsilon'.*(relaxationRateSlope*oxygenExtractionFraction*echoTime);
coefficientK3  = 1 - epsilon';

% get the Euler indices
Indices = DCM.M.idx;

% BOLD signal time course
y = restingVenousVolume*( bsxfun(@times,coefficientK1,(1 - (q(Indices,:)))) +...
            bsxfun(@times,coefficientK2,(1 - (q(Indices,:)./v(Indices,:)))) +...
            bsxfun(@times,coefficientK3,(1-v(Indices,:))));

end

author: translationalneuromodeling

repo: tapas


text:
function res = verification( des, varargin)
%VERIFICATION Run the patch verification task
%  RES = VERIFICATION(DESCRIPTOR) Compute the results for descriptor
%  DESCRIPTOR. Descriptor can either be an object of `desc.memdesc` or a
%  string, in which case a descriptor from
%
%     `<HBROOT>/data/descriptors/DESCRIPTOR`
%
%  are loaded.
%  Stores the results in csv files in:
%
%      `<HBROOT>/matlab/results/default/DESCRIPTOR/verification.csv`
%
%  If the scores file already exists, skips the computation. Use
%  `'override', true` to change this behaviour.
%
%  Additionally accepts the following 'OptionName', OptionValue arguments:
%
%  'scoresroot' :: '<HBROOT>/matlab/results/default/'
%   Change for a different target path for scores.
%
%  'taskpath' :: '<HBROOT>/tasks'
%   Location of te tasks files.
%
%  'geom_noise' :: {'easy', 'hard', 'tough'}
%  Limit the geometry noise used for positives.
%
%  'negs' :: {'inter', 'intra'}
%  Used negative sets.
%
%  'split' :: {'a', 'b', 'c', 'illum', 'view', 'full'}
%  Limit the splits on which the scores are computed.
%
%  'override' :: false
%   If true, do override existing score file.

% Copyright (C) 2017 Karel Lenc
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).
if ischar(des), [des, varargin] = desc.memdesc(des, varargin{:}); end;
opts.scoresroot = fullfile(hb_path, 'matlab', 'scores', 'default');
[opts, varargin] = vl_argparse(opts, varargin);
opts.taskspath = fullfile(hb_path, 'tasks');
opts.scorespath = fullfile(opts.scoresroot, des.name, 'verification.csv');
opts.split = {'a', 'b', 'c', 'illum', 'view', 'full'};
opts.negs = {'inter', 'intra'};
opts.geom_noise = {'easy', 'hard', 'tough'};
opts.methods = struct('name', {'balanced', 'imbalanced'}, ...
  'args', {{}, {'posneg_ratio', 0.2}});
opts.override = false;
opts.verbose = false;
[opts, ~] = vl_argparse(opts, varargin);

if opts.verbose, display(opts); end;
if ~iscell(opts.split), opts.split = {opts.split}; end;
if ~iscell(opts.negs), opts.negs = {opts.negs}; end;
if ~iscell(opts.geom_noise), opts.geom_noise  = {opts.geom_noise}; end;
vl_xmkdir(fileparts(opts.scorespath));
if ~opts.override && exist(opts.scorespath, 'file')
  res = readtable(opts.scorespath);
  fprintf('Results loaded from %s.\n', opts.scorespath);
  return;
end;

pospairs = @(split) fullfile(opts.taskspath, ...
  sprintf('verif_pos_split-%s.csv', split));
negpairs = @(split, neg) fullfile(opts.taskspath, ...
  sprintf('verif_neg_%s_split-%s.csv', neg, split));

numtasks = numel(opts.split)*numel(opts.negs)*numel(opts.geom_noise);
status = utls.textprogressbar(numtasks, 'updatestep', 1, ...
  'startmsg', sprintf('Evaluating "verification" for %s ', des.name));
stepi = 1;
res = {};
descs_n = des;
for gni = 1:numel(opts.geom_noise)
  gnoise = opts.geom_noise{gni};
  for si = 1:numel(opts.split)
    split = opts.split{si};
    pos_dists = bench.verification.compute(pospairs(split), descs_n, gnoise);
    for ni = 1:numel(opts.negs)
      neg = opts.negs{ni};
      neg_dists = bench.verification.compute(negpairs(split, neg), descs_n, gnoise);
      for mi = 1:numel(opts.methods)
        method = opts.methods(mi);
        res_s = bench.verification.eval(pos_dists, neg_dists, method.args{:});
        res{end+1} = struct(...
          'descriptor', descs_n.name, 'split', split, 'negs', neg, ...
          'geom_noise', gnoise, 'method', method.name, ...
          'pr_auc', res_s.pr_auc, 'pr_ap', res_s.pr_ap, 'roc_auc', res_s.roc_auc);
      end
      status(stepi); stepi = stepi+1;
    end
  end
end
res = struct2table(cell2mat(res), 'AsArray', true);
writetable(res, opts.scorespath);
fprintf('\n');

author: hpatches

repo: hpatches-benchmark


text:
function [xi,w]=quadraturePoints1D(n,algorithm,c1)
%%QUADRATUREPOINTS1D Obtain quadrature points and weights to efficiently
%           numerically evaluate 1D integrals involving various weighting
%           functions. The quadrature points and weights are based off
%           properties of orthogonal polynomials. The evaluation of a
%           continuous integral using cubature points is
%           integral_lowL^upL w(x)*f(x) dx=sum_{i=1}^n w_i*f(xi(i))
%           where the equality holds up for all polynomials up to a
%           particular degree. For high-order polynomials and other
%           functions, quadrature integration is just an approximation.
%
%INPUTS: n A positive integer such that 2n-1 is the highest degree to
%          which the quadrature points are accurate. n points are returned
%          by the function.
% algorithm This specifies the type of weighting function and the range of
%          the integral. These points can generally be transformed for
%          integrals over other regions/ with other parameters. Possible
%          values are
%          0 (The default if omitted or an empty matrix is passed) The
%            weighting function is w(x)=1/(sqrt(2*pi))*exp(-x^2/2), the
%            integration interval is (-Inf,Inf). The algorithm of [2] is
%            used with parameters from Table 22.7 in Ch. 22 of [1].
%          1 w(x)=exp(-x^2) on the interval (-Inf,Inf). The algorithm of
%            [2] is used with parameters from Table 22.7 in Ch. 22 of [1].
%          2 w(x)=1 on (-1,1). The function GaussLegendrePoints1D is used.
%            Note that the transformation
%            xiNew=xi*(b-a)/2+(b+a)/2;
%            wNew=w*(b-a)/2;
%            can be used to transform the points and weights to the
%            weighting function w(y)=1 on the range (a,b).
%          3 w(x)=(1-x^2)^(c1-1/2) on (-1,1) with c1>-1/2. The algorithm of
%            [2] is used with the values from Table 22.7 of [1]. Formula 10
%            is used for the c1=0 case. This becomes numerically unstable
%            for c1 close to but not equal to zero.
%          4 w(x)=exp(-x) on (0,Inf). The algorithm of [2] is used with the
%            values from Table 22.7 of [1].
%          5 w(x)=x^c1*exp(-x) on (0,Inf). The algorithm of [2] is used
%            with the values from Table 22.7 of [1] for c1>-1.
%          6 w(x)=(1+x)^c1 on (-1,1) for c1>-1. This is a special case of
%            the Jacobi polynomials. The algorithm of [2] is used with
%            values from Table 22.7 of [1].
%          7 w(x)=x^c1 on (0,1), c1>-1. The algorithm of [2] is used. The
%            three-term recursion was derived by explicitly evaluating
%            the integrals with the desired orthogonality constraints
%            until a pattern could be identified for an arbitrary order.
%          8 w(x)=|x|^c1 on (-1,1), c1>=0 and c1 is an integer. The
%            algorithm of [2] is used. The three-term recursion was derived
%            by explicitly evaluating the integrals with the desired
%            orthogonality constraints until a pattern could be identified
%            for an arbitrary order. Separate patterns for even and odd
%            integers were found.
%          9 w(x)=|x|^c1*exp(-x^2) on (-Inf,Inf), c1>=0 and c1 is an
%            integer. The algorithm of [2] is used. The three-term
%            recursion was derived by explicitly evaluating the integrals
%            with the desired orthogonality constraints until a pattern
%            could be identified for an arbitrary order. Separate patterns
%            for even and odd integers were found.
%         10 w(x)=1/sqrt(1-x^2) on (-1,1). The formula in 25.4.38 in [1] is
%            used. Note that the transformation
%            xiNew=(b+a)/2+xi*(b-a)/2; can be used to change the
%            quadrature points to the weighting function
%            w(y)=1/sqrt((y-a)*(b-y)) on (a,b).
%         11 w(x)=sqrt(1-x^2) on (-1,1). The formula in 25.4.40 in [1] is
%            used. Note that the transformation xiNew=(b+a)/2+xi*(b-a)/2
%            can be used to change the quadrature points to the weighting
%            function w(y)=sqrt((y-a)*(b-y)) on (a,b).
%         12 w(x)=sqrt(x/(1-x)) on (0,1). The formula in 25.4.42 in [1] is
%            used. Note that the transformation xiNew=a+(b-a)*xi can be
%            used to change the quadrature points to the weighting
%            function w(x)=sqrt((x-a)/(b-x)) on (a,b).
%      c1 This is a parameter that is only used with certain algorithms
%         described above.
%
%OUTPUTS: xi A 1 X n vector containing the quadrature points.
%          w An nX1 vector of the weights associated with the quadrature
%            points. Depending on the algorithm, these may or may not sum
%            to one.
%
%In the algorithms implemented in this function, with the exception of 11
%and 12, this function obtains points by passing parameters for a three-
%point recursion to the orthoPolyZerosFromRecur function.
%
%Note that the function linCubPoints2MultiDim can be given a handle to this
%function to produce multi-dimensional cubature formula.
%
%REFERENCES:
%[1] Abramowitz, M. and Stegun, I. A. (Eds.). "Orthogonal Polynomials."
%    in Ch. 22, 25 in Handbook of Mathematical Functions with Formulas,
%    Graphs, and Mathematical Tables, 9th printing. New York: Dover, 1972.
%[2] G. H. Golub and J. H. Welsh, "Calculation of Gauss quadrature rules,"
%    Mathematics of Computation, vol. 23, pp. 221-230, 1969.
%
%August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(algorithm))
    algorithm=0;
end

switch(algorithm)
    case 0%Weighting function w(x)=1/(sqrt(2*pi))*exp(-x^2/2)
          %on (-Inf,Inf)
        mu0=1;
        a=@(i)1;
        b=@(i)0;
        c=@(i)i-1;
        [xi,w]=orthoPolyZerosFromRecur(n,a,b,c,mu0);
    case 1%Weighting function w(x)=exp(-x^2) on (-Inf,Inf)
        mu0=sqrt(pi);
        a=@(i)2;
        b=@(i)0;
        c=@(i)2*(i-1);
        [xi,w]=orthoPolyZerosFromRecur(n,a,b,c,mu0);
    case 2%Weighting function w(x)=1 on (-1,1)
        [xi,w]=GaussLegendrePoints1D(n);
    case 3%Weighting function w(x)=(1-x^2)^(c1-1/2) on (-1,1) with c1>-1/2
          %and c1~=0. Given c1=0, Formula 10 is used.
        if(c1==0)
            [xi,w]=quadraturePoints1D(n,10);
            return;
        end

        mu0=(sqrt(pi)*gamma((1/2)+c1))/gamma(1+c1);
        a=@(i)2*(i+c1-1)/i;
        b=@(i)0;
        c=@(i)(i+2*c1-2)/i;
        [xi,w]=orthoPolyZerosFromRecur(n,a,b,c,mu0);
    case 4%Weighting function w(x)=exp(-x) on (0,Inf)
        mu0=1;
        a=@(i)-1/i;
        b=@(i)(2*(i-1)+1)/i;
        c=@(i)(i-1)/i;
        [xi,w]=orthoPolyZerosFromRecur(n,a,b,c,mu0);
    case 5%Weighting function w(x)=x^c1*exp(-x) on (0,Inf) for c1>-1.
        mu0=gamma(1+c1);
        a=@(i)(-1/i);
        b=@(i)(2*i-1+c1)/i;
        c=@(i)(i-1+c1)/i;
        [xi,w]=orthoPolyZerosFromRecur(n,a,b,c,mu0);
    case 6%Weighting function w(x)=(1+x)^c1 on (-1,1) for c1>-1.
          %This uses a special case of the Jacobi polynomials.
        mu0=2^(c1+1)/(c1+1);
        a=@(i)(2*i+c1-1)*(2*i+c1)*(2*i+c1-2)/(2*i*(i+c1)*(2*i+c1-2));
        b=@(i)(2*i+c1-1)*(-c1^2)/(2*i*(i+c1)*(2*i+c1-2));
        c=@(i)2*(i-1)*(i+c1-1)*(2*i+c1)/(2*i*(i+c1)*(2*i+c1-2));
        [xi,w]=orthoPolyZerosFromRecur(n,a,b,c,mu0);
    case 7%Weighting function w(x)=x^c1 on (0,1) for c1>-1
        mu0=1/(1+c1);
        a=@(k)1;
        b=@(k)(-(c1^2+(2*(k-1)+1)*c1+2*k*(k-1))/((2*(k-1)+c1)*(2*k+c1)));
        c=@(k)((k-1)^2*(c1+k-1)^2/((2*k-3+c1)*(2*k-2+c1)^2*(2*k-1+c1)));
        [xi,w]=orthoPolyZerosFromRecur(n,a,b,c,mu0);
    case 8%Weighting function w(x)=|x|^c1 on (-1,1), c1>=0 and c1 is an
          %integer.
        if(mod(c1,2)==0)%Even integer
            c1=c1/2;
            mu0=2/(1+2*c1);
            a=@(k)1;
            b=@(k)0;
            c=@(k)((k-1+2*c1*mod(k+1,2))^2/((2*k-3+2*c1)*(2*k-1+2*c1)));
        else%Odd integer
            c1=(c1-1)/2;
            mu0=1/(1+c1);
            a=@(k)1;
            b=@(k)0;
            c=@(k)(((k-mod(k,2))/2+mod(k+1,2)*c1)^2/((k-1+c1)*(k+c1)));
        end
        [xi,w]=orthoPolyZerosFromRecur(n,a,b,c,mu0);
    case 9%Weighting function w(x)=|x|^c1*exp(-x^2) on (-Inf,Inf), c1>=0
           %and c1 is an integer.
        if(mod(c1,2)==0)%Even integer
            c1=c1/2;
            mu0=gamma(c1+1/2);
            a=@(k)1;
            b=@(k)0;
            c=@(k)((k-1)/2+c1*mod(k+1,2));
        else%Odd integer
            c1=(c1-1)/2;
            mu0=factorial(c1);
            a=@(k)1;
            b=@(k)0;
            c=@(k)((k-mod(k,2))/2+c1*mod(k+1,2));
        end
        [xi,w]=orthoPolyZerosFromRecur(n,a,b,c,mu0);
    case 10%Weighting function w(x)=1/sqrt(1-x^2) on (-1,1) from 25.4.38 in
          %[1].
        xi=cos((2*(1:n)-1)*pi/(2*n));
        w=repmat(pi/n,[n,1]);
    case 11%Weighting function w(x)=sqrt(1-x^2) on (-1,1) from 25.4.40 in
          %[1].
        xi=cos((1:n)*pi/(n+1));
        w=(pi/(n+1))*sin((1:n)'*pi/(n+1)).^2;
    case 12%Weighting function w(x)=sqrt(x/(1-x)) on (0,1) from 25.4.42 in
          %[1].
        xi=cos((pi/2)*((2*(1:n)-1)/(2*n+1))).^2;
        w=(2*pi/(2*n+1))*xi';
    otherwise
        error('Unknown algorithm specified');
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function hermite_poly_phys_values_test ( )

%*****************************************************************************80
%
%% HERMITE_POLY_PHYS_VALUES_TEST demonstrates the use of HERMITE_POLY_PHYS_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_POLY_PHYS_VALUES_TEST:\n' );
  fprintf ( 1, '  HERMITE_POLY_PHYS_VALUES stores values of\n' );
  fprintf ( 1, '  the physicist''s Hermite polynomials.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N      X            H(N,X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, x, fx ] = hermite_poly_phys_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %4d  %12f  %24.16e\n', n, x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Averaging Quaternions
% Since quaternions are not regular vectors, but rather representations
% of orientation, an average quaternion cannot just be obtained by taking
% a weighted mean. This function implements the work done by paper by
% F. Landis Merkley to calculate the average quaternion. The algorithm
% explained by F. Landis Markley at:
% http://www.acsu.buffalo.edu/~johnc/ave_quat07.pdf
% For this particular implementation, I would also like to reference Mandar
% Harshe:
% http://www-sop.inria.fr/members/Mandar.Harshe/knee-joint/html/index.html
%
% This algorithm is compared by rotqrmean from VoiceBox and found to
% produce quite similar results, yet it is more elegant, much simpler to
% implement and follow. (Though, there might be difference in signs)
%
% Usage :
% Q is an Mx4 matrix, where each row stores a quaternion to be averaged.
% In return, the function outputs Qavg, which is a single quaternion
% corresponding to the average.
%
% Tolga Birdal

function [Qwavg]=weighted_avg_quat(Q, weights)

% Form the symmetric accumulator matrix
A=zeros(4,4);
M=size(Q,1);

for i=1:M
    q=Q(i,:)';
    A=q*q'*weights(i)+A; % rank 1 update
end

% scale
A=A/sum(weights);

% Get the eigenvector corresponding to largest eigen value
[Qwavg, Eval]=eigs(A,1);

end

repo: ekfmonoslam


text:
function ShuZiZiMu = DuQuSZZM(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,...
                              s11,s12,s13,s14,s15,s16,s17,s18,s19,s20,...
                              s21,s22,s23,s24,s25,s26,s27,s28,s29,s30,...
                              s31,s32,s33,s34)
%0-9
%10-14 ABCDE 15-19FGHJK 20-24 LMNPQ 25-29 RSTUV 30-33 WXYZ
ShuZiZiMu(:,:,1)   =imresize(im2bw(s1, graythresh(rgb2gray(s1))), [110 55],'bilinear');
ShuZiZiMu(:,:,2)   =imresize(im2bw(s2, graythresh(rgb2gray(s2))), [110 55],'bilinear');
ShuZiZiMu(:,:,3)   =imresize(im2bw(s3, graythresh(rgb2gray(s3))), [110 55],'bilinear');
ShuZiZiMu(:,:,4)   =imresize(im2bw(s4, graythresh(rgb2gray(s4))), [110 55],'bilinear');
ShuZiZiMu(:,:,5)   =imresize(im2bw(s5, graythresh(rgb2gray(s5))), [110 55],'bilinear');
ShuZiZiMu(:,:,6)   =imresize(im2bw(s6, graythresh(rgb2gray(s6))), [110 55],'bilinear');
ShuZiZiMu(:,:,7)   =imresize(im2bw(s7, graythresh(rgb2gray(s7))), [110 55],'bilinear');
ShuZiZiMu(:,:,8)   =imresize(im2bw(s8, graythresh(rgb2gray(s8))), [110 55],'bilinear');
ShuZiZiMu(:,:,9)   =imresize(im2bw(s9, graythresh(rgb2gray(s9))), [110 55],'bilinear');
ShuZiZiMu(:,:,10)  =imresize(im2bw(s10,graythresh(rgb2gray(s10))),[110 55],'bilinear');
ShuZiZiMu(:,:,11)  =imresize(im2bw(s11,graythresh(rgb2gray(s11))),[110 55],'bilinear');
ShuZiZiMu(:,:,12)  =imresize(im2bw(s12,graythresh(rgb2gray(s12))),[110 55],'bilinear');
ShuZiZiMu(:,:,13)  =imresize(im2bw(s13,graythresh(rgb2gray(s13))),[110 55],'bilinear');
ShuZiZiMu(:,:,14)  =imresize(im2bw(s14,graythresh(rgb2gray(s14))),[110 55],'bilinear');
ShuZiZiMu(:,:,15)  =imresize(im2bw(s15,graythresh(rgb2gray(s15))),[110 55],'bilinear');
ShuZiZiMu(:,:,16)  =imresize(im2bw(s16,graythresh(rgb2gray(s16))),[110 55],'bilinear');
ShuZiZiMu(:,:,17)  =imresize(im2bw(s17,graythresh(rgb2gray(s17))),[110 55],'bilinear');
ShuZiZiMu(:,:,18)  =imresize(im2bw(s18,graythresh(rgb2gray(s18))),[110 55],'bilinear');
ShuZiZiMu(:,:,19)  =imresize(im2bw(s19,graythresh(rgb2gray(s19))),[110 55],'bilinear');
ShuZiZiMu(:,:,20)  =imresize(im2bw(s20,graythresh(rgb2gray(s20))),[110 55],'bilinear');
ShuZiZiMu(:,:,21)  =imresize(im2bw(s21,graythresh(rgb2gray(s21))),[110 55],'bilinear');
ShuZiZiMu(:,:,22)  =imresize(im2bw(s22,graythresh(rgb2gray(s22))),[110 55],'bilinear');
ShuZiZiMu(:,:,23)  =imresize(im2bw(s23,graythresh(rgb2gray(s23))),[110 55],'bilinear');
ShuZiZiMu(:,:,24)  =imresize(im2bw(s24,graythresh(rgb2gray(s24))),[110 55],'bilinear');
ShuZiZiMu(:,:,25)  =imresize(im2bw(s25,graythresh(rgb2gray(s25))),[110 55],'bilinear');
ShuZiZiMu(:,:,26)  =imresize(im2bw(s26,graythresh(rgb2gray(s26))),[110 55],'bilinear');
ShuZiZiMu(:,:,27)  =imresize(im2bw(s27,graythresh(rgb2gray(s27))),[110 55],'bilinear');
ShuZiZiMu(:,:,28)  =imresize(im2bw(s28,graythresh(rgb2gray(s28))),[110 55],'bilinear');
ShuZiZiMu(:,:,29)  =imresize(im2bw(s29,graythresh(rgb2gray(s29))),[110 55],'bilinear');
ShuZiZiMu(:,:,30)  =imresize(im2bw(s30,graythresh(rgb2gray(s30))),[110 55],'bilinear');
ShuZiZiMu(:,:,31)  =imresize(im2bw(s31,graythresh(rgb2gray(s31))),[110 55],'bilinear');
ShuZiZiMu(:,:,32)  =imresize(im2bw(s32,graythresh(rgb2gray(s32))),[110 55],'bilinear');
ShuZiZiMu(:,:,33)  =imresize(im2bw(s33,graythresh(rgb2gray(s33))),[110 55],'bilinear');
ShuZiZiMu(:,:,34)  =imresize(im2bw(s34,graythresh(rgb2gray(s34))),[110 55],'bilinear');

author: hangxyz

repo: License-Plate-Recognition-by-MATLAB


text:
function [M, H, noiseFractions] = hyperNacp(M, h, w)
% HYPERMNF Performs the noise adjusted principal component transform (NACP)
%  hyperMnf performs the noise adjust principal component transform on the
% data and uses spatial (row) offsets of the data to estimate the
% covariance matrix of the data.
%
% Usage
%   M = hyperNacp(M, h, w)
% Inputs
%   M - 2D matrix (p x N)
%   h - height of image in pixels
%   w - width of image in pixels
% Outputs
%   M - 2D transformed data
%   H - 2D transformation matrix
%   noiseFractions - Estimates of the noise fraction for each band
%
% References
%   C-I Change and Q Du, "Interference and Noise-Adjusted Principal
% Components Analysis," IEEE TGRS, Vol 36, No 5, September 1999.

[p, N] = size(M);

% Remove mean from data
u = mean(M.').';
for k=1:N
    M(:,k) = M(:,k) - u;
end

% Compute to rotation of the signal+noise
sigmaZ = hyperCov(M);
M = hyperConvert3d(M, h, w, p);

% Estimate the covariance of the noise.
dX = zeros(h-1, w, p);
for i=1:(h-1)
    dX(i, :, :) = M(i, :, :) - M(i+1, :, :);
end
dX = hyperConvert2d(dX);

% Compute the covariance of the noise signal estimate.
sigmaN = hyperCov(dX);

% Orthonormalize the noise subspace.
[U,deltaN,E] = svd(sigmaN);
F = E*inv(sqrt(deltaN));  % Rotation components of noise orthonormalized
% F now whitens the noise.

% Rotates the signal+noise cov so that the noise is whitened (all noise
% powers are equal)
sigmaAdj = F'*sigmaZ*F;

[U,gammaAdj,G] = svd(sigmaAdj);
H = G*F;

% Compute noise fractions
noiseFractions = diag(gammaAdj);

% Perform transform
M = H*hyperConvert2d(M);

author: isaacgerg

repo: matlabHyperspectralToolbox


text:
clear all; close all;
I=imread('onion.png');
J=rgb2gray(I);
K=fft2(J);
L=fftshift(K);
M=ifft2(K);
figure;
subplot(121);
imshow(uint8(abs(L)/198));
subplot(122);
imshow(uint8(M));

author: HuangCongQing


text:
function C = times(A,B)
%TIMES Element-wise multiplication for ktensor.
%
%   TIMES(A,B) denotes element-by-element multiplication.
%
%   C = TIMES(A,B) is called for the syntax 'A .* B' when A or B is a
%   tensor.
%
%   See also KTENSOR, SPTENSOR/TIMES.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if ~isequal(size(A),size(B))
    error('Must be two tensors of the same size');
end

switch class(B)
    case {'sptensor','tensor'}
        % Call back to sptensor version.
        C = times(B,A);
        return;
    otherwise
        error('Invalid second argument for ktensor/times');
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm