text:
% QP test problem. Demo file
%
clc;

disp('----------------------------------');
disp('Problem #1:');
disp(' ');
disp('min x1^2 + x^2  - 3*x1 - 1/4*x2');
disp(' ');
disp('subject to');
disp(' |x1|+|x2| <= 1');
disp('----------------------------------');

H=2*eye(2);
q=[-3,-1/4]';
A=[1 1;1 -1;-1 1;-1 -1];
b=ones(4,1);

[x, obj, lambda, info] = qpng (H, q, A, b)

pause;
clc

disp('----------------------------------');
disp('Problem #2:');
disp(' ');
disp('min 4*x1^2 + 1 x1*x2 + 4*x2^2 + 3*x1 ? 4*x2');
disp(' ');
disp('s.t. x1+x2 <= 5');
disp('     x1-x2 == 0');
disp('     x1 >= 0');
disp('     x2 >= 0');
disp('----------------------------------');

H=[8 1;1 8];
q=[3 -4]';
A=[1 1; 1 -1];
b=[5 0]';
ctype=char(['U', 'E']');
lb=[0 0]';
ub=[];

[x, obj, lambda, info] = qpng (H, q, A, b,ctype,lb,ub)

pause;
clc

disp('----------------------------------');
disp('Problem #3:');
disp(' ');
disp('min .5*x1''*H*x + q''*x');
disp(' ');
disp('s.t. -2.5 <= x1 <= 2.5');
disp('----------------------------------');

H=[.16 -1.2 2.4 -1.4;
   -1.2 12.0 -27.0 16.8;
   2.4 -27.0 64.8 -42.0;
   -1.4 16.8 -42.0 28.0];
q=[5.04 -59.4 146.4 -96.6]';
A=[-1 0 0 0; 1 0 0 0];
b=[2.5; 2.5];

[x, obj, lambda, info] = qpng (H, q, A, b)

pause;
clc

disp('----------------------------------');
disp('Problem #4:');
disp(' Solve a portfolio example. ');
disp(' See Section 13.7 of "Applications of optimization with Xpress-MP.');
disp('----------------------------------');

H=[4 3 -1 0;
   3 6 1 0;
   -1 1 10 0;
   0 0 0 10];
q=[];
A=[1 1 1 1;
   8 9 12 7];
b=[1; 7];
ctype=char(['E', 'L']');
lb=zeros(4,1);
ub=[];

[x, obj, lambda, info] = qpng (H, q, A, b, ctype, lb, ub)

pause;
clc

disp('----------------------------------');
disp('Problem #5:');
disp(' An indefinite QP problem is not solved x*=(3, 25/4). ');
disp(' ');
disp('min .5*x1^2 - .5*x2^2');
disp(' ');
disp('s.t. x2 >= 0');
disp('s.t. x1+x2 >= 2');
disp('s.t. -5x1+4x2 <= 10');
disp('s.t. x1 <= 3');
disp('----------------------------------');

H=[1 0;0 -1];
q=[];
A=[0 1;1 2;-5 4;1 0];
b=[0 2 10 3]';
ctype=char(['L', 'L', 'U', 'U']);

[x, obj, lambda, info]=qpng(H, q, A, b)

author: opencobra

repo: cobratoolbox


text:
function [U_final, V_final, nIter_final, objhistory_final] = GNMF(X, k, W, options, U, V)
% Graph regularized Non-negative Matrix Factorization (GNMF)
%
% where
%   X
% Notation:
% X ... (mFea x nSmp) data matrix
%       mFea  ... number of words (vocabulary size)
%       nSmp  ... number of documents
% k ... number of hidden factors
% W ... weight matrix of the affinity graph
%
% options ... Structure holding all settings
%               options.alpha ... the regularization parameter.
%                                 [default: 100]
%                                 alpha = 0, GNMF boils down to the ordinary NMF.
%
%
% You only need to provide the above four inputs.
%
% X = U*V'
%
% References:
% [1] Deng Cai, Xiaofei He, Xiaoyun Wu, and Jiawei Han. "Non-negative
% Matrix Factorization on Manifold", Proc. 2008 Int. Conf. on Data Mining
% (ICDM'08), Pisa, Italy, Dec. 2008.
%
% [2] Deng Cai, Xiaofei He, Jiawei Han, Thomas Huang. "Graph Regularized
% Non-negative Matrix Factorization for Data Representation", IEEE
% Transactions on Pattern Analysis and Machine Intelligence, , Vol. 33, No.
% 8, pp. 1548-1560, 2011.
%
%
%   version 2.0 --April/2009
%   version 1.0 --April/2008
%
%   Written by Deng Cai (dengcai AT gmail.com)
%

if min(min(X)) < 0
    error('Input should be nonnegative!');
end

if ~isfield(options,'error')
    options.error = 1e-5;
end
if ~isfield(options, 'maxIter')
    options.maxIter = [];
end

if ~isfield(options,'nRepeat')
    options.nRepeat = 10;
end

if ~isfield(options,'minIter')
    options.minIter = 30;
end

if ~isfield(options,'meanFitRatio')
    options.meanFitRatio = 0.1;
end

if ~isfield(options,'alpha')
    options.alpha = 100;
end

nSmp = size(X,2);

if isfield(options,'alpha_nSmp') && options.alpha_nSmp
    options.alpha = options.alpha*nSmp;
end

if isfield(options,'weight') && strcmpi(options.weight,'NCW')
    feaSum = full(sum(X,2));
    D_half = X'*feaSum;
    X = X*spdiags(D_half.^-.5,0,nSmp,nSmp);
end

if ~isfield(options,'Optimization')
    options.Optimization = 'Multiplicative';
end

if ~exist('U','var')
    U = [];
    V = [];
end

switch lower(options.Optimization)
    case {lower('Multiplicative')}
        [U_final, V_final, nIter_final, objhistory_final] = GNMF_Multi(X, k, W, options, U, V);
    otherwise
        error('optimization method does not exist!');
end

author: ZJULearning

repo: MatlabFunc


text:
function regular_training(video, method, frames)

opts.expDir = ['net/' method '/' num2str(frames) '/' video] ;

opts.train.batchSize = 5 ;
opts.train.numEpochs = 20;
opts.train.continue = false ;
opts.train.useGpu = true ;
opts.train.learningRate    = 1e-3;
opts.train.expDir = opts.expDir ;

% --------------------------------------------------------------------
%                                                         Prepare data
% --------------------------------------------------------------------

imgDir = ['../CDNetDataset/' method '/' num2str(frames) 'frames/' video '/input'];
labelDir = ['../CDNetDataset/' method '/' num2str(frames) 'frames/' video '/GT'];
imdb = getImdb(imgDir,labelDir);

mask = imread(['../CDNetDataset/' method '/' num2str(frames) 'frames/' video '/ROI.bmp']);
mask = mask(:,:,1);
A = max(max(mask));
mask(mask == A) = 1;
if size(mask,1) > 400 || size(mask,2) >400
    mask = imresize(mask, 0.5, 'nearest');
end
imdb.mask = single(double(mask));

imdb.half_size = 15;

%%%%%%Yi%%%%%% redefined the net
load('net');

net.layers{end-1} = struct('type', 'conv', ...
    'filters', 0.1*randn(1,1,64,1, 'single'), ...
    'biases', zeros(1, 1, 'single'), ...
    'stride', 1, ...
    'pad', 0) ;
net.layers{end} = struct('type', 'sigmoidcrossentropyloss');

load('meanPixel.mat');
imdb.meanPixel = meanPixel;

[net,info] = cnn_train_adagrad(net, imdb, @getBatch,...
    opts.train,'errorType','euclideanloss',...
    'conserveMemory', true);

end

function [im, labels] = getBatch(imdb, batch)
% --------------------------------------------------------------------

half_size = imdb.half_size;
meanPixel = imdb.meanPixel;

for ii = 1:numel(batch)
    imagename = imdb.images.name{batch(ii)};
    im_ii = single(imread(imagename));

    labelname = imdb.images.labels{batch(ii)};
    roi = imread(labelname);
    labels_ii = zeros(size(roi,1),size(roi,2));
    labels_ii( roi == 50 )  = 0.25;       %shade
    labels_ii( roi == 170 ) = 0.75;       %object boundary
    labels_ii( roi == 255 ) = 1;          %foreground

    % resize the image to half size
    if size(im_ii,1) > 400 || size(im_ii,2) >400
        im_ii = imresize(im_ii, 0.5, 'nearest');
        labels_ii = imresize(labels_ii, 0.5, 'nearest');
    end
    im_large = padarray(im_ii,[half_size,half_size],'symmetric');
    im_ii = bsxfun(@minus, im_large, meanPixel);

    im(:,:,:,ii) = im_ii;
    labels(:,:,1,ii) = labels_ii;
    labels(:,:,2,ii) = double(imdb.mask);

end
end

function imdb = getImdb(imgDir, labelDir)

files = dir([imgDir '/*.jpg']);
label_files = dir([labelDir '/*.png']);
names = {};labels = {};

for ii = 1:numel(files)
    names{end+1} = [imgDir '/' files(ii).name];
    labels{end+1} = [labelDir '/' label_files(ii).name];
end

imdb.images.set = ones(1,numel(names));
imdb.images.name = names ;
imdb.images.labels = labels;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function r8_erf_test ( )

%*****************************************************************************80
%
%% R8_ERF_TEST tests R8_ERF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 August 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_ERF_TEST:\n' );
  fprintf ( 1, '  R8_ERF evaluates the error function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     X      Exact F   R8_ERF(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = erf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r8_erf ( x );

    fprintf ( 1, '  %6f  %12f  %12f\n', x, fx, fx2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function geometry_test047 ( )

%*****************************************************************************80
%
%% TEST047 tests MINQUAD.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST047\n' );
  fprintf ( 1, '  MINQUAD finds the minimum of a function\n' );
  fprintf ( 1, '    F(X) = A * X * X + B * X + C\n' );
  fprintf ( 1, '    within an interval, given three data points.\n' );
%
%  Case 1: a minimum is in the interval.
%  y = ( x - 1 ).^2 + 4
%
  x1 = 0.0;
  y1 = ( x1 - 1.0 ).^2 + 4.0;

  x2 = 2.0;
  y2 = ( x2 - 1.0 ).^2 + 4.0;

  x3 = 3.0;
  y3 = ( x3 - 1.0 ).^2 + 4.0;

  [ xmin, ymin ] = minquad ( x1, y1, x2, y2, x3, y3 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The minimum lies in the interval.\n' );
  fprintf ( 1, '  X1,   Y1   = %8f  %8f\n', x1, y1 );
  fprintf ( 1, '  X2,   Y2   = %8f  %8f\n', x2, y2 );
  fprintf ( 1, '  X3,   Y3   = %8f  %8f\n', x3, y3 );
  fprintf ( 1, '  XMIN, YMIN = %8f  %8f\n', xmin, ymin );
%
%  Case 2: the minimum is to the left of the interval.
%  y = ( x - 1 ).^2 + 4
%
  x1 = 2.0;
  y1 = ( x1 - 1.0 ).^2 + 4.0;

  x2 = 4.0;
  y2 = ( x2 - 1.0 ).^2 + 4.0;

  x3 = 5.0;
  y3 = ( x3 - 1.0 ).^2 + 4.0;

  [ xmin, ymin ] = minquad ( x1, y1, x2, y2, x3, y3 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The minimum is to the left of the interval\n' );
  fprintf ( 1, '  X1,   Y1   = %8f  %8f\n', x1, y1 );
  fprintf ( 1, '  X2,   Y2   = %8f  %8f\n', x2, y2 );
  fprintf ( 1, '  X3,   Y3   = %8f  %8f\n', x3, y3 );
  fprintf ( 1, '  XMIN, YMIN = %8f  %8f\n', xmin, ymin );
%
%  Case 3: the function is flat.
%
  x1 = 11.0;
  y1 = 6.0;

  x2 = 6.0;
  y2 = 6.0;

  x3 = 2.0;
  y3 = 6.0;

  [ xmin, ymin ] = minquad ( x1, y1, x2, y2, x3, y3 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The function is flat.\n' );
  fprintf ( 1, '  X1,   Y1   = %8f  %8f\n', x1, y1 );
  fprintf ( 1, '  X2,   Y2   = %8f  %8f\n', x2, y2 );
  fprintf ( 1, '  X3,   Y3   = %8f  %8f\n', x3, y3 );
  fprintf ( 1, '  XMIN, YMIN = %8f  %8f\n', xmin, ymin );
%
%  Case 4: the function has a maximum.
%  y = - ( x - 1 ).^2 + 4
%
  x1 = 0.0;
  y1 = - ( x1 - 1.0 ).^2 + 4.0;

  x2 = 2.0;
  y2 = - ( x2 - 1.0 ).^2 + 4.0;

  x3 = 3.0;
  y3 = - ( x3 - 1.0 ).^2 + 4.0;

  [ xmin, ymin ] = minquad ( x1, y1, x2, y2, x3, y3 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The function has a maximum.\n' );
  fprintf ( 1, '  X1,   Y1   = %8f  %8f\n', x1, y1 );
  fprintf ( 1, '  X2,   Y2   = %8f  %8f\n', x2, y2 );
  fprintf ( 1, '  X3,   Y3   = %8f  %8f\n', x3, y3 );
  fprintf ( 1, '  XMIN, YMIN = %8f  %8f\n', xmin, ymin );

  return
end


text:
function pass = test_spherefun( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e2*pref.techPrefs.chebfuneps;

% Example 1
f = ballfun(@(x,y,z)cos(x.*y));
g = spherefun(f);
h = spherefun(@(x,y,z)cos(x.*y));
pass(1) = norm(g-h) < tol;

% Example 2
f = ballfun(@(x,y,z)sin(z));
g = spherefun(f, 0.5);
h = spherefun(@(x,y,z)sin(0.5*z));
pass(2) = norm(g-h) < tol;

% Example 3
f = ballfun(@(x,y,z)x);
g = f(10,:,:,'spherical');
h = spherefun(@(x,y,z)10*x);
pass(3) = norm(g-h) < tol;

if (nargout > 0)
    pass = all(pass(:));
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%kmexchg 'Exchange Rows or Columns of a Matrix'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros mexchg.pane file
%
% Parameters:
% InputFile: i 'Input Matrix', required: 'Input matrix'
% OutputFile: o 'Output Matrix', required: 'Output matrix'
% Integer: r1 'Row #1 ', default: 1: 'row #1'
% Integer: r2 'Row #2 ', default: 1: 'row #2'
% Integer: c1 'Column #1 ', default: 1: 'Column #1'
% Integer: c2 'Column #2 ', default: 1: 'Column #2'
%
% Example: o = kmexchg(i, {'i','';'o','';'r1',1;'r2',1;'c1',1;'c2',1})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% mexchg - Exchange Rows or Columns of a Matrix
%
%  DESCRIPTION
% .I mexchg
% is used to interchage a pair of columns and/or a pair of rows of a matrix.
% If only a row exchange is desired, set the column indices equal to each other.
% Likewise, if only a column exhange is desired, set the row indices equal to each
% other.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [ref,listRxnColour] = addColour(parsed,listRxn_Color,colorScheme)
% Changes colour attributes of the reaction links in a parsed `CellDesigner`
% model structure given a list of reaction IDs
%
% USAGE:
%
%    [ref, listRxnColour] = addColour(parsed, listRxn_Color, colorScheme)
%
% INPUTS:
%    parsed:          A parsed model structure generated by `parseCD` function
%    listRxn_Color:   A list of reaction IDs that need to be highilighted by
%                     changing the colour attributes of the reaciton links in
%                     the `CellDesigner` model. The first column stores a list
%                     of reaction IDs whose reaction links need to be
%                     highlighted, whereas the second column saves a list of
%                     html Colours.
%
% OUTPUT:
%
%    ref:             An updated parsed `CellDesigner` model structure
%    listRxnColour:
%
% .. Author: - Longfei Mao October 2014

if nargin<3

    colorScheme=2;

end

ref=parsed;
listRxn=listRxn_Color(:,1); % coloumn 1: reaction IDs;
a=size(listRxn_Color)
if length(a)<2
    listColor=[];
    listColor=listRxn_Color(:,2); % coloumn 2: html colour codes;
    p=1;
else
    p=0;
end

cScheme={{'#f16359';'#f18d59';'#f9ea76';'#acf3be';'#117e9a'}; % 'Strawberry' % Strawberry Orchard Color Palette
    {'#def595';'#cae45b';'#81c828';'#3dc800';'#279700'};  % 'Spring' Spring Growth Color Palette
    {'#ffedad';'#ffd599';'#ffc471';'#ffb93f';'#ffac37'};  % 'Tangerine' Tangerine Color Palette
    {'#00f9ff';'#009fff';'#0078ff';'#005eff';'#0900ff'};  % Blue Vibes
    {'#fe8181';'#fe5757';'#fe2e2e';'#cb2424';'#b62020'}}  % Red pt1

hexCode={};

for i=1:size(cScheme,1);
    for m=1:size(cScheme{1},1)
        if length(cScheme{i}{m})<9;
            hexCode{i}{m}=[cScheme{i}{m}(1),'ff',cScheme{i}{m}(2:end)]
        end
    end
end

num=0;

%results=[];

[ID_row,ID_Column]=size(ref.r_info.ID);

for m=1:ID_row;
    for n=1:ID_Column;
        r=iscellstr(ref.r_info.ID(m,n));
        if ~r;
            %results(or,1)=m;
            %results(or,2)=n;
            ref.r_info.ID{m,n}=' '
        end;
    end;
end
% listRxnColour=listRxn;
column=[];
for r=1:length(listRxn);

%     if numCol>0;
%
%         indRxns=find(ismember(ref.r_info.ID(:,numCol),listRxn{r}));
%         newRxnName=ref.r_info.ID(indRxns,2);
%
%     else

        newRxnName=listRxn{r};
        id=find(ismember(ref.r_info.ID,listRxn(r)))
        if numel(id)>1;
            id=id(1);
        end
% %         if id
% %             [m,n]=size(ref.r_info.ID);
% %             if id>m*(n-1)  % the third column of ref.r_info.ID contains reaction ID; for example: {'re5160','re8','DESAT16_2'}
% %                 newRxnName=ref.r_info.ID{id-2*m}; % newRxnName is defined to be the ID in the first column of ref.r_info.ID.
% %                 column=3;
% %             end
% %         end
        found=0;
        if id

            for rr=1:n;
                if id<m
                    newRxnName=ref.r_info.ID{id};
                    break;
                elseif and(id>rr*m,id<(rr+1)*m);
                    newRxnName=ref.r_info.ID{id-rr*m};
                    break;
                end
            end
        else
            if ~isfield(ref,listRxn{r})
                newRxnName=strcat('R_',listRxn{r});
                if  isempty(strfind(newRxnName,'(e)'))
                    newRxnName=strrep(newRxnName,'(e)','_e');
                end
            end
        end

    if ~isfield(ref,newRxnName)
        disp(listRxn{r});
        fprintf('error ! the listRxn{%d}',r);
        r=r+1
    else
        [rw,cw]=size(ref.(newRxnName).color)
        for  ddr=1:rw
            if ischar(ref.(newRxnName).width{ddr,1}) % if it is char such as ('1.0');then convert it into double;
                ref.(newRxnName).width{ddr,1}=str2double(ref.(newRxnName).width{ddr,1})
            end
            %                 try
            w(1)=ref.(newRxnName).width{ddr,1}
            %             catch
            % %                 disp(w(1))
            %                  disp(newRxnName);
            %                 disp(ref.(newRxnName).width{ddr,1});
            %             end

            %disp(ref.(newRxnName).width(ddr,1));disp('dddd');
            %disp(w(1));
            %              if isempty(w)
            %                  w=0;
            %              end            %
            fprintf('w value is %d\n',w(1));
            if w(1)>=0&&w(1)<=10;  % flux ranges from 2 to 10 will be highligthed
                if p==0;
                    %                     the colours (Hex triplet, e.g.,
                    %                     https://closedxml.codeplex.com/wikipage?title=Excel%20Indexed%20Colors)
                    %% 'FF' as in '#FF0000FF' is the prefix

                    if  w(1)==0
                        colorStr='#ffaed5fc';
                    elseif w(1)>0&&w(1)<=3;
                        colorStr=hexCode{colorScheme}{1} % '#FF0000FF' % flux between 1 and 2 will be highlighted in blue.
                    elseif w(1)>2&&w(1)<=4
                        colorStr=hexCode{colorScheme}{2} % flux between 2 and 5 will be highlighted in pink.

                    elseif w(1)>3&&w(1)<=5;
                        colorStr=hexCode{colorScheme}{3}

                    elseif w(1)>4&&w(1)<=6;
                        colorStr=hexCode{colorScheme}{4}
                    elseif w(1)>6;
                        colorStr=hexCode{colorScheme}{5} % flux between 2 and 5 will be highlighted in red.
                    end
                    if w(1)~=0
                    num=num+1;
                    listRxnColour(num,1)={newRxnName};
                    listRxnColour(num,2)={colorStr};
                    end

                elseif p==1
                    colorStr=listColor{ddr};  % ddr is the row number for each reaction;
                end
                fprintf('p value is %d\n',p);
                colorStr=strrep(colorStr,'#','');
                fprintf('colorStr is %s\n',colorStr);
                if strcmp(newRxnName,'re9103');
                    disp('goodl');
                end

                for  ddc=1:cw
                    ref.(newRxnName).color{ddr,ddc}=colorStr;
                    fprintf('set %s ''s colour to %s \n',newRxnName,ref.(newRxnName).color{ddr,ddc});
                end
            end
        end

    end
end


text:
%% EXAMPLE MONTE-CARLO SETUP CONFIGURATION FILE (monteCarlo_example.m) %%%%
% This function is intended to demonstrate how the monte-carlo aspect of
% OMAS may be used to run numerous cycles of a series of different
% algorithms, agents and scenarios.

% Author: James A. Douthwaite 30/04/18

% Notes*
% - The monte-carlo object accepts a cell array of objects [studies,objects]
%   (where 'studies' are seen as independant sets of cycles/analytics).
% - It is currently more efficient to hand smaller sets of object arrays
%   per Monte-Carlo instance rather than run many studies using one.
% - The object array can be defined in same manner as OMAS, using any
%   scenario or object set, the Monte-Carlo instance will then perturb the
%   global positions of all objects before the cycle begins.

% ADD THE PROGRAM PATHS
clear all; close all;
addpath('environment');
addpath('objects');
addpath('scenarios');                                                      % Required system/toolbox paths

fprintf('[SETUP]\tInitialising Monte-Carlo example.\n');
% /////////////////// AGENT & SCENARIO PARAMETERS /////////////////////////
sim_agentOrbit  = 30;
sim_agentVelocity = 0;
sim_waypointOrbit = 35;

% //////////////////////// OMAS CONFIGURATION /////////////////////////////
sim_warningDistance = 10;
sim_maxDuration = 120;
sim_timeStep = 0.25;
sim_verbosity = 0;

% ///////////////////// MONTE-CARLO CONFIGURATION /////////////////////////
monteCarlo_algorithms = {'agent_2D_VO',...
                         'agent_2D_RVO',...
                         'agent_2D_HRVO'};
                         %'agent_2D_RVO2'};                                % The types of agents simulated
monteCarlo_agentPopulations = [2,4];

% ///////////////////// MONTE-CARLO CONFIGURATION /////////////////////////
[~, userdir] = system('echo %USERPROFILE%');    % Get desktop path
monteCarlo_outputPath = strcat(userdir,'\desktop\US18_data');
monteCarlo_positionSigma = 0.2;
monteCarlo_velocitySigma = 0;
monteCarlo_cycles = 5;                        % [FIXED]
monteCarlo_parallel = 0;
monteCarlo_offOnComplete = 0;

% ///////////////////// DEFINE THE MC OBJECT MATRIX ///////////////////////
studyVector = cell(1,2*max(monteCarlo_agentPopulations));
studyObjects = [];
for n = 1:numel(monteCarlo_agentPopulations)                               % Differing agent numbers
    % FOR EACH ALGORITHM TO BE EVALUATED
    for index = 1:numel(monteCarlo_algorithms)                             % Differing agent algorithms
        % CREATE THE AGENT SET
        agentVector = cell(1,monteCarlo_agentPopulations(n));              % Agent container
        for i = 1:monteCarlo_agentPopulations(n)
            agentVector{1,i} = eval(monteCarlo_algorithms{index});         % Assign agents
        end
        % ///////////////////// SCENARIO CONFIGURATION ////////////////////
        % Concentric scenario
%         initialisedObjects = GetScenario_concentric('agents',agentVector,...
%                                                     'agentOrbit',sim_agentOrbit,...
%                                                     'agentVelocity',sim_agentVelocity,...
%                                                     'waypointOrbit',sim_waypointOrbit);
        % Opposing lines scenarios
        initialisedObjects = GetScenario_twoLines('agents',agentVector,...
                                                  'agentVelocity',sim_agentVelocity);
        % /////////////////////////////////////////////////////////////////
        sessionSample = studyVector;
        sessionSample(1:numel(initialisedObjects)) = initialisedObjects;
        studyObjects = vertcat(studyObjects,sessionSample);
    end
end

%% ///////////// BEGIN SIMULATING THE MONTE-CARLO INSTANCES ///////////////
fprintf('[SETUP]\t Defining Monte-Carlo Simulation Series...\n');
[monteCarlo] = OMAS_monteCarlo('objects',studyObjects,...
    'duration',sim_maxDuration,...
    'dt',sim_timeStep,...
    'warningDistance',sim_warningDistance,...
    'verbosity',sim_verbosity,...
    'positionSigma',monteCarlo_positionSigma,...
    'velocitySigma',monteCarlo_velocitySigma,...
    'threadPool',monteCarlo_parallel,...
    'shutDownOnComplete',monteCarlo_offOnComplete,...
    'cycles',monteCarlo_cycles,...
    'outputPath',monteCarlo_outputPath);

% EVALUATE ALL OF THE PROPOSED CONDITIONS
monteCarlo.evaluateAllCycles();

repo: OpenMAS


text:
function [ r, seed ] = r8vec_uniform_01 ( n, seed )

%*****************************************************************************80
%
%% R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real R(N,1), the vector of pseudorandom values.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  i4_huge = 2147483647;

  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_UNIFORM_01 - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'R8VEC_UNIFORM_01 - Fatal error!' );
  end

  r = zeros ( n, 1 );

  for i = 1 : n

    k = floor ( seed / 127773 );

    seed = 16807 * ( seed - k * 127773 ) - k * 2836;

    if ( seed < 0 )
      seed = seed + i4_huge;
    end

    r(i,1) = seed * 4.656612875E-10;

  end

  return
end


text:
function r = mtimes(a,b,takeRight)
% orientation times Miller and orientation times orientation
%
% Syntax
%   o = o1 * o2
%   r = o * h
%   h = inv(o) * r
%
% Input
%  o - @orientation
%  h - @Miller indice
%  r - @vector3d
%
% See also
% orientation/times

% this is some shortcut for internal use
if nargin == 3
  r = mtimes@rotation(a,b,takeRight);
  return
end

% special case multiplication with +-1
if isnumeric(a) || isnumeric(b)
  r = mtimes@rotation(a,b);
  return
end

% orientation times symmetry
if isa(b,'symmetry')
  r = mtimes@quaternion(a,b.rot,0);
  return
end

% ensure inner symmetries coincide
[a, left, right] = ensureSym(a,b);

% orientation times object
if ~isa(b,'quaternion')
  r = rotate_outer(b,a);
  return
end

% rotation multiplication
r = mtimes@quaternion(a,b,isa(b,'orientation'));

% convert back to orientation
if isa(right,'crystalSymmetry') || isa(left,'crystalSymmetry')
  r.CS = right;
  r.SS = left;
else % otherwise it is only a rotation anymore
  r = rotation(r);
end

end

author: mtex-toolbox

repo: mtex

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


text:
function loaded_image = MimLoadImageFromDicomFiles(image_path, filenames, reporting)
    % MimLoadImageFromDicomFiles. Loads a series of DICOM files into a 3D volume
    %
    %     Syntax
    %     ------
    %
    %         loaded_image = MimLoadImageFromDicomFiles(image_path, filenames, reporting)
    %
    %             loaded_image    a PTKImage containing the 3D volume
    %
    %             image_path      location of the DICOM files to load
    %
    %             filenames       a filename string, or cell array of filenames
    %
    %             reporting (optional) - an object implementing CoreReportingInterface
    %                             for reporting progress and warnings
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    % Create a reporting object if none was provided
    if nargin < 3
        reporting = CoreReportingDefault();
    end

    % A single filename canbe specified as a string
    if ischar(filenames)
        filenames = {filenames};
    end

    dicomLibrary = DMFallbackDicomLibrary.getLibrary();

    [image_volume_wrapper, representative_metadata, slice_thickness, global_origin_mm] = DMLoadMainImageFromDicomFiles(image_path, filenames, dicomLibrary, reporting);

    if ~isempty(representative_metadata) && isfield(representative_metadata, 'Modality') && ~isempty(representative_metadata.Modality)
        if ~MimModalityIsSupported(representative_metadata.Modality)
            reporting.Error('MimLoadImageFromDicomFiles:ModalityNotSupported', ['The ' representative_metadata.Modality ' modality is not supported']);
        end
    end

    if isempty(image_volume_wrapper.RawImage)
        reporting.Error('MimLoadImageFromDicomFiles:NoPixelData', 'The DICOM file contains no pixel data');
    end

    if isempty(slice_thickness)
        reporting.ShowWarning('MimLoadImageFromDicomFiles:NoSliceThickness', 'No information found about the slice thickness. Arbitrarily setting slice thickness to 1');
        slice_thickness = 1;
    end

    % Detect and remove padding values
    padding_value = MimRemovePaddingValues(image_volume_wrapper, representative_metadata, reporting);

    % Construct a PTKDicomImage from the loaded image volume
    loaded_image = PTKDicomImage.CreateDicomImageFromMetadata(image_volume_wrapper, representative_metadata, slice_thickness, global_origin_mm, padding_value, reporting);

end

author: tomdoel

repo: pulmonarytoolkit


text:
function compspec
% Competing species phase diagram
%    using MATLAB ode
%
%   $Ekkehard Holzbecher  $Date: 2006/09/04 $
%--------------------------------------------------------------------------
T = 1000;                % maximum time
r = [1; 1];              % rates
e = [1; 1];              % equilibria
lambda = 1;              % lambda parameter
c0 = [1; 1];             % initial concentrations

gtraj = 1;               % trajectory plot
gquiv = 20;              % arrow field plot; value for no. of arrows in 1D
xmin = 0; xmax = 1;      % x- interval for arrow field plot
ymin = 0; ymax = 1;      % y-     "     "    "     "     "
scale = 2;               % scaling factor for arrows
%----------------------execution-------------------------------------------

options = odeset('AbsTol',1e-20);
[~,c] = ode15s(@CS,[0 T],c0,options,r,e,lambda);

%---------------------- graphical output ----------------------------------

if (gtraj)
    plot (c(:,1)',c(:,2)'); hold on;
    plot (e(1),0,'s'); plot (0,e(2),'s');
    legend ('trajectory');
    xlabel ('specie 1'); ylabel ('specie 2');
    title ('Competing species');
end
if (gquiv)
    [x,y] = meshgrid (linspace(xmin,xmax,gquiv),linspace(ymin,ymax,gquiv));
    dy = zeros(gquiv,gquiv,2);
    for i = 1:gquiv
        for j = 1:gquiv
            dy(i,j,:) = CS(0,[x(i,j);y(i,j)],r,e,lambda);
        end
    end
    quiver (x,y,dy(:,:,1),dy(:,:,2),scale);
end

%---------------------- function ------------------------------------------
function dydt = CS(~,y,r,e,lambda)
k = [e(1)/(1+lambda*y(2)/y(1)); e(2)/(1+y(1)/y(2)/lambda)];
dydt = r.*y.*(1-y./k);

author: Sable

repo: mcbench-benchmarks


text:
function I = mi_gd2(x, y, Ym, biascorrect, demeaned, class_means)
% MI_GD Mutual information (MI) between a Gaussian and a discrete
%         variable in bits
%   I = mi_gd(x,y,Ym) returns the MI between the (possibly multidimensional)
%   Gaussian variable x and the discrete variable y.
%   Rows of x correspond to samples, columns to dimensions/variables.
%   (Samples first axis)
%   y should contain integer values in the range [0 Ym-1] (inclusive).
%
%   biascorrect : true / false option (default true) which specifies whether
%   bias correction should be applied to the esimtated MI.
%   demeaned : false / true option (default false) which specifies whether the
%   input data already has zero mean (true if it has been copula-normalized)

persistent previous_nvar previous_y bias_unc bias_cond;

% ensure samples first axis for vectors
if isvector(x)
    x = x(:);
end
if ndims(x)~=2
    error('mi_gd: input arrays should be 2d')
end
if isvector(y)
    y = y(:);
else
%    error('mi_gd: only univariate discrete variable supported');
end

Ntrl = size(x,1);
Nvar = size(x,2);

if isequal(previous_y, y)
  computebias = false;
else
  computebias = true;
end

if ~isequal(previous_nvar, Nvar)
  computebias = true;
end

if size(y,1) ~= Ntrl
    error('mi_gd: number of trials do not match');
end

% default option values
if nargin<4
    biascorrect = true;
end
if nargin<5
    demeaned = false;
end

if ~demeaned
    x = bsxfun(@minus,x,sum(x,1)/Ntrl);
end

% class-conditional entropies
Ntrl_y = sum(y);
Hcond = zeros(1,Ym);
for yi=1:Ym
    xm = x(y(:,yi),:);
    %Ntrl_y(yi) = size(xm,1);
    %xm = bsxfun(@minus,xm,sum(xm,1)/Ntrl_y(yi));
    Cm = (xm'*xm) / (Ntrl_y(yi) - 1);
    chCm = chol(Cm);
    Hcond(yi) = sum(log(diag(chCm)));% + c*Nvar;
end
% class weights
w = Ntrl_y ./ Ntrl;

% input data is class-demeaned, this needs to be accounted for in the
% unconditional entropies
c = diag(sqrt(Ntrl_y))*class_means;%*diag(Ntrl_y);

% unconditional entropy from unconditional Gaussian fit
Cx = (x'*x + c'*c) / (Ntrl-1);
chC = chol(Cx);
Hunc = sum(log(diag(chC)));% + c*Nvar; % the commented out bit drops out in the subtraction below

% apply bias corrections
ln2 = log(2);
if biascorrect

    if computebias,
      vars = 1:Nvar;

      psiterms_unc = psi((Ntrl - vars)/2) / 2;
      dterm_unc    = (ln2 - log(Ntrl-1)) / 2;
      bias_unc     = Nvar*dterm_unc + sum(psiterms_unc);

      dterm_cond    = (ln2 - log(Ntrl_y-1)) / 2;
      psiterms_cond = zeros(1,Ym);
      for vi=vars
        idx = (Ntrl_y-vi);
        psiterms_cond = psiterms_cond + psi(idx/2);
      end
      bias_cond = Nvar*dterm_cond + (psiterms_cond/2);
    end

    Hunc  = Hunc  - bias_unc;
    Hcond = Hcond - bias_cond;
end

I = Hunc - w*Hcond(:);% sum(w .* Hcond);
% convert to bits
I = I / ln2;

previous_y = y;
previous_nvar = Nvar;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function roi = roiGrayDisk(seg, mr, radius, startPt);
% Create an ROI as a disk of gray matter with the specified startPt
% and radius.
%
% roi = roiGrayDisk(seg, mr, [radius=dialog], [startPt=mesh cursor]);
%
%
% ras, 02/2007.
if nargin<2, error('Not enough input args.');	end
if notDefined('mapMethod'),     mapMethod = 'all'; end

msh = segGet(seg, 'SelectedMesh');
[nodes edges] = segGet(seg, 'gray');

if isempty(msh.vertexGrayMap) | all(msh.vertexGrayMap(:)==0)
    msh.vertexGrayMap = mrmMapVerticesToGray(msh.initVertices, nodes, ...
							msh.mmPerVox, edges);
end

if notDefined('startPt')
	vtx = mrmGet(msh, 'cursorvertex');

	if vtx<1
		error('Mesh Cursor not pointing to a vertex on the mesh.')
	end

	I = msh.vertexGrayMap(1,vtx); % layer 1 node for this vertex
	startPt = nodes([2 1 3],I);
end

if notDefined('radius')
	r = inputdlg({'Enter Gray Disk Radius (mm):'}, mfilename, 1, {'3'});
	radius = r{1};
end

%% initialize empty ROI based on the same anatomy coords as the segmentation
roi = roiCreate('I|P|R');
roi.voxelSize = msh.mmPerVox;
if isstruct(mr)
	roi.referenceMR = mr.name;
else
	roi.referenceMR = mr;
end

%% find index of gray node closest to start point
startNode = segGet(seg, 'NearestNode', startPt);

%% compute the distance between each gray node and the start point
dist = mrManDist(nodes, edges, startNode, msh.mmPerVox, -1, radius);

%% find those nodes whose distance is within the radius
inside = find(dist >= 0); % only searches within radius; invalid is -1

%% get the coords of those nodes, set as ROI coords
roi.definedCoords = nodes([2 1 3],inside);
roi.coords = roi.definedCoords;
roi = roiCheckCoords(roi, mr);
roi.name = sprintf('Disk ROI (%i mm about %s)', radius, num2str(startPt'));
roi.comments = [sprintf('Created by %s %s \n', mfilename, datestr(now)) ...
				sprintf('Start Point %s \n', num2str(startPt)) ...
				sprintf('Start Node %i \n', startNode) ...
				sprintf('Radius %i mm', radius)];

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function hist_han = histogram(obj, varargin)
% Create a histogram of image values or a series of histograms for each
% image in an image_vector (e.g., fmri_data) object
%
% :Usage:
% ::
%
%     hist_handles = histogram(obj, [optional inputs])
%
% For objects: Type methods(object_name) for a list of special commands
%              Type help object_name.method_name for help on specific
%              methods.
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2016 Tor Wager
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
% ..
%
% :Inputs:
%
%   **obj:**
%       an image_vector object
%
% :Optional Inputs:
%
% **'byimage', 'separate'**
%       Plot histogram for each image separately.
%       The default is to aggregate across all images in the object.
%
% **'by_tissue_type'**
%       Separate into gray, white, CSF tissue compartments and plot
%       histograms for each compartment.
%       This uses the extract_gray_white_csf method, which in turn
%       currently uses the images
%       'gray_matter_mask.img' 'canonical_white_matter.img' 'canonical_ventricles.img'
%       These images are based on the SPM8 a priori tissue probability
%       maps, but they have been cleaned up and made symmetrical and/or eroded
%       so that the white and CSF compartments are unlikely to contain very
%       much gray matter.  The gray compartment is currently more
%       inclusive. The potential value of this is that signal in the CSF/white
%       compartments may be removed from images prior to/during analysis
%
% % **'singleaxis'**
%       If (only if) 'byimage', plot all histograms on single axis in
%       different random colors. Replaces default behavior of one axis per
%       image.
%
% **'plot'**
%       No function right now.
%
% **'nbins'**
%       Followed by number of bins to use (default is 'auto', n_values/1000)
%
% **'color'**
%       Followed by color [r g b] triplet or text, e.g., 'b'
%
% **'mask'**
%       Followed by mask image name or fmri_data object
%
% **'doline'**
%       Plot lines for histograms
%
% **'noline'**
%       Suppress lines for histograms
%
% **'nofill'**
%       Suppress histogram fill
%
% **'nofigure'**
%       Suppress figure creation
%
% :Outputs:
%
%   **hist_han:**
%        Graphics handles to histogram lines or bars
%
%   **out2:**
%        description of out2
%
% :Examples:
% ::
%
% hist_han = histogram(wp_alone, 'byimage', 'color', 'b');
% hist_han = histogram(wp_alone, 'byimage', 'by_tissue_type');
% hist_han = histogram(wp_alone, 'by_tissue_type');
%
%
% :References:
%   None. These are just histograms.
%
% :See also:
%   - image_histogram, image_intensity_histograms
%
%
% ..
%    Programmers' notes:
%   Documentation and update, July 2016, Tor Wager
%   Added line to maximize figure window, May 2022, Lukas Van Oudenhove
% ..

% -------------------------------------------------------------------------
% DEFAULTS AND INPUTS
% -------------------------------------------------------------------------

doplot = 1;
doline = 1;
dofill = 1;
dofigure = 1;
do_by_image = 0;
nbins = [];
color = [.3 .3 .3];
mask = [];
by_tissue_type = 0;
singleaxis = false;

% initalize optional variables to default values here.

% optional inputs with default values
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case {'byimage', 'by_image', 'separate'}, do_by_image = 1;

            case 'singleaxis', singleaxis = true;

            case 'plot', doplot = 1;

            case 'nbins', nbins = varargin{i+1}; varargin{i+1} = [];
            case 'color', color = varargin{i+1}; varargin{i+1} = [];
            case 'mask', mask = varargin{i+1}; varargin{i+1} = [];

            case 'doline', doline = 1;
            case 'noline', doline = 0;
            case 'nofill', dofill = 0;
            case {'nofigure', 'nofig'}, dofigure = 0;

            case 'by_tissue_type', by_tissue_type = 1; varargin{i} = [];

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

% -------------------------------------------------------------------------
% SPECIAL FUNCTION: BY TISSUE TYPE
% -------------------------------------------------------------------------
if by_tissue_type

    [values, components, full_data_objects] = extract_gray_white_csf(obj);

    hist_han = cell(1, 3);
    colors = {[.2 .2 .2] [1 .6 .6] [.3 .3 1]};

    for i = 3:-1:1

        hist_han{i} = histogram(full_data_objects{i}, varargin{:}, 'color', colors{i}, 'nofigure', 'nbins', 100, 'nofill');

        full_data_objects{i} = remove_empty(full_data_objects{i});  % just in case

        stdevs(:, i) = nanstd(full_data_objects{i}.dat);

    end

    % Plot locations of each and relationships

    create_figure('relationships', 1, 3);
    set(gcf,'WindowState','maximized'); %@lukasvo76 added to improve figure quality when publishing
    disp('------------------------------------------------------------')
    disp('Plotting mean values by tissue type')
    fprintf('Gray = gray, red = white, blue = CSF\n');

    barplot_columns(values, 'nofig', 'colors', colors, 'names', {'gray' 'white' 'CSF'});
    %set(gca, 'XTickLabel', {'gray' 'white' 'CSF'});
    title('Mean values by tissue type');

    subplot(1, 3, 2);
    disp('------------------------------------------------------------')
    disp('Plotting spatial standard deviations across voxels by tissue type')
    fprintf('Gray = gray, red = white, blue = CSF\n');

    barplot_columns(stdevs, 'nofig', 'colors', colors, 'names', {'gray' 'white' 'CSF'});
    %set(gca, 'XTickLabel', {'gray' 'white' 'CSF'});
    title('Std. Dev. by tissue type');

    subplot(1, 3, 3);
    plotmatrix([values stdevs]);
    title('Relations [means stdevs]')

    return

end

% -------------------------------------------------------------------------
% MASKING
% -------------------------------------------------------------------------

if ~isempty(mask)

    if isa(mask, 'image_vector')
        obj = apply_mask(obj, mask);
    else
        obj = apply_mask(obj, fmri_data(mask));
    end

end

% -------------------------------------------------------------------------
% HISTOGRAMS
% -------------------------------------------------------------------------

if do_by_image
    % By image: one figure per image in object
    % -----------------------------------------------------

    % Set scale range
    Xtmp = obj.dat(:);
    Xtmp(Xtmp == 0 | isnan(Xtmp)) = [];
    mylim = prctile(Xtmp, [.1 99.9]);

    % Set n bins if auto-select
    if isempty(nbins), nbins = ceil(size(obj.dat, 1) ./ 1000); end

    % Set up subplots
    nimgs = size(obj.dat, 2);
    [rows, cols] = deal(floor(nimgs .^ .5));
    while rows * cols < nimgs, cols = cols + 1; end
    naxes = rows * cols;

    if dofigure
        if singleaxis
            create_figure('histogram');
            set(gcf,'WindowState','maximized'); %@lukasvo76 added to improve figure quality when publishing
        else
            create_figure('histogram', rows, cols);
            set(gcf,'WindowState','maximized'); %@lukasvo76 added to improve figure quality when publishing
        end
    end

    if singleaxis
        % Set systematic colors according to index number
        if nimgs < 24
            rgbcell = seaborn_colors(nimgs);
        else
            rgbcell = seaborn_colors(24);

            while length(rgbcell) < nimgs
                rgbcell = [rgbcell; seaborn_colors(min(24, nimgs - length(rgbcell)))];
            end
        end

    end

    % Do for each image
    for i = 1:nimgs

        if singleaxis

            hist_han(i) = create_hist(obj.dat(:, i), nbins, doline, false, rgbcell{i});

        else

            subplot(rows, cols, i);
            hist_han(i) = create_hist(obj.dat(:, i), nbins, doline, dofill, color);

        end

        % set axes
        set(gca, 'XLim', mylim);

    end

    if ~singleaxis
        % turn off remaining axes
        for i = nimgs + 1:naxes
            subplot(rows, cols, i)
            axis off
        end
    else
        % single axis legend
        if nimgs < 40
            case_id = cellfun(@num2str, num2cell(1:nimgs), 'UniformOutput', false);
            legend(hist_han, case_id)
        end

        xlabel('Value'); ylabel('Frequency')
        title('Histogram of values');

    end

else
    % Overall: one histogram
    % -----------------------------------------------------

    % Set n bins if auto-select
    if isempty(nbins), nbins = ceil(size(obj.dat, 1) ./ 500); end

    if dofigure
        create_figure('histogram');
        set(gcf,'WindowState','maximized'); %@lukasvo76 added to improve figure quality when publishing
    end

    Xtmp = obj.dat(:);
    hist_han = create_hist(Xtmp, nbins, doline, dofill, color);

    % Line at 0
    hh = plot_vertical_line(0);
    set(hh, 'LineStyle', '--');

    xlabel('Value'); ylabel('Frequency');
    title('Histogram of values');

end

end % FUNCTION

function han = create_hist(Xtmp, nbins, doline, dofill, color)

% remove zeros
Xtmp(Xtmp == 0 | isnan(Xtmp)) = [];

% Get bins
x = linspace(prctile(Xtmp, .01), prctile(Xtmp, 99.9), nbins);

[h, x] = hist(Xtmp, x);

% convert to PDF
h = h ./ sum(h);

if doline
    han = plot(x, h, 'LineWidth', 3, 'color', color);

else
    han = bar(x, h);
    set(han, 'FaceColor', color, 'EdgeColor', 'none');

end

if dofill
    fill([x(1)-eps x x(end)+eps], [0 h 0], color, 'FaceAlpha', .2);
end

hold on; % still trouble with automatic turn-off

axis tight

% set x-axis
% should not need to do this if we've chosen bins carefully above
% wh = find(h > .0001);
% xlims = [x(min(wh)) x(max(wh))];
% set(gca, 'XLim', xlims);

drawnow

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function X = my_fftn(x)

X = fft(x);
dd = ndims(x);
for id=2:dd
	ord = 1:dd;
	ord([1 id]) = [id 1];
	x = permute(X, ord);
	X = fft(x);
	X = ipermute(X, ord);
end

author: thomaskuestner


text:
function [V, policy, iter, cpu_time] = mdp_policy_iteration(P, R, discount, policy0, max_iter, eval_type)

% mdp_policy_iteration  Resolution of discounted MDP
%                       with policy iteration algorithm
% Arguments ---------------------------------------------------------------
% Let S = number of states, A = number of actions
%   P(SxSxA) = transition matrix
%              P could be an array with 3 dimensions or
%              a cell array (1xA), each cell containing a matrix (SxS) possibly sparse
%   R(SxSxA) or (SxA) = reward matrix
%              R could be an array with 3 dimensions (SxSxA) or
%              a cell array (1xA), each cell containing a sparse matrix (SxS) or
%              a 2D array(SxA) possibly sparse
%   discount = discount rate, in ]0, 1[
%   policy0(S) = starting policy, optional
%   max_iter = maximum number of iteration to be done, upper than 0,
%              optional (default 1000)
%   eval_type = type of function used to evaluate policy:
%              0 for mdp_eval_policy_matrix, else mdp_eval_policy_iterative
%              optional (default 0)
% Evaluation --------------------------------------------------------------
%   V(S)   = value function
%   policy(S) = optimal policy
%   iter     = number of done iterations
%   cpu_time = used CPU time
%--------------------------------------------------------------------------
% In verbose mode, at each iteration, displays the number
% of differents actions between policy n-1 and n

% MDPtoolbox: Markov Decision Processes Toolbox
% Copyright (C) 2009  INRA
% Redistribution and use in source and binary forms, with or without modification,
% are permitted provided that the following conditions are met:
%    * Redistributions of source code must retain the above copyright notice,
%      this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright notice,
%      this list of conditions and the following disclaimer in the documentation
%      and/or other materials provided with the distribution.
%    * Neither the name of the <ORGANIZATION> nor the names of its contributors
%      may be used to endorse or promote products derived from this software
%      without specific prior written permission.
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
% INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
% OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
% OF THE POSSIBILITY OF SUCH DAMAGE.

cpu_time = cputime;

global mdp_VERBOSE;

% check of arguments
if iscell(P); S = size(P{1},1); else S = size(P,1); end;
if discount <= 0 || discount >= 1
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: Discount rate must be in ]0; 1[')
    disp('--------------------------------------------------------')
elseif nargin > 3 && (size(policy0,1)~=S || any(mod(policy0,1)) || any(policy0<1) || any(policy0>S))
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: policy0 must a (1xS) vector with integer from 1 to S')
    disp('--------------------------------------------------------')
elseif nargin > 4 && max_iter <= 0
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: The maximum number of iteration must be upper than 0')
    disp('--------------------------------------------------------')
else

    PR = mdp_computePR(P,R);

    % initialization of optional arguments
    if nargin < 6; eval_type = 0; end;
    if nargin < 5; max_iter = 1000; end;
    if nargin < 4;
        % initialization of policy:
        % the one wich maximizes the expected immediate reward
        [nil, policy0] = mdp_bellman_operator(P,PR,discount,zeros(S,1));
    end;

    if mdp_VERBOSE; disp('  Iteration  Number_of_different_actions'); end;

    iter = 0;
    policy = policy0;
    is_done = false;
    while ~is_done
        iter = iter + 1;
        if  (eval_type==0)
            V = mdp_eval_policy_matrix(P,PR,discount,policy);
        else
            V = mdp_eval_policy_iterative(P,PR,discount,policy);
        end;
        [nil, policy_next] = mdp_bellman_operator(P,PR,discount,V);

	n_different = sum(policy_next ~= policy);
        if mdp_VERBOSE; disp(['       ' num2str(iter) '                 '  num2str(n_different)]); end;

        if all(policy_next==policy) || iter == max_iter
            is_done = true;
        else
            policy = policy_next;
        end;
    end;

end;

cpu_time = cputime - cpu_time;


text:
%% Demonstration of the 2D Discrete Gram Transformation
%
% (C) 2013 Matthew harker and Paul O'Leary
% Institute for Automation
% University of Leoben
% A-8700 Leoben
% Austria
%
% email: office@harkeroleary.org
%
% This file produces a demonstration of the 2D discrete gram polynomial
% transformation of an image and the corresponding inverse transformation.
% It used the function |dgt| to compute the discrete Gram transformation
% and the function |idgt| to compute the inverse transformation.
%
% It is also demonstarted that the error in the reconstrcuted image ins
% negligible, i.e., approaching the numerical accuracy available in MATLAB.
% This demonstrates that the polynomial basis functions are free from error
% evant at very high degrees in this example d = 1024.
%
% This tool requires the discrete orthogonal toolbox: dopBox:
%
%   http://www.mathworks.com/matlabcentral/fileexchange/41250
%
%%
close all;
clear all;
%
% Set some defaults
%
FontSize = 12;
set(0,'DefaultaxesFontName','Times');
set(0,'DefaultaxesFontSize',FontSize);
set(0,'DefaulttextFontName','Times');
set(0,'DefaulttextFontSize',FontSize);
set(0,'DefaultTextInterpreter', 'latex');
%
%% Load an image and convert to a double gray data set
%
% This is an image of an art object. The snails are modelled in modelling
% material, they are not real animals.
%
[pic, map] = imread('snails1.jpg');
%
% Make the data mean free so that the spectrum is better viwible
%
D = double( rgb2gray( pic ));
D = D - mean(D(:));
%% View the Test Image
%
% View the image
%
fig1 = figure;
imagesc( pic );
axis image;
%
title('Photo Finish ($$\copyright$$ 2012 Paul O''Leary and Wolfgang Trettnak)');
xlabel('These are artificial snails and not real animals');
%
%% Compute the Gram Spectrum
%
% Call the |dgt2| function to compute the two dimensional Gram polynomial
% transformation. This is akin to the fft2 for the Fourier basis.
%
[S, Bx, By] = dgt2( D );
%
[ny,nx] = size( S );
%
% Compute the degrees in x and y directions.
%
xScale = 0:(nx-1);
yScale = 0:(ny-1);
%% View the 2D Polynomial Spectrum
%
% This is the complete 2D spectrum for the complete image.
%
fig2 = figure;
imagesc( xScale, yScale, S );
axis image;
colorbar
title('2D Gram Spectrum (information is at lower degrees)');
xlabel('Polynomial degree in $$x$$');
ylabel('Polynomial degree in $$y$$');
%
%% Zoomed Section of Spectrum
%
% Virtually all of the Information of the test image is concentrated at
% lower degrees. This can be shown by zooming in on the relevant portion
% of the spectrum.
%
fig3 = figure;
imagesc( xScale, yScale, S);
axis image;
colorbar
title('2D Gram Spectrum');
xlabel('Polynomial degree in $$x$$');
ylabel('Polynomial degree in $$y$$');
axis([0,20,0,20]);
%
%
%% Reconstruct the Image and Test for Errors
%
% The aim here is to reconstruct the image from its spectrum and then to
% compute a norm for the total error. This gives a measure for the quality
% of the transform inverse-transform pair.
%
% Note: the basis functions computed for the transform are passed to the
% inverse transform. this avoids the necessity to recompute the polynomial
% basis functions.
%
% As can be seen in the final figure, the Frobeneus norm of the error is in
% the order of $10^{-15}$, consequently, the polynomial transformation can
% be regarded as free from error.
%
Dr = idgt2( S, Bx, By );
%
% Calculate the difference between the original gray scale image and its
% reconstruction.
%
E = D - Dr;
%
% Compute the ratio of the Frobenius of the error to the norm of the data.
%
normE = norm( E, 'fro');
normD = norm( D, 'fro');
normEtoD = normE / normD;
%
%% View the Reconstructed Image.
%
fig3 = figure;
imagesc( Dr );
axis image;
colormap gray;
axis off;
title(['Reconstructed Gray Image ($$ \epsilon = ',num2str(normEtoD),' $$)']);
xlabel('These are artificial snails and not real animals');


text:
close all
clear all

N=5;
names={'A','B','C','D','E'};
dag=zeros(N);
dag(1,2)=1;
dag(2,[3 4])=1;
dag([3 4],5)=1;
[xx yy] = draw_graph(dag,names,ones(N,1));
title('Cheng example.');

node_sizes=2*ones(1,N);
bnet=mk_bnet(dag,node_sizes);

if 1
  disp('Generating DataSet');
  bnet.CPD{1} = tabular_CPD(bnet, 1, [0.4 0.6]);
  bnet.CPD{2} = tabular_CPD(bnet, 2, [0.2 0.3 0.8 0.7]);
  bnet.CPD{3} = tabular_CPD(bnet, 3, [0.1 0.2 0.9 0.8]);
  bnet.CPD{4} = tabular_CPD(bnet, 4, [0.6 0.8 0.4 0.2]);
  bnet.CPD{5} = tabular_CPD(bnet, 5, [0.9 0.8 0.7 0.6 0.1 0.2 0.3 0.4]);
  m=10000;
  cheng = cell(N,m);
  for i=1:m
    cheng(:,i)=sample_bnet(bnet);
  end
  cheng = cell2num(cheng);
  %save -ascii cheng cheng
else
  load -ascii cheng.mat
end

%  profile clear
%  profile on
  [Phase_3, Phase_2, Phase_1, UPhase_3] = learn_struct_bnpc(cheng, node_sizes, 0.05, 0)
%  profile off
%  profile report report_cheng

figure
draw_graph(Phase_1,names,ones(N,1),xx,yy);
title('PhaseI');
figure
draw_graph(Phase_2,names,ones(N,1),xx,yy);
title('PhaseII');
%figure
%draw_graph(UPhase_3,names,ones(N,1),xx,yy);
%title('undirected PhaseIII');
figure
draw_graph(Phase_3,names,ones(N,1),xx,yy);
title('PhaseIII');

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [ o, x, w ] = en_r2_05_2 ( n )

%*****************************************************************************80
%
%% EN_R2_05_2 implements the Stroud rule 5.2 for region EN_R2.
%
%  Discussion:
%
%    The rule has order O = 2 * N^2 + 1.
%
%    The rule has precision P = 5.
%
%    EN_R2 is the entire N-dimensional space with weight function
%
%      w(x) = exp ( - x1^2 - x2^2 ... - xn^2 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971,
%    ISBN: 0130438936,
%    LC: QA311.S85.
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Output, integer O, the order.
%
%    Output, real X(N,O), the abscissas.
%
%    Output, real W(O), the weights.
%
  o = 2 * n * n + 1;
  volume = sqrt ( pi^n );

  a = 2 * volume / ( n + 2 );
  b = ( 4 - n ) * volume / 2 / ( n + 2 )^2;
  c = volume / ( n + 2 )^2;

  r = sqrt ( ( n + 2 ) / 2 );
  s = sqrt ( ( n + 2 ) / 4 );

  x = zeros ( n, o );
  w = zeros ( o, 1 );

  k = 0;
%
%  1 point.
%
  k = k + 1;
% x(1:n,k) = 0;
  w(k) = a;
%
%  2 * N points.
%
  for i = 1 : n
    k = k + 1;
    x(i,k) = - r;
    w(k) = b;
    k = k + 1;
    x(i,k) = + r;
    w(k) = b;
  end
%
%  4 * ( N * ( N - 1 ) / 2 ) points.
%
  for i = 1 : n - 1
    for j = i + 1 : n
      k = k + 1;
      x(i,k) = - s;
      x(j,k) = - s;
      w(k) = c;
      k = k + 1;
      x(i,k) = - s;
      x(j,k) = + s;
      w(k) = c;
      k = k + 1;
      x(i,k) = + s;
      x(j,k) = - s;
      w(k) = c;
      k = k + 1;
      x(i,k) = + s;
      x(j,k) = + s;
      w(k) = c;
    end
  end

  return
end


text:
Ie=imread('background05_empty.jpg');
%imshow(Ie);

I=imread('background05.jpg');
% imshow(I);
dI=double(rgb2gray(Ie))-double(rgb2gray(I));

%imshow(dI,[min(dI(:))  max(dI(:))]);

tr=0.45*max(max(abs(dI(:))));
% tr=0.1*max(abs(dI(:)));

bw=abs(dI)>tr;

%imshow(bw);

bw(226,426:428)=false;
bw(227,328)=false;
% imshow(bw);

L = bwlabel(bw,4);
map=rand(max(L(:))+1,3);
%imshow(L+1,map);

ns=max(L(:)); % number of initial segments

% find rectangles of initial segments:
six=zeros(ns,1);
siy=zeros(ns,1); % positions
six2=zeros(ns,1);
siy2=zeros(ns,1); % positions
for c=1:ns
    bwt=(L==c);
    [yt xt]=find(bwt);
    x1=min(xt);
    x2=max(xt);
    y1=min(yt);
    y2=max(yt);
    six(c)=x1;
    siy(c)=y1;
    six2(c)=x2;
    siy2(c)=y2;
end

[sy sx tmp3]=size(Ie);

%in_rect=false(sy,sx,ns);
in_rect=false(sy*sx,ns);
waitbar(0,'in rect');
for yc=1:sy
    for xc=1:sx
        for c=1:ns
            if (siy(c)<=yc)&&(yc<=siy2(c))&&(six(c)<=xc)&&(xc<=six2(c))
                %in_rect(yc,xc,c)=true;
                in_rect(yc+(xc-1)*sy,c)=true;
            end
        end
    end
    waitbar(yc/sy);
end

in_rect_u=unique(in_rect,'rows'); % unicue pixels class
% in_rect_u(class_counter,:) - what segments used

nsr0=size(in_rect_u,1); % new number of segments
nsr=nsr0-1;

sixr=zeros(nsr,1);
siyr=zeros(nsr,1); % positions
six2r=zeros(nsr,1);
siy2r=zeros(nsr,1); % positions

% find limits:
%in_rect(pixel_number,segment_number)
%for cr=1:nsr
tic
waitbar(0,'new limits');
for cr1=1:nsr % exlude big
    cr=cr1+1;
    in_rect_u_t=in_rect_u(cr,:); % in_rect_u_t(segment_number)
    sgu=find(in_rect_u_t); % what segments used

    rect=false(sy,sx);

    for yc=1:sy
        for xc=1:sx
            %if all(in_rect_u_t==())
            pn=yc+(xc-1)*sy; % pixel number
            %in_rect(pn,:) - in wat seggments curent pixel
            if all(in_rect_u_t==in_rect(pn,:))
                rect(yc,xc)=true;
            end
        end
    end

    [py px]=find(rect);
    siyr(cr1)=min(py);
    siy2r(cr1)=max(py);
    sixr(cr1)=min(px);
    six2r(cr1)=max(px);
    waitbar(cr1/nsr);
end

sir=cell(nsr,1); % small images
sibwr=cell(nsr,1); % small images, masks
for cr1=1:nsr % exlude big
    cr=cr1+1;
    y1=siyr(cr1);
    y2=siy2r(cr1);
    x1=sixr(cr1);
    x2=six2r(cr1);
    rect=false(sy,sx);
    rect(y1:y2,x1:x2)=true;
    bwt=false(sy,sx);
    bwt(rect)=bw(rect);
    It=trim_with_mask(I,bwt);
    sir{cr1}=It(y1:y2,x1:x2,:);
    sibwr{cr1}=bwt(y1:y2,x1:x2);
end

% rename back:
six=sixr;
siy=siyr;
six2=six2r;
siy2=siy2r;

si=sir;
sibw=sibwr;

save('segments','si','sibw','six','siy','six2','siy2');


text:
% generates a rounded rectangle as a test track for tracking control applications

l1 = 100;
l2 = 20;
lstep = 2;
r = 12.5;
phistep = 0.1;
phi = phistep:phistep:(pi/2);
b0 = 10;

% build track snippets
trackparts{1} = [lstep:lstep:l1; zeros(1, l1/lstep)];
trackparts{2} = [lstep:lstep:l2; zeros(1, l2/lstep)];
trackparts{3} = [r*cos(phi+3/2*pi); r*sin(phi+3/2*pi)+r];

% order of track snippets
trackgenerator = [1, 3, 2, 3, 1, 3, 2, 3];
% rotation of track snippets
rotgenerator = [0, 0, pi/2, pi/2, pi, pi, 3*pi/2, 3*pi/2] + pi/2;

% build up track
x_m = 0;
y_m = 0;
for i = 1:1:length(trackgenerator)
  R = [cos(rotgenerator(i)), -sin(rotgenerator(i));...
    sin(rotgenerator(i)), cos(rotgenerator(i))];
  % rotate track part to the right position
  x_loc = R(1, :)*trackparts{trackgenerator(i)};
  y_loc = R(2, :)*trackparts{trackgenerator(i)};
  % transform trackpart origin to the last point
  x_m = [x_m, x_loc+x_m(end)];
  y_m = [y_m, y_loc+y_m(end)];
end
% generate width vector
b_m = ones(1, length(x_m))*b0;

% write to csv
dlmwrite('roundedRectangle.csv', [x_m(1:end-1)',y_m(1:end-1)',b_m(1:end-1)'], ';');

figure;
plot(x_m(1:end-1), y_m(1:end-1));
axis equal;
grid on;

author: TUMFTM


text:
function out = ne( a,b )
%ne (~=) for nested cell arrays
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

% determine input types
inType = {class(a), class(b)};
if(strcmp(inType{1},'TRAFO') && strcmp(inType{2},'TRAFO'))
    [outA, idx] = flattenCellMatrix(a.data);
    [outB, idxB] = flattenCellMatrix(b.data);
    meta = a.meta;
    meta_b = b.meta;
    clear 'a' 'b'
    out = false;

    if(~all(cellfun(@(x,y) isequal(x,y), idx, idxB)))
        error('TRAFO::ne: Nested cells must have the same size');
    end
    if (iscell(meta) && iscell(meta_b))
        if(~isempty(meta) && ~isempty(meta_b)) % unequal reconstruction information can happen -> this is what we are interested in to check
            if(iscell(meta{1}) && iscell(meta_b{1}))
                if(~all(cellfun(@(x,y) isequal(x,y), flattenCellMatrix(meta), flattenCellMatrix(meta_b))))
                    out = true;
    %                 error('TRAFO::ne: Unequal reconstruction information');
                end
            else
                if(~all(cellfun(@(x,y) isequal(x,y), meta, meta_b)))
                    out = true;
    %                 error('TRAFO::ne: Unequal reconstruction information');
                end
            end
        end
        out = all(cellfun(@(x,y) ~isequal(x,y), outA, outB)) | out;
    else
        if(~isequal(meta, meta_b))
            out = true; % can not be the same, because the size of the original image stored as meta info is different
        else
            out = all(cellfun(@(x,y) ~isequal(x,y), outA, outB)) | out; %check if entries are the same
        end
    end

elseif(strcmp(inType{1},'TRAFO')) % assume b is double/uint/int
    [out, ~] = flattenCellMatrix(a.data);
%     meta = a.meta;
    clear 'a'

    out = all(cellfun(@(x) ~isequal(x,b), out));

elseif(strcmp(inType{2},'TRAFO'))
    [out, ~] = flattenCellMatrix(b.data);
%     meta = b.meta;
    clear 'b'

    out = all(cellfun(@(x) ~isequal(a,x), out));

else
    error('TRAFO::ne: Impossible constellation!');
end

end


text:
function texim = getTextonImage(im, textonNodes)

if size(im, 3)==3
    im = rgb2gray(im);
end
feat = single(MRS4fast(imfilter(im, fspecial('gaussian', 3, 1))));
%idx = getNearestHierarchy(feat, textonNodes);
% leafnum = zeros(numel(textonNodes), 1);
% leafnum([textonNodes.isleaf]) = 1:sum([textonNodes.isleaf]);
% idx = leafnum(idx);

idx = getNearest(feat, textonNodes.centers);

texim = reshape(idx, size(im));

author: Cloud-CV

repo: object-proposals


text:
function [fmin,xmin]=ConjugateGradientMethod(x0)

% initialization
xk=x0;
gk=grad_obj(xk);
dk=-gk;

% iteration
for i=1:length(x0)
    % line search
    alphak=fminbnd(@(alpha) phi(alpha,xk,dk),0,10);
    % update xk
    tempgk=gk;
    tempxk=xk;
    tempdk=dk;
    xk=tempxk+alphak*tempdk;
    gk=grad_obj(xk);
    tempbetak=(gk'*(gk-tempgk))/(dk'*(gk-tempgk));
    dk=-gk+tempbetak*tempdk;
end
xmin=xk;
fmin=obj(xk);

repo: Optimization-Theory-and-Algorithm


text:
function axesSDG = InitAxesSDG(handles)

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

if isempty(handles)
    hAxes = axes();
else
    hAxes = handles.axesSDG;
end

default_colors = [...
    0.2  0.6  0.1;
    1.0  0.5  0.0;
    1.0  0.0  1.0;
    0.5  0.5  1.0;
    0.0  1.0  1.0;
    1.0  0.0  0.0;
    0.2  0.3  0.1;
    0.8  0.6  0.6;
    0.5  1.0  0.5;
    0.5  1.0  1.0;
    0.0  0.0  0.0;
    0.2  0.2  0.2;
    0.4  0.4  0.4;
    0.6  0.6  0.6;
    0.8  0.8  0.8 ...
];

p = [getAppDir(), 'SDGcolors.csv'];
try
    rgb01 = csvread(p) / 255;
    if size(rgb01, 1) < 6
        rgb01 = default_colors;
        fprintf('SDGcolors.csv is not an appropriate table of more than 6 RGB values. Using default colors.\n')
    else
        fprintf('Loaded %s colors from SDGcolors.csv.\n', num2str(size(rgb01, 1)))
    end
catch  % File not found
    rgb01 = default_colors;
    fprintf('Failed to load SDGcolors.csv. Using default colors.\n')
end

axesSDG = struct(...
                 'handles', struct(...
                                   'axes',hAxes, ...
                                   'SD', [], ...
                                   'ch', [] ...
                                  ), ...
                 'iCh', [], ...
                 'iSrcDet', [], ...
                 'SDPairColors', rgb01...
               );

author: BUNPC


text:
%% The Class vector3d
%
%%
%


text:
function p00_nvar_test ( problem_num )

%*****************************************************************************80
%
%% P00_NVAR_TEST prints the problem size.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM_NUM, the number of problems.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'P00_NVAR_TEST\n' );
  fprintf ( 1, '  List the problem size.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   Problem    Option      Size\n' );
  fprintf ( 1, '\n' );

  for problem = 1 : problem_num
    option_num = p00_option_num ( problem );
    fprintf ( 1, '\n' );
    for option = 1 : option_num
      nvar = p00_nvar ( problem, option );
      fprintf ( 1, '  %8d  %8d  %8d\n', problem, option, nvar );
    end
  end

  return
end


text:
function Vi = spm_get_vc(I,factor)
% Generate error covariance components for factorial designs
% FORMAT Vi = spm_get_vc(I,factor)
% I         - n x m matrix of factor level indicators
%             I(n,i) is the level of factor i for observation n
% factor(i) - structure array of sphericity assumptions for each factor
% .variance - 1 for different variance among levels of factor i
% .dept     - 1 for dependencies within levels of factor i
%
% Vi        - cell vector of covariance components
%__________________________________________________________________________
%
% spm_get_vc generates variance components for a given design. For each
% factor, the user specifies whether its levels have identical variances
% and are independent. The individual components for each factor are
% combined into covariance components by using the Kronecker tensor
% product. If there are unequal number of observations at different levels,
% the function specifies covariance components for a full factorial design
% first and subsequently removes unwanted rows and columns from the
% covariance matrices.
%
% The functionality of spm_get_vc is similar to that of spm_non_sphericity.
% The difference is that spm_get_vc can accommodate any number of factors
% and is more general, because it can cope with different number of
% observations under different levels of a factor.
%__________________________________________________________________________
% Copyright (C) 2006 Freiburg Brain Imaging
% Copyright (C) 2008-2013 Wellcome Trust Centre for Neuroimaging

% Volkmar Glauche
% $Id: spm_get_vc.m 5293 2013-03-01 16:41:46Z guillaume $

%-Numbers of scans and factors
%--------------------------------------------------------------------------
[nscan,nfactor] = size(I);

%-Make sure each row of Iin is unique
%==========================================================================
[Iu,Ii,Ij]  = unique(I,'rows');
if size(Iu,1) < nscan
    nfactor = nfactor + 1;
    uf      = zeros(nscan, 1);
    for k = 1:max(Ij)
        uf(Ij==k) = 1:sum(Ij==k);
    end
    I       = [I uf];
end
nlevel      = max(I);

%-Non-sphericity assumptions
%--------------------------------------------------------------------------
% First factor is replications, assume identical variance and independence.
% Pad with zeroes in case there are less than nfactor factors specified.
variance    = [0 cat(2, factor.variance) zeros(1,nfactor)];
dept        = [0 cat(2, factor.dept) zeros(1,nfactor)];

% (i) generate generic index
%==========================================================================
Igen = zeros(prod(nlevel), nfactor);
Igen(:,1) = kron(ones(1,prod(nlevel(2:end))),1:nlevel(1))';
for cf = 2:(nfactor-1)
    Igen(:,cf) = kron(ones(1,prod(nlevel((cf+1):end))),kron(1:nlevel(cf),ones(1,prod(nlevel(1:(cf-1))))))';
end
Igen(:,nfactor) = kron(1:nlevel(nfactor),ones(1,prod(nlevel(1:(nfactor-1)))))';

% (ii) generate error variance components
%==========================================================================
Vi = {};
for f=1:nfactor

    % identical/non-identical variances
    % for each factor, create a single variance component if variances are
    % identical across levels, and level specific variance components if
    % variances are non-identical
    %----------------------------------------------------------------------
    nVi = {};
    if ~variance(f)
        nVi{1} = speye(nlevel(f),nlevel(f));
    else
        for l1=1:nlevel(f)
            nVi{l1} = sparse(l1,l1,1,nlevel(f),nlevel(f));
        end
    end
    if dept(f)
        for l1 = 1:nlevel(f)
            for l2 = 1:(l1-1)
                nVi{end+1} = sparse([l1 l2],[l2 l1],1,nlevel(f), ...
                                         nlevel(f));
            end
        end
    end

    % combine current factor components with previous ones, thus building
    % up covariance components block by block
    %----------------------------------------------------------------------
    if isempty(Vi)
        Vi = nVi;
    else
        oVi = Vi;
        Vi = {};
        for nv = 1:numel(nVi)
            for ov = 1:numel(oVi)
                Vi{end+1} = kron(nVi{nv}, oVi{ov});
            end
        end
    end
end

% (iii) sort out rows/columns & remove all-zero variance components
%==========================================================================
[unused,ind] = ismember(I,Igen,'rows');
az = false(size(Vi));

for cVi = 1:numel(Vi)
    Vi{cVi} = Vi{cVi}(ind,ind);
    az(cVi) = full(all(Vi{cVi}(:) == 0));
end
Vi = Vi(~az);

dupl = false(size(Vi));
for cVi = 1:numel(Vi)
    if ~dupl(cVi)
        for cVi1 = (cVi+1):numel(Vi)
            dupl(cVi1) = dupl(cVi1)||full(all(Vi{cVi}(:) == Vi{cVi1}(:)));
        end
    end
end
Vi = Vi(~dupl);

author: spm


text:
function huffman_demo2
%HUFFMAN_DEMO2
%   This demo shows compression-decompression of different data.

%   $Author: Giuseppe Ridino' $
%   $Revision: 1.0 $  $Date: 02-Jul-2004 16:47:25 $

fprintf('\n')
disp('This is an example of function using Huffman compression-decompression')
disp('Press a key to continue...')
pause
fprintf('\n')

fprintf('\n')
disp('The string below will be compressed and decompressed')
disp('using norm2huff and huff2norm functions')
disp('Press a key to continue...')
pause
fprintf('\n')

data1 = 'This is a string demo. It will be compressed using norm2huff function.'

fprintf('\n')
disp('Now we will compress and decompress it.')
disp('But pay attention that data mut be converted to uint8!')
disp('Press a key to continue...')
pause
fprintf('\n')

% string data
data1 = uint8(data1);
[zipped1,info1] = norm2huff(data1)
unzipped1 = char(huff2norm(zipped1,info1))

isOK = isequal(data1,unzipped1)

whos data1 zipped1 unzipped1

bar(frequency(data1)); axis tight

%  random data1
fprintf('\n')
disp('Now we will try to compress a random sequence of 1000 bytes')
disp('Press a key to continue...')
pause
fprintf('\n')

data3 = uint8(256*rand(1,1000));
[zipped3,info3] = norm2huff(data3);
unzipped3 = huff2norm(zipped3,info3);

isOK = isequal(data3,unzipped3)

whos data3 zipped3 unzipped3

bar(frequency(data3)); axis tight

fprintf('\n')
disp('It is not well compressed. Why?')
disp('This is because the sequence is normally random!')
disp('Press a key to continue...')
pause
fprintf('\n')

%  random data2
fprintf('\n')
disp('Now we will try to compress a random sequence of 1000 bytes but not normally distributed')
disp('Press a key to continue...')
pause
fprintf('\n')

data4 = uint8(256*rand(1,1000).^3);
[zipped4,info4] = norm2huff(data4);
unzipped4 = huff2norm(zipped4,info4);

isOK = isequal(data4,unzipped4)

whos data4 zipped4 unzipped4

bar(frequency(data4)); axis tight

fprintf('\n')
disp('As the distribution is not uniform, the compression ratio is a bit more significant')
disp('Press a key to continue...')
pause
fprintf('\n')

% TERMINATED
fprintf('\n')
disp('#### FINISH ####')


text:
function dy = scopy ( n, dx, incx, dy, incy )

%*****************************************************************************80
%
%% SCOPY copies a vector X to a vector Y.
%
%  Discussion:
%
%    The routine uses unrolled loops for increments equal to one.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 May 2006
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer N, the number of elements in DX and DY.
%
%    Input, real DX(*), the first vector.
%
%    Input, integer INCX, the increment between successive entries of DX.
%
%    Input, real DY(*), the second vector, into which elements are copied.
%
%    Input, integer INCY, the increment between successive entries of DY.
%
%    Output, real DY(*), the second vector, with elements copied from DX.
%
  if ( n <= 0 )
    dy = [];
    return
  end

  if ( incx == 1 & incy == 1 )

    m = mod ( n, 7 );

    if ( m ~= 0 )

      dy(1:m) = dx(1:m);

    end

    for i = m+1 : 7 : n
      dy(i) = dx(i);
      dy(i + 1) = dx(i + 1);
      dy(i + 2) = dx(i + 2);
      dy(i + 3) = dx(i + 3);
      dy(i + 4) = dx(i + 4);
      dy(i + 5) = dx(i + 5);
      dy(i + 6) = dx(i + 6);
    end

  else

    if ( 0 <= incx )
      ix = 1;
    else
      ix = ( -n + 1 ) * incx + 1;
    end

    if ( 0 <= incy )
      iy = 1;
    else
      iy = ( -n + 1 ) * incy + 1;
    end

    for i = 1 : n
      dy(iy) = dx(ix);
      ix = ix + incx;
      iy = iy + incy;
    end

  end

  return
end


text:
function [ fea, out ] = ex_heattransfer7( varargin )
%EX_HEATTRANSFER7 1D Transient heat diffusion with analytic solution.
%
%   [ FEA, OUT ] = EX_HEATTRANSFER7( VARARGIN ) Transient heat
%   diffusion problem with analytic solution. A 1 m rod is kept at
%   fixed temperature on one end and constant outward heat flux at the
%   other end as in the following illustration.
%
%              +---------- L=1m ----------+ T = 25
%            q_n = 1       T(t=0) = 25
%
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       hmax        scalar {0.1}           Grid cell size
%       sfun        string {sflag1}        Finite element shape function
%       solver      string fenics/{}       Use FEniCS or default solver
%       ischeme     scalar {2}/1/3         Time stepping scheme
%       tmax        scalar {0.2}           Maximum time
%       tstep       scalar {0.01}          Time step size
%       iplot       scalar {1}/0           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct
%
%   See also EX_HEATTRANSFER8.

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { 'hmax',     0.1;
            'sfun',     'sflag1';
            'solver',   '';
            'ischeme',  2;
            'tmax',     0.2;
            'tstep',    0.01;
            'nstbwe',   0;
            'iplot',    1;
            'tol',      1e-3;
            'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});

% Grid generation.
fea.grid = linegrid( round(1/opt.hmax), 0, 1 );

% Problem definition.
fea.sdim  = { 'x' };                      % Space coordinate name.
fea = addphys( fea, @heattransfer );      % Add heat transfer physics mode.
fea.phys.ht.sfun = { opt.sfun };          % Set shape function.

% Equation coefficients.
fea.phys.ht.eqn.coef{1,end} = 1;          % Density.
fea.phys.ht.eqn.coef{2,end} = 1;          % Heat capacity.
fea.phys.ht.eqn.coef{3,end} = 1;          % Thermal conductivity.
fea.phys.ht.eqn.coef{6,end} = { 25 };     % Initial temperature.

% Boundary conditions.
fea.phys.ht.bdr.sel = [ 4 1 ];
fea.phys.ht.bdr.coef{1,end} = { [] 25 };
fea.phys.ht.bdr.coef{4,end}{1}{1} = -1;

% Parse physics modes and problem struct.
fea = parsephys(fea);
fea = parseprob(fea);

% Compute solution.
if( strcmp(opt.solver,'fenics') )
  fea = fenics( fea, 'fid', opt.fid, ...
                'tstep', opt.tstep, 'tmax', opt.tmax, 'ischeme', opt.ischeme );
  tlist = fea.sol.t;
else
  [fea.sol.u, tlist] = solvetime( fea, 'fid', opt.fid, 'init', {'T0_ht'}, 'ischeme', opt.ischeme, ...
                                  'tmax', opt.tmax, 'tstep', opt.tstep, 'nstbwe', opt.nstbwe );
end

% Postprocessing.
T_ref = refsol( fea.grid.p', tlist(end) );
if( opt.iplot>0 )
  postplot( fea, 'surfexpr', 'T', 'axequal', 0 )
  title(['Temperature at t=',num2str(tlist(end))])
  xlabel('x')
  ylabel('T')

  hold on
  plot( fea.grid.p, T_ref, 'r--' )
end

% Error checking.
T_sol = evalexpr( 'T', fea.grid.p, fea );
out.err  = norm( abs(T_sol-T_ref)/T_ref );
out.pass = out.err<opt.tol;

if( nargout==0 )
  clear fea out
end

% -----------------------------------
function [ u ] = refsol( x, t, dtol )

if( nargin<3 )
  dtol = 1e-7;
end

u0  = x + 24;
bdo = true;
n   = 0;
while( bdo )

  n   = n + 1;
  u   = u0 + 8/(1-2*n)^2/pi^2*cos((n-1/2)*pi*x)*exp(-((n-1/2)^2*pi^2)*t);
  bdo = any( max( u - u0 ) > dtol );
  u0  = u;

  if( n>1e3 )
    warning( ['Reference solution did not converge to tolerance ',num2str(dtol)] )
    break
  end
end

author: precise-simulation

repo: featool-multiphysics


text:
%
% output = bilateralFilter( data, edge, ...
%                          edgeMin, edgeMax, ...
%                          sigmaSpatial, sigmaRange, ...
%                          samplingSpatial, samplingRange )
%
% Bilateral and Cross-Bilateral Filter using the Bilateral Grid.
%
% Bilaterally filters the image 'data' using the edges in the image 'edge'.
% If 'data' == 'edge', then it the standard bilateral filter.
% Otherwise, it is the 'cross' or 'joint' bilateral filter.
% For convenience, you can also pass in [] for 'edge' for the normal
% bilateral filter.
%
% Note that for the cross bilateral filter, data does not need to be
% defined everywhere.  Undefined values can be set to 'NaN'.  However, edge
% *does* need to be defined everywhere.
%
% data and edge should be of the greyscale, double-precision floating point
% matrices of the same size (i.e. they should be [ height x width ])
%
% data is the only required argument
%
% edgeMin and edgeMax specifies the min and max values of 'edge' (or 'data'
% for the normal bilateral filter) and is useful when the input is in a
% range that's not between 0 and 1.  For instance, if you are filtering the
% L channel of an image that ranges between 0 and 100, set edgeMin to 0 and
% edgeMax to 100.
%
% edgeMin defaults to min( edge( : ) ) and edgeMax defaults to max( edge( : ) ).
% This is probably *not* what you want, since the input may not span the
% entire range.
%
% sigmaSpatial and sigmaRange specifies the standard deviation of the space
% and range gaussians, respectively.
% sigmaSpatial defaults to min( width, height ) / 16
% sigmaRange defaults to ( edgeMax - edgeMin ) / 10.
%
% samplingSpatial and samplingRange specifies the amount of downsampling
% used for the approximation.  Higher values use less memory but are also
% less accurate.  The default and recommended values are:
%
% samplingSpatial = sigmaSpatial
% samplingRange = sigmaRange
%

% Copyright (c) 2007 Jiawen Chen, Sylvain Paris, and Fredo Durand
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THE SOFTWARE.

function output = bilateralFilter( data, edge, edgeMin, edgeMax, sigmaSpatial, sigmaRange, ...
    samplingSpatial, samplingRange )

if( ndims( data ) > 2 ),
    error( 'data must be a greyscale image with size [ height, width ]' );
end

if( ~isa( data, 'double' ) ),
    error( 'data must be of class "double"' );
end

if ~exist( 'edge', 'var' ),
    edge = data;
elseif isempty( edge ),
    edge = data;
end

if( ndims( edge ) > 2 ),
    error( 'edge must be a greyscale image with size [ height, width ]' );
end

if( ~isa( edge, 'double' ) ),
    error( 'edge must be of class "double"' );
end

inputHeight = size( data, 1 );
inputWidth = size( data, 2 );

if ~exist( 'edgeMin', 'var' ),
    edgeMin = min( edge( : ) );
    warning( 'edgeMin not set!  Defaulting to: %f\n', edgeMin );
end

if ~exist( 'edgeMax', 'var' ),
    edgeMax = max( edge( : ) );
    warning( 'edgeMax not set!  Defaulting to: %f\n', edgeMax );
end

edgeDelta = edgeMax - edgeMin;

if ~exist( 'sigmaSpatial', 'var' ),
    sigmaSpatial = min( inputWidth, inputHeight ) / 16;
    fprintf( 'Using default sigmaSpatial of: %f\n', sigmaSpatial );
end

if ~exist( 'sigmaRange', 'var' ),
    sigmaRange = 0.1 * edgeDelta;
    fprintf( 'Using default sigmaRange of: %f\n', sigmaRange );
end

if ~exist( 'samplingSpatial', 'var' ),
    samplingSpatial = sigmaSpatial;
end

if ~exist( 'samplingRange', 'var' ),
    samplingRange = sigmaRange;
end

if size( data ) ~= size( edge ),
    error( 'data and edge must be of the same size' );
end

% parameters
derivedSigmaSpatial = sigmaSpatial / samplingSpatial;
derivedSigmaRange = sigmaRange / samplingRange;

paddingXY = floor( 2 * derivedSigmaSpatial ) + 1;
paddingZ = floor( 2 * derivedSigmaRange ) + 1;

% allocate 3D grid
downsampledWidth = floor( ( inputWidth - 1 ) / samplingSpatial ) + 1 + 2 * paddingXY;
downsampledHeight = floor( ( inputHeight - 1 ) / samplingSpatial ) + 1 + 2 * paddingXY;
downsampledDepth = floor( edgeDelta / samplingRange ) + 1 + 2 * paddingZ;

gridData = zeros( downsampledHeight, downsampledWidth, downsampledDepth );
gridWeights = zeros( downsampledHeight, downsampledWidth, downsampledDepth );

% compute downsampled indices
[ jj, ii ] = meshgrid( 0 : inputWidth - 1, 0 : inputHeight - 1 );

% ii =
% 0 0 0 0 0
% 1 1 1 1 1
% 2 2 2 2 2

% jj =
% 0 1 2 3 4
% 0 1 2 3 4
% 0 1 2 3 4

% so when iterating over ii( k ), jj( k )
% get: ( 0, 0 ), ( 1, 0 ), ( 2, 0 ), ... (down columns first)

di = round( ii / samplingSpatial ) + paddingXY + 1;
dj = round( jj / samplingSpatial ) + paddingXY + 1;
dz = round( ( edge - edgeMin ) / samplingRange ) + paddingZ + 1;

% perform scatter (there's probably a faster way than this)
% normally would do downsampledWeights( di, dj, dk ) = 1, but we have to
% perform a summation to do box downsampling
for k = 1 : numel( dz ),

    dataZ = data( k ); % traverses the image column wise, same as di( k )
    if ~isnan( dataZ  ),

        dik = di( k );
        djk = dj( k );
        dzk = dz( k );

        gridData( dik, djk, dzk ) = gridData( dik, djk, dzk ) + dataZ;
        gridWeights( dik, djk, dzk ) = gridWeights( dik, djk, dzk ) + 1;

    end
end

% make gaussian kernel
kernelWidth = 2 * derivedSigmaSpatial + 1;
kernelHeight = kernelWidth;
kernelDepth = 2 * derivedSigmaRange + 1;

halfKernelWidth = floor( kernelWidth / 2 );
halfKernelHeight = floor( kernelHeight / 2 );
halfKernelDepth = floor( kernelDepth / 2 );

[gridX, gridY, gridZ] = meshgrid( 0 : kernelWidth - 1, 0 : kernelHeight - 1, 0 : kernelDepth - 1 );
gridX = gridX - halfKernelWidth;
gridY = gridY - halfKernelHeight;
gridZ = gridZ - halfKernelDepth;
gridRSquared = ( gridX .* gridX + gridY .* gridY ) / ( derivedSigmaSpatial * derivedSigmaSpatial ) + ( gridZ .* gridZ ) / ( derivedSigmaRange * derivedSigmaRange );
kernel = exp( -0.5 * gridRSquared );

% convolve
blurredGridData = convn( gridData, kernel, 'same' );
blurredGridWeights = convn( gridWeights, kernel, 'same' );

% divide
blurredGridWeights( blurredGridWeights == 0 ) = -2; % avoid divide by 0, won't read there anyway
normalizedBlurredGrid = blurredGridData ./ blurredGridWeights;
normalizedBlurredGrid( blurredGridWeights < -1 ) = 0; % put 0s where it's undefined

% for debugging
% blurredGridWeights( blurredGridWeights < -1 ) = 0; % put zeros back

% upsample
[ jj, ii ] = meshgrid( 0 : inputWidth - 1, 0 : inputHeight - 1 ); % meshgrid does x, then y, so output arguments need to be reversed
% no rounding
di = ( ii / samplingSpatial ) + paddingXY + 1;
dj = ( jj / samplingSpatial ) + paddingXY + 1;
dz = ( edge - edgeMin ) / samplingRange + paddingZ + 1;

% interpn takes rows, then cols, etc
% i.e. size(v,1), then size(v,2), ...
output = interpn( normalizedBlurredGrid, di, dj, dz );

author: vitorsr

repo: SIHR

save_path: github-repos/MATLAB/vitorsr-SIHR


text:
function suborder = tetrahedron_ncc_suborder ( rule, suborder_num )

%*****************************************************************************80
%
%% TETRAHEDRON_NCC_SUBORDER returns the suborders for an NCC rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Peter Silvester,
%    Symmetric Quadrature Formulae for Simplexes,
%    Mathematics of Computation,
%    Volume 24, Number 109, January 1970, pages 95-100.
%
%  Parameters:
%
%    Input, integer RULE, the index of the rule.
%
%    Input, integer SUBORDER_NUM, the number of suborders of the rule.
%
%    Output, integer SUBORDER(SUBORDER_NUM), the suborders of the rule.
%
  if ( rule == 1 )
    suborder(1:suborder_num) = [ ...
      1 ];
  elseif ( rule == 2 )
    suborder(1:suborder_num) = [ ...
      4 ];
  elseif ( rule == 3 )
    suborder(1:suborder_num) = [ ...
      4, 6 ];
  elseif ( rule == 4 )
    suborder(1:suborder_num) = [ ...
      4, 12, 4 ];
  elseif ( rule == 5 )
    suborder(1:suborder_num) = [ ...
      4, 12, 6, 12, 1 ];
  elseif ( rule == 6 )
    suborder(1:suborder_num) = [ ...
      4, 12, 12, 12, 12, 4 ];
  elseif ( rule == 7 )
    suborder(1:suborder_num) = [ ...
      4, 12, 12, 12, 6, 24, 4, 4, 6 ];

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'TETRAHEDRON_NCC_SUBORDER - Fatal error!\n' );
    fprintf ( 1, '  Illegal RULE = %d\n', rule );
    error ( 'TETRAHEDRON_NCC_SUBORDER - Fatal error!' )

  end

  return
end

author: johannesgerer


text:
function drawCamera()

close all

addpath ../icosahedron2sphere

points = icosahedron2sphere(1);

points = points(points(:,3)>=0,:);

%{
plot3(points(:,1),points(:,2),points(:,3),'.')
title(sprintf('Level %d with %d points',1,size(points,1)))
axis equal
axis tight
xlabel('x');
ylabel('y');
zlabel('z');
%}

aspect_ratio = 4/3;
focal_length = 0.12;
h_fov = 54.4/180*pi;

w = tan(h_fov/2)*focal_length;
h = w/aspect_ratio;

camera=[...
0 -w           +w            +w            -w
0 -h           -h            +h            +h
0 focal_length focal_length  focal_length  focal_length];

%plotCamera(camera);

for i=1:size(points,1)
    center_ray = -points(i,:);

    projLen = sqrt(center_ray(1)^2+center_ray(2)^2);

    height = -projLen^2 / center_ray(3);

    upVector = [center_ray(1) center_ray(2) height];

    upVector = upVector /norm(upVector);

    %sinTheta = center_ray(3);

    leftVector = cross(center_ray,upVector);
    leftVector = leftVector/norm(leftVector);

    R = [leftVector' upVector' center_ray'];

    currentCamera = R * camera + repmat(points(i,:)',1,5);

    plotCamera(currentCamera);
end

axis equal
axis([-1 1 -1 1 -1 1]);
xlabel('x');
ylabel('y');
zlabel('z');

function plotCamera(camera)

%mid_ray = 1;
%side_rays = 1.2;

% frame
plot3(camera(1,[2 3 4 5 2]),camera(2,[2 3 4 5 2]),camera(3,[2 3 4 5 2]),'-k'); hold on;

% side rays
for i=2:5
    plot3(camera(1,[1 i]),camera(2,[1 i]),camera(3,[1 i]),'-k'); hold on;
end

author: jianxiongxiao

repo: ProfXkit