text:
function iNewDipoles = import_dipoles(iStudy, DataFile, DipolesFiles)
% IMPORT_DIPOLES: Imports a dipoles file (Neuromag xfit or CTF).
%
% USAGE:
%         iNewDipoles = import_dipoles(iStudy, DataFile, DipolesFiles) : Read files and save them in brainstorm database
%         iNewDipoles = import_dipoles(iStudy, DataFile)               : Ask files to user, read it, and save it in database
%         iNewDipoles = import_dipoles(iStudy)                         : Same wihtout attaching dipoles to any data file
%
% INPUT:
%    - iStudy       : Index of the study where to import the DipolesFiles
%    - DipolesFiles : Full filename, or cell array of filenames, of the dipoles files to import
%                     => if not specified : file to import is asked to the user

% @=============================================================================
% 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, 2010-2012

%% ===== PARSE INPUTS =====
% Argument: DipolesFiles
if (nargin < 3)
    DipolesFiles = [];
elseif ~iscell(DipolesFiles)
    DipolesFiles = {DipolesFiles};
end
% Argument: DataFile
if (nargin < 3) || isempty(DataFile)
    DataFile = '';
end

% Detect file format
FileFormat = '';
if ~isempty(DipolesFiles)
    % Get the file extenstion
    [fPath, fBase, fExt] = bst_fileparts(DipolesFiles{1});
    if ~isempty(fExt)
        fExt = lower(fExt(2:end));
        % Detect file format by extension
        switch lower(fExt)
            case 'bdip', FileFormat = 'BDIP';
            case 'mat',  FileFormat = 'BST';
            case 'dip',  FileFormat = 'CTFDIP';
            otherwise,   error('Unsupported file format.');
        end
        % Display assumed file format
        disp(['Default file format for this extension: ' FileFormat]);
        disp('If you want to specify the extension, please run this function without arguments.');
    end
end

%% ===== SELECT DIPOLES FILE =====
% If file to load was not defined : open a dialog box to select it
if isempty(DipolesFiles)
    % Get default import directory and formats
    LastUsedDirs   = bst_get('LastUsedDirs');
    DefaultFormats = bst_get('DefaultFormats');
    % Get Dipoles file
    [DipolesFiles, FileFormat] = java_getfile('open', ...
            'Import dipoles...', ...   % Window title
            LastUsedDirs.ImportData, ...    % Last used directory
            'multiple', 'files', ...        % Selection mode
            {{'.bdip'},    'Xfit binary (*.bdip)',      'BDIP'; ...
             {'_dipoles'}, 'Brainstorm (dipoles*.mat)', 'BST'; ...
             {'.dip'},     'CTF dipoles (*.dip)',       'CTFDIP' ...
            }, DefaultFormats.DipolesIn);
    % If no file was selected: exit
    if isempty(DipolesFiles)
        return
    % Make sure it's an array of files
    elseif ~iscell(DipolesFiles)
        DipolesFiles = {DipolesFiles};
    end
    % Save default import directory
    LastUsedDirs.ImportData = bst_fileparts(DipolesFiles{1});
    bst_set('LastUsedDirs', LastUsedDirs);
    % Save default import format
    DefaultFormats.DipolesIn = FileFormat;
    bst_set('DefaultFormats',  DefaultFormats);
end

%% ===== GET INFORMATION =====
% Get various information
isProgressBar = bst_progress('isVisible');
ProtocolInfo = bst_get('ProtocolInfo');
iNewDipoles = [];
% Get output study
sStudy = bst_get('Study', iStudy);
% Get channel file
ChannelFile = bst_get('ChannelFileForStudy', sStudy.FileName);

%% ===== GET COORDINATES TRANSFORMATION =====
TransfMeg = [];
% Get coordinates transformations
if ~isempty(ChannelFile)
    % Load channel file transformations
    ChannelMat = in_bst_channel(ChannelFile, 'TransfMeg');
    % If transformation Neuromag Head CS -> SCS transformation is defined
    if isfield(ChannelMat, 'TransfMeg') && ~isempty(ChannelMat.TransfMeg) && (length(ChannelMat.TransfMeg) >= 2)
        TransfMeg = ChannelMat.TransfMeg{2};
        % Add other transformations if there are some
        for i = 3:length(ChannelMat.TransfMeg)
            TransfMeg = ChannelMat.TransfMeg{i} * TransfMeg;
        end
    else
        errMsg = ['Warning: No transformation defined in this channel file.' 10 ...
                  'The dipoles locations might not be registered with the anatomy.'];
    end
else
    errMsg = ['Warning: No channel file defined for this study.' 10 10 ...
              'You should import the sensors definition before importing the dipoles, ' 10 ...
              'if not the dipoles locations might not be registered with the anatomy.'];
end
% Warning if no transformation
if isempty(TransfMeg)
    % Ask user confirmation
    isConfirm = java_dialog('confirm', [errMsg 10 10 'Import dipoles anyway ?' 10 10], 'Import dipoles');
    % If cancelled: exit
    if ~isConfirm
        return
    end
end

% Loop on each file
for iFile = 1:length(DipolesFiles)
    %% ===== LOAD DIPOLES FILE =====
    % Progress bar
    if ~isProgressBar
        bst_progress('start', 'Import dipoles', ['Loading file "' DipolesFiles{iFile} '"...']);
    end
    % Load file
    switch FileFormat
        case 'BDIP'
            DipolesMat = in_dipoles_bdip(DipolesFiles{iFile});
        case 'BST'
            DipolesMat = load(DipolesFiles{iFile});
        case 'CTFDIP'
            DipolesMat = in_dipoles_ctfdip(DipolesFiles{iFile});
        otherwise
            error('Unsupported file format.');
    end
    % Add data file
    DipolesMat.DataFile = DataFile;
    % Check that something was read
    if isempty(DipolesMat)
        error('Dipoles file could not be imported.');
    end
    % Add History field
    DipolesMat = bst_history('add', DipolesMat, 'import', ['Import from: ' DipolesFiles{iFile}]);

    %% ===== CHANGE COORDINATE SYSTEM =====
    if ~isempty(TransfMeg)
        R = TransfMeg(1:3,1:3); % Rotation
        T = TransfMeg(1:3,4);   % Translation
        % Add History field
        DipolesMat = bst_history('add', DipolesMat, 'import', sprintf('Rotation: [%1.3f,%1.3f,%1.3f; %1.3f,%1.3f,%1.3f; %1.3f,%1.3f,%1.3f]', R'));
        DipolesMat = bst_history('add', DipolesMat, 'import', sprintf('Translation: [%1.3f,%1.3f,%1.3f]', T));
        % Process each dipole
        for iDip = 1:length(DipolesMat.Dipole)
            % Position
            DipolesMat.Dipole(iDip).Loc = R * DipolesMat.Dipole(iDip).Loc + T;
            % Orientation (amplitude in each direction)
            DipolesMat.Dipole(iDip).Amplitude = R * DipolesMat.Dipole(iDip).Amplitude;
        end
    end

    %% ===== SAVE NEW FILE =====
    % Get imported base name
    [tmp__, importedBaseName, importedExt] = bst_fileparts(DipolesFiles{iFile});
    importedBaseName = strrep(importedBaseName, 'dipoles_', '');
    importedBaseName = strrep(importedBaseName, '_dipoles', '');
    importedBaseName = strrep(importedBaseName, 'dipoles', '');
    % Limit number of chars
    if (length(importedBaseName) > 15)
        importedBaseName = importedBaseName(1:15);
    end
    % Create output filename
    OutputFile = bst_fullfile(ProtocolInfo.STUDIES, bst_fileparts(sStudy.FileName), ['dipoles_', importedBaseName, '.mat']);
    OutputFile = file_unique(OutputFile);
    % Save new file in Brainstorm format
    bst_save(OutputFile, DipolesMat, 'v7');

    %% ===== UPDATE DATABASE =====
    % Create structure
    BstDipolesMat = db_template('Dipoles');
    BstDipolesMat.FileName = file_short(OutputFile);
    BstDipolesMat.Comment  = DipolesMat.Comment;
    BstDipolesMat.DataFile = DataFile;
    % Add to study
    iDipole = length(sStudy.Dipoles) + 1;
    sStudy.Dipoles(iDipole) = BstDipolesMat;
    iNewDipoles = [iNewDipoles, iDipole];
end

% Save study
bst_set('Study', iStudy, sStudy);
% Update tree
panel_protocols('UpdateNode', 'Study', iStudy);
% Save database
db_save();

% Progress bar
if ~isProgressBar
    bst_progress('stop');
end

author: brainstorm-tools


text:
function assert_grid(theta_start,theta_delta,theta_end,rho_start,rho_delta,rho_end,varargin)
% check for valid grid parameters

if check_option(varargin,'degree')
  stop = 180;
else
  stop = pi;
end

assert(theta_start>=0 && theta_start <= theta_end && theta_end<=stop &&...
  theta_delta > 0 && theta_delta <= theta_end - theta_start && ...
  rho_start>=0 && rho_start < rho_end && rho_end<=2*stop &&...
  rho_delta > 0 && rho_delta <= rho_end - rho_start);

author: mtex-toolbox

repo: mtex

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


text:
function plot_data(frame_array,D,X,S,n1,n2)

    T=size(frame_array,2);
    for k=1:T
        subplot('Position',[(k-1)/T 0.68 1/T 0.3]), imshow(reshape(D(:,frame_array(k)),n1,n2),[])
        subplot('Position',[(k-1)/T 0.35 1/T 0.3]), imshow(reshape(X(:,frame_array(k)),n1,n2),[])
        subplot('Position',[(k-1)/T 0.02 1/T 0.3]), imshow(reshape(S(:,frame_array(k)),n1,n2),[])
    end

% N=size(frame_array,2);
% for k=1:N
%     subplot(3,N,k), imshow(reshape(D(:,frame_array(k)),144,176),[])
%     subplot(3,N,k+N), imshow(reshape(X(:,frame_array(k)),144,176),[])
%     subplot(3,N,k+2*N), imshow(reshape(S(:,frame_array(k)),144,176),[])
% end

author: andrewssobral

repo: lrslibrary


text:
%kiconvolve 'Perform convolution or correlation on image data'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros iconvolve.pane file
%
% Parameters:
% InputFile: i1 'Input image', required: 'input image'
% InputFile: i2 'Kernel ', required: 'kernel'
% Toggle: upcast 'Upcast data', default: 0: 'Upcast data to double or double complex for processing'
% Integer: wc 'Width ', default: 0: 'Width coordinate of kernel hotspot'
% Integer: hc 'Height ', default: 0: 'Height coordinate of kernel hotspot'
% OutputFile: o 'Output image', required: 'output image'
%
% Example: o = kiconvolve({i1, i2}, {'i1','';'i2','';'upcast',0;'wc',0;'hc',0;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% iconvolve - Perform convolution or correlation on image data
%
%  DESCRIPTION
% \fBiconvolve" is a simple interface to the more general \fBklinearop\fP
% routine in the DATAMANIP toolbox. \fBiconvolve" is set up to perform
% convolution or correlation on image data (WxHx?x?x?) using a WxH kernel,
% in zero-padded, linear mode (not circular).
%
% The kernel origin, center, or "hotspot" can be set using the -wc and -hc
% arguments; otherwise the KPDS_KERNEL_ORIGIN attribute for the kernel object
% is queried for that information.
%
%
%
%  EXAMPLES
% iconvolve -i1 ball.xv -i2 laplacian -corr 0 -o edges
%
%  "SEE ALSO"
% kfft(1), kconvolve(1)
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kiconvolve(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kiconvolve(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'i2', '__input';'upcast', 0;'wc', 0;'hc', 0;'o', '__output'};
maxval={0,0,0,1,1,0};
minval={0,0,0,1,1,0};
istoggle=[0,0,1,1,1,0];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','Toggle','Integer','Integer','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kconvolve"  -doff'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [ magSpect_oct, freqVec_oct ] = octaveBandMean( magSpect, freqVec, octSpace, centerFreq )
% Given a magnitude spectrum this function will calculate the average (single, third, nth) octave band magnitude
%
% Syntax:	[ MAGSPECT_OCT, FREQVEC_OCT ] = OCTAVEBANDMEAN( MAGSPECT, FREQVEC, OCTSPACE, CENTERFREQ )
%
% Inputs:
% 	magSpect - An arbitrary magnitude spectrum as a vector
% 	freqVec - The corresponding frequencies for each magnitude value
% 	octSpace - A single value between 0 and 1 specifying the octave spacing
% 	centerFreq - The center frequency for the octave band
%
% Outputs:
% 	magSpect_oct - The magnitude spectrum averaged per octave band
% 	freqVec_oct - The corresponding frequencies for the output magnitude vector
%
% See also:

% Author: Jacob Donley
% University of Wollongong
% Email: jrd089@uowmail.edu.au
% Copyright: Jacob Donley 2016, 2017
% Date: 8 July 2016
% Revision: 0.2 (29 March 2017)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 4
    centerFreq=1e3;
end
if nargin < 3
    octSpace = 1/3;
end

if size(magSpect,2) > size(magSpect,1), magSpect = magSpect.'; end
if size(freqVec,2) > size(freqVec,1), freqVec = freqVec.'; end

freqVec_ = nonzeros(freqVec); %remove zero values

n = floor(log(freqVec_( 1 )/centerFreq)/log(2)/octSpace): ...
    floor(log(freqVec_(end)/centerFreq)/log(2)/octSpace);

freqVec_oct = centerFreq * (2 .^ (n*octSpace));
if freqVec_oct(end) ~= freqVec_(end)
    % This may happen if there is no octave band that is centered at the
    % Nyquist frequency (fs/2)
    freqVec_oct(end+1) = freqVec_(end); % So we find an average at the
    % Nyquist frequency that is of the same width as the other octave bands
    warning(['The sampling frequency used does not have a Nyquist ' ...
    'frequency that, after division by the centre frequency, is a ' ...
    'power of two. ' ...
    'i.e. (' num2str(freqVec_(end)*2) '/2)/' num2str(centerFreq) ' is ' ...
    'not a power of two.']);
end
fd = 2^(octSpace/2);
fupper = freqVec_oct * fd;
flower = freqVec_oct / fd;

octBands = [flower' fupper'];

[~, oBandInds] = min(abs( ...
    repmat(freqVec,1,2,length(octBands)) ...
    - repmat(permute(octBands,[3 2 1]),length(freqVec),1,1)  ));

oBandInds = permute(oBandInds, [3 2 1]);

magSpect_oct = zeros(1,length(freqVec_oct));
for band = 1:length(freqVec_oct)
    magSpect_oct(band+1) = mean(  ...
        magSpect( oBandInds(band,1):oBandInds(band,2) ) );
end

magSpect_oct(1) = magSpect(1);
freqVec_oct  = [0 freqVec_oct];
end

author: jtkim-kaist

repo: Speech-enhancement


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: EXAMPLE DATA ANALYSIS CODE TO LOOK AT FLOW IN EMPTY CHANNEL
%
%     Note:
%           (1) This code analyzes viz_IB2d code for channel flow in
%               /data_analysis/Example_For_Data_Analysis/Example_Flow_In_Channel/viz_IB2d
%           (2) Produces a plot of cross-sectional mag. of velocity for
%               different points along the channel at three times.
%           (3) USER-DEFINED functions are functions that users should make to
%               analyze their specific data sets
%           (4) MUST make sure to 'addpath' to where DA_Blackbox is, i.e.,
%               line 61
%           (5) MUST make sure to set path to desired dataset, i.e., in line
%               56
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Example_Channel_Flow_Analysis()

% TEMPORAL INFO FROM input2d %
dt = 1e-4;      % Time-step
Tfinal = 0.015; % Final time in simulation
pDump=50;       % Note: 'print_dump' should match from input2d

% DATA ANALYSIS INFO %
start=1;                             % 1ST interval # included in data analysis
finish=3;                            % LAST interval # included in data analysis
dump_Times = (start:1:finish)*pDump; % Time vector when data was printed in analysis

% SET PATH TO DESIRED viz_IB2d DATA and hier_IB2d DATA %
pathViz = 'viz_IB2d';
pathForce='hier_IB2d_data';

% SET PATH TO DA_BLACKBOX %
addpath('../../DA_Blackbox');

for i=start:1:finish

    % Points to desired data viz_IB2d data file
    if i<10
       numSim = ['000', num2str(i)];
    elseif i<100
       numSim = ['00', num2str(i)];
    elseif i<1000
       numSim = ['0', num2str(i)];
    else
       numSim = num2str(i);
    end

    % Imports immersed boundary positions %
    [xLag,yLag] = give_Lag_Positions(pathViz,numSim);

    % Imports (x,y) grid values and ALL EULERIAN DATA %
    %                      DEFINITIONS
    %          x: x-grid                y: y-grid
    %       Omega: vorticity           P: pressure
    %    uMag: mag. of velocity
    %    uX: mag. of x-Velocity   uY: mag. of y-Velocity
    %    U: x-directed velocity   V: y-directed velocity
    %    Fx: x-directed Force     Fy: y-directed Force
    %
    %  Note: U(j,i): j-corresponds to y-index, i to the x-index
    %
    %
    Eulerian_Flags(1) = 0;   % OMEGA
    Eulerian_Flags(2) = 0;   % PRESSURE
    Eulerian_Flags(3) = 1;   % uMAG
    Eulerian_Flags(4) = 0;   % uX (mag. x-component of velocity)
    Eulerian_Flags(5) = 0;   % uY (mag. x-component of velocity)
    Eulerian_Flags(6) = 0;   % uVEC (vector components of velocity: U,V)
    Eulerian_Flags(7) = 0;   % Fx (x-component of force )
    Eulerian_Flags(8) = 0;   % Fy (y-component of force)
    Eulerian_Flags(9) = 0;   % C (concentration)
    %
    [x,y,Omega,P,uMag,uX,uY,U,V,Fx,Fy,C] = import_Eulerian_Data(pathViz,numSim,Eulerian_Flags);

    % Imports Lagrangian Pt. FORCE (magnitude) DATA %
    %                      DEFINITIONS
    %
    %      fX_Lag: forces in x-direction on boundary
    %      fY_Lag: forces in y-direction on boundary
    %       fLagMag: magnitude of force at boundary
    %   fLagNorm: magnitude of NORMAL force at boundary
    %   fLagTan: magnitude of TANGENT force at boundary
    %
    [fX_Lag,fY_Lag,fLagMag,fLagNorm,fLagTan] = import_Lagrangian_Force_Data(pathForce,numSim);

    %
    %
    % *** USER DEFINED FUNCTIONS TO GET DESIRED ANALYSIS PT. INDICES *** %
    %                                                                    %
    if i==start
        %xPts = [0.125 0.225 0.325 0.425];
        xPts = [0.25 0.375 0.5 0.675];
        yPts = [0.405 0.595];
        [xInds,yInds] = give_Desired_Analysis_Points(x,y,xPts,yPts);
        vel_data = zeros(length(yInds),length(xInds),finish);
        %Inds = put_x_y_Indices_Together_For_Analysis(xInds,yInds);
    end

    %                                                                    %
    % ***** USER DEFINED FUNCTION TO SAVE DESIRED VELOCITY DATA *****    %
    %                                                                    %
    vel_data = store_Desired_Magnitude_Velocity_Data(uMag,vel_data,xInds,yInds,i);

end % END OF LOOP OVER ANALYZED TIME-PTS

%                                                                    %
% ***** USER DEFINED FUNCTION TO PLOT DESIRED VELOCITY DATA *****   %
%                                                                    %
yVals = y(yInds);
plot_Desired_Data(yVals,vel_data);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% USER-FUNCTION: Finds desired analysis point indices
%
%        INPUTS: x: x-grid pts (row vector)
%                y: y-grid pts (column vector)
%                xPts: desired x-pts to analyze
%                yPts: desired y-pts to analyze
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xInds,yInds] = give_Desired_Analysis_Points(x,y,xPts,yPts)

% Get x-Indices
xInds = zeros(length(xPts),1);
for i=1:length(xPts)
   xPt = xPts(i);   % Get x-Pts
   k = 1;           % Initialize loop iteration variable
   while x(k) < xPt
      xInds(i) = k;
      k = k+1;
   end
end

% Get y-Indices
yIndsAux = zeros(length(yPts),1);
for i=1:length(yPts)
   yPt = yPts(i);   % Get x-Pts
   k = 1;           % Initialize loop iteration variable
   while y(k) < yPt
      yIndsAux(i) = k;
      k = k+1;
   end
end
yInds = yIndsAux(1):1:yIndsAux(2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% USER-FUNCTION: Rearranges desired analysis indices
%
%        INPUTS: xInds: indices for x-region of interest
%                yInds: indices for y-region of interest
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Inds = put_x_y_Indices_Together_For_Analysis(xInds,yInds)

% Order new index matrix %
Inds = zeros(length(xInds)*length(yInds),2);
count = 1;
for i=1:length(xInds)
    for j=1:length(yInds)
        Inds(count,1) = xInds(i);
        Inds(count,2) = yInds(j);
        count=count+1;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% USER-FUNCTION: Stores desired magnitude of velocity data
%
%        INPUTS: uMag: magnitude of velocity from simulation
%                vel_data: stored mag. of velocity data in 3D-matrix
%                xInds/yInds: indices of where to save data
%                i: ith time-step to store data from
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function vel_data = store_Desired_Magnitude_Velocity_Data(uMag,vel_data,xInds,yInds,i)

% NOTE: vel_data(k,i): : k: magnitude of velocity data  (row)
%                      : j: xPt you're storing data for (column)
%                      : i: ith level you're storing in (level)

for j=1:length(xInds)
    for k=1:length(yInds)
       vel_data(k,j,i) = uMag( yInds(k),xInds(j) );
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% USER-FUNCTION: Plots magnitude of velocity data
%
%        INPUTS: vel_data: stored mag. of velocity data in 3D-matrix
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function plot_Desired_Data(yVals,vel_data)

fs = 16; % Font Size
lw = 5;  % Line Width
ms = 32; % MarkerSize

% Set Figure Size
FigHand = figure(1);
set(FigHand,'Position',[100,100,1024,895]);
%
% Make Figure!
%
subplot(3,1,1)
mat = vel_data(:,:,1);
maxVal = max(max(mat));
plot(yVals(2:end),vel_data(2:end,1,1),'.-','LineWidth',lw,'MarkerSize',ms); hold on;
plot(yVals(2:end),vel_data(2:end,2,1),'r.-','LineWidth',lw,'MarkerSize',ms); hold on;
plot(yVals(2:end),vel_data(2:end,3,1),'g.-','LineWidth',lw,'MarkerSize',ms); hold on;
plot(yVals(2:end),vel_data(2:end,4,1),'k.-','LineWidth',lw,'MarkerSize',ms); hold on;
axis([0.4 0.6 0 1.1*maxVal]);
leg=legend('x=0.25','x=0.375','x=0.50','x=0.625');
title('t=0.005s');
ylabel('Mag. Velocity','FontSize',fs); xlabel('y','FontSize',fs);
set(leg,'FontSize',fs);
set(gca,'FontSize',fs-1);
%
subplot(3,1,2)
mat = vel_data(:,:,2);
maxVal = max(max(mat));
plot(yVals(2:end),vel_data(2:end,1,2),'.-','LineWidth',lw,'MarkerSize',ms);  hold on;
plot(yVals(2:end),vel_data(2:end,2,2),'r.-','LineWidth',lw,'MarkerSize',ms); hold on;
plot(yVals(2:end),vel_data(2:end,3,2),'g.-','LineWidth',lw,'MarkerSize',ms); hold on;
plot(yVals(2:end),vel_data(2:end,4,2),'k.-','LineWidth',lw,'MarkerSize',ms); hold on;
axis([0.4 0.6 0 1.1*maxVal]);
leg=legend('x=0.25','x=0.375','x=0.50','x=0.625');
ylabel('Mag. Velocity','FontSize',fs); xlabel('y','FontSize',fs);
set(leg,'FontSize',fs);
set(gca,'FontSize',fs-1);
title('t=0.01s');
%
subplot(3,1,3)
mat = vel_data(:,:,3);
maxVal = max(max(mat));
plot(yVals(2:end),vel_data(2:end,1,3),'.-','LineWidth',lw,'MarkerSize',ms); hold on;
plot(yVals(2:end),vel_data(2:end,2,3),'r.-','LineWidth',lw,'MarkerSize',ms); hold on;
plot(yVals(2:end),vel_data(2:end,3,3),'g.-','LineWidth',lw,'MarkerSize',ms); hold on;
plot(yVals(2:end),vel_data(2:end,4,3),'k.-','LineWidth',lw,'MarkerSize',ms); hold on;
axis([0.4 0.6 0 1.1*maxVal]);
leg=legend('x=0.25','x=0.375','x=0.50','x=0.625');
ylabel('Mag. Velocity','FontSize',fs); xlabel('y','FontSize',fs);
set(leg,'FontSize',fs);
set(gca,'FontSize',fs-1);
title('t=0.015s');

author: nickabattista


text:
function b = isPerpendicular3d(v1, v2, varargin)
%ISPERPENDICULAR3D Check orthogonality of two 3D vectors.
%
%   B = isPerpendicular3d(V1, V2)
%   where V1 and V2 are 2 [1x3] arrays, returns 1 if the vectors are
%   orthogonal, and 0 otherwise.
%
%   Also works when V1 and V2 are two [Nx3] arrays with same number of
%   rows. In this case, return a [Nx1] array containing 1 at the positions
%   of parallel vectors.
%
%   Also works when one of V1 or V2 is scalar and the other one is [Nx3]
%   array, in this case return [Nx1] results.
%
%   B = isPerpendicular3d(V1, V2, TOL)
%   Specifies the absolute tolerance (default is 1e-14).
%
%
%   Example
%   isPerpendicular3d([1 0 0], [0 1 0])
%   ans =
%       1
%
%   isPerpendicular3d([1 0 1], [1 0 0])
%   ans =
%       0
%
%   See also
%   vectors3d, isParallel3d

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2006-04-25
% Copyright 2006-2022 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas)

% check if tolerance is specified
tol = 1e-14;
if ~isempty(varargin)
    tol = varargin{1};
end

% compute perpendicularity test
b = abs(sum(bsxfun(@times, v1, v2), 2)) < tol;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [Gx_org,Gy_org,Gz_org]=compute_gradients(img_pca,ins_msk,Subject,Figures,Streamlines)

% Compute eigenmap's gradient magnitude

% INPUT
% img_pca is the eigenmap image in which gradients are computed
% ins_msk is a mask. Only gradients for voxels in the mask are shown in
% figure
% Subject: prefix of output
% Figures: 0->suppress figures; 1->print figures
% slices
% Streamlines:1-> write out vector file; 0->do not write out vector file

% OUTPUT
% Local gradient orientations for each voxel

Erode=1; %erosion kernel
N=size(img_pca);
ind_ins=find(~~ins_msk);

% Compute gradient
% updated script to clean up edges
img_new=zeros(N);
img_new(ind_ins)=img_pca(ind_ins);
img_pca=img_new; %Remove gradients computed in boundary voxels added by dilation

[Gy,Gx,Gz,mag]=compute_grads_local(img_pca);
% %Swap order of Gx and Gy because the first output of the function is the
% %gradient along direction of increasing column subscripts.
% %With this re-ordering
% %Gx is gradient along rows    (1st dimension)
% %Gy is gradient along columns (2nd dimension)
% %Gz is gradeint in 3rd dimensions

img_mag=mag;

Gx_org=Gx; Gy_org=Gy; Gz_org=Gz;

%Crop for a better visulization
Gx=Gx.*~~ins_msk;
Gy=Gy.*~~ins_msk;
Gz=Gz.*~~ins_msk;

%Bounding box
tmp=squeeze(sum(ins_msk,1));
[y,z]=find(tmp);
max_y=max(y); min_y=min(y);
max_z=max(z); min_z=min(z);
tmp=squeeze(sum(ins_msk,2));
[x,z]=find(tmp);
max_x=max(x); min_x=min(x);
Nx=max_x-min_x+1;
Ny=max_y-min_y+1;
Nz=max_z-min_z+1;

%Write out img_pca and img_mag after cropping
if ~isempty(Subject)
    fprintf('Write out %s\n',[Subject,'magnitude.nii'])
    mat2nii(img_mag(min_x:max_x,min_y:max_y,min_z:max_z),[Subject,'magnitude.nii']);

    fprintf('Write out %s\n',[Subject,'eigenvector.nii'])
    mat2nii(img_pca(min_x:max_x,min_y:max_y,min_z:max_z),[Subject,'eigenvector.nii']);
end

if Streamlines
    u=repmat([1:Nx]',1,Ny); u=repmat(u,1,1,Nz);
    v=repmat([1:Ny],Nx,1); v=repmat(v,1,1,Nz);
    w=zeros(Nx,Ny,Nz);
    for i=1:Nz
        w(:,:,i)=ones(Nx,Ny)*i;
    end
    slice_msk=imerode(~~squeeze(ins_msk(min_x:max_x,min_y:max_y,min_z:max_z)),ones(Erode,Erode));
    gx=squeeze(Gx(min_x:max_x,min_y:max_y,min_z:max_z)).*slice_msk;
    gy=squeeze(Gy(min_x:max_x,min_y:max_y,min_z:max_z)).*slice_msk;
    gz=squeeze(Gz(min_x:max_x,min_y:max_y,min_z:max_z)).*slice_msk;

    save ([Subject,'VectorFile.mat'],'gy','gx','gz','v','u','w','Gx_org','Gy_org','Gz_org');
end

if Figures
    Mag=0; % Mag=0:print out eigenmap slices;
           % Mag=1:print out gradient magnitude slices

    % Colormap consistent with Trackvis
    mycolormap=dlmread('trackvis_jet.txt');
    mycolormap(1,:)=1; % White background

    % Rescale according to the value in the magnitude map
    img_pca(ind_ins)=img_pca(ind_ins)/0.03*63+2;
    img_mag(ind_ins)=img_mag(ind_ins)/0.03*63+2;

%     %Use default colormap
%     mycolormap=parula(256); mycolormap(1,:)=1; % White background
%     img_pca(ind_ins)=img_pca(ind_ins)/max(img_pca(ind_ins))*255+2;
%     img_mag(ind_ins)=img_mag(ind_ins)/max(img_mag(ind_ins))*255+2;
%     %Rescale to 256 colors:
%     %Bin 1->2 is white
%     %Bin 2->3 is first color
%     %Bin 256->257 is last color

    %Sagittal
    for i=1:N(1)
        sz(i)=sum(sum(~~ins_msk(i,:,:)));
    end
    [~,ind_srt]=sort(sz,'descend'); %find slices with greatest coverage

    Slices=ind_srt(1:4);
    hf=figure; hf.Position=[100,300,2000,300]; hf.Color='w';
    for i=1:length(Slices)
        subplot(1,length(Slices),i)
        if Mag
            %1st dimension of image is placed on rows
            %2nd dimension of image is placed on columns
            im=image(squeeze(img_mag(Slices(i),min_y:max_y,min_z:max_z)));
        else
            im=image(squeeze(img_pca(Slices(i),min_y:max_y,min_z:max_z)));
        end
        slice_msk=imerode(~~squeeze(ins_msk(Slices(i),min_y:max_y,min_z:max_z)),ones(Erode,Erode));
        colormap(mycolormap);
        hold on;
        gx=squeeze(Gx(Slices(i),min_y:max_y,min_z:max_z)).*slice_msk;
        gy=squeeze(Gy(Slices(i),min_y:max_y,min_z:max_z)).*slice_msk;
        gz=squeeze(Gz(Slices(i),min_y:max_y,min_z:max_z)).*slice_msk;

        % Even the magnitude for figure
        mag=sqrt(gy.^2+gx.^2+gz.^2);
        gy=gy./mag; gx=gx./mag; gz=gz./mag;

        u=repmat([1:Ny]',1,Nz); v=repmat([1:Nz],Ny,1); w=ones(Ny,Nz);
        q=quiver3(v,u,w,gz,gy,gx);

        %First input is gradient in direction of increasing columns
        %Second input is gradient in direction of increasing rows
        q.LineWidth=0.5;
        q.Color=[192,192,192]/255;
        q.AutoScaleFactor=2;
        view(-90, 90);
        axis off;
        ax=gca; ax.XLabel.Visible='on'; ax.YLabel.Visible='on';
        axis equal;
    end

    % Coronal
    for i=1:N(2)
        sz(i)=sum(sum(~~ins_msk(:,i,:)));
    end
    [~,ind_srt]=sort(sz,'descend'); %find slices with greatest coverage
    Slices=ind_srt(1:4);
    hf=figure; hf.Position=[100,800,2000,300]; hf.Color='w';
    for i=1:length(Slices)
        subplot(1,length(Slices),i)
        if Mag
            %1st dimension of image is placed on rows
            %2nd dimension of image is placed on columns
            tmp=squeeze(img_mag(min_x:max_x,Slices(i),min_z:max_z));
            im=image(tmp);
        else
            tmp=squeeze(img_pca(min_x:max_x,Slices(i),min_z:max_z));
            im=image(tmp);
        end
        slice_msk=imerode(~~squeeze(ins_msk(min_x:max_x,Slices(i),min_z:max_z)),ones(Erode,Erode));
        colormap(mycolormap);
        hold on;
        gx=squeeze(Gx(min_x:max_x,Slices(i),min_z:max_z)).*slice_msk;
        gy=squeeze(Gy(min_x:max_x,Slices(i),min_z:max_z)).*slice_msk;
        gz=squeeze(Gz(min_x:max_x,Slices(i),min_z:max_z)).*slice_msk;

        % Even the magnitude
        mag=sqrt(gy.^2+gx.^2+gz.^2);
        gy=gy./mag; gx=gx./mag; gz=gz./mag;

        u=repmat([1:Nx]',1,Nz); v=repmat([1:Nz],Nx,1); w=ones(Nx,Nz);
        q=quiver3(v,u,w,gz,gx,gy);

        %First input is gradient in direction of increasing columns
        %Second input is gradient in direction of increasing rows
        q.LineWidth=0.5;
        q.Color=[192,192,192]/255;
        q.AutoScaleFactor=2;
        view(-90, 90);
        axis off;
        ax=gca; ax.XLabel.Visible='on'; ax.YLabel.Visible='on';
        axis equal;
    end

    %Axial
    for i=1:N(3)
        sz(i)=sum(sum(~~ins_msk(:,:,i)));
    end
    [~,ind_srt]=sort(sz,'descend'); %find slices with greatest coverage
    Slices=ind_srt(1:4);
    hf=figure; hf.Position=[100,1200,2000,300]; hf.Color='w';
    for i=1:length(Slices)
        subplot(1,length(Slices),i)
        if Mag
            %1st dimension of image is placed on rows
            %2nd dimension of image is placed on columns
            tmp=squeeze(img_mag(min_x:max_x,min_y:max_y,Slices(i)));
            im=image(tmp);
        else
            tmp=squeeze(img_pca(min_x:max_x,min_y:max_y,Slices(i)));
            im=image(tmp);
        end
        slice_msk=imerode(~~squeeze(ins_msk(min_x:max_x,min_y:max_y,Slices(i))),ones(Erode,Erode));
        colormap(mycolormap);
        hold on;
        gx=squeeze(Gx(min_x:max_x,min_y:max_y,Slices(i))).*slice_msk;
        gy=squeeze(Gy(min_x:max_x,min_y:max_y,Slices(i))).*slice_msk;
        gz=squeeze(Gz(min_x:max_x,min_y:max_y,Slices(i))).*slice_msk;

        % Even the magnitude
        mag=sqrt(gy.^2+gx.^2+gz.^2);
        gy=gy./mag; gx=gx./mag; gz=gz./mag;

        u=repmat([1:Nx]',1,Ny); v=repmat([1:Ny],Nx,1); w=ones(Nx,Ny);

        q=quiver3(v,u,w,gy,gx,gz);

        %First input is gradient in direction of increasing columns
        %Second input is gradient in direction of increasing rows
        q.LineWidth=0.5;
        q.Color=[192,192,192]/255;
        q.AutoScaleFactor=2;

        view(-90, 90);
        axis off;
        ax=gca; ax.XLabel.Visible='on'; ax.YLabel.Visible='on';
        axis equal;
    end
end

author: yetianmed

repo: subcortex

save_path: github-repos/MATLAB/yetianmed-subcortex


text:
function B = adjugate(A)

% This finds the adjugate (adjoint) of square matrix A,
% and is valid even if A is singular or complex-valued.
% With u, s, and v obtained from [u,s,v] = svd(A), it
% makes use of the identity adj(A) = det(u*v')*v*adj(s)*u',
% which holds even if A and s are singular.  The expression,
% diag(prod(reshape(s0(ix),n-1,n),1)), accomplishes the
% evaluation of adj(s), each of whose diagonal elements
% is the product of all but one of the diagonal elements
% of s.  This requires that A be n x n where n >= 2.
% Roger Stafford - 10/18/06

[m,n] = size(A);
if (m ~= n) | (n < 2)
 error('Matrix A should be size n x n with n >= 2.')
end
[u,s,v] = svd(A);
s0 = diag(s);
ix = toeplitz(ones(n-1,1),[1 zeros(1,n-1)]) ...
     + repmat((1:n-1)',1,n);
B = det(u*v')*v*diag(prod(reshape(s0(ix),n-1,n),1))*u';

% ---

author: Sable

repo: mcbench-benchmarks


text:
function [V,F] = schwarz_lantern(m,n)
  % SCHWARZ_LANTERN Construct a Schwarz lattern with m axial slices and n radial
  % vertices.
  %
  % Inputs:
  %   m  number of axial "slices"
  %   n  number of vertices around each slice
  % Outputs:
  %   V  mn by 3 list of vertex positions
  %   F  #F by 3 list of triangle indices into V
  %
  [X,Y,Z] = cylinder(ones(m,1),n);
  %[F,V] = surf2patch(X,Y,Z,'triangles');
  [Q,V] = surf2patch(X,Y,Z);
  R = axisangle2matrix([0 0 1],-2*pi/n*0.5);
  other = ((1:n+1)-1)*m+((2:2:m-mod(m,2)+1)');
  V(other,:) = V(other,:)*R;
  other = ((1:n)-1)*(m-1)+((1:2:m-1)');
  %tsurf(Q,V,'Tets',false,'FaceIndices',1,'CData',sparse(other,1,1,size(Q,1),1));
  %view(90,0);
  other = ismember(1:size(Q,1),other);
  F = [Q(~other,[1 2 3]);Q(~other,[1 3 4]);Q(other,[1 2 4]);Q(other,[4 2 3])];
  %tsurf(F,V);
  %view(90,0);
  % Heavy-handed way of removing duplicates...
  [V,~,IM] = remove_duplicate_vertices(V,eps);
  F = IM(F);
end

author: alecjacobson

repo: gptoolbox


text:
function test_bug2186

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

data = [];
for i=1:10
  data.label{i} = num2str(i);
end
for i=1:20
  data.time{i} = (1:1000)/1000;
  data.trial{i} = randn(10, 1000);
end

cfg = [];
cfg.trials = 1:10;
data1 = ft_selectdata(cfg, data);
assert(data1.trial{1}(1)==data.trial{1}(1));

cfg = [];
cfg.trials = 11:20;
data2 = ft_selectdata(cfg, data);
assert(data2.trial{1}(1)==data.trial{11}(1));

%%

cfg = [];
cfg.keeptrials = 'no';
cfg.covariance = 'yes';
timelock1 = ft_timelockanalysis(cfg, data1);
timelock2 = ft_timelockanalysis(cfg, data2);

append12 = ft_appendtimelock([], timelock1, timelock2);
assert(isfield(append12, 'trial')); % renamed from avg
assert(isfield(append12, 'dof'));
assert(isfield(append12, 'var'));
assert(isfield(append12, 'cov')); % this one is new

%%

timelock1a = timelock1;
for i=1:10
  timelock1a.label{i} = num2str(i+10);
end
append1a = ft_appendtimelock([], timelock1, timelock1a);
assert(isequal(size(append1a.cov), [20 20]));
assert(isnan(append1a.cov(20, 1))); % it should be a block diagonal matrix with nan off-diagonal

%%

cfg = [];
cfg.keeptrials = 'yes';
cfg.covariance = 'yes';
timelock1 = ft_timelockanalysis(cfg, data1);
timelock2 = ft_timelockanalysis(cfg, data2);

append12 = ft_appendtimelock([], timelock1, timelock2);
assert(isfield(append12, 'trial')); % avg, dof and var are not present
assert(isfield(append12, 'cov')); % this one is new

%%

timelock1a = timelock1;
for i=1:10
  timelock1a.label{i} = num2str(i+10);
end
append1a = ft_appendtimelock([], timelock1, timelock1a);
assert(isequal(size(append1a.cov), [10 20 20]));
assert(isnan(append1a.cov(1, 20, 1))); % it should be a block diagonal matrix with nan off-diagonal

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% +MSG
%
% Files
%   dbdisp    -
%   dbfprintf -
%   errordisp -
%   infodisp  -
%   warndisp  - display a warning message to command line

author: CelsoReyes


text:
function bgdiff = bgtdiff(ll,llex)
    %bgtdiff.m                                         A.Allmann
    %calculates the time difference between jth and biggest event
    %in an aftershock sequence
    %works with eqtime from function clustime.m
    %gives the indices ac of the eqs not already related to cluster k1
    %Last modification 6/95
    global  eqtime k1 bg mbg

    z=ll(llex);                               %indices of the eqs to examine
    %mbg(k1)=bg(k1);                   %bgdiff not calculated from the last
    %  biggest eq but the first biggest eq
    bgdiff=zeros(length(z),1);
    if bg(k1)~=mbg(k1)                          %if more eqs with biggest magnitude
        tm1 = max(find(z<=mbg(k1)));                 %position of mbg

        if size(tm1)~=0                          %if eqs before mbg(k1)
            bgdiff(1:tm1,1)=eqtime(z(1:tm1),1)-eqtime(bg(k1));

            if tm1~=length(z)                    %if mbg(k1) is not the last eq to examine
                bgdiff((tm1+1):length(z),1)=eqtime(z((tm1+1):length(z)),1)-eqtime(mbg(k1));
            end

        else                          %bg and mbg are smaller than all indices in z
            bgdiff=eqtime(z)-eqtime(mbg(k1));
        end

    else                              %only one eq with biggest magnitude
        bgdiff=eqtime(z)-eqtime(bg(k1));

    end


text:
function [dimensions,x1,x2,x3] = xyz_axes_selection(x,y,z)
%XYZ_AXES_SELECTION returns the first non-singleton axes and a vector
%indicating which axes are selected
%
%   Usage: [dimensions,x1,x2,x3] = xyz_axes_selection(x,y,z)
%
%   Input parameters:
%       x,y,z      - vectors/matrices containing the x-, y- and z-axis values / m
%
%   Output parameters:
%       dimensions - 1x3 vector containing 1 or 0 to indicate the activity
%                    of the single dimensions in the order [x y z]
%       x1         - vector/matrix containing the first axis / m
%       x2         - vector/matrix containing the second axis / m
%       x3         - vector/matrix containing the third axis / m
%
%   XYZ_AXES_SELECTION(x,y,z) returns an indication vector for the x-, y- and
%   z-axis if we have any activity on this axis or if it is a singleton axis.
%   In addition, the axes are reordered starting first with the non-singleton
%   axes.
%
%   See also: plot_sound_field, xyz_grid, is_dim_custom

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input parameters ====================================
nargmin = 3;
nargmax = 3;
narginchk(nargmin,nargmax);
isargnumeric(x,y,z);

%% ===== Computation =====================================================
dims = {x,y,z};
dimensions = ~is_dim_singleton(dims{:});
Nd = sum(dimensions);

newdims = {x(1),y(1),z(1)};  % default case, if all dimensions are singleton
newdims(Nd+1:end) = newdims(~dimensions); % move singleton dimensions to the end
newdims(1:Nd) = dims(dimensions);  % move non-singleton dimensions to the front

[x1, x2, x3] = newdims{:};

author: sfstoolbox

repo: sfs-matlab


text:
% Tow Thomas Biquad Bandpass Filter - Leverriers Algorithm
% File:  TowThomasLev.m
% 9/22/02
clc;clear;format short g;
K=1e3;uF=1e-6;
R1=200*K;R2=10*K;R3=20*K;R4=10*K;R5=20*K;R6=20*K;
C1=0.0159*uF;C2=C1;
X=[R1 R2 R3 R4 R5 R6 C1 C2]; % Put components in vector form
%
[A,B,D,E]=tow(X);
%
% Uncomment following two lines to output arrays to text file qbout.txt
%fid=fopen('c:\M_files\qbout.txt','w'); % Use local directory
%diary c:\M_Files\qbout.txt; % Use local directory
A
B
D
E
%
% Display results of Leverrier's Algorithm
%
F1=eye(2);T1=-trace(A*F1)/1
F0=A*F1+T1*F1
T0=-trace(A*F0)/2
Y1=D*F1*B+E*T1
Y0=D*F0*B+E*T0
%
% Uncomment the following two lines to close text file.
%diary off
%status=fclose(fid);
%
% Get frequency sweep; BF = Beginning (Log) Frequency, ND = Number of Decades
% PD = Points per Decade
BF=2;ND=2;PD=50;Lit=ND*PD+1;L=linspace(BF,BF+ND,Lit);
for i=1:Lit
   F=10^L(i);s=2*pi*F*j;
   for k=1:3 % Get all three transfer functions from transfer matrix
      num=E(k)*s^2+Y1(k)*s+Y0(k);
      den=s^2+T1*s+T0;
      Vo(k,i)=20*log10(abs(num/den));
   end
end
%
% Plot Vo
%
h=plot(L,Vo(1,:),'k',L,Vo(2,:),'r',L,Vo(3,:),'b');
set(h,'LineWidth',2);
grid on
xlabel('Log Freq(Hz)');
ylabel('dBV');
title('Transfer Matrix Outputs');
legend('V2','V4','V6');
figure(1);


text:
classdef nndagnn < nntest

  properties
    net;
    x;
    class;
    inputs;
    outputs;
  end

  methods (TestClassSetup)
    function initNet(test, device)
      test.net = [];
      test.net.layers = {} ;
      test.net.layers{end+1} = struct('type', 'conv', ...
        'weights', {{test.randn(5,5,1,20), test.zeros(1, 20)}}, ...
        'stride', 1, ...
        'pad', 0) ;
      test.net.layers{end+1} = struct('type', 'pool', ...
        'method', 'max', ...
        'pool', [2 2], ...
        'stride', 2, ...
        'pad', 0) ;
      test.net.layers{end+1} = struct('type', 'conv', ...
        'weights', {{test.randn(5,5,20,50),test.zeros(1,50)}}, ...
        'stride', 1, ...
        'pad', 0) ;
      test.net.layers{end+1} = struct('type', 'pool', ...
        'method', 'max', ...
        'pool', [2 2], ...
        'stride', 2, ...
        'pad', 0) ;
      test.net.layers{end+1} = struct('type', 'conv', ...
        'weights', {{test.randn(4,4,50,500),test.zeros(1,500)}}, ...
        'stride', 1, ...
        'pad', 0) ;
      test.net.layers{end+1} = struct('type', 'relu') ;
      test.net.layers{end+1} = struct('type', 'conv', ...
        'weights', {{test.randn(1,1,500,10),test.zeros(1,10)}}, ...
        'stride', 1, ...
        'pad', 0) ;
      test.net.layers{end+1} = struct('type', 'softmaxloss') ;
      % Fill the missing values
      test.net = dagnn.DagNN.fromSimpleNN(vl_simplenn_tidy(test.net));
      test.inputs = test.net.getInputs();
      test.outputs = test.net.getOutputs();

      test.x = test.randn(32, 32, 1, 20) ;
      test.class = test.toDevice(randi(10, 20, 1));
      test.net.move(device) ;
    end
  end

  methods (Test)
    function simpleRun(test)
      % Verify the forget functionality for the forward pass
      test.net.conserveMemory = true;
      test.forward();
      for ri = 1:numel(test.net.vars)
        if ismember(test.net.vars(ri).name, test.outputs)
          test.verifyNotEmpty(test.net.vars(ri).value);
        else
          test.verifyEmpty(test.net.vars(ri).value);
        end
      end
      % Verify the forget functionality for the backward pass
      test.backward();
      for ri = 1:numel(test.net.vars)
        if ismember(test.net.vars(ri).name, test.inputs)
          test.verifyNotEmpty(test.net.vars(ri).value);
        else
          test.verifyEmpty(test.net.vars(ri).value);
        end
      end
      for pi = 1:numel(test.net.params)
        test.verifyNotEmpty(test.net.params(pi).der);
      end
    end

    function conserveMemory(test)
      % Verify the forget functionality for the forward pass
      test.net.conserveMemory = false;
      test.forward();
      for ri = 1:numel(test.net.vars)
          test.verifyNotEmpty(test.net.vars(ri).value);
      end
      % Verify the forget functionality for the backward pass
      test.backward();
      for ri = 1:numel(test.net.vars)
        test.verifyNotEmpty(test.net.vars(ri).value);
        if ~ismember(test.net.vars(ri).name, test.inputs)
          test.verifyNotEmpty(test.net.vars(ri).der);
        end
      end
    end

    function precious(test)
      % Verify that the precious argument works
      selLayer = 3;
      outputIdx = test.net.layers(selLayer).outputIndexes;
      test.net.vars(outputIdx).precious = true;
      test.net.conserveMemory = true;
      test.forward();
      for ri = 1:numel(test.net.vars)
        if ismember(test.net.vars(ri).name, test.outputs) || ri == outputIdx
          test.verifyNotEmpty(test.net.vars(ri).value);
        else
          test.verifyEmpty(test.net.vars(ri).value);
        end
      end
      test.net.vars(outputIdx).precious = false;
    end
  end

  methods
    function forward(test)
      test.net.reset();
      test.net.eval({'x0', test.x, 'label', test.class});
    end
    function backward(test)
      test.net.reset();
      test.net.eval({'x0', test.x, 'label', test.class}, {'x8', 1});
    end
  end
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function varargout = clipPolygonHP(poly, line, varargin)
%CLIPPOLYGONHP Clip a polygon with a half-plane defined by a directed line.
%
%   POLY2 = clipPolygonHP(POLY, LINE)
%   POLY is a [Nx2] array of points, and LINE is given as [x0 y0 dx dy].
%   The result POLY2 is also an array of points, sometimes smaller than
%   poly, and that can be [0x2] (empty polygon). POLY2 contains the part of
%   POLY on the left side of the directed line.
%
%   [POLY_L, POLY_R] = clipPolygonHP(POLY, LINE, 'method', 'polyshape')
%   Uses MATLAB polyshape objects and functions to clip the polygon by the
%   line. Returns the right part POLY_R in addition to the left part
%   POLY_L in the polygon cell format.
%
%   Example
%     line = [0.4 0 1 1];
%     r = [2.5, 2, 1];
%     poly = flipud(circleToPolygon([0 0 r(1)], round(2*pi*r(1))));
%     poly2 = clipPolygonHP(poly, line);
%     figure('color','w','numbertitle','off','name','Method: legland')
%     axis equal tight; hold on; xlabel('x'); ylabel('y')
%     fillPolygon(poly2)
%     poly2_centroid = polygonCentroid(poly2);
%     drawLabels(poly2_centroid,'L ','HorizontalAlignment','Right')
%     scatter(poly2_centroid(1), poly2_centroid(2),[],'k','filled')
%     midCircle = circleToPolygon([0 0 r(2)], round(2*pi*r(2)));
%     innerCircle = flipud(circleToPolygon([0 0 r(3)], round(2*pi*r(3))));
%     poly = {poly, midCircle, innerCircle};
%     clipPolygonHP(poly, line, 'method','polyshape','debug',1);
%
%   See also
%   clipPolygon

% ------
% Author: David Legland, oqilipo
% E-mail: david.legland@inrae.fr
% Created: 2005-07-31
% Copyright 2005-2022 INRA - Cepia Software Platform

% Parsing
p = inputParser;
logParValidFunc=@(x) (islogical(x) || isequal(x,1) || isequal(x,0));
addParameter(p,'method','legland',@(x) any(validatestring(x,{'legland','polyshape'})));
addParameter(p,'debugVisu',false,logParValidFunc);
parse(p,varargin{:});

method = p.Results.method;
debugVisu = p.Results.debugVisu;

if iscell(poly) && numel(poly)==1
    poly = poly{1};
elseif iscell(poly) && length(poly) > 1
    method = 'polyshape';
end

if ~iscell(poly) && any(isnan(poly(:)))
    poly = splitPolygons(poly);
    method = 'polyshape';
end

switch method
    case 'legland'
        % avoid to process empty polygons
        if size(poly, 1)<3
            varargout{1} = zeros([0 2]);
            return;
        end

        % ensure the last point is the same as the first one
        if sum(poly(end, :)==poly(1,:))~=2
            poly = [poly; poly(1,:)];
        end

        N = size(poly, 1);
        edges = [poly([N 1:N-1], :) poly];

        b = isLeftOriented(poly, line);

        % case of totally clipped polygon
        if sum(b)==0
            varargout{1} = zeros(0, 2);
            return;
        end

        poly2 = zeros(0, 2);

        i=1;
        while i<=N
            if isLeftOriented(poly(i,:), line)
                % keep all points located on the right side of line
                poly2 = [poly2; poly(i,:)]; %#ok<AGROW>
            else
                % compute of preceeding edge with line
                if i>1
                    poly2 = [poly2; intersectLineEdge(line, edges(i, :))]; %#ok<AGROW>
                end
                % go to the next point on the left side
                i=i+1;
                while i<=N

                    % find the next point on the right side
                    if isLeftOriented(poly(i,:), line)
                        % add intersection of previous edge
                        poly2 = [poly2; intersectLineEdge(line, edges(i, :))]; %#ok<AGROW>

                        % add current point
                        poly2 = [poly2; poly(i,:)]; %#ok<AGROW>

                        % exit the second loop
                        break;
                    end
                    i=i+1;
                end
            end
            i=i+1;
        end

        % remove last point if it is the same as the first one
        if sum(poly2(end, :)==poly(1,:))==2
            poly2 = poly2(1:end-1, :);
        end

        if debugVisu
            figure('color','w','numbertitle','off', ...
                'name', ['Debug Figure: ' mfilename ...
                '.m, Method: ' method]);
            axis equal tight; hold on; xlabel('x'); ylabel('y')
            fillPolygon(poly2,'g')
            poly2_centroid = polygonCentroid(poly2);
            drawLabels(poly2_centroid,'L ','HorizontalAlignment','Right')
            scatter(poly2_centroid(1), poly2_centroid(2),[],'k','filled')
        end

        varargout{1} = poly2;

    case 'polyshape'
        warning('off','MATLAB:polyshape:repairedBySimplify')
        polyShape = polygonToPolyshape(poly, 'debugVisu',0);
        warning('on','MATLAB:polyshape:repairedBySimplify')
        [bblim(1:2), bblim(3:4)] = boundingbox(polyShape);
        bbXdist = bblim(2)-bblim(1);
        bbYdist = bblim(4)-bblim(3);
        % Increase the bounding box a little bit
        bblim(1) = bblim(1)-0.01*bbXdist;
        bblim(2) = bblim(2)+0.01*bbXdist;
        bblim(3) = bblim(3)-0.01*bbYdist;
        bblim(4) = bblim(4)+0.01*bbYdist;
        % Bounding box
        BB = [bblim(1) bblim(3); bblim(2) bblim(3); ...
            bblim(2) bblim(4); bblim(1) bblim(4)];
        % Clip the bounding box by the line
        BB_L = clipPolygonHP(BB, line);
        PS_R = subtract(polyShape, polyshape(BB_L));
        lineRev = [line(1:2) -line(3:4)];
        BB_R = clipPolygonHP(BB, lineRev);
        PS_L = subtract(polyShape, polyshape(BB_R));
        lineSeg = clipLine(line, bblim);
        lineSeg = [lineSeg(1:2); lineSeg(3:4)];
        % Intersection edges
        itsEdges = intersect(polyShape,lineSeg);
        itsEdges(2:end+1, 3:4) = itsEdges;
        itsEdges(1,3:4) = nan;
        itsEdges(end,1:2) = nan;
        itsEdges(any(isnan(itsEdges),2),:) = [];

        if debugVisu
            figure('color','w','numbertitle','off', ...
                'name', ['Debug Figure: ' mfilename ...
                '.m, Method: ' method]);
            axis equal tight; hold on; xlabel('x'); ylabel('y')
            plot(PS_R,'FaceColor','g','EdgeColor','k')
            plot(PS_L,'FaceColor','b','EdgeColor','k')
            [PS_R_centroid(1), PS_R_centroid(2)] = centroid(PS_R);
            drawLabels(PS_R_centroid,' R','HorizontalAlignment','Left')
            scatter(PS_R_centroid(1), PS_R_centroid(2),[],'g','filled')
            [PS_L_centroid(1), PS_L_centroid(2)] = centroid(PS_L);
            drawLabels(PS_L_centroid,'L ','HorizontalAlignment','Right')
            scatter(PS_L_centroid(1), PS_L_centroid(2),[],'b','filled')
            drawEdge(itsEdges,'LineStyle','-','LineWidth',2,'Color','k')
        end

        varargout{1} = splitPolygons(PS_L.Vertices);
        varargout{2} = splitPolygons(PS_R.Vertices);
        varargout{3} = itsEdges;

end

end


text:
function [universal_path, is_valid] = checkPath(path)

% SYNTAX:
%   universal_path = checkPath(path)
%
% INPUT:
%   path
%
% OUTPUT:
%   universal_path
%   < is_valid >            optional, contains the status of existence
%
% DESCRIPTION:
%   Conversion of path OS specific to a universal one: "/" or "\" are converted to filesep
%   if the second output parameter is present is_valid contains the status of existence
%       2 => is a file
%       7 => is a folder

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Andrea Gatti
%  Contributors:     Andrea Gatti
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

if nargout == 2
    [universal_path, is_valid] = File_Name_Processor.checkPath(path);
else
    universal_path = File_Name_Processor.checkPath(path);
end

author: goGPS-Project


text:
function p = getcustompwd

% GETCUSTOMPWD determines the present working directory, ensuring that
% it is not the private directory of the distributed computing toolbox.

% Copyright (C) 2011-2012, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are for faster processing on subsequent calls
persistent previous_pwd previous_argout

t = pwd;

if isequal(t, previous_pwd)
  % don't do the processing again, but return the previous values from cache
  p = previous_argout;
  return
end

% don't use the present directory if it contains the peer code
% it will confuse the worker with a potentially different mex file
if strcmp(pwd, fileparts(mfilename('fullpath')))
  warning('will not change directory to %s', t);
  p = [];
else
  p = t;
end

% remember the current input and output arguments, so that they can be
% reused on a subsequent call in case the same input argument is given
previous_pwd    = t;
previous_argout = p;


function Y = vl_nnsoftmaxloss(X,c,dzdy)
% VL_NNSOFTMAXLOSS  CNN combined softmax and logistic loss
%    Y = VL_NNSOFTMAX(X, C) applies the softmax operator followed by
%    the logistic loss the data X. X has dimension H x W x D x N,
%    packing N arrays of W x H D-dimensional vectors.
%
%    C contains the class labels, which should be integer in the range
%    1 to D.  C can be an array with either N elements or with H x W x
%    1 x N dimensions. In the fist case, a given class label is
%    applied at all spatial locations; in the second case, different
%    class labels can be specified for different locations.
%
%    D can be thought of as the number of possible classes and the
%    function computes the softmax along the D dimension. Often W=H=1,
%    but this is not a requirement, as the operator is applied
%    convolutionally at all spatial locations.
%
%    DZDX = VL_NNSOFTMAXLOSS(X, C, DZDY) computes the derivative DZDX
%    of the CNN with respect to the input X given the derivative DZDY
%    with respect to the block output Y. DZDX has the same dimension
%    as X.

% Copyright (C) 2014 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

%X = X + 1e-6 ;
sz = [size(X,1) size(X,2) size(X,3) size(X,4)] ;

% index from 0
c = c - 1 ;

if numel(c) == sz(4)
  % one label per image
  c = reshape(c, [1 1 1 sz(4)]) ;
  c = repmat(c, [sz(1) sz(2)]) ;
else
  % one label per spatial location
  sz_ = [size(c,1) size(c,2) size(c,3) size(c,4)] ;
  assert(isequal(sz_, [sz(1) sz(2) 1 sz(4)])) ;
end

% convert to indeces
c_ = 0:numel(c)-1 ;
c_ = 1 + ...
  mod(c_, sz(1)*sz(2)) + ...
  (sz(1)*sz(2)) * c(:)' + ...
  (sz(1)*sz(2)*sz(3)) * floor(c_/(sz(1)*sz(2))) ;

% compute softmaxloss
Xmax = max(X,[],3) ;
ex = exp(bsxfun(@minus, X, Xmax)) ;

n = sz(1)*sz(2) ;
if nargin <= 2
  t = Xmax + log(sum(ex,3)) - reshape(X(c_), [sz(1:2) 1 sz(4)]) ;
  Y = sum(t(:)) / n ;
else
  Y = bsxfun(@rdivide, ex, sum(ex,3)) ;
  Y(c_) = Y(c_) - 1;
  Y = Y * (dzdy / n) ;
end


text:
function determ = creation_determinant ( n )

%*****************************************************************************80
%
%% CREATION_DETERMINANT returns the determinant of the CREATION matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real DETERM, the determinant.
%
  determ = 0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [s] = loadLTPLLogs(Ts, Ts_slow)

% Author: Thomas Herrmann     Date: 19-01-2021
%
% Description: This function loads the logged LTPL data, which arrived at
% the control module and prepares these LTPL information for a log data
% replay in Simulink. Using this funcitonality, the controller can be
% developed further.
%
% Input:    Ts - logging file
%           Ts_slow - slow logging file
%
% Output: s - struct with data from timeseries input oject

% specify signals you want to extract, 'd' must be a fixed string
sig_names = struct(...
    'd1', 'debug_slow_tartraj_x_m', ...
    'd2', 'debug_slow_tartraj_y_m', ...
    'd3', 'debug_slow_tartraj_psi_rad', ...
    'd4', 'debug_slow_tartraj_kappa_radpm', ...
    'd5', 'debug_slow_tartraj_v_mps', ...
    'd6', 'debug_slow_tartraj_ax_mps2', ...
    'd7', 'debug_slow_tartraj_s_loc_m', ...
    'd8', 'debug_slow_tartraj_LapCnt', ...
    'd9', 'debug_slow_tartraj_ax_lim_mps2', ...
    'd10', 'debug_slow_tartraj_ay_lim_mps2', ...
    'd11', 'debug_slow_tartraj_banking_rad');

% specify integer of how much faster the fast logging file was recorded
slow_fast_factor = 10;

%% position of vehicle (fast logging file)

% get time of first entry in slow debug file
t_zero = Ts_slow.debug_slow_ax_traj_mps2.Time(1);
% find synchronisation point in fast debug file
t_zero_fast_idx = find(t_zero==Ts.debug_mvdc_path_matching_debug_ActualTrajPoint_s_glob_m.Time);

% get global s coordinate from fast logging file which matches every entry
% in the slow logging file
s_global_slow_match_m = Ts.debug_mvdc_path_matching_debug_ActualTrajPoint_s_glob_m.Data(t_zero_fast_idx:slow_fast_factor:end);

%%
% initialize empty struct
s = struct();

s.s_global_m = s_global_slow_match_m;

% remove all entries with only zeros at beginning of race
zero_counter = 1;
for k=1:length(Ts_slow.(sig_names.d1).Data)
    if sum(abs(Ts_slow.(sig_names.d1).Data(k, :))) == 0
        zero_counter = zero_counter + 1;
    end
end

s.s_global_m = s.s_global_m(zero_counter:end);
for k=1:length( fields(sig_names) )
    en = ['d', num2str(k)];
    s.(sig_names.(en)) = Ts_slow.(sig_names.(en)).Data(zero_counter:end, :);
end

%% Write to data dictionary
DDObj = ...
    Simulink.data.dictionary.open('TrajectoryPlanningEmulation.sldd');
dataSectObj = getSection(DDObj, 'Design Data');
try addEntry(dataSectObj, 'ltpl_log', s);
% if entry was there already modify it
catch e
    if isa(e, 'MSLException')
        ltplLogObj = getEntry(dataSectObj, 'ltpl_log');
        setValue(ltplLogObj, s);
    else
        print('Something went wrong.')
    end
end

% set switch to enable local log replay
switchObj = getEntry(dataSectObj, 'P_VDC_TrajEmulation_mode');
setValue(switchObj, 1);
% save & close
saveChanges(DDObj);
close(DDObj);

% set start position
Scenario_DD = Simulink.data.dictionary.open('raceline.sldd');
dataSectObj = getSection(Scenario_DD, 'Design Data');
x0_pose = getEntry(dataSectObj, 'x0_vehiclepose_stm');
setValue(x0_pose, ...
    [s.debug_slow_tartraj_x_m(1, 1), ...
     s.debug_slow_tartraj_y_m(1, 1), ...
     s.debug_slow_tartraj_psi_rad(1, 1)]);
saveChanges(Scenario_DD);
close(Scenario_DD);

end

author: TUMFTM


text:
function cut = disjunctivecut(varargin)
%
%
% x = sdpvar(2,1);
% A1 = randn(8,2);
% b1 = rand(8,1)*2-A1*[3;3];
% A2 = randn(8,2);
% b2 = rand(8,1)*2-A2*[-3;3];
% A3 = randn(8,2);
% b3 = rand(8,1)*2-A3*[3;-3];
% C1 = A1*x < b1;
% C2 = A2*x < b2;
% C3 = A3*x < b3;
% cut = disjunctivecut(C1,C2,C3,[-8;0]);cut = disjunctivecut(C1,C2,C3,[-8;0]);
% plot([cut,-10<x<10]);plot(hull(C1,C2,C3),x);plot(C1,[],'y');plot(C2,[],'y');plot(C3,[],'y')
xstar = varargin{end};
x = [];
beta = sdpvar(1);
alpha = sdpvar(length(xstar),1);
for i = 1:nargin-1
    [Imodel,Iax1,Iax2,p{i}] = export(varargin{i},[],sdpsettings,[],[],0);
    neq(i) = p{i}.K.f;
    b{i} = -p{i}.F_struc(:,1);
    A{i} = p{i}.F_struc(:,2:end);
    mu{i} = sdpvar(length(b{i}),1);
end

Objective  = alpha'*xstar-beta;
Constraints = [-1<alpha<1];

summu = 0;
for i = 1:nargin-1
    summu = summu + sum(mu{i});
    Constraints = [Constraints,alpha' == mu{i}'*A{i}];
    Constraints = [Constraints,beta  <= mu{i}'*b{i}];
    Constraints = [Constraints,mu{i}(p{i}.K.f+1:end)>0];
end

solvesdp(Constraints,Objective,sdpsettings('verbose',0));

x = recover(p{1}.used_variables);
cut = (-double(beta)+double(alpha)'*x >= 0);

return

%x = recover(sdpvar(C1));

%[Imodel,Iax1,Iax2,p1] = export(C1,[],sdpsettings,[],[],0);
%[Omodel,Oax1,Oax2,p2] = export(C2,[],sdpsettings,[],[],0);

%neq1 = p1.K.f;
%neq2 = p2.K.f;

%e1 = p1.F_struc*[1;x];
%e2 = p2.F_struc*[1;x];
%Model1 = [e2(1+p1.K.f:end)>=0];
%Model2 = [e2(1+p2.K.f:end)>=0];
if 0
Ab1 = getbase(sdpvar(Model1));
Ab2 = getbase(sdpvar(Model2));
b1 = -Ab1(:,1);
A1 =  Ab1(:,2:end);
b2 = -Ab2(:,1);
A2 = Ab2(:,2:end);

b1 = -p1.F_struc(:,1);
A1 = p1.F_struc(:,2:end);
b2 = -p2.F_struc(:,1);
A2 = p2.F_struc(:,2:end);
end

alpha = sdpvar(length(xstar),1);
beta = sdpvar(1);
mu1 = sdpvar(length(b1),1);
mu2 = sdpvar(length(b2),1);

Objective  = alpha'*xstar-beta;
Constraint = [alpha' == mu1'*A1,alpha' == mu2'*A2, beta <= mu1'*b1, beta <= mu2'*b2,mu1(neq1+1:end)>0,mu2(neq2+1:end)>0];
Constraint = [Constraint,-1<alpha<1,sum(mu1)+sum(mu2)<10];

solvesdp(Constraint,Objective,sdpsettings('verbose',0));

cut = (-double(beta)+double(alpha)'*x >= 0);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [cData3M, xLim, yLim] = CERRDoseColorWash(hAxis, dose2M, doseXVals, doseYVals, offset, CT2M, CTXVals, CTYVals, scanSet)
% function CERRDoseColorWash"
% Create dose display as a colorwash in axis hAxis, using dose2M as
% defined at doseXVals, doseYVals. CT2M is a B&W CT image defined at
% CTXVals, CTYVals that will have dose interpolated to it and displayed.
%
% CT2M, CTXvals, and CTYVals are optional--if they do not exist, the
% dose is displayed using its own x,y values.
%
% Created: 24 Nov 02, JOD.
%  LM: 03 Dec 02.
%      29 Dec 02, test for 'trans' view before assigning old dose matrix;
%                  also added 'axes commands' to assure correct GUI behavior; JOD.
%      05 Jan 03, JOD.
%      06 Jan 03, CZ (fixed problem of having colorbar marks showing which are larger than maximum dose.)
%      13 Jan 03, CZ visualRef changes
%      13 Jan 03, JOD, minor colorbar marks off.
%      16 Jan 03, JOD, don't use dose size to get CT size ever.
%      20 Jan 03, JOD, visual ref dose isn't show if over dose max.
%      23 Feb 03, JOD, fixed colorbar labels over 100 Gy.
%      07 Apr 03, JOD, adjusted display of strings.
%      09 Apr 03, JOD, fix bug associated with empty levelsV.
%      23 May 03, JRA, added code to draw washes with offset values, and deal with colorbar+labels:May 25, also added hatchpattern
%                       For further changes, see CVS logs.
%
%Usage:
%   [cData3M, xLim, yLim] = CERRDoseColorWash(hAxis, dose2M, doseXVals, doseYVals, offset, CT2M, CTXVals, CTYVals)
%
% 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/>.

global stateS planC
indexS = planC{end};

if strcmpi(get(hAxis,'tag'),'doseCompareAxes')
    colorbarFrameMin = stateS.colorbarFrameMinCompare;
    colorbarFrameMax = stateS.colorbarFrameMaxCompare;
    colorbarRange    = stateS.colorbarRangeCompare;
    doseDisplayRange = stateS.doseDisplayRangeCompare;
else
    colorbarFrameMin = stateS.colorbarFrameMin;
    colorbarFrameMax =  stateS.colorbarFrameMax;
    colorbarRange    = stateS.colorbarRange;
    doseDisplayRange = stateS.doseDisplayRange;
end

if stateS.imageRegistration
    if ~isempty(stateS.doseFusionColormap)
        c = CERRColorMap(stateS.doseFusionColormap);
    else
        c = CERRColorMap(stateS.optS.doseColormap);
    end
else
    c = CERRColorMap(stateS.optS.doseColormap);
end

if stateS.optS.doubleSidedColorbar
    c = [flipud(c);c];
end

[n, m] = size(dose2M);

noCT = 0;

%Get CT data
if ~exist('CT2M','var') || (isempty(CT2M) && isempty(CTXVals))
    CT2M = zeros(size(dose2M));
    CTXVals = doseXVals;
    CTYVals = doseYVals;
    minCol = 1;
    maxCol = m;
    maxRow = n;
    minRow = 1;
    noCT = 1;
else
    dose2M = finterp2(doseXVals, doseYVals, dose2M, CTXVals, CTYVals, 1, 0);

    % get row/col for clipping the dose grid
    minCol = findnearest(CTXVals,min(doseXVals));
    maxCol = findnearest(CTXVals,max(doseXVals));
    maxRow = findnearest(CTYVals,min(doseYVals));
    minRow = findnearest(CTYVals,max(doseYVals));

    % Swap min/max based on the order of x/v vectors for trans/sag/cor views
    minColTmp = min(minCol,maxCol);
    maxCol = max(minCol,maxCol);
    minCol = minColTmp;
    minRowTmp = min(minRow,maxRow);
    maxRow = max(minRow,maxRow);
    minRow = minRowTmp;

%     %% for DDM : Use Kriging to fit dose
%     [CTXValsM, CTYValsM] = meshgrid(CTXVals, CTYVals);
%     [xUnifV, yUnifV, jnk] = getUniformScanXYZVals(planC{indexS.scan}(1));
%     [jnk1, jnk2, zCTV] = getScanXYZVals(planC{indexS.scan}(1));
%     %structureNameAnalyzeC =  {'MMR', 'MML', 'TMR', 'TML', 'PMR', 'PML', 'PLR', 'PLL'};
%     structureNameAnalyzeC =  {'All_Muscles_plus_2cm'};
%     allStructureNames = {planC{indexS.structures}.structureName};
%     structNum = [];
%     for strIndex = 1:length(structureNameAnalyzeC)
%         structNum = [structNum find(strcmpi(structureNameAnalyzeC{strIndex}, allStructureNames))];
%     end
%
%     viewStr = getAxisInfo(hAxis,'view');
%     coord = getAxisInfo(hAxis,'coord');
%     switch lower(viewStr)
%         case 'transverse'
%             dim = 3;
%             sliceNum = findnearest(zCTV,coord);
%
%             [segs, planC, isError] = getRasterSegments(structNum(1), planC, sliceNum);
%             for structNumMask = structNum(2:end)
%                 segsTmp = getRasterSegments(structNumMask, planC, sliceNum);
%                 segs = [segs; segsTmp];
%             end
%
%             maskM = rasterToMask(segs, 1,planC);
%
%         case 'sagittal'
%             dim = 1;
%             sliceNum = findnearest(xUnifV,coord);
%             maskM = getStructureMask(structNum(1), sliceNum, dim, planC);
%             for structNumMask = structNum(2:end)
%                 maskM = maskM | getStructureMask(structNumMask, sliceNum, dim, planC);
%             end
%
%         case 'coronal'
%             dim = 2;
%             sliceNum = findnearest(yUnifV,coord);
%             maskM = getStructureMask(structNum(1), sliceNum, dim, planC);
%             for structNumMask = structNum(2:end)
%                 maskM = maskM | getStructureMask(structNumMask, sliceNum, dim, planC);
%             end
%     end
%     %[iV,jV] =  find(~isnan(dose2M));
%     %dmodelDose = dacefit([doseXVals(jV)' doseYVals(iV)'], dose2M(sub2ind(size(dose2M),iV,jV)), @regpoly0, @correxp, 10, 1e-1, 20);
%     % DoseKrigV = predictor([CTXValsM(:) CTYValsM(:)], dmodelDose);
%
%     [xVals, yVals, zVals] = getDoseXYZVals(planC{indexS.dose}(2));
%     dose3M = planC{indexS.dose}(2).doseArray;
%     [iV,jV,kV] =  find3d(~isnan(dose3M));
%     dataV = dose3M(sub2ind(size(dose3M),iV,jV,kV));
%     dmodelDose = dacefit([xVals(jV)' yVals(iV)' zVals(kV)'], double(dataV(:)), @regpoly0, @correxp, 10, 1e-1, 20);
%
%     CTXValsEval = CTXValsM(maskM);
%     CTYValsEval = CTYValsM(maskM);
%     if dim == 3
%         DoseKrigV = predictor([CTXValsEval(:) CTYValsEval(:) coord*ones(length(CTXValsEval(:)),1)], dmodelDose);
%     elseif dim == 1
%         DoseKrigV = predictor([coord*ones(length(CTXValsEval(:)),1) CTXValsEval(:) CTYValsEval(:)], dmodelDose);
%     else % dim = 2
%         DoseKrigV = predictor([CTXValsEval(:) coord*ones(length(CTXValsEval(:)),1) CTYValsEval(:)], dmodelDose);
%     end
%     dose2M = zeros(size(CTXValsM));
%     dose2M(maskM) = DoseKrigV;

end

ctSize = [size(CT2M,1), size(CT2M,2)];
%
%Mask of pixels w/dose above the cutoff.
lowerBound = doseDisplayRange(1);
upperBound = doseDisplayRange(2);

if(offset)
    negativeMask = [dose2M < offset];
    maskM = [dose2M <= (upperBound+offset) & dose2M >= (lowerBound+offset)];
else
    maskM = [dose2M <= upperBound & dose2M >= lowerBound];
end

if stateS.optS.transparentZeroDose
    %     if offset == 0  %cuts-off dose=0
    %         maskM(dose2M <= offset) = 0;
    %     end
    if offset
        maskM(dose2M == offset) = 0;
    else
        maskM(dose2M == 0) = 0;
    end
end

if stateS.optS.calcDoseInsideSkinOnly
    [structNum,isSkinUniform] = getStructureIndex(hAxis,'skin',planC);
    if isSkinUniform
        skinMaskM = getSkinMask(hAxis,structNum,planC);
        maskM = maskM & skinMaskM;
        %dose2M(skinMaskM) = 0;
    end
end

% Clip dose grid
maskTmpM = maskM;
maskM = false(size(maskM));
maskM(minRow:maxRow,minCol:maxCol) = maskTmpM(minRow:maxRow,minCol:maxCol);

%Fit dose to colormap, but first replace full dose with only dose values that will be displayed.
lowerBound = colorbarRange(1);
upperBound = colorbarRange(2);
if ~(lowerBound == 0 && upperBound == 0) && ~(lowerBound == upperBound)
    percentBelow = (colorbarRange(1) - colorbarFrameMin) / (colorbarRange(2) - colorbarRange(1));
    percentAbove = (colorbarFrameMax - colorbarRange(2)) / (colorbarRange(2) - colorbarRange(1));
    nElements = size(c, 1);
    lowVal = c(1,:);
    hiVal = c(end,:);

    c = [repmat(lowVal, [round(percentBelow*nElements),1]);c;repmat(hiVal, [round(percentAbove*nElements),1])];
end

partialDose2M = dose2M(maskM);

partialDose =((partialDose2M-(colorbarFrameMin+offset))/(colorbarFrameMax - colorbarFrameMin)) * (size(c,1) + 0.5);

roundPartialDose = round(partialDose);
partialDoseClip = clip(roundPartialDose,1,size(c,1),'limits');

%build RGB Matrix by indexing into colormap.
cData3M = c(partialDoseClip, 1:3);

%Temporarily out of commission: doesnt matter since we no longer use texture heavily.
if(offset) && stateS.optS.negativeTexture
    textureMap = [];
    if isempty(textureMap)
        thin = 1; %intensity values for thin and thick bars
        thick = 0.95;
        [sizeX,sizeY] = size(dose2M);
        textureMap = repmat([thin, thin, thick, thick;
            thick, thin, thin, thick;
            thick, thick, thin, thin;
            thin, thick, thick, thin], ceil(sizeX/4), ceil(sizeY/4));

        %resize texturemap to fit the size of the dose.
        textureMap = textureMap(1:sizeX, 1:sizeY);
    end

    negativeValues = negativeMask;
    negativeValuesM = repmat(negativeValues, [1 1 3]);
    textureM = repmat(textureMap, [1 1 3]);
    mask3M = repmat(maskM, [1 1 3]);
    cData3M(negativeValuesM(mask3M)) = cData3M(negativeValuesM(mask3M)) .* (textureM(negativeValuesM & mask3M));
end

CTBackground3M = [];

if stateS.CTToggle == 1 && ~noCT %Don't show very low doses

    CTOffset = planC{indexS.scan}(scanSet(1)).scanInfo(1).CTOffset;

    scanUID = ['c',repSpaceHyp(planC{indexS.scan}(scanSet).scanUID(max(1,end-61):end))];

    colorCT = CERRColorMap(stateS.scanStats.Colormap.(scanUID));
    CTLevel     = stateS.scanStats.CTLevel.(scanUID) + CTOffset;
    CTWidth     = stateS.scanStats.CTWidth.(scanUID);
    CTLow       = CTLevel - CTWidth/2;
    CTHigh      = CTLevel + CTWidth/2;

    %CT2M = clip(CT2M, CTLow, CTHigh, 'limits');

    % Get Min/max of CT2M scan to scale accordingly.
    %minCT = min(CT2M(:));
    %maxCT = max(CT2M(:));
    %scanMin = stateS.scanStats.minScanVal.(scanUID);
    %scanMax = stateS.scanStats.maxScanVal.(scanUID);
    %CTLow = max(CTLow,scanMin);
    %CTHigh = max(CTLow,min(CTHigh,scanMax));

    CT2M = clip(CT2M, CTLow, CTHigh, 'limits');

    %This is a trick for speed.  Map the CT data from 1...N+1 bins, which
    %results (only) the maxValue exceeding N after floored.  Replicate
    %the colorCT last element to display the maxValue correctly.
    if CTLow ~= CTHigh
        ctScaled = (CT2M - CTLow) / ((CTHigh - CTLow) / size(colorCT,1)) + 1;
        %ctScaled = (CT2M - minCT) / ((maxCT - minCT) / size(colorCT,1)) + 1; % bug
        ctClip = uint16(ctScaled);
        colorCT(end+1,:) = colorCT(end,:);
    else
        ctClip = ones(size(CT2M));
    end
    %ctClip = ctClip(1:ctSize(1),1:ctSize(2));

    %colorCT = (1-stateS.doseAlphaValue.trans)*colorCT; % APA comment

    %build CT background by indexing into CTcolormap. Optimal mtd for speed.
    %CTBackground3M = reshape((1-stateS.doseAlphaValue.trans)*colorCT(ctClip,1:3),ctSize(1),ctSize(2),3);
    CTBackground3M = colorCT(ctClip(:),1:3);

    %Check status of transparency value.
    if ~isfield(stateS, 'doseAlphaValue')
        stateS.doseAlphaValue.trans = stateS.optS.initialTransparency;
    end

    %Build 3 layer mask of dose pixel locations.
    %mask3M = repmat(maskM, [1 1 3]);
    mask3M = repmat(maskM(:), [1 3]);

    %Add dose to the CT, merging it with CT values @mask3M based on alpha.
    if ~isempty(cData3M(:))
        % CTBackground3M(mask3M) = (cData3M(:) .* stateS.doseAlphaValue.trans) + (CTBackground3M(mask3M) .* (1-stateS.doseAlphaValue.trans));
        if stateS.imageRegistration
            CTBackground3M(mask3M) = [cData3M(:) CTBackground3M(mask3M)] * [stateS.doseFusionAlpha 1-stateS.doseFusionAlpha]';
        else
            %CTBackground3M(mask3M) = [cData3M(:) CTBackground3M(mask3M)] * [stateS.doseAlphaValue.trans 1-stateS.doseAlphaValue.trans]';
            CTBackground3M(mask3M) = (1-stateS.doseAlphaValue.trans)*CTBackground3M(mask3M) + cData3M(:)*stateS.doseAlphaValue.trans;
        end
    end
    CTBackground3M = reshape(CTBackground3M,ctSize(1),ctSize(2),3);
    cData3M = CTBackground3M;
end

if ~isempty(CTBackground3M)
    cData3M = CTBackground3M;
else
    CTBackground3M = zeros(ctSize(1), ctSize(2), 3);
    mask3M = repmat(maskM, [1 1 3]);
    CTBackground3M(mask3M) = cData3M(:);
    cData3M = CTBackground3M;
end

% If in print mode, replace area outside of skin with white. OR dose
% outside skin is set to zero. the code repeats so just take the same code
if stateS.printMode % | stateS.optS.calcDoseInsideSkinOnly

    %     structNum = strmatch('skin',lower({planC{indexS.structures}.structureName}),'exact');
    %     [assocScan, relStructNum] = getStructureAssociatedScan(structNum, planC);
    %     isSkinUniform = any(bitget(planC{indexS.structureArray}(assocScan).bitsArray,relStructNum));

    %%%%%%%%%%% DK replace above two lines with the following code to get
    %%%%%%%%%%% the actual structure number in case there are multiple
    %%%%%%%%%%% structures with the same name.
%     aI = get(hAxis,'userdata');
%
%     scanSet = aI.scanSets;
%
%     assocScansV = getStructureAssociatedScan(1:length(planC{indexS.structures}), planC);
%
%     indXStr = find(assocScansV == scanSet);
%
%     allStructureNames = {planC{indexS.structures}(indXStr).structureName};
%
%     structNum = find(strcmpi('skin', allStructureNames));
%
%     if structNum <= 52
%         bitsArray = planC{indexS.structureArray}(scanSet(1)).bitsArray;
%         isSkinUniform = any(bitget(bitsArray,structNum));
%     else
%         cellNum = ceil((structNum-52)/8)+1;
%         bitsArray = planC{indexS.structureArrayMore}(scanSet(1)).bitsArray{cellNum-1};
%         isSkinUniform = any(bitget(bitsArray,structNum-52-(cellNum-2)*8));
%     end
%
%     if ~isSkinUniform
%         if ~isempty(structNum)
%             CERRStatusString('Skin structure is not uniformized. Unable to turn dose off outside skin.')
%         end
%         structNum = [];
%     end

    [structNum,isSkinUniform] = getStructureIndex(hAxis,'skin',planC);

    %%%%%%%%%%%%++++++++++++++++++++++++++++++++++++++++%%%%%%%%%%%%%%%%%%%

%     %obtain x,y and z-vals
%     [xUnifV, yUnifV, jnk] = getUniformScanXYZVals(planC{indexS.scan}(scanSet));
%     [jnk1, jnk2, zCTV] = getScanXYZVals(planC{indexS.scan}(scanSet));

    if ~isempty(structNum)

%         viewStr = getAxisInfo(hAxis,'view');
%         coord = getAxisInfo(hAxis,'coord');
%         switch lower(viewStr)
%             case 'transverse'
%                 dim = 3;
%                 sliceNum = findnearest(zCTV,coord);
%
%                 [segs, planC, isError] = getRasterSegments(structNum, planC, sliceNum);
%
%                 maskM = rasterToMask(segs, scanSet,planC);
%
%             case 'sagittal'
%                 dim = 1;
%                 sliceNum = findnearest(xUnifV,coord);
%                 maskM = getStructureMask(structNum, sliceNum, dim, planC);
%
%
%             case 'coronal'
%                 dim = 2;
%                 sliceNum = findnearest(yUnifV,coord);
%                 maskM = getStructureMask(structNum, sliceNum, dim, planC);
%
%         end

        maskM = getSkinMask(hAxis,structNum,planC);

        mask3M = repmat(maskM, [1 1 3]);

        if ((isSkinUniform && dim~=3) || dim==3) & stateS.printMode
            set(hAxis,'color',[1 1 1])

            if strcmp(class(cData3M),'uint8')
                cData3M(~mask3M) = 255;
            elseif stateS.printMode
                cData3M(~mask3M) = 1;
            end
        end

%         if stateS.optS.calcDoseInsideSkinOnly && ~stateS.printMode && ((isSkinUniform && dim~=3) || dim==3)
%             cData3M(~mask3M) = 0;
%         end
    end
end

%hold on
set(hAxis, 'nextPlot', 'add');

xLim = CTXVals([1 end]);
yLim = CTYVals([1 end]);

return;

function maskM = getSkinMask(hAxis,structNum,planC)

indexS = planC{end};

%aI = get(hAxis,'userdata');
scanSet = getAxisInfo(hAxis,'scanSets');
viewStr = getAxisInfo(hAxis,'view');
coord = getAxisInfo(hAxis,'coord');

%obtain x,y and z-vals
[xUnifV, yUnifV, jnk] = getUniformScanXYZVals(planC{indexS.scan}(scanSet));
[jnk1, jnk2, zCTV] = getScanXYZVals(planC{indexS.scan}(scanSet));

switch lower(viewStr)
    case 'transverse'
        dim = 3;
        sliceNum = findnearest(zCTV,coord);

        [segs, planC, isError] = getRasterSegments(structNum, planC, sliceNum);

        maskM = rasterToMask(segs, scanSet,planC);

    case 'sagittal'
        dim = 1;
        sliceNum = findnearest(xUnifV,coord);
        maskM = getStructureMask(structNum, sliceNum, dim, planC);

    case 'coronal'
        dim = 2;
        sliceNum = findnearest(yUnifV,coord);
        maskM = getStructureMask(structNum, sliceNum, dim, planC);

end
return;

function [structNum,isSkinUniform] = getStructureIndex(hAxis,structName,planC)

indexS = planC{end};

aI = getAxisInfo(hAxis);

scanSet = aI.scanSets;
assocScansV = getStructureAssociatedScan(1:length(planC{indexS.structures}), planC);

indXStr = find(assocScansV == scanSet);

allStructureNames = {planC{indexS.structures}(indXStr).structureName};

structNum = find(strcmpi(structName, allStructureNames));

if structNum <= 52
    bitsArray = planC{indexS.structureArray}(scanSet(1)).bitsArray;
    isSkinUniform = any(bitget(bitsArray,structNum));
else
    cellNum = ceil((structNum-52)/8)+1;
    bitsArray = planC{indexS.structureArrayMore}(scanSet(1)).bitsArray{cellNum-1};
    isSkinUniform = any(bitget(bitsArray,structNum-52-(cellNum-2)*8));
end

if ~isSkinUniform
    if ~isempty(structNum)
        CERRStatusString([structName, ' structure is not uniformized. Unable to turn dose off outside ',structName])
    end
    structNum = [];
end
return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function x = r8vec_chebyspace ( n, a, b )

%*****************************************************************************80
%
%% R8VEC_CHEBYSPACE creates a vector of Chebyshev spaced values in [A,B].
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 June 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, real A, B, the first and last entries.
%
%    Output, real X(N,1), a vector of Chebyshev spaced data.
%
  x = zeros ( n, 1 );

  if ( n == 1 )

    x(1) = ( a + b ) / 2.0;

  else

    for i = 1 : n

      theta = ( n - i ) * pi / ( n - 1 );

      c = cos ( theta );

      if ( mod ( n, 2 ) == 1 )
        if ( 2 * i - 1 == n )
          c = 0.0;
        end
      end

      x(i) = ( ( 1.0 - c ) * a  ...
             + ( 1.0 + c ) * b ) ...
             /   2.0;

    end

  end

  return
end


text:
% ------------------------------------------------------
% SwarmOps - Heuristic optimization for Matlab
% Copyright (C) 2003-2010 Magnus Erik Hvass Pedersen.
% Please see the file license.txt for license details.
% SwarmOps on the internet: http://www.Hvass-Labs.org/
% ------------------------------------------------------

% Example optimization problem. You may use this as
% a starting point for custom problems.
% Parameters:
%     x; position in the search-space.
%     data; data-struct for optimization problem.
% Returns:
%     fitness; the measure to be minimized.
function fitness = myproblem(x, data)
    % Retrieve data from struct.
    r = data.MyExtraData;

    % Displace position.
    t = x-r;

    % Compute and return fitness.
    fitness = (2*t(1)-t(2))^2 + ...
              (3*t(2)-2*t(3))^2 + ...
              (4*t(3)-3*t(4))^2 + ...
              (t(4)-4*t(1))^2 + ...
              sum(t.^2);
end

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


text:
function marginal = marginal_family(engine, i, t)
% MARGINAL_FAMILY Compute the marginal on the specified family (ff)
% marginal = marginal_family(engine, i, t)

if nargin < 3, t = 1; end

% The method is similar to the following HMM equation:
% xi(i,j,t) = normalise( alpha(i,t) * transmat(i,j) * obsmat(j,t+1) * beta(j,t+1) )
% where xi(i,j,t) = Pr(Q(t)=i, Q(t+1)=j | y(1:T))

bnet = bnet_from_engine(engine);

if myismember(i, engine.onodes)
  ps = parents(bnet.dag, i);
  p = ps(1);
  marginal = pot_to_marginal(engine.marginals{p,t});
  marginal.domain = [p i];
  return;
end

if t==1
  marginal = pot_to_marginal(engine.marginals{i,t});
  return;
end

bnet = bnet_from_engine(engine);
ss = length(bnet.intra);
pot = engine.CPDpot{i,t};
c = engine.obschild(i);
pot = multiply_by_pot(pot, engine.CPDpot{c,t});
pot = multiply_by_pot(pot, engine.back{i,t});
ps = parents(bnet.dag, i+ss);
for p=ps(:)'
  pot = multiply_by_pot(pot, engine.fwd{p,t-1});
end
marginal = pot_to_marginal(normalize_pot(pot));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function out=hasimag(in)

out=false;
if any(in~='0')
 out=true;
end

author: opencobra

repo: cobratoolbox


text:
function obj = coeftq1_pnp_func_new(in1)
coef_J11 = in1(:,11);
coef_J20 = in1(:,20);
coef_J30 = in1(:,30);
coef_J26 = in1(:,26);
coef_J45 = in1(:,45);
coef_J63 = in1(:,63);
coef_J55 = in1(:,55);
coef_J39 = in1(:,39);
coef_J49 = in1(:,49);
coef_J59 = in1(:,59);
obj = [coef_J11,coef_J20,coef_J26,coef_J30,coef_J39,coef_J45,coef_J49,coef_J55,coef_J59,coef_J63];