text:
classdef MrUnitTest < matlab.unittest.TestCase
    % Implements unit testing for all MrClasses
    %
    % EXAMPLE
    %   MrUnitTest
    %
    %   See also

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

    properties (TestParameter)
        % MrDimInfo
        testVariantsDimInfo = {'1', '2', '3', '4', '5', '6', '7', '8', '9', ...
            '10', '11'};
        testFile = {'3DNifti', '4DNifti', 'Folder', 'ParRec'};
        testVariantsDimInfoSplit = {'singleDim', 'twoDims', 'nonExistingDim', ...
            'charSplitDim', 'differentIndex', 'cellofChars', 'emptySplitDim'};
        testVariantsDimInfoSelect = {'singleDim', 'multipleDimsWithSelection',...
            'type', 'invert', 'removeDims', 'unusedVarargin'};
        testCaseLoadMat = {'checkTempDir', 'oneVar', 'objectAsStruct', ...
            'className', 'noMatch', 'tooManyMatch', 'withVarName'};
        % MrAffineTransformation
        testVariantsAffineTrafo = {'propVal', 'matrix'};
        testFileAffineTrafo = {'3DNifti', '4DNifti', 'ParRec'};
        % MrImageGeometry
        testVariantsImageGeom = {'makeReference', ...
            'dimInfoAndaffineTransformation', 'matrix', 'dimInfo', ...
            'affineTransformation', 'timing_info'};
        % MrDataNd
        testVariantsDataNd = {'matrix', 'matrixWithDimInfo', 'matrixWithPropVal'};
        testVariantsArithmeticOperation = {'minus', 'plus', 'power', 'rdivide', 'times'};
        testVariantsDimensionOperation = {'circshift', 'flip', 'fliplr', 'flipud', ...
            'resize', 'rot90', 'select', 'split'};
        testVariantsValueOperation = {'cumsum', 'diff', 'fft', 'hist', 'ifft', ...
            'isreal', 'max', 'maxip', 'mean', 'power', 'prctile', 'real', ...
            'rms', 'rmse', 'unwrap'};
        testVariantsShiftTimeseries = {0}; % verbosity level (plots)
        testVariantsSelect = {'multipleDims', 'invert', 'removeDims', 'unusedVarargin'};
        % MrImage
        MrImageLoadConditions = {'4DNifti', 'FilePlusDimLabelsUnits', ...
            'FilePlusResolutions', 'FilePlussamplingWidths', ...
            'FilePlusSamplingPoints', 'FilePlusShearRotation', ...
            'FilePlusSelect', 'FilePlusDimInfoPropVals', ...
            'FilePlusAffineTransformation', 'FilePlusFirstSamplingPoint', ...
            'FilePlusLastSamplingPoint', 'FilePlusArrayIndex', ...
            'FilePlusOriginIndex'};
    end
    %% MrDimInfo
    methods (Test, TestTags = {'Constructor', 'MrDimInfo'})
        this = MrDimInfo_constructor(this, testVariantsDimInfo)
        this = MrDimInfo_load_from_file(this, testFile)
        this = MrDimInfo_load_from_mat(this, testCaseLoadMat)
    end

    methods (Test, TestTags = {'Methods', 'MrDimInfo'})
        this = MrDimInfo_get_add_remove(this)
        this = MrDimInfo_index2sample(this)
        this = MrDimInfo_permute(this)
        this = MrDimInfo_split(this, testVariantsDimInfoSplit)
        this = MrDimInfo_select(this, testVariantsDimInfoSelect)
        this = MrDimInfo_update_and_validate_properties_from(this)
    end

    %% MrAffineTransformation
    methods (Test, TestTags = {'Constructor', 'MrAffineTransformation'})
        this = MrAffineTransformation_constructor(this, testVariantsAffineTrafo)
        this = MrAffineTransformation_load_from_file(this, testFileAffineTrafo)
    end

    methods (Test, TestTags = {'Methods', 'MrAffineTransformation'})
        this = MrAffineTransformation_transformation(this)
        this = MrAffineTransformation_affineMatrix(this)
    end
    %% MrImageGeometry
    methods (Test, TestTags = {'Constructor', 'MrImageGeometry'})
        this = MrImageGeometry_constructor(this, testVariantsImageGeom)
        this = MrImageGeometry_load_from_file(this, testFile)
    end

    methods (Test, TestTags = {'Methods', 'MrImageGeometry'})
        this = MrImageGeometry_create_empty_image(this)
    end
    %% MrDataNd
    methods (Test, TestTags = {'Constructor', 'MrDataNd'})
        this = MrDataNd_constructor(this, testVariantsDataNd)
        % loading of nifti data will be tested in MrImage (since there,
        % also the affineTransformation is created)
    end

    methods (Test, TestTags = {'Methods', 'MrDataNd'})
        this = MrDataNd_arithmetic_operation(this, testVariantsArithmeticOperation);
        this = MrDataNd_permute(this);
        this = MrDataNd_select(this, testVariantsSelect);
        % this = MrDataNd_dimension_operation(this, testDimensionOperation);
        this = MrDataNd_value_operation(this, testVariantsValueOperation);
        this = MrDataNd_shift_timeseries(this, testVariantsShiftTimeseries);
        % these take too long, need to be shortened
        % this = MrDataNd_split_epoch(this);
    end

    %% MrImage
    methods (Test, TestTags = {'Constructor', 'MrImage'})
        this = MrImage_load_from_file(this, MrImageLoadConditions)
    end

end

author: translationalneuromodeling

repo: tapas


text:
function Q_gpfa = test_gpfa(seed)

if nargin >= 1
  randn('state', seed);
  rand('state', seed);
end

%
% Generate some data
%

D = 3;

in_x = 1:200;
N_x = length(in_x);
X = zeros(D,N_x);
X(1,:) = in_x/N_x;%sin(2*pi*in_w/10);
X(2,:) = cos(2*pi*in_x/7);
X(3,:) = randn(N_x,1);

in_w = 1:100;
N_w = length(in_w);
W = zeros(D,N_w);
W(1,:) = cos(2*pi*in_w/20);
W(2,:) = cos(2*pi*in_w/10);
W(3,:) = randn(N_w,1);

s = 0.5;
Y_noiseless = W'*X;
Y = Y_noiseless + s*randn(N_w,N_x);

Imv = (rand(size(Y)) < 0.4);
Imv(:,15) = true;
Y(Imv) = NaN;

%
% GP module for X
%

N_p = 10;
pseudo_x = linspace(min(in_x),max(in_x), N_p);
is_pseudo_x = false(D,1);

D2_xx = sq_dist(in_x);
D2_pp = sq_dist(pseudo_x);
D2_px = sq_dist(pseudo_x, in_x);

% Covariance functions
covfunc = @(D2) gp_cov_se(D2);
covfunc_x = cell(D,1);
% $$$ covfunc_x{1} = gp_cov_jitter(gp_cov_se(D2_xx));
covfunc_x{1} = gp_cov_pseudo(gp_cov_jitter(gp_cov_se(D2_pp)), ...
                             gp_cov_se(D2_px), ...
                             gp_cov_se(diag(D2_xx)));
is_pseudo_x(1) = true;
theta_x{1} = 30;

covfunc_x{2} = gp_cov_jitter(covfunc(D2_xx));
theta_x{2} = 3;

covfunc_x{3} = gp_cov_jitter(gp_cov_pp(sqrt(D2_xx),1));
theta_x{3} = 1.1;

% $$$ figure
% $$$ imagesc(covfunc_x{1}(theta_x{1}));
% $$$ figure
% $$$ imagesc(covfunc_x{2}(theta_x{2}));
% $$$ return

%
% GP module for W
%

D2_ww = sq_dist(in_w);

% Covariance functions
covfunc_w = cell(D,1);
covfunc_w{1} = gp_cov_scale(gp_cov_jitter(gp_cov_se(D2_ww)));
theta_w{1} = [1; 3];
covfunc_w{2} = gp_cov_scale(gp_cov_jitter(gp_cov_se(D2_ww)));
theta_w{2} = [1; 2];
covfunc_w{3} = gp_cov_scale(gp_cov_jitter(gp_cov_pp(sqrt(D2_ww),1)));
theta_w{3} = [1; 1.1];

%
% Isotropic noise module
%

%noise_module = noise_module_fixed(1/s^2 * ones(N_w, N_x));
noise_module = noise_module_isotropic(N_w, N_x, 1e-3, 1e-3, 'init', 100);

%
% VB inference
%

%
% TODO:
%
% - pseudo inputs
%
% - put the noise to Q(W)
%
% - test more components
%
% - rotation, optimize the hyperparameters jointly?
%
% - weighted noise
%

%
% GPFA
%

% GP modules with component-wise factorization

% NOTES:
%
% The algorithm can be quite sensitive to hyperparameter initialization and
% update schedule. Also, it seems that it's best to do rotations rarely and
% after the hyperparameters have been updated at least once. The rotations
% in GPFA are approximately optimized, so some problems can be caused by
% that. So it might be a good idea NOT to use rotations..

% Update only the smooth components at the beginning
maxiter = 30;
update_schedule = cell(maxiter,1);
update_schedule(:) = {1:D};
update_schedule(1:10) = {1};
%update_schedule(6:10) = {1:2};

X_module = factor_module_gp_factorized(N_x, covfunc_x, theta_x, ...
                                       'is_pseudo', is_pseudo_x, ...
                                       'update_hyperparameters', [5:5:1000], ...
                                       'update_schedule', {update_schedule});
W_module = factor_module_gp_factorized(N_w, covfunc_w, theta_w, ...
                                       'update_hyperparameters', [5:5:1000], ...
                                       'update_schedule', {update_schedule});

Q_gpfa = vbfa(D, Y, W_module, X_module, noise_module, ...
              'maxiter', maxiter, ...
              'update_noise', 1, ...
              'rotate', false, ...
              'rotation_checkgrad', false);

figure
plot(Q_gpfa.loglikelihood)

%
% PCA
%

X_module = factor_module_iid(D, N_x);
W_module = factor_module_ard(D, N_w);

Q_pca = vbfa(D, Y, W_module, X_module, noise_module, ...
         'maxiter', 30, ...
         'update_noise', 1, ...
         'rotate', 1);

%
% Results
%

hax = tsplot(X);
title(hax(1), 'True X');
hax = tsplot(W);
title(hax(1), 'True W');

hax = tsplot(Q_gpfa.X);
title(hax(1), 'GPFA X');
hax = tsplot(Q_gpfa.W);
title(hax(1), 'GPFA W');

hax = tsplot(Q_pca.X);
title(hax(1), 'PCA X');
hax = tsplot(Q_pca.W);
title(hax(1), 'PCA W');

Yh_gpfa = Q_gpfa.W'*Q_gpfa.X;
Yh_pca = Q_pca.W'*Q_pca.X;
hax = tsplot([Y_noiseless(:,15)';
              Yh_gpfa(:,15)';
              Yh_pca(:,15)']);
title(hax(1), 'Reconstructions of observations');

% Comparison
%noise_std_gpfa = Q_gpfa.Tau(1)^(-0.5);
%noise_std_pca = Q_pca.Tau(1)^(-0.5);
recon_error_gpfa = rmse(Y_noiseless, Yh_gpfa)
recon_error_pca = rmse(Y_noiseless, Yh_pca)
test_error_gpfa = rmse(Y_noiseless(Imv), Yh_gpfa(Imv))
test_error_pca = rmse(Y_noiseless(Imv), Yh_pca(Imv))

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function applanet_prt(jdtdb)

% print aphelion/perihelion conditions

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global au pnames drsaved iap_flg iplanet1

% convert julian day to calendar date and time strings

[cdstr, utstr] = jd2str(jdtdb);

if (iap_flg == 1)

    % perihelion

    fprintf('\ntime and conditions at perihelion of ');
    disp(deblank(pnames(iplanet1, 1:7)));
    fprintf('============================================\n\n');

else

    % aphelion

    fprintf('\ntime and conditions at aphelion of ');
    disp(deblank(pnames(iplanet1, 1:7)));
    fprintf('==========================================\n\n');

end

fprintf('calendar date                ');

disp(cdstr);

fprintf('\nTDB time                     ');

disp(utstr);

fprintf('\nTDB Julian date              %16.8f \n', jdtdb);

jdutc = tdb2utc (jdtdb);

[cdstr, utstr] = jd2str(jdutc);

fprintf('\nUTC time                     ');

disp(utstr);

fprintf('\nUTC Julian date              %16.8f \n\n', jdutc);

% display heliocentric distance (kilometers and AU)

fprintf('heliocentric distance        %16.6f kilometers\n', drsaved);

fprintf('                             %16.12f AU', drsaved / au);

% toggle minima flag

iap_flg = -iap_flg;

author: Sable

repo: mcbench-benchmarks


text:
function chrpak_test011 ( )

%*****************************************************************************80
%
%% TEST011 tests CH_CAP.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 June 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST011\n' );
  fprintf ( 1, '  CH_CAP uppercases a character.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  C  CH_CAP(C)\n' );
  fprintf ( 1, '\n' );

  c = 'F';
  fprintf ( 1, '  %c      %c\n', c, ch_cap ( c ) );
  c = 'f';
  fprintf ( 1, '  %c      %c\n', c, ch_cap ( c ) );
  c = '1';
  fprintf ( 1, '  %c      %c\n', c, ch_cap ( c ) );
  c = 'b';
  fprintf ( 1, '  %c      %c\n', c, ch_cap ( c ) );
  c = '&';
  fprintf ( 1, '  %c      %c\n', c, ch_cap ( c ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [rpm] = kHz2rpm(kHz)
% Convert frequency from kilohertz to revolutions per minute.
% Chad A. Greene 2012
rpm = kHz*60*1e+3;


text:
function matRad_progress(currentIndex, totalNumberOfEvaluations)
% matRad progress bar
%
% call
%   matRad_progress(currentIndex, totalNumberOfEvaluations)
%
% input
%   currentIndex:               current iteration index
%   totalNumberOfEvaluations:   maximum iteration index
%
% output
%   graphical display of progess. make sure there is no other output
%   written during the loop to prevent confusion
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 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.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% If it's not the first step, erase the stuff printed before
if (currentIndex == 1)
    fprintf('Progress: ');
end

if (currentIndex > 1)
  Length = numel(sprintf('%3.2f %%',(currentIndex-1)/totalNumberOfEvaluations*100));
  fprintf(repmat('\b',1,Length));
end

% Print the progress tool
fprintf('%3.2f %%',currentIndex/totalNumberOfEvaluations*100);

% After the last iteration print a newline command
if (currentIndex == totalNumberOfEvaluations)
    fprintf('\n');
end

end

repo: matRad


text:
function [] = main(file, k, iterations)
    %=======================================
    %           Initialising Data
    %=======================================
    % k = 3;
    % iterations = 20;
    % file = 'yeast_test.txt';
    data = load(file);
    rows = size(data, 1);
    cols = size(data, 2);
    data = data(1:rows, 1:cols-1);
    clusters = randi([1 k], rows, 1);
    clustered_data = [data clusters];
    mean_matrix = zeros(k, cols-1);

    %=======================================
    %           Calculating Mean
    %=======================================
    for i = 1:k
       index = clustered_data(:, end) == i;
       indexed_data = clustered_data(index, 1:end-1);
       mean_matrix(i, :) = mean(indexed_data);

    end

    %=======================================
    %           Computing Error
    %=======================================
    error = get_error(clustered_data, mean_matrix);
    fprintf('After initialization: error = %.4f \n', error);

    %=======================================
    %           Starting Iterations
    %=======================================
    for p = 1:iterations
       for q = 1:rows
          %=======================================
          % Deciding which Cluster data belongs
          %=======================================
          dist = get_euclidean(data(q, :), mean_matrix);
          [minimum_row, minimum_col] = min(dist);
          clusters(q) = minimum_col;
       end
       clustered_data = [data, clusters];

       %=======================================
       %          Calculating Mean
       %=======================================
       for i = 1:k
          index = clustered_data(:, end) == i;
          indexed_data = clustered_data(index, 1:end-1);
          mean_matrix(i, :) = mean(indexed_data);
       end
       %=======================================
       %           Computing Error
       %=======================================
       error = get_error(clustered_data, mean_matrix);
       fprintf('After iteration %d: error = %.4f \n', p, error);
    end
end

function [error] = get_error(data, mean_matrix)
    %=======================================
    %           Computing Error
    %=======================================
    error = 0;
    for j = 1: size(data, 1)
       c = data(j, end);
       dist = get_euclidean(data(j, 1:end-1), mean_matrix(c, 1:end));
       error = error + dist;
    end
end

function [distance_matrix] = get_euclidean(data, mean_matrix)
    %=======================================
    %           Calculating Euclidean
    %=======================================
    dist = data - mean_matrix;
    dist = dist.^2;
    dist = sum(dist, 2);
    distance_matrix = sqrt(dist);
end

repo: Machine-Learning-Models


text:
function [ y, m, d, f ] = frac_carry_republican ( y, m, d, f )

%*****************************************************************************80
%
%% FRAC_CARRY_REPUBLICAN carries fractions from days in a Republican YMDF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, D, real F,
%    a YMDF date.
%
  if ( f < 1.0 )
    return
  end

  days = floor ( f );

  f = f - days;
  d = d + days;

  [ y, m, d ] = day_carry_republican ( y, m, d );

  return
end


text:
%% Copyright (C) 2014-2016, 2019, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym potential (@var{v})
%% @defmethodx @@sym potential (@var{v}, @var{x})
%% @defmethodx @@sym potential (@var{v}, @var{x}, @var{y})
%% Symbolic potential of a vector field.
%%
%% Finds the potential of the vector field @var{v} with respect to
%% the variables @var{x}$.  The potential is defined up to an
%% additive constant, unless the third argument is given; in which
%% case the potential is such that @var{p} is zero at the point
%% @var{y}.
%%
%% Example:
%% @example
%% @group
%% syms x y z
%% f = x*y*z;
%% g = gradient (f)
%%  @result{} g = (sym 3×1 matrix)
%%      ⎡y⋅z⎤
%%      ⎢   ⎥
%%      ⎢x⋅z⎥
%%      ⎢   ⎥
%%      ⎣x⋅y⎦
%% potential (g)
%%  @result{} (sym) x⋅y⋅z
%% @end group
%% @end example
%%
%% Return symbolic @code{nan} if the field has no potential (based
%% on checking if the Jacobian matrix of the field is
%% nonsymmetric).  For example:
%% @example
%% @group
%% syms x y
%% a = [x; x*y^2];
%% potential (a)
%%  @result{} (sym) nan
%% @end group
%% @end example
%%
%% @seealso{@@sym/gradient, @@sym/jacobian}
%% @end defmethod

function p = potential(v, x, y)

  if (nargin > 3)
    print_usage ();
  end

  assert (isvector(v), 'potential: defined for vector fields')

  if (nargin == 1)
    x = symvar(v);
  end

  % orient same as vec field
  x = reshape(x, size(v));

  if (nargin < 3)
    y = 0*x;
  end

  assert ((length(v) == length(x)) && (length(x) == length(y)), ...
          'potential: num vars must match vec length')

  cmd = { '(v, x, y) = _ins'
          'if not v.is_Matrix:'
          '    v = Matrix([v])'
          '    x = Matrix([x])'
          '    y = Matrix([y])'
          'G = v.jacobian(x)'
          'if not G.is_symmetric():'
          '    return S.NaN,'
          '_lambda = sympy.Dummy("lambda", real=True)'
          'q = y + _lambda*(x - y)'
          'vlx = v.subs(list(zip(list(x), list(q))), simultaneous=True)'
          'p = integrate((x-y).dot(vlx), (_lambda, 0, 1))'
          'return p.simplify(),' };

  p = pycall_sympy__ (cmd, sym(v), x, sym(y));

end

%!error potential (sym(1), 2, 3, 4)

%!shared x,y,z
%! syms x y z

%!test
%! % 1D
%! f = 3*x^2;
%! F = x^3;
%! assert (isequal (potential(f), F))
%! assert (isequal (potential(f, x), F))
%! assert (isequal (potential(f, x, 0), F))
%! assert (isequal (potential(f, x, 2), F - 8))

%!test
%! F = x*exp(y) + (z-1)^2;
%! f = gradient(F);
%! G = potential(f, [x;y;z], [0;1;1]);
%! assert (isAlways (G == F))

%!test
%! F = x*exp(y);
%! f = gradient(F);
%! G = potential(f);
%! assert (isAlways (G == F))

%!test
%! % no potential exists
%! syms x y
%! a = [x; x*y^2];
%! assert (isnan (potential (a)))

%!shared

%!xtest
%! % fails b/c of sympy #8458 (piecewise expr that should simplify)
%! syms x
%! f = cos(x);
%! assert (isequal (potential(f, x), sin(x)))

repo: octsympy


text:
function [p] = get_line_xyz(h)
% File:      get_line_xyz.m
% Author:    Ioannis Filippidis, jfilippidis@gmail.com
% Date:      2012.06.15
% Language:  MATLAB R2012a
% Purpose:   get line point coordinates and fix z if empty
% Copyright: Ioannis Filippidis, 2012-

% get defined data-points
x = get(h, 'XData');
y = get(h, 'YData');
z = get(h, 'ZData');

% 2d line ?
if isempty(z)
    z = zeros(size(x) );
end

p = [x; y; z];


text:
% ACOS   Inverse cosine, result in radians.
%    ACOS(X) is the arccosine of the elements of X. Complex
%    results are obtained if ABS(x) > 1.0 for some element.
%
%    See also COS, ACOSD.
%
%    Reference page in Doc Center
%       doc acos
%
%    Other functions named acos
%
%       codistributed/acos    gpuArray/acos    sym/acos    ts/acos
%

author: jmaih


text:
function [x,fval,exitflag,info] = opti_nomad(fun,x0,lb,ub,nlcon,nlrhs,xtype,opts)
%OPTI_NOMAD Solve a Global MINLP / NLP using NOMAD (Nonlinear Optimization using the MADS Algorithm)
%
%   min f(x)       subject to:    nlcon(x) <= nlrhs
%    x                            lb <= x <= ub
%                                 xi in Z
%                                 xj in {0,1}
%
%   x = opti_nomad(fun,x0) solves a Global NLP where fun is the
%   objective function and x0 is a starting guess.
%
%   x = opti_nomad(fun,x0,lb,ub) solves subject to decision variable bounds
%   lb <= x <= ub. Infinite bounds will be ignored by the solver.
%
%   x = opti_nomad(fun,x0,lb,ub,nlcon,nlrhs) solves subject to the
%   nonlinear inequality constraints nlcon(x) <= nlrhs.
%
%   x = opti_nomad(fun,...,nlrhs,xtype) solves subject to integer and
%   binary constraints specified in xtype. xtype is a string of characters
%   ('C', 'I', 'B') for continuous, integer, or binary constraints.
%
%   x = opti_nomad(fun,...,xtype,opts) uses opts to pass optiset options to
%   the solver.
%
%   [x,fval,exitflag,info] = opti_nomad(...) returns the objective value
%   at the solution, together with the solver exitflag, and an information
%   structure.
%
%   THIS IS A WRAPPER FOR NOMAD
%   See supplied Lesser GNU Public License

%   Copyright (C) 2012 Jonathan Currie (I2C2)

if(nargin < 8), opts = optiset; else opts = optiset(opts); end
if(nargin < 7), xtype = []; end
if(nargin < 6), nlrhs = []; end
if(nargin < 5), nlcon = []; end
if(nargin < 4), ub = []; end
if(nargin < 3), lb = []; end
if(nargin < 2), error('NOMAD requires at least 2 arguments'); end

%Addin nomad settings if specified
if(isfield(opts,'solverOpts') && ~isempty(opts.solverOpts))
    nopts = nomadset(opts.solverOpts);
else
    nopts = [];
end
%Override with optiset settings (note NOMAD does not accept unknown (optiset) settings)
nopts = nomadset(nopts,'max_iterations',opts.maxiter,'max_bb_eval',opts.maxfeval,...
                'max_time',opts.maxtime,'display_degree',dispLevel(opts.display),...
                'min_mesh_size',num2str(opts.tolrfun),'iterfun',opts.iterfun);

%If we have any inf bounds, set initial mesh size (C.Tribes idea)
if(any(isinf(lb)) || any(isinf(ub)) && isempty(nopts.initial_mesh_size))
    nopts = nomadset('initial_mesh_size','10'); %better heuristic on size here?
end
%Check we have a valid x0
if(isempty(x0) || any(isnan(x0)))
    error('NOMAD requires an initial guess, x0!');
end

t = tic;
% Run NOMAD
[x, fval, exitflag, iter, nfval] = nomad(fun,x0,lb,ub,nlcon,nlrhs,xtype,nopts);

%Collect Results
info.Iterations = iter;
info.FuncEvals = nfval;
info.Time = toc(t);
info.Algorithm = 'NOMAD: Nonlinear Optimization using the MADS Algorithm';

switch(exitflag)
    case 1
        info.Status = 'Converged / Target Reached';
    case 0
        info.Status = 'Exceeded Iterations / Function Evaluations / Time';
    case -1
        info.Status = 'Infeasible / Mesh Limit Reached';
    case -2
        info.Status = 'Initialization Error';
    case -3
        info.Status = 'NOMAD Error';
    case -5
        info.Status = 'User Exited';
    otherwise
        info.Status = 'NOMAD Error';
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [ a, ipvt, info ] = dgefa ( a, lda, n )

%*****************************************************************************80
%
%% DGEFA factors a real matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 June 2005
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real A(LDA,N), the matrix to be factored.
%
%    Input, integer LDA, the leading dimension of A.
%
%    Input, integer N, the order of the matrix A.
%
%    Output, real A(LDA,N), an upper triangular matrix and the multipliers
%    used to obtain it.  The factorization can be written A=L*U, where L is
%    a product of permutation and unit lower triangular matrices, and U is
%    upper triangular.
%
%    Output, integer IPVT(N), the pivot indices.
%
%    Output, integer INFO, singularity indicator.
%    0, normal value.
%    K, if U(K,K) == 0.  This is not an error condition for this subroutine,
%    but it does indicate that DGESL or DGEDI will divide by zero if called.
%    Use RCOND in DGECO for a reliable indication of singularity.
%

%
%  Gaussian elimination with partial pivoting.
%
  info = 0;

  for k = 1 : n - 1
%
%  Find L = pivot index.
%
    l = idamax ( n-k+1, a(k:n,k), 1 ) + k - 1;
    ipvt(k) = l;
%
%  Zero pivot implies this column already triangularized.
%
    if ( a(l,k) == 0.0 )
      info = k;
      continue
    end
%
%  Interchange if necessary.
%
    if ( l ~= k )
      [ a(l,k), a(k,k) ] = r8_swap ( a(l,k), a(k,k) );
    end
%
%  Compute multipliers.
%
    a(k+1:n,k) = - a(k+1:n,k) / a(k,k);
%
%  Row elimination with column indexing.
%
    for j = k+1 : n
      t = a(l,j);
      if ( l ~= k )
        a(l,j) = a(k,j);
        a(k,j) = t;
      end
      a(k+1:n,j) = daxpy ( n-k, t, a(k+1:n,k), 1, a(k+1:n,j), 1 );
    end

  end

  ipvt(n) = n;

  if ( a(n,n) == 0.0 )
    info = n;
  end

  return
end


text:
function frame2gif(volume,videoname,delaytime,resizeRate,height)
if exist('delaytime','var')~=1
    delaytime=0.5;
end
if exist('resizeRate','var')~=1
    resizeRate=1;
end
if exist('height','var')~=1
    height=0;
end

nframes=size(volume,4);

for i=1:nframes
    im=volume(:,:,:,i);
    if height~=0
        [h,w,nchannels]=size(im);
        resizeRate=height/h;
    end
    if resizeRate<1
        im=imresize(imfilter(im,fspecial('gaussian',5,0.7),'same','replicate'),resizeRate,'bicubic');
    else if resizeRate>1
            im=imresize(im,resizeRate,'bicubic');
        end
    end
    [X,map]=rgb2ind(im,256);
    if i==1
        imwrite(X,map,videoname,'DelayTime',delaytime,'LoopCount',Inf);
    else
        imwrite(X,map,videoname,'WriteMode','append','DelayTime',delaytime);
    end
end

author: qingsenyangit

repo: AHDRNet


text:
function [a]=minus(b,c)
%A=B-C
%   [A]=MINUS(B,C) Subtraction of two tensors
%
%Adds two TT-tensors B and C
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

a=b+(-1.0)*c;
return
end

author: oseledets

repo: TT-Toolbox


text:
function circle = makeCartCircle(radius, num_points, center_pos, arc_angle, plot_circle)
%MAKECARTCIRCLE     Create a 2D Cartesian circle or arc.
%
% DESCRIPTION:
%       MakeCartCircle creates a 2 x num_points array of the Cartesian
%       coordinates of points evenly distributed over a circle or arc (if
%       arc_angle is given).
%
% USAGE:
%       circle = makeCartCircle(radius, num_points)
%       circle = makeCartCircle(radius, num_points, center_pos)
%       circle = makeCartCircle(radius, num_points, center_pos, arc_angle)
%       circle = makeCartCircle(radius, num_points, center_pos, arc_angle, plot_circle)
%
% INPUTS:
%       radius          - circle radius [m]
%       num_points      - number of points in the circle
%
% OPTIONAL INPUTS:
%       center_pos      - [x, y] position of the circle center [m]
%                         (default = [0, 0])
%       arc_angle       - arc angle for incomplete circle [radians]
%                         (default = 2*pi)
%       plot_circle     - Boolean controlling whether the Cartesian points
%                         are plotted (default = false)
%
% OUTPUTS:
%       circle          - 2 x num_points array of Cartesian coordinates
%
% ABOUT:
%       author          - Bradley Treeby
%       date            - 5th June 2009
%       last update     - 20th March 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also cart2grid, makeCartSphere, makeCircle

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

% check for plot_circle input
if nargin < 5 || isempty(plot_circle)
    plot_circle = false;
end

% check for arc_angle input
if nargin < 4 || isempty(arc_angle)
    arc_angle = 2*pi;
    full_circle = true;
elseif arc_angle == 2*pi;
    full_circle = true;
else
    full_circle = false;
end

% check for center_pos input
if nargin < 3 || isempty(center_pos)
    cx = 0;
    cy = 0;
else
    cx = center_pos(1);
    cy = center_pos(2);
end

% ensure there is only a total of num_points including the endpoints when
% arc_angle is not equal to 2*pi
if ~full_circle
    num_points = num_points - 1;
end

% create angles
angles = (0:(num_points))*arc_angle/(num_points) + pi/2;

% discard repeated final point if arc_angle is equal to 2*pi
if full_circle
    angles = angles(1:end-1);
end

% create cartesian grid
% circle = flipud([radius*cos(angles); radius*sin(-angles)]);   % B.0.3
circle = ([radius*cos(angles); radius*sin(-angles)]);           % B.0.4

% offset if needed
circle(1, :) = circle(1, :) + cx;
circle(2, :) = circle(2, :) + cy;

% plot results
if plot_circle

    % select suitable axis scaling factor
    [x_sc, scale, prefix] = scaleSI(max(abs(circle(:))));

    % create the figure
    figure;
    plot(circle(2,:)*scale, circle(1,:)*scale, 'b.');
    set(gca, 'YDir', 'reverse');
    xlabel(['y-position [' prefix 'm]']);
    ylabel(['x-position [' prefix 'm]']);
    axis equal;

end

repo: Aero-matlab


text:
method = 'manual';  % random, uniform, manual
frames = 200;       % 50 100 150 200

ms_regular_training('baseline/highway', method, frames);
ms_regular_training('baseline/office', method, frames);
ms_regular_training('baseline/pedestrians', method, frames);
ms_regular_training('baseline/PETS2006', method, frames);

ms_regular_training('badWeather/blizzard', method, frames);
ms_regular_training('badWeather/skating', method, frames);
ms_regular_training('badWeather/snowFall', method, frames);
ms_regular_training('badWeather/wetSnow', method, frames);

ms_regular_training('cameraJitter/badminton', method, frames);
ms_regular_training('cameraJitter/boulevard', method, frames);
ms_regular_training('cameraJitter/sidewalk', method, frames);
ms_regular_training('cameraJitter/traffic', method, frames);

ms_regular_training('dynamicBackground/boats', method, frames);
ms_regular_training('dynamicBackground/canoe', method, frames);
ms_regular_training('dynamicBackground/fall', method, frames);
ms_regular_training('dynamicBackground/fountain01', method, frames);
ms_regular_training('dynamicBackground/fountain02', method, frames);
ms_regular_training('dynamicBackground/overpass', method, frames);

ms_regular_training('intermittentObjectMotion/abandonedBox', method, frames);
ms_regular_training('intermittentObjectMotion/parking', method, frames);
ms_regular_training('intermittentObjectMotion/sofa', method, frames);
ms_regular_training('intermittentObjectMotion/streetLight', method, frames);
ms_regular_training('intermittentObjectMotion/tramstop', method, frames);
ms_regular_training('intermittentObjectMotion/winterDriveway', method, frames);

ms_regular_training('lowFramerate/port_0_17fps', method, frames);
ms_regular_training('lowFramerate/tramCrossroad_1fps', method, frames);
ms_regular_training('lowFramerate/tunnelExit_0_35fps', method, frames);
ms_regular_training('lowFramerate/turnpike_0_5fps', method, frames);

ms_regular_training('nightVideos/bridgeEntry', method, frames);
ms_regular_training('nightVideos/busyBoulvard', method, frames);
ms_regular_training('nightVideos/fluidHighway', method, frames);
ms_regular_training('nightVideos/streetCornerAtNight', method, frames);
ms_regular_training('nightVideos/tramStation', method, frames);
ms_regular_training('nightVideos/winterStreet', method, frames);

ms_regular_training('PTZ/continuousPan', method, frames);
ms_regular_training('PTZ/intermittentPan', method, frames);
ms_regular_training('PTZ/twoPositionPTZCam', method, frames);
ms_regular_training('PTZ/zoomInZoomOut', method, frames);

ms_regular_training('shadow/backdoor', method, frames);
ms_regular_training('shadow/bungalows', method, frames);
ms_regular_training('shadow/busStation', method, frames);
ms_regular_training('shadow/copyMachine', method, frames);
ms_regular_training('shadow/cubicle', method, frames);
ms_regular_training('shadow/peopleInShade', method, frames);

ms_regular_training('thermal/corridor', method, frames);
ms_regular_training('thermal/diningRoom', method, frames);
ms_regular_training('thermal/lakeSide', method, frames);
ms_regular_training('thermal/library', method, frames);
ms_regular_training('thermal/park', method, frames);

ms_regular_training('turbulence/turbulence0', method, frames);
ms_regular_training('turbulence/turbulence1', method, frames);
ms_regular_training('turbulence/turbulence2', method, frames);
ms_regular_training('turbulence/turbulence3', method, frames);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function index_test01 ( )

%*****************************************************************************80
%
%% INDEX_TEST01 tests INDEX0 and INDEX1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 November 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'INDEX_TEST01\n' );
  fprintf ( 1, '  INDEX0 indexes a 1D array with zero base,\n' );
  fprintf ( 1, '  INDEX1 indexes a 1D array with  unit base.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             Min Index   Max\n' );
  fprintf ( 1, '\n' );

  i_min = 1;
  i = 3;
  i_max = 5;
  fprintf ( 1, '  1D Index  %4d  %4d  %4d\n', i_min,     i,     i_max );

  value = index0 ( i_min, i, i_max );
  index_min = 0;
  index_max = index_min + i_max - i_min;
  fprintf ( 1, '  Index0:   %4d  %4d  %4d\n', index_min, value, index_max );

  value = index1 ( i_min, i, i_max );
  index_min = 1;
  index_max = index_min + i_max - i_min;
  fprintf ( 1, '  Index1:   %4d  %4d  %4d\n', index_min, value, index_max );

  return
end


text:
classdef LineSegmentDetector < handle
    %LINESEGMENTDETECTOR  Line segment detector class
    %
    % Following the algorithm described at [Rafael12].
    %
    % ## References
    % [Rafael12]:
    % > Rafael Grompone von Gioi, Jeremie Jakubowicz, Jean-Michel Morel,
    % > and Gregory Randall. Lsd: a line segment detector. 2012.
    %
    % See also: cv.LineSegmentDetector.LineSegmentDetector, cv.HoughLines,
    %  cv.FastLineDetector, houghlines
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    methods
        function this = LineSegmentDetector(varargin)
            %LINESEGMENTDETECTOR  Creates a LineSegmentDetector object and initializes it
            %
            %     lsd = cv.LineSegmentDetector()
            %     lsd = cv.LineSegmentDetector('OptionName', optionValue, ...)
            %
            % ## Options
            % * __Refine__ The way found lines will be refined, one of:
            %   * __None__ No refinement applied.
            %   * __Standard__ (default) Standard refinement is applied. E.g.
            %     breaking arches into smaller straighter line approximations.
            %   * __Advanced__ Advanced refinement. Number of false alarms is
            %     calculated, lines are refined through increase of precision,
            %     decrement in size, etc.
            % * __Scale__ The scale of the image that will be used to find the
            %   lines. In the range `[0,1)`. default 0.8
            % * __SigmaScale__ Sigma for Gaussian filter. It is computed as
            %   `sigma = SigmaScale/Scale`. default 0.6
            % * __QuantError__ Bound to the quantization error on the gradient
            %   norm. default 2.0
            % * __AngleTol__ Gradient angle tolerance in degrees. default 22.5
            % * __DetectionThreshold__ Detection threshold:
            %   `-log10(NFA) > DetectionThreshold`. Used only when advanced
            %   refinement is chosen. default 0
            % * __MinDensity__ Minimal density of aligned region points in the
            %   enclosing rectangle. default 0.7
            % * __NBins__ Number of bins in pseudo-ordering of gradient
            %   modulus. default 1024
            %
            % The cv.LineSegmentDetector algorithm is defined using the
            % standard values. Only advanced users may want to edit those, as
            % to tailor it for their own application.
            %
            % See also: cv.LineSegmentDetector.detect
            %
            this.id = LineSegmentDetector_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.LineSegmentDetector
            %
            if isempty(this.id), return; end
            LineSegmentDetector_(this.id, 'delete');
        end
    end

    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.LineSegmentDetector.empty
            %
            LineSegmentDetector_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Checks if algorithm object is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the algorithm object is empty
            %   (e.g. in the very beginning or after unsuccessful read).
            %
            % See also: cv.LineSegmentDetector.clear
            %
            b = LineSegmentDetector_(this.id, 'empty');
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.LineSegmentDetector.save, cv.LineSegmentDetector.load
            %
            name = LineSegmentDetector_(this.id, 'getDefaultName');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in a file storage.
            %
            % See also: cv.LineSegmentDetector.load
            %
            LineSegmentDetector_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from a file storage.
            % The previous model state is discarded.
            %
            % See also: cv.LineSegmentDetector.save
            %
            LineSegmentDetector_(this.id, 'load', fname_or_str, varargin{:});
        end
    end

    methods
        function [lines, width, prec, nfa] = detect(this, img)
            %DETECT  Finds lines in the input image
            %
            %     lines = lsd.detect(img)
            %     [lines, width, prec, nfa] = lsd.detect(img)
            %
            % ## Input
            % * __img__ A grayscale (`uint8`) input image.
            %
            % ## Output
            % * __lines__ A cell-array of 4-element vectors of the form
            %   `{[x1, y1, x2, y2], ..}` specifying the beginning and ending
            %   point of a line. Where point 1 `[x1,y1]` is the start, point 2
            %   `[x2,y2] the end. Returned lines are strictly oriented
            %   depending on the gradient.
            % * __width__ Vector of widths of the regions, where the lines are
            %   found. E.g. Width of line.
            % * __prec__ Vector of precisions with which the lines are found.
            % * __nfa__ Vector containing number of false alarms in the line
            %   region, with precision of 10%. The bigger the value,
            %   logarithmically better the detection. This vector will be
            %   calculated only when the object refine type is 'Advanced',
            %   empty otherwise.
            %   * -1 corresponds to 10 mean false alarms
            %   * 0 corresponds to 1 mean false alarm
            %   * 1 corresponds to 0.1 mean false alarms
            %
            % See also: cv.LineSegmentDetector.drawSegments
            %
            [lines, width, prec, nfa] = LineSegmentDetector_(this.id, 'detect', img);
        end

        function img = drawSegments(this, img, lines)
            %DRAWSEGMENTS  Draws the line segments on a given image
            %
            %     img = lsd.drawSegments(img, lines)
            %
            % ## Input
            % * __img__ The image, where the lines will be drawn. Should be
            %   bigger or equal to the image, where the lines were found.
            % * __lines__ A vector of the lines that needed to be drawn.
            %
            % ## Output
            % * __img__ Output image with drawn line segments.
            %
            % See also: cv.LineSegmentDetector.compareSegments
            %
            img = LineSegmentDetector_(this.id, 'drawSegments', img, lines);
        end

        function [img, count] = compareSegments(this, sz, lines1, lines2, varargin)
            %COMPARESEGMENTS  Draws two groups of lines in blue and red, counting the non overlapping (mismatching) pixels
            %
            %     [img, count] = lsd.compareSegments(sz, lines1, lines2)
            %     [...] = lsd.compareSegments(..., 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __sz__ The size of the image, where `lines1` and `lines2` were
            %   found `[w,h]`.
            % * __lines1__ The first group of lines that needs to be drawn. It
            %   is visualized in blue color.
            % * __lines2__ The second group of lines. They visualized in red
            %   color.
            %
            % ## Output
            % * __img__ color image with the two groups of lines drawn.
            % * __count__ count of non overlapping (mismatching) pixels.
            %
            % ## Options
            % * __Image__ Optional image, where the lines will be drawn. The
            %   image should be color (3-channel) in order for `lines1` and
            %   `lines2` to be drawn in the above mentioned colors.
            %
            % See also: cv.LineSegmentDetector.drawSegments
            %
            [img, count] = LineSegmentDetector_(this.id, 'compareSegments', ...
                sz, lines1, lines2, varargin{:});
        end
    end
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function machine_test ( )

%*****************************************************************************80
%
%% MAIN_TEST tests the MACHINE library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 April 2007
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MACHINE_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the MACHINE library.\n' );

  d1mach_test ( );
  i1mach_test ( );
  r1mach_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MACHINE_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function [image_roidb, bbox_means, bbox_stds] = mask_prepare_image_roidb_voc(conf, imdbs, roidbs, bbox_means, bbox_stds)
% [image_roidb, bbox_means, bbox_stds] = fast_rcnn_prepare_image_roidb(conf, imdbs, roidbs, cache_img, bbox_means, bbox_stds)
%   Gather useful information from imdb and roidb
%   pre-calculate mean (bbox_means) and std (bbox_stds) of the regression
%   term for normalization
% --------------------------------------------------------
% Fast R-CNN
% Reimplementation based on Python Fast R-CNN (https://github.com/rbgirshick/fast-rcnn)
% Copyright (c) 2015, Shaoqing Ren
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

    if ~exist('bbox_means', 'var')
        bbox_means = [];
        bbox_stds = [];
    end

    if ~iscell(imdbs)
        imdbs = {imdbs};
        roidbs = {roidbs};
    end

    imdbs = imdbs(:);
    roidbs = roidbs(:);

    image_roidb = ...
        cellfun(@(x, y) ... // @(imdbs, roidbs)
               arrayfun(@(z) ... //@([1:length(x.image_ids)])
                        struct('image_path', x.image_at(z), 'image_seg_path', x.image_seg_at(z), 'image_id', x.image_ids{z}, ...
                        'im_size', x.sizes(z, :), 'imdb_name', x.name, 'overlap', y.rois(z).overlap, 'boxes', y.rois(z).boxes, ...
                        'class', y.rois(z).class, 'image', [], 'bbox_targets', [], 'mask_targets', []), ...
                        [1:length(x.image_ids)]', 'UniformOutput', true),...
        imdbs, roidbs, 'UniformOutput', false);

    image_roidb = cat(1, image_roidb{:});

    % enhance roidb to contain bounding-box regression targets
    [image_roidb, bbox_means, bbox_stds] = append_bbox_regression_targets(conf, image_roidb, bbox_means, bbox_stds);
end

function [image_roidb, means, stds] = append_bbox_regression_targets(conf, image_roidb, means, stds)
    % means and stds -- (k+1) * 4, include background class

    num_images = length(image_roidb);
    % Infer number of classes from the number of columns in gt_overlaps
    num_classes = size(image_roidb(1).overlap, 2);
    valid_imgs = true(num_images, 1);
    for i = 1:num_images
       rois = image_roidb(i).boxes;
       [image_roidb(i).bbox_targets, image_roidb(i).mask_targets, valid_imgs(i)] = ...
           compute_targets(conf, rois, image_roidb(i).overlap, image_roidb(i).image_seg_path);
    end
    if ~all(valid_imgs)
        image_roidb = image_roidb(valid_imgs);
        num_images = length(image_roidb);
        fprintf('Warning: fast_rcnn_prepare_image_roidb: filter out %d images, which contains zero valid samples\n', sum(~valid_imgs));
    end
    if ~(exist('means', 'var') && ~isempty(means) && exist('stds', 'var') && ~isempty(stds))
        % Compute values needed for means and stds
        % var(x) = E(x^2) - E(x)^2
        class_counts = zeros(num_classes, 1) + eps;
        sums = zeros(num_classes, 4);
        squared_sums = zeros(num_classes, 4);
        for i = 1:num_images
           targets = image_roidb(i).bbox_targets;
           for cls = 1:num_classes
              cls_inds = find(targets(:, 1) == cls);
              if ~isempty(cls_inds)
                 class_counts(cls) = class_counts(cls) + length(cls_inds);
                 sums(cls, :) = sums(cls, :) + sum(targets(cls_inds, 2:end), 1);
                 squared_sums(cls, :) = squared_sums(cls, :) + sum(targets(cls_inds, 2:end).^2, 1);
              end
           end
        end

        means = bsxfun(@rdivide, sums, class_counts);
        stds = (bsxfun(@minus, bsxfun(@rdivide, squared_sums, class_counts), means.^2)).^0.5;
        % add background class
        means = [0, 0, 0, 0; means];
        stds = [0, 0, 0, 0; stds];

    end
    % Normalize targets
    for i = 1:num_images
        targets = image_roidb(i).bbox_targets;
        for cls = 1:num_classes
            cls_inds = find(targets(:, 1) == cls);
            if ~isempty(cls_inds)
                image_roidb(i).bbox_targets(cls_inds, 2:end) = ...
                    bsxfun(@minus, image_roidb(i).bbox_targets(cls_inds, 2:end), means(cls+1, :));
                image_roidb(i).bbox_targets(cls_inds, 2:end) = ...
                    bsxfun(@rdivide, image_roidb(i).bbox_targets(cls_inds, 2:end), stds(cls+1, :));
            end
        end
    end
end

function [bbox_targets, mask_targets, is_valid] = compute_targets(conf, rois, overlap, img_seg_path)
%function [bbox_targets, is_valid] = compute_targets(conf, rois, overlap)
    overlap = full(overlap);

    [max_overlaps, max_labels] = max(overlap, [], 2);

    % ensure ROIs are floats
    rois = single(rois);

    bbox_targets = zeros(size(rois, 1), 5, 'single');
    mask_targets = zeros(size(rois,1), 28, 28, 'single');
    % Indices of ground-truth ROIs
    gt_inds = find(max_overlaps == 1);

    if ~isempty(gt_inds)
        % Indices of examples for which we try to make predictions
        ex_inds = find(max_overlaps >= conf.bbox_thresh);

        % Get IoU overlap between each ex ROI and gt ROI
        ex_gt_overlaps = boxoverlap(rois(ex_inds, :), rois(gt_inds, :));

        %Small error between gt overlap and calculated overlap
        assert(all(abs(max(ex_gt_overlaps, [], 2) - max_overlaps(ex_inds)) < 10^-4));

        % Find which gt ROI each ex ROI has max overlap with:
        % this will be the ex ROI's gt target
        [~, gt_assignment] = max(ex_gt_overlaps, [], 2);
        gt_rois = rois(gt_inds(gt_assignment), :);
        ex_rois = rois(ex_inds, :);

        [regression_label] = fast_rcnn_bbox_transform(ex_rois, gt_rois);

        bbox_targets(ex_inds, :) = [max_labels(ex_inds), regression_label];
        seg_image = imread(img_seg_path);
        [max_h, max_w] = size(seg_image);

        %Get mask targets
        for i=1:size(ex_rois,1)
          roi_orig = ex_rois(round(i),:); %convert to int before accessing indices
          mask = seg_image(max(1, roi_orig(2)): min(max_h, roi_orig(4)), max(1, roi_orig(1)):min(max_w, roi_orig(3)));
          mask = imresize(mask, [28, 28], 'nearest');
          mask(mask == 255) = 0; % We consider void boundaries as background
          mask_targets(i,:,:) = mask;
        end
    end
    % Select foreground ROIs as those with >= fg_thresh overlap
    is_fg = max_overlaps >= conf.fg_thresh;
    % Select background ROIs as those within [bg_thresh_lo, bg_thresh_hi)
    is_bg = max_overlaps < conf.bg_thresh_hi & max_overlaps >= conf.bg_thresh_lo;

    % check if there is any fg or bg sample. If no, filter out this image
    is_valid = true;
    if ~any(is_fg | is_bg)
        is_valid = false;
    end
end

author: jasjeetIM

repo: Mask-RCNN

save_path: github-repos/MATLAB/jasjeetIM-Mask-RCNN


text:
function prepro = spm_cfg_eeg_inv_prepro
% Configuration file for configuring imaging source inversion reconstruction
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_cfg_eeg_inv_prepro.m 6444 2015-05-21 11:15:48Z guillaume $

D = cfg_files;
D.tag = 'D';
D.name = 'M/EEG datasets';
D.filter = 'mat';
D.num = [1 Inf];
D.help = {'Select the M/EEG mat files.'};

val = cfg_entry;
val.tag = 'val';
val.name = 'Inversion index';
val.strtype = 'n';
val.help = {'Index of the cell in D.inv where the forward model can be found and the results will be stored.'};
val.val = {1};

all = cfg_const;
all.tag = 'all';
all.name = 'All';
all.val  = {1};

condlabel = cfg_entry;
condlabel.tag = 'condlabel';
condlabel.name = 'Condition label';
condlabel.strtype = 's';
condlabel.val = {''};

conditions = cfg_repeat;
conditions.tag = 'conditions';
conditions.name = 'Conditions';
conditions.help = {'Specify the labels of the conditions to be included in the inversion'};
conditions.num  = [1 Inf];
conditions.values  = {condlabel};
conditions.val = {condlabel};

whatconditions = cfg_choice;
whatconditions.tag = 'whatconditions';
whatconditions.name = 'What conditions to include?';
whatconditions.values = {all, conditions};
whatconditions.val = {all};

standard = cfg_const;
standard.tag = 'standard';
standard.name = 'Standard';
standard.help = {'Use default settings for the inversion'};
standard.val  = {1};

invtype = cfg_menu;
invtype.tag = 'invtype';
invtype.name = 'Inversion type';
invtype.help = {'Select the desired inversion type'};
invtype.labels = {'GS', 'ARD', 'MSP (GS+ARD)' 'COH', 'IID', 'EBB'};
invtype.values = {'GS', 'ARD', 'MSP', 'LOR', 'IID', 'EBB'};
invtype.val = {'GS'};

woi = cfg_entry;
woi.tag = 'woi';
woi.name = 'Time window of interest';
woi.strtype = 'r';
woi.num = [1 2];
woi.val = {[-Inf Inf]};
woi.help = {'Time window to include in the inversion (ms)'};

foi = cfg_entry;
foi.tag = 'foi';
foi.name = 'Frequency window of interest';
foi.strtype = 'r';
foi.num = [1 2];
foi.val = {[0 256]};
foi.help = {'Frequency window (the same as high-pass and low-pass in the GUI)'};

hanning = cfg_menu;
hanning.tag = 'hanning';
hanning.name = 'PST Hanning window';
hanning.help = {'Multiply the time series by a Hanning taper to emphasize the central part of the response.'};
hanning.labels = {'yes', 'no'};
hanning.values = {1, 0};
hanning.val = {1};

priorsmask  = cfg_files;
priorsmask.tag = 'priorsmask';
priorsmask.name = 'Priors file';
priorsmask.filter = '(.*\.gii$)|(.*\.mat$)|(.*\.nii(,\d+)?$)|(.*\.img(,\d+)?$)';
priorsmask.num = [0 1];
priorsmask.help = {'Select a mask or a mat file with priors.'};
priorsmask.val = {{''}};

space = cfg_menu;
space.tag = 'space';
space.name = 'Prior image space';
space.help = {'Space of the mask image.'};
space.labels = {'MNI', 'Native'};
space.values = {1, 0};
space.val = {1};

priors = cfg_branch;
priors.tag = 'priors';
priors.name = 'Source priors';
priors.help = {'Restrict solutions to pre-specified VOIs'};
priors.val  = {priorsmask, space};

locs  = cfg_entry;
locs.tag = 'locs';
locs.name = 'Source locations';
locs.strtype = 'r';
locs.num = [Inf 3];
locs.help = {'Input source locations as n x 3 matrix'};
locs.val = {zeros(0, 3)};

radius = cfg_entry;
radius.tag = 'radius';
radius.name = 'Radius of VOI (mm)';
radius.strtype = 'r';
radius.num = [1 1];
radius.val = {32};

restrict = cfg_branch;
restrict.tag = 'restrict';
restrict.name = 'Restrict solutions';
restrict.help = {'Restrict solutions to pre-specified VOIs'};
restrict.val  = {locs, radius};

custom = cfg_branch;
custom.tag = 'custom';
custom.name = 'Custom';
custom.help = {'Define custom settings for the inversion'};
custom.val  = {woi, foi, hanning};

isstandard = cfg_choice;
isstandard.tag = 'isstandard';
isstandard.name = 'Inversion parameters';
isstandard.help = {'Choose whether to use standard or custom inversion parameters.'};
isstandard.values = {standard, custom};
isstandard.val = {standard};

modality = cfg_menu;
modality.tag = 'modality';
modality.name = 'Select modalities';
modality.help = {'Select modalities for the inversion (only relevant for multimodal datasets).'};
modality.labels = {'All', 'EEG', 'MEG', 'MEGPLANAR', 'EEG+MEG', 'MEG+MEGPLANAR', 'EEG+MEGPLANAR'};
modality.values = {
    {'All'}
    {'EEG'}
    {'MEG'}
    {'MEGPLANAR'}
    {'EEG', 'MEG'}
    {'MEG', 'MEGPLANAR'}
    {'EEG', 'MEGPLANAR'}
    }';
modality.val = {{'All'}};

prepro = cfg_exbranch;
prepro.tag = 'prepro';
prepro.name = 'Inversion preprocessing';
prepro.val = {D, val, whatconditions, isstandard, modality};
prepro.help = {'Run imaging source reconstruction'};
prepro.prog = @run_prepro;
prepro.vout = @vout_prepro;
prepro.modality = {'EEG'};

function  out = run_prepro(job)

D = spm_eeg_load(job.D{1});

inverse = [];
if isfield(job.whatconditions, 'condlabel')
    inverse.trials = job.whatconditions.condlabel;
end

if isfield(job.isstandard, 'custom')
    %inverse.type = job.isstandard.custom.invtype;
    inverse.woi  = fix([max(min(job.isstandard.custom.woi), 1000*D.time(1)) min(max(job.isstandard.custom.woi), 1000*D.time(end))]);
    inverse.Han  = job.isstandard.custom.hanning;
    inverse.lpf  =  fix(min(job.isstandard.custom.foi));
    inverse.hpf  =  fix(max(job.isstandard.custom.foi));

end

[mod, list] = modality(D, 1, 1);
if strcmp(job.modality{1}, 'All')
    inverse.modality  = list;
else
    inverse.modality  = intersect(list, job.modality);
end

if numel(inverse.modality) == 1
    inverse.modality = inverse.modality{1};
end

D = {};

for i = 1:numel(job.D)
    D{i} = spm_eeg_load(job.D{i});

    D{i}.val = job.val;

    D{i}.con = 1;

    if ~isfield(D{i}, 'inv')
        error('Forward model is missing for subject %d.', i);
    elseif  numel(D{i}.inv)<D{i}.val || ~isfield(D{i}.inv{D{i}.val}, 'forward')
        if D{i}.val>1 && isfield(D{i}.inv{D{i}.val-1}, 'forward')
            D{i}.inv{D{i}.val} = D{i}.inv{D{i}.val-1};
            warning('Duplicating the last forward model for subject %d.', i);
        else
            error('Forward model is missing for subject %d.', i);
        end
    end

    D{i}.inv{D{i}.val}.inverse = inverse;
end

D = spm_eeg_invert_prepro(D);

if ~iscell(D)
    D = {D};
end

for i = 1:numel(D)
    save(D{i});
end

out.D = job.D;

function dep = vout_prepro(job)
% Output is always in field "D", no matter how job is structured
dep = cfg_dep;
dep.sname = 'M/EEG dataset(s) after imaging source reconstruction';
% reference field "D" from output
dep.src_output = substruct('.','D');
% this can be entered into any evaluated input
dep.tgt_spec   = cfg_findspec({{'filter','mat'}});

author: spm


text:
function marginal = marginal_nodes(engine, nodes, t, fam)
% MARGINAL_NODES Compute the marginal on the specified query nodes (bk)
%
%   marginal = marginal_nodes(engine, i, t)
% returns Pr(X(i,t) | Y(1:T)), where X(i,t) is the i'th node in the t'th slice.
%
%   marginal = marginal_nodes(engine, query, t)
% returns Pr(X(query(1),t), ... X(query(end),t) | Y(1:T)),
% where 't' specifies the time slice of the earliest node in the query.
% 'query' cannot span more than 2 time slices.
%
% Example:
% Consider a DBN with 2 nodes per slice.
% Then t=2, nodes=[1 3] refers to node 1 in slice 2 and node 1 in slice 3.

if nargin < 3, t = 1; end
if nargin < 4, fam = 0; else fam = 1; end

% clpot{t} contains slice t-1 and t
% Example
% clpot #: 1    2    3
% slices:  1  1,2  2,3
% For filtering, we must take care not to take future evidence into account.
% For smoothing, clpot{1} does not exist.

bnet = bnet_from_engine(engine);
ss = length(bnet.intra);

if t < engine.T
  slice = t+1;
  nodes2 = nodes;
else % earliest t is T, so all nodes fit in one slice
  slice = engine.T;
  nodes2 = nodes + ss;
end

c = clq_containing_nodes(engine.jtree_engine, nodes2, fam);
assert(c >= 1);

%disp(['computing marginal on ' num2str(nodes) ' t = ' num2str(t)]);
%disp(['using ' num2str(nodes2) ' slice = ' num2str(slice) 'clq = ' num2str(c)]);

bigpot = engine.clpot{c, slice};

pot = marginalize_pot(bigpot, nodes2, engine.maximize);
marginal = pot_to_marginal(pot);

% we convert the domain to the unrolled numbering system
% so that update_ess extracts the right evidence.
marginal.domain = nodes+(t-1)*ss;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%% Copyright (C) 2014, 2016, 2019 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym null (@var{A})
%% Basis for the nullspace of a symbolic matrix.
%%
%% Return a matrix whose columns are a basis for the nullspace of
%% the matrix.
%% Examples:
%% @example
%% @group
%% A = sym([1 1; 2 0]);
%% null (A)
%%   @result{} (sym) []  (empty 2×0 matrix)
%% @end group
%%
%% @group
%% A = sym([1 2; 1 2]);
%% null (A)
%%   @result{} (sym 2×1 matrix)
%%
%%       ⎡-2⎤
%%       ⎢  ⎥
%%       ⎣1 ⎦
%% @end group
%%
%% @group
%% A = sym(zeros(2,2));
%% null (A)
%%   @result{} (sym 2×2 matrix)
%%
%%       ⎡1  0⎤
%%       ⎢    ⎥
%%       ⎣0  1⎦
%% @end group
%% @end example
%%
%% @seealso{@@sym/rank, @@sym/orth}
%% @end defmethod

function r = null(A)

  cmd = { 'A = _ins[0]'
          'if not A.is_Matrix:'
          '    A = sympy.Matrix([A])'
          'ns = A.nullspace()'
          'if len(ns) == 0:'
          '    return sympy.zeros(A.cols, 0),'
          'return sympy.Matrix.hstack(*ns),' };

  r = pycall_sympy__ (cmd, A);

end

%!test
%! A = sym([1 2; 3 4]);
%! assert (isempty (null (A)))

%!assert (isempty (null (sym(4))))

%!test
%! A = sym([1 2 3; 3 4 5]);
%! assert (isequal (null(A), sym([1;-2;1])))


text:
function varargout = filetype_check_uri(filename, ftyp)

% FILETYPE_CHECK_URI
%
% Supported URIs are
%   buffer://<host>:<port>
%   fifo://<filename>
%   global://<varname>
%   mysql://<user>:<password>@<host>:<port>
%   rfb://<password>@<host>:<port>
%   serial:<port>?key1=value1&key2=value2&...
%   shm://<filename>
%   tcp://<host>:<port>
%   udp://<host>:<port>
%
% The URI schemes supproted by these function are not the official schemes.
% See the documentation included inside this function for more details.
% RFC4395 defines an IANA-maintained registry of URI Schemes. See also
% http://www.iana.org/assignments/uri-schemes.html and
% http://en.wikipedia.org/wiki/URI_scheme#Generic_syntax.

% Copyright (C) 2007, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% 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: filetype_check_uri.m 2885 2011-02-16 09:41:58Z roboos $

% these are for remembering the type on subsequent calls with the same input arguments
persistent previous_argin previous_argout

if nargin<2
  % ensure that all input arguments are defined
  ftyp = [];
end

current_argin = {filename, ftyp};
if isequal(current_argin, previous_argin)
  % don't do the detection again, but return the previous value from cache
  varargout = previous_argout;
  return
end

sep = find(filename==':');
if ~isempty(sep)
  scheme = filename(1:(sep(1)-1));
else
  scheme = [];
end

if nargin>1
  % only return true or false
  varargout{1} = strcmp(scheme, ftyp);
else
  % return the full details of this URI scheme
  switch scheme
    case 'shm'
      % shm://<filename>
      % the filename is optional, usually it can and will be read from shared memory
      if length(filename)>6
        varargout{1} = filename(7:end);
      else
        varargout{1} = [];
      end

    case 'fifo'
      % fifo://<filename>
      varargout{1} = filename(8:end);

    case 'global'
      % global://<varname>
      varargout{1} = filename(10:end);

    case 'buffer'
      % buffer://<host>:<port>
      tok = tokenize(filename(10:end), ':');
      varargout{1} = tok{1};
      varargout{2} = str2num(tok{2});

    case 'tcp'
      % tcp://<host>:<port>
      tok = tokenize(filename(7:end), ':');
      varargout{1} = tok{1};
      varargout{2} = str2num(tok{2});

    case 'udp'
      % udp://<host>:<port>
      tok = tokenize(filename(7:end), ':');
      varargout{1} = tok{1};
      varargout{2} = str2num(tok{2});

    case 'rfb'
      % rfb://<password>@<host>:<port>
      tok0 = tokenize(filename(7:end), '@');
      tok1 = tokenize(tok0{2}, ':');
      varargout{1} = tok0{1};
      varargout{2} = tok1{1};
      varargout{3} = str2num(tok1{2});

    case 'mysql'
      % mysql://<user>:<password>@<host>:<port>
      tok0 = tokenize(filename(9:end), '@');
      tok1 = tokenize(tok0{1}, ':');
      tok2 = tokenize(tok0{2}, ':');
      varargout{1} = tok1{1};
      varargout{2} = tok1{2};
      varargout{3} = tok2{1};
      if length(tok2)>1
        varargout{4} = str2num(tok2{2});
      else
        varargout{4} = [];
      end

    case 'serial'
      % serial:<Port>?key1=value1&key2=value2&...
      % the supported optional arguments are
      %   BaudRate
      %   DataBits
      %   DataTerminalReady
      %   FlowControl
      %   Parity
      %   Port
      %   ReadAsyncMode
      %   RequestToSend
      %   StopBits
      %   Terminator
      tok0 = tokenize(filename(8:end), '?');
      varargout{1} = tok0{1};
      varargout{2} = [];
      if length(tok0)>1
        tok1 = tokenize(tok0{2}, '&');
        for i=1:length(tok1)
          tok2 = tokenize(tok1{i}, '=');
          opt{(i-1)*2+1} = tok2{1};
          opt{(i-1)*2+2} = tok2{2};
        end
        varargout{2} = opt;
      end

    otherwise
      error('unsupported scheme in URI')
  end
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
current_argout  = varargout;
previous_argin  = current_argin;
previous_argout = current_argout;

return % main()

% RFC4395 defines an IANA-maintained registry of URI Schemes.
% see also http://www.iana.org/assignments/uri-schemes.html
% and http://en.wikipedia.org/wiki/URI_scheme#Generic_syntax
%
% Scheme Description Reference
% -------------------------------------------------------------
% aaa   Diameter Protocol   [RFC3588]
% aaas  Diameter Protocol with Secure Transport [RFC3588]
% acap  application configuration access protocol   [RFC2244]
% cap   Calendar Access Protocol    [RFC4324]
% cid   content identifier  [RFC2392]
% crid  TV-Anytime Content Reference Identifier [RFC4078]
% data  data    [RFC2397]
% dav   dav [RFC-ietf-webdav-rfc2518bis-18.txt]
% dict  dictionary service protocol [RFC2229]
% dns   Domain Name System  [RFC4501]
% fax   fax [RFC3966]
% file  Host-specific file names    [RFC1738]
% ftp   File Transfer Protocol  [RFC1738]
% go    go  [RFC3368]
% gopher    The Gopher Protocol [RFC4266]
% h323  H.323   [RFC3508]
% http  Hypertext Transfer Protocol [RFC2616]
% https Hypertext Transfer Protocol Secure  [RFC2818]
% icap  Internet Content Adaptation Protocol    [RFC3507]
% im    Instant Messaging   [RFC3860]
% imap  internet message access protocol    [RFC2192]
% info  Information Assets with Identifiers in Public Namespaces    [RFC4452]
% ipp   Internet Printing Protocol  [RFC3510]
% iris  Internet Registry Information Service   [RFC3981]
% iris.beep iris.beep   [RFC3983]
% iris.xpc  iris.xpc    [RFC-ietf-crisp-iris-xpc-06.txt]
% iris.xpcs iris.xpcs   [RFC-ietf-crisp-iris-xpc-06.txt]
% iris.lwz  iris.lwz    [RFC-ietf-crisp-iris-lwz-08.txt]
% ldap  Lightweight Directory Access Protocol   [RFC4516]
% mailto    Electronic mail address [RFC2368]
% mid   message identifier  [RFC2392]
% modem modem   [RFC3966]
% msrp  Message Session Relay Protocol  [RFC-ietf-simple-message-sessions-19.txt]
% msrps Message Session Relay Protocol Secure   [RFC-ietf-simple-message-sessions-19.txt]
% mtqp  Message Tracking Query Protocol [RFC3887]
% mupdate   Mailbox Update (MUPDATE) Protocol   [RFC3656]
% news  USENET news [RFC1738]
% nfs   network file system protocol    [RFC2224]
% nntp  USENET news using NNTP access   [RFC1738]
% opaquelocktoken   opaquelocktokent    [RFC-ietf-webdav-rfc2518bis-18.txt]
% pop   Post Office Protocol v3 [RFC2384]
% pres  Presence    [RFC3859]
% rtsp  real time streaming protocol    [RFC2326]
% service   service location    [RFC2609]
% shttp Secure Hypertext Transfer Protocol  [RFC2660]
% sip   session initiation protocol [RFC3261]
% sips  secure session initiation protocol  [RFC3261]
% snmp  Simple Network Management Protocol  [RFC4088]
% soap.beep soap.beep   [RFC3288]
% soap.beeps    soap.beeps  [RFC3288]
% tag   tag [RFC4151]
% tel   telephone   [RFC3966]
% telnet    Reference to interactive sessions   [RFC4248]
% tftp  Trivial File Transfer Protocol  [RFC3617]
% thismessage   multipart/related relative reference resolution [RFC2557]
% tip   Transaction Internet Protocol   [RFC2371]
% tv    TV Broadcasts   [RFC2838]
% urn   Uniform Resource Names (click for registry) [RFC2141]
% vemmi versatile multimedia interface  [RFC2122]
% xmlrpc.beep   xmlrpc.beep [RFC3529]
% xmlrpc.beeps  xmlrpc.beeps    [RFC3529]
% xmpp  Extensible Messaging and Presence Protocol  [RFC4622]
% z39.50r   Z39.50 Retrieval    [RFC2056]
% z39.50s   Z39.50 Session  [RFC2056]

author: PatternRecognition

repo: OpenBMI


text:
function [post,model] = spm_vb_roi_anova (VOI_fname,SPM,factor)
% Bayesian ANOVA for a region of interest
% FORMAT [post,model] = spm_vb_roi_anova (VOI_fname,SPM,factor)
%
% VOI_fname   - VOI filename
% SPM         - SPM data structure
% factor      - data structure relating conditions to levels of factors
%
% model       - data structure describing models
%               (m).F             model evidence
%               (m).X             design matrix
% post        - Posterior probabilities of
%               .factor1        main effect of factor 1
%               .factor2        main effect of factor 2
%               .interaction    interaction
%               .average        average
%__________________________________________________________________________
% Copyright (C) 2005-2014 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_vb_roi_anova.m 6079 2014-06-30 18:25:37Z spm $

if nargin < 2
    [Pf, sts] = spm_select(1,'^SPM\.mat$','Select SPM.mat');
    if ~sts, return; end
    swd = spm_file(Pf,'fpath');
    load(fullfile(swd,'SPM.mat'))
end

load(VOI_fname);

Y    = xY.y;
xyz  = xY.XYZmm;
N    = size(xyz,2);
M    = diag(SPM.xVol.M);
m    = 1./abs(M(1:3));
xyz  = (m*ones(1,N)).*xyz;
vxyz = spm_vb_neighbors(xyz',1);

%-Set number of AR coefficients
%--------------------------------------------------------------------------
try
    SPM.PPM.AR_P;
catch
    SPM.PPM.AR_P = 0;
end

%-Specify type of prior for regression coefficients
%--------------------------------------------------------------------------
try
    SPM.PPM.priors.W;
catch
    if N==1
        SPM.PPM.priors.W = 'Voxel - Shrinkage';
    else
        SPM.PPM.priors.W = 'Spatial - LORETA';
    end
end

%-Specify type of prior for AR coefficients
%--------------------------------------------------------------------------
try
    SPM.PPM.priors.A;
catch
    if N==1
        SPM.PPM.priors.A = 'Voxel - Shrinkage';
    else
        SPM.PPM.priors.A = 'Spatial - LORETA';
    end
end

% Get matrices that will remove low-frequency drifts
% if high pass filters have been specified
%--------------------------------------------------------------------------
s=1;
sess_nScan=length(SPM.xX.K(s).row);
if size(SPM.xX.K(s).X0,2) > 0
    X0=SPM.xX.K(s).X0;
    hpf(s).R0=eye(sess_nScan)-X0*pinv(X0);
else
    hpf(s).R0=eye(sess_nScan);
end

%-Filter data to remove low frequencies
R0Y = hpf(s).R0*Y(SPM.Sess(s).row,:);

%-Set optimisation parameters
%--------------------------------------------------------------------------
try
    SPM.PPM.maxits;
catch
    SPM.PPM.maxits = 16;
end
try
    SPM.PPM.tol;
catch
    SPM.PPM.tol = 0.00001;
end

%-Specify basis functions
%--------------------------------------------------------------------------
% SPM.xBF.name='hrf';
% SPM.xBF.order=1;
SPM.xBF.name='hrf (with time derivative)';
SPM.xBF.order=2;
SPM.xBF.length=32;
SPM.xBF = spm_get_bf(SPM.xBF);

nf=length(factor);

model = spm_vb_models (SPM,factor);

original_SPM = SPM;

%-Fit models
%--------------------------------------------------------------------------
for m=1:6

    if nf==2 || (nf==1 && (m==1 || m==2 || m==6))
        % fit model
        SPM = original_SPM;

        if ~(m==1 || m==6)
            % Get design matrix for relevant input set
            SPM.Sess(1).U=model(m).U;
            SPM.Sess(1).U=spm_get_ons(SPM,1);
            SPM=spm_fmri_design(SPM,0);     % Don't write SPM.mat file
            model(m).X=SPM.xX.X;
        end
        slice = spm_vb_init_volume (model(m).X,SPM.PPM.AR_P);

        slice.maxits=SPM.PPM.maxits;
        slice.tol=SPM.PPM.tol;
        slice.compute_det_D=1;
        slice.verbose=1;
        slice.update_w=1;
        slice.update_lambda=1;
        slice.update_F=1;
        slice = spm_vb_set_priors(slice,SPM.PPM.priors,vxyz);
        slice = spm_vb_glmar(R0Y,slice);

        model(m).F=slice.F;

        model(m).slice=slice;
    end
end

if nf==2

    F=[model(3).F,model(2).F];
    F=F-mean(F);
    post.factor1=exp(F(1))/sum(exp(F));

    F=[model(4).F,model(2).F];
    F=F-mean(F);
    post.factor2=exp(F(1))/sum(exp(F));

    F=[model(6).F,model(5).F];
    F=F-mean(F);
    post.interaction=exp(F(1))/sum(exp(F));

    F=[model(2).F,model(1).F];
    F=F-mean(F);
    post.average=exp(F(1))/sum(exp(F));

elseif nf==1

    F=[model(2).F,model(1).F];
    F=F-mean(F);
    post.average=exp(F(1))/sum(exp(F));

    F=[model(6).F,model(2).F];
    F=F-mean(F);
    post.factor1=exp(F(1))/sum(exp(F));

end


text:
function readAllDICOM(pathRead,pathSave,nBatch,nameSaveOption)
% -------------------------------------------------------------------------
% function readAllDICOM(pathRead,pathSave,batchFlag,nameSaveOption)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function reads the DICOM content of all the sub-folder tree of a
% starting directory defined by 'pathRead'. It then organizes the data
% (files throughout the starting directory are associated by
% 'SeriesInstanceUID') in a cell of structures called 'sData' (stands for
% STAMP data), and it finally computes the region of interest (ROI) defined
% by an associated RTstruct, and save as well all the  REG, RTdose and
% RTplan structures(if present in the sub-folder tree of the starting
% directory).
% All sData structures hereby created are saved in 'pathSave' with a name
% varying with the variable 'nameSaveOption'.
%
% DIFFERENTIATION/ASSOCIATION OF DICOM FILES: 1)imaging, 2)RTstruct, 3)REG, 4)RTdose, 5)RTplan.
% 1) Imaging volumes are differentiated by the 'SeriesInstanceUID' field
% 2) Association between a RTstruct and an imaging volume is performed with:
%    - 'ReferencedFrameOfReferenceSequence.Item_1.FrameOfReferenceUID' field, or
%    - 'FrameOfReferenceUID' field.
% 3) Association between a REG and an imaging volume is performed with:
%    - 'RegistrationSequence.Item_2.FrameOfReferenceUID' field, or
%    - 'FrameOfReferenceUID' field.
% 4) Association between a RTdose and an imaging volume is performed with:
%    - 'FrameOfReferenceUID' field.
% 5) Association between a RTplan and an imaging volume is performed with:
%    - 'FrameOfReferenceUID' field.
%
% If multiple imaging volumes have the same 'FrameOfReferenceUID' field,
% they will all be assigned a RTstruct, REG, RTdose and RTplan possessing
% that same field.
%
% IMPORTANT: This function has only been tested on Ubuntu 14.04. May work
%            using MAC, will fail using Windows.
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathRead: String specifying the full path to the starting directory
%              where the DICOM path to read are located.
%              --> Ex: '/home/myStudy/DICOM'
% 2. pathSave: String specifying the full path to the directory where to
%              save all the sData structures created by the current function.
%              --> Ex: '/home/myStudy/DATA'
% 3. nBatch: Numerical value specifying the number of batch to use in the
%            parallel computations (use 0 for serial).
%            --> Ex: 8
% 4. nameSaveOption: (optional). If this argument is not present (default)
%                    all sData structures will be saved with a filename
%                    defined by the 'Modality' field present in the
%                    DICOM header of their corresponding imaging volume. If
%                    multiple volumes of the same modality are present for
%                    the same 'PatientID', the different volumes will be
%                    numerated (e.g. CT1, CT2, CT3, etc.). If this argument
%                    is present and set to 'folder', all sData structures
%                    will be saved with a filename defined by the name of
%                    the folder in which all DICOM files of a given imaging
%                    volume are located (useful to keep track of meaningful
%                    names). The 'folder' option assumes that all imaging
%                    DICOM files of a given imaging volume are present in
%                    the same directory. Note: the 'folder' argument may
%                    lead to overwritting of data if not organized properly.
%                    Finally, if this argument is present and set to
%                    'modality',  all sData structures will be saved with a
%                    filename defined by the 'Modality' field in the DICOM
%                    header, but without enumerating multiple volumes of
%                    the same modality and same 'PatientID' (this may lead
%                    to overwritting, but is faster if the user is sure
%                    that only one scan of each modality for each patient
%                    is present in 'pathRead'.
%                    --> Options: - No argument
%                                 - 'folder'
%                                 - 'modality'
% -------------------------------------------------------------------------
% OUTPUTS: Multiple sData files
% - sData: Cell of structures organizing the content of the volume data,
%          DICOM headers, DICOM RTstruct* (used to compute the ROI) and
%          DICOM REGstruct* (used to register a given volume to another),
%          DICOM RTdose* and DICOM RTplan*.
%          * If present in the sub-folder tree
%    --> sData{1}: Explanation of cell content
%    --> sData{2}: Imaging data and ROI defintion
%    --> sData{3}: DICOM headers of imaging data
%    --> sData{4}: DICOM RTstruct (if applicable)
%    --> sData{5}: DICOM REGstruct (if applicable)
%    --> sData{6}: DICOM RTdose (if applicable)
%    --> sData{7}: DICOM RTplan (if applicable)
% -------------------------------------------------------------------------
% AUTHOR(S):
% - Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: February 2016
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
warning off

% INITIALIZATION
cd(pathRead)
stackPathRS = {}; stackFrameRS = {}; % Full path and FrameOfReferenceUID of the RTstruct files.
stackPathREG = {}; stackFrameREG = {}; % Full path and FrameOfReferenceUID of the REG files.
stackPathRD = {}; stackFrameRD = {}; % Full path and FrameOfReferenceUID of the RTdosefiles.
stackPathRP = {}; stackFrameRP = {}; % Full path and FrameOfReferenceUID of the RTplan files.
cellSeriesID = {}; % Cell of 'SeriesInstanceUID' of different imaging volumes (string).
cellFrameID = {}; % Cell of 'FrameOfReferenceUID' of different imaging volumes (string). Cell index is associated to the index of cellSeriesID.
cellPathImages = {}; % Cell of paths to the different imaging volumes. Cell index is associated to the index of cellSeriesID. Each cell in turn contains all the different path to the dicom images of a given volume.
cellPathRS = {}; % Cell of paths to the different RTstruct file (struct). Cell index is associated to the index of cellSeriesID.
cellPathREG = {}; % Cell of paths to the different REG file (struct). Cell index is associated to the index of cellSeriesID.
cellPathRD = {}; % Cell of paths to the different RTdose file (struct). Cell index is associated to the index of cellSeriesID.
cellPathRP = {}; % Cell of paths to the different RTplan file (struct). Cell index is associated to the index of cellSeriesID.
nameSave = {}; % Cell of filenames for all created sData structures (string). Cell index is associated to the index of cellSeriesID.
fidErr = fopen('readAllDICOM_err','w');
fidRead = fopen('readAllDICOM_log','w');

% SCANNING ALL FOLDERS IN INITIAL DIRECTORY
stackFolder = {pathRead}; % Stack of folder paths to scan
ind = strfind(stackFolder{1},'/'); nameDir = stackFolder{1}((ind(end)+1):end);
indUnder = strfind(nameDir,'_');
if ~isempty(indUnder)
    for under = numel(indUnder)
        nameDir = [nameDir(1:indUnder(under)-1),'\',nameDir(indUnder(under):end)];
    end
end
while ~isempty(stackFolder)
    cd(stackFolder{1})
    fprintf(fidRead,['\nScanning files in directory: ',stackFolder{1},' ... ']);
    list = dir; nElem = numel(list);
    for i = 1:nElem
        if ~strcmp(list(i).name,'.') && ~strcmp(list(i).name,'..') && isempty(strfind(list(i).name,'DS_Store')) % Good enough for linux, maybe add other conditions for MAC and Windows
            if list(i).isdir
                cd(list(i).name), stackFolder = [stackFolder;pwd]; cd ..
            elseif isdicom(list(i).name)
                try
                    info = dicominfo(list(i).name);
                    if strcmp(info.Modality,'MR') || strcmp(info.Modality,'PT') || strcmp(info.Modality,'CT')
                        [indSeriesID] = findUIDcellIndex(info.SeriesInstanceUID,cellSeriesID);
                        if indSeriesID > numel(cellSeriesID) % New volume
                            cellSeriesID = [cellSeriesID cell(1)]; cellSeriesID{indSeriesID} = info.SeriesInstanceUID;
                            cellFrameID = [cellFrameID cell(1)]; cellFrameID{indSeriesID} = info.FrameOfReferenceUID;
                            cellPathImages = [cellPathImages cell(1)]; cellPathImages{indSeriesID} = {};
                            cellPathRS = [cellPathRS cell(1)]; cellPathRS{indSeriesID} = {};
                            cellPathREG = [cellPathREG cell(1)]; cellPathREG{indSeriesID} = {};
                            cellPathRD = [cellPathRD cell(1)]; cellPathRD{indSeriesID} = {};
                            cellPathRP = [cellPathRP cell(1)]; cellPathRP{indSeriesID} = {};
                            nameSave = [nameSave cell(1)];
                            if nargin > 3 && strcmp(nameSaveOption,'folder')
                                path = pwd;
                                ind = strfind(path,'/'); % NOT GOOD FOR WINDOWS
                                nameSave{indSeriesID} = path((ind(end)+1):end);
                            elseif nargin > 3 && strcmp(nameSaveOption,'modality')
                                nameSave{indSeriesID} = info.Modality;
                            else
                                try
                                    nameSave{indSeriesID} = info.SeriesDescription;
                                catch
                                    nameSave{indSeriesID} = info.SeriesInstanceUID;
                                end
                            end
                        end
                        cellPathImages{indSeriesID} = [cellPathImages{indSeriesID};fullfile(pwd,list(i).name)];
                    elseif strcmp(info.Modality,'RTSTRUCT')
                        stackPathRS = [stackPathRS;fullfile(pwd,list(i).name)];
                        try
                            frameUID = info.ReferencedFrameOfReferenceSequence.Item_1.FrameOfReferenceUID;
                        catch
                            frameUID = info.FrameOfReferenceUID;
                        end
                        stackFrameRS = [stackFrameRS;frameUID];
                    elseif strcmp(info.Modality,'REG')
                        stackPathREG = [stackPathREG;fullfile(pwd,list(i).name)];
                        try
                            frameUID = info.RegistrationSequence.Item_2.FrameOfReferenceUID;
                        catch
                            frameUID = info.FrameOfReferenceUID;
                        end
                        stackFrameREG = [stackFrameREG;frameUID];
                    elseif strcmp(info.Modality,'RTDOSE')
                        stackPathRD = [stackPathRD;fullfile(pwd,list(i).name)];
                        frameUID = info.FrameOfReferenceUID;
                        stackFrameRD = [stackFrameRD;frameUID];
                    elseif strcmp(info.Modality,'RTPLAN')
                        stackPathRP = [stackPathRP;fullfile(pwd,list(i).name)];
                        frameUID = info.FrameOfReferenceUID;
                        stackFrameRP = [stackFrameRP;frameUID];
                    end
                catch
                    fprintf(fidErr,'Error while reading: %s\n',fullfile(pwd,list(i).name));
                end
            end
        end
    end
    stackFolder(1) = [];
    fprintf(fidRead,'DONE!');
end
fclose(fidErr);
fprintf(fidRead,'\n');

% ASSOCIATING ALL RTSTRUCT TO IMAGING VOLUMES
nRS = numel(stackPathRS);
if nRS
    fprintf(fidRead,'\nASSOCIATING ALL RTSTRUCT TO IMAGING VOLUMES ... ');
    for i = 1:nRS
        [indSeriesID] = findUIDcellIndex(stackFrameRS{i},cellFrameID);
        for n = 1:numel(indSeriesID)
            cellPathRS{indSeriesID(n)} = [cellPathRS{indSeriesID(n)};stackPathRS{i}];
        end
    end
    fprintf(fidRead,'DONE!');
end

% ASSOCIATING ALL REG TO IMAGING VOLUMES
nREG = numel(stackPathREG);
if nREG
    fprintf(fidRead,'\nASSOCIATING ALL REG TO IMAGING VOLUMES ... ');
    for i = 1:nREG
        [indSeriesID] = findUIDcellIndex(stackFrameREG{i},cellFrameID);
        for n = 1:numel(indSeriesID)
            cellPathREG{indSeriesID(n)} = [cellPathREG{indSeriesID(n)};stackPathREG{i}];
        end
    end
    fprintf(fidRead,'DONE!');
end

% ASSOCIATING ALL RTdose TO IMAGING VOLUMES
nRD = numel(stackPathRD);
if nRD
    fprintf(fidRead,'\nASSOCIATING ALL RTdose TO IMAGING VOLUMES ... ');
    for i = 1:nRD
        [indSeriesID] = findUIDcellIndex(stackFrameRD{i},cellFrameID);
        for n = 1:numel(indSeriesID)
            cellPathRD{indSeriesID(n)} = [cellPathRD{indSeriesID(n)};stackPathRD{i}];
        end
    end
    fprintf(fidRead,'DONE!');
end

% ASSOCIATING ALL RTplan TO IMAGING VOLUMES
nRP = numel(stackPathRP);
if nRP
    fprintf(fidRead,'\nASSOCIATING ALL RTplan TO IMAGING VOLUMES ... ');
    for i = 1:nRP
        [indSeriesID] = findUIDcellIndex(stackFrameRP{i},cellFrameID);
        for n = 1:numel(indSeriesID)
            cellPathRP{indSeriesID(n)} = [cellPathRP{indSeriesID(n)};stackPathRP{i}];
        end
    end
    fprintf(fidRead,'DONE!');
end

% READING ALL IMAGES TO CREATE ALL sData STRUCTURES
nScans = numel(cellPathImages);
if nBatch
    if nScans < nBatch
        nBatch = nScans;
    end
    fprintf(fidRead,['\n\nREADING AND PROCESSING ALL DICOM DATA USING ',num2str(nBatch),' PARALLEL BATCH ... ']);
    cd(pathRead), mkdir('BATCH_LOG'), cd('BATCH_LOG'), pathBatch = pwd;
    [scans] = batchPatients(nScans,nBatch); % Scans for each batch
    save('workspace','scans','pathSave','cellPathImages','cellPathRS','cellPathREG','cellPathRD','cellPathRP','nameSave'), pause(2)
    for i = 1:nBatch
        nScan = numel(scans{i});
        nameScript = ['batch',num2str(i),'_script.m'];
        fid = fopen(nameScript,'w');
        fprintf(fid,'dummy = 0;\n');
        fprintf(fid,'load(''workspace'')\n');
        for j = 1:nScan
            fprintf(fid,['sDataCreation_FromDICOMpaths(pathSave,cellPathImages{scans{',num2str(i),'}(',num2str(j),')},cellPathRS{scans{',num2str(i),'}(',num2str(j),')},cellPathREG{scans{',num2str(i),'}(',num2str(j),')},cellPathRD{scans{',num2str(i),'}(',num2str(j),')},cellPathRP{scans{',num2str(i),'}(',num2str(j),')},nameSave{scans{',num2str(i),'}(',num2str(j),')})\n']);
        end
        fprintf(fid,['system(''touch batch',num2str(i),'_end'');\n']);
        fprintf(fid,'clear all');
        fclose(fid);
        system(['matlab -nojvm -nodisplay -nodesktop -nosplash < ',nameScript,' >& ',nameScript(1:end-1),'log &']);
    end
    waitBatch(pathBatch,30,nBatch)
    cd ..
    %system('rm -r BATCH_LOG');
else
    fprintf(fidRead,'\n\nREADING AND PROCESSING ALL DICOM DATA ... ');
    for i = 1:nScans
        sDataCreation_FromDICOMpaths(pathSave,cellPathImages{i},cellPathRS{i},cellPathREG{i},cellPathRD{i},cellPathRP{i},nameSave{i});
    end
end
fprintf(fidRead,'DONE!\n');

% POST-PROCESSING OF SAVED FILES IF 4th ARGUMENT IS SET TO 'modality'
if nargin < 4
    fprintf(fidRead,'\nCORRECTING NAME OF SAVED FILES ... ');
    type = {'.CTscan','.PTscan','.MRscan'}; nType = numel(type);
    cd(pathSave)
    list = dir; nElem = numel(list);
    cellNames = {};
    for i = 1:nElem
        if ~list(i).isdir
            cellNames = [cellNames,list(i).name];
        end
    end
    nNames = numel(cellNames);
    temp = '';
    for i = 1:nNames
        try
            load(cellNames{i})
            idString = sData{3}(1).PatientID;
            if strcmp(idString,temp)
                ok = 0;
            else
                ok = 1;
            end
        catch
            ok = 0;
        end
        if ok
            cellTemp = strfind(cellNames,idString);
            for c = 1:numel(cellTemp)
                if isempty(cellTemp{c})
                    cellTemp{c} = [0];
                end
            end
            idIndex = find(cell2mat(cellTemp));
            nIndex = numel(idIndex); nameComp = cell(1,nIndex);
            for j = 1:nIndex
                nameComp{j} = cellNames{idIndex(j)};
            end
            for j = 1:nType
                cellTemp = strfind(nameComp,type{j});
                for c = 1:numel(cellTemp)
                    if isempty(cellTemp{c})
                        cellTemp{c} = [0];
                    end
                end
                scanIndex = find(cell2mat(cellTemp));
                nScanIndex = numel(scanIndex);
                if nScanIndex > 0
                    for k = 1:nScanIndex
                        ind1 = numel(idString) + 1; % Getting the '_' position after the PatientID
                        ind2 = numel(nameComp{scanIndex(k)}) - 10; % Getting the '.' position
                        name = nameComp{scanIndex(k)};
                        indWrong = strfind(name,'\');
                        if ~isempty(indWrong)
                            nInd = length(indWrong);
                            for n = 1:nInd
                                name = [name(1:indWrong(n)+n-2),'\',name(indWrong(n)+n-1:end)];
                            end
                        end
                        indWrong = strfind(name,'(');
                        if ~isempty(indWrong)
                            nInd = length(indWrong);
                            for n = 1:nInd
                                name = [name(1:indWrong(n)+n-2),'\',name(indWrong(n)+n-1:end)];
                            end
                        end
                        indWrong = strfind(name,')');
                        if ~isempty(indWrong)
                            nInd = length(indWrong);
                            for n = 1:nInd
                                name = [name(1:indWrong(n)+n-2),'\',name(indWrong(n)+n-1:end)];
                            end
                        end
                        indWrong = strfind(name,'&');
                        if ~isempty(indWrong)
                            nInd = length(indWrong);
                            for n = 1:nInd
                                name = [name(1:indWrong(n)+n-2),'\',name(indWrong(n)+n-1:end)];
                            end
                        end
                        indWrong = strfind(name,' ');
                        if ~isempty(indWrong)
                            nInd = length(indWrong);
                            for n = 1:nInd
                                name = [name(1:indWrong(n)+n-2),'\',name(indWrong(n)+n-1:end)];
                            end
                        end
                        if nScanIndex == 1
                            system(['mv ',name,' ',nameComp{scanIndex(k)}(1:ind1(1)),type{j}(2:3),nameComp{scanIndex(k)}(ind2(1):end)]);
                        else
                            system(['mv ',name,' ',nameComp{scanIndex(k)}(1:ind1(1)),type{j}(2:3),num2str(k),nameComp{scanIndex(k)}(ind2(1):end)]);
                        end
                    end
                end
            end
            temp = idString;
        end
    end
    fprintf(fidRead,'DONE!\n');
end
fclose(fidRead);

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function input_info=my_init_input_info(input_info)

if input_info.is_group_data
   input_info_child_groups=input_info.data_child_groups;
   child_num=length(input_info_child_groups);
   for c_idx=1:child_num
       input_info_child_groups{c_idx}=my_init_input_info(input_info_child_groups{c_idx});
   end
   input_info.data_child_groups=input_info_child_groups;
   input_info.child_valid_flags=true(child_num, 1);
end

input_info.bp_finished=false;
input_info.forward_finished=true;

input_info.aux=[];
input_info.dzdx=[];
input_info.dzdw=[];

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
% This function calculates the ground truth images for a given dataset by
% averaging multiple consecutive frames at the stop-motion time steps.
function groundTruthSequence = createGroundTruthImages(datasetDir, type)

    % Get all files in the dataset directory.
    files = dir([datasetDir, '/*.', type]);
    numFiles = length(files);
    fileIdx = 1;
    groundTruthIdx = 1;
    while fileIdx <= numFiles
        % Read the first frame of the k-th time step.
        filename = files(fileIdx).name;
        frames(:,:,1) = im2double( imread([datasetDir, '/', filename]) );

        % Read all remaining frames of the k-th time step.
        for frameIdx = 1:30
            % Get filename and extract the frame number.
            if fileIdx + frameIdx > numFiles
                % We reached the end of the file list.
                break;
            end
            filename = files(fileIdx + frameIdx).name;
            [~, filenameWithoutEx, ~] = fileparts(filename);
            frameNumber = str2num( filenameWithoutEx(end-1:end) ); %#ok<ST2NM>

            if frameNumber ~= 0
                % New frame for the current time step.
                frames(:,:,frameIdx+1) = im2double( imread([datasetDir, '/', filename]) ); %#ok<AGROW>
            else
                % We reached the last frame of the current time step.
                break;
            end
        end

        % Calculate ground truth image of the k-th time step by the
        % temporal mean (or median) of all corresponding frames.
        groundTruthSequence(:,:,groundTruthIdx) = mean(frames, 3); %#ok<AGROW>

        fileIdx = fileIdx + frameIdx;
        groundTruthIdx = groundTruthIdx + 1;
    end

author: thomas-koehler

repo: SupER