text:
function CM = confMatrix( IDXtrue, IDXpred, ntypes )
% Generates a confusion matrix according to true and predicted data labels.
%
% CM(i,j) denotes the number of elements of class i that were given label
% j.  In other words, each row i contains the predictions for elements whos
% actual class was i.  If IDXpred is perfect, then CM is a diagonal matrix
% with CM(i,i) equal to the number of instances of class i.
%
% To normalize CM to [0,1], divide each row by sum of that row:
%  CMnorm = CM ./ repmat( sum(CM,2), [1 size(CM,2)] );
%
% USAGE
%  CM = confMatrix( IDXtrue, IDXpred, ntypes )
%
% INPUTS
%  IDXtrue     - [nx1] array of true labels [int values in 1-ntypes]
%  IDXpred     - [nx1] array of predicted labels [int values in 1-ntypes]
%  ntypes      - maximum number of types (should be > max(IDX))
%
% OUTPUTS
%  CM          - ntypes x ntypes confusion array with integer values
%
% EXAMPLE
%  IDXtrue = [ones(1,25) ones(1,25)*2];
%  IDXpred = [ones(1,10) randint2(1,30,[1 2]) ones(1,10)*2];
%  CM = confMatrix( IDXtrue, IDXpred, 2 )
%  confMatrixShow( CM, {'class-A','class-B'}, {'FontSize',20} )
%
% See also CONFMATRIXSHOW
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.12
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

IDXtrue=IDXtrue(:); IDXpred=IDXpred(:);

%%% convert common binary labels [-1/+1] or [0/1] to [1/2]
if( ntypes==2 )
  IDX = [IDXtrue;IDXpred];
  if( min(IDX)>=-1 && max(IDX)<=1 && all(IDX~=0))
    IDXtrue=IDXtrue+2;  IDXpred=IDXpred+2;
    IDXtrue(IDXtrue==3) = 2;  IDXpred(IDXpred==3) = 2;
  elseif( min(IDX)>=0 && max(IDX)<=1 )
    IDXtrue=IDXtrue+1;  IDXpred=IDXpred+1;
  end
end

%%% error check
[IDXtrue,er] = checkNumArgs( IDXtrue, [], 0, 2 ); error(er);
[IDXpred,er] = checkNumArgs( IDXpred, [], 0, 2 ); error(er);
if( length(IDXtrue)~=length(IDXpred) )
  error('Lengths of IDXs must match up.'); end
if( max([IDXtrue;IDXpred])>ntypes )
  error(['ntypes = ' int2str(ntypes) ' not large enough']); end

%%% generate CM
CM = zeros(ntypes);
for i=1:ntypes
  vals = IDXpred( IDXtrue==i );
  for j=1:ntypes; CM(i,j) = sum(vals==j); end
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function paths = getPaths(doMkdir)
	%% Data dir
	paths.dataDir = fullfile('/work3/', 'sgupta', 'cvpr13Release', 'data');
		paths.pcDir = fullfile(paths.dataDir, 'pointCloud');
		paths.colorImageDir = fullfile(paths.dataDir, 'colorImage');

	%% Cache dir
	paths.cacheDir = fullfile('/work3', 'sgupta', 'cvpr13Release', 'cachedir');

	RUNNAME = 'release';
	paths.runDir = fullfile(paths.cacheDir, RUNNAME);

	paths.cacheDir = fullfile(paths.runDir, 'cache');
		paths.gradientsForRecognition = fullfile(paths.cacheDir, 'gradientsRecog');
		paths.colorCues = fullfile(paths.cacheDir, 'colorCues');
		paths.depthCues = fullfile(paths.cacheDir, 'depthCues');
		paths.featuresDir = fullfile(paths.cacheDir, 'features');
		paths.mapDir = fullfile(paths.cacheDir, 'map');
		paths.featureCache = fullfile(paths.cacheDir, 'featureCache');
		paths.outsDir = fullfile(paths.cacheDir, 'outs');
		paths.sceneOutsDir = fullfile(paths.cacheDir, 'sceneOuts');
			paths.ucmFDir = fullfile(paths.cacheDir, 'ucmFeatures', 'features');
			paths.ucmFCacheDir = fullfile(paths.cacheDir, 'ucmFeatures', 'cache');
			paths.ucmGTDir = fullfile(paths.cacheDir, 'ucmFeatures', 'ucm_gto');

	paths.modelDir = fullfile(paths.runDir, 'model');
		paths.categroySpecificModels = fullfile(paths.modelDir, 'categorySpecific');
		paths.ucmModels = fullfile(paths.modelDir, 'ucm');

	paths.outDir = fullfile(paths.runDir, 'output');
		paths.ssOutDir = fullfile(paths.outDir, 'semanticSegmentation');
		paths.amodalDir = fullfile(paths.outDir, 'amodal');
		paths.ucmDir = fullfile(paths.outDir, 'ucm');

	%% Path for the image stack library ..
	pathstr = fileparts(mfilename('fullpath'));
		paths.siftLib = fullfile(pathstr, '..', 'external', 'colorDescriptor');
		paths.imageStackLib = fullfile(pathstr, '..', 'external', 'ImageStack');

	if(exist('doMkdir', 'var') && doMkdir)
		f = fieldnames(paths);
		for i = 1:length(f)
			d = paths.(f{i});
			if(~exist(d, 'dir'))
				mkdir(d);
			end
		end
	end
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
%GETSIZE Dataset size and number of classes
%
%  [M,K,C] = GETSIZE(A,DIM)
%
% INPUT
%   A    Dataset
%   DIM  1,2 or 3 : the number of the output argument to be returned
%
% OUTPUT
%   M    Number of objects
%   K    Number of features
%   C    Number of classes
%
% DESCRIPTION
% Returns size of the dataset A and the number of classes. C is determined
% from the number of labels stored in A.LABLIST. If DIM = 1,2 or 3, just
% one of these numbers is returned, e.g. C = GETSIZE(A,3).

author: marianux

repo: ecg-kit

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


text:
% Chapter 8 - Planar Systems.
% Program_8a - Program to Plot a Vector Field.
% Save M-file as vectorfield.m. Do NOT run this function file.
% Copyright Birkhauser 2013. Stephen Lynch.

% See the phase protraits in the book.
function vectorfield(deqns,xval,yval,t)
if nargin==3;
    t=0;
end
m=length(xval);
n=length(yval);
x1=zeros(n,m);
y1=zeros(n,m);
for a=1:m
  for b=1:n
    pts = feval(deqns,t,[xval(a);yval(b)]);
    x1(b,a) = pts(1);
    y1(b,a) = pts(2);
  end
end
arrow=sqrt(x1.^2+y1.^2);
quiver(xval,yval,x1./arrow,y1./arrow,.5,'r');
axis tight;

% End of Program_8a.

author: Sable

repo: mcbench-benchmarks


text:
%% Get data from an HTML table straight into MATLAB
% Copyright 2008 - 2010 The MathWorks, Inc.

%% Navigate to a webpage with a table
% You have to use the MATLAB browser because this takes advantage of the
% matlabcolon (matlab:) protocol.
web('http://finance.yahoo.com/q/ks?s=GOOG')
% give it time to load the page
pause(3)
%% Call the function to grab the data
getTableFromWeb
%% Updated Browser
% This shows what the browser looks like.  This is just a screenshot but it
% would be whatever page you are on with MATLAB icons next to tables.
im = imread('mlPage.jpg');
imshow(im);
%% Select the table
% You can select the table interactively by clicking on the MATLAB icon or
% passing a number to the getTableFromWeb function.  I have selected to get
% the Valuation Measures
myTableData = getTableFromWeb(8)


text:
function [RIR_cell] = ISM_RIR_bank_GPU(setupstruc,RIRFileName,varargin)
%ISM_RIR_bank  Bank of RIRs using Lehmann & Johansson's image-source method
%
% [RIR_CELL] = ISM_RIR_bank(SETUP_STRUC,RIR_FILE_NAME)
% [RIR_CELL] = ISM_RIR_bank( ... ,'arg1',val1,'arg2',val2,...)
%
% This function generates a bank of room impulse responses (RIRs) for a
% particular user-defined room setup, using Lehmann and Johansson's
% implementation of the image-source method (see: "Prediction of energy
% decay in room impulse responses simulated with an image-source model", J.
% Acoust. Soc. Am., vol. 124(1), pp. 269-277, July 2008). The input
% SETUP_STRUC is a structure of enviromental parameters containing the
% following fields:
%
%          Fs: sampling frequency (in Hz).
%        room: 1-by-3 vector of enclosure dimensions (in m),
%              [x_length y_length z_length].
%     mic_pos: N-by-3 matrix, [x1 y1 z1; x2 y2 z2; ...] positions of N
%              microphones (in m).
%    src_traj: M-by-3 matrix, [x1 y1 z1; x2 y2 z2; ...] positions of M
%              source trajectory points (in m).
%  T20 or T60: scalar value (in s), desired reverberation time.
%           c: (optional) sound velocity (in m/s).
% abs_weights: (optional) 1-by-6 vector of absorption coefficients weights,
%              [w_x1 w_x2 w_y1 w_y2 w_z1 w_z2].
%
% If the field SETUP_STRUC.c is undefined, the function assumes a default
% value of sound velocity of 343 m/s.
%
% The field 'abs_weight' corresponds to the relative weights of each of the
% six absorption coefficients resulting from the desired reverberation time.
% For instance, defining 'abs_weights' as [1 1 0.8 0.8 0.6 0.6] will result
% in the absorption coefficients (alpha) for the walls in the y-dimension
% being 20% smaller compared to the x-dimension walls, whereas the floor
% and ceiling will end up with absorption coefficients 40% smaller (e.g.,
% to simulate the effects of a concrete floor and ceiling). If this field
% is omitted, the parameter 'abs_weight' will default to [1 1 1 1 1 1],
% which leads to uniform absorption coefficients for all room boundaries.
%
% The structure SETUP_STRUC may contain one of the two fields 'T60' or
% 'T20'. This function will automatically determine which reverberation
% type is used and compute the desired room absorption coefficients
% accordingly. T20 is defined as the time required for the impulse response
% energy to decay from -5 to -25dB, whereas T60 corresponds to the time
% required by the impulse response energy to decay by 60dB. Setting the
% corresponding field value to 0 achieves anechoic impulse responses
% (direct path only).
%
% In addition, a number of other (optional) parameters can be set using a
% series of 'argument'--value pairs. The following parameters (arguments)
% can be used:
%
%   'Delta_dB': scalar (in dB), parameter determining how much the resulting
%               impulse response is cropped: the impulse response is
%               computed until the time index where its overall energy
%               content has decreased by 'Delta_dB' decibels, after which
%               the computations stop. Not relevant if the reverberation
%               time is set to 0 (anechoic case). Defaults to 50.
% 'SilentFlag': set to 1 to disable this function's on-screen messages.
%               Defaults to 0.
%
% This function returns a 2-dimensional cell array RIR_CELL containing the
% RIRs for each source trajectory point and each microphone, organised as
% follows: RIR_CELL{mic_index,traj_index}. The resulting filter length
% may differ slightly in each computed RIR.
%
% This function also saves the computation results on file. The argument
% RIR_FILE_NAME determines the name of the .mat file where the variable
% RIR_CELL is to be saved. If a file already exists with the same name as
% the input argument, the user will be prompted to determine whether the
% file is to be overwritten or not. The given parameter RIR_FILE_NAME can
% be a full access path to the desired file. If no access path is given,
% the file is saved in the current working directory.

% Release date: November 2009
% Author: Eric A. Lehmann, Perth, Australia (www.eric-lehmann.com)
%
% Copyright (C) 2009 Eric A. Lehmann
%
% 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/>.

VarList = {'Delta_dB'       50;         % maximum attenuation in RIR computation
           'SilentFlag'     0;          % set to 1 to disable on-screen messages
           'useGPU'         0};         % whether to use GPU
eval(SetUserVars(VarList,varargin));    % set user-definable variables

if isempty(RIRFileName)
    DoNotSave = 1;
else
    DoNotSave = 0;
    if length(RIRFileName)<=4 || ~strcmpi(RIRFileName(end-3:end),'.mat'),
        RIRFileName = [RIRFileName '.mat'];
    end
end
if exist(RIRFileName,'file')==2,
    foo = input(' [ISM_RIR_bank] The file name passed as input already exists. Overwrite? [yes/no]: ','s');
    if ~strcmpi(foo,'yes');
        fprintf(' [ISM_RIR_bank] Terminating execution now (no data was saved).\n');
        return
    end
end

Fs = setupstruc.Fs;
room = setupstruc.room;
micpos = setupstruc.mic_pos;
straj = setupstruc.src_traj;

if isfield(setupstruc,'abs_weights'),
    weights = setupstruc.abs_weights;
else
    weights = ones(1,6);
end
if isfield(setupstruc,'c'),
    cc = setupstruc.c;
else
    cc = 343;
end
if isfield(setupstruc,'T60'),
    if isfield(setupstruc, 'reflect_weights')
        beta = setupstruc.reflect_weights;
    else
        alpha = ISM_AbsCoeff('t60',setupstruc.T60,room,weights,'LehmannJohansson','c',cc);
    end
    rttype = 'T60'; rtval = setupstruc.T60;
elseif isfield(setupstruc,'T20'),
    alpha = ISM_AbsCoeff('t20',setupstruc.T20,room,weights,'LehmannJohansson','c',cc);
    rttype = 'T20'; rtval = setupstruc.T20;
else
    error('Missing T60 or T20 field.');
end
if isfield(setupstruc, 'reflect_weights')==0
    beta = sqrt(1-alpha);
end

nMics = size(micpos,1);     % number of microphones
nSPts = size(straj,1);      % number of source trajectory points

%-=:=- Compute RIR bank -=:=-
RIR_cell = cell(nMics,nSPts); % pre-allocate cell array
if ~SilentFlag, PrintLoopPCw(' [ISM_RIR_bank] Computing room impulse responses. '); end;
for mm=1:nMics,
    X_rcv = micpos(mm,:);
    for tt=1:nSPts,         % compute ISM room impulse response for each source-receiver combinations
        if ~SilentFlag, PrintLoopPCw((mm-1)*nSPts+tt,nMics*nSPts); end;
        X_src = straj(tt,:);
        if strcmpi(rttype, 'T60')       % modified by Xiong Xiao: if we use T60, simply set RIR length to T60 length
            RIR_cell{mm,tt} = ISM_RoomResp_GPU(Fs,beta,rttype,rtval,X_src,X_rcv,room,'SilentFlag',1,'c',cc,'Delta_dB',Delta_dB, 'MaxDelay', rtval, 'useGPU', useGPU);
        else
            RIR_cell{mm,tt} = ISM_RoomResp_GPU(Fs,beta,rttype,rtval,X_src,X_rcv,room,'SilentFlag',1,'c',cc,'Delta_dB',Delta_dB, 'useGPU', useGPU);
        end
    end
end

%-=:=- Save results into .mat file -=:=-
if DoNotSave==0
    save(RIRFileName,'RIR_cell');
    if ~SilentFlag, fprintf(' [ISM_RIR_bank] RIR bank parameter ''RIR_cell'' saved in file ''%s''\n',RIRFileName); end;
end

author: singaxiong

repo: SignalGraph


text:
% POP_CROSSF - Return estimates and plots of event-related spectral coherence
%
% Usage:
%       >> pop_crossf(EEG, typeproc, num1, num2, tlimits,cycles,
%                                 'key1',value1,'key2',value2, ... );
% Inputs:
%   INEEG    - Input EEG dataset
%   typeproc - Type of processing:
%                  1 = process two raw-data channels,
%                  0 = process two ICA components
%   num1     - First component or channel number
%   num2     - Second component or channel number
%   tlimits  - [mintime maxtime] Sub-epoch time limits in ms
%   cycles   -   >0 -> Number of cycles in each analysis wavelet
%                 0 -> Use FFTs (with constant window length)
%
% Optional inputs: As for CROSSF.  See >> help crossf
%
% Outputs: Same as CROSSF. No outputs are returned when a
%          window pops-up to ask for additional arguments
%
% Author: Arnaud Delorme, CNL / Salk Institute, 11 March 2002
%
% See also: TIMEF, EEGLAB

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

% 03-18-02 added title -ad & sm
% 04-04-02 added outputs -ad & sm

function varargout = pop_crossf(EEG, typeproc, num1, num2, tlimits, cycles, varargin );

varargout{1} = '';
% display help if not enough arguments
% ------------------------------------
if nargin < 2
	help pop_crossf;
	return;
end;
lastcom = [];
if nargin < 3
	popup = 1;
else
	popup = ischar(num1) | isempty(num1);
	if ischar(num1)
		lastcom = num1;
	end
end

