text:
function nfgWriteConTrackParams(ctrparamsFile, binDir, wmFile, roiFile, pdfFile, nMinNodes, nMaxNodes, fStepSizeMm, nSamples)
%Write the parameters file for ConTrack programs
%
%   nfgWriteConTrackParams(ctrparamsFile, binDir, wmFile, roiFile, pdfFile,
%   nMinNodes, nMaxNodes, fStepSizeMm, nSamples)
%
%   contrack_gen and contrack_score require this file to run.
%   ctrparamsFile : Filename to write.
%   binDir : Directory to find images.
%   wmFile : White matter mask relative to binDir.
%   roiFile : ROI mask relative to binDir.
%   pdfFile : PDF image relative to binDir.
%
% AUTHORS:
%   2009.08.05 : AJS wrote it
%
% NOTES:

if ieNotDefined('nSamples'); nSamples = 50000; end

ctr = ctrCreate();
ctr=ctrSet(ctr,'image_directory',binDir);
ctr=ctrSet(ctr,'fa_filename',wmFile);
ctr=ctrSet(ctr,'mask_filename',roiFile);
ctr=ctrSet(ctr,'pdf_filename',pdfFile);
ctr=ctrSet(ctr,'max_nodes', nMaxNodes);
ctr=ctrSet(ctr,'min_nodes', nMinNodes);
ctr=ctrSet(ctr,'step_size', fStepSizeMm);
ctr=ctrSet(ctr,'desired_samples', nSamples);
ctrSave(ctr,ctrparamsFile);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [voxAmps, tS] = twEventPrediction(mv, rmParams, anal);
% twEventPrediction - compute a set of predicted voxel amplitudes
% and time series for a  set of voxels, given retinotopy parameters estimated
% from a traveling wave analysis.
%
%  [voxAmps tS] = twEventPrediction(mv, rmParams, anal);
%
% The 'mv' struct is a multivoxel struct, created by mv_init or
% multiVoxelUI.
%
% The anal struct is the output from rmVisualizeRFs. The voxels flag
% specifies the voxels for which to make predictions (defaulting to all
% voxels in the current ROI), and the roi flag specifies which ROI to make
% predictions for in the anal struct (default is the first one).
%
% ras, 09/2006.
if nargin<3, error('Not enough input args.');               end
if notDefined('roi'),    roi = 1;                           end
if notDefined('voxels'), voxels = 1:length(anal.x0{roi});   end

nConds = sum(mv.trials.condNums>0);
nVoxels = length(voxels);

% to properly normalize the voxel amplitudes, we need
% to divide by a measure of the number of samples per square degree:
% the idea is that the density per deg^2 is constant over time, but
% sampling may change.
dX = unique(diff(anal.X(:))); dY = unique(diff(anal.Y(:)));
dX = dX(2); dY = dY(2); % should be at least 2 unique vals, 0 and step size
samplesPerDeg2 = ceil( [1/dX] * [1/dY] );

% convert stimuli in rmParams to a more useable format
xRange = unique(anal.X); yRange = unique(anal.Y);
[rmParams images] = rmStimulusMatrix(rmParams, xRange, yRange, 0);

% compute voxAmps
fprintf('[%s] Computing predicted voxel amplitudes ', mfilename);
for c = 1:nConds
    mask = logical( round(images(:,:,c)) );

    for v = 1:nVoxels
        ii = voxels(v);
        sigma = anal.sigma{roi}(ii) .* 2;
        x0 = anal.x0{roi}(ii);
        y0 = anal.y0{roi}(ii);
        beta = anal.beta{roi}(ii);
%         sig = anal.log10p{roi}(ii);
        RF = beta * rfGaussian2D(anal.X, anal.Y, sigma, sigma, 0, x0, y0);

        % multiply by the two stimulus representations
        voxAmps(v,c) = sum( RF(mask) ); % / [sum(mask(:))/samplesPerDeg2];
    end

    fprintf('.');
end
fprintf(' done.\n');

% using voxAmps as scaling betas, compute predicted tSeries
erParams = mv.params; erParams.glmHRF = -1; % delta function only
X = glm_createDesMtx(mv.trials, erParams);
X = X(:,2:end); % remove null condition onsets
nFrames = size(X, 1);
betas = ones(1, size(X,2)); % note this is beta for the prediction, not the
                            % retinotopy model
for v = 1:nVoxels
    betas(1:nConds) = voxAmps(voxels(v),:);
%     timeSeries = conv2(X * betas', rmParams.analysis.Hrf{1}, 'full');
    timeSeries = filter(rmParams.analysis.Hrf{1}, 1, X * betas');
    tS(:,v) = timeSeries(1:nFrames);
end

% remove NaNs
voxAmps( isnan(voxAmps) ) = 0;
tS( isnan(tS) ) = 0;

return


text:
function cdf = inverse_gaussian_cdf ( x, a, b )

%*****************************************************************************80
%
%% INVERSE_GAUSSIAN_CDF evaluates the Inverse Gaussian CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the CDF.
%    0.0 < X.
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < A,
%    0.0 < B.
%
%    Output, real CDF, the value of the CDF.
%
  if ( x <= 0.0 )

    cdf = 0.0;

  else

    x1 = sqrt ( b / x ) * ( x - a ) / a;
    cdf1 = normal_01_cdf ( x1 );

    x2 = - sqrt ( b / x ) * ( x + a ) / a;
    cdf2 = normal_01_cdf ( x2 );

    cdf = cdf1 + exp ( 2.0 * b / a ) * cdf2;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [theta,cmethod,interpolation,preconditioner] = amgoption(option)
%% AMGOPTIONS default options for algebraic multigrid solver
%
%   theta = 0.025;	cmethod = 'C'; interpolation = 'S'; preconditioner = 'W';
%
%   Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if isfield(option,'theta')
    theta = option.theta;
else
    theta = 0.025;
end

if isfield(option,'cmethod')
    cmethod = option.cmethod;
else
    cmethod = 'C';
end

if isfield(option,'interpolation')
    interpolation = option.interpolation;
else
    interpolation = 'T';
end
if isfield(option,'preconditioner')
    preconditioner = option.preconditioner;
else
    preconditioner = 'W';
end

repo: Aero-matlab


text:
clear();

numRows     = 3;
numVectors  = 6;

%% Load / Generate Data

mW = randn(numRows, numRows);
mW = (mW.' * mW) + (0.00005 + eye(numRows));

mXX = randn(numRows, numVectors); %<! The set of {x}_{i} (Each column)

mX = zeros(numVectors, numRows * numRows);
for ii = 1:numVectors
    mX(ii, :) = kron(mXX(:, ii).', mXX(:, ii).');
end

vY = randn(numVectors, 1);

objVal = 0;

for ii = 1:numVectors
    objVal = objVal + (((mXX(:, ii).' * mW * mXX(:, ii)) - vY(ii)) ^ 2);
end

objVal
sum((mX * mW(:) - vY) .^ 2)

% ii = 8;
% mXX(:, ii).' * mW * mXX(:, ii)
% mX(ii, :) * mW(:)

author: RoyiAvital

repo: StackExchangeCodes


text:
function [varargout] = bst_simmeeg(varargin)
% BST_SIMMEEG:  Calls SimMEEG GUI to simulate signals and imports results in database
%
% USAGE:  NewFiles = bst_simmeeg('GUI', iStudy)

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

eval(macro_method);
end

%% ===== SIMMEEG GUI =====
function NewFiles = GUI(iStudy)
    NewFiles = {};
    % Install/load SimMEEG plugin
    [isOk, errMsg] = bst_plugin('InstallInteractive', 'simmeeg');
    if ~isOk
        return;
    end
    % Progress bar
    bst_progress('start', 'SimMEEG', 'Initializing...');
    bst_plugin('SetProgressLogo', 'simmeeg');
    % Load anatomy + sensors + headmodel
    bst = LoadInputs(iStudy);
    % Reset SimMEEG global variable h
    global h;
    h = [];
    % Call SimMEEG
    try
        SimMEEG_GUI(bst);
        isStarted = 1;
    catch
        bst_error();
        isStarted = 0;
    end
    % Close progress bar
    bst_progress('stop');
    bst_plugin('SetProgressLogo', []);
    % If the main figure was created
    if isfield(h, 'main_fig') && ~isempty(h.main_fig) && ishandle(h.main_fig)
        % Initialization successful: Wait for it to be closed
        if isStarted
            waitfor(h.main_fig);
        % Initialization failed: Close it
        else
            close(h.main_fig);
            return;
        end
    end
    % Import simulated signals
    NewFiles = ImportSimulations(iStudy);
    % Clear memory
    h = [];
end

%% ===== LOAD INPUTS =====
% Set all fields requested in sm_bst2ft_anatomy_from_bst_files
function bst = LoadInputs(iStudy)
    % Get study and subject info
    ProtocolInfo = bst_get('ProtocolInfo');
    sStudy = bst_get('Study', iStudy);
    sSubject = bst_get('Subject', sStudy.BrainStormSubject);
    % Initialize FieldTrip
    [isInstalled, errMsg, PlugFt] = bst_plugin('Install', 'fieldtrip', 1, '20200911');
    if ~isInstalled
        bst_error(errMsg, 'Plugin manager', 0);
        return;
    end
    % Folders
    bst.FieldTrip_dir = bst_fileparts(which(PlugFt.TestFile));  % (only "fieldtrip-20200911" tested)
    bst.subj_anat_dir = bst_fullfile(ProtocolInfo.SUBJECTS, bst_fileparts(sStudy.BrainStormSubject));
    bst.subj_data_dir = bst_fullfile(ProtocolInfo.STUDIES, bst_fileparts(sStudy.FileName));
    % MRI file
    if ~isempty(sSubject.iAnatomy)
        bst.subj_MriFile = file_fullpath(sSubject.Anatomy(sSubject.iAnatomy).FileName);
    else
        bst.subj_MriFile = [];
    end
    % Scalp
    if ~isempty(sSubject.iScalp)
        bst.subj_scalpFile = file_fullpath(sSubject.Surface(sSubject.iScalp).FileName);
    else
        bst.subj_scalpFile = [];
    end
    % Outer skull
    if ~isempty(sSubject.iOuterSkull)
        bst.subj_skullFile = file_fullpath(sSubject.Surface(sSubject.iOuterSkull).FileName);
    else
        bst.subj_skullFile = [];
    end
    % Inner skull
    if ~isempty(sSubject.iInnerSkull)
        bst.subj_brainFile = file_fullpath(sSubject.Surface(sSubject.iInnerSkull).FileName);
    else
        bst.subj_brainFile = [];
    end
    % Cortex => The one used for the computation of the selected forward model
    if ~isempty(sSubject.iCortex)
        bst.subj_cortexFile = file_fullpath(sSubject.Surface(sSubject.iCortex).FileName);
    else
        bst.subj_cortexFile = [];
    end

    bst.subj_wmFile = bst.subj_cortexFile;
    bst.subj_pialFile = bst.subj_cortexFile;

    % Get channel files
    bst.subj_sens_meg_file = [];
    bst.subj_sens_eeg_file = [];
    if ~isempty(sStudy.Channel)
        if ismember('MEG', sStudy.Channel.Modalities)
            bst.subj_sens_meg_file = file_fullpath(sStudy.Channel.FileName);
        end
        if ismember('EEG', sStudy.Channel.Modalities)
            bst.subj_sens_eeg_file = file_fullpath(sStudy.Channel.FileName);
        end
    end

    % Get head model
    bst.subj_hdm_meg_vol_file = [];
    bst.subj_hdm_meg_cortex_file = [];
    bst.subj_hdm_eeg_vol_file = [];
    bst.subj_hdm_eeg_cortex_file = [];
    if ~isempty(sStudy.iHeadModel)
        hm = sStudy.HeadModel(sStudy.iHeadModel);
        % Surface
        if strcmpi(hm.HeadModelType, 'surface')
            % MEG
            if ~isempty(hm.MEGMethod)
                bst.subj_hdm_meg_cortex_file = file_fullpath(hm.FileName);
            end
            % EEG
            if ~isempty(hm.EEGMethod)
                bst.subj_hdm_eeg_cortex_file = file_fullpath(hm.FileName);
            end
        % Volume
        elseif strcmpi(hm.HeadModelType, 'volume')
            % MEG
            if ~isempty(hm.MEGMethod)
                bst.subj_hdm_meg_vol_file = file_fullpath(hm.FileName);
            end
            % EEG
            if ~isempty(hm.EEGMethod)
                bst.subj_hdm_eeg_vol_file = file_fullpath(hm.FileName);
            end
        end
        % If a reference cortex surface is defined in the file: overwrite the default cortex
        hmMat = in_bst_headmodel(hm.FileName, 0, 'SurfaceFile');
        if ~isempty(hmMat.SurfaceFile) && file_exist(file_fullpath(hmMat.SurfaceFile))
            bst.subj_cortexFile = file_fullpath(hmMat.SurfaceFile);
        end
    end
end

%% ===== IMPORT SIMULATIONS =====
% Get the signals simulated by SimMEEG from the global data h, and import them in the Brainstorm DB
function NewFiles = ImportSimulations(iStudy)
    % Global SimMEEG variable
    global h;
    % Returned files
    NewFiles = {};
    % No simulated data
    if isempty(h) || ~isfield(h, 'sim_data') || isempty(h.sim_data)
        return;
    end
    % Get available fields
    allFields = {'sens_final', 'sig_final', 'sig_wav', 'prepost_wav', 'noise_wav', 'prepost_win', 'sig_win', ...
        'sens_noise', 'sens_noise_scaled', 'sens_sig_data', 'sens_noise_final', 'sens_final_org', 'sens_noise_final_org', 'sens_sig_data_org'};
    availFields = {};
    for i = 1:length(allFields)
        if isfield(h.sim_data, allFields{i}) && ~isempty(h.sim_data.(allFields{i}))
            availFields{end+1} = allFields{i};
        end
    end
    if isempty(availFields)
        return;
    end
    % Ask fields to import
    isSelect = strcmpi(availFields, 'sens_final');
    if ~any(isSelect)
        isSelect = strcmpi(availFields, 'sig_final');
    end
    isSelect = java_dialog('checkbox', 'Select variables to import:', 'Import SimMEEG output', [], availFields, isSelect);
    if isempty(isSelect)
        return;
    end
    importGroups = availFields(isSelect == 1);
    % Get study
    sStudy = bst_get('Study', iStudy);
    % Timestamp
    c = clock;
    strTime = sprintf('_%02.0f%02.0f%02.0f_%02.0f%02.0f', c(1)-2000, c(2:5));
    % Progress bar
    nTrials = length(importGroups) * size(h.sim_data.(importGroups{1}),3);
    bst_progress('start', 'SimMEEG', 'Importing simulated files...', 0, nTrials);
    % Detect previous executions of simmeeg in this study
    iRun = 1;
    if ~isempty(sStudy.Matrix)
        while any(cellfun(@(c)and(length(c)>3, strcmpi(c(1:3),sprintf('%02d-',iRun))), {sStudy.Matrix.Comment}))
            iRun = iRun + 1;
        end
    end
    % Load channel file
    if ~isempty(sStudy.Channel)
        ChannelMat = in_bst_channel(sStudy.Channel.FileName);
        nChannels = length(ChannelMat.Channel);
        iChannels = cellfun(@(c)find(strcmpi(c,{ChannelMat.Channel.Name}),1), h.anatomy.sens.label);
    end
    % Loop to import groups of trials
    for iGroup = 1:length(importGroups)
        % List of trials to import
        trials = h.sim_data.(importGroups{iGroup});
        nSources = size(trials,2);
        % Save matrix
        if ismember(importGroups{iGroup}, {'sig_final', 'sig_wav', 'prepost_wav', 'noise_wav', 'prepost_win', 'sig_win'})
            % Create a "matrix" structure
            sMat = db_template('matrixmat');
            sMat.Time        = h.sim_data.cfg.study.lat_sim;
            sMat.Description = cell(nSources,1);
            for iSource = 1:nSources
                sMat.Description{iSource} = sprintf('Source %d', iSource);
            end
            % Add history entry
            sMat = bst_history('add', sMat, 'process', 'Generated with SimMEEG');
            % Add extra cfg structure
            sMat.cfg = h.sim_data.cfg;
            % Loop on each trial
            for iTrial = 1:size(trials,3)
                % Trial values
                sMat.Comment = sprintf('%02d-%s (#%d)', iRun, importGroups{iGroup}, iTrial);
                sMat.Value   = trials(:,:,iTrial)';
                % Output filename
                FileName = bst_process('GetNewFilename', bst_fileparts(sStudy.FileName), ['matrix_simmeeg', strTime, '_', strrep(importGroups{iGroup}, '_', '-'), sprintf('_trial%03d', iTrial)], 0);
                % Save on disk
                bst_save(FileName, sMat, 'v6');
                NewFiles{end+1} = FileName;
                % Add structure to database
                sNew = db_template('Matrix');
                sNew.FileName = file_short(FileName);
                sNew.Comment  = sMat.Comment;
                iItem = length(sStudy.Matrix) + 1;
                sStudy.Matrix(iItem) = sNew;
                % Increment progress bar
                bst_progress('inc',1);
            end
        % Save MEG/EEG data
        else
            % Create a "matrix" structure
            sMat = db_template('datamat');
            sMat.Time = h.sim_data.cfg.study.lat_sim;
            sMat.ChannelFlag = ones(nChannels,1);
            sMat.Device = 'SimMEEG';
            % Add history entry
            sMat = bst_history('add', sMat, 'process', 'Generated with SimMEEG');
            % Add extra cfg structure
            sMat.cfg = h.sim_data.cfg;
            % Loop on each trial
            for iTrial = 1:size(trials,3)
                % Trial values
                sMat.Comment = sprintf('%02d-%s (#%d)', iRun, importGroups{iGroup}, iTrial);
                sMat.F = zeros(nChannels, size(trials,1));
                sMat.F(iChannels,:) = trials(:,:,iTrial)';
                % Output filename
                FileName = bst_process('GetNewFilename', bst_fileparts(sStudy.FileName), ['data_simmeeg', strTime, '_', strrep(importGroups{iGroup}, '_', '-'), sprintf('_trial%03d', iTrial)], 0);
                % Save on disk
                bst_save(FileName, sMat, 'v6');
                NewFiles{end+1} = FileName;
                % Add structure to database
                sNew = db_template('Data');
                sNew.FileName = file_short(FileName);
                sNew.Comment  = sMat.Comment;
                iItem = length(sStudy.Data) + 1;
                sStudy.Data(iItem) = sNew;
                % Increment progress bar
                bst_progress('inc',1);
            end
        end
    end
    % Update database
    bst_set('Study', iStudy, sStudy);
    panel_protocols('UpdateNode', 'Study', iStudy);
    % Close progress bar
    bst_progress('stop');
end

author: brainstorm-tools


text:
function pde = Helmholtzdata3
%% HELMHOTLZDATA4 trigonometric  data for PML Helmholtz equation
% -laplace u -k^2 u = f;
% The Intresting domain is [0,1]^2; The computing domain including the PML
%domain is [-0.1,1.1]^2;
% boundary condition
% u =0 on the PML boundary;
% Created by Jie Zhou.
% The data come from
% Advances in Iterative Methods and Preconditioners for the Helmholtz Equation
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

pde = struct('f',@f,'k2',@k2,'g_D',@g_D,'d',@d);

    % load data (right hand side function)
    function rhs =  f(p)

      x = p(:,1); y = p(:,2);
     k0 = size(p,1);
    rhs = zeros(k0,1);
     ii = find(abs(x-0.5)<=1.0e-2&abs(y-0.5)<=1.0e-2);
    rhs(ii) = 10000;
    end

    function PMLp = d(p)
    global omegal k  c
    c = 340;
    omegal = k*340;
          i = sqrt(-1);
          x = p(:,1); y = p(:,2);
      size0 = size(x,1);
     sigma1 = ones(size0,1);
     sigma2 = ones(size0,1);
    index = find((x<0)|(x>1)|(y>1)|(y<0));
    sigma1(index) = x(index);
    sigma2(index) = y(index);

    s1 = 1 + sigma1/((i)*omegal);
    s2 = 1 + sigma2/((i)*omegal);

    PMLp = [s2./s1 s1./s2];
    end
    function s =  g_D(p)
    s = zeros(size(p,1),1);
    end
    function wavenumber = k2(p)
    global k omegal  c
    omegal = k*340;
          i = sqrt(-1);
          x = p(:,1); y = p(:,2);
      size0 = size(x,1);
     sigma1 = ones(size0,1);
     sigma2 = ones(size0,1);
    index = find((x<0)|(x>1)|(y>1)|(y<0));
    sigma1(index) = x(index);
    sigma2(index) = y(index);

    s1 = 1 + sigma1/((i)*omegal);
    s2 = 1 + sigma2/((i)*omegal);
    wavenumber = k.^2.*s1.*s2;
    end

end

repo: ifem


text:
function r8mat_norm_fro_test ( )

%*****************************************************************************80
%
%% R8MAT_NORM_FRO_TEST tests R8MAT_NORM_FRO.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2014
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = 4;
  a = zeros ( m, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8MAT_NORM_FRO_TEST\n' );
  fprintf ( 1, '  R8MAT_NORM_FRO computes a Frobenius norm of an R8MAT;\n' );

  t1 = 0.0;
  k = 0;
  for i = 1 : m
    for j = 1 : n
      k = k + 1;
      a(i,j) = k;
      t1 = t1 + k * k;
    end
  end

  t1 = sqrt ( t1 );

  r8mat_print ( m, n, a, '  A:' );

  t2 = r8mat_norm_fro ( m, n, a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Expected norm = %g\n', t1 );
  fprintf ( 1, '  Computed norm = %g\n', t2 );

  return
end


text:
function [c, r, t, cpulse] = tapas_physio_read_physlogfiles_GE(log_files, ...
verbose)
% reads out physiological time series and timing vector depending on the
% MR scanner vendor and the modality of peripheral cardiac monitoring (ECG
% or pulse oximetry)
%
%   [cpulse, rpulse, t, c] = tapas_physio_read_physlogfiles_GE(logfiles, ...
%                               verbose)
%
%   NOTE: if one
%
% IN    log_files
%       .log_cardiac        contains ECG or pulse oximeter time course
%                           for GE: ECGData...
%       .log_respiration    contains breathing belt amplitude time course
%                           for GE: RespData...
%       .sampling_interval  1 entry: sampling interval (seconds)
%                           for both log files
%                           2 entries: 1st entry sampling interval (seconds)
%                           for cardiac logfile, 2nd entry for respiratory
%                           logfile
%       verbose
%       .level              debugging plots are created if level >=3
%       .fig_handles        appended by handle to output figure
%
% OUT
%   cpulse              time events of R-wave peak in cardiac time series (seconds)
%                       for GE: usually empty
%   r                   respiratory time series
%   t                   vector of time points (in seconds)
%                       NOTE: This assumes the default sampling rate of 40
%                       Hz
%   c                   cardiac time series (ECG or pulse oximetry)
%
% EXAMPLE
%   [ons_secs.cpulse, ons_secs.rpulse, ons_secs.t, ons_secs.c] =
%       tapas_physio_read_physlogfiles_GE(logfiles);
%
%   See also tapas_physio_main_create_regressors

% Author: Lars Kasper
% Created: 2013-02-16
% Copyright (C) 2013 Institute for Biomedical Engineering, ETH/Uni Zurich.
%
% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

%% read out values
DEBUG = verbose.level >= 3;

hasRespirationFile = ~isempty(log_files.respiration);
hasCardiacFile = ~isempty(log_files.cardiac);

if hasRespirationFile
    r = load(log_files.respiration, 'ascii');
else
    r = [];
end

if hasCardiacFile
    c = load(log_files.cardiac, 'ascii');
else
    c = [];
end

%% interpolate to greater precision, if 2 different sampling rates are given

dt = log_files.sampling_interval;

hasDifferentSamplingRates = numel(dt) > 1;

if hasDifferentSamplingRates && hasCardiacFile && hasRespirationFile
    dtCardiac = dt(1);
    dtRespiration = dt(2);
    isHigherSamplingCardiac = dtCardiac < dtRespiration;

    nSamplesRespiration = size(r,1);
    nSamplesCardiac = size(c,1);

    tCardiac = -log_files.relative_start_acquisition + ...
        ((0:(nSamplesCardiac-1))*dtCardiac)';

    tRespiration = -log_files.relative_start_acquisition + ...
        ((0:(nSamplesRespiration-1))*dtRespiration)';

    if isHigherSamplingCardiac
        t = tCardiac;
        rInterp = interp1(tRespiration, r, t);

        if DEBUG
            fh = plot_interpolation(tRespiration, r, t, rInterp, ...
                {'respiratory', 'cardiac'});
            verbose.fig_handles(end+1) = fh;
        end
        r = rInterp;

    else
        t = tRespiration;
        cInterp = interp1(tCardiac, c, t);

        if DEBUG
            fh = plot_interpolation(tCardiac, c, t, cInterp, ...
                {'cardiac', 'respiratory'});
            verbose.fig_handles(end+1) = fh;
        end
        c = cInterp;

    end

else
    nSamples = size(c,1);
    t = -log_files.relative_start_acquisition + ((0:(nSamples-1))*dt)';
end

cpulse = [];
end

% local function to plot interpolation result
function fh = plot_interpolation(tOrig, yOrig, tInterp, yInterp, ...
    stringOrigInterp)
fh = tapas_physio_get_default_fig_params;
stringTitle = sprintf('Read-In: Interpolation of %s signal', stringOrigInterp{1});
set(fh, 'Name', stringTitle);
plot(tOrig, yOrig, 'go--');  hold all;
plot(tInterp, yInterp,'r.');
legend({
    sprintf('after interpolation to %s timing', ...
    stringOrigInterp{1}), ...
    sprintf('original %s time series', stringOrigInterp{2}) });
title(stringTitle);
xlabel('time (seconds');
end

author: translationalneuromodeling

repo: tapas


text:
function MakeTableOfFunctions()
% MakeTableOfFunctions Outputs a csv summarizing INP_mops.txt
% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

% Read in INP_mops.txt to get a list of functions
fid = fopen('INP_mops.txt');
Mops = textscan(fid,'%s %s','CommentStyle','#','CollectOutput',1);
fclose(fid);
Mops = Mops{1};
Mops = Mops(:,1);

% Split at the parenthesis:
MopsSplit = regexp(Mops,'(','split','ignorecase');
mFiles = cellfun(@(x)x{1},MopsSplit,'UniformOutput',false);
mFiles = unique(mFiles);
numMFiles = length(mFiles);

% Get help lines as descriptions:
helpLines = cell(numMFiles,1);
for i = 1:numMFiles
    helpMe = help(mFiles{i});
    newLines = regexp(helpMe,'\n');
    theTopLine = helpMe(1:newLines(1)-1);
    notSpaces = regexp(theTopLine,'\w*');
    helpLines{i} = theTopLine(notSpaces(2):end);
end

% Make a table
T = table(mFiles,helpLines);

% Output to csv:
writetable(T,'codeSummary.csv','Delimiter',',','QuoteStrings',true);

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [a1, b1, c1, a2, b2, c2] = aaad(A, B, C)
%AAAD  gives both solutions to the angle-angle-angle problem, in degrees.
%
%   AAAD(A, B, C) will result in NaNs if the existence condition
%   |pi - |A|-|B|| <= |C| <= pi - ||A| - |B||
%   is not met.
%
%   See also AAAD.

% Rody P.S. Oldenhuis
% Delft University of Technology
% oldenhuis@gmail.com
%
% Crated     : 23/Feb/2009
% Last edited: 30/Nov/2012

    % find both solutions by calling aaa directly
    r2d = 180/pi;
    d2r = 1/r2d;
    [a1, b1, c1, a2, b2, c2] = aaa(A*d2r, B*d2r, C*d2r);
    [a1, b1, c1, a2, b2, c2] = deal(a1*r2d, b1*r2d, c1*r2d, a2*r2d, b2*r2d, c2*r2d);

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% Discrete All-Pole (DAP) envelope (AR model)
%
% Estimate an AR model based on given spectral peaks following the method in [1].
%
% Input
%  sins    : [Hz;amp] [2xN] Spectral peaks with frequency and linear amplitudes.
%            (as provided by sin_analysis.m)
%  fs      : [Hz] Sampling frequency
%  order   : Order of the AR model
% [dftlen] : DFT length of the envelope (if requested by the output)
% [opt]    : Additional options (see code below)
%
% Output
%  g  : gain factor
%  a  : AR coefficients
%  E  : Amplitude envelope (if dftlen is not empty)
%
% Reference
%  [1] El-Jaroudi, A., Makhoul, J.: Discrete All-Pole Modeling, IEEE Transactions
%      on Signal Processing 39(2), 411–423, 1991.
%
% Copyright (c) 2011 University of Crete - Computer Science Department
%
% License
%  This file is under the LGPL license,  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. This file 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.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function [g a E opt] = env_dap(sins, fs, order, dftlen, opt)

    if nargin<5
        opt.maxit     = 50;     % Maximum number of iterations allowed
        opt.alpha     = 0.5;    % Convergence speed
        opt.dISthresh = 10e-7;  % Stopping threshold for Itakura-Saito err diff
        opt.minbw     = [];     % [Hz] Minimum bandwidth of the poles.
                                %    Slow down the computation significantly !
                                %    It should be at least zero in order to ensure
                                %    the stability of the AR filter.
                                %    Advised to set it to 50Hz

        opt.debug = 0;          % 1: show only final solution
                                % 2: show each step
    end
    if nargin==0; g=opt; return; end

    if ~isempty(opt.minbw); minrr=bandwidth2rootradius(opt.minbw, fs); end

    wm = 2*pi*sins(1,:)'/fs; % Radian frequencies
    P = abs(sins(2,:));      % Magnitues
    N = size(sins,2);        % Number of Peaks

    B = exp(-j*wm*double(0:order));  % The Fourier basis
    Binv = exp(j*wm*double(0:order));% The inverse Fourier basis

    if opt.debug>1
        subplot(211);
            hold off;
            plot(fs*wm/(2*pi), log(P), 'xk');
            hold on;
    end

    % Compute the target AR from the power-spectrum P^2, eq. (5) in [1]
    % Surely a mistake in [1] because the autocorr is the inv Fourier transform
    % of the POWER spectrum and not that of the magnitude.
    R = (1/N)*real(P.^2*Binv);
    Rmatinv = inv(toeplitz(R)); % Prepare the inverse of the AR matrix

    % Initial guess using the Levinson-Durbin recursion
    [a,e] = levinson(R,order);
    a = (1/sqrt(e))*a'; % Put the gain factor into the AR coefficients

    Phat = 1./abs(B*a); % Get the model amplitude response
    if opt.debug>1; subplot(211); plot(fs*wm/(2*pi), log(Phat), 'g'); end

    err = nan(1,opt.maxit);
    it=1;
    cont = true;
    while cont && it<opt.maxit
        if opt.debug>1
            subplot(211);
                plot(fs*wm/(2*pi), log(abs(Phat)));
                F = fs*(0:dftlen/2)/dftlen;
                E = gba2hspec(1, 1, a, dftlen);
                plot(F, log(abs(E)), 'r');
                xlim([0 fs/2]);
                xlabel('Frequency [Hz]');
                ylabel('Amplitude [log]');

            subplot(212);
                plot(log10(err));
                xlabel('Iterations');
                ylabel('Itakura-Saito error [log10]');
            pause
        end

        A = B*a; % Compute the frequency response of the model

        h = (1/N)*real(B.'*(1./A)); % Compute the reversed impulse response

        anew = a*(1-opt.alpha) + opt.alpha * (Rmatinv*h); % Update the AR coeffs

        % If asked, ensure stability of the AR filter
        % by limiting the bandwidth of the poles
        if ~isempty(opt.minbw); anew=polystab2(anew, minrr); end

        % Compute the Itakura-Saito error
        Phat = 1./abs(B*anew); % Get the model amplitude response
        err(it)  =  1/N *sum(P(:)./Phat-log(P(:)./Phat)-1); % eq. (14) in [1]

        if it>2 % Start at 3 in order to skip the influence of the initial guess
            cont = err(it-1)-err(it)>opt.dISthresh;
        end

        a = anew;

        it = it + 1;
    end

    % Put the gain back to a separate parameter in order to have a0=1
    g = 1./a(1);
    a = a./a(1);

    if nargout>2
        if isempty(dftlen); dftlen=4096; end
        E = gba2hspec(g, 1, a, dftlen);
    end

    if opt.debug==1
        hold off;
        plot(sins(1,:), log(sins(2,:)), 'xk');
        hold on;
        F = fs*(0:dftlen/2)/dftlen;
        E = gba2hspec(g, 1, a, dftlen);
        plot(F, log(abs(E)));
        xlim([0 fs/2]);
        keyboard
    end

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function voxels = mask2voxel(mask)
% convert from 3-D mask to voxel list in canonical orientation
%
% :Usage:
% ::
%
%     function voxels = mask2voxel(mask)
%
% [i j k] = row, column, slice
%
% [x y z] in brain if brain is in analyze format
%
% (x is rows, y is columns, z is slices)
%
% ..
%    Tor Wager, 10/17/01
%
%    Refactored to be orders of magnitude faster
%    Matthew Davidson, 6/1/2006
% ..

% ..
%     Loop below is very slow. The find() function is much faster.
%
%     voxels = [];
%     index = 1;
%     for i = 1:size(mask,1)
%         for j = 1:size(mask,2)
%             for k = 1:size(mask,3)
%                 if mask(i,j,k) > 0, voxels(index,:) = [i j k];,index = index+1;,end
%             end
%         end
%     end
% ..

    [X Y Z] = ind2sub(size(mask), find(mask > 0));
    voxels = [X Y Z];
    voxels = sortrows(voxels);

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function fig = scopegui()
% SCOPEGUI is GUI-part for oscilloscope SCOPE.M function
% This is the machine-generated representation of a Handle Graphics object
% and its children.  Note that handle values may change when these objects
% are re-created. This may cause problems with any callbacks written to
% depend on the value of the handle at the time the object was saved.
%
% To reopen this object, just type the name of the M-file at the MATLAB
% prompt. The M-file and its associated MAT-file must be on your path.

load scopegui

h0 = figure('Units','points', ...
	'BackingStore','off', ...
	'Color',[0.8 0.8 0.8], ...
	'Colormap',mat0, ...
	'DeleteFcn','scope h_Close;', ...
	'MenuBar','none', ...
	'Name','DAPScope', ...
	'NumberTitle','off', ...
	'PointerShapeCData',mat1, ...
	'Position',[87 137.25 546.75 324], ...
	'RendererMode','manual', ...
	'Tag','scope', ...
	'WindowButtonDownFcn','zoom xdown', ...
	'WindowButtonUpFcn','ones;');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'Callback','scope(''h_Start'');', ...
	'ListboxTop',0, ...
	'Position',[468.375 69.75 63.75 24.75], ...
	'String','Start', ...
	'Tag','h_Start');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'Callback','scope(''h_Stop'');', ...
	'Enable','off', ...
	'ListboxTop',0, ...
	'Position',[468.375 43.5 63.75 24], ...
	'String','Stop', ...
	'Tag','h_Stop');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'Callback','scope h_Close;', ...
	'ListboxTop',0, ...
	'Position',[468.375 17.25 63.75 24], ...
	'String','Close', ...
	'Tag','h_Close');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0 0.5 0.5], ...
	'Enable','inactive', ...
	'FontWeight','demi', ...
	'ForegroundColor',[1 1 1], ...
	'ListboxTop',0, ...
	'Position',[468 132 64.5 16.5], ...
	'String','360mks/cnt', ...
	'Style','edit', ...
	'Tag','h_Dwell');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0 0.5 0.5], ...
	'Enable','inactive', ...
	'FontWeight','demi', ...
	'ForegroundColor',[1 1 1], ...
	'ListboxTop',0, ...
	'Position',[466.875 246.75 66.75 15], ...
	'String','0', ...
	'Style','edit', ...
	'Tag','h_SpNum');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'Callback','scope(''h_DwTime'');', ...
	'ListboxTop',0, ...
	'Max',333, ...
	'Min',20, ...
	'Position',[470.25 114 60 14.25], ...
	'String','DwTime', ...
	'Style','slider', ...
	'Tag','h_DwTime', ...
	'Value',100);
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0 0 1], ...
	'Enable','inactive', ...
	'FontWeight','demi', ...
	'ForegroundColor',[1 1 0], ...
	'ListboxTop',0, ...
	'Position',[459 289.5 82.5 15.75], ...
	'String','LanscView 1.0', ...
	'Style','edit', ...
	'Tag','Label1Text');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'FontWeight','bold', ...
	'ListboxTop',0, ...
	'Position',[467.625 264 65.25 15.75], ...
	'String','Trace Num', ...
	'Style','text', ...
	'Tag','Label2Text');
h1 = axes('Parent',h0, ...
	'Units','points', ...
	'Box','on', ...
	'CameraUpVector',[0 1 0], ...
	'CameraUpVectorMode','manual', ...
	'Color',[0 0 0], ...
	'ColorOrder',mat2, ...
	'DrawMode','fast', ...
	'FontSize',8, ...
	'FontWeight','demi', ...
	'NextPlot','replacechildren', ...
	'Position',[142.5 27 312.75 268.5], ...
	'Tag','h_Axes', ...
	'XColor',[0 0 0], ...
	'XLim',mat3, ...
	'XLimMode','manual', ...
	'YColor',[0 0 0], ...
	'YLim',[-5000 5000], ...
	'YLimMode','manual', ...
	'ZColor',[0 0 0]);
h2 = text('Parent',h1, ...
	'Color',[0 0 0], ...
	'FontSize',8, ...
	'FontWeight','demi', ...
	'HandleVisibility','off', ...
	'HorizontalAlignment','center', ...
	'Position',[235.2568512568512 5140.056022408962 17.32050807568877], ...
	'String','Signal', ...
	'Tag','h_AxesText4', ...
	'VerticalAlignment','bottom');
set(get(h2,'Parent'),'Title',h2);
h2 = text('Parent',h1, ...
	'Color',[0 0 0], ...
	'FontSize',8, ...
	'FontWeight','demi', ...
	'HandleVisibility','off', ...
	'HorizontalAlignment','center', ...
	'Position',[235.2568512568512 -5560.224089635854 17.32050807568877], ...
	'String','count', ...
	'Tag','h_AxesText2', ...
	'VerticalAlignment','cap');
set(get(h2,'Parent'),'XLabel',h2);
h2 = text('Parent',h1, ...
	'Color',[0 0 0], ...
	'HandleVisibility','off', ...
	'HorizontalAlignment','right', ...
	'Position',[-9.666225666225643 6036.414565826328 17.32050807568877], ...
	'Tag','h_AxesText1', ...
	'UserData',[0 1024 -5000 5000], ...
	'Visible','off');
set(get(h2,'Parent'),'ZLabel',h2);
h2 = line('Parent',h1, ...
	'Color',[1 1 0], ...
	'EraseMode','background', ...
	'LineWidth',0.2, ...
	'Tag','h_AxesLine2', ...
	'XData',mat4, ...
	'YData',mat5);
h2 = line('Parent',h1, ...
	'Color',[0 1 0], ...
	'EraseMode','background', ...
	'LineWidth',0.2, ...
	'Tag','h_AxesLine1', ...
	'XData',mat6, ...
	'YData',mat7);
h2 = text('Parent',h1, ...
	'Color',[0 0 0], ...
	'HandleVisibility','off', ...
	'HorizontalAlignment','center', ...
	'Position',[83.25685125685125 -42.01680672269049 17.32050807568877], ...
	'Rotation',90, ...
	'Tag','h_AxesText3', ...
	'VerticalAlignment','baseline');
set(get(h2,'Parent'),'YLabel',h2);
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'Callback','scope h_Mode;', ...
	'ListboxTop',0, ...
	'Position',[468.375 185.25 63.75 15.75], ...
	'String','Switch View', ...
	'Tag','h_Mode');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0 0.5 0.5], ...
	'Enable','inactive', ...
	'FontWeight','demi', ...
	'ForegroundColor',[1 1 1], ...
	'ListboxTop',0, ...
	'Position',[16.5 250 60 15.5], ...
	'String','0', ...
	'Style','edit', ...
	'Tag','h_Inp1_ed');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'Callback','scope h_Inp1_sl;', ...
	'ListboxTop',0, ...
	'Max',100, ...
	'Position',[80.25 251.25 14.25 30.75], ...
	'String','Shift, mV', ...
	'Style','slider', ...
	'Tag','h_Inp1_sl');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.501960784313725 0.501960784313725 0.501960784313725], ...
	'FontWeight','bold', ...
	'ForegroundColor',[0 1 0], ...
	'ListboxTop',0, ...
	'Position',[16.5 285.25 75 15.75], ...
	'String','Channel 1', ...
	'Style','text', ...
	'Tag','Label2Text');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'FontWeight','bold', ...
	'ListboxTop',0, ...
	'Position',[15.75 267.75 60 15.75], ...
	'String','Input', ...
	'Style','text', ...
	'Tag','h_Inp1_lb');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'Enable','inactive', ...
	'FontWeight','bold', ...
	'ListboxTop',0, ...
	'Position',[16.5 232.25 60 15.75], ...
	'String','Scale', ...
	'Style','text', ...
	'Tag','h_Scale1_lb');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'Callback','scope h_Scale1_sl;', ...
	'ListboxTop',0, ...
	'Max',101, ...
	'Min',1, ...
	'Position',[79.5 215.25 14.25 30.75], ...
	'String','Shift, mV', ...
	'Style','slider', ...
	'Tag','h_Scale1_sl', ...
	'Value',1);
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0 0.5 0.5], ...
	'Enable','inactive', ...
	'FontWeight','demi', ...
	'ForegroundColor',[1 1 1], ...
	'ListboxTop',0, ...
	'Position',[16.5 214.75 60 15.5], ...
	'String','1', ...
	'Style','edit', ...
	'Tag','h_Scale1_ed');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'FontWeight','bold', ...
	'ListboxTop',0, ...
	'Position',[466.875 150.75 66.75 15.75], ...
	'String','Dwell Time', ...
	'Style','text', ...
	'Tag','Label2Text');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0 0.5 0.5], ...
	'Enable','inactive', ...
	'FontWeight','demi', ...
	'ForegroundColor',[1 1 1], ...
	'ListboxTop',0, ...
	'Position',[16.5 179.25 60 15.75], ...
	'String','0', ...
	'Style','edit', ...
	'Tag','h_Shift1_ed');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'Callback','scope h_Shift1_sl;', ...
	'ListboxTop',0, ...
	'Max',5000, ...
	'Min',-5000, ...
	'Position',[79.5 179.25 14.25 30.75], ...
	'String','Shift, mV', ...
	'Style','slider', ...
	'Tag','h_Shift1_sl');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'Enable','inactive', ...
	'FontWeight','bold', ...
	'ListboxTop',0, ...
	'Position',[16.5 197 60 15.75], ...
	'String','Shift, mV', ...
	'Style','text', ...
	'Tag','h_Shift1_lb');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'FontWeight','bold', ...
	'ListboxTop',0, ...
	'Position',[16.5 52.75 60 15.75], ...
	'String','Shift, mV', ...
	'Style','text', ...
	'Tag','h_Shift2_lb');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'Callback','scope h_Shift2_sl;', ...
	'ListboxTop',0, ...
	'Max',5000, ...
	'Min',-5000, ...
	'Position',[79.5 35.25 14.25 30.75], ...
	'String','Shift, mV', ...
	'Style','slider', ...
	'Tag','h_Shift2_sl');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0 0.5 0.5], ...
	'FontWeight','demi', ...
	'ForegroundColor',[1 1 1], ...
	'ListboxTop',0, ...
	'Position',[16.5 35.25 60 15.5], ...
	'String','0', ...
	'Style','edit', ...
	'Tag','h_Shift2_ed');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0 0.5 0.5], ...
	'FontWeight','demi', ...
	'ForegroundColor',[1 1 1], ...
	'ListboxTop',0, ...
	'Position',[16.5 70.5 60 15.75], ...
	'String','1', ...
	'Style','edit', ...
	'Tag','h_Scale2_ed');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'Callback','scope h_Scale2_sl;', ...
	'ListboxTop',0, ...
	'Max',101, ...
	'Min',1, ...
	'Position',[79.5 70.5 14.25 30.75], ...
	'String','Shift, mV', ...
	'Style','slider', ...
	'Tag','h_Scale2_sl', ...
	'Value',1);
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'FontWeight','bold', ...
	'ListboxTop',0, ...
	'Position',[16.5 88.25 60 15.75], ...
	'String','Scale', ...
	'Style','text', ...
	'Tag','h_Scale2_lb');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'FontWeight','bold', ...
	'ListboxTop',0, ...
	'Position',[16.5 123.75 60 15.75], ...
	'String','Input', ...
	'Style','text', ...
	'Tag','h_Inp2_lb');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.501960784313725 0.501960784313725 0.501960784313725], ...
	'FontWeight','bold', ...
	'ForegroundColor',[1 1 0], ...
	'ListboxTop',0, ...
	'Position',[16.5 141.5 75 15.75], ...
	'String','Channel 2', ...
	'Style','text', ...
	'Tag','Label2Text');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'Callback','scope h_Inp2_sl;', ...
	'ListboxTop',0, ...
	'Max',100, ...
	'Position',[79.5 107.25 14.25 30.75], ...
	'String','Shift, mV', ...
	'Style','slider', ...
	'Tag','h_Inp2_sl');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0 0.5 0.5], ...
	'FontWeight','demi', ...
	'ForegroundColor',[1 1 1], ...
	'ListboxTop',0, ...
	'Position',[17.25 106 60 15.75], ...
	'String','0', ...
	'Style','edit', ...
	'Tag','h_Inp2_ed');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.752941176470588 0.752941176470588 0.752941176470588], ...
	'FontWeight','bold', ...
	'ListboxTop',0, ...
	'Position',[467.625 221.25 65.25 15.75], ...
	'String','View Mode', ...
	'Style','text', ...
	'Tag','Label2Text');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0 0.5 0.5], ...
	'Enable','inactive', ...
	'FontWeight','demi', ...
	'ForegroundColor',[1 1 1], ...
	'ListboxTop',0, ...
	'Position',[466.875 204 66.75 15], ...
	'String','X-T', ...
	'Style','edit', ...
	'Tag','h_View_ed');
if nargout > 0, fig = h0; end

author: Sable

repo: mcbench-benchmarks


text:
% fbp2_example.m
% Example of how to use fbp2.m
% This compares "ordinary" FBP with FBP based on Mojette sampling.
% todo: this should eventually supercede fbp_fan_arc_example.m
% Copyright 2005-12-16, Jeff Fessler, University of Michigan

if ~isvar('sino'), printm 'sino'
	down = 2;
	ig = image_geom('nx', 512, 'ny', 504, 'fov', 500);
	ig = ig.downsample(down);

	% parallel-beam
	sg = sino_geom('par', 'nb', 888, 'na', 984, ...
		'orbit', 180, ... % usual case
		'strip_width', 'dr', ... % not important for FBP
		'dr', 541/949, 'offset_r', 0.25);
%		'orbit', 360, ... % 2008-10-14 just for testing
	sg = sg.downsample(down);

%	ell = [0 0 200 200 0 1; 0 0 10 10 0 0];
%	ell = [30 0 120 100 0 1; 0 0 10 10 0 0];
%	ell = [20 -15 225 150 30 10; 40 0 10 15 0 1];
	ell = []; clim = (1 + [-1 1] * 0.05) * 1000;
%	ell = [100 70 f.ds/2 f.ds/2 0 100]; % point source

	[xtrue, ell] = ellipse_im(ig, ell, 'oversample', 4, 'rot', 90, ...
		'hu_scale', 1000);
	sino = ellipse_sino(sg, ell, 'oversample', 4);

	im plc 2 4
	im(1, xtrue, 'xtrue', clim), cbar
	im(5, sg.s, sg.ad, sino, 'sino'), cbar

	t = floor(min(sg.nb * sg.d, ig.fov)/2);
	ig.mask = ellipse_im(ig, [0 0 t t 0 1]) > 0;
	im(8, ig.mask, 'mask'), cbar
prompt
end

%
% conventional FBP reconstruction
%
if ~isvar('r_std'), printm 'fbp std'
	fg.std = fbp2(sg, ig);
	cpu etic
	r_std = fbp2(sino, fg.std);
	cpu etoc 'fbp std recon time'

	im(2, r_std, 'FBP matlab std', clim), cbar
	im(6, r_std - xtrue, 'error'), cbar
prompt
end

if 0 % compare single thread time
	fg.std1 = fbp2(sg, ig, 'nthread', 1);
	cpu etic
	r_std = fbp2(sino, fg.std1);
	cpu etoc 'fbp std recon time, nthread=1'
return
end

% dsc is too slow to be worth it
if 0 && ~isvar('r_dsc'), printm 'fbp dsc'
	fg.dsc = fbp2(sg, ig, 'type', 'dsc');
	cpu etic
	r_dsc = fbp2(sino, fg.dsc);
	cpu etoc 'fbp dsc recon time'

	im(3, r_dsc, 'FBP matlab dsc', clim), cbar
	im(6, r_dsc - xtrue, 'error'), cbar
prompt
end

%
% Mojette FBP reconstruction
% (an experimental approach - not essential)
%
if has_mex_jf
	if ~isvar('r_moj'), printm 'fbp moj'
		cpu etic
		fg.moj = fbp2(sg, ig, 'type', 'mojette', 'nthread', 1);
		cpu etoc 'moj setup time'

		cpu etic
		r_moj = fbp2(sino, fg.moj);
		cpu etoc 'moj recon time'

		im(3, r_moj, 'FBP mojette', clim), cbar
		im(7, r_moj - xtrue, 'error'), cbar
	prompt
	end

	max_percent_diff(xtrue, r_std)
	max_percent_diff(xtrue, r_moj)
	printm('sums: %g %g', [sum(r_std(:)) sum(r_moj(:))] / sum(xtrue(:)))
	nrms(r_std, xtrue)
	nrms(r_moj, xtrue)

	if 0 % profiles
		ix=1:ig.nx; iy=1:ig.ny; ix=round(ig.nx*0.1973); ii=iy;
		clf, plot(ii, xtrue(ix,iy), 'y:', ...
			ii, r_std(ix,iy), 'c-', ii, r_moj(ix,iy), 'g--')
		axis([[0.14 0.7]*ig.nx clim])
		legend('true', 'FBP', 'Mojette', 'location', 'south')
	return
	end
end

if ~has_aspire, return, end % check consistency with aspire

if ~isvar('r_asp')
	dir = test_dir;
	f.sino = [dir 'sino.fld'];
	f.image = [dir 'image.fld'];
	f.dsc = [dir 't.dsc'];
	fld_write(f.sino, sino)

	t = aspire_pair(sg, ig);
	char_array_write(f.dsc, t)
	f.win = 'boxcar,1,0,1';
	com = sprintf('echo y | i -chat 0 fbp dsc %s %s %s %s', ...
		f.image, f.sino, f.dsc, f.win);
%	eval(['!' com])
	cpu etic
	disp(os_run(com))
	cpu etoc 'aspire time'

	if 0 % compare filters
		tmp = fld_read('fft_filt.fld');
		sum(tmp) / sum(test)
%		clf, plot([tmp test])
		max_percent_diff(test, tmp)
	return
	end

	if 0 % compare filtered projections
		% seem to match except for slight shift??
		tmp = fld_read('proj_filt.fld');
		clf, plot([tmp(:,1) test(:,1)])
		plot(tmp(:,1)-test(:,1))
		max_percent_diff(test, tmp)
		sum(tmp(:)) / sum(test(:))
		minmax(test)
		minmax(tmp)
		minmax(test-tmp)
	return
	end

	r_asp = fld_read(f.image);
	im(4, r_asp, 'aspire', clim), cbar
	im(8, r_asp - xtrue, 'error'), cbar
	max_percent_diff(r_std, r_asp)
prompt
end

if 1 % show all images
	im plc 2 4; clim = [0.99 1.04]*1000; elim = [-0.1 0.1]*1000;
	im(1, xtrue, 'xtrue', clim), cbar
	im(5, sg.s, sg.ad, sino, 'sino'), cbar
	im(2, r_std, 'FBP matlab', clim), cbar
	im(6, r_std - xtrue, 'error', elim), cbar
	if has_mex_jf
	im(3, r_moj, 'FBP mojette', clim), cbar
	im(7, r_moj - xtrue, 'error', elim), cbar
	im(4, r_asp, 'FBP aspire', clim), cbar
	im(8, r_asp - xtrue, 'error', elim), cbar
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%% Example
% real time control of the KUKA iiwa 7 R 800
% Moving EEF of the robot in Z direction using a sinisoidal function

% An example script, it is used to show how to use the different
% functions of the KUKA Sunrise matlab toolbox

% First run the following script in Matlab
% Then start the client on the KUKA iiwa controller

% Mohammad SAFEEA, 24th of April 2018

close all,clear all;clc;

warning('off')

ip='172.31.1.147'; % The IP of the controller
% start a connection with the server
global t_Kuka;
t_Kuka=net_establishConnection( ip );

if ~exist('t_Kuka','var') || isempty(t_Kuka) || strcmp(t_Kuka.Status,'closed')
  warning('Connection could not be establised, script aborted');
  return;
end
% Joints positions
jPos={0,0,0,-pi/2,0,pi/2,0};
%setBlueOn(t); % turn on blue light
relVel=0.15;
movePTPJointSpace( t_Kuka , jPos, relVel); % move to initial configuration
%% Get Cartesian position of EEF
fprintf('Cartesian position')
eefpos=getEEFPos( t_Kuka );
eefposDist=eefpos;
disp(eefpos)
%% Start direct servo in Cartesian space
realTime_startDirectServoCartesian(t_Kuka);
disp('Starting direct servo in Cartesian space');
w=2; % motion constants, frequency rad/sec
A=75; % motion constants, amplitude of motion (mm)
tic;
deltaT=0;
counter=0;
initiationFlag=0;
disp('Enter control loop, stream EEF positions')
try
    %% Control loop
    while(deltaT<(6*pi/w))
        if(initiationFlag==0)
            initiationFlag=1;
            t_0=toc;
            t0=t_0;
        else
            time=toc;
            deltaT=time-t0;
            %%%%%%%%%%%%%%%%%%%%%%%%
            %% Perform trajectory calculation here
            eefposDist{3}=eefpos{3}-A*(1-cos(w*deltaT));
            %%%%%%%%%%%%%%%%%%%%%%%%
            % Send EEF position to robot
            if(toc-t_0>0.003)
                counter=counter+1;
                sendEEfPositionf( t_Kuka ,eefposDist);
                t_0=toc;
            end
        end
    end
    tstart=t0;
    tend=time;
    rate=counter/(tend-tstart);
    %% Stop the direct servo motion
    realTime_stopDirectServoCartesian( t_Kuka );
    fprintf('\nThe rate of update per second is: \n');
    disp(rate);
    fprintf('\n')
    pause(2);
    %% turn off light
    %setBlueOff(t);
    %% turn off the server
    net_turnOffServer( t_Kuka );
    fclose(t_Kuka);
    warning('on')
catch
    %% turn off the server
    net_turnOffServer( t_Kuka );
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function fieldnameIndex = compileFieldnameIndex(matFiles)
% Traverse the extracts of all matfiles, compiling a list of all fields and
% indexing from fieldnames to the matfile num.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

disp('planDB>> Compiling fieldname Index...');
fieldnameIndex = struct('fieldname', {},'planIndices',{});
for i=1:length(matFiles)
    if matFiles(i).isPlanC
        newfieldnames = listAllFieldNames(matFiles(i).extract);
        oldfieldnames = [fieldnameIndex.fieldname]';
        [C, INew, IOld] = intersect(newfieldnames, oldfieldnames);

        for j=1:length(IOld)
            fieldnameIndex(IOld(j)).planIndices(end+1) = i;
        end

        if isempty(oldfieldnames)
            oldfieldnames={};
        end

        [C, INew] = setdiff(newfieldnames, oldfieldnames);

        for j=1:length(INew)
            fieldnameIndex(end+1).fieldname = newfieldnames(INew(j));
            fieldnameIndex(end).planIndices = i;
        end

    end
end

for i=1:length(fieldnameIndex)
    fieldnameIndex(i).fieldname = fieldnameIndex(i).fieldname{:};
end

disp('planDB>> Fieldname Index Complete.');

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function p = minkowski_sum(a, b)

if size(a, 2) == 1 || size(b, 2) == 1
  p = bsxfun(@plus, a, b);
else
  p = zeros(2, size(a, 2) * size(b, 2));
  idx = 0;
  for j = 1:size(a, 2)
    p(:,idx+(1:size(b,2))) = bsxfun(@plus, a(:,j), b);
    idx = idx + size(b,2);
  end
  k = convhull(p(1,:), p(2,:), 'simplify', true);
  p = p(:,k(1:end-1));
end

author: rdeits

repo: iris-distro

save_path: github-repos/MATLAB/rdeits-iris-distro


text:
function c8_cube_root_test ( )

%*****************************************************************************80
%
%% C8_CUBE_ROOT_TEST tests C8_CUBE_ROOT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2013
%
%  Author:
%
%    John Burkardt
%
  test_num = 10;
  seed = 123456678;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8_CUBE_ROOT_TEST\n' );
  fprintf ( 1, '  C8_CUBE_ROOT computes the principal cube root of a C8.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '            C1=random            C2=C8_CUBE_ROOT(C1)         C3=C2*C2*C2\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    [ c1, seed ] = c8_uniform_01 ( seed );
    c2 = c8_cube_root ( c1 );
    c3 = c1 * c1 * c1;

    fprintf ( 1, '  %10f  %10f    %10f  %10f    %10f  %10f\n', ...
      real ( c1 ), imag ( c1 ), ...
      real ( c2 ), imag ( c2 ), ...
      real ( c3 ), imag ( c3 ) );

  end

  return
end


text:
function [lossValue, gradients, predictions] = vl_structuredNetwork_pairwiseModel(net, x, gradients, labels, dzdy, varargin)
%vl_structuredNetwork_pairwiseModel implements the evaluation of the structured network together with its gradient comuptation

opts = struct;
opts.conserveMemory = true ;
opts.sync = true ;
opts.disableDropout = false ;
opts.backPropagateType = 'all'; % 'all', 'unaryAndPairwise', 'onlyUnary', 'onlyPairwise'
opts.computeMaxMarginals = false;
opts.disableDropoutFeatureExtractor = false;
opts.lossNormalization = 1;
opts = vl_argparse(opts, varargin);

if (nargin <= 4) || isempty(dzdy)
    doder = false ;
else
    doder = true ;
end

%% preparation
numImages = length(labels);
numNodes = zeros(numImages, 1);
numEdges = zeros(numImages, 1);
for iImage = 1 : numImages
    numNodes( iImage ) = length( labels{iImage}.candidateBatchIds );
    numEdges( iImage ) = size(labels{ iImage }.clusteredEdges.bbIds, 1);
end
numCandidates = sum(numNodes);
if size(x, 4) ~= numCandidates
    error('vl_structuredNetwork_pairwiseModel:badSizeInput', 'Size of data input "x" is incompatible with "labels"');
end

n = numel(net.layers) ;
if doder
    if isempty(gradients)
        gradients = cell( length(net.layers), 1);
    end
    for i = 1 : n
        if isfield( net.layers{i}, 'weights' )
            for j = 1 : length( net.layers{i}.weights )
                gradients{i}.dzdw{j} = zeros( size(net.layers{i}.weights{j}), 'like', net.layers{i}.weights{j} );
            end
        end
    end
end

%% forward pass
% forward pass on the feature extraction network
disableFeatureExtractorDropout = opts.disableDropout || opts.disableDropoutFeatureExtractor;
numLayers_features = length(net.featureExtractor.layers);
res_features = initRes( numLayers_features );
res_features = vl_simplenn_pairwiseModel_forwardPass(net.featureExtractor, x, net.layers, res_features, ...
    'conserveMemory', opts.conserveMemory, ...
    'sync', opts.sync, ...
    'disableDropout', disableFeatureExtractorDropout);

extractedFeatures = res_features(end).x;
numFeatures = size( extractedFeatures, 3 );

% forward pass on the unary potential network
numLayers_unary = length(net.unaryNetwork.layers);
res_unary = initRes( numLayers_unary );
res_unary = vl_simplenn_pairwiseModel_forwardPass( net.unaryNetwork, extractedFeatures, net.layers, res_unary, ...
    'conserveMemory', opts.conserveMemory, ...
    'sync', opts.sync, ...
    'disableDropout', opts.disableDropout);

% forward pass on the pairwise potential network
x_pairwise = zeros( 1, 1, numFeatures * 2, sum(numEdges), 'like', x);
startId = 0;
batchId_inEdge = zeros( sum(numEdges), 2 );
for iImage = 1 : numImages
    curIds = startId + (1 : numEdges(iImage));
    curCandidateBatchIds = labels{ iImage }.candidateBatchIds( labels{ iImage }.clusteredEdges.bbIds );
    x_pairwise(:,:,1:numFeatures,curIds) = extractedFeatures(:,:,:, curCandidateBatchIds(:,1) );
    x_pairwise(:,:,numFeatures + 1 : 2 * numFeatures,curIds) = extractedFeatures(:,:,:, curCandidateBatchIds(:,2) );

    batchId_inEdge( curIds, :) = curCandidateBatchIds;

    startId = startId + numEdges(iImage);
end

numLayers_pairwise = length(net.pairwiseNetwork.layers);
res_pairwise = initRes( numLayers_pairwise );
res_pairwise = vl_simplenn_pairwiseModel_forwardPass( net.pairwiseNetwork, x_pairwise, net.layers, res_pairwise, ...
    'conserveMemory', opts.conserveMemory, ...
    'sync', opts.sync, ...
    'disableDropout', opts.disableDropout);

%% loss evaluation
unaryPotentials = gather(res_unary(end).x);
pairwisePotentials = gather(res_pairwise(end).x);

switch net.lossLayer.type
    case 'logisticScoresCompact'
        [lossValue, unaryDerivative, pairwiseDerivative, predictions] = vl_logisticScoreLoss_pairwiseCompactModel( unaryPotentials, pairwisePotentials, labels, dzdy, opts.computeMaxMarginals);
    case 'svmStructCompact'
        [lossValue, unaryDerivative, pairwiseDerivative, predictions] = vl_svmStructLoss_pairwiseCompactModel( unaryPotentials, pairwisePotentials, labels, dzdy, opts.computeMaxMarginals, opts.lossNormalization);
    otherwise
        error( ['Unknown structured loss: ', net.lossLayer.type] );
end

%% backward pass
if doder
    if isequal( opts.backPropagateType, 'all' ) || isequal( opts.backPropagateType, 'onlyUnary' ) || isequal( opts.backPropagateType, 'unaryAndPairwise' )
        % backward pass on the unary potential network
        [res_unary, gradients] = vl_simplenn_pairwiseModel_backwardPass( net.unaryNetwork, extractedFeatures, net.layers, res_unary, gradients, unaryDerivative, ...
            'conserveMemory', opts.conserveMemory, ...
            'sync', opts.sync, ...
            'disableDropout', opts.disableDropout );
    end

    if isequal( opts.backPropagateType, 'all' )
        % backward pass of the unary signal through the feature extracting network
        % start accumulating data to backprop through the feature extraction network:
        featureBackPropSignal = res_unary(1).dzdx;
    end

    if isequal( opts.backPropagateType, 'all' ) || isequal( opts.backPropagateType, 'onlyPairwise' ) || isequal( opts.backPropagateType, 'unaryAndPairwise' )
        % backward pass through pairwise potentials network
        [res_pairwise, gradients] = vl_simplenn_pairwiseModel_backwardPass( net.pairwiseNetwork, x_pairwise, net.layers, res_pairwise, gradients, pairwiseDerivative, ...
            'conserveMemory', opts.conserveMemory, ...
            'sync', opts.sync, ...
            'disableDropout', opts.disableDropout );
    end

    % prepare data for passing backward pairwise signal
    if isequal( opts.backPropagateType, 'all' )
        I = eye( numCandidates, 'like', x );
        ids1 = I(batchId_inEdge(:,1), :);
        ids2 = I(batchId_inEdge(:,2), :);

        data1 = res_pairwise(1).dzdx(:,:, 1 : numFeatures, :);
        data2 = res_pairwise(1).dzdx(:,:, numFeatures + 1 : 2 * numFeatures, :);

        data1 = reshape(data1, [], sum( numEdges ));
        data2 = reshape(data2, [], sum( numEdges ));

        featureBackPropSignal = featureBackPropSignal ...
            + reshape( data1 * ids1, size(featureBackPropSignal) ) ;
        featureBackPropSignal = featureBackPropSignal ...
            + reshape( data2 * ids2, size(featureBackPropSignal) ) ;

%         % SLOW VERSION
%         edgeOffset = 0;
%         for iImage = 1 : numImages
%             for iEdge = 1 : numEdges( iImage )
%                 instance1 = labels{iImage}.clusteredEdges.bbIds(iEdge, 1);
%                 instance1 = labels{iImage}.candidateBatchIds( instance1 );
%
%                 featureBackPropSignal(:,:,:,instance1) = featureBackPropSignal(:,:,:,instance1) ...
%                     + res_pairwise(1).dzdx(:,:, 1 : numFeatures, edgeOffset + iEdge);
%
%                 instance2 = labels{iImage}.clusteredEdges.bbIds(iEdge, 2);
%                 instance2 = labels{iImage}.candidateBatchIds( instance2 );
%
%                 featureBackPropSignal(:,:,:,instance2) = featureBackPropSignal(:,:,:,instance2) ...
%                     + res_pairwise(1).dzdx(:,:, numFeatures + 1 : 2 * numFeatures, edgeOffset + iEdge);
%             end
%             edgeOffset = edgeOffset + numEdges( iImage ) ;
%         end

        % backward pass of the pairwise signal through the feature extracting network
        [res_features, gradients] = vl_simplenn_pairwiseModel_backwardPass( net.featureExtractor, x, net.layers, res_features, gradients, featureBackPropSignal, ...
            'conserveMemory', opts.conserveMemory, ...
            'sync', opts.sync, ...
            'disableDropout', disableFeatureExtractorDropout );
    end
end

end

function res = initRes(n)
res = struct(...
    'x', cell(1,n+1), ...
    'dzdx', cell(1,n+1), ...
    'dzdw', cell(1,n+1), ...
    'aux', cell(1,n+1), ...
    'time', num2cell(zeros(1,n+1)), ...
    'backwardTime', num2cell(zeros(1,n+1))) ;
end

author: aosokin


text:
function mpe = find_mpe(engine, local_evidence)
% FIND_MPE Find the most probable explanation of the data
% function mpe = find_mpe(engine, local_evidence
%
% local_evidence{i}(j) = Pr(observation at node i | S(i)=j)
%
% This finds the marginally most likely value for each hidden node.
% It may give inconsistent results if there are ties.

[mpe, niter] = bp_mpe_mrf2(engine.mrf2.adj_mat, engine.mrf2.pot, local_evidence, ...
			   'max_iter', engine.max_iter, 'momentum', engine.momentum, ...
			   'tol', engine.tol);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [g,dgdx,dgdP] = g_GLM(x,P,u,in)

if isempty(u) || isequal(u,0)
    g = in.X*P;
    dgdx = [];
    dgdP = in.X';
else
    g = in.X(u,:)*P;
    dgdx = [];
    dgdP = in.X(u,:)';
end

author: MBB-team

repo: VBA-toolbox


text:
function value = isCategorical(tbl,name)
% Determine which columns in a table contain categorical variables.
%
% tbl  - The table
% name - Optional. Specify a specific column or cell array of columns to
%           check.
% OUTPUT
% value - logical matching the columns in tbl or the name cell-array
isCatColumn = @(x) (isa(x,'categorical') || iscellstr(x) || isstring(x) || ischar(x) || islogical(x));
value = logical(varfun(isCatColumn,tbl,'OutputFormat','Uniform'));
if nargin>1
    stay = ismember(tbl.Properties.VariableNames,name);
    value= value(stay);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


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

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

% Last Modified by GUIDE v2.5 17-Oct-2006 10:33:23

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

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

% --- Executes just before cluster_tool_table is made visible.
function cluster_tool_table_OpeningFcn(hObject, eventdata, handles, varargin)
movegui(hObject,'center');
set(handles.title,'string',[cluster_tool_getver]);
handles.pwd=pwd;handles.fdir=pwd;
for k=1:size(varargin(:),1)
    if strcmp(varargin{k},'cl'),handles.cl=varargin{k+1};
    elseif strcmp(varargin{k},'fdir'),handles.fdir=varargin{k+1};
    end
end
% s={'Print Labels (Talairach atlas--coordinates converted','from MNI using mni2tal)'};
% set(handles.tal_labels,'String',s);

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

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

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes cluster_tool_table wait for user response (see UIRESUME)
% uiwait(handles.figure1);

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

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in quit.
function quit_Callback(hObject, eventdata, handles)
delete(handles.figure1);

% --- Executes on button press in ret.
function ret_Callback(hObject, eventdata, handles)
if isfield(handles,'cl'),cluster_tool_tools('cl',handles.cl,'fdir',handles.fdir);
else cluster_tool_tools('fdir',handles.fdir);end
delete(handles.figure1);

% --- Executes on button press in write.
function write_Callback(hObject, eventdata, handles)

% --- Executes on button press in do_subcl.
function do_subcl_Callback(hObject, eventdata, handles)

% --- Executes on button press in do_labels.
function do_labels_Callback(hObject, eventdata, handles)
if get(handles.do_labels,'Value')
    set(handles.tal_labels,'Value',0)
end

% --- Executes on button press in done.
function done_Callback(hObject, eventdata, handles)
if get(handles.tal_labels,'Value')
    for i=length(handles.cl)
        handles.cl(i).XYZmm=mni2tal(handles.cl(i).XYZmm);
    end
end
s='cluster_table(handles.cl,get(handles.do_subcl,''Value''),any([get(handles.do_labels,''Value'') get(handles.tal_labels,''Value'')])';
if get(handles.write,'Value')
    handles.pwd=pwd;cd(handles.fdir);
    [FileName,PathName]=uiputfile('*.tsv','Save as tab seperated values...');
    cd(handles.pwd);
    if isempty(FileName)||~ischar(FileName),beep,disp('No file found'),return,end
    handles.fdir=PathName;
    s=[s ',''writefile'',[PathName FileName]'];
end
if get(handles.do_labels,'Value')||get(handles.tal_labels,'Value')
    try
        xyz=evalin('base','xyz');
        L3=evalin('base','L3');
        L5=evalin('base','L5');
        s=[s ',''tal_info'',xyz,L3,L5'];
    catch
        if get(handles.tal_labels,'Value')
            s=[s ',''talairach'''];
        end
    end
end
s=[s ');'];
eval(s)
beep,disp('Done printing table.')

function tal_labels_CreateFcn(hObject, eventdata,handles)
% s={'Print Labels (Talairach atlas--coordinates converted','from MNI using mni2tal)'};
% set(handles.tal_labels,'String',s)

% --- Executes on button press in tal_labels.
function tal_labels_Callback(hObject, eventdata, handles)
if get(handles.tal_labels,'Value')
    set(handles.do_labels,'Value',0)
end


text:
function [V,F3,F4]=loadawobj(modelname,opts)
% loadawobj
% Load an Wavefront/Alias obj style model. Will only consider polygons with
% 3 or 4 vertices.
% Programme will also ignore normal and texture data. It will also ignore any
% part of the obj specification that is not a polygon mesh, ie nurbs, these
% deficiencies can probably be remedied relatively easily
%
% [v,f3,f4]=loadawobj(modelname)
%
% To come, s=loadawobj(modelname) where s is a structure
% s should contain containing texture coords, vertex norms, vn, sub
% ranges g3 and g4, and possibly even faces > 4
%
% See also Anders Sandberg's vertface2obj.m and saveobjmesh.m
%
% there is still a good chance that obj files will not load. I
% would be grateful for any reports and examples of those that fail.
%
% W.S. Harwin, University Reading, 2006,2010.
% Matlab BSD license
% thanks also to Doug Hackett

version=0.2;
if nargin <1
  disp('specify model name')
end

fid = fopen(modelname,'r');
if (fid<0)
  error(['can not open file: ' modelname]);
  return ;
end

vnum=1;
f3num=1;
f4num=1;
vtnum=1;
vnnum=1;
g3num=1;
g4num=1;
Vtmp=[];

% Line by line passing of the obj file

while ~feof(fid)
  Ln=fgets(fid);
  Ln=removespace(Ln);
  objtype=sscanf(Ln,'%s',1);
  l=length(Ln);
  if l==0  % isempty(s) ;
%    disp(['empty' Ln]);
    continue
  end
%  disp(Ln)
  Lyn=Ln; %temp hack
  switch objtype
    case '#' % comment
      disp(Ln);
    case 'v' % vertex
      v=sscanf(Ln(2:end),'%f');
      Vtmp(:,vnum)=v;
      vnum=vnum+1;
    case 'vt'	% textures
      if vtnum==1
        disp(Ln);
        vtnum=vtnum+1;
      end
    case 'g' % sub mesh
      disp(Ln);
%      if f3num > g3num(end);
      g3num=[g3num f3num];
%      if f4num> g4num(end)
      g4num=[g4num f4num];
    case 'mtllib' % material library
        disp(Ln);
    case 'usemtl' % use this material name
        disp(Ln);
    case 'l' % Line
        disp(Ln);
    case 's' %smooth shading across polygons
        disp(Ln);
    case 'vn' % normals
      if vnnum==1
        disp(Ln);
        vnnum=vnnum+1;
      end
    case 'f' % faces
      nvrts=length(findstr(Ln,' ')); % spaces as a predictor of n vertices
      slashpat=findstr(Ln,'/');
      nslash=length(slashpat);
      if nslash >1 % dblslash can be 0, 1 or >1
        dblslash=slashpat(2)-slashpat(1); else dblslash=0;
      end
      Ln=Ln(3:end); % get rid of the f
      if nslash == 0 % Face = vertex
%        disp('-v');
        f1=sscanf(Ln,'%f');
      elseif nslash == nvrts && dblslash>1 % Face = v/tc
%        disp('-v/tc');
        data1=sscanf(Ln,'%f/%f');
        if nvrts == 3
          f1=data1([1 3 5]);
          tc1=data1([2 4 6]);
        end
        if nvrts == 4;
          f1=data1([1 3 5 7]);
          tc1=data1([2 4 6 8]);
        end
      elseif nslash == 2*nvrts && dblslash==1 % v//n
%        disp('-v//n');
        data1=sscanf(Ln,'%f//%f');
        if nvrts == 3
          f1=data1([1 3 5]);
          vn1=data1([2 4 6]);
          Vn3(:,f3num)=f1;
        end
        if nvrts == 4;
          f1=data1([1 3 5 7]);
          vn11=data1([2 4 6 8]);
          Vn4(:,f4num)=f1;
        end
      elseif nslash == 2*nvrts && dblslash>1 % v/tc/n
%        disp('-v/tc/n');
        data1=sscanf(Ln,'%f/%f/%f');
        if nvrts == 3
          f1=data1([1 4 7]);
          tc1=data1([2 5 8]);
          vn1=data1([3 6 9]);
          Vn3(:,f3num)=f1;
        end
        if nvrts == 4;
          f1=data1([1 4 7 10]);
          tc1=data1([2 5 8 11]);
          vn1=data1([3 6 9 12]);
          Vn4(:,f4num)=f1;
        end
      end
% Now put the data into the array(s)
      if nvrts == 3
          F3(:,f3num)=f1;
          f3num=f3num+1;
      elseif nvrts ==4
          F4(:,f4num)=f1;
          f4num=f4num+1;
      else
          warning(sprintf('v nvrts=%d %s',nvrts, Ln));
      end
    otherwise
%      if ~strcmp(Lyn,char([13 10])) % carriage return
%        if ~ all(Lyn == ' ') % ignore lines with only space
          disp(['unprocessed-' Ln '-']); % see what has not been processed
%          double(Ln)
          %        end
%      end
    end

end

fclose(fid);

% plot if no output arguments are given
if nargout ==0
  if exist('F3','var')
    patch('Vertices',Vtmp','Faces',F3','FaceColor','g');
  end
  if exist('F4','var')
    patch('Vertices',Vtmp','Faces',F4','FaceColor','b');
  end
  axis('equal')
  clear Vtmp F3 F4
end

if nargout >=2
  V=Vtmp;
  if ~ exist('F3','var')
    warning('No 3 element faces')
    F3=[];
  end
  if nargout ==3
    if ~ exist('F4','var')
      warning('No 4 element faces')
      F4=[];
    end
  end
end

if nargout ==1
  V.v=Vtmp;
  V.f3=F3;
  if exist('F4','var')
    V.f4=F4;end
  V.g3=g3num;
  V.g4=g4num;
  if exist('Vn3','var')
    V.Vn3=Vn3;  end
  if exist('Vn4','var')
    V.Vn4=Vn4;  end
end

function Lyn=removespace(Lyn)
% A not an elegant way to remove
% surplus space
Lyn=strtrim(Lyn);
Lyn=strrep(Lyn,'       ',' '); % 8-2 .. 12-6
Lyn=strrep(Lyn,'    ',' '); % 5-2 6-3 4-1
Lyn=strrep(Lyn,'  ',' '); % 3-2 2-1
Lyn=strrep(Lyn,'  ',' ');
Lyn=strrep(Lyn,char([13 10]),''); % remove cr/lf
Lyn=strrep(Lyn,char([10]),''); % remove lf


text:
function aux_MLS(params, hParentFigure)
% function aux_MLS(params, hParentFigure);
%-------------------------------------------
% Function to plot the best fitting model of seismicity variation parameters using the maximum likelihood score determination
%
% Incoming variables:
% params        : all variables
% hParentFigure : Handle of the parent figure
%
% J.Woessner, woessner@seismo.ifg.ethz.ch
% last update: 04.11.02

% Get the axes handle of the plotwindow
axes(sv_result('GetAxesHandle', hParentFigure, [], guidata(hParentFigure)));
hold on;
% Select a point in the plot window with the mouse
[fX, fY] = ginput(1);
disp(['X: ' num2str(fX) ' Y: ' num2str(fY)]);
% Plot a small circle at the chosen place
plot(fX,fY,'ok');

% Get closest gridnode for the chosen point on the map
[fXGridNode fYGridNode,  nNodeGridPoint] = calc_ClosestGridNode(params.mPolygon, fX, fY);
plot(fXGridNode, fYGridNode, '*r');
hold off;

% Get the data for the grid node
mNodeCatalog_ = params.mCatalog(params.caNodeIndices{nNodeGridPoint}, :);
%%%% Doe: Determine next grid point and earthquakes associated with it %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

params.fTimePeriod = params.fTimePeriod/365;
% Split the gridpoint catalog according to the defined Splittime
[mFirstCatalog, mSecondCatalog, fFirstPeriodExact, fSecondPeriodExact, fFirstPeriod,...
        result.fSecondPeriod] = ex_SplitCatalog(mNodeCatalog_, params.fSplitTime, params.bTimePeriod,...
    params.fTimePeriod, params.bTimePeriod, params.fTimePeriod);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Start First Plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nModel = params.mValueGrid(nNodeGridPoint,6)
switch nModel
case 1
    [fProb_nochange, fBic_nochange] = calc_loglikelihood_nochange(mFirstCatalog, mSecondCatalog);
    sTxt_nochange = ['No change!']
case 2
    [fdM, fProb_dM, fBic_dM, mLikeli_dM] = calc_loglikelihood_dM(mFirstCatalog, mSecondCatalog);
    plot_mls1p(mLikeli_dM(:,1), mLikeli_dM(:,2), 1);
    fdM = fdM;
    sTxt_dM = ['Simple magnitude shift: dM = ' num2str(fdM) ' , '...
            ' Max. likelihood score = ' num2str(fProb_dM) ' , BIC = ' num2str(fBic_dM)]
case 3
    [fS, fProb_stretch, fBic_stretch, mLikeli_dS] = calc_loglikelihood_stretch(mFirstCatalog, mSecondCatalog);
    plot_mls1p(mLikeli_dS(:,1), mLikeli_dS(:,2), 2);
    fdS = fS;
    sTxt_stretch = ['Magnitude stretch: c = ' num2str(fS) ' , '...
            ' Max. likelihood score = ' num2str(fProb_stretch) ' , BIC = ' num2str(fBic_stretch)]
case 4
    [fFac, fProb_Rate, fBic_Rate, mLikeli_Rate] = calc_loglikelihood_rate(mFirstCatalog, mSecondCatalog);
    plot_mls1p(mLikeli_Rate(:,1), mLikeli_Rate(:,2), 3);
    fRf = fFac;
    sTxt_rate = ['Rate change: R_f = ' num2str(fFac) ' , '...
            ' Max. likelihood score = ' num2str(fProb_Rate) ' , BIC = ' num2str(fBic_Rate)]
case 5
    [fdM_rate, fdM_Fac, fProb_dMrate, fBic_dMrate, mLikeli_dMrate] = calc_loglikelihood_dM_rate(mFirstCatalog, mSecondCatalog);
    plot_mls2p(mLikeli_dMrate(:,1), mLikeli_dMrate(:,2), mLikeli_dMrate(:,3),3);
    fdM = fdM_rate;
    fRf = fdM_Fac;
    sTxt_dMrate = ['Magnitude shift and rate change: dM = ' num2str(fdM_rate) ' , R_f =' num2str(fdM_Fac) ' , '...
            ' Max. likelihood score = ' num2str(fProb_dM) ' , BIC = ' num2str(fBic_dMrate)]
case 6
    [fdM_st, fStretch, fProb_Trans, fBic_Trans, mLikeli_Trans] = calc_loglikelihood_Trans(mFirstCatalog, mSecondCatalog);
    plot_mls2p(mLikeli_Trans(:,1), mLikeli_Trans(:,2), mLikeli_Trans(:,3), 1);
    fdM = fdM_st;
    fdS = fStretch;
    sTxt_Trans = ['Magnitude transformation: c = ' num2str(fStretch) ' , dM_st = ' num2str(fdM_st) ' , '...
            ' Max. likelihood score = ' num2str(fProb_Trans) ' , BIC = ' num2str(fBic_Trans)]
case 7
    [fdS_rate, fdS_Fac, fProb_dSrate, fBic_dSrate mLikeli_dSrate] = calc_loglikelihood_stretch_rate(mFirstCatalog, mSecondCatalog);
    plot_mls2p(mLikeli_dSrate(:,1), mLikeli_dSrate(:,2), mLikeli_dSrate(:,3), 2);
    fdS = fdS_rate;
    fRf = fdS_Fac;
    sTxt_dSrate = ['Stretch and rate change: c = ' num2str(fdS_rate) ' , R_f =' num2str(fdS_Fac) ' , '...
            ' Max. likelihood score = ' num2str(fProb_dSrate) ' , BIC = ' num2str(fBic_dSrate)]
case 8
    [fdM_all, fdS_all, fFac_all, fProb_all, fBic_all, mLikeli_all] = calc_loglikelihood_dMdSrate(mFirstCatalog, mSecondCatalog);
    plot_mls3p(mLikeli_all(:,1), mLikeli_all(:,2), mLikeli_all(:,3), mLikeli_all(:,4));
    fdM = fdM_all;
    fdS = fdS_all;
    fRf = fFac_all;
    sTxt_dSdMrate = ['Shift, Stretch and Rate change: dM = ' num2str(fdM_all) ' , c = ' num2str(fdS_all) ' , R_f =' num2str(fFac_all) ' , '...
            ' Max. likelihood score = ' num2str(fProb_all) ' , BIC = ' num2str(fBic_all)]
otherwise
    errordlg('No model chosen','Model error');
    return;
end
% Plot model
plot_svmodel(fdM, fdS, fRf, mFirstCatalog, mSecondCatalog);

% Determine all Bics and show
[fProb_nochange, fBic_nochange] = calc_loglikelihood_nochange(mFirstCatalog, mSecondCatalog);
[fdM, fProb_dM, fBic_dM, mLikeli_dM] = calc_loglikelihood_dM(mFirstCatalog, mSecondCatalog);
[fS, fProb_stretch, fBic_stretch, mLikeli_dS] = calc_loglikelihood_stretch(mFirstCatalog, mSecondCatalog);
[fFac, fProb_Rate, fBic_Rate, mLikeli_Rate] = calc_loglikelihood_rate(mFirstCatalog, mSecondCatalog);
[fdM_rate, fdM_Fac, fProb_dMrate, fBic_dMrate, mLikeli_dMrate] = calc_loglikelihood_dM_rate(mFirstCatalog, mSecondCatalog);
[fdS_rate, fdS_Fac, fProb_dSrate, fBic_dSrate mLikeli_dSrate] = calc_loglikelihood_stretch_rate(mFirstCatalog, mSecondCatalog);
[fdM_st, fStretch, fProb_Trans, fBic_Trans, mLikeli_Trans] = calc_loglikelihood_Trans(mFirstCatalog, mSecondCatalog);
[fdM_all, fdS_all, fFac_all, fProb_all, fBic_all, mLikeli_all] = calc_loglikelihood_dMdSrate(mFirstCatalog, mSecondCatalog);

vBic = [fBic_nochange; fBic_dM; fBic_stretch; fBic_Rate; fBic_dMrate;  fBic_Trans; fBic_dSrate; fBic_all]

author: CelsoReyes


text:
classdef TestTonemap
    %TestTonemap

    methods (Static)
        function test_1
            hdr = cv.imread(fullfile(mexopencv.root(),'test','memorial.hdr'), 'Flags',-1);
            hdr = cv.resize(hdr, 0.5, 0.5);

            obj = cv.Tonemap();
            obj.Gamma = 2.2;

            ldr = obj.process(hdr);
            validateattributes(ldr, {'single'}, ...
                {'size',size(hdr), 'nonnegative'});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function duration=getDuration(obj)
% get the duration of the test
%
% duration=getDuration(obj)
%
% Input:
%   obj                 MOxUnitTestReport instance
%
% Output:
%   duration            total time in seconds that it took to run all tests
%                       in the obj instance

    get_outcome_duration=@(i)getDuration(getTestOutcome(obj,i));

    all_durations=arrayfun(get_outcome_duration,...
                            1:countTestOutcomes(obj));

    duration=sum(all_durations);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% SP_PRECOMPUTE: compute the basis functions in the whole mesh. The output
%  can consume a lot of memory. Use with care.
%
%     space = sp_precompute (space, msh);
%     space = sp_precompute (space, msh, 'option', value);
%
% INPUT:
%
%    space: object representing the discrete function space (see sp_vector)
%    msh: mesh object containing the quadrature information (see msh_cartesian)
%   'option', value: additional optional parameters, currently available options are:
%
%              Name     |   Default value |  Meaning
%           ------------+-----------------+----------------------------------
%            value      |      true       |  compute shape_functions
%            gradient   |      false      |  compute shape_function_gradients
%            divergence |      false      |  compute shape_function_divs
%            curl       |      false      |  compute shape_function_curls
%
% OUTPUT:
%
%    space: object containing the information of the input object, plus the
%            fields of the old structure, that are listed below. If no option
%            is given all the fields are computed. If an option is given,
%            only the selected fields will be computed.
%
%    FIELD_NAME      (SIZE)                               DESCRIPTION
%    ncomp           (scalar)                              number of components of the functions of the space
%    ndof            (scalar)                              total number of degrees of freedom
%    ndof_dir        (ncomp x ndim matrix)                 for each component, number of degrees of freedom along each direction
%    nsh_max         (scalar)                              maximum number of shape functions per element
%    nsh             (1 x msh.nel vector)                  actual number of shape functions per each element
%    connectivity    (nsh_max x msh.nel vector)            indices of basis functions that do not vanish in each element
%    shape_functions (ncomp x msh.nqn x nsh_max x msh.nel) basis functions evaluated at each quadrature node in each element
%    shape_function_gradients
%             (ncomp x rdim x msh.nqn x nsh_max x msh.nel) basis function gradients evaluated at each quadrature node in each element
%    shape_function_divs (msh.nqn x nsh_max x msh.nel)     basis function divergence evaluated at each quadrature node in each element
%    shape_function_curls
%         2D:  (msh.nqn x nsh_max x msh.nel)               basis function curl evaluated at each quadrature node in each element
%         3D:  (3 x msh.nqn x nsh_max x msh.nel)
%
% Copyright (C) 2009, 2010 Carlo de Falco
% Copyright (C) 2015, 2019 Rafael Vazquez
%
%    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/>.

function sp_out = sp_precompute (sp, msh, varargin)

value = true;
gradient = false;
divergence = false;
curl = false;
hessian = false;

if (~isempty (varargin))
  if (~rem (length (varargin), 2) == 0)
    error ('sp_precompute: options must be passed in the [option, value] format');
  end
  for ii=1:2:length(varargin)-1
    if (strcmpi (varargin {ii}, 'value'))
      value = varargin {ii+1};
    elseif (strcmpi (varargin {ii}, 'gradient'))
      gradient = varargin {ii+1};
    elseif (strcmpi (varargin {ii}, 'curl'))
      curl = varargin {ii+1};
    elseif (strcmpi (varargin {ii}, 'divergence'))
      divergence = varargin {ii+1};
    elseif (strcmpi (varargin {ii}, 'hessian'))
      hessian = varargin {ii+1};
    else
      warning ('Ignoring unknown option %s', varargin {ii});
    end
  end
end

if (~isstruct (msh))
  msh = msh_precompute (msh);
end

grad_param = gradient || divergence || curl || hessian;
value_param = value || grad_param;
div_param = false; curl_param = false;
switch (lower (sp.transform))
  case {'curl-preserving'}
    curl_param = curl;
  case {'div-preserving'}
    div_param = divergence;
end

sp_out = sp_precompute_param (sp, msh, 'value', value_param, 'gradient', grad_param, 'divergence', div_param, 'curl', curl_param, 'hessian', hessian);

switch (lower (sp.transform))
  case {'grad-preserving'}
    sp_out = sp_vector_grad_preserving_transform (sp_out, msh, value, gradient, curl, divergence, hessian);
  case {'curl-preserving'}
    sp_out = sp_vector_curl_preserving_transform (sp_out, msh, value, curl);
    if (gradient || divergence || hessian)
      warning ('Gradient, divergence and hessian not implemented for curl-preserving transformation')
    end
  case {'div-preserving'}
    sp_out = sp_vector_div_preserving_transform (sp_out, msh, value, gradient, curl, divergence);
    if (hessian)
      warning ('Hessian not implemented for div-preserving transformation')
    end
end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [W] = spm_Volt_W(u)
% returns basis functions used for Volterra expansion
% FORMAT [W] = spm_Volt_W(u);
% u  - times {seconds}
% W  - basis functions (mixture of Gammas)
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_Volt_W.m 1143 2008-02-07 19:33:33Z spm $

u     = u(:);
W     = [];
for i = 2:4
    m   = (2^i);
    s   = sqrt(m);
    W   = [W spm_Gpdf(u,(m/s)^2,m/s^2)];
end

author: spm


text:
%DECOLOR  Transforms a color image to a grayscale image
%
%     [grayscale, color_boost] = cv.decolor(src)
%     [grayscale, color_boost] = cv.decolor(src, 'OptionName',optionValue, ...)
%
% ## Input
% * __src__ Input 8-bit 3-channel color image.
%
% ## Output
% * __grayscale__ Output 8-bit 1-channel grayscale image.
% * **color_boost** Output 8-bit 3-channel contrast-boosted color image. It is
%   obtained by converting `src` to Lab color space, replacing the lightness
%   channel with `grayscale`, and converting back to RGB color space.
%
% ## Options
% * __FlipChannels__ whether to flip the order of color channels in input
%   `src` and output `color_boost`, between MATLAB's RGB order and OpenCV's
%   BGR (input: RGB->BGR, output: BGR->RGB). default false
%
% It is a basic tool in digital printing, stylized black-and-white photograph
% rendering, and in many single channel image processing applications [CL12].
%
% This function is to be applied on color images.
%
% ## References
% [CL12]:
% > Cewu Lu, Li Xu, and Jiaya Jia. "Contrast preserving decolorization".
% > In Computational Photography (ICCP), 2012 IEEE International Conference
% > on, pages 1-7. IEEE, 2012.
%
% See also: cv.cvtColor
%


text:
function cl_ext_make_resid(con_files, varargin)
% This function will create residual images (4d) and mask image in
% a current or assigned directory in order to use them in estimating smoothness
% (relevant functions: spm_est_smoothness (SPM), 3dFWHMx (AFNI), smoothest (FSL).
%
% :Usage:
% ::
%
%     function cl_ext_make_resid(conimgs, varargin)
%
% :Inputs:
%
%   **con_files:**
%        contrast image file names; This could be a cell array or
%        strings. This could be 4d images.
%
%        Best: Input a cell string. e.g., for a string matrix:
%
%        Use cl_ext_make_resid(cellstr(imgs)); % save residual images
%
%   * If you are not providing the absolute paths of the images, you need to
%   be in the directory that has the image files.
%
% :Outputs:
%
%   **Res4d.nii:**
%        residual images saved by SPM.
%
%   **mask.nii:**
%        the mask image that was used.
%
% :Options for varargin:
%
%   **'mask'**
%        This option can be used to estimate a cluster size for the correction for multiple
%        comparisons "within the mask". You can put in a ROI mask or gray matter,
%        whatever. If you don't specify a mask image, brainmask.nii (default) will be
%        used, but the image has to be in your path.
%         e.g.)
%         ::
%
%              mask = fullfile(basedir, 'ROI_image.img');
%              mask = which('scalped_avg152T1_graymatter_smoothed.img'); % limited to gray matter
%
%   **'outputdir'**
%        With this option, this will save residual and mask images and in the
%        outputdir directory. If you don't give outputdir, the current directory
%        will be used (default).
%
% This function calls cl_ext_spm_spm.m, which is a modified spm_spm not to
% delete residual images.
%
% ..
%    Choong-Wan (Wani) Woo, 01/22/2013
%    modified by Wani, 05/18/2013
% ..

outputdir = pwd;
mask = [];

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            case 'mask', mask = varargin{i+1};
            case 'outputdir', outputdir = varargin{i+1};
        end
    end
end

temp_dir = fullfile(outputdir, 'cl_extent_estimation');
mkdir(temp_dir);

if isempty(mask)
    try
        mask = which('brainmask.nii');
    catch
        warning('We tried to use a default mask (brainmask.nii), but we couldn''t use it. Maybe you need to add "brainmask.nii" into your path. We are going to use the first contrast image as a mask.');
        mask = conimgs{1};
    end
end

matlabbatch{1}.spm.stats.factorial_design.dir{1} = temp_dir; % result dir

% file names should be in cell arrays.

if ~iscell(con_files)
    if size(con_files, 1) == 1
        con_files = filenames(con_files, 'absolute', 'char');
        con_files = expand_4d_filenames(con_files);
        for ii = 1:size(con_files,1), conimgs{ii} = deblank(con_files(ii,:)); end
    else
        for ii = 1:size(con_files,1)
            conimgs(ii) = filenames(deblank(con_files(ii,:)), 'absolute'); % wani added this line for efficiency
%             con_files_temp(ii,:) = filenames(con_files(ii,:), 'absolute', 'char');
        end
%         con_files = con_files_temp;
    end
%     for i = 1:size(con_files,1)
%         conimgs{i} = deblank(con_files(i,:));
%     end
    conimgs = conimgs';
else
    if length(con_files) == 1
        con_files = expand_4d_filenames(con_files);
        for i = 1:size(con_files,1)
            conimgs{i} = deblank(con_files(i,:));
        end
        conimgs = conimgs';
    else
        conimgs = con_files;
    end
end

con_num = length(conimgs);
if con_num == 1
    error('The number of contrast images is one. Please check the image names.');
end

matlabbatch{1}.spm.stats.factorial_design.des.t1.scans = conimgs; % each con image for 2nd level
% matlabbatch{1}.spm.stats.factorial_design.cov = struct('c', [], 'cname',[],'iCFI',[],'iCC',[]);

matlabbatch{1}.spm.stats.factorial_design.masking.tm.tm_none = 1;
matlabbatch{1}.spm.stats.factorial_design.masking.im = 0;

matlabbatch{1}.spm.stats.factorial_design.masking.em{1} = mask;

matlabbatch{1}.spm.stats.factorial_design.globalc.g_omit = 1;
matlabbatch{1}.spm.stats.factorial_design.globalm.gmsca.gmsca_no = 1;
matlabbatch{1}.spm.stats.factorial_design.globalm.glonorm = 1;

cd(temp_dir);
save spm_2nd matlabbatch;

ndf = [1 con_num-1];

if exist('SPM.mat', 'file'), delete('SPM.mat'); end
spm_jobman('run', 'spm_2nd.mat');

load SPM.mat;

try
    cl_ext_spm_spm(SPM); % you need cl_ext_spm_spm.m in your path
catch
    error('You need "cl_ext_spm_spm.m" in your path.');
end

mask_img = filenames('mask*.img','char');
Res_imgs = filenames('ResI_*.img', 'char');

dat_mask = fmri_data(mask_img);
dat_mask.fullpath = fullfile(outputdir, 'mask.nii');
write(dat_mask);

dat_res = fmri_data(Res_imgs);
dat_res.fullpath = fullfile(outputdir, 'Res4d.nii');
write(dat_res);

% rmdir(temp_dir, 's');

cd(outputdir);

end


text:
function h=m_range_ring(long,lat,range,varargin)
% M_RANGE_RING Creates range rings on a map
%    M_RANGE_RING(LONG,LAT,RANGE) creates a range ring of range RANGE
%    km centered at the position specified by LONG and LAT. Range rings
%    will generally appear as small (almost) circles for short ranges,
%    but will be distorted at longer ranges.
%
%    If RANGE is a vector, concentric rings at the specified ranges
%    are drawn. If LONG,LAT are vectors, rings are drawn around
%    all specified locations.
%
%    The appearance of lines can be modified using the usual
%    line properties thus:
%    M_RANGE_RING(LONG,LAT,RANGE, <line property/value pairs>)
%
%    Sometimes you may need to adjust the number of points plotted
%    in each range ring (this can happen if the ring is at the extreme
%    edge of certian projections). THis can be done using
%    M_RANGE_RING(LONG,LAT,RANGE,NPTS, <line property/value pairs>)
%
%    NB: Earth radius is assumed to 6378.137km (WGS-84 value), and
%    calculations are for spherical geometry.

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

% 6/Nov/00 - eliminate returned stuff if ';' neglected (thx to D Byrne)

global MAP_VAR_LIST

pi180=pi/180;
earth_radius=6378.137;
n=72;

if ~isempty(varargin)  &&  ~ischar(varargin{1})
 n=varargin{1};varargin(1)=[];
end

c=range(:)'/earth_radius;

h=[];
for k=1:length(long)
  rlat=lat(k)*pi180;
  rlong=long(k)*pi180;
  if long(k)<MAP_VAR_LIST.longs(1), rlong=rlong+2*pi; end
  if long(k)>MAP_VAR_LIST.longs(2), rlong=rlong-2*pi; end

  x=sin([0:n-1]'/(n-1)*2*pi)*c;
  y=cos([0:n-1]'/(n-1)*2*pi)*c;
  on=ones(n,1);

  Y=(asin(on*cos(c)*sin(rlat) + (on*cos(rlat)*(sin(c)./c)).*y))/pi180;
  switch lat(k)
    case 90
      X=(rlong+atan2(x,-y))/pi180;
    case -90
      X=(rlong+atan2(x,y))/pi180;
    otherwise
      X=(rlong+atan2(x.*(on*sin(c)),on*(cos(rlat)*cos(c).*c) - (on*sin(rlat)*sin(c)).*y ) )/pi180;
  end

  nz=zeros(1,length(range(:)));
  X=X+cumsum([nz;diff(X)<-300]-[nz;diff(X)>300])*360;

  kk=find(X(1,:)~=X(end,:));
  X2=X(:,kk)+360;X2(X2>MAP_VAR_LIST.longs(2))=NaN;
  X3=X(:,kk)-360;X3(X3<MAP_VAR_LIST.longs(1))=NaN;

  [XX,YY]=m_ll2xy([X,X2,X3],[Y,Y(:,kk),Y(:,kk)],'clip','on');

  % Get rid of 2-point lines (these are probably clipped lines spanning the window)
  fk=finite(XX(:));
  st=find(diff(fk)==1)+1;
  ed=find(diff(fk)==-1);
  if length(st)<length(ed), st=[1;st]; end
  if length(ed)<length(st), ed=[ed;length(fk)]; end
  k=find((ed-st)==1);
  XX(st(k))=NaN;

  h=[h;line(XX,YY,varargin{:},'tag','m_range_ring')];

end

if nargout==0
 clear h
end