text:
%Contents of the AMD sparse matrix ordering package:
%
% amd2          - p = amd2 (A), the approximate minimum degree ordering of A
% amd_demo      - a demo of amd2, using the can_24 matrix
% amd_make      - to compile amd2 for use in MATLAB
% amd_install   - compile and install amd2 for use in MATLAB
%
% See also:  amd, amd2, colamd, symamd, colmmd, symmmd, umfpack
%
% Note that amd2 and the built-in amd function in MATLAB 7.3 and later are one
% and the same.
%
% Example:
%   p = amd2 (A) ;

% Copyright 1994-2007, Tim Davis, University of Florida,
% Patrick R. Amestoy, and Iain S. Duff.

help Contents


text:
function out=eps(x)
%eps: Floating point relative accuracy (as defined in matlab R14).

one=mp(1);%automatically inherit default precision
two=one+one;
if isempty(x)
    out=one*power(two,-one.precision*one);%automatically inherit data type
else
    out=abs(x)*power(two,-x(1).precision*one);
end

author: opencobra

repo: cobratoolbox


text:
function naca_test02 ( m, p, t )

%*****************************************************************************80
%
%% NACA_TEST02 tests NACA4_CAMBERED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real M, the maximum camber.
%    0 < M.
%
%    Input, real P, the location of maximum camber.
%    0.0 < P < 1.0.
%
%    Input, real T, the maximum relative thickness.
%    0.0 < T <= 1.0.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NACA_TEST02\n' );
  fprintf ( 1, '  NACA4_CAMBERED evaluates (xu,yu) and (xl,yl) for a NACA\n' );
  fprintf ( 1, '  cambered airfoil defined by a 4-digit code.\n' );

  if ( nargin < 1 )
    m = 0.02;
  end

  if ( nargin < 2 )
    p = 0.4;
  end

  if ( nargin < 3 )
    t = 0.12;
  end

  c = 10.0;
  n = 51;

  xc = linspace ( 0.0, c, n );
  [ xu, yu, xl, yl ] = naca4_cambered ( m, p, t, c, xc );
%
%  Plot the wing surface.
%
  plot ( xu, yu, 'b-', ...
         xl(n:-1:1), yl(n:-1:1), 'b-', ...
         'Linewidth', 3 );

  axis equal
  grid on
  xlabel ( '<---X--->', 'Fontsize', 16 );
  ylabel ( '<---Y--->', 'Fontsize', 16 );
  title ( 'NACA 4-digit cambered airfoil', 'Fontsize', 24 );

  filename = 'cambered.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Graphics saved in file "%s"\n', filename );
%
%  Save it to a file.  List the points in counter clockwise order.
%
  xy = [ xl, xu(n:-1:1); yl, yu(n:-1:1) ];
  filename = 'cambered_data.txt';
  r8mat_write ( filename, 2, 2 * n, xy );
  fprintf ( 1, '  Data saved in file "%s"\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function  test_suite = test_intersectLineCircle
%TESTINTERSECTLINECIRCLE  One-line description here, please.
%
%   output = testIntersectLineCircle(input)
%
%   Example
%   testIntersectLineCircle
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-06-06,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testIntersect(testCase) %#ok<*DEFNU>
% Should find two distinct points.

center = [10 0];
l1 = [center 0 1];
c1 = [center 5];

pts = intersectLineCircle(l1, c1);

exp = [10 -5; 10 5];
assertEqual(testCase, exp, pts, 'AbsTol', .01);

function testTangentLine(testCase)
% Should find twice the same point.

center = [10 0];
l1 = [15 0 0 1];
c1 = [center 5];

pts = intersectLineCircle(l1, c1);

exp = [15 0;15 0];
assertEqual(testCase, exp, pts, 'AbsTol', .01);

function testNoIntersect(testCase)
% Should return a 2-by-2 array full of NaN.

center = [10 0];
l1 = [16 0 0 1];
c1 = [center 5];

pts = intersectLineCircle(l1, c1);

exp = [NaN NaN;NaN NaN];
assertEqual(testCase, exp, pts);

function test_ManyLines_ManyCircles(testCase) %#ok<*DEFNU>

lines    = [ 0 0 1 0; 0 0 0 1];
circles = [ 0 0 1 ; 0 0 2];
intersectLineCircle (lines, circles);

pts = intersectLineCircle(lines, circles);

assertEqual(testCase, size(pts), [2 2 2]);

function test_ManyLines_ManyCircles_Tangents(testCase) %#ok<*DEFNU>

lines    = [ 0 0 1 0; 2 2 1 0];
circles = [ 0 0 1 ; 0 0 1];
intersectLineCircle (lines, circles);

pts = intersectLineCircle(lines, circles);

assertEqual(testCase, size(pts), [2 2]);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% [d,xi,yi,thetai] = ellipsedist_hack(xc,yc,a,b,theta,u,v)
% compute distance from points (u,v) to ellipse (xc,yc,a,b,theta)
function [d,xi,yi,thetai] = ellipsedist_hack(xc,yc,a,b,theta,u,v,npoints)

if nargin < 8,
  npoints = 100;
end

[x,y,theta] = ellipsepoints(a,b,xc,yc,theta,npoints);
[d,i] = min(dist2([x(:),y(:)],[u(:),v(:)]),[],1);
d = sqrt(d);
xi = x(i); yi = y(i); thetai = theta(i);
d = reshape(d,size(u));

author: kristinbranson

repo: JAABA


text:
function [firstCrossing, pointOfCrossing] = BF_PointOfCrossing(x,threshold)
% BF_PointOfCrossing  Linearly interpolate to the point of crossing a threshold
%
%---INPUT:
% x, a vector
% threshold, a threshold x crosses.
%
%---OUTPUTS:
% firstCrossing, the first discrete value after which a crossing event has occurred
% pointOfCrossing, the (linearly) interpolated point of crossing

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

% Find index of x at which the first crossing event occurs:
if x(1) > threshold
    firstCrossing = find((x - threshold < 0),1,'first');
else
    firstCrossing = find((x - threshold > 0),1,'first');
end

if isempty(firstCrossing)
    % Never crosses
    N = length(x);
    firstCrossing = N;
    pointOfCrossing = N;
else
    % Continuous version---the point of crossing
    valueBeforeCrossing = x(firstCrossing - 1);
    valueAfterCrossing = x(firstCrossing);
    pointOfCrossing = firstCrossing - 1 + (threshold - valueBeforeCrossing)/(valueAfterCrossing - valueBeforeCrossing);
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function d = mahalanobis(varargin)
%MAHALANOBIS Computes the Mahalanobis distance.
%   D = MAHALANOBIS(Y, X) computes the Mahalanobis distance between
%   each vector in Y to the mean (centroid) of the vectors in X, and
%   outputs the result in vector D, whose length is size(Y, 1).  The
%   vectors in X and Y are assumed to be organized as rows.  The
%   input data can be real or complex. The outputs are real
%   quantities.
%
%   D = MAHALANOBIS(Y, CX, MX) computes the Mahalanobis distance
%   between each vector in Y and the given mean vector, MX. The
%   results are output in vector D, whose length is size(Y, 1).  The
%   vectors in Y are assumed to be organized as the rows of this
%   array. The input data can be real or complex. The outputs are
%   real quantities. In addition to the mean vector MX, the
%   covariance matrix CX of a population of vectors X also must be
%   provided. Use function COVMATRIX (Section 11.5) to compute MX and
%   CX.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.6 $  $Date: 2005/01/18 13:44:47 $

% Reference: Acklam, P. J. [2002]. "MATLAB Array Manipulation Tips
% and Tricks." Available at
%     home.online.no/~pjacklam/matlab/doc/mtt/index.html
% or at
%     www.prenhall.com/gonzalezwoodseddins

param = varargin; % Keep in mind that param is a cell array.
Y = param{1};
ny = size(Y, 1); % Number of vectors in Y.

if length(param) == 2
   X = param{2};
   % Compute the mean vector and covariance matrix of the vectors
   % in X.
   [Cx, mx] = covmatrix(X);
elseif length(param) == 3 % Cov. matrix and mean vector provided.
   Cx = param{2};
   mx = param{3};
else
   error('Wrong number of inputs.')
end
mx = mx(:)'; % Make sure that mx is a row vector.

% Subtract the mean vector from each vector in Y.
Yc = Y - mx(ones(ny, 1), :);

% Compute the Mahalanobis distances.
d = real(sum(Yc/Cx.*conj(Yc), 2));

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
% EEGPLOT2EVENT - convert EEGPLOT rejections into events
%                   compatible with the eeg_eegrej function for rejecting
%                   continuous portions of datasets.
%
% Usage:
%   >> [events] = eegplot2event( eegplotrej, type, colorin, colorout );
%
% Inputs:
%   eegplotrej - EEGPLOT output (TMPREJ; see eegplot for more details)
%   type       - type of the event. Default -1.
%   colorin    - only extract rejection of specific colors (here a n x 3
%                array must be given). Default: extract all rejections.
%   colorout   - do not extract rejection of specified colors.
%
% Outputs:
%   events     - array of events compatible with the eeg_eegrej function
%                for rejecting continuous portions of datasets.
%
% Example:
% NEWEEG = eeg_eegrej(EEG,eegplot2event(TMPREJ, -1));
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: EEGPLOT, EEG_MULTIEEGPLOT, EEGPLOT2TRIAL, EEGLAB

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function events = eegplot2event( TMPREJ, type, color, colorout )

if nargin < 1
   help eegplot2event;
   return;
end;

if nargin < 2
   type = -1;
end;

% only take into account specific colors
% --------------------------------------
if exist('color') == 1
	for index = 1:size(color,1)
   		tmpcol1 = TMPREJ(:,3) + 255*TMPREJ(:,4) + 255*255*TMPREJ(:,5);
   		tmpcol2 = color(index,1)+255*color(index,2)+255*255*color(index,3);
   		I = find( tmpcol1 == tmpcol2);
   		if isempty(I)
      		fprintf('Warning: color [%d %d %d] not found\n', ...
      			color(index,1), color(index,2), color(index,3));
   		end
   		TMPREJ = TMPREJ(I,:);
   	end;
end

% remove other specific colors
% ----------------------------
if exist('colorout') == 1
	for index = 1:size(colorout,1)
   		tmpcol1 = TMPREJ(:,3) + 255*TMPREJ(:,4) + 255*255*TMPREJ(:,5);
   		tmpcol2 = colorout(index,1)+255*colorout(index,2)+255*255*colorout(index,3);
   		I = find( tmpcol1 ~= tmpcol2);
   		TMPREJ = TMPREJ(I,:);
   	end;
end

events = [];
if ~isempty(TMPREJ)
   events = TMPREJ(:,1:5);
   events = [ type*ones(size(events,1), 1) ones(size(events,1), 1) round(events(:,1:2)) events(:,3:5)];
end

return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function h = p31_h ( n, x )

%*****************************************************************************80
%
%% P31_H evaluates the Hessian for problem 31.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real H(N,N), the N by N Hessian matrix.
%
  h = zeros ( n, n );

  m = 5;
  a = [ 4.0, 4.0, 4.0, 4.0; ...
        1.0, 1.0, 1.0, 1.0; ...
        8.0, 8.0, 8.0, 8.0; ...
        6.0, 6.0, 6.0, 6.0; ...
        3.0, 7.0, 3.0, 7.0 ]';

  c = [ 0.1, 0.2, 0.2, 0.4, 0.6 ]';

  for ii = 1 : n
    for jj = 1 : n
      for j = 1 : m
        d = c(j) + sum ( ( x(1:n) - a(1:n,j) ).^2 );
        h(ii,jj) = h(ii,jj) - 8.0 * ( x(ii) - a(ii,j) ) * ( x(jj) - a(jj,j) ) / d^3;
        if ( ii == jj )
          h(ii,jj) = h(ii,jj) + 2.0 / d^2;
        end
      end
    end
  end

  return
end


text:
function [TruePositives, FalseNegatives] = testFermentationProducts(model, microbeID, biomassReaction, database, inputDataFolder)
% Performs an FVA and reports those fermentation products (exchange reactions)
% that can be secreted by the model and should be secreted according to
% data (true positives) and those fermentation products that cannot be secreted by
% the model but should be secreted according to in vitro data (false
% negatives).
%
% INPUT
% model             COBRA model structure
% microbeID         Microbe ID in carbon source data file
% biomassReaction   Biomass objective functions (low flux through BOF
%                   required in analysis)
% database          Structure containing rBioNet reaction and metabolite
%                   database
% inputDataFolder   Folder with experimental data and database files
%                   to load
%
% OUTPUT
% TruePositives     Cell array of strings listing all fermentation products
%                   (exchange reactions) that can be secreted by the model
%                   and in in vitro data.
% FalseNegatives    Cell array of strings listing all fermentation products
%                   (exchange reactions) that cannot be secreted by the model
%                   but should be secreted according to in vitro data.
%
% .. Author:
%       Stefania Magnusdottir, Nov 2017
%       Almut Heinken, Jan 2018-reduced number of reactions minimized and
%                      maximized to speed up the computation
%                      March 2022 - changed code to string-matching to make it
%                      more robust

global CBT_LP_SOLVER
if isempty(CBT_LP_SOLVER)
    initCobraToolbox
end

% read fermentation product table
dataTable = readInputTableForPipeline([inputDataFolder filesep 'FermentationTable.txt']);

% remove the reference columns
dataTable(:,find(strncmp(dataTable(1,:),'Ref',3))) = [];

corrRxns = {'Acetate kinase (acetate producer or consumer)','EX_ac(e)';'Bifid shunt','EX_ac(e)';'Acetogen pathway','EX_ac(e)';'Formate producer','EX_for(e)';'D-lactate producer','EX_lac_D(e)';'L-lactate producer','EX_lac_L(e)';'Ethanol producer or consumer','EX_etoh(e)';'Succinate producer','EX_succ(e)';'Propionate from succinate','EX_ppa(e)';'Propionate from propane-1,2-diol','EX_ppa(e)';'Propionate from lactate (acrylate pathway)','EX_ppa(e)';'Propionate from threonine','EX_ppa(e)';'Wood-Werkman cycle','EX_ppa(e)';'Butyrate via butyryl-CoA: acetate CoA transferase','EX_but(e)';'Butyrate via butyrate kinase','EX_but(e)';'Butyrate from lysine via butyrate-acetoacetate CoA-transferase','EX_but(e)';'Butyrate from glutarate or glutamate','EX_but(e)';'Butyrate from 4-hydroxybutyrate or succinate','EX_but(e)';'Hydrogen from ferredoxin oxidoreductase','EX_h2(e)';'Hydrogen from formate hydrogen lyase','EX_h2(e)';'Methanogenesis','EX_ch4(e)';'Sulfate reducer','EX_h2s(e)';'Isobutyrate producer','EX_isobut(e)';'Isovalerate producer','EX_isoval(e)';'Acetoin producer','EX_actn_R(e)';'2,3-butanediol producer','EX_btd_RR(e)';'Indole producer','EX_indole(e)';'Phenylacetate producer','EX_pac(e)';'Butanol producer','EX_btoh(e)';'Valerate producer','EX_M03134(e)'};

TruePositives = {};  % true positives (uptake in vitro and in silico)
FalseNegatives = {};  % false negatives (uptake in vitro not in silico)

% find microbe index in fermentation table
mInd = find(strcmp(dataTable(:,1), microbeID));
if isempty(mInd)
    warning(['Microbe "', microbeID, '" not found in fermentation product data file.'])
else
    % perform FVA to identify uptake metabolites
    % set BOF
    if ~any(ismember(model.rxns, biomassReaction)) || nargin < 3
        error(['Biomass reaction "', biomassReaction, '" not found in model.'])
    end
    model = changeObjective(model, biomassReaction);
    % set a low lower bound for biomass
    %     model = changeRxnBounds(model, biomassReaction, 1e-3, 'l');
    % list exchange reactions
    exchanges = model.rxns(strncmp('EX_', model.rxns, 3));
    % open all exchanges
    model = changeRxnBounds(model, exchanges, -1000, 'l');
    model = changeRxnBounds(model, exchanges, 1000, 'u');

    % get the reactions to test
    rxns = {};
    for i=2:size(dataTable,2)
        if contains(version,'(R202') % for Matlab R2020a and newer
            if dataTable{mInd,i}==1
                findCorrRxns = find(strcmp(corrRxns(:,1),dataTable{1,i}));
                rxns = union(rxns,corrRxns(findCorrRxns,2:end));
            end
        else
            if strcmp(dataTable{mInd,i},'1')
                findCorrRxns = find(strcmp(corrRxns(:,1),dataTable{1,i}));
                rxns = union(rxns,corrRxns(findCorrRxns,2:end));
            end
        end
    end

    % flux variability analysis on reactions of interest
    rxns = unique(rxns);
    rxns = rxns(~cellfun('isempty', rxns));
    rxnsInModel=intersect(rxns,model.rxns);
    if ~isempty(rxnsInModel)
        currentDir=pwd;
        try
            [minFlux, maxFlux, ~, ~] = fastFVA(model, 0, 'max', 'ibm_cplex', ...
                rxnsInModel, 'S');
        catch
            warning('fastFVA could not run, so fluxVariability is instead used. Consider installing fastFVA for shorter computation times.');
            cd(currentDir)
            [minFlux, maxFlux] = fluxVariability(model, 0, 'max', rxnsInModel);
        end

        % active flux
        flux = rxnsInModel(maxFlux > 1e-6);
    else
        flux = {};
    end

    % which reaction should carry flux according to in vitro data
    for i=2:size(dataTable,2)
        rxn={};
        if contains(version,'(R202') % for Matlab R2020a and newer
            if dataTable{mInd,i}==1
                rxn = corrRxns{find(strcmp(corrRxns(:,1),dataTable{1,i})),2};
            end
        else
            if strcmp(dataTable{mInd,i},'1')
                rxn = corrRxns{find(strcmp(corrRxns(:,1),dataTable{1,i})),2};
            end
        end
        if ~isempty(rxn)
            % add any that are not in model/not carrying flux to the false negatives
            if ~isempty(intersect(rxn,flux))
                TruePositives = union(TruePositives,rxn);
            else
                FalseNegatives=union(FalseNegatives,rxn);
            end
        end
    end
end

% replace reaction IDs with metabolite names
if ~isempty(TruePositives)
    TruePositives = TruePositives(~cellfun(@isempty, TruePositives));
    TruePositives=strrep(TruePositives,'EX_','');
    TruePositives=strrep(TruePositives,'(e)','');

    for i=1:length(TruePositives)
        TruePositives{i}=database.metabolites{find(strcmp(database.metabolites(:,1),TruePositives{i})),2};
    end
end

% warn about false negatives
if ~isempty(FalseNegatives)
    FalseNegatives = FalseNegatives(~cellfun(@isempty, FalseNegatives));
    FalseNegatives=strrep(FalseNegatives,'EX_','');
    FalseNegatives=strrep(FalseNegatives,'(e)','');
    for i = 1:length(FalseNegatives)
        FalseNegatives{i}=database.metabolites{find(strcmp(database.metabolites(:,1),FalseNegatives{i})),2};
    end
end

end


text:
function element_num = grid_element_num ( code, nelemx, nelemy )

%*****************************************************************************80
%
%% GRID_ELEMENT_NUM returns the number of elements in a grid.
%
%  Discussion:
%
%    The number of elements generated will be NELEMX * NELEMY for
%    quadrilaterals, or 2 * NELEMX * NELEMY for triangles.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string CODE, identifies the element desired.
%    Legal values include 'Q4', 'Q8', 'Q9', 'Q12', 'Q16', 'QL', 'T3',
%    'T4', 'T6' and 'T10'.
%
%    Input, integer NELEMX, NELEMY, the number of quadrilaterals along the
%    X and Y directions.
%
%    Output, integer ELEMENT_NUM, the number of elements in the grid.
%
  if ( s_eqi ( code, 'Q4' ) )
    element_num = grid_q4_element_num ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'Q8' ) )
    element_num = grid_q8_element_num ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'Q9' ) )
    element_num = grid_q9_element_num ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'Q12' ) )
    element_num = grid_q12_element_num ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'Q16' ) )
    element_num = grid_q16_element_num ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'QL' ) )
    element_num = grid_ql_element_num ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'T3' ) )
    element_num = grid_t3_element_num ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'T4' ) )
    element_num = grid_t4_element_num ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'T6' ) )
    element_num = grid_t6_element_num ( nelemx, nelemy );
  elseif ( s_eqi ( code, 'T10' ) )
    element_num = grid_t10_element_num ( nelemx, nelemy );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GRID_ELEMENT_NUM - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of CODE = "%s".', code );
    element_num = -1;
    error ( 'GRID_ELEMENT_NUM - Fatal error!' )
  end

  return
end


text:
function [FcellZ, FcellNeuZ] = extractZSignals(ops, zpos, Fcell, FcellNeu, subpixel)

Nk = size(Fcell{1},1);

F = [];
Fneu = [];
for k = 1:numel(Fcell)
    F = [F Fcell{k}];
    Fneu = [Fneu FcellNeu{k}];
end

% compute each cell's profile in Z
zspread  = 10;
iZ = [-zspread : 1/subpixel : zspread];
zind = zpos * subpixel + find(iZ==0);
zF = zeros(Nk, numel(iZ));
zFneu = zeros(Nk, numel(iZ));
for j = 1:numel(iZ)
    zF(:,j)    = mean(F(:, zind == j), 2);
    zFneu(:,j) = mean(Fneu(:, zind == j), 2);
end

if sum(zind==1) == 0
    jm = min(zind);
    zF(:, 1:jm-1) = repmat(zF(:, jm), 1, jm-1);
    zFneu(:, 1:jm-1) = repmat(zFneu(:, jm), 1, jm-1);
end

if sum(zind==numel(iZ)) == 0
    jm = max(zind);
    zF(:, jm+1:end) = repmat(zF(:, jm), 1, numel(iZ) - jm);
    zFneu(:, jm+1:end) = repmat(zFneu(:, jm), 1, numel(iZ) - jm);
end

zF= fixnangaps(zF);
zFneu= fixnangaps(zFneu);

zF = my_conv2(zF, 0.5, 2);
zFneu = my_conv2(zFneu, 0.5, 2);
plot(zF(1,:));

%%
Fz      = zF(:, zind);
Fzneu    = zFneu(:, zind);

csumNframes = [0 cumsum(ops.Nframes)];
FcellZ       = cell(1, length(ops.Nframes));
FcellNeuZ    = cell(1, length(ops.Nframes));
for i = 1:length(ops.Nframes)
    FcellZ{i}     = Fz(:, csumNframes(i) + (1:ops.Nframes(i)));
    FcellNeuZ{i}  = Fzneu(:, csumNframes(i) + (1:ops.Nframes(i)));
end

author: cortex-lab


text:
function centroid = findBallFcn(greenBall1, thresh, imageType, axH)

% Copyright 2011 The MathWorks, Inc.

error(nargchk(3, 4, nargin, 'struct'));

if isempty(greenBall1)
    return;
end
if ischar(greenBall1)
    greenBall1 = imread(greenBall1);
end

%% Find Green Object
% This script reads in an image file and then attempts to find a green
% object in the image. It is designed to find one green ball and highlight
% that ball on the original image

%% Extract each color
% Next we using indexing to extract three 2D matrices from the 3D image
% data corresponding to the red, green, and blue components of the image.

