text:
function x = jpeg2im(y)
%JPEG2IM Decodes an IM2JPEG compressed image.
%   X = JPEG2IM(Y) decodes compressed image Y, generating
%   reconstructed approximation X.  Y is a structure generated by
%   IM2JPEG.
%
%   See also IM2JPEG.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.4 $  $Date: 2003/10/26 18:39:08 $

error(nargchk(1, 1, nargin));             % Check input arguments

m = [16 11  10  16  24  40  51  61        % JPEG normalizing array
     12  12  14  19  26  58  60  55       % and zig-zag reordering
     14  13  16  24  40  57  69  56       % pattern.
     14  17  22  29  51  87  80  62
     18  22  37  56  68  109 103 77
     24  35  55  64  81  104 113 92
     49  64  78  87  103 121 120 101
     72  92  95  98  112 100 103 99];

order = [1 9  2  3  10 17 25 18 11 4  5  12 19 26 33  ...
        41 34 27 20 13 6  7  14 21 28 35 42 49 57 50  ...
        43 36 29 22 15 8  16 23 30 37 44 51 58 59 52  ...
        45 38 31 24 32 39 46 53 60 61 54 47 40 48 55  ...
        62 63 56 64];
rev = order;                          % Compute inverse ordering
for k = 1:length(order)
   rev(k) = find(order == k);
end

m = double(y.quality) / 100 * m;      % Get encoding quality.
xb = double(y.numblocks);             % Get x blocks.
sz = double(y.size);
xn = sz(2);                           % Get x columns.
xm = sz(1);                           % Get x rows.
x = huff2mat(y.huffman);              % Huffman decode.
eob = max(x(:));                      % Get end-of-block symbol

z = zeros(64, xb);   k = 1;           % Form block columns by copying
for j = 1:xb                          % successive values from x into
   for i = 1:64                       % columns of z, while changing
      if x(k) == eob                  % to the next column whenever
         k = k + 1;   break;          % an EOB symbol is found.
      else
         z(i, j) = x(k);
         k = k + 1;
      end
   end
end

z = z(rev, :);                                 % Restore order
x = col2im(z, [8 8], [xm xn], 'distinct');     % Form matrix blocks
x = blkproc(x, [8 8], 'x .* P1', m);           % Denormalize DCT
t = dctmtx(8);                                 % Get 8 x 8 DCT matrix
x = blkproc(x, [8 8], 'P1 * x * P2', t', t);   % Compute block DCT-1
x = uint8(x + 128);                            % Level shift

repo: weiyanmin


text:
function [A, B] = convertHypergraphToBipartiteGraph(S, printLevel)
% Converts a hypergraph into an undirected bipartite graph
%
% USAGE:
%
%    [A, B] = convertHypergraphToBipartiteGraph(S, printLevel)
%
% INPUT:
%    S:             `m x n` matrix with `m` nodes and `n` hyperedges
%    printLevel:    timer
%
% OUTPUTS:
%    A:             `(m+n) x (m+n)` adjacency matrix for an undirected graph (symmetric)
%    B:             `(m+n) x nnz(S)` incidence matrix for a directed graph
%
% .. Author: - Ronan Fleming 2013

if ~exist('printLevel', 'var')
    printLevel = 0;
end

[nMet, nRxn] = size(S);

if printLevel
    tic
end

S = sparse(S);
nnzS = nnz(S);

[row, col, v] = find(S);
numbers = 1:nnzS;
rowIndices = zeros(2 * nnzS, 1);
colIndices = zeros(2 * nnzS, 1);

values = zeros(2 * nnzS, 1);
rowIndices(1:2:end) = row;
colIndices(1:2:end) = numbers;
values(1:2:end) = v;

rowIndices(2:2:end) = nMet + col;
colIndices(2:2:end) = numbers;
values(2:2:end) = sign(v) * -1;

% incidence matrix for a bipartite graph
B = sparse(rowIndices, colIndices, values);
if printLevel
    toc
end

if printLevel
    tic
end

% create the adjacency matrix from undirected incidence matrix
A = inc2adj(B ~= 0);
if printLevel
    toc
end

% sanity check
assert(all(sum(B ~= 0, 1) == 2))
end

% old code
% %converts to bipartite hypergraph
% A=sparse(nMet+nRxn,nMet+nRxn);
%
% for j=1:nRxn
%     for i=1:nMet
%         if S(i,j)~=0
%             A(i,nMet+j)=1;
%             A(nMet+j,i)=1;
%         end
%     end
% end

% old methods
% switch method
% case 1
%     %incidence matrix for a bipartite graph
%     B=sparse(nMet+nRxn,nnzS);
%     k=1;
%     for j=1:nRxn
%         for i=1:nMet
%             if S(i,j)~=0
%                 if S(i,j)<0
%                     B(i,k)=S(i,j);
%                     B(nMet+j,k)=1;
%                 else
%                     B(i,k)=S(i,j);
%                     B(nMet+j,k)=-1;
%                 end
%                 k=k+1;
%
%             end
%         end
%     end
% case 2
%     rowIndices=zeros(2*nnzS,1);
%     colIndices=zeros(2*nnzS,1);
%     values=zeros(2*nnzS,1);
%     k=1;
%     p=1;
%     for j=1:nRxn
%         for i=1:nMet
%             if S(i,j)~=0
%                 if S(i,j)<0
%                     rowIndices(p)=i;
%                     colIndices(p)=k;
%                     values(p)=S(i,j);
%                     p=p+1;
%                     rowIndices(p)=nMet+j;
%                     colIndices(p)=k;
%                     values(p)=1;
%                     p=p+1;
%                 else
%                     rowIndices(p)=i;
%                     colIndices(p)=k;
%                     values(p)=S(i,j);
%                     p=p+1;
%                     rowIndices(p)=nMet+j;
%                     colIndices(p)=k;
%                     values(p)=-1;
%                     p=p+1;
%                 end
%                 k=k+1;
%             end
%         end
%     end

author: opencobra

repo: cobratoolbox


text:
function x = stretch_r ( r, i )

%*****************************************************************************80
%
%% STRETCH_R returns a data component given (R).
%
%  Discussion:
%
%    This routine shifts by 1 and stretches by 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the coordinate of a point that lies on the
%    boundary of the cube.
%
%    Input, integer I, the component of the data which is to be returned.
%
%    Output, real X, the I-th component of the data vector X, evaluated
%    at the point (R), which is on an endpoint of the unit line segment.
%
  if ( i == 1 )
    x = 2.0 * r + 1.0;
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'STRETCH_R - Fatal error!\n' );
    fprintf ( 1, '  Illegal component index I = %d\n', i );
    error ( 'STRETCH_R - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%                                 ADA = getada1(ADA, A,Ajc2,perm, d, blkstart)
% GETADA1  Compute ADA(i,j) = (D(d^2; LP,Lorentz)*A.t(:,i))' *A.t(:,j),
%   and exploit sparsity as much as possible.
%   Ajc2 points just beyond LP/Lorentz nonzeros for each column
%   blkstart = K.qblkstart partitions into Lorentz blocks.
%
%   IMPORTANT 1: only LP and sparse Lorentz part. PSD part ignored altogether.
%   For Lorentz, it uses only det(dk) * ai[k]'*aj[k].
%   IMPORTANT 2: Computes ADA only on triu(ADA(Aord.lqperm,Aord.lqperm)).
%     Remaining entries are set to 0. (CAUTION: sparse(ADA) will therefore
%     destroy the sparsity structure !).
%
% ******************** INTERNAL FUNCTION OF SEDUMI ********************
%
% See also sedumi, getada2, getada3

function ADA = getada1(ADA, A,Ajc2,perm, d, blkstart)
%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% 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 2 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, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%

disp('The SeDuMi binaries are not installed.')
disp('In Matlab, launch "install_sedumi" in the folder you put the SeDuMi files.')
disp('For more information see the file Install.txt.')
error(' ')

repo: REST


text:
% Description:
%
%     Create a design matrix for groups.
%
% Syntax:
%
%     [ X, terms ] = mG2X(groups, offset)
%
% Inputs:
%
%     groups  - [ N x G ] (int)  - columns of qualitative variables
%     offset  - [ 1 x 1 ] (int)  - group offset index for terms
%     options - [ 1 x P ] (cell) - see Details and Options
%     columns - [ 1 x T ] (int)  - indices into columns of groups
%
% Outputs:
%
%     X     - [ N x M ] (double)
%     terms - [ 1 x M ] (cell)
%
% Details:
%
%     By default, mG2X.m uses a canonical coding of the design matrix that
%     includes the first level of each grouping variable in the intercept.
%     This default coding may be replaced either by over-determined coding,
%     which includes a column of ones and zeros for every level of every
%     grouping variable, or by sigma-restricted coding, which estimates an
%     overall mean and offsets from the overall mean for every level of every
%     grouping variable, excluding the last level of each grouping variable.
%
% Options:
%
%     'over-determined'  - use over-determined coding for the design matrix
%     'sigma-restricted' - use sigma-restricted coding for the design matrix
%     'verbose'          - display extra information to the command window
%
% Examples:
%
% Notes:
%
% Author(s):
%
%     William Gruner (williamgruner@gmail.com)
%
% References:
%
% Acknowledgements:
%
%     Many thanks to Dr. Erik Erhardt and Dr. Elena Allen of the Mind Research
%     Network (www.mrn.org) for their continued collaboration.
%
% Version:
%
%     $Author: williamgruner $
%     $Date: 2010-04-12 07:14:07 -0600 (Mon, 12 Apr 2010) $
%     $Revision: 494 $

function [ X, terms ] = mG2X(groups, offset, options, columns)

    if ~exist('offset', 'var') || isempty(offset)
        offset = 0;
    end

    if ~exist('options', 'var')
        options = cell(0);
    end

    if ~exist('columns', 'var')
        columns = 1 : size(groups, 2);
    end

    terms = {};
    X     = [];

    if ~isempty(strmatch('verbose', options, 'exact'))
        fprintf('\n')
    end

    for i = columns

        [ B, I, J ] = unique(groups(:, i));

        if ~isempty(strmatch('verbose', options, 'exact'))
            fprintf('Factor %d represents column %d of groups and has %d levels.\n', ...
                i + offset, i, length(B))
        end

        if ~isempty(strmatch('over-determined', options, 'exact'))

            for j = 1 : length(B)
                X(J ~= j, length(terms) + j) =  0;
                X(J == j, length(terms) + j) =  1;
            end

            for j = 1 : length(B)
                terms{end + 1} = i + offset;
            end

        elseif ~isempty(strmatch('sigma-restricted', options, 'exact'))

            for j = 1 : length(B) - 1
                X(J ~=         j, length(terms) + j) =  0;
                X(J ==         j, length(terms) + j) =  1;
                X(J == length(B), length(terms) + j) = -1;
            end

            for j = 1 : length(B) - 1
                terms{end + 1} = i + offset;
            end

        else

            for j = 2 : length(B)
                X(J ~= j, length(terms) + j - 1) = 0;
                X(J == j, length(terms) + j - 1) = 1;
            end

            for j = 1 : length(B) - 1
                terms{end + 1} = i + offset;
            end

        end

    end

author: Sable

repo: mcbench-benchmarks


text:
function [W V] = genhull(U, wtype)
%GENHULL Generate "hull" matrices from orthonormal matrices
%	W = GENHULL(U)
%	W = GENHULL(U, wtype)
%	[W V] = GENHULL(U)
%	[W V] = GENHULL(U, wtype)
%
%	U     - matrices in a cell with orthonormal column vectors (U{i}' * U{i} = I)
%	wtype - weight function type for each matrix (default: 'snnn')
%
%	W     - "hull" matrices (weight functions)
%	V     - vertex matrices (W{i}*V{i} = U{i})
%
%	wtype parameter controls the type of the column vectors in the resulting
%	W matrices.
%	"hull" matrix means that the row sums are 1 and each element is positive
%		ie.: W{i}*ones == ones && all(all(W{i} >= 0)) &&
%		     exists V{i} such that W{i}*V{i} = U{i}
%
%	possible values of wtype:
%		'eye': return eye(size(U{i},1))
%		'ortho': keep as is: W{i} = U{i} (not hull matrix!)
%		'snnn': general non-negative sum-normalized (convex combination)
%		'cno': 'snnn' and each weight function comes close to 1 (tight hull)
%		'irno': 'snnn' and each weight function comes close to 0
%		'box': perform box_decomp
%
%	See also HOSVD, DECOMP.

if nargin <= 1
	% TODO
	wtype = 'close';
end

if ~iscell(U)
	% only one matrix
	if nargout < 2
		W = decomp(U, wtype);
	else
		[W V] = decomp(U, wtype);
	end
	return
end

if ~iscell(wtype)
	tmp = wtype;
	wtype = cell(1,length(U));
	for i = 1:length(U)
		wtype{i} = tmp;
	end
end

% TODO: allow skipping last few U?
%assert(length(U) == length(wtype));

W = cell(1,length(wtype));
if nargout < 2
	for i = 1:length(wtype)
		if isempty(U{i})
			W{i} = [];
		else
			W{i} = decomp(U{i}, wtype{i});
		end
	end
else
	V = cell(1,length(wtype));
	for i = 1:length(wtype)
		if isempty(U{i})
			W{i} = [];
			V{i} = [];
		else
			[W{i} V{i}] = decomp(U{i}, wtype{i});
		end
	end
end


text:
function [p,q,D,sc] = dpfast(M,C,T,G)
% [p,q,D,sc] = dpfast(M,C,T,G)
%    Use dynamic programming to find a min-cost path through matrix M.
%    Return state sequence in p,q; full min cost matrix as D and
%    local costs along best path in sc.
%    This version gives the same results as dp.m, but uses dpcore.mex
%    to run ~200x faster.
%    C is a step matrix, with rows (i step, j step, cost factor)
%    Default is [1 1 1.0;0 1 1.0;1 0 1.0];
%    Another good one is [1 1 1;1 0 1;0 1 1;1 2 2;2 1 2]
%    T selects traceback origin: 0 is to any edge; 1 is top right (default);
%    T > 1 finds path to min of anti-diagonal T points away from top-right.
%    Optional G defines length of 'gulleys' for T=0 mode; default 0.5
%    (i.e. accept path to only 50% of edge nearest top-right)
% 2003-04-04,2005-04-04 dpwe@ee.columbia.edu $Header: /Users/dpwe/projects/dtw/RCS/dpfast.m,v 1.6 2008/03/14 14:40:50 dpwe Exp dpwe $

% Copyright (c) 2003 Dan Ellis <dpwe@ee.columbia.edu>
% released under GPL - see file COPYRIGHT

if nargin < 2
  % Default step / cost matrix
  C = [1 1 1.0;0 1 1.0;1 0 1.0];
end

if nargin < 3
  % Default: path to top-right
  T = 1;
end

if nargin < 4
  % how big are gulleys?
  G = 0.5;  % half the extent
end

if sum(isnan(M(:)))>0
  error('dpwe:dpfast:NAN','Error: Cost matrix includes NaNs');
end

if min(M(:)) < 0
  disp('Warning: cost matrix includes negative values; results may not be what you expect');
end

[r,c] = size(M);

% Core cumulative cost calculation coded as mex
[D,phi] = dpcore(M,C);

p = [];
q = [];

%% Traceback from top left?
%i = r;
%j = c;

if T == 0
  % Traceback from lowest cost "to edge" (gulleys)
  TE = D(r,:);
  RE = D(:,c);
  % eliminate points not in gulleys
  TE(1:round((1-G)*c)) = max(max(D));
  RE(1:round((1-G)*r)) = max(max(D));
  if (min(TE) < min(RE))
    i = r;
    j = max(find(TE==min(TE)));
  else
    i = max(find(RE==min(RE)));
    j = c;
  end
else
  if min(size(D)) == 1
    % degenerate D has only one row or one column - messes up diag
    i = r;
    j = c;
  else
    % Traceback from min of antidiagonal
    %stepback = floor(0.1*c);
    stepback = T;
    slice = diag(fliplr(D),-(r-stepback));
    [mm,ii] = min(slice);
    i = r - stepback + ii;
    j = c + 1 - ii;
  end
end

p=i;
q=j;

sc = M(p,q);

while i > 1 & j > 1
%  disp(['i=',num2str(i),' j=',num2str(j)]);
  tb = phi(i,j);
  i = i - C(tb,1);
  j = j - C(tb,2);
  p = [i,p];
  q = [j,q];
  sc = [M(i,j),sc];
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function v_cheb = legcoeffs2chebvals(c_leg, varargin)
%LEGCOEFFS2CHEBVALS  Convert Legendre coefficients to Chebyshev values.
%   V_CHEB = LEGCOEFFS2CHEBVALS(C_LEG) converts the vector C_LEG of Legendre
%   coefficients to a vector V_CHEB of values at second-kind Chebyshev points.
%
%   V_CHEB = LEGCOEFFS2CHEBVALS(C_LEG, 1) is similar, but evaluates at
%   first-kind Chebyshev points.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Parse inputs:
v1 = {}; v2 = {};
for k = 1:numel(varargin)
    if ( isnumeric(varargin{k}) )
        v1{1} = varargin{k};
    elseif ( ischar(varargin{k}) )
        v2{1} = varargin{k};
    end
end

c_cheb = leg2cheb(c_leg, v2{:}); % Convert to Chebyshev coefficients.
v_cheb = chebcoeffs2chebvals(c_cheb, v1{:});  % Convert to Chebyshev values.

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function strategy = beatTheBookie(dat, bet, marg)

money = 0; % initial amount of money
m = 1; % counter for money
accuracy = [];
acc = 1;
nGames = size(dat, 1);
bet_games = [];
nValidOdds = 3;

for gm = 1 : nGames

    score_home = dat(gm, 1);
    score_away = dat(gm, 2);
    if score_home > score_away
        result = 1;
    elseif score_home == score_away
        result = 2;
    else
        result = 3;
    end

    averages = dat(gm, 3:5); % [avg_home avg_draw avg_away];
    maximums = dat(gm, 6:8); % [max_home max_draw max_away];
    counts = dat(gm, 9:11);

    % check that there is at least one valid bet. If the bets are still not
    % online for this time, continue to the next game
    if any(isnan(averages)) || any(isnan(maximums))
        continue
    end

    % Apply formula and decide whether to bet
    earn_margin(1) = ((1 ./ averages(1) - marg) * maximums(1) - 1) * (counts(1) > nValidOdds) ;
    earn_margin(2) = ((1 ./ averages(2) - marg) * maximums(2) - 1) * (counts(2) > nValidOdds);
    earn_margin(3) = ((1 ./ averages(3) - marg) * maximums(3) - 1) * (counts(3) > nValidOdds);

    if sum(earn_margin > 0) >= 1

        [~, bet_result] = max(earn_margin);
        possible_earn = bet  * (maximums(bet_result) - 1);

        max_odds(m) = maximums(bet_result); %#ok
        mean_odds(m) = averages(bet_result); %#ok
        ids(m) = bet_result;

        % calculate loss / earning
        if isequal(bet_result, result)
            money(m + 1) = money(m) + possible_earn;
            accuracy(acc) = 1;
        else
            money(m + 1) = money(m) - bet;
            accuracy(acc) = 0;

        end

        bet_games(m) = gm;
        m = m + 1;
        acc = acc + 1;
    end

end

strategy.money = money;
strategy.odds = max_odds;
strategy.mean_odds = mean_odds;
strategy.accuracy = accuracy;
strategy.ids = ids;

end

repo: BeatTheBookie


text:
% SATGLOBE4 - Draw an idealized satellite view of earth
%
% This file renders a fully manipulatable satellite view
% of earth at a resolution of four pixels per degree, with added
% international political boundaries and gridlines.
% The imagery was obtained from NASA, then postprocessed.
% The globe was rendered using the Matlab Mapping Toolbox.
%
% The Mapping Toolbox is not needed to use this file; however,
% if you have the toolbox, you will be able to use the plot3m
% command to add your own graphics. If not, you can simply use plot3.
%
% In order to save storage space, this m-file loads image
% data from the file satglobe.mat, and then creates the
% graticule mesh itself. This process allows users who
% do not have the Matlab Mapping Toolbox to render the
% figure, but it does take a few moments to compute the
% mesh. Using this trick, the data storage is reduced
% considerably; however, once the figure is rendered, you
% may wish to save it as a regular Matlab figure file
% to increase speed.
%
% Michael Kleder, 2004

function satglobe4
load satglobe4
th=repmat((0:.25:180)'*pi/180,[1 1441]);
ph=repmat((-180:.25:180)*pi/180,[721 1]);
s.children(1).properties.XData = sin(th).*cos(ph);
s.children(1).properties.YData = sin(th).*sin(ph);
s.children(1).properties.ZData = cos(th);
s.children(1).properties.CData = double(c)/255;
figure;
struct2handle(s,gcf);
set(gcf,'color','k','renderer','zbuffer','inverthardcopy','off','name',...
   'Earth at 4 Pixels per Degree, by Michael Kleder');
return


text:
function R=Rxyz(t,flag)

if flag==1
    R=[1       0         0;
       0     cos(t)   -sin(t);
       0     sin(t)    cos(t)]';
elseif flag==2
     R=[cos(t)       0         sin(t);
           0         1           0;
       -sin(t)       0         cos(t)]';
elseif flag==3
    R=[   cos(t)      -sin(t)       0;
          sin(t)       cos(t)       0;
           0            0           1]';
end

return;

repo: GINav


text:
% function [xpad, nup, n1, n2] = padsignal(x, padtype, padlength)
%
% Pads signal and returns indices of original signal
%
% Input:
%	x: original signal
%	padtype (optional): either 'symmetric' (default) or 'replicate'
%   padlength (optional): number of samples to pad on each side; default is nearest power of 2
% Output:
%	x: padded signal
%   nup: next power of 2
%   n1: length on left
%   n2: length on right
%
%---------------------------------------------------------------------------------
%    Synchrosqueezing Toolbox
%    Authors: Eugene Brevdo, Gaurav Thakur
%---------------------------------------------------------------------------------
function [xpad, nup, n1, n2] = padsignal(x, padtype, padlength)

[nup, n1, n2] = p2up(length(x));
if nargin<3
	[nup, n1, n2] = p2up(length(x));	%if padlength not given, pad up to nearest power of 2
else
	nup = length(x)+2*padlength;
	n1 = padlength-1;
	n2 = padlength;
end

%xl = padarray(x(:), n1, padtype, 'pre');
%xr = padarray(x(:), n2, padtype, 'post');
%xpad = [xl(1:n1); x(:); xr(end-n2+1:end)];
%padarray needs image processing toolbox; below is equivalent code to avoid that

x=x(:);
n=length(x);
if strcmpi(padtype,'symmetric')
	xl = repmat([x;flipud(x)],[ceil(n1/(2*n)),1]);
	xr = repmat([flipud(x);x],[ceil(n2/(2*n)),1]);
elseif strcmpi(padtype,'replicate')
	xl = x(1)*ones(n1,1);
	xr = x(end)*ones(n2,1);
end

xpad = [xl(end-n1+1:end); x; xr(1:n2)];

author: ebrevdo

repo: synchrosqueezing


text:
function [A,d]=affineTransBetweenTriangles(v1,v2)
%%AFFINETRANSBETWEENTRIANGLES Given the vertices of one triangle in 2D and
%       the vertices of a second triangle in 2D, obtain a matrix M and a
%       vector d such that p2=M*p1+d transforms a point p1 in the first
%       triangle into a point p2 in the second triangle.
%
%INPUTS: v1 The 2X3 set of vertices in the first triangle.
%        v2 The 2X3 set of vertices in the second triangle (into which we
%           wish to transform the first triangle).
%
%OUTPUTS: A A 2X2 scale/rotation matrix.
%         d A 2X1 origin shift vector.
%
%Equations for the transformation are given in Chapter 5 of [1].
%
%EXAMPLE:
%This just find the transformation between two sets of vertices and
%transforms the first set of vertices into the second set. The error is on
%the order of typical finite precision limitations.
% v1=[ 5, -22,   3;
%     18,   8, -13];
% v2=[-4, 35, -13;
%      3, 27,  30];
% [A,d]=affineTransBetweenTriangles(v1,v2);
% v2Trans=bsxfun(@plus,A*v1,d);
% RelErr=max(abs((v2Trans-v2)./v2))
%
%REFERENCES:
%[1] K. Anjyo and H. Ochiai, Mathematical Basics of Motion and Deformation
%    in Computer Graphics, 2nd ed. Morgan and Claypool Publishers, 2017.
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

A1=[v1;
   1,1,1];
A2=[v2;
   1,1,1];

M=A2/A1;

A=M(1:2,1:2);
d=M(1:2,3);

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function dotmfolders = findDotMFolders(subdir, exclList)
global MAXPATHLENGTH
MAXPATHLENGTH = 8;

dotmfolders = {};

if ~exist('subdir','var')
    subdir = pwd;
end
if ~exist('exclList','var')
    exclList = {};
end

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

subdirFullpath = filesepStandard_startup(subdir,'full');

if ~ispathvalid_startup(subdirFullpath, 'dir')
    fprintf('Warning: folder %s doesn''t exist\n', subdirFullpath);
    return;
end

% If current subjdir is in the exclList then go back to curr dir and exit
if isExcluded(subdirFullpath, exclList)
    return;
end

dirs = dir([subdirFullpath, '*']);
if isempty(dirs)
    return;
end

if isdotmfolder(subdirFullpath)
    dotmfolders = {filesepStandard_startup(subdirFullpath, 'nameonly')};
end

for ii = 1:length(dirs)
    if ~dirs(ii).isdir
        continue;
    end
    if dirs(ii).name(1) == '.'
        continue;
    end
    dotmfolders = [dotmfolders; findDotMFolders([subdirFullpath, dirs(ii).name], exclList)]; %#ok<AGROW>
end

% -------------------------------------------------------------------------
function b = isdotmfolder(folder)
global MAXPATHLENGTH
b = false;
if ~ispathvalid_startup(folder, 'dir')
    return
end
if isempty(dir([folder,'/*.m']))
    % Exceptions to rule that 'dotm' folder must have at least one '.m' file:
    % it is a an executable folder (i.e. '/bin')
    if ~isempty(strfind(folder, '/bin/')) %#ok<*STREMP>
        b = true;
        return
    end
    return;
else
    rootdir = which('findDotMFolders');
    rootdir = fileparts(rootdir);
    rootdir = pathsubtract_startup(rootdir, 'Utils/submodules','nochange');
    p = pathsubtract_startup(folder, rootdir);
    if length(find(p=='/')) > MAXPATHLENGTH
        return
    end
end
b = true;

% -------------------------------------------------------------------------
function b = isExcluded(pname, exclList)
b = true;
if pname(end)=='/'
    pname(end) = '';
end
if ~ispathvalid_startup(pname,'dir')
    return;
end
[~,f,e] = fileparts(pname);
for ii = 1:length(exclList)
    [c,d] = str2cell_startup(exclList{ii},'*');
    if isempty(d) && strcmp(c{1}, [f,e])
        return;
    end

    % Get list of all folders matching exclList{ii} pattern, whether it be
    % a single folder name or a wildcard pattern
    for jj = 1:length(c)
        k = strfind(c{jj}, [f,e]);
        if isempty(k)
            break;
        end
    end
    if ~isempty(k)
        return;
    end
end
b = false;

author: BUNPC


text:
%kvlocmin 'Regional Min'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vlocmin.pane file
%
% Parameters:
% InputFile: i 'Input Image ', required: 'Input image'
% InputFile: str 'Struct Element', required: 'Structuring Element'
% OutputFile: o 'Output Image', required: 'resulting output image'
%
% Example: o = kvlocmin({i, str}, {'i','';'str','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vlocmin - Regional Min
%
%  DESCRIPTION
% .I vlocmin
% will find the local minimum of an image.
%
%
%
%
%  EXAMPLES
% vlocmin -i feath.xv -str b.str -o outimage.xv
%
% Will find the local minimum of the input image "feath.xv". The connectivity is defined by the structuring element b. The transformed image is written to "outimage.xv".
%
%  "SEE ALSO"
% vcast1(1)
%
%  RESTRICTIONS
% .I vlocmin
% can be defined for all data types supported by Khoros (except bit), but at the moment it has been implemented just for the unsigned char type.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993-1997 Junior Barrera, Roberto Lotufo.  All rights reserved.
%

function varargout = kvlocmin(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,..] = kvlocmin(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'str', '__input';'o', '__output'};
maxval={0,0,0};
minval={0,0,0};
istoggle=[0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','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 'vlocmin"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [z0,deltaZ]=addAstroRefrac(algorithm,plhObs,zTrue,Rh,P,T,wl)
%%ADDASTROREFRAC Add the effects of refraction to the true zenith angle of
%                an object outside the atmosphere using low-precision
%                atmospheric models for an observer near the surface of
%                the Earth viewing an object using a narrowband sensor.
%
%INPUTS: algorithm This specified the algorithm used. The possible values
%                  are:
%                  0: Use a numerical integration method from Chapter 7.2
%                     of [1] and from [2]. This algorithm uses the inputs
%                     plhObs,z0, Rh, P, T, and wl and makes use of the
%                     Sinclair atmospheric model. This algorithm is the
%                     most precise of all of the methods. If a point is too
%                     far below the horizon, z0 and deltaZ will be empty
%                     matrices.
%                  1: Use the simple formula of [3] for the refraction
%                     experienced by an observer at sea level viewing light
%                     with a wavelength of 0.574 micrometers (yellow). This
%                     algorithm uses the inputs z0, Rh, P, T. This
%                     algorithm is only valid for refraction-corrupted
%                     zenith distances below 70 degrees. Values that are
%                     too high will lead to z0 and deltaZ being empty
%                     matrices.
%                  2: Use the algorithm from the International Astronomical
%                     Union's (IAU) standards of fundamental astronomy
%                     library. This algorithm uses the inputs z0,Rh,P,T,wl
%                     for an observer at sea level. This algorithm will
%                     produce results for all positive zenith distances,
%                     though the results might not be very good for large
%                     values.
%           plhObs The WGS-84 ellipsoidal latitude, longitude and height
%                  of the observer in meters. Only algorithm 0 uses this
%                  parameter, and only the latitude and height above the
%                  reference ellipsoid are used (the latitude does not
%                  matter). The height is treated as an approximate height
%                  above mean sea level. Given the precision of the
%                  algorithm, the fact that geoid undulations are ignored
%                  probably does not matter.For algorithms 1 and 2, an
%                  empty matrix can be passed.
%            zTrue A vector or matrix of true (refraction-free) positive
%                  zenith "distances" (in radians) of the celestial object
%                  being observed. A zenith distance is the observed angle
%                  of an object down from the gravitational vertical (the
%                  zenith) at the observer's location. This value is
%                  inaccurate as one approaches pi/2 radians (the horizon).
%                  zTrue must be greater than zero. Given the precision of
%                  the algorithms that are available, an angle with respect
%                  to the vertical defined by the WGS-84 reference
%                  ellipsoid could probably be substituted for an angle
%                  with respect to the true gravitational vertical.
%               Rh The relative humidity at the observer (between 0 and 1).
%                  If this parameter is omitted or an empty matrix is
%                  passed, then Constants.standardRelHumid is used.
%                P The atmospheric pressure at the observer in Pascals
%                  (N/m^2). If this parameter is omitted or an empty matrix
%                  is passed, then Constants.standardAtmosphericPressure is
%                  used.
%                T The air temperature at the observer in degrees Kelvin.
%                  If this parameter is omitted or an empty matrix is
%                  passed, then Constants.standardTemp is used.
%               wl The wavelength at which the observation is made in units
%                  of meters. If this parameter is omitted or an empty
%                  matrix is passed, then a wavelength of 0.574 micrometers
%                  is used, which is in the visible spectrum (a rather
%                  yellow color).
%
%OUTPUTS: z0  An NX1 vector of the zTrue values in radians with atmospheric
%             refraction added.
%      deltaZ The refraction value that was applied. z0=zTrue-deltaZ.
%
%The function removeAstroRefrac computes deltaZ for the problem where z0 is
%known and zTrue is unknown. This function iterates the solution a few
%times to try to solve the inverse problem. A fixed 20 iterations are used,
%which is generally sufficient to ensure convergence to within working
%precision limits for all of the algorithms.
%
%REFERENCES:
%[1] S. E. Urban and K. P.Seidelmann, Eds.,Explanatory Supplement to the
%    Astronomical Almanac, 3rd ed. Mill Valley, CA: University Science
%    Books, 2013.
%[2] C. Y. Hohenkerk and A. T. Sinclair, "The computation of angular
%    atmospheric refraction at large zenith angles," United Kingdom
%    Hydrographic Office, HM Nautical Almanac, Tech. Rep. 63, Apr. 1985.
%    http://astro.ukho.gov.uk/data/tn/naotn63.pdf
%
%April 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if(nargin<7||isempty(wl))
       wl=0.574e-6;
    end

    if(nargin<6||isempty(T))
       T=Constants.standardTemp;
    end

    if(nargin<5||isempty(P))
        P=Constants.standardAtmosphericPressure;
    end

    if(nargin<4||isempty(Rh))
        Rh=Constants.standardRelHumid;
    end

    %The initial estimate of deltaZ is given by solving the inverse problem
    %at zTrue
    [~,deltaZ]=removeAstroRefrac(algorithm,plhObs,zTrue,Rh,P,T,wl);

    if(~isempty(deltaZ))
        numIter=20;
        for curIter=1:numIter
            [~,deltaZ]=removeAstroRefrac(algorithm,plhObs,zTrue-deltaZ,Rh,P,T,wl);
            if(isempty(deltaZ))
                %This can arise if the observation ends up too far
                %underground.
                break
            end
        end
    end
    if(isempty(deltaZ))
        z0=[];
    else
        z0=zTrue-deltaZ;
    end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
function [xUpdate,LUpdate,PUpdate,innov,Pzz,W]=separatedCovUpdate(xPred,LPred,TPred,z,R,H,c)
%%SEPARATEDCOVUPDATE Perform the measurement update step in the separated
%                    covariance filter. This filter takes a maximum assumed
%                    acceleration (or other moment) for the target and
%                    provides the optimum estimate in terms of a cost
%                    function that trades off between estimation accuracy
%                    and estimator delay. The filter assumes a linear
%                    measurement model.
%
%INPUTS: xPred The xDimX1 predicted target state.
%        LPred The xDimXzDim predicted delay vector (defined before
%              Equation 12 in [1]). The use of multiple columns represents
%              a choice in how the algorithm was generalized to multiple
%              dimensions.
%        TPred The xDimXxDim total error matrix. This is a combination of
%              errors due to measurement noise and filter lag.
%            z The zDimX1 vector measurement. zDim should be the same as
%              the number of position components in the state.
%            R The zDimXzDim measurement covariance matrix.
%            H The zDimXxDim measurement matrix for a linear measurement
%              model. That is z=H*x+w, where w is measurement noise having
%              covariance matrix R.
%            c The confidence region under consideration by the filter.
%              0<c<1. If this parameter is omitted or an empty matrix is
%              passed, the default value of c=0.99 is used.
%
%OUTPUTS: xUpdate The xDimX1 updated state vector.
%         LUpdate The xDimXzDim updated delay vector.
%         PUpdate The xDimXxDim covariance matrix of the state estimate.
%                 This is a combination of LUpdate and PUpdate.
%      innov, Pzz The zDimX1 innovation and a zDimXzDim matrix S that is
%                 akin to an innovation covariance matrix are returned in
%                 case one wishes to analyze the consistency of the
%                 estimator or use those values in gating or likelihood
%                 evaluation.
%               W The gain used in the update. This can be useful when
%                 gating and using the function calcMissedGateCov.
%
%The equations for the algorithm are given in Table 1 of [1]. The algorithm
%is presented in 1D for a state consisting of position and velocity.
%However, the equations are given in vector form and can thus be used in
%multiple dimensions, which is done here.
%
%The filter in 1 is 1D. To generalize the filter to 3D, L is redefined so
%that each column contains the lag for one particular dimension of motion.
%To keep the solution the same as in Table 1, L becomes a matrix where the
%elements in each column that do not correspond to components for that
%dimensions of motion are zero.
%
%In [1], no clear method of initializing this type of tracking filter is
%provided. A simple way to initialize the filter would be to use two
%Cartesian converted measurements to obtain a state estimate and covariance
%PInit as one would do with a normal Kalman filter (one could, for example,
%use the KalmanFIRSmoother function) and then set TUpdate=c^2*PInit;
%LUpdate=zeros(xDim,zDim); and PUpdate=PInit.
%
%REFERENCES:
%[1] G. J. Portmann, J. R. Moore, and W. G. Bath, "Separated covariance
%    filtering," in Proceedings of the IEEE International Radar Conference,
%    Arlington, VA, 7-10 May 1990, pp. 456-460.
%
%September 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<7||isempty(c))
   c=0.99;
end

Pzz=H*TPred*H'+c^2*R;
%Ensure symmetry
Pzz=(Pzz+Pzz')/2;

W=TPred*H'/Pzz;%The gain

innov=z-H*xPred;%The innovation
xUpdate=xPred+W*innov;

xDim=size(xPred,1);
diff=eye(xDim,xDim)-W*H;

LUpdate=diff*LPred;
TUpdate=diff*TPred;

%Ensure symmetry
TUpdate=(TUpdate+TUpdate')/2;

PUpdate=(TUpdate-LUpdate*LUpdate')/c^2;

%Ensure symmetry
PUpdate=(PUpdate+PUpdate')/2;
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
function test_bug27

% MEM 8gb
% WALLTIME 00:10:00
% DEPENDENCY ft_preproc_lowpassfilter ft_resampledata ft_resampledata ft_datatype_raw fixtimeaxes

% this script tests for bug 27 and for the solution
% the 'bug' is that if ft_resampledata is called on data with different time axes
% after resampling the time axes of the trials may be shifted with respect to one
% another with a fractional amount (governed by the input fsample)

%% generate some data
% shifted by one sample
data1          = [];
data1.trial{1} = randn(2,1000);
data1.trial{2} = randn(2,1000);
data1.time{1}  = [0:0.001:0.999] - 0.5;
data1.time{2}  = [0.001:0.001:1] - 0.5;
data1.label    = {'chan1';'chan2'};
data1.fsample  = 1000;
data1.cfg.trl  = zeros(2,3);

for k = 1:2
  data1.trial{k} = ft_preproc_lowpassfilter(data1.trial{k}, 1000, 40);
end

%% resample with a factor 1x
cfg            = [];
cfg.resamplefs = 1000;
cfg.detrend    = 'no';
data2          = ft_resampledata(cfg, data1);

[ok, msg] = isalmostequal(data1.time, data2.time, 'abstol', 1e-10);
assert(ok, msg);

%% resample
cfg            = [];
cfg.resamplefs = 256;
cfg.detrend    = 'no';
data2          = ft_resampledata(cfg, data1);

figure
subplot(2,2,1), plot(data1.time{1}, data1.trial{1}, 'b'); xlim([data1.time{1}(1) data1.time{1}(end)])
subplot(2,2,2), plot(data1.time{2}, data1.trial{2}, 'b'); xlim([data1.time{2}(1) data1.time{2}(end)])
subplot(2,2,3), plot(data2.time{1}, data2.trial{1}, 'r'); xlim([data1.time{1}(1) data1.time{1}(end)])
subplot(2,2,4), plot(data2.time{2}, data2.trial{2}, 'r'); xlim([data1.time{2}(1) data1.time{2}(end)])

%% make it a bit more extreme
data1.time{2} = data1.time{2} + 3.1152;
data2 = ft_resampledata(cfg, data1);

figure
subplot(2,2,1), plot(data1.time{1}, data1.trial{1}, 'b'); xlim([data1.time{1}(1) data1.time{1}(end)])
subplot(2,2,2), plot(data1.time{2}, data1.trial{2}, 'b'); xlim([data1.time{2}(1) data1.time{2}(end)])
subplot(2,2,3), plot(data2.time{1}, data2.trial{1}, 'r'); xlim([data1.time{1}(1) data1.time{1}(end)])
subplot(2,2,4), plot(data2.time{2}, data2.trial{2}, 'r'); xlim([data1.time{2}(1) data1.time{2}(end)])

%% make it a bit more extreme
cfg.resamplefs = 17;
data2 = ft_resampledata(cfg, data1);

figure
subplot(2,2,1), plot(data1.time{1}, data1.trial{1}, 'b'); xlim([data1.time{1}(1) data1.time{1}(end)])
subplot(2,2,2), plot(data1.time{2}, data1.trial{2}, 'b'); xlim([data1.time{2}(1) data1.time{2}(end)])
subplot(2,2,3), plot(data2.time{1}, data2.trial{1}, 'r'); xlim([data1.time{1}(1) data1.time{1}(end)])
subplot(2,2,4), plot(data2.time{2}, data2.trial{2}, 'r'); xlim([data1.time{2}(1) data1.time{2}(end)])

%%

% this contains a problematic data structure from Mats
load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug27.mat'), 'data');

clear sample0 time0
for i=1:numel(data.time)
  % find the sample corresponding with time zero
  sample0(i)  = nearest(data.time{i}, 0);
  time0(i)    = data.time{i}(sample0(i));
end
figure
plot(time0, '.-'); % this looks ok, no jitter
assert(std(time0)<1e-6, 'there is too much jitter in time0');

cfg = [];
cfg.resamplefs = 600;
data_resampled = ft_resampledata(cfg, data);

clear sample0 time0
for i=1:numel(data.time)
  % find the sample corresponding with time zero
  sample0(i)  = nearest(data_resampled.time{i}, 0);
  time0(i)    = data_resampled.time{i}(sample0(i));
end
figure
plot(time0, '.-'); % this showed quite some jitter over trials
assert(std(time0)<1e-6, 'there is too much jitter in time0');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = get_soln(om, set_type, tags, name, idx)
%GET_SOLN  Fetch solution values for specific named/indexed sets.
%   VALS = OM.GET_SOLN(SET_TYPE, NAME)
%   VALS = OM.GET_SOLN(SET_TYPE, NAME, IDX)
%   VALS = OM.GET_SOLN(SET_TYPE, TAGS, NAME)
%   VALS = OM.GET_SOLN(SET_TYPE, TAGS, NAME, IDX)
%
%   Returns named/indexed results for a solved model, evaluated at
%   the solution found.
%
%   Inputs:
%       SET_TYPE - one of the following, specifying the type of set:
%           'var' - variables
%           'lin' - linear constraints
%           'nle' - nonlinear equality constraints
%           'nli' - nonlinear inequality constraints
%           'nlc' - nonlinear costs
%           'qdc' - quadratic costs
%       TAGS - char array or cell array of char arrays specifying the
%           desired output(s). Valid tags vary by SET_TYPE as follows:
%           'var' - default is {'x', 'mu_l', 'mu_u'}
%               'x' - value of solution variable
%               'mu_l' - shadow price on variable lower bound
%               'mu_u' - shadow price on variable upper bound
%           'lin' - default is {'g', 'mu_l', 'mu_u'}
%               'g' - 1 x 2 cell array of upper and lower constraint
%                   values, {A*x - u, l - A*x}
%               'Ax_u' - upper constraint value, A*x - u
%               'l_Ax' - lower constraint value, l - A*x
%               'mu_l' - shadow price on constraint lower bound
%               'mu_u' - shadow price on constraint upper bound
%           'nle' - default is {'g', 'lam', 'dg'}
%               'g' - constraint value g(x)
%               'lam' - shadow price on constraint
%               'dg' - Jacobian of constraint
%           'nli' - default is {'h', 'mu', 'dh'}
%               'h' - constraint value h(x)
%               'mu' - shadow price on constraint
%               'dh' - Jacobian of constraint
%           'nlc' and 'qdc' - default is {'f', 'df', 'd2f'}
%               'f' - cost function value f(x) (for 'qdc' can return a vector)
%               'df' - gradient of cost function
%               'd2f' - Hessian of cost function
%       NAME - char array specifying the name of the set
%       IDX  - cell array specifying the indices of the set
%
%   Outputs:
%       Variable number of outputs corresponding to TAGS input. If TAGS
%       is empty or not specified, the calling context will define the
%       number of outputs, returned in order of default tags for the
%       specified SET_TYPE.
%
%   Examples:
%       [P, muPmin, muPmax] = om.get_soln('var', 'P');
%       [mu_u, mu_l] = om.get_soln('lin', {'mu_u', 'mu_l'}, 'lin_con_1');
%       dg_b_2_3 = om.get_soln('nle', 'dg', 'nle_con_b', {2,3});
%
%   For a complete set of solution vector values and shadow prices, using
%   the PARSE_SOLN method may be more efficient.
%
%   See also PARSE_SOLN.

%   MP-Opt-Model
%   Copyright (c) 2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

%% input arg handling
if nargin == 3              %% om.get_soln(set_type, name)
    idx = [];
    name = tags;
    tags = {};
elseif nargin == 4
    if ischar(name)         %% om.get_soln(set_type, tags, name)
        idx = [];
    else                    %% om.get_soln(set_type, name, idx)
        idx = name;
        name = tags;
        tags = {};
    end
end

%% set up tags for default outputs
if isempty(tags)
    switch set_type
        case 'var'
            tags = {'x', 'mu_l', 'mu_u'};
        case 'lin'
            if strcmp(om.problem_type(), 'LEQ')
                tags = {'f'};   %% 'Ax_u', 'l_Ax' are also options
            else
                tags = {'g', 'mu_l', 'mu_u'};   %% 'Ax_u', 'l_Ax' are also options
            end
        case 'nle'
            tags = {'g', 'lam', 'dg'};
        case 'nli'
            tags = {'h', 'mu', 'dh'};
        case 'nlc'
            tags = {'f', 'df', 'd2f'};
        case 'qdc'
            tags = {'f', 'df', 'd2f'};
    end
elseif ~iscell(tags)
    tags = { tags };
end

%% set up indexing
om_ff = om.(set_type);
if isempty(idx)         %% simple named set
    N = om_ff.idx.N.(name);
    i1 = om_ff.idx.i1.(name);           %% starting row index
    iN = om_ff.idx.iN.(name);           %% ending row index
else                    %% indexed named set
    %% calls to substruct() are relatively expensive, so we pre-build the
    %% structs for addressing cell and numeric array fields, updating only
    %% the subscripts before use
    sc = struct('type', {'.', '{}'}, 'subs', {name, idx});  %% cell array field
    sn = sc; sn(2).type = '()';                         %% num array field
    N = subsref(om_ff.idx.N, sn);
    i1 = subsref(om_ff.idx.i1, sn);     %% starting row index
    iN = subsref(om_ff.idx.iN, sn);     %% ending row index
end

%% get outputs
varargout = cell(1, nargout);
s = om.soln;
if N && ~isempty(s.eflag)
    switch set_type
        case 'var'
            for k = 1:nargout
                switch tags{k}
                    case 'x'
                        varargout{k} = s.x(i1:iN);
                    case 'mu_l'
                        varargout{k} = s.lambda.lower(i1:iN);
                    case 'mu_u'
                        varargout{k} = s.lambda.upper(i1:iN);
                    otherwise
                        error('opt_model/get_soln: unknown tag ''%s''', tags{k});
                end
            end
        case 'lin'
            if strcmp(om.problem_type(), 'LEQ') %% tag must be 'f'
                varargout{1} = s.f(i1:iN);
            else
                if any(ismember({'g', 'Ax_u', 'l_Ax'}, tags(1:nargout)))
                    g = cell(1,2);
                    [g{:}] = om.eval_lin_constraint(s.x, name, idx);
                end
                for k = 1:nargout
                    switch tags{k}
                        case 'g'
                            varargout{k} = g;
                        case 'Ax_u'
                            varargout{k} = g{1};
                        case 'l_Ax'
                            varargout{k} = g{2};
                        case 'mu_l'
                            varargout{k} = s.lambda.mu_l(i1:iN);
                        case 'mu_u'
                            varargout{k} = s.lambda.mu_u(i1:iN);
                        otherwise
                            error('opt_model/get_soln: unknown tag ''%s''', tags{k});
                    end
                end
            end
        case 'nle'
            if ismember('dg', tags(1:nargout))
                [g, dg] = om.eval_nln_constraint(s.x, 1, name, idx);
            elseif ismember('g', tags(1:nargout))
                g = om.eval_nln_constraint(s.x, 1, name, idx);
            end
            for k = 1:nargout
                switch tags{k}
                    case 'g'
                        varargout{k} = g;
                    case 'dg'
                        varargout{k} = dg;
                    case 'lam'
                        varargout{k} = s.lambda.eqnonlin(i1:iN);
                    otherwise
                        error('opt_model/get_soln: unknown tag ''%s''', tags{k});
                end
            end
        case 'nli'
            if ismember('dh', tags(1:nargout))
                [h, dh] = om.eval_nln_constraint(s.x, 0, name, idx);
            elseif ismember('h', tags(1:nargout))
                h = om.eval_nln_constraint(s.x, 0, name, idx);
            end
            for k = 1:nargout
                switch tags{k}
                    case 'h'
                        varargout{k} = h;
                    case 'dh'
                        varargout{k} = dh;
                    case 'mu'
                        varargout{k} = s.lambda.ineqnonlin(i1:iN);
                    otherwise
                        error('opt_model/get_soln: unknown tag ''%s''', tags{k});
                end
            end
        case 'nlc'
            if ismember('d2f', tags(1:nargout))
                [f, df, d2f] = om.eval_nln_cost(s.x, name, idx);
            elseif ismember('df', tags(1:nargout))
                [f, df] = om.eval_nln_cost(s.x, name, idx);
            else
                f = om.eval_nln_cost(s.x, name, idx);
            end
            for k = 1:nargout
                switch tags{k}
                    case 'f'
                        varargout{k} = f;
                    case 'df'
                        varargout{k} = df;
                    case 'd2f'
                        varargout{k} = d2f;
                    otherwise
                        error('opt_model/get_soln: unknown tag ''%s''', tags{k});
                end
            end
        case 'qdc'
            if ismember('d2f', tags(1:nargout))
                [f, df, d2f] = om.eval_quad_cost(s.x, name, idx);
            elseif ismember('df', tags(1:nargout))
                [f, df] = om.eval_quad_cost(s.x, name, idx);
            else
                f = om.eval_quad_cost(s.x, name, idx);
            end
            for k = 1:nargout
                switch tags{k}
                    case 'f'
                        varargout{k} = f;
                    case 'df'
                        varargout{k} = df;
                    case 'd2f'
                        varargout{k} = d2f;
                    otherwise
                        error('opt_model/get_soln: unknown tag ''%s''', tags{k});
                end
            end
        otherwise
            error('opt_model/get_soln: unknown set_type ''%s''', set_type);
    end     %% switch set_type
    end
end     %% if N

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function K = anisotropic_rbf(kern,dat1,dat2,ind1,ind2,kerParam),

% K = rbf(d1,d2,ind1,ind2,param), compute the kernel
%     matrix between d1 and d2
% for a rbf kernel exp(-||x-z||^2/(2*param^2))
%     where x is from d1 and z from d2

w = [];
for i = 1:length(kerParam)-1
    w = [w,kerParam{i}];
end
v = kerParam{end};

X2 = get_x(dat2,ind2);
X1 = get_x(dat1,ind1);

for i  =1:size(X1,1)
    for j = 1:size(X2,1)
         s = sum(w.*((X1(i,:) - X2(j,:)).^2));
         K(j,i) = v*exp(-0.5*s);
     end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% geopdes_inv: the inverse matrix or, for codimension > 0 (rdim > ndim)
% the Moore-Penrose pseudo-inverse.
%
% To be used with msh.geo_map_jac, to compute DF^{-1}.
% In the case of codimension > 0, it computes
%   J = (DF^t * DF)^{-1} * DF^t
%   det = sqrt ( det(DF^t * DF)^{-1} )
%
%  [JinvT, det] = geopdes_invT__ (geo_map_jac)
%
% OUTPUT:
%   Jinv: the computed matrix evaluated at every point. Size (ndim x rdim x nqn x nel)
%   det:  the determinant evaluated at every point. Size (nqn x nel)
%
% Copyright (C) 2010 Carlo de Falco
% Copyright (C) 2015 Rafael Vazquez
%
% 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 Octave; see the file COPYING.  If not, see
% <http://www.gnu.org/licenses/>.

function [Jinv, det] = geopdes_inv__ (v)

  [JinvT, det] = geopdes_invT__ (v);
  Jinv = permute (JinvT, [2, 1, 3:numel(size(v))]);

end

%!test
%! A = [1 2; 3 4];
%! Ainv = geopdes_inv__ (A);
%! assert (Ainv * A, eye(2), 1e-14)
%!
%!test
%! A = [1 2 3; 4 5 6; 7 8 10];
%! Ainv = geopdes_inv__ (A);
%! assert (Ainv * A, eye(3), 1e-14)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function r8poly_lagrange_0_test ( )

%*****************************************************************************80
%
%% R8POLY_LAGRANGE_0_TEST tests R8POLY_LAGRANGE_0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  npol = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8POLY_LAGRANGE_0_TEST\n' );
  fprintf ( 1, '  R8POLY_LAGRANGE_0 evaluates the Lagrange\n' );
  fprintf ( 1, '  factor W(X) at a point.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of data points is %d\n', npol );
%
%  Set the abscissas of the polynomials.
%
  xlo = 0.0;
  xhi = npol - 1;

  xpol = r8vec_even ( npol, xlo, xhi );

  r8vec_print ( npol, xpol, '  Abscissas:' );
%
%  Evaluate W(X).
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X          W(X)\n' );
  fprintf ( 1, '\n' );

  nx = 4 * npol - 1;

  for ival = 1 : nx

    xval = r8vec_even_select ( nx, xlo, xhi, ival );

    w = r8poly_lagrange_0 ( npol, xpol, xval );

    fprintf ( 1, '%12f  %12e\n', xval, w );

  end

  return
end


text:
function prediction = NIN(varargin)
%NIN Returns a Network-in-Network model for CIFAR10
%   M = models.NIN() returns the model proposed in:
%
%     Lin et al, "Network in network", arXiv technical report 2013.
%
%   models.NIN(..., 'option', value, ...) accepts the following options:
%
%   `input`:: default input
%     Specifies an input (images) layer for the network. If unspecified, a
%     new one is created.
%
%   `numClasses`:: 10
%     Number of output classes.
%
%   `batchNorm`:: true
%     Whether to use batch normalization.
%
%   Any other options will be passed to models.ConvBlock(), and can be used
%   to change the activation function, weight initialization, etc.
%
%   Suggested SGD training options are also returned in the struct M.meta.

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

  % parse options. unknown arguments will be passed to ConvBlock (e.g.
  % batchNorm).
  opts.input = Input('name', 'images', 'gpu', true) ;  % default input layer
  opts.numClasses = 10 ;  % number of predicted classes
  opts.batchNorm = true ;  % whether to use batch normalization
  [opts, convBlockArgs] = vl_argparse(opts, varargin, 'nonrecursive') ;

  % build network
  images = opts.input ;

  % first NIN block
  channels = [192 160 96] ;  % number of outputs channels per conv layer
  ker = [5 5] ;  % conv kernel
  poolKer = [3 3] ;  % pooling kernel
  poolMethod = 'max' ;  % pooling method
  pad = 2 ;  % input padding
  m1 = ninBlock(images, 3, channels, ker, pad, ...
    poolKer, poolMethod, convBlockArgs, opts.batchNorm, false) ;
  outChannels = channels(3) ;  % output channels of the NIN block

  % second NIN block
  channels = [192 192 192] ;
  ker = [5 5] ;
  poolKer = [3 3] ;
  poolMethod = 'avg' ;
  pad = 2 ;
  m2 = ninBlock(m1, outChannels, channels, ker, pad, ...
    poolKer, poolMethod, convBlockArgs, opts.batchNorm, false) ;
  outChannels = channels(3) ;

  % third NIN block
  channels = [192 192 opts.numClasses] ;
  ker = [3 3] ;
  poolKer = [7 7] ;
  poolMethod = 'avg' ;
  pad = 1 ;
  prediction = ninBlock(m2, outChannels, channels, ker, pad, ...
    poolKer, poolMethod, convBlockArgs, opts.batchNorm, true) ;

  % default training options for this network
  defaults.batchSize = 100 ;
  defaults.weightDecay = 0.0005 ;
  defaults.imageSize = [32, 32, 3] ;
  % the default learning rate schedule
  if ~opts.batchNorm
    defaults.learningRate = [0.002, 0.01, 0.02, 0.04 * ones(1,80), 0.004 * ones(1,10), 0.0004 * ones(1,10)] ;
    defaults.numEpochs = numel(defaults.learningRate) ;
  else
    defaults.learningRate = 0.01 ;
    defaults.numEpochs = 40 ;
  end
  prediction.meta = defaults ;

end

function block = ninBlock(in, inChannels, outChannels, ...
  ker, pad, poolKer, poolMethod, convBlockArgs, bnorm, final)

  % get conv block generator with the given options.
  % default activation is ReLU, with post-activation batch normalization.
  conv = models.ConvBlock('batchNorm', bnorm, convBlockArgs{:}) ;

  % if it's the final layer (prediction), the last conv block has no
  % activation or batch-norm
  finalArgs = {} ;
  if final
    finalArgs = {'activation', 'none', 'batchNorm', false} ;
  end

  % 3 conv blocks
  c1 = conv(in, 'size', [ker(1:2), inChannels, outChannels(1)], 'pad', pad) ;
  c2 = conv(c1, 'size', [1, 1, outChannels(1), outChannels(2)]) ;
  c3 = conv(c2, 'size', [1, 1, outChannels(2), outChannels(3)], finalArgs{:}) ;

  % pooling
  p1 = vl_nnpool(c3, poolKer, 'method', poolMethod, 'stride', 2) ;

  % dropout, skipped if it's the final layer (prediction)
  if ~final
    block = vl_nndropout(p1, 'rate', 0.5) ;
  else
    block = p1 ;
  end
end

repo: MFT


text:
function status = cbarSave(cbar, pth);
% Save a colorbar to a MATLAB file.
%
% status = cbarSave(cbar, [pth=dialog]);
%
% Returns a 1 if the save was successful and a 0 otherwise.
%
% ras, 02/24/07.
status = 0;

if notDefined('cbar'), error('Need a colorbar.'); end

if notDefined('pth') % dialog
	pth = mrvSelectFile('w', 'mat', [], 'Save a colorbar...');
end

save(pth, '-struct', 'cbar');

if exist(pth, 'file')
	status = 1;
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


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

    methods
        function obj = Ty(varargin)
            obj@ETS2('Ty', varargin{:});

        end

        % function s = char(ets)
        %     s = ets.what;
        % end
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function [s] = find_first_path(guesses,quiet)
  si = find(cellfun(@(guess) exist(guess,'file'),guesses),1,'first');
  if isempty(si)
    if nargin>1 && quiet
      s = [];
      return;
    else
      error('Could not find path');
    end
  end
  s = guesses{si};
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [nearest_neighbours] = find_nearest_neighbours( database, desc, dist_ratio, max_dist )

% [nearest_neighbours] = find_nearest_neighbours( database, desc, max_dist )
%
% Find the indices of the nearest neighbours of the given desriptors in the
% specified database.  Uses euclidean distance.
%
% Input:
% database - descriptor database created by add_descriptors_to_database.
% desc - descriptors from the SIFT function.
% dist_ratio - maximum ratio between distances of nearest and second closest
%   neighbour for a match to be allowed.
%
% Output:
% nearest_neighbours - indices of the nearest neighbours for the descriptors
%   (descriptors with no neighbour closer than max_dist will have index 0).
%
% Thomas F. El-Maraghi
% May 2004
%%
%if ~exist( 'dist_ratio' )
%   dist_ratio = 0.8;
%end
%%
nearest_neighbours = zeros(size(desc,1),1);
for k = 1:size(desc,1)

   dist = sqrt(sum((database.desc - repmat(desc(k,:),size(database.desc,1),1)).^2,2));  % thid is a column vector of euclidean distances (one element for each database descriptor)
   [nn1_dist idx] = min(dist);
   dist(idx) = max(dist);
   nn2_dist = min(dist);
   if nn1_dist/nn2_dist >= dist_ratio || nn1_dist>max_dist %%%%% I added the last bit (&& ...)
      idx = 0;
   end
   nearest_neighbours(k,1:2) = [idx nn1_dist];
end


text:
function res = imCropBorder(img, pad, varargin)
%IMCROPBORDER Crop the border around a 2D or 3D image
%
%   IMG = imCropBorder(IMG, BORDER)
%   Crop the image by removing BORDER pixels around the image. The size of
%   new image is given by size(IMG)-2*BORDER.
%   The function works for 2D or 3D images, grayscale or color.
%
%   IMG = imCropBorder(IMG, [BORDER1 BORDER2])
%   IMG = imCropBorder(IMG, [BORDER1 BORDER2 BORDER3])
%   Specifies different sizes for borders in first and second directions
%   (and third direction for 3D images)
%
%   IMG = imCropBorder(IMG, [BORDER10 BORDER11 BORDER20 BORDER21])
%   IMG = imCropBorder(IMG, [BORDER10 BORDER11 BORDER20 BORDER21 BORDER30 BORDER31])
%   Specifies different sizes for borders in first and second directions,
%   and in the beginning and at the end in each direction.
%
%
%   Example
%   % Add 20 pixel in each side of rice image
%     img = imread('rice.png');
%     img2 = imAddBorder(img, 20);
%     imshow(img2);
%
%   % add white borders with different sizes
%     img = imread('rice.png');
%     img2 = imAddBorder(img, [10 20 30 40], 255);
%     imshow(img2);
%
%   % add cyan border around a color image
%     img = imread('peppers.png');
%     img2 = imAddBorder(img, [30 20], [0 1 1]);
%     imshow(img2);
%
%   See also
%     imCrop, imAddBorder
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2018-11-05,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

%   History

%% Initialisations

% size of input image
dim = size(img);
isColor = size(img, 3) == 3;
isPlanar = length(dim) == 2 || (length(dim) == 3 && isColor);

if isPlanar
    %% Planar case

    % convert pad to have 4 values
    if length(pad) == 1
        pad = [pad pad pad pad];
    elseif length(pad)==2
        pad = [pad(1) pad(1) pad(2) pad(2)];
    end

    dim1 = (pad(1)+1):(size(img,1)-pad(2));
    dim2 = (pad(3)+1):(size(img,2)-pad(4));

    % fillup result image with initial image
    res = img(dim1, dim2, :);

elseif length(dim) - isColor == 3
    %% Case of 3D image

    % convert pad to have 6 values
    if length(pad) == 1
        pad = [pad pad pad pad pad pad];
    elseif length(pad) == 3
        pad = [pad(1) pad(1) pad(2) pad(2) pad(3) pad(3)];
    end

    % compute index of voxels to keep along each dimension
    dim1 = (pad(1)+1):(size(img,1)-pad(2));
    dim2 = (pad(3)+1):(size(img,2)-pad(4));
    dim3 = (pad(5)+1):(size(img,3+isColor)-pad(6));

    % select voxels to keep
    if isColor
        res = img(dim1, dim2, :, dim3);
    else
        res = img(dim1, dim2, dim3);
    end

else
    error('Image dimension not managed');
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [biomassValues, targetValues] = runProductionEnvelope(model, targetRxn, biomassRxn, nPts)
% runProductionEnvelope
%   Calculates the byproduct secretion envelope
%
% Input:
%   model            a model structure
%   targetRxn        identifier of target metabolite production reaction
%   biomassRxn       identifier of biomass reaction
%   nPts             number of points in the plot (opt, default 20)
%
% Output:
%   biomassValues    Biomass values for plotting
%   targetValues     Target upper and lower bounds for plotting
%
% Modified from COBRA Toolbox productionEnvelope.m
%
% Usage: [biomassValues, targetValues] = runProductionEnvelope(model, targetRxn, biomassRxn, nPts)

if nargin < 4
    nPts = 20;
end

% Run FBA to get upper bound for biomass
model   = setParam(model,'obj',biomassRxn,1);
[solMax,hsSol]  = solveLP(model);
solMax  = solMax.x(logical(model.c));
model.c = -model.c;
solMin  = solveLP(model,0,[],hsSol);
solMin  = solMin.x(logical(model.c));

% Create biomass range vector
biomassValues = linspace(solMin,solMax,nPts);
targetUpperBound = nan(1,numel(biomassValues));
targetLowerBound = nan(1,numel(biomassValues));

fprintf('Creating production envelope...   0%% complete');
% Max/min for target production
model = setParam(model,'obj',targetRxn,1);
for i = 1:length(biomassValues)
    progress=pad(num2str(floor(i/numel(biomassValues)*100)),3,'left');
    fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b%s%% complete',progress);
    model1 = setParam(model,'eq',biomassRxn,biomassValues(i));
    sol = solveLP(model1,0,[],hsSol);
    if (sol.stat > 0)
        targetUpperBound(i) = sol.x(logical(model.c));
    else
        targetUpperBound(i) = NaN;
    end
    model1.c = -model1.c;
    sol = solveLP(model1,0,[],hsSol);
    if (sol.stat > 0)
        targetLowerBound(i) = sol.x(logical(model1.c));
    else
        targetLowerBound(i) = NaN;
    end
end
fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\bCOMPLETE\n');

% Plot results
plot([biomassValues fliplr(biomassValues)],[targetUpperBound fliplr(targetLowerBound)],'blue','LineWidth',2);
axis tight;
ylabel([strrep(targetRxn,'_','-') ' (mmol/gDW h)']);
xlabel('Growth rate (1/h)');
targetValues = [targetLowerBound; targetUpperBound];
end

author: SysBioChalmers

repo: RAVEN


%train
%dataset_dir = '/home/gyeongsikmoon/workspace/Data/Hand_pose_estimation/NYU/dataset/train/';
%save_dir = '/home/gyeongsikmoon/workspace/Data/Hand_pose_estimation/NYU/dataset/train/parsed/';
%tot_frame_num = 72757;

%test
dataset_dir = '/home/gyeongsikmoon/workspace/Data/Hand_pose_estimation/NYU/dataset/test/';
save_dir = '/home/gyeongsikmoon/workspace/Data/Hand_pose_estimation/NYU/dataset/test/parsed/';
tot_frame_num = 8252;

kinect_index = 1;

for image_index = 1:tot_frame_num
    filename_prefix = sprintf('%d_%07d', kinect_index, image_index);

    if exist([dataset_dir, 'depth_', filename_prefix, '.png'], 'file')

        %% Load and display a depth example
        % The top 8 bits of depth are packed into green and the lower 8 bits into blue.
        depth = imread([dataset_dir, 'depth_', filename_prefix, '.png']);
        depth = uint16(depth(:,:,3)) + bitsll(uint16(depth(:,:,2)), 8);

        fp_save = fopen([save_dir, 'depth_', filename_prefix, '.bin'],'w');
        fwrite(fp_save,permute(depth,[2,1,3]),'float');
        fclose(fp_save);

        %delete(strcat(folderpath,img_name));
    end

end


text:
function varargout = real(varargin)
%REAL      Real part of a CHEBFUN2.
%
% See also IMAG.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

[varargout{1:nargout}] = real@separableApprox(varargin{:});

end


text:
function line_enhanced = EnhanceLine(band)
%ENHANCELINE enhance NDBI to light urban/built-up areas and dark other
%bright surface, such as desert, rock. ref Guindon et. RSE 2004

% also can see the details at
% https://homepages.inf.ed.ac.uk/rbf/HIPR2/linedet.htm.
%     band=data_toabt.BandGreen;
%
%     band=ndbi;

%% line with a length of three pixels
%     template1=[-1 1 0;
%                -1 1 0;
%                -1 1 0;];
%     template2=[0 1 -1;
%                0 1 -1;
%                0 1 -1;];
%     line_enhanced1 = imfilter(band,template1);
%     line_enhanced2 = imfilter(band,template2);
%     line_enhanced=max(line_enhanced1,line_enhanced2);
%
%     template1=[-1 -1 -1;
%                 1  1  1;
%                 0  0  0;];
%     template2=[0  0  0;
%                1  1  1;
%               -1 -1 -1;];
%     line_enhanced1 = imfilter(band,template1);
%     line_enhanced2 = imfilter(band,template2);
%     line_enhanced=max(line_enhanced1,line_enhanced);
%     line_enhanced=max(line_enhanced2,line_enhanced);
%
%     template1=[1 -1 -1;
%                0  1 -1;
%                0  0  1;];
%     template2=[1  0  0;
%               -1  1  0;
%               -1 -1  1;];
%     line_enhanced1 = imfilter(band,template1);
%     line_enhanced2 = imfilter(band,template2);
%     line_enhanced=max(line_enhanced1,line_enhanced);
%     line_enhanced=max(line_enhanced2,line_enhanced);
%
%     template1=[-1 -1  1;
%                -1  1  0;
%                 1  0  0;];
%     template2=[0  0  1;
%                0  1 -1;
%                1 -1 -1;];
%
%     line_enhanced1 = imfilter(band,template1);
%     line_enhanced2 = imfilter(band,template2);
%     line_enhanced=max(line_enhanced1,line_enhanced);
%     line_enhanced=max(line_enhanced2,line_enhanced);
% 	line_enhanced = line_enhanced./3;

    template=[-1 2 -1;
              -1 2 -1;
              -1 2 -1;];
	template = template./6;
    line_enhanced = imfilter(band,template);

    template=[-1 -1 -1;
               2  2  2;
              -1 -1 -1;];
	template = template./6;
    line_enhanced_new = imfilter(band,template);
    line_enhanced=max(line_enhanced_new,line_enhanced);

    template =[2 -1 -1;
               -1  2 -1;
               -1  -1  2;];
	template = template./6;
    line_enhanced_new = imfilter(band,template);
    line_enhanced=max(line_enhanced_new,line_enhanced);

    template =[-1  -1  2;
               -1  2 -1;
               2 -1 -1;];
	template = template./6;
    line_enhanced_new = imfilter(band,template);
    line_enhanced=max(line_enhanced_new,line_enhanced);
end

author: GERSL

repo: Fmask

save_path: github-repos/MATLAB/GERSL-Fmask


text:
function e = calcError(odf1,odf2,varargin)
% calculate approximation error between two ODFs
%
% Syntax
%   e = calcError(odf1,odf2)
%   e = calcError(odf,pf,'RP')
%
% Input
%  odf1, odf2 - @SO3Fun
%  pf   - @PoleFigure
%  S3G  - @SO3Grid of quadrature nodes (optional)
%
% Options
%  L0 - measure of the orientation space where $|odf1 -- odf2|>\epsilon|
%  L1 - L^1 error
%  L2 - L^2 error
%  RP - RP  error (default)
%  resolution - resolution used for calculation of the error
%
% See also
% PoleFigure/calcODF PoleFigure/calcError

% compare with a pole figure
if isa(odf2,'PoleFigure'), e = calcError(odf2,odf1,varargin{:}); return;end

evaluated = ~check_option([{odf2} varargin],'evaluated');

% compare two odfs
if evaluated
  assert(odf1.CS.Laue == odf2.CS.Laue && odf1.SS.Laue == odf2.SS.Laue,...
    'Input ODFs does not have same symmetry.');
end

% TODO
% Fourier based algorithm
if check_option(varargin,'L2')

  e = norm(odf1 - odf2)./norm(odf2);

% quadrature  rule based algorithm
else

  % get approximation grid
  S3G = extract_SO3grid(odf1,'resolution',5*degree,varargin{:});

  % eval ODFs
  if evaluated %second ODF allready evaluated
    d1 = eval(odf2,orientation(S3G),varargin{:});
  else
    d1 = get_option([{odf2} varargin],'evaluated',[],'double');
  end

  if isa(odf2,'double')
    d2 = odf1;
  else
    d2 = eval(odf1,orientation(S3G),varargin{:});
  end

  % calculate the error
  if check_option(varargin,'L0')
    epsilon = get_option(varargin,'L0',1);
    for i = 1:length(epsilon)
      e(i) = sum(abs(d1(:)-d2(:)) > epsilon(i))/numel(d1);
    end
  elseif check_option(varargin,'L2')
    e = norm(d1(:)-d2(:)) / norm(d2(:));
  else
    e = sum(abs(d1(:)-d2(:))) / numel(d1) /2;
  end
end

author: mtex-toolbox

repo: mtex

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


text:
function stat = test_ReadYaml()
% this function tests reading in the yaml file

stat.ok = 1;
stat.desc = '';
try
    %stat.test_ReadYaml_SimpleStructure = test_ReadYaml_SimpleStructure();
    %stat.test_ReadYaml_DateTime = test_ReadYaml_DateTime();
    fprintf('Testing read ');
    stat.test_RY_Matrices = test_RY_Matrices();
    fprintf('.');
    stat.test_RY_Whitespaces = test_RY_Whitespaces();
    fprintf('.');
    stat.test_RY_FloatingPoints = test_RY_FloatingPoints();
    fprintf('.');
    stat.test_RY_Indentation = test_RY_Indentation();
    fprintf('.');
    stat.test_RY_SequenceMapping = test_RY_SequenceMapping();
    fprintf('.');
    stat.test_RY_Simple = test_RY_Simple();
    fprintf('.');
    stat.test_RY_Time = test_RY_Time();
    fprintf('.');
    stat.test_RY_TimeVariants = test_RY_TimeVariants();
    fprintf('.');
    stat.test_RY_Import = test_RY_Import();
    fprintf('.');
    stat.test_RY_ImportDef = test_RY_ImportDef();
    fprintf('.');
    stat.test_RY_ImportNonex = test_RY_ImportNonex();
    fprintf('.');
    stat.test_RY_Inheritance = test_RY_Inheritance();
    fprintf('.');
    stat.test_RY_InheritanceMultiple = test_RY_InheritanceMultiple();
    fprintf('.');
    stat.test_RY_InheritanceLoop = test_RY_InheritanceLoop();
    fprintf('.');
    stat.test_RY_usecase_01 = test_RY_usecase_01();
    fprintf('.\n');

catch
    stat.ok = 0;
    stat.desc  = 'Program crash';
end

end

function result = PTH_PRIMITIVES()
    result = sprintf('Data%stest_primitives%s',filesep,filesep);
end

function result = PTH_IMPORT()
    result = sprintf('Data%stest_import%s',filesep,filesep);
end

function result = PTH_INHERITANCE()
    result = sprintf('Data%stest_inheritance%s',filesep,filesep);
end

function stat = test_RY_Matrices()
    stat.ok = 1;
    stat.desc = '';
    try
        ry = ReadYaml([PTH_PRIMITIVES() 'matrices.yaml']);
        tv = load([PTH_PRIMITIVES() 'matrices.mat']);
        if ~isequal(ry, tv.testval)
            stat.desc  = 'Wrong values loaded';
            stat.ok = 0;
        end;
    catch
        stat.ok = 0;
        stat.desc = 'Crash';
    end;
end

function stat = test_RY_FloatingPoints()
    stat.ok = 1;
    stat.desc = '';
    try
        ry = ReadYaml([PTH_PRIMITIVES() 'floating_points.yaml']);
        tv = load([PTH_PRIMITIVES() 'floating_points.mat']);
        if ~isequalwithequalnans(ry, tv.testval)
            stat.desc  = 'Wrong values loaded';
            stat.ok = 0;
        end;
    catch
        stat.ok = 0;
        stat.desc = 'Crash';
    end;
end

function stat = test_RY_Indentation()
    stat.ok = 1;
    stat.desc = '';
    try
        ry = ReadYaml([PTH_PRIMITIVES() 'indentation.yaml']);
        tv = load([PTH_PRIMITIVES() 'indentation.mat']);
        if ~isequal(ry, tv.testval)
            stat.desc  = 'Wrong values loaded';
            stat.ok = 0;
        end;
    catch
        stat.ok = 0;
        stat.desc = 'Crash';
    end;
end

function stat = test_RY_SequenceMapping()
    stat.ok = 1;
    stat.desc = '';
    try
        ry = ReadYaml([PTH_PRIMITIVES() 'sequence_mapping.yaml']);
        tv = load([PTH_PRIMITIVES() 'sequence_mapping.mat']);
        if ~isequal(ry, tv.testval)
            stat.desc  = 'Wrong values loaded';
            stat.ok = 0;
        end;
    catch
        stat.ok = 0;
        stat.desc = 'Crash';
    end;
end

function stat = test_RY_Simple()
    stat.ok = 1;
    stat.desc = '';
    try
        ry = ReadYaml([PTH_PRIMITIVES() 'simple.yaml']);
        tv = load([PTH_PRIMITIVES() 'simple.mat']);
        if ~isequal(ry, tv.testval)
            stat.desc  = 'Wrong values loaded';
            stat.ok = 0;
        end;
    catch
        stat.ok = 0;
        stat.desc = 'Crash';
    end;
end

function stat = test_RY_Time()
    stat.ok = 1;
    stat.desc = '';
    try
        ry = ReadYaml([PTH_PRIMITIVES() 'time.yaml']);
        tv = load([PTH_PRIMITIVES() 'time.mat']);
        if ~isequal(ry, tv.testval)
            stat.desc  = 'Wrong values loaded';
            stat.ok = 0;
        end;
    catch
        stat.ok = 0;
        stat.desc = 'Crash';
    end;
end

function stat = test_RY_TimeVariants()
    stat.ok = 1;
    stat.desc = '';
    try
        ry = ReadYaml([PTH_PRIMITIVES() 'time_variants.yaml']);
        tv = load([PTH_PRIMITIVES() 'time_variants.mat']);
        if ~isequal(ry, tv.testval)
            stat.desc  = 'Wrong values loaded';
            stat.ok = 0;
        end;
    catch
        stat.ok = 0;
        stat.desc = 'Crash';
    end;
end

function stat = test_RY_Import()
    stat.ok = 1;
    stat.desc = '';
    try
        ry = ReadYaml([PTH_IMPORT() 'import.yaml']);
        tv = load([PTH_IMPORT() 'import.mat']);
        if ~isequal(ry, tv.testval)
            stat.desc  = 'Wrong values loaded';
            stat.ok = 0;
        end;
    catch
        stat.ok = 0;
        stat.desc = 'Crash';
    end;
end

function stat = test_RY_ImportDef()
    stat.ok = 1;
    stat.desc = '';
    try
        ry = ReadYaml([PTH_IMPORT() 'import_def.yaml']);
        tv = load([PTH_IMPORT() 'import_def.mat']);
        if ~isequal(ry, tv.testval)
            stat.desc  = 'Wrong values loaded';
            stat.ok = 0;
        end;
    catch
        stat.ok = 0;
        stat.desc = 'Crash';
    end;
end

function stat = test_RY_ImportNonex()
    stat.ok = 0;
    stat.desc = 'Did not end with any exception.';
    try
        try
            ry = ReadYaml([PTH_IMPORT() 'import_nonex.yaml'],1);
        catch ex
            if strcmp(ex.identifier, 'MATLAB:MATYAML:FileNotFound')
                stat.desc = '';
                stat.ok = 1;
            else
                rethrow(ex);
            end;
        end;

    catch
        stat.ok = 0;
        stat.desc = 'Crash';
    end;
end

function stat = test_RY_Inheritance()
    stat.ok = 1;
    stat.desc = '';
    try
        ry = ReadYaml([PTH_INHERITANCE() 'inheritance.yaml']);
        tv = load([PTH_INHERITANCE() 'inheritance.mat']);
        if ~isequal(ry, tv.testval)
            stat.desc  = 'Wrong values loaded';
            stat.ok = 0;
        end;
    catch
        stat.ok = 0;
        stat.desc = 'Crash';
    end;
end

function stat = test_RY_InheritanceMultiple()
    stat.ok = 1;
    stat.desc = '';
    try
        ry = ReadYaml([PTH_INHERITANCE() 'inheritance_multiple.yaml']);
        tv = load([PTH_INHERITANCE() 'inheritance_multiple.mat']);
        if ~isequal(ry, tv.testval)
            stat.desc  = 'Wrong values loaded';
            stat.ok = 0;
        end;
    catch
        stat.ok = 0;
        stat.desc = 'Crash';
    end;
end

function stat = test_RY_InheritanceLoop()
    stat.ok = 0;
    stat.desc = 'Did not end with any exception.';
    try
        try
            ry = ReadYaml([PTH_INHERITANCE() 'inheritance_loop.yaml']);
        catch ex
            if strcmp(ex.identifier, 'MATLAB:MATYAML:inheritedtwice')
                stat.desc = '';
                stat.ok = 1;
            else
                rethrow(ex);
            end;
        end;

    catch
        stat.ok = 0;
        stat.desc = 'Crash';
    end;
end

function stat = test_RY_Whitespaces()
    stat.ok = 1;
    stat.desc = '';
    try
        ry = ReadYaml([PTH_PRIMITIVES() 'whitespaces.yaml']);
        if ~isfield(ry,'ImageFile') || ~isfield(ry,'ContoursCount')
            stat.desc  = 'Wrong values loaded';
            stat.ok = 0;
        end;
    catch
        stat.ok = 0;
        stat.desc = 'Crash';
    end;
end

function stat = test_RY_usecase_01()
    stat.ok = 1;
    stat.desc = '';
    try
        ry = ReadYaml([PTH_PRIMITIVES() 'usecase_struct_01.yaml']);
        tv = load([PTH_PRIMITIVES() 'usecase_struct_01.mat']);
        if ~isequalwithequalnans(ry, tv.testval)
            stat.desc  = 'Wrong values loaded';
            stat.ok = 0;
        end;
    catch
        stat.ok = 0;
        stat.desc = 'Crash';
    end;
end

function stat = test_ReadYaml_SimpleStructure()

stat.ok = 1;
stat.desc = '';
try
    s = ReadYaml('simple.yaml');

    ages = [s.age];

    if not(isequal([33 27], ages))  || not(all(ismember({'John Smith', 'Mary Smith'}, {s.name}) ))
        stat.desc  = ' Wrong values loaded';
        stat.ok = 0;
    end

catch
       stat.desc  = 'Program crash';
       stat.ok = 0;
end

end

function stat = test_ReadYaml_DateTime()

stat.ok = 1;
stat.desc = '';
try
    s = ReadYaml('time.yaml');

    if ~isa(s.Data.B1_S_SW{1},'DateTime')
        stat.desc  = ' Wrong data type of datetick';
        stat.ok = 0;
    end
    if isa(s.Data.B1_S_SW{2},'DateTime')
        stat.desc  = ' Wrong data type of datetick';
        stat.ok = 0;
    end
catch
       stat.desc  = 'Program crash';
       stat.ok = 0;
end
end

author: ziyuw

repo: rembo

save_path: github-repos/MATLAB/ziyuw-rembo


text:
classdef ipfSpotKey < ipfColorKey
%
% Colorizes single spots in the inverse pole figure with individual colors.
%
% Syntax
%   ipfKey = ipfSpotKey
%   ipfKey.inversePoleFigureDirection = zvector;
%   ipfKey.center = Miller(1,0,0,cs); % the centers of the spots in the inverse pole figure
%   ipfKey.color = [0 0 1];           % the color of the spots
%   ipfKey.psi = S2DeLaValleePoussinKernel('halfwidth',7.5*degree);
%
%   color = ipfKey.orientation2color(ori)
%
% Class Properties
%  center  - list of crystal directions @Miller
%  color   - n-by-3 list representing RGB values, one for each center
%  psi     - @S2Kernel providing the width and brightness for colored fibre
%  inversePoleFigureDirection - specimen direction @vector3d
%
% See also
% EBSDAdvancedMapping

  properties
    center % list of crystal directions @Miller
    color  % list of RGB values, one for each center
    psi    % @2Kernel providing the width and brightness for colored fibre
  end

  methods
    function oM = ipfSpotKey(varargin)
      oM = oM@ipfColorKey(varargin{:});

      oM.center = get_option(varargin,'center',Miller(0,0,1,oM.CS1));
      oM.CS1 = oM.center.CS;
      oM.color = get_option(varargin,'color',[1 0 0]);
      oM.psi = get_option(varargin,'S2Kernel',...
        S2DeLaValleePoussinKernel('halfwidth',get_option(varargin,'halfwidth',10*degree)));

      oM.dirMap = directionColorKey(oM.CS1,'dir2color',@(varargin) oM.dir2color(varargin{:}));

    end

    function rgb = dir2color(oM,h)

      h = normalize(h);
      rgb = ones(length(h),3);

      for k=1:length(oM.center)

        w = oM.psi.eval(dot(h,normalize(oM.center(k)))) ./ oM.psi.eval(1);

        if ~any(oM.color(k,:))% fix in case of black
          cdata = repmat([0 1 1],length(h),1);
          cdata(:,2) = w(:).*cdata(:,2);
          cdata = reshape(hsv2rgb(cdata),[],3);
          cdata(:,1)=cdata(:,2);
          rgb = rgb.*cdata;
        else
          cdata = rgb2hsv(repmat(oM.color(k,:),length(h),1));
          cdata(:,2) = w(:).*cdata(:,2);
          cdata(:,3) = 1;
          cdata = reshape(hsv2rgb(cdata),[],3);
          rgb = rgb.*cdata;
        end
      end
    end
  end
end


text:
function [LD,p,q] = ldlchol (A,beta)					    %#ok
%LDLCHOL sparse A=LDL' factorization
%   Note that L*L' (LCHOL) and L*D*L' (LDLCHOL) factorizations are faster than
%   R'*R (CHOL2 and CHOL) and use less memory.  The LL' and LDL' factorization
%   methods use tril(A).  A must be sparse.
%
%   Example:
%   LD = ldlchol (A)            return the LDL' factorization of A
%   [LD,p] = ldlchol (A)        similar [R,p] = chol(A), but for L*D*L'
%   [LD,p,q] = ldlchol (A)      factorizes A(q,q) into L*D*L', where q is a
%                               fill-reducing ordering
%
%   LD = ldlchol (A,beta)       return the LDL' factorization of A*A'+beta*I
%   [LD,p] = ldlchol (A,beta)   like [R,p] = chol(A*A'+beta+I)
%   [LD,p,q] = ldlchol (A,beta) factorizes A(q,:)*A(q,:)'+beta*I into L*D*L'
%
%   The output matrix LD contains both L and D.  D is on the diagonal of LD, and
%   L is contained in the strictly lower triangular part of LD.  The unit-
%   diagonal of L is not stored.  You can obtain the L and D matrices with
%   [L,D] = ldlsplit (LD).  LD is in the form needed by ldlupdate.
%
%   Explicit zeros may appear in the LD matrix.  The pattern of LD matches the
%   pattern of L as computed by symbfact2, even if some entries in LD are
%   explicitly zero.  This is to ensure that ldlupdate and ldlsolve work
%   properly.  You must NOT modify LD in MATLAB itself and then use ldlupdate
%   or ldlsolve if LD contains explicit zero entries; ldlupdate and ldlsolve
%   will fail catastrophically in this case.
%
%   You MAY modify LD in MATLAB if you do not pass it back to ldlupdate or
%   ldlsolve.  Just be aware that LD contains explicit zero entries, contrary
%   to the standard practice in MATLAB of removing those entries from all
%   sparse matrices.  LD = sparse2 (LD) will remove any zero entries in LD.
%
%   See also LDLUPDATE, LDLSOLVE, LDLSPLIT, CHOL2, LCHOL, CHOL, SPARSE2

%   Copyright 2006-2007, Timothy A. Davis
%   http://www.cise.ufl.edu/research/sparse

error ('ldlchol mexFunction not found') ;


text:
function  numTSeries = countTSeries(view)
%numTSeries = countTSeries(view)
%
%Counts the tSeries
%
% 12.22.98 RFD: changed countFiles to work with new tSeries format ('.dat')
%				It should still work with old-style '.mat' tSeries files too.

if exist('view','var')
    subdir = view.subdir;
else
    subdir = 'Inplane';
end

[numScans,scanDirList] = countDirs('Scan*',tSeriesDir(view));
for scan=1:numScans
    dirPathStr = fullfile(tSeriesDir(view),scanDirList{scan});
    numTSeries(scan) = countFiles('tSeries*',dirPathStr);
end


text:
function [surf, labels] = imSurfaceAreaEstimate(img, varargin)
% Estimate surface area of a binary 3D structure.
%
%   Usage
%   S = imSurfaceAreaEstimate(IMG)
%   Estimate the surface area of the structure within the image, without
%   measuring surface area of borders.
%   The aim of this function is to be called by the "imSurfaceAreaDensity"
%   function, for providing an estimate of surface area density within a
%   representative volume of interest.
%
%
%   Example
%   imSurfaceAreaEstimate
%
%   See also
%   imSurfaceArea, imSurfaceAreaDensity
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2010-07-26,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRAE - Cepia Software Platform.

%% Process input arguments

% check image dimension
if ndims(img) ~= 3
    error('first argument should be a 3D image');
end

% in case of a label image, return a vector with a set of results
if ~islogical(img)
    labels = unique(img);
    labels(labels==0) = [];
    surf = zeros(length(labels), 1);
    for i = 1:length(labels)
        surf(i) = imSurfaceAreaEstimate(img==labels(i), varargin{:});
    end
    return;
end

%% Process input arguments

% in case of binary image, compute only one label...
labels = 1;

% default number of directions
nDirs = 13;

% default image resolution
delta = [1 1 1];

% Process user input arguments
while ~isempty(varargin)
    var = varargin{1};
    if ~isnumeric(var)
        error('option should be numeric');
    end

    % option is either connectivity or resolution
    if isscalar(var)
        nDirs = var;
    else
        delta = var;
    end
    varargin(1) = [];
end

%% Use the LUT to estimate surface area

lut = imSurfaceAreaLut(delta, nDirs);
bch = imBinaryConfigHisto(img);
surf = sum(bch .* lut);