text:
function export_mesh(fileName,meshStruct)

% function export_mesh(fileName,meshStruct)
% ------------------------------------------------------------------------
% This function exports a .MESH file based on the meshStruct which should
% contain the vertices and elements.
%
% Elements can be specified as either:
% meshStruct.elements or meshStruct.E
%
% Vertices can be specified as either:
% meshStruct.vertices or meshStruct.V or meshStruct.nodes
%
% This function currently only supporte tetrahedral and hexahedral element
% exporting.
%
% The element type is either inferred from the size of the element matrix
% or can be set using:
% meshStruct.elementType
% Use meshStruct.elementType='Hexahedral', for 8-noded hexahedral elements
% Use meshStruct.elementType='Tetrahedral', for 4-noded tetrahedral elements
%
% This implementation was created using examples of .MESH files contained
% in: https://github.com/cnr-isti-vclab/HexaLab
%
%
% Change log:
% 2021/07/20 Kevin M. Moerman Created for use with https://www.hexalab.net/
%
% To do:
% - Proper treatment of vertex/element "attributes"
% - Support more .MESH options (surface meshes, FEA types,....)
% - Test for tetrahedral meshes
% - Expand and test for surface meshes
% ------------------------------------------------------------------------

%% Parse input

%Get elements
if isfield(meshStruct,'E')
    E=meshStruct.E;
elseif isfield(meshStruct,'elements')
    E=meshStruct.elements;
end

%Get vertices
if isfield(meshStruct,'V')
    V=meshStruct.V;
elseif isfield(meshStruct,'nodes')
    V=meshStruct.nodes;
elseif isfield(meshStruct,'vertices')
    V=meshStruct.vertices;
end

%Get/set element type
if isfield(meshStruct,'elementType')
    elementType=meshStruct.elementType;
else
    switch size(E,2)
        case 4 % 4 noded tetrahedral elements
            elementType='Tetrahedra';
        case 8 % 8 noded tetrahedral elements
            elementType='Hexahedra';
    end
end

%Replace "FEBio-style" element type designations by .MESH types
if strcmp(elementType,'tet4')
    elementType='Tetrahedra';
elseif strcmp(elementType,'hex8')
    elementType='Hexahedra';
end

%% Open file for writing

fileID = fopen(fileName,'w'); %Open file for writing

%% Write .MESH file version line
% Example:
%
% MeshVersionFormatted 1

fprintf(fileID,'%s \n','MeshVersionFormatted 1');

%% Write dimensionality line
% Example:
%
% Dimension 3

fprintf(fileID,'Dimension %d \n',size(V,2));

%% Write line specifying number of vertices
% Example:
%
% Vertices 12

fprintf(fileID,'Vertices %d \n',size(V,1));

%% Write vertex data field
% Example:
%
% 0 0 0 1
% -1.0000000000000000e+00 -1.0000000000000000e+00 -1.0000000000000000e+00 1

%Prepare data
V_id=ones(size(V,1),1);
V_field=[V V_id(:)];

%Write field
tForm=repmat('%0.16e ',1,size(V,2)); %Replicated text format
tForm=[tForm(1:end-1),' %i \n']; %Text format with attribute and end of line
fprintf(fileID,tForm,V_field');

%% Write line specifying number of elements and element type
% Example:
%
% Hexahedra 50

fprintf(fileID,[elementType,' %d \n'],size(E,1));

%% Write element data field

%Prepare data
E_id=ones(size(E,1),1); %1:1:size(E,1);
E_field=[E E_id(:)];

%Write field
tForm=repmat('%d ',1,size(E,2)); %Replicated text format
tForm=[tForm(1:end-1),' %i \n']; %Text format with attribute and end of line
fprintf(fileID,tForm,E_field');

%% Write "End" line
fprintf(fileID,'%s \n','End');

%% Close file for writing
fclose(fileID);

end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [calscores,params] = parallel_cal_augm(w,scores)
%

if nargin==0
    test_this();
    return;
end

if ~exist('scores','var') || isempty(scores)
    calscores = sprintf(['parallel calibration:',repmat(' %g',1,length(w))],w);
    return;
end

[m,n] = size(scores);
scores = [scores;zeros(1,n)];
wsz = 2*m;

[whead,wtail] = splitvec_fh(wsz,w);
params.get_w0 = @(wfuse) init_w0(wfuse);
params.tail = wtail;

waugm = augmentmatrix_fh(m,0,whead);
calscores = linTrans(waugm,@(w)map_this(w),@(w)transmap_this(w));

    function w0 = init_w0(wfuse)
        scal = wfuse(1:end-1);
        offs = wfuse(end);
        W = [scal*(m+1);((m+1)/m)*offs*ones(m,1)];
        w0 = W(:);
    end

    function y = map_this(w)
        w = reshape(w,m+1,2);
        y = bsxfun(@times,scores,w(:,1));
        y = bsxfun(@plus,y,w(:,2));
    end

    function w = transmap_this(y)
        y = reshape(y,m+1,n);
        w = [sum(y.*scores,2),sum(y,2)];
    end

end

function test_this()

    scores = randn(4,10);
    [sys,params] = parallel_cal_augm([],scores);
    w0 = params.get_w0();
    test_MV2DF(sys,w0);

end

author: nesl


text:
function [result] = oneMatrix(size)
    result =ones(size, 'double', 'gpuArray');
end

author: jiweil

repo: Hierarchical-Neural-Autoencoder


text:
function s = signum (x)
%SIGNUM compute and display the sign of a column vector x
% Example
%   s = signum(x)
% See also: testall

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

s = ones (length (x),1) ;
s (find (x < 0)) = -1 ;     %#ok
disp ('s =') ;
disp (s) ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ output_args ] = bz_LFPPowerDist( LFP,varargin )
% bz_LFPPowerDist(LFP) calculates the power distribution of an LFP signal.
%   NOTE: THIS FUNCTION IS UNDER DEVELOPMENT. Feel free to chip in.

%
%INPUTS
%   LFP     [nt x 1] vector of the LFP signal -or- basename
%   (optional)
%   sf          sampling frequency of the LFP (default 1250Hz)
%   int         restrict analysis to specific time intervals
%   spectype    'wavelet' or 'FFT'
%   frange      [lowf highf]
%   nfreqs      number of frequencies to look at
%   SHOWFIG     true/false (default: false)
%   figfolder   folder to save output figures
%
%% DEV
%datasetfolder = '/Users/dlevenstein/Dropbox/Research/Datasets/BWData/';
%figfolder = '/Users/dlevenstein/Dropbox/Research/Current Projects/misc/PowerDistribution/';
%recname = '20140526_277um';

%load(fullfile(datasetfolder,recname,[recname,'_LFP.mat']))
%load(fullfile(datasetfolder,recname,[recname,'_SleepScore.mat']))
%%
sf = 1250;
%int = StateIntervals.NREMpacket;

%LFP = LFP.CTX;
%LFP = NormToInt(LFP,int,sf,'modZ');

spectype = 'FFT';

%%
warning('This function is still under development. No promises')

frange = [1 128];
nfreqs = 100;

switch spectype
    case 'FFT'
        freqlist = logspace(log10(frange(1)),log10(frange(2)),nfreqs);
        window = 1;
        noverlap = 0.8;
        window = window*sf;
        noverlap = noverlap*sf;
        [spec,freqs,t_FFT] = spectrogram(LFP,window,noverlap,freqlist,sf);
        spec = abs(spec)';

        [~,inintIDX] = RestrictInts(t_FFT',int);
        intspec = spec(inintIDX,:);
    case 'wavelet'
        downsamplefactor = 2;
        sf_down = sf./downsamplefactor;
        intLFP = IsolateEpochs2(downsample(LFP,2),int,0,sf_down);
        ncyc = 5;
        [freqlist,t,intspec] = bz_WaveSpec(intLFP,frange,nfreqs,ncyc,1/sf_down,'log');
        intspec = cat(2,intspec{:});
        intspec = abs(intspec)';
end

%%
numpowerbins = 200;

%minpower = min(intspec(:)); maxpower = max(intspec(:));
switch spectype
    case 'FFT'
        minpower = 0.5;maxpower = 2.85e3;
    case 'wavelet'
        minpower = -2;maxpower = 2;
end

powerbins = linspace(log10(minpower),log10(maxpower),numpowerbins);
[powerdist_mean] = hist(log10(intspec),powerbins);

%%
figure
%subplot(2,2,1)
imagesc(log2(freqlist),powerbins,powerdist_mean)
axis xy
LogScale('x',2)

%%
alldists = {'birnbaumsaunders','exponential','extreme value','gamma',...
    'generalized extreme value','generalized pareto','inverse gaussian',...
    'logistic','loglogistic','lognormal','nakagami','normal','rayleigh',...
    'rician','tlocationscale','weibull'};

%Distirbutions were removed that consistently showed bad fit to cortical
%LFP data during NREM
testdists = {'gamma','loglogistic','lognormal',...
    'rayleigh','weibull'};

showexamples = [];

D = {};PF = {};bestfit={};
for ff = 1:nfreqs
    ff
    if ismember(ff,showexamples)
        [distfits] = allfitdist(intspec(:,ff),'PDF');
    else
        [distfits] = allfitdist(intspec(:,ff));
    end

    %Keep only the distributions tested
    distnames = {distfits(:).DistName};
    keepdists = ismember(distnames,testdists);
    distfits = distfits(keepdists);

    %Check the best and worst-fitting distribution
    bestfit{ff} = distfits(1).DistName;
    worstfit{ff} = distfits(end).DistName;

    %Sort alphabetically
    distnames = {distfits(:).DistName};
    [distnames,sortdist] = sort(distnames);

    D{ff} = distfits(sortdist);
end

%%
AICs = cellfun(@(X) [X(:).AIC],D,'UniformOutput',false);
AICs = cat(1,AICs{:});

ndists = length(distnames);
bestfitmat = zeros(nfreqs,ndists);
for nn = 1:ndists
    bestfitmat(:,nn) = nn.*strcmp(distnames{nn},bestfit);
end
bestfitmat(bestfitmat==0)=nan;
%%
figure
subplot(2,2,1)
    imagesc(log2(freqlist),powerbins,powerdist_mean)
    axis xy
    LogScale('x',2);LogScale('y',10)
    xlabel('f (Hz)')
    ylabel('Power (AU)')
    title([spectype,' Power Distribution'])

subplot(4,2,5)
    plot(log2(freqlist),log10(AICs),'LineWidth',1)

    legend(distnames,'location','southwest')
    axis tight
    LogScale('x',2);%LogScale('y',10)
    xlabel('f (Hz)')
    box off
    ylabel('log(AIC)')

subplot(6,2,11)
        plot(log2(freqlist),bestfitmat,'.','markersize',20)
        set(gca,'YTick',1:ndists)
        set(gca,'YTickLabels',distnames)
    LogScale('x',2);%LogScale('y',10)
    xlabel('f (Hz)')
    box off

    NiceSave('LFPPowerDist',figfolder,recname)

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ remove ] = selectHierarchies( data, NODETYPES )

%% select the hierarchies which meet our requirement

labellist = data.labellist;
kids = data.kids;
remove = 0;

% there should be one floor
ind = strmatch('floor',labellist,'exact');
if(length(ind)~=1)
    remove = 1;
    fprintf('don not have a floor\n');
end

% there should be one bed
%ind = strmatch('bed',labellist,'exact');
%if(length(ind)~=1)
%    remove = 1;
%    fprintf("don't have one bed\n");
%end

% the number of objects should be between 2~10, 8~15 including the walls
objnum = length(labellist);
%if(objnum<14||objnum>34)
if(objnum<8||objnum>15)
    remove = 1;
    fprintf('obj num:%d\n', objnum);
end

% check kids
for i = 1:length(kids)
    k = kids{i};
    if(length(k)>0)
        flag = k(1);
        if(flag~=0&&flag~=NODETYPES.COOCCUR&&flag~=NODETYPES.SUPPORT&&flag~=NODETYPES.SURROUND&&flag~=NODETYPES.ROOM&&flag~=NODETYPES.WALL)
            remove = 1;
            fprintf('node type error\n', objnum);
        end
        if(flag==NODETYPES.COOCCUR||flag==NODETYPES.SUPPORT||flag==NODETYPES.WALL)
            if(length(k)~=3)
                remove = 1;
                fprintf('group/support node error\n', objnum);
            end
        end
        if(flag==NODETYPES.SURROUND)
            if(length(k)~=4)
                remove = 1;
                fprintf('surround node error\n', objnum);
            end
        end
        if(flag==NODETYPES.ROOM)
            if(length(k)~=6)
                remove = 1;
                fprintf('room node error\n', objnum);
            end
        end
    else
        remove = 1;
    end
end

% set the max room size
MAXROOMSIZE = 7;
obblist = data.obblist;
for i = 1:size(obblist,2)
    cent = obblist(1:3,i);
    sizes = obblist(10:12,i);
    ind = find(abs(cent)>MAXROOMSIZE);
    ind2 = find(abs(sizes)>MAXROOMSIZE);
    if(length(ind)>0||length(ind2)>0)
        remove = 1;
        fprintf('exceed maxroom size\n', objnum);
    end
end

end

repo: GRAINS


text:
function [sz] = size(tt,p)
%Mode sizes of the TT-matrix
%   [SZ]=SIZE(TT) Returns the mode sizes of the TT-matrix as a dx2 array of
%   integers
%
%   [SZ]=SIZE(TT,P) Returns the mode-P sizes of the TT-matrix as an array
%   of d integers
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

if (nargin==1) || isempty(p)
    n=tt.n; m=tt.m;
    sz=[n,m];
else
    switch p
        case 1
            sz=tt.n;
        case 2
            sz=tt.m;
        otherwise
            error('tt_matrix/size : illegal value of parameter p.\n');
    end
end

return
end

author: oseledets

repo: TT-Toolbox


text:
function [ know, x ] = p37_sol ( n )

%*****************************************************************************80
%
%% P37_SOL returns the solution for problem 37.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the problem.  This value
%    is only needed for those problems with variable N.
%
%    Output, integer KNOW.
%    If KNOW is 0, then the solution is not known.
%    If KNOW is positive, then the solution is known, and is returned in X.
%
%    Output, real X(N), the solution, if known.
%
  know = 1;

  x = pi * ones ( n, 1 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [done] = optMaxGenTerm(ops,bPop,endPop)
% function [done] = optMaxGenTerm(ops,bPop,endPop)
%
% Returns 1, i.e. terminates the GA, when either the maximal_generation is
% reached or when the optimal function val is found.
%
% ops    - a vector of options [maximum_generation optimal epsilon]
% bPop   - a matrix of best solutions [generation_found solution_string]
% endPop - the current generation of solutions

% Binary and Real-Valued Simulation Evolution for Matlab
% Copyright (C) 1996 C.R. Houck, J.A. Joines, M.G. Kay
%
% C.R. Houck, J.Joines, and M.Kay. A genetic algorithm for function
% optimization: A Matlab implementation. ACM Transactions on Mathmatical
% Software, Submitted 1996.
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 1, or (at your option)
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. A copy of the GNU
% General Public License can be obtained from the
% Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

currentGen = ops(1);
maxGen     = ops(2);
optimal    = ops(3);
epsilon    = ops(4);
fitIndex   = size(endPop,2);
bestSolVal = max(endPop(:,fitIndex));
done       = (currentGen >= maxGen) | ((optimal - bestSolVal) <= epsilon);

author: HuangCongQing


text:
function simp_plot_2d ( p, t )

%*****************************************************************************80
%
%% SIMP_PLOT_2D displays a plot of the triangles that form a mesh in 2D.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, real P(NP,2), the coordinates of a set of nodes in 2D.
%
%    Input, integer T(NT,3), a list of the nodes which make up each triangle in the mesh.
%
  h = trimesh ( t, p(:,1), p(:,2), 0*p(:,1), ...
    'facecolor', 'none', 'edgecolor', 'k' );

  view ( 2 )
  axis equal
  axis off

  return
end


text:
function result = p00_rat_transform ( problem, order )

%*****************************************************************************80
%
%% P00_RAT_TRANSFORM applies a rational transform and Gauss-Legendre rule.
%
%  Discussion:
%
%    To approximate:
%
%      Integral ( alpha <= x < Infinity ) f(x) dx
%
%    Transform:
%
%      u = 1 / ( 1 + x )
%      du = - dx / ( 1 + x )^2
%
%      x = ( 1 - u ) / u
%      dx = - du / u^2
%
%      x = alpha    => u = 1 / ( 1 + alpha )
%      x = Infinity => u = 0
%
%    Transformed integral:
%
%      Integral ( 0 < u <= 1 / ( 1 + alpha ) ) f ( ( 1 - u ) / u ) du / u^2
%
%    We apply a Gauss-Legendre rule here, but we could easily use any rule
%    that avoids evaluation at U = 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur Stroud, Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966,
%    LC: QA299.4G3S7.
%
%  Parameters:
%
%    Input, integer PROBLEM, the index of the problem.
%
%    Input, integer ORDER, the order of the Gauss-Legendre rule
%    to apply.
%
%    Output, real RESULT, the approximate integral.
%
  alpha = p00_alpha ( problem );
%
%  Get the abscissas and weights for Gauss-Legendre quadrature.
%
  [ u, weight ] = legendre_compute ( order );
%
%  Modify the weights from [-1,1] to [0,1/(1+alpha)].
%
  weight(1:order) = weight(1:order) / 2.0 / ( 1.0 + alpha );
%
%  Linear transform of abscissas from [-1,1] to [0,1/(1+alpha)].
%
  u(1:order) = ( ( 1.0 + u(1:order) ) / ( 1.0 + alpha ) ...
               + ( 1.0 - u(1:order) ) * 0.0 )           ...
               / ( 2.0              );
%
%  Define U_RAT = ( 1 - U ) / U.
%
  u_rat(1:order) = ( 1.0 - u(1:order) ) ./ u(1:order);
%
%  Evaluate F ( (1-U)/U ).
%
  f_vec = p00_fun ( problem, order, u_rat );
%
%  The integrand is F ( (1-U)/U ) / U^2
%
  f_vec(1:order) = f_vec(1:order) ./ u(1:order).^2;
%
%  Sum.
%
  result = weight(1:order) * f_vec(1:order)';

  return
end


text:
function rs_imwrite_bands(varargin)
% rs_imwrite - Write an ENVI image with differnt bands with names (Copy Right - Zhe Zhu, Sept. 4th 2018)
% Examples
% rs_imwrite(image,'filename',info,n_names);
% Input
%   image: image matrix
%   filename: File names
%   info: Image properties and cartographic information for ENVI or GeoTIFF file
%   n_names   : Band names 1985:2015 or 'disturb type';

if nargin < 4
    fprintf('Failed! Less Than Three Variables! \n');
    return
elseif nargin == 4
    % Create input variables
    image = varargin{1};
    fname = varargin{2};
    info = varargin{3};
    n_names = varargin{4};

    im_size=size(image);
    im_size(3)=size(image,3);

    % ENVI data format numbers
    d = [1 2 3 4 5 12];

    % Check user input
    if ~ischar(fname)
        error('fname should be a char string');
    end

    cl=class(image);
    switch cl
        case 'uint8'
            t = d(1);
        case 'int16'
            t = d(2);
        case 'int32'
            t = d(3);
        case 'single'
            t = d(4);
        case 'double'
            t = d(5);
        case 'uint16'
            t = d(6);
        otherwise
            error('Data type not recognized');
    end

    % update info
    info.samples = im_size(2);
    info.lines = im_size(1);
    info.bands = im_size(3);
    info.data_type = t;

    try multibandwrite(image,fname,info.interleave);
        fprintf('Write ENVI Images %s',fname);
        fprintf((' . '));
    catch
        fprintf('Images Failed! Not ENVI Variables! \n');
        return
    end

    % Write header file
    fprintf((' . '));
    fid = fopen(strcat(fname,'.hdr'),'w');

    fprintf(fid,'%s \n','ENVI');
    n_field = fieldnames(info);
    n_value = struct2cell(info);

    for i = 1:length(n_field)-1
        n_field{i} = strrep(n_field{i},'_',' ');
        if ischar(n_value{i})
            fprintf(fid,'%s = %s \n',n_field{i},n_value{i});
        else
            fprintf(fid,'%s = %d \n',n_field{i},n_value{i});
        end
    end

    % band names
    fprintf(fid,'band names = {');
    if ischar(n_names)
        fprintf(fid,'%s',n_names);
    else
        for i=1:length(n_names)
            str_n = num2str(n_names(i));
            fprintf(fid,'%s',str_n);
            if i < length(n_names)
                fprintf(fid,',');
            end
        end
    end
    fprintf(fid,'}\n');

        elements={'samples =','lines   =','bands   =',...
        'header offset =','file type =', 'data type =',...
        'interleave =','sensor type =','byte order =',...
        'map info =','coordinate system string =','wavelength units =',...
        'band names ='};

    fprintf(fid,'%s \n','ENVI');

    try
        fprintf(fid,'%s %s \n','description =',info.description);
    catch
        fprintf(fid,'%s %s \n','description =',info.descirption); % to fix the inccorect info made by the old data-preparesion funciton.
    end
    fprintf(fid,'%s %d \n',elements{1,1},im_size(2));
    fprintf(fid,'%s %d \n',elements{1,2},im_size(1));
    fprintf(fid,'%s %d \n',elements{1,3},im_size(3));
    fprintf(fid,'%s %d \n',elements{1,4},info.header_offset);
    fprintf(fid,'%s %s \n',elements{1,5},info.file_type);
    fprintf(fid,'%s %d \n',elements{1,6},t);
    fprintf(fid,'%s %s \n',elements{1,7},info.interleave);
    fprintf(fid,'%s %s \n',elements{1,8},info.sensor_type);
    fprintf(fid,'%s %d \n',elements{1,9},info.byte_order);
    fprintf(fid,'%s %s \n',elements{1,10},info.map_info);
% % %     fprintf(fid,'%s %s \n',elements{1,11},info.coordinate_system_string);
% % %     fprintf(fid,'%s %s \n',elements{1,12},info.wavelength_units);
% % %     fprintf(fid,'%s %s \n',elements{1,13},info.band_names);

    fclose(fid);
    fprintf((' . \n'));
else
    fprintf('Failed! Too Many Variables! \n');
    return
end

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
% Read all images and extract point coordinates.
%
% All information needed are stored and retrieved
% from the function CONFIGDATA

% $Author: svoboda $
% $Revision: 2.6 $
% $Id: im2points.m,v 2.6 2005/05/23 16:26:03 svoboda Exp $
% $State: Exp $

clear all;

% add path to config data
addpath ../../CommonCfgAndIO
% add path for graphical output if needed
addpath ../OutputFunctions

SHOWFIG	  = 0; % show images during point extraction
STEP4STAT = 1; % step for computing average and std images, if 1 then all images taken

% Read configuration from whatever is specified on command-line (via --config=FILENAME)
config = read_configuration();

im.dir = config.paths.img;
im.ext = config.files.imgext;

NoCams = size(config.files.idxcams,2);	% number of cameras

% load image names
for i=1:NoCams,
  seq(i).camId = config.files.idxcams(i);
  if seq(i).camId > -1
	if strfind(config.expname,'oscar')
	  seq(i).data = dir([sprintf(im.dir,seq(i).camId),config.files.imnames,'*.',im.ext]);
	else
	  seq(i).data = dir([sprintf(im.dir,seq(i).camId),sprintf(config.files.imnames,seq(i).camId),im.ext]);
	end
  else
	seq(i).data = dir([im.dir,sprintf(config.files.imnames),im.ext]);
  end
  seq(i).size = size(seq(i).data,1);
  if seq(i).size<4
	error('Not enough images found. Wrong image path or name pattern?');
  end
end

% create an occupancy matrix for image frames
occmat=1; try config.files.maxid; catch occmat=0; end
if occmat,
	NoPoints = config.files.maxid;
	FrameMat = zeros(config.files.maxid,NoCams);
	for i=1:NoCams,
		seq(i).imgidx = zeros(size(1:NoPoints));
		for j=1:size(seq(i).data,1)
			FrameMat(str2num(seq(i).data(j).name(config.files.posid)),i)=j;
		end
	end
else
	NoPoints = min([seq.size]);
	FrameMat = zeros(NoPoints,NoCams);
	for i=1:NoCams,
		FrameMat(:,i) = [1:NoPoints]';
	end
end

% In fact, some frames might be without any calibration point

% Becouse of non-consistent stopping of capturing, the sequences might
% have different number of images, select the minimal value as the right one

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% beginning of the findings

t = cputime;
for i=1:NoCams,
  if ~exist(sprintf(config.files.avIM,seq(i).camId)),
	disp(sprintf('The average image of the camera %d is being computed',seq(i).camId));
	avIM = zeros(size(imread([sprintf(im.dir,seq(i).camId),seq(i).data(1).name])));
	pointIdx = 1:STEP4STAT:seq(i).size;
	for j=pointIdx,
	  IM = imread([sprintf(im.dir,seq(i).camId),seq(i).data(j).name]);
	  avIM = avIM + double(IM);
	end
	avIM = uint8(round(avIM./size(pointIdx,2)));
	imwrite(avIM,sprintf(config.files.avIM,seq(i).camId));
  else	disp('Average file already exists');
  end
end
disp(sprintf('Elapsed time for computation of average images: %4.2f [sec]',cputime-t))
% compute the standard deviations images that will be used for finding LEDs
% if not already computed
t = cputime;
for i=1:NoCams,
  if ~exist(sprintf(config.files.stdIM,seq(i).camId)),
	avIM = double(imread(sprintf(config.files.avIM,seq(i).camId)));
	disp(sprintf('The image of standard deviations of the camera %d is being computed',seq(i).camId));
	stdIM = zeros(size(imread([sprintf(im.dir,seq(i).camId),seq(i).data(1).name])));
	pointIdx = 1:STEP4STAT:seq(i).size;
	for j=pointIdx,
	  IM = imread([sprintf(im.dir,seq(i).camId),seq(i).data(j).name]);
	  stdIM = stdIM + (double(IM)-avIM).^2;
	end
	stdIM = uint8(round(sqrt(stdIM./(size(pointIdx,2)-1))));
	imwrite(stdIM,sprintf(config.files.stdIM,seq(i).camId));
  else
	disp('Image of standard deviations already exists')
  end
end

disp(sprintf('Elapsed time for computation of variance images: %4.2f [sec]',cputime-t))

% find points in the images
Ws    = [];	  % joint image matrix
Res	  = [];	  % resolution of cameras
% UsableFramesIdx = find(sum(FrameMat')>2);
IdMat = ones(NoCams,NoPoints);
% IdMat is very important for Martinec&Pajdla filling [ECCV2002]
% it is a NoCams x NoPoints matrix,
% IdMat(i,j) = 0 -> no j-th point in i-th
% IdMat(i,j) = 1 -> point successfully detected

disp('*********************************************')
disp('Finding points (laser projections) in cameras')
disp(sprintf('Totally %d cameras, %d images for each cam', NoCams, NoPoints'))
disp('*********************************************')
for i=1:NoCams,
  t1 = cputime;
  disp(sprintf('Finding points in camera No: %0.2d',config.files.idxcams(i)))
  Points = [];
  avIM  = imread(sprintf(config.files.avIM,seq(i).camId));
  stdIM	= imread(sprintf(config.files.stdIM,seq(i).camId));
  for j=1:NoPoints,
	  fprintf(1,'\b\b\b\b\b\b %5d',j);
	  idx2data = FrameMat(j,i);
	  if idx2data
		  [pos,err] = getpoint([sprintf(im.dir,seq(i).camId),seq(i).data(idx2data).name], SHOWFIG, config.imgs, avIM, stdIM);
	  else
		  err = 1;
	  end
	if err
	  IdMat(i,j) = 0;
	  Points = [Points, [NaN; NaN; NaN]];
	else
	  Points = [Points, [pos; 1]];
	end
  end
  Ws = [Ws; Points];
  Res= [Res; size(avIM,2), size(avIM,1)];
  t2 = cputime;
  disp(sprintf('\nElapsed time for finding points in one camera: %d minutes %d seconds',floor((t2-t1)/60), round(mod((t2-t1),60))))
  disp(sprintf('%4d points found in camera No: %0.2d',sum(Points(3,:)>0),config.files.idxcams(i)));
end

%%% End of the findings
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strfind(config.expname,'oscar')
  % needs special care for handling projector data
  ProjPoints = load(config.files.projdata,'-ASCII');
  Ws = [Ws; ProjPoints(:,end-1:end)'; ones(size(ProjPoints(:,1)'))];
  IdMat = [IdMat; ones(size(ProjPoints(:,1)'))];
  Res	= [Res; config.imgs.projres];
end

save(config.files.points, 'Ws','-ASCII')
save(config.files.Res, 'Res', '-ASCII')
save(config.files.IdMat, 'IdMat', '-ASCII')

% display the overall statistics
disp('Overall statistics from im2points:  ************************  ')
disp(sprintf('Total number of frames (possible 3D points): %d',NoPoints))
disp(sprintf('Total number of cameras %d', NoCams))
disp('More important statistics: *********************************  ')
disp(sprintf('Detected 3D points:                    %d', sum(sum(IdMat)>0)))
disp(sprintf('Detected 3D points in at least 3 cams: %d', sum(sum(IdMat)>2)))
disp(sprintf('Detected 3D points in ALL cameras:     %d', sum(sum(IdMat)==NoCams)))

author: strawlab

repo: MultiCamSelfCal


text:
function [val1,val2,val3,val4,val5,val6,val7] = dtiGetValFromTensors(dt6, coords, xform, valName, interpMethod)
% Interpolates dt6 tensor field and computes stats at each of the coords
%
%  [val1,val2,val3,val4,val5,val6,val7] = ...
%     dtiGetValFromTensors(dt6, coords, [xform], [valName], [interpMethod])
%
% This  also works for a scalar image in place of the dt6. In that case,
% 'valName' is ignored. (OK, but what does it do in that case? BW)
%
% Inputs:
%  xform: the transform that converts coords to dt6 indices. Default is
%    eye(4) (ie. no xform).  As an example,
%    dtiFiberGroupPropertyWeightedAverage sets the xform to
%    inv(dt.xformToAcpc).  This is the inverse of the transform from coords
%    to acpc, suggesting that the coords in there are acpc and the xform
%    maps from acpc to image space.
%  coords: Nx3 matrix of coords. The returned val will be of length N.  The
%    coords are often in ACPC space while the dt6 indices are in image
%    space.  So xform would move from ACPC to image, typically.
%  interpMethod: 'nearest', 'trilin' (default), 'spline'
%
% The values returned (val1, val2, etc.) differ depending on the string in
% valName.  The current valnaem options are:
%    - 'fa' (fractional anisotropy) (DEFAULT)
%    - 'md' (mean diffusivity)
%    - 'eigvals' (triplet of values for 1st, 2nd and 3rd eigenvalues)
%    - 'shapes' (triplet of values indicating linearity, planarity and
%              spherisity)
%    - 'dt6' (the full tensor in [Dxx Dyy Dzz Dxy Dxz Dyz] format
%    - 'pdd' (principal diffusion direction)
%    - 'linearity'
%    - 'fa md pdd', 'fa md ad rd', 'fa md pdd dt6'
%    - 'fa md ad rd shape'
%
% HISTORY:
% 2005.03.18 RFD (bob@white.stanford.edu) wrote it.
% 2006.08.07 RFD: we no longer set NaNs to 0. If there are missing
%   data, the caller should know about it and deal with as they wish.
% 2012.09 BW - Coments and consistency with other methods
%
% Bob (c) Stanford VISTA Team

%% Parameter initialization
if(~exist('xform','var') || isempty(xform)),  xform = eye(4); end
if(~exist('valName','var') || isempty(valName)),  valName = 'fa'; end
if(~exist('interpMethod','var') || isempty(interpMethod))
    interpMethod = 'trilin';
end
if(size(coords,2)~=3), coords = coords'; end
if(size(coords,2)~=3), error('coords must be an Nx3 array!'); end

% The coordinates are in YYY space. The data are in XXX space. We convert
% the coordinates from YYY space to XXX space. The parameter xfrom is the
% way to do the conversion, and this is sent in. Hence, this routine has no
% way to know what the coordinate frames are. Typically, the coordinates
% are in ACPC and we convert to image. A comment from someone who really
% knows would be nice (BW).
if(~all(all(xform==eye(4))))
    coords = mrAnatXformCoords(xform, coords);
end

%% Data volumetric interpolation
% We use the SPM interpolation interpolating the measurements to the
% sampling resolution of the coordinates.This requires some parameters that
% depend on the method.  Those parameters are set here.  A comment would be
% nice (BW).
switch lower(interpMethod)
    case 'nearest'
        interpParams = [0 0 0 0 0 0];
    case 'trilin'
        interpParams = [1 1 1 0 0 0];
    case 'spline'
        interpParams = [7 7 7 0 0 0];
    otherwise
        error(['Unknown interpMethod "' interpMethod '".']);
end
val_dt6 = zeros(size(coords,1), size(dt6,4));
for ii=1:size(dt6,4)
    bsplineCoefs  = spm_bsplinc(dt6(:,:,:,ii), interpParams);
    val_dt6(:,ii) = spm_bsplins(bsplineCoefs, coords(:,1), coords(:,2), coords(:,3), interpParams);
end
clear dt6 coords;

%% Initialize the returns.  Not all are used.
val1 = []; val2 = []; val3 = []; val4 = []; val5 = []; val6 = [];
if(size(val_dt6,2)~=6)
    % This might be the scalar case?  (BW)
    val1 = val_dt6(:,1);
    return;
end

% The user passes in with spaces and upper/lower variation.  We squeeze out
% the spaces and force to lower here.  That makes the switch easier.
valName = mrvParamFormat(valName);
switch lower(valName)
    case 'dt6'
        val1 = val_dt6(:,1);
        val2 = val_dt6(:,2);
        val3 = val_dt6(:,3);
        val4 = val_dt6(:,4);
        val5 = val_dt6(:,5);
        val6 = val_dt6(:,6);
    case 'eigvals'
        val_dt6 = reshape(val_dt6, [size(val_dt6,1) 1 1 6]);
        [nil, eigVal] = dtiSplitTensor(val_dt6);
        val1 = squeeze(eigVal(:,:,:,1));
        val2 = squeeze(eigVal(:,:,:,2));
        val3 = squeeze(eigVal(:,:,:,3));
    case 'fa'
        val_dt6 = reshape(val_dt6, [size(val_dt6,1) 1 1 6]);
        [nil, eigVal] = dtiSplitTensor(val_dt6);
        val1 = dtiComputeFA(eigVal);
        %val(isnan(val)) = 0;
    case 'md'
        % mean diffusivity: trace/3, where trace is the sum of the diagonal
        % elements (ie. the first three dt6 values)
        val1 = sum(val_dt6(:,1:3),2)./3;
    case {'shapes','linearity'}
        val_dt6 = reshape(val_dt6, [size(val_dt6,1) 1 1 6]);
        [eigVec, eigVal] = dtiSplitTensor(val_dt6);
        [val1, val2, val3] = dtiComputeWestinShapes(eigVal);
    case 'pdd'
        % principal diffusion direction
        val_dt6 = reshape(val_dt6, [size(val_dt6,1) 1 1 6]);
        [eigVec, eigVal] = dtiSplitTensor(val_dt6);
        val1 = squeeze(eigVec(:,:,:,[1 2 3],1)); % Should be [1 3 2]?
        %val(isnan(val)) = 0;
    case 'famdpdd'
        % FA, Mean diffusivity, PDD
        val_dt6 = reshape(val_dt6, [size(val_dt6,1) 1 1 6]);
        [eigVec, eigVal] = dtiSplitTensor(val_dt6);
        val1 = dtiComputeFA(eigVal);
        %val(isnan(val)) = 0;
        val2 = sum(val_dt6(:,1:3),2)./3;
        val3 = squeeze(eigVec(:,:,:,[1 2 3],1)); % Should be [1 3 2]?
        %val3(isnan(val3)) = 0;
    case 'famdadrd'
        % FA, mean diffusivity, axial diffusivity, radial diffusivity
        val_dt6 = reshape(val_dt6, [size(val_dt6,1) 1 1 6]);
        [eigVec, eigVal] = dtiSplitTensor(val_dt6);
        val1 = dtiComputeFA(eigVal);  % FA
        %val(isnan(val)) = 0;
        val2 = sum(val_dt6(:,1:3),2)./3;  % Mean diffusivity
        val3 = squeeze(eigVal(:,:,:,1));  % Axial diffusivity
        val4 = squeeze(eigVal(:,:,:,2)+eigVal(:,:,:,3))./2; % Radial
    case 'famdpdddt6'
        % FA MD PDD DT6
        val_dt6 = reshape(val_dt6, [size(val_dt6,1) 1 1 6]);
        [eigVec, eigVal] = dtiSplitTensor(val_dt6);
        val1 = dtiComputeFA(eigVal);
        %val(isnan(val)) = 0;
        val2 = sum(val_dt6(:,1:3),2)./3;
        val3 = squeeze(eigVec(:,:,:,[1 2 3],1)); % Should be [1 3 2]?
        %val3(isnan(val3)) = 0;
        val4 = val_dt6;
    case 'famdadrdshape'
        % FA MD AD RD shape
        val_dt6 = reshape(val_dt6, [size(val_dt6,1) 1 1 6]);
        [eigVec, eigVal] = dtiSplitTensor(val_dt6);
        val1 = dtiComputeFA(eigVal);
        %val(isnan(val)) = 0;
        val2 = sum(val_dt6(:,1:3),2)./3;
        val3 = squeeze(eigVal(:,:,:,1));
        val4 = squeeze(eigVal(:,:,:,2)+eigVal(:,:,:,3))./2;
        %val5=linearity (cp); val6=planarity (cp); val7=sphericity (cs)
        [val5, val6, val7] = dtiComputeWestinShapes(eigVal);

    otherwise
        error(['Unknown tensor value "' valName '".']);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [yd] = km2yd(km)
% Convert length from kilometers to yards.
% Chad A. Greene 2012
yd = km*1093.613298338;

author: Sable

repo: mcbench-benchmarks


text:
function [ x, y, z, w ] = ld0350 ( )

%*****************************************************************************80
%
%% LD0350 computes the 350 point Lebedev angular grid.
%
%  Modified:
%
%    14 September 2010
%
%  Author:
%
%    Dmitri Laikov
%
%  Reference:
%
%    Vyacheslav Lebedev, Dmitri Laikov,
%    A quadrature formula for the sphere of the 131st
%    algebraic order of accuracy,
%    Russian Academy of Sciences Doklady Mathematics,
%    Volume 59, Number 3, 1999, pages 477-481.
%
%  Parameters:
%
%    Output, real X(N), Y(N), Z(N), W(N), the coordinates
%    and weights of the points.
%
  n = 0;
  x = zeros(350,1);
  y = zeros(350,1);
  z = zeros(350,1);
  w = zeros(350,1);
  a = 0.0;
  b = 0.0;
  v = 0.3006796749453936E-02;
  [ n, x, y, z, w ] = gen_oh ( 1, n, a, b, v, x, y, z, w );
  v = 0.3050627745650771E-02;
  [ n, x, y, z, w ] = gen_oh ( 3, n, a, b, v, x, y, z, w );
  a = 0.7068965463912316;
  v = 0.1621104600288991E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.4794682625712025;
  v = 0.3005701484901752E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.1927533154878019;
  v = 0.2990992529653774E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6930357961327123;
  v = 0.2982170644107595E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.3608302115520091;
  v = 0.2721564237310992E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6498486161496169;
  v = 0.3033513795811141E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.1932945013230339;
  v = 0.3007949555218533E-02;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.3800494919899303;
  v = 0.2881964603055307E-02;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.2899558825499574;
  b = 0.7934537856582316;
  v = 0.2958357626535696E-02;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.9684121455103957E-01;
  b = 0.8280801506686862;
  v = 0.3036020026407088E-02;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.1833434647041659;
  b = 0.9074658265305127;
  v = 0.2832187403926303E-02;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );

  return
end


text:
function FitResultsSave_nii(FitResults,fname_copyheader,folder)
% Save FitResutls with headers copied from nifti file
% FitResultsSave_nii(FitResults,fname_copyheader)
% FitResultsSave_nii(FitResults,fname_copyheader,folder)
%
% Example:
%   FitResultsSave_nii(FitResults,'merged_crop_eddy_moco.nii')
if ~exist('folder','var'), folder = 'FitResults'; end

if length(strfind(folder,filesep)) <= 1 && length(strfind(folder,['.' filesep]))~=1
% Normally users should pass folder variable in the ./myDir format to create folder
% myDir at their current directory. This condition guards against users who pass
% folder without './' expression preceding the folder name (e.g. myDir or /myDir).
% The conditional statement checks for the existence of this case, if present,
% it prepends absolute path of the pwd to the folder name to ensure that
% the folder is created without possibly conflicting with another folder that may
% be saved in the MATLAB/Octave's search path (e.g. Test).

    folder = [pwd filesep folder];
    mkdir(folder);
end

if ~exist(folder,'dir')
% Avoids warnings to amazing users who would like to save somehting to an
% existing directory using absolute path.
    mkdir(folder);
end

for i = 1:length(FitResults.fields)
    map = FitResults.fields{i};
    file = strcat(map,'.nii.gz');
    if ~exist('fname_copyheader','var') || isempty(fname_copyheader)
        save_nii_v2(make_nii(FitResults.(map)),fullfile(folder,file),[],64);
    else
    save_nii_v2(FitResults.(map),fullfile(folder,file),fname_copyheader,64);
    end
end
if isfield(FitResults,'Model')
    FitResults.Model = objProps2struct(FitResults.Model);
end
save(fullfile(folder,'FitResults.mat'),'-struct','FitResults')

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
classdef HandleObj < handle
  properties
    data
  end
end

author: kristinbranson

repo: JAABA


text:
%% sparseSphericalMask
% Below is a demonstration of the features of the |sparseSphericalMask| function

%%
clear; close all; clc;

%% Syntax
% |[subIndMask,patchDataOut]=sparseSphericalMask(rSphere,nSubTri,optPatchOut);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.


text:
% batch run full-clustering on all fish

% f.pushbutton_autoclus_Callback

global VAR;

data_masterdir = GetCurrentDataDir();

% range_fish = [5,6,7];
% M_ClusGroup = [2,2,2,2];
% M_Cluster = [1,1,1,1];
range_fish = 2:2;
M_ClusGroup = 1;
M_Cluster = 1;
M_stim = [1];
% M_fish_set = [1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2];

%%
M_param = 5:5:20%0.3:0.1:0.8;

ParamScores = cell(length(M_param),length(range_fish));
ParamScores_raw = cell(length(M_param),length(range_fish));

for k_param = 1:length(M_param),
    numK2 = M_param(k_param);
    %thres_split = M_param(k_param);
    %setappdata(hfig,'thres_split',thres_split);

    for k_fish = 1:length(range_fish),
        i_fish = range_fish(k_fish);
        disp(i_fish);
        LoadFullFish(hfig,i_fish);
        absIX = getappdata(hfig,'absIX');

        %% partitions for CV
        timelists = getappdata(hfig,'timelists');
        timelists_names = getappdata(hfig,'timelists_names');
        periods = getappdata(hfig,'periods');
        %if length(periods)>1,
            timelistsCV = cell(length(M_stim),2);
            for k_stim = 1:length(M_stim), % :3
                i_stim = M_stim(k_stim);
                TL = timelists{i_stim};
                period = periods(i_stim);
                nrep = size(TL,2)/periods(i_stim); % integer
                n = floor(nrep/2);
                timelistsCV{k_stim,1} = TL(1):TL(n*period);
                timelistsCV{k_stim,2} = TL(1+n*period):TL(2*n*period);
            end
        %end

        %% CV loop: auto-clustering with the partitions
        Score = zeros(length(M_stim),2);

        for k_stim = 1:length(M_stim), % :3
            i_stim = M_stim(k_stim);
            NumClus = zeros(1,2);
            CIX = cell(1,2);
            GIX = cell(1,2);
            for k = 1:2,
                %% Cluster to start auto-clustering
                i_ClusGroup = M_ClusGroup(k_fish);
                i_Cluster = M_Cluster(k_fish);
                ClusGroup = VAR(i_fish).ClusGroup{i_ClusGroup};
                numK = ClusGroup(i_Cluster).numK;
                gIX = ClusGroup(i_Cluster).gIX;
                cIX_abs = ClusGroup(i_Cluster).cIX_abs; % convert absolute index to index used for this dataset
                [~,cIX] = ismember(cIX_abs,absIX);

                % ~UpdateTimeIndex
                tIX = timelistsCV{k_stim,k};
                M_0 = GetTimeIndexedData_Default_Direct(hfig,cIX,tIX,'isAllCells');

                isWkmeans = 1;
                [cIX,gIX] = AutoClustering(cIX,gIX,absIX,i_fish,M_0,isWkmeans,numK2);

                NumClus(k) = length(unique(gIX));
                CIX{k} = cIX;
                GIX{k} = gIX;
            end
            % plot cell-matching figure
            Score(k_stim,1) = HungarianCV(NumClus(1),NumClus(2),CIX{1},CIX{2},GIX{1},GIX{2});
            Score(k_stim,2) = HungarianCV(NumClus(2),NumClus(1),CIX{2},CIX{1},GIX{2},GIX{1});
        end

        ParamScores{k_param,k_fish} = mean(mean(Score));
        ParamScores_raw{k_param,k_fish} = Score;
    end

end

author: xiuyechen

repo: FishExplorer


text:
function cellInfo = Optical_Flow_Bus(varargin)
% OPTICAL_FLOW_BUS returns a cell array containing bus object information
%
% Optional Input: 'false' will suppress a call to Simulink.Bus.cellToObject
%                 when the MATLAB file is executed.
% The order of bus element attributes is as follows:
%   ElementName, Dimensions, DataType, SampleTime, Complexity, SamplingMode, DimensionsMode, Min, Max, DocUnits, Description

suppressObject = false;
if nargin == 1 && islogical(varargin{1}) && varargin{1} == false
    suppressObject = true;
elseif nargin > 1
    error('Invalid input argument(s) encountered');
end

cellInfo = { ...
  { ...
    'Optical_Flow_Bus', ...
    '', ...
    '', ...
    'Auto', ...
    '-1', {...
{'timestamp', 1, 'uint32', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('ms'), ''}; ...
{'vx', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('m/s'), ''}; ...
{'vy', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('m/s'), ''}; ...
{'quality', 1, 'uint8', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'reserved1', 1, 'uint8', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'reserved2', 1, 'uint16', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
    } ...
  } ...
}';

if ~suppressObject
    % Create bus objects in the MATLAB base workspace
    Simulink.Bus.cellToObject(cellInfo)
end

author: Firmament-Autopilot

repo: FMT-Model


text:
function response = Craby(n,prev)

if n>4 && n<9 || n==1 || n>18
    response = 'cooperate';
else
    response = 'defect';
end


text:
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Compile the Simulink model through Matlab command line.

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Compiling the Simulink model...')

pause(1.5)

try
    positionControlBalancing([], [], [], 'compile')
    positionControlBalancing([], [], [], 'term')

catch ME

    errorMessages = ME;
end

clc

disp('Compilation done.')

% warning about Simulink timing
warning('The model will anyways start with FEW SECONDS OF DELAY after pressing the ''Start Model'' button.')

author: robotology

repo: whole-body-controllers


text:
function [ y, m, d, h, n ] = minute_borrow_common ( y, m, d, h, n )

%*****************************************************************************80
%
%% MINUTE_BORROW_COMMON "borrows" an hour of minutes in a Common date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, D, H, N, the year,
%    month, day, hour and minute representing a date.  On input, N
%    might be negative.
%    On output, H should have decreased by one, and N gone up by 60.
%
  while ( n < 0 )

    n = n + 60;
    h = h - 1;

    [ y, m, d, h ] = hour_borrow_common ( y, m, d, h );

  end

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ============================================
% patch_light.m
%
% Goal: surface rendering with lighting
%       need to open figure or subplot before use
%
% Li Shen
% 02/21/2006 - create

function h = patch_light(vertices,faces)

% different color options
fcolor = [0.55 0.5 0.55];
fcolor = [0.3 0.2 0.2]*2;
fcolor = [0.6 0.4 0.4];
fcolor = [0.5 0.2 0.2];
fcolor = [0.5 0.5 0.5];
fcolor = [0.6 0.3 0.3];
h = patch('faces', faces, 'vertices', vertices,'FaceColor', fcolor, 'EdgeAlpha', 0, 'FaceLighting', 'phong', ...
      'AmbientStrength', 0.9, 'SpecularColorReflectance', 0.5, 'DiffuseStrength',0.5,'SpecularExponent',5);

c = [1 1 1]*0.5;
light('Position',[ 1  0 0],'Style','infinite', 'Color', c);
light('Position',[-1  0 0],'Style','infinite', 'Color', c);
light('Position',[ 0  1 0],'Style','infinite', 'Color', c);
light('Position',[ 0 -1 0],'Style','infinite', 'Color', c);
lighting gouraud;

axis image; box on; view(3);

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function model = modelTieParam(model, paramsList)

% MODELTIEPARAM Tie parameters of a model together.
% FORMAT
% DESC groups of parameters of a model to be seen as one parameter
% during optimisation of the model.
% ARG model : the model for which parameters are being tied
% together.
% ARG paramsList : indices of parameteres to group together. The
% indices are provided in a cell array. Each cell in the array
% contains a vector of indices of parameters that should be
% considered as one parameter. Each group of parameters in each
% cell should obviously be mutually exclusive.
% Alternatively each element of the cell array can be a string which
% is interpreted as a regular expression of names of parameters
% (as returned by modelExtractParam) to be tied.
% RETURN model : the model with the parameters grouped together.
%
% SEEALSO : modelExtractParam, modeExpandParam, modelLogLikeGradients
%
% COPYRIGHT : Neil D. Lawrence, 2003, 2006, 2008
%
% MODIFICATIONS: Antti Honkela, 2009

% MLTOOLS

if ~isfield(model, 'paramGroups')
  if isfield(model, 'nParams')
    model.paramGroups = speye(model.nParams);
  elseif isfield(model, 'numParams')
    model.paramGroups = speye(model.numParams);
  else
    error('Model does not list number of parameters.');
  end
end
colToDelete = [];
try,
  [params, names] = modelExtractParam(model);
catch
  try,
    [params, names] = kernExtractParam(model);
  catch
    names = {};
  end
end

for i = 1:length(paramsList)
  if ischar(paramsList{i}),
    paramIndices=find(~cellfun('isempty', regexp(names, paramsList{i})));
    if isempty(paramIndices),
      warning(['No matches for parameter tie spec: ', paramsList{i}])
    end
  else
    paramIndices=sort(paramsList{i});
  end
  if any(paramIndices(1) == colToDelete)
    error('Parameter is already being tied')
  end

  for j = 2:length(paramIndices)
    model.paramGroups(paramIndices(j), paramIndices(1)) = 1;
    if any(paramIndices(j) == colToDelete)
      error('Parameter has already been tied')
    end
    colToDelete = [colToDelete paramIndices(j)];
  end
end

model.paramGroups(:, colToDelete) = [];
if isfield(model, 'nParams')
  % Update to the new number of parameters.
  model.nParams = size(model.paramGroups, 2);
elseif isfield(model, 'numParams')
  model.numParams = size(model.paramGroups, 2);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [featureLexiconNameList, ...
          featureLexiconFileNameList, ...
          featureLexiconAnimalTypeList] = ...
  getFeatureLexiconListsFromXML(jaabarootdir)

% Figure out the name of the XML file holding the list of the possible
% feature lexicons
if isdeployed,
  filename = deployedRelative2Global('params/featureConfigList.xml');
else
  filename = 'featureConfigList.xml';
  if nargin >= 1,
    filename = fullfile(jaabarootdir,'perframe','params',filename);
  end

end
% Read in the XML file
featureLexiconListThing = ReadXMLParams(filename);
% break out that information into three cell arrays, each of length equal
% to the number of feature lexicons
%   featureLexiconNameList{i} holds the name of the ith lexicon, e.g.
%   'flies', 'larvae', 'vivek_mice'.
%   featureLexiconFileNameList{i} holds the file name of the XML file
%     describing the ith lexicon
%   featureLexiconAnimalTypeList{i} holds the name of the animal that the
%     ith lexicon is designed to be used with.  E.g. 'fly', or 'mouse', or
%     'larva'.  This should be singular, not plural.
featureLexiconNameList=fieldnames(featureLexiconListThing);
nFeatureLexicons=length(featureLexiconNameList);
featureLexiconFileNameList=cell(nFeatureLexicons,1);
featureLexiconAnimalTypeList=cell(nFeatureLexicons,1);
for i=1:nFeatureLexicons
  featureLexiconName=featureLexiconNameList{i};
  featureLexiconFileNameList{i}=featureLexiconListThing.(featureLexiconName).file;
  if nargin >= 1,
    featureLexiconFileNameList{i} = fullfile(jaabarootdir,'perframe',featureLexiconFileNameList{i});
  end
  featureLexiconAnimalTypeList{i}=featureLexiconListThing.(featureLexiconName).animal;
end

end


text:
function sum = cfg_example_sum
% Example script that creates an cfg_exbranch to sum two numbers. The
% inputs are entered as vector, the output is just a single
% number. This function differs from cfg_example_add2 (except from names)
% only in the specification of input1.num.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_example_sum.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

%% Input Item
input1         = cfg_entry; % This is the generic data entry item
input1.name    = 'Input a'; % The displayed name
input1.tag     = 'a';       % The name appearing in the harvested job structure. This name must be unique among all items in the val field of the superior node
input1.strtype = 'e';       % No restriction on what type of data is entered. This could be used to restrict input to real numbers, integers ...
input1.num     = [1 Inf];   % Number of inputs required (2D-array with exactly one row and two column)
input1.help    = {'This is the input vector.','The elements will be added together.'}; % help text displayed

%% Executable Branch
sum      = cfg_exbranch;       % This is the branch that has information about how to run this module
sum.name = 'sum';             % The display name
sum.tag  = 'cfg_example_sum'; % The name appearing in the harvested job structure. This name must be unique among all items in the val field of the superior node
sum.val  = {input1};    % The items that belong to this branch. All items must be filled before this branch can run or produce virtual outputs
sum.prog = @cfg_example_run_sum;  % A function handle that will be called with the harvested job to run the computation
sum.vout = @cfg_example_vout_sum; % A function handle that will be called with the harvested job to determine virtual outputs
sum.help = {'Add two numbers.'};

%% Local Functions
% The cfg_example_vout_sum function can go here, it is not useful outside
% the batch environment.
function vout = cfg_example_vout_sum(job)
% Determine what outputs will be present if this job is run. In this case,
% the structure of the inputs is fixed, and the output is always a single
% number. Note that input items may not be numbers, they can also be
% dependencies.

vout = cfg_dep;                        % The dependency object
vout.sname      = 'sum(a)';       % Displayed dependency name
vout.src_output = substruct('()',{1}); % The output subscript reference. This could be any reference into the output variable created during computation

author: spm


text:
function planC = divideStructureAntPostLeftRight(structNum, planC)
% function divideStructureAntPostLeftRight(structNum, planC)
%
% APA, 08/23/2012

if ~exist('planC','var')
    global planC
end

global stateS

indexS = planC{end};

scanNum = getStructureAssociatedScan(structNum,planC);

newStructNumT = length(planC{indexS.structures}) + 1;
newStructNumL = length(planC{indexS.structures}) + 2;
newStructNumB = length(planC{indexS.structures}) + 3;
newStructNumR = length(planC{indexS.structures}) + 4;

newStructTS = newCERRStructure(scanNum, planC, newStructNumT);
newStructLS = newCERRStructure(scanNum, planC, newStructNumL);
newStructBS = newCERRStructure(scanNum, planC, newStructNumB);
newStructRS = newCERRStructure(scanNum, planC, newStructNumR);

for slcNum = 1:length(planC{indexS.scan}(scanNum).scanInfo)

    % Calculate centroid for this slice
    [rasterSegments, planC, isError]    = getRasterSegments(structNum,planC);
    rasterIndices = find(rasterSegments(:,6) == slcNum);

    zValue = planC{indexS.scan}(scanNum).scanInfo(slcNum).zValue;

    if ~isempty(rasterIndices)

        maskM = rasterToMask(rasterSegments(rasterIndices,:), scanNum, planC);

        %Get r,c,s list of all points in mask.
        [rV,cV] = find(maskM);

        %Take the mean of all points... unweighted as this is a mask.
        rowCOM = mean(rV);
        colCOM = mean(cV);
        slcCOM = 1; % dummy value

        %Convert from rcs, to xyz coordinates.
        [xc,yc] = mtoxyz(rowCOM, colCOM, slcCOM, scanNum, planC, 'uniform');

        for segNum = 1:length(planC{indexS.structures}(structNum).contour(slcNum).segments)
            xV = planC{indexS.structures}(structNum).contour(slcNum).segments(segNum).points(:,1);
            yV = planC{indexS.structures}(structNum).contour(slcNum).segments(segNum).points(:,2);
            [xyTV,xyLV,xyBV,xyRV] = dividePolygon(xV,yV,xc,yc);
            newStructTS.contour(slcNum).segments(segNum).points = [xyTV xyTV(:,1).^0*zValue];
            newStructLS.contour(slcNum).segments(segNum).points = [xyLV  xyLV(:,1).^0*zValue];
            newStructBS.contour(slcNum).segments(segNum).points = [xyBV  xyBV(:,1).^0*zValue];
            newStructRS.contour(slcNum).segments(segNum).points = [xyRV  xyRV(:,1).^0*zValue];
        end

    else
        newStructTS.contour(slcNum).segments(1).points = [];
        newStructLS.contour(slcNum).segments(1).points = [];
        newStructBS.contour(slcNum).segments(1).points = [];
        newStructRS.contour(slcNum).segments(1).points = [];
    end

end

stateS.structsChanged = 1;

newStructLS.structureName = [planC{indexS.structures}(structNum).structureName '_RIGHT'];
newStructRS.structureName = [planC{indexS.structures}(structNum).structureName '_LEFT'];
patientPosition = planC{indexS.scan}.scanInfo(1).DICOMHeaders.PatientPosition;
if isequal(patientPosition,'HFP') || isequal(patientPosition,'FFP')
    newStructTS.structureName = [planC{indexS.structures}(structNum).structureName '_POSTERIOR'];
    newStructBS.structureName = [planC{indexS.structures}(structNum).structureName '_ANTERIOR'];
else
    newStructTS.structureName = [planC{indexS.structures}(structNum).structureName '_ANTERIOR'];
    newStructBS.structureName = [planC{indexS.structures}(structNum).structureName '_POSTERIOR'];
end
planC{indexS.structures} = dissimilarInsert(planC{indexS.structures}, newStructTS, newStructNumT);
planC{indexS.structures} = dissimilarInsert(planC{indexS.structures}, newStructLS, newStructNumL);
planC{indexS.structures} = dissimilarInsert(planC{indexS.structures}, newStructBS, newStructNumB);
planC{indexS.structures} = dissimilarInsert(planC{indexS.structures}, newStructRS, newStructNumR);

planC = getRasterSegs(planC, [newStructNumT newStructNumL newStructNumB newStructNumR]);
planC = updateStructureMatrices(planC, newStructNumT);
planC = updateStructureMatrices(planC, newStructNumL);
planC = updateStructureMatrices(planC, newStructNumB);
planC = updateStructureMatrices(planC, newStructNumR);

% Refresh View
if ~isempty(stateS) && isfield(stateS,'handle') && isfield(stateS.handle,'CERRSliceViewer') && isnumeric(stateS.handle.CERRSliceViewer)
    stateS.structsChanged = 1;
    CERRRefresh
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% eeg_getica() - get ICA component activation. Recompute if necessary.
%
% >> mergelocs = eeg_getica(EEG, comp);
%
% Inputs:
%     EEG     - EEGLAB dataset structure
%     comp    - component index
%
% Output:
%     icaact  - ICA component activity
%
% Author: Arnaud Delorme, 2006

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

function icaact = eeg_getica(EEG, comp)

  if nargin < 1
    help eeg_getica;
    return;
  end;
  if nargin < 2
    comp = 1:size(EEG.icaweights,1);
  end;

  if ~isempty(EEG.icaact)
    icaact = EEG.icaact(comp,:,:);
  else
    disp('Recomputing ICA activations');
    if isempty(EEG.icachansind)
        EEG.icachansind = 1:EEG.nbchan;
        disp('Channels indices are assumed to be in regular order and arranged accordingly');
    end
    icaact = (EEG.icaweights(comp,:)*EEG.icasphere)*reshape(EEG.data(EEG.icachansind,:,:), length(EEG.icachansind), EEG.trials*EEG.pnts);
    icaact = reshape( icaact, size(icaact,1), EEG.pnts, EEG.trials);
  end;

author: PatternRecognition

repo: OpenBMI


text:
% SIMSAMPLE Sample from SIM kernel

% KERN

colordef white
randn('seed', 1e8)
rand('seed', 1e8)

numSamp = 100;
t = linspace(0, 5, numSamp)';
kern = kernCreate(t, {'multi', 'rbf', 'sim', 'sim', 'sim'});
for i = 1:length(kern.comp)
  kern.comp{i}.inverseWidth = 1/(0.75*0.75);
end
kern.comp{2}.decay = 5;
kern.comp{2}.variance = 25;
kern.comp{3}.decay = 1;
kern.comp{3}.variance = 1;
kern.comp{4}.decay = 0.5;
kern.comp{4}.variance = 0.25;

params = kernExtractParam(kern);
kern = kernExpandParam(kern, params);

K = kernCompute(kern, t);

imagesc(K, [-1.1 1.1]);
handle = [];
fontSize = 24;
set(gca, 'fontname', 'times')
set(gca, 'fontsize', fontSize)
handle = [handle; text('Interpreter', 'latex', 'string', '$$f(t)$$', 'position', ...
            [50 410], 'fontsize', fontSize, 'horizontalalignment', 'center')];
handle = [handle; text('Interpreter', 'latex', 'string', '$$x_1(t)$$', 'position', ...
            [150 410], 'fontsize', fontSize, 'horizontalalignment', 'center')];
handle = [handle; text('Interpreter', 'latex', 'string', '$$x_2(t)$$', 'position', ...
            [250 410], 'fontsize', fontSize, 'horizontalalignment', 'center')];
handle = [handle; text('Interpreter', 'latex', 'string', '$$x_3(t)$$', 'position', ...
            [350 410], 'fontsize', fontSize, 'horizontalalignment', 'center')];
set(handle, 'horizontalalignment', 'right')
handle = [handle; text('Interpreter', 'latex', 'string', '$$f(t)$$', 'position', ...
            [-10 50], 'fontsize', fontSize, 'horizontalalignment', 'center')];
handle = [handle; text('Interpreter', 'latex', 'string', '$$x_1(t)$$', 'position', ...
            [-10 150], 'fontsize', fontSize, 'horizontalalignment', 'center')];
handle = [handle; text('Interpreter', 'latex', 'string', '$$x_2(t)$$', 'position', ...
            [-10 250], 'fontsize', fontSize, 'horizontalalignment', 'center')];
handle = [handle; text('Interpreter', 'latex', 'string', '$$x_3(t)$$', 'position', ...
            [-10 350], 'fontsize', fontSize, 'horizontalalignment', 'center')];
axis off
colorbar
if exist('printDiagram') && printDiagram
  printPlot('simTestKernelImage', '../tex/diagrams', '../html');
end

numSamp2 = 150;
y = gsamp(zeros(size(K), 1), 0.5*(real(K) + real(K')), numSamp2)';
t = t - 0;
counter = 0;
for i = 1:numSamp2
  x = reshape(y(:, i), numSamp, length(kern.comp));
  if any(x(:, 1)<0)
    continue
  end
  figure
  counter = counter + 1;
  a = plot(t, x(:, 1), 'k-');
  hold on
  a = [a plot(t, x(:, 2), 'r-')];
  a = [a plot(t, x(:, 3), 'g-')];
  a = [a plot(t, x(:, 4), 'b-')];
  set(a, 'linewidth', 2);
  set(gca, 'fontname', 'arial');
  set(gca, 'fontsize', 14);
  if exist('printDiagram') && printDiagram
    printPlot(['simSample'  num2str(counter)], '../tex/diagrams', '../html');
  end
end


text:
function [tuples,gain]=kCard2DAssign(C,k,maximize)
%%KCARD2DASSIGN Solve the two-dimensional k-cardinality assignment problem
%          with a rectangular cost matrix C. The problem being solved can
%          be formulated as minimize (or maximize)
%          \sum_{i=1}^{numRow}\sum_{j=1}^{numCol}C_{i,j}*x_{i,j}
%          subject to
%          \sum_{j=1}^{numCol}x_{i,j}<=1 for all i
%          \sum_{i=1}^{numRow}x_{i,j}<=1 for all j
%          \sum_{i=1}^{numRow}\sum_{j=1}^{numCol}x_{i,j}=k
%          x_{i,j}=0 or 1.
%
%INPUTS: C A numRowXnumCol cost matrix that does not contain any NaNs and
%          where the largest finite element minus the smallest element is a
%          finite quantity (does not overflow) when performing minimization
%          and where the smallest finite element minus the largest
%          element is finite when performing maximization. Forbidden
%          assignments can be given costs of +Inf for minimization and -Inf
%          for maximization.
%        k The integer number of assignments to make.
%          k<=min(numRow,numCol).
% maximize If true, the minimization problem is transformed into a
%          maximization problem. The default if this parameter is omitted
%          or an empty matrix is passed is false.
%
%OUTPUTS: tuples A 2XnumRow set of assignment values. This is ordered
%                [rowIndex;columnIndex]. If no feasible solution exists,
%                then an empty matrix will be returned.
%          gain This is the value of the cost. This is the sum of the
%               values in C corresponding to the tuples.
%
%As noted in [1], the k-cardinality 2D assignment problem can be
%transformed into a standard rectangular 2D assignment problem. This
%function implements that transformation and calls assign2D.
%
%EXAMPLE:
% C=[7,   51,  52,  87;
%    50,  12,   0,  64;
%    27,  77,   0,  18;
%    62,   0,   3,   8];
% k=4;
% [tuples4,gain4]=kCard2DAssign(C,4)
% [tuples3,gain3]=kCard2DAssign(C,3)
% [tuples2,gain2]=kCard2DAssign(C,2)
% [tuples1,gain1]=kCard2DAssign(C,1)
%One will see optimal gains of 25, 7, 0, and 0.
%
%REFERENCES:
%[1] A. Volgenant, "Solving the k-cardinality assignment problem by
%    transformation," European Journal of Operational Research, vol. 157,
%    no. 2, pp. 322-331, 1 Sep. 2004.
%
%June 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if(nargin<3||isempty(maximize))
        maximize=false;
    end

    numRow=size(C,1);
    numCol=size(C,2);

    if(k==0)
       tuples=[];
       gain=0;
       return;
    end

    if(k<0||k~=fix(k))
        error('k is invalid.')
    end

    if(k>numRow||k>numCol)
        %The problem is not feasible.
        tuples=[];
        gain=-1;
        return
    end

    COrig=C;

    if(maximize==true)
        CDelta=max(C(:));

        %If C is all negative, do not shift.
        if(CDelta<0)
            CDelta=0;
        end

        C=-C+CDelta;
    else
        CDelta=min(C(:));

        %If C is all positive, do not shift.
        if(CDelta>0)
            CDelta=0;
        end

        C=C-CDelta;
    end

    %The minimum value in C is now 0 or something positive. We shall shift
    %the values to guarantee that 0 is less than every value in C.
    %Find the maximum finite value in C.
    CMax=max(max(C(isfinite(C))));

    %We shall now offset everything by a small fraction of the maximum
    %value. This is so that we guarantee that values will be assigned
    %to the zero padded rows that shall be added to C.
    COffset=CMax*1e-10;
    C=C+COffset;

    %Perform 2D assignment on the augmented matrix. An algorithm that does
    %not explicitely construct that assignment matrix is also possible.
    [~,row4Col]=assign2D([C;zeros(numCol-k,numCol)],false);

    if(isempty(row4Col))
        %If the problem is infeasible.
        tuples=[];
        gain=-1;
        return
    end

    tuples=zeros(2,k);
    curTuple=1;
    gain=0;
    for curCol=1:numCol

        curRow=row4Col(curCol);

        if(curRow<=numRow)
            gain=gain+COrig(curRow,curCol);
            tuples(:,curTuple)=[curRow;curCol];
            curTuple=curTuple+1;

           if(curTuple>k)
               break;
           end
        end
    end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function a = mag(a)
%MAG          Absolute value for non-intervals, same as abs (for completeness)
%
%  c = mag(a);
%
%Replaces abss. Thanks to Arnold for pointing to better name.
%

% written  01/04/05     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  a = abs(a);

repo: OpenMAS


text:
function kern = velotransKernParamInit(kern)

% VELOTRANSKERNPARAMINIT VELOTRANS kernel parameter initialisation.
% This kernel allows you to take any other kernel in the tool box and
% translate with respect to the input space with respect to time. The translation allows a moving field such as clouds moving over a land surface.
%
% SEEALSO : cmpndKernCreate
%
% FORMAT
% DESC initialises the velocity translate
%  kernel structure with some default parameters.
% ARG kern : the kernel structure which requires initialisation.
% RETURN kern : the kernel structure with the default parameters placed in.
%
% SEEALSO : kernCreate, kernParamInit
%
% COPYRIGHT : Neil D. Lawrence, 2011

% KERN

kern = cmpndKernParamInit(kern);

% Add in the translation parameters
kern.nParams = kern.nParams + kern.inputDimension-1;
kern.velocity = zeros(1, kern.inputDimension-1);


text:
function mySubplot(figureHandle, subplotWidth, subplotHeight, subplotPos, img, imgTitle, newMap)
% MYSUBPLOT creates a matrix of subplots, each with a custom colormap
    changeMap = sprintf('colormap %s', newMap);
    figure(figureHandle)
    subplot(subplotWidth, subplotHeight, subplotPos)
	imagesc(img)
    eval(changeMap);
    freezeColors
	pbaspect([size(img,2),size(img,1),1]);
	title(imgTitle)
    axis off;
end

author: bertinetto

repo: staple

save_path: github-repos/MATLAB/bertinetto-staple


text:
function compareMetaboliteFormulae(modelT)
% Prints out a tab delimited file with the abbreviations, reconstruction
% metabolite formluae, and group contribution metabolite formulae.
%
% USAGE:
%
%    compareMetaboliteFormulae(modelT)
%
% INPUT:
%    modelT:    output of `setupThermoModel`

fid=fopen('metaboliteFormulae.txt','w');

[nMet,nRxn]=size(modelT.S);

for m=1:nMet
    fprintf(fid,'%20s%20s\t%20s\n',modelT.mets{m},modelT.metFormulas{m},modelT.met(m).formulaMarvin);
end
fclose(fid);

author: opencobra

repo: cobratoolbox


text:
function shiftTSeries(vw,shift,scanList)
%
% function shiftTSeries(view,[shift],[scanList])
%
% Shifts the time series (circular shift), to advance the phase of the response.
% This is useful for approximately compensating for the hemodynamic delay before
% averaging across scans, especially for CW and CCW retinotopies.
%
% scanList: vector of scan numbers to flip. Default: prompt user.
% shift: number of frames to shift. Default: 1.
%
% If you change this function make parallel changes in:
%    flipTSeries
%
% djh, 6/2001
if ieNotDefined('scanList'),  scanList = er_selectScans(vw); end
if ieNotDefined('shift'),     shift = 1;                       end

waitHandle = mrvWaitbar(0,'Shifting tSeries.  Please wait...');

for iScan=1:length(scanList)

    scan = scanList(iScan);
    tSeriesFull = [];
    dimNum = 0;
    for slice=sliceList(vw,scan);
        tSeries = loadtSeries(vw,scan,slice);
        tSeries = circularShift(tSeries,0,-shift);
        dimNum = numel(size(tSeries));
        tSeriesFull = cat(dimNum + 1, tSeriesFull, tSeries); %Combine together
        mrvWaitbar(iScan/length(scanList))
    end

    if dimNum == 3
        tSeriesFull = reshape(tSeriesFull,[1,2,4,3]);
    end %if

    savetSeries(tSeriesFull,vw,scan);
end

close(waitHandle);

if length(scanList)==1
    fprintf('Shifted tSeries for scan %i by %i frames.\n',scanList,shift);
else
    fprintf('Shifted tSeries for scans %s by %i frames.\n',...
        num2str(scanList),shift);
end

return


text:
function [energyinband_array,sqi_array,hr_array] = ComputeMVM(ecg,ann_ecg,Fs,segment_size,normalize)

% [energyinband_array,sqi_array] = ComputeMVM(ecg,ann,Fs,segment_size,normalize)
%   OVERVIEW:   This function returns the QRS MVM measured in non overlapping
%   windows of size 'segment_size' in minutes. It also returns the
%   signal quality and mean heart rate for each window.
%
%   INPUT:      MANDATORY:
%               ecg             : a single row of ECG data in samples. Must be sampled at 1 KHz.
%
%               ann             : annotation for the ecg passed as input. The ann variable
%                                 is a struct with Q,R,S fiduical points in sample number as fields
%                                 for the struct. i.e. ann.R, ann.Q, ...
%
%               Fs              : sampling frequency for the ecg signal (Hz)
%
%               segment_size        : the length of non overlapping windows
%                                 being analyzed in minutes. Default size
%                                 is 5 minutes
%
%               normalize       : boolean value (0 or 1). When set to 1
%                               will normalize each analysis window by the median R
%                               amplitude. Default is 1.
%
%
%   OUTPUT:
%               energyinband_array     : the energy for QRS complex morphological variability measured in
%                                        the every 2-7 beats region for each analysis window.
%
%               sqi_array        : the signal quality index measured for each analysis window.
%
%               hr_array          : the median heart rate measured for analysis
%                                 window
%
%   REF:
%   The script is based on the algorithm given in the following paper,
%   Liu Y, Syed Z, Scirica BM, Morrow DA, Guttag JV, Stultz CM.
%   ECG morphological variability in beat space for risk stratification after acute coronary syndrome.
%   J Am Heart Assoc. 2014;3(3):e000981. Published 2014 Jun 24. doi:10.1161/JAHA.114.000981
%	REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Written by Ismail Sadiq
%	COPYRIGHT (C) 2019
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information. The license may be found in
%       the Documents folder of the Physionet-Cardiovascular-Signal-Toolbox.
%

% divide into 5 minute segments
% initialize variables
segment_size_samples = Fs*60*segment_size; % size of analysis window in samples
Nsegments = floor(length(ecg)/segment_size_samples); % number of segments analyzed
energyinband_array = NaN(1, Nsegments); % array for storing the MVM energy measured in each analysis window
sqi_array = NaN(1, Nsegments); % array of storing the signal quality for each window analyzed
hr_array = NaN(1,Nsegments); % array for storing the heart rate estimate for each window analyzed
lengthflag = 0;

sqi_threshold = 0.7;

% convert ecg to row vector
if (size(ecg,1) > size(ecg,2))
   ecg = ecg';
end

if (Nsegments == 0)   % length of data < segment_size_samples, analyze available data
    Nsegments = 1; lengthflag = 1;
end

for segmentidx = 1:Nsegments

    try

        % isolate segment
        if(lengthflag)
            % if the length of signal is < segment_size in samples,
            % determine annotations till the end of the record.
            fiveminsegment = ecg(segmentidx*segment_size_samples+1-segment_size_samples:end);
            ann.R = ann_ecg.R(ann_ecg.R < length(ecg));
            ann.R = ann.R(ann.R > segmentidx*segment_size_samples+1-segment_size_samples) - (segmentidx-1)*segment_size_samples;% + 1;
            ann.Q = ann_ecg.Q(ann_ecg.Q < length(ecg));
            ann.Q = ann.Q(ann.Q > segmentidx*segment_size_samples+1-segment_size_samples) - (segmentidx-1)*segment_size_samples;% + 1;
            ann.S = ann_ecg.S(ann_ecg.S < length(ecg));
            ann.S = ann.S(ann.S > segmentidx*segment_size_samples+1-segment_size_samples) - (segmentidx-1)*segment_size_samples;% + 1;
            ann.QRSon = ann_ecg.QRSon(ann_ecg.QRSon < length(ecg));
            ann.QRSon = ann.QRSon(ann.QRSon > segmentidx*segment_size_samples+1-segment_size_samples) - (segmentidx-1)*segment_size_samples;
            ann.QRSoff = ann_ecg.QRSoff(ann_ecg.QRSoff < length(ecg));
            ann.QRSoff = ann.QRSoff(ann.QRSoff > segmentidx*segment_size_samples+1-segment_size_samples) - (segmentidx-1)*segment_size_samples;
            lengthflag = 0;
        else
            % get ecg signal and corresponding annotations in the current
            % window
            fiveminsegment = ecg(segmentidx*segment_size_samples+1-segment_size_samples:segmentidx*segment_size_samples);
            ann = ann_ecg;
            ann.R = ann_ecg.R(ann_ecg.R < segmentidx*segment_size_samples);
            ann.R = ann.R(ann.R > segmentidx*segment_size_samples+1-segment_size_samples) - (segmentidx-1)*segment_size_samples;% + 1;
            ann.Q = ann_ecg.Q(ann_ecg.Q < segmentidx*segment_size_samples);
            ann.Q = ann.Q(ann.Q > segmentidx*segment_size_samples+1-segment_size_samples) - (segmentidx-1)*segment_size_samples;% + 1;
            ann.S = ann_ecg.S(ann_ecg.S < segmentidx*segment_size_samples);
            ann.S = ann.S(ann.S > segmentidx*segment_size_samples+1-segment_size_samples) - (segmentidx-1)*segment_size_samples;% + 1;
            ann.QRSon = ann_ecg.QRSon(ann_ecg.QRSon < segmentidx*segment_size_samples);
            ann.QRSon = ann.QRSon(ann.QRSon > segmentidx*segment_size_samples+1-segment_size_samples) - (segmentidx-1)*segment_size_samples;
            ann.QRSoff = ann_ecg.QRSoff(ann_ecg.QRSoff < segmentidx*segment_size_samples);
            ann.QRSoff = ann.QRSoff(ann.QRSoff > segmentidx*segment_size_samples+1-segment_size_samples) - (segmentidx-1)*segment_size_samples;
            %figure(2); plot(fiveminsegment); hold on; % Check annotations after loading
            %scatter(ann.R(~isnan(ann.R)), fiveminsegment(ann.R(~isnan(ann.R)))); hold off;
        end

        ecg_mv = fiveminsegment; clear fiveminsegment;

        % remove any nan values from the annotations
        R = round(ann.R);
        sum(isnan(R)); R = R(~isnan(R)); sum(isnan(R));
        Q = round(ann.Q);
        sum(isnan(Q)); Q = Q(~isnan(Q)); sum(isnan(Q));
        S = round(ann.S);
        sum(isnan(S)); S = S(~isnan(S)); sum(isnan(S));

        % compute the median Q,S and R amplitude, used to determine
        % normalization factor
        median_Q_amp = nanmedian(ecg_mv(1,Q));
        median_R_amp = nanmedian(ecg_mv(1,R));
        median_S_amp = nanmedian(ecg_mv(1,S));

        % compute normalization factor
        if (median_S_amp <= median_Q_amp || isnan(median_Q_amp))
            norm_factor = median_R_amp - median_S_amp;
        else
            if (~isempty(median_Q_amp))
                norm_factor = median_R_amp - median_Q_amp;
            else
                % if unable to determine, set to 1
                disp('unable to compute norm factor, setting to 1')
                norm_factor = 1;
            end
        end

        % normalize signal
        if (normalize)
            ecg_mv = ecg_mv / norm_factor;
        end

        % generate second set of annotations
        refqrs = ann.R;
        % mean impute missing val
        ecg_mv(isnan(ecg_mv)) = nanmean(ecg_mv);
        testqrs = wqrsm(ecg_mv, Fs);
        thres = 0.1; margin = 0; windowlen = 60*segment_size;
        % Determine signal quality
        refqrs = refqrs./Fs; testqrs = testqrs./Fs; % convert to time (s)
        [current_win_sqi,Se,PPV,Nb] = run_sqi(refqrs,testqrs,thres,margin,windowlen,Fs);
        sqi_array(segmentidx) = current_win_sqi;
        if (current_win_sqi > sqi_threshold)   % may need to be set higher

            % perform morphological variability on QRS

            % determine qrs onset and qrs offset annotations over 5 minute segment
            qrson = ann.QRSon;
            qrsoff = ann.QRSoff;
            if (qrsoff(1) < qrson(1))
                startidx = find(qrsoff < qrson(1));
                qrsoff = qrsoff(startidx+1:end);
                qrson = qrson(1:length(qrsoff));
            else
                qrson = qrson(1:length(qrsoff));
            end

            % make sure sufficient annotations to estimate median S
            % amplitude
            if (-1*median_S_amp > median_R_amp && (length(S) > 200))
                alignmentpoint = 'S';
            else
                alignmentpoint = 'R';
            end

            % re-compute R peak
            radj = zeros(1,length(qrson)); % r-adjusted, used in ectopic beat removal
            % setting max value between qrson and qrsoff as r-pk
            for qrsonidx = 1:length(qrson)
                %[pks,locs] = findpeaks(ecg(qrson(qrsonidx):qrsoff(qrsonidx)));
                [pks,locs] = findpeaks(ecg_mv(qrson(qrsonidx):qrsoff(qrsonidx)));
                maxpk = max(pks); maxloc = locs((maxpk == pks));
                if (isempty(maxloc))
                    radj(qrsonidx) = qrson(qrsonidx);
                else
                    radj(qrsonidx) = qrson(qrsonidx)-1+maxloc(1);   % 1 in case of duplicates
                end
            end
            % compute alignemnt point
            % find optimal alignment points max r-pk or s-valley
            alignpoint = zeros(1,length(qrson)); % used in qrs complex isolation
            for qrsonidx = 1:length(qrson)

                if (strcmp(alignmentpoint, 'R'))
                    [pks,locs] = findpeaks((ecg_mv(qrson(qrsonidx):qrsoff(qrsonidx))));
                else
                    [pks,locs] = findpeaks(-1*(ecg_mv(qrson(qrsonidx):qrsoff(qrsonidx))));
                end

                if (isempty(pks))   % if no positive peaks define mid-point between qrson and qrsoff as alignment point
                    locs = round(mean(qrson(qrsonidx), qrsoff(qrsonidx)));
                    pks = ecg_mv(locs);
                    alignpoint(qrsonidx) = locs; continue;
                end
                maxpk = max(pks); maxloc = locs((maxpk == pks));
                alignpoint(qrsonidx) = qrson(qrsonidx)-1+maxloc(1);   % 1 in case of duplicates
            end

            % detect pvc beats
            signal = ecg_mv; th = 0.1;
            [qrs_time pvc_output] = detectpvc2(signal',Fs,th);
            qrs_time = qrs_time(logical(pvc_output));

            % only use 5 minute signals with complete annotations to maintane validity
            % of fourier analysis (may want to replace with interpolation techhnique)
            if (sum(isnan(qrson)) || sum(isnan(qrsoff)))
                disp(['annotations incomplete ' num2str(fileidx)]);
                continue;
            end

            % ectopic beat removal
            qrsonectopicfree = qrson(1);
            rectopicfree = radj(1);
            qrsoffectopicfree = qrsoff(1);
            alignpointsectopicfree = alignpoint(1);
            for ridx = 2:length(radj)-1   % start from idx 2 and goto length(r)-1
                if (ridx <= length(radj)-40)
                    %meanRR = mean(diff(radj(ridx-1:ridx+38)));   % compute mean over 40 beat segment
                    meanRR = median(diff(radj(ridx-1:ridx+38)));   % compute median over 40 beat segment
                else
                    %meanRR = mean(diff(radj(end-39:end)));   % compute mean over last 40 beat segment
                    meanRR = median(diff(radj(end-39:end)));   % compute median over last 40 beat segment
                end
                preRR = radj(ridx) - radj(ridx-1);    % pre-RR interval
                postRR = radj(ridx+1) - radj(ridx);   % post-RR interval
                % if the pre-RR and post-RR vary from the mean RR by more than
                % 20% consider as ectopic
                if ((((meanRR-preRR) > 0.2*meanRR) && ((postRR-meanRR) > 0.2*meanRR)) || (((meanRR-preRR) > 0.2*meanRR) && ((meanRR-postRR) > 0.2*meanRR))) % && ((meanRR-postRR) > 0.2*meanRR)))  % 326, 367
                    continue;
                else % include the QRS complex, may also want to include t-wave onset/offset/amplitude/QT-interval
                    qrsonectopicfree = [qrsonectopicfree qrson(ridx)];
                    rectopicfree = [rectopicfree radj(ridx)];
                    qrsoffectopicfree = [qrsoffectopicfree qrsoff(ridx)];
                    alignpointsectopicfree = [alignpointsectopicfree alignpoint(ridx)];
                end
            end

            % remove ectopic beats
            [C,ia,ib] = intersect(rectopicfree, qrs_time);
            rectopicfree(ia) = [];
            qrsonectopicfree(ia) = [];
            qrsoffectopicfree(ia) = [];
            alignpointsectopicfree(ia) = [];

            % calculate NTWDseries and TWDseries
            % Isolate complexes
            [Complexes] = IsolateQRSComplexes(ecg_mv, qrsoffectopicfree, alignpointsectopicfree, qrsonectopicfree);

            % NTWD and TWD
            [NTWDseries TWDseries] = SquaredDiffComputation_beatbybeat(Complexes);

            % compute spectrum
            nfft = 256; % standardize the length of the fft to 256
            pxx = pwelch(NTWDseries,[],[],nfft);
            %             pxx = pwelch(TWDseries,[],[],nfft);
            N = nfft;
            stepsize = 0.5/length(pxx);

            dl = 2; du = 7;
            indexeverytwobeats = round(N/dl, 0);    % corresponding to every 2 beats
            indexeverysevenbeats = floor(N/du); % corresponding to every 7 beats

            % compute energy in every 2-7 beats region
            energyinband = trapz((indexeverysevenbeats:indexeverytwobeats) * stepsize, pxx(indexeverysevenbeats:indexeverytwobeats));
            energyinband_array(segmentidx) = energyinband; disp(['MVM for current window is ' num2str(energyinband)]);

            % hr est
            hr_array(segmentidx) = 60/(nanmedian(diff(alignpointsectopicfree))/Fs);
        else
            disp('unclean segment')
        end

    catch
        disp('error in segment')
        energyinband_array(segmentidx) = -inf;
        sqi_array(segmentidx) = -inf;
        hr_array(segmentidx) = -inf;
    end

end

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
% script to clean the workspace after a branch change
warning('off');
GUI_Edit_Settings.closeGUI
clear all
clear Logger Prj_Settings Go_Settings Settings_Interface Command_Interpreter Command_Settings Core_UI GUI_Edit_Settings Core Remote_Resource_Manager Parallel_Manager Go_Slave Go_Wait_Bar Zernike
clear read_asc.m
% delete last_settings.ini
warning('off');

author: goGPS-Project


text:
function nfgSetupMRDiffusion(phantomDir)
%Setup mrDiffusion files for BlueMatter NFG tests.
%
%   nfgSetupMRDiffusion(phantomDir)
%
%
% AUTHORS:
%   2009.08.05 : AJS wrote it
%
% NOTES:
%   * We are automatically dividing the NFG b-value by 1000, but there
%   should be a more appropriate dynamic way to do this.
%

% Directories
dtDir = nfgGetName('dtDir',phantomDir);
% Input Files
nfgGradFile = nfgGetName('nfgGradFile',phantomDir);
noisyImg = nfgGetName('noisyImg',phantomDir);
brainMaskFile = nfgGetName('brainMaskFile',phantomDir);
% Output Files
bvalsFile = nfgGetName('bvalsFile',phantomDir);
bvecsFile = nfgGetName('bvecsFile',phantomDir);

% XXX b-value factor to adjust into our units
bValFactor = 1000;

% Convert NFG grad file to ours
disp(' '); disp('Converting NFG gradient file to mrDiffusion format ...');
grad = load(nfgGradFile,'-ascii');
bvals = grad(:,4)/bValFactor;
bvecs = grad(:,1:3)';
fid = fopen(bvalsFile,'wt');
fprintf(fid, '%1.3f ', bvals);
fclose(fid);
fid = fopen(bvecsFile,'wt');
fprintf(fid, '%1.4f ', bvecs(1,:)); fprintf(fid, '\n');
fprintf(fid, '%1.4f ', bvecs(2,:)); fprintf(fid, '\n');
fprintf(fid, '%1.4f ', bvecs(3,:));
fclose(fid);

% Do tensor fitting
numBootstraps=300;
disp(' '); disp(['Tensor fitting with ' num2str(numBootstraps) ' bootstraps ...']);
dtiRawFitTensorMex(noisyImg, bvecsFile, bvalsFile, dtDir, numBootstraps);

% Fix automatically generated brain mask
disp(' '); disp('Fixing automatically generated brain mask ...');
ni = niftiRead(brainMaskFile);
ni.data(:) = 1;
writeFileNifti(ni);

% Need to produce pdf image now
disp(' '); disp('Creating ConTrack PDF file ...');
mtrCreateConTrackOptionsFromROIs(0,1,dtDir);

return;

author: vistalab


text:
function [ output_args ] = vVecA2vVecB( vxA, vyA, pxA, pxB )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

end

author: TUMFTM


text:
function [ res_map ] = siam_eval( net_match, patches, feat_groups )
%DSIAM_EVAL This function evaluates the score of each positions on the
%search map
% Input:
%       net_match   - matching net
%       patches     - feature groups of the target
%       feat_groups - feature groups of the search image patch
%       sw_sz       - size of the search window [width hight]
% Output:
%       res_map     - response map
% By Xin Li, 2017-12

    for feat_i =1 : numel(feat_groups)

        [x_h, x_w, ~, num_scale] = size(feat_groups{feat_i});
        [z_h, z_w,~, num_scale_1] = size(patches);
        num_scale = max([num_scale,num_scale_1]);
        z_bh = floor(z_h/2); z_bw = floor(z_w/2);
        res_tmp = zeros(x_h, x_w, num_scale,'like',feat_groups{feat_i});

        net_match.eval({'w_feat', patches,'x_feat', feat_groups{feat_i}});
        res_tmp(z_bh+1:z_bh+x_h-z_h+1, z_bw+1:z_bw+x_w-z_w+1,:) = ...
            squeeze( net_match.vars(net_match.getVarIndex('m_score')).value);
        res_ = reshape(gather(res_tmp),x_h,x_w,[]);

        res_map = res_/max(res_(:));

    end

end

author: XinLi-zn

repo: TADT

save_path: github-repos/MATLAB/XinLi-zn-TADT


text:
function [stFiles] = ex_GetAllFiles(sBaseDir)
% function [stFiles] = ex_GetAllFiles(sBaseDir)
% ---------------------------------------------
% Collects all filenames form a given directory and its subdirectories
%
% Input parameter:
%   sBaseDir    Base directory for filesearch. Function looks into this directory and its subdirectories
%
% Output parameter:
%   stFiles     Vector of structure containing filenames. Filenames are stored in .name element
%
% Danijel Schorlemmer
% October 7, 2002

% Init return values
stFiles = [];

% Get the list of the base directory
files = dir(sBaseDir);

% Iterate through this list store the filenames
for nCnt = 1:(length(files))
  % If element is a directory, analyze it
  if files(nCnt).isdir == 1
    % If it doesn't point upwards, recurse into the subdirectory
    if ~strcmp(files(nCnt).name, '.')  &&  ~strcmp(files(nCnt).name, '..')
      sDir = fullfile(fileparts(sBaseDir), files(nCnt).name);
      stFiles = [stFiles; GetAllFiles(sDir)];
    end
  else
    stFiles = [stFiles; struct('name', fullfile(fileparts(sBaseDir), files(nCnt).name))];
  end
end

author: CelsoReyes


text:
function [cst,cstJac] = autoGen_cst_footVel(q1p,q2p,q4p,q5p,q1m,q2m,q4m,q5m,dq1p,dq2p,dq4p,dq5p,dq1m,dq2m,dq4m,dq5m,l1,l2,l4,l5)
%AUTOGEN_CST_FOOTVEL
%    [CST,CSTJAC] = AUTOGEN_CST_FOOTVEL(Q1P,Q2P,Q4P,Q5P,Q1M,Q2M,Q4M,Q5M,DQ1P,DQ2P,DQ4P,DQ5P,DQ1M,DQ2M,DQ4M,DQ5M,L1,L2,L4,L5)

%    This function was generated by the Symbolic Math Toolbox version 6.3.
%    25-Oct-2015 18:36:52

t2 = sin(q1p);
t3 = sin(q2p);
t4 = sin(q4p);
t5 = sin(q5p);
t6 = sin(q1m);
t7 = sin(q2m);
t8 = sin(q4m);
t9 = sin(q5m);
cst = [dq1p.*l1.*t2+dq2p.*l2.*t3-dq4p.*l4.*t4-dq5p.*l5.*t5;-dq1m.*l1.*t6-dq2m.*l2.*t7+dq4m.*l4.*t8+dq5m.*l5.*t9];
if nargout > 1
    cstJac = reshape([0.0,0.0,dq1p.*l1.*cos(q1p),0.0,dq2p.*l2.*cos(q2p),0.0,0.0,0.0,-dq4p.*l4.*cos(q4p),0.0,-dq5p.*l5.*cos(q5p),0.0,l1.*t2,0.0,l2.*t3,0.0,0.0,0.0,-l4.*t4,0.0,-l5.*t5,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,-dq1m.*l1.*cos(q1m),0.0,-dq2m.*l2.*cos(q2m),0.0,0.0,0.0,dq4m.*l4.*cos(q4m),0.0,dq5m.*l5.*cos(q5m),0.0,-l1.*t6,0.0,-l2.*t7,0.0,0.0,0.0,l4.*t8,0.0,l5.*t9,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],[2,32]);
end

author: MatthewPeterKelly

repo: OptimTraj