% pop up window
% -------------
if popup
	[txt vars] = gethelpvar('timef.m');

	geometry = { [1 0.5 0.5] [1 0.5 0.5] [1 0.5 0.5] [1 0.5 0.5] [0.92 0.1 0.78] [1 0.5 0.5] [1 0.8 0.2] [1] [1 1]};
    uilist = { { 'Style', 'text', 'string', fastif(typeproc, 'First channel number', 'First component number'), 'fontweight', 'bold'  } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,3,[],'1') } {} ...
			   { 'Style', 'text', 'string', fastif(typeproc, 'Second channel number', 'Second component number'), 'fontweight', 'bold'  } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,4,[],'2') } {} ...
			   { 'Style', 'text', 'string', 'Epoch time range [min max] (msec)', 'fontweight', 'bold', ...
				 'tooltipstring', 'Sub epoch time limits' } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,5,[],[ int2str(EEG.xmin*1000) ' ' int2str(EEG.xmax*1000) ]) } {} ...
			   { 'Style', 'text', 'string', 'Wavelet cycles (0->FFT, see >> help timef)', 'fontweight', 'bold', ...
				 'tooltipstring', context('cycles',vars,txt) } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,6,[],'3 0.5') } {} ...
			   { 'Style', 'text', 'string',  '[set]->Linear coher / [unset]->Phase coher', 'fontweight', 'bold', ...
				 'tooltipstring', ['Compute linear inter-trial coherence (coher)' 10 ...
					'OR Amplitude-normalized inter-trial phase coherence (phasecoher)'] } ...
			   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'phasecoher','present',1) } { } ...
			   { 'Style', 'text', 'string', 'Bootstrap significance level (Ex: 0.01 -> 1%)', 'fontweight', 'bold', ...
				 'tooltipstring', context('alpha',vars,txt) } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,'alpha') } {} ...
			   { 'Style', 'text', 'string', 'Optional timef() arguments (see Help)', 'fontweight', 'bold', ...
				 'tooltipstring', 'See crossf() help via the Help button on the right...' } ...
			   { 'Style', 'edit', 'string', '''padratio'', 4' } ...
			   { 'Style', 'pushbutton', 'string', 'Help', 'callback',  'pophelp(''crossf'');' } ...
			   {} ...
			   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'plotersp','present',0), 'string', ...
				 'Plot coherence amplitude', 'tooltipstring', ...
				 'Plot coherence ampltitude image in the upper panel' } ...
			   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'plotphase','present',0), 'string', ...
				 'Plot coherence phase', 'tooltipstring', ...
				 'Plot coherence phase image in the lower panel' } ...
			 };

	result = inputgui( geometry, uilist, 'pophelp(''pop_crossf'');', ...
					   fastif(typeproc, 'Plot channel cross-coherence -- pop_crossf()', ...
							  'Plot component cross-coherence -- pop_crossf()'));
	if length( result ) == 0 return; end

	num1     = eval( [ '[' result{1} ']' ] );
	num2     = eval( [ '[' result{2} ']' ] );
	tlimits	 = eval( [ '[' result{3} ']' ] );
	cycles	 = eval( [ '[' result{4} ']' ] );
    if result{5}
    	options = [',''type'', ''coher''' ];
    else
		options = [',''type'', ''phasecoher''' ];
    end;

    % add topoplot
    % ------------
	if isfield(EEG.chanlocs, 'theta')
        if ~isfield(EEG, 'chaninfo'), EEG.chaninfo = []; end
		if typeproc == 1
			options = [options ', ''topovec'', [' int2str([num1 num2]) ...
                       '], ''elocs'', EEG.chanlocs, ''chaninfo'', EEG.chaninfo' ];
		else % typeproc == 0
			options = [options ', ''topovec'', EEG.icawinv(:, [' int2str([num1 num2]) ...
                       '])'', ''elocs'', EEG.chanlocs, ''chaninfo'', EEG.chaninfo' ];
		end
	end

    % add title
    % ---------
	if isempty( findstr(  'title', result{7}))
        if ~isempty(EEG.chanlocs) && typeproc
            chanlabel1 = EEG.chanlocs(num1).labels;
            chanlabel2 = EEG.chanlocs(num2).labels;
        else
            chanlabel1 = int2str(num1);
            chanlabel2 = int2str(num2);
        end
		if result{5}
            options = [options ', ''title'',' fastif(typeproc, '''Channel ', '''Component ') chanlabel1 '-' chanlabel2 ...
					' Coherence'''];
        else
            options = [options ', ''title'',' fastif(typeproc, '''Channel ', '''Component ') chanlabel1 '-' chanlabel2 ...
					' Phase Coherence''' ];
		end
	end
	if ~isempty( result{6} )
		options      = [ options ', ''alpha'',' result{6} ];
	end
	if ~isempty( result{7} )
		  options = [ options ',' result{7} ];
	end
	if ~result{8}
		options = [ options ', ''plotersp'', ''off''' ];
	end
	if ~result{9}
		options = [ options ', ''plotphase'', ''off''' ];
	end
    figure; try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end;
else
	options = [ ',' vararg2str(varargin) ];
end

% compute epoch limits
% --------------------
if isempty(tlimits)
	tlimits = [EEG.xmin, EEG.xmax]*1000;
end;
pointrange1 = round(max((tlimits(1)/1000-EEG.xmin)*EEG.srate, 1));
pointrange2 = round(min((tlimits(2)/1000-EEG.xmin)*EEG.srate, EEG.pnts));
pointrange = [pointrange1:pointrange2];

% call function sample either on raw data or ICA data
% ---------------------------------------------------
if typeproc == 1
	tmpsig1 = EEG.data(num1,pointrange,:);
	tmpsig2 = EEG.data(num2,pointrange,:);
else
	if ~isempty( EEG.icasphere )
        tmpsig1 = eeg_getdatact(EEG, 'component', num1, 'samples', pointrange);
        tmpsig2 = eeg_getdatact(EEG, 'component', num2, 'samples', pointrange);
	else
		error('You must run ICA first');
	end;
end;
tmpsig1 = reshape( tmpsig1, 1, size(tmpsig1,2)*size(tmpsig1,3));
tmpsig2 = reshape( tmpsig2, 1, size(tmpsig2,2)*size(tmpsig2,3));

% outputs
% -------
outstr = '';
if ~popup
    for io = 1:nargout, outstr = [outstr 'varargout{' int2str(io) '},' ]; end
    if ~isempty(outstr), outstr = [ '[' outstr(1:end-1) '] =' ]; end
end

% plot the datas and generate output command
% --------------------------------------------
if length( options ) < 2
    options = '';
end
varargout{1} = sprintf('figure; pop_crossf( EEG, %d, %d, %d, [%s], [%s] %s);', ...
          typeproc, num1, num2, int2str(tlimits), num2str(cycles), options);

%options = [ options ', ''ydir'', ''norm''' ];
com = sprintf( '%s crossf( tmpsig1, tmpsig2, length(pointrange), [tlimits(1) tlimits(2)], EEG.srate, cycles %s);', outstr, options);
eval(com)

return;

% get contextual help
% -------------------
function txt = context(var, allvars, alltext);
	loc = strmatch( var, allvars);
	if ~isempty(loc)
		txt= alltext{loc(1)};
	else
		disp([ 'warning: variable ''' var ''' not found']);
		txt = '';
	end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function a = cheby_t_inverse ( n )

%*****************************************************************************80
%
%% CHEBY_T_INVERSE returns the inverse of the CHEBY_T matrix.
%
%
%  Example:
%
%    N = 11
%
%      1   .   .  .   .  .  .  .  .  .  .
%      .   1   .  .   .  .  .  .  .  .  .
%      1   .   1  .   .  .  .  .  .  .  .  /   2
%      .   3   .  1   .  .  .  .  .  .  .  /   4
%      3   .   4  .   1  .  .  .  .  .  .  /   8
%      .  10   .  5   .  1  .  .  .  .  .  /  16
%     10   .  15  .   6  .  1  .  .  .  .  /  32
%      .  35   . 21   .  7  .  1  .  .  .  /  64
%     35   .  56  .  28  .  8  .  1  .  .  / 128
%      . 126   . 84   . 36  .  9  .  1  .  / 256
%    126   . 210  . 120  . 45  . 10  .  1  / 512
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 20007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  a(1,1) = 1.0;

  if ( n == 1 )
    return
  end

  a(2,2) = 1.0;

  if ( n == 2 )
    return
  end

  for i = 3 : n
    for j = 1 : n
      if ( j == 1 )
        a(i,j) =                      a(i-1,j+1)   / 2.0;
      elseif ( j == 2 )
        a(i,j) = ( 2.0 * a(i-1,j-1) + a(i-1,j+1) ) / 2.0;
      elseif ( j < n )
        a(i,j) = (       a(i-1,j-1) + a(i-1,j+1) ) / 2.0;
      else
        a(i,j) =         a(i-1,j-1)                / 2.0;
      end
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Animate_Double_Pendulum(t,K,P)

%FUNCTION:
%   This function is used to generate an animation of the double pendulum,
%   using precomputed values.
%
%INPUTS:
%   t = a vector of time, (tout, returned by ode45)
%   K = Kinematics struct, returned by Double_Pendulum_Kinematics.m
%   P = Parameters struct, returned by Set_Parameters.m
%
%OUTPUTS:
%   Figure 1  --  Animation of double pendulum
%

%Which figure to plopt on?
figure(1);

%Get the various position vectors
P1_x = K.r_P1_O(1,:);
P1_y = K.r_P1_O(2,:);
P2_x = K.r_P2_O(1,:);
P2_y = K.r_P2_O(2,:);

%Get the font sizes
    TitleFontSize = P.plot.TitleFontSize;
    LabelFontSize = P.plot.LabelFontSize;

%Display a few things:
Parameters_to_Display = {...
        ['Slow Motion Factor: ' num2str(P.sim.slow_motion_factor)];...
        '';...
        ['Mass 1:  ' num2str(P.dyn.m1,4) '  kg'];...
        ['Mass 2:  ' num2str(P.dyn.m2,4) '  kg'];...
        ['Gravity:  ' num2str(P.dyn.g,4) '  m/s^2'];...
        ['Length:  ' num2str(P.dyn.L,4) '  m'];...
};

tic;    %Start a timer
Loop_Time = 0;    %store how long has the simulation been running
i=2;   %Start at second index (for interpolation purposes)
Max_i = length(t);
T_end = t(end);   %Ending time of one step
L = P.dyn.L;   %Leg length
Bounds = 1.1 * 2*L * [-1,1,-1,1];

while Loop_Time < T_end;  %Loop while the CPU time is less than the end of the simulation's time
    %The next few lines pull the time step that is closest to the real time
    Loop_Time = toc/P.sim.slow_motion_factor;   %Get the current time  (Taking slow motion into accunt if desired)
    while (i<Max_i) && (Loop_Time > (t(i))) %While we are not at the end of the data and the CPU time is ahead of the simulation time
        i=i+1;   %Go to the next time frame
    end
    %Now t(i-1) < Loop_Time < t(i) should be true

    %Linear interpolation scheme
        t_m = t(i-1);
        t_p = t(i);

        Pend_X_m = [0;P1_x(i-1); P2_x(i-1)];
        Pend_Y_m = [0;P1_y(i-1); P2_y(i-1)];

        Pend_X_p = [0;P1_x(i); P2_x(i)];
        Pend_Y_p = [0;P1_y(i); P2_y(i)];

        dPend_X = Pend_X_p - Pend_X_m;
        dPend_Y = Pend_Y_p - Pend_Y_m;

        dt = t_p - t_m;
        dt_cpu = Loop_Time - t_m;
        m = dt_cpu/dt;

        Pend_X = m*dPend_X + Pend_X_p;
        Pend_Y = m*dPend_Y + Pend_Y_p;

    clf
    plot(Pend_X,Pend_Y,'k-','LineWidth',6)
    hold on
    plot(Pend_X,Pend_Y,'b.','MarkerSize',50)
    plot(0,0,'ko','MarkerSize',40,'LineWidth',3)
    title(['Simulation Time: ' num2str(t(i)) ' seconds'],'FontSize',TitleFontSize)

    text(Bounds(2),mean(Bounds(3:4)),Parameters_to_Display,'FontSize',LabelFontSize);

    axis(Bounds); axis equal, axis manual; axis off;
    drawnow

end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [TB] = bit2TB(bit)
% Convert computery things from bits to terabytes.
% Chad A. Greene 2012
TB = bit*2^-43;


text:
function assertNotEqual(a, b, message)
% assert that two inputs are not equal
%
% assertNotEqual(a,b,[msg])
%
% Inputs:
%   a           first input  } of any
%   b           second input } type
%   msg             optional custom message
%
% Raises:
%   'assertNotEqual:equal'                a and b are equal, and are either
%                                         both sparse or both not sparse
%
% Examples:
%   assertNotEqual('foo','bar');
%   %|| % passes without output
%
%   assertNotEqual(1,sparse(1));
%   %|| % passes without output
%
%   assertNotEqual('foo','foo');
%   %|| error('elements are equal');
%
% Notes:
%   - If a custom message is provided, then any error message is prefixed
%     by this custom message
%   - In this function, NaN values are considered equal.
%   - This function attempts to show similar behaviour as in
%     Steve Eddins' MATLAB xUnit Test Framework (2009-2012)
%     URL: http://www.mathworks.com/matlabcentral/fileexchange/
%                           22846-matlab-xunit-test-framework
%

    % Note: although it may seem more logical to compare class before size,
    % for compatibility reasons the order of tests matches that of the
    % MATLAB xUnit framework

    if isequaln_wrapper(a,b) && issparse(a)==issparse(b)
        whatswrong='inputs are equal';
        error_id='assertNotEqual:equal';

    else
        % elements are equal
        return;
    end

    if nargin<3
        message='';
    end

    full_message=moxunit_util_input2str(message,whatswrong,a,b);

    if moxunit_util_platform_is_octave()
        error(error_id,'%s',full_message);
    else
        throwAsCaller(MException(error_id, '%s', full_message));
    end

function tf=isequaln_wrapper(a,b)
% wrapper to support old versions of Matlab
    persistent has_equaln;

    if isempty(has_equaln)
        has_equaln=~isempty(which('isequaln'));
    end

    if has_equaln
        tf=isequaln(a,b);
    else
        tf=isequalwithequalnans(a,b);
    end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% CAMCALD Camera calibration from data points
%
% C = CAMCALD(D) is the camera matrix (3x4) determined by least squares
% from corresponding world and image-plane points.  D is a table
% of points with rows of the form [X Y Z U V] where (X,Y,Z) is the
% coordinate of a world point and [U,V] is the corresponding image
% plane coordinate.
%
% [C,E] = CAMCALD(D) as above but E is the maximum residual error after
% back substitution [pixels].
%
% Notes:
% - This method assumes no lense distortion affecting the image plane
%   coordinates.
%
% See also CentralCamera.

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

function [C, resid] = camcald(XYZ, uv)

    if numcols(XYZ) ~= numcols(uv)
        error('must have same number of world and image-plane points');
    end

    n = numcols(XYZ);
    if n < 6,
        error('At least 6 points required for calibration');
    end
%
% build the matrix as per Ballard and Brown p.482
%
% the row pair are one row at this point
%
    A = [ XYZ' ones(n,1) zeros(n,4) -repmat(uv(1,:)', 1,3).*XYZ' ...
          zeros(n,4) XYZ' ones(n,1) -repmat(uv(2,:)', 1,3).*XYZ'  ];
%
% reshape the matrix, so that the rows interleave
%
    A = reshape(A',11, n*2)';
    if rank(A) < 11,
        error('Rank deficient,  perhaps points are coplanar or collinear');
    end

    B = reshape( uv, 1, n*2)';

    C = A\B;    % least squares solution
    resid = max(max(abs(A * C - B)));
    if resid > 1,
        warning('Residual greater than 1 pixel');
    end
    fprintf('maxm residual %f pixels.\n', resid);
    C = reshape([C;1]',4,3)';

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [point,sc] = harris_strongest(im,sigma,mrg,edgerej)
% HARRIS_STRONGEST Extract strongest Harris point.
%   HARRIS_STRONGEST(IM) extracts the strongest harris point in image IM.
%
%   HARRIS_STRONGEST(IM, SG, MRG, EDGEREJ) accepts additional inputs as
%   follows:
%
%   INPUT
%   =====
%   im     : the graylevel image
%   sigma  : STD of the smoothing mask in pixels
%   mrg    : inner margin to be ignored
%   edgerej: maximum ratio allowed between largest and smallest eigenvalues
%
%   [P,SC] = HARRIS_STRONGEST(...) returns also the point's strength as a
%   scalar score (non-normalized).
%
%   OUTPUT
%   ======
%   point  : the interest point extracted
%   sc     : strength of corner point
%
%   EXAMPLE
%   =======
%   [point,sc] = harris_strongest(im)

% Author :: Vincent Garcia - multiple point detector
% Date   :: 05/12/2007
%
% Author :: Joan Sola - strongest point detector, margin ignoring, edge
% rejection
% Date   :: 01/08/2008
%
% REFERENCES
% ==========
% C.G. Harris and M.J. Stephens. "A combined corner and edge detector",
% Proceedings Fourth Alvey Vision Conference, Manchester.
% pp 147-151, 1988.
%
% Alison Noble, "Descriptions of Image Surfaces", PhD thesis, Department
% of Engineering Science, Oxford University 1989, p45.
%
% C. Schmid, R. Mohrand and C. Bauckhage, "Evaluation of Interest Point Detectors",
% Int. Journal of Computer Vision, 37(2), 151-172, 2000.
%

persistent x dx dy ex g; % gaussian filter

% Input options management
switch nargin
    case 1
        sigma = 2;
        mrg = 0;
        edgerej = 20;
    case 2
        mrg = 0;
        edgerej = 20;
    case 3
        edgerej = 20;
end

% only luminance value
im = single(im(:,:,1));

% image size
[sx,sy] = size(im);

% derivative masks
s_D = 0.5*sigma;
if isempty(x)
    x  = -round(2.0*s_D):round(2.0*s_D);
    ex = exp(-x.*x/(2*s_D^2)) ./ (s_D*sqrt(2*pi));
    dx = s_D^2 * x .* ex;
    dy = dx';
    g = conv2(ex,ex');
    g = g/sum(sum(g));
end

% image derivatives
Ix = conv2(im, dx, 'same');
Iy = conv2(im, dy, 'same');

Ix2 = conv2(Ix.^2, g, 'same'); % Smoothed squared image derivatives
Iy2 = conv2(Iy.^2, g, 'same');
Ixy = conv2(Ix.*Iy, g, 'same');

% interest point response
% Original Harris measure.
% k = 0.06; cim = (Ix2.*Iy2 - Ixy.^2) - k*(Ix2 + Iy2).^2;

% Alison Noble measure.
% cim = (Ix2.*Iy2 - Ixy.^2)./(Ix2 + Iy2 + eps);

% Shi and Tomasi -- smallest eigenvalue
a = Ix2+Iy2;
b = sqrt((Ix2-Iy2).^2+4*Ixy.^2);
cim = a-b;  % Smallest eigenvalue.

% supress corner measures close to margin
cim([1:mrg sx-mrg+1:sx],:) = -9e9;
cim(:,[1:mrg sy-mrg+1:sy]) = -9e9;

% Find strongest corner [u;v] and score
[scv,v] = max(cim);
[sc,u]  = max(scv);     % score and u-coordinate
v       = v(u);         % v-coordinate

% reject edges
l = sc;                 % smallest eigenvalue is the score
L = a(v,u) + b(v,u);    %  largest eigenvalue
if L/l > edgerej        % reject edge detections
    sc = 0;             % return nothing if winner is an edge
end

point = [u;v];

% comment out this part for displaying purposes
% figure(98)
% imshow(im,[0 255]);colormap(gray(255))
% axis on
% hold on;
% sz = 3;
% rectangle('Position',[u-sz,v-sz,2*sz,2*sz],'Curvature',[0 0],'EdgeColor','b','LineWidth',1);
% figure(99)
% imshow(cim,[0 max(max(cim))]); axis on;colormap(gray(255))
% rectangle('Position',[u-sz,v-sz,2*sz,2*sz],'Curvature',[0 0],'EdgeColor','b','LineWidth',1);
% point;
% sc;

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [dSP,sp]=grShortPath(E,s,t)
% Function dSP=grShortPath(E) solve the task about
% the shortest path between any vertexes of digraph.
% Input parameter:
%   E(m,2) or (m,3) - the arrows of digraph and their weight;
%     1st and 2nd elements of each row is numbers of vertexes;
%     3rd elements of each row is weight of arrow;
%     m - number of arrows.
%     If we set the array E(m,2), then all weights is 1.
% Output parameter:
%   dSP(n,n) - the matrix of shortest path.
%   Each element dSP(i,j) is the shortest path
%   from vertex i to vertex j (may be inf,
%   if vertex j is not accessible from vertex i).
% Uses the algorithm of R.W.Floyd, S.A.Warshall.
% [dSP,sp]=grShortPath(E,s,t) - find also
% the shortest paths from vertex s (source) to vertex t (tail).
% In this case output parameter sp is vector with numbers
% of vertexes, included to shortest path from s to t.
% Author: Sergiy Iglin
% e-mail: siglin@yandex.ru
% personal page: http://iglin.exponenta.ru
% Acknowledgements to Prof. Gerard Biau
% (Universite Montpellier II, France)
% for testing of this algorithm.

% ============= Input data validation ==================
if nargin<1,
  error('There are no input data!')
end
[m,n,E] = grValidation(E); % E data validation

% ================ Initial values ===============
dSP=ones(n)*inf; % initial distances
dSP((E(:,2)-1)*n+E(:,1))=E(:,3);

% ========= The main cycle of Floyd-Warshall algorithm =========
for j=1:n,
  i=setdiff((1:n),j);
  dSP(i,i)=min(dSP(i,i),repmat(dSP(i,j),1,n-1)+repmat(dSP(j,i),n-1,1));
end
sp=[];
if (nargin<3)|(isempty(s))|(isempty(t)),
  return
end
s=s(1);
t=t(1);
if (~(s==round(s)))|(~(t==round(t)))|(s<1)|(s>n)|(t<1)|(t>n),
  error(['s and t must be integer from 1 to ' num2str(n)])
end
if isinf(dSP(s,t)), % t is not accessible from s
  return
end
dSP1=dSP;
dSP1(1:n+1:n^2)=0; % modified dSP
l=ones(m,1); % label for each arrow
sp=t; % final vertex
while ~(sp(1)==s),
  nv=find((E(:,2)==sp(1))&l); % all labeled arrows to sp(1)
  vnv=abs((dSP1(s,sp(1))-dSP1(s,E(nv,1)))'-E(nv,3))<eps*1e3; % valided arrows
  l(nv(~vnv))=0; % labels of not valided arrows
  if all(~vnv), % invalided arrows
    l(find((E(:,1)==sp(1))&(E(:,2)==sp(2))))=0;
    sp=sp(2:end); % one step back
  else
    nv=nv(vnv); % rested vaded arrows
    sp=[E(nv(1),1) sp]; % add one vertex to shortest path
  end
end
return

author: HuangCongQing


text:
function [hFig, iDS, iFig] = view_timeseries(DataFile, Modality, RowNames, hFig)
% VIEW_TIMESERIES: Display times series in a new figure.
%
% USAGE: [hFig, iDS, iFig] = view_timeseries(DataFile, Modality=[], RowNames=[], hFig=[])
%        [hFig, iDS, iFig] = view_timeseries(DataFile, Modality=[], RowNames=[], 'NewFigure')
%        [hFig, iDS, iFig] = view_timeseries(DataFiles, ...)
%
% INPUT:
%     - DataFile  : Path to data file to visualize
%     - DataFiles : Cell-array of file paths (simply calls this function recursively for each file)
%     - Modality  : Modality to display with the input Data file
%     - RowNames  : Cell array of channel names to plot in this figure
%     - "NewFigure" : force new figure creation (do not re-use a previously created figure)
%     - hFig        : Specify the figure in which to display the MRI
%
% OUTPUT :
%     - hFig : Matlab handle to the 3DViz figure that was created or updated
%     - iDS  : DataSet index in the GlobalData variable
%     - iFig : Indice of returned figure in the GlobalData(iDS).Figure array

% @=============================================================================
% 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-2023

%% ===== INITIALIZATION =====
global GlobalData;
% Parse inputs
if (nargin < 3) || isempty(RowNames)
    RowNames = {};
elseif ischar(RowNames)
    RowNames = {RowNames};
end
if (nargin < 2) || isempty(Modality)
    % Get default modality
    [AllMod,DispMod,Modality] = bst_get('ChannelModalities', DataFile);
    % Replace SEEG and ECOG with SEEG+ECOG
    if ~isempty(DispMod) && all(ismember({'SEEG','ECOG'}, DispMod))
        DispMod = cat(2, {'ECOG+SEEG'}, setdiff(DispMod, {'SEEG','ECOG'}));
        if ismember(Modality, {'SEEG','ECOG'})
            Modality = 'ECOG+SEEG';
        end
    elseif ~isempty(AllMod) && ismember(Modality, {'SEEG','ECOG'}) && all(ismember({'SEEG','ECOG'}, AllMod))
        Modality = 'ECOG+SEEG';
        DispMod = union(DispMod, {'ECOG+SEEG','SEEG','ECOG'});
    end
else
    DispMod = [];
end
% Get target figure
if (nargin < 4) || isempty(hFig)
    hFig = [];
    iFig = [];
    NewFigure = 0;
elseif ischar(hFig) && strcmpi(hFig, 'NewFigure')
    hFig = [];
    iFig = [];
    NewFigure = 1;
elseif ishandle(hFig)
    [hFig,iFig,iDS] = bst_figures('GetFigure', hFig);
    NewFigure = 0;
else
    error('Invalid figure handle.');
end
% Multiple files in input: call function recursively
if iscell(DataFile)
    hFigInput = hFig;
    hFig = zeros(1,length(DataFile));
    iDS = zeros(1,length(DataFile));
    iFig = zeros(1,length(DataFile));
    for iFile = 1:length(DataFile)
        [hFig(iFile), iDS(iFile), iFig(iFile)] = view_timeseries(DataFile{iFile}, Modality, RowNames, hFigInput);
        % Remove the borders in the figure
        hAxes = findobj(hFig(iFile), '-depth', 1, 'tag', 'AxesGraph');
        set(hFig(iFile), bst_get('ResizeFunction'), []);
        set(hAxes, 'Units', 'normalized', 'Position', [0 0 1 1]);
    end
    return;
end

%% ===== GET A DATASET AND LOAD DATA =====
% Get DataFile information
[sStudy, iData, ChannelFile] = bst_memory('GetFileInfo', DataFile);
% If Channel is not defined
if isempty(ChannelFile)
    Modality = 'EEG';
end
% If not loaded yet
if isempty(hFig)
    % Load file
    iDS = bst_memory('LoadDataFile', DataFile);
    % If no DataSet is accessible : error
    if isempty(iDS)
        return
    end
end
% Check that the selected modality can be displayed, if not select another one
if ~isempty(DispMod)
    % List of expected modalities, in order of preference
    AllMod = {'MEG', 'MEG MAG', 'MEG GRAD', 'EEG', 'ECOG', 'SEEG', 'ECOG+SEEG'};
    % Get the list of good modalities
    GoodMod = unique({GlobalData.DataSet(iDS).Channel(GlobalData.DataSet(iDS).Measures.ChannelFlag == 1).Type});
    GoodMod = intersect(GoodMod, DispMod);
    % Add MEG for Elekta systems
    if any(ismember({'MEG MAG', 'MEG GRAD'}, GoodMod))
        GoodMod{end+1} = 'MEG';
    end
    % Add combined ECOG+SEEG
    if all(ismember({'ECOG', 'SEEG'}, GoodMod))
        GoodMod{end+1} = 'ECOG+SEEG';
    end
    % Get the preferred modalities
    iMod = find(ismember(AllMod, GoodMod));
    GoodMod = AllMod(iMod);
    % If the selected modality is not in the list of the good ones: select the first of the good ones.
    if ~isempty(iMod) && ~ismember(Modality, GoodMod)
        Modality = GoodMod{1};
    end
end

%% ===== CREATE A NEW FIGURE =====
bst_progress('start', 'View time series', 'Loading data...');
if isempty(hFig)
    % Prepare FigureId structure
    FigureId.Type     = 'DataTimeSeries';
    FigureId.SubType  = '';
    FigureId.Modality = Modality;
    % Create TimeSeries figure
    if NewFigure
        [hFig, iFig, isNewFig] = bst_figures('CreateFigure', iDS, FigureId, 'AlwaysCreate', RowNames);
    else
        [hFig, iFig, isNewFig] = bst_figures('CreateFigure', iDS, FigureId, [], RowNames);
    end
    if isempty(hFig)
        bst_error('Could not create figure', 'View time series', 0);
        return;
    end
else
    isNewFig = 0;
end
% Add DataFile to figure appdata
setappdata(hFig, 'DataFile', DataFile);
setappdata(hFig, 'StudyFile', GlobalData.DataSet(iDS).StudyFile);
setappdata(hFig, 'SubjectFile', GlobalData.DataSet(iDS).SubjectFile);

%% ===== SELECT ROWS =====
LinesColor = {};
% Select only the channels that we need to plot
if ~isempty(RowNames)
    % Get the channels normally displayed in this figure
    iSelChanMod = GlobalData.DataSet(iDS).Figure(iFig).SelectedChannels;
    AllChannels = {GlobalData.DataSet(iDS).Channel.Name};
    % Remove all spaces
    AllChannels = cellfun(@(c)strrep(c,' ',''), AllChannels, 'UniformOutput', 0);
    RowNames    = cellfun(@(c)strrep(c,' ',''), RowNames,    'UniformOutput', 0);
    % Get the channels that are requested from the command line call (RowNames argument)
    iSelChanCall = [];
    for i = 1:length(RowNames)
        iSelChanCall = [iSelChanCall, find(strcmpi(RowNames{i}, AllChannels))];
        % NIRS: color based on the name
        if isequal(Modality, 'NIRS') && (length(iSelChanCall) == i)
            if ~isempty(strfind(lower(RowNames{i}), 'hbo'))
                LinesColor{i} = [1,0,0];
            elseif ~isempty(strfind(lower(RowNames{i}), 'hbr'))
                LinesColor{i} = [0,0,1];
            elseif ~isempty(strfind(lower(RowNames{i}), 'hbt'))
                LinesColor{i} = [0,1,0];
            end
        else
            LinesColor = {};
        end
    end
    % Keep only the intersection of the two selections (if non-empty)
    if ~isempty(iSelChanCall) && ~isempty(intersect(iSelChanMod, iSelChanCall))
        GlobalData.DataSet(iDS).Figure(iFig).SelectedChannels = intersect(iSelChanMod, iSelChanCall);
    end
    % Redraw position of figures
    if (iFig > 1)
        gui_layout('Update');
    end
end

%% ===== CONFIGURE FIGURE =====
% Static dataset ?
setappdata(hFig, 'isStatic', (GlobalData.DataSet(iDS).Measures.NumberOfSamples <= 2));
% Raw file
isRaw = strcmpi(GlobalData.DataSet(iDS).Measures.DataType, 'raw');
% Statistics?
% isStat = strcmpi(GlobalData.DataSet(iDS).Measures.DataType, 'stat');
isStat = ~ismember(GlobalData.DataSet(iDS).Measures.DataType, {'recordings', 'raw'});

% Create time-series information structure
if isNewFig
    % Create figure structure
    TsInfo = db_template('TsInfo');
    TsInfo.FileName      = DataFile;
    TsInfo.Modality      = Modality;
    TsInfo.DisplayMode   = bst_get('TSDisplayMode');
    TsInfo.LinesLabels   = {};
    TsInfo.AxesLabels    = {};
    TsInfo.LinesColor    = LinesColor;
    TsInfo.RowNames      = RowNames;
    TsInfo.MontageName   = [];
    TsInfo.DefaultFactor = figure_timeseries('GetDefaultFactor', Modality);
    TsInfo.DisplayUnits  = GlobalData.DataSet(iDS).Measures.DisplayUnits;
    TsInfo.FlipYAxis     = ~isempty(Modality) && ismember(Modality, {'EEG','MEG','MEG GRAD','MEG MAG','SEEG','ECOG','NIRS'}) && ~isStat && bst_get('FlipYAxis');
    TsInfo.AutoScaleY    = bst_get('AutoScaleY');
    TsInfo.NormalizeAmp  = 0;
    TsInfo.Resolution    = [0 0];
    TsInfo.ShowXGrid     = bst_get('ShowXGrid');
    TsInfo.ShowYGrid     = bst_get('ShowYGrid');
    TsInfo.ShowZeroLines = bst_get('ShowZeroLines');
    TsInfo.ShowEventsMode = bst_get('ShowEventsMode');
    % Hide events only for multiple RAW figures
    allId = [GlobalData.DataSet(iDS).Figure(1:iFig-1).Id];
    TsInfo.ShowEvents = ~isRaw || ((iFig == 1) || ~any(strcmpi({allId.Type}, 'DataTimeSeries')));
else
    TsInfo = getappdata(hFig, 'TsInfo');
    TsInfo.FileName = DataFile;
end
setappdata(hFig, 'TsInfo', TsInfo);
% Get default montage
if isNewFig
    sMontage = panel_montage('GetCurrentMontage', Modality);
    % If displaying a SEEG for which a bipolar montage has already been applied: ignore current montage
    if ~isempty(sMontage) && isempty(RowNames) && strcmpi(Modality, 'SEEG') && all(cellfun(@(c)any(c=='-'), {GlobalData.DataSet(iDS).Channel.Name}))
        TsInfo.MontageName = [];
    % Use previous montage
    elseif ~isempty(sMontage) && isempty(RowNames) && (~isStat || strcmpi(sMontage.Type, 'selection')) && ~ismember(sMontage.Name, {'ICA components[tmp]', 'SSP components[tmp]'})
        TsInfo.MontageName = sMontage.Name;
    else
        TsInfo.MontageName = [];
    end
    setappdata(hFig, 'TsInfo', TsInfo);
end
% If no bad channels are available: do not accept bad channels montage
if isequal(TsInfo.MontageName, 'Bad channels') && ~any(GlobalData.DataSet(iDS).Measures.ChannelFlag == -1)
    TsInfo.MontageName = [];
    setappdata(hFig, 'TsInfo', TsInfo);
end
% Default display mode: Force to 'butterfly' if there is only one channel to display
if isequal(TsInfo.MontageName, 'Bad channels')
    nChan = nnz(GlobalData.DataSet(iDS).Measures.ChannelFlag == -1);
else
    nChan = length(GlobalData.DataSet(iDS).Figure(iFig).SelectedChannels);
end
if (nChan == 1) && strcmpi(TsInfo.DisplayMode, 'column')
    TsInfo.DisplayMode = 'butterfly';
    setappdata(hFig, 'TsInfo', TsInfo);
end

%% ===== PLOT TIME SERIES =====
% Update figure selection
bst_figures('SetCurrentFigure', hFig, '2D');
% Plot figure
isOk = figure_timeseries('PlotFigure', iDS, iFig);
% Some error occured during the display procedure: close the window and stop process
if ~isOk
    close(hFig);
    return;
end

%% ===== UPDATE ENVIRONMENT =====
% Uniformize time series scales if required
isUniform = bst_get('UniformizeTimeSeriesScales');
if ~isempty(isUniform) && (isUniform == 1)
    figure_timeseries('UniformizeTimeSeriesScales', 1);
end
% Set the Y scale factor if necessary
if isRaw && ~isempty(Modality) && (Modality(1) ~= '$')
    newScale = bst_get('FixedScaleY', [Modality, TsInfo.DisplayMode]);
    if ~isempty(newScale)
        figure_timeseries('SetScaleY', iDS, iFig, newScale);
    end
end
% Set figure visible
if strcmpi(get(hFig,'Visible'), 'off')
    set(hFig, 'Visible', 'on');
end
% Set the time label visible
figure_timeseries('SetTimeVisible', hFig, 1);
% Last updates
if isNewFig
    % Update record tab
    panel_record('UpdateDisplayOptions', hFig);
    % Update stat tab
    if isStat
        panel_stat('CurrentFigureChanged_Callback', hFig);
    % Select Record tab
    else
        gui_brainstorm('SetSelectedTab', 'Record');
    end
else
    % Update figure name
    bst_figures('UpdateFigureName', hFig);
end
% Open Cluster tab
if ~isempty(GlobalData.DataSet(iDS).Clusters)
    gui_brainstorm('ShowToolTab', 'Cluster');
end
% Close progress bar
drawnow;
bst_progress('stop');

author: brainstorm-tools


text:
function layer2colors=riverplot_recolor_layer2(sim_matrix,layer1colors)

% ..
%    Programmers' notes:
%    List dates and changes here, and author of changes
%    Created July 2016 by Tor Wager
%
%    8/21/2017 Stephan Geuter
%    changed normalization to used absolute similarity values (l21). Before, negative
%    values could distort the normalization factor resulting RGB values
%    outside the [0-1] interval.
% ..

[n2, n1]=size(sim_matrix);
for i=1:n2
    temp_color=[0 0 0];
    for j=1:n1
        temp_color= temp_color + layer1colors{j}*abs(sim_matrix(i,j))/sum(abs(sim_matrix(i,:)));
    end

    if any(isnan(temp_color));
       temp_color=[.5 .5 .5];
    end
    layer2colors{i}=temp_color;
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function out = echo_enc_mirror(signal, text, d0, d1, alpha, L)
%ECHO_ENC_MIRRORED Echo Hiding with Mirrored Echo Kernels
%
%   INPUT VARIABLES
%       signal : Cover signal
%       text   : Message to hide
%       d0     : Delay rate for bit0
%       d1     : Delay rate for bit1
%       alpha  : Echo amplitude
%       L      : Length of frames
%
%   OUTPUT VARIABLES
%       out    : Stego signal
%
%   Kadir Tekeli (kadir.tekeli@outlook.com)

if nargin < 4
	d0 = 150;     %Delay rate for bit0
	d1 = 200;     %Delay rate for bit1
end

if nargin < 5
	alpha = 0.5;  %Echo amplitude
end

if nargin < 6
	L = 8*1024;   %Length of frames
end

[s.len, s.ch] = size(signal);
bit = getBits(text);
nframe = floor(s.len/L);
N = nframe - mod(nframe,8);      %Number of frames (for 8 bit)

if (length(bit) > N)
	warning('Message is too long, being cropped!');
	bits = bit(1:N);
else
	warning('Message is being zero padded...');
	bits = [bit, num2str(zeros(N-length(bit), 1))'];
end

[echo_zro,echo_one] = mirror_echo(signal,d0,d1,alpha);  %Echo signals
mix = mixer(L,bits,0,1,256)*ones(1, s.ch);              %Mixer signal

%%%%%%%%%%%%%%%%%%%%%%% EMBEDDING MESSAGE... %%%%%%%%%%%%%%%%%%%%%%%
out = signal(1:N*L, :) + echo_zro(1:N*L, :) .* abs(mix-1) ...
                       + echo_one(1:N*L, :) .* mix;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
out = [out; signal(N*L+1:s.len, :)];   %Rest of the signal
end

author: ktekeli

repo: audio-steganography-algorithms


text:
function [c,mx,my] = covc(x, y, dim)

if nargin==2,
  dim = y;
  y   = x;
end

if dim==1,
  c  = x'*y;
elseif dim==2,
  c = x*y';
end

mx = sum(x,dim);
my = sum(y,dim);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function D = kn2sd(K)
% Transform a kernel matrix (or inner product matrix) to a squared distance matrix
% Input:
%   K: n x n kernel matrix
% Ouput:
%   D: n x n squared distance matrix
% Written by Mo Chen (sth4nth@gmail.com).
d = diag(K);
D = -2*K+bsxfun(@plus,d,d');

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
%ANIMATE Create an animation
%
% Helper class for creating animations.  Saves snapshots of a figture as a
% folder of individual PNG format frames numbered 0000.png, 0001.png and so
% on.
%
% Example::
%
%          anim = Animate('movie');
%
%          for i=1:100
%              plot(...);
%              anim.add();
%          end
%
% To convert the image files to a movie you could use a tool like ffmpeg
%           % ffmpeg -r 10 -i movie/*.png out.mp4

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

classdef Animate < handle
    properties
        frame
        dir
        resolution
    end

    methods
        function a = Animate(name, res)
            %ANIMATE.ANIMATE Create an animation class
            %
            % A = ANIMATE(NAME, OPTIONS) initializes an animation, and creates a folder
            % called NAME to hold the individual frames.
            %
            % Options::
            % 'resolution',R    Set the resolution of the saved image to R pixels per
            % inch.

            a.frame = 0;
            a.dir = name;
            mkdir(name);
            if nargin > 1
                a.resolution = res;
            else
                a.resolution = [];
            end
            delete( fullfile(name, '*.png') );

        end

        function add(a, fh)
            %ANIMATE.ADD Adds current plot to the animation
            %
            % A.ADD() adds the current figure in PNG format to the animation
            % folder with a unique sequential filename.
            %
            % A.ADD(FIG) as above but captures the figure FIG.
            %
            % See also print.

            if nargin < 2
                fh = gcf;
            end

            if isempty(a.resolution)
                print(fh, '-dpng', fullfile(a.dir, sprintf('%04d.png', a.frame)));
            else
                print(fh, '-dpng', sprintf('-r%d', a.resolution), fullfile(a.dir, sprintf('%04d.png', a.frame)));
            end
            a.frame = a.frame + 1;
        end

        function close(a)
        end
    end
end

author: Allopart

repo: rbpf-gmapping


text:
% Uncomment below if "cannot figure out class error"
% http://in.mathworks.com/matlabcentral/answers/137332-error-using-matlab-compiler-app-applicationcompiler-dependency-analysis

% f = dir(fullfile(matlabroot, 'toolbox', 'slvnv'));
% for i=1:length(f)
% rmpath(fullfile(matlabroot, 'toolbox', 'slvnv', f(i).name));
% end
%

mcc -o computeFeaturesCompiled -W main:computeFeaturesCompiled ...
  -T link:exe -d /home/mayank/work/FlySpaceTime/code/ ...
  -v /home/mayank/work/FlySpaceTime/code/computeFeaturesCompiled.m ...
  -a /home/mayank/work/FlySpaceTime/code/deepmatching/deepmex.mexa64 ...
  -a /home/mayank/work/FlySpaceTime/code/filehandling ...
  -a /home/mayank/work/FlySpaceTime/code/misc ...
  -a /home/mayank/work/FlySpaceTime/code/optflow_deqing ...
  -a /home/mayank/work/FlySpaceTime/code/toolbox

delete('readme.txt')
delete('mccExcludedFiles.log')
delete('requiredMCRProducts.txt')
delete('PackagingLog.txt')

author: kristinbranson

repo: JAABA


text:
function [yPred, PInvPred]=infoFilterDiscPred(yPrev,PInvPrev,F,Q,u)
%%INFOFILTERDISCPRED Perform the discrete-time prediction step that comes
%                    with the standard linear information filter with
%                    additive process noise.
%
%INPUTS: yPrev The xDimX1 information state at the previous time-step. The
%              information state is the inverse covariance matrix times the
%              target state.
%     PInvPrev The xDimXxDim inverse of the state covariance matrix at the
%              previous time-step.
%            F An xDim X xDim state transition matrix.
%            Q The xDimX xDim process noise covariance matrix. This can be
%              singular.
%            u An optional xDimX1 vector that is the control input. If
%              omitted, no control input is used.
%
%OUTPUTS: yPred The xDim X 1 predicted information state vector.
%      PInvPred The predicted xDim X xDim inverse state covariance matrix.
%
%The implementation of the prediction step given here is from the flow
%chart given in [1]. The matrix G in the paper is omitted, since any
%control input can be pre-multipled by the matrix. The implied discrete-
%time dynamic model is
%x(k)=F(k-1)*x(k-1)+u(k-1)+noise
%where x(k) is the state at time k. However, an information filter
%propagates
%y=inv(P)*x
%and
%PInv=inv(P)
%instead of propagating P and x, where P is the covariance matrix
%associated with the Gaussian state x. This allows the filter to be used
%even when PInv is singular.
%
%More information on information filtering is given in Chapter 7.2 of [2].
%
%REFERENCES:
%[1] D. F. Crouse, P. Willett, and Y. Bar-Shalom, "A low-complexity
%    sliding-window Kalman FIR smoother for discrete-time models," IEEE
%    Signal Processing Letters, vol. 17, no. 2, pp. 177-180, Feb. 2009.
%[2] Y. Bar-Shalom, X. R. Li, and T. Kirubarajan, Estimation with
%    Applications to Tracking and Navigation. New York: John Wiley and
%    Sons, Inc, 2001.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    DInv=F'+PInvPrev/(F)*Q;
    PInvPred=DInv\PInvPrev/(F);

    %Handle possible loss of symmetry due to order of operations and finite
    %precision limitations.
    PInvPred=(PInvPred+PInvPred')/2;

    if(nargin<5||isempty(u))
        yPred=DInv\yPrev;
    else
        yPred=DInv\yPrev+PInvPred*u;
    end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function RR = SqueezeAngularFT(R);
% Usage:
%   RR = SqueezeAngularFT(R);
% Inputs:
%   R    Squared array of Fourier samples. Samples inside an
%        interior square are a priori zero.
%        a-priori zero.
% Outputs:
%   RR   Same as R but with the entries corresponding to the
%        interior square deleted
% Description:
%    R is an area of coefficients obtained after scale and angular
%    separation. Values inside an interior square are identically
%    zero. SqueezeAngularFT essentially removes this interior
%    square.

  nn = size(R);
  n2 = 2*nn(3);
  n = 2*n2;
  boxcnt = nn(2);
  boxlen = nn(4)/2;

  [ix,w] = DetailMeyerWindow([n2/4 n2/2],3);
  alpha_max = max(ix)./n2;
  Lmax = ceil(alpha_max*boxlen);
  mid = (boxlen - Lmax + 1):(boxlen + Lmax);

  RR = zeros(nn(1),nn(2),nn(3),2*Lmax);

  RR = R(:,:,:,mid);

author: thomaskuestner


text:
function BoschettiEnc(img, name, bos_rateE, bos_rateRGB, nBit, tmo_operator)
%
%
%       BoschettiEnc(img, name, bos_rateE, bos_rateRGB, nBit)
%
%
%       Input:
%           -img: input HDR image
%           -name: is output name of the image. If img is empty
%                  an HDR image with filename 'name' is loaded
%           -bos_rateE: JPEG2000 compression rate for the E layer
%           -bos_rateRGB: JPEG2000 compression rate for the RGB layer
%           -nBit: number of bit of the encoding. The maximum is 16.
%           -tmo_operator: an handle to a tone mapping operator
%
%
%     Copyright (C) 2012  Francesco Banterle
%
%     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/>.
%

if(~exist('name', 'var'))
    name = 'bosc_enc';
end

if(~exist('bos_rateRGB', 'var'))
    rateRGB = 15;
else
    rateRGB = bos_rateRGB;
end

if(~exist('bos_rateE', 'var'))
    rateE = 15;
else
    rateE = bos_rateE;
end

if(~exist('nBit','var'))
    nBit = 16;
end

if(~exist('tmo_operator','var'))
    tmo_operator = @ReinhardTMO;
end

%Tone mapping
imgTMO = tmo_operator(img);
%LDR Encoding
imgTMO = GammaTMO(imgTMO, 2.2, 0.0, 0);

%Quantization
maxVal = 2^nBit - 1;
imgTMO = round(imgTMO * maxVal) / maxVal;

%Computing E
epsilon = 1e-4;
epi = 1.0 / maxVal;
E = log2(img ./ (imgTMO + epi) + epsilon);
E = mean(E, 3);

%Encoding E
maxE = max(E(:));
minE = min(E(:));
Eq = (E - minE) / (maxE - minE);

%metadata string
metatadata = [num2str(nBit), ' ', num2str(maxE), ' ', num2str(minE)];

if(nBit == 16)
    Eq = uint16(Eq * maxVal);
end

imwrite(Eq,[name,'_bos_E.jp2'], 'Mode', 'lossy', 'CompressionRatio', rateE);

%Decoding E
EqDec = double(imread([name, '_bos_E.jp2'])) / maxVal;
EDec = EqDec * (maxE - minE) + minE;

%Computing RGB
RGB = zeros(size(img));
div = 2.^EDec;
for i=1:size(img, 3)
    RGB(:,:,i) = (img(:,:,i) ./ div);
end

%Encoding RGB
if(nBit == 16)
    RGB = uint16(RGB * maxVal);
end

nameOut = [name,'_bos_RGB.jp2'];
imwrite(RGB, nameOut, 'Mode', 'lossy', 'CompressionRatio', rateRGB, 'Comment', metatadata);
end

author: banterle


text:
% Matlab Script: cluster_win.m
% ----------------------------
%
% Plots the distance window length and the time window length against magnitude
%
%

%% Magntiude areas
vMagnitude = (0:0.1:10);
vMagnitudea = (0:0.1:6.5);
vMagnitudeb = (6.5:0.1:10);

% Gardener and Knopoff, 1974
vSpaceGaKn74 = 10.^(0.1238*vMagnitude+0.983);
vTimeGaKn74 = 10.^(0.5409*vMagnitudea-0.547);
vTimeGaKn74b = 10.^(0.032*vMagnitudeb+2.7389); % M>=6.5

% Uhrhammer 1976
vSpaceUr = exp(-1.024+0.804*vMagnitude);
vTimeUr = exp(-2.87+1.235*vMagnitude);

% Gruenthal, 1985
vSpaceGr85 = 10.^(0.1060*vMagnitude+1.0982);
vTimeGr85 = 10.^(0.5055*vMagnitude-0.1329);
% Gruenthal, pers. communication
vSpaceGr = exp(1.77+sqrt(0.037+1.02*vMagnitude));
vTimeGra = exp(-3.95+sqrt(0.62+17.32*vMagnitudea));
vTimeGrb = 10.^(2.8+0.024*vMagnitudeb); % M >= 6.5

% Youngs, 1987
%% Maximum window limits
%% Space
vMagY1 = (0:0.01:2.43);
vMagY2 = (2.43:0.01:5.87);
vMagY3 = (5.87:0.01:10);
vSpaceY1 = 20+0*vMagY1;
vSpaceY2 = 10.^(0.1159*vMagY2+1.0197);
vSpaceY3 = 10.^(0.5281*vMagY3-1.3937);
%% Time
vMagTY1 = (0:0.01:3.89);
vMagTY2 = (3.89:0.01:10);
vTimeY1 = 30+0*vMagTY1;
vTimeY2 = 10.^(0.4916*vMagTY2-0.4317);

%% Minimum window limits
%% Space
vMinY1 = (0:0.01:4.41);
vMinY2 = (4.41:0.01:4.98);
vMinY3 = (4.98:0.01:6.42);
vMinY4 = (6.42:0.01:10);
vSpaceMinY1 = 10+0*vMinY1;
vSpaceMinY2 = 10.^(0.5281*vMinY2-1.3290);
vSpaceMinY3 = 10.^(0.3313*vMinY3-0.3490);
vSpaceMinY4 = 10.^(0.1154*vMinY4+1.0371);
%% Time
vMinTY1 = (0:0.01:5);
vMinTY2 = (5:0.01:10);
vTimeMinY1 = 5+0*vMinTY1;
vTimeMinY2 = 10.^(1.0526*vMinTY2-4.561);

% Plotting time window
hPlot = figure_w_normalized_uicontrolunits('tag','Timewindow','Numbertitle','off','Name','Time window length','Units','normalized');
hAxe = axes('tag','ax_timewindow','Nextplot','replace','box','on');

set(gca,'tag','ax_timewindow','Nextplot','replace','box','on','Xticklabel', [0 10 100]);
axs1=findobj('tag','ax_timewindow');
axes(axs1(1));
semilogy(vMagnitudea,vTimeGaKn74,'Color',[1 0 0],'Linewidth', 2);
hold on;
semilogy(vMagnitude,vTimeGr85,'Color',[0 0.5 0],'Linewidth', 2);
semilogy(vMagnitudea,vTimeGra,'Color',[0 0.8 0],'Linewidth', 2);
semilogy(vMagnitude,vTimeUr,'Color',[0.5 0 0],'Linewidth', 2);
semilogy(vMagTY1,vTimeY1,vMagTY2,vTimeY2,'Color',[0 0 1],'Linewidth', 2);
semilogy(vMagnitudeb,vTimeGaKn74b,'Color',[1 0 0],'Linewidth', 2);
semilogy(vMinTY1,vTimeMinY1,vMinTY2,vTimeMinY2,'Color',[0 0 1],'Linewidth', 2);
semilogy(vMagnitudeb,vTimeGrb,'Color',[0 0.8 0],'Linewidth', 2);
legend('Gardner & Knopoff (1974)','Gruenthal (1985)','Gruenthal (pers.)','Urhammer (1976)','Mod. Young (1987)');
xlabel('Magnitude');
ylabel('Time / [days]');
set(gca,'Xlim',[0 9]);%,'Ylim',[0.1 3000],'Yticklabel', [0.1 1 10 100 1000]);
grid on;

% Plotting space window
hPlot2 = figure_w_normalized_uicontrolunits('tag','Spacewindow','Numbertitle','off','Name','Spacial window length','Units','normalized');
hAxe2 = axes('tag','ax_spacewindow','Nextplot','replace','box','on');

set(gca,'tag','ax_spacewindow','Nextplot','replace','box','on');
axs1=findobj('tag','ax_spacewindow');
axes(axs1(1));
semilogy(vMagnitude,vSpaceGaKn74,'Color',[1 0 0],'Linewidth', 2);
hold on;
semilogy(vMagnitude,vSpaceGr85,'Color',[0 0.5 0],'Linewidth', 2);
semilogy(vMagnitude,vSpaceGr,'Color',[0 0.8 0],'Linewidth', 2);
semilogy(vMagnitude,vSpaceUr,'Color',[0.5 0 0],'Linewidth', 2);
semilogy(vMagY1,vSpaceY1,vMagY2,vSpaceY2,vMagY3,vSpaceY3,'Color',[0 0 1],'Linewidth', 2);
semilogy(vMinY1,vSpaceMinY1,vMinY2,vSpaceMinY2,vMinY3,vSpaceMinY3,vMinY4,vSpaceMinY4,'Color',[0 0 1],'Linewidth', 2);
legend('Gardner & Knopoff (1974)','Gruenthal (1985)','Gruenthal (pers. comm)','Urhammer (1976)','Mod. Young (1987)');
ylabel('Distance / [km]')
xlabel('Magnitude');
set(gca,'Xlim',[0 9],'Ylim', [5 300],'Yticklabel', [10 100]);
grid on;

author: CelsoReyes


text:
function cor_n = coriolis(lat, vel, h)
% coriolis: calculates Coriolis forces in the navigation frame.
%
% INPUT
%       lat, Mx1 latitude (radians).
%       vel, Mx3 NED velocities (m/s).
%         h, Mx1 altitude (m).
%
% OUTPUT
%		cor_n, Mx3 Coriolis forces vector in the nav-frame (m/s^2).
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% Reference:
%
%	R. Gonzalez, J. Giribet, and H. Patiño. NaveGo: a
% simulation framework for low-cost integrated navigation systems,
% Journal of Control Engineering and Applied Informatics}, vol. 17,
% issue 2, pp. 110-120, 2015. Eq. 11.
%
% Version: 002
% Date:    2022/01/25
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

M = max(size(lat));

cor_n = zeros(M, 3);

for i = 1:M

    omega_en_n = transport_rate(lat(i), vel(i,1), vel(i,2), h(i));
    omega_ie_n = earth_rate(lat(i));

    cor_n(i,:) = (omega_en_n + 2*omega_ie_n) * vel(i,:)';
end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
%   FactorMarginalization Sums given variables out of a factor in log space.
%   B = FactorMarginalization(A,V) computes the factor with the variables
%   in V summed out. The factor data structure has the following fields:
%       .var    Vector of variables in the factor, e.g. [1 2 3]
%       .card   Vector of cardinalities corresponding to .var, e.g. [2 2 2]
%       .val    Value table of size prod(.card)
%
%   The resultant factor should have at least one variable remaining or this
%   function will throw an error.
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

function B = FactorMarginalization(A,V)
B.var = A.var(2);
B.card = A.card(1);
Val = reshape(A.val,B.card,B.card);

if(V==A.var(2))
    Val = Val';
    B.var = A.var(1);
end

B.val = log(sum(exp(bsxfun(@minus, Val, max(Val)))))+max(Val);

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function variance = r8row_variance ( m, n, x )

%*****************************************************************************80
%
%% R8ROW_VARIANCE returns the variances of an R8ROW.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns in the array.
%
%    Input, real X(M,N), the R8ROW whose row means are desired.
%
%    Output, real VARIANCE(M), the variances of the rows of X.
%
  for i = 1 : m

    mean = sum ( x(i,1:n) ) / n;

    variance(i) = sum ( ( x(i,1:n) - mean ).^2 );

    if ( 1 < n )
      variance(i) = variance(i) / ( n - 1 );
    else
      variance(i) = 0.0;
    end

  end

  return
end

author: johannesgerer


text:
% eegthresh() -  reject trials with out-of-bounds channel values within a
%                specified epoch time range.
%
% Usage:
%   >> [Iin, Iout, newsignal, elec] = eegthresh( signal, frames, ...
%                elecs, negthresh, posthresh, timerange, starttime,endtime);
%
% Required inputs:
%   signal     - 2-D data matrix [channels, frames*sweeps],
%                or 3-D data matrix [channels, frames, sweeps]
%   frames     - number of points per epoch
%   elecs      - [int vector] electrode indices to reject on
%   negthresh  - minimum rejection threshold(s) in uV. This can be an array
%                of values for individual electrodes. If fewer values than
%                electrodes, the last value is used for the remaining
%                electrodes.
%   posthresh  - maximum rejection threshold(s) in uV (same syntax as for
%                negthresh)
%   timerange  - [mintime maxtime] time range limits of the signal
%   starttime  - Starting limit (in seconds or Hz) of range to perform
%                rejection (same syntax  as negthresh)
%   endtime    - Ending limit (in seconds or Hz) of range to perform
%                rejection (same syntax  as negthresh)
%
% Outputs:
%   Iin        - Indexes of epochs accepted
%   Iout       - Indexes of epochs rejected
%   newsignal  - input data after epoch rejection
%   elec       - electrode that triggered the rejections (array of 0s
%                and 1s with the same number of columns as Iout
%                and number of rows = number of electrodes).
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: pop_eegthresh(), eeglab()

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

function [Iin, Iout, newsignal, elec] = eegthresh( signal, pnts, electrodes, negthresh, posthresh, timerange, starttime, endtime);

if nargin < 7
	help eegthresh;
	return;
end;

if starttime < timerange(1)
	disp('eegthresh: starting point out of range, adjusted');
    starttime = timerange(1);
end;
if endtime > timerange(2)
	disp('eegthresh: ending point out of range, adjusted');
    endtime = timerange(2);
end;

% complete thresholds values if necessary
%----------------------------------------
if size(posthresh,2) < size(electrodes,2)
	posthresh = [ posthresh posthresh(end)*ones(1,size(electrodes,2)-size(posthresh,2))];
end;
if size(negthresh,2) < size(electrodes,2)
	negthresh = [ negthresh negthresh(end)*ones(1,size(electrodes,2)-size(negthresh,2))];
end;

% complete timeselect values if necessary
%----------------------------------------
if size(starttime,2) < size(electrodes,2)
	starttime = [ starttime starttime(end)*ones(1,size(electrodes,2)-size(starttime,2))];
end;
if size(endtime,2) < size(electrodes,2)
	endtime = [ endtime endtime(end)*ones(1,size(electrodes,2)-size(endtime,2))];
end;

% find the maximum for each trial
%--------------------------------
sweeps = size(signal(1,:),2)/pnts;
signal = reshape(signal(electrodes,:), size(electrodes(:),1), pnts, sweeps);

% reject the selected trials
%---------------------------
elec = zeros(size(electrodes(:),1), sweeps);
allelec = zeros(1, sweeps);
for indexe = 1:size(electrodes(:),1)
	% transform the time range
	% ------------------------
	framelowlimit  = max(1,floor((starttime(indexe)-timerange(1))/(timerange(2)-timerange(1))*(pnts-1))+1);
	framehighlimit = floor((endtime(indexe)  -timerange(1))/(timerange(2)-timerange(1))*(pnts-1))+1;

	% remove outliers
	% ---------------
	sigtmp = squeeze(signal(indexe,framelowlimit:framehighlimit,:));
    if size(signal,3) == 1, sigtmp = sigtmp'; end;
	sigmax = max(sigtmp, [], 1);
	sigmin = min(sigtmp, [], 1);
	elec(indexe,:) = ( sigmin < negthresh(indexe) ) | ( sigmax > posthresh(indexe) );
	allelec = allelec | elec(indexe,:);
end;
Iout = find( allelec == 1 );
Iin  = find( allelec == 0 );
elec = elec(:,Iout);

% reject the selected trials
%---------------------------
newsignal = reshape(signal, size(signal,1), pnts, sweeps);
if ~isempty(Iin)
	newsignal = newsignal(:,:,Iin);
	if ndims(signal) == 2
		newsignal = newsignal(:,:);
	end;
else
	newsignal = [];
end;
return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [f]=comp_isepdgtreal(coef,g,L,a,M,phasetype)
%COMP_ISEPDGTREAL  Separable IDGT.
%   Usage:  f=comp_isepdgtreal(c,g,L,a,M);
%
%   This is a computational routine. Do not call it directly.
%
%   Input must be in the M x N x W format, so the N and W dimension is
%   combined.
%
%   See also: idgt

%   AUTHOR : Peter L. Søndergaard.
%   TESTING: OK
%   REFERENCE: OK

Lwindow=size(g,1);

if phasetype==1
    % Change from time-invariant phase convention to a
    % frequency-invariant one
    b=L/M;
    M2=floor(M/2)+1;
    N=size(coef,2);
    %M2short=ceil(M/2);

    TimeInd = (0:(N-1))/N;
    FreqInd = (0:(M2-1))*b;

    phase = FreqInd'*TimeInd;
    phase = exp(-2*1i*pi*phase);

    % Handle multisignals
    coef = bsxfun(@times,coef,phase);
end;

if L==Lwindow
    % Do full-window algorithm.

    % Call the computational subroutine.
    f = comp_idgtreal_long(coef,g,L,a,M);

else
    % Do filter bank algorithm.
    % Call the computational subroutine.
    f = comp_idgtreal_fb(coef,g,L,a,M);
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function v = isnaturalnumber(x)

v = isnumeric(x) & x > 0 & round(x) == x;


text:
function Hrefine2D(refineflag)

% function Hrefine2D(refineflag)
% purpose:  apply non-conforming refinement to elements labelled in refineflag

Globals2D;

% 1.1 Count vertices
Nv    = length(VX(:));

% 1.2 Find and count elements to be refined
ref = sort(find(refineflag));
Nrefine = length(ref);

% 1.3 Extract vertex numbers of elements to refine
v1 = EToV(ref, 1); v2 = EToV(ref, 2); v3 = EToV(ref, 3);

% 1.4 Uniquely number all face centers
v4 = max( 1 + Nfaces*(0:K-1)', EToF(:,1) + Nfaces*(EToE(:,1)-1) );
v5 = max( 2 + Nfaces*(0:K-1)', EToF(:,2) + Nfaces*(EToE(:,2)-1) );
v6 = max( 3 + Nfaces*(0:K-1)', EToF(:,3) + Nfaces*(EToE(:,3)-1) );

% 2.0 Extract face center vertices for elements to refine
v4 = v4(ref);      v5 = v5(ref);      v6 = v6(ref);

% 2.1 Renumber face centers contiguously from Nv+1
ids = unique([v4;v5;v6]);
newids(ids) = (1:length(ids))';
v4 = Nv+newids(v4)'; v5 = Nv+newids(v5)'; v6 = Nv+newids(v6)';

% 2.2 Replace original triangle with triangle connecting edge centers
EToV(ref,:) = [v4,v5,v6];

% 3.0 Add extra triangles to EToV
EToV(K+1:K+3*Nrefine,1) = [v1;v2;v3]; % first  vertices of new elements
EToV(K+1:K+3*Nrefine,2) = [v4;v5;v6]; % second vertices of new elements
EToV(K+1:K+3*Nrefine,3) = [v6;v4;v5]; % third  vertices of new elements

% 3.1 Create boundary condition type for refined elements
bcsave = BCType(ref,:);
BCType(ref, :) = 0; % now internal faces

BCType(K+1:K+Nrefine, 1) = bcsave(:, 1);
BCType(K+1:K+Nrefine, 3) = bcsave(:, 3);

BCType(K+Nrefine+1:K+2*Nrefine, 1) = bcsave(:, 2);
BCType(K+Nrefine+1:K+2*Nrefine, 3) = bcsave(:, 1);

BCType(K+2*Nrefine+1:K+3*Nrefine, 1) = bcsave(:, 3);
BCType(K+2*Nrefine+1:K+3*Nrefine, 3) = bcsave(:, 2);

% 3.2 Find vertex locations of elements to be refined
x1 = VX(v1);  x2 = VX(v2);  x3 = VX(v3);
y1 = VY(v1);  y2 = VY(v2);  y3 = VY(v3);

% 3.3 Add coordinates for refined edge centers
VX(v4) = 0.5*(x1+x2); VX(v5) = 0.5*(x2+x3); VX(v6) = 0.5*(x3+x1);
VY(v4) = 0.5*(y1+y2); VY(v5) = 0.5*(y2+y3); VY(v6) = 0.5*(y3+y1);

% 3.4 Increase element count
K = K+3*Nrefine;
return;

repo: Aero-matlab


text:
function patcht(FF,VV,TF,VT,I,Options)
% This function PATCHT, will show a triangulated mesh like Matlab function
% Patch but then with a texture.
%
% Alec: Essentially this calls "surface" for each triangle of the triangle mesh
% with the appropriate texture coordinates
%
% patcht(FF,VV,TF,VT,I,Options);
%
% inputs,
%   FF : Face list #F x 3 with vertex indices
%   VV : Vertices #V x 3
%   TF : Texture list #F x 3 with texture vertex indices
%   VT : Texture Coordinates s 2 x #VT, range must be [0..1] or real pixel postions
%   I : The texture-image RGB [O x P x 3] or Grayscale [O x P]
%   Options : Structure with options for the textured patch such as
%           EdgeColor, EdgeAlpha see help "Surface Properties :: Functions"
%
%   Options.PSize : Special option, defines the image texturesize for each
%           individual  polygon, a low number gives a more block
%           like texture, defaults to 64;
%
% note:
%   On a normal PC displaying 10,000 faces will take about 6 sec.
%
% Example,
%
%  % Load Data;
%   load testdata;
%  % Show the textured patch
%   figure, patcht(FF,VV,TF,VT,I);
%  % Allow Camera Control (with left, right and center mouse button)
%   mouse3d
%
%  % file containing image
%  image_source = 'woody.png';
%  % load image with alpha mask
%  [im,map,alpha] = imread(image_source);
%  % mesh based on alpha values with approximately 100 boundary vertices
%  [V,F] = png2mesh(image_source,0,100);
%  % treat alpha as white in display
%  im = over(im,alpha,ones(size(im)),1);
%  % dummy display white image the same size as input image
%  imshow(ones([max(V(:,2))-min(V(:,2)) max(V(:,1))-min(V(:,1))]));
%  % display original mesh
%  % reverse y-direction so mesh display is upright (we're displaying the mesh
%  % ontop of the imshow() display to get pixel perfect images)
%  patcht(F,[V(:,1) size(im,1)-V(:,2)],F,V,permute(flipdim(im,1),[2 1 3]));
%  % create some trivial deformation
%  new_V = [V(:,1) 2*V(:,2)];
%  figure;
%  % dummy display white image the same size as input image
%  imshow(ones([max(new_V(:,2))-min(new_V(:,2)) max(new_V(:,1))-min(new_V(:,1))]));
%  % display deformed mesh
%  patcht(F,[new_V(:,1) max(new_V(:,2)+0.5)-new_V(:,2)],F,V,permute(flipdim(im,1),[2 1 3]));
%
% Function is written by D.Kroon University of Twente (July 2010)

% FaceColor is a texture
Options.FaceColor='texturemap';
% noe edges
Options.EdgeColor='none';

% Size of texture image used for every triangle
if(isfield(Options,'PSize'))
    sizep=round(Options.PSize(1));
    Options=rmfield(Options,'PSize');
else
    sizep=64;
end

% Check input sizes
if(size(FF,2)~=size(TF,2))
    error('patcht:inputs','Face list must be equal in size to texture-index list');
end

if((ndims(I)~=2)&&(ndims(I)~=3))
    error('patcht:inputs','No valid Input texture image');
end

% Detect if grayscale or color image
switch(size(I,3))
    case 1
        iscolor=false;
    case 3
        iscolor=true;
    otherwise
        error('patcht:inputs','No valid Input texture image');
end

if(max(VT(:))<2)
    % Remap texture coordinates to image coordinates
    VT2(:,1)=(size(I,1)-1)*(VT(:,1))+1;
    VT2(:,2)=(size(I,2)-1)*(VT(:,2))+1;
else
    VT2=VT;
end

% Calculate the texture interpolation values
[lambda1 lambda2 lambda3 jind]=calculateBarycentricInterpolationValues(sizep);

% Split texture-image in r,g,b to allow fast 1D index
Ir=I(:,:,1); if(iscolor), Ig=I(:,:,2); Ib=I(:,:,3); end

% The Patch used for every triangle (rgb)
Jr=zeros([(sizep+1) (sizep+1) 1],class(I));
if(iscolor)
    Jg=zeros([(sizep+1) (sizep+1) 1],class(I));
    Jb=zeros([(sizep+1) (sizep+1) 1],class(I));
end

% if no z coordinates are given then use 0s
if size(VV,2) < 3
  VV = [VV zeros(size(VV,1),1)];
end

hold on;
% Loop through all triangles of the mesh
for i=1:size(FF,1)
    % Get current triangle vertices and current texture-vertices
    V=VV(FF(i,:),:);
    Vt=VT2(TF(i,:),:);

    % Define the triangle as a surface
    x=[V(1,1) V(2,1); V(3,1) V(3,1)];
    y=[V(1,2) V(2,2); V(3,2) V(3,2)];
    z=[V(1,3) V(2,3); V(3,3) V(3,3)];

    % Define the texture coordinates of the surface
    tx=[Vt(1,1) Vt(2,1) Vt(3,1) Vt(3,1)];
    ty=[Vt(1,2) Vt(2,2) Vt(3,2) Vt(3,2)] ;
    xy=[tx(1) ty(1); tx(2) ty(2); tx(3) ty(3); tx(3) ty(3)];

    % Calculate texture interpolation coordinates
    pos(:,1)=xy(1,1)*lambda1+xy(2,1)*lambda2+xy(3,1)*lambda3;
    pos(:,2)=xy(1,2)*lambda1+xy(2,2)*lambda2+xy(3,2)*lambda3;
    pos=round(pos); pos=max(pos,1); pos(:,1)=min(pos(:,1),size(I,1)); pos(:,2)=min(pos(:,2),size(I,2));
    posind=(pos(:,1)-1)+(pos(:,2)-1)*size(I,1)+1;

    % Map texture to surface image
    Jr(jind)=Ir(posind);
    J(:,:,1)=Jr;
    if(iscolor)
        Jg(jind)=Ig(posind);
        Jb(jind)=Ib(posind);
        J(:,:,2)=Jg;
        J(:,:,3)=Jb;
    end

    % Show the surface
    surface(x,y,z,J,Options)
end
hold off;

function [lambda1 lambda2 lambda3 jind]=calculateBarycentricInterpolationValues(sizep)
% Define a triangle in the upperpart of an square, because only that
% part is used by the surface function
x1=sizep; y1=sizep; x2=sizep; y2=0; x3=0 ;y3=0;
% Calculate the bary centric coordinates (instead of creating a 2D image
% with the interpolation values, we map them directly to an 1D vector)
detT = (x1-x3)*(y2-y3) - (x2-x3)*(y1-y3);
[x,y]=ndgrid(0:sizep,0:sizep); x=x(:); y=y(:);
lambda1=((y2-y3).*(x-x3)+(x3-x2).*(y-y3))/detT;
lambda2=((y3-y1).*(x-x3)+(x1-x3).*(y-y3))/detT;
lambda3=1-lambda1-lambda2;
% Make from 2D (surface)image indices 1D image indices
[jx jy]=ndgrid(sizep-(0:sizep)+1,sizep-(0:sizep)+1);
jind=(jx(:)-1)+(jy(:)-1)*(sizep+1)+1;

author: alecjacobson

repo: gptoolbox


text:
% Comparing the SPM beta weights and the ones derived from mrVista
%
% We ran SPM version 2 in this directory.  The time series data are stored
% in the imageSlice sub-directory. The data here is a toy set with only one
% slice that allows us to debug the beta weight computations in the mrVista
% GLM and in the SPM analysis.
%
% When we ran the SPM model we created an SPM.mat file that is in this
% directory that contains all of the relevant information for the two types
% of events that take place in the experiment: rhymes and lines.  You can
% view the design matrix by running SPM2 | fMRI-time series | Review (or
% something like that). The design matrix is also just read in below.
%

% We could only run the SPM2 code on slightly older versions of Matlab, say
% Matlab 7.0.4.  This doesn't run properly on Matlab 2007a, for example.

%% The first part is about properly computing the model weights
% This part confirmed that our weights and SPM are the same up to a single
% global scale factor.

% Current directory.
chdir(fullfile(mrvRootPath,'EventRelated','glmTest'))

% This is the original SPM Michal and I created, but I cleared the .xM
% field so that there is no mask.  This makes the data comparable to Rory's
% calculation, point by point.
load('SPM-noMask');

% The spm calculation of the betas is in spm_spm_local
spm_defaults
[SPM,pKX] = spm_spm_local(SPM);

% We set X, the design matrix, to be the one we built using the SPM gui.
% this design matrix has 1 predictor per condition, plus a column of ones.
% it uses the hrf from SPM
X = SPM.xX.X;
% figure(1); plot(X)
% tmp = pinv(pKX); max(abs(X(:) - tmp(:)))

% We created the beta files using the SPM interface (Estimate).  We read
% the files here using the builtin Matlab analyze reader.
% chdir('SPM-UI')
tmp = analyzeRead('beta_0001');  spmBeta(:,1) = tmp(:);
tmp = analyzeRead('beta_0002');  spmBeta(:,2) = tmp(:);
tmp = analyzeRead('beta_0003');  spmBeta(:,3) = tmp(:);
[r,c] = size(tmp);
% imtool(reshape(spmBeta(:,1),r,c))
% figure(1); imagesc(reshape(spmBeta(:,1),r,c)); axis image;

% We convert the time series in analyze format into our own time series.
inFileRoot = ...
    fullfile(mrvRootPath,'EventRelated',...
    'glmTest',...
    'imageFiles',...
    'Scan1-Slice6');
outFileRoot = inFileRoot;
nVols = 78;
firstVolIndex = 0;
doRotate = 0;
scaleFact = [1,1];
flipudFlag = 0;
fliplrFlag = 0;
tSeries = analyze2mrLoadRet3TSeries(inFileRoot,outFileRoot,nVols,...
    firstVolIndex,doRotate,scaleFact,flipudFlag,fliplrFlag);

% figure(1); plot(tSeries)
%   max(abs(KWY(:) - tSeries(:)))

% Compute the beta weights from mrVista
tr = 2; nh = 1;
model    = glm(tSeries, X, tr, nh);
mrvBetas = squeeze(model.betas);

% Compare the beta weights for the different model terms
ii = 1;
mrvB = reshape(mrvBetas(ii,:),r,c);
figure(2); imagesc(mrvB); axis image; colormap(gray(256))
figure(1); hist(mrvB(:),100)

% Compare the values directly.
figure(1); plot(mrvBetas(ii,:),spmBeta(:,ii),'.')
grid on

% We can check the images read through the SPM call and our call based on
% analyze Read. They are the same.
% vName = 'imageFiles\Scan1-Slice6\000.img';
vName = 'C:\u\brian\Matlab\VISTASOFT\mrLoadRet-3.0\EventRelated\glmTest\imageFiles\Scan1-Slice6\077.img';

tmp = analyzeRead(vName);
figure(1); hist(tmp(:),100)
imtool(tmp/max(tmp(:)));
mean(tmp(:))

% V = VY(i)
V.fname = vName;
V.pinfo = [0.0066, 0, 0]';
r = 124; c= 99;
[X,Y] = meshgrid(1:r,1:c);
Z = ones(size(X(:)));
XYZ = [X(:),Y(:),Z(:)]';
spmY = spm_get_data(V,XYZ);
figure(2);hist(spmY(:),100)

spmY = reshape(spmY,c,r)';
imtool(spmY/max(spmY(:)));
mean(spmY(:))

%% We analyze the t-contrast images using mrVista.  We created a
% t-contrast using SPM via its GUI.  These are in SPMT004 and con_004.
chdir(fullfile(mrvRootPath,'EventRelated','glmTest'))

% For the contrasts, we used this SPM
load('SPM');
X = SPM.xX.X;

% Get the time series
inFileRoot = ...
    fullfile(mrvRootPath,'EventRelated',...
    'glmTest',...
    'imageFiles',...
    'Scan1-Slice6');
outFileRoot = inFileRoot;
nVols = 78;
firstVolIndex = 0;
doRotate = 0;
scaleFact = [1,1];
flipudFlag = 0;
fliplrFlag = 0;
tSeries = analyze2mrLoadRet3TSeries(inFileRoot,outFileRoot,nVols,...
    firstVolIndex,doRotate,scaleFact,flipudFlag,fliplrFlag);

% Build the GLM model
tr = 2; nh = 1;
model    = glm(tSeries, X, tr, nh);

active = 1;   % Rhyme condition
control= 2;   % Line condition
[stat, ces, vSig, units] = glm_contrast(model,active,control);
[statT, ces, vSig, units] = glm_contrast(model,active,control,'t');

r = 124; c= 99;
mrvCON = reshape(stat,r,c);
figure(1); imagesc(mrvCON); axis image

spmCON = analyzeRead('con_0004');
spmCON = reshape(spmCON,r,c);
figure(2); imagesc(spmCON); axis image
imtool(spmCON)

plot(spmCON(:),mrvCON(:),'.'); axis equal; grid on

spmT = analyzeRead('SPMT_0004');
spmT = reshape(spmT,r,c);
figure(3); imagesc(spmT); axis image

mrvCONT = reshape(statT,r,c);
figure(1); imagesc(mrvCONT); axis image

% Notice that there are a lot of mrvCONT values that are high were the SPMT
% values are zero.  We think these are from places that SPM masks and mrV
% doesn't.  We need to figure this out.
plot(spmT(:),mrvCONT(:),'.');axis equal; grid on

% hist(mrvCONT(:),20)

l = (spmT(:) == 0);
tmp = zeros(size(spmT));
tmp(l) = 1;
tmp = reshape(tmp,r,c);
figure(4); subplot(1,2,1), imagesc(tmp)
colorbar('vert'); axis image

l = (mrvCONT(:) == 0);
tmp = zeros(size(mrvCONT));
tmp(l) = 1;
tmp = reshape(tmp,r,c);
figure(4); subplot(1,2,2), imagesc(tmp)
colorbar('vert'); axis image

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function test_bug3195

% MEM 8gb
% WALLTIME 00:20:00
% DEPENDENCY ft_prepare_sourcemodel ft_inside_headmodel

%%
% load the data

cd(dccnpath('/home/common/matlab/fieldtrip/data/test/bug3195'));

load('template_FEM.mat')
load('sourcemodel.mat')

% the system matrix has not been computed yet, hence ft_headmodeltype fails
template_FEM.type = 'simbio';

%%
% take a subset of the cortical sheet vertices

sel = randperm(8004);
sel = sel(1:1000);

cfg = [];
cfg.headmodel = template_FEM;
cfg.sourcemodel.pos = sourcemodel.pos(sel,:);
sourcemodel1 = ft_prepare_sourcemodel(cfg);

figure
ft_plot_mesh(sourcemodel1.pos(sourcemodel1.inside,:))

%%
% construct a coarse 3D grid

cfg = [];
cfg.headmodel = template_FEM;
cfg.resolution = 20;
sourcemodel2 = ft_prepare_sourcemodel(cfg);

figure
ft_plot_mesh(sourcemodel2.pos(sourcemodel2.inside,:))


text:
function test_bug1571

% MEM 5gb
% WALLTIME 00:10:00
% DEPENDENCY ft_preprocessing ft_preproc_dftfilter dftfilter ft_freqanalysis ft_singleplotER

% Philipp Hintze wrote:
% The issue with the continuous data is the following: I tried using
% [this script]
% but looking at averaged timelocked data for the subject, the data looks
% identical to the result I get using no filter at all, i.e., contaminated by
% strong line noise. Including padding makes no difference either, therefore
% the brackets.

dataset = dccnpath('/home/common/matlab/fieldtrip/data/test/bug1571/20101230_1010.cnt');

cfg                     = [];
cfg.channel      = {'EEG'};
cfg.datafile     = dataset;
cfg.headerfile   = dataset;
cfg.dataset      = dataset;
cfg.continuous   = 'yes';

unfilteredContinuousData  = ft_preprocessing(cfg);

cfg.dftfilter           = 'yes';
filteredContinuousData  = ft_preprocessing(cfg);

assert(~isequal(filteredContinuousData.trial{1}(:), unfilteredContinuousData.trial{1}(:)))

cfg = [];
cfg.channel = {'Fp1'};
cfg.method = 'mtmfft';
cfg.taper = 'hanning';
cfg.foi = 1:100;
fftFilt = ft_freqanalysis(cfg, filteredContinuousData);
figure;
subplot(2,1,1)
ft_singleplotER(cfg, fftFilt)

cfg = [];
cfg.channel = {'Fp1'};
cfg.method = 'mtmfft';
cfg.taper = 'hanning';
cfg.foi = 1:100;
fftUnfilt = ft_freqanalysis(cfg, unfilteredContinuousData);
subplot(2,1,2)
ft_singleplotER(cfg, fftUnfilt)

%% synthetic signal with wrong sampling rate

fsample = 512;

data = {};
data.fsample = fsample;
data.label = {'1', '2', '3'};
nchans = numel(data.label);
real_fsample = 512.01;
data.time{1} = 0:1/real_fsample:100;
data.trial{1} = rand(nchans, numel(data.time{1}));

cfg                     = [];
cfg.continuous          = 'yes';
cfg.dftfilter           = 'yes';
filteredData  = ft_preprocessing(cfg, data);

assert(~isequal(data.trial{1}, filteredData.trial{1}(:)))

fprintf('done.\n');