text:
function [prior,binSize] = part_location_det(expidx)

p = rcnn_exp_params(expidx);
exp_dir = [p.expDir '/' p.shortName];

conf = rcnn_config('sub_dir', '/cachedir/test', 'exp_dir', exp_dir);
% saveTo = [conf.cache_dir '/gtHist/'];
saveTo = [conf.cache_dir '/detHist/'];
if (~exist(saveTo,'dir'))
    mkdir(saveTo);
end

% sc = p.refHeight/200;
sc = 2;
binSize = p.prior_bin_size*sc;
delta = 300*sc;
nBins = round(2*delta/binSize);

try
    assert(false);
%     load(fname, 'prior');
catch
%     load(p.testGT);
    load(p.trainGT);
%     if (exist('single_person_annolist','var'))
%         annolist = single_person_annolist;
%     end

    fnameDist = [fileparts(p.evalTest) '/distAll'];
    load(fnameDist, 'keypointsAll');
    assert(length(keypointsAll) == length(annolist));

    [~,parts] = util_get_parts24();
    nParts = length(p.pidxs);
    prior = zeros(2*delta,2*delta,nParts);
    histLoc = zeros(nBins,nBins,nParts);
    nEx = zeros(nParts,1);
    X0 = size(prior,1)/2;
    Y0 = size(prior,2)/2;

    for imgidx = 1:length(annolist)
        fprintf('.');
        rect = annolist(imgidx).annorect;
        for ridx = 1:length(rect)
            objpos = rect(ridx).objpos;
            points = rect(ridx).annopoints.point;
            for pidx = 1:length(parts)
                pidxLin = find(p.pidxs == parts(pidx).id);
                if (~isempty(pidxLin))
                    assert(parts(pidx).pos(1) == parts(pidx).pos(2));
%                     pp = util_get_annopoint_by_id(points,parts(pidx).pos(1));
                    pp = keypointsAll(imgidx).det(parts(pidx).pos(1)+1,1:2);
                    if (~isempty(pp))
                        iy = round(Y0 + sc*(pp(:,2) - objpos.y));
                        ix = round(X0 + sc*(pp(:,1) - objpos.x));
%                         iy = round(Y0 + sc*(pp.y - objpos.y));
%                         ix = round(X0 + sc*(pp.x - objpos.x));
                        iy = max(iy,1);
                        iy = min(iy,2*delta);
                        ix = max(ix,1);
                        ix = min(ix,2*delta);
                        iy_hist = ceil(iy/binSize);
                        ix_hist = ceil(ix/binSize);
                        histLoc(iy_hist,ix_hist,pidxLin) = histLoc(iy_hist,ix_hist,pidxLin) + 1;
                        iy_prior = (iy_hist-1)*binSize+1:iy_hist*binSize;
                        ix_prior = (ix_hist-1)*binSize+1:ix_hist*binSize;
                        prior(iy_prior,ix_prior,pidxLin) = prior(iy_prior,ix_prior,pidxLin) + ones(binSize,binSize);
                        nEx(pidxLin) = nEx(pidxLin) + 1;
                    end
                end
            end
        end
        if (~mod(imgidx, 100))
            fprintf(' %d/%d\n',imgidx,length(annolist));
        end
    end
    fprintf(' done\n');
    if (isfield(p,'sigma_gauss'))
        sigma_gauss = p.sigma_gauss;
    else
        sigma_gauss = binSize;
    end
    for pidx = 1:nParts
%         prior(:,:,pidx) = prior(:,:,pidx)/nEx(pidx);
        prior(:,:,pidx) = prior(:,:,pidx)/max(max(prior(:,:,pidx)));
%         pr = prior(:,:,pidx);
%         pr(pr > 0) = 1;
%         filterMask = gausswin(2*sigma_gauss+1)*gausswin(2*sigma_gauss+1)';
%         prior(:,:,pidx) = filter2(filterMask, pr);%/nEx(pidx);
%         prior(:,:,pidx) = prior(:,:,pidx)/max(max(prior(:,:,pidx)));
    end
%     save(fname, 'prior');
end

labels = {'rankle','rknee','rhip','lhip','lknee','lankle','rwrist','relbow','rshoulder','lshoulder','lelbow','lwrist','neck','tophead'};
tics = [1 size(prior,1)/6:size(prior,1)/6:size(prior,1)];
for pidx = 1:length(parts)
    figure(101);clf;%colormap gray;
    pidxLin = find(p.pidxs == parts(pidx).id);
    if (~isempty(pidxLin))
        pr = prior(:,:,pidxLin);
        imagesc(log(log(pr+1))); hold on; axis equal;
        plot([size(prior,2)/2; size(prior,2)/2], [1; size(prior,1)]', 'k-', 'lineWidth', 2);
        plot([1; size(prior,2)], [size(prior,1)/2; size(prior,1)/2]', 'k-', 'lineWidth', 2);
        title([labels{pidxLin}]);
        set(gca,'XTick',tics,'YTick',tics);
        fImgName = [saveTo '/partDetHist_' labels{pidxLin}];
        print(gcf, '-dpng', [fImgName '.png']);
    end
end
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [H] = lpfilter(type,M,N,D0,n)
%LPFILTER Computes freq. domain lowpass filters
%		THIS IS NOT A STANDARD MATLAB FUNCTION
%		H = lpfilter (type,M,N,D0,n) creates the
%		transfer function of a lowpass filter, H, of
%		the specified type and size MxN. Acceptable
%		values for type, D0, and n are:
%
%				'ideal'				Ideal lowpass filter with
%									cutoff frequency D0. If
%									supplied, n is ignored.
%				'btw'				Butterworth lowpass filter
%									of order n, and cutoff D0.
%				'gauss'			    Gaussian lowpass filter with
%									cutoff (standard deviation)D0.
%									If supplied, n is ignored.
%		M and N should be even numbers for DFT
%		filtering.
%
%		Class support: double, uint8, uint16
%		The output is of class double

%		Verify the number of inputs

error(nargchk(4,5,nargin));

%		Issue warning if M or N not even

L1 = M/2;
L2 = N/2;
if((L1-round(L1)) | (L2-round(L2)))
		disp('Warning: M & N must be even for DFT filtering')
end

%		Set up meshgrid to compute filter functions
%		running from 0 to M-1 and 0 to N-1. But,
%		recall that U and V run from 1 to M and 1 to N,
%		repectively.  See Section 2.10.4 for an
%		explanation of meshgrid. Recall also that
%		meshgrid reverses the order of rows and
%		columns, so the final result must be transposed
%		in order to preserve the original order of
%		of the coordinates.

[U,V] = meshgrid(0:M-1,0:N-1);

%   Compute the filter center

U0=L1;
V0=L2;

%		Compute distances D(U,V)

D=((U-U0).^2 + (V-V0).^2).^0.5;

%		Begin filter computations

switch type
	case 'ideal'
		H = D<=D0;
        H = H';
	case 'btw'
		H = 1./(1 + (D./D0).^(2*n));
        H = H';
	case 'gauss'
		H = exp(-(D.^2)./(2*(D0^2)));
        H = H';
	otherwise
		error('Unknown filter type')
end

%		End of function

repo: weiyanmin


text:
function [dist_matrix] = pos2dist(positions)
% pos2dist - Compute the relative distances matrix from the positions
%            matrix
%
% Inputs:
%   positions - matrix of size (3,nb_agents)
%
% Ouputs:
%   dist_matrix - symmetric matrix with pairwise distances
%
%
    distances = pdist(positions');
    dist_matrix = squareform(distances);

end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
%compute velang_head
function [data,units]=compute_velanghead(trx,n)
larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
velanghead=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    velanghead{1,i}=bsxfun(@atan2,trx(larva).dyhead_mm,trx(larva).dxhead_mm);
end

units=parseunits('rad');
data=velanghead;

author: kristinbranson

repo: JAABA


text:
function hermite_rule ( n, a, b, scale, filename )

%*****************************************************************************80
%
%% HERMITE_RULE generates a Gauss-Hermite rule.
%
%  Discussion:
%
%    This program computes a Gauss-Hermite quadrature rule
%    and writes it to a file.
%
%    The integral to be approximated has the form
%
%      C * integral ( -oo < x < +oo ) f(x) rho(x) dx
%
%    where the weight rho(x) is:
%
%      rho(x) = exp ( - b * ( x - a )^2 ) * sqrt ( b / pi ) dx
%
%    and A and B are parameters.
%
%    The user specifies:
%    * N, the number of points in the rule;
%    * A, the center point;
%    * B, the scale factor in the exponential;
%    * SCALE, is 1 if the weights are to be normalized;
%    * FILENAME, the root name of the output files.
%
%    If SCALE = 0, then the factor C in front of the integrand is 1.
%    If SCALE is nonzero, then the factor C is sqrt ( B ) / sqrt ( PI ).
%    which means that the function f(x)=1 will integrate to 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_RULE\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compute a Gauss-Hermite rule for approximating\n' );
  fprintf ( 1, '    integral ( -oo < x < +oo ) f(x) rho(x) dx\n' );
  fprintf ( 1, '  where the weight rho(x) is:\n' );
  fprintf ( 1, '    exp ( - b * ( x - a )^2 ) * sqrt ( b / pi ) dx\n' );
  fprintf ( 1, '  using N points.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The user specifies N, A, B, SCALE, FILENAME.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  N is the number of points;\n' );
  fprintf ( 1, '  A is the center point (typically 0).\n' );
  fprintf ( 1, '  B is the exponential scale factor (typically 1).\n' );
  fprintf ( 1, '  SCALE is 1 if the weights are to be normalized.\n');
  fprintf ( 1, '  FILENAME is used to generate 3 files:\n' );
  fprintf ( 1, '    filename_w.txt - the weight file\n' );
  fprintf ( 1, '    filename_x.txt - the abscissa file.\n' );
  fprintf ( 1, '    filename_r.txt - the region file.\n' );
%
%  Initialize the parameters.
%
  alpha = 0.0;
  beta = 0.0;
%
%  Get N.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  N is the number of points in the rule, such as 5.\n' );
    n = input ( '  Enter N: ' );
  elseif ( ischar ( n ) )
    n = str2num ( n );
  end
%
%  Get A.
%
  if ( nargin < 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  A is the center point, typically 0.\n' );
    fprintf ( 1, '\n' );
    a = input ( '  Enter A: ' );
  elseif ( ischar ( a ) )
    a = str2num ( a );
  end
%
%  Get B.
%
  if ( nargin < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  B is the scale factor, typically 0.5 or 1.0.\n' );
    fprintf ( 1, '\n' );
    b = input ( '  Enter B: ' );
  elseif ( ischar ( b ) )
    b = str2num ( b );
  end
%
%  Get SCALE.
%
  if ( nargin < 4 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  SCALE is 1 to normalize the weights, 0 otherwise.\n' );
    fprintf ( 1, '\n' );
    b = input ( '  Enter SCALE (0/1): ' );
  elseif ( ischar ( scale ) )
    scale = str2num ( scale );
  end
%
%  Get FILENAME.
%
  if ( nargin < 5 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  FILENAME is the ''root name'' of the quadrature files).\n' );
    filename = input ( '  Enter FILENAME as a quoted string:' );
  end
%
%  Input summary.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Input summary:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  N = %d\n', n );
  fprintf ( 1, '  A = %f\n', a );
  fprintf ( 1, '  B = %f\n', b );
  fprintf ( 1, '  SCALE = %d\n', scale );
  fprintf ( 1, '  FILENAME = "%s".\n', filename );
%
%  Construct the rule.
%
  kind = 6;
  [ x, w ] = cgqf ( n, kind, alpha, beta, a, b );
%
%  Normalize the rule.
%  This way, the weights add up to 1.
%
  if ( scale == 1.0 )
    w(1:n) = w(1:n) * sqrt ( b ) / sqrt ( pi );
  end
%
%  Set the R values to suggest an infinite interval.
%
  r = zeros ( 2, 1 );
  r(1) = - r8_huge ( );
  r(2) =   r8_huge ( );
%
%  Write the rule.
%
  rule_write ( n, filename, x, w, r );
%
%  Terminate.
%
  fprintf ( 1,  '\n' );
  fprintf ( 1,  'HERMITE_RULE:\n' );
  fprintf ( 1,  '  Normal end of execution.\n' );
  fprintf ( 1,  '\n' );
  timestamp ( );

  return
end
function [ t, wts ] = cdgqf ( nt, kind, alpha, beta )

%*****************************************************************************80
%
%% CDGQF computes a Gauss quadrature formula with default A, B and simple knots.
%
%  Discussion:
%
%    This routine computes all the knots and weights of a Gauss quadrature
%    formula with a classical weight function with default values for A and B,
%    and only simple knots.
%
%    There are no moments checks and no printing is done.
%
%    Use routine EIQFS to evaluate a quadrature computed by CGQFS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev,            (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Output, real T(NT), the knots.
%
%    Output, real WTS(NT), the weights.
%
  parchk ( kind, 2 * nt, alpha, beta );
%
%  Get the Jacobi matrix and zero-th moment.
%
  [ aj, bj, zemu ] = class_matrix ( kind, nt, alpha, beta );
%
%  Compute the knots and weights.
%
  [ t, wts ] = sgqf ( nt, aj, bj, zemu );

  return
end
function [ t, wts ] = cgqf ( nt, kind, alpha, beta, a, b )

%*****************************************************************************80
%
%% CGQF computes knots and weights of a Gauss quadrature formula.
%
%  Discussion:
%
%    The user may specify the interval (A,B).
%
%    Only simple knots are produced.
%
%    The user may request that the routine print the knots and weights,
%    and perform a moment check.
%
%    Use routine EIQFS to evaluate this quadrature formula.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,+oo)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-oo,+oo)   |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,+oo)     (x-a)^alpha*(x+b)^beta
%    9, Chebyshev Type 2,     (a,b)       ((b-x)*(x-a))^(+0.5)
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Input, real A, B, the interval endpoints.
%
%    Output, real T(NT), the knots.
%
%    Output, real WTS(NT), the weights.
%

%
%  Compute the Gauss quadrature formula for default values of A and B.
%
  [ t, wts ] = cdgqf ( nt, kind, alpha, beta );
%
%  All knots have multiplicity = 1.
%
  mlt = zeros(nt,1);
  mlt(1:nt) = 1;
%
%  NDX(I) = I.
%
  ndx = ( 1 : nt );
%
%  Scale the quadrature rule.
%
  [ t, wts ] = scqf ( nt, t, mlt, wts, nt, ndx, kind, alpha, beta, a, b );

  return
end
function [ aj, bj, zemu ] = class_matrix ( kind, m, alpha, beta )

%*****************************************************************************80
%
%% CLASS_MATRIX computes the Jacobi matrix for a quadrature rule.
%
%  Discussion:
%
%    This routine computes the diagonal AJ and subdiagonal BJ
%    elements of the order M tridiagonal symmetric Jacobi matrix
%    associated with the polynomials orthogonal with respect to
%    the weight function specified by KIND.
%
%    For weight functions 1-7, M elements are defined in BJ even
%    though only M-1 are needed.  For weight function 8, BJ(M) is
%    set to zero.
%
%    The zero-th moment of the weight function is returned in ZEMU.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev,            (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
%
%    Input, integer M, the order of the Jacobi matrix.
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Output, real AJ(M), BJ(M), the diagonal and subdiagonal
%    of the Jacobi matrix.
%
%    Output, real ZEMU, the zero-th moment.
%
  temp = eps;

  parchk ( kind, 2 * m - 1, alpha, beta );

  temp2 = 0.5;

  if ( 500.0 * temp < abs ( ( gamma ( temp2 ) )^2 - pi ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CLASS - Fatal error!\n' );
    fprintf ( 1, '  Gamma function does not match machine parameters.\n' );
    error ( 'CLASS - Fatal error!' );
  end

  bj = zeros(m,1);
  aj = zeros(m,1);

  if ( kind == 1 )

    ab = 0.0;

    zemu = 2.0 / ( ab + 1.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      abi = i + ab * mod ( i, 2 );
      abj = 2 * i + ab;
      bj(i) = abi * abi / ( abj * abj - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 2 )

    zemu = pi;

    aj(1:m) = 0.0;

    bj(1) =  sqrt ( 0.5 );
    bj(2:m) = 0.5;

  elseif ( kind == 3 )

    ab = alpha * 2.0;
    zemu = 2.0^( ab + 1.0 ) * gamma ( alpha + 1.0 )^2 ...
      / gamma ( ab + 2.0 );

    aj(1:m) = 0.0;
    bj(1) = 1.0 / ( 2.0 * alpha + 3.0 );
    for i = 2 : m
      bj(i) = i * ( i + ab ) / ( 4.0 * ( i + alpha )^2 - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 4 )

    ab = alpha + beta;
    abi = 2.0 + ab;
    zemu = 2.0^( ab + 1.0 ) * gamma ( alpha + 1.0 ) ...
      * gamma ( beta + 1.0 ) / gamma ( abi );
    aj(1) = ( beta - alpha ) / abi;
    bj(1) = 4.0 * ( 1.0 + alpha ) * ( 1.0 + beta ) ...
      / ( ( abi + 1.0 ) * abi * abi );
    a2b2 = beta * beta - alpha * alpha;

    for i = 2 : m
      abi = 2.0 * i + ab;
      aj(i) = a2b2 / ( ( abi - 2.0 ) * abi );
      abi = abi^2;
      bj(i) = 4.0 * i * ( i + alpha ) * ( i + beta ) * ( i + ab ) ...
        / ( ( abi - 1.0 ) * abi );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 5 )

    zemu = gamma ( alpha + 1.0 );

    for i = 1 : m
      aj(i) = 2.0 * i - 1.0 + alpha;
      bj(i) = i * ( i + alpha );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 6 )

    zemu = gamma ( ( alpha + 1.0 ) / 2.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      bj(i) = ( i + alpha * mod ( i, 2 ) ) / 2.0;
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 7 )

    ab = alpha;
    zemu = 2.0 / ( ab + 1.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      abi = i + ab * mod(i,2);
      abj = 2 * i + ab;
      bj(i) = abi * abi / ( abj * abj - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 8 )

    ab = alpha + beta;
    zemu = gamma ( alpha + 1.0 ) * gamma ( - ( ab + 1.0 ) ) ...
      / gamma ( - beta );
    apone = alpha + 1.0;
    aba = ab * apone;
    aj(1) = - apone / ( ab + 2.0 );
    bj(1) = - aj(1) * ( beta + 1.0 ) / ( ab + 2.0 ) / ( ab + 3.0 );
    for i = 2 : m
      abti = ab + 2.0 * i;
      aj(i) = aba + 2.0 * ( ab + i ) * ( i - 1 );
      aj(i) = - aj(i) / abti / ( abti - 2.0 );
    end

    for i = 2 : m - 1
      abti = ab + 2.0 * i;
      bj(i) = i * ( alpha + i ) / ( abti - 1.0 ) * ( beta + i ) ...
        / ( abti^2 ) * ( ab + i ) / ( abti + 1.0 );
    end

    bj(m) = 0.0;
    bj(1:m) =  sqrt ( bj(1:m) );

  end

  return
end
function [ d, z ] = imtqlx ( n, d, e, z )

%*****************************************************************************80
%
%% IMTQLX diagonalizes a symmetric tridiagonal matrix.
%
%  Discussion:
%
%    This routine is a slightly modified version of the EISPACK routine to
%    perform the implicit QL algorithm on a symmetric tridiagonal matrix.
%
%    The authors thank the authors of EISPACK for permission to use this
%    routine.
%
%    It has been modified to produce the product Q' * Z, where Z is an input
%    vector and Q is the orthogonal matrix diagonalizing the input matrix.
%    The changes consist (essentialy) of applying the orthogonal transformations
%    directly to Z as they are generated.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%    Roger Martin, James Wilkinson,
%    The Implicit QL Algorithm,
%    Numerische Mathematik,
%    Volume 12, Number 5, December 1968, pages 377-383.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real D(N), the diagonal entries of the matrix.
%
%    Input, real E(N), the subdiagonal entries of the
%    matrix, in entries E(1) through E(N-1).
%
%    Input, real Z(N), a vector to be operated on.
%
%    Output, real D(N), the diagonal entries of the diagonalized matrix.
%
%    Output, real Z(N), the value of Q' * Z, where Q is the matrix that
%    diagonalizes the input symmetric tridiagonal matrix.
%
  itn = 30;

  prec = eps;

  if ( n == 1 )
    return
  end

  e(n) = 0.0;

  for l = 1 : n

    j = 0;

    while ( 1 )

      for m = l : n

        if ( m == n )
          break
        end

        if ( abs ( e(m) ) <= prec * ( abs ( d(m) ) + abs ( d(m+1) ) ) )
          break
        end

      end

      p = d(l);

      if ( m == l )
        break
      end

      if ( j == itn )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'IMTQLX - Fatal error!\n' );
        fprintf ( 1, '  Iteration limit exceeded.\n' );
        error ( 'IMTQLX - Fatal error!' );
      end

      j = j + 1;
      g = ( d(l+1) - p ) / ( 2.0 * e(l) );
      r =  sqrt ( g * g + 1.0 );
      g = d(m) - p + e(l) / ( g + r8_sign ( g ) * abs ( r ) );
      s = 1.0;
      c = 1.0;
      p = 0.0;
      mml = m - l;

      for ii = 1 : mml

        i = m - ii;
        f = s * e(i);
        b = c * e(i);

        if ( abs ( f ) >= abs ( g ) )
          c = g / f;
          r =  sqrt ( c * c + 1.0 );
          e(i+1) = f * r;
          s = 1.0 / r;
          c = c * s;
        else
          s = f / g;
          r =  sqrt ( s * s + 1.0 );
          e(i+1) = g * r;
          c = 1.0 / r;
          s = s * c;
        end

        g = d(i+1) - p;
        r = ( d(i) - g ) * s + 2.0 * c * b;
        p = s * r;
        d(i+1) = g + p;
        g = c * r - b;
        f = z(i+1);
        z(i+1) = s * z(i) + c * f;
        z(i) = c * z(i) - s * f;

      end

      d(l) = d(l) - p;
      e(l) = g;
      e(m) = 0.0;

    end

  end

  for ii = 2 : n

     i = ii - 1;
     k = i;
     p = d(i);

     for j = ii : n
       if ( d(j) < p )
         k = j;
         p = d(j);
       end
     end

     if ( k ~= i )
       d(k) = d(i);
       d(i) = p;
       p = z(i);
       z(i) = z(k);
       z(k) = p;
     end

  end

  return
end
function parchk ( kind, m, alpha, beta )

%*****************************************************************************80
%
%% PARCHK checks parameters ALPHA and BETA for classical weight functions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev,            (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
%
%    Input, integer M, the order of the highest moment to
%    be calculated.  This value is only needed when KIND = 8.
%
%    Input, real ALPHA, BETA, the parameters, if required
%    by the value of KIND.
%
  if ( kind <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARCHK - Fatal error!\n' );
    fprintf ( 1, '  KIND <= 0.\n' );
    error ( 'PARCHK - Fatal error!' );
  end
%
%  Check ALPHA for Gegenbauer, Jacobi, Laguerre, Hermite, Exponential.
%
  if ( 3 <= kind && alpha <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARCHK - Fatal error!\n' );
    fprintf ( 1, '  3 <= KIND and ALPHA <= -1.\n' );
    error ( 'PARCHK - Fatal error!' );
  end
%
%  Check BETA for Jacobi.
%
  if ( kind == 4 && beta <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARCHK - Fatal error!\n' );
    fprintf ( 1, '  KIND == 4 and BETA <= -1.0.\n' );
    error ( 'PARCHK - Fatal error!' );
  end
%
%  Check ALPHA and BETA for rational.
%
  if ( kind == 8 )
    tmp = alpha + beta + m + 1.0;
    if ( 0.0 <= tmp || tmp <= beta )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PARCHK - Fatal error!\n' );
      fprintf ( 1, '  KIND == 8 but condition on ALPHA and BETA fails.\n' );
      error ( 'PARCHK - Fatal error!' );
    end
  end

  return
end
function value = r8_huge ( )

%*****************************************************************************80
%
%% R8_HUGE returns a "huge" real number.
%
%  Discussion:
%
%    The value returned by this function is NOT required to be the
%    maximum representable R8.  This value varies from machine to machine,
%    from compiler to compiler, and may cause problems when being printed.
%    We simply want a "very large" but non-infinite number.
%
%    MATLAB provides a built-in symbolic constant "inf" that can be used
%    if a huge number is really what you want!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VALUE, a huge number.
%
  value = 1.0E+30;

  return
end
function value = r8_sign ( x )

%*****************************************************************************80
%
%% R8_SIGN returns the sign of an R8.
%
%  Discussion:
%
%    The value is +1 if the number is positive or zero, and it is -1 otherwise.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the number whose sign is desired.
%
%    Output, real VALUE, the sign of X.
%
  if ( 0 <= x )
    value = +1.0;
  else
    value = -1.0;
  end

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
%  For smaller data files, and less precision, try:
%
%     fprintf ( output_unit, '  %14.6f', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %24.16f', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function rule_write ( order, filename, x, w, r )

%*****************************************************************************80
%
%% RULE_WRITE writes a quadrature rule to a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ORDER, the order of the rule.
%
%    Input, string FILENAME, specifies the output files.
%    write files 'filename_w.txt', 'filename_x.txt', 'filename_r.txt' defining
%    weights, abscissas, and region.
%
%    Input, real X(ORDER), the abscissas.
%
%    Input, real W(ORDER), the weights.
%
%    Input, real R(2), the region.
%
  filename_x = strcat ( filename, '_x.txt' );
  filename_w = strcat ( filename, '_w.txt' );
  filename_r = strcat ( filename, '_r.txt' );

  fprintf ( 1, '\n' );
  fprintf ( 1,'  Creating quadrature files.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  "Root" file name is   "%s".\n', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Weight file will be   "%s".\n', filename_w );
  fprintf ( 1, '  Abscissa file will be "%s".\n', filename_x );
  fprintf ( 1, '  Region file will be   "%s".\n', filename_r );

  r8mat_write ( filename_w, 1, order, w' );
  r8mat_write ( filename_x, 1, order, x' );
  r8mat_write ( filename_r, 1, 2,     r' );

  return
end
function [ t, wts ] = scqf ( nt, t, mlt, wts, nwts, ndx, kind, alpha, ...
  beta, a, b )

%*****************************************************************************80
%
%% SCQF scales a quadrature formula to a nonstandard interval.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, real T(NT), the original knots.
%
%    Input, integer MLT(NT), the multiplicity of the knots.
%
%    Input, real WTS(NWTS), the weights.
%
%    Input, integer NWTS, the number of weights.
%
%    Input, integer NDX(NT), used to index the array WTS.
%    For more details see the comments in CAWIQ.
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,+oo)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-oo,+oo)   |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,+oo)     (x-a)^alpha*(x+b)^beta
%    9, Chebyshev Type 2,     (a,b)       ((b-x)*(x-a))^(+0.5)
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Input, real A, B, the interval endpoints.
%
%    Output, real T(NT), the scaled knots.
%
%    Output, real WTS(NWTS), the scaled weights.
%
  temp = eps;

  parchk ( kind, 1, alpha, beta )

  if ( kind == 1 )

    al = 0.0;
    be = 0.0;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 2 )

    al = -0.5;
    be = -0.5;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 3 )

    al = alpha;
    be = alpha;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 4 )

    al = alpha;
    be = beta;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 5 )

    if ( b <= 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  B <= 0.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = a;
    slp = 1.0 / b;
    al = alpha;
    be = 0.0;

  elseif ( kind == 6 )

    if ( b <= 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  B <= 0.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = a;
    slp = 1.0 / sqrt ( b );
    al = alpha;
    be = 0.0;

  elseif ( kind == 7 )

    al = alpha;
    be = 0.0;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 8 )

    if ( a + b <= 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  A + B <= 0.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = a;
    slp = a + b;
    al = alpha;
    be = beta;

  elseif ( kind == 9 )

    al = 0.5;
    be = 0.5;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  end

  p = slp^( al + be + 1.0 );

  for k = 1 : nt

    t(k) = shft + slp * t(k);
    l = abs ( ndx(k) );

    if ( l ~= 0 )
      tmp = p;
      for i = l : l + mlt(k) - 1
        wts(i) = wts(i) * tmp;
        tmp = tmp * slp;
      end
    end

  end

  return
end
function [ t, wts ] = sgqf ( nt, aj, bj, zemu )

%*****************************************************************************80
%
%% SGQF computes knots and weights of a Gauss Quadrature formula.
%
%  Discussion:
%
%    This routine computes all the knots and weights of a Gauss quadrature
%    formula with simple knots from the Jacobi matrix and the zero-th
%    moment of the weight function, using the Golub-Welsch technique.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, real AJ(NT), the diagonal of the Jacobi matrix.
%
%    Input, real BJ(NT), the subdiagonal of the Jacobi
%    matrix, in entries 1 through NT-1.  On output, BJ has been overwritten.
%
%    Input, real ZEMU, the zero-th moment of the weight function.
%
%    Output, real T(NT), the knots.
%
%    Output, real WTS(NT), the weights.
%

%
%  Exit if the zero-th moment is not positive.
%
  if ( zemu <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SGQF - Fatal error!\n' );
    fprintf ( 1, '  ZEMU <= 0.\n' );
    error ( 'SGQF - Fatal error!' );
  end
%
%  Set up vectors for IMTQLX.
%
  wts = zeros ( nt, 1 );

  wts(1) = sqrt ( zemu );
  wts(2:nt) = 0.0;
%
%  Diagonalize the Jacobi matrix.
%
  [ t, wts ] = imtqlx ( nt, aj, bj, wts );

  wts(1:nt) = wts(1:nt).^2;

  return
end
function timestamp ( )

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

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function poly = selectPolygon(varargin)
% select a polygon by mouse
%
% Syntax
%   poly = selectPolygon
%   ebsd = ebsd(inpolygon(ebsd,poly))

disp('Use the mouse to select the vertices of a polygon.');
disp('Click right to finish your selection.')

poly = [];
button = 1;
hold on

h = line('xdata',[],'ydata',[],'Color','k','MarkerFaceColor','k','MarkerEdgeColor',...
  'w','MarkerSize',10,'LineWidth',2,'Marker','s');

while button == 1

    [xi,yi,button] = ginput(1);
    poly(end+1,:) = [xi,yi]; %#ok<AGROW>

    set(h,'XData',poly(:,1),'YData',poly(:,2));

end

% ensure the polygon is always closed
poly(end+1,:) = poly(1,:);
set(h,'XData',poly(:,1),'YData',poly(:,2));

hold off

author: mtex-toolbox

repo: mtex

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


text:
function [HDR] = leadidcodexyz(arg1)
% LeadIdCodeXYZ uses the Label information for computing the
%       LeadIdCode and the XYZ position of the EEG Electrodes
%       according to Annex A of FEF Vital Signs Format [1]
%
%   HDR = leadidcodexyz(HDR);
%	adds HDR.LeadIdCode and HDR.ELEC.XYZ, if needed.
%
% see also: SLOAD, SOPEN, PHYSICALUNITS, doc/leadidtable_scpecg.txt, doc/elecpos.txt
%
% Reference(s):
% [1] CEN/TC251/PT40 (2001)
% 	File Exchange Format for Vital Signs - Annex A
%
% Birbaumer, N. (2006). Brain-computer-interface research: Coming of age. Clinical Neurophysiology, 117:479??3.
%  http://www.acns.org/pdfs/ACFDD46.pdf.
% ACNS (2006). Guidelines for standard electrode position nomenclature. American Clinical
% Neurophysiology Society. http://www.acns.org/pdfs/ACFDD46.pdf.

%
% 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.

%	$Id: leadidcodexyz.m 2602 2011-02-03 13:05:03Z schloegl $
%	Copyright (C) 2006,2007,2008,2009 by Alois Schloegl <a.schloegl@ieee.org>
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/

global BIOSIG_GLOBAL;
% BIOSIG_GLOBAL=[]; %%% used for debugging, only.

if ~isfield(BIOSIG_GLOBAL,'Phi');
	BIOSIG_GLOBAL.ISLOADED_XYZ = 0 ;
end;
if ~isfield(BIOSIG_GLOBAL,'ISLOADED_XYZ')
	BIOSIG_GLOBAL.ISLOADED_XYZ = 0 ;
end;
if ~BIOSIG_GLOBAL.ISLOADED_XYZ;
    f = which('getfiletype.m'); 	% identify path to biosig
    [p,f,e] = fileparts(f);
    [p,f,e] = fileparts(p);

    BIOSIG_GLOBAL.ISLOADED_XYZ = 0 ;

    N = 0;
    fid = fopen(fullfile(p,f,'leadidtable_scpecg.txt'),'r');
    s = char(fread(fid,[1,inf],'uint8'));
    fclose(fid);

    Code = repmat(NaN, 200, 1); Phi = Code; Theta = Code;
    while ~isempty(s),
        [t,s] = strtok(s,[10,13]);
        if ~length(t)
        elseif ~strncmp(t,'#',1)
            ix3 = strfind(t,'MDC_ECG_LEAD_');
            if isempty(ix3)
                ix3 = length(t)+1;
            end
            [t1,t2] = strtok(t(1:ix3-1),[9,32]);
            [t2,t3] = strtok(t2,[9,32]);
            id = str2double(t2);
            N  = N + 1;
            Labels{N,1}	  = t1;
            Code(N,1)         = id;
            Description{N,1}  = deblank(t3);
            %MDC_ECG_LEAD{N,1} = t(ix3+13:end)
            MDC_ECG_LEAD{N,1} = t(ix3:end);
        end;
    end;
    N1 = N;

        % load table
        fid = fopen(fullfile(p,f,'elecpos.txt'),'r');
        t = char(fread(fid,[1,inf],'uint8'));
        fclose(fid);

        % extract table information
        while ~isempty(t)
                [x,t] = strtok(t,[10,13]);
                if isempty(x)
                elseif strncmp(x,'#',1)
                else
                        N = N + 1;
                        [num,status,strarray] = str2double(x);
                        Code(N,1)   = num(1);
                        Labels{N,1} = upper(strarray{2});
                        Phi(N,1)    = num(3);
                        Theta(N,1)  = num(4);
                end;
        end;
        Phi   = Phi(:)  *pi/180;
        Theta = Theta(:)*pi/180;

        % loading is done only once.
        BIOSIG_GLOBAL.XYZ = [sin(Theta).*cos(Phi), sin(Theta).*sin(Phi), cos(Theta)];
        BIOSIG_GLOBAL.Phi          = Phi*180/pi;
        BIOSIG_GLOBAL.Theta        = Theta*180/pi;
        BIOSIG_GLOBAL.LeadIdCode   = Code;
        BIOSIG_GLOBAL.Label        = Labels;
        BIOSIG_GLOBAL.Description  = Description;
        BIOSIG_GLOBAL.MDC_ECG_LEAD = MDC_ECG_LEAD;

        BIOSIG_GLOBAL.ISLOADED_XYZ = 1;
end;

if nargin<1,
        HDR.LeadIdCode = BIOSIG_GLOBAL.LeadIdCode;
        HDR.Label      = BIOSIG_GLOBAL.Label;
        HDR.ELEC.XYZ   = BIOSIG_GLOBAL.XYZ;
        HDR.TYPE       = 'ELPOS';

else    % electrode code and position

        if isstruct(arg1)
                HDR = arg1;
        elseif isnumeric(arg1),
                HDR.LeadIdCode = arg1;
        else
                HDR.Label = arg1;
        end;

        tmp.flag1 = isfield(HDR,'ELEC');
        if tmp.flag1,
                tmp.flag1 = isfield(HDR.ELEC,'XYZ');
        end;
        if tmp.flag1,
                tmp.flag1 = any(HDR.ELEC.XYZ(:));
        end;
        tmp.flag2 = isfield(HDR,'LeadIdCode');
        tmp.flag3 = isfield(HDR,'Label');

        if (~tmp.flag1 || ~tmp.flag2 || ~tmp.flag3),
        	if 0,
        	elseif tmp.flag3,
                        if ischar(HDR.Label)
                                HDR.Label = cellstr(HDR.Label);
                        end;
	                NS = length(HDR.Label);
	        elseif tmp.flag2,
	        	NS = length(HDR.LeadIdCode);
        	elseif isfield(HDR,'NS')
        		NS = HDR.NS;
        		HDR.LeadIdCode = zeros(1,HDR.NS);
	        end;

                if tmp.flag3,
	                if ~tmp.flag1,
				HDR.ELEC.XYZ   = repmat(NaN,NS,3);
				HDR.ELEC.Phi   = repmat(NaN,NS,1);
				HDR.ELEC.Theta = repmat(NaN,NS,1);
	        	end;
                	if ~tmp.flag2,
                       		HDR.LeadIdCode = repmat(NaN,NS,1);
	        	end;
                	for k = 1:NS;
				Label = upper(deblank(HDR.Label{k}));
				pos = find(Label==':');
				if ~isempty(pos)
					Label = Label(pos+1:end);
                            	end;
	                        ix = strmatch(Label,BIOSIG_GLOBAL.Label,'exact');

	                        if length(ix)==2,
	                        	%%%%% THIS IS A HACK %%%%%
	                        	%% solve ambiguity for 'A1','A2'; could be EEG or ECG
	                        	if sum(HDR.LeadIdCode(1:k)>=996)>sum(HDR.LeadIdCode(1:k)<996)
	                        		%% majority are EEG electrodes,
	                        		ix = ix(find(BIOSIG_GLOBAL.LeadIdCode(ix)>996));
	                        	else
	                        		%% majority are ECG electrodes,
	                        		ix = ix(find(BIOSIG_GLOBAL.LeadIdCode(ix)<996));
	                        	end;
	                        elseif isempty(ix)
		                        ix = strmatch(deblank(HDR.Label{k}),BIOSIG_GLOBAL.MDC_ECG_LEAD,'exact');
	                        end;

        	                if (length(ix)==1),
	                	        if ~tmp.flag1,
        	                	        HDR.ELEC.XYZ(k,1:3) = BIOSIG_GLOBAL.XYZ(ix,:);
                	        	        HDR.ELEC.Phi(k)   = BIOSIG_GLOBAL.Phi(ix);
                        		        HDR.ELEC.Theta(k) = BIOSIG_GLOBAL.Theta(ix);
                        		end;
                	        	if ~tmp.flag2,
                                		HDR.LeadIdCode(k,1) = BIOSIG_GLOBAL.LeadIdCode(ix);
	                        	end;
        	                end;
                        end;
                else
			HDR.Label = cell(NS,1);
			for k = 1:NS;
				ix = find(BIOSIG_GLOBAL.LeadIdCode==HDR.LeadIdCode(k));
				if (length(ix)>=1),
					ix=ix(1);	%
	                                HDR.Label{k} = BIOSIG_GLOBAL.Label{ix};
					if ~tmp.flag1,
        		                        HDR.ELEC.XYZ(k,1:3) = BIOSIG_GLOBAL.XYZ(ix,1:3);
        		                        HDR.ELEC.Phi(k,1)   = BIOSIG_GLOBAL.Phi(ix);
        		                        HDR.ELEC.Theta(k,1) = BIOSIG_GLOBAL.Theta(ix);
					end;
				else
					HDR.Label{k} = ['#',int2str(k)];
				end;
			end;
                end;
	end;
	if tmp.flag3 && ~any(HDR.LeadIdCode),
               	for k = 1:HDR.NS;
			ix = strmatch(upper(HDR.Label{k}),BIOSIG_GLOBAL.Label,'exact');
			if length(ix)==1,
				HDR.LeadIdCode(k) = BIOSIG_GLOBAL.LeadIdCode(ix);
			end;
		end;
        end;

end;

author: PatternRecognition

repo: OpenBMI


text:
function [acre] = m22acre(m2)
% Convert area from square meters to acres.
% Chad A. Greene 2012
acre = m2*0.0002471053814672;

author: Sable

repo: mcbench-benchmarks


text:
% Information-Flow Matte Refinement
% This function implements the matte refinement approach described in
% Yagiz Aksoy, Tunc Ozan Aydin, Marc Pollefeys, "Designing Effective
% Inter-Pixel Information Flow for Natural Image Matting", CVPR, 2017.
% 'alphaHat' and 'confidences' parameters are typically obtained by a
% sampling-based natural matting algorithm. 'confidences' is filled by
% ones if not provided. Optional input parameter 'params' can be
% customized by editing the default values in the struct returned
% by 'getMattingParams('IFM').
% - **_K parameters represent the number of nonlocal neighbors found
%   for color mixture,  and intra-U flows, while **_xyw define the effect of
%   spatial proximity.
% - **_mult define the weight of each information flow.
% - loc_*** define the parameters for the matting Laplacian.
% - refinement_mult determines how much trust is given to the initial
%   alpha estimation

function alpha = informationFlowMatteRefinement(image, trimap, alphaHat, confidences, params, suppressMessages)
    abmtSetup
    tic;
    if ~exist('confidences', 'var') || isempty(confidences)
        confidences = ones(size(alphaHat(:,:,1)));
    end
    if ~exist('params', 'var') || isempty(params)
        params = getMattingParams('IFM');
    end
    if ~exist('suppressMessages', 'var') || isempty(suppressMessages)
        suppressMessages = false;
    end
    if(~suppressMessages) display('Matte refinement via Information-Flow Matting...'); end

    image = im2double(image);
    trimap = im2double(trimap(:,:,1));
    alphaHat = im2double(alphaHat(:,:,1));

    % Compute L_IFM
    unk = trimap < 0.8 & trimap > 0.2;
    dilUnk = imdilate(unk, ones(3, 3));
    if(~suppressMessages) display('     Computing color mixture flow...'); end
    Lap = affinityMatrixToLaplacian(colorMixtureAffinities(image, params.cm_K, dilUnk, [], params.cm_xyw));
    Lap = params.cm_mult * (Lap' * Lap);
    if(~suppressMessages) display('     Computing matting Laplacian...'); end
    Lap = Lap + params.loc_mult * affinityMatrixToLaplacian(mattingAffinity(image, dilUnk, params.loc_win, params.loc_eps));
    if(~suppressMessages) display('     Computing intra-U flow...'); end
    Lap = Lap + params.iu_mult * affinityMatrixToLaplacian(colorSimilarityAffinities(image, params.iu_K, unk, unk, params.iu_xyw));

    if(~suppressMessages) display('     Solving for alphas...'); end
    alpha = solveForAlphas(Lap, trimap, params.lambda, params.usePCGtoSolve, alphaHat, confidences, params.refinement_mult);

    alpha = reshape(alpha, [size(image, 1), size(image, 2)]);

    dur = toc;
    if(~suppressMessages) display(['Done. It took ' num2str(dur) ' seconds.']); end
end

author: yaksoy

repo: AffinityBasedMattingToolbox


text:
classdef SetPolyThrottleModelActionOptimVar < AbstractOptimizationVariable
    %SetRPYSteeringModelActionOptimVar Summary of this class goes here
    %   Detailed explanation goes here

    properties
        varObj = ThrottlePolyModel.getDefaultThrottleModel()

        lb(1,4) double
        ub(1,4) double

        varConst(1,1) logical = false;
        varLin(1,1) logical   = false;
        varAccel(1,1) logical = false;

        varTimeOffset(1,1) logical = false;
    end

    methods
        function obj = SetPolyThrottleModelActionOptimVar(varObj)
            obj.varObj = varObj;
            obj.varObj.optVar = obj;

            obj.id = rand();
        end

        function x = getXsForVariable(obj)
            x = NaN(1,4);

            if(obj.varConst)
                x(1) = obj.varObj.throttleModel.constTerm;
            end
            if(obj.varLin)
                x(2) = obj.varObj.throttleModel.linearTerm;
            end
            if(obj.varAccel)
                x(3) = obj.varObj.throttleModel.accelTerm;
            end

            if(obj.varTimeOffset)
                x(4) = obj.varObj.throttleModel.tOffset;
            end

            x(isnan(x)) = [];
        end

        function [lb, ub] = getBndsForVariable(obj)
            useTf = obj.getUseTfForVariable();

            lb = obj.lb(useTf);
            ub = obj.ub(useTf);
        end

        function [lb, ub] = getAllBndsForVariable(obj)
            lb = obj.lb;
            ub = obj.lb;
        end

        function setBndsForVariable(obj, lb, ub)
            if(length(lb) == 4 && length(ub) == 4)
                obj.lb = lb;
                obj.ub = ub;
            else
                useTf = obj.getUseTfForVariable();

                obj.lb(useTf) = lb;
                obj.ub(useTf) = ub;
            end
        end

        function useTf = getUseTfForVariable(obj)
            useTf = [obj.varConst obj.varLin obj.varAccel ...
                     obj.varTimeOffset];
        end

        function setUseTfForVariable(obj, useTf)
            obj.varConst = useTf(1);
            obj.varLin = useTf(2);
            obj.varAccel = useTf(3);
            obj.varTimeOffset = useTf(4);
        end

        function updateObjWithVarValue(obj, x)
            ind = 1;

            if(obj.varConst)
                obj.varObj.throttleModel.constTerm = x(ind);
                ind = ind + 1;
            end
            if(obj.varLin)
                obj.varObj.throttleModel.linearTerm = x(ind);
                ind = ind + 1;
            end
            if(obj.varAccel)
                obj.varObj.throttleModel.accelTerm = x(ind);
                ind = ind + 1;
            end

            if(obj.varTimeOffset)
                obj.varObj.throttleModel.tOffset = x(ind);
                ind = ind + 1; %#ok<NASGU>
            end
        end

        function nameStrs = getStrNamesOfVars(obj, evtNum, varLocType)
           if(evtNum > 0)
                subStr = sprintf('Event %i',evtNum);
            else
                subStr = varLocType;
            end

            nameStrs = {sprintf('%s Throttle Constant', subStr), ...
                        sprintf('%s Throttle Rate', subStr), ...
                        sprintf('%s Throttle Acceleration', subStr), ...
                        sprintf('%s Throttle Time Offset', subStr)};

            nameStrs = nameStrs(obj.getUseTfForVariable());
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ c, e ] = polynomial_sort ( o, c, e )

%*****************************************************************************80
%
%% POLYNOMIAL_SORT sorts the information in a polynomial.
%
%  Discussion;
%
%    The coefficients C and exponents E are rearranged so that
%    the elements of E are in ascending order.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 November 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer O, the "order" of the polynomial.
%
%    Input/output, real C(O), the coefficients of the polynomial.
%
%    Input/output, integer E(O), the indices of the exponents of the polynomial.
%
  [ e, i ] = sort ( e );
  c = c(i);

  return
end


text:
function parcel_clusters(clpos_data, clneg_data)
% :Usage:
% ::
%
%     parcel_clusters(clpos_data, clneg_data)
%
% No outputs. Saves all output in separate directory.
%
% First, get eigenvectors for each subject.
%
% We're interested obtaining PARCELS of voxels that tend to co-activate,
% or have the same activation profile.
%
% We can find these by using clustering algorithms to group voxels with
% similar profiles.
%
% Because we have a many voxel x many voxel covariance matrix for each
% subject (lots of data!), it's important to reduce the dimensionality of
% the problem and peform clustering on a REDUCED_DIMENSIONAL space.
%
% We use PCA to do this.  Instead of clustering activation profiles (e.g.,
% time-courses) directly, we cluster eigenvector loadings for each voxel on
% a reduced set of components that explains most of the variance in the
% data.
%
% Similar voxels will have similar loadings across the set of
% eigenvectors.  e.g., two voxels may load high on components [1 3 and 5],
% and low on components [10 and 13].  If they have the same pattern of
% loadings, they should be considered part of the same CLASS.  Groups
% of voxels that are contiguous in space and are members of the same CLASS
% are called parcels.
%
% Images and outputs are saved in their own subdirectory called
% Parcellation_info
%
% The main outputs are:
% parcel_cl             % parcels, one cell per subject, one parcel per
%                       element within cells. same format as clpos_data
%
% parcel_cl_avgs        % parcels, one parcel per element within cells.
%                         same format as clpos_data2
%
% parcel_cl_avgs(x).timeseries contains one cell per subject, with data
% averaged across voxels within that parcel for that subject
% This kind of output is useful, because you can input it directly into
% other mediation analyses.
% [paths, stats2] = mediation(SETUP.data.X, SETUP.data.Y, parcel_cl_avgs(1).timeseries, 'plots', 'verbose', 'names', {'Hi-Low Cue' 'Pain Report' 'Parcel'}, 'boot');
% cluster_orthviews(parcel_cl_avgs(1), {[0 1 0]}, 'add');
% ::
%
%     cd('/Volumes/SCNAlpha/Data_and_Tools/SpeechTask/analysis/wb_multisubject_correl_HR_corrected/mediation_Xprepvsb_Mbrain_Yhr')
%     load cl_b_fdr05_002_01_k3_1_1_prune
%
% then run
%
% There are 2 dimension-reduction steps:
%   1. within-subjects
%   2. is on eigenvectors concatenated across subjects

initial_eigval_limit = 15;      % Number of eigenvalues to save initially for each subject
                                % Need this to reduce computational burden

n_eigs = 7;                     % Number of eigenvectors to use per subject in the clustering algorithm.

% We form a matrix of [voxels x subjects*eigenvectors] (called group_eigenvectors)
% For example, 7 eigenvectors x 10 subjects means clustering will be done
% in a 70 dimensional space.
% This matrix is subjected to a second step of PCA data reduction.  A
% smaller number of scores are saved, and these scores represent sets of
% eigenvector loadings in the original 70 dimensional space.

n_eigs_across = 12;             % Number of dimensions to save out of the original subjects*eigenvectors set
                                % Clustering is done on component scores.
                                % These are canonical eigenvector patterns
                                % that capture regular variations across
                                % subjects. Think of them as compressed
                                % eigenvectors.

% More n_eigs_across  means that clustering will be done in a more complex space.
% Increasing this will tend to split the voxels into smaller parcels, and
% decreasing it will tend to lump them into larger parcels.

n_class_range = [2:20];       % test solutions from this many to this many CLASSES

nclasses = 20;

% The final number of classes you want to use in clustering.
% Increasing this will tend to split the voxels into smaller parcels, and
% decreasing it will tend to lump them into larger parcels.

do_nclasses_search = 0;         % search over n_class_range?

doprune = 1;                    % prune voxels and parcels that are too small or whose voxels don't inter-correlate

mysavedir = 'Parcellation_info';
if ~exist(mysavedir, 'dir'), mkdir(mysavedir), end

create_figure('eigenvalues');

nsubjects = length(clpos_data);
fprintf('Subjects: %3.0f\n', nsubjects);

if ~isempty(clneg_data)
    test_subj_dat = [cat(2, clpos_data{1}(:).all_data) cat(2, clneg_data{1}(:).all_data)];
else
    test_subj_dat = [cat(2, clpos_data{1}(:).all_data)];
end

nvox = size(test_subj_dat, 2);
fprintf('Voxels in mask area: %3.0f\n', nvox);

group_eigenvalues = zeros(nsubjects, initial_eigval_limit);

group_scores = cell(1, nsubjects);

fprintf('PCA: Subject ')

%Clustering of multivariate data is most stable when the data is not sparse,
% i.e., the dimensionality is low relative to the number of observations.
% To limit the dimensionality of the data, a spatio-temporal dimension-reduction
% step is first performed on the [n x v x N] data matrix of  AUC data for
% n trials x v voxels x N participants.  A temporal data reduction is first performed
% to identify components with correlated AUC trial time series within each participant,
% followed by a spatial reduction to identify components with correlated spatial patterns
% across subjects.  First, the [n x v] matrix of AUC data for each participant was
% subjected to PCA, using the [v x v] correlation matrices.  Based on the scree plots
% across subjects, we saved the first [n_eigs] eigenvectors.
% These eigenvectors explained xxx +- xxx% (st. dev. across subjects) of the variance
% in the full dataset. These eigenvectors were scaled by their variances (eigenvalues)
% and concatenated across subjects to form an [v x N*7] matrix of eigenvectors,
% where N=27 subjects.   This matrix was subjected to another (spatial) PCA step to
% identify components with similar spatial maps across participants.
% We retained [n_eigs_across] eigenvectors based on the scree plot, which explained xx% of the
% variance across individuals. This is a data reduction step, and the results are not
% expected to depend strongly on the number of eigenvectors retained at either step,
% as long as most of the variance in the data is explained.

% Temporal reduction
% -------------------------------------------------------------------------
for s = 1:nsubjects

    fprintf('%3.0f', s)

    if ~isempty(clneg_data)
        subj_dat = [cat(2, clpos_data{s}(:).all_data) cat(2, clneg_data{s}(:).all_data)];
    else
        subj_dat = [cat(2, clpos_data{s}(:).all_data)];
    end

    nanvec = any(isnan(subj_dat)) | all(subj_dat == 0);
    wh_bad = find(nanvec);
    if any(wh_bad)
        fprintf('Warning! Subject %3.0f s has missing data (0 or NaN) for these voxels: ', s)
        fprintf('%3.0f ', wh_bad);
        fprintf('\n')

        subj_dat(:, wh_bad) = [];
    end

    % to do PCA on correlation matrix rather than cov
    subj_dat = zscore(subj_dat);

    %[U, eigenvalues, eigenvectors] = svd(subj_dat, 'econ'); % almost, but
    %scaling isn't right, so just use princomp, which does it all

    [eigenvectors, score, eigenvalues] = princomp(subj_dat, 'econ');

    clear subj_dat

    % insert bad voxels back in
    if any(wh_bad)
        for i = 1:size(eigenvectors, 2)
        ev(:, i) = naninsert(nanvec, eigenvectors(:, i));
        end

        ev(isnan(ev)) = 0;
        eigenvectors = ev;
    end

    % The first [initial_eigval_limit]
    group_eigenvalues(s, :) = eigenvalues(1:initial_eigval_limit)';

    % we want to scale the eigenvectors by their variances (eigenvalues),
    % so that components that account for more variation in the data are weighted more heavily.
    eigenvectors = eigenvectors(:, 1:initial_eigval_limit) * diag(eigenvalues(1:initial_eigval_limit));

    group_eigenvectors{s} = eigenvectors;

    plot(group_eigenvalues(s, :), 'ko-');
    drawnow

    % Calculate variance explained for first [initial_eigval_limit]
    % eigenvalues.
    ev = cumsum(group_eigenvalues');
    ev = ev ./ repmat(sum(group_eigenvalues', 1), size(ev, 1), 1);
    evn = ev(n_eigs, :);  % explained variance for these eigs; approximate as it is only proportion of first n eigs
    fprintf('Temporal reduction: Eigenvectors explain %3.2f%% +- %3.2f%% of variance.', 100*mean(evn), 100*std(evn));

end

clear eigenvalues eigenvectors

%%

%n_eigs = input('Enter number of eigenvectors to save: ');

plot_vertical_line(n_eigs);
drawnow

if ~exist(mysavedir, 'dir'), mkdir(mysavedir); end

saveas(gcf, fullfile(mysavedir, 'Eigenvalues'), 'png');

for i = 1:nsubjects
    group_eigenvectors{s}(:, n_eigs + 1 : end) = [];
end

group_eigenvectors = cat(2, group_eigenvectors{:});

%%

%input('Enter range of clusters: ');
%niter = 100;

names = [];

% c = nmdsfig_tools('cluster_solution', [], group_eigenvectors, n_class_range, niter, names);

%% 2nd Dimension redution step : to get scores in lower-dim space for
% clustering

% Spatial reduction
% -------------------------------------------------------------------------

[across_eigenvectors, across_scores, across_eigenvalues] = princomp(group_eigenvectors, 'econ');

create_figure('group_eigenvalues', 1, 3);
plot(across_eigenvalues(1:initial_eigval_limit), 'ko-');
title('Across subjects eigenvalues');

% Calculate variance explained for first [n_eigs_across] eigenvalues.
ev = cumsum(across_eigenvalues);
ev = ev ./ sum(across_eigenvalues);
evn = ev(n_eigs_across, :);  % explained variance for these eigs; approximate as it is only proportion of first n eigs
fprintf('Spatial reduction: Eigenvectors explain %3.2f%% of variance.', 100*evn);

%n_eigs_across = input('Enter number of eigenvectors to save: ');

scores_to_cluster = across_scores(:, 1:n_eigs_across);

plot_vertical_line(n_eigs_across);

subplot(1, 3, 2)
imagesc(scores_to_cluster);

drawnow
%% CLUSTER voxels

classes = [];
sil_vals = {};
mean_sil = [];

if do_nclasses_search

    fprintf('Clustering : ')

    % problem with silhouette is that it really finds natural break-point in
    % data; so looks good for 2 clusters with pos/neg groups in data usually.

    clear s

    for i = n_class_range

        fprintf('%3.0f ', i);

        classes = clusterdata(scores_to_cluster, 'linkage', 'average', 'maxclust', i);
        sil_vals{i} = silhouette(scores_to_cluster, classes);
        mean_sil(i) = mean(sil_vals{i});

    end

    fprintf('\n');

    subplot(1, 3, 3)
    plot(n_class_range, mean_sil(n_class_range), 'ko-', 'MarkerFaceColor', [.2 .6 1], 'LineWidth', 2);

    saveas(gcf, fullfile(mysavedir, 'Group_eigs_and_clustering'), 'png');

end

cl = [clpos_data clneg_data];

disp(['Saving data file: ' mysavedir filesep 'parcellation.mat'])
save(fullfile(mysavedir, 'parcellation'), 'cl', 'n*', '*eig*', '*score*', 'classes', '*sil*')

%% % Get parcels of contiguous regions
% Save averages over voxels for each subject, within each region

fprintf('Getting parcels and associated data: Clustering with %3.0f classes\n', nclasses)

classes = clusterdata(scores_to_cluster, 'linkage', 'average', 'maxclust', nclasses);

clear parcel_cl
disp('Getting contiguous regions: These become parcels');

fprintf('Subject ');
for s = 1:nsubjects

    fprintf('%3.0f', s)

    if ~isempty(clneg_data)
        cl = [clpos_data{s} clneg_data{s}];
    else
        cl = clpos_data{s};
    end

    CLU = clusters2CLU(cl);

    for i = 1:max(classes)

        my_cl = CLU;
        my_cl.XYZmm = my_cl.XYZmm(:, classes == i);
        my_cl.XYZ = my_cl.XYZ(:, classes == i);
        my_cl.Z = my_cl.Z(:, classes == i);

        my_cl.all_data = my_cl.all_data(:, classes == i);

        class_clusters{i} = tor_extract_rois([], my_cl, my_cl);

        [class_clusters{i}.from_class] = deal(i);
    end

    parcel_cl{s} = cat(2, class_clusters{:});

end
fprintf('\n')

% Prune parcels here
if doprune
    [parcel_cl, meanpval, meancor] = prune_parcels(parcel_cl);

    disp('Saving meanpval and meancor for pruned parcels in parcellation.mat')
    save(fullfile(mysavedir, 'parcellation.mat'), '-append', 'meanpval', 'meancor')
end

parcel_cl_avgs = parcel_cl{1};

% another convenient format
for i = 1:length(parcel_cl_avgs)
    parcel_cl_avgs(i).all_data = cell(1, nsubjects);
    parcel_cl_avgs(i).timeseries = cell(1, nsubjects);
    for s = 1:nsubjects
        parcel_cl_avgs(i).all_data{s} = parcel_cl{s}(i).all_data;
        parcel_cl_avgs(i).timeseries{s} = parcel_cl{s}(i).timeseries;
    end
end

fprintf('\n')

disp('Saving parcel_cl and parcel_cl_avgs in parcellation.mat')
save(fullfile(mysavedir, 'parcellation.mat'), '-append', 'class_clusters', 'parcel*')

%volInfo = iimg_read_img('mask.img', 2);

%% RE-do networks on these parcels
% re-define class clusters
% refine class (network) membership
% ---------------------------
[parcel_cl_avgs, NMDS, class_clusters] = parcel_cl_nmds(parcel_cl_avgs);

disp('Saving NMDS structure and final class clusters in parcellation.mat')
save(fullfile(mysavedir, 'parcellation.mat'), '-append', 'NMDS', 'class_clusters')

%% Plots
% ---------------------------
% Plot: data panel
% Orthviews of parcels
% Montages of parcels
parcel_cl_nmds_plots(parcel_cl_avgs, NMDS, 'save', 'savedir', 'Parcellation_info')

end

%% Sub-functions

function [parcel_cl, meanpval, meancor] = prune_parcels(parcel_cl)

    nsubjects = length(parcel_cl);

    %% remove parcels with too few voxels
    vcutoff = 3;
    p_cutoff = .002;

    whomit = false(size(parcel_cl{1}));
    nvox = cat(1, parcel_cl{1}.numVox);
    whomit(nvox < vcutoff) = 1;
    fprintf('Eliminated %3.0f parcels smaller than %3.0f voxels\n', sum(whomit), vcutoff)
    fprintf('Keeping %3.0f parcels\n', sum(~whomit))
    for s = 1:nsubjects
        parcel_cl{s}(whomit) = [];
    end

    nparcels = length(parcel_cl{s});
    [meancor, meanpval, vox_removed] = deal(zeros(nparcels, 1));
    omit_parcels = false(nparcels, 1);

    %
    % Get data across subjects for one parcel
    for p = 1:nparcels

        fprintf('Parcel %3.0f ', p);

        % Get data across subjects for one parcel (p)
        % ---------------------------------------------
        dat = cell(nsubjects, 1);
        for s = 1:nsubjects
            dat{s} = parcel_cl{s}(p).all_data;
        end
        %dat = cat(1, dat{:});

        % Correlate, and get stats
        % ---------------------------------------------
        %[c, pvals] = corrcoef(dat);
        clear c
        for i = 1:nsubjects
            c(:, :, i) = corrcoef(dat{i});
        end

        % Note: p vals will not be correct if there are NaNs!

        mc = nanmean(c, 3);
        se = nanstd(c, 0, 3) ./ sqrt(nsubjects);

        mc = mc .* (1 - eye(size(mc)));
        mc = squareform(mc);

        se = se .* (1 - eye(size(se)));
        se = squareform(se);

        t = mc ./ se;
        pvals = 2 .* (1 - tcdf(abs(t), nsubjects - 1));

        meanpval(p) = mean(pvals);
        meancor(p) = mean(mc);

        pvals = squareform(pvals);

        % Omit the bad (unrelated) voxels
        % ---------------------------------------------
        pv = mean(pvals, 2);

        wh = pv > p_cutoff;

        vox_removed(p) = sum(wh);

        if sum(wh) > length(wh) - vcutoff + 1
            % we have zero or 1 valid voxels left; omit the whole parcel
            omit_parcels(p) = 1;
        elseif any(wh)
            % omit bad voxels

            for s = 1:nsubjects
                parcel_cl{s}(p).XYZmm(:, wh) = [];
                parcel_cl{s}(p).XYZ(:, wh) = [];
                parcel_cl{s}(p).all_data(:, wh) = [];

                parcel_cl{s}(p).timeseries = nanmean(parcel_cl{s}(p).all_data, 2);
                parcel_cl{s}(p).numVox = sum(~wh);

            end

            pvals(wh, :) = [];
            pvals(:, wh) = [];
            mc = squareform(mc);
            mc(wh, :) = [];
            mc(:, wh) = [];

            meanpval(p) = mean(squareform(pvals));
            meancor(p) = mean(squareform(mc));
        end

        fprintf(' mean inter-voxel corr: %3.2f, mean p = %3.4f, vox excluded = %3.0f', meancor(p), meanpval(p), vox_removed(p))
        if omit_parcels(p), fprintf(' OMITTED'); end
        fprintf('\n')

    end

    % remove whole bad parcels
    for s = 1:nsubjects
        parcel_cl{s}(omit_parcels) = [];
    end

    meanpval(omit_parcels) = [];
    meancor(omit_parcels) = [];

    % remove parcels that are now too small
    whomit = false(size(parcel_cl{1}));
    nvox = cat(1, parcel_cl{1}.numVox);
    whomit(nvox < vcutoff) = 1;
    fprintf('Eliminated %3.0f parcels smaller than %3.0f voxels\n', sum(whomit), vcutoff)
    fprintf('Keeping %3.0f parcels\n', sum(~whomit))
    for s = 1:nsubjects
        parcel_cl{s}(whomit) = [];
    end

    nvox = cat(1, parcel_cl{1}.numVox);
    create_figure('Pruned parcels', 2, 1);
    plot(nvox, 'ko', 'MarkerSize', 8, 'MarkerFaceColor', [.5 .5 1]);
    title('Number of voxels in each parcel');
    xlabel('Parcel index number')
    plot_horizontal_line(vcutoff)

    subplot(2, 1, 2)
    plot(meancor, 'ko', 'MarkerSize', 8, 'MarkerFaceColor', [.5 .5 1]);
    title('Mean within-subject correlation value with other voxels in parcel')
    drawnow

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
close all; clear all; clc;
% Create the directory for storing results
spx.fs.ensure_dir('bin');

h = spx.data.motion.Hopkins155;
ne = h.num_examples;
fprintf('Number of sequences: %d\n', ne);

motion3 = 0;
motion2 = 0;
tstart = tic;
% pre-load all examples
h.load_all_examples();
h.describe();
examples = h.get_2_3_motions();
ne = length(examples);
fprintf('Number of 2, 3 motions: %d\n', ne);
for i=1:ne
    fprintf('Example (%d): ', i);
    example = examples{i};
    fprintf('%s, %d motions, %d points, %d frames;', example.name, example.num_motions, example.num_points, example.num_frames);
    switch example.num_motions
        case 2
            motion2 = motion2 + 1;
        case 3
            motion3 = motion3 + 1;
        otherwise
            error('Invalid number of motions.');
    end
    fprintf(' counts: ');
    fprintf('%d ', example.counts);
    [start_indices, end_indices] = spx.cluster.start_end_indices(example.counts);
    % The dataset
    X = example.X;
    fprintf(' effective rank: ');
    for k=1:example.num_motions
        ss = start_indices(k);
        ee = end_indices(k);
        XX = X(:, ss:ee);
        singular_values  = svd(XX, 'econ')';
        r = spx.la.svd.mahdi_rank(singular_values);
        fprintf('%d,  ', r);
    end
    fprintf('\n');

end
elapsed = toc(tstart);
fprintf('Time taken: %.2f seconds \n', elapsed);
fprintf('3 motions: %d\n', motion3);
fprintf('2 motions: %d\n', motion2);

author: indigits

repo: sparse-plex


text:
function V = se3ToVec(se3mat)
% *** CHAPTER 3: RIGID-BODY MOTIONS ***
% Takes se3mat a 4x4 se(3) matrix
% Returns the corresponding 6-vector (representing spatial velocity).
% Example Input:
%
% clear; clc;
% se3mat = [[0, -3, 2, 4]; [3, 0, -1, 5]; [-2, 1, 0, 6]; [0, 0, 0, 0]];
% V = se3ToVec(se3mat)
%
% Output:
% V =
%     1
%     2
%     3
%     4
%     5
%     6

V = [se3mat(3, 2); se3mat(1, 3); se3mat(2, 1); se3mat(1: 3, 4)];
end

author: ShuoYangRobotics

repo: QuadrupedSim


text:
% make. Make an Ot.dll

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 2008-03-18 Created (Johan Meijdam, VORtech)

function make()
   % Omnitrans include path:
   incpath = '"C:\Program Files\Omnitrans International\OmniTRANS 5.1.5\ToolKit"';

   % Debug:
   % eval(['mex -argcheck -g -compatibleArrayDims -output Ot.dll -I' incpath ' Ot.cpp']);
   % Release:
   eval(['mex -O -compatibleArrayDims -output Ot.dll -I' incpath ' Ot.cpp']);

end


text:
function new_name = fname_concat(name, str)

    split = regexp(name, '.mat', 'split');
    new_name = strcat(split{1}, str, '.mat');
end

author: layumi

repo: Image-Text-Embedding


text:
function varargout = exp_eval(x,iters)
% Evaluate the given expression data structure (or short: expression).
% Out-Args = exp_eval(Expression, Iterations)
%
% exp_eval operates on data structures that represent expressions, and evaluates them. The canonical
% expression data structure is of the form:
%  struct
%    .head: function_handle
%    .parts: cell array
%
% This structure represents the expression head(parts{:}). Expression data structures can be nested,
% and can have their value attached, if that value is a struct (they can have a "payload", in a
% sense). In this case, they are structured in the form:
%  struct
%    .whatever_fields_belong_to_the_value
%    .tracking: struct
%      .expression: the actual expression data structure
%
% This way, EEG and STUDY data sets can be regular value types and expressions at the same time,
% which allows to annotate EEGLAB data sets with an expression which (uniquely) identifies them, and
% which reproduces them. Expressions are created by some of the expressions/exp_* functions, and by
% functions using the exp_beginfun/exp_endfun contract (i.e., BCILAB filters and data set editing
% components). For this reason, a call like io_quickload('/data/test.set') or
% flt_resample(eeg,[],200) do by default return an unevaluated expression data structure, which must
% be evaluated first (using exp_eval()), to get access to its value. In other words,
% exp_eval(io_quickload('/data/test.set')) gives the EEGLAB data set. This can be viewed as "lazy"
% evaluation. BCILABs framework functions (bci_train, bci_predict, bci_preproc, etc.) do this
% implicitly and automatically.
%
%
% In:
%   Expression : an expression
%   Iterations : maximum number of repeated evaluations applied (default: 1)
%                for real symbolic computations, pass Inf here
%
% Out:
%   Out-Args   : the result of the evaluation, may still be unevaluated if the expression could not
%                be completely evaluated
%
% Notes:
%   For the purposes of exp_eval(), every data structure that can be expressed in MATLAB is an
%   expression.
%   * Almost all data structures, however, are considered to be 'literal' (and atomic) expressions,
%     such as numbers, strings, function handles, but also recursive data structures, such as cell
%     arrays and structs. These literals cannot be further evaluated. Therefore, almost all MATLAB
%     expressions produce literal expressions, such as, for example, fft(randn(100,1)).
%
%   * A subset of expressions, however, yields by default data that represents unevaluated (or
%     partially evaluated) expressions. Such expressions are formed using either functions that
%     follow the exp_beginfun/exp_endfun protocol, or symbols introduced using exp_symbol or
%     exp_declare. The bulk of these expressions "evaluate" into structures with the two fields
%     'head' and 'parts', where 'head' is a function_handle and 'parts' is a cell array. The purpose
%     of exp_eval() is to fully evaluate these unevaluated expressions, by applying the head to the
%     parts ("calling" the function head with parts as its arguments, in MATLAB terms).
%
%   * To ease interaction with other MATLAB code, some expressions may retain, after having been
%     evaluated (into a MATLAB structure), a subfield that holds the original unevaluated expression
%     that reproduces the data when re-evaluated. These expressions are called "impure" expressions,
%     because they can be viewed both as an unevaluated expression and an evaluated expression. They
%     can thus be processed both by expression-oriented code, and by code that expects values.
%
% Usage Notes:
%   * any MATLAB function can be converted into one that is not immediately evaluated, by adding the
%     two lines "if ~exp_beginfun() return; end", and "exp_endfun;" to its definition, before and
%     after the remaining code.
%   * using exp_symbol, exp_declare_symbols, or @name_of_no_defined_function, arbitrary identifiers
%     can be introduced which are not immediately evaluated (but some of which could be evaluated by
%     exp_eval, given that function code and/or values are defined for them).
%
% Further Reading:
%   A good way to understand symbolic expressions and the operations on them is to read the
%   Mathematica online documentation. The expression system shares some commonalities with
%   Mathematica, with the following considerations:
%    * A subset of functions is implemented with analogous semantics, where SomeFunction is called
%      exp_somefunction in MATLAB
%    * MATLAB has a larger palette of atomic expressions (in addition to String, Number, Symbol, it
%      has Value, which reflects structs, matrices, etc.) In the future, more builtin data structure
%      may be reflected as List, SparseArray, etc.
%    * () are used instead of [], and most operators, such as //., :> cannot be expressed in MATLAB,
%      so that their functional forms must be used.
%    * symbols must be formally introduced before use (with exp_declare_symbols), or expressed
%      either as exp_symbol('x') or as @x
%    * expressions of the form f(a)(b) cannot be interpreted by MATLAB. The workaround is
%      x=f(a);x(b).
%    * expressions are not automatically evaluated; this must be done explicitly with exp_eval().
%      also, the evaluator implemented here is less aggressive, and can be invoked for a limited
%      number of steps.
%    * for the implementation of functions, the exp_beginfun/exp_endfun prologue-epilogue pattern is
%      available.
%
% Examples:
%   % evaluate a dataset-loading expression
%   eeg = exp_eval(io_loadset('/data/projects/test.set'))
%
%   % write up a sequence of filter steps, and finally evalate it; note that, if this particular
%   % sequence had been executed previously, and was cached on disk, the result will be directly
%   % retrieved from the disk cache
%   eeg = io_loadset('/data/projects/test.set')
%   eeg = flt_resample(eeg,200)
%   eeg = flt_iir(eeg,[0.5 2],'highpass')
%   eeg = exp_eval(eeg);
%
% See also:
%   exp_beginfun, exp_endfun, exp_eval_optimized
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-15
dp;

if nargin < 2
    iters = 1;
elseif ~(isnumeric(iters) && isreal(iters) && isscalar(iters) && round(iters) == iters) || iters < 1
    error('The given Iterations argument must be a positive integer.');
end

varargout = {x};
while iters > 0
    % check if we got a canonically represented expression
    if isscalar(x) && all(isfield(x,{'head','parts'}))
        chead = char(x.head);
        % we did, check if it's an application of a lambda or builtin function
        if chead(1) == '@' || (exist(char(x.head),'builtin') && ~exist(char(x.head),'file'))
            % if so, need to pre-evaluate the parts
            x.parts = cellfun(@(v)exp_eval(v,iters),x.parts,'UniformOutput',false);
            % if not all of the args could be evaluated: return unevaluated
            if ~all(cellfun(@is_evaluated,x.parts))
                return; end
        end
        % apply the function
        try
            varargout = hlp_wrapresults(x.head,x.parts{:});
        catch e
            if ~(strcmp(e.identifier,'MATLAB:UndefinedFunction') && strcmp(e.stack(1).name,'hlp_wrapresults'))
                % got a legitimate error, check if we're in dbstop-if-error mode
                settings = dbstatus;
                if any(strcmp({settings.cond},'error')) && ~hlp_resolve('disable_dbstop_if_error_msg',false)
                    % if in dbstop if error mode, we're re-running the line to get the debugger to stop at the right place
                    disp_once('exp_eval: caught error while in "dbstop if error" mode; re-running offending code to get break point...\n  error traceback was: %s\n',hlp_handleerror(e));
                    varargout = hlp_wrapresults(x.head,x.parts{:});
                else
                    rethrow(e);
                end
            else
                % otherwise we just return unevaluated
                return;
            end
        end
    elseif isa(x,'function_handle')
        str = char(x);
        if str(1) ~= '@'
            % regular function handle and undefined?: do a lookup
            if is_undefined_function(x)
                varargout = {hlp_resolve(x,x)}; end
        else
            % lambda function
            try
                % try to obtain the function symbol, if it's a symbolic lambda function
                varargout = {hlp_resolve(get_function_symbol(x),x)};
            catch
            end
        end
    end

    iters = iters-1;

    if iters == 0 || isempty(varargout) || utl_same(x,varargout{1}) || (length(varargout) > 1 && is_evaluated(varargout{1}))
        % evaluation terminates only when the expression no longer changes, or the maximum number of
        % iterations is exceeded (or when multiple evaluated args are being returned and the first
        % one is fully evaluated... - because we want to return these additional outputs)
        break;
    else
        % otherwise we start over...
        x = varargout{1};
    end
end

repo: REST


text:
function [C, S] = QLiftDec2Nevill(X, N, filtername)
%-----------------------------------------------------------------------------
% QLiftDec2Nevill
% Multilevel 2-D decomposition by the lifting scheme and using quincunx grids
%
% Calls for: NevilleR2Q, stencilR2Q, stencilCrop, stencilxgridfRVC, QLmaxlev,
%            storeQ1001, storeR,
%            getcolor01, getcolor10, getcolor00, getcolor11,
%            putcolor01, putcolor10, putcolor00, putcolor11.
% See also: QLiftRec2NV
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: May 16, 2002.
% (c) 1999-2002 Stichting CWI, Amsterdam
%-----------------------------------------------------------------------------
%Firstly, check input data
%
if  isempty(X)
  error(' QLiftDec2Nevill - empty matrix ');
else
  if mod(N, 2) == 1
    error(' QLiftDec2Nevill - only an even number of levels is accepted ');
  end
  if QLmaxlev(size(X), filtername) < N
    error(' QLiftDec2Nevill - too many levels requested ');
  end
  if N < 2
    disp([' QLiftDec2Nevill - WARNING too few levels requested ' ...
          '-> empty decomposition ']);
  end
end
C = []; S = [];
%
%Secondly, some initializing of filters
%
switch lower(filtername)
 case 'neville2'
   [Pa, centerPa] = NevilleR2Q(2);
 case 'neville4'
   [Pa, centerPa] = NevilleR2Q(4);
 case 'neville6'
   [Pa, centerPa] = NevilleR2Q(6);
 case 'neville8'
   [Pa, centerPa] = NevilleR2Q(8);
 otherwise
   error([' QLiftDec2Nevill - unknown filter ' filtername]);
end
%
% Here Pa is the stencil of a prediction step,
% e.g. Pa = 0.250*[0 1 0; 1 0 1; 0 1 0]; centerPa = [2 2];
%
Ua = 0.5 * Pa; centerUa = centerPa;
% Ua is the stencil of the update step and as such determined by Pa.
%
[Pb, centerPb] = stencilR2Q(Pa, centerPa);
[Pb, centerPb] = stencilCrop(Pb, centerPb);
% where Pb is the stencil of a prediction step,
% e.g. Pb = 0.250*[1 0 1; 0 0 0; 1 0 1]; centerPb = [2 2];
%
Ub = 0.5 * Pb; centerUb = centerPb; % Compare to V11!!,
%                                     what about center of stencil??
% Ub is the stencil of the update step and as such determined by Pb.
%
%Thirdly, start decomposition
%
O = X; % For the sake of efficient use of memory this could be improved upon.
% We descend to coarser grids, integer lev indicates number of scale.

for lev=1:2:N
%
   [nO, mO] = size(O);
   if ( nO < 3 ) || ( mO < 3)
     error(' QLiftDec2Nevill - too many levels ');
   else
     sizeO = size(O);
   end
%
%  The Lifting Scheme proceeds from a rectangular grid
%  towards a quincunx grid.
%
%  Stage: predict
%  Convolution of prediction stencil Pa with WHOLE gridfunction O.
%  Note: the operations count is twice as much as it could be.
   PaO = stencilxgridfRVC(O, Pa, centerPa);
%
%  Quincunx grid Q0011 is the union of the values at .00 and .11: "even slots"
%  Quincunx grid Q1001 is the union of the values at .10 and .01: "odd slots"
   Q1001D01 = getcolor01(O) - getcolor01(PaO);
   Q1001D10 = getcolor10(O) - getcolor10(PaO);
   clear PaO;
%  At this point the union (quincunx) of Q1001D01 & Q1001D10
%  contains the DETAILS of O.
%  Note: this is not an "in place" implementation.
%  For the inverse transform Q1001D01 and Q1001D10 have to be stored:
   [C, S] = storeQ1001( Q1001D10, Q1001D01, lev, 'd', C, S);
%
%  Stage: update
   RectDetail = putcolor10(Q1001D10, sizeO) + putcolor01(Q1001D01, sizeO);
   clear Q1001D01 Q1001D10;
   UaD = stencilxgridfRVC(RectDetail, Ua, centerUa);
   clear RectDetail;
%
   Q0011A00 = getcolor00(O) + getcolor00(UaD);
   Q0011A11 = getcolor11(O) + getcolor11(UaD);
   clear UaD O;
%
%  At this point the union (quincunx) of Q0011A00 & Q0011A11
%  contains the updated APPROXIMATION of O, the DETAILS of O
%  were in the union (quincunx) of Q1001D01 & Q1001D10 (see above).
%
%  The Lifting Scheme proceeds by a subsequent step from quincunx
%  to rectangular grid.
%
%  Q0011 is split into the 11 colour with the "odd slots" and
%  the 00 colour with the "even slots".
%
%  Stage: predict
   DETAIL11 = Q0011A11 - ...
      getcolor11(stencilxgridfRVC(putcolor00(Q0011A00, sizeO), Pb, centerPb));
   clear Q0011A11;
%  Note: the operations count is twice as much as it could be
%  DETAIL11 presents the detail gridfunction w.r.t. Q0011
%
%  For the inverse transform DETAIL11 has to be stored:
   [C, S] = storeR( DETAIL11, lev+1, 'd', C, S);
%
%  Stage: update
%  Break up next line in parts in case checking is desired
   APPROX00 = Q0011A00 + ...
      getcolor00(stencilxgridfRVC(putcolor11(DETAIL11, sizeO), Ub, centerUb));
   clear Q0011A00 DETAIL11;
%  APPROX00 now represents the updated version of the approximation of Q0011
%
%  At this point gridfunction DETAIL11 containing the DETAILS has been stored,
%  gridfunction APPROX00 contains the updated APPROXIMATION, on the (down-
%  sampled) rectangular grid and has to be stored as well if at the highest
%  scale.
%  Note that APPROX00 is downsampled onto a rectangular grid with dimensions of
%  half size of the original O.
   if lev+1 >= N
     [C, S] = storeR(APPROX00, lev+1, 'a', C, S);
%    It is obligatory that at least at one scale the Approximation has to be
%    stored or else the scheme cannot be inverted.
     clear APPROX00;
%    In the Lifting Scheme all scales have now been processed!
   else
%    We proceed to the next scale.
     O = APPROX00; clear APPROX00;
   end
end
%-----------------------------------------------------------------------------


text:
%%
% Copyright (c) 2019-present, Mahmoud Afifi
%
% This source code is licensed under the license found in the
% LICENSE file in the root directory of this source tree.
%
% Please, cite the following paper if you use this code:
%
% Mahmoud Afifi and Michael S. Brown. What else can fool deep learning?
% Addressing color constancy errors on deep neural network performance.
% ICCV, 2019
%
% Email: mafifi@eecs.yorku.ca | m.3afifi@gmail.com
%%

addpath('src');
addpath('models');
savepath
disp('Done!');

author: mahmoudnafifi


text:
function labeled_output = PTKWatershedMeyerFromStartingPoints( ~, ~ )
    % PTKWatershedMeyerFromStartingPoints Implementation of Meyer flooding algorithm
    %
    %     This is a Matlab mex file and must be compiled before use.
    %
    %     To compile, type
    %
    %         mex PTKWatershedMeyerFromStartingPoints
    %
    %     in the Matlab command window.

    error('PTKWatershedMeyerFromStartingPoints has not been compiled. You must compile using mex PTKWatershedMeyerFromStartingPoints. Alternatively, use the slower Matlab equivalent function mex PTKWatershedMeyerFromStartingPointsMatlab().');
end

author: tomdoel

repo: pulmonarytoolkit


text:
function fitness = trace_fit_gaussian(C,fr,t_int,fac)

if nargin < 4 || isempty(fac); fac = 3; end
if nargin < 3 || isempty(t_int); t_int = 0.25; end
if nargin < 2 || isempty(fr); fr = 30; end

Np = round(t_int*fr);
[K,T] = size(C);
bas = zeros(K,1);
sn = zeros(K,1);
for i = 1:K
    [~,density,xmesh] = kde(C(i,:));
    [~,ind] = max(density);
    bas(i) = xmesh(ind);
    sn(i) = std(C(i,C(i,:)<bas(i)))/sqrt(1-2/pi);
end

z = bsxfun(@times, 1./(fac*sn(:)), bsxfun(@minus, C, bas));  % normalized z scores
z_pr = 1 - normcdf(z);

filt_z = filter(ones(1,Np),1,log(z_pr),[],2);
fitness = min(filt_z,[],2);

author: flatironinstitute

repo: CaImAn-MATLAB


text:
classdef TestSuperpixelSEEDS
    %TestSuperpixelSEEDS

    methods (Static)
        function test_1
            img = cv.imread(fullfile(mexopencv.root(),'test','fruits.jpg'), ...
                'Color',true, 'ReduceScale',2);
            lab = cv.cvtColor(img, 'RGB2Lab');

            superpix = cv.SuperpixelSEEDS(size(lab), 150, 4);
            superpix.iterate(lab);

            num = superpix.getNumberOfSuperpixels();
            validateattributes(num, {'numeric'}, ...
                {'scalar', 'integer', 'nonnegative'});

            L = superpix.getLabels();
            validateattributes(L, {'int32'}, ...
                {'size',[size(img,1) size(img,2)], '>=',0, '<',num});

            mask = superpix.getLabelContourMask();
            validateattributes(mask, {'logical'}, ...
                {'size',[size(img,1) size(img,2)]});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function mb_lbfgs

% 2010 m.bangert@dkfz.de
% fitting a gaussian distribution to data with a quasi newton method using
% a lbfgs update for the inverse hessian matrix
%
% example : mb_bfgs (function does not take any argument)

% options
visLineSearch = 1;
visBool       = 1;

% sample data from a gaussian distribution
randn('seed',12345); %set the seed of the random number generator - if you do not want to do that comment...
numOfSamples = 250;
data(:,1)    = linspace(-2,4,numOfSamples) + 6/numOfSamples*randn(1,numOfSamples);
data(:,2)    = exp(-(linspace(-2,4,250)-1.5).^2) + .05*randn(1,250);

% plot data
if visBool
    close all
    figure
    hold on
    plot(data(:,1),data(:,2),'b+')
    xlabel('X')
    ylabel('Y')
    visHandle = [];
end

% specify input
mem          = 10;        % number of past gradients and function values used for inverse hessian contruction
x            = NaN*ones(3,mem);
x(:,1)       = [1 1 1]';  % starting value for parameters
objFuncValue = NaN*ones(1,mem);
dx           = NaN*ones(3,mem);

s_k = ones(3,mem-1);
y_k = ones(3,mem-1);
r_k = ones(mem-1,1);
a_k = ones(1,mem-1);

fprintf(['\nFitting a gaussian distribution to data with an L-BFGS-B algorithm...\n\n']);

% 1st calculation of objective function and gradient
objFunc         = @(x) sum((x(3,1).*exp(-(data(:,1)-x(1,1)).^2./x(2,1))-data(:,2)).^2);
objFuncValue(1) = objFunc(x(:,1));
objFuncValue(2) = objFuncValue(1) + 1;
dx(:,1)          = mb_numDiff(objFunc,x(:,1));

% iterate
iter      = 0;
numOfIter = 100;
prec      = 1e-5;

% convergence if gradient smaller than prec, change in objective function
% smaller than prec or maximum number of iteration reached...
while iter < numOfIter && abs((objFuncValue(2)-objFuncValue(1))/objFuncValue(1))>prec && norm(dx(:,1))>prec

    % implementation of L-BFGS according to
    % http://en.wikipedia.org/wiki/L-BFGS

    % inverse hessian update
    q = dx(:,1);
    for i = 1:min(iter,mem-1)
        a_k(i) = r_k(i)*s_k(:,i)'*q;
        q      = q - a_k(i)*y_k(:,i);
    end
    z = s_k(:,1)'*y_k(:,1)/(y_k(:,1)'*y_k(:,1))*q; % this corresponds to H*q where H is approximated as described in Nocedal 9.1
    for i = min(iter,mem-1):-1:1
        b = r_k(i)*y_k(:,i)'*z;
        z = z + s_k(:,i)*(a_k(i)-b);
    end

    % obtain search direction
    dir = -z;

    % 2.1 armijo linesearch to to find acceptable stepsize alpha
    alpha = mb_backtrackingLineSearch(objFunc,objFuncValue(1),x(:,1),dx(:,1),dir);

    % 3
    p = alpha*dir;

        % visualize the line search (optional) used for debugging
        if visLineSearch
            figure
            hold on
            alphas   = linspace(0,2*alpha,100);
            alphaVis = NaN*alphas;
            for i = 1:100
                alphaVis(i) = objFunc(x(:,1)+alphas(i)*dir);
            end
            plot(alphas,alphaVis,'b') % function in direction of line search
            plot(alpha,objFunc(x(:,1)+alpha*dir),'r.') % the steplength found
            plot(alphas,objFuncValue(1) + 1e-1*alphas*(dir'*dx(:,1)),'g') % constraint
            pause(.5)
            close
        end

    % 2.2 update x
    x(:,2:end) = x(:,1:end-1);
    x(:,1)     = x(:,1) + p;
    s_k        = -diff(x,[],2);

    % update objective function (and remember old objective function to
    % check for convergence)
    objFuncValue(2) = objFuncValue(1);
    objFuncValue(1) = objFunc(x(:,1));

    % increment iteration counter
    iter = iter + 1;

    fprintf(1,'Iteration %d: alpha_min=%f, OF=%f\n',iter,alpha,objFuncValue(1));

    % calculate difference of gradients
    dx(:,2:end) = dx(:,1:end-1);
    dx(:,1)     = mb_numDiff(objFunc,x(:,1));
    y_k         = -diff(dx,[],2);

    r_k = 1./diag(y_k'*s_k);

    % update function with new x and plot
    if visBool
        pause(0.5)
        fitFunc = @(t) x(3).*exp(-(t-x(1)).^2./x(2));
        delete(visHandle);
        visHandle = plot(data(:,1),fitFunc(data(:,1)),'r','LineWidth',4);
        drawnow;
    end

end

fprintf(['\n' num2str(iter) ' iteration(s) performed to converge\n'])

fprintf(1,'Final solution: x=(%f,%f,%f)\n',x(1),x(2),x(3));

end


text:
function s = y_to_s_persian ( y )

%*****************************************************************************80
%
%% Y_TO_S_PERSIAN writes a Persian year into a string.
%
%  Format:
%
%    AP YearNumber
%
%  Discussion:
%
%    "AP" stands for "anno Persico".
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year.
%
%    Output, string S, a representation of the year.
%
  s = sprintf ( 'AP %d', y );

  return
end


text:
function disp_(c)
%DISP_        Display of interval gradients in "_" notation
%
%   disp_(c)
%

% written  06/04/09     S.M. Rump
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/05/12     S.M. Rump  complete redesign
%

  INTLAB_GRADIENT_NUMVAR = getappdata(0,'INTLAB_GRADIENT_NUMVAR');

  loose = strcmp(get(0,'FormatSpacing'),'loose');

  numvar = size(c.dx,2);
  if numvar~=INTLAB_GRADIENT_NUMVAR
    warning('**** number of dependent variables and partial derivatives do not coincide')
  end

  name = inputname(1);
  if isempty(name)                    % happens for display(gradientinit(random))
    name = 'ans';
  end

  INTLAB_INTVAL_DISPLAY = getappdata(0,'INTLAB_INTVAL_DISPLAY');
  setappdata(0,'INTLAB_INTVAL_DISPLAY','Display_');
  display(c,name)
  setappdata(0,'INTLAB_INTVAL_DISPLAY',INTLAB_INTVAL_DISPLAY);

repo: OpenMAS