text:
function area = voronoi_areas_direct ( d_num, d_xyz, face_num, face_d, v_xyz )

%*****************************************************************************80
%
%% VORONOI_AREAS_DIRECT computes the area of each polygon in a Voronoi diagram.
%
%  Discussion:
%
%    This algorithm works directly from the Delaunay triangulation and the
%    Voronoi vertex locations, without having to determine the explicit
%    description of each Voronoi polygon as a list of Voronoi vertices.
%    This means this routine should be significantly faster than
%    VORONOI_AREAS, at least for larger problems.
%
%
%    The Delaunay triangle has vertices D1, D2, D3, and a Voronoi vertex
%    V which we imagine is inside the triangle, although in fact it need
%    not be.
%
%    Include now the Delaunay vertex averages, D12, D23, and D31, and
%    the Delaunay triangle can be divided into 3 pieces, contributing to
%    the Voronoi polygons associated with D1, D2 and D3 as follows:
%
%      D1 gets D1:V:D31 + D1:D12:V
%      D2 gets D2:V:D12 + D2:D23:V
%      D3 gets D3:V:D23 + D3:D31:V
%
%    The advantage of computing things this way is that it can be done
%    directly from the Delaunay triangle information.  If we must first
%    determine the explicit list of Voronoi vertices that form a particular
%    Voronoi polygon, then compute the area, it is cumbersome and expensive
%    to determine this connectivity information.
%
%    Presumably, if things have been done correctly, the cases where
%    the Voronoi vertex is outside of the Delaunay triangle will correspond
%    to situations where some of the areas come out negative, and
%    sum to the correct final result.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 May 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer D_NUM, the number of data points and Voronoi polygons.
%
%    Input, real D_XYZ(3,D_NUM), the coodinates of the data points.
%
%    Input, integer FACE_NUM, the number of faces, and Voronoi vertices.
%
%    Input, integer FACE_D(3,FACE_NUM), the data points that form
%    each triangle.
%
%    Input, real V_XYZ(3,FACE_NUM), the coordinates of the Voronoi vertices.
%
%    Output, real AREA(D_NUM), the area of each Voronoi polygon.
%
  area = zeros ( d_num, 1 );

  r = 1.0;
%
%  Consider each triangle.
%
  for face = 1 : face_num
%
%  Define some things.
%
    p1 = face_d(1,face);
    p2 = face_d(2,face);
    p3 = face_d(3,face);

    d1 = d_xyz ( 1:3, p1 );
    d2 = d_xyz ( 1:3, p2 );
    d3 = d_xyz ( 1:3, p3 );

    d12 = 0.5 * ( d1 + d2 );
    d23 = 0.5 * ( d2 + d3 );
    d31 = 0.5 * ( d3 + d1 );
    v = v_xyz(1:3,face);
%
%  Force all points to lie on the sphere.
%
    d12 = d12 / norm ( d12 );
    d23 = d23 / norm ( d23 );
    d31 = d31 / norm ( d31 );
    v = v / norm ( v );
%
%  Compute orientation and area of the 6 subtriangles.
%
    o1a = stri_vertices_to_orientation ( d1, v, d31 );
    a1a = stri_vertices_to_area ( r, d1, v, d31 );
    o1b = stri_vertices_to_orientation ( d1, d12, v );
    a1b = stri_vertices_to_area ( r, d1, d12, v );

    o2a = stri_vertices_to_orientation ( d2, v, d12 );
    a2a = stri_vertices_to_area ( r, d2, v, d12 );
    o2b = stri_vertices_to_orientation ( d2, d23, v );
    a2b = stri_vertices_to_area ( r, d2, d23, v );

    o3a = stri_vertices_to_orientation ( d3, v, d23 );
    a3a = stri_vertices_to_area ( r, d3, v, d23 );
    o3b = stri_vertices_to_orientation ( d3, d31, v );
    a3b = stri_vertices_to_area ( r, d3, d31, v );
%
%  Contribute to the Voronoi areas.
%
    area(p1) = area(p1) + o1a * a1a + o1b * a1b;
    area(p2) = area(p2) + o2a * a2a + o2b * a2b;
    area(p3) = area(p3) + o3a * a3a + o3b * a3b;

  end

  return
end


text:
function r8vec_sort_heap_index_d_test ( )

%*****************************************************************************80
%
%% R8VEC_SORT_HEAP_INDEX_D_TEST tests R8VEC_SORT_HEAP_INDEX_D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  n = 20;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_SORT_HEAP_INDEX_D_TEST\n' );
  fprintf ( 1, '  R8VEC_SORT_HEAP_INDEX_D creates a descending\n' );
  fprintf ( 1, '  sort index for an R8VEC.\n' );

  b = 0.0;
  c = 3.0 * n;
  seed = 123456789;

  [ a, seed ] = r8vec_uniform_ab ( n, b, c, seed );

  r8vec_print_some ( n, a, 1, 10, '  Unsorted array:' );

  indx = r8vec_sort_heap_index_d ( n, a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  After indexed descending sort:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I, INDX(I), A(I)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %6d  %6d  %12f\n', i, indx(i), a(i) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now use the index array to carry out the\n' );
  fprintf ( 1, '  permutation implicitly.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  INDX(I), ARRAY(INDX(I))\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %6d  %12f\n', indx(i), a(indx(i)) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function b = hp(fs)

% Equiripple Highpass filter designed using the FIRPM function.
% All frequency values are in Hz.

Fstop = 0.125;           % Stopband Frequency
Fpass = 0.25;            % Passband Frequency
Dstop = 0.001;           % Stopband Attenuation
Dpass = 0.057501127785;  % Passband Ripple
dens  = 20;              % Density Factor

% Calculate the order from the parameters using FIRPMORD.
[N, Fo, Ao, W] = firpmord([Fstop, Fpass]/(fs/2), [0 1], [Dstop, Dpass]);

% Calculate the coefficients using the FIRPM function.
b  = firpm(N, Fo, Ao, W, {dens});

end

author: nesl


text:
function fn = fieldnames(item)

% function fn = fieldnames(item)
% Return a list of all (inherited and non-inherited) field names.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: fieldnames.m 4867 2012-08-30 13:04:51Z volkmar $

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

fn1 = fieldnames(item.cfg_item);
fn2 = mysubs_fields;

fn = [fn1(:); fn2(:)]';

author: spm


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%

% problem 2- graph of
% u(t+1)*u(t-1)

 t1=-5:.1:1;
 t2=1:.1:10;
 x1=zeros(size(t1));
 x2=ones(size(t2));
 t=[t1 t2];
 x=[x1 x2];
 plot(t,x);
 ylim([-0.1 1.1]);

 figure
 t=-5:.1:10;
 x=heaviside(t+1).*heaviside(t-1) ;
 plot(t,x);
 ylim([-0.1 1.1]);

author: Sable

repo: mcbench-benchmarks


text:
function addQuitMenuItem()
    % ADDQUITMENUITEM add quit menu to main file menu
    hQuit=findall(gcf,'Label','Quit Zmap');
    if isempty(hQuit)
        mainfile=findall(gcf,'Tag','figMenuFile');
        uimenu(mainfile,'Label','Quit Zmap','Separator','on','MenuSelectedFcn',@(~,~)restartZmap);
    end
end

author: CelsoReyes


text:
function ori = discreteSample(SO3F,npoints,varargin)
% draw a random sample
%

% spread points over different centers
ic = discretesample([SO3F.weights(:);SO3F.c0], npoints).';

isUniform = ic == length(SO3F.weights)+1;

% some uniform random orientations
ori = orientation.rand(length(ic),SO3F.CS,SO3F.SS);

% the remaining points
npoints = nnz(~isUniform);
if npoints == 0, return; end

% take random rotational axes for the remaining samples
axis = vector3d.rand(npoints);

% random rotational angles
M = 1000000; % discretisation parameter

hw = min(4*SO3F.psi.halfwidth,90*degree);
t = linspace(cos(hw),1,M);

% compute cummulative distribution function
c = 4 / pi * cumsum(sqrt(1-t.^2) .* SO3F.psi.eval(t)) / M;
c = c ./ c(end);

% random sample with respect to the CDF
r = rand(npoints,1);
[~,id] = histc(r,c);
angle = 2 * acos(t(id)).';

% set up random orientations
ori(~isUniform) = times(reshape(SO3F.center(ic(~isUniform)),[],1), ...
  rotation.byAxisAngle(axis,angle),false);

% random symmetry elements
ori = ori .* SO3F.CS.rot(randi(SO3F.CS.numSym,npoints,1));
if SO3F.SS.numSym>1
  ori = SO3F.SS.rot(randi(SO3F.SS.numSym,npoints,1)) .* ori;
end

author: mtex-toolbox

repo: mtex

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


text:
% textsc() - places text in screen coordinates and places
%            a title at the top of the figure.
%
% Usage:
%   H = textsc(X,Y,TXT) places the text string, TXT
%   at the normalized coordinates X and Y.  H is the
%   handle to the text object.
%
%   H = textsc(TXT,'title') places a title at the top
%   of the figure window.  This is useful when you
%   want to place a single title above multiple
%   subplots.
%
% Notes: textsc creates an invisible AXES which occupies
% the entire FIGURE.  The units of the AXES are
% normalized (range from 0 to 1).  textsc checks
% all the children of the current FIGURE to determine
% if an AXES exist which meets these criteria already
% exist.  If one does, then it places the text relative
% to that AXES.
%
% Author: John L. Galenski, January 21, 1994

% Written by John L. Galenski III
% All Rights Reserved  January 21, 1994
% LDM031695jlg

function H = textsc(x,y,txt);

% Basic error checking
if nargin < 2
    y = 'title';
end

% Check to see if AXES already exist
ch = get(gcf,'Children');
ax = findobj(gcf,'Type','axes','Tag','TEXTSC');
if isempty(ax)
  ax = axes('Units','Normal','Position',[0 0 1 1], ...
            'Visible','Off','Tag','TEXTSC');
else
  axes(ax)
end

% Place the text
if isstr(y) && isstr(x) && strcmp(lower(y),'title')  % Subplot title
  txt = x;
  x = .5;
  tmp = text('Units','normal','String','tmp','Position',[0 0 0]);
  ext = get(tmp,'Extent');
  delete(tmp)
  H = ext(4);
  y = 1 - .60*H;
end
h = text(x,y,txt,'VerticalAlignment','Middle', ...
         'HorizontalAlignment','Center','interpreter', 'none' );

% Make the original AXES current
if ~isempty(ch)
  set(gcf,'CurrentAxes',min(ch))
end

% Check for output
if nargout == 1
  H = h;
end

author: PatternRecognition

repo: OpenBMI


text:
function J = image2Cell(I)

    numChannels = size(I,3);
    for channelIdx = 1:numChannels
        J{channelIdx} = squeeze(I(:,:,channelIdx));
    end

author: thomas-koehler

repo: SupER


text:
function dtr = sv_clock_bias(t, toc, a0, a1, a2, e, sqrtA, toe, Delta_n, M0)
% 输入:
% a0 a1 a2 toc: 卫星时钟校正模型方程中3个参数，  toc: 第一数据块参考时间, 被用作时钟校正模型中时间参考点

% dtr： 卫星时钟偏差

dtr = a0+a1*(t-toc)+a2*(t-toc).^2;%t为未做钟差改正的观测时刻

F = -4.442807633e-10;
mu = 3.986005e14;
A = sqrtA^2;
cmm = sqrt(mu/A^3); % computed mean motion
tk = t - toe;
% account for beginning or end of week crossover
if (tk > 302400)
    tk = tk-604800;
end
if (tk < -302400)
    tk = tk+604800;
end
% apply mean motion correction
n = cmm + Delta_n;

% Mean anomaly
mk = M0 + n*tk;

% solve for eccentric anomaly
Ek = mk;
Ek = mk + e*sin(Ek);
Ek = mk + e*sin(Ek);
Ek = mk + e*sin(Ek);

% dsv 时间为s
dtr = dtr + F*e*sqrtA*sin(Ek);

end

author: yandld


text:
classdef NMPCDCBF2 < handle
    % MPC with distance constraints
    properties
        system
        params
        x0
        x_curr
        time_curr = 0.0
        xlog = []
        ulog = []
        tt = 0
        distlog = []
        solvertime = []
        xopenloop = {}
        uopenloop = {}
        u_cost = 0
        obs
    end
    methods
        function self = NMPCDCBF2(x0, system, params)
            % Define MPC_CBF controller
            self.x0 = x0;
            self.x_curr = x0;
            self.system = system;
            self.params = params;
        end

        function sim(self, time)
            % Simulate the system until a given time
            tic;
            xk = self.x_curr;
            while self.time_curr < time
                [~, uk, tt] = self.solveNMPCDCBF1(self.x_curr);%if infeasiblility happens, stop
                if tt == -1
                    self.tt = tt;%record computing time for each time-step
                    return
                end
                xk = self.system.A * xk + [xk(4,1) * cos(xk(3,1))*self.system.dt;xk(4,1) * sin(xk(3,1))*self.system.dt;0;0] + self.system.C * uk;
                % update system
                self.x_curr = xk;
                self.time_curr = self.time_curr + self.system.dt;
                self.xlog = [self.xlog, xk];
                self.ulog = [self.ulog, uk];
                self.u_cost = self.u_cost + uk'*uk;
                tt = tt + toc;
                self.tt = tt;
                return
            end

        end

        function [xopt, uopt,tt] = solveNMPCDCBF1(self, xk)
            % Solve NMPC-DCBF
            [feas, x, u, J] = self.solve_cftoc1(xk);
            if ~feas
                xopt = [];
                uopt = [];
                tt = -1;
                return
            else
                xopt = x(:,2);
                uopt = u(:,1);
                tt = 0;
            end
        end

        function [feas, xopt, uopt, Jopt] = solve_cftoc1(self, xk)
            % Solve CFTOC
            % extract variables
            N = self.params.N;
            % define variables and cost
            x = sdpvar(4, N+1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            u = sdpvar(4, N);
            constraints = [];
            cost = 0;
            % initial constraint
            constraints = [constraints; x(:,1) == xk];
            % add constraints and costs
            AA=self.system.A;
            BB=self.system.B;
            CC=self.system.C;
            constraints = [constraints;
                    self.system.xl <= x(:,1) <= self.system.xu;
                    self.system.ul <= u(:,1) <= self.system.uu;
                    x(:,2) == AA * x(:,1) + [x(4,1)*cos(x(3,1))*self.system.dt;x(4,1)*sin(x(3,1))*self.system.dt;0; 0] + CC * u(:,1)];
                cost = cost + (x(:,1)-[3;0.01;0;0])'*self.params.Q*(x(:,1)-[3;0.01;0;0]) + (u(:,1)-[0;0;1;1])'*self.params.R*(u(:,1)-[0;0;1;1]); % self.params.S*(w(:,1)-1)^2;
            for i = 2:1:N
                constraints = [constraints;
                    self.system.xl <= x(:,i) <= self.system.xu;
                    self.system.ul <= u(:,i) <= self.system.uu;
                    x(:,i+1) == AA * x(:,i) + [x(4,i)*cos(x(3,i))*self.system.dt;x(4,i)*sin(x(3,i))*self.system.dt;0; 0] + CC * u(:,i)];
                cost = cost + (x(:,i)-[3;0.01;0;0])'*self.params.Q*(x(:,i)-[3;0.01;0;0]) + (u(:,i)-[0;0;1;1])'*self.params.R*(u(:,i)-[0;0;1;1]);% self.params.S*(w(:,i)-1)^2;
            end
            % add CBF constraints
            for i = 1:N-1
                pos = self.obs.pos1;
                r = self.obs.r1 ;
                b = (x([1:2],i)-pos)'*((x([1:2],i)-pos))-r^2;
                b_next = (x([1:2],i+1)-pos)'*((x([1:2],i+1)-pos))-r^2;
                b_next_next = (x([1:2],i+2)-pos)'*((x([1:2],i+2)-pos))-r^2;
                b1 = (b_next - b)/self.system.dt + self.params.gamma1/self.system.dt * (b);
                b1_next = (b_next_next - b_next)/self.system.dt + self.params.gamma1/self.system.dt * (b_next);
                constraints = [constraints; b1_next  >=u(4,i)*(1-self.params.gamma2)*b1;b_next  >= u(3,i)*(1-self.params.gamma1)*b];%Highest order mcbf=2
            end
            % add terminal cost
            cost = cost + (x(:,N+1)-[3;0.01;0;0])'*self.params.P*(x(:,N+1)-[3;0.01;0;0]);
            ops = sdpsettings('solver','ipopt','verbose',0);
            % solve optimization
            diagnostics = optimize(constraints, cost, ops);
            if diagnostics.problem == 0
                feas = true;
                xopt = value(x);
                uopt = value(u);
                Jopt = value(cost);
            else
                feas = false;
                xopt = [];
                uopt = [];
                Jopt = [];
            end
            pos = self.obs.pos1;
            r = self.obs.r1 ;
            self.distlog = [self.distlog, (r^2-(xk(1:2)-pos)'*(xk(1:2)-pos))];
            self.xopenloop{size(self.xopenloop,2)+1} = xopt;
            self.uopenloop{size(self.uopenloop,2)+1} = uopt;
            self.solvertime = [self.solvertime, diagnostics.solvertime];
            fprintf('solver time: %f\n', diagnostics.solvertime);
        end
    end
end

author: HybridRobotics

repo: NMPC-DCLF-DCBF


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

% CIFAR-10 model from
% M. Lin, Q. Chen, and S. Yan. Network in network. CoRR, abs/1312.4400, 2013.

net.layers = {} ;
b=0 ;

% Block 1
net.layers{end+1} = struct('type', 'conv', ...
                           'name', 'conv1', ...
                           'weights', {{0.01*randn(5,5,3,192,'single'), b*ones(1,192,'single')}}, ...
                           'learningRate', [.1 2], ...
                           'stride', 1, ...
                           'pad', 2) ;
net.layers{end+1} = struct('type', 'relu', 'name', 'relu1') ;
net.layers{end+1} = struct('type', 'conv', ...
                           'name', 'cccp1', ...
                           'weights', {{0.05*randn(1,1,192,160, 'single'), b*ones(1,160,'single')}}, ...
                           'learningRate', [.1 2], ...
                           'stride', 1, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'relu', 'name', 'relu_cccp1') ;
net.layers{end+1} = struct('type', 'conv', ...
                           'name', 'cccp2', ...
                           'weights', {{0.05*randn(1,1,160,96,'single'), b*ones(1,96,'single')}}, ...
                           'learningRate', [.1 2], ...
                           'stride', 1, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'relu', 'name', 'relu_cccp2') ;
net.layers{end+1} = struct('name', 'pool1', ...
                           'type', 'pool', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'dropout', 'name', 'dropout1', 'rate', 0.5) ;

% Block 2
net.layers{end+1} = struct('type', 'conv', ...
                           'name', 'conv2', ...
                           'weights', {{0.05*randn(5,5,96,192,'single'), b*ones(1,192,'single')}}, ...
                           'learningRate', [.1 2], ...
                           'stride', 1, ...
                           'pad', 2) ;
net.layers{end+1} = struct('type', 'relu', 'name', 'relu2') ;
net.layers{end+1} = struct('type', 'conv', ...
                           'name', 'cccp3', ...
                           'weights', {{0.05*randn(1,1,192,192,'single'), b*ones(1,192,'single')}}, ...
                           'learningRate', [.1 2], ...
                           'stride', 1, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'relu', 'name', 'relu_cccp3') ;
net.layers{end+1} = struct('type', 'conv', ...
                           'name', 'cccp4', ...
                           'weights', {{0.05*randn(1,1,192,192, 'single'), b*ones(1,192,'single')}}, ...
                           'learningRate', [.1 2], ...
                           'stride', 1, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'relu', 'name', 'relu_cccp4') ;
net.layers{end+1} = struct('name', 'pool2', ...
                           'type', 'pool', ...
                           'method', 'avg', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'dropout', 'name', 'dropout2', 'rate', 0.5) ;

% Block 3
net.layers{end+1} = struct('type', 'conv', ...
                           'name', 'conv3', ...
                           'weights', {{0.05*randn(3,3,192,192,'single'), b*ones(1, 192, 'single')}}, ...
                           'learningRate', [.1 2], ...
                           'stride', 1, ...
                           'pad', 1) ;
net.layers{end+1} = struct('type', 'relu', 'name', 'relu3') ;
net.layers{end+1} = struct('type', 'conv', ...
                           'name', 'cccp5', ...
                           'weights', {{0.05*randn(1,1,192,192,'single'), b*ones(1,192,'single')}}, ...
                           'learningRate', [.1 2], ...
                           'stride', 1, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'relu', 'name', 'relu_cccp5') ;
net.layers{end+1} = struct('type', 'conv', ...
                           'name', 'cccp6', ...
                           'weights', {{0.05*randn(1,1,192,10, 'single'), b*ones(1,10,'single')}}, ...
                           'learningRate', 0.1*[.1 2], ...
                           'stride', 1, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'relu', 'name', 'relu_cccp6') ;
net.layers{end+1} = struct('type', 'pool', ...
                           'name', 'pool3', ...
                           'method', 'avg', ...
                           'pool', [7 7], ...
                           '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.5*ones(1,30) 0.1*ones(1,10) 0.02*ones(1,5)]  ;
net.meta.trainOpts.weightDecay = 0.0005 ;
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:
function x = f2_abscissas ( n )

%*****************************************************************************80
%
%% F2_ABSCISSAS computes Fejer Type 2 abscissas.
%
%  Discussion:
%
%    The interval is [-1,+1].
%
%    The abscissas are the cosines of equally spaced angles.
%    The angles are computed as N+2 equally spaced values between 0 and PI,
%    but with the first and last angle omitted.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Philip Davis, Philip Rabinowitz,
%    Methods of Numerical Integration,
%    Second Edition,
%    Dover, 2007,
%    ISBN: 0486453391,
%    LC: QA299.3.D28.
%
%    Walter Gautschi,
%    Numerical Quadrature in the Presence of a Singularity,
%    SIAM Journal on Numerical Analysis,
%    Volume 4, Number 3, 1967, pages 357-362.
%
%    Joerg Waldvogel,
%    Fast Construction of the Fejer and Clenshaw-Curtis Quadrature Rules,
%    BIT Numerical Mathematics,
%    Volume 43, Number 1, 2003, pages 1-18.
%
%  Parameters:
%
%    Input, integer N, the order of the rule.
%
%    Output, real X(N), the abscissas.
%
  if ( n == 1 )
    x(1) = 0.0;
    return
  elseif ( n == 2 )
    x(1) = -0.5;
    x(2) =  0.5;
    return
  end

  for i = 1 : n
    theta(i) = ( n + 1 - i ) * pi ...
             / ( n + 1     );
  end

  x(1:n) = cos ( theta(1:n) );

  return
end


text:
function r = getResidual(SR, LR, W, photometricParams)
% GETRESIDUAL Get residual error for low-resolution and super-resolved
% data.
%   GETRESIDUAL computes the the residual error caused by a super-resolved
%   image with the associated low-resolution frames and the model
%   parameters (system matrix and photometric parameters).

    if nargin < 4 || isempty(photometricParams)
        r = LR - W*SR;
    else
        numFrames = size(photometricParams.mult, 3);
        numLRPixel = length(LR)/numFrames;
        if isvector(photometricParams.mult(:,:,1))
            bm = zeros(size(LR));
            ba = zeros(size(LR));
            for k = 1:numFrames
                bm( ((k-1)*numLRPixel + 1):(k*numLRPixel) ) = repmat(photometricParams.mult(k), numLRPixel, 1);
                ba( ((k-1)*numLRPixel + 1):(k*numLRPixel) ) = repmat(photometricParams.add(k), numLRPixel, 1);
            end
        else
            bm = zeros(size(LR));
            ba = zeros(size(LR));
            for k = 1:numFrames
                bm( ((k-1)*numLRPixel + 1):(k*numLRPixel) ) = imageToVector(photometricParams.mult(:,:,k));
                ba( ((k-1)*numLRPixel + 1):(k*numLRPixel) ) = imageToVector(photometricParams.add(:,:,k));
            end
        end
        r = LR - bm .* (W*SR) - ba;
    end


text:
%-fanDTasia ToolBox------------------------------------------------------------------
% This Matlab script is part of the fanDTasia ToolBox: a Matlab library for Diffusion
% Weighted MRI (DW-MRI) Processing, Diffusion Tensor (DTI) Estimation, High-order
% Diffusion Tensor Analysis, Diffusion Kurtosis Imaging (DKI) Estimation,
% Tensor ODF estimation, Visualization and more.
%
% A Matlab Tutorial on DW-MRI can be found in:
% http://www.cise.ufl.edu/~abarmpou/lab/fanDTasia/tutorial.php
%
%-CITATION---------------------------------------------------------------------------
% If you use this software please cite the following papers on 4th-order tensors:
% 1) A. Barmpoutis et al. "Diffusion Kurtosis Imaging: Robust Estimation from DW-MRI
%    using Homogeneous Polynomials", In the Proceedings of ISBI, 2011, pp. 262-265.
% 2) A. Barmpoutis and B.C. Vemuri, "A Unified Framework for Estimating Diffusion
%    Tensors of any order with Symmetric Positive-Definite Constraints",
%    In the Proceedings of ISBI, 2010, pp. 1385-1388.
%
%-DESCRIPTION------------------------------------------------------------------------
% This demo script shows how to compute the Diffusion Kurtosis Coefficients from a given
% Diffusion-Weighted MRI dataset. The method guarantees that the estimated diffusivity is
% positive semi-definite, and computes the diffusion and kurtosis tensors using the method
% in Sec. 4.1 of the ISBI'11 paper. Here the given demo dataset consists of 5 voxels,
% 30 gradient directions x 2 b-values from the real brain dataset used in the ISBI'11 paper.
%
%-USE--------------------------------------------------------------------------------
% [D,W]=DEMO_DKI_Estimation_Method1_v2;
%
% D: is a vector with the computed Unique Coefficients of the 2nd-order Diffusion Tensor
% W: is a vector with the computed Unique Coefficients of the 4th-order Kurtosis Tensor
%
%-DISCLAIMER-------------------------------------------------------------------------
% You can use this source code for non commercial research and educational purposes
% only without licensing fees and is provided without guarantee or warrantee expressed
% or implied. You cannot repost this file without prior written permission from the
% authors. If you use this software please cite the following papers:
% 1) A. Barmpoutis et al. "Diffusion Kurtosis Imaging: Robust Estimation from DW-MRI
%    using Homogeneous Polynomials", In the Proceedings of ISBI, 2011.
% 2) A. Barmpoutis and B.C. Vemuri, "A Unified Framework for Estimating Diffusion
%    Tensors of any order with Symmetric Positive-Definite Constraints",
%    In the Proceedings of ISBI, 2010, pp. 1385-1388.
%
%-AUTHOR-----------------------------------------------------------------------------
% Angelos Barmpoutis, PhD
% Digital Worlds Institute
% University of Florida, Gainesville, FL 32611, USA
% angelbar at ufl dot edu
%------------------------------------------------------------------------------------
function [DKI_D,DKI_W]=DEMO_DKI_Estimation_Method1_v2

%%% DATA OPENING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%open data files
[S,B]=openFDT('real_data_5voxels.fdt');

%S0 signal, no diffusion weighting
S0=S(:,:,:,1);

%acquisition shell 1, 30 orientations
S_1real=S(:,:,:,[2:31]);
GradientOrientations_1=B([2:31],[1:3]);
BValue_1=B(2,4);

%acquisition shell 2, 30 orientations
S_2real=S(:,:,:,[32:61]);
GradientOrientations_2=B([32:61],[1:3]);
BValue_2=B(32,4);

%%% OPTIONAL: ADD RICIAN NOISE TO THE DATA FOR QUANTITATIVE EVALUATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%

stdv=input('Do you want to add Rician noise to the data for validation? \n If yes, then give the Std. Dev. (e.g: 0.1), otherwise type 0.\n STD.DEV.:');
S_1=S_1real;
for i=1:size(S_1real,4)
    S_1(:,:,:,i)=sqrt((S_1real(:,:,:,i)+stdv*S0.*randn(size(S_1real,1),size(S_1real,2),size(S_1real,3))).^2+(stdv*S0.*randn(size(S_1real,1),size(S_1real,2),size(S_1real,3))).^2);
end
S_2=S_2real;
for i=1:size(S_2real,4)
    S_2(:,:,:,i)=sqrt((S_2real(:,:,:,i)+stdv*S0.*randn(size(S_2real,1),size(S_2real,2),size(S_2real,3))).^2+(stdv*S0.*randn(size(S_2real,1),size(S_2real,2),size(S_2real,3))).^2);
end
%your data can have as many shells and bvalues you want

%%% INITIALIZATION - COMPUTING AUXILIARY DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Construct a constant set of polynomial coefficients C
C_order4=constructSetOf321Polynomials(4)'; %computes C from section 5.1 (ISBI'10)
A_1=D4toDKImatrix(BValue_1); %Computes the matrix A from Table 1 (ISBI'11)
A_2=D4toDKImatrix(BValue_2); %Computes the matrix A from Table 1 (ISBI'11)

%shell 1
G_1_order2=constructMatrixOfMonomials(GradientOrientations_1, 2); %computes G from section 5.1 (ISBI'10)
G_1_order4=constructMatrixOfMonomials(GradientOrientations_1, 4); %computes G from section 5.1 (ISBI'10)
Gbig_1=[-BValue_1*G_1_order2 BValue_1*BValue_1/6*G_1_order4];  %all the monomials for orders 2 and 4.

%shell 2
G_2_order2=constructMatrixOfMonomials(GradientOrientations_2, 2); %computes G from section 5.1 (ISBI'10)
G_2_order4=constructMatrixOfMonomials(GradientOrientations_2, 4); %computes G from section 5.1 (ISBI'10)
Gbig_2=[-BValue_2*G_2_order2 BValue_2*BValue_2/6*G_2_order4];  %all the monomials for orders 2 and 4.

%your data can have as many shells and bvalues you want
Gbig=[Gbig_1; Gbig_2]; %all the monomials for orders 2 and 4 and bvalues b1 and b2.

%%%%%% MAIN LOOP - METHOD: LINEAR FITTING - NO CONSTRAINTS %%%%%%%%%%%%%%%%%%%%%%%%%
for x=1:size(S,1)
    for y=1:size(S,2)
        for z=1:size(S,3)

            logS_1=log(squeeze(S_1(x,y,z,:))/S0(x,y,z));
            logS_2=log(squeeze(S_2(x,y,z,:))/S0(x,y,z));

            %The following 2 steps implement the method in Sec. 4.1 of the ISBI'11 paper.
            %Step 1: Compute a positive-definite 4th-order tensor for each b-value
            D4_1=C_order4*lsqnonneg(-G_1_order4*C_order4, logS_1);%computes a positive-definite tensor for b1
            D4_2=C_order4*lsqnonneg(-G_2_order4*C_order4, logS_2);%computes a positive-definite tensor for b2

            %Step 2: Compute DKI from the positive definite 4th-order tensors.
            dki=pinv([A_1;A_2])*[D4_1;D4_2];

            %Optional Validation if user adds noise to the data
            if stdv>0
                logS_1=log(squeeze(S_1real(x,y,z,:))/S0(x,y,z));
                logS_2=log(squeeze(S_2real(x,y,z,:))/S0(x,y,z));
                err(:,x,y,z)=abs(Gbig*dki-[logS_1; logS_2]);
            end

            %Store the data
            DKI_D(:,x,y,z)=dki([1:6]); %The 6 unique coefficients of the diffusion tensor D
            %If you want you can put the result in the form of a 3x3 matrix
            %D=[dki(6) dki(5)/2 dki(4)/2
            %   dki(5)/2 dki(3) dki(2)/2
            %   dki(4)/2 dki(2)/2 dki(1)];

            DKI_W(:,x,y,z)=dki([7:21]); %The 15 unique coefficients of the kurtosis tensor W
            %You can see which coefficient is which you can use the function: printTensor(DKI_W(:,x,y,z),4)
            % or if you want to plot a tensor or a tensor field as spherical functions
            % you have to download the plotTensors.m function developed by Angelos Barmpoutis, Ph.D.
            % and then uncomment the following lines.
            %
            % plotTensors(DKI_D(:,1,1,1),1,[321 1]); %3D ellipsoidal plot of D
            % plotTensors(DKI_W(:,1,1,1),1,[321 1]); %3D ellipsoidal plot of W

            %Optional Calculation of Dapp and Kapp
            Dapp(:,x,y,z) = G_1_order2*dki(1:6);
            Kapp(:,x,y,z) = (G_1_order4*dki(7:21))./((G_1_order2*dki(1:6)).^2);

        end
    end
end

%%%%%% OPTIONAL: PRINT RESULTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
all_err=[];
all_Dapp=[];
all_Kapp=[];
counter=0;
meanD=zeros(6,1);
meanW=zeros(15,1);
for x=1:size(S,1)
    for y=1:size(S,2)
        for z=1:size(S,3)
            if stdv>0
               all_err=[all_err; err(:,x,y,z)];
            end
            all_Dapp=[all_Dapp; Dapp(:,x,y,z)];
            all_Kapp=[all_Kapp; Kapp(:,x,y,z)];
            meanD=meanD+DKI_D(:,x,y,z);
            meanW=meanW+DKI_W(:,x,y,z);
            counter=counter+1;
        end
    end
end
meanD=meanD/counter;
meanW=meanW/counter;

fprintf(1,'\n-----RESULTS:-----\n\n');
fprintf('Number of fitted voxels: %d\n',counter);
if stdv>0
    fprintf(1,'Fitting Error: %.4f (std. dev. %.4f)\n',mean(all_err),std(all_err));
end
fprintf(1,'Mean Dapp: %.4f (std. dev. %.4f)\n',mean(all_Dapp),std(all_Dapp));
fprintf(1,'Mean Kapp: %.4f (std. dev. %.4f)\n',mean(all_Kapp),std(all_Kapp));
fprintf(1,'\nMean Diffusion Tensor D:\n');
printTensor(meanD,2);
fprintf(1,'\nMean Kurtosis Tensor W:\n');
printTensor(meanW,4);

% If you want to plot a Diffusion or Kurtosis tensor as spherical functions
% you have to download the plotTensors.m function developed by Angelos Barmpoutis, Ph.D.
% and then uncomment the following lines.
%
% subplot(1,2,1)
% plotTensors(meanD,1,[321 1]);
% title('Mean Diffusion Tensor');
% subplot(1,2,2);
% plotTensors(meanW,1,[321 1]);
% title('Mean Kurtosis Tensor');

fprintf(1,'\nIf you use this software please cite the following papers on DKI and DTI estimation:\n');
fprintf(1,'1) A. Barmpoutis et al. "Diffusion Kurtosis Imaging: Robust Estimation from DW-MRI\n');
fprintf(1,'   using Homogeneous Polynomials", In the Proceedings of ISBI, 2011, pp. 262-265.\n');
fprintf(1,'2) A. Barmpoutis and B.C. Vemuri, "A Unified Framework for Estimating Diffusion Tensors\n');
fprintf(1,'   of any order with Symmetric Positive-Definite Constraints",\n');
fprintf(1,'   In the Proceedings of ISBI, 2010, pp. 1385-1388.\n');


text:
%io_CSIload_rda.m
%Brenden Kadota, SunnyBrook Hosptial 2021.
%
% USAGE:
% [out] = io_CSIload_rda(rda_filename);
%
% DESCRIPTION:
% Reads in siemens rda data (.rda file).
%
% io_CSIload_rda outputs the data in CSI structure format, with fields corresponding to time
% scale, fids, frequency scale, spectra, and header fields containing
% information about the acquisition.  The resulting matlab structure can be
% operated on by the other functions in this CSI toolbox.
%
% INPUTS:
% rda_filename   = filename of Siemens rda data to load.
%
% OUTPUTS:
% out        = Input dataset in FID-A structure format.

function [out] = io_CSIload_rda(rda_filename)

    file = fopen(rda_filename);

    head_start_text = '>>> Begin of header <<<';
    head_end_text   = '>>> End of header <<<';

    tline = fgets(file);

    while (~contains(tline , head_end_text))
        tline = fgets(file);
        if (~contains(tline , head_start_text) + ~contains(tline , head_end_text ) == 2)

            % Store this data in the appropriate format

            occurence_of_colon = strfind(tline, ':');
            variable = tline(1:occurence_of_colon-1) ;
            value    = tline(occurence_of_colon+1 : length(tline)) ;

            switch variable
                case { 'PatientID' , 'PatientName' , 'StudyDescription' , 'PatientBirthDate' , 'StudyDate' , 'StudyTime' , 'PatientAge' , 'SeriesDate' , ...
                        'SeriesTime' , 'SeriesDescription' , 'ProtocolName' , 'PatientPosition' , 'ModelName' , 'StationName' , 'InstitutionName' , ...
                        'DeviceSerialNumber', 'InstanceDate' , 'InstanceTime' , 'InstanceComments' , 'SequenceName' , 'SequenceDescription' , 'Nucleus' ,...
                        'TransmitCoil' }
                    rda.(variable) = value;
                case { 'PatientSex' }
                    % Sex converter! (int to M,F,U)
                    switch value
                        case 0
                            rda.sex = 'Unknown';
                        case 1
                            rda.sex = 'Male';
                        case 2

                            rda.sex = 'Female';
                    end

                case {  'SeriesNumber' , 'InstanceNumber' , 'AcquisitionNumber' , 'NumOfPhaseEncodingSteps' , 'NumberOfRows' , 'NumberOfColumns' , 'VectorSize' }
                    %Integers
                    rda.(variable) = str2double(value);
                case { 'PatientWeight' , 'TR' , 'TE' , 'TM' , 'DwellTime' , 'NumberOfAverages' , 'MRFrequency' , 'MagneticFieldStrength' , 'FlipAngle' , ...
                        'SliceThickness' ,  'FoVHeight' , 'FoVWidth', 'FoV3D' , 'PercentOfRectFoV' , 'PixelSpacingRow' , 'PixelSpacingCol', 'PixelSpacing3D'}
                    %Floats
                    rda.(variable) = str2double(value);
                case {'SoftwareVersion[0]' }
                    rda.software_version = value;
                case {'CSIMatrixSize[0]' }
                    rda.CSIMatrix_Size(1) = str2double(value);
                case {'CSIMatrixSize[1]' }
                    rda.CSIMatrix_Size(2) = str2double(value);
                case {'CSIMatrixSize[2]' }
                    rda.CSIMatrix_Size(3) = str2double(value);
                case {'CSIMatrixSizeOfScan[0]'}
                    rda.CSIMatrixSizeOfScan(1) = str2double(value);
                case {'CSIMatrixSizeOfScan[1]'}
                    rda.CSIMatrixSizeOfScan(2) = str2double(value);
                case {'CSIMatrixSizeOfScan[2]'}
                    rda.CSIMatrixSizeOfScan(3) = str2double(value);
                case {'PositionVector[0]' }
                    rda.PositionVector(1) = str2double(value);
                case {'PositionVector[1]' }
                    rda.PositionVector(2) = str2double(value);
                case {'PositionVector[2]' }
                    rda.PositionVector(3) = str2double(value);
                case {'RowVector[0]' }
                    rda.RowVector(1) = str2double(value);
                case {'RowVector[1]' }
                    rda.RowVector(2) = str2double(value);
                case {'RowVector[2]' }
                    rda.RowVector(3) = str2double(value);
                case {'ColumnVector[0]' }
                    rda.ColumnVector(1) = str2double(value);
                case {'ColumnVector[1]' }
                    rda.ColumnVector(2) = str2double(value);
                case {'ColumnVector[2]' }
                    rda.ColumnVector(3) = str2double(value);
                case {'VOIPositionSag' }
                    rda.imageOrigin(1) = str2double(value);
                case {'VOIPositionCor' }
                    rda.imageOrigin(2) = str2double(value);
                case {'VOIPositionTra' }
                    rda.imageOrigin(3) = str2double(value);
                case {'VOINormalSag' }
                    rda.sliceVector(1) = str2double(value);
                case {'VOINormalCor' }
                    rda.sliceVector(2) = str2double(value);
                case {'VOINormalTra' }
                    rda.sliceVector(3) = str2double(value);

                otherwise
                    % We don't know what this variable is.  Report this just to keep things clear
                    %disp(['Unrecognised variable ' , variable ]);
            end

        else
            % Don't bother storing this bit of the output
        end

    end

    %
    % So now we should have got to the point after the header text
    %
    % Siemens documentation suggests that the data should be in a double complex format (8bytes for real, and 8 for imaginary?)
    %

    bytes_per_point = 16;
    complex_data = fread(file , rda.CSIMatrix_Size(1) * rda.CSIMatrix_Size(2) * rda.CSIMatrix_Size(3) *rda.VectorSize * 2 , 'double');
    %fread(fid , 1, 'double');  %This was a check to confirm that we had read all the data (it passed!)
    fclose(file);

    % Now convert this data into something meaningful

    %Reshape so that we can get the real and imaginary separated
    hmm = reshape(complex_data,  2 , rda.VectorSize , rda.CSIMatrix_Size(1) ,  rda.CSIMatrix_Size(2) ,  rda.CSIMatrix_Size(3) );

    %Combine the real and imaginary into the complex matrix
    fids = complex(hmm(1,:,:,:,:),hmm(2,:,:,:,:));
    fids = squeeze(fids);

    %fids = conj(fids);

    % make calculations for the output mrs structure
    sz = size(fids);
    dwelltime = rda.DwellTime/1000000;
    spectralwidth = 1/dwelltime;
    txfrq = rda.MRFrequency*1000000;
    dims.t = 1;
    dims.subSpecs = 0;
    dims.x = 2;
    fids = flip(fids, dims.x);
    dims.y = 3;
    dims.z = 0;
    dims.coils = 0;
    dims.averages = 0;
    dims.extras = 0;
    dims.ky = 0;
    dims.kx = 0;
    Bo = rda.MagneticFieldStrength;
    rawAverages = rda.NumberOfAverages;
    averages = 1;
    subspecs =1;
    rawSubspecs = 'na';
    date = strip(rda.StudyDate);
    seq = rda.SequenceDescription;
    TE = rda.TE;
    TR = rda.TR;
    pointsToLeftShift = 'N/A';

    %Calculate t and ppm arrays using the calculated parameters:
    f = createCoordinates(spectralwidth/2, spectralwidth/sz(1));
    ppm = -f/(Bo*42.577);
    ppm = ppm+4.6082;

    t = 0:dwelltime:(sz(1)-1)*dwelltime;

    out = struct();
    %FILLING IN DATA STRUCTURE
    out = setData(out, fids);
    out = setSize(out, sz);
    out = setPPM(out, ppm);
    out = setAdcTime(out, t);
    out = setSpectralTime(out, t);
    out = setSpectralDwellTime(out, dwelltime);
    out = setSpectralWidth(out, spectralwidth);
    out.txfrq = txfrq;
    out.date = datetime(str2double(date(1:4)), str2double(date(5:6)), str2double(date(7:8)));
    out.dims = dims;
    out.Bo = Bo;
    out.averages = averages;
    out.rawAverages = rawAverages;
    out.subspecs = subspecs;
    out.rawSubspecs = rawSubspecs;
    out.seq = seq;
    out.te = TE;
    out.tr = TR;
    out.pointsToLeftshift = pointsToLeftShift;
    out = setFov(out, 'x', rda.FoVHeight);
    out = setFov(out, 'y', rda.FoVWidth);
    out = setFov(out, 'z', rda.FoV3D);
    out = setVoxelSize(out, 'x', getFov(out, 'x')/getSizeFromDimensions(out, {'x'}));
    out = setVoxelSize(out, 'y', getFov(out, 'y')/getSizeFromDimensions(out, {'y'}));
    if(getDimension(out, 'z') ~= 0)
        out = setVoxelSize(out, 'z', getFov(out, 'z')/getSizeFromDimensions(out, {'z'}));
    else
        out = setVoxelSize(out, 'z', getFov(out, 'z'));
    end

    affine_matrix = calculateAffineMatrix(rda);
    out.affineMatrix = affine_matrix;
    out = setImageOrigin(out, rda.imageOrigin);
    out = findCoordinates(out);
    out = setFlags(out);
end

function affine_matrix = calculateAffineMatrix(rda)
    %may have to update x and y orientations
    affine_matrix_rotation = [rda.RowVector', rda.ColumnVector', rda.sliceVector', [0 0 0]'];
    affine_matrix_rotation = cat(1, affine_matrix_rotation,  [0, 0, 0, 1]);

    affine_matrix_translate = eye(4);
    affine_matrix_translate(1,4) = rda.PositionVector(1);
    affine_matrix_translate(2,4) = rda.PositionVector(2);
    affine_matrix_translate(3,4) = rda.PositionVector(3);

    affine_matrix_scale = eye(4);
    affine_matrix_scale(1,1) = rda.PixelSpacingCol;
    affine_matrix_scale(2,2) = rda.PixelSpacingRow;
    affine_matrix_scale(3,3) = rda.PixelSpacing3D;

    % for some reason rotation seems to be in the wrong direction. Fix here
    affine_reverse_rotation = affine_matrix_rotation';
    middleVoxelY = rda.CSIMatrix_Size(2)/2;
    middleVoxelX = rda.CSIMatrix_Size(1)/2;
    middleCoordinate = affine_matrix_scale * [middleVoxelX; middleVoxelY; 0; 1];
    middleCoordinateRoate = affine_reverse_rotation * middleCoordinate;
    adjacentTriangleEdge = middleCoordinateRoate(2);
    hypotenuseLength = middleCoordinate(2);
    rotationDistanceAxial = acos(adjacentTriangleEdge/hypotenuseLength);
    oppositeTriangleEdge = sin(rotationDistanceAxial) * hypotenuseLength;
    affine_matrix_coorection = [1 0 0 0
                                0 1 0 0
                                0 0 1 -abs(2*oppositeTriangleEdge)
                                0 0 0 1];

    affine_matrix = affine_matrix_translate * affine_matrix_coorection * affine_reverse_rotation * affine_matrix_scale;
end

function out = setFlags(out)
    %FILLING IN THE FLAGS
    out.flags.writtentostruct = 1;
    out.flags.gotparams = 1;
    out.flags.leftshifted = 0;
    out.flags.filtered = 0;
    out.flags.zeropadded = 0;
    out.flags.freqcorrected = 0;
    out.flags.phasecorrected = 0;
    out.flags.averaged = 1;
    out.flags.addedrcvrs = 1;
    out.flags.subtracted = 0;
    out.flags.writtentotext = 0;
    out.flags.downsampled = 0;
    out.flags.spatialFT = 1;
    out.flags.spectralFT = 0;
    out.flags.isCartesian = 1;
    if out.dims.subSpecs == 0
        out.flags.isISIS = 0;
    else
        out.flags.isISIS = (out.sz(out.dims.subSpecs)==4);
    end
end

function out = findCoordinates(out)
    xCoordinate = createCoordinates(getFov(out, 'x')/2, getVoxSize(out, 'x'));
    yCoordinate = createCoordinates(getFov(out, 'y')/2, getVoxSize(out, 'y'));
    zCoordinate = createCoordinates(getFov(out, 'z')/2, getVoxSize(out, 'z'));
    out = setCoordinates(out, 'x', xCoordinate);
    out = setCoordinates(out, 'y', yCoordinate);
    out = setCoordinates(out, 'z', zCoordinate);
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [final_sound] = im2sound(filename, ext, f_sample, f_low, ...
    f_high, amp_mod, sample_t)

%INPUTS:
%'filename' - Name of the image to be encoded (not including extension
%ext' - Extension of the image (not including "." at the beginning).
%'f_sample' - Sampling frequency (Hz)
%'f_low' - Lowest frequency (Hz) (e.g. 40)
%'f_high' - Highest frequency (Hz) (e.g. 6000)
%'amp_mod' - Multiplication factor for the amplitude.  Decrease until
%image is clear. Too high and the waveform clips.  Too low and the image
%is very dark (e.g. 0.00002)
%'sample_t' - Duration of the sample in seconds.  Longer samples have
%better quality (e.g. 10)

%OUTPUTS:
%'final_sound' - the final sound containing the image.  This is
%automatically saved to a .wav file with the original image filename

%INITIALISING VARIABLES:
%The waveform at each time point.  This is reset at the beginning of each
%time point
temp_sound = 0;
%The final waveform
final_sound = 0;

%MAIN BODY
%Loading the sample image and calculating the image size
raw_im = imread(strcat(filename,'.',ext));
size_raw_im = size(raw_im);

%Making a frequency table for the height of the image.  Each row of the
%image is assigned a particular frequency from the corresponding row of
%this table.  The frequencies are linearly distributed between the highest and
%lowest user-definied frequencies.  "f_step" is the increment between each
%adjacent frequency
f_step = (f_high - f_low)/size_raw_im(1,1);
f_table = (f_high:-f_step:f_low);

%The final sound will dwell on each column of the image for a specific
%time.  This time is defined by "t_start" and "t_end".  It depends on how
%long the user determined the sound-clip should be and how wide (how many
%columns) the image is.
t_step = (sample_t/size_raw_im(1,2));

%Initial values for the start and end times.  These will be increased at
%the end of each loop iteration (when the script moves onto the next column
%of the image).
t_start = 0;
t_end = t_step;

%The loop which generates the sound file.  At each iteration it generates a
%segment of the final sound file, which is temporarily saved to
%"temp_sound".  This segment is built up of frequencies from that
%particular column of the image.
for j = 1:size_raw_im(1,2)
    %Initialising the variable (the sound for each frequency (row) is added
    %to the existing sound)
    temp_sound = 0;

    %Setting the time in matrix format
    t = t_start:1/f_sample:(t_end);

    %For each iteration of this loop, the script goes down the current
    %column of the image and generates a waveform of the frequency
    %specified in "f_table".  The amplitude of the waveform is determined
    %by the pixel intensity.  This generated waveform is added to all the
    %previously generated waveforms in that particular column
    for i = size_raw_im(1,1):-1:1
        temp_sound = temp_sound+ sin(2*pi*t*f_table(i))*...
            double(raw_im(i,j))*amp_mod;

    end

    %At the end of each column the segment of sound generated is added to
    %the end of the existing sound file ("final_sound").
    final_sound = cat(2,final_sound,temp_sound);

    %The temporary sound is cleared ready for the start of the next column
    clear temp_sound

    %Moving to the next time frame
    t_start = t_start + t_step;
    t_end = t_end + t_step;

end

%This saves "final_sound" to the '.wav' file of the same name as the input
%file
wavwrite(final_sound, f_sample, strcat(filename, '.wav'));


text:
function fem2d_pack_test09 ( )

%*****************************************************************************80
%
%% TEST09 demonstrates ELEMENTS_EPS with T4 elements.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  element_order = 4;
  nelemx = 7;
  nelemy = 5;

  element_num = grid_t4_element_num ( nelemx, nelemy );
  node_num = grid_t4_node_num ( nelemx, nelemy );

  code = 'T4';
  file_name = 'fem2d_pack_test_t4.eps';
  title = 'Grid of T4 Elements';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST09\n' );
  fprintf ( 1, '  ELEMENTS_EPS creates an Encapsulated PostScript\n' );
  fprintf ( 1, '  file containing an image of a T4 mesh.\n' );

  element_node = grid_t4_element ( nelemx, nelemy );

  node = 0;

  for row = 0 : nelemy

    y = ( ( 3*nelemy - row ) *  0.0   ...
        + (          + row ) *  6.0 ) ...
        / ( 3*nelemy       );

    for col = 0 : nelemx

      node = node + 1;

      x = ( ( 2*nelemx - col ) * 0.0   ...
          + (          + col ) * 6.0 ) ...
          / ( 2*nelemx       );

      node_x(node) = x;
      node_y(node) = y;

    end
%
%  Skip over the two rows of interior nodes.
%
    node = node + nelemx;
    node = node + nelemx;

  end
%
%  The coordinates of interior nodes are the average of the vertices.
%
  for element = 1 : element_num
    node = element_node(4,element);
    node_x(node) = sum ( node_x(element_node(1:3,element)) ) / 3.0;
    node_y(node) = sum ( node_y(element_node(1:3,element)) ) / 3.0;
  end

  for element = 1 : element_num
    element_mask(element) = 1;
  end

  elements_eps ( file_name, node_num, node_x, node_y, code, ...
    element_order, element_num, element_mask, element_node, title );

  return
end


text:
function [FaceArea, VertArea] = tess_area(Vertices, Faces)
% TESS_AREA: Compute the surface area associated with each face and each vertex.

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

% Compute the area of all the faces
r12 = Vertices(Faces(:,1),:);        % temporary holding
r13 = Vertices(Faces(:,3),:) - r12;  % negative of r31
r12 = Vertices(Faces(:,2),:) - r12;  % from 1 to 2
FaceArea = sqrt(sum(bst_cross(r12,r13,2).^2, 2)) / 2;

% Compute the triangle area only if needed
if (nargout >= 2)
    % Build vertex-face connectivity matrix, with the area information
    nFaces = size(Faces,1);
    rowno = double([Faces(:,1); Faces(:,2); Faces(:,3)]);
    colno = [1:nFaces, 1:nFaces, 1:nFaces]';
    data  = [FaceArea; FaceArea; FaceArea];
    VertFacesArea = sparse(rowno,colno,data);

    % Compute the vertex area: 1/3 of each triangle involving this vertex
    VertArea = 1/3 * full(sum(VertFacesArea,2));
end

author: brainstorm-tools


text:
function classifier= bbci_apply_evalClassifier(fv, bbci_classifier)
%BBCI_APPLY_EVALCLASSIFIER - Apply classifier to feature vector
%
%Synopsis:
%  CLASSIFIER= bbci_apply_evalClassifier(FV, BBCI_CLASSIFIER)

% 02-2011 Benjamin Blankertz

classifier.x= bbci_classifier.fcn(bbci_classifier.C, fv);
classifier.x= classifier.x(:);

author: bbci


text:
% Profiler extension for Kernel Normalized Least-Mean-Square algorithm with
% Coherence Criterion
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

classdef knlms_profiler < knlms

    properties (GetAccess = 'public', SetAccess = 'private')
        elapsed = 0; % elapsed time
        prev_dict_size = 0; % previous dictionary size for growth check
    end

    methods

        function kaf = knlms_profiler(parameters) % constructor
            if nargin<1, parameters = struct(); end
            kaf = kaf@knlms(parameters);
        end

        function flops = lastflops(kaf) % flops for last iteration
            m = size(kaf.dict,1);
            if kaf.prev_dict_size < m % growing
                m1 = m - 1;
                n2 = 1;
            else
                m1 = m;
                n2 = 0;
            end
            m3 = m;
            floptions = struct(...
                'sum', 3*m3, ...
                'mult', 3*m3 + 1, ...
                sprintf('%s_kernel',kaf.kerneltype), [m1+n2, 1, size(kaf.dict,2)]);

            flops = kflops(floptions);
        end

        %% flops breakdown

        % k = kernel(x,kaf.dict,kaf.kerneltype,kaf.kernelpar);
        % kernel: m1

        % h = kernel(x,kaf.dict,kaf.kerneltype,kaf.kernelpar);
        % kernel: n2 % 1 element when growing, 0 when not

        % kaf.alpha = kaf.alpha + kaf.eta / (kaf.eps + h*h') * (y - h*kaf.alpha) * h';
        % sum: 3*m3
        % mult: 3*m3 + 1
        % div: 1

        %%

        function train_profiled(kaf,x,y)
            kaf.prev_dict_size = size(kaf.dict,1);
            t1 = tic;
            kaf.train(x,y);
            t2 = toc(t1);
            kaf.elapsed = kaf.elapsed + t2;
        end

        function bytes = lastbytes(kaf) % bytes used in last iteration
            m = size(kaf.dict,1);
            bytes = 8*(m + m*size(kaf.dict,2)); % 8 bytes for double precision
            % alpha, dict
        end

    end
end

repo: kafbox


text:
function options = defaultOptions;

% DEFAULTOPTIONS The default options for optimisation.
% FORMAT
% DESC returns a default options vector for optimisation.
% RETURN options : the default options vector.
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006
%
% SEEALSO : scg, conjgrad, quasinew

% NDLUTIL

options = [0,  1e-4, 1e-4, 1e-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1e-8, 0.1, 0];

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni

%Fourier Transfrom properties

%	convolution

%x1(t)=u(t)-u(t-2)   and x2(t)=u(t)-u(t-4)

% F^-1{X1(w)X2(w)}
syms t w
x1=heaviside(t)-heaviside(t-2);
x2=heaviside(t)-heaviside(t-4);
X1=fourier(x1,w);
X2=fourier(x2,w);
right =ifourier(X1*X2,t);
ezplot(right,[0 8]);

% convolution of x1(t) with x2(t)
figure
t1=0:.01:2;
t2=2.01:.01:4;
x1=[ones(size(t1)) zeros(size(t2))];
x2=ones(size([t1 t2]));
y=conv(x1,x2)*.01;
plot(0:.01:8,y);


text:
function [E,Ev,Fe,Ve]=patch2edgeIND(F,V)

%% DERIVE EDGE MATRIX

% E=[]; %Can be improved through memory allocation and indexing in loop
% for i=1:1:size(F,2)-1;
%   E=[E; F(:,i) F(:,i+1)];
% end
% E=[E; F(:,end) F(:,1)];

%Format of column index in F
EColumnInd=[(1:size(F,2)); (1:size(F,2))];
EColumnInd=[EColumnInd(2:end) EColumnInd(1)];

%Derive edges matrix
E=F(:,EColumnInd)'; %Use index into F to create edges matrix
E=reshape(E,2,numel(E)/2)';

%%
E=sort(E,2); %Sort edge order
[E,~,ind2] = unique(E,'rows'); %Removing double edges, i.e. [1  4] = [4  1]

%%

Fe=reshape(1:numel(F),size(F,2),size(F,1))';
Fe=ind2(Fe);

ind_E=(1:size(E,1))'*ones(1,2);

Ve=sort((sparse(E(:),ind_E(:),ind_E(:),size(V,1),size(E,1))),2);
[~,J]=find(Ve);
Ve=full(Ve(:,min(J):end));

A=E(Ve(Ve>0),:);
B=zeros(size(Ve));
B(Ve>0)=A(:,1);
C=zeros(size(Ve));
C(Ve>0)=A(:,2);
D=(1:size(Ve,1))'*ones(1,size(Ve,2));
B(B==D)=0;
C(C==D)=0;
Ev=B+C;

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function chantype = ft_chantype(input, desired)

% FT_CHANTYPE determines for each individual channel what chantype of data it
% represents, e.g. a planar gradiometer, axial gradiometer, magnetometer,
% trigger channel, etc. If you want to know what the acquisition system is
% (e.g. ctf151 or neuromag306), you should not use this function but
% FT_SENSTYPE instead.
%
% Use as
%   type = ft_chantype(hdr)
%   type = ft_chantype(sens)
%   type = ft_chantype(label)
% or as
%   type = ft_chantype(hdr,   desired)
%   type = ft_chantype(sens,  desired)
%   type = ft_chantype(label, desired)
%
% If the desired unit is not specified as second input argument, this
% function returns a Nchan*1 cell-array with a string describing the type
% of each channel.
%
% If the desired unit is specified as second input argument, this function
% returns a Nchan*1 boolean vector with "true" for the channels of the
% desired type and "false" for the ones that do not match.
%
% The specification of the channel types depends on the acquisition system,
% for example the ctf275 system includes the following type of channels:
% meggrad, refmag, refgrad, adc, trigger, eeg, headloc, headloc_gof.
%
% See also FT_READ_HEADER, FT_SENSTYPE, FT_CHANUNIT

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

% these are for remembering the type on subsequent calls with the same input arguments
persistent previous_argin previous_argout

% this is to avoid a recursion loop
persistent recursion
if isempty(recursion)
  recursion = false;
end

if nargin<2
  desired = [];
end

% determine the type of input, this is handled similarly as in FT_CHANUNIT
isheader = isa(input, 'struct')  && isfield(input, 'label') && isfield(input, 'Fs');
isdata   = isa(input, 'struct')  && ~isheader && (isfield(input, 'hdr') || isfield(input, 'grad') || isfield(input, 'elec') || isfield(input, 'opto'));
isgrad   = isa(input, 'struct')  && isfield(input, 'label') && isfield(input, 'pnt')  &&  isfield(input, 'ori'); % old style
iselec   = isa(input, 'struct')  && isfield(input, 'label') && isfield(input, 'pnt')  && ~isfield(input, 'ori'); % old style
isgrad   = (isa(input, 'struct') && isfield(input, 'label') && isfield(input, 'coilpos')) || isgrad;             % new style
iselec   = (isa(input, 'struct') && isfield(input, 'label') && isfield(input, 'elecpos')) || iselec;             % new style
isopto   = isa(input, 'struct')  && isfield(input, 'label') && isfield(input, 'transceiver');
islabel  = isa(input, 'cell')    && ~isempty(input) && isa(input{1}, 'char');

if isheader
  % this speeds up the caching in real-time applications
  input.nSamples = 0;
end

current_argin = {input, desired};
if isequal(current_argin, previous_argin)
  % don't do the chantype detection again, but return the previous output from cache
  chantype = previous_argout{1};
  return
end

if isdata
  % the hdr, grad, elec or opto structure might have a different set of channels
  origlabel = input.label;

  if isfield(input, 'hdr')
    input = input.hdr;
    isheader = true;
  elseif isfield(input, 'grad')
    input = input.grad;
    isgrad = true;
  elseif isfield(input, 'elec')
    input = input.elec;
    iselec = true;
  elseif isfield(input, 'opto')
    input = input.opto;
    isopto = true;
  else
    % at least it contains channel labels
    islabel = true;
  end
end

if isheader
  label = input.label;
  numchan = length(label);
elseif isgrad
  label   = input.label;
  numchan = length(label);
elseif iselec
  label   = input.label;
  numchan = length(label);
elseif islabel
  label   = input;
  numchan = length(label);
elseif isfield(input, 'label')
  % this is a last resort: I don't know what it is, but perhaps the labels are informative
  label   = input.label;
  numchan = length(label);
else
  ft_error('the input that was provided to this function cannot be deciphered');
end

if isfield(input, 'chantype')
  % start with the provided channel types
  chantype = input.chantype(:);
else
  % start with unknown chantype for all channels
  chantype = repmat({'unknown'}, numchan, 1);
end

if ~any(strcmp(chantype, 'unknown'))
  % all channels are known, don't bother doing any further heuristics

elseif ft_senstype(input, 'unknown')
  % don't bother doing subsequent checks to determine the chantype

elseif isheader && (ft_senstype(input, 'neuromag') || ft_senstype(input, 'babysquid74'))
  % channames-KI is the channel kind, 1=meg, 202=eog, 2=eeg, 3=trigger (I am not sure, but have inferred this from a single test file)
  % chaninfo-TY is the Coil chantype (0=magnetometer, 1=planar gradiometer)
  if isfield(input, 'orig') && isfield(input.orig, 'channames')
    for sel=find(input.orig.channames.KI(:)==202)'
      chantype{sel} = 'eog';
    end
    for sel=find(input.orig.channames.KI(:)==2)'
      chantype{sel} = 'eeg';
    end
    for sel=find(input.orig.channames.KI(:)==3)'
      chantype{sel} = 'digital trigger';
    end
    % determine the MEG channel subtype
    selmeg=find(input.orig.channames.KI(:)==1)';
    for i=1:length(selmeg)
      if input.orig.chaninfo.TY(i)==0
        chantype{selmeg(i)} = 'megmag';
      elseif input.orig.chaninfo.TY(i)==1
        % FIXME this might also be a axial gradiometer in case the BabySQUID data is read with the old reading routines
        chantype{selmeg(i)} = 'megplanar';
      end
    end

  elseif isfield(input, 'orig') && isfield(input.orig, 'chs') && isfield(input.orig.chs, 'coil_type')
    % all the chs.kinds and chs.coil_types are obtained from the MNE manual, p.210-211
    for sel=find([input.orig.chs.kind]==1 & [input.orig.chs.coil_type]==2)' % planar gradiometers
      chantype(sel) = {'megplanar'}; % Neuromag-122 planar gradiometer
    end
    for sel=find([input.orig.chs.kind]==1 & [input.orig.chs.coil_type]==3012)' %planar gradiometers
      chantype(sel) = {'megplanar'}; % Type T1 planar grad
    end
    for sel=find([input.orig.chs.kind]==1 & [input.orig.chs.coil_type]==3013)' %planar gradiometers
      chantype(sel) = {'megplanar'}; % Type T2 planar grad
    end
    for sel=find([input.orig.chs.kind]==1 & [input.orig.chs.coil_type]==3014)' %planar gradiometers
      chantype(sel) = {'megplanar'}; % Type T3 planar grad
    end
    for sel=find([input.orig.chs.kind]==1 & [input.orig.chs.coil_type]==3022)' %magnetometers
      chantype(sel) = {'megmag'};    % Type T1 magenetometer
    end
    for sel=find([input.orig.chs.kind]==1 & [input.orig.chs.coil_type]==3023)' %magnetometers
      chantype(sel) = {'megmag'};    % Type T2 magenetometer
    end
    for sel=find([input.orig.chs.kind]==1 & [input.orig.chs.coil_type]==3024)' % magnetometers
      chantype(sel) = {'megmag'};    % Type T3 magenetometer
    end
    for sel=find([input.orig.chs.kind]==1 & [input.orig.chs.coil_type]==7001)' % axial gradiometer
      chantype(sel) = {'megaxial'};
    end
    for sel=find([input.orig.chs.kind]==301)' % MEG reference channel, located far from head
      chantype(sel) = {'ref'};
    end
    for sel=find([input.orig.chs.kind]==2)'   % EEG channels
      chantype(sel) = {'eeg'};
    end
    for sel=find([input.orig.chs.kind]==201)' % MCG channels
      chantype(sel) = {'mcg'};
    end
    for sel=find([input.orig.chs.kind]==3)' % Stim channels
      if any([input.orig.chs(sel).logno] == 101) % new systems: 101 (and 102, if enabled) are digital; low numbers are 'pseudo-analog' (if enabled)
        chantype(sel([input.orig.chs(sel).logno] == 101)) = {'digital trigger'};
        chantype(sel([input.orig.chs(sel).logno] == 102)) = {'digital trigger'};
        chantype(sel([input.orig.chs(sel).logno] <= 32))  = {'analog trigger'};
        others = [input.orig.chs(sel).logno] > 32 & [input.orig.chs(sel).logno] ~= 101 & ...
          [input.orig.chs(sel).logno] ~= 102;
        chantype(sel(others)) = {'other trigger'};
      elseif any(ismember([input.orig.chs(sel).logno], [14 15 16])) % older systems: STI 014/015/016 are digital; lower numbers 'pseudo-analog'(if enabled)
        chantype(sel([input.orig.chs(sel).logno] == 14)) = {'digital trigger'};
        chantype(sel([input.orig.chs(sel).logno] == 15)) = {'digital trigger'};
        chantype(sel([input.orig.chs(sel).logno] == 16)) = {'digital trigger'};
        chantype(sel([input.orig.chs(sel).logno] <= 13)) = {'analog trigger'};
        others = [input.orig.chs(sel).logno] > 16;
        chantype(sel(others)) = {'other trigger'};
      else
        ft_warning('There does not seem to be a suitable trigger channel.');
        chantype(sel) = {'other trigger'};
      end
    end
    for sel=find([input.orig.chs.kind]==202)' % EOG
      chantype(sel) = {'eog'};
    end
    for sel=find([input.orig.chs.kind]==302)' % EMG
      chantype(sel) = {'emg'};
    end
    for sel=find([input.orig.chs.kind]==402)' % ECG
      chantype(sel) = {'ecg'};
    end
    for sel=find([input.orig.chs.kind]==502)' % MISC
      chantype(sel) = {'misc'};
    end
    for sel=find([input.orig.chs.kind]==602)' % Resp
      chantype(sel) = {'respiration'};
    end
  end

elseif ft_senstype(input, 'babysquid74')
  % the name can be something like "MEG 001" or "MEG001" or "MEG 0113" or "MEG0113"
  % i.e. with two or three digits and with or without a space
  sel = myregexp('^MEG', label);
  chantype(sel) = {'megaxial'};

elseif ft_senstype(input, 'neuromag122')
  % the name can be something like "MEG 001" or "MEG001" or "MEG 0113" or "MEG0113"
  % i.e. with two or three digits and with or without a space
  sel = myregexp('^MEG', label);
  chantype(sel) = {'megplanar'};

elseif ft_senstype(input, 'neuromag306') && isgrad
  % there should be 204 planar gradiometers and 102 axial magnetometers
  if isfield(input, 'tra')
    tmp = sum(abs(input.tra)>0,2);
    sel = (tmp==1);
    chantype(sel) = {'megmag'};
    sel = (tmp==2);
    chantype(sel) = {'megplanar'};
  end

elseif ft_senstype(input, 'neuromag306') && islabel
  sel = myregexp('^MEG.*1$', label);
  chantype(sel) = {'megmag'};
  sel = myregexp('^MEG.*2$', label);
  chantype(sel) = {'megplanar'};
  sel = myregexp('^MEG.*3$', label);
  chantype(sel) = {'megplanar'};

elseif ft_senstype(input, 'neuromag306_combined') && islabel
  % the magnetometers are detected, the combined channels remain unknown
  sel = myregexp('^MEG.*1$', label);
  chantype(sel) = {'megmag'};

elseif ft_senstype(input, 'ctf') && isheader
  % The following is according to "CTF MEG(TM) File Formats" pdf, Release 5.2.1
  %
  % eMEGReference      0 Reference magnetometer channel
  % eMEGReference1     1 Reference 1st-order gradiometer channel
  % eMEGReference2     2 Reference 2nd-order gradiometer channel
  % eMEGReference3     3 Reference 3rd-order gradiometer channel
  % eMEGSensor         4 Sensor magnetometer channel located in head shell
  % eMEGSensor1        5 Sensor 1st-order gradiometer channel located in head shell
  % eMEGSensor2        6 Sensor 2nd-order gradiometer channel located in head shell
  % eMEGSensor3        7 Sensor 3rd-order gradiometer channel located in head shell
  % eEEGRef            8 EEG unipolar sensors not on the scalp
  % eEEGSensor         9 EEG unipolar sensors on the scalp
  % eADCRef           10 (see eADCAmpRef below)
  % eADCAmpRef        10 ADC amp channels from HLU or PIU (old electronics)
  % eStimRef          11 Stimulus channel for MEG41
  % eTimeRef          12 Time reference coming from video channel
  % ePositionRef      13 Not used
  % eDACRef           14 DAC channel from ECC or HLU
  % eSAMSensor        15 SAM channel derived through data analysis
  % eVirtualSensor    16 Virtual channel derived by combining two or more physical channels
  % eSystemTimeRef    17 System time showing elapsed time since trial started
  % eADCVoltRef       18 ADC volt channels from ECC
  % eStimAnalog       19 Analog trigger channels
  % eStimDigital      20 Digital trigger channels
  % eEEGBipolar       21 EEG bipolar sensor not on the scalp
  % eEEGAflg          22 EEG ADC over range flags
  % eMEGReset         23 MEG resets (counts sensor jumps for crosstalk purposes)
  % eDipSrc           24 Dipole source
  % eSAMSensorNorm    25 Normalized SAM channel derived through data analy- sis
  % eAngleRef         26 Orientation of head localization field
  % eExtractionRef    27 Extracted signal from each sensor of field generated by each localization coil
  % eFitErr           28 Fit error from each head localization coil
  % eOtherRef         29 Any other type of sensor not mentioned but still valid
  % eInvalidType      30 An invalid sensor

  % start with an empty one
  origSensType = [];
  if isfield(input, 'orig')
    if isfield(input.orig, 'sensType') && isfield(input.orig, 'Chan')
      % the header was read using the open-source MATLAB code that originates from CTF and that was modified by the FCDC
      origSensType = input.orig.sensType;
    elseif isfield(input.orig, 'res4') && isfield(input.orig.res4, 'senres')
      % the header was read using the CTF p-files, i.e. readCTFds
      origSensType =  [input.orig.res4.senres.sensorTypeIndex];
    elseif isfield(input.orig, 'sensor') && isfield(input.orig.sensor, 'info')
      % the header was read using the CTF importer from the NIH and Daren Weber
      origSensType = [input.orig.sensor.info.index];
    end
  end

  if isempty(origSensType)
    ft_warning('could not determine channel chantype from the CTF header');
  end

  for sel=find(origSensType(:)==0)'
    chantype{sel} = 'refmag';
  end
  for sel=find(origSensType(:)==1)'
    chantype{sel} = 'refgrad';
  end
  for sel=find(origSensType(:)==5)'
    chantype{sel} = 'meggrad';
  end
  for sel=find(origSensType(:)==9)'
    chantype{sel} = 'eeg';
  end
  for sel=find(origSensType(:)==11)'
    % Stimulus channel for MEG41
    chantype{sel} = 'trigger';
  end
  for sel=find(origSensType(:)==13)'
    chantype{sel} = 'headloc'; % these represent the x, y, z position of the head coils
  end
  for sel=find(origSensType(:)==17)'
    chantype{sel} = 'clock';
  end
  for sel=find(origSensType(:)==18)'
    chantype{sel} = 'adc';
  end
  for sel=find(origSensType(:)==20)'
    % Digital trigger channels
    chantype{sel} = 'trigger';
  end
  for sel=find(origSensType(:)==28)'
    chantype{sel} = 'headloc_gof'; % these represent the goodness of fit for the head coils
  end
  for sel=find(origSensType(:)==29)'
    chantype{sel} = 'reserved'; % these are "reserved for future use", but relate to head localization
  end

elseif ft_senstype(input, 'ctf') && isgrad
  % in principle it is possible to look at the number of coils, but here the channels are identified based on their name
  sel = myregexp('^M[ZLR][A-Z][0-9][0-9]$', input.label);
  chantype(sel) = {'meggrad'};            % normal gradiometer channels
  sel = myregexp('^S[LR][0-9][0-9]$', input.label);
  chantype(sel) = {'meggrad'};            % normal gradiometer channels in the 64 channel CTF system
  sel = myregexp('^B[GPQR][0-9]$', input.label);
  chantype(sel) = {'refmag'};             % reference magnetometers
  sel = myregexp('^[GPQR][0-9][0-9]$', input.label);
  chantype(sel) = {'refgrad'};            % reference gradiometers

elseif ft_senstype(input, 'ctf') && islabel
  % the channels have to be identified based on their name alone
  sel = myregexp('^M[ZLR][A-Z][0-9][0-9]$', label);
  chantype(sel) = {'meggrad'};            % normal gradiometer channels
  sel = myregexp('^S[LR][0-9][0-9]$', label);
  chantype(sel) = {'meggrad'};            % normal gradiometer channels in the 64 channel CTF system
  sel = myregexp('^B[GPR][0-9]$', label);
  chantype(sel) = {'refmag'};             % reference magnetometers
  sel = myregexp('^[GPQR][0-9][0-9]$', label);
  chantype(sel) = {'refgrad'};            % reference gradiometers

elseif ft_senstype(input, 'bti')
  if isfield(input, 'orig') && isfield(input.orig, 'config')
    configname = {input.orig.config.channel_data.name};
    configtype = [input.orig.config.channel_data.type];

    if ~isequal(configname(:), input.label(:))
      % reorder the channels according to the order in input.label
      [sel1, sel2] = match_str(input.label, configname);
      configname = configname(sel2);
      configtype = configtype(sel2);
      configdata = input.orig.config.channel_data(sel2);
    end
    numloops = zeros(size(configdata));
    for i=1:length(configdata)
      if isfield(configdata(i).device_data, 'total_loops')
        numloops(i) = configdata(i).device_data.total_loops;
      end
    end

    % these are taken from bti2grad
    chantype(configtype==1 & numloops==1) = {'megmag'};
    chantype(configtype==1 & numloops==2) = {'meggrad'};
    chantype(configtype==2) = {'eeg'};
    chantype(configtype==3) = {'ref'}; % not known if mag or grad
    chantype(configtype==4) = {'aux'};
    chantype(configtype==5) = {'trigger'};

    % refine the distinction between refmag and refgrad to make the types
    % in grad and header consistent
    sel = myregexp('^M[CLR][xyz][aA]*$', label);
    chantype(sel) = {'refmag'};
    sel = myregexp('^G[xyz][xyz]A$', label);
    chantype(sel) = {'refgrad'};
  else
    % determine the chantype on the basis of the channel labels
    % all 4D-BTi MEG channels start with "A" followed by a number
    % all 4D-BTi reference channels start with M or G
    % all 4D-BTi EEG channels start with E, except for the 248-MEG/32-EEG system in Warsaw where they end with -1
    sel = myregexp('^A[0-9]+$', label);
    chantype(sel) = {'meg'};
    sel = myregexp('^M[CLR][xyz][aA]*$', label);
    chantype(sel) = {'refmag'};
    sel = myregexp('^G[xyz][xyz]A$', label);
    chantype(sel) = {'refgrad'};

    if isgrad && isfield(input, 'tra')
      gradtype = repmat({'unknown'}, size(input.label));
      gradtype(strncmp('A', input.label, 1)) = {'meg'};
      gradtype(strncmp('M', input.label, 1)) = {'refmag'};
      gradtype(strncmp('G', input.label, 1)) = {'refgrad'};
      % look at the number of coils of the meg channels
      selchan = find(strcmp('meg', gradtype));
      for k = 1:length(selchan)
        ncoils = length(find(input.tra(selchan(k),:)==1));
        if ncoils==1
          gradtype{selchan(k)} = 'megmag';
        elseif ncoils==2
          gradtype{selchan(k)} = 'meggrad';
        end
      end
      [selchan, selgrad] = match_str(label, input.label);
      chantype(selchan) = gradtype(selgrad);
    end

    % deal with additional channel types based on the names
    if isheader && issubfield(input, 'orig.channel_data.chan_label')
      tmplabel = {input.orig.channel_data.chan_label};
      tmplabel = tmplabel(:);
    else
      tmplabel = label; % might work
    end
    sel = find(strcmp('unknown', chantype));
    if ~isempty(sel)
      chantype(sel) = ft_chantype(tmplabel(sel));
      sel       = find(strcmp('unknown', chantype));
      if ~isempty(sel)
        % channels that start with E are assumed to be EEG
        % channels that end with -1 are also assumed to be EEG, see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2389
        chantype(sel(cellfun(@(x) strcmp(x(end-1:end),'-1') || strcmp(x(1),'E'), label(sel)))) = {'eeg'};
      end
    end
  end

elseif ft_senstype(input, 'itab') && isheader
  origtype = [input.orig.ch.type];
  chantype(origtype==0) = {'unknown'};
  chantype(origtype==1) = {'ele'};
  chantype(origtype==2) = {'mag'}; % might be magnetometer or gradiometer, look at the number of coils
  chantype(origtype==4) = {'ele ref'};
  chantype(origtype==8) = {'mag ref'};
  chantype(origtype==16) = {'aux'};
  chantype(origtype==32) = {'param'};
  chantype(origtype==64) = {'digit'};
  chantype(origtype==128) = {'flag'};
  % these are the channels that are visible to FieldTrip
  chansel = 1:input.orig.nchan;
  chantype = chantype(chansel);

elseif ft_senstype(input, 'yokogawa') && isheader
  % This is to recognize Yokogawa channel types from the original header
  % This is from the original documentation
  NullChannel                   = 0;
  MagnetoMeter                  = 1;
  AxialGradioMeter              = 2;
  PlannerGradioMeter            = 3;
  RefferenceChannelMark         = hex2dec('0100');
  RefferenceMagnetoMeter        = bitor( RefferenceChannelMark, MagnetoMeter      );
  RefferenceAxialGradioMeter    = bitor( RefferenceChannelMark, AxialGradioMeter  );
  RefferencePlannerGradioMeter  = bitor( RefferenceChannelMark, PlannerGradioMeter);
  TriggerChannel                = -1;
  EegChannel                    = -2;
  EcgChannel                    = -3;
  EtcChannel                    = -4;
  if ft_hastoolbox('yokogawa_meg_reader')
    label = input.label;
    sel = myregexp('[0-9][0-9][0-9]$', label);
    chantype(sel) = {'null'};
    sel = myregexp('^M[0-9][0-9][0-9]$', label);
    chantype(sel) = {'megmag'};
    sel = myregexp('^AG[0-9][0-9][0-9]$', label);
    chantype(sel) = {'meggrad'};
    sel = myregexp('^PG[0-9][0-9][0-9]$', label);
    chantype(sel) = {'megplanar'};
    sel = myregexp('^RM[0-9][0-9][0-9]$', label);
    chantype(sel) = {'refmag'};
    sel = myregexp('^RAG[0-9][0-9][0-9]$', label);
    chantype(sel) = {'refgrad'};
    sel = myregexp('^RPG[0-9][0-9][0-9]$', label);
    chantype(sel) = {'refplanar'};
    sel = myregexp('^TRIG[0-9][0-9][0-9]$', label);
    chantype(sel) = {'trigger'};
    %% Possible labels categorized in "eeg"
    sel_A = myregexp('^A[^G]*[0-9hzZ]$', label);
    sel_P = myregexp('^P[^G]*[0-9hzZ]$', label);
    sel_T = myregexp('^T[^R]*[0-9hzZ]$', label);
    sel_E = myregexp('^E$', label);
    sel_Z = myregexp('^[zZ]$', label);
    sel_M = myregexp('^M[0-9]$', label);
    sel_O = myregexp('^[BCFION]\w*[0-9hzZ]$', label);
    sel_EEG = myregexp('^EEG[0-9][0-9][0-9]$', label);
    sel = logical( sel_A + sel_P + sel_T + sel_E + sel_Z + sel_M + sel_O + sel_EEG );
    clear sel_A sel_P sel_T sel_E sel_Z sel_M sel_O sel_EEG
    chantype(sel) = {'eeg'};
    %% Additional EOG, ECG labels
    sel = myregexp('^EO[0-9]$', label); % EO
    chantype(sel) = {'eog'};
%    sel = myregexp('^ECG[0-9][0-9][0-9]$', label);
    sel_X = myregexp('^X[0-9]$', label); % X
    sel_ECG = myregexp('^ECG[0-9][0-9][0-9]$', label);
    sel = logical( sel_X + sel_ECG );
    clear sel_X sel_ECG
    chantype(sel) = {'ecg'};
    sel = myregexp('^ETC[0-9][0-9][0-9]$', label);
    chantype(sel) = {'etc'};

%   % shorten names
%    ch_info = input.orig.channel_info.channel;
%    type_orig = [ch_info.type];
%    sel = (type_orig == NullChannel);
%    chantype(sel) = {'null'};
%    sel = (type_orig == MagnetoMeter);
%    chantype(sel) = {'megmag'};
%    sel = (type_orig == AxialGradioMeter);
%    chantype(sel) = {'meggrad'};
%    sel = (type_orig == PlannerGradioMeter);
%    chantype(sel) = {'megplanar'};
%    sel = (type_orig == RefferenceMagnetoMeter);
%    chantype(sel) = {'refmag'};
%    sel = (type_orig == RefferenceAxialGradioMeter);
%    chantype(sel) = {'refgrad'};
%    sel = (type_orig == RefferencePlannerGradioMeter);
%    chantype(sel) = {'refplanar'};
%    sel = (type_orig == TriggerChannel);
%    chantype(sel) = {'trigger'};
%    sel = (type_orig == EegChannel);
%    chantype(sel) = {'eeg'};
%    sel = (type_orig == EcgChannel);
%    chantype(sel) = {'ecg'};
%    sel = (type_orig == EtcChannel);
%    chantype(sel) = {'etc'};

  elseif ft_hastoolbox('yokogawa')
    sel = (input.orig.channel_info(:, 2) == NullChannel);
    chantype(sel) = {'null'};
    sel = (input.orig.channel_info(:, 2) == MagnetoMeter);
    chantype(sel) = {'megmag'};
    sel = (input.orig.channel_info(:, 2) == AxialGradioMeter);
    chantype(sel) = {'meggrad'};
    sel = (input.orig.channel_info(:, 2) == PlannerGradioMeter);
    chantype(sel) = {'megplanar'};
    sel = (input.orig.channel_info(:, 2) == RefferenceMagnetoMeter);
    chantype(sel) = {'refmag'};
    sel = (input.orig.channel_info(:, 2) == RefferenceAxialGradioMeter);
    chantype(sel) = {'refgrad'};
    sel = (input.orig.channel_info(:, 2) == RefferencePlannerGradioMeter);
    chantype(sel) = {'refplanar'};
    sel = (input.orig.channel_info(:, 2) == TriggerChannel);
    chantype(sel) = {'trigger'};
    sel = (input.orig.channel_info(:, 2) == EegChannel);
    chantype(sel) = {'eeg'};
    sel = (input.orig.channel_info(:, 2) == EcgChannel);
    chantype(sel) = {'ecg'};
    sel = (input.orig.channel_info(:, 2) == EtcChannel);
    chantype(sel) = {'etc'};
  end

elseif ft_senstype(input, 'yokogawa') && isgrad
  % all channels in the gradiometer definition are meg
  % chantype(1:end) = {'meg'};
  % channels are identified based on their name: only magnetic as isgrad==1
  sel = myregexp('^M[0-9][0-9][0-9]$', input.label);
  chantype(sel) = {'megmag'};
  sel = myregexp('^AG[0-9][0-9][0-9]$', input.label);
  chantype(sel) = {'meggrad'};
  sel = myregexp('^PG[0-9][0-9][0-9]$', input.label);
  chantype(sel) = {'megplanar'};
  sel = myregexp('^RM[0-9][0-9][0-9]$', input.label);
  chantype(sel) = {'refmag'};
  sel = myregexp('^RAG[0-9][0-9][0-9]$', input.label);
  chantype(sel) = {'refgrad'};
  sel = myregexp('^RPG[0-9][0-9][0-9]$', input.label);
  chantype(sel) = {'refplanar'};

elseif ft_senstype(input, 'yokogawa') && islabel
  % the yokogawa channel labels are a mess, so autodetection is not possible
  % chantype(1:end) = {'meg'};
  sel = myregexp('[0-9][0-9][0-9]$', label);
  chantype(sel) = {'null'};
  sel = myregexp('^M[0-9][0-9][0-9]$', label);
  chantype(sel) = {'megmag'};
  sel = myregexp('^AG[0-9][0-9][0-9]$', label);
  chantype(sel) = {'meggrad'};
  sel = myregexp('^PG[0-9][0-9][0-9]$', label);
  chantype(sel) = {'megplanar'};
  sel = myregexp('^RM[0-9][0-9][0-9]$', label);
  chantype(sel) = {'refmag'};
  sel = myregexp('^RAG[0-9][0-9][0-9]$', label);
  chantype(sel) = {'refgrad'};
  sel = myregexp('^RPG[0-9][0-9][0-9]$', label);
  chantype(sel) = {'refplanar'};
  sel = myregexp('^TRIG[0-9][0-9][0-9]$', label);
  chantype(sel) = {'trigger'};
  %% Possible labels categorized in "eeg"
  sel_A = myregexp('^A[^G]*[0-9hzZ]$', label);
  sel_P = myregexp('^P[^G]*[0-9hzZ]$', label);
  sel_T = myregexp('^T[^R]*[0-9hzZ]$', label);
  sel_E = myregexp('^E$', label);
  sel_Z = myregexp('^[zZ]$', label);
  sel_M = myregexp('^M[0-9]$', label);
  sel_O = myregexp('^[BCFION]\w*[0-9hzZ]$', label);
  sel_EEG = myregexp('^EEG[0-9][0-9][0-9]$', label);
  sel = logical( sel_A + sel_P + sel_T + sel_E + sel_Z + sel_M + sel_O + sel_EEG );
  clear sel_A sel_P sel_T sel_E sel_Z sel_M sel_O sel_EEG
  chantype(sel) = {'eeg'};
  %% Additional EOG, ECG labels
  sel = myregexp('^EO[0-9]$', label); % EO
  chantype(sel) = {'eog'};
% sel = myregexp('^ECG[0-9][0-9][0-9]$', label);
  sel_X = myregexp('^X[0-9]$', label); % X
  sel_ECG = myregexp('^ECG[0-9][0-9][0-9]$', label);
  sel = logical( sel_X + sel_ECG );
  clear sel_X sel_ECG
  chantype(sel) = {'ecg'};
  sel = myregexp('^ETC[0-9][0-9][0-9]$', label);
  chantype(sel) = {'etc'};

elseif ft_senstype(input, 'itab') && isheader
  sel = ([input.orig.ch.type]==0);
  chantype(sel) = {'unknown'};
  sel = ([input.orig.ch.type]==1);
  chantype(sel) = {'unknown'};
  sel = ([input.orig.ch.type]==2);
  chantype(sel) = {'megmag'};
  sel = ([input.orig.ch.type]==8);
  chantype(sel) = {'megref'};
  sel = ([input.orig.ch.type]==16);
  chantype(sel) = {'aux'};
  sel = ([input.orig.ch.type]==64);
  chantype(sel) = {'digital'};
  % not all channels are actually processed by FieldTrip, so only return
  % the types fopr the ones that read_header and read_data return
  chantype = chantype(input.orig.chansel);

elseif ft_senstype(input, 'itab') && isgrad
  % the channels have to be identified based on their name alone
  sel = myregexp('^MAG_[0-9][0-9][0-9]$', label);
  chantype(sel) = {'megmag'};
  sel = myregexp('^MAG_[0-9][0-9]$', label); % for the itab28 system
  chantype(sel) = {'megmag'};
  sel = myregexp('^MAG_[0-9]$', label); % for the itab28 system
  chantype(sel) = {'megmag'};
  sel = myregexp('^REF_[0-9][0-9][0-9]$', label);
  chantype(sel) = {'megref'};
  sel = myregexp('^AUX.*$', label);
  chantype(sel) = {'aux'};

elseif ft_senstype(input, 'itab') && islabel
  % the channels have to be identified based on their name alone
  sel = myregexp('^MAG_[0-9][0-9][0-9]$', label);
  chantype(sel) = {'megmag'};
  sel = myregexp('^REF_[0-9][0-9][0-9]$', label);
  chantype(sel) = {'megref'};
  sel = myregexp('^AUX.*$', label);
  chantype(sel) = {'aux'};

elseif ft_senstype(input, 'eeg') && islabel
  % use an external helper function to define the list with EEG channel names
  chantype(match_str(label, ft_senslabel('eeg1005'))) = {'eeg'};          % this includes all channels from the 1010 and 1020 arrangement
  chantype(match_str(label, ft_senslabel(ft_senstype(input)))) = {'eeg'}; % this will work for biosemi, egi and other detected channel arrangements

elseif ft_senstype(input, 'eeg') && iselec
  % all channels in an electrode definition must be eeg channels
  chantype(:) = {'eeg'};

elseif ft_senstype(input, 'eeg') && isheader
  % use an external helper function to define the list with EEG channel names
  chantype(match_str(input.label, ft_senslabel(ft_senstype(input)))) = {'eeg'};

elseif ft_senstype(input, 'plexon') && isheader
  % this is a complete header that was read from a Plexon *.nex file using read_plexon_nex
  for i=1:numchan
    switch input.orig.VarHeader(i).Type
      case 0
        chantype{i} = 'spike';
      case 1
        chantype{i} = 'event';
      case 2
        chantype{i} = 'interval';  % Interval variables?
      case 3
        chantype{i} = 'waveform';
      case 4
        chantype{i} = 'population'; % Population variables ?
      case 5
        chantype{i} = 'analog';
      otherwise
        % keep the default 'unknown' chantype
    end
  end

end % ft_senstype

% if possible, set additional types based on channel labels
label2type = {
  {'ecg', 'ekg'};
  {'emg'};
  {'eog', 'heog', 'veog'};
  {'lfp'};
  {'eeg'};
  {'trigger', 'trig', 'dtrig'};
  };
for i = 1:numel(label2type)
  for j = 1:numel(label2type{i})
    chantype(intersect(strmatch(label2type{i}{j}, lower(label)), find(strcmp(chantype, 'unknown')))) = label2type{i}(1);
  end
end

if isdata
  % the input was replaced by one of hdr, grad, elec, opto
  [sel1, sel2] = match_str(origlabel, input.label);
  origtype = repmat({'unknown'}, size(origlabel));
  origtype(sel1) = chantype(sel2);
  % the hdr, grad, elec or opto structure might have a different set of channels
  chantype = origtype;
end

if all(strcmp(chantype, 'unknown')) && ~recursion
  % try whether only lowercase channel labels makes a difference
  if islabel
    recursion = true;
    chantype = ft_chantype(lower(input));
    recursion = false;
  elseif isfield(input, 'label')
    input.label = lower(input.label);
    recursion = true;
    chantype = ft_chantype(input);
    recursion = false;
  end
end

if all(strcmp(chantype, 'unknown')) && ~recursion
  % try whether only uppercase channel labels makes a difference
  if islabel
    recursion = true;
    chantype = ft_chantype(upper(input));
    recursion = false;
  elseif isfield(input, 'label')
    input.label = upper(input.label);
    recursion = true;
    chantype = ft_chantype(input);
    recursion = false;
  end
end

if nargin>1
  % return a boolean vector
  if isequal(desired, 'meg') || isequal(desired, 'ref')
    % only compare the first three characters, i.e. meggrad or megmag should match
    chantype = strncmp(desired, chantype, 3);
  else
    chantype = strcmp(desired, chantype);
  end
end

% remember the current input and output arguments, so that they can be
% reused on a subsequent call in case the same input argument is given
current_argout = {chantype};
previous_argin  = current_argin;
previous_argout = current_argout;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function match = myregexp(pat, list)
match = false(size(list));
for i=1:numel(list)
  match(i) = ~isempty(regexp(list{i}, pat, 'once'));
end


text:
function Mh = hmxPlusMtimes(Mh,alpha,Ml,Mr)
%+========================================================================+
%|                                                                        |
%|         OPENHMX - LIBRARY FOR H-MATRIX COMPRESSION AND ALGEBRA         |
%|           openHmx is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : hmxPlusMtimes.m                               |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Sum and product of H-Matrix with the rule     |
%|  `---'  |                                                              |
%+========================================================================+

% Check dimensions
if (size(Mh,1) ~= size(Ml,1))
    error('hmxPlusMtimes.m : matrix dimensions must agree.')
end
if (size(Mh,2) ~= size(Mr,2))
    error('hmxPlusMtimes.m : matrix dimensions must agree.')
end
if (size(Ml,2) ~= size(Mr,1))
    error('hmxPlusMtimes.m : matrix dimensions must agree.')
end

%%% H-Matrix + alpha * H-Matrix * H-Matrix --> H-Matrix
if isa(Mh,'hmx') && isa(Ml,'hmx') && isa(Mr,'hmx')

    % H-Matrix + H-Matrix * H-Matrix --> H-Matrix   (recursion)
    if (Mh.typ==0) && (Ml.typ==0) && (Mr.typ==0)
        % First bloc : M1 -> M1 + alpha * (Ml1 * Mr1 + Ml2 * Mr3)
        Mh.chd{1} = hmxPlusMtimes(Mh.chd{1},alpha,Ml.chd{1},Mr.chd{1});
        Mh.chd{1} = hmxPlusMtimes(Mh.chd{1},alpha,Ml.chd{2},Mr.chd{3});

        % Second bloc : M2 -> M2 + alpha * (Ml1 * Mr2 + Ml2 * Mr4)
        Mh.chd{2} = hmxPlusMtimes(Mh.chd{2},alpha,Ml.chd{1},Mr.chd{2});
        Mh.chd{2} = hmxPlusMtimes(Mh.chd{2},alpha,Ml.chd{2},Mr.chd{4});

        % Third bloc : M3 -> M3 + alpha * (Ml3 * Mr1 + Ml4 * Mr3)
        Mh.chd{3} = hmxPlusMtimes(Mh.chd{3},alpha,Ml.chd{3},Mr.chd{1});
        Mh.chd{3} = hmxPlusMtimes(Mh.chd{3},alpha,Ml.chd{4},Mr.chd{3});

        % Fourth bloc : M4 -> M4 + alpha * (Ml3 * Mr2  + Ml4 * Mr4)
        Mh.chd{4} = hmxPlusMtimes(Mh.chd{4},alpha,Ml.chd{3},Mr.chd{2});
        Mh.chd{4} = hmxPlusMtimes(Mh.chd{4},alpha,Ml.chd{4},Mr.chd{4});

        % Fusion
        Mh = hmxFusion(Mh);

    % H-Matrix + H-Matrix * Compr --> H-Matrix
    elseif (Mh.typ==0) && (Ml.typ==0) && (Mr.typ==1)
        A  = alpha * ( Ml * Mr.dat{1} );
        B  = Mr.dat{2};
        Mh = hmxPlusAB(Mh,A,B);

    % H-Matrix + H-Matrix * Full --> unavailable
    elseif (Mh.typ==0) && (Ml.typ==0) && (Mr.typ==2)
        error('hmxPlusMtimes : unvailable case')

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

    % H-Matrix + Compr * H-Matrix --> H-Matrix
    elseif (Mh.typ==0) && (Ml.typ==1) && (Mr.typ==0)
        A  = alpha * Ml.dat{1};
        B  = Ml.dat{2} * Mr;
        Mh = hmxPlusAB(Mh,A,B);

    % H-Matrix + Compr * Compr --> H-Matrix
    elseif (Mh.typ==0) && (Ml.typ==1) && (Mr.typ==1)
        A  = alpha * Ml.dat{1};
        B  = (Ml.dat{2} * Mr.dat{1}) * Mr.dat{2};
        Mh = hmxPlusAB(Mh,A,B);

    % H-Matrix + Compr * Full --> H-Matrix
    elseif (Mh.typ==0) && (Ml.typ==1) && (Mr.typ==2)
        error('hmxPlusMtimes : unvailable case')

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

    % H-Matrix + Full * H-Matrix --> Full
    elseif (Mh.typ==0) && (Ml.typ==2)
        error('hmxPlusMtimes : unvailable case')

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

    % Compr + H-Matrix * H-Matrix --> H-Matrix
    elseif (Mh.typ==1) && (Ml.typ==0) && (Mr.typ==0)
        Mh.typ = 0;
        Mh.row = Ml.row;
        Mh.col = Mr.col;
        for i = 1:4
            Mh.chd{i} = hmx(Ml.chd{i}.pos{1},Mr.chd{i}.pos{2},...
                Mh.dat{1}(Ml.row{i},:) , Mh.dat{2}(:,Mr.col{i}) , ...
                Ml.tol ) ;
        end
        Mh.dat = [];
        Mh     = hmxPlusMtimes(Mh,alpha,Ml,Mr);

    % Compr + H-Matrix * Compr --> Compr
    elseif (Mh.typ==1) && (Ml.typ==0) && (Mr.typ==1)
        A  = alpha * (Ml * Mr.dat{1});
        B  = Mr.dat{2};
        Mh = hmxPlusAB(Mh,A,B);

    % Compr + H-Matrix * Full --> Full
    elseif (Mh.typ==1) && (Ml.typ==0) && (Mr.typ==2)
        error('hmxPlusMtimes : unvailable case')

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

    % Compr + Compr * H-Matrix --> Compr
    elseif (Mh.typ==1) && (Ml.typ==1) && (Mr.typ==0)
        A  = alpha * Ml.dat{1};
        B  = Ml.dat{2} * Mr;
        Mh = hmxPlusAB(Mh,A,B);

    % Compr + Compr * Compr --> Compr
    elseif (Mh.typ==1) && (Ml.typ==1) && (Mr.typ==1)
        A  = alpha * Ml.dat{1};
        B  = (Ml.dat{2} * Mr.dat{1}) * Mr.dat{2};
        Mh = hmxPlusAB(Mh,A,B);

    % Compr + Compr * Full --> Compr
    elseif (Mh.typ==1) && (Ml.typ==1) && (Mr.typ==2)
        A  = alpha * Ml.dat{1};
        B  = Ml.dat{2} * Mr.dat;
        Mh = hmxPlusAB(Mh,A,B);

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

    % Compr + Full * H-Matrix --> Full
    elseif (Mh.typ==1) && (Ml.typ==2) && (Mr.typ==0)
        error('hmxPlusMtimes : unvailable case')

    % Compr + Full * Compr --> Compr
    elseif (Mh.typ==1) && (Ml.typ==2) && (Mr.typ==1)
        A  = alpha * (Ml.dat * Mr.dat{1});
        B  = Mr.dat{2};
        Mh = hmxPlusAB(Mh,A,B);

    % Compr + Full * Full --> Full
    elseif (Mh.typ==1) && (Ml.typ==2) && (Mr.typ==2)
        Mh.dat = full(Mh) + alpha * (Ml.dat * Mr.dat);
        Mh.typ = 2;

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

    % Full + --- --> Full
    elseif (Mh.typ==2)
        Mh.dat = full(Mh) + alpha * (full(Ml) * full(Mr));
        Mh.typ = 2;

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

    % Unknown type
    else
        error('hmxPlusMtimes : unvailable case')
    end

%%% Unavailable
else
    error('hmxPlusMtimes : unvailable case')
end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan