text:
function gmsh_mesh3d_write ( gmsh_filename, m, node_num, node_x, ...
  element_order, element_num, element_node )

%*****************************************************************************80
%
%% GMSH_MESH3D_WRITE writes 3D mesh data as a Gmsh mesh file.
%
%  Discussion:
%
%    The node ordering for the 20 node element is not standard.
%
%    Assuming the vertices are A, B, C and D, Gmsh uses the following ordering:
%
%    1:    a
%    2:        b
%    3:            c
%    4:                d
%    5: (2*a  +b        )/3
%    6: (  a+2*b        )/3
%    7: (    2*b+  c    )/3
%    8: (      b+2*c    )/3
%    9: (  a    +2*c    )/3
%   10: (2*a    +  c    )/3
%   11: (2*a        +  d)/3
%   12: (  a        +2*d)/3
%   13: (     b     +2*d)/3
%   14: (   2*b     +  d)/3
%   15: (       +  c+2*d)/3
%   16: (       +2*c+  d)/3
%   17: (  a+  b+  c    )/3
%   18: (  a+  b    +  d)/3
%   19: (      b+  c+  d)/3
%   20: (  a+      c+  d)/3
%
%    Leo Rebholz used the following ordering:
%
%    1:    a
%    2:        b
%    3:            c
%    4:                d
%    5: (2*a  +b        )/3
%    6: (2*a    +  c    )/3
%    7: (  a+2*b        )/3
%    8: (  a    +2*c    )/3
%    9: (  a+  b+  c    )/3
%   10: (    2*b+  c    )/3
%   11: (      b+2*c    )/3
%   12: (2*a        +  d)/3
%   13: (   2*b     +  d)/3
%   14: (       +2*c+  d)/3
%   15: (  a+  b    +  d)/3
%   16: (      b+  c+  d)/3
%   17: (  a+      c+  d)/3
%   18: (  a        +2*d)/3
%   19: (     b     +2*d)/3
%   20: (       +  c+2*d)/3
%
%    Since the only 20 node data we have is from Leo, we will assume that
%    all 20 node input data is in Leo's format, and needs to be converted
%    to the Gmsh convention.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 October 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Christophe Geuzaine, Jean-Francois Remacle,
%    Gmsh: a three-dimensional finite element mesh generator with
%    built-in pre- and post-processing facilities,
%    International Journal for Numerical Methods in Engineering,
%    Volume 79, Number 11, pages 1309-1331, 2009.
%
%  Parameters:
%
%    Input, string GMSH_FILENAME, the name of the Gmsh file.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_X(M,NODE_NUM), the node coordinates.
%
%    Input, integer ELEMENT_ORDER, the order of the elements.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM), the nodes
%    that make up each element.
%
  leo_to_gmsh = [...
     1,  2,  3,  4,  5, ...
     7, 10, 11,  8,  6, ...
    12, 18, 19, 13, 20, ...
    14,  9, 15, 16, 17 ];
%
%  Enforce 1-based indexing.
%
  element_node = mesh_base_one ( node_num, element_order, element_num, ...
    element_node );
%
%  Open the file.
%
  gmsh = fopen ( gmsh_filename, 'wt' );

  if ( gmsh < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GMSH_MESH3D_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'GMSH_MESH3D_WRITE - Error!' );
  end
%
%  Write the data.
%
  fprintf ( gmsh, '$MeshFormat\n' );
  fprintf ( gmsh, '2.2 0 8\n' );
  fprintf ( gmsh, '$EndMeshFormat\n' );

  fprintf ( gmsh, '$Nodes\n' );
  fprintf ( gmsh, '%d\n', node_num );
  for node = 1 : node_num
    fprintf ( gmsh, '%d', node );
    for dim = 1 : 3
      if ( dim <= m )
        fprintf ( gmsh, '  %g', node_x(dim,node) );
      else
        fprintf ( gmsh, '  %g', 0.0 );
      end
    end
    fprintf ( gmsh, '\n' );
  end
  fprintf ( gmsh, '$EndNodes\n' );
%
%  These are the Gmsh codes for 4, 10 and 20 node tetrahedral elements.
%
  if ( element_order == 4 )
    element_type = 4;
  elseif ( element_order == 10 )
    element_type = 11;
  elseif ( element_order == 20 )
    element_type = 29;
  end

  tag_num = 2;
  tag1 = 0;
  fprintf ( gmsh, '$Elements\n' );
  fprintf ( gmsh, '%d\n', element_num );
  for element = 1 : element_num
    fprintf ( gmsh, '%d  %d  %d  %d  %d', ...
      element, element_type, tag_num, tag1, element );
    for vertex = 1 : element_order
      if ( element_order == 20 )
        v = leo_to_gmsh(vertex);
      else
        v = vertex;
      end
      fprintf ( gmsh, '  %d', element_node(v,element) );
    end
    fprintf ( gmsh, '\n' );
  end
  fprintf ( gmsh, '$EndElements\n' );

  fclose ( gmsh );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% EX_LAPLACE_1D: solve the Poisson problem in the unit interval with a B-spline discretization.

% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = nrbline ([0 0], [1 0]);

% Type of boundary conditions for each side of the domain
problem_data.nmnn_sides   = [];
problem_data.drchlt_sides = [1 2];

% Physical parameters
problem_data.c_diff  = @(x) ones(size(x));

% Source and boundary terms
n = 10;
problem_data.f = @(x) n^2 * pi^2 * sin (n*pi*x);
problem_data.g = @(x, ind) n * pi * cos (n*pi*x) * (2*ind-3); % The last part fixes the sign
problem_data.h = @(x, ind) sin (n*pi*x);

% Exact solution (optional)
problem_data.uex     = @(x) sin (n*pi*x);
problem_data.graduex = @(x) n * pi * cos (n*pi*x);

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.degree     = 3;     % Degree of the splines
method_data.regularity = 2;     % Regularity of the splines
method_data.nsub       = 20;    % Number of subdivisions
method_data.nquad      = 4;     % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
[geometry, msh, space, u] = solve_laplace (problem_data, method_data);

% 4) POST-PROCESSING
% PLOT IN MATLAB. COMPARISON WITH THE EXACT SOLUTION

vtk_pts = {linspace(0, 1, 150)};
[eu, F] = sp_eval (u, space, geometry, vtk_pts);
plot (F, eu, F, problem_data.uex(F))
legend ('Numerical solution', 'Exact solution'); axis tight

% Display errors of the computed solution in the L2 and H1 norm
[error_h1, error_l2] = ...
           sp_h1_error (space, msh, u, problem_data.uex, problem_data.graduex)

%!demo
%! ex_laplace_1d

%!test
%! problem_data.geo_name = nrbline ([0 0], [1 0]);
%! problem_data.nmnn_sides   = [];
%! problem_data.drchlt_sides = [1 2];
%! problem_data.c_diff  = @(x) ones(size(x));
%! n = 10;
%! problem_data.f = @(x) n^2 * pi^2 * sin (n*pi*x);
%! problem_data.g = @(x, ind) n * pi * cos (n*pi*x) * (2*ind-3); % The last part fixes the sign
%! problem_data.h = @(x, ind) sin (n*pi*x);
%! problem_data.uex     = @(x) sin (n*pi*x);
%! problem_data.graduex = @(x) n * pi * cos (n*pi*x);
%! method_data.degree     = 3;     % Degree of the splines
%! method_data.regularity = 2;     % Regularity of the splines
%! method_data.nsub       = 20;    % Number of subdivisions
%! method_data.nquad      = 4;     % Points for the Gaussian quadrature rule
%! [geometry, msh, space, u] = solve_laplace (problem_data, method_data);
%! [error_h1, error_l2] = sp_h1_error (space, msh, u, problem_data.uex, problem_data.graduex);
%! assert (msh.nel, 20)
%! assert (space.ndof, 23)
%! assert (error_h1, 0.824973349002851, 1e-14)
%! assert (error_l2, 0.00870630412762487, 1e-15)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function Distance = ComputeDistance( obj, pos )

meta = obj.Metadata{1};

%compute line distance and heading in the ground plane
StartPos = obj.axescoords2native(pos(1,:));
StopPos = obj.axescoords2native(pos(2,:));

DeltaX = (StartPos(1)-StopPos(1))*meta.Grid.Col.SS/cosd(meta.SCPCOA.TwistAng);
DeltaY = (StartPos(2)-StopPos(2))*meta.Grid.Row.SS/cosd(meta.SCPCOA.GrazeAng);

Distance = sqrt(DeltaX*DeltaX+DeltaY*DeltaY);

%alternate method...get lat/lon for each point and then compute the
%distance between the points
lla1 = point_slant_to_ground([StartPos(2), StartPos(1)]',meta);
lla2 = point_slant_to_ground([StopPos(2), StopPos(1)]',meta);

[d1km d2km]=lldistkm(lla1(1:2),lla2(1:2));

Distance = d1km*1000;

end

author: ngageoint


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%
% Transformations of the time variable

t=-1:.1:3;
 x=t.*exp(-t);
 plot(t,x);
 legend('x(t)')

 figure
 plot(-t,x);
 legend('x(-t)')

 figure
 a=2;
 plot((1/a)*t,x)
 legend('x(2t)') ;

 figure
  a=1/2;
 plot((1/a)*t,x)
 legend('x(1/2 t)') ;

 figure
 t0=2;
 plot(t+t0,x)
 legend('x(t-2)') ;

 figure
  t0=-3;
 plot(t+t0,x);
 legend('x(t+3)') ;

 figure
 plot(t-1,x)
 legend('x(t+ 1)')

 figure
  plot(0.5*(t-1),x)
 legend('x(2t-1)')

 figure
 plot(-0.5*(t-1),x)
 legend('x(1-2t)')

author: Sable

repo: mcbench-benchmarks


text:
function value = mean(sF, varargin)
% calculates the mean value for an univariate S2Fun or calculates the mean along a specified dimension fo a multimodal S2Fun
%
% Syntax
%   value = mean(sF)
%   sF = mean(sF, d)
%
% Input
%  sF - @S2FunHarmonic
%  d - dimension to take the mean value over
%
% Output
%  sF - S2FunHarmonic
%  value - double
%
% Description
%
% If sF is a 3x3 S2Fun then
% mean(sF) returns a 3x3 matrix with the mean values of each function
% mean(sF, 1) returns a 1x3 S2Fun wich contains the pointwise means values along the first dimension
%

s = size(sF);
if nargin == 1
  sF = sF.subSet(':');
  value = real(sF.fhat(1, :))/sqrt(4*pi);
  value = reshape(value, s);
else
  s = size(sF);
  value = sum(sF, varargin{1})./s(varargin{1});
end

end

author: mtex-toolbox

repo: mtex

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


text:
function pass = test_composition_operators( pref )
% Check that composition operations are working.

if ( nargin == 0 )
    pref = chebfunpref;
end

tol = 1e3*pref.cheb2Prefs.chebfun2eps;

f = diskfun(@(x,y) cos(x.*y) + sin(x.*(y)) );

% Multiplication
x = diskfun(@(x,y) x);
y = diskfun(@(x,y) y);
exact = @(x,y) (cos(x.*y) + sin(x.*y)).*sin((x-.1).*(y+.4));
g = diskfun(@(x,y) exact(x,y));
pass(1) = ( norm( g - f.*sin((x-.1).*(y+.4)) ) < tol );

% Cosine
exact = @(x,y) cos(cos(x.*y) + sin(x.*y));
g = diskfun(@(x,y) exact(x,y));
pass(2) = ( norm( g - cos(f) ) < tol );

% Cosh
exact = @(x,y) cosh( cos(x.*y) + sin(x.*y) );
g = diskfun(@(x,y) exact(x,y));
pass(3) = ( norm( g - cosh(f) ) < tol );

% Sine
exact = @(x,y)  sin(cos(x.*y) + sin(x.*y));
g = diskfun(@(x,y,z) exact(x,y));
pass(4) = ( norm( g - sin(f) ) < tol );

% Sinh
exact = @(x,y) sinh( cos(x.*y) + sin(x.*y));
g = diskfun(@(x,y) exact(x,y));
pass(5) = ( norm( g - sinh(f) ) < tol );

% Multiple operations:
f = diskfun(@(x,y)  sin(pi*x.*y));
pass(6) = (norm(f+f+f-3*f) < 100*tol);
pass(7) = (norm(f.*f-f.^2) < tol);

% Composition with a CHEBFUN with one column:
f = diskfun(@(x,y) x + y);
g = chebfun(@(t) t.^2, [-2, 2]);
h = compose(f, g);
h_true = diskfun(@(x,y) (x + y).^2);
pass(8) = ( norm(h - h_true) < tol );

% ... and a CHEBFUN with two columns:
G = chebfun(@(t) [t.^2, exp(t)], [-2, 2]);
H = compose(f, G);
H_true = diskfunv(@(x,y) (x + y).^2, @(x,y) exp(x + y));
pass(9) = ( norm(H - H_true) < tol );

% Composition with a CHEBFUN2 and a CHEBFUN2V:
f = diskfun(@(x,y) x + y);
g = chebfun2(@(x,y) x.^2 + y.^2, [-2, 2, -2, 2]);
h = compose(f, g);
h_true = diskfun(@(x,y) (x + y).^2);
pass(10) = ( norm(h - h_true) < tol );

H = compose(f, [g; g]);
H_true = [h_true; h_true];
pass(11) = ( norm(H - H_true) < tol );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function ft_write_data(filename, dat, varargin)

% FT_WRITE_DATA exports electrophysiological data such as EEG to a file.
%
% Use as
%   ft_write_data(filename, dat, ...)
%
% The specified filename can contain the filename extension. If it has no filename
% extension not, it will be added automatically.
%
% Additional options should be specified in key-value pairs and can be
%   'header'         header structure that describes the data, see FT_READ_HEADER
%   'event'          event structure that corresponds to the data, see FT_READ_EVENT
%   'chanindx'       1xN array, for selecting a subset of channels from header and data
%   'dataformat'     string, see below
%   'append'         boolean, not supported for all formats
%
% The supported dataformats for writing are
%   edf
%   gdf
%   anywave_ades
%   brainvision_eeg
%   neuralynx_ncs
%   neuralynx_sdma
%   plexon_nex
%   fcdc_matbin
%   fcdc_mysql
%   fcdc_buffer
%   flac, m4a, mp4, oga, ogg, wav (audio formats)
%   matlab
%
% For EEG data formats, the input data is assumed to be scaled in microvolt.
%
% See also FT_READ_HEADER, FT_READ_DATA, FT_READ_EVENT, FT_WRITE_EVENT

% Copyright (C) 2007-2014, 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$

global data_queue    % for fcdc_global
global header_queue  % for fcdc_global
global db_blob       % for fcdc_mysql
if isempty(db_blob)
  db_blob = 0;
end

% get the options
append        = ft_getopt(varargin, 'append', false);
nbits         = ft_getopt(varargin, 'nbits', 16); % for audio
chanindx      = ft_getopt(varargin, 'chanindx');
hdr           = ft_getopt(varargin, 'header');
evt           = ft_getopt(varargin, 'event');
dataformat    = ft_getopt(varargin, 'dataformat');

if isempty(dataformat)
  % only do the autodetection if the format was not specified
  dataformat = ft_filetype(filename);
end

if startsWith(dataformat, 'audio_')
  % support for  audio formats is implemented in a generic fashion
  dataformat = dataformat(7:end);
end

% convert 'yes' or 'no' string into boolean
append = istrue(append);

% determine the data size
[nchans, nsamples] = size(dat);

% ensure that the directory exists
isdir_or_mkdir(fileparts(filename));

switch dataformat

  case 'empty'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % just pretend that we are writing the data, this is only for debugging
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [numC, numS] = size(dat);
    fprintf(1,'Pretending to write %i samples from %i channels...\n',numS,numC);
    % Insert a small delay to make this more realitic for testing purposes
    % The time for writing to an actual location will differ and depend on
    % the amount of data
    pause(0.001);

  case 'fcdc_global'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % store it in a global variable, this is only for debugging
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isempty(hdr)
      header_queue = hdr;
    end
    if isempty(data_queue) || ~append
      data_queue = dat;
    else
      data_queue = cat(2, data_queue, dat);
    end

  case 'fcdc_buffer'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % network transparent buffer
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isempty(chanindx)
      % assume that the header corresponds to the original multichannel
      % file and that the data represents a subset of channels
      hdr.label  = hdr.label(chanindx);
      hdr.nChans = length(chanindx);
    end

    [host, port] = filetype_check_uri(filename);

    type = {
      'char'
      'uint8'
      'uint16'
      'uint32'
      'uint64'
      'int8'
      'int16'
      'int32'
      'int64'
      'single'
      'double'
      };

    wordsize = {
      1 % 'char'
      1 % 'uint8'
      2 % 'uint16'
      4 % 'uint32'
      8 % 'uint64'
      1 % 'int8'
      2 % 'int16'
      4 % 'int32'
      8 % 'int64'
      4 % 'single'
      8 % 'double'
      };

    % this should only be done the first time
    if ~append && ~isempty(hdr)
      % reformat the header into a buffer-compatible format
      packet.fsample   = hdr.Fs;
      packet.nchans    = hdr.nChans;
      packet.nsamples  = 0;
      packet.nevents   = 0;
      packet.data_type = find(strcmp(type, class(dat))) - 1; % zero-offset
      if isfield(hdr,'label') && iscell(hdr.label)
        packet.channel_names = hdr.label;
      end
      if isfield(hdr,'siemensap')
        if isa(hdr.siemensap, 'uint8')
          packet.siemensap = hdr.siemensap;
        else
          %          try
          %            packet.siemensap = matlab2sap(hdr.siemensap);
          %          catch
          warning 'Ignoring field "siemensap"';
          %          end
        end
      end
      if isfield(hdr,'nifti_1')
        if isa(hdr.nifti_1, 'uint8')
          packet.nifti_1 = hdr.nifti_1;
        else
          try
            packet.nifti_1 = encode_nifti1(hdr.nifti_1);
          catch
            warning 'Ignoring field "nifti_1"';
          end
        end
      end
      if isfield(hdr,'ctf_res4')
        if isa(hdr.ctf_res4, 'uint8')
          packet.ctf_res4 = hdr.ctf_res4;
        else
          warning 'Ignoring non-uint8 field "ctf_res4"';
        end
      end

      % try to put_hdr and initialize if necessary
      try
        % try writing the packet
        buffer('put_hdr', packet, host, port);

      catch
        if contains(lasterr, 'Buffer size N must be an integer-valued scalar double.')
          % this happens if the MATLAB75/toolbox/signal/signal/buffer
          % function is used instead of the FieldTrip buffer
          ft_error('the FieldTrip buffer mex file was not found on your path, it should be in fieldtrip/fileio/private');

        elseif contains(lasterr, 'failed to create socket') && (strcmp(host, 'localhost') || strcmp(host, '127.0.0.1'))

          % start a local instance of the TCP server
          ft_warning('starting FieldTrip buffer on %s:%d', host, port);
          buffer('tcpserver', 'init', host, port);
          pause(1);

          % rewrite the packet until success
          success = false;
          while ~success
            try
              % it may take some time before the TCP server is fully initialized
              % try writing the packet again
              buffer('put_hdr', packet, host, port);
              success = true;
            catch
              success = false;
            end
          end
        end % if strfind...

      end % try

    end % writing header

    if ~isempty(dat)
      max_nsamples = 32556;
      if size(dat,2)>max_nsamples
        % FIXME this is a hack to split large writes into multiple smaller writes
        % this is to work around a problem observed in the neuralynx proxy
        % when sampling 32 channels at 32KHz
        begsample = 1;
        while begsample<=size(dat,2)
          endsample = begsample - 1 + max_nsamples;
          endsample = min(endsample, size(dat,2));
          % if append is already one of the arguments, remove it from varargin
          indx = find(strcmp(varargin, 'append')); % find the "append" key
          if ~isempty(indx)
            indx = [indx indx+1];                  % remove the key and the value
            varargin(indx) = [];
          end
          ft_write_data(filename, dat(:,begsample:endsample), varargin{:}, 'append', false);
          begsample = endsample + 1;
        end
      else
        % FIXME this is the normal code, which will also be used recursively
        % reformat the data into a buffer-compatible format
        packet.nchans    = size(dat,1);
        packet.nsamples  = size(dat,2);
        packet.data_type = find(strcmp(type, class(dat))) - 1; % zero-offset
        packet.bufsize   = numel(dat) * wordsize{strcmp(type, class(dat))};
        packet.buf       = dat;
        buffer('put_dat', packet, host, port);
      end % if data larger than chuncksize
    end

  case {'brainvision_eeg', 'brainvision_vhdr'}
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % combination of *.eeg and *.vhdr file
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if append
      ft_error('appending data is not yet supported for this data format');
    end

    if nchans~=hdr.nChans && length(chanindx)==nchans
      % assume that the header corresponds to the original multichannel
      % file and that the data represents a subset of channels
      hdr.label  = hdr.label(chanindx);
      hdr.nChans = length(chanindx);
    end
    % the header should at least contain the following fields
    %   hdr.label
    %   hdr.nChans
    %   hdr.Fs
    write_brainvision_eeg(filename, hdr, dat, evt);

  case 'fcdc_matbin'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % multiplexed data in a *.bin file (ieee-le, 64 bit floating point values),
    % accompanied by a MATLAB V6 file containing the header
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [path, file, ext] = fileparts(filename);
    headerfile = fullfile(path, [file '.mat']);
    datafile   = fullfile(path, [file '.bin']);

    if append && exist(headerfile, 'file') && exist(datafile, 'file')
      % read the existing header and perform a sanity check
      old = load(headerfile);
      assert(old.hdr.nChans==size(dat,1));

      % update the existing header
      hdr          = old.hdr;
      hdr.nSamples = hdr.nSamples + nsamples;

      % there are no new events
      if isfield(old, 'event')
        event = old.event;
      else
        event = [];
      end

      save(headerfile, 'hdr', 'event', '-v6');

      % update the data file
      fid = fopen_or_error(datafile,'ab','ieee-le');
      fwrite(fid, dat, hdr.precision);
      fclose(fid);

    else
      hdr.nSamples = nsamples;
      hdr.nTrials  = 1;
      if nchans~=hdr.nChans && length(chanindx)==nchans
        % assume that the header corresponds to the original multichannel
        % file and that the data represents a subset of channels
        hdr.label     = hdr.label(chanindx);
        hdr.nChans    = length(chanindx);
      end
      if ~isfield(hdr, 'precision')
        hdr.precision = 'double';
      end
      % there are no events
      event = [];
      % write the header file
      save(headerfile, 'hdr', 'event', '-v6');

      % write the data file
      fid = fopen_or_error(datafile,'wb','ieee-le');
      fwrite(fid, dat, hdr.precision);
      fclose(fid);
    end

  case 'fcdc_mysql'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % write to a MySQL server listening somewhere else on the network
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % check that the required low-level toolbox is available
    ft_hastoolbox('mysql', 1);
    db_open(filename);
    if ~isempty(hdr) && isempty(dat)
      % insert the header information into the database
      if db_blob
        % insert the structure into the database table as a binary blob
        db_insert_blob('fieldtrip.header', 'msg', hdr);
      else
        % make a structure with the same elements as the fields in the database table
        s             = struct;
        s.Fs          = hdr.Fs;           % sampling frequency
        s.nChans      = hdr.nChans;       % number of channels
        s.nSamples    = hdr.nSamples;     % number of samples per trial
        s.nSamplesPre = hdr.nSamplesPre;  % number of pre-trigger samples in each trial
        s.nTrials     = hdr.nTrials;      % number of trials
        s.label       = mxSerialize(hdr.label);
        try
          s.msg = mxSerialize(hdr);
        catch
          ft_warning(lasterr);
        end
        db_insert('fieldtrip.header', s);
      end

    elseif isempty(hdr) && ~isempty(dat)
      dim = size(dat);
      if numel(dim)==2
        % ensure that the data dimensions correspond to ntrials X nchans X samples
        dim = [1 dim];
        dat = reshape(dat, dim);
      end
      ntrials = dim(1);
      for i=1:ntrials
        if db_blob
          % insert the data into the database table as a binary blob
          db_insert_blob('fieldtrip.data', 'msg', reshape(dat(i,:,:), dim(2:end)));
        else
          % create a structure with the same fields as the database table
          s = struct;
          s.nChans   = dim(2);
          s.nSamples = dim(3);
          try
            s.data = mxSerialize(reshape(dat(i,:,:), dim(2:end)));
          catch
            ft_warning(lasterr);
          end
          % insert the structure into the database
          db_insert('fieldtrip.data', s);
        end
      end

    else
      ft_error('you should specify either the header or the data when writing to a MySQL database');
    end

  case 'matlab'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % plain MATLAB file
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [path, file, ext] = fileparts(filename);
    filename = fullfile(path, [file '.mat']);
    if      append &&  exist(filename, 'file')
      % read the previous header and data from MATLAB file
      prev = load(filename);
      if ~isempty(hdr) && ~isequal(hdr, prev.hdr)
        ft_error('inconsistent header');
      else
        % append the new data to that from the MATLAB file
        dat = cat(2, prev.dat, dat);
      end
    elseif  append && ~exist(filename, 'file')
      % file does not yet exist, which is not a problem
    elseif ~append &&  exist(filename, 'file')
      ft_warning('deleting existing file ''%s''', filename);
      delete(filename);
    elseif ~append && ~exist(filename, 'file')
      % file does not yet exist, which is not a problem
    end
    save(filename, 'dat', 'hdr');

  case 'mff'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % MFF files using Phillips plugin
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ft_hastoolbox('mffmatlabio', 1);
    mff_fileio_write(filename, hdr, dat, evt);

  case 'neuralynx_sdma'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % The first version of this file format contained in the first 8 bytes the
    % channel label as string. Subsequently it contained 32 bit integer values.
    %
    % The second version of this file format starts with 8 bytes describing (as
    % a space-padded string) the data type. The channel label is contained in
    % the filename as dataset.chanlabel.bin.
    %
    % The third version of this file format starts with 7 bytes describing (as
    % a zero-padded string) the data type, followed by the 8th byte which
    % describes the downscaling for the 8 and 16 bit integer representations.
    % The downscaling itself is represented as uint8 and should be interpreted as
    % the number of bits to shift. The channel label is contained in the
    % filename as dataset.chanlabel.bin.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    statuschannel = {
      'stx'
      'pid'
      'siz'
      'tsh'
      'tsl'
      'cpu'
      'ttl'
      'x01'
      'x02'
      'x03'
      'x04'
      'x05'
      'x06'
      'x07'
      'x08'
      'x09'
      'x10'
      'crc'
      };

    dirname = filename;
    clear filename
    [path, file] = fileparts(dirname);
    for i=1:hdr.nChans
      downscale(i) = 0;
      if ~isempty(strmatch(hdr.label{i}, statuschannel))
        format{i} = 'uint32';
      else
        format{i} = 'int32';
      end
      filename{i} = fullfile(dirname, [file '.' hdr.label{i} '.bin']);
    end

    if ~isfolder(dirname)
      mkdir(dirname);
    end

    % open and write to the output files, one for each selected channel
    fid = zeros(hdr.nChans,1);
    for j=1:hdr.nChans

      if append==false
        fid(j) = fopen_or_error(filename{j}, 'wb', 'ieee-le'); % open the file
        magic = format{j};                               % this used to be the channel name
        magic((end+1):8) = 0;                            % pad with zeros
        magic(8) = downscale(j);                         % number of bits to shift
        fwrite(fid(j), magic(1:8));                      % write the 8-byte file header
      else
        fid(j) = fopen_or_error(filename{j}, 'ab', 'ieee-le');    % open the file for appending
      end % if append

      % convert the data into the correct class
      buf = dat(j,:);
      if ~strcmp(class(buf), format{j})
        switch format{j}
          case 'int16'
            buf = int16(buf);
          case 'int32'
            buf = int32(buf);
          case 'single'
            buf = single(buf);
          case 'double'
            buf = double(buf);
          case 'uint32'
            buf = uint32(buf);
          otherwise
            ft_error('unsupported format conversion');
        end
      end

      % apply the scaling, this corresponds to bit shifting
      buf = buf ./ (2^downscale(j));

      % write the segment of data to the output file
      fwrite(fid(j), buf, format{j}, 'ieee-le');

      fclose(fid(j));
    end % for each channel

  case {'flac' 'm4a' 'mp4' 'oga' 'ogg' 'wav'}
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % This writes data Y to a Windows WAVE file specified by the file name
    % WAVEFILE, with a sample rate of FS Hz and with NBITS number of bits.
    % NBITS must be 8, 16, 24, or 32.  For NBITS < 32, amplitude values
    % outside the range [-1,+1] are clipped
    %
    % Supported extensions for AUDIOWRITE are:
    %   .flac
    % 	.m4a
    % 	.mp4
    % 	.oga
    % 	.ogg
    % 	.wav
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if append
      ft_error('appending data is not supported for this data format');
    end

    if nchans~=hdr.nChans && length(chanindx)==nchans
      % assume that the header corresponds to the original multichannel
      % file and that the data represents a subset of channels
      hdr.label  = hdr.label(chanindx);
      hdr.nChans = length(chanindx);
    end

    if nchans~=1
      ft_error('this format only supports single channel continuous data');
    end

    [p, f, x] = fileparts(filename);
    if isempty(x)
      % append the format as extension
      filename = [filename '.' dataformat];
    end

    audiowrite(filename, dat, hdr.Fs, 'BitsPerSample', nbits);

  case 'plexon_nex'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % single or mulitple channel Plexon NEX file
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if append
      ft_error('appending data is not yet supported for this data format');
    end

    [path, file] = fileparts(filename);
    filename = fullfile(path, [file, '.nex']);
    if nchans~=1
      ft_error('only supported for single-channel data');
    end
    % construct a NEX structure with  the required parts of the header
    nex.hdr.VarHeader.Type       = 5; % continuous
    nex.hdr.VarHeader.Name       = hdr.label{1};
    nex.hdr.VarHeader.WFrequency = hdr.Fs;
    if isfield(hdr, 'FirstTimeStamp')
      nex.hdr.FileHeader.Frequency = hdr.Fs * hdr.TimeStampPerSample;
      nex.var.ts = hdr.FirstTimeStamp;
    else
      ft_warning('no timestamp information available');
      nex.hdr.FileHeader.Frequency  = nan;
      nex.var.ts = nan;
    end
    nex.var.indx = 0;
    nex.var.dat  = dat;

    write_plexon_nex(filename, nex);

    if 0
      % the following code snippet can be used for testing
      [nex2.var, nex2.hdr] = read_plexon_nex(filename, 'channel', 1);
    end

  case 'neuralynx_ncs'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % single channel Neuralynx NCS file
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if append
      ft_error('appending data is not yet supported for this data format');
    end

    if nchans>1
      ft_error('only supported for single-channel data');
    end

    [path, file, ext] = fileparts(filename);
    filename = fullfile(path, [file, '.ncs']);

    if nchans~=hdr.nChans && length(chanindx)==nchans
      % assume that the header corresponds to the original multichannel
      % file and that the data represents a subset of channels
      % WARNING the AD channel index assumes that the data was read from a DMA or SDMA file
      % the first 17 channels contain status info, this number is zero-offset
      ADCHANNEL  = chanindx - 17 - 1;
      LABEL      = hdr.label{chanindx};
    elseif hdr.nChans==1
      ADCHANNEL  = -1;            % unknown
      LABEL      = hdr.label{1};  % single channel
    else
      ft_error('cannot determine channel label');
    end

    FSAMPLE    = hdr.Fs;
    RECORDNSMP = 512;
    RECORDSIZE = 1044;

    % cut the downsampled LFP data into record-size pieces
    nrecords = ceil(nsamples/RECORDNSMP);
    fprintf('construct ncs with %d records\n', nrecords);

    % construct a ncs structure with all header details and the data in it
    ncs                = [];
    ncs.NumValidSamp   = ones(1,nrecords) * RECORDNSMP;   % except for the last block
    ncs.ChanNumber     = ones(1,nrecords) * ADCHANNEL;
    ncs.SampFreq       = ones(1,nrecords) * FSAMPLE;
    ncs.TimeStamp      = zeros(1,nrecords,'uint64');

    if rem(nsamples, RECORDNSMP)>0
      % the data length is not an integer number of records, pad the last record with zeros
      dat = cat(2, dat, zeros(nchans, nrecords*RECORDNSMP-nsamples));
      ncs.NumValidSamp(end) = rem(nsamples, RECORDNSMP);
    end

    ncs.dat = reshape(dat, RECORDNSMP, nrecords);

    for i=1:nrecords
      % timestamps should be 64 bit unsigned integers
      ncs.TimeStamp(i) = uint64(hdr.FirstTimeStamp) + uint64((i-1)*RECORDNSMP*hdr.TimeStampPerSample);
    end

    % add the elements that will go into the ascii header
    ncs.hdr.CheetahRev            = '4.23.0';
    ncs.hdr.NLX_Base_Class_Type   = 'CscAcqEnt';
    ncs.hdr.NLX_Base_Class_Name   = LABEL;
    ncs.hdr.RecordSize            = RECORDSIZE;
    ncs.hdr.ADChannel             = ADCHANNEL;
    ncs.hdr.SamplingFrequency     = FSAMPLE;

    % write it to a file
    fprintf('writing to %s\n', filename);
    write_neuralynx_ncs(filename, ncs);

    if 0
      % the following code snippet can be used for testing
      ncs2 = read_neuralynx_ncs(filename, 1, inf);
    end

  case 'gdf'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % multiple channel GDF file
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if append
      ft_error('appending data is not yet supported for this data format');
    end
    if ~isempty(chanindx)
      % assume that the header corresponds to the original multichannel
      % file and that the data represents a subset of channels
      hdr.label  = hdr.label(chanindx);
      hdr.nChans = length(chanindx);
    end
    write_gdf(filename, hdr, dat);

  case 'edf'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % multiple channel European Data Format file
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if append
      ft_error('appending data is not yet supported for this data format');
    end
    if ~isempty(chanindx)
      % assume that the header corresponds to the original multichannel
      % file and that the data represents a subset of channels
      hdr.label  = hdr.label(chanindx);
      hdr.nChans = length(chanindx);
    end
    write_edf(filename, hdr, dat);

  case 'anywave_ades'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % see http://meg.univ-amu.fr/wiki/AnyWave:ADES
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if append
      ft_error('appending data is not yet supported for this data format');
    end
    if ~isempty(chanindx)
      % assume that the header corresponds to the original multichannel
      % file and that the data represents a subset of channels
      hdr.label     = hdr.label(chanindx);
      hdr.chantype  = hdr.chantype(chanindx);
      hdr.chanunit  = hdr.chanunit(chanindx);
      hdr.nChans    = length(chanindx);
    end

    dattype = unique(hdr.chantype);
    datunit = cell(size(dattype));
    for i=1:numel(dattype)
      unit = hdr.chanunit(strcmp(hdr.chantype, dattype{i}));
      if ~all(strcmp(unit, unit{1}))
        ft_error('channels of the same type with different units are not supported');
      end
      datunit{i} = unit{1};
    end

    % only change these after checking channel types and units
    chantype = adestype(hdr.chantype);
    dattype  = adestype(dattype);

    % ensure that all channels have the right scaling
    for i=1:size(dat,1)
      switch chantype{i}
        case 'MEG'
          dat(i,:) = dat(i,:) * ft_scalingfactor(hdr.chanunit{i}, 'T');
        case 'Reference'
          dat(i,:) = dat(i,:) * ft_scalingfactor(hdr.chanunit{i}, 'T');
        case 'GRAD'
          dat(i,:) = dat(i,:) * ft_scalingfactor(hdr.chanunit{i}, 'T/m');
        case 'EEG'
          dat(i,:) = dat(i,:) * ft_scalingfactor(hdr.chanunit{i}, 'V');
        case 'SEEG'
          dat(i,:) = dat(i,:) * ft_scalingfactor(hdr.chanunit{i}, 'V');
        case 'EMG'
          dat(i,:) = dat(i,:) * ft_scalingfactor(hdr.chanunit{i}, 'V');
        case 'ECG'
          dat(i,:) = dat(i,:) * ft_scalingfactor(hdr.chanunit{i}, 'V');
        otherwise
          % FIXME I am not sure what scaling to apply
      end
    end

    [p, f, x] = fileparts(filename);
    filename = fullfile(p, f); % without extension
    mat2ades(dat, filename, hdr.Fs, hdr.label, chantype, dattype, datunit);

  otherwise
    ft_error('unsupported data format');
end % switch dataformat

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function type = adestype(type)
for i=1:numel(type)
  switch lower(type{i})
    case 'meggrad'
      type{i} = 'MEG'; % this is for CTF and BTi/4D
    case 'megmag'
      type{i} = 'MEG'; % this is for Neuromag and BTi/4D
    case 'megplanar'
      type{i} = 'GRAD'; % this is for Neuromag
    case {'refmag' 'refgrad'}
      type{i} = 'Reference';
    case 'eeg'
      type{i} = 'EEG';
    case {'seeg' 'ecog' 'ieeg'}
      type{i} = 'SEEG'; % all intracranial channels
    case 'ecg'
      type{i} = 'ECG';
    case 'emg'
      type{i} = 'EMG';
    case 'trigger'
      type{i} = 'Trigger';
    case 'source'
      type{i} = 'Source'; % virtual channel
    otherwise
      type{i} = 'Other';
  end
end

author: spm


text:
function y = cvx_isconstant( x, full ) %#ok
narginchk(1,2);
if nargin == 1,
    y = true;
else
    y = true( size( x ) );
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [U, S] = pca(X)
%PCA Run principal component analysis on the dataset X
%   [U, S, X] = pca(X) computes eigenvectors of the covariance matrix of X
%   Returns the eigenvectors U, the eigenvalues (on diagonal) in S
%

% Useful values
[m, n] = size(X);

% You need to return the following variables correctly.
U = zeros(n);
S = zeros(n);

% ====================== YOUR CODE HERE ======================
% Instructions: You should first compute the covariance matrix. Then, you
%               should use the "svd" function to compute the eigenvectors
%               and eigenvalues of the covariance matrix.
%
% Note: When computing the covariance matrix, remember to divide by m (the
%       number of examples).
%

Sigma = 1.0/m .* X' * X;

[U, S, V] = svd(Sigma);

% =========================================================================

end

author: atinesh-s

repo: Coursera-Machine-Learning-Stanford


text:
function q=roteu2qr(m,t)
%ROTEU2QR converts a sequence of euler angles to a real unit quaternion
% Inputs:
%
%     M(1,n)   a string of n characters from the set {'x','y','z'}
%              or, equivalently, a vector whose elements are 1, 2, or 3
%     T(1,n)   n rotation angles
%
% Outputs:
%
%     Q(1,4)   output quaternion
%
% The string M specifies the axes about which the rotations are performed.
% You cannot have the same axis in adjacent positions and so there are 12
% possibilities. Common ones are "ZXZ" and "ZYX". A positive rotation is clockwise
% if looking along the axis away from the origin.

% Suggestions:
%   (1) Should allow 1,2,3 as well as x,y,z to specify the axes

%
%      Copyright (C) Mike Brookes 2007
%      Version: $Id: roteu2qr.m,v 1.2 2007/11/21 12:42:36 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y=[2 4 1 3 1 3 2 4; 3 2 1 4 1 4 3 2; 3 4 2 1 1 2 4 3];
% m consists of a sequence of axes e.g. 'zxy'
% and t gives the rotation angles in radians
q=[1 0 0 0]';
if ischar(m)
    m=lower(m)-'w';
end
if any(abs(m-2)>1), error('Euler axis must be x,y or z'); end
for i=1:length(m)
    x=y(m(i),:);
    b=0.5*t(i);
    c=cos(b);
    s=sin(b);
    r=zeros(4,1);
    r(x(1:2))=q(x(3:4));
    r(x(5:6))=-q(x(7:8));
    q=c*q+s*r;
end
f=find(q~=0);
if (q(f(1))<0), q=-q; end

author: decouples


text:
function volROIfiles2mrGray(roiFilename, writeDir)
%
% volROIfiles2mrGray(roiFilename, [writeDir])
%
% PURPOSE:  Read ROI.mat files (mrLoadRet format) and write
% out the ROIs in a format that mrGray can read. The new files
% will have the same file name as the ROI files and will be
% written to writeDir (defaults to pwd).
%
% roiFilename can either be a single filename string, a cell
% array of file names, or a single directory name (in which
% case all .mat files in the dir will be converted).
%
%
% HISTORY:
%   2002.01.11 RFD (bob@white.stanford.edu) wrote it.
%

if(~iscell(roiFilename))
		 d = dir(roiFilename);
		 if(length(d)==1 & d.isdir==0)
		     roiFilename = {roiFilename};
		 else
		     % It's a directory- get all .mat files
         p = roiFilename;
         roiFilename = {};
         for(ii=1:length(d))
             if(~d(ii).isdir & length(d(ii).name)>3 ...
                   & strcmpi(d(ii).name(end-3:end), '.mat'))
               roiFilename{end+1} = fullfile(p,d(ii).name);
             end
         end
     end
end

for(ii=1:length(roiFilename))
    roi = load(roiFilename{ii});
    [p, n, e] = fileparts(roiFilename{ii});
    if(exist('writeDir', 'var'))
        p = writeDir;
    end
    outName = [n,'.roi'];
    disp(['Writing ',outName,' to ',p,'.']);
    volROI2mrGray(outName, roi.ROI, p);
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [dat_out,clindx,keepit] = iimg_cluster_prune(dat,datsig,volInfo)
% Prunes an image data vector (dat) assumed to contain suprathreshold
% contiguous clusters ("blobs") by saving only those blobs that have one or
% more significant (nonzero) elements in another image, datsig
%
% :Usage:
% ::
%
%     [dat_out,clindx,keepit] = iimg_cluster_prune(dat,datsig,volInfo)
%
% One intended use is to define an FWE-corrected significance map in
% datsig, and report blobs at some lower threshold (in dat) that have at
% least one corrected voxel.
%
% Try iimg_threshold to create dat vectors from Analyze images (e.g.,
% statistic images)
%
% ..
%    tor wager, july 06
% ..

% check data type and reduce to in-mask only if necessary
[dattype,dat] = iimg_check_indx(dat,volInfo,'masked');
[dattype2,datsig] = iimg_check_indx(datsig,volInfo,'masked');

% output in same format as dat input
switch dattype
    case 'full'
        dat_out = zeros(volInfo.nvox,1);          % in full image space
    case 'masked'
        dat_out = zeros(volInfo.n_inmask,1);      % in-mask space
    otherwise
        error('Internal error.  dattype should have been checked earlier.')
end

clindx = iimg_cluster_index(dat,volInfo.xyzlist');

n = max(clindx);
keepit = zeros(1,n);

for i = 1:n
    wh = clindx == i;
    keepit(i) = any(datsig(wh));

    if keepit(i)
        switch dattype
            case 'full'
                dat_out(volInfo.wh_inmask(wh)) = dat(wh);
            case 'masked'
                dat_out(wh) = dat(wh);
        end
    end
end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [modelMin, modelPruned, Ex_Rxns] = generateCompactExchModel(model,minGrowth,biomassRxn,prune,fastFVA)
% This function identifies a subnetwork with the least number of possible exchange
% reactions given the model and the applied constraints. It returns the resulting pruned model.
%
% USAGE:
%
%    [modelMin, modelPruned, Ex_Rxns] = generateCompactExchModel(model, minGrowth, biomassRxn, prune, fastFVA)
%
% INPUTS:
%    model:         model structure
%    minGrowth:     minimal Growth rate to be set on biomass reaction
%    biomassRxn:    biomass reaction name (default: 'biomass_reaction2')
%    prune:         optional: to prune the model based on exchange reactions
%                   (default: 1)
%    fastFVA:       optional: to use fastFVA instead of fluxvariability for
%                   computing FVA results (default: 0)
%    medium:        (default: {})
%
% OUTPUTS:
%    modelUpdated:  same as input model but constraints on blocked reactions
%                   are set to be 0
%    modelPruned:   pruned model, where all blocked reactions are removed
%                   (attention this seems to cause issues with GPRs)
%    Ex_Rxns:       List of exchange reactions in pruned model
%
% .. Author: - Ines Thiele, 02/2014

medium={};

if ~exist('biomassRxn','var')
    biomassRxn = 'biomass_reaction2';
end
if ~exist('prune','var')
    prune = 1;
end
if ~exist('fastFVA','var')
    fastFVA = 0;
end

% find exchange reactions
cnt=1;
for t=1:length(model.rxns)
    if  strfind(model.rxns{t}, 'EX_')
        Ex_Rxns1All(cnt,1) =model.rxns(t); %make exchange reaction list
        cnt=cnt+1;
    elseif  strfind(model.rxns{t}, 'Ex_')
        Ex_Rxns1All(cnt,1) =model.rxns(t); %make exchange reaction list
        cnt=cnt+1;
    end
end
% exclude all exchanges that have been set
Ex_Rxns2Min = Ex_Rxns1All;

model.lb(find(ismember(model.rxns,biomassRxn)))=minGrowth;% based on slowlest cell line in data

% exclude ions
Ex_Rxns2Min(ismember(Ex_Rxns2Min,medium))=[];

% identify iteratively the minimal exchange reaction network
OptExchRxns = Ex_Rxns2Min;
LO = length(OptExchRxns);
OptExchRxnsLast = OptExchRxns;
modelMin = model;
first = 1;
while LO>0
    if length(OptExchRxns)>0
        [modelMin,AddedExchange] = findMinCardModel(modelMin,OptExchRxns);
        if first == 0
            [OptExchRxns] = findOptExchRxns(modelMin,OptExchRxns);
        elseif first == 1
           [OptExchRxns] = findOptExchRxns(modelMin,AddedExchange);
            first =1;
        end
        LO = length(OptExchRxnsLast) - length(OptExchRxns)
        OptExchRxnsLast = OptExchRxns;
    else
        LO = 0;
    end
end

if prune == 1
    % prune the model. modelMin has same dimension but all blockedRxns
    % have bounds set to 0
    [modelMin,modelPruned, Ex_Rxns] = pruneModel(modelMin,minGrowth,biomassRxn);
else
    modelPruned = struct();
    Ex_Rxns = '';
end

author: opencobra

repo: cobratoolbox


text:
function [planC,allLabelNamesC,dcmExportOptS] = ...
    processAndImportAIOutput(planC,userOptS,origScanNumV,scanNumV,...
    outputScanNum,algorithm,hashID,sessionPath,cmdFlag,hWait)
%[planC,,origScanNumV,allLabelNamesC,dcmExportOptS] = ...
%processAndImportAIOutput(planC,userOptS,origScanNumV,scanNumV,...
% outputScanNum,algorithm,hashID,sessionPath,cmdFlag,hWait)
%--------------------------------------------------------------------------
% AI 08/29/22

outputS = userOptS.output;

%Loop over model outputs

outputC = fieldnames(outputS);
for nOut = 1:length(outputC)

    outType = outputC{nOut};

    switch(lower(outType))

        case 'labelmap'
            %Segmentations

            % Import segmentations
            if ishandle(hWait)
                waitbar(0.9,hWait,'Importing segmentation results to CERR');
            end
            [planC,allLabelNamesC,dcmExportOptS] = ...
                processAndImportSeg(planC,origScanNumV,scanNumV,...
                outputScanNum,sessionPath,userOptS);

        case 'dvf'
            %Deformation vector field

            outFmt = outputS.DVF.outputFormat;
            DVFpath = fullfile(sessionPath,'outputH5','DVF');
            DVFfile = dir([DVFpath,filesep,'*.h5']);
            outFile = fullfile(DVFpath,DVFfile.name);
            switch(lower(outFmt))
                case 'h5'
                    DVF4M = h5read(outFile,'/dvf');
                otherwise
                    error('Invalid model output format %s.',outFmt)
            end

            % Convert to CERR coordinate sytem

            if ~iscell(planC)
                cerrDir = planC;
                cerrDirS = dir(cerrDir);
                cerrFile = cerrDirS(3).name;
                planC = loadPlanC(fullfile(cerrDir,cerrFile),tempdir);
            else
                cerrFile = '';
            end
            indexS = planC{end};

            % Get associated scan num
            idS = userOptS.outputAssocScan.identifier;
            assocScan = getScanNumFromIdentifiers(idS,planC);

            tempOptS = userOptS;
            outTypesC = fieldnames(userOptS.output);
            matchIdx = strcmpi(outTypesC,'DVF');
            outTypesC = outTypesC(~matchIdx);
            tempOptS.output = rmfield(tempOptS.output,outTypesC);
            niiOutDir = tempOptS.output.DVF.outputDir;
            DVFfilename = strrep(DVFfile.name,'.h5','');
            dimsC = {'dx','dy','dz'};
            niiFileNameC = cell(1,length(dimsC));
            for nDim = 1:size(DVF4M,1)
                DVF3M = squeeze(DVF4M(nDim,:,:,:));
                DVF3M = permute(DVF3M,[2,3,1]);
                [DVF3M,planC] = joinH5planC(assocScan,DVF3M,[DVFfilename,'_'...
                    dimsC{nDim}],tempOptS,planC);
                niiFileNameC{nDim} = fullfile(niiOutDir,[DVFfilename,'_'...
                    dimsC{nDim},'.nii.gz']);
                fprintf('\n Writing DVF to file %s\n',niiFileNameC{nDim});
                DVF3M_nii = make_nii(DVF3M);
                save_nii(DVF3M_nii, niiFileNameC{nDim}, 0);
            end

            %Calc. deformation magnitude
            DVFmag3M = zeros(size(DVF3M));
            assocScanUID = planC{indexS.scan}(assocScan).scanUID;
            for nDim = 1:size(DVF4M,1)
                doseNum = length(planC{indexS.dose})-nDim+1;
                doseArray3M = double(getDoseArray(doseNum,planC));
                DVFmag3M = DVFmag3M + doseArray3M.^2;
            end
            DVFmag3M = sqrt(DVFmag3M);
            description = 'Deformation magnitude';
            planC = dose2CERR(DVFmag3M,[],description,'',description,...
                'CT',[],'no',assocScanUID, planC);

            % Store to deformS
            indexS = planC{end};
            idS = userOptS.register.baseScan.identifier;
            baseScanNum = getScanNumFromIdentifiers(idS,planC,1);
            idS = userOptS.register.movingScan.identifier;
            movScanNum = getScanNumFromIdentifiers(idS,planC,1);

            planC{indexS.deform}(end+1).baseScanUID = ...
                planC{indexS.scan}(baseScanNum).scanUID;
            planC{indexS.deform}(end+1).movScanUID = ...
                planC{indexS.scan}(movScanNum).scanUID;

            planC{indexS.deform}(end+1).algorithm = algorithm;
            planC{indexS.deform}(end+1).registrationTool = 'CNN';
            planC{indexS.deform}(end+1).algorithmParamsS.singContainerHash = ...
                hashID;
            planC{indexS.deform}(end+1).DVFfileName = niiFileNameC;

            if ~isempty(cerrFile)
                save_planC(planC,[],'PASSED',cerrFile);
                planC = cerrFile;
            end

        case 'derivedimage'
            %Read output image
            outFmt = outputS.derivedImage.outputFormat;
            outputImgType = outputS.derivedImage.imageType;
            imgPath = fullfile(sessionPath,'outputH5');
            imgFile = dir([imgPath,filesep,'*.h5']);
            outFile = fullfile(imgFile,DVFfile.name);
            switch(lower(outFmt))
                case 'h5'
                    img3M = h5read(outFile,'/scan'); %Dataset name for derived images?
                otherwise
                    error('Invalid model output format %s.',outFmt)
            end

            indexS = planC{end};
            initTextureS = initializeCERR('texture');
            initTextureS(end+1).textureUID = createUID('texture');
            planC{indexS.texture} = ...
                dissimilarInsert(planC{indexS.texture},initTextureS);
            currentTexture = length(planC{indexS.texture});
            planC{indexS.texture}(currentTexture).textureUID = ...
                createUID('TEXTURE');

            %Get associated scan index
            identifierS = userOptS.outputAssocScan.identifier;
            idS = rmfield(identifierS,{'warped','filtered'});
            idC = fieldnames(idS);
            if ~isempty(idC)
                assocScanNum = getScanNumFromIdentifiers(identifierS,planC);
            else
                assocScanNum = 1; %Assoc with first scan by default
            end
            assocScanUID = planC{indexS.scan}(assocScanNum).scanUID;
            planC{indexS.texture}(currentTexture).assocScanUID = assocScanUID;

            %Get associated structure index
            sizeV = size(getScanArray(assocScanNum,planC));
            minc = 1;
            maxr = sizeV(1);
            uniqueSlicesV = 1:sizeV(3);
            if isfield(userOptS.input,'structure')
                strC = {planC{indexS.structures}.structureName};
                if isfield(userOptS.input.structure,'name')
                    strName =  userOptS.input.structure.name;
                else
                    if isfield(userOptS.input.structure,'strNameToLabelMap')
                        strName =  userOptS.input.structure.strNameToLabelMap.structureName;
                    end
                end
                strIdx = getMatchingIndex(strName,strC,'EXACT');
                if ~isempty(strIdx)
                    assocStrUID = planC{indexS.structures}(strIdx).strUID;
                    planC{indexS.texture}(currentTexture).assocStructUID = assocStrUID;
                    mask3M = getStrMask(strIdx,planC);
                    [~,maxr,minc,~,~,~] = compute_boundingbox(mask3M);
                    uniqueSlicesV = find(sum(sum(mask3M))>0);
                end
            end

            % Assign parameters based on category of texture
            planC{indexS.texture}(currentTexture).parameters = ...
                userOptS;
            planC{indexS.texture}(currentTexture).description = algorithm;

            % Create Texture Scans
            [xValsV, yValsV] = ...
                getScanXYZVals(planC{indexS.scan}(assocScanNum));
            dx = median(diff(xValsV));
            dy = median(diff(yValsV));
            zV = zVals(uniqueSlicesV);
            regParamsS.horizontalGridInterval = dx;
            regParamsS.verticalGridInterval = dy;
            regParamsS.coord1OFFirstPoint = xVals(minc);
            regParamsS.coord2OFFirstPoint   = yVals(maxr);
            regParamsS.zValues  = zV;
            regParamsS.sliceThickness = ...
            [planC{indexS.scan}(assocScanNum).scanInfo(uniqueSlicesV).sliceThickness];
            assocTextureUID = planC{indexS.texture}(currentTexture).textureUID;

            %Save to planC
            planC = scan2CERR(img3M,outputImgType,'Passed',regParamsS,...
                assocTextureUID,planC);

        otherwise
            error('Invalid output type '' %s ''.',outType)

    end
    userOptS.output.(outType) =  outputS;

end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function outputImage = ifft(this, applicationDimension)
% IFFT (including IFFT-shift and normalization) using image2k
%
%   Y = MrImage()
%   Y.ifft(applicationDimension)
%
% This is a method of class MrImage.
%
% IN
%
% OUT
%
% EXAMPLE
%   Y = MrImage();
%   Y.ifft(4); % convert each voxel's time series to frequency space by
%               applying iFFT to 4th dimension
%   Y.ifft('2D') % slice-wise iFFT for transversal slices
%   Y.ifft([1 2]) % slice-wise iFFT, same as previous
%   Y.ifft([2 3]) % slice-wise iFFT, but for sagittal slices (dim 2 and 3 = 1
%                %  slice)
%   Y.ifft('3D') % volume-wise iFFT
%   Y.ifft([1 2 3]) % volume-wise iFFT, same as previous
%   Y.ifft([1 2 4]) % volume-wise iFFT, but for time series of a k-space
%                  % voxel, acquired in a 2D k-space
%
%
%   See also MrImage MrImage.image2k

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2015-12-12
% Copyright (C) 2015 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

outputImage = image2k(this, applicationDimension);

author: translationalneuromodeling

repo: tapas


text:
function [FitResult,Spectrum] = multi_comp_fit_v2(data_vol, EchoTimes, DecayMatrix, T2, Opt, varargin)
%
% multi_comp_fit_v2(data_file_name, EchoTimes, DecayMatrix, T2, Opt, ['ROI', roi_vol, 'tissue', mask_vol])
%
%**************************************************************************
% DESCRIPTION:
%
% A script to compute the multi-component T2 or T2* spectrum for an
% echo train data set using regularized (and unregulatrized) NNLS.
%
% If an ROI mask if given, the user will be prompted to proceed with a
% voxel-wise or ROI analysis.
%   voxel-wise analysis: The script will output superimposed plots of the T2/T2*
%                        spectrum for each voxel, and calculate the average
%                        MWF and gmT2/T2* for the ROI.
%   ROI analysis: The script will output plots of the average T2/T2* spectrum and
%                 the average raw data + fits for that ROI.
%
% If no mask is given then the script will output maps of the geometric
% mean T2/T2* (<T2/T2*>), the myelin water fraction (MWF) and the mono
% exponential T2/T2*.
%
% If a tissue classification mask is given (cls_file_name), the user will
% be prompted as to which tissue he/she wishes to analyze, and the maps
% will be created for that tissue only.
%
%**************************************************************************
% INPUTS:
%
%   * data_vol     = echo train data set to be analyzed
%   * EchoTimes    = vector with all echo times
%   * DecayMatrix  = matrix
%   * Opt          = struct that must contains:
%       1) RelaxationType   : 'T2' or 'T2star'
%       2) Sigma            : noise's sigma
%       3) lower_cutoff_MW  : lower cutoff on Myelin Water T2
%       4) upper_cutoff_MW  : upper cutoff on Myelin Water T2
%       5) upper_cutoff_IEW : upper cutoff on Intra/Extracellular Water T2
%   * 'ROI'        = ROI processing flag
%   * roi_vol      = ROI mask
%   * 'tissue'     = tissue processing flag
%   * mask_vol     = tissue mask, can be used to limit processing time
%   * varargin     = optional file limiting the number of echoes to be used
%                    for analysis. This can be:
%                  1) max_echoes.mnc file obtained from'in_plane_correction.m'
%                     to perform in-plane field inhomogeneity correction
%                  2) correction_mask.mnc, obtained from 'Gz_correction.m'.
%                     This binary file has voxels where the field gradient (Gz)
%                     was > 2mG/cm set to 0 (no field inhomog. correction performed*),
%                     and areas where Gz was < 2mG/cm set to 1.
%                     *The correction fails where Gz > 2mG/cm. In order to limit
%                     effect of field inhomogeneities in these voxels,
%                     processing is limited to 45 echoes, instead of 64.
%
%**************************************************************************
% EXAMPLE USES:
% multi_comp_fit_v2(EchoTimes, DecayMatrix, T2, Opt, 'tissue', cls_mask)
% --> Analyze vol T2 relaxation data, using a
% tissue classification mask: cls_mask.
%
% AUTHOR: Eva Alonso Ortiz (eva.alonso.ortiz@gmail.com)
% DATE LAST MODIFIED:
% March 2016 - WIP: complex data anaylis is not currently working, do not
% attempt to use it!
%
%**************************************************************************
% SPECIAL VERSION adapted for qMRLab
% By Ian Gagnon, 2017
% For the original version, see multi_comp_fit.m
% *************************************************************************

%%-------------------------------------------------------------------------
%% check existence of data files
%%-------------------------------------------------------------------------

% set processing labels
ROI_flag           = 0;
voxelwise_ROI_flag = 0;
tissue_flag        = 0;
mask_flag          = 0;

% default number of inputs
ndef_inputs = 5;

if nargin > ndef_inputs
    if nargin < ndef_inputs+3
        mask_opts = nargin - ndef_inputs;
    else
        mask_opts = nargin - ndef_inputs - 1;
    end

    for counter = 2:2:(mask_opts)
        switch varargin{counter-1}
            case{'ROI'}
                roi_vol = varargin{counter};
                ROI_flag = 1;
                voxelwise_ROI_flag = input('Voxel-wise ROI analysis (1) or average ROI analysis (0)?: ');
            case{'tissue'}
                mask_vol = varargin{counter};
                tissue_flag = 1;
        end
    end
end

%%-------------------------------------------------------------------------
%% Data informations
%%-------------------------------------------------------------------------

data_dim    = size(data_vol);
data_slices = data_dim(1,3);
data_height = data_dim(1,1);
data_width  = data_dim(1,2);
data_voxels = data_height*data_width;
num_echoes  = data_dim(1,4);

%%-------------------------------------------------------------------------
%% Compatibility mask/data
%%-------------------------------------------------------------------------

if ROI_flag == 1
    roi_voxels = size(roi_vol,1)*size(roi_vol,2);
    if roi_voxels ~= data_voxels
        error(sprintf('\nError in multi_comp_fit_v2: Mask file dimensions do not match data image file.\n'));
    end
end

if tissue_flag ~= 0
    mask_voxels = size(mask_vol,1)*size(mask_vol,2);
    if mask_voxels ~= data_voxels
        error(sprintf('\nError in multi_comp_fit_v2: Mask file dimensions do not match data image file.\n'));
    end
end

%%-------------------------------------------------------------------------
%% default values for NNLS fitting
%%-------------------------------------------------------------------------

T2.num = 120;

% set default values for reg-NNLS (taken from C. Chia)
if ~moxunit_util_platform_is_octave
    set(0,'RecursionLimit',5000)
end
mu = 0.25;
chi2range = [2 2.5];
chi2_min  = chi2range(1);
chi2_max  = chi2range(2);

%%-------------------------------------------------------------------------
%% apply ROI mask to data volumes
%%-------------------------------------------------------------------------

if (ROI_flag == 1 && voxelwise_ROI_flag == 0)
    mean_roi_data = get_avgROI_data(roi_vol, data_vol, data_dim);
end

%%-------------------------------------------------------------------------
%% NNLS fitting routine presets
%%-------------------------------------------------------------------------

% find cutoff indices
lower_cutoff_MW_index  = find_cutoff_index(Opt.lower_cutoff_MW, T2.vals);
upper_cutoff_MW_index  = find_cutoff_index(Opt.upper_cutoff_MW, T2.vals);
upper_cutoff_IEW_index = find_cutoff_index(Opt.upper_cutoff_IEW, T2.vals);

%%-------------------------------------------------------------------------
%% data fitting and analysis
%%-------------------------------------------------------------------------
if (ROI_flag == 1 && voxelwise_ROI_flag == 0)

    for slice = 1:data_slices

            %------------------------------
            % Do multi-exponential fitting
            %------------------------------

            % Do non-regularized NNLS
            [spectrum_NNLS(slice,:), chi2_NNLS(slice)] = do_NNLS(DecayMatrix, double(squeeze(mean_roi_data(slice,:))), Opt.Sigma(slice));

            ssq_res_NNLS(slice) = sum((DecayMatrix*squeeze(spectrum_NNLS(slice,:)') - squeeze(mean_roi_data(slice,:))').^2);
            s0_NNLS(slice)      = sum(spectrum_NNLS(slice,:));
            mwf_NNLS(slice)     = sum(spectrum_NNLS(slice,lower_cutoff_MW_index:upper_cutoff_MW_index))/s0_NNLS(slice);

            if isnan(ssq_res_NNLS(slice))
                ssq_res_NNLS(slice) = 0;
            end
            if isnan(s0_NNLS(slice))
                s0_NNLS(slice) = 0;
            end
            if isnan(mwf_NNLS(slice))
                mwf_NNLS(slice) = 0;
            end

            % Calculate the geometric mean of the T2 distribition for the non-reg NNLS
            gm_t2_NNLS(slice)     = exp(sum(squeeze(spectrum_NNLS(slice,:))'.*log(T2.vals))/sum(spectrum_NNLS(slice,:)));
            gm_IEW_t2_NNLS(slice) = exp(sum(squeeze(spectrum_NNLS(slice,upper_cutoff_MW_index:upper_cutoff_IEW_index))'.*log(T2.vals(upper_cutoff_MW_index:upper_cutoff_IEW_index)))/sum(spectrum_NNLS(slice,upper_cutoff_MW_index:upper_cutoff_IEW_index)));

            if isnan(gm_t2_NNLS(slice))
                gm_t2_NNLS(slice) = 0;
            end
            if isnan(gm_IEW_t2_NNLS(slice))
                gm_IEW_t2_NNLS(slice) = 0;
            end

            % Do regulaized NNLS
            [spectrum_regNNLS(slice,:), chi2_regNNLS(slice)] = ...
            iterate_NNLS(mu,chi2_min,chi2_max,T2.num,double(squeeze(mean_roi_data(slice,:))),DecayMatrix,chi2_NNLS(slice),Opt.Sigma(slice));

            ssq_res_regNNLS(slice) = sum((DecayMatrix*squeeze(spectrum_regNNLS(slice,:)') - squeeze(mean_roi_data(slice,:))').^2);
            s0_regNNLS(slice)      = sum(spectrum_regNNLS(slice,:));
            mwf_regNNLS(slice)     = sum(spectrum_regNNLS(slice,1:upper_cutoff_MW_index))/s0_regNNLS(slice);

            if isnan(ssq_res_regNNLS(slice))
                ssq_res_regNNLS(slice) = 0;
            end
            if isnan(s0_regNNLS(slice))
                s0_regNNLS(slice) = 0;
            end
            if isnan(mwf_regNNLS(slice))
                mwf_regNNLS(slice) = 0;
            end

            % Calculate the geometric mean of the T2 distribition for the reg NNLS
            gm_t2_regNNLS(slice)     = exp(sum(squeeze(spectrum_regNNLS(slice,:))'.*log(T2.vals))/sum(spectrum_regNNLS(slice,:)));
            gm_IEW_t2_regNNLS(slice) = exp(sum(squeeze(spectrum_regNNLS(slice,upper_cutoff_MW_index:upper_cutoff_IEW_index))'.*log(T2.vals(upper_cutoff_MW_index:upper_cutoff_IEW_index)))/sum(spectrum_regNNLS(slice,upper_cutoff_MW_index:upper_cutoff_IEW_index)));

            if isnan(gm_t2_regNNLS(slice))
                gm_t2_regNNLS(slice) = 0;
            end
            if isnan(gm_IEW_t2_regNNLS(slice))
                gm_IEW_t2_regNNLS(slice) = 0;
            end

            %-----------------------------
            % Do mono-exponential fitting
            %-----------------------------

            % Fit to single T2 decay component for comparison
            t2_guess = 100;
            guess = [t2_guess, mean_roi_data(slice,1)];

            [s0_fit(slice), t2_fit(slice)] = mono_exp_fit(double(squeeze(mean_roi_data(slice,:))), EchoTimes, double(guess));

            if isnan(s0_fit(slice))
                s0_fit(slice) = 0;
            end
            if isnan(s0_fit(slice))
                s0_fit(slice) = 0;
            end
    end

else
    if (ROI_flag == 1 && voxelwise_ROI_flag == 1)
        mask_flag = 1;
        mask_vol  = roi_vol;
    else
        if (tissue_flag == 0 && voxelwise_ROI_flag == 0)

            % create brain mask to limit processing time (and mask out noise in
            % images!)
%             create_brainmask(data_file_name);
            mask_flag = 1;
            % open brain mask
        else
            mask_flag = tissue_flag;
        end
    end

    % initialize all data vectors to zero
    spectrum_NNLS  = zeros(data_height,data_width,data_slices,T2.num);
    chi2_NNLS      = zeros(data_height,data_width,data_slices);
    gm_t2_NNLS     = zeros(data_height,data_width,data_slices);
    gm_IEW_t2_NNLS = zeros(data_height,data_width,data_slices);
    ssq_res_NNLS   = zeros(data_height,data_width,data_slices);
    s0_NNLS        = zeros(data_height,data_width,data_slices);
    mwf_NNLS       = zeros(data_height,data_width,data_slices);

    spectrum_regNNLS        = zeros(data_height,data_width,data_slices, T2.num);
    spectrum_regNNLS_A      = zeros(data_height,data_width,data_slices, T2.num);
    spectrum_regNNLS_B      = zeros(data_height,data_width,data_slices, T2.num);
    amp_spectrum_regNNLS    = zeros(data_height,data_width,data_slices, T2.num);
    phase_spectrum_regNNLS  = zeros(data_height,data_width,data_slices, T2.num);
    chi2_regNNLS            = zeros(data_height,data_width,data_slices);
    gm_t2_regNNLS           = zeros(data_height,data_width,data_slices);
    gm_IEW_t2_regNNLS       = zeros(data_height,data_width,data_slices);
    gm_MW_t2_regNNLS        = zeros(data_height,data_width,data_slices);
    ssq_res_regNNLS         = zeros(data_height,data_width,data_slices);
    s0_regNNLS              = zeros(data_height,data_width,data_slices);
    mwf_regNNLS             = zeros(data_height,data_width,data_slices);
    mw_delf_regNNLS         = zeros(data_height,data_width,data_slices);

    s0_fit = zeros(data_height,data_width,data_slices);
    t2_fit = zeros(data_height,data_width,data_slices);

    if (ROI_flag == 1 && voxelwise_ROI_flag == 1)

        %-------------------------------------------
        % Prepare figure for superimposed spectrum
        %-------------------------------------------
        figure; hold on;
        set(gca,'xscale','log','FontSize',20,'XMinorTick','on')
        xlim([T2.range(1) T2.range(2)])
        ylabel('Normalized Signal');
        if strcmp(Opt.RelaxationType,'T2')
            xlabel('T2 (ms)');
            title('T_2 Spectrum','FontSize',20);
        elseif strcmp(Opt.RelaxationType,'T2star')
            xlabel('T2* (ms)');
            title('T_2^* Spectrum','FontSize',20);
        end
        x = [Opt.upper_cutoff_MW, Opt.upper_cutoff_MW];
        y = [0, 1];
        plot(x,y,'r-','LineWidth', 4);
        grid minor;

    end

    for slice = 1:data_slices

        for i = 1:data_height
            for j = 1:data_width

                % only fill in what corresponds to the mask
                if mask_vol(i,j,slice) == mask_flag
                    this_echo_times             = EchoTimes;
                    this_num_echoes             = num_echoes;
                    this_DecayMatrix            = DecayMatrix;
                    this_lower_cutoff_MW_index  = lower_cutoff_MW_index;
                    this_upper_cutoff_MW_index  = upper_cutoff_MW_index;
                    this_upper_cutoff_IEW_index = upper_cutoff_IEW_index;

                    % Do non-regularized NNLS
                    [spectrum_NNLS(i,j,slice,:), chi2_NNLS(i,j,slice)] = do_NNLS(this_DecayMatrix, double(squeeze(data_vol(i,j,slice,1:this_num_echoes)))', Opt.Sigma(slice));

                    % Do regulaized NNLS
                    [spectrum_regNNLS(i,j,slice,:), chi2_regNNLS(i,j,slice)] = ...
                    iterate_NNLS(mu,chi2_min,chi2_max,double(squeeze(data_vol(i,j,slice,1:this_num_echoes)))',this_DecayMatrix,chi2_NNLS(i,j,slice),Opt.Sigma(slice));
                    %------------------------------
                    % Do multi-exponential fitting
                    %------------------------------

                    s0_regNNLS(i,j,slice)  = sum(spectrum_regNNLS(i,j,slice,:));
                    mwf_regNNLS(i,j,slice) = sum(spectrum_regNNLS(i,j,slice,1:this_upper_cutoff_MW_index))/s0_regNNLS(i,j,slice);

                    % Calculate the geometric mean of the T2 distribition for the reg NNLS
                    gm_t2_regNNLS(i,j,slice)     = exp(sum(squeeze(spectrum_regNNLS(i,j,slice,:)).*log(T2.vals))/sum(spectrum_regNNLS(i,j,slice,:)));
                    gm_IEW_t2_regNNLS(i,j,slice) = exp(sum(squeeze(spectrum_regNNLS(i,j,slice,this_upper_cutoff_MW_index:this_upper_cutoff_IEW_index)).*log(T2.vals(this_upper_cutoff_MW_index:this_upper_cutoff_IEW_index)))/sum(spectrum_regNNLS(i,j,slice,this_upper_cutoff_MW_index:this_upper_cutoff_IEW_index)));
                    gm_MW_t2_regNNLS(i,j,slice)  = exp(sum(squeeze(spectrum_regNNLS(i,j,slice,this_lower_cutoff_MW_index:this_upper_cutoff_MW_index)).*log(T2.vals(this_lower_cutoff_MW_index:this_upper_cutoff_MW_index)))/sum(spectrum_regNNLS(i,j,slice,this_lower_cutoff_MW_index:this_upper_cutoff_MW_index)));

                    if isnan(gm_t2_regNNLS(i,j,slice))
                        gm_t2_regNNLS(i,j,slice) = 0;
                    end
                    if isnan(gm_IEW_t2_regNNLS(i,j,slice))
                        gm_IEW_t2_regNNLS(i,j,slice) = 0;
                    end
                    if isnan(gm_MW_t2_regNNLS(i,j,slice))
                        gm_MW_t2_regNNLS(i,j,slice) = 0;
                    end
                    if isnan(s0_regNNLS(i,j,slice))
                        s0_regNNLS(i,j,slice) = 0;
                    end
                    if isnan(mwf_regNNLS(i,j,slice))
                        mwf_regNNLS(i,j,slice) = 0;
                    end

                    %-----------------------------
                    % Do mono-exponential fitting
                    %-----------------------------

                    if (ROI_flag == 1 && voxelwise_ROI_flag == 1)

                        %-----------------------------
                        %  Plot the T2 spectrum
                        %-----------------------------

                        % Normalize distributions
                        norm_spectrum_NNLS(i,j,slice,:)    = spectrum_NNLS(i,j,slice,:)/max(spectrum_NNLS(i,j,slice,:));
                        norm_spectrum_regNNLS(i,j,slice,:) = spectrum_regNNLS(i,j,slice,:)/max(spectrum_regNNLS(i,j,slice,:));

                        plot(T2.vals,squeeze(norm_spectrum_regNNLS(i,j,slice,:)), 'b-','LineWidth',0.1);
                        print('-djpeg','superimposed_spectrum.jpeg');

                    end
                end
            end
        end
    end
end

%--------------------------------------------------------------------------

if ( ROI_flag == 1 && voxelwise_ROI_flag == 0 )

    for slice = 1:data_slices

        %% Plot the T2 spectrum

        % Normalize distributions
        norm_spectrum_NNLS(slice,:)    = spectrum_NNLS(slice,:)/max(spectrum_NNLS(slice,:));
        norm_spectrum_regNNLS(slice,:) = spectrum_regNNLS(slice,:)/max(spectrum_regNNLS(slice,:));

        figure; hold on;
        grid minor;
        %bar1 = bar(T2.vals,norm_spectrum_NNLS(slice,:),'FaceColor', 'g', 'EdgeColor', 'g','LineWidth',1);
        %set(bar1,'BarWidth',0.1);

%         %--------------------
%         % If the data is simulated: Plot "true" spectrum (otherwise,
%         % comment this out)
%         %--------------------
%         T2 = [5 50]';
%         frac = [0.12 0.88]';
%
%         true_spectrum = zeros(size(T2.vals));
%         for j = 1: size(T2)
%             for i = 1:(size(T2.vals)-1)
%                 if i == 1
%                     if T2(j) == T2.vals(i)
%                         true_spectrum(i) = frac(j);
%                     end
%                 end
%
%                 if ((T2.vals(i+1) >= T2(j)) && (T2.vals(i) <= T2(j)))
%                     true_spectrum(i) = frac(j);
%                 end
%             end
%         end
%         bar2 = bar(T2.vals,true_spectrum(:),'FaceColor', 'm', 'EdgeColor', 'm','LineWidth',1);
%         %set(bar2,'BarWidth',0.1);
%         %--------------------

        plot(T2.vals,norm_spectrum_regNNLS(slice,:), 'b-','LineWidth',1.5);
        set(gca,'xscale','log','FontSize',20)
        xlim([T2.range(1) T2.range(2)])
        %xlim([0.002*1e3 4*1e3]) %keep T2* and T2 ranges the same for plot
        ylabel('Normalized Signal');
        if strcmp(Opt.RelaxationType,'T2')
            xlabel('T2 (ms)');
            title('T_2 Spectrum','FontSize',20);
        elseif strcmp(Opt.RelaxationType,'T2star')
            xlabel('T2* (ms)');
            title('T_2^* Spectrum','FontSize',20);
        end
        x = [Opt.lower_cutoff_MW, Opt.upper_cutoff_MW];
        y = [0, 1];
        plot(x,y,'r-','LineWidth', 2);
%         legend('location', 'NorthEast','NNLS multi-exp fit','True Spectrum','regularized NNLS multi-exp fit')
%         legend('location', 'NorthEast','NNLS multi-exp fit','regularized NNLS multi-exp fit')
%         legend('location', 'NorthEast','regularized NNLS multi-exp fit')
        print('-djpeg','spectrum.jpeg')

        %% Plot the decay curve and fitted line

        % Reconstruct the signal intensities and plot back on
        % top of the original data.

        single_sig_recon(slice,:)    = s0_fit(slice).*exp(-EchoTimes/t2_fit(slice));
        multi_sig_recon(slice,:)     = DecayMatrix * spectrum_NNLS(slice,:)';
        reg_multi_sig_recon(slice,:) = DecayMatrix * spectrum_regNNLS(slice,:)';

        figure; hold on;
        plot(EchoTimes, mean_roi_data(slice,:), 'kx', 'MarkerSize',10);
%         plot(echo_times, single_sig_recon(slice,:), 'g-.','LineWidth',2);
%         plot(echo_times, multi_sig_recon(slice,:), 'b-','LineWidth',2);
        plot(EchoTimes, reg_multi_sig_recon(slice,:), 'r--','LineWidth',2);
        if strcmp(Opt.RelaxationType,'T2')
%             legend('location', 'best', 'SE data','mono-exponential fit','NNLS multi-exp fit','regularized NNLS multi-exp fit');
            legend('location', 'best', 'SE data','regularized NNLS multi-exp fit');
        elseif strcmp(Opt.RelaxationType,'T2star')
%             legend('location', 'best','GRE data','mono-exponential fit','NNLS multi-exp fit','regularized NNLS multi-exp fit');
            legend('location', 'best', 'SE data','regularized NNLS multi-exp fit');
        end
        xlabel('TE (ms)');
        ylabel('Signal (arb)');
        grid on;
        print('-djpeg','signal_decay.jpeg');

        if strcmp(Opt.RelaxationType,'T2')
            %disp(fprintf('\nmulti-exp NNLS (slice %d) => <T2>: %3.2f ms, MWF:  %3.2f', slice, gm_t2_NNLS(slice), mwf_NNLS(slice)));
            disp(fprintf('multi-exp regNNLS (slice %d) => IE water <T2>: %3.2f ms, MWF:  %3.2f', slice, gm_IEW_t2_regNNLS(slice), 100*mwf_regNNLS(slice)));
            %disp(fprintf('Mono-exponential T2 fit (slice %d): %3.2f ms', slice, t2_fit(slice)));
        elseif strcmp(Opt.RelaxationType,'T2star')
            %disp(fprintf('\nmulti-exp NNLS (slice %d) => <T2*>: %3.2f ms, MWF:  %3.2f', slice, gm_t2_NNLS(slice), mwf_NNLS(slice)));
            disp(fprintf('multi-exp regNNLS (slice %d) => IE water <T2*>: %3.2f ms, MWF:  %3.2f', slice, gm_IEW_t2_regNNLS(slice), 100*mwf_regNNLS(slice)));
            %disp(fprintf('Mono-exponential T2* fit (slice %d): %3.2f ms', slice, t2_fit(slice)));
        end

   end

elseif ROI_flag == 0
    FitResult.MWF   = 100*mwf_regNNLS;
    FitResult.T2IEW = gm_IEW_t2_regNNLS;
    FitResult.T2MW  = gm_MW_t2_regNNLS;
    Spectrum        = squeeze(spectrum_regNNLS);
end

if ( tissue_flag ~= 0 || voxelwise_ROI_flag == 1 )
    mwf_regNNLS       = reshape(mwf_regNNLS,data_voxels,data_slices);
    gm_IEW_t2_regNNLS = reshape(gm_IEW_t2_regNNLS,data_voxels,data_slices);
    for slice = 1:data_slices
        mROI_mwf_regNNLS(slice)         = mean(nonzeros(mwf_regNNLS(:,slice)));
        valtmp           = std(nonzeros(mwf_regNNLS(:,slice)));
        if isempty(valtmp), stdROI_regNNLS(slice) = nan; else stdROI_regNNLS(slice) = valtmp; end
        mROI_gm_IEW_t2_regNNLS(slice)   = mean(nonzeros(gm_IEW_t2_regNNLS(:,slice)));
        valtmp = std(nonzeros(gm_IEW_t2_regNNLS(:,slice)));
        if isempty(valtmp), stdROI_gm_IEW_t2_regNNLS(slice) = nan; else stdROI_gm_IEW_t2_regNNLS(slice) = valtmp; end

%         mROI_gm_MW_t2_regNNLS(slice) = mean(nonzeros(gm_MW_t2_regNNLS(:,slice)));
%         stdROI_gm_MW_t2_regNNLS(slice) = std(nonzeros(gm_MW_t2_regNNLS(:,slice)));
%         disp(fprintf('\nAverage regMWF over tissue mask (slice %d): %3.2f +/- %3.2f %%', slice, 100*mROI_mwf_regNNLS(slice), 100*stdROI_regNNLS(slice)));
%         disp(fprintf('\nAverage IE water gm_T2 over tissue mask (slice %d): %3.2f +/- %3.2f ms', slice, mROI_gm_IEW_t2_regNNLS(slice), stdROI_gm_IEW_t2_regNNLS(slice)));
%         disp(fprintf('\nAverage Myelin water gm_T2 over tissue mask (slice %d): %3.2f +/- %3.2f ms', slice, mROI_gm_MW_t2_regNNLS(slice), stdROI_gm_MW_t2_regNNLS(slice)));
    end
end

%--------------------------------------------------------------------------

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [countV,ktransAll,veAll,kepAll,vpAll] = calcToftsFit(planC)
% Function to generate and write out DCE parameter maps, compute and save statistics to file.
% AI 5/26/16
% AI 7/06/16  Updated to provide option to scale AIF using reference region time course.
% AI 9/08/16  Removed philips scaling by SS to match Kristen's code.
% ------------------------------------------------------------------------------------------------
% INPUTS
% planC : Patient CERR file.
%         This function uses CERR structures if available.
%         Alternatively, user can select DICOM mask files.
%         Masks/ROIs are assumed to be on sequential slices (slStart to slEnd).
%
% Patient parameters are read from user-selected .txt file.
% For a list of valid input parameters, see ReadDCEParameterFile.m
% -------------------------------------------------------------------------------------------------

%% Get inputs

% Get CERR archive
if nargin==0
    [fname,fpath] = uigetfile('*.mat');
    fprintf('\n Loading patient CERR file ....\n');
    load(fullfile(fpath,fname));
    fprintf('\n Load complete.\n');
end
indexS = planC{end};
% Get parameter file input
paramDir = fullfile(getCERRPath,'PlanAnalysis','DCE-MR analysis','sample_param_and_coeff_files');
[inputS,filePath] = ReadDCEParameterFile(paramDir);
outPath = fileparts(filePath);
[paramS,shiftS] = getInputs(inputS);
% Get DCE scan set when multiple series' are available
planD = planC;
scanListC = unique({planD{indexS.scan}.scanType});
if numel(scanListC)>1
    [series,ok] = listdlg('ListString',scanListC,'Name','Select series','ListSize',[300,250],...
        'PromptString','Select DCE series','SelectionMode','Single');
    if ~ok
        return
    end
    scanNumV = strcmp({planD{indexS.scan}.scanType},scanListC{series});
    planD{indexS.scan} = planD{indexS.scan}(scanNumV);
    %Keep structures associated with selected scan
    structNumV = ismember([planD{indexS.structures}.associatedScan],find(scanNumV));
    planD{indexS.structures} =  planD{indexS.structures}(structNumV);
    %Discard empty structures
    emptyStructV = cellfun(@isempty,{planD{indexS.structures}.rasterSegments});
    planD{indexS.structures} =  planD{indexS.structures}(~emptyStructV);
end

%% Display select header info
headerS.manufacturer = planD{indexS.scan}(1).scanInfo(1).DICOMHeaders.Manufacturer;
headerS.xSize = planD{indexS.scan}(1).scanInfo(1).DICOMHeaders.Rows;
headerS.ySize = planD{indexS.scan}(1).scanInfo(1).DICOMHeaders.Columns;
headerS.TRms = planD{indexS.scan}(1).scanInfo(1).DICOMHeaders.RepetitionTime;
tDel = getTDel(planD);
headerS.tDel = double(tDel)/1000 ;                                       % convert time step to seconds from ms
headerS.B0 = planD{indexS.scan}(1).scanInfo(1).DICOMHeaders.MagneticFieldStrength;
headerS.nSlices = size(planD{indexS.scan}(1).scanArray,3);
FA = planD{indexS.scan}(1).scanInfo(1).DICOMHeaders.FlipAngle;
headerS.FA = shiftS.FAF*FA;                                              % correct flip angle by multiplicative factor FAF
if isfield(paramS,'outfileBase')
    outfileBase = paramS.outfileBase;
else
    outfileBase = '';
end

%Get no. frames
switch(headerS.manufacturer)
    case 'GE MEDICAL SYSTEMS'
        userSel = questdlg('Manufacturer: GE MEDICAL SYSTEMS','Scanner check','Continue','Cancel');
        if strcmp(userSel,'Cancel')
            return;
        end
        totalImages = planD{indexS.scan}(1).scanInfo(1).DICOMHeaders.ImagesinAcquisition;
        headerS.frames = totalImages/headerS.nSlices;
    case 'Philips Healthcare'
        headerS.frames = planD{indexS.scan}(1).scanInfo(1).DICOMHeaders.NumberofTemporalPositions;
end

fprintf(['\nVendor: %s\nNo. slices: %d\nNo. frames: %d\nTR(ms): %g\nDel_T(s): %g\n',...
    'Flip angle: %g\nImage Size(rows x cols): %d x %d\nB0: %g\n'],headerS.manufacturer,...
    headerS.nSlices,headerS.frames,headerS.TRms,headerS.tDel,headerS.FA,headerS.xSize,headerS.ySize,headerS.B0);

%% Compute model parameters
% ** Time shift for ROI time course (from input) **
ROIPointShift = shiftS.TROIPointShift * strcmp(paramS.model,'T') + ...
    shiftS.ETROIPointShift * strcmp(paramS.model,'ET');
shiftS.ROIPointShift = ROIPointShift;
% ** r1  (Using relaxivities from Pintaske, et. al. Invest Radiol 2006) **
%     Gd-DTPA, 1.5T   =  3.9 L/(mmol x s)
%     Gd-DTPA, 3.0T   =  3.3 L/(mmol x s)
%     Gd-BOPTA (multihance) 1.5T  = 8.1 L/(mmol x s)
%     Gd-BOPTA (multihance) 3.0T = 6.3 L/(mmol x s)

r1 = 3.9*(headerS.B0 < 2.0)*strcmp(paramS.cAgent,'gd')+ 3.3*(headerS.B0 > 2.0)*strcmp(paramS.cAgent,'gd') + ...
    8.1*(headerS.B0 < 2.0)*strcmp(paramS.cAgent,'mh') + 6.3*(headerS.B0 > 2.0)*strcmp(paramS.cAgent,'mh');
r1 = r1(1);

%** Vector of sampling times (1 x frames size) **
ind = double(0:(headerS.frames-1));                      % ind: vector with elements 0,1,2,3...AIF_frames-1
tsec = ind*headerS.tDel;
tmin = tsec./60;                                         % ktrans, etc. are in /min.  Parker coefs are in  min or /min.
TR = headerS.TRms/1000;                                  % TR from DICOM header is in ms. Change to s.

%**  Parameter array for the concentrarion calculation:
%    pCalcConc(1) = R1 = relaxivity of contrast agent (units are mmolxs)
%    pCalcConc(2) = TR (s)
%    pCalcConc(3) = flip angle (degrees)
%    pCalcConc(4) = pre-contrast T1 for tissue (s)
pCalcConc = [r1, TR, headerS.FA, paramS.T10];

%** Get Parker coefficient file **
[coefFile, fPath] = uigetfile([paramDir '/*.txt'],'Select coefficient file');
coefFile = fullfile(fPath,coefFile);

%% Get 3D mask
mask3M = ones(size(planD{indexS.scan}(1).scanArray));
roiSlicesV = [];
%if strcmp (paramS.maskSource,'MM') % Use CERR ROI (for Kristen)
structListC = cellfun(@lower,{planD{indexS.structures}.structureName},'un',0);
%---- For Kristen Zakian ---
%To Skip reference structures (labeled 'muscle') :
isMuscleC = strfind(structListC,'muscle');
strIdxC = cellfun(@isempty,isMuscleC,'un',0);
strNumV = find([strIdxC{:}]);
structListC = structListC(strNumV);
%--------------------------------
if numel(structListC)>1
    [strSel,ok] = listdlg('ListString',structListC,'Name','Select structure','ListSize',[300,250],...
        'PromptString','Select structure','SelectionMode','Single');
    if ~ok
        return
    end
else
    strSel = 1;
end
strNumV = strNumV(strSel);
structListC = structListC(strSel);
for l  = 1:length(structListC)
    rasterSegments = getRasterSegments(strNumV(l),planD);
    maskSlices = unique(rasterSegments(:,6));
    mask3M(:,:,maskSlices) = rasterToMask(rasterSegments,1, planD);
    roiSlicesV = [roiSlicesV;maskSlices];
end
% For Kristen (to use DICOM masks )
%else
%     % Get mask from file
%     [maskName,maskPathname] = getMaskFiles(paramS.slStart,pathBase);
%     maskFileS = dir([maskPathname,maskName,'*']);
%     numMasks = paramS.slEnd-paramS.slStart+1;
%     readMasks3M = zeros(headerS.xSize,headerS.ySize,numMasks);
%     for k = 1:numMasks
%         readMasks3M(:,:,k) = dicomread(fullfile(maskPathname,maskFileS(k).name));
%     end
%     %Get corresponding cerr file slice nos
%     [~,sliceNoV] = sort([planD{indexS.scan}(1).scanInfo.imageNumber]);
%     roiSlicesV = sliceNoV(paramS.slStart:paramS.slEnd);
%     mask3M(:,:,roiSlicesV) = readMasks3M;
%     % Copy to planD object
%     planD = maskToCERRStructure(mask3M, 0, 1, 'Imported structure', planD);
%     roiSlicesV = sort(roiSlicesV);
%end
%fprintf('ROIslices: %d',roiSlicesV);
paramS.slStart = roiSlicesV(1);
paramS.slEnd = roiSlicesV(end);
workSlices = paramS.slEnd-paramS.slStart+1;

%If using T1 maps, load T1 maps for all slices into a 3D array
if strcmp(paramS.T1Map , 'y')
    T13M = getT1Maps(headerS.xSize, headerS.ySize, headerS.nSlices, paramS.slStart, workSlices);
else
    T13M = [];
end

%% Read scans grouped by slice location & time point
sliceC = cell(1,workSlices);
infoC = cell(1,workSlices);
for k = 1:workSlices
    timeSlices = arrayfun(@(x) x.scanArray(:,:,roiSlicesV(k)),planD{indexS.scan},'un',0);
    sliceC{k} = cat(3,timeSlices{:});
    infoC{k} = planD{indexS.scan}(1).scanInfo(roiSlicesV(k)).DICOMHeaders;
end
if isfield(planD{indexS.scan}(1).scanInfo(1),'scaleSlope') & ~isempty(planD{indexS.scan}(1).scanInfo(1).scaleSlope)
    scaleSlope = planD{indexS.scan}(1).scanInfo(1).scaleSlope;
else
    scaleSlope = 1;
end
sliceTimeCourse4M = cat(4,sliceC{:});
sliceTimeCourse4M = sliceTimeCourse4M *scaleSlope; %To match Kristen's code
sliceTimeCourse4M = uint16(sliceTimeCourse4M);     %To match Kristen's code

%% Generate AIFP
fprintf('\n Generating AIF from Parker Coefficient file...');
AIFP = genAIF(tsec,coefFile,headerS.frames,tmin,paramS,shiftS);
fprintf('\n AIF generation complete.\n');

%% Get reference region (muscle) mask

useRefRegion = questdlg('Use muscle reference to scale AIF?','Reference region','Yes','No','No') ;

if strcmp(useRefRegion, 'No')
    AIFPScaled = AIFP;
else
    %if strcmp (paramS.refMaskSource,'MM') %For Kristen
    structListC = {planD{indexS.structures}.structureName};
    [strNum,ok] = listdlg('ListString',structListC,'SelectionMode','Single',...
                  'PromptString','Select reference structure');
    if ~ok
        return
    end
    muscleRasterSegS = getRasterSegments(strNum,planD);
    [muscleMaskM,muscleROISlice] = rasterToMask(muscleRasterSegS,1, planD);
    %     else
    %         %For Kristen: to read DICOM masks
    %         [refMaskName,refMaskPath] = uigetfile('.dcm','Select muscle mask file');
    %         slNum = strfind(refMaskName,'_SL');
    %         muscleROISlice = str2num(refMaskName(slNum+3:slNum+4));
    %         muscleMaskM = dicomread(fullfile(refMaskPath,refMaskName));
    %         %Copy reference region mask to CERR plan
    %         refMask3M  = zeros(size(planD{indexS.scan}(1).scanArray));
    %         muscleROISlice = sliceNoV(muscleROISlice);
    %         refMask3M(:,:,muscleROISlice) = muscleMaskM;
    %         planD = maskToCERRStructure(refMask3M, 0, 1, 'Muscle', planD);
    %     end
    % Scale AIF using reference region time course
    muscleSliceC = arrayfun(@(x) x.scanArray(:,:,muscleROISlice),planD{indexS.scan},'un',0);
    muscleSliceTimeCourse3M  = cat(3,muscleSliceC{:});
    muscleTimeCourse3M = ROIMaskedSlices(muscleSliceTimeCourse3M,muscleMaskM);
    AIFPScaled = refRegionScale(AIFP,inputS,shiftS,filePath,headerS,muscleTimeCourse3M,muscleMaskM);
end

%% Write DICOM Maps to file
getVals = 'Yes';
nRuns = 1;
while(strcmp(getVals,'Yes'))
    if nRuns~=1                %Runs with value from parameter file if checkShift='n' the first time
     shiftS.checkShift = 'y';  %Automatically set checkShift to 'y' if user chooses to run again with different shift
     [countV,ktransAll,veAll,kepAll,vpAll,dceRsqAll,UserInShift,UserInBase] = ...
     getDCEToftsParams(paramS,shiftS,headerS,roiSlicesV,sliceTimeCourse4M,pCalcConc,...
     mask3M,T13M,tmin,AIFPScaled,outFolder);
    else
     [countV,ktransAll,veAll,kepAll,vpAll,dceRsqAll,UserInShift,UserInBase,outFolder] = ...
     getDCEToftsParams(paramS,shiftS,headerS,roiSlicesV,sliceTimeCourse4M,pCalcConc,...
     mask3M,T13M,tmin,AIFPScaled,outPath);
    end
    getVals = questdlg(sprintf('Current shift = %d.\n Run again with different shift?',UserInShift),...
        'Change shift?','Yes','No','No');
    nRuns = nRuns + 1;
end

%% Store shifts to file
fid = fopen(filePath);
inputC = textscan(fid, '%s %s %n','endofline','\r\n');
fclose(fid);
fieldsC = fieldnames(paramS);
if strcmp(paramS.model,'T') %Tofts
shiftIdx = strcmp(fieldsC,'TROIPointShift');
else  %Extended tofts
shiftIdx = strcmp(fieldsC,'ETROIPointShift');
end
baseIdx = strcmp(fieldsC,'basepts');
inputC{3}(shiftIdx) = UserInShift;
inputC{3}(baseIdx) = UserInBase;
fid = fopen(filePath,'w+');
fmt = '\r\n%s\t%s\t%d';
for lineNum = 1:size(inputC{1},1)
        col1 = inputC{1}(lineNum);
        col2 = inputC{2}(lineNum);
        col3 = inputC{3}(lineNum);
        fprintf(fid,fmt,col1{1},col2{1},col3);
end
fclose(fid);

%% Write statistics to file
keepIdx = dceRsqAll> paramS.cutOff;
ktransAll(~keepIdx) = 0;
veAll(~keepIdx) = 0;
kepAll(~keepIdx) = 0;
vpAll(~keepIdx) = 0;
writeStats(outFolder,outfileBase,paramS.model,paramS.cutOff,countV,ktransAll,veAll,kepAll,vpAll);

ktransAll(ktransAll==100) = nan;
ktransAll(~keepIdx) = nan;
veAll(veAll==100) = nan;
veAll(~keepIdx) = nan;
kepAll(kepAll==100) = nan;
kepAll(~keepIdx) = nan;
vpAll(vpAll==100) = nan;
vpAll(~keepIdx) = nan;

end


text:
function [chirpF, Freqs] = bst_chirplet(sRate, nTime, chirpCenterFreqs)
% BST_CHIRPLET: Compute the Phase-Amplitude Coupling in one of several time series (directPAC)
%
% INPUTS:
%    - sRate  : Signal sampling rate (in Hz)
%    - nTime  : Number of time points of the signal to filter
%    - chirpCenterFreqs: Center frequencies of the chirplets to calculate
%
% DOCUMENTATION:
%    - The current code is inspired from Ryan Canolty's code provided originally with the article:
%         Canolty RT, Edwards E, Dalal SS, Soltani M, Nagarajan SS, Kirsch HE, Berger MS, Barbaro NM, Knight RT,
%         "High gamma power is phase-locked to theta oscillations in human neocortex",
%         Science, 2006 Sep 15;313(5793):1626-8.

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

% ===== PARSE INPUTS =====
% To avoid out-of-memeory issues: Check this on your machine, machine-specific threshold
% ***** TODO: EVALUALATE THOSE LINES *****
if (nTime > 2^23)
    nFreq = nTime;
else
    % Fixed parameter for computational ease
    nFreq = 2^ceil(log2(nTime));
end
% Raw time_support
Freqs = (sRate/nFreq) * (0:nFreq-1);
inds = Freqs > (sRate/2);
Freqs(inds) = Freqs(inds) - sRate;
% Reduce storage space
% % ***** TODO: EVALUALATE THIS LINES *****
% Freqs = single(Freqs);

% ===== CALCULATE CHIRPLETS =====
% Initialize returned matrix
chirpF = zeros(1, nFreq, length(chirpCenterFreqs));
% Make set of chirplets
fbw = 0.15;
for iif = 1:length(chirpCenterFreqs)
    % Assign or compute duration parameter
    v0 = chirpCenterFreqs(iif);  % center_frequency
    c0 = 0;  % chirp_rate
    s0 = log((2*log(2)) / (fbw^2*pi*v0^2));
    % Frequency support
    std_multiple = 6;
    vstd = sqrt((exp(-s0) + c0^2*exp(s0)) / (4*pi));
    v = Freqs; % in Hz
    iFreq = find(...
        (v0 - std_multiple * vstd <= v) & ...
        (v <= v0 + std_multiple * vstd));
    % Shorten to include only chirplet support
    v = v(iFreq);
    % Chirplet in frequency domain:
    Gk = 2^(1/4)*sqrt(-1i*c0+exp(-s0))^-1 * exp(-s0/4 + (exp(s0)*pi*(v-v0).^2)/(-1+1i*c0*exp(s0)));
    n1 = sqrt(length(Freqs)) / norm(Gk);
    % Because of discrete sampling and different time/freq sample numbers
    Gk = n1 * Gk;  % filter
    % Report in returned structures
    chirpF(1, iFreq, iif) = Gk;
end

author: brainstorm-tools