text:
function [feat,idxs_bbox_pair,cb1,cb2] = get_spatial_features_diff_img_dx_dy(boxes,idxs_bbox_pair,rot_offset,scores,class_id,bVis,lab)

if (nargin < 3)
    rot_offset = 0;
end

if (nargin < 6)
    bVis = false;
end

if (nargin < 7)
    lab = ones(size(boxes,1),1);
end

cb1 = [mean(boxes(idxs_bbox_pair(:,1),[1 3]),2) mean(boxes(idxs_bbox_pair(:,1),[2 4]),2)];
cb2 = [mean(boxes(idxs_bbox_pair(:,2),[1 3]),2) mean(boxes(idxs_bbox_pair(:,2),[2 4]),2)];

deltaX = abs(cb1(:,1)-cb2(:,1));
deltaY = abs(cb1(:,2)-cb2(:,2));
angle = atan2(deltaY,deltaX);
angle = angle - rot_offset;
angle = wrapMinusPiPifast(angle);

% feat = cat(2,dist,angle,scores(idxs_bbox_pair(:,1),:),scores(idxs_bbox_pair(:,2),:));
feat = cat(2,deltaX,deltaY,angle,scores(idxs_bbox_pair(:,1),:),scores(idxs_bbox_pair(:,2),:));

if (bVis)
    idxs = 1:min(size(cb1,1),100);
%     figure(100);clf;

    plot(cb1(idxs,1),cb1(idxs,2),'b+','MarkerSize',10);
    plot(cb2(idxs,1),cb2(idxs,2),'g+','MarkerSize',10);

    for i = 1:length(idxs)
        if (lab(idxs(i)) == 1)
            plot([cb1(idxs(i),1); cb2(idxs(i),1)],[cb1(idxs(i),2); cb2(idxs(i),2)],'r-','lineWidth',1);
        else
            plot([cb1(idxs(i),1); cb2(idxs(i),1)],[cb1(idxs(i),2); cb2(idxs(i),2)],'b-','lineWidth',1);
        end
    end

    legendName = {['pidx ' num2str(class_id(1)) '-' num2str(class_id(2))]};
    legend(legendName);
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function test_polygon_recovery(record)
import iris.inflate_region;

if nargin < 1
  record = false;
end

r = 1;

lb = [-r;-r];
ub = [r;r];

n_obs = 7;
th = linspace(0,2*pi,n_obs+1);
% th = [0,pi/4,pi/2,pi,5*pi/4,3*pi/2,2*pi];

% th = cumsum(rand(6, 1) * 3*pi./2);
% th = th * (2*pi / th(end));

obstacles = zeros(2, 2, length(th)-1);
for j = 1:length(th)-1
  obstacles(:,:,j) = [[r*cos(th(j));r*sin(th(j))],[r*cos(th(j+1));r*sin(th(j+1))]];
end

A_bounds = [-1,0;0,-1;1,0;0,1];
b_bounds = [-lb; ub];

start = 0.5 .* (rand(2,1) .* (ub - lb) + lb);
% start = [0;0];
% start = [0;0.98];

% profile on
[A,b,C,d,results] = inflate_region(obstacles, A_bounds, b_bounds, start);
% profile viewer
iris.drawing.animate_results(results, record);

end

author: rdeits

repo: iris-distro

save_path: github-repos/MATLAB/rdeits-iris-distro


text:
function S=inputdlgStruct(varargin)

%% Parse input
switch nargin
    case 1
        nameCell=varargin{1};
        defaultOptions=[];
        dialogTitle=[];
    case 2
        nameCell=varargin{1};
        defaultOptions=varargin{2};
        dialogTitle=[];
    case 3
        nameCell=varargin{1};
        defaultOptions=varargin{2};
        dialogTitle=varargin{3};
end

if isempty(dialogTitle)
   dialogTitle='Input dialog';
end

if isempty(defaultOptions)
   defaultOptions=repmat({''},size(nameCell));
end

%%
s=25+max([cellfun(@numel,nameCell) cellfun(@numel,defaultOptions)]); %Set sizes
entryCell = inputdlg(nameCell,dialogTitle,[1 s],defaultOptions); %Open dialog box
[S]=cellPair2struct(nameCell,entryCell,1); %Convert output to structure

%%
% _*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 [accuracy, frames] = estimate_accuracy(trajectory, sequence, varargin)
% estimate_accuracy Calculate accuracy score
%
% Calculate accuracy score as average overlap over the entire sequence.
%
% Input:
% - trajectory (cell): A trajectory as a cell array of regions.
% - sequence (cell or structure): Can be another trajectory or a valid sequence descriptor.
% - varargin[Burnin] (integer): Number of frames that have to be ignored after the failure.
% - varargin[IgnoreUnknown] (boolean): Ignore frames where the overlap is
% unknown.
% - varargin[BindWithin] (boolean or vector): Bind the overlap calculation to the region
% within the image. If the sequence variable is a sequence structure then a boolean value
% is sufficient to establishe bounding region. Otherwise a bounding region has to be specified
% manually.
%
% Output:
% - accuracy (double): Average overlap.
% - frames (double vector): Per-frame overlaps.
%

ignore_unknown = true;
burnin = 0;
bind_within = get_global_variable('bounded_overlap', true);

for j=1:2:length(varargin)
    switch lower(varargin{j})
        case 'burnin', burnin = max(0, varargin{j+1});
        case 'ignoreunknown', ignore_unknown = varargin{j+1};
        case 'bindwithin', bind_within = varargin{j+1};
        otherwise, error(['unrecognized argument ' varargin{j}]);
    end
end

if burnin > 0

    mask = cellfun(@(r) numel(r) == 1 && r == 1, trajectory, 'UniformOutput', true);

    if is_octave()
        se = logical([zeros(burnin - 1, 1); ones(burnin, 1)]);
    else
        se = strel('arbitrary', [zeros(burnin - 1, 1); ones(burnin, 1)]);
    end;

    % ignore the next 'burnin' frames
    mask = imdilate(mask, se);

else

    mask = false(size(trajectory, 1), 1);

end;

if ~ignore_unknown
    unknown = cellfun(@(r) numel(r) == 1 && r == 0, trajectory, 'UniformOutput', true);
end;

trajectory(mask) = {0};

if islogical(bind_within)
    if bind_within && isstruct(sequence)
        bounds = [sequence.width, sequence.height] - 1;
    else
        bounds = [];
    end;
else
    bounds = bind_within;
end;

if isstruct(sequence)
    frames = calculate_overlap(trajectory, sequence_get_region(sequence, 1:sequence.length), bounds);
else
    frames = calculate_overlap(trajectory, sequence, bounds);
end;

if ~ignore_unknown
    frames(unknown) = 0;
end;

overlap = frames(~isnan(frames)); % filter-out illegal values

% Handle cases, where no overlap is available
if isempty(overlap)
    accuracy = 0;
else
    accuracy = mean(overlap);
end;

author: votchallenge

repo: toolkit-legacy


text:
function [sigma,bnd] = bdsqr(alpha,beta)

% BDSQR: Compute the singular values and bottom element of
%        the left singular vectors of a (k+1) x k lower bidiagonal
%        matrix with diagonal alpha(1:k) and lower bidiagonal beta(1:k),
%        where length(alpha) = length(beta) = k.
%
% [sigma,bnd] = bdsqr(alpha,beta)
%
% Input parameters:
%   alpha(1:k)   : Diagonal elements.
%   beta(1:k)    : Sub-diagonal elements.
% Output parameters:
%   sigma(1:k)  : Computed eigenvalues.
%   bnd(1:k)    : Bottom elements in left singular vectors.

% Below is a very slow replacement for the BDSQR MEX-file.

% warning('PROPACK:NotUsingMex','Using slow matlab code for bdsqr.')
k = length(alpha);
if min(size(alpha)') ~= 1  | min(size(beta)') ~= 1
  error('alpha and beta must be vectors')
elseif length(beta) ~= k
  error('alpha and beta must have the same lenght')
end
B = spdiags([alpha(:),beta(:)],[0,-1],k+1,k);
[U,S,V] = svd(full(B),0);
sigma = diag(S);
bnd = U(end,1:k)';

author: andrewssobral

repo: lrslibrary


text:
function h = fast_viewvol(vol,volres,frameres)
% h = fast_viewvol(vol,volres,frameres)

h = [];

if(nargin < 1 | nargin > 3)
  fprintf('USAGE: h = fast_viewvol(vol,volres,frameres)\n');
  return;
end

if(nargin < 2)      volres = [1 1 1]; end
if(isempty(volres)) volres = [1 1 1]; end

if(nargin < 3)        frameres = 1; end
if(isempty(frameres)) frameres = 1; end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function dx = pendulumSystem(x,A,P,grid,dyn)
% dx = pendulumSystem(x,A,P,pendulum)
%
% This function computes the closed loop dynamics of the pendulum, using
% the optimal controller from the MDP solve
%

th = x(1,:);
w = x(2,:);

u = pendulumController(x,A,P,grid);
dw = pendulumDynamics(th,w,u,dyn);

dx = [w;dw];

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [out] = percolation_1(p1,S,dt)
%percolation_1

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

% Flux function
% ------------------
% Description:  Percolation at a constant rate
% Constraints:  f <= S/dt
% @(Inputs):    p1   - base percolation rate [mm/d]
%               S    - current storage [mm]
%               dt   - time step size [d]

out = min(p1,S/dt);

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function u = poisson(varargin)
%POISSON   Poisson solver in the unit ball with Dirichlet or Neumann
%          boundary conditions.
%   U = POISSON(F, BC, M, N, P) is the solution to the Poisson equation in
%   the unit ball with right-hand side F and Dirichlet boundary data
%   U(1,lambda,theta) = @(lambda,theta) BC(lambda, theta). It uses a
%   discretization size of M x N x P.
%
%   U = POISSON(F, BC, M) is the same as POISSON(F, BC, M, M, M).
%
%   U = POISSON(F, BC, M, N, P, 'neumann') is the solution to the Poisson
%   equation with right-hand side F and Neuamnn boundary data U(1,lambda,theta) =
%   @(lambda,theta) BC(lambda, theta). It uses a discretization size of M x N x P.
%
%   U = POISSON(F, BC, M, 'neumann') is the same as POISSON(F, BC, M, M, M, 'neumann').
%
%   Also see HELMHOLTZ.

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

% Call Helmholtz command with zero frequency:
u = helmholtz(varargin{1}, 0, varargin{2:end});
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [B1, C1, c1, B2, C2, c2] = ssad(a, b, A)
%SSAD  gives both solutions to the side-side-angle problem, in degrees.
%
%   SSAD(a, b, A) will result in NaNs if the existence condition
%   |sin b * sin A| <= | sin a |  is not met.
%
%   See also SSA, MALD, MSLD.

% Rody P.S. Oldenhuis
% Delft University of Technology
% oldenhuis@gmail.com
%
% Crated     : 23/Feb/2009
% Last edited: 30/Nov/2012

    % find both solutions by calling ssa directly
    r2d = 180/pi;
    d2r = 1/r2d;
    [B1, C1, c1, B2, C2, c2] = ssa(a*d2r, b*d2r, A*d2r);
    [B1, C1, c1, B2, C2, c2] = deal(B1*r2d, C1*r2d, c1*r2d, B2*r2d, C2*r2d, c2*r2d);

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%kpamgenpattern 'shifts a pattern to yield a full scan '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros pamgenpattern.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'First Input data object'
% Integer: numX 'Number of Steps X', default: 5: 'number of steps into X direction'
% Integer: numY 'Number of Steps Y', default: 5: 'number of steps into Y direction'
% Integer: stepX 'Step X', default: 1: 'steps into X direction to jump to while stepping'
% Integer: stepY 'Step Y', default: 1: 'steps into Y direction to jump to while stepping'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kpamgenpattern(i, {'i','';'numX',5;'numY',5;'stepX',1;'stepY',1;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% pamgenpattern - shifts a pattern to yield a full scan
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1996-2003, Rainer Heintzmann,  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function result = circularShift(matrix, colshift, rowshift)
% CIRCULARSHIFT: Circular shifting of a matrix/image, i.e., pixels that get
% shifted off one side of the image are put back on the other side.
%
% result = circularShift(matrix, colshift, rowshift)
%
% EPS, DJH '96

lastrow = size(matrix, 1);
lastcol = size(matrix, 2);

result = matrix;

% Shift the cols
if (colshift>0)
  result = [result(:,[lastcol-colshift+1:lastcol]) ...
	    result(:,[1:lastcol-colshift])];
else
  colshift = -colshift;
  result = [result(:,[colshift+1:lastcol]) ...
	    result(:,[1:colshift])];
end

% Shift the rows
if (rowshift>0)
  result = [result([lastrow-rowshift+1:lastrow],:) ; ...
	    result([1:lastrow-rowshift],:)];
else
  rowshift = -rowshift;
  result = [result([rowshift+1:lastrow],:) ; ...
	    result([1:rowshift],:)];
end


text:
function cg_test ( )

%*****************************************************************************80
%
%% CG_TEST tests CG.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CG_TEST:\n' );
  fprintf ( 1, '  CG uses the Conjugate Gradient \n' );
  fprintf ( 1, '  iterative method to approximate the solution \n' );
  fprintf ( 1, '  of a linear system A * x = b.\n' );

  n = 10;

  A = zeros ( n, n );

  for i = 1 : n
    A(i,i) = 2.0;
  end
  for i = 1 : n-1
    A(i,i+1) = -1;
  end
  for i = 2 : n
    A(i-1,i) = -1;
  end
  x = [ 1 : n ]';
  b = A * x;
  x = ones ( n, 1 );

  M = 0;
  max_it = 10;
  tol = 0.0001;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  For this example, the order of the system is N = %d\n', n );
  fprintf ( 1, '  The matrix A is the simple tridiagonal -1, 2, -1.\n' );
  fprintf ( 1, '  The correct solution is x = [ 1, 2, ..., n].\n' );
  fprintf ( 1, '  The right hand side b is determined by computing A * x.\n' );
  fprintf ( 1, '  The exact x is then replaced by a vector of all 1''s for\n' );
  fprintf ( 1, '  use as a starting guess.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Other parameters are set as follows:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The maximum number of steps is %d.\n', max_it );
  fprintf ( 1, '  The error tolerance is %f\n', tol );

  [ x, error_norm, iter, flag ] = cg ( A, x, b, M, max_it, tol );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The CG routine has returned with FLAG = %d\n', flag );
  if ( flag == 0 )
    fprintf ( 1, '  This indicates that the iteration has converged.\n' );
  elseif ( flag == 1 )
    fprintf ( 1, '  This indicates that the iteration has NOT converged.\n' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of iterations taken was %d\n', iter );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The L2 norm of the error per iteration:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : iter
    fprintf ( 1, '  %4d  %f\n', i, error_norm(i) );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The computed solution vector X:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %4d  %f\n', i, x(i) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CG_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function sF = plus(sF1, sF2)
% overloads sF1 + sF2

if isa(sF2,'S2FunHarmonic')
  sF = sF2 + sF1;
  return
end

sF = S2FunHandle(@(v) sF1.eval(v) + sF2.eval(v));

end

author: mtex-toolbox

repo: mtex

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


text:
% Post processing the Finite Element Results
% Plotting the profile of components on Finite Element mesh / deformed mesh
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Warning : On running this the workspace memory will be deleted. Save if
% any data present before running the code !!
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%--------------------------------------------------------------------------
% Code written by : Siva Srinivas Kolukula                                |
%                   Senior Research Fellow                                |
%                   Structural Mechanics Laboratory                       |
%                   Indira Gandhi Center for Atomic Research              |
%                   India                                                 |
% E-mail : allwayzitzme@gmail.com                                         |
%          http://sites.google.com/site/kolukulasivasrinivas/            |
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%
% Variable descriptions
%      coordinates.dat - data file of nodal coordinates of the nodes
%      nodes.dat       - data file of elemental nodal connectivity
%      displacements   - displacement results obtained from FEA
%      factor          - amplication factor for deformed mesh
%
% NOTE : Please note that in coordinates ,displacements first column is
%        node number and in nodes forst column is element number .
%--------------------------------------------------------------------------

clear all;clc ; close all;
%--------------------------------------------------------------------------
%  input data
%--------------------------------------------------------------------------
load coordinates.dat ;
load nodes.dat ;
load displacements ;
%--------------------------------------------------------------------------
% Sorting the values accordingly
%--------------------------------------------------------------------------
UX = displacements(:,2) ;
UY = displacements(:,3) ;
UZ = displacements(:,4) ;
 U = sqrt(UX.^2+UY.^2+UZ.^2) ;
RX = displacements(:,5) ;
RY = displacements(:,6) ;
RZ = displacements(:,7) ;
%--------------------------------------------------------------------------
% Plotting the mesh and profiles
%--------------------------------------------------------------------------
PlotMesh(coordinates,nodes) ;   %  Plot the FEM mesh
%
component = UZ ;
PlotFieldonMesh(coordinates,nodes,component) ; % Plot the component profile on mesh
%
depl = [UX UY UZ] ;
factor = 90 ;                % Plot the component profile on deformed mesh
PlotFieldonDefoMesh(coordinates,nodes,factor,depl,component) ;

author: Sable

repo: mcbench-benchmarks


text:
function varargout = drawMesh(vertex, face, varargin)
% Plot polygonal 2D or 3D mesh.
%
% Usage:    drawMesh(vertex, face)
%           drawMesh(vertex, face, 'wire')
%           hMesh = drawMesh(vertex, face, color, 'surf', alpha)
%
% INPUT:
% vertex        - nv-by-2 matrix of [x y] or
%                 nv-by-3 matrix of [x y z] vertex coordinates
% face          - nf-by-? faces matrix
% (optional)
% 'wire' or
% 'surf'        - plot wire mesh or surface (default) mesh;
% color         - string, defining the face color (only in 'surf' mode);
%                 default: 'g' (green color);
% alpha         - a scalar between 0 and 1, defining transparency:
%                 0 - transparent, 1 - opaque (default).
%
% OUTPUT:
% (optional)
% hMesh         - handle to the patch object.
%
% Examples:
% 3D:
%  load('queen.mat'); figure(1); clf; drawMesh(vertex, face);
%  load('queen.mat'); figure(1); clf; drawMesh(vertex, face, 'wire');
%  load('queen.mat'); figure(1); clf; drawMesh(vertex, face, 'r', .5);
% 2D:
%  load('home.mat');  figure(1); clf; drawMesh(vertex, face, 'wire');
%
% See also: drawVector, drawPlane, drawSpan, drawLine.

% Copyright (c) 2009, Dr. Vladimir Bondarenko <http://sites.google.com/site/bondsite>

% Check input:
error(nargchk(2,4,nargin));
[mv,nv] = size(vertex);
[mf,nf] = size(face);
if (nv~=3)&&(nv~=2) , error('Wrong dimensions of the vertex matrix.'); end;

% Defaults:
alfa     = 1;
meshType = 'surf';
fcolor = 'g';

% Parse input:
if nargin > 2
    for ii=1:nargin-2
        if ischar(varargin{ii})
           if strcmpi(varargin{ii}, 'surf')||strcmpi(varargin{ii}, 'wire')
               meshType = varargin{ii};
           else
                 fcolor = varargin{ii};
           end;
        elseif isnumeric(varargin{ii})
            alfa = varargin{ii};
        end
    end
end
if (alfa > 1)||(alfa < 0), error('The alpha parameter must be in [0, 1].');end;

% Plot the mesh
switch meshType
    case 'surf'
        hMesh = patch('vertices', vertex,'faces', face,...
                      'facecolor', fcolor, 'FaceAlpha', alfa);
        camlight('headlight');
    case 'wire'
        hMesh = patch('vertices', vertex,'faces', face,...
                      'facecolor', 'none');
end
%axis equal tight off
if nargout==1, varargout{1} = hMesh; end


text:
function pink_noise_test ( )

%*****************************************************************************80
%
%% PINK_NOISE_TEST tests the PINK_NOISE library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 May 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PINK_NOISE_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the PINK_NOISE library.\n' );

  pink_noise_test01 ( );
  pink_noise_test02 ( );
  pink_noise_test03 ( );
  pink_noise_test04 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PINK_NOISE_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function alpha=solveAlpha(I,consts_map,consts_vals,varargin)

  [h,w,c]=size(I);
  img_size=w*h;

  A=getLaplacian1(I,consts_map,varargin{:});

  D=spdiags(consts_map(:),0,img_size,img_size);
  lambda=100;
  x=(A+lambda*D)\(lambda*consts_map(:).*consts_vals(:));

  alpha=max(min(reshape(x,h,w),1),0);

author: luanfujun

repo: deep-photo-styletransfer


text:
function varargout = holdertable(X)
% Holder table function
%
%   HOLDERTABLE([x1, x2]) returns the value of the Holder table
%   function at the specified points. [x1] and [x2] may be vectors.
%   The search domain is
%
%               -10 < x_i < 10
%
%   The four global minima are near the edges of the interval, and have a
%   function value of
%
%       f(x*) = -1.92085026.

% Author: Rody P.S. Oldenhuis
% Delft University of Technology
% E-mail: oldenhuis@dds.nl
% Last edited 20/Jul/2009

    % if no input is given, return dimensions, bounds and minimum
    if (nargin == 0)
        varargout{1} = 2;  % # dims
        varargout{2} = [-10, -10]; % LB
        varargout{3} = [+10, +10]; % UB
        varargout{4} = [+8.055023472141116e+000   +9.664590028909654e+000
                        -8.055023472141116e+000   +9.664590028909654e+000
                        +8.055023472141116e+000   -9.664590028909654e+000
                        -8.055023472141116e+000   -9.664590028909654e+000]; % solution
        varargout{5} = -1.920850256788675e+001; % function value at solution

    % otherwise, output function value
    else

        % keep values within the search interval
        X(X < -10) = inf;      X(X > 10) = inf;

        % split input vector X into x1, x2
        if size(X, 1) == 2
            x1 = X(1, :);        x2 = X(2, :);
        else
            x1 = X(:, 1);        x2 = X(:, 2);
        end

        % output function value
        varargout{1} = -abs(sin(x1).*cos(x2).*exp(abs(1 - sqrt(x1.^2 + x2.^2)/pi)));

    end

end


text:
function [x, y] = theta2xy(theta1, theta2, l1, l2)
% Angle to Position Conversion for 2 link SCARA robot
%
% [Inputs]
%	theta1      : angle1 [deg]
%   theta2      : angle2 [deg]
%	l1          : link1 length [m]
%	l2          : link2 length [m]
% [Outputs]
%	x		    : x coordinate of the edge of link2 [m]
%   y			: y coordinate of the edge of link2 [m]

% deg to rad conversion
theta1 = theta1 * pi / 180;
theta2 = theta2 * pi / 180;

% position
x = l1 * cos(theta1) + l2 * cos(theta1 + theta2);
y = l1 * sin(theta1) + l2 * sin(theta1 + theta2);


text:
function ft_realtime_fmriproxy(cfg)

% FT_REALTIME_FMRIPROXY simulates an fMRI acquisition system by writing volumes in a
% cycle of about 2 seconds. The voxel data is written as a column vector with X*Y*Z
% channels, where X and Y are the readout and phase-encoding resolution, and Z is the
% number of slices. The voxel data consists of a ellipsoid (a bit like a head) with
% added lateralized activation (20 scan cycle) and noise.
%
% This function also writes out events of type='scan' and value='pulse' when the
% simulated scanner initiates a scan, and value='ready' when a hypothetical
% processing pipeline is finished with that scan, just after writing out the volume
% data itself. There is an artificial delay of 1.3*TR between the two events.
%
% Use as
%   ft_realtime_fmriproxy(cfg)
%
% The target to write the data to is configured as
%   cfg.target.datafile      = string, target destination for the data (default = 'buffer://localhost:1972')
%
% You can also select a resolution of the simulated volumes (default = [64,64,20]) like
%   cfg.voxels = [64 64 32]
% and the repetition time (TR, default = 0.08*number of slices) in seconds using
%   cfg.TR = 2.0
%
% To stop this realtime function, you have to press Ctrl-C
%
% See also FT_REALTIME_SIGNALPROXY, FT_REALTIME_SIGNALVIEWER

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

% set the defaults
if isempty(cfg) || ~isfield(cfg, 'target') || ~isfield(cfg.target, 'datafile')
  cfg.target.datafile = 'buffer://localhost:1972';
end
cfg.target.dataformat = [];

if isfield(cfg, 'voxels')
   voxels = cfg.voxels(1:3);
else
   voxels = [64 64 20];
end

if isfield(cfg,'TR')
   TR = cfg.TR;
   Tslice = TR / voxels(3);
else
   Tslice = 0.08;
   TR = Tslice * voxels(3);
end

nifti = [];
nifti.dim = voxels;
nifti.pixdim = [3.5 3.5 3.5]; % size of voxel in mm
nifti.slice_duration = Tslice;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create a FieldTrip compatible header structure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hdr = [];
hdr.nChans = prod(voxels);
hdr.nSamples = 0;
hdr.Fs = 1/TR;
hdr.nSamplesPre        = 0;
hdr.nTrials            = 1;
hdr.nifti_1 = encode_nifti1(nifti);

% prepare two event structure that will go along with every sample/scan
evp = [];
evp.type = 'scan';
evp.value = 'pulse';
evp.offset = 0;
evp.duration = 0;
evp.sample = 0;

evr = evp;
evr.value = 'ready';

% assume acquisition + processing takes 30% of a TR cycle
TP = 0.3 * TR;

% make some fake MRI data
% Sphere = volume containing a spherical "head"
% Left/Right = volume of same size, adds a bit of fake activation
x = linspace(-1,1,voxels(1));
y = linspace(-1,1,voxels(2));
z = linspace(-1,1,voxels(3));
[X,Y,Z] = meshgrid(x,y,z);

R = sqrt(X.^2 + Y.^2 + Z.^2);
Sphere = single(1200 * (1./(1 + exp(40*(R-0.8)))));
R = sqrt((X-0.4).^2 + Y.^2 + Z.^2);
Left  = single(60 * (1./(1 + exp(40*(R-0.4)))));
R = sqrt((X+0.4).^2 + Y.^2 + Z.^2);
Right = single(60 * (1./(1 + exp(40*(R-0.4)))));

stopwatch = tic;

% write header to FieldTrip buffer
ft_write_data(cfg.target.datafile, single([]), 'header', hdr, 'dataformat', cfg.target.dataformat, 'append', false);

numPulse = 0;
numScan = 0;

while true
	% simulate acquisition of pixeldata
	t0 = toc(stopwatch);

	% write 'scan pulse' event
	numPulse = numPulse + 1;
	evp.sample = numPulse;
	fprintf('Sending pulse: %3i, clock time = %f\n', numPulse, t0);
	ft_write_event(cfg.target.datafile, evp);

	if numPulse > 1
		% fake acquisition + processing
		lateral = 0.5 + 0.5*sin(2*pi*numScan/20);

		Noise = single(30*randn(voxels));
		scan = single(Sphere + lateral*Left + (1-lateral)*Right + Noise);

		tp = toc(stopwatch);
		pause(TP - (tp-t0));

		ft_write_data(cfg.target.datafile, scan(:), 'append', true);
		% write 'scan ready' event
		numScan = numScan + 1;
		evr.sample = numScan;
		fprintf('Scans written: %3i, clock time = %f\n', numScan, toc(stopwatch));
		ft_write_event(cfg.target.datafile, evr);
	end

	tr = toc(stopwatch);
	pause(numPulse*TR - tr);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef ALGORITHM < handle & matlab.mixin.Heterogeneous
%ALGORITHM - The superclass of algorithms.
%
%   This is the superclass of algorithms. An object of ALGORITHM stores the
%   settings of the algorithm and the data generated in current execution.
%
% ALGORITHM properties:
%   parameter       <any>       parameters of the algorithm
%   save            <scalar>    number of populations saved in an execution
%   outputFcn       <function>	function called after each generation
%   pro             <class>     problem solved in current execution
%   result          <cell>      populations saved in current execution
%   metric          <struct>    metric values of current populations
%
% ALGORITHM methods:
%   ALGORITHM       <protected> the constructor setting all the properties specified by user
%   Solve           <public>    use the algorithm to solve a problem
%   main            <public>	the main function of the algorithm
%   NotTerminated 	<protected>	the function called after each generation of the execution
%   ParameterSet	<protected>	obtain the parameters of the algorithm

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties(SetAccess = protected)
        parameter = {};                 % Parameters of the algorithm
        save      = -10;            	% Number of populations saved in an execution
        outputFcn = @DefaultOutput;     % Function called after each generation
        pro;                            % Problem solved in current execution
        result;                         % Populations saved in current execution
        metric;                         % Metric values of current populations
    end
    methods(Access = protected)
        function obj = ALGORITHM(varargin)
        %ALGORITHM - The constructor of ALGORITHM.
        %
        %   Alg = algName('Name',Value,'Name',Value,...) generates an
        %   object with the properties specified by the inputs. algName is
        %   a subclass of ALGORITHM, while ALGORITHM cannot be instantiated
        %   directly.
        %
        %   Example:
        %        Algorithm = MOEAD('parameter',4,'save',1)

            isStr = find(cellfun(@ischar,varargin(1:end-1))&~cellfun(@isempty,varargin(2:end)));
            for i = isStr(ismember(varargin(isStr),{'parameter','save','outputFcn'}))
                obj.(varargin{i}) = varargin{i+1};
            end
        end
    end
    methods(Sealed)
        function Solve(obj,Problem)
        %Solve - Use the algorithm to solve a problem.
        %
        %   obj.Solve(Pro) uses the algorithm to solve a problem, where Pro
        %   is a PROBLEM object.
        %
        %   In terms of the default obj.outputFcn, the result will be
        %   displayed when obj.save = 0 and saved when obj.save > 0.
        %
        %   Example:
        %       Algorithm.Solve(Problem)

            try
                obj.result = {};
                obj.metric = struct('runtime',0);
                obj.pro    = Problem;
                obj.pro.FE = 0;
                addpath(fileparts(which(class(obj))));
                addpath(fileparts(which(class(obj.pro))));
                tic; obj.main(obj.pro);
            catch err
                if ~strcmp(err.identifier,'PlatEMO:Termination')
                    rethrow(err);
                end
            end
        end
    end
    methods
        function main(obj,Problem)
        %main - The main function of the algorithm.
        %
        %   This function is expected to be implemented in each subclass of
        %   ALGORITHM, which is usually called by ALGORITHM.Solve.
        end
    end
    methods(Access = protected, Sealed)
        function nofinish = NotTerminated(obj,Population)
        %NotTerminated - The function called after each generation of the
        %execution.
        %
        %   obj.NotTerminated(P) stores the population P as the result of
        %   the current execution, and returns true if the algorithm should
        %   be terminated, i.e., the number of function evaluations or
        %   runtime exceeds.
        %
        %   obj.outputFcn is called here, whose runtime will not be counted
        %   in the runtime of current execution.
        %
        %   Example:
        %       while Algorithm.NotTerminated(Population)
        %           ... ...
        %       end

            obj.metric.runtime = obj.metric.runtime + toc;
            if obj.pro.maxRuntime < inf
                obj.pro.maxFE = obj.pro.FE*obj.pro.maxRuntime/obj.metric.runtime;
            end
            num   = max(1,abs(obj.save));
            index = max(1,min(min(num,size(obj.result,1)+1),ceil(num*obj.pro.FE/obj.pro.maxFE)));
            obj.result(index,:) = {obj.pro.FE,Population};
            drawnow('limitrate');
            obj.outputFcn(obj,obj.pro);
            nofinish = obj.pro.FE < obj.pro.maxFE;
            assert(nofinish,'PlatEMO:Termination',''); tic;
        end
        function varargout = ParameterSet(obj,varargin)
        %ParameterSet - Obtain the parameters of the algorithm.
        %
        %   [p1,p2,...] = obj.ParameterSet(v1,v2,...) sets the values of
        %   parameters p1, p2, ..., where each parameter is set to the
        %   value given in obj.parameter if obj.parameter is specified, and
        %   set to the value given in v1, v2, ... otherwise.
        %
        %   Example:
        %       [p1,p2,p3] = Algorithm.ParameterSet(1,2,3)

            varargout = varargin;
            specified = ~cellfun(@isempty,obj.parameter);
            varargout(specified) = obj.parameter(specified);
        end
    end
    methods(Sealed)
        function Scores = CalMetric(obj,metName)
        % Calculate metric values

            if ~isfield(obj.metric,metName)
                obj.metric.(metName) = [cell2mat(obj.result(:,1)),cellfun(@(S)obj.pro.CalMetric(metName,S),obj.result(:,2))];
            end
            Scores = obj.metric.(metName);
        end
    end
end

function DefaultOutput(Algorithm,Problem)
% The default output function of ALGORITHM

    clc; fprintf('%s on %d-objective %d-variable %s (%6.2f%%), %.2fs passed...\n',class(Algorithm),Problem.M,Problem.D,class(Problem),Problem.FE/Problem.maxFE*100,Algorithm.metric.runtime);
    if Problem.FE >= Problem.maxFE
        if Algorithm.save < 0
            if Problem.M > 1
                Population = Algorithm.result{end};
                if length(Population) >= size(Problem.optimum,1); name = 'HV'; else; name = 'IGD'; end
                value = Algorithm.CalMetric(name);
                figure('NumberTitle','off','Name',sprintf('%s : %.4e  Runtime : %.2fs',name,value(end),Algorithm.CalMetric('runtime')));
                title(sprintf('%s on %s',class(Algorithm),class(Problem)),'Interpreter','none');
                top = uimenu(gcf,'Label','Data source');
                g   = uimenu(top,'Label','Population (obj.)','CallBack',{@(h,~,Pro,P)eval('Draw(gca);Pro.DrawObj(P);cb_menu(h);'),Problem,Population});
                uimenu(top,'Label','Population (dec.)','CallBack',{@(h,~,Pro,P)eval('Draw(gca);Pro.DrawDec(P);cb_menu(h);'),Problem,Population});
                uimenu(top,'Label','True Pareto front','CallBack',{@(h,~,P)eval('Draw(gca);Draw(P,{''\it f\rm_1'',''\it f\rm_2'',''\it f\rm_3''});cb_menu(h);'),Problem.optimum});
                cellfun(@(s)uimenu(top,'Label',s,'CallBack',{@(h,~,A)eval('Draw(gca);Draw(A.CalMetric(h.Label),''-k.'',''LineWidth'',1.5,''MarkerSize'',10,{''Number of function evaluations'',strrep(h.Label,''_'','' ''),[]});cb_menu(h);'),Algorithm}),{'IGD','HV','GD','Feasible_rate'});
                set(top.Children(4),'Separator','on');
                g.Callback{1}(g,[],Problem,Population);
            else
                best = Algorithm.CalMetric('Min_value');
                if isempty(best); best = nan; end
                figure('NumberTitle','off','Name',sprintf('Min value : %.4e  Runtime : %.2fs',best(end),Algorithm.CalMetric('runtime')));
                title(sprintf('%s on %s',class(Algorithm),class(Problem)),'Interpreter','none');
                top = uimenu(gcf,'Label','Data source');
                uimenu(top,'Label','Population (dec.)','CallBack',{@(h,~,Pro,P)eval('Draw(gca);Pro.DrawDec(P);cb_menu(h);'),Problem,Algorithm.result{end}});
                cellfun(@(s)uimenu(top,'Label',s,'CallBack',{@(h,~,A)eval('Draw(gca);Draw(A.CalMetric(h.Label),''-k.'',''LineWidth'',1.5,''MarkerSize'',10,{''Number of function evaluations'',strrep(h.Label,''_'','' ''),[]});cb_menu(h);'),Algorithm}),{'Min_value','Feasible_rate'});
                set(top.Children(2),'Separator','on');
                top.Children(2).Callback{1}(top.Children(2),[],Algorithm);
            end
        elseif Algorithm.save > 0
            folder = fullfile('Data',class(Algorithm));
            [~,~]  = mkdir(folder);
            file   = fullfile(folder,sprintf('%s_%s_M%d_D%d_',class(Algorithm),class(Problem),Problem.M,Problem.D));
            runNo  = 1;
            while exist([file,num2str(runNo),'.mat'],'file') == 2
                runNo = runNo + 1;
            end
            result = Algorithm.result;
            metric = Algorithm.metric;
            save([file,num2str(runNo),'.mat'],'result','metric');
        end
    end
end

function cb_menu(h)
% Switch between the selected menu

    set(get(get(h,'Parent'),'Children'),'Checked','off');
    set(h,'Checked','on');
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function tfrscalt
%TFRSCALT Unit test for the time-frequency representation TFRSCALO.

%	O. Lemoine - June 1996.

% We test each property of the corresponding TFR :

N=128;

% Covariance by translation in time
t1=60; t2=70; f=0.3; W=0;
sig1=amgauss(N,t1).*fmconst(N,f,t1);
sig2=amgauss(N,t2).*fmconst(N,f,t2);
tfr1=tfrscalo(sig1,1:N,W,0.1,0.4,128);
tfr2=tfrscalo(sig2,1:N,W,0.1,0.4,128);
[tr,tc]=size(tfr1);
nu=round(f*(tc-1)*2)+1;
tfr=tfr1-tfr2(:,modulo((1:tc)-t1+t2,tc));
if any(any(abs(tfr)>sqrt(eps))),
 error('tfrscalo test 1 failed');
end

% Covariance by dilation
t=N/2; f=0.2; T=2*sqrt(N); a=2; W=8;
sig1=amgauss(N,t,T).*fmconst(N,f,t);
sig2=amgauss(a*N,a*t,T*a).*fmconst(a*N,f/a,a*t);
[tfr1,t1,f1]=tfrscalo(sig1,1:N  ,W,0.01,0.49,N);
[tfr2,t2,f2]=tfrscalo(sig2,1:a*N,W,0.01,0.49,N);
Max1=max(max(tfr1)); Max2=max(max(tfr2));
[I1,J1]=find(tfr1==Max1); [I2,J2]=find(tfr2==Max2);
if abs(f1(I1)-a*f2(I2))>1e-2 | J2~=a*J1,
 error('tfrscalo test 2 failed');
end

% Reality of the TFR
sig=noisecg(N); W=5;
tfr=tfrscalo(sig,1:N,W,0.01,0.5,N);
if sum(any(abs(imag(tfr))>sqrt(eps)))~=0,
 error('tfrscalo test 3 failed');
end

% Energy conservation
sig=fmsin(N,.1,.4); W=6; Nf=2*N ;
[tfr,t,f]=tfrscalo(sig,1:N,W,0.01,0.49,Nf);
Es=norm(sig)^2/Nf;
Etfr=integ2d(tfr,t,f)/N;
if abs(Es-Etfr)>sqrt(eps),
 error('tfrscalo test 4 failed');
end

% Positivity
if any(any(tfr<0)),
 error('tfrscalo test 5 failed');
end

% Same energy in the time-scale plane for 2 gaussian atoms at different scales
sig=amgauss(256).*(fmconst(256,.15)+fmconst(256,.35));
[tfr,t,f]=tfrscalo(sig,1:256,12,.01,.49,512);
int1=integ2d(tfr(1:430,:),t,f(1:430));
int2=integ2d(tfr(431:512,:),t,f(431:512));
if abs(int1-int2)>1e-4,
 error('tfrscalo test 6 failed');
end

N=127;

% Covariance by dilation
t=round(N/2); f=0.2; T=2*sqrt(N); a=2; W=8;
sig1=amgauss(N,t,T).*fmconst(N,f,t);
sig2=amgauss(a*N,a*t,T*a).*fmconst(a*N,f/a,a*t);
[tfr1,t1,f1]=tfrscalo(sig1,1:N  ,W,0.01,0.49,N);
[tfr2,t2,f2]=tfrscalo(sig2,1:a*N,W,0.01,0.49,N);
Max1=max(max(tfr1)); Max2=max(max(tfr2));
[I1,J1]=find(tfr1==Max1); [I2,J2]=find(tfr2==Max2);
if abs(f1(I1(1))-a*f2(I2))>1e-2 | J2~=a*J1(1)+1,
 error('tfrscalo test 7 failed');
end

% Reality of the TFR
sig=noisecg(N); W=5;
tfr=tfrscalo(sig,1:N,W,0.01,0.5,N);
if sum(any(abs(imag(tfr))>sqrt(eps)))~=0,
 error('tfrscalo test 8 failed');
end

% Energy conservation
sig=fmsin(N,.1,.4); W=6; Nf=2*N+1 ;
[tfr,t,f]=tfrscalo(sig,1:N,W,0.01,0.49,Nf);
SP = fft(hilbert(sig));
indmin = 1+round(0.01*(N-2));
indmax = 1+round(0.49*(N-2));
SPana = SP(indmin:indmax);
Es=SPana'*SPana/Nf;
Etfr=integ2d(tfr,t,f);
if abs(Es-Etfr)>1e-3,
 error('tfrscalo test 9 failed');
end

% Positivity
if any(any(tfr<0)),
 error('tfrscalo test 10 failed');
end

% Same energy in the time-scale plane for 2 gaussian atoms at different scales
sig=amgauss(N).*(fmconst(N,.15)+fmconst(N,.35));
[tfr,t,f]=tfrscalo(sig,1:N,12,.01,.49,2*N+1);
int1=integ2d(tfr(1:210,:),t,f(1:210));
int2=integ2d(tfr(211:255,:),t,f(211:255));
if abs(int1-int2)>5e-4,
 error('tfrscalo test 11 failed');
end

author: HeLiangHIT


text:
% Reindexes whole cell simulation.
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 6/30/2011
function runReindexing(outputDirectory)
setWarnings();
setPath();
setPreferences();

import edu.stanford.covert.cell.sim.util.DiskLogger;

DiskLogger.reindexTimeCourses(outputDirectory);

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function bnet = bnet_from_engine(engine)
% BNET_FROM_ENGINE Return the bnet structure stored inside the engine (smoother_engine)
% bnet = bnet_from_engine(engine)

bnet = bnet_from_engine(engine.tbn_engine);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%% CHUTE  Chutes and ladders analysis.
%
%  Discussion:
%
%    This script determines the probability of finishing a one-player
%    game of chutes and ladders in exactly N moves (PDF) or in no more
%    than N moves (CDF).

%
%    The game starts at square 0, and ends when square 100 is reached.
%
%  Modified:
%
%    19 September 2014
%
%  Author:

%

%    Desmond Higham, Nicholas Higham

%

%  Reference:

%

%    Desmond Higham, Nicholas Higham,

%    MATLAB Guide,

%    SIAM, 2005,

%    ISBN13: 9780898717891.

%
  N = 100;
%
% "+1" translates square to state.
%
  top = [ 1  4  9 16 21 28 36 47 49 51 56 62 64 71 80  87 93 95 98] + 1;
  bot = [38 14 31  6 42 84 44 26 11 67 53 19 60 91 100 24 73 75 78] + 1;

  P = toeplitz(zeros(1,N+1),[0 ones(1,6) zeros(1,N-6)]);

  for k = N-4:N+1
    P(k,k) = k-N+5;
  end
  P = P/6;

  for k = 1:length(top)
    r = top(k); s = bot(k);     % Chute or ladder from r to s.
    P(:,s) = P(:,s) + P(:,r);   % Add column r to column s.
  end
  P(top,:) = [];
  P(:,top) = [];  % Remove starts of chutes and ladders.

  figure(1)
  spy(P)

  M = 200;
  cumprob = zeros(M,1);
  cumprob(1) = P(1,end);
  v = P(1,:);
  for n = 2:M,
    v = v*P;
    cumprob(n) = v(end);
  end

  figure(2)
  colormap([0.8,0.4,0.4])
  bar(diff([0;cumprob]))
  title('Probability for Game Length','FontSize',12,'FontWeight','Bold')
  grid on
  xlim([0 M])

figure(3)

colormap([0.8,0.4,0.4])
  bar(cumprob)
  title('Cumulative Probability for Game Length',...
      'FontSize',12,'FontWeight','Bold')
  grid on
  xlim([0 M])


text:
function a = chow_inverse ( alpha, beta, n )

%*****************************************************************************80
%
%% CHOW_INVERSE returns the inverse of the CHOW matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, the ALPHA value.  A typical value is 1.0.
%
%    Input, real BETA, the BETA value.  A typical value is 0.0.
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  if ( 0.0 == alpha & beta == 0.0 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'CHOW_INVERSE - Fatal error!\n' );
    fprintf ( 1, '  The Chow matrix is not invertible, because\n' );
    fprintf ( 1, '  ALPHA = 0 and BETA = 0.\n' );
    error ( 'CHOW_INVERSE - Fatal error!' );

  elseif ( 0.0 == alpha & beta ~= 0.0 )

    for i = 1 : n
      for j = 1 : n

        if ( i <= j )
          a(i,j) = (-1)^(j-i) / beta^(j+1-i);
        else
          a(i,j) = 0.0;
        end

      end
    end

    return

  elseif ( 0.0 ~= alpha & beta == 0.0 )

    if ( 1 < n )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CHOW_INVERSE - Fatal error!\n' );
      fprintf ( 1, '  The Chow matrix is not invertible, because\n' );
      fprintf ( 1, '  BETA = 0 and 1 < N.\n' );
      error ( 'CHOW_INVERSE - Fatal error!' );
    end

    a(1,1) = 1.0 / alpha;

    return

  end

  d(0+1) = 1.0;
  d(1+1) = beta;
  for i = 2 : n
    d(i+1) = beta * d(i-1+1) + alpha * beta * d(i-2+1);
  end

  dp(-1+2) = 1.0 / beta;
  dp(0+2) = 1.0;
  dp(1+2) = alpha + beta;
  for i = 2 : n
    dp(i+2) = d(i+1) + alpha * d(i-1+1);
  end

  for i = 1 : n
    for j = 1 : i-1
      a(i,j) = - alpha * ( alpha * beta )^(i-j) * dp(j-2+2) * d(n-i+1) ...
        / dp(n+2);
    end
    for j = i : n
      a(i,j) = (-1)^(i+j) * dp(i-1+2) * d(n+1-j+1) / ( beta * dp(n+2) );
    end
  end

  return
end


text:
function c8mat_identity_test ( )

%*****************************************************************************80
%
%% C8MAT_IDENTITY_TEST tests C8MAT_IDENTITY.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2015
%
%  Author:
%
%    John Burkardt
%
  m = 4;
  n = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8MAT_IDENTITY_TEST\n' );
  fprintf ( 1, '  C8MAT_IDENTITY returns the complex identity matrix.\n' );

  a = c8mat_identity ( n );

  c8mat_print ( m, n, a, '  The identity matrix:' );

  return
end


text:
function source = ft_read_cifti(filename, varargin)

% FT_READ_CIFTI read functional data or functional connectivity from a cifti-1 or
% cifti-2 file. The functional data can consist of a dense or a parcellated
% representation. The geometrical description of the brainordinates can consist of
% triangulated surfaces or voxels in a regular 3-D volumetric grid. If available,
% it also reads the geometrical description of the surfaces from the accompanying
% gifti files.
%
% Use as
%   data = ft_read_cifti(filename, ...)
%
% If the file contains a dense representation of functional data, the output data
% structure is organized according to the FT_DATATYPE_SOURCE or FT_DATATYPE_VOLUME
% definition.
%
% If the contains a parcellated representation of functional data, the output data
% structure is organized according to the FT_DATATYPE_TIMELOCK or FT_DATATYPE_FREQ
% definition. In addition, the description of the geometry wil be represented in a
% data.brainordinate field, which is organized according to the FT_DATATYPE_SOURCE
% or FT_DATATYPE_VOLUME definition.
%
% Any optional input arguments should come in key-value pairs and may include
%   'readdata'         = boolean, can be false or true (default depends on file size)
%   'readsurface'      = boolean, can be false or true (default = true)
%   'cortexleft'       = string, filename with left cortex (optional, default is automatic)
%   'cortexright'      = string, filename with right cortex (optional, default is automatic)
%   'hemisphereoffset' = number, amount in milimeter to move the hemispheres apart from each other (default = 0)
%   'mapname'          = string, 'field' to represent multiple maps separately, or 'array' to represent as array (default = 'field')
%   'debug'            = boolean, write a debug.xml file (default = false)
%
% See also FT_WRITE_CIFTI, FT_READ_MRI, FT_WRITE_MRI

% Known limitations and bugs
% - it will fail if multiple MatrixIndicesMap contain BrainModels
% - fibers (i.e. dfan and dfibersamp) are unsupported/untested
% - metadata is unsupported

% Copyright (C) 2013-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$

readdata         = ft_getopt(varargin, 'readdata', []);   % the default depends on file size, see below
readsurface      = ft_getopt(varargin, 'readsurface', true);
cortexleft       = ft_getopt(varargin, 'cortexleft', {});
cortexright      = ft_getopt(varargin, 'cortexright', {});
hemisphereoffset = ft_getopt(varargin, 'hemisphereoffset', 0); % in mm, move the two hemispheres apart from each other
debug            = ft_getopt(varargin, 'debug', false);
mapname          = ft_getopt(varargin, 'mapname', 'field');
dataformat       = ft_getopt(varargin, 'dataformat', []);

if ft_filetype(filename, 'compressed')
  % the file is compressed, unzip on the fly
  inflated = true;
  origfile = filename;
  filename = inflate_file(filename);
else
  inflated = false;
  origfile = filename;
end

% read the header section
hdr = read_nifti2_hdr(filename);

% xml_offset = 540+12;
% xml_size   = hdr.vox_offset-xml_offset-8;

fid = fopen_or_error(filename, 'rb', hdr.endian);

% determine the file size, this is used to catch endian errors
fseek(fid, 0, 'eof');
filesize = ftell(fid);
fseek(fid, 0, 'bof');

if isempty(readdata)
  % set the default for readdata
  if filesize>1e9
    ft_warning('Not reading data by default in case filesize>1GB. Please specify the ''readdata'' option.');
    readdata = false;
  else
    readdata = true;
  end
else
  % convert 'yes'/'no' into boolean
  readdata = istrue(readdata);
end

fseek(fid, 540, 'bof');
hdrext = fread(fid, [1 4], 'int8');
if hdrext(1)~=1
  ft_error('cifti requires a header extension');
end

% determine the size of the header extension
esize = fread(fid, 1, 'int32=>int32');
etype = fread(fid, 1, 'int32=>int32');

hdrsize = 540;
voxsize = filesize-hdr.vox_offset;
if esize>(filesize-hdrsize-voxsize)
  ft_warning('the endianness of the header extension is inconsistent with the nifti-2 header');
  esize = swapbytes(esize);
  etype = swapbytes(etype);
end

if etype~=32 && etype~=swapbytes(int32(32)) % FIXME there is an endian problem
  ft_error('the header extension type is not cifti');
end

% read the extension content, subtract the 8 bytes from esize and etype
xmldata = fread(fid, [1 esize-8], 'uint8=>char');

% the size of the extension must be an integer multiple of 16 bytes according to http://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/extension.html
% consequently there might be some zero-padding at the end of the XML section
if any(xmldata==0)
  xmldata = xmldata(xmldata>0);
end

if debug
  try
    % write the xml section to a temporary file
    xmlfile = 'debug.xml';
    tmp = fopen_or_error(xmlfile, 'w');
    fwrite(tmp, xmldata);
    fclose(tmp);
  end
end

% ensure that the external toolbox is present, this adds gifti/@xmltree
ft_hastoolbox('gifti', 1);

% convert the character data to an xmltree object
tree = xmltree(xmldata);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert the xmltree object to a generic structure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

numericAttributeTypes = {'NumberOfMatrices', 'AppliesToMatrixDimension', 'IndexOffset', 'IndexCount', 'SurfaceNumberOfNodes', 'VolumeDimensions', 'SurfaceNumberOfVertices', 'SeriesStart', 'SeriesStep', 'NumberOfSeriesPoints', 'SeriesExponent', 'Vertices', 'MeterExponent'};

Cifti             = struct(); % the parent of the XML tree, it only contains version info
MatrixIndicesMap  = [];
Parcel            = [];
NamedMap          = [];
Surface           = [];
Volume            = [];
BrainModel        = struct([]); % this will be constructed both for dense and parcellated files

attr = attributes(tree, 'get', 1);
if ~iscell(attr), attr = {attr}; end % treat one attribute just like multiple attributes
for j=1:numel(attr)
  if any(strcmp(attr{j}.key, numericAttributeTypes))
    Cifti.(attr{j}.key) = str2num(attr{j}.val);
  else
    Cifti.(attr{j}.key) = attr{j}.val;
  end
end

uid_MatrixIndicesMap = find(tree,'/CIFTI/Matrix/MatrixIndicesMap');
for i=1:length(uid_MatrixIndicesMap)
  map = branch(tree, uid_MatrixIndicesMap(i));

  % get the attributes of each map
  attr = attributes(map, 'get', 1);
  for j=1:numel(attr)
    if any(strcmp(attr{j}.key, numericAttributeTypes))
      MatrixIndicesMap(i).(attr{j}.key) = str2num(attr{j}.val);
    else
      MatrixIndicesMap(i).(attr{j}.key) = attr{j}.val;
    end
  end

  switch Cifti.Version
    case {'1' '1.0'}
      uid_Volume = find(tree,'/CIFTI/Matrix/Volume');
    case {'2' '2.0'}
      uid_Volume = find(map,'/MatrixIndicesMap/Volume');
  end
  % the following will fail if there are multiple volumes
  if ~isempty(uid_Volume)
    switch Cifti.Version
      case {'1' '1.0'}
        volume = branch(tree, uid_Volume);
      case {'2' '2.0'}
        volume = branch(map, uid_Volume);
    end
    attr = attributes(volume, 'get', 1); % there should only be one attribute here
    if ~iscell(attr), attr = {attr}; end % treat one attribute just like multiple attributes
    for j=1:numel(attr)
      if any(strcmp(attr{j}.key, numericAttributeTypes))
        Volume.(attr{j}.key) = str2num(attr{j}.val);
      else
        Volume.(attr{j}.key) = attr{j}.val;
      end
    end
    uid_Transform = find(volume,'/Volume/TransformationMatrixVoxelIndicesIJKtoXYZ');
    if ~isempty(uid_Transform)
      transform = branch(volume, uid_Transform);
      attr = attributes(transform, 'get', 1);
      if isstruct(attr), attr = {attr}; end % treat one attribute just like multiple attributes
      for j=1:numel(attr)
        if any(strcmp(attr{j}.key, numericAttributeTypes))
          Volume.(attr{j}.key) = str2num(attr{j}.val);
        else
          Volume.(attr{j}.key) = attr{j}.val;
        end
      end
      Volume.Transform = str2num(get(transform, 2, 'value'));
      Volume.Transform = reshape(Volume.Transform, [4 4])'; % it needs to be transposed
    end
  end

  uid_Surface = find(map,'/MatrixIndicesMap/Surface');
  if ~isempty(uid_Surface)
    for j=1:length(uid_Surface)
      surface = branch(map, uid_Surface(j));

      % get the attributes of each surface
      attr = attributes(surface, 'get', 1);
      if isstruct(attr), attr = {attr}; end % treat one attribute just like multiple attributes
      for k=1:numel(attr)
        if any(strcmp(attr{k}.key, numericAttributeTypes))
          Surface(j).(attr{k}.key) = str2num(attr{k}.val);
        else
          Surface(j).(attr{k}.key) = attr{k}.val;
        end
      end % for

    end
  end

  uid_Parcel = find(map,'/MatrixIndicesMap/Parcel');
  if ~isempty(uid_Parcel)
    for j=1:length(uid_Parcel)
      parcel = branch(map, uid_Parcel(j));

      % get the attributes of each parcel
      attr = attributes(parcel, 'get', 1);
      if isstruct(attr), attr = {attr}; end % treat one attribute just like multiple attributes
      for k=1:numel(attr)
        if any(strcmp(attr{k}.key, numericAttributeTypes))
          Parcel(j).(attr{k}.key) = str2num(attr{k}.val);
        else
          Parcel(j).(attr{k}.key) = attr{k}.val;
        end
      end % for

      uid_VoxelIndicesIJK = find(parcel, '/Parcel/VoxelIndicesIJK');
      if ~isempty(uid_VoxelIndicesIJK)
        tmp = str2num(get(parcel, children(parcel, uid_VoxelIndicesIJK), 'value'));
        Parcel(j).VoxelIndicesIJK = reshape(tmp, 3, [])' + 1; % transpose, one offset
      else
        Parcel(j).VoxelIndicesIJK = [];
      end

      Parcel(j).BrainStructure = {};

      uid_Vertices = find(parcel, '/Parcel/Vertices');
      for k=1:length(uid_Vertices)
        vertices = branch(parcel, uid_Vertices(k));

        attr = attributes(vertices, 'get', 1);
        if isstruct(attr), attr = {attr}; end % treat one attribute just like multiple attributes
        for l=1:numel(attr)
          switch attr{l}.key
            case 'BrainStructure'
              Parcel(j).BrainStructure{k} = attr{l}.val;
          end
        end % for

        Parcel(j).Vertices{k} = str2num(get(vertices, children(vertices, find(vertices, 'Vertices')), 'value')) + 1; % one offset
      end
    end
  end

  uid_NamedMap = find(map, '/MatrixIndicesMap/NamedMap');
  if ~isempty(uid_NamedMap)
    for j=1:length(uid_NamedMap)

      namedmap = branch(map, uid_NamedMap(j));
      NamedMap(j).MapName = get(namedmap, children(namedmap, find(namedmap, '/NamedMap/MapName')), 'value');

      uid_LabelTable = find(namedmap, '/NamedMap/LabelTable');
      for k=1:length(uid_LabelTable);
        labeltable = branch(namedmap, uid_LabelTable(k));
        uid_Label = find(labeltable, '/LabelTable/Label');
        for l=1:length(uid_Label)
          % there are also potentially intersting atributes here, but I don't know what to do with them
          NamedMap(j).LabelTable.Label{l} = get(labeltable, children(labeltable, uid_Label(l)), 'value');
          attr = attributes(branch(labeltable, uid_Label(l)), 'get', 1);
          if isstruct(attr), attr = {attr}; end % treat one attribute just like multiple attributes
          for m=1:numel(attr)
            switch attr{m}.key
              case 'Key'
                NamedMap(j).LabelTable.Key(l)   = str2double(attr{m}.val);
              case 'Red'
                NamedMap(j).LabelTable.Red(l)   = str2double(attr{m}.val);
              case 'Green'
                NamedMap(j).LabelTable.Green(l) = str2double(attr{m}.val);
              case 'Blue'
                NamedMap(j).LabelTable.Blue(l)  = str2double(attr{m}.val);
              case 'Alpha'
                NamedMap(j).LabelTable.Alpha(l) = str2double(attr{m}.val);
            end
          end
        end
      end % for each LabelTable
    end % for each NamedMap
  end % if NamedMap

  uid_BrainModel = find(map, '/MatrixIndicesMap/BrainModel');
  for j=1:length(uid_BrainModel)
    brainmodel = branch(map, uid_BrainModel(j));

    % get the attributes of each model
    attr = attributes(brainmodel, 'get', 1);
    if isstruct(attr), attr = {attr}; end % treat one attribute just like multiple attributes
    for k=1:numel(attr)
      if any(strcmp(attr{k}.key, numericAttributeTypes))
        MatrixIndicesMap(i).BrainModel(j).(attr{k}.key) = str2num(attr{k}.val);
      else
        MatrixIndicesMap(i).BrainModel(j).(attr{k}.key) = attr{k}.val;
      end
    end % for each attr

    switch MatrixIndicesMap(i).BrainModel(j).ModelType
      case 'CIFTI_MODEL_TYPE_SURFACE'
        switch Cifti.Version
          case {'1' '1.0'}
            uid = find(brainmodel, '/BrainModel/NodeIndices');
            % use VertexIndices rather than NodeIndices
            try, MatrixIndicesMap(i).BrainModel(j).VertexIndices = str2num(get(brainmodel, children(brainmodel, uid), 'value')) + 1; end
            % also copy the global surface information to a higher level
            Surface(end+1).BrainStructure          = MatrixIndicesMap(i).BrainModel(j).BrainStructure;
            Surface(end  ).SurfaceNumberOfVertices = MatrixIndicesMap(i).BrainModel(j).SurfaceNumberOfNodes;

          case {'2' '2.0'}
            uid = find(brainmodel, '/BrainModel/VertexIndices');
            try, MatrixIndicesMap(i).BrainModel(j).VertexIndices = str2num(get(brainmodel, children(brainmodel, uid), 'value')) + 1; end
            % also copy the global surface information to a higher level
            Surface(end+1).BrainStructure          = MatrixIndicesMap(i).BrainModel(j).BrainStructure;
            Surface(end  ).SurfaceNumberOfVertices = MatrixIndicesMap(i).BrainModel(j).SurfaceNumberOfVertices;

          otherwise
            ft_error('unsupported cifti version');
        end % switch version

      case 'CIFTI_MODEL_TYPE_VOXELS'
        tmp = str2num(get(brainmodel, children(brainmodel, find(brainmodel, '/BrainModel/VoxelIndicesIJK')), 'value'));
        MatrixIndicesMap(i).BrainModel(j).VoxelIndicesIJK = reshape(tmp, 3, [])' + 1; % transpose, one offset

      otherwise
        ft_error('unsupported ModelType');
    end % switch ModelType
  end % for each BrainModel
end % for each MatrixIndicesMap

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read the voxel data section
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if readdata
  fseek(fid, hdr.vox_offset, 'bof');
  switch hdr.datatype
    case   2, [voxdata, nitemsread] = fread(fid, inf, 'uchar');   assert(nitemsread==prod(hdr.dim(2:end)), 'could not read all data');
    case   4, [voxdata, nitemsread] = fread(fid, inf, 'short');   assert(nitemsread==prod(hdr.dim(2:end)), 'could not read all data');
    case   8, [voxdata, nitemsread] = fread(fid, inf, 'int');     assert(nitemsread==prod(hdr.dim(2:end)), 'could not read all data');
    case  16, [voxdata, nitemsread] = fread(fid, inf, 'float');   assert(nitemsread==prod(hdr.dim(2:end)), 'could not read all data');
    case  64, [voxdata, nitemsread] = fread(fid, inf, 'double');  assert(nitemsread==prod(hdr.dim(2:end)), 'could not read all data');
    case 512, [voxdata, nitemsread] = fread(fid, inf, 'ushort');  assert(nitemsread==prod(hdr.dim(2:end)), 'could not read all data');
    case 768, [voxdata, nitemsread] = fread(fid, inf, 'uint');    assert(nitemsread==prod(hdr.dim(2:end)), 'could not read all data');
    otherwise, ft_error('unsupported datatype');
  end
  % hdr.dim(1) is the number of dimensions
  % hdr.dim(2) is reserved for the x-dimension
  % hdr.dim(3) is reserved for the y-dimension
  % hdr.dim(4) is reserved for the z-dimension
  % hdr.dim(5) is reserved for the time-dimension
  % hdr.dim(6:8) are used for CIFTI
  voxdata = reshape(voxdata, hdr.dim(6:end));
end
fclose(fid);

if inflated
  % compressed file has been unzipped on the fly, clean up
  delete(filename);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert to FieldTrip source representation, i.e. according to FT_DATATYPE_SOURCE and FT_DATATYPE_PARCELLATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dimord = cell(size(MatrixIndicesMap));

% MatrixIndicesMap.IndicesMapToDataType can be
% CIFTI_INDEX_TYPE_BRAIN_MODELS The dimension represents one or more brain models.
% CIFTI_INDEX_TYPE_PARCELS      The dimension represents a parcellation scheme.
% CIFTI_INDEX_TYPE_SERIES       The dimension represents a series of regular samples.
% CIFTI_INDEX_TYPE_SCALARS      The dimension represents named scalar maps.
% CIFTI_INDEX_TYPE_LABELS       The dimension represents named label maps.

for i=1:length(MatrixIndicesMap)
  switch MatrixIndicesMap(i).IndicesMapToDataType
    case 'CIFTI_INDEX_TYPE_BRAIN_MODELS'
      dimord(MatrixIndicesMap(i).AppliesToMatrixDimension+1) = {'pos'};

      % concatenate all into an array of brain models
      for j=1:length(MatrixIndicesMap(i).BrainModel)
        BrainModel = cat(1, BrainModel, MatrixIndicesMap(i).BrainModel(j));
      end

    case 'CIFTI_INDEX_TYPE_PARCELS'
      dimord(MatrixIndicesMap(i).AppliesToMatrixDimension+1) = {'chan'};

      % construct an array of brain models
      IndexOffset = 0;
      BrainModelParcelName = {};

      for j=1:numel(Parcel)

        for k=1:numel(Parcel(j).BrainStructure)
          sel = strcmp({Surface(:).BrainStructure}, Parcel(j).BrainStructure{k});
          BrainModel(end+1).IndexOffset             = IndexOffset;
          BrainModel(end  ).IndexCount              = numel(Parcel(j).Vertices{k});
          BrainModel(end  ).BrainStructure          = Parcel(j).BrainStructure{k};
          BrainModel(end  ).ModelType               = 'CIFTI_MODEL_TYPE_SURFACE';
          BrainModel(end  ).SurfaceNumberOfVertices = Surface(sel).SurfaceNumberOfVertices;
          BrainModel(end  ).VertexIndices           = Parcel(j).Vertices{k};
          BrainModel(end  ).VoxelIndicesIJK         = [];
          IndexOffset = IndexOffset + numel(Parcel(j).Vertices{k});
          BrainModelParcelName{end+1} = Parcel(j).Name;
        end % for each BrainStructure

        if ~isempty(Parcel(j).VoxelIndicesIJK)
          BrainModel(end+1).IndexOffset             = IndexOffset;
          BrainModel(end  ).IndexCount              = size(Parcel(j).VoxelIndicesIJK,1);
          BrainModel(end  ).BrainStructure          = 'CIFTI_STRUCTURE_INVALID';
          BrainModel(end  ).ModelType               = 'CIFTI_MODEL_TYPE_VOXELS';
          BrainModel(end  ).SurfaceNumberOfVertices = [];
          BrainModel(end  ).VertexIndices           = [];
          BrainModel(end  ).VoxelIndicesIJK         = Parcel(j).VoxelIndicesIJK;
          IndexOffset = IndexOffset + size(Parcel(j).VoxelIndicesIJK,1);
          BrainModelParcelName{end+1} = Parcel(j).Name;
        end

      end % for each Parcel

    case 'CIFTI_INDEX_TYPE_SERIES'
      % this only applies to cifti version 2
      switch MatrixIndicesMap(i).SeriesUnit
        case 'SECOND'
          dimord(MatrixIndicesMap(i).AppliesToMatrixDimension+1) = {'time'};
          Cifti.time = (((1:MatrixIndicesMap(i).NumberOfSeriesPoints)-1) * MatrixIndicesMap(i).SeriesStep + MatrixIndicesMap(i).SeriesStart) * 10^MatrixIndicesMap(i).SeriesExponent;
        case 'HERTZ'
          dimord(MatrixIndicesMap(i).AppliesToMatrixDimension+1) = {'freq'};
          Cifti.freq = (((1:MatrixIndicesMap(i).NumberOfSeriesPoints)-1) * MatrixIndicesMap(i).SeriesStep + MatrixIndicesMap(i).SeriesStart) * 10^MatrixIndicesMap(i).SeriesExponent;
          % case 'METER'
          % case 'RADIAN'
        otherwise
          ft_error('unsupported SeriesUnit');
      end % switch SeriesUnit

    case 'CIFTI_INDEX_TYPE_SCALARS'
      dimord{MatrixIndicesMap(i).AppliesToMatrixDimension+1} = []; % scalars are not explicitly represented
      for j=1:length(NamedMap)
        Cifti.mapname{j} = fixname(NamedMap(j).MapName);
      end

    case 'CIFTI_INDEX_TYPE_LABELS'
      dimord{MatrixIndicesMap(i).AppliesToMatrixDimension+1} = []; % labels are not explicitly represented
      for j=1:length(NamedMap)
        key = NamedMap(j).LabelTable.Key;
        lab = NamedMap(j).LabelTable.Label;
        sel = key>0;
        if isfield(NamedMap(j).LabelTable, 'Red')
          % assume rgba to be also specified
          rgba = [NamedMap(j).LabelTable.Red(:) ...
                  NamedMap(j).LabelTable.Green(:) ...
                  NamedMap(j).LabelTable.Blue(:) ...
                  NamedMap(j).LabelTable.Alpha(:)];
          Cifti.rgba{j} = rgba(sel,:);
        end
        Cifti.labeltable{j}(key(sel)) = lab(sel);
        Cifti.mapname{j} = fixname(NamedMap(j).MapName);
      end

    case 'CIFTI_INDEX_TYPE_TIME_POINTS'
      % this only applies to cifti-1, in cifti-2 this has been replaced by CIFTI_INDEX_TYPE_SERIES
      dimord(MatrixIndicesMap(i).AppliesToMatrixDimension+1) = {'time'};
      switch MatrixIndicesMap(i).TimeStepUnits
        case 'NIFTI_UNITS_SEC'
          Cifti.fsample = 1/str2double(MatrixIndicesMap(i).TimeStep);
        otherwise
          % other units should be trivial to implement
          ft_error('unsupported TimeStepUnits "%s"', MatrixIndicesMap(i).TimeStepUnits);
      end % switch TimeStepUnits

    otherwise
      ft_error('unsupported IndicesMapToDataType');
  end % switch IndicesMapToDataType
end % for each MatrixIndicesMap

dimord = dimord(~cellfun(@isempty, dimord));
source.dimord = sprintf('%s_', dimord{:});
source.dimord(end) = [];

if ~isempty(BrainModel)
  % the BrainModel is constructed both for dense and parcellated data

  dataIndex           = cell(size(BrainModel));
  greynodeIndex       = cell(size(BrainModel));
  brainstructureIndex = cell(size(BrainModel));
  surfaceIndex        = nan(size(BrainModel)); % remains nan if it maps into the volume

  geombeg = [BrainModel.IndexOffset]+1;
  geomend = geombeg + [BrainModel.IndexCount] - 1;

  for i=1:numel(BrainModel)
    dataIndex{i} = geombeg(i):geomend(i);
    switch BrainModel(i).ModelType
      case 'CIFTI_MODEL_TYPE_SURFACE'

        switch Cifti.Version
          case {'1' '1.0'}
            try
              greynodeIndex{i}     = BrainModel(i).VertexIndices;
            catch
              greynodeIndex{i}     = 1:BrainModel(i).SurfaceNumberOfNodes;
            end
            brainstructureIndex{i} = 1:BrainModel(i).SurfaceNumberOfNodes;
          case {'2', '2.0'}
            greynodeIndex{i}       = BrainModel(i).VertexIndices;
            brainstructureIndex{i} = 1:BrainModel(i).SurfaceNumberOfVertices;
          otherwise
            ft_error('unsupported cifti version');
        end % switch

        sel = find(strcmp({Surface(:).BrainStructure}, BrainModel(i).BrainStructure));
        assert(numel(sel)==1);
        surfaceIndex(i) = sel;

      case 'CIFTI_MODEL_TYPE_VOXELS'
        greynodeIndex{i}       = 1:BrainModel(i).IndexCount;
        brainstructureIndex{i} = 1:BrainModel(i).IndexCount;
        surfaceIndex(i) = nan; % does not map onto surface

      otherwise
        ft_error('unsupported ModelType "%s"', BrainModel(i).ModelType);
    end
  end % for each BrainModel

  pos      = zeros(0,3);
  posIndex = zeros(0,1);

  % concatenate all vertices of all surfaces, including vertices that do not have data
  for i=1:numel(Surface)
    pos      = cat(1, pos, nan(Surface(i).SurfaceNumberOfVertices, 3));
    posIndex = cat(1, posIndex, i*ones(Surface(i).SurfaceNumberOfVertices, 1));
  end

  % it would be possible to represent all voxels, but for efficiency we only include voxel positions with data
  if ~isempty(Volume) && any(isnan(surfaceIndex))
    tmp       = ft_warp_apply(Volume.Transform, cat(1, BrainModel(isnan(surfaceIndex)).VoxelIndicesIJK));
    pos       = cat(1, pos, tmp);
    posIndex  = cat(1, posIndex, nan(size(tmp,1),1));
  end

  % the surface vertices come before the voxels
  if ~isempty(Surface)
    voxeloffset = sum([Surface.SurfaceNumberOfVertices]);
  else
    voxeloffset = 0;
  end

  greynodeOffset = nan(size(BrainModel));
  for i=1:numel(BrainModel)
    if strcmp(BrainModel(i).ModelType, 'CIFTI_MODEL_TYPE_SURFACE')
      sel = find(strcmp({Surface(:).BrainStructure}, BrainModel(i).BrainStructure));
      greynodeOffset(i) = find(posIndex==sel, 1, 'first') - 1;
    else
      sel = strcmp({BrainModel(1:i-1).ModelType}, 'CIFTI_MODEL_TYPE_VOXELS');
      greynodeOffset(i) = voxeloffset + sum([BrainModel(sel).IndexCount]);
    end
    % shift the greynodes to become consistent with the voxel data
    greynodeIndex{i} = greynodeIndex{i} + greynodeOffset(i);
    % shift the brainstructures to become consistent with the voxel data
    brainstructureIndex{i} = brainstructureIndex{i} + greynodeOffset(i);
  end

end % if BrainModel

brainordinate.brainstructure      = zeros(size(pos,1),1);
brainordinate.brainstructurelabel = {};

for i=1:numel(BrainModel)
  indx = find(strcmp(brainordinate.brainstructurelabel, BrainModel(i).BrainStructure));
  if isempty(indx)
    indx = length(brainordinate.brainstructurelabel)+1;
  end
  brainordinate.brainstructure(brainstructureIndex{i}) = indx;
  brainordinate.brainstructurelabel{indx} = BrainModel(i).BrainStructure;
end

list1 = {
  'CIFTI_STRUCTURE_CORTEX_LEFT'
  'CIFTI_STRUCTURE_CORTEX_RIGHT'
  'CIFTI_STRUCTURE_CEREBELLUM'
  'CIFTI_STRUCTURE_ACCUMBENS_LEFT'
  'CIFTI_STRUCTURE_ACCUMBENS_RIGHT'
  'CIFTI_STRUCTURE_ALL_GREY_MATTER'
  'CIFTI_STRUCTURE_ALL_WHITE_MATTER'
  'CIFTI_STRUCTURE_AMYGDALA_LEFT'
  'CIFTI_STRUCTURE_AMYGDALA_RIGHT'
  'CIFTI_STRUCTURE_BRAIN_STEM'
  'CIFTI_STRUCTURE_CAUDATE_LEFT'
  'CIFTI_STRUCTURE_CAUDATE_RIGHT'
  'CIFTI_STRUCTURE_CEREBELLAR_WHITE_MATTER_LEFT'
  'CIFTI_STRUCTURE_CEREBELLAR_WHITE_MATTER_RIGHT'
  'CIFTI_STRUCTURE_CEREBELLUM_LEFT'
  'CIFTI_STRUCTURE_CEREBELLUM_RIGHT'
  'CIFTI_STRUCTURE_CEREBRAL_WHITE_MATTER_LEFT'
  'CIFTI_STRUCTURE_CEREBRAL_WHITE_MATTER_RIGHT'
  'CIFTI_STRUCTURE_CORTEX'
  'CIFTI_STRUCTURE_DIENCEPHALON_VENTRAL_LEFT'
  'CIFTI_STRUCTURE_DIENCEPHALON_VENTRAL_RIGHT'
  'CIFTI_STRUCTURE_HIPPOCAMPUS_LEFT'
  'CIFTI_STRUCTURE_HIPPOCAMPUS_RIGHT'
  'CIFTI_STRUCTURE_INVALID'
  'CIFTI_STRUCTURE_OTHER'
  'CIFTI_STRUCTURE_OTHER_GREY_MATTER'
  'CIFTI_STRUCTURE_OTHER_WHITE_MATTER'
  'CIFTI_STRUCTURE_PALLIDUM_LEFT'
  'CIFTI_STRUCTURE_PALLIDUM_RIGHT'
  'CIFTI_STRUCTURE_PUTAMEN_LEFT'
  'CIFTI_STRUCTURE_PUTAMEN_RIGHT'
  'CIFTI_STRUCTURE_THALAMUS_LEFT'
  'CIFTI_STRUCTURE_THALAMUS_RIGHT'
  };

list2 = {
  'CORTEX_LEFT'
  'CORTEX_RIGHT'
  'CEREBELLUM'
  'ACCUMBENS_LEFT'
  'ACCUMBENS_RIGHT'
  'ALL_GREY_MATTER'
  'ALL_WHITE_MATTER'
  'AMYGDALA_LEFT'
  'AMYGDALA_RIGHT'
  'BRAIN_STEM'
  'CAUDATE_LEFT'
  'CAUDATE_RIGHT'
  'CEREBELLAR_WHITE_MATTER_LEFT'
  'CEREBELLAR_WHITE_MATTER_RIGHT'
  'CEREBELLUM_LEFT'
  'CEREBELLUM_RIGHT'
  'CEREBRAL_WHITE_MATTER_LEFT'
  'CEREBRAL_WHITE_MATTER_RIGHT'
  'CORTEX'
  'DIENCEPHALON_VENTRAL_LEFT'
  'DIENCEPHALON_VENTRAL_RIGHT'
  'HIPPOCAMPUS_LEFT'
  'HIPPOCAMPUS_RIGHT'
  'INVALID'
  'OTHER'
  'OTHER_GREY_MATTER'
  'OTHER_WHITE_MATTER'
  'PALLIDUM_LEFT'
  'PALLIDUM_RIGHT'
  'PUTAMEN_LEFT'
  'PUTAMEN_RIGHT'
  'THALAMUS_LEFT'
  'THALAMUS_RIGHT'
  };

% replace the long name with the short name, i.e remove 'CIFTI_STRUCTURE_' where applicable
[dum, indx1, indx2] = intersect(brainordinate.brainstructurelabel, list1);
brainordinate.brainstructurelabel(indx1) = list2(indx2);

if ~isempty(Parcel)
  brainordinate.parcellation = zeros(size(pos,1),1);
  brainordinate.parcellationlabel = {};
  for i=1:numel(Parcel)
    brainordinate.parcellationlabel{i} = Parcel(i).Name;
    sel = find(strcmp(BrainModelParcelName, Parcel(i).Name));
    for j=1:numel(sel)
      brainordinate.parcellation(greynodeIndex{sel(j)}) = i; % FIXME should this be greynodeIndex or brainstructureIndex?
    end
  end
end

if readdata
  if isempty(Parcel)
    % the data is dense, make it consistent with the graynode positions
    dataIndex     = [dataIndex{:}];
    greynodeIndex = [greynodeIndex{:}];
    Ngreynodes    = size(pos,1);
  else
    % the data is defined on parcels
    Ngreynodes    = length(Parcel);
    dataIndex     = 1:Ngreynodes;
    greynodeIndex = 1:Ngreynodes;
  end

  switch source.dimord
    % the following representations are directly consistent with FieldTrip
    case {'pos' 'chan'}
      [m, n] = size(voxdata);
      if m>n
        dat = nan(Ngreynodes,n);
        dat(greynodeIndex(dataIndex),:) = voxdata;
      else
        dat = nan(Ngreynodes,m);
        dat(greynodeIndex(dataIndex),:) = transpose(voxdata);
      end
    case {'pos_time' 'chan_time'}
      Ntime = size(voxdata,2);
      dat = nan(Ngreynodes,Ntime);
      dat(greynodeIndex(dataIndex),:) = voxdata;
    case {'pos_pos' 'chan_chan'}
      dat = nan(Ngreynodes,Ngreynodes);
      dat(greynodeIndex(dataIndex),greynodeIndex(dataIndex)) = voxdata;
    case {'pos_pos_time' 'chan_chan_time'}
      Ntime = size(voxdata,3);
      dat = nan(Ngreynodes,Ngreynodes,Ntime);
      dat(greynodeIndex(dataIndex),greynodeIndex(dataIndex),:) = voxdata;
    case {'pos_pos_freq' 'chan_chan_freq'}
      Ntime = size(voxdata,3);
      dat = nan(Ngreynodes,Ngreynodes,Ntime);
      dat(greynodeIndex(dataIndex),greynodeIndex(dataIndex),:) = voxdata;

      % the following representations need to be transposed to be consistent with FieldTrip
    case 'time_pos'
      Ntime = size(voxdata,1);
      dat = nan(Ngreynodes,Ntime);
      dat(greynodeIndex(dataIndex),:) = transpose(voxdata);
      source.dimord = 'pos_time';
    case 'time_chan'
      Ntime = size(voxdata,1);
      dat = nan(Ngreynodes,Ntime);
      dat(greynodeIndex(dataIndex),:) = transpose(voxdata);
      source.dimord = 'chan_time';
    case 'time_pos_pos'
      Ntime = size(voxdata,1);
      dat = nan(Ngreynodes,Ngreynodes,Ntime);
      dat(greynodeIndex(dataIndex),greynodeIndex(dataIndex),:) = permute(voxdata, [2 3 1]);
      source.dimord = 'pos_pos_time';
    case 'time_chan_chan'
      Ntime = size(voxdata,1);
      dat = nan(Ngreynodes,Ngreynodes,Ntime);
      dat(greynodeIndex(dataIndex),greynodeIndex(dataIndex),:) = permute(voxdata, [2 3 1]);
      source.dimord = 'chan_chan_time';
      case 'freq_pos'
      Nfreq = size(voxdata,1);
      dat = nan(Ngreynodes,Nfreq);
      dat(greynodeIndex(dataIndex),:) = transpose(voxdata);
      source.dimord = 'pos_freq';
    case 'freq_chan'
      Nfreq = size(voxdata,1);
      dat = nan(Ngreynodes,Nfreq);
      dat(greynodeIndex(dataIndex),:) = transpose(voxdata);
      source.dimord = 'chan_freq';
    case 'freq_pos_pos'
      Nfreq = size(voxdata,1);
      dat = nan(Ngreynodes,Ngreynodes,Nfreq);
      dat(greynodeIndex(dataIndex),greynodeIndex(dataIndex),:) = permute(voxdata, [2 3 1]);
      source.dimord = 'pos_pos_freq';
    case 'freq_chan_chan'
      Nfreq = size(voxdata,1);
      dat = nan(Ngreynodes,Ngreynodes,Nfreq);
      dat(greynodeIndex(dataIndex),greynodeIndex(dataIndex),:) = permute(voxdata, [2 3 1]);
      source.dimord = 'chan_chan_freq';

    otherwise
      ft_error('unsupported dimord %s', source.dimord);
  end % switch

  if isfield(Cifti, 'mapname') && isfield(Cifti, 'labeltable') && strcmp(mapname, 'array')
    allthesame = true;
    for i=2:length(Cifti.labeltable)
      allthesame = allthesame && isequal(Cifti.labeltable{1}, Cifti.labeltable{i});
    end
    if allthesame
      ft_warning('using the same labels for all maps in the array');
      source.datalabel = Cifti.labeltable{1};
      Cifti = rmfield(Cifti, 'labeltable');
    else
      ft_error('multiple maps cannot be represented as array in the presence of different labeltables');
    end
  end

  if isfield(Cifti, 'mapname') && (length(Cifti.mapname)>1 || isfield(Cifti, 'labeltable'))
    switch mapname
      case 'field'
        % use distict names if there are multiple scalars or labels
        for i=1:length(Cifti.mapname)
          fieldname = Cifti.mapname{i};
          if isfield(Cifti, 'labeltable')
            if length(fieldname)>58
              % truncate it, needed to be able to append 'label' to the end
              fieldname = fieldname(1:58);
            end
            % append 'label' to the end
            source.([fieldname 'label']) = Cifti.labeltable{i}(:);
          end
          if isfield(Cifti, 'rgba')
            source.([fieldname 'rgba']) = Cifti.rgba{i};
          end
          source.(fieldname) = dat(:,i);
        end
      case 'array'
        source.mapname = {NamedMap.MapName}; % keep the original names, not the field names
        source.data    = dat;
        source.dimord  = [source.dimord '_mapname'];
      otherwise
        ft_error('incorrect specification of mapname "%s"', mapname);
    end % switch mapname
  else
    % the name of the data will be based on the filename
    source.data = dat;
  end
end % if readdata

source = copyfields(Cifti, source, {'time', 'freq'});
source.hdr = hdr;
source.unit = 'mm'; % per definition

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% try to get the geometrical information from the corresponding gifti files
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% use the filename prior to decompression
filename = origfile;

[p, f, x] = fileparts(filename);
t = tokenize(f, '.');

subject  = 'unknown';
dataname = 'unknown';
geomodel = '';

% the following assumes HCP/WorkBench/Caret file naming conventions
if length(t)==2
  subject  = t{1};
  dataname = t{2};
elseif length(t)==3
  subject  = t{1};
  dataname = t{2};
  content  = t{3};
elseif length(t)==4
  subject  = t{1};
  dataname = t{2};
  geomodel = t{3};
  content  = t{4};
elseif length(t)==5
  subject  = t{1};
  dataname = [t{2} '.' t{3}];
  geomodel = t{4};
  content  = t{5};
else
  ft_warning('cannot parse file name');
end

if readsurface
  % construct a list of possible file names for the surface geometry
  Lfilelist = {
    [subject '.L' '.midthickness'  '.' geomodel '.surf.gii']
    [subject '.L' '.pial'          '.' geomodel '.surf.gii']
    [subject '.L' '.white'         '.' geomodel '.surf.gii']
    [subject '.L' '.inflated'      '.' geomodel '.surf.gii']
    [subject '.L' '.very_inflated' '.' geomodel '.surf.gii']
    [subject '.L' '.sphere'        '.' geomodel '.surf.gii']
    [subject '.L' '.'              '.' geomodel '.surf.gii']
    [subject '.L' '.midthickness'               '.surf.gii']
    [subject '.L' '.pial'                       '.surf.gii']
    [subject '.L' '.white'                      '.surf.gii']
    [subject '.L' '.inflated'                   '.surf.gii']
    [subject '.L' '.very_inflated'              '.surf.gii']
    [subject '.L' '.sphere'                     '.surf.gii']
    [subject '.L'                               '.surf.gii']
    [subject '.CORTEX_LEFT'                     '.surf.gii']
    };

  Rfilelist = {
    [subject '.R' '.midthickness'  '.' geomodel  '.surf.gii']
    [subject '.R' '.pial'          '.' geomodel  '.surf.gii']
    [subject '.R' '.white'         '.' geomodel  '.surf.gii']
    [subject '.R' '.inflated'      '.' geomodel  '.surf.gii']
    [subject '.R' '.very_inflated' '.' geomodel  '.surf.gii']
    [subject '.R' '.sphere'        '.' geomodel  '.surf.gii']
    [subject '.R' '.'              '.' geomodel  '.surf.gii']
    [subject '.R' '.midthickness'                '.surf.gii']
    [subject '.R' '.pial'                        '.surf.gii']
    [subject '.R' '.white'                       '.surf.gii']
    [subject '.R' '.inflated'                    '.surf.gii']
    [subject '.R' '.very_inflated'               '.surf.gii']
    [subject '.R' '.sphere'                      '.surf.gii']
    [subject '.R'                                '.surf.gii']
    [subject '.CORTEX_RIGHT'                     '.surf.gii']
    };

  Bfilelist = {
    [subject '.midthickness'  '.' geomodel '.surf.gii']
    [subject '.pial'          '.' geomodel '.surf.gii']
    [subject '.white'         '.' geomodel '.surf.gii']
    [subject '.inflated'      '.' geomodel '.surf.gii']
    [subject '.very_inflated' '.' geomodel '.surf.gii']
    [subject '.sphere'        '.' geomodel '.surf.gii']
    [subject                  '.' geomodel '.surf.gii']
    [subject '.midthickness'               '.surf.gii']
    [subject '.pial'                       '.surf.gii']
    [subject '.white'                      '.surf.gii']
    [subject '.inflated'                   '.surf.gii']
    [subject '.very_inflated'              '.surf.gii']
    [subject '.sphere'                     '.surf.gii']
    [subject                               '.surf.gii']
    [subject '.CORTEX'                     '.surf.gii']
    };

  % assume that the surface files are in the same directory as the cifti file
  for i=1:numel(Lfilelist)
    Lfilelist{i} = fullfile(p, Lfilelist{i});
  end
  for i=1:numel(Rfilelist)
    Rfilelist{i} = fullfile(p, Rfilelist{i});
  end
  for i=1:numel(Bfilelist)
    Bfilelist{i} = fullfile(p, Bfilelist{i});
  end

  Lfilelist = cat(1, cortexleft, Lfilelist);
  Rfilelist = cat(1, cortexright, Rfilelist);

  tri = zeros(0,3);
  for i=1:length(Surface)

    switch Surface(i).BrainStructure
      case 'CIFTI_STRUCTURE_CORTEX_LEFT'
        for j=1:length(Lfilelist)
          if exist(Lfilelist{j}, 'file')
            fprintf('reading CORTEX_LEFT surface from %s\n', Lfilelist{j});
            mesh = ft_read_headshape(Lfilelist{j}, 'unit', 'mm'); % volume and surface should be in consistent units, gifti is defined in mm, wb_view also expects mm
            mesh.pos(:,1) = mesh.pos(:,1) - hemisphereoffset;
            pos(posIndex==i,:) = mesh.pos;
            tri = cat(1, tri, mesh.tri + find(posIndex==i, 1, 'first') - 1);
            break
          end
        end % for each Lfilelist

      case 'CIFTI_STRUCTURE_CORTEX_RIGHT'
        for j=1:length(Rfilelist)
          if exist(Rfilelist{j}, 'file')
            fprintf('reading CORTEX_RIGHT surface from %s\n', Rfilelist{j});
            mesh = ft_read_headshape(Rfilelist{j}, 'unit', 'mm'); % volume and surface should be in consistent units, gifti is defined in mm, wb_view also expects mm
            mesh.pos(:,1) = mesh.pos(:,1) + hemisphereoffset;
            pos(posIndex==i,:) = mesh.pos;
            tri = cat(1, tri, mesh.tri + find(posIndex==i, 1, 'first') - 1);
            break
          end
        end % for each Rfilelist

      otherwise
        for j=1:length(Bfilelist)
          if exist(Bfilelist{j}, 'file')
            fprintf('reading %s surface from %s\n', Surface(i).BrainStructure(17:end), Bfilelist{j});
            mesh = ft_read_headshape(Bfilelist{j}, 'unit', 'mm'); % volume and surface should be in consistent units, gifti is defined in mm, wb_view also expects mm
            pos(posIndex==i,:) = mesh.pos;
            tri = cat(1, tri, mesh.tri + find(posIndex==i, 1, 'first') - 1);
            break
          end
        end % for each Bfilelist

    end % switch BrainStructure
  end
else
  tri = [];
end % if readsurface

% add the vertex and voxel positions
brainordinate.pos = pos;
if ~isempty(tri)
  % add the surface triangulations
  brainordinate.tri = tri;
end

if ~isempty(Volume)
  brainordinate.dim        = Volume.VolumeDimensions;
  brainordinate.transform  = Volume.Transform;
end

if isempty(Parcel)
  % copy the geometrical description of the brainordinates into the main structure
  source = copyfields(brainordinate, source, fieldnames(brainordinate));
else
  % it is a parcellated source structure, i.e. represented by one channel per parcel
  % copy the geometrical description of the brainordinates into a sub-structure
  source.brainordinate = brainordinate;
  source.label = {Parcel(:).Name}';
end

% haslabeltable = false;
% if ~isempty(NamedMap)
%   % the following assumes a single NamedMap
%   if isfield(NamedMap, 'LabelTable')
%     % use the key-label combination in the label table
%     haslabeltable    = true;
%     key              = NamedMap.LabelTable.Key;
%     source.datalabel = NamedMap.LabelTable.Label(:);
%   end
% end

if readdata
  if isfield(source, 'data')
    % rename the data field
    source.(fixname(dataname)) = source.data;
    source = rmfield(source, 'data');
%     % adopt FT convention for parcel-to-label mapping
%     if haslabeltable
%       tempdata = nan+zeros(size(source.data));
%       for k = 1:numel(key)
%         tempdata(source.data==key(k)) = k;
%       end
%       source.data = tempdata;
%     end
%     source = rmfield(source, 'data');
  end

  % rename the datalabel field
  if isfield(source, 'datalabel')
    source.(fixname([dataname 'label'])) = source.datalabel;
    source = rmfield(source, 'datalabel');
  end

end % if readdata


text:
function f = sin(f)
%SIN   Sine of a CHEBFUN3.
%   SIN(F) returns the sine of a CHEBFUN3 object F.
%
% See also CHEBFUN3/COMPOSE.

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

% Check for empty:
if ( isempty(f) )
    return
end

op = @(x,y,z) sin(feval(f, x, y, z));   % Resample.
f = chebfun3(op, f.domain);             % Call constructor.

end


text:
%io_loadjmrui.m
%Jamie Near, McGill University 2014.
%
% USAGE:
% out=io_loadjmrui(filename);
%
% DESCRIPTION:
% Load a jMRUI text file into matlab structure format.
%
% INPUTS:
% filename   = filename of the jMRUI txt file.
%
% OUTPUTS:
% out = Input dataset in FID-A structure format.

function out=io_loadjmrui(filename);

%LOAD IN JMRUI .txt FILE
[RF,info]=io_readjmrui(filename);

txfrq=str2num(info.TransmitterFrequency);
sz=[str2num(info.PointsInDataset) 1];
dwelltime=str2num(info.SamplingInterval)*1e-3;
spectralwidth=1/dwelltime;
date=str2num(info.DateOfExperiment);
Bo=txfrq/42577000;  %JMRUI HEADER INCORRECTLY SAYS 3.0

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

dims.t=1;
dims.coils=0;
dims.averages=0;
dims.subSpecs=0;
dims.extras=0;

fids=RF(:,1);

specs=fftshift(ifft(fids,[],dims.t),dims.t);

f=[(-spectralwidth/2)+(spectralwidth/(2*sz(1))):spectralwidth/(sz(1)):(spectralwidth/2)-(spectralwidth/(2*sz(1)))];
ppm=-f/(Bo*42.577);
ppm=ppm+4.65;

out.fids=fids;
out.specs=specs;
out.sz=sz;
out.ppm=ppm;
out.t=t;
out.spectralwidth=spectralwidth;
out.dwelltime=dwelltime;
out.txfrq=txfrq;
out.date=date;
out.dims=dims;
out.Bo=Bo;
out.seq='';
out.te=[];
out.tr=[];
out.pointsToLeftshift=0;

%Write 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.writtentotext=0;
out.flags.downsampled=0;
out.flags.isFourSteps=0;
out.flags.subtracted=1;

author: CIC-methods

repo: FID-A

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


text:
function L = add_mex_source(L, directory, relName, matchPlatform, excludePlatform, extras)
% function L = add_mex_source(L, directory, relName, matchPlatform, excludePlatform, extras)
%
% Input + output argument L is a structure array of directory names, source file names,
% and extra arguments required for the compilation of MEX files. This function will
% create a new element of this structure and append it to L.
%
% Further inputs:
%   directory
%      target directory of the mex-file
%   relName
%      source file relative to 'directory'
%   matchPlatform
%      list of platforms this MEX file should only be compiled for.
%      use an empty matrix [] to compile for all platforms
%   excludePlatform
%      list of platforms this MEX file should NOT be compiled for.
%   extras
%      extra arguments to the MEX command, e.g. additional source files

% (C) 2010 S. Klanke

% Check if this file only needs compilation on certain platforms (including this one)
if nargin>3 && ~isempty(matchPlatform)
   ok = false;
   for k=1:numel(matchPlatform)
       if strcmp(matchPlatform{k}, computer)
		  ok = true;
		  break;
	   end
	end
	if ~ok
	   return
	end
end

% Check if this file cannot be compiled on certain platforms (including this one)
if nargin>4 && ~isempty(excludePlatform)
   ok = true;
   for k=1:numel(excludePlatform)
      if strcmp(excludePlatform{k}, computer)
         return;
      end
   end
end

L(end+1).dir   = directory;
L(end).relName = relName;
if nargin>5
  L(end).extras = extras;
end