text:
function [y,nbr_bits] = perform_spiht_coding(x,options)

% perform_spiht_coding - SPIHT coding of wavelet coefficients
%
% Coding :
%   options.Jmin = ??;      % minimum scale of the transform
%   options.nb_bits = ??;   % target number of bits
%   [stream,nbr_bits] = perform_spiht_coding(MW,options);
% Decoding :
%   MW = perform_spiht_coding(stream);
%
% This is a simple wrapper of the code of Jing Tian
% <scuteejtian at hotmail.com>

options.null = 0;

if isfield(options, 'arithmetic_coding')
    arithmetic_coding = options.arithmetic_coding;
else
    arithmetic_coding = 1;
end

if size(x,1)>1 && size(x,2)>1
    if isfield(options, 'Jmin')
        Jmin = options.Jmin;
    else
        warning('You should provide options.Jmin');
        Jmin = 4;
    end
    if isfield(options, 'nb_bits')
        nb_bits = options.nb_bits;
    else
        nb_bits = floor( prod(size(x))*0.1 );
    end
    %-----------   Coding   ----------------
    Jmax = log2(size(x,1))-1;
    level = Jmax-Jmin+1;
    y = func_SPIHT_Enc(x, nb_bits, level); y = y(:);
    if arithmetic_coding
        % remove trailing 2
        I = find(y==2); y(I) = [];
        % the 3 first entry are [size nb_bitsplanes level]
        [z,nbr_bits] = perform_arithmetic_coding(y(4:end), +1, options);
        y = [y(1:3); z];
        nbr_bits = nbr_bits + 16; % approx 16 bits
    else
        nbr_bits = nb_bits;
    end
else
    %-----------   Decoding   ----------------
    if arithmetic_coding
        x = x(:);
        % the 3 first entry are [size nb_bitsplanes level]
        z = perform_arithmetic_coding(x(4:end), -1, options);
        x = [x(1:3); z];
        nbr_bits = -1;
    end
    y = func_SPIHT_Dec(x(:)');
end

function out = func_SPIHT_Enc(m, max_bits, level)
% Matlab implementation of SPIHT (without Arithmatic coding stage)
%
% Encoder
%
% input:    m : input image in wavelet domain
%           max_bits : maximum bits can be used
%           level : wavelet decomposition level
%
% output:   out : bit stream
%
% Jing Tian
% Contact me : scuteejtian@hotmail.com
% This program is part of my undergraduate project in GuangZhou, P. R. China.
% April - July 1999

%-----------   Initialization  -----------------
bitctr = 0;
out = 2*ones(1,max_bits - 14);
n_max = floor(log2(abs(max(max(m)'))));
Bits_Header = 0;
Bits_LSP = 0;
Bits_LIP = 0;
Bits_LIS = 0;

%-----------   output bit stream header   ----------------
% image size, number of bit plane, wavelet decomposition level should be
% written as bit stream header.
out(1,[1 2 3]) = [size(m,1) n_max level]; bitctr = bitctr + 24;
index = 4;
Bits_Header = Bits_Header + 24;

%-----------   Initialize LIP, LSP, LIS   ----------------
temp = [];
bandsize = 2.^(log2(size(m, 1)) - level + 1);
temp1 = 1 : bandsize;
for i = 1 : bandsize
    temp = [temp; temp1];
end
LIP(:, 1) = temp(:);
temp = temp';
LIP(:, 2) = temp(:);
LIS(:, 1) = LIP(:, 1);
LIS(:, 2) = LIP(:, 2);
LIS(:, 3) = zeros(length(LIP(:, 1)), 1);
pstart = 1;
pend = bandsize / 2;
for i = 1 : bandsize / 2
    LIS(pstart : pend, :) = [];
    pdel = pend - pstart + 1;
    pstart = pstart + bandsize - pdel;
    pend = pend + bandsize - pdel;
end
LSP = [];

n = n_max;

%-----------   coding   ----------------
while(bitctr < max_bits)

    % Sorting Pass
    LIPtemp = LIP; temp = 0;
    for i = 1:size(LIPtemp,1)
        temp = temp+1;
        if (bitctr + 1) >= max_bits
            if (bitctr < max_bits)
                out(length(out))=[];
            end
            return
        end
        if abs(m(LIPtemp(i,1),LIPtemp(i,2))) >= 2^n % 1: positive; 0: negative
            out(index) = 1; bitctr = bitctr + 1;
            index = index +1; Bits_LIP = Bits_LIP + 1;
            sgn = m(LIPtemp(i,1),LIPtemp(i,2))>=0;
            out(index) = sgn; bitctr = bitctr + 1;
            index = index +1; Bits_LIP = Bits_LIP + 1;
            LSP = [LSP; LIPtemp(i,:)];
            LIP(temp,:) = []; temp = temp - 1;
        else
            out(index) = 0; bitctr = bitctr + 1;
            index = index +1;
            Bits_LIP = Bits_LIP + 1;
        end
    end

    LIStemp = LIS; temp = 0; i = 1;
    while ( i <= size(LIStemp,1))
        temp = temp + 1;
        if LIStemp(i,3) == 0
            if bitctr >= max_bits
                return
            end
            max_d = func_MyDescendant(LIStemp(i,1),LIStemp(i,2),LIStemp(i,3),m);
            if max_d >= 2^n
                out(index) = 1; bitctr = bitctr + 1;
                index = index +1; Bits_LIS = Bits_LIS + 1;
                x = LIStemp(i,1); y = LIStemp(i,2);

                if (bitctr + 1) >= max_bits
                    if (bitctr < max_bits)
                        out(length(out))=[];
                    end
                    return
                end
                if abs(m(2*x-1,2*y-1)) >= 2^n
                    LSP = [LSP; 2*x-1 2*y-1];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    sgn = m(2*x-1,2*y-1)>=0;
                    out(index) = sgn; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                else
                    out(index) = 0; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    LIP = [LIP; 2*x-1 2*y-1];
                end

                if (bitctr + 1) >= max_bits
                    if (bitctr < max_bits)
                        out(length(out))=[];
                    end
                    return
                end
                if abs(m(2*x-1,2*y)) >= 2^n
                    LSP = [LSP; 2*x-1 2*y];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    sgn = m(2*x-1,2*y)>=0;
                    out(index) = sgn; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                else
                    out(index) = 0; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    LIP = [LIP; 2*x-1 2*y];
                end

                if (bitctr + 1) >= max_bits
                    if (bitctr < max_bits)
                        out(length(out))=[];
                    end
                    return
                end
                if abs(m(2*x,2*y-1)) >= 2^n
                    LSP = [LSP; 2*x 2*y-1];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    sgn = m(2*x,2*y-1)>=0;
                    out(index) = sgn; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                else
                    out(index) = 0; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    LIP = [LIP; 2*x 2*y-1];
                end

                if (bitctr + 1) >= max_bits
                    if (bitctr < max_bits)
                        out(length(out))=[];
                    end
                    return
                end
                if abs(m(2*x,2*y)) >= 2^n
                    LSP = [LSP; 2*x 2*y];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    sgn = m(2*x,2*y)>=0;
                    out(index) = sgn; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                else
                    out(index) = 0; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    LIP = [LIP; 2*x 2*y];
                end

                if ((2*(2*x)-1) < size(m) & (2*(2*y)-1) < size(m))
                    LIS = [LIS; LIStemp(i,1) LIStemp(i,2) 1];
                    LIStemp = [LIStemp; LIStemp(i,1) LIStemp(i,2) 1];
                end
                LIS(temp,:) = []; temp = temp-1;

            else
                out(index) = 0; bitctr = bitctr + 1;
                index = index +1; Bits_LIS = Bits_LIS + 1;
            end
        else
            if bitctr >= max_bits
                return
            end
            max_d = func_MyDescendant(LIStemp(i,1),LIStemp(i,2),LIStemp(i,3),m);
            if max_d >= 2^n
                out(index) = 1; bitctr = bitctr + 1;
                index = index +1;
                x = LIStemp(i,1); y = LIStemp(i,2);
                LIS = [LIS; 2*x-1 2*y-1 0; 2*x-1 2*y 0; 2*x 2*y-1 0; 2*x 2*y 0];
                LIStemp = [LIStemp; 2*x-1 2*y-1 0; 2*x-1 2*y 0; 2*x 2*y-1 0; 2*x 2*y 0];
                LIS(temp,:) = []; temp = temp - 1;
            else
                out(index) = 0; bitctr = bitctr + 1;
                index = index +1; Bits_LIS = Bits_LIS + 1;
            end
        end
        i = i+1;
    end

    % Refinement Pass
    temp = 1;
    value = floor(abs(2^(n_max-n+1)*m(LSP(temp,1),LSP(temp,2))));
    while (value >= 2^(n_max+2) & (temp <= size(LSP,1)))
        if bitctr >= max_bits
            return
        end
        s = bitget(value,n_max+2);
        out(index) = s; bitctr = bitctr + 1;
        index = index +1; Bits_LSP = Bits_LSP + 1;
        temp = temp + 1;
        if temp <= size(LSP,1)
            value = floor(abs(2^(n_max-n+1)*m(LSP(temp,1),LSP(temp,2))));
        end
    end

    n = n - 1;
end

function m = func_SPIHT_Dec(in)
% Matlab implementation of SPIHT (without Arithmatic coding stage)
%
% Decoder
%
% input:    in : bit stream
%
% output:   m : reconstructed image in wavelet domain
%
% Jing Tian
% Contact me : scuteejtian@hotmail.com
% This program is part of my undergraduate project in GuangZhou, P. R. China.
% April - July 1999

%-----------   Initialization  -----------------
% image size, number of bit plane, wavelet decomposition level should be
% written as bit stream header.
m = zeros(in(1,1));
n_max = in(1,2);
level = in(1,3);
ctr = 4;

%-----------   Initialize LIP, LSP, LIS   ----------------
temp = [];
bandsize = 2.^(log2(in(1,1)) - level + 1);
temp1 = 1 : bandsize;
for i = 1 : bandsize
    temp = [temp; temp1];
end
LIP(:, 1) = temp(:);
temp = temp';
LIP(:, 2) = temp(:);

LIS(:, 1) = LIP(:, 1);
LIS(:, 2) = LIP(:, 2);
LIS(:, 3) = zeros(length(LIP(:, 1)), 1);
pstart = 1;
pend = bandsize / 2;
for i = 1 : bandsize / 2
    LIS(pstart : pend, :) = [];
    pdel = pend - pstart + 1;
    pstart = pstart + bandsize - pdel;
    pend = pend + bandsize - pdel;
end
LSP = [];

%-----------   coding   ----------------
n = n_max;
while (ctr <= size(in,2))

    %Sorting Pass
    LIPtemp = LIP; temp = 0;
    for i = 1:size(LIPtemp,1)
        temp = temp+1;
        if ctr > size(in,2)
            return
        end
        if in(1,ctr) == 1
            ctr = ctr + 1;
            if in(1,ctr) > 0
                m(LIPtemp(i,1),LIPtemp(i,2)) = 2^n + 2^(n-1);
            else
                m(LIPtemp(i,1),LIPtemp(i,2)) = -2^n  - 2^(n-1);
            end
            LSP = [LSP; LIPtemp(i,:)];
            LIP(temp,:) = []; temp = temp - 1;
        end
        ctr = ctr + 1;
    end

    LIStemp = LIS; temp = 0; i = 1;
    while ( i <= size(LIStemp,1))
        temp = temp + 1;
        if ctr > size(in,2)
            return
        end
        if LIStemp(i,3) == 0
            if in(1,ctr) == 1
                ctr = ctr + 1;
                x = LIStemp(i,1); y = LIStemp(i,2);

                if ctr > size(in,2)
                    return
                end
                if in(1,ctr) == 1
                    LSP = [LSP; 2*x-1 2*y-1];
                    ctr = ctr + 1;
                    if in(1,ctr) == 1
                        m(2*x-1,2*y-1) = 2^n + 2^(n-1);
                    else
                        m(2*x-1,2*y-1) = -2^n  - 2^(n-1);
                    end
                    ctr = ctr + 1;
                else
                    LIP = [LIP; 2*x-1 2*y-1];
                    ctr = ctr + 1;
                end

                if ctr > size(in,2)
                    return
                end
                if in(1,ctr) == 1
                    ctr = ctr + 1;
                    LSP = [LSP; 2*x-1 2*y];
                    if in(1,ctr) == 1;
                        m(2*x-1,2*y) = 2^n + 2^(n-1);
                    else
                        m(2*x-1,2*y) = -2^n  - 2^(n-1);
                    end
                    ctr = ctr + 1;
                else
                    LIP = [LIP; 2*x-1 2*y];
                    ctr = ctr + 1;
                end

                if ctr > size(in,2)
                    return
                end
                if in(1,ctr) == 1
                    ctr = ctr + 1;
                    LSP = [LSP; 2*x 2*y-1];
                    if in(1,ctr) == 1
                        m(2*x,2*y-1) = 2^n + 2^(n-1);
                    else
                        m(2*x,2*y-1) = -2^n  - 2^(n-1);
                    end
                    ctr = ctr + 1;
                else
                    LIP = [LIP; 2*x 2*y-1];
                    ctr = ctr + 1;
                end

                if ctr > size(in,2)
                    return
                end
                if in(1,ctr) == 1
                    ctr = ctr + 1;
                    LSP = [LSP; 2*x 2*y];
                    if in(1,ctr) == 1
                        m(2*x,2*y) = 2^n + 2^(n-1);
                    else
                        m(2*x,2*y) = -2^n  - 2^(n-1);
                    end
                    ctr = ctr + 1;
                else
                    LIP = [LIP; 2*x 2*y];
                    ctr = ctr + 1;
                end

                if ((2*(2*x)-1) < size(m) & (2*(2*y)-1) < size(m))
                    LIS = [LIS; LIStemp(i,1) LIStemp(i,2) 1];
                    LIStemp = [LIStemp; LIStemp(i,1) LIStemp(i,2) 1];
                end
                LIS(temp,:) = []; temp = temp-1;

            else
                ctr = ctr + 1;
            end
        else
            if in(1,ctr) == 1
                x = LIStemp(i,1); y = LIStemp(i,2);
                LIS = [LIS; 2*x-1 2*y-1 0; 2*x-1 2*y 0; 2*x 2*y-1 0; 2*x 2*y 0];
                LIStemp = [LIStemp; 2*x-1 2*y-1 0; 2*x-1 2*y 0; 2*x 2*y-1 0; 2*x 2*y 0];
                LIS(temp,:) = []; temp = temp - 1;
            end
            ctr = ctr + 1;
        end
        i = i+1;
    end

    % Refinement Pass
    temp = 1;
    value = m(LSP(temp,1), LSP(temp,2));
    while (abs(value) >= 2^(n+1) & (temp <= size(LSP,1)))
        if ctr > size(in,2)
            return
        end

        value = value + ((-1)^(in(1,ctr) + 1)) * (2^(n-1))*sign(m(LSP(temp,1),LSP(temp,2)));
        m(LSP(temp,1),LSP(temp,2)) = value;
        ctr = ctr + 1;
        temp = temp + 1;
        if temp <= size(LSP,1)
            value = m(LSP(temp,1),LSP(temp,2));
        end
    end

    n = n-1;
end

function value = func_MyDescendant(i, j, type, m)
% Matlab implementation of SPIHT (without Arithmatic coding stage)
%
% Find the descendant with largest absolute value of pixel (i,j)
%
% input:    i : row coordinate
%           j : column coordinate
%           type : type of descendant
%           m : whole image
%
% output:   value : largest absolute value
%
% Jing Tian
% Contact me : scuteejtian@hotmail.com
% This program is part of my undergraduate project in GuangZhou, P. R. China.
% April - July 1999

s = size(m,1);

S = [];

index = 0; a = 0; b = 0;

while ((2*i-1)<s & (2*j-1)<s)
    a = i-1; b = j-1;

    mind = [2*(a+1)-1:2*(a+2^index)];
    nind = [2*(b+1)-1:2*(b+2^index)];

    chk = mind <= s;
    len = sum(chk);
    if len < length(mind)
        mind(len+1:length(mind)) = [];
    end

    chk = nind <= s;
    len = sum(chk);
    if len < length(nind)
        nind(len+1:length(nind)) = [];
    end

    S = [S reshape(m(mind,nind),1,[])];

    index = index + 1;
    i = 2*a+1; j = 2*b+1;
end

if type == 1
    S(:,1:4) = [];;
end

value = max(abs(S));

author: gpeyre

repo: matlab-toolboxes


text:
function [ p, t ] = cvt_triangle_uniform ( n, sample_num, delaunay_display )

%*****************************************************************************80
%
%% CVT_TRIANGLE_UNIFORM demonstrates how a CVT can be computed and displayed in MATLAB.
%
%  Discussion:
%
%    This simple example carries out an iterative CVT calculation in a
%    triangle, with a uniform density.  The initial placement of the
%    generators is random.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Qiang Du, Vance Faber, Max Gunzburger,
%    Centroidal Voronoi Tessellations: Applications and Algorithms,
%    SIAM Review,
%    Volume 41, 1999, pages 637-676.
%
%  Parameters:
%
%    Input, integer N, the number of generators.
%
%    Input, integer SAMPLE_NUM, the number of sample points.
%
%    Input, logical DELAUNAY_DISPLAY, is TRUE (nonzero) if the Delaunay
%    triangulation is to be displayed.
%
%    Output, real P(N,2), the location of the generators.
%
%    Output, integer T(NT,3), information defining the Delaunay
%    triangulation of the generators.  NT is the number of triangles,
%    which varies depending on the arrangement of the generators.
%

%
%  Here's the triangle we will use.
%
  triangle = [ 0.50, 1.00; ...
               0.00, 0.75; ...
               0.95, 0.00 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_TRIANGLE_UNIFORM:\n' );
  fprintf ( 1, '  A simple demonstration of a CVT computation\n' );
  fprintf ( 1, '  (Centroidal Voronoi Tessellation)\n' );
  fprintf ( 1, '  in a triangle, with a uniform density.\n' );

  if ( nargin < 1 )
    n = 100;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CVT_TRIANGLE_UNIFORM - Note:\n' );
    fprintf ( 1, '  No value of N was supplied.\n' );
    fprintf ( 1, '  N is the number of generators.\n' );
    fprintf ( 1, '  A default value N = %d will be used.\n', n );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  User specified number of generators = %d\n',  n );
  end

  if ( nargin < 2 )
    sample_num = 1000 * n;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CVT_TRIANGLE_UNIFORM - Note:\n' );
    fprintf ( 1, '  No value of SAMPLE_NUM was supplied.\n' );
    fprintf ( 1, '  SAMPLE_NUM is the number of sample points.\n' );
    fprintf ( 1, '  A default value SAMPLE_NUM = %d will be used.\n', ...
      sample_num );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  User specified number of sample points = %d\n', ...
      sample_num );
  end

  if ( nargin < 3 )
    delaunay_display = 0;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CVT_TRIANGLE_UNIFORM - Note:\n' );
    fprintf ( 1, '  No value of DELAUNAY_DISPLAY was supplied.\n' );
    fprintf ( 1, '  DELAUNAY_DISPLAY is TRUE (nonzero) if the\n' );
    fprintf ( 1, '  Delaunay triangulation is also to be displayed.\n' );
    fprintf ( 1, '  A default value DELAUNAY_DISPLAY = %d will be used.\n', ...
      delaunay_display );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  User specified DELAUNAY_DISPLAY = %d\n', ...
      delaunay_display );
  end
%
%  This switch is set to 1 (TRUE) if the ACCUMARRAY command is available.
%  That speeds up the calculation a lot.  If you don't have the ACCUMARRAY
%  command, just set this to 0.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_TRIANGLE_UNIFORM:\n' );
  fprintf ( 1, '  MATLAB''s ACCUMARRAY command can be used for faster\n' );
  fprintf ( 1, '  computation.  This command is not available in\n' );
  fprintf ( 1, '  some versions of MATLAB.  If ACCUMARRAY is available,\n' );
  fprintf ( 1, '  simply make sure that the ACCUMARAY_AVAILABLE variable\n' );
  fprintf ( 1, '  is set to 1!\n' );

  accumarray_available = 1;

  if ( accumarray_available )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The ACCUMARRAY command will be used.\n' );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The ACCUMARRAY command will NOT be used.\n' );
  end
%
%  Clear the figure screen, if already open.
%
  clf
%
%  Randomize the initial locations of the generators.
%
  p(1:2,1:n) = triangle_uniform ( triangle, n );

  it = 0;

  while ( 1 )
%
%  Compute the Delaunay triangle information T for the current nodes.
%
    triangulation = delaunay ( p(1,:), p(2,:) );
%
%  Display the Delaunay triangulation, if requested.
%
    if ( delaunay_display )
      subplot ( 1, 2, 2 )
      trimesh ( triangulation, p(:,1), p(:,2), zeros(n,1) )
      axis ( [ -0.1, 1.1, -0.1, 1.1 ] )
      line ( [ triangle(1,1), triangle(1,2), triangle(1,3), triangle(1,1) ], ...
           [ triangle(2,1), triangle(2,2), triangle(2,3), triangle(2,1) ], ...
        'Color', 'r' );
      title_string = sprintf ( 'Delaunay, step %d', it );
      title ( title_string );
      axis equal
      view ( 2 )
    end
%
%  Display the CVT generators, and the associated Voronoi diagram.
%
    if ( delaunay_display )
      subplot ( 1, 2, 1 )
    end

    voronoi ( p(1,:), p(2,:), triangulation );

    axis ( [ -0.1, 1.1, -0.1, 1.1 ] )
    line ( [ triangle(1,1), triangle(1,2), triangle(1,3), triangle(1,1) ], ...
           [ triangle(2,1), triangle(2,2), triangle(2,3), triangle(2,1) ], ...
      'Color', 'r' );
    title_string = sprintf ( 'Voronoi, step %d', it );
    title ( title_string );
    axis equal
    drawnow
%
%  Generate sample points.
%
%  These sample points implicitly define the geometry of the region.
%  If the region is not a unit square, then the range of the sample
%  data must be changed.
%
%  The data is sampled uniformly.  If a nonuniform density is desired,
%  then the sampling must be done in a biased way.
%
    s = triangle_uniform ( triangle, sample_num );
%
%  For each sample point, find K, the index of the nearest generator.
%  We do this efficiently by using the Delaunay information with
%  Matlab's DSEARCH command, rather than a brute force nearest neighbor
%  computation.
%
    k(1:sample_num,1) = dsearch ( p(1,:), p(2,:), triangulation, s(1,:), s(2,:) );
%
%  The centroid of the Voronoi region associated with each generator
%  is approximated by the average of the sample points it was closest to.
%
    if ( accumarray_available )

      count(1:n) = accumarray ( k, ones(sample_num,1) );
      centroid(1,1:n) = accumarray ( k, s(1,:) );
      centroid(2,1:n) = accumarray ( k, s(2,:) );

    else

      count(1:n) = 0;
      centroid(1,1:n) = 0.0;
      centroid(2,1:n) = 0.0;

      for i = 1 : sample_num
        j = k(i);
        count(j) = count(j) + 1;
        centroid(1,j) = centroid(1,j) + s(1,i);
        centroid(2,j) = centroid(2,j) + s(2,i);
      end

    end
%
%  Replace the generators by the centroids.
%
    p(1,1:n) = ( centroid(1,1:n) ./ count(1:n) )';
    p(2,1:n) = ( centroid(2,1:n) ./ count(1:n) )';

    string = input ( 'RETURN, or Q to quit: ', 's' );

    if ( string == 'q' | string == 'Q' )
        break
    end

    it = it + 1;

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_TRIANGLE_UNIFORM:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function p = triangle_uniform ( t, n )

%*****************************************************************************80
%
%% TRIANGLE_UNIFORM returns random points in a triangle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the vertices of the triangle.
%
%    Input, integer N, the number of points to generate.
%
%    Output, real P(2,N), random points in the triangle.
%
  dim_num = 2;
  alpha = rand ( 1, n );
%
%  Interpret R as a percentage of the triangle's area.
%
%  Imagine a line L, parallel to side 1, so that the area between
%  vertex 1 and line L is R percent of the full triangle's area.
%
%  The line L will intersect sides 2 and 3 at a fraction
%  ALPHA = SQRT ( R ) of the distance from vertex 1 to vertices 2 and 3.
%
  alpha(1:n) = sqrt ( alpha(1:n) );
%
%  Determine the coordinates of the points on sides 2 and 3 intersected
%  by line L.
%
  for dim = 1 : dim_num
    p12(dim,1:n) = ( 1.0 - alpha(1:n) ) * t(dim,1) ...
                         + alpha(1:n)   * t(dim,2);

    p13(dim,1:n) = ( 1.0 - alpha(1:n) ) * t(dim,1) ...
                         + alpha(1:n)   * t(dim,3);
  end
%
%  Now choose, uniformly at random, a point on the line L.
%
  alpha = rand ( 1, n );

  for dim = 1 : dim_num
    p(dim,1:n) = ( 1.0 - alpha(1:n) ) .* p12(dim,1:n) ...
                       + alpha(1:n)   .* p13(dim,1:n);
  end

 return
end

author: johannesgerer

repo: jburkardt-m


text:
function setTableData(hTable, tableData)
%setTableData Summary of this function goes here
%   Detailed explanation goes here

    try
%         jScrollPane = findjobj(hTable);
%         jscroll = jScrollPane.getVerticalScrollBar();
%         sVal = get(jscroll,'value');
        set(hTable,'Data',tableData);
%         drawnow;
%         set(jscroll, 'value', sVal);
%         jScrollPane.repaint;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         jScrollPane = findjobj(hTable);
%         jscroll = jScrollPane.getVerticalScrollBar();
%         sVal = get(jscroll,'value');
%
%         jTable = jScrollPane.getViewport.getView;
%         model = jTable.getModel();
%         if(model.getRowCount() == size(tableData,1))
%             for(i=1:size(tableData,1))
%                 for(j=1:size(tableData,2))
%                     value = tableData{i,j};
%                     if(isnumeric(value))
%                         value = num2str(value);
%                     end
% %                     javaMethodEDT('setValueAt', jTable, value,i-1,j-1);
%                     jTable.setValueAt(value,i-1,j-1);
%                 end
%             end
%             drawnow;
%         else
%             set(hTable,'Data',tableData);
%             drawnow;
%             set(jscroll, 'value', sVal);
%             jScrollPane.repaint;
%         end

    catch ME
        disp(ME.message);
        set(hTable,'Data',tableData);
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [values]=readparm(fname,parm,numval,log_flag)
%READPARM read a parameter from a file
%   [values]=readparm(fname,parm,numval,log_flag)
%   Andy Hooper, June 2011
%
%   ======================================================================
%   09/2006 AH: 'v' option added

%
if nargin<3
    numval=1;
end

if nargin<4
    log_flag=1;
end

fid=fopen(fname);
parms=textscan(fid,'%s');
fclose(fid);

parms=parms{1};
ix=strmatch(parm,parms);
msg=[parms{ix},'='];
for i=1:numval
    parmval=parms{ix+i};
    msg=[msg,parmval,' '];
    if isempty(str2num(parmval))
        values{i}=parmval;
    else
        values{i}=str2num(parmval);
    end
end

if log_flag
    logit(msg)
end

if numval==1
    values=values{1};
end

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function distMat = BagDistMat(data,bagModels)

distMat = zeros(size(data,1),length(bagModels)*length(bagModels{1}),'single');
count = 1;
for mno = 1:length(bagModels)
  curModel = bagModels{mno};
  for j = 1:length(curModel)
    curWk = curModel(j);
    dd = data(:,curWk.dim)*curWk.dir;
    tt = curWk.tr*curWk.dir;
    distMat(:,count) = sign( (dd>tt)-0.5)*curWk.alpha;
    count = count+1;
  end
end

author: kristinbranson

repo: JAABA


text:
function c = ch_cap ( c )

%% CH_CAP capitalizes a single character.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, the character to capitalize.
%
%    Output, character C, the capitalized character.
%
  if ( 'a' <= c & c <= 'z' )
    c = c + 'A' - 'a';
  end

  return
end


text:
function f  = replace_chromosome(intermediate_chromosome, M, V,pop)

%% function f  = replace_chromosome(intermediate_chromosome,pro,pop)
% This function replaces the chromosomes based on rank and crowding
% distance. Initially until the population size is reached each front is
% added one by one until addition of a complete front which results in
% exceeding the population size. At this point the chromosomes in that
% front is added subsequently to the population based on crowding distance.

%  Copyright (c) 2009, Aravind Seshadri
%  All rights reserved.

[N, m] = size(intermediate_chromosome);
f = zeros(pop,m);                                            % modified by zzb
% Get the index for the population sort based on the rank
% [temp,index] = sort(intermediate_chromosome(:,M + V + 1)); % modified by zzb
[~,index] = sort(intermediate_chromosome(:,M + V + 1));      % modified by zzb
% clear temp m;                                              % modified by zzb
sorted_chromosome = zeros(N,m);                              % modified by zzb
% Now sort the individuals based on the index
for i = 1 : N
    sorted_chromosome(i,:) = intermediate_chromosome(index(i),:);
end

% Find the maximum rank in the current population
max_rank = max(intermediate_chromosome(:,M + V + 1));

% Start adding each front based on rank and crowing distance until the
% whole population is filled.
previous_index = 0;
for i = 1 : max_rank
    % Get the index for current rank i.e the last the last element in the
    % sorted_chromosome with rank i.
    % current_index = max(find(sorted_chromosome(:,M + V + 1) == i));      % modified by zzb
    current_index = find(sorted_chromosome(:,M + V + 1) == i, 1, 'last' ); % modified by zzb
    % Check to see if the population is filled if all the individuals with
    % rank i is added to the population.
    if current_index > pop
        % If so then find the number of individuals with current rank i.
        remaining = pop - previous_index;
        % Get information about the individuals in the current rank i.
        temp_pop = sorted_chromosome(previous_index + 1 : current_index, :);
        % Sort the individuals with rank i in the descending order based on
        % the crowding distance.
        % [temp_sort,temp_sort_index] = sort(temp_pop(:, M + V + 2),'descend'); % modified by zzb
        [~,temp_sort_index] = sort(temp_pop(:, M + V + 2),'descend');           % modified by zzb
        % Start filling individuals into the population in descending order
        % until the population is filled.
        for j = 1 : remaining
            f(previous_index + j,:) = temp_pop(temp_sort_index(j),:);
        end
        return;
    elseif current_index < pop
        % Add all the individuals with rank i into the population.
        f(previous_index + 1 : current_index, :) = ...
            sorted_chromosome(previous_index + 1 : current_index, :);
    else
        % Add all the individuals with rank i into the population.
        f(previous_index + 1 : current_index, :) = ...
            sorted_chromosome(previous_index + 1 : current_index, :);
        return;
    end
    % Get the index for the last added individual.
    previous_index = current_index;
end

author: sfvsfv

repo: Mathematical-modeling


text:
function g = ndsimKernGradient(kern, t1, varargin)

% NDSIMKERNGRADIENT Gradient of SIM kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the
% single input motif
% kernel's parameters. As well as the kernel structure and the
% input positions, the user provides a matrix PARTIAL which gives
% the partial derivatives of the function with respect to the
% relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG t : the input locations for which the gradients are being
% computed.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in X.
%
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG t1 : the input locations associated with the rows of the
% kernel matrix.
% ARG t2 : the input locations associated with the columns of the
% kernel matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as X1 and the same number of columns
% as X2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% That is, this function computes sum(sum(deriv(K,a).*partial)) for
% each scalar parameter a, where deriv(K,a) is the matrix-valued
% derivative of the kernel K (computed between times t and t2) with
% respect to the parameter a.
%
% SEEALSO simKernParamInit, kernGradient, simKernDiagGradient, kernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% COPYRIGHT : Jaakko Peltonen, 2011

% KERN

if length(varargin)<2
  t2 = t1;
  covGrad = varargin{1};
else
  t2 = varargin{1};
  covGrad = varargin{2};
end

sigma = sqrt(2/kern.inverseWidth);
if isfield(kern, 'isNegativeS') && (kern.isNegativeS == true)
    variancemultiplier = (kern.sensitivity*kern.sensitivity);
else
    variancemultiplier = kern.variance;
end

dim1 = size(t1, 1);
dim2 = size(t2, 1);
t1Mat = t1(:, ones(1, dim2));
t2Mat = t2(:, ones(1, dim1))';
diffT = (t1Mat - t2Mat);

k1a = sqrt(pi)/2*( t1Mat.*erf(t1Mat/sigma) + t2Mat.*erf(t2Mat/sigma) - diffT.*erf(diffT/sigma) );
k1b = sigma/2*( exp(-(t1Mat/sigma).^2) + exp(-(t2Mat/sigma).^2) - exp(-(diffT/sigma).^2) - 1 );
K = (k1a+k1b)*sigma*variancemultiplier;

% k = ndsimKernCompute(kern, t1, t2);
% deriv_variancemultiplier = sum(sum((k/variancemultiplier).*covGrad));
% deriv_sigma = k/sigma + variancemultiplier*sigma/2*(exp(-(t1Mat/sigma).^2)+exp(-(t2Mat/sigma).^2)-exp(-(diffT/sigma).^2)-1);
% deriv_sigma = sum(sum(deriv_sigma.*covGrad));
% % sigma=sqrt(2/inversewidth)
% % --> df/dinversewidth = (df/dsigma)*(dsigma/dinversewidth)
% % = (df/dsigma)*((-1/2)*sqrt(2)*(inversewidth^(-3/2)))
% % = (df/dsigma)*(-1/2*sigma/inversewidth)
% deriv_inversewidth=deriv_sigma*(-0.5*sigma/kern.inverseWidth);

deriv_variancemultiplier = sum(sum(sigma*(k1a+k1b).*covGrad));
deriv_sigma = (k1a+2*k1b)*variancemultiplier;
% deriv_sigma = (k1a+k1b)*variancemultiplier ...
%     + variancemultiplier*sigma/2*(exp(-(t1Mat/sigma).^2)+exp(-(t2Mat/sigma).^2)-exp(-(diffT/sigma).^2)-1) ...
%     + variancemultiplier/sigma*(-(t1Mat.^2).*exp(-(t1Mat/sigma).^2)-(t2Mat.^2).*exp(-(t2Mat/sigma).^2)+(diffT.^2).*exp(-(diffT/sigma).^2)) ...
%    + variancemultiplier/sigma*(-(t1Mat.^2).*exp(-(t1Mat/sigma).^2)-(t2Mat.^2).*exp(-(t2Mat/sigma).^2)+(diffT.^2).*exp(-(diffT/sigma).^2)) ;

deriv_sigma = sum(sum(deriv_sigma.*covGrad));
deriv_inversewidth=deriv_sigma*(-0.5*sigma/kern.inverseWidth);

if isfield(kern, 'isNegativeS') && (kern.isNegativeS == true)
  % variancemultiplier=sensitivity^2
  % --> df/dsensitivity = (df/dvarmult)*(dvarmult/dsensitivity)
  % = (df/dvarmult)*(2*sensitivity)
  deriv_sensitivity=deriv_variancemultiplier*2*kern.sensitivity;
  g = [deriv_inversewidth deriv_sensitivity];
else
  g = [deriv_inversewidth deriv_variancemultiplier];
end

% gaussianinitial currently unsupported

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [respWave] = Wfm(data, up)
%WFM extracts a respiratory signal corresponding to FM using the continuous
%wavelet transform (Morlet wavelet).
%
%	[respWave] = Wfm(data, up)
%
%	Inputs:
%       data            raw signal data
%       up              universal parameters structure
%
%	Outputs:
%       respWave        a respiratory signal
%

%% Setup
data.t = (1/data.fs)*(1:length(data.v));
data.v = data.v;

%% Downsample to make processing possible
up.paramSet.filt_resample_fs = 25;   % temporarily change downsample freq (perhaps to 50) to increase resolution?
d_s = downsample_data(data, up);
up.paramSet.filt_resample_fs = 25;

%% CWT
% Specify characteristics
s0  = 6/d_s.fs;  % smallest scale
ds = 0.001; % spacing between scales
NbSc = 3000; % number of scales
SCA = {s0,ds,NbSc, 'lin'}; % specify scales
cwtstruct = cwtft({d_s.v, 1/d_s.fs},'scales',SCA);
scales = cwtstruct.scales;
F = scal2frq(scales,cwtstruct.wav,1/d_s.fs); F = d_s.fs./F;
%contour(d_s.t,F,real(cwtstruct.cfs));
%xlabel('Seconds'); ylabel('Pseudo-frequency');
%hold on
mag_mat = cwtstruct.cfs;
rel_rows = F >= up.paramSet.hr_range(1)/60 & F <= up.paramSet.hr_range(2)/60;
rel_mag_mat = mag_mat(rel_rows, :);
[rel_mags, rel_els] = max(rel_mag_mat);
%plot(d_s.t,F(rel_els), 'r', 'LineWidth', 3)
[fm_sig, am_sig] = deal(d_s);
fm_sig.v = F(rel_els);
am_sig.v = abs(rel_mags);

%% Downsample result to usual downsample freq
fm_sig = downsample_data(fm_sig, up);
am_sig = downsample_data(am_sig, up);

%% Store
respWave = fm_sig;

end

author: peterhcharlton

repo: RRest


text:
%% CELL_JOB_ITHACA sets up a job of tasks making animation frames.
%
%  Discussion:
%
%    There are 10 TIF files, constituting a time sequence of images
%    of cells under a microscope.
%
%    Each file is to be processed in such a way that the cells are
%    detected and indicated by a white boundary.
%
%    The processed data is then to be saved in a new TIF file.
%
%    The processing of each file is independent of the others.
%    This means that MATLAB's batch processing facility can be used,
%    involving one job of 10 tasks.
%
%    We assume that this script will be submitted from a desktop machine,
%    for execution on a remote machine.  This raises the question of
%    where the input files are at the beginning, where the output files
%    are at the end, and how the particular task is able to find the
%    input and create the output!
%
%    Unlike the BATCH command, the CREATE_JOB and CREATE_TASK commands
%    do not seem to allow for any way to choose the directory where
%    execution begins.
%
%    One awkward way to try to deal with these problems is as follows:
%    a) "by hand", create a directory on the remote machine, and move the input
%       files there before submitting the job.
%    b) inserting a CD command in the task code that moves to the directory.
%    c) "by hand", collect the processed files from the directory
%       once the tasks are completed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2010
%
%  Author:
%
%    John Burkardt
%
  clear

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CELL_JOB_ITHACA\n' );
  fprintf ( 1, '  Set up and execute a job of independent tasks.\n' );
  fprintf ( 1, '  Each task opens a TIF file that is an image of cells,\n' );
  fprintf ( 1, '  marks the border of each cell in white,\n' );
  fprintf ( 1, '  and writes the modified image to a new TIF file.\n' );
%
%  There are exactly 10 frames.
%
%  We could set this value explicitly:
%
  if ( 0 )

    frame_total = 10;
%
%  but here is a way to count the files automatically.
%
  else
    p = which ( 'AT3_1m4_01.tif' );
    filelist = dir ( [ fileparts(p) filesep 'AT3_1m4_*.tif' ] );
    frame_total = length ( { filelist.name } );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  There are %d files to process.\n', frame_total );
  fprintf ( 1, '  The first file is called          "AT3_1m4_01.tif",\n' );
  fprintf ( 1, '  and its processed version will be "BT3_1m4_01.tif".\n' );
%
%  Define the job.
%  Our image files are on the local machine, so our processing must also
%  be here.
%
  work_directory = '/home/burkardt/public_html/m_src/cell_detection_tasks';

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Define the job:\n' );

  job = createJob ( ...
    'configuration', 'ithaca_2009b', ...
    'FileDependencies', { 'cell_task.m' } );
%
%  Define the tasks that constitute the job.
%  The value of TASK is enough to identify the file to be processed.
%
  task_num = frame_total;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Define %d tasks:\n', task_num );
  fprintf ( 1, '  Working directory is "%s".\n', work_directory );

  for task = 1 : task_num
    task_id = createTask ( job, @cell_task, 0, { work_directory, task } );
  end
%
%  Submit the job.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Submit the job:\n' );

  submit ( job );
%
%  Wait for results.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Wait for job to finish:\n' );

  timestamp ( );
  wait ( job );
  timestamp ( );
%
%  The following commands will capture and print error messages
%  from the tasks, if something bad occurs.
%
  errmsgs = get ( job.Tasks, {'ErrorMessage'} );
  nonempty = ~cellfun ( @isempty, errmsgs );
  celldisp ( errmsgs(nonempty) );
%
%  There are no output arguments from the tasks.  Otherwise, we
%  would gather and inspect them here.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The individual frame files have been created.\n' );
%
%  Clean up by freeing memory associated with the job.
%
  destroy ( job );
%
%  Terminate;
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CELL_JOB_ITHACA:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );


text:
function setup()
    %  setup allows the user to setup the earthquake datafile, overlaying faults, mainshocks

    % turned into function by Celso G Reyes 2017

    report_this_filefun();
    ZG=ZmapGlobal.Data;
    % make the interface
    %

    if ~exist('ZG.bin_dur', 'var')
        ZG.bin_dur = days(14);
    end

    % This is the info window text
    %
    ttlStr='ZMAP Data Import                               ';
    hlpStr1= ...
        ['                                                '
        'For fast and easy data access ZMAP uses the     '
        'internal matlab data format (*.mat) data files).'
        'Earthquakes catalogs as well as additional      '
        'information must first be loaded as and ASCII   '
        'file. Once this data is saved as a matlab *.mat '
        'file it can be reloaded directly into ZMAP.     '
        '                                                '
        'Several Data types can be integrated:           '
        ' EARTHQUAKE CATALOGS: You will need to supply   '
        ' formated ASCII file (e.g hypoellispe format) or'
        ' an ASCII files with the following information  '
        ' in columns sperated by at least one blanck:    '
        '                                                '
        '  lon lat year month day mag depth hour min     '
        '                                                '
        'Example: (California)                           '
        '-116.86 34.35 86 03 27 4.2 15.0 10 25           '
        'Please note the minus sign for W longitudes!    '
        'You chose the magnitude that you would like to  '
        'work with. Zmap works with decimales (0-100) in '
        'the minute position, not with degrees (0-60)!   '];
    hlpStr2= ...
        ['COASTLINE DATA: ZMAP will plot coastlines and/or'
        'state borders on maps. You need to supply an    '
        'ascii datafile with columns seperated by at     '
        'least one blanck:                               '
        'lon  lat  (e.g. -116.86  34.34)                 '
        'A "lift pen" command can be initiade by the line'
        'Inf  Inf (Therefore you avoid connecting islands'
        ' etc.)                                          '
        '                                                '
        'FAULTS DATA: Faults data is imported in the same'
        'way as coastlines.                              '
        'SYMBOLS: Two type of symbols can be displayed on'
        'the seismicity maps: (1) Epicenters of large    '
        'earthquakes as + signs and (2) main faults as   '
        'thick lines. The input format is identical to   '
        'the above - note that you can seperate multiple '
        'main faults by "Inf Inf".                       ' ];
    hlpStr3= ...
        ['The Clear button will remove existing data      '
        'or overlay symbols.                             '
        '                                                ' ];
    % Find out if figure already exists
    %
    loda=findobj('Type','Figure','-and','Name','Import Data into ZMAP');

    % Set up the setup window Enviroment
    %
    if isempty(loda)
        loda = figure_w_normalized_uicontrolunits( ...
            'Name','Import Data into ZMAP',...
            'NumberTitle','off', ...
            'NextPlot','new', ...
            'Visible','off', ...
            'Position',[ (ZG.fipo(3:4) - [500 400]) (ZG.map_len - [100 100]) ]);

        te = text(0.1,0.99,'Load all the available data as ASCII files.\newlinePress <Save> to save the new catalog \newlineafter you have loaded all ASCII files!');
        set(te,'FontSize',ZG.fontsz.m);

        te1 = text(0.40,0.80,'Data ') ;
        set(te1,'FontSize',ZG.fontsz.m,'Color','r','FontWeight','bold')

        te2 = text(0.25,0.58,'Overlay Symbols') ;
        set(te2,'FontSize',ZG.fontsz.m,'Color','r','FontWeight','bold');

        uicontrol('BackGroundColor',[0.8 0.8 0.8]','Style','Pushbutton',...
            'Position',[.01 .65 .30 .08 ],...
            'Units','normalized',...
            'callback',@loaddb_callback,...
            'String','EQ Datafile');

        uicontrol('BackGroundColor',[0.8 0.8 0.8]','Style','Pushbutton',...
            'Position',[.35 .65 .35 .08 ],...
            'Units','normalized',...
            'callback',@loaddb_focal_callback,...
            'String','EQ Datafile (+focal)');

        uicontrol('BackGroundColor',[0.8 0.8 0.8]','Style','Pushbutton',...
            'Position',[.80 .65 .15 .08 ],...
            'Units','normalized',...
            'callback',@clearevents_callback,...
            'String','Clear ');

        uicontrol('BackGroundColor',[0.8 0.8 0.8]','Style','Pushbutton',...
            'Position',[.25 .45 .45 .08 ],...
            'Units','normalized',...
            'callback',@loadmainshocks_callback,...
            'String','Mainshocks ');
        uicontrol('BackGroundColor',[0.8 0.8 0.8]','Style','Pushbutton',...
            'Position',[.75 .45 .15 .08 ],...
            'Units','normalized',...
            'callback',@clearmainshocks_callback,...
            'String','Clear ');

        uicontrol('BackGroundColor',[0.8 0.8 0.8]','Style','Pushbutton',...
            'Position',[.25 .35 .45 .08 ],...
            'Units','normalized',...
            'callback',@loadfaults_callback,...
            'String','Faults ');
        uicontrol('BackGroundColor',[0.8 0.8 0.8]','Style','Pushbutton',...
            'Position',[.75 .35 .15 .08 ],...
            'Units','normalized',...
            'callback',@clearfaults_callback,...
            'String','Clear ');

        uicontrol('BackGroundColor',[0.8 0.8 0.8]','Style','Pushbutton',...
            'Position',[.25 .25 .45 .08 ],...
            'Units','normalized',...
            'callback',@load_mainrupt_callback,...
            'String','Mainrupture');
        uicontrol('BackGroundColor',[0.8 0.8 0.8]','Style','Pushbutton',...
            'Position',[.75 .25 .15 .08 ],...
            'Units','normalized',...
            'callback',@clear_mainrupt_callback,...
            'String','Clear ');

        uicontrol('BackGroundColor',[0.8 0.8 0.8]','Style','Pushbutton',...
            'Position',[.25 .15 .45 .08 ],...
            'Units','normalized',...
            'callback',@loadcoast_callback,...
            'String','Coastline/Borders');
        uicontrol('BackGroundColor',[0.8 0.8 0.8]','Style','Pushbutton',...
            'Position',[.75 .15 .15 .08 ],...
            'Units','normalized',...
            'callback',@clearcoast_callback,...
            'String','Clear ');

        uicontrol('Style','Pushbutton',...
            'Position',[.75 .03 .15 .08 ],...
            'Units','normalized', 'Callback', @close_callback,'String','Cancel');

        uicontrol('Style','Pushbutton',...
            'Position',[.05 .03 .25 .08 ],...
            'Units','normalized',...
            'callback',@save_callback,...
            'String','Save as *.mat');

        uicontrol('Style','Pushbutton',...
            'Position',[.45 .03 .15 .08 ],...
            'Units','normalized',...
            'callback',@info_callback,...
            'String','Info');

        watchoff

    end   %if figure exist
    figure(loda)
    set(gca,'box','off',...
        'SortMethod','childorder','TickDir','out','FontWeight','bold',...
        'visible','off','FontSize',ZG.fontsz.m,'Linewidth',1.2)

    set(loda,'Visible','on');

    function loaddb_callback(mysrc,myevt)
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        dataimpo('noini',1);
    end

    function loaddb_focal_callback(mysrc,myevt)
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        loadasci('focal','of');
    end

    function clearevents_callback(mysrc,myevt)
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        a = [];
        zmap_update_displays();
    end

    function loadmainshocks_callback(mysrc,myevt)
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        loadasci('mainshock','of');
    end

    function clearmainshocks_callback(mysrc,myevt)
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        main  = [];
        zmap_update_displays();
    end

    function loadfaults_callback(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        loadasci('faults','of');
    end

    function clearfaults_callback(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        faults  = [];
        zmap_update_displays();
    end

    function load_mainrupt_callback(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        loadasci('mainfaults','of');
    end

    function clear_mainrupt_callback(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        mainfault  = [];
        zmap_update_displays();
    end

    function loadcoast_callback(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        loadasci('coastline','of');
    end

    function clearcoast_callback(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        coastline  = [];
        zmap_update_displays();
    end

    function close_callback(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        close;
    end

    function save_callback(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        msg.infodisp('  ','Save Data');

        [file1,path1] = uiputfile(fullfile(ZmapGlobal.Data.Directories.data, '*.mat'), 'Filename?');
        save(fullfile(path1, file1),'a','faults','mainfault','coastline','main','infstri');
        close(loda);
        zmap_update_displays();

    end

    function info_callback(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        zmaphelp(ttlStr,hlpStr1,hlpStr2,hlpStr3);
    end

end

author: CelsoReyes


text:
% test for landmark heuristic computation

n = 300;

rep = 'images/landmark_distance/';
if exist(rep)~=7
    mkdir(rep);
end
rep_eps = [rep 'eps/'];
if exist(rep_eps)~=7
    mkdir(rep_eps);
end

name = 'bump';
name = 'map';
name = 'mountain';
name = 'stephanodiscusniagarae';
name = 'cavern';
name = 'gaussian';
name = 'road2';
name = 'constant';

save_images = 0;

nbr_landmarks = 10;
landmark_init = 'rand';
landmark_init = 'farthest';

[M,W] = load_potential_map(name, n);

% just a string representing the number of landmarks points
num_str = num2string_fixeddigit(nbr_landmarks, 3);
marker_size = 20;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% pick starting point
if ~exist('start_points')
    [start_points,end_points] = pick_start_end_point(M);
end

disp('Performing front propagation.');
[D,S] = perform_fast_marching(W, start_points);

% compute landmark
disp('Initializing Landmarks');
switch lower(landmark_init)
    case 'rand'
        landmark = floor( rand(2, nbr_landmarks)*(n-1) )+1;
    case 'farthest'
        landmark = perform_farthest_point_sampling( W, [], nbr_landmarks );
    case 'farthestboundary'
        landmark = perform_farthest_point_sampling_boundary( W, landmark, 1 );
end

% compute distance to landmark
DL = zeros(n,n,nbr_landmarks);
for i=1:nbr_landmarks
    disp('Performing front propagation.');
    [DL(:,:,i),S] = perform_fast_marching(W, landmark(:,i));
end

% compute the approximated distance
[D1,Z] = compute_heuristic_landmark(DL,start_points);

clf;
if ~save_images
    subplot(2,2,1);
end
hold on;
imagesc(W');
plot_scattered(landmark);
axis tight; axis square; axis off;
hold off;
title('Potential with landmarks');
colormap jet(256);
if save_images
    axis image; axis off;
    str = [name '_' num_str 'landmarks_' landmark_init '_potential'];
    saveas(gcf, [rep str '.png'], 'png');
    saveas(gcf, [rep_eps str '.eps'], 'eps');
end

if ~save_images
    subplot(2,2,3);
end
hold on;
imagesc(D');
contour(D', 'k');
axis tight; axis square; axis off;
plot(start_points(1,:), start_points(2,:), 'kx', 'MarkerSize', 10);
hold off;
title('Distance to point');
colormap jet(256);
if save_images
    axis image; axis off;
    str = [name '_' num_str 'landmarks_' landmark_init '_realdist'];
    saveas(gcf, [rep str '.png'], 'png');
    saveas(gcf, [rep_eps str '.eps'], 'eps');
end

if ~save_images
    subplot(2,2,4);
end
hold on;
imagesc(D1');
contour(D1', 'k');
axis tight; axis square; axis off;
plot(start_points(1,:), start_points(2,:), 'kx', 'MarkerSize', 10);
plot(landmark(1,:), landmark(2,:), 'k.', 'MarkerSize', marker_size);
hold off;
title('Approximate distance to point');
colormap jet(256);
if save_images
    axis image; axis off;
    str = [name '_' num_str 'landmarks_' landmark_init '_approxdist'];
    saveas(gcf, [rep str '.png'], 'png');
    saveas(gcf, [rep_eps str '.eps'], 'eps');
end

% str = [rep name '_' num2str(nbr_landmarks) 'landmarks_' landmark_init];
% saveas(gcf, [str '.png'], 'png');

return;

% plot the influence zones
clf;
hold on;
imagesc(Z');
plot_scattered(start_points);
axis off;
hold off;
title('Approximate distance to point');
colormap jet(256);

return;

%%% ANIMATION! %%%
r = 0.3;
nb = 200;
i = 0;
for theta = (0:1/nb:1-1/nb)*2*pi
    i = i+1;
    start_points = 0.5 + r * [cos(theta);sin(theta)];
    start_points = round( start_points*(n-1) )+1;
    [D1,Z] = compute_distance_landmark(DL,start_points);
    clf;
    imagesc(Z);
    axis image;
    axis off;
    colormap jet(256);

    num_str = num2str(i);
    if i<10
        num_str = ['0' num_str];
    end
    if i<100
        num_str = ['0' num_str];
    end
    str = ['anim/' name '_' num_str];
    saveas(gcf, [str '.png'], 'png');
end

author: alecjacobson

repo: gptoolbox


text:
function sources_to_keep = update_sources(qs,sources,addhull_sources);

p = size(sources,2);
% put the sources to add at the end of the list of sources
sources_diff = 1:size(sources,1);
sources_diff(addhull_sources) = [];
sources_to_keep = sources([ sources_diff ],:);
sources_to_add = sources([ addhull_sources ],:);
ltoadd = size(sources_to_add,1);
ltokeep = size(sources_to_keep,1);

while ~isempty(sources_to_add)
    % at each iteration, take one of the sources to add, remove it and add
    % all parents, then check descendants
    rest_of_sources = [sources_to_keep; sources_to_add(2:end,:)];
    local_to_add = sources_to_add(1,:) ;
    toadd  = add_one_source_c(int32(rest_of_sources'),int32(local_to_add),int32(qs));

%     toadd_old = [];
%     for i=1:p
%         temp  = local_to_add;
%         if (temp(i)<qs(i))   % has not reached the end of the dimension
%             temp(i) = temp(i)+1;
%             % check that temp is not a descendant of any sources
%             if ~any(all( rest_of_sources - repmat(temp,ltoadd+ltokeep-1,1)<=0,2))
%                 toadd_old = [ toadd_old;  temp];
%             end
%         end
%     end
%
%     if (size(toadd,1)~=size(toadd_old,1)), keyboard; end
%         toadd
%         toadd_old
%         pause
        sources_to_keep = [ sources_to_keep; toadd ];
        sources_to_add(1,:) = [];
        ltoadd = size(sources_to_add,1);
        ltokeep = size(sources_to_keep,1);
    end

repo: REST


text:
function jed = ymdhms_to_jed_common ( y, m, d, h, n, s )

%*****************************************************************************80
%
%% YMDHMS_TO_JED_COMMON converts a Common YMDHMS date to a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, D, H, N, S, the YMDHMS date.
%
%    Output, real JED, the Julian Ephemeris Date.
%
  [ y1, m1, d1, f1 ] = ymdhms_to_ymdf_common ( y, m, d, h, n, s );

  jed = ymdf_to_jed_common ( y1, m1, d1, f1 );

  return
end


text:
%KLMS Karhunen Loeve Mapping, followed by scaling
%
%   [W,FRAC] = KLMS(A,N)
%   [W,N]    = KLMS(A,FRAC)
%
% INPUT
%   A            Dataset
%   N  or FRAC   Number of dimensions (>= 1) or fraction of variance (< 1)
%                to retain; if > 0, perform PCA; otherwise MCA. Default: N = inf.
%
% OUTPUT
%   W            Affine Karhunen-Loeve mapping
%   FRAC or N    Fraction of variance or number of dimensions retained.
%
% DESCRIPTION
% First a Karhunen Loeve Mapping is performed (i.e. PCA or MCA on the average
% prior-weighted class covariance matrix). The result is scaled by the mean
% class standard deviations. For N and FRAC, see KLM.
%
% Default N: select all ('pre-whiten' the average covariance matrix, i.e.
% orthogonalize and scale). The resulting mapping has a unit average
% covariance matrix.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, KLM, PCA

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: klms.m,v 1.2 2006/03/08 22:06:58 duin Exp $

function [w,truefrac] = klms(a,n)

		if (nargin < 2), n = []; end;
	if (nargin < 1) | (isempty(a))
		w = prmapping('klms',n);
		w = setname(w,'Scaled KL Mapping');
		return
	end

	[w,truefrac] = klm(a,n);        % Calculate KL mapping
	b = a*w;                        % Combine KL mapping with scaling on
	w = w*scalem(b,'c-variance');   % KL-mapped data
	w = setname(w,'Scaled KL Mapping');

	return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [metabolite_structure,IDAdded] = addInfoFromMolFiles(metabolite_structure,folderName,startSearch,endSearch)
% This function creates inchiStrings, smiles, and inchiKeys from provided mol files,
% in the case that these fields are empty (NaN) in the structure.
%
% INPUT
% metabolite_structure  metabolite structure
% folderName            name of folder that contains the mol structures
% startSearch           specify where the search should start in the
%                       metabolite structure. Must be numeric (optional, default: all metabolites
%                       in the structure will be search for)
% endSearch             specify where the search should end in the
%                       metabolite structure. Must be numeric (optional, default: all metabolites
%                       in the structure will be search for)
%
% OUTPUT
% metabolite_structure  Updated metabolite structure
%
%
% Ines Thiele, 09/2021

F = fieldnames(metabolite_structure);
if ~exist('startSearch','var')
    startSearch = 1;
end
if ~exist('endSearch','var')
    endSearch = length(F);
end

if ~exist('folderName','var')
    folderName= 'ctf-main/mets/molFiles';
end
files = dir(folderName);
for i = 1 : length(files)
    fileNames{i,1} = files(i).name;
end
cnt =1;

if exist('fileNames','var')
    for i = startSearch : endSearch
        VMHId = metabolite_structure.(F{i}).VMHId ;

        if isempty(metabolite_structure.(F{i}).inchiString) || length(find(isnan(metabolite_structure.(F{i}).inchiString))) >0 % no inchistring exists

            if strmatch([metabolite_structure.(F{i}).VMHId '.mol'],fileNames,'exact') % but mol file
                % generate inchistring from mol file
                [~,tmp] = system(['obabel ' folderName filesep VMHId  '.mol -o inchi']);
                tmp2=split(tmp,'==');
                inchiString = strtrim(tmp2{1});
                inchiString = regexprep(inchiString,'1 molecule converted','');
                inchiString = regexprep(inchiString,'\n','');
                if ~isempty(inchiString)
                    metabolite_structure.(F{i}).inchiString = inchiString;
                    metabolite_structure.(F{i}).inchiString_source = ['Mol File in ' folderName];
                    IDAdded{cnt,1} = VMHId;
                    IDAdded{cnt,2} = inchiString;
                    IDAdded{cnt,3} = VMHId;cnt = cnt + 1;
                end
                metabolite_structure.(F{i}).hasmolfile =1;
                metabolite_structure.(F{i}).hasmolfile_source = ['Mol File in ' folderName];
            end
        end

        if ~isempty(metabolite_structure.(F{i}).inchiString) && length(find(isnan(metabolite_structure.(F{i}).inchiString))) ==0 % no inchistring exists
            if  isempty(metabolite_structure.(F{i}).inchiKey) || length(find(isnan(metabolite_structure.(F{i}).inchiKey))) >0 % no inchistring exists
                inchiString =metabolite_structure.(F{i}).inchiString;
                [result] = convertInchiString2format(inchiString,'inchiKey');

                if ~isnan(result)

                    metabolite_structure.(F{i}).inchiKey = (result);
                    metabolite_structure.(F{i}).inchiKey_source = ['Mol File in ' folderName];
                    IDAdded{cnt,1} = VMHId;
                    IDAdded{cnt,2} = metabolite_structure.(F{i}).inchiKey;
                    IDAdded{cnt,3} = VMHId;cnt = cnt + 1;
                end
            end
            if  isempty(metabolite_structure.(F{i}).smile) || length(find(isnan(metabolite_structure.(F{i}).smile))) >0 % no inchistring exists

                inchiString =metabolite_structure.(F{i}).inchiString;
                [result] = convertInchiString2format(inchiString,'smiles');

                if ~isnan(result)

                    metabolite_structure.(F{i}).smile_source = ['Mol File in ' folderName];
                    metabolite_structure.(F{i}).smile = (result);
                    IDAdded{cnt,1} = VMHId;
                    IDAdded{cnt,2} =    metabolite_structure.(F{i}).smile;
                    IDAdded{cnt,3} = VMHId;cnt = cnt + 1;
                end
            end
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function ops = convertOpenEphysToRawBInary(ops)

fname       = fullfile(ops.root, sprintf('%s.dat', ops.fbinary));
fidout      = fopen(fname, 'w');
%
clear fs
for j = 1:ops.Nchan
   fs{j} = dir(fullfile(ops.root, sprintf('*CH%d_*.continuous', j) ));
end
nblocks = cellfun(@(x) numel(x), fs);
if numel(unique(nblocks))>1
   error('different number of blocks for different channels!')
end
%
nBlocks     = unique(nblocks);
nSamples    = 1024;  % fixed to 1024 for now!

fid = cell(ops.Nchan, 1);

tic
for k = 1:nBlocks
    for j = 1:ops.Nchan
        fid{j}             = fopen(fullfile(ops.root, fs{j}(k).name));
        % discard header information
        fseek(fid{j}, 1024, 0);
    end
    %
    nsamps = 0;
    flag = 1;
    while 1
        samples = zeros(nSamples * 1000, ops.Nchan, 'int16');
        for j = 1:ops.Nchan
            collectSamps    = zeros(nSamples * 1000, 1, 'int16');

            rawData         = fread(fid{j}, 1000 * (nSamples + 6), '1030*int16', 10, 'b');

            nbatches        = ceil(numel(rawData)/(nSamples+6));
            for s = 1:nbatches
                rawSamps = rawData((s-1) * (nSamples + 6) +6+ [1:nSamples]);
                collectSamps((s-1)*nSamples + [1:nSamples]) = rawSamps;
            end
            samples(:,j)         = collectSamps;
        end

        if nbatches<1000
            flag = 0;
        end
        if flag==0
            samples = samples(1:s*nSamples, :);
        end

        samples         = samples';
        fwrite(fidout, samples, 'int16');

        nsamps = nsamps + size(samples,2);

        if flag==0
            break;
        end
    end
    ops.nSamplesBlocks(k) = nsamps;

    for j = 1:ops.Nchan
       fclose(fid{j});
    end

end

fclose(fidout);

toc

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
% test_mecv2_reconstruction.m -
% image reconstruction using mirror-extended curvelets

if(1)
  N = 512;
  H = N/2;
  [a,b] = ndgrid(0:N-1);
  x = a+b-N;
  x = (x-min(min(x)))/max(max(x-min(min(x))));
  ns = ceil(log2(N) - 3);
  nag = 16;
  ci = 128;

  c = mefcv2(x,N,N,ns,nag);

  cfs = [];
  for s=1:length(c);    for w=1:length(c{s});      cfs = [cfs; abs(c{s}{w}(:))];    end;  end
  cfs = sort(abs(cfs)); cfs = cfs(end:-1:1);
  val = cfs(ci);
  d = c;
  cnt = 0;
  for s=1:length(c)
    for w=1:length(c{s})
      cnt = cnt + sum(sum(abs(c{s}{w})>=val+eps));
      d{s}{w} = c{s}{w}.*(abs(c{s}{w})>=val+eps);
    end
  end
  fprintf(1, 'number of coefficients used = %d\n', cnt);
  y = meicv2(d,N,N,ns,nag);
  figure;  imagesc(real(y)); axis equal; axis tight; colormap gray; colorbar;
  title('Partial reconstruction with mirror-extended curvelets');
end

author: thomaskuestner


text:
function prepareFsetVASARI(pathWORK,pathSTUDY)
% -------------------------------------------------------------------------
% function prepareFsetVASARI(pathWORK,pathSTUDY)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function prepares feature sets for multivariable modeling
% experiments with VASARI features
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathWORK: Full path to the WORKSPACE of the LGG study. Defined in
%              masterScript_LGG.m.
%              --> Ex: '/myProject/WORKSPACE/'
% 2. pathSTUDY: Full path to where the organized data will be saved.
%               Defined in masterScript_LGG.m
%               --> Ex: '/myProject/WORKSPACE/STUDY_DATA'
% -------------------------------------------------------------------------
% OUTPUTS: Organized 'feature sets data in
% '/myProject/WORKSPACE/VASARI/FSET
% -------------------------------------------------------------------------
% AUTHOR(S):
% - Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2017
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015-2017  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;

cd(pathSTUDY), load('outcomes')
nameOutcomes = fieldnames(outcomes); nOutcomes = numel(nameOutcomes);

cd(pathWORK), mkdir('VASARI'), cd('VASARI'), save('outcomes','outcomes'), mkdir('FSET'), cd('FSET'), pathFset = pwd;

for o = 1:nOutcomes
    outcomeName = nameOutcomes{o};
    cd(pathSTUDY), load(['vasari_',outcomeName]) % Variable 'vasari' now in MATLAB workspace
    outcome = outcomes.(outcomeName);
    vasari(isnan(vasari)) = 0; nVas = size(vasari,2); nameVas = cell(nVas,1);
    for i = 1:nVas
        nameVas{i} = ['F',num2str(i)];
    end
    delete = [];
    for i = 1:nVas
        if numel(unique(vasari(:,i))) < 3 % Deleted, not enough variation
            delete = [delete,i];
        end
    end
    vasari(:,delete) = []; nameVas(delete) = [];

    fSet = struct;
    fSet.Data = vasari;
    fSet.Info = nameVas;
    cd(pathFset), save(['FSET_VASARI_',outcomeName],'fSet')
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
%% Reads ECG header in AHA format
% Reads ECG annotations in AHA format. Implements the documentation
% available in the help of the application provided with the database.
%
% Reads ECG headers in AHA format. Implements the documentation
% available in the help of the application provided with the database.
%
% Arguments:
%   + filename: recording to extract header block.
%
% Output:
%   + heasig: header with the ECG properties.
%
%
% See also read_AHA_format, read_AHA_ann, read_ECG, ECGwrapper
%
% Author: Mariano Llamedo Soria
% <matlab:web('mailto:llamedom@electron.frba.utn.edu.ar','-browser') (email)>
% Version: 0.1 beta
% Birthdate: 17/12/2010
% Last update: 19/11/2014
% Copyright 2008-2015
%
function heasig = read_AHA_header(filename)

heasig = [];

%default settings for the whole DB
fidECG = fopen( filename, 'r');

if( fidECG > 0 )

    try
        fseek(fidECG, 0, 'eof');
        bytes_totales = ftell(fidECG);
        %esto es as� en la DB, si fuera necesario se podria cambiar.
        heasig.nsamp = (bytes_totales / 5);
        fclose(fidECG);
    catch ME
        fclose(fidECG);
        rethrow(ME)
    end

    heasig.nsig = 2;
    heasig.freq = 250;
    heasig.gain = repmat(400,1,heasig.nsig);
    heasig.adczero = zeros(1,heasig.nsig);
    % hardcoded fields
    [~, heasig.recname ] = fileparts(filename);
    heasig.btime = '00:00:00';
    heasig.bdate = '01/01/2000';
    heasig.desc = char(strcat(repmat({'ECG'}, heasig.nsig,1), cellstr(num2str(colvec(1:heasig.nsig))) ));

end


text:
function [sourcemodel, cfg] = ft_prepare_leadfield(cfg, data)

% FT_PREPARE_LEADFIELD computes the forward model for many dipole locations
% on a regular 2D or 3D sourcemodel and stores it for efficient inverse modelling
%
% Use as
%   [sourcemodel] = ft_prepare_leadfield(cfg, data)
%
% It is necessary to input the data on which you want to perform the
% inverse computations, since that data generally contain the gradiometer
% information and information about the channels that should be included in
% the forward model computation. The data structure can be either obtained
% from FT_PREPROCESSING, FT_FREQANALYSIS or FT_TIMELOCKANALYSIS. If the data is empty,
% all channels will be included in the forward model.
%
% The configuration should contain
%   cfg.channel            = Nx1 cell-array with selection of channels (default = 'all'),
%                            see FT_CHANNELSELECTION for details
%
% The positions of the sources can be specified as a regular 3-D
% sourcemodel that is aligned with the axes of the head coordinate system
%   cfg.xgrid      = vector (e.g. -20:1:20) or 'auto' (default = 'auto')
%   cfg.ygrid      = vector (e.g. -20:1:20) or 'auto' (default = 'auto')
%   cfg.zgrid      = vector (e.g.   0:1:20) or 'auto' (default = 'auto')
%   cfg.resolution = number (e.g. 1 cm) for automatic sourcemodel generation
% Alternatively the position of a few sources at locations of interest can
% be specified, for example obtained from an anatomical or functional MRI
%   cfg.sourcemodel.pos        = N*3 matrix with position of each source
%   cfg.sourcemodel.inside     = N*1 vector with boolean value whether sourcemodel point is inside brain (optional)
%   cfg.sourcemodel.dim        = [Nx Ny Nz] vector with dimensions in case of 3-D sourcemodel (optional)
%
% The volume conduction model of the head should be specified as
%   cfg.headmodel     = structure with volume conduction model, see FT_PREPARE_HEADMODEL
%
% The EEG or MEG sensor positions can be present in the data or can be specified as
%   cfg.elec          = structure with electrode positions or filename, see FT_READ_SENS
%   cfg.grad          = structure with gradiometer definition or filename, see FT_READ_SENS
%
% Optionally, you can modify the leadfields by reducing the rank (i.e.
% remove the weakest orientation), or by normalizing each column.
%   cfg.reducerank      = 'no', or number (default = 3 for EEG, 2 for MEG)
%   cfg.normalize       = 'yes' or 'no' (default = 'no')
%   cfg.normalizeparam  = depth normalization parameter (default = 0.5)
%   cfg.backproject     = 'yes' or 'no' (default = 'yes') determines when reducerank is applied
%                         whether the lower rank leadfield is projected back onto the original
%                         linear subspace, or not.
%
% Depending on the type of headmodel, some additional options may be
% specified.
%
% For OPENMEEG based headmodels:
%   cfg.openmeeg.batchsize    = scalar (default 100e3), number of dipoles
%                               for which the leadfield is computed in a
%                               single call to the low-level code. Trades off
%                               memory efficiency for speed.
%   cfg.openmeeg.dsm          = 'no'/'yes', reuse existing DSM if provided
%   cfg.openmeeg.keepdsm      = 'no'/'yes', option to retain DSM (no by default)
%   cfg.openmeeg.nonadaptive  = 'no'/'yes'
%
% For SINGLESHELL based headmodels:
%   cfg.singleshell.batchsize = scalar or 'all' (default 1), number of dipoles
%                               for which the leadfield is computed in a
%                               single call to the low-level code. Trades off
%                               memory efficiency for speed.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% If you specify this option the input data will be read from a *.mat
% file on disk. This mat files should contain only a single variable named 'data',
% corresponding to the input structure.
%
% See also FT_SOURCEANALYSIS, FT_DIPOLEFITTING, FT_PREPARE_HEADMODEL, FT_PREPARE_SOURCEMODEL

% Undocumented local options:
% cfg.feedback
% cfg.sel50p      = 'no' (default) or 'yes'
% cfg.lbex        = 'no' (default) or a number that corresponds with the radius
% cfg.mollify     = 'no' (default) or a number that corresponds with the FWHM

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar data
ft_preamble provenance data
ft_preamble trackconfig

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% the data can be passed as input arguments or can be read from disk
hasdata = exist('data', 'var');

if ~hasdata
  % the data variable will be passed to the prepare_headmodel function below
  % where it would be used for channel selection
  data = [];
else
  % check if the input data is valid for this function
  data = ft_checkdata(data);
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed', {'hdmfile', 'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed', {'vol',     'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed', {'grid',    'sourcemodel'});
cfg = ft_checkconfig(cfg, 'renamed', {'om',      'openmeeg'});
cfg = ft_checkconfig(cfg, 'renamed', {'elecfile', 'elec'});
cfg = ft_checkconfig(cfg, 'renamed', {'gradfile', 'grad'});
cfg = ft_checkconfig(cfg, 'renamed', {'optofile', 'opto'});

% set the defaults
cfg.normalize      = ft_getopt(cfg, 'normalize',      'no');
cfg.normalizeparam = ft_getopt(cfg, 'normalizeparam', 0.5);
cfg.lbex           = ft_getopt(cfg, 'lbex',           'no');
cfg.sel50p         = ft_getopt(cfg, 'sel50p',         'no');
cfg.feedback       = ft_getopt(cfg, 'feedback',       'text');
cfg.mollify        = ft_getopt(cfg, 'mollify',        'no');
cfg.patchsvd       = ft_getopt(cfg, 'patchsvd',       'no');
cfg.backproject    = ft_getopt(cfg, 'backproject',    'yes'); % determines whether after rank reduction the subspace projected leadfield is backprojected onto the original space
% cfg.reducerank   = ft_getopt(cfg, 'reducerank', 'no');      % the default for this depends on EEG/MEG and is set below

cfg = ft_checkconfig(cfg, 'renamed', {'tightgrid', 'tight'});  % this is moved to cfg.sourcemodel.tight by the subsequent createsubcfg
cfg = ft_checkconfig(cfg, 'renamed', {'sourceunits', 'unit'}); % this is moved to cfg.sourcemodel.unit by the subsequent createsubcfg

% put the low-level options pertaining to the sourcemodel in their own field
cfg = ft_checkconfig(cfg, 'createsubcfg', {'sourcemodel'});
% move some fields from cfg.sourcemodel back to the top-level configuration
cfg = ft_checkconfig(cfg, 'createtopcfg', {'sourcemodel'});

% this code expects the inside to be represented as a logical array
cfg.sourcemodel = ft_checkconfig(cfg.sourcemodel, 'renamed',  {'pnt' 'pos'});
cfg = ft_checkconfig(cfg, 'inside2logical', 'yes');

if strcmp(cfg.sel50p, 'yes') && strcmp(cfg.lbex, 'yes')
  ft_error('subspace projection with either lbex or sel50p is mutually exclusive');
end

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

% collect and preprocess the electrodes/gradiometer and head model
[headmodel, sens, cfg] = prepare_headmodel(cfg, data);

% set the default for reducing the rank of the leadfields
if ft_senstype(sens, 'eeg')
  cfg.reducerank = ft_getopt(cfg, 'reducerank', 3);
else
  cfg.reducerank = ft_getopt(cfg, 'reducerank', 2);
end

% construct the sourcemodel for which the leadfield will be computed
tmpcfg           = keepfields(cfg, {'sourcemodel', 'mri', 'headshape', 'symmetry', 'smooth', 'threshold', 'spheremesh', 'inwardshift', 'xgrid' 'ygrid', 'zgrid', 'resolution', 'tight', 'warpmni', 'template', 'showcallinfo'});
tmpcfg.headmodel = headmodel;
if ft_senstype(sens, 'eeg')
  tmpcfg.elec = sens;
elseif ft_senstype(sens, 'meg')
  tmpcfg.grad = sens;
end
sourcemodel = ft_prepare_sourcemodel(tmpcfg);

% check whether units are equal (NOTE: this was previously not required,
% this check can be removed if the underlying bug is resolved. See
% http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2387
if ~isfield(headmodel, 'unit') || ~isfield(sourcemodel, 'unit') || ~isfield(sens, 'unit')
  ft_warning('cannot determine the units of all geometric objects required for leadfield computation (headmodel, sourcemodel, sensor configuration). THIS CAN LEAD TO WRONG RESULTS! (refer to http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2387)');
else
  if ~strcmp(headmodel.unit, sourcemodel.unit) || ~strcmp(sourcemodel.unit, sens.unit)
    ft_error('geometric objects (headmodel, sourcemodel, sensor configuration) are not expressed in the same units (this used to be allowed, and will be again in the future, but for now there is a bug which prevents a correct leadfield from being computed; see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2387)');
  end
end

% find the indices of all sourcemodel points that are inside the brain
insideindx = find(sourcemodel.inside);

if ft_headmodeltype(headmodel, 'openmeeg')

  ft_hastoolbox('openmeeg', 1);  % add to path (if not yet on path)

  % repeated system calls to the openmeeg executable makes it rather slow
  % calling it once is much more efficient
  fprintf('calculating leadfield for all positions at once, this may take a while...\n');

  if(~isfield(cfg,'openmeeg'))
    cfg.openmeeg = [];
  end
  batchsize   = ft_getopt(cfg.openmeeg, 'batchsize',100e3); % number of voxels per DSM batch; set to e.g. 1000 if not much RAM available
  dsm         = ft_getopt(cfg.openmeeg, 'dsm');             % reuse existing DSM if provided
  keepdsm     = ft_getopt(cfg.openmeeg, 'keepdsm', 'no');   % retain DSM
  nonadaptive = ft_getopt(cfg.openmeeg, 'nonadaptive', 'no');

  ndip       = length(insideindx);
  numchunks  = ceil(ndip/batchsize);
  if(numchunks > 1)
    if istrue(keepdsm)
      ft_warning('Keeping DSM output not supported when the computation is split into batches')
    end
    keepdsm = false;
  end

  try
    % DSM computation is computationally intensive:
    % As it can be reused with same voxel sourcemodel (i.e. if voxels are defined in
    % MRI coordinates rather than MEG coordinates), optionally save result.
    % Dense voxel grids may require several gigabytes of RAM, so optionally
    % split into smaller batches

    [h2sens,ds2sens] = ft_sensinterp_openmeeg(sourcemodel.pos(insideindx,:), headmodel, sens);

    % use pre-existing DSM if present
    if(~isempty(dsm))
      lf = ds2sens + h2sens*headmodel.mat*dsm;
    else
      lf = zeros(size(ds2sens)); % pre-allocate Msensors x Nvoxels

      for ii = 1:numchunks
        % select sourcemodel positions for this batch
        diprange = (((ii-1)*batchsize + 1):(min((ii)*batchsize,ndip)));
        % remap with 3 orientations per position
        diprangeori = [((ii-1)*3*batchsize + 1):(min((ii)*3*batchsize,3*ndip))];
        dsm = ft_sysmat_openmeeg(sourcemodel.pos(insideindx(diprange),:), headmodel, sens, nonadaptive);
        lf(:,diprangeori) = ds2sens(:,diprangeori) + h2sens*headmodel.mat*dsm;

        if istrue(keepdsm)
          % retain DSM in cfg if desired
          cfg.openmeeg.dsm = dsm;
        end

        dipindx = insideindx(diprange);
      end
    end
  catch
    me = lasterror;
    rethrow(me);
  end

  % apply montage, if applicable
  if isfield(sens, 'tra')
    lf = sens.tra * lf;
  end

  % lead field computation already done, but pass to ft_compute_leadfield so that
  % any post-computation options can be applied (e.g., normalization, etc.)
  lf = ft_compute_leadfield(sourcemodel.pos(diprange,:), sens, headmodel, 'lf', lf, 'reducerank', cfg.reducerank, 'normalize', cfg.normalize, 'normalizeparam', cfg.normalizeparam, 'backproject', cfg.backproject);

  % reshape result into sourcemodel.leadfield cell-array
  for i=1:ndip
    sourcemodel.leadfield{insideindx(i)} = lf(:,3*(i-1) + [1:3]);
  end
  clear lf

elseif ft_headmodeltype(headmodel, 'singleshell')
  cfg.singleshell = ft_getopt(cfg, 'singleshell', []);
  batchsize       = ft_getopt(cfg.singleshell, 'batchsize', 1);
  if ischar(batchsize) && strcmp(batchsize, 'all')
    batchsize = length(insideindx);
  end

  dippos     = sourcemodel.pos(insideindx,:);
  ndip       = length(insideindx);
  numchunks  = ceil(ndip/batchsize);

  ft_progress('init', cfg.feedback, 'computing leadfield');
  for k = 1:numchunks
    ft_progress(k/numchunks, 'computing leadfield %d/%d\n', k, numchunks);
    diprange = (((k-1)*batchsize + 1):(min(k*batchsize,ndip)));
    tmp      = ft_compute_leadfield(dippos(diprange,:), sens, headmodel, 'reducerank', cfg.reducerank, 'normalize', cfg.normalize, 'normalizeparam', cfg.normalizeparam, 'backproject', cfg.backproject);
    for i=1:length(diprange)
      thisindx = insideindx(diprange(i));
      if istrue(cfg.backproject)
        sourcemodel.leadfield{thisindx} = tmp(:,(i-1)*3+(1:3));
      else
        sourcemodel.leadfield{thisindx} = tmp(:,(i-1)*cfg.reducerank+(1:cfg.reducerank));
      end

      if isfield(cfg, 'sourcemodel') && isfield(cfg.sourcemodel, 'mom')
        % multiply with the normalized dipole moment to get the leadfield in the desired orientation
        sourcemodel.leadfield{thisindx} = sourcemodel.leadfield{thisindx} * sourcemodel.mom(:,thisindx);
      end
    end
  end
  ft_progress('close');

else
  ft_progress('init', cfg.feedback, 'computing leadfield');
  for i=1:length(insideindx)
    % compute the leadfield on all sourcemodel positions inside the brain
    ft_progress(i/length(insideindx), 'computing leadfield %d/%d\n', i, length(insideindx));
    thisindx = insideindx(i);
    sourcemodel.leadfield{thisindx} = ft_compute_leadfield(sourcemodel.pos(thisindx,:), sens, headmodel, 'reducerank', cfg.reducerank, 'normalize', cfg.normalize, 'normalizeparam', cfg.normalizeparam, 'backproject', cfg.backproject);

    if isfield(cfg, 'sourcemodel') && isfield(cfg.sourcemodel, 'mom')
      % multiply with the normalized dipole moment to get the leadfield in the desired orientation
      sourcemodel.leadfield{thisindx} = sourcemodel.leadfield{thisindx} * sourcemodel.mom(:,thisindx);
    end
  end % for all sourcemodel locations inside the brain
  ft_progress('close');
end

% represent the leadfield for positions outside the brain as empty array
sourcemodel.leadfield(~sourcemodel.inside) = {[]};

% add the label of the channels
sourcemodel.label           = sens.label;
sourcemodel.leadfielddimord = '{pos}_chan_ori';

% mollify the leadfields
if ~strcmp(cfg.mollify, 'no')
  sourcemodel = mollify(cfg, sourcemodel);
end

% combine leadfields in patches and do an SVD on them
if ~strcmp(cfg.patchsvd, 'no')
  sourcemodel = patchsvd(cfg, sourcemodel);
end

% compute the 50 percent channel selection subspace projection
if ~strcmp(cfg.sel50p, 'no')
  sourcemodel = sel50p(cfg, sourcemodel, sens);
end

% compute the local basis function expansion (LBEX) subspace projection
if ~strcmp(cfg.lbex, 'no')
  sourcemodel = lbex(cfg, sourcemodel);
end

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

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble previous   data
ft_postamble provenance sourcemodel
ft_postamble history    sourcemodel
ft_postamble savevar    sourcemodel

author: spm


text:
function residual ( )

%*****************************************************************************80
%
%% RESIDUAL assembles the residual vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 April 2006
%
  global area
  global beta
  global difeta
  global fzero
  global indx
  global kappa
  global lambda1
  global lambda2
  global m
  global nel
  global neqn
  global nnodes
  global node
  global nq
  global nu1
  global nu2
  global nunk
  global xc
  global xl
  global xr

  resid = zeros(nunk*neqn,1 );
  xq = zeros(3);
  wq = zeros(3);
%
%  Set up quadrature information on [-1,1]
%
  [ xq, wq ] = gauss3pt;

  for it = 1 : nel

    for iq = 1 : nq

      x = xc(node(it,1)) + area(it) / 2.0 * ( xq(iq) + 1.0 );
      ar = area(it) * wq(iq);
%
%  Evaluate solutions at quadrature points.
%
      [ vh, vhx ] = eval_pt ( x, it, vcur );
      [ ch, chx ] = eval_pt ( x, it, ccur );
      [ fh, fhx ] = eval_pt ( x, it, fcur );
      [ etah, etahx ] = eval_pt ( x, it, etacur );
      [ vhold, vhxold ] = eval_pt ( x, it, vold );
      [ chold, chxold ] = eval_pt ( x, it, cold );
      [ fhold, fhxold ] = eval_pt ( x, it, fold );
      [ etahold, etahxold ] = eval_pt ( x, it, etaold );

      tauxovertau = tau_prob ( fh, fhx, ch, chx );

      for nuk = 1 : 4

        for in = 1 : nnodes

          ip = node(it,in);
          i = indx(ip);

          if ( 0 < i )

            [ bb, bx ] = quadbf ( x, it, in, xc, node );

            term_linear = 0;
            term_nonlin = 0;

            if ( nuk == 1 )
              term_time = ( vh - vhold ) / dt;
              term_nonlin = ( lambda1 * vh * etah ) / ( 1.0 + nu1 * vh );
            elseif ( nuk == 2 )
              term_time = ( ch - chold ) / dt;
              term_nonlin = -lambda1 * vh * etah / ( 1.0 + nu1 * vh );
            elseif ( nuk == 3 )
              term_time = ( fh - fhold ) / dt;
              term_nonlin = -beta * fh * ( fzero - fh ) * etah + ...
                lambda2 * ch * fh / ( 1.0 + nu2 * ch);
            elseif ( nuk == 4 )
              term_time = ( etah - etahold ) / dt;
              term_linear = difeta * bx * etah * ( etahx / etah - tauxovertau );
            end

            k = i + ( nuk - 1 ) * nunk;
            [rhs_source] = rhsfun_source ( x, time, bb, nuk );

            resid(k) = resid(k) + ( rhs_source  ...
              - ( term_time + term_nonlin ) * bb - term_linear ) * ar;

          end

        end

      end

    end

  end
%
%  Calculate norm of residual vector.
%
  resid_norm = sqrt ( resid' * resid ) / ( neqn * nunk )

  return
end


text:
function msh = meshCreate(mshType)
% mrMesh creation routine
%
%   msh = meshCreate(mshType);
%
% We only create a vistaMesh type.  In the future we may design additional
% mesh structures.  See notes below about the properties of the msh fields.
%
% See also:  meshSet/Get  and mrmSet/Get
%
% Example:
%    msh = meshCreate;
%    msh = meshCreate('vista mesh');
%
% Stanford VISTA team

if ieNotDefined('mshType'), mshType = 'vistaMesh'; end

mshType = mrvParamFormat(mshType);

switch lower(mshType)
    case 'vistamesh'
        msh = vistaMeshCreate;
    otherwise
        error('Unknown mesh type %s\n',mshType);
end

return;

%------------------------------
function msh = vistaMeshCreate
%Default settings for a mrVista mesh
%
%   msh = vistaMeshCreate;
%
%  * The triangles are triplets of vertices.
%  * The vertices are numbered from [0,n-1], consistent with C numbering,
%  but not Matlab numbering.  This is necessary to work with mrMeshSrv.
%
%
% fields = {'name', 'host', 'id', 'actor', 'mmPerVox', 'lights', 'origin', ...
%     'initialvertices', 'vertices', 'triangles', 'colors', 'normals', 'curvature',...
%     'ngraylayers', 'vertexGrayMap', 'fibers',...
%     'smooth_sinc_method', 'smooth_relaxation', 'smooth_iterations', 'mod_depth'};
%
% (c) Stanford VISTA Team

msh.name = '';
msh.type = 'vistaMesh';
msh.host = 'localhost';
msh.id   = -1;    % Figure this out as soon as possible.
msh.filename = [];
msh.path = [];
msh.actor = 33;   % This is default.  But we should be able to change.
msh.mmPerVox = [1 1 1];
msh.lights   = {};
msh.origin   = [];
msh.initVertices = [];   % Initial vertices, without smoothing.
msh.vertices = [];
msh.triangles = [];

% Surface shading typically for pseudo-color or for showing curvature
msh.colors = [];
msh.mod_depth = 0.25;

msh.normals = [];   % Normals to the patches.  Can be computed using Matlab

msh.curvature = []; % Not sure how we compute this.  Uh oh.

% Not sure why these are here.  They relate the vertices to gray map.
% Probably essential for VISTASOFT/mrBOLD
msh.grayLayers = [];
msh.vertexGrayMap =[];

% Not sure why these are here too, but probably related to mrDiffusion
msh.fibers = [];

% Mesh smoothing related
msh.smooth_sinc_method = 0;
msh.smooth_relaxation = 0.5;
msh.smooth_iterations = 32;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function AnimateTrajectory(x,y,z)

ax = newplot([]);
axis(ax,[min(x) max(x) min(y) max(y) min(z) max(z)])
%axis(ax,[0 max(x) 0 max(y) 0 max(z)])
head = line('parent',ax,'color','r','marker','.','markersize',15,'erase','xor', ...
    'xdata',x(1),'ydata',y(1),'zdata',z(1));
body = line('parent',ax,'color','b','linestyle','-','erase','none', ...
    'xdata',[],'ydata',[],'zdata',[]);
box on
%set(gca,'yDir','reverse')
% Grow the body
for i = 1:length(x)
    pause(.01)
    set(head,'xdata',x(i),'ydata',y(i),'zdata',z(i))
    set(body,'xdata',x(1:i),'ydata',y(1:i),'zdata',z(1:i))
end

author: Sable

repo: mcbench-benchmarks


text:
function exact = p00_exact ( problem, dim_num )

%*****************************************************************************80
%
%% P00_EXACT returns the exact integral for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the number of the desired test problem.
%
%    Input, integer DIM_NUM, the dimension of the problem.
%
%    Output, real EXACT, the exact value of the integral.
%
  if ( problem == 1 )
    exact = p01_exact ( dim_num );
  elseif ( problem == 2 )
    exact = p02_exact ( dim_num );
  elseif ( problem == 3 )
    exact = p03_exact ( dim_num );
  elseif ( problem == 4 )
    exact = p04_exact ( dim_num );
  elseif ( problem == 5 )
    exact = p05_exact ( dim_num );
  elseif ( problem == 6 )
    exact = p06_exact ( dim_num );
  elseif ( problem == 7 )
    exact = p07_exact ( dim_num );
  elseif ( problem == 8 )
    exact = p08_exact ( dim_num );
  elseif ( problem == 9 )
    exact = p09_exact ( dim_num );
  elseif ( problem == 10 )
    exact = p10_exact ( dim_num );
  elseif ( problem == 11 )
    exact = p11_exact ( dim_num );
  elseif ( problem == 12 )
    exact = p12_exact ( dim_num );
  elseif ( problem == 13 )
    exact = p13_exact ( dim_num );
  elseif ( problem == 14 )
    exact = p14_exact ( dim_num );
  elseif ( problem == 15 )
    exact = p15_exact ( dim_num );
  elseif ( problem == 16 )
    exact = p16_exact ( dim_num );
  elseif ( problem == 17 )
    exact = p17_exact ( dim_num );
  elseif ( problem == 18 )
    exact = p18_exact ( dim_num );
  elseif ( problem == 19 )
    exact = p19_exact ( dim_num );
  elseif ( problem == 20 )
    exact = p20_exact ( dim_num );
  elseif ( problem == 21 )
    exact = p21_exact ( dim_num );
  elseif ( problem == 22 )
    exact = p22_exact ( dim_num );
  elseif ( problem == 23 )
    exact = p23_exact ( dim_num );
  elseif ( problem == 24 )
    exact = p24_exact ( dim_num );
  elseif ( problem == 25 )
    exact = p25_exact ( dim_num );
  elseif ( problem == 26 )
    exact = p26_exact ( dim_num );
  elseif ( problem == 27 )
    exact = p27_exact ( dim_num );
  elseif ( problem == 28 )
    exact = p28_exact ( dim_num );
  elseif ( problem == 29 )
    exact = p29_exact ( dim_num );
  elseif ( problem == 30 )
    exact = p30_exact ( dim_num );
  elseif ( problem == 31 )
    exact = p31_exact ( dim_num );
  elseif ( problem == 32 )
    exact = p32_exact ( dim_num );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_EXACT - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem number = %d\n', problem );
    error ( 'P00_EXACT - Fatal error!' );
  end

  return
end

author: johannesgerer


text:
function [H] = globalrescale(f);

% GLOBALRESCALE creates the homogenous spatial transformation matrix
% for a 7 parameter rigid-body transformation with global rescaling
%
% Use as
%   [H] = globalrescale(f)
%
% The transformation vector f should contain the
%   x-shift
%   y-shift
%   z-shift
% followed by the
%   pitch (rotation around x-axis)
%   roll  (rotation around y-axis)
%   yaw   (rotation around z-axis)
% followed by the
%   global rescaling factor

% Copyright (C) 2000-2005, Robert Oostenveld
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% $Log: globalrescale.m,v $
% Revision 1.1  2009/01/30 04:02:06  arno
% *** empty log message ***
%
% Revision 1.3  2005/08/15 08:15:32  roboos
% reimplemented the rotate function, which contained an error (the error is in the AIR technical reference)
% changed all functions to be dependent on the rotate, translate and scale function
% all functions now behave consistenly, which also means that they are not compleetly backward compatible w.r.t. the order of the rotations
%
% Revision 1.2  2004/05/19 09:57:07  roberto
% added GPL copyright statement, added CVS log item
%

% compute the homogenous transformation matrix for the translation
T = translate(f([1 2 3]));

% compute the homogenous transformation matrix for the rotation
R = rotate_new(f([4 5 6]));

% compute the homogenous transformation matrix for the global scaling
S = scale(f([7 7 7]));

% compute the homogenous transformation matrix for the combination
H = T*R*S;