text:
% Standard text for function documentation
%
% First line: One-line summary description of function
%
% :Usage:
% ::
%
%     [list outputs here] = function_name(list inputs here, [optional inputs])
%
% For objects: Type methods(object_name) for a list of special commands
%              Type help object_name.method_name for help on specific
%              methods.
%
% ..
%     Author and copyright information:
%
%     Copyright (C) <year>  <name of author>
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
% ..
%
% :Inputs:
%
%   **param1:**
%        description of param1
%
%   **param2:**
%        description of param2
%
% :Optional Inputs:
%   **param1:**
%        description of param1
%
%   **param2:**
%        description of param2
%
% :Outputs:
%
%   **out1:**
%        description of out1
%
%   **out2:**
%        description of out2
%
%   **'doplot', [logical flag]:**
%        Create plots; default = true. 'noplot' to turn off.
%
%   **'doverbose', [logical flag]:**
%        Verbose output; default = true. 'noverbose' to turn off.
%
%   **'dooptimize', [logical flag]:**
%        Optimize hyperparameters; default = true. 'nooptimize' to turn off.
%
%   **'doprepeats', [num repeats]:**
%        Repeat cross-val with different partitions; default = 10.
%        Enter number of repeats, 'norepeats' to turn off.
%
% :Examples:
% ::
%
%    % give examples of code here
%    param1 = abc();
%    param2 = xyz();
%    [out1,out2] = func_call(param1, param2)
%
% :References:
%   CITATION(s) HERE
%
% :See also:
%   - list other functions related to this one, and alternatives*
%

% ..
%    Programmers' notes:
%    List dates and changes here, and author of changes
% ..

% BELOW IS A STANDARD TEMPLATE FOR DEFINING VARIABLE (OPTIONAL) INPUT
% ARGUMENTS. MANY FUNCTIONS NEED TO PARSE OPTIONAL ARGS, SO THIS MAY BE
% USEFUL. I (TOR) RECOMMEND COPYING AND PASTING THE DOCUMENTATION ABOVE AND
% ALL THE STUFF THROUGH ARGUMENT VALIDATION, BELOW, AND MODIFYING IT FOR
% YOUR FUNCTION.

% SEE BELOW FOR TEMPLATE FOR OBJECT CONSTRUCTORS. THIS HAS STRUCTURED
% DOCUMENTATION WITH SOME TEXT FORMATTING RULES THAT WILL MAKE IT
% COMPATIBLE WITH MATLAB'S DOC COMMAND, WHICH WILL THEN CREATE AN HTML PAGE
% LISTING METHODS AND ATTRIBUTES. I (TOR) RECOMMEND COPYING AND PASTING THE
% DOCUMENTATION BELOW AND MODIFYING IT FOR YOUR FUNCTION.

% ..
%    DEFAULTS AND INPUTS
% ..

% Some useful display helper functions
% --------------------------------------------------------

n_cols = 80;
sep_str = repmat('_', 1, n_cols);  % see textwrap

dashes = '----------------------------------------------';
printstr = @(dashes) disp(dashes);
printhdr = @(str) fprintf('%s\n%s\n%s\n', dashes, str, dashes);

% Enforce valid variable names in a cell array of strings: Eliminate special characters and leading numbers
[variable_names, namewarnings] = format_text_letters_only(variable_names, 'numbers', 'cleanup', 'squeeze', 'underscore_ok');

% ----------------------------------------------------------------------
% Parse inputs
% ----------------------------------------------------------------------
% This 2019 version uses the inputParser object. Older schemes are below.
% Note: With this, you can pass in EITHER keyword, value pairs OR a
% structure with fields that define keywords.
% Note: You can also make this a subfunction at the end of your function,
% e.g., function ARGS = parse_inputs(varargin)...end, and call it using:

ARGS = parse_inputs(varargin{:});

% If you want to distribute arguments back out to variables, use this:

fn = fieldnames(ARGS);

for i = 1:length(fn)
    str = sprintf('%s = ARGS.(''%s'');', fn{i}, fn{i});
    eval(str)
end

% Logical flags
% ----------------------------------------------------------------------
if any(strcmp(varargin, 'noverbose')), doverbose = false; end
if any(strcmp(varargin, 'noplots')), doplots = false; end

% ======== This goes at the end of the file / after the main function ========
function ARGS = parse_inputs(varargin)

p = inputParser;

% Validation functions - customized for each type of input
% ----------------------------------------------------------------------

valfcn_scalar = @(x) validateattributes(x, {'numeric'}, {'nonempty', 'scalar'});

valfcn_number = @(x) validateattributes(x, {'numeric'}, {'nonempty'}); % scalar or vector

% Validation: Region object, structure, or [x1 x2 x3] triplet
valfcn_custom = @(x) isstruct(x) || isa(x, 'region') || (~isempty(x) && all(size(x) - [1 3] == 0) && all(isnumeric(x)));

% Validation: [x1 x2 x3] triplet
valfcn_xyz = @(x) validateattributes(x, {'numeric'}, {'nonempty', 'size', [1 3]});

valfcn_logical = @(x) validateattributes(x, {'numeric'}, {'nonempty', 'scalar', '>=', 0, '<=', 1}); % could enter numeric 0,1 or logical

% Required inputs
% ----------------------------------------------------------------------
p.addRequired('x', valfcn_custom);
p.addRequired('y', valfcn_custom);

% Optional inputs
% ----------------------------------------------------------------------
% Pattern: keyword, value, validation function handle

p.addParameter('color', [.9 .2 0], valfcn_xyz);
p.addParameter('bendpercent', .15, valfcn_number); % can be scalar or vector
p.addParameter('thickness', .1, valfcn_scalar);
p.addParameter('nstreamlines', 30, valfcn_scalar);

% Parse inputs and distribute out to variable names in workspace
% ----------------------------------------------------------------------
% e.g., p.parse([30 1 0], [-40 0 10], 'bendpercent', .1);
p.parse(varargin{:});

ARGS = p.Results;

end % parse_inputs

% ==================== end pattern =====================================

% Next is the previously used scheme for CANlab tools for several years

% -------------------------------------------------------------------------
% DEFAULT ARGUMENT VALUES
% -------------------------------------------------------------------------

rowsz = [];
doplot = 0;
basistype = 'spm+disp';
% initalize optional variables to default values here.

% -------------------------------------------------------------------------
% OPTIONAL INPUTS
% -------------------------------------------------------------------------

% This is a compact way to assign multiple variables. The input argument
% names and variable names must match, however:

allowable_inputs = {'var_names' 'p_thr' 'dospearman' 'dopartial' 'dofdr' 'dofigure' 'doimage' 'docircles' 'dotext' 'colorlimit' 'text_x_offset' 'text_y_offset' 'text_fsize' 'text_nonsig_color' 'text_sig_color' 'partitions' 'partitioncolors' 'partitionlabels'};

keyword_inputs = {'noorthviews'};

% optional inputs with default values - each keyword entered will create a variable of the same name

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case allowable_inputs

                eval([varargin{i} ' = varargin{i+1}; varargin{i+1} = [];']);

            case keyword_inputs
                % Skip deal with this below

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

% 2nd pass: Keyword inputs. These supersede earlier inputs
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case 'writetofile'
                writetofile = true;
                if isempty(movieoutfile)
                    movieoutfile = fullfile(pwd, 'rmssd_movie_file');
                    fprint('Writing file with default name:\n%s\n', movieoutfile);
                end

            case {'nomovie' 'nodisplay'}
                showmovie = false;

        end
    end
end

% 2nd pass: Keyword inputs. These supersede earlier inputs

if any(strcmp(varargin{i}, 'writetofile'))

    writetofile = true;
    if isempty(movieoutfile)
        movieoutfile = fullfile(pwd, 'rmssd_movie_file');
        fprint('Writing file with default name:\n%s\n', movieoutfile);
    end

end

if any(strcmp(varargin{i}, 'nomovie')) || any(strcmp(varargin{i}, 'nodisplay'))

    showmovie = false;

end

% This pattern will flexibly assign arguments based on keywords.
% The names of the input keyword and variable created do not need to match.
% Multiple input keywords can be mapped to the same variable.

% optional inputs with default values
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case {'rows', 'rowsz'}, rowsz = varargin{i+1}; varargin{i+1} = [];
            case 'plot', doplot = 1;
            case 'basistype', basistype = varargin{i+1}; varargin{i+1} = [];

            case allowable_inputs
                % skip - handled above

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

% -------------------------------------------------------------------------
% VALIDATE ATTRIBUTES OF INPUTS
% -------------------------------------------------------------------------

validateattributes(X,{'numeric'},{'2d'},'plot_correlation_matrix','X', 1);

logical_args = {'dofdr' 'false' 'dospearman' 'dofigure' 'doimage' 'docircles' 'dotext'};
for i = 1:length(logical_args)

    my_arg = eval([logical_args{i} ';']);
    validateattributes(my_arg,{'logical'},{'scalar'},'plot_correlation_matrix',logical_args{i});

end

cell_args = {'partitioncolors' 'partitionlabels' 'var_names'};
for i = 1:length(cell_args)

    my_arg = eval([cell_args{i} ';']);
    validateattributes(my_arg,{'cell'},{},'plot_correlation_matrix',cell_args{i});

end

vector_args = {'text_sig_color' 'text_nonsig_color'};
for i = 1:length(vector_args)

    my_arg = eval([vector_args{i} ';']);
    validateattributes(my_arg,{'numeric' 'vector' 'nonnegative' 'nonnan'},{},'plot_correlation_matrix',vector_args{i});

end

supportedImageClasses    = {'int8','uint8','int16','uint16','int32','uint32','single','double','logical'};
supportedImageAttributes = {'real','nonsparse','nonempty'};
validateattributes(V,supportedImageClasses,supportedImageAttributes,mfilename,'V');

% SEE:
%  validateattributes Check validity of array.
%     validateattributes(A,CLASSES,ATTRIBUTES) validates that array A belongs
%     to at least one of the specified CLASSES and has all of the specified
%     ATTRIBUTES.
%
% See also inputParser.m

% ==================== end pattern =====================================

% THIS IS ANOTHER WAY THAT IS MORE STREAMLINED, BUT ACTUALLY LESS INTUITIVE

% -------------------------------------------------------------------------
% DEFAULTS AND INPUTS
% -------------------------------------------------------------------------

% Defaults
% -----------------------------------
% Set color maps for + / - values
poscm = colormap_tor([1 0 .5], [1 1 0], [.9 .6 .1]);  %reddish-purple to orange to yellow
negcm = colormap_tor([0 0 1], [0 1 1], [.5 0 1]);  % cyan to purple to dark blue

% optional inputs with default values
% -----------------------------------
% - allowable_args is a cell array of argument names
% - avoid spaces, special characters, and names of existing functions
% - variables will be assigned based on these names
%   i.e., if you use an arg named 'cl', a variable called cl will be
%   created in the workspace

allowable_args = {'cl', 'ycut_mm', 'pos_colormap', 'neg_colormap', ...
    'surface_handles', 'existingfig'};

default_values = {[], [], poscm, negcm, ...
    [], 0};

% define actions for each input
% -----------------------------------
% - cell array with one cell for each allowable argument
% - these have special meanings in the code below
% - allowable actions for inputs in the code below are: 'assign_next_input' or 'flag_on'

actions = {'assign_next_input', 'assign_next_input', 'assign_next_input', 'assign_next_input', ...
    'assign_next_input', 'flag_on'};

% logical vector and indices of which inputs are text
textargs = cellfun(@ischar, varargin);
whtextargs = find(textargs);

for i = 1:length(allowable_args)

    % assign default
    % -------------------------------------------------------------------------

    eval([allowable_args{i} ' =  default_values{i};']);

    wh = strcmp(allowable_args{i}, varargin(textargs));

    if any(wh)
        % Optional argument has been entered
        % -------------------------------------------------------------------------

        wh = whtextargs(wh);
        if length(wh) > 1, warning(['input ' allowable_args{i} ' is duplicated.']); end

        switch actions{i}
            case 'assign_next_input'
                eval([allowable_args{i} ' = varargin{wh(1) + 1};']);

            case 'flag_on'
                eval([allowable_args{i} ' = 1;']);

            otherwise
                error(['Coding bug: Illegal action for argument ' allowable_args{i}])
        end

    end % argument is input
end

% END DEFAULTS AND INPUTS
% -------------------------------------------------------------------------

% Below are some default headers

% -------------------------------------------------------------------------
% -------------------------------------------------------------------------
%
% Sub-functions
%
% -------------------------------------------------------------------------
% -------------------------------------------------------------------------

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

% OBJECT FUNCTION DEFINITION - followed by one-line description
% image_vector: An object that allows for storage and manipulation of neuroimaging data
%
% -------------------------------------------------------------------------
% Features and philosophy:
% -------------------------------------------------------------------------
%
% - Store image data in a flat (2-d), space-efficient voxels x images matrix
% - Save space: Tools to remove out-of-image voxels and empty (zero/NaN) voxels, store data in single-precision format
% - Analysis-friendly: 2-d matrices can be read and analyzed in multiple packages/algorithms
% - Meta-data included to convert back to 3-d image volume space,
%   with easy tools (methods) to reconstruct and visualize images
%   Built-in resampling makes it easy to compare/combine datasets with different voxel sizes and image bounding boxes
%   Reduces overhead for statisticians/data scientists unfamiliar with neuroimaging to apply their algorithms.
% - Multiple images can be stored in a single object
% - Methods have short, intuitive names, and perform high-level functions specialized for neuroimaging:
% - Visualization (plot, orthviews, surface, montage, histogram, isosurface methods)
% - Image manipulation (apply_mask, get_wh_image, resample_space, compare_space, flip, threshold methods)
% - Data extraction (apply_atlas, apply_parcellation, extract_gray_white_csf, extract_roi_averages)
% - Analysis (ica, mahal, image_math, and many more in the fmri_data subclass)
% - Provenance: Ability to track and update history of changes to objects
%
% image_vector is a superclass of several other object types:
% fmri_data
% statistic_image
% atlas
%
% These object types have additional analysis methods and properties, and
% are most often used. You will rarely create an image_vector object
% directly. Typically you would create instances of one of the classes above.
%
% -------------------------------------------------------------------------
% To construct/create a new instance of an object:
% -------------------------------------------------------------------------
%
% As with any neuroimaging-data object class in this toolbox
% (image_vector, fmri_data, statistic_image, atlas), you can create an object by
% specifying the names of a set of images you want to load in and passing them into the
% function with the same name as the object class (e.g., fmri_data). This
% will create an object with image data and image space information attached.
%
% Use fmri_data instead of image_vector to load image data:
% data_obj = fmri_data(image_names);
%
% ...where image_names is a character array or 1 x n cell array of strings
%
% -------------------------------------------------------------------------
% Properties and methods
% -----------------------------------------------------------------------
% Properties are data fields associated with an object.
% Type the name of an object (class instance) you create to see its
% properties, and a link to its methods (things you can run specifically
% with this object type). For example: After creating an fmri_data object
% called fmri_dat, as above, type fmri_dat to see its properties.
%
% There are many other methods that you can apply to fmri_data objects to
% do different things.
% - Try typing methods(fmri_data) for a list.
% - You always pass in an fmri_data object as the first argument.
% - Methods include utilities for many functions - e.g.,:
% - resample_space(fmri_dat) resamples the voxels
% - write(fmri_dat) writes an image file to disk (careful not to overwrite by accident!)
% - regress(fmri_dat) runs multiple regression
% - predict(fmri_dat) runs cross-validated machine learning/prediction algorithms
%
% Key properties and methods (a partial list; type doc image_vector for more):
% -------------------------------------------------------------------------
% image_vector Properties:
%   dat                     - Image data, a [voxels x images] matrix, single-format
%   fullpath                - List of image names loaded into object with full paths
%   history                 - History of object processing, for provenance
%   image_names             - List of image names loaded into object, no paths
%   removed_images          - Vector of images that have been removed (saves space; see remove_empty.m, replace_empty.m)
%   removed_voxels          - Vector of empty in-mask voxels that have been removed (saves space; see remove_empty.m, replace_empty.m)
%   volInfo                 - Structure with info on brain mask (saves space) and mapping voxels to brain space
%
% image_vector Methods:
%   General:
% . 	descriptives        -  Get descriptives for an fmri_data or other image_vector object
%       enforce_variable_types	-  Re-casts variables in objects into standard data types, which can save
%       flip                - Flips images stored in an object left-to-right
%     	history             - Display history for image_vector object
%       write               - Write an image_vector object to hard drive as an Analyze image (uses .fullpath field for image names)
%   Data extraction:
%       apply_atlas         - Computes the mean value or pattern expression for each reference region specified in an atlas object
%       apply_mask          - Apply a mask image (image filename or fmri_mask_image object) to an image_vector object
%       apply_parcellation  - Computes the mean value or pattern expression for each parcel specified in a data object
%       extract_gray_white_csf	- Extracts mean values (values) and top 5 component scores (components)
%       extract_roi_averages	- This image_vector method a extracts and averages data stored in an fmri_data object
%
%   Handling brain space and image selection:
%       compare_space           - Compare spaces of two image_vector objects
%       get_wh_image            - For an image_vector with multiple images (cases, contrasts, etc.), select a subset.
%       reconstruct_image       - Reconstruct a 3-D or 4-D image from image_vector object obj
%       remove_empty            - remove vox: logical vector of custom voxels to remove, VOX x 1
%       reparse_contiguous      - Re-construct list of contiguous voxels in an image based on in-image
%       replace_empty           - Replace empty/missing values in an image data object
%       resample_space          - Resample the images in an fmri_data object (obj) to the space of another
%
%   Display and visualization:
%   	display_slices      - Creates 3 separate montage views - ax, cor, sagg in a special figure window
%       histogram           - Create a histogram of image values or a series of histograms for each
%       image_similarity_plot - Associations between images in object and set of 'spatial basis function' images (e.g., 'signatures' or pre-defined maps)
%       isosurface          - Create and visualize an isosurface created from the boundaries in an image object.
%       montage             - Create a montage of an image_vector (or statistic_image or fmri_data)
%       orthviews               - display SPM orthviews for CANlab image_vector (or fmri_data, statistic_image) object
%       pattern_surf_plot_mip	- axial maximum intensity projection pattern surface plot
%   	sagg_slice_movie	- Movie of successive differences (sagittal slice)
%       slices              - Create a montage of single-slice results for every image in an image_vector object
%       surface             - Render image data on brain surfaces; options for cutaways and canonical surfaces
%       wedge_plot_by_atlas	- Plot a data object or 'signature' pattern divided into local regions
%
%   Data processing and analysis:
%   	ica                 - Spatial ICA of an fmri_data object
%       image_math          - Perform simple mathematical and boolean operations on image objects (see also plus, minus, power)
%       mahal               - Mahalanobis distance for each image in a set compared to others in the set
%       mean                - Mean across a set of images. Returns a new image_vector object.
%       preprocess          - Preprocesses data in an image_vector (e.g., fmri_data) object; many options for filtering and outlier id
%       qc_metrics_second_level	- Quality metrics for a 2nd-level analysis (set of images from different subjects)
%       searchlight         - Run searchlight multivariate prediction/classification on an image_vector
%       threshold           - Threshold image_vector (or fmri_data or fmri_obj_image) object based on raw threshold values
%       union               - ...and intersection masks for two image_vector objects
%
% -------------------------------------------------------------------------
% Examples and help:
% -------------------------------------------------------------------------
%
% To list properties and methods for this object, type:
% doc image_vector, methods(image_vector)
%
% Example 1:
% Load a sample dataset into an fmri_data object (subclass of image_vector)
% This loads one of a set of named image collections used in demos/help:
% data_obj = load_image_set('emotionreg');
%
% You can load the same images manually, by locating the files, listing
% their names in a character array (or 1 x n cell array of strings), and
% then passing those into fmri_data:
%
% filedir = what(fullfile('CanlabCore', 'Sample_datasets', 'Wager_et_al_2008_Neuron_EmotionReg'));
% image_names = filenames(fullfile(filedir.path, '*img'));
% data_obj = fmri_data(image_names);
%
% Now you can interact with the object.  Try, e.g.,:
% methods(data_obj)                               % List methods for object type
% plot(data_obj)                                  % Custom fmri_data specific plots
% t = ttest(data_obj);                            % Perform a voxel-wise one-sample t-test across images
% t = threshold(t, .005, 'unc', 'k', 10);         % Re-threshold with extent threshold of 10 contiguous voxels
% r = region(t);                                  % Turn t-map into a region object with one element per contig region
%
% For more examples and walkthroughs, see the CANlab_help_examples
% repository at https://github.com/canlab/CANlab_help_examples
%
% Some example tutorials:
% canlab_help_1_installing_tools
% canlab_help_2_load_a_sample_dataset
% canlab_help_3_voxelwise_t_test_walkthrough
% canlab_help_4_write_data_to_image_file_format
% canlab_help_5_regression_walkthrough
%
% License and authors:
% -------------------------------------------------------------------------
% By Tor Wager and CANlab members and collaborators
% GNU General Public License; see http://www.gnu.org/licenses/
% See Programmers' Notes in code for more info
%
%
% -------------------------------------------------------------------------
% See also:
% fmri_data
% statistic_image
% atlas
% region


text:
function str=civilizedStringFromCellArrayOfStrings(c,separator)

if ~exist('separator','var')
  separator=',';
end

n=length(c);
if n==0
  str='';
elseif n==1
  str=c{1};
elseif n==2
  str=[c{1} ' and ' c{2}];
else
  str='';
  for i=1:n-1
    str=[str c{i} separator ' '];
  end
  str=[str 'and ' c{n}];
end

end

author: kristinbranson

repo: JAABA


text:
function [hr] = yr2hr(yr)
% Convert time from Julian years (365.25 days) to hours.
% Chad Greene 2012
hr = yr*8766;

author: Sable

repo: mcbench-benchmarks


text:
% add/append a Data Quality struct to an ERPset
% part of the ERPLAB Toolbox
% axs May 2019
%
%
function ERP = add_dq_measure(ERP, DQ_struct, name)

% check input
assert(isstruct(ERP));

% If DQ is currently absent
if isfield(ERP,'dataquality') == 0

    ERP.dataquality(1).type = 'empty';
    ERP.dataquality(1).times = [];
    ERP.dataquality(1).data = [];
    ERP.dataquality(1).time_window_labels = {};
    ERP.dataquality(1).comments = [];
end

% with no DQ_struct specified, populate with empty
if exist('DQ_struct','var') == 0 || isempty(DQ_struct)
    DQ_struct.type = 'empty';
    DQ_struct.times = [];
    DQ_struct.data = [];
    DQ_struct.time_window_labels = {};
    DQ_struct.comments = [];
end

if exist('name','var') == 0 || isempty(name)
    name = DQ_struct.type;
end

%

n_existing_dq = length(ERP.dataquality);
was_empty = strcmpi(ERP.dataquality(n_existing_dq).type,'empty');

if was_empty
    dq_slot = n_existing_dq;
else
    dq_slot = n_existing_dq + 1;
end

ERP.dataquality(dq_slot) = DQ_struct;

[ERP, serror] = sorterpstruct(ERP);

if serror
    warning('ERP sorting error in making DQ struct')
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
addpath(genpath('BasicFuncPano'));
addpath(genpath('DDSampling'));
addpath(genpath('GroundTruth'));
addpath(genpath('ObjectHypothesisGeneration'));
addpath(genpath('RoomHypothesisSampling'));
addpath(genpath('Toolbox'));
addpath(genpath('VpEstimation'));
addpath(genpath('ResultVisualization'));

author: zouchuhang

repo: LayoutNet


text:
% readneurodat() - read neuroscan location file (.dat)
%
% Usage:
%   >> [ CHANLOCS labels theta theta ] = readneurodat( filename );
%
% Inputs:
%   filename       - file name or matlab cell array { names x_coord y_coord }
%
% Outputs:
%   CHANLOCS       - [structure] EEGLAB channel location data structure. See
%                    help readlocs()
%   labels         - [cell arrya] channel labels
%   theta          - [float array]array containing 3-D theta angle electrode
%                    position (in degree)
%   phi            - [float array]array containing 3-D phi angle electrode
%                    position (in degree)
%
% Author: Arnaud Delorme, CNL / Salk Institute, 28 Nov 2003
%
% See also: readlocs(), readneurolocs()

% Copyright (C) Arnaud Delorme, CNL / Salk Institute, 15 March 2003, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [chanlocs, labels, theta, phi] = readneurodat(filename)

    if nargin < 1
        help readneurodat;
        return;
    end

    % enter file name here
    % --------------------
    %tmp = loadtxt('/home/ftp/pub/locfiles/neuroscan/cap128.dat');
    %tmp = loadtxt('/home/arno/temp/quik128.DAT');
    tmp = loadtxt(filename);

    % resort electrodes
    % -----------------
    [tmp2, tmpind] = sort(celltomat(tmp(:,1))');
    tmp = tmp(tmpind,:);

    % convert to polar coordinates
    % ----------------------------
    %figure; plot(celltomat(tmp(:,2)), celltomat(tmp(:,3)), '.');
    [phi,theta] = cart2pol(celltomat(tmp(:,end-1)), celltomat(tmp(:,end)));
    theta = theta/513.1617*44;
    phi   = phi/pi*180;

    % convert to other types of coordinates
    % -------------------------------------
    labels = tmp(:,end-2)';
    numLabel = labels{1};
    if isnumeric(numLabel)
        if numLabel ~= round(numLabel) % floating point label
            error('Wrong file format')
        end
    end
    chanlocs = struct('labels', labels, 'sph_theta_besa', mattocell(theta)', 'sph_phi_besa', mattocell(phi)');
    chanlocs = convertlocs( chanlocs, 'sphbesa2all');
    for index = 1:length(chanlocs)
        chanlocs(index).labels = num2str(chanlocs(index).labels);
    end
    theta = theta/pi*180;
    fprintf('Note: .dat file contain polar 2-D coordinates. EEGLAB will use these coordinates\n');
    fprintf('      to recreated 3-D coordinates.\n');
    fprintf('\n');
    fprintf('Warning: if the channels are clustered in the middle or oddly spaced removed\n');
    fprintf('         the peripheral channels and then used the optimize head function\n');
    fprintf('         in the Channel Locations GUI. It seems that sometimes the peripherals\n');
    fprintf('         are throwing off the spacing.\n');

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
classdef S2Grid < vector3d
%
% The class S2Grid represent spherical grids. The central difference to a
% simple list of @vector3d is that S2Grid provides more efficient methods
% for finding the closes points in the grid with respect to a given
% direction.
%
% Syntax
%
%   % regular grid with fixed polar and azimuthal spacing
%   S2G = regularS2Grid('resolution',5*degree)
%   S2G = regularS2Grid('theta',(0:5:80)*degree,'rho',(0:5:355)*degree)
%
%   % regular grid optimized for plotting
%   S2G = plotS2Grid('resolution',1.5*degree,'upper')
%
%   % equispaced nodes  with given resolution
%   S2G = equispacedS2Grid('resolution',5*degree)
%
% Options
%  points     - [nrho,ntheta] number of points
%  resolution - resolution of a equispaced grid
%  theta      - polar angle
%  rho        - azimuthal angle
%  minRho     - starting rho angle (default 0)
%  maxRho     - maximum rho angle (default 2*pi)
%  minTheta   - starting theta angle (default 0)
%  maxTheta   - maximum theta angle (default pi)
%
% Flags
%  lower, uper, complete - restrict hemisphere
%  antipodal  - include <VectorsAxes.html antipodal symmetry>
%  no_center  - ommit point at center
%  restrict2minmax - restrict margins to min / max
%
% Examples
%
%   S2G = equispacedS2Grid('resolution',5*degree,'maxTheta',70*degree)
%   plot(S2G)
%
% See also
% vector3d.vector3d plotS2Grid regularS2Grid equispaceS2Grid

properties

  thetaGrid = S1Grid([],0,pi);
  rhoGrid = S1Grid([],0,2*pi);
  res = 2*pi;

end

methods

  function S2G = S2Grid(thetaGrid,rhoGrid,varargin)

    % call superclass method
    v = calcGrid(thetaGrid,rhoGrid);
    [S2G.x,S2G.y,S2G.z] = double(v);
    S2G.thetaGrid = thetaGrid;
    S2G.rhoGrid = rhoGrid;
    S2G.res = get_option(varargin,'resolution',2*pi);
    S2G.antipodal = check_option(varargin,'antipodal');
  end

  function v = vector3d(S2G)
    v = vector3d(S2G.x,S2G.y,S2G.z);
    v.opt = S2G.opt;
    v.antipodal = S2G.antipodal;
  end

end

end

author: mtex-toolbox

repo: mtex

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


text:
function A = makehatch(hatch)
%MAKEHATCH Predefined hatch patterns
%  MAKEHATCH(HATCH) returns a matrix with the hatch pattern for HATCH
%   according to the following table:
%      HATCH        pattern
%     -------      ---------
%        /          right-slanted lines
%        \          left-slanted lines
%        |          vertical lines
%        -          horizontal lines
%        +          crossing vertical and horizontal lines
%        x          criss-crossing lines
%        .          single dots
%
%  See also: APPLYHATCH

%  Copyright 2009 The MathWorks, Inc.

n = 6;
A=zeros(n);
switch (hatch)
 case '/'
  A = fliplr(eye(n));
 case '\'
  A = eye(n);
 case '|'
  A(:,1) = 1;
 case '-'
  A(1,:) = 1;
 case '+'
  A(:,1) = 1;
  A(1,:) = 1;
 case 'x'
  A = eye(n) | fliplr(diag(ones(n-1,1),-1));
 case '.'
  A(1:2,1:2)=1;
 otherwise
  error(['Undefined hatch pattern "' hatch '".']);
end


text:
function n = hex_grid_01_n ( nodes_per_layer )

%*****************************************************************************80
%
%% HEX_GRID_01_N computes the number of unit square hex grid points.
%
%  Discussion:
%
%    This routine determines the value of N, the number of
%    hex grid points, from the fundamental hexagonal grid
%    parameter NODES_PER_LAYER.
%
%    A hexagonal grid is defined in the unit square [0,1] x [0,1].
%
%    All nodes of the grid lie on one of LAYERS horizontal lines.
%    The first of these lines is the X axis, and each successive
%    line is HY units higher.
%
%    On all the odd numbered lines, there are NODES_PER_LAYER points,
%    equally spaced from 0 to 1, with a spacing of HX.
%
%    On the even numbered lines, there are NODES_PER_LAYER-1 points,
%    whose values are the midpoints of successive intervals on
%    an odd numbered line.  (The grid is staggered).
%
%    HY = HX * sqrt ( 3 ) / 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODES_PER_LAYER, the number of grid points on the first
%    horizontal layer of points.
%
%    Output, integer N, the number of hex grid points.
%
  if ( nodes_per_layer < 1 )

    n = 0;

  elseif ( nodes_per_layer == 1 )

    n = 1;

  else

    layers = hex_grid_01_layers ( nodes_per_layer );

    n = nodes_per_layer       * floor ( ( layers + 1 ) / 2 ) + ...
      ( nodes_per_layer - 1 ) * floor ( ( layers     ) / 2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function M=convertPF(P,Pp,isProj)
% Canonical case. Converts from P to F
% I should remove that function
%
% Vincent's Structure From Motion Toolbox      Version 1.1
% Copyright (C) 2008-2011 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

F=Pp; Pp=P;

if isProj
  if isempty(F)
    % Reference: HZ2, p246, Table 9.1
    ep=Pp(:,4);
    M = skew(ep)*Pp(:,1:3);
    return
  else
    % Reference: HZ2, p256, Result 9.14
    [U,S,V] = svd(F); %#ok<NASGU>
    ep = U(:,end);
    M = [ skew(ep)*F ep ];
    return
  end
else    % F has form 14.1, HZ2, p345 :  [ 0 0 a; 0 0 b; c d e ]
  if isempty(F)
    % Reference: HZ2, p348, table 14.1
    M=[ 0 0 Pp(2,3); 0 0 -Pp(1,3); Pp(1,3)*Pp(2,1)-Pp(1,1)*Pp(2,3) ...
      Pp(1,3)*Pp(2,2)-Pp(1,2)*Pp(2,3) Pp(1,3)*Pp(2,4)-Pp(2,3)*P(1,4) ];
    return
  else
    % Reference: HZ2, p348, table 14.1
    M=zeros(3,4); M(3,4)=1;
    M(1,3)=-F(2,3);
    M(2,3)=F(1,3);
    M(1,1)=-(F(3,1)/M(1,3)*M(2,3)/M(1,3)-1)/(1+(M(2,3)/M(1,3))^2);
    M(2,1)=(F(3,1)+M(1,1)*M(2,3))/M(1,3);

    M(1,2)=-((F(3,2)/M(1,3)-1)*M(2,3)/M(1,3))/(1+(M(2,3)/M(1,3))^2);
    M(2,2)=(F(3,2)+M(1,2)*M(2,3))/M(1,3);
    M(1,4)=1; M(2,4)=(F(3,3)+M(2,3)*M(1,4))/M(1,3);
    return
  end
end
error('Bad input');

author: vrabaud


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2010, John T. Ramshur, jramshur@gmail.com
%
% This file is part of HRVAS
%
% HRVAS 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.
%
% HRVAS 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 HRVAS.  If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function output = freqDomainHRV(ibi,VLF,LF,HF,AR_order,window, ...
    noverlap,nfft,fs,methods,flagPlot)
%freqDomainHRV - calculates freq domain HRV using FFT, AR, and Lomb-Scargle
%methods
%
%Inputs:    ibi = 2Dim array of time (s) and inter-beat interval (s)
%           AR_order = order of AR model
%           window = # of samples in window
%           noverlap = # of samples to overlap
%           fs = cubic spline interpolation rate / resample rate (Hz)
%           nfft = # of points in the frequency axis
%           methods = cell array of strings that defines the methods used to
%               calculate freqDomain. The default is all to use
%               all three methods.
%               methods={'welch','ar','lomb'}
%           flagPlot = flag to tell function to plot PSD. 1=plot,
%           0=don't plot, default is 0.
%Outputs:   output is a structure containg all HRV. One field for each
%           PSD method.
%           Output units include:
%               peakHF,LF,VLF (Hz)
%               aHF,aLF,aVLF (ms^2)
%               pHF,pLF,pVLF (%)
%               nHF,nLF,nVLF (%)
%               PSD (ms^2/Hz)
%               F (Hz)
%Usage:  (1) To compute freq. domain HRV on a ibi data set named dIBI
%        using VLF=[0.0-0.16], LF =[0.16-0.6], HF=[0.6 3],
%        AR model order = 16, welch window width = 256,
%        # of overlap pnts in welch window (50%) = 128, # of pnts in fft = 512,
%        IBI resample rate = 10Hz
%
%        Use: output = freqDomainHRV(sampledata,[0 .16],[.16 .6],[.6 3], ...
%                       16, 256, 128, 512, 10);
%
%        (2) To do the above and also plot all three power
%        spectrum densities (PSD)
%
%        Use: output = freqDomainHRV(sampledata,[0 .16],[.16 .6],[.6 3], ...
%                       16,256,128,512,10,{'welch','ar','lomb'},1);

    %check input
    if nargin<9
        error('Not enough input arguments!')
    elseif nargin<10
        methods={'welch','ar','lomb'};
        flagPlot=false;
    elseif nargin<11
        flagPlot=false;
    end

    flagWelch=false; flagAR=false; flagLomb=false;
    for m=1:length(methods)
        if strcmpi(methods{m},'welch')
            flagWelch=true;
        elseif strcmpi(methods{m},'ar')
            flagAR=true;
        elseif strcmpi(methods{m},'lomb')
            flagLomb=true;
        end
    end

    t=ibi(:,1); %time (s)
    y=ibi(:,2); %ibi (s)

    y=y.*1000; %convert ibi to ms
    %assumes ibi units are seconds

    maxF=fs/2;

    %prepare y
    y=detrend(y,'linear');
    y=y-mean(y);

    %Welch FFT
    if flagWelch
        [output.welch.psd,output.welch.f] = ...
            calcWelch(t,y,window,noverlap,nfft,fs);
        output.welch.hrv = ...
            calcAreas(output.welch.f,output.welch.psd,VLF,LF,HF);
    else
        output.welch=emptyData(nfft,maxF);
    end

    %AR
    if flagAR
        [output.ar.psd,output.ar.f]=calcAR(t,y,fs,nfft,AR_order);
        output.ar.hrv=calcAreas(output.ar.f,output.ar.psd,VLF,LF,HF);
    else
        output.ar=emptyData(nfft,maxF);
    end

    %Lomb
    if flagLomb
        [output.lomb.psd,output.lomb.f]=calcLomb(t,y,nfft,maxF);
        output.lomb.hrv = ...
            calcAreas(output.lomb.f,output.lomb.psd,VLF,LF,HF,true);
    else
        output.lomb=emptyData(nfft,maxF);
    end

    %plot all three psd
    if flagPlot
    figure;
    h1=subplot(3,1,1);
    plotPSD(h1,output.welch.f,output.welch.psd,VLF,LF,HF,[0 0.6],[]);
    legend('welch')
    h2=subplot(3,1,2);
    plotPSD(h2,output.ar.f,output.ar.psd,VLF,LF,HF,[0 0.6],[]);
    legend('AR')
    h3=subplot(3,1,3);
    plotPSD(h3,output.lomb.f,output.lomb.psd,VLF,LF,HF,[0 0.6],[]);
    legend('Lomb-Scargle')
    end
end

function [PSD,F]=calcWelch(t,y,window,noverlap,nfft,fs)
%calFFT - Calculates the PSD using Welch method.
%
%Inputs:
%Outputs:

    %Prepare y
    t2 = t(1):1/fs:t(length(t));%time values for interp.
    y=interp1(t,y,t2','spline')'; %cubic spline interpolation
    y=y-mean(y); %remove mean

    %Calculate Welch PSD using hamming windowing
    [PSD,F] = pwelch(y,window,noverlap,(nfft*2)-1,fs,'onesided');

end

function [PSD,F]=calcAR(t,y,fs,nfft,AR_order)
%calAR - Calculates the PSD using Auto Regression model.
%
%Inputs:
%Outputs:

    %Prepare y
    t2 = t(1):1/fs:t(length(t)); %time values for interp.
    y=interp1(t,y,t2,'spline')'; %cubic spline interpolation
    y=y-mean(y); %remove mean
    y = y.*hamming(length(y)); %hamming window

    %Calculate PSD
    %Method 1
%     [A, variance] = arburg(y,AR_order); %AR using Burg method
%     [H,F] = freqz(1,A,nfft,fs);
%     PSD=(abs(H).^2).*(variance/fs); %malik, p.67
    %Method 2
    [PSD,F]=pburg(y,AR_order,(nfft*2)-1,fs,'onesided');
    %Method 3
%      h=spectrum.burg;
%      hpsd = psd(h, y, 'NFFT', nfft, 'Fs', 2);
%      F=hpsd.Frequencies;
%      PSD=hpsd.Data;

end

function [PSD,F]=calcLomb(t,y,nfft,maxF)
%calLomb - Calculates the PSD using Lomb-Scargle method.
%
%Inputs:
%Outputs:

    %Calculate PSD
    deltaF=maxF/nfft;
    F = linspace(0.0,maxF-deltaF,nfft);
    PSD=lomb2(y,t,F,false); %calc lomb psd
end

function output=calcAreas(F,PSD,VLF,LF,HF,flagNorm)
%calcAreas - Calulates areas/energy under the PSD curve within the freq
%bands defined by VLF, LF, and HF. Returns areas/energies as ms^2,
%percentage, and normalized units. Also returns LF/HF ratio.
%
%Inputs:
%   PSD: PSD vector
%   F: Freq vector
%   VLF, LF, HF: array containing VLF, LF, and HF freq limits
%   flagNormalize: option to normalize PSD to max(PSD)
%Output:
%
%Usage:
%
%
%   Modified from Gary Clifford's ECG Toolbox: calc_lfhf.m

    if nargin<6
       flagNorm=false;
    end

    %normalize PSD if needed
    if flagNorm
        PSD=PSD/max(PSD);
    end

    % find the indexes corresponding to the VLF, LF, and HF bands
    iVLF= (F>=VLF(1)) & (F<=VLF(2));
    iLF = (F>=LF(1)) & (F<=LF(2));
    iHF = (F>=HF(1)) & (F<=HF(2));

    %Find peaks
      %VLF Peak
      tmpF=F(iVLF);
      tmppsd=PSD(iVLF);
      [pks,ipks] = zipeaks(tmppsd);
      if ~isempty(pks)
        [tmpMax i]=max(pks);
        peakVLF=tmpF(ipks(i));
      else
        [tmpMax i]=max(tmppsd);
        peakVLF=tmpF(i);
      end
      %LF Peak
      tmpF=F(iLF);
      tmppsd=PSD(iLF);
      [pks,ipks] = zipeaks(tmppsd);
      if ~isempty(pks)
        [tmpMax i]=max(pks);
        peakLF=tmpF(ipks(i));
      else
        [tmpMax i]=max(tmppsd);
        peakLF=tmpF(i);
      end
      %HF Peak
      tmpF=F(iHF);
      tmppsd=PSD(iHF);
      [pks,ipks] = zipeaks(tmppsd);
      if ~isempty(pks)
        [tmpMax i]=max(pks);
        peakHF=tmpF(ipks(i));
      else
        [tmpMax i]=max(tmppsd);
        peakHF=tmpF(i);
      end

    % calculate raw areas (power under curve), within the freq bands (ms^2)
    aVLF=trapz(F(iVLF),PSD(iVLF));
    aLF=trapz(F(iLF),PSD(iLF));
    aHF=trapz(F(iHF),PSD(iHF));
    aTotal=aVLF+aLF+aHF;

    %calculate areas relative to the total area (%)
    pVLF=(aVLF/aTotal)*100;
    pLF=(aLF/aTotal)*100;
    pHF=(aHF/aTotal)*100;

    %calculate normalized areas (relative to HF+LF, n.u.)
    nLF=aLF/(aLF+aHF);
    nHF=aHF/(aLF+aHF);

    %calculate LF/HF ratio
    lfhf =aLF/aHF;

    %create output structure
    if flagNorm
        output.aVLF=round(aVLF*1000)/1000;
        output.aLF=round(aLF*1000)/1000;
        output.aHF=round(aHF*1000)/1000;
        output.aTotal=round(aTotal*1000)/1000;
    else
        output.aVLF=round(aVLF*100)/100; % round
        output.aLF=round(aLF*100)/100;
        output.aHF=round(aHF*100)/100;
        output.aTotal=round(aTotal*100)/100;
    end
    output.pVLF=round(pVLF*10)/10;
    output.pLF=round(pLF*10)/10;
    output.pHF=round(pHF*10)/10;
    output.nLF=round(nLF*1000)/1000;
    output.nHF=round(nHF*1000)/1000;
    output.LFHF=round(lfhf*1000)/1000;
    output.peakVLF=round(peakVLF(1)*100)/100;
    output.peakLF=round(peakLF(1)*100)/100;
    output.peakHF=round(peakHF(1)*100)/100;
end

function plotPSD(aH,F,PSD,VLF,LF,HF,limX,limY)

    color.vlf=[.5 .5 1];    %vlf color
    color.lf=[.7 .5 1];     %lf color
    color.hf=[.5 1 1];      %hf color

    % find the indexes corresponding to the VLF, LF, and HF bands
    iVLF= find( (F>=VLF(1)) & (F<VLF(2)) );
    iLF = find( (F>=LF(1)) & (F<LF(2)) );
    iHF = find( (F>=HF(1)) & (F<HF(2)) );

    %plot area under PSD curve
    area(aH,F(:),PSD(:),'FaceColor',[.8 .8 .8]);
    hold(aH);
    area(aH,F(iVLF(1):iVLF(end)+1),PSD(iVLF(1):iVLF(end)+1), ...
        'FaceColor',color.vlf);
    area(aH,F(iLF(1):iLF(end)+1),PSD(iLF(1):iLF(end)+1), ...
        'FaceColor',color.lf);
    area(aH,F(iHF(1):iHF(end)+1),PSD(iHF(1):iHF(end)+1), ...
        'FaceColor',color.hf);

    if ~isempty(limX)
        set(aH,'xlim',limX)
    else
        limX=[min(F) max(F)];
    end
    if ~isempty(limY)
        set(aH,'ylim',limY)
    else
        limY=[min(PSD) max(PSD)];
    end

    %draw vertical lines around freq bands
    line1=line([VLF(2) VLF(2)],[limY(1) limY(2)]);
    set(line1,'color',[1 0 0],'parent',aH);
    line2=line([LF(2) LF(2)],[limY(1) limY(2)]);
    set(line2,'color',[1 0 0],'parent',aH);
    line3=line([HF(2) HF(2)],[limY(1) limY(2)]);
    set(line3,'color',[1 0 0],'parent',aH);

    hold(aH)

end

function output=emptyData(nfft,maxF)
%create output structure of zeros

    output.hrv.aVLF=0;
    output.hrv.aLF=0;
    output.hrv.aHF=0;
    output.hrv.aTotal=0;
    output.hrv.pVLF=0;
    output.hrv.pLF=0;
    output.hrv.pHF=0;
    output.hrv.nLF=0;
    output.hrv.nHF=0;
    output.hrv.LFHF=0;
    output.hrv.peakVLF=0;
    output.hrv.peakLF=0;
    output.hrv.peakHF=0;

    %PSD with all zeros
    deltaF=maxF/nfft;
    output.f = linspace(0.0,maxF-deltaF,nfft);
    output.psd=zeros(length(output.f),1);

end

function [pks locs]=zipeaks(y)
%zippeaks: finds local maxima of input signal y
%Usage:  peak=zipeaks(y);
%Returns 2x(number of maxima) array
%pks = value at maximum
%locs = index value for maximum
%
%Reference:  2009, George Zipfel (Mathworks File Exchange #24797)

%check dimentions
if isempty(y)
    Warning('Empty input array')
    pks=[]; locs=[];
    return
end
[rows cols] = size(y);
if cols==1 && rows>1 %all data in 1st col
    y=y';
elseif cols==1 && rows==1
    Warning('Short input array')
    pks=[]; locs=[];
    return
end

%Find locations of local maxima
%yD=1 at maxima, yD=0 otherwise, end point maxima excluded
    N=length(y)-2;
    yD=[0 (sign(sign(y(2:N+1)-y(3:N+2))-sign(y(1:N)-y(2:N+1))-.1)+1) 0];
%Indices of maxima and corresponding values of y
    Y=logical(yD);
    I=1:length(Y);
    locs=I(Y);
    pks=y(Y);
end

author: jramshur

repo: HRVAS

save_path: github-repos/MATLAB/jramshur-HRVAS


text:
% mlf2pdf_ex is an example for how to use mlf2pdf.
% It creates two buttons (Exit and PDF), which you may use either to exit the
% script or to create the pdf file from the current figure. So you can edit
% the plot "online" and create the pdf from the figure the way it is.

function mlf2pdf_ex

fpl=figure('units','centimeters',...%plotfigure
           'Position',[10 10 10 5]); % position: [llx lly width height]
f=figure('units','centimeters','Position',[5 10 2 2]);     %Buttons
uicontrol(f,'Style','pushbutton','String','-> PDF',...
            'Units','normalized','Pos',[0 0 1 .5],'Callback',@callback_mlf2pdf);
uicontrol(f,'Style','pushbutton','String','Exit',...
            'Units','normalized','Pos',[0 0.5 1 .5],'Callback','close all');
figure(fpl);                %activate plot figure

filename='yourfilename';    %enter your filename here
yourplot;                   %place your plot commands here

 function callback_mlf2pdf(dummy1,dummy2)
    mlf2pdf(fpl,filename,'SIUnits');
 end

 function yourplot
        X=0:pi/10:2*pi;
        Y=sin(X);
        plot(X,Y);
        axis([0,2*pi,-1.5,1.5]);
        grid on;
end
end


text:
%==========================================================================
%
%   Active contour with Chen-Vese Method
%   for image segementation
%
%   Implemented by Yue Wu (yue.wu@tufts.edu)
%   Tufts University
%   Feb 2009
%   http://sites.google.com/site/rexstribeofimageprocessing/
%
%   all rights reserved
%   Last update 02/26/2009
%--------------------------------------------------------------------------
%   Usage of varibles:
%   input:
%       I           = any gray/double/RGB input image
%       mask        = initial mask, either customerlized or built-in
%       num_iter    = total number of iterations
%       mu          = weight of length term
%       method      = submethods pick from ('chen','vector','multiphase')
%
%   Types of built-in mask functions
%       'small'     = create a small circular mask
%       'medium'    = create a medium circular mask
%       'large'     = create a large circular mask
%       'whole'     = create a mask with holes around
%       'whole+small' = create a two layer mask with one layer small
%                       circular mask and the other layer with holes around
%                       (only work for method 'multiphase')
%   Types of methods
%       'chen'      = general CV method
%       'vector'    = CV method for vector image
%       'multiphase'= CV method for multiphase (2 phases applied here)
%
%   output:
%       phi0        = updated level set function
%
%--------------------------------------------------------------------------
%
% Description: This code implements the paper: "Active Contours Without
% Edges" by Chan and Vese for method 'chen', the paper:"Active Contours Without
% Edges for vector image" by Chan and Vese for method 'vector', and the paper
% "A Multiphase Level Set Framework for Image Segmentation Using the
% Mumford and Shah Model" by Chan and Vese.
%
%--------------------------------------------------------------------------
% Deomo: Please see HELP file for details
%==========================================================================

function [seg, phi0]= chenvese(I,mask,num_iter,mu,method)

%%
%-- Default settings
%   length term mu = 0.2 and default method = 'chan'
if(~exist('mu','var'))
    mu=0.2;
end

if(~exist('method','var'))
    method = 'chan';
end

%-- End default settings

%%
%-- Initializations on input image I and mask
%  resize original image
row=size(I,1);
col=size(I,2);
s = 200./min(size(I,1),size(I,2)); % resize scale
if s<1
    I = imresize(I,s);
end

%   auto mask settings
if ischar(mask)
    switch lower (mask)
        case 'small'
            mask = maskcircle2(I,'small');
        case 'medium'
            mask = maskcircle2(I,'medium');
        case 'large'
            mask = maskcircle2(I,'large');
        case 'whole'
            mask = maskcircle2(I,'whole');
            %mask = init_mask(I,30);
        case 'whole+small'
            m1 = maskcircle2(I,'whole');
            m2 = maskcircle2(I,'small');
            mask = zeros(size(I,1),size(I,2),2);
            mask(:,:,1) = m1(:,:,1);
            mask(:,:,2) = m2(:,:,2);
        otherwise
            error('unrecognized mask shape name (MASK).');
    end
else
    if s<1
        mask = imresize(mask,s);
    end
    if size(mask,1)>size(I,1) || size(mask,2)>size(I,2)
        error('dimensions of mask unmathch those of the image.')
    end
    switch lower(method)
        case 'multiphase'
            if  (size(mask,3) == 1)
                error('multiphase requires two masks but only gets one.')
            end
    end

end

switch lower(method)
    case 'chan'
        if size(I,3)== 3
            P = rgb2gray(uint8(I));
            P = double(P);
        elseif size(I,3) == 2
            P = 0.5.*(double(I(:,:,1))+double(I(:,:,2)));
        else
            P = double(I);
        end
        layer = 1;

    case 'vector'
        s = 200./min(size(I,1),size(I,2)); % resize scale
        I = imresize(I,s);
        mask = imresize(mask,s);
        layer = size(I,3);
        if layer == 1
            display('only one image component for vector image')
        end
        P = double(I);

    case 'multiphase'
        layer = size(I,3);
        if size(I,1)*size(I,2)>200^2
            s = 200./min(size(I,1),size(I,2)); % resize scale
            I = imresize(I,s);
            mask = imresize(mask,s);
        end

        P = double(I);  %P store the original image
    otherwise
        error('!invalid method')
end
%-- End Initializations on input image I and mask

%%
%--   Core function
switch lower(method)
    case {'chan','vector'}
        %-- SDF
        %   Get the distance map of the initial mask

        mask = mask(:,:,1);
        phi0 = bwdist(mask)-bwdist(1-mask)+im2double(mask)-.5;
        %   initial force, set to eps to avoid division by zeros
        force = eps;
        %-- End Initialization

        %-- Display settings
        figure();
        subplot(2,2,1); imshow(I); title('Input Image');
        subplot(2,2,2); contour(flipud(phi0), [0 0], 'r','LineWidth',1); title('initial contour');
        subplot(2,2,3); title('Segmentation');
        %-- End Display original image and mask

        %-- Main loop
        for n=1:num_iter
            inidx = find(phi0>=0); % frontground index
            outidx = find(phi0<0); % background index
            force_image = 0; % initial image force for each layer
            for i=1:layer
                L = im2double(P(:,:,i)); % get one image component
                c1 = sum(sum(L.*Heaviside(phi0)))/(length(inidx)+eps); % average inside of Phi0
                c2 = sum(sum(L.*(1-Heaviside(phi0))))/(length(outidx)+eps); % verage outside of Phi0
                force_image=-(L-c1).^2+(L-c2).^2+force_image;
                % sum Image Force on all components (used for vector image)
                % if 'chan' is applied, this loop become one sigle code as a
                % result of layer = 1
            end

            % calculate the external force of the image
            force = mu*kappa(phi0)./max(max(abs(kappa(phi0))))+1/layer.*force_image;

            % normalized the force
            force = force./(max(max(abs(force))));

            % get stepsize dt
            dt=0.5;

            % get parameters for checking whether to stop
            old = phi0;
            phi0 = phi0+dt.*force;
            new = phi0;
            indicator = checkstop(old,new,dt);
%             if n==1
%                 figure;
%                 showphi(I,phi0,n);
%             end
%

            % intermediate output

            if(mod(n,20) == 0)
                showphi(I,phi0,n);
            end;
            if indicator % decide to stop or continue
                showphi(I,phi0,n);

                %make mask from SDF
                seg = phi0<=0; %-- Get mask from levelset
                if s<1
                    seg = imresize(seg,[row col]);
                end

                subplot(2,2,4); imshow(seg); title('Global Region-Based Segmentation');

                return;
            end
        end;
        showphi(I,phi0,n);

        %make mask from SDF
        seg = phi0<=0; %-- Get mask from levelset

        subplot(2,2,4); imshow(seg); title('Global Region-Based Segmentation');
    case 'multiphase'
        %-- Initializations
        %   Get the distance map of the initial masks
        mask1 = mask(:,:,1);
        mask2 = mask(:,:,2);
        phi1=bwdist(mask1)-bwdist(1-mask1)+im2double(mask1)-.5;%Get phi1 from the initial mask 1
        phi2=bwdist(mask2)-bwdist(1-mask2)+im2double(mask2)-.5;%Get phi1 from the initial mask 2

        %-- Display settings
        figure();
        subplot(2,2,1);
        if layer ~= 1
            imshow(I); title('Input Image');
        else
            imagesc(P); axis image; colormap(gray);title('Input Image');
        end
        subplot(2,2,2);
        hold on
        contour(flipud(mask1),[0,0],'r','LineWidth',2.5);
        contour(flipud(mask1),[0,0],'x','LineWidth',1);
        contour(flipud(mask2),[0,0],'g','LineWidth',2.5);
        contour(flipud(mask2),[0,0],'x','LineWidth',1);
        title('initial contour');
        hold off
        subplot(2,2,3); title('Segmentation');
        %-- End display settings

        %Main loop
        for n=1:num_iter
            %-- Narrow band for each phase
            nb1 = find(phi1<1.2 & phi1>=-1.2); %narrow band of phi1
            inidx1 = find(phi1>=0); %phi1 frontground index
            outidx1 = find(phi1<0); %phi1 background index

            nb2 = find(phi2<1.2 & phi2>=-1.2); %narrow band of phi2
            inidx2 = find(phi2>=0); %phi2 frontground index
            outidx2 = find(phi2<0); %phi2 background index
            %-- End initiliazaions on narrow band

            %-- Mean calculations for different partitions
            %c11 = mean (phi1>0 & phi2>0)
            %c12 = mean (phi1>0 & phi2<0)
            %c21 = mean (phi1<0 & phi2>0)
            %c22 = mean (phi1<0 & phi2<0)

            cc11 = intersect(inidx1,inidx2); %index belong to (phi1>0 & phi2>0)
            cc12 = intersect(inidx1,outidx2); %index belong to (phi1>0 & phi2<0)
            cc21 = intersect(outidx1,inidx2); %index belong to (phi1<0 & phi2>0)
            cc22 = intersect(outidx1,outidx2); %index belong to (phi1<0 & phi2<0)

            f_image11 = 0;
            f_image12 = 0;
            f_image21 = 0;
            f_image22 = 0; % initial image force for each layer

            for i=1:layer
                L = im2double(P(:,:,i)); % get one image component

                if isempty(cc11)
                    c11 = eps;
                else
                    c11 = mean(L(cc11));
                end

                if isempty(cc12)
                    c12 = eps;
                else
                    c12 = mean(L(cc12));
                end

                if isempty(cc21)
                    c21 = eps;
                else
                    c21 = mean(L(cc21));
                end

                if isempty(cc22)
                    c22 = eps;
                else
                    c22 = mean(L(cc22));
                end

                %-- End mean calculation

                %-- Force calculation and normalization
                % force on each partition

                f_image11=(L-c11).^2.*Heaviside(phi1).*Heaviside(phi2)+f_image11;
                f_image12=(L-c12).^2.*Heaviside(phi1).*(1-Heaviside(phi2))+f_image12;
                f_image21=(L-c21).^2.*(1-Heaviside(phi1)).*Heaviside(phi2)+f_image21;
                f_image22=(L-c22).^2.*(1-Heaviside(phi1)).*(1-Heaviside(phi2))+f_image22;
            end

            % sum Image Force on all components (used for vector image)
            % if 'chan' is applied, this loop become one sigle code as a
            % result of layer = 1

            % calculate the external force of the image

            % curvature on phi1
            curvature1 = mu*kappa(phi1);
            curvature1 = curvature1(nb1);
            % image force on phi1
            fim1 = 1/layer.*(-f_image11(nb1)+f_image21(nb1)-f_image12(nb1)+f_image22(nb1));
            fim1 = fim1./max(abs(fim1)+eps);

            % curvature on phi2
            curvature2 = mu*kappa(phi2);
            curvature2 = curvature2(nb2);
            % image force on phi2
            fim2 = 1/layer.*(-f_image11(nb2)+f_image12(nb2)-f_image21(nb2)+f_image22(nb2));
            fim2 = fim2./max(abs(fim2)+eps);

            % force on phi1 and phi2
            force1 = curvature1+fim1;
            force2 = curvature2+fim2;
            %-- End force calculation

            % detal t
            dt = 1.5;

            old(:,:,1) = phi1;
            old(:,:,2) = phi2;

            %update of phi1 and phi2
            phi1(nb1) = phi1(nb1)+dt.*force1;
            phi2(nb2) = phi2(nb2)+dt.*force2;

            new(:,:,1) = phi1;
            new(:,:,2) = phi2;

            indicator = checkstop(old,new,dt);

            if indicator
                showphi(I, new, n);
                %make mask from SDF
                seg11 = (phi1>0 & phi2>0); %-- Get mask from levelset
                seg12 = (phi1>0 & phi2<0);
                seg21 = (phi1<0 & phi2>0);
                seg22 = (phi1<0 & phi2<0);

                se = strel('disk',1);
                aa1 = imerode(seg11,se);
                aa2 = imerode(seg12,se);
                aa3 = imerode(seg21,se);
                aa4 = imerode(seg22,se);
                seg = aa1+2*aa2+3*aa3+4*aa4;
                if s<1
                    seg = imresize(seg,[row col]);
                end
                subplot(2,2,4); imagesc(seg);axis image;title('Global Region-Based Segmentation');

                return
            end
            % re-initializations
            phi1 = reinitialization(phi1, 0.6);%sussman(phi1, 0.6);%
            phi2 = reinitialization(phi2, 0.6);%sussman(phi2,0.6);

            %intermediate output
            if(mod(n,20) == 0)
                phi(:,:,1) = phi1;
                phi(:,:,2) = phi2;
                showphi(I, phi, n);
            end;
        end;
        phi(:,:,1) = phi1;
        phi(:,:,2) = phi2;
        showphi(I, phi, n);
        %make mask from SDF
        seg11 = (phi1>0 & phi2>0); %-- Get mask from levelset
        seg12 = (phi1>0 & phi2<0);
        seg21 = (phi1<0 & phi2>0);
        seg22 = (phi1<0 & phi2<0);

        se = strel('disk',1);
        aa1 = imerode(seg11,se);
        aa2 = imerode(seg12,se);
        aa3 = imerode(seg21,se);
        aa4 = imerode(seg22,se);
        seg = aa1+2*aa2+3*aa3+4*aa4;
        %seg = bwlabel(seg);
        subplot(2,2,4); imagesc(seg);axis image;title('Global Region-Based Segmentation');

end
if s<1
    seg = imresize(seg,[row col]);
end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%% Toolbox differencial calculus - A toolbox to handle differential operators
%
% Copyright (c) 2008 Gabriel Peyre
%

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

%%
% First add to the path additional scripts.
path(path, 'toolbox/');
clear options;

%% Vector Fields
% A vector field |v| is an array of dimension |(n,m,2)|, and |v(:,:,1)| is
% the |x| component of the field, while |v(:,:,2)| is the |y| component

%%
% We can create synthetic fields with sources and vortices.

n = 150; % size of the field
v = load_flow('mixed', n);
% normalize it
v = perform_vf_normalization(v);
% display it
options.subsampling = 5;
clf;
subplot(1,2,1);
options.display_streamlines = 0;
plot_vf( v, [], options );
title('Display with arrows'); axis tight;
subplot(1,2,2);
options.display_streamlines = 1;
plot_vf( v, [], options );
axis xy; axis equal;
title('Display with streamlines'); axis tight;

%%
% One can separate a vector field into a sum of incompressible flow and
% irrotional flow.

% create a random flow
n = 100;
options.bound = 'per';
v = randn(n,n,2);
v = perform_vf_normalization( perform_blurring(v,50) );
% compute the Hodge decomposition
[v1,v2] = compute_hodge_decompositon(v,options);
% display
options.subsampling = 4;
options.display_streamlines = 0;
options.normalize_flow = 1;
options.lgd = { 'Original=Irrot+Incomp' 'Irrotational' 'Incompressible'  };
plot_vf( {v v1 v2}, [], options );

%% Tensor Fields
% A tensor field |T| is an array either of size |(n,m,2,2)|, then
% |T(:,:,i,j)| are the component along axis |i| and |j| of the tensor, or
% an array of size |(n,m,3)|, where we have exploited the fact that the
% tensor is symmetric.

%%
% We can compute the structure tensor of an image, which is derived by
% avergaging the tensor product of the gradient with it self.

%%
% First we load an image.

n = 128;
M = load_image('lena',256);
M = rescale(crop(M,n));

%%
% Then we compute the structure tensor, and display it.
% Each tensor is displayed with a small ellipsoid, whose direction matches
% the direction of the tensor, and color (blue to red) matches its size.

% compute the structure tensor, using an averaging of 8
T = compute_structure_tensor(M,1.5,8);
% we cheat a little for the display, we rescale the tensor field
U = perform_tensor_mapping(T,+1);
U(:,:,1) = perform_histogram_equalization(U(:,:,1), 'linear');
T1 = perform_tensor_mapping(U,-1);
% Display the tensor fields. The color is proportional to the size of the tensor.
clf;
options.sub = 5;
plot_tensor_field(T1, M, options);

%% Image Advection and Fluid Dyamics
% You can compute various things from an image along a vector field.

%%
% Load an image and create a random smooth vector field.

n = 256;
M = load_image('lena',n);
options.bound = 'per';
v = perform_blurring(randn(n,n,2), 40, options);
[tmp,v] = compute_hodge_decompositon(v,options);
v = perform_vf_normalization(v);

%%
% One can perform iterative image advection.

dt = .7; niter = 30;
Ma = M; Msvg = {};
for i=1:niter
    Ma = perform_image_advection(Ma, dt*v, options);
    if mod(i,5)==1
        Msvg{end+1} = Ma;
    end
end
% display
clf;
imageplot(Msvg);

%%
% One can solve the Navier Stockes equations for incompressible fluids
% starting from an image.

n = 128;
M = load_image('lena',256);
M = rescale(crop(M,n));
% initial flow
v = perform_blurring(randn(n,n,2), 40, options);
[tmp,v] = compute_hodge_decompositon(v,options);
v = perform_vf_normalization(v);
% options for the PDE solver
dt = .3;
options.viscosity = 2*dt; % diffusion per frame
options.advspeed = 1*dt; % advection per frame
options.viscosity_texture = .3*dt; % diffusion of the texture
options.texture_histo = 'linear'; % fix the contrast
options.display = 0;
options.niter_fluid = 100;
% solve the PDE
[vlist,A] = perform_fluid_dynamics(v,M,options);
% display
sel = round( linspace(1,size(A,3),6) );
B = mat2cell(A(:,:,sel),n,n,ones(6,1));
clf;
imageplot(B);

%%
% One can perform _Line Integral Convolution_ (LIC), see my image processing toolbox for more
% information.

author: gpeyre

repo: matlab-toolboxes


text:
function [U] = spm_dcm_eeg_channelmodes(dipfit,Nm,xY)
% Returns the channel eigenmodes
% FORMAT [U] = spm_dcm_eeg_channelmodes(dipfit,Nm)
% FORMAT [U] = spm_dcm_eeg_channelmodes(dipfit,Nm,xY)
% dipfit  - spatial model specification
% Nm      - number of modes required (upper bound)
% xY      - data structure
% U       - channel eigenmodes
%__________________________________________________________________________
%
% Uses SVD (an eigensolution) to identify the patterns with the greatest
% prior covariance; assuming independent source activity in the specified
% spatial (forward) model.
%
% if xY is specifed a CVA (a generalised eigensolution) will be used to
% find the spatial modes that are best by the spatial model
%
% U is scaled to ensure trace(U'*L*L'*U) = Nm
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_eeg_channelmodes.m 6806 2016-06-08 20:04:30Z karl $

% number of channels and modes
%--------------------------------------------------------------------------
if nargin < 2, Nm = 8; end

% Spatial modes
%--------------------------------------------------------------------------
[pE pC] = spm_L_priors(dipfit);

% evaluate eigenmodes of gain of covariance in sensor space
%--------------------------------------------------------------------------
dGdg  = spm_diff('spm_erp_L',pE,dipfit,1,'nocat');
L     = spm_cat(dGdg(find(spm_vec(pC))));

% eigen-mode reduction
%--------------------------------------------------------------------------
[U S] = spm_svd(L*L',exp(-8));
S     = diag(S);

if nargin > 2

    % reduce lead field to 32 or less d.f.
    %----------------------------------------------------------------------
    n    = min(size(U,2),32);
    L    = U(:,1:n)*diag(sqrt(S(1:n)));

    % response variable
    %----------------------------------------------------------------------
    for i = 1:length(xY.y)
        Y{i} = xY.y{i};
        Y{i} = Y{i}';
    end
    Y     = spm_cat(Y);
    CVA   = spm_cva(Y,L);
    U     = spm_en(CVA.w);
    S     = U'*(L*L')*U;
    S     = diag(S);
end

% eigen-mode reduction
%--------------------------------------------------------------------------
try
    U = U(:,1:Nm);
    S = S(  1:Nm);
end

% re-scale spatial projector
%--------------------------------------------------------------------------
U     = U/sqrt(mean(S));

author: spm


text:
function pdf = quasigeometric_pdf ( x, a, b )

%*****************************************************************************80
%
%% QUASIGEOMETRIC_PDF evaluates the Quasigeometric PDF.
%
%  Discussion:
%
%    PDF(A,B;X) =    A                     if 0  = X;
%               = (1-A) * (1-B) * B^(X-1)  if 1 <= X.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Darren Glass, Philip Lowry,
%    Quasiquasigeometric Distributions and Extra Inning Baseball Games,
%    Mathematics Magazine,
%    Volume 81, Number 2, April 2008, pages 127-137.
%
%    Paul Nahin,
%    Digital Dice: Computational Solutions to Practical Probability Problems,
%    Princeton University Press, 2008,
%    ISBN13: 978-0-691-12698-2,
%    LC: QA273.25.N34.
%
%  Parameters:
%
%    Input, integer X, the independent variable.
%    0 <= X
%
%    Input, real A, the probability of 0 successes.
%    0.0 <= A <= 1.0.
%
%    Input, real B, the depreciation constant.
%    0.0 <= B < 1.0.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x < 0 )

    pdf = 0.0;

  elseif ( x == 0 )

    pdf = a;

  elseif ( b == 0.0 )

    if ( x == 1 )
      pdf = 1.0;
    else
      pdf = 0.0;
    end

  else

    pdf = ( 1.0 - a ) * ( 1.0 - b ) * b^( x - 1 );

  end

  return
end


text:
function serverStarted = quenchCheckServer(host)
%
%     serverStarted = quenchCheckServer(host)
%
% Author: dakers
% Purpose:
%   Test whether the Quench application has been started.
%
% mrmCheckServer

if(~exist('host','var')||isempty(host)) host = 'localhost:4001'; end

stat = mrMesh(host, 999, 'ping');

% If the status is not -1000, the window opened and the server is there.
% So, close the window and return.
if stat ~= -1000
    serverStarted = 1;
    return;
else
    serverStarted = 0;
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function blas1_c_test09 ( )

%*****************************************************************************80
%
%% TEST09 tests CSCAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2006
%
%  Author:
%
%    John Burkardt
%
  n = 6;

  for j = 1 : n
    x(j) = 10 * j + j * i;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST09\n' );
  fprintf ( 1, '  CSCAL multiplies a complex scalar times a vector.\n' );

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

  da = 5.0 + 0.0 * i;
  x = cscal ( n, da, x, 1 );

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

  for j = 1 : n
    x(j) = 10 * j + j * i;
  end

  da = -2.0 + 1.0 * i;
  x = cscal ( 3, da, x, 2 );

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

  return
end


text:
% Return the total number of edges given the adjacency matrix
% INPUTs: adjacency matrix, nxn
% OUTPUTs: m - total number of edges/links
%
% Note: Valid for both directed and undirected, simple or general graph
% Other routines used: selfLoops.m, isSymmetric.m
% GB, last updated Sep 19, 2012

function m = numEdges(adj)

sl=selfLoops(adj); % counting the number of self-loops

if isSymmetric(adj) && sl==0    % undirected simple graph
    m=sum(sum(adj))/2;

elseif isSymmetric(adj) && sl>0
    m=(sum(sum(adj))-sl)/2+sl; % counting the self-loops only once

elseif not(isSymmetric(adj))   % directed graph (not necessarily simple)
    m=sum(sum(adj));

end

author: aeolianine

repo: octave-networks-toolbox


text:
function [phi,dphi]=pt2TetrahedronCoords3D(x,v)
%%PT2TETRAHEDRONCOORDS3D Convert a 3D point to barycentric coordinates in
%      a tetrahedron using the relative areas of subtetrahedra, as
%      discussed at the beginning of [1]. A signed tetrahedral volume area
%      is used so that points outside of the tetrahedron can also be
%      represented (albeit with some negative weights).
%
%INPUTS: x The 3XnumPts set of points to convert into tetrahedral
%          coordinates.
%        v The 3X4 set of vertices of the tetrahedron.
%
%OUTPUTS: phi A 4XnumPts set of the coordinates in the triangle.
%        dphi A 4X3XnumPts set of derivatives of the elements of phi
%             (rows) taken with respect to the elements of x (columns) for
%             each input measurement (3rd dimension).
%
%In [1], the coordinates are defined as the ratio of the volume of a number
%of sub-tetrahedra (with x as a vertex) and the tetrahedron given in v.
%Here, things differ a little bit in that signed tetrahedral volumes are
%used with a very specific ordering. This allows this function to work
%with points inside the tetrahedron (with all positive weights) as well as
%with points outside the tetrahdron (with some negative weights).
%
%EXAMPLE 1:
%In this example, the coordinates of three points are found and the points
%are recreated from the phi values. The residual error of the reverse
%conversion is seen to be on the order of finite precision errors.
% x=sqrt(3)/3;
% h=sqrt(6)/3;
% d=sqrt(3)/6;
% v=[x,  -d,  -d, 0;
%    0, 1/2,-1/2, 0;
%    0,   0,   0, h];
% xPts=[[0;0;0.5],[x;0;0],[1;1;1]];
% phi=pt2TetrahedronCoords3D(xPts,v);
% xBack=barycentricCoords2Pt(phi,v);
% ResidualErr=max(max(abs(xBack-xPts)))
%
%EXAMPLE 2:
%In this example, we verify that the derivative is consistent with
%numerical differentiation. The absolute error is consistent with what one
%might expect due to finite precision limitations.
% v=randn(3,4);
% x=randn(3,1);
% [~,dphi]=pt2TetrahedronCoords3D(x,v);
% f=@(y)pt2TetrahedronCoords3D(y,v);
% dPhiNumDiff=numDiff(x,f,4);
% RelErr=max(max(abs((dphi-dPhiNumDiff)./dPhiNumDiff)))
%
%REFERENCES:
%[1] M. S. Floater, "Generalized barycentric coordinates and applications,"
%    Acta Numerica, vol. 24, pp. 161-214, 1 May 2015.
%
%December 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

a=v(:,1);
b=v(:,2);
c=v(:,3);
d=v(:,4);

numX=size(x,2);
phi=zeros(4,numX);
dphi=zeros(4,3,numX);
for curX=1:numX
    [Va,dVa]=tetrahedronVolume([b,x(:,curX),d,c],false,2);
    [Vb,dVb]=tetrahedronVolume([a,x(:,curX),c,d],false,2);
    [Vc,dVc]=tetrahedronVolume([a,x(:,curX),d,b],false,2);
    [Vd,dVd]=tetrahedronVolume([a,x(:,curX),b,c],false,2);

    %Nominally, V=tetrahedronVolume(v,false); should equal Va+Vb+Vc+Vd.
    %However, to ensure that phi is normalized to more significant digits,
    %we compute V directly from the sum rather than computing it once
    %outside this loop and reusing it.
    V=Va+Vb+Vc+Vd;

    phi(:,curX)=[Va;Vb;Vc;Vd]/V;

    dphi(:,:,curX)=[dVa;
                    dVb;
                    dVc;
                    dVd]/V;
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [dcenter,i0,i1,j0,j1] = dcenter_pair(trx,fly1,fly2)
% modified AR 3/8/2018 add i0,i1,j0,j1 to output
% initialize
dcenter = nan(1,trx(fly1).nframes);

% get start and end frames of overlap
t0 = max(trx(fly1).firstframe,trx(fly2).firstframe);
t1 = min(trx(fly1).endframe,trx(fly2).endframe);
i0 = nan;
i1 = nan;
j0 = nan;
j1 = nan;
% no overlap
if t1 < t0,
  return;
end

% indices for these frames
i0 = t0 + trx(fly1).off;
i1 = t1 + trx(fly1).off;
j0 = t0 + trx(fly2).off;
j1 = t1 + trx(fly2).off;

% centroid distance
dx = trx(fly2).x_mm(j0:j1)-trx(fly1).x_mm(i0:i1);
dy = trx(fly2).y_mm(j0:j1)-trx(fly1).y_mm(i0:i1);
z = sqrt(dx.^2 + dy.^2);
dcenter(i0:i1) = z;


text:
function [ w, seed ] = direction_uniform_nd ( dim_num, seed )

%*****************************************************************************80
%
%% DIRECTION_UNIFORM_ND generates a random direction vector in ND.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the space.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real W(DIM_NUM), a random direction vector,
%    with unit norm.
%
%    Output, integer SEED, a seed for the random number generator.
%

%
%  Get N values from a standard normal distribution.
%
  [ w, seed ] = r8vec_normal_01 ( dim_num, seed );
%
%  Compute the length of the vector.
%
  norm = sqrt ( sum ( w(1:dim_num).^2 ) );
%
%  Normalize the vector.
%
  w(1:dim_num) = w(1:dim_num) / norm;

  return
end


text:
function gid = slguid()
%SLGUID Generates a GUID (Global Unique Identifier)
%
% $ Syntax $
%   - gid = slguid()
%
% $ Description $
%   - gid = slguid() generates a 128-bit GUID, which is represented by
%     a 16 x 1 UINT8 array. If it succeeds, an 1 x 16 array is returned
%     otherwise an error is reported.
%
% $ Remarks $
%   - In current implementation, it is based on win32guid_core mex
%     cpp function, thus only Win32 or compatible platform is supported.
%
% $ History $
%   - Created by Dahua Lin, on Aug 12, 2006
%

if ~strcmpi(computer(), 'PCWIN') && ~strcmpi(computer(), 'PCWIN64')
    error('sltoolbox:notsupportplatform', ...
        'Only win32 or x64 platform is supported for slguid');
end

[b, gid] = win32guid_core();

if ~b
    error('sltoolbox:rterror', ...
        'Fail to generate GUID');
end

author: lmthang


text:
function	Xcov = embed_covariance(X,D,tau)
% Covariance matrix in time delay embedding space
%  Xcov = embed_covariance(X,D,tau)
% --- Input
%  X : input signal : [Xdim x Tsample x Ntrial]
%  D : time embedding dimension  (integer)
%  tau : delay time in sample number (integer)
% --- Output
%  Xcov([i,t1],[j,t2]) = Xtau([i,t1],:) * Xtau([j,t2],:)'
%  Xtau([i,t0],t) = X(i ,t + tau*(D - t0)) : [(Xdim*D) x T x Ntrial]
%  T = Tsample - tau*(D-1)
%
% 2008-5-22 Masa-aki Sato

[N ,Tall, Nr] = size(X);

Tdelay = tau*(D-1);

T = Tall - Tdelay;

% Time delayed embedding index
Tbgn = Tdelay + 1 - (0:tau:Tdelay);
Tend = Tbgn + T - 1;

% X-dim index
Iend = (1:D)*N;
Ibgn = Iend - N + 1;

Xcov = zeros(N*D,N*D);

for i=1:D
  for j=1:i
  	XX = reshape(X(:,Tbgn(i):Tend(i),:), N,T*Nr) ...
  	   * reshape(X(:,Tbgn(j):Tend(j),:), N,T*Nr)';

  	Xcov(Ibgn(i):Iend(i), Ibgn(j):Iend(j)) = XX;
  	Xcov(Ibgn(j):Iend(j), Ibgn(i):Iend(i)) = XX';
  end
end

Xcov = (Xcov + Xcov')/2;

author: KamitaniLab

repo: GenericObjectDecoding


text:
function test_bug950

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_megrealign test_bug950

% the issue explored here is a reputed crash in megrealign due to a problem
% in the channelposition function.
% in addition balancing is not properly taken into account when creating
% the headmodel for the inverse/forward steps

% load in some data
load(dccnpath('/home/common/matlab/fieldtrip/data/test/latest/raw/meg/preproc_ctf151.mat'));

cfg = [];
cfg.gradient = 'G3BR';
data = ft_denoise_synthetic(cfg, data);

template = data.grad;
template.chanpos(:,3) = template.chanpos(:,3)+1;
template.coilpos(:,3) = template.coilpos(:,3)+1;

cfg = [];
cfg.template{1}    = template;
cfg.inwardshift    = 1;
cfg.headmodel.o    = [0 0 4];
cfg.headmodel.r    = 8;
cfg.headmodel.unit = 'cm';
data2 = ft_megrealign(cfg, data);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


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

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

% Last Modified by GUIDE v2.5 06-Mar-2020 17:08:32

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

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

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

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

    centerUIFigure(hObject);

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

    handles = populateGUI(optionsObj, handles);

    % Update handles structure
    guidata(hObject, handles);

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

function handles = populateGUI(options, handles)
    setOptsDoubleValueStrInGUI(handles, options, 'integratorStepSize', 'integratorOutputStepText');

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

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

        setOptsDoubleValueInObject(handles, options, 'integratorStepSize', 'integratorOutputStepText');

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

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

    errMsg = validateDoubleValue(handles, errMsg, 'integratorOutputStepText', 'Integrator Output Step', 1E-6, Inf, false);

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

    if(isempty(errMsg))
        uiresume(handles.lvd_builtInIntegratorOptionsGUI);
    else
        msgbox(errMsg,'Invalid Integrator Options Inputs','error');
    end

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

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

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

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

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% load an image from uiuc texture database
all_images = retrieve_uiuc(1,1);
x = all_images{1}{1};

% compute its 2d scattering
tic;
options.feat = 'scatt2d_uiuc';
[feat_f,outmeta,outprecomp] = gfeat(x,options);
toc; % about 20 seconds on a core i7 2.4Ghz

% check number of coefficients :
% there should be
% 9 (scales) * 8 (orientations) coefficients of order 1
% and
% (7+5+3+1) (scales) * 8^2 (orientations) coefficients of order 2
theoritical_size = 9*8 + (7+5+3+1)*8^2;

assert(numel(feat_f) == theoritical_size);

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [idx, row, col] = extidx2original(idx, original_size, extended_size)
%EXTIDX2ORIGINAL Convert from extended indexes to original indexes
    [row, col] = ind2sub(extended_size,idx);

    offright = col > (extended_size(2) + original_size(2)) / 2;
    offleft = col < (extended_size(2) - original_size(2)) / 2 + 1;
    notoff = ~(offleft | offright);
    col(offright)=col(offright) - (extended_size(2) + original_size(2)) / 2;
    col(offleft) = col(offleft) + original_size(2) - (extended_size(2) - original_size(2)) / 2;
    col(notoff)=col(notoff) - (extended_size(2) - original_size(2)) / 2;

    idx = sub2ind(original_size,row,col);
end

author: jfaghm

repo: OceanEddies

save_path: github-repos/MATLAB/jfaghm-OceanEddies


text:
function pass = test_univariate( pref )
% Rank-1 PDEs can be solved by univariate spectral methods.  Do the results
% match.
% Alex Townsend, March 2013.

if ( nargin < 1 )
    pref = chebfunpref();
end
tol = 1000*pref.cheb2Prefs.chebfun2eps;

% Simple example in y-variable.
N = chebop2(@(u) diff(u,2,1) + u); N.dbc = 1; N.ubc = 1; u = N \ 0;
L = chebop(@(u) diff(u,2) + u); L.lbc = 1; L.rbc = 1; v = L \ 0;

pass(1) = (length(u) == 1);
pass(2) = (norm(u(0,:) - v) < tol);

% Simple example in x-variable.
N = chebop2(@(u) diff(u,2,2) + u); N.lbc = 1; N.rbc = 1; u = N \ 0;
L = chebop(@(u) diff(u,2) + u); L.lbc = 1; L.rbc = 1; v = L \ 0;

pass(3) = (length(u) == 1);
pass(4) = (norm(u(:,0) - v.') < tol);

%
% % Simple example in x-variable.
% N = chebop2(@(u) diff(diff(u,1,1),1,2) + diff(u,1,1));
% N.lbc = 1; N.dbc = 1; u = N \ 0;
%
% L = chebop(@(u) diff(u,2) + u); L.lbc = 1; L.rbc = 1; v = L \ 0;
%
% pass(j) = (length(u) == 1);  j = j+1;
% pass(j) = (norm(u(:,0) - v) < tol); j = j+1;

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ mK ] = CreateConvMtx2D( mH, numRows, numCols, convShape )
% ----------------------------------------------------------------------------------------------- %
% [ mK ] = CreateConvMtx2D( mH, numRows, numCols, convShape )
% Generates a Convolution Matrix for the 2D Kernel (The Matrix mH) with
% support for different convolution shapes (Full / Same / Valid).
% Input:
%   - mH                -   Input 2D Convolution Kernel.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - numRows           -   Number of Rows.
%                           Number of rows of the image to be convolved.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, ...}.
%   - numCols           -   Number of Columns.
%                           Number of columns of the image to be convolved.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, ...}.
%   - convShape         -   Convolution Shape.
%                           The shape of the convolution which the output
%                           convolution matrix should represent. The
%                           options should match MATLAB's conv2() function
%                           - Full / Same / Valid.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3}.
% Output:
%   - mK                -   Convolution Matrix.
%                           The output convolution matrix. The product of
%                           the matrix 'mK' and and image 'mI' in its
%                           column stack form ('mK * mI(:)') is equivalent
%                           to the convolution of 'mI' with the kernel 'mH'
%                           using the corresponding convolution shape
%                           ('conv2(mI, mH, convShapeString)').
%                           Structure: Matrix (Sparse).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References:
%   1.  MATLAB's 'convmtx2()' - https://www.mathworks.com/help/images/ref/convmtx2.html.
% Remarks:
%   1.  In caes the same convolution is applied on many images, stacking
%       them into a matrix (Each image as columns stacked vector) and
%       applying convolution on each column by matrix multiplication might
%       be more efficient than applying classic convolution per image.
%   2.  The output matrux has the form of doubly block toeplitz matrix. The
%       convolution shape sets where the diagonal of the first column to
%       appear.
% TODO:
%   1.
%   Release Notes:
%   -   1.0.001     22/01/2018  Royi Avital
%       *   Fixed issue with the creatioon of the sparse diagonal matrix.
%           The vector to initialize the diagonal was much bigger than
%           needed. Improved performance in the Unit Test by factor of 10.
%   -   1.0.000     16/01/2018  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

CONVOLUTION_SHAPE_FULL  = 1;
CONVOLUTION_SHAPE_SAME  = 2;
CONVOLUTION_SHAPE_VALID = 3;

numColsKernel   = size(mH, 2);
numBlockMtx     = numColsKernel;

cBlockMtx = cell(numBlockMtx, 1);

for ii = 1:numBlockMtx
    cBlockMtx{ii} = CreateConvMtx1D(mH(:, ii), numRows, convShape);
end

switch(convShape)
    case(CONVOLUTION_SHAPE_FULL)
        % For convolution shape - 'full' the Doubly Block Toeplitz Matrix
        % has the first column as its main diagonal.
        diagIdx     = 0;
        numRowsKron = numCols + numColsKernel - 1;
    case(CONVOLUTION_SHAPE_SAME)
        % For convolution shape - 'same' the Doubly Block Toeplitz Matrix
        % has the first column shifted by the kernel horizontal radius.
        diagIdx     = floor(numColsKernel / 2);
        numRowsKron = numCols;
    case(CONVOLUTION_SHAPE_VALID)
        % For convolution shape - 'valid' the Doubly Block Toeplitz Matrix
        % has the first column shifted by the kernel horizontal length.
        diagIdx     = numColsKernel - 1;
        numRowsKron = numCols - numColsKernel + 1;
end

vI = ones(min(numRowsKron, numCols), 1);
mK = kron(spdiags(vI, diagIdx, numRowsKron, numCols), cBlockMtx{1});
for ii = 2:numBlockMtx
    diagIdx = diagIdx - 1;
    mK = mK + kron(spdiags(vI, diagIdx, numRowsKron, numCols), cBlockMtx{ii});
end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
% test_all_fbp

% todo:
% cuboid_im test
% cuboid_proj test

list = {
'cbct_back test'
'ct_geom test'
'image_geom test'
'sino_geom test'
'cylinder_proj test'
'df_example1'
'ellipse_im test'
'ellipse_sino test'
'ellipsoid_proj test'
'ellipsoid_im test'
'fbp_fan_arc_example'
'fbp_fan_arc_point'
'fbp_fan_flat_example'
'fbp_ramp test'
'fbp2_sino_filter test'
'fbp2_example'
'feldkamp_example'
'jaszczak1 test'
'ir_radon_zwart_powell test'
'rebin_helix test' % helix_example
'rect_im test'
'rect_sino test'
%'sphere_proj test'
};

im nan-fail
run_mfile_local(list)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
classdef gnonomicProjection < sphericalProjection
  % gnonomic projection

  methods

     function proj = gnonomicProjection(varargin)
      proj = proj@sphericalProjection(varargin{:});
    end

    function [x,y] = project(sP,v,varargin)
      % compute polar angles

      [rho,theta] = project@sphericalProjection(sP,v,varargin{:});

      % map to upper hemisphere
      ind = find(theta > pi/2+10^(-10));
      theta(ind)  = pi - theta(ind);

      % turn around antipodal vectors
      sP.sR.antipodal = false; v.antipodal = false;
      ind = ~sP.sR.checkInside(v);
      rho(ind) = rho(ind) + pi;

      % formula for stereographic projection
      r =  tan(theta);

      % compute coordinates
      x = reshape(cos(rho) .* r,size(v));
      y = reshape(sin(rho) .* r,size(v));

    end

    function v = iproject(sP,x,y)
      rho = atan2(y,x);
      theta = atan(sqrt(x.^2 + y.^2));
      v = vector3d('theta',theta,'rho',rho);
    end

  end

end

author: mtex-toolbox