text:
function M = load_image(type, n, options)

% load_image - load benchmark images.
%
%   M = load_image(name, n, options);
%
%   name can be:
%   Synthetic images:
%       'chessboard1', 'chessboard', 'square', 'squareregular', 'disk', 'diskregular', 'quaterdisk', '3contours', 'line',
%       'line_vertical', 'line_horizontal', 'line_diagonal', 'line_circle',
%       'parabola', 'sin', 'phantom', 'circ_oscil',
%       'fnoise' (1/f^alpha noise).
%   Natural images:
%       'boat', 'lena', 'goldhill', 'mandrill', 'maurice', 'polygons_blurred', or your own.
%
%   Copyright (c) 2004 Gabriel Peyre

if nargin<2
    n = 512;
end
options.null = 0;

type = lower(type);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% parameters for geometric objects
eta = 0.1;              % translation
gamma = 1/sqrt(2);      % slope
if isfield( options, 'eta' )
    eta = options.eta;
end
if isfield( options, 'gamma' )
    eta = options.gamma;
end
if isfield( options, 'radius' )
    radius = options.radius;
end
if isfield( options, 'center' )
    center = options.center;
end
if isfield( options, 'center1' )
    center1 = options.center1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% for the line, can be vertical / horizontal / diagonal / any
if strcmp(type, 'line_vertical')
    eta = 0.5;              % translation
    gamma = 0;      % slope
elseif strcmp(type, 'line_horizontal')
    eta = 0.5;              % translation
    gamma = Inf;      % slope
elseif strcmp(type, 'line_diagonal')
    eta = 0;              % translation
    gamma = 1;      % slope
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% for some blurring
sigma = 0;
if isfield(options, 'sigma')
    sigma = options.sigma;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch lower(type)
    case {'letter-x' 'letter-v' 'letter-z' 'letter-y'}
        if isfield(options, 'radius')
            r = options.radius;
        else
            r = 10;
        end
        M = create_letter(type(8), r, n);

    case 'l'
        r1 = [.1 .1  .3 .9];
        r2 = [.1 .1 .9 .4];
        M = double( draw_rectangle(r1,n) | draw_rectangle(r2,n) );

    case 'ellipse'
        c1 = [0.15 0.5];
        c2 = [0.85 0.5];
        if isfield(options, 'eccentricity')
            eccentricity = options.eccentricity;
        else
            eccentricity = 1.3;
        end
        x = linspace(0,1,n);
        [Y,X] = meshgrid(x,x);
        d = sqrt((X-c1(1)).^2 + (Y-c1(2)).^2) + sqrt((X-c2(1)).^2 + (Y-c2(2)).^2);
        M = double( d<=eccentricity*sqrt( sum((c1-c2).^2) ) );
    case 'ellipse-thin'
        options.eccentricity = 1.06;
        M = load_image('ellipse', n, options);
    case 'ellipse-fat'
        options.eccentricity = 1.3;
        M = load_image('ellipse', n, options);

    case 'square-tube'
        if isfield(options, 'tube_width')
            w = options.tube_width;
        else
            w = 0.06;
        end
        c1 = [.25 .5];
        c2 = [.75 .5];
        r1 = [c1 c1] + .18*[-1 -1 1 1];
        r2 = [c2 c2] + .18*[-1 -1 1 1];
        r3 = [c1(1)-w c1(2)-w c2(1)+w c2(2)+w];
        M = double( draw_rectangle(r1,n) | draw_rectangle(r2,n) | draw_rectangle(r3,n) );
    case 'square-tube-1'
        options.tube_width = 0.03;
        M = load_image('square-tube', n, options);
    case 'square-tube-2'
        options.tube_width = 0.06;
        M = load_image('square-tube', n, options);
    case 'square-tube-3'
        options.tube_width = 0.09;
        M = load_image('square-tube', n, options);
    case 'polygon'
        if isfield(options, 'nb_points')
            nb_points = options.nb_points;
        else
            nb_points = 9;
        end
        if isfield(options, 'scaling')
            scaling = options.scaling;
        else
            scaling = 1;
        end
        theta = sort( rand(nb_points,1)*2*pi );
        radius = scaling*rescale(rand(nb_points,1), 0.1, 0.93);
        points = [cos(theta) sin(theta)] .* repmat(radius, 1,2);
        points = (points+1)/2*(n-1)+1; points(end+1,:) = points(1,:);
        M = draw_polygons(zeros(n),0.8,{points'});
        [x,y] = ind2sub(size(M),find(M));
        p = 100; m = length(x);
        lambda = linspace(0,1,p);
        X = n/2 + repmat(x-n/2, [1 p]) .* repmat(lambda, [m 1]);
        Y = n/2 + repmat(y-n/2, [1 p]) .* repmat(lambda, [m 1]);
        I = round(X) + (round(Y)-1)*n;
        M = zeros(n); M(I) = 1;
    case 'polygon-8'
        options.nb_points = 8;
        M = load_image('polygon', n, options);
    case 'polygon-10'
        options.nb_points = 8;
        M = load_image('polygon', n, options);
    case 'polygon-12'
        options.nb_points = 8;
        M = load_image('polygon', n, options);
    case 'pacman'
        if isfield(options, 'theta')
            theta = options.theta;
        else
            theta = 30 * 2*pi/360;
        end
        options.radius = 0.45;
        M = load_image('disk', n, options);
        x = linspace(-1,1,n);
        [Y,X] = meshgrid(x,x);
        T =atan2(Y,X);
        M = M .* (1-(abs(T)<theta/2));
    case 'square-hole'
        options.radius = 0.45;
        M = load_image('disk', n, options);
        options.scaling = 0.5;
        M = M - load_image('polygon-10', n, options);

    case 'grid-circles'
        if isempty(n)
            n = 256;
        end
        if isfield(options, 'frequency')
            f = options.frequency;
        else
            f = 30;
        end
        if isfield(options, 'width')
            eta = options.width;
        else
            eta = 0.3;
        end
        x = linspace(-n/2,n/2,n) - round(n*0.03);
        y = linspace(0,n,n);
        [Y,X] = meshgrid(y,x);
        R = sqrt(X.^2+Y.^2);
        theta = 0.05*pi/2;
        X1 = cos(theta)*X+sin(theta)*Y;
        Y1 = -sin(theta)*X+cos(theta)*Y;
        A1 = abs(cos(2*pi*R/f))<eta;
        A2 = max( abs(cos(2*pi*X1/f))<eta, abs(cos(2*pi*Y1/f))<eta );

        M = A1;
        M(X1>0) = A2(X1>0);

    case 'chessboard1'
        x = -1:2/(n-1):1;
        [Y,X] = meshgrid(x,x);
        M = (2*(Y>=0)-1).*(2*(X>=0)-1);

    case 'chessboard'
        if ~isfield( options, 'width' )
            width = round(n/16);
        else
            width = options.width;
        end
        [Y,X] = meshgrid(0:n-1,0:n-1);
        M = mod( floor(X/width)+floor(Y/width), 2 ) == 0;

    case 'square'
        if ~isfield( options, 'radius' )
            radius = 0.6;
        end
        x = -1:2/(n-1):1;
        [Y,X] = meshgrid(x,x);
        M = max( abs(X),abs(Y) )<radius;

    case 'squareregular'
        M = rescale(load_image('square',n,options));
        if not(isfield(options, 'alpha'))
            options.alpha = 3;
        end
        S = load_image('fnoise',n,options);
        M = M + rescale(S,-0.3,0.3);

    case 'regular1'
        options.alpha = 1;
        M = load_image('fnoise',n,options);
    case 'regular2'
        options.alpha = 2;
        M = load_image('fnoise',n,options);
    case 'regular3'
        options.alpha = 3;
        M = load_image('fnoise',n,options);

    case 'sparsecurves'
        options.alpha = 3;
        M = load_image('fnoise',n,options);
        M = rescale(M);
        ncurves = 3;
        M = cos(2*pi*ncurves);

    case 'square_texture'
        M = load_image('square',n);
        M = rescale(M);
        % make a texture patch
        x = linspace(0,1,n);
        [Y,X] = meshgrid(x,x);
        theta = pi/3;
        x = cos(theta)*X + sin(theta)*Y;
        c = [0.3,0.4]; r = 0.2;
        I = find( (X-c(1)).^2 + (Y-c(2)).^2 < r^2 );
        eta = 3/n; lambda = 0.3;
        M(I) = M(I) + lambda * sin( x(I) * 2*pi / eta );

    case 'oscillatory_texture'
        x = linspace(0,1,n);
        [Y,X] = meshgrid(x,x);
        theta = pi/3;
        x = cos(theta)*X + sin(theta)*Y;
        c = [0.3,0.4]; r = 0.2;
        I = find( (X-c(1)).^2 + (Y-c(2)).^2 < r^2 );
        eta = 3/n; lambda = 0.3;
        M = sin( x * 2*pi / eta );

    case {'line', 'line_vertical', 'line_horizontal', 'line_diagonal'}
        x = 0:1/(n-1):1;
        [Y,X] = meshgrid(x,x);
        if gamma~=Inf
            M = (X-eta) - gamma*Y < 0;
        else
            M = (Y-eta) < 0;
        end

    case 'grating'
        x = linspace(-1,1,n);
        [Y,X] = meshgrid(x,x);
        if isfield(options, 'theta')
            theta = options.theta;
        else
            theta = 0.2;
        end
        if isfield(options, 'freq')
            freq = options.freq;
        else
            freq = 0.2;
        end
        X = cos(theta)*X + sin(theta)*Y;
        M = sin(2*pi*X/freq);

    case 'disk'
        if ~isfield( options, 'radius' )
            radius = 0.35;
        end
        if ~isfield( options, 'center' )
            center = [0.5, 0.5];    % center of the circle
        end
        x = 0:1/(n-1):1;
        [Y,X] = meshgrid(x,x);
        M = (X-center(1)).^2 + (Y-center(2)).^2 < radius^2;

    case 'diskregular'
        M = rescale(load_image('disk',n,options));
        if not(isfield(options, 'alpha'))
            options.alpha = 3;
        end
        S = load_image('fnoise',n,options);
        M = M + rescale(S,-0.3,0.3);

    case 'quarterdisk'
        if ~isfield( options, 'radius' )
            radius = 0.95;
        end
        if ~isfield( options, 'center' )
            center = -[0.1, 0.1];    % center of the circle
        end
        x = 0:1/(n-1):1;
        [Y,X] = meshgrid(x,x);
        M = (X-center(1)).^2 + (Y-center(2)).^2 < radius^2;

    case 'fading_contour'
        if ~isfield( options, 'radius' )
            radius = 0.95;
        end
        if ~isfield( options, 'center' )
            center = -[0.1, 0.1];    % center of the circle
        end
        x = 0:1/(n-1):1;
        [Y,X] = meshgrid(x,x);
        M = (X-center(1)).^2 + (Y-center(2)).^2 < radius^2;
        theta = 2/pi*atan2(Y,X);
        h = 0.5;
        M = exp(-(1-theta).^2/h^2).*M;

    case '3contours'
        radius = 1.3;
        center = [-1, 1];
        radius1 = 0.8;
        center1 = [0, 0];
        x = 0:1/(n-1):1;
        [Y,X] = meshgrid(x,x);
        f1 = (X-center(1)).^2 + (Y-center(2)).^2 < radius^2;
        f2 = (X-center1(1)).^2 + (Y-center1(2)).^2 < radius1^2;
        M = f1 + 0.5*f2.*(1-f1);

    case 'line_circle'

        gamma = 1/sqrt(2);

        x = linspace(-1,1,n);
        [Y,X] = meshgrid(x,x);
        M1 = double( X>gamma*Y+0.25 );
        M2 = X.^2 + Y.^2 < 0.6^2;
        M = 20 + max(0.5*M1,M2) * 216;

    case 'fnoise'

        % generate an image M whose Fourier spectrum amplitude is
        %   |M^(omega)| = 1/f^{omega}
        if isfield(options, 'alpha')
            alpha = options.alpha;
        else
            alpha = 1;
        end
        M = gen_noisy_image(n,alpha);

    case 'gaussiannoise'
        % generate an image of filtered noise with gaussian
        if isfield(options, 'sigma')
            sigma = options.sigma;
        else
            sigma = 10;
        end
        M = randn(n);
        m = 51;
        h = compute_gaussian_filter([m m],sigma/(4*n),[n n]);
        M = perform_convolution(M,h);
        return;

    case {'bwhorizontal','bwvertical','bwcircle'}

        [Y,X] = meshgrid(0:n-1,0:n-1);
        if strcmp(type, 'bwhorizontal')
            d = X;
        elseif strcmp(type, 'bwvertical')
            d = Y;
        elseif strcmp(type, 'bwcircle')
            d = sqrt( (X-(n-1)/2).^2 + (Y-(n-1)/2).^2 );
        end
        if isfield(options, 'stripe_width')
            stripe_width = options.stripe_width;
        else
            stripe_width = 5;
        end
        if isfield(options, 'black_prop')
            black_prop = options.black_prop;
        else
            black_prop = 0.5;
        end
        M = double( mod( d/(2*stripe_width),1 )>=black_prop );

    case 'parabola'

        % curvature
        if isfield(options, 'c')
            c = options.c;
        else
            c = 0.1;
        end
        % angle
        if isfield(options, 'theta');
            theta = options.theta;
        else
            theta = pi/sqrt(2);
        end
        x = -0.5:1/(n-1):0.5;
        [Y,X] = meshgrid(x,x);
        Xs = X*cos(theta) + Y*sin(theta);
        Y =-X*sin(theta) + Y*cos(theta); X = Xs;
        M = Y>c*X.^2;

    case 'sin'

        [Y,X] = meshgrid(-1:2/(n-1):1, -1:2/(n-1):1);
        M = Y >= 0.6*cos(pi*X);
        M = double(M);

    case 'circ_oscil'

        x = linspace(-1,1,n);
        [Y,X] = meshgrid(x,x);
        R = sqrt(X.^2+Y.^2);
        M = cos(R.^3*50);

    case 'phantom'

        M = phantom(n);

    case 'periodic_bumps'

        if isfield(options, 'nbr_periods')
            nbr_periods = options.nbr_periods;
        else
            nbr_periods = 8;
        end
        if isfield(options, 'theta')
            theta = options.theta;
        else
            theta = 1/sqrt(2);
        end
        if isfield(options, 'skew')
            skew = options.skew;
        else
            skew = 1/sqrt(2);
        end

        A = [cos(theta), -sin(theta); sin(theta), cos(theta)];
        B = [1 skew; 0 1];
        T = B*A;
        x = (0:n-1)*2*pi*nbr_periods/(n-1);
        [Y,X] = meshgrid(x,x);
        pos = [X(:)'; Y(:)'];
        pos = T*pos;
        X = reshape(pos(1,:), n,n);
        Y = reshape(pos(2,:), n,n);
        M = cos(X).*sin(Y);

    case 'noise'
        if isfield(options, 'sigma')
            sigma = options.sigma;
        else
            sigma = 1;
        end
        M = randn(n);

    otherwise
        ext = {'gif', 'png', 'jpg', 'bmp', 'tiff', 'pgm', 'ppm'};
        for i=1:length(ext)
            name = [type '.' ext{i}];
            if( exist(name) )
                M = imread( name );
                M = double(M);
                if not(isempty(n)) && (n~=size(M, 1) || n~=size(M, 2)) && nargin>=2
                    M = image_resize(M,n,n);
                end
                return;
            end
        end
        error( ['Image ' type ' does not exists.'] );
end

M = double(M);

if sigma>0
    h = compute_gaussian_filter( [9 9], sigma/(2*n), [n n]);
    M = perform_convolution(M,h);
end

M = rescale(M) * 256;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function M = create_letter(a, r, n)

c = 0.2;
p1 = [c;c];
p2 = [c; 1-c];
p3 = [1-c; 1-c];
p4 = [1-c; c];
p4 = [1-c; c];
pc = [0.5;0.5];
pu = [0.5; c];

switch a
    case 'x'
        point_list = { [p1 p3] [p2 p4] };
    case 'z'
        point_list = { [p2 p3 p1 p4] };
    case 'v'
        point_list = { [p2 pu p3] };
    case 'y'
        point_list = { [p2 pc pu] [pc p3] };

end
% fit image
for i=1:length(point_list)
    a = point_list{i}(2:-1:1,:);
    a(1,:) = 1-a(1,:);
    point_list{i} = round( a*(n-1)+1 );
end
M = draw_polygons(zeros(n),r,point_list);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function sk = draw_polygons(mask,r,point_list)

sk = mask*0;
for i=1:length(point_list)
    pl = point_list{i};
    for k=2:length(pl)
        sk = draw_line(sk,pl(1,k-1),pl(2,k-1),pl(1,k),pl(2,k),r);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function sk = draw_line(sk,x1,y1,x2,y2,r)

n = size(sk,1);
[Y,X] = meshgrid(1:n,1:n);
q = 100;
t = linspace(0,1,q);
x = x1*t+x2*(1-t); y = y1*t+y2*(1-t);
if r==0
    x = round( x ); y = round( y );
    sk( x+(y-1)*n ) = 1;
else
    for k=1:q
        I = find((X-x(k)).^2 + (Y-y(k)).^2 <= r^2 );
        sk(I) = 1;
    end
end

function M = gen_noisy_image(n,alpha)

% gen_noisy_image - generate a noisy cloud-like image.
%
%   M = gen_noisy_image(n,alpha);
%
% generate an image M whose Fourier spectrum amplitude is
%   |M^(omega)| = 1/f^{omega}
%
%   Copyright (c) 2004 Gabriel Peyr?

if nargin<1
    n = 128;
end
if nargin<2
    alpha = 1.5;
end

if mod(n(1),2)==0
    x = -n/2:n/2-1;
else
    x = -(n-1)/2:(n-1)/2;
end

[Y,X] = meshgrid(x,x);
d = sqrt(X.^2 + Y.^2) + 0.1;
f = rand(n)*2*pi;

M = (d.^(-alpha)) .* exp(f*1i);
% M = real(ifft2(fftshift(M)));

M = ifftshift(M);
M = real( ifft2(M) );

function y = gen_signal_2d(n,alpha)

% gen_signal_2d -  generate a 2D C^\alpha signal of length n x n.
%   gen_signal_2d(n,alpha) generate a 2D signal C^alpha.
%
%   The signal is scale in [0,1].
%
%   Copyright (c) 2003 Gabriel Peyr?

% new new method

[Y,X] = meshgrid(0:n-1, 0:n-1);

A = X+Y+1;
B = X-Y+n+1;

a = gen_signal(2*n+1, alpha);
b = gen_signal(2*n+1, alpha);
y = a(A).*b(B);
% M = a(1:n)*b(1:n)';

return;

% new method
h = (-n/2+1):(n/2); h(n/2)=1;
[X,Y] = meshgrid(h,h);
h = sqrt(X.^2+Y.^2+1).^(-alpha-1/2);
h = h .* exp( 2i*pi*rand(n,n) );
h = fftshift(h);
y = real( ifft2(h) );

m1 = min(min(y));
m2 = max(max(y));
y = (y-m1)/(m2-m1);

return;

%% old code

y = rand(n,n);
y = y - mean(mean(y));
for i=1:alpha
    y = cumsum(cumsum(y)')';
    y = y - mean(mean(y));
end
m1 = min(min(y));
m2 = max(max(y));
y = (y-m1)/(m2-m1);

function newimg = image_resize(img,p1,q1,r1)

% image_resize - resize an image using bicubic interpolation
%
%   newimg = image_resize(img,nx,ny,nz);
% or
%   newimg = image_resize(img,newsize);
%
%   Works for 2D, 2D 2 or 3 channels, 3D images.
%
%   Copyright (c) 2004 Gabriel Peyr?

if nargin==2
    % size specified as an array
    q1 = p1(2);
    if length(p1)>2
        r1 = p1(3);
    else
        r1 = size(img,3);
    end
    p1 = p1(1);
end

if nargin<4
    r1 = size(img,3);
end

if ndims(img)<2 || ndims(img)>3
    error('Works only for grayscale or color images');
end

if ndims(img)==3 && size(img,3)<4
    % RVB image
    newimg = zeros(p1,q1, size(img,3));
    for m=1:size(img,3)
        newimg(:,:,m) = image_resize(img(:,:,m), p1, q1);
    end
    return;
elseif ndims(img)==3
    p = size(img,1);
    q = size(img,2);
    r = size(img,3);
    [Y,X,Z] = meshgrid( (0:q-1)/(q-1), (0:p-1)/(p-1), (0:r-1)/(r-1)  );
    [YI,XI,ZI] = meshgrid( (0:q1-1)/(q1-1), (0:p1-1)/(p1-1), (0:r1-1)/(r1-1) );
    newimg = interp3( Y,X,Z, img, YI,XI,ZI ,'cubic');
    return;
end

p = size(img,1);
q = size(img,2);
[Y,X] = meshgrid( (0:q-1)/(q-1), (0:p-1)/(p-1) );
[YI,XI] = meshgrid( (0:q1-1)/(q1-1), (0:p1-1)/(p1-1) );
newimg = interp2( Y,X, img, YI,XI ,'cubic');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function M = draw_rectangle(r,n)

x = linspace(0,1,n);
[Y,X] = meshgrid(x,x);
M = double( (X>=r(1)) & (X<=r(3)) & (Y>=r(2)) & (Y<=r(4)) ) ;


text:
function ft_realtime_fmriviewer(cfg)

% FT_REALTIME_FMRIVIEWER allows for realtime visualization of the fMRI data stream
%
% Use as
%   ft_realtime_fmriviewer(cfg)
% with the following configuration options
%   cfg.bufferdata = whether to start on the 'first or 'last' data that is available (default = 'last')
%
% The source of the data is configured as
%   cfg.dataset       = string
% or alternatively to obtain more low-level control as
%   cfg.datafile      = string
%   cfg.headerfile    = string
%   cfg.eventfile     = string
%   cfg.dataformat    = string, default is determined automatic
%   cfg.headerformat  = string, default is determined automatic
%   cfg.eventformat   = string, default is determined automatic

% Copyright (C) 2012, Stefan Klanke
%
% 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$

if ~isfield(cfg, 'headerfile'),    	cfg.headerfile = 'buffer://localhost:1972'; end
if ~isfield(cfg, 'datafile'),      	cfg.datafile = cfg.headerfile; end
if ~isfield(cfg, 'eventfile'),    	cfg.eventfile = cfg.headerfile; end
if ~isfield(cfg, 'headerformat'),   cfg.headerformat = []; end
if ~isfield(cfg, 'dataformat'),   	cfg.dataformat = []; end
if ~isfield(cfg, 'bufferdata'),     cfg.bufferdata = 'last';  end % first or last

% translate dataset into datafile+headerfile
cfg = ft_checkconfig(cfg, 'dataset2files', 'yes');
cfg = ft_checkconfig(cfg, 'required', {'datafile' 'headerfile'});

% ensure that the persistent variables related to caching are cleared
clear ft_read_header
% start by reading the header from the realtime buffer
hdr = ft_read_header(cfg.headerfile);

disp(hdr);

if isfield(hdr,'nifti_1')
	ddim = double(hdr.nifti_1.dim);
    width  = ddim(1);
    height = ddim(2);
    numSlices = ddim(3);
elseif isfield(hdr,'siemensap') && isstruct(hdr.siemensap)
    width = siemensap.sKSpace.lBaseResolution;
    phaseFOV = siemensap.sSliceArray.asSlice{1}.dPhaseFOV;
    readoutFOV = siemensap.sSliceArray.asSlice{1}.dReadoutFOV;
    height = width * phaseFOV / readoutFOV;
    numSlices = siemensap.sSliceArray.lSize;
else
  ft_warning('No protocol information found!')
  width = sqrt(hdr.nChans);
  height = width;
  numSlices = 1;
end

% open the figure and set the colormap
h = figure;
colormap(gray);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is the general BCI loop where realtime incoming data is handled
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

blocksize  = 1;
prevSample = 0;

% simple smoothing kernel
kern1d = [1 2 1];
kern = convn(kern1d'*kern1d,reshape(kern1d,1,1,3))
kern = kern/sum(kern(:));

while true

  % determine number of samples available in buffer
  hdr = ft_read_header(cfg.headerfile, 'headerformat', cfg.headerformat, 'cache', true);

  if hdr.nSamples > prevSample
    % determine the samples to process
    if strcmp(cfg.bufferdata, 'last')
      begsample  = hdr.nSamples*hdr.nTrials - blocksize + 1;
      endsample  = hdr.nSamples*hdr.nTrials;
    elseif strcmp(cfg.bufferdata, 'first')
      begsample  = prevSample+1;
      endsample  = prevSample+blocksize;
    else
      ft_error('unsupported value for cfg.bufferdata');
    end

    prevSample  = endsample;
    fprintf('\nprocessing scan %d\n', begsample);

    % read data from buffer (only the last scan)
    dat = ft_read_data(cfg.datafile, 'header', hdr, 'dataformat', cfg.dataformat, 'begsample', begsample, 'endsample', endsample);

    % read events from buffer
    ev = ft_read_event(cfg.eventfile, 'header', hdr);
    if ~isempty(ev)
      ev_sam = [ev.sample];
      ind = find(ev_sam==begsample-1);
      if ~isempty(ind)
        disp([ev(ind).type ev(ind).value])
      end
    end

    % go to the correct figure
    figure(h);

    % actually, you should reshape to 3D (remove the *)
    S = reshape(dat, [width  height*numSlices]);

    V = reshape(dat, [width  height numSlices]);
    tic;
    Vn = convn(V,kern);
    toc
    Vn = Vn(2:end-1,2:end-1,2:end-1);
    size(Vn)
    S = reshape(Vn, [width height*numSlices]);

    imagesc(S(:,:)',[0 2048]);

    % force Matlab to update the figure
    drawnow
  else
    pause(0.01);
  end % if enough new samples
end % while true

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ solutionsValid ] = validateOptForceSol(model,posOptForceSets,typeRegOptForceSets, target)
relreacs = unique(model.rxns(posOptForceSets(:)));
[minRelFluxes,maxRelFluxes] = fluxVariability(model,0,'max',relreacs);

solutionsValid = true;
sol = optimizeCbModel(model);
tempmodel = model;
tempmodel.lb(find(model.c)) = sol.x(find(model.c));
tempmodel.c = double(ismember(model.rxns,target));
tempmodel.osenseStr = 'min';
sol = optimizeCbModel(tempmodel);
ores = sol.x(ismember(model.rxns,target));

%These changes are extreme, and they should only be used for the actual
%test Case. In general they might not be applicable.
changes = {'knockout',@(cmodel,treac) changeRxnBounds(cmodel,treac,0,'b');...
           'upregulation',@(cmodel,treac) changeRxnBounds(cmodel,treac,min(sol.x(ismember(cmodel.rxns,treac))+100,maxRelFluxes(ismember(relreacs,treac))),'l');...
           'downregulation',@(cmodel,treac) changeRxnBounds(cmodel,treac,max(sol.x(ismember(cmodel.rxns,treac))-100,minRelFluxes(ismember(relreacs,treac))),'u')};
for i = 1:size(posOptForceSets,1)
    cmodel = model;
    for j = 1:size(posOptForceSets,2)
        if posOptForceSets(i,j) > 0
            cfunc = changes{ismember(changes(:,1),typeRegOptForceSets{i,j}),2};
            cmodel = cfunc(cmodel,cmodel.rxns(posOptForceSets(i,j)));
        end
    end
    csol = optimizeCbModel(cmodel);
    cmodel.lb(find(cmodel.c)) = csol.x(find(cmodel.c));
    cmodel.c = double(ismember(model.rxns,target));
    cmodel.osenseStr = 'min';
    csol = optimizeCbModel(cmodel);
    if ~(csol.x(ismember(model.rxns,target)) > ores)
        solutionsValid = false;
        return
    end
end

author: opencobra

repo: cobratoolbox


text:
function bf10 = corrbf(r,n)
%
% bf10 = corrbf(r,n)
%
% Calculates JZS Bayes factor for correlation r and sample size n.
% This quantifies the evidence in favour of the alternative hypothesis.
% See Wetzels & Wagemakers, 2012, Psychon Bull Rev for details.
%

% Function to be integrated - updated after personal communication from EJW (old one didn't work with large N)
F = @(g,r,n) exp(((n-2)./2).*log(1+g)+(-(n-1)./2).*log(1+(1-r.^2).*g)+(-3./2).*log(g)+-n./(2.*g));

% Bayes factor calculation
bf10 = sqrt((n/2)) / gamma(1/2) * integral(@(g) F(g,r,n),0,Inf);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function close_gripper(coppelia, robot_number)
    robot_name = strcat('R', int2str(robot_number));
    gripper_name = strcat(robot_name, '_gripper');

    if gripper_name =='R2_gripper';
        gripper_name='R1_gripper#0';
    end

    [res retInts retFloats retStrings retBuffer]=coppelia.sim.simxCallScriptFunction(coppelia.clientID,gripper_name,coppelia.sim.sim_scripttype_childscript,'open_close_gripper_function',[0],[],[],[],coppelia.sim.simx_opmode_blocking);
    coppelia_wait(coppelia, 10)
end

repo: ARTE


text:
function f = p06_f ( x )

%*****************************************************************************80
%
%% P06_F evaluates the objective function for problem 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Richard Brent,
%    Algorithms for Minimization Without Derivatives,
%    Prentice Hall 1973,
%    Reprinted Dover, 2002
%
%  Parameters:
%
%    Input, real X, the argument of the objective function.
%
%    Output, real F, the value of the objective function.
%
  f = 2.0 - x;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function T = permute(T,order)
%PERMUTE Permute tensor dimensions.
%
%   B = PERMUTE(A,ORDER) rearranges the dimensions of A so that they
%   are in the order specified by the vector ORDER. The result has the
%   same values of A, but the order of the subscripts needed to access
%   any particular element are rearranged as specified by ORDER.
%
%   See also TENSOR, TENSOR/SIZE, TENSOR/NDIMS, PERMUTE.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if ndims(T) ~= numel(order)
  error('Invalid permutation order');
end

% Check for special case of permuting an order-1 object (which has
% no effect but confuses MATLAB's permute command which doesn't
% think that there is such a thing as a 1D-array).
if isequal(order,1)
    return;
end

% Check for special case of empty object (which has
% no effect but confuses MATLAB's permute command which doesn't
% think that there is such a thing as an empty array).
if isempty(order)
    return;
end

% Note that permute does error checking on order, so we don't worry
% about it.
T.data = permute(T.data,order);
T.size = T.size(order);

return;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


function [ err ] = residual_KR_robust( X1, X2, imsize1, imsize2, paras, sigma )

% parameretes sigma indicates the distance scope for inliers with homopraphy matrix, in pixels

k1 = paras(1);
k2 = paras(2);
theta = paras(3:5);
% yaw = paras(3);
% pitch = paras(4);
% roll = paras(5);

K1 = [k1, 0, imsize1(2)/2;
     0, k1, imsize1(1)/2;
     0,  0, 1];
K2 = [k2, 0, imsize2(2)/2;
      0, k2, imsize2(1)/2;
      0,  0, 1];
theta_m = [0         -theta(3) theta(2)
           theta(3)  0         -theta(1)
           -theta(2) theta(1)  0];
R = expm(theta_m);
% Ry = [cos(yaw),     0,              -sin(yaw);
%       0,            1,              0;
%       sin(yaw),     0,              cos(yaw)] ;
% Rp = [1,            0               0;
%       0,            cos(pitch),     sin(pitch);
%       0,            -sin(pitch),    cos(pitch)] ;
% Rr = [cos(roll),    sin(roll),      0;
%       -sin(roll),   cos(roll),      0;
%       0,            0,              1] ;
% R = Rr * Rp * Ry;

err = residual_H(X1, X2, K1, K2, R);

outlier = (abs(err) > sigma);
err(outlier) = sign(err(outlier)) .* (sigma + sigma * log(abs(err(outlier))/sigma));
% err(outlier) = sign(err(outlier)) .* sqrt(2*sigma*abs(err(outlier)) - sigma*sigma);

end


text:
function N = spm_dicom_metadata(N,hdr)
% Export image metadata as side-car JSON file
% FORMAT N = spm_dicom_metadata(N,hdr)
% N(input)  - nifti object
% hdr       - a single header from spm_dicom_headers
% N(output) - unchanged nifti object (for potential future use)
%
% This function creates JSON-encoded metadata during DICOM to NIfTI
% conversion, including all acquisition parameters, and saves them as a
% JSON side-car file.
%
% See also: spm_dicom_convert
%__________________________________________________________________________
% Copyright (C) 2017-2018 Wellcome Trust Centre for Neuroimaging

% Evelyne Balteau, Cyclotron Research Centre, University of Liege, Belgium
% $Id: spm_dicom_metadata.m 7482 2018-11-12 12:18:08Z guillaume $

%-Provenance and general description of the image (hMRI toolbox)
dicom_convert_version = {sprintf('%s %s', spm_check_version, version),sprintf('spm_dicom_convert.m - %s', spm('Version'))};
metadata.history.procstep = struct('descrip','dicom to nifti import', 'version', {dicom_convert_version}, 'procpar', []);
metadata.history.input(1) = struct('filename','AnonymousFileName', 'history',[]);
if isfield(hdr,'ImageType')
    metadata.history.output = struct('imtype',hdr.ImageType, 'units','a.u.');
else
    metadata.history.output = struct('imtype','Unprocessed MR image', 'units','a.u.');
end

%-Acquisition parameters with complete DICOM header (hMRI toolbox)
hdr = reformat_spm_dicom_header(hdr);
hdr = anonymise_metadata(hdr); % default is basic anonymisation
metadata.acqpar = hdr;

%-Save metadata as side-car JSON file
spm_jsonwrite(spm_file(N.dat.fname,'ext','json'), metadata, struct('indent','\t'));

%==========================================================================
function hdr = reformat_spm_dicom_header(hdr)
% To tidy up and rearrange CSA fields in the header, including formatting
% the ASCII part into a proper MATLAB structure (Note: this is specific to
% Siemens DICOM format but could be extended to other cases where needed)
%
% FORMAT hdr = reformat_spm_dicom_header(hdr)
% where hdr is a dicom header structure as output by spm_dicom_headers.

% list of CSA-like fields:
CSAlist = {'CSAImageHeaderInfo', 'CSASeriesHeaderInfo','CSANonImageHeaderInfoVA','CSAMiscProtocolHeaderInfoVA','CSANonImageHeaderInfoVB','CSAMiscProtocolHeaderInfoVB'};
% NB: might be necessary to add cases 'Private_0029_1110' and
% 'Private_0029_1210' for spectroscopic data (see spm_dicom_essentials.m)

for ccsa = 1:length(CSAlist)
    if isfield(hdr,CSAlist{ccsa})
        tdyhdr = tidy_CSA(hdr.(CSAlist{ccsa}));
        if isfield(tdyhdr,'MrPhoenixProtocol')
            % works for Siemens VB, VD & VE DICOM format:
            tdyhdr.MrPhoenixProtocol = read_ASCII(tdyhdr.MrPhoenixProtocol);
        elseif isfield(tdyhdr,'MrProtocol')
            % works for Siemens VA DICOM format:
            tdyhdr.MrProtocol = read_ASCII(tdyhdr.MrProtocol);
        end
        hdr.(CSAlist{ccsa}) = tdyhdr;
    end
end

%==========================================================================
function hdrout = anonymise_metadata(hdr, opts)
% FORMAT hdrout = anonymise_metadata(hdr, opts)
% hdr is a DICOM header read with spm_dicom_header
% opts are anonymisation options:
%       opts.anonym = 'none': no anonymisation, all patient data kept in
%                       the metadata.
%                     'full': no patient information is kept at all
%                     'basic': patient ID (presumably not containing his
%                       name), age (years at the time of the data
%                       acquisition), sex, size and weight are kept,
%                       patient name, date of birth and DICOM filename
%                       (often containing the patient name) are removed.
%
% !!! IMPORTANT WARNING: EFFECTIVE ANONYMISATION IS NOT GUARANTEED !!!
% The anonymisation implemented here depends on the structure and content
% of the DICOM header and might not be effective in many cases.

if nargin<2
    opts.anonym = 'basic';
end
if ~iscell(hdr)
    hdrout{1} = hdr;
else
    hdrout = hdr;
end

if ~strcmp(opts.anonym,'none')
    for i=1:length(hdrout)
        if isfield(hdrout{i},'PatientName')
            hdrout{i}.PatientName = 'Anonymous';
        end
        if isfield(hdrout{i},'PatientBirthDate')
            t1 = hdrout{i}.PatientBirthDate;
            t2 = hdrout{i}.StudyDate;
            hdrout{i}.PatientAge = round((t2-t1)*10/365.25)/10;
            hdrout{i} = rmfield(hdrout{i},'PatientBirthDate');
        end
        if isfield(hdrout{i},'Filename')
            hdrout{i} = rmfield(hdrout{i},'Filename');
        end

        if strcmp(opts.anonym, 'full')
            try, hdrout{i} = rmfield(hdrout{i},'PatientID'); end
            try, hdrout{i} = rmfield(hdrout{i},'PatientSex'); end
            try, hdrout{i} = rmfield(hdrout{i},'PatientAge'); end
            try, hdrout{i} = rmfield(hdrout{i},'PatientSize'); end
            try, hdrout{i} = rmfield(hdrout{i},'PatientWeight'); end
            try, hdrout{i}.CSASeriesHeaderInfo = rmfield(hdrout{i}.CSASeriesHeaderInfo,'UsedPatientWeight'); end
        end
    end
end

%==========================================================================
function tdyhdr = tidy_CSA(csahdr)
% FORMAT tdyhdr = tidy_CSA(csahdr)
% DESCRIPTION:
% To rearrange CSAImageHeaderInfo, CSASeriesHeaderInfo, ... structures
% and make it a simplified, easier to browse structure (Siemens specific).
% USAGE:
% tdyhdr = tidy_CSA(csahdr)
% where csahdr is a structure, content of the CSA field.

tdyhdr = [];
for ccsa = 1:length(csahdr)
    val = get_numaris4_val(csahdr,csahdr(ccsa).name);
    % if val is empty, let's not waste disk space with empty fields
    if ~isempty(val)
        % for elements having a value representation corresponding to a
        % numerical value (or an array of numerical values), convert char
        % into numbers. Here is a list of VR corresponding to numerical
        % values:
        % - DS (Decimal String)
        % - FL (Floating Point Single)
        % - FD (Floating Point Double)
        % - IS (Integer String)
        % - OD (Other Double String)
        % - OF (Other Float String)
        % - SL (Signed Long)
        % - SS (Signed Short)
        % - UL (Unsigned Long)
        % - US (Unsigned Short)
        switch deblank(csahdr(ccsa).vr)
            case {'DS','FL','FD','IS','OD','OF','SL','SS','UL','US'}
                try
                    tmp = zeros(size(val,1),1);
                    for k = 1:size(val,1)
                        tmp(k) = str2num(val(k,:));
                    end
                    val = tmp;
                catch
                    fprintf('Trouble reading CSA header %s (%s) = %s\n',csahdr(ccsa).name, deblank(csahdr(ccsa).vr), val(1,:));
                    val = [];
                end
            otherwise
                val = deblank(val);
        end
        % make sure no "-" in field name (invalid otherwise)
        csahdr(ccsa).name = strrep(csahdr(ccsa).name,'-','');
        tdyhdr.(csahdr(ccsa).name) = val;
    end
end

%==========================================================================
function val = get_numaris4_val(str,name)
% copied from spm_dicom_convert
name = deblank(name);
val  = {};
for i=1:length(str)
    if strcmp(deblank(str(i).name),name)
        if isfield(str(i),'nitems')
            for j=1:str(i).nitems
                if  str(i).item(j).xx(1)
                    val = [val {str(i).item(j).val}];
                end
            end
            %         else
            %             fprintf(1,'Found empty CSA header with name %s\n',name);
        end
        break;
    end
end
val = strvcat(val{:});

%==========================================================================
function ascout = read_ASCII(ascin)
% FORMAT ascout = read_ASCII(ascin)
%
% DESCRIPTION:
% To parse the content of the ASCII part of the DICOM header (Siemens
% specific) and convert it into a matlab structure. Compatible with VA, VB,
% VD and also reading *.SR files (PhoenixProtocols).
%
% USAGE:
% ascout = read_ASCII(ascin)
% where ascin is the char content of either
% hdr.CSASeriesHeaderInfo.MrPhoenixProtocol or
% hdr.CSASeriesHeaderInfo.MrProtocol according to Siemens software version,
% and asc is a matlab structure containing the information contained in the
% ASCCONV BEGIN - ASCCONV END part of the DICOM header.
%
% NOTE: slightly different formatting from read_ascconv in spm_dicom_convert

% only for debugging purpose
ENABLE_DEBUG = false;

% extract the portion between ### ASCCONV BEGIN ### and ### ASCCONV END ###
ascin = regexprep(ascin,'^.*### ASCCONV BEGIN [^#]*###(.*)### ASCCONV END ###.*$','$1');

% split ascin into lines
% asclines = strsplit(ascin,'\n'); % ok for Matlab 2013 and later versions
% replaced by textscan (with delimiter '\n') for compatibility with
% Matlab 2012 and earlier :/...
tmp = textscan(ascin,'%s','delimiter','\n');
asclines = tmp{1};

% do a first cleaning pass over the data:
% - replace indexes [] by () + increment index +1 (C>Matlab indexing)
% - replace double "" and single " by single '
% - replace hexadecimal values (e.g. "0x01") by decimal value
% - delete lines where a field starts or ends by "_"
% - delete end of lines after "#"
asclines = regexprep(asclines,{'\[([0-9]*)\]','["]+','^([^"]*)0x([0-9a-fA-F]*)',' #.*','^.*\._.*$'},{'($1+1)','''','$1hex2dec(''$2'')','',''});

% initialise variables
ascout = [];
clinenum = 1;

while clinenum<length(asclines)

    if ENABLE_DEBUG, fprintf('Line #%d - [%s]\n', clinenum, asclines{clinenum}); end
    if ~isempty(asclines{clinenum})
        try
            [tlhrh,tmp] = regexp(asclines{clinenum}, '(?:=)+', 'split', 'match'); %#ok<*ASGLU>
            % first process the name of the parameter (hdrnam)
            % split into fields
            [hdrnam,tmp] = regexp(tlhrh{1}, '(?:\.)+', 'split', 'match');
            % check whether any field starts with a number and if so,
            % replace it by "index<number>..."
            isfirstdigit = cellfun(@(x) isstrprop(x(1),'digit'),hdrnam);
            if any(isfirstdigit)
                for cfield=1:length(hdrnam)
                    if isfirstdigit(cfield)
                        hdrnam{cfield} = ['index' hdrnam{cfield}];
                    end
                end
            end
            % concatenate back the fields
            hdrnam = sprintf('.%s', hdrnam{:});
            % now process the value (hdrval)
            hdrval = strtrim(tlhrh{2});
            % and eval the whole line
            eval(sprintf('ascout%s = %s;', hdrnam, hdrval));
            if ENABLE_DEBUG, fprintf('ascout%s = %s;\n', hdrnam, hdrval); end
        catch
            fprintf('AscConv: Error evaluating [ %s; ]\n', asclines{clinenum});
        end
    end
    clinenum = clinenum+1;
end

author: spm


text:
% Copyright 2009 - 2010 The MathWorks, Inc.
function ObjTrack01(position)
%#eml
% Figure setup
numPts = 300;
figure;hold;grid;

% Kalman filter loop
for idx = 1: numPts
    % Get the input data
    z = position(:,idx);

    % Use Kalman filter to estimate the location
    y = kalman01(z);

    % Plot the results
    plot_trajectory(z,y);
end
hold;
end

author: Sable

repo: mcbench-benchmarks


text:
classdef MetaboliteUtil
    methods (Static)
        function localizedMetabolites = localizeMetabolites(sim, metabolites)
            mets = sim.state('Metabolite');
            nMets = numel(mets.wholeCellModelIDs);
            nComp = sim.compartment.count;
            cIdx = sim.compartment.cytosolIndexs;
            mIdx = sim.compartment.membraneIndexs;
            hMetIdxs = mets.hydrophobicIndexs;

            localizedMetabolites = zeros(nMets, nComp);
            localizedMetabolites(:, cIdx) = metabolites;
            localizedMetabolites(hMetIdxs, mIdx) = localizedMetabolites(hMetIdxs, cIdx);
            localizedMetabolites(hMetIdxs, cIdx) = 0;
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function test_suite = test_polygonSecondAreaMoments
%TEST_POLYGONSECONDAREAMOMENTS Tests for polygonSecondAreaMoments.
%
%   Example
%       runtests('test_polygonSecondAreaMoments.m')
%
%   References
%       Expected (_exp) values are taken from:
%       https://en.wikipedia.org/wiki/List_of_second_moments_of_area
%       2020 - Bedford & Liechti - Mechanics of Materials
%           DOI: 10.1007/978-3-030-22082-2
%

% ------
% Author: oqilipo
% E-mail: N/A
% Created: 2022-12-21, using MATLAB 9.13.0.2080170 (R2022b) Update 1
% Copyright 2022

test_suite = functiontests(localfunctions);

function testCircle(testCase)
r = randi([1 100]);
circle = circleToPolygon([0 0 r], round(2*pi*r*1e2));

[Ixx, Iyy, Ixy] = polygonSecondAreaMoments(circle);
Ixx_exp = pi/4*r^4;
Iyy_exp = Ixx_exp;
Ixy_exp = 0;

testCase.assertEqual(Ixx, Ixx_exp, 'AbsTol',3*1e-1);
testCase.assertEqual(Iyy, Iyy_exp, 'AbsTol',3*1e-1);
testCase.assertEqual(Ixy, Ixy_exp, 'AbsTol',1e-1);

function testAnnulus(testCase)
r2 = randi([50 100]);
r1 = r2 - randi([1 49]);
outerCircle = circleToPolygon([0 0 r2], round(2*pi*r2*1e2));
innerCircle = flipud(circleToPolygon([0 0 r1], round(2*pi*r1*1e2)));
annulus = {outerCircle; innerCircle};

[Ixx, Iyy, Ixy] = polygonSecondAreaMoments(annulus);
Ixx_exp = pi/4*(r2^4-r1^4);
Iyy_exp = Ixx_exp;
Ixy_exp = 0;

testCase.assertEqual(Ixx, Ixx_exp, 'AbsTol',2*1e-1);
testCase.assertEqual(Iyy, Iyy_exp, 'AbsTol',2*1e-1);
testCase.assertEqual(Ixy, Ixy_exp, 'AbsTol',1e-1);

function testCircularSector(testCase)
r = randi([1 100]);
circle = circleToPolygon([0 0 r] , round(2*pi*r*1e2));
theta = (pi/2-0.2)*rand+0.1;
xPosLine = reverseLine(createLine(theta/2));
xNegLine = createLine(-theta/2);

semiCircle = clipPolygonHP(circle, xPosLine);
circularSector = clipPolygonHP(semiCircle, xNegLine);

[Ixx, Iyy, Ixy] = polygonSecondAreaMoments(circularSector);
Ixx_exp = (theta-sin(theta))*r^4/8;
Iyy_exp = (theta+sin(theta))*r^4/8;
Ixy_exp = 0;

testCase.assertEqual(Ixx, Ixx_exp, 'AbsTol',1e-1);
testCase.assertEqual(Iyy, Iyy_exp, 'AbsTol',1e-1);
testCase.assertEqual(Ixy, Ixy_exp, 'AbsTol',1e-7);

function testCenteredSemiCircle(testCase)
r = randi([1 100]);
circle = circleToPolygon([0 0 r], round(2*pi*r*1e2));
semiCircle = clipPolygonHP(circle, [0 0 1 0]);
centeredSemiCircle = semiCircle - polygonCentroid(semiCircle);

[Ixx, Iyy, Ixy] = polygonSecondAreaMoments(centeredSemiCircle);
Ixx_exp = (pi/8-8/(9*pi))*r^4;
Iyy_exp = pi/8*r^4;
Ixy_exp = 0;

testCase.assertEqual(Ixx, Ixx_exp, 'AbsTol',2*1e-1);
testCase.assertEqual(Iyy, Iyy_exp, 'AbsTol',2*1e-1);
testCase.assertEqual(Ixy, Ixy_exp, 'AbsTol',1e-7);

function testSemiCircle(testCase)
r = randi([1 100]);
circle = circleToPolygon([0 0 r], round(2*pi*r*1e2));
semiCircle = clipPolygonHP(circle, [0 0 1 0]);

[Ixx, Iyy, Ixy] = polygonSecondAreaMoments(semiCircle);
Ixx_exp = pi/8*r^4;
Iyy_exp = Ixx_exp;
Ixy_exp = 0;

testCase.assertEqual(Ixx, Ixx_exp, 'AbsTol',2*1e-1);
testCase.assertEqual(Iyy, Iyy_exp, 'AbsTol',2*1e-1);
testCase.assertEqual(Ixy, Ixy_exp, 'AbsTol',1e-7);

function testQuarterCircle(testCase)
r = randi([1 100]);
circle = circleToPolygon([0 0 r], round(2*pi*r*1e2));
semiCircle = clipPolygonHP(circle, [0 0 1 0]);
quarterCircle = clipPolygonHP(semiCircle, [0 0 0 -1]);

[Ixx, Iyy, Ixy] = polygonSecondAreaMoments(quarterCircle);
Ixx_exp = pi/16*r^4;
Iyy_exp = pi/16*r^4;
Ixy_exp = r^4/8;

testCase.assertEqual(Ixx, Ixx_exp, 'AbsTol',1e-1);
testCase.assertEqual(Iyy, Iyy_exp, 'AbsTol',1e-1);
testCase.assertEqual(Ixy, Ixy_exp, 'AbsTol',1e-1);

function testCenteredQuarterCircle(testCase)
r = randi([1 100]);
circle = circleToPolygon([0 0 r], round(2*pi*r*1e2));
semiCircle = clipPolygonHP(circle, [0 0 1 0]);
quarterCircle = clipPolygonHP(semiCircle, [0 0 0 -1]);
centeredQuarterCircle = quarterCircle - polygonCentroid(quarterCircle);

[Ixx, Iyy, Ixy] = polygonSecondAreaMoments(centeredQuarterCircle); %#ok<ASGLU>
Ixx_exp = (pi/16-4/(9*pi))*r^4;
Iyy_exp = Ixx_exp;
% Ixy_exp = ?;

testCase.assertEqual(Ixx, Ixx_exp, 'AbsTol',1e-1);
testCase.assertEqual(Iyy, Iyy_exp, 'AbsTol',1e-1);

function testEllipse(testCase)
a = randi([51 100]);
b = randi([1 50]);
ellipse = ellipseToPolygon([0 0 a b], 4*a*sqrt(1-b^2/a^2)*1e2);

[Ixx, Iyy, Ixy] = polygonSecondAreaMoments(ellipse);
Ixx_exp = pi/4*a*b^3;
Iyy_exp = pi/4*a^3*b;
Ixy_exp = 0;

testCase.assertEqual(Ixx, Ixx_exp, 'AbsTol',6*1e-1);
testCase.assertEqual(Iyy, Iyy_exp, 'AbsTol',6*1e-1);
testCase.assertEqual(Ixy, Ixy_exp, 'AbsTol',1e-1);

function testRectangle(testCase)
b = randi([51 100]);
h = randi([1 50]);
rectangle = [0 0; b 0; b h; 0 h];

[Ixx, Iyy, Ixy] = polygonSecondAreaMoments(rectangle);
Ixx_exp = b*h^3/3;
Iyy_exp = b^3*h/3;
Ixy_exp = b^2*h^2/4;

testCase.assertEqual(Ixx, Ixx_exp);
testCase.assertEqual(Iyy, Iyy_exp);
testCase.assertEqual(Ixy, Ixy_exp);

function testCenteredRectangle(testCase)
b = randi([51 100]);
h = randi([1 50]);
centeredRectangle = [b/2 -h/2; b/2 h/2; -b/2 h/2; -b/2 -h/2];

[Ixx, Iyy, Ixy] = polygonSecondAreaMoments(centeredRectangle);
Ixx_exp = b*h^3/12;
Iyy_exp = b^3*h/12;
Ixy_exp = 0;

testCase.assertEqual(Ixx, Ixx_exp);
testCase.assertEqual(Iyy, Iyy_exp);
testCase.assertEqual(Ixy, Ixy_exp);

function testHollowRectangle(testCase)
b = randi([51 100]);
h = randi([10 50]);
b1 = b - randi([10 20]);
h1 = h - randi([1 8]);
outerRectangle = [b/2 -h/2; b/2 h/2; -b/2 h/2; -b/2 -h/2];
innerRectangle = [-b1/2 -h1/2; -b1/2 h1/2; b1/2 h1/2; b1/2 -h1/2];
hollowRectangle = {outerRectangle; innerRectangle};

[Ixx, Iyy, Ixy] = polygonSecondAreaMoments(hollowRectangle);
Ixx_exp = (b*h^3 - b1*h1^3)/12;
Iyy_exp = (b^3*h - b1^3*h1)/12;
Ixy_exp = 0;

testCase.assertEqual(Ixx, Ixx_exp);
testCase.assertEqual(Iyy, Iyy_exp);
testCase.assertEqual(Ixy, Ixy_exp);

function testCenteredTriangle(testCase)
b = randi([51 100]);
h = randi([51 100]);
a = randi([1 25]);
triangle = [0 0; b 0; a h];
centeredTriangle = triangle - polygonCentroid(triangle);

[Ixx, Iyy, Ixy] = polygonSecondAreaMoments(centeredTriangle);
Ixx_exp = b*h^3/36;
Iyy_exp = (b^3*h-b^2*h*a+b*h*a^2)/36;
Ixy_exp = -b*h^2/72*(b-2*a);

testCase.assertEqual(Ixx, Ixx_exp, 'AbsTol',1e-8);
testCase.assertEqual(Iyy, Iyy_exp, 'AbsTol',1e-8);
testCase.assertEqual(Ixy, Ixy_exp, 'AbsTol',1e-8);

function testTriangle(testCase)
b = randi([51 100]);
h = randi([51 100]);
a = randi([1 25]);
triangle = [0 0; b 0; a h];

[Ixx, Iyy, Ixy] = polygonSecondAreaMoments(triangle);
Ixx_exp = b*h^3/12;
Iyy_exp = (b^3*h+b^2*h*a+b*h*a^2)/12;
Ixy_exp = b*h^2/24*(b+2*a);

testCase.assertEqual(Ixx, Ixx_exp, 'AbsTol',1e-8);
testCase.assertEqual(Iyy, Iyy_exp, 'AbsTol',1e-8);
testCase.assertEqual(Ixy, Ixy_exp, 'AbsTol',1e-8);

function testlSection(testCase)
L = randi([51 100]);
t = randi([1 25]);
lSection = [0 0; L 0; L t; t t; t L; 0 L];
lSection = lSection - polygonCentroid(lSection);

[Ixx, Iyy, Ixy] = polygonSecondAreaMoments(lSection);
Ixx_exp = (t*(5*L^2-5*L*t+t^2)*(L^2-L*t+t^2))/(12*(2*L-t));
Iyy_exp = Ixx_exp;
Ixy_exp = (L^2*t*(L-t)^2)/(4*(t-2*L));

testCase.assertEqual(Ixx, Ixx_exp, 'AbsTol',1e-8);
testCase.assertEqual(Iyy, Iyy_exp, 'AbsTol',1e-8);
testCase.assertEqual(Ixy, Ixy_exp, 'AbsTol',1e-8);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% GRADIENT DESCENT function.
% Iteratively optimizes theta model parameters.
function [theta, J, exit_flag] = gradient_descent(X, y, theta, lambda)
    % X - training set.
    % y - training output values.
    % theta - model parameters.
    % lambda - regularization parameter.

    % Set Options
    options = optimset('GradObj', 'on', 'MaxIter', 400);

    % Optimize
    [theta, J, exit_flag] = fminunc(@(t)(gradient_callback(X, y, t, lambda)), theta, options);
end

author: trekhleb

repo: machine-learning-octave


text:
function varargout = ma_InsertMassDumpGUI(varargin)
% MA_INSERTMASSDUMPGUI MATLAB code for ma_InsertMassDumpGUI.fig
%      MA_INSERTMASSDUMPGUI, by itself, creates a new MA_INSERTMASSDUMPGUI or raises the existing
%      singleton*.
%
%      H = MA_INSERTMASSDUMPGUI returns the handle to a new MA_INSERTMASSDUMPGUI or the handle to
%      the existing singleton*.
%
%      MA_INSERTMASSDUMPGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MA_INSERTMASSDUMPGUI.M with the given input arguments.
%
%      MA_INSERTMASSDUMPGUI('Property','Value',...) creates a new MA_INSERTMASSDUMPGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before ma_InsertMassDumpGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to ma_InsertMassDumpGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help ma_InsertMassDumpGUI

% Last Modified by GUIDE v2.5 19-Jun-2014 19:45:21

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @ma_InsertMassDumpGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @ma_InsertMassDumpGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before ma_InsertMassDumpGUI is made visible.
function ma_InsertMassDumpGUI_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to ma_InsertMassDumpGUI (see VARARGIN)

% Choose default command line output for ma_InsertMassDumpGUI
handles.output = hObject;
handles.ma_MainGUI = varargin{1};

% Update handles structure
guidata(hObject, handles);

if(length(varargin)>1)
    event = varargin{2};
    populateGUIWithEvent(handles, event);
end

% GUI Setup
setPropellantNameLabels(handles);
populateThrusterCombo(handles, handles.thrustersCombo);
massDumpTypeCombo_Callback(handles.massDumpTypeCombo, [], handles);
thrustersCombo_Callback(handles.thrustersCombo, [], handles);

% UIWAIT makes ma_InsertMassDumpGUI wait for user response (see UIRESUME)
uiwait(handles.ma_InsertMassDumpGUI);

function setPropellantNameLabels(handles)
    maData = getappdata(handles.ma_MainGUI,'ma_data');
    propNames = maData.spacecraft.propellant.names;

    set(handles.propName1Text,'String', propNames{1}, 'TooltipString', propNames{1});
    set(handles.propName2Text,'String', propNames{2}, 'TooltipString', propNames{2});
    set(handles.propName3Text,'String', propNames{3}, 'TooltipString', propNames{3});

function populateThrusterCombo(handles, hThrustersCombo)
    maData = getappdata(handles.ma_MainGUI,'ma_data');
    thrusters = maData.spacecraft.thrusters;

    thrusterComboStr = cell(length(thrusters),1);
    for(i=1:length(thrusters))
        thrusterComboStr{i} = thrusters{i}.name;
    end

    set(hThrustersCombo,'String',thrusterComboStr);

function populateGUIWithEvent(handles, event)
    maData = getappdata(handles.ma_MainGUI,'ma_data');

    set(handles.titleLabel, 'String', 'Edit Mass Dump');
    set(handles.ma_InsertMassDumpGUI, 'Name', 'Edit Mass Dump');
    set(handles.nameText, 'String', event.name);

    type = event.subType;
    switch type
        case 'basic'
            value = findValueFromComboBox('Mass Dump', handles.massDumpTypeCombo);
            set(handles.massDumpTypeCombo,'value',value);

            dmasses = event.dumpValue;

            set(handles.dryMassText,'String', fullAccNum2Str(dmasses(1)));
            set(handles.fuelOxMassText,'String', fullAccNum2Str(dmasses(2)));
            set(handles.monopropMassText,'String', fullAccNum2Str(dmasses(3)));
            set(handles.xenonMassText,'String', fullAccNum2Str(dmasses(4)));
        case 'delta-v'
            value = findValueFromComboBox('Delta-V Dump', handles.massDumpTypeCombo);
            set(handles.massDumpTypeCombo,'value',value);

            deltav = event.dumpValue;
            set(handles.deltavText,'String',fullAccNum2Str(deltav));

            thrusters = maData.spacecraft.thrusters;
            thrusterVal = 1;
            for(i=1:length(thrusters)) %#ok<*NO4LP>
                thrusterFromArr = thrusters{i};
                if(thrusterFromArr.id == event.thruster.id)
                    thrusterVal = i;
                    break;
                end
            end

            set(handles.thrustersCombo,'value',thrusterVal);
    end

% --- Outputs from this function are returned to the command line.
function varargout = ma_InsertMassDumpGUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
if(isempty(handles))
    varargout{1} = [];
else
    maData = getappdata(handles.ma_MainGUI,'ma_data');

    name = get(handles.nameText, 'String');

    contents = cellstr(get(handles.massDumpTypeCombo,'String'));
    selected = contents{get(handles.massDumpTypeCombo,'Value')};

    switch selected
        case 'Mass Dump'
            subType = 'basic';

            dryMass = str2double(get(handles.dryMassText,'String'));
            fuelOxMass = str2double(get(handles.fuelOxMassText,'String'));
            monopropMass = str2double(get(handles.monopropMassText,'String'));
            xenonMass = str2double(get(handles.xenonMassText,'String'));

            massDumpValue = [dryMass fuelOxMass monopropMass xenonMass];

            thruster = [];
        case 'Delta-V Dump'
            subType = 'delta-v';

            deltav = str2double(get(handles.deltavText,'String'));

            massDumpValue = [deltav];

            thrusterNum = get(handles.thrustersCombo,'value');
            thruster = maData.spacecraft.thrusters{thrusterNum};
    end

    varargout{1} = ma_createMassDump(name, subType, massDumpValue, thruster);
    close(handles.ma_InsertMassDumpGUI);
end

% --- Executes on button press in saveCloseButton.
function saveCloseButton_Callback(hObject, eventdata, handles)
% hObject    handle to saveCloseButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    errMsg = validateInputs(handles);

    if(isempty(errMsg))
        uiresume(handles.ma_InsertMassDumpGUI);
    else
        msgbox(errMsg,'Errors were found while inserting a mass dump.','error');
    end

function errMsg = validateInputs(handles)
    contents = cellstr(get(handles.massDumpTypeCombo,'String'));
    selected = contents{get(handles.massDumpTypeCombo,'Value')};

    errMsg = {};

    switch selected
        case 'Mass Dump'
            dryMass = str2double(get(handles.dryMassText,'String'));
            enteredStr = get(handles.dryMassText,'String');
            numberName = 'Dry Mass to Dump';
            lb = 0;
            ub = Inf;
            isInt = false;
            errMsg = validateNumber(dryMass, numberName, lb, ub, isInt, errMsg, enteredStr);

            fuelOxMass = str2double(get(handles.fuelOxMassText,'String'));
            enteredStr = get(handles.fuelOxMassText,'String');
            numberName = 'Fuel/Ox Mass to Dump';
            lb = 0;
            ub = Inf;
            isInt = false;
            errMsg = validateNumber(fuelOxMass, numberName, lb, ub, isInt, errMsg, enteredStr);

            monopropMass = str2double(get(handles.monopropMassText,'String'));
            enteredStr = get(handles.monopropMassText,'String');
            numberName = 'Monoprop Mass to Dump';
            lb = 0;
            ub = Inf;
            isInt = false;
            errMsg = validateNumber(monopropMass, numberName, lb, ub, isInt, errMsg, enteredStr);

            xenonMass = str2double(get(handles.xenonMassText,'String'));
            enteredStr = get(handles.xenonMassText,'String');
            numberName = 'Xenon Mass to Dump';
            lb = 0;
            ub = Inf;
            isInt = false;
            errMsg = validateNumber(xenonMass, numberName, lb, ub, isInt, errMsg, enteredStr);

        case 'Delta-V Dump'
            deltav = str2double(get(handles.deltavText,'String'));
            enteredStr = get(handles.deltavText,'String');
            numberName = 'Delta-V to Dump';
            lb = 0;
            ub = Inf;
            isInt = false;
            errMsg = validateNumber(deltav, numberName, lb, ub, isInt, errMsg, enteredStr);
    end

% --- Executes on button press in cancelButton.
function cancelButton_Callback(hObject, eventdata, handles)
% hObject    handle to cancelButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    close(handles.ma_InsertMassDumpGUI);

function nameText_Callback(hObject, eventdata, handles)
% hObject    handle to nameText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of nameText as text
%        str2double(get(hObject,'String')) returns contents of nameText as a double

% --- Executes during object creation, after setting all properties.
function nameText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to nameText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in massDumpTypeCombo.
function massDumpTypeCombo_Callback(hObject, eventdata, handles)
% hObject    handle to massDumpTypeCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns massDumpTypeCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from massDumpTypeCombo
    contents = cellstr(get(hObject,'String'));
    selected = contents{get(hObject,'Value')};

    switch selected
        case 'Mass Dump'
            set(handles.dryMassText,'Enable','on');
            set(handles.fuelOxMassText,'Enable','on');
            set(handles.monopropMassText,'Enable','on');
            set(handles.xenonMassText,'Enable','on');
            set(handles.deltavText,'Enable','off');
            set(handles.thrustersCombo,'Enable','off');
        case 'Delta-V Dump'
            set(handles.dryMassText,'Enable','off');
            set(handles.fuelOxMassText,'Enable','off');
            set(handles.monopropMassText,'Enable','off');
            set(handles.xenonMassText,'Enable','off');
            set(handles.deltavText,'Enable','on');
            set(handles.thrustersCombo,'Enable','on');
    end

% --- Executes during object creation, after setting all properties.
function massDumpTypeCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to massDumpTypeCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function fuelOxMassText_Callback(hObject, eventdata, handles)
% hObject    handle to fuelOxMassText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of fuelOxMassText as text
%        str2double(get(hObject,'String')) returns contents of fuelOxMassText as a double
newInput = get(hObject,'String');
newInput = attemptStrEval(newInput);
set(hObject,'String', newInput);

% --- Executes during object creation, after setting all properties.
function fuelOxMassText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to fuelOxMassText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function monopropMassText_Callback(hObject, eventdata, handles)
% hObject    handle to monopropMassText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of monopropMassText as text
%        str2double(get(hObject,'String')) returns contents of monopropMassText as a double
newInput = get(hObject,'String');
newInput = attemptStrEval(newInput);
set(hObject,'String', newInput);

% --- Executes during object creation, after setting all properties.
function monopropMassText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to monopropMassText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function xenonMassText_Callback(hObject, eventdata, handles)
% hObject    handle to xenonMassText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of xenonMassText as text
%        str2double(get(hObject,'String')) returns contents of xenonMassText as a double
newInput = get(hObject,'String');
newInput = attemptStrEval(newInput);
set(hObject,'String', newInput);

% --- Executes during object creation, after setting all properties.
function xenonMassText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to xenonMassText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function dryMassText_Callback(hObject, eventdata, handles)
% hObject    handle to dryMassText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of dryMassText as text
%        str2double(get(hObject,'String')) returns contents of dryMassText as a double
newInput = get(hObject,'String');
newInput = attemptStrEval(newInput);
set(hObject,'String', newInput);

% --- Executes during object creation, after setting all properties.
function dryMassText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to dryMassText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function deltavText_Callback(hObject, eventdata, handles)
% hObject    handle to deltavText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of deltavText as text
%        str2double(get(hObject,'String')) returns contents of deltavText as a double
newInput = get(hObject,'String');
newInput = attemptStrEval(newInput);
set(hObject,'String', newInput);

% --- Executes during object creation, after setting all properties.
function deltavText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to deltavText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in thrustersCombo.
function thrustersCombo_Callback(hObject, eventdata, handles)
% hObject    handle to thrustersCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns thrustersCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from thrustersCombo
    maData = getappdata(handles.ma_MainGUI,'ma_data');
    propNames = maData.spacecraft.propellant.names;
    thrusters = maData.spacecraft.thrusters;

    thruster = thrusters{get(handles.thrustersCombo,'Value')};
    propTypeIndex = ma_getPropTypeIndexFromInternalPropName(thruster.propType);
    set(handles.fuelTypeBeingDumpedLabel,'String',propNames{propTypeIndex});

% --- Executes during object creation, after setting all properties.
function thrustersCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to thrustersCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on key press with focus on ma_InsertMassDumpGUI or any of its controls.
function ma_InsertMassDumpGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to ma_InsertMassDumpGUI (see GCBO)
% eventdata  structure with the following fields (see FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            saveCloseButton_Callback(handles.saveCloseButton, [], handles);
        case 'enter'
            saveCloseButton_Callback(handles.saveCloseButton, [], handles);
        case 'escape'
            close(handles.ma_InsertMassDumpGUI);
    end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [hFig, iDS, iFig] = view_timefreq(TimefreqFile, DisplayMode, RowName, isNewFigure, Function)
% VIEW_TIMEFREQ: Display times frequency maps in a new figure.
%
% USAGE: [hFig, iDS, iFig] = view_timefreq(TimefreqFile, DisplayMode='SingleSensor', RowName=[], isNewFigure=0, Function=[])
%
% INPUT:
%     - TimefreqFile : Path to time-frequency file to visualize
%     - DisplayMode  : {'SingleSensor', 'AllSensors', '2DLayout', '2DLayoutOpt'}
%     - RowName      : Name of the row to display from the input timefreq file
%     - isNewFigure  : If 1, force the creation of a new figure
%
% OUTPUT :
%     - hFig : Matlab handle to the figure that was created or updated
%     - iDS  : DataSet index in the GlobalData variable
%     - iFig : Indice of returned figure in the GlobalData(iDS).Figure array

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2010-2014

%% ===== INITIALIZATION =====
% GlobalData : create if not existing yet
global GlobalData;
% Parse inputs
if (nargin < 2)
    DisplayMode = 'SingleSensor';
end
if (nargin < 3) || isempty(RowName)
    RowName = [];
elseif iscell(RowName)
    RowName = RowName{1};
end
if (nargin < 4) || isempty(isNewFigure) || (isNewFigure == 0)
    CreateMode = '';
else
    CreateMode = 'AlwaysCreate';
end
if (nargin < 5) || isempty(Function)
    Function = [];
end

if ~isempty(strfind(lower(TimefreqFile), 'spike_field_coherence')) ...
        || ~isempty(strfind(lower(TimefreqFile), 'noise_correlation')) ...
        || ~isempty(strfind(lower(TimefreqFile), 'rasterplot'))...
        || ~isempty(strfind(lower(TimefreqFile), 'spiking_phase_locking'))
    isEphysFile = 1;
    GlobalData.UserFrequencies.HideFreqPanel = 1;
else
    isEphysFile = 0;
end

%% ===== GET ALL ACCESSIBLE DATA =====
% Get study
[sStudy, iStudy, iItem, DataType, sTimefreq] = bst_get('AnyFile', TimefreqFile);
if isempty(sStudy)
    error('File is not registered in database.');
end

%% ===== LOAD TIME-FREQUENCY FILE =====
bst_progress('start', 'View time-frequency map', 'Loading data...');
% Load file
[iDS, iTimefreq] = bst_memory('LoadTimefreqFile', TimefreqFile);
if isempty(iDS)
    % error('Cannot load timefreq file.');
    hFig = [];
    iFig = [];
    return
end
% If there is only one row available: for DisplayMode to 'SingleSensor'
if (length(GlobalData.DataSet(iDS).Timefreq(iTimefreq).RowNames) == 1)
    DisplayMode = 'SingleSensor';
end
% Detect modality
Modality = GlobalData.DataSet(iDS).Timefreq(iTimefreq).Modality;
AllModalities = GlobalData.DataSet(iDS).Timefreq(iTimefreq).AllModalities;
LayoutRows = [];
% Get the all row names for this file
AllRows = figure_timefreq('GetRowNames', GlobalData.DataSet(iDS).Timefreq(iTimefreq).RefRowNames, GlobalData.DataSet(iDS).Timefreq(iTimefreq).RowNames);
% Check display type
if any(strcmpi(DisplayMode, {'2DLayout', '2DLayoutOpt'}))
    % Sensor types that cannot be displayed
    if ~isempty(Modality) && ~ismember(Modality, {'MEG','MEG GRAD','MEG MAG','EEG','NIRS'})
        disp(['TIMEFREQ> Error: Cannot display 2DLayout for modality "' Modality '", using AllSensors instead.']);
        DisplayMode = 'AllSensors';
    % If there are multiple modalities available in the file
    elseif isempty(Modality) && ~isempty(AllModalities)
        % Check for main sensor categories
        DispMod = intersect(AllModalities, {'MEG','MEG GRAD','MEG MAG','EEG','ECOG','NIRS'});
        % No available display types
        if isempty(DispMod)
            Modality = [];
        % If there are multiple modalities but one only that can be displayed as a topography
        elseif (length(DispMod) == 1)
            Modality = DispMod{1};
        % Else: ask user what to display
        elseif ~isempty(DispMod)
            res = java_dialog('question', ['This file contains multiple sensor types.' 10 'Which modality would you like to display?'], 'Select sensor type', [], DispMod);
            if isempty(res) || strcmpi(res, 'Cancel')
                return;
            else
                Modality = res;
            end
        end
        % Get the displayed rows only
        if ~isempty(Modality)
            iChan = channel_find(GlobalData.DataSet(iDS).Channel, Modality);
            LayoutRows = intersect(AllRows, {GlobalData.DataSet(iDS).Channel(iChan).Name});
        end
    end
    % No solution found
    if isempty(Modality)
        disp('TIMEFREQ> Error: Cannot display 2DLayout for this file, using AllSensors instead.');
        DisplayMode = 'AllSensors';
    end
end
% Make sure that the required row names are available
if ~isempty(RowName) && ischar(RowName) && ~ismember(RowName, AllRows)
    error(['There is no entry "' RowName '" in this file.']);
end

%% ===== CREATE A NEW FIGURE =====
% Prepare FigureId structure
FigureId.Type     = 'Timefreq';
FigureId.SubType  = DisplayMode;
FigureId.Modality = Modality;
% Create TimeSeries figure
[hFig, iFig] = bst_figures('CreateFigure', iDS, FigureId, CreateMode, sTimefreq.FileName);
if isempty(hFig)
    error('Cannot create figure');
end

%% ===== INITIALIZE FIGURE =====
% Configure app data
setappdata(hFig, 'DataFile',     GlobalData.DataSet(iDS).DataFile);
setappdata(hFig, 'StudyFile',    GlobalData.DataSet(iDS).StudyFile);
setappdata(hFig, 'SubjectFile',  GlobalData.DataSet(iDS).SubjectFile);
% Static dataset
setappdata(hFig, 'isStatic',     (GlobalData.DataSet(iDS).Timefreq(iTimefreq).NumberOfSamples <= 2));
%setappdata(hFig, 'isStaticFreq', (size(GlobalData.DataSet(iDS).Timefreq(iTimefreq).TF,3) <= 1));
setappdata(hFig, 'isStaticFreq', 0);
% Get figure data
TfInfo = getappdata(hFig, 'Timefreq');
% Create options structure
TfInfo.FileName    = sTimefreq.FileName;
TfInfo.Comment     = sTimefreq.Comment;
TfInfo.DisplayMode = DisplayMode;
TfInfo.iFreqs      = [];
if ismember(TfInfo.DisplayMode, {'2DLayout', '2DLayoutOpt', 'AllSensors'})
    TfInfo.RowName = LayoutRows;
elseif ~isempty(RowName)
    TfInfo.RowName = RowName;
elseif iscell(AllRows)
    TfInfo.RowName = AllRows{1};
else
    TfInfo.RowName = AllRows(1);
end
% Default function
if isEphysFile
    TfInfo.Function = 'power';
    TfInfo.DisplayMeasure = 0;
    TfInfo.DisableHideEdgeEffects = 1;
elseif ~isempty(Function)
    TfInfo.Function = Function;
else
    TfInfo.Function = process_tf_measure('GetDefaultFunction', GlobalData.DataSet(iDS).Timefreq(iTimefreq));
end
% EPhys specific fields
if ~isempty(strfind(lower(TimefreqFile), 'noise_correlation'))
    DataMat = in_bst_data(TimefreqFile, 'NeuronNames');
    TfInfo.NeuronNames = DataMat.NeuronNames;
elseif ~isempty(strfind(lower(TimefreqFile), 'rasterplot'))
    TfInfo.DisplayAsDots = 1;
    TfInfo.DisableSmoothDisplay = 1;
elseif ~isempty(strfind(lower(TimefreqFile), 'spiking_phase_locking'))
    TfInfo.DisplayAsPhase = 1;
end
% Set figure data
setappdata(hFig, 'Timefreq', TfInfo);
% Add colormap (stat or timefreq)
if ~isempty(GlobalData.DataSet(iDS).Timefreq(iTimefreq).ColormapType)
    ColormapType = GlobalData.DataSet(iDS).Timefreq(iTimefreq).ColormapType;
elseif strcmpi(file_gettype(TimefreqFile), 'ptimefreq') || ~isempty(strfind(lower(TimefreqFile), 'zscore')) || ~isempty(strfind(lower(TimefreqFile), 'ersd'))
    ColormapType = 'stat2';
else
    ColormapType = 'timefreq';
end
% Display units
DisplayUnits = GlobalData.DataSet(iDS).Timefreq(iTimefreq).DisplayUnits;
% Add colormap to figure
bst_colormaps('AddColormapToFigure', hFig, ColormapType, DisplayUnits);
% Display options panel
isDisplayTab = ~strcmpi(TfInfo.Function, 'other') || (~isempty(TfInfo.RowName) && (ischar(TfInfo.RowName) || iscell(TfInfo.RowName)));
if isDisplayTab
    gui_brainstorm('ShowToolTab', 'Display');
end

%% ===== PLOT DATA =====
% Plot time-freq map
figure_timefreq('UpdateFigurePlot', hFig);
% Update figure name
bst_figures('UpdateFigureName', hFig);

%% =====  DISPLAY FIGURE =====
% Reset selected figure, so that the following figure selection can update the timefreq panel
bst_figures('SetCurrentFigure', [], 'TF');
% Set selected figure
bst_figures('SetCurrentFigure', hFig, 'TF');
% Select display options
if isDisplayTab
    panel_display('UpdatePanel', hFig);
end
% Set figure visible
set(hFig, 'Visible', 'on');
bst_progress('stop');

end

author: brainstorm-tools


text:
function [cfg, movement] = ft_detect_movement(cfg, data)

% FT_SACCADE_DETECTION performs detection of movements such as saccades and
% microsaccades, but also joystick movements, from time series data over multiple
% trials. Different methods for detecting movements are implemented, which are
% described in detail below:
%
% VELOCITY2D - detects micro/saccades using a two-dimensional (2D) velocity according
% to "Engbert R, Kliegl R (2003) Vision Res 43:1035-1045". The vertical and the
% horizontal eyetracker time series (for one eye) are transformed into velocities and
% microsaccades are indentified as "outlier" eye movements that exceed a given
% threshold for velocity and duration. This method has the additional options
%     cfg.velocity2D.kernel   = vector 1 x nsamples, kernel to compute velocity (default = [1 1 0 -1 -1].*(data.fsample/6);
%     cfg.velocity2D.demean   = 'no' or 'yes', whether to apply centering correction (default = 'yes')
%     cfg.velocity2D.mindur   = minimum microsaccade durantion in samples (default = 3);
%     cfg.velocity2D.velthres = threshold for velocity outlier detection (default = 6);
%
% CLUSTERING - detects movements according to "Otero-Millan et al., (2014) J Vis 14".
%
% Use as
%   [cfg, movement] = ft_detect_movement(cfg, data)
% where the input data should be organised in a structure as obtained from the
% FT_PREPROCESSING function.
%
% The configuration can contain the following options
%   cfg.method  = string representing the method for movement detection
%                 'velocity2D' detects microsaccades using the 2D velocity
%                 'clustering' use unsupervised clustering method to detect microsaccades
%   cfg.channel = Nx1 cell-array with selection of channels, see FT_CHANNELSELECTION for details, (default = 'all')
%   cfg.trials  = 'all' or a selection given as a 1xN vector (default = 'all')
%
% The output argument "movement" is a Nx3 matrix. The first and second columns
% specify the begining and end samples of a movement period (saccade, joystick, ...),
% and the third column contains the peak velocity/acceleration movement. The thrid
% column allows to convert movements into spike data representation, making it
% compatible with the spike toolbox functions.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_DATABROWSER, FT_DATATYPE_SPIKE

% Copyright (C) 2014, Diego Lozano-Soldevilla
%
% 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$

% FIXME the help mentioned the

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the initial part deals with parsing the input options and data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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

% ensure that the input data is valid for this function, this will also do
% backward-compatibility conversions of old data that for example was
% read from an old *.mat file
data = ft_checkdata(data, 'datatype', {'raw'}, 'feedback', 'yes', 'hassampleinfo', 'yes');

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels', 'trial'}); % prevent accidental typos, see issue 1729

% set the defaults
cfg.method   = ft_getopt(cfg, 'method',   'velocity2D');
cfg.feedback = ft_getopt(cfg, 'feedback', 'yes');

if isfield(data, 'fsample')
  fsample = getsubfield(data, 'fsample');
else
  fsample = 1./(mean(diff(data.time{1})));
end

% set the defaults for the various microsaccade detection methods
switch cfg.method
  case 'velocity2D'
    % Engbert R, Kliegl R (2003) Microsaccades uncover the orientation of
    % covert attention. Vision Res 43:1035-1045.
    kernel = [1 1 0 -1 -1].*(fsample/6); % this is equivalent to Engbert et al (2003) Vis Res, eqn. (1)
    if ~isfield(cfg.velocity2D, 'kernel'),   cfg.velocity2D.kernel  = kernel; end
    if ~isfield(cfg.velocity2D, 'demean'),   cfg.velocity2D.demean  = 'yes';  end
    if ~isfield(cfg.velocity2D, 'mindur'),   cfg.velocity2D.mindur  =  3;     end % minimum microsaccade duration in samples
    if ~isfield(cfg.velocity2D, 'velthres'), cfg.velocity2D.velthres = 6;     end
  case 'clustering'
    ft_error('not implemented yet');
    % Otero-Millan J, Castro JLA, Macknik SL, Martinez-Conde S (2014)
    % Unsupervised clustering method to detect microsaccades. J Vis 14.
  otherwise
    ft_error('unsupported option for cfg.method');
end % switch method

% select channels and trials of interest, by default this will select all channels and trials
tmpcfg = keepfields(cfg, {'trials', 'channel', 'tolerance', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
data = ft_selectdata(tmpcfg, data);
[cfg, data] = rollback_provenance(cfg, data);

% determine the size of the data
ntrial = length(data.trial);
nchan  = length(data.label);   % number of channels

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the actual computation is done in the middle part
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
movement = [];

ft_progress('init', cfg.feedback, 'processing trials');
% do all the computations
for i=1:ntrial
  ft_progress(i/ntrial, 'finding microsaccades trial %d of %d\n', i, ntrial);

  dat = data.trial{i};
  nsample = size(dat,2);

  switch cfg.method
    case 'velocity2D'

      % demean horizontal and vertical time courses
      if strcmp(cfg.velocity2D.demean, 'yes')
        dat = ft_preproc_polyremoval(dat, 0, 1, nsample);
      end

      %% eye velocity computation
      % deal with padding
      n = size(cfg.velocity2D.kernel,2);
      pad = ceil(n/2);
      dat = ft_preproc_padding(dat, 'localmean', pad);

      % convolution. See Engbert et al (2003) Vis Res, eqn. (1)
      if n<100
        % heuristic: for large kernel the convolution is faster when done along
        % the columns, weighing against the costs of doing the transposition.
        % the threshold of 100 is a bit ad hoc.
        vel = convn(dat,   cfg.velocity2D.kernel,   'same');
      else
        vel = convn(dat.', cfg.velocity2D.kernel.', 'same').';
      end
      % cut the eges
      vel = ft_preproc_padding(vel, 'remove', pad);

      %% microsaccade detection
      % compute velocity thresholds as in Engbert et al (2003) Vis Res, eqn. (2)
      medianstd = sqrt( median(vel.^2,2) - (median(vel,2)).^2 );

      % Engbert et al (2003) Vis Res, eqn. (3)
      radius = cfg.velocity2D.velthres*medianstd;

      % compute test criterion: ellipse equation
      test = sum((vel./radius(:,ones(1,nsample))).^2,1);
      sacsmp = find(test>1); % microsaccade's indexing

      %% determine microsaccades per trial
      % first find eye movements of n-consecutive time points
      j = find(diff(sacsmp)==1);
      j1 = [j; j+1];
      com = intersect(j,j+1);
      cut = ~ismember(j1,com);
      sacidx = reshape(j1(cut),2,[]);

      for k=1:size(sacidx,2)
        duration = sacidx(1,k):sacidx(2,k);
        if size(duration,2) >= cfg.velocity2D.mindur
          % finding peak velocity by Pitagoras
          begtrl = sacsmp(duration(1,1));
          endtrl = sacsmp(duration(1,end));

          [peakvel, smptrl] = max(sqrt(sum(vel(:,begtrl:endtrl).^2,1)));
          veltrl = sacsmp(duration(1,smptrl)); % peak velocity microsaccade sample -> important for spike conversion

          trlsmp = data.sampleinfo(i,1):data.sampleinfo(i,2);
          begsample = trlsmp(1, begtrl); % begining microsaccade sample
          endsample = trlsmp(1, endtrl); % end microsaccade sample
          velsample = trlsmp(1, veltrl); % velocity peak microsaccade sample
          movement(end+1,:) = [begsample endsample velsample];
        end
      end

    case 'clustering'
      ft_error('not implemented yet');
      % Otero-Millan J, Castro JLA, Macknik SL, Martinez-Conde S (2014)
      % Unsupervised clustering method to detect microsaccades. J Vis 14.

    otherwise
      ft_error('unsupported option for cfg.method');
  end % switch method
end % for each trial
ft_progress('close');

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


text:
function ellipsLat=reducedLat2EllipsLat(reducedLat,f)
%%REDUCEDLAT2ELLIPSLAT Convert reduced (parameteric) latitude into an
%                      ellipsoidal (geodetic) latitude. A parameteric
%                      latitude is an angle of elevation drawn from the
%                      center of a sphere whose radius is equal to the
%                      semi-major axis of the ellipsoid. The parametric
%                      latitude is the angle of the point on the sphere
%                      where a vertical line drawn from the major axis of
%                      the ellipse through the point on the ellipse
%                      intersects the sphere. This is a way of mapping an
%                      ellipsoid onto a sphere.
%
%INPUTS: reducedLat A vector or matrix of reduced latitudes in radians.
%                   These typically range from -pi/2 to pi/2.
%                 f The flattening factor of the reference ellipsoid. If
%                   this argument is omitted, the value in
%                   Constants.WGS84Flattening is used.
%
%OUTPUTS: reducedLat The ellipsoidal (geodetic) latitude ranging from
%                   -pi/2 to pi/2, for each reduced latitude.
%
%%The definition of the reduced (parameteric) latitude from Chapter 3 (pg.
%18) of [1] is used. The conversion relation between the latitudes is in
%Section 3.4 of [2].
%
%The reduced latitude has been used for solving geodescic problems on an
%ellipsoid by transforming the problem to an equivalent problem on a
%sphere.
%
%The inverse of this function is ellipsLat2ReducedLat.
%
%REFERENCES:
%[1] J. P. Snyder, "Map projections- a working manual," U.S. Geological
%    Survey, Tech. Rep. 1395, 1987.
%[2] R. H. Rapp, "Geometric geodesy, part I," Ohio State University
%    Department of Geodetic Science and Surveying, Tech. Rep., Apr. 1991.
%    [Online]. Available: http://hdl.handle.net/1811/24333
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2)
    f=Constants.WGS84Flattening;
end

%The first numerical eccentricity of the ellipsoid.
e=sqrt(2*f-f^2);

%The wrapping helps deal with precision limitations that might push the
%input value slightly above or below +/-pi/2. If the wrapping were not
%done, then instead of getting, for example, pi/2, one will get pi/2 as a
%return value.
reducedLat=wrapRange(reducedLat,-pi/2,pi/2,true);

ellipsLat=atan(tan(reducedLat)/sqrt(1-e^2));
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function varargout = pdftops(cmd)
%PDFTOPS  Calls a local pdftops executable with the input command
%
% Example:
%   [status result] = pdftops(cmd)
%
% Attempts to locate a pdftops executable, finally asking the user to
% specify the directory pdftops was installed into. The resulting path is
% stored for future reference.
%
% Once found, the executable is called with the input command string.
%
% This function requires that you have pdftops (from the Xpdf package)
% installed on your system. You can download this from:
% http://www.foolabs.com/xpdf
%
% IN:
%   cmd - Command string to be passed into pdftops.
%
% OUT:
%   status - 0 iff command ran without problem.
%   result - Output from pdftops.

% Copyright: Oliver Woodford, 2009-2010

% Thanks to Jonas Dorn for the fix for the title of the uigetdir window on
% Mac OS.
% Thanks to Christoph Hertel for pointing out a bug in check_xpdf_path
% under linux.
% 23/01/2014 - Add full path to pdftops.txt in warning.

% Call pdftops
[varargout{1:nargout}] = system(sprintf('"%s" %s', xpdf_path, cmd));
return

function path_ = xpdf_path
% Return a valid path
% Start with the currently set path
path_ = user_string('pdftops');
% Check the path works
if check_xpdf_path(path_)
    return
end
% Check whether the binary is on the path
if ispc
    bin = 'pdftops.exe';
else
    bin = 'pdftops';
end
if check_store_xpdf_path(bin)
    path_ = bin;
    return
end
% Search the obvious places
if ispc
    path_ = 'C:\Program Files\xpdf\pdftops.exe';
else
    path_ = '/usr/local/bin/pdftops';
end
if check_store_xpdf_path(path_)
    return
end
% Ask the user to enter the path
while 1
    if strncmp(computer,'MAC',3) % Is a Mac
        % Give separate warning as the uigetdir dialogue box doesn't have a
        % title
        uiwait(warndlg('Pdftops not found. Please locate the program, or install xpdf-tools from http://users.phg-online.de/tk/MOSXS/.'))
    end
    base = uigetdir('/', 'Pdftops not found. Please locate the program.');
    if isequal(base, 0)
        % User hit cancel or closed window
        break;
    end
    base = [base filesep];
    bin_dir = {'', ['bin' filesep], ['lib' filesep]};
    for a = 1:numel(bin_dir)
        path_ = [base bin_dir{a} bin];
        if exist(path_, 'file') == 2
            break;
        end
    end
    if check_store_xpdf_path(path_)
        return
    end
end
error('pdftops executable not found.');

function good = check_store_xpdf_path(path_)
% Check the path is valid
good = check_xpdf_path(path_);
if ~good
    return
end
% Update the current default path to the path found
if ~user_string('pdftops', path_)
    warning('Path to pdftops executable could not be saved. Enter it manually in %s.', fullfile(fileparts(which('user_string.m')), '.ignore', 'pdftops.txt'));
    return
end
return

function good = check_xpdf_path(path_)
% Check the path is valid
[good, message] = system(sprintf('"%s" -h', path_));
% system returns good = 1 even when the command runs
% Look for something distinct in the help text
good = ~isempty(strfind(message, 'PostScript'));
return

author: Cloud-CV

repo: object-proposals