text:
function taser_about( )
%TASER_ABOUT Brings up the "About" window in the TASER menu
%
% Author: Tim Cox, NRL, and Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

fig_hand=figure('Name', 'Taser: About', 'NumberTitle', 'off', ...
    'MenuBar', 'none', 'Toolbar', 'none');

% Draw title
text_ax = axes('Parent',fig_hand,'Position',[0.05 .7 0.5 .3]);
axis(text_ax, 'off');
title_text = text(.5, 0.6, 'TASER', 'FontSize', 20, 'FontWeight', 'bold', 'FontAngle', 'italic', 'HorizontalAlignment', 'center', 'Parent', text_ax);
if isdeployed()&&exist('taser_compile_date','file')
    set(title_text,'Position',[.5 .7]); % Make room for version text
    text(.5, 0.4, ['Version ' taser_compile_date], 'FontSize', 12, 'HorizontalAlignment', 'center', 'Parent', text_ax)
end
text(.5, 0.2, 'Naval Research Laboratory', 'FontSize', 10, 'HorizontalAlignment', 'center', 'Parent', text_ax)
text(.5, 0.05, 'National Geospatial-Intelligence Agency', 'FontSize', 10, 'HorizontalAlignment', 'center', 'Parent', text_ax)

% Draw scroll bar with text info
license_ui = uicontrol(fig_hand, 'Style', 'edit',...
    'Units', 'normalized', ...
    'Position', [.05 .05 .55 .6], ...
    'min',0,'max',5,... % No effect, except max must be greater than min
    'HorizontalAlign','left',...
    'BackgroundColor',[1 1 1],...
    'enable','inactive');
scrollbar_string = '';
credits_path = fullfile(fileparts(mfilename('fullpath')),'taser_credits.txt');
if exist(credits_path,'file')
    fid = fopen(credits_path);
    scrollbar_string = fread(fid,Inf,'*char').';
    fclose(fid);
end
% Assumes this file is in Tools/MatlabWorkbench/about path of MATLAB SAR Toolbox
license_path = fullfile([fileparts(mfilename('fullpath')) filesep '..' filesep '..' filesep '..'],'license.txt');
if exist(license_path,'file')
    fid = fopen(license_path);
    fgetl(fid);
    license_strings = fread(fid,Inf,'*char').';
    fclose(fid);
    % Concatenate so that uicontrol will wrap properly.
    license_strings = regexprep(license_strings, '\r\n(?!\r\n)(?<!\r\n\r\n)', ''); % Single occurence only
    scrollbar_string = [scrollbar_string sprintf('\r\n\r\nLicense:\r\n\r\n') license_strings];
end
set(license_ui,'String',scrollbar_string);

% Draw seals
nga_ax = axes('Parent', fig_hand,'Position', [.6 .55 .4 .4]);
nrl_ax = axes('Parent', fig_hand,'Position', [.6 .05 .4 .4]);
image(imread(fullfile(fileparts(mfilename('fullpath')),'NGA_Seal.png'),...
    'BackgroundColor',get(fig_hand,'Color')),'Parent',nga_ax);
image(imread(fullfile(fileparts(mfilename('fullpath')),'NRL_Seal.png'),...
    'BackgroundColor',get(fig_hand,'Color')),'Parent',nrl_ax);
set([nga_ax nrl_ax], 'DataAspectRatio', [1 1 1]);
axis([nga_ax nrl_ax], 'off');

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
close all
clear

%% Setup and Parameters
x0 = [-5; -5; 0; 0];
time_total = 20.0;
dt = 0.2;
P = 100*eye(4);
Q = 10*eye(4);
R = eye(2);
N = 8;
xmin = [-5; -5; -5; -5];
xmax = [5; 5; 5; 5];
umin = [-1; -1];
umax = [1; 1];

%% Discrete-time double integrator 2D
system.dt = dt;
system.A = [1 0 dt 0;
    0 1 0 dt;
    0 0 1 0;
    0 0 0 1];
system.B = [0.5*dt^2 0;
    0 0.5*dt^2;
    dt 0;
    0 dt];
system.xl = xmin;
system.xu = xmax;
system.ul = umin;
system.uu = umax;

%% MPC-CBF parameters
params.Q = Q;
params.R = R;
params.P = P;
params.N = N;
params.gamma = 0.5;

%% Obstacle
obs.pos = [-2; -2.25];
obs.r = 1.5;

%% Simulate MPC-CBF
gamma_list = linspace(0.1, 0.5, 5);
controller_mpc_cbf_list = {};
for i = 1:size(gamma_list, 2)
    new_params = params;
    new_params.N = 5;
    new_params.gamma = gamma_list(i);
    controller_mpc_cbf = MPCCBF(x0, system, new_params);
    controller_mpc_cbf.obs = obs;
    controller_mpc_cbf.sim(time_total);
    controller_mpc_cbf_list{i} = controller_mpc_cbf;
end

%% Computational time benchmark
for i = 1:size(gamma_list, 2)
    controller_mpc_cbf = controller_mpc_cbf_list{i};
    fprintf('Computational time for MPC-CBF3: mean %.3f, std %.3f, min %.3f, max %.3f, input cost %.3f, min dist %f\n',...
    [mean(controller_mpc_cbf.solvertime),...
    std(controller_mpc_cbf.solvertime),...
    min(controller_mpc_cbf.solvertime),...
    max(controller_mpc_cbf.solvertime),...
    controller_mpc_cbf.u_cost,...
    min(controller_mpc_cbf.distlog)]);
end

author: HybridRobotics

repo: NMPC-DCLF-DCBF


text:
function value = ddot ( n, dx, incx, dy, incy )

%*****************************************************************************80
%
%% DDOT forms the dot product of two vectors.
%
%  Discussion:
%
%    This routine uses unrolled loops for increments equal to one.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 June 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vectors.
%
%    Input, real DX(*), the first vector.
%
%    Input, integer INCX, the increment between successive entries in DX.
%
%    Input, real DY(*), the second vector.
%
%    Input, integer INCY, the increment between successive entries in DY.
%
%    Output, real VALUE, the sum of the product of the
%    corresponding entries of DX and DY.
%
  ddot = 0.0;
  value = 0.0;

  if ( n <= 0 )
    return
  end
%
%  Code for unequal increments or equal increments
%  not equal to 1.
%
  if ( incx ~= 1 | incy ~= 1 )

    if ( 0 <= incx )
      ix = 1;
    else
      ix = ( - n + 1 ) * incx + 1;
    end

    if ( 0 <= incy )
      iy = 1;
    else
      iy = ( - n + 1 ) * incy + 1;
    end

    for i = 1 : n
      value = value + dx(ix) * dy(iy);
      ix = ix + incx;
      iy = iy + incy;
    end
%
%  Code for both increments equal to 1.
%
  else

    m = mod ( n, 5 );

    for i = 1 : m
      value = value + dx(i) * dy(i);
    end

    for i = m+1 : 5 : n

      value = value + dx(i  ) * dy(i  ) ...
                    + dx(i+1) * dy(i+1) ...
                    + dx(i+2) * dy(i+2) ...
                    + dx(i+3) * dy(i+3) ...
                    + dx(i+4) * dy(i+4);
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ n_data, n, k, x, b ] = bernstein_poly_01_values ( n_data )

%*****************************************************************************80
%
%% BERNSTEIN_POLY_01_VALUES returns some values of the Bernstein polynomials.
%
%  Discussion:
%
%    The Bernstein polynomials are assumed to be based on [0,1].
%
%    The formula for the Bernstein polynomials is
%
%      B(N,I)(X) = [N!/(I!*(N-I)!)] * (1-X)**(N-I) * X**I
%
%    In Mathematica, the function can be evaluated by:
%
%      Binomial[n,i] * (1-x)^(n-i) * x^i
%
%  First values:
%
%    B(0,0)(X) = 1
%
%    B(1,0)(X) =      1-X
%    B(1,1)(X) =                X
%
%    B(2,0)(X) =     (1-X)**2
%    B(2,1)(X) = 2 * (1-X)    * X
%    B(2,2)(X) =                X**2
%
%    B(3,0)(X) =     (1-X)**3
%    B(3,1)(X) = 3 * (1-X)**2 * X
%    B(3,2)(X) = 3 * (1-X)    * X**2
%    B(3,3)(X) =                X**3
%
%    B(4,0)(X) =     (1-X)**4
%    B(4,1)(X) = 4 * (1-X)**3 * X
%    B(4,2)(X) = 6 * (1-X)**2 * X**2
%    B(4,3)(X) = 4 * (1-X)    * X**3
%    B(4,4)(X) =                X**4
%
%  Special values:
%
%    B(N,I)(X) has a unique maximum value at X = I/N.
%
%    B(N,I)(X) has an I-fold zero at 0 and and N-I fold zero at 1.
%
%    B(N,I)(1/2) = C(N,K) / 2**N
%
%    For a fixed X and N, the polynomials add up to 1:
%
%      Sum ( 0 <= I <= N ) B(N,I)(X) = 1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, the degree of the polynomial.
%
%    Output, integer K, the index of the polynomial.
%
%    Output, real X, the argument of the polynomial.
%
%    Output, real B, the value of the polynomial B(N,K)(X).
%
  n_max = 15;

  b_vec = [ ...
     0.1000000000000000E+01, ...
     0.7500000000000000E+00, ...
     0.2500000000000000E+00, ...
     0.5625000000000000E+00, ...
     0.3750000000000000E+00, ...
     0.6250000000000000E-01, ...
     0.4218750000000000E+00, ...
     0.4218750000000000E+00, ...
     0.1406250000000000E+00, ...
     0.1562500000000000E-01, ...
     0.3164062500000000E+00, ...
     0.4218750000000000E+00, ...
     0.2109375000000000E+00, ...
     0.4687500000000000E-01, ...
     0.3906250000000000E-02 ];

  k_vec = [ ...
    0, ...
    0, 1, ...
    0, 1, 2, ...
    0, 1, 2, 3, ...
    0, 1, 2, 3, 4 ];

  n_vec = [ ...
    0, ...
    1, 1, ...
    2, 2, 2, ...
    3, 3, 3, 3, ...
    4, 4, 4, 4, 4 ];

  x_vec = [ ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    k = 0;
    x = 0.0;
    b = 0.0;
  else
    n = n_vec(n_data);
    k = k_vec(n_data);
    x = x_vec(n_data);
    b = b_vec(n_data);
  end

  return
end


text:
function writeMesh_off(fileName, vertices, faces)
%WRITEMESH_OFF Writes a mesh into a text file in OFF format.
%
%   writeMesh_off(FNAME, V, F)
%
%   Example
%   writeMesh_off
%
%   See also
%      meshes3d, readMesh_off, writeMesh_ply

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2018-04-26,    using Matlab 9.4.0.813654 (R2018a)
% Copyright 2018 INRA - Cepia Software Platform.

if ~ischar(fileName)
    error('First argument must contain the name of the file');
end

% optionnaly parses data
if isstruct(vertices)
    faces = vertices.faces;
    vertices = vertices.vertices;
end

% open file for writing text
f = fopen(fileName, 'wt');
if (f == -1)
	error('Couldn''t open the file %s', fileName);
end

% write the header line
fprintf(f, 'OFF\n');

% write number of vertices and of faces
nVertices = size(vertices, 1);
nFaces = size(faces, 1);
if iscell(faces)
    nFaces = length(faces);
end
fprintf(f, '%d %d\n', nVertices, nFaces);

% Write vertex info
format = '%g %g %g\n';
for iv = 1:nVertices
    fprintf(f, format, vertices(iv, :));
end

% Write face info
if isnumeric(faces)
    % simply write face vertex indices
    ns = size(faces, 2);
    format = ['%d' repmat(' %d', 1, ns) '\n'];
    for iFace = 1:nFaces
        fprintf(f, format, ns, faces(iFace, :)-1);
    end
else
    % if faces are stored in a cell array, the number of vertices in each
    % face may be different, and we need to process each face individually
    for iFace = 1:nFaces
        ns = length(faces{iFace});
        format = ['%d' repmat(' %d', 1, ns) '\n'];
        fprintf(f, format, ns, faces{iFace}-1);
    end
end

% close the file
fclose(f);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% The COBRAToolbox: testFindExtremePathway.m
%
% Purpose:
%     - testFindExtremePathway tests the functionality of findExtremePathway.
%
% Authors:
%     - Sylvain Arreckx March 2017
%     - Laurent Heirendt March 2017
%
% Test problem from
%     "System Modeling in Cellular Biology: From Concepts to Nuts and Bolts",
%     section "Stoichiometric and Constraint Based Modeling", MIT Press

% save the current path
currentDir = pwd;

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

model.S = [1, 0, 0, 0,-1,-1,-1, 0, 0, 0
           0, 1, 0, 0, 1, 0, 0,-1,-1, 0
           0, 0, 0, 0, 0, 1, 0, 1, 0,-1
           0, 0, 0, 0, 0, 0, 1, 0, 0,-1
           0, 0, 0,-1, 0, 0, 0, 0, 0, 1
           0, 0,-1, 0, 0, 0, 0, 0, 1, 1];

% define list of solver packages
solverPkgs = {'gurobi', 'glpk'};

obj = [0.5985; 0.4709; 0.6959; 0.6999; 0.6385; 0.0336; 0.0688; 0.3196; 0.5309; 0.6544; 0.4076; 0.8200];

refV = [0.2500; 0; 0.2500; 0; 0.2500; 0; 0; 0; 0.2500; 0];
tol = 1e-14;

for k = 1:length(solverPkgs)

    % print success message
    fprintf('   Testing findExtremePathway using %s ... ', solverPkgs{k});

    % change the solver
    solverOK = changeCobraSolver(solverPkgs{k}, 'LP', 0);

    if solverOK == 1
%
%         % testing findExtremePathway with different arguments
%         if isfield(model, 'revRxns')
%             model = rmfield(model, 'revRxns');
%         end
        %testing missing lb/ub
        try
            v = findExtremePathway(model, obj);
        catch ME
            assert(length(ME.message) > 0)
        end

        model.ub = [-1,  1, -1, -1, -1, -1, -1,  1, -1, -1];
        model.lb = [ 0, -1,  0,  0,  0,  0,  0, -1,  0,  0];
        v = findExtremePathway(model);
        assert(all(abs(model.S * v) < tol))

        v = findExtremePathway(model, obj);
        assert(all(refV == v))
        assert(all(model.S * v == zeros(size(model.S, 1), 1)))

        model2.S = [-1,  0,  1;
                     1, -1,  0;
                     0,  1, -1];

        %solverOK = changeCobraSolver('gurobi', 'LP', 0);

        [x, output] = findExtremePool(model2);

        assert(all(abs(model2.S * x) < tol))

        % delete generated files
        delete('*.ine');
        delete('*.ext');
    end

    % print success message
    fprintf('Done\n');
end

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
%% efw
% Below is a demonstration of the features of the |efw| function

%% Syntax
% |efw(hf);|

%% Description
% The |efw| function, the export figure widget, adds a push button to a
% figure toolbar to link with the |export_fig| function. Press the button
% to start exporting a figure. Users can specify file names, formats,
% resolution and also additional export_fig options.
%%
% The Export Figure Widget requires the external function |export_fig|
% created by Oliver Woodford and Yair Altman. It can be obtained from the
% Mathworks Central File Exchange:
%%
% <http://www.mathworks.com/matlabcentral/fileexchange/23629-export-fig>

%% Examples

%% Using the preloaded Export Figure Widget (efw) in |cFigure|
% A |cFigure| window has the Export Figure Widget loaded by default. To use
% |efw| click on the icon in the toolbar.
%
% <<efw_icon_menubar.png>>

%%
% Pressing the |efw| button will open a basic |inputdlg| allowing users to
% specify all the usual export_fig options. Hints are given in brackets
% behind the input labels. The default entries are altered according to the
% previous usage within the current figure.

%%
% <<efw_inputdlg.png>>

%% Use with |figure|
% The default MATLAB |figure| does not contain the Export Figure Widget. To
% load it here the user must enter |efw;| either before or after axes are
% created. e.g.:
%%
% |figure; surf(peaks(25)); axis equal; axis tight; efw;|
%%
% which is equivalent to:
%%
% |figure; efw; surf(peaks(25)); axis equal; axis tight;|

%% Tip
% You can create your own figure function that simply contains |figure;
% efw;| to create a standard MATLAB figure containing the Export Figure
% Widget.

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function z = h2z(h);

% Z = h2z(H)
%
% Hybrid to Admittance transformation
% only for 2-by-2 matrices
%
% martie 27

if h(2,2) == 0
  disp('correspondent admittance matrix non-existent');
else
z(1,1) = h(1,1) - h(1,2)*h(2,1)/h(2,2);
z(1,2) = h(1,2)/h(2,2);
z(2,1) = -h(2,1)/h(2,2);
z(2,2) = 1/h(2,2);
end;

author: Sable

repo: mcbench-benchmarks


text:
function planner(rows,cols,obs,start,goal)
%
% path planner simulator
%
density = size(obs,2)/(rows*cols);% find obstacle density
ndensity = 9;
for d = 9:10:89
    if d <= density
        ndensity = d;
    else
        break;
    end
end
csetname = sprintf('curveset_%d_%d_%d',rows,cols,ndensity);
if(evalin('base',sprintf('exist(''%s'',''var'')',csetname))==0)% try to load curveset from matlab runtime memory
    if(evalin('base',sprintf('exist(''%s.mat'',''file'')',csetname))==0)% try to load curveset from hard drive
        disp('Generating curve set');
        evalin('base',sprintf('gen_curves(%d,%d,%d)',rows,cols,ndensity));% generate curveset
    else
        disp('Loading curve set from file');
        evalin('base',sprintf('load %s.mat',csetname));
    end
    curveset = evalin('base', csetname);
else
    disp('Loading curve set from memory');
    curveset = evalin('base', csetname);
end
clear csetname;
sol1 = forward_chain(rows,cols,obs,goal,start,100);% find time to reach shallowest solution
if (sol1 == -1)
    disp('Run program after reloading world');% Dynamics changed
    return;
end
la = find_good_la(rows,cols,obs,goal,start,sol1,curveset);% estimate lookahead to use
world = extend_world(rows,cols,obs,goal,la);% extend world in time
pot_vals = calc_pot_values(world);% Calculate potential values
[res path] = traverse_world(start,pot_vals);% Follow negative potential gradient
if (res == 0)
    disp('Run program after reloading world');% Dynamics changed or Goal became unreachable
    return;
end
show_path(world,path)% display path traversed
end


text:
%SerialLink.plot3d Graphical display and animation of solid model robot
%
% R.plot3d(Q, options) displays and animates a solid model of the robot.
% The robot is displayed at the joint angle Q (1xN), or
% if a matrix (MxN) it is animated as the robot moves along the M-point trajectory.
%
% Options::
%
% 'color',C         A cell array of color names, one per link.  These are
%                   mapped to RGB using colorname().  If not given, colors
%                   come from the axis ColorOrder property.
% 'alpha',A         Set alpha for all links, 0 is transparant, 1 is opaque
%                   (default 1)
% 'path',P          Overide path to folder containing STL model files
% 'workspace', W    Size of robot 3D workspace, W = [xmn, xmx ymn ymx zmn zmx]
% 'floorlevel',L    Z-coordinate of floor (default -1)
%-
% 'delay',D         Delay betwen frames for animation (s)
% 'fps',fps         Number of frames per second for display, inverse of 'delay' option
% '[no]loop'        Loop over the trajectory forever
% '[no]raise'       Autoraise the figure
% 'movie',M         Save frames as files in the folder M
%-
% 'scale',S         Annotation scale factor
% 'ortho'           Orthographic view (default)
% 'perspective'     Perspective view
% 'view',V          Specify view V='x', 'y', 'top' or [az el] for side elevations,
%                   plan view, or general view by azimuth and elevation
%                   angle.
%-
% '[no]wrist'       Enable display of wrist coordinate frame
% 'xyz'             Wrist axis label is XYZ
% 'noa'             Wrist axis label is NOA
% '[no]arrow'       Display wrist frame with 3D arrows
%-
% '[no]tiles'       Enable tiled floor (default true)
% 'tilesize',S      Side length of square tiles on the floor (default 0.2)
% 'tile1color',C   Color of even tiles [r g b] (default [0.5 1 0.5]  light green)
% 'tile2color',C   Color of odd tiles [r g b] (default [1 1 1] white)
%-
% '[no]jaxes'       Enable display of joint axes (default true)
% '[no]joints'      Enable display of joints
%-
% '[no]base'        Enable display of base shape
%
% Notes::
% - Solid models of the robot links are required as STL ascii format files,
%   with extensions .stl
% - Suitable STL files can be found in the package ARTE: A ROBOTICS TOOLBOX
%   FOR EDUCATION by Arturo Gil, https://arvc.umh.es/arte
% - The root of the solid models is an installation of ARTE with an empty
%   file called arte.m at the top level
% - Each STL model is called 'linkN'.stl where N is the link number 0 to N
% - The specific folder to use comes from the SerialLink.model3d property
% - The path of the folder containing the STL files can be specified using
%   the 'path' option
% - The height of the floor is set in decreasing priority order by:
%   - 'workspace' option, the fifth element of the passed vector
%   - 'floorlevel' option
%   - the lowest z-coordinate in the link1.stl object
%
% Authors::
% - Peter Corke, based on existing code for plot()
% - Bryan Moutrie, demo code on the Google Group for connecting ARTE and RTB
% - Don Riley, function rndread() extracted from cad2matdemo (MATLAB
%   File Exchange)
%
% See also SerialLink.plot, plotbotopt3d, SerialLink.animate, SerialLink.teach, SerialLink.fkine.

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

function plot3d(robot, q, varargin)

    if robot.mdh
        error('RTB:plot3d:badmodel', '3D models are defined for standard, not modified, DH parameters');
    end

    clf
    opt = plot_options(robot, varargin);

    %-- load the shape if need be

    nshapes = robot.n+1;

    if isempty(robot.faces)
        % no 3d model defined, let's try to load one

        if isempty(opt.path)
            % first find the path to the models
            pth = which('arte.m');
            if ~pth
                error('RTB:plot3d:nomodel', 'no 3D model found, install the RTB contrib zip file');
            end

            % find the path to this specific model
            pth = fullfile(fileparts(pth), 'robots', robot.model3d);
        else
            pth = opt.path;
        end

        % now load the STL files
        robot.points = cell(1, robot.n+1);
        robot.faces = cell(1, robot.n+1);
        for i=1:nshapes
            [F,P] = rndread( fullfile(pth, sprintf('link%d.stl', i-1)) );
            robot.points{i} = P;
            robot.faces{i} = F;
        end
    end

    % if a base is specified set the floor height to this
    if isempty(opt.workspace)
        % workspace not provided, fall through the options for setting floor level
        if ~isempty(opt.floorlevel)
            opt.ws(5) = opt.floorlevel;
        elseif opt.base
            mn = min( robot.points{1} );
            opt.ws(5) = mn(3);
        end
    end
    opt.floorlevel = opt.ws(5);

% TODO
% should test if the plot exists, pinch the logic from plot()

    %-- set up to plot
    % create an axis
    ish = ishold();
    if ~ishold
        % if hold is off, set the axis dimensions
        axis(opt.ws);
        set(gca, 'ZLimMode', 'manual');
        axis(opt.ws);
        hold on
    end

    if opt.raise
        % note this is a very time consuming operation
        figure(gcf);
    end

    if strcmp(opt.projection, 'perspective')
        set(gca, 'Projection', 'perspective');
    end

    grid on
    xlabel('X'); ylabel('Y'); zlabel('Z');

    %--- create floor if required
    if opt.tiles
        create_tiled_floor(opt);
    end

    %--- configure view and lighting
    if isstr(opt.view)
        switch opt.view
            case 'top'
                view(0, 90);
            case 'x'
                view(0, 0);
            case 'y'
                view(90, 0)
            otherwise
                error('rtb:plot3d:badarg', 'view must be: x, y, top')
        end
    elseif isnumeric(opt.view) && length(opt.view) == 2
        view(opt.view)
    else
        campos([2 2 1]);
    end

    daspect([1 1 1]);
    light('Position', [0 0 opt.reach*2]);
    light('Position', [1 0.5 1]);

    %-- figure the colors for each shape
    if isempty(opt.color)
        % if not given, use the axis color order
        C = get(gca,'ColorOrder');
    else
        C = [];
        for c=opt.color
            C = [C; colorname(c{1})];
        end
    end

    %--- create the robot
    %  one patch per shape, use hgtransform to animate them later
    group = hggroup('Tag', robot.name);
    ncolors = numrows(C);

    h = [];
    for link=0:robot.n
        if link == 0
            if ~opt.base
                continue;
            end

            patch('Faces', robot.faces{link+1}, 'Vertices', robot.points{link+1}, ...
                'FaceColor', C(mod(link,ncolors)+1,:), 'EdgeAlpha', 0, 'FaceAlpha', opt.alpha);
        else
            h.link(link) = hgtransform('Tag', sprintf('link%d', link), 'Parent', group);

            patch('Faces', robot.faces{link+1}, 'Vertices', robot.points{link+1}, ...
                'FaceColor', C(mod(link,ncolors)+1,:), 'EdgeAlpha', 0, 'FaceAlpha', opt.alpha, ...
            'Parent', h.link(link));
        end
    end

    % enable mouse-based 3D rotation
    rotate3d on

    h.wrist = [];  % HACK, should be trplot
    h.robot = robot;
    h.link = [0 h.link];
    set(group, 'UserData', h);

    robot.animate(q);

    if ~ish
        hold off
    end
end

function opt = plot_options(robot, optin)
    opt.color = [];
    opt.path = [];  % override path
    opt.alpha = 1;

        % timing/looping
    opt.delay = 0.1;
    opt.fps = [];
    opt.loop = false;

    opt.raise = false;

    % general appearance
    opt.scale = 1;

    opt.workspace = [];
    opt.floorlevel = [];

    opt.name = true;
    opt.projection = {'ortho', 'perspective'};
    opt.view = [];

    % tiled floor
    opt.tiles = true;
    opt.tile1color = [0.5 1 0.5];  % light green
    opt.tile2color = [1 1 1];  % white
    opt.tilesize = 0.2;

    % the base or pedestal
    opt.base = true;

    % wrist
    opt.wrist = true;
    opt.wristlabel = {'xyz', 'noa'};
    opt.arrow = true;

    % joint rotation axes
    opt.jaxes = false;

    % misc
    opt.movie = [];

    % build a list of options from all sources
    %   1. the M-file plotbotopt if it exists
    %   2. robot.plotopt
    %   3. command line arguments
    if exist('plotbotopt3d', 'file') == 2
        options = [plotbotopt3d robot.plotopt3d optin];
    else
        options = [robot.plotopt3d optin];
    end

    % parse the options
    [opt,args] = tb_optparse(opt, options);
    if ~isempty(args)
        error(['unknown option: ' args{1}]);
    end

    if ~isempty(opt.fps)
        opt.delay = 1/opt.fps;
    end
    % figure the size of the figure

    if isempty(opt.workspace)
        %
        % simple heuristic to figure the maximum reach of the robot
        %
        L = robot.links;
        if any(L.isprismatic)
            error('Prismatic joint(s) present: requires the ''workspace'' option');
        end
        reach = 0;
        for i=1:robot.n
            reach = reach + abs(L(i).a) + abs(L(i).d);
        end

        % if we have a floor, quantize the reach to a tile size
        if opt.tiles
            reach = opt.tilesize * ceil(reach/opt.tilesize);
        end

        % now create a 3D volume based on this reach
        opt.ws = [-reach reach -reach reach -reach reach];

        % if a floorlevel has been given, ammend the 3D volume
        if ~isempty(opt.floorlevel)
            opt.ws(5) = opt.floorlevel;
        end
    else
        % workspace is provided
        reach = min(abs(opt.workspace));
        if opt.tiles
            % set xy limits to be integer multiple of tilesize
            opt.ws(1:4) = opt.tilesize * round(opt.workspace(1:4)/opt.tilesize);
            opt.ws(5:6) = opt.workspace(5:6);
        end
    end

    opt.reach = reach;

    % update the fundamental scale factor (given by the user as a multiplier) by a length derived from
    % the overall workspace dimension
    %  we need that a lot when creating the robot model
    opt.scale = opt.scale * reach/40;

    % deal with a few options that need to be stashed in the SerialLink object
    % movie mode has not already been flagged
    if opt.movie
        robot.framenum = 0;
        robot.moviepath = opt.movie;
    else
        robot.framenum = [];
    end

    if ~isempty(opt.movie)
        mkdir(opt.movie);
        framenum = 1;
    end

    robot.delay = opt.delay;
    robot.loop = opt.loop;
end

% Extracted from cad2matdemo
% http://www.mathworks.com/matlabcentral/fileexchange/3642-cad2matdemo-m/content/cad2matR2.zip
%
% Copyright (c) 2003, Don Riley
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of the Walla Walla University nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

function [fout, vout, cout] = rndread(filename)

    % Reads CAD STL ASCII files, which most CAD programs can export.
    % Used to create Matlab patches of CAD 3D data.
    % Returns a vertex list and face list, for Matlab patch command.
    %
    % filename = 'hook.stl';  % Example file.
    %
    fid=fopen(filename, 'r'); %Open the file, assumes STL ASCII format.
    if fid == -1
        error('File could not be opened, check name or path.')
    end
    %
    % Render files take the form:
    %
    %solid BLOCK
    %  color 1.000 1.000 1.000
    %  facet
    %      normal 0.000000e+00 0.000000e+00 -1.000000e+00
    %      normal 0.000000e+00 0.000000e+00 -1.000000e+00
    %      normal 0.000000e+00 0.000000e+00 -1.000000e+00
    %    outer loop
    %      vertex 5.000000e-01 -5.000000e-01 -5.000000e-01
    %      vertex -5.000000e-01 -5.000000e-01 -5.000000e-01
    %      vertex -5.000000e-01 5.000000e-01 -5.000000e-01
    %    endloop
    % endfacet
    %
    % The first line is object name, then comes multiple facet and vertex lines.
    % A color specifier is next, followed by those faces of that color, until
    % next color line.
    %
    CAD_object_name = sscanf(fgetl(fid), '%*s %s');  %CAD object name, if needed.
    %                                                %Some STLs have it, some don't.
    vnum=0;       %Vertex number counter.
    report_num=0; %Report the status as we go.
    VColor = 0;
    %
    while feof(fid) == 0                    % test for end of file, if not then do stuff
        tline = fgetl(fid);                 % reads a line of data from file.
        fword = sscanf(tline, '%s ');       % make the line a character string
        % Check for color
        if strncmpi(fword, 'c',1) == 1;    % Checking if a "C"olor line, as "C" is 1st char.
            VColor = sscanf(tline, '%*s %f %f %f'); % & if a C, get the RGB color data of the face.
        end                                % Keep this color, until the next color is used.
        if strncmpi(fword, 'v',1) == 1;    % Checking if a "V"ertex line, as "V" is 1st char.
            vnum = vnum + 1;                % If a V we count the # of V's
            report_num = report_num + 1;    % Report a counter, so long files show status
            if report_num > 249;
                %disp(sprintf('Reading vertix num: %d.',vnum));
                report_num = 0;
            end
            v(:,vnum) = sscanf(tline, '%*s %f %f %f'); % & if a V, get the XYZ data of it.
            c(:,vnum) = VColor;              % A color for each vertex, which will color the faces.
        end                                 % we "*s" skip the name "color" and get the data.
    end
    %   Build face list; The vertices are in order, so just number them.
    %
    fnum = vnum/3;      %Number of faces, vnum is number of vertices.  STL is triangles.
    flist = 1:vnum;     %Face list of vertices, all in order.
    F = reshape(flist, 3,fnum); %Make a "3 by fnum" matrix of face list data.
    %
    %   Return the faces and vertexs.
    %
    fout = F';  %Orients the array for direct use in patch.
    vout = v';  % "
    cout = c';
    %
    fclose(fid);
end

% draw a tiled floor in the current axes
function create_tiled_floor(opt)

    xmin = opt.ws(1);
    xmax = opt.ws(2);
    ymin = opt.ws(3);
    ymax = opt.ws(4);

    % create a colored tiled floor
    xt = xmin:opt.tilesize:xmax;
    yt = ymin:opt.tilesize:ymax;
    Z = opt.floorlevel*ones( numel(yt), numel(xt));
    C = zeros(size(Z));
    [r,c] = ind2sub(size(C), 1:numel(C));
    C = bitand(r+c,1);
    C = reshape(C, size(Z));
    C = cat(3, opt.tile1color(1)*C+opt.tile2color(1)*(1-C), ...
        opt.tile1color(2)*C+opt.tile2color(2)*(1-C), ...
        opt.tile1color(3)*C+opt.tile2color(3)*(1-C));
    [X,Y] = meshgrid(xt, yt);
    surface(X, Y, Z, C, ...
        'FaceColor','texturemap',...
        'EdgeColor','none',...
        'CDataMapping','direct');
end

author: Allopart

repo: rbpf-gmapping


text:
function triangle_points_plot ( file_name, node_xy, node_show, point_num, ...
  point_xy, point_show )

%*****************************************************************************80
%
%% TRIANGLE_POINTS_PLOT plots a triangle and some points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 October 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILE_NAME, the name of the output file.
%
%    Input, real NODE_XY(2,3), the coordinates of the nodes
%    of the triangle.
%
%    Input, integer NODE_SHOW,
%   -1, do not show the triangle, or the nodes.
%    0, show the triangle, do not show the nodes;
%    1, show the triangle and the nodes;
%    2, show the triangle, the nodes and number them.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, real POINT_XY(2,POINT_NUM), the coordinates of the
%    points.
%
%    Input, integer POINT_SHOW,
%    0, do not show the points;
%    1, show the points;
%    2, show the points and number them.
%
  node_num = 3;
  x_ps_max = 576;
  x_ps_max_clip = 594;
  x_ps_min = 36;
  x_ps_min_clip = 18;
  y_ps_max = 666;
  y_ps_max_clip = 684;
  y_ps_min = 126;
  y_ps_min_clip = 108;
%
%  We need to do some figuring here, so that we can determine
%  the range of the data, and hence the height and width
%  of the piece of paper.
%
  x_max = max ( max ( node_xy(1,1:node_num) ), ...
                max ( point_xy(1,1:point_num) ) );
  x_min = min ( min ( node_xy(1,1:node_num) ), ...
                min ( point_xy(1,1:point_num) ) );
  x_scale = x_max - x_min;

  x_max = x_max + 0.05 * x_scale;
  x_min = x_min - 0.05 * x_scale;
  x_scale = x_max - x_min;

  y_max = max ( max ( node_xy(2,1:node_num) ), ...
                max ( point_xy(2,1:point_num) ) );
  y_min = min ( min ( node_xy(2,1:node_num) ), ...
                min ( point_xy(2,1:point_num) ) );
  y_scale = y_max - y_min;

  y_max = y_max + 0.05 * y_scale;
  y_min = y_min - 0.05 * y_scale;
  y_scale = y_max - y_min;

  if ( x_scale < y_scale )

    delta = round ( ( x_ps_max - x_ps_min ) ...
      * ( y_scale - x_scale ) / ( 2.0 * y_scale ) );

    x_ps_max = x_ps_max - delta;
    x_ps_min = x_ps_min + delta;

    x_ps_max_clip = x_ps_max_clip - delta;
    x_ps_min_clip = x_ps_min_clip + delta;

    x_scale = y_scale;

  elseif ( y_scale < x_scale )

    delta = round ( ( y_ps_max - y_ps_min ) ...
      * ( x_scale - y_scale ) / ( 2.0 * x_scale ) );

    y_ps_max      = y_ps_max - delta;
    y_ps_min      = y_ps_min + delta;

    y_ps_max_clip = y_ps_max_clip - delta;
    y_ps_min_clip = y_ps_min_clip + delta;

    y_scale = x_scale;

  end

  file_unit = fopen ( file_name, 'wt' );

  if ( file_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGLE_POINTS_PLOT - Fatal error!\n' );
    fprintf ( 1, '  Can not open output file.\n' );
    return
  end

  fprintf ( file_unit, '%!PS-Adobe-3.0 EPSF-3.0\n' );
  fprintf ( file_unit, '%%Creator: triangulation_order3_plot.f90\n' );
  fprintf ( file_unit, '%%Title: %s\n', file_name );
  fprintf ( file_unit, '%%Pages: 1\n' );

  fprintf ( file_unit, '%%BoundingBox:  %4d  %4d  %4d  %4d\n', ...
    x_ps_min, y_ps_min, x_ps_max, y_ps_max );
  fprintf ( file_unit, '%%Document-Fonts: Times-Roman\n' );
  fprintf ( file_unit, '%%LanguageLevel: 1\n' );
  fprintf ( file_unit, '%%EndComments\n' );
  fprintf ( file_unit, '%%BeginProlog\n' );
  fprintf ( file_unit, '/inch {72 mul} def\n' );
  fprintf ( file_unit, '%%EndProlog\n' );
  fprintf ( file_unit, '%%Page: 1 1\n' );
  fprintf ( file_unit, 'save\n' );
  fprintf ( file_unit, '%\n' );
  fprintf ( file_unit, '%  Set the RGB line color to very light gray.\n' );
  fprintf ( file_unit, '%\n' );
  fprintf ( file_unit, '0.900  0.900  0.900 setrgbcolor\n' );
  fprintf ( file_unit, '%\n' );
  fprintf ( file_unit, '%  Draw a gray border around the page.\n' );
  fprintf ( file_unit, '%\n' );
  fprintf ( file_unit, 'newpath\n' );
  fprintf ( file_unit, '  %4d  %4d  moveto\n', x_ps_min, y_ps_min );
  fprintf ( file_unit, '  %4d  %4d  lineto\n', x_ps_max, y_ps_min );
  fprintf ( file_unit, '  %4d  %4d  lineto\n', x_ps_max, y_ps_max );
  fprintf ( file_unit, '  %4d  %4d  lineto\n', x_ps_min, y_ps_max );
  fprintf ( file_unit, '  %4d  %4d  lineto\n', x_ps_min, y_ps_min );
  fprintf ( file_unit, 'stroke\n' );
  fprintf ( file_unit, '%\n' );
  fprintf ( file_unit, '%  Set the RGB color to black.\n' );
  fprintf ( file_unit, '%\n' );
  fprintf ( file_unit, '0.000  0.000  0.000 setrgbcolor\n' );
  fprintf ( file_unit, '%\n' );
  fprintf ( file_unit, '%  Set the font and its size.\n' );
  fprintf ( file_unit, '%\n' );
  fprintf ( file_unit, '/Times-Roman findfont\n' );
  fprintf ( file_unit, '0.50 inch scalefont\n' );
  fprintf ( file_unit, 'setfont\n' );
  fprintf ( file_unit, '%\n' );
  fprintf ( file_unit, '%  Print a title.\n' );
  fprintf ( file_unit, '%\n' );
  fprintf ( file_unit, '%  210  702  moveto\n' );
  fprintf ( file_unit, '%  (Triangulation)  show\n' );
  fprintf ( file_unit, '%\n' );
  fprintf ( file_unit, '%  Define a clipping polygon.\n' );
  fprintf ( file_unit, '%\n' );
  fprintf ( file_unit, 'newpath\n' );
  fprintf ( file_unit, '  %4d  %4d  moveto\n', x_ps_min_clip, y_ps_min_clip );
  fprintf ( file_unit, '  %4d  %4d  lineto\n', x_ps_max_clip, y_ps_min_clip );
  fprintf ( file_unit, '  %4d  %4d  lineto\n', x_ps_max_clip, y_ps_max_clip );
  fprintf ( file_unit, '  %4d  %4d  lineto\n', x_ps_min_clip, y_ps_max_clip );
  fprintf ( file_unit, '  %4d  %4d  lineto\n', x_ps_min_clip, y_ps_min_clip );
  fprintf ( file_unit, 'clip newpath\n' );
%
%  Draw the nodes.
%
  if ( 1 <= node_show )

    circle_size = 5;

    fprintf ( file_unit, '%\n' );
    fprintf ( file_unit, '%  Draw filled dots at the nodes.\n' );
    fprintf ( file_unit, '%\n' );
    fprintf ( file_unit, '%  Set the RGB color to blue.\n' );
    fprintf ( file_unit, '%\n' );
    fprintf ( file_unit, '0.000  0.150  0.750 setrgbcolor\n' );
    fprintf ( file_unit, '%\n' );

    for node = 1 : 3

      x_ps = round ( ...
        ( ( x_max - node_xy(1,node)         ) * x_ps_min   ...
        + (         node_xy(1,node) - x_min ) * x_ps_max ) ...
        / ( x_max                   - x_min ) );

      y_ps = round ( ...
        ( ( y_max - node_xy(2,node)         ) * y_ps_min   ...
        + (         node_xy(2,node) - y_min ) * y_ps_max ) ...
        / ( y_max                   - y_min ) );

      fprintf ( file_unit, 'newpath  %4d  %4d  %4d  0 360 arc closepath fill\n',...
        x_ps, y_ps, circle_size );

    end

  end
%
%  Label the nodes.
%
  if ( 2 <= node_show )

    fprintf ( file_unit, '%\n' );
    fprintf ( file_unit, '%  Label the nodes:\n' );
    fprintf ( file_unit, '%\n' );
    fprintf ( file_unit, '%  Set the RGB color to darker blue.\n' );
    fprintf ( file_unit, '%\n' );
    fprintf ( file_unit, '0.000  0.250  0.850 setrgbcolor\n' );
    fprintf ( file_unit, '/Times-Roman findfont\n' );
    fprintf ( file_unit, '0.20 inch scalefont\n' );
    fprintf ( file_unit, 'setfont\n' );
    fprintf ( file_unit, '%\n' );

    for node = 1 : node_num

      x_ps = round ( ...
        ( ( x_max - node_xy(1,node)         ) * x_ps_min   ...
        + (       + node_xy(1,node) - x_min ) * x_ps_max ) ...
        / ( x_max                   - x_min ) );

      y_ps = round ( ...
        ( ( y_max - node_xy(2,node)         ) * y_ps_min   ...
        + (         node_xy(2,node) - y_min ) * y_ps_max ) ...
        / ( y_max                   - y_min ) );

      fprintf ( file_unit, '%4d  %4d  moveto ( %4d ) show\n', x_ps, y_ps+5, node );

    end

  end
%
%  Draw the points.
%
  if ( point_num <= 200 )
    circle_size = 5;
  elseif ( point_num <= 500 )
    circle_size = 4;
  elseif ( point_num <= 1000 )
    circle_size = 3;
  elseif ( point_num <= 5000 )
    circle_size = 2;
  else
    circle_size = 1;
  end

  if ( 1 <= point_show )

    fprintf ( file_unit, '%\n' );
    fprintf ( file_unit, '%  Draw filled dots at the points.\n' );
    fprintf ( file_unit, '%\n' );
    fprintf ( file_unit, '%  Set the RGB color to green.\n' );
    fprintf ( file_unit, '%\n' );
    fprintf ( file_unit, '0.150  0.750  0.000 setrgbcolor\n' );
    fprintf ( file_unit, '%\n' );

    for point = 1 : point_num

      x_ps = round ( ...
        ( ( x_max - point_xy(1,point)         ) * x_ps_min   ...
        + (         point_xy(1,point) - x_min ) * x_ps_max ) ...
        / ( x_max                     - x_min ) );

      y_ps = round ( ...
        ( ( y_max - point_xy(2,point)         ) * y_ps_min   ...
        + (         point_xy(2,point) - y_min ) * y_ps_max ) ...
        / ( y_max                     - y_min ) );

      fprintf ( file_unit, 'newpath  %4d  %4d  %4d  0 360 arc closepath fill\n', ...
        x_ps, y_ps, circle_size );

    end

  end
%
%  Label the points.
%
  if ( 2 <= point_show )

    fprintf ( file_unit, '%\n' );
    fprintf ( file_unit, '%  Label the point:\n' );
    fprintf ( file_unit, '%\n' );
    fprintf ( file_unit, '%  Set the RGB color to darker green.\n' );
    fprintf ( file_unit, '%\n' );
    fprintf ( file_unit, '0.250  0.850  0.000 setrgbcolor\n' );
    fprintf ( file_unit, '/Times-Roman findfont\n' );
    fprintf ( file_unit, '0.20 inch scalefont\n' );
    fprintf ( file_unit, 'setfont\n' );
    fprintf ( file_unit, '%\n' );

    for point = 1 : point_num

      x_ps = round ( ...
        ( ( x_max - point_xy(1,point)         ) * x_ps_min   ...
        + (       + point_xy(1,point) - x_min ) * x_ps_max ) ...
        / ( x_max                     - x_min ) );

      y_ps = round ( ...
        ( ( y_max - point_xy(2,point)         ) * y_ps_min   ...
        + (         point_xy(2,point) - y_min ) * y_ps_max ) ...
        / ( y_max                     - y_min ) );

      fprintf ( file_unit, '%4d  %4d  moveto ( %4d ) show\n', x_ps, y_ps+5, point );

    end

  end
%
%  Draw the triangle.
%
  if ( 0 <= node_show )

    fprintf ( file_unit, '%\n' );
    fprintf ( file_unit, '%  Set the RGB color to red.\n' );
    fprintf ( file_unit, '%\n' );
    fprintf ( file_unit, '0.900  0.200  0.100 setrgbcolor\n' );
    fprintf ( file_unit, '%\n' );
    fprintf ( file_unit, '%  Draw the triangle.\n' );
    fprintf ( file_unit, '%\n' );

    fprintf ( file_unit, 'newpath\n' );

    for i = 1 : 4

      node = i4_wrap ( i, 1, 3 );

      x_ps = round ( ...
        ( ( x_max - node_xy(1,node)         ) * x_ps_min   ...
        + (         node_xy(1,node) - x_min ) * x_ps_max ) ...
        / ( x_max                   - x_min ) );

      y_ps = round ( ...
        ( ( y_max - node_xy(2,node)         ) * y_ps_min   ...
        + (         node_xy(2,node) - y_min ) * y_ps_max ) ...
        / ( y_max                   - y_min ) );

      if ( i == 1 )
        fprintf ( file_unit, '%3d  %3d  moveto\n', x_ps, y_ps );
      else
        fprintf ( file_unit, '%3d  %3d  lineto\n', x_ps, y_ps );
      end

    end

    fprintf ( file_unit, 'stroke\n' );

  end

  fprintf ( file_unit, '%\n' );
  fprintf ( file_unit, 'restore  showpage\n' );
  fprintf ( file_unit, '%\n' );
  fprintf ( file_unit, '%  End of page.\n' );
  fprintf ( file_unit, '%\n' );
  fprintf ( file_unit, '%%Trailer\n' );
  fprintf ( file_unit, '%%EOF\n' );

  fclose ( file_unit );

  return
end


text:
function this = parse_input_parameter(this, param)
%PARSE_INPUT_PARAMETER parses and set parameters in the struct PARAM
%
%   This is a member function of the class 'alt_ba_optical_flow'.
%
% Authors: Deqing Sun, Department of Computer Science, Brown University
% Contact: dqsun@cs.brown.edu
% $Date: 2009 $
% $Revision: $
%
% Copyright 2009-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

if length(param) ==1
    param = param{1};
end;

if mod(length(param), 2) ~=0
    error('Parse_input_parameter: Input parameters must be given in pairs (name and value)!');
end;

i = 1;

while (i <= length(param))

    if ischar(param{i+1})
        param{i+1} = str2num(param{i+1});
    end;

    switch lower(param{i})

        case 'lambda'
            this.lambda         = param{i+1};
            this.lambda_q       = param{i+1};

        case 'pyramid_levels'
            this.pyramid_levels = param{i+1};
            this.auto_level     = false;

        case 'pyramid_spacing'
            this.pyramid_spacing        = param{i+1};

        case 'gnc_pyramid_levels'
            this.gnc_pyramid_levels     = param{i+1};

        case 'gnc_pyramid_spacing'
            this.gnc_pyramid_spacing    = param{i+1};

        case 'sigma_d'
            for j = 1:length(this.spatial_filters);
                this.rho_spatial_u{j}.param   = param{i+1};
                this.rho_spatial_v{j}.param   = param{i+1};
            end;

        case 'sigma_s'
            this.rho_data.param         = param{i+1};

    end;

    i = i+2;
end;

author: kristinbranson

repo: JAABA


text:
% goes through all matlab files in the folder
% loads file
% saves as CSV file

base_folder_smartmeter = 'data/eco/smartmeter/';
base_folder_plugs = 'data/eco/plugs/';
target_folder_smartmeter = 'data/eco/smartmeter_csv/';
target_folder_plugs = 'data/eco/plugs_csv/';
process_smart_meter = 1;
process_plugs = 0;

% household = [ 2:6 ];
% plugs = [ 1:12 ];
household = [ 4 ];
plugs = [ 3, 8 ];

for h = household

    houseStr = num2str(h, '%02d');

	%% smart meter
    if process_smart_meter == 1
        source_folder = [base_folder_smartmeter, houseStr, '/'];
        target_folder = [target_folder_smartmeter, houseStr, '/'];
        if ~exist(target_folder, 'dir')
            mkdir(target_folder);
        end
        files = dir([source_folder, '*.mat']);
        for day = 1:size(files,1)
            fprintf('Processing smart meter household %s, day %d of %d\n', houseStr, day, size(files,1));
            filename = files(day,1).name;
            source_file = [source_folder, filename];
            vars = whos('-file',source_file);
            load(source_file);
            eval(['smartmeter_data=' vars.name ';']);
            eval(['clear ' vars.name ';']);
            data_to_write = [ ...
                smartmeter_data.powerallphases ...
                smartmeter_data.powerl1 ...
                smartmeter_data.powerl2 ...
                smartmeter_data.powerl3 ...
                smartmeter_data.currentneutral ...
                smartmeter_data.currentl1 ...
                smartmeter_data.currentl2 ...
                smartmeter_data.currentl3 ...
                smartmeter_data.voltagel1 ...
                smartmeter_data.voltagel2 ...
                smartmeter_data.voltagel3 ...
                smartmeter_data.phaseanglevoltagel2l1 ...
                smartmeter_data.phaseanglevoltagel3l1 ...
                smartmeter_data.phaseanglecurrentvoltagel1 ...
                smartmeter_data.phaseanglecurrentvoltagel2 ...
                smartmeter_data.phaseanglecurrentvoltagel3 ...
            ];
            dlmwrite([target_folder, filename(1:end-4), '.csv'], data_to_write, 'precision', 9);
        end
    end

    if process_plugs == 1
        for plug = plugs
            plugStr = num2str(plug, '%02d');
            plug_dir = [ base_folder_plugs, houseStr, '/', plugStr, '/'];
            if ~exist(plug_dir, 'dir')
                continue;
            end
            target_folder = [target_folder_plugs, houseStr, '/', plugStr, '/'];
            if ~exist(target_folder, 'dir')
                mkdir(target_folder)
            end
            files = dir([plug_dir, '*.mat']);
            for day = 1:size(files,1)
                fprintf('Processing plug %s for household %s, day %d of %d\n', plugStr, houseStr, day, size(files,1));
                filename = files(day,1).name;
                source_file = [plug_dir, filename];
                vars = whos('-file',source_file);
                load(source_file);
                eval(['plug_data=' vars.name ';']);
                eval(['clear ' vars.name ';']);
                data_to_write = plug_data.consumption;
                dlmwrite([target_folder, filename(1:end-4), '.csv'], data_to_write, 'precision', 9);
            end
        end
    end
end

author: beckel

repo: nilm-eval

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


text:
function [img,mp] = wedgeMap(mp,diameter,totAngle)
%
%  [img,mp] = wedgeMap(mp,diameter,totAngle)
%
%AUTHOR:  Baseler, Wandell
%DATE:    07.23.96
%PURPOSE:
%  Make one of those darn circular color maps to
% represent the polar angle
%
%ARGUMENTS:
%
% mp:  A polar angle color map, the last value should
%      be the background color
%            (Default = [ hsv(64) ; .5 .5 .5])
% diameter:  The size of the returned image.
%            (Default 128)
% totAngle:  The angle around the circle the map
%    should sweep out.  [0 is the noon position, and
%    angles run clockwise]
%            (Default = 180)
%RETURNS:
%
% img:  The resulting image
% mp:   The colormap of the image
%    To display the result use image(img),colormap(mp)
%

% Test parameter
% diameter = 256;
% mp = [hsv(128); 0.5 0.5 0.5];
% totAngle = 360;

if nargin < 3,  totAngle = 180; end
if nargin < 2,  diameter = 128; end
if nargin < 1,   mp = [hsv(64); .5 .5 .5]; end

% Input arguments
%
nX = diameter;
nY = diameter(1);
nMap = size(mp,1);

% Create a grid of (X,Y) values
%
[X Y] = meshgrid(1:nX,1:nY);

% Center the grid around (0,0)
%
X = X - (nX/2); Y = Y - (nY/2);

%  Find the angle for each of the X,Y points
%
ang = zeros(size(X));
posX = (X > 0);
ang(posX) = atan(Y(posX) ./ (X(posX) + eps))*360/(2*pi);
ang = ang + 90;
tmp = rot90(ang + 180,2);
ang(~posX) = tmp(~posX);

dist = sqrt(X.^2 + Y.^2);
radius = diameter/2;

% Pick out those locations that are within a radius
% and that are on the right side of the axis (i.e., X>0)
%
inRadius = (dist < radius);
inAngle = ( ang < totAngle);

% Scale angles to sweep out the color map except
% for the last entry (which is used for background)
%
ang(inAngle) = mrScale(ang(inAngle),1,(nMap - 1));

% Find positions NOT in the selected range
%
ang(~inRadius) = nMap*ones(size(ang(~inRadius)));
ang(~inAngle) =  nMap*ones(size(ang(~inAngle)));
%imagesc(ang)
img = ang;

% Make an image
%
% figure(1)
% clf
% hold off
% colormap(mp)
% imagesc(img), axis image, axis off
%
% tiffwrite(ang,mp,'rotateMap.tif');
% unix('xv rotateMap.tif &')

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%%% Note: run the 'GenerateData_model_64_25_Res_Bnorm_Adam.m' to generate
%%% training data first.

%%%-------------------------------------------------------------------------
%%% configuration
%%%-------------------------------------------------------------------------
opts.modelName        = 'model_64_25_Res_Bnorm_Adam'; %%% model name
opts.learningRate     = [logspace(-3,-3,30) logspace(-4,-4,20)];%%% you can change the learning rate
opts.batchSize        = 128; %%% default
opts.gpus             = [1]; %%% this code can only support one GPU!

%%% solver
opts.solver           = 'Adam';

opts.gradientClipping = false; %%% Set 'true' to prevent exploding gradients in the beginning.
opts.expDir      = fullfile('data', opts.modelName);
opts.imdbPath    = fullfile(opts.expDir, 'imdb.mat');

%%%-------------------------------------------------------------------------
%%%   Initialize model and load data
%%%-------------------------------------------------------------------------
%%%  model
net  = feval(['DnCNN_init_',opts.modelName]);

%%%  load data
imdb = load(opts.imdbPath) ;

%%%-------------------------------------------------------------------------
%%%   Train
%%%-------------------------------------------------------------------------

[net, info] = DnCNN_train(net, imdb, ...
    'expDir', opts.expDir, ...
    'learningRate',opts.learningRate, ...
    'solver',opts.solver, ...
    'gradientClipping',opts.gradientClipping, ...
    'batchSize', opts.batchSize, ...
    'modelname', opts.modelName, ...
    'gpus',opts.gpus) ;

author: cszn

repo: DnCNN

save_path: github-repos/MATLAB/cszn-DnCNN


text:
clear all; clc; close all;

OpenBMI('C:\Users\CVPR\Desktop\Open_Github') % Edit the variable BMI if necessary
global BMI;
BMI.EEG_DIR=['G:\data2'];

%% DATA LOAD MODULE
file=fullfile(BMI.EEG_DIR, '\2016_08_05_hkkim_training');
marker= {'1','right';'2','left';'3','foot';'4','rest'};
[EEG.data, EEG.marker, EEG.info]=Load_EEG(file,{'device','brainVision';'marker', marker;'fs', 500});

field={'x','t','fs','y_dec','y_logic','y_class','class', 'chan'};
CNT=opt_eegStruct({EEG.data, EEG.marker, EEG.info}, field);
CNT=prep_selectClass(CNT,{'class',{'right', 'rest'}});

CNT2 = prep_laplacian(CNT, {'Channel', {'C3', 'Cz', 'C4'}})

%% CROSS-VALIDATION MODULE
CV.var.band=[7 20];
CV.var.interval=[750 3500];
CV.prep={ % commoly applied to training and test data before data split
    'CNT=prep_filter(CNT, {"frequency", band})'
    'SMT=prep_segmentation(CNT, {"interval", interval})'
    };
CV.train={
    '[SMT, CSP_W, CSP_D]=func_csp(SMT,{"nPatterns", [3]})'
    'FT=func_featureExtraction(SMT, {"feature","logvar"})'
    '[CF_PARAM]=func_train(FT,{"classifier","LDA"})'
    };
CV.test={
    'SMT=func_projection(SMT, CSP_W)'
    'FT=func_featureExtraction(SMT, {"feature","logvar"})'
    '[cf_out]=func_predict(FT, CF_PARAM)'
    };
CV.option={
'KFold','7'
% 'leaveout'
};

[loss]=eval_crossValidation(CNT, CV); % input : eeg, or eeg_epo

author: PatternRecognition

repo: OpenBMI


text:
function Datacube = Pooling(ktensor)
% input
%   ktensor(Cell Array): Each cell contains a ktensor caculated by a given
%   channel of convolutionized Datacube
% output
%   Datacube(Cell Array): For each cell of ktensor, discarding the third
%   dimension, outproduct the first and the second dimension with
%   corresponding lambda to form a Datacube with equal size of space
%   Maxpooling ?
num_channels = size(ktensor,1);
Datacube = cell(num_channels,1);
for i=1:num_channels
    num_of_U = size(ktensor{i}.U{1},2);
    combined_basis = [];
    for j=1:num_of_U
        temp = ktensor{i}.lambda(j) * ktensor{i}.U{1}(:,j) * ktensor{i}.U{2}(:,j)';
        combined_basis(:,:,j) = temp;
    end
    Datacube{i} = combined_basis;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [z,dz,ymu,ys,fmu,fs,fpi] = acqNegSqEI(xi,target,gpstruct,optimState,grad_flag)
%ACQNEGEI Acquisition function for (negative) expected squared improvement.

if nargin < 5 || isempty(grad_flag); grad_flag = false; end

n = size(xi,1);

if grad_flag && n > 1
    error('acqNegSqEI:gradient', ...
        'Gradient of acquisition function is provided only at one test point XI (row vector).');
end

if grad_flag
    [ymu,ys2,fmu,fs2,hypw,dymu,dys2,dfmu,dfs2] = gppred(xi,gpstruct,'central');
else
    [ymu,ys2,fmu,fs2,hypw] = gppred(xi,gpstruct);
end
fs = sqrt(fs2);
ys = sqrt(ys2);

% Probability of improvement
gammaz = real((target - fmu)./fs);
fpi = 0.5*erfc(-gammaz/sqrt(2));

% Expected squared improvement
z = -fs.^2.*((gammaz.^2+1).*fpi + gammaz.*exp(-0.5*(gammaz.^2))/sqrt(2*pi));

try
    z = sum(bsxfun(@times,hypw(~isnan(hypw)),z(~isnan(hypw),:)),1);
catch
    z = Inf(1,n);
    dz = NaN(n,size(xi,2));
    return;
end

if grad_flag
    % Gradient of probability of improvement
    dfs = 0.5*dfs2./fs;
    dgammaz = -(dfmu.*fs + (target - fmu).*dfs)./fs2;
    dfpi = -0.5*dgammaz/sqrt(2)*(-2*exp(-gammaz.^2/2)/sqrt(pi));

    % Gradient of expected improvement
    %dz = -(dfs.*gammaz.*fpi + dgammaz.*fs.*fpi + dfpi.*gammaz.*fs) ...
    %    + (fs.*gammaz.*dgammaz - dfs).*exp(-0.5*gammaz.^2)/sqrt(2*pi);
    %dz = sum(bsxfun(@times,hypw,dz(~isnan(hypw),:)),1);
else
    dz = NaN(n,size(xi,2));     % Gradient not estimated
end

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function y = pow_p( x, p )

%POW_P   Positive branch of the power function.
%   POW_P(X,P) computes a convex or concave branch of the power function:
%           P < 0: POW_P(X,P) = X.^P if X >  0, +Inf otherwise
%      0 <= P < 1: POW_P(X,P) = X.^P if X >= 0, -Inf otherwise
%      1 <= P    : POW_P(X,P) = X.^P if X >= 0, +Inf otherwise
%   Both P and X must be real.
%
%   Disciplined convex programming information:
%       The geometry of POW_P(X,P) depends on the precise value of P,
%       which must be a real constant:
%                P < 0: convex  and nonincreasing; X must be concave.
%           0 <= P < 1: concave and nondecreasing; X must be concave.
%           1 <= P    : convex  and nonmonotonic;  X must be affine.
%       In all cases, X must be real.

y = power( pdom( x ), p );

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
%{
 * Copyright (C) 2013-2020, The Regents of The University of Michigan.
 * All rights reserved.
 * This software was developed in the Biped Lab (https://www.biped.solutions/)
 * under the direction of Jessy Grizzle, grizzle@umich.edu. This software may
 * be available under alternative licensing terms; contact the address above.
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 * The views and conclusions contained in the software and documentation are those
 * of the authors and should not be interpreted as representing official policies,
 * either expressed or implied, of the Regents of The University of Michigan.
 *
 * AUTHOR: Bruce JK Huang (bjhuang[at]umich.edu)
 * WEBSITE: https://www.brucerobot.com/
%}

function cost = verifyCornerAccuracy(X_verification, Y_verification, P)
    C_X_transformed = P * X_verification;
    C_X_transformed = C_X_transformed ./ C_X_transformed(3,:);
    cost = (norm(C_X_transformed(1:2,:) - Y_verification(1:2,:), 'fro'))^2;
end

author: UMich-BipedLab


text:
function [xH, H, HNMax, HNMin] = hopkins_main (dataXY_all, xlim1, xlim2, ylim1, ylim2, m, N, Nbins, envelopes, Nsimul, filename)
% HOPKINS_MAIN computes hopkins statistics
% [xH, H, HNMax, HNMin] = hopkins_main (dataXY_all, xlim1, xlim2, ylim1, ylim2, m, N, Nbins, envelopes, Nsimul, filename)
% dataXY_all - whole data
% [xlim1, xlim2, ylim1, ylim2] - selected ROI
% m - number of random events and points for Hopkins computation
% N - number of iterations for histogram computation
% Nbins - numbero of bins in histogram
% envelopes - if set to 1 maximum and minimum envelopes for Poisson random
% process are computed
% Nsimul - number of simulations for envelopes computation
% filename - name of the txt file where results are written

Hopp.xlim1=xlim1; Hopp.xlim2=xlim2; Hopp.ylim1=ylim1; Hopp.ylim2=ylim2;
Hopp.Nbins = Nbins; Hopp.m = m; Hopp.N = N;
if envelopes, Hopp.Nsimul = Nsimul; end

figure
dataXY = ROIdata(dataXY_all, Hopp.xlim1, Hopp.xlim2, Hopp.ylim1, Hopp.ylim2, 1);
Hopp.Npoints = length(dataXY);

fprintf('Computing Hopkins statistics for selected ROI... \n');
h = waitbar(0,'Computing Hopkins...');
for ii=1:Hopp.N
    HopD(ii) = hopkins(dataXY, Hopp.xlim1, Hopp.xlim2, Hopp.ylim1, Hopp.ylim2, Hopp.m);
    waitbar(ii/Hopp.N,h)
end
close (h)

[H, xH] = hist(HopD, Nbins);
H = H/trapz(xH, H); % normalization

if envelopes
    [HNMin, HNMax, xH] = ...
        hopkinsEnv(Hopp.Npoints, Hopp.xlim1, Hopp.xlim2, Hopp.ylim1, Hopp.ylim2, Hopp.m,Hopp.N, xH, Hopp.Nsimul);
    if nargin>10
        writedata(xH, HNMax, Hopp,  [filename '_envelope_Max'], ['Hopkins statistics - Maximum of ' num2str(Hopp.Nsimul) ' simulations of Poisson process']);
        writedata(xH, HNMin, Hopp,  [filename '_envelope_Min'], ['Hopkins statistics - Mainimum of ' num2str(Hopp.Nsimul) ' simulations of Poisson process'])
    end
else
    HNMin = [];
    HNMax = [];
end

if nargin>10
    writedata(xH, H, Hopp,  filename, 'Hopkins statistics')
    save (filename)
end

function [HhistMin, HhistMax, xHhist] = hopkinsEnv(Npoints, xlim1, xlim2, ylim1, ylim2, m,N, bins, Nsimul)
% function [HhistMin, HhistMax, xHhist] = hopkinsEnv(Npoints, xlim1, xlim2,
% ylim1, ylim2, m,N, bins, Nsimul)

h = waitbar(0,'Computing Hopkins Envelope...');
for ii=1 : Nsimul
    waitbar(ii/Nsimul,h)
    simNoise = generateNoise(Npoints, [xlim1,xlim2,ylim1,ylim2]);
    for jj=1 : N
        H(jj) = hopkins (simNoise,xlim1, xlim2, ylim1, ylim2, m);
    end
    [Hhist, xHhist] = hist(H, bins);
    Hhist = Hhist/trapz(xHhist, Hhist); % normalization
    if ii==1
        HhistMin = Hhist;
        HhistMax = Hhist;
    else
        HhistMin = min(HhistMin, Hhist);
        HhistMax = max(HhistMax, Hhist);
    end
end
close (h)

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function get_sofea()
    disp( [' The restrictive policy of the mathworks File Exchange prevents me from posting the actual source'])
    disp( [' Please visit:'])
    disp( ['     hogwarts.ucsd.edu/~pkrysl/sofea']);
    disp( [' Petr Krysl'])
end


text:
function pde = jumpmgdata2
%% Data of JUMPMG1
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

pde = struct('f',@f,'g_D',@g_D,'d',@omega);

    function s = f(p) % load data (right hand side function)
    s = zeros(size(p,1),1);
    end

    function s = g_D(p) % Dirichlet boundary condition
    s = zeros(size(p,1),1);
    x = p(:,1);
    idx = (abs(x-1)<eps);
    s(idx) = 1;
    s(~idx) = 0;
    end

    function c = omega(p) % diffusion constant
    global epsilon
    c = zeros(size(p,1),1);
    x = p(:,1); y = p(:,2); z = p(:,3);
    idx = ((x>0) & (x<0.5) & (y>0) & (y<0.5) & (z>0) & (z<0.5)) | ...
          ((x<0) & (x>-0.5) & (y<0) & (y>-0.5) & (z<0) & (z>-0.5));
    c(idx) = 1;
    c(~idx) = epsilon;
    end
end

repo: ifem


text:
 function out = cpu(arg, varargin)
%function out = cpu(arg, varargin)
%|
%| cpu etic
%| cpu etoc printarg
%|	tic/toc elapsed time
%|
%| cpu etic
%| cpu tic
%| cpu toc
%| cpu toc printarg
%| work like tic/toc except they use cpu time instead of wall time.
%|
%| also:
%| If printarg begins with a ':' then caller_name preceeds when printing
%|
%| for multiple levels of clock use:
%|	c = cpu('etic');
%|	(run some slow commands)
%|	cpu(c, printarg)

if nargin < 1, ir_usage, end
if streq(arg, 'test'), cpu_test, return, end

% todo: add default behaviour

persistent t
if ~isvar('t') || isempty(t)
	t.tic = [];
	t.clock = [];
end

if isnumeric(arg) && numel(arg) == 6
	out = etime(clock, arg);

elseif streq(arg, 'tic')
	t.tic = cputime;
	if length(varargin), error 'tic takes no option', end
	return

elseif streq(arg, 'etic')
	if length(varargin), error 'tic takes no option', end
	if nargout == 1
		out = clock;
	else
		t.clock = clock;
	end
	return

elseif streq(arg, 'toc')
	if isempty(t.tic)
		error 'must initialize cpu with tic first'
	end
	out = cputime - t.tic;

elseif streq(arg, 'etoc')
	if isempty(t.clock)
		error 'must initialize cpu with etic first'
	end
	out = etime(clock, t.clock);

else
	error 'bad argument'
end

if length(varargin) == 1
	arg = varargin{1};
	if ischar(arg)
		if arg(1) == ':'
			[name line] = caller_name;
			if ~isempty(line) && line ~= 0
				name = [name sprintf(' %d', line)];
			end
			printf('%s%s %g', name, arg, out)
		else
			printf('%s %g', arg, out)
		end
	end

	if ~nargout
		clear out
	end
elseif length(varargin) > 1
	error 'too many arguments'
end

function cpu_test
cpu tic
cpu('toc');
cpu toc cpu_test:toc
cpu etic
cpu('etoc');
cpu etoc cpu_test:etoc
out = cpu('etoc', 'cpu_test:out');

tmp = cpu('etic');
pause(0.1)
cpu(tmp, 'ok')

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function pb = bsdsNonMaxSuppression(pbnew, pbold)

if size(pbold)~=size(pbnew)
    pbold = imresize(pbold, size(pbnew), 'nearest');
end
pb = pbnew.*double(pbold>0);

% [h, w] = size(pbnew);
% norient = size(pball, 3);
% theta = (0:norient-1)/norient*pi;
%
% [h2, w2, tmp] = size(pball);
% if h2~=h || w~=w2
%     pball = imresize(pball, [h w], 'nearest');
% end
%
% % nonmax suppression and max over orientations
% [unused,maxo] = max(pball,[],3);
% pb = zeros(h,w);
% %theta = zeros(h,w);
% r = 2.5;
% for i = 1:norient,
%   mask = (maxo == i);
%   %a = fitparab(pball(:,:,i),r,r,theta(i));
%   %pbi = nonmax(max(0,a),gtheta(i));
%   pbi = nonmax(pbnew,theta(i));
%   pb = max(pb,pbi.*mask);
% end
% pb = max(0,min(1,pb));
%
% % mask out 1-pixel border where nonmax suppression fails
% pb(1,:) = 0;
% pb(end,:) = 0;
% pb(:,1) = 0;
% pb(:,end) = 0;

author: Cloud-CV

repo: object-proposals


text:
function [varargout]=edgeVec(E,V)

% function [N,Vp,Nv]=edgeVec(E,V)
% ------------------------------------------------------------------------
%
% ------------------------------------------------------------------------

%% Compute edge vectors
Vp=V(E(:,1),:); %Edge vector origin
N=V(E(:,2),:)-Vp; %Edge vector

%% Collect output
varargout{1}=N;
varargout{2}=Vp;

if nargout==3
    Nv=faceToVertexMeasure(E,V,N); %Edge vectors at vertices
    varargout{3}=Nv;
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.


text:
function [ r, g, b ] = ppma_read ( file_name )

%*****************************************************************************80
%
%% PPMA_READ opens an ASCII PPM file and reads the data.
%
%  The PPMA file is assumed to have the format:
%
%    P3
%    # feep.ppma
%    ncol nrow
%    maxrgb, the largest legal value for the data.
%    row 1, a list of ncol triples between 0 and maxrgb.
%    row 2, a list of ncol triples between 0 and maxrgb.
%    ...
%    row nrow, a list of ncol triples between 0 and maxrgb.
%
%    Lines beginning with '#' are comments.
%
%  Example:
%
%    P3
%    # feep.ppma
%    4 4
%    15
%     0  0  0  0  0  0    0  0  0   15  0 15
%     0  0  0  0 15  7    0  0  0    0  0  0
%     0  0  0  0  0  0    0 15  7    0  0  0
%    15  0 15  0  0  0    0  0  0    0  0  0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILE_NAME, the name of the PPMA file to read.
%
%    Output, integer R{NROW,NCOL], G{NROW,NCOL], B{NROW,NCOL] is the RGB data read from the file.
%
  FALSE = 0;
  TRUE = 1;

  r = [];
  g = [];
  b = [];

  fid = fopen ( file_name, 'r' );

  if ( fid < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PPMA_READ - Fatal error!\n' );
    fprintf ( 1, '  Could not open the input file.\n' );
    error ( 'PPMA_READ - Fatal error!' );
    return;
  end
%
%  Read the first line.
%
  line = fgets ( fid );
%
%  Verify that the first two characters are the "magic number".
%  Matlab strncmp returns 1 for equality, and 0 for inequality.
%
  if ( strncmp ( line, 'P3', 2 ) == 0 )
    return;
  end
%
%  Move to the next noncomment line.
%
  while ( 1 )
    line = fgets ( fid );
    if ( line(1) ~= '#' )
      break;
    end
  end
%
%  Extract NCOL and NROW.
%
  [ array, count ] = sscanf ( line, '%d' );
  ncol = array(1);
  nrow = array(2);
%
%  Move to the next noncomment line.
%
  while ( 1 )
    line = fgets ( fid );
    if ( line(1) ~= '#' )
      break;
    end
  end
%
%  Extract MAXRGB, and ignore it.
%
  [ array, count ] = sscanf ( line, '%d' );
  maxRGB = array(1);
%
%  Set aside space for R, G< B.
%
  r = zeros ( nrow, ncol );
  g = zeros ( nrow, ncol );
  b = zeros ( nrow, ncol );

  i = 1;
  j = 0;
  p = 0;
  done = FALSE;

  while ( done == FALSE )
%
%  Move to the next noncomment line.
%
    while ( 1 )
      line = fgets ( fid );
      if ( line(1) ~= '#' )
        break;
      end
    end
    [ array, count ] = sscanf ( line, '%d' );
%
%  Each value that you read goes into the "next" open entry.
%  When reading R, we update the row and column indices.
%
    for k = 1 : count

      if ( p == 0 )
        j = j + 1;
        if ( ncol < j )
          j = 1;
          i = i + 1;
        end
        r(i,j) = array(k);
      elseif ( p == 1 )
        g(i,j) = array(k);
      elseif ( p == 2 )
        b(i,j) = array(k);
        if ( i == nrow & j == ncol )
          done = TRUE;
        end
      end

      p = mod ( p+1, 3 );

    end

  end

  fclose ( fid );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function asa058_test01 ( )

%*****************************************************************************80
%
%% TEST01 tries out the ASA058 routine.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 February 2008
%
%  Author:
%
%    John Burkardt
%
  k = 5;
  m = 2;
  n = 100;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  Test the CLUSTR algorithm.\n' );
  fprintf ( 1, '  Applied Statistics Algorithm 58\n' );
%
%  Read the data.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Reading the data.\n' );

  input_unit = fopen ( 'points_100.txt', 'rt' );

  for i = 1 : n
    x(i,1:m) = fscanf ( input_unit, '%f', m );
  end

  fclose ( input_unit );
%
%  Print a few data values.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  First 5 data values:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : 5
    fprintf ( 1, '  %8d', i );
    for j = 1 : m
      fprintf ( 1, '  %14f', x(i,j) );
    end
    fprintf ( 1, '\n' );
  end
%
%  Initialize the cluster centers arbitrarily.
%
  for i = 1 : k
    for j = 1 : m
      d(i,j) = x(i,j);
    end
  end
%
%  Compute the clusters.
%
  nz = 1;
  k2 = k;

  [ d, dev, b, e ] = clustr ( x, d, n, m, k, nz, k2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Cluster  Population  Energy\n' );
  fprintf ( 1, '\n' );

  for i = 1 : k
    fprintf ( 1, '  %8d  %8d  %14f\n', i, e(i), dev(i) );
  end

  e_sum = sum ( e(1:k) );
  dev_sum = sum ( dev(1:k) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     Total  %8d  %14f\n', e_sum, dev_sum );

  return
end


text:
function gradfun = approxgradientFD(problem, options)
% Gradient approx. fnctn handle based on finite differences of the cost.
%
% function gradfun = approxgradientFD(problem)
% function gradfun = approxgradientFD(problem, options)
%
% Input:
%
% A Manopt problem structure (already containing the manifold and enough
% information to compute the cost) and an options structure (optional),
% containing one option:
%    options.stepsize (positive double; default: 2^-23).
%    options.subspacedim (positive integer; default: [], for M.dim()).
%
% If the cost cannot be computed on 'problem', a warning is issued.
%
% Output:
%
% Returns a function handle, encapsulating a generic finite difference
% approximation of the gradient of the problem cost. The finite difference
% is based on M.dim()+1 computations of the cost.
%
% The returned gradfun has this calling pattern:
%
%   function gradfd = gradfun(x)
%   function gradfd = gradfun(x, storedb)
%   function gradfd = gradfun(x, storedb, key)
%
% x is a point on the manifold problem.M, storedb is a StoreDB object,
% and key is the StoreDB key to point x.
%
% Usage:
%
% Typically, the user will set problem.M and other fields to define the
% cost (typically, problem.cost). Then, to use this generic purpose
% gradient approximation:
%
%   problem.approxgrad = approxgradientFD(problem, options);
%
% See also: steepestdescent conjugategradient

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Nov. 1, 2016.
% Contributors:
% Change log:

    % This gradient approximation is based on the cost:
    % check availability.
    if ~canGetCost(problem)
        warning('manopt:approxgradFD:nocost', ...
                'approxgradFD requires the cost to be computable.');
    end

    % Set local defaults here, and merge with user options, if any.
    localdefaults.stepsize = 2^-23;
    localdefaults.subspacedim = [];
    if ~exist('options', 'var') || isempty(options)
        options = struct();
    end
    options = mergeOptions(localdefaults, options);

    % % Finite-difference parameters
    % How far do we look?
    stepsize = options.stepsize;
    % Approximate the projection of the gradient on a random subspace of
    % what dimension? If [], uses full tangent space.
    subspacedim = options.subspacedim;

    % Build and return the function handle here. This extra construct via
    % funhandle makes it possible to make storedb and key optional.
    gradfun = @funhandle;
    function gradfd = funhandle(x, storedb, key)
        % Allow omission of the key, and even of storedb.
        if ~exist('key', 'var')
            if ~exist('storedb', 'var')
                storedb = StoreDB();
            end
            key = storedb.getNewKey();
        end
        gradfd = gradientFD(stepsize, subspacedim, problem, x, storedb, key);
    end

end

function gradfd = gradientFD(stepsize, subspacedim, problem, x, storedb, key)
% This function does the actual work.
%
% Original code: Nov. 1, 2016 (NB).

    % Evaluate the cost at the root point
    fx = getCost(problem, x, storedb, key);

    % Pick an orthonormal basis for the tangent space at x, or a subspace
    % thereof. The default is a full subspace. If a strict subspace is
    % picked, the returned vector approximates the orthogonal projection of
    % the gradient to that subspace.
    B = tangentorthobasis(problem.M, x, subspacedim);

    % Use finite differences to approximate the directional derivative
    % along each direction in the basis B.
    df = zeros(size(B));
    for k = 1 : numel(B)
        % Move in the B{k} direction
        xk = problem.M.retr(x, B{k}, stepsize);
        keyk = storedb.getNewKey();
        % Evaluate the cost there
        fxk = getCost(problem, xk, storedb, keyk);
        % Don't keep this point in cache
        storedb.remove(keyk);
        % Finite difference
        df(k) = (fxk - fx)/stepsize;
    end

    % Build the gradient approximation.
    gradfd = lincomb(problem.M, x, B, df);

end

author: NicolasBoumal

repo: manopt


text:
function [OutputError] = PIVlab_capture_OPTOcam_save(OPTOcam_vid,nr_of_images,ImagePath,frame_nr_display,bitmode)
if bitmode==8
	bitmultiplicator=1;
elseif bitmode==12
	bitmultiplicator = 16; %bring 12bit data to 16 bits full histogram, otherwise images outside Matlab are not displayed correctly (too dark).
end
hgui=getappdata(0,'hgui');
OutputError=0;
OPTOcam_frames_to_capture = nr_of_images*2;
if getappdata(hgui,'cancel_capture') ~=1 %capture was not cancelled --> save images from RAM to disk
	OPTOcam_data = getdata(OPTOcam_vid,OPTOcam_frames_to_capture); %ruft alle Frames in RAM ab. Frame 1,2,3 sind müll
	cntr=0;
	for image_save_number=1:2:size(OPTOcam_data,4)
		imgA_path=fullfile(ImagePath,['PIVlab_' sprintf('%4.4d',cntr) '_A.tif']);
		imgB_path=fullfile(ImagePath,['PIVlab_' sprintf('%4.4d',cntr) '_B.tif']);
		imwrite(OPTOcam_data(:,:,:,image_save_number)*bitmultiplicator,imgA_path,'compression','none'); %tif file saving seems to be the fastest method for saving data...
		imwrite(OPTOcam_data(:,:,:,image_save_number+1)*bitmultiplicator,imgB_path,'compression','none');
		cntr=cntr+1;
		set(frame_nr_display,'String',['Saving images to disk: Image pair ' num2str(cntr) ' of ' num2str(size(OPTOcam_data,4)/2)]);
		drawnow limitrate;
	end
end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function [output, validFrameMask] = F_hadamard(input_layers)
input = input_layers{1}.a;
input2 = input_layers{2}.a;

output = input .* input2;

[D,M,N] = size(input);
if N==1
    validFrameMask = [];
else
    validFrameMask = getValidFrameMask(input_layers{1});
end

end

author: singaxiong

repo: SignalGraph


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: SSD and Force Fields
%
%==============================================================================

clear, close all, help(mfilename);
sdiag = @(a) spdiags(reshape(a,[],1),0,length(a),length(a));
setup2DHNSPData;
level = 7; omega = ML{level}.omega; m = ML{level}.m; % load data
imgModel('reset','imgModel','splineInter')
viewImage('set','axis','off');
[T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega);
X = reshape(getCellCenteredGrid(omega,m),[],2);

[Tc,dT] = imgModel(T,omega,X);
[Rc,dR] = imgModel(R,omega,X);

axT   = [0.55,0.85,0.32,0.47];
axR   = [0.6,0.9,0.20,0.35];
axTF  = [0.55,0.70,0.32,0.395];

FAIRfigure(1,'color','w');
viewImage(Tc,omega,m); hold on; axis off;
ph = plot(axT([1,1,2,2,1]),axT([3,4,4,3,3]),'w-','linewidth',4);

set(ph,'visible','off');  axis(axT);

clf;
viewImage(Rc,omega,m); hold on; axis off;
ph = plot(axR([1,1,2,2,1]),axR([3,4,4,3,3]),'w-','linewidth',4);

set(ph,'visible','off');  axis(axR);

clf;
viewImage(128+0.5*(Tc-Rc),omega,m); hold on; axis off;
% ph = plot(axT([1,1,2,2,1]),axT([3,4,4,3,3]),'w-','linewidth',4);
% set(ph,'visible','off');  axis(axR);

X = getCellCenteredGrid(omega,m); X = reshape(X,[],2);
n = prod(m);
gradT     = spdiags(dT,[0,n]);
maxGradT  = norm(gradT,'inf');
ngradT    = 128+128/maxGradT*gradT;
normGradT = sqrt(sum(gradT.^2,2));
J     = find(normGradT>5e1);
F     = sdiag(Tc-Rc)*gradT;      maxF = norm(F,'inf');
nF    = 128+128/maxF*F;
normF = sqrt(sum(F.^2,2));
K     = find(normF>1e2);

figure(1); clf; set(1,'position',FAIRposition(800),'color','w');
clf; viewImage(ngradT(:,1),omega,m); axis(axT);
clf; viewImage(ngradT(:,2),omega,m); axis(axT);

clf; viewImage(128+0.5*Tc,omega,m); hold on;
qh = quiver(X(J,1),X(J,2),gradT(J,1),gradT(J,2),2);
set(qh,'linewidth',3,'color','k'); axis(axTF);

clf; viewImage(128+0.5*Tc,omega,m); hold on;
qh = quiver(X(K,1),X(K,2),F(K,1),F(K,2),2);
set(qh,'linewidth',3,'color','k'); axis(axTF)

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


text:
% Fast Inter-Harmonic Reconstruction as a pre-process for LPC-based
% spectral envelope estimation.
%
%
% Description
%  This technique reconstructs inter-harmonics of the voice signal,
% as a pre-process for an efficient spectral envelope estimation in
% high-pitched voices. The technique is fully described in [1].

%
% Inputs
%  wave             : [samples] [Nx1] input signal (speech signal)
%  Fs               : [Hz]      [1x1] sampling frequency
%  f0               : F0 estimate (a 0 value is for an unvoiced frame)
%  order            : Order used for the LP analysis
%
%
% Outputs
%  LP               : LP coefficients
%  Energy           : energy of the frame
%
%
% Example
%  Please see the HOWTO_envelope.m example file.
%  Please see http://tcts.fpms.ac.be/~drugman/Toolbox/ for more details.
%
% References
%  [1] T.Drugman, Y. Stylianou, "Fast Inter-Harmonic Reconstruction for
%      Spectral Envelope Estimation in High-Pitched Voices", vol. 21, pp.
%      1418-1422, IEEE Signal Processing Letters, 2014.
%      http://tcts.fpms.ac.be/~drugman/files/SPL-FIHR.pdf
%
% Copyright (c) 2014 Toshiba Cambridge Research Laboratory
%
% License
%  This code will be part of the GLOAT toolbox with the following
%  licence:
%  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.
%
% This function will also be part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Thomas Drugman thomas.drugman@umons.ac.be

function [LP,Energy] = env_fihr(Seg,Fs,F0_local,order)

F0_target=100;
% F0* in the paper. Here set to 100 Hz.

Energy=sum(Seg.^2);

if F0_local>0

    TmpSignal=ones(1,length(Seg));
    t_tmp=1:length(TmpSignal);

    I=ceil(F0_local/F0_target);
    for k=1:I-1
        W=(I-k)/I;
        % We use here a linear weighting function. Other functions are
        % possible, as long as they meet the properties mentioned in
        % the paper.
        TmpSignal=TmpSignal+2*W*cos(2*pi*(k/I)*F0_local/Fs*t_tmp);

        % Equation (2) in the paper
    end

    Seg2=Seg.*TmpSignal';
else
    Seg2=Seg;
end

LP=lpc(Seg2,order);

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep