text:
function dist = line_par_point_dist_3d ( f, g, h, x0, y0, z0, p )

%*****************************************************************************80
%
%% LINE_PAR_POINT_DIST_3D: distance ( parametric line, point ) in 3D.
%
%  Discussion:
%
%    The parametric form of a line in 3D is:
%
%      X = X0 + F * T
%      Y = Y0 + G * T
%      Z = Z0 + H * T
%
%    We normalize by choosing F*F+G*G+H*H=1 and 0 <= F.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Adrian Bowyer and John Woodwark,
%    A Programmer's Geometry,
%    Butterworths, 1983.
%
%  Parameters:
%
%    Input, real F, G, H, X0, Y0, Z0, the parametric line
%    parameters.
%
%    Input, real P(3,1), the point whose distance from the line is
%    to be measured.
%
%    Output, real DIST, the distance from the point to the line.
%
  dx =   g * ( f * ( p(2,1) - y0 ) - g * ( p(1,1) - x0 ) ) ...
       + h * ( f * ( p(3,1) - z0 ) - h * ( p(1,1) - x0 ) );

  dy =   h * ( g * ( p(3,1) - z0 ) - h * ( p(2,1) - y0 ) ) ...
       - f * ( f * ( p(2,1) - y0 ) - g * ( p(1,1) - x0 ) );

  dz = - f * ( f * ( p(3,1) - z0 ) - h * ( p(1,1) - x0 ) ) ...
       - g * ( g * ( p(3,1) - z0 ) - h * ( p(2,1) - y0 ) );

  dist = sqrt ( dx * dx + dy * dy + dz * dz ) ...
    / ( f * f + g * g + h * h );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function STResNet_stage2(varargin)

if ~isempty(gcp('nocreate'))
    delete(gcp)
end

opts.train.gpus = [ 1 ] ;

% addpath('../network_surgery');
opts = cnn_setup_environment(opts);

% opts.dataSet = 'hmdb51';
opts.dataSet = 'ucf101';

opts.dropOutRatio = NaN;
opts.backpropFuseFrom = 1;
opts.nSplit =  1 ;
opts.addPool3D = 1 ;
opts.singleSoftMax = 0;
opts.nFrames = 11 ;
% opts.train.backpropDepth = 'pool5';
opts.train.learningRate =  1*[ 1e-4*ones(1,1) 1e-5*ones(1,1)]  ;
opts.train.augmentation = 'f25noCtr';

opts.train.epochFactor = 1 ;

opts.train.batchSize = 128  ;
opts.train.numSubBatches = 32 ;
opts.train.cheapResize = 0 ;
opts.poolMethod = 'max';
opts.poolSz=5 ;
opts.poolStride=2 ;

model = ['ST-ResNet50-final-split=' num2str(opts.nSplit)];

opts.train.numSubBatches =  ceil(opts.train.numSubBatches / max(numel(opts.train.gpus),1));

opts.train.memoryMapFile = fullfile(tempdir, 'ramdisk', ['matconvnet' num2str(opts.nSplit ) '.bin']) ;
opts.dataDir = fullfile(opts.dataPath, opts.dataSet) ;
opts.splitDir = 'ucf101_splits'; nClasses = 101;
opts.imdbPath = fullfile(opts.dataDir, [opts.dataSet '_split' num2str(opts.nSplit) 'imdb.mat']);

opts.model = fullfile(opts.modelPath, [opts.dataSet '-ST-ResNet50-split=' num2str(opts.nSplit) '.mat']) ;

if strcmp(opts.dataSet, 'hmdb51')
  opts.splitDir = 'hmdb51_splits'; nClasses = 51;
  opts.flowDir = strrep(opts.flowDir, 'ucf101','hmdb51');
  opts.imdbPath = fullfile(opts.dataDir, ['hmdb_split' num2str(opts.nSplit) 'imdb.mat']);
end

opts.expDir = fullfile(opts.dataDir, [opts.dataSet '-' model]) ;
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.train.saveAllPredScores = 1;
opts.train.denseEval = 1;
opts.train.plotDiagnostics = 0 ;
opts.train.continue = 1 ;
opts.train.prefetch = 1 ;
opts.train.expDir = opts.expDir ;
opts.train.numAugments = 1;
opts.train.frameSample = 'random';
opts.train.nFramesPerVid = 1;

opts = vl_argparse(opts, varargin) ;

% -------------------------------------------------------------------------
%                                                   Database initialization
% -------------------------------------------------------------------------
if exist(opts.imdbPath)
  imdb = load(opts.imdbPath) ;
  imdb.flowDir = opts.flowDir;
else
  imdb = cnn_setup_data(opts) ;
  save(opts.imdbPath, '-struct', 'imdb', '-v6') ;
end

% -------------------------------------------------------------------------
%                                                    Network initialization
% -------------------------------------------------------------------------
if ~exist(opts.model,'file')
  [~, baseModel] = fileparts(opts.model);
  fprintf('Downloading base model file: %s ...\n', baseModel);
  mkdir(fileparts(opts.model)) ;
  urlwrite(...
  ['http://ftp.tugraz.at/pub/feichtenhofer/st-res/stage1/' baseModel '.mat'], ...
    opts.model) ;
end
net = load(opts.model) ;

if isfield(net, 'net'), net=net.net;end
net = dagnn.DagNN.loadobj(net);

if opts.addPool3D

%   poolLayers = {'res2a', 'res3a',  'res4a', 'res5a'; };
%   poolLayers = {'pool5'; };
  poolLayers = {'res5c_relu'; };

  for j=1:numel(poolLayers)
    for s = {'spatial', 'temporal'}
      i_pool = find(strcmp({net.layers.name},[poolLayers{j} '_' char(s)]));

      block = dagnn.PoolTime() ;
      block.poolSize = [1 opts.poolSz];
      block.pad = [0 0 0 0];
      block.stride = [1 opts.poolStride];
      block.method = opts.poolMethod;

      name = [poolLayers{j} '_pool_' char(s)];

      disp(['injecting ' name ' as PoolTime'])

      net.addLayerAt(i_pool, name, block, ...
                    [net.layers(i_pool).outputs], {name}) ;

      % chain input of l that has layer as input
      for l = 1:numel(net.layers)
          if ~strcmp(net.layers(l).name, name)
            sel = find(strcmp(net.layers(l).inputs, net.layers(i_pool).outputs{1})) ;
            if any(sel)
                net.layers(l).inputs{sel} = name;
            end;
          end
      end

    end
  end
end % add pool3d

if opts.addPool3D
  opts.train.augmentation = 'f25noCtr';
  opts.train.frameSample = 'temporalStrideRandom';
  opts.train.nFramesPerVid = opts.nFrames * 1;
  opts.train.temporalStride = 1:15;
  opts.train.valmode = 'temporalStrideRandom';
  opts.train.numValFrames = 25 ;
  opts.train.saveAllPredScores = 1 ;
  opts.train.denseEval = 1;
  opts.train.temporalFullConvTest = 1;
end

net.meta.normalization.rgbVariance = [];
opts.train.train = find(ismember(imdb.images.set, [1])) ;
opts.train.train = repmat(opts.train.train,1,opts.train.epochFactor);

zero_drs =  find(arrayfun(@(x) isa(x.block,'dagnn.DropOut') && x.block.rate == 0, net.layers)) ;
net.removeLayer({net.layers(zero_drs).name});

if ~isnan(opts.dropOutRatio)
  dr_layers = find(arrayfun(@(x) isa(x.block,'dagnn.DropOut'), net.layers)) ;
  if opts.dropOutRatio > 0
    net.layers(dr_layers).block.rate = opts.dropOutRatio;
  else
    net.removeLayer({net.layers(dr_layers).name});
  end
end

if opts.singleSoftMax
  pred_layers = [];
  for l=1:numel(net.layers)
    if isempty( net.layers(l).params ), continue; end;
    if size(net.params(net.getParamIndex(net.layers(l).params{1})).value,4) == nClasses || ...
        size(net.params(net.getParamIndex(net.layers(l).params{1})).value,5) == nClasses % 3D FC layer
          pred_layers = [pred_layers l];
          net.vars(net.layers(l).outputIndexes).precious = 1;
    end
  end
  pred_layers = fliplr(pred_layers) ; % remove the spatial layer
  paramsIdx1 = net.getParamIndex(net.layers(pred_layers(1)).params);
  paramsIdx2 = net.getParamIndex(net.layers(pred_layers(2)).params);
  for p = 1:numel(paramsIdx1)
    sz = size(net.params(paramsIdx1(p)).value);
    if numel(sz) > 2
      net.params(paramsIdx1(p)).value = cat(3,net.params(paramsIdx1(p)).value, net.params(paramsIdx2(p)).value);
    else
      net.params(paramsIdx1(p)).value = net.params(paramsIdx1(p)).value + net.params(paramsIdx2(p)).value;
    end
  end
  block = dagnn.Concat() ;
  newName = ['singleSoftMaxConcat'];
  net.addLayer(newName, block, ...
               [net.layers(pred_layers).inputs], ...
               newName) ;
  net.layers(pred_layers(1)).inputs =    newName ;
  % remove layers of the other prediction
  for l = numel(net.layers):-1:1
    for f = net.layers(pred_layers(2)).outputs
       sel = find(strcmp(f, net.layers(l).inputs )) ;
       if ~isempty(sel)
         fprintf('removing ayer %s \n', net.layers(l).name);
         net.removeLayer({net.layers(l).name});
       end
    end
  end
  net.removeLayer({net.layers(pred_layers(2)).name});
end

net.layers(~cellfun('isempty', strfind({net.layers(:).name}, 'err'))) = [] ;

net.rebuild() ;

opts.train.derOutputs = {} ;
for l=1:numel(net.layers)
  if isa(net.layers(l).block, 'dagnn.Loss') && isempty(strfind(net.layers(l).block.loss, 'err'))
    if opts.backpropFuseFrom || ~isempty(strfind(net.layers(l).name, opts.train.fuseInto ))
      fprintf('setting derivative for layer %s \n', net.layers(l).name);
      opts.train.derOutputs = [opts.train.derOutputs, net.layers(l).outputs, {1}] ;
    end
    net.addLayer(['err1_' net.layers(l).name(end-7:end) ], dagnn.Loss('loss', 'classerror'), ...
             net.layers(l).inputs, 'error') ;
  end
end

net.conserveMemory = 1 ;
fn = getBatchWrapper_rgbflow(net.meta.normalization, opts.numFetchThreads, opts.train) ;
[info] = cnn_train_dag(net, imdb, fn, opts.train) ;

author: feichtenhofer

repo: st-resnet


text:
classdef matRad_MaxDVH < DoseObjectives.matRad_DoseObjective
% matRad_MaxDVH Implements a penalized maximum DVH objective
%   See matRad_DoseObjective for interface description
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2020 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    properties (Constant)
        name = 'Max DVH';
        parameterNames = {'d', 'V^{max}'};
        parameterTypes = {'dose','numeric'};
    end

    properties
        parameters = {30,95};
        penalty = 1;
    end

    methods
        function obj = matRad_MaxDVH(penalty,dRef,vMaxPercent)
            %If we have a struct in first argument
            if nargin == 1 && isstruct(penalty)
                inputStruct = penalty;
                initFromStruct = true;
            else
                initFromStruct = false;
                inputStruct = [];
            end

            %Call Superclass Constructor (for struct initialization)
            obj@DoseObjectives.matRad_DoseObjective(inputStruct);

            %now handle initialization from other parameters
            if ~initFromStruct
                if nargin >= 3 && isscalar(vMaxPercent)
                    obj.parameters{2} = vMaxPercent;
                end

                if nargin >= 2 && isscalar(dRef)
                    obj.parameters{1} = dRef;
                end

                if nargin >= 1 && isscalar(penalty)
                    obj.penalty = penalty;
                end
            end
        end

        %% Calculates the Objective Function value
        function fDose = computeDoseObjectiveFunction(obj,dose)
            % get reference Volume
            refVol = obj.parameters{2}/100;

            % calc deviation
            deviation = dose - obj.parameters{1};

            % calc d_ref2: V(d_ref2) = refVol
            d_ref2 = matRad_calcInversDVH(refVol,dose);

            deviation(dose < obj.parameters{1} | dose > d_ref2) = 0;

            % claculate objective function
            fDose = (obj.penalty/numel(dose))*(deviation'*deviation);
        end

        %% Calculates the Objective Function gradient
        function fDoseGrad   = computeDoseObjectiveGradient(obj,dose)
            % get reference Volume
            refVol = obj.parameters{2}/100;

            % calc deviation
            deviation = dose - obj.parameters{1};

            % calc d_ref2: V(d_ref2) = refVol
            d_ref2 = matRad_calcInversDVH(refVol,dose);

            deviation(dose < obj.parameters{1} | dose > d_ref2) = 0;

            % calculate delta
            fDoseGrad = 2 * (obj.penalty/numel(dose))*deviation;
        end
    end

end

repo: matRad


text:
%% Ex. 8 Another example of elementary functions with a vectorial variable

a = [2 3 5];
b = 2*a.^2+3*a+4

%Output:
%        b = 18 31 69
%   Remark: The content of b is
%    b = [2*(a(1))^2+3*a(1)+4 2*(a(2))^2+3*a(2)+4 2*(a(3))^2+3*a(3)+4].

author: TheAlgorithms

repo: MATLAB-Octave


text:
function [ element_node, element_att ] = triangle_element_data_example ( ...
  element_num, element_order, element_att_num )

%*****************************************************************************80
%
%% TRIANGLE_ELEMENT_DATA_EXAMPLE returns the element information for the example.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 October 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_ORDER, the order of the elements.
%
%    Input, integer ELEMENT_ATT_NUM, the number of element
%    attributes.
%
%    Output, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM),
%    the indices of the nodes that make up each element.
%
%    Output, real ELEMENT_ATT(ELEMENT_ATT_NUM,ELEMENT_NUM),
%    the attributes of each element.
%
  element_att = [];

  element_node = [ ...
    1,  2,  6; ...
    7,  6,  2; ...
    2,  3,  7; ...
    8,  7,  3; ...
    3,  4,  8; ...
    9,  8,  4; ...
    4,  5,  9; ...
   10,  9,  5; ...
    6,  7, 11; ...
   12, 11,  7; ...
    7,  8, 12; ...
   13, 12,  8; ...
    8,  9, 13; ...
   14, 13,  9; ...
    9, 10, 14; ...
   15, 14, 10; ...
   11, 12, 16; ...
   17, 16, 12; ...
   12, 13, 17; ...
   18, 17, 13; ...
   16, 17, 19; ...
   20, 19, 17; ...
   17, 18, 20; ...
   21, 20, 18 ]';

  return
end


function [ err ] = residual_KR( X1, X2, K1, K2, R)

% homopraphy matrix
H = K2*R/K1;
n = size(X1,2);

x1 = X1(1,:); y1 = X1(2,:);
x2 = X2(1,:); y2 = X2(2,:);

w = H(3,1)*x1 + H(3,2)*y1 + H(3,3);
x1_2 = (H(1,1)*x1 + H(1,2)*y1 + H(1,3)) ./ w;
y1_2 = (H(2,1)*x1 + H(2,2)*y1 + H(2,3)) ./ w;
err1_2 = [x1_2 - x2; y1_2 - y2];

err = reshape(err1_2,[2*n,1]);

end


text:
function F = filmStrip( I, overlap, delta, border )
% Used to display R stacks of T images as a "filmstrip".
%
% See examples below to see what is meant by "filmstrip".
%
% USAGE
%  F = filmStrip( I, overlap, delta, border )
%
% INPUTS
%  I          - MxNxTxR or MxNx1xTxR or MxNx3xTxR array
%               (of bw or color images). R can equal 1.
%  overlap    - amount of overlap between successive frames
%  delta      - amount to shift each successive frame upward
%  border     - width of black border around each frame
%
% OUTPUTS
%  F       - filmstrip
%
% EXAMPLE - one filmstrip
%  load images;
%  F1 = filmStrip( video(:,:,1:15), 10, 2, 5 );   figure(1); im(F1); % one
%  F2 = filmStrip( videos(:,:,:,1:10), 5, 2, 3 ); figure(2); im(F2); % many
%
% See also MONTAGE2
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

% convert I to be of type double, and have dimensions MxNxCxTxR
I = double(I); I = I/max(I(:)); sizI = size(I);
if(~any(sizI(3)==[1 3])); I=reshape(I,[sizI(1:2),1,sizI(3:end)]); end
I = padarray( I, [border border 0 0 0], 0, 'both' );
[mRows, nCols, nColor, nFrame, nStrip] = size(I);

% size of final filmstip object
sizF1 = [mRows+delta*(nFrame-1), nFrame*nCols-overlap*(nFrame-1), nColor];
sizF=sizF1;  sizF(1)=sizF1(1)*nStrip - ((nFrame-1-2)*delta)*(nStrip-1);
mRowsF1 = sizF1(1);

for i=1:nStrip

  % Create i-th film strip
  Fi = -ones( sizF1 );  row = 1;  col = sizF1(2);
  for f=nFrame:-1:1
    Fi( row:(row+mRows-1), (col-nCols+1):col, : ) = I(:,:,:,f,i);
    row = row + delta;  col = col - nCols + overlap;
  end

  % stop if creating single filmstrip
  if( nStrip==1 ); F=Fi; break; end

  % merge with the previous film strips
  if( i==1 ); F = -ones( sizF );  row2=1;  end
  Fc = F( row2:(row2+mRowsF1-1), : );
  locs=(Fc<0);  Fc(locs) = Fi(locs);
  F( row2:(row2+mRowsF1-1), :  ) = Fc;
  row2 = row2 + mRowsF1 - ((nFrame-1-2)*delta);
end
F(F<0)=1;

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
%% INCREMENTAL TENSOR LEARNING
%%% Incremental tensor analysis: Theory and applications (2008)
%%% http://dl.acm.org/citation.cfm?id=1409621
close all; clear all; clc;

%% LIBRARIES
addpath('libs/poblano_toolbox_1.1');
addpath('libs/tensor_toolbox_2.5');
addpath('libs/nway331');
addpath('libs/itl');

%% TENSORS
A = full(ttensor(tenrand([2,3,4]),{rand(10,2),rand(30,3), rand(40,4)}));
B = full(ttensor(tenrand([2,3,4]),{rand(10,2),rand(30,3), rand(40,4)}));
D = {A,A,A,A,A,B,B,B,B,B};

%% Dynamic Tensor Decomposition (DTA)
%%% Compute DTA with no forgetting factor
[T,C] = DTA(D{1},[2,3,4]);
for i = 2:10
  [T,C] = DTA(D{i},[2,3,4],C);
  err = norm(full(T)-D{i})/norm(D{i});
  fprintf('tensor #%d has error %f\n',i,err);
end
%%% Compute DTA with forgetting factor alpha = 0.1
[T,C] = DTA(D{1},[2,3,4]);
alpha = 0.1;
for i = 2:10
  [T,C] = DTA(D{i},[2,3,4],C, alpha);
  err = norm(full(T)-D{i})/norm(D{i});
  fprintf('tensor #%d has error %f\n',i,err);
end

%% Streaming Tensor Decomposition (STA)
%%% Compute STA with no forgetting factor
[T,S] = STA(D{1},[2,3,4]);
for i = 2:10
  [T,S] = STA(D{i},[2,3,4],T,S);
  err = norm(full(T)-D{i})/norm(D{i});
  fprintf('tensor #%d has error %f\n',i,err);
end
%%% Compute STA with forgetting factor alpha = 0.9
[T,S] = STA(D{1},[2,3,4]);
alpha = 0.9;
for i = 2:10
  [T,S] = STA(D{i},[2,3,4],T,S, alpha);
  err = norm(full(T)-D{i})/norm(D{i});
  fprintf('tensor #%d has error %f\n',i,err);
end

%% Window-based Tensor Decomposition (WTA)
[T,C] = WTA(D{1},[2,3,4]);
for i = 2:10
  [T,C] = WTA(D{i},[2,3,4],tensor(T));
  err = norm(full(T)-D{i})/norm(D{i});
  fprintf('tensor #%d has error %e\n',i,err);
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
%% backup the current results before merging
neuron_bk = neuron.copy();
ctr = neuron.estCenter();
K = size(neuron.A,2);   %neuron numbers
ind_del = false(K, 1);
Coor = neuron.get_contours();
[img, col0, AA] = neuron.overlapA([], 0.1);  % overlap neurons
color_all = zeros(K,3);
col = col0;
for m=1:3
color_all(:,m) = mod(col,2);
col = floor(col/2);
end

%% UI CONTROL
IDs = [];
fig_merge = figure('position', [10, 300, 1400, 500], 'papersize', [14, 5]);
ax_all = axes('position', [0.005 0.1, 0.4, 0.8]);
axes(ax_all);
imagesc(img, 'parent', ax_all);  hold on;
axis equal off tight;
pb_select = uicontrol(fig_merge, 'style', 'pushbutton', ...
    'unit', 'normalized',  'position', [0.05, 0.01, 0.1, 0.08], ...
    'string', 'add neurons', 'fontweight', 'bold', ...
    'fontsize', 16,     'callback', 'cnmfe_manual_merge_pick_neuron;');
pb_zoom = uicontrol(fig_merge, 'style', 'pushbutton', ...
    'unit', 'normalized',  'position', [0.16, 0.01, 0.1, 0.08], ...
    'string', 'zoom out', 'fontweight', 'bold', ...
    'fontsize', 16,     'callback', 'axes(ax_all); cla; imagesc(img); axis tight;');
pb_finish = uicontrol(fig_merge, 'style', 'pushbutton', ...
    'unit', 'normalized',  'position', [0.27, 0.01, 0.1, 0.08], ...
    'string', 'FINISH', 'fontweight', 'bold', ...
    'fontsize', 16,     'callback', 'neuron.delete(ind_del); close(gcf); ');

ax_selected = axes('position', [0.41, 0.6, 0.15, 0.35]);
ax_trace = axes('position', [0.6, 0.55, 0.35, 0.4]);
ax_merged = axes('position', [0.41, 0.05, 0.15, 0.35]);
ax_merged_trace = axes('position', [0.6, 0.05, 0.35, 0.4]);

%%
pb_remove = uicontrol(fig_merge, 'style', 'pushbutton', ...
    'unit', 'normalized',  'position', [.41, 0.5, 0.08, 0.08], ...
    'string', 'remove', 'fontweight', 'bold', ...
    'fontsize', 16, ...
    'callback', 'cnmfe_manual_merge_remove_neuron;');
pb_merge = uicontrol(fig_merge, 'style', 'pushbutton', ...
    'unit', 'normalized',  'position', [.49, 0.5, 0.08, 0.08], ...
    'string', 'merge', 'fontweight', 'bold', ...
    'fontsize', 16, ...
    'callback', 'cnmfe_manual_merge_merge_neurons;');
pb_reverse = uicontrol(fig_merge, 'style', 'pushbutton', ...
    'unit', 'normalized',  'position', [.41, 0.42, 0.08, 0.08], ...
    'string', 'cancel', 'fontweight', 'bold', ...
    'fontsize', 16, ...
    'callback', 'cnmfe_manual_merge_cancel;');
pb_reverse = uicontrol(fig_merge, 'style', 'pushbutton', ...
    'unit', 'normalized',  'position', [.49, 0.42, 0.08, 0.08], ...
    'string', 'NEW', 'fontweight', 'bold', ...
    'fontsize', 16, ...
    'callback', 'IDs=[]; axes(ax_selected); cla; axes(ax_trace); cla; ');

author: zhoupc


text:
function b = ParseBatch(bfile)

%ParseBatch - Parse batch file.
%
%  USAGE
%
%    b = ParseBatch(bfile)
%
%    bfile          batch file listing the parameters for each iteration
%
%  SEE
%
%    See also StartBatch, ShowBatch, InitBatch.
%

% Copyright (C) 2007-2013 by Michaël Zugaro
%
% 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.

% Check number of parameters
if nargin < 1,
	error('Incorrect number of parameters (type ''help ParseBatch'' for details).');
end

% Make sure batch file is valid
s = dbstack;
if ~isstring_FMAT(bfile) || ~exist(bfile,'file'),
	error(['Batch file not found (type ''help ' s.name ''' for details).']);
end

% Open batch file
f = fopen(bfile,'r');
if f == -1, error(['Could not open file ''' bfile '''.']); end

item = 1;
l = 1;
nFields = 0;
b.field = {};
% Read batch file
while ~feof(f),
	field = 0;
	% Read next line
	line = fgetl(f);
	while ~isempty(line),
		% Get next field
		[token,line] = strtok(line);
		% Skip rest of line if this is a comment
		if isempty(token) | token(1) == '%', break; end
		field = field + 1;
		% Determine if this is a number, vector or string
		n = str2num(token);
		if ~isempty(n),
			% It is a number, convert it to numerical format
			b.field{item,field} = n;
		elseif token(1) ~= '[',
			% It is a string, keep it as it is
			b.field{item,field} = token;
		else
			% It is a vector (or matrix)
			[token,line] = strtok([token line],'[]');
			n = str2num(token);
			if isempty(n),
				error(['Incorrect vector or matrix in batch file (line ' int2str(l) ', element ' int2str(field) ').']);
			end
			b.field{item,field} = n;
			line = line(2:end);
		end
	end
	if field > 0,
		% Make sure all lines have the same number of elements
		if item == 1,
			nFields = field;
		elseif field ~= nFields,
			if l == 2,
				error(['Incoherent numbers of elements in batch file (' int2str(nFields)  ' elements on line 1, but ' int2str(field) ' elements on line 2).']);
			else
				error(['Incoherent numbers of elements in batch file (' int2str(nFields)  ' elements on lines 1-' int2str(l-1) ', but ' int2str(field) ' elements on line ' int2str(l) ').']);
			end
		end
		% Next item
		b.line(item) = l;
		item = item + 1;
	end
	% Update line counter (for error messages)
	l = l + 1;
end

% Close file
fclose(f);

% Reset iterator
b.currentItem = 0;
b.currentField = 0;

% Empty fields
b.log = -1;
b.mfile = [];

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function out = initializeCERR(cellName)
%"initializeCERR"
%   Return a skeleton planC structure with each cell's struct array having
%   zero length and properly set fieldnames.
%
%   An optional argument allows the user to specifiy a single empty planC
%   cell to be returned rather than the full empty planC, ie 'scan',
%   'structures', 'dose' etc.
%
%Usage:
%   planInitC = initializeCERR;
%   doseInitS = initializeCERR('dose')
%   scanInitS = initializeCERR('scan')
%
%   etc.
%
%Author: J. O. Deasy, deasy@radonc.wustl.edu
%
%Rerence:
%W. Harms, Specifications for Tape/Network Format for Exchange of Treatment Planning Information,
%version 3.22., RTOG 3D QA Center, (http://rtog3dqa.wustl.edu), 1997.
%
%Last Modified:  22 Apr 02, by JOD.
%                30 Apr 03, JOD, set doseScale default to 1 as per RTOG spec. 4.00.
%                28 Jun 06, JRA, Modified function prototype to return a
%                                single variable, planInitC, rather than
%                                individual structS.  Also modified the
%                                base scanInitS to reflect planC's current
%                                structure.
%
%
% 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

if nargout > 1
    error('Invalid number of outputs requested from initializeCERR.  Check for the new single output prototype!');
end

headerInitS = struct(...
    'archive'                  ,  '',...
    'tapeStandardNumber'       ,  '',...
    'intercomparisonStandard'  ,  '',...
    'institution'              ,  '',...
    'dateCreated'              ,   '',...
    'writer'                   ,   '',...
    'sponsorID'                ,   '',...
    'protocolID'               ,   '',...
    'subjectID'                ,   '',...
    'submissionID'             ,   '',...
    'timeSaved'                ,   '',...
    'lastSavedInVer'           ,   '',...
    'CERRImportVersion'        ,   '',...
    'anonymizedID'             ,   '');

commentInitS = struct( ...
    'imageNumber'             ,   '', ...
    'imageType'               ,   '', ...
    'caseNumber'              ,   '', ...
    'patientName'             ,   '', ...
    'writer'                  ,   '', ...
    'commentDescription'      ,   '', ...
    'dateWritten'             ,   '', ...
    'fileOfOrigin'            ,   '', ...
    'unitNumber'              ,   '', ...
    'headInOut'               ,   '', ...
    'comment'                 ,   '');

scanInitS = struct(...
    'scanArray'               ,   '', ...
    'scanType'                ,   '', ...
    'scanInfo'                ,   '', ...
    'uniformScanInfo'         ,   '', ...
    'scanArraySuperior'       ,   '', ...
    'scanArrayInferior'       ,   '', ...
    'thumbnails'              ,   '', ...
    'transM'                  ,   '', ...
    'scanUID'                 ,   '', ...
    'assocDeformUID'          ,   '', ...
    'assocTextureUID'         ,   '', ...
    'assocBaseScanUID'        ,   '', ...
    'assocMovingScanUID'      ,   ''...
    );

structureInitS = struct( ...
    'imageNumber'               ,   '', ...
    'imageType'                 ,   '', ...
    'caseNumber'                ,   '', ...
    'roiNumber'                 ,   '', ...
    'patientName'               ,   '', ...
    'structureName'             ,   '', ...
    'ROIInterpretedType'        ,   '', ... % DK %% New addition to match DICOM standard
    'numberRepresentation'      ,   '', ...
    'structureFormat'           ,   '', ...
    'numberOfScans'             ,   '', ...
    'maximumNumberScans'        ,   '', ...
    'maximumPointsPerSegment'   ,   '', ...
    'maximumSegmentsPerScan'    ,   '', ...
    'structureEdition'          ,   '', ...
    'unitNumber'                ,   '', ...
    'writer'                    ,   '', ...
    'dateWritten'               ,   '', ...
    'structureColor'            ,   '', ...
    'structureDescription'      ,   '', ...
    'roiGenerationAlgorithm'    ,   '', ...
    'roiGenerationDescription'  ,   '', ...
    'studyNumberOfOrigin'       ,   '', ...
    'contour'                   ,   '', ...   %CERR addition
    'rasterSegments'            ,   '', ...   %CERR addition
    'DSHPoints'                 ,   '', ...   %CERR addition
    'orientationOfStructure'    ,   '', ...   %ORIENTATIONOFSTRUCTURE for BAT compatibility (non-RTOG keyword).
    'transferProtocol'          ,   '', ...
    'DICOMHeaders'              ,   '', ...
    'visible'                   ,   '', ...   %CERR addition
    'associatedScan'            ,   '', ...
    'strUID'                    ,   '', ...
    'assocScanUID'              ,   '', ...
    'structSetSopInstanceUID'   ,   '', ...
    'rasterized'                ,   '', ...
    'referencedFrameOfReferenceUID','', ...
    'referencedSeriesUID',''...
    );      %CERR addition

structureArrayInitS        = struct( ...
    'indicesArray'              ,   '', ...
    'bitsArray'                 ,   '', ...
    'assocScanUID'              ,   '', ...
    'structureSetUID'           ,   ''  ...
    );

beamGeometryInitS          = struct( ...
    'imageNumber'                ,   '', ...
    'imageType'                  ,   '', ...
    'caseNumber'                 ,   '', ...
    'patientName'                ,   '', ...
    'beamNumber'                 ,   '', ...
    'beamModality'               ,   '', ...
    'beamEnergyMeV'              ,   '', ...
    'beamDescription'            ,   '', ...
    'RxDosePerTxGy'              ,   '', ...
    'MonitorUnitsPerTx'          ,   '', ...   %KU added 7Oct07
    'numberOfTx'                 ,   '', ...
    'fractionGroupID'            ,   '', ...
    'beamType'                   ,   '', ...
    'planIDOfOrigin'             ,   '', ...
    'collimatorType'             ,   '', ...
    'apertureType'               ,   '', ...
    'apertureDescription'        ,   '', ...
    'collimatorAngle'            ,   '', ...
    'gantryAngle'                ,   '', ...
    'couchAngle'                 ,   '', ...
    'headInOut'                  ,   '', ...
    'nominalIsocenterDistance'   ,   '', ...
    'numberRepresentation'       ,   '', ...
    'apertureID'                 ,   '', ...
    'wedgeAngle'                 ,   '', ...
    'wedgeRotationAngle'         ,   '', ...
    'arcAngle'                   ,   '', ...
    'machineID'                  ,   '', ...
    'beamWeight'                 ,   '', ...
    'weightUnits'                ,   '', ...
    'compensator'                ,   '', ...
    'compensatorFormat'          ,   '', ...
    'file'                       ,   '', ...
    'beamUID'                    ,   ''  ...
    );       %CERR addition

beamsInitS = struct(...
    'PatientName'                     ,   '', ...
    'PatientID'                       ,   '', ...
    'PatientBirthDate'                ,   '', ...
    'PatientSex'                      ,   '', ...
    'AcquisitionGroupLength'          ,   '', ...
    'RelationshipGroupLength'         ,   '', ...
    'ImagePresentationGroupLength'    ,   '', ...
    'PixelPaddingValue'               ,   '', ...
    'PlanGroupLength'                 ,   '', ...
    'RTPlanLabel'                     ,   '', ...
    'RTPlanDate'                      ,   '', ...
    'RTPlanTime'                      ,   '', ...
    'RTPlanGeometry'                  ,   '', ...
    'TreatmentSites'                  ,   '', ...
    'PrescriptionDescription'         ,   '', ...
    'DoseReferenceSequence'           ,   '', ...
    'FractionGroupSequence'           ,   '', ...
    'BeamSequence'                    ,   '', ...
    'PatientSetupSequence'            ,   '', ...
    'ReferencedRTGroupLength'         ,   '', ...
    'ReferencedStructureSetSequence'  ,   '', ...
    'ReferencedDoseSequence'          ,   '', ...
    'ReviewGroupLength'               ,   '', ...
    'ApprovalStatus'                  ,   '', ...
    'ReviewDate'                      ,   '', ...
    'ReviewTime'                      ,   '', ...
    'ReviewerName'                    ,   '', ...
    'SOPInstanceUID'                  ,   '', ...
    'BeamUID'                         ,   '' ...
    );

digitalFilmInitS = struct( ...
    'imageNumber'                ,   '', ...
    'imageType'                  ,   '', ...
    'caseNumber'                 ,   '', ...
    'patientName'                ,   '', ...
    'beamNumber'                 ,   '', ...
    'beamDescription'            ,   '', ...
    'filmNumber'                 ,   '', ...
    'filmDate'                   ,   '', ...
    'filmType'                   ,   '', ...
    'filmSource'                 ,   '', ...
    'numberOfDimensions'         ,   '', ...
    'sizeOfDimension1'           ,   '', ...
    'sizeOfDimension2'           ,   '', ...
    'grid1Units'                 ,   '', ...
    'grid2Units'                 ,   '', ...
    'sourceImageDistance'        ,   '', ...
    'xOffset'                    ,   '', ...
    'yOffset'                    ,   '', ...
    'bytesPerPixel'              ,   '', ...
    'numberRepresentation'       ,   '', ...
    'filmDescription'            ,   '', ...
    'unitNumber'                 ,   '', ...
    'ODScale'                    ,   '', ...
    'bitsPerPixel'               ,   '', ...
    'collimatorAngle'            ,   '', ...
    'image'                      ,   '');

doseInitS  = struct( ....
    'imageNumber'                ,   '', ...
    'imageType'                  ,   '', ...
    'caseNumber'                 ,   '', ...
    'patientName'                ,   '', ...
    'doseNumber'                 ,   '', ...
    'doseType'                   ,   '', ...
    'doseSummationType'          ,   '', ...
    'refBeamNumber'              ,   '', ...
    'refFractionGroupNumber'     ,   '', ...
    'numberMultiFrameImages'     ,   '', ...
    'doseUnits'                  ,   '', ...
    'doseScale'                  ,   '', ...
    'fractionGroupID'            ,   '', ...
    'numberOfTx'                 ,   '', ...
    'orientationOfDose'          ,   '', ...
    'imagePositionPatient'       ,   '', ...
    'imageOrientationPatient'    ,   '', ...
    'numberRepresentation'       ,   '', ...
    'numberOfDimensions'         ,   '', ...
    'sizeOfDimension1'           ,   '', ...
    'sizeOfDimension2'           ,   '', ...
    'sizeOfDimension3'           ,   '', ...
    'coord1OFFirstPoint'         ,   '', ...
    'coord2OFFirstPoint'         ,   '', ...
    'horizontalGridInterval'     ,   '', ...
    'verticalGridInterval'       ,   '', ...
    'doseDescription'            ,   '', ...
    'doseEdition'                ,   '', ...
    'unitNumber'                 ,   '', ...
    'writer'                     ,   '', ...
    'dateWritten'                ,   '', ...
    'planNumberOfOrigin'         ,   '', ...
    'planEditionOfOrigin'        ,   '', ...
    'studyNumberOfOrigin'        ,   '', ...
    'studyInstanceUID'           ,   '', ...
    'versionNumberOfProgram'     ,   '', ...
    'xcoordOfNormaliznPoint'     ,   '', ...
    'ycoordOfNormaliznPoint'     ,   '', ...
    'zcoordOfNormaliznPoint'     ,   '', ...
    'doseAtNormaliznPoint'       ,   '', ...
    'doseError'                  ,   '', ...
    'coord3OfFirstPoint'         ,   '', ...
    'depthGridInterval'          ,   '', ...
    'planIDOfOrigin'             ,   '', ...
    'doseArray'                  ,   '', ...  %CERR addition
    'zValues'                    ,   '', ...  %CERR addition
    'delivered'                  ,   '', ...  %CERR addition
    'cachedColor'                ,   '', ...  %CERR addition
    'cachedTime'                 ,   '', ...  %CERR addition
    'numCachedSlices'            ,   '', ...  %CERR addition
    'transferProtocol'           ,   '', ...
    'DICOMHeaders'               ,   '', ...
    'associatedScan'             ,   '', ...  %CERR addition
    'transM'                     ,   '', ...
    'doseUID'                    ,   '', ...
    'assocScanUID'               ,   '', ...
    'assocBeamUID'               ,   '', ...
    'frameOfReferenceUID'        ,   '', ...
    'refStructSetSopInstanceUID' ,   '', ...
    'Rx'                         ,   '', ...  % Priscription Dose
    'doseOffset'                 ,   '' ...   %CERR addition
    );

DVHInitS = struct( ...
    'imageNumber'                ,   '', ...
    'imageType'                  ,   '', ...
    'caseNumber'                 ,   '', ...
    'patientName'                ,   '', ...
    'structureName'              ,   '', ...
    'doseUnits'                  ,   '', ...
    'doseType'                   ,   '', ...
    'volumeType'                 ,   '', ...
    'numberOfPairs'              ,   '', ...
    'maximumNumberPairs'         ,   '', ...
    'numberRepresentation'       ,   '', ...
    'planIDOfOrigin'             ,   '', ...
    'fractionIDOfOrigin'         ,   '', ...  %CERR addition
    'doseScale'                  ,   '', ...
    'volumeScale'                ,   '', ...
    'dateOfDVH'                  ,   '', ...
    'DVHMatrix'                  ,   '', ...  %CERR addition
    'doseIndex'                  ,   '', ...  %CERR addition, index in planC of the dose distribution used to generate this DVH
    'dvhUID'                     ,   '', ...
    'assocStrUID'                ,   '', ...
    'assocDoseUID'               ,   ''  ...
    );

IVHInitS = struct( ...
    'imageNumber'                ,   '', ...
    'imageType'                  ,   '', ...
    'caseNumber'                 ,   '', ...
    'patientName'                ,   '', ...
    'structureName'              ,   '', ...
    'scanUnits'                  ,   '', ...
    'scanType'                   ,   '', ...
    'volumeType'                 ,   '', ...
    'numberOfPairs'              ,   '', ...
    'maximumNumberPairs'         ,   '', ...
    'numberRepresentation'       ,   '', ...
    'scanScale'                  ,   '', ...
    'volumeScale'                ,   '', ...
    'dateOfIVH'                  ,   '', ...
    'IVHMatrix'                  ,   '', ...  %CERR addition
    'scanIndex'                  ,   '', ...  %CERR addition, index in planC of the dose distribution used to generate this DVH
    'ivhUID'                     ,   '', ...
    'assocStrUID'                ,   '', ...
    'assocScanUID'               ,   ''  ...
    );

planParamS = struct(...
    'doseFractionID'        ,'',...
    'structureName'         ,'',...
    'lh_goalType'           ,'',...
    'lh_goal'               ,'',...
    'lh_unit'               ,'',...
    'operator'              ,'',...
    'rh_goalType'           ,'',...
    'rh_goal'               ,'',...
    'rh_unit'               ,'',...
    'category'              ,'',...
    'failPercent'           ,'',...
    'volume'                ,'',...
    'minDose'               ,'',...
    'maxDose'               ,'',...
    'meanDose'              ,'',...
    'result'                ,'',...
    'assocDoseUID'          ,'',...
    'assocStructUID'        ,''...
    );

seedGeometryInitS = struct(...
    'imageNumber'                ,   '', ...
    'imageType'                  ,   '', ...
    'caseNumber'                 ,   '', ...
    'patientName'                ,   '', ...
    'seedModel'                  ,   '', ...
    'isotope'                    ,   '', ...
    'seedStrength'               ,   '', ...
    'strengthUnits'              ,   '', ...
    'dateOfImplant'              ,   '', ...
    'numberOfSeeds'              ,   '', ...
    'numberRepresentation'       ,   '', ...
    'planIDOfOrigin'             ,   '', ...
    'pointsM'                    ,   '' ...  %CERR addition;
    );

deformS         = struct(...
 'baseScanUID'                  ,   '', ...
 'movScanUID'                   ,   '', ...
 'algorithm'                    ,   '', ...
 'algorithmParamsS'             ,   '', ...
 'deformUID'                    ,   '',...
 'registrationTool'             ,   '',...
 'DVFfileName'                  ,   ''...
  );

registrationS         = struct(...
 'rigidS'                       ,   '', ...
 'deformS'                      ,   '', ...
 'regUID'                       ,   '' ...
  );

textureInitS           = struct(...
 'category'                   ,   '', ...
 'description'                ,   '', ...
 'patchSize'                  ,   '', ...
 'patchUnit'                  ,   '', ...
 'paramS'                     ,   '', ...
 'assocScanUID'               ,   '', ...
 'assocStructUID'             ,   '', ...
 'textureUID'                 ,   '' ...
 );

featureSetInitS           = struct(...
 'type'                       ,   '', ...
 'name'                       ,   '', ...
 'description'                ,   '', ...
 'valuesS'                    ,   '', ...
 'paramS'                     ,   '', ...
 'assocScanUID'               ,   '', ...
 'assocDoseUID'               ,   '', ...
 'assocStrUID'                ,   '', ...
 'featureSetUID'              ,   '' ...
 );

RTTreatmentInitS = struct('');          %CERR addition

importLogInitS = struct(...             %CERR addition
    'importLog'                 , '', ...
    'startTime'                 , '', ...
    'endTime'                   , '');

IMInitS = struct(...                    %CERR addition
    'IMSetup'                   , '', ...
    'IMDosimetry'               , '', ...
    'IMUID'                     , '');

GSPSInitS = struct(...
    'SOPInstanceUID'                            ,   '', ...
    'graphicAnnotationS'                        ,   '', ...
    'textAnnotationS'                           ,   '', ...
    'presentLabel'                              ,   '', ...
    'presentOrder'                              ,   '', ...
    'presentRGBColor'                           ,   '', ...
    'presentDescription'                        ,   '', ...
    'presentCreationDate'                       ,   '', ...
    'annotUID'                                  ,   '');

segmentLabelInitS = struct(...             %CERR addition
    'name'                  ,'', ...
    'labelC'                ,'', ...
    'labelColorM'           ,'', ...
    'valueS'                ,'',...
    'segLabelUID'              ,'',...
    'assocStructUID'        ,'');

planInitC = cell(1,16);

indexS.header               = 1;
indexS.comment              = 2;
indexS.scan                 = 3;
indexS.structures           = 4;
indexS.structureArray       = 5;
indexS.structureArrayMore   = 6;
indexS.beams                = 7;
indexS.beamGeometry         = 8;
indexS.dose                 = 9;
indexS.DVH                  = 10;
indexS.IVH                  = 11;
indexS.planParam            = 12;
indexS.seedGeometry         = 13;
indexS.digitalFilm          = 14;
indexS.RTTreatment          = 15;
indexS.IM                   = 16;
indexS.GSPS                 = 17;
indexS.deform               = 18;
indexS.registration         = 19;
indexS.texture              = 20;
indexS.featureSet           = 21;
indexS.importLog            = 22;
indexS.segmentLabel         = 23;
indexS.CERROptions          = 24;

% check if this is microRT and/or RPC film study and change the study
% accordingly
try
    if stateS.optS.chkMicroRT
        indexS.microRTP = 25;
        indexS.indexS   = 26;
        microRTPInitS = initializeMicroRTP(planInitC);
    else
        indexS.indexS = 25;
    end
catch
    indexS.indexS = 25;
end

indexInitS = indexS;

% Initialize the fields of each struct to be empty
headerInitS(1) = []; commentInitS(1) = []; scanInitS(1) = []; doseInitS(1) = [];
planInitC(1) = [];  DVHInitS(1) = []; IVHInitS(1) = [];
digitalFilmInitS(1) = []; importLogInitS(1) = []; IMInitS(1) = [];structureInitS(1) = [];
beamGeometryInitS(1) = []; structureArrayInitS(1) = []; seedGeometryInitS(1) = [];
deformS(1) = []; GSPSInitS(1) = []; textureInitS(1) = []; featureSetInitS(1) = []; segmentLabelInitS(1) = [];
% seedGeometryInitS(1) = []; This is never used commented by DK

planInitC{indexS.header}            = headerInitS;
planInitC{indexS.comment}           = commentInitS;
planInitC{indexS.scan}              = scanInitS;
planInitC{indexS.structures}        = structureInitS;
planInitC{indexS.structureArray}    = structureArrayInitS;
planInitC{indexS.structureArrayMore}= structureArrayInitS;
planInitC{indexS.beamGeometry}      = beamGeometryInitS;
planInitC{indexS.beams}             = beamsInitS;
planInitC{indexS.dose}              = doseInitS;
planInitC{indexS.DVH}               = DVHInitS;
planInitC{indexS.IVH}               = IVHInitS;
planInitC{indexS.planParam}         = planParamS;
planInitC{indexS.seedGeometry}      = seedGeometryInitS;
planInitC{indexS.digitalFilm}       = digitalFilmInitS;
planInitC{indexS.RTTreatment}       = RTTreatmentInitS;
planInitC{indexS.deform}            = deformS;
planInitC{indexS.importLog}         = importLogInitS;
planInitC{indexS.IM}                = IMInitS;
planInitC{indexS.GSPS}              = GSPSInitS;
planInitC{indexS.registration}      = registrationS;
planInitC{indexS.texture}           = textureInitS;
planInitC{indexS.featureSet}        = featureSetInitS;
planInitC{indexS.segmentLabel}      = segmentLabelInitS;
planInitC{indexS.CERROptions}       = struct([]); %Currently blank.  Set on planC creation.

%Store the indexS so requesting functions can extract the prototype they
%want, ie dose, structures etc.
planInitC{indexS.indexS}            = indexInitS;

try
    planInitC{indexS.microRTP}      = microRTPInitS;
end

%Handle various output options.
if nargin == 0
    out = planInitC;
elseif nargin == 1

    if isfield(indexS, cellName)
        cellNum = getfield(indexS, cellName);
        out = planInitC{cellNum};
    else
        error('Cell name requested from initializeCERR does not exist in planInitC.');
    end

    clear headerInitS commentInitS scanInitS doseInitS planInitC DVHInitS IVHInitS
    clear digitalFilmInitS importLogInitS IMInitS structureInitS beamGeometryInitS
    clear structureArrayInitS

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
classdef MimCoordinateSystem
    % MimCoordinateSystem. An enumeration used to specify the coordinate system used when
    %     importing or exporting coordinates
    %
    % Different imaging programs use different coordinate systems. Therefore,
    % if you import or export image coordinates from/to another imaging program,
    % you need to tell MIM what origin to use, otherwise the coordinates will
    % not match.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    enumeration
        PTK                % Coordinates in mm relatve to the top-left corner of the image volume
        Dicom              % Coordinates in mm relative to the scanner origin
        DicomUntranslated  % Origin in the centre of the first voxel of the image
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [ E ] = getSketchTokenEdgemap( img )
%GETSKETCHTOKENEDGEMAP Summary of this function goes here
%   Detailed explanation goes here
load('./Toolbox/SketchTokens-master/models/forest/modelSmall.mat');
st = stDetect( img, model );
E = stToEdges( st, 1 );

end

author: zouchuhang

repo: LayoutNet


text:
function [p sres sres_ns] = ResidScan(res, FWHM)
% function [p sres sres_ns] = ResidScan(res, FWHM)
%
% Calculates P(M>=t) where M is the max value of the smoothed residuals.
% In this implementation the residuals are smoothed using a Gaussian
% kernel.
%
% INPUT:
%
% res - residual time course
% FWHM - Full Width Half Maximum (in time units)
%
% OUTPUT:
%
% p - pvalues
% sres - smoothed residuals
% sres_ns - smoothed residuals (non standardized)
%
% By Martin Lindquist & Ji-Meng Loh, July 2007
%
% Edited by ML on 10/02/09

res_ns = res;
res = res./std(res);
len = length(res);

% Create Gaussian Kernel
sig = ceil(FWHM/(2*sqrt(2*log(2))));
klen = 3*sig;
kern = normpdf((-klen:klen),0,sig);
kern = kern./sqrt(sum(kern.^2));

% Convolve
x = conv(res,kern);
sres = x((klen + 1):(end-klen));

x = conv(res_ns,kern/sum(kern));
sres_ns = x((klen + 1):(end-klen));

% Find Max value
[a,location] = max(abs(sres));

% Find p-values using Gaussian Random Field theory
z = Euler_p(1, a, len, FWHM);
z = 2*z;        %Two-sided test
p = min(1, z);

end

% END MAIN FUNCTION

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Subfunctions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function pval = Euler_p(myDim, value, N, fwhm)
% function z = Euler_p(myDim, value, N, fwhm)
%
% Finds the p value using the expected Euler characteristic.
%
% This function returns P(M \ge value) using the approximation
% \sum_{d=0}^D R_d(V) \rho_d(value) following Worsley et al's "A Unified
% Statistical Approach for Determining Significant Signals in Images of
% Cerebral Activation".
%
% INPUTS:
%
% myDim - the number of dimensions in the data
% value - the value of the maximum.
% N     - the number of (time) points in that 1 dimension
% fwhm  - the full width half maximum
%
% OUTPUTS:
%
% pval - the p-value

% NOTE: CURRENTLY THIS FUNCTION IS ONLY IMPLEMENTED FOR THE 1D CASE

  % Constants
  myfactor = 4*log(2);
  pi2 = 2*pi;
  exptsq = exp(-(value^2)/2);

  % Euler Characteristc Densties
  rho = zeros(5,1);
  rho(1) = 1-normcdf(value);
  rho(2) = myfactor^(0.5)*exptsq/pi2;
  rho(3) = myfactor * exptsq * value / (pi2 ^ (1.5));
  rho(4) = myfactor ^ (1.5) * exptsq * (value^2-1) / (pi2 ^2);
  rho(5) = myfactor ^2 * exptsq * (value^3-3*value) / (pi2 ^ (5/2));

  % Resel Count
  R0 = 1;
  R1 = N/fwhm;

  % P-value
  pval = R0 * rho(1) + R1 * rho(2);

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function loc=RLSloc(X)
%Syntax: loc=RLSloc(X)
%_____________________
%
% Calculates the Reweighterd Least Squares (RLS) location parameter
% of the columns of a matrix X. If X is a vector, it returns the RLS
% location parameter of its components. If X is a scalar, it returns
% X.
%
% loc is the RLS estimated vector of locations.
% X is the matrix with the data sets.
%
% Reference:
% Rousseeuw PJ, Leroy AM (1987): Robust regression and outlier detection. Wiley.
%
%
% Alexandros Leontitsis
% Institute of Mathematics and Statistics
% University of Kent at Canterbury
% Canterbury
% Kent, CT2 7NF
% U.K.
%
% University e-mail: al10@ukc.ac.uk (until December 2002)
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% Sep 3, 2001.

if nargin<1 | isempty(X)==1
   error('Not enough input arguments.');
else
   % If X is a row vector make it a column vector
   if size(X,1)==1
      X=X';
   end
   % n is the length of the data set
   n=size(X,1);
end

% For a single data point there is no need to proceed
if n==1
   loc=X;
else
   % Find the LMS location parameter
   l=LMSloc(X);
   % Find the LMS scale parameter
   s=LMSsca(X);

   % Calculate the location parameter for every column of X
   for i=1:size(X,2);
      r=X(:,i)-l(i);
      if s(i)==0
         w=1:n;
      else
         w=find(abs(r)/s(i)<=2.5);
      end
      loc(i)=mean(X(w,i));
   end
end

author: Sable

repo: mcbench-benchmarks


text:
function AlignTsp2Whl_All(varargin)

% USAGE
%     AlignTsp2Whl_All('parameters','values')
%
% Inputs:
%    SIMPLEST THING TO DO: call AlignTsp2Whl_All()
%
%    WARNING: it is assumed that the root folder name and the final merge name are the same. If not, use
%    the 'mergename' option. Typically the root folder is something like
%    /xxx/yyy/AnimalZZZ-DayA/ and the merged data are in /xxx/yyy/AnimalZZZ-DayA/AnimalZZZ-DayA/
%
% - options:
%     'mergename': if the merged data are not in 'fbasename/fbasename' but in
%     'fbasename/mergename', then indicate with this option the name of the
%     merged data directory
%     'merge': if set to 1 (default 1) the function will merge the whl
%     dat into one single whl file 'fbasename/mergename/mergename.whl'.
%     'colorIx': a 2 value vector [colFront colRear] which defines which
%     color from [R G B] is the front and rear LEDs. Default: [1 3]*
%     'folderNames': a cell array of filenames in the order they will be
%     merged. If omitted, the program looks at the merged files in the xml
%     file contained in '/dir1/dir2/.../fbasename/mergename/mergename.xml'
%     'fromMpg': if 1 (default 0) will reconstruct the tracking from the
%     video
%
% Dependencies: xml_toolbox (easy to find on the web)
%
% Adrien Peyrache, 2012

merge=1;
colorIx=[1 3];
tspFiles = {};
tspFilesIx = [];
frommpg = 0;
[dummy fbasename dummy] = fileparts(pwd);

mergename = fbasename;

for i = 1:2:length(varargin),
  if ~isa(varargin{i},'char'),
    error(['Parameter ' num2str(i+3) ' is not a property (type ''help bz_LoadBinary'' for details).']);
  end
  switch(lower(varargin{i})),
    case 'mergename',
      mergename = varargin{i+1};
      if ~ischar(mergename)
        error('Incorrect value for property ''mergename''. Should be char array.');
      end
    case 'merge',
      merge = varargin{i+1};
      if merge ~= 1 || merge ~= 0
        error('Incorrect value for property ''merge''. Should be logical.');
      end
    case 'colorix',
      colorIx = varargin{i+1};
      if length(colorIx)~=2
        error('Incorrect value for property ''colorIx''. Should be 3 element vector of logical values.');
      end
    case 'folderNames',
      tspFiles = varargin{i+1};
      isachar = 1;
      for fix=1:length(tspFiles),isachar=isachar*isa(tspFiles{fix},'char');end
      if ~isa(tspFiles,'cell') || ~isachar
        error('Incorrect value for property ''tspFiles''. Should be a cell array of filenames');
      end
     case 'frommpg',
         frommpg = varargin{i+1};
         if ~isnumeric(frommpg)
            error('Incorrect value for property ''frommpg''. Should be 0 pr 1');
         end
     case 'tspFilesix',
         tspFilesIx = varargin{i+1};
  end
end

% if isempty(tspFiles)
%     xmldata = xml_load([fbasename '.xml']);
%     segnames = xmldata.multiFileProcessing.files;
%
%     if isempty(tspFilesIx)
%         tspFilesIx = (1:length(segnames));
%     else
%         tspFilesIx = tspFilesIx(:)';
%     end
%     tspFiles = cell(length(tspFilesIx),1);
%     for ii=1:length(tspFilesIx)
%         tspFiles{ii} = segnames(tspFilesIx(ii)).fileBaseName;
%     end
% end

tsptemp = dir('*.tsp');%get list of tsp files (to structure);
for a = 1:length(tsptemp)
    tspFiles{a}=tsptemp(a).name(1:end-4);
end

fprintf('Creating whl files...\n')
for fix=1:length(tspFiles)
    fname = tspFiles{fix};
    fprintf('\t%s\n',fname)
%     try
        AlignTsp2Whl(fname,colorIx,frommpg);
%     catch
%         keyboard
%     end
end

if merge
    fprintf('Merging files...\n')
    whlall = [];
    for fix=1:length(tspFiles)
        fname = tspFiles{fix};
        whl = load([fname '.whl']);
        whlall = [whlall;whl];
    end
    dlmwrite([mergename '.whl'],whlall,'delimiter','\t');
end


text:
function title = p04_title ( )

%*****************************************************************************80
%
%% P04_TITLE returns a title for problem 04.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = '#4: The unit hexagon with hexagonal hole.';

  return
end


text:
function [curr_segments] = filter_segments(seg_obj, gp_obj, graph_idx, ...
                                           curr_segments)
%FILTER_SEGMENTS filters (discards) the segments based on different
% criterion
%
% @authors:     Ahmad Humayun,  Fuxin Li
% @contact:     ahumayun@cc.gatech.edu
% @affiliation: Georgia Institute of Technology
% @date:        Fall 2013 - Summer 2014

    % first check if we want to filter segments for the current graph
    if ~seg_obj.segm_params.graph_filter_segs{graph_idx}
        curr_segments.segs_meta_info.energies = [];
        return;
    end

    fprintf('\tFiltering segments\n');

    t_fltr = tic;
    % if no segments, return
    if isempty(curr_segments.cut_segs)
        return;
    end

    % separate each segment into its connected components. Also remove the
    % components (split segments) which are below min_npixels size
    [curr_segments] = ...
        separate_conn_comp(curr_segments, seg_obj, ...
                           seg_obj.segm_params.filter_min_seg_pixels);

    % just computes the pairwise edges cut value divided by the number of
    % pairwise edges cut (this is called the cut ratio)
    [curr_segments, t_energy] = compute_energies(curr_segments, ...
                                                 seg_obj, gp_obj);

    % filter segment which are above the filter_max_energy cut ratio
    [curr_segments] = filter_energies(curr_segments, seg_obj, ...
                                seg_obj.segm_params.filter_max_energy, ...
                                t_energy);

    % randomly throw away segments if more than filter_max_rand
    [curr_segments] = filter_rand(curr_segments, seg_obj, ...
                                  seg_obj.segm_params.filter_max_rand);

    % if multiple segments overlap, only keep the one with the lowest cut
    % ratio
    [curr_segments] = remove_repeated_segments(curr_segments, seg_obj, ...
                                               gp_obj);

    % if user specified specified seeds, discard anything that doesn't
    % include the seed location itself (which can happen with color unaries
    % for instance)
    if strcmp(seg_obj.segm_params.graph_seed_gen_method,'gen_user_seeds')
        % find segments which cover all the seed superpixels (supposes that
        % there is only one seed in seed_sets)
        seed_cover = ...
            curr_segments.cut_segs(seg_obj.precomputed_seeds.seed_sets,:);
        conn = (sum(seed_cover,1) ...
            == sum(seg_obj.precomputed_seeds.seed_sets));
        [curr_segments] = cherry_pick_segments(curr_segments, conn);
    end

    seg_obj.timings.seg_filtering_time = ...
        [seg_obj.timings.seg_filtering_time, toc(t_fltr)];
end

function [curr_segments] = separate_conn_comp(curr_segments, seg_obj, ...
                                              min_npixels)
    fprintf('\t\tSeparating connected components ... ');

    t_conn = tic;

    if(nargin == 1)
        min_npixels = 25;
    end

    % separate each segment into its connected components by
    % finding connected components on the induced graph. Also remove the
    % components which are below min_npixels size
    [segms, num_comps] = sp_conncomp_mex(curr_segments.cut_segs, ...
                                         seg_obj.sp_data.edgelet_sp, ...
                                         seg_obj.sp_data.sp_seg_szs, ...
                                         min_npixels);

    % make the mapping which is an array of the length of the total
    % number of new segments (each of which forms a connected
    % component) where each array location gives the index/id of
    % the original segment returned by the min-cut
    temp = duplicateElems(1:length(num_comps), num_comps);
    curr_segments.segs_meta_info.seg_mapping_final_to_orig = ...
        curr_segments.segs_meta_info.seg_mapping_final_to_orig(temp);

    curr_segments.cut_segs = segms;

    seg_obj.num_segs.after_splitting_conncomp = ...
        [seg_obj.num_segs.after_splitting_conncomp, sum(num_comps)];

    time_util(seg_obj, 'init_filter_time', t_conn, 1, 1);
end

function [curr_segments, t_energy] = compute_energies(curr_segments, ...
                                                      seg_obj, gp_obj)
    fprintf('\t\tComputing energies ... ');

    t_energy = tic;

    cut_ratio = zeros(1,size(curr_segments.cut_segs,2));
    for i = 1:length(cut_ratio)
        links_across = ...
            gp_obj.pairwise_graph(~curr_segments.cut_segs(:,i), ...
                                  curr_segments.cut_segs(:,i));
        cut = full(sum(links_across(:)));
        % cut ratio (minimizing it is called the sparsest cut problem)
        n_edges_across = nnz(links_across);
        cut_ratio(i) = cut / n_edges_across;
    end
    curr_segments.segs_meta_info.energies = cut_ratio;

    fprintf('%.2fs\n', toc(t_energy));
end

function [curr_segments] = filter_energies(curr_segments, seg_obj, ...
                                           max_energy, t_energy)
    fprintf('\t\tRemoving high energy solutions ... ');

    t_energy_fltr = tic;

    % minimum number of segments to keep
    min_n_segms = 5;

    [sorted_cut_ratio, sorted_ind] = ...
        sort(curr_segments.segs_meta_info.energies, 'ascend');
    last_acceptable = find(sorted_cut_ratio <= max_energy, 1, 'last');
    reject_segs_ind = sorted_ind(last_acceptable+1:end);
    % if the remaining number of segments would be less
    num_remain_segs = length(sorted_cut_ratio) - length(reject_segs_ind);
    if num_remain_segs < min_n_segms
        num_keep_more = min_n_segms - num_remain_segs;
        num_keep_more = min(num_keep_more, length(reject_segs_ind));
        reject_segs_ind(1:num_keep_more) = [];
    end

    % remove segments not wanted from curr_segments
    [curr_segments] = cherry_pick_segments(curr_segments, ...
                                           setdiff(1:length(sorted_ind), ...
                                                   reject_segs_ind));

    seg_obj.num_segs.after_energy_filtering = ...
        [seg_obj.num_segs.after_energy_filtering, ...
         size(curr_segments.cut_segs,2)];

    fprintf('%.2fs\n', toc(t_energy_fltr));

    time_util(seg_obj, 'energy_filter_time', t_energy, 1, 0);
end

function [curr_segments] = filter_rand(curr_segments, seg_obj, ...
                                       max_rand_pick)
    fprintf('\t\tPick %d random segments ... ', max_rand_pick);

    t_rand = tic;

    if size(curr_segments.cut_segs,2) > max_rand_pick
        randn = randperm(size(curr_segments.cut_segs,2));
        curr_segments = cherry_pick_segments(curr_segments, ...
                                             randn(1:max_rand_pick));
    end

    seg_obj.num_segs.after_random_picking = ...
        [seg_obj.num_segs.after_random_picking, ...
         size(curr_segments.cut_segs,2)];

    time_util(seg_obj, 'rand_filter_time', t_rand, 1, 1);
end

function [curr_segments] = remove_repeated_segments(curr_segments, ...
                                                    seg_obj, gp_obj)
    fprintf('\t\tRemoving repeated segments ... ');

    OVERLAP_THRESH = 0.95;

    t_repeat = tic;

    seg_sel = [];

    % divide segments by their graph sub-methods
    ub = cumsum(gp_obj.graph_sets_per_method);
    lb = [1, ub(1:end-1)+1];
%     % if you don't want division by graph sub-methods (comment out above)
%     ub = inf;
%     lb = 1;

    % find sols_to_unary_mapping for the segments remaining
    sols_to_unary_map = curr_segments.segs_meta_info.sols_to_unary_mapping;
    seg_mapping = curr_segments.segs_meta_info.seg_mapping_final_to_orig;
    curr_unary_mapping = sols_to_unary_map(seg_mapping);
    for gsubm = 1:length(ub)
        subm_idx = curr_unary_mapping >= lb(gsubm) & ...
                   curr_unary_mapping <= ub(gsubm);
        cut_segs = curr_segments.cut_segs(:,subm_idx);

        % smart way to compute pairwise overlaps btw segments if above 0.95
        overlap_mat = overlap_over_threshold(cut_segs, OVERLAP_THRESH);

        % find segment sets which are stringed together with high overlap
        bw_mat = overlap_mat >= OVERLAP_THRESH;
        [num_conn, conncomps_t] = graphconncomp(sparse(bw_mat));
        conncomps = zeros(size(subm_idx));
        conncomps(subm_idx) = conncomps_t;
        this_sel = zeros(1,num_conn);
        % for each segment set, select one segment with lowest energy
        for i = 1:num_conn
            s1 = find(conncomps == i);
            [~,b] = min(curr_segments.segs_meta_info.energies(conncomps == i));
            this_sel(i) = s1(b);
        end
        % collate to the list of segments to be cherry picked
        seg_sel = [seg_sel, this_sel];
    end

    % keep only segments which have lowest energy amongst similar segments
    [curr_segments] = cherry_pick_segments(curr_segments, seg_sel);

    seg_obj.num_segs.after_repeat_remove = ...
        [seg_obj.num_segs.after_repeat_remove, length(seg_sel)];

    time_util(seg_obj, 'seg_similar_filter_time', t_repeat, 1, 1);
end

function [curr_segments] = cherry_pick_segments(curr_segments, to_pick_idxs)
% this function removes segments from curr_segments structure

    curr_segments.cut_segs = curr_segments.cut_segs(:, to_pick_idxs);
    curr_segments.segs_meta_info.seg_mapping_final_to_orig = ...
        curr_segments.segs_meta_info.seg_mapping_final_to_orig(:, to_pick_idxs);

    if isfield(curr_segments.segs_meta_info, 'energies') && ...
            ~isempty(curr_segments.segs_meta_info.energies)
        curr_segments.segs_meta_info.energies = ...
            curr_segments.segs_meta_info.energies(to_pick_idxs);
    end
end

author: Cloud-CV

repo: object-proposals


text:
%% Rotational Approximation and Interpolation
%
%%
% On this page, we want to cover the topic of function approximation from
% discrete values on the Rotation group. To simulate this, we have stored some
% nodes and corresponding function values which we can load. The csv-file
% contains the Euler angles $\phi_1$, $\Phi$ and $\phi_2$ of the nodes and the function
% value in the fourth column. Lets import these data using the function
% <orientation.load.html |load|>

fname = fullfile(mtexDataPath, 'orientation', 'dubna.csv');
[nodes, S] = orientation.load(fname,'columnNames',{'phi1','Phi','phi2','values'});

%%
% The second output |S| is a struct that contains a field |S.values| with
% the function values from the fourth column. Next, we can make a section
% plot to see, what we are dealing with

plotSection(nodes, S.values,'all');

%%
% Now, we want to find a function which coincides with the given function
% values in the nodes reasonably well.

%% Interpolation
%
%%
% Interpolation is done by the <SO3Fun.interpolate |interpolate|> command
% of class <SO3Fun.SO3Fun |SO3Fun|>

SO3F = SO3Fun.interpolate(nodes, S.values,'exact');
plot(SO3F)

%%
% The interpolation is done by lsqr. Hence the error is not in machine
% precision.
norm(SO3F.eval(nodes) - S.values)

%%
% If we don't restrict ourselfs to the given function values in the nodes, we have more
% freedom, which can be seen in the case of approximation.

%% Approximation
%
% In contrast to interpolation we are now not restricted to the function
% values in the nodes but still want to keep the error reasonably small.
%
%%
% One way is to interpolate the function similary as before, without the
% option |'exact'|.
%
%%
% Another way is to approximate the rotational function with a series of
% <WignerFunctions.html Wigner-D functions> (Harmonic series).
% We don't take as many Wigner-D functions as there are nodes,
% such that we are in the overdetermined case. In that way we don't have a
% chance of getting the error in the nodes zero but hope for a smoother
% approximation. This can be achieved by the <SO3FunHarmonic.approximation |approximation|>
% command of the class <SO3FunHarmonic.SO3FunHarmonic |SO3FunHarmonic|>

SO3F2 = SO3FunHarmonic.approximation(nodes, S.values);
plot(SO3F2)

%%
% Plotting this function, we can immidiately see, that we have a much
% smoother function. But one has to keep in mind that the error in the data
% nodes is not zero as in the case of interpolation.

norm(eval(SO3F, nodes) - S.values)

%%
% But this may not be of great importance like in the case of function
% approximation from noisy function values, where we don't know the exact
% function values anyways.

%%
%
% The strategy underlying the |approximation|-command
% to obtain such an approximation works via Wigner-D functions
% (<SO3FunHarmonicSeries Basics of rotational harmonics>). For that,
% we seek for so-called Fourier-coefficients ${\bf \hat f} = (\hat
% f^{0,0}_0,\dots,\hat f^{N,N}_N)^T$ such that
%
% $$ g(x) = \sum_{n=0}^N\sum_{k,l = -n}^n \hat f_n^{k,l} D_n^{k,l}(x) $$
%
% approximates our function. A basic strategy to achieve this is through
% least squares, where we minimize the functional
%
% $$ \sum_{m=1}^M|f(x_m)-g(x_m)|^2 $$
%
% for the data nodes $x_m$, $m=1,\dots,M$, $f(x_m)$ the target function
% values and $g(x_m)$ our approximation evaluated in the given data nodes.
%
% This can be done by the |lsqr| function of Matlab, which efficiently
% seeks for roots of the derivative of the given functional (also known as
% normal equation). In the process we compute the matrix-vector product
% with the Fourier-matrix multible times, where the Fourier-matrix is given
% by
%
% $$ F = [D_n^{k,l}(x_m)]_{m = 1,\dots,M;~n = 0,\dots,N,\,k,l = -n,\dots,n}. $$
%
% This matrix-vector product can be computed efficiently with the use of
% the nonequispaced SO(3) Fourier transform
% <https://www-user.tu-chemnitz.de/~potts/nfft/nfsoft.php NSOFT>
% or faster by the combination of an Wigner-transform together with a
% <https://www-user.tu-chemnitz.de/~potts/nfft/index.php NFFT>.
%
% We end up with the Fourier-coefficients of our approximation $g$, which
% describe our approximation.
%

author: mtex-toolbox

repo: mtex

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


text:
function fem1d_bvp_linear_test ( )

%*****************************************************************************80
%
%% FEM1D_BVP_LINEAR_TEST tests the FEM1D_BVP_LINEAR library.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/fem1d_bvp_linear/fem1d_bvp_linear_test.m
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 June 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_BVP_LINEAR_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the FEM1D_BVP_LINEAR library.\n' );

  fem1d_bvp_linear_test01 ( );
  fem1d_bvp_linear_test02 ( );
  fem1d_bvp_linear_test03 ( );
  fem1d_bvp_linear_test04 ( );
  fem1d_bvp_linear_test05 ( );
  fem1d_bvp_linear_test06 ( );
  fem1d_bvp_linear_test07 ( );
  fem1d_bvp_linear_test08 ( );
  fem1d_bvp_linear_test09 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_BVP_LINEAR_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function matrixIndices = blkIndices(blockIndex, blockSize)
% The function returns the indices corresponding to the blockIndex-th block
% for a matrix or vector with blocks of size blockIndex X blockIndex
%
% Luca Carlone
% Georgia Institute of Technology
% 1/5/2013

numBlocks = length(blockIndex);
matrixIndices = zeros(numBlocks*blockSize,1);

for i=1:numBlocks
  matrixIndices(blockSize*i-(blockSize-1) : blockSize*i) = [blockSize*blockIndex(i)-(blockSize-1) : blockSize*blockIndex(i)];
end

end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled
% 	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:  nickabattista[at]gmail.com
% IB2d Created: May 27th, 2015
% Institution: TCNJ
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs or non-invariant beams*)
% 	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 (nickabattista[at]gmail.co) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: updates the beam attributes!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function beams_info = update_nonInv_Beams(dt,current_time,beams_info)

% beams_info:   col 1: 1ST PT.
%               col 2: MIDDLE PT. (where force is exerted)
%               col 3: 3RD PT.
%               col 4: beam stiffness
%               col 5: x-curvature
%               col 6: y-curvature

%IDs = beams_info(:,1);   % Gives Middle Pt.

% Time initialization
t1 = 0.25;                          % Period A->B (down-stroke)
t2 = 0.25;                          % Period B->A (up-stroke)
tR = 0.0;                           % Resting time between A->B and B->A
period = (t1+t2+tR);                % Time it takes to move to the right
t = rem(current_time,period);       % Recalculate Time (using modular arithmetic)

% Cubic Interpolation Information
p1 = 0.10;
p2 = 0.90;

% a COEFFICIENTS
a0 = 0;
a1 = 0;
a2 = 0;
a3 = 11.111111111111088;

% b COEFFICIENTS
b0 = 0.013888888888889;
b1 = -0.416666666666666;
b2 = 4.166666666666658;
b3 = -2.777777777777770;

% c COEFFICIENTS
c0 =  -10.111111111111073;
c1 =  33.333333333333222;
c2 = -33.333333333333222;
c3 =  11.111111111111073;

% Read In y_Pts for two Phases!
[xP1,yP1,yP2] = read_File_In('swimmer.phases'); % NOTE xP1 = xP2
xP2 = xP1;

%
% FIRST WE COMPUTE THE INTERPOLATE GEOMETRY BETWEEN BOTH PHASES
%

%
% START THE INTERPOLATING BETWEEN STATES!
%
if t <= t1 % STATE A -> STATE B

    % Scaling time for appropriate use in interp. function so tTilde\in[0,1]
    tTilde = (t/t1);

    % Evaluate Pieceise Cubic Interpolation Poly
    if ( tTilde<=p1 )
        gFUNC = a0 + a1*tTilde + a2*tTilde^2 + a3*tTilde^3;
    elseif ( (tTilde>p1) && (tTilde<=p2) )
        gFUNC = b0 + b1*tTilde + b2*tTilde^2 + b3*tTilde^3;
    else
        gFUNC = c0 + c1*tTilde + c2*tTilde^2 + c3*tTilde^3;
    end

    %xPts = xP1 + gFUNC*( xP2 - xP1 );
    yPts = yP1 + gFUNC*( yP2 - yP1 );

elseif ( t >= t1+tR ) % STATE B -> A

    % Scaling time for appropriate use in interp. function so tTilde\in[0,1]
    tTilde = (t-t1-tR)/(t2);

    % Evaluate Pieceise Cubic Interpolation Poly
    if tTilde<=p1
        gFUNC = a0 + a1*tTilde + a2*tTilde^2 + a3*tTilde^3;
    elseif tTilde<=p2
        gFUNC = b0 + b1*tTilde + b2*tTilde^2 + b3*tTilde^3;
    else
        gFUNC = c0 + c1*tTilde + c2*tTilde^2 + c3*tTilde^3;
    end

    %xPts = xP2 + gFUNC*( xP1 - xP2 );
    yPts = yP2 + gFUNC*( yP1 - yP2 );

end

%
% NOW WE UPDATE THE CURAVTURES APPROPRIATELY
%
%beams_info(:,5) = xPts(1:end-2)+xPts(3:end)-2*xPts(2:end-1);
beams_info(:,6) = yPts(1:end-2)+yPts(3:end)-2*yPts(2:end-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: reads in info from file
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x1,y1,y2] = read_File_In(file_name)

filename = file_name;  %Name of file to read in

fileID = fopen(filename);

    % Read in the file, use 'CollectOutput' to gather all similar data together
    % and 'CommentStyle' to to end and be able to skip lines in file.
    C = textscan(fileID,'%f %f %f','CollectOutput',1);

fclose(fileID);        %Close the data file.

mat_info = C{1};   %Stores all read in data

%Store all elements in matrix
mat = mat_info(1:end,1:end);

x1 =  mat(:,1); %store xVals1/2
y1 =  mat(:,2); %store yVals1
y2 =  mat(:,3); %store yVals2

author: nickabattista


text:
function varargout = ParseArgs(Args, varargin)

% Check input arguments
assert(nargin >= 2, ...
    'Function takes at least two arguments.');
assert(nargout == nargin-1, ...
    'Number of output arguments must be equal with the number of input arguments minus 1.');
assert(iscell(Args) && (isempty(Args) || isvector(Args)), ...
    'Args must be a cell vector.');
for i = 1:length(varargin)
    assert(iscellstr(varargin{i}) && isvector(varargin{i}), ...
        'All input arguments except the first one must be cell vectors of strings.');
end

% Convert Args to structure
if numel(Args) == 1
    if isstruct(Args{1})
        ArgStruct = Args{1};
    else
        error('If Args has only one element, that element must be a structure.');
    end
else
    ArgStruct = NameValue2Struct(Args{:});
end

% Make a cell array with all names in varargin
NumNames = sum(cellfun(@length,varargin));
AllNames = cell([1 NumNames]);
k = 0;
for i = 1:length(varargin)
    for j = 1:length(varargin{i})
        k = k + 1;
        AllNames{k} = varargin{i}{j};
    end
end
assert(k == NumNames);

% Check if Args contains only arguments with names that are not among those in varargin
ArgNames = fieldnames(ArgStruct);
for i = 1:length(ArgNames)
    if ~any(strcmp(ArgNames{i},AllNames))
        error('Unknown parameter: ''%s''',ArgNames{i});
    end
end

% Create output structures
for i = 1:nargout
    ParamStruct = struct();
    for j = 1:length(varargin{i})
        ParamName = varargin{i}{j};
        if isfield(ArgStruct,ParamName)
            ParamStruct.(ParamName) = ArgStruct.(ParamName);
        else
            ParamStruct.(ParamName) = [];
        end
    end
    varargout{i} = ParamStruct;
end

    function Struct = NameValue2Struct(varargin)
        for n = 2:2:nargin
            if iscell(varargin{n})
                % Put cell arrays in another cell array so that we don't get a
                % structure vector after the conversion
                varargin{n} = {varargin{n}};
            end
        end
        Struct = struct(varargin{:});
    end

end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
%% axisLim
% Below is a demonstration of the features of the |axisLim| function

%%
clear; close all; clc;

%% Syntax
% |axLim=axisLim(V);|

%% Description
% This function computes appropiate axis limits for the input vertices V.
% The vertices may be an k x l x m array, where by by k is the number of
% vertices, l is the number of dimensions (e.g. 2 or 3), and m is for
% instance a time (or other) dimension. The function returs axLim which are
% appropriate axis limits such that the coordinates in V can be displayed
% appropriately (and in a tight fashion).

%% Examples

%% Example 1: setting axis limits tightly around a coordinate set in 2D

%%
% Create example data
t=linspace(0,2*pi,250)';
r=5+2*sin(5*t);
V=r.*[cos(t) sin(t)];

%%
% Compute appropriate limits
axLim=axisLim(V) %Axis limits for vertices

%%
% Assign axis limits using |axisLim|

cFigure;
plot(V(:,1),V(:,2),'b-');
axis equal; axis(axisLim(V));
drawnow;

%% Example 2: setting axis limits tightly around a coordinate set in 3D

%%
% Create example data
[F,V]=stanford_bunny;

%%
% Compute appropriate limits
axLim=axisLim(V) %Axis limits for vertices

%%
% Assign axis limits using |axisLim|

cFigure;
gpatch(F,V);
axisGeom; camlight headlight;
axis(axisLim(V));
gdrawnow;

%% Example 3: setting axis limits tightly around multiple coordinate sets in 3D

%%
% Create example data
[F1,V1]=graphicsModels(3);
[F2,V2]=graphicsModels(4);
V2=V2+1;
%%
% Compute appropriate limits
axLim=axisLim(V1,V2) %Axis limits for vertices

%%
% Assign axis limits using |axisLim|

cFigure; hold on;
gpatch(F1,V1,'gw');
gpatch(F2,V2,'rw');
axisGeom; camlight headlight;
axis(axisLim(V1,V2));
gdrawnow;

%% Example 3: setting axis limits tightly around multiple coordinate sets with varying dimensions

%%
% Create example data
V1=r.*[cos(t) sin(t)];
[F2,V2]=geoSphere(2,1);

%%
% Compute appropriate limits
axLim=axisLim(V1,V2) %Axis limits for vertices

%%
% Assign axis limits using |axisLim|

cFigure; hold on;
plot(V1(:,1),V1(:,2),'b-');
gpatch(F2,V2,'rw');
axisGeom; camlight headlight;
axis(axisLim(V1,V2));
gdrawnow;

%% Example 4: setting axis limits tightly for coordinates over time

[F,V]=stanford_bunny; %Some graphics data

nSteps=50;
t=linspace(0,pi,nSteps)';
U=[-100*t/pi zeros(nSteps,1) 50*sin(t)];

V_DEF=zeros(size(V,1),3,nSteps);
for q=1:1:nSteps
    V_DEF(:,:,q)=V+U(q*ones(size(V,1),1),:);
end

%%
% Compute appropriate limits
axLim=axisLim(V_DEF) %Axis limits for vertices

%%
% Assign axis limits using |axisLim|

hf=cFigure;
hp=gpatch(F,V,'gw');
axisGeom; camlight headlight;
axis(axisLim(V_DEF)); %Set limits to be suitable across time
drawnow;

%%
% Use |anim8| to animate the scene

animStruct.Time=linspace(0,1,nSteps); %Time vector

for q=1:1:nSteps
    %Set entries in animation structure
    animStruct.Handles{q}=hp; %Handles of objects to animate
    animStruct.Props{q}={'Vertices'}; %Properties of objects to animate
    animStruct.Set{q}={V_DEF(:,:,q)}; %Property values for to set in order to animate
end

% Start |anim8| gui
anim8(hf,animStruct);

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON