text:
% VL_DEMO_SIFT_BASIC  Demo: SIFT: basic functionality

pfx = fullfile(vl_root,'figures','demo') ;
randn('state',0) ;
rand('state',0) ;
figure(1) ; clf ;

% --------------------------------------------------------------------
%                                                        Load a figure
% --------------------------------------------------------------------
I = imread(fullfile(vl_root,'data','roofs1.jpg')) ;

image(I) ; colormap gray ;
axis equal ; axis off ; axis tight ;
vl_demo_print('sift_basic_0') ;

% --------------------------------------------------------------------
%                                       Convert the to required format
% --------------------------------------------------------------------
I = single(rgb2gray(I)) ;

clf ; imagesc(I)
axis equal ; axis off ; axis tight ;
vl_demo_print('sift_basic_1') ;

% --------------------------------------------------------------------
%                                                             Run SIFT
% --------------------------------------------------------------------
[f,d] = vl_sift(I) ;

hold on ;
perm = randperm(size(f,2)) ;
sel  = perm(1:50) ;
h1   = vl_plotframe(f(:,sel)) ; set(h1,'color','k','linewidth',3) ;
h2   = vl_plotframe(f(:,sel)) ; set(h2,'color','y','linewidth',2) ;

vl_demo_print('sift_basic_2') ;

delete([h1 h2]);

h3 = vl_plotsiftdescriptor(d(:,sel),f(:,sel)) ;
set(h3,'color','k','linewidth',2) ;
h4 = vl_plotsiftdescriptor(d(:,sel),f(:,sel)) ;
set(h4,'color','g','linewidth',1) ;
h1   = vl_plotframe(f(:,sel)) ; set(h1,'color','k','linewidth',3) ;
h2   = vl_plotframe(f(:,sel)) ; set(h2,'color','y','linewidth',2) ;

vl_demo_print('sift_basic_3') ;

% --------------------------------------------------------------------
%                                                      Custom keypoint
% --------------------------------------------------------------------
clf ; imagesc(I(1:200,1:200,:))
axis equal ; axis off ; axis tight ;
hold on;

fc = [100;100;10;-pi/8] ;
[f,d] = vl_sift(I,'frames',fc) ;

h3   = vl_plotsiftdescriptor(d,f) ;  set(h3,'color','k','linewidth',3) ;
h4   = vl_plotsiftdescriptor(d,f) ;  set(h4,'color','g','linewidth',2) ;
h1   = vl_plotframe(f) ; set(h1,'color','k','linewidth',4) ;
h2   = vl_plotframe(f) ; set(h2,'color','y','linewidth',2) ;

vl_demo_print('sift_basic_4') ;

% --------------------------------------------------------------------
%                                   Custom keypoints with orientations
% --------------------------------------------------------------------

fc = [100;100;10;0] ;
[f,d] = vl_sift(I,'frames',fc,'orientations') ;

delete([h1 h2 h3 h4]) ;

h3   = vl_plotsiftdescriptor(d,f) ;  set(h3,'color','k', 'linewidth',3) ;
h4   = vl_plotsiftdescriptor(d,f) ;  set(h4,'color','g', 'linewidth',2) ;
h1   = vl_plotframe(f) ; set(h1,'color','k','linewidth',4) ;
h2   = vl_plotframe(f) ; set(h2,'color','y','linewidth',2) ;

vl_demo_print('sift_basic_5') ;

author: Cloud-CV

repo: object-proposals


text:
% Copyright Mohammad SAFEEA, 17th-Aug-2017
function setJointAngleToText(anglesTextHandlesCellArray, index,angle )
%% Set the joint angle in degrees for the
    angle=180*angle/pi;
    set(anglesTextHandlesCellArray{index},'String',num2str(angle));
    %set(anglesTextHandlesCellArray{1},'BackgroundColor','red');

end

repo: KST-Kuka-Sunrise-Toolbox


text:
function [solvers,keep,allsolvers] = getavailablesolvers(findallsolvers,options);

solvers = definesolvers;allsolvers = solvers;
keep = ones(length(solvers),1);

if ~findallsolvers
    for i = 1:length(solvers)
        isavailable = 1;
        j = 1;

        while (j <= length(solvers(i).checkfor)) && isavailable
            s = exist(solvers(i).checkfor{j},'file');
            s = (s~=0) & (s~=7);
            isavailable = isavailable & s;
            j = j + 1;
        end
        if ~isavailable
            keep(i)=0;
        end
    end
end

if nargout == 1 || nargout==3
    solvers = solvers(find(keep));
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function L = boundary_length(V,F)
%BOUNDARY_LENGTH Compute the length of a mesh's boundary
%
% Input:
%  V,F  input mesh
% Output:
%  L  the length of the boundary

%Compute the indices of all boundary edges.
O = outline(F);

%Compute a matrix whose rows are equal to edgeEnd-edgeStart for each
%boundary edge.
BE = V(O(:,2),:) - V(O(:,1),:);

%Compute the total length of all boundary edges.
L = sum(normrow(BE));

end

author: odedstein

repo: sgi-introduction-course


text:
function EVAL = Evaluate(ACTUAL,PREDICTED)
% This fucntion evaluates the performance of a classification model by
% calculating the common performance measures: Accuracy, Sensitivity,
% Specificity, Precision, Recall, F-Measure, G-mean.
% Input: ACTUAL = Column matrix with actual class labels of the training
%                 examples
%        PREDICTED = Column matrix with predicted class labels by the
%                    classification model
% Output: EVAL = Row matrix with all the performance measures

idx = (ACTUAL()==1);

p = length(ACTUAL(idx));
n = length(ACTUAL(~idx));
N = p+n;

tp = sum(ACTUAL(idx)==PREDICTED(idx));
tn = sum(ACTUAL(~idx)==PREDICTED(~idx));
fp = n-tn;
fn = p-tp;

tp_rate = tp/p;
tn_rate = tn/n;

accuracy = (tp+tn)/N;
sensitivity = tp_rate;
specificity = tn_rate;
precision = tp/(tp+fp);
recall = sensitivity;
f_measure = 2*((precision*recall)/(precision + recall));
gmean = sqrt(tp_rate*tn_rate);

EVAL = [accuracy sensitivity specificity precision recall f_measure gmean];

author: jindongwang

repo: activityrecognition


text:
% test for drawing a spherical triangulation

rep = '../toolbox_graph_data/wrl/';
name = 'pawn';
name = 'test2';
name = 'test1';
[vertex,face] = read_wrl([rep name '.wrl']);

% project on sphere [should use here a spherical parameterization instead]
d = sqrt(sum(vertex.^2));
svertex = vertex./repmat(d, [3 1]);

% draw simply the projected mesh
clf;
plot_mesh(svertex,face);
shading faceted;

% now draw nice curve on the sphere
options.target_face = round(size(svertex,2)/2);
clf;
plot_spherical_triangulation(svertex,face, options);

author: gpeyre

repo: matlab-toolboxes


text:
% TEST_HANNING_STANDALONE Test case for HANNING_STANDALONE
%
% See also
%   HANNING (Signal Processing Toolbox)

classdef test_hanning_standalone < matlab.unittest.TestCase
    methods (Test)
        function testBasic(testcase)
            maxLength = 8191;
            for window_length = 1:maxLength % lasts a few seconds
                hanning_old = hanning(window_length); % requires DSP toolbox
                hanning_new = hanning_standalone(window_length);
                testcase.assertEqual(all(hanning_old==hanning_new));
            end
        end
    end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [] = test_shrink_mex()
% Test shrink_mex.c and shrink_mex2.cc

    check_correctness()
    time_single_run()

end

function [] = check_correctness()
    rng(271828);

    m = 700;
    n = 600;
    X = randn(m,n);
    lambda = rand();
    %NOTE: shrink_mex doesn't properly handle vector lambda; shrink_mex2 fixes this
    %lambda = rand(1,n); % shrink_mex doesn't properly handle vector lambda; shrink_mex2 fixes this
    soffset = randn();
    voffset = randn(1,n);

    shrink_ref = @(X, lambda) sign(X) .* max( abs(X) - lambda, 0);
    shrink_test = @(X, lambda) shrink_mex(X, lambda);
    shrink_mex2(struct('num_threads', 1));
    shrink_test2 = @(X, lambda) shrink_mex2(X, lambda);

    fprintf('Testing shrink(X, lambda):\n');
    ref = shrink_ref(X, lambda);
    test = shrink_test(X, lambda);
    test2 = shrink_test2(X, lambda);
    fprintf('  rel_err  = %1.5e\n', norm(test - ref, 'fro') / norm(ref, 'fro'));
    fprintf('  rel_err2 = %1.5e\n', norm(test2 - ref, 'fro') / norm(ref, 'fro'));

    shrink_ref = @(X, lambda, offset) sign(X - offset) .* max( abs(X - offset) - lambda, 0);
    shrink_test = @(X, lambda, offset) shrink_mex(X, lambda, offset);
    shrink_mex2(struct('num_threads', 1));
    shrink_test2 = @(X, lambda, offset) shrink_mex2(X, lambda, offset);

    fprintf('Testing shrink(X, lambda, scalar offset):\n');
    ref = shrink_ref(X, lambda, soffset);
    test = shrink_test(X, lambda, soffset);
    test2 = shrink_test2(X, lambda, soffset);
    fprintf('  rel_err  = %1.5e\n', norm(test - ref, 'fro') / norm(ref, 'fro'));
    fprintf('  rel_err2 = %1.5e\n', norm(test2 - ref, 'fro') / norm(ref, 'fro'));

    fprintf('Testing shrink(X, lambda, vector offset):\n');
    ref = shrink_ref(X, lambda, voffset);
    test = shrink_test(X, lambda, voffset);
    test2 = shrink_test2(X, lambda, voffset);
    fprintf('  rel_err  = %1.5e\n', norm(test - ref, 'fro') / norm(ref, 'fro'));
    fprintf('  rel_err2 = %1.5e\n', norm(test2 - ref, 'fro') / norm(ref, 'fro'));

end

function [] = time_single_run()
% Time shrinkage implementations

    minimum_runtime = 1;

    % An example from sparse subspace clustering (SSC)
    n_rows = 6000;
    n_cols = 6000;
    x = randn(n_rows, n_cols);
    lambda = 1;
    offset = randn(1, n_cols);

    opt = struct('num_threads', 4);
    shrink_mex2(opt);

    % Make the shrinkage operator
    %shrink = @(X, lambda) sign(X) .* max( abs(X) - lambda, 0);
    %shrink = @(X, lambda) shrink_mex(X, lambda);
    shrink = @(X, lambda) shrink_mex2(X, lambda);

    %shrink = @(X, lambda, offset) sign(X - offset) .* max( abs(X - offset) - lambda, 0);
    %shrink = @(X, lambda, offset) shrink_mex(X, lambda, offset);
    %shrink = @(X, lambda, offset) shrink_mex2(X, lambda, offset);

    % Warm up
    n_done = 0;
    t_ = tic();
    while true
        if toc(t_) >= minimum_runtime
            break
        end

        y = shrink(x, lambda);
        %y = shrink(x, lambda, offset);
        n_done = n_done + 1;
    end

    % Measure runtime
    times = zeros(n_done,1);
    for n=1:n_done
        t_ = tic();
        y = shrink(x, lambda);
        %y = shrink(x, lambda, offset);
        times(n) = toc(t_);
    end

    fprintf('shrinkage min/mean/max runtime = %1.5e  %1.5e  %1.5e  seconds\n', min(times), mean(times), max(times));

end

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function f = tan( f )
%TAN   Tangent of a SEPARABLEAPPROX.

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

% Empty check:
if ( isempty( f ) )
    return
end

f = compose( f, @tan );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function i4_log_10_test ( )

%*****************************************************************************80
%
%% I4_LOG_10_TEST tests I4_LOG_10.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 13;

  x = [ 0, 1, 2, 3, 9, 10, 11, 99, 101, -1, -2, -3, -9 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_LOG_10_TEST\n' );
  fprintf ( 1, '  I4_LOG_10: whole part of log base 10,\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  X, I4_LOG_10\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n

    fprintf ( 1, '%6d  %12d\n', x(i), i4_log_10 ( x(i) ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef TestRotate
    %TestRotate

    methods (Static)
        function test_1
            src = magic(5);

            codes = {'90CW', '90CCW', '180'};
            for i=1:numel(codes)
                dst = cv.rotate(src, codes{i});
                validateattributes(dst, {class(src)}, {'size',size(src)});
            end
        end

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

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function monogrid_poisson_1d_test01_mono ( )

%*****************************************************************************80
%
%% MONOGRID_POISSON_1D_TEST01_MONO tests MONOGRID_POISSON_1D on test case 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 December 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MONOGRID_POISSON_1D_TEST01_MONO\n' );
  fprintf ( 1, '  MONOGRID_POISSON_1D solves a 1D Poisson BVP\n' );
  fprintf ( 1, '  using the Gauss-Seidel method.\n' );

  a = 0.0;
  b = 1.0;
  ua = 0.0;
  ub = 0.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  -u"(x) = 1, for %g < x < %g\n', a, b );
  fprintf ( 1, '  u(%g) = %g, u(%g) = %g.\n', a, ua, b, ub );
  fprintf ( 1, '  Solution is u(x) = ( -x^2 + x ) / 2\n' );

  for k = 5 : 5

    n = 2^k;

    x = ( linspace ( a, b, n + 1 ) )';

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Mesh index K = %d\n', k );
    fprintf ( 1, '  Number of intervals N=2^K = %d\n', n );
    fprintf ( 1, '  Number of nodes = 2^K+1 =   %d\n', n + 1 );

    [ u, it_num ] = monogrid_poisson_1d ( n, a, b, ua, ub, @force1, @exact1 );

    fprintf ( 1, '\n' );
    fprintf ( 1, '     I        X(I)      U(I)         U Exact(X(I))\n' );
    fprintf ( 1, '\n' );
    for i = 1 : n + 1
      fprintf ( 1, '  %4d  %10f  %14g  %14g\n', i, x(i), u(i), exact1 ( x(i) ) );
    end

    fprintf ( 1, '\n' );

    difmax = 0.0;
    for i = 1 : n + 1
      difmax = max ( difmax, abs ( u(i) - exact1 ( x(i) ) ) );
    end
    fprintf ( 1, '  Maximum error = %g\n', difmax );
    fprintf ( 1, '  Number of iterations = %d\n', it_num );

  end

  return
end


text:
function [hPatch,outP,FV, cl, myLight] = mask2surface(P,varargin)
% ::
%
%    [hPatch,outP,FV, cl, myLight] = mask2surface(P)
%
% :Input:
%
%   **P:**
%        a file name for a mask file, OR a clusters structure
%
% :Inputs:
%
%   **hPatch:**
%        handle to surface patch
%
%   **outP:**
%        file name, with path, of mat file containing faces and vertices
%
%   **FV:**
%        isosurface
%
%   **cl:**
%        clusters structure with coords and meshgrid info
%
%   **varargin:**
%        suppress lighting (0) or do lighting (1)
%
%   **varargin{2}:**
%        color for patch objects
%
% ..
%    tor wager 9/23/02
% ..
%
% uses get_cluster_volume, adapted from Sergey Pakhomov
%
% use with cluster_surf.m to map activations onto surfaces:
% cluster_surf(clusters1,clusters2,outP,10,{[0 1 0] [1 0 0]})
%
% or use getvertexcolors.m to map colors onto surface using hPatch.

hPatch = []; outP = [];  cl = []; myLight = [];
mycolor = [.5 .5 .5];
myalpha = 1;

P = which(P);
if isempty(P), error('Cannot find image.');, end

if isstr(P), cl = mask2clusters(P);, else, cl = P;, P=cl(1).P;, end

height_threshold = .9;
dolight = 1;

if length(varargin) > 0, dolight = varargin{1};, end
if length(varargin) > 1, mycolor = varargin{2};, end

set(0,'CurrentFigure', gcf)

omit = cat(1,cl.numVox);
cl(find(omit==1)) = [];

for i = 1:length(cl)

    % added for SPM2 compatibility
    cl(i).voxSize = abs(cl(i).voxSize);

    cl2(i) = get_cluster_volume(cl(i));
    V = smooth3(cl2(i).vTal,'gaussian',[5 5 5]);

	FV(i) = isosurface(cl2(i).xTal, cl2(i).yTal, cl2(i).zTal, V, height_threshold);
	hPatch(i) = patch(FV);
	isonormals(cl2(i).xTal, cl2(i).yTal, cl2(i).zTal, cl2(i).vTal, hPatch(i));

    try
		set(hPatch, 'Tag', P, 'FaceColor', mycolor, 'EdgeColor', 'none', 'FaceAlpha', myalpha);
	catch
		set(hPatch, 'Tag', P, 'FaceColor', mycolor, 'EdgeColor', 'none');
	end

end

cl = cl2;  clear cl2

if dolight
    view(135,30);
    myLight = standardMRIlighting('full',[hPatch hPatch]);
    material dull
end

[d f e] = fileparts(P);
outP = fullfile(d, filesep,['surf_' f '.mat']);
%if outP(1) == filesep, outP = outP(2:end);end
save(outP, 'FV');

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% Maximum volume inscribed ellipsoid in a polyhedron
% Section 8.4.1, Boyd & Vandenberghe "Convex Optimization"
% Original version by Lieven Vandenberghe
% Updated for CVX by Almir Mutapcic - Jan 2006
% (a figure is generated)
%
% We find the ellipsoid E of maximum volume that lies inside of
% a polyhedra C described by a set of linear inequalities.
%
% C = { x | a_i^T x <= b_i, i = 1,...,m } (polyhedra)
% E = { Bu + d | || u || <= 1 } (ellipsoid)
%
% This problem can be formulated as a log det maximization
% which can then be computed using the det_rootn function, ie,
%     maximize     log det B
%     subject to   || B a_i || + a_i^T d <= b,  for i = 1,...,m

% problem data
n = 2;
px = [0 .5 2 3 1];
py = [0 1 1.5 .5 -.5];
m = size(px,2);
pxint = sum(px)/m; pyint = sum(py)/m;
px = [px px(1)];
py = [py py(1)];

% generate A,b
A = zeros(m,n); b = zeros(m,1);
for i=1:m
  A(i,:) = null([px(i+1)-px(i) py(i+1)-py(i)])';
  b(i) = A(i,:)*.5*[px(i+1)+px(i); py(i+1)+py(i)];
  if A(i,:)*[pxint; pyint]-b(i)>0
    A(i,:) = -A(i,:);
    b(i) = -b(i);
  end
end

% formulate and solve the problem
cvx_begin
    variable B(n,n) symmetric
    variable d(n)
    maximize( det_rootn( B ) )
    subject to
       for i = 1:m
           norm( B*A(i,:)', 2 ) + A(i,:)*d <= b(i);
       end
cvx_end

% make the plots
noangles = 200;
angles   = linspace( 0, 2 * pi, noangles );
ellipse_inner  = B * [ cos(angles) ; sin(angles) ] + d * ones( 1, noangles );
ellipse_outer  = 2*B * [ cos(angles) ; sin(angles) ] + d * ones( 1, noangles );

clf
plot(px,py)
hold on
plot( ellipse_inner(1,:), ellipse_inner(2,:), 'r--' );
plot( ellipse_outer(1,:), ellipse_outer(2,:), 'r--' );
axis square
axis off
hold off

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [out] = saturation_9(In,S,St,varargin)
%saturation_9

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Deficit store: Saturation excess from a store that has
%               reached maximum capacity
% Constraints:  -
% @(Inputs):    In   - incoming flux [mm/d]
%               S    - current storage [mm]
%               St   - threshold for flow generation [mm], 0 for deficit
%               store
%               varargin(1) - smoothing variable r (default 0.01)
%               varargin(2) - smoothing variable e (default 5.00)

if size(varargin,2) == 0
    out = In.*smoothThreshold_storage_logistic(S,St);
elseif size(varargin,2) == 1
    out = In.*smoothThreshold_storage_logistic(S,St,varargin(1));
elseif size(varargin,2) == 2
    out = In.*smoothThreshold_storage_logistic(S,St,varargin(1),varargin(2));
end

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function [ dose ] = matRad_interpDicomDoseCube( ct, currDose )
% matRad function to interpolate a given Dicom Dose Cube dicom RTDOSE data
%
% call
%   [ dose ] = matRad_interpDicomDoseCube( ct, currDose )
%
% input
%   ct:             ct imported by the matRad_importDicomCt function
%   currDose:       one (of several) dose cubes which should be interpolated
%
% output
%   dose:           struct with different actual current dose cube and several
%                   meta data
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% read information out of the RT file
dosefile = currDose{1};
if verLessThan('matlab','9')
    doseInfo = dicominfo(dosefile);
else
    doseInfo = dicominfo(dosefile,'UseDictionaryVR',true);
end

% read the dosefile itself
dosedata = dicomread(dosefile);
dose.cube = double(dosedata);
% dose.cube = single(dosedata);

% give it an internal name
dose.internalName = currDose{12};

% read out the resolution
dose.resolution.x = doseInfo.PixelSpacing(1);
dose.resolution.y = doseInfo.PixelSpacing(2);
dose.resolution.z = doseInfo.SliceThickness;

% target resolution is ct.resolution
target_resolution = ct.resolution;

% convert dosedata to 3-D cube
dose.cube = squeeze(dose.cube(:,:,1,:));

% ct resolution is target resolution, now convert to new cube;

% generating grid vectors
x = doseInfo.ImagePositionPatient(1) + doseInfo.ImageOrientationPatient(1) * ...
                                       doseInfo.PixelSpacing(1) * double([0:doseInfo.Columns - 1]);
y = doseInfo.ImagePositionPatient(2) + doseInfo.ImageOrientationPatient(5) * ...
                                       doseInfo.PixelSpacing(2) * double([0:doseInfo.Rows - 1]);
z = [doseInfo.ImagePositionPatient(3) + doseInfo.GridFrameOffsetVector];

% set up grid matrices - implicit dimension permuation (X Y Z-> Y X Z)
% Matlab represents internally in the first matrix dimension the
% ordinate axis and in the second matrix dimension the abscissas axis
[ X,  Y,  Z] = meshgrid(x,y,z);
[Xq, Yq, Zq] = meshgrid(ct.x,ct.y,ct.z);

% get GridScalingFactor
gridScale = double(doseInfo.DoseGridScaling);
% rescale dose.cube
dose.cube = gridScale * dose.cube;

% interpolation to ct grid - cube is now stored in Y X Z
dose.cube = interp3(X,Y,Z,dose.cube,Xq,Yq,Zq,'linear',0);

% write new parameters
dose.resolution = ct.resolution;
dose.x = ct.x;
dose.y = ct.y;
dose.z = ct.z;

% write Dicom-Tags
dose.dicomInfo.PixelSpacing            = [target_resolution.x; ...
                                                target_resolution.y];
dose.dicomInfo.ImagePositionPatient    = [min(dose.x); min(dose.y); min(dose.z)];
dose.dicomInfo.SliceThickness          = target_resolution.z;
dose.dicomInfo.ImageOrientationPatient = doseInfo.ImageOrientationPatient;
dose.dicomInfo.DoseType                = doseInfo.DoseType;
dose.dicomInfo.DoseSummationType       = doseInfo.DoseSummationType;
%dose.dicomInfo.InstanceNumber          = doseInfo.InstanceNumber; %Not
%always given
dose.dicomInfo.SOPClassUID             = doseInfo.SOPClassUID;
dose.dicomInfo.SOPInstanceUID          = doseInfo.SOPInstanceUID;
dose.dicomInfo.ReferencedRTPlanSequence = doseInfo.ReferencedRTPlanSequence;

end

repo: matRad


text:
function asa266_test07 ( )

%*****************************************************************************80
%
%% TEST07 tests PPCHI2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 August 2010
%
%  Author:
%
%    John Burkardt
%
  nitest = 9;
  njtest = 9;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST07\n' );
  fprintf ( 1, '  PPCHI2 computes the percentage points\n' );
  fprintf ( 1, '  of the chi squared distribution.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  CDF, PPCHI2(CDF)\n' );
  fprintf ( 1, '\n' );

  for j = 1 : njtest

    v = j;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  For Chi^2 parameter value %f\n', v );
    fprintf ( 1, '\n' );

    for i = 1 : nitest

      cdf = i / ( nitest + 1 );
      [ x1, ifault ] = ppchi2 ( cdf, v );

      fprintf ( 1, '  %12f  %12f\n', cdf, x1 );

    end

  end

  return
end


text:
function status = isrealvec(x)

% ISREALVEC returns true for a real row or column vector
%
% Use as
%   status = isrealvec(x)
%
% See also ISNUMERIC, ISREAL, ISVECTOR, ISREALMAT

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

% This is a drop-in replacement for MATLAB/toolbox/ident/idutils/isrealvec.m
% which requires the Ident toolbox. The following can be used instead to
% replicate the same functionality.
%
% built-in (/Applications/MATLAB_R2010a.app/toolbox/matlab/datatypes/isnumeric)
% built-in (/Applications/MATLAB_R2010a.app/toolbox/matlab/elfun/isreal)
% built-in (/Applications/MATLAB_R2010a.app/toolbox/matlab/elmat/isvector)
% /Applications/MATLAB_R2010a.app/toolbox/matlab/elmat/ndims.m

status = isnumeric(x) && isreal(x) && length(size(x))==2 && any(size(x)==1);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = process_ersd_ab( varargin )
% PROCESS_ERSD_AB: Compute event related perturbation (synchrnonization / desynchrnonization)
%
% DESCRIPTION:
%    This function calculates event related perturbation (ERS/ERD) as the percentage of a
%    decrease or increase during a test interval (T), as compared to a reference interval (R).
%    The following formula is used: ERSP = (R-T)/R x 100.

% @=============================================================================
% 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: Nikola Vukovic, University of Cambridge, 2013
%          Francois Tadel, 2013-2015

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Event related perturbation (ERS/ERD) [DEPRECATED]';
    sProcess.FileTag     = 'ersd';
    sProcess.Category    = 'Filter2';
    sProcess.SubGroup    = 'Standardize';
    sProcess.Index       = 205;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/TimeFrequency#Normalized_time-frequency_maps';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'results', 'timefreq', 'matrix'};
    sProcess.OutputTypes = {'data', 'results', 'timefreq', 'matrix'};
    sProcess.nInputs     = 2;
    sProcess.nMinFiles   = 1;
    % Default values for some options
    sProcess.isSourceAbsolute = 1;
    sProcess.processDim       = 1;    % Process channel by channel
    sProcess.isPaired         = 1;

    % Definition of the options
    sProcess.options.description.Comment = ['For each signal in input:
' ...
                                            '1) FilesA: Compute the mean m for the baseline
' ...
                                            '2) FilesB: For each time sample, calculates a percentage of increase/decrease
' ...
                                            'ERSD = (Data-m)/m x 100

'];
    sProcess.options.description.Type    = 'label';
    % === Baseline time window
    sProcess.options.baseline.Comment = 'Baseline (Files A):';
    sProcess.options.baseline.Type    = 'baseline';
    sProcess.options.baseline.Value   = [];
    % === Sensor types
    sProcess.options.sensortypes.Comment = 'Sensor types or names (empty=all): ';
    sProcess.options.sensortypes.Type    = 'text';
    sProcess.options.sensortypes.Value   = 'MEG, EEG';
    sProcess.options.sensortypes.InputTypes = {'data'};
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    % Get baseline
    if isfield(sProcess.options, 'baseline') && isfield(sProcess.options.baseline, 'Value') && iscell(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value{1})
        Time = sProcess.options.baseline.Value{1};
    else
        Time = [];
    end
    % Format comment
    if isempty(Time)
        Comment = 'Event related perturbation: [All file]';
    elseif any(abs(Time) > 2)
        Comment = sprintf('Event related perturbation: [%1.3fs,%1.3fs]', Time(1), Time(2));
    else
        Comment = sprintf('Event related perturbation: [%dms,%dms]', round(Time(1)*1000), round(Time(2)*1000));
    end
end

%% ===== RUN =====
function sInputB = Run(sProcess, sInputA, sInputB) %#ok<DEFNU>
    % Get options
    if isfield(sProcess.options, 'baseline') && isfield(sProcess.options.baseline, 'Value') && iscell(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value{1})
        BaselineBounds = sProcess.options.baseline.Value{1};
    else
        BaselineBounds = [];
    end
    % Get baseline indices
    if ~isempty(BaselineBounds)
        iBaseline = bst_closest(BaselineBounds, sInputA.TimeVector);
        if (iBaseline(1) == iBaseline(2)) && any(iBaseline(1) == sInputA.TimeVector)
            error('Invalid baseline definition.');
        end
        iBaseline = iBaseline(1):iBaseline(2);
    % Get all file
    else
        iBaseline = 1:size(sInputA.A,2);
    end
    % Compute ERS/ERD
    sInputB.A = Compute(sInputA.A(:,iBaseline,:), sInputB.A);
    % Change DataType
    if ~strcmpi(sInputB.FileType, 'timefreq')
        sInputB.DataType = 'zscore';
    end
    % Default colormap
    if strcmpi(sInputB.FileType, 'results')
        sInputB.ColormapType = 'stat1';
    else
        sInputB.ColormapType = 'stat2';
    end
    % Add new tag to comment
    if isfield(sProcess.options, 'source_abs') && sProcess.options.source_abs.Value
        sInputB.Comment = [sInputB.Comment, ' | abs'];
    end
end

%% ===== COMPUTE =====
function B_data = Compute(A_baseline, B_data)
    disp('BST> process_ersd_ab.m is deprecated, use "Standardize > Baseline normalization" instead.');
    % Compute baseline statistics
    meanBaseline = mean(A_baseline, 2);
    % Remove null variance values
    meanBaseline(meanBaseline == 0) = 1e-12;
    % Compute event related perturbation
    B_data = bst_bsxfun(@minus, B_data, meanBaseline);
    B_data = bst_bsxfun(@rdivide, B_data, meanBaseline) .* 100;
end

author: brainstorm-tools


text:
function str = VBA_summary(out,newlines)
% writes a summary string from standard output of VBA model inversion
% function str = VBA_summary(out)
% IN:
%   - out: the 'out' structure of VBA inversion routine
%   - newlines: flag for inserting line separators (\n)
% OUT:
%   - str: a cell array of strings, which summarize the VBA inversion

try;newlines;catch,newlines=0;end

[LLH0] = VBA_LMEH0(out.y,out.options);
try F = out.F(end); catch, F = '?'; end
many = length(out.options.sources)>1; % more than one source?

str{1} = sprintf(['Date: ',datestr(out.date)]);
if ~out.options.OnLine
    s0 = ['VB converged in ',num2str(out.it),' iterations'];
else
    s0 = ['Online VB algorithm'];
end
try
    if floor(out.dt./60) == 0
        timeString = [num2str(floor(out.dt)),' sec'];
    else
        timeString = [num2str(floor(out.dt./60)),' min'];
    end
    str{2} = sprintf([s0,' (took ~',timeString,')']);
catch
    str{2} = sprintf(s0);
end
if many
    datastr = [' (number of sources=',num2str(length(out.options.sources)),')'];
else
    datastr = [];
end
str{3} = sprintf(['Dimensions of the model:','\n ',...
    '    - data: p=',num2str(out.dim.p),datastr,'\n ',...
    '    - time samples: t=',num2str(out.dim.n_t),'\n ',...
    '    - hidden states: n=',num2str(out.dim.n),'\n ',...
    '    - evolution parameters: n_theta=',num2str(out.dim.n_theta),'\n ',...
    '    - observation parameters: n_phi=',num2str(out.dim.n_phi),'\n ',...
    '    - inputs: n_u=',num2str(out.dim.u)]);
if numel(out.options.sources)>1
    tmp = ' (multisource)';
else
    switch out.options.sources.type
        case 0
            tmp = ' (gaussian data)';
        case 1
            tmp = ' (binomial data)';
        case 2
            tmp = ' (multinomial data)';
    end
end
if out.options.UNL
    so = 'un-normalized likelihood';
else
    so = 'observation';
end
if out.dim.n >= 1
    if isinf(out.options.priors.a_alpha) && isequal(out.options.priors.b_alpha,0)
        str{4} = sprintf('This was a deterministic dynamical system');
    else
        str{4} = sprintf('This was a stochastic dynamical system');
    end
    if isa(out.options.g_fname,'function_handle')
        gfn = func2str(out.options.g_fname);
    else
        gfn = out.options.g_fname;
    end
    if isequal(gfn,'g_embed')
        gfn0 = out.options.inG.g_fname;
        if isa(gfn0,'function_handle')
            gfn0 = func2str(gfn0);
        end
        gfn = [gfn,' (',gfn0,')'];
        str{4} = [str{4},' (with delay embedding)'];
    end
    if isa(out.options.f_fname,'function_handle')
        ffn = func2str(out.options.f_fname);
    else
        ffn = out.options.f_fname;
    end
    if isequal(ffn,'f_embed')
        ffn0 = out.options.inF.f_fname;
        if isa(ffn0,'function_handle')
            ffn0 = func2str(ffn0);
        end
        ffn = [ffn,' (',ffn0,')'];
    end
    str{4} = sprintf([str{4},'\n ',...
        '    - ',so,' function: ',gfn,tmp,'\n ',...
        '    - evolution function: ',ffn]);
else
    str{4} = ['The model was static (no hidden states)','\n '];
    if isa(out.options.g_fname,'function_handle')
        gfn = func2str(out.options.g_fname);
    else
        gfn = out.options.g_fname;
    end
    str{4} = sprintf([str{4},'    - ',so,' function: ',gfn,tmp]);
end
str{5} = sprintf(['Bayesian log model evidences:','\n',...
    '     - full model: log p(y|m) > ',num2str(F,'%4.3e'),'\n',...
    '     - null hypothesis: log p(y|H0) = ',num2str(LLH0,'%4.3e')]);
if ~out.options.OnLine && out.dim.n >= 1 && ~isinf(out.options.priors.a_alpha) && ~isequal(out.options.priors.b_alpha,0)
    Fd = out.options.init.out.F;
    str{5} = sprintf([str{5},'\n ',...
        '    - deterministic variant: log p(y|m,eta=0) > ',num2str(Fd,'%4.3e') ]);
end
% str{5} = [str{5}, '\n'];

gsi = find([out.options.sources.type]==0);
bsi = find([out.options.sources.type]~=0);
many = length(out.options.sources)>1;
R2str = ['     - determin. coeff. (R2):  '];
LLstr = ['     - log-likelihood: '];
AICstr = ['     - AIC: '];
BICstr = ['     - BIC: '];
if ~isempty(gsi)
    R2str = [R2str,catnum2str(out.fit.R2,gsi,many,'%')];
    CAstr = [];
    LLstr = [LLstr,catnum2str(out.fit.LL,gsi,many)];
    AICstr = [AICstr,catnum2str(out.fit.AIC,gsi,many)];
    BICstr = [BICstr,catnum2str(out.fit.BIC,gsi,many)];
    separator = ', ';
else
    separator = [];
end
if  ~isempty(bsi)
    R2str = [R2str,separator,catnum2str(out.fit.R2,bsi,many,'%')];
    CAstr = ['     - balanced classif. acc.:  '];
    CAstr = [CAstr,catnum2str(out.fit.bacc,bsi,many,'%'),'\n'];
    LLstr = [LLstr,separator,catnum2str(out.fit.LL,bsi,many)];
    AICstr = [AICstr,separator,catnum2str(out.fit.AIC,bsi,many)];
    BICstr = [BICstr,separator,catnum2str(out.fit.BIC,bsi,many)];
end
R2str = [R2str,'\n'];
LLstr = [LLstr,'\n'];
AICstr = [AICstr,'\n'];
% BICstr = [BICstr,'\n'];
str{6} = ['Classical goodness-of-fit metrics:','\n',...
    R2str,...
    CAstr,...
    LLstr,...
    AICstr,...
    BICstr];

if newlines
    for i=1:length(str)
        str{i} = [str{i},'\n'];
    end
end

function str = catnum2str(x,ind,many,flag)
try,flag;catch,flag=[];end
str = [];
for i=1:length(ind)
    si=ind(i);
    if isequal(flag,'%')
        str = [str,', ',num2str(100*x(si),'%2.1f%%'),'%'];
    else
        str = [str,', ',num2str(x(si),'%4.3e')];
    end
    if many
        str = [str,' (source #',num2str(si),')'];
    end
end
str(1:2) = [];

author: MBB-team

repo: VBA-toolbox


text:
function [y] = spm_int_J(P,M,U)
% integrates a MIMO nonlinear system using the Jacobian
% FORMAT [y] = spm_int_J(P,M,U)
% P  - model parameters
% M  - model structure
% U  - input structure or matrix
%
% y  - (v x l)  response y = g(x,u,P)
%__________________________________________________________________________
% Integrates the MIMO system described by
%
%        dx/dt = f(x,u,P,M)
%        y     = g(x,u,P,M)
% or
%        dx/dt = f(x,u,P)
%        y     = g(x,u,P)
%
% using the update scheme:
%
%    x(t + dt) = x(t) + U*dx(t)/dt
%
%            U = (expm(dt*J) - I)*inv(J)
%            J = df/dx
%
% at input times.  This integration scheme evaluates the update matrix (Q)
% at each time point
%
%--------------------------------------------------------------------------
%
% SPM solvers or integrators
%
% spm_int_ode:  uses ode45 (or ode113) which are one and multi-step solvers
% respectively.  They can be used for any ODEs, where the Jacobian is
% unknown or difficult to compute; however, they may be slow.
%
% spm_int_J: uses an explicit Jacobian-based update scheme that preserves
% nonlinearities in the ODE: dx = (expm(dt*J) - I)*inv(J)*f.  If the
% equations of motion return J = df/dx, it will be used; otherwise it is
% evaluated numerically, using spm_diff at each time point.  This scheme is
% infallible but potentially slow, if the Jacobian is not available (calls
% spm_dx).
%
% spm_int_E: As for spm_int_J but uses the eigensystem of J(x(0)) to eschew
% matrix exponentials and inversion during the integration. It is probably
% the best compromise, if the Jacobian is not available explicitly.
%
% spm_int_B: As for spm_int_J but uses a first-order approximation to J
% based on J(x(t)) = J(x(0)) + dJdx*x(t).
%
% spm_int_L: As for spm_int_B but uses J(x(0)).
%
% spm_int_U: like spm_int_J but only evaluates J when the input changes.
% This can be useful if input changes are sparse (e.g., boxcar functions).
% It is used primarily for integrating EEG models
%
% spm_int:  Fast integrator that uses a bilinear approximation to the
% Jacobian evaluated using spm_bireduce. This routine will also allow for
% sparse sampling of the solution and delays in observing outputs. It is
% used primarily for integrating fMRI models
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_int_J.m 6801 2016-05-29 19:18:06Z karl $

% convert U to U.u if necessary and M(1) to M
%--------------------------------------------------------------------------
if ~isstruct(U), u.u = U; U = u; end
try, dt = U.dt;  catch, dt = 1;  end
M       = M(1);

% state equation; add [0] states if not specified
%--------------------------------------------------------------------------
try
    f   = fcnchk(M.f,'x','u','P','M');
catch
    f   = @(x,v,P,M)x;sparse(0,1);
    M.n = 0;
    M.x = sparse(0,0);
end

% and output nonlinearity
%--------------------------------------------------------------------------
try
    g   = fcnchk(M.g,'x','u','P','M');
catch
    g   = @(x,v,P,M)x;
end

% Initial states and inputs
%--------------------------------------------------------------------------
try
    u   = U.u(1,:);
catch
    u   = sparse(1,M.m);
end
try
    x   = M.x;
catch
    x   = sparse(0,1);
    M.x = x;
end

% check function format
%--------------------------------------------------------------------------
if ~isa(f,'function_handle')
    try
        f(x,u,P,M);
    catch
        f = inline(char(f),'x','v','P','M');
    end
end
if ~isa(g,'function_handle')
    try
        g(x,u,P,M);
    catch
        g = inline(char(g),'x','v','P','M');
    end
end

% default delay operator
%--------------------------------------------------------------------------
D = 1;

% integrate
%==========================================================================
for i = 1:size(U.u,1)

    % input
    %----------------------------------------------------------------------
    try
        u = U.u(i,:);
    end

    % dx(t)/dt and Jacobian df/dx
    %----------------------------------------------------------------------
    if nargout(f) >= 3
        [fx,dfdx,D] = f(x,u,P,M);

    elseif nargout(f) == 2
        [fx,dfdx]   = f(x,u,P,M);

    else
        fx          = f(x,u,P,M);
        dfdx        = spm_cat(spm_diff(f,x,u,P,M,1));
    end

    % update dx = (expm(dt*J) - I)*inv(J)*fx
    %----------------------------------------------------------------------
    x      = spm_unvec(spm_vec(x) + spm_dx(D*dfdx,D*fx,dt),x);

    % output - implement g(x)
    %----------------------------------------------------------------------
    if nargin(g) > 3
        y(:,i) = spm_vec(g(x,u,P,M));
    else
        y(:,i) = spm_vec(g(x,u,P));
    end

end

% transpose
%--------------------------------------------------------------------------
y      = real(y');

author: spm


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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym dawson (@var{x})
%% Symbolic Dawson (scaled imaginary error) function.
%%
%% Example:
%% @example
%% @group
%% syms x
%% dawson (x)
%%   @result{} ans = (sym)
%%             2
%%           -x
%%       √π⋅ℯ   ⋅erfi(x)
%%       ───────────────
%%              2
%% @end group
%% @end example
%% @seealso{dawson, @@sym/erfc, @@sym/erf, @@sym/erfcx, @@sym/erfi, @@sym/erfinv, @@sym/erfcinv}
%% @end defmethod

function y = dawson(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('lambda a: exp(-a**2)*erfi(a)*(sqrt(S(pi))/2)', x);
end

%!test
%! if (exist ('OCTAVE_VERSION', 'builtin'))
%! % dawson missing on Matlab, Issue #742
%! A = dawson([1 2]);
%! B = double(dawson(sym([1 2])));
%! assert(A, B, -eps)
%! end

repo: octsympy


text:
function [fig] = plot_trajectories_offline(pos_history, N, ...
                colors, fontsize, map)
% plot_trajectpries_offline - Plot trajectories of the agents.

x0 = 10;
y0 = 10;
width = 400;
height = 400;

fig = figure('Name','Offline swarm trajectories','NumberTitle','off');
hold on;
grid on;
box on;

% Plot environment
if ~isempty(map)
    fig = draw_cylinders(fig,map);
end

% Plot trajetcories
for agent = 1:N
    hold on;

    if ~isempty(colors)
        plot3(pos_history(:,(agent-1)*3+2), pos_history(:,(agent-1)*3+1), ...
            - pos_history(:,(agent-1)*3+3), 'Color', colors(:,agent));
    else
        plot3(pos_history(:,(agent-1)*3+2), pos_history(:,(agent-1)*3+1), ...
            - pos_history(:,(agent-1)*3+3));
    end

end

xlabel('Y Position [m]','fontsize',fontsize);
ylabel('X Position [m]','fontsize',fontsize);
zlabel('Z Position [m]','fontsize',fontsize);
view(2);

set(fig,'units','pixels','position',[x0,y0,width,height]);
set(fig,'PaperPositionMode','auto');

end

author: lis-epfl

repo: swarmlab

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


text:
function [M,padpre] = padImage(MData, pof2, padval)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

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

dim = size(MData);
s = 2^pof2;
newdim = ceil(dim/s)*s;

paddim = newdim-dim;
padpre = floor(paddim/2);
padpost = paddim-padpre;

if exist('padarray.m','file')
    M = padarray(MData,padpre,padval,'pre');
    M = padarray(M,padpost,padval,'post');
else
    M = padarray_oct(MData,padpre,padval,'pre');
    M = padarray_oct(M,padpost,padval,'post');
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [out] = meanangle(in,dim,sens);

% MEANANGLE will calculate the mean of a set of angles (in degrees) based
% on polar considerations.
%
% Usage: [out] = meanangle(in,dim)
%
% in is a vector or matrix of angles (in degrees)
% out is the mean of these angles along the dimension dim
%
% If dim is not specified, the first non-singleton dimension is used.
%
% A sensitivity factor is used to determine oppositeness, and is how close
% the mean of the complex representations of the angles can be to zero
% before being called zero.  For nearly all cases, this parameter is fine
% at its default (1e-12), but it can be readjusted as a third parameter if
% necessary:
%
% [out] = meanangle(in,dim,sensitivity)
%
% Written by J.A. Dunne, 10-20-05
%

if nargin<3
    sens = 1e-12;
end

if nargin<2
    ind = min(find(size(in)>1));
    if isempty(ind)
        %This is a scalar
        out = in;
        return
    end
    dim = ind;
end

in = in * pi/180;

in = exp(i*in);
mid = mean(in,dim);
out = atan2(imag(mid),real(mid))*180/pi;
out(abs(mid)<sens) = nan;

author: Sable

repo: mcbench-benchmarks


text:
%% This file is the main file of camparing the noise sensitivity of iSINDy
% and SINDy-PI method. This file takes a long time to run.
%
% Last Update: 2019/07/17
% Coded By: K

%% Close all, clear all, clc
close all;clear all; clc;
addpath('./Functions')
addpath('Datas')
set(0,'defaulttextInterpreter','latex')
%% Simulate the budworm population growth and gather the simulation data

% Define the system parameters
jx=0.6;Vmax=1.5;Km=0.3;

% Determine the simulation time step and time span
dt=0.1; T=5; tspan=0:dt:T;

% Define noise level and add gaussian noise to the data
noise=0;

% Define whehter you have control, if you have it, please define it
Control=0;u=0;

%Define whether you want to shuffel the final data
Shuffle=0;

%% Set up some parameters
% Get the number of states we have
n_state=1;

% Define the control input(Should be zero in our example)
n_control=0;

% Choose whether you want to display actual ODE or not
disp_actual_ode=1;

% If the ODEs you want to display is the actual underlyting dynamics of the
% system, please set actual as 1
actual=1;

% Print the actual ODE we try to discover
Print_ODEs(@(t,y)MMK_ODE(t,y,jx,Vmax,Km),n_state,n_control,disp_actual_ode,actual);

% Create symbolic states
dz=sym('dz',[n_state,1]);

% Now we first create the parameters of the function right hand side
Highest_Poly_Order_Guess=1;
Highest_Trig_Order_Guess=0;
Highest_U_Order_Guess=0;

% Then create the right hand side library parameters
Highest_Poly_Order=4;
Highest_Trig_Order=0;
Highest_U_Order=0;
Highest_dPoly_Order=1;

% Determine whether you want to normalize your libary or not. 1 is yes and
% 0 is no.
NormalizeLib=0;

% Determine the model selection method
Model_Selection_Method=1;
% Determine the prediction step
Prediction_Steps=0;

% The follwoing parameters will determine the for loop and loop through the
% different noise level
d_percent=1;
percent=0;
percent_start=1;
% Must be smaller or equal to 24
percent_end=24;

% Determine a vector to store whether SINDy-PI identification is correct
is_Right=zeros(percent_end-percent_start+1,1);

% Define a matrix to store the null space vector
Xi_ns=cell(percent_end-percent_start+1,1);

% Turn off the rank deficient warning
MSGID='MATLAB:rankDeficientMatrix';
warning('off',MSGID);

% Determine how many iteration you need
Final_pin=30;

% Choose which denoise method you want:
% 1 is adding Guassian noise to actual derivative, 2 is direct
% finite difference, 3 is TVRegDiff.
Noisy_Data_Generation_Method=3;

% Determine how many initial conditions you need
init_num=2400;

% Create the new directory to save the function files
FolderName=strcat('Result_i_SINDy_Derivative_Method',num2str(Noisy_Data_Generation_Method),'_PridictionStep_',num2str(Prediction_Steps),'_IniNum_',num2str(init_num));
[fld_status, fld_msg, fld_msgID]=mkdir('TempFunctions_iSINDy');
[fld_status, fld_msg, fld_msgID]=mkdir(FolderName);
addpath('TempFunctions_iSINDy')

%% Now calculate the sensitivity of the iSINDy and SINDy-PI to the noise.

% Print the start process
fprintf('Start calculating ....\n\n\n')

for total=1:Final_pin
    percent=0;

    for percent_iter=percent_start:d_percent:percent_end
        % Set the pin one step forward
        percent=percent+1;

        % Set the noise level
        noise_level(percent,1)=Determine_Noise_Level(percent_iter);

        % Print the current noise level
        fprintf('Using the noise level as %i on implicit-SINDy.\n',noise_level(percent,1))

        % Read files, load the training data and testing data
        File_Name=strcat('Data_Iter_',num2str(total),'_NoiseLevel_',num2str(noise_level(percent,1)),'_NoiseMethod_',num2str(Noisy_Data_Generation_Method),'_GoodLuck.mat');
        load(File_Name)
        [size1,size2]=size(xt_test);

        %% Run implicit-SINDy (some function in this section is obtained from the github code iSINDy )
        fprintf('\n\n\n\n USing implicit-SINDy now...\n')

        % Sweep through the states
        %tic
        for iter=1:n_state
            fprintf('\n \n Calculating the %i expression...\n',iter)

            % Build library data
            [iSINDy_Data,iSINDy_Struct]=SINDyLib(Data,dData(:,iter),u,Highest_Poly_Order,Highest_Trig_Order,Highest_U_Order,Highest_dPoly_Order);

            % Use Null-Space method and set the tolerance
            tol = 1e-5;

            % Plotting option
            pflag=0;

            % Print process
            fprintf('\n \n \t Calculating the null space...\n')
            try
                % Sparse Regression
                [Xi_ns{percent,1}, indTheta, lambdavec, numterms, errorv] = ADMpareto(iSINDy_Data, tol, pflag);
                Xi_dum=zeros(10,28);
                Xi_dum=Xi_ns{percent,1};
                %Define the variable to store the ODE expression and score of each
                % expression. This is necessary to perform the parallel for loop.
                ODEs_iS=cell(n_state,size(Xi_dum,2)-1);
                if percent==1 && iter==1
                    Score_iS=NaN*ones(percent_end-percent_start+1,size(Xi_dum,2)-1);
                end
                dz_dum=dz(iter);

                % Sweep the null space and store the score
                parfor j=1:(size(Xi_dum,2)-1)

                    % Print which vector we are working on
                    fprintf('\t\t Sweeping the %i null sapce vector...\n',j)

                    % Validate each sparse vector
                    try
                        % Solve for the ODE expression
                        ODEs_iS{iter,j}=solve(vpa(cell2sym(iSINDy_Struct)*Xi_dum(:,j))==0,dz_dum);

                        % Store the previous result into a matrix
                        Generate_ODE_RHS(ODEs_iS{iter,j},n_state,n_control,strcat('TempFunctions_iSINDy/ParForiS',num2str(j)));

                        % Calculate the accuracy of the file
                        Test_Score=0;
                        Test_Score=Get_Score(dxt_test,xt_test,u,Control,tspan,Shuffle,strcat('ParForiS',num2str(j)),Prediction_Steps,dt,size1,size2,Model_Selection_Method);
                        Score_iS(percent,j)=Test_Score;
                    catch
                        ODEs_iS{iter,j}=NaN;
                        Score_iS(percent,j)=NaN;
                    end
                end
            catch
                fprintf('\n\t\t\tSomething is wrong....\n')
            end

        end

    end

    % Store the result for this iteration

    % Store the minimum of implicit SINDy
    for pin1=1:size(Score_iS,1)
        [min_iS_Val(total,pin1),Index1(total,pin1)]=min(Score_iS(pin1,:));

    end

    % Test whether the correct structure is identified
    for pin1=1:size(Score_iS,1)
        Xi_dum=Xi_ns{pin1,1};
        Coffs=Xi_dum(:,Index1(total,pin1));
        Coffs_dum=Coffs~=0;
        if Coffs_dum==[1;1;0;0;0;1;1;0;0;0]
            is_Right(pin1,1)=is_Right(pin1,1)+1;
            Coffs_norm=Coffs/Coffs(7);
            Coffs_error(total,pin1)=norm(Coffs_norm-[-0.18;0.9;0;0;0;0.3;1;0;0;0]);
        elseif Coffs_dum==[0;1;1;0;0;0;1;1;0;0]
            is_Right(pin1,1)=is_Right(pin1,1)+1;
            Coffs_norm=Coffs/Coffs(8);
            Coffs_error(total,pin1)=norm(Coffs_norm-[0;-0.18;0.9;0;0;0;0.3;1;0;0]);
        elseif Coffs_dum==[0;0;1;1;0;0;0;1;1;0]
            is_Right(pin1,1)=is_Right(pin1,1)+1;
            Coffs_norm=Coffs/Coffs(9);
            Coffs_error(total,pin1)=norm(Coffs_norm-[0;0;-0.18;0.9;0;0;0;0.3;1;0]);
        elseif Coffs_dum==[0;0;0;1;1;0;0;0;1;1]
            is_Right(pin1,1)=is_Right(pin1,1)+1;
            Coffs_norm=Coffs/Coffs(10);
            Coffs_error(total,pin1)=norm(Coffs_norm-[0;0;0;-0.18;0.9;0;0;0;0.3;1]);
        end
    end

    for pin1=1:size(Score_iS,1)
        try
            Num=sum(Coffs_error(:,pin1)~=0);
            Ave_Coff_Error(pin1,1)=sum(Coffs_error(:,pin1))/Num;
        catch
            Ave_Coff_Error(pin1,1)=NaN;
        end
    end

    %% Plot the result
    figure(1)
    plot(noise_level,log(1+mean(min_iS_Val,1)),'o','linewidth',2,'color','blue')
    drawnow
    grid on
    title('Peformance Comparison','FontSize',18)
    xlabel('Noise Level $\sigma$','FontSize', 18)
    ylabel('Average L2 Norm','FontSize', 18)
    set(gca,'FontSize',18);
    set(gcf,'Position',[100 100 600 400]);
    set(gcf,'PaperPositionMode','auto');

    figure(2)
    plot(noise_level,is_Right/total,'o')
    drawnow
    grid on
    title('Success Rate','FontSize',18)
    xlabel('Noise Level $\sigma$','FontSize', 18)
    ylabel('Success Rate','FontSize', 18)
    ylim([0 1])
    set(gca,'FontSize',18);
    set(gcf,'Position',[100 100 600 400]);
    set(gcf,'PaperPositionMode','auto');
    set(gca, 'XScale', 'log')

    figure(3)
    plot(noise_level,Ave_Coff_Error,'o')
    drawnow
    grid on
    title('Average Parameter Error','FontSize',18)
    xlabel('Noise Level $\sigma$','FontSize', 18)
    ylabel('L2 norm','FontSize', 18)
    ylim([0 total])
    set(gca,'FontSize',18);
    set(gcf,'Position',[100 100 600 400]);
    set(gcf,'PaperPositionMode','auto');
    set(gca, 'XScale', 'log')

    % Save the calculation result of current iteration
    cc=clock;
    ResultName=strcat(FolderName,'/implicit_SINDY_Result',num2str(total),'__',num2str(cc(3)),'_',num2str(cc(4)),'_',num2str(cc(5)),'_',num2str(round(cc(6))),'.mat');
    save(ResultName,'Score_iS','Xi_ns')

end

author: dynamicslab

repo: SINDy-PI


text:
function poly2 = clipConvexPolygon3dHP(poly, plane)
%CLIPCONVEXPOLYGON3DHP Clip a convex 3D polygon with Half-space
%
%   POLY2 = clipConvexPolygon3dHP(POLY, PLANE)
%   POLY is a N-by-3 array of points, and PLANE is given as:
%   [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2].
%   The result POLY2 is also an array of 3d points, sometimes smaller than
%   poly, and that can be 0-by-3 (empty polygon).
%
%   POLY2 = clipConvexPolygon3dHP(POLY, PT0, NORMAL)
%   uses plane with normal NORMAL and containing point PT0.
%
%
%   See also:
%   polygons3d, polyhedra
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2007-01-05
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

%   HISTORY
%   2007/14/09 fix postprocessing of last point

% ensure last point is the same as the first one
if sum(poly(end, :) == poly(1,:)) ~= 3
    poly = [poly; poly(1,:)];
end

% initialize empty polygon
poly2 = zeros(0, 2);

% compute visible points
below = isBelowPlane(poly, plane);

% case of empty polygon
if sum(below) == 0
    return;
end

% case of totally clipped polygon
if sum(below) == length(below)
    poly2 = poly;
    return;
end

% indices of edges intersecting the plane
ind = find(below ~= below([2:end 1]));

% compute intersection points: they are 2 for a convex polygon
lines = createLine3d(poly(ind, :), poly(ind+1, :));
pInt = intersectLinePlane(lines, plane);

% insert intersection points and remove invisible points
if below(1)
    poly2 = [poly(1:ind(1), :); pInt; poly(ind(2)+1:end, :)];
else
    poly2 = [pInt(1, :); poly(ind(1)+1:ind(2), :); pInt(2, :)];
end

% remove last point if it is the same as the first one
if sum(poly2(end, :) == poly2(1,:)) == 3
    poly2(end, :) = [];
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [success]=checksignres_favar(Sj,fj)
 % check here for success only

% if there are no sign restrictions on this shock, don't do anything and automatically count as success
if isempty(Sj)
success=1;
% if there are sign restriction on this shock, check them
else
% check if the restrictions hold
   % if yes, count as success
   if all(Sj*fj>=0)
   success=1;
%    % if the restrictions do not hold, there may still be a possibility by switching the sign of qj
%    elseif all(Sj*(-fj)>=0)
%    qj=-qj;
%    success=1;
   % else, if there is no way to have qj succesful, count as a fail
   else
   success=0;
   end
end

author: european-central-bank

repo: BEAR-toolbox


text:
function K = hole_correlation ( s, t )

%*****************************************************************************80
%
%% HOLE_CORRELATION evaluates the hole correlation function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real S(*), T(*), pairs of argument values.
%
%    Output, real K(*), the correlation function values
%
  d = abs ( s - t );

  K = ( 1 - d ) .* exp ( - d );

  return
end


text:
function [p] = linortfit2(xdata, ydata)
% LINORTFIT2  Fit a line to data by ORTHOGONAL least-squares.
%    P = LINORTFIT2(X,Y) finds the coefficients of a 1st-order polynomial
%    that best fits the data (X,Y) in an ORTHOGONAL least-squares sense.
%    Consider the line P(1)*t + P(2), and the minimum (Euclidean) distance
%    between this line and each datapoint [X(i) Y(i)] -- LINORTFIT2 finds
%    P(1) and P(2) such that the sum of squared distances is minimized.

if ~isequal(size(xdata), size(ydata))
    error('linortfit2:XYSizeMismatch',...
          'X and Y vectors must be the same size.');
end

[N,C] = linortfitn([xdata(:) ydata(:)]);

% The hyperplane given by N * [x; y] + C == 0 is optimal.
% Convert to the form  y = p(1)*x + p(2), just like polyfit(xdata,ydata,1).
p = - [N(1)  C] / N(2);


text:
%% Efficient subpixel image registration by cross-correlation.
% Registers two images (2-D rigid translation) within a  fraction
% of a pixel specified by the user. Instead of computing a zero-padded FFT
% (fast Fourier transform), this code uses selective upsampling by a
% matrix-multiply DFT (discrete FT) to dramatically reduce computation time and memory
% without sacrificing accuracy. With this procedure all the image points are used to
% compute the upsampled cross-correlation in a very small neighborhood around its peak. This
% algorithm is referred to as the single-step DFT algorithm in [1].
%
% [1] Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup,
% "Efficient subpixel image registration algorithms," Opt. Lett. 33,
% 156-158 (2008).

%% Syntax
% The code receives the FFT of the reference and the shifted images, and an
% (integer) upsampling factor. The code expects FFTs with DC in (1,1) so do not use
% fftshift.
%
%    output = dftregistration(fft2(f),fft2(g),usfac);
%
% The images are registered to within 1/usfac of a pixel.
%
% output(1) is the normalized root-mean-squared error (NRMSE) [1] between f and
% g.
%
% output(2) is the global phase difference between the two images (should be
% zero if images are real-valued and non-negative).
%
% output(3) and output(4) are the row and column shifts between f and g respectively.
%
%    [output Greg] = dftregistration(fft2(f),fft2(g),usfac);
%
% Greg is an optional output, it returns the Fourier transform of the registered version of g,
% where the global phase difference [output(2)] is also compensated.

%% Obtain a reference and shifted images
% To illustrate the use of the algorithm, lets obtain a reference and a
% shifted image. First we read the reference image f(x,y)
f = im2double(imread('cameraman.tif'));
%%
% Define g(x,y) as a version of f(x,y) shifted by fractional values of a
% pixel and multiplied by a global phase.
deltar = 3.48574;
deltac = 8.73837;
phase = 2;
[nr,nc]=size(f);
Nr = ifftshift([-fix(nr/2):ceil(nr/2)-1]);
Nc = ifftshift([-fix(nc/2):ceil(nc/2)-1]);
[Nc,Nr] = meshgrid(Nc,Nr);
g = ifft2(fft2(f).*exp(i*2*pi*(deltar*Nr/nr+deltac*Nc/nc))).*exp(-i*phase);
figure(1);
subplot(1,2,1);
imshow(abs(f));
title('Reference image, f(x,y)')
subplot(1,2,2);
imshow(abs(g));
title('Shifted image, g(x,y)')
%%
% We have shifted the image by 8.73837 and 3.48574 pixels in the x and
% y direction, respectively, and added a phase of 2 radians to g(x,y). The shift
% was implemented by applying a linear phase on its
% FT, thus we have assumed that the images wrap around (features leaving one
% side of the window reappear on the opposite side) and that the image is band-limited
% (interpolated by a sinc function). Cross-correlation image registration by DFTs
% (both the matrix-multiply DFT and the zero-padded FFT) share these assumptions.
%
% This registration technique is well suited to compare images that are captured
% in Fourier domain (i.e. to evaluate an image reconstruction by holography
% or phase retrieval) which are strictly band-limited and exhibit the
% wrap-around effect.
%
% Even though the registration code assumes band-limited images that wrap around, we
% have obtained very good results when applying it to
% band-limited microscope images, and aliased imagery. That is when shifting
% the image brings in new content instead of wrapping it around or when the
% images are not band-limited.

%% Sample Image Registration
% dftregistration.m receives the FT of f and g and the upsampling factor.
% The code expects DC of the FTs at (1,1) so don't use fftshift.
%
% We now use the image registration code to register f and g within 0.01
% pixels by specifying an upsampling parameter of 100
[output Greg] = dftregistration(fft2(f),fft2(g),100);
display(output),
%%
% The pixel shift error (difference between the true and obtained shifts)
% is 0.0016 and 0.0043 in the x and y directions respectively. Well within
% the expected accuracy of 0.01. Notice that using the conventional zero-padded
% FFT approach with the same accuracy, would
% require computation of a 25,600x25,600 FFT, which would require more than
% 19 Gbytes of RAM and a very comfortable chair.
%
% The following plot shows the reference image and the registered image.
figure(1);
subplot(1,2,1);
imshow(abs(f));
title('Reference image, f(x,y)')
subplot(1,2,2);
imshow(abs(ifft2(Greg)));
title('Registered image, gr(x,y)')
%% Disclaimer
% I have made every effort to evaluate the proper working of this code
% under many different conditions. However, it is the responsibility of
% the user to ensure that this registration code is adequate for their
% application.
%
% Feel free to e-mail me with questions or comments.


text:
function sinogram = myRadon(image,thetas)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% radon transformation -> schlegel & bille 9.1.1
% written by Mark Bangert
% m.bangert@dkfz.de 2011

numOfAngularProjections  = length(thetas);
numOfParallelProjections = size(image,1);

sinogram = zeros(numOfParallelProjections,numOfAngularProjections);

% loop over the number of angles
for i = 1:length(thetas)

   % rotate image
   tmpImage      = imrotate(image,-thetas(i),'bilinear','crop');

   % fill sinogram
   sinogram(:,i) = sum(tmpImage,2);

   % visualization on the fly
   imagesc(sinogram);
   drawnow

end


text:
function results = vl_test_hog(varargin)
% VL_TEST_HOG
vl_test_init ;

function s = setup()
s.im = im2single(vl_impattern('roofs1')) ;
[x,y]= meshgrid(linspace(-1,1,128)) ;
s.round = single(x.^2+y.^2);
s.imSmall = s.im(1:128,1:128,:) ;
s.imSmall = s.im ;
s.imSmallFlipped = s.imSmall(:,end:-1:1,:) ;

function test_basic_call(s)
cellSize = 8 ;
hog = vl_hog(s.im, cellSize) ;

function test_bilinear_orientations(s)
cellSize = 8 ;
vl_hog(s.im, cellSize, 'bilinearOrientations') ;

function test_variants_and_flipping(s)
variants = {'uoctti', 'dalaltriggs'} ;
numOrientationsRange = 3:9 ;
cellSize = 8 ;

for cellSize = [4 8 16]
  for i=1:numel(variants)
    for j=1:numel(numOrientationsRange)
      args = {'bilinearOrientations', ...
              'variant', variants{i}, ...
              'numOrientations', numOrientationsRange(j)} ;
      hog = vl_hog(s.imSmall, cellSize, args{:}) ;
      perm = vl_hog('permutation', args{:}) ;
      hog1 = vl_hog(s.imSmallFlipped, cellSize, args{:}) ;
      hog2 = hog(:,end:-1:1,perm) ;
      %norm(hog1(:)-hog2(:))
      vl_assert_almost_equal(hog1,hog2,1e-3) ;
    end
  end
end

function test_polar(s)
cellSize = 8 ;
im = s.round ;
for b = [0 1]
  if b
    args = {'bilinearOrientations'} ;
  else
    args = {} ;
  end
  hog1 = vl_hog(im, cellSize, args{:}) ;
  [ix,iy] = vl_grad(im) ;
  m = sqrt(ix.^2 + iy.^2) ;
  a = atan2(iy,ix) ;
  m(:,[1 end]) = 0 ;
  m([1 end],:) = 0 ;
  hog2 = vl_hog(cat(3,m,a), cellSize, 'DirectedPolarField', args{:}) ;
  vl_assert_almost_equal(hog1,hog2,norm(hog1(:))/1000) ;
end

author: Cloud-CV

repo: object-proposals


text:
function [status, results] = mrtrix_mrconvert(orig_data_file, mif_filename, bkgrnd, verbose)
%
% Convert a nifti image to a mrtrix .mif image.
% Example
% mrtrix_mrconvert('t1.nii.gz', 't1.mif')

if notDefined('bkgrnd'),  bkgrnd  = false;end
if notDefined('verbose'), verbose = true; end

% First, deal with the possibility that the original data is compressed:
if strcmp(  orig_data_file(end-1:end), 'gz')
    % Try uncompressing it with gunzip:
    [status, results] = system(sprintf('gunzip %s', orig_data_file));
    % Need to change the target for conversion:
    convert_file = orig_data_file(1:end-3);
    if status ~= 0 && ~exist(convert_file,'file')
        error('There was some problem extracting the raw data file. \n UNIX Message: %s \n',results);
    end
else
    convert_file = orig_data_file;
end

cmd_str = sprintf('mrconvert %s %s', convert_file, mif_filename);

[status,results] = mrtrix_cmd(cmd_str, bkgrnd, verbose);

% Rezip it:
if orig_data_file(end-1:end) == 'gz'

    [s,r] = system(sprintf('gzip %s', convert_file));
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function nksub = ksubset_enum ( k, n )

%*****************************************************************************80
%
%% KSUBSET_ENUM enumerates the K element subsets of an N set.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer K, the number of elements each K subset must
%    have. 0 <= K <= N.
%
%    Input, integer N, the number of elements in the master set.
%    0 <= N.
%
%    Output, integer NKSUB, the number of distinct elements.
%
  nksub = i4_choose ( n, k );

  return
end

repo: jburkardt-m


text:
function view = stimPlusTSeriesMovie(view,stim,scans,slices,varargin);
%  view = stimPlusTSeriesMovie(view,stim,[scans,slices,options]);
%
% Similar to makeTSeriesMovie, makes a movie with the tSeries
% overlayed over anatomies, but with some additional features:
%
%   1) Shows next to the movie the stimuli that were showing
%   at each frame;
%
%   2) Has a simple interface allowing the threshold to be
%   set for what (normalized) values of the tSeries are
%   overlayed over the anatomies.
%
% stim must be a struct containing the following fields:
%   images: cell specifying which images were shown. Each entry
%           in the cell can be a matrix containing the image, or
%           a string specifying the path to load the image;
%           [If left empty, will try to load these from the third
%           column of the parfiles assigned to the selected scans]
%   frames: array specifying which frames from the tSeries to use
%           [default is 1:numFrames, but if, e.g., many images appear
%           within one TR you may 'double up' those entries of the
%           frames field]
%   order:  array specifying which image (index in the images cell)
%           to show for each frame. Should be the same length as
%           the frames field. [If empty, will use the conditions
%           specified in the parfiles assigned to the selected scans]
%
%
% Defaults to the current scan and slices, though multiple
% scans/slices may be specified.
%
%
% 05/04 ras.

% the callbacks from the UI controls give as the
% first argument a string specifying the sort of
% callback. Check:
if ischar(view)
    % it's a callback, not a call to open a new movie
    switch view
        case 'play', playStimMovie; return;
        case 'rebuild', rebuildStimMovie;
        case 'export', exportStimMovie; return;
        otherwise, return;
    end
elseif ~isstruct(view)
    % if it's not a struct, something's wrong
    help stimPlusTSeriesMovie;
    error('stimPlusTSeriesMovie called wrong.');
end

% --------- Everything below is code to build the movie ----------- %

if ieNotDefined('scans')    scans = viewGet(view,'curScan');   end
if ieNotDefined('slices')   slices = viewGet(view,'curSlice'); end

%%%%% default params
thresh = 0.6;       % normalized func intensity val
detrend = 0;        % detrending takes a long time on a whole tSeries
cmap = hot(256);    % color map for tSeries overlay
movieRate = 6;     % rate in frames per second to play movie
TR = viewGet(view,'TR');

%%%%% parse options
% put code here

% check that the stim struct has the proper fields
if ~isfield(stim,'images') | ~isfield(stim,'frames') | ~isfield(stim,'order')
    help stimPlusTSeriesMovie;
    error('Second argument does not have proper fields.');
end

% allow for some fields to be empty; this means get the information
% from the parfiles:
if isempty(stim.images)
    trials = exp5_concatParfiles(view,scans);
    stim.images = trials.image;
end

if isempty(stim.frames)
    global dataTYPES;
    dt = viewGet(view,'curdt');
    nFrames = 0;
    for s = scans
        nFrames = nFrames + dataTYPES(dt).scanParams(s).nFrames;
    end
    stim.frames = 1:nFrames;
end

if isempty(stim.order)
   % assume the conditions in the parfiles point to the index
   % of which image to use in trials.images -- interpolate
   % to make order a vector showing which condition at each frame
   trials = exp5_concatParfiles(view,scans);
   ok = find(trials.onset < nFrames);
   stim.order = NaN*ones(1,nFrames);
   stim.order(trials.onset(ok)) = trials.cond(ok);
   % we've set the onset frames of each cond; now fill in the intervals
   fillin = find(isnan(stim.order));
   for j = fillin
       lastOnset = find((1:nFrames < j) & ~isnan(stim.order));
       stim.order(j) = stim.order(lastOnset(end));
   end
end

% load any images specified by paths in stim.images
for i = 1:length(stim.images)
    if ischar(stim.images{i})
        if exist(stim.images{i},'file')
            stim.images{i} = imread(stim.images{i},'jpg');
        else
            error(sprintf('File %s not found.',stim.images{i}));
        end
    end
end

% open the movie window; add handles to stim struct
stim.cmap = cmap;
stim = openMovieWin(view,stim,scans,slices);

% Get, rescale anatomy image
anatIm = makeMontage(view.anat,slices);
anatClip = viewGet(view,'anatClip');
minVal = double(min(anatIm(:)));
maxVal = double(max(anatIm(:)));
anatClipMin = min(anatClip)*(maxVal-minVal) + minVal;
anatClipMax = max(anatClip)*(maxVal-minVal) + minVal;
anatIm = (rescale2(double(anatIm),[anatClipMin,anatClipMax],[1,256]));
anatIm = normalize(anatIm,0,1);
imSize = size(anatIm);

% make all the stimulus images squares with the edge equal to
% the y-size of the fnctional image (this will ensure everything
% stays the same relative size):
for i = 1:length(stim.images)
    if ~isequal(size(stim.images{i}),[imSize(1) imSize(1)])
        stim.images{i} = imresize(stim.images{i},[imSize(1) imSize(1)]);
    end
end

% make a spacer b/w the functional img and stimulus img
spacer = uint8(205*ones(imSize(1),50,3));

% load the tSeries for the selected scans and slices
fprintf('Loading tSeries...');
for s = 1:length(scans)
    for i = 1:length(slices)
        tS = loadtSeries(view,scans(s),slices(i));
        if s==1
            tSeriesAll{i} = tS;
        else
            tSeriesAll{i} = [tSeriesAll{i} tS];
        end
    end
end
fprintf('done.\n');

% ----- build the movie ----- %
nFrames = length(stim.frames);
dims = sliceDims(view,scans(1));
for f = 1:nFrames

    % construct the functional img (anat + overlay) for this frame
    func = zeros(imSize(1),imSize(2),3); % truecolor image
    clear tSeries;
    for ii = 1:length(slices)
        tSeries(:,:,ii) = reshape(tSeriesAll{ii}(f,:),dims);
    end
    tSeries = makeMontage(tSeries);
    overlay = upSampleRep(tSeries,size(anatIm));
    overlay = normalize(overlay,0,1);
    mask = find(overlay >= thresh);
    overlay(mask) = normalize(overlay(mask),1,256); % indices into cmap
    overlay = round(overlay);
    for col = 1:3
        tmp = anatIm;
        tmp(mask) = cmap(overlay(mask),col);
        func(:,:,col) = tmp;
    end
    func = uint8(256 .* func);

    % get stimulus img for this frame
    whichStim = stim.order(f);
    if whichStim==0  % null condition, maybe a better way to do this
        stimulus = zeros(imSize(1),imSize(1));
    else
        stimulus = stim.images{whichStim};
    end
    stimulus = repmat(stimulus,[1 1 3]);
    stimulus = uint8(stimulus);

    % merge the images together, display and grab a movie frame
    axes(stim.handles.axes);
    img = [stimulus spacer func];
    imshow(img);
    htmp = text(20,20,num2str(f)); set(htmp,'Color','r');
    title(sprintf('%i secs',f*TR),'FontSize',16);
    stim.movie(f) = getframe;

    % might also be good to have this be attached to
    % the view like makeTSeriesMovie, but I have to
    % do some other things to ensure full compatibility:
    view.ui.movie.movie(f) = getframe;
end

colorbar horiz;
colormap(stim.cmap);

% set some important params on the stim struct; stash in fig's UserData
stim.viewName = view.name;
stim.thresh = thresh;
stim.detrend = detrend;
stim.movieRate = movieRate;
set(stim.handles.figure,'UserData',stim);

return
% /----------------------------------------------------------------------/ %

% /----------------------------------------------------------------------/ %
function stim = openMovieWin(view,stim,scans,slices);
% opens the movie window and adds the handles
stim.handles.figure = figure('Units','Normalized','Position',[.2 .2 .6 .6]);
set(gcf,'Name',['tSeries movie, scan ',num2str(scans),', slices ',num2str(slices)]);
% set(gcf,'DefaultTextColor','blue');
stim.handles.axes = axes('Position',[.05 .3 .75 .6]);
cbstr = 's = get(gcf,''UserData''); movie(s.movie,1,s.movieRate);';
stim.handles.play = uicontrol('Style','pushbutton','Units','Normalized',...
                              'Position',[.83 .7 .12 .08],'String','Play',...
                              'Callback',cbstr);
stim.handles.rebuild = uicontrol('Style','pushbutton','Units','Normalized',...
                              'Position',[.83 .5 .12 .08],'String','Rebuild',...
                              'Callback','stimPlusTSeriesMovie(''rebuild'');');
stim.handles.export = uicontrol('Style','pushbutton','Units','Normalized',...
                              'Position',[.83 .3 .12 .08],'String','Export',...
                              'Callback','stimPlusTSeriesMovie(''export'');');
return
% /----------------------------------------------------------------------/ %

% /----------------------------------------------------------------------/ %
function playStimMovie(nTimes);
% plays the movie once it's made.
if ieNotDefined('nTimes')   nTimes = 1;   end
stim = get(gcf,'UserData');
movie(stim.movie,nTimes,stim.movieRate);
return
% /----------------------------------------------------------------------/ %

% /----------------------------------------------------------------------/ %
function rebuildStimMovie;
% rebuild the stimulus movie, with user-modified params.
return
% /----------------------------------------------------------------------/ %

% /----------------------------------------------------------------------/ %
function exportStimMovie;
% export the movie to an AVI file.
if isunix
    compression = 'None';
else
    compression = 'Indeo5';
end
stim = get(gcf,'UserData');
[fname, pth] = myUiPutFile(pwd,'*.avi','Name your saved movie...');
savePath = fullfile(pth,fname);
movie2avi(stim.movie,savePath,'FPS',stim.movieRate,'Compression',compression);
fprintf('Saved %s.\n',savePath);
return


text:
function cross_time=MeanCrossingRate(data)
%	Acc feas Computation
%   INPUT:
%   data NO*1
%   OUTPUT:
%   cross_time: 1*1 mean crossing rate

mean_series = data - mean(data);
cross_time = 0;
for i = 2:length(data)
    if mean_series(i-1)*mean_series(i) < 0
        cross_time = cross_time + 1;
    end
end

author: jindongwang

repo: activityrecognition


text:
function [A,B,h,k,phi]=polyEllipsForm2ScalRotTrans(coeffs)
%%POLYELLIPSFORM2SCALROTTRANS Given a 2D ellipse represented as a quadratic
%   polynomial such that f(x,y)=0 is points on the ellipse, this function
%   obtains the parameterization for a parameteric representation of the
%   points in terms of scale factors (A, and B), a rotation angle (phi) and
%   translation coordinates (h, and k). Thus, if t is a parameteric
%   parameter that runs from 0 to 2*pi, x and y coordinates of points on
%   the ellipse are given by
%   x=A*cos(phi)*cos(t)-B*sin(phi)*sin(t)+h;
%   y=A*sin(phi)*cos(t)+B*cos(phi)*sin(t)+k;
%
%INPUTS: coeffs A 3X3 matrix of the coefficients for the bivariate
%               quadratic polynomial. These are arranged such
%               that coeffs(a1,a2) corresponds to the coefficient
%               of an x1^(a1-1)*x2^(a2-1) term. This is the format used in
%               polyValMultiDim.
%
%OUTPUTS: A, B The scalar positive scale factors.
%         h, k The scalar offsets in the x and y coordinates.
%          phi The scalar rotation angle in radians.
%
%The relations here come from inverting Equation 5 in [1]. Using the
%notation of [1], we note that we can write BB=(1/A2-1/B2)*sin(2*phi)
%and also AA-CC==(1/A^2-1/B^2)*cos(2*phi), so we use an inverse transgent to
%find phi. After that, the other terms are found by substituting and
%solving each of the equations. In the final equation for FF, we repalce
%the 1 with a c that we solve for. That is a scale factor (since one can
%multiply coeffs by any constant without changing the result). Given the
%scale factor, we then replace A with A/sqrt(c) and B with B/sqrt(c) to
%undo the scaling.
%
%EXAMPLE:
%To show that this produces correct results, we plot an ellipse using the
%drawEllipse function, which uses a different parameterization. We then
%convert that parameterization into polynomial coefficients using
%quadEllipsForm2Poly. Those coefficients are then transformed into the
%parameters for the scaled, rotated, and translated parametric form of
%this function. Being in parametric form, we plot the ellipse over the
%other one. One can see that they are the same.
% A=inv([4, -1.5;
%       -1.5,1]);
% x0=[2;3]*10;
% c=2.5;
% %Plot the original ellipse
% figure(1)
% clf
% drawEllipse(x0,A,c,'-k','linewidth',4)
% hold on
% %Convert the ellipse into polynomial form.
% coeffs=quadEllipsForm2Poly(A,x0,c)/1000;
% %Convert the polynomial form into the scaled, rotated, translated form.
% [A,B,h,k,phi]=polyEllipsForm2ScalRotTrans(coeffs);
% %This is a parameteric form, so plot the points on the ellipse:
% numPoints=500;
% t=linspace(0,2*pi,numPoints);
% x=A*cos(phi)*cos(t)-B*sin(phi)*sin(t)+h;
% y=A*sin(phi)*cos(t)+B*cos(phi)*sin(t)+k;
% plot(x,y,'-y','linewidth',1)
% legend('Original Ellipse','Parametric Ellipse')
% xlabel('x')
% ylabel('y')
%
%REFERENCES:
%[1] G. B. Hughes and M. Chraibi, "Calculating ellipse overlap areas,"
%    Computing and Visualization in Science, vol. 15, no. 5, pp. 291-301,
%    Oct. 2012.
%
%January 2023 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

AA=coeffs(3,1);
BB=coeffs(2,2);
CC=coeffs(1,3);
DD=coeffs(2,1);
EE=coeffs(1,2);
FF=coeffs(1,1);

phi=atan2(BB,AA-CC)/2;
sinPhi=sin(phi);
cosPhi=cos(phi);
%Using a double angle identity:
cos2Phi=cosPhi^2-sinPhi^2;

A2=2*cos2Phi/(AA-CC+(AA+CC)*cos2Phi);
B2=2*cos2Phi/(CC-AA+(AA+CC)*cos2Phi);

h=(1/2)*(-A2*DD+(B2-A2)*cosPhi*EE*sinPhi+(A2-B2)*DD*sinPhi^2);
k=(1/2)*(-B2*EE+(B2-A2)*cosPhi*DD*sinPhi-(A2-B2)*EE*sinPhi^2);

c=(h*cosPhi+k*sinPhi)^2/A2+(h*sinPhi-k*cosPhi)^2/B2-FF;

%Adjust the scaling of everything.
A2=A2*c;
B2=B2*c;

A=sqrt(A2);
B=sqrt(B2);

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 [img, permInds, flipInds] = stackRotate90(img, axis, varargin)
%STACKROTATE90 Rotate a 3D image by 90 degrees around one image axis
%
%   RES = stackRotate90(IMG, AXIS);
%   IMG is a 3D image (either gray scale or color), and AXIS is the axis
%   number, in XYZ convention: 1-> X-axis, 2->Y-axis, 3->Z-axis.
%   AXIS can also be specified as letter: 'x', 'y' or 'z'.
%
%   RES = stackRotate90(IMG, AXIS, NUMBER);
%   Apply NUMBER rotation around the axis. NUMBER is the number of
%   rotations to apply, between 1 and 3. NUMER can also be negative, in
%   this case the rotation is performed in reverse direction.
%
%   [RES PERM] = stackRotate90(...);
%   [RES PERM FLIP] = stackRotate90(...);
%   Returns also informations about rotation:
%   PERM is the permutation of indices from the first image, such that
%   stackSize(RES) = stackSize(permute(IMG, PERM));
%   FLIP is the set of dimension indices that were flipped after dimension
%   permutation.
%
%   Example
%     img = discreteEllipsoid(1:90, 1:100, 1:110, [50 50 50 30 20 10]);
%     img2 = stackRotate90(img, 'x');
%     figure(1); clf;
%     subplot(121); imshow(rot90(img(:,:,50)));
%     title('rotated slice');
%     subplot(122); imshow(squeeze(img2(50,:,:)));
%     title('slice of rotated image');
%
%   See also
%   imStacks, rotateStack90
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2010-05-18,    using Matlab 7.9.0.529 (R2009b)
% http://www.pfl-cepia.inra.fr/index.php?page=slicer
% Copyright 2010 INRA - Cepia Software Platform.

% parse axis, and check bounds
axis = parseAxisIndex(axis);

% get rotation parameters in xyz ordering
[permInds, flipInds] = getStackRotate90Params(axis, varargin{:});

% check if image is color
colorImage = length(size(img)) > 3;

% convert indices in xyz ordering to ijk ordering
permInds2 = xyz2ijk(permInds([2 1 3]), colorImage);
flipInds2 = xyz2ijk(flipInds, colorImage);

% in the case of color image, adds the channel coordinate after spatial
% coordinates
if colorImage
    permInds2 = [permInds2([1 2]) 3 permInds2(3)];
end

% apply matrix dimension permutation
img = permute(img, permInds2);

% depending on rotation, some dimensions must be fliped
for i=1:length(flipInds2)
    img = flip(img, flipInds2(i));
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function c = butter(c,varargin)

% C = BUTTER(C,[TYPE],CUTOFF,[POLES])
% This function creates and applies a butterworth filter to each waveform
% in the correlation object. The function returns a correlation object.
% Filtering is performed by an underlying call to the filterobject/filtfilt
% routine. Note that filtfilt applies the filter once in each direction to
% minimize phase distortion (effectively a zero-phase filter). This also
% has the effect of doubling the order of the filter. 2 or 4 poles should
% be sufficient for most applications.
%
% EXAMPLES:
%  c = BUTTER(c,[1 5])            band pass filter on 1-5 Hz (4 poles)
%  c = BUTTER(c,'B',[1 5])        same as previous example
%  c = BUTTER(c,'L',5)            low  pass filter below 5 Hz (4 poles)
%  c = BUTTER(c,'H',1)            high pass filter above 5 Hz (4 poles)
%  c = BUTTER(c,...,2)            use 2 poles

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$
% TODO: should check to ensure low cutoff period is >> trace length

% GET INPUTS
if ischar(varargin{1})       % set filter type
    type = varargin{1};
    varargin = varargin(2:end);
else
    type = 'B';
end

cutoff = varargin{1};            % set cutoff frequencies

if length(varargin)>1
   poles = varargin{2};
else
    poles = 4;                  % set number of poles
end

% CHECK FIRST LETTER OF FILTER TYPE
type = upper(type(1));
if (type~='B') && (type~='H') && (type~='L')
    error('Filter type not recognized')
end;

% CHECK NUMBER OF CUTOFF FREQUENCIES
if type=='B'
    if length(cutoff)~=2
        error('Two cutoff frequecies needed for bandpass filter')
    end;
else
    if length(cutoff)~=1
        error('High and lowpass filters require one cutoff frequency')
    end;
end;

% CHECK FREQUENCY RANGE
if cutoff(end) >= get(c,'NYQ')
    error('Frequency cutoffs exceed the Nyquist frequency')
end;

% duration = get(c,'DURATION_EPOCH');
% if 1/cutoff(1) >= duration(1)
%     warning('Frequency cutoff is very low relative to trace length');
% end;
% uration(1)
% 1/cutoff(1)

% APPLY FILTER
f = filterobject(type,cutoff,poles);
c.W = filtfilt(f,c.W);

author: geoscience-community-codes

repo: GISMO


text:
classdef MimSubjectList < MimModel
    methods (Access = protected)
        function value = run(obj)
            database = obj.Callback.getMim().GetImageDatabase();
            [projectNames, projectIds] = database.GetListOfProjects();
            subjectList = {};
            for projectIndex = 1 : numel(projectNames)
                projectName = projectNames{projectIndex};
                projectId = projectIds{projectIndex};
                [names, ids, short_visible_names, num_series, num_patients_combined, patient_id_map] = database.GetListOfPatientNamesAndSeriesCount(projectId, true);
                for nameIndex = 1 : numel(names)
                    subjectId = ids{nameIndex};
                    name = names{nameIndex};
                    parameters = {};
                    parameters.subjectId = subjectId;
                    parameters.projectName = projectName;
                    parameters.projectId = projectId;
                    parameters.subjectName = name;

                    modelId = obj.Callback.buildModelId('MimWSSubject', parameters);
                    subjectList{end + 1} = MimSubjectList.SubjectListEntry(modelId, name, subjectId, projectName, [], []);

                end
            end
            value = subjectList;
        end
    end

    methods (Static, Access = private)
        function subjectListEntry = SubjectListEntry(modelId, label, id, project, uri, insert_date)
            subjectListEntry = struct();
            subjectListEntry.modelId = modelId;
            subjectListEntry.label = label;
            subjectListEntry.ID = id;
            subjectListEntry.project = project;
            subjectListEntry.URI = uri;
            subjectListEntry.insert_date = insert_date;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function  [Ypred, ps] = evaluateHist(model, X)

grids = model.grids;
ntypes = size(model.hist,3);

logps = zeros(size(X,1), ntypes);
for j = 1:size(model.grids,2)
    xp = X(:, j+1);
    xp(xp<grids(1,j))   = grids(1,j);
    xp(xp>grids(end,j)) = grids(end,j);

    [~, ~, ibin]    = histcounts(xp, grids(:,j));

    logps           = logps + log(sq(model.hist(ibin, j, :)));
end

ps = ones(1,ntypes)/ntypes;
for j = 1:10
    L    = bsxfun(@plus, logps, log(ps));
    L    = bsxfun(@minus, L, max(L, [], 2));
    rs   = exp(L) + 1e-5;
    rs   = bsxfun(@rdivide, rs, sum(rs,2));
    ps   = mean(rs,1);
end

Ypred = rs(:,1);

author: cortex-lab


text:
function [gid, fid] = HDF5_GroupOpen(fileobj, location)

% Find data location
if ischar(fileobj)
    fid = H5F.open(fileobj);
else
    fid = fileobj;
end
try
    gid = H5G.open(fid, location);
catch
    gid.double = -1;
end

author: fieldtrip

repo: fieldtrip


text:
% -------------------------------------------------------------------------------------------------
function rect  = get_axis_aligned_rect(region)
%GETAXISALIGNEDRECT computes axis-aligned rect with same area as the rotated one (REGION)
% -------------------------------------------------------------------------------------------------
if size(region,2) == 8
    cx = mean(region(:,1:2:end),2);
    cy = mean(region(:,2:2:end),2);
    x1 = min(region(:,1:2:end),[],2);
    x2 = max(region(:,1:2:end),[],2);
    y1 = min(region(:,2:2:end),[],2);
    y2 = max(region(:,2:2:end),[],2);
    x1y1x2y2 = region(:,1:2) - region(:,3:4);
    x2y2x3y3 = region(:,3:4) - region(:,5:6);
    A1 = sqrt(sum(x1y1x2y2.*x1y1x2y2,2)).* sqrt(sum(x2y2x3y3.*x2y2x3y3,2));
    A2 = (x2 - x1) .* (y2 - y1);
    s = sqrt(A1./A2);
    w = s .* (x2 - x1) + 1;
    h = s .* (y2 - y1) + 1;
    rect = [cx-w/2,cy-h/2,w,h];
else
    rect = region;
end
end

author: foolwood

repo: DCFNet

save_path: github-repos/MATLAB/foolwood-DCFNet