text:
function varargout = meshVertexClustering(vertices, faces, spacing, varargin)
%MESHVERTEXCLUSTERING Simplifies a mesh using vertex clustering.
%
%   [V2, F2] = meshVertexClustering(V, F, SPACING)
%   [V2, F2] = meshVertexClustering(MESH, SPACING)
%   MESH2 = meshVertexClustering(...)
%
%   Simplifies a mesh using vertex clustering. Input mesh is specified
%   either by a pair V, F containing the vertex coordinates and the faces
%   informations, or by a structure with fields 'vertices' and 'faces'.
%
%   The SPACING input defines the size of the grid. It can be either a
%   scalar (uniform grid) or a 1-by-3 row vector.
%
%   The output is specified either in two outputs, or in a structure with
%   fields 'vertices' and 'faces'.
%
%   Example
%     [x, y, z]  = meshgrid(1:100, 1:100, 1:100);
%     img = hypot3(x-51.12, y-52.23, z-53.34);
%     [faces, vertices] = isosurface(img, 45);
%     [v2, f2] = meshVertexClustering(vertices, faces, 10);
%     figure; axis equal; axis([0 100 0 100 0 100]);
%     drawMesh(v2, f2);
%
%   See also
%     reducepatch, smoothMesh

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2019-01-28, using Matlab 9.5.0.944444 (R2018b)
% Copyright 2019-2022 INRA - Cepia Software Platform

%% Initialisation

if isstruct(vertices)
    if nargin > 2
        varargin = [{spacing} varargin(:)];
    end
    spacing = faces;

    mesh = vertices;
    vertices = mesh.vertices;
    faces = mesh.faces;
end

% ensure input mesh is a triangulation
faces = triangulateFaces(faces);

% ensure spacing is a 1-by-3 array
if isscalar(spacing)
    spacing = [spacing spacing spacing];
end

% extract grid origin
origin = [0 0 0];
if ~isempty(varargin)
    origin = varargin{1};
end

%% Apply grid simplification

% identify the vertices belonging to the same grid
[v2, I, J] = unique(round(bsxfun(@rdivide, bsxfun(@minus, vertices, origin), spacing)), 'rows');

%% compute reduced vertex coordinates

% compute coordinates of new vertices
for iVertex = 1:length(I)
    gridVertices = vertices(J == iVertex, :);
    v2(iVertex, :) = mean(gridVertices, 1);
end

%% Compute new faces

% create empty array
faces2 = zeros(0, 3);

% iterate over old faces, and keep only faces whose vertices belong to
% different cell grids
nFaces = size(faces, 1);
for iFace = 1:nFaces
    % current face
    face = faces(iFace, :);

    % equivalent face with new vertices
    face2 = J(face)';

    % some vertices may belong to same cell, so we need to adjust
    % processing
    nInds = length(unique(face2));
    if nInds == 3
        % vertices belong to three different cells -> create a new face

        % keep smaller vertex at first position
        [tmp, indMin] = min(face2); %#ok<ASGLU>
        face2 = circshift(face2, [1-indMin 0]);

        % append the new face to the array
        faces2 = [faces2 ; face2]; %#ok<AGROW>
    end
end

% remove duplicate faces
faces2 = unique(faces2, 'rows');

if nargout == 1
    varargout{1} = struct('vertices', v2, 'faces', faces2);
else
    varargout = {v2, faces2};
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function header = sofa_get_header(sofa)
%SOFA_GET_HEADER returns the header of a SOFA file or struct
%
%   Usage: header = sofa_get_header(sofa)
%
%   Input parameters:
%       sofa    - impulse response data set (SOFA struct/file)
%
%   Output parameters:
%       header  - SOFA header
%
%   SOFA_GET_HEADER(sofa) returns the header of the given SOFA file or struct.
%   For the struct the SOFA file has to loaded before with SOFAload().
%   For a description of the SOFA file format see: https://sofaconventions.org
%
%   See also: sofa_get_data, sofa_is_file, get_ir, SOFAload

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 1;
nargmax = 1;
narginchk(nargmin,nargmax)

%% ===== Computation ====================================================
% Get only the metadata of the SOFA data set
if sofa_is_file(sofa)
    header = SOFAload(sofa,'nodata');
else
    header = sofa;
    if isfield(sofa.Data,'IR')
        header.Data = rmfield(sofa.Data,'IR');
    end
end

author: sfstoolbox

repo: sfs-matlab


text:
function [cl] = gal2cl(gal)
% Convert volume from US liquid gallons to centiliters.
% Chad Greene 2012
cl = gal*378.5411784;

author: Sable

repo: mcbench-benchmarks


text:
%% Parameters
Parameters_struct.CenterFrequency = 5e9; % 5 GHz
Parameters_struct.Bandwidth = 20e6; % 20 MHz
Parameters_struct.Ts = 1/Parameters_struct.Bandwidth; % 50 ns
%% Load Given data
load('Long_preamble_slot_Frequency'); % [1x64]
load('data_Payload_1'); % [1x48]
load('data_Payload_2'); % [1x48]
Parameters_struct.Long_preamble_slot_Frequency = Long_preamble_slot_Frequency;
Parameters_struct.data_Payload_1 = data_Payload_1;
Parameters_struct.data_Payload_2 = data_Payload_2;

author: MeowLucian


text:
function vw = buildGrayCoords(vw, coordsPth, keepAll, classPth, numGrayLayers)
% vw = buildGrayCoords(vw, [coordsPth='Gray/coords.mat'], [keepAll=0], [classPth], [numGrayLayers]):
%
% Build a coords.mat file for a gray view. This was split
% off from getGrayCoords when I found the design was getting
% unwieldy -- even when a coords file had been built, there were
% still frequent checks on external files, and I was asked to
% rebuild this nearly every time. I consolidated all of the looking
% for external files into this function -- anything required for
% a gray view that's not already contained in a session directory
% should be loaded up here, and saved to coords.mat. This should
% happen a maximum of once per session/segmentation.
%
% INPUTS:
%	vw: gray or volume view [default: initialize a hidden volume view,
%	which gets converted to gray as the coords are built.]
%
%	coordsPth: path to the coords file to save. [default: Gray/coords.mat']
%
%	keepAll: flag to keep all gray matter nodes, rather than just the
%	intersection with the inplane data. [default: 0, don't do this.]
%	The advantage to keeping all nodes is that you can bring in data from
%	other sessions, and do analyses involving the whole gray matter graph.
%	The disadvantage is it can take more disk space / memory, and you can't
%	map back from the volume to the inplane (though this should be
%	fixable?)
%
%	classPth: optional argument to specify which class file to use
%
%
% ras 08/04.
% ras 06/09: allows user to provide the file paths externally, rather than
% being prompted.
global HOMEDIR

if notDefined('keepAll'),	keepAll=0;			end
if notDefined('classPth'),	classPth = {};		end
if notDefined('vw'),
    % convert from a volume view.
    vw = initHiddenVolume;
    vw = viewSet(vw, 'viewType','Gray');
    vw = viewSet(vw, 'subdir', 'Gray');
end
if notDefined('coordsPth')
    coordsPth = fullfile(HOMEDIR, 'Gray', 'coords.mat');
end

if ~isfield(vw,'leftPath'),  vw = viewSet(vw,'leftpath','(none assigned)');    end
if ~isfield(vw,'rightPath'), vw = viewSet(vw,'rightPath','(none assigned)');   end

% if external file paths are provided, attach them to the view, so that
% they'll be the default files accessed when building the
if isempty(classPth)
    startDir = getAnatomyPath();
    filterspec = {'*.nii.gz;*.nii', 'NIFTI class file'; '*.*', 'All files'};
    classPth = getPathStrDialog(startDir, 'Select class file', filterspec);
end

if iscell(classPth), classPth = classPth{1}; end

% Check class type. Should be nifti.
classType = mrGrayCheckClassType(classPth);

% Since we moved to NIFTI, left and right class files are the same
leftClassFile   = classPth;
rightClassFile  = classPth;
vw.leftPath     = classPth;
vw.rightPath    = classPth;

%% load the gray coords from the class/gray files
% Load the gray nodes and edges
[allLeftNodes,  allLeftEdges,  leftPath]  = loadGrayNodes('left', classPth, numGrayLayers);
[allRightNodes, allRightEdges, rightPath] = loadGrayNodes('right',leftPath, numGrayLayers);

% Load the inplane coords (loaded temporarily into vw.coords)
vw = getVolCoords(vw);
ipCoords = viewGet(vw, 'coords');

% Here, we have the option to include all the nodes in the gray view.
% However, we can't map back to the inplane any more.
if keepAll
    disp(['Warning: You are keeping all the gray nodes ' ...
        '(not computing an intersection with the inplane). ' ...
        'This can be useful in some circumstances. But it  ' ...
        'also means that you can never map things back to ' ...
        'the INPLANE view']);

    leftNodes = allLeftNodes;
    leftEdges = allLeftEdges;
    if ~isempty(allLeftNodes)
        leftCoords = uint16(allLeftNodes([2 1 3],:));
    else
        leftCoords = [];
    end
    keepLeft = 1:length(leftNodes); %#ok<*NASGU>

    rightNodes = allRightNodes;
    rightEdges = allRightEdges;
    if ~isempty(allRightNodes)
        rightCoords = uint16(allRightNodes([2 1 3],:));
    else
        rightCoords = [];
    end
    keepRight = 1:length(rightNodes);
else

    waitHandle = mrvWaitbar(0,'Computing gray coordinates.  Please wait...');
    if ~isempty(allLeftNodes)
        % Find gray nodes that are in the inplanes. Note: nodes are
        % (x,y,z) not (y,x,z), unlike everything else in mrLoadRet.
        [leftCoords,~,keepLeft] = intersectCols(uint16(ipCoords),uint16(allLeftNodes([2 1 3],:)));

        % Keep nodes that intersect with the inplanes
        [leftNodes,leftEdges,nLeftNotReached] = keepNodes(allLeftNodes,allLeftEdges,double(keepLeft));
    else
        leftNodes=[];
        leftEdges=[];
        leftCoords=[];
        keepLeft = [];
    end

    % Repeat for right hemisphere
    mrvWaitbar(1/2)
    if ~isempty(allRightNodes)
        [rightCoords, ~, keepRight] = ...
            intersectCols(uint16(ipCoords), uint16(allRightNodes([2 1 3],:)));
        [rightNodes, rightEdges, nRightNotReached] = ...
            keepNodes(allRightNodes, allRightEdges, keepRight);
    else
        rightNodes=[];
        rightEdges=[];
        rightCoords=[];
        keepRight = [];
    end
    close(waitHandle);
end % End check on whether we intersect with inplane

% Concatenate coords from the two hemispheres
coords = [double(leftCoords) double(rightCoords)];

% Concantenate the left and right gray graphs
if ~isempty(rightNodes)
    rightNodes(5,:) = rightNodes(5,:) + length(leftEdges);
    rightEdges = rightEdges + size(leftNodes,2);
end

nodes = [leftNodes rightNodes];
edges = [leftEdges rightEdges];

% make path local, if HOMEDIR is contained within it
% (this happens if, for instance, the vAnat directory is
% linked to from the session directory):
if strfind(leftPath,HOMEDIR);
    leftPath = leftPath(length(HOMEDIR)+2:end); % +2 b/c of filesep
    leftPath(leftPath=='\') = '/';  % \ win filesep -> / generic filesep
end

rightPath = leftPath; leftClassFile = leftPath; rightClassFile = leftPath;

% Verify save path exists - if not, make it
if ~exist(fileparts(coordsPth), 'dir'), mkdir(fileparts(coordsPth)); end

% Save to coords file
save(coordsPth,'coords','nodes','edges',...
    'allLeftNodes','allLeftEdges','allRightNodes','allRightEdges',...
    'leftPath','rightPath','keepLeft','keepRight');

if exist('leftClassFile','var'), save(coordsPth,'leftClassFile','-append'); end
if exist('rightClassFile','var'),save(coordsPth,'rightClassFile','-append');end

fprintf('Saved gray coords info to %s. \n',coordsPth);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function spm_srender(job)
% A function for rendering surfaces
% FORMAT spm_srender(job)
% job - a job structure (see tbx_cfg_render.m)
%__________________________________________________________________________
% Copyright (C) 2008-2018 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_srender.m 7381 2018-07-25 10:27:54Z guillaume $

fg  = spm_figure('GetWin','Graphics');
ren = get(fg,'Renderer');
clf(fg);
set(fg,'Renderer','OpenGL');
ax = axes('Parent',fg,'DeleteFcn',['rotate3d off; set(gcf,''Renderer'',''' ren ''');']);

try
    set(fg,'CurrentAxes',ax);
    cameratoolbar(fg);
    drawnow;
end

for i=1:numel(job.Object)
     obj = job.Object(i);
     for j=1:numel(obj.SurfaceFile)
         FVo = struct(gifti(obj.SurfaceFile{j}));
         FV  = struct('faces',FVo.faces,'vertices',FVo.vertices);
         p   = patch(FV, 'Parent',ax,...
             'FaceColor', [obj.Color.Red,obj.Color.Green, obj.Color.Blue],...
             'FaceVertexCData',  [],...
             'EdgeColor',        'none',...
             'FaceLighting',     'gouraud',...
             'EdgeLighting',     'gouraud',...
             'SpecularStrength', obj.SpecularStrength,...
             'AmbientStrength',  obj.AmbientStrength,...
             'DiffuseStrength',  obj.DiffuseStrength,...
             'SpecularExponent', obj.SpecularExponent,...
             'FaceAlpha',        obj.FaceAlpha);
    end
end
for i=1:numel(job.Light)
    obj = job.Light(i);
    l   = light('Parent',ax,...
        'Position',obj.Position,...
        'Color',[obj.Color.Red,obj.Color.Green, obj.Color.Blue]);
end

%set(0,'CurrentFigure',fg);
set(fg,'CurrentAxes',ax);
axis image equal off;
drawnow;

author: spm


text:
%ISOBJIM test if the dataset contains objects that are images
%
%  N = ISOBJIM(A)
%      ISOBJIM(A)
%
% INPUT
%  A  input dataset
%
% OUTPUT
%  N  logical value
%
% DESCRIPTION
% True if the dataset contains objects that are images. If no output is required,
% false outputs are turned into errors. This may be used for assertion.
%
% In case the objects in A should be interpreted as 1-dimensional object
% images, use  A = SETFEATSIZE(A,[1 GETFEATSIZE(A)]) for conversion.
%
% SEE ALSO (PRTools Guide)
% ISDATASET, ISDATAIM, SETFEATSIZE, GETFEATSIZE

% $Id: isobjim.m,v 1.5 2009/01/31 15:43:10 duin Exp $

function n = isobjim(a)

	n = (isa(a,'prdataset') & length(a.featsize) > 1) | isdatafile(a);

	% generate error if input is not a dataset with image data with
   % pixels being features AND no output is requested (assertion)

	if nargout == 0 & n == 0
		error([newline '---- Dataset with object images expected -----'])
	end

return

author: marianux

repo: ecg-kit

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


text:
function output = callsdpa(interfacedata)

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
x0      = interfacedata.x0;
ub      = interfacedata.ub;
lb      = interfacedata.lb;

% Bounded variables converted to constraints
if ~isempty(ub)
    [F_struc,K] = addStructureBounds(F_struc,K,ub,lb);
end

% Convert from internal (sedumi) format
[mDIM,nBLOCK,bLOCKsTRUCT,c,F] = sedumi2sdpa(F_struc,c,K);

if options.verbose==0
    options.sdpa.print = 'no';
else
    options.sdpa.print = 'display';
end

if options.savedebug
    ops = options.sdpa;
    save sdpadebug mDIM nBLOCK bLOCKsTRUCT c F ops
end

if options.showprogress;showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);end

solvertime = tic;
[objVal,x,X,Y,INFO]=sdpam(mDIM,nBLOCK,bLOCKsTRUCT,c,F,[],[],[],options.sdpa);
solvertime = toc(solvertime);

% Create variables in YALMIP internal format
Primal = x;

Dual = [];
for i = 1:length(Y)
    Dual = [Dual;Y{i}(:)];
end

Slack = [];
if options.saveduals
    for i = 1:length(X)
        Slack = [Slack;X{i}(:)];
    end
end

switch (INFO.phasevalue)
    case 'pdOPT'
        problem = 0;
    case {'noINFO','pFEAS','dFEAS'}
        problem = 3;
    case {'pdFEAS'}
        problem = 4;
    case 'pFEAS_dINF'
        problem = 2;
    case 'pINF_dFEAS'
        problem = 1;
    case 'pUNBD'
        problem = 2;
    case 'dUNBD'
        problem = 1;
    case 'pdINF'
        problem = 12;
    otherwise
        problem = -1;
end
infostr = yalmiperror(problem,interfacedata.solver.tag);

if options.savesolveroutput
    solveroutput.objVal = objVal;
    solveroutput.x = x;
    solveroutput.X = X;
    solveroutput.Y = Y;
    solveroutput.INFO = INFO;
else
    solveroutput = [];
end

if options.savesolverinput
    solverinput.mDIM = mDIM;
    solverinput.nBLOCK=nBLOCK;
    solverinput.bLOCKsTRUCT=bLOCKsTRUCT;
    solverinput.c=c;
    solverinput.F=F;
else
    solverinput = [];
end

% Standard interface
output = createOutputStructure(Primal,Dual,[],problem,infostr,solverinput,solveroutput,solvertime);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function edgenb=edgeneighbors(t,opt)
%
% edgenb=edgeneighbors(t,opt)
%
% to find neighboring triangular elements in a triangule surface
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%     t: a triangular surface element list, 3 columns of integers
%     opt: if opt='general', return the edge neighbors for a general
%          triangular surface: each edge can be shared by more than 2
%          triangles; if ignored, we assume all triangles are shared by no
%          more than 2 triangles.
%
% output:
%     edgenb: if opt is not supplied, edgenb is a size(t,1) by 3 array with
%     each element being the triangle ID of the edge neighbor of that
%     triangle. For each row, the order of the neighbors is listed as those
%     sharing edges [1 2], [2 3] and [3 1] between the triangle nodes.
%
%     when opt='general', edgenb is a cell array with a length of size(t).
%     each member of the cell array is a list of edge neighbors (the order
%     is not defined).
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

edges=[t(:,[1,2]);
       t(:,[2,3]);
       t(:,[3 1])];
edges=sort(edges,2);
[foo,ix,jx]=unique(edges,'rows');

if(nargin==2)
  if(strcmp(opt,'general'))
        ne=size(t,1);
        edgenb=cell(ne,1);
        for i=1:ne
            % this is very slow, need to be optimized
            nb=unique(mod([find(jx==jx(i) | jx==jx(i+ne) | jx==jx(i+2*ne))]',ne),'first');
            nb(nb==0)=ne;
            edgenb{i}=nb(nb~=i);
        end
        return;
  else
        error(['supplied option "' opt '" is not supported.']);
  end
end

if(isoctavemesh)
        u=unique(jx);
        qx=u(hist(jx,u)==2);
else
        vec=histc(jx,1:max(jx));
        qx=find(vec==2);
end

nn=max(t(:));
ne=size(t,1);
edgenb=zeros(size(t));

% now I need to find all repeatitive elements
% that share a face, to do this, unique('first')
% will give me the 1st element, and 'last' will
% give me the second. There will be no more than 2

% doing this is 60 times faster than doing find(jx==qx(i))
% inside a loop

[ujx,ii]=unique(jx,'first');
[ujx,ii2]=unique(jx,'last');

% iddup is the list of all pairs that share a common face

iddup=[ii(qx) ii2(qx)];
faceid=ceil(iddup/ne);
eid=mod(iddup,ne);
eid(eid==0)=ne;

% now rearrange this list into an element format

for i=1:length(qx)
	edgenb(eid(i,1),faceid(i,1))=eid(i,2);
	edgenb(eid(i,2),faceid(i,2))=eid(i,1);
end

% edgenb may contain 0s, that just means the corresponding
% face is a boundary face and has no neighbor.

% if the second option is 'surface', I am going to find
% and return surface patches only

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [V,E] = eig(T,varargin)
% compute the eigenvalues and eigenvectors of a tensor
%
% Syntax
%   E = eig(T)
%   [V,E] = eig(T)
%
% Input
%  T - list of M rank 2 @tensor
%
% Output
%  E - 3xM list of eigen values
%  V - 3xM list eigen @vector3d
%

switch T.rank

  case 1
  case 2
    if all(T.isSymmetric(:))
      [V,E] = eig3(T.M,varargin{:});
    else
      E = zeros(3,length(T));
      V = zeros(3,3,length(T));
      for i = 1:length(T)
        [V(:,:,i),E(:,i)] = eig(T.M(:,:,i),'vector');
      end
    end
  case 3
  case 4
    E = zeros(6,length(T));
    V = zeros(6,6,length(T));
    M = tensor42(T.M,T.doubleConvention);
    for i = 1:length(T)
      [V(:,:,i),E(:,i)] = eig(M(:,:,i),'vector');
    end
end

if nargout <= 1
  V = E;
elseif isa(V,'double')
  V = reshape(vector3d(V(1,:),V(2,:),V(3,:),'antipodal'),3,[]);
end

author: mtex-toolbox

repo: mtex

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


text:
function A = add_hhmm_end_state(transprob, termprob)
% ADD_HMM_END_STATE Combine trans and term probs into transmat for automaton with an end state
% function A = add_hhmm_end_state(transprob, termprob)
%
% A(i,k,j) = Pr( i->j | Qps=k), where i in 1:Q, j in 1:(Q+1), and Q+1 is the end state
% This implements the equation in sec 4.6 of my tech report, where
% transprob(i,k,j) = \tilde{A}_k(i,j), termprob(k,j) = \tau_k(j)
%
% For the top level, the k index is missing.

Q = size(transprob,1);
toplevel = (ndims(transprob)==2);
if toplevel
  Qk = 1;
  transprob = reshape(transprob, [Q 1 Q]);
  termprob = reshape(termprob, [1 Q]);
else
  Qk = size(transprob, 2);
end

A = zeros(Q, Qk, Q+1);
A(:,:,Q+1) = termprob';

for k=1:Qk
  for i=1:Q
    for j=1:Q
      A(i,k,j) = transprob(i,k,j) * (1-termprob(k,i));
    end
  end
end

if toplevel
  A = squeeze(A);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function imgOut = ExponentialTMO(img, exp_k)
%
%       imgOut = ExponentialTMO(img, exp_k)
%
%
%       Input:
%           -img: input HDR image
%           -exp_k: appearance value [1, +inf)
%
%       Output
%           -imgOut: tone mapped image
%
%     Copyright (C) 2010-15 Francesco Banterle
%
%     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/>.
%

%is it a three color channels image?
check13Color(img);

checkNegative(img);

if(~exist('exp_k', 'var'))
    exp_k = 1;
end

if(exp_k <= 0)
    exp_k = 1;
end

%Luminance channel
L = lum(img);

Lwa = logMean(L); %geometric mean

%dynamic range reduction
Ld = 1 - exp(-exp_k * ( L / Lwa));

%change luminance in img
imgOut = ChangeLuminance(img, L, Ld);

end

author: banterle


text:
function polynomial_dif_test ( )

%*****************************************************************************80
%
%% POLYNOMIAL_DIF_TEST tests POLYNOMIAL_DIF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 November 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYNOMIAL_DIF_TEST\n' );
  fprintf ( 1, '  POLYNOMIAL_DIF computes derivatives of a polynomial.\n' );

  m = 2;

  o1 = 4;
  c1 = [ 2.0, 3.0, 4.0, 5.0 ];
  e1 = [ 1, 10, 12, 32 ];
  title1 = '  P(X) =';
  fprintf ( 1, '\n' );
  polynomial_print ( m, o1, c1, e1, title1 );

  dif = [ 2, 1 ];

  [ o, c, e ] = polynomial_dif ( m, o1, c1, e1, dif );
  title = '  d3 P(X) dx1 dx1 dx2 =';
  fprintf ( 1, '\n' );
  polynomial_print ( m, o, c, e, title );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%kkedit ' '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kedit.pane file
%
% Parameters:
% InputFile: i 'Input File', optional: 'Input ASCII file'
% OutputFile: o 'Output File', optional: 'Resulting output ASCII File'
%
% Example: o = kkedit(i, {'i','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kedit - ASCII file editor
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc.  All rights reserved.
%

function varargout = kkedit(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkedit(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output'};
maxval={1,1};
minval={1,1};
istoggle=[1,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kedit"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function varargout = quantumstates(varargin)
%QUANTUMSTATES    Compute and plot Schroedinger eigenstates.
%   This program computes and plots eigenvalues lambda and eigenfunctions u
%   for the equation Lu = lambda*u, where L is the Schroedinger operator
%   defined by Lu(x) = -h^2*u"(x) + V(x)*u(x).  Here h is a small parameter
%   and the potential function V is given as a Chebfun. The domain of the
%   problem is the domain of V, with boundary conditions u=0 at both ends.
%
%   Inputs:
%
%       QUANTUMSTATES(V) plots 10 eigenstates for h=0.1
%       QUANTUMSTATES(V, n) plots n eigenstates for h=0.1
%       QUANTUMSTATES(V, h), h noninteger, plots 10 eigenstates for given h
%       QUANTUMSTATES(V, n, h) plots n eigenstates for given h
%       QUANTUMSTATES(..., 'noplot') produces no plot
%
%   Outputs:
%
%       D = QUANTUMSTATES(...) returns a vector D of eigenvalues
%       [U, D] = QUANTUMSTATES(...) returns a quasimatrix U of eigenfunctions
%       and a diagonal matrix of eigenvalues
%
%   Examples:
%
%       x = chebfun('x', [-3, 3]);
%       V = x.^2;                 % harmonic oscillator, or
%       V = abs(x);               % absolute value, or
%       V = (x.^2-1).^4;          % double well

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

%% Parsing of inputs:
noplot = strcmpi(varargin{nargin}, 'noplot');      % check for no plot
nargin1 = nargin;                                  % no of input args
if ( noplot )
    nargin1 = nargin1 - 1;
end
V = varargin{1};                                   % potential function
n = 10;                                            % default no of states
h = 0.1;                                           % default constant
if ( nargin1 == 3 )
    n = varargin{2};
    h = varargin{3};
end
if ( nargin1 == 2 )
    v2 = varargin{2};
    if ( v2 == round(v2) )
        n = v2;
    else
        h = v2;
    end
end

%% Eigenvalue computation:
[xmin, xmax] = domain(V);                          % domain of problem
% Create a CHEBOP with Dirichlet BCs
L = chebop([xmin, xmax]);
L.lbc = 0; L.rbc = 0;
L.op = @(x,u) -h^2*diff(u,2) + V.*u;               % Schroedinger operator
[U, D] = eigs(L, n, 'sr');                         % compute evals/efuns
d = diag(D);                                       % vector of evals
[d, ii] = sort(d);                                 % sort them
U = U(:,ii);

%% Outputs:
if ( nargout == 2 )
    varargout = {U, diag(d)};
else
    varargout = {d};
end

if ( noplot )
    % If we're not plotting, then we return.
    return
end

%% Plot:

holdState = ishold;

% Plot the potential function:
LW = 'linewidth'; lw = 1;
plot(V, 'k', LW, lw, 'jumpline', '-k'), hold on
s = sprintf('h = %4g      %d eigenstates', h, n);
title(s, 'fontsize', 12)

% Vertical limits:
ymax = max(d);
ymin = min(V);
ydiff = ymax - ymin;
ymax = ymax + .2*ydiff;
ymin = ymin - 0*ydiff;

% V values at endpoints:
Vxmin = feval(V, xmin);
Vxmax = feval(V, xmax);
if ( ymax > Vxmin )                              % The potential
    plot(xmin*[1 1], [ymax, Vxmin], 'k', LW, lw)  %   V(x) effectively
end                                              %   goes to infinity
if ( ymax > Vxmax )                              %   at the endpoints,
    plot(xmax*[1 1], [ymax Vxmax], 'k', LW, lw)   %   so we make the
end                                              %   plot show this.
dx = .05*(xmax - xmin);
dy = .25*ydiff/max(5, n);

% Plot the eigenfunction, lifted by the corresponding eigenvalue:
W = dy*U;
W = num2cell(W);
for j = 1:n
    umm = minandmax(W{j});
    if ( umm(2) < -umm(1) )
        W{j} = -W{j};
    end
    W{j} = W{j} + d(j);
end
W = horzcat(W{:});
plot(W, LW, lw)

% Plot V(x) again (so that black ends up on top):
plot(V, 'k', LW, lw, 'jumpline', '-k')
if ( ymax > Vxmin )
    plot(xmin*[1, 1], [ymax, Vxmin], 'k', LW, lw)
end
if ( ymax > Vxmax )
    plot(xmax*[1, 1], [ymax Vxmax], 'k', LW, lw)
end

% Set axis:
axis([xmin - dx, xmax + dx, ymin - dy, ymax]), drawnow

if ( ~holdState )
    % Stop holding:
    hold off
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% explore balloon model stability

clear all
close all
clc

xg1 = [-5:0.2:5];
xg2 = xg1;
n = length(xg1);
lm = zeros(n,n);
kas = [-4:0.2:4];
kaf = [-3:0.2:3];

try

    load demo_stability_hrf.mat
    % get max eigenvalue
    P = [0;0;0;0;0;0];
    for i=1:n
        for j=1:n
            J = VBA_numericDiff(@f_HRF,1,[xg1(i);xg2(j);0;0],P,0,[]);
            lm(i,j) = max(real(eig(J)));
        end
    end

    % get separatrix
    stab = zeros(n,n);
    stab(lm>=0) = -1;
    L = del2(stab);
    sep = zeros(n,n);
    sep(L<0)=1;
    p = length(find(L==1));
    y = [];
    for i=1:n
        for j=1:n
            if sep(i,j)==1
                y = [y,[xg1(i);xg2(j)]];
            end
        end
    end
    g_fname = @g_exp;
    options.inG.x = y(2,:)';
    options.priors.muPhi = [1;1;0];
    options.priors.SigmaPhi = eye(3);
    options.priors.SigmaPhi(1,1) = 0;
    options.priors.SigmaPhi(3,3) = 0;
    options.verbose = 0;
    options.DisplayWin = 0;
    dim.n_theta         = 0;
    dim.n_phi           = 3;
    dim.n               = 0;
    [posterior,out] = VBA_NLStateSpaceModel(y(1,:)',[],[],g_fname,dim,options);

catch

    % loop over different kas/kaf
    phi = zeros(length(kaf),length(kas));
    for k=1:length(kas)
        for l=1:length(kaf)
            k,l
            % get max eigenvalue
            P = [0;0;kaf(l);kas(k);0;0];
            for i=1:n
                for j=1:n
                    J = VBA_numericDiff(@f_HRF,1,[xg1(i);xg2(j);0;0],P,0,[]);
                    lm(i,j) = max(real(eig(J)));
                end
            end

            % get separatrix
            stab = zeros(n,n);
            stab(lm>=0) = -1;
            L = del2(stab);
            sep = zeros(n,n);
            sep(L<0)=1;
            p = length(find(L==1));
            y = [];
            for i=1:n
                for j=1:n
                    if sep(i,j)==1
                        y = [y,[xg1(i);xg2(j)]];
                    end
                end
            end
            g_fname = @g_exp;
            options.inG.x = y(2,:)';
            options.priors.muPhi = [1;1;0];
            options.priors.SigmaPhi = eye(3);
            options.priors.SigmaPhi(1,1) = 0;
            options.priors.SigmaPhi(3,3) = 0;
            options.verbose = 0;
            options.DisplayWin = 0;
            dim.n_theta         = 0;
            dim.n_phi           = 3;
            dim.n               = 0;
            [posterior,out] = VBA_NLStateSpaceModel(y(1,:)',[],[],g_fname,dim,options);
            phi(l,k) = posterior.muPhi(2);

        end
    end

end

% store for display purposes
x = options.inG.x;
gx = out.suffStat.gx;
PHI = posterior.muPhi;

gkaf = kron(ones(1,length(kas)),kaf)';
gkaf = gkaf - mean(gkaf);
gkas = kron(kas,ones(1,length(kaf)))';
gkas = gkas - mean(gkas);
options = [];
g_fname = @g_exp2d;
options.inG.gkaf = gkaf;
options.inG.gkas = gkas;
options.priors.muPhi = ones(1,1);
options.priors.SigmaPhi = 1e4*eye(1);
options.verbose = 1;
options.DisplayWin = 1;
dim.n_theta         = 0;
dim.n_phi           = 1;
dim.n               = 0;
[posterior,out] = VBA_NLStateSpaceModel(...
    phi(:),[],[],g_fname,dim,options);
I = reshape(out.suffStat.gx,length(kaf),length(kas));

% display results

hf = figure('color',[1 1 1],'position',[649,40,767,1105],'menubar','none');
ha = subplot(3,2,1,'parent',hf);
hi = imagesc(flipud(lm),'parent',ha);
xlabel(ha,'x2: blood inflow')
ylabel(ha,'x1: vasodilatory signal')
title(ha,'max eigenvalue')
colorbar('peer',ha);

ha = subplot(3,2,2,'parent',hf);
hi = plot(y(2,:),y(1,:),'.','parent',ha);
xlabel(ha,'x2: blood inflow')
ylabel(ha,'x1: vasodilatory signal')
set(ha,...
    'xlim',[min(xg1),max(xg1)],...
    'ylim',[min(xg1),max(xg1)],...
    'nextplot','add')
axis(ha,'square')
grid(ha,'on')
title(ha,...
    ['separatrix: kas=',...
    num2str(P(4)),...
    ' ; kaf=',...
    num2str(P(3))])
set(ha,'nextplot','add')
[xs,is] = sort(x);
ys = gx(is);
plot(ha,xs,ys,'k','parent',ha)
plot(0,0,'r+','parent',ha);
str = ['x(1) = ',...
    'theta',...
    '*exp(',...
    'x(2)) : theta = ',...
    num2str(PHI(2))];
legend(ha,{'numerical',str,'equilibrium'})

ha = subplot(3,2,3,'parent',hf);
hi = imagesc(phi,'parent',ha);
colorbar('peer',ha);
xlabel(ha,'kas')
ylabel(ha,'kaf')
title('separatrix parameter (theta)')

ha = subplot(3,2,4,'parent',hf);
hi = plot(kas,phi','parent',ha);
xlabel(ha,'kas')
ylabel(ha,'separatrix parameter (theta)')
grid(ha,'on');
axis(ha,'tight')
for i=1:length(kaf)
    leg{i} = ['kaf=',num2str(kaf(i))];
end
legend(ha,leg)
title(ha,'separatrix parameter (theta)')

ha = subplot(3,2,5,'parent',hf);
hi = imagesc(I,'parent',ha);
colorbar('peer',ha);
xlabel(ha,'kas')
ylabel(ha,'kaf')
title('fitted separatrix parameter (theta)')

ha = subplot(3,2,6,'parent',hf);
hi = plot(phi(:),I(:),'k.','parent',ha);
mi = min([phi(:);I(:)]);
ma = max([phi(:);I(:)]);
set(ha,'nextplot','add')
plot(ha,[mi ma],[mi ma],'r')
xlabel(ha,'sep param')
ylabel(ha,'fitted sep param (theta)')
title('sep param model fit')

str = ['theta = ',...
    num2str(posterior.muPhi(1)),...
    '*exp(kaf/2 - |kas-kaf/2|)'];
legend(ha,{'numerical',str})
grid(ha,'on');
set(ha,'xlim',[mi ma],'ylim',[mi ma])

author: MBB-team

repo: VBA-toolbox


text:
function Acc = bssfo_off2on(CNT,general_initparam,bssfo_param)
opt = opt_cellToStruct(general_initparam);
bssfo_opt = opt_cellToStruct(bssfo_param);

% Pre-processing
%Train
CNT_off = prep_selectChannels(CNT{1}, {'Index', opt.channel_index});
CNT_off =prep_filter(CNT_off , {'frequency', bssfo_opt.init_band});
SMT_off = prep_segmentation(CNT_off, {'interval', opt.time_interval});

%Test
CNT_on= prep_selectChannels(CNT{2}, {'Index', opt.channel_index});
CNT_on =prep_filter(CNT_on , {'frequency', bssfo_opt.init_band});
SMT_on = prep_segmentation(CNT_on, {'interval', opt.time_interval});
%% BSSFO - optimal band and weight
[FilterBand]=func_bssfo(SMT_off, {'classes', {'right', 'left'}; ...
    'frequency', {bssfo_opt.mu_band,bssfo_opt.beta_band}; 'std', {5, 25}; 'numBands', bssfo_opt.numBands; ...
    'numCSPPatterns', opt.CSPFilter; 'numIteration', bssfo_opt.numIteration});

%% BSSFO - CSP filter from optimal band
for iii=1:bssfo_opt.numIteration
    CNTtr=prep_filter(CNT{1}, {'frequency', FilterBand.sample(:,iii)'});
    CNTtr= prep_selectChannels(CNTtr, {'Index', opt.channel_index});
    SMT=prep_segmentation(CNTtr, {'interval', opt.time_interval});
    [SMT, CSP_W, CSP_D]=func_csp(SMT,{'nPatterns', opt.CSPFilter});
    FT=func_featureExtraction(SMT, {'feature','logvar'});
    [CF_PARAM]=func_train(FT,{'classifier','LDA'});
    clear CNTclass CNTtr SMT FT CSP_D

    CNTte=prep_filter(CNT{2}, {'frequency', FilterBand.sample(:,iii)'});
    CNTte= prep_selectChannels(CNTte, {'Index', opt.channel_index});
    SMTte=prep_segmentation(CNTte, {'interval', opt.time_interval});
    SMTfb=func_projection(SMTte, CSP_W);
    FTfb=func_featureExtraction(SMTfb, {'feature','logvar'});
    [cf_out]=func_predict(FTfb, CF_PARAM);
    ensemble(:,iii)=cf_out.*FilterBand.weight(iii); %
    clear CNTte SMTte FTfb CSP_W
end
sum_ensemble = sum(ensemble,2); %
%%
[loss out]=eval_calLoss(CNT{2}.y_dec, sum_ensemble); %
Acc=1-loss;
end

author: PatternRecognition

repo: OpenBMI


text:
%| Gtomo3_3s_test.m
%| Test the Gtomo3 object for 3d (SPECT) with many threads vs 1 thread
%|
%| 2017-07-11 updated to test new capability of multiple objects

% f3d_mex('chat', int32(2)) % for debugging

if 0 || ~isvar('x'), printm 'x'
	ig = image_geom('nx', 16, 'ny', 16, 'nz', 10, 'dx', 4, 'dz', 4);
	ig.mask = ig.circ > 0;
	x = single(convn(double(ig.mask), ones(3,3,1)/3^2, 'same') >= 1);
	x(end/4, end/2, 3) = 2;
	im plc 2 3
	im(1, ig.mask)
	im(4, x)
end

if 0 || ~isvar('A1'), printm 'A1'
%	f.chat = 998; % debug
	f.chat = 0;
	f.option = {};

	f3d_mex('free:all') % for testing, want only these two in play

	if 1
		if 1
			f.fwhm_collimator = 1;
			f.fwhm_iso = 2; % depth-dependent gaussian blur
			f.psfs = '-';
			f.blur = sprintf(',gauss,%g,%g', ...
				f.fwhm_collimator, f.fwhm_iso);
		elseif 0
			f.psfs = '-';
			f.blur = ',none';
		else % stress fftw
			f.psfs = '/tmp/t,psfs.fld';
			psfs = make_3s_psfs(ny, 1, 1.2*nx, 0, 2/nx);
			f.blur = ',fft'; % fails due to fftw issues?
			fld_write(f.psfs, psfs)
		end
%		mask = []; % for 3s object!
		f.na = 6;
		f.mumap = '-';
		f.sfilter = 1;
		dx = 4;
		f.sys_type = '3s@%g,%g,%g,360,0,%d%s@%s@%s@-%d,%d,%d';
		f.sys_type = sprintf(f.sys_type, ig.dx, ig.dx, ig.dz, ...
			f.sfilter, f.blur, f.mumap, f.psfs, ig.nx, ig.nz, f.na);

%		3s@[-|sx,sy,sz,orbit,orbit_start,spline_filter[,blur_method]]
%			@mumap.file@filter.file@-nu,nv,nview
	end

	f.option = {f.option{:}, 'chat', f.chat, 'checkmask', 0};

	A1 = Gtomo3(f.sys_type, ig.mask, ig.nx, ig.ny, ig.nz, ...
		f.option{:}, 'nthread', 1);

	tester_tomo2(A1, ig.mask) % put it through paces
end

if 1 || ~isvar('y1'), printm 'y1'
	cpu etic
	y1 = A1 * x;
	t1 = cpu('etoc', sprintf('%d threads', A1.arg.nthread));
	im(2, y1)
end

if 0 || ~isvar('A2'), printm 'A2'
	A2 = Gtomo3(f.sys_type, ig.mask, ig.nx, ig.ny, ig.nz, ...
		f.option{:}, ...
			'nthread', 2+2*jf('ncore'), ...
			'nthread_max', 100);

	tester_tomo2(A2, ig.mask, 'A2', A1) % put it through paces / compare
end

if 0 || ~isvar('y2'), printm 'y2'
	y2 = A2 * x;
	t2 = cpu('etoc', sprintf('%d threads', A2.arg.nthread));

	im(5, y2)
	im(6, y2-y1)

        jf_equal(y1, y2)
end

f3d_mex('show:all')
if 1 % clean up after testing
	clear A1 A2
	f3d_mex('free:all')
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function volumeStruct = readLeksellVolumeFile(filename)
%"readLeksellVolumeFile"
%   Reads a Leksell volume file using decodeLeksellData and places the
%   fields into a datastructure with meaningful names.  This file contains
%   all contouring information for CERR.
%
%JRA 6/13/05
%
%LM: KRK, 8/11/07,  Added scaling on the xy contour positions to change mm
%                   to cm.
%
%Usage:
%   function volumeStruct = readLeksellVolumeFile(filename)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

fid = fopen(filename, 'r', 'b');

data = decodeLeksellData(fid);

fclose(fid);

%First iterate over structures included in volume file.
for i = 1:length(data) - 1

    rawStruct = data{i};

    structureS(i).structName = rawStruct{2};
    structureS(i).registeredTo = rawStruct{1};

    allContours = rawStruct{4};

    contourS = [];

    %Iterate over contours for this structure
    for j=1:length(allContours)-1

        oneContour = allContours{j};

        contourS(j).sliceNum = oneContour{1};
        contourS(j).nVertices = oneContour{2};

        cData     = oneContour{3};
        xyMat = [];
        for k = 1:length(cData)
            xyMat = [xyMat;cData{k}{1}'];
        end

        contourS(j).contour = xyMat;
        contourS(j).rowExtent = [oneContour{6} oneContour{7}];
        contourS(j).colExtent = [oneContour{8} oneContour{9}];
        contourS(j).slcExtent = [oneContour{10} oneContour{11}];

    end

    structureS(i).contour = contourS;

end

volumeStruct = structureS;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [acc, cm, classcount, passthresh] = ...
    msAnalyzeResult(imsegs, labels, pg, DO_SP, ignore, thresh, niter)

nclasses = size(pg{1}, 2);
cm = zeros(nclasses, nclasses);
total = 0;
acc = 0;

passthresh = 0;
pttotal = 0;

if ~exist('thresh', 'var') || isempty(thresh)
    thresh = 0;
end

for f = 1:numel(pg)

    if exist('niter', 'var') && ~isempty(niter)
        pg{f} = mean(pg{f}(:, :, 1:niter), 3);
    end

    lab = labels{f}(:);

    if ~isempty(imsegs)
        npixels = imsegs(f).npixels(:) / sum(imsegs(f).npixels);
    end

    pg{f}(:, ignore) = 0;

    ind = false(size(lab));
    for k = 1:numel(ignore)
        ind(lab==ignore(k)) = true;
    end
    lab(ind) = 0;
    [maxval, maxlab] = max(pg{f}, [], 2);

    if DO_SP

        pttotal = pttotal + sum(npixels.*(lab>0));
        lab(maxval<thresh) = 0; %ignore labels assigned below threshold confidence
        passthresh = passthresh + sum(npixels.*(lab>0));

        acc = acc + sum((lab==maxlab).*npixels);
        total = total + sum(npixels((lab>0)));

        for s = 1:numel(maxlab)
            if lab(s)~=0
                cm(lab(s), maxlab(s)) = cm(lab(s), maxlab(s)) + npixels(s);
            end
        end
    else

        if ~isempty(imsegs)
            maxlab = maxlab(imsegs(f).segimage);
            maxval = maxval(imsegs(f).segimage);
        end

        pttotal = pttotal + sum((lab>0)/numel(lab));
        lab(maxval<thresh) = 0; %ignore labels assigned below threshold confidence
        passthresh = passthresh + sum((lab>0)/numel(lab));

%           maxlab2 = maxlab; maxlab2(maxval<0.5) = 0;
%         figure(2), hold off, imagesc(maxlab2, [0 21]), axis image
%         figure(3), imagesc(maxlab, [0 21]), axis image
%         figure(4), imagesc(reshape(lab, size(maxlab)), [0 21]), axis image, colormap jet
%         pause

        acc = acc + sum(lab==maxlab(:))/numel(lab);
        total = total + sum(lab(:)>0)/numel(lab);
        for y1 = 1:nclasses
            ind = lab==y1;
            if any(ind)
                for y2 = 1:nclasses
                    ind2 = ind & (maxlab(:)==y2);
                    cm(y1, y2) = cm(y1, y2) + sum(ind2)/numel(lab);
                end
            end
        end
    end
end

acc = acc / total;

classcount = sum(cm, 2)';
classcount = classcount / sum(classcount);

cm = cm ./ max(repmat(sum(cm, 2), 1, size(cm, 2)), 1E-9);

passthresh = passthresh / pttotal;

%disp(['vacc: ' num2str(vacc) '   hacc: ' num2str(hacc)])
%disp(num2str(sum(ny, 1) / sum(totalv)))

author: zouchuhang

repo: LayoutNet


text:
function Z = ECGModel(X,Phasemn)
% Z = ECGModel(X,Phaemn)
% Overview: Generate Synthetic ECG beat using parameters defined in X over phase
% defined in the Phasemn variable.
%
% Inputs:
%           X - contains the amplitude (alphai), standard deviation (bi) and phase (tetai) for
%           each Gaussian function used for generating the synthetic ECG beat.
%
%           Phasemn - array contains the phase points over which synthetic
%           ECG beat is generated.
%
% Outputs:
%           Z - synthetic ECG beat generated using Gaussian parameter estimates.
%
% Open Source ECG Toolbox, version 1.0, November 2006
% Released under the GNU General Public License
% Copyright (C) 2006  Reza Sameni
% Sharif University of Technology, Tehran, Iran -- LIS-INPG, Grenoble, France
% reza.sameni@gmail.com
% Last modified:
% on 11/30/2020 by Ismail Sadiq
%
% TODO: Return the analytic Jacobians to accelerate the convergence of the
% nonlinear least squares solver. Added on Feb. 2019
%
% 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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

L = (length(X)/3);

alphai = X(1:L);
bi = X(L+1:2*L);
tetai = X(2*L+1:3*L);

Z = zeros(size(Phasemn));
for j = 1:length(alphai),
    dtetai = rem(Phasemn - tetai(j) + pi,2*pi)-pi;
    Z = Z + alphai(j) .* exp(-dtetai .^2 ./ (2*bi(j) .^ 2));
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
%GRIDBDRE Generate grid edge boundary information.
%
%   [ BE, E, EV ] = GRIDBDRE( B, C ) Given arrays B and C with cell
%   connectivities and face boundary information generates an array
%   with edge boundary information.
%
%       Input       Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       b           (3/6,n_b)              Array with boundary information, the first row
%                                          indicates the cell number, followed by cell
%                                          face, boundary number, and optionally the
%                                          normal vectors for each boundary face
%       c           (n_pcell,n_cells)      Cell connectivity pointer array where the
%                                          entries in each column correspond to the
%                                          cell vertex numbers.
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       be          (5,n_b)                Array with boundary information, the first row
%                                          indicates the global cell number, followed by
%                                          local cell face, local edge, and global edge
%                                          and edge boundary numbers.
%       e           (n_e,n_c)              Array with global edge numbers for all cells.
%       ev          (n_e,2)                Array with indices to vertices for each edge.
%
%   See also GRIDBDR, GRIDEDGE

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function scanQuery(command,varargin)
% scanQuery
%
% Queries the scan Intensity values and displays in the CERR Status Bar
%
% written DK
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC stateS

switch upper(command)
    case 'SCANQUERYSTART'
        cP = get(gcbo, 'CurrentPoint');
        %delete([findobj('tag', 'scanQueryPoint')]);
        delete(stateS.handle.scanQueryPoint)
        stateS.handle.scanQueryPoint = [];
        %line([cP(1,1) cP(1,1)], [cP(2,2) cP(2,2)], 'tag', 'scanQueryPoint', 'userdata', gcbo, 'eraseMode', 'xor', 'parent', gcbo, 'marker', '+', 'color', [1 1 1], 'hittest', 'off');
        stateS.handle.scanQueryPoint = line([cP(1,1) cP(1,1)], [cP(2,2) cP(2,2)], 'tag', 'scanQueryPoint', 'userdata', gcbo, 'parent', gcbo, 'marker', '+', 'color', [1 1 1], 'hittest', 'off');
        return;

    case 'SCANQUERYMOTION'
        %dQP = findobj('tag', 'scanQueryPoint');
        hAxis = get(stateS.handle.scanQueryPoint, 'userdata');
        [view, coord, scanSets] = getAxisInfo(hAxis, 'view', 'coord', 'scanSets');

        if isempty(scanSets)
            CERRStatusString('Cannot query scan in this axis: no scan is being displayed.')
            return;
        end

        scanSet = scanSets(1);
        if isempty(varargin)
            cP = get(hAxis, 'CurrentPoint');
            set(stateS.handle.scanQueryPoint, 'XData', [cP(1,1) cP(1,1)]);
            set(stateS.handle.scanQueryPoint, 'YData', [cP(2,2) cP(2,2)]);
        else
            xd = get(stateS.handle.scanQueryPoint, 'XData');
            yd = get(stateS.handle.scanQueryPoint, 'YData');
            cP = [xd(:) yd(:)];
        end

        switch lower(view)
            case 'transverse'
                x = cP(1,1); y = cP(2,2); z = coord;
            case 'sagittal'
                y = cP(1,1); z = cP(2,2); x = coord;
            case 'coronal'
                x = cP(1,1); z = cP(2,2); y = coord;
            otherwise
                return;
        end

        %Get scan's transM, and convert requested point to scan coords.
        transM = getTransM('scan', scanSet, planC);
        [xD, yD, zD] = applyTransM(inv(transM), x, y, z);

        %Get the actual scan value using the converted point.
        scan = getScanAt(scanSet,xD,yD,zD,planC);
        indexS = planC{end};
        %imageType = planC{indexS.scan}(scanSet).scanInfo(1).imageType;
        if ~isempty(planC{indexS.scan}(scanSet).scanInfo(1).CTOffset)
            scan = scan - planC{indexS.scan}(scanSet).scanInfo(1).CTOffset;
        end
%         if strfind(upper(imageType), 'CT')
%             CTOffset = planC{indexS.scan}(scanSet).scanInfo(1).CTOffset;
%             scan = scan - CTOffset;
%         end

        CERRStatusString(['x = ' num2str(x) ', y = ' num2str(y) ', z = ' num2str(z) ' scan: ' num2str(scan)], 'gui');
        return;

    case 'SCANQUERYMOTIONDONE'
        hFig = gcbo;
        set(hFig, 'WindowButtonMotionFcn', '', 'WindowButtonUpFcn', '');
        return;
end


text:
classdef meanshift_G_is_Gaussian < dagnn.ElementWise
    properties
        delta=0.1
    end
    properties % (Transient)
        numInputs
        SIZE_=[]
    end

    % [TODO]: current version only supports batchSize=1; need to extend to
    % multiple input images
    methods
        function outputs = forward(obj, inputs, params)
            obj.numInputs = numel(inputs);
            % obj.SIZE_ = inputs{2};
            %outputs{1} = exp((inputs{1}-1)/(obj.delta^2));

            outputs{1} = exp( -0.5*(inputs{1}.^2)/(obj.delta^2));
        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
            derInputs = cell(1, numel(inputs));
            dzdy = derOutputs{1};
            S = inputs{1};
            %derInputs{1} = dzdy.* ((1/obj.delta^2)* exp((S-1)/(obj.delta^2)));

            derInputs{1} = dzdy.* (-1/(obj.delta^2)*inputs{1}).*exp(-0.5*(inputs{1}.^2)/(obj.delta^2));

            derParams = {} ;
        end

        function obj = meanshift_G_is_Gaussian(varargin)
            obj.load(varargin) ;
        end
    end
end

author: aimerykong


text:
function fibredir = GetFibreOrientation(modelname, fittedpars)
%
% function fibredir = GetFibreOrientation(modelname, fittedpars)
%
% Returns the fibre orientation as a unit vector from fitted parameters.
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

theta = fittedpars(:, GetParameterIndex(modelname, 'theta'));
phi = fittedpars(:, GetParameterIndex(modelname, 'phi'));

fibredir = [cos(phi).*sin(theta) sin(phi).*sin(theta) cos(theta)]';

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
%**************************************************************
%* mex interface to Andy Liaw et al.'s C code (used in R package randomForest)
%* Added by Abhishek Jaiantilal ( abhishek.jaiantilal@colorado.edu )
%* License: GPLv2
%* Version: 0.02
%
% Calls Regression Random Forest
% A wrapper matlab file that calls the mex file
% This does prediction given the data and the model file
%**************************************************************

function Y_hat = regRF_predict(X,model)
    %function Y_hat = regRF_predict(X,model)
    %requires 2 arguments
    %X: data matrix
    %model: generated via regRF_train function
	if nargin~=2
		error('need atleast 2 parameters,X matrix and model');
	end

	Y_hat = mexRF_predict(X',model.lDau,model.rDau,model.nodestatus,model.nrnodes,model.upper,model.avnode,model.mbest,model.ndtree,model.ntree);

    if ~isempty(find(model.coef)) %for bias corr
        Y_hat = model.coef(1) + model.coef(2)*Y_hat;
    end
	clear mexRF_predict

repo: sr-metric


text:
function A = repmatmatch(a,B)
%Replicate and tile an array to match the size of a given N-D array
%
% Syntax:	A = repmatmatch(a,B)
%
% Inputs:
% 	a - Input array to tile
% 	B - N-D array to match the size of
%
% Outputs:
% 	A - The replicated and tiled copy of input 'a'
%
% Example:
% 	a = [1 2 3];
% 	B = rand(2,3);
% 	A = repmatmatch(a,B);
%
% See also: repmat

% Author: Jacob Donley
% University of Wollongong
% Email: jrd089@uowmail.edu.au
% Copyright: Jacob Donley 2017
% Date: 28 August 2017
% Version: 0.1 (28 August 2017)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

aSz = size(a);
BSz = size(B);

if any(rem(BSz ./ aSz,1))
    error(['The size of each dimension of ''a'' should match ''B'' ' ...
       'or be evenly divisible by the corresponding dimension of ''B''.'])
end

A = repmat(a, BSz ./ aSz);

end

author: jtkim-kaist

repo: Speech-enhancement


text:
function [ft3] = gal2ft3(gal)
% Convert volume from US liquid gallons to cubic feet.
% Chad Greene 2012
ft3 = gal*0.13368055556;

author: Sable

repo: mcbench-benchmarks


text:
function class_info=gen_class_info_ade()

class_info_ADE=load('class_info_ADE.mat');
class_info_ADE=class_info_ADE.class_info_ADE;

class_names=class_info_ADE.Name;
assert(size(class_names, 2)==1);
class_names=cat(1, {'void'}, class_names);

class_label_values=uint8([0 1:150]);

class_info=[];

class_info.class_names = class_names;

class_info.class_label_values=class_label_values;
class_info.background_label_value=uint8(1);
class_info.void_label_values=uint8(0);

% addpath ../libs/VOCdevkit_2012/VOCcode
class_info.mask_cmap = VOClabelcolormap(256);

class_info=process_class_info(class_info);

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
%% How do I know whether AR1 noise is an appropriate null hypothesis to test against?
% It is usually an appropriate null hypothesis if the theoretical AR1
% spectrum is  degreesa good model degrees for the power decay in the observed spectrum.
% I recommend to simply visually compare the two power spectra:

X=rednoise(200,.8);
[P,freq]=pburg(zscore(X),7,[],1);
aa=ar1(X);
Ptheoretical=(1-aa.^2)./(abs(1-aa.*exp(-2*pi*i*freq))).^2;
semilogy(freq,P/sum(P),freq,Ptheoretical/sum(Ptheoretical),'r');
legend('observed',sprintf('Theoretical AR1=%.2f',aa),'location','best')

author: grinsted

repo: wavelet-coherence


text:
function net = cnn_cifar_init(varargin)
opts.networkType = 'simplenn' ;
opts = vl_argparse(opts, varargin) ;

lr = [.1 2] ;

% Define network CIFAR10-quick
net.layers = {} ;

% Block 1
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{0.01*randn(5,5,3,32, 'single'), zeros(1, 32, 'single')}}, ...
                           'learningRate', lr, ...
                           'stride', 1, ...
                           'pad', 2) ;
net.layers{end+1} = struct('type', 'pool', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', [0 1 0 1]) ;
net.layers{end+1} = struct('type', 'relu') ;

% Block 2
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{0.05*randn(5,5,32,32, 'single'), zeros(1,32,'single')}}, ...
                           'learningRate', lr, ...
                           'stride', 1, ...
                           'pad', 2) ;
net.layers{end+1} = struct('type', 'relu') ;
net.layers{end+1} = struct('type', 'pool', ...
                           'method', 'avg', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', [0 1 0 1]) ; % Emulate caffe

% Block 3
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{0.05*randn(5,5,32,64, 'single'), zeros(1,64,'single')}}, ...
                           'learningRate', lr, ...
                           'stride', 1, ...
                           'pad', 2) ;
net.layers{end+1} = struct('type', 'relu') ;
net.layers{end+1} = struct('type', 'pool', ...
                           'method', 'avg', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', [0 1 0 1]) ; % Emulate caffe

% Block 4
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{0.05*randn(4,4,64,64, 'single'), zeros(1,64,'single')}}, ...
                           'learningRate', lr, ...
                           'stride', 1, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'relu') ;

% Block 5
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{0.05*randn(1,1,64,10, 'single'), zeros(1,10,'single')}}, ...
                           'learningRate', .1*lr, ...
                           'stride', 1, ...
                           'pad', 0) ;

% Loss layer
net.layers{end+1} = struct('type', 'softmaxloss') ;

% Meta parameters
net.meta.inputSize = [32 32 3] ;
net.meta.trainOpts.learningRate = [0.05*ones(1,30) 0.005*ones(1,10) 0.0005*ones(1,5)] ;
net.meta.trainOpts.weightDecay = 0.0001 ;
net.meta.trainOpts.batchSize = 100 ;
net.meta.trainOpts.numEpochs = numel(net.meta.trainOpts.learningRate) ;

% Fill in default values
net = vl_simplenn_tidy(net) ;

% Switch to DagNN if requested
switch lower(opts.networkType)
  case 'simplenn'
    % done
  case 'dagnn'
    net = dagnn.DagNN.fromSimpleNN(net, 'canonicalNames', true) ;
    net.addLayer('error', dagnn.Loss('loss', 'classerror'), ...
             {'prediction','label'}, 'error') ;
  otherwise
    assert(false) ;
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
%CHECK12 Decide between monadic and dyadic datafile operations
%
% [CHECK,A,B] = CHECK12(A,B)

author: marianux

repo: ecg-kit


text:
function determ = r8mat_determinant ( n, a )

%*****************************************************************************80
%
%% R8MAT_DETERMINANT computes the determinant of a matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the orderof the matrix A.
%
%    Input, real A(N,N), the matrix.
%
%    Output, real DETERM, the determinant.
%
  b(1:n,1:n) = a(1:n,1:n);

  [ b, pivot, info ] = r8mat_gefa ( b, n );

  determ = r8mat_gedet ( b, n, pivot );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Sfi = GetSf(angles, offsets, T2f, SfTable)
%GetSf interpolate Sf values from precomputed table

Sfi = zeros(length(angles),1);
printed = false;

% This is a little workspace trick to limit the number of
% Sf table interpolation warnings to the command window/console.
% If the counterSfMiss variable exists in the base workspace,
% its value will nbe subjected to evaluation. If not, will be assigned
% with 0 and broadcasted to the base workspace.
% This is not a go-to *.m practice, especially if done without enough
% comments.

for ii = 1:length(angles)

[xi, yi, zi] = meshgrid(offsets(ii), angles(ii), T2f);
Sfi(ii) = interp3(SfTable.offsets, SfTable.angles, SfTable.T2f, SfTable.values, xi, yi, zi);

    if (isnan(Sfi(ii)))

        if ~evalin('base','exist(''counterSfMiss'')')
          counterSfMiss = 0;
          assignin('base','counterSfMiss',counterSfMiss);
         else
         counterSfMiss  = evalin('base','counterSfMiss');

        % Print this once for all the angles. Allow 10 global prints in total.
        % Fetch the variable from workspace. Doing this here instead of L14
        % means one less condition, which matters when there are thousands.

        if ~printed && counterSfMiss < 11
         % Get value from base kspace.
          cprintf('magenta','Cannot interpolate value from current Sf table : angle: %f; offset: %f; T2f: %f\n',angles(ii), offsets(ii), T2f);
          cprintf('blue','%s','Calculating the missing Sf value...');
          if counterSfMiss==10
            cprintf('blue','%s','Remaining warnings for missing Sf value interpolations have been silenced for this processs');
          end
          counterSfMiss = counterSfMiss + 1;
          assignin('base','counterSfMiss',counterSfMiss);
          printed = true;
        end

        end

        MTpulse = GetPulse(angles(ii),offsets(ii),SfTable.PulseTrf,SfTable.PulseShape,SfTable.PulseOpt);
        Sfi(ii) = computeSf(T2f, MTpulse);
    end
end


text:
function results = PTKGetRadiusForAirways(centreline_results, lung_image, radius_approximation, reporting, figure_airways_3d)
    % PTKGetRadiusForAirways. Computes the radius for a segmented airway tree.
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    reporting.ShowProgress('Computing radius for each branch');

    [radius_results, skeleton_tree] = GetRadius(lung_image, centreline_results.IntermediateAirwaySkeleton, radius_approximation, reporting, figure_airways_3d);
    centreline_tree_model = PTKTreeModel.CreateFromSkeletonTree(skeleton_tree, lung_image);

    results = rmfield(centreline_results, 'IntermediateAirwaySkeleton');
    results.AirwayCentrelineTree = centreline_tree_model;
end

function [results, airway_skeleton] = GetRadius(lung_image, airway_skeleton, radius_approximation, reporting, figure_airways_3d)
    results = {};

    number_of_segments = airway_skeleton.CountBranches;
    segments_done = 0;

    lung_image_as_double = lung_image.BlankCopy;
    lung_image_as_double.ChangeRawImage(double(lung_image.RawImage));

    segments = airway_skeleton.GetBranchesAsListByGeneration;

    for next_segment = segments
        reporting.UpdateProgressValue(round(100*segments_done/number_of_segments));
        segments_done = segments_done + 1;

        next_result = PTKComputeRadiusForBranch(next_segment, lung_image_as_double, radius_approximation, figure_airways_3d, reporting);
        results{end+1} = next_result;

        next_segment.Radius = next_result.Radius;
        next_segment.WallThickness = next_result.WallThickness;

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [gal] = oz2gal(oz)
% Convert volume from US liquid ounces to US liquid gallons.
% Chad Greene 2012
gal = oz/128;


text:
function varargout = xyz
%XYZ   Three chebfun3 obejcts of the identity on [-1, 1, -1, 1, -1, 1].
%   CHEB.XYZ is shorthand for the expressions
%   X = CHEBFUN3(@(X,Y,Z) X),
%   Y = CHEBFUN3(@(X,Y,Z) Y), and
%   Z = CHEBFUN3(@(X,Y,Z) Z).

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

x = chebfun3(@(x,y,z) x);
y = chebfun3(@(x,y,z) y);
z = chebfun3(@(x,y,z) z);

if nargout>1
    % For the syntax [x,y,z] = cheb.xyz
    varargout{1} = x;
    varargout{2} = y;
    varargout{3} = z;
else
    % For the syntax cheb.xyz we still want to put these variables into the
    % workspace:
    assignin('base', 'x', x);
    assignin('base', 'y', y);
    assignin('base', 'z', z);
end

end

author: chebfun

repo: chebfun


text:
function [L, D, P, rho] = ldlt_skew(A)
%LDLT_SKEW  Block LDL^T factorization for a skew-symmetric matrix.
%           Given a real, skew-symmetric A,
%           [L, D, P, RHO] = LDLT_SKEW(A) computes a permutation P,
%           a unit lower triangular L, and a block diagonal D
%           with 1x1 and 2x2 diagonal blocks, such that P*A*P' = L*D*L'.
%           A partial pivoting strategy of Bunch is used.
%           RHO is the growth factor.

%           Reference:
%           J. R. Bunch, A note on the stable decomposition of skew-symmetric
%              matrices. Math. Comp., 38(158):475-479, 1982.
%           N. J. Higham, Accuracy and Stability of Numerical Algorithms,
%              Second edition, Society for Industrial and Applied Mathematics,
%              Philadelphia, PA, 2002; chap. 11.

%           This routine does not exploit skew-symmetry and is not designed to
%           be efficient.

if ~isreal(A) | ~isequal(triu(A,1)',-tril(A,-1))
    error('Must supply real, skew-symmetric matrix.')
end

n = length(A);
k = 1;
D = zeros(n);
L = eye(n);
pp = 1:n;
if nargout >= 4
   maxA = norm(A(:), inf);
   rho = maxA;
end

while k < n

      if max( abs(A(k+1:n,k)) ) == 0

         s = 1;
         % Nothing to do.

      else

         s = 2;

         if k < n-1
            [colmaxima, rowindices] = max( abs(A(k+1:n, k:k+1)) );
            [biggest, colindex] = max(colmaxima);
            row = rowindices(colindex)+k; col = colindex+k-1;

            % Permute largest element into (k+1,k) position.
            % NB: k<->col permutation must be done before k+1<->row one.
            A( [k, col], : ) = A( [col, k], : );
            A( :, [k, col] ) = A( :, [col, k] );
            A( [k+1, row], : ) = A( [row, k+1], : );
            A( :, [k+1, row] ) = A( :, [row, k+1] );
            L( [k, col], : ) = L( [col, k], : );
            L( :, [k, col] ) = L( :, [col, k] );
            L( [k+1, row], : ) = L( [row, k+1], : );
            L( :, [k+1, row] ) = L( :, [row, k+1] );
            pp( [k, col] ) = pp( [col, k] );
            pp( [k+1, row] ) = pp( [row, k+1] );
         end

         E = A(k:k+1,k:k+1);
         D(k:k+1,k:k+1) = E;
         C = A(k+2:n,k:k+1);
         temp = C/E;
         L(k+2:n,k:k+1) = temp;
         A(k+2:n,k+2:n) = A(k+2:n,k+2:n) + temp*C';  % Note the plus sign.
         % Restore skew-symmetry.
         A(k+2:n,k+2:n) = 0.5 * (A(k+2:n,k+2:n) - A(k+2:n,k+2:n)');

         if nargout >= 4, rho = max(rho, max(max(abs(A(k+2:n,k+2:n)))) ); end

     end

     k = k + s;
     if k >= n-2, D(k:n,k:n) = A(k:n,k:n); break, end;

end

if nargout >= 3, P = eye(n); P = P(pp,:); end
if nargout >= 4, rho = rho/maxA; end

author: opencobra

repo: cobratoolbox


text:
function dx = not2dx ( x1, y1, x2, y2, x3, y3, x4, y4 )

%*****************************************************************************80
%
%% NOT2DX evaluates a factor for serendipity basis functions.
%
%  Discussion:
%
%    not2(x1,y1,x2,y2,x3,y3,x4,y4) evaluates at the point (x1,y1), the basis
%    factor that is 0 at (x2,y2) and (x3,y3) and 1 at (x4,y4):
%
%          ( ( x1 - x2 ) * ( y3 - y2 ) - ( x3 - x2 ) * ( y1 - y2 ) )
%        / ( ( x4 - x2 ) * ( y3 - y2 ) - ( x3 - x2 ) * ( y4 - y2 ) )
%
%    not2dx returns the derivative of this function with respect to X1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X1, Y2, the evaluation point.
%
%    Input, real X2, Y2, X3,Y3, values that define the factor.
%
%    Output, real DX, the derivative of the basis function factor
%    with respect to X1.
%
  dx = (   1.0       * ( y3 - y2 ) +                 0.0       ) ...
     / ( ( x4 - x2 ) * ( y3 - y2 ) - ( x3 - x2 ) * ( y4 - y2 ) );

  return
end


text:
function [ is_ok, text_error, total_dur ] = checkTiming( system, varargin )
%checkTiming(sys, objects, ...)
%   Function checks whether timing of the specified evets is aligned to
%   the corresponding raster

    if (isempty(varargin))
        text_error=['empty or damaged block detected' ];
        is_ok=false;
        return;
    end
    total_dur=mr.calcDuration(varargin{:});
    is_ok=div_check(total_dur,system.blockDurationRaster);
    if (is_ok)
        text_error=[];
    else
        text_error=['total duration:' num2str(total_dur*1e6) 'us' ];
    end
    for i=1:length(varargin)
        e=varargin{i};
        if isnumeric(e) % special handling for blockDuration
            continue;
        end
        assert(isstruct(e), 'wrong format of the variable aguments, list of structures is expected');
        ok=true;
        if length(e)>1
            % for now this is only the case for arrays of extensions, but
            % we actually cannot check extensons anyway...
            continue;
        end
        if isfield(e, 'type') && (strcmp(e.type,'adc') || strcmp(e.type,'rf'))
            raster=system.rfRasterTime;
        else
            raster=system.gradRasterTime;
        end
        if isfield(e, 'delay')
            if e.delay<-eps
                ok=false;
            end
            if ~div_check(e.delay,raster)
                ok=false;
            end
        end
        if isfield(e, 'duration')
            if ~div_check(e.duration,raster)
                ok=false;
            end
        end
        if isfield(e, 'dwell') % special case ADC
            if e.dwell<system.adcRasterTime || abs(round(e.dwell/system.adcRasterTime)*system.adcRasterTime - e.dwell)>1e-10
                ok=false;
            end
        end
        if isfield(e, 'type') && strcmp(e.type,'trap')
            if ~div_check(e.riseTime, system.gradRasterTime) || ~div_check(e.flatTime, system.gradRasterTime) || ~div_check(e.fallTime, system.gradRasterTime)
                ok=false;
            end
        end
        if ~ok
            is_ok=false;
            if ~isempty(text_error)
                text_error = [text_error ' '];
            end
            text_error = [text_error '[ '];
            if isfield(e, 'type')
                text_error = [text_error 'type:' e.type ' ' ];
            end
            if isfield(e, 'delay')
                text_error = [text_error 'delay:' num2str(e.delay*1e6) 'us ' ];
            end
            if isfield(e, 'duration')
                text_error = [text_error 'duration:' num2str(e.duration*1e6) 'us ' ];
            end
            if isfield(e, 'dwell')
                text_error = [text_error 'dwell:' num2str(e.dwell*1e9) 'ns ' ];
            end
            if isfield(e, 'type') && strcmp(e.type,'trap')
                text_error = [text_error 'riseTime:' num2str(e.riseTime*1e6) 'us flatTime:' num2str(e.flatTime*1e6) 'us fallTime:' num2str(e.fallTime*1e6) 'us '];
            end
            text_error = [text_error ']'];
        end
    end
end

function out = div_check(a, b)
%  checks wheher a can be divided by b to an accuracy of 1e-9
    c = a / b;
    out = (abs( c - round(c) ) < 1e-9);
end

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
% [INPUT]
% data = A float t-by-n matrix (-Inf,Inf) representing the model input.
% a = A float [0.01,0.10] representing the target quantile.
% k = An integer [1,60] representing the target lag.
% cis = A float (0.0,0.1] representing the significance level of confidence intervals (optional, default=0.050).
% cib = An integer [10,1000] representing the number of bootstrap iterations of confidence intervals (optional, default=100).
%
% [OUTPUT]
% cq = A float (-Inf,Inf) representing the cross-quantilogram.
% ci = A row vector of floats (-Inf,Inf) of length 2 representing the lower and upper confidence intervals.
%
% [NOTES]
% The model computes partial cross-quantilograms when n is greater than 2 using exogenous variables from 3 to n.

function [cq,ci] = cross_quantilograms_sb(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('data',@(x)validateattributes(x,{'double'},{'real' 'finite' '2d' 'nonempty'}));
        ip.addRequired('a',@(x)validateattributes(x,{'double'},{'real' 'finite' '>=' 0.01 '<=' 0.10 'scalar'}));
        ip.addRequired('k',@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 1 '<=' 60 'scalar'}));
        ip.addOptional('cis',0.050,@(x)validateattributes(x,{'double'},{'real' 'finite' '>' 0 '<=' 0.1 'scalar'}));
        ip.addOptional('cib',100,@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 10 '<=' 1000 'scalar'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    data = validate_input(ipr.data);
    a = ipr.a;
    k = ipr.k;
    cis = ipr.cis;
    cib = ipr.cib;

    nargoutchk(2,2);

    [cq,ci] = cross_quantilograms_sb_internal(data,a,k,cis,cib);

end

function [cq,ci] = cross_quantilograms_sb_internal(data,a,k,cis,cib)

    [t,n] = size(data);
    len = t - k;
    partial = n > 2;

    cis = cis / 2;

    d = zeros(len,n);
    d(:,1) = data(k+1:t,1);
    d(:,2:n) = data(1:len,2:n);

    block_length = ppw_optimal_block_length(d);
    g = mean(block_length(:,1));

    a_sb = ones(len,n) .* a;
    cq_sb = zeros(cib,1);

    if (partial)
        for i = 1:cib
            indices = indices_bootstrap(len,g);

            d_sb = d(indices,:);
            q_sb = (d_sb <= repmat(gumbel_quantile(d_sb,a),len,1)) - a_sb;

            h_sb = q_sb.' * q_sb;

            if (det(h_sb) <= 1e-08)
                hi_sb = pinv(h_sb);
            else
                hi_sb = inv(h_sb);
            end

            cq_sb(i) = -hi_sb(1,2) / sqrt(hi_sb(1,1) * hi_sb(2,2));
        end
    else
        for i = 1:cib
            indices = indices_bootstrap(len,g);

            d_sb = d(indices,:);
            q_sb = (d_sb <= repmat(gumbel_quantile(d_sb,a),len,1)) - a_sb;

            h_sb = q_sb.' * q_sb;

            cq_sb(i) = h_sb(1,2) / sqrt(h_sb(1,1) * h_sb(2,2));
        end
    end

    q = (data <= repmat(gumbel_quantile(data,a),t,1)) - (ones(t,n) .* a);

    d = zeros(len,n);
    d(:,1) = q(k+1:t,1);
    d(:,2:n) = q(1:len,2:n);

    h = d.' * d;

    if (partial)
        if (det(h) <= 1e-08)
            hi = pinv(h);
        else
            hi = inv(h);
        end

        cq = -hi(1,2) / sqrt(hi(1,1) * hi(2,2));
    else
        cq = h(1,2) / sqrt(h(1,1) * h(2,2));
    end

    cqc = cq_sb - cq;
    ci = [min(0,gumbel_quantile(cqc,cis)) max(0,gumbel_quantile(cqc,1 - cis))];

end

function q = gumbel_quantile(x,p)

    index = 1 + ((size(x,1) - 1) * p);
    low = floor(index);
    high = ceil(index);

    x = sort(x);
    x_low = x(low,:);
    x_high = x(high,:);

    h = max(index - low,0);
    q = (h .* x_high) + ((1 - h) .* x_low);

end

function indices = indices_bootstrap(n,g)

    indices = [ceil(n * rand()); zeros(n - 1,1)];

    u = rand(n,1) < g;
    indices(u) = ceil(n .* rand(sum(u),1));

    zi = find(~u(2:n));
    indices(zi + 1) = indices(zi) + 1;

    fi = indices > n;
    indices(fi) = indices(fi) - n;

end

function bl = ppw_optimal_block_length(x)

    [t,n] = size(x);

    k = max(sqrt(log10(t)),5);
    c = 2 * sqrt(log10(t) / t);

    b_max = ceil(min(3 * sqrt(t),t / 3));
    m_max = ceil(sqrt(t)) + k;

    bl = zeros(n,2);

    for ppw_i = 1:n
        x_i = x(:,ppw_i);

        p1 = m_lag(x_i,m_max);
        p1 = p1(m_max+1:end,:);
        p1 = corr([x_i(m_max+1:end) p1]);
        p1 = p1(2:end,1);

        p2 = [m_lag(p1,k).' p1(end-k+1:end)];
        p2 = p2(:,k+1:end);
        p2 = sum((abs(p2) < (ones(k,m_max - k + 1) .* c))).';

        p3 = [(1:length(p2)).' p2];
        p3 = p3(p2 == k,:);

        if (isempty(p3))
            m_hat = find(abs(p1) > c,1,'last');
        else
            m_hat = p3(1,1);
        end

        m = min(2 * m_hat,m_max);

        if (m > 0)
            mm = (-m:m).';

            p1 = m_lag(x_i,m);
            p1 = p1(m+1:end,:);
            p1 = cov([x_i(m+1:end),p1]);

            act = sortrows([-(1:m).' p1(2:end,1)],1);
            ac = [act(:,2); p1(:,1)];

            mmn = mm ./ m;
            kernel_weights = ((abs(mmn) >= 0) .* (abs(mmn) < 0.5)) + (2 .* (1 - abs(mmn)) .* (abs(mmn) >= 0.5) .* (abs(mmn) <= 1));

            acw = kernel_weights .* ac;
            acw_ss = sum(acw)^2;

            g_hat = sum(acw .* abs(mm));
            dcb_hat = (4/3) * acw_ss;
            dsb_hat = 2 * acw_ss;

            b_comp1 = 2 * g_hat^2;
            b_comp2 = t^(1 / 3);
            bl_vl = min((b_comp1 / dsb_hat)^(1/3) * b_comp2,b_max);
            bl_cb = min((b_comp1 / dcb_hat)^(1/3) * b_comp2,b_max);

            bl(ppw_i,:) = [bl_vl bl_cb];
        else
            bl(ppw_i,:) = 1;
        end
    end

    function l = m_lag(x,n)

        mn = numel(x);
        l = ones(mn,n);

        for ml_i = 1:n
            l(ml_i+1:mn,ml_i) = x(1:mn-ml_i,1);
        end

    end

end

function data = validate_input(data)

    n = size(data,2);

    if (n < 2)
        error('The value of ''data'' is invalid. Expected input to be a matrix with at least 2 columns.');
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
%INPAINT2  The function implements different single-image inpainting algorithms
%
%     dst = cv.inpaint2(src, mask)
%     dst = cv.inpaint2(src, mask, 'OptionName', optionValue, ...)
%
% ## Input
% * __src__ source image, it could be of any type (8/16/32-bit integers or
%   32/64-bit floating-points) and any number of channels from 1 to 4. In case
%   of 3- and 4-channels images the function expect them in CIELab colorspace
%   or similar one, where first color component shows intensity, while second
%   and third shows colors. Nonetheless you can try any colorspaces.
% * __mask__ mask (8-bit 1-channel of same size as `src`), where non-zero
%   pixels indicate valid image area, while zero pixels indicate area to be
%   inpainted.
%
% ## Output
% * __dst__ Output image with the same size and type as `src`.
%
% ## Options
% * __Method__ Inpainting algorithms, one of:
%   * __ShiftMap__ (default) This algorithm searches for dominant
%     correspondences (transformations) of image patches and tries to
%     seamlessly fill-in the area to be inpainted using this transformations.
%
% The function reconstructs the selected image area from known area.
% See the original paper [He2012] for details.
%
% ## References
% [He2012]:
% > Kaiming He, Jian Sun. "Statistics of patch offsets for image completion".
% > In Computer Vision-ECCV 2012, pages 16-29. Springer, 2012.
%
% See also: cv.inpaint
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function MDDot=rotMatSecondDeriv2D(theta,thetaDot,thetaDDot,isClockwise)
%%ROTMATSECONDDERIV2D Given an initial angle of rotation, and the first and
%   second derivatives of the angle with respect to time, find the
%   derivative of the 2D rotation matrix. This is the second time
%   derivative of rotMat2D assuming that the rotation angle changes as a
%   function of time.
%
%INPUTS: theta The initial rotation angle radians about which the
%              time-varying 2D rotation matrix would rotate a 2X1 vector at
%              the current time.
%     thetaDot The derivative of theta with respect to time.
%    thetaDDot The second derivative of theta with respect to time.
%  isClockwise An optional boolean parameter indication whether the
%              rotation implied by theta is clockwise. The default if
%              omitted or an empty matrix is passed is false.
%
%OUTPUTS: MDDot The 2X2 second derivative of the rotation matrix (as one
%               could obtain using rotMat2D) with respect to time.
%
%EXAMPLE:
%This example shows that this function is consistent with numeric
%differentiation. The relative error of the finite differenced solution
%typically implies for 6-8 digits of precision.
% theta=2*pi*rand(1);
% thetaDot=randn(1);
% thetaDDot=randn(1);
% thetaState=[theta;thetaDot;thetaDDot];
% isClockwise=false;
% MDDot=rotMatSecondDeriv2D(theta,thetaDot,thetaDDot,isClockwise);
%
% MDot=rotMatDeriv2D(theta,thetaDot,isClockwise);
% %Propagate deltaT into the future for finite differencing. This is just a
% %second order linear model applied to the thetaState.
% deltaT=1e-8;
% order=2;
% F=FPolyKal(deltaT,3,order);
% thetaStateDelta=F*thetaState;
% MDotDelta=rotMatDeriv2D(thetaStateDelta(1),thetaStateDelta(2),isClockwise);
% MDDotNumDiff=(MDotDelta-MDot)/deltaT;
% RelErr=max(abs((MDDotNumDiff(:)-MDDot(:))./MDDot(:)))
%
%December 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(isClockwise))
    isClockwise=false;
end

sinTheta=sin(theta);
cosTheta=cos(theta);

if(isClockwise)
    MDDot=[-cosTheta*thetaDot^2-sinTheta*thetaDDot,-sinTheta*thetaDot^2+cosTheta*thetaDDot;
            sinTheta*thetaDot^2-cosTheta*thetaDDot,-cosTheta*thetaDot^2-sinTheta*thetaDDot];
else
    MDDot=[-cosTheta*thetaDot^2-sinTheta*thetaDDot, sinTheta*thetaDot^2-cosTheta*thetaDDot;
           -sinTheta*thetaDot^2+cosTheta*thetaDDot,-cosTheta*thetaDot^2-sinTheta*thetaDDot];
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 [status, results] = mrtrix_dwi2tensor (in_file, out_file, b_file, verbose)

%
% Calculate diffusion tensors.
%
% Parameters
% ----------
% in_file: The name of a diffusion file in .mif format
% out_file: The name of the resulting dti file in .mif format
% b_file: The name of a mrtrix format gradient vector file (default: 'encoding.b')
%
% Returns
% -------
% status: whether (0) or not (1) the operation succeeded
% results: the results of the operation in the terminal
%
% Notes
% -----
% http://www.brain.org.au/software/mrtrix/tractography/preprocess.html

if notDefined('verbose')
   verbose = true;
end

if notDefined('bkgrnd')
   bkgrnd = false;
end

if notDefined('b_file')
    b_file = 'encoding.b';
end

% This command generates  tensors:
cmd_str = sprintf('dwi2tensor %s %s -grad %s',in_file, out_file, b_file);

% Send it to mrtrix:
[status,results] = mrtrix_cmd(cmd_str, bkgrnd, verbose);

author: vistalab

repo: vistasoft


text:
function components_test02 ( )

%*****************************************************************************80
%
%% COMPONENTS_TEST02 tests I4MAT_COMPONENTS on a simple case.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 March 2011
%
%  Author:
%
%    John Burkardt
%
  m = 9;
  n = 17;

  a = [ ...
    0, 0, 0, 0, 0, 0, 0, 0, 0; ...
    0, 0, 1, 0, 0, 1, 0, 0, 0; ...
    0, 1, 1, 0, 1, 1, 1, 0, 0; ...
    0, 1, 1, 1, 1, 1, 1, 0, 0; ...
    0, 0, 1, 1, 1, 0, 0, 0, 0; ...
    0, 0, 1, 1, 1, 0, 0, 0, 0; ...
    0, 1, 1, 1, 0, 1, 0, 1, 0; ...
    0, 1, 1, 0, 0, 1, 0, 1, 0; ...
    0, 0, 1, 0, 0, 0, 0, 1, 0; ...
    0, 0, 0, 0, 1, 0, 1, 1, 0; ...
    0, 1, 0, 1, 1, 0, 1, 0, 0; ...
    0, 1, 1, 1, 1, 1, 0, 0, 0; ...
    0, 0, 1, 1, 0, 1, 0, 1, 0; ...
    0, 0, 1, 1, 0, 1, 0, 1, 0; ...
    0, 1, 1, 0, 1, 0, 1, 1, 0; ...
    0, 1, 0, 0, 1, 0, 1, 1, 0; ...
    0, 0, 0, 0, 0, 0, 0, 0, 0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'COMPONENTS_TEST02\n' );
  fprintf ( 1, '  I4MAT_COMPONENTS finds and labels connected\n' );
  fprintf ( 1, '  components in a 2D integer array.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  A:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : m
    fprintf ( 1, '    ' );
    for j = 1 : n
      fprintf ( 1, '%d', a(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  [ component_num, c ] = i4mat_components ( m, n, a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of components = %d\n', component_num );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  C:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : m
    fprintf ( 1, '    ' );
    for j = 1 : n
      fprintf ( 1, '%d', c(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
%GRABFG Framegrabber grabbing image function.
%
%SYNTAX imgmatrix = grabfg(devinfo) or
%       moviematrix = grabfg(devinfo,imgmatrix) or
%       [imgmatrix, moviematrix] = grabfg(devinfo) or
%       grabfg(devinfo, imgmatrix) or
%       grabfg(devinfo,imgmatrix,moviematrix) or
%       grabfg(devinfo,imgmatrix,moviematrix, scaling factor)
%
%DESCRIPTION This function will let the user grab an
%image from a opened framegrabber device. The user has
%to give the framegrabbers device information structure
%as input argument. Executing the openfg function has
%retrieved this device information structure
%
%INPUT The device information structure (devinfo) of
%an opened framegrabber device and the eventually the
%preallocated image matrix  (imgmatrix);
%
%REMARKS It will be verified if the input is a valid
%device information structure. If it's valid, it will
%try to get the deviceid from the device information
%structure. If the second input is given, the function
%will check the size of the second input and checks if
%there's enough room to store the grabbed image. If
%there's enough room the variable (of type imgmatrix)
%will be overwritten without creating a new matrix,
%else a new matrix will be created and stored in that
%variable.
%
%OUTPUT If the input isn't valid or the framegrabber
%isn't able to grab an image, some error message will
%be shown to the user. If everything's ok it will return
%the grabbed image stored in a Matlab matrix (imgmatrix).
%
%The image is stored in transposed way. To show the image
%correctly in Matlab you can use the accent ' on the
%imgmatrix to correct it (=> image(imgmatrix') )
%
%Grabfg also generates true matlab movie matrices. The
%resolution reduction factor defaults to 2.

%EXAMPLE There are two ways to use this function. Both
%will be shown here. Variable m is the device info
%structure which was returned from the function openfg.
%
%1. Just grab an image and store it in the variable im,
%   type:
%               im = grabfg(m);
%
%   To show this image in Matlab
%   (in the right way, not transposed) type:
%               image(im');
%
%2. Just grab an image and movie matrix and store it in the
%   variable im and movie, type:
%               [im,movie] = grabfg(m);
%
%   To show the movie matrix in Matlab
%   (the movie matrix does not need to be transposed) type:
%               movie(movie);
%
%3. First create a matrix of the appropriate size and
%   type where the image will be stored:
%               im = uint8(zeros(768,576));
%
%   Then grab the image:
%               grabfg(m,im);
%
%   To show the image in Matlab type:
%               image(im');

disp('Error: grabfg not found.')


text:
function g = netderiv(w, net, x)
%NETDERIV Evaluate derivatives of network outputs by weights generically.
%
%	Description
%
%	G = NETDERIV(W, NET, X) takes a weight vector W and a network data
%	structure NET, together with the matrix X of input vectors, and
%	returns the gradient of the outputs with respect to the weights
%	evaluated at W.
%
%	See also
%	NETEVFWD, NETOPT
%

%	Copyright (c) Ian T Nabney (1996-2001)

fstr = [net.type, 'deriv'];
net = netunpak(net, w);
g = feval(fstr, net, x);

author: bayesnet

repo: bnt


text:
% Constructing a random graph based on a given degree sequence.
% Idea source: Molloy M. & Reed, B. (1995) Random Structures and Algorithms 6, 161-179
%
% INPUTs: a graphic sequence of numbers, 1xn
% OUTPUTs: adjacency matrix of resulting graph, nxn
%
% Note: The simple version of this algorithm gets stuck about half
%       of the time, so in this implementation the last problematic
%       edge is rewired.
%
% Other routines used: adj2edgeL.m, rewireThisEdge.m, edgeL2adj.m
% GB: last updated, Oct 25 2012

function adj= randomGraphFromDegreeSequence(Nseq)

stubs=Nseq;                % assign degrees to stubs
adj = zeros(length(Nseq)); % initialize adjacency matrix

old_sum = 0;
cnt=0;

while sum(stubs)>0   % while no more stubs are left to connect

  if cnt>5                       % if rewiring did not work 5 times

    el = adj2edgeL(adj);
    ind = find(stubs>0);

    if length(ind) == 1; elr = rewireThisEdge([el; ind(1) ind(1) 1],ind(1),ind(1));  end
    if length(ind) == 2;  elr = rewireThisEdge([el; ind(1) ind(2) 1; ind(2) ind(1) 1],ind(1),ind(2)); end

    if length(ind)>2 || isempty(elr)           % restart algorithm
      printf('randomGraphFromDegreeSequence(): restarting ...\n')
      stubs = Nseq;
      adj = zeros(length(Nseq));
      old_sum = 0;
      cnt=0;

    else
      adj = edgeL2adj(elr);  % return matrix with last edge rewired
      return

    end

  end

  new_sum = sum(stubs);

  if old_sum==new_sum; cnt = cnt+1; end       % no new nodes have been connected, counter+1
  if old_sum~=new_sum; cnt=0; end             % new connections, restart count

  [~,n1] = max(stubs);                % pick the node with highest number of remaining stubs

  old_sum = sum(stubs);

  ind = find(stubs>0);
  n2 = ind(randi(length(ind)));

  if n1==n2; continue; end            % no self-loops

  if adj(n1,n2)>0; continue; end      % no double edges
  adj(n1,n2)=1; adj(n2,n1)=1;
  stubs(n1) = stubs(n1) - 1;
  stubs(n2) = stubs(n2) - 1;

end

author: aeolianine

repo: octave-networks-toolbox