r = greenBall1(:, :, 1);
g = greenBall1(:, :, 2);
b = greenBall1(:, :, 3);

%% Calculate Green
% Then we perform an arithmetic operation on the matrices as a whole to try
% to create one matrix that represents an intensity of green.

justGreen = g - r/2 - b/2;

%% Threshold the image
% Now we can set a threshold to separate the parts of the image that we
% consider to be green from the rest.

if nargin == 4
    bw = justGreen > thresh;
else
    bw = justGreen > 80;
end

%% Remove small groups
% We can use special functions provided by the Image Processing toolbox to
% quickly perform common image processing tasks. Here we are using
% BWAREAOPEN to remove groups of pixels less than 30.

ball1 = bwareaopen(bw, 30);

%% Find center
% Now we are using REGIONPROPS to extract the centroid of the group of
% pixels representing the ball.

s  = regionprops(ball1, {'centroid','area'});
if isempty(s)
    centroid = [];
else
    [maxArea, id] = max([s.Area]); %#ok<ASGLU>
    centroid = s(id).Centroid;
end
switch imageType
   case 'video'
      imshow(greenBall1, 'Parent', axH);
      if ~isempty(centroid)
         line(centroid(1), centroid(2), 'Parent', axH, 'Color', 'w', 'Marker', 'p', 'MarkerSize', 20, 'MarkerFaceColor', 'r')
      end
   case 'bw'
      imshow(ball1, 'Parent', axH);
end

author: Sable

repo: mcbench-benchmarks


text:
function [ y, m, d, f ] = jed_to_ymdf_republican ( jed )

%*****************************************************************************80
%
%% JED_TO_YMDF_REPUBLICAN converts a JED to a Republican YMDF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Richards,
%    Algorithm F,
%    Mapping Time, The Calendar and Its History,
%    Oxford, 1999, pages 324-325.
%
%  Parameters:
%
%    Input, real JED, the Julian Ephemeris Date.
%
%    Output, integer Y, M, D, real F,
%    the YMDF date.
%

%
%  Determine the computational date (Y'/M'/D').
%
  j = floor ( jed + 0.5 );
  f = ( jed + 0.5 ) - j;

  g =floor ( ( 4 * j + 578797 ) / 146097 );
  g = floor ( ( 3 * g ) / 4 ) - 51;
  j_prime = j + 111 + g;

  y_prime =     floor ( ( 4 * j_prime + 3 ) / 1461 );
  t_prime = floor ( mod ( 4 * j_prime + 3, 1461 ) / 4 );
  m_prime = floor ( t_prime / 30 );
  d_prime =   mod ( t_prime, 30 );
%
%  Convert the computational date to a calendar date.
%
  d = d_prime + 1;
  m = mod ( m_prime, 13 ) + 1;
  y = y_prime - 6504 + floor ( ( 13 - m ) / 13 );

  return
end


text:
function [center,radius] = minboundsphere(xyz)
% minboundsphere: Compute the minimum radius enclosing sphere of a set of (x,y,z) triplets
% usage: [center,radius] = minboundsphere(xyz)
%
% arguments: (input)
%  xyz - nx3 array of (x,y,z) triples, describing points in R^3
%        as rows of this array.
%
%
% arguments: (output)
%  center - 1x3 vector, contains the (x,y,z) coordinates of
%        the center of the minimum radius enclosing sphere
%
%  radius - scalar - denotes the radius of the minimum
%        enclosing sphere
%
%
% Example usage:
% Sample uniformly from the interior of a unit sphere.
% As the sample size increases, the enclosing sphere
% should asymptotically approach center = [0 0 0], and
% radius = 1.
%
%   xyz = rand(10000,3)*2-1;
%   r = sqrt(sum(xyz.^2,2));
%   xyz(r>1,:) = [];          % 5156 points retained
%   tic,[center,radius] = minboundsphere(xyz);toc
%
%   Elapsed time is 0.199467 seconds.
%
%   center = [0.00017275   8.5006e-05   0.00012015]
%
%   radius = 0.9999
%
% Example usage:
% Sample from the surface of a unit sphere. Within eps
% or so, the result should be center = [0 0 0], and radius = 1.
%
%   xyz = randn(10000,3);
%   xyz = xyz./repmat(sqrt(sum(xyz.^2,2)),1,3);
%   tic,[center,radius] = minboundsphere(xyz);toc
%
%   Elapsed time is 0.614762 seconds.
%
%   center =
%      4.6127e-17  -2.5584e-17   7.2711e-17
%
%   radius =
%       1
%
%
% See also: minboundrect, minboundcircle
%
%
% Author: John D'Errico
% E-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 1/23/07

% not many error checks to worry about
sxyz = size(xyz);
if (length(sxyz)~=2) || (sxyz(2)~=3)
  error 'xyz must be an nx3 array of points'
end
n = sxyz(1);

% start out with the convex hull of the points to
% reduce the problem dramatically. Note that any
% points in the interior of the convex hull are
% never needed.
if n>4
  tri = convhulln(xyz,{'QJ' 'Pp'});

  % list of the unique points on the convex hull itself
  hlist = unique(tri(:));

  % exclude those points inside the hull as not relevant
  xyz = xyz(hlist,:);

end

% now we must find the enclosing sphere of those that
% remain.
n = size(xyz,1);

% special case small numbers of points. If we trip any
% of these cases, then we are done, so return.
switch n
  case 0
    % empty begets empty
    center = [];
    radius = [];
    return
  case 1
    % with one point, the center has radius zero
    center = xyz;
    radius = 0;
    return
  case 2
    % only two points. center is at the midpoint
    center = mean(xyz,1);
    radius = norm(xyz(1,:) - center);
    return
  case 3
    % exactly 3 points. For this odd case, just use enc4,
    % appending a new point at the centroid. This is simpler
    % than other solutions that would have reduced the
    % problem to 2-d. enc4 will do that anyway.
    [center,radius] = enc4([xyz;mean(xyz,1)]);
    return
  case 4
    % exactly 4 points
    [center,radius] = enc4(xyz);
    return
end

% pick a tolerance
tol = 10*eps*max(max(abs(xyz),[],1) - min(abs(xyz),[],1));

% more than 4 points. for no more than 15 points in the hull,
% just do an exhaustive search.
if n <= 15
  % for 15 points, there are only nchoosek(15,4) = 1365
  % sets to look through. this is only about a second.
  asets = nchoosek(1:n,4);

  center = inf(1,3);
  radius = inf;
  for i = 1:size(asets,1)
    aset = asets(i,:);
    iset = setdiff(1:n,aset);

    % get the enclosing sphere for the current set
    [centeri,radiusi] = enc4(xyz(aset,:));

    % are all the inactive set points inside the circle?
    ri = sqrt(sum((xyz(iset,:) - repmat(centeri,n-4,1)).^2,2));

    [rmax,k] = max(ri);
    if ((rmax - radiusi) <= tol) && (radiusi < radius)
      center = centeri;
      radius = radiusi;
    end
  end

else
  % Use an active set strategy, on many different
  % random starting sets.
  center = inf(1,3);
  radius = inf;

  for i = 1:250
    aset = randperm(n); % a random start, but quite adequate
    iset = aset(5:n);
    aset = aset(1:4);

    flag = true;
    iter = 0;
    centeri = inf(1,3);
    radiusi = inf;
    while flag && (iter < 12)
      iter = iter + 1;

      % get the enclosing sphere for the current set
      [centeri,radiusi] = enc4(xyz(aset,:));

      % are all the inactive set points inside the circle?
      ri = sqrt(sum((xyz(iset,:) - repmat(centeri,n-4,1)).^2,2));

      [rmax,k] = max(ri);
      if (rmax - radiusi) <= tol
        % the active set enclosing sphere also enclosed
        % all of the inactive points. We are done.
        flag = false;
      else
        % it must be true that we can replace one member of aset
        % with iset(k). That k'th element was farthest out, so
        % it seems best (a greedy algorithm) to swap it in. The
        % problem with the greedy algorithm, is it gets trapped
        % in a cycle at times. but since we are restarting the
        % algorithm multiple times, this will work.
        s1 = [aset([2 3 4]),iset(k)];
        [c1,r1] = enc4(xyz(s1,:));
        if (norm(c1 - xyz(aset(1),:)) <= r1)
          centeri = c1;
          radiusi = r1;

          % update the active/inactive sets
          swap = aset(1);
          aset = [iset(k),aset([2 3 4])];
          iset(k) = swap;

          % bounce out to the while loop
          continue
        end
        s1 = [aset([1 3 4]),iset(k)];
        [c1,r1] = enc4(xyz(s1,:));
        if (norm(c1 - xyz(aset(2),:)) <= r1)
          centeri = c1;
          radiusi = r1;

          % update the active/inactive sets
          swap = aset(2);
          aset = [iset(k),aset([1 3 4])];
          iset(k) = swap;

          % bounce out to the while loop
          continue
        end
        s1 = [aset([1 2 4]),iset(k)];
        [c1,r1] = enc4(xyz(s1,:));
        if (norm(c1 - xyz(aset(3),:)) <= r1)
          centeri = c1;
          radiusi = r1;

          % update the active/inactive sets
          swap = aset(3);
          aset = [iset(k),aset([1 2 4])];
          iset(k) = swap;

          % bounce out to the while loop
          continue
        end
        s1 = [aset([1 2 3]),iset(k)];
        [c1,r1] = enc4(xyz(s1,:));
        if (norm(c1 - xyz(aset(4),:)) <= r1)
          centeri = c1;
          radiusi = r1;

          % update the active/inactive sets
          swap = aset(4);
          aset = [iset(k),aset([1 2 3])];
          iset(k) = swap;

          % bounce out to the while loop
          continue
        end

        % if we get through to this point, then something went wrong.
        % Active set problem. Increase tol, then try again.
        tol = 2*tol;

      end
    end

    % have we improved over the best set so far?
    if radiusi < radius
      center = centeri;
      radius = radiusi;
    end
  end
end

% =======================================
%  begin subfunctions
% =======================================
function [center,radius] = enc4(xyz)
% minimum radius enclosing sphere for exactly 4 points in R^3
%
% xyz is a 4x3 array
%
% Note that enc4 will attempt to pass a sphere through all
% 4 of the supplied points. When the set of points proves to
% be degenerate, perhaps because of collinearity of 3 or
% more of the points, or because the 4 points are coplanar,
% then the sphere would nominally have infinite radius. Since
% there must be a finite radius sphere to enclose any set of
% finite valued points, enc4 will provide that sphere instead.
%
% In addition, there are some non-degenerate sets of points
% for which the circum-sphere is not minimal. enc4 will always
% try to find the minimum radius enclosing sphere.

% interpoint distance matrix D
% dfun = @(A) (A(:,[1 1 1 1]) - A(:,[1 1 1 1])').^2;
dfun = inline('(A(:,[1 1 1 1]) - A(:,[1 1 1 1])'').^2','A');
D = sqrt(dfun(xyz(:,1)) + dfun(xyz(:,2)) + dfun(xyz(:,3)));

% Find the most distant pair. Test if their circum-sphere
% also encloses the other points. If it does, then we are
% done.
[dij,ij] = max(D(:));
[i,j] = ind2sub([4 4],ij);
others = setdiff(1:4,[i,j]);
radius = dij/2;
center = (xyz(i,:) + xyz(j,:))/2;
if (norm(center - xyz(others(1),:))<=radius) && ...
   (norm(center - xyz(others(2),:))<=radius)
  % we can stop here.
  return
end

% next, we need to test each triplet of points, finding their
% enclosing sphere. If the 4th point is also inside, then we
% are done.
ind = 1:3;
[center,radius,isin] = enc3_4(xyz(ind,:),xyz(4,:),D(ind,ind));
if isin
  % the 4th point was inside this enclosing sphere.
  return
end

ind = [1 2 4];
[center,radius,isin] = enc3_4(xyz(ind,:),xyz(3,:),D(ind,ind));
if isin
  % the 3rd point was inside this enclosing sphere.
  return
end

ind = [1 3 4];
[center,radius,isin] = enc3_4(xyz(ind,:),xyz(2,:),D(ind,ind));
if isin
  % the second point was inside this enclosing sphere.
  return
end

ind = [2 3 4];
[center,radius,isin] = enc3_4(xyz(ind,:),xyz(1,:),D(ind,ind));
if isin
  % the first point was inside this enclosing sphere.
  return
end

% find the circum-sphere that passes through all 4 points
% since we have passed all the other tests, we need not
% worry here about singularities in the system of
% equations.
A = 2*(xyz(2:4,:)-repmat(xyz(1,:),3,1));
rhs = sum(xyz(2:4,:).^2 - repmat(xyz(1,:).^2,3,1),2);
center = (A\rhs)';
radius = norm(center - xyz(1,:));

% =======================================
function [center,radius,isin] = enc3_4(xyz,xyztest,Di)
% minimum radius enclosing sphere for exactly 3 points in R^3
%
% xyz - a 3x3 array, with each row as a point in R^3
%
% xyztest - 1x3 vector, a point to be tested if it is
%       inside the generated enclosing sphere.
%
% Di - 3x3 array of interpoint distances

% test the farthest pair of points. do they form a diameter
% of the sphere?
if Di(1,2)>=max(Di(1,3),Di(2,3))
  center = mean(xyz([1 2],:),1);
  radius = Di(1,2)/2;
  isin = (norm(xyz(3,:) - center)<=radius) && (norm(xyztest - center)<=radius);
elseif Di(1,3)>=max(Di(1,2),Di(2,3))
  center = mean(xyz([1 3],:),1);
  radius = Di(1,3)/2;
  isin = (norm(xyz(2,:) - center)<=radius) && (norm(xyztest - center)<=radius);
elseif Di(2,3)>=max(Di(1,2),Di(1,3))
  center = mean(xyz([2 3],:),1);
  radius = Di(2,3)/2;
  isin = (norm(xyz(1,:) - center)<=radius) && (norm(xyztest - center)<=radius);
end
if isin
  % we found the minimal enclosing sphere already
  return
end

% If we drop down to here, no singularities should
% happen (I've already caught any degeneracies.)

% We transform the three points into a plane, then
% compute the enclosing sphere in that plane.

% translate to the origin
xyz0 = xyz(1,:);
xyzt = xyz(2:3,:) - [xyz0;xyz0];

rot = orth(xyzt');

% uv is composed of 2 points, in 2-d, plus we
% have the origin (after the translation)
uv = xyzt*rot;

A = 2*uv;
rhs = sum(uv.^2,2);
center = (A\rhs)';
radius = norm(center - uv(1,:));

% rotate and translate back
center = center*rot' + xyz0;

% test if the 4th point is enclosed also
isin = (norm(xyztest - center)<=radius);


text:
function classString = classIDToString(classID)

labels = getClassLabels;

if classID<1 || classID>length(labels)
    classString='unknown';
else
    classString = char(labels{classID});
end

end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function [refBody, scName] = orbitPanelGetOrbitFromSFSContextCallBack(mainGUIHandle, hSMA, hECC, hINC, hRAAN, hARG, varargin)
%orbitPanelGetOrbitFromSFSContextCallBack Summary of this function goes here
%   Detailed explanation goes here
    if(not(isempty(mainGUIHandle)))
        mainUserData = get(mainGUIHandle,'UserData');
        prevPathName = mainUserData{3,1};
    else
        mainUserData = [];
        prevPathName = [];
    end

    refBody = [];
%     [orbit,PathName,scName] = importOrbitGUI(1, prevPathName);

    output = AppDesignerGUIOutput({[],[],''});
    importOrbitGUI_App(1, prevPathName, output);
    orbit = output.output{1};
    PathName = output.output{3};
    scName = output.output{3};

    if(not(isempty(orbit)))
        if(not(isempty(hSMA)))
            set(hSMA, 'String', fullAccNum2Str(orbit{3}));
        end

        if(not(isempty(hECC)))
            set(hECC, 'String', fullAccNum2Str(orbit{4}));
        end

        if(not(isempty(hINC)))
            set(hINC, 'String', fullAccNum2Str(orbit{5}));
        end

        if(not(isempty(hRAAN)))
            set(hRAAN, 'String', fullAccNum2Str(orbit{6}));
        end

        if(not(isempty(hARG)))
            set(hARG, 'String', fullAccNum2Str(orbit{7}));
        end

        refBody = orbit{10};

        if(nargin == 8)
            hMA = varargin{1};
            hEpoch = varargin{2};
            set(hMA, 'String', fullAccNum2Str(orbit{8}));
            set(hEpoch, 'String', fullAccNum2Str(orbit{9}));
        end

        if(not(isempty(mainUserData)))
            mainUserData{3,1} = PathName;
            set(mainGUIHandle,'UserData',mainUserData);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function varargout = welcome(varargin)
clc;

disp('                     HELLO ;)');
disp('                     @)=--   ')
% WELCOME M-file for welcome.fig
%      WELCOME, by itself, creates a new WELCOME or raises the existing
%      singleton*.
%
%      H = WELCOME returns the handle to a new WELCOME or the handle to
%      the existing singleton*.
%
%      WELCOME('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in WELCOME.M with the given input arguments.
%
%      WELCOME('Property','Value',...) creates a new WELCOME or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before welcome_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to welcome_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help welcome

% Last Modified by GUIDE v2.5 23-Jan-2008 08:50:01

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @welcome_OpeningFcn, ...
                   'gui_OutputFcn',  @welcome_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before welcome is made visible.
function welcome_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to welcome (see VARARGIN)

% Choose default command line output for welcome
handles.output = hObject;

% printing the top text
t=clock;
if t(4)>1 && t(4)<=12
    temp = 1;
elseif t(4)>12 && t(4)<=15
    temp = 2;
else
    temp = 3;
end
switch temp
    case 1

        set(handles.text1,'string','GOOD MORNING. WELCOME TO THE WORLD OF MATLAB ');
        backgroundImage = importdata('good_morning.jpg');
        %select the axes
        axes(handles.axes1);
        %place image onto the axes
        image(backgroundImage);
        %remove the axis tick marks
        axis off

    case 2
        set(handles.text1,'string','GOOD AFTERNOON. WELCOME TO THE WORLD OF MATLAB ');
        backgroundImage = importdata('good_afternoon.jpg');
        %select the axes
        axes(handles.axes1);
        %place image onto the axes
        image(backgroundImage);
        %remove the axis tick marks
        axis off

    case 3
       set(handles.text1,'string','GOOD EVENING. WELCOME TO THE WORLD OF MATLAB ');
       backgroundImage = importdata('good_evening.jpg');
        %select the axes
        axes(handles.axes1);
        %place image onto the axes
        image(backgroundImage);
        %remove the axis tick marks
        axis off
end

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes welcome wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = welcome_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
close all;
disp('Have a nice Day :)');

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
simulink
new_system('New','Model');
open_system('New');
close all;


text:
function range = p17_range ( )

%*****************************************************************************80
%
%% P17_RANGE returns an interval bounding the root for problem 17.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real RANGE(2), the minimum and maximum values of
%    an interval containing the root.
%
  range(1) = 2.0;
  range(2) = 3.0;

  return
end


text:
classdef ScrollingPanel < uix.Container & uix.mixin.Panel
    %uix.ScrollingPanel  Scrolling panel
    %
    %  p = uix.ScrollingPanel(p1,v1,p2,v2,...) constructs a scrolling panel
    %  and sets parameter p1 to value v1, etc.
    %
    %  A scrolling panel is a standard container (uicontainer) that shows
    %  one its contents and hides the others.
    %
    %  See also: uix.Panel, uix.BoxPanel, uix.TabPanel, uicontainer

    %  Copyright 2009-2016 The MathWorks, Inc.
    %  $Revision: 1775 $ $Date: 2018-12-06 06:07:40 +0000 (Thu, 06 Dec 2018) $

    properties( Dependent )
        Heights % heights of contents, in pixels and/or weights
        MinimumHeights % minimum heights of contents, in pixels
        VerticalOffsets % vertical offsets of contents, in pixels
        VerticalSteps % vertical slider steps, in pixels
        Widths % widths of contents, in pixels and/or weights
        MinimumWidths % minimum widths of contents, in pixels
        HorizontalOffsets % horizontal offsets of contents, in pixels
        HorizontalSteps % horizontal slider steps, in pixels
        MouseWheelEnabled % mouse wheel scrolling enabled [on|off]
    end

    properties( Access = protected )
        Heights_ = zeros( [0 1] ) % backing for Heights
        MinimumHeights_ = zeros( [0 1] ) % backing for MinimumHeights
        Widths_ = zeros( [0 1] ) % backing for Widths
        MinimumWidths_ = zeros( [0 1] ) % backing for MinimumWidths
        HorizontalSliders = matlab.ui.control.UIControl.empty( [0 1] ) % sliders
        VerticalSliders = matlab.ui.control.UIControl.empty( [0 1] ) % sliders
        BlankingPlates = matlab.ui.control.UIControl.empty( [0 1] ) % blanking plates
        HorizontalSteps_ = zeros( [0 1] ) % steps
        VerticalSteps_ = zeros( [0 1] ) % steps
    end

    properties( Access = private )
        MouseWheelListener = [] % mouse listener
        MouseWheelEnabled_ = 'on' %
        ScrollingListener = [] % slider listener
        ScrolledListener = [] % slider listener
    end

    properties( Constant, Access = protected )
        SliderSize = 20 % slider size, in pixels
        SliderStep = 10 % slider step, in pixels
    end

    events( NotifyAccess = private )
        Scrolling
        Scrolled
    end

    methods

        function obj = ScrollingPanel( varargin )
            %uix.ScrollingPanel  Scrolling panel constructor
            %
            %  p = uix.ScrollingPanel() constructs a scrolling panel.
            %
            %  p = uix.ScrollingPanel(p1,v1,p2,v2,...) sets parameter p1 to
            %  value v1, etc.

            % Set properties
            try
                uix.set( obj, varargin{:} )
            catch e
                delete( obj )
                e.throwAsCaller()
            end

        end % constructor

    end % structors

    methods

        function value = get.Heights( obj )

            value = obj.Heights_;

        end % get.Heights

        function set.Heights( obj, value )

            % For those who can't tell a column from a row...
            if isrow( value )
                value = transpose( value );
            end

            % Check
            assert( isa( value, 'double' ), 'uix:InvalidPropertyValue', ...
                'Property ''Heights'' must be of type double.' )
            assert( all( isreal( value ) ) && ~any( isinf( value ) ) && ...
                ~any( isnan( value ) ), 'uix:InvalidPropertyValue', ...
                'Elements of property ''Heights'' must be real and finite.' )
            assert( isequal( size( value ), size( obj.Contents_ ) ), ...
                'uix:InvalidPropertyValue', ...
                'Size of property ''Heights'' must match size of contents.' )

            % Set
            obj.Heights_ = value;

            % Mark as dirty
            obj.Dirty = true;

        end % set.Heights

        function value = get.MinimumHeights( obj )

            value = obj.MinimumHeights_;

        end % get.MinimumHeights

        function set.MinimumHeights( obj, value )

            % For those who can't tell a column from a row...
            if isrow( value )
                value = transpose( value );
            end

            % Check
            assert( isa( value, 'double' ), 'uix:InvalidPropertyValue', ...
                'Property ''MinimumHeights'' must be of type double.' )
            assert( all( isreal( value ) ) && ~any( isinf( value ) ) && ...
                all( value >= 0 ), 'uix:InvalidPropertyValue', ...
                'Elements of property ''MinimumHeights'' must be non-negative.' )
            assert( isequal( size( value ), size( obj.Heights_ ) ), ...
                'uix:InvalidPropertyValue', ...
                'Size of property ''MinimumHeights'' must match size of contents.' )

            % Set
            obj.MinimumHeights_ = value;

            % Mark as dirty
            obj.Dirty = true;

        end % set.MinimumHeights

        function value = get.VerticalOffsets( obj )

            sliders = obj.VerticalSliders;
            if isempty( sliders )
                value = zeros( size( sliders ) );
            else
                value = -vertcat( sliders.Value ) - 1;
                value(value<0) = 0;
            end

        end % get.VerticalOffsets

        function set.VerticalOffsets( obj, value )

            % Check
            assert( isa( value, 'double' ), 'uix:InvalidPropertyValue', ...
                'Property ''VerticalOffsets'' must be of type double.' )
            assert( all( isreal( value ) ) && ~any( isinf( value ) ) && ...
                ~any( isnan( value ) ), 'uix:InvalidPropertyValue', ...
                'Elements of property ''VerticalOffsets'' must be real and finite.' )
            assert( isequal( size( value ), size( obj.Contents_ ) ), ...
                'uix:InvalidPropertyValue', ...
                'Size of property ''VerticalOffsets'' must match size of contents.' )

            % Set
            sliders = obj.VerticalSliders;
            for ii = 1:numel( sliders )
                sliders(ii).Value = -value(ii) - 1;
            end

            % Mark as dirty
            obj.Dirty = true;

        end % set.VerticalOffsets

        function value = get.VerticalSteps( obj )

            value = obj.VerticalSteps_;

        end % get.VerticalSteps

        function set.VerticalSteps( obj, value )

            % For those who can't tell a column from a row...
            if isrow( value )
                value = transpose( value );
            end

            % Check
            assert( isa( value, 'double' ), 'uix:InvalidPropertyValue', ...
                'Property ''VerticalSteps'' must be of type double.' )
            assert( all( isreal( value ) ) && ~any( isinf( value ) ) && ...
                ~any( isnan( value ) ) && all( value > 0 ), ...
                'uix:InvalidPropertyValue', ...
                'Elements of property ''VerticalSteps'' must be real, finite and positive.' )
            assert( isequal( size( value ), size( obj.Contents_ ) ), ...
                'uix:InvalidPropertyValue', ...
                'Size of property ''VerticalSteps'' must match size of contents.' )

            % Set
            obj.VerticalSteps_ = value;

            % Mark as dirty
            obj.Dirty = true;

        end % set.VerticalSteps

        function value = get.Widths( obj )

            value = obj.Widths_;

        end % get.Widths

        function set.Widths( obj, value )

            % For those who can't tell a column from a row...
            if isrow( value )
                value = transpose( value );
            end

            % Check
            assert( isa( value, 'double' ), 'uix:InvalidPropertyValue', ...
                'Property ''Widths'' must be of type double.' )
            assert( all( isreal( value ) ) && ~any( isinf( value ) ) && ...
                ~any( isnan( value ) ), 'uix:InvalidPropertyValue', ...
                'Elements of property ''Widths'' must be real and finite.' )
            assert( isequal( size( value ), size( obj.Contents_ ) ), ...
                'uix:InvalidPropertyValue', ...
                'Size of property ''Widths'' must match size of contents.' )

            % Set
            obj.Widths_ = value;

            % Mark as dirty
            obj.Dirty = true;

        end % set.Widths

        function value = get.MinimumWidths( obj )

            value = obj.MinimumWidths_;

        end % get.MinimumWidths

        function set.MinimumWidths( obj, value )

            % For those who can't tell a column from a row...
            if isrow( value )
                value = transpose( value );
            end

            % Check
            assert( isa( value, 'double' ), 'uix:InvalidPropertyValue', ...
                'Property ''MinimumWidths'' must be of type double.' )
            assert( all( isreal( value ) ) && ~any( isinf( value ) ) && ...
                all( value >= 0 ), 'uix:InvalidPropertyValue', ...
                'Elements of property ''MinimumWidths'' must be non-negative.' )
            assert( isequal( size( value ), size( obj.Widths_ ) ), ...
                'uix:InvalidPropertyValue', ...
                'Size of property ''MinimumWidths'' must match size of contents.' )

            % Set
            obj.MinimumWidths_ = value;

            % Mark as dirty
            obj.Dirty = true;

        end % set.MinimumWidths

        function value = get.HorizontalOffsets( obj )

            sliders = obj.HorizontalSliders;
            if isempty( sliders )
                value = zeros( size( sliders ) );
            else
                value = vertcat( sliders.Value );
                value(value<0) = 0;
            end

        end % get.HorizontalOffsets

        function set.HorizontalOffsets( obj, value )

            % Check
            assert( isa( value, 'double' ), 'uix:InvalidPropertyValue', ...
                'Property ''HorizontalOffsets'' must be of type double.' )
            assert( all( isreal( value ) ) && ~any( isinf( value ) ) && ...
                ~any( isnan( value ) ), 'uix:InvalidPropertyValue', ...
                'Elements of property ''HorizontalOffsets'' must be real and finite.' )
            assert( isequal( size( value ), size( obj.Contents_ ) ), ...
                'uix:InvalidPropertyValue', ...
                'Size of property ''HorizontalOffsets'' must match size of contents.' )

            % Set
            sliders = obj.HorizontalSliders;
            for ii = 1:numel( sliders )
                sliders(ii).Value = value(ii);
            end

            % Mark as dirty
            obj.Dirty = true;

        end % set.HorizontalOffsets

        function value = get.HorizontalSteps( obj )

            value = obj.HorizontalSteps_;

        end % get.HorizontalSteps

        function set.HorizontalSteps( obj, value )

            % For those who can't tell a column from a row...
            if isrow( value )
                value = transpose( value );
            end

            % Check
            assert( isa( value, 'double' ), 'uix:InvalidPropertyValue', ...
                'Property ''HorizontalSteps'' must be of type double.' )
            assert( all( isreal( value ) ) && ~any( isinf( value ) ) && ...
                ~any( isnan( value ) ) && all( value > 0 ), ...
                'uix:InvalidPropertyValue', ...
                'Elements of property ''HorizontalSteps'' must be real, finite and positive.' )
            assert( isequal( size( value ), size( obj.Contents_ ) ), ...
                'uix:InvalidPropertyValue', ...
                'Size of property ''HorizontalSteps'' must match size of contents.' )

            % Set
            obj.HorizontalSteps_ = value;

            % Mark as dirty
            obj.Dirty = true;

        end % set.HorizontalSteps

        function value = get.MouseWheelEnabled( obj )

            value = obj.MouseWheelEnabled_;

        end % get.MouseWheelEnabled

        function set.MouseWheelEnabled( obj, value )

            assert( ischar( value ) && any( strcmp( value, {'on','off'} ) ), ...
                'uix:InvalidArgument', ...
                'Property ''MouseWheelEnabled'' must ''on'' or ''off''.' )
            listener = obj.MouseWheelListener;
            if ~isempty( listener )
                listener.Enabled = strcmp( value, 'on' );
            end
            obj.MouseWheelEnabled_ = value;

        end % set.MouseWheelEnabled

    end % accessors

    methods( Access = protected )

        function redraw( obj )
            %redraw  Redraw

            % Return if no contents
            selection = obj.Selection_;
            if selection == 0, return, end

            % Retrieve width and height of selected contents
            contentsWidth = obj.Widths_(selection);
            minimumWidth = obj.MinimumWidths_(selection);
            contentsHeight = obj.Heights_(selection);
            minimumHeight = obj.MinimumHeights_(selection);

            % Retrieve selected contents and corresponding decorations
            child = obj.Contents_(selection);
            vSlider = obj.VerticalSliders(selection);
            hSlider = obj.HorizontalSliders(selection);
            plate = obj.BlankingPlates(selection);

            % Compute dimensions
            bounds = hgconvertunits( ancestor( obj, 'figure' ), ...
                [0 0 1 1], 'normalized', 'pixels', obj );
            width = bounds(3);
            height = bounds(4);
            sliderSize = obj.SliderSize; % slider size
            vSliderWidth = sliderSize * ...
                (contentsHeight > height | ...
                minimumHeight > height); % first pass
            hSliderHeight = sliderSize * ...
                (contentsWidth > width - vSliderWidth | ...
                minimumWidth > width - vSliderWidth);
            vSliderWidth = sliderSize * ...
                (contentsHeight > height - hSliderHeight | ...
                minimumHeight > height - hSliderHeight); % second pass
            vSliderWidth = min( vSliderWidth, width ); % limit
            hSliderHeight = min( hSliderHeight, height ); % limit
            vSliderHeight = height - hSliderHeight;
            hSliderWidth = width - vSliderWidth;
            widths = uix.calcPixelSizes( width, ...
                [contentsWidth;vSliderWidth], ...
                [minimumWidth;vSliderWidth], 0, 0 );
            contentsWidth = widths(1); % to be offset
            heights = uix.calcPixelSizes( height, ...
                [contentsHeight;hSliderHeight], ...
                [minimumHeight;hSliderHeight], 0, 0 );
            contentsHeight = heights(1); % to be offset

            % Compute positions
            contentsPosition = [1 1+hSliderHeight+vSliderHeight-contentsHeight contentsWidth contentsHeight];
            vSliderPosition = [1+hSliderWidth 1+hSliderHeight vSliderWidth vSliderHeight];
            hSliderPosition = [1 1 hSliderWidth hSliderHeight];
            platePosition = [1+hSliderWidth 1 vSliderWidth hSliderHeight];

            % Compute and set vertical slider properties
            if vSliderWidth == 0 || vSliderHeight == 0 || vSliderHeight <= vSliderWidth
                % Slider is invisible or incorrectly oriented
                set( vSlider, 'Style', 'text', 'Enable', 'inactive', ...
                    'Position', vSliderPosition, ...
                    'Min', 0, 'Max', 1, 'Value', 1 )
            else
                % Compute properties
                vSliderMin = 0;
                vSliderMax = contentsHeight - vSliderHeight;
                vSliderValue = -vSlider.Value; % negative sign convention
                vSliderValue = max( vSliderValue, vSliderMin ); % limit
                vSliderValue = min( vSliderValue, vSliderMax ); % limit
                vStep = obj.VerticalSteps_(selection);
                vSliderStep(1) = min( vStep / vSliderMax, 1 );
                vSliderStep(2) = max( vSliderHeight / vSliderMax, vSliderStep(1) );
                contentsPosition(2) = contentsPosition(2) + vSliderValue;
                % Set properties
                set( vSlider, 'Style', 'slider', 'Enable', 'on', ...
                    'Position', vSliderPosition, ...
                    'Min', -vSliderMax, 'Max', -vSliderMin, ...
                    'Value', -vSliderValue, 'SliderStep', vSliderStep )
            end

            % Compute and set horizontal slider properties
            if hSliderHeight == 0 || hSliderWidth == 0 || hSliderWidth <= hSliderHeight
                % Slider is invisible or incorrectly oriented
                set( hSlider, 'Style', 'text', 'Enable', 'inactive', ...
                    'Position', hSliderPosition, ...
                    'Min', -1, 'Max', 0, 'Value', -1 )
            else
                % Compute properties
                hSliderMin = 0;
                hSliderMax = contentsWidth - hSliderWidth;
                hSliderValue = hSlider.Value; % positive sign convention
                hSliderValue = max( hSliderValue, hSliderMin ); % limit
                hSliderValue = min( hSliderValue, hSliderMax ); % limit
                hStep = obj.HorizontalSteps_(selection);
                hSliderStep(1) = min( hStep / hSliderMax, 1 );
                hSliderStep(2) = max( hSliderWidth / hSliderMax, hSliderStep(1) );
                contentsPosition(1) = contentsPosition(1) - hSliderValue;
                % Set properties
                set( hSlider, 'Style', 'slider', 'Enable', 'on', ...
                    'Position', hSliderPosition, ...
                    'Min', hSliderMin, 'Max', hSliderMax, ...
                    'Value', hSliderValue, 'SliderStep', hSliderStep )
            end

            % Set contents and blanking plate positions
            uix.setPosition( child, contentsPosition, 'pixels' )
            set( plate, 'Position', platePosition )

        end % redraw

        function addChild( obj, child )
            %addChild  Add child
            %
            %  c.addChild(d) adds the child d to the container c.

            % Add to sizes
            obj.Widths_(end+1,:) = -1;
            obj.MinimumWidths_(end+1,:) = 1;
            obj.Heights_(end+1,:) = -1;
            obj.MinimumHeights_(end+1,:) = 1;
            obj.VerticalSliders(end+1,:) = uicontrol( ...
                'Internal', true, 'Parent', obj, 'Units', 'pixels', ...
                'Style', 'slider' );
            obj.HorizontalSliders(end+1,:) = uicontrol( ...
                'Internal', true, 'Parent', obj, 'Units', 'pixels', ...
                'Style', 'slider' );
            obj.BlankingPlates(end+1,:) = uicontrol( ...
                'Internal', true, 'Parent', obj, 'Units', 'pixels', ...
                'Style', 'text', 'Enable', 'inactive' );
            obj.VerticalSteps_(end+1,:) = obj.SliderStep;
            obj.HorizontalSteps_(end+1,:) = obj.SliderStep;
            obj.updateSliderListeners()

            % Call superclass method
            addChild@uix.mixin.Panel( obj, child )

        end % addChild

        function removeChild( obj, child )
            %removeChild  Remove child
            %
            %  c.removeChild(d) removes the child d from the container c.

            % Remove from sizes
            tf = obj.Contents_ == child;
            obj.Widths_(tf,:) = [];
            obj.MinimumWidths_(tf,:) = [];
            obj.Heights_(tf,:) = [];
            obj.MinimumHeights_(tf,:) = [];
            obj.VerticalSliders(tf,:) = [];
            obj.HorizontalSliders(tf,:) = [];
            obj.BlankingPlates(tf,:) = [];
            obj.VerticalSteps_(tf,:) = [];
            obj.HorizontalSteps_(tf,:) = [];
            obj.updateSliderListeners()

            % Call superclass method
            removeChild@uix.mixin.Panel( obj, child )

        end % removeChild

        function reparent( obj, ~, newFigure )
            %reparent  Reparent container
            %
            %  c.reparent(a,b) reparents the container c from the figure a
            %  to the figure b.

            if isempty( newFigure )
                obj.MouseWheelListener = [];
            else
                listener = event.listener( newFigure, ...
                    'WindowScrollWheel', @obj.onMouseScrolled );
                listener.Enabled = strcmp( obj.MouseWheelEnabled_, 'on' );
                obj.MouseWheelListener = listener;
            end

        end % reparent

        function reorder( obj, indices )
            %reorder  Reorder contents
            %
            %  c.reorder(i) reorders the container contents using indices
            %  i, c.Contents = c.Contents(i).

            % Reorder
            obj.Widths_ = obj.Widths_(indices,:);
            obj.MinimumWidths_ = obj.MinimumWidths_(indices,:);
            obj.Heights_ = obj.Heights_(indices,:);
            obj.MinimumHeights_ = obj.MinimumWidths_(indices,:);
            obj.VerticalSliders = obj.VerticalSliders(indices,:);
            obj.HorizontalSliders = obj.HorizontalSliders(indices,:);
            obj.BlankingPlates = obj.BlankingPlates(indices,:);
            obj.VerticalSteps_ = obj.VerticalSteps_(indices,:);
            obj.HorizontalSteps_ = obj.HorizontalSteps_(indices,:);

            % Call superclass method
            reorder@uix.mixin.Panel( obj, indices )

        end % reorder

        function showSelection( obj )
            %showSelection  Show selected child, hide the others
            %
            %  c.showSelection() shows the selected child of the container
            %  c, and hides the others.

            % Call superclass method
            showSelection@uix.mixin.Panel( obj )

            % Show and hide sliders based on selection
            selection = obj.Selection_;
            for ii = 1:numel( obj.Contents_ )
                if ii == selection
                    obj.VerticalSliders(ii).Visible = 'on';
                    obj.HorizontalSliders(ii).Visible = 'on';
                    obj.BlankingPlates(ii).Visible = 'on';
                else
                    obj.VerticalSliders(ii).Visible = 'off';
                    obj.HorizontalSliders(ii).Visible = 'off';
                    obj.BlankingPlates(ii).Visible = 'off';
                end
            end

        end % showSelection

    end % template methods

    methods( Access = private )

        function onSliderScrolling( obj, ~, ~ )
            %onSliderScrolling  Event handler

            % Mark as dirty
            obj.Dirty = true;

            % Raise event
            notify( obj, 'Scrolling' )

        end % onSliderScrolling

        function onSliderScrolled( obj, ~, ~ )
            %onSliderScrolled  Event handler

            % Mark as dirty
            obj.Dirty = true;

            % Raise event
            notify( obj, 'Scrolled' )

        end % onSliderScrolled

        function onMouseScrolled( obj, ~, eventData )
            %onMouseScrolled  Event handler

            sel = obj.Selection_;
            if sel == 0
                return
            else
                % Get pointer position and panel bounds
                pp = getpixelposition( obj, true );
                f = ancestor( obj, 'figure' );
                cpu = f.CurrentPoint; % figure Units
                cpwhu = [cpu 0 0]; % [x y] to [x y w h] for hgconvertunits
                cpwh = hgconvertunits( f, cpwhu, f.Units, 'pixels', obj ); % pixels
                cp = cpwh(1:2); % [x y w h] to [x y]

                % Check that pointer is over panel
                if cp(1) < pp(1) || cp(1) > pp(1) + pp(3) || ...
                        cp(2) < pp(2) || cp(2) > pp(2) + pp(4), return, end
                % Scroll
                if strcmp( obj.VerticalSliders(sel).Enable, 'on' ) % scroll vertically
                    delta = eventData.VerticalScrollCount * ...
                        eventData.VerticalScrollAmount * obj.VerticalSteps(sel);
                    obj.VerticalOffsets(sel) = obj.VerticalOffsets(sel) + delta;
                elseif strcmp( obj.HorizontalSliders(sel).Enable, 'on' ) % scroll horizontally
                    delta = eventData.VerticalScrollCount * ...
                        eventData.VerticalScrollAmount * obj.HorizontalSteps(sel);
                    obj.HorizontalOffsets(sel) = obj.HorizontalOffsets(sel) + delta;
                end
                % Raise event
                notify( obj, 'Scrolled' )
            end

        end % onMouseScrolled

    end % event handlers

    methods( Access = private )

        function updateSliderListeners( obj )
            %updateSliderListeners  Update listeners to slider events

            if isempty( obj.VerticalSliders )
                obj.ScrollingListener = [];
                obj.ScrolledListener = [];
            else
                obj.ScrollingListener = event.listener( ...
                    [obj.VerticalSliders; obj.HorizontalSliders], ...
                    'ContinuousValueChange', @obj.onSliderScrolling );
                obj.ScrolledListener = event.listener( ...
                    [obj.VerticalSliders; obj.HorizontalSliders], ...
                    'Action', @obj.onSliderScrolled );
            end

        end % updateSliderListeners

    end % helpers

end % classdef

author: goGPS-Project


text:
classdef inversion_recovery < AbstractModel
%inversion_recovery: Compute a T1 map using Inversion Recovery data
%
% Assumptions:
% (1) Gold standard for T1 mapping
% (2) Infinite TR
%
% Inputs:
%   IRData      Inversion Recovery data (4D)
%   (Mask)      Binary mask to accelerate the fitting (OPTIONAL)
%
% Outputs:
%   T1          transverse relaxation time [ms]
%   b           arbitrary fit parameter (S=a + b*exp(-TI/T1))
%   a           arbitrary fit parameter (S=a + b*exp(-TI/T1))
%   idx         index of last polarity restored datapoint (only used for magnitude data)
%   res         Fitting residual
%
%
% Protocol:
%	IRData  [TI1 TI2...TIn] inversion times [ms]
%   TimingTable   [TR] repetition time [ms]
%
% Options:
%   method          Method to use in order to fit the data, based on whether complex or only magnitude data acquired.
%     'complex'     Complex dataset.
%     'magnitude'   Magnitude dataset.
%
%
%   fitModel        T1 fitting moddel.
%     'Barral'      Fitting equation: a+bexp(-TI/T1)
%     'General'     Fitting equation: c(1-2exp(-TI/T1)+exp(-TR/T1))
%
% Example of command line usage (see also showdemo inversion_recovery_batch):
%   Model = inversion_recovery;  % Create class from model
%   Model.Prot.IRData.Mat=[350.0000; 500.0000; 650.0000; 800.0000; 950.0000; 1100.0000; 1250.0000; 1400.0000; 1700.0000];
%   data = struct;  % Create data structure
%   data.MET2data ='IRData.mat';  % Load data
%   data.Mask = 'Mask.mat';
%   FitResults = FitData(data,Model); %fit data
%   FitResultsSave_mat(FitResults);
%
%       For more examples: qMRusage(inversion_recovery)
%
% Author: Ilana Leppert, 2017
%
% References:
%   Please cite the following if you use this module:
%       A robust methodology for in vivo T1 mapping. Barral JK, Gudmundson E, Stikov N, Etezadi-Amoli M, Stoica P, Nishimura DG. Magn Reson Med. 2010 Oct;64(4):1057-67. doi: 10.1002/mrm.22497.
%   In addition to citing the package:
%     Karakuzu A., Boudreau M., Duval T.,Boshkovski T., Leppert I.R., Cabana J.F.,
%     Gagnon I., Beliveau P., Pike G.B., Cohen-Adad J., Stikov N. (2020), qMRLab:
%     Quantitative MRI analysis, under one umbrella doi: 10.21105/joss.02343

properties (Hidden=true)
    onlineData_url = 'https://osf.io/cmg9z/download?version=3';
end

	properties
        MRIinputs = {'IRData','Mask'}; % input data required
        xnames = {'T1','rb','ra'}; % name of the fitted parameters
        voxelwise = 1; % voxel by voxel fitting?

        % fitting options
        st           = [  600    -1000      500 ]; % starting point
        lb           = [    0.0001   -10000        0.0001 ]; % lower bound
        ub           = [ 5000        0    10000 ]; % upper bound
        fx           = [    0        0        0 ]; % fix parameters

        % Protocol
        Prot = struct('IRData', struct('Format',{'TI(ms)'},'Mat',[350 500 650 800 950 1100 1250 1400 1700]'),...
                      'TimingTable', struct('Format',{{'TR(ms)'}},'Mat',2500)); %default protocol
        % Model options
        buttons = {'method',{'Magnitude','Complex'}, 'fitModel',{'Barral','General'}}; %selection buttons
        options = struct(); % structure filled by the buttons. Leave empty in the code

        % Simulation Options
        Sim_Single_Voxel_Curve_buttons = {'SNR',50,'T1',600,'M0',1000,'TR',3000,'FAinv',180,'FAexcite',90,'Update input variables','pushbutton'};%'FArefocus',180
        Sim_Optimize_Protocol_buttons = {'# of volumes',5,'Population size',100,'# of migrations',100};

    end

methods (Hidden=true)
% Hidden methods goes here.
end

    methods
        % -------------CONSTRUCTOR-------------------------------------------------------------------------
        function  obj = inversion_recovery()
            obj.options = button2opts(obj.buttons);
        end

        function obj = UpdateFields(obj)
            obj.Prot.IRData.Mat = sort(obj.Prot.IRData.Mat);
        end
        function xnew = SimOpt(obj,x,Opt)
            [ra,rb] = ComputeRaRb(obj,x,Opt);
            xnew = [Opt.T1 rb ra];

        end

        % -------------IR EQUATION-------------------------------------------------------------------------
        function Smodel = equation(obj, x)
            % Generates an IR signal based on fit parameters
            x = mat2struct(x,obj.xnames); % if x is a structure, convert to vector

            % equation
            Smodel = x.ra + x.rb * exp(-obj.Prot.IRData.Mat./x.T1);
            if (strcmp(obj.options.method, 'Magnitude'))
                Smodel = abs(Smodel);
            end
        end

        % -------------EXPLICIT IR EQUATION-------------------------------------------------------------------------
        function [ra,rb] = ComputeRaRb(obj,x,Opt)

            % Some sanity checks
            [ErrMsg]=[];

            for brkloop=1:1
                if Opt.TR < max(obj.Prot.IRData.Mat) %TR can't be less than max TI
                    txt=['The TR (' num2str(Opt.TR) ') cannot be less than max TI (' num2str(max(obj.Prot.IRData.Mat)),')'];
                    ErrMsg = txt; break
                end
                if Opt.T1 < 0 || Opt.T1 > 10000
                    txt='Choose a reasonable value for T1 (0-10000 s)';
                    ErrMsg = txt; break
                end
                if Opt.FAinv < 120 || Opt.FAinv > 220
                    txt='Choose a reasonable value for the inversion FA (120-220 deg)';
                    ErrMsg = txt; break
                end
                 if Opt.FAexcite < 50 || Opt.FAexcite > 120
                    txt='Choose a reasonable value for the excitation FA (50-120 deg)';
                    ErrMsg = txt; break
                end
            end
            if ~isempty(ErrMsg)
                if moxunit_util_platform_is_octave
                    errordlg(ErrMsg,'Input Error');
                else
                    Mode = struct('WindowStyle','modal','Interpreter','tex');
                    errordlg(ErrMsg,'Input Error', Mode);
                    error(ErrMsg);
                end
            end

            % equation for GRE-IR
            ra = Opt.M0 * (1-cos(Opt.FAinv*pi/180)*exp(-Opt.TR/Opt.T1))/(1-cos(Opt.FAinv*pi/180)*cos(Opt.FAexcite*pi/180)*exp(-Opt.TR/Opt.T1));
            rb = -Opt.M0 * (1-cos(Opt.FAinv*pi/180))/(1-cos(Opt.FAinv*pi/180)*cos(Opt.FAexcite*pi/180)*exp(-Opt.TR/Opt.T1));
            %Smodel = ra + rb * exp(-obj.Prot.IRData.Mat./x.T1);
            %if (strcmp(obj.options.method, 'Magnitude'))
            %    Smodel = abs(Smodel);
            %end
        end

        % -------------DATA FITTING-------------------------------------------------------------------------
        function FitResults = fit(obj,data)
            % Fits the data
            %

            data = data.IRData;

            switch obj.options.fitModel
                case 'Barral'
                    [T1,rb,ra,res,idx] = fitT1_IR(data,obj.Prot.IRData.Mat,obj.options.method);
                    FitResults.T1  = T1;
                    FitResults.rb  = rb;
                    FitResults.ra  = ra;
                    FitResults.res = res;
                    if (strcmp(obj.options.method, 'Magnitude'))
                        FitResults.idx = idx;
                    end
                case 'General'
                    approxFlag = 3; % Selects the equation c(1-2exp(-TI/T1)+exp(TR/T1))

                    params.TI = obj.Prot.IRData.Mat';
                    params.TR = obj.Prot.TimingTable.Mat;

                    if strcmp(obj.options.method, 'Magnitude')
                        % Make sure data vector is a column vector
                        data = data(:);

                        % Find the min of the data (which is nearest to
                        % signal null to flip
                        [~, minInd] = min(data);

                        % Signal inversion algorithm to fir the T1 curve
                        for ii = 1:2
                            % Fit the data by inverting the sign of the
                            % datapoints up until the TI where signal is
                            % null, then also without that point. The
                            % fit with the smallest residual is our best
                            % guess for up to where to flip the sign of the
                            % signal.
                            if ii == 1
                                % First, we set all elements up to and including
                                % the smallest element to minus
                                dataTmp = data.*[-ones(minInd,1); ones(length(data) - minInd,1)];
                            elseif ii == 2
                                % Second, we set all elements up to (not including)
                                % the smallest element to minus
                                dataTmp = data.*[-ones(minInd-1,1); ones(length(data) - (minInd-1),1)];
                            end
                            [fitVals{ii}, resnorm(ii)] = inversion_recovery.fit_lm(dataTmp, params, approxFlag);
                        end
                        [~,ind] = min(resnorm); % Index of the minimum residual will be which signal fit results to choose.
                        FitResults.T1 = fitVals{ind}.T1;
                        FitResults.ra = fitVals{ind}.ra;
                        FitResults.rb = fitVals{ind}.rb;
                        FitResults.res = resnorm(ind);
                        FitResults.idx = ind;
                    elseif strcmp(obj.options.method, 'Complex')
                        params.dataType = 'complex';
                        [fitVals, resnorm] = inversion_recovery.fit_lm(data(:), params, 3);
                        FitResults.T1 = fitVals.T1;
                        FitResults.ra = fitVals.ra;
                        FitResults.rb = fitVals.rb;
                        FitResults.res = resnorm;
                    end
            end
        end

        function plotModel(obj, FitResults, data)
            % Plots the fit
            %
            % :param FitResults: [struct] Fitting parameters
            % :param data: [struct] input data
            if nargin<2 || isempty(FitResults), FitResults = obj.st; end
            if exist('data','var')
                data = data.IRData;
                % plot
                plot(obj.Prot.IRData.Mat,data,'.','MarkerSize',15)
                hold on
                if (strcmp(obj.options.method, 'Magnitude'))
                    % plot the polarity restored data points
                    data_rest = -1.*data(1:FitResults.idx);
                    plot(obj.Prot.IRData.Mat(1:FitResults.idx),data_rest,'o','MarkerSize',5,'MarkerEdgeColor','b','MarkerFaceColor',[1 0 0])
                end
            end

            % compute model
            obj.Prot.IRData.Mat = linspace(min(obj.Prot.IRData.Mat),max(obj.Prot.IRData.Mat),100);
            Smodel = equation(obj, FitResults);

            % plot fitting curve
            plot(obj.Prot.IRData.Mat,Smodel,'Linewidth',3)
            hold off
            xlabel('Inversion Time [ms]','FontSize',15);
            ylabel('Signal','FontSize',15);
            legend('data', 'polarity restored', 'fit','Location','best')
            set(gca,'FontSize',15)
        end

        function [FitResults, data] = Sim_Single_Voxel_Curve(obj, x, Opt,display)
            % Simulates Single Voxel
            %
            % :param x: [struct] fit parameters
            % :param Opt.SNR: [struct] signal to noise ratio to use
            % :param display: 1=display, 0=nodisplay
            % :returns: [struct] FitResults, data (noisy dataset)

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

            Smodel = equation(obj, x);
            sigma = max(abs(Smodel))/Opt.SNR;
            if (strcmp(obj.options.method, 'Magnitude'))
                data.IRData = ricernd(Smodel,sigma);
            else
                data.IRData = random('normal',Smodel,sigma);
            end
            FitResults = fit(obj,data);
            if display
                plotModel(obj, FitResults, data);
            end
        end

        function SimVaryResults = Sim_Sensitivity_Analysis(obj, OptTable, Opt)
            % SimVaryGUI
            SimVaryResults = SimVary(obj, Opt.Nofrun, OptTable, Opt);
        end

        function SimRndResults = Sim_Multi_Voxel_Distribution(obj, RndParam, Opt)
            % SimRndGUI
            SimRndResults = SimRnd(obj, RndParam, Opt);
        end

%         function schemeLEADER = Sim_Optimize_Protocol(obj,xvalues,Opt)
%             % schemeLEADER = Sim_Optimize_Protocol(obj,xvalues,nV,popSize,migrations)
%             % schemeLEADER = Sim_Optimize_Protocol(obj,obj.st,30,100,100)
%             % Optimize Inversion times
%             nV         = Opt.Nofvolumes;
%             popSize    = Opt.Populationsize;
%             migrations = Opt.Nofmigrations;
%
%             sigma  = .05;
%             TImax = 5000;
%             TImin = 50;
%             GenerateRandFunction = @() rand(nV,1)*(TImax-TImin)+TImin; % do not sort TI values... or you might fall in a local minima
%             CheckProtInBoundFunc = @(Prot) min(max(50,Prot),TImax);
%             % Optimize Protocol
%             [retVal] = soma_all_to_one(@(Prot) mean(SimCRLB(obj,Prot,xvalues,sigma)), GenerateRandFunction, CheckProtInBoundFunc, migrations, popSize, nV, obj.Prot.IRData.Mat(:,1));
%
%             % Generate Rest
%             schemeLEADER = retVal.schemeLEADER;
%
%             fprintf('SOMA HAS FINISHED \n')
%
%         end

    end

    % CLI-only implemented static methods. Can be called directly from
    % class - no object needed.
    methods(Static)
        function Mz = analytical_solution(params, seqFlag, approxFlag)
            %ANALYTICAL_SOLUTION  Analytical equations for the longitudinal magnetization of
            %steady-state inversion recovery experiments with either a gradient echo
            %(GRE-IR) or spin-echo (SE-IR) readouts.
            %   Reference: Barral, J. K., Gudmundson, E. , Stikov, N. , Etezadi?Amoli,
            %   M. , Stoica, P. and Nishimura, D. G. (2010), A robust methodology for
            %   in vivo T1 mapping. Magn. Reson. Med., 64: 1057-1067.
            %   doi:10.1002/mrm.22497
            %
            %   params: struct with the required parameters for the sequence and
            %   approximation. See below for list.
            %
            %   seqFlag: String. Either 'GRE-IR' or 'SE-IR'
            %   approxFlag: Integer between 1 and 4.
            %       1: General equation (no approximation).
            %       2: Ideal 180 degree pulse approximation of case 1.
            %       3: Ideal 90 degree pulse approximation of case 2, and readout term
            %          absorbed into constant.
            %       4: Long TR (TR >> T1) approximation of case 3.
            %
            %   **PARAMS PROPERTIES**
            %   All times in ms, all angles in degrees.
            %  'GRE-IR'
            %       case 1: T1, TR, TI, EXC_FA, INV_FA, constant (optional)
            %       case 2: T1, TR, TI, EXC_FA, constant (optional)
            %       case 3: T1, TR, TI, constant (optional)
            %       case 4: T1, TI, constant (optional)
            %
            %  'SE-IR'
            %       case 1: Same as 'GRE-IR' case + SE_FA, TE
            %       case 2: Same as 'GRE-IR' case + TE
            %       case 3: Same as 'GRE-IR' case + TE
            %       case 4: Same as 'GRE-IR' case
            %

            Mz = ir_equations(params, seqFlag, approxFlag);

        end

        function [Mz, Msig] = bloch_sim(params)
            %BLOCH_SIM Bloch simulations of the GRE-IR pulse sequence.
            % Simulates 100 spins params.Nex repetitions of the IR pulse
            % sequences.
            %
            % params: Struct with the following fields:
            %   INV_FA: Inversion pulse flip angle in degrees.
            %   EXC_FA: Excitation pulse flip angle in degrees.
            %   TI: Inversion time (ms).
            %   TR: Repetition time (ms).
            %   TE: Echo time (ms).
            %   T1: Longitudinal relaxation time (ms).
            %   T2: Transverse relaxation time (ms).
            %   Nex: Number of excitations
            %
            %   (optional
            %       df: Off-resonance frequency of spins relative to excitation pulse (in Hz)
            %       crushFlag: Numeric flag for perfect spoiling (1) or partial spoiling (2).
            %       partialDephasing: Partial dephasing fraction (between [0, 1]). 1 = no dephasing, 0 = complete dephasing (sele
            %       inc: Phase spoiling increment in degrees.
            %
            % Outputs:
            %   Mz: Longitudinal magnetization at time TI (prior to excitation pulse).
            %   Msig: Complex signal produced by the transverse magnetization at time TE after excitation.
            %

            %% Setup parameters
            %

            alpha = deg2rad(params.INV_FA);
            beta  = deg2rad(params.EXC_FA);
            TI = params.TI;
            TR = params.TR;
            T1 = params.T1;

            TE = params.TE;
            T2 = params.T2;

            Nex = params.Nex;

            %% Optional parameers

            if isfield(params, 'df')
                df = params.df;
            else
                df = 0;
            end

            if isfield(params, 'crushFlag')
                crushFlag = params.crushFlag;
            else
                crushFlag = 1;
            end

            if isfield(params, 'partialDephasing')
                partialDephasing = params.partialDephasing;
            else
                partialDephasing = 1;
            end

            if isfield(params, 'inc')
                inc = deg2rad(params.inc);
            else
                inc = 0;
            end

            %% Simulate for every TI's
            %

            for ii = 1:length(TI)

                [Msig(ii), Mz(ii)] = ir_blochsim(                  ...
                                                 alpha,            ...
                                                 beta,             ...
                                                 TI(ii),           ...
                                                 T1,               ...
                                                 T2,               ...
                                                 TE,               ...
                                                 TR,               ...
                                                 crushFlag,        ...
                                                 partialDephasing, ...
                                                 df,               ...
                                                 Nex,              ...
                                                 inc               ...
                                                 );

            end

        end

        function [fitVals, resnorm] = fit_lm(data, params, approxFlag)
            % FIT_LM Levenberg-Marquardt fitting of GRE-IR data.
            %
            % data: Array for a single voxel (length = #TI)
            % params: Properties TI and TR (except for approxFlag = 4,
            %         where only TI is needed)
            % approxFlag: Same flags numbering & equations as for the
            %             analytical_solution class method.

            switch approxFlag
                case 1
                    TI = params.TI;
                    TR = params.TR;

                    %    [constant, T1, EXC_FA, INV_FA]
                    x0 = [1, 1000, 90, 180];

                    options.Algorithm = 'levenberg-marquardt';
                    options.Display = 'off';

                    [x, resnorm] = lsqnonlin(@(x)ir_loss_func_1(x, TR, TI, data), x0, [], [], options);

                    fitVals.INV_FA = x(4);
                    fitVals.EXC_FA = x(3);
                    fitVals.T1 = x(2);
                    fitVals.c = x(1);
                case 2
                    TI = params.TI;
                    TR = params.TR;

                    %    [constant, T1, EXC_FA]
                    x0 = [1, 1000, 90];

                    options.Algorithm = 'levenberg-marquardt';
                    options.Display = 'off';

                    [x, resnorm] = lsqnonlin(@(x)ir_loss_func_2(x, TR, TI, data), x0, [], [], options);

                    fitVals.EXC_FA = x(3);
                    fitVals.T1 = x(2);
                    fitVals.c = x(1);

                case 3
                    TI = params.TI;
                    TR = params.TR;

                    % Find the min of the data
                    [~, minInd] = min(abs(data));

                    T1est = -TI(minInd)/log(1/2); % Estimate from null (or min) of data and the equation in Case 4. Used as the first guess for the fitting point.

                    if isfield(params, 'dataType') && strcmp(params.dataType,'complex')
                        % Fit comlex data
                        %    [constant, T1]
                        if max(abs(data)) > 0
                            dataNorm = data/max(abs(data)); % Normalize the data to fit with the initial constant guess of 1.
                        else
                            dataNorm = data;
                        end

                        % Initial fit estimates
                        %    [Re(constant), Im(constant), T1)]
                        x0 = [1, 0, T1est];

                        options.Algorithm = 'trust-region-reflective';
                        options.Display = 'off';

                        % Bound the fitting variables to -2 to 2
                        % (constant) and 0 to 5000 (T1 in ms).
                        [x, resnorm] = lsqnonlin(@(x)ir_loss_func_3(x, TR, TI, dataNorm(:)', params.dataType), x0, [-2, -2, 0], [2, 2, 5000], options);

                        fitVals.T1 = x(3);

                        % ra and rb calculation from Equation 3
                        if max(abs(data)) > 0
                            fitVals.ra = (x(1)+1i*x(2))*max(abs(data)) * (1 + exp(-TR/fitVals.T1));
                            fitVals.rb = -2*(x(1)+1i*x(2))*max(abs(data));
                        else
                            fitVals.ra = (x(1)+1i*x(2)) * (1 + exp(-TR/fitVals.T1));
                            fitVals.rb = -2*(x(1)+1i*x(2));
                        end
                    else
                        % Fit magnitude data
                        %    [constant, T1]
                        x0 = [1, T1est];
                        if max(abs(data)) > 0
                            dataNorm = data/max(abs(data)); % Normalize the data to fit with the initial constant guess of 1.
                        else
                            dataNorm = data;
                        end

                        options.Algorithm = 'trust-region-reflective';
                        options.Display = 'off';

                        % Bound the fitting variables to -2 to 2
                        % (constant) and 0 to 5000 (T1 in ms).
                        [x, resnorm] = lsqnonlin(@(x)ir_loss_func_3(x, TR, TI, dataNorm(:)'), x0, [0, 0], [2, 5000], options);

                        fitVals.T1 = x(2);

                        % ra and rb calculation from Equation 3
                        if max(abs(data)) > 0
                            fitVals.ra = x(1)*max(abs(data)) * (1 + exp(-TR/fitVals.T1));
                            fitVals.rb = -2*x(1)*max(abs(data));
                        else
                            fitVals.ra = x(2) * (1 + exp(-TR/fitVals.T1));
                            fitVals.rb = -2*x(2);
                        end

                    end

                case 4
                    TI = params.TI;
                    %    [constant, T1]
                    x0 = [1, 1000];

                    options.Algorithm = 'levenberg-marquardt';
                    options.Display = 'off';

                    [x, resnorm] = lsqnonlin(@(x)ir_loss_func_4(x, TI, data), x0, [], [], options);

                    fitVals.T1 = x(2);
                    fitVals.c = x(1);
            end

        end

    end

    methods(Access = protected)
        function obj = qMRpatch(obj,loadedStruct, version)
            obj = qMRpatch@AbstractModel(obj,loadedStruct, version);
            % 2.0.10
            if checkanteriorver(version,[2 0 10])
                % Update buttons for simulation
                snrValue = obj.Sim_Single_Voxel_Curve_buttons{2};
                obj.Sim_Single_Voxel_Curve_buttons = {'SNR' snrValue 'T1' 600 'M0' 1000 'TR' 3000 'FAinv' 180 'FAexcite' 90 'Update input variables' 'pushbutton'};
            end
             if checkanteriorver(version,[2 4 1])
                % Update buttons for simulation
                obj.Prot = struct('IRData', struct('Format',{'TI(ms)'},'Mat',[350 500 650 800 950 1100 1250 1400 1700]'),...
                                  'TimingTable', struct('Format',{{'TR(ms)'}},'Mat',2500)); %default protocol
                % Model options
                obj.buttons = {'method',{'Magnitude','Complex'}, 'fitModel',{'Barral','General'}}; %selection buttons
                obj.options = button2opts(obj.buttons);
            end
        end
    end

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% function failed_tutorial_fem

% WALLTIME 01:00:00
% MEM 6gb
% DEPENDENCY ft_volumesegment ft_prepare_mesh ft_prepare_headmodel ft_headmodel_simbio

% this test script is based on the tutorial under development at
% http://www.fieldtriptoolbox.org/development/simbio

cd(dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf'))

mri = ft_read_mri('Subject01.mri');

figure
cfg = [];
ft_sourceplot(cfg, mri);

cfg     = [];
cfg.dim = mri.dim;
mri     = ft_volumereslice(cfg, mri);

figure
cfg = [];
ft_sourceplot(cfg,mri);

cfg           = [];
cfg.output    = {'gray','white','csf','skull','scalp'};
segmentedmri  = ft_volumesegment(cfg, mri);

seg_i            = ft_datatype_segmentation(segmentedmri, 'segmentationstyle', 'indexed');

cfg              = [];
cfg.funparameter = 'seg';
cfg.funcolormap  = lines(6);
cfg.location     = 'center';
cfg.atlas        = seg_i;
ft_sourceplot(cfg, seg_i);

cfg        = [];
cfg.shift  = 0.3;
cfg.method = 'hexahedral';
mesh = ft_prepare_mesh(cfg, segmentedmri);

cfg              = [];
cfg.method       = 'simbio';
cfg.conductivity = [0.33 0.14 1.79 0.01 0.43];           % order follows mesh.tissuelabel
vol = ft_prepare_headmodel(cfg, mesh);
% the call to "ft_prepare_headmodel" took 320 seconds and required the additional allocation of an estimated 501 MB

% you may need to specify the full path to the file
elec = ft_read_sens('standard_1020.elc');

figure
hold on
ft_plot_mesh(mesh,'surfaceonly','yes','vertexcolor','none','edgecolor','none','facecolor',[0.5 0.5 0.5],'face alpha',0.5)
camlight
ft_plot_sens(elec,'style', 'sk');

nas = mri.hdr.fiducial.mri.nas;
lpa = mri.hdr.fiducial.mri.lpa;
rpa = mri.hdr.fiducial.mri.rpa;

vox2head = mri.transform;

nas = ft_warp_apply(vox2head, nas, 'homogenous');
lpa = ft_warp_apply(vox2head, lpa, 'homogenous');
rpa = ft_warp_apply(vox2head, rpa, 'homogenous');

% create a structure similar to a template set of electrodes
fid.chanpos       = [nas; lpa; rpa];       % CTF head coordinates of fiducials
fid.label         = {'Nz','LPA','RPA'};    % use the same labels as those in elec
fid.unit          = 'mm';                  % use the same units as those in mri

% alignment
cfg               = [];
cfg.method        = 'fiducial';
cfg.elec          = elec;                  % the electrodes we want to align
cfg.template      = fid;                   % the template we want to align to
cfg.fiducial      = {'Nz', 'LPA', 'RPA'};  % labels of fiducials in fid and in elec
elec_aligned      = ft_electroderealign(cfg);

figure;
hold on;
ft_plot_mesh(mesh,'surfaceonly','yes','vertexcolor','none','edgecolor','none','facecolor',[0.5 0.5 0.5],'face alpha',0.5)
camlight
ft_plot_sens(elec_aligned,'style','sr');

if false
  % this should only run in interactive mode
  cfg = [];
  cfg.method = 'interactive';
  cfg.elec = elec_aligned;
  cfg.headshape = vol;
  elec_aligned2 = ft_electroderealign(cfg);
else
  % this is more or less the same as the interactive alignment
  transform = [
    1 0 0 15
    0 1 0 0
    0 0 1 -10
    0 0 0 1
    ];
  elec_aligned2 = ft_transform_geometry(transform, elec_aligned);
end

%%%%%%%%%

cfg = [];
cfg.headmodel = vol;
cfg.sourcemodel.resolution = 10;
cfg.elec = elec_aligned2;
cfg.sourcemodel.unit = 'mm';
source = ft_prepare_sourcemodel(cfg);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip