text:
classdef SubcellsMesher < handle

    properties (GetAccess = public, SetAccess = protected)
        subcells
    end

    properties (GetAccess = protected, SetAccess = protected)
        coord_iso
        coord_global
        levelSet
        connec
    end

    properties (Access = private)
    end

    properties (GetAccess = protected, SetAccess = private)
        mesh
        interior_coord_iso
        cutPoints
        nCutPoints
        cell_nodes
        cell_connec
        cell_levelSet
        boundary_levelSet
        coordsBackground
        posNodes
        levelSetBackground

    end

    methods (Access = protected, Abstract)

        computeCoordinates(obj)
        computeConnectivities(obj)
        computeLevelSet(obj)

    end

    methods (Access = public, Static)

        function obj = create(cParams)
            f = SubcellsMesherFactory();
            obj = f.create(cParams);
        end

    end

    methods (Access = public)

        function computeSubcells(obj,cParams)

            cellConnec = cParams.cellConnec;
            cutPoints  = cParams.cutPoints;

            obj.saveInputs(cellConnec,cutPoints);
            obj.computeLevelSet();
            obj.computeCoordinates();
            obj.computeConnectivities();
            obj.packResultsInStruct();
        end

    end

    methods (Access = protected)

        function init(obj,cParams)
            obj.posNodes = cParams.posNodes;
            obj.levelSetBackground = cParams.levelSetBackground;
            obj.coordsBackground = cParams.coordsBackground;
        end
    end

    methods (Access = protected, Static)

        function connectivities = computeDelaunay(coordinates)
            DT = delaunayTriangulation(coordinates);
            connectivities = DT.ConnectivityList;
        end

        function w = patch(u,v)
            w = [u;v];
        end

    end

    methods (Access = private)

        function saveInputs(obj,cell_connec,cutPoints)
            obj.cell_connec = cell_connec;
            obj.cutPoints = cutPoints;
        end

        function packResultsInStruct(obj)
            obj.subcells.levelSet     = obj.levelSet;
            obj.subcells.coord_iso    = obj.coord_iso;
            obj.subcells.coord_global = obj.coord_global;
            obj.subcells.connec       = obj.connec;

            obj.subcells.nNodes       = size(obj.coord_iso,1);
            obj.subcells.nSubcells    = size(obj.connec,1);
        end

    end

    methods

        function coord = get.interior_coord_iso(obj)
            coord = obj.patch(obj.cell_nodes,obj.cutPoints.ISO);
        end

        function n = get.nCutPoints(obj)
            n = size(obj.cutPoints.ISO,1);
        end

        function nodes = get.cell_nodes(obj)
            nodes = obj.posNodes;
        end

        function lvlSet = get.cell_levelSet(obj)
            lvlSet = obj.levelSetBackground(obj.cell_connec);
        end

        function lvlSet = get.boundary_levelSet(obj)
            lvlSet = zeros(obj.nCutPoints,1);
        end

    end
end

author: SwanLab

repo: Swan


text:
%kcpRGBNTSCtoI1I2I3 'converts an image from RGB NTSC to I1I2I3 '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros cpRGBNTSCtoI1I2I3.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'First Input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kcpRGBNTSCtoI1I2I3(i, {'i','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% cpRGBNTSCtoI1I2I3 - converts an image from RGB NTSC to I1I2I3
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1997 Rafael Santos. Khoros (C) Khoral Research, Inc.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function item = clearval(item, dflag)

% function item = clearval(item, dflag)
% This is a generic function to clear the contents of the val field of a
% cfg_item. It is usable for all leaf cfg_item classes (cfg_entry,
% cfg_files, cfg_menu).
% dflag is ignored for leaf entries - the val field is cleared
% unconditionally.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: clearval.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

item.val = {};

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function out=bgsubtractbyhand(in,peval)
% out=bgsubtractbyhand(in,peval.bg)
mfprintf(peval.fid,'Background (%g) subtracted and clipped by hand.\n',peval.bg)
out=max(in-peval.bg,eps);


text:
clear ra
for i = 1 : length(list)
    s = size(imread(list{i}));
    s = max(s);
    t=gt{i};
    t = t(:,3)-t(:,1);
    t=t/s*2048;
    ra{i} = t;
end
t = cellfun(@(x) unique(11-floor(log2(x))), ra, 'uni', 0);
 class = cellfun(@(x) max(1, min(ceil(log2(x   )*10)-50, 60)), ra, 'uni', 0);

author: liuyuisanai

repo: RSA-for-object-detection


text:
function d =  generate(a)

if a.seed>0 rand('seed',a.seed);   end;
x=(rand(a.l,a.n)-0.5)*2;
w=ones(1,a.n);
prem=a.n-min(a.relevant_n,a.n);
w(1:prem)=0;
y=sign(x*w');

d=data([get_name(a) ' '] ,x,y);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function teacher_stats(varargin)
%TEACHER_STATS - compute the distribution of teacher predictions
%   TEACHER_STATS(VARARGIN) computes the distribution of emotion predictions
%   made by the teahcer model on the EmoVoxCeleb dataset and the Afew6.0
%   dataset [1] and combines them in a figure.
%
%   TEACHER_STATS(..'name', value) accepts the following options:
%
%   `figurePath` :: fullfile(vl_rootnn, 'data/xEmo18/emovoxceleb-figure.pdf')
%    Path to location where figure will be stored.
%
%    References:
%     [1] Dhall, Abhinav, et al. "Collecting large, richly annotated
%       facial-expression databases from movies." IEEE multimedia 19.3%
%       (2012): 34-41.
%
% Copyright (C) 2018 Samuel Albanie, Arsha Nagrani
% Licensed under The MIT License [see LICENSE.md for details]

  opts.figurePath = fullfile(vl_rootnn, ...
                       'data/emoVoxCeleb/emovoxceleb-figure.pdf') ;
  opts.afewLogits = fullfile(vl_rootnn, 'data/emoVoxCeleb/afew-logits.mat') ;
  opts = vl_argparse(opts, varargin) ;

  fprintf('loading imdb of predictions on EmoVoxCeleb...') ; tic ;
  opts.teacher = 'senet50-ferplus' ;
  imdb = fetch_emovoxceleb_imdb(opts.teacher) ;
  allLogits = vertcat(imdb.wavLogits{:}) ;
  [~, emoVoxPreds] = max(allLogits, [], 2) ;
  fprintf('done in %g (s) \n', toc) ;

  fprintf('loading imdb of teacher predictions on AFEW...') ; tic ;
  % Note: These are stored online as a convenience but can be recomputed
  % if desired.
  if ~exist(fileparts(opts.afewLogits), 'dir')
		mkdir(fileparts(opts.afewLogits)) ;
	end
	fetchLogitsFromInternet(opts.afewLogits, 'afew-logits') ;
  xx = load(opts.afewLogits) ;
  emos = vertcat(xx.faceLogits{:}) ;
  [~, afewPreds] = max(emos, [], 2) ;
  fprintf('done in %g (s) \n', toc) ;
  plotHistogram(emoVoxPreds, afewPreds, opts.figurePath) ;
end

% --------------------------------------------------------------------
function plotHistogram(emoVoxPreds, afewPreds, figurePath)
% --------------------------------------------------------------------

  ferPlusEmotions = {'Neutral', 'Happiness', 'Surprise', 'Sadness', ...
                    'Anger', 'Disgust', 'Fear', 'Contempt'} ;

  olive = [ 125 ; 168 ; 50 ] / 255 ; % pick a colour scheme
  set(0,'defaulttextinterpreter','latex') ;

  h = figure(1) ;
  [emoCeleb, ~] = histcounts(emoVoxPreds, 0.5:8.5) ;
  [compared, edges] = histcounts(afewPreds, 0.5:8.5) ;
  label1 = 'EmoVoxCeleb' ;
  label2 = 'Afew 6.0' ;
	ctrs = edges(1)+(1:length(edges)-1).*diff(edges);   % Create Centres
	b = bar(ctrs, [emoCeleb ; compared]', 1, 'BarWidth', 1) ;
  b(1).EdgeColor = olive ; b(1).FaceColor = olive ; b(1).FaceAlpha = 0.8 ;
  b(2).EdgeColor = 'red'; b(2).FaceColor = 'red'; b(2).FaceAlpha = 0.8 ;

	set(gca,'YScale','log') ; grid on ;
	set(gca, 'xtick', ctrs, 'xticklabel', ferPlusEmotions) ;
  ylabel('Number of frames') ;
  set(groot, 'defaultAxesTickLabelInterpreter','latex') ;
  set(groot, 'defaultLegendInterpreter','latex') ;
  legend({label1, label2}) ;
  a = get(gca,'XTickLabel') ;
  set(gca,'XTickLabel', a, 'fontsize', 15) ;
  xlim([1 9]) ;
	xtickangle(60) ;

  if exist('zs_dispFig', 'file'), zs_dispFig ; end

  set(h,'Units','Inches');
  pos = get(h,'Position');
  set(h,'PaperPositionMode', 'Auto', ...
         'PaperUnits', 'Inches', 'PaperSize',[pos(3), pos(4)]) ;
  print(figurePath, '-dpdf','-r0') ;
end

% -------------------------------------------------
function fetchLogitsFromInternet(destPath, imdbName)
% -------------------------------------------------

  waiting = true ;
  prompt = sprintf(...
        strcat('Logits were not found at %s\nWould you like to ', ...
        ' download them from THE INTERNET (y/n)?\n'), destPath) ;

  baseUrl = 'http://www.robots.ox.ac.uk/~albanie/data/cross-modal-emotions' ;

  switch imdbName
    case 'afew-logits'
      url = [baseUrl '/afew-logits.mat'] ;
    otherwise
      error('did not recognise imdb name %s\n', imdbName) ;
  end

  if exist(destPath, 'file')
    fprintf('file already exists at destination, skipping...\n') ;
    return
  end

  while waiting
    str = input(prompt,'s') ;
    switch str
      case 'y'
        if ~exist(fileparts(destPath), 'dir'), mkdir(fileparts(destPath)) ; end
        fprintf('VoxCelebImdb ... \n') ;
        urlwrite(url, destPath) ;
        return ;
      case 'n', throw(exception) ;
      otherwise, fprintf('input %s not recognised, please use `y/n`\n', str) ;
    end
  end
end

author: albanie

repo: mcnCrossModalEmotions


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%
%
% problem 2- Evaluate the Laplace Transform of  y''(t) and then replace y(t) by sin(t)u(t)

sym 'y(t)' ;
syms t s
z=laplace( diff('y(t)',2),s)

z=subs(z,'y(t)',sin(t))

z=subs(z,'y(0)',0)
z=subs(z,'D(y)(0)',1)
simplify(z)

%verification
x=sin(t);
x2=diff(x,2,t)
laplace(x2,s)

author: Sable

repo: mcbench-benchmarks


text:
% Support only 2 dim array
classdef CMatrix
   properties (Constant)
      mex = CMatrix.mexSelector();
   end

   methods (Static)
      function func = mexSelector()
         if isarm()
            func = @CMatrixArmMex;
         else
            func = @CMatrixMex;
         end
      end
   end

   properties
      x;
   end

   methods
      function o = CMatrix(a)
         if iscell(a) || isa(a, 'uint8')
            o.x = a;
         else
            o.x = CMatrix.toMex(a);
         end
      end

      %% Sizes
      function [m, n] = size(a, dim)
         a = a.x;
         if (iscell(a))
            s = double(size(a{1}));
         else
            s = [size(a, 2) size(a, 3)];
         end

         if (nargin == 1)
            if (nargout == 2)
               m = s(1);
               n = s(2);
            else
               m = s;
            end
         else
            m = s(dim);
         end
      end

      function r = length(a)
         s = size(a);
         if (min(s) == 0)
            r = 0;
         else
            r = max(s);
         end
      end

      function r = numel(a)
         r = prod(size(a)); %#ok<*PSIZE>
      end

      function r = isscalar(a)
         r = (numel(a) == 1);
      end

      function r = isvector(a)
         r = (size(a,1) == 1 || size(a,2) == 1);
      end

      function r = isempty(a)
         r = (numel(a) == 0);
      end

      function r = isrow(a)
         r = (size(a,1) == 1);
      end

      function r = iscolumn(a)
         r = (size(a,2) == 1);
      end

      function r = ismatrix(~)
         r = true;
      end

      function r = end(a, k, n)
         if n == 2
            r = size(a,k);
         else
            r = size(a,1) * size(a,2);
         end
      end

      %% Conversions
      function disp(a)
         disp(double(a))
      end

      function r = double(a)
         r = CMatrix.UnaryOp('double', a);
      end

      function r = logical(a)
         r = CMatrix.UnaryOp('logical', a);
      end

      function r = sparse(a,b,c,d,e,f)
         if nargin == 1
            if ~issparse(a)
               A = logical(sparse(double(a)));
               s = struct('type', '()');
               s.subs = {find(A)};
               s = CMatrix(subsref(a, s));
               s.x = reshape(s.x, [size(s.x,1) numel(s.x)/size(s.x,1)]);
               r = CMatrix({A;s.x});
            else
               r = CMatrix(a);
            end
         else
            assert(nargin == 3 || nargin == 5 || nargin == 6, "unsupported number of arguments for sparse");
            if nargin == 3
               i = a; j = b; s = c; m = max(i); n = max(j); nzmax = length(s);
               if isempty(i)
                  m = 0; n = 0;
               end
            elseif nargin == 5
               i = a; j = b; s = c; m = d; n = e; nzmax = length(s);
            elseif nargin == 6
               i = a; j = b; s = c; m = d; n = e; nzmax = f;
            end
            A = sparse(i, j, ones(nzmax,1,'logical'),m,n,nzmax);
            assert(length(s) == nzmax);
            s = CMatrix(s(:));
            r = CMatrix({A;s.x});
         end
      end

      function r = full(a)
         if issparse(a)
            A = a.x{1};
            s = a.x{2};
            r = zeros(size(s,1), size(A,1), size(A,2), 'uint8');
            idx = find(a.x{1});
            r(:, idx) = s;
            r = CMatrix(r);
         else
            r = CMatrix(a);
         end
      end

      function r = issparse(a)
         r = iscell(a.x);
      end

      function r = issymmetric(a)
         if size(a,1) ~= size(a,2)
            r = false;
         else
            r = full(all(a == a', 'all'));
         end
      end

      function r = diag(a)
         assert(nargin == 1, "Support only one parameters")
         if isvector(a)
            r = CMatrix(sparse(a));
            if size(a,1) == 1
               r = r';
            end
            r.x{1} = diag(double(a));
            if ~issparse(a)
               r = full(r);
            end
         else
            k = min(size(a));
            if k == 0
               r = CMatrix([]);
            else
               idx = (1:k) + size(a,1) * (0:(k-1));
               r =  subsref(a, struct('type', '()', 'subs', {{idx}}));
            end
            if issparse(a)
               r = sparse(r);
            end
            if size(a,2) < size(a,1)
               r = r';
            end
         end
      end

      function r = reshape(a, s)
         if issparse(a)
            a.x{1} = reshape(a.x{1}, s);
         else
            a.x = reshape(a.x, [size(a.x,1) s]);
         end
         r = CMatrix(a.x);
      end

      function r = eps(~)
         r = CMatrix(CMatrix.mex('eps'));
      end

      %% Comparisons
      function r = lt(a, b)
         r = CMatrix.BinaryOp('lt', a, b);
         if ((issparse(a) || issparse(b)) && ~issparse(r)), r = sparse(r); end
      end

      function r = gt(a, b)
         r = CMatrix.BinaryOp('gt', a, b);
         if ((issparse(a) || issparse(b)) && ~issparse(r)), r = sparse(r); end
      end

      function r = le(a, b)
         r = ~gt(a,b);
      end

      function r = ge(a, b)
         r = ~lt(a,b);
      end

      function r = ne(a, b)
         r = CMatrix.BinaryOp('ne', a, b);
         if ((issparse(a) || issparse(b)) && ~issparse(r)), r = sparse(r); end
      end

      function r = eq(a, b)
         r = ~ne(a,b);
      end

      function r = and(a, b)
         r = CMatrix.BinaryOp('and', a, b);
      end

      function r = or(a, b)
         r = CMatrix.BinaryOp('or', a, b);
      end

      function r = not(a)
         r = ~logical(a);
      end

      %% Arithmetic
      function r = plus(a, b)
         r = CMatrix.BinaryOp('plus', a, b);
      end

      function r = minus(a, b)
         r = CMatrix.BinaryOp('minus', a, b);
      end

      function r = times(a, b)
         r = CMatrix.BinaryOp('times', a, b);
      end

      function r = power(a, b)
         assert(isscalar(b) && b == 2, "support only power(a, 2)");
         r = a .* a;
      end

      function r = rdivide(a, b)
         if issparse(b)
            b = full(b);
         end
         r = CMatrix.BinaryOp('rdivide', a, b);

         if issparse(a) % To fix the bug regarding x / 0
            r2 = double(a) ./ double(b);
            r(~isfinite(r2)) = r2(~isfinite(r2));
         end
      end

      function r = ldivide(a, b)
         r = b./a;
      end

      function r = uminus(a)
         r = CMatrix.UnaryOp('uminus', a);
      end

      function r = uplus(a)
         r = a;
      end

      function r = abs(a)
         r = CMatrix.UnaryOp('abs', a);
      end

      function r = sqrt(a)
         r = CMatrix.UnaryOp('sqrt', a);
      end

      function r = sum(varargin)
         r = CMatrix.ReductionOp('sum', varargin{:});
      end

      function r = prod(varargin)
         r = CMatrix.ReductionOp('prod', varargin{:});
      end

      function r = all(a, dim)
         if nargin == 1
            r = all(logical(a));
         else
            r = all(logical(a), dim);
         end
      end

      function r = any(a, dim)
         if nargin == 1
            r = any(logical(a));
         else
            r = any(logical(a), dim);
         end
      end

      function r = max(a, b, c)
         if nargin == 1
            r = CMatrix.ReductionOp('max', a);
         elseif nargin == 3 && size(b,1) == 0 && size(b,2) == 0
            r = CMatrix.ReductionOp('max', a, c);
         else
            assert(nargin == 2);
            r = CMatrix.BinaryOp('max2', a, b);
            if issparse(a) || issparse(b)
               r = sparse(r);
            end
         end
      end

      function r = min(a, b, c)
         if nargin == 1
            r = CMatrix.ReductionOp('min', a);
         elseif nargin == 3 && size(b,1) == 0 && size(b,2) == 0
            r = CMatrix.ReductionOp('min', a, c);
         else
            assert(nargin == 2);
            r = CMatrix.BinaryOp('min2', a, b);
            if issparse(a) || issparse(b)
               r = sparse(r);
            end
         end
      end

      function r = norm(a)
         assert(isvector(a), "norm supports only vectors.");
         r = sqrt(full(sum(a.^2)));
      end

      function r = nnz(a)
         r = full(sum(logical(a),'all'));
      end

      %% Access
      function [i,j,v] = find(a)
         if nargout == 1
            i = find(logical(a));
         elseif nargout == 2
            [i,j] = find(logical(a));
         else
            a = sparse(CMatrix(a));
            [i,j] = find(a.x{1});
            v = CMatrix(a.x{2});
            if size(a,1) == 1
               v = v';
            end
            if numel(v) == 0
               v = CMatrix([]);
            end
         end
      end

      %% Matrix operations
      function r = ctranspose(a)
         r = transpose(a);
      end

      function r = transpose(a)
         r = CMatrix.UnaryOp('transpose', a);
      end

      function r = mtimes(a, b)
         if isscalar(a) || isscalar(b)
            r = a .* b;
         else
            r = CMatrix.BinaryOp('mtimes', a, b);
         end
      end

      function r = mldivide(a, b)
         assert(size(a,1) == size(b,1), 'Incompatible Size.');
         s = [size(a,2), size(b,2)];
         if size(a,2) == 0 || size(a,1) == 0 || size(b,2) == 0
            r = CMatrix(zeros(s));
         elseif isscalar(a)
            r = reshape(b ./ a, s);
         else
            b_full = full(b); % sparsity usually does not help. Do this for simplicity
            r = CMatrix.BinaryOp('mldivide', a, b_full);
         end
         if (issparse(a) && issparse(b))
            r = sparse(r);
         end
      end

      function r = mrdivide(a, b)
         r = (b'\a')';
      end

      function r = horzcat(varargin)
         allDense = true;
         for k = 1:length(varargin)
            allDense = allDense & ~issparse(varargin{k});
         end

         if (allDense)
            cat_input = cell(length(varargin), 1);
            for k = 1:length(varargin)
               cat_input{k} = varargin{k}.x;
            end
            r = CMatrix(cat(3, cat_input{:}));
         else
            A_bool = [];
            m = []; v = [];
            for k = 1:length(varargin)
               A = sparse(CMatrix(varargin{k}));
               if k == 1
                  m = size(A,1);
                  A_bool = A.x{1};
                  v = A.x{2};
               else
                  assert(m == size(A,1), 'Incompatible size.');
                  A_bool = [A_bool A.x{1}];
                  v = [v A.x{2}];
               end
            end
            r = CMatrix({A_bool;v});
         end
      end

      function r = vertcat(varargin)
         for i = 1:length(varargin)
            varargin{i} = varargin{i}';
         end
         r = horzcat(varargin{:})';
      end

      function r = subsref(a,s)
         if strcmp(s(1).type, '()')
            if issparse(a)
               idx = subsref(CMatrix.ToIndex(a),s);
               A = logical(idx);
               [~,~,idx] = find(idx);
               v = a.x{2};
               v = v(:, idx);
               r = CMatrix({A;v});
            else
               idx = subsref(CMatrix.ToIndex(a),s);
               r = CMatrix(reshape(a.x(:,idx), [size(a.x,1) size(idx)]));
            end
         else
            r = builtin('subsref',a,s);
         end
      end

      function r = subsasgn(a,s,b)
         if strcmp(s(1).type, '()')
            [a_idx, a_v] = CMatrix.ToIndex(a);
            [b_idx, b_v] = CMatrix.ToIndex(b);
            b_idx = max(a_idx,[],'all') + b_idx;
            v = [a_v b_v];

            r_idx = subsasgn(a_idx,s,b_idx);

            if issparse(r_idx)
               A = logical(r_idx);
               [~,~,idx] = find(r_idx);
               v = v(:, idx);
               r = {A;v};
            else
               r = reshape(v(:, r_idx), [size(v,1) size(r_idx)]);
            end
            r = CMatrix(r);
         else
            r = builtin('subsasgn',a,s,b);
         end
      end

      %% cholesky related
      function r = dissect(a)
         r = dissect(double(a));
      end

      function r = amd(a)
         r = amd(double(a));
      end

      function r = colamd(a)
         r = colamd(double(a));
      end

      function r = symbfact(a)
         r = symbfact(double(a));
      end

      function [r1, r2] = etree(a)
         [r1, r2] = etree(double(a));
      end

      function r = chol(a)
         r = CMatrix.UnaryOp('chol', a);
      end
   end

   methods (Static)
      function r = ones(varargin)
         r = CMatrix(ones(varargin{:}));
      end

      function r = zeros(varargin)
         r = CMatrix(zeros(varargin{:}));
      end

      function r = eye(varargin)
         r = CMatrix(eye(varargin{:}));
      end

      function r = rand(varargin)
         r = CMatrix(rand(varargin{:}));
      end

      function r = randn(varargin)
         r = CMatrix(randn(varargin{:}));
      end

      function r = randi(varargin)
         r = CMatrix(randi(varargin{:}));
      end

      function r = toMex(a)
         if isa(a, 'CMatrix')
            r = a.x;
         elseif isobject(a)
            r = a.x;
            r = CMatrix.mex('toMex', r);
         else
            r = CMatrix.mex('toMex', double(a));
         end
      end

      function r = UnaryOp(cmd, a)
         a = CMatrix.toMex(a);
         r = CMatrix.mex(cmd, a);

         if iscell(r) || isa(r, 'uint8')
            r = CMatrix(r);
         end
      end

      function r = BinaryOp(cmd, a, b)
         a = CMatrix.toMex(a);
         b = CMatrix.toMex(b);
         r = CMatrix.mex(cmd, a, b);

         if iscell(r) || isa(r, 'uint8')
            r = CMatrix(r);
         end
      end

      function [r, v] = ToIndex(a)
         if nargout == 2 && ~isa(a, 'CMatrix')
            a = CMatrix(a);
         end

         if issparse(a)
            [i,j] = find(double(a));
            r = sparse(i,j,1:length(i),size(a,1),size(a,2));

            if nargout == 2
               v = a.x{2};
            end
         else
            r = reshape(1:numel(a),[size(a,1) size(a,2)]);

            if nargout == 2
               v = reshape(a.x, [size(a.x,1) numel(a.x)/size(a.x,1)]);
            end
         end
      end

      function r = ReductionOp(cmd, a, dim_)
         a_size = size(a);
         a = CMatrix(a);

         if nargin == 2
            if (size(a,1) > 1)
               dim = 1;
            elseif (size(a,2) > 1)
               dim = 2;
            else
               dim = 1;
               a = subsref(a, struct('type', '()', 'subs', {{':'}}));
            end
         else
            assert(strcmp(dim_, 'all') || isscalar(dim_), 'dim can either be all or a scalar');
            if strcmp(dim_, 'all')
               dim = 1;
               a = subsref(a, struct('type', '()', 'subs', {{':'}}));
            else
               dim = dim_;
            end
         end

         if size(a,1) == 0
            switch cmd
               case 'sum'
                  r = CMatrix(0.0 * ones(1, size(a,2)));
               case 'prod'
                  r = CMatrix(1.0 * ones(1, size(a,2)));
               case 'max'
                  r = CMatrix(ones(0, a_size(2)));
               case 'min'
                  r = CMatrix(ones(0, a_size(2)));
            end
            if dim == 2
               r = r';
            end
         elseif dim == 1
            r = CMatrix.mex(cmd, a.x);
            r = CMatrix(r);
         elseif dim == 2
            a = a';
            r = CMatrix.mex(cmd, a.x);
            r = CMatrix(r)';
         end

         if issparse(a)
            r = sparse(r);
         end
      end
   end
end

author: opencobra

repo: cobratoolbox


text:
function M = ptv_params_to_matrix(params, transform_type)
    if isnumeric(transform_type)
        if transform_type == -1
            transform_type = 'translate';
        elseif transform_type == -2
            transform_type = 'rigid';
        elseif transform_type == -3
            transform_type = 'rigid_scale';
        elseif transform_type == -4
            transform_type = 'affine';
        end
    end

    if strcmp(transform_type, 'affine')
        M = [1 + params(3), params(4), params(1); ...
          params(5), 1 + params(6), params(2)];
    elseif strcmp(transform_type, 'rigid')
        theta = params(3);
        M = [cos(theta), -sin(theta), params(1); ...
             sin(theta), cos(theta), params(2)];
    elseif strcmp(transform_type, 'rigid_scale')
        theta = params(3);
        scale = 1 + params(4);
        M = [cos(theta) * scale, -sin(theta) * scale, params(1); ...
             sin(theta) * scale, cos(theta) * scale, params(2)];
    elseif strcmp(transform_type, 'translate')
        M = [1, 0, params(1); ...
             0, 1, params(2)];
    end
end

repo: pTVreg


text:
function [mg] = kg2mg(kg)
% Convert mass from kilograms to milligrams.
% Chad Greene 2012
mg = kg*1000000;


text:
classdef stats

methods(Static)

    function [frequencies, labels] = relative_frequencies(data)
        % Returns the frequencies of values of a discrete random variable in the dataset
        n = length(data);
        % Reshape data to a row vector
        data = spx.vector.reshape_as_row_vec(data);
        % Identify unique labels in the data set and their frequencies
        [labels, ~, label_map] = unique(data);
        % Lets get the number of unique labels
        nl = numel(labels);
        % Let's construct an n x nl sparse matrix in which
        % each column identifies c-th label
        % each row identifies the r-th data point
        % each row contains only one 1 and rest of the entries are 0.
        % The column in which 1 is stored identifies the label of the data entry
        % We store this information as the sparse matrix.
        % We already know that the matrix will have n non-zero entries
        non_zero_max = n;
        rows = 1:n;
        columns = label_map;
        index_to_label_matrix = sparse(rows, columns,  1,    n, nl,    non_zero_max);
        % Take mean over each column to obtain the frequencies
        frequencies = full(mean(index_to_label_matrix,1));
    end

    function [ statistic ] = compute_statistic_per_vector(...
        receivedSequence, N, statisticFunc)
        %COMPUTE_STATISTIC_PER_VECTOR Computes a statistic for each vector
        %   Detailed explanation goes here

        % Number of samples
        L = length(receivedSequence);
        % Number of vectors
        M = round(L / N);
        % Each row now represents one vector
        receivedSequence = reshape(receivedSequence, N, M)';
        statistic = zeros(M, 1);
        for i=1:M
            data = receivedSequence(i, :);
            statistic(i) = statisticFunc(data);
        end
    end

    function result = format_descriptive_statistics(x)
        mu = mean(x);
        med = median(x);
        rng = range(x);
        r = iqr(x);
        sigma = std(x);
        [m1, min_idx] = min(x);
        [m2, max_idx] = max(x);
        m3 = mad(x);
        result = sprintf('Min: %.2f (%d), Max: %.2f (%d), Mean: %.2f, Median: %.2f, Range : %.2f, IQR: %.2f, Deviation: %.2f, MAD: %.2f', m1, min_idx, m2, max_idx, mu, med, rng, r, sigma, m3);
    end

    function y = rand_subset(n, k)
        % Returns y has a column vector of k values sampled uniformly
        % at random without replacement from the integers 1:n
        rp = randperm(n);
        y = rp(1:k);
        y = y(:);
    end

end

end

author: indigits

repo: sparse-plex


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

%% -*- texinfo -*-
%% @defun  show_system_info ()
%% Show information about the current system.
%%
%% @seealso{computer, ver}
%% @end defun

function show_system_info ()
  [platform, array_maxsize, endian] = computer ();
  arch = computer ('arch');

  disp ('System info');
  disp ('-----------');
  disp ('');
  fprintf ('Platform: %s\n', platform);
  fprintf ('Array Maximum Size: %d\n', array_maxsize);
  fprintf ('Endianness: %s\n', endian);
  fprintf ('Architecture: %s\n', arch);
  ver ();
  disp ('');
end

repo: octsympy


text:
function M = spm_meanm(A)
% Compute barycentre of matrix exponentials
% FORMAT M = spm_meanm(A)
% A - A 3D array, where each slice is a matrix
% M - the resulting mean
%
% Note that matrices should not be too dissimilar to each other or the
% procedure fails.
% See http://hal.archives-ouvertes.fr/hal-00699361/
%__________________________________________________________________________
% Copyright (C) 2012-2019 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_meanm.m 7563 2019-04-01 10:39:24Z guillaume $

N = size(A,3);
M = eye(size(A,1),size(A,2));

for iter = 1:1024
    S = zeros(size(M));
    for i=1:N
        L = real(logm(M\A(:,:,i)));
        S = S + L;
    end
    S = S/N;
    M = M*expm(S);
    %imagesc(M); drawnow
    %fprintf('%d\t%g\n', iter,sum(S(:).^2));
    if sum(S(:).^2)<1e-20
        break;
    end
end

author: spm


text:
clear; close all; clc;

%% Fig 3. 한 점 (x,y)에서 순간 변화율과 평균 변화율 (동영상)
for x1 = [1, 0.5, 0.1]+1
    plotXY(-1,3,-1,3,true,20)
    x = linspace(-1,3, 100);
    y = 1/3*x.^2 + 1;
    hold on;
    plot(x,y,'k','linewidth',2)

    x = linspace(-0.5, 2.5, 100);

    plot(x, 2/3*x + 2/3,'b','linewidth',3)
    plot(1, 4/3 ,'o','markerfacecolor','b','markeredgecolor','none','markersize',10)

    y1 = 1/3 * x1^2  + 1;

    y = (y1-4/3)/(x1-1) * (x-1) + 4/3;
    h(1) = plot(x1, y1, 'o', 'markerfacecolor','r','markeredgecolor','none','markersize',10);
    h(2) = plot(x, y, 'r','linewidth',3);

end

%% Fig 3. 한 점 (x,y)에서 순간 변화율과 평균 변화율 (동영상)
clear v
v = VideoWriter('fig3.mp4','MPEG-4');
v.FrameRate = 30;
v.Quality = 100;
open(v);

clear F

plotXY(-1,3,-1,3,true,20)
x = linspace(-1,3, 100);
y = 1/3*x.^2 + 1;
hold on;
plot(x,y,'k','linewidth',2)

x = linspace(-0.5, 2.5, 100);

plot(x, 2/3*x + 2/3,'b','linewidth',3)
plot(1, 4/3 ,'o','markerfacecolor','b','markeredgecolor','none','markersize',10)

i=1;
for x1 = linspace(3,1,100)

    y1 = 1/3 * x1^2  + 1;

    y = (y1-4/3)/(x1-1) * (x-1) + 4/3;
    h(1) = plot(x1, y1, 'o', 'markerfacecolor','r','markeredgecolor','none','markersize',10);
    h(2) = plot(x, y, 'r','linewidth',3);
    drawnow;

    F(i)=getframe(gcf);
    if x1>1
        delete(h)
    end
    i=i+1;
end

for i=1:length(F)
    % convert the image to a frame
    frame = F(i) ;
    writeVideo(v, frame);
end

close(v)

%% 3d plot

close all;
[X,Y] = meshgrid(linspace(-3,3,20));

surf(X, Y, -X.^2-Y.^2+1)
xlabel('$$x$$','interpreter','latex'); ylabel('$$y$$','interpreter','latex'); zlabel('$$z=-x^2-y^2+1$$','interpreter','latex');
hold on;
xx = -1 * ones(1,100);
yy = linspace(-3,3,100);
zz = -xx.^2 - yy.^2  + 1;
plot3(xx,yy,zz,'r','linewidth',3)

% tangential line
% (-1, -1, -1)에서...?
% round z round y = - 2y
t = linspace(-0.5,1,100);

x = -1 * ones(1, 100);
y = -2 + 4 * t;
z = -4 + 16*t

plot3(x,y,z,'k','linewidth',3);
plot3(-1,-2,-4,'o','markersize',7,'markerfacecolor','k','markeredgecolor','none')

axis square

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
classdef TestSimpleBlobDetector
    %TestSimpleBlobDetector

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','tsukuba_l.png');
        kfields = {'pt', 'size', 'angle', 'response', 'octave', 'class_id'};
    end

    methods (Static)
        function test_detect_img
            obj = cv.SimpleBlobDetector('MinThreshold',50, 'MaxThreshold',220);
            typename = obj.typeid();

            img = imread(TestSimpleBlobDetector.im);
            kpts = obj.detect(img);
            validateattributes(kpts, {'struct'}, {'vector'});
            assert(all(ismember(TestSimpleBlobDetector.kfields, fieldnames(kpts))));
        end

        function test_detect_imgset
            img = imread(TestSimpleBlobDetector.im);
            imgs = {img, img};
            obj = cv.SimpleBlobDetector();

            kpts = obj.detect(imgs);
            validateattributes(kpts, {'cell'}, {'vector', 'numel',numel(imgs)});
            cellfun(@(kpt) validateattributes(kpt, {'struct'}, {'vector'}), kpts);
            cellfun(@(kpt) assert(all(ismember(TestSimpleBlobDetector.kfields, fieldnames(kpt)))), kpts);
        end

        function test_detect_mask
            img = imread(TestSimpleBlobDetector.im);
            mask = zeros(size(img), 'uint8');
            mask(:,1:end/2) = 255;  % only search left half of the image
            obj = cv.SimpleBlobDetector();
            kpts = obj.detect(img, 'Mask',mask);
            validateattributes(kpts, {'struct'}, {'vector'});
            assert(all(ismember(TestSimpleBlobDetector.kfields, fieldnames(kpts))));
            xy = cat(1, kpts.pt);
            assert(all(xy(:,1) <= ceil(size(img,2)/2)));
        end

        function test_error_1
            try
                cv.SimpleBlobDetector('foobar');
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function tno_persistence_test ( )

%*****************************************************************************80
%
%% TNO_PERSISTENCE_TEST tests the persistence of TNO data.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TNO_PERSISTENCE_TEST:\n' );
  fprintf ( 1, '  TNO is given a level L, and returns points and weights\n' );
  fprintf ( 1, '  of a Truncated Normal Odd (TNO) quadrature rule.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The first time TNO is called for a particular order N\n' );
  fprintf ( 1, '  it must compute the corresponding rule.  It prints a\n' );
  fprintf ( 1, '  message, and saves the data in persistent arrays.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  To verify this, compute a sequence of levels L\n' );
  fprintf ( 1, '  for which the values of N repeat.  The message should\n' );
  fprintf ( 1, '  only show up once for each value of N.\n' );

  fprintf ( 1, '\n' );
  for l = 1 : 6
    n = tno_order ( l );
    fprintf ( 1, '  L = %d, N = %d\n', l, n );
    [ x, w ] = tno ( l );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ycorr,yfit] = bf(y,varargin)
% Baseline Fit each column in "x".
% Syntax: [ycorr,yfit] = bf(y,pts,avgpts,method,confirm);
%   ycorr = bf(y);  ycorr = bf(y,method); ycorr = bf(y,avgpts);
%   ycorr = bf(y,pts);  ycorr = bf(y,pts,avgpts);
% A baseline fit is interpolated from selected points and then applied to the data.
%   "y" is a vector or array.
%       If an array, the baseline fit is performed for each column of data (dim 1).
%   Arguments following "y" may be in any order.
%   "pts" is vector specifying the indices of the points for a baseline fit.
%       If not specified then a plot is displayed and the user is instructed to
%       interactively select points for the baseline fit.
%       End points are always automatically included for interactive "pts" selection,
%       and do not need to explicitly selected. It is recommended that the end points
%       be included in any list of "pts".
%       It is not necessary to order or sort "pts" before use.
%   "avgpts" determines the width in points for the calculation of the mean y(x)
%       value, where x is a selected point in "pts". (Default = 3).
%       This can be helpful for noisy data.
%   "method" controls the algorithm applied for the baseline fit. The routine uses
%       Matlab's interp1 command. "method" must be one of the methods supported by
%       interp1. (Default is 'spline').
%   "confirm", if specified as the string 'confirm', will allow the user to see the
%       result and to confirm it is acceptable. If not the user can reslect "pts".
%   "ycorr" is the baseline corrected data in the same format as "y".
%   "yfit" is a vector or array with the interpolated baseline fit.
%
% Examples:
%   [y,yfit] = bf(y,'confirm','linear');
%       "y" will be plotted and the user is instructed to select points for the fit.
%       A baseline will be linearly interpolated from the selected points and will be
%       plotted together with "y". The user is prompted as to whether to redo the
%       baseline selection. Upon completion, the corrected data "y" and the fitted
%       baseline "yfit" are output.
%   ycorr = bf(y,5);
%      "y" is plotted and the user is instructed to select points for the fit.
%       The baseline fit is based on the mean value of "y" over 5 points centered on
%       the selected points. Cubic spline interpolation is used for the baseline fit.
%       The corrected data "ycorr" is output.
%   ycorr = bf(y,[5,10,15,30,35,40],'pchip');
%       Points with the specified indices are used to calculate a baseline fit using
%       the piecewise cubic Hermite interpolation method. No data is plotted.
%       The baseline fit is based on the mean value of "y" over 3 points centered on
%       the selected points. The corrected data "ycorr" is output.
%
% See Also:   interp1, spline, ginput

% Copyright 2009 Mirtech, Inc.
% Created by    Mirko Hrovat    08/01/2009  contact:mhrovat@email.com

def_method  = 'spline';
def_avgpts  = 3;

method = [];
avgpts = [];
pts    = [];
confirm = false;
for n = 2:nargin,
    f = varargin{n-1};
    if ischar(f),
        if strcmpi(f,'confirm'),
            confirm = true;
        else
            method = f;
        end
    elseif isnumeric(f) && numel(f) == 1,
        avgpts = f;
    elseif isnumeric(f) && numel(f) > 1,
        pts = f;
    elseif isempty(f),
        continue
    else
        error ('  Invalid input argument!')
    end
end
if isempty(method),     method = def_method;        end
if isempty(avgpts),     avgpts = def_avgpts;        end
dimy = size(y);
lst = dimy(1);
newdimy = [dimy(1),prod(dimy(2:end))];
y = reshape(y,newdimy);
x = 1:lst;
if isempty(pts),
    interactive = true;
else
    interactive = false;
end
if interactive || confirm,
    bffig = figure;
else
    bffig = 0;
end
ok = false;
while ~ok,
    if interactive,
        plot(x,real(y(:,1)))
        set(bffig,'Name','Baseline Fit - Select points')
        fprintf(['\n Now select baseline points to fit by positioning cursor,',...
            '\n   and selecting points with mouse button or key press.',...
            '\n Press Enter key when done.\n'])
        [a,b] = ginput;                                 %#ok
        pts = round(a.');
    end
    pts = sort(pts);
    pts(diff(pts)==0) = [];         % delete duplicate points
    if pts(1)~=1,       pts = [1,pts];          end     %#ok
    if pts(end)~=lst,   pts = [pts,lst];        end     %#ok
    npts = numel(pts);
    pss = zeros(npts,2);
    pss(:,1) = pts - floor(avgpts/2);
    pss(:,2) = pss(:,1) + avgpts;
    pss(pss < 1) = 1;
    pss(pss > lst) = lst;
    yavg = zeros([npts,newdimy(2)]);
    for n = 1:npts,
        yavg(n,:) = mean(y(pss(n,1):pss(n,2),:),1);
    end
    yfit = interp1(pts,yavg,x,method);
    if size(yfit,1) ==1,
        yfit = shiftdim(yfit,1);    % make yfit a column if it is a row vector
    end
    if confirm,
        interactive = true;
        figure(bffig)
        plot(x,real(y(:,1)),'b',x,real(yfit(:,1)),'r',pts,real(yavg(:,1)),'ob')
        set(bffig,'Name','Baseline Fit - Verify baseline')
        answer = input('  Do you to redo fit and reselect baseline points?[N] ','s');
        if isempty(answer),     answer = 'n';   end
        if strcmpi(answer,'y'),
            ok = false;
        else
            ok = true;
        end
    else
        ok = true;
    end
end
if any(findobj('Type','figure')==bffig),
    close(bffig),                   % close figure if it exists
end
ycorr = y - yfit;
ycorr = reshape(ycorr,dimy);
yfit = reshape(yfit,dimy);


text:
% select frames from the input stream
%
function grad = B_frame_select(input_layer, future_layers, curr_layer)
input = input_layer.a;
[D,T,N] = size(input);

future_grad = GetFutureGrad(future_layers, curr_layer);

words = ExtractWordsFromString_v2(curr_layer.frameSelect);
selectionType = words{1};
switch selectionType
    case 'last'
        if length(words)>1; nFrameSelect = str2num(words{2}); else nFrameSelect = 1; end
        if N>1; [mask, variableLength] = getValidFrameMask(input_layer); else variableLength = 0; end

        precision = class(gather(future_grad(1,1,1)));
        if strcmpi(class(future_grad), 'gpuArray')
            grad = gpuArray.zeros(D,T, N, precision);
        else
            grad = zeros(D,T, N, precision);
        end
        if variableLength
            last_idx = gather(GetLastValidFrameIndex(mask));
            for i=1:nSeg
                grad(:,(last_idx(i)-nFrameSelect+1):last_idx(i)) = future_grad;
            end
        else
            grad(:,max(1,end-nFrameSelect+1):end,:) = future_grad;
        end
end

end

author: singaxiong

repo: SignalGraph


text:
function atlas_out = split_atlas_into_contiguous_regions(atlas_obj)
% Divide regions with multiple contiguous blobs into separate labeled regions for each blob
%
% Take an atlas object whose labeled regions contain multiple contiguous blobs and divide each
% contiguous blob into a separate labeled region.
% - Note: This version eliminates probability maps - handling them not implemented yet
%
% July 2018, Tor Wager
%
% "split_atlas" method use cases:
% split_atlas_by_hemisphere: We have a defined set of bilateral regions that
% we want to "hard-split" into left and right. Multiple discontiguous regions
% with the same label will be kept together.
%
% split_atlas_into_contiguous_regions: We want to (1) keep contiguous blobs together
% that may cross the midline, and (2) separate contiguous blobs with the
% same label into separate labeled regions.

[n_regions, n_regions_with_data, missing_regions] = num_regions(atlas_obj);

atlas_out = atlas_obj;

atlas_out.labels = {};
atlas_out.label_descriptions = {};

atlas_out.probability_maps = []; % Eliminate probability maps - handling them not implemented yet

atlas_out.dat = zeros(size(atlas_out.dat)); % avoid 0 for resampling

fprintf('Splitting %d regions: 000', n_regions)

for i = 1:n_regions

    fprintf('%3.0f', i);

    subatlas = select_atlas_subset(atlas_obj, i);

    %     wh = logical(subatlas.dat);                      % which voxels in region
    %
    %     xyz = subatlas.volInfo.xyzlist(wh, :);           % x-coordinates in voxels
    %
    %     XYZmm = voxel2mm(xyz',subatlas.volInfo.mat);

    % parse into contiguous regions. split if there are multiple.
    r = atlas2region(subatlas);
    r = reparse_continguous(r);

    % clean up left-out voxels
    wh_omit = cat(1, r.numVox) < 3;
    r(wh_omit) = [];

    % newr is atlas object for this original region, with contig regions
    % divided

    new_subatlas = region2atlas(r, subatlas);

    new_subatlas.probability_maps = []; % Eliminate probability maps - handling them not implemented yet

    new_subatlas.labels = get_new_labels(r, subatlas.labels{1});

    atlas_out = merge_atlases(atlas_out, new_subatlas);

end % region

%atlas_out.dat(atlas_out.dat == 9999) = 0;

atlas_out.references = unique(atlas_out.references, 'rows');

disp('Done.');

end % function

function new_labels = get_new_labels(r, old_label)
% later: could merge all L, all R, all M

k = length(r);

for j = 1:k

    voxsign = sign(r(j).XYZmm(1, :)); % sign of each x coord. - = left
    modal_x = mode(voxsign);

    switch modal_x
        case -1
            labelstr = '_L';
        case 1
            labelstr = '_R';
        case 0
            labelstr = '_M';
    end

    % proportional asymmetry. 1 is vary lateralized. 0 is balanced
    % across L/R.  If symmetrical, then label _M
    prop_asym = abs(sum(voxsign == -1) - sum(voxsign == 1)) ./ length(voxsign);

    if prop_asym < .5, labelstr = '_M'; end

    new_labels{j} = [old_label labelstr];

end % j

% add left/right/mid

%     % split region into left, right, and midline within 10 mm
%     isleft = XYZmm(1, :) < 0;
%     isright = XYZmm(1, :) > 0;   % what to do with midline?
%     ismid = XYZmm == 0;

end % subfunction


text:
function [A,M] = assemblematrix3(node,elem,lumpflag)
%% ASSEMBLEMATRIX3 matrix for diffusion and reaction
%
% [A,M] = ASSEMBLEMATRIX3(node,elem) return the stiffness matrix and the
% mass matrix.
%
% [A,M] = ASSEMBLEMATRIX3(node,elem,1) return the stiffness matrix and the
% lumped mass matrix. Note that in the output M is a vector not a matrix. A
% sparse diagonal matrix using M as diaongal can be obtained by
% spdiags(M,0,N,N);
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

N = size(node,1);
A = sparse(N,N);
M = sparse(N,N);

%% Compute geometric quantities and gradient of local basis
[Dphi,volume] = gradbasis3(node,elem);

%% Assemble stiffness matrix
for i = 1:4
    for j = i:4
        Aij = dot(Dphi(:,:,i),Dphi(:,:,j),2).*volume;
        ii = double(elem(:,i));
        jj = double(elem(:,j));
        if (j==i)
            A = A + sparse(ii,jj,Aij,N,N);
        else
            A = A + sparse([ii;jj],[jj;ii],[Aij; Aij],N,N);
        end
        if (nargout > 1) && (~exist('lumpflag','var') || lumpflag == 0 )
            if (j==i)
                M = M + sparse(ii,jj,volume/10,N,N);
            else
                M = M + sparse([ii;jj],[jj;ii],[volume/20; volume/20],N,N);
            end
        end
    end
end

%% Assemble mass matrix by mass lumping
if exist('lumpflag','var') && lumpflag==1 && nargout > 1
    M = accumarray(elem(:),[volume;volume;volume;volume]/4,[N,1]);
end

repo: Aero-matlab


text:
function [ point_coord, face_order, face_point ] = ...
  tetrahedron_rhombic_shape_3d ( point_num, face_num, face_order_max )

%*****************************************************************************80
%
%% TETRAHEDRON_RHOMBIC_SHAPE_3D describes a rhombic tetrahedron in 3D.
%
%  Discussion:
%
%    Call TETRAHEDRON_RHOMBIC_SIZE_3D first, to get dimension information.
%
%    The tetrahedron is described using 10 nodes.  If we label the vertices
%    P0, P1, P2 and P3, then the extra nodes lie halfway between vertices,
%    and have the labels P01, P02, P03, P12, P13 and P23.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Anwei Liu, Barry Joe,
%    Quality Local Refinement of Tetrahedral Meshes Based
%    on 8-Subtetrahedron Subdivision,
%    Mathematics of Computation,
%    Volume 65, Number 215, July 1996, pages 1183-1200.
%
%  Parameters:
%
%    Input, integer POINT_NUM, the number of points in the shape.
%
%    Input, integer FACE_NUM, the number of faces in the shape.
%
%    Input, integer FACE_ORDER_MAX, the maximum number of vertices per face.
%
%    Output, real POINT_COORD(3,POINT_NUM), the vertices.
%
%    Output, integer FACE_ORDER(FACE_NUM), the number of vertices
%    for each face.
%
%    Output, integer FACE_POINT(FACE_ORDER_MAX,FACE_NUM); FACE_POINT(I,J)
%    contains the index of the I-th point in the J-th face.  The
%    points are listed in the counter clockwise direction defined
%    by the outward normal at the face.
%
  dim_num = 3;

  a =        1.0   / sqrt ( 3.0 );
  b = sqrt ( 2.0 ) / sqrt ( 3.0 );
  c = sqrt ( 3.0 ) /        6.0;
  d =        1.0   / sqrt ( 6.0 );
  z = 0.0;
%
%  Set the point coordinates.
%
  point_coord(1:dim_num,1)  = [ -b,  z,  z ]';
  point_coord(1:dim_num,2)  = [  z, -a,  z ]';
  point_coord(1:dim_num,3)  = [  z,  a,  z ]';
  point_coord(1:dim_num,4)  = [  z,  z,  b ]';
  point_coord(1:dim_num,5)  = [ -d, -c,  z ]';
  point_coord(1:dim_num,6)  = [ -d,  c,  z ]';
  point_coord(1:dim_num,7)  = [ -d,  z,  d ]';
  point_coord(1:dim_num,8)  = [  z,  z,  z ]';
  point_coord(1:dim_num,9)  = [  z, -c,  d ]';
  point_coord(1:dim_num,10) = [  z,  c,  d ]';
%
%  Set the face orders.
%
  face_order(1:face_num) = [ 6, 6, 6, 6 ];
%
%  Set faces.
%
  face_point(1:face_order_max,1:face_num) = [ ...
     1,  5,  2,  9,  4,  7; ...
     2,  8,  3, 10,  4,  9; ...
     3,  6,  1,  7,  4, 10; ...
     1,  6,  3,  8,  2,  5 ]';

  return
end


text:
function test_bug2170

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

% the 3rd dataset is different in the sense that it is a *.fif file accompanied
% by an *.eve file. Those files are also used in the Tristan babysquid74 system.

filename1 = dccnpath('/home/common/matlab/fieldtrip/data/test/original/meg/neuromag306/raw.fif');
filename2 = dccnpath('/home/common/matlab/fieldtrip/data/test/original/meg/neuromag306/run_01_raw.fif');
filename3 = dccnpath('/home/common/matlab/fieldtrip/data/test/bug2170/an05a4_ss.fif');

assert(ft_filetype(filename1, 'neuromag_fif'));
assert(ft_filetype(filename2, 'neuromag_fif'));
assert(ft_filetype(filename3, 'neuromag_fif'));

event1 = ft_read_event(filename1);
event2 = ft_read_event(filename2);
event3 = ft_read_event(filename3);

% Note: on 7 June 2017 Robert made a change to ft_read_trigger, causing
% chantype='other trigger' also to be returned. In the 1st and 2nd dataset that
% corresponds to STI201 and STI301.

type = unique({event1.type});
for i=1:numel(type)
  figure
  tmp = event1(strcmp({event1.type}, type{i}));
  plot([tmp.sample], [tmp.value], '.')
  title(type{i})
end

assert(length(unique({event1.type}))==7+3+1); % 'STI001' 'STI002' 'STI003' 'STI004' 'STI009' 'STI010' 'STI011' 'STI101' 'STI201' 'STI301' and an additional 'Trigger'
assert(length(unique({event2.type}))==2);     % 'STI101' and 'STI201'
assert(length(unique({event3.type}))==2);     % 'STI101' and 'trigger'


text:
function v = volume(sR,S2G)
% volume of a spherical region
%
%
% TODO: make this better!!
% formula: v = S - (n-2)*pi
% where S is the sum of the interior angles and n is the number of vertices

if nargin == 1, S2G = equispacedS2Grid(0.5*degree); end

v = nnz(sR.checkInside(S2G))/length(S2G);

author: mtex-toolbox

repo: mtex

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


text:
function y =Fault_decision_b2m2_t2(u)
global st17

% Using model defined by the structure st17 to make decision about type2 fault in the pitch beta2m2

res=u(1);u=u(2:end);
if res>0.5
sker=st17.x2sup+(abs(u))'*abs(u)*ones(st17.Nlsup,1)-2*st17.xsup*abs(u);
y=(st17.w)'*exp(-sker./(2*(st17.sigma).^2))+st17.b;
else
    y=0;
end
if y>=0.5
    y=1;
else
    y=0;
end

end


text:
function sgmga_size_table

%*****************************************************************************80
%
%% MAIN is the main program for SGMGA_SIZE_TABLE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2010
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../sandia_rules' );

  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SGMGA_SIZE_TABLE:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Make tables of point counts.\n' );
  fprintf ( 1, '  Measure the CPU time for each table.\n' );
%
%  Clenshaw-Curtis Grid (1), slow exponential growth (4).
%
  rule_1d = 1;
  growth_1d = 4;
  np_1d = 0;
  p_1d = [];
  dim_min = 1;
  dim_max = 5;
  level_max_min = 0;
  level_max_max = 5;
  tic ( );
  sgmga_size_tabulate ( rule_1d, growth_1d, np_1d, p_1d, dim_min, dim_max, ...
    level_max_min, level_max_max );
  ctime = toc ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  CPU Time = %f\n', ctime );
%
%  Clenshaw Curtis Grid (1), exponential growth (6).
%
  rule_1d = 1;
  growth_1d = 6;
  np_1d = 0;
  p_1d = [];
  dim_min = 1;
  dim_max = 5;
  level_max_min = 0;
  level_max_max = 5;
  tic ( );
  sgmga_size_tabulate ( rule_1d, growth_1d, np_1d, p_1d, dim_min, dim_max, ...
    level_max_min, level_max_max );
  ctime = toc ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  CPU Time = %f\n', ctime );
%
%  Clenshaw Curtis Grid (1), exponential growth (6).
%
  if ( 0 )
    rule_1d = 1;
    growth_1d = 6;
    np_1d = 0;
    p_1d = [];
    dim_min = 6;
    dim_max = 10;
    level_max_min = 0;
    level_max_max = 5;
    sgmga_size_tabulate ( rule_1d, growth_1d, np_1d, p_1d, dim_min, dim_max, ...
      level_max_min, level_max_max );
    ctime = toc ( );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  CPU Time = %f\n', ctime );
  end
%
%  Clenshaw Curtis Grid (1), exponential growth (6).
%
  rule_1d = 1;
  growth_1d = 6;
  np_1d = 0;
  p_1d = [];
  dim_min = 100;
  dim_max = 100;
  level_max_min = 0;
  level_max_max = 2;
  tic ( );
  sgmga_size_tabulate ( rule_1d, growth_1d, np_1d, p_1d, dim_min, dim_max, ...
    level_max_min, level_max_max );
  ctime = toc ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  CPU Time = %f\n', ctime );
%
%  Gauss Patterson Grid (3), slow exponential growth (4).
%
  rule_1d = 3;
  growth_1d = 4;
  np_1d = 0;
  p_1d = [];
  dim_min = 1;
  dim_max = 5;
  level_max_min = 0;
  level_max_max = 5;
  tic ( );
  sgmga_size_tabulate ( rule_1d, growth_1d, np_1d, p_1d, dim_min, dim_max, ...
    level_max_min, level_max_max );
  ctime = toc ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  CPU Time = %f\n', ctime );
%
%  Gauss Patterson Grid (3), moderate exponential growth (5).
%
  rule_1d = 3;
  growth_1d = 5;
  np_1d = 0;
  p_1d = [];
  dim_min = 1;
  dim_max = 5;
  level_max_min = 0;
  level_max_max = 5;
  tic ( );
  sgmga_size_tabulate ( rule_1d, growth_1d, np_1d, p_1d, dim_min, dim_max, ...
    level_max_min, level_max_max );
  ctime = toc ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  CPU Time = %f\n', ctime );
%
%  Gauss Patterson Grid (3), exponential growth (6).
%
  rule_1d = 3;
  growth_1d = 6;
  np_1d = 0;
  p_1d = [];
  dim_min = 1;
  dim_max = 5;
  level_max_min = 0;
  level_max_max = 5;
  tic ( );
  sgmga_size_tabulate ( rule_1d, growth_1d, np_1d, p_1d, dim_min, dim_max, ...
    level_max_min, level_max_max );
  ctime = toc ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  CPU Time = %f\n', ctime );
%
%  Gauss Legendre Grid (4), slow linear odd growth (2).
%
  rule_1d = 4;
  growth_1d = 2;
  np_1d = 0;
  p_1d = [];
  dim_min = 1;
  dim_max = 5;
  level_max_min = 0;
  level_max_max = 5;
  tic ( );
  sgmga_size_tabulate ( rule_1d, growth_1d, np_1d, p_1d, dim_min, dim_max, ...
    level_max_min, level_max_max );
  ctime = toc ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  CPU Time = %f\n', ctime );
%
%  Gauss Legendre Grid (4), moderate linear growth (3).
%
  rule_1d = 4;
  growth_1d = 3;
  np_1d = 0;
  p_1d = [];
  dim_min = 1;
  dim_max = 5;
  level_max_min = 0;
  level_max_max = 5;
  tic ( );
  sgmga_size_tabulate ( rule_1d, growth_1d, np_1d, p_1d, dim_min, dim_max, ...
    level_max_min, level_max_max );
  ctime = toc ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  CPU Time = %f\n', ctime );
%
%  Gauss Laguerre Grid (7), moderate linear growth (3).
%
  rule_1d = 7;
  growth_1d = 3;
  np_1d = 0;
  p_1d = [];
  dim_min = 1;
  dim_max = 5;
  level_max_min = 0;
  level_max_max = 5;
  tic ( );
  sgmga_size_tabulate ( rule_1d, growth_1d, np_1d, p_1d, dim_min, dim_max, ...
    level_max_min, level_max_max );
  ctime = toc ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  CPU Time = %f\n', ctime );
%
%  Hermite Genz Keister (10), slow exponential growth (4).
%
  rule_1d = 10;
  growth_1d = 4;
  np_1d = 0;
  p_1d = [];
  dim_min = 1;
  dim_max = 5;
  level_max_min = 0;
  level_max_max = 5;
  tic ( );
  sgmga_size_tabulate ( rule_1d, growth_1d, np_1d, p_1d, dim_min, dim_max, ...
    level_max_min, level_max_max );
  ctime = toc ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  CPU Time = %f\n', ctime );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SGMGA_SIZE_TABLE:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../sandia_rules' );

  return
end


text:
function detail_layer_out = StevensonDetailEnhancement(detail_layer, F_L)
%
%       detail_layer_out = StevensonDetailEnhancement(detail_layer, F_L)
%
%       This function adjusts the detail layer for taking into account the
%       Stevenson effect.
%
%       input:
%           -detail_layer: the detail layer of an image
%           -F_L: a luminance dependent factor
%
%       output:
%           -detail_layer_out: processed detail layer
%
%     Copyright (C) 2014-15  Francesco Banterle
%
%     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/>.
%

exponent = (F_L + 0.8).^0.25;

[~, ~, col] = size(detail_layer);

detail_layer_out = zeros(size(detail_layer));

for i=1:col
    detail_layer_out(:,:,i) = detail_layer(:,:,i).^exponent;
end

end

author: banterle


text:
function feat = ExtractSpatialCovFeat(X, nCh, context_size, shift, binStep, rowOnly, nChLogMag)
if nargin<5;    binStep = 1; end
if nargin<6;    rowOnly = 0; end
if nargin<7;    nChLogMag = 7; end

[DC, T] = size(X);
nFFT = DC/nCh;
X2 = reshape(X, nFFT, nCh, T);

if 0    % apply Mel filterbank
    nFbank = 40;
    MelWindow = mel_window_FE(nFbank, nFFT-1, 16000)';
    MelWindow(:,end+1) = 0;
    X_tmp = reshape(X2, nFFT, nCh*T);
    X_fbank = MelWindow * X_tmp;
    data = reshape(X_fbank, nFbank, nCh, T);
else
    data = X2;
end

nBin = size(data,1);
bin_selector = binStep+1:binStep:nBin-1;    % we don't want the first and last bin as they are always real values
nSelected = length(bin_selector);

% compute the spatial covariance using context
if strcmpi(class(data), 'gpuArray')
    useGPU = 1;
else
    useGPU = 0;
end
spatialCov = ComplexSpectrum2SpatialCov(data(bin_selector,:,:), context_size, shift, useGPU);

% normalize the covariance matrix by their diagonal elements. This step
% removes the effect of spectral power and only retains the phase
% information (ideally).
nFrCov = size(spatialCov,4);
if 0
    for i=1:nSelected
        for j=1:nFrCov
            meanPower(i,j) = mean(diag(spatialCov(:,:,i,j)));
        end
    end
elseif 0    % still slow
    spatialCovCell = num2cell(gather(spatialCov), [1 2]);
    meanPowerCell = cellfun(@(x) mean(diag(x)), spatialCovCell, 'UniformOutput', 0);
    meanPower = squeeze(cell2mat(meanPowerCell));
else    % use indexing
    dimSelectMask = zeros(nCh,nCh,nSelected);
    for i=1:nSelected
        dimSelectMask(:,:,i) = eye(nCh);
    end
    dimSelectIdx = find(dimSelectMask(:) == 1);
    spatialCov2 = reshape(spatialCov, nCh^2*nSelected,nFrCov);
    diag_part = spatialCov2(dimSelectIdx,:);
    meanPower = squeeze(mean(reshape(diag_part, nCh, nSelected, nFrCov)));
end
spatialCovNorm = bsxfun(@times, permute(spatialCov, [3 4 1 2]), 1./meanPower);
spatialCovNorm = permute(spatialCovNorm, [3 4 1 2]);

% get the upper triangle off-diagonal elements which are complex-valued
dimSelectMask = zeros(nCh,nCh,nSelected);
if rowOnly
    dimSelectMask(1,2:end,:) = 1;
else
    for i=1:nCh
        dimSelectMask(i,i+1:end,:) = 1;
    end
end
dimSelectIdx = find(dimSelectMask(:) == 1);

spatialCovNorm2 = reshape(spatialCovNorm, nCh^2*nSelected,nFrCov);
real_part = real(spatialCovNorm2(dimSelectIdx,:));
imag_part = imag(spatialCovNorm2(dimSelectIdx,:));

% get the diagoanl elements which are real values
dimSelectMask = zeros(nCh,nCh,nSelected);
for i=1:nSelected
    for j=1:nChLogMag
        dimSelectMask(j,j,i) = 1;
    end
end
dimSelectIdx = find(dimSelectMask(:) == 1);
spatialCov2 = reshape(spatialCov, nCh^2*nSelected,nFrCov);
diag_part = spatialCov2(dimSelectIdx,:);
diag_part = log(max(eps,abs(diag_part)));

% get the final feature vector
feat = [real_part; imag_part; diag_part];

end


text:
%InterfaceUtil
% Checks if a class (or class instance) properly implements a specificed
% interface.
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 9/12/2010
classdef InterfaceUtil
    methods (Static = true)
        %Throws an exception is the class/class instance doesn't properly
        %implement the specified interface.
        %
        %Inputs:
        %- instanceOrClassName is an instance of a class, or a char equal to the
        %  name of a class
        %- metaclass is a struct of duck type meta.class with properties
        %  - Properties: cell array of objects of duck type meta.property. Each
        %    element of Properties must have the field Name,
        %    corresponding to the name of the property whose attributes should
        %    be validated. The other field names of Properties must equal
        %    names of property attributes, and the field values represent
        %    acceptable values of the corresponding property attribute. A
        %    Properties field value of null means that any value of the
        %    property attribute is acceptable.
        %  - Methods: cell array of objects of duck type meta.method. Each
        %    element of Methods must have the field Name, corresponding
        %    to the name of the method whose attributes should be validated. The
        %    other field names of Methods must equal names of method
        %    attributes, and the field values represent acceptable values of the
        %    corresponding method attribute. A Methods field value of
        %    null means that any value of the method attribute is acceptable.
        function assertInterface(instanceOrClassName, instanceOrMetaclass)
            import edu.stanford.covert.util.InterfaceUtil;

            if isa(instanceOrMetaclass, 'meta.class') || isstruct(instanceOrMetaclass)
                metadata = instanceOrMetaclass;
            elseif ischar(instanceOrMetaclass)
                metadata = meta.class.fromName(instanceOrMetaclass);
            else
                metadata = metaclass(instanceOrMetaclass);
            end

            InterfaceUtil.assertPropertyInterface(instanceOrClassName, metadata.Properties);
            InterfaceUtil.assertMethodInterface(instanceOrClassName, metadata.Methods);
        end
    end

    methods (Static = true, Access = protected)
        function assertPropertyInterface(instanceOrClassName, interface)
            import edu.stanford.covert.util.InterfaceUtil;

            if isempty(interface)
                return;
            end

            %validate input
            if ~isvector(interface) || ~iscell(interface) || ~all(cellfun(@InterfaceUtil.hasDuckTypeMetaProperty, interface))
                throw(MException('Interface:invalidInput', 'interface must be a cell array of objects of duck type meta.property'))
            end

            %get metadata
            if ischar(instanceOrClassName)
                metadata = meta.class.fromName(instanceOrClassName);
            else
                metadata = metaclass(instanceOrClassName);
            end

            %check property interface
            propertyNames = cellfun(@(property) property.Name, metadata.Properties, 'UniformOutput', false);
            for i = 1:numel(interface)
                iProperty = find(strcmp(propertyNames, interface{i}.Name), 1, 'first');
                interfaceFieldNames = setdiff(fieldnames(interface{i}), {'Name', 'Description', 'DetailedDescription', 'DefiningClass'});

                if isempty(iProperty)
                    throw(MException('Interface:incompleteInterface', ...
                        'The class %s must implement a %s property', ...
                        metadata.Name, interface{i}.Name));
                end

                property = metadata.Properties{iProperty};
                for j = 1:numel(interfaceFieldNames)
                    if isempty(interface{i}.(interfaceFieldNames{j}))
                        continue;
                    end

                    interfaceValue = interface{i}.(interfaceFieldNames{j});
                    interfaceValueStr = interface{i}.(interfaceFieldNames{j});
                    if ischar(interfaceValue)
                        interfaceValueStr = {interfaceValue};
                        interfaceValue = {interfaceValue};
                    elseif ~iscell(interfaceValue)
                        interfaceValueStr = cellfun(@num2str, num2cell(interfaceValue), 'UniformOutput', false);
                    end

                    propertyValue = property.(interfaceFieldNames{j});
                    if ischar(propertyValue)
                        %GetAccess, SetAccess
                        propertyValue = {propertyValue};
                    elseif isa(propertyValue, 'function_handle')
                        %GetMethod, SetMethod
                        propertyValue = ~isempty(propertyValue);
                    end

                    if ~ismember(propertyValue, interfaceValue)
                        throw(MException('Interface:incompatibleInterface', ...
                            'The class %s %s property must have attribute %s with value {%s}', ...
                            metadata.Name, interface{i}.Name, interfaceFieldNames{j}, ['''' strjoin(''', ''',interfaceValueStr{:}) '''']));
                    end
                end
            end
        end

        function assertMethodInterface(instanceOrClassName, interface)
            import edu.stanford.covert.util.InterfaceUtil;

            if isempty(interface)
                return;
            end

            %validate input
            if ~isvector(interface) || ~iscell(interface) || ~all(cellfun(@InterfaceUtil.hasDuckTypeMetaMethod, interface))
                throw(MException('Interface:invalidInput', 'interface must be a cell array of objects of duck type meta.method'))
            end

            %get metadata
            if ischar(instanceOrClassName)
                metadata = meta.class.fromName(instanceOrClassName);
            else
                metadata = metaclass(instanceOrClassName);
            end

            %check method interface
            methodNames = cellfun(@(method) method.Name, metadata.Methods, 'UniformOutput', false);
            for i = 1:numel(interface)
                iMethod = find(strcmp(methodNames, interface{i}.Name), 1, 'first');
                interfaceFieldNames = setdiff(fieldnames(interface{i}), {'Name', 'Description', 'DetailedDescription', 'DefiningClass'});

                if isempty(iMethod)
                    throw(MException('Interface:incompleteInterface', ...
                        'The class %s must implement a %s method', ...
                        metadata.Name, interface{i}.Name));
                end

                method = metadata.Methods{iMethod};
                for j = 1:numel(interfaceFieldNames)
                    if isempty(interface{i}.(interfaceFieldNames{j}))
                        continue;
                    end

                    interfaceValue = interface{i}.(interfaceFieldNames{j});
                    interfaceValueStr = interface{i}.(interfaceFieldNames{j});
                    if ischar(interfaceValue)
                        interfaceValueStr = {interfaceValue};
                        interfaceValue = {interfaceValue};
                    elseif ~iscell(interfaceValue)
                        interfaceValueStr = cellfun(@num2str, num2cell(interfaceValue), 'UniformOutput', false);
                    end

                    methodValue = method.(interfaceFieldNames{j});
                    if ischar(methodValue)
                        %Access
                        methodValue = {methodValue};

                        if ~ismember(methodValue, interfaceValue)
                            throw(MException('Interface:incompatibleInterface', ...
                                'The class %s %s method must have attribute %s with value {%s}', ...
                                metadata.Name, interface{i}.Name, interfaceFieldNames{j}, ['''' strjoin(''', ''',interfaceValueStr{:}) '''']));
                        end
                    elseif iscell(methodValue)
                        %InputNames, OutputNames
                        methodValue = numel(methodValue);
                        if iscell(interfaceValue)
                            interfaceValue = numel(interfaceValue);
                        elseif numel(interfaceValue)~=1
                            throw(MException('Interface:invalidInput', 'interface must specify number or names of %s',interfaceFieldNames{j}));
                        end

                        if methodValue < interfaceValue
                            throw(MException('Interface:incompatibleInterface', ...
                                'The class %s %s method must have at at least %d %s', ...
                                metadata.Name, interface{i}.Name, interfaceValue, interfaceFieldNames{j}));
                        end
                    else
                        if ~ismember(methodValue, interfaceValue)
                            throw(MException('Interface:incompatibleInterface', ...
                                'The class %s %s method must have attribute %s with value {%s}', ...
                                metadata.Name, interface{i}.Name, interfaceFieldNames{j}, ['''' strjoin(''', ''',interfaceValueStr{:}) '''']));
                        end
                    end
                end
            end
        end

        function tf = hasDuckTypeMetaProperty(prop)
            tf = false;

            if isa(prop, 'meta.property')
                tf = true;
                return;
            end

            metaPropertyFieldNames = {'Name'; 'Description'; 'DetailedDescription';
                'GetAccess'; 'SetAccess'; 'Dependent'; 'Constant'; 'Abstract'; 'Transient';
                'Hidden'; 'GetObservable'; 'SetObservable'; 'AbortSet'; 'GetMethod';
                'SetMethod'; 'HasDefault'; 'DefaultValue'; 'DefiningClass'};

            if isstruct(prop) && ismember('Name',fieldnames(prop)) && all(ismember(fieldnames(prop), metaPropertyFieldNames))
                tf = true;
                return
            end
        end

        function tf = hasDuckTypeMetaMethod(method)
            tf = false;

            if isa(method, 'meta.method')
                tf = true;
                return;
            end

            metaMethodFieldNames = {'Name'; 'Description'; 'DetailedDescription';
                'Access'; 'Static'; 'Abstract'; 'Sealed'; 'Hidden'; 'InputNames';
                'OutputNames'; 'DefiningClass'};

            if isstruct(method) && ismember('Name',fieldnames(method)) && all(ismember(fieldnames(method), metaMethodFieldNames))
                tf = true;
                return
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [gray_enhanced]=gray_level_images(current_image);
 %% trying dc coefficients
[mrows,ncolumns]=size(current_image);
 current_image=mat2gray(current_image);
        H = fspecial('disk',10);
        gray_filtering=current_image;
        blurred_gray = imfilter(gray_filtering,H,'replicate');
        k=2;
    for i=1:1:mrows
        for j=1:1:ncolumns
gray_enhanced(i,j)=mat2gray((gray_filtering(i,j)+(k.*(gray_filtering(i,j)-blurred_gray(i,j)))));
        end
    end
%figure(10001);
imshow(gray_enhanced);


text:
function [coord_norm] = warp_fsaverage_sym(cfg, elec)

% WARP_FSAVERAGE_SYM maps left or right hemisphere electrodes onto
% FreeSurfer's fsaverage_sym's left hemisphere. To perform this mapping,
% you first need to have processed the subject's MRI with FreeSurfer's
% recon-all functionality and additionaly have registered the subject's resulting
% surfaces to freesurfer fsaverage_sym template using surfreg as described
% in section 1.2 of https://surfer.nmr.mgh.harvard.edu/fswiki/Xhemi
%
% The configuration must contain the following options
%   cfg.headshape      = string, filename containing subject headshape
%                      (e.g. <path to freesurfer/surf/lh.pial>)
%   cfg.fshome         = string, path to freesurfer
%
% See also FT_ELECTRODEREALIGN, WARP_FSAVERAGE

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

subj_pial = ft_read_headshape(cfg.headshape);
[PATHSTR, NAME] = fileparts(cfg.headshape); % lh or rh
if strcmp(NAME, 'lh')
  subj_reg = ft_read_headshape([PATHSTR filesep 'lh.fsaverage_sym.sphere.reg']);
elseif strcmp(NAME, 'rh')
  subj_reg = ft_read_headshape([PATHSTR(1:strfind(PATHSTR, [filesep 'surf'])-1) filesep 'xhemi' filesep 'surf' filesep 'lh.fsaverage_sym.sphere.reg']);
end
if ~isfolder([cfg.fshome filesep 'subjects' filesep 'fsaverage_sym']) || ~isfolder([PATHSTR(1:strfind(PATHSTR, [filesep 'surf'])-1) filesep 'xhemi'])
  ft_error(['fsaverage_sym and/or xhemi folders cannot be found'])
end
fsavg_pial = ft_read_headshape([cfg.fshome filesep 'subjects' filesep 'fsaverage_sym' filesep 'surf' filesep 'lh.pial']);
fsavg_reg = ft_read_headshape([cfg.fshome filesep 'subjects' filesep 'fsaverage_sym' filesep 'surf' filesep 'lh.sphere.reg']); % always map onto the left hemi

for e = 1:numel(elec.label)
  % subject space (3D surface): electrode pos -> vertex index
  dist = sqrt(sum(((subj_pial.pos - repmat(elec.elecpos(e,:), size(subj_pial.pos,1), 1)).^2),2));
  [dum, minidx] = min(dist);

  % intersubject space (3D sphere): vertex index -> vertex pos -> template vertex index
  dist2 = sqrt(sum(((fsavg_reg.pos - repmat(subj_reg.pos(minidx,:), size(fsavg_reg.pos,1), 1)).^2),2));
  [dum, minidx2] = min(dist2);
  clear minidx

  % template space (3D surface): template vertex index -> template electrode pos
  coord_norm(e,:) = fsavg_pial.pos(minidx2,:);
  clear minidx2
end
clear subj_pial subj_reg fsavg_pial fsavg_reg


text:
function [ c0, esterr ] = padua2 ( deg, degmax, npd, wpd, fpd )

%*****************************************************************************80
%
%% PADUA2 computes the Padua interpolation coefficient matrix.
%
%  Discussion:
%
%    This subroutine computes the coefficient matrix C0, in the
%    orthonormal Chebyshev basis T_j(x)T_{k-j}(y), 0 <= j <= k <= DEG,
%    T_0(x)=1, T_j(x) = sqrt(2) * cos(j * acos(x)), of the
%    interpolation polynomial of degree DEG of the function values FPD
%    at the set of NPD Padua points (PD1,PD2) in the square [-1,1]^2.
%
%    The interpolant may be evaluated at an arbitrary point by the
%    function PD2VAL. PD1, PD2 and WPD are the Padua points and weights
%    computed by PDPTS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2014
%
%  Author:
%
%    Original FORTRAN77 version by Marco Caliari, Stefano De Marchi,
%    Marco Vianello.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Marco Caliari, Stefano de Marchi, Marco Vianello,
%    Algorithm 886:
%    Padua2D: Lagrange Interpolation at Padua Points on Bivariate Domains,
%    ACM Transactions on Mathematical Software,
%    Volume 35, Number 3, October 2008, Article 21, 11 pages.
%
%  Parameters:
%
%    Input, integer DEG, the degree of approximation.
%
%    Input, integer DEGMAX, the maximum degree allowed.
%
%    Input, integer NPD, the number of Padua points.
%
%    Input, real WPD(NPD), the weights.
%
%    Input, real FPD(NPD), the value at the Padua points
%    of the function to be interpolated.
%
%    Output, real C0(0+1:DEG+1,0+1:DEG+1), the coefficient matrix.
%
%    Output, real ESTERR, the estimated error.
%
  BASE = 1;
%
%  Build the matrix P_2 and store it in RAUX2.
%
  raux2 = zeros ( deg + 1, deg + 2 );

  for i = 0 : deg + 1
    angle = i * pi / ( deg + 1 );
    pt = - cos ( angle );
    raux2(0+BASE:deg+BASE,i+1) = cheb ( deg, pt );
  end
%
%  Build the matrix G(f) and store it in C0.
%
  k = 0;
  for j = 0 : deg + 1
    for i = 0 : deg
      if ( mod ( i + j, 2 ) == 0 )
        k = k + 1;
        c0(i+BASE,j+BASE) = fpd(k) * wpd(k);
      else
        c0(i+BASE,j+BASE) = 0.0;
      end
    end
  end
%
%  Compute the matrix-matrix product G(f)*P_2' and store it in RAUX1.
%
  raux1 = zeros(degmax+1,deg+2);
  raux1 = dgemm ( 'n', 't', deg + 1, deg + 1, deg + 2, 1.0, ...
    c0, degmax + 2, raux2, degmax + 1, 0.0, raux1, degmax + 1 );
%
%  Build the matrix P_1 and store it in RAUX2.
%
  for i = 0 : deg
    angle = i * pi / deg;
    pt = - cos ( angle );
    raux2(0+BASE:deg+BASE,i+1) = cheb ( deg, pt );
  end
%
%  Compute the matrix-matrix product C(f) = P_1 * ( G(f) * P_2' )
%  and store it in C0.
%
  c0 = dgemm ( 'n', 'n', deg + 1, deg + 1, deg + 1, 1.0, ...
    raux2, degmax + 1, raux1, degmax + 1, 0.0, c0, degmax + 2 );

  c0(deg+BASE,0+BASE) = c0(deg+BASE,0+BASE) / 2.0;
%
%  Estimate the error.
%
  esterr = 0.0;
  for j = 0 : 2
    for i = 0 : deg - j
      esterr = esterr + abs ( c0(i+BASE,deg-i-j+BASE) );
    end
  end
  esterr = 2.0 * esterr;

  return
end


text:
function all_parts = add_sp_parts(sp, sp_list)
% Sp are made of parts of original superpixelation. This function gives the
% list of combined parts of a collection of sp.

all_parts = [];

for i = sp_list
    all_parts = [all_parts, sp{i}.parts];
end

author: Cloud-CV

repo: object-proposals


function [  ] = test2( vis, compressIndex, totalNum )
%TEST2 Summary of this function goes here
%   Detailed explanation goes here

    row = zeros(1, totalNum);
    for i=1:length(vis)
        row(1, vis(i)+1) = 1;
    end

    column = zeros(totalNum, 1);
    for i=1:length(compressIndex)
        column(compressIndex(i),1) = 1;
    end

    b = row*column
end


text:
function box_behnken_test01 ( )

%*****************************************************************************80
%
%% BOX_BEHNKEN_TEST01 tests BOX_BEHNKEN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 October 2008
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;

  range = [ ...
    0.0, 10.0,  5.0; ...
    1.0, 11.0, 15.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BOX_BEHNKEN_TEST01\n' );
  fprintf ( 1, '  BOX_BEHNKEN computes a Box-Behnken dataset.\n' );

  r8mat_transpose_print ( dim_num, 2, range, '  The ranges:' );

  x_num = box_behnken_size ( dim_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  For dimension DIM_NUM = %d\n', dim_num );
  fprintf ( 1, ' the Box-Behnken design is of size %d\n', x_num );

  x = box_behnken ( dim_num, x_num, range );

  r8mat_transpose_print ( dim_num, x_num, x, '  The Box-Behnken design:' );

  return
end


text:
function dt6_new = dtiResliceTensorAffine(dt6, xform, mmPerVoxIn, boundingBox, mmPerVoxOut)
% dt6_new = dtiResliceTensorAffine(dt6, xform, mmPerVoxIn, [boundingBox], [mmPerVoxOut])
%
% Applies an affine transform to a tensor volume in dt6 format.
%
% REQUIRES:
% spm99 or spm2 in the path.
%
% HISTORY:
% 2004.02.01 RFD & ASH Wrote it.
% 2004.08.06 RFD & ASH- we now use the Finite Strain method, which means
% that we no longer ignore shears.
% 2005.03.25 RFD- fixed the bug that was introduced with the last edit
% (2004.08.06). The tensor rotation was the inverse of what it should have
% been! I also think that applying the tensor rotation after the
% interpolation was introducing some subtle artifacts, so I switched back
% to rotating the tensors before interpolating (the order really shouldn't
% matter, so we should look into this more to understand why it does.)
%
% example:
% sn3dParams = load('dti_analyses/may021126_B0_sn3d.mat');
% xformToB0 = computeCannonicalXformFromIfile('dti/B0.001');
% xformToTensor = inv(xformToB0);
% [eigVec,eigVal,mm] = dtiLoadTensor('/snarp/u1/dti/adultData/may021126/dti/Vectors.float');
% dt6 = dtiRebuildTensor(eigVec, eigVal);
% dt6_new = dtiResliceTensor(dt6, sn3dParams, mmPerVoxOut, xformToTensor);
% [newVec, newVal] = dtiSplitTensor(dt6_new);
% fa = dtiComputeFA(newVal);
% figure; imagesc(fa(:,:,23)); axis equal; colormap gray; axis xy;

if(~exist('mmPerVoxOut') | isempty(mmPerVoxOut))
    mmPerVoxOut = [1 1 1];
end

if(~exist('boundingBox') | isempty(boundingBox))
    % create a default bounding box for the resliced data.
    % The following is the spm default bounding box. Note that the box is
    % defined in Talairach space (units = mm).
    boundingBox = [-78 -120 -60;
                    78  80   85];
end

x = (boundingBox(1,1):mmPerVoxOut(1):boundingBox(2,1));
y = (boundingBox(1,2):mmPerVoxOut(2):boundingBox(2,2));
z = (boundingBox(1,3):mmPerVoxOut(3):boundingBox(2,3));

[X,Y,Z] = meshgrid(x, y, z);
clear x y z;

talCoords = [X(:) Y(:) Z(:)];
newSize = size(X);
clear X Y Z;

% Does it matter if tensor rotation is done before or after rotation?
% ROTATE THE TENSORS
% NOTE: we want to apply inv(xform) to the tensors, since xform maps from
% the new space to the old (the correct mapping for the interpolation,
% since we interpolate by creating a grid in the new space and fill it by
% pulling data from the old space.)
rigidXform = dtiFiniteStrainDecompose(inv(xform));
% Apply just the rigid rotation component of the xform to the tensors
dt6 = dtiXformTensors(dt6, rigidXform);

imgCoords = mrAnatXformCoords(xform, talCoords)';
%imgCoords = imgCoords./repmat(mmPerVoxIn',1,size(imgCoords,2));

% the imgCoords are specified in mm, as an offset from the origin. So, we
% convert the 1-indexed imgCoords from mrAnatXformCoords to zero-indexed and
% scale to real mm.
imgCoords = imgCoords - 1;
imgCoords(1,:) = imgCoords(1,:).*mmPerVoxIn(1);
imgCoords(2,:) = imgCoords(2,:).*mmPerVoxIn(2);
imgCoords(3,:) = imgCoords(3,:).*mmPerVoxIn(3);

% We need to mask away the crazy values outside the head.
% Also, the Pajevic algorithm will fold the data over to deal with the
% edges. That might be nice for computing edge values, but it makes for
% ugly images and crazy fiber tacings. The following will build a mask of
% the good image values (ie. values inside the head), interpolate it to the
% new volume size, and then invert it so it marks bad values, those
% outside the head.
dt6_mask = sum(abs(dt6),4)>0.001;
dt6_mask = dtiResliceVolume(double(dt6_mask), xform, boundingBox, mmPerVoxOut);
dt6_mask = dt6_mask>0.5;
dt6_mask = dtiCleanImageMask(dt6_mask);
%figure; imagesc(makeMontage(dt6_mask)); axis equal off; colormap gray;
dt6_mask = ~dt6_mask;

% *** HACK so that matlab can find the shared library
old = pwd;
cd(fileparts(which('dtiTensorInterp_Pajevic')));
dt6_new = dtiTensorInterp_Pajevic(dt6, [imgCoords(1,:);imgCoords(2,:);imgCoords(3,:)]', mmPerVoxIn, 1, mmPerVoxIn./2);
cd(old);
dt6_new(repmat(dt6_mask(:),1,6)) = 0;
dt6_new = reshape(dt6_new, [newSize, 6]);

return;

%
% DEBUGGING:
[vec,val] = dtiSplitTensor(dt6_new);
val(val<100) = 0;
fa_interp = dtiComputeFA(val);
figure; imagesc(fa_interp(:,:,20)); axis equal xy tight off; colorbar;
b0 = loadAnalyze('dti_analyses/may021126_nB0.img');
b0 = permute(b0, [2 1 3]);
figure; imagesc(b0(:,:,20)); axis equal xy tight off; colorbar;
mask = b0>500;
fa_interp(~mask) = 0;
figure; imagesc(fa_interp(:,:,20)); axis equal xy tight off; colorbar;
m = makeMontage(fa_interp);
figure; imagesc(m); axis equal xy tight off; colorbar;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function IC = clusterMontage( I, IDX, nElem, pad )
% Used for visualization of clusters of images and videos.
%
% Organizes the matrix I, which is either MxNxR for images or MxNxTxR for
% videos, into a [M x N x T x nElem x nCluster] array IC, where each
% element IC(:,:,:,:,i) is the set of objects belonging to cluster i.
% If not all clusters have the same size, if pad==1 blank elements are
% added to pad the clusters so they do in fact have the same size, and if
% pad is 0 then IC is a cell vector where IC{i} is the set of objects
% belonging to cluster i.
%
% To display only a subset of clusters, given in a vector cs, use IDXb:
%  IDXb = zeros(size(IDX)); for i=1:length(cs) IDXb(IDX==cs(i))=i; end;
%
% USAGE
%  IC = clusterMontage( I, IDX, nElem, [pad] )
%
% INPUTS
%  I       - MxNxR array of images or MxNxTxR array of videos
%  IDX     - cluster membership (Rx1 integer vector) [see kmeans2.m]
%  nElem   - max number of instances to show of each cluster
%  pad     - [1] pads each cluster w blanks so it has exactly nElem
%
% OUTPUTS
%  IC      - if pad==1  [M x N x T x nElem x nCluster] array
%          - if pad==0  nCluster cell of [M x N x T x nElem_i] arrays
%
% EXAMPLE - images
%  load( 'images.mat' );
%  keep=randSample(144,80); IDXi=IDXi(keep); images=images(:,:,keep);
%  IC = clusterMontage( images, IDXi, 9, 0 );
%  figure(1); montage2( IC )
%
% EXAMPLE - videos
%  load( 'images.mat' );
%  IC = clusterMontage( videos, IDXv, 9, 0 );
%  figure(1); playMovie( IC )
%
% See also KMEANS2, MONTAGE2, PLAYMOVIE, CELL2ARRAY
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<4 || isempty(pad) ); pad=1; end

% error check
siz = size(I); nd = ndims(I);
if(nd~=3 && nd~=4); error('I must be 3 or 4 dimensional array'); end;
inds = {':'}; inds = inds(:,ones(1,nd-1));

% discard outliers
keepLocs = find( IDX>0 ); IDX = IDX(keepLocs); I=I(inds{:},keepLocs);

% sample both I and IDX so have nElem per cluster
uIDX=unique(IDX)';
for i=uIDX
  locs = find(IDX==i);  nLocs = length(locs);
  if( nLocs>nElem )
    keepLocs = [find(IDX~=i); locs(randSample(nLocs,nElem))];
    IDX = IDX(keepLocs); I=I(inds{:},keepLocs);
  elseif( nLocs<nElem && pad )
    nAdd = nElem-nLocs;
    IDX = [IDX; repmat(i,[nAdd,1])]; %#ok<AGROW>
    I = cat( nd, I, repmat(uint8(0),[siz(1:nd-1) nAdd]));
  end;
end;

% string out I
if( pad )
  IC = repmat( uint8(0), [siz(1:nd-1), nElem, length(uIDX)] );
  for i=uIDX; IC(inds{:},:,i) = I(inds{:},IDX==i); end
else
  IC = cell(1,max(IDX));
  for i=uIDX; IC{i} = I(inds{:},IDX==i); end
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function obj = setVertexProperties(obj, vertex, varargin)
% Sets the values of vertex properties
%
% Parameters:
%  vertex: the names of the vertices to be configured. It must be a cell string or a
%  character array.
%  @type cellstr
%  PropName: a character vector specifiy the name of a vertex property
%  for the added vertices. It must be given as a ''Name''-''Value'' pairs.
%  PropValue: the values of the vertex properties specified by the
%  ''PropName''. It must be scalar or a cell array with the the same length
%  of the number of the vertices to be added. If it is a scalar, all edges
%  will be assign this scalar value.
%  Use syntex
%  @verbatim sys =  addVertex(sys, Vert, PropName,PropValue,...) @endverbatim
%  @type char
%
% @note A struct can be used to replace the (PropName, PropValue) inputs.
% To be valid input, the structure must have fields with the same name as
% vertex properties. For mutilple nodes, use the array of structures instead of
% the structure of arrays. Invalid fields in the structure will be ignored,
% and missing properties will be added with their default values.

% check if the specified vertices exist in the graph
vert = findnode(obj.Gamma, vertex);
if any(find(~vert))
    error('The following vertices are NOT found in the graph: %s\n',...
        implode(vertex(~vert),', '));
end

nNode = numel(vert);

if ~ischar(varargin{1})
    error(['The thrid argument is expected to be a character vector.\n',...
        'Instead it is a %s'], class(varargin{1}));
end

% create a struct from the input ''Name-Value'' pair arguments
% create a struct array from the input
props = struct(varargin{:});
% if not a column array, convert to column array
if ~iscolumn(props)
    props = props';
end

if length(props) == 1 && nNode > 1
    % if the property value is a scalar, but the number of edges is greater
    % than one, create a structure array to replicate the property value
    for i=2:nNode
        props(i) = props(1);
    end
end

assert(length(props)==nNode,...
    'The number of property values ''%d'' must match the number of nodes ''%d''.\n',...
    length(props), nNode);
propNames = fieldnames(props);

valid_props = obj.VertexProperties;
% first check if these fields are valid property names
check_flags = cellfun(@(x)any(strcmp(x, valid_props.Name)), propNames);
if ~ all(check_flags)
    warning('The following input parameters are not valid properties, they will be not added to the graph: %s\n',...
        implode(propNames(~check_flags),', '));
end
% remove invalid fields
props = rmfield(props, propNames(~check_flags));

% get updated field names
propNames = fieldnames(props);

prop_pos = cellfun(@(x)find(strcmp(x, valid_props.Name)), propNames);

% go through all vertex properties to check if the value of a certain
% property is given
for i=1:numel(propNames)

    % if a property exists in the input argument, validate its
    % value
    HybridSystem.validatePropAttribute(propNames{i},...
        {props.(propNames{i})},...
        valid_props.Type{prop_pos(i)},valid_props.Attribute{prop_pos(i)});

    % if the value is not numeric or a cell, convert it to cell to
    % prevent future concatenation goes wrong
    for j=1:numel(vert)
        k = vert(j);
        if iscell(obj.Gamma.Nodes.(propNames{i})(j))
            obj.Gamma.Nodes.(propNames{i}){k} = props(j).(propNames{i});
        else
            obj.Gamma.Nodes.(propNames{i})(k) = props(j).(propNames{i});
        end
    end

end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev