text:
function res = saveastiff(data, path, options)
% options.color
%   : true or FALSE
%   : If this is true, third dimension should be 3 and the data is saved as a color image.
% options.compress
%   : 'no', 'lzw', 'jpeg' or 'adobe'.
%     Compression type.
%       'no'    : Uncompressed(Default)
%       'lzw'   : lossless LZW
%       'jpeg'  : lossy JPEG (When using JPEG compression, ImageWidth,
%                 ImageLength, and RowsPerStrip must be multiples of 16.)
%       'adobe' : lossless Adobe-style
% options.message
%   : TRUE or false.
%     If this is false, all messages are skipped.
% options.append
%   : true or FALSE
%     If path is exist, the data is appended to an existing file.
%     If path is not exist, this options is ignored.
% options.overwrite
%   : true or FALSE
%     Overwrite to an existing file.
% options.big
%   : true or FALSE,
%     Use 64 bit addressing and allows for files > 4GB
%
% Defalut value of 'options' is
%     options.color     = false;
%     options.compress  = 'no';
%     options.message   = true;
%     options.append    = false;
%     options.overwrite = false;
%     options.big       = false;
%
% res : Return value. It is 0 when the function is finished with no error.
%       If an error is occured in the function, it will have a positive
%       number (error code).
%
% Copyright (c) 2012, YoonOh Tak
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of the Gwangju Institute of Science and Technology (GIST), Republic of Korea nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

tStart = tic;
errcode = 0;
try
%% Init options parameter
if nargin < 3 % Use default options
    options.color = false;
    options.compress = 'no';
    options.message = true;
    options.append = false;
    options.overwrite = false;
end
if ~isfield(options, 'message'),   options.message   = true; end
if ~isfield(options, 'append'),    options.append    = false; end
if ~isfield(options, 'compress'),  options.compress  = 'no';  end
if ~isfield(options, 'color'),     options.color     = false; end
if ~isfield(options, 'overwrite'), options.overwrite = false; end
if  isfield(options, 'big') == 0,  options.big       = false; end

if isempty(data), errcode = 1; assert(false); end
if (options.color == false && ndims(data) > 3) || ...
   (options.color == true && ndims(data) > 4)
    % Maximum dimension of a grayscale image is 3 of [height, width, frame]
    % Maximum dimension of a color image is 4 of [height, width, color, frame]
    errcode = 2; assert(false);
end

%% Get image informations
% http://www.awaresystems.be/imaging/tiff/tifftags/photometricinterpretation.html
if ~options.color
    if ndims(data) >= 4, errcode = 2; assert(false); end;
    [height, width, depth] = size(data);
    tagstruct.Photometric = Tiff.Photometric.MinIsBlack;
%     tagstruct.Photometric = Tiff.Photometric.MinIsWhite;
%     tagstruct.Photometric = Tiff.Photometric.Mask;
%     tagstruct.Photometric = Tiff.Photometric.Separated;
else
    if ndims(data) >= 5, errcode = 2; assert(false); end;
    [height, width, cc, depth] = size(data); % cc: color channels. 3: rgb, 4: rgb with alpha channel
    if cc ~= 3 && cc ~= 4, errcode = 3; assert(false); end;
    tagstruct.Photometric = Tiff.Photometric.RGB;
%     tagstruct.Photometric = Tiff.Photometric.CIELab;
%     tagstruct.Photometric = Tiff.Photometric.ICCLab;
%     tagstruct.Photometric = Tiff.Photometric.ITULab;
%     (Unsupported)tagstruct.Photometric = Tiff.Photometric.Palette;
%     (Unsupported)tagstruct.Photometric = Tiff.Photometric.YCbCr;
end
tagstruct.ImageLength = height;
tagstruct.ImageWidth = width;
tagstruct.PlanarConfiguration = Tiff.PlanarConfiguration.Chunky; % (RGB RGB,RGB RGB,RGB RGB), http://www.awaresystems.be/imaging/tiff/tifftags/planarconfiguration.html
% tagstruct.PlanarConfiguration = Tiff.PlanarConfiguration.Separate; % (RRR RRR, GGG GGG, BBB BBB), % http://www.awaresystems.be/imaging/tiff/tifftags/planarconfiguration.html

%% Complex number
% http://www.awaresystems.be/imaging/tiff/tifftags/samplesperpixel.html
if ~options.color && isreal(data) % Grayscale image with real numbers
    tagstruct.SamplesPerPixel = 1;
    data = reshape(data, height, width, 1, depth);
elseif ~options.color && ~isreal(data) % Grayscale image with complex numbers
    tagstruct.SamplesPerPixel = 2;
    data = reshape([real(data) imag(data)], height, width, 2, depth);
elseif options.color && isreal(data) % Color image with real numbers
    tagstruct.SamplesPerPixel = cc;
    if cc == 4
        tagstruct.ExtraSamples = Tiff.ExtraSamples.AssociatedAlpha; % The forth channel is alpha channel
    end
    data = reshape(data, height, width, cc, depth);
elseif options.color && ~isreal(data) % Color image with complex numbers
    tagstruct.SamplesPerPixel = cc * 2;
    if cc == 3
        tagstruct.ExtraSamples = repmat(Tiff.ExtraSamples.Unspecified, 1, 3); % 3(real)+3(imag) = 6 = 3(rgb) + 3(Extra)
    else
        tagstruct.ExtraSamples = repmat(Tiff.ExtraSamples.Unspecified, 1, 5); % 4(real)+4(imag) = 8 = 3(rgb) + 5(Extra)
    end
    data = reshape([real(data) imag(data)], height, width, cc*2, depth);
end

%% Image compression
% http://www.awaresystems.be/imaging/tiff/tifftags/compression.html
switch lower(options.compress)
    case 'no'
        tagstruct.Compression = Tiff.Compression.None;
    case 'lzw'
        tagstruct.Compression = Tiff.Compression.LZW;
    case 'jpeg'
        tagstruct.Compression = Tiff.Compression.JPEG;
    case 'adobe'
        tagstruct.Compression = Tiff.Compression.AdobeDeflate;
    otherwise
        % Use tag nubmer in http://www.awaresystems.be/imaging/tiff/tifftags/compression.html
        tagstruct.Compression = options.compress;
end

%% Sample format
% http://www.awaresystems.be/imaging/tiff/tifftags/sampleformat.html
switch class(data)
    % Unsupported Matlab data type: char, logical, cell, struct, function_handle, class.
    case {'uint8', 'uint16', 'uint32'}
        tagstruct.SampleFormat = Tiff.SampleFormat.UInt;
    case {'int8', 'int16', 'int32'}
        tagstruct.SampleFormat = Tiff.SampleFormat.Int;
        if options.color
            errcode = 4; assert(false);
        end
    case {'single', 'double', 'uint64', 'int64'}
        tagstruct.SampleFormat = Tiff.SampleFormat.IEEEFP;
    otherwise
        % (Unsupported)Void, ComplexInt, ComplexIEEEFP
        errcode = 5; assert(false);
end

%% Bits per sample
% http://www.awaresystems.be/imaging/tiff/tifftags/bitspersample.html
switch class(data)
    case {'uint8', 'int8'}
        tagstruct.BitsPerSample = 8;
    case {'uint16', 'int16'}
        tagstruct.BitsPerSample = 16;
    case {'uint32', 'int32'}
        tagstruct.BitsPerSample = 32;
    case {'single'}
        tagstruct.BitsPerSample = 32;
    case {'double', 'uint64', 'int64'}
        tagstruct.BitsPerSample = 64;
    otherwise
        errcode = 5; assert(false);
end

%% Rows per strip
maxstripsize = 8*1024;
tagstruct.RowsPerStrip = ceil(maxstripsize/(width*(tagstruct.BitsPerSample/8)*size(data,3))); % http://www.awaresystems.be/imaging/tiff/tifftags/rowsperstrip.html
if tagstruct.Compression == Tiff.Compression.JPEG
    tagstruct.RowsPerStrip = max(16,round(tagstruct.RowsPerStrip/16)*16);
end

%% Overwrite check
if exist(path, 'file') && ~options.append
    if ~options.overwrite
        errcode = 6; assert(false);
    end
end

%% Save path configuration
path_parent = pwd;
[pathstr, fname, fext] = fileparts(path);
if ~isempty(pathstr)
    if ~exist(pathstr, 'dir')
        mkdir(pathstr);
    end
    cd(pathstr);
end

%% Write image data to a file
file_opening_error_count = 0;
while ~exist('tfile', 'var')
    try
        if ~options.append % Make a new file
            s=whos('data');
            if s.bytes > 2^32-1 || options.big
                tfile = Tiff([fname, fext], 'w8'); % Big Tiff file
            else
                tfile = Tiff([fname, fext], 'w');
            end
        else
            if ~exist([fname, fext], 'file') % Make a new file
                s=whos('data');
                if s.bytes > 2^32-1 || options.big
                    tfile = Tiff([fname, fext], 'w8'); % Big Tiff file
                else
                    tfile = Tiff([fname, fext], 'w');
                end
            else % Append to an existing file
                tfile = Tiff([fname, fext], 'r+');
                while ~tfile.lastDirectory(); % Append a new image to the last directory of an exiting file
                    tfile.nextDirectory();
                end
                tfile.writeDirectory();
            end
        end
    catch
        file_opening_error_count = file_opening_error_count + 1;
        pause(0.1);
        if file_opening_error_count > 5 % automatically retry to open for 5 times.
            reply = input('Failed to open the file. Do you wish to retry? Y/n: ', 's');
            if isempty(reply) || any(upper(reply) == 'Y')
                file_opening_error_count = 0;
            else
                errcode = 7;
                assert(false);
            end
        end
    end
end

for d = 1:depth
    tfile.setTag(tagstruct);
    tfile.write(data(:, :, :, d));
    if d ~= depth
       tfile.writeDirectory();
    end
end

tfile.close();
if exist('path_parent', 'var'), cd(path_parent); end

tElapsed = toc(tStart);
if options.message
    display(sprintf('The file was saved successfully. Elapsed time : %.3f s.', tElapsed));
end

catch exception
%% Exception management
    if exist('tfile', 'var'), tfile.close(); end
    switch errcode
        case 1
            if options.message, error '''data'' is empty.'; end;
        case 2
            if options.message, error 'Data dimension is too large.'; end;
        case 3
            if options.message, error 'Third dimesion (color depth) should be 3 or 4.'; end;
        case 4
            if options.message, error 'Color image cannot have int8, int16 or int32 format.'; end;
        case 5
            if options.message, error 'Unsupported Matlab data type. (char, logical, cell, struct, function_handle, class)'; end;
        case 6
            if options.message, error 'File already exists.'; end;
        case 7
            if options.message, error(['Failed to open the file ''' path '''.']); end;
        otherwise
            if exist('fname', 'var') && exist('fext', 'var')
                delete([fname fext]);
            end
            if exist('path_parent', 'var'), cd(path_parent); end
            rethrow(exception);
    end
    if exist('path_parent', 'var'), cd(path_parent); end
end
res = errcode;
end

author: flatironinstitute

repo: NoRMCorre


text:
function panel = BstPanel( varargin )
% Constructor for object BstPanel.
% A BstPanel object holds mainly a java Swing container to be displayed in the
% Brainstorm main window.
%
% Constructor call :
%     BstPanel(name, jHandle, sControls)
%     BstPanel() : just to have a data template
%
% Data structure
%   - jHandle    : java handle to a javax.swing.JComponent object (panel handle)
%   - name       : string identifier to index and display this panel
%   - container  : string identifier of the panel container in which it should be displayed
%   - sControls  : list of the useful controls of the panel

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2008-2010

% Parse inputs
if (nargin == 0)
    % Create an empty object (template to create structures)
    panel.jHandle    = [];
    panel.name       = '';
    panel.container  = struct('type',   '', ...
                              'handle', {});
    panel.sControls  = [];
    panel = class(panel, 'BstPanel');
    panel = repmat(panel, 0);
    return;
elseif (nargin == 3)
    name       = varargin{1};
    jHandle    = varargin{2};
    sControls  = varargin{3};
else
    error('Usage : BstPanel(name, jHandle, sControls)');
end

% Check inputs
if (~isa(jHandle, 'javax.swing.JComponent'))
    error('First argument ''jHandle'' must be a javax.swing.JComponent object');
elseif (~ischar(name))
    error('Second argument ''name'' must be a Matlab string');
elseif (~isstruct(sControls) && ~isempty(sControls))
    error('Third argument ''sControls'' must be a structure of Java Swing controls');
end
% Check intergrity of the structure sControls
if (isempty(sControls))
    sControls = struct();
end

% Create data structure
panel.jHandle          = jHandle;
panel.name             = name;
panel.container.type   = '';
panel.container.handle = {};
panel.sControls        = sControls;
% Set object class
panel = class(panel, 'BstPanel');

author: brainstorm-tools


text:
function c = tapas_squared_pe_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuraton for the squared-prediction error optimization of perceptual parameters
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The squared-prediction error optimization infers the perceptual parameter values that lead to the
% best predictions of input according to the (conditional) criterion of the sum of squared
% errors. The criterion is conditional in the sense that the priors on the perceptual parameters
% retain a certain weight determined by the parameter zeta, whose prior is defined below. Zeta can
% be interpreted as an inverse weight on the prediction errors: greater values of zeta lead to less
% influence of prediction errors as opposed to priors. One would usually leave zeta fixed.
%
% Usage:
%     tapas_fitModel([], inputs, '<perceptual_model>', 'tapas_squared_pe_config', ...)
%
% Note that the first argument (responses) is empty.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Model name
c.model = 'tapas_gaussian_obs';

% Sufficient statistics of Gaussian parameter priors
%
% Zeta
c.logzemu = log(0.05);
c.logzesa = 0;

% Gather prior settings in vectors
c.priormus = [
    c.logzemu,...
         ];

c.priorsas = [
    c.logzesa,...
         ];

% Model filehandle
c.obs_fun = @tapas_squared_pe;

% Handle to function that transforms observation parameters to their native space
% from the space they are estimated in
c.transp_obs_fun = @tapas_squared_pe_transp;

return;

author: translationalneuromodeling

repo: tapas


text:
function dimNumber = getDimension(MRSIStruct, dimLabel)
    arguments
        MRSIStruct (1,1) struct
        dimLabel char
    end
    switch(lower(dimLabel))
        case 'x'
            dimNumber = MRSIStruct.dims.x;
        case 'y'
            dimNumber = MRSIStruct.dims.y;
        case 'z'
            dimNumber = MRSIStruct.dims.z;
        case 'coils'
            dimNumber = MRSIStruct.dims.coils;
        case 'averages'
            dimNumber = MRSIStruct.dims.averages;
        case 'average'
            dimNumber = MRSIStruct.dims.averages;
        case 'kx'
            dimNumber = MRSIStruct.dims.kx;
        case 'ky'
            dimNumber = MRSIStruct.dims.ky;
        case 'kz'
            dimNumber = MRSIStruct.dims.kz;
        case 'subspecs'
            dimNumber = MRSIStruct.dims.subSpecs;
        case 'subspec'
            dimNumber = MRSIStruct.dims.subSpecs;
        case 'extras'
            dimNumber = MRSIStruct.dims.extras;
        case 'extra'
            dimNumber = MRSIStruct.dims.extras;
        case 't'
            dimNumber = MRSIStruct.dims.t;
        otherwise
            error('Dimension not found %s', dimLabel);
    end
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function F = restrict(F, newDomain)
%RESTRICT   Restrict a CHEBFUN object to a subinterval.
%   G = RESTRICT(F, [S1, S2]) returns a CHEBFUN G defined on the interval [S1,
%   S2] which agrees with F on that interval. Any interior breakpoints in
%   F.DOMAIN within [S1, S2] are kept in G.DOMAIN.
%
%   G = RESTRICT(F, S), where S is a row vector, will introduce additional
%   interior breakpoints at S(2:end-2).
%
%   In both cases, if S(1) > S(end), S(1) < F.domain(1), or S(end) >
%   F.domain(end), then an error is returned. If S is empty or a scalar, then an
%   empty CHEBFUN G is returned.
%
%   Note that G will not be 'simplified'. If this is required, call G =
%   SIMPLIFY(RESTRICT(F)), or G = F{S}.
%
% See also OVERLAP, SUBSREF, DEFINE, SIMPLIFY.

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

% Tweak the domain of a quasimatrix input:
[F, newDomain] = tweakDomain(F, newDomain);

% Loop over the columns:
for k = 1:numel(F)
    F(k) = columnRestrict(F(k), newDomain);
end

end

function f = columnRestrict(f, newDomain)

% Ignore duplicate entries in newDomain:
newDomain([false, ~diff(newDomain)]) = [];

% Empty case:
if ( isempty(f) )
    return
elseif ( isempty(newDomain) || numel(newDomain) == 1 )
    f = chebfun();
    return
end

% Cast a periodic chebfun to a regualr chebfun:
if ( isPeriodicTech(f) )
    f = chebfun(f);
end

% Grab domain from f:
oldDomain = f.domain;

% Trivial case:
if ( (numel(newDomain) == 2) && domainCheck(f, newDomain([1, end])) )
    % The domains are the same!
    return
end

% Check the domain is valid:
if ( (newDomain(1) < oldDomain(1)) || (newDomain(end) > oldDomain(end)) || ...
        any(diff(newDomain) < 0) )
    % newDom is not a valid subinterval of oldDom!
    error('CHEBFUN:CHEBFUN:restrict:subdom', 'Not a valid subdomain.');
end

% Obtain FUN cell and pointValues from f:
funs = f.funs;
pointValues = f.pointValues;

% Discard intervals to the left:
discardIntsLeft = oldDomain(2:end) < newDomain(1);
oldDomain([discardIntsLeft, false]) = [];
funs(discardIntsLeft) = [];
pointValues([discardIntsLeft, false],:,:) = [];

% Discard intervals to the right:
discardIntsRright = oldDomain(1:end-1) > newDomain(end);
oldDomain([false, discardIntsRright]) = [];
funs(discardIntsRright) = [];
pointValues([false, discardIntsRright],:,:) = [];

% Take the union of the new and old domains:
if ( ~isempty( oldDomain(2:end-1) ) ) % Required due to Matlab union() behavior.
    % NB: the old endpoints will either be in newDomain or are not required.
    newDomain = union(oldDomain(2:end-1), newDomain);
end
numFuns = numel(funs);

% Initialise storage for new FUN objects and pointValues:
newFuns = cell(1, numel(newDomain)-1);

% Loop through each fun and restrict as required:
l = 0;
for k = 1:numFuns
    % Find the breaks which correspond to the kth fun:
    subsIdx = (newDomain >= oldDomain(k)) & (newDomain <= oldDomain(k+1));
    if ( sum(subsIdx) == 2 )
        % This interval is already a FUN: (i.e., no new breaks to introduce)
        newFuns{l+1} = restrict(funs{k}, newDomain(subsIdx));
        l = l + 1;
    else
        numSubs = sum(subsIdx)-1;
        if ( numSubs > 0 )
            % Restrict the FUN at the corresponding break points:
            newFuns(l+(1:numSubs)) = restrict(funs{k}, newDomain(subsIdx));
            l = l + numSubs;
        end
    end
end

% Update the pointValues:
newPointValues = chebfun.getValuesAtBreakpoints(newFuns);
% Restore existing pointValues:
[mask, locB] = ismember(oldDomain, newDomain);
locB(~logical(locB)) = [];
newPointValues(locB,:) = pointValues(mask,:);

% Attach data to CHEBFUN to return as output:
f.domain = newDomain;
f.funs = newFuns;
f.pointValues = newPointValues;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [tfr,t,f] = tfrmh(x,t,N,trace);
%TFRMH	Margenau-Hill time-frequency distribution.
%	[TFR,T,F]=TFRMH(X,T,N,TRACE) computes the Margenau-Hill
%	distribution of a discrete-time signal X,
%	or the cross Margenau-Hill representation between two signals.
%
%	X     : signal if auto-MH, or [X1,X2] if cross-MH.
%	T     : time instant(s)          (default : 1:length(X)).
%	N     : number of frequency bins (default : length(X)).
%	TRACE : if nonzero, the progression of the algorithm is shown
%                                        (default : 0).
%	TFR   : time-frequency representation. When called without
%               output arguments, TFRMH runs TFRQVIEW.
%	F     : vector of normalized frequencies.
%
%	Example :
%	 sig=fmlin(128,0.1,0.4); tfrmh(sig,1:128,128,1);
%
%	See also  all the time-frequency representations listed in
%	 the file CONTENTS (TFR*)

%	F. Auger, May-August 1994, July 1995.
%	Copyright (c) 1996 by CNRS (France).
%
%  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., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin == 0),
 error('At least 1 parameter required');
end;
[xrow,xcol] = size(x);
if (nargin == 1),
 t=1:xrow; N=xrow; trace=0;
elseif (nargin == 2),
 N=xrow; trace=0;
elseif (nargin == 3),
 trace = 0;
end;

if (N<0),
 error('N must be greater than zero');
end;
[trow,tcol] = size(t);
if (xcol==0) | (xcol>2),
 error('X must have one or two columns');
elseif (trow~=1),
 error('T must only have one row');
elseif (2^nextpow2(N)~=N),
 fprintf('For a faster computation, N should be a power of two\n');
end;

tfr= zeros (N,tcol) ;
if trace, disp('Margenau-Hill distribution'); end;
for icol=1:tcol,
 ti= t(icol); tau=-min([N-ti,xrow-ti]):(ti-1);
 indices= rem(N+tau,N)+1;
 if trace, disprog(icol,tcol,10); end;
 tfr(indices,icol)=x(ti,1)*conj(x(ti-tau,xcol));
end;

tfr= real(fft(tfr));

if (nargout==0),
 tfrqview(tfr,x,t,'tfrmh');
elseif (nargout==3),
 if rem(N,2)==0,
  f=[0:N/2-1 -N/2:-1]'/N;
 else
  f=[0:(N-1)/2 -(N-1)/2:-1]'/N;
 end;
end;

author: HeLiangHIT


text:
function SetJointVelocities(idx, dq)
global uLINK

for n=1:length(idx)
    j = idx(n);
    uLINK(j).dq = dq(n);
end

author: s-kajita

repo: IntroductionToHumanoidRobotics


text:
% function [face,skin_region]=face(I);
%
% skin_region=skin(I);
%
% se = strel('disk',3);
% dil = imdilate(skin_region,se);        % morphologic dilation
% d2 = imfill(dil, 'holes');   % morphologic fill
% face = bwdist(~d2);            % computing minimal euclidean distance to non-white pixel
% figure;imshow(face,[]);

function [face_a,skin_region]=face(I);

skin_region=skin(I);

se = strel('disk',5);
se2 = strel('disk',3);
er = imerode(skin_region,se2);
cl = imclose(er,se);
dil = imdilate(cl,se);        % morphologic dilation
dil = imdilate(dil,se);
cl2 = imclose(dil,se);
d2 = imfill(cl2, 'holes');   % morphologic fill
facearea = bwdist(~d2);            % computing minimal euclidean distance to non-white pixel
% figure;imshow(facearea,[]);

% imshow(d2);
face(:,:,1)=double(I(:,:,1)).*d2;
face(:,:,2)=double(I(:,:,2)).*d2;
face(:,:,3)=double(I(:,:,3)).*d2;
face_a=uint8(face);
% figure;imshow(face_a);

author: Sable

repo: mcbench-benchmarks


text:
function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
    = grw_start(fid, data, varargins, opts)
% This file is an entry for the EG strategy.
%
% function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
%            = grw_start(fid, data, varargins, opts)
% cum_ret: cumulative wealth achived at the end of a period.
% cumprod_ret: cumulative wealth achieved till the end each period.
% daily_ret: daily return achieved by a strategy.
% daily_portfolio: daily portfolios
%
% data: market sequence vectors
% fid: handle for write log file
% varargins: variable parameters
% opts: option parameter for behvaioral control
%
% Example: [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
%          = grw_start(fid, data, {0.00005, 0}, opts)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extract the parameters
sigma =varargins{1};    % Switching parameter
tc = varargins{2};      % transaction cost fee rate

% Run the grw algorithm multiple runs
N=[10 100 200];

fw = ones(3, 1);

for i = 1:3,
    [fw(i, 1)] = grw_run(fid, data, N(i), sigma, tc, opts);
end

A=[ones(1, 3); 1./N; 1./(N.^2)]';
beta = A\fw;

cum_ret = beta(1);

fprintf(1, 'cum_ret:%f\n', cum_ret);

% multiple times, ignore
cumprod_ret = 0;
daily_ret = 0;
daily_portfolio = 0;

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Eugenio Alcala Baselga
% Date: 02/06/2018
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ u_opt, u_predicted, x_predicted, x_err] = NLMPC_Kinematic_Computation...
    ( controller, Hp, KC, oldu, u_predicted, x_predicted, x_Real,...
    kin_states, kin_inputs, x_Ref, y_Ref, theta_Ref, vel_Ref, omega_Ref,...
    longitud_vect, car)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\\ Vehicle kinematic control:
%   controller: object provided by the MPC design
%   Hp: Prediction horizon
%   KC: discrete kinematic time
%   oldu: previous control action
%   u_predicted: vector of predicted control actions in the last iteration
%   x_predicted: vector of predicted states in the last iteration
%   x_Real: vector of measured vehicle states [x y theta v alpha w]
%   kin_states: number of kinematic states: 3
%   kin_inputs: number of kinematic control actions: 2
%   x_Ref: x reference (planner)
%   y_Ref: y reference (planner)
%   theta_Ref: theta reference (planner)
%   vel_Ref: velocity reference (planner)
%   omega_Ref: angular velocity reference (planner)
%   longitud_vect: length of vectors for simulation
%   car: vehicle object
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if KC==1
    % Initializing vectors:
    theta_cos           = ones(1,Hp+1);
    theta_sin           = zeros(1,Hp+1);
    theta_err_cos       = ones(1,Hp);
    U_OPT_VECTOR        = zeros(kin_inputs,Hp,longitud_vect);
    X_OPT_VECTOR        = zeros(kin_states,Hp,longitud_vect);
end

%%%% Constraints %%%%
% Control action constraints:
    % CA
        u_min_ini   = vehicle.u_min_ini;
        u_max_ini   = vehicle.u_max_ini;
        u_min       = vehicle.u_min;
        u_max       = vehicle.u_max;
    % derivative of CA
        du_min_ini  = vehicle.du_min_ini;
        du_max_ini  = vehicle.du_max_ini;
        du_min      = vehicle.du_min;
        du_max      = vehicle.du_max;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Terminal Set LMIs
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %[ PP, WW ] = TerminalSet_Continuous_Knemtic_Cmputtion();
    %[ PP, WW ] = TerminalSet_Discrete_Knemtic_Cmputtion(Ts_Kcontrol);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Errors de verdad (no forecasted from the LPV model):
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    x_err(1)    = (x_Ref(1)-x_Real(1))*cos(x_Real(3)) + (y_Ref(1)-x_Real(2))*sin(x_Real(3));
    x_err(2)    = -(x_Ref(1)-x_Real(1))*sin(x_Real(3)) + (y_Ref(1)-x_Real(2))*cos(x_Real(3));
    x_err(3)    = theta_Ref(1) - x_Real(3);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Moving Horizon Constraints:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if KC>1
        for j=1:Hp+1
            theta_cos(j) = cos(theta_Ref(1));
            theta_sin(j) = sin(theta_Ref(1));
            if(j>1)
                theta_err_cos(j-1) = cos(x_predicted(3,1)); % cos (theta_ERROR forecasted)
            end
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Computing PREDICTIONS of KIN-DYN LPV Vehicle Model. Note that
        % this A and B matrices are in continuous time due to the
        % discretization is carried out inside the controller object.
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        [Ac_pred, Bc_pred, Bref] = KinemError_LPV_Model_predicted(KC,Hp,x_predicted(3,:),...
            u_predicted(2,:), vel_Ref );
        for i=1:Hp
            Ac_pred(:,:,i) = Ac_pred(:,:,1);
            Bc_pred(:,:,i) = Bc_pred(:,:,1);
        end
    else
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Computing KIN-DYN LPV Vehicle Model for First iteration. Note that
        % this A and B matrices are in continuous time due to the
        % discretization is carried out inside the controller object.
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        [Ac_pred, Bc_pred, Bref] = KinemError_LPV_Model(KC,Hp,x_err(3),...
            x_Real(6), vel_Ref(1) );
    end

    %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Optimization stage:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%     if(KC > car.t_ini)
        inputs = {x_err', oldu, vel_Ref, omega_Ref,...
            theta_cos, theta_sin, u_min, u_max, du_min, du_max};
%     else
%         inputs = {x_err', oldu, vel_Ref, omega_Ref,...
%             theta_cos, theta_sin, u_min_ini, u_max_ini,...
%             du_min_ini, du_max_ini};
%     end

    [solutions,diagnostics] = controller{inputs};

    if diagnostics == 1
        u_opt = oldu;
        U_OPT_VECTOR(:,:,KC) = U_OPT_VECTOR(:,:,KC-1);
        X_OPT_VECTOR(:,:,KC) = X_OPT_VECTOR(:,:,KC-1);
    else
        U_OPT_VECTOR(:,:,KC) = solutions{1};
        u_opt = U_OPT_VECTOR(:,1,KC);
        X_OPT_VECTOR(:,:,KC) = solutions{2}; % Kinematic errors
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Predictions obtained from the solver:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    x_predicted = X_OPT_VECTOR(:,:,KC);
    u_predicted = U_OPT_VECTOR(:,:,KC);

end


text:
% Fig. 9.8  Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
 n=[1 1];
 d = [1 0 0];
 rlocus(n,d)
 axis([-6 2 -3 3])
 hold on
 r=roots([1 1 1]);
 plot(r,'*')
  z=0:.1:.9;
 wn= 1:6;
 sgrid(z, wn)
 hold off


text:
function [forecastmatrix]=forecastsim(data_endo_a,data_exo_p,beta,n,p,k,horizon)

% function [forecastmatrix]=forecastsim(data_endo_a,data_exo_p,beta,n,p,k,horizon)
% computes the matrix of unconditional forecasts, using the chain rule of forecasts (see technical guide p38)
% inputs:  - matrix 'data_endo_a': matrix of pre-forecast endogenous data
%          - matrix 'data_exo_p': predicted values for the exogenous variables over the forecast periods
%          - vector 'beta': vectorised form of VAR coefficients (definined in 1.1.12)
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - integer 'k': number of coefficients to estimate for each equation in the BVAR model (defined p 7 of technical guide)
%          - integer 'horizon': number of forecast periods
% outputs: - matrix 'forecast_matrix': matrix recording the forecast values

% uses the chain rule of forecasting, p38 of the technical guide

% compute the reduced matrix Y
Y=data_endo_a(end-p+1:end,:);

% reshape beta to obtain B
B=reshape(beta,k,n);

% repeat the process for periods T+1 to T+h
for jj=1:horizon

% step 1
% use the function lagx to obtain the matrix temp
temp=bear.lagx(Y,p-1);

% step 2
% define the reduced regressor matrix X
% if no exogenous variable is present at all in the model (neither constant nor other exogenous), define X only from the endogenous variables
if isempty(data_exo_p)==1
X=[temp(end,:)];
% if there are exogenous variables, concatenate them next to the endogenous
else
X=[temp(end,:) data_exo_p(jj,:)];
end

% step 3
% obtain the predicted value for T+jj
yp=X*B;

% step 4
% concatenate yp to the top of Y
Y=[Y;yp];

% repeat until values are obtained for T+h
end

% record the values in the matrix forecastmatrix
forecastmatrix=Y(p+1:end,:);

author: european-central-bank

repo: BEAR-toolbox


text:
% -*- INTERNAL UNDOCUMENTED FUNCTION -*-
function [geometry, msh, space_v, vel, space_p, press] = ...
                          solve_stokes_3d (problem_data, method_data)

warning ('geopdes:obsolete', 'Function SOLVE_STOKES_3D is obsolete. Using SOLVE_STOKES instead')

[geometry, msh, space_v, vel, space_p, press] = solve_stokes (problem_data, method_data);

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
DSP_MIN_SCALE = 1 / 6;
DSP_MAX_SCALE = 3;
DSP_NUM_SCALES = 10;

pool = gcp('nocreate');
if isempty(pool)
    pool = parpool(maxNumCompThreads());
end

parfor i = 1:num_images
    fprintf('Computing features for %s [%d/%d]', ...
            image_names{i}, i, num_images);

    if exist(keypoint_paths{i}, 'file') ...
            && exist(descriptor_paths{i}, 'file')
        fprintf(' -> skipping, already exist\n');
        continue;
    end

    tic;

    % Read the image for keypoint detection, patch extraction and
    % descriptor computation.
    image = imread(image_paths{i});
    if ismatrix(image)
        image = single(image);
    else
        image = single(rgb2gray(image));
    end

    % Read the pre-computed SIFT keypoints.
    keypoints = read_keypoints(keypoint_paths{i});

    % Compute the descriptors for the detected keypoints.
    if size(keypoints, 1) == 0
        descriptors = zeros(0, 128);
    else
        % Create DSP keypoints, sampled around the standard SIFT detections.
        num_keypoints = size(keypoints, 1);
        dsp_keypoints = zeros(4, num_keypoints, DSP_NUM_SCALES);
        dsp_scale_idx = 1;
        for dsp_scale = linspace(DSP_MIN_SCALE, DSP_MAX_SCALE, DSP_NUM_SCALES)
            dsp_keypoints([1 2 4],:,dsp_scale_idx) = keypoints(:,[1 2 4])';
            dsp_keypoints(3,:,dsp_scale_idx) = dsp_scale * keypoints(:,3)';
            dsp_scale_idx = dsp_scale_idx + 1;
        end

        % Transpose to VLFeat format.
        dsp_keypoints = reshape(dsp_keypoints, ...
                                [4, num_keypoints * DSP_NUM_SCALES]);

        % Extract the descriptors from the DSP keypoints.
        [~, descriptors] = vl_covdet(image, 'Frames', dsp_keypoints, ...
                                     'Descriptor', 'SIFT');

        % Aggregate the descriptors across all scales.
        descriptors = reshape(descriptors, ...
                              [128, num_keypoints, DSP_NUM_SCALES]);
        descriptors = mean(double(descriptors), 3)';
    end

    % Make sure that each keypoint has one descriptor.
    assert(size(keypoints, 1) == size(descriptors, 1));

    % Write the descriptors to disk for matching.
    write_descriptors(descriptor_paths{i}, descriptors);

    fprintf(' in %.3fs\n', toc);
end

author: ahojnnes

repo: local-feature-evaluation


text:
function mean = r8col_mean ( m, n, a )

%*****************************************************************************80
%
%% R8COL_MEAN returns the column means of an R8COL.
%
%  Example:
%
%    A =
%      1  2  3
%      2  6  7
%
%    MEAN =
%      1.5  4.0  5.0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), the array to be examined.
%
%    Output, real MEAN(N), the means, or averages, of the columns.
%
  for j = 1 : n
    mean(j) = sum ( a(1:m,j) );
  end

  mean(1:n) = mean(1:n) / m;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% L1 Filtering (Liu et al. 2011)
% process_video('RPCA', 'L1F', 'dataset/demo.avi', 'output/demo_L1F.avi');
[L,S] = rpca_l1f(M);

author: andrewssobral

repo: lrslibrary


text:
function [varargout]=pillowHex(E,V,C,shrinkFactor)

%Defive core element
[E_core,V_core]=scalePatch(E,V,shrinkFactor);

Vp=[V;V_core]; %Collect nodes
E_core=E_core+size(V,1); %Fix node indices in E_core

%% FORMAT FOR FACES
indTop=1:1:4;
indBottom=5:1:8;
indFront=[2 3 7 6];
indBack=[5 8 4 1];
indSide1=[1 2 6 5];
indSide2=[3 4 8 7];

%% THE TOP ELEMENT
%   top-of-core top-of-outer
E1=[E_core(:,indTop) E(:,indTop)];

%% THE BOTTOM ELEMENT
%   bottom-of-core bottom-of-outer
E2=[E_core(:,indBottom) E(:,indBottom)];

%% THE FRONT ELEMENT
%   bottom=front of E2 top=front of E1
E3=[E2(:,indFront) E1(:,indFront)];

%% THE BACK ELEMENT
%   bottom=back of E2 top=back of E1
E4=[E2(:,indBack) E1(:,indBack)];

%% THE SIDE1 ELEMENT
%   bottom=side1 of E2 top=side1 of E1
E5=[E2(:,indSide1) E1(:,indSide1)];

%% THE SIDE2 ELEMENT
%   bottom=side2 of E2 top=side2 of E1
E6=[E2(:,indSide2) E1(:,indSide2)];

%% Gather element sets
Ep=[E_core; E1; E2; E3; E4; E5; E6];
Cp=repmat(C,[7,1]);

%% Collect output
varargout{1}=Ep;
varargout{2}=Vp;
varargout{3}=Cp;

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


function  plot_shcks_dcmp_(pplotvar,ex_names_,leg,input,BVAR,options)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Filippo Ferroni, 6/1/2015
% Revised, 2/15/2017
% Revised, 3/21/2018

% pplotvar:     cell with the variable names for the decomp
% ex_names_:    shocks grouping
% leg:          shocks grouping names
% input:        outpu tof 'histdecomposition.m'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nshocks = BVAR.N;

for v = 1 : nshocks
    eval(['namesshock{' num2str(v) '} =  ''Shck ' num2str(v) ''';'])
end

TT   = 1:1:size(input,1);
Tlim   = [TT(1) TT(end)];

[~,positions] = ismember(pplotvar,BVAR.varnames);
% Setting the time
% i1     = length(TT);
% nplots = [1 1];
% T      = TT;

pplotvarname   = pplotvar;
dcmp_group_yes = 0;
dmcp_type      = 'stacked';
tags           = [];
colors_decomp_yes = 0;
savefig_yes = 0;
% initial_state_dcmp = 0;
% addplot_=0;
% addplot0_=0;

if nargin > 5
    if isfield(options,'time') ==1
        TT = options.time;
        if isfield(options,'Tlim') ==1;
            Tlim=options.Tlim;
            if Tlim(1) < TT(1)
                warning('You have set a intitial date that starts earlier than the first obs');
                warning('I change it with the frist obs');
                Tlim(1) = TT(1);
            end
            if Tlim(2) > TT(end)
                warning('You have set a final date that exceeds the forecast horizon');
                warning('I change it with the endo of forecast');
                Tlim(2) = TT(end);
            end
        end

    end

    if isfield(options,'tags') ==1
        tags = options.tags;
    end
    if isfield(options,'dcmp_grouped') ==1 && options.dcmp_grouped ==1
        dcmp_group_yes = 1;
        dmcp_type      = 'grouped';
    end
    if isfield(options,'plotvarnames') ==1
        pplotvarname = options.plotvarnames;
        if length(pplotvarname)~=length(pplotvar)
            error('The number of plot titles (pplotvarname) and of plot variables needs to coincide')
        end
    end
    if isfield(options,'colors_decomp') ==1 && options.colors_decomp==1
        colors_decomp_yes =1;
    end
    if isfield(options,'colors_decomp') ==1 && options.colors_decomp==2
        colors_decomp_yes =2;
    end
    if isfield(options,'colors_decomp') ==1 && options.colors_decomp==3
        colors_decomp_yes =3;
    end
    if isfield(options,'saveas_dir') ==1
        savefig_yes = 1;
        %   setting the folder where to save the figure
        fnam_dir = options.saveas_dir;
        if exist(fnam_dir,'dir') == 0
            mkdir(fnam_dir)
        end
    end

    %     if isfield(options,'addplots_yes') ==1  && options.addplots_yes==1,
    %         tags = [ tags '_memo'];
    %         if isfield(input,'frcsts') ==1
    %             addplot_=1;
    %             frcsts = input.frcsts.states.Mean(:,positions);
    %             s_s = input.frcsts.states.steady(positions);
    %         end
    %         if isfield(input,'frcsts0') ==1
    %             addplot0_=1;
    %             frcsts0 = input.frcsts0.states.Mean(:,positions);
    %
    %         end
    %     end
end

deco  = input;

% setting the names of the figure to save
fnam_suffix = [tags '_shcks_dcmp'];

ngroups0 = size(ex_names_,1);
% ngroups  = ngroups0+1+no_initial_effect;
if colors_decomp_yes == 0
    func = @(x) colorspace('RGB->Lab',x);
    MAP = distinguishable_colors(ngroups0+1,'w',func);
    % MAP = CreateColorMap(ngroups0+1);
    MAP(end,:) = [0.7 0.7 0.7];
elseif colors_decomp_yes == 1
    MAP = zeros(4,3);
    MAP(end,:) = [0.7 0.7 0.7]; % gray
    MAP(end-1,:) = [0.2 0.5 0.99];      % blue
    MAP(end-2,:) =  [1 1 0];        % yellow
    MAP(end-3,:) = [0.8 0 0.8] ;            % purple
    %     MAP(end-4,:) = [0 0.7 0];          % green
    %     MAP(end-5,:) = [1 0 0];             % red
    %     MAP(end-6,:) = [0 0 1];             % blue
    %     MAP(end-7,:) = [.5 .5 0];           % light green
elseif colors_decomp_yes == 2
    MAP = zeros(10,3);
    MAP(end,:) = [0.7 0.7 0.7]; % gray
    MAP(end-1,:) = [0.2 0.5 0.99];       % blue
    MAP(end-2,:) = [0.8 0 0.8] ;        % purple
    MAP(end-3,:) = [1 1 0] ;            % yellow
    MAP(end-4,:) = [0 0.7 0];          % green
    MAP(end-5,:) = [1 0 0];             % red
    MAP(end-6,:) = [0 0 1];             % blue
    MAP(end-7,:) = [.5 .5 0];          % light green
    MAP(end-8,:) = [0 0.25 0]; %brown
    MAP(end-9,:) = [0 0.9 0.9]; % violet
elseif colors_decomp_yes == 3
    MAP = zeros(11,3);
    MAP(end,:) = [0.7 0.7 0.7]; % gray
    MAP(end-1,:) = [0.2 0.5 0.99];       % blue
    MAP(end-2,:) = [0.8 0 0.8] ;        % purple
    MAP(end-3,:) = [1 1 0] ;            % yellow
    MAP(end-4,:) = [0 0.7 0];          % green
    MAP(end-5,:) = [1 0 0];             % red
    MAP(end-6,:) = [0 0 1];             % blue
    MAP(end-7,:) = [.5 .5 0];          % light green
    MAP(end-8,:) = [0 0.25 0]; %brown
    MAP(end-9,:) = [0 0.9 0.9]; % violet
    MAP(end-10,:) = [0.5 0.1 0.1]; % violet
end

st = find(Tlim(1)==TT);
en = find(Tlim(2)==TT);
if savefig_yes == 1
    fidTxt = fopen([fnam_dir '\legenda_' fnam_suffix '_plots.txt'],'w');
    fprintf(fidTxt,['LEGENDA ' tags ' SHOCK DECOMPOSITION PLOTS\n']);
    fprintf(fidTxt,['\n']);
end

for j = 1 : size(pplotvar,2)
    clear sdec sdec_tot,
    indx = positions(j);
    sdec0 = squeeze(deco(:,indx,:));
    for i=1:ngroups0
        clear index,
        for ii=1:size(ex_names_{i},2)
            indbuf = strmatch(ex_names_{i}{ii},namesshock,'exact');
            if ~isempty(indbuf)
                index(ii) = indbuf;
            elseif ~isempty(ex_names_{i}{ii})
                error(['Shock name ',ex_names_{i}{ii}, ' not found.' ]);
            end
        end
        sdec(:,i)=sum(sdec0(:,index),2);
        sdec0(:,index)=0;
    end

    h= figure('Name',['Shocks Decomposition for ' pplotvarname{j}]);
    sdec_tot=[sdec,  sum(sdec0,2)];
    if dcmp_group_yes == 0
        ind_pos = (sdec_tot>0);
        ind_neg = (sdec_tot<0);
        temp_neg  = cumsum(sdec_tot.*ind_neg ,2).*ind_neg;
        temp_pos  = cumsum(sdec_tot.*ind_pos ,2).*ind_pos;
        temp      = temp_neg + temp_pos;
        for kk = size(temp,2) : -1 : 1
            hold on
            bbar = bar(TT(st:en),temp(st:en,kk),dmcp_type,'EdgeColor',[0 0 0]);
            set(bbar, 'FaceColor', MAP(kk,:))
            shading faceted; hold on;
        end
        leg0 = leg(end:-1:1);

    else
        temp      = sdec_tot;
        bbar = bar(TT(st:en),temp(st:en,:),dmcp_type,'EdgeColor',[0 0 0]); colormap(MAP);
        %         hleg=legend(leg(1:1:end),'interpreter','none','location','Best');
        %         shading faceted;
        %
    end
    %     set(hleg,'position',[0.5 0.15 0.4 0.2],'units','normalized')
    hold on
    axis tight
    fillips = sum(sdec_tot,2);
    hold on, h1=plot(TT(st:en),fillips(st:en),'k-d');
    set(h1,'MarkerFaceColor', 'k')

    %     if addplot_ ==1
    %         hold on, h1=plot(TT(st:en),frcsts(st:en,j),'k-*','LineWidth',2);
    %         hold on, h1=plot(TT(st:en),s_s(j)*ones(length(TT(st:en)),1),'k-.','LineWidth',2);
    %         leg0{end+1} = 'Current Forecast';
    %         leg0{end+1} = 'steady state';
    %     end
    %
    %     if addplot0_ ==1
    %         hold on, h1=plot(TT(st:en),frcsts0(st:en,j),'r-o','LineWidth',2);
    %         leg0{end+1} = 'Previous Forecast';
    %     end

    set(gcf,'position' ,[50 50 800 650])
    hleg=legend(leg0,'interpreter','none','location','Best');
    shading faceted;

    %     define_.timestart + (define_.nobs-1)/4
    %plot the last in-+sample obs
    %     hold on; plot([TT(define_.nobs) TT(define_.nobs)],[low up],'color',[1 0 0],'LineWidth',2)
    %     h=vline(define_.timestart + (define_.nobs-1)/4, 'k', 'Last Obs');
    %     set(h,'Linewidth',2)

    title(pplotvarname{j})
    %     set(gca,'Xtick',TT(st:6:en))
    %     tmp_str= sample2date(TT(st:6:en));
    set(gca,'Xtick',TT(st:6:en))
    %     tmp_str= sample2date(TT(st:2:en));
    %     set(gca,'Xticklabel',tmp_str)
    %     tmp_str= sample2date(TT);
    %     STR_RECAP = [ 'model_' fnam_suffix '_' tmp_str{st} '_'  tmp_str{end} '_' int2str(j)];

    if savefig_yes == 1
        STR_RECAP = [ fnam_dir '\svar_' fnam_suffix '_' int2str(j)];

        saveas(gcf,STR_RECAP,'fig');
        if strcmp(version('-release'),'2022b') == 0
            saveas(gcf,STR_RECAP,'pdf');
        end
        fprintf(fidTxt,['The figure sdcmp_' fnam_suffix '_' int2str(j) ' contains the following variable:\n' ]);
        tmp = strrep(char(pplotvarname{j} ),'\',' ');
        fprintf(fidTxt,[tmp '; ']);

        fprintf(fidTxt,['\n' ]);
        fprintf(fidTxt,['\n']);
        %close all
    end

end

if savefig_yes == 1
    fclose(fidTxt);
end

end


text:
function [xi,w]=fifthOrderTriangleCubPoints()
%%FIFTHORDERTRIANGLECUBPOINTS Obtain fifth-order cubature points for
%   integration over a triangle in 2D. The points and weights are for the
%   triangle with vertices (1,0), (0,1), (0,0), but can be transformed to
%   any triangle using transformSimplexTriPoints.
%
%INPUTS: None
%
%OUTPUTS: xi A 2XnumCubPoints set of points for the standard triangle.
%          w A 1XnumCubPoints set of cubature weights. This sums to the
%            volume of the triangle (1/2).
%
%This function implements the points given in [1] (7 points).
%
%EXAMPLE:
%Given the vertices of the simplex, we compare a fifth-order moment
%computed using these cubature points to one computed using
%monomialIntSimplex. The results are the same within typical finite
%precision limits.
% [xi,w]=fifthOrderTriangleCubPoints();
% alpha=[3;2];
% theMoment=findMomentFromSamp(alpha,xi,w)
% intVal=monomialIntSimplex(alpha)
%
%REFERENCES:
%[1] F. D. Witherden and P. E. Vincent, "On the identification of symmetric
%    quadrature rules for finite element methods," Computer and Mathematics
%    with Applications, vol. 69, no. 10, pp. 1232-1241, May 2015.
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

M=[-0.33333333333333333333333333333333333333,   -0.33333333333333333333333333333333333333,                                        0.45;
   -0.79742698535308732239802527616975234389,    0.59485397070617464479605055233950468778,    0.25187836108965430519136789100036266732;
    0.59485397070617464479605055233950468778,   -0.79742698535308732239802527616975234389,    0.25187836108965430519136789100036266732;
   -0.79742698535308732239802527616975234389,   -0.79742698535308732239802527616975234389,    0.25187836108965430519136789100036266732;
  -0.059715871789769820459117580973104798968,   -0.88056825642046035908176483805379040206,    0.26478830557701236147529877566630399935;
   -0.88056825642046035908176483805379040206,  -0.059715871789769820459117580973104798968,    0.26478830557701236147529877566630399935;
  -0.059715871789769820459117580973104798968,  -0.059715871789769820459117580973104798968,    0.26478830557701236147529877566630399935];
w=M(:,3);
xi=M(:,1:2)';
%Transform the points to the standard triangle.
v1=[-1,-1, 1;
    -1, 1,-1];
v2=[1,0,0;
    0,1,0];
[A,d]=affineTransBetweenTriangles(v1,v2);
xi=bsxfun(@plus,A*xi,d);
w=w/4;
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% vgg_wedge  Wedge product of N-1 N-vectors (generalization of cross product).
%
% Y = vgg_wedge(X)  Wedge product of columns/rows of X.
%   Y ... double (1,N).
%   X ... double (N,N-1).
% It is Y = X(:,1) \wedge X(:,2) \wedge ... \wedge X(:,N-1). For N=3,
% wedge product is the same as cross (vector) product. E.g., for N=4,
% wedge product of three 4-vectors is in fact computation of a plane in
% 3-D projective space from 3 points in the plane.
% The sign is chosen so that for any square matrix X it is: det(X) == wedge(X(:,1:end-1))*X(:,end)
% Works also dually for
%   Y ... double (N,1).
%   X ... double (N-1,N).
%
% Y = vgg_wedge(X_1,X_2,...,X_{N-1})  Wedge product for each (N-1)-tuple of corresponding columns of X_n.
%   X_n ... double (N,K)
%   Y ... double (K,N)
% Equivalent to
%   for k = 1:K, Y(k,:) = wedge([X_1(:,k) ... X_{N-1}(:,k)]); end
% E.g.: wedge(X1,X2) is the same as cross(X1,X2)' but faster.
% Dual form is not available.

function Y = vgg_wedge(varargin)

if nargin == 1

  X = varargin{1};

  [N,Nm1] = size(X);
  if Nm1>N
    Y = vgg_wedge(X')';
    return
  end

  switch N
   case 3 % make it faster for special case N==3
    Y = [X(2,1).*X(3,2)-X(3,1).*X(2,2),...
         X(3,1).*X(1,2)-X(1,1).*X(3,2),...
         X(1,1).*X(2,2)-X(2,1).*X(1,2)];
   otherwise
    for n = 1:N
      Y(n) = (-1)^(n+N)*det(X([1:n-1 n+1:N],:));
    end
  end

else

  N = nargin + 1;
  switch N
   case 3 % make it faster for special case N==3
    X1 = varargin{1};
    X2 = varargin{2};
    Y = [(X1(2,:).*X2(3,:)-X1(3,:).*X2(2,:))',...
         (X1(3,:).*X2(1,:)-X1(1,:).*X2(3,:))',...
         (X1(1,:).*X2(2,:)-X1(2,:).*X2(1,:))'];
   otherwise
    for k = 1:size(varargin{1},2)
      for n = 1:N-1
	X(:,n) = varargin{n}(:,k);
      end
      Y(k,:) = vgg_wedge(X);
    end
  end

end

return

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
function [fDeltaProbability, fProbabilityN, fProbabilityH, mPredictionFMD, vObservedFMD, vMagnitudeBins, vProbH, vProbN] ...
        = pt_poissonian(mLearningCatalog, fLearningTime, mObservedCatalog, fObservedTime, nMinimumNumber, fBValueH, fBValueN, ...
        fMcH, fMcN, fMinMag, fMaxMag)
    % Calculation of likelihoods for two models based on a- and b-values
    %
    % [fDeltaProbability, fProbabilityN, fProbabilityH, mPredictionFMD, vObservedFMD, vMagnitudeBins, vProbH, vProbN]
    %   = pt_poissonian(mLearningCatalog, fLearningTime, mObservedCatalog, fObservedTime, nMinimumNumber, fBValueH, fBValueN,
    %                   fMcH, fMcN, fMinMag, fMaxMag);
    %
    %
    % Input parameters:
    %   mLearningCatalog    Earthquake catalog of the learning period
    %   fLearningTime       Length of learning period (can be different than the exact length of mLearningCatalog)
    %   mObservedCatalog    Earthquake catalog of the observed period
    %   fObservedTime       Length of observed period (can be different than the exact length of mObservedCatalog)
    %   nMinimumNumber      Minimum number of earthquakes in the catalog for calculating the output values
    %   fBValueH            b-value for test hypothesis
    %   fBValueN            b-value for null hypothesis
    %   fMcH                Magnitude of completeness for the test hypothesis
    %   fMcN                Magnitude of completeness for the null hypothesis
    %   fMinMag             Minimum magnitude for testing
    %   fMaxMag             Maximum magnitude for testing
    %
    % Output parameters:
    %   fDeltaProbability   Difference of the log-likelihood between both of the models
    %   fProbabilityN       Log-likelihood of the null hypothesis
    %   fProbabilityH       Log-likelihood of the test hypothesis
    %   mPredictionFMD      Forecasted FMD of both models (mPredictionFMD(:,1) test hypothesis, (:,2) null hypothesis)
    %   vObservedFMD        Vector with observations corresponding to mPredictionFMD
    %   vMAgnitudeBins      Vector containing the tested magnitude bins
    %   vProbH              log-likelihoods per magnitude bin of the test hypothesis
    %   vProbN              log-likelihoods per magnitude bin of the null hypothesis
    %
    % Danijel Schorlemmer
    % April 30, 2003

    report_this_filefun();

    % Index description
    % xxxxN : Variable with a value for the null hypothesis model
    % xxxxH : Variable with a value for the hypothesis model

    if ((length(mLearningCatalog(:,1)) < nMinimumNumber) | (isnan(fBValueH)) | (isnan(fBValueN)))
        fDeltaProbability = nan;
        fProbabilityH = nan;
        fProbabilityN = nan;
        mPredictionFMD = [];
        vObservedFMD = [];
        vMagnitudeBins = [];
        vProbH = [];
        vProbN =[];
    else
        % Select the higher magnitude of completeness
        fMc = max([fMcH fMcN]);
        % Cut the first catalog at magnitude of completeness
        vSel = mLearningCatalog(:,6) >= fMc;
        mCalcACatalog = mLearningCatalog(vSel,:);
        if length(mCalcACatalog(:,1)) >= nMinimumNumber
            % Compute the maximum likelihood a-values for both models
            fAValueH = calc_MaxLikelihoodA(mCalcACatalog, fBValueH);
            fAValueN = calc_MaxLikelihoodA(mCalcACatalog, fBValueN);
            % Create the predicted FMD
            vCnt = (fMinMag:0.1:fMaxMag+0.1)'; % Add one more magnitude bin for later use of diff()
            vNumberH = 10.^(fAValueH - (fBValueH * vCnt))/fLearningTime * fObservedTime;
            vNumberN = 10.^(fAValueN - (fBValueN * vCnt))/fLearningTime * fObservedTime;
            mNumbers = [vNumberH vNumberN];
            % Determine the number of events in each magnitude bin
            mPredictionFMD = -diff(mNumbers);
            % Cut the second catalog at the minimum magnitude for testing
            vSel = mObservedCatalog(:,6) >= fMinMag;
            mCalcCatalog = mObservedCatalog(vSel,:);
            % Create the FMD for the period of observation
            vObservedFMD = histogram(mCalcCatalog(:,6), fMinMag:0.1:fMaxMag);
            % Calculate the likelihoods for both of the models
            vProbH = calc_logpoisspdf(vObservedFMD', mPredictionFMD(:,1));
            vProbN = calc_logpoisspdf(vObservedFMD', mPredictionFMD(:,2));
            % Return the values
            fProbabilityH = sum(vProbH);
            fProbabilityN = sum(vProbN);
            fDeltaProbability = fProbabilityN - fProbabilityH;
            vObservedFMD = vObservedFMD';   % Simply for returning the vector in the same way as vPredictionFMD
            vMagnitudeBins = (fMinMag:0.1:fMaxMag)';
        else
            fDeltaProbability = nan;
            fProbabilityH = nan;
            fProbabilityN = nan;
            mPredictionFMD = [];
            vObservedFMD = [];
            vMagnitudeBins = [];
            vProbH = [];
            vProbN =[];
        end
    end
end

author: CelsoReyes


text:
function pass = test_roots08( pref )
% Check that the marching squares and Bezoutian agree with each other.

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

%%
f = chebfun2(@(x,y)sin(10*x-y/10));
g = chebfun2(@(x,y)cos(3*x.*y));
r1 = roots([f;g],'ms');
r2 = roots([f;g],'resultant');
pass(j) = ( norm(sort(r1(:,1))-sort(r2(:,1))) < tol ); j = j + 1;
pass(j) = ( norm(sort(r1(:,2))-sort(r2(:,2))) < tol ); j = j + 1;

%%
f = chebfun2(@(x,y)sin(10*x-y/10) + y);
g = chebfun2(@(x,y)cos(10*y-x/10) - x);
r1 = roots([f;g],'ms');
r2 = roots([f;g],'resultant');
pass(j) = ( norm(sort(r1(:,1))-sort(r2(:,1))) < tol ); j = j + 1;
pass(j) = ( norm(sort(r1(:,2))-sort(r2(:,2))) < tol ); j = j + 1;

end


text:
function y = r8vec_sqctf ( n, x )

%*****************************************************************************80
%
%% R8VEC_SQCTF computes a "slow" quarter cosine transform forward of an R8VEC.
%
%  Discussion:
%
%    This routine is provided for illustration and testing.  It is inefficient
%    relative to optimized routines that use fast Fourier techniques.
%
%    For 0 <= I <= N-1,
%
%      Y(I) = (1/N) Sum ( 0 <= J <= N-1 ) X(J) * cos ( PI * I * (J+1/2) / N )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    William Briggs, Van Emden Henson,
%    The Discrete Fourier Transform,
%    SIAM, 1995,
%    QA403.5 B75
%
%  Parameters:
%
%    Input, integer N, the number of data values.
%
%    Input, real X(N), the data sequence.
%
%    Output, real Y(N), the transformed data.
%
  y(1:n) = 0.0;

  for i = 1 : n
    for j = 1 : n
      theta = 0.5 * pi * ( i - 1 ) * ( 2 * j - 1 ) / n;
      y(i) = y(i) + x(j) * cos ( theta  );
    end
  end

  y(1:n) = y(1:n) / n;

  return
end


text:
%typ==0 -> nav=geodetic_frame, typ==1 -> nav=wander_frame
%although this implementation can be used as wander frame mechanization,
%this is still singular as it explicitly computes the wander angle for
%transport rate computation.
%For non-singular implementation see strapdown_wander_quat (which uses
%geowander to compute earth curvature)

function [Cbn_new, Vn_new, Cen_new, h_new]=strapdown_Cen_dcm(Cbn, Vn, Cen, h, a, w, dt, typ)

%%Geo Params
[ll, wander]=dcm2llh_v000(Cen);
Cgn=[wander(2) wander(1) 0;-wander(1) wander(2) 0;0 0 1];
llh=[ll;h];
[Rn, Re, g, sL, cL, WIE_E]=geoparam_v000(llh);

wie_n=Cen*[0; 0; WIE_E];
wen_n=Cgn*[0 1/(Re+llh(3)) 0; -1/(Rn+llh(3)) 0 0;0 -sL/cL/(Re+llh(3)) 0]*Cgn'*Vn;

if (typ==1) %wander transport rate
    wen_n(3)=0;
end

%%%Update attitude
%Part I:Body frame update
rot=w*dt;
rot_norm=norm(rot);
sr_a=1-(rot_norm^2/6)+(rot_norm^4/120);
sr_b=(1/2)-(rot_norm^2/24)+(rot_norm^4/720);
mx_a=eye(3)+sr_a*skew(rot)+sr_b*skew(rot)*skew(rot);

%Part II:Nav Frame update
rot=-(wen_n+wie_n)*dt;
rot_norm=norm(rot);
sr_a=1-(rot_norm^2/6)+(rot_norm^4/120);
sr_b=(1/2)-(rot_norm^2/24)+(rot_norm^4/720);
mx_b=eye(3)+sr_a*skew(rot)+sr_b*skew(rot)*skew(rot);

Cbn_new=mx_b*Cbn*mx_a;

%%%Update Velocity
vel_inc1=(Cbn*(a*dt))+[0;0;g]*dt;
vel_inc2=(cross(Vn,2*wie_n+wen_n))*dt;
Vn_new=Vn+vel_inc1+vel_inc2;

%%%Update Cen (position + wander_angle)
rot=-wen_n*dt;
rot_norm=norm(rot);
sr_a=1-(rot_norm^2/6)+(rot_norm^4/120);
sr_b=(1/2)-(rot_norm^2/24)+(rot_norm^4/720);
mx_b=eye(3)+sr_a*skew(rot)+sr_b*skew(rot)*skew(rot);
Cen_new=mx_b*Cen;

%update height
h_new=h-Vn(3)*dt;

repo: ekfmonoslam


text:
function Path = TestPathInGraph(DAG,Start,End)
% Input:DAG is the graph, p the start point, q the end point.
% Output:if there is a path from p to q, return Path = 1,otherwise Path = 0;

% test the cycle by DFS algorithm
Path = 0;
Visited = zeros( 1,size( DAG,1 ) );
[ List,Visited ] = SearchPath( DAG,Start,[],Visited );

    while isempty( List ) == 0
         %Node = List(1)
         Descendant = List(1);
         if Descendant == End
             Path = 1;
             break;
         end
         List( 1 )=[];
         if Visited( Descendant ) ==0
            [List,Visited] = SearchPath(DAG,Descendant,List,Visited);
         end
    end
end

function [ List,Visited ] = SearchPath( DAG,Start,L,V )
 V( Start ) = 1;
 Visited=V;
 Descendant = find( DAG(Start,:) ~= 0 );
 List = [L Descendant];
end


text:
function b_n = beamWeightsHypercardioid2Spherical(N)
%BEAMWEIGHTSHYPERCARDIOID2SPHERICAL Generate spherical beamweights for hypercardioids
%
%   The hypercardioid is the pattern that maximizes the directivity-factor
%   for a certain SH order N. The hypercardioid is also the plane-wave
%   decomposition beamformer in the SHD, also called 'regular' because the
%   beamweights are just the SH values on the beam-direction.
%   Since the pattern is axisymmetric only the N+1 coefficients of m=0 are
%   returned.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% BEAMWEIGHTSHYPERCARDIOID2SPHERICAL.M - 11/7/2013
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c_n = 4*pi/(N+1)^2 * getSH(N,[0 0],'real');
b_n = zeros(N+1,1);
for n=0:N
    b_n(n+1) = c_n((n+1)^2 -n);
end

author: polarch

repo: Spherical-Array-Processing


text:
function [A] = spm_matrix(P, order)
% returns an affine transformation matrix
% FORMAT [A] = spm_matrix(P, order)
% P(1)  - x translation
% P(2)  - y translation
% P(3)  - z translation
% P(4)  - x rotation about - {pitch} (radians)
% P(5)  - y rotation about - {roll}  (radians)
% P(6)  - z rotation about - {yaw}   (radians)
% P(7)  - x scaling
% P(8)  - y scaling
% P(9)  - z scaling
% P(10) - x affine
% P(11) - y affine
% P(12) - z affine
%
% order (optional) application order of transformations.
%
% A     - affine transformation matrix
%___________________________________________________________________________
%
% spm_matrix returns a matrix defining an orthogonal linear (translation,
% rotation, scaling or affine) transformation given a vector of
% parameters (P).  By default, the transformations are applied in the
% following order (i.e., the opposite to which they are specified):
%
% 1) shear
% 2) scale (zoom)
% 3) rotation - yaw, roll & pitch
% 4) translation
%
% This order can be changed by calling spm_matrix with a string as a
% second argument. This string may contain any valid MATLAB expression
% that returns a 4x4 matrix after evaluation. The special characters 'S',
% 'Z', 'R', 'T' can be used to reference the transformations 1)-4)
% above. The default order is 'T*R*Z*S', as described above.
%
% SPM uses a PRE-multiplication format i.e. Y = A*X where X and Y are 4 x n
% matrices of n coordinates.
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id$

% pad P with 'null' parameters
%---------------------------------------------------------------------------
q  = [0 0 0 0 0 0 1 1 1 0 0 0];
P  = [P q((length(P) + 1):12)];

% default multiplication order if not specified
%---------------------------------------------------------------------------
if nargin < 2
    order = 'T*R*Z*S';
end;

T  =   [1   0   0   P(1);
        0   1   0   P(2);
        0   0   1   P(3);
        0   0   0   1];

R1  =  [1    0      0          0;
        0    cos(P(4))  sin(P(4))  0;
        0   -sin(P(4))  cos(P(4))  0;
        0    0      0          1];

R2  =  [cos(P(5))  0    sin(P(5))  0;
        0          1    0      0;
       -sin(P(5))  0    cos(P(5))  0;
        0          0    0          1];

R3  =  [cos(P(6))   sin(P(6))   0  0;
       -sin(P(6))   cos(P(6))   0  0;
        0           0           1  0;
        0           0       0  1];

R   = R1*R2*R3;

Z   =  [P(7)    0       0       0;
        0       P(8)    0       0;
        0       0       P(9)    0;
        0       0       0       1];

S   =  [1       P(10)   P(11)   0;
        0       1   P(12)   0;
        0       0       1   0;
        0       0       0       1];

A = eval(sprintf('%s;', order));
if ~isnumeric(A) || ndims(A) ~= 2 || any(size(A) ~= 4)
    error('Order expression ''%s'' did not return a valid 4x4 matrix.', ...
          order);
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function out=errorplot(x,y,dy,plottypes,dx);
%
%  2 dimensional plot with error bars
%
%	By: P.S.Basran		6/17/03
%
%       Toronto-Sunnybrook Cancer Centre
%       Dept. Medical Physics
%       Parminder.Basran@sw.ca
%
%
%
%	Function:
%		This function produces a two dimensional plot with error bars.
%
%	Syntax:
% 		out=errorplot(x,y,dy,plottypes,dx)
%
%	Inputs:
%		    x - is a [1 x n] dimensional vector
%           y - is a [m x n] dimensional matrix for multidimensional plots
%          dy - is a [m x n] dimensional matrix whose dimensions must match those of y
%   plottypes - is a [1 x 4] dimensional vector whose entries specify the plottype in the
%               same fashion as the variable 's' in the function plot(x,y,s).
%               NOTE: current version of this requires the plottype to be specified in
%               all instances, and to be 4 characters in length.
%          dx - is a real number that specifices the dimensions of the horizontal tics
%               for each error bar. The default value is 0.25 the dimension of x.
%
%	Outputs:
%		   - plot with legend, labeled '1', '2', etc.
%
%   Example:
%
%           %For a single plot:
%
%           x=[-10:1:10];
%           y1=x.^2;
%           dy1=0.3*y1.*rand(size(x)); %
%           plottype1=[':   '];
%           errorplot(x,y1,dy1,plottype1);
%
%           % For multiple plots:
%
%           y2=0.5*x.^2;
%           dy2=0.2*y2.*rand(size(x)); %
%           plottype2=['o-- '];
%           y=[y1; y2];
%           dy=[dy1; dy2];
%           errorplot(x,y,dy,[plottype1 plottype2],1);
%

%   Modification Log:
%       June 17 - 2003: alpha version
%       March 11 - 2004: fix for the plottypes ....thanks to Sean Verret
%

error(nargchk(2,5,nargin));

if nargin <5, dx=(x(2)-x(1))/4;
    if nargin < 4, plottypes='';
        if nargin < 3, dy='';
            [a,b]=size(y);
            for i=1:a, plot(x,y(i,:)); end % default to a regular plot for 2 inputs
        end, end, end

%Provide plot and legend
[a,b]=size(y);
if a>b,
    y=y';
    [a,b]=size(y);
end

tx=[];
for i=1:a,
    plot(x,y(i,:),plottypes(4*(i-1)+1:4*i));;
    tx=[tx; 'plot ' num2str(i)];
    hold on;
end
legend(tx);

%Now the error bars
for i=1:a,
    for j=1:b,
        %horizontal lines
        line([x(j) x(j)], [y(i,j)-dy(i,j) y(i,j)+dy(i,j)]);
        %top vertical line
        line([x(j)-dx x(j)+dx], [y(i,j)+dy(i,j) y(i,j)+dy(i,j)]);
        %bottom vertical line
        line([x(j)-dx x(j)+dx], [y(i,j)-dy(i,j) y(i,j)-dy(i,j)]);
     end
end

hold off


text:
function [timeZoneInfo,displayName]=getTimeZoneInfo(timeZoneName,Jul1,Jul2)
%GETTIMEZONEINFO Get information about a time zone. If a date is provided,
%                one can obtain an offset from UTC taking daylight savings
%                time into account, if relevant. Note, however, that many
%                countries only add a leap second at their local midnight
%                and not when UTC does, so the offsets might be off by a
%                second around when a leap second is added.
%
%INPUTS: timeZoneName The name of the time zone. This is one of the short
%                   names returned by getTimeZoneList(). A meaningless
%                   input will cause UTC to be returned.
%           Jul1,Jul2 Optionally, the time as a pseudo-Julian date in UTC.
%                   The units of the date are days. The full date is the
%                   sum of both terms. The date is broken into two parts to
%                   provide more bits of precision. It does not matter how
%                   the date is split. If this input is provided, then
%                   daylight savings time information can be provided on
%                   the output.
%
%OUTPUTS: timeZoneInfo A 2X1 or 4X1 (if a UTC date is provided) providing
%                   information on the time zone. timeZoneInfo(1) is the
%                   offset from UTC in seconds for time time zone, without
%                   daylight savings time. timeZoneInfo(2) is true if the
%                   time zone is listed as having ever supported daylight
%                   savings time or if it is listed as ever supporting it
%                   in the future. timeZoneInfo(3) is the offset in seconds
%                   from UTC at the given date and timeZoneInfo(4) is true
%                   if the time zone is in daylight savings time at the
%                   given date.
%       displayName A more detailed version of timeZoneName.
%
%This function just calls the appropriate Java commands to obtain time zone
%information.
%
%EXAMPLE:
%For the Easten Time Zone of the United States on 11 January 2016, one can
%use
% [UTC1,UTC2]=Cal2UTC(2016,1,11,0,0,0);
% [timeZoneInfo,displayName]=getTimeZoneInfo('EST5EDT',UTC1,UTC2)
%One will see that daylight savings time is honored, but is not active at
%the given time, so both offsets are -18000 second (=-5 hours).
%
%May 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

timeZone=java.util.TimeZone.getTimeZone(timeZoneName);

displayName=char(timeZone.getDisplayName());

if(nargin>1&&~isempty(Jul1))
    timeZoneInfo=zeros(4,1);
    timeZoneInfo(1)=timeZone.getRawOffset()/1e3;
    timeZoneInfo(2)=timeZone.observesDaylightTime() || timeZone.useDaylightTime();

    [TT1,TT2]=UTC2TT(Jul1,Jul2);
    [UTCRef1,UTCRef2]=Cal2UTC(1970,1,1,0,0,0);
    [TTRef1,TTRef2]=UTC2TT(UTCRef1,UTCRef2);

    %It should be that TT1>TT2 and TTRef1>TTRef2 due to how the conversion
    %functions work.
    %There are 86400 seconds in all Julian days in terrestrial time. The
    %multiplication by 1000 turns it into milliseconds.
    timeDiffMilliSeconds=int64((TT1-TTRef1)*86400*1000)+int64((TT2-TTRef2)*86400*1000);

    timeZoneInfo(3)=timeZone.getOffset(timeDiffMilliSeconds)/1e3;
    timeZoneInfo(4)=timeZoneInfo(1)~=timeZoneInfo(3);%Must currently be in DST.
else
    timeZoneInfo=zeros(2,1);
    timeZoneInfo(1)=timeZone.getRawOffset()/1e3;
    timeZoneInfo(2)=timeZone.observesDaylightTime() || timeZone.useDaylightTime();
end
end

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


text:
function [bbci, data]= bbci_calibrate_csp_tiny(bbci, data)
%BBCI_CALIBRATE_CSP_TINY - Calibrate for SMR Modulations with CSP (Tiny ver.)
%
%This function is called by bbci_calibrate
%(if BBCI.calibate.fcn is set to @bbci_calibrate_csp_tiny).
%Via BBCI.calibrate.settings, the details can be specified, see below.
%
%Synopsis:
% [BBCI, DATA]= bbci_calibrate_csp_tiny(BBCI, DATA)
%
%Arguments:
%  BBCI -  the field 'calibrate.settings' holds parameters specific to
%          calibrate CSP-based BCI processing.
%  DATA -  holds the calibration data
%
%Output:
%  BBCI - Updated BBCI structure in which all necessary fields for
%     online operation are set, see bbci_apply_structures.
%  DATA - As input.
%
%BBCI.calibrate.settings may include the following parameters:
%  ival: [1x2 DOUBLE] interval on which CSP is performed.
%  band: [1x2 DOUBLE] frequency band on which CSP is performed.
%  clab: [CELL] Labels of the channels that are used for classification,
%     default {'not','E*','Fp*','AF*','OI*','I*','*9','*10'}.
%  nPatters: [INT>0] number of CSP patterns which are considered from each
%     side of the eigenvalue spectrum. Note, that not neccessarily all of
%     these are used for classification, see settings.pattern.
%     Default is 3.
%  model: [CHAR or CELL] Classification model.
%     Default {'RLDAshrink', 'gamma',0, store_means',1, 'scaling',1}.

% 11-2011 Benjamin Blankertz

default_clab=  {'not','E*','Fp*','AF*','OI*','I*','*9','*10'};
default_model= {@train_RLDAshrink, 'Gamma',0, 'StoreMeans',1, 'Scaling',1};

props= {'clab'         default_clab   'CELL{CHAR}'
        'ival'         [750 3750]     '!DOUBLE[1 2]'
        'band'         [8 33]         '!DOUBLE[1 2]'
        'nPatterns'    3              '!INT'
        'model'        default_model  'FUNC|CELL'
        'filtOrder'    5              '!INT'
       };
opt= opt_setDefaults('bbci.calibrate.settings', props);

[filt_b,filt_a]= butter(opt.filtOrder, opt.band/data.cnt.fs*2);
cnt_flt= proc_filt(data.cnt, filt_b, filt_a);

bbci.signal.clab= data.cnt.clab(util_chanind(data.cnt, opt.clab));

fv= proc_segmentation(cnt_flt, data.mrk, opt.ival, 'clab',bbci.signal.clab);
[fv_csp, csp_w, A, la]= proc_csp(fv, 'SelectFcn',...
                                 {@cspselect_equalPerClass, opt.nPatterns});

bbci.signal.proc= {{@online_linearDerivation, csp_w}, ...
                   {@online_filt, filt_b, filt_a}};

bbci.feature.ival= [-750 0];
bbci.feature.proc= {@proc_variance, @proc_logarithm};

fv_csp= bbci_calibrate_evalFeature(fv_csp, bbci.feature);
bbci.classifier.C= trainClassifier(fv_csp, opt.model);

bbci.quit_condition.marker= 255;

author: bbci


text:
% Initialization code for running Garmin serial
% input model for xPC

Ts=1/75; % sample time in seconds


text:
function [xh,x,t]=simobsv(G,H)
[y,t,x]=step(G);
G=ss(G); A=G.a; B=G.b; C=G.c; D=G.d;
[y1,xh1]=step((A-H*C),(B-H*D),C,D,1,t);
[y2,xh2]=lsim((A-H*C),H,C,D,y,t);
xh=xh1+xh2;


text:
function denoised = KSVD_WRAP(image_path, sigma)

    Noisy_Image = im2double(imread(image_path));

    redChannel = Noisy_Image(:,:,1); % Red channel
    greenChannel = Noisy_Image(:,:,2); % Green channel
    blueChannel = Noisy_Image(:,:,3); % Blue channel

    redChannel=im2double(redChannel);
    if (length(size(redChannel))>2)
        redChannel = rgb2gray(redChannel);
    end
    if (max(redChannel(:))<2)
        redChannel = redChannel*255;
    end

    greenChannel=im2double(greenChannel);
    if (length(size(greenChannel))>2)
        greenChannel = rgb2gray(greenChannel);
    end
    if (max(greenChannel(:))<2)
        greenChannel = greenChannel*255;
    end

    blueChannel=im2double(blueChannel);
    if (length(size(blueChannel))>2)
        blueChannel = rgb2gray(blueChannel);
    end
    if (max(blueChannel(:))<2)
        blueChannel = blueChannel*255;
    end

    bb=8; % block size
    RR=2; % redundancy factor
    K=RR*bb^2; % number of atoms in the dictionary

    [Denoised_Image_Red,~] = denoiseImageKSVD(redChannel, sigma,K);
    [Denoised_Image_Green,~] = denoiseImageKSVD(greenChannel, sigma,K);
    [Denoised_Image_Blue,~] = denoiseImageKSVD(blueChannel, sigma,K);

    recombinedRGBImage = cat(3, Denoised_Image_Red, Denoised_Image_Green, Denoised_Image_Blue);

    denoised = uint8(255 * mat2gray(recombinedRGBImage));

end

author: lbasek

repo: image-denoising-benchmark


text:
function [cl] = l2cl(l)
% Convert volume from liters to centiliters.
% Chad Greene 2012
cl = l*100;


text:
% Script demo_load_fileformats - to be excluded
% Shows versatile file format loading capabilities of MrImage.load
%
% demo_load_fileformats
%
%
% See also

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2017-02-14
% Copyright (C) 2017 Institute for Biomedical Engineering
% University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
% <http://www.gnu.org/licenses/>.

pathExamples    = tapas_uniqc_get_path('data');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1. Load different types of nifti
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% a) load data from fileName and updates both name and parameters.save of
% nifti files, header is read to update MrImage.parameters
fileName = fullfile(pathExamples, 'nifti', 'rest', 'fmri_short.nii');

Y1 = MrImage(fileName);
disp(Y1);
disp(Y1.dimInfo);

%% b) load data from fileName with different load options

% selectedVolumes: Note that selectedVolumes is specific to loading nifti,
% cpx and par/rec files and refers to the 4rd dimension of the loaded data;
% the advantage is it avoids loading the full data and then selecting a subset.
Y2 = MrImage(fileName, 'selectedVolumes', 3:2:9);
disp(Y2.dimInfo);

% updateProperties: Note that per default the name of the MrImage object is
% set to the file. For more information see MrDataNd.read_single_file. In
% the example here, the save parameters are set to the load path. Per
% default, a new path depending on the pwd is created to prevent accidental
% overwrite. Be careful with this option!
Y3 = MrImage(fileName, 'updateProperties', 'save');
disp(Y3.parameters.save);

%% c) load data from fileName with additinal dimInfo information

% default: The dimInfo information is gathered from the header (same as Y1,
% nothing to do here).
Y4 = MrImage(fileName);

% fileName_dimInfo.mat: An additional dimInfo.mat file exisits.
Y4.dimInfo.units = {'this', 'and', 'that', 'too'};
Y4.save;

% load file, dimInfo is automatically added
Y5 = MrImage(Y4.get_filename);
disp(Y5.dimInfo);

% compare to only loading the nifti file
[fp, fn] = fileparts(Y4.get_filename);
delete(fullfile(fp, [fn, '_dimInfo.mat']));
Y6 = MrImage(Y4.get_filename);
disp(Y6.dimInfo);

% restore correct dimInfo via dimInfo argument
Y7 = MrImage(Y4.get_filename, 'dimInfo', Y5.dimInfo);

disp(Y7.dimInfo);

% or, alternatively, via prop/val pair
Y8 = MrImage(Y4.get_filename, 'units', Y4.dimInfo.units);
disp(Y8.dimInfo)

%% d) load filename from data with additional affine geometry
Y9 = MrImage(fileName, 'shear', [0.2 -0.1 0]);
disp(Y9.geometry);

%% 2. Load multiple files in folder
% a) load multiple .nii files in folder with filenames containing additional
% dimension information but no additinal _dimInfo.mat files
fileNameSplitSubset = fullfile(pathExamples, 'nifti', 'split', 'subset');
YSplitSubset = MrImage(fileNameSplitSubset);
disp(YSplitSubset.dimInfo);

% b) load multiple nifti files in folder without the filenames containing
% dimension information
fileNameSplitRes = fullfile(pathExamples, 'nifti', 'split', 'split_residual_images');
YSplitRes = MrImage(fileNameSplitRes);
disp(YSplitRes.dimInfo);

% c) load multiple nifti files in folder with filenames containing
% additional dimension information and select
fileNameSplitFull = fullfile(pathExamples, 'nifti', 'split', 'subset');
select.TE = 2;
select.t = 3;
YSplitSelect = MrImage(fileNameSplitFull, 'select', select);

% d) load multiple nifti files in folder with filenames containing
% additional dimension information and _dimInfo file
fileNameSplitDimInfo = fullfile(pathExamples, 'nifti', 'split', 'with_dim_info');
YSplitDimInfo = MrImage(fileNameSplitDimInfo);

% e) load multiple nifti files in folder with filenames containing
% additional dimension information and select which is not a dimension of
% the files
selectError.TE = 2;
selectError.t = 7;
selectError.doesNotExist = 3;
YSplitError = MrImage(fileNameSplitFull, 'select', selectError);

% f) select for MrSeries
YMrSeries = MrSeries(fileName, 'select', {'t', 4:2:8});
disp(YMrSeries.data.dimInfo);

author: translationalneuromodeling

repo: tapas


text:
function V=mrLoadRet3StatsToAnalyze(outFileName,scanNum,mapName)
%
% V=mrLoadRet3StatsToAnalyze(outFileName,scanNum,mapName)
%
% AUTHOR:  Wade
% DATE: 06.10.03
% PURPOSE:
% Converts a statistical map in a mrLoadRet 3 Gray view into a 8 bit
% analyze file. Data are scaled from 0 to 2^8-1
% Uses the current scan if no scanNum is passed
% mapName can be 'co', 'amp', 'ph' or any other similar statistical map in
% the VOLUME structure.
% If no mapName is passed, it defaults to writing out 'co'.
% If embedDimensions are passed, the data are embedded into an even larger
% volume. This is useful when
% dealing with BV volume anatomies that have been embedded into 256x256x256
%
% RETURNS:
% Number of bytes written in main data block.
%
% EXAMPLE:
% 1) V=mrLoadRet3StatsToAnalyze('temp');
% 2) V=mrLoadRet3StatsToAnalyze('temp_ph','ph');
%
% NOTES
% Originall coded to transfer mrLoadRet gray stats into SSI's EMSE / mrViewer
% $Author: wade $
% $Date: 2003/09/09 21:18:59 $
mrGlobals;
if (isempty(selectedVOLUME))
    error('You must select (click within) a volume window before proceeding');
end

if ((~exist('scanNum','var')) | (isempty(scanNum)))
    scanNum=getCurScan(VOLUME{selectedVOLUME});
end
% Get a coherence threshold
coThresh=getCothresh(VOLUME{selectedVOLUME});

if ((~exist('outFileName','var')) | (isempty(outFileName)))
    a=dataTYPES(VOLUME{selectedVOLUME}.curDataType).name;
    outFileName=[a,'_',int2str(scanNum),'_co_',int2str(coThresh*100)];
end

if ((~exist('mapName','var')) | (isempty(mapName)))
    mapName='co';
end
volAnatSize=size(VOLUME{selectedVOLUME}.anat);
% We need to permute this into the correct format.
volAnatSize=volAnatSize([2 1 3]);
if ((~exist('embedDimensions','var')) | (isempty(embedDimensions)))
    embedDimensions=volAnatSize;
end
if ((~exist('VOLUME','var')) | (isempty(VOLUME{selectedVOLUME}.co)))
    error('mrLoadRet must be running and the GRAY corAnal must be loaded');
end
if (isempty(VOLUME{selectedVOLUME}.anat))
    error('The anatomy must be loaded');
end

% % Find the size of the full anatomy
coords=VOLUME{selectedVOLUME}.coords;
coords=coords';
%
% minVals=min(coords);
% maxVals=max(coords);
% bbSize=maxVals-minVals+1
%
% if(find(bbSize<=0))
%     error('The GRAY data are in an invalid VOI');
% end

dataVolume=uint8(ones(size(VOLUME{selectedVOLUME}.anat)));
[ySiz xSiz zSiz]=size(VOLUME{selectedVOLUME}.anat);

% Do the mother of all sub2inds to get a linear index into dataVolume
coords=sub2ind([ySiz xSiz zSiz],coords(:,1),coords(:,2),coords(:,3));
if (~strcmp(mapName,'co'))
    error('Only co maps allowed for now');
end

co=(VOLUME{selectedVOLUME}.co(scanNum));
co=co{1};

% Now threshold according to the current cothresh
co(co<coThresh)=0;

blurKern=ones(3,3,3);
blurKern=blurKern./sum(blurKern(:));
blurKern=blurKern*20;

co=convn(co,blurKern,'same');
co(co>=1)=0.9999; % Just to make sure
co(co<=0)=0;

% Scale co from 0 to 2^8-1
max(co);
min(co);
% co=log10(co*10);
co(co<=0)=0;
co=co*((2^8)-1);

max(co)
min(co)

% Send it into the large data volume
dataVolume(coords)=(co);

dataVolume=permute(dataVolume,[3 2 1]);
dataVolume=flipdim(dataVolume,3);
%dataVolume=flipdim(dataVolume,2);
dataVolume=flipdim(dataVolume,1);
% Call SPM writevol routine to write out the data.

    s=spm_hwrite(outFileName,[ySiz xSiz zSiz],[1 1 1],1,spm_type('uint8'),0);
    V=spm_vol(outFileName);

    V.descrip=['Converted from tSeries file in session',mrSESSION.sessionCode,' : ',mrSESSION.subject,':  on ',datestr(now)];

    s=spm_write_vol(V,double(dataVolume));

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% step1:
% use
%     pipeline = pipeline_init('upper');
% for initiating a fashion landmark detection pipeline
%
% step2:
% use
%     prediction = pipeline_forword(img,pipeline);
% for getting landmarks of images
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initiate pipeline
addpath(genpath('../caffe/matlab'));
pipeline = pipeline_init('upper');

image_path = '../data/FLD_upper/';

name_list = dir(image_path);

for i = 1:length(name_list)

    if name_list(i).isdir
        continue;
    end
    img_name = [image_path, name_list(i).name];
    img = imread(img_name);

    % forward
    prediction = pipeline_forword(img, pipeline);

    % show result
    pipeline_show_results(img, prediction);

end

% release pipeline
pipeline_release;

repo: fashion-landmarks


text:
function p = emptyNumericalModel
p.F_struc = [];
p.K.f = 0;
p.K.l = 0;
p.K.q = [];
p.K.p = [];
p.K.e = 0;
p.K.s = 0;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function determ = daub10_determinant ( n )

%*****************************************************************************80
%
%% DAUB10_DETERMINANT returns the determinant of the DAUB10 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real DETERM, the determinant.
%
  determ = + 1.0;

  return
end


text:
function [p,r] = cs_scc (A)                                                 %#ok
%CS_SCC strongly-connected components of a square sparse matrix.
%   [p,r] = cs_scc(A) finds a permutation p so that A(p,p) is permuted into
%   block upper triangular form.  The diagonal of A is ignored.  The kth block
%   is given by A (s,s) where s = r(k):r(k+1)-1.  A must be square.
%   For bipartite or rectangular graphs, use cs_scc2.
%
%   Example:
%       Prob = UFget ('HB/arc130') ; A = Prob.A ; [p r] = cs_scc (A) ;
%       cspy (A (p,p)) ;
%
%   See also CS_DMPERM, DMPERM, CS_SCC2.

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

error ('cs_scc mexFunction not found') ;


text:
function Qdec = num2Qdec(M,Vars)
% Qdec = num2Qdec(M,Vars)
%
% Converts numeric matrix M to two dimensional cell string array Qdec.
%
% Input
% M: Numeric matrix.
% Vars: One dimensional cell string array with the name of the Qdec
% variables corresponding to the columns of M.
%
% Output
% Qdec: Two dimensional cell string array of Qdec data .
%
% $Revision: 1.2 $  $Date: 2015/01/06 17:14:53 $
% Original Author: Jorge Luis Bernal Rusiel
% CVS Revision Info:
%    $Author: mreuter $
%    $Date: 2015/01/06 17:14:53 $
%    $Revision: 1.2 $
%
if nargin < 2
    error('Too few inputs');
end;
szM = size(M);
Qdec = cell(szM(1)+1,szM(2));
Qdec(1,1:end) = Vars;
for i=2:szM(1)+1
    for j=1:szM(2)
        Qdec{i,j} = num2str(M(i-1,j));
    end;
end;


text:
% M = MEAN2(MTX)
%
% Sample mean of a matrix.

function res = mean2(mtx)

res = mean(mean(mtx));

repo: SelfExSR


text:
function y = Evaluate_FT(x,shift,boxlen,center,w)
% Evaluate_FT -- 1d unequispaced Fourier transform followed by windowing
% Usage:
%   y = Evaluate_FT(x,shift,boxlen,center,w)
% Inputs:
%   x	      vector of length n
%   shift     vector of shifts
%   boxlen    half length of interval around each shift
%   center    boolean variable: 0/1 = unbiased/biased FT
%   w         tapering window of size 2*boxlen
% Outputs:
%   y         vector of length: (# shifts) * (2*boxlen)
% Description:
%  Evaluates the FT at the frequency points
%                   omega(j,k) = shift(j) + k, -boxlen <= k <boxlen
%  If center = 0,
%
%  y = sum_{-n/2 <= t < n/2} exp(-i 2pi omega t/n) x(t), -n/2 <= t < n/2
%
%  If center = 1,
%
%  y = sum_{0 <= t < n} exp(-i 2pi omega t/n) x(t), 0 <= t < n
%
%  For each shift parameter (shift(j)), the outcome y_j is then
%  tapered with the window w
%
%  Evaluate_FT uses two different methods for computing the
%  unesquispaced Fourier transform. If the number of shifts is less
%  than 8, it uses an exact method, otherwise it uses the
%  accelerated USFFT.
%  See Also
%    Adj_Evaluate_FT, USFFT, USFT_simple
%
% By Emmanuel candes, 2003-2004

  if nargin < 5,
    w = ones(1,2*boxlen);
  end

  if nargin < 4,
    center = 0;
  end

  n = length(x);
  boxcnt = length(shift);

  if (boxcnt <= 16)
	  y = USFT_simple(x,shift,boxlen,center,w); % nonequispaced FT
  else
    % Make frequency grid
    k = (-boxlen):(boxlen -1);
    omega = meshgrid(shift,k) + meshgrid(k,shift).';
    omega = 2*pi*(omega(:))./n;

    y = USFFT(x,omega,16,6,center)./sqrt(n); % Take USFFT
    y = repmat(w(:),boxcnt,1) .* y; % Tapering
  end

author: thomaskuestner


text:
function labeled_output = PTKWatershedFromStartingPoints( ~, ~ )
    % PTKWatershedFromStartingPoints Implementation of watershed flooding algorithm
    %
    %     This is a Matlab mex file and must be compiled before use.
    %
    %     To compile, type
    %
    %         mex PTKWatershedFromStartingPoints
    %
    %     in the Matlab command window.

    error('PTKWatershedMeyerFromStartingPoints has not been compiled. You must compile using mex PTKWatershedFromStartingPoints. Alternatively, use the slower Matlab equivalent function mex PTKWatershedFromStartingPointsMatlab().');
end

author: tomdoel

repo: pulmonarytoolkit


text:
%SFLAG4 Shape function driver (fourth order P4/Q4 Lagrange polynomials).
%
%   [ VBASE, NLDOF, XLDOF, SFUN ] = SFLAG4( I_EVAL, N_SDIM, N_VERT, I_DOF, XI, AINVJAC, VBASE )
%   Evaluates conforming quartic Lagrange shape functions with values
%   defined in the nodes, and edges, (and also faces, and cell centers for
%   quadrilaterals, tetrahedra, and hexahedrals).
%
%       Input       Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       i_eval      scalar:  1             Evaluate function values
%                           >1             Evaluate values of derivatives
%       n_sdim      scalar: 1-3            Number of space dimensions
%       n_vert      scalar: 2-8            Number of vertices per cell
%       i_dof       scalar: 1-n_ldof       Local basis function to evaluate
%       xi          [n_sdim(+1)]           Local coordinates of evaluation point
%       aInvJac     [n,n_sdim(+1)*n_sdim]  Inverse of transformation Jacobian
%       vBase       [n]                    Preallocated output vector
%                                                                                .
%       Output      Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       vBase       [n]                    Evaluated function values
%       nLDof       [4]                    Number of local degrees of freedom on
%                                          vertices, edges, faces, and cell interiors
%       xLDof       [n_sdim,n_ldof]        Local coordinates of local dofs
%       sfun        string                 Function name of called shape function
%
%   See also SF_LINE_P4, SF_TRI_P4, SF_TET_P4, SF_QUAD_Q4, SF_HEX_Q4

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function test3
%TEST3 test cs_lsolve, cs_ltsolve, cs_usolve, cs_chol
%
% Example:
%   test3
% See also: testall

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

index = UFget ;
[ignore f] = sort (max (index.nrows, index.ncols)) ;
f = f (1:100) ;

clf
% f = f(1)

for i = f
    Prob = UFget (i) ;
    disp (Prob) ;
    A = Prob.A ;
    [m n] = size (A) ;
    if (~isreal (A) | m ~= n)                                               %#ok
        continue
    end

    A = A*A' + 2*n*speye (n) ;
    try
        p = amd (A) ;
    catch
        p = symamd (A) ;
    end
    try
        L0 = chol (A)' ;
    catch
        continue
    end
    b = rand (n,1) ;

    C = A(p,p) ;
    c = condest (C) ;
    fprintf ('condest: %g\n', c) ;

    x1 = L0\b ;
    x2 = cs_lsolve (L0,b) ;
    err = norm (x1-x2,1) ;
    if (err > 1e-12 * c)
        error ('!') ;
    end

    x1 = L0'\b ;
    x2 = cs_ltsolve (L0,b) ;
    err = norm (x1-x2,1) ;
    if (err > 1e-10 * c)
        error ('!') ;
    end

    U = L0' ;

    x1 = U\b ;
    x2 = cs_usolve (U,b) ;
    err = norm (x1-x2,1) ;
    if (err > 1e-10 * c)
        error ('!') ;
    end

    L2 = cs_chol (A) ;
    subplot (2,3,1) ; spy (L0) ;
    subplot (2,3,4) ; spy (L2) ;
    err = norm (L0-L2,1) ;
    if (err > 1e-8 * c)
        error ('!') ;
    end

    L1 = chol (C)' ;
    L2 = cs_chol (C) ;
    subplot (2,3,2) ; spy (L1) ;
    subplot (2,3,5) ; spy (L2) ;
    err = norm (L1-L2,1) ;
    if (err > 1e-8 * c)
        error ('!') ;
    end

    [L3,p] = cs_chol (A) ;
    C = A(p,p) ;
    L4 = chol (C)' ;
    subplot (2,3,3) ; spy (L4) ;
    subplot (2,3,6) ; spy (L3) ;
    err = norm (L4-L3,1) ;
    if (err > 1e-8 * c)
        error ('!') ;
    end

    drawnow

end


text:
function x=anasing(N,t0,h) ;
%ANASING Lipschitz singularity.
%	X=ANASING(N,T0,H) generates the N-points Lipschitz singularity
%	centered around T=T0 : X(T) = |T-T0|^H.
%
%	N  : number of points in time
%	T0 : time localization of the singularity  (default : N/2)
%	H  : strenght of the Lipschitz singularity (positive or
%	     negative)                             (default : 0.0)
%	X  : the time row vector containing the signal samples
%
%	Example :
%	 x=anasing(128); plot(real(x));
%
%	See also anastep, anapulse, anabpsk, doppler.

%	P. Goncalves - September 1995
%	Copyright (c) 1995 Rice University
%
%  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., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin == 0),
 error ( 'The number of parameters must be at least 1.' );
elseif (nargin == 1),
 t0=N/2; h=0.0;
elseif (nargin ==2),
 h=0.0;
end;

if (N<=0),
 error('N must be greater or equal to 1.');
elseif (t0<=0) | (t0>N),
 error('t0 must be between 1 and N.');
else
 if h <= 0
   f = (1/N:1/N:0.5-1/N) ;
   y = zeros(1,N/2);
   y(2:N/2) = (f.^(-1-h)).*exp(-i*2*pi*f.*(t0-1));
   x = real(ifft(y,N)) ;
   x = x./max(x);
   x = x.' - sign(min(x))*abs(min(x)) ;
 else
   t = 1:N;
   x = abs(t-t0).^h;
   x = max(x)-x.' ;
 end
end

x=hilbert(x);

author: HeLiangHIT