text:
function groupExperiments_STS(pathWORK,fSetNameType,freedomMat,maxOrder)
% -------------------------------------------------------------------------
% function groupExperiments_STS(pathWORK,fSetNameType,freedomMat,maxOrder)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function finds the best combination of model order and texture
% extraction parameter degree of freedom in terms of [AUC]0.632+ for all
% experiments performed for a given feature set. See ref. [1] for more
% details.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% - pathWORK: Full path to the STS WORKSPACE directory.
% - fSetNameType: String specifying the name of the type of feature set
%                 (e.g., 'PET', 'SEPARATE', 'FUSED', etc.)
% - freedomMat:  Matrix of row vectors of 1's and 0's to specify the degree
%                of freedom on texture extraction parameters for all
%                experiments. For example, for an ith experiment where
%                extraction parameters 1, 2 and 4 in paramAll are allowed
%                to vary, use freedomMat(i,:) = [1,1,0,1].
% - maxOrder: Integer specifying the maximal model order to construct.
% -------------------------------------------------------------------------
% OUTPUTS: The optimal prediction performance results for every model order
%          (in terms of etexture extraction parameter degree of freedom)
%          are saved in a folder named 'RESULTS' in the STS WORKSPACE.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: May 2015
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
cd([pathWORK,'/RESULTS'])

% LOADING ALL RESULTS FOR THE DIFFERENT EXPERIMENTS
nFreedom = size(freedomMat,1);
nParamType = size(freedomMat,2);
cellResults = cell(1,nFreedom);
for i = 1:nFreedom
    nameOpen = ['RESULTS_',fSetNameType,'_'];
    for j = 1:nParamType
        nameOpen = [nameOpen,num2str(freedomMat(i,j))];
    end
    temp = load(nameOpen); temp = struct2cell(temp); temp = temp{1};
    cellResults{i} = temp;
end

% FINDING THE BEST COMBINATIONS (model order/degree of freedom)
results = struct;
for i = 1:maxOrder
    orderName = ['Order',num2str(i)];
    auc = zeros(nFreedom,1);
    for j = 1:nFreedom
        auc(j) = cellResults{j}.(orderName).AUC632;
    end
    [~,indMax] = max(auc);
    results.(orderName) = cellResults{indMax}.(orderName);
    results.(orderName).freedom = freedomMat(indMax,:);
end

% SAVING RESULTS
save(['RESULTS_',fSetNameType,'_BEST'],'results')

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [ a, b ] = p30_lim ( )

%*****************************************************************************80
%
%% P30_LIM returns the integration limits for problem 30.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the lower and upper limits of integration.
%
  a = 2.0;
  b = 7.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [gx,dG_dX,dG_dPhi,d2G_dXdPhi] = g_ERP(x,Phi,u,inG)
% neural mass observation function (DCM for ERPs)

n = size(x,1);          % should be 9
nPhi = length(Phi);     % should be 1

g = exp(Phi(1)).*[0  0   0   0   0   0   0   0   1
                  0  0   0   0   0   0   0   0   0.3
                  0  0   0   0   0   0   0   0   -0.2
                  0  0   0   0   0   0   0   0   -0.7];

% state observation
gx = g*x;

%------ gradients evaluations ------%

% wrt the hidden states
dG_dX = g';

% wrt the observation parameters
dG_dPhi = (g*x)';

% mixed partial derivatives
d2G_dXdPhi(:,:,1) = zeros(n,nPhi);
d2G_dXdPhi(9,1,1) = exp(Phi(1));
d2G_dXdPhi(:,:,2) = 0.3*d2G_dXdPhi(:,:,1);
d2G_dXdPhi(:,:,3) = -0.2*d2G_dXdPhi(:,:,1);
d2G_dXdPhi(:,:,4) = -0.7*d2G_dXdPhi(:,:,1);

author: MBB-team

repo: VBA-toolbox


text:
function [H] = quaternion(q)

% QUATERNION returns the homogenous coordinate transformation matrix corresponding to
% a coordinate transformation described by 7 quaternion parameters.
%
% Use as
%   [H] = quaternion(Q)
% where
%   Q       [q0, q1, q2, q3, q4, q5, q6] vector with parameters
%   H       corresponding homogenous transformation matrix
%
% If the input vector has length 6, it is assumed to represent a unit quaternion without scaling.
%
% See Neuromag/Elekta MaxFilter manual version 2.2, section "D2 Coordinate Matching", page 77 for more details and
% https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Conversion_to_and_from_the_matrix_representation

% Copyright (C) 2016-2017, Robert Oostenveld
%
% 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 2 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% 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 numel(q)==6
  % this is used a lot by the Neuromag/Elekta software, where the first element is left out and a rigid body transformation wothout scaling is used.
  % see also https://github.com/mne-tools/mne-python/blob/maint/0.15/mne/transforms.py#L1137
  q0 = sqrt(1 - q(1)^2 - q(2)^2 - q(3)^2);
  q = [q0 q];
end

if numel(q)~=7
  ft_error('incorrect input vector');
end

% all of these quaternions are zero-offset in the original equation, but one-offset in the MATLAB vector
q0 = q(0+1);
q1 = q(1+1);
q2 = q(2+1);
q3 = q(3+1);
q4 = q(4+1);
q5 = q(5+1);
q6 = q(6+1);

R = [
  q0^2+q1^2-q2^2-q3^2  2*(q1*q2-q0*q3)      2*(q1*q3+q0*q2)
  2*(q1*q2+q0*q3)      q0^2+q2^2-q1^2-q3^2  2*(q2*q3-q0*q1)
  2*(q1*q3-q0*q2)      2*(q2*q3+q0*q1)      q0^2+q3^2-q1^2-q2^2
  ];

T = [q4 q5 q6]';

H = eye(4,4);
H(1:3,1:3) = R;
H(1:3,4)   = T;

author: spm


text:
function [cfg] = ft_topoplotER(cfg, varargin)

% FT_TOPOPLOTER plots the topographic distribution over the head of a 2-dimensional
% data representations such as the event-related potentials or fields, or a power
% or connectivity spectrum.
%
% Use as
%   ft_topoplotER(cfg, timelock)
% or
%   ft_topoplotER(cfg, freq)
%
% The data can be an ERP/ERF produced by FT_TIMELOCKANALYSIS, a power spectrum
% (without time dimension) produced by FT_FREQANALYSIS or a connectivity spectrum
% produced by FT_CONNECTIVITYANALYSIS. Also, the output to FT_FREQSTATISTICS and
% FT_TIMELOCKSTATISTICS can be visualised.
%
% The configuration can have the following parameters
%   cfg.parameter          = field that contains the data to be plotted as color, for example 'avg', 'powspctrm' or 'cohspctrm' (default is automatic)
%   cfg.maskparameter      = field in the data to be used for masking of data. It should have alues between 0 and 1, where 0 corresponds to transparent.
%   cfg.xlim               = limit for 1st dimension in data (e.g., time), can be 'maxmin' or [xmin xmax] (default = 'maxmin')
%   cfg.zlim               = limits for color dimension, 'maxmin', 'maxabs', 'zeromax', 'minzero', or [zmin zmax] (default = 'maxmin')
%   cfg.channel            = Nx1 cell-array with selection of channels (default = 'all'), see FT_CHANNELSELECTION for details
%   cfg.refchannel         = name of reference channel for visualising connectivity, can be 'gui'
%   cfg.baseline           = 'yes','no' or [time1 time2] (default = 'no'), see FT_TIMELOCKBASELINE or FT_FREQBASELINE
%   cfg.baselinetype       = 'absolute' or 'relative' (default = 'absolute')
%   cfg.trials             = 'all' or a selection given as a 1xN vector (default = 'all')
%   cfg.colormap           = string, or Nx3 matrix, see FT_COLORMAP
%   cfg.marker             = 'on', 'labels', 'numbers', 'off'
%   cfg.markersymbol       = channel marker symbol (default = 'o')
%   cfg.markercolor        = channel marker color (default = [0 0 0] (black))
%   cfg.markersize         = channel marker size (default = 2)
%   cfg.markerfontsize     = font size of channel labels (default = 8 pt)
%   cfg.highlight          = 'off', 'on', 'labels', 'numbers'
%   cfg.highlightchannel   =  Nx1 cell-array with selection of channels, or vector containing channel indices see FT_CHANNELSELECTION
%   cfg.highlightsymbol    = highlight marker symbol (default = 'o')
%   cfg.highlightcolor     = highlight marker color (default = [0 0 0] (black))
%   cfg.highlightsize      = highlight marker size (default = 6)
%   cfg.highlightfontsize  = highlight marker size (default = 8)
%   cfg.hotkeys            = enables hotkeys (pageup/pagedown/m) for dynamic zoom and translation (ctrl+) of the color limits
%   cfg.colorbar           = whether to show a colorbar alongside the figure (default = 'no')
%                            'no'                 do not show a colorbar
%                            'yes'                at the default MATLAB location
%                            'North'              inside plot box near top
%                            'South'              inside bottom
%                            'East'               inside right
%                            'West'               inside left
%                            'NorthOutside'       outside plot box near top
%                            'SouthOutside'       outside bottom
%                            'EastOutside'        outside right
%                            'WestOutside'        outside left
%   cfg.colorbartext       =  string indicating the text next to colorbar
%   cfg.interplimits       = limits for interpolation (default = 'head')
%                            'electrodes'         to furthest electrode
%                            'head'               to edge of head
%   cfg.interpolation      = 'linear','cubic','nearest','v4' (default = 'v4') see GRIDDATA
%   cfg.style              = plot style (default = 'both')
%                            'straight'           colormap only
%                            'contour'            contour lines only
%                            'both'               both colormap and contour lines
%                            'fill'               constant color between lines
%                            'blank'               only the head shape
%   cfg.gridscale          = scaling grid size that determines resolution of figure (default = 67)
%   cfg.shading            = 'flat' or 'interp' (default = 'flat')
%   cfg.comment            = 'no', 'auto' or 'xlim' (default = 'auto')
%                            'auto': date, xparam and zparam limits are printed
%                            'xlim': only xparam limits are printed
%   cfg.commentpos         = string or two numbers, position of the comment (default = 'leftbottom')
%                            'lefttop' 'leftbottom' 'middletop' 'middlebottom' 'righttop' 'rightbottom'
%                            'title' to place comment as title
%                            'layout' to place comment as specified for COMNT in layout
%                            [x y] coordinates
%   cfg.interactive        = Interactive plot 'yes' or 'no' (default = 'yes')
%                            In a interactive plot you can select areas and produce a new
%                            interactive plot when a selected area is clicked. Multiple areas
%                            can be selected by holding down the SHIFT key.
%   cfg.directionality     = '', 'inflow' or 'outflow' specifies for
%                            connectivity measures whether the inflow into a
%                            node, or the outflow from a node is plotted. The
%                            (default) behavior of this option depends on the dimord
%                            of the input data (see below).
%   cfg.layout             = specify the channel layout for plotting using one of the supported ways (see below).
%   cfg.interpolatenan     = 'yes' or 'no', whether to interpolate over channels containing NaNs (default = 'yes')
%   cfg.figure             = 'yes' or 'no', whether to open a new figure. You can also specify a figure handle from FIGURE, GCF or SUBPLOT. (default = 'yes')
%
% For the plotting of directional connectivity data the cfg.directionality option
% determines what is plotted. The default value and the supported functionality
% depend on the dimord of the input data. If the input data is of dimord
% 'chan_chan_XXX', the value of directionality determines whether, given the
% reference channel(s), the columns (inflow), or rows (outflow) are selected for
% plotting. In this situation the default is 'inflow'. Note that for undirected
% measures, inflow and outflow should give the same output. If the input data is of
% dimord 'chancmb_XXX', the value of directionality determines whether the rows in
% data.labelcmb are selected. With 'inflow' the rows are selected if the
% refchannel(s) occur in the right column, with 'outflow' the rows are selected if
% the refchannel(s) occur in the left column of the labelcmb-field. Default in this
% case is '', which means that all rows are selected in which the refchannel(s)
% occur. This is to robustly support linearly indexed undirected connectivity
% metrics. In the situation where undirected connectivity measures are linearly
% indexed, specifying 'inflow' or 'outflow' can result in unexpected behavior.
%
% The layout defines how the channels are arranged. You can specify the
% layout in a variety of ways:
%  - you can provide a pre-computed layout structure, see FT_PREPARE_LAYOUT
%  - you can give the name of an ascii layout file with extension *.lay
%  - you can give the name of an electrode file
%  - you can give an electrode definition, i.e. "elec" structure
%  - you can give a gradiometer definition, i.e. "grad" structure
% If you do not specify any of these and the data structure contains an
% electrode or gradiometer structure, that will be used for creating a
% layout. If you want to have more fine-grained control over the layout
% of the subplots, you should create your own layout file.
%
% See also FT_SINGLEPLOTER, FT_MULTIPLOTER, FT_SINGLEPLOTTFR, FT_MULTIPLOTTFR,
% FT_TOPOPLOTTFR, FT_PREPARE_LAYOUT

% Undocumented options:
%
% It is possible to use multiple highlight-selections (e.g.: multiple
% statistical clusters of channels) To do this, all the content of
% the highlight-options (including cfg.highlight) should be placed
% in a cell-array (even if the normal content was already in a
% cell-array). Specific marker settings (e.g. color, size) are defaulted
% when not present.
%
% Example (3 selections):
% cfg.highlight          = {'labels', 'labels', 'numbers'}
% cfg.highlightchannel   = {{'MZF03','MZC01','MRT54'}, [1:5], 'C*'}
% cfg.highlightsymbol    = {'o',[],'+'}        % the empty option will be defaulted
% cfg.highlightcolor     = {'r',[0 0 1]};      % the missing option will be defaulted
% cfg.highlightsize      = [];                 % will be set to default, as will the missing cfg.highlightfontsize
%
% Other options:
% cfg.labeloffset (offset of labels to their marker, default = 0.005)

% Copyright (C) 2005-2020, F.C. Donders Centre
%
% 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$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEVELOPERS NOTE: This code is organized in a similar fashion for multiplot/singleplot/topoplot
% and for ER/TFR and should remain consistent over those 6 functions.
% Section 1: general cfg handling that is independent from the data
% Section 2: data handling, this also includes converting bivariate (chan_chan and chancmb) into univariate data
% Section 3: cfg handling that depends on the data
% Section 4: actual plotting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Section 1: general cfg handling that is independent from the data

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar varargin
ft_preamble provenance varargin

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% this is needed for the figure title
if isfield(cfg, 'dataname') && ~isempty(cfg.dataname)
  dataname = cfg.dataname;
elseif isfield(cfg, 'inputfile') && ~isempty(cfg.inputfile)
  dataname = cfg.inputfile;
elseif nargin>1
  dataname = arrayfun(@inputname, 2:nargin, 'UniformOutput', false);
else
  dataname = {};
end

% make sure figure window titles are labeled appropriately, pass this onto the actual plotting function
cfg.funcname = mfilename;
cfg.dataname = dataname;

% prepare the layout, this should be done only once
tmpcfg = keepfields(cfg, {'layout', 'channel', 'rows', 'columns', 'commentpos', 'skipcomnt', 'scalepos', 'skipscale', 'projection', 'viewpoint', 'rotate', 'width', 'height', 'elec', 'grad', 'opto', 'linecolor', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
cfg.layout = ft_prepare_layout(tmpcfg, varargin{1});

% call the common function that is shared between ft_topoplotER and ft_topoplotTFR
cfg = topoplot_common(cfg, varargin{:});

% remove this field again, it is only used for figure labels
cfg = removefields(cfg, 'funcname');

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous varargin
ft_postamble provenance
ft_postamble savefig

% add a menu to the figure, but only if the current figure does not have subplots
menu_fieldtrip(gcf, cfg, false);

if ~ft_nargout
  % don't return anything
  clear cfg
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function obj_to_tri_surface ( obj_filename, node_filename, element_filename )

%*****************************************************************************80
%
%% OBJ_TO_TRI_SURFACE transfers 3D graphics from an OBJ file to a TRI_SURFACE file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OBJ_FILENAME, the name of the OBJ file to be read.
%
%    Input, string NODE_FILENAME, the name of the TRI_SURFACE node file
%    to be written.
%
%    Input, string ELEMENT_FILENAME, the name of the TRI_SURFACE element
%    file to be written.
%

%
%  Get the size of the objects in the OBJ file.
%
  [ node_num, face_num, normal_num, order_max ] = obj_size ( obj_filename );

  obj_size_print ( obj_filename, node_num, face_num, normal_num, order_max );
%
%  Read the OBJ file.
%
  [ node_xyz, face_order, face_node, normal_vector, vertex_normal ] = ...
    obj_read ( obj_filename, node_num, face_num, normal_num, order_max );
%
%  Convert faces to triangles.
%
  element_num = sum ( face_order(1:face_num) ) - 2 * face_num;

  element_node = faces_to_triangles ( order_max, face_num, face_order, ...
    face_node, element_num );
%
%  Write the data to the TRI_SURFACE files.
%
  dim_num = size ( node_xyz, 1 );
  order_num = 3;

  tri_surface_size_print ( node_filename, element_filename, dim_num, ...
    node_num, order_num, element_num )

  tri_surface_write ( node_filename, element_filename, dim_num, node_num, ...
    order_num, element_num, node_xyz, element_node );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  TRI_SURFACE data written to:\n' );
  fprintf ( 1, '  node file    "%s".\n', node_filename );
  fprintf ( 1, '  element file "%s".\n', element_filename );

  return
end
function c = ch_cap ( c )

%*****************************************************************************80
%
%% CH_CAP capitalizes a single character.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, the character to capitalize.
%
%    Output, character C, the capitalized character.
%
  if ( 'a' <= c & c <= 'z' )
    c = c + 'A' - 'a';
  end

  return
end
function truefalse = ch_eqi ( c1, c2 )

%*****************************************************************************80
%
%% CH_EQI is a case insensitive comparison of two characters for equality.
%
%  Example:
%
%    CH_EQI ( 'A', 'a' ) is TRUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C1, C2, the characters to compare.
%
%    Output, logical TRUEFALSE, is TRUE (1) if the characters are equal.
%
  FALSE = 0;
  TRUE = 1;

  if ( ch_cap ( c1 ) == ch_cap ( c2 ) )
    truefalse = TRUE;
  else
    truefalse = FALSE;
  end

  return
end
function value = ch_index ( s, c )

%*****************************************************************************80
%
%% CH_INDEX is the first occurrence of a character in a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be searched.
%
%    Input, character C, the character to be searched for.
%
%    Output, integer VALUE, the location of the first occurrence of C
%    in the string, or 0 if C does not occur.
%
  value = 0;

  for i = 1 : length ( s )

    if ( s(i:i) == c )
      value = i;
      return
    end

  end

  return
end
function value = ch_is_control ( ch )

%*****************************************************************************80
%
%% CH_IS_CONTROL is TRUE if a character is a control character.
%
%  Discussion:
%
%    A "control character" has ASCII code <= 31 or 127 <= ASCII code.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character CH, the character to be tested.
%
%    Output, integer CH_IS_CONTROL, TRUE if the character is a control
%    character, and FALSE otherwise.
%
  if ( ch <= 31 | 127 <= ch )
    value = 1;
  else
    value = 0;
  end

  return
end
function truefalse = ch_is_digit ( c )

%*****************************************************************************80
%
% CH_IS_DIGIT returns TRUE if the character C is a digit.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 December 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, a character.
%
%    Output, integer TRUEFALSE, is TRUE (1) if C is a digit, FALSE (0) otherwise.
%
  TRUE = 1;
  FALSE = 0;

  if ( '0' <= c & c <= '9' )
    truefalse = TRUE;
  else
    truefalse = FALSE;
  end

  return
end
function digit = ch_to_digit ( c )

%*****************************************************************************80
%
%% CH_TO_DIGIT returns the integer value of a base 10 digit.
%
%  Example:
%
%     C   DIGIT
%    ---  -----
%    '0'    0
%    '1'    1
%    ...  ...
%    '9'    9
%    ' '    0
%    'X'   -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, the decimal digit, '0' through '9' or blank
%    are legal.
%
%    Output, integer DIGIT, the corresponding integer value.  If C was
%    'illegal', then DIGIT is -1.
%
  if ( '0' <= c & c <= '9' )

    digit = c - '0';

  elseif ( c == ' ' )

    digit = 0;

  else

    digit = -1;

  end

  return
end
function triangle_node = faces_to_triangles ( order_max, face_num, ...
  face_order, face_node, triangle_num )

%*****************************************************************************80
%
%% FACES_TO_TRIANGLES replaces higher order faces by triangles.
%
%  Modified:
%
%    28 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ORDER_MAX, the maximum face order.
%
%    Input, integer FACE_NUM, the number of faces.
%
%    Input, integer FACE_ORDER(FACE_NUM), the order of each face.
%
%    Input, integer FACE_NODE(ORDER_MAX,FACE_NUM), the nodes that make up each face.
%
%    Input, integer TRIANGLE_NUM, the number of triangles that will be generated.
%
%    Output, integer TRIANGLE_NODE(3,TRIANGLE_NUM), the triangles.
%
  triangle_node = zeros ( 3, triangle_num );

  triangle = 0;

  for face = 1 : face_num
    for order = 3 : face_order ( face )
      triangle = triangle + 1;
      triangle_node(1,triangle) = face_node(1,face);
      triangle_node(2,triangle) = face_node(order-1,face);
      triangle_node(3,triangle) = face_node(order,face);
    end
  end

  return
end
function i4mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Discussion:
%
%    An I4MAT is an array of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 June 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %12d', round ( table(i,j) ) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function [ node_xyz, face_order, face_node, normal_vector, vertex_normal ] = ...
  obj_read ( input_file_name, node_num, face_num, normal_num, order_max )

%*****************************************************************************80
%
%% OBJ_READ reads graphics information from a Wavefront OBJ file.
%
%  Discussion:
%
%    It is intended that the information read from the file can
%    either start a whole new graphics object, or simply be added
%    to a current graphics object via the '<<' command.
%
%    This is controlled by whether the input values have been zeroed
%    out or not.  This routine simply tacks on the information it
%    finds to the current graphics object.
%
%  Example:
%
%    #  magnolia.obj
%
%    v -3.269770 -39.572201 0.876128
%    v -3.263720 -39.507999 2.160890
%    ...
%    v 0.000000 -9.988540 0.000000
%    vn 1.0 0.0 0.0
%    ...
%    vn 0.0 1.0 0.0
%
%    f 8 9 11 10
%    f 12 13 15 14
%    ...
%    f 788 806 774
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer NODE_NUM, the number of points.
%
%    Input, integer FACE_NUM, the number of faces.
%
%    Input, integer NORMAL_NUM, the number of normal vectors.
%
%    Input, integer ORDER_MAX, the maximum number of vertices per face.
%
%    Output, real ( kind = 8 ) NODE_XYZ(3,NODE_NUM), the coordinates of points.
%
%    Output, integer FACE_ORDER(FACE_NUM), the number of vertices per face.
%
%    Output, integer FACE_NODE(ORDER_MAX,FACE_NUM), the nodes making faces.
%
%    Output, real ( kind = 8 ) NORMAL_VECTOR(3,NORMAL_NUM), normal vectors.
%
%    Output, integer VERTEX_NORMAL(ORDER_MAX,FACE_NUM), the indices of normal
%    vectors per vertex.
%
  ierror = 0;

  face = 0;
  node = 0;
  normal = 0;
  text_num = 0;

  face_node = zeros ( order_max, face_num );
  face_order = zeros ( face_num, 1 );
  node_xyz = zeros ( 3, node_num );
  normal_vector = zeros ( 3, normal_num );
  vertex_normal = zeros ( order_max, face_num );
%
%  If no file input, try to get one from the user.
%
  if ( nargin < 1 )
    input_file_name = input ( 'Enter the name of the ASCII OBJ file.' );
    if ( isempty ( input_file_name ) )
      return
    end
  end
%
%  Open the file.
%
  input_file_unit = fopen ( input_file_name, 'r' );

  if ( input_file_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'OBJ_READ - Fatal error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'OBJ_READ - Fatal error!' );
    return
  end
%
%  Read a line of text from the file.
%
  while ( 1 )

    text = fgetl ( input_file_unit );

    if ( text == -1 )
      break
    end

    text_num = text_num + 1;
%
%  Replace any control characters (in particular, TAB's) by blanks.
%
    s_control_blank ( text );

    done = 1;
    word_index = 0;
%
%  Read a word from the line.
%
    [ word, done ] = word_next_read ( text, done );
%
%  If no more words in this line, read a new line.
%
    if ( done )
      continue
    end
%
%  If this word begins with '#' or '$', then it's a comment.  Read a new line.
%
    if ( word(1) == '#' | word(1) == '$' )
      continue
    end

    word_index = word_index + 1;

    if ( word_index == 1 )
      word_one = word;
    end
%
%  BEVEL
%  Bevel interpolation.
%
    if ( s_eqi ( word_one, 'BEVEL' ) )
%
%  BMAT
%  Basis matrix.
%
    elseif ( s_eqi ( word_one, 'BMAT' ) )
%
%  C_INTERP
%  Color interpolation.
%
    elseif ( s_eqi ( word_one, 'C_INTERP' ) )
%
%  CON
%  Connectivity between free form surfaces.
%
    elseif ( s_eqi ( word_one, 'CON' ) )
%
%  CSTYPE
%  Curve or surface type.
%
    elseif ( s_eqi ( word_one, 'CSTYPE' ) )
%
%  CTECH
%  Curve approximation technique.
%
    elseif ( s_eqi ( word_one, 'CTECH' ) )
%
%  CURV
%  Curve.
%
    elseif ( s_eqi ( word_one, 'CURV' ) )
%
%  CURV2
%  2D curve.
%
    elseif ( s_eqi ( word_one, 'CURV2' ) )
%
%  D_INTERP
%  Dissolve interpolation.
%
    elseif ( s_eqi ( word_one, 'D_INTERP' ) )
%
%  DEG
%  Degree.
%
    elseif ( s_eqi ( word_one, 'DEG' ) )
%
%  END
%  End statement.
%
    elseif ( s_eqi ( word_one, 'END' ) )
%
%  F V1 V2 V3 ...
%    or
%  F V1/VT1/VN1 V2/VT2/VN2 ...
%    or
%  F V1//VN1 V2//VN2 ...
%
%  Face.
%  A face is defined by the vertices.
%  Optionally, slashes may be used to include the texture vertex
%  and vertex normal indices.
%
    elseif ( s_eqi ( word_one, 'F' ) )

      face = face + 1;

      vertex = 0;

      while ( 1 )

        [ word, done ] = word_next_read ( text, done );

        if ( done )
          break
        end

        vertex = vertex + 1;
        order_max = max ( order_max, vertex );
%
%  Locate the slash characters in the word, if any.
%
        i1 = ch_index ( word, '/' );
        if ( 0 < i1 )
          i2 = ch_index ( word(i1+1), '/' ) + i1;
        else
          i2 = 0;
        end
%
%  Read the vertex index.
%
        itemp = s_to_i4 ( word );

        face_node(vertex,face) = itemp;
        face_order(face) = face_order(face) + 1;
%
%  If there are two slashes, then read the data following the second one.
%
        if ( 0 < i2 )

          itemp = s_to_i4 ( word(i2+1) );

          vertex_normal(vertex,face) = itemp;

        end

      end
%
%  G
%  Group name.
%
    elseif ( s_eqi ( word_one, 'G' ) )
%
%  HOLE
%  Inner trimming loop.
%
    elseif ( s_eqi ( word_one, 'HOLE' ) )
%
%  L
%  A line, described by a sequence of vertex indices.
%  Are the vertex indices 0 based or 1 based?
%
    elseif ( s_eqi ( word_one, 'L' ) )
%
%  LOD
%  Level of detail.
%
    elseif ( s_eqi ( word_one, 'LOD' ) )
%
%  MG
%  Merging group.
%
    elseif ( s_eqi ( word_one, 'MG' ) )
%
%  MTLLIB
%  Material library.
%
    elseif ( s_eqi ( word_one, 'MTLLIB' ) )
%
%  O
%  Object name.
%
    elseif ( s_eqi ( word_one, 'O' ) )
%
%  P
%  Point.
%
    elseif ( s_eqi ( word_one, 'P' ) )
%
%  PARM
%  Parameter values.
%
    elseif ( s_eqi ( word_one, 'PARM' ) )
%
%  S
%  Smoothing group.
%
    elseif ( s_eqi ( word_one, 'S' ) )
%
%  SCRV
%  Special curve.
%
    elseif ( s_eqi ( word_one, 'SCRV' ) )
%
%  SHADOW_OBJ
%  Shadow casting.
%
    elseif ( s_eqi ( word_one, 'SHADOW_OBJ' ) )
%
%  SP
%  Special point.
%
    elseif ( s_eqi ( word_one, 'SP' ) )
%
%  STECH
%  Surface approximation technique.
%
    elseif ( s_eqi ( word_one, 'STECH' ) )
%
%  STEP
%  Stepsize.
%
    elseif ( s_eqi ( word_one, 'STEP' ) )
%
%  SURF
%  Surface.
%
    elseif ( s_eqi ( word_one, 'SURF' ) )
%
%  TRACE_OBJ
%  Ray tracing.
%
    elseif ( s_eqi ( word_one, 'TRACE_OBJ' ) )
%
%  TRIM
%  Outer trimming loop.
%
    elseif ( s_eqi ( word_one, 'TRIM' ) )
%
%  USEMTL
%  Material name.
%
    elseif ( s_eqi ( word_one, 'USEMTL' ) )
%
%  V X Y Z
%  Geometric vertex.
%
    elseif ( s_eqi ( word_one, 'V' ) )

      node = node + 1;

      for i = 1 : 3
        [ word, done ] = word_next_read ( text, done );
        temp = s_to_r8 ( word );
        node_xyz(i,node) = temp;
      end
%
%  VN
%  Vertex normals.
%
    elseif ( s_eqi ( word_one, 'VN' ) )

      normal = normal + 1;

      for i = 1 : 3
        [ word, done ] = word_next_read ( text, done );
        temp = s_to_r8 ( word );
        normal_vector(i,normal) = temp;
      end
%
%  VT
%  Vertex texture.
%
    elseif ( s_eqi ( word_one, 'VT' ) )
%
%  VP
%  Parameter space vertices.
%
    elseif ( s_eqi ( word_one, 'VP' ) )
%
%  Unrecognized keyword.
%
    else

    end

  end

  fclose ( input_file_unit );

  return
end
function [ node_num, face_num, normal_num, order_max ] = obj_size ( ...
  input_file_name )

%*****************************************************************************80
%
%% OBJ_SIZE determines sizes of graphics objects in an Alias OBJ file.
%
%  Discussion:
%
%    The only items of interest to this routine are vertices,
%    faces, and normal vectors.
%
%  Example:
%
%    #  magnolia.obj
%
%    v -3.269770 -39.572201 0.876128
%    v -3.263720 -39.507999 2.160890
%    ...
%    v 0.000000 -9.988540 0.000000
%
%    vn 1.0 0.0 0.0
%    ...
%    vn 0.0 1.0 0.0
%
%    f 8 9 11 10
%    f 12 13 15 14
%    ...
%    f 788 806 774
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the input file name.
%
%    Output, integer NODE_NUM, the number of points.
%
%    Output, integer FACE_NUM, the number of faces.
%
%    Output, integer NORMAL_NUM, the number of normal vectors.
%
%    Output, integer ORDER_MAX, the maximum face order.
%
  ierror = 0;

  face_num = 0;
  node_num = 0;
  normal_num = 0;
  order_max = 0;
  text_num = 0;
%
%  If no file input, try to get one from the user.
%
  if ( nargin < 1 )
    input_file_name = input ( 'Enter the name of the ASCII OBJ file.' );
    if ( isempty ( input_file_name ) )
      return
    end
  end
%
%  Open the file.
%
  input_file_unit = fopen ( input_file_name, 'r' );

  if ( input_file_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'OBJ_SIZE - Fatal error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'OBJ_SIZE - Fatal error!' );
    return
  end
%
%  Read a line of text from the file.
%
  while ( 1 )

    text = fgetl ( input_file_unit );

    if ( text == -1 )
      break
    end

    text_num = text_num + 1;
%
%  Replace any control characters (in particular, TABs) by blanks.
%
    s_control_blank ( text );

    done = 1;
    word_index = 0;
%
%  Read a word from the line.
%
    [ word, done ] = word_next_read ( text, done );
%
%  If no more words in this line, read a new line.
%
    if ( done )
      continue
    end
%
%  If this word begins with '#' or '$', then it is a comment.  Read a new line.
%
    if ( word(1) == '#' | word(1) == '$' )
      continue
    end

    word_index = word_index + 1;

    if ( word_index == 1 )
      word_one = word;
    end
%
%  F V1 V2 V3 ...
%    or
%  F V1/VT1/VN1 V2/VT2/VN2 ...
%    or
%  F V1//VN1 V2//VN2 ...
%
%  Face.
%  A face is defined by the vertices.
%  Optionally, slashes may be used to include the texture vertex
%  and vertex normal indices.
%
    if ( s_eqi ( word_one, 'F' ) )

      face_num = face_num + 1;

      vertex = 0;

      while ( 1 )

        [ word, done ] = word_next_read ( text, done );

        if ( done )
          break
        end

        vertex = vertex + 1;
        order_max = max ( order_max, vertex );
%
%  Locate the slash characters in the word, if any.
%
        i1 = ch_index ( word, '/' );
        if ( 0 < i1 )
          i2 = ch_index ( word(i1+1), '/' ) + i1;
        else
          i2 = 0;
        end
%
%  Read the vertex index.
%
        itemp = s_to_i4 ( word );
%
%  If there are two slashes, then read the data following the second one.
%
        if ( 0 < i2 )
          itemp = s_to_i4 ( word(i2+1) );
        end

      end
%
%  V X Y Z W
%  Geometric vertex.
%
    elseif ( s_eqi ( word_one, 'V' ) )

      node_num = node_num + 1;
      continue
%
%  VN
%  Vertex normals.
%
    elseif ( s_eqi ( word_one, 'VN' ) )

      normal_num = normal_num + 1;
      continue

    end

  end

  fclose ( input_file_unit );

  return
end
function obj_size_print ( input_file_name, node_num, face_num, normal_num, ...
  order_max )

%*****************************************************************************80
%
%% OBJ_SIZE_PRINT prints sizes associated with an OBJ file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer NODE_NUM, the number of vertices defined.
%
%    Input, integer FACE_NUM, the number of faces defined.
%
%    Input, integer NORMAL_NUM, the number of normal vectors defined.
%
%    Input, integer ORDER_MAX, the maximum number of vertices per face.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Object sizes for OBJ file "%s":\n', input_file_name );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Nodes              = %d\n', node_num );
  fprintf ( 1, '  Faces              = %d\n', face_num );
  fprintf ( 1, '  Maximum face order = %d\n', order_max );
  fprintf ( 1, '  Normal vectors     = %d\n', normal_num );

  return
end
function value = r8_epsilon ( )

%*****************************************************************************80
%
%% R8_EPSILON returns the R8 roundoff unit.
%
%  Discussion:
%
%    The roundoff unit is a number R which is a power of 2 with the
%    property that, to the precision of the computer's arithmetic,
%      1 < 1 + R
%    but
%      1 = ( 1 + R / 2 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VALUE, the roundoff unit.
%
  value = eps;

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Discussion:
%
%    An R8MAT is an array of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
%  For smaller data files, and less precision, try:
%
%     fprintf ( output_unit, '  %14.6e', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %24.16e', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function s = s_control_blank ( s )

%*****************************************************************************80
%
%% S_CONTROL_BLANK replaces control characters with blanks.
%
%  Discussion:
%
%    A "control character" has ASCII code <= 31 or 127 <= ASCII code.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, string S, the string to be transformed.
%
  s_length = s_len_trim ( s );

  for i = 1 : s_length
    if ( ch_is_control ( s(i) ) )
      s(i) = ' ';
    end
  end

  return
end
function value = s_eqi ( s1, s2 )

%*****************************************************************************80
%
%% S_EQI is a case insensitive comparison of two strings for equality.
%
%  Example:
%
%    S_EQI ( 'Anjana', 'ANJANA' ) is TRUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S1, S2, the strings to compare.
%
%    Output, logical VALUE, is TRUE if the strings are equal.
%
  FALSE = 0;
  TRUE = 1;

  len1 = length ( s1 );
  len2 = length ( s2 );
  lenc = min ( len1, len2 );

  value = FALSE;

  for i = 1 : lenc

    c1 = ch_cap ( s1(i) );
    c2 = ch_cap ( s2(i) );

    if ( c1 ~= c2 )
      value = FALSE;
      return
    end

  end

  for i = lenc + 1 : len1
    if ( s1(i) ~= ' ' )
      value = FALSE;
      return
    end
  end

  for i = lenc + 1 : len2
    if ( s2(i) ~= ' ' )
      value = FALSE;
      return
    end
  end

  value = TRUE;

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
%% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LEN, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
function ival = s_to_i4 ( s )

%*****************************************************************************80
%
%% S_TO_I4 reads an integer value from a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, a string to be examined.
%
%    Output, integer IVAL, the integer value read from the string.
%
  sgn = 1;
  state = 0;
  ival = 0;

  i = 0;

  while ( i < s_len_trim ( s ) )

    i = i + 1;
    c = s(i);

    if ( state == 0 )

      if ( c == ' ' )

      elseif ( c == '-' )
        state = 1;
        sgn = -1;
      elseif ( c == '+' )
        state = 1;
        sgn = +1;
      elseif ( '0' <= c & c <= '9' )
        state = 2;
        ival = c - '0';
      else
        fprintf ( '\n' );
        fprintf ( 'S_TO_I4 - Fatal error!\n' );
        fprintf ( '  Illegal character %c while in state %d.\n', c, state );
        return;
      end
%
%  Have read the sign, now expecting the first digit.
%
    elseif ( state == 1 )

      if ( c == ' ' )

      elseif ( '0' <= c & c <= '9' )
        state = 2;
        ival = c - '0';
      else
        fprintf ( '\n' );
        fprintf ( 'S_TO_I4 - Fatal error!\n' );
        fprintf ( '  Illegal character %c while in state %d.\n', c, state );
        return
      end
%
%  Have read at least one digit, expecting more.
%
    elseif ( state == 2 )

      if ( '0' <= c & c <= '9' )
        ival = 10 * ival + c - '0';
      else
        ival = sgn * ival;
        return;
      end

    end

  end
%
%  If we read all the characters in the string, see if we're OK.
%
  if ( state ~= 2 )
    fprintf ( '\n' );
    fprintf ( 'S_TO_I4 - Fatal error!\n' );
    fprintf ( '  Did not read enough information to define an integer!\n' );
    return;
  end

  ival = sgn * ival;

  return
end
function [ r, lchar, ierror ] = s_to_r8 ( s )

%*****************************************************************************80
%
%% S_TO_R8 reads an R8 from a string.
%
%  Discussion:
%
%    This routine will read as many characters as possible until it reaches
%    the end of the string, or encounters a character which cannot be
%    part of the real number.
%
%    Legal input is:
%
%       1 blanks,
%       2 '+' or '-' sign,
%       2.5 spaces
%       3 integer part,
%       4 decimal point,
%       5 fraction part,
%       6 'E' or 'e' or 'D' or 'd', exponent marker,
%       7 exponent sign,
%       8 exponent integer part,
%       9 exponent decimal point,
%      10 exponent fraction part,
%      11 blanks,
%      12 final comma or semicolon.
%
%    with most quantities optional.
%
%  Example:
%
%    S                 R
%
%    '1'               1.0
%    '     1   '       1.0
%    '1A'              1.0
%    '12,34,56'        12.0
%    '  34 7'          34.0
%    '-1E2ABCD'        -100.0
%    '-1X2ABCD'        -1.0
%    ' 2E-1'           0.2
%    '23.45'           23.45
%    '-4.2E+2'         -420.0
%    '17d2'            1700.0
%    '-14e-2'         -0.14
%    'e2'              100.0
%    '-12.73e-9.23'   -12.73 * 10.0**(-9.23)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string containing the
%    data to be read.  Reading will begin at position 1 and
%    terminate at the end of the string, or when no more
%    characters can be read to form a legal real.  Blanks,
%    commas, or other nonnumeric data will, in particular,
%    cause the conversion to halt.
%
%    Output, real R, the value that was read from the string.
%
%    Output, integer LCHAR, the number of characters of S that were used to form R.
%
%    Output, integer IERROR, is 0 if no error occurred.
%
  s_length = s_len_trim ( s );
  ierror = 0;
  r = 0.0;
  lchar = -1;
  isgn = 1;
  rtop = 0.0;
  rbot = 1.0;
  jsgn = 1;
  jtop = 0;
  jbot = 1;
  ihave = 1;
  iterm = 0;

  while ( 1 )

    lchar = lchar + 1;
    c = s(lchar+1);
%
%  Blank character.
%
    if ( c == ' ' )

      if ( ihave == 2 )

      elseif ( ihave == 6 | ihave == 7 )
        iterm = 1;
      elseif ( 1 < ihave )
        ihave = 11;
      end
%
%  Comma.
%
    elseif ( c == ',' | c == ';' )

      if ( ihave ~= 1 )
        iterm = 1;
        ihave = 12;
        lchar = lchar + 1;
      end
%
%  Minus sign.
%
    elseif ( c == '-' )

      if ( ihave == 1 );
        ihave = 2;
        isgn = -1;
      elseif ( ihave == 6 )
        ihave = 7;
        jsgn = -1;
      else
        iterm = 1;
      end
%
%  Plus sign.
%
    elseif ( c == '+' )

      if ( ihave == 1 )
        ihave = 2;
      elseif ( ihave == 6 )
        ihave = 7;
      else
        iterm = 1;
      end
%
%  Decimal point.
%
    elseif ( c == '.' )

      if ( ihave < 4 )
        ihave = 4;
      elseif ( 6 <= ihave & ihave <= 8 )
        ihave = 9;
      else
        iterm = 1;
      end
%
%  Exponent marker.
%
    elseif ( ch_eqi ( c, 'E' ) | ch_eqi ( c, 'D' ) )

      if ( ihave < 6 )
        ihave = 6;
      else
        iterm = 1;
      end
%
%  Digit.
%
    elseif ( ihave < 11 & ch_is_digit ( c ) )

      if ( ihave <= 2 )
        ihave = 3;
      elseif ( ihave == 4 )
        ihave = 5;
      elseif ( ihave == 6 | ihave == 7 )
        ihave = 8;
      elseif ( ihave == 9 )
        ihave = 10;
      end

      d = ch_to_digit ( c );

      if ( ihave == 3 )
        rtop = 10.0 * rtop + d;
      elseif ( ihave == 5 )
        rtop = 10.0 * rtop + d;
        rbot = 10.0 * rbot;
      elseif ( ihave == 8 )
        jtop = 10 * jtop + d;
      elseif ( ihave == 10 )
        jtop = 10 * jtop + d;
        jbot = 10 * jbot;
      end
%
%  Anything else is regarded as a terminator.
%
    else
      iterm = 1;
    end
%
%  If we haven't seen a terminator, and we haven't examined the
%  entire string, go get the next character.
%
    if ( iterm == 1 | s_length <= lchar + 1 )
      break;
    end

  end
%
%  If we haven't seen a terminator, and we have examined the
%  entire string, then we're done, and LCHAR is equal to S_LENGTH.
%
  if ( iterm ~= 1 & lchar + 1 == s_length )
    lchar = s_length;
  end
%
%  Number seems to have terminated.  Have we got a legal number?
%  Not if we terminated in states 1, 2, 6 or 7!
%
  if ( ihave == 1 | ihave == 2 | ihave == 6 | ihave == 7 )
    ierror = ihave;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'S_TO_R8 - Fatal error!\n' );
    fprintf ( 1, '  IHAVE = %d\n', ihave );
    error ( 'S_TO_R8 - Fatal error!' );
  end
%
%  Number seems OK.  Form it.
%
  if ( jtop == 0 )
    rexp = 1.0;
  else

    if ( jbot == 1 )
      rexp = 10.0^( jsgn * jtop );
    else
      rexp = jsgn * jtop;
      rexp = rexp / jbot;
      rexp = 10.0^rexp;
    end

  end

  r = isgn * rexp * rtop / rbot;

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
function tri_surface_size_print ( node_filename, element_filename, dim_num, ...
  node_num, order_num, triangle_num )

%*****************************************************************************80
%
%% TRI_SURFACE_SIZE_PRINT prints sizes associated with a TRI_SURFACE file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string NODE_FILENAME, the name of the node file.
%
%    Input, string TRIANGLE_FILENAME, the name of the triangle file.
%
%    Input, integer DIM_NUM, the number of spatial dimensions.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer ORDER_NUM, the order of the triangles.
%
%    Input, integer TRIANGLE_NUM, the number of triangles.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRI_SURFACE_SIZE_PRINT:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Node file     "%s".\n', node_filename );
  fprintf ( 1, '  Triangle file "%s".\n', element_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension  = %d\n', dim_num );
  fprintf ( 1, '  Nodes              = %d\n', node_num );
  fprintf ( 1, '  Triangle order     = %d\n', order_num );
  fprintf ( 1, '  Triangles          = %d\n', triangle_num );

  return
end
function tri_surface_write ( node_filename, element_filename, dim_num, ...
  node_num, order_num, triangle_num, node_xyz, triangle_node )

%*****************************************************************************80
%
%% TRI_SURFACE_WRITE writes graphics information to a pair of TRI_SURFACE files.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string NODE_FILENAME, the name of the node file.
%
%    Input, string ELEMENT_FILENAME, the name of the element file.
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer NODE_NUM, the number of points.
%
%    Input, integer ORDER_NUM, the order of the triangles.
%
%    Input, integer TRIANGLE_NUM, the number of triangles.
%
%    Input, real NODE_XYZ(DIM_NUM,NODE_NUM), the node coordinates.
%
%    Input, integer TRIANGLE_NODE(ORDER_NUM,TRIANGLE_NUM),
%    the nodes that form the triangles.
%
  r8mat_write ( node_filename, dim_num, node_num, node_xyz );

  i4mat_write ( element_filename, order_num, triangle_num, ...
    triangle_node );

  return
end
function [ word, done ] = word_next_read ( s, done )

%*****************************************************************************80
%
%% WORD_NEXT_READ "reads" words from a string, one at a time.
%
%  Special cases:
%
%    The following characters are considered to be a single word,
%    whether surrounded by spaces or not:
%
%      " ( ) { } [ ]
%
%    Also, if there is a trailing comma on the word, it is stripped off.
%    This is to facilitate the reading of lists.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 December 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, a string, presumably containing words
%    separated by spaces.
%
%    Input, logical DONE.
%    TRUE, if we are beginning a new string;
%    FALSE, if we are continuing to process the current string.
%
%    Output, string WORD.
%    If DONE is FALSE, then WORD contains the "next" word read.
%    If DONE is TRUE, then WORD is blank, because there was no more to read.
%
%    Output, logical DONE.
%      FALSE if another word was read,
%      TRUE if no more words could be read.
%
  persistent lenc;
  persistent next;

  tab = char ( 9 );
%
%  We "remember" LENC and NEXT from the previous call.
%
%  An input value of DONE = TRUE signals a new line of text to examine.
%
  if ( done )

    next = 1;
    done = 0;
    lenc = s_len_trim ( s );

    if ( lenc <= 0 )
      done = 1;
      word = ' ';
      return
    end

  end
%
%  Beginning at index NEXT, search the string for the next nonblank,
%  which signals the beginning of a word.
%
  ilo = next;
%
%  ...S(NEXT:) is blank.  Return with WORD = ' ' and DONE = TRUE.
%
  while ( 1 )

    if ( lenc < ilo )
      word = ' ';
      done = 1;
      next = lenc + 1;
      return
    end
%
%  If the current character is blank, skip to the next one.
%
    if ( s(ilo) ~= ' ' & s(ilo) ~= tab )
      break
    end

    ilo = ilo + 1;

  end
%
%  ILO is the index of the next nonblank character in the string.
%
%  If this initial nonblank is a special character,
%  then that's the whole word as far as we're concerned,
%  so return immediately.
%
  if ( s(ilo) == '"' | ...
       s(ilo) == '(' | ...
       s(ilo) == ')' | ...
       s(ilo) == '{' | ...
       s(ilo) == '}' | ...
       s(ilo) == '[' | ...
       s(ilo) == ']' )

    word = s(ilo);
    next = ilo + 1;
    return

  end
%
%  Now search for the last contiguous character that is not a
%  blank, TAB, or special character.
%
  next = ilo + 1;

  while ( next <= lenc )

    if ( s(next) == ' ' )
      break;
    elseif ( s(next) == tab )
      break;
    elseif ( s(next) == '"' )
      break;
    elseif ( s(next) == '(' )
      break;
    elseif ( s(next) == ')' )
      break;
    elseif ( s(next) == '{' )
      break;
    elseif ( s(next) == '}' )
      break;
    elseif ( s(next) == '[' )
      break;
    elseif ( s(next) == ']' )
      break;
    end

    next = next + 1;

  end

  if ( s(next-1) == ',' )
    word = s(ilo:next-2);
  else
    word = s(ilo:next-1);
  end

  return
end


text:
function er4OC
%EG4OC    Example 4 of optimal control tutorial.
%    This example is from the following reference:
%    S.N.Avvakumov and Yu.N.Kiselev, "Boundary value problem for ordinary
%    differential equations with applications to optimal control".
%    Example 5 on page 8

global mu;
mu = 0.5;
solinit = bvpinit(linspace(0,1,10),[2;2;0;-1],4);

sol = bvp4c(@ode,@bc,solinit);

figure(1)
lines = {'k-.','r--','b-'};
ut = sol.y(4,:)./sqrt(mu*sol.y(3,:).^2+sol.y(4,:).^2);
t = sol.parameters*sol.x;
plot(t,ut,lines{1});
% axis([-0.1 1.1 0 1.7]);
title('Smoothed control');
xlabel('t');
ylabel('u');
drawnow; hold on;

figure(2)
plot(sol.y(1,:),sol.y(2,:),lines{1});
% axis([-0.1 1.1 0 1.7]);
title('Admissible trajectories');
xlabel('x_1');
ylabel('x_2');
drawnow; hold on;

% the solution for one value of mu is used as guess for the next.
for i=2:3
  if i==2
      mu = 0.3;
  else
      mu = .1;
  end;
  % After creating function handles, the new value of mu
  % will be used in nested functions.
  sol = bvp4c(@ode,@bc,sol);
  figure(2)
  plot(sol.y(1,:),sol.y(2,:),lines{i});
  drawnow;

  figure(1)
  ut = sol.y(4,:)./sqrt(mu*sol.y(3,:).^2+sol.y(4,:).^2);
  t = sol.parameters*sol.x;
  plot(t,ut,lines{i});
  drawnow;

end
figure(1); legend('mu = 0.5', 'mu = 0.3', 'mu = 0.1', ...
                  'Location', 'NorthWest'); hold off;
figure(2); legend('mu = 0.5', 'mu = 0.3', 'mu = 0.1', ...
                  'Location', 'NorthWest'); hold off;

% print -djpeg90 -f1 -r300 eg4_trajectory.jpg;
% print -djpeg90 -f2 -r300 eg4_control.jpg

% --------------------------------------------------------------------------
function dydt = ode(t,y,T)
global mu;
term = sqrt(mu*y(3)^2+y(4)^2);
dydt = T*[ y(2) + mu*y(3)/term
           y(4)/term
           0
           -y(3)];

% --------------------------------------------------------------------------
% boundary conditions, with 4 states and 1 parameters, 5 conditions are
% needed: x1(0) =2, x2(0) = 2; x1(1) = 0; x2(1) = 0; x3(1)^2+x4(1)^2 = 1;
function res = bc(ya,yb,T)
res = [ya(1)-2
       ya(2)-2
       yb(1)
       yb(2)
       yb(3)^2+yb(4)^2-1];

author: Sable

repo: mcbench-benchmarks


text:
function varargout = drawLine(lin, varargin)
%DRAWLINE Draw a straight line clipped by the current axis
%
%   drawLine(LINE);
%   Draws the line LINE on the current axis, by using current axis to clip
%   the line.
%
%   drawLine(LINE, PARAM, VALUE);
%   Specifies drawing options.
%
%   H = drawLine(...)
%   Returns a handle to the created line object. If clipped line is not
%   contained in the axis, the function returns -1.
%
%   Example
%     figure; hold on; axis equal;
%     axis([0 100 0 100]);
%     drawLine([30 40 10 20]);
%     drawLine([30 40 20 -10], 'Color', 'm', 'LineWidth', 2);
%     drawLine([-30 140 10 20]);
%
%   See also:
%   lines2d, createLine, drawEdge
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 31/10/2003.
%

%   HISTORY
%   25/05/2004 add support for multiple lines (loop)
%   23/05/2005 add support for arguments
%   03/08/2010 bug for lines outside box (thanks to Reto Zingg)
%   04/08/2010 rewrite using clipLine
%   2011-10-11 add management of axes handle

% extract handle of axis to draw in
if isAxisHandle(lin)
    ax = lin;
    lin = varargin{1};
    varargin(1) = [];
else
    ax = gca;
end

% default style for drawing lines
if length(varargin) ~= 1
    varargin = [{'color', 'b'}, varargin];
end

% extract bounding box of the current axis
xlim = get(ax, 'xlim');
ylim = get(ax, 'ylim');

% clip lines with current axis box
clip = clipLine(lin, [xlim ylim]);
ok   = isfinite(clip(:,1));

% initialize result array to invalide handles
h = -1 * ones(size(lin, 1), 1);

% draw valid lines
h(ok) = plot(ax, clip(ok, [1 3])', clip(ok, [2 4])', varargin{:});

% return line handle if needed
if nargout > 0
    varargout = {h};
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
% Author : FUAT COGUN
% Date   : 16.07.2009
%
%
% Find the Covariance Matrix of given region.
%
%
% Inputs  : -I                                                  (uint8 Matrix)
% ======    -x position                                         (Integer)
%           -y position                                         (Integer)
%           -size                                               (Integer)
%
% Outputs : -Covariance Matrix                                  (7x7 Matrix)
% =======

function CR = findCovarianceMatrix(I, positionX, positionY, size)

    % First and second derivatives
    d = [-1 0 1];
    dd = [-1 2 -1];
    dI = double(I);

    % Ix, Iy, Ixx, Iyy
    Ix = conv2(d, dI);
    Iy = conv2(d,1,dI);
    Ixx = conv2(dd, dI);
    Iyy = conv2(dd,1,dI);

    for j = 1:size
        for i = 1:size

            f(i+(j-1)*size,:) = [positionX+i-1 positionY+j-1 dI(positionY+j-1, positionX+i-1) ...
                         Ix(positionY+j-1, positionX+i-1) Iy(positionY+j-1, positionX+i-1) ...
                        Ixx(positionY+j-1, positionX+i-1) Iyy(positionY+j-1,positionX+i-1)];

        end
    end

    % vector of means of features in region R
    uR = mean(f);

    T = zeros(7);
    for k = 1:size^2
        temp = (f(k,:)-uR)'*(f(k,:)-uR);
        T = T + temp;
    end

    CR = (1/size^2)*T;


text:
function D = addparts(D, mode)
% This function looks at the parts of the object and adds them inside the
% object struct. This useful when perfoming queries.
%
% D = addparts(D)
%
% The objects will have this new fields:
% D(i).annotation.object(j).parts.object(n).name
% D(i).annotation.object(j).parts.object(n).polygon
% D(i).annotation.object(j).parts.object(n).crop
% D(i).annotation.object(j).parts.object(n).occluded
%
% This function also allows inserting the object as a part of itself. This
% might be useful in some cases...
%
% D = addparts(D, 'addobjectasparts')

if nargin<2
    mode = '';
end

for i = 1:length(D)
    D(i).annotation = dealparts(D(i).annotation);

    if isfield(D(i).annotation, 'object')
        for j = 1:length(D(i).annotation.object)
            Nparts=0;
            % check that the field object is not inside already
            if ~isfield(D(i).annotation.object(j), 'parts') || ~isfield(D(i).annotation.object(j).parts, 'object')
                % chech that there are parts
                if isfield(D(i).annotation.object(j), 'parts') && isfield(D(i).annotation.object(j).parts, 'hasparts')
                    parts = D(i).annotation.object(j).parts.hasparts;
                    Nparts = length(parts);

                    % insert parts
                    for n = 1:Nparts
                        D(i).annotation.object(j).parts.object(n) = D(i).annotation.object(parts(n));
                    end
                end

                % if option 'addobjectasparts' then add object as part
                if strcmp(mode, 'addobjectasparts')
                    D(i).annotation.object(j).parts.object(Nparts+1).name = D(i).annotation.object(j).name;
                    D(i).annotation.object(j).parts.object(Nparts+1).occluded = D(i).annotation.object(j).occluded;
                    D(i).annotation.object(j).parts.object(Nparts+1).polygon = D(i).annotation.object(j).polygon;
                    D(i).annotation.object(j).parts.object(Nparts+1).crop = D(i).annotation.object(j).crop;
                end
            end
        end
    end
end

function annotation = dealparts(annotation)

if isfield(annotation, 'object') && isfield(annotation.object(1), 'parts')
    Nobjects = length(annotation.object);

    ids = [annotation.object.id];

    for n = 1:Nobjects
        parts = annotation.object(n).parts;
        if isfield(parts, 'ispartof')
            if isnumeric(annotation.object(n).parts.ispartof)
                [foo, ispartof] = ismember(annotation.object(n).parts.ispartof, ids);
            else
                [foo, ispartof] = ismember(str2num(annotation.object(n).parts.ispartof), ids);
            end
            annotation.object(n).parts.ispartof = ispartof;
        end

        if isfield(parts, 'hasparts')
            if isnumeric(annotation.object(n).parts.hasparts)
                [foo, hasparts] = ismember(annotation.object(n).parts.hasparts, ids);
            else
                [foo, hasparts] = ismember(str2num(annotation.object(n).parts.hasparts), ids);
            end
            hasparts = setdiff(hasparts, 0);
            annotation.object(n).parts.hasparts = hasparts;
        end
    end

end

author: CSAILVision

repo: LabelMeToolbox


text:
function varargout = conj(varargin)
%CONJ   Complex conjugate of a DISKFUN.
%   CONJ(F) returns the complex conjugate of F.  For a complex F, CONJ(F) =
%   REAL(F) - i*IMAG(F).
%
%   Since only real-valued DISKFUNS are presently supported, this
%   function is trivial.
%
% See also DISKFUN/REAL, DISKFUN/IMAG.

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

[varargout{1:nargout}] = conj@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [flow,labels] = test1()

% TEST1 Shows how to use the library to compute
%   a minimum cut on the following graph:
%
%                SOURCE
%		       /       \
%		     1/         \2
%		     /      3    \
%		   node0 -----> node1
%		     |   <-----   |
%		     |      4     |
%		     \            /
%		     5\          /6
%		       \        /
%		          SINK
%
%   (c) 2008 Michael Rubinstein, WDI R&D and IDC
%   $Revision: 140 $
%   $Date: 2008-09-15 15:35:01 -0700 (Mon, 15 Sep 2008) $
%

A = sparse(2,2);
A(1,2)=3;
A(2,1)=4;
T = sparse(2,2);
T(1,1)=1;
T(2,1)=2;
T(1,2)=5;
T(2,2)=6;

[flow,labels] = maxflow(A,T)


text:
function h = errorpatch(x,y,l,u,varargin)

% make inputs uniform
x = x(:);
y = y(:);
if numel(l) == 1,
  l = repmat(l,size(x));
else
  l = l(:);
end
if ~isnumeric(u),
  varargin = [{u},varargin];
  u = l;
elseif nargin < 4,
  u = l;
elseif numel(u) == 1,
  u = repmat(u,size(x));
else
  u = u(:);
end

% sort by x
[x,order] = sort(x);
y = y(order);
l = l(order);
u = u(order);

h = patch([x;flipud(x)],[y-l;flipud(y+u)],[.7,.7,.7],...
  'EdgeColor','none',varargin{:});

author: kristinbranson

repo: JAABA


text:
function k = multiKernDiagCompute(kern, x)

% MULTIKERNDIAGCOMPUTE Compute diagonal of MULTI kernel.
% FORMAT
% DESC computes the diagonal of the kernel matrix for the multiple output block kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : a vector containing the diagonal of the kernel matrix
% computed at the given points.
%
% SEEALSO : multiKernParamInit, kernDiagCompute, kernCreate, multiKernCompute
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% COPYRIGHT : Pei Gao, 2007

% MODIFICATIONS : Mauricio Alvarez, 2008

% KERN

if iscell(x)
  dim = 0;
  for i = 1:length(x)
    dim = dim + size(x{i},1);
  end
  k = zeros(dim, 1);
  startVal = 1;
  endVal = size(x{1},1);
  for i = 1:length(kern.comp)
      if ~isempty(x{i})
          k(startVal:endVal) = kernDiagCompute(kern.comp{i}, x{i});
      end
      startVal = endVal + 1;
      if i+1 <= length(kern.comp)
          endVal = endVal + size(x{i+1},1);
      end
  end
else
  k = zeros(size(x, 1)*kern.numBlocks, 1);
  endVal = size(x, 1);
  startVal = 1;
  for i = 1:length(kern.comp)
    k(startVal:endVal, 1)  = kernDiagCompute(kern.comp{i}, x);
    startVal = endVal + 1;
    endVal = endVal + size(x, 1);
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef TestBitwiseNot
    %TestBitwiseNot

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','img001.jpg');
    end

    methods (Static)
        function test_rgb_image
            img = cv.imread(TestBitwiseNot.im, 'ReduceScale',2);

            % rectangular mask
            [h,w,~] = size(img);
            mask = false([h w]);
            mask(100:h-100,100:w-100) = true;

            out = cv.bitwise_not(img);
            validateattributes(out, {class(img)}, {'size',size(img)});
            expected = my_bitwise_not(img);
            assert(isequal(out, expected));

            out = cv.bitwise_not(img, 'Mask',mask);
            validateattributes(out, {class(img)}, {'size',size(img)});
            expected = my_bitwise_not(img, mask);
            assert(isequal(out, expected));

            out = cv.bitwise_not(img, 'Mask',mask, 'Dest',img);
            validateattributes(out, {class(img)}, {'size',size(img)});
            expected = my_bitwise_not(img, mask, img);
            assert(isequal(out, expected));
        end

        function test_float_image
            img = cv.imread(TestBitwiseNot.im, 'ReduceScale',2);
            img = single(img) ./ 255;

            % circular mask
            [h,w,~] = size(img);
            [X,Y] = meshgrid(1:w,1:h);
            c = fix([w h]/2); r = 50;
            mask = ((X-c(1)).^2 + (Y-c(2)).^2) < r^2;

            out = cv.bitwise_not(img);
            validateattributes(out, {class(img)}, {'size',size(img)});
            expected = my_bitwise_not(img);
            assert(isequaln(out, expected));

            out = cv.bitwise_not(img, 'Mask',mask);
            validateattributes(out, {class(img)}, {'size',size(img)});
            expected = my_bitwise_not(img, mask);
            assert(isequaln(out, expected));

            out = cv.bitwise_not(img, 'Mask',mask, 'Dest',img);
            validateattributes(out, {class(img)}, {'size',size(img)});
            expected = my_bitwise_not(img, mask, img);
            assert(isequaln(out, expected));
        end

        function test_vector
            src = uint8([1 2 3]);
            mask = [true false true];

            % uint8([254 253 252])
            out = cv.bitwise_not(src);
            validateattributes(out, {class(src)}, {'size',size(src)});
            expected = bitcmp(src);
            assert(isequal(out, expected));

            % uint8([254 0 252])
            out = cv.bitwise_not(src, 'Mask',mask);
            validateattributes(out, {class(src)}, {'size',size(src)});
            expected = bitcmp(src);
            expected(~mask) = 0;
            assert(isequal(out, expected));

            % uint8([254 2 252])
            out = cv.bitwise_not(src, 'Mask',mask, 'Dest',src);
            validateattributes(out, {class(src)}, {'size',size(src)});
            expected = bitcmp(src);
            expected(~mask) = src(~mask);
            assert(isequal(out, expected));
        end

        function test_error_argnum
            try
                cv.bitwise_not();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

function out = my_bitwise_not(src, mask, dst)
    %MY_BITWISE_NOT  Similar to cv.bitwise_not using core MATLAB functions

    if nargin < 2, mask = true(size(src,1), size(src,2)); end
    if nargin < 3, dst = zeros(size(src), class(src)); end

    % calculate bitwise NOT
    if isinteger(src)
        out = bitcmp(src);
    elseif isfloat(src)
        if isa(src, 'double')
            klass = 'uint64';
        elseif isa(src, 'single')
            klass = 'uint32';
        end
        out = bitcmp(typecast(src(:), klass));
        out = reshape(typecast(out, class(src)), size(src));
    end

    % apply masking
    for k=1:size(out,3)
        out_slice = out(:,:,k);
        dst_slice = dst(:,:,k);
        out_slice(~mask) = dst_slice(~mask);
        out(:,:,k) = out_slice;
    end
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function SDUpdater(fhndl)

persistent app cleanup;

if ~isempty(cleanup)
    cleanup = [];
end

if ~ispc
    error('The SD Updater is only supported on Windows platforms');
end

app = feval(fhndl);
cleanup = onCleanup(@()delete(app));

author: analogdevicesinc


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

function snippets = generateSnippetsFromPlevels(plevel)

    % generate snippets using power levels

    snippets = struct;
    snippets.duration = {};
    snippets.mean = {};
    snippets.std = {};
    snippets.start = [];
    snippets.end = [];
    snippetStart = 1;
    while snippetStart < length(plevel.startidx) - 4
        snippetFound = 0;
        for snippetEnd = snippetStart + 1: snippetStart + 3
           if min(plevel.mean(snippetStart:snippetEnd)) + 10 < plevel.mean(snippetStart)
               break;
           end
           if abs(plevel.mean(snippetStart) - plevel.mean(snippetEnd+1)) < 10
               snippets.mean{end+1} = plevel.mean(snippetStart+1:snippetEnd) - plevel.mean(snippetStart);
               snippets.std{end+1} = plevel.std(snippetStart+1:snippetEnd);
               snippets.duration{end+1} = plevel.duration(snippetStart+1:snippetEnd);
               snippets.start(end+1) = snippetStart+1;
               snippets.end(end+1) = snippetEnd;
               snippetFound = 1;
               snippetStart = snippetEnd + 1;
               break;
           end
        end

        if snippetFound == 0
           snippetStart = snippetStart + 1;
        end
    end

end

author: beckel

repo: nilm-eval

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


text:
% ECGKIT
% Version 0.1 01-Oct-2014
% Files

%   InstallECGkit   - Installation script of ecg-kit
%   UnInstallECGkit - Uninstallation script of ecg-kit
%   ECGkit_examples - Example of how to use the ecg-kit

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function c=ref_dft_1(f)
%REF_DFT_1  Reference DFT by doubling
%   Usage:  c=ref_dft_1(f);
%
%   REF_DFT_1(f) computes a DFT of f by upsampling f and inserting zeros
%   at the odd positions.
%
%   This is not an efficient method, it is just meant to illustrate a
%   symmetry of the DFT.

L=size(f,1);
W=size(f,2);

flong=zeros(2*L,W,assert_classname(f));
flong(1:2:end-1)=f;

fflong=fft(flong)/sqrt(L);

c=fflong(1:L,:);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function logmmse(filename,outfile)

%
%  Implements the logMMSE algorithm [1].
%
%  Usage:  logmmse(noisyFile, outputFile)
%
%         infile - noisy speech file in .wav format
%         outputFile - enhanced output file in .wav format
%
%
%  Example call:  logmmse('sp04_babble_sn10.wav','out_log.wav');
%
%  References:
%   [1] Ephraim, Y. and Malah, D. (1985). Speech enhancement using a minimum
%       mean-square error log-spectral amplitude estimator. IEEE Trans. Acoust.,
%       Speech, Signal Process., ASSP-23(2), 443-445.
%
% Authors: Philipos C. Loizou
%
% Copyright (c) 2006 by Philipos C. Loizou
% $Revision: 0.0 $  $Date: 10/09/2006 $
%-------------------------------------------------------------------------

if nargin<2
   fprintf('Usage: logmmse(noisyfile.wav,outFile.wav) \n\n');
   return;
end

[x, Srate, bits]= wavread( filename);	%nsdata is a column vector

% =============== Initialize variables ===============

len=floor(20*Srate/1000); % Frame size in samples
if rem(len,2)==1, len=len+1; end;
PERC=50; % window overlap in percent of frame size
len1=floor(len*PERC/100);
len2=len-len1; % update rate in samples

win=hanning(len);  % define window
win = win*len2/sum(win);  % normalize window for equal level output

% Noise magnitude calculations - assuming that the first 6 frames is
% noise/silence

nFFT=2*len;
noise_mean=zeros(nFFT,1);
j=1;
for m=1:6
    noise_mean=noise_mean+abs(fft(win.*x(j:j+len-1),nFFT));
    j=j+len;
end
noise_mu=noise_mean/6;
noise_mu2=noise_mu.^2;

%--- allocate memory and initialize various variables

x_old=zeros(len1,1);
Nframes=floor(length(x)/len2)-floor(len/len2);
xfinal=zeros(Nframes*len2,1);

%===============================  Start Processing =======================================================
%
k=1;
aa=0.98;
mu=0.98;
eta=0.15;

ksi_min=10^(-25/10);

for n=1:Nframes

    insign=win.*x(k:k+len-1);

    spec=fft(insign,nFFT);
    sig=abs(spec); % compute the magnitude
    sig2=sig.^2;

    gammak=min(sig2./noise_mu2,40);  % limit post SNR to avoid overflows
    if n==1
        ksi=aa+(1-aa)*max(gammak-1,0);
    else
        ksi=aa*Xk_prev./noise_mu2 + (1-aa)*max(gammak-1,0);     % a priori SNR
        ksi=max(ksi_min,ksi);  % limit ksi to -25 dB
    end

    log_sigma_k= gammak.* ksi./ (1+ ksi)- log(1+ ksi);
    vad_decision= sum(log_sigma_k)/ len;
    if (vad_decision< eta)
        % noise only frame found
        noise_mu2= mu* noise_mu2+ (1- mu)* sig2;
    end
    % ===end of vad===

    A=ksi./(1+ksi);  % Log-MMSE estimator
    vk=A.*gammak;
    ei_vk=0.5*expint(vk);
    hw=A.*exp(ei_vk);

    sig=sig.*hw;
    Xk_prev=sig.^2;

    xi_w= ifft( hw .* spec,nFFT);
    xi_w= real( xi_w);

    xfinal(k:k+ len2-1)= x_old+ xi_w(1:len1);
    x_old= xi_w(len1+ 1: len);

    k=k+len2;

end

wavwrite(xfinal,Srate,16,outfile);

author: jtkim-kaist

repo: Speech-enhancement


text:
% ****************************************************************
% This is a function used for robotics project
%
% MAE 505: Robotic Final Project.
% Created By: Gan Tao
% Date: 10 Dec 2003.
% ****************************************************************

function [X_dot] = good_lang2(t, q)

% Xc                ***X coordinate of center mass of the mobile platform
% Yc                ***Y coordinate of center mass of the mobile platform
% fai               ***heading angle of the platform measured from the X axis of the world coordinates
% theta_R           ***angular displacement of the right wheel
% theta_L           ***angular displacement of the left wheel
% theta_dot_R       ***angular velocity of the right wheel??
% theta_dot_L       ***
% Xr                ***X coordinate of world coordinates of the reference point
% Yr                ***Y coordinate of world coordinates of the reference point
% Xr_c              ***
% Yr_c              ***
% c                 ***c=(r/2*b)
% b                 ***the distance between the driving wheels and the axis of symmetry
% d                 ***the distance from Po(the intersection of the axis of symmetry with the driving wheel axis) to Pc(the center
%                   ***of mass of the platform

Xr_c=5;
Yr_c=0;
c   =0.5;
b   =0.75;
d   =1;
Kp  =2;
Kd  =5;

Xc          = q(1); % xc
Yc          = q(2); % yc
fai         = q(3); % fai
theta_R     = q(4); % theta right
theta_L     = q(5); % theta left
theta_dot_R = q(6); % theta dot right
theta_dot_L = q(7); % theta dot left

% ****************** The path to trace:**************************************************************************

y_desire = t;       %line case ii x=y;
x_desire = t;

xy_dot_desire = [1; 1]; %2x1 matrix         line x=y case ii;

xy_2_dot_desire = [0; 0]; %2x1 matrix

xy_desire =[x_desire y_desire]';

%****************** To get the new input V(y2_dot) **************************************************************

Xr = Xc+Xr_c*cos(fai)-Yr_c*sin(fai);
Yr = Yc+Xr_c*sin(fai)+Yr_c*cos(fai);

xy = [Xr;Yr];

S = [c*(b*cos(fai)-d*sin(fai)) c*(b*cos(fai)+d*sin(fai));
     c*(b*sin(fai)+d*cos(fai)) c*(b*sin(fai)-d*cos(fai));
     c                         -c                       ;
     1                         0                        ;
     0                         1                        ]; %5*2;

faifunction = [ c*(b*cos(fai)-d*sin(fai))+(-Xr_c*sin(fai)-Yr_c*cos(fai))*c,c*(b*cos(fai)+d*sin(fai))-(-Xr_c*sin(fai)-Yr_c*cos(fai))*c;
                c*(b*sin(fai)+d*cos(fai))+(Xr_c*cos(fai)-Yr_c*sin(fai))*c,c*(b*sin(fai)-d*cos(fai))-(Xr_c*cos(fai)-Yr_c*sin(fai))*c ]; %2*2;

%*** faifunction_dot = diff(faifunction,fai);*****
niu = [theta_dot_R;theta_dot_L];

faifunction_dot =[ c*(-b*sin(fai)-d*cos(fai))+(-Xr_c*cos(fai)+Yr_c*sin(fai))*c,c*(-b*sin(fai)+d*cos(fai))-(-Xr_c*cos(fai)+Yr_c*sin(fai))*c;
             c*(b*cos(fai)-d*sin(fai))+(-Xr_c*sin(fai)-Yr_c*cos(fai))*c,c*(b*cos(fai)+d*sin(fai))-(-Xr_c*sin(fai)-Yr_c*cos(fai))*c ]; %2*2;

xy_dot  = faifunction*niu;%2*1;

xy2_dot = xy_2_dot_desire+Kd*(xy_dot_desire - xy_dot)+Kp*(xy_desire - xy);%2*1;

%**************** Using V(xy2_dot) to get the U input ************************************************************

U = inv(faifunction)*(xy2_dot-faifunction_dot*niu); %2*1

%**************** Using U to get the X_dot ***********************************************************************

f1 = S*niu;%5*1;

X_dot = [f1(1);f1(2);f1(3);f1(4);f1(5); 0; 0] + [0 0;0 0;0 0;0 0;0 0;1 0;0 1]*[U(1);U(2)];


text:
function [ f,qualMeasOut ] = ASD_POCS(proj,geo,angles,maxiter,varargin)
%ASD_POCS Solves the ASD_POCS total variation constrained image in 3D
% tomography.
%
%   ASD_POCS(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem
%   using the projection data PROJ taken over ALPHA angles, corresponding
%   to the geometry descrived in GEO, using NITER iterations.
%
%   ASD_POCS(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The
%   possible options in OPT are:
%
%
%   'lambda':      Sets the value of the hyperparameter for the SART iterations.
%                  Default is 1
%
%   'lambdared':   Reduction of lambda.Every iteration
%                  lambda=lambdared*lambda. Default is 0.99
%
%   'init':        Describes diferent initialization techniques.
%                   •  'none'     : Initializes the image to zeros (default)
%
%                   •  'FDK'      : intializes image to FDK reconstrucition
%   'TViter':      Defines the amount of TV iterations performed per SART
%                  iteration. Default is 20
%
%   'alpha':       Defines the TV hyperparameter. default is 0.002
%
%   'alpha_red':   Defines the reduction rate of the TV hyperparameter
%
%   'Ratio':       The maximum allowed image/TV update ration. If the TV
%                  update changes the image more than this, the parameter
%                  will be reduced.default is 0.95
%   'maxL2err'     Maximum L2 error to accept an image as valid. This
%                  parameter is crucial for the algorithm, determines at
%                  what point an image should not be updated further.
%                  Default is 20% of the FDK L2 norm.
%   'Verbose'      1 or 0. Default is 1. Gives information about the
%                  progress of the algorithm.
%
% 'OrderStrategy'  Chooses the subset ordering strategy. Options are
%                  'ordered' :uses them in the input order, but divided
%                  'random'  : orders them randomply
%                  'angularDistance': chooses the next subset with the
%                                     biggest angular distance with the ones used.
% 'redundancy_weighting': true or false. Default is true. Applies data
%                         redundancy weighting to projections in the update step
%                         (relevant for offset detector geometry)
%  'groundTruth'  an image as grounf truth, to be used if quality measures
%                 are requested, to plot their change w.r.t. this known
%                 data.
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri and Manasavee Lohvithee
%--------------------------------------------------------------------------

%% parse inputs
blocksize=1;
[beta,beta_red,f,ng,verbose,alpha,alpha_red,rmax,epsilon,OrderStrategy,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,varargin);

measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:)));
if ~any(isnan(gt(:)))
    QualMeasOpts{end+1}='error_norm';
    res_prev=gt;
    clear gt
end
if nargout<2 && measurequality
    warning("Image metrics requested but none catched as output. Call the algorithm with 3 outputs to store them")
    measurequality=false;
end
qualMeasOut=zeros(length(QualMeasOpts),maxiter);

[alphablocks,orig_index]=order_subsets(angles,blocksize,OrderStrategy);
angles_reorder=cell2mat(alphablocks);
index_angles=cell2mat(orig_index);

% does detector rotation exists?
if ~isfield(geo,'rotDetector')
    geo.rotDetector=[0;0;0];
end
%% Create weigthing matrices for the SART step
% the reason we do this, instead of calling the SART fucntion is not to
% recompute the weigths every ASD-POCS iteration, thus effectively doubling
% the computational time
% Projection weight, W
W=computeW(geo,angles,gpuids);

% Back-Projection weigth, V
V=computeV(geo,angles,alphablocks,orig_index,'gpuids',gpuids);

if redundancy_weights
    % Data redundancy weighting, W_r implemented using Wang weighting
    % reference: https://iopscience.iop.org/article/10.1088/1361-6560/ac16bc

    num_frames = size(proj,3);
    W_r = redundancy_weighting(geo);
    W_r = repmat(W_r,[1,1,num_frames]);
    % disp('Size of redundancy weighting matrix');
    % disp(size(W_r));
    W = W.*W_r; % include redundancy weighting in W
end

%%
stop_criteria=0;
iter=0;
offOrigin=geo.offOrigin;
offDetector=geo.offDetector;
rotDetector=geo.rotDetector;
DSD=geo.DSD;
DSO=geo.DSO;
while ~stop_criteria %POCS
    % If quality is going to be measured, then we need to save previous image
    if measurequality && ~strcmp(QualMeasOpts,'error_norm')
        res_prev = f; % only store if necesary
    end

    f0=f;
    if (iter==0 && verbose==1);tic;end
    iter=iter+1;

    for jj=1:size(angles,2)
        if size(offOrigin,2)==size(angles,2)
            geo.offOrigin=offOrigin(:,index_angles(:,jj));
        end
        if size(offDetector,2)==size(angles,2)
            geo.offDetector=offDetector(:,index_angles(:,jj));
        end
        if size(rotDetector,2)==size(angles,2)
            geo.rotDetector=rotDetector(:,index_angles(:,jj));
        end
        if size(DSD,2)==size(angles,2)
            geo.DSD=DSD(jj);
        end
        if size(DSO,2)==size(angles,2)
            geo.DSO=DSO(jj);
        end
        %         proj_err=proj(:,:,jj)-Ax(f,geo,angles(:,jj));          %                                 (b-Ax)
        %         weighted_err=W(:,:,jj).*proj_err;                   %                          W^-1 * (b-Ax)
        %         backprj=Atb(weighted_err,geo,angles(:,jj));            %                     At * W^-1 * (b-Ax)
        %         weigth_backprj=bsxfun(@times,1./V(:,:,jj),backprj); %                 V * At * W^-1 * (b-Ax)
        %         f=f+beta*weigth_backprj;                          % x= x + lambda * V * At * W^-1 * (b-Ax)
        % Enforce positivity
        f=f+beta* bsxfun(@times,1./V(:,:,jj),Atb(W(:,:,index_angles(:,jj)).*(proj(:,:,index_angles(:,jj))-Ax(f,geo,angles_reorder(:,jj),'gpuids',gpuids)),geo,angles_reorder(:,jj),'gpuids',gpuids));
        % non-negativity constrain
        if nonneg
            f=max(f,0);
        end
    end

    geo.offDetector=offDetector;
    geo.offOrigin=offOrigin;
    geo.DSD=DSD;
    geo.DSO=DSO;
    geo.rotDetector=rotDetector;
    % Save copy of image.
    if measurequality
        qualMeasOut(:,iter)=Measure_Quality(res_prev,f,QualMeasOpts);
    end
    % compute L2 error of actual image. Ax-b
    dd=im3Dnorm(Ax(f,geo,angles,'gpuids',gpuids)-proj,'L2');
    % compute change in the image after last SART iteration
    dp_vec=(f-f0);
    dp=im3Dnorm(dp_vec,'L2');

    if iter==1
        dtvg=alpha*dp;
        %Convert the steepest-descent step-size from a fraction of a
        %step-size to an absolute image distance on the first iteration.
    end
    f0=f;

    %  TV MINIMIZATION
    % =========================================================================
    %  Call GPU to minimize TV
    f=minimizeTV(f0,dtvg,ng,'gpuids',gpuids);    %   This is the MATLAB CODE, the functions are sill in the library, but CUDA is used nowadays
    %                                   for ii=1:ng
    % %                                 Steepest descend of TV norm
    %                                      tv(ng*(iter-1)+ii)=im3Dnorm(f,'TV','forward');
    %                                      df=weighted_gradientTVnorm2(f,0.002);
    %                                      df=df./im3Dnorm(df,'L2');
    %                                      f=f-dtvg.*df;
    %                                    end

    % update parameters
    % ==========================================================================

    % compute change by TV min
    dg_vec=(f-f0);
    dg=im3Dnorm(dg_vec,'L2');
    % if change in TV is bigger than the change in SART AND image error is still bigger than acceptable
    if dg>rmax*dp && dd>epsilon
        dtvg=dtvg*alpha_red;
    end
    % reduce SART step
    beta=beta*beta_red;
    % Check convergence criteria
    % ==========================================================================

    %Define c_alpha as in equation 21 in the journal
    c=dot(dg_vec(:),dp_vec(:))/(norm(dg_vec(:),2)*norm(dp_vec(:),2));
    %This c is examined to see if it is close to -1.0

    if (c<-0.99 && dd<=epsilon) || beta<0.005|| iter>=maxiter
        if verbose
            disp(['Stopping criteria met']);
            disp(['   c    = ' num2str(c),      '(Desired: c<-0.99)']);
            disp(['   beta = ' num2str(beta),   '(Desired: beta<0.005)']);
            disp(['   iter = ' num2str(iter), ]);
        end
        stop_criteria=true;
    end
    if (iter==1 && verbose==1)
        expected_time=toc*maxiter;
        disp('ADS_POCS');
        disp(['Expected duration   :    ',secs2hms(expected_time)]);
        disp(['Expected finish time:    ',datestr(datetime('now')+seconds(expected_time))]);
        disp('');
    end

end

end

function [beta,beta_red,f0,ng,verbose,alpha,alpha_red,rmax,epsilon,OrderStrategy,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,argin)

opts=     {'lambda','lambda_red','init','tviter','verbose','alpha','alpha_red','ratio','maxl2err','orderstrategy','qualmeas','nonneg','gpuids','redundancy_weighting','groundtruth'};
defaults=ones(length(opts),1);
% Check inputs
nVarargs = length(argin);
if mod(nVarargs,2)
    error('TIGRE:ASD_POCS:InvalidInput','Invalid number of inputs')
end

% check if option has been passed as input
for ii=1:2:nVarargs
    ind=find(ismember(opts,lower(argin{ii})));
    if ~isempty(ind)
        defaults(ind)=0;
    else
        error('TIGRE:ASD_POCS:InvalidInput',['Optional parameter "' argin{ii} '" does not exist' ]);
    end
end

for ii=1:length(opts)
    opt=opts{ii};
    default=defaults(ii);
    % if one option isnot default, then extract value from input
    if default==0
        ind=double.empty(0,1);jj=1;
        while isempty(ind)
            ind=find(isequal(opt,lower(argin{jj})));
            jj=jj+1;
        end
        if isempty(ind)
            error('TIGRE:ASD_POCS:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]);
        end
        val=argin{jj};
    end
    % parse inputs
    switch opt
        % Verbose
        %  =========================================================================
        case 'verbose'
            if default
                verbose=1;
            else
                verbose=val;
            end
            if ~is2014bOrNewer
                warning('TIGRE:Verbose mode not available for older versions than MATLAB R2014b');
                verbose=false;
            end
            % Lambda
            %  =========================================================================
            % Its called beta in ASD-POCS
        case 'lambda'
            if default
                beta=1;
            else
                if length(val)>1 || ~isnumeric( val)
                    error('TIGRE:ASD_POCS:InvalidInput','Invalid lambda')
                end
                beta=val;
            end
            % Lambda reduction
            %  =========================================================================
        case 'lambda_red'
            if default
                beta_red=0.99;
            else
                if length(val)>1 || ~isnumeric( val)
                    error('TIGRE:ASD_POCS:InvalidInput','Invalid lambda')
                end
                beta_red=val;
            end
            % Initial image
            %  =========================================================================
        case 'init'
            if default || strcmp(val,'none')
                f0=zeros(geo.nVoxel','single');

            else
                if strcmp(val,'FDK')
                    f0=FDK(proj, geo, angles);
                else
                    error('TIGRE:ASD_POCS:InvalidInput','Invalid init')

                end
            end
            % Number of iterations of TV
            %  =========================================================================
        case 'tviter'
            if default
                ng=20;
            else
                ng=val;
            end
            %  TV hyperparameter
            %  =========================================================================
        case 'alpha'
            if default
                alpha=0.002; % 0.2
            else
                alpha=val;
            end
            %  TV hyperparameter redution
            %  =========================================================================
        case 'alpha_red'
            if default
                alpha_red=0.95;
            else
                alpha_red=val;
            end
            %  Maximum update ratio
            %  =========================================================================
        case 'ratio'
            if default
                rmax=0.95;
            else
                rmax=val;
            end
            %  Maximum L2 error to have a "good image"
            %  =========================================================================
        case 'maxl2err'
            if default
                epsilon=im3Dnorm(Ax(FDK(proj,geo,angles),geo,angles)-proj,'L2')*0.2; %heuristic
            else
                epsilon=val;
            end
            %  Order strategy
            %  =========================================================================
        case 'orderstrategy'
            if default
                OrderStrategy='random';
            else
                OrderStrategy=val;
            end
            % Image Quality Measure
            %  =========================================================================
        case 'qualmeas'
            if default
                QualMeasOpts={};
            else
                if iscellstr(val)
                    QualMeasOpts=val;
                else
                    error('TIGRE:ASD_POCS:InvalidInput','Invalid quality measurement parameters');
                end
            end
            %  Non negative
            %  =========================================================================
        case 'nonneg'
            if default
                nonneg=true;
            else
                nonneg=val;
            end
            %  GPU Ids
            %  =========================================================================
        case 'gpuids'
            if default
                gpuids = GpuIds();
            else
                gpuids = val;
            end
            % Data Redundancy weighting
        case 'redundancy_weighting'
            if default
                redundancy_weights = true;
            else
                redundancy_weights = val;
            end
        case 'groundtruth'
            if default
                gt=nan;
            else
                gt=val;
            end
        otherwise
            error('TIGRE:ASD_POCS:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in ASD_POCS()']);

    end
end

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE