text:
function fList = GetRecon(rawDir, fList, rotFlag, shifts, littleEndian)
% Reads recon data from Glover's recon files
%
%  fList = GetRecon(rawDir, fList, rotFlag, shifts, littleEndian)
%
% Output data are saved in Inplane directory under TSeries and optional
% complexTSeries subdirectories.
%
% rawDir:  Directory with the data
% fList:   0 means you need to find them all.  This should probably be
% empty, not 0
%
% rotFlag: Optional 90-degree rotation can be obtained by setting
%
% shifts:  Optional functional-image shift can also be obtained (to
%    deal with fallback errors or otherwise) by specifying a
%    2D shifts vector. Format is [yShift, xShift] in mm. Positive
%    shifts are to the right and upward.
%
% littleEndian:  Byte format
%
% Returns the list of selected Pfiles.  This may be a subset of all the
% Pfiles.
%
% DBR 8/00
%
% djh & dbr, 7/18/02, Fixed major bug in shifting. It was shifting
% only some of the images because of confusion about junkframes.
% kgs 08/04: added rotFlag
% ras 12/04: added scanList
% ras 03/05: added littleEndian (for scanner upgrade), got
% rid of obsolete compFlag

global mrSESSION HOMEDIR

if isempty(mrSESSION),    load mrSESSION;                       end
if isempty(HOMEDIR),     HOMEDIR = pwd;                         end
if notDefined('rawDir'), rawDir = fullfile(HOMEDIR, 'Raw');     end
if notDefined('littleEndian'),    littleEndian = 1;             end
if notDefined('rotFlag'),     rotFlag = 0;                      end
if notDefined('shifts')
    shifts = [0, 0];
    shiftFlag = 0;
else
    shiftFlag = 1;
end

rawDir = fullfile(rawDir,'Pfiles');

% User selects Pfiles:
if notDefined('fList') || fList==0
    pFileList = SelectPfiles(mrSESSION.functionals);
    fList = find(pFileList);
end

% Set up directories:
parentDir = fullfile(HOMEDIR, 'Inplane', 'Original');
tSeriesDir = fullfile(parentDir, 'TSeries');

if ~exist(tSeriesDir,'dir')
    fprintf('\ntSeries directory %s does not exist.\n',tSeriesDir);
    mkdir(tSeriesDir);
end

if fList==0, fList = 1:length(mrSESSION.functionals); end

% determine byte format (Endian Flag)
% (this depends on when the scan was run: scans run before 3/2005 at the
% Stanford Lucas Center use big-endian format; newer ones use little-endian
% format)
if littleEndian==1,     endianFlag = 'ieee-le';
else                   endianFlag = 'ieee-be';
end

nScans = length(fList);
for iScan=1:nScans
    % Scan loop:
    scan = fList(iScan);
    scanParams = mrSESSION.functionals(scan);

    % if no crop, the first row may accidentally be set to 0
    if isequal(scanParams.crop(1,:),[0 0])
        scanParams.crop(1,:) = [1 1];
    end

    pixelShifts = round(shifts./scanParams.voxelSize(1:2));
    % check whether scanParams.PfileName is absolute path
    if isempty(fileparts(scanParams.PfileName))
        fName = fullfile(rawDir, scanParams.PfileName);
    else
        fName = scanParams.PfileName;
    end

    disp(['Scan ', int2str(scan), ' [', fName, ']'])
    % sName = ['Scan', int2str(scan)];
    % We assign the output name based on how many we recon, not how many
    % possible Pfiles we could have reconned.
    sName = ['Scan', int2str(iScan)];
    MakeDir(sName,tSeriesDir);
    scanDir = fullfile(tSeriesDir, sName);

%     if isfield(scanParams,'cropSize')
%         nRows = scanParams.cropSize(1);
%         nCols = scanParams.cropSize(2);
%     else
%         nRows = scanParams.fullSize(1);
%         nCols = scanParams.fullSize(2);
%     end

    for iSlice=1:length(scanParams.slices)
        % Slice loop:
        disp(['   slice ', int2str(scanParams.slices(iSlice))]);

        % Load fullsize images
        tSeries = LoadRecon(scanParams,fName,iSlice,rotFlag,endianFlag);

        % Remove junk frames
        f0 = scanParams.junkFirstFrames+1;
        nFrames = scanParams.nFrames;
        fEnd = f0 + nFrames - 1;
        tSeries = tSeries(f0:fEnd, :, :);

        % Shift
        if shiftFlag
            for iF=1:nFrames
                tSeries(iF, :, :) = shift(squeeze(tSeries(iF, :, :)), pixelShifts);
            end
        end

        % Crop
        x0 = scanParams.crop(1, 1);
        xN = scanParams.crop(2, 1);
        y0 = scanParams.crop(1, 2);
        yN = scanParams.crop(2, 2);
        tSeries = tSeries(:, y0:yN, x0:xN);

        % Reshape to standard t-series shape
        tSeries = reshape(tSeries, nFrames, (yN-y0+1)*(xN-x0+1)); %#ok<NASGU>

        % Save tSeries file
        tsName = ['tSeries', int2str(scanParams.slices(iSlice))];
        tSeriesFile = fullfile(scanDir, tsName);
        save(tSeriesFile,'tSeries');
    end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [tgrad,fgrad,c]=gabphasegrad(method,varargin)
%GABPHASEGRAD   Phase gradient of the DGT
%   Usage:  [tgrad,fgrad,c] = gabphasegrad('dgt',f,g,a,M);
%           [tgrad,fgrad]   = gabphasegrad('phase',cphase,a);
%           [tgrad,fgrad]   = gabphasegrad('abs',s,g,a);
%
%   `[tgrad,fgrad]=gabphasegrad(method,...)` computes the relative
%   time-frequency gradient of the phase of the |dgt| of a signal.
%   The derivative in time *tgrad* is the relative instantaneous
%   frequency while the frequency derivative *fgrad* is the negative
%   of the local group delay.
%
%   *tgrad* is a measure the deviation from the current channel frequency,
%   so a value of zero means that the instantaneous frequency is equal to
%   the center frequency of the considered channel, a positive value means
%   the true absolute intantaneous frequency is higher than the current
%   channel frequency and vice versa.
%   Similarly, *fgrad* is a measure of deviation from the current time
%   positions.
%
%   *fgrad* is scaled such that distances are measured in samples. Similarly,
%   *tgrad* is scaled such that the Nyquist frequency (the highest possible
%   frequency) corresponds to a value of L/2. The absolute time and
%   frequency positions can be obtained as
%
%      tgradabs = bsxfun(@plus,tgrad,fftindex(M)*L/M);
%      fgradabs = bsxfun(@plus,fgrad,(0:L/a-1)*a);
%
%   Please note that neither *tgrad* and *fgrad* nor *tgradabs* and
%   *fgradabs* are true derivatives of the |dgt| phase. To obtain the true
%   phase derivatives, one has to explicitly pass either 'freqinv' or
%   'timeinv' flags and scale both *tgrad* and *fgrad* by 2*pi/L.
%
%   The computation of *tgrad* and *fgrad* is inaccurate when the absolute
%   value of the Gabor coefficients is low. This is due to the fact the the
%   phase of complex numbers close to the machine precision is almost
%   random. Therefore, *tgrad* and *fgrad* may attain very large random values
%   when `abs(c)` is close to zero.
%
%   The computation can be done using four different methods.
%
%     'dgt'    Directly from the signal using algorithm by Auger and
%              Flandrin.
%
%     'phase'  From the phase of a DGT of the signal. This is the
%              classic method used in the phase vocoder.
%
%     'abs'    From the absolute value of the DGT. Currently this
%              method works only for Gaussian windows.
%
%     'cross'  Directly from the signal using algorithm by Nelson.
%
%   `[tgrad,fgrad]=gabphasegrad('dgt',f,g,a,M)` computes the time-frequency
%   gradient using a DGT of the signal *f*. The DGT is computed using the
%   window *g* on the lattice specified by the time shift *a* and the number
%   of channels *M*. The algorithm used to perform this calculation computes
%   several DGTs, and therefore this routine takes the exact same input
%   parameters as |dgt|.
%
%   The window *g* may be specified as in |dgt|. If the window used is
%   'gauss', the computation will be done by a faster algorithm.
%
%   `[tgrad,fgrad,c]=gabphasegrad('dgt',f,g,a,M)` additionally returns the
%   Gabor coefficients *c*, as they are always computed as a byproduct of the
%   algorithm.
%
%   `[tgrad,fgrad]=gabphasegrad('cross',f,g,a,M)` does the same as above
%   but this time using algorithm by Nelson which is based on computing
%   several DGTs.
%
%   `[tgrad,fgrad]=gabphasegrad('phase',cphase,a)` computes the phase
%   gradient from the phase *cphase* of a DGT of the signal. The original DGT
%   from which the phase is obtained must have been computed using a
%   time-shift of *a* using the default phase convention (`'freqinv'`) e.g.::
%
%        [tgrad,fgrad]=gabphasegrad('phase',angle(dgt(f,g,a,M)),a)
%
%   `[tgrad,fgrad]=gabphasegrad('abs',s,g,a)` computes the phase gradient
%   from the spectrogram *s*. The spectrogram must have been computed using
%   the window *g* and time-shift *a* e.g.::
%
%        [tgrad,fgrad]=gabphasegrad('abs',abs(dgt(f,g,a,M)),g,a)
%
%   `[tgrad,fgrad]=gabphasegrad('abs',s,g,a,difforder)` uses a centered finite
%   diffence scheme of order *difforder* to perform the needed numerical
%   differentiation. Default is to use a 4th order scheme.
%
%   Currently the `'abs'` method only works if the window *g* is a Gaussian
%   window specified as a string or cell array.
%
%   See also: resgram, gabreassign, dgt
%
%   References: aufl95 cmdaaufl97 fl65 ltfatnote042

% AUTHOR: Peter L. Søndergaard, 2008; Zdenek Průša 2015

%narginchk(4,6);

% If no phaseconv flag was passed, add 'relative'
definput = arg_gabphasederivconv;
if ~any(cellfun(@(el) any(strcmpi(el,varargin)),definput.flags.phaseconv))
    varargin{end+1} = 'relative';
end

if nargout<3
    phased = gabphasederiv({'t','f'},method,varargin{:});
else
    [phased,c] = gabphasederiv({'t','f'},method,varargin{:});
end

[tgrad,fgrad] = deal(phased{:});

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [status, errmsg] = ensureDirExists(pth)
% Make sure a directory exists, creating parent directories if needed
%
%  [status, errmsg] = ensureDirExists(pth);
%
% Status = 1 if the path exists or is successfully created, 0 otherwise.
%
% ras 06/06.

if ~ischar(pth), help(mfilename); error('Need a path string'); end

status = 0;
errmsg = '';

if exist(fullpath(pth),'dir') & ~isequal(pth(2),'\')
    % Person sent in a full path that exists.  We are done.
    status = 1;
    return
end

% The person sent in a path that doesn't fully exist.  So we are going to
% live dangerously and create the path for them.  Part of the path might
% exist, and we follow that down, until we start creating the new elements
% of the path.
%

% find all the directory strings
clear subdirs
if isequal(pth(2),'\')
    % Trap the case when we have a mount: \\white.stanford.edu, for example
    tmpPath = pth(3:end);
    tmpPath = ['C:\',tmpPath];  % Cheap hack to make explode work
    tmpSubDirs = explode(filesep,tmpPath);
    subdirs = cell(1,length(tmpSubDirs) - 1);
    for ii=2:length(tmpSubDirs)
        subdirs{ii-1} = tmpSubDirs{ii};
    end
    % Append the \\ in front of the computer name
    tmp = subdirs{1}; tmp = ['\\',tmp]; subdirs{1} = tmp;
else
    subdirs = explode(filesep, fullpath(pth));
end

% step through each subdir, progressively checking
% child directories within parent dirs
if isunix,  par = [filesep subdirs{1}];
else        par = subdirs{1};
end

for ii = 2:length(subdirs)
    if ~exist(fullfile(par,subdirs{ii}),'dir')
        try
            [status, errmsg] = mkdir(par,subdirs{ii});
        catch
            fprintf('ensureDir Exists: Creation of intermediate dir failed: ')
            fprintf('%s\n',errmsg)
            return
        end
        fprintf('Created directory %s%s%s\n',par,filesep,subdirs{ii});
    end

    par = fullfile(par,subdirs{ii});
end

status = 1;

return


text:
function [ ALL_SCENE ] = packupScene( ALL_HYPS, CAN_POOL )
%PACKUPSCENE Summary of this function goes here
%   Detailed explanation goes here
global config

[H,W] = size(ALL_HYPS);
ALL_SCENE = cell(H,W);
all_objects = CAN_POOL.sel_hyps.objects;
room_obj.out_points_w = CAN_POOL.room.out_points_w;
room_obj.name = CAN_POOL.room.name;
room_obj.type = CAN_POOL.room.type;
room_obj.points = CAN_POOL.room.points;
room_obj.x_w = CAN_POOL.room.x_w;
room_obj = addBoundingBox(room_obj,'both');
room_obj.objtype = config.roomtypeid; % change this to be room id
for i = 1:H*W
    local = room_obj;
    sel_hyps = ALL_HYPS(i).sel_hyps;
    for j = 1:length(sel_hyps)
        if sel_hyps(j).fixed
            newObj = all_objects(sel_hyps(j).selID);
%             name = get_object_type_livingroom(j); % change this to be function
%             name = get_object_type(j); % change this to be function
            name = feval(config.typefunc, j);
            for k = 1:length(newObj)
                newObj(k).name = name{1};
                newObj(k).objtype = j;
            end
            local = [local; newObj];
        end
    end
    ALL_SCENE{i} = local;
end

end

author: zouchuhang

repo: LayoutNet


text:
 function st = newfft(om_in, Nd_in, varargin)
%function st = newfft(om, Nd, [options])
%|
%| New version of NUFFT (pun intended) that uses real interpolation kernels.
%| (The original NUFFT code used complex interpolation needlessly.)
%|
%| This returns a "strum" object with methods for both forward and adjoint
%| d-dimensional NUFFT operations.
%| The forward operation is:
%| X(om_m) = \sum_{n=0}^{N-1} x[n] exp(-1i * om_m * n) for m=1,...,M
%| The adjoint operation is:
%| x_adj[n] = \sum_{m=1}^M X(om_m) exp(+1i * om_m * n) for n=0,...,N-1
%| Note that the adjoint is not the "inverse" NUFFT in general.
%|
%| This routine has numerous options for investigative purposes, but is
%| designed so that the default options should be very good choices.
%| Providing the frequencies and the image size and should suffice.
%| Reducing the neighborhood size 'Jd' (default 6) and/or reducing the
%| over-sampled DFT size 'Kd' (default 2*Nd) may be useful for acceleration.
%|
%| in
%|	om [M,d]	"digital" frequencies in radians (can be empty!)
%|			(if empty, then user must provide 'om' to methods)
%|	Nd [d]		image dimensions (N1,N2,...,Nd)
%|
%| options
%|	'Jd' [d]	# of neighbors used (in each direction). (default: 6)
%|	'Kd' [d]	FFT sizes (should be >= Nd). (default: 2*Nd)
%|	n_shift [d]	n = 0-n_shift to N-1-n_shift (default: 0)
%|		Like fft(), the NUFFT expects the signals to be x(0,0), ...
%|		Use n_shift = [N1/2, N2/2, ...] for x(-N1/2,-N2/2,...), ...
%|
%|	'mode'	char	how to compute the NUFFT (default: 'table1')
%|			'exact'	slow FT
%|			'table0' table with nearest neighbor interpolation
%|			'table1' table with linear interpolation (default)
%|			'gg' gaussian factorization of Greengard & Lee 04 (todo)
%|			'sparse' precompute large sparse matrix
%|				caution: this option may require lots of memory!
%|				(this option used internally to make tables too)
%|			recommended: table0 or table1 to save memory.
%|
%|	'oversample' int table oversampling factor (default: 2^9 for table1)
%|	'gram' 0|1	if 1, precompute additional terms needed for gram matrix
%|	'printmem' 0|1	if 1, print memory usage
%|	'phasing' char	'real' : new real kernels (strongly recommended default)
%|			'complex' : original nufft.m complex kernels
%|			'none' : no phase (straw man; table uses it internally)
%|			'flipreal' : real kernels with sign flips (unsupported)
%|	'dotime' char	report cpu time? (default: '')
%|
%|	'ktype'	char	type of interpolation kernel (default: 'minmax:kb')
%|
%| ktype options:
%|	'diric'		Dirichlet interpolator (exact only if Jd = Kd)
%|	'linear'	linear interpolator (a terrible straw man)
%|	'minmax:kb'	minmax interpolator with excellent KB scaling!
%|	'minmax:tuned'	minmax interpolator, somewhat numerically tuned scaling
%|	'minmax:unif'	minmax with uniform scaling factors (not recommended)
%|	'minmax:user'	minmax interpolator with user parameters required:
%|				'alpha', {alpha}, 'beta', {beta}
%|	'kb:minmax'	kaiser-bessel (KB) interpolator (minmax best alpha, m)
%|	'kb:beatty'	KB with parameters from Beatty et al T-MI Jun 2005
%|	'kb:user'	KB with user-specified KB parameters required:
%|				'kb_m', [m] 'kb_alf', [alpha]
%|	@kernel		user-provided inline interpolation kernel(k,J)
%|			(or a cell array of kernels, one for each dimension)
%|			example ..., 'table', 2^11, 'minmax:kb'
%|			where 2^11 is the table over-sampling factor.
%|
%| out
%|	st	strum object with several methods including the following:
%|		st.fft(x, [om])
%|		st.adj(Xo, [om])
%|		st.p		[M, *Kd]	sparse interpolation matrix
%|						(or empty if table-based)
%|		st.sn		[(Nd)]		scaling factors
%|		st.Nd,Jd,Kd,om	copies of inputs
%|
%| *Nd is shorthand for prod(Nd).
%| (Nd) is shorthand for (N1,N2,...,Nd)
%|
%| Copyright 2007-6-3, Jeff Fessler, The University of Michigan

if nargin == 1 && streq(om_in, 'test'), newfft_test, return, end
if nargin < 2, ir_usage, end

cpu etic

% inputs
st.om = om_in; % [M,d] frequency samples
st.Nd = Nd_in; % [d] signal dimentions
st.dd = length(st.Nd); % dimensionality of input space (usually 2 or 3)

% defaults
st.gram = false;
st.Jd = 6 * ones(1, st.dd);
st.Kd = 2 * st.Nd;
st.n_shift = zeros(1, st.dd);
st.mode = 'table1';
st.oversample = []; % aka Ld for table mode
st.order = []; % for table mode
st.ktype = 'kb:minmax';
st.kb_m = []; % [dd] KB parameters
st.kb_alf = [];
st.alpha = {}; % [dd] minmax parameters
st.beta = {};
st.tol = 0;
st.printmem = false;
st.is_kaiser_scale = false;
st.phasing = 'real'; % use new real table by default
st.dotime = '';

% options
st = vararg_pair(st, varargin);
dotime = st.dotime; st = rmfield(st, 'dotime');

st.phase_before = []; % place holders
st.phase_after = [];
st.flips = [];

% special cases of input sampling pattern
if ischar(st.om)
	st.om = nufft_samples(st.om, st.Nd);
end

% checks
if st.dd ~= length(st.Jd) || st.dd ~= length(st.Kd)
	error 'inconsistent dim'
end
if st.dd ~= length(st.n_shift)
	fail('n_shift needs %d columns', st.dd)
end

if ~isempty(st.om) && st.dd ~= size(st.om,2)
	fail('omega needs %d columns', st.dd)
end

if st.gram, fail 'todo: gram not done', end

if any(st.Kd < st.Nd), warning 'Kd < Nd unlikely to work.  Try Kd=2*Nd', end

%
% "midpoint" of scaling factors
%
switch st.phasing
case 'real'
	st.Nmid = floor(st.Nd / 2); % new
otherwise
	st.Nmid = (st.Nd - 1) / 2; % old
end

%
% different interpolation modes
%
switch st.mode

case 'exact' % exact interpolation for testing
	st = strum(st, { ...
		'fft', @newfft_exact_for, '(x, [om])';
		'adj', @newfft_exact_adj, '(X, [om])';
		'p', @newfft_exact_p, '([om])';
		'sn', @(st) ones(st.Nd), '()';
		});

case {'table0', 'table1'} % precomuted interpolator table
	st = newfft_table_init(st);

	% set up phase corrections not included in table initialization
	if streq(st.phasing, 'real') || streq(st.phasing, 'flipreal')
		st.phase_before = newfft_phase_before(st.Kd, st.Nmid);
		st.phase_after = @(om) newfft_phase_after(om, st.Nmid, st.n_shift);
		if ~isempty(st.om) % phase that goes after interpolation:
			st.phase_after = st.phase_after(st.om);
		end
	end

case 'gg' % greengard's gaussian
	error 'todo: gg'
	st = strum(st, { ...
		'fft', @newfft_gg_for, '(x, [om])';
		'adj', @newfft_gg_adj, '(X, [om])';
		'p', @newfft_gg_p, '([om])';
		'sn', @newfft_gg_sn, '()';
		});

case 'sparse' % interpolator based on a sparse matrix
	if isempty(st.om), error 'sparse mode requires "om"', end
	st = newfft_init_sparse(st);

otherwise
	fail('unknown mode %s', st.mode)
end

if ~isempty(dotime)
	tmp = whos('st');
	tmp = num2str(tmp.bytes);
	tmp = ['newfft setup ' st.mode ' ' st.phasing(1) ' ' tmp ' ' dotime];
	cpu('etoc', tmp)
end

%
% newfft_init_sparse()
%
% create an interpolator based on a sparse matrix.
% this matrix will be large for intersting problem sizes, so this
% mode is not recommended.  but it is supported in part because it
% is needed for generating samples of the interpolator for the table mode.
%
function st = newfft_init_sparse(st)

om = st.om;

%
% different interpolation kernel mechanisms
%
ktype = st.ktype;

switch class(ktype)
case 'cell' % cell array of kernel functions: {kernel1, kernel2, ..., kernelD}
	if isa(ktype{1}, 'inline') || isa(ktype{1}, 'function_handle')
		if length(ktype) ~= dd, error 'wrong # of kernels', end
		st.kernel = ktype;
		ktype = 'inline';
	else
		error 'cell array should be inline kernels!?'
	end

case {'inline', 'function_handle'} % single inline kernel for all dimensions
	for id = 1:st.dd
		st.kernel{id} = ktype; % all same
	end
	ktype = 'inline';

case 'char'
	% a string that describes the type of interpolator, see below

otherwise
	fail('unknown kernel type %s', class(ktype))
end

%
% interpolator set up
%
Nd = st.Nd;
Jd = st.Jd;
Kd = st.Kd;
dd = st.dd;
is_kaiser_scale = false;

switch ktype
case 'inline'
	% already did it above

case 'diric' % exact interpolator
	if any(Jd ~= Kd), warn 'diric inexact unless Jd=Kd', end
	ktype = 'inline';
	for id = 1:dd
		N = Nd(id);
		K = Kd(id);
		if 1 && streq(st.phasing, 'real')
			N = 2 * floor((K+1)/2) - 1; % trick
		end
		st.kernel{id} = @(k,J) N / K * nufft_diric(k, N, K, true);
	end

case 'linear' % linear interpolator straw man
	ktype = 'inline';
	kernel = inline('(1 - abs(k/(J/2))) .* (abs(k) < J/2)', 'k', 'J');
	for id = 1:dd
		st.kernel{id} = kernel;
	end

case 'kb:beatty' % KB with Beatty et al parameters
	is_kaiser_scale = true;
	if ~isempty(st.kb_alf) || ~isempty(st.kb_m)
		warn 'kb_alf and kb_m ignored'
	end
	K_N = Kd ./ Nd;
	st.kb_alf = pi * sqrt( Jd.^2 ./ K_N.^2 .* (K_N - 1/2).^2 - 0.8 );
%	pr st.kb_alf ./ Jd % approximately 2.34 for K_N = 2
	st.kb_m = zeros(1,dd);
	for id = 1:dd
		st.kernel{id} = kaiser_bessel('inline', Jd(id), ...
				st.kb_alf(id), st.kb_m(id));
	end

case 'kb:minmax' % KB with minmax-optimized parameters
	is_kaiser_scale = true;

	if ~isempty(st.kb_alf) || ~isempty(st.kb_m)
		warn 'kb_alf and kb_m ignored'
	end
	for id = 1:dd
		[st.kernel{id} st.kb_alf(id) st.kb_m(id)] = ...
			kaiser_bessel('inline', Jd(id));
	end

case 'kb:user' % KB with user-defined parameters
	is_kaiser_scale = true;

	if isempty(st.kb_alf) || isempty(st.kb_m)
		fail 'kb_alf and kb_m required'
	end
	if (length(st.kb_alf) ~= dd) || (length(st.kb_m) ~= dd)
		fail('#alpha=%d #m=%d vs dd=%d', ...
			length(st.kb_alf), length(st.kb_m), dd)
	end
	for id = 1:dd
		st.kernel{id} = kaiser_bessel('inline', Jd(id), ...
				st.kb_alf(id), st.kb_m(id));
	end

case 'minmax:kb' % minmax interpolator with KB scaling factors
	for id = 1:dd
		[st.alpha{id}, st.beta{id}] = ...
			nufft_alpha_kb_fit(Nd(id), Jd(id), Kd(id), ...
				'Nmid', st.Nmid(id));
	end

case 'minmax:tuned' % minmax with numerically "tuned" scaling factors
	for id = 1:dd
		[st.alpha{id}, st.beta{id}, ok] = ...
			nufft_best_alpha(Jd(id), 0, Kd(id)/Nd(id));
		if ~ok, error 'unknown J,K/N', end
	end

case 'minmax:user' % minmax interpolator with user-provided scaling factors
	if isempty(st.alpha) || isempty(st.beta)
		error 'user must provide alpha/beta'
	end
	if length(st.alpha) ~= dd || length(st.beta) ~= dd
		error 'alpha/beta size mismatch'
	end

case 'minmax:unif' % minmax with straw man uniform scaling factors
	for id = 1:dd
		st.alpha{id} = 1;
		st.beta{id} = 0;
	end

otherwise
	fail('unknown kernel type %s', ktype)
end

%
% scaling factors: "outer product" of 1D vectors
%
st.sn = 1;
for id=1:dd
	if 1 && streq(st.ktype, 'linear')
		tmp = newfft_scale_tri(Nd(id), Jd(id), Kd(id), st.Nmid);
	elseif streq(st.ktype, 'diric')
		tmp = ones(Nd(id),1);
	elseif is_kaiser_scale
		nc = [0:Nd(id)-1]' - st.Nmid(id);
		tmp = 1 ./ kaiser_bessel_ft(...
			nc/Kd(id), Jd(id), st.kb_alf(id), st.kb_m(id), 1);
	elseif streq(ktype, 'inline')
		tmp = 1 ./ nufft_interp_zn(0, Nd(id), Jd(id), Kd(id), ...
			st.kernel{id}, st.Nmid(id));
	else
		tmp = nufft_scale(Nd(id), Kd(id), ...
			st.alpha{id}, st.beta{id}, st.Nmid(id));
	end
	tmp = reale(tmp);
	st.sn = st.sn(:) * tmp';
end
if length(Nd) > 1
	st.sn = reshape(st.sn, Nd);	% [(Nd)]
else
	st.sn = st.sn(:);	% [(Nd)]
end

%
% [J?,M] interpolation coefficient vectors.  will need kron of these later
%
for id=1:dd
	N = Nd(id);
	J = Jd(id);
	K = Kd(id);
	if isfield(st, 'kernel')
		[c, arg] = ...
		nufft_coef(om(:,id), J, K, st.kernel{id});	% [J?,M]
	else
		alpha = st.alpha{id};
		beta = st.beta{id};
		T = nufft_T(N, J, K, st.tol, alpha, beta);	% [J?,J?]
		[r, arg] = ...
		nufft_r(om(:,id), N, J, K, alpha, beta);	% [J?,M]
		c = T * r;
		clear T r
	end

	%
	% indices into oversampled FFT components
	%
	koff = nufft_offset(om(:,id), J, K);	% [M,1] to leftmost near nbr
	k0 = outer_sum([1:J]', koff');	% [J?,M] arbitrary integers
	kd{id} = mod(k0, K);	% [J?,M] {0,...,K?-1} (DFT indices)

	gam = 2*pi/K;

	switch st.phasing
	case {'real', 'none'}
		phase = 1;
	case 'complex'
		phase_scale = 1i * gam * (N-1)/2;
		phase = exp(phase_scale * arg);	% [J?,M] linear phase
	case 'flipreal'
		isodd = @(n) mod(n,2) == 1;
		phase = ones(size(k0)); % [J?,M]
		flip = isodd((kd{id} - k0) / K * (N-1)); % sign flip every K
		phase(flip) = -1; % sign flip (if N is even)
	otherwise
		fail('unknown phasing %s', st.phasing)
	end
	ud{id} = phase .* c;		% [J?,M]

end % id
clear c arg gam phase phase_scale koff k0 N J K

% st.M = size(om,1);
M = size(om,1);

%
% build sparse matrix that is [M,*Kd]
% with *Jd nonzero entries per frequency point
%
if st.printmem
	printm('Needs at least %g Gbyte RAM', prod(Jd)*M*8/2^30*2)
end

kk = kd{1};	% [J1,M]
uu = ud{1};	% [J1,M]
for id = 2:dd
	Jprod = prod(Jd(1:id));
	tmp = kd{id} * prod(Kd(1:(id-1)));
	kk = block_outer_sum(kk, tmp);		% outer sum of indices
	kk = reshape(kk, Jprod, M);
	uu = block_outer_prod(uu, ud{id});	% outer product of coefficients
	uu = reshape(uu, Jprod, M);
end % now kk and uu are [*Jd, M]

%
% handle phase shifts
%
uu = conj(uu); % [*Jd,M] ala Hermitian transpose of interpolation coefficients
switch st.phasing
case 'complex'
	phase = exp(1i * (om * st.n_shift(:))).';	% [1,M]
	uu = uu .* phase(ones(1,prod(Jd)),:);	% [*Jd,M]
	if streq(st.mode, 'table', 5) % moved from newfft_table_init.m to here
		st.phase_after = @(om) exp(1i * (om * col(st.n_shift))); % [M,1]
	end
case {'real', 'flipreal'}
	st.phase_before = newfft_phase_before(Kd, st.Nmid);
	if ~isempty(st.om) % precompute phase that goes after interpolation
		st.phase_after = newfft_phase_after(st.om, st.Nmid, st.n_shift);
	else
		st.phase_after = @(om) newfft_phase_after(om, st.Nmid, st.n_shift);
	end
case 'none'
	% do nothing
otherwise
	error 'bug'
end

mm = repmat(1:M, prod(Jd), 1); % [*Jd,M]
st.p = sparse(mm(:), 1+kk(:), uu(:), M, prod(Kd)); % [M, *Kd] sparse matrix
% sparse object, to better handle single precision operations!
st.p = Gsparse(st.p, 'odim', [M 1], 'idim', [prod(Kd) 1]);

st = strum(st, { ...
	'fft', @newfft_approx_for, '(x, [om])';
	'adj', @newfft_approx_adj, '(X, [om])';
	});

%
% in
%	x1	[J1,M]
%	x2	[J2,M]
% out
%	y	[J1,J2,M]	y(i1,i2,m) = x1(i1,m) + x2(i2,m)
%
function y = block_outer_sum(x1, x2)
[J1 M] = size(x1);
[J2 M] = size(x2);
xx1 = reshape(x1, [J1 1 M]);	% [J1,1,M] from [J1,M]
xx1 = xx1(:,ones(J2,1),:);	% [J1,J2,M], emulating ndgrid
xx2 = reshape(x2, [1 J2 M]);	% [1,J2,M] from [J2,M]
xx2 = xx2(ones(J1,1),:,:);	% [J1,J2,M], emulating ndgrid
y = xx1 + xx2;			% [J1,J2,M]

function y = block_outer_prod(x1, x2)
[J1 M] = size(x1);
[J2 M] = size(x2);
xx1 = reshape(x1, [J1 1 M]);	% [J1,1,M] from [J1,M]
xx1 = xx1(:,ones(J2,1),:);	% [J1,J2,M], emulating ndgrid
xx2 = reshape(x2, [1 J2 M]);	% [1,J2,M] from [J2,M]
xx2 = xx2(ones(J1,1),:,:);	% [J1,J2,M], emulating ndgrid
y = xx1 .* xx2;			% [J1,J2,M]

%
% newfft_phase_before()
% phase factor that gets multiplied by DFT (before interpolation)
%
function phase = newfft_phase_before(Kd, Nmid)

phase = 0;
for id = 1:length(Kd)
	tmp = 2 * pi * [0:Kd(id)-1] / Kd(id) * Nmid(id);
	phase = outer_sum(phase, tmp); % [(Kd)] when done
end
phase = exp(1i * phase);

%
% newfft_phase_after()
% phase factor that multiplies the DTFT (after interpolation)
%
function phase = newfft_phase_after(om, Nmid, n_shift)

phase = exp(1i * (om * col(n_shift - Nmid))); % [M,1]

%
% newfft_scale_tri()
% scale factors when kernel is 'linear'
% tri(u/J) <-> J sinc^2(J x)
%
function sn = newfft_scale_tri(N, J, K, Nmid)
nc = [0:N-1] - Nmid;
fun = @(x) J * nufft_sinc(J * x / K).^2;
cent = fun(nc);
sn = 1 ./ cent;

% try the optimal formula
tmp = 0;
LL = 3;
for ll=-LL:LL
	tmp = tmp + abs(fun(nc - ll*K)).^2;
end
sn = cent ./ tmp;

%
% newfft_test_time()
% compare compute times of real vs complex, sparse vs table
%
function newfft_test_time

Nd = [1 1] * 2^8;
[tmp om] = mri_trajectory('radial', {}, Nd, Nd);
pr length(om)
rng(0)
x0 = rand([Nd 1]);

arg = {om, Nd, 'dotime', ' '};
s0r = newfft(arg{:}, 'mode', 'table0');
s0c = newfft(arg{:}, 'mode', 'table0', 'phasing', 'complex');
s1r = newfft(arg{:}, 'mode', 'table1');
ssr = newfft(arg{:}, 'mode', 'sparse');
ssc = newfft(arg{:}, 'mode', 'sparse', 'phasing', 'complex');

tmp = @(st) newfft_test_time_one(st, x0, 2);
tmp(s0r)
tmp(s0c)
tmp(s1r)
tmp(ssr)
tmp(ssc)

function newfft_test_time_one(st, x0, nn)
st.fft(x0); % warm up
tmp = [st.mode ' ' st.phasing(1)];
cpu etic
for ii=1:nn
	st.fft(x0); % trial
end
cpu('etoc', tmp)

%
% newfft_test()
%
function newfft_test

newfft_test_time
prompt

modes = {'sparse', 'table0', 'table1'};

phasings = {'complex', 'real'};
% 'flipreal' no longer needed thanks to floor(N/2)
% 'none' is for internal table only

Nd_list = [20 10 8];
for id=0:3
	if id == 0
		Nd = 1 + Nd_list(1); % test odd case
	else
		Nd = Nd_list(1:id);
	end
	dd = length(Nd);

	rng(0)
	om = 3 * 2 * pi * (rand(100,length(Nd)) - 0.5);
%	om = sort(om);
%	om = linspace(-1,1,601)' * 3*pi;
%	om = [-8:8]'/2 * pi;
	%om = 'epi';
	x0 = rand([Nd 1]);
%	x0 = zeros([Nd 1]); x0(1) = 1; % unit vector for testing
%	x0 = ones([Nd 1]);

	% exact
	st_e = newfft(om, Nd, 'mode', 'exact');
	Xe = st_e.fft(x0);
	xe = st_e.adj(Xe);
	if 0
		pe = st_e.p;
		equivs(pe * x0(:), Xe)
		equivs(reshape(pe' * Xe, [Nd 1]), xe)
	end

	ktypes = {{'linear', 'Jd', 2*ones(1,dd)}, ...
		'minmax:unif', ...
		{'minmax:user', 'alpha', num2cell(ones(1,dd)), ...
			'beta', num2cell(0.5 * ones(1,dd))}, ...
		{'diric', 'Jd', 2*Nd-0, 'oversample', []}, ...
		'minmax:kb', 'minmax:tuned', ...
		'kb:minmax', 'kb:beatty', ...
		{'kb:user', 'kb_m', 0*Nd, 'kb_alf', 2.34 * 6 + 0*Nd}
		};

	for ii=4:length(ktypes) % skip poor ones
		ktype = ktypes{ii};
		if ~iscell(ktype), ktype = {ktype}; end

	for jj=1:length(modes)
	for ip=1:length(phasings)
		sc = newfft(st_e.om, st_e.Nd, 'phasing', 'complex', ...
			'mode', 'table0', 'ktype', ktype{:});

		st = newfft(st_e.om, st_e.Nd, 'phasing', phasings{ip}, ...
			'mode', modes{jj}, 'ktype', ktype{:});

		if streq(st.phasing, 'complex') && streq(st.mode, 'table1')
			continue
		end

%		pr minmax(st.sn)
		pad = @(s,n) [s blanks(n-length(s))];
		key = [sprintf('%2d ', st.Jd(1)) st.ktype];
		key = [st.mode ' ' num2str(id) st.phasing(1) ' ' pad(key,16)];
		Xs = st.fft(x0);
		max_percent_diff(Xe, Xs, key)
%		plot(abs(Xe), abs(Xs), 'o'), prompt

		xs = st.adj(Xe);
		max_percent_diff(xe, xs, key)
	end % ip
	end % jj
	end % ii
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [M,R] = ballatsq(N)
% BALLATSQ - Balanced Latin Square
%    M = BALLATSQ(N) creates a balanced latin square of size N containing
%    the numbers 1 to N. N should be an even positive integer.
%
%    [M, R] = BALLATSQ(N) also returns a randomized balanced latin square in R.
%
%    A latin square of size M is a MxM matrix filled with the M different
%    numbers in such a way that each number occurs exactly once in each row
%    and exactly once in each column. A balanced latin square has the further
%    restriction that each sequence of two values in a row does not occur
%    more than once in the whole matrix. They have applications in the
%    design of experiments.
%    More information: http://en.wikipedia.org/wiki/Latin_square
%
%    Example:
%      M = ballatsq(6) % ->
%      %      1     2     6     3     5     4
%      %      2     3     1     4     6     5
%      %      3     4     2     5     1     6
%      %      4     5     3     6     2     1
%      %      5     6     4     1     3     2
%      %      6     1     5     2     4     3
%
%     Note that "sort(ballatsq(N),1-2)" will return 1:N in each column-row.
%
%    See also MAGIC, GALLERY, CIRCSHIFT
%             LATSQ, CIRCULANT (FEX)

% for Matlab R13 and up
% version 2.2 (mar 2009)
% (c) Jos van der Geest
% email: jos@jasen.nl

% History
% 1.0 (feb 2006) created
% 2.0 (apr 2006) previous algorithm was faulty. Thanks to Jen M. (UK)
% 2.1 (sep 2006) can now also return a randomize balanced latin square
% 2.2 (mar 2009) added internal comments, changed variable names, improved
%     randomization

if nargin ~= 1 || ~isnumeric(N) || numel(N)~=1 || N<2 || rem(N,2) || fix(N)~=N,
    error('Single argument should be a positive even integer') ;
end

V = 1:N ;        % values in each column
CSI = [V ; -V] ; % circular shift index
M = zeros(N) ;   % pre-allocation

% the first colum just contains the numbers 1 to N
M(:,1) = V(:) ;
for i=2:N,
    % every column contains the numbers 1 to N shifted circularly according
    % to CSI
    M(:,i) =  rem(V(:)+CSI(i-1)+(N-1),N) + 1 ;
end

if nargout==2,
    % randomize the values 1:N
    rN = randperm(N) ;
    R = rN(M) ;
end

author: Sable

repo: mcbench-benchmarks


text:
clear;close all;clc;j=1i;
Global_Parameters;
%% Hardware Parameters
Mode='transmitRepeat'; % Select Mode
tx_object = sdrtx('ZedBoard and FMCOMMS2/3/4', ...
           'IPAddress',            '192.168.3.2', ...
           'CenterFrequency',      Parameters_struct.CenterFrequency, ...
           'BasebandSampleRate',   Parameters_struct.Bandwidth, ...  % Bandwidth
           'Gain',                 -10, ...
           'ChannelMapping',       1);
%          'EnableBurstMode',1,...

%% Button Setting
figure('Name','TX','NumberTitle','off');
TransmittingDisplay = uicontrol('Style', 'text', 'Position',[55,150,155,35],'String', 'Transmitting','FontSize',20,'HorizontalAlignment','left','BackgroundColor',[0.937 0.867 0.867]);
button = uicontrol; % Generate GUI button
set(button,'String','Stop !','Position',[80 50 100 60]); % Add "Stop !" text
set(gcf,'Units','centimeters','position',[3 3 7 6]); % Set the postion of GUI
%% TX Load
load('TX_signal'); % [1x972]
% transmitRepeat Mode
TX_Hardware = repmat(TX_signal.',5,1); % Transmit Data must be >= 4096 % [4860x1]
state = 1;
%% Main
switch Mode
    case 'step'
        while(state == 1)
           step(tx_object,TX_Hardware);
           % ----- Button Behavior -----%
           set(button,'Callback','setstate0_TX'); % Set the reaction of pushing button
           drawnow;
        end
        release(tx_object);

    case 'transmitRepeat'
        transmitRepeat(tx_object,TX_Hardware);
        % ----- Button Behavior -----%
        set(button,'Callback','setstate0_TX'); % Set the reaction of pushing button
end

author: MeowLucian


text:
function [X,Zf] = moving_average(N,X,Zi,dim)
% Like filter() for the special case of moving-average kernels.
% [X,Zf] = moving_average(N,X,Zi,Dim)
%
% This is an overall very fast implementation whose running time does now grow with N (beyond
% N=100). The algorithm does not run into numerical problems for large data sizes unlike the usual
% cumsum-based implementations.
%
% In:
%   N : filter length in samples
%
%   X : data matrix
%
%   Zi : initial filter conditions (default: [])
%
%   Dim : dimension along which to filter (default: first non-singleton dimension)
%
% Out:
%   X : the filtered data
%
%   Zf : final filter conditions
%
% See also:
%   filter
%
%                           Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                           2012-01-10

% Copyright (C) Christian Kothe, SCCN, 2012, christian@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify it under the terms of the GNU
% General Public License as published by the Free Software Foundation; either version 2 of the
% License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without
% even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License along with this program; if not,
% write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
% USA

% determine the dimension along which to filter
if nargin <= 3
    if isscalar(X)
        dim = 1;
    else
        dim = find(size(X)~=1,1);
    end
end

% empty initial state
if nargin <= 2
    Zi = []; end

lenx = size(X,dim);
if lenx == 0
    % empty X
    Zf = Zi;
else
    if N < 100
        % small N: use filter
        [X,Zf] = filter(ones(N,1)/N,1,X,Zi,dim);
    else
        % we try to avoid permuting dimensions below as this would increase the running time by ~3x
        if ndims(X) == 2
            if dim == 1
                % --- process along 1st dimension ---
                if isempty(Zi)
                    % zero initial state
                    Zi = zeros(N,size(X,2));
                elseif size(Zi,1) == N-1
                    % reverse engineer filter's initial state (assuming a moving average)
                    tmp = diff(Zi(end:-1:1,:),1,1);
                    Zi = [tmp(end:-1:1,:); Zi(end,:)]*N;
                    Zi = [-sum(Zi,1); Zi];
                elseif ~isequal(size(Zi),[N,size(X,2)])
                    error('These initial conditions do not have the correct format.');
                end

                % pre-pend initial state & get dimensions
                Y = [Zi; X]; M = size(Y,1);
                % get alternating index vector (for additions & subtractions)
                I = [1:M-N; 1+N:M];
                % get sign vector (also alternating, and includes the scaling)
                S = [-ones(1,M-N); ones(1,M-N)]/N;
                % run moving average
                X = cumsum(bsxfun(@times,Y(I(:),:),S(:)),1);
                % read out result
                X = X(2:2:end,:);

                % construct final state
                if nargout > 1
                    Zf = [-(X(end,:)*N-Y(end-N+1,:)); Y(end-N+2:end,:)]; end
            else
                % --- process along 2nd dimension ---
                if isempty(Zi)
                    % zero initial state
                    Zi = zeros(N,size(X,1));
                elseif size(Zi,1) == N-1
                    % reverse engineer filter's initial state (assuming a moving average)
                    tmp = diff(Zi(end:-1:1,:),1,1);
                    Zi = [tmp(end:-1:1,:); Zi(end,:)]*N;
                    Zi = [-sum(Zi,1); Zi];
                elseif ~isequal(size(Zi),[N,size(X,1)])
                    error('These initial conditions do not have the correct format.');
                end

                % pre-pend initial state & get dimensions
                Y = [Zi' X]; M = size(Y,2);
                % get alternating index vector (for additions & subtractions)
                I = [1:M-N; 1+N:M];
                % get sign vector (also alternating, and includes the scaling)
                S = [-ones(1,M-N); ones(1,M-N)]/N;
                % run moving average
                X = cumsum(bsxfun(@times,Y(:,I(:)),S(:)'),2);
                % read out result
                X = X(:,2:2:end);

                % construct final state
                if nargout > 1
                    Zf = [-(X(:,end)*N-Y(:,end-N+1)) Y(:,end-N+2:end)]'; end
            end
        else
            % --- ND array ---
            [X,nshifts] = shiftdim(X,dim-1);
            shape = size(X); X = reshape(X,size(X,1),[]);

            if isempty(Zi)
                % zero initial state
                Zi = zeros(N,size(X,2));
            elseif size(Zi,1) == N-1
                % reverse engineer filter's initial state (assuming a moving average)
                tmp = diff(Zi(end:-1:1,:),1,1);
                Zi = [tmp(end:-1:1,:); Zi(end,:)]*N;
                Zi = [-sum(Zi,1); Zi];
            elseif ~isequal(size(Zi),[N,size(X,2)])
                error('These initial conditions do not have the correct format.');
            end

            % pre-pend initial state & get dimensions
            Y = [Zi; X]; M = size(Y,1);
            % get alternating index vector (for additions & subtractions)
            I = [1:M-N; 1+N:M];
            % get sign vector (also alternating, and includes the scaling)
            S = [-ones(1,M-N); ones(1,M-N)]/N;
            % run moving average
            X = cumsum(bsxfun(@times,Y(I(:),:),S(:)),1);
            % read out result
            X = X(2:2:end,:);

            % construct final state
            if nargout > 1
                Zf = [-(X(end,:)*N-Y(end-N+1,:)); Y(end-N+2:end,:)]; end

            X = reshape(X,shape);
            X = shiftdim(X,ndims(X)-nshifts);
        end
    end
end

repo: REST


text:
function [cs,h]=m_contourf(long,lat,data,varargin)
%  M_CONTOURF Adds filled contours to a map
%    M_CONTOURF(LONG,LAT,DATA,...) is the same as M_CONTOUR except
%    that contours are filled. Areas of data above a given level are
%    filled, areas below are left blank or are filled by a lower level.
%    NaN's in the data leave holes in the filled plot/
%
%    [CS,H] = M_CONTOURF(...) returns contour matrix C as described in
%    CONTOURC and a vector H of handles to PATCH objects (for use by
%    CLABEL).
%
%    See also M_CONTOUR, CONTOURF

% Rich Pawlowicz (rich@ocgy.ubc.ca) 17/Jan/1998
%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.

% 19/02/98 - type - should have been 'clip','patch', rather than 'off'.
%  9/12/98 - handle all-NaN plots without letting contour crash.
% 6/Nov/00 - eliminate returned stuff if ';' neglected (thx to D Byrne)

global MAP_PROJECTION

% Have to have initialized a map first

if isempty(MAP_PROJECTION)
  disp('No Map Projection initialized - call M_PROJ first!');
  return;
end

if min(size(long))==1  &&  min(size(lat))==1
 [long,lat]=meshgrid(long,lat);
end

[X,Y]=m_ll2xy(long,lat,'clip','on');  %First find the points outside

i=isnan(X);      % For these we set the *data* to NaN...
data(i)=NaN;

                 % And then recompute positions without clipping. THis
                 % is necessary otherwise contouring fails (X/Y with NaN
                 % is a no-no. Note that this only clips properly down
                 % columns of long/lat - not across rows. In general this
                 % means patches may nto line up properly a right/left edges.
if any(i(:)), [X,Y]=m_ll2xy(long,lat,'clip','patch'); end

if any(~i(:))
 [cs,h]=contourf(X,Y,data,varargin{:});
 set(h,'tag','m_contourf');
else
  cs=[];h=[];
end

if nargout==0
 clear cs h
end

author: CelsoReyes


text: any(cellfun(@exist,{'mdwt', 'midwt', 'mirdwt', 'mrdwt'})==3)


text:
%MDL_P8 Create model of Puma robot on an XY base
%
% MDL_P8 is a script that creates the workspace variable p8 which is an
% 8-axis robot comprising a Puma 560 robot on an XY base.  Joints 1 and 2
% are the base, joints 3-8 are the robot arm.
%
% Also define the workspace vectors:
%   qz         zero joint angle configuration
%   qr         vertical 'READY' configuration
%   qstretch   arm is stretched out in the X direction
%   qn         arm is at a nominal non-singular configuration
%
% Notes::
% - SI units of metres are used.
%
% References::
% - Robotics, Vision & Control, 1st edn,
%   P. Corke, Springer 2011.  Sec 7.3.4.
%
% See also mdl_puma560, SerialLink.

% MODEL: Unimation, Puma560, on XY base, redundant, 8DOF, standard_DH

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

% create the base
platform = SerialLink( [0 0 0 -pi/2 1; -pi/2 0 0 pi/2 1], ...
    'base', troty(pi/2), 'name', 'platform' );

% load the standard Puma
mdl_puma560;
% add the pedestal height to link 1
p560.links(1).d = 30 * 0.0254;

% compose the two robots
p8 = SerialLink( [platform, p560], 'name', 'P8');

% set the range for the prismatic links
p8.links(1).qlim=[-1 1];
p8.links(2).qlim=[-1 1];

author: petercorke

repo: robotics-toolbox-matlab


text:
%% test modulate for all oasis functions.
col = {[0 114 178],[0 158 115], [213 94 0],[230 159 0],...
    [86 180 233], [204 121 167], [64 224 208], [240 228 66]}; % colors
plot_cvx = false;

%% example 3: foopsi, convolution kernel
g = [1.7, -0.712];         % AR coefficient
noise = 1;
T = 3000;
framerate = 30;
firerate = 0.5;
b = 0;              % baseline
N = 20;              % number of trials
seed = 3;          % seed for genrating random variables
[Y, trueC, trueS] = gen_data(g, noise, T, framerate, firerate, b, N, seed);
y = Y(1,:);
true_c = trueC(1,:);  %#ok<*NASGU>
true_s = trueS(1,:);
taus = ar2exp(g);
w = 200;
taus = ar2exp(g);
ht = exp2kernel(taus, w);
% case 1: use the difference of two exponential functions to construct a
% kernel
lambda = 0.1;
[c_oasis, s_oasis] = deconvolveCa(y, 'exp2', taus, 'foopsi', 'lambda', lambda, ...
    'shift', 100, 'window', 200);  %#ok<*ASGLU>

figure('name', 'FOOPSI, exp2, known: g, lambda', 'papersize', [15, 4]);
show_results;

% case 2: use the kernel directly
lambda = 0.1;
[c_oasis, s_oasis] = deconvolveCa(y, 'kernel', ht, 'foopsi', 'lambda', ...
    lambda, 'shift', 100, 'window', 200);  %#ok<*ASGLU>

figure('name', 'FOOPSI, kernel, known: g, lambda', 'papersize', [15, 4]);
show_results;

%% case 3: estimate the time constants
lambda = 0;
taus = ar2exp(g);
[c_oasis, s_oasis, options] = deconvolveCa(y, 'exp2', 'foopsi', 'lambda', lambda, ...
    'shift', 100, 'window', 200, 'smin', 0.5);  %#ok<*ASGLU>

figure('name', 'FOOPSI, exp2, known: g, lambda', 'papersize', [15, 4]);
show_results;

author: zhoupc


text:
function [V,E,I] = spline_to_poly(P,C,tol)
  % SPLINE_TO_POLY Evaluate a cubic Bezier spline as a polyline where each
  % segment corresponds to a locally flat segment of the curve up to given
  % tolerance
  %
  % [V,E] = spline_to_poly(P,C,tol)
  %
  % Inputs:
  %   P  #P by dim list of control point locations
  %   C  #C by 4 list of indices into P of cubic Bezier curves
  %   tol  tolerance
  % Outputs:
  %   V  #V by dim list of vertex locations
  %   E  #E by dim list of edge indices into V
  %   I  #E list of indices into 1:#C
  %

  V = P;
  E = [];
  I = [];
  % consider each cubic
  for c = 1:size(C,1)
    Pc = cubic_flat_eval(P(C(c,:),:),tol);
    J = [C(c,1) size(V,1)+(1:size(Pc,1)-2) C(c,4)];
    Ec = [J(1:end-1);J(2:end)]';
    E = [E; Ec];
    V = [V;Pc(2:end-1,:)];
    I = [I;repmat(c,size(Ec,1),1)];
  end
  [V,~,~,E] = remove_unreferenced(V,E);

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [targets,expdirs_chosen,framestarts,nframesperseg] = ...
  ChooseExampleIntervals(expdirs,classifierparamsfile,varargin)

if ~iscell(expdirs)
  expdirs = {expdirs};
end
nexps = numel(expdirs);

score_norm_prctile = 80;

% hysteresis for counting bouts
min_score_close = -.1;
max_dist_close = 5;
max_dist_open = 5;
max_score_open = .25;
nframesperseg = 500;
behaviorsuse = [];

% for weighing number of bouts of different behaviors
weightlogoff = .1;

% for weighting the same fly, same sex
downweight_samefly = 5;
downweight_samesex = 2;

% scheme for weighting
weightingscheme = 'sumlogsumscore';

score_thresh = .25;
[nintervals,score_norm_prctile,...
  min_score_close,max_dist_close,...
  max_dist_open,max_score_open,...
  nframesperseg,...
  downweight_samefly,...
  downweight_samesex,...
  weightlogoff,...
  behaviorsuse,...
  weightingscheme] = ...
  myparse(varargin,'nintervals',1,...
  'score_norm_prctile',score_norm_prctile,...
  'min_score_close',min_score_close,...
  'max_dist_close',max_dist_close,...
  'max_dist_open',max_dist_open,...
  'max_score_open',max_score_open,...
  'nframesperseg',nframesperseg,...
  'downweight_samefly',downweight_samefly,...
  'downweight_samesex',downweight_samesex,...
  'weightlogoff',weightlogoff,...
  'behaviorsuse',behaviorsuse,...
  'weightingscheme',weightingscheme);

classifierparams = ReadClassifierParamsFile(classifierparamsfile);
nbehaviors = numel(classifierparams);
behaviors = cell(1,nbehaviors);
scorefns = cell(1,nbehaviors);
for i = 1:nbehaviors,
  if iscell(classifierparams(i).behaviors.names),
    behaviors{i} = sprintf('%s_',classifierparams(i).behaviors.names{:});
    behaviors{i} = behaviors{i}(1:end-1);
  else
    behaviors{i} = classifierparams(i).behaviors.names;
  end
  behaviors{i} = regexprep(behaviors{i},'(_|^)([a-z])','${upper($2)}');
  scorefn = classifierparams(i).file.scorefilename;
  scorefns{i} = regexprep(scorefn,'\.mat$','');
end

%% create trx structure

%global STORED_TRX;
%if isempty(STORED_TRX),

trx = Trx('trxfilestr',classifierparams(1).file.trxfilename,...
  'moviefilestr',classifierparams(1).file.moviefilename,...
  'perframedir',classifierparams(1).file.perframedir);
if isfield(classifierparams,'perframe') && isfield(classifierparams(1).perframe,'params'),
  trx.SetPerFrameParams(classifierparams(1).perframe.params);
end
if isfield(classifierparams,'perframe') && isfield(classifierparams(1).perframe,'landmark_params'),
  trx.SetLandmarkParams(classifierparams(1).perframe.landmark_params);
end

for i = 1:nexps,
  trx.AddExpDir(expdirs{i},'openmovie',false);
end

% STORED_TRX = trx;
% else
%   trx = STORED_TRX;
% end

%% process scores

scores = cell(1,trx.nflies);
for fly = 1:trx.nflies,
  scores{fly} = nan(nbehaviors,trx(fly).nframes);
  for behi = 1:nbehaviors,
    scores{fly}(behi,:) = trx(fly).(scorefns{behi});
  end
end

% normalize
score_norms = nan(nbehaviors,1);
for behi = 1:nbehaviors,
  tmp = cell2mat(cellfun(@(x) x(behi,:),scores,'UniformOutput',false));
  score_norms(behi) = prctile(abs(tmp),score_norm_prctile);
end
for fly = 1:trx.nflies,
  for behi = 1:nbehaviors,
    scores{fly}(behi,:) = max(-1,min(1,scores{fly}(behi,:)/score_norms(behi,:)));
  end
end

labels = cell(1,trx.nflies);
for fly = 1:trx.nflies,

  labels{fly} = false(nbehaviors,trx(fly).nframes);
  for behi = 1:nbehaviors,
    label0 = scores{fly}(behi,:) > score_thresh;

    % close small holes where score >= min_score_close
    label_lenient = scores{fly}(behi,:) >= min_score_close;
    [t0s,t1s] = get_interval_ends(label0); t1s = t1s-1;
    gapdist = t0s(2:end)-t1s(1:end-1)-1;
    issmallgap = gapdist <= max_dist_close;
    label1 = label0;
    for i = find(issmallgap),
      doclose = all(label_lenient(t1s(i)+1:t0s(i+1)-1));
      if doclose,
        label1(t1s(i)+1:t0s(i+1)-1) = true;
      end
    end

    % remove short bouts where max score <= max_score_open
    [t0s,t1s] = get_interval_ends(label1); t1s = t1s-1;
    label_strict = scores{fly}(behi,:) > max_score_open;
    gapdist = t1s-t0s+1;
    issmallgap = gapdist <= max_dist_open;
    label2 = label1;
    for i = find(issmallgap),
      doopen = ~any(label_strict(t0s(i):t1s(i)));
      if doopen,
        label2(t0s(i):t1s(i)) = false;
      end
    end
    labels{fly}(behi,:) = label2;

  end
end

%%

% count the number of bout starts and ends in each interval
fil = ones(1,nframesperseg);
nbouts = cell(1,trx.nflies);
w0 = cell(1,trx.nflies);
for fly = 1:trx.nflies,

  isstart = double(~labels{fly}(:,1:end-1) & labels{fly}(:,2:end));
  isend = double(labels{fly}(:,1:end-1) & ~labels{fly}(:,2:end));
  nbouts{fly} = .5*imfilter(isstart,fil,0,'same') + ...
    .5*imfilter(isend,fil,0,'same');

  switch weightingscheme,
    case 'sumlogsumscore',
      logsumscores = nan(size(labels{fly}));
      minv = log(weightlogoff);
      logsumscores(:) = minv;
      for behi = 1:nbehaviors,
        [t0s,t1s] = get_interval_ends(labels{fly}(behi,:)); t1s=t1s-1;
        % csscores(i) is the sum of scores from 1:i-1
        csscores = [0,cumsum(scores{fly}(behi,:))];
        sumscores = csscores(t1s+1)-csscores(t0s);
        lsscurr = log(max(0,sumscores)+weightlogoff)./(t1s-t0s+1);
        for i = 1:numel(t0s),
          logsumscores(behi,t0s(i):t1s(i)) = lsscurr(i);
        end
      end
      logsumscores = logsumscores-minv;
      sumlogsumscore = imfilter(logsumscores,fil,0,'same');
      w0{fly} = sumlogsumscore;
    case 'lognbouts',
      w0{fly} = log(nbouts{fly}+weightlogoff);
    case 'isbout',
      w0{fly} = double(nbouts{fly} >= 1);
    case 'none',
      w0{fly} = ones(size(nbouts{fly}));
    case 'random',
      w0{fly} = rand(size(nbouts{fly}));
  end
end

% make sure weights are non-negative
if ~strcmpi(weightingscheme,'none'),
  tmp = min(min(cell2mat(w0),[],1),[],2);
  w0 = cellfun(@(x) x-tmp,w0,'UniformOutput',false);
end
%
% % normalize per-behavior
% Z = zeros(nbehaviors,1);
% for behi = 1:nbehaviors,
%   Z(behi) = sum(cell2mat(cellfun(@(x) x(behi,:),w0,'UniformOutput',false)));
% end
% Z = Z / max(Z);
% w0 = cellfun(@(x) bsxfun(@rdivide,x,Z),w0,'UniformOutput',false);

% make sure sex classification is uniform for each bout
ismale = cellfun(@(x) strcmpi(x,'M'),trx.sex,'UniformOutput',false);
isallowed = cell(1,trx.nflies);
for fly = 1:trx.nflies,
  nmale = imfilter(double(ismale{fly}),fil,nan,'same');
  isallowed{fly} = nmale == 0 | nmale == nframesperseg;
end

w0mat = cell2mat(w0);
w0mat(:,~cell2mat(isallowed)) = nan;
nboutsmat = cell2mat(nbouts);
ismalemat = cell2mat(ismale);
flyidx = zeros(1,numel(ismalemat));
flyoffs = cumsum([0,trx.nfliespermovie(1:end-1)]);
flystarts = cumsum([1,trx(1:trx.nflies-1).nframes]);
flyidx(flystarts) = 1;
flyidx = cumsum(flyidx);

framestarts = nan(1,nintervals);
targets = nan(1,nintervals);
expis = nan(1,nintervals);

% only care about some behaviors
if ~isempty(behaviorsuse)
  behaviorsuseidx = ismember(behaviors,behaviorsuse);
else
  behaviorsuseidx = true(1,nbehaviors);
end
for segi = 1:nintervals,

  % choose the highest weight interval
  [~,i] = max(sum(w0mat(behaviorsuseidx,:),1));
  ismalecurr = ismalemat(i);
  fly = flyidx(i);
  i0 = i-flystarts(fly)+1-floor(nframesperseg/2);
  t0 = i0-trx(fly).off;
  framestarts(segi) = t0;
  expi = trx.fly2exp(fly);
  flyi = fly-flyoffs(expi);
  targets(segi) = flyi;
  expis(segi) = expi;

  % downweight this fly
  idx = flyidx == fly;
  w0mat(:,idx) = w0mat(:,idx) / downweight_samefly;

  % downweight this fly's sex
  idx = ismalemat == ismalecurr;
  w0mat(:,idx) = w0mat(:,idx) / downweight_samesex;

  % downweight behaviors selected
  w0mat = bsxfun(@rdivide,w0mat,(nboutsmat(:,i)+1));

  % don't allow overlapping bouts
  i00 = max(i0-nframesperseg,1)+flystarts(fly)-1;
  i10 = min(i0+2*nframesperseg,trx(fly).nframes)+flystarts(fly)-1;
  w0mat(:,i00:i10) = nan;

end

expdirs_chosen = expdirs(expis);

author: kristinbranson

repo: JAABA


text:
function [Ta,Ka,Ba,Tg,Kg,Bg,Tm2a,Bm,Vm,mag_strength]=ImuCalibration_Gesture(data)
% input data raw IMU data from mpu9250
% data :time accelerometer  gyroscope   magnetometer
%  cal_acc=Ta*Ka*(raw_acc+Ba)
%  cal_gyro=Tg*Kg*(raw_gyro+Bg)
%  cal_mag=Tm2a*(raw_mag+Bm)
%
% author  Zhang Xin

[~,fix_point,rotation]=FindFixData(data,30);

[Ta,Ka,Ba]=ICRA2014_acc(fix_point);

Bg=-mean(fix_point(:,4:6),1)';

n=size(rotation,1);

rotation{n+1}=Ta;
rotation{n+2}=Ka;
rotation{n+3}=Ba;
rotation{n+4}=Bg;

[Tg,Kg]=ICRA_2014_gyro(rotation);

%[Tm2a,Bm,Vm]=mag2acc_matrix(fix_point,Ta,Ka,Ba);

[Tm2a,Bm,Vm,mag_strength]=Cal_mag4acc_frame(rotation,fix_point,Tg,Kg);

Set_Bias_Gyro=[0.1,-0.4,1.5];

See_Gesture( data,Ta,Ka,Ba,Tg,Kg,Bg,Tm2a,Bm,Vm,Set_Bias_Gyro);

end

author: shenshikexmu

repo: IMUCalibration-Gesture


text:
function r4vec_uniform_ab_test ( )

%*****************************************************************************80
%
%% R4VEC_UNIFORM_AB_TEST tests  R4VEC_UNIFORM_AB.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2014
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  a = -1.0;
  b = +5.0;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R4VEC_UNIFORM_AB_TEST\n' );
  fprintf ( 1, '  R4VEC_UNIFORM_AB computes a random R4VEC.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  %g <= x <= %g\n', a, b );
  fprintf ( 1, '  Initial seed is %d\n', seed );

  [ v, seed ] = r4vec_uniform_ab ( n, a, b, seed );

  r4vec_print ( n, v, '  Uniform R4VEC:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
filename='Gripping_triangle_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'nonadjoint_compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
constraint_case = 'INEQUALITY';
optimizer = 'SLERP'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [r,bb,f, han] = timeseries_prplot(y,X,cols,varargin)
% :Usage:
% ::
%
% Partial residual plot of a subset [cols] of columns regressing y on X
%
%    [r,bb,f] = timeseries_prplot(y,X,cols,varargin)
%
% - Plots timeseries data (y') removing other regressors (~cols) against partial fitted response (X[cols])
% - Makes sure there is an intercept at the end, whether one was entered or not.
% - This is *not* a partial regression plot
%
% See also: https://www.itl.nist.gov/div898/software/dataplot/refman1/auxillar/partregr.htm#:~:text=Partial%20regression%20plots%20are%20most,independent%20variables%20in%20the%20model).
% "Partial regression plots are most commonly used to identify leverage
% points and influential data points that might not be leverage points.
% Partial residual plots are most commonly used to identify the nature of the
% relationship between Y and Xi (given the effect of the other independent variables
% in the model)"
%
% :Inputs:
%
%   **y:**
%        y is n x 1 data points, X is n x k model matrix of predictors
%
%        y is adjusted to remove all effects OTHER THAN those columns of X
%        specified in cols, creating a partial residual vector y'.
%
%   **X:**
%        The fitted response to X(:,cols) is plotted against y' to graphically assess the
%        effect of particular columns of X on y.
%
%   varargin includes two optional arguments (enter both *together*):
%   These will create "boxcars" showing trials onsets, e.g., for fMRI
%     1) a vector of trial onsets (to shade in plot)
%     2) length of elements to shade after trial onset.
%
% :Outputs:
%
%   **r:**
%        partial residuals
%
%   **bb:**
%        betas
%
%   **f:**
%        partial fitted response
%
% :Examples:
% ::
%
%    timeseries_prplot(Yfla,X,[2 4],x2,18);
%
% to average across sessions 1 and 2:
%
% ..
%    Tor Wager
% ..

% Make sure there is an intercept at the end, whether one was entered or not.
X = intercept(X, 'end');

% get betas from full model
b = pinv(X) * y;
bb = b;

% get fit from subset of predictors of interest here

f = X(:,cols) * b(cols);

% get partial residuals, excluding the subset of interest

X(:,cols) = []; b(cols) = [];
r = y - X * b;

% plot
tor_fig; hold on;
han = plot(r,'k','LineWidth', 2);
han(2) = plot(f,'b','LineWidth', 2);
legend({'Partial residual (data)' 'Partial fit'});

if length(varargin) > 0
    x2 = varargin{1}; len = varargin{2};
    sc = 2 * max(abs(r));

    for i = 1:length(x2), hh(i) = fill([x2(i) x2(i)+len x2(i)+len, x2(i)],sc*([0 0 1 1]-.5),[.5 .5 .5]); end

    set(hh,'EdgeColor','none','FaceAlpha',.5)
end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [ room ] = convert_room_rep( room, labelset )

% convert the objects of each room to the [obb+label] format
% INPUT: room - the room extract from .json file
%        labelset - the considered label set. Labels in this set use
%        one-hot vector, others use zero vectors
% OUTPUT: the room with updated obblist, labellist, objlist

obblist = [];
objlist = room.objlist;
for i = 1:length(objlist)
    obj = objlist{i};

    % local directions
    lup = obj.up;
    lfront = obj.front;
    laxis = cross(lup,lfront);
    laxis = laxis/norm(laxis,2);
    localbbox = obj.localbbox;

    transform = reshape(obj.transform,4,4);
    for j = 1:3
        transform(1:3,j) = transform(1:3,j)/norm(transform(1:3,j),2);
    end

    % local and global max&min coordinates
    globbox = obj.bbox;
    scale = norm(globbox.max - globbox.min,2)/norm(localbbox.max-localbbox.min,2);
    lminp = localbbox.min;
    lmaxp = localbbox.max;
    ldia = lmaxp - lminp;
    gminp = globbox.min;
    gmaxp = globbox.max;

    % compute the global directions
    gfront = transform*[lfront(:);0];
    gfront = gfront(1:3);
    gfront = gfront/norm(gfront,2);
    gup = transform*[lup(:);0];
    gup = gup(1:3);
    gup = gup/norm(gup,2);
    gaxis = cross(gup,gfront);
    gaxis = gaxis/norm(gaxis,2);

    % compute the sizes along local axes
    fsize = abs(ldia*lfront(:))*scale;
    usize = abs(ldia*lup(:))*scale;
    asize = abs(ldia*laxis(:))*scale;
    cent = (gmaxp+gminp)/2;

    % label vector
    label = zeros(length(labelset),1);
    ind = strmatch(obj.label,labelset,'exact');
    if(length(ind)>0)
        label(ind)=1;
    end

    vec = [cent(:);gfront(:);gup(:);fsize;usize;asize;label(:)];
    obblist(:,i) = vec;
end

iskept = sum(obblist(end-length(labelset)+1:end,:));
ind = find(iskept>0);
room.obblist = obblist(:,ind);
room.objlist = objlist(ind);
room.labellist = room.labellist(ind);

end

repo: GRAINS


text:
function [result] = oneMatrix(size, isGPU, type)
%%%
%
% Initialize zero matrix with GPU supprot
%
% Thang Luong @ 2015, <lmthang@stanford.edu>
%
%%%
  if isGPU
    result = ones(size, type, 'gpuArray');
  else
    result = ones(size, type);
  end
end

author: lmthang


text:
function b = onRay(point, ray)
%ONRAY test if a point belongs to a ray
%
%   B = onRay(PT, RAY);
%   Returns 1 if point PT belongs to the ray RAY.
%   PT is given by [x y] and RAY by [x0 y0 dx dy].
%
%   See also:
%   rays2d, points2d, onLine
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 31/10/2003.
%

%   HISTORY
%   07/07/2005 : normalize condition to test if on the line
%       and add support of multiple rays or points
%   22/05/2009 deprecate

% deprecation warning
warning('geom2d:deprecated', ...
    '''onRay'' is deprecated, use ''isPointOnRay'' instead');

% number of rays and points
Nr = size(line, 1);
Np = size(point, 1);

% if several rays or several points, adapt sizes of arrays
x0 = repmat(ray(:,1)', Np, 1);
y0 = repmat(ray(:,2)', Np, 1);
dx = repmat(ray(:,3)', Np, 1);
dy = repmat(ray(:,4)', Np, 1);
xp = repmat(point(:,1), 1, Nr);
yp = repmat(point(:,2), 1, Nr);

% test if points belongs to the ray
b1 = abs((xp-x0).*dy-(yp-y0).*dx)./sqrt(dx.*dx+dy.*dy) < 1e-13;

% check if points lie the good direction on the rays
ind = abs(dx)>abs(dy);
t = zeros(size(b1));
t(ind) = (xp(ind)-x0(ind))./dx(ind);
t(~ind) = (yp(~ind)-y0(~ind))./dy(~ind);

% combine the two tests
b = b1 && (t>0);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
classdef BoundaryMeshCreator < handle

    methods (Access = public, Static)

        function obj = create(cParams)
            f = BoundaryMeshCreatorFactory();
            obj = f.create(cParams);
        end

    end

end


text:
function [epoch_std, year] = mydatestd (epoch)
    temp = mydatevec(epoch);
    year = temp(:,1);
    [year_len, year_start] = mydatestd_aux(year);

    std_year_len = mydatestd_aux();

    epoch_std = (epoch - year_start) ./ year_len .* std_year_len;
end

%!test
%! % mydatestd()
%! test('mydatestd_aux')


text:
function points = mergeClosePoints(points, varargin)
%MERGECLOSEPOINTS Merge points that are closer than a given distance.
%
%   PTS2 = mergeClosePoints(PTS, DIST)
%   Remove points in the array PTS such that no points closer than the
%   distance DIST remain in the array.
%
%   PTS2 = mergeClosePoints(PTS)
%   If the distance is not specified, the default value 1e-14 is used.
%
%
%   Example
%     pts = rand(200, 2);
%     pts2 = mergeClosePoints(pts, .1);
%     figure; drawPoint(pts, '.');
%     hold on; drawPoint(pts2, 'mo');
%
%   See also
%     points2d, removeMultipleVertices

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2013-10-04, using Matlab 7.9.0.529 (R2009b)
% Copyright 2013-2022 INRA - Cepia Software Platform

% default values
minDist = 1e-14;
if ~isempty(varargin)
    minDist = varargin{1};
end

i = 1;
while i < size(points, 1)
    dist = distancePoints(points(i,:), points);
    inds = dist < minDist;
    inds(i) = 0;

    points(inds, :) = [];

    % switch to next point
    i = i + 1;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function c = multiprod(a, b, idA, idB)
%MULTIPROD  Multiplying 1-D or 2-D subarrays contained in two N-D arrays.
%   C = MULTIPROD(A,B) is equivalent  to C = MULTIPROD(A,B,[1 2],[1 2])
%   C = MULTIPROD(A,B,[D1 D2]) is eq. to C = MULTIPROD(A,B,[D1 D2],[D1 D2])
%   C = MULTIPROD(A,B,D1) is equival. to C = MULTIPROD(A,B,D1,D1)
%
%   MULTIPROD performs multiple matrix products, with array expansion (AX)
%   enabled. Its first two arguments A and B are "block arrays" of any
%   size, containing one or more 1-D or 2-D subarrays, called "blocks" (*).
%   For instance, a 5?6?3 array may be viewed as an array containing five
%   6?3 blocks. In this case, its size is denoted by 5?(6?3). The 1 or 2
%   adjacent dimensions along which the blocks are contained are called the
%   "internal dimensions" (IDs) of the array (?).
%
%   1) 2-D by 2-D BLOCK(S) (*)
%         C = MULTIPROD(A, B, [DA1 DA2], [DB1 DB2]) contains the products
%         of the P?Q matrices in A by the R?S matrices in B. [DA1 DA2] are
%         the IDs of A; [DB1 DB2] are the IDs of B.
%
%   2) 2-D by 1-D BLOCK(S) (*)
%         C = MULTIPROD(A, B, [DA1 DA2], DB1) contains the products of the
%         P?Q matrices in A by the R-element vectors in B. The latter are
%         considered to be R?1 matrices. [DA1 DA2] are the IDs of A; DB1 is
%         the ID of B.
%
%   3) 1-D by 2-D BLOCK(S) (*)
%         C = MULTIPROD(A, B, DA1, [DB1 DB2]) contains the products of the
%         Q-element vectors in A by the R?S matrices in B. The vectors in A
%         are considered to be 1?Q matrices. DA1 is the ID of A; [DB1 DB2]
%         are the IDs of B.
%
%   4) 1-D BY 1-D BLOCK(S) (*)
%      (a) If either SIZE(A, DA1) == 1 or SIZE(B, DB1) == 1, or both,
%             C = MULTIPROD(A, B, DA1, DB1) returns products of scalars by
%             vectors, or vectors by scalars or scalars by scalars.
%      (b) If SIZE(A, DA1) == SIZE(B, DB1),
%             C = MULTIPROD(A, B, [0 DA1], [DB1 0]) or
%             C = MULTIPROD(A, B, DA1, DB1) virtually turns the vectors
%             contained in A and B into 1?P and P?1 matrices, respectively,
%             then returns their products, similar to scalar products.
%             Namely, C = DOT2(A, B, DA1, DB1) is equivalent to
%             C = MULTIPROD(CONJ(A), B, [0 DA1], [DB1 0]).
%      (c) Without limitations on the length of the vectors in A and B,
%             C = MULTIPROD(A, B, [DA1 0], [0 DB1]) turns the vectors
%             contained in A and B into P?1 and 1?Q matrices, respectively,
%             then returns their products, similar to outer products.
%             Namely, C = OUTER(A, B, DA1, DB1) is equivalent to
%             C = MULTIPROD(CONJ(A), B, [DA1 0], [0 DB1]).
%
%   Common constraints for all syntaxes:
%      The external dimensions of A and B must either be identical or
%      compatible with AX rules. The internal dimensions of each block
%      array must be adjacent (DA2 == DA1 + 1 and DB2 == DB1 + 1 are
%      required). DA1 and DB1 are allowed to be larger than NDIMS(A) and
%      NDIMS(B). In syntaxes 1, 2, and 3, Q == R is required, unless the
%      blocks in A or B are scalars.
%
%   Array expansion (AX):
%      AX is a powerful generalization to N-D of the concept of scalar
%      expansion. Indeed, A and B may be scalars, vectors, matrices or
%      multi-dimensional arrays. Scalar expansion is the virtual
%      replication or annihilation of a scalar which allows you to combine
%      it, element by element, with an array X of any size (e.g. X+10,
%      X*10, or []-10). Similarly, in MULTIPROD, the purpose of AX is to
%      automatically match the size of the external dimensions (EDs) of A
%      and B, so that block-by-block products can be performed. ED matching
%      is achieved by means of a dimension shift followed by a singleton
%      expansion:
%      1) Dimension shift (see SHIFTDIM).
%            Whenever DA1 ~= DB1, a shift is applied to impose DA1 == DB1.
%            If DA1 > DB1, B is shifted to the right by DA1 - DB1 steps.
%            If DB1 > DA1, A is shifted to the right by DB1 - DA1 steps.
%      2) Singleton expansion (SX).
%            Whenever an ED of either A or B is singleton and the
%            corresponding ED of the other array is not, the mismatch is
%            fixed by virtually replicating the array (or diminishing it to
%            length 0) along that dimension.
%
%   MULTIPROD is a generalization for N-D arrays of the matrix
%   multiplication function MTIMES, with AX enabled. Vector inner, outer,
%   and cross products generalized for N-D arrays and with AX enabled are
%   performed by DOT2, OUTER, and CROSS2 (MATLAB Central, file #8782).
%   Elementwise multiplications (see TIMES) and other elementwise binary
%   operations with AX enabled are performed by BAXFUN (MATLAB Central,
%   file #23084). Together, these functions make up the ?ARRAYLAB toolbox?.
%
%   Input and output format:
%      The size of the EDs of C is determined by AX. Block size is
%      determined as follows, for each of the above-listed syntaxes:
%      1) C contains P?S matrices along IDs MAX([DA1 DA2], [DB1 DB2]).
%      2) Array     Block size     ID(s)
%         ----------------------------------------------------
%         A         P?Q  (2-D)     [DA1 DA2]
%         B         R    (1-D)     DB1
%         C (a)     P    (1-D)     MAX(DA1, DB1)
%         C (b)     P?Q  (2-D)     MAX([DA1 DA2], [DB1 DB1+1])
%         ----------------------------------------------------
%         (a) The 1-D blocks in B are not scalars (R > 1).
%         (b) The 1-D blocks in B are scalars (R = 1).
%      3) Array     Block size     ID(s)
%         ----------------------------------------------------
%         A           Q  (1-D)     DA1
%         B         R?S  (2-D)     [DB1 DB2]
%         C (a)       S  (1-D)     MAX(DA1, DB1)
%         C (b)     R?S  (2-D)     MAX([DA1 DA1+1], [DB1 DB2])
%         ----------------------------------------------------
%         (a) The 1-D blocks in A are not scalars (Q > 1).
%         (b) The 1-D blocks in A are scalars (Q = 1).
%      4)     Array     Block size         ID(s)
%         --------------------------------------------------------------
%         (a) A         P        (1-D)     DA1
%             B         Q        (1-D)     DB1
%             C         MAX(P,Q) (1-D)     MAX(DA1, DB1)
%         --------------------------------------------------------------
%         (b) A         P        (1-D)     DA1
%             B         P        (1-D)     DB1
%             C         1        (1-D)     MAX(DA1, DB1)
%         --------------------------------------------------------------
%         (c) A         P        (1-D)     DA1
%             B         Q        (1-D)     DB1
%             C         P?Q      (2-D)     MAX([DA1 DA1+1], [DB1 DB1+1])
%         --------------------------------------------------------------
%
%   Terminological notes:
%   (*) 1-D and 2-D blocks are generically referred to as "vectors" and
%       "matrices", respectively. However, both may be also called
%       ?scalars? if they have a single element. Moreover, matrices with a
%       single row or column (e.g. 1?3 or 3?1) may be also called ?row
%       vectors? or ?column vectors?.
%   (?) Not to be confused with the "inner dimensions" of the two matrices
%       involved in a product X * Y, defined as the 2nd dimension of X and
%       the 1st of Y (DA2 and DB1 in syntaxes 1, 2, 3).
%
%   Examples:
%    1) If  A is .................... a 5?(6?3)?2 array,
%       and B is .................... a 5?(3?4)?2 array,
%       C = MULTIPROD(A, B, [2 3]) is a 5?(6?4)?2 array.
%
%       A single matrix A pre-multiplies each matrix in B
%       If  A is ........................... a (1?3)    single matrix,
%       and B is ........................... a 10?(3?4) 3-D array,
%       C = MULTIPROD(A, B, [1 2], [3 4]) is a 10?(1?4) 3-D array.
%
%       Each matrix in A pre-multiplies each matrix in B (all possible
%       combinations)
%       If  A is .................... a (6?3)?5   array,
%       and B is .................... a (3?4)?1?2 array,
%       C = MULTIPROD(A, B, [1 2]) is a (6?4)?5?2 array.
%
%   2a) If  A is ........................... a 5?(6?3)?2 4-D array,
%       and B is ........................... a 5?(3)?2   3-D array,
%       C = MULTIPROD(A, B, [2 3], [2]) is   a 5?(6)?2   3-D array.
%
%   2b) If  A is ........................... a 5?(6?3)?2 4-D array,
%       and B is ........................... a 5?(1)?2   3-D array,
%       C = MULTIPROD(A, B, [2 3], [2]) is   a 5?(6?3)?2 4-D array.
%
%   4a) If both A and B are .................. 5?(6)?2   3-D arrays,
%       C = MULTIPROD(A, B, 2) is .......... a 5?(1)?2   3-D array, while
%   4b) C = MULTIPROD(A, B, [2 0], [0 2]) is a 5?(6?6)?2 4-D array
%
%   See also DOT2, OUTER, CROSS2, BAXFUN, MULTITRANSP, MULTITRACE, MULTISCALE.

% $ Version: 2.1 $
% CODE      by:            Paolo de Leva
%                          (Univ. of Rome, Foro Italico, IT)    2009 Jan 24
%           optimized by:  Paolo de Leva
%                          Jinhui Bai (Georgetown Univ., D.C.)  2009 Jan 24
% COMMENTS  by:            Paolo de Leva                        2009 Feb 24
% OUTPUT    tested by:     Paolo de Leva                        2009 Feb 24
% -------------------------------------------------------------------------

%error( nargchk(2, 4, nargin) ); % Allow 2 to 4 input arguments % HK
narginchk(2, 4); % added by HK
switch nargin % Setting IDA and/or IDB
    case 2, idA = [1 2]; idB = [1 2];
    case 3, idB = idA;
end

% ESC 1 - Special simple case (both A and B are 2D), solved using C = A * B

     if ndims(a)==2 && ndims(b)==2 && ...
         isequal(idA,[1 2]) && isequal(idB,[1 2])
         c = a * b; return
     end

% MAIN 0 - Checking and evaluating array size, block size, and IDs

     sizeA0 = size(a);
     sizeB0 = size(b);
     [sizeA, sizeB, shiftC, delC, sizeisnew, idA, idB, ...
     squashOK, sxtimesOK, timesOK, mtimesOK, sumOK] = ...
                                           sizeval(idA,idB, sizeA0,sizeB0);

% MAIN 1 - Applying dimension shift (first step of AX) and
%          turning both A and B into arrays of either 1-D or 2-D blocks

     if sizeisnew(1), a = reshape(a, sizeA); end
     if sizeisnew(2), b = reshape(b, sizeB); end

% MAIN 2 - Performing products with or without SX (second step of AX)

     if squashOK % SQUASH + MTIMES (fastest engine)
         c = squash2D_mtimes(a,b, idA,idB, sizeA,sizeB, squashOK);
     elseif timesOK % TIMES (preferred w.r. to SX + TIMES)
         if sumOK, c = sum(a .* b, sumOK);
         else      c =     a .* b; end
     elseif sxtimesOK % SX + TIMES
         if sumOK, c = sum(bsxfun(@times, a, b), sumOK);
         else      c =     bsxfun(@times, a, b); end
     elseif mtimesOK % MTIMES (rarely used)
         c = a * b;
     end

% MAIN 3 - Reshaping C (by inserting or removing singleton dimensions)

     [sizeC sizeCisnew] = adjustsize(size(c), shiftC, false, delC, false);
     if sizeCisnew, c = reshape(c, sizeC); end

function c = squash2D_mtimes(a, b, idA, idB, sizeA, sizeB, squashOK)
% SQUASH2D_MTIMES  Multiproduct with single-block expansion (SBX).
%    Actually, no expansion is performed. The multi-block array is
%    rearranged from N-D to 2-D, then MTIMES is applied, and eventually the
%    result is rearranged back to N-D. No additional memory is required.
%    One and only one of the two arrays must be single-block, and its IDs
%    must be [1 2] (MAIN 1 removes leading singletons). Both arrays
%    must contain 2-D blocks (MAIN 1 expands 1-D blocks to 2-D).

    if squashOK == 1 % A is multi-block, B is single-block (squashing A)

        % STEP 1 - Moving IDA(2) to last dimension
        nd = length(sizeA);
        d2 = idA(2);
        order = [1:(d2-1) (d2+1):nd d2]; % Partial shifting
        a = permute(a, order); % ...?Q

        % STEP 2 - Squashing A from N-D to 2-D
        q = sizeB(1);
        s = sizeB(2);
        lengthorder = length(order);
        collapsedsize = sizeA(order(1:lengthorder-1));
        n = prod(collapsedsize);
        a = reshape(a, [n, q]); % N?Q
        fullsize = [collapsedsize s]; % Size to reshape C back to N-D

    else % B is multi-block, A is single-block (squashing B)

        % STEP 1 - Moving IDB(1) to first dimension
        nd = length(sizeB);
        d1 = idB(1);
        order = [d1 1:(d1-1) (d1+1):nd]; % Partial shifting
        b = permute(b, order); % Q?...

        % STEP 2 - Squashing B from N-D to 2-D
        p = sizeA(1);
        q = sizeA(2);
        lengthorder = length(order);
        collapsedsize = sizeB(order(2:lengthorder));
        n = prod(collapsedsize);
        b = reshape(b, [q, n]); % Q?N
        fullsize = [p collapsedsize]; % Size to reshape C back to N-D

    end

    % FINAL STEPS - Multiplication, reshape to N-D, inverse permutation
    invorder(order) = 1 : lengthorder;
    c = permute (reshape(a*b, fullsize), invorder);

function [sizeA, sizeB, shiftC, delC, sizeisnew, idA, idB, ...
          squashOK, sxtimesOK, timesOK, mtimesOK, sumOK] = ...
                                          sizeval(idA0,idB0, sizeA0,sizeB0)
%SIZEVAL   Evaluation of array size, block size, and IDs
%    Possible values for IDA and IDB:
%        [DA1 DA2], [DB1 DB2]
%        [DA1 DA2], [DB1]
%        [DA1],     [DB1 DB2]
%        [DA1],     [DB1]
%        [DA1 0],   [0 DB1]
%        [0 DA1],   [DB1 0]
%
%    sizeA/B     Equal to sizeA0/B0 if RESHAPE is not needed in MAIN 1
%    shiftC, delC    Variables controlling MAIN 3.
%    sizeisnew   1x2 logical array; activates reshaping of A and B.
%    idA/B       May change only if squashOK ~= 0
%    squashOK    If only A or B is a multi-block array (M-B) and the other
%                is single-block (1-B), it will be rearranged from N-D to
%                2-D. If both A and B are 1-B or M-B arrays, squashOK = 0.
%                If only A (or B) is a M-B array, squashOK = 1 (or 2).
%    sxtimesOK, timesOK, mtimesOK    Flags controlling MAIN 2 (TRUE/FALSE).
%    sumOK       Dimension along which SUM is performed. If SUM is not
%                needed, sumOK = 0.

% Initializing output arguments

    idA = idA0;
    idB = idB0;
     squashOK = 0;
    sxtimesOK = false;
      timesOK = false;
     mtimesOK = false;
        sumOK = 0;
    shiftC = 0;
    delC = 0;

% Checking for gross input errors

    NidA = numel(idA);
    NidB = numel(idB);
    idA1 = idA(1);
    idB1 = idB(1);
    if  NidA>2 || NidB>2 || NidA==0 || NidB==0 || ...
           ~isreal(idA1) ||    ~isreal(idB1)   || ...
        ~isnumeric(idA1) || ~isnumeric(idB1)   || ...
                 0>idA1  ||          0>idB1    || ... % negative
         idA1~=fix(idA1) ||  idB1~=fix(idB1)   || ... % non-integer
         ~isfinite(idA1) ||  ~isfinite(idB1) % Inf or NaN
        error('MULTIPROD:InvalidDimensionArgument', ...
        ['Internal-dimension arguments (e.g., [IDA1 IDA2]) must\n', ...
         'contain only one or two non-negative finite integers']);
    end

% Checking Syntaxes containing zeros (4b/c)

    declared_outer = false;
    idA2 = idA(NidA); % It may be IDA1 = IDA2 (1-D block)
    idB2 = idB(NidB);

    if any(idA==0) || any(idB==0)

        % "Inner products": C = MULTIPROD(A, B, [0 DA1], [DB1 0])
        if idA1==0 && idA2>0 && idB1>0 && idB2==0
            idA1 = idA2;
            idB2 = idB1;
        % "Outer products": C = MULTIPROD(A, B, [DA1 0], [0 DB1])
        elseif idA1>0 && idA2==0 && idB1==0 && idB2>0
            declared_outer = true;
            idA2 = idA1;
            idB1 = idB2;
        else
            error('MULTIPROD:InvalidDimensionArgument', ...
            ['Misused zeros in the internal-dimension arguments\n', ...
            '(see help heads 4b and 4c)']);
        end
        NidA = 1;
        NidB = 1;
        idA = idA1;
        idB = idB1;

    elseif (NidA==2 && idA2~=idA1+1) || ...  % Non-adjacent IDs
           (NidB==2 && idB2~=idB1+1)
        error('MULTIPROD:InvalidDimensionArgument', ...
        ['If an array contains 2-D blocks, its two internal dimensions', ...
        'must be adjacent (e.g. IDA2 == IDA1+1)']);
    end

% ESC - Case for which no reshaping is needed (both A and B are scalars)

    scalarA = isequal(sizeA0, [1 1]);
    scalarB = isequal(sizeB0, [1 1]);
    if scalarA && scalarB
        sizeA = sizeA0;
        sizeB = sizeB0;
        sizeisnew = [false false];
        timesOK = true; return
    end

% Computing and checking adjusted sizes
% The lengths of ADJSIZEA and ADJSIZEB must be >= IDA(END) and IDB(END)

    NsA = idA2 - length(sizeA0); % Number of added trailing singletons
    NsB = idB2 - length(sizeB0);
    adjsizeA = [sizeA0 ones(1,NsA)];
    adjsizeB = [sizeB0 ones(1,NsB)];
    extsizeA = adjsizeA([1:idA1-1, idA2+1:end]); % Size of EDs
    extsizeB = adjsizeB([1:idB1-1, idB2+1:end]);
    p = adjsizeA(idA1);
    q = adjsizeA(idA2);
    r = adjsizeB(idB1);
    s = adjsizeB(idB2);
    scalarsinA = (p==1 && q==1);
    scalarsinB = (r==1 && s==1);
    singleA = all(extsizeA==1);
    singleB = all(extsizeB==1);
    if q~=r && ~scalarsinA && ~scalarsinB && ~declared_outer
       error('MULTIPROD:InnerDimensionsMismatch', ...
             'Inner matrix dimensions must agree.');
    end

% STEP 1/3 - DIMENSION SHIFTING (FIRST STEP OF AX)
%   Pipeline 1 (using TIMES) never needs left, and may need right shifting.
%   Pipeline 2 (using MTIMES) may need left shifting of A and right of B.

    shiftA = 0;
    shiftB = 0;
    diffBA = idB1 - idA1;
    if scalarA % Do nothing
    elseif singleA && ~scalarsinB, shiftA = -idA1 + 1; %  Left shifting A
    elseif idB1 > idA1,            shiftA = diffBA;    % Right shifting A
    end
    if scalarB % Do nothing
    elseif singleB && ~scalarsinA, shiftB = -idB1 + 1; %  Left shifting B
    elseif idA1 > idB1,            shiftB = -diffBA;   % Right shifting B
    end

% STEP 2/3 - SELECTION OF PROPER ENGINE AND BLOCK SIZE ADJUSTMENTS

    addA  = 0; addB  = 0;
    delA  = 0; delB  = 0;
    swapA = 0; swapB = 0;
    idC1 = max(idA1, idB1);
    idC2 = idC1 + 1;
    checktimes = false;

    if (singleA||singleB) &&~scalarsinA &&~scalarsinB % Engine using MTIMES

        if singleA && singleB
            mtimesOK = true;
            shiftC=idC1-1; % Right shifting C
            idC1=1; idC2=2;
        elseif singleA
            squashOK = 2;
            idB = [idB1, idB1+1] + shiftB;
        else % singleB
            squashOK = 1;
            idA = [idA1, idA1+1] + shiftA;
        end

        if NidA==2 && NidB==2 % 1) 2-D BLOCKS BY 2-D BLOCKS
            % OK
        elseif NidA==2        % 2) 2-D BLOCKS BY 1-D BLOCKS
            addB=idB1+1; delC=idC2;
        elseif NidB==2        % 3) 1-D BLOCKS BY 2-D BLOCKS
            addA=idA1; delC=idC1;
        else                  % 4) 1-D BLOCKS BY 1-D BLOCKS
            if declared_outer
                addA=idA1+1; addB=idB1;
            else
                addA=idA1; addB=idB1+1; delC=idC2;
            end
        end

    else % Engine using TIMES (also used if SCALARA || SCALARB)

        sxtimesOK = true;

        if NidA==2 && NidB==2 % 1) 2-D BLOCKS BY 2-D BLOCKS

            if scalarA || scalarB
                timesOK=true;
            elseif scalarsinA && scalarsinB % scal-by-scal
                checktimes=true;
            elseif scalarsinA || scalarsinB || ... % scal-by-mat
                (q==1 && r==1)  % vec-by-vec ("outer")
            elseif p==1 && s==1 % vec-by-vec ("inner")
                swapA=idA1; sumOK=idC1; checktimes=true;
            elseif s==1 % mat-by-vec
                swapB=idB1; sumOK=idC2;
            elseif p==1 % vec-by-mat
                swapA=idA1; sumOK=idC1;
            else % mat-by-mat
                addA=idA2+1; addB=idB1; sumOK=idC2; delC=idC2;
            end

        elseif NidA==2 % 2) 2-D BLOCKS BY 1-D BLOCKS

            if scalarA || scalarB
                timesOK=true;
            elseif scalarsinA && scalarsinB % scal-by-scal
                addB=idB1; checktimes=true;
            elseif scalarsinA % scal-by-vec
                delA=idA1;
            elseif scalarsinB % mat-by-scal
                addB=idB1;
            elseif p==1 % vec-by-vec ("inner")
                delA=idA1; sumOK=idC1; checktimes=true;
            else % mat-by-vec
                addB=idB1; sumOK=idC2; delC=idC2;
            end

        elseif NidB==2 % 3) 1-D BLOCKS BY 2-D BLOCKS

            if scalarA || scalarB
                timesOK=true;
            elseif scalarsinA && scalarsinB % scal-by-scal
                addA=idA1+1; checktimes=true;
            elseif scalarsinB % vec-by-scal
                delB=idB2;
            elseif scalarsinA % scal-by-mat
                addA=idA1+1;
            elseif s==1 % vec-by-vec ("inner")
                delB=idB2; sumOK=idC1; checktimes=true;
            else % vec-by-mat
                addA=idA1+1; sumOK=idC1; delC=idC1;
            end

        else % 4) 1-D BLOCKS BY 1-D BLOCKS

            if scalarA || scalarB
                timesOK=true;
            elseif declared_outer % vec-by-vec ("outer")
                addA=idA1+1; addB=idB1;
            elseif scalarsinA && scalarsinB % scal-by-scal
                checktimes=true;
            elseif scalarsinA || scalarsinB % vec-by-scal
            else % vec-by-vec
                sumOK=idC1; checktimes=true;
            end
        end
    end

% STEP 3/3 - Adjusting the size of A and B. The size of C is adjusted
%            later, because it is not known yet.

    [sizeA, sizeisnew(1)] = adjustsize(sizeA0, shiftA, addA, delA, swapA);
    [sizeB, sizeisnew(2)] = adjustsize(sizeB0, shiftB, addB, delB, swapB);

    if checktimes % Faster than calling BBXFUN
        diff = length(sizeB) - length(sizeA);
        if isequal([sizeA ones(1,diff)], [sizeB ones(1,-diff)])
            timesOK = true;
        end
    end

function [sizeA, sizeisnew] = adjustsize(sizeA0, shiftA, addA, delA, swapA)
% ADJUSTSIZE  Adjusting size of a block array.

    % Dimension shifting (by adding or deleting trailing singleton dim.)
    if     shiftA>0, [sizeA,newA1] = addsing(sizeA0, 1, shiftA);
    elseif shiftA<0, [sizeA,newA1] = delsing(sizeA0, 1,-shiftA);
    else   sizeA = sizeA0;  newA1  = false;
    end
    % Modifying block size (by adding, deleting, or moving singleton dim.)
    if      addA, [sizeA,newA2] = addsing(sizeA, addA+shiftA, 1); % 1D-->2D
    elseif  delA, [sizeA,newA2] = delsing(sizeA, delA+shiftA, 1); % 2D-->1D
    elseif swapA, [sizeA,newA2] = swapdim(sizeA,swapA+shiftA); % ID Swapping
    else                 newA2  = false;
    end
    sizeisnew = newA1 || newA2;

function [newsize, flag] = addsing(size0, dim, ns)
%ADDSING   Adding NS singleton dimensions to the size of an array.
%   Warning: NS is assumed to be a positive integer.
%   Example: If the size of A is ..... SIZE0 = [5 9 3]
%            NEWSIZE = ADDSING(SIZE0, 3, 2) is [5 9 1 1 3]

    if dim > length(size0)
        newsize = size0;
        flag = false;
    else
        newsize = [size0(1:dim-1), ones(1,ns), size0(dim:end)];
        flag = true;
    end

function [newsize, flag] = delsing(size0, dim, ns)
%DELSING   Removing NS singleton dimensions from the size of an array.
%   Warning: Trailing singletons are not removed
%   Example: If the size of A is SIZE0 = [1 1 1 5 9 3]
%            NEWSIZE = DELSING(SIZE, 1, 3) is  [5 9 3]

    if dim > length(size0)-ns % Trailing singletons are not removed
        newsize = size0;
        flag = false;
    else % Trailing singl. added, so NEWSIZE is guaranteed to be 2D or more
        newsize = size0([1:dim-1, dim+ns:end, dim]);
        flag = true;
    end

function [newsize, flag] = swapdim(size0, dim)
%SWAPDIM   Swapping two adjacent dimensions of an array (DIM and DIM+1).
%   Used only when both A and B are multi-block arrays with 2-D blocks.
%   Example: If the size of A is .......... 5?(6?3)
%            NEWSIZE = SWAPIDS(SIZE0, 2) is 5?(3?6)

    newsize = [size0 1]; % Guarantees that dimension DIM+1 exists.
    newsize = newsize([1:dim-1, dim+1, dim, dim+2:end]);
    flag = true;

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function [cellmap] = createObjMap(inputImages,varargin)
	% Creates a cellmap from a ZxXxY input matrix of input images.
	% Biafra Ahanonu
	% started: 2013.10.12
	% inputs
		%
	% outputs
		%

	[cellmap] = ciapkg.image.createObjMap(inputImages,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function rd = jed_to_rd ( jed )

%*****************************************************************************80
%
%% JED_TO_RD converts a JED to an RD.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Reingold, Nachum Dershowitz,
%    Calendrical Calculations, the Millennium Edition,
%    Cambridge, 2002.
%
%  Parameters:
%
%    Input, real JED, the Julian Ephemeris Date.
%
%    Output, real RD, the RD date.
%
  rd_epoch = epoch_to_jed_rd ( );

  rd = jed - rd_epoch;

  return
end


text:
function distance_map_in_meters = distance_in_meters_cityscapes(...
    depth_map_in_meters, camera_parameters_file)
%DISTANCE_IN_METERS_CITYSCAPES  Compute air thickness, i.e. distance of depicted
%object from camera center, expressed in meters, as a dense map with same
%resolution as the image, using a dense depth map and intrinsic camera
%parameters as inputs.

% Retrieve relevant intrinsic camera parameters: focal length and optical
% center, both expressed in pixel coordinates.
[~, f_x, c_x, c_y] = camera_parameters_cityscapes(camera_parameters_file);

% Compute medium (i.e. air) thickness in meters from depth. The derivation of
% the formula in the final lines of code is based on similar triangles.

[height, width] = size(depth_map_in_meters);
[X, Y] = meshgrid(1:width, 1:height);

distance_map_in_meters = depth_map_in_meters .*...
    sqrt((f_x ^ 2 + (X - c_x) .^ 2 + (Y - c_y) .^ 2) / f_x ^ 2);

end

author: sakaridis


text:
function g = whiteblockKernGradient(kern, x, x2, covGrad)

% WHITEBLOCKKERNGRADIENT Gradient of WHITEBLOCK kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the
% white noise block kernel's parameters. As well as the kernel structure
% and the input positions, the user provides a matrix PARTIAL which gives
% the partial derivatives of the function with respect to the
% relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x : the input locations for which the gradients are being
% computed.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in X.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x1 : the input locations associated with the rows of the
% kernel matrix.
% ARG x2 : the input locations associated with the columns of the
% kernel matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as X1 and the same number of columns
% as X2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% SEEALSO whiteblockKernParamInit, whiteblockKernDiagGradient
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if nargin < 4
    covGrad = x2;
end

g = zeros(1, kern.nout);
if nargin < 4
    startOne = 1;
    endOne = 0;
    for i=1:kern.nout
        if iscell(x)
            endOne = endOne + size(x{i},1);
        else
            endOne = endOne + size(x,1);
        end
        g(1, i) = trace(covGrad(startOne:endOne, startOne:endOne));
        startOne = endOne + 1;
    end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function pass = test_Poisson( )
% Check correctness of Poisson solver on the sphere:

% Discretization sizes:
m = 40;
n = 40;
tol = 1e3*chebfunpref().cheb2Prefs.chebfun2eps;

% Example 1:
f = spherefun(@(lam,th) -6*cos(lam).*cos(th).*sin(th));
exact = spherefun(@(lam,th) sin(th).*cos(th).*cos(lam));
u = spherefun.poisson(f, 0, m, n);
pass(1) = ( norm(u - exact, inf) < tol );

% Example 2:
f = spherefun(@(lam,th) -4*(3*cos(th)+5*cos(3*th)).*sin(lam).*sin(th));
exact = spherefun(@(lam,th) -2*sin(lam).*sin(2*th).*sin(th).^2);
u = spherefun.poisson(f, 0, m, n);
pass(2) = ( norm(u - exact, inf) < tol );

% Inline example:
f = spherefun(@(lam,th) -6*(-1+5*cos(2*th)).*sin(lam).*sin(2*th));
exact = spherefun(@(lam,th) -2*sin(lam).*sin(2*th).*sin(th).^2 -...
            sin(lam).*sin(th).*cos(th) + .5*sin(lam).*sin(2*th).*cos(2*th));
u = spherefun.poisson(f, 0, m, n);
pass(3) = ( norm(u - exact, inf) < tol );

% Check that the mean is zero.
pass(4) = ( abs(mean2(u)) < tol );

% Check that the code properly deals with a right hand side without a
% mean of zero.
f = spherefun(@(x,y,z) 1 + x);
warning('off','CHEBFUN:SPHEREFUN:POISSON:meanRHS');
u = spherefun.poisson(f, 0, 10);
warning('on','CHEBFUN:SPHEREFUN:POISSON:meanRHS');
% This f - mean2(f)
g = spherefun(@(x,y,z) x);
% Solution with the mean of f set to zero
v = spherefun.poisson(g, 0, 10);
pass(5) = ( norm(u - v) < tol );

% Check that the code allows the mean of the solution to be set to
% something other than zero.
f = spherefun(@(x,y,z) x.*y.*z );
u = spherefun.poisson(f, 1, 10);
pass(6) = ( abs(mean2(u)-1) < tol );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [stop, reason] = stoppingcriterion(problem, x, options, info, last)
% Checks for standard stopping criteria, as a helper to solvers.
%
% function [stop, reason] = stoppingcriterion(problem, x, options, info, last)
%
% Executes standard stopping criterion checks, based on what is defined in
% the info(last) stats structure and in the options structure.
%
% The returned number 'stop' is 0 if none of the stopping criteria
% triggered, and a (strictly) positive integer otherwise. The integer
% identifies which criterion triggered:
%  0 : Nothing triggered;
%  1 : Cost tolerance reached;
%  2 : Gradient norm tolerance reached;
%  3 : Max time exceeded;
%  4 : Max iteration count reached;
%  5 : Maximum number of cost evaluations reached;
%  6 : User defined stopfun criterion triggered.
%
% The output 'reason' is a string describing the triggered event.

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%
%   April 2, 2015 (NB):
%       'reason' now contains the option (name and value) that triggered.

    stop = 0;
    reason = '';

    stats = info(last);

    % Target cost attained
    if isfield(stats, 'cost') && isfield(options, 'tolcost') && ...
       stats.cost <= options.tolcost
        reason = sprintf('Cost tolerance reached; options.tolcost = %g.', options.tolcost);
        stop = 1;
        return;
    end

    % Target gradient norm attained
    if isfield(stats, 'gradnorm') && isfield(options, 'tolgradnorm') && ...
       stats.gradnorm < options.tolgradnorm
        reason = sprintf('Gradient norm tolerance reached; options.tolgradnorm = %g.', options.tolgradnorm);
        stop = 2;
        return;
    end

    % Allotted time exceeded
    if isfield(stats, 'time') && isfield(options, 'maxtime') && ...
       stats.time >= options.maxtime
        reason = sprintf('Max time exceeded; options.maxtime = %g.', options.maxtime);
        stop = 3;
        return;
    end

    % Allotted iteration count exceeded
    if isfield(stats, 'iter') && isfield(options, 'maxiter') && ...
       stats.iter >= options.maxiter
        reason = sprintf('Max iteration count reached; options.maxiter = %g.', options.maxiter);
        stop = 4;
        return;
    end

    % Allotted function evaluation count exceeded
    if isfield(stats, 'costevals') && isfield(options, 'maxcostevals') && ...
       stats.costevals >= options.maxcostevals
        reason = sprintf('Maximum number of cost evaluations reached; options.maxcostevals = %g.', options.maxcostevals);
        stop = 5;
    end

    % Check whether the possibly user defined stopping criterion
    % triggers or not.
    if isfield(options, 'stopfun')
        userstop = options.stopfun(problem, x, info, last);
        if userstop
            reason = 'User defined stopfun criterion triggered; see options.stopfun.';
            stop = 6;
            return;
        end
    end

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function plotdecisionboundary(Theta, X, Y)

plotdata(X(:,2:3), Y);

hold on;

if size(X, 2) <= 3

    plot_X = [min(X(:,2))-2,  max(X(:,2))+2];

    plot_Y = (-1./Theta(3)).*(Theta(2).*plot_X + Theta(1));

    plot(plot_X, plot_Y);

    legend('Admitted', 'Not admitted', 'Decision Boundary');

    axis([30, 100, 30, 100]);

else

    U = linspace(-1, 1.5, 50);

    V = linspace(-1, 1.5, 50);

    Z = zeros(length(U), length(V));

    for i = 1:length(U)

        for j = 1:length(V)

            Z(i,j) = mapFeature(U(i), V(j))*Theta;

        end

    end

    Z = Z';

end

hold off;

end

author: TheAlgorithms

repo: MATLAB-Octave


text:
function [xm, ym, xfit,yfit] = companal(icapix, rs)

sica = size(icapix,3);
fprintf('maxima: \n')
for ii=1:sica
    aim = abs(imresize(icapix(:,:,ii),1/rs));
    m(ii) = max(max(aim));
    [ym(ii), xm(ii)] = find(m(ii)==aim);
    [ycog(ii), xcog(ii)] = cog(imresize(icapix(:,:,ii),1/rs)); %center of gravity
    [xfit(ii), yfit(ii), sig(ii)] = fitgauss2d(imresize(icapix(:,:,ii),1/rs));
    fprintf('%g . maximum: %g %g \tcog: %g %g \tfit: %g %g \tsig: %g\n',ii , [xm(ii), ym(ii)], [xcog(ii) ycog(ii)], [xfit(ii), yfit(ii), sig(ii)] );
